FV3 Bundle
cloud.F90
Go to the documentation of this file.
1 MODULE cloud
2 
3 IMPLICIT NONE
4 
5 PRIVATE
6 
7 !Make subroutines public for adjoint checkpointing
8 PUBLIC cloud_driver, cloud_tidy, meltfreeze, get_ice_fraction, convec_src, pdf_width, ls_cloud, &
9  autoconversion_ls, autoconversion_cnv, ice_settlefall_ls, ice_settlefall_cnv, precipandevap, &
10  pdffrac, pdfcondensate, cons_sundq3, cons_alhx, cons_microphys, marshpalm, evap_cnv, subl_cnv, ldradius, &
11  dqsat_bac, dqsats_bac
12 
13 CONTAINS
14 
15 subroutine cloud_driver ( DT, IM, JM, LM, th, q, ple, &
16  CNV_DQLDT, CNV_MFD, CNV_PRC3, CNV_UPDF, &
17  QI_ls, QL_ls, QI_con, QL_con, CF_LS, CF_con, &
18  FRLAND, PHYSPARAMS, ESTBLX, KHu, KHl, &
19  CONS_RUNIV, CONS_KAPPA, CONS_AIRMW, &
20  CONS_H2OMW, CONS_GRAV, CONS_ALHL, &
21  CONS_ALHF, CONS_PI, CONS_RGAS, &
22  CONS_CP, CONS_VIREPS, CONS_ALHS, &
23  CONS_TICE, CONS_RVAP, CONS_P00, do_moist_physics )
24 
25 !NAME LIST
26 !DT - Physics timestep
27 !IM - Number of points in x-dir (per processor)
28 !JM - Number of points in y-dir
29 !LM - Number of vertical levels
30 !th (d)- Potential temperature (K), defined with p0 = 1000hPa (th=T at the surface)
31 !q (d)- Specific humidity (kg kg^-1)
32 !ple - Pressure at level edge (Pa)
33 !CNV_DQLDT (d)- Convective condensate source from RAS (kg m^-2 s^-1)
34 !CNV_MFD (d)- Convective detraining mass flux from RAS (kg m^-2 s^-1)
35 !CNV_PRC3 (d)- Convective precipitation from RAS (kg m^-2 s^-1)
36 !CNV_UPDF (d)- Convective updraft areal fraction from RAS (1)
37 !QI_ls (d)- Mass fraction of large scale cloud ice water (kg kg^-1)
38 !QL_ls (d)- Mass fraction of large scale cloud liquid water (kg kg^-1)
39 !QI_con (d)- Mass fraction of convective cloud ice water (kg kg^-1)
40 !QL_con (d)- Mass fraction of convective cloud liquid water (kg kg^-1)
41 !CF_LS (d)- Large scale cloud area fraction (1)
42 !CF_con (d)- Convective cloud area fraction (1)
43 !FRLAND - Fraction of land, 0 to 1
44 !PHYSPARAMS - Set of constants
45 !ESTBLX - Table of values for calculating qsat and dqsatdT
46 
47 !(d) = TLM Dependent variable
48 
49 IMPLICIT NONE
50 
51 !INPUTS
52 integer, intent(in) :: im, jm, lm, do_moist_physics
53 real(8), intent(in) :: dt, frland(im,jm), physparams(:)
54 real(8), intent(in), dimension(IM,JM,LM) :: cnv_dqldt, cnv_mfd, cnv_updf, cnv_prc3
55 
56 real(8), intent(in) :: estblx(:)
57 
58 real(8), intent(in), dimension(IM,JM,0:LM) :: ple
59 integer, intent(in), dimension(IM,JM) :: khu, khl
60 
61 !MAPL_CONSTANTS REDEFINED FOR USE IN AUTODIFF TOOL
62 real(8), intent(in) :: cons_runiv, cons_kappa, cons_airmw
63 real(8), intent(in) :: cons_h2omw, cons_grav, cons_alhl
64 real(8), intent(in) :: cons_alhf, cons_pi, cons_rgas
65 real(8), intent(in) :: cons_cp, cons_vireps, cons_alhs
66 real(8), intent(in) :: cons_tice, cons_rvap, cons_p00
67 
68 !PROGNOSTICS
69 real(8), intent(inout), dimension(IM,JM,LM) :: th, q
70 real(8), intent(inout), dimension(IM,JM,LM) :: qi_ls, ql_ls, qi_con, ql_con, cf_con, cf_ls
71 
72 !OUTPUTS (DIAGNOSTICS)
73 
74 !LOCALS
75 integer :: i, j, k, l, ktop
76 
77 real(8), dimension (IM,JM,LM) :: ph, pih, mass, imass, t, dzet, qddf3, rh, dp, dm
78 real(8), dimension (IM,JM,LM+1) :: zet
79 real(8), dimension (IM,JM,0:LM) :: p, pi
80 real(8), dimension (IM,JM,LM) :: qs, dqsdt, dqs
81 
82 real(8), dimension(IM,JM) :: vmip
83 
84 real(8) :: cf_tot !Precip amounts and fall rate
85 
86 real(8) :: alpha, alhx3, rhcrit
87 real(8) :: aa, bb !Microphyiscal constants
88 
89 real(8), parameter :: pi_0 = 4.*atan(1.)
90 real(8), parameter :: rho_w = 1.0e3 !Density of liquid water in kg/m^3
91 real(8) :: t_ice_max
92 
93 !PHYSPARAMS constants
94 real(8) :: cnv_beta,anv_beta,ls_beta,rh00,c_00,lwcrit,c_acc,c_ev_r,c_ev_s,cldvol2frc
95 real(8) :: rhsup_ice,shr_evap_fac,min_cld_water,cld_evp_eff,ls_sdqv2,ls_sdqv3,ls_sdqvt1
96 real(8) :: anv_sdqv2,anv_sdqv3,anv_sdqvt1,anv_to_ls,n_warm,n_ice,n_anvil,n_pbl
97 real(8) :: anv_icefall_c,ls_icefall_c,revap_off_p,cnvenvfc,wrhodep,t_ice_all,cnviceparam
98 real(8) :: cnvddrfc,anvddrfc,lsddrfc,minrhcrit,maxrhcrit,turnrhcrit,maxrhcritland
99 real(8) :: min_rl,min_ri,max_rl,max_ri,ri_anv
100 integer :: nsmax,disable_rad,icefrpwr,tanhrhcrit,fr_ls_wat,fr_ls_ice,fr_an_wat,fr_an_ice,pdfflag
101 
102 real(8) :: lsenvfc, anvenvfc
103 real(8) :: qrn_cu, qsn_cu, qrn_an, qsn_an, qrn_ls, qsn_ls, qrn_cu_1d
104 real(8) :: qt_tmpi_1, qt_tmpi_2, qlt_tmp, qit_tmp
105 
106 real(8) :: prn_above_cu_new, prn_above_an_new, prn_above_ls_new
107 real(8) :: prn_above_cu_old, prn_above_an_old, prn_above_ls_old
108 real(8) :: psn_above_cu_new, psn_above_an_new, psn_above_ls_new
109 real(8) :: psn_above_cu_old, psn_above_an_old, psn_above_ls_old
110 real(8) :: evap_dd_cu_above_new, evap_dd_an_above_new, evap_dd_ls_above_new
111 real(8) :: evap_dd_cu_above_old, evap_dd_an_above_old, evap_dd_ls_above_old
112 real(8) :: subl_dd_cu_above_new, subl_dd_an_above_new, subl_dd_ls_above_new
113 real(8) :: subl_dd_cu_above_old, subl_dd_an_above_old, subl_dd_ls_above_old
114 
115 real(8) :: area_ls_prc1, area_upd_prc1, area_anv_prc1
116 real(8) :: tot_prec_upd, tot_prec_anv, tot_prec_ls, area_upd_prc, area_anv_prc, area_ls_prc
117 real(8) :: qtmp2
118 
119 real(8) :: rhexcess, tpw, negtpw
120 
121  !LS_CLOUD FILTERING
122  integer :: ii, cloud_pertmod
123  real(8) :: xx(8)
124  real(8) :: ttraj, qtraj, qi_lstraj, qi_contraj, ql_lstraj, ql_contraj, cf_lstraj, cf_contraj, phtraj
125  real(8) :: tpert, qpert, qi_lspert, qi_conpert, ql_lspert, ql_conpert, cf_lspert, cf_conpert
126  real(8) :: jacobian(8,8), a(8,8)
127 
128 
129 
130  !Total Filtering
131  real(8) :: totfilt_t, totfilt_ql, totfilt_qi
132  real(8) :: t_p_preall, ql_ls_p_preall, ql_con_p_preall, qi_ls_p_preall, qi_con_p_preall
133 
134  !Sink Filtering
135  real(8) :: sinkfilt_ql, sinkfilt_qi, sinkfilt_cf
136  real(8) :: t_p_presink, q_p_presink
137  real(8) :: ql_ls_p_presink, ql_con_p_presink
138  real(8) :: qi_ls_p_presink, qi_con_p_presink
139  real(8) :: cf_con_p_presink
140 
141 
142 !Highest level of calculations
143 ktop = 30
144 
145 !Get Constants from CLOUDPARAMS
146 cnv_beta = physparams(1) ! Area factor for convective rain showers (non-dim)
147 anv_beta = physparams(2) ! Area factor for anvil rain showers (non-dim)
148 ls_beta = physparams(3) ! Area factor for Large Scale rain showers (non-dim)
149 rh00 = physparams(4) ! Critical relative humidity
150 c_00 = physparams(5)
151 lwcrit = physparams(6)
152 c_acc = physparams(7)
153 c_ev_r = physparams(8)
154 c_ev_s = physparams(56)
155 cldvol2frc = physparams(9)
156 rhsup_ice = physparams(10)
157 shr_evap_fac = physparams(11)
158 min_cld_water = physparams(12)
159 cld_evp_eff = physparams(13)
160 nsmax = int( physparams(14) )
161 ls_sdqv2 = physparams(15)
162 ls_sdqv3 = physparams(16)
163 ls_sdqvt1 = physparams(17)
164 anv_sdqv2 = physparams(18)
165 anv_sdqv3 = physparams(19)
166 anv_sdqvt1 = physparams(20)
167 anv_to_ls = physparams(21)
168 n_warm = physparams(22)
169 n_ice = physparams(23)
170 n_anvil = physparams(24)
171 n_pbl = physparams(25)
172 disable_rad = int( physparams(26) )
173 anv_icefall_c = physparams(28)
174 ls_icefall_c = physparams(29)
175 revap_off_p = physparams(30)
176 cnvenvfc = physparams(31)
177 wrhodep = physparams(32)
178 t_ice_all = physparams(33) + cons_tice
179 cnviceparam = physparams(34)
180 icefrpwr = int( physparams(35) + .001 )
181 cnvddrfc = physparams(36)
182 anvddrfc = physparams(37)
183 lsddrfc = physparams(38)
184 tanhrhcrit = int( physparams(41) )
185 minrhcrit = physparams(42)
186 maxrhcrit = physparams(43)
187 turnrhcrit = physparams(45)
188 maxrhcritland = physparams(46)
189 fr_ls_wat = int( physparams(47) )
190 fr_ls_ice = int( physparams(48) )
191 fr_an_wat = int( physparams(49) )
192 fr_an_ice = int( physparams(50) )
193 min_rl = physparams(51)
194 min_ri = physparams(52)
195 max_rl = physparams(53)
196 max_ri = physparams(54)
197 ri_anv = physparams(55)
198 pdfflag = int(physparams(57))
199 
200 t_ice_max = cons_tice
201 
202 !Initialize the saving of downdraft values.
203 prn_above_cu_new = 0.
204 prn_above_an_new = 0.
205 prn_above_ls_new = 0.
206 psn_above_cu_new = 0.
207 psn_above_an_new = 0.
208 psn_above_ls_new = 0.
209 evap_dd_cu_above_new = 0.
210 evap_dd_an_above_new = 0.
211 evap_dd_ls_above_new = 0.
212 subl_dd_cu_above_new = 0.
213 subl_dd_an_above_new = 0.
214 subl_dd_ls_above_new = 0.
215 
216 !Convert to hPa and average pressure to temperature levels
217 p = ple*0.01
218 ph = 0.5*(p(:,:,0:lm-1) + p(:,:,1:lm ) )
219 
220 !Calculate Exner Pressure at temperature levels
221 pi = (p /1000.)**(cons_rgas/cons_cp)
222 pih = (ph/1000.)**(cons_rgas/cons_cp)
223 
224 !Calculate temperature
225 t = th*pih
226 
227 !Compute QS and DQSDT
228  call dqsat_bac(dqsdt,qs,t,ph,im,jm,lm,estblx,cons_h2omw,cons_airmw)
229 
230 !Relative humidity
231 rh = q/qs
232 
233 !Compute layer mass and 1/mass
234 mass = ( p(:,:,1:lm) - p(:,:,0:lm-1) )*100./cons_grav
235 imass = 1 / mass
236 
237 !Level thickness
238 dzet(:,:,1:lm) = th(:,:,1:lm) * (pi(:,:,1:lm) - pi(:,:,0:lm-1)) * cons_cp/cons_grav
239 
240 !Level heights
241 zet(:,:,lm+1) = 0.0
242 DO k = lm, 1, -1
243  zet(:,:,k) = zet(:,:,k+1)+dzet(:,:,k)
244 END DO
245 
246 WHERE ( zet(:,:,1:lm) < 3000. )
247  qddf3 = -( zet(:,:,1:lm)-3000. ) * zet(:,:,1:lm) * mass
248 ELSEWHERE
249  qddf3 = 0.
250 END WHERE
251 
252 DO i = 1,im
253  DO j = 1,jm
254  vmip(i,j) = sum( qddf3(i,j,:) )
255  END DO
256 END DO
257 DO k = 1,lm
258  qddf3(:,:,k) = qddf3(:,:,k) / vmip
259 END DO
260 
261 !Pressure and mass thickness for use in cleanup.
262 dp = ( ple(:,:,1:lm)-ple(:,:,0:lm-1) )
263 dm = dp*(1./cons_grav)
264 
265 !Begin loop over all grid boxes.
266 DO i = 1,im
267  DO j = 1,jm
268  DO k = ktop,lm
269 
270  !Save the inputs to the scheme for filtering
271  t_p_preall = t(i, j, k)
272  ql_ls_p_preall = ql_ls(i, j, k)
273  ql_con_p_preall = ql_con(i, j, k)
274  qi_ls_p_preall = qi_ls(i, j, k)
275  qi_con_p_preall = qi_con(i, j, k)
276 
277  if (k == ktop) then
278  tot_prec_upd = 0.
279  tot_prec_anv = 0.
280  tot_prec_ls = 0.
281 
282  area_upd_prc = 0.
283  area_anv_prc = 0.
284  area_ls_prc = 0.
285  end if
286 
287  !Initialize precips, except QRN_CU which comes from RAS
288  qrn_ls = 0.
289  qrn_an = 0.
290  qrn_cu_1d = 0.
291  qsn_ls = 0.
292  qsn_an = 0.
293  qsn_cu = 0.
294 
295  qrn_cu_1d = cnv_prc3(i,j,k) !Ras Rain
296 
297  !Tidy up where fractions or cloud is too low
298  call cloud_tidy( q(i,j,k), &
299  t(i,j,k), &
300  ql_ls(i,j,k), &
301  qi_ls(i,j,k), &
302  cf_ls(i,j,k), &
303  ql_con(i,j,k), &
304  qi_con(i,j,k), &
305  cf_con(i,j,k), &
306  cons_alhl, &
307  cons_alhs, &
308  cons_cp )
309 
310  !Phase changes for large scale cloud.
311  call meltfreeze ( dt, &
312  t(i,j,k), &
313  ql_ls(i,j,k), &
314  qi_ls(i,j,k), &
315  t_ice_all, &
316  t_ice_max, &
317  icefrpwr, &
318  cons_alhl, &
319  cons_alhs, &
320  cons_cp )
321 
322  !Phase changes for convective cloud.
323  call meltfreeze ( dt, &
324  t(i,j,k), &
325  ql_con(i,j,k), &
326  qi_con(i,j,k), &
327  t_ice_all, &
328  t_ice_max, &
329  icefrpwr, &
330  cons_alhl, &
331  cons_alhs, &
332  cons_cp )
333 
334 
335 
336 
337  !STAGE 1 - Compute convective clouds from RAS diagnostics
338  call convec_src( dt, &
339  mass(i,j,k), &
340  imass(i,j,k), &
341  t(i,j,k), &
342  q(i,j,k), &
343  cnv_dqldt(i,j,k), &
344  cnv_mfd(i,j,k), &
345  ql_con(i,j,k), &
346  qi_con(i,j,k), &
347  cf_con(i,j,k), &
348  qs(i,j,k), &
349  cons_alhs, &
350  cons_alhl, &
351  cons_cp, &
352  t_ice_all, &
353  t_ice_max, &
354  icefrpwr )
355 
356 
357  !STAGE 2a - Get PDF attributes
358  call pdf_width ( ph(i,j,k), &
359  frland(i,j), &
360  maxrhcrit, &
361  maxrhcritland, &
362  turnrhcrit, &
363  minrhcrit, &
364  pi_0, &
365  alpha )
366 
367  alpha = max( alpha , 1.0 - rh00 )
368  rhcrit = 1.0 - alpha
369 
370 
371  cloud_pertmod = 1
372 
373  !STAGE 2b - Use PDF to compute large scale cloud effects and diagnostics,
374  ! also update convection clouds
375  call ls_cloud( dt, &
376  alpha, &
377  pdfflag, &
378  ph(i,j,k), &
379  t(i,j,k), &
380  q(i,j,k), &
381  ql_ls(i,j,k), &
382  ql_con(i,j,k), &
383  qi_ls(i,j,k), &
384  qi_con(i,j,k), &
385  cf_ls(i,j,k), &
386  cf_con(i,j,k), &
387  cons_alhl, &
388  cons_alhf, &
389  cons_alhs, &
390  cons_cp, &
391  cons_h2omw, &
392  cons_airmw, &
393  t_ice_all, &
394  t_ice_max, &
395  icefrpwr, &
396  estblx, &
397  cloud_pertmod, &
398  do_moist_physics )
399 
400  !SAVE PRESINKS INPUTS
401  t_p_presink = t(i,j,k)
402  q_p_presink = q(i,j,k)
403  qi_ls_p_presink = qi_ls(i,j,k)
404  qi_con_p_presink = qi_con(i,j,k)
405  ql_ls_p_presink = ql_ls(i,j,k)
406  ql_con_p_presink = ql_con(i,j,k)
407  cf_con_p_presink = cf_con(i,j,k)
408 
409 
410  !Clean up where too much overall cloud.
411  cf_tot = cf_ls(i,j,k) + cf_con(i,j,k)
412  IF ( cf_tot > 1.00 ) THEN
413  cf_ls(i,j,k) = cf_ls(i,j,k)*(1.00 / cf_tot )
414  cf_con(i,j,k) = cf_con(i,j,k)*(1.00 / cf_tot )
415  END IF
416  cf_tot = cf_ls(i,j,k) + cf_con(i,j,k)
417 
418 
419  !STAGE 3 - Evap, Sublimation and Autoconversion
420 
421  !Evaporation and sublimation of anvil cloud
422  call evap_cnv( dt, &
423  rhcrit, &
424  ph(i,j,k), &
425  t(i,j,k), &
426  q(i,j,k), &
427  ql_con(i,j,k), &
428  qi_con(i,j,k), &
429  cf_con(i,j,k), &
430  cf_ls(i,j,k), &
431  qs(i,j,k), &
432  rho_w, &
433  cld_evp_eff, &
434  cons_h2omw, &
435  cons_airmw, &
436  cons_alhl, &
437  cons_rvap, &
438  cons_rgas, &
439  cons_pi, &
440  cons_cp )
441 
442  call subl_cnv( dt, &
443  rhcrit, &
444  ph(i,j,k), &
445  t(i,j,k), &
446  q(i,j,k), &
447  ql_con(i,j,k), &
448  qi_con(i,j,k), &
449  cf_con(i,j,k), &
450  cf_ls(i,j,k), &
451  qs(i,j,k), &
452  rho_w, &
453  cld_evp_eff, &
454  cons_h2omw, &
455  cons_airmw, &
456  cons_alhl, &
457  cons_rvap, &
458  cons_rgas, &
459  cons_pi, &
460  cons_cp, &
461  cons_alhs )
462 
463  !Autoconversion
464  call autoconversion_ls ( dt, &
465  ql_ls(i,j,k), &
466  qrn_ls, &
467  t(i,j,k), &
468  ph(i,j,k), &
469  cf_ls(i,j,k), &
470  ls_sdqv2, &
471  ls_sdqv3, &
472  ls_sdqvt1, &
473  c_00, &
474  lwcrit, &
475  dzet(i,j,k) )
476 
477  call autoconversion_cnv ( dt, &
478  ql_con(i,j,k), &
479  qrn_an, &
480  t(i,j,k), &
481  ph(i,j,k), &
482  cf_con(i,j,k), &
483  anv_sdqv2, &
484  anv_sdqv3, &
485  anv_sdqvt1, &
486  c_00, &
487  lwcrit, &
488  dzet(i,j,k) )
489 
490 
491  !STAGE 4 - Fall and Re-evaporation of precip
492  call ice_settlefall_cnv( wrhodep, &
493  qi_con(i,j,k), &
494  ph(i,j,k), &
495  t(i,j,k), &
496  cf_con(i,j,k), &
497  cons_rgas, &
498  khu(i,j), &
499  khl(i,j), &
500  k, &
501  dt, &
502  dzet(i,j,k), &
503  qsn_an, &
504  anv_icefall_c )
505 
506  call ice_settlefall_ls( wrhodep, &
507  qi_ls(i,j,k), &
508  ph(i,j,k), &
509  t(i,j,k), &
510  cf_ls(i,j,k), &
511  cons_rgas, &
512  khu(i,j), &
513  khl(i,j), &
514  k, &
515  dt, &
516  dzet(i,j,k), &
517  qsn_ls, &
518  ls_icefall_c )
519 
520  !"Freeze" out any conv. precip, not done in RAS. This is
521  ! precip w/ large particles, so freezing is strict.
522  qtmp2 = 0.
523  if ( t(i,j,k) < cons_tice ) then
524  qtmp2 = qrn_cu_1d
525  qsn_cu = qrn_cu_1d
526  qrn_cu_1d = 0.
527  t(i,j,k) = t(i,j,k) + qsn_cu*(cons_alhs-cons_alhl)/cons_cp
528  end if
529 
530  !Area
531  area_ls_prc1 = 0.0
532  area_upd_prc1 = 0.0
533  area_anv_prc1 = 0.0
534 
535  tot_prec_upd = tot_prec_upd + ( ( qrn_cu_1d + qsn_cu ) * mass(i,j,k) )
536  area_upd_prc = area_upd_prc + ( cnv_updf(i,j,k)* ( qrn_cu_1d + qsn_cu )* mass(i,j,k) )
537 
538  tot_prec_anv = tot_prec_anv + ( ( qrn_an + qsn_an ) * mass(i,j,k) )
539  area_anv_prc = area_anv_prc + ( cf_con(i,j,k)* ( qrn_an + qsn_an )* mass(i,j,k) )
540 
541  tot_prec_ls = tot_prec_ls + ( ( qrn_ls + qsn_ls ) * mass(i,j,k) )
542  area_ls_prc = area_ls_prc + ( cf_ls(i,j,k)* ( qrn_ls + qsn_ls )* mass(i,j,k) )
543 
544  if ( tot_prec_anv > 0.0 ) area_anv_prc1 = max( area_anv_prc/tot_prec_anv, 1.e-6 )
545  if ( tot_prec_upd > 0.0 ) area_upd_prc1 = max( area_upd_prc/tot_prec_upd, 1.e-6 )
546  if ( tot_prec_ls > 0.0 ) area_ls_prc1 = max( area_ls_prc/tot_prec_ls, 1.e-6 )
547 
548  area_ls_prc1 = ls_beta * area_ls_prc1
549  area_upd_prc1 = cnv_beta * area_upd_prc1
550  area_anv_prc1 = anv_beta * area_anv_prc1
551 
552  IF (k == lm) THEN ! We have accumulated over the whole column
553  if ( tot_prec_anv > 0.0 ) area_anv_prc = max( area_anv_prc/tot_prec_anv, 1.e-6 )
554  if ( tot_prec_upd > 0.0 ) area_upd_prc = max( area_upd_prc/tot_prec_upd, 1.e-6 )
555  if ( tot_prec_ls > 0.0 ) area_ls_prc = max( area_ls_prc/tot_prec_ls, 1.e-6 )
556 
557  area_ls_prc = ls_beta * area_ls_prc
558  area_upd_prc = cnv_beta * area_upd_prc
559  area_anv_prc = anv_beta * area_anv_prc
560  END IF
561 
562  !Get micro-physical constants
563  call cons_alhx( t(i,j,k), &
564  alhx3, &
565  t_ice_max, &
566  t_ice_all, &
567  cons_alhs, &
568  cons_alhl )
569 
570  call cons_microphys( t(i,j,k), &
571  ph(i,j,k), &
572  qs(i,j,k), &
573  aa, &
574  bb, &
575  cons_h2omw, &
576  cons_airmw, &
577  cons_rvap, &
578  alhx3 )
579 
580 
581  !Precip Scheme Expects Total Cloud Liquid
582  qlt_tmp = ql_ls(i,j,k) + ql_con(i,j,k)
583  qit_tmp = qi_ls(i,j,k) + qi_con(i,j,k)
584 
585  prn_above_cu_old = prn_above_cu_new
586  psn_above_cu_old = psn_above_cu_new
587  evap_dd_cu_above_old = evap_dd_cu_above_new
588  subl_dd_cu_above_old = subl_dd_cu_above_new
589 
590  !Precip and Evap for Convection
591  call precipandevap( k , &
592  ktop , &
593  lm , &
594  dt , &
595  frland(i,j) , &
596  rhcrit , &
597  qrn_cu_1d , &
598  qsn_cu , &
599  qlt_tmp , &
600  qit_tmp , &
601  t(i,j,k) , &
602  q(i,j,k) , &
603  mass(i,j,k) , &
604  imass(i,j,k) , &
605  ph(i,j,k) , &
606  dzet(i,j,k) , &
607  qddf3(i,j,k) , &
608  aa , &
609  bb , &
610  area_upd_prc1 , &
611  prn_above_cu_old , &
612  prn_above_cu_new , &
613  psn_above_cu_old , &
614  psn_above_cu_new , &
615  evap_dd_cu_above_old , &
616  evap_dd_cu_above_new , &
617  subl_dd_cu_above_old , &
618  subl_dd_cu_above_new , &
619  cnvenvfc , &
620  cnvddrfc , &
621  cons_alhf , &
622  cons_alhs , &
623  cons_alhl , &
624  cons_cp , &
625  cons_tice , &
626  cons_h2omw , &
627  cons_airmw , &
628  revap_off_p , &
629  c_acc , &
630  c_ev_r , &
631  c_ev_s , &
632  rho_w , &
633  estblx )
634 
635  prn_above_an_old = prn_above_an_new
636  psn_above_an_old = psn_above_an_new
637  evap_dd_an_above_old = evap_dd_an_above_new
638  subl_dd_an_above_old = subl_dd_an_above_new
639 
640  !Precip and Evap for Anvil
641  anvenvfc = 1.0
642  call precipandevap( k , &
643  ktop , &
644  lm , &
645  dt , &
646  frland(i,j) , &
647  rhcrit , &
648  qrn_an , &
649  qsn_an , &
650  qlt_tmp , &
651  qit_tmp , &
652  t(i,j,k) , &
653  q(i,j,k) , &
654  mass(i,j,k) , &
655  imass(i,j,k) , &
656  ph(i,j,k) , &
657  dzet(i,j,k) , &
658  qddf3(i,j,k) , &
659  aa , &
660  bb , &
661  area_anv_prc1 , &
662  prn_above_an_old , &
663  prn_above_an_new , &
664  psn_above_an_old , &
665  psn_above_an_new , &
666  evap_dd_an_above_old , &
667  evap_dd_an_above_new , &
668  subl_dd_an_above_old , &
669  subl_dd_an_above_new , &
670  anvenvfc , &
671  anvddrfc , &
672  cons_alhf , &
673  cons_alhs , &
674  cons_alhl , &
675  cons_cp , &
676  cons_tice , &
677  cons_h2omw , &
678  cons_airmw , &
679  revap_off_p , &
680  c_acc , &
681  c_ev_r , &
682  c_ev_s , &
683  rho_w , &
684  estblx )
685 
686  prn_above_ls_old = prn_above_ls_new
687  psn_above_ls_old = psn_above_ls_new
688  evap_dd_ls_above_old = evap_dd_ls_above_new
689  subl_dd_ls_above_old = subl_dd_ls_above_new
690 
691  !Precip and Evap for Large Scale
692  lsenvfc = 1.0
693  call precipandevap( k , &
694  ktop , &
695  lm , &
696  dt , &
697  frland(i,j) , &
698  rhcrit , &
699  qrn_ls , &
700  qsn_ls , &
701  qlt_tmp , &
702  qit_tmp , &
703  t(i,j,k) , &
704  q(i,j,k) , &
705  mass(i,j,k) , &
706  imass(i,j,k) , &
707  ph(i,j,k) , &
708  dzet(i,j,k) , &
709  qddf3(i,j,k) , &
710  aa , &
711  bb , &
712  area_ls_prc1 , &
713  prn_above_ls_old , &
714  prn_above_ls_new , &
715  psn_above_ls_old , &
716  psn_above_ls_new , &
717  evap_dd_ls_above_old , &
718  evap_dd_ls_above_new , &
719  subl_dd_ls_above_old , &
720  subl_dd_ls_above_new , &
721  lsenvfc , &
722  lsddrfc , &
723  cons_alhf , &
724  cons_alhs , &
725  cons_alhl , &
726  cons_cp , &
727  cons_tice , &
728  cons_h2omw , &
729  cons_airmw , &
730  revap_off_p , &
731  c_acc , &
732  c_ev_r , &
733  c_ev_s , &
734  rho_w , &
735  estblx )
736 
737 
738  if ( (ql_ls(i,j,k) + ql_con(i,j,k)) > 0.00 ) then
739  qt_tmpi_1 = 1./(ql_ls(i,j,k) + ql_con(i,j,k))
740  else
741  qt_tmpi_1 = 0.0
742  endif
743  ql_ls(i,j,k) = ql_ls(i,j,k) * qlt_tmp * qt_tmpi_1
744  ql_con(i,j,k) = ql_con(i,j,k) * qlt_tmp * qt_tmpi_1
745 
746  if ( (qi_ls(i,j,k) + qi_con(i,j,k)) > 0.00 ) then
747  qt_tmpi_2 = 1./(qi_ls(i,j,k) + qi_con(i,j,k))
748  else
749  qt_tmpi_2 = 0.0
750  endif
751  qi_ls(i,j,k) = qi_ls(i,j,k) * qit_tmp * qt_tmpi_2
752  qi_con(i,j,k) = qi_con(i,j,k) * qit_tmp * qt_tmpi_2
753 
754 
755 
756  !SINK FILTERING
757 ! if (do_moist_physics == 1) then
758 ! SINKfilt_ql = 0.65
759 ! SINKfilt_qi = 0.65
760 ! SINKfilt_CF = 1.0
761 ! elseif (do_moist_physics == 2) then
762 ! SINKfilt_ql = 0.9
763 ! SINKfilt_qi = 0.9
764 ! SINKfilt_CF = 1.0
765 ! endif
766 !
767 ! if (k < 50) then
768 ! qi_ls(i,j,k) = SINKfilt_qi * qi_ls(i,j,k) + (1.0-SINKfilt_qi) * qi_ls_p_presink
769 ! qi_con(i,j,k) = SINKfilt_qi * qi_con(i,j,k) + (1.0-SINKfilt_qi) * qi_con_p_presink
770 ! q(i,j,k) = SINKfilt_qi * q(i,j,k) + (1.0-SINKfilt_qi) * q_p_presink
771 ! endif
772 !
773 ! if ( abs(k-62) .le. 2) then
774 ! ql_ls(i,j,k) = SINKfilt_ql * ql_ls(i,j,k) + (1.0-SINKfilt_ql) * ql_ls_p_presink
775 ! ql_con(i,j,k) = SINKfilt_ql * ql_con(i,j,k) + (1.0-SINKfilt_ql) * ql_con_p_presink
776 ! endif
777 !
778 ! cf_con(i,j,k) = SINKfilt_CF * cf_con(i,j,k) + (1.0-SINKfilt_CF) * cf_con_p_presink
779 !
780 !
781 !
782 ! !TOTAL FILTERING
783 ! TOTfilt_T = 0.25
784 ! t(i,j,k) = TOTfilt_T * t(i,j,k) + (1.0-TOTfilt_T) * t_p_preall
785 !
786 ! if (do_moist_physics == 1) then
787 ! TOTfilt_ql = 0.75
788 ! TOTfilt_qi = 1.0
789 ! elseif (do_moist_physics == 2) then
790 ! TOTfilt_ql = 0.5
791 ! TOTfilt_qi = 0.5
792 ! endif
793 !
794 ! ql_ls(i,j,k) = TOTfilt_ql * ql_ls(i,j,k) + (1.0-TOTfilt_ql) * ql_ls_p_preall
795 ! ql_con(i,j,k) = TOTfilt_ql * ql_con(i,j,k) + (1.0-TOTfilt_ql) * ql_con_p_preall
796 !
797 ! qi_ls(i,j,k) = TOTfilt_qi * qi_ls(i,j,k) + (1.0-TOTfilt_qi) * qi_ls_p_preall
798 ! qi_con(i,j,k) = TOTfilt_qi * qi_con(i,j,k) + (1.0-TOTfilt_qi) * qi_con_p_preall
799 
800 
801 
802  endDO
803  endDO
804 endDO
805 
806 !Clean up of excess relative humidity
807 rhexcess = 1.1
808 call dqsat_bac(dqsdt,qs,t,ph,im,jm,lm,estblx,cons_h2omw,cons_airmw)
809 
810 where ( q > rhexcess*qs )
811  dqs = (q - rhexcess*qs)/( 1.0 + rhexcess*dqsdt*cons_alhl/cons_cp )
812 elsewhere
813  dqs = 0.0
814 endwhere
815 
816 q = q - dqs
817 t = t + (cons_alhl/cons_cp)*dqs
818 
819 !Clean up Q<0
820 do j=1,jm
821  do i=1,im
822 
823  !Total precipitable water
824  tpw = sum( q(i,j,:)*dm(i,j,:) )
825 
826  negtpw = 0.
827  do l=1,lm
828  if ( q(i,j,l) < 0.0 ) then
829  negtpw = negtpw + ( q(i,j,l)*dm( i,j,l ) )
830  q(i,j,l) = 0.0
831  endif
832  enddo
833 
834  do l=1,lm
835  if ( q(i,j,l) >= 0.0 ) then
836  q(i,j,l) = q(i,j,l)*( 1.0+negtpw/(tpw-negtpw) )
837  endif
838  enddo
839 
840  end do
841 end do
842 
843 !Convert temperature back to potential temperature
844 th = t/pih
845 
846 end subroutine cloud_driver
847 
848 
849 
850 !SUBROUTINES
851 subroutine cloud_tidy( QV, TE, QLC, QIC, CF, QLA, QIA, AF, CONS_ALHL, CONS_ALHS, CONS_CP )
852 
853 Implicit None
854 
855 real(8), intent(inout) :: te,qv,qlc,cf,qla,af,qic,qia
856 real(8), intent(in) :: cons_alhl, cons_alhs, cons_cp
857 
858  !Fix if Anvil cloud fraction too small
859  if (af < 1.e-5) then
860  qv = qv + qla + qia
861  te = te - (cons_alhl/cons_cp)*qla - (cons_alhs/cons_cp)*qia
862  af = 0.
863  qla = 0.
864  qia = 0.
865  end if
866 
867  !Fix if LS cloud fraction too small
868 ! if ( CF < 1.E-5 ) then
869 ! QV = QV + QLC + QIC
870 ! TE = TE - (CONS_ALHL/CONS_CP)*QLC - (CONS_ALHS/CONS_CP)*QIC
871 ! CF = 0.
872 ! QLC = 0.
873 ! QIC = 0.
874 ! end if
875 
876  !LS LIQUID too small
877  if ( qlc < 1.e-8 ) then
878  qv = qv + qlc
879  te = te - (cons_alhl/cons_cp)*qlc
880  qlc = 0.
881  end if
882  !LS ICE too small
883  if ( qic < 1.e-8 ) then
884  qv = qv + qic
885  te = te - (cons_alhs/cons_cp)*qic
886  qic = 0.
887  end if
888 
889  !Anvil LIQUID too small
890  if ( qla < 1.e-8 ) then
891  qv = qv + qla
892  te = te - (cons_alhl/cons_cp)*qla
893  qla = 0.
894  end if
895  !Anvil ICE too small
896  if ( qia < 1.e-8 ) then
897  qv = qv + qia
898  te = te - (cons_alhs/cons_cp)*qia
899  qia = 0.
900  end if
901 
902  !Fix ALL cloud quants if Anvil cloud LIQUID+ICE too small
903  if ( ( qla + qia ) < 1.e-8 ) then
904  qv = qv + qla + qia
905  te = te - (cons_alhl/cons_cp)*qla - (cons_alhs/cons_cp)*qia
906  af = 0.
907  qla = 0.
908  qia = 0.
909  end if
910  !Ditto if LS cloud LIQUID+ICE too small
911  if ( ( qlc + qic ) < 1.e-8 ) then
912  qv = qv + qlc + qic
913  te = te - (cons_alhl/cons_cp)*qlc - (cons_alhs/cons_cp)*qic
914  cf = 0.
915  qlc = 0.
916  qic = 0.
917  end if
918 
919 end subroutine cloud_tidy
920 
921 subroutine meltfreeze( DT, TE, QL, QI, T_ICE_ALL, T_ICE_MAX, ICEFRPWR, &
922  CONS_ALHL, CONS_ALHS, CONS_CP )
923 
924 Implicit None
925 
926 !Inputs
927 real(8), intent(in) :: dt, t_ice_all, t_ice_max
928 integer, intent(in) :: icefrpwr
929 real(8), intent(in) :: cons_alhl, cons_alhs, cons_cp
930 
931 !Prognostic
932 real(8), intent(inout) :: te,ql,qi
933 
934 !Locals
935 real(8) :: fqi, dqil
936 real(8), parameter :: taufrz = 1000.
937 
938  fqi = 0.0
939  dqil = 0.0
940 
941  call get_ice_fraction( te, t_ice_all, t_ice_max, icefrpwr, fqi )
942 
943  !Freeze liquid
944  if ( te <= t_ice_max ) then
945  dqil = ql *(1.0 - exp( -dt * fqi / taufrz ) )
946  end if
947 
948  dqil = max( 0., dqil )
949  qi = qi + dqil
950  ql = ql - dqil
951  te = te + (cons_alhs-cons_alhl)*dqil/cons_cp
952 
953  dqil = 0.
954  !Melt ice instantly above 0^C
955  if ( te > t_ice_max ) then
956  dqil = -qi
957  end if
958 
959  dqil = min( 0., dqil )
960  qi = qi + dqil
961  ql = ql - dqil
962  te = te + (cons_alhs-cons_alhl)*dqil/cons_cp
963 
964 end subroutine meltfreeze
965 
966 
967 subroutine convec_src( DT, MASS, iMASS, TE, QV, DCF, DMF, QLA, QIA, AF, QS, CONS_ALHS, &
968  CONS_ALHL, CONS_CP, T_ICE_ALL, T_ICE_MAX, ICEFRPWR )
969 
970 IMPLICIT NONE
971 
972 !Inputs
973 real(8), intent(IN) :: dt, t_ice_all, t_ice_max
974 integer, intent(in) :: icefrpwr
975 real(8), intent(in) :: mass, imass, qs
976 real(8), intent(in) :: dmf, dcf
977 real(8), intent(in) :: cons_alhs, cons_alhl, cons_cp
978 
979 !Prognostic
980 real(8), intent(inout) :: te, qv
981 real(8), intent(inout) :: qla, qia, af
982 
983 !Locals
984 real(8), parameter :: minrhx = 0.001 !Minimum allowed env RH
985 real(8) :: tend, qvx, fqi
986 
987  !Namelist
988  !DT - Timestep
989  !MASS - Level Mass
990  !iMASS - 1/Level Mass
991  !TE - Temperature
992  !QV - Specific Humidity
993  !DCF - CNV_DQL from RAS
994  !DMF - CNV_MFD from RAS
995  !QLA - Convective cloud liquid water
996  !QIA - Convective cloud liquid ice
997  !AF - Convective cloud fraction
998  !QS - Qsat
999 
1000  !Zero out locals
1001  tend = 0.0
1002  qvx = 0.0
1003  fqi = 0.0
1004 
1005  !Addition of condensate from RAS
1006  tend = dcf*imass
1007  call get_ice_fraction( te, t_ice_all, t_ice_max, icefrpwr, fqi )
1008  qla = qla + (1.0-fqi)* tend*dt
1009  qia = qia + fqi * tend*dt
1010 
1011  !Convective condensation has never frozen so latent heat of fusion
1012  te = te + (cons_alhs-cons_alhl) * fqi * tend*dt/cons_cp
1013 
1014  !Compute Tiedtke-style anvil fraction
1015  tend=dmf*imass
1016  af = af + tend*dt
1017  af = min( af , 0.99 )
1018 
1019  ! Check for funny (tiny, negative) external QV, resulting from assumed QV=QSAT within anvil.
1020  ! Simply constrain AF assume condensate just gets "packed" in
1021 
1022  if ( af < 1.0 ) then
1023  qvx = ( qv - qs * af )/(1.-af)
1024  else
1025  qvx = qs
1026  end if
1027 
1028  !If saturated over critial value and there is Anvil
1029  if ( (( qvx - minrhx*qs ) < 0.0 ) .and. (af > 0.) ) then
1030  af = (qv - minrhx*qs )/( qs*(1.0-minrhx) )
1031  end if
1032 
1033  if ( af < 0. ) then ! If still cant make suitable env RH then destroy anvil
1034  af = 0.0
1035  qv = qv + qla + qia
1036  te = te - (cons_alhl*qla + cons_alhs*qia)/cons_cp
1037  qla = 0.0
1038  qia = 0.0
1039  end if
1040 
1041 end subroutine convec_src
1042 
1043 
1044 
1045 subroutine pdf_width (PP,FRLAND,maxrhcrit,maxrhcritland,turnrhcrit,minrhcrit,pi_0,ALPHA )
1046 
1047  Implicit none
1048 
1049  !Inputs
1050  real(8), intent(in) :: pp, frland
1051 
1052  real(8), intent(in) :: maxrhcrit, maxrhcritland, turnrhcrit, minrhcrit, pi_0
1053 
1054  !Prognostic
1055  real(8), intent(inout) :: alpha
1056 
1057  !Locals
1058  real(8) :: a1, tempmaxrh
1059 
1060  !Namelist
1061  !PP - Pressure (hPa)
1062  !FRLAND - Fraction of land
1063  !maxrhcrit - Constant
1064  !maxrhcritland - Constant
1065  !turnrhcrit - Constant
1066  !minrhcrit - Constant
1067  !pi_0 - Constant
1068  !ALPHA - 1/2*PDFwidth so RH_crit=1.0-alpha
1069 
1070  ! Use Slingo-Ritter (1985) formulation for critical relative humidity
1071  ! array a1 holds the critical rh, ranges from 0.8 to 1
1072 
1073  tempmaxrh = maxrhcrit
1074  if (frland > 0.05) then
1075  tempmaxrh = maxrhcritland
1076  endif
1077 
1078  a1 = 1.0
1079  if (pp .le. turnrhcrit) then
1080 
1081  a1 = minrhcrit
1082 
1083  else
1084 
1085  a1 = minrhcrit + (tempmaxrh-minrhcrit)/(19.) * &
1086  ((atan( (2.*(pp- turnrhcrit)/(1020.-turnrhcrit)-1.) * &
1087  tan(20.*pi_0/21.-0.5*pi_0) ) + 0.5*pi_0) * 21./pi_0 - 1.)
1088 
1089  end if
1090 
1091  a1 = min(a1,1.)
1092  alpha = 1. - a1
1093 
1094  alpha = min( alpha , 0.25 ) ! restrict RHcrit to > 75%
1095 
1096 end subroutine pdf_width
1097 
1098 
1099 
1100 
1101 
1102 
1103 
1104 subroutine ls_cloud( DT, ALPHA, PDFSHAPE, PL, TE, QV, QCl, QAl, QCi, QAi, CF, AF, &
1105  CONS_ALHL, CONS_ALHF, CONS_ALHS, CONS_CP, CONS_H2OMW, CONS_AIRMW, T_ICE_ALL, &
1106  T_ICE_MAX, ICEFRPWR, ESTBLX, cloud_pertmod, dmp )
1107 
1108 IMPLICIT NONE
1109 
1110 !Inputs
1111 real(8), intent(in) :: dt, alpha, pl, t_ice_all, t_ice_max
1112 integer, intent(in) :: pdfshape, cloud_pertmod, dmp
1113 integer, intent(in) :: icefrpwr
1114 real(8), intent(in) :: cons_alhl, cons_alhf, cons_alhs, cons_cp, cons_h2omw, cons_airmw
1115 real(8), intent(in) :: estblx(:)
1116 
1117 !Prognostic
1118 real(8), intent(inout) :: te, qv, qcl, qci, qal, qai, cf, af
1119 
1120 !Locals
1121 integer :: n
1122 
1123 real(8) :: qco, cfo, qao, qt, qmx, qmn, dq
1124 real(8) :: teo, qsx, dqsx, qs, dqs, tmparr
1125 real(8) :: qcx, qvx, cfx, qax, qc, qa, fqi, fqi_a, dqai, dqal, dqci, dqcl
1126 
1127 real(8) :: ten, qsp, cfp, qvp, qcp
1128 real(8) :: tep, qsn, cfn, qvn, qcn
1129 real(8) :: alhx, sigmaqt1, sigmaqt2
1130 
1131 real(8), dimension(1) :: dqsx1,qsx1,teo1,pl1
1132 
1133  !Namelist
1134  !DT - Timestep
1135  !ALPHA - PDF half width
1136  !PL - Pressure (hPa)
1137  !TE - Temperature
1138  !QV - Specific humidity
1139  !QCl - Convective cloud liquid water
1140  !QAl - Large scale cloud liquid water
1141  !QCi - Convective cloud liquid ice
1142  !QAi - Large scale cloud liquid ice
1143  !CF - Large scale cloud fraction
1144  !AF - Convective cloud fraction
1145 
1146  qc = qcl + qci
1147  qa = qal + qai
1148 
1149  if ( qa > 0.0 ) then
1150  fqi_a = qai / qa
1151  else
1152  fqi_a = 0.0
1153  end if
1154 
1155  teo = te
1156 
1157  call dqsats_bac(dqsx, qsx,teo,pl,estblx,cons_h2omw, cons_airmw)
1158 
1159  if ( af < 1.0 ) then
1160  if (dmp == 1) then
1161  if ( (1.-af) .gt. 0.02) then
1162  tmparr = 1./(1.-af)
1163  else
1164  tmparr = 1./(1.-af)
1165  endif
1166  elseif (dmp == 2) then
1167  tmparr = 1./(1.-af)
1168  endif
1169  else
1170  tmparr = 0.0
1171  end if
1172 
1173  cfx = cf*tmparr
1174  qcx = qc*tmparr
1175  qvx = ( qv - qsx*af )*tmparr
1176 
1177  if ( af >= 1.0 ) then
1178  qvx = qsx*1.e-4
1179  end if
1180 
1181  if ( af > 0. ) then
1182  qax = qa/af
1183  else
1184  qax = 0.
1185  end if
1186 
1187  qt = qcx + qvx
1188 
1189  tep = teo
1190  qsn = qsx
1191  ten = teo
1192  cfn = cfx
1193  qvn = qvx
1194  qcn = qcx
1195 
1196  dqs = dqsx
1197 
1198  !Begin iteration
1199  !do n=1,4
1200  n = 1
1201 
1202  qsp = qsn
1203  qvp = qvn
1204  qcp = qcn
1205  cfp = cfn
1206 
1207  !Dont call again as not looping
1208  dqs = dqsx
1209  qsn = qsx
1210  !call DQSATs_BAC(DQS, QSn, TEn, PL, ESTBLX, CONS_H2OMW, CONS_AIRMW)
1211 
1212  tep = ten
1213  call get_ice_fraction( tep, t_ice_all, t_ice_max, icefrpwr, fqi )
1214 
1215  sigmaqt1 = alpha*qsn
1216  sigmaqt2 = alpha*qsn
1217 
1218  if (pdfshape .eq. 2) then
1219  !For triangular, symmetric: sigmaqt1 = sigmaqt2 = alpha*qsn (alpha is half width)
1220  !For triangular, skewed r : sigmaqt1 < sigmaqt2
1221  sigmaqt1 = alpha*qsn
1222  sigmaqt2 = alpha*qsn
1223  endif
1224 
1225  !Compute cloud fraction
1226  if (cloud_pertmod == 0) then
1227  call pdffrac(1,qt,sigmaqt1,sigmaqt2,qsn,cfn)
1228  elseif (cloud_pertmod == 1) then
1229  call pdffrac(4,qt,sigmaqt1,sigmaqt2,qsn,cfn)
1230  endif
1231 
1232  !Compute cloud condensate
1233  call pdfcondensate(pdfshape,qt,sigmaqt1,sigmaqt2,qsn,qcn)
1234 
1235  !Adjustments to anvil condensate due to the assumption of a stationary TOTAL
1236  !water PDF subject to a varying QSAT value during the iteration
1237  if ( af > 0. ) then
1238  qao = qax ! + QSx - QS
1239  else
1240  qao = 0.
1241  end if
1242 
1243  alhx = (1.0-fqi)*cons_alhl + fqi*cons_alhs
1244 
1245  if (pdfshape .eq. 1) then
1246  qcn = qcp + ( qcn - qcp ) / ( 1. - (cfn * (alpha-1.) - (qcn/qsn))*dqs*alhx/cons_cp)
1247  elseif (pdfshape .eq. 2) then
1248  !This next line needs correcting - need proper d(del qc)/dT derivative for triangular
1249  !for now, just use relaxation of 1/2.
1250  if (n.ne.4) then
1251  qcn = qcp + ( qcn - qcp ) *0.5
1252  endif
1253  endif
1254 
1255  qvn = qvp - (qcn - qcp)
1256  ten = tep + (1.0-fqi)*(cons_alhl/cons_cp)*( (qcn - qcp)*(1.-af) + (qao-qax)*af ) &
1257  + fqi* (cons_alhs/cons_cp)*( (qcn - qcp)*(1.-af) + (qao-qax)*af )
1258 
1259  !enddo ! qsat iteration
1260 
1261  cfo = cfn
1262  cf = cfn
1263  qco = qcn
1264  teo = ten
1265 
1266  ! Update prognostic variables. QCo, QAo become updated grid means.
1267  if ( af < 1.0 ) then
1268 
1269  cf = cfo * ( 1.-af)
1270  qco = qco * ( 1.-af)
1271  qao = qao * af
1272 
1273  else !Grid box filled with anvil
1274 
1275  ! Special case AF=1, i.e., box filled with anvil. Note: no guarantee QV_box > QS_box
1276 
1277  cf = 0. ! Remove any other cloud
1278  qao = qa + qc ! Add any LS condensate to anvil type
1279  qco = 0. ! Remove same from LS
1280 
1281  qt = qao + qv ! Total water
1282 
1283  ! Now set anvil condensate to any excess of total water over QSx (saturation value at top)
1284  qao = max( qt - qsx, 0. )
1285 
1286  end if
1287 
1288  !Partition new condensate into ice and liquid taking care to keep both >=0 separately.
1289  !New condensate can be less than old, so Delta can be < 0.
1290 
1291  qcx = qco - qc
1292  dqcl = (1.0-fqi)*qcx
1293  dqci = fqi *qcx
1294 
1295  !Large Scale Partition
1296  if ((qcl+dqcl)<0.) then
1297  dqci = dqci + (qcl+dqcl)
1298  dqcl = -qcl !== dQCl - (QCl+dQCl)
1299  end if
1300  if ((qci+dqci)<0.) then
1301  dqcl = dqcl + (qci+dqci)
1302  dqci = -qci !== dQCi - (QCi+dQCi)
1303  end if
1304 
1305  qax = qao - qa
1306  dqal = qax ! (1.0-fQi)*QAx
1307  dqai = 0. ! fQi * QAx
1308 
1309  !Convective partition
1310  if ((qal+dqal)<0.) then
1311  dqai = dqai + (qal+dqal)
1312  dqal = -qal
1313  end if
1314  if ((qai+dqai)<0.) then
1315  dqal = dqal + (qai+dqai)
1316  dqai = -qai
1317  end if
1318 
1319  ! Clean-up cloud if fractions are too small
1320  if ( af < 1.e-5 ) then
1321  dqai = -qai
1322  dqal = -qal
1323  end if
1324  if ( cf < 1.e-5 ) then
1325  dqci = -qci
1326  dqcl = -qcl
1327  end if
1328 
1329  qai = qai + dqai
1330  qal = qal + dqal
1331  qci = qci + dqci
1332  qcl = qcl + dqcl
1333 
1334  !Update specific humidity
1335  qv = qv - ( dqai+dqci+dqal+dqcl)
1336 
1337  !Update temperature
1338  te = te + (cons_alhl*( dqai+dqci+dqal+dqcl)+cons_alhf*(dqai+dqci))/ cons_cp
1339 
1340  !Take care of situations where QS moves past QA during QSAT iteration (QAo negative).
1341  !"Evaporate" offending QA
1342  if ( qao <= 0. ) then
1343  qv = qv + qai + qal
1344  te = te - (cons_alhs/cons_cp)*qai - (cons_alhl/cons_cp)*qal
1345  qai = 0.
1346  qal = 0.
1347  af = 0.
1348  end if
1349 
1350 end subroutine ls_cloud
1351 
1352 subroutine pdffrac (flag,qtmean,sigmaqt1,sigmaqt2,qstar,clfrac)
1353 
1354  IMPLICIT NONE
1355 
1356  !Inputs
1357  INTEGER, INTENT(IN) :: flag
1358  real(8), INTENT(IN) :: qtmean, sigmaqt1, sigmaqt2, qstar
1359 
1360  !Prognostic
1361  real(8), INTENT(INOUT) :: clfrac
1362 
1363  !LOCALS
1364  REAL(8) :: qtmode, qtmin, qtmax
1365  REAL(8) :: rh, rhd, q1, q2
1366 
1367 
1368  if (flag.eq.1) then !Tophat PDF
1369 
1370  if ((qtmean+sigmaqt1).lt.qstar) then
1371  clfrac = 0.
1372  else
1373  if (sigmaqt1.gt.0.) then
1374  clfrac = min((qtmean + sigmaqt1 - qstar),2.*sigmaqt1)/(2.*sigmaqt1)
1375  else
1376  clfrac = 1.
1377  endif
1378  endif
1379 
1380  elseif(flag.eq.2) then !Triangular PDF
1381 
1382  qtmode = qtmean + (sigmaqt1-sigmaqt2)/3.
1383  qtmin = min(qtmode-sigmaqt1,0.)
1384  qtmax = qtmode + sigmaqt2
1385 
1386  if (qtmax.lt.qstar) then
1387  clfrac = 0.
1388  elseif ( (qtmode.le.qstar).and.(qstar.lt.qtmax) ) then
1389  clfrac = (qtmax-qstar)*(qtmax-qstar) / ( (qtmax-qtmin)*(qtmax-qtmode) )
1390  elseif ( (qtmin.le.qstar).and.(qstar.lt.qtmode) ) then
1391  clfrac = 1. - ( (qstar-qtmin)*(qstar-qtmin) / ( (qtmax-qtmin)*(qtmode-qtmin) ) )
1392  elseif ( qstar.le.qtmin ) then
1393  clfrac = 1.
1394  endif
1395 
1396  elseif(flag.eq.3) then !TANH function for the perturabtions
1397 
1398  !Tophat PDF for the reference part
1399  if ((qtmean+sigmaqt1).lt.qstar) then
1400  clfrac = 0.
1401  else
1402  if (sigmaqt1.gt.0.) then
1403  clfrac = min((qtmean + sigmaqt1 - qstar),2.*sigmaqt1)/(2.*sigmaqt1)
1404  else
1405  clfrac = 1.
1406  endif
1407  endif
1408 
1409  elseif(flag.eq.4) then !Linear function for the perturabtions
1410 
1411  !Tophat PDF for the reference part
1412  if ((qtmean+sigmaqt1).lt.qstar) then
1413  clfrac = 0.
1414  else
1415  if (sigmaqt1.gt.0.) then
1416  clfrac = min((qtmean + sigmaqt1 - qstar),2.*sigmaqt1)/(2.*sigmaqt1)
1417  else
1418  clfrac = 1.
1419  endif
1420  endif
1421 
1422  endif
1423 
1424 end subroutine pdffrac
1425 
1426 
1427 subroutine pdfcondensate (flag,qtmean4,sigmaqt14,sigmaqt24,qstar4,condensate4)
1428 
1429  IMPLICIT NONE
1430 
1431  !Inputs
1432  INTEGER, INTENT(IN) :: flag
1433  real(8), INTENT(IN) :: qtmean4, sigmaqt14, sigmaqt24, qstar4
1434 
1435  !Prognostic
1436  real(8), INTENT(INOUT) :: condensate4
1437 
1438  !Locals
1439  real(8) :: qtmode, qtmin, qtmax, consta, constb, cloudf
1440  real(8) :: term1, term2, term3
1441  real(8) :: qtmean, sigmaqt1, sigmaqt2, qstar, condensate
1442 
1443  qtmean = qtmean4
1444  sigmaqt1 = sigmaqt14
1445  sigmaqt2 = sigmaqt24
1446  qstar = qstar4
1447 
1448  if (flag.eq.1) then
1449 
1450  if (qtmean+sigmaqt1.lt.qstar) then
1451  condensate = 0.d0
1452  elseif (qstar.gt.qtmean-sigmaqt1) then
1453  if (sigmaqt1.gt.0.d0) then
1454  condensate = (min(qtmean + sigmaqt1 - qstar,2.d0*sigmaqt1)**2)/ (4.d0*sigmaqt1)
1455  else
1456  condensate = qtmean-qstar
1457  endif
1458  else
1459  condensate = qtmean-qstar
1460  endif
1461 
1462  elseif (flag.eq.2) then
1463 
1464  qtmode = qtmean + (sigmaqt1-sigmaqt2)/3.d0
1465  qtmin = min(qtmode-sigmaqt1,0.d0)
1466  qtmax = qtmode + sigmaqt2
1467 
1468  if ( qtmax.lt.qstar ) then
1469  condensate = 0.d0
1470  elseif ( (qtmode.le.qstar).and.(qstar.lt.qtmax) ) then
1471  constb = 2.d0 / ( (qtmax - qtmin)*(qtmax-qtmode) )
1472  cloudf = (qtmax-qstar)*(qtmax-qstar) / ( (qtmax-qtmin)*(qtmax-qtmode) )
1473  term1 = (qstar*qstar*qstar)/3.d0
1474  term2 = (qtmax*qstar*qstar)/2.d0
1475  term3 = (qtmax*qtmax*qtmax)/6.d0
1476  condensate = constb * (term1-term2+term3) - qstar*cloudf
1477  elseif ( (qtmin.le.qstar).and.(qstar.lt.qtmode) ) then
1478  consta = 2.d0 / ( (qtmax - qtmin)*(qtmode-qtmin) )
1479  cloudf = 1.d0 - ( (qstar-qtmin)*(qstar-qtmin) / ( (qtmax-qtmin)*(qtmode-qtmin) ) )
1480  term1 = (qstar*qstar*qstar)/3.d0
1481  term2 = (qtmin*qstar*qstar)/2.d0
1482  term3 = (qtmin*qtmin*qtmin)/6.d0
1483  condensate = qtmean - ( consta * (term1-term2+term3) ) - qstar*cloudf
1484  elseif ( qstar.le.qtmin ) then
1485  condensate = qtmean-qstar
1486  endif
1487 
1488  elseif (flag.eq.3) then
1489 
1490  !Reference part done normally
1491  !IF (qtmean + sigmaqt1 .LT. qstar) THEN
1492  ! condensate = 0.d0
1493  !ELSE IF (qstar .GT. qtmean - sigmaqt1) THEN
1494  ! IF (sigmaqt1 .GT. 0.d0) THEN
1495  ! IF (qtmean + sigmaqt1 - qstar .GT. 2.d0*sigmaqt1) THEN
1496  ! min1 = 2.d0*sigmaqt1
1497  ! ELSE
1498  ! min1 = qtmean + sigmaqt1 - qstar
1499  ! END IF
1500  ! condensate = min1**2/(4.d0*sigmaqt1)
1501  ! ELSE
1502  ! condensate = qtmean - qstar
1503  ! END IF
1504  !ELSE
1505  ! condensate = qtmean - qstar
1506  !END IF
1507 
1508  !Perturbation part from linear
1509  if (qtmean - qstar > -0.5e-3) then
1510  condensate = qtmean - qstar
1511  else
1512  condensate = 0.0
1513  endif
1514 
1515  endif
1516 
1517  condensate4 = condensate
1518 
1519 end subroutine pdfcondensate
1520 
1521 
1522 
1523 subroutine evap_cnv( DT, RHCR, PL, TE, QV, QL, QI, F, XF, QS, RHO_W, CLD_EVP_EFF, &
1524  CONS_H2OMW, CONS_AIRMW, CONS_ALHL, CONS_RVAP, CONS_RGAS, CONS_PI, CONS_CP)
1525 
1526 IMPLICIT NONE
1527 
1528 !Inputs
1529 real(8), intent(in) :: dt, rhcr, pl, xf, qs, rho_w, cld_evp_eff
1530 real(8), intent(in) :: cons_h2omw, cons_airmw, cons_alhl, cons_rvap, cons_rgas, cons_pi, cons_cp
1531 
1532 !Prognostics
1533 real(8), intent(inout) :: te, qv, ql, qi, f
1534 
1535 !Locals
1536 real(8) :: es, radius, k1, k2, teff, qcm, evap, rhx, qc, a_eff, epsilon
1537 
1538 real(8), parameter :: k_cond = 2.4e-2
1539 real(8), parameter :: diffu = 2.2e-5
1540 real(8), parameter :: nn = 50.*1.0e6
1541 
1542  epsilon = cons_h2omw/cons_airmw
1543  a_eff = cld_evp_eff
1544 
1545  !EVAPORATION OF CLOUD WATER.
1546  es = 100.* pl * qs / ( (epsilon) + (1.0-(epsilon))*qs ) ! (100 <-^ convert from mbar to Pa)
1547 
1548  rhx = min( qv/qs , 1.00 )
1549 
1550  k1 = (cons_alhl**2) * rho_w / ( k_cond*cons_rvap*(te**2))
1551  k2 = cons_rvap * te * rho_w / ( diffu * (1000./pl) * es )
1552  !Here DIFFU is given for 1000 mb so 1000./PR accounts for increased diffusivity at lower pressure.
1553 
1554  if ( ( f > 0.) .and. ( ql > 0. ) ) then
1555  qcm=ql/f
1556  else
1557  qcm=0.
1558  end if
1559 
1560  call ldradius(pl,te,qcm,nn,rho_w,radius,cons_rgas,cons_pi)
1561 
1562  if ( (rhx < rhcr ) .and.(radius > 0.0) ) then
1563  teff = (rhcr - rhx) / ((k1+k2)*radius**2) ! / (1.00 - RHx)
1564  else
1565  teff = 0.0 ! -999.
1566  end if
1567 
1568  evap = a_eff*ql*dt*teff
1569  evap = min( evap , ql )
1570 
1571  qc=ql+qi
1572  if (qc > 0.) then
1573  f = f * ( qc - evap ) / qc
1574  end if
1575 
1576  qv = qv + evap
1577  ql = ql - evap
1578  te = te - (cons_alhl/cons_cp)*evap
1579 
1580 end subroutine evap_cnv
1581 
1582 subroutine subl_cnv( DT, RHCR, PL, TE, QV, QL, QI, F, XF, QS, RHO_W, CLD_EVP_EFF, &
1583  CONS_H2OMW, CONS_AIRMW, CONS_ALHL, CONS_RVAP, CONS_RGAS, CONS_PI, CONS_CP, CONS_ALHS)
1584 
1585 IMPLICIT NONE
1586 
1587 !INPUTS
1588 real(8), intent(in) :: dt, rhcr, pl, xf, qs, rho_w, cld_evp_eff
1589 real(8), intent(in) :: cons_h2omw, cons_airmw, cons_alhl, cons_rvap, cons_rgas, cons_pi, cons_cp, cons_alhs
1590 
1591 !PROGNOSTIC
1592 real(8), intent(inout) :: te, qv, ql, qi, f
1593 
1594 !LOCALS
1595 real(8) :: es, radius, k1, k2, teff, qcm, subl, rhx, qc, a_eff, nn, epsilon
1596 
1597 real(8), parameter :: k_cond = 2.4e-2
1598 real(8), parameter :: diffu = 2.2e-5
1599 
1600  epsilon = cons_h2omw/cons_airmw
1601 
1602  a_eff = cld_evp_eff
1603 
1604  nn = 5.*1.0e6
1605 
1606  es = 100.* pl * qs / ( (epsilon) + (1.0-(epsilon))*qs ) ! (100 s <-^ convert from mbar to Pa)
1607 
1608  rhx = min( qv/qs , 1.00 )
1609 
1610  k1 = (cons_alhl**2) * rho_w / ( k_cond*cons_rvap*(te**2))
1611  k2 = cons_rvap * te * rho_w / ( diffu * (1000./pl) * es )
1612 
1613  !Here DIFFU is given for 1000 mb so 1000./PR accounts for increased diffusivity at lower pressure.
1614  if ( ( f > 0.) .and. ( qi > 0. ) ) then
1615  qcm=qi/f
1616  else
1617  qcm=0.
1618  end if
1619 
1620  call ldradius(pl,te,qcm,nn,rho_w,radius,cons_rgas,cons_pi)
1621 
1622  if ( (rhx < rhcr ) .and.(radius > 0.0) ) then
1623  teff = ( rhcr - rhx) / ((k1+k2)*radius**2) ! / (1.00 - RHx)
1624  else
1625  teff = 0.0 ! -999.
1626  end if
1627 
1628  subl = a_eff*qi*dt*teff
1629  subl = min( subl , qi )
1630 
1631  qc=ql+qi
1632  if (qc > 0.) then
1633  f = f * ( qc - subl ) / qc
1634  end if
1635 
1636  qv = qv + subl
1637  qi = qi - subl
1638  te = te - (cons_alhs/cons_cp)*subl
1639 
1640 end subroutine subl_cnv
1641 
1642 subroutine ldradius(PL,TE,QCL,NN,RHO_W,RADIUS,CONS_RGAS,CONS_PI)
1643 
1644 IMPLICIT NONE
1645 
1646 !Inputs
1647 real(8), intent(in) :: te,pl,nn,qcl,rho_w
1648 real(8), intent(in) :: cons_rgas,cons_pi
1649 
1650 !Outputs
1651 real(8), intent(out) :: radius
1652 
1653 !Equiv. Spherical Cloud Particle Radius in m
1654 radius = ((qcl * (100.*pl / (cons_rgas*te )))/(nn*rho_w*(4./3.)*cons_pi))**(1./3.)
1655 
1656 end subroutine ldradius
1657 
1658 subroutine autoconversion_ls( DT, QC, QP, TE, PL, F, SUNDQV2, SUNDQV3, SUNDQT1, &
1659  C_00, LWCRIT, DZET)
1660 
1661 IMPLICIT NONE
1662 
1663 !Inputs
1664 real(8), intent(in) :: dt, te, pl, dzet, sundqv2, sundqv3, sundqt1, c_00, lwcrit
1665 
1666 !Prognostic
1667 real(8), intent(inout) :: qc, qp, f
1668 
1669 !Locals
1670 real(8) :: acf0, acf, c00x, iqccrx, f2, f3, rate, dqp, qcm, dqfac
1671 
1672  !Zero Locals
1673  acf0 = 0.0
1674  acf = 0.0
1675  c00x = 0.0
1676  iqccrx = 0.0
1677  f2 = 0.0
1678  f3 = 0.0
1679  rate = 0.0
1680  dqp = 0.0
1681  qcm = 0.0
1682  dqfac = 0.0
1683 
1684  CALL cons_sundq3(te, sundqv2, sundqv3, sundqt1, f2, f3 )
1685 
1686  c00x = c_00 * f2 * f3
1687  iqccrx = f2 * f3 / lwcrit
1688 
1689  if ( ( f > 0.) .and. ( qc > 0. ) )then
1690  qcm = qc/f
1691  else
1692  qcm = 0.
1693  end if
1694 
1695  rate = c00x * ( 1.0 - exp( - ( qcm * iqccrx )**2 ) )
1696 
1697  !Temporary kluge until we can figure a better to make thicker low clouds.
1698  f2 = 1.0
1699  f3 = 1.0
1700 
1701  !Implement ramps for gradual change in autoconv
1702 
1703  !Thicken low high lat clouds
1704  if ( pl .GE. 775. .AND. te .LE. 275. ) then
1705  !F3 = max(-0.016 * PL + 13.4, 0.2)
1706  f3 = 0.2
1707  end if
1708  if ( pl .GE. 825. .AND. te .LE. 282. ) then
1709  !F3 = max(0.11 * TE - 30.02, 0.2)
1710  f3 = 0.2
1711  end if
1712  if ( pl .GE. 775. .AND. pl .LT. 825. .AND. te .LE. 282. .AND. te .GT. 275.) then
1713  !F3 = min(max(-0.016*PL + 0.11 * TE - 16.85, 0.2),1.)
1714  f3 = 0.2
1715  end if
1716  if ( pl .GE. 825. .AND. te .LE. 275. ) then
1717  f3 = 0.2
1718  end if
1719  if ( pl .LE. 775. .OR. te .GT. 282. ) then
1720  f3 = 1.
1721  end if
1722 
1723  !Thin-out low tropical clouds
1724  if ( pl .GE. 950. .AND. te .GE. 285. ) then
1725  f3 = min(0.2 * te - 56, 2.)
1726  end if
1727  if ( pl .GE. 925. .AND. te .GE. 290. ) then
1728  f3 = min(0.04 * pl - 36., 2.)
1729  end if
1730  if ( pl .GE. 925. .AND. pl .LT. 950. .AND. te .GT. 285. .AND. te .LT. 290.) then
1731  f3 = max(min(0.04*pl + 0.2 * te - 94., 2.),1.)
1732  end if
1733  if ( pl .GE. 950. .AND. te .GE. 290. ) then
1734  f3 = 2.
1735  end if
1736 
1737  f3 = max( f3, 0.1 )
1738  rate = f3 * rate
1739 
1740  dqp = qc*( 1.0 - exp( -rate * dt ) )
1741  dqp = max( dqp , 0.0 ) ! Protects against floating point problems for tiny RATE
1742 
1743 
1744  !Wipe-out warm fogs
1745  dqfac = 0.
1746  if ( pl .GE. 975. .AND. te .GE. 280. ) then
1747  dqfac = max(min(0.2 * te - 56., 1.),0.)
1748  end if
1749  if ( pl .GE. 950. .AND. te .GE. 285. ) then
1750  dqfac = max(min(0.04 * pl - 38., 1.),0.)
1751  end if
1752  if ( pl .GE. 950. .AND. pl .LT. 975. .AND. te .GT. 280. .AND. te .LT. 285.) then
1753  dqfac = max(min(0.04*pl + 0.2 * te - 95., 1.),0.)
1754  end if
1755  if ( ( pl >= 975. ) .AND. (te >= 285. ) ) then
1756  dqfac = 1.
1757  end if
1758 
1759  dqp = max(dqp, dqfac*qc)
1760 
1761  qc = qc - dqp
1762  qp = qp + dqp
1763 
1764  !IF LARGE SCALE THEN
1765  if ( ((qc + dqp) > 0.) ) then
1766  f = qc * f / (qc + dqp )
1767  end if
1768 
1769 
1770 END SUBROUTINE autoconversion_ls
1771 
1772 subroutine autoconversion_cnv( DT, QC, QP, TE, PL, F, SUNDQV2, SUNDQV3, SUNDQT1, &
1773  C_00, LWCRIT, DZET)
1774 
1775 IMPLICIT NONE
1776 
1777 !Inputs
1778 real(8), intent(in) :: dt, te, pl, dzet, sundqv2, sundqv3, sundqt1, c_00, lwcrit
1779 
1780 !Prognostic
1781 real(8), intent(inout) :: qc, qp, f
1782 
1783 !Locals
1784 real(8) :: acf0, acf, c00x, iqccrx, f2, f3, rate, dqp, qcm, dqfac
1785 
1786  !Zero Locals
1787  acf0 = 0.0
1788  acf = 0.0
1789  c00x = 0.0
1790  iqccrx = 0.0
1791  f2 = 0.0
1792  f3 = 0.0
1793  rate = 0.0
1794  dqp = 0.0
1795  qcm = 0.0
1796  dqfac = 0.0
1797 
1798  CALL cons_sundq3(te, sundqv2, sundqv3, sundqt1, f2, f3 )
1799 
1800  c00x = c_00 * f2 * f3
1801  iqccrx = f2 * f3 / lwcrit
1802 
1803  if ( ( f > 0.) .and. ( qc > 0. ) )then
1804  qcm = qc/f
1805  else
1806  qcm = 0.
1807  end if
1808 
1809  rate = c00x * ( 1.0 - exp( - ( qcm * iqccrx )**2 ) )
1810 
1811  !Temporary kluge until we can figure a better to make thicker low clouds.
1812  f2 = 1.0
1813  f3 = 1.0
1814 
1815  !Implement ramps for gradual change in autoconv
1816 
1817  !Thicken low high lat clouds
1818  if ( pl .GE. 775. .AND. te .LE. 275. ) then
1819  !F3 = max(-0.016 * PL + 13.4, 0.2)
1820  f3 = 0.2
1821  end if
1822  if ( pl .GE. 825. .AND. te .LE. 282. ) then
1823  !F3 = max(0.11 * TE - 30.02, 0.2)
1824  f3 = 0.2
1825  end if
1826  if ( pl .GE. 775. .AND. pl .LT. 825. .AND. te .LE. 282. .AND. te .GT. 275.) then
1827  !F3 = min(max(-0.016*PL + 0.11 * TE - 16.85, 0.2),1.)
1828  f3 = 0.2
1829  end if
1830  if ( pl .GE. 825. .AND. te .LE. 275. ) then
1831  f3 = 0.2
1832  end if
1833  if ( pl .LE. 775. .OR. te .GT. 282. ) then
1834  f3 = 1.
1835  end if
1836 
1837  !Thin-out low tropical clouds
1838  if ( pl .GE. 950. .AND. te .GE. 285. ) then
1839  f3 = min(0.2 * te - 56, 2.)
1840  end if
1841  if ( pl .GE. 925. .AND. te .GE. 290. ) then
1842  f3 = min(0.04 * pl - 36., 2.)
1843  end if
1844  if ( pl .GE. 925. .AND. pl .LT. 950. .AND. te .GT. 285. .AND. te .LT. 290.) then
1845  f3 = max(min(0.04*pl + 0.2 * te - 94., 2.),1.)
1846  end if
1847  if ( pl .GE. 950. .AND. te .GE. 290. ) then
1848  f3 = 2.
1849  end if
1850 
1851  f3 = max( f3, 0.1 )
1852  rate = f3 * rate
1853 
1854  dqp = qc*( 1.0 - exp( -rate * dt ) )
1855  dqp = max( dqp , 0.0 ) ! Protects against floating point problems for tiny RATE
1856 
1857 
1858  !Wipe-out warm fogs
1859  dqfac = 0.
1860  if ( pl .GE. 975. .AND. te .GE. 280. ) then
1861  dqfac = max(min(0.2 * te - 56., 1.),0.)
1862  end if
1863  if ( pl .GE. 950. .AND. te .GE. 285. ) then
1864  dqfac = max(min(0.04 * pl - 38., 1.),0.)
1865  end if
1866  if ( pl .GE. 950. .AND. pl .LT. 975. .AND. te .GT. 280. .AND. te .LT. 285.) then
1867  dqfac = max(min(0.04*pl + 0.2 * te - 95., 1.),0.)
1868  end if
1869  if ( ( pl >= 975. ) .AND. (te >= 285. ) ) then
1870  dqfac = 1.
1871  end if
1872 
1873  dqp = max(dqp, dqfac*qc)
1874 
1875  qc = qc - dqp
1876  qp = qp + dqp
1877 
1878 END SUBROUTINE autoconversion_cnv
1879 
1880 subroutine get_ice_fraction(TEMP, T_ICE_ALL, T_ICE_MAX, ICEFRPWR, ICEFRCT)
1881 
1882 IMPLICIT NONE
1883 
1884 !Inputs
1885 real(8), intent(in) :: temp, t_ice_all, t_ice_max
1886 integer, intent(in) :: icefrpwr
1887 
1888 !Outputs
1889 real(8), intent(out) :: icefrct
1890 
1891 
1892  icefrct = 0.00
1893  if ( temp <= t_ice_all ) then
1894  icefrct = 1.000
1895  else if ( (temp > t_ice_all) .AND. (temp <= t_ice_max) ) then
1896  icefrct = 1.00 - ( temp - t_ice_all ) / ( t_ice_max - t_ice_all )
1897  end if
1898 
1899  icefrct = min(icefrct,1.00)
1900  icefrct = max(icefrct,0.00)
1901 
1902  icefrct = icefrct**icefrpwr
1903 
1904 
1905 end subroutine get_ice_fraction
1906 
1907 
1908 SUBROUTINE cons_sundq3(TEMP,RATE2,RATE3,TE1, F2, F3)
1909 
1910 IMPLICIT NONE
1911 
1912 !Inputs
1913 real(8), intent(in) :: rate2, rate3, te1, temp
1914 
1915 !Outputs
1916 real(8), intent(out) :: f2, f3
1917 
1918 !Locals
1919 real(8), parameter :: te0 = 273.
1920 real(8), parameter :: te2 = 200.
1921 real(8) :: jump1
1922 
1923 
1924  jump1= (rate2-1.0) / ( ( te0-te1 )**0.333 )
1925 
1926  !Ice - phase treatment
1927  IF ( temp .GE. te0 ) THEN
1928  f2 = 1.0
1929  f3 = 1.0
1930  END IF
1931  IF ( ( temp .GE. te1 ) .AND. ( temp .LT. te0 ) ) THEN
1932  if (abs(te0 - temp) .gt. 0.0) then !Linearisation security
1933  f2 = 1.0 + jump1 * (( te0 - temp )**0.3333)
1934  else
1935  f2 = 1.0
1936  endif
1937  f3 = 1.0
1938  END IF
1939  IF ( temp .LT. te1 ) THEN
1940  f2 = rate2 + (rate3-rate2)*(te1-temp)/(te1-te2)
1941  f3 = 1.0
1942  END IF
1943 
1944  f2 = min(f2,27.0)
1945 
1946 
1947 end subroutine cons_sundq3
1948 
1949 
1950 subroutine cons_microphys(TEMP,PR,Q_SAT,AA,BB,CONS_H2OMW,CONS_AIRMW,CONS_RVAP,ALHX3)
1951 
1952 IMPLICIT NONE
1953 
1954 !Inputs
1955 real(8), intent(in) :: temp, q_sat, pr, alhx3
1956 real(8), intent(in) :: cons_h2omw, cons_airmw, cons_rvap
1957 
1958 !Outputs
1959 real(8), intent(out) :: aa, bb
1960 
1961 !Locals
1962 real(8), parameter :: k_cond = 2.4e-2
1963 real(8), parameter :: diffu = 2.2e-5
1964 real(8) :: e_sat, epsi
1965 
1966  epsi = cons_h2omw/cons_airmw
1967 
1968  e_sat = 100.* pr * q_sat /( (epsi) + (1.0-(epsi))*q_sat ) ! (100 converts from mbar to Pa)
1969 
1970  aa = ( alhx3**2 ) / ( k_cond*cons_rvap*(temp**2) )
1971  bb = cons_rvap*temp / ( diffu*(1000./pr)*e_sat )
1972 
1973 end subroutine cons_microphys
1974 
1975 
1976 
1977 subroutine cons_alhx(T,ALHX3,T_ICE_MAX,T_ICE_ALL,CONS_ALHS,CONS_ALHL)
1978 
1979 IMPLICIT NONE
1980 
1981 !Inputs
1982 real(8), intent(in) :: t, t_ice_max, t_ice_all
1983 real(8), intent(in) :: cons_alhs, cons_alhl
1984 
1985 !Outputs
1986 real(8), intent(out) :: alhx3
1987 
1988  if ( t < t_ice_all ) then
1989  alhx3 = cons_alhs
1990  end if
1991 
1992  if ( t > t_ice_max ) then
1993  alhx3 = cons_alhl
1994  end if
1995 
1996  if ( (t <= t_ice_max) .and. (t >= t_ice_all) ) then
1997  alhx3 = cons_alhs + (cons_alhl - cons_alhs)*( t - t_ice_all ) /( t_ice_max - t_ice_all )
1998  end if
1999 
2000 end subroutine cons_alhx
2001 
2002 subroutine marshpalm(RAIN,PR,DIAM3,NTOTAL,W,VE)
2003 
2004 Implicit None
2005 
2006 !Inputs
2007 real(8), intent(in ) :: rain, pr ! in kg m^-2 s^-1, mbar
2008 
2009 !Outputs
2010 real(8), intent(out) :: diam3, ntotal, w, ve
2011 
2012 !Locals
2013 integer :: iqd
2014 real(8), parameter :: n0 = 0.08 !cm^-3
2015 real(8) :: rain_day,slopr,diam1
2016 
2017 real(8) :: rx(8) , d3x(8)
2018 
2019  rain_day = 0.0
2020  slopr = 0.0
2021  diam1 = 0.0
2022 
2023  !Marshall-Palmer sizes at different rain-rates: avg(D^3)
2024  !RX = (/ 0. , 5. , 20. , 80. , 320. , 1280., 5120., 20480. /) ! rain per in mm/day
2025  rx(1) = 0.
2026  rx(2) = 5.
2027  rx(3) = 20.
2028  rx(4) = 80.
2029  rx(5) = 320.
2030  rx(6) = 1280.
2031  rx(7) = 5120.
2032  rx(8) = 20480.
2033 
2034  !D3X= (/ 0.019, 0.032, 0.043, 0.057, 0.076, 0.102, 0.137, 0.183 /)
2035  d3x(1) = 0.019
2036  d3x(2) = 0.032
2037  d3x(3) = 0.043
2038  d3x(4) = 0.057
2039  d3x(5) = 0.076
2040  d3x(6) = 0.102
2041  d3x(7) = 0.137
2042  d3x(8) = 0.183
2043 
2044  rain_day = rain * 3600. *24.
2045 
2046  IF ( rain_day <= 0.00 ) THEN
2047  diam1 = 0.00
2048  diam3 = 0.00
2049  ntotal= 0.00
2050  w = 0.00
2051  END IF
2052 
2053  DO iqd = 1,7
2054  IF ( (rain_day <= rx(iqd+1)) .AND. (rain_day > rx(iqd) ) ) THEN
2055  slopr =( d3x(iqd+1)-d3x(iqd) ) / ( rx(iqd+1)-rx(iqd) )
2056  diam3 = d3x(iqd) + (rain_day-rx(iqd))*slopr
2057  END IF
2058  END DO
2059 
2060  IF ( rain_day >= rx(8) ) THEN
2061  diam3=d3x(8)
2062  END IF
2063 
2064  ntotal = 0.019*diam3
2065 
2066  diam3 = 0.664 * diam3
2067 
2068  w = (2483.8 * diam3 + 80.)*sqrt(1000./pr)
2069 
2070  ve = max( 0.99*w/100. , 1.000 )
2071 
2072  diam1 = 3.0*diam3
2073 
2074  diam1 = diam1/100.
2075  diam3 = diam3/100.
2076  w = w/100.
2077  ntotal = ntotal*1.0e6
2078 
2079 end subroutine marshpalm
2080 
2081 subroutine ice_settlefall_cnv( WXR, QI, PL, TE, F, CONS_RGAS, KHu, KHl, k, DT, DZ, QP, ANV_ICEFALL_C )
2082 
2083 IMPLICIT NONE
2084 
2085 !Inputs
2086 real(8), intent(in) :: wxr, pl, te, dz, dt, anv_icefall_c
2087 real(8), intent(in) :: cons_rgas
2088 integer, intent(in) :: khu, khl, k
2089 
2090 real(8), intent(inout) :: qi, f, qp
2091 
2092 !Locals
2093 real(8) :: rho, xim, lxim, qixp, vf
2094 
2095 
2096  rho = 1000.*100.*pl/(cons_rgas*te) ! 1000 TAKES TO g m^-3 ; 100 takes mb TO Pa
2097 
2098  if ( ( f > 0.) .and. ( qi > 0. ) ) then
2099  xim = (qi/f)*rho
2100  else
2101  xim = 0.
2102  end if
2103 
2104  if ( xim > 0.) then
2105  lxim = log10(xim)
2106  else
2107  lxim = 0.0
2108  end if
2109 
2110  vf = 128.6 + 53.2*lxim + 5.5*lxim**2
2111 
2112  !VF = VF*100./MAX(PL,10.) ! Reduce/increase fall speeds for high/low pressure (NOT in LC98!!! )
2113  ! Assume unmodified they represent situation at 100 mb
2114 
2115  if (wxr > 0.) then
2116  vf = vf * ( 100./max(pl,10.) )**wxr
2117  endif
2118 
2119  vf = vf/100.
2120 
2121  if (khu .gt. 0 .and. khl .gt. 0) then
2122  if ( (k-1 .ge. khu) .and. (k-1 .le. khl) ) then
2123  vf = 0.01 * vf
2124  end if
2125  endif
2126 
2127  vf = anv_icefall_c * vf
2128 
2129  qixp = 0.0
2130 
2131  qixp = qi * ( vf * dt / dz )
2132  qixp = min( qixp , qi )
2133 
2134  qixp = max( qixp, 0.0 ) ! protects against precision problem
2135 
2136  qp = qp + qixp
2137  qi = qi - qixp
2138 
2139 
2140 
2141 end SUBROUTINE ice_settlefall_cnv
2142 
2143 
2144 subroutine ice_settlefall_ls( WXR, QI, PL, TE, F, CONS_RGAS, KHu, KHl, k, DT, DZ, QP, LS_ICEFALL_C )
2145 
2146 IMPLICIT NONE
2147 
2148 !Inputs
2149 real(8), intent(in) :: wxr, pl, te, dz, dt, ls_icefall_c
2150 real(8), intent(in) :: cons_rgas
2151 integer, intent(in) :: khu, khl, k
2152 
2153 real(8), intent(inout) :: qi, f, qp
2154 
2155 !Locals
2156 real(8) :: rho, xim, lxim, qixp, vf
2157 
2158  rho = 1000.*100.*pl/(cons_rgas*te) ! 1000 TAKES TO g m^-3 ; 100 takes mb TO Pa
2159 
2160  if ( ( f > 0.) .and. ( qi > 0. ) ) then
2161  xim = (qi/f)*rho
2162  else
2163  xim = 0.
2164  end if
2165 
2166  if ( xim > 0.) then
2167  lxim = log10(xim)
2168  else
2169  lxim = 0.0
2170  end if
2171 
2172  if (abs(xim) .gt. 0.0) then !Linearisation security
2173  vf = 109.0*(xim**0.16)
2174  else
2175  vf = 0.0
2176  endif
2177 
2178  !VF = VF*100./MAX(PL,10.) ! Reduce/increase fall speeds for high/low pressure (NOT in LC98!!! )
2179  ! Assume unmodified they represent situation at 100 mb
2180 
2181  if (wxr > 0.) then
2182  vf = vf * ( 100./max(pl,10.) )**wxr
2183  endif
2184 
2185  vf = vf/100.
2186 
2187  if (khu .gt. 0 .and. khl .gt. 0) then
2188  if ( (k-1 .ge. khu) .and. (k-1 .le. khl) ) then
2189  vf = 0.01 * vf
2190  end if
2191  endif
2192 
2193  vf = ls_icefall_c * vf
2194 
2195  qixp = 0.0
2196 
2197  qixp = qi * ( vf * dt / dz )
2198  qixp = min( qixp , qi )
2199 
2200  qixp = max( qixp, 0.0 ) ! protects against precision problem
2201 
2202  qp = qp + qixp
2203  qi = qi - qixp
2204 
2205  if ( ((qi + qixp) > 0.) ) then
2206  f = qi * f / (qi + qixp )
2207  end if
2208 
2209 
2210 end SUBROUTINE ice_settlefall_ls
2211 
2212 
2213 subroutine precipandevap(K, KTOP, LM, DT, FRLAND, RHCR3, QPl, QPi, QCl, QCi, TE, QV, mass, imass, &
2214  PL, dZE, QDDF3, AA, BB, AREA, PFl_above_in, PFl_above_out, PFi_above_in, PFi_above_out, &
2215  EVAP_DD_above_in, EVAP_DD_above_out, SUBL_DD_above_in, SUBL_DD_above_out, &
2216  ENVFC, DDRFC, &
2217  CONS_ALHF, CONS_ALHS, CONS_ALHL, CONS_CP, CONS_TICE,CONS_H2OMW,CONS_AIRMW, REVAP_OFF_P, &
2218  C_ACC, C_EV_R, C_EV_S, RHO_W, ESTBLX )
2219 
2220 IMPLICIT NONE
2221 
2222 !Inputs
2223 integer, intent(in) :: k, lm, ktop
2224 real(8), intent(in) :: dt, mass, imass, pl, aa, bb, rhcr3, dze, qddf3, area, frland, envfc, ddrfc
2225 real(8), intent(in) :: cons_alhf, cons_alhs, cons_alhl, cons_cp, cons_tice, cons_h2omw, cons_airmw
2226 real(8), intent(in) :: revap_off_p
2227 real(8), intent(in) :: c_acc, c_ev_r, c_ev_s, rho_w
2228 real(8), intent(in) :: estblx(:)
2229 
2230 !Prognostics
2231 real(8), intent(inout) :: qv, qpl, qpi, qcl, qci, te
2232 real(8), intent(inout) :: pfl_above_in, pfl_above_out, pfi_above_in, pfi_above_out
2233 real(8), intent(inout) :: evap_dd_above_in, evap_dd_above_out, subl_dd_above_in, subl_dd_above_out
2234 
2235 !Locals
2236 integer :: ns, nsmx, itr,l
2237 
2238 real(8) :: pfi, pfl, qs, dqs, envfrac, tko, qko, qstko, dqstko, rh_box, t_ed, qplko, qpiko
2239 real(8) :: ifactor, rainrat0, snowrat0, fallrn, fallsn, vesn, vern, nrain, nsnow, efactor
2240 real(8) :: tinlayerrn, diamrn, droprad, tinlayersn, diamsn, flakrad
2241 real(8) :: evap, subl, accr, mltfrz, evapx, sublx, evap_dd,subl_dd,ddfract, landseaf
2242 real(8) :: tau_frz, tau_mlt
2243 
2244 real(8), parameter :: trmv_l = 1.0 !m/s
2245 logical, parameter :: taneff = .false.
2246 
2247 !Fraction of precip falling through "environment" vs through cloud
2248 real(8), parameter :: b_sub = 1.00
2249 
2250  envfrac = envfc
2251 
2252  IF ( area > 0. ) THEN
2253  ifactor = 1./ ( area )
2254  ELSE
2255  ifactor = 1.00
2256  END if
2257 
2258  ifactor = max( ifactor, 1.) !
2259 
2260  !Start at top of precip column:
2261  !
2262  ! a) Accrete
2263  ! b) Evaporate/Sublimate
2264  ! c) Rain/Snow-out to next level down
2265  ! d) return to (a)
2266 
2267  !Update saturated humidity
2268  call dqsats_bac(dqs,qs,te,pl,estblx,cons_h2omw,cons_airmw)
2269 
2270  ddfract = ddrfc
2271 
2272  IF (k == ktop) THEN
2273 
2274  pfl=qpl*mass
2275  pfi=qpi*mass
2276 
2277  evap_dd = 0.
2278  subl_dd = 0.
2279 
2280  ELSE
2281 
2282  qpl = qpl + pfl_above_in * imass
2283  pfl = 0.00
2284 
2285  qpi = qpi + pfi_above_in * imass
2286  pfi = 0.00
2287 
2288  accr = b_sub * c_acc * ( qpl*mass ) *qcl
2289 
2290  accr = min( accr , qcl )
2291 
2292  qpl = qpl + accr
2293  qcl = qcl - accr
2294 
2295  !Accretion of liquid condensate by falling ice/snow
2296  accr = b_sub * c_acc * ( qpi*mass ) *qcl
2297 
2298  accr = min( accr , qcl )
2299 
2300  qpi = qpi + accr
2301  qcl = qcl - accr
2302 
2303  !! Liquid freezes when accreted by snow
2304  te = te + cons_alhf*accr/cons_cp
2305 
2306  rainrat0 = ifactor*qpl*mass/dt
2307  snowrat0 = ifactor*qpi*mass/dt
2308 
2309  call marshpalm(rainrat0,pl,diamrn,nrain,fallrn,vern)
2310  call marshpalm(snowrat0,pl,diamsn,nsnow,fallsn,vesn)
2311 
2312  tinlayerrn = dze / ( fallrn+0.01 )
2313  tinlayersn = dze / ( fallsn+0.01 )
2314 
2315  !Melting of Frozen precipitation
2316  tau_frz = 5000. ! time scale for freezing (s).
2317 
2318  mltfrz = 0.0
2319  IF ( (te > cons_tice ) .and.(te <= cons_tice+5. ) ) THEN
2320  mltfrz = tinlayersn * qpi *( te - cons_tice ) / tau_frz
2321  mltfrz = min( qpi , mltfrz )
2322  te = te - cons_alhf*mltfrz/cons_cp
2323  qpl = qpl + mltfrz
2324  qpi = qpi - mltfrz
2325  END IF
2326 
2327  mltfrz = 0.0
2328  IF ( te > cons_tice+5. ) THEN ! Go Ahead and melt any snow/hail left above 5 C
2329  mltfrz= qpi
2330  te = te - cons_alhf*mltfrz/cons_cp
2331  qpl = qpl + mltfrz
2332  qpi = qpi - mltfrz
2333  END IF
2334 
2335  mltfrz = 0.0
2336  if ( k >= lm-1 ) THEN
2337  IF ( te > cons_tice+0. ) THEN ! Go Ahead and melt any snow/hail left above 0 C in lowest layers
2338  mltfrz= qpi
2339  te = te - cons_alhf*mltfrz/cons_cp
2340  qpl = qpl + mltfrz
2341  qpi = qpi - mltfrz
2342  END IF
2343  endif
2344 
2345  !Freezing of liquid precipitation
2346  mltfrz = 0.0
2347  IF ( te <= cons_tice ) THEN
2348  te = te + cons_alhf*qpl/cons_cp
2349  qpi = qpl + qpi
2350  mltfrz = qpl
2351  qpl = 0.
2352  END IF
2353 
2354  !In the exp below, evaporation time scale is determined "microphysically" from temp,
2355  !press, and drop size. In this context C_EV becomes a dimensionless fudge-fraction.
2356  !Also remember that these microphysics are still only for liquid.
2357 
2358  qko = qv
2359  tko = te
2360  qplko = qpl
2361  qpiko = qpi
2362 
2363  !do itr = 1,1
2364  itr = 1
2365 
2366  dqstko = dqs
2367  qstko = qs + dqstko * ( tko - te )
2368  qstko = max( qstko , 1.0e-7 )
2369 
2370  rh_box = qko/qstko
2371 
2372  qko = qv
2373  tko = te
2374 
2375  IF ( rh_box < rhcr3 ) THEN
2376  efactor = rho_w * ( aa + bb ) / (rhcr3 - rh_box )
2377  else
2378  efactor = 9.99e9
2379  end if
2380 
2381  if ( frland < 0.1 ) then
2382  landseaf = 0.5 ! Over Ocean
2383  else
2384  landseaf = 0.5 ! Over Land
2385  end if
2386 
2387  landseaf = 1.00
2388 
2389  !Rain falling
2390  if ( (rh_box < rhcr3) .AND. (diamrn > 0.00) .AND. (pl > 100.) .AND. (pl < revap_off_p) ) then
2391  droprad=0.5*diamrn
2392  t_ed = efactor * droprad**2
2393  t_ed = t_ed * ( 1.0 + dqstko*cons_alhl/cons_cp )
2394  evap = qpl*(1.0 - exp( -c_ev_r * vern * landseaf * envfrac * tinlayerrn / t_ed ) )
2395  ELSE
2396  evap = 0.0
2397  END if
2398 
2399  !Snow falling
2400  if ( (rh_box < rhcr3) .AND. (diamsn > 0.00) .AND. (pl > 100.) .AND. (pl < revap_off_p) ) then
2401  flakrad=0.5*diamsn
2402  t_ed = efactor * flakrad**2
2403  t_ed = t_ed * ( 1.0 + dqstko*cons_alhs/cons_cp )
2404  subl = qpi*(1.0 - exp( -c_ev_s * vesn * landseaf * envfrac * tinlayersn / t_ed ) )
2405  ELSE
2406  subl = 0.0
2407  END IF
2408 
2409  !if (itr == 1) then
2410  ! EVAPx = EVAP
2411  ! SUBLx = SUBL
2412  !else
2413  ! EVAP = (EVAP+EVAPx) /2.0
2414  ! SUBL = (SUBL+SUBLx) /2.0
2415  !endif
2416 
2417  qko = qv + evap + subl
2418  tko = te - evap * cons_alhl / cons_cp - subl * cons_alhs / cons_cp
2419 
2420  !enddo
2421 
2422  qpi = qpi - subl
2423  qpl = qpl - evap
2424 
2425  !Put some re-evap/re-subl precip in to a \quote{downdraft} to be applied later
2426  evap_dd = evap_dd_above_in + ddfract*evap*mass
2427  evap = evap - ddfract*evap
2428  subl_dd = subl_dd_above_in + ddfract*subl*mass
2429  subl = subl - ddfract*subl
2430 
2431  qv = qv + evap + subl
2432  te = te - evap * cons_alhl / cons_cp - subl * cons_alhs / cons_cp
2433 
2434  pfl = qpl*mass
2435  pfi = qpi*mass
2436 
2437  end if
2438 
2439  evap = qddf3*evap_dd/mass
2440  subl = qddf3*subl_dd/mass
2441  qv = qv + evap + subl
2442  te = te - evap * cons_alhl / cons_cp - subl * cons_alhs / cons_cp
2443 
2444  qpi = 0.
2445  qpl = 0.
2446 
2447  pfl_above_out = pfl
2448  pfi_above_out = pfi
2449 
2450  evap_dd_above_out = evap_dd
2451  subl_dd_above_out = subl_dd
2452 
2453 end subroutine precipandevap
2454 
2455 
2456 subroutine dqsat_bac(DQSi,QSSi,TEMP,PLO,im,jm,lm,ESTBLX,CONS_H2OMW,CONS_AIRMW)
2457 !COMPUTES SATURATION VAPOUR PRESSURE QSSi AND GRADIENT w.r.t TEMPERATURE DQSi.
2458 !INPUTS ARE TEMPERATURE AND PLO (PRESSURE AT T-LEVELS)
2459 !VALES ARE COMPUTED FROM LOOK-UP TALBE (PIECEWISE LINEAR)
2460 
2461  IMPLICIT NONE
2462 
2463  !Inputs
2464  integer :: im,jm,lm
2465  real(8), dimension(im,jm,lm) :: temp, plo
2466  real(8) :: estblx(:)
2467  real(8) :: cons_h2omw, cons_airmw
2468 
2469  !Outputs
2470  real(8), dimension(im,jm,lm) :: dqsi, qssi
2471 
2472  !Locals
2473  real(8), parameter :: max_mixing_ratio = 1.0
2474  real(8) :: esfac
2475 
2476  integer :: i, j, k
2477 
2478  real(8) :: tl, tt, ti, dqsat, qsat, dqq, qq, pl, pp, dd
2479  integer :: it
2480 
2481  integer, parameter :: degsubs = 100
2482  real(8), parameter :: tmintbl = 150.0, tmaxtbl = 333.0
2483  integer, parameter :: tablesize = nint(tmaxtbl-tmintbl)*degsubs + 1
2484 
2485  esfac = cons_h2omw/cons_airmw
2486 
2487  do k=1,lm
2488  do j=1,jm
2489  do i=1,im
2490 
2491  tl = temp(i,j,k)
2492  pl = plo(i,j,k)
2493 
2494  pp = pl*100.0
2495 
2496  if (tl<=tmintbl) then
2497  ti = tmintbl
2498  elseif(tl>=tmaxtbl-.001) then
2499  ti = tmaxtbl-.001
2500  else
2501  ti = tl
2502  end if
2503 
2504  tt = (ti - tmintbl)*degsubs+1
2505  it = int(tt)
2506 
2507  dqq = estblx(it+1) - estblx(it)
2508  qq = (tt-it)*dqq + estblx(it)
2509 
2510  if (pp <= qq) then
2511  qsat = max_mixing_ratio
2512  dqsat = 0.0
2513  else
2514  dd = 1.0/(pp - (1.0-esfac)*qq)
2515  qsat = esfac*qq*dd
2516  dqsat = (esfac*degsubs)*dqq*pp*(dd*dd)
2517  end if
2518 
2519  dqsi(i,j,k) = dqsat
2520  qssi(i,j,k) = qsat
2521 
2522  end do
2523  end do
2524  end do
2525 
2526 end subroutine dqsat_bac
2527 
2528 subroutine dqsats_bac(DQSi,QSSi,TEMP,PLO,ESTBLX,CONS_H2OMW,CONS_AIRMW)
2529 !COMPUTES SATURATION VAPOUR PRESSURE QSSi AND GRADIENT w.r.t TEMPERATURE DQSi.
2530 !INPUTS ARE TEMPERATURE AND PLO (PRESSURE AT T-LEVELS)
2531 !VALES ARE COMPUTED FROM LOOK-UP TALBE (PIECEWISE LINEAR)
2532 
2533  IMPLICIT NONE
2534 
2535  !Inputs
2536  real(8) :: temp, plo
2537  real(8) :: estblx(:)
2538  real(8) :: cons_h2omw,cons_airmw
2539 
2540  !Outputs
2541  real(8) :: dqsi, qssi
2542 
2543  !Locals
2544  real(8), parameter :: max_mixing_ratio = 1.0
2545  real(8) :: esfac
2546 
2547  real(8) :: tl, tt, ti, dqsat, qsat, dqq, qq, pl, pp, dd
2548  integer :: it
2549 
2550  integer, parameter :: degsubs = 100
2551  real(8), parameter :: tmintbl = 150.0, tmaxtbl = 333.0
2552  integer, parameter :: tablesize = nint(tmaxtbl-tmintbl)*degsubs + 1
2553 
2554  esfac = cons_h2omw/cons_airmw
2555 
2556  tl = temp
2557  pl = plo
2558 
2559  pp = pl*100.0
2560 
2561  if (tl<=tmintbl) then
2562  ti = tmintbl
2563  elseif(tl>=tmaxtbl-.001) then
2564  ti = tmaxtbl-.001
2565  else
2566  ti = tl
2567  end if
2568 
2569  tt = (ti - tmintbl)*degsubs+1
2570  it = int(tt)
2571 
2572  dqq = estblx(it+1) - estblx(it)
2573  qq = (tt-it)*dqq + estblx(it)
2574 
2575  if (pp <= qq) then
2576  qsat = max_mixing_ratio
2577  dqsat = 0.0
2578  else
2579  dd = 1.0/(pp - (1.0-esfac)*qq)
2580  qsat = esfac*qq*dd
2581  dqsat = (esfac*degsubs)*dqq*pp*(dd*dd)
2582  end if
2583 
2584  dqsi = dqsat
2585  qssi = qsat
2586 
2587 end subroutine dqsats_bac
2588 
2589 END MODULE cloud
cloud::autoconversion_cnv
subroutine, public autoconversion_cnv(DT, QC, QP, TE, PL, F, SUNDQV2, SUNDQV3, SUNDQT1, C_00, LWCRIT, DZET)
Definition: cloud.F90:1774
cloud::ldradius
subroutine, public ldradius(PL, TE, QCL, NN, RHO_W, RADIUS, CONS_RGAS, CONS_PI)
Definition: cloud.F90:1643
cloud::pdfcondensate
subroutine, public pdfcondensate(flag, qtmean4, sigmaqt14, sigmaqt24, qstar4, condensate4)
Definition: cloud.F90:1428
cloud::dqsat_bac
subroutine, public dqsat_bac(DQSi, QSSi, TEMP, PLO, im, jm, lm, ESTBLX, CONS_H2OMW, CONS_AIRMW)
Definition: cloud.F90:2457
cloud::precipandevap
subroutine, public precipandevap(K, KTOP, LM, DT, FRLAND, RHCR3, QPl, QPi, QCl, QCi, TE, QV, mass, imass, PL, dZE, QDDF3, AA, BB, AREA, PFl_above_in, PFl_above_out, PFi_above_in, PFi_above_out, EVAP_DD_above_in, EVAP_DD_above_out, SUBL_DD_above_in, SUBL_DD_above_out, ENVFC, DDRFC, CONS_ALHF, CONS_ALHS, CONS_ALHL, CONS_CP, CONS_TICE, CONS_H2OMW, CONS_AIRMW, REVAP_OFF_P, C_ACC, C_EV_R, C_EV_S, RHO_W, ESTBLX)
Definition: cloud.F90:2219
cloud::evap_cnv
subroutine, public evap_cnv(DT, RHCR, PL, TE, QV, QL, QI, F, XF, QS, RHO_W, CLD_EVP_EFF, CONS_H2OMW, CONS_AIRMW, CONS_ALHL, CONS_RVAP, CONS_RGAS, CONS_PI, CONS_CP)
Definition: cloud.F90:1525
cloud::ice_settlefall_ls
subroutine, public ice_settlefall_ls(WXR, QI, PL, TE, F, CONS_RGAS, KHu, KHl, k, DT, DZ, QP, LS_ICEFALL_C)
Definition: cloud.F90:2145
cloud::ice_settlefall_cnv
subroutine, public ice_settlefall_cnv(WXR, QI, PL, TE, F, CONS_RGAS, KHu, KHl, k, DT, DZ, QP, ANV_ICEFALL_C)
Definition: cloud.F90:2082
cloud::pdffrac
subroutine, public pdffrac(flag, qtmean, sigmaqt1, sigmaqt2, qstar, clfrac)
Definition: cloud.F90:1353
cloud
Definition: cloud.F90:1
cloud::cons_alhx
subroutine, public cons_alhx(T, ALHX3, T_ICE_MAX, T_ICE_ALL, CONS_ALHS, CONS_ALHL)
Definition: cloud.F90:1978
cloud::cloud_driver
subroutine, public cloud_driver(DT, IM, JM, LM, th, q, ple, CNV_DQLDT, CNV_MFD, CNV_PRC3, CNV_UPDF, QI_ls, QL_ls, QI_con, QL_con, CF_LS, CF_con, FRLAND, PHYSPARAMS, ESTBLX, KHu, KHl, CONS_RUNIV, CONS_KAPPA, CONS_AIRMW, CONS_H2OMW, CONS_GRAV, CONS_ALHL, CONS_ALHF, CONS_PI, CONS_RGAS, CONS_CP, CONS_VIREPS, CONS_ALHS, CONS_TICE, CONS_RVAP, CONS_P00, do_moist_physics)
Definition: cloud.F90:24
cloud::pdf_width
subroutine, public pdf_width(PP, FRLAND, maxrhcrit, maxrhcritland, turnrhcrit, minrhcrit, pi_0, ALPHA)
Definition: cloud.F90:1046
cloud::cons_microphys
subroutine, public cons_microphys(TEMP, PR, Q_SAT, AA, BB, CONS_H2OMW, CONS_AIRMW, CONS_RVAP, ALHX3)
Definition: cloud.F90:1951
cloud::meltfreeze
subroutine, public meltfreeze(DT, TE, QL, QI, T_ICE_ALL, T_ICE_MAX, ICEFRPWR, CONS_ALHL, CONS_ALHS, CONS_CP)
Definition: cloud.F90:923
cloud::convec_src
subroutine, public convec_src(DT, MASS, iMASS, TE, QV, DCF, DMF, QLA, QIA, AF, QS, CONS_ALHS, CONS_ALHL, CONS_CP, T_ICE_ALL, T_ICE_MAX, ICEFRPWR)
Definition: cloud.F90:969
cloud::cloud_tidy
subroutine, public cloud_tidy(QV, TE, QLC, QIC, CF, QLA, QIA, AF, CONS_ALHL, CONS_ALHS, CONS_CP)
Definition: cloud.F90:852
cloud::cons_sundq3
subroutine, public cons_sundq3(TEMP, RATE2, RATE3, TE1, F2, F3)
Definition: cloud.F90:1909
cloud::subl_cnv
subroutine, public subl_cnv(DT, RHCR, PL, TE, QV, QL, QI, F, XF, QS, RHO_W, CLD_EVP_EFF, CONS_H2OMW, CONS_AIRMW, CONS_ALHL, CONS_RVAP, CONS_RGAS, CONS_PI, CONS_CP, CONS_ALHS)
Definition: cloud.F90:1584
cloud::dqsats_bac
subroutine, public dqsats_bac(DQSi, QSSi, TEMP, PLO, ESTBLX, CONS_H2OMW, CONS_AIRMW)
Definition: cloud.F90:2529
cloud::autoconversion_ls
subroutine, public autoconversion_ls(DT, QC, QP, TE, PL, F, SUNDQV2, SUNDQV3, SUNDQT1, C_00, LWCRIT, DZET)
Definition: cloud.F90:1660
max
#define max(a, b)
Definition: mosaic_util.h:33
cloud::marshpalm
subroutine, public marshpalm(RAIN, PR, DIAM3, NTOTAL, W, VE)
Definition: cloud.F90:2003
cloud::get_ice_fraction
subroutine, public get_ice_fraction(TEMP, T_ICE_ALL, T_ICE_MAX, ICEFRPWR, ICEFRCT)
Definition: cloud.F90:1881
min
#define min(a, b)
Definition: mosaic_util.h:32
cloud::ls_cloud
subroutine, public ls_cloud(DT, ALPHA, PDFSHAPE, PL, TE, QV, QCl, QAl, QCi, QAi, CF, AF, CONS_ALHL, CONS_ALHF, CONS_ALHS, CONS_CP, CONS_H2OMW, CONS_AIRMW, T_ICE_ALL, T_ICE_MAX, ICEFRPWR, ESTBLX, cloud_pertmod, dmp)
Definition: cloud.F90:1107