178 subroutine shoc_work (ix, nx, nzm, nz, dtn, &
179 prsl, delp, phii, phil, u, v, omega, tabs, &
180 qwv, qi, qc, qpi, qpl, rhc, supice, &
181 pcrit, cefac, cesfac, tkef1, dis_opt, &
182 cld_sgs, tke, hflx, evap, prnum, tkh, &
183 wthv_sec, ntlnc, ncpl, ncpi, &
184 cp, ggr, lcond, lfus, rv, rgas, pi, epsv, eps)
190 real,
intent(in) :: cp, ggr, lcond, lfus, rv, rgas, pi, epsv, eps
191 integer,
intent(in) :: ix
192 integer,
intent(in) :: nx
194 integer,
intent(in) :: nzm
195 integer,
intent(in) :: nz
196 integer,
intent(in) :: ntlnc
197 real,
intent(in) :: dtn
199 real,
intent(in) :: pcrit
200 real,
intent(in) :: cefac
201 real,
intent(in) :: cesfac
202 real,
intent(in) :: tkef1
203 real,
intent(in) :: dis_opt
205 real,
intent(in) :: hflx(nx)
206 real,
intent(in) :: evap(nx)
210 real,
intent(in) :: prsl (ix,nzm)
211 real,
intent(in) :: delp (ix,nzm)
212 real,
intent(in) :: phii (ix,nz )
213 real,
intent(in) :: phil (ix,nzm)
214 real,
intent(in) :: u (ix,nzm)
215 real,
intent(in) :: v (ix,nzm)
216 real,
intent(in) :: omega (ix,nzm)
217 real,
intent(inout) :: tabs (ix,nzm)
218 real,
intent(inout) :: qwv (ix,nzm)
219 real,
intent(inout) :: qc (ix,nzm)
220 real,
intent(inout) :: qi (ix,nzm)
222 real,
intent(inout) :: ncpl (nx,nzm)
223 real,
intent(inout) :: ncpi (nx,nzm)
224 real,
intent(in) :: qpl (nx,nzm)
225 real,
intent(in) :: qpi (nx,nzm)
227 real,
intent(in) :: rhc (nx,nzm)
228 real,
intent(in) :: supice
229 real,
intent(out) :: cld_sgs(ix,nzm)
231 real,
intent(inout) :: tke (ix,nzm)
233 real,
intent(inout) :: tkh (ix,nzm)
234 real,
intent(in) :: prnum (nx,nzm)
235 real,
intent(inout) :: wthv_sec (ix,nzm)
237 real,
parameter :: zero=0.0_kp, one=1.0_kp, half=0.5_kp, two=2.0_kp, &
238 three=3.0_kp, oneb3=one/three, twoby3=two/three, fourb3=twoby3+twoby3
239 real,
parameter :: sqrt2 = sqrt(two), twoby15 = two / 15.0_kp, &
240 nmin = 1.0_kp, ri_cub = 6.4e-14_kp, rl_cub = 1.0e-15_kp, &
241 skew_facw=1.2_kp, skew_fact=0.0_kp, &
242 tkhmax=300.0_kp, qcmin=1.0e-9_kp
243 real :: lsub, fac_cond, fac_fus, cpolv, fac_sub, ggri, kapa, gocp, &
244 rog, sqrtpii, epsterm, onebeps, onebrvcp
248 real,
parameter :: lambda = 0.04_kp
250 real,
parameter :: min_tke = 1.0e-4_kp
252 real,
parameter :: max_tke = 40.0_kp
255 real,
parameter :: max_eddy_length_scale = 1000.0_kp
257 real,
parameter :: max_eddy_dissipation_time_scale = 2000.0_kp
258 real,
parameter :: Pr = 1.0_kp
261 real,
parameter :: pt19=0.19_kp, pt51=0.51_kp, pt01=0.01_kp, atmin=0.01_kp, atmax=one-atmin
262 real,
parameter :: Cs = 0.15_kp, epsln=1.0e-6_kp
264 real,
parameter :: Ck = 0.1_kp
274 real,
parameter :: Ce = ck**3/cs**4, ces = ce
278 real,
parameter :: vonk=0.4_kp
279 real,
parameter :: tscale=400.0_kp
280 real,
parameter :: w_tol_sqd = 4.0e-04_kp
282 real,
parameter :: w_thresh = 0.0_kp, thresh = 0.0_kp
283 real,
parameter :: w3_tol = 1.0e-20_kp
288 real,
parameter :: tbgmin = 233.16_kp
291 real,
parameter :: tbgmax = 273.16_kp
292 real,
parameter :: a_bg = one/(tbgmax-tbgmin)
297 real,
parameter :: thl2tune = 1.0_kp, qw2tune = 1.0_kp, qwthl2tune = 1.0_kp, &
299 thl_tol = 1.0e-2_kp, rt_tol = 1.0e-4_kp
301 integer,
parameter :: nitr=6
322 real thl_sec (nx,nzm)
323 real qwthl_sec(nx,nzm)
324 real wqw_sec (nx,nzm)
325 real wthl_sec (nx,nzm)
328 real wqp_sec (nx,nzm)
332 real isotropy (nx,nzm)
335 real conv_vel2(nx,nzm)
340 real diag_frac, diag_qn, diag_qi, diag_ql
361 real,
dimension(nx,nzm) :: total_water, brunt2, thv, tkesbdiss
362 real,
dimension(nx,nzm) :: def2
363 real,
dimension(nx) :: denom, numer, l_inf, cldarr, thedz, thedz2
365 real lstarn, depth, omn, betdz, bbb, term, qsatt, dqsat, &
366 conv_var, tkes, skew_w, skew_qw, aterm, w1_1, w1_2, w2_1, &
367 w2_2, w3var, thl1_1, thl1_2, thl2_1, thl2_2, qw1_1, qw1_2, qw2_1, &
368 qw2_2, ql1, ql2, w_ql1, w_ql2, &
369 r_qwthl_1, r_wqw_1, r_wthl_1, testvar, s1, s2, std_s1, std_s2, C1, C2, &
370 thl_first, qw_first, w_first, Tl1_1, Tl1_2, betatest, pval, pkap, &
371 w2thl, w2qw,w2ql, w2ql_1, w2ql_2, &
372 thec, thlsec, qwsec, qwthlsec, wqwsec, wthlsec, thestd,dum, &
373 cqt1, cthl1, cqt2, cthl2, qn1, qn2, qi1, qi2, omn1, omn2, &
374 basetemp2, beta1, beta2, qs1, qs2, &
375 esval, esval2, om1, om2, epss, &
376 lstarn1, lstarn2, sqrtw2, sqrtthl, sqrtqt, &
377 sqrtstd1, sqrtstd2, tsign, tvar, sqrtw2t, wqls, wqis, &
378 sqrtqw2_1, sqrtqw2_2, sqrtthl2_1, sqrtthl2_2, sm, prespot, &
379 corrtest1, corrtest2, wrk, wrk1, wrk2, wrk3, onema, pfac
382 integer i,k,km1,ku,kd,ka,kb
395 sqrtpii = one/sqrt(pi+pi)
397 onebeps = one/epsterm
398 onebrvcp = one/(rv*cp)
406 zi(i,k) = phii(i,k) * ggri
415 if (qc(i,k) < zero)
then
416 qwv(i,k) = qwv(i,k) + qc(i,k)
417 tabs(i,k) = tabs(i,k) - fac_cond * qc(i,k)
420 if (qi(i,k) < zero)
then
421 qwv(i,k) = qwv(i,k) + qi(i,k)
422 tabs(i,k) = tabs(i,k) - fac_sub * qi(i,k)
448 if (qwv(i,k) < zero)
then
449 qwv(i,k) = qwv(i,km1) + qwv(i,k) * delp(i,k) / delp(i,km1)
456 zl(i,k) = phil(i,k) * ggri
457 wrk = one / prsl(i,k)
458 qv(i,k) = max(qwv(i,k), zero)
459 thv(i,k) = tabs(i,k) * (one+epsv*qv(i,k))
460 w(i,k) = - rog * omega(i,k) * thv(i,k) * wrk
461 qcl(i,k) = max(qc(i,k), zero)
462 qci(i,k) = max(qi(i,k), zero)
469 total_water(i,k) = qcl(i,k) + qci(i,k) + qv(i,k)
471 prespot = (100000.0_kp*wrk) ** kapa
472 bet(i,k) = ggr/(tabs(i,k)*prespot)
473 thv(i,k) = thv(i,k)*prespot
476 gamaz(i,k) = gocp * zl(i,k)
479 hl(i,k) = tabs(i,k) + gamaz(i,k) - fac_cond*(qcl(i,k)+qpl(i,k)) &
480 - fac_sub *(qci(i,k)+qpi(i,k))
490 adzi(i,k) = zl(i,k) - zl(i,km1)
491 adzl(i,km1) = zi(i,k) - zi(i,km1)
495 adzi(i,1) = (zl(i,1)-zi(i,1))
496 adzi(i,nz) = adzi(i,nzm)
497 adzl(i,nzm) = zi(i,nz) - zi(i,nzm)
499 wthl_sec(i,1) = hflx(i)
500 wqw_sec(i,1) = evap(i)
525 elseif (k == nzm)
then
530 if (tke(i,k) > zero)
then
532 wrk = half*(tkh(i,ka)*prnum(i,ka)+tkh(i,kb)*prnum(i,kb))*(w(i,ku) - w(i,kd)) &
533 * sqrt(tke(i,k)) / (zl(i,ku) - zl(i,kd))
534 w_sec(i,k) = max(twoby3 * tke(i,k) - twoby15 * wrk, zero)
550 wrk1 = one / adzi(i,k)
552 wrk3 = max(tkh(i,k),epsln) * wrk1
554 sm = half*(isotropy(i,k)+isotropy(i,km1))*wrk1*wrk3
559 wrk1 = hl(i,k) - hl(i,km1) &
560 + (qpl(i,k) - qpl(i,km1)) * fac_cond &
561 + (qpi(i,k) - qpi(i,km1)) * fac_sub
562 wthl_sec(i,k) = - wrk3 * wrk1
566 wrk2 = total_water(i,k) - total_water(i,km1)
567 wqw_sec(i,k) = - wrk3 * wrk2
571 thl_sec(i,k) = thl2tune * sm * wrk1 * wrk1
575 qw_sec(i,k) = qw2tune * sm * wrk2 * wrk2
580 qwthl_sec(i,k) = qwthl2tune * sm * wrk1 * wrk2
589 thl_sec(i,1) = thl_sec(i,2)
590 qw_sec(i,1) = qw_sec(i,2)
591 qwthl_sec(i,1) = qwthl_sec(i,2)
606 subroutine tke_shoc()
611 real grd,betdz,Cee,lstarn, lstarp, bbb, omn, omp,qsatt,dqsat, smix, &
612 buoy_sgs,ratio,a_prod_sh,a_prod_bu,a_diss,a_prod_bu_debug, buoy_sgs_debug, &
613 tscale1, wrk, wrk1, wtke, wtk2, rdtn, tkef2
614 integer i,k,ku,kd,itr,k1
618 call tke_shear_prod(def2)
624 tke(i,k) = max(min_tke,tke(i,k))
625 tkesbdiss(i,k) = zero
645 elseif(k == nzm)
then
651 if (dis_opt > 0)
then
653 wrk = (zl(i,k)-zi(i,1)) / adzl(i,1) + 1.5_kp
654 cek(i) = (one + two / max((wrk*wrk - 3.3_kp), 0.5_kp)) * cefac
671 a_prod_bu = ggr / thv(i,k) * wthv_sec(i,k)
678 buoy_sgs = - (a_prod_bu+a_prod_bu) / (tkh(i,ku)+tkh(i,kd) + 0.0001_kp)
683 if (buoy_sgs <= zero)
then
686 smix = min(grd,max(0.1_kp*grd, 0.76_kp*sqrt(tke(i,k)/(buoy_sgs+1.0e-10_kp))))
690 cee = cek(i) * (pt19 + pt51*ratio) * max(one, sqrt(pcrit/prsl(i,k)))
693 a_prod_sh = half*(def2(i,ku)*tkh(i,ku)*prnum(i,ku) &
694 + def2(i,kd)*tkh(i,kd)*prnum(i,kd))
707 wrk = (dtn*cee) / smixt(i,k)
708 wrk1 = wtke + dtn*(a_prod_sh+a_prod_bu)
711 wtke = min(max(min_tke, wtke), max_tke)
712 a_diss = wrk*sqrt(wtke)
713 wtke = wrk1 / (one+a_diss)
714 wtke = tkef1*wtke + tkef2*wtk2
718 tke(i,k) = min(max(min_tke, wtke), max_tke)
719 a_diss = wrk*sqrt(tke(i,k))
721 tscale1 = (dtn+dtn) / a_diss
723 tkesbdiss(i,k) = rdtn*a_diss*tke(i,k)
728 if (buoy_sgs <= zero)
then
729 isotropy(i,k) = min(max_eddy_dissipation_time_scale, tscale1)
731 isotropy(i,k) = min(max_eddy_dissipation_time_scale, &
732 tscale1/(one+lambda*buoy_sgs*tscale1*tscale1))
749 tkh(i,k) = min(tkhmax, wrk * (isotropy(i,k) * tke(i,k) &
750 + isotropy(i,k1) * tke(i,k1)))
755 end subroutine tke_shoc
758 subroutine tke_shear_prod(def2)
762 real,
intent(out) :: def2(nx,nzm)
764 real rdzw, wrku, wrkv, wrkw
772 rdzw = one / adzi(i,k)
773 wrku = (u(i,k)-u(i,k1)) * rdzw
774 wrkv = (v(i,k)-v(i,k1)) * rdzw
776 def2(i,k) = wrku*wrku + wrkv*wrkv
781 def2(i,1) = (u(i,1)*u(i,1) + v(i,1)*v(i,1)) / (zl(i,1)*zl(i,1))
784 end subroutine tke_shear_prod
786 subroutine eddy_length()
792 real wrk, wrk1, wrk2, wrk3
793 integer i, k, kk, kl, ku, kb, kc, kli, kui
815 if (qcl(i,k)+qci(i,k) <= qcmin)
then
816 tkes = sqrt(tke(i,k)) * adzl(i,k)
817 numer(i) = numer(i) + tkes*zl(i,k)
818 denom(i) = denom(i) + tkes
827 if (denom(i) > zero .and. numer(i) > zero)
then
828 l_inf(i) = min(0.1_kp * (numer(i)/denom(i)), 100.0_kp)
842 thedz(:) = adzi(:,kc)
843 elseif (k == nzm)
then
848 thedz(:) = adzi(:,kc) + adzi(:,k)
855 betdz = bet(i,k) / thedz(i)
857 tkes = sqrt(tke(i,k))
861 wrk = qcl(i,k) + qci(i,k)
866 omn = qcl(i,k) / (wrk+1.0e-20_kp)
871 lstarn = fac_cond + (one-omn)*fac_fus
874 dqsat = omn * dtqsatw(tabs(i,k),prsl(i,k)) &
875 + (one-omn) * dtqsati(tabs(i,k),prsl(i,k))
879 qsatt = omn * qsatw(tabs(i,k),prsl(i,k)) &
880 + (one-omn) * qsati(tabs(i,k),prsl(i,k))
884 bbb = (one + epsv*qsatt-wrk-qpl(i,k)-qpi(i,k) &
885 + 1.61_kp*tabs(i,k)*dqsat) / (one+lstarn*dqsat)
889 brunt(i,k) = betdz*(bbb*(hl(i,kc)-hl(i,kb)) &
890 + (bbb*lstarn - (one+lstarn*dqsat)*tabs(i,k)) &
891 * (total_water(i,kc)-total_water(i,kb)) &
892 + (bbb*fac_cond - (one+fac_cond*dqsat)*tabs(i,k))*(qpl(i,kc)-qpl(i,kb)) &
893 + (bbb*fac_sub - (one+fac_sub*dqsat)*tabs(i,k))*(qpi(i,kc)-qpi(i,kb)) )
901 bbb = one + epsv*qv(i,k) - qpl(i,k) - qpi(i,k)
902 brunt(i,k) = betdz*( bbb*(hl(i,kc)-hl(i,kb)) &
903 + epsv*tabs(i,k)*(total_water(i,kc)-total_water(i,kb)) &
904 + (bbb*fac_cond-tabs(i,k))*(qpl(i,kc)-qpl(i,kb)) &
905 + (bbb*fac_sub -tabs(i,k))*(qpi(i,kc)-qpi(i,kb)) )
911 if (brunt(i,k) >= zero)
then
912 brunt2(i,k) = brunt(i,k)
931 if (tkes > zero .and. l_inf(i) > zero)
then
932 wrk1 = one / (tscale*tkes*vonk*zl(i,k))
933 wrk2 = one / (tscale*tkes*l_inf(i))
934 wrk1 = wrk1 + wrk2 + pt01 * brunt2(i,k) / tke(i,k)
935 wrk1 = sqrt(one / max(wrk1,1.0e-8_kp)) * (one/0.3_kp)
937 smixt(i,k) = min(max_eddy_length_scale, wrk1)
976 if (cldarr(i) == 1)
then
984 wrk = qcl(i,k) + qci(i,k)
985 if (wrk > qcmin)
then
991 if (qcl(i,k+1)+qci(i,k+1) <= qcmin)
then
1002 if (kl > 0 .and. ku >= kl)
then
1006 conv_var = conv_var+ 2.5_kp*adzi(i,kk)*bet(i,kk)*wthv_sec(i,kk)
1008 conv_var = conv_var ** oneb3
1010 if (conv_var > zero)
then
1012 depth = (zl(i,ku)-zl(i,kl)) + adzl(i,kl)
1020 wrk = conv_var/(depth*depth*sqrt(tke(i,kk))) &
1021 + pt01*brunt2(i,kk)/tke(i,kk)
1023 smixt(i,kk) = min(max_eddy_length_scale, (one/0.3_kp)*sqrt(one/wrk))
1037 end subroutine eddy_length
1040 subroutine conv_scale()
1055 conv_vel2(i,1) = zero
1069 conv_vel2(i,k) = conv_vel2(i,k-1) &
1070 + 2.5_kp*adzi(i,k)*bet(i,k)*wthv_sec(i,k)
1074 end subroutine conv_scale
1077 subroutine check_eddy()
1081 integer i, k, kb, ks, zend
1101 wrk = 0.1_kp*adzl(i,k)
1103 smixt(i,k) = max(wrk, min(max_eddy_length_scale,smixt(i,k)))
1109 if (qcl(i,kb) == zero .and. qcl(i,k) > zero .and. brunt(i,k) > 1.0e-4_kp)
then
1117 end subroutine check_eddy
1127 integer i, k, kb, kc
1129 real bet2, f0, f1, f2, f3, f4, f5, iso, isosqr, &
1130 omega0, omega1, omega2, X0, Y0, X1, Y1, AA0, AA1, buoy_sgs2, &
1131 wrk, wrk1, wrk2, wrk3, avew
1135 real,
parameter :: c=7.0_kp, a0=0.52_kp/(c*c*(c-2.0_kp)), a1=0.87_kp/(c*c), &
1136 a2=0.5_kp/c, a3=0.6_kp/(c*(c-2.0_kp)), a4=2.4_kp/(3.0_kp*c+5.0_kp), &
1137 a5=0.6_kp/(c*(3.0_kp*c+5.0_kp))
1158 thedz(i) = one / adzi(i,k)
1159 thedz2(i) = one / adzl(i,kb)
1163 thedz(i) = one / adzi(i,k)
1164 thedz2(i) = one / (adzl(i,k)+adzl(i,kb))
1170 iso = half*(isotropy(i,k)+isotropy(i,kb))
1172 buoy_sgs2 = isosqr*half*(brunt(i,k)+brunt(i,kb))
1173 bet2 = half*(bet(i,k)+bet(i,kb))
1179 avew = half*(w_sec(i,k)+w_sec(i,kb))
1203 wrk2 = thedz2(i)*wrk1*wrk1*iso
1204 wrk3 = thl_sec(i,kc) - thl_sec(i,kb)
1206 f0 = wrk2 * wrk1 * wthl_sec(i,k) * wrk3
1208 wrk = wthl_sec(i,kc) - wthl_sec(i,kb)
1210 f1 = wrk2 * (wrk*wthl_sec(i,k) + half*avew*wrk3)
1213 f2 = thedz(i)*wrk1*wthl_sec(i,k)*(w_sec(i,k)-w_sec(i,kb)) &
1214 + (thedz2(i)+thedz2(i))*bet(i,k)*isosqr*wrk
1216 f3 = thedz2(i)*wrk1*wrk + thedz(i)*bet2*isosqr*(wthl_sec(i,k)*(tke(i,k)-tke(i,kb)))
1218 wrk1 = thedz(i)*iso*avew
1219 f4 = wrk1*(w_sec(i,k)-w_sec(i,kb) + tke(i,k)-tke(i,kb))
1221 f5 = wrk1*(w_sec(i,k)-w_sec(i,kb))
1226 omega0 = a4 / (one-a5*buoy_sgs2)
1227 omega1 = omega0 / (c+c)
1228 omega2 = omega1*f3+(5.0_kp/4.0_kp)*omega0*f4
1232 wrk1 = one / (one-(a1+a3)*buoy_sgs2)
1233 wrk2 = one / (one-a3*buoy_sgs2)
1234 x0 = wrk1 * (a2*buoy_sgs2*(one-a3*buoy_sgs2))
1235 y0 = wrk2 * (two*a2*buoy_sgs2*x0)
1236 x1 = wrk1 * (a0*f0+a1*f1+a2*(one-a3*buoy_sgs2)*f2)
1237 y1 = wrk2 * (two*a2*(buoy_sgs2*x1+(a0/a1)*f0+f1))
1241 aa0 = omega0*x0 + omega1*y0
1242 aa1 = omega0*x1 + omega1*y1 + omega2
1251 w3(i,k) = (aa1-1.2_kp*x1-1.5_kp*f5)/(c-1.2_kp*x0+aa0)
1264 end subroutine canuto
1266 subroutine assumed_pdf()
1279 real wrk, wrk1, wrk2, wrk3, wrk4, bastoeps, eps_ss1, eps_ss2, cond_w
1305 pfac = pval * 1.0e-5_kp
1310 thl_first = hl(i,k) + fac_cond*qpl(i,k) + fac_sub*qpi(i,k)
1311 qw_first = total_water(i,k)
1321 w3var = half*(w3(i,kd)+w3(i,ku))
1322 thlsec = max(zero, half*(thl_sec(i,kd)+thl_sec(i,ku)) )
1323 qwsec = max(zero, half*(qw_sec(i,kd)+qw_sec(i,ku)) )
1324 qwthlsec = half * (qwthl_sec(i,kd) + qwthl_sec(i,ku))
1325 wqwsec = half * (wqw_sec(i,kd) + wqw_sec(i,ku))
1326 wthlsec = half * (wthl_sec(i,kd) + wthl_sec(i,ku))
1328 w3var = half*w3(i,k)
1329 thlsec = max(zero, half*thl_sec(i,k))
1330 qwsec = max(zero, half*qw_sec(i,k))
1331 qwthlsec = half * qwthl_sec(i,k)
1332 wqwsec = half * wqw_sec(i,k)
1333 wthlsec = half * wthl_sec(i,k)
1344 if (w_sec(i,k) > zero)
then
1345 sqrtw2 = sqrt(w_sec(i,k))
1349 if (thlsec > zero)
then
1350 sqrtthl = sqrt(thlsec)
1354 if (qwsec > zero)
then
1355 sqrtqt = sqrt(qwsec)
1367 IF (w_sec(i,k) <= w_tol_sqd)
THEN
1379 cond_w = 1.2_kp*sqrt2*max(w3_tol, sqrtw2*sqrtw2*sqrtw2)
1380 w3var = max(-cond_w, min(cond_w, w3var))
1383 skew_w = w3var / (sqrtw2*sqrtw2*sqrtw2)
1393 aterm = max(atmin,min(half*(one-skew_w*sqrt(one/(4.0_kp*wrk*wrk*wrk+skew_w*skew_w))),atmax))
1399 wrk = sqrt(onema/aterm)
1400 w1_1 = sqrtw2t * wrk
1401 w1_2 = - sqrtw2t / wrk
1403 w2_1 = w2_1 * w_sec(i,k)
1404 w2_2 = w2_2 * w_sec(i,k)
1410 IF (thlsec <= thl_tol*thl_tol .or. abs(w1_2-w1_1) <= w_thresh)
THEN
1419 corrtest1 = max(-one,min(one,wthlsec/(sqrtw2*sqrtthl)))
1421 thl1_1 = -corrtest1 / w1_2
1422 thl1_2 = -corrtest1 / w1_1
1424 wrk1 = thl1_1 * thl1_1
1425 wrk2 = thl1_2 * thl1_2
1426 wrk3 = three * (one - aterm*wrk1 - onema*wrk2)
1427 wrk4 = -skew_facw*skew_w - aterm*wrk1*thl1_1 - onema*wrk2*thl1_2
1430 wrk = three * (thl1_2-thl1_1)
1431 if (wrk /= zero)
then
1432 thl2_1 = thlsec * min(100.0_kp,max(zero,(thl1_2*wrk3-wrk4)/(aterm*wrk)))
1433 thl2_2 = thlsec * min(100.0_kp,max(zero,(-thl1_1*wrk3+wrk4)/(onema*wrk)))
1439 thl1_1 = thl1_1*sqrtthl + thl_first
1440 thl1_2 = thl1_2*sqrtthl + thl_first
1442 sqrtthl2_1 = sqrt(thl2_1)
1443 sqrtthl2_2 = sqrt(thl2_2)
1449 IF (qwsec <= rt_tol*rt_tol .or. abs(w1_2-w1_1) <= w_thresh)
THEN
1458 corrtest2 = max(-one,min(one,wqwsec/(sqrtw2*sqrtqt)))
1460 qw1_1 = - corrtest2 / w1_2
1461 qw1_2 = - corrtest2 / w1_1
1463 tsign = abs(qw1_2-qw1_1)
1467 IF (tsign > 0.4_kp)
THEN
1468 skew_qw = skew_facw*skew_w
1469 ELSEIF (tsign <= 0.2_kp)
THEN
1472 skew_qw = (skew_facw/0.2_kp) * skew_w * (tsign-0.2_kp)
1475 wrk1 = qw1_1 * qw1_1
1476 wrk2 = qw1_2 * qw1_2
1477 wrk3 = three * (one - aterm*wrk1 - onema*wrk2)
1478 wrk4 = skew_qw - aterm*wrk1*qw1_1 - onema*wrk2*qw1_2
1479 wrk = three * (qw1_2-qw1_1)
1481 if (wrk /= zero)
then
1482 qw2_1 = qwsec * min(100.0_kp,max(zero,( qw1_2*wrk3-wrk4)/(aterm*wrk)))
1483 qw2_2 = qwsec * min(100.0_kp,max(zero,(-qw1_1*wrk3+wrk4)/(onema*wrk)))
1489 qw1_1 = qw1_1*sqrtqt + qw_first
1490 qw1_2 = qw1_2*sqrtqt + qw_first
1492 sqrtqw2_1 = sqrt(qw2_1)
1493 sqrtqw2_2 = sqrt(qw2_2)
1499 w1_1 = w1_1*sqrtw2 + w_first
1500 w1_2 = w1_2*sqrtw2 + w_first
1504 testvar = aterm*sqrtqw2_1*sqrtthl2_1 + onema*sqrtqw2_2*sqrtthl2_2
1506 IF (testvar == zero)
THEN
1509 r_qwthl_1 = max(-one,min(one,(qwthlsec-aterm*(qw1_1-qw_first)*(thl1_1-thl_first) &
1510 -onema*(qw1_2-qw_first)*(thl1_2-thl_first))/testvar))
1519 tl1_1 = thl1_1 - gamaz(i,k)
1520 tl1_2 = thl1_2 - gamaz(i,k)
1526 IF (tl1_1 >= tbgmax)
THEN
1528 esval = min(
fpvsl(tl1_1), pval)
1529 qs1 = eps * esval / (pval-0.378_kp*esval)
1530 ELSE IF (tl1_1 <= tbgmin)
THEN
1532 esval = min(
fpvsi(tl1_1), pval)
1533 qs1 = epss * esval / (pval-0.378_kp*esval)
1535 om1 = max(zero, min(one, a_bg*(tl1_1-tbgmin)))
1536 lstarn1 = lcond + (one-om1)*lfus
1537 esval = min(
fpvsl(tl1_1), pval)
1538 esval2 = min(
fpvsi(tl1_1), pval)
1539 qs1 = om1 * eps * esval / (pval-0.378_kp*esval) &
1540 + (one-om1) * epss * esval2 / (pval-0.378_kp*esval2)
1546 beta1 = lstarn1 / tl1_1
1547 beta1 = beta1 * beta1 * onebrvcp
1552 IF (tl1_1 == tl1_2)
THEN
1556 IF (tl1_2 >= tbgmax)
THEN
1558 esval = min(
fpvsl(tl1_2), pval)
1559 qs2 = eps * esval / (pval-0.378_kp*esval)
1560 ELSE IF (tl1_2 <= tbgmin)
THEN
1562 esval = min(
fpvsi(tl1_2), pval)
1563 qs2 = epss * esval / (pval-0.378_kp*esval)
1565 om2 = max(zero, min(one, a_bg*(tl1_2-tbgmin)))
1566 lstarn2 = lcond + (one-om2)*lfus
1567 esval = min(
fpvsl(tl1_2), pval)
1568 esval2 = min(
fpvsi(tl1_2), pval)
1569 qs2 = om2 * eps * esval / (pval-0.378_kp*esval) &
1570 + (one-om2) * epss * esval2 / (pval-0.378_kp*esval2)
1576 beta2 = lstarn2 / tl1_2
1577 beta2 = beta2 * beta2 * onebrvcp
1582 qs1 = qs1 * rhc(i,k)
1583 qs2 = qs2 * rhc(i,k)
1587 cqt1 = one / (one+beta1*qs1)
1588 wrk = qs1 * (one+beta1*qw1_1) * cqt1
1590 cthl1 = cqt1*wrk*cpolv*beta1*pkap
1592 wrk1 = cthl1 * cthl1
1595 std_s1 = sqrt(max(zero, wrk1*thl2_1+wrk2*qw2_1 &
1596 - two*cthl1*sqrtthl2_1*cqt1*sqrtqw2_1*r_qwthl_1))
1601 IF (std_s1 > zero)
THEN
1602 wrk = s1 / (std_s1*sqrt2)
1603 c1 = max(zero, min(one, half*(one+erf(wrk))))
1605 IF (c1 > zero) qn1 = s1*c1 + (std_s1*sqrtpii)*exp(-wrk*wrk)
1606 ELSEIF (s1 >= qcmin)
THEN
1615 IF (qw1_1 == qw1_2 .and. thl2_1 == thl2_2 .and. qs1 == qs2)
THEN
1624 cqt2 = one / (one+beta2*qs2)
1625 wrk = qs2 * (one+beta2*qw1_2) * cqt2
1627 cthl2 = wrk*cqt2*cpolv*beta2*pkap
1628 wrk1 = cthl2 * cthl2
1631 std_s2 = sqrt(max(zero, wrk1*thl2_2+wrk2*qw2_2 &
1632 - two*cthl2*sqrtthl2_2*cqt2*sqrtqw2_2*r_qwthl_1))
1637 IF (std_s2 > zero)
THEN
1638 wrk = s2 / (std_s2*sqrt2)
1639 c2 = max(zero, min(one, half*(one+erf(wrk))))
1640 IF (c2 > zero) qn2 = s2*c2 + (std_s2*sqrtpii)*exp(-wrk*wrk)
1641 ELSEIF (s2 >= qcmin)
THEN
1649 diag_frac = aterm*c1 + onema*c2
1651 om1 = max(zero, min(one, (tl1_1-tbgmin)*a_bg))
1652 om2 = max(zero, min(one, (tl1_2-tbgmin)*a_bg))
1654 qn1 = min(qn1,qw1_1)
1655 qn2 = min(qn2,qw1_2)
1663 diag_qn = min(max(zero, aterm*qn1 + onema*qn2), total_water(i,k))
1664 diag_ql = min(max(zero, aterm*ql1 + onema*ql2), diag_qn)
1665 diag_qi = max(zero, diag_qn - diag_ql)
1668 om1 = max(zero, min(one, (tabs(i,k)-tbgmin)*a_bg))
1669 lstarn1 = lcond + (one-om1)*lfus
1670 tabs(i,k) = hl(i,k) - gamaz(i,k) + fac_cond*(diag_ql+qpl(i,k)) &
1671 + fac_sub *(diag_qi+qpi(i,k)) &
1672 + tkesbdiss(i,k) * (dtn/cp)
1679 if (ncpl(i,k) > nmin)
then
1680 ncpl(i,k) = diag_ql/max(qc(i,k),1.0e-10_kp)*ncpl(i,k)
1682 ncpl(i,k) = max(diag_ql/(fourb3*pi*rl_cub*997.0_kp), nmin)
1684 if (ncpi(i,k) > nmin)
then
1685 ncpi(i,k) = diag_qi/max(qi(i,k),1.0e-10_kp)*ncpi(i,k)
1687 ncpi(i,k) = max(diag_qi/(fourb3*pi*ri_cub*500.0_kp), nmin)
1694 qwv(i,k) = max(zero, total_water(i,k) - diag_qn)
1695 cld_sgs(i,k) = diag_frac
1698 wqls = aterm * ((w1_1-w_first)*ql1) + onema * ((w1_2-w_first)*ql2)
1699 wqis = aterm * ((w1_1-w_first)*qi1) + onema * ((w1_2-w_first)*qi2)
1702 wqlsb(k) = wqlsb(k) + wqls
1703 wqisb(k) = wqisb(k) + wqis
1708 wrk = epsv * thv(i,k)
1710 bastoeps = onebeps * thv(i,k)
1713 wthv_sec(i,k) = wthlsec + wrk*wqwsec &
1714 + (fac_cond-bastoeps)*wqls &
1715 + (fac_sub-bastoeps) *wqis &
1716 + ((lstarn1/cp)-thv(i,k))*half*(wqp_sec(i,kd)+wqp_sec(i,ku))
1718 wthv_sec(i,k) = wthlsec + wrk*wqwsec &
1719 + (fac_cond-bastoeps)*wqls &
1720 + (fac_sub-bastoeps) *wqis &
1721 + ((lstarn1/cp)-thv(i,k))*half*wqp_sec(i,k)
1733 end subroutine assumed_pdf
1741 real function esatw(t)
1743 real a0,a1,a2,a3,a4,a5,a6,a7,a8
1744 data a0,a1,a2,a3,a4,a5,a6,a7,a8 / &
1745 6.11239921, 0.443987641, 0.142986287e-1, &
1746 0.264847430e-3, 0.302950461e-5, 0.206739458e-7, &
1747 0.640689451e-10, -0.952447341e-13,-0.976195544e-15/
1749 dt = max(-80.,t-273.16)
1750 esatw = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
1753 real function qsatw(t,p)
1760 qsatw = 0.622 * esat/max(esat,p-0.378*esat)
1766 real function esati(t)
1768 real a0,a1,a2,a3,a4,a5,a6,a7,a8
1769 data a0,a1,a2,a3,a4,a5,a6,a7,a8 / &
1770 6.11147274, 0.503160820, 0.188439774e-1, &
1771 0.420895665e-3, 0.615021634e-5, 0.602588177e-7, &
1772 0.385852041e-9, 0.146898966e-11, 0.252751365e-14/
1777 else if(t.gt.185.)
then
1779 esati = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
1781 dt = max(-100.,t-273.16)
1782 esati = 0.00763685 + dt*(0.000151069+dt*7.48215e-07)
1786 real function qsati(t,p)
1797 real function dtesatw(t)
1799 real a0,a1,a2,a3,a4,a5,a6,a7,a8
1800 data a0,a1,a2,a3,a4,a5,a6,a7,a8 / &
1801 0.443956472, 0.285976452e-1, 0.794747212e-3, &
1802 0.121167162e-4, 0.103167413e-6, 0.385208005e-9, &
1803 -0.604119582e-12, -0.792933209e-14, -0.599634321e-17/
1805 dt = max(-80.,t-273.16)
1806 dtesatw = a0 + dt* (a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
1807 end function dtesatw
1809 real function dtqsatw(t,p)
1813 dtqsatw = 100.0*0.622*dtesatw(t)/p
1814 end function dtqsatw
1816 real function dtesati(t)
1818 real a0,a1,a2,a3,a4,a5,a6,a7,a8
1819 data a0,a1,a2,a3,a4,a5,a6,a7,a8 / &
1820 0.503223089, 0.377174432e-1, 0.126710138e-2, &
1821 0.249065913e-4, 0.312668753e-6, 0.255653718e-8, &
1822 0.132073448e-10, 0.390204672e-13, 0.497275778e-16/
1826 dtesati = dtesatw(t)
1827 else if(t > 185.)
then
1829 dtesati = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
1831 dt = max(-100.,t-273.16)
1832 dtesati = 0.0013186 + dt*(2.60269e-05+dt*1.28676e-07)
1834 end function dtesati
1837 real function dtqsati(t,p)
1841 dtqsati = 100.0*0.622*dtesati(t)/p
1842 end function dtqsati