# 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. 

# 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 

# 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; }