CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
cldmacro.F
1
3
7 module cldmacro
8!=======================================================================
9! Anning Cheng 2/18/2016 replaced GEO condensation scheme
10! with those from 2M microphysics
11 use wv_saturation, only : epsqs, ttrice, hlatv, hlatf, pcf, rgasv
12! & ,vqsatd2_water_single,
13! & vqsatd2_ice_single,vqsatd2_single
14 use funcphys, only : fpvs, fpvsl, fpvsi
15! use GEOS_UtilsMod, only:QSAT=>GEOS_Qsat, DQSAT=>GEOS_DQsat,
16! & QSATLQ=>GEOS_QsatLQU, QSATIC=>GEOS_QsatICE
17#ifdef GEOS5
18 use mapl_constantsmod, only: mapl_tice , mapl_cp , mapl_grav ,
19 & mapl_alhs , mapl_alhl , mapl_alhf , mapl_rgas , mapl_h2omw,
20 & mapl_airmw, mapl_rvap , mapl_pi , mapl_r8 , mapl_r4
21 use mapl_basemod, only: mapl_undef
22#endif
23#ifdef NEMS_GSM
24 use physcons, mapl_tice => con_t0c, mapl_grav => con_g,
25 & mapl_cp => con_cp, mapl_alhl => con_hvap,
26 & mapl_alhf => con_hfus, mapl_pi => con_pi,
27 & mapl_rgas => con_rd, mapl_rvap => con_rv
28#endif
29
30
31 implicit none
32
33
34! save
35 private
36
37 PUBLIC macro_cloud
38 PUBLIC update_cld
39 public meltfrz_inst
40 public fix_up_clouds_2m
41 PUBLIC cloud_ptr_stubs
42
43!! Some parameters set by PHYSPARAMS
44
45! integer :: NSMAX, DISABLE_RAD, ICEFRPWR
46 integer :: nsmax, disable_rad, icefrpwr
47 &, fr_ls_wat, fr_ls_ice, fr_an_wat, fr_an_ice
48
49 real :: cnv_beta
50 real :: anv_beta
51 real :: ls_beta
52 real :: rh00
53 real :: c_00
54 real :: lwcrit
55 real :: c_acc
56 real :: c_ev_r
57 real :: c_ev_s
58 real :: cldvol2frc
59 real :: rhsup_ice
60 real :: shr_evap_fac
61 real :: min_cld_water
62 real :: cld_evp_eff
63 real :: ls_sdqv2
64 real :: ls_sdqv3
65 real :: ls_sdqvt1
66 real :: anv_sdqv2
67 real :: anv_sdqv3
68 real :: anv_sdqvt1
69 real :: anv_to_ls
70 real :: n_warm
71 real :: n_ice
72 real :: n_anvil
73 real :: n_pbl
74 real :: anv_icefall_c
75 real :: ls_icefall_c
76 real :: revap_off_p
77 real :: cnvenvfc
78 real :: wrhodep
79 real :: t_ice_all
80 real :: cnviceparam
81 real :: cnvddrfc
82 real :: anvddrfc
83 real :: lsddrfc
84! integer :: tanhrhcrit
85! real :: minrhcrit
86! real :: maxrhcrit
87! real :: maxrhcritland
88 real :: turnrhcrit
89 real :: turnrhcrit_upper
90 real :: min_ri, max_ri, min_rl, max_rl, ri_anv
91
92
93 real, parameter :: t_ice_max = mapl_tice
94 real, parameter :: rho_w = 1.0e3
95 real, parameter :: min_cld_frac = 1.0e-8
96 real, parameter :: mapl_alhs = mapl_alhl+mapl_alhf
97
98 real, parameter :: alhlbcp = mapl_alhl/mapl_cp
99 &, alhfbcp = mapl_alhf/mapl_cp
100 &, alhsbcp = alhlbcp+alhfbcp
101
102
103! real, parameter :: PI_0 = 4.*atan(1.)
104 real :: omeps, trinv, t_ice_denom
105
106!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
107
108 contains
109
112 subroutine macro_cloud(IRUN, LM, DT, alf_fac, PP_dev, PPE_dev &
113 &, QLWDTR_dev &
114 &, TH_dev, Q_dev &
115 &, QLW_LS_dev, QLW_AN_dev, QIW_LS_dev &
116 &, QIW_AN_dev, ANVFRC_dev, CLDFRC_dev &
117 &, PHYSPARAMS, SCLMFDFR &
118 &, ALPHT_dev &
119 &, CNV_FICE_dev &
120 &, CNV_NDROP_dev, CNV_NICE_dev, SCICE_dev &
121 &, NCPL_dev, NCPI_dev, PFRZ_dev &
122 &, lprnt, ipr, rhc, pdfflag, qc_min )
123! &, KCBL, lprnt, ipr, rhc )
124
125 integer, intent(in ) :: irun, lm, pdfflag
126 real, intent(in ) :: dt, alf_fac, qc_min(2)
127 real, intent(in ), dimension(IRUN, LM) :: pp_dev
128 real, intent(in ), dimension(IRUN,0:LM) :: ppe_dev
129! real, intent(in ), dimension(IRUN ) :: FRLAND_dev
130! real, intent(in ), dimension(IRUN, LM) :: RMFDTR_dev
131 real, intent(in ), dimension(IRUN, LM) :: qlwdtr_dev
132! real, intent(inout), dimension(IRUN, LM) :: QRN_CU_dev
133! real, intent(inout), dimension(IRUN, LM) :: CNV_UPDFRC_dev
134! real, intent(in ), dimension(IRUN, LM) :: U_dev
135! real, intent(in ), dimension(IRUN, LM) :: V_dev
136 real, intent(in ), dimension(IRUN, LM) :: rhc
137 real, intent(inout), dimension(IRUN, LM) :: th_dev
138 real, intent(inout), dimension(IRUN, LM) :: q_dev
139 real, intent(inout), dimension(IRUN, LM) :: qlw_ls_dev
140 real, intent(inout), dimension(IRUN, LM) :: qlw_an_dev
141 real, intent(inout), dimension(IRUN, LM) :: qiw_ls_dev
142 real, intent(inout), dimension(IRUN, LM) :: qiw_an_dev
143 real, intent(inout), dimension(IRUN, LM) :: anvfrc_dev
144 real, intent(inout), dimension(IRUN, LM) :: cldfrc_dev
145! real, intent( out), dimension(IRUN ) :: PRECU_dev
146! real, intent( out), dimension(IRUN ) :: CUARF_dev
147! real, intent( out), dimension(IRUN ) :: SNRCU_dev
148 real, intent(in ), dimension(58 ) :: physparams
149 real, intent(in ) :: sclmfdfr
150! real, intent(in ), dimension(IRUN, LM) :: QST3_dev
151! real, intent(in ), dimension(IRUN, LM) :: DZET_dev
152! real, intent(in ), dimension(IRUN, LM) :: QDDF3_dev
153! real, intent( out), dimension(IRUN, LM) :: RHX_dev
154! real, intent( out), dimension(IRUN, LM) :: REV_CN_dev
155! real, intent( out), dimension(IRUN, LM) :: RSU_CN_dev
156! real, intent( out), dimension(IRUN, LM) :: ACLL_CN_dev
157! real, intent( out), dimension(IRUN, LM) :: ACIL_CN_dev
158! real, intent( out), dimension(IRUN,0:LM) :: PFL_CN_dev
159! real, intent( out), dimension(IRUN,0:LM) :: PFI_CN_dev
160! real, intent( out), dimension(IRUN, LM) :: PDFL_dev
161! real, intent( out), dimension(IRUN, LM) :: PDFI_dev
162 real, intent( out), dimension(IRUN, LM) :: alpht_dev
163! real, intent( out), dimension(IRUN, LM) :: CFPDF_dev
164! real, intent( out), dimension(IRUN, LM) :: DQRL_dev
165! real, intent( out), dimension(IRUN, LM) :: VFALLSN_CN_dev
166! real, intent( out), dimension(IRUN, LM) :: VFALLRN_CN_dev
167 real, intent(inout), dimension(IRUN, LM) :: cnv_fice_dev
168 real, intent(inout), dimension(IRUN, LM) :: cnv_ndrop_dev
169 real, intent(inout), dimension(IRUN, LM) :: cnv_nice_dev
170 real, intent(inout), dimension(IRUN, LM) :: scice_dev
171 real, intent(inout), dimension(IRUN, LM) :: ncpl_dev
172 real, intent(inout), dimension(IRUN, LM) :: ncpi_dev
173 real, intent(out), dimension(IRUN, LM) :: pfrz_dev
174! real, intent(out), dimension(IRUN, LM) :: QRAIN_CN
175! real, intent(out), dimension(IRUN, LM) :: QSNOW_CN
176
177! real, dimension(IRUN, LM) :: FRZ_PP_dev
178! integer, intent(in), dimension(IRUN) :: KCBL
179 logical lprnt
180 integer ipr
181
182
183! GPU The GPUs need to know how big local arrays are during compile-time
184! as the GPUs cannot allocate memory themselves. This command resets
185! this a priori size to LM for the CPU.
186
187
188 integer :: i , j , k , l
189
190! integer :: FRACTION_REMOVAL
191
192 real :: mass, imass, totfrc, temp, alpha
193 &, dti, tx1, tend, fqi
194
195! real :: MASS, iMASS, TOTFRC, QRN_CU_1D, QSN_CU, QRN_ALL, QSN_ALL
196! &, QTMP1, QTMP2, QTMP3, QTOT, TEMP, RHCRIT, AA3, BB3, ALPHA
197! &, VFALL, VFALLRN, VFALLSN, TOT_PREC_UPD, AREA_UPD_PRC
198! &, AREA_UPD_PRC_tolayer
199! &, PRN_CU_above, PSN_CU_above
200! &, AREA_UPD_PRC_tolayer, U_above,U_below, V_above,V_below
201! &, DZET_above,DZET_below, PRN_CU_above, PSN_CU_above
202! &, EVAP_DD_CU_above, SUBL_DD_CU_above
203! &, NIX, TOTAL_WATER, dti, tx1, tend, fqi, psinv, pops
204
205 logical :: use_autoconv_timescale
206!
207 real, parameter :: rl_cub = 1.0e-15, ri_cub = 6.4e-14
208!
209
210 omeps = 1.0 - epsqs
211 trinv = 1.0 / ttrice
212 dti = 1.0 / dt
213
214 cnv_beta = physparams(1)
215 anv_beta = physparams(2)
216 ls_beta = physparams(3)
217 rh00 = physparams(4)
218 c_00 = physparams(5)
219 lwcrit = physparams(6)
220 c_acc = physparams(7)
221 c_ev_r = physparams(8)
222 c_ev_s = physparams(56)
223 cldvol2frc = physparams(9)
224 rhsup_ice = physparams(10)
225 shr_evap_fac = physparams(11)
226 min_cld_water = physparams(12)
227 cld_evp_eff = physparams(13)
228 nsmax = int( physparams(14) )
229 ls_sdqv2 = physparams(15)
230 ls_sdqv3 = physparams(16)
231 ls_sdqvt1 = physparams(17)
232 anv_sdqv2 = physparams(18)
233 anv_sdqv3 = physparams(19)
234 anv_sdqvt1 = physparams(20)
235 anv_to_ls = physparams(21)
236 n_warm = physparams(22)
237 n_ice = physparams(23)
238 n_anvil = physparams(24)
239 n_pbl = physparams(25)
240 disable_rad = int( physparams(26) )
241 anv_icefall_c = physparams(28)
242 ls_icefall_c = physparams(29)
243 revap_off_p = physparams(30)
244 cnvenvfc = physparams(31)
245 wrhodep = physparams(32)
246 t_ice_all = physparams(33) + mapl_tice
247 cnviceparam = physparams(34)
248 icefrpwr = int( physparams(35) + .001 )
249 cnvddrfc = physparams(36)
250 anvddrfc = physparams(37)
251 lsddrfc = physparams(38)
252! tanhrhcrit = INT( PHYSPARAMS(41) )
253! minrhcrit = PHYSPARAMS(42)
254! maxrhcrit = PHYSPARAMS(43)
255 turnrhcrit = physparams(45) * 0.001
256! maxrhcritland = PHYSPARAMS(46)
257 fr_ls_wat = int( physparams(47) )
258 fr_ls_ice = int( physparams(48) )
259 fr_an_wat = int( physparams(49) )
260 fr_an_ice = int( physparams(50) )
261 min_rl = physparams(51)
262 min_ri = physparams(52)
263 max_rl = physparams(53)
264 max_ri = physparams(54)
265 ri_anv = physparams(55)
266! pdfflag = INT(PHYSPARAMS(57))
267
268
269 turnrhcrit_upper = physparams(58) * 0.001
270
271 use_autoconv_timescale = .false.
272
273 t_ice_denom = 1.0 / (t_ice_max-t_ice_all)
274
275 run_loop: DO i = 1, irun
276! Anning initialization here
277! PRN_CU_above = 0.
278! PSN_CU_above = 0.
279! EVAP_DD_CU_above = 0.
280! SUBL_DD_CU_above = 0.
281! psinv = 1.0 / ppe_dev(i,lm)
282
283 k_loop: DO k = 1, lm
284
285! if (K == 1) then
286! TOT_PREC_UPD = 0.
287! AREA_UPD_PRC = 0.
288! end if
289
290! if (K == LM ) then
291! PRECU_dev(I) = 0.
292! SNRCU_dev(I) = 0.
293! CUARF_dev(I) = 0.
294! end if
295
296! QRN_CU_1D = 0.
297! QSN_CU = 0.
298! VFALL = 0.
299
300! PFL_CN_dev(I,K) = 0.
301! PFI_CN_dev(I,K) = 0.
302
303! IF (K == 1) THEN
304! PFL_CN_dev(I,0) = 0.
305! PFI_CN_dev(I,0) = 0.
306! END IF
307
308! RHX_dev(I,K) = 0.0
309! REV_CN_dev(I,K) = 0.0
310! RSU_CN_dev(I,K) = 0.0
311! ACLL_CN_dev(I,K) = 0.0
312! ACIL_CN_dev(I,K) = 0.0
313! PDFL_dev(I,K) = 0.0
314! PDFI_dev(I,K) = 0.0
315! ALPHT_dev(I,K) = 0.0
316! CFPDF_dev(I,K) = 0.0
317! DQRL_dev(I,K) = 0.0
318! VFALLSN_CN_dev(I,K) = 0.0
319! VFALLRN_CN_dev(I,K) = 0.0
320! VFALLSN = 0.0
321! VFALLRN = 0.0
322
323! DNDCNV_dev(I, K) = 0.0
324! DNCCNV_dev(I, K) = 0.0
325! RAS_DT_dev(I, K) = 0.0
326
327! QRAIN_CN(I,K) = 0.0
328! QSNOW_CN(I,K) = 0.0
329! NIX = 0.0
330
331! QRN_CU_1D = QRN_CU_dev(I,K)
332
333! MASS = (PPE_dev(I,K) - PPE_dev(I,K-1))
334! & * (100./MAPL_GRAV)
335! iMASS = 1.0 / MASS
336! TEMP = TH_dev(I,K)
337! FRZ_PP_dev(I,K) = 0.00
338
339 alpht_dev(i,k) = 0.0
340 mass = (ppe_dev(i,k) - ppe_dev(i,k-1)) * (100./mapl_grav)
341 imass = 1.0 / mass
342 temp = th_dev(i,k)
343
344
345! NOT USED??? - Moorthi
346! TOTAL_WATER = (QIW_AN_dev(I,K) + QLW_AN_dev(I,K)
347! & + QIW_LS_dev(I,K) + QLW_LS_dev(I,K))*MASS
348! & + QLWDTR_dev(I,K)*DT
349
350
352
353 if (temp < t_ice_all) then
354 fqi = 1.0
355 elseif (temp > t_ice_max) then
356 fqi = 0.0
357 else
358 fqi = cnv_fice_dev(i,k)
359 end if
360 tx1 = (1.0-fqi)*qlwdtr_dev(i,k)
361 if (tx1 > 0.0 .and. cnv_ndrop_dev(i,k) <= 0.0) then
362 cnv_ndrop_dev(i,k) = tx1 / ( 1.333 * mapl_pi *rl_cub*997.0)
363 end if
364
365 tx1 = fqi*qlwdtr_dev(i,k)
366 if (tx1 > 0.0 .and. cnv_nice_dev(i,k) <= 0.0) then
367 cnv_nice_dev(i,k) = tx1 / ( 1.333 * mapl_pi *ri_cub*500.0)
368 end if
369
370 tx1 = imass*dt
371 ncpl_dev(i,k) = max(ncpl_dev(i,k)+cnv_ndrop_dev(i,k)*tx1,0.0)
372 ncpi_dev(i,k) = max(ncpi_dev(i,k)+cnv_nice_dev(i,k)*tx1,0.0)
373
374
375! TEND = RMFDTR_dev(I,K)*iMASS * SCLMFDFR
376! ANVFRC_dev(I,K) = min(ANVFRC_dev(I,K) + TEND*DT, 1.0)
377
378!
379! DCNVi_dev(I,K) = (QIW_AN_dev(I,K) - DCNVi_dev(I,K) ) * DTi
380! DCNVL_dev(I,K) = (QLW_AN_dev(I,K) - DCNVL_dev(I,K) ) * DTi
381! DNDCNV_dev(I,K) = (NCPL_dev(I,K) - DNDCNV_dev(I,K)) * DTi
382! DNCCNV_dev(I,K) = (NCPI_dev(I,K) - DNCCNV_dev(I,K)) * DTi
383
384
385! if (k == 1 .or. k == lm) then
386! U_above = 0.0
387! U_below = 0.0
388! V_above = 0.0
389! V_below = 0.0
390! DZET_above = 0.0
391! DZET_below = 0.0
392! else
393! U_above = U_dev(i,k-1)
394! U_below = U_dev(i,k+1)
395! V_above = V_dev(i,k-1)
396! V_below = V_dev(i,k+1)
397! DZET_above = DZET_dev(i,k-1)
398! DZET_below = DZET_dev(i,k+1)
399! end if
400
401! call pdf_spread (K, LM, U_dev(I,K), U_above, U_below,
402! & V_dev(I,K), V_above, V_below,
403! & DZET_above, DZET_below, CNV_UPDFRC_dev(I,K),
404! & PP_dev(I,K), ALPHA, ALPHT_dev(I,K),
405! & FRLAND_dev(I) )
406! pops = PP_dev(I,K) * psinv
407
408! call pdf_spread (K, LM, PP_dev(I,K), ALPHA, ALPHT_dev(I,K),
409! call pdf_spread (K, LM, pops, ALPHA, ALPHT_dev(I,K),
410! & FRLAND_dev(I), rhc(i) )
411
412! ALPHA = max(1.0e-4, 1.0-rhc(i,k))
413! ALPHT_dev(I,K) = ALPHA * alf_fac
414!
415 alpha = max(1.0e-4, 1.0-rhc(i,k)) * alf_fac
416 alpht_dev(i,k) = alpha
417
418! RHCRIT = 1.0 - ALPHA
419
420!================================
421
423 call pfreezing (alpha , pp_dev(i,k) , temp , q_dev(i,k),
424 & qlw_ls_dev(i,k), qlw_an_dev(i,k),
425 & qiw_ls_dev(i,k), qiw_an_dev(i,k),
426 & scice_dev(i,k) , cldfrc_dev(i,k),
427 & anvfrc_dev(i,k), pfrz_dev(i,k), pdfflag)
428
429
430!=============Collect convective precip==============
431
432!*********************** begin of if(false)********************************
433! if(.false.) then
434! QTMP1 = 0.
435! QTMP2 = 0.
436! QTMP3 = 0.
437! QRN_ALL = 0.
438! QSN_ALL = 0.
439
440! if ( TEMP < MAPL_TICE ) then
441!! QTMP2 = QRN_CU_1D
442! QSN_CU = QRN_CU_1D
443! QRN_CU_1D = 0.
444! TEMP = TEMP + QSN_CU * ALHFbCP
445! end if
446
447! AREA_UPD_PRC_tolayer = 0.0
448
449
450! TOT_PREC_UPD = TOT_PREC_UPD + ((QRN_CU_1D + QSN_CU) * MASS)
451! AREA_UPD_PRC = AREA_UPD_PRC + (CNV_UPDFRC_dev(I,K)*
452! & (QRN_CU_1D + QSN_CU )* MASS)
453
454! if ( TOT_PREC_UPD > 0.0 ) AREA_UPD_PRC_tolayer =
455! & MAX( AREA_UPD_PRC/TOT_PREC_UPD, 1.E-6 )
456
457! AREA_UPD_PRC_tolayer = CNV_BETA * AREA_UPD_PRC_tolayer
458
459! IF (K == LM) THEN
460! if (TOT_PREC_UPD > 0.0) AREA_UPD_PRC = MAX( AREA_UPD_PRC/
461! & TOT_PREC_UPD, 1.E-6 )
462! AREA_UPD_PRC = CNV_BETA * AREA_UPD_PRC
463! CUARF_dev(I) = MIN( AREA_UPD_PRC, 1.0 )
464! END IF
465
466
467! CALL MICRO_AA_BB_3 (TEMP,PP_dev(I,K),QST3_dev(I,K),AA3,BB3)
468
469
470! QTMP1 = QLW_LS_dev(I,K) + QLW_AN_dev(I,K)
471! QTMP2 = QIW_LS_dev(I,K) + QIW_AN_dev(I,K)
472! QTOT = QTMP1 + QTMP2
473
474! call PRECIP3 (K, LM, DT, FRLAND_dev(I), RHCRIT, QRN_CU_1D,
475! & QSN_CU, QTMP1, QTMP2, TEMP, Q_dev(I,K), mass,
476! & imass, PP_dev(I,K), DZET_dev(I,K),
477! & QDDF3_dev(I,K), AA3,BB3,AREA_UPD_PRC_tolayer,
478! & PRECU_dev(I), SNRCU_dev(I), PRN_CU_above,
479! & PSN_CU_above, EVAP_DD_CU_above,
480! & SUBL_DD_CU_above, REV_CN_dev(I,K),
481! & RSU_CN_dev(I,K), ACLL_CN_dev(I,K),
482! & ACIL_CN_dev(I,K), PFL_CN_dev(I,K),
483! & PFI_CN_dev(I,K), VFALLRN, VFALLSN,
484! & FRZ_PP_dev(I,K), CNVENVFC, CNVDDRFC,
485! & ANVFRC_dev(I,k), CLDFRC_dev(I,k),
486! & PP_dev(I,KCBL(I)),i)
487
488! VFALLSN_CN_dev(I,K) = VFALLSN
489! VFALLRN_CN_dev(I,K) = VFALLRN
490
491! if (.not. use_autoconv_timescale) then
492! if (VFALLSN .NE. 0.) then
493! QSN_ALL = QSN_ALL + PFI_CN_dev(I,K)/VFALLSN
494! end if
495! if (VFALLRN .NE. 0.) then
496! QRN_ALL = QRN_ALL + PFL_CN_dev(I,K)/VFALLRN
497! end if
498! end if
499
500!! if (.true.) then
501
502! tx1 = QLW_LS_dev(I,K) + QLW_AN_dev(I,K)
503! IF (tx1 > 1.e-20 ) THEN
504! QTMP3 = 1.0 / tx1
505! ELSE
506! QTMP3 = 0.0
507! END IF
508! tx1 = QTMP1 * QTMP3
509! QLW_LS_dev(I,K) = QLW_LS_dev(I,K) * tx1
510! QLW_AN_dev(I,K) = QLW_AN_dev(I,K) * tx1
511! NCPL_dev(I, K) = NCPL_dev(I,K) * tx1
512
513! tx1 = QIW_LS_dev(I,K) + QIW_AN_dev(I,K)
514! IF (tx1 > 1.0e-20 ) THEN
515! QTMP3 = 1.0 / tx1
516! ELSE
517! QTMP3 = 0.0
518! END IF
519! tx1 = QTMP2 * QTMP3
520! QIW_LS_dev(I,K) = QIW_LS_dev(I,K) * tx1
521! QIW_AN_dev(I,K) = QIW_AN_dev(I,K) * tx1
522! NCPI_dev(I, K) = NCPI_dev(I,K) * tx1
523
524
525! QTMP3 = QIW_LS_dev(I,K) + QIW_AN_dev(I,K)
526! & + QLW_LS_dev(I,K) + QLW_AN_dev(I,K)
527
528! If (QTOT > 0.0) then
529! tx1 = QTMP3/QTOT
530! CLDFRC_dev(I,k) = CLDFRC_dev(I,k)*tx1
531! ANVFRC_dev(I,k) = ANVFRC_dev(I,k)*tx1
532! end if
533
534!! end if
535
536
537! tx1 = (MAPL_RGAS*0.01) * temp / PP_dev(I,K)
538
539! QRAIN_CN(I,K) = QRN_ALL * tx1
540! QSNOW_CN(I,K) = QSN_ALL * tx1
541! QRN_CU_dev(I,K) = QRN_CU_1D
542
543! TOTFRC = CLDFRC_dev(I,K) + ANVFRC_dev(I,K)
544
545! IF ( TOTFRC > 1.00 ) THEN
546! tx1 = 1.0 / TOTFRC
547! CLDFRC_dev(I,k) = CLDFRC_dev(I,k) * tx1
548! ANVFRC_dev(I,k) = ANVFRC_dev(I,k) * tx1
549! END IF
550
551! TOTFRC = CLDFRC_dev(I,K) + ANVFRC_dev(I,K)
552
553! end if
554!*********************** end of if(false)********************************
555! if (lprnt .and. i== ipr) write(0,*)'in macrocld1 clffrc=',
556! & CLDFRC_dev(I,K) ,' k=',k
557
559 CALL fix_up_clouds_2m( q_dev(i,k) , temp , qlw_ls_dev(i,k),
560 & qiw_ls_dev(i,k), cldfrc_dev(i,k), qlw_an_dev(i,k),
561 & qiw_an_dev(i,k), anvfrc_dev(i,k), ncpl_dev(i, k),
562 & ncpi_dev(i, k), qc_min)
563
564! if (lprnt .and. i== ipr) write(0,*)'in macrocld1 clffrc=',
565! & CLDFRC_dev(I,K) , ' k=',k
566
567 th_dev(i,k) = temp
568
569 end do k_loop
570
571
572 end do run_loop
573
574 END SUBROUTINE macro_cloud
575
576
577!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
578!! P R O C E S S S U B R O U T I N E S !!
579!! * * * * * !!
580!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
581
582
584 subroutine pdf_spread (K, LM, PP, ALPHA, ALPHT_DIAG, FRLAND, rhc)
585
586 integer, intent(in) :: k,lm
587 real, intent(in) :: PP, FRLAND, rhc
588
589 real, intent(out) :: ALPHA, ALPHT_DIAG
590
591! real, parameter :: slope = 20.0, slope_up = 20.0
592 real, parameter :: slope = 0.02, slope_up = 0.02
593
594 real :: aux1, aux2, maxalpha
595
596! maxalpha = 1.0 - minrhcrit
597 maxalpha = 1.0 - rhc
598
599 aux1 = min(max((pp - turnrhcrit)/slope, -20.0), 20.0)
600 aux2 = min(max((turnrhcrit_upper - pp)/slope_up, -20.0), 20.0)
601
602 if (frland > 0.05) then
603! aux1 = 1.0
604 aux1 = 1.0 / (1.0+exp(aux1+aux1))
605 else
606 aux1 = 2.0 / (1.0+exp(aux1+aux1))
607 end if
608
609 aux2 = 1.0 / (1.0+exp(aux2))
610
611 alpha = max(1.0e-4, min(0.3, maxalpha*aux1*aux2))
612! alpha = min(0.3, maxalpha*aux1*aux2) !Anning
613
614 alpht_diag = alpha
615
616 end subroutine pdf_spread
617
618!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
619
621 subroutine fix_up_clouds_2m(QV, TE, QLC, QIC, CF, QLA, QIA, AF, &
622 & NL, NI, qc_min)
623
624 real, intent(in) :: qc_min(2)
625 real, intent(inout) :: te,qv,qlc,cf,qla,af,qic,qia, nl, ni
626
627! real, parameter :: qmin = 1.0e-8, qmini = 1.0e-7
628! real, parameter :: nmin = 1.0e-3, cfmin = 1.0e-5
629 real, parameter :: nmin = 1.0, cfmin = 1.0e-5
630 &, ri_cub = 6.4e-14, rl_cub = 1.0e-15
631 &, fourb3 = 4.0/3.0
632
633 if (af <= cfmin) then ! Fix if Anvil cloud fraction too small
634 qv = qv + qla + qia
635 te = te - alhlbcp*qla - alhsbcp*qia
636 af = 0.
637 qla = 0.
638 qia = 0.
639
640 if ( cf <= cfmin) then ! Fix if LS cloud fraction too small
641 qv = qv + qlc + qic
642 te = te - alhlbcp*qlc - alhsbcp*qic
643 cf = 0.
644 qlc = 0.
645 qic = 0.
646 endif
647 endif
648
649 if (qlc <= qc_min(1)) then ! LS LIQUID too small
650 qv = qv + qlc
651 te = te - alhlbcp*qlc
652 qlc = 0.
653 endif
654
655 if (qic <= qc_min(2)) then ! LS ICE too small
656 qv = qv + qic
657 te = te - alhsbcp*qic
658 qic = 0.
659 endif
660
661 if (qla <= qc_min(1)) then ! Anvil LIQUID too small
662 qv = qv + qla
663 te = te - alhlbcp*qla
664 qla = 0.
665 endif
666
667 if (qia <= qc_min(2)) then ! Anvil ICE too small
668 qv = qv + qia
669 te = te - alhsbcp*qia
670 qia = 0.
671 endif
672
673 if (qla+qia <= qc_min(1)) then ! Fix ALL cloud quants if Anvil cloud LIQUID+ICE too small
674 qv = qv + qla + qia
675 te = te - alhlbcp*qla - alhsbcp*qia
676 af = 0.
677 qla = 0.
678 qia = 0.
679 endif
680
681 if (qlc+qic <= qc_min(1)) then ! Ditto if LS cloud LIQUID+ICE too small
682 qv = qv + qlc + qic
683 te = te - alhlbcp*qlc - alhsbcp*qic
684 cf = 0.
685 qlc = 0.
686 qic = 0.
687 endif
688
689 if (qla+qlc <= qc_min(1)) then
690 nl = 0.0
691 elseif (nl <= nmin) then ! make sure NL > 0 if Q >0
692 nl = max((qla+qlc)/( fourb3 * mapl_pi *rl_cub*997.0), nmin)
693 endif
694
695 if (qia+qic <= qc_min(2)) then
696 ni = 0.0
697 elseif (ni <= nmin) then ! make sure NI > 0 if Q >0
698 ni = max((qia+qic)/( fourb3 * mapl_pi *ri_cub*500.0), nmin)
699 endif
700
701 end subroutine fix_up_clouds_2m
702
703
704
705!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
709 subroutine update_cld( irun, lm, DT, ALPHA, qc_min, &
710 & PDFSHAPE, PL, QV, QCl, QAl, &
711 & QCi, QAi, TE, CF, AF, &
712 & SCICE, NI, NL)
713! & SCICE, NI, NL, NCnuc)
714
715 integer, intent(in) :: irun, lm, pdfshape
716 real, intent(in) :: dt, qc_min(2)
717 real, intent(in), dimension(irun,lm) :: alpha, pl
718! real, intent(in), dimension(irun,lm) :: ALPHA, PL, NCnuc
719 real, intent(inout), dimension(irun,lm) :: te, qv, qcl, qci &
720 &, CF, QAl, QAi, AF, NI, NL, SCICE
721
722! real :: CFO, pl100, QT, DQ, QSx, DQsx, QCx, QC, QA
723 real :: cfo, pl100, qt, dq, qsx, qcx, qc, qa &
724 &, QX, QSLIQ, QSICE, CFALL, DQx, FQA, tem
725
726 real :: esl, esi, esn !temp use only Anning
727
728 integer :: i,k
729
730 do k=1,lm
731 do i=1,irun
732 if (qv(i,k) > 1.0e-6) then
733 qc = qcl(i,k) + qci(i,k)
734 qa = qal(i,k) + qai(i,k)
735 !Anning do not let empty cloud exist
736 if(qc <= 0.) cf(i,k) = 0.
737 if(qa <= 0.) af(i,k) = 0.
738 qcx = qc + qa
739 qt = qcx + qv(i,k)
740 cfall = af(i,k) + cf(i,k)
741
742!================================================
743! Find the cloud fraction that would correspond to the current condensate
744
745 pl100 = pl(i,k)*100
746
747 if (qcx > 0.0) then
748 tem = 1.0 / qcx
749 fqa = qa *tem
750 esl = min(fpvsl(te(i,k)), pl100)
751 qsliq = min(epsqs*esl/(pl100-omeps*esl), 1.)
752 esi = min(fpvsi(te(i,k)), pl100)
753 qsice = min(epsqs*esi/(pl100-omeps*esi), 1.)
754
755 qsx = ( (qcl(i,k)+qal(i,k))*qsliq
756 & + (qci(i,k)+qai(i,k))*qsice ) *tem
757 else
758 fqa = 0.0
759 esn = min(fpvs(te(i,k)), pl100)
760 qsx = min(epsqs*esn/(pl100-omeps*esn), 1.)
761 endif
762
763! if (TE(i,k) > T_ICE_ALL) SCICE(i,k) = 1.0
764
765 qx = qt - qsx*scice(i,k)
766
767! recalculate QX if too low and SCICE<SHOM
768
769 if (qcx > qc_min(1)) then
770 if (qx <= qcx) then
771 cfo = 1.0 / (1.0 + sqrt(1.0-qx/qcx) )
772! DQ = (Qcx+QCx) / (CFo*CFo)
773 else
774 cfo = 1.0 !Outside of distribution but still with condensate
775! DQ = (QSx+QSx) * ALPHA(i,k)
776 endif
777 else
778 cfo = 0.
779 endif
780
781 cfall = min(1.0, max(cfo, 0.0))
782
783 af(i,k) = cfall * fqa
784 cf(i,k) = cfall - af(i,k)
785
786! if (TE(i,k) > T_ICE_ALL) then ! don't do anything else for cirrus
787
788 call hystpdf( dt, alpha(i,k), pdfshape, qc_min, pl(i,k)
789 &, qv(i,k), qcl(i,k), qal(i,k), qci(i,k)
790 &, qai(i,k), te(i,k), cf(i,k), af(i,k)
791 &, scice(i,k), ni(i,k), nl(i,k))
792
793! endif
794 !Anning do not let empty cloud exist
795 if(qcl(i,k)+qci(i,k) <= 0.0) cf(i,k) = 0.0
796 if(qal(i,k)+qai(i,k) <= 0.0) af(i,k) = 0.0
797 else
798 cf(i,k) = 0.0
799 af(i,k) = 0.0
800 endif
801 enddo
802 enddo
803
804
805 end subroutine update_cld
806
807!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
808!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
811 subroutine hystpdf( DT, ALPHA, PDFSHAPE, qc_min, PL, QV, QCl, QAl&
812 &, QCi, QAi, TE, CF, AF, SCICE, NI, NL)
813! &, QCi, QAi, TE, CF, AF, SCICE, NI, NL, i, k)
814
815 real, intent(in) :: DT, ALPHA, PL, qc_min(2)
816 integer, intent(in) :: pdfshape
817 real, intent(inout) :: TE, QV, QCl, QCi, CF, QAl, QAi, AF, &
818 & NI, NL, SCICE
819
820 integer, parameter :: nmax=10
821
822 real :: QCO, QVO, CFO, QAO, TAU
823 real :: QT, QMX, QMN, DQ, QVtop, sigmaqt1, sigmaqt2, qsnx
824
825! real :: TEO, QSx, DQsx, QS, DQs
826 real :: QSx, DQSx, QS, DQs
827 &, tep, qsp, cfp, qvp, qcp
828 &, ten, qsn, cfn, qvn, qcn
829
830! real :: QCx, QVx, CFx, QAx, QC, QA, fQi, fQi_A
831 real :: QCx, QVx, CFx, QAx, QC, QA, fQi
832 &, dqai, dqal, dqci, dqcl
833
834! real :: QX, QSLIQ, QSICE, CFALL, DQx, FQA, pl100, tmpARR
835 real :: QX, QSLIQ, QSICE, DQx, pl100, tmpARR
836 &, alhx, dqcall, esn, desdt, tc, hltalt, tterm
837
838! integer :: N, i, k
839 integer :: N
840
841 qc = qcl + qci
842 qa = qal + qai
843
844! QT = QC + QA + QV
845! CFALL = AF + CF
846! FQA = 0.0
847! fQi = 0.0
848
849! if (QA+QC > 0.0) FQA = QA / (QA+QC)
850! if (QA > 0.0) fQi_A = QAi / QA
851! if (QT > 0.0) fQi = (QCI+QAI) / QT
852
853! TEo = TE
854!
855 if (te <= t_ice_all) then
856 fqi = 1.0
857 elseif (te >= t_ice_max) then
858 fqi = 0.0
859 else
860 fqi = (1.0 - (te-t_ice_all)*t_ice_denom) ** icefrpwr
861 endif
862
863 pl100 = pl*100
864
865 esn = min(fpvs(te), pl100)
866 qsx = min(epsqs*esn/(pl100-omeps*esn), 1.)
867
868 if (qsx < 1.0) then
869 tc = te - mapl_tice
870 if (te < mapl_tice) then
871 hltalt = hlatv + hlatf * min(-tc*trinv,1.0)
872 else
873 hltalt = hlatv - 2369.0*tc
874 end if
875 if (tc >= -ttrice .and. tc < 0.0) then
876 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3)
877 & + tc*(pcf(4) + tc*pcf(5))))
878 else
879 tterm = 0.0
880 endif
881 desdt = hltalt*esn/(rgasv*te*te) + tterm*trinv
882 dqsx = qsx*pl100*desdt/(esn*(pl100-omeps*esn))
883 else
884 dqsx = 0.0
885 endif
886
887 if (af < 1.0) then
888 tmparr = 1.0 / (1.0-af)
889 else
890 tmparr = 0.0
891 endif
892
893 cfx = cf*tmparr
894 qcx = qc*tmparr
895 qvx = (qv - qsx*af) * tmparr
896
897! if ( AF >= 1.0 ) QVx = QSx*1.e-4
898 if (af > 0.0) then
899 qax = qa/af
900 else
901 qax = 0.0
902 endif
903
904 qt = qcx + qvx
905
906! TEp = TEo
907 qsn = qsx
908 ten = te
909 cfn = cfx
910 qvn = qvx
911 qcn = qcx
912 dqs = dqsx
913
914 do n=1,nmax
915
916 qvp = qvn
917 qcp = qcn
918 cfp = cfn
919 tep = ten
920
921 if(pdfshape < 2) then
922 sigmaqt1 = alpha*qsn
923! sigmaqt1 = ALPHA*QSn
924 sigmaqt2 = sigmaqt1
925 elseif(pdfshape == 2) then
926 sigmaqt1 = alpha*qsn
927 sigmaqt2 = sigmaqt1
928 elseif(pdfshape == 4) then
929 sigmaqt1 = max(alpha/sqrt(3.0), 0.001)
930 else
931 write(0,*)' Aborting : invalid pdfshape=',pdfshape
932 stop
933 endif
934
935 qsnx = qsn*scice
936! if (QCI >= 0.0 .and. qsn > qt) qsnx = qsn
937
938 call pdffrac(pdfshape,qt,sigmaqt1,sigmaqt2,qsnx,cfn)
939 call pdfcondensate(pdfshape,qt,sigmaqt1,sigmaqt2,qsnx,qcn,cfn)
940
941 dqcall = qcn - qcp
942 cf = cfn * ( 1.0-af)
943
944! call Bergeron_iter (DT, PL, TEp, QT, QCi, QAi, QCl, QAl,
945! & CF, AF, NL, NI, DQCALL, fQi)
946
947 if (af > 0.) then
948 qao = qax
949 else
950 qao = 0.
951 end if
952
953
954 alhx = (1.0-fqi)*alhlbcp + fqi*alhsbcp
955
956 if(pdfshape == 1) then
957 qcn = qcp + (qcn- qcp)
958 & / (1.0 - (cfn*(alpha-1.0) - qcn/qsn) *dqs*alhx)
959 elseif(pdfshape == 2) then
960 if (n < nmax) qcn = qcp + ( qcn - qcp ) * 0.5
961 endif
962
963 qvn = qvp - (qcn - qcp)
964 ten = tep + alhx * ((qcn-qcp)*(1.0-af) + (qao-qax)*af)
965
966 if (abs(ten-tep) < 0.00001) exit
967
968 if (ten <= t_ice_all) then
969 fqi = 1.0
970 elseif (ten >= t_ice_max) then
971 fqi = 0.0
972 else
973 fqi = (1.0 - (te-t_ice_all)*t_ice_denom) ** icefrpwr
974 endif
975
976 dqs = 0.0
977
978 if (n < nmax) then
979 esn = min(fpvs(ten), pl100)
980 qsn = min(epsqs*esn/(pl100-omeps*esn), 1.0)
981
982 if (qsn < 1.0) then
983 tc = ten - mapl_tice
984 if (ten < mapl_tice) then
985 hltalt = hlatv + hlatf * min(-tc*trinv,1.0)
986 else
987 hltalt = hlatv - 2369.0*tc
988 end if
989 if (tc >= -ttrice .and. tc < 0.0) then
990 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3)
991 & + tc*(pcf(4) + tc*pcf(5))))
992 else
993 tterm = 0.0
994 end if
995 desdt = hltalt*esn/(rgasv*ten*ten) + tterm*trinv
996 dqs = qsn*pl100*desdt / (esn*(pl100-omeps*esn))
997! else
998! DQS = 0.0
999 endif
1000 endif
1001
1002 enddo
1003
1004 cfo = cfn
1005 cf = cfn
1006 qco = qcn
1007! QVo = QVn
1008! TEo = TEn
1009! TE = TEn
1010
1011! Update prognostic variables. Deal with special case of AF=1
1012! Temporary variables QCo, QAo become updated grid means.
1013
1014 if (af < 1.0) then
1015 cf = cfo * (1.0-af)
1016 qco = qco * (1.0-af)
1017 qao = qao * af
1018 else
1019! Special case AF=1, i.e., box filled with anvil.
1020! - Note: no guarantee QV_box > QS_box
1021
1022 cf = 0. ! Remove any other cloud
1023 qao = qa + qc ! Add any LS condensate to anvil type
1024 qco = 0. ! Remove same from LS
1025 qt = qao + qv ! Total water
1026
1027! Now set anvil condensate to any excess of total water
1028! over QSx (saturation value at top)
1029 qao = max(qt-qsx, 0.)
1030 endif
1031
1032 dqcl = 0.0
1033 dqci = 0.0
1034 dqal = 0.0
1035 dqai = 0.0
1036
1037!large scale QCx is not in envi
1038
1039 qcx = qco - qc
1040! Anning Cheng prevented unstable here
1041! if (QCx < -1.e-3) QCx = -1.e-3
1042 if (qcx < 0.0) then
1043 dqcl = max(qcx, -qcl)
1044 dqci = max(qcx-dqcl, -qci)
1045 else
1046 dqci = qcx * fqi
1047 dqcl = qcx - dqci
1048 end if
1049
1050!Anvil QAx is not in anvil
1051 qax = qao - qa
1052! Anning Cheng prevented unstable here
1053! if(QAx < -1.e-3) QAx = -1.e-3
1054
1055 if (qax < 0.0) then
1056 dqal = max(qax, -qal)
1057 dqai = max(qax-dqal, -qai)
1058 else
1059 dqai = qax * fqi
1060 dqal = qax - dqai
1061 end if
1062
1063! if(.false.) then !Anning turn it off causing unstable
1064 if ( af < 1.e-5 ) then
1065 dqai = -qai
1066 dqal = -qal
1067 af = 0.0
1068 end if
1069 if ( cf < 1.e-5 ) then
1070 dqci = -qci
1071 dqcl = -qcl
1072 cf = 0.0
1073 end if
1074! end if
1075
1076 qai = qai + dqai
1077 qal = qal + dqal
1078 qci = qci + dqci
1079 qcl = qcl + dqcl
1080 qv = qv - ( dqai+dqci+dqal+dqcl)
1081
1082
1083 te = te + (alhlbcp * (dqai+dqci+dqal+dqcl)
1084 & + alhfbcp * (dqai+dqci))
1085
1086
1087
1088 if ( qao <= 0. ) then
1089 qv = qv + qai + qal
1090 te = te - alhsbcp*qai - alhlbcp*qal
1091 qai = 0.
1092 qal = 0.
1093 af = 0.
1094 end if
1095
1096 CALL fix_up_clouds_2m(qv, te, qcl, qci, cf, qal, qai, af, nl, ni
1097 &, qc_min)
1098
1099 end subroutine hystpdf
1100
1101!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1104 subroutine pdffrac (flag,qtmean,sigmaqt1,sigmaqt2,qstar,clfrac)
1105 implicit none
1106
1107 integer flag
1108
1109 real :: qtmean, sigmaqt1, sigmaqt2, qstar, clfrac
1110
1111 real :: qtmode, qtmin, qtmax, qtmedian, aux
1112
1113 if(flag == 1) then
1114 aux = qtmean + sigmaqt1 - qstar
1115 if (aux < 0.0) then
1116 clfrac = 0.
1117 else
1118 if(sigmaqt1 > 0.0) then
1119 clfrac = min(0.5*aux/sigmaqt1, 1.0)
1120 else
1121 clfrac = 1.
1122 endif
1123 endif
1124 elseif(flag == 2) then
1125 qtmode = qtmean + (sigmaqt1-sigmaqt2)/3.
1126 qtmin = min(qtmode-sigmaqt1,0.)
1127 qtmax = qtmode + sigmaqt2
1128 if(qtmax < qstar) then
1129 clfrac = 0.
1130 elseif ( (qtmode <= qstar).and.(qstar < qtmax) ) then
1131 clfrac = (qtmax-qstar)*(qtmax-qstar) /
1132 & ((qtmax-qtmin)*(qtmax-qtmode))
1133 elseif ( (qtmin <= qstar).and.(qstar < qtmode) ) then
1134 clfrac = 1. - ((qstar-qtmin)*(qstar-qtmin)
1135 & /( (qtmax-qtmin)*(qtmode-qtmin)))
1136 elseif ( qstar <= qtmin ) then
1137 clfrac = 1.
1138 endif
1139 elseif(flag == 4) then
1140 if (qtmean > 1.0e-20) then
1141 qtmedian = qtmean*exp(-0.5*sigmaqt1*sigmaqt1)
1142 aux = log(qtmedian/qstar)/sqrt(2.0)/sigmaqt1
1143 aux = min(max(aux, -20.0), 20.0)
1144 clfrac = 0.5*(1.0+erf_app(aux))
1145 else
1146 clfrac = 0.0
1147 end if
1148 endif
1149
1150 return
1151 end subroutine pdffrac
1152
1153!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1154
1157 subroutine pdfcondensate (flag,qtmean4,sigmaqt14,sigmaqt24, &
1158 & qstar4,condensate4, clfrac4)
1159 implicit none
1160
1161 integer flag
1162
1163 real qtmean4, sigmaqt14, sigmaqt24, qstar4, condensate4, clfrac4
1164
1165 real *8 :: qtmode, qtmin, qtmax, constA, constB, cloudf
1166 &, term1, term2, term3
1167 &, qtmean, sigmaqt1, sigmaqt2, qstar, condensate
1168 &, qtmedian, aux, clfrac, tx1
1169
1170 qtmean = dble(qtmean4)
1171 sigmaqt1 = dble(sigmaqt14)
1172 sigmaqt2 = dble(sigmaqt24)
1173 qstar = dble(qstar4)
1174 clfrac = dble(clfrac4)
1175
1176 if(flag == 1) then
1177 if(qtmean+sigmaqt1 < qstar) then
1178 condensate = 0.d0
1179 elseif(qstar > qtmean-sigmaqt1) then
1180 if(sigmaqt1 > 0.d0) then
1181 tx1 = min(qtmean+sigmaqt1-qstar, 2.d0*sigmaqt1)
1182 condensate = tx1*tx1 / (4.d0*sigmaqt1)
1183 else
1184 condensate = qtmean - qstar
1185 endif
1186 else
1187 condensate = qtmean - qstar
1188 endif
1189 elseif(flag == 2) then
1190 qtmode = qtmean + (sigmaqt1-sigmaqt2)/3.d0
1191 qtmin = min(qtmode-sigmaqt1,0.d0)
1192 qtmax = qtmode + sigmaqt2
1193 if ( qtmax < qstar ) then
1194 condensate = 0.d0
1195 elseif ( qtmode <= qstar .and. qstar < qtmax ) then
1196 constb = 2.d0 / ( (qtmax - qtmin)*(qtmax-qtmode) )
1197 cloudf = (qtmax-qstar)*(qtmax-qstar) * 0.5d0 * constb
1198 term1 = (qstar*qstar*qstar)/3.d0
1199 term2 = (qtmax*qstar*qstar)/2.d0
1200 term3 = (qtmax*qtmax*qtmax)/6.d0
1201 condensate = constb * (term1-term2+term3) - qstar*cloudf
1202 elseif ( qtmin <= qstar .and. qstar < qtmode ) then
1203 consta = 2.d0 / ((qtmax-qtmin)*(qtmode-qtmin))
1204 cloudf = 1.d0 - (qstar-qtmin)*(qstar-qtmin)*0.5d0*consta
1205 term1 = qstar*qstar*qstar/3.d0
1206 term2 = qtmin*qstar*qstar/2.d0
1207 term3 = qtmin*qtmin*qtmin/6.d0
1208 condensate = qtmean - (consta*(term1-term2+term3))
1209 & - qstar*cloudf
1210 elseif ( qstar <= qtmin ) then
1211 condensate = qtmean - qstar
1212 endif
1213
1214 elseif(flag == 4) then
1215
1216 if (qtmean > 1.0e-20) then
1217 aux = 0.5*sigmaqt1*sigmaqt1
1218 qtmedian = qtmean*exp(-aux)
1219 aux = (aux + log(qtmedian/qstar))/(sqrt(2.0)*sigmaqt1)
1220 aux = min(max(aux, -20.0), 20.0)
1221 aux = 1.0 + dble(erf_app(sngl(aux)))
1222 condensate = (0.5*qtmean*aux/qstar - clfrac) * qstar
1223 condensate= min(condensate, qtmean)
1224 else
1225 condensate = 0.0
1226 end if
1227
1228 endif
1229 condensate4 = real(condensate)
1230
1231 return
1232 end subroutine pdfcondensate
1233
1234!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1237 subroutine cnvsrc( DT, ICEPARAM, SCLMFDFR, MASS, iMASS, PL, &
1238 & TE, QV, DCF, DMF, QLA, QIA, CF, AF, QS, &
1239 & NL, NI, CNVFICE, CNVNDROP, CNVNICE)
1240
1241 real, intent(in) :: DT, ICEPARAM, SCLMFDFR, MASS, iMASS, QS &
1242 &, DMF,PL, DCF, CF
1243 real, intent(inout) :: TE, AF,QV, QLA, QIA, NI, NL &
1244 &, CNVFICE, CNVNDROP, CNVNICE
1245
1246 real :: TEND,QVx,QCA,fQi
1247
1248 integer, parameter :: STRATEGY = 3
1249 real, parameter :: RL_cub = 1.0e-15, ri_cub = 6.4e-14
1250 &, minrhx = 0.001
1251
1252 fqi = 0.0
1253
1254 if (te < t_ice_all) then
1255 fqi = 1.0
1256 elseif (te > t_ice_max) then
1257 fqi = 0.0
1258 else
1259 fqi = cnvfice
1260 end if
1261
1262
1263! TEND = DCF*iMASS
1264
1265! QLA = QLA + (1.0-fQi)* TEND*DT
1266! QIA = QIA + fQi * TEND*DT
1267
1268
1269 if ( ( (1.0-fqi)*dcf > 0.0) .and. (cnvndrop <= 0.0)) then
1270 cnvndrop = (1.0-fqi)*dcf/( 1.333 * mapl_pi *rl_cub*997.0)
1271 end if
1272
1273 if ((fqi*dcf > 0.0) .and. (cnvnice <= 0.0)) then
1274 cnvnice = fqi*dcf/( 1.333 * mapl_pi *ri_cub*500.0)
1275 end if
1276
1277 nl = max(nl + cnvndrop*imass*dt, 0.0)
1278 ni = max(ni + cnvnice*imass*dt, 0.0)
1279
1280! TE = TE + (alhsbcp-alhlbcp) * fQi * TEND * DT
1281
1282! QCA = QLA + QIA
1283
1284 tend = dmf*imass * sclmfdfr
1285 af = min(af + tend*dt, 1.0)
1286
1287 if ( af < 1.0 ) then
1288 qvx = ( qv - qs * af )/(1.-af)
1289 else
1290 qvx = qs
1291 end if
1292
1293! if (STRATEGY == 1) then
1294! if ( (( QVx - minrhx*QS ) < 0.0 ) .and. (AF > 0.) ) then
1295! AF = (QV - minrhx*QS )/( QS*(1.0-minrhx) )
1296! end if
1297! if ( AF < 0. ) then
1298! AF = 0.
1299! QV = QV + QLA + QIA
1300! TE = TE - (alhlbcp*QLA + alhsbcp*QIA)
1301! QLA = 0.
1302! QIA = 0.
1303! end if
1304! else if (STRATEGY == 2) then
1305! if ( (( QVx - minrhx*QS ) < 0.0 ) .and. (AF > 0.) ) then
1306! QV = QV + (1.-AF)*( minrhx*QS - QVx )
1307! QCA = QCA - (1.-AF)*( minrhx*QS - QVx )
1308! TE = TE - (1.-AF)*( minrhx*QS - QVx )* alhlbcp
1309! end if
1310! end if
1311
1312 end subroutine cnvsrc
1313
1314!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1317 subroutine precip3( K,LM , DT , FRLAND , RHCR3 , QPl , QPi , &
1318 & QCl , QCi , TE , QV , mass , imass , PL , dZE , QDDF3 , AA , BB ,&
1319 & AREA , RAIN , SNOW , PFl_above , PFi_above , EVAP_DD_above, &
1320 & SUBL_DD_above, REVAP_DIAG , RSUBL_DIAG , ACRLL_DIAG , &
1321 & ACRIL_DIAG , PFL_DIAG , PFI_DIAG , VFALLRN , VFALLSN , FRZ_DIAG ,&
1322 & ENVFC,DDRFC, AF, CF, PCBL,i )
1323
1324
1325 integer, intent(in) :: K,LM,i
1326
1327 real, intent(in ) :: DT
1328
1329 real, intent(inout) :: QV,QPl,QPi,QCl,QCi,TE
1330
1331 real, intent(in ) :: mass,imass
1332 real, intent(in ) :: PL
1333 real, intent(in ) :: AA,BB
1334 real, intent(in ) :: RHCR3
1335 real, intent(in ) :: dZE
1336 real, intent(in ) :: QDDF3
1337 real, intent( out) :: RAIN,SNOW
1338 real, intent(in ) :: AREA
1339 real, intent(in ) :: FRLAND
1340
1341 real, intent(inout) :: PFl_above, PFi_above
1342 real, intent(inout) :: EVAP_DD_above, SUBL_DD_above
1343
1344 real, intent( out) :: REVAP_DIAG
1345 real, intent( out) :: RSUBL_DIAG
1346 real, intent( out) :: ACRLL_DIAG,ACRIL_DIAG
1347 real, intent( out) :: PFL_DIAG, PFI_DIAG
1348 real, intent(inout) :: FRZ_DIAG
1349 real, intent( out) :: VFALLSN, VFALLRN
1350
1351 real, intent(in ) :: ENVFC,DDRFC
1352
1353 real, intent(in ) :: AF,CF, PCBL
1354
1355
1356 real :: PFi,PFl,QS,dQS,ENVFRAC
1357 real,save :: TKo,QKo,QSTKo,DQSTKo,RH_BOX,T_ED,QPlKo,QPiKo
1358 real :: Ifactor,RAINRAT0,SNOWRAT0
1359 real :: FALLRN,FALLSN,VEsn,VErn,NRAIN,NSNOW,Efactor
1360
1361 real :: TinLAYERrn,DIAMrn,DROPRAD
1362 real :: TinLAYERsn,DIAMsn,FLAKRAD,pl100
1363
1364 real :: EVAP,SUBL,ACCR,MLTFRZ,EVAPx,SUBLx
1365 real :: EVAP_DD,SUBL_DD,DDFRACT
1366 real :: LANDSEAF
1367
1368 real :: tmpARR, CFR, aux
1369
1370 real, parameter :: TRMV_L = 1.0
1371
1372 real :: TAU_FRZ, TAU_MLT, QSICE, DQSI
1373
1374 integer :: NS, NSMX, itr,L
1375
1376 logical, parameter :: taneff = .true.
1377
1378 real, parameter :: B_SUB = 1.00
1379 real esl,esi, esn,desdt,weight,tc,hlatsb,hlatvp,hltalt,tterm,
1380 & gam
1381 logical lflg
1382
1383
1384 pl100=pl*100
1385 if(taneff) then
1386 aux = min(max((pl- pcbl)/10.0, -20.0), 20.0)
1387 aux = 1.0/(1.0+exp(-aux))
1388 envfrac = envfc + (1.0-envfc)*aux
1389 envfrac = min(envfrac,1.)
1390 else
1391 envfrac = envfc
1392 endif
1393
1394 cfr= af+cf
1395 if ( cfr < 0.99) then
1396 tmparr = 1./(1.-cfr)
1397 else
1398 tmparr = 0.0
1399 end if
1400
1401
1402 IF ( area > 0. ) THEN
1403 ifactor = 1./ ( area )
1404 ELSE
1405 ifactor = 1.00
1406 END if
1407
1408 ifactor = max( ifactor, 1.)
1409
1410 pfl_diag = 0.
1411 pfi_diag = 0.
1412 acril_diag = 0.
1413 acrll_diag = 0.
1414 revap_diag = 0.
1415 rsubl_diag = 0.
1416
1417
1418! dQS = DQSAT( TE, PL, QSAT = QS )
1419! call vqsatd2_single( TE, pl*100., esl,QS,DQS)
1420 esn=min(fpvs(te),pl100)
1421 qs= min(epsqs*esn/(pl100-omeps*esn),1.)
1422! TKO is not defined yet here Anning Cheng
1423 tko = te
1424! QSICE = QSATIC( min(TKo, T_ICE_MAX), PL*100.0 , DQ=DQSI )
1425! call vqsatd2_ice_single(min(TKo, T_ICE_MAX),PL*100.0,
1426! & esl,QSICE,DQSI)
1427 esi=min(fpvsi(tko),pl100)
1428 qsice= min(epsqs*esi/(pl100-omeps*esi),1.)
1429 hltalt = hlatv + hlatf
1430 desdt = hltalt*esi/(rgasv*tko*tko)
1431 if (qsice < 1.0) then
1432 gam = hltalt*qsice*pl100*desdt/(mapl_cp*esi
1433 & * (pl100 - omeps*esi))
1434 else
1435 gam = 0.0
1436 endif
1437 dqsi = (mapl_cp/hltalt)*gam
1438
1439
1440 ddfract = ddrfc
1441
1442 IF (k == 1) THEN
1443 pfl = qpl*mass
1444 pfi = qpi*mass
1445
1446 evap_dd = 0.
1447 subl_dd = 0.
1448
1449 vfallrn = 0.0
1450 vfallsn = 0.0
1451 ELSE
1452 qpl = qpl + pfl_above * imass
1453 pfl = 0.00
1454 qpi = qpi + pfi_above * imass
1455 pfi = 0.00
1456
1457
1458 accr = b_sub * c_acc * ( qpl*mass ) *qcl
1459
1460 accr = min( accr , qcl )
1461
1462 qpl = qpl + accr
1463 qcl = qcl - accr
1464
1465 acrll_diag = accr / dt
1466
1467
1468 accr = b_sub * c_acc * ( qpi*mass ) *qcl
1469
1470 accr = min( accr , qcl )
1471
1472 qpi = qpi + accr
1473 qcl = qcl - accr
1474
1475 te = te + alhfbcp*accr
1476
1477 acril_diag = accr / dt
1478
1479 rainrat0 = ifactor*qpl*mass/dt
1480 snowrat0 = ifactor*qpi*mass/dt
1481
1482 call marshpalmq2(rainrat0,pl,diamrn,nrain,fallrn,vern)
1483 call marshpalmq2(snowrat0,pl,diamsn,nsnow,fallsn,vesn)
1484
1485 IF ( frland < 0.1 ) THEN
1486
1487 END IF
1488
1489 vfallrn = fallrn
1490 vfallsn = fallsn
1491
1492 tinlayerrn = dze / ( max(fallrn,0.)+0.01 )
1493 tinlayersn = dze / ( max(fallsn,0.)+0.01 )
1494
1495 tau_frz = 5000.
1496
1497 mltfrz = 0.0
1498 IF ( (te > mapl_tice ) .and. (te <= mapl_tice+5. ) ) THEN
1499 mltfrz = tinlayersn * qpi *( te - mapl_tice ) / tau_frz
1500 mltfrz = min( qpi , mltfrz )
1501 te = te - alhfbcp*mltfrz
1502 qpl = qpl + mltfrz
1503 qpi = qpi - mltfrz
1504 END IF
1505 frz_diag = frz_diag - mltfrz / dt
1506
1507 mltfrz = 0.0
1508 IF ( te > mapl_tice+5. ) THEN
1509 mltfrz = qpi
1510 te = te - alhfbcp*mltfrz
1511 qpl = qpl + mltfrz
1512 qpi = qpi - mltfrz
1513 END IF
1514 frz_diag = frz_diag - mltfrz / dt
1515
1516 mltfrz = 0.0
1517 if ( k >= lm-1 ) THEN
1518 IF ( te > mapl_tice+0. ) THEN
1519 mltfrz = qpi
1520 te = te - alhfbcp*mltfrz
1521 qpl = qpl + mltfrz
1522 qpi = qpi - mltfrz
1523 END IF
1524 endif
1525 frz_diag = frz_diag - mltfrz / dt
1526
1527
1528 mltfrz = 0.0
1529 IF ( te <= mapl_tice ) THEN
1530 te = te + alhfbcp*qpl
1531 qpi = qpl + qpi
1532 mltfrz = qpl
1533 qpl = 0.
1534 END IF
1535 frz_diag = frz_diag + mltfrz / dt
1536
1537
1538 qko = qv
1539 tko = te
1540 qplko = qpl
1541 qpiko = qpi
1542
1543 do itr = 1,3
1544
1545! DQSTKo = DQSAT ( TKo , PL,QSAT=QSTko )
1546! call vqsatd2_single(TKo, pl*100., esl,QSTko,DQSTKo)
1547 esn=min(fpvs(tko),pl100)
1548 qstko= min(epsqs*esn/(pl100-omeps*esn),1.)
1549 tc = tko - mapl_tice
1550 lflg = (tc >= -ttrice .and. tc < 0.)
1551 weight = min(-tc*trinv,1.0)
1552 hlatsb = hlatv + weight*hlatf
1553 hlatvp = hlatv - 2369.0*tc
1554 if (tko < mapl_tice) then
1555 hltalt = hlatsb
1556 else
1557 hltalt = hlatvp
1558 end if
1559 if (lflg) then
1560 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3)
1561 & +tc*(pcf(4) + tc*pcf(5))))
1562 else
1563 tterm = 0.
1564 end if
1565 desdt = hltalt*esn/(rgasv*tko*tko) + tterm*trinv
1566 dqstko=(epsqs + omeps*qstko)/(pl100 - omeps*esn)*desdt
1567
1568! QSICE = QSATIC( min(TKo, T_ICE_MAX), PL*100.0 , DQ=DQSI )
1569! call vqsatd2_ice_single(min(TKo, T_ICE_MAX),PL*100.0,
1570! & esl,QSICE,DQSI)
1571 esi=min(fpvsi(tko),pl100)
1572 qsice= min(epsqs*esi/(pl100-omeps*esi),1.)
1573 hltalt = hlatv + hlatf
1574 desdt = hltalt*esi/(rgasv*tko*tko)
1575 if (qsice < 1.0) then
1576 gam = hltalt*qsice*pl100*desdt/(mapl_cp*esi
1577 & * (pl100 - omeps*esi))
1578 else
1579 gam = 0.0
1580 endif
1581 dqsi = (mapl_cp/hltalt)*gam
1582
1583 qstko = max( qstko , 1.0e-7 )
1584 qsice = max( qsice , 1.0e-7 )
1585
1586 if (tmparr > 0.0) then
1587 qko =(qko -qstko*cfr)*tmparr
1588 rh_box = qko/qstko
1589 else
1590 rh_box = qko/qstko
1591 end if
1592
1593 IF ( rh_box < rhcr3 ) THEN
1594 efactor = rho_w * ( aa + bb ) / (rhcr3 - rh_box )
1595 else
1596 efactor = 9.99e9
1597 end if
1598
1599
1600 landseaf = 1.00
1601
1602
1603 if ( ( rh_box < rhcr3 ) .AND. ( diamrn > 0.00 ) .AND.
1604 & ( pl > 100.) .AND. ( pl < revap_off_p ) ) then
1605 droprad = 0.5*diamrn
1606 t_ed = efactor * droprad**2
1607 t_ed = t_ed * ( 1.0 + dqstko*alhlbcp )
1608
1609 evap = qpl*(1.0 - exp( -c_ev_r * vern * landseaf *envfrac*
1610 & tinlayerrn / t_ed ) )
1611 ELSE
1612 evap = 0.0
1613 END if
1614
1615
1616
1617 if (tmparr > 0.0) then
1618 qko = (qko -qsice*cfr)*tmparr
1619 rh_box = qko/qsice
1620 else
1621 rh_box = qko/qsice
1622 end if
1623 IF ( rh_box < rhcr3 ) THEN
1624 efactor = 0.5*rho_w * ( aa+bb ) / (rhcr3-rh_box )
1625 else
1626 efactor = 9.99e9
1627 end if
1628
1629
1630 if ( ( rh_box < rhcr3 ) .AND. ( diamsn > 0.00 ) .AND.
1631 & ( pl > 100. ) .AND. ( pl < revap_off_p ) ) then
1632 flakrad = 0.5*diamsn
1633 t_ed = efactor * flakrad**2
1634 t_ed = t_ed * ( 1.0 + dqsi*alhsbcp)
1635 subl = qpi*(1.0 - exp( -c_ev_s * vesn * landseaf * envfrac
1636 & * tinlayersn / t_ed ) )
1637
1638 ELSE
1639 subl = 0.0
1640 END IF
1641
1642 if (itr == 1) then
1643 evapx = evap
1644 sublx = subl
1645 else
1646 evap = (evap+evapx) /2.0
1647 subl = (subl+sublx) /2.0
1648 endif
1649
1650 evap = evap*(1.-cfr)
1651 subl = subl*(1.-cfr)
1652! Anning prevent negative QPi and QPl
1653 subl = min(qpi, max(subl,0.))
1654 evap = min(qpl, max(evap,0.))
1655
1656
1657 qko =qko + evap + subl
1658 tko = tko - evap * alhlbcp - subl * alhsbcp
1659
1660 enddo
1661 qpi = qpi - subl
1662 qpl = qpl - evap
1663
1664
1665 evap_dd = evap_dd_above + ddfract*evap*mass
1666 evap = evap - ddfract*evap
1667 subl_dd = subl_dd_above + ddfract*subl*mass
1668 subl = subl - ddfract*subl
1669
1670
1671 qv = qv + evap + subl
1672 te = te - evap * alhlbcp - subl * alhsbcp
1673
1674 revap_diag = evap / dt
1675 rsubl_diag = subl / dt
1676
1677 pfl = qpl*mass
1678 pfi = qpi*mass
1679
1680 pfl_diag = pfl/dt
1681 pfi_diag = pfi/dt
1682 end if
1683
1684
1685
1686 evap = qddf3*evap_dd/mass
1687 subl = qddf3*subl_dd/mass
1688! Anning prevent negative QPi and QPl
1689 subl = min(qpi, max(subl,0.))
1690 evap = min(qpl, max(evap,0.))
1691 qv = qv + evap + subl
1692 te = te - evap * alhlbcp - subl * alhsbcp
1693 revap_diag = revap_diag + evap / dt
1694 rsubl_diag = rsubl_diag + subl / dt
1695
1696 IF (k == lm) THEN
1697 rain = pfl/dt
1698 snow = pfi/dt
1699 END IF
1700
1701 qpi = 0.
1702 qpl = 0.
1703
1704 pfl_above = pfl
1705 pfi_above = pfi
1706
1707 evap_dd_above = evap_dd
1708 subl_dd_above = subl_dd
1709
1710 end subroutine precip3
1711
1712!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1715 subroutine marshpalmq2(RAIN,PR,DIAM3,NTOTAL,W,VE)
1716
1717 real, intent(in ) :: RAIN,PR
1718 real, intent(out) :: DIAM3,NTOTAL,W,VE
1719
1720 real :: RAIN_DAY,LAMBDA,A,B,SLOPR,DIAM1
1721
1722 real, parameter :: N0 = 0.08
1723
1724 INTEGER :: IQD
1725
1726 real :: RX(8) , D3X(8)
1727
1728
1729
1730 rx = (/ 0. , 5. , 20. , 80. , 320. , 1280., 4*1280., 16*1280. /)
1731 d3x = (/ 0.019, 0.032, 0.043, 0.057, 0.076, 0.102, 0.137 ,
1732 & 0.183 /)
1733
1734 rain_day = rain * 3600. *24.
1735
1736 IF ( rain_day <= 0.00 ) THEN
1737 diam1 = 0.00
1738 diam3 = 0.00
1739 ntotal= 0.00
1740 w = 0.00
1741 END IF
1742
1743 DO iqd = 1,7
1744 IF ( (rain_day <= rx(iqd+1)) .AND. (rain_day > rx(iqd))) THEN
1745 slopr =( d3x(iqd+1)-d3x(iqd) ) / ( rx(iqd+1)-rx(iqd))
1746 diam3 = d3x(iqd) + (rain_day-rx(iqd))*slopr
1747 END IF
1748 END DO
1749
1750 IF ( rain_day >= rx(8) ) THEN
1751 diam3=d3x(8)
1752 END IF
1753
1754 ntotal = 0.019*diam3
1755
1756 diam3 = 0.664 * diam3
1757
1758 w = (2483.8 * diam3 + 80.)*sqrt(1000./pr)
1759
1760 ve = max( 0.99*w/100. , 1.000 )
1761
1762 diam1 = 3.0*diam3
1763
1764 diam1 = diam1/100.
1765 diam3 = diam3/100.
1766 w = w/100.
1767 ntotal = ntotal*1.0e6
1768
1769 end subroutine marshpalmq2
1770!==========================================================
1773 subroutine micro_aa_bb_3(TEMP,PR,Q_SAT,AA,BB)
1774
1775 real, intent(in ) :: TEMP,Q_SAT
1776 real, intent(in ) :: PR
1777 real, intent(out) :: AA,BB
1778
1779 real :: E_SAT
1780
1781 real, parameter :: EPSILON = 0.622
1782 real, parameter :: K_COND = 2.4e-2
1783 real, parameter :: DIFFU = 2.2e-5
1784
1785 e_sat = 100.* pr * q_sat /( (epsilon) + (1.0-(epsilon))*q_sat )
1786
1787 aa = ( get_alhx3(temp)**2 ) / ( k_cond*mapl_rvap*(temp**2) )
1788
1789
1790 bb = mapl_rvap*temp / ( diffu*(1000./pr)*e_sat )
1791
1792 end subroutine micro_aa_bb_3
1793
1794!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1795
1797 function ldradius3(PL,TE,QCL,NN) RESULT(RADIUS)
1798
1799 real, intent(in) :: te,pl,nn,qcl
1800 real :: radius
1801
1802 real :: muu,rho
1803
1804
1805 rho = 100.*pl / (mapl_rgas*te )
1806 muu = qcl * rho
1807 radius = muu/(nn*rho_w*(4./3.)*mapl_pi)
1808 radius = radius**(1./3.)
1809
1810
1811 end function ldradius3
1812
1813!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1814
1816 function ice_fraction (TEMP) RESULT(ICEFRCT)
1817 real, intent(in) :: temp
1818 real :: icefrct
1819
1820 icefrct = 0.00
1821 if ( temp <= t_ice_all ) then
1822 icefrct = 1.000
1823 else if ( (temp > t_ice_all) .AND. (temp <= t_ice_max) ) then
1824 icefrct = 1.00 - ( temp - t_ice_all ) / ( t_ice_max -
1825 & t_ice_all )
1826 end if
1827 icefrct = min(icefrct,1.00)
1828 icefrct = max(icefrct,0.00)
1829
1830 icefrct = icefrct**icefrpwr
1831
1832 end function ice_fraction
1833
1834!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1836 function get_alhx3(T) RESULT(ALHX3)
1837
1838 real, intent(in) :: t
1839 real :: alhx3
1840
1841 real :: t_x
1842
1843 t_x = t_ice_max
1844
1845 if ( t < t_ice_all ) then
1846 alhx3=mapl_alhs
1847 end if
1848
1849 if ( t > t_x ) then
1850 alhx3=mapl_alhl
1851 end if
1852
1853 if ( (t <= t_x) .and. (t >= t_ice_all) ) then
1854 alhx3 = mapl_alhs + (mapl_alhl-mapl_alhs)
1855 & *( t - t_ice_all )/( t_x - t_ice_all )
1856 end if
1857
1858 end function get_alhx3
1859
1860
1861!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1862
1864 real function icefrac(t,t_trans,t_freez)
1865
1866 real, intent(in) :: t
1867 real, intent(in),optional :: t_trans
1868 real, intent(in),optional :: t_freez
1869
1870 real :: t_x,t_f
1871
1872 if (present( t_trans )) then
1873 t_x = t_trans
1874 else
1875 t_x = t_ice_max
1876 endif
1877 if (present( t_freez )) then
1878 t_f = t_freez
1879 else
1880 t_f = t_ice_all
1881 endif
1882
1883
1884 if ( t < t_f ) icefrac=1.000
1885
1886 if ( t > t_x ) icefrac=0.000
1887
1888 if ( t <= t_x .and. t >= t_f ) then
1889 icefrac = 1.00 - ( t - t_f ) /( t_x - t_f )
1890 endif
1891
1892 end function icefrac
1893
1894
1895
1896!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1897!Parititions DQ into ice and liquid. Follows Morrison and Gettelman, 2008
1898!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1901 subroutine bergeron_iter ( DTIME , PL , TE , QV , QILS , QICN ,
1902 & QLLS , QLCN , CF , AF , NL , NI , DQALL , FQI )
1903
1904 real , intent(in ) :: DTIME, PL, TE
1905 real , intent(inout ) :: DQALL
1906 real , intent(in) :: QV, QLLS, QLCN, QICN, QILS
1907 real , intent(in) :: CF, AF, NL, NI
1908 real, intent (out) :: FQI
1909 real :: DC, TEFF,QCm,DEP, QC, QS, RHCR, DQSL, DQSI, QI, TC, DIFF,
1910 & denair, denice, aux, dcf, qtot, lhcorr, ql, dqi, dql, qvinc,
1911 & qsliq, cfall, new_qi, new_ql, qsice, fqi_0, qs_0, dqs_0, fqa,
1912 & nix
1913 real esl,esi, esn,desdt,weight,hlatsb,hlatvp,hltalt,tterm,
1914 & gam,pl100
1915
1916
1917 pl100=pl*100
1918 diff = 0.0
1919 dep=0.0
1920 qi = qils + qicn
1921 ql = qlls +qlcn
1922 qtot=qi+ql
1923 fqa = 0.0
1924 if (qtot .gt. 0.0) fqa = (qicn+qils)/qtot
1925 nix= (1.0-fqa)*ni
1926
1927 dqall=dqall/dtime
1928 cfall= min(cf+af, 1.0)
1929 tc=te-273.0
1930 fqi_0 = fqi
1931
1932
1933
1934 if (te .ge. t_ice_max) then
1935 fqi = 0.0
1936 elseif(te .le. t_ice_all) then
1937 fqi = 1.0
1938 else
1939
1940
1941 fqi = 0.0
1942 if (qils .le. 0.0) return
1943
1944 qvinc= qv
1945! QSLIQ = QSATLQ( TE , PL*100.0 , DQ=DQSL )
1946! call vqsatd2_water_single(TE,PL*100.0,esl,QSLIQ,DQSL)
1947 esl=min(fpvsl(te),pl100)
1948 qsliq= min(epsqs*esl/(pl100-omeps*esl),1.)
1949
1950! QSICE = QSATIC( TE , PL*100.0 , DQ=DQSI )
1951! call vqsatd2_ice_single(TE,PL*100.0,esl,QSICE,DQSI)
1952 esi=min(fpvsi(te),pl100)
1953 qsice= min(epsqs*esi/(pl100-omeps*esi),1.)
1954 hltalt = hlatv + hlatf
1955 desdt = hltalt*esi/(rgasv*te*te)
1956 if (qsice < 1.0) then
1957 gam = hltalt*qsice*pl100*desdt/(mapl_cp*esi
1958 & * (pl100 - omeps*esi))
1959 else
1960 gam = 0.0
1961 endif
1962 dqsi = (mapl_cp/hltalt)*gam
1963
1964
1965
1966
1967 qvinc =min(qvinc, qsliq)
1968
1969 diff=(0.211*1013.25/(pl+0.1))*(((te+0.1)/273.0)**1.94)*1e-4
1970 denair = pl*100.0/mapl_rgas/te
1971 denice = 1000.0*(0.9167 - 1.75e-4*tc -5.0e-7*tc*tc)
1972 lhcorr = ( 1.0 + dqsi*alhsbcp)
1973
1974 if ((nix .gt. 1.0) .and. (qils .gt. 1.0e-10)) then
1975 dc = max((qils/(nix*denice*mapl_pi))**(0.333), 20.0e-6)
1976 else
1977 dc = 20.0e-6
1978 end if
1979
1980
1981 teff= nix*denair*2.0*mapl_pi*diff*dc/lhcorr
1982
1983 dep=0.0
1984 if ((teff .gt. 0.0) .and. (qils .gt. 1.0e-14)) then
1985 aux =max(min(dtime*teff, 20.0), 0.0)
1986 dep=(qvinc-qsice)*(1.0-exp(-aux))/dtime
1987 end if
1988
1989 dep=max(dep, -qils/dtime)
1990
1991 dqi = 0.0
1992 dql = 0.0
1993 fqi=0.0
1994
1995 if (dqall .ge. 0.0) then
1996
1997 if (dep .gt. 0.0) then
1998 dqi = min(dep, dqall + qlls/dtime)
1999 dql = dqall - dqi
2000 else
2001 dql=dqall
2002 dqi = 0.0
2003 end if
2004 end if
2005
2006
2007 if (dqall .lt. 0.0) then
2008 dql = max(dqall, -qlls/dtime)
2009 dqi = max(dqall - dql, -qils/dtime)
2010 end if
2011
2012 if (dqall .ne. 0.0) fqi=max(min(dqi/dqall, 1.0), 0.0)
2013 end if
2014 end subroutine bergeron_iter
2015
2016
2017
2018!=============================================================================
2023 subroutine pfreezing ( ALPHA , PL , TE , QV , QCl , QAl , QCi , &
2024 & QAi , SC_ICE , CF , AF , PF, pdfflag)
2025
2026 integer, intent(in) :: pdfflag
2027 real , intent(in) :: PL, ALPHA, QV, SC_ICE, AF, TE, &
2028 & QCl, QCi, QAl, QAi, CF
2029 real , intent(out) :: PF
2030
2031 real :: qt, QCx, QSn, tmpARR, CFALL, QVx, CFio, QA, QAx, QC, QI, &
2032 & QL, DQSx, sigmaqt1, sigmaqt2, qsnx, esl, esi,pl100
2033
2034 pl100 = pl*100
2035
2036 qa = qal + qai
2037 qc = qcl + qci
2038 cfall = af
2039
2040 if ( cfall >= 1.0 ) then
2041 pf = 0.0
2042 return
2043 end if
2044
2045! QSn = QSATIC( TE , PL*100.0 , DQ=DQSx )
2046! call vqsatd2_ice_single(TE,PL*100.0,esl,QSn,DQSx)
2047
2048 esi = min(fpvsi(te),pl100)
2049 qsn = max(min(epsqs*esi/(pl100-omeps*esi), 1.0), 1.0e-9)
2050
2051 tmparr = 0.0
2052 if ( cfall < 0.99 ) then
2053 tmparr = 1./(1.0-cfall)
2054 end if
2055
2056 qcx = qc*tmparr
2057 qvx = ( qv - qsn*cfall )*tmparr
2058! QVx = QV*tmpARR
2059
2060 qt = qcx + qvx
2061
2062 cfio = 0.0
2063
2064 qsn = qsn*sc_ice
2065
2066 if(pdfflag < 2) then
2067 sigmaqt1 = max(alpha, 0.1) * qsn
2068 sigmaqt2 = max(alpha, 0.1) * qsn
2069 elseif(pdfflag == 2) then
2070! for triangular, symmetric: sigmaqt1 = sigmaqt2 = alpha*qsn (alpha is half width)
2071! for triangular, skewed r : sigmaqt1 < sigmaqt2
2072! try: skewed right below 500 mb
2073!!! if(pl.lt.500.) then
2074 sigmaqt1 = alpha * qsn
2075 sigmaqt2 = alpha * qsn
2076 elseif(pdfflag == 4) then
2077 sigmaqt1 = max(alpha/sqrt(3.0), 0.001)
2078 else
2079 write(0,*)' Aborting : invalid pdfflag=',pdfflag
2080 stop
2081 endif
2082
2083 call pdffrac(pdfflag,qt,sigmaqt1,sigmaqt2,qsn,cfio)
2084
2085 pf = min(max(cfio*(1.0-cfall), 0.0), 0.999)
2086
2087
2088 end subroutine pfreezing
2089
2090
2091
2092
2093!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Instantaneous freezing of condensate!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2094
2097 subroutine meltfrz_inst(IM, LM, TE, QCL, QAL, QCI, QAI, NL, NI)
2098
2099 integer, intent(in) :: im, lm
2100 real , intent(inout), dimension(:,:) :: te,qcl,qci, qal, qai, &
2101 & NI, NL
2102
2103 real , dimension(im,lm) :: fqi,dqil, dqmax, qltot, qitot, fqal,
2104 & fqai, dnil, fqa
2105
2106 qitot = qci+qai
2107 qltot = qcl + qal
2108 fqa = 0.0
2109
2110 where (qitot+qltot > 0.0)
2111 fqa= (qai+qal)/(qitot+qltot)
2112 end where
2113
2114 dqil = 0.0
2115 dnil = 0.0
2116 dqmax = 0.0
2117
2118 where( te <= t_ice_all )
2119 dqmax = (t_ice_all - te)/(alhsbcp-alhlbcp)
2120 dqil = min(qltot , dqmax)
2121 end where
2122
2123 where ((dqil <= dqmax) .and. (dqil > 0.0))
2124 dnil = nl
2125 end where
2126
2127 where ((dqil > dqmax) .and. (dqil > 0.0))
2128 dnil = nl*dqmax/dqil
2129 end where
2130
2131 dqil = max( 0., dqil )
2132! Anning for moisture conservation 11/22/2016
2133! QITOT = max(QITOT + dQil, 0.0)
2134! QLTOT = max(QLTOT - dQil, 0.0)
2135 dqil = min(qltot,dqil)
2136 qitot = qitot + dqil
2137 qltot = qltot - dqil
2138 nl = nl - dnil
2139 ni = ni + dnil
2140 te = te + (alhsbcp-alhlbcp)*dqil
2141
2142 dqil = 0.0
2143 dnil = 0.0
2144 dqmax = 0.0
2145
2146
2147 where( te > t_ice_max )
2148 dqmax = (te-t_ice_max) / (alhsbcp-alhlbcp)
2149 dqil = min(qitot, dqmax)
2150 endwhere
2151
2152 where ((dqil .le. dqmax) .and. (dqil .gt. 0.0))
2153 dnil = ni
2154 end where
2155 where ((dqil .gt. dqmax) .and. (dqil .gt. 0.0))
2156 dnil = ni*dqmax/dqil
2157 end where
2158 dqil = max( 0., dqil )
2159! Anning for moisture conservation 11/22/2016
2160! QLTOT = max(QLTOT+ dQil, 0.)
2161! QITOT = max(QITOT - dQil, 0.)
2162 dqil = min(qitot,dqil)
2163 qitot = qitot - dqil
2164 qltot = qltot + dqil
2165 nl = nl + dnil
2166 ni = ni - dnil
2167
2168 te = te - (alhsbcp-alhlbcp)*dqil
2169
2170 qci = qitot*(1.0-fqa)
2171 qai = qitot*fqa
2172 qcl = qltot*(1.0-fqa)
2173 qal = qltot*fqa
2174
2175 end subroutine meltfrz_inst
2176
2177
2178
2179!======================================
2181 subroutine cloud_ptr_stubs ( &
2182 & SMAXL, SMAXI, WSUB, CCN01, CCN04, CCN1, NHET_NUC, NLIM_NUC, SO4, &
2183 & ORG, BCARBON, DUST, SEASALT, NCPL_VOL, NCPI_VOL, NRAIN, NSNOW, &
2184 & CDNC_NUC, INC_NUC, SAT_RAT, QSTOT, QRTOT, CLDREFFS, CLDREFFR, &
2185 & DQVDT_micro,DQIDT_micro, DQLDT_micro, DTDT_micro, RL_MASK, &
2186 & RI_MASK, KAPPA, SC_ICE, CFICE, CFLIQ, RHICE, RHLIQ, RAD_CF, &
2187 & RAD_QL, RAD_QI, RAD_QS, RAD_QR, RAD_QV, CLDREFFI, CLDREFFL, &
2188 & NHET_IMM, NHET_DEP, NHET_DHF, DUST_IMM, DUST_DEP, DUST_DHF, SCF, &
2189 & SCF_ALL, SIGW_GW, SIGW_CNV, SIGW_TURB, SIGW_RC, RHCmicro, &
2190 & DNHET_IMM, NONDUST_IMM, NONDUST_DEP, BERG, BERGSO, MELT, &
2191 & DNHET_CT, DTDT_macro, QCRES, DT_RASP, FRZPP_LS, SNOWMELT_LS, &
2192 & QIRES, AUTICE, PFRZ, DNCNUC, DNCHMSPLIT, DNCSUBL, DNCAUTICE, &
2193 & DNCACRIS, DNDCCN, DNDACRLS, DNDEVAPC, DNDACRLR, DNDAUTLIQ)
2194! & DNDCNV, DNCCNV)
2195
2196
2197
2198 real , pointer , dimension(:,:,:) :: smaxl,smaxi, wsub, ccn01, &
2199 & CCN04, CCN1, NHET_NUC, NLIM_NUC, SO4, ORG, BCARBON, DUST, &
2200 & SEASALT, NCPL_VOL, NCPI_VOL, NRAIN, NSNOW, CDNC_NUC, INC_NUC, &
2201 & SAT_RAT, QSTOT, QRTOT, CLDREFFS, CLDREFFR, DQVDT_micro, &
2202 &DQIDT_micro, DQLDT_micro, DTDT_micro, RL_MASK, RI_MASK, KAPPA, &
2203 & SC_ICE, CFICE, CFLIQ, RHICE, RHLIQ, ALPH, RAD_CF, RAD_QL, RAD_QI,&
2204 & RAD_QS, RAD_QR, RAD_QV, CLDREFFI, CLDREFFL, NHET_IMM, NHET_DEP, &
2205 & NHET_DHF, DUST_IMM, DUST_DEP, DUST_DHF, SCF, SCF_ALL, SIGW_GW, &
2206 & SIGW_CNV, SIGW_TURB, SIGW_RC, RHCmicro, DNHET_IMM, NONDUST_IMM, &
2207 & NONDUST_DEP, BERG, BERGSO, MELT, DNHET_CT, DTDT_macro, QCRES, &
2208 & DT_RASP, FRZPP_LS, SNOWMELT_LS, QIRES, AUTICE, PFRZ, DNCNUC, &
2209 & DNCHMSPLIT, DNCSUBL, DNCAUTICE, DNCACRIS, DNDCCN, DNDACRLS, &
2210 & DNDEVAPC, DNDACRLR, DNDAUTLIQ
2211! & DNDEVAPC, DNDACRLR, DNDAUTLIQ, DNDCNV, DNCCNV
2212
2213
2214!DONIF
2215
2216 IF( ASSOCIATED(smaxl) ) smaxl = 0.
2217 IF( ASSOCIATED(smaxi) ) smaxi = 0.
2218 IF( ASSOCIATED(wsub) ) wsub = 0.
2219 IF( ASSOCIATED(ccn01) ) ccn01 = 0.
2220 IF( ASSOCIATED(ccn04) ) ccn04 = 0.
2221 IF( ASSOCIATED(ccn1) ) ccn1 = 0.
2222 IF( ASSOCIATED(nhet_nuc) ) nhet_nuc = 0.
2223 IF( ASSOCIATED(nlim_nuc) ) nlim_nuc = 0.
2224 IF( ASSOCIATED(so4) ) so4 = 0.
2225 IF( ASSOCIATED(org) ) org = 0.
2226 IF( ASSOCIATED(bcarbon) ) bcarbon = 0.
2227 IF( ASSOCIATED(dust) ) dust = 0.
2228 IF( ASSOCIATED(seasalt) ) seasalt = 0.
2229 IF( ASSOCIATED(ncpl_vol) ) ncpl_vol = 0.
2230 IF( ASSOCIATED(ncpi_vol) ) ncpi_vol = 0.
2231
2232 IF( ASSOCIATED(nrain) ) nrain = 0.
2233 IF( ASSOCIATED(nsnow) ) nsnow = 0.
2234 IF( ASSOCIATED(cdnc_nuc) ) cdnc_nuc = 0.
2235 IF( ASSOCIATED(inc_nuc) ) inc_nuc = 0.
2236 IF( ASSOCIATED(sat_rat) ) sat_rat = 0.
2237 IF( ASSOCIATED(qstot) ) qstot = 0.
2238 IF( ASSOCIATED(qrtot) ) qrtot = 0.
2239
2240 IF( ASSOCIATED(dqvdt_micro) ) dqvdt_micro = 0.
2241 IF( ASSOCIATED(dqidt_micro) ) dqidt_micro = 0.
2242 IF( ASSOCIATED(dqldt_micro) ) dqldt_micro = 0.
2243 IF( ASSOCIATED(dtdt_micro) ) dtdt_micro = 0.
2244 IF( ASSOCIATED(dtdt_macro) ) dtdt_macro = 0.
2245
2246 IF( ASSOCIATED(rl_mask) ) rl_mask = 0.
2247 IF( ASSOCIATED(ri_mask) ) ri_mask = 0.
2248 IF( ASSOCIATED(kappa) ) kappa = 0.
2249 IF( ASSOCIATED(sc_ice)) sc_ice = 0.
2250 IF( ASSOCIATED(rhice) ) rhice = 0.
2251 IF( ASSOCIATED(rhliq) ) rhliq = 0.
2252 IF( ASSOCIATED(cfice) ) cfice = 0.
2253 IF( ASSOCIATED(cfliq) ) cfliq = 0.
2254 IF( ASSOCIATED(alph) ) alph = 0.
2255
2256
2257 IF( ASSOCIATED(rad_cf) ) rad_cf = 0.
2258 IF( ASSOCIATED(rad_ql) ) rad_ql = 0.
2259 IF( ASSOCIATED(rad_qi) ) rad_qi = 0.
2260 IF( ASSOCIATED(rad_qs) ) rad_qs = 0.
2261 IF( ASSOCIATED(rad_qr) ) rad_qr = 0.
2262 IF( ASSOCIATED(rad_qv) ) rad_qv = 0.
2263 IF( ASSOCIATED(cldreffi) ) cldreffi = 0.
2264 IF( ASSOCIATED(cldreffl) ) cldreffl = 0.
2265 IF( ASSOCIATED(cldreffs) ) cldreffs = 0.
2266 IF( ASSOCIATED(cldreffr) ) cldreffr = 0.
2267
2268 IF( ASSOCIATED(nhet_imm) ) nhet_imm = 0.
2269 IF( ASSOCIATED(nhet_dep) ) nhet_dep = 0.
2270 IF( ASSOCIATED(nhet_dhf) ) nhet_dhf = 0.
2271 IF( ASSOCIATED(dust_imm) ) dust_imm = 0.
2272 IF( ASSOCIATED(dust_dep) ) dust_dep = 0.
2273 IF( ASSOCIATED(dust_dhf) ) dust_dhf = 0.
2274 IF( ASSOCIATED(nondust_imm) ) nondust_imm = 0.
2275 IF( ASSOCIATED(nondust_dep) ) nondust_dep = 0.
2276
2277
2278 IF( ASSOCIATED(scf) ) scf = 0.
2279 IF( ASSOCIATED(scf_all) ) scf_all = 0.
2280 IF( ASSOCIATED(sigw_gw) ) sigw_gw = 0.
2281 IF( ASSOCIATED(sigw_cnv) ) sigw_cnv = 0.
2282 IF( ASSOCIATED(sigw_turb) ) sigw_turb = 0.
2283 IF( ASSOCIATED(sigw_rc) ) sigw_rc = 0.
2284 IF( ASSOCIATED(rhcmicro) ) rhcmicro = 0.
2285 IF( ASSOCIATED(dnhet_imm) ) dnhet_imm = 0.
2286 IF( ASSOCIATED(berg) ) berg = 0.
2287 IF( ASSOCIATED(bergso)) bergso = 0.
2288 IF( ASSOCIATED(melt) ) melt = 0.
2289 IF( ASSOCIATED(dnhet_ct) ) dnhet_ct = 0.
2290 IF( ASSOCIATED(dt_rasp) ) dt_rasp = 0.
2291
2292 IF( ASSOCIATED(qcres) ) qcres = 0.
2293 IF( ASSOCIATED(qires) ) qires = 0.
2294 IF( ASSOCIATED(autice) ) autice = 0.
2295 IF( ASSOCIATED(frzpp_ls) ) frzpp_ls = 0.
2296 IF( ASSOCIATED(snowmelt_ls) ) snowmelt_ls = 0.
2297 IF( ASSOCIATED(pfrz) ) pfrz = 0.
2298
2299 IF( ASSOCIATED(dncnuc) ) dncnuc = 0.
2300 IF( ASSOCIATED(dncsubl) ) dncsubl = 0.
2301 IF( ASSOCIATED(dnchmsplit) ) dnchmsplit = 0.
2302 IF( ASSOCIATED(dncautice) ) dncautice = 0.
2303 IF( ASSOCIATED(dncacris) ) dncacris = 0.
2304 IF( ASSOCIATED(dndccn) ) dndccn = 0.
2305 IF( ASSOCIATED(dndacrls) ) dndacrls = 0.
2306 IF( ASSOCIATED(dndacrlr) ) dndacrlr = 0.
2307 IF( ASSOCIATED(dndevapc) ) dndevapc = 0.
2308 IF( ASSOCIATED(dndautliq) ) dndautliq = 0.
2309! IF( ASSOCIATED(DNDCNV) ) DNDCNV = 0.
2310! IF( ASSOCIATED(DNCCNV) ) DNCCNV = 0.
2311
2312 end subroutine cloud_ptr_stubs
2313
2314!C=======================================================================
2319 REAL function erf_app(x)
2320 REAL :: x
2321 real*8:: aa(4), axx, y
2322 DATA aa /0.278393d0,0.230389d0,0.000972d0,0.078108d0/
2323
2324 y = dabs(dble(x))
2325 axx = 1.d0 + y*(aa(1)+y*(aa(2)+y*(aa(3)+y*aa(4))))
2326 axx = axx*axx
2327 axx = axx*axx
2328 axx = 1.d0 - (1.d0/axx)
2329 if(x.le.0.) then
2330 erf_app = sngl(-axx)
2331 else
2332 erf_app = sngl(axx)
2333 endif
2334 RETURN
2335 END FUNCTION
2336
2337 end module cldmacro
real function erf_app(x)
REAL FUNCTION erf (overwrites previous versions) THIS SUBROUTINE CALCULATES THE ERROR FUNCTION USING ...
Definition cldmacro.F:2320
subroutine pdffrac(flag, qtmean, sigmaqt1, sigmaqt2, qstar, clfrac)
This subroutine calculates.
Definition cldmacro.F:1105
subroutine, public fix_up_clouds_2m(qv, te, qlc, qic, cf, qla, qia, af, nl, ni, qc_min)
Definition cldmacro.F:623
subroutine pdfcondensate(flag, qtmean4, sigmaqt14, sigmaqt24, qstar4, condensate4, clfrac4)
This subroutine calculates.
Definition cldmacro.F:1159
real function icefrac(t, t_trans, t_freez)
Definition cldmacro.F:1865
subroutine micro_aa_bb_3(temp, pr, q_sat, aa, bb)
This subroutine calculates.
Definition cldmacro.F:1774
subroutine, public update_cld(irun, lm, dt, alpha, qc_min, pdfshape, pl, qv, qcl, qal, qci, qai, te, cf, af, scice, ni, nl)
This subroutine calculates the cloud fraction that would correspond to the current condensate.
Definition cldmacro.F:713
subroutine, public cloud_ptr_stubs(smaxl, smaxi, wsub, ccn01, ccn04, ccn1, nhet_nuc, nlim_nuc, so4, org, bcarbon, dust, seasalt, ncpl_vol, ncpi_vol, nrain, nsnow, cdnc_nuc, inc_nuc, sat_rat, qstot, qrtot, cldreffs, cldreffr, dqvdt_micro, dqidt_micro, dqldt_micro, dtdt_micro, rl_mask, ri_mask, kappa, sc_ice, cfice, cfliq, rhice, rhliq, rad_cf, rad_ql, rad_qi, rad_qs, rad_qr, rad_qv, cldreffi, cldreffl, nhet_imm, nhet_dep, nhet_dhf, dust_imm, dust_dep, dust_dhf, scf, scf_all, sigw_gw, sigw_cnv, sigw_turb, sigw_rc, rhcmicro, dnhet_imm, nondust_imm, nondust_dep, berg, bergso, melt, dnhet_ct, dtdt_macro, qcres, dt_rasp, frzpp_ls, snowmelt_ls, qires, autice, pfrz, dncnuc, dnchmsplit, dncsubl, dncautice, dncacris, dndccn, dndacrls, dndevapc, dndacrlr, dndautliq)
Definition cldmacro.F:2194
real function get_alhx3(t)
Definition cldmacro.F:1837
subroutine, public macro_cloud(irun, lm, dt, alf_fac, pp_dev, ppe_dev, qlwdtr_dev, th_dev, q_dev, qlw_ls_dev, qlw_an_dev, qiw_ls_dev, qiw_an_dev, anvfrc_dev, cldfrc_dev, physparams, sclmfdfr, alpht_dev, cnv_fice_dev, cnv_ndrop_dev, cnv_nice_dev, scice_dev, ncpl_dev, ncpi_dev, pfrz_dev, lprnt, ipr, rhc, pdfflag, qc_min)
This subroutine is the cloud macrophysics scheme in MG micriphysics.
Definition cldmacro.F:123
subroutine marshpalmq2(rain, pr, diam3, ntotal, w, ve)
This subroutine calculates.
Definition cldmacro.F:1716
subroutine pdf_spread(k, lm, pp, alpha, alpht_diag, frland, rhc)
Definition cldmacro.F:585
subroutine hystpdf(dt, alpha, pdfshape, qc_min, pl, qv, qcl, qal, qci, qai, te, cf, af, scice, ni, nl)
This subroutine calculates.
Definition cldmacro.F:813
subroutine bergeron_iter(dtime, pl, te, qv, qils, qicn,
This subroutine calculates.
Definition cldmacro.F:1902
subroutine cnvsrc(dt, iceparam, sclmfdfr, mass, imass, pl, te, qv, dcf, dmf, qla, qia, cf, af, qs, nl, ni, cnvfice, cnvndrop, cnvnice)
This subroutine calculates.
Definition cldmacro.F:1240
subroutine, public meltfrz_inst(im, lm, te, qcl, qal, qci, qai, nl, ni)
This subroutine calculates instantaneous freezing of condensate.
Definition cldmacro.F:2098
subroutine pfreezing(alpha, pl, te, qv, qcl, qal, qci, qai, sc_ice, cf, af, pf, pdfflag)
This calculates the probability of finding a supersaturated parcel in the grid cell....
Definition cldmacro.F:2025
real function ice_fraction(temp)
Definition cldmacro.F:1817
real function ldradius3(pl, te, qcl, nn)
Definition cldmacro.F:1798
subroutine precip3(k, lm, dt, frland, rhcr3, qpl, qpi, qcl, qci, te, qv, mass, imass, pl, dze, qddf3, aa, bb, area, rain, snow, pfl_above, pfi_above, evap_dd_above, subl_dd_above, revap_diag, rsubl_diag, acrll_diag, acril_diag, pfl_diag, pfi_diag, vfallrn, vfallsn, frz_diag, envfc, ddrfc, af, cf, pcbl, i)
This subroutine calculates.
Definition cldmacro.F:1323
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26
This module contains some of the most frequently used math and physics constants for GCM models.
Definition physcons.F90:36