Basics

Diagnostic equations to calculate necessary parameters for microphysics from prognostic variables.

\( \rho_{d} = \rho – \rho_{v} – \rho_{c} – \rho_{r} – \rho_{i} – \rho_{s}\)
\(C_{pml} = (\rho_{d} \cdot 1004. + \rho_{v} \cdot 1885. + (\rho_{c}+\rho_{r}) \cdot 4186. + (\rho_{i} + \rho_{s}) \cdot 2110.\)
\(R_{ml}= (/rho_{d} \cdot 287. + rho_{v} \cdot * 461.)\)
\(\theta = exp(( \mu \theta + R_{ml}*log(p_{0}))/C_{pml}\)
\(T_{s} = \frac{\mu \theta}{(C_{pml} – R_{ml})}+ \frac{1.}{C_{pml/rml} -1.} \cdot log(R_{ml})\)
\(p_{s}= \frac{\mu \theta}{C_{pml-rml}}+\frac{1.}{(1.-R_{ml}/C_{pml})} \cdot log(R_{ml})\)
\(T = exp^{T_{s}}\)	temperature
\(p = exp^{p_{s}}\)	pressure
\(p_{v} = T_{s} + log(461.0 \cdot \rho–{v})\)	vaporpressure
\(p_{vs} = 59.260894042828049 + \frac{1885.0-4186.0}{461.0} \cdot T_{s} – \frac{3128518.15}{461.0} \cdot exp^{-T_s}\)	saturation vapor pressure
\(L_{v} = 3128518.15 + (1885.0-4186.0) \cdot exp^{T_{s}}\)	latent heat for condensation/evaporation
\(L_{vi} = 3.330e+5\)	latent heat for freezing/melting
\(S = 100.0 \cdot (exp^{p_{v}}/exp^{p_{vs}}-1.)\)	supersaturation in %
\(RelCloud = 0.1\)	relaxation parameters for smoothing
\(RelCloud2 = 1.0\)

Parameters for different distribution functions and particle properties for the classes cloud water, rain water, ice and snow.

p_cloud.nu    = 0.333333;
p_cloud.mu    = 0.666666;
p_cloud.x_max = 2.60e-10;
p_cloud.x_min = 4.20e-15; 
p_cloud.a_geo = 1.24e-01;
p_cloud.b_geo = 0.333333;
p_cloud.a_vel = 3.75e+05;
p_cloud.b_vel = 0.666667;
p_cloud.a_ven = 0.780000;
p_cloud.b_ven = 0.308000;
p_cloud.cap   = 2.0;
p_cloud.k_au  = 0.0;
p_cloud.k_sc  = 0.0;
p_cloud.k_c   = 0.0;
p_cloud.k_1   = 0.0;
p_cloud.k_2   = 0.0;

p_rain.nu    = -0.666666;
p_rain.mu    =  0.333333;
p_rain.x_max = 3.00e-06;
p_rain.x_min = 2.60e-10;
p_rain.a_geo = 1.24e-01;
p_rain.b_geo = 0.333333;
p_rain.a_vel = 1.59e+02;
p_rain.b_vel = 0.266667;
p_rain.a_ven = 0.780000;
p_rain.b_ven = 0.308000;
p_rain.cap   = 2.0;
p_rain.k_au  = 0.0;
p_rain.k_sc  = 0.0;
p_rain.k_c   = 0.0;
p_rain.k_1   = 0.0;
p_rain.k_2   = 2.0;

// HK87 'hex plates'
p_ice.nu    = -0.333333;
p_ice.mu    =  0.333333;
p_ice.x_max = 1.00e-07;
p_ice.x_min = 1.00e-12;
p_ice.a_geo = 2.17e-01;
p_ice.b_geo = 0.302115;
p_ice.a_vel = 4.19e+01;
p_ice.b_vel = 0.260000;
p_ice.a_ven = 0.860000;
p_ice.b_ven = 0.280000;
p_ice.cap   = 3.141593;
p_ice.k_au  = 0.0;
p_ice.k_sc  = 0.0;
p_ice.k_c   = 0.0;
p_ice.k_1   = 0.0;
p_ice.k_2   = 2.0; 

 //Locatelli and Hobbs
p_snow.nu    = 0.500000;
p_snow.mu    = 0.50000;
p_snow.x_max = 2.00e-06;
p_snow.x_min = 1.73e-09;
p_snow.a_geo = 8.16e-00;
p_snow.b_geo = 0.526316;
p_snow.a_vel = 2.77e+01;
p_snow.b_vel = 0.215790;
p_snow.a_ven = 0.780000;
p_snow.b_ven = 0.308000;
p_snow.cap   = 2.0;
p_snow.k_au  = 0.0;
p_snow.k_sc  = 0.0;
p_snow.k_c   = 0.0;
p_snow.k_1   = 0.0;
p_snow.k_2   = 2.0;

further parameters needed in autoconversion and self collection kernels

p_cloud.k_c  = 4.44e+9;     //..long-Kernel
p_cloud.k_1  = 4.00e+2;     //..parameter for Phi-function.
p_cloud.k_2  = 0.7e+0;      //..parameter for Phi-function.
p_cloud.k_au = 1.99047e+19; //..parameter for autoconversion.
p_cloud.k_sc = 2.05131e+10; //..parameter for selfcollection.

Warm Microphysics

Nucleation

kccn = 0.462
upper cell
Tszr =mythetazr/(cpmlzr-rmlzr)+1./(cpmlzr/rmlzr-1.)*log(rmlzr);
pszr =mythetazr/(cpmlzr-rmlzr)+1./(1.-rmlzr/cpmlzr)*log(rmlzr);
pvzr =Tszr+log(461.0*rhovzr);
pvszr=59.260894042828049+(1885.0-4186.0)/461.0*Tszr-3128518.15/461.0*exp(-Tszr);  
S2   =100.0*(exp(pvzr)/exp(pvszr)-1.0);
dSdzw=(S2-S)/dz*wLoc
calculate nucleation rate
if(S>0.0 && dSdzw>0.0 && exp(Ts)>233.15 && wLoc>0.0)
nuc=nv*kccn/S*dSdzw;
forcing=nuc;
limiter=nv;
dnc    =RelCloud2*(forcing+limiter-sqrt(forcing*forcing+limiter*limiter));
drhoc  =dnc*1e-12; // 1e-12 is smallest drop mass
drhoTh =Lv/T*drhoc+(Ts*(4186.0-1885.0)+ps*461.0)*drhoc

Condensation/Evaporation cloud water

forcing= rhov-(exp(pvs-Ts)/461.0);
limiter= texture2D(t1,pos).a;
drhoc  = RelCloud*(forcing-limiter+sqrt(forcing*forcing+limiter*limiter));
drhoTh = Lv/T*drhoc+(Ts*(4186.0-1885.0)+ps*461.0)*drhoc

Autoconversion Cloudwater to Rain

x_s=p_cloud.x_max;
k_rr=7.12; 
 if(rhoc>1e-9)
 {
  x_c=min(max(rhoc/(nc+1e-33),p_cloud.x_min),p_cloud.x_max); //..mean mass in SI
  au = p_cloud.k_au*rhoc*rhoc*x_c*x_c*rho0/(rho+1e-33);      //..autoconversion rate (SB2000)
 if(rhoc>1e-6)
 {
  tau=min(max(1.0-rhoc/(rhoc+rhor+1e-33),1e-33),1.0e0-1.0e-1);
  phi=400.*pow(tau,0.7)*pow(1.-pow(tau,0.7),3.);
  au =au*(1.0+phi/pow(1.0-tau,2.0))  } }
  
forcing= au/x_s;
limiter= nc;
dnr    = RelCloud2*(forcing+limiter-sqrt(forcing*forcing+limiter*limiter));
forcing= au; 
limiter= texture2D(t1,pos).a; //rhoc;
drhor  = RelCloud2*(forcing+limiter-sqrt(forcing*forcing+limiter*limiter))

Selfcollection cloud droplets

sc     = p_cloud.k_sc*rhoc*rhoc*rho0/(rho+1e-33);
forcing= -sc;
limiter= nc; 
dnc    = RelCloud2*(forcing-limiter+sqrt(forcing*forcing+limiter*limiter))

Selfcollection Rain

x_r=min(max(rhor/(nr+1e-33),p_rain.x_min),p_rain.x_max);
D_r=p_rain.a_geo*pow(x_r,p_rain.b_geo);
sc =k_rr*nr*rhor*sqrt(rho0/(rho+1e-33));
forcing= -sc;
limiter= nr; 
dnrsc  = RelCloud2*(forcing-limiter+sqrt(forcing*forcing+limiter*limiter))

Sedimentation Rain

a_qc=8420.0; 
 b_qc=0.8;
 rho_qc=1000.0;
 N_qc  =8.0;
 
 alf=9.65e+0;
 bet=1.03e+1;
 gam=6.00e+2;
   
 if(rhor>1.0e-10)
 {
  xr = rhor/(nr+1e-33);                                      //..mean mass in SI
  xr = min(max(xr,p_rain.x_min),p_rain.x_max);
  Dm = pow(6./(1000.0*3.14)*xr,1./3.);
  mue= (p_rain.nu+1.0e0)/p_rain.b_geo - 1.0e0;
  Dr = pow(pow(Dm,3.0)/((mue+3.)*(mue+2.)*(mue+1.)),1./3.);
  vn = alf-bet/pow(1.0+gam*Dr,mue+1.);
  vq = alf-bet/pow(1.0+gam*Dr,mue+4.);
  vn = vn*pow(rho0/(rho+1e-33),0.5);
  vq = vq*pow(rho0/(rho+1e-33),0.5);
  vn = max(vn,1.e-1);
  vq = max(vq,1.e-1);
  vn = min(vn,2.e+1);
  vq = min(vq,2.e+1);
  
  Flux   = min(rhor*vq/dz,rhor/2./dt);
  drhor -= Flux;
  drhoTh-= Ts*4186.0*Flux;
 
  Flux   = min(nr*vn/dz,nr/2./dt);
  dnr   -= Flux;
 }

Subsidence

h=float(zpos)*dz;
  //if(h<2260.)Subsidence=vsub/2260.*h; 
  Subsidence=vsub;
 
  thetar = exp((mythetazr+rmlzr*log(p0))/cpmlzr);
  thetal = exp((mythetazl+rmlzl*log(p0))/cpmlzl);
 
  Flux   = Subsidence*(thetar-thetal)/2./dz;
  drhoTh+= cpml/theta*Flux;
 
  Flux   = Subsidence*(rhorzr-rhorzl)/2./dz;
  drhor += Flux;
  drhoTh+= Ts*4186.0*Flux;
 
  Flux   = Subsidence*(rhovzr-rhovzl)/2./dz;
  drhov += Flux;
  drhoTh+= (Ts*1885.0+ps*461.0)*Flux; 
 
  Flux   = Subsidence*(rhoczr-rhoczl)/2./dz;
  drhoc += Flux;
  drhoTh+= Ts*4186.0*Flux;
 
  Flux   = Subsidence*(nczr-nczl)/2./dz;
  dnc   += Flux;
 
  Flux   = Subsidence*(nrzr-nrzl)/2./dz;
  dnr   += Flux

Evaporation Rain

K_T  = 2.500e-2;
 nu_l = 1.460e-5 ;    //..kinematic viscosity of air
 N_sc = 0.710;        //..Schmidt-Number (PK, S.541)
 n_f  = 0.333;        //..exponent of N_sc in ventilation coefficient (PK, S.541)
 m_f  = 0.500;        //..exponent of N_re in ventilation coefficient (PK, S.541)
 q_krit = 1.000e-9;       
  
 a_q = 0.429250753972163723;
 b_q = 0.180944125457962318;
 c_r = 1.0 / p_rain.cap;
 
 if(S<0.0 && rhor>0.0)
 {
  D_vtp = pow(8.7602e-5*T,1.81)/p;
  g_d   = 4.0*3.14/(Lv*Lv/(K_T*287.0*T*T)+287.0*T/(D_vtp*exp(pvs)));
 
  x_r   = min(max(rhor/(nr+1e-10),p_rain.x_min),p_rain.x_max);
  //x_r   = max(5e-5*rhor,min(rain.x_min,0.05*rhor));
 
  d_r   = p_rain.a_geo * pow(x_r,p_rain.b_geo);                 
  v_r   = p_rain.a_vel * pow(x_r,p_rain.b_vel) * pow(rho0/(rho+1e-33),0.5); 
               
  N_re  = v_r*d_r/nu_l;
  f_v   = pow(N_sc,n_f)*pow(N_re,m_f);
              
  f_q   = a_q+b_q*f_v;
             
  eva   = g_d*max(nr,1e3)*c_r*d_r*S;
  eva_q = min(-5e-9,100.*f_q*eva);
  eva_n = eva_q/x_r;
 }  
 
 forcing= eva_q;
 limiter= texture2D(t2,pos).r; //rhor;
 drhor  = RelCloud2*(forcing-limiter+sqrt(forcing*forcing+limiter*limiter));
 drhov  = -drhor; //RelCloud*(forcing+limiter-sqrt(forcing*forcing+limiter*limiter));
 
 forcing= eva_n;
 limiter= nr; //nr;
 dnr    = RelCloud2*(forcing-limiter+sqrt(forcing*forcing+limiter*limiter));
 dnv    = -dnr;
 
 drhoTh = -Lv/T*drhov-(Ts*(4186.0-1885.0)+ps*461.0)*drhov

Cold Microphysics

Freezing of cloud droplets

 if(T<273.15)
 {
  T_c=T-273.15;
  facg=2.173;
 
  if(T_c< -50.0)
  {
   fr_q= rhoc;
   fr_n= nc;
  }
  else
  {
   x_c = min(max(rhoc/(nc+1e-33),p_cloud.x_min),p_cloud.x_max);    //..mean mass
 
   //..homogeneous freezing (Jeffrey und Austin (1997), see Cotton und Field (2001) as well)
   if(T_c>-30.0) 
    j_hom = 1.0e6 * exp(-7.63-2.996*(T_c+30.0));                  //..J in 1/(m3 s) 
   else  
    j_hom = 1.0e6 * exp(-243.4-14.75*T_c-0.307*pow(T_c,2.0)-0.00287*pow(T_c,3.0)-0.0000102*pow(T_c,4.0));
 
   fr_q= j_tot*rhoc*x_c*facg;
   fr_n= j_tot*rhoc;
  }
 }
 
 forcing= fr_q;
 limiter= texture2D(t1,pos).a; //rhoc;
 drhoi  = RelCloud*(forcing+limiter-sqrt(forcing*forcing+limiter*limiter));
 
 forcing= fr_n;
 limiter= texture2D(t3,pos).r; //nc;
 dni    = RelCloud*(forcing+limiter-sqrt(forcing*forcing+limiter*limiter));
 
 drhoTh = Lvi/T*drhoi+(Ts*(2110.0-4186.0))*drhoi

Selfcollection ice and snow

ice_s_vel=0.25; // dispersion der fallges.
 snow_s_vel=0.25; // dispersion der fallges.
 delta_n=3.56212;
 theta_n=0.100138;
  
 if(rhoi>1e-9)
 {
  x_i = min(max(rhoi/(ni+1e-33),p_ice.x_min),p_ice.x_max); 
  D_i = pow(p_ice.a_geo*x_i,p_ice.b_geo);                             //..mean diameter
 
  if(T>273.15) e_coll = 1.0;                                          //..temperature dependent efficiency (Cotton et al., 1986) 
  else         e_coll = min(pow(10.0,0.035*(T-273.15)-0.7),0.2e0);    //  (see Straka, 1989; S. 53)
   
  v_i = pow(p_ice.a_vel*x_i,p_ice.b_vel)*rho0/(rho+1e-33);            //..mean sedimentationvelocity
  sn = 3.14*0.25e0*e_coll*ni*  ni*D_i*D_i*sqrt(1.0*v_i*v_i+2.0*ice_s_vel*ice_s_vel);
  sq = 3.14*0.25e0*e_coll*ni*rhoi*D_i*D_i*sqrt(1.0*v_i*v_i+2.0*ice_s_vel*ice_s_vel); 
 }
 
 forcing  = -sq;
 limiter  = texture2D(t2,pos).g; //rhoi;
 drhoi    = RelCloud*(forcing-limiter+sqrt(forcing*forcing+limiter*limiter));
 
 forcing  = -sn;
 limiter  = texture2D(t3,pos).b; //ni;
 dni      = RelCloud*(forcing-limiter+sqrt(forcing*forcing+limiter*limiter));
 
 sn=0.0;
 
 if(rhos>0.0)
 {
  if(T>273.15) e_coll = 1.0;                                        //..temperature dependent efficiency (Cotton et al., 1986) 
  else         e_coll = max(0.1e0,min(exp(0.09*(T-273.15)),1.0e0)); //  (see Straka, 1989; S. 53)
 
  x_s = min(max(rhos/(ns+1e-33),p_snow.x_min),p_snow.x_max);        //..mean mass in SI     
  D_s = pow(p_snow.a_geo*x_s,p_snow.b_geo);                         //..mean diameter
  v_s = pow(p_snow.a_vel*x_s,p_snow.b_vel)*rho0/(rho+1e-33);        //..mean sedimentationvelocity 
 
  sn = 3.14*0.125e0*e_coll*ns*ns*delta_n*D_s*D_s*sqrt(theta_n*v_s*v_s+2.0*snow_s_vel*snow_s_vel);
 }
 
 forcing  = -sn;
 limiter  = texture2D(t3,pos).a; //ns;
 dns      = RelCloud*(forcing-limiter+sqrt(forcing*forcing+limiter*limiter))

Melting of snow

a_melt_q=0.67625;
 b_melt_q=0.293159;
 
 if(T>273.15 && rhos>0.0)
 {
  e_a = 6.1078e2*exp(1.72693882e1*(T-273.15)/(T-3.586e1));
  x_s = min(max(rhos/(ns+1e-33),p_snow.x_min),p_snow.x_max);       //..mean mass in SI     
  D_s = pow(p_snow.a_geo*x_s,p_snow.b_geo);                        //..mean diameter
  v_s = pow(p_snow.a_vel*x_s,p_snow.b_vel)*rho0/(rho+1e-33);       //..mean sedimentation velocity
 
  N_re= v_s*D_s/1.460e-5;                                          //..mean Reynolds number
  fv_q= a_melt_q + b_melt_q * pow(0.710,0.333) * pow(N_re,0.500);  //..mean ventilation coefficient for water vapor
 
  D_T  = 2.5e-2/(1005.7*rho);
  fh_q = D_T/3e-5*fv_q;
 
  melt = 2.0*3.14/Lvi*D_s*ns; 
  melt_h = melt * 2.5e-2*(T-273.15);
  melt_v = melt * 3e-5*Lv/461*(e_a/T-6.10780000e2/273.15);
  melt_q = (melt_h*fh_q + melt_v*fv_q);
  melt_n = min(max((melt_q-rhos)/x_s+ns,0.0e0),ns);
  
  melt_q = min(rhos,max(melt_q,0.0e0));
  melt_n = min(ns,max(melt_n,0.0e0));
 
  if(T>283.15)
  {
   melt_q= rhos;
   melt_n= ns;
  }
 }
 
 forcing= melt_q;
 limiter= texture2D(t2,pos).b; //rhos;
 drhor  = RelCloud*(forcing+limiter-sqrt(forcing*forcing+limiter*limiter));
  
 forcing= melt_n;
 limiter= texture2D(t3,pos).a; //ns;
 dnr    = RelCloud*(forcing+limiter-sqrt(forcing*forcing+limiter*limiter));
 
 drhoTh = -Lvi/T*drhor-(Ts*(2110.0-4186.0))*drhor

snow evaporation

K_T  = 2.500e-2;
 nu_l = 1.460e-5 ;    //..Kinem. Visc. von Luft
 N_sc = 0.710;        //..Schmidt-Zahl (PK, S.541)
 n_f  = 0.333;        //..Exponent von N_sc im Vent-koeff. (PK, S.541)
 m_f  = 0.500;        //..Exponent von N_re im Vent-koeff. (PK, S.541)
 q_krit = 1.000e-9;       
  
 a_q = 0.429250753972163723;
 b_q = 0.180944125457962318;
 c_s = 1.0 / p_snow.cap;
 
 if(S<0.0 && rhos>0.0)
 {
  D_vtp = pow(8.7602e-5*T,1.81)/p;
  g_d   = 4.0*3.14/(Lv*Lv/(K_T*287.0*T*T)+287.0*T/(D_vtp*exp(pvs)));
 
  x_s   = min(max(rhos/(ns+1e-20),p_snow.x_min),p_snow.x_max);
 
  d_s   = p_snow.a_geo * pow(x_s,p_snow.b_geo);                 
  v_s   = p_snow.a_vel * pow(x_s,p_snow.b_vel) * sqrt(rho0/(rho+1e-33)); // corr-factor missing
               
  N_re  = v_s*d_s/nu_l;
  f_v   = pow(N_sc,n_f)*pow(N_re,m_f);
              
  f_q   = a_q+b_q*f_v;
        
  eva   = g_d*ns*c_s*d_s*f_v*S;
  eva_q = f_q*eva;
  eva_n = eva_q/x_s;
 }
 
 forcing= eva_q;
 limiter= texture2D(t2,pos).b; //rhos;
 drhov  = RelCloud*(forcing+limiter-sqrt(forcing*forcing+limiter*limiter));
 
 forcing= eva_n;
 limiter= texture2D(t3,pos).a; //ns;
 dnv    = RelCloud*(forcing+limiter-sqrt(forcing*forcing+limiter*limiter));
 
 drhoTh = -2.830e+6/T*drhov-(Ts*(2110.0-1885.0)+ps*461.0)*drhov

Sedimentation Snow

if(rhos>1.0e-10)
 {
  xs = rhos/ns;
  xs = min(max(xs,p_snow.x_min),p_snow.x_max);
  lam= pow(0.083333*xs,p_snow.b_vel)*sqrt(rho0/(rho+1e-33));
  vn = 42.7154*lam;
  vq = 54.1323*lam;
  vn = max(vn,0.1e+0);
  vq = max(vq,0.1e+0);
  vn = min(vn,3.0e+0);
  vq = min(vq,3.0e+0);
 
  Flux   = min(rhos*vq/dz,rhos/2./dt);
  drhos -= Flux;
  drhoTh-= Ts*2110.0*Flux;
 
  Flux   = min(ns*vn/dz,ns/2./dt);
  dns   -= Flux;
 }