Skip to content
Snippets Groups Projects
Commit c0ab8db1 authored by Harald.Sodemann's avatar Harald.Sodemann
Browse files

Merge branch 'correct_ventilation_factor' into 'master'

Correct heavy isotope mass tendency equations

See merge request !2
parents c5f8352f 3710a9f3
Branches master
Tags V1.1
1 merge request!2Correct heavy isotope mass tendency equations
......@@ -91,22 +91,22 @@ end
% -----------------------------------------------------------------------------
if strcmp(namelist.isocalc,'stewart') % Calculation according to Stewart (1975)
% relative humidity of the environment with respect to the drop temperature
rh=profiles.equih;
% BO and BH are used below
BO=alphaO(Tint).*DDo.*(1-rh);
BH=alphaH(Tint).*DDh.*(1-rh);
% (3e) in Stewart (1975)
gammaO=alphaO(Tint).*rh./(1-BO);
gammaH=alphaH(Tint).*rh./(1-BH);
% (3d) in Stewart (1975)
betaO=(1-BO)./(BO);
betaH=(1-BH)./(BH);
% Start calculation from top to bottom
for i=length(outarray):-1:2
if Tint(i)>0 && f(i)<=1 % Temperature positive and drop is shrinking (change of composition by growth is neglected)
......@@ -117,16 +117,16 @@ if strcmp(namelist.isocalc,'stewart') % Calculation according to Stewart (1975)
ROp(i-1)=ROp(i);
end
end
elseif strcmp(namelist.isocalc,'explicit') % Calculation using masses of all isotopic species rather than ratios
% Calculate initial mass of all isotopic species from the isotopic
% ratio
mh2o(length(outarray))=(1-RHp(end)-ROp(end))*profiles.m(end); % mass of HHO
mhdo(length(outarray))=RHp(end)*profiles.m(end); % mass of HDO
m18o(length(outarray))=ROp(end)*profiles.m(end); % mass of HH18O
mtot(length(outarray))=profiles.m(end); % total mass
mvaptot(length(outarray))=0; % evaporated water mass between current and previous altitude
RHvap(length(outarray))=0; % composition of the evaporated water
ROvap(length(outarray))=0;
......@@ -139,15 +139,11 @@ elseif strcmp(namelist.isocalc,'explicit') % Calculation using masses of all iso
Tint(i)=mean([Tint(i),Tint(i-1)]);
% Determine the time that the hydrometeor uses to dexcend from z(i) to z(i-1)
dt=(outarray(i)-outarray(i-1))/terminal_velocity(profiles.diam(i),p(i),T(i));
if Tint(i)>0 % Tint is positive: fractionation occurs
% Determine some coefficients
diff(i)=diffusion_coefficient(T(i),p(i)); % Diffusivity of the light isotope (H2O)
fmo=mass_vent_coeff(profiles.diam(i),T(i),p(i),rh(i),'heavyO')/...
mass_vent_coeff(profiles.diam(i),T(i),p(i),rh(i)); % Ratio of mass ventilation coefficients for 18O
fmh=mass_vent_coeff(profiles.diam(i),T(i),p(i),rh(i),'heavyH')/...
mass_vent_coeff(profiles.diam(i),T(i),p(i),rh(i)); % Ratio of mass ventilation coefficients for 2H
% mass rate of change of the light isotope ( dm(H2O)/dt );
% See Equations (3.11 and 3.22) in Graf 2017 (Diss ETH No. 24777)
dmldt(i-1)=2*pi*profiles.diam(i)*diff(i)/const.Rw*...
......@@ -156,17 +152,17 @@ elseif strcmp(namelist.isocalc,'explicit') % Calculation using masses of all iso
(1-RHv(i)-ROv(i))-... % ratio of light isotope in vapour
satvappress(Tint(i))./(Tint(i)+273.15)*...
(1-RHp(i)/alphaH(Tint(i))-ROp(i)/alphaO(Tint(i))));% ratio of light isotope in hydrometeor
mh2o(i-1)=mh2o(i)+dt*dmldt(i-1); % new mh2o
% Define factor with all constants (for readability)
F=2*pi*profiles.diam(i)*diff(i)*mass_vent_coeff(profiles.diam(i),T(i),p(i),rh(i))/const.Rw;
% Same as above, but for m18o and mhdo (see Eq. (3.22) in Graf 2017)
m18o(i-1)=m18o(i)+dt*F*(fmo/DDova)^nexp*(rh(i)*ROv(i)*satvappress(T(i))./(T(i)+273.15)-...
m18o(i-1)=m18o(i)+dt*F*(1./DDova)^nexp*(rh(i)*ROv(i)*satvappress(T(i))./(T(i)+273.15)-...
ROp(i)/alphaO(Tint(i))*satvappress(Tint(i))./(Tint(i)+273.15));
mhdo(i-1)=mhdo(i)+dt*F*(fmh/DDhva)^nexp*(rh(i)*RHv(i)*satvappress(T(i))./(T(i)+273.15)-...
mhdo(i-1)=mhdo(i)+dt*F*(1./DDhva)^nexp*(rh(i)*RHv(i)*satvappress(T(i))./(T(i)+273.15)-...
RHp(i)/alphaH(Tint(i))*satvappress(Tint(i))./(Tint(i)+273.15));
% Total change of mass (Eq. 3.11)
mtot(i-1)=mtot(i)+dt*2*pi*profiles.diam(i)*diff(i)/const.Rw*...
mass_vent_coeff(profiles.diam(i),T(i),p(i),rh(i))*...
......@@ -175,13 +171,13 @@ elseif strcmp(namelist.isocalc,'explicit') % Calculation using masses of all iso
% Calculate new isotopic ratios
RHp(i-1)=mhdo(i-1)/(m18o(i-1)+mhdo(i-1)+mh2o(i-1));
ROp(i-1)=m18o(i-1)/(m18o(i-1)+mhdo(i-1)+mh2o(i-1));
% Calculate evaporated mass
mvaptot(i-1)=mtot(i)-mtot(i-1);
% Calculate isotopic ratio of the evaporated mass
RHvap(i-1)=(mhdo(i)-mhdo(i-1))/(m18o(i)-m18o(i-1)+mhdo(i)-mhdo(i-1)+mh2o(i)-mh2o(i-1));
ROvap(i-1)=(m18o(i)-m18o(i-1))/(m18o(i)-m18o(i-1)+mhdo(i)-mhdo(i-1)+mh2o(i)-mh2o(i-1));
else % Tint is not positive
if profiles.m(i-1)-profiles.m(i)<=0 % mass decreases:
% Relative change of mass is the same for all isotope species (no fractionation)
......@@ -253,4 +249,4 @@ if namelist.fulloutput==1
profiles.RHvap=RHvap;
profiles.ROvap=ROvap;
end
end
\ No newline at end of file
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment