CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
rascnv.F90
1
4
6 module rascnv
7
8 USE machine , ONLY : kind_phys
9 implicit none
10 public :: rascnv_init, rascnv_run
11 private
12 logical :: is_initialized = .false.
13!
14 integer, parameter :: kp = kind_phys
15 integer, parameter :: nrcmax=32 ! Maximum # of random clouds per 1200s
16
17 integer, parameter :: idnmax=999
18 real (kind=kind_phys), parameter :: delt_c=1800.0_kp/3600.0_kp &
19! Adjustment time scales in hrs for deep and shallow clouds
20! &, adjts_d=3.0, adjts_s=0.5
21! &, adjts_d=2.5, adjts_s=0.5
22 &, adjts_d=2.0_kp, adjts_s=0.5_kp
23!
24 logical, parameter :: fix_ncld_hr=.true.
25
26!
27 real (kind=kind_phys), parameter :: zero=0.0_kp, half=0.5_kp &
28 &, pt25=0.25_kp, one=1.0_kp &
29 &, two=2.0_kp, four=4.0_kp &
30 &, twoo3=two/3.0_kp &
31 &, four_p2=4.0e2_kp, one_m10=1.0e-10_kp&
32 &, one_m6=1.0e-6_kp, one_m5=1.0e-5_kp &
33 &, one_m2=1.0e-2_kp, one_m1=1.0e-1_kp &
34 &, oneolog10=one/log(10.0_kp) &
35 &, rain_min=1.0e-13_kp &
36 &, facmb=0.01_kp & ! conversion factor from Pa to hPa (or mb)
37 &, cmb2pa=100.0_kp ! Conversion from hPa to Pa
38!
39! real (kind=kind_phys), parameter :: frac=0.5_kp, crtmsf=0.0_kp &
40 real (kind=kind_phys), parameter :: frac=0.1_kp, crtmsf=0.0_kp &
41 &, tfrac_max=0.15_kp &
42 &, rhfacs=0.75_kp, rhfacl=0.75_kp &
43 &, face=5.0_kp, delx=10000.0_kp &
44 &, ddfac=face*delx*0.001_kp &
45 &, max_neg_bouy=0.15_kp &
46! &, max_neg_bouy=pt25_kp &
47 &, testmb=0.1_kp, testmbi=one/testmb &
48 &, dpd=0.5_kp, rknob=1.0_kp, eknob=1.0_kp
49
50!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
51! logical, parameter :: aw_scal=.false., cumfrc=.true. &
52 logical, parameter :: aw_scal=.true., cumfrc=.true. &
53 &, updret=.false., vsmooth=.false. &
54 &, wrkfun=.false., crtfun=.true. &
55 &, calkbl=.true., botop=.true., revap=.true. &
56 &, advcld=.true., advups=.false.,advtvd=.true.
57! &, advcld=.true., advups=.true., advtvd=.false.
58! &, advcld=.true., advups=.false.,advtvd=.false.
59
60
61 real(kind=kind_phys), parameter :: tf=233.16_kp, tcr=273.16_kp &
62 &, tcrf=one/(tcr-tf), tcl=2.0_kp
63
64!
65! For pressure gradient force in momentum mixing
66! real (kind=kind_phys), parameter :: pgftop=0.80, pgfbot=0.30 &
67! No pressure gradient force in momentum mixing
68 real (kind=kind_phys), parameter :: pgftop=0.0_kp, pgfbot=0.0_kp &
69! real (kind=kind_phys), parameter :: pgftop=0.55, pgfbot=0.55 &
70 &, pgfgrad=(pgfbot-pgftop)*0.001_kp&
71 &, cfmax=0.1_kp
72!
73! For Tilting Angle Specification
74!
75 real(kind=kind_phys) :: refp(6), refr(6), tlac(8), plac(8), &
76 tlbpl(7), drdp(5)
77!
78 DATA plac/100.0, 200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0/
79 DATA tlac/ 35.0, 25.0, 20.0, 17.5, 15.0, 12.5, 10.0, 7.5/
80 DATA refp/500.0, 300.0, 250.0, 200.0, 150.0, 100.0/
81 DATA refr/ 1.0, 2.0, 3.0, 4.0, 6.0, 8.0/
82!
83 real(kind=kind_phys) :: ac(16), ad(16)
84!
85 integer, parameter :: nqrp=500001
86 real(kind=kind_phys) :: c1xqrp, c2xqrp, tbqrp(nqrp), &
87 tbqra(nqrp), tbqrb(nqrp)
88!
89 integer, parameter :: nvtp=10001
90 real(kind=kind_phys) :: c1xvtp, c2xvtp, tbvtp(nvtp)
91!
92 real(kind=kind_phys) :: afc, facdt, &
93 grav, cp, alhl, alhf, rgas, rkap, nu, pi, &
94 t0c, rv, cvap, cliq, csol, ttp, eps, epsm1,&
95!
96 onebg, gravcon, onebcp, gravfac, elocp, &
97 elfocp, oneoalhl, cmpor, picon, zfac, &
98 deg2rad, piinv, testmboalhl, &
99 rvi, facw, faci, hsub, tmix, den
100
101 contains
102
103! -----------------------------------------------------------------------
104! CCPP entry points for gfdl cloud microphysics
105! -----------------------------------------------------------------------
106
113 subroutine rascnv_init(me, dt, con_g, con_cp, con_rd, &
114 con_rv, con_hvap, con_hfus, con_fvirt, &
115 con_t0c, con_ttp, con_cvap, con_cliq, &
116 con_csol, con_eps, con_epsm1, &
117 errmsg, errflg)
118!
119 Implicit none
120!
121 integer, intent(in) :: me
122 real(kind=kind_phys), intent(in) :: dt, &
123 con_g, con_cp, con_rd, con_rv, con_hvap, &
124 con_hfus, con_fvirt, con_t0c, con_cvap, con_cliq, &
125 con_csol, con_ttp, con_eps, con_epsm1
126
127 character(len=*), intent(out) :: errmsg
128 integer, intent(out) :: errflg
129!
130 real(kind=kind_phys), parameter :: actp=1.7_kp, facm=1.00_kp
131!
132 real(kind=kind_phys) :: ph(15), a(15)
133!
134 DATA ph/150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0, 500.0 &
135 &, 550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0/
136!
137 DATA a/ 1.6851, 1.1686, 0.7663, 0.5255, 0.4100, 0.3677 &
138 &, 0.3151, 0.2216, 0.1521, 0.1082, 0.0750, 0.0664 &
139 &, 0.0553, 0.0445, 0.0633/
140!
141 real(kind=kind_phys) tem, actop, tem1, tem2
142 integer i, l
143!
144! Initialize CCPP error handling variables
145 errmsg = ''
146 errflg = 0
147 if (is_initialized) return
148! set critical workfunction arrays
149 actop = actp*facm
150 DO l=1,15
151 a(l) = a(l)*facm
152 ENDDO
153 DO l=2,15
154 tem = one / (ph(l) - ph(l-1))
155 ac(l) = (ph(l)*a(l-1) - ph(l-1)*a(l)) * tem
156 ad(l) = (a(l) - a(l-1)) * tem
157 ENDDO
158 ac(1) = actop
159 ac(16) = a(15)
160 ad(1) = zero
161 ad(16) = zero
162!
163 CALL setqrp
164 CALL setvtp
165!
166 do i=1,7
167 tlbpl(i) = (tlac(i)-tlac(i+1)) / (plac(i)-plac(i+1))
168 enddo
169 do i=1,5
170 drdp(i) = (refr(i+1)-refr(i)) / (refp(i+1)-refp(i))
171 enddo
172!
173! VTP = 36.34*SQRT(1.2)* (0.001)**0.1364
174!
175 afc = -(1.01097e-4_kp*dt)*(3600.0_kp/dt)**0.57777778_kp
176!
177 if (fix_ncld_hr) then
178 facdt = delt_c / dt
179 else
180 facdt = one / 3600.0_kp
181 endif
182!
183 grav = con_g ; cp = con_cp ; alhl = con_hvap
184 alhf = con_hfus ; rgas = con_rd
185 nu = con_fvirt ; t0c = con_t0c
186 rv = con_rv ; cvap = con_cvap
187 cliq = con_cliq ; csol = con_csol ; ttp = con_ttp
188 eps = con_eps ; epsm1 = con_epsm1
189!
190 pi = four*atan(one) ; piinv = one/pi
191 onebg = one / grav ; gravcon = cmb2pa * onebg
192 onebcp = one / cp ; gravfac = grav / cmb2pa
193 rkap = rgas * onebcp ; deg2rad = pi/180.0_kp
194 elocp = alhl * onebcp ; elfocp = (alhl+alhf) * onebcp
195 oneoalhl = one/alhl ; cmpor = cmb2pa / rgas
196 picon = half*pi*onebg ; zfac = 0.28888889e-4_kp * onebg
197 testmboalhl = testmb/alhl
198!
199 rvi = one / rv ; facw = cvap - cliq
200 faci = cvap - csol ; hsub = alhl + alhf
201 tmix = ttp - 20.0_kp ; den = one / (ttp-tmix)
202!
203
204 if (me == 0) write(0,*) ' NO DOWNDRAFT FOR CLOUD TYPES' &
205 &, ' DETRAINING AT NORMALIZED PRESSURE ABOVE ',dpd
206!
207 is_initialized = .true.
208
209!
210 end subroutine rascnv_init
211!
212!!
213!!===================================================================== !
214!! rascnv_run: !
215!! !
216!! program history log: !
217!! Oct 2019 -- shrinivas moorthi !
218!! !
219!! !
220!! ==================== defination of variables ====================
221!! !
222!! !
223!! inputs: size
224!! !
225!! im - integer, horiz dimension and num of used pts 1 !
226!! k - integer, vertical dimension 1 !
227!! dt - real, time step in seconds 1 !
228!! dtf - real, dynamics time step in seconds 1 !
229!! rannum - real, array holding random numbers between 0 an 1 (im,nrcm) !
230!! tin - real, input temperature (K)
231!! qin - real, input specific humidity (kg/kg)
232!! uin - real, input zonal wind component
233!! vin - real, input meridional wind component
234!! ccin - real, input condensates+tracers
235!! fscav - real
236!! prsi - real, layer interface pressure
237!! prsl - real, layer mid pressure
238!! prsik - real, layer interface Exner function
239!! prslk - real, layer mid Exner function
240!! phil - real, layer mid geopotential height
241!! phii - real, layer interface geopotential height
242!! kpbl - integer pbl top index
243!! cdrag - real, drag coefficient
244!! rainc - real, convectinve rain (m/sec)
245!! kbot - integer, cloud bottom index
246!! ktop - integer, cloud top index
247!! knv - integer, 0 - no convvection; 1 - convection
248!! ddvel - downdraft induced surface wind
249!! flipv - logical, true if input data from bottom to top
250!! me - integer, current pe number
251!! area - real, grid area
252!! ccwf - real, multiplication factor for critical workfunction
253!! nrcm - integer, number of random numbers at each grid point
254!! rhc - real, critical relative humidity
255!! ud_mf - real, updraft mass flux
256!! dd_mf - real, downdraft mass flux
257!! dt_mf - real, detrained mass flux
258!! qw0 - real, min cloud water before autoconversion
259!! qi0 - real, min cloud ice before autoconversion
260!! dlqfac - real,fraction of condensated detrained in layers
261!! kdt - integer, current teime step
262!! revap - logial, when true reevaporate falling rain/snow
263!! qlcn - real
264!! qicn - real
265!! w_upi - real
266!! cf_upi - real
267!! cnv_mfd - real
268!! cnv_dqldt- real
269!! clcn - real
270!! cnv_fice - real
271!! cnv_ndrop- real
272!! cnv_nice - real
273!! mp_phys - integer, microphysics option
274!! mp_phys_mg - integer, flag for MG microphysics option
275!! trcmin - real, floor value for tracers
276!! ntk - integer, index representing TKE in the tracer array
277!!
278!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
279
284 subroutine rascnv_run(IM, k, itc, ntc, ntr, dt, dtf &
285 &, ccwf, area, dxmin, dxinv &
286 &, psauras, prauras, wminras, dlqf, flipv &
287 &, me, rannum, nrcm, mp_phys, mp_phys_mg &
288 &, ntk, kdt, rhc &
289 &, tin, qin, uin, vin, ccin, fscav &
290 &, prsi, prsl, prsik, prslk, phil, phii &
291 &, KPBL, CDRAG, RAINC, kbot, ktop, kcnv &
292 &, DDVEL, ud_mf, dd_mf, dt_mf &
293 &, QLCN, QICN, w_upi, cf_upi, CNV_MFD &
294 &, CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE &
295 &, errmsg, errflg)
296!
297!*********************************************************************
298!*********************************************************************
299!************ Relaxed Arakawa-Schubert ******************
300!************ Parameterization ******************
301!************ Plug Compatible Driver ******************
302!************ 23 May 2002 ******************
303!************ ******************
304!************ Developed By ******************
305!************ ******************
306!************ Shrinivas Moorthi ******************
307!************ ******************
308!************ EMC/NCEP ******************
309!*********************************************************************
310!*********************************************************************
311!
312!
313 Implicit none
314!
315! input
316!
317 logical, intent(in) :: flipv
318!
319 integer, intent(in) :: im, k, itc, ntc, ntr, me, nrcm, ntk, kdt &
320 &, mp_phys, mp_phys_mg
321 integer, dimension(:), intent(out) :: kbot, ktop
322 integer, dimension(:), intent(inout) :: kcnv
323 integer, dimension(:), intent(in) :: kpbl
324!
325 real(kind=kind_phys), intent(in) :: dxmin, dxinv, ccwf(:) &
326 &, psauras(:), prauras(:) &
327 &, wminras(:), dlqf(:)
328!
329 real(kind=kind_phys), dimension(:,:), intent(in) :: prsi, prsik, phii
330
331 real(kind=kind_phys), dimension(:,:), intent(inout) :: tin, qin, uin, vin
332 real(kind=kind_phys), dimension(:,:), intent(in) :: prsl, prslk, phil &
333 &, rhc
334 real(kind=kind_phys), dimension(:,:), intent(out), optional :: ud_mf
335 real(kind=kind_phys), dimension(:,:), intent(out) :: dd_mf, dt_mf
336 real(kind=kind_phys), dimension(:,:), intent(inout), optional :: qlcn, qicn, w_upi &
337 &, cnv_mfd &
338 &, cnv_dqldt, clcn &
339 &, cnv_fice, cnv_ndrop &
340 &, cnv_nice, cf_upi
341 real(kind=kind_phys), dimension(:) , intent(in) :: area, cdrag
342 real(kind=kind_phys), dimension(:) , intent(out) :: rainc
343 real(kind=kind_phys), dimension(:) , intent(out), optional :: ddvel
344 real(kind=kind_phys), dimension(:,:), intent(in) :: rannum
345 real(kind=kind_phys), intent(inout) :: ccin(:,:,:)
346 real(kind=kind_phys), intent(in) :: dt, dtf
347!
348! Added for aerosol scavenging for GOCART
349!
350 real(kind=kind_phys), intent(in) :: fscav(:)
351
352! &, ctei_r(im), ctei_rm
353 character(len=*), intent(out) :: errmsg
354 integer, intent(out) :: errflg
355!
356! locals
357!
358 real(kind=kind_phys) :: trcmin(ntr+2)
359 real(kind=kind_phys), dimension(k) :: toi, qoi, tcu, qcu &
360 &, pcu, clw, cli, qii, qli&
361 &, phi_l, prsm,psjm &
362 &, alfind, rhc_l &
363! &, alfinq, alfind, rhc_l &
364 &, qoi_l, qli_l, qii_l
365 real(kind=kind_phys), dimension(k+1) :: prs, psj, phi_h, flx, flxd
366
367
368 integer, dimension(100) :: ic
369 real(kind=kind_phys), parameter :: clwmin=1.0e-10_kp
370!
371 real(kind=kind_phys), allocatable :: alfint(:,:), uvi(:,:) &
372 &, trcfac(:,:), rcu(:,:)
373 real(kind=kind_phys) dtvd(2,4)
374! &, DPI(K)
375 real(kind=kind_phys) cfac, tem, sgc, ccwfac, tem1, tem2, rain &
376 &, wfnc,tla,pl,qiid,qlid, c0, c0i, dlq_fac, sumq&
377 &, rainp
378! integer :: nrcmax ! Maximum # of random clouds per 1200s
379!
380 Integer kcr, kfx, ncmx, nc, ktem, i, ii, lm1, l &
381 &, ntrc, ll, km1, kp1, ipt, lv, KBL, n &
382 &, KRMIN, KRMAX, KFMAX, kblmx, irnd,ib &
383 &, kblmn, ksfc, ncrnd
384 real(kind=kind_phys) sgcs(k)
385!
386! Scavenging related parameters
387!
388 real fscav_(ntr+2) ! Fraction scavenged per km
389!
390 fscav_ = -999.0_kp ! By default no scavenging
391 if (itc > 0 .and. ntc > 0) then
392 n = itc + ntc - 1
393 if (n <= ntr + 2) then
394 fscav_(itc:n) = fscav
395 else
396 errmsg = 'Error in rascnv_run: test ntr >= itc + ntc - 3 FAILED'
397 errflg = 1
398 return
399 end if
400 end if
401 trcmin = -99999.0_kp
402 if (ntk-2 > 0) trcmin(ntk-2) = 1.0e-4_kp
403
405
406 errmsg = ''
407 errflg = 0
408!
409 km1 = k - 1
410 kp1 = k + 1
411 if (flipv) then
412 ksfc = 1
413 else
414 ksfc = kp1
415 endif
416!
417 ntrc = ntr
418 IF (cumfrc) THEN
419 ntrc = ntrc + 2
420 ENDIF
421 if (ntrc > 0) then
422 if (.not. allocated(trcfac)) allocate (trcfac(k,ntrc))
423 if (.not. allocated(uvi)) allocate (uvi(k,ntrc))
424 if (.not. allocated(rcu)) allocate (rcu(k,ntrc))
425 do n=1, ntrc
426 do l=1,k
427 trcfac(l,n) = one ! For other tracers
428 rcu(l,n) = zero
429 enddo
430 enddo
431 endif
432!
433!!!!! initialization for microphysics ACheng
434 if(mp_phys == mp_phys_mg) then
435 do l=1,k
436 do i=1,im
437 qlcn(i,l) = zero
438 qicn(i,l) = zero
439 w_upi(i,l) = zero
440 cf_upi(i,l) = zero
441 cnv_mfd(i,l) = zero
442! CNV_PRC3(i,l) = zero
443 cnv_dqldt(i,l) = zero
444 clcn(i,l) = zero
445 cnv_fice(i,l) = zero
446 cnv_ndrop(i,l) = zero
447 cnv_nice(i,l) = zero
448 enddo
449 enddo
450 endif
451!
452 if (.not. allocated(alfint)) allocate(alfint(k,ntrc+4))
453!
454! call set_ras_afc(dt)
455! AFC = -(1.04E-4*DT)*(3600./DT)**0.578
456! AFC = -(1.01097E-4*DT)*(3600./DT)**0.57777778
457!
458 do l=1,k
459 do i=1,im
460 ud_mf(i,l) = zero
461 dd_mf(i,l) = zero
462 dt_mf(i,l) = zero
463 enddo
464 enddo
465 DO ipt=1,im
466
467 tem1 = max(zero, min(one, (log(area(ipt)) - dxmin) * dxinv))
468 tem2 = one - tem1
469 ccwfac = ccwf(1)*tem1 + ccwf(2)*tem2
470 dlq_fac = dlqf(1)*tem1 + dlqf(2)*tem2
471 tem = one + dlq_fac
472 c0i = (psauras(1)*tem1 + psauras(2)*tem2) * tem
473 c0 = (prauras(1)*tem1 + prauras(2)*tem2) * tem
474 if (ccwfac == zero) ccwfac = half
475!
476! ctei = .false.
477! if (ctei_r(ipt) > ctei_rm) ctei = .true.
478!
479! Compute NCRND :
480! if flipv is true, then input variables are from bottom
481! to top while RAS goes top to bottom
482!
483 tem = one / prsi(ipt,ksfc)
484
485 krmin = 1
486 krmax = km1
487 kfmax = krmax
488 kblmx = 1
489 kblmn = 1
490 sgcs(k) = one
491 DO l=1,km1
492 ll = l
493 if (flipv) ll = kp1 -l ! Input variables are bottom to top!
494 sgc = prsl(ipt,ll) * tem
495 sgcs(l) = sgc
496 IF (sgc <= 0.050_kp) krmin = l
497! IF (SGC <= 0.700_kp) KRMAX = L
498! IF (SGC <= 0.800_kp) KRMAX = L
499 IF (sgc <= 0.760_kp) krmax = l
500! IF (SGC <= 0.930_kp) KFMAX = L
501 IF (sgc <= 0.970_kp) kfmax = l ! Commented on 20060202
502! IF (SGC <= 0.700_kp) kblmx = L ! Commented on 20101015
503 IF (sgc <= 0.600_kp) kblmx = l !
504! IF (SGC <= 0.650_kp) kblmx = L ! Commented on 20060202
505 IF (sgc <= 0.980_kp) kblmn = l !
506 ENDDO
507 krmin = max(krmin,2)
508
509!
510 if (fix_ncld_hr) then
511!!! NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/1200) + 0.50001
512 ncrnd = min(nrcmax, (krmax-krmin+1)) * (dtf/1800) + 0.10001_kp
513! NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/1200) + 0.10001
514! NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/900) + 0.50001
515! NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/600) + 0.50001
516! NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/360) + 0.50001
517! & + 0.50001
518! NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * min(1.0,DTF/360) + 0.1
519 else
520 ncrnd = min(nrcmax, (krmax-krmin+1))
521 endif
522 ncrnd = min(nrcm,max(ncrnd, 1))
523!
524 kcr = min(k,krmax)
525 ktem = min(k,kfmax)
526 kfx = ktem - kcr
527
528 IF (kfx > 0) THEN
529 IF (botop) THEN
530 DO nc=1,kfx
531 ic(nc) = ktem + 1 - nc
532 ENDDO
533 ELSE
534 DO nc=kfx,1,-1
535 ic(nc) = ktem + 1 - nc
536 ENDDO
537 ENDIF
538 ENDIF
539!
540 ncmx = kfx + ncrnd
541 IF (ncrnd > 0) THEN
542 DO i=1,ncrnd
543 ii = mod(i-1,nrcm) + 1
544 irnd = (rannum(ipt,ii)-0.0005_kp)*(kcr-krmin+1)
545 ic(kfx+i) = irnd + krmin
546 ENDDO
547 ENDIF
548!
549 do l=1,k
550 clw(l) = zero
551 cli(l) = zero
552 ! to be zero i.e. no environmental condensate!!!
553 qii(l) = zero
554 qli(l) = zero
555! Initialize heating, drying, cloudiness etc.
556 tcu(l) = zero
557 qcu(l) = zero
558 pcu(l) = zero
559 flx(l) = zero
560 flxd(l) = zero
561 do n=1,ntrc
562 rcu(l,n) = zero
563 enddo
564 enddo
565 flx(kp1) = zero
566 flxd(kp1) = zero
567 rain = zero
568!
569 if (flipv) then ! Input variables are bottom to top!
570 do l=1,k
571 ll = kp1 - l
572 ! Transfer input prognostic data into local variable
573 toi(l) = tin(ipt,ll)
574 qoi(l) = qin(ipt,ll)
575
576 prsm(l) = prsl(ipt,ll) * facmb
577 psjm(l) = prslk(ipt,ll)
578 phi_l(l) = phil(ipt,ll)
579 rhc_l(l) = rhc(ipt,ll)
580!
581 if (ntrc > ntr) then ! CUMFRC is true
582 uvi(l,ntr+1) = uin(ipt,ll)
583 uvi(l,ntr+2) = vin(ipt,ll)
584 endif
585!
586 if (ntr > 0) then ! tracers such as O3, dust etc
587 do n=1,ntr
588 uvi(l,n) = ccin(ipt,ll,n+2)
589 if (abs(uvi(l,n)) < 1.0e-20_kp) uvi(l,n) = zero
590 enddo
591 endif
592 enddo
593 do l=1,kp1
594 ll = kp1 + 1 - l ! Input variables are bottom to top!
595 prs(ll) = prsi(ipt,l) * facmb
596 psj(ll) = prsik(ipt,l)
597 phi_h(ll) = phii(ipt,l)
598 enddo
599!
600 if (ccin(ipt,1,2) <= -998.0_kp) then ! input ice/water are together
601 do l=1,k
602 ll = kp1 -l
603 tem = ccin(ipt,ll,1) &
604 & * max(zero, min(one, (tcr-toi(l))*tcrf))
605 ccin(ipt,ll,2) = ccin(ipt,ll,1) - tem
606 ccin(ipt,ll,1) = tem
607 enddo
608 endif
609 if (advcld) then
610 do l=1,k
611 ll = kp1 -l ! Input variables are bottom to top!
612 qii(l) = ccin(ipt,ll,1)
613 qli(l) = ccin(ipt,ll,2)
614 enddo
615 endif
616 kbl = max(min(k, kp1-kpbl(ipt)), k/2)
617!
618 else ! Input variables are top to bottom!
619
620 do l=1,k
621 ! Transfer input prognostic data into local variable
622 toi(l) = tin(ipt,l)
623 qoi(l) = qin(ipt,l)
624
625 prsm(l) = prsl(ipt, l) * facmb
626 psjm(l) = prslk(ipt,l)
627 phi_l(l) = phil(ipt,l)
628 rhc_l(l) = rhc(ipt,l)
629!
630 if (ntrc > ntr) then ! CUMFRC is true
631 uvi(l,ntr+1) = uin(ipt,l)
632 uvi(l,ntr+2) = vin(ipt,l)
633 endif
634!
635 if (ntr > 0) then ! tracers such as O3, dust etc
636 do n=1,ntr
637 uvi(l,n) = ccin(ipt,l,n+2)
638 if (abs(uvi(l,n)) < 1.0e-20_kp) uvi(l,n) = zero
639 enddo
640 endif
641 enddo
642 DO l=1,kp1
643 prs(l) = prsi(ipt,l) * facmb
644 psj(l) = prsik(ipt,l)
645 phi_h(l) = phii(ipt,l)
646 ENDDO
647!
648 if (ccin(ipt,1,2) <= -998.0_kp) then ! input ice/water are together
649 do l=1,k
650 tem = ccin(ipt,l,1) &
651 & * max(zero, min(one, (tcr-toi(l))*tcrf))
652 ccin(ipt,l,2) = ccin(ipt,l,1) - tem
653 ccin(ipt,l,1) = tem
654 enddo
655 endif
656 if (advcld) then
657 do l=1,k
658 qii(l) = ccin(ipt,l,1)
659 qli(l) = ccin(ipt,l,2)
660 enddo
661 endif
662!
663 kbl = kpbl(ipt)
664!
665 endif ! end of if (flipv) then
666!
667! do l=k,kctop(1),-1
668!! DPI(L) = 1.0 / (PRS(L+1) - PRS(L))
669! enddo
670!
671! print *,' ipt=',ipt
672
673 if (advups) then ! For first order upstream for updraft
674 alfint(:,:) = one
675 elseif (advtvd) then ! TVD flux limiter scheme for updraft
676! alfint(:,:) = one
677 alfint(:,:) = half
678 l = krmin
679 lm1 = l - 1
680 dtvd(1,1) = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1) &
681 & + alhl*(qoi(l)-qoi(lm1))
682 dtvd(1,2) = qoi(l) - qoi(lm1)
683 dtvd(1,3) = qli(l) - qli(lm1)
684 dtvd(1,4) = qii(l) - qii(lm1)
685 do l=krmin+1,k
686 lm1 = l - 1
687
688! write(0,*)' toi=',toi(l),toi(lm1),' phi_l=',phi_l(l),phi_l(lm1)
689! &,' qoi=',qoi(l),qoi(lm1),' cp=',cp,' alhl=',alhl
690
691 dtvd(2,1) = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1) &
692 & + alhl*(qoi(l)-qoi(lm1))
693
694! write(0,*)' l=',l,' dtvd=',dtvd(:,1)
695
696 if (abs(dtvd(2,1)) > 1.0e-10_kp) then
697 tem1 = dtvd(1,1) / dtvd(2,1)
698 tem2 = abs(tem1)
699 alfint(l,1) = one - half*(tem1 + tem2)/(one + tem2) ! for h
700 endif
701
702! write(0,*)' alfint=',alfint(l,1),' l=',l,' ipt=',ipt
703
704 dtvd(1,1) = dtvd(2,1)
705!
706 dtvd(2,2) = qoi(l) - qoi(lm1)
707
708! write(0,*)' l=',l,' dtvd2=',dtvd(:,2)
709
710 if (abs(dtvd(2,2)) > 1.0e-10_kp) then
711 tem1 = dtvd(1,2) / dtvd(2,2)
712 tem2 = abs(tem1)
713 alfint(l,2) = one - half*(tem1 + tem2)/(one + tem2) ! for q
714 endif
715 dtvd(1,2) = dtvd(2,2)
716!
717 dtvd(2,3) = qli(l) - qli(lm1)
718
719! write(0,*)' l=',l,' dtvd3=',dtvd(:,3)
720
721 if (abs(dtvd(2,3)) > 1.0e-10_kp) then
722 tem1 = dtvd(1,3) / dtvd(2,3)
723 tem2 = abs(tem1)
724 alfint(l,3) = one - half*(tem1 + tem2)/(one + tem2) ! for ql
725 endif
726 dtvd(1,3) = dtvd(2,3)
727!
728 dtvd(2,4) = qii(l) - qii(lm1)
729
730! write(0,*)' l=',l,' dtvd4=',dtvd(:,4)
731
732 if (abs(dtvd(2,4)) > 1.0e-10_kp) then
733 tem1 = dtvd(1,4) / dtvd(2,4)
734 tem2 = abs(tem1)
735 alfint(l,4) = one - half*(tem1 + tem2)/(one + tem2) ! for qi
736 endif
737 dtvd(1,4) = dtvd(2,4)
738 enddo
739!
740 if (ntrc > 0) then
741 do n=1,ntrc
742 l = krmin
743 dtvd(1,1) = uvi(l,n) - uvi(l-1,n)
744 do l=krmin+1,k
745 dtvd(2,1) = uvi(l,n) - uvi(l-1,n)
746
747! write(0,*)' l=',l,' dtvdn=',dtvd(:,1),' n=',n,' l=',l
748
749 if (abs(dtvd(2,1)) > 1.0e-10_kp) then
750 tem1 = dtvd(1,1) / dtvd(2,1)
751 tem2 = abs(tem1)
752 alfint(l,n+4) = one - half*(tem1 + tem2)/(one + tem2) ! for tracers
753 endif
754 dtvd(1,1) = dtvd(2,1)
755 enddo
756 enddo
757 endif
758 else
759 alfint(:,:) = half ! For second order scheme
760 endif
761 alfind(:) = half
762!
763! write(0,*)' after alfint for ipt=',ipt
764
765! Resolution dependent press grad correction momentum mixing
766
767 if (cumfrc) then
768 do l=krmin,k
769 tem = one - max(pgfbot, min(pgftop, pgftop+pgfgrad*prsm(l)))
770 trcfac(l,ntr+1) = tem
771 trcfac(l,ntr+2) = tem
772 enddo
773 endif
774!
775! if (calkbl) kbl = k
776
777 if (calkbl) then
778 kbl = kblmn
779 else
780 kbl = min(kbl, kblmn)
781 endif
782!
783 DO nc=1,ncmx ! multi cloud loop
784!
785 ib = ic(nc) ! cloud top level index
786 if (ib > kbl-1) cycle
787!
788!****************************************************************************
789! if (advtvd) then ! TVD flux limiter scheme for updraft
790! l = ib
791! lm1 = l - 1
792! dtvd(1,1) = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1)
793! & + alhl*(qoi(l)-qoi(lm1))
794! dtvd(1,2) = qoi(l) - qoi(lm1)
795! dtvd(1,3) = qli(l) - qli(lm1)
796! dtvd(1,4) = qii(l) - qii(lm1)
797! do l=ib+1,k
798! lm1 = l - 1
799! dtvd(2,1) = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1)
800! & + alhl*(qoi(l)-qoi(lm1))
801! if (abs(dtvd(2,1)) > 1.0e-10) then
802! tem1 = dtvd(1,1) / dtvd(2,1)
803! tem2 = abs(tem1)
804! alfint(l,1) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2) ! for h
805! endif
806! dtvd(1,1) = dtvd(2,1)
807!
808! dtvd(2,2) = qoi(l) - qoi(lm1)
809! if (abs(dtvd(2,2)) > 1.0e-10) then
810! tem1 = dtvd(1,2) / dtvd(2,2)
811! tem2 = abs(tem1)
812! alfint(l,2) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2) ! for q
813! endif
814! dtvd(1,2) = dtvd(2,2)
815!
816! dtvd(2,3) = qli(l) - qli(lm1)
817! if (abs(dtvd(2,3)) > 1.0e-10) then
818! tem1 = dtvd(1,3) / dtvd(2,3)
819! tem2 = abs(tem1)
820! alfint(l,3) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2) ! for ql
821! endif
822! dtvd(1,3) = dtvd(2,3)
823!
824! dtvd(2,4) = qii(l) - qii(lm1)
825! if (abs(dtvd(2,4)) > 1.0e-10) then
826! tem1 = dtvd(1,4) / dtvd(2,4)
827! tem2 = abs(tem1)
828! alfint(l,4) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2) ! for qi
829! endif
830! dtvd(1,4) = dtvd(2,4)
831! enddo
832!
833! if (ntrc > 0) then
834! do n=1,ntrc
835! l = ib
836! dtvd(1,1) = uvi(l,n) - uvi(l-1,n)
837! do l=ib+1,k
838! dtvd(2,1) = uvi(l,n) - uvi(l-1,n)
839! if (abs(dtvd(2,1)) > 1.0e-10) then
840! tem1 = dtvd(1,1) / dtvd(2,1)
841! tem2 = abs(tem1)
842! alfint(l,n+4) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2) ! for tracers
843! endif
844! dtvd(1,1) = dtvd(2,1)
845! enddo
846! enddo
847! endif
848! endif
849!****************************************************************************
850!
851 wfnc = zero
852 do l=ib,kp1
853 flx(l) = zero
854 flxd(l) = zero
855 enddo
856!
857 tla = -10.0_kp
858!
859 qiid = qii(ib) ! cloud top level ice before convection
860 qlid = qli(ib) ! cloud top level water before convection
861!
862 rainp = rain
863
864 CALL cloud(k, kp1, ib, ntrc, kblmx, kblmn &
865 &, frac, max_neg_bouy, vsmooth, aw_scal &
866 &, revap, wrkfun, calkbl, crtfun &
867 &, dt, kdt, tla, dpd &
868 &, alfint, rhfacl, rhfacs, area(ipt) &
869 &, ccwfac, cdrag(ipt), trcfac &
870 &, alfind, rhc_l, phi_l, phi_h, prs, prsm,sgcs &
871 &, toi, qoi, uvi, qli, qii, kbl, ddvel(ipt) &
872 &, tcu, qcu, rcu, pcu, flx, flxd, rain, wfnc, fscav_ &
873 &, trcmin, ntk-2, c0, wminras(1), c0i, wminras(2) &
874 &, dlq_fac)
875! &, ctei)
876
877!
878 if (flipv) then
879 do l=ib,k
880 ll = kp1 -l ! Input variables are bottom to top!
881 ud_mf(ipt,ll) = ud_mf(ipt,ll) + flx(l+1)
882 dd_mf(ipt,ll) = dd_mf(ipt,ll) + flxd(l+1)
883 enddo
884 ll = kp1 - ib
885 dt_mf(ipt,ll) = dt_mf(ipt,ll) + flx(ib)
886
887 if (mp_phys == mp_phys_mg) then ! Anning Cheng for microphysics 11/14/2015
888
889 cnv_mfd(ipt,ll) = cnv_mfd(ipt,ll) + flx(ib)/dt
890
891!! CNV_DQLDT(ipt,ll) = CNV_DQLDT(ipt,ll)
892!! & + max(0.,(QLI(ib)+QII(ib)-qiid-qlid))/dt
893 cnv_dqldt(ipt,ll) = cnv_dqldt(ipt,ll) + flx(ib)* &
894 & max(0.,(qli(ib)+qii(ib)-qiid-qlid))/dt
895!! & max(0.,(QLI(ib)+QII(ib)))/dt/3.
896 if(flx(ib)<0) write(*,*)"AAA666", flx(ib),qli(ib),qii(ib) &
897 & ,ipt,ll
898 endif
899
900 else
901
902 do l=ib,k
903 ud_mf(ipt,l) = ud_mf(ipt,l) + flx(l+1)
904 dd_mf(ipt,l) = dd_mf(ipt,l) + flxd(l+1)
905 enddo
906 dt_mf(ipt,ib) = dt_mf(ipt,ib) + flx(ib)
907
908 if (mp_phys == mp_phys_mg) then ! Anning Cheng for microphysics 11/14/2015
909 cnv_mfd(ipt,ib) = cnv_mfd(ipt,ib) + flx(ib)/dt
910!! CNV_DQLDT(ipt,ib) = CNV_DQLDT(ipt,ib)
911!! & + max(0.,(QLI(ib)+QII(ib)-qiid-qlid))/dt
912 cnv_dqldt(ipt,ib) = cnv_dqldt(ipt,ib) + flx(ib)* &
913 & max(zero,(qli(ib)+qii(ib)-qiid-qlid))/dt
914!! & max(0.,(QLI(ib)+QII(ib)))/dt/3.
915 if(flx(ib)<0) write(*,*)"AAA666", flx(ib),qli(ib),qii(ib) &
916 & ,ipt,ib
917 endif
918
919 endif
920!
921!
922! Warning!!!!
923! ------------
924! By doing the following, CLOUD does not contain environmental
925! condensate!
926!
927 if (.not. advcld) then
928 do l=1,k
929 clw(l) = clw(l) + qli(l)
930 cli(l) = cli(l) + qii(l)
931 qli(l) = zero
932 qii(l) = zero
933 enddo
934 endif
935!
936 ENDDO ! End of the NC loop!
937!
938 rainc(ipt) = rain * 0.001_kp ! Output rain is in meters
939 if (rainc(ipt) < rain_min) rainc(ipt) = zero
940
941 ktop(ipt) = kp1
942 kbot(ipt) = 0
943
944 kcnv(ipt) = 0
945
946
947 do l=k,1,-1
948! qli(l) = max(qli(l), zero)
949! qii(l) = max(qii(l), zero)
950! clw(i) = max(clw(i), zero)
951! cli(i) = max(cli(i), zero)
952
953 if (sgcs(l) < 0.93_kp .and. abs(tcu(l)) > one_m10) then
954! if (sgcs(l) < 0.90_kp .and. tcu(l) .ne. zero) then
955! if (sgcs(l) < 0.85_kp .and. tcu(l) .ne. zero) then
956 kcnv(ipt) = 1
957 endif
958! New test for convective clouds ! added in 08/21/96
959 if (clw(l)+cli(l) > zero .OR. &
960 & qli(l)+qii(l) > clwmin) ktop(ipt) = l
961 enddo
962 do l=1,km1
963 if (clw(l)+cli(l) > zero .OR. &
964 & qli(l)+qii(l) > clwmin) kbot(ipt) = l
965 enddo
966!
967 if (flipv) then
968 do l=1,k
969 ll = kp1 - l
970 tin(ipt,ll) = toi(l) ! Temperature
971 qin(ipt,ll) = qoi(l) ! Specific humidity
972 uin(ipt,ll) = uvi(l,ntr+1) ! U momentum
973 vin(ipt,ll) = uvi(l,ntr+2) ! V momentum
974
975!! for 2M microphysics, always output these variables
976 if (mp_phys == mp_phys_mg) then
977 if (advcld) then
978 qlcn(ipt,ll) = max(qli(l)-ccin(ipt,ll,2), zero)
979 qicn(ipt,ll) = max(qii(l)-ccin(ipt,ll,1), zero)
980 cnv_fice(ipt,ll) = qicn(ipt,ll) &
981 & / max(1.0e-10_kp,qlcn(ipt,ll)+qicn(ipt,ll))
982 else
983 qlcn(ipt,ll) = qli(l)
984 qicn(ipt,ll) = qii(l)
985 cnv_fice(ipt,ll) = qii(l)/max(1.0e-10_kp,qii(l)+qli(l))
986 endif
987 cf_upi(ipt,ll) = max(zero,min(0.02_kp*log(one+ &
988 & 500.0_kp*ud_mf(ipt,ll)/dt), cfmax))
989! & 500*ud_mf(ipt,ll)/dt), 0.60))
990 clcn(ipt,ll) = cf_upi(ipt,ll) !downdraft is below updraft
991 w_upi(ipt,ll) = ud_mf(ipt,ll)*toi(l)*rgas / &
992 & (dt*max(cf_upi(ipt,ll),1.0e-12_kp)*prsl(ipt,ll))
993 endif
994
995 if (ntr > 0) then
996 do n=1,ntr
997 ccin(ipt,ll,n+2) = uvi(l,n) ! Tracers
998 enddo
999 endif
1000 enddo
1001 if (advcld) then
1002 do l=1,k
1003 ll = kp1 - l
1004 ccin(ipt,ll,1) = qii(l) ! Cloud ice
1005 ccin(ipt,ll,2) = qli(l) ! Cloud water
1006 enddo
1007 else
1008 do l=1,k
1009 ll = kp1 - l
1010 ccin(ipt,ll,1) = ccin(ipt,ll,1) + cli(l)
1011 ccin(ipt,ll,2) = ccin(ipt,ll,2) + clw(l)
1012 enddo
1013 endif
1014!
1015 ktop(ipt) = kp1 - ktop(ipt)
1016 kbot(ipt) = kp1 - kbot(ipt)
1017!
1018 else
1019
1020 do l=1,k
1021 tin(ipt,l) = toi(l) ! Temperature
1022 qin(ipt,l) = qoi(l) ! Specific humidity
1023 uin(ipt,l) = uvi(l,ntr+1) ! U momentum
1024 vin(ipt,l) = uvi(l,ntr+2) ! V momentum
1025
1026!! for 2M microphysics, always output these variables
1027 if (mp_phys == mp_phys_mg) then
1028 if (advcld) then
1029 qlcn(ipt,l) = max(qli(l)-ccin(ipt,l,2), zero)
1030 qicn(ipt,l) = max(qii(l)-ccin(ipt,l,1), zero)
1031 cnv_fice(ipt,l) = qicn(ipt,l) &
1032 & / max(1.0e-10_kp,qlcn(ipt,l)+qicn(ipt,l))
1033 else
1034 qlcn(ipt,l) = qli(l)
1035 qicn(ipt,l) = qii(l)
1036 cnv_fice(ipt,l) = qii(l)/max(1.0e-10_kp,qii(l)+qli(l))
1037 endif
1038!! CNV_PRC3(ipt,l) = PCU(l)/dt
1039! CNV_PRC3(ipt,l) = zero
1040! if(PCU(l) < zero) write(*,*)"AAA777",PCU(l),ipt,l
1041 cf_upi(ipt,l) = max(zero,min(0.02_kp*log(one+ &
1042 & 500.0_kp*ud_mf(ipt,l)/dt), cfmax))
1043! & 500*ud_mf(ipt,l)/dt), 0.60))
1044 clcn(ipt,l) = cf_upi(ipt,l) !downdraft is below updraft
1045 w_upi(ipt,l) = ud_mf(ipt,l)*toi(l)*rgas / &
1046 & (dt*max(cf_upi(ipt,l),1.0e-12_kp)*prsl(ipt,l))
1047 endif
1048
1049 if (ntr > 0) then
1050 do n=1,ntr
1051 ccin(ipt,l,n+2) = uvi(l,n) ! Tracers
1052 enddo
1053 endif
1054 enddo
1055 if (advcld) then
1056 do l=1,k
1057 ccin(ipt,l,1) = qii(l) ! Cloud ice
1058 ccin(ipt,l,2) = qli(l) ! Cloud water
1059 enddo
1060 else
1061 do l=1,k
1062 ccin(ipt,l,1) = ccin(ipt,l,1) + cli(l)
1063 ccin(ipt,l,2) = ccin(ipt,l,2) + clw(l)
1064 enddo
1065 endif
1066 endif
1067!
1068! Velocity scale from the downdraft!
1069!
1070 ddvel(ipt) = ddvel(ipt) * ddfac * grav / (prs(kp1)-prs(k))
1071!
1072 ENDDO ! End of the IPT Loop!
1073
1074 deallocate (alfint, uvi, trcfac, rcu)
1075!
1076 RETURN
1077 end subroutine rascnv_run
1078
1080 SUBROUTINE cloud( &
1081 & K, KP1, KD, NTRC, KBLMX, kblmn &
1082 &, FRACBL, MAX_NEG_BOUY, vsmooth, aw_scal &
1083 &, REVAP, WRKFUN, CALKBL, CRTFUN &
1084 &, DT, KDT, TLA, DPD &
1085 &, ALFINT, RHFACL, RHFACS, area, ccwf, cd, trcfac &
1086 &, alfind, rhc_ls, phil, phih, prs, prsm, sgcs &
1087 &, TOI, QOI, ROI, QLI, QII, KPBL, DSFC &
1088 &, TCU, QCU, RCU, PCU, FLX, FLXD, CUP, WFNC,fscav_ &
1089 &, trcmin, ntk, c0, qw0, c0i, qi0, dlq_fac)
1090! &, ctei)
1091
1092!
1093!***********************************************************************
1094!******************** Relaxed Arakawa-Schubert ************************
1095!****************** Plug Compatible Scalar Version *********************
1096!************************ SUBROUTINE CLOUD ****************************
1097!************************ October 2004 ****************************
1098!******************** VERSION 2.0 (modified) *************************
1099!************* Shrinivas.Moorthi@noaa.gov (301) 683-3718 ***** ********
1100!***********************************************************************
1101!*References:
1102!-----------
1103! NOAA Technical Report NWS/NCEP 99-01:
1104! Documentation of Version 2 of Relaxed-Arakawa-Schubert
1105! Cumulus Parameterization with Convective Downdrafts, June 1999.
1106! by S. Moorthi and M. J. Suarez.
1107!
1108! Relaxed Arakawa-Schubert Cumulus Parameterization (Version 2)
1109! with Convective Downdrafts - Unpublished Manuscript (2002)
1110! by Shrinivas Moorthi and Max J. Suarez.
1111!
1112!***********************************************************************
1113!
1114!===> UPDATES CLOUD TENDENCIES DUE TO A SINGLE CLOUD
1115!===> DETRAINING AT LEVEL KD.
1116!
1117!***********************************************************************
1118!
1119!===> TOI(K) INOUT TEMPERATURE KELVIN
1120!===> QOI(K) INOUT SPECIFIC HUMIDITY NON-DIMENSIONAL
1121!===> ROI(K,NTRC)INOUT TRACER ARBITRARY
1122!===> QLI(K) INOUT LIQUID WATER NON-DIMENSIONAL
1123!===> QII(K) INOUT ICE NON-DIMENSIONAL
1124
1125!===> PRS(KP1) INPUT PRESSURE @ EDGES MB
1126!===> PRSM(K) INPUT PRESSURE @ LAYERS MB
1127!===> SGCS(K) INPUT Local sigma
1128!===> PHIH(KP1) INPUT GEOPOTENTIAL @ EDGES IN MKS units
1129!===> PHIL(K) INPUT GEOPOTENTIAL @ LAYERS IN MKS units
1130!===> PRJ(KP1) INPUT (P/P0)^KAPPA @ EDGES NON-DIMENSIONAL
1131!===> PRJM(K) INPUT (P/P0)^KAPPA @ LAYERS NON-DIMENSIONAL
1132
1133!===> K INPUT THE RISE & THE INDEX OF THE SUBCLOUD LAYER
1134!===> KD INPUT DETRAINMENT LEVEL ( 1<= KD < K )
1135!===> NTRC INPUT NUMBER OF TRACERS. MAY BE ZERO.
1136!===> kblmx INPUT highest level the pbl can take
1137!===> kblmn INPUT lowest level the pbl can take
1138!===> DPD INPUT Critical normalized pressure (i.e. sigma) at the cloud top
1139! No downdraft calculation if the cloud top pressure is higher
1140! than DPD*PRS(KP1)
1141!
1142!===> TCU(K ) UPDATE TEMPERATURE TENDENCY DEG
1143!===> QCU(K ) UPDATE WATER VAPOR TENDENCY (G/G)
1144!===> RCU(K,NTRC)UPDATE TRACER TENDENCIES ND
1145!===> PCU(K) UPDATE PRECIP @ BASE OF LAYER KG/M^2
1146!===> FLX(K ) UPDATE MASS FLUX @ TOP OF LAYER KG/M^2
1147!===> CUP UPDATE PRECIPITATION AT THE SURFACE KG/M^2
1148!
1149 IMPLICIT NONE
1150!
1151 real (kind=kind_phys), parameter :: rhmax=1.0_kp & ! MAX RELATIVE HUMIDITY
1152 &, quad_lam=1.0_kp & ! MASK FOR QUADRATIC LAMBDA
1153 &, rhram=0.05_kp & ! PBL RELATIVE HUMIDITY RAMP
1154! &, RHRAM=0.15_kp !& ! PBL RELATIVE HUMIDITY RAMP
1155 &, hcritd=4000.0_kp & ! Critical Moist Static Energy for Deep clouds
1156! &, HCRITS=2000.0_kp & ! Critical Moist Static Energy for Shallow clouds
1157 &, hcrits=2500.0_kp & ! Critical Moist Static Energy for Shallow clouds
1158 &, pcrit_lcl=250.0_kp & ! Critical pressure difference between boundary layer top
1159 ! layer top and lifting condensation level (hPa)
1160! &, hpert_fac=1.01_kp !& ! Perturbation on hbl when ctei=.true.
1161! &, hpert_fac=1.005_kp !& ! Perturbation on hbl when ctei=.true.
1162 &, qudfac=quad_lam*half &
1163 &, shalfac=3.0_kp &
1164! &, qudfac=quad_lam*pt25, shalfac=3.0_kp !& ! Yogesh's
1165! &, c0ifac=0.07_kp & ! following Han et al, 2016 MWR
1166! &, c0ifac=0.001_kp & ! following Han et al, 2017 Weather and Forecasting
1167 &, c0ifac=0.0_kp &
1168 &, dpnegcr = 150.0_kp
1169! &, dpnegcr = 100.0_kp
1170! &, dpnegcr = 200.0_kp
1171!
1172 real(kind=kind_phys), parameter :: errmin=0.0001_kp &
1173 &, errmi2=0.1_kp*errmin &
1174! &, rainmin=1.0e-9_kp !&
1175 &, rainmin=1.0e-8_kp &
1176 &, oneopt9=one/0.09_kp &
1177 &, oneopt4=one/0.04_kp
1178 real(kind=kind_phys), parameter :: almax=1.0e-2_kp &
1179 &, almin1=0.0_kp, almin2=0.0_kp
1180 real(kind=kind_phys), parameter :: bldmax=300.0_kp, bldmin=25.0_kp
1181!
1182! INPUT ARGUMENTS
1183
1184! LOGICAL REVAP, WRKFUN, CALKBL, CRTFUN, CALCUP, ctei
1185 LOGICAL REVAP, WRKFUN, CALKBL, CRTFUN, CALCUP
1186 logical vsmooth, aw_scal
1187 INTEGER K, KP1, KD, NTRC, kblmx, kblmn, ntk
1188
1189
1190 real(kind=kind_phys), dimension(K) :: toi, qoi, prsm, qli, qii&
1191 &, phil, sgcs, rhc_ls &
1192 &, alfind
1193 real(kind=kind_phys), dimension(KP1) :: prs, phih
1194 real(kind=kind_phys), dimension(K,NTRC) :: roi, trcfac
1195 real(kind=kind_phys), dimension(ntrc) :: trcmin
1196 real(kind=kind_phys) :: cd, dsfc
1197 INTEGER :: KPBL, KBL, KB1, kdt
1198
1199 real(kind=kind_phys) alfint(k,ntrc+4)
1200 real(kind=kind_phys) fracbl, max_neg_bouy, dpd &
1201 &, rhfacl, rhfacs, area, ccwf &
1202 &, c0, qw0, c0i, qi0, dlq_fac
1203
1204! UPDATE ARGUMENTS
1205
1206 real(kind=kind_phys), dimension(K) :: tcu, qcu, tcd, qcd, pcu
1207 real(kind=kind_phys), dimension(KP1) :: flx, flxd
1208 real(kind=kind_phys), dimension(K,NTRC) :: rcu
1209 real(kind=kind_phys) :: cup
1210!
1211! TEMPORARY WORK SPACE
1212
1213 real(kind=kind_phys), dimension(KD:K) :: hol, qol, hst, qst &
1214 &, tol, gmh, akt, akc, bkc, ltl, rnn &
1215 &, fco, pri, qil, qll, zet, xi, rns &
1216 &, q0u, q0d, vtf, cil, cll, etai, dlq &
1217 &, wrk1, wrk2, dhdp, qrb, qrt, evp &
1218 &, ghd, gsd, etz, cldfr, sigf, rho
1219
1220 real(kind=kind_phys), dimension(KD:KP1) :: gaf, gms, gam, dlb &
1221 &, dlt, eta, prl, buy, etd, hod, qod, wvl
1222 real(kind=kind_phys), dimension(KD:K-1) :: etzi
1223
1224 real(kind=kind_phys) fscav_(ntrc)
1225
1226 LOGICAL ep_wfn, cnvflg, LOWEST, DDFT, UPDRET
1227
1228 real(kind=kind_phys) alm, det, hcc, clp &
1229 &, hsu, hsd, qtl, qtv &
1230 &, akm, wfn, hos, qos &
1231 &, amb, tx1, tx2, tx3 &
1232 &, tx4, tx5, qis, qls &
1233 &, hbl, qbl, rbl(ntrc), wcbase &
1234 &, qlb, qib, pris &
1235 &, wfnc, tx6, acr &
1236 &, tx7, tx8, tx9, rhc &
1237 &, hstkd, qstkd, ltlkd, q0ukd, q0dkd, dlbkd &
1238 &, qtp, qw00, qi00, qrbkd &
1239 &, hstold, rel_fac, prism &
1240 &, tl, pl, ql, qs, dqs, st1, sgn, tau, &
1241 & qtvp, hb, qb, tb, qqq, &
1242 & hccp, ds, dh, ambmax, x00, epp, qtlp, &
1243 & dpi, dphib, dphit, del_eta, detp, &
1244 & tem, tem1, tem2, tem3, tem4, &
1245 & st2, st3, st4, st5, &
1246 & errh, errw, erre, tem5, &
1247 & tem6, hbd, qbd, st1s, shal_fac, hmax, hmin, &
1248 & dhdpmn, avt, avq, avr, avh &
1249 &, train, dof, cldfrd, tla, gmf &
1250 &, fac, rsum1, rsum2, rsum3, dpneg, hcrit &
1251 &, actevap,arearat,deltaq,mass,massinv,potevap &
1252 &, teq,qsteq,dqdt,qeq &
1253 &, clfrac, dt, clvfr, delzkm, fnoscav, delp
1254! &, CLFRAC, DT, clf, clvfr, delzkm, fnoscav, delp
1255! &, almin1, almin2
1256
1257 INTEGER I, L, N, KD1, II, iwk, idh, lcon &
1258 &, IT, KM1, KTEM, KK, KK1, LM1, LL, LP1, kbls, kmxh &
1259 &, kblh, kblm, kblpmn, kmax, kmaxm1, kmaxp1, klcl, kmin, kmxb
1260!
1261!***********************************************************************
1262!
1263! almin2 = 0.2 * sqrt(pi/area)
1264! almin1 = almin2
1265
1266 km1 = k - 1
1267 kd1 = kd + 1
1268
1269 do l=1,k
1270 tcd(l) = zero
1271 qcd(l) = zero
1272 enddo
1273!
1274 cldfrd = zero
1275 dof = zero
1276 prl(kp1) = prs(kp1)
1277!
1278 DO l=kd,k
1279 rnn(l) = zero
1280 zet(l) = zero
1281 xi(l) = zero
1282!
1283 tol(l) = toi(l)
1284 qol(l) = qoi(l)
1285 prl(l) = prs(l)
1286 cll(l) = qli(l)
1287 cil(l) = qii(l)
1288 buy(l) = zero
1289
1290 wvl(l) = zero
1291 ENDDO
1292 wvl(kp1) = zero
1293!
1294 if (vsmooth) then
1295 do l=kd,k
1296 wrk1(l) = tol(l)
1297 wrk2(l) = qol(l)
1298 enddo
1299 do l=kd1,km1
1300 tol(l) = pt25*wrk1(l-1) + half*wrk1(l) + pt25*wrk1(l+1)
1301 qol(l) = pt25*wrk2(l-1) + half*wrk2(l) + pt25*wrk2(l+1)
1302 enddo
1303 endif
1304!
1305 DO l=kd, k
1306 dpi = one / (prl(l+1) - prl(l))
1307 pri(l) = gravfac * dpi
1308!
1309 pl = prsm(l)
1310 tl = tol(l)
1311
1312 rho(l) = cmb2pa * pl / (rgas*tl*(one+nu*qol(l)))
1313
1314 akt(l) = (prl(l+1) - pl) * dpi
1315!
1316 CALL qsatcn(tl, pl, qs, dqs)
1317!
1318 qst(l) = qs
1319 gam(l) = dqs * elocp
1320 st1 = one + gam(l)
1321 gaf(l) = oneoalhl * gam(l) / st1
1322
1323 ql = max(min(qs*rhmax,qol(l)), one_m10)
1324 qol(l) = ql
1325
1326 tem = cp * tl
1327 ltl(l) = tem * st1 / (one+nu*(qst(l)+tl*dqs))
1328 vtf(l) = one + nu * ql
1329 eta(l) = one / (ltl(l) * vtf(l))
1330
1331 hol(l) = tem + ql * alhl
1332 hst(l) = tem + qs * alhl
1333!
1334 ENDDO
1335!
1336 eta(kp1) = zero
1337 gms(k) = zero
1338!
1339 akt(kd) = half
1340 gms(kd) = zero
1341!
1342 clp = zero
1343!
1344 gam(kp1) = gam(k)
1345 gaf(kp1) = gaf(k)
1346!
1347 DO l=k,kd1,-1
1348 dphib = phil(l) - phih(l+1)
1349 dphit = phih(l) - phil(l)
1350!
1351 dlb(l) = dphib * eta(l) ! here eta contains 1/(L*(1+nu*q))
1352 dlt(l) = dphit * eta(l)
1353!
1354 qrb(l) = dphib
1355 qrt(l) = dphit
1356!
1357 eta(l) = eta(l+1) + dphib
1358
1359 hol(l) = hol(l) + eta(l)
1360 hstold = hst(l)
1361 hst(l) = hst(l) + eta(l)
1362!
1363 eta(l) = eta(l) + dphit
1364 ENDDO
1365!
1366! For the cloud top layer
1367!
1368 l = kd
1369
1370 dphib = phil(l) - phih(l+1)
1371!
1372 dlb(l) = dphib * eta(l)
1373!
1374 qrb(l) = dphib
1375 qrt(l) = dphib
1376!
1377 eta(l) = eta(l+1) + dphib
1378
1379 hol(l) = hol(l) + eta(l)
1380 hst(l) = hst(l) + eta(l)
1381!
1382! To determine KBL internally -- If KBL is defined externally
1383! the following two loop should be skipped
1384!
1385 if (sgcs(kd) < 0.5_kp) then
1386 hcrit = hcritd
1387 elseif (sgcs(kd) > 0.65_kp) then
1388 hcrit = hcrits
1389 else
1390 hcrit = (hcrits*(sgcs(kd)-0.5_kp) + hcritd*(0.65_kp-sgcs(kd)))&
1391 & * (one/0.15_kp)
1392 endif
1393 IF (calkbl) THEN
1394 ktem = max(kd+1, kblmx)
1395 hmin = hol(k)
1396 kmin = k
1397 do l=km1,kd,-1
1398 if (hmin > hol(l)) then
1399 hmin = hol(l)
1400 kmin = l
1401 endif
1402 enddo
1403 if (kmin == k) return
1404 hmax = hol(k)
1405 kmax = k
1406 do l=km1,ktem,-1
1407 if (hmax < hol(l)) then
1408 hmax = hol(l)
1409 kmax = l
1410 endif
1411 enddo
1412 kmxb = kmax
1413 if (kmax < kmin) then
1414 kmax = k
1415 kmxb = k
1416 hmax = hol(kmax)
1417 elseif (kmax < k) then
1418 do l=kmax+1,k
1419! if (abs(hol(kmax)-hol(l)) > half * hcrit) then
1420 if (abs(hol(kmax)-hol(l)) > hcrit) then
1421 kmxb = l - 1
1422 exit
1423 endif
1424 enddo
1425 endif
1426 kmaxm1 = kmax - 1
1427 kmaxp1 = kmax + 1
1428 kblpmn = kmax
1429!
1430 dhdp(kmax:k) = zero
1431 dhdpmn = dhdp(kmax)
1432 do l=kmaxm1,ktem,-1
1433 dhdp(l) = (hol(l)-hol(l+1)) / (prl(l+2)-prl(l))
1434 if (dhdp(l) < dhdpmn) then
1435 dhdpmn = dhdp(l)
1436 kblpmn = l + 1
1437 elseif (dhdp(l) > zero .and. l <= kmin) then
1438 exit
1439 endif
1440 enddo
1441 kbl = kmax
1442 if (kblpmn < kmax) then
1443 do l=kblpmn,kmaxm1
1444 if (hmax-hol(l) < half*hcrit) then
1445 kbl = l
1446 exit
1447 endif
1448 enddo
1449 endif
1450!
1451 klcl = kd1
1452 if (kmax > kd1) then
1453 do l=kmaxm1,kd1,-1
1454 if (hmax > hst(l)) then
1455 klcl = l+1
1456 exit
1457 endif
1458 enddo
1459 endif
1460
1461! if (klcl == kd .or. klcl < ktem) return
1462
1463! This is to handle mid-level convection from quasi-uniform h
1464
1465 if (kmax < kmxb) then
1466 kmax = max(kd1, min(kmxb,k))
1467 kmaxm1 = kmax - 1
1468 kmaxp1 = kmax + 1
1469 endif
1470
1471
1472! if (prl(Kmaxp1) - prl(klcl) > 250.0 ) return
1473
1474 ii = max(kbl,kd1)
1475 kbl = max(klcl,kd1)
1476 tem = min(50.0_kp,max(10.0_kp,(prl(kmaxp1)-prl(kd))*0.10_kp))
1477 if (prl(kmaxp1) - prl(ii) > tem .and. ii > kbl) kbl = ii
1478
1479 if (kbl .ne. ii) then
1480 if (prl(kmaxp1)-prl(kbl) > bldmax) kbl = max(kbl,ii)
1481 endif
1482 if (kbl < ii) then
1483 if (hol(ii)-hol(ii-1) > half*hcrit) kbl = ii
1484 endif
1485
1486 if (prl(kbl) - prl(klcl) > pcrit_lcl) return
1487!
1488! KBL = min(kmax, MAX(KBL,KBLMX))
1489 kbl = min(kblmn, max(kbl,kblmx))
1490! kbl = min(kblh,kbl)
1491!!!
1492! tem1 = max(prl(kP1)-prl(k), &
1493! & min((prl(kbl) - prl(kd))*0.05, 10.0))
1494!! & min((prl(kbl) - prl(kd))*0.05, 20.0))
1495!! & min((prl(kbl) - prl(kd))*0.05, 30.0))
1496! if (prl(kp1)-prl(kbl) < tem1) then
1497! KTEM = MAX(KD+1, KBLMX)
1498! do l=k,KTEM,-1
1499! tem = prl(kp1) - prl(l)
1500! if (tem > tem1) then
1501! kbl = min(kbl,l)
1502! exit
1503! endif
1504! enddo
1505! endif
1506! if (kbl == kblmx .and. kmax >= km1) kbl = k - 1
1507!!!
1508
1509 kpbl = kbl
1510
1511 ELSE
1512 kbl = kpbl
1513 ENDIF
1514!
1515 kbl = min(kmax,max(kbl,kd+2))
1516 kb1 = kbl - 1
1517!
1518 if (prl(kmaxp1)-prl(kbl) > bldmax .or. kb1 <= kd ) then
1519! & .or. PRL(Kmaxp1)-PRL(KBL) < bldmin) then
1520 return
1521 endif
1522!
1523 pris = one / (prl(kp1)-prl(kbl))
1524 prism = one / (prl(kmaxp1)-prl(kbl))
1525 tx1 = eta(kbl) ! geopotential height at KBL
1526!
1527 gms(kbl) = zero
1528 xi(kbl) = zero
1529 zet(kbl) = zero
1530!
1531 shal_fac = one
1532! if (prl(kbl)-prl(kd) < 300.0 .and. kmax == k) shal_fac = shalfac
1533 if (prl(kbl)-prl(kd) < 350.0_kp .and. kmax == k) shal_fac = shalfac
1534 DO l=kmax,kd,-1
1535 IF (l >= kbl) THEN
1536 eta(l) = (prl(kmaxp1)-prl(l)) * prism
1537 ELSE
1538 zet(l) = (eta(l) - tx1) * onebg
1539 xi(l) = zet(l) * zet(l) * (qudfac*shal_fac)
1540 eta(l) = zet(l) - zet(l+1)
1541 gms(l) = xi(l) - xi(l+1)
1542 ENDIF
1543 ENDDO
1544 if (kmax < k) then
1545 do l=kmaxp1,kp1
1546 eta(l) = zero
1547 enddo
1548 endif
1549!
1550 hbl = hol(kmax) * eta(kmax)
1551 qbl = qol(kmax) * eta(kmax)
1552 qlb = cll(kmax) * eta(kmax)
1553 qib = cil(kmax) * eta(kmax)
1554 tx1 = qst(kmax) * eta(kmax)
1555!
1556 DO l=kmaxm1,kbl,-1
1557 tem = eta(l) - eta(l+1)
1558 hbl = hbl + hol(l) * tem
1559 qbl = qbl + qol(l) * tem
1560 qlb = qlb + cll(l) * tem
1561 qib = qib + cil(l) * tem
1562 tx1 = tx1 + qst(l) * tem
1563 ENDDO
1564
1565! if (ctei .and. sgcs(kd) > 0.65) then
1566! hbl = hbl * hpert_fac
1567! qbl = qbl * hpert_fac
1568! endif
1569
1570! Find Min value of HOL in TX2
1571 tx2 = hol(kd)
1572 idh = kd1
1573 DO l=kd1,kb1
1574 IF (hol(l) < tx2) THEN
1575 tx2 = hol(l)
1576 idh = l ! Level of minimum moist static energy!
1577 ENDIF
1578 ENDDO
1579 idh = 1
1580! IDH = MAX(KD1, IDH)
1581 idh = max(kd, idh) ! Moorthi May, 31, 2019
1582!
1583 tem1 = hbl - hol(kd)
1584 tem = hbl - hst(kd1) - ltl(kd1) * nu *(qol(kd1)-qst(kd1))
1585 lowest = kd == kb1
1586
1587 lcon = kd
1588 do l=kb1,kd1,-1
1589 if (hbl >= hst(l)) then
1590 lcon = l
1591 exit
1592 endif
1593 enddo
1594!
1595 if (lcon == kd .or. kbl <= kd .or. prl(kbl)-prsm(lcon) > 150.0_kp) &
1596 & return
1597!
1598 tx1 = rhfacs - qbl / tx1 ! Average RH
1599
1600 cnvflg = (tem > zero .OR. (lowest .AND. tem1 >= zero)) &
1601 & .AND. tx1 < rhram
1602
1603 IF (.NOT. cnvflg) RETURN
1604!
1605 rhc = max(zero, min(one, exp(-20.0_kp*tx1) ))
1606!
1607 wcbase = 0.1_kp
1608 if (ntrc > 0) then
1609 DO n=1,ntrc
1610 rbl(n) = roi(kmax,n) * eta(kmax)
1611 ENDDO
1612 DO n=1,ntrc
1613 DO l=kmaxm1,kbl,-1
1614 rbl(n) = rbl(n) + roi(l,n)*(eta(l)-eta(l+1))
1615 ENDDO
1616 ENDDO
1617!
1618! if (ntk > 0 .and. aw_scal) then
1619 if (ntk > 0) then
1620 if (rbl(ntk) > zero) then
1621 wcbase = min(two, max(wcbase, sqrt(twoo3*rbl(ntk))))
1622! wcbase = min(one, max(wcbase, sqrt(twoo3*rbl(ntk))))
1623 endif
1624 endif
1625
1626 endif
1627!
1628 tx4 = zero
1629 tx5 = zero
1630!
1631 tx3 = qst(kbl) - gaf(kbl) * hst(kbl)
1632 DO l=kbl,k
1633 qil(l) = max(zero, min(one, (tcr-tcl-tol(l))*tcrf))
1634 ENDDO
1635!
1636 DO l=kb1,kd1,-1
1637 lp1 = l + 1
1638 tem = qst(l) - gaf(l) * hst(l)
1639 tem1 = (tx3 + tem) * half
1640 st2 = (gaf(l)+gaf(lp1)) * half
1641!
1642 fco(lp1) = tem1 + st2 * hbl
1643
1644 rnn(lp1) = zet(lp1) * tem1 + st2 * tx4
1645 gmh(lp1) = xi(lp1) * tem1 + st2 * tx5
1646!
1647 tx3 = tem
1648 tx4 = tx4 + eta(l) * hol(l)
1649 tx5 = tx5 + gms(l) * hol(l)
1650!
1651 qil(l) = max(zero, min(one, (tcr-tcl-tol(l))*tcrf))
1652 qll(lp1) = (half*alhf) * st2 * (qil(l)+qil(lp1)) + one
1653 ENDDO
1654!
1655! FOR THE CLOUD TOP -- L=KD
1656!
1657 l = kd
1658!
1659 lp1 = l + 1
1660 tem = qst(l) - gaf(l) * hst(l)
1661 tem1 = (tx3 + tem) * half
1662 st2 = (gaf(l)+gaf(lp1)) * half
1663!
1664 fco(lp1) = tem1 + st2 * hbl
1665 rnn(lp1) = zet(lp1) * tem1 + st2 * tx4
1666 gmh(lp1) = xi(lp1) * tem1 + st2 * tx5
1667!
1668 fco(l) = tem + gaf(l) * hbl
1669 rnn(l) = tem * zet(l) + (tx4 + eta(l)*hol(l)) * gaf(l)
1670 gmh(l) = tem * xi(l) + (tx5 + gms(l)*hol(l)) * gaf(l)
1671!
1672! Replace FCO for the Bottom
1673!
1674 fco(kbl) = qbl
1675 rnn(kbl) = zero
1676 gmh(kbl) = zero
1677!
1678 qil(kd) = max(zero, min(one, (tcr-tcl-tol(kd))*tcrf))
1679 qll(kd1) = (half*alhf) * st2 * (qil(kd) + qil(kd1)) + one
1680 qll(kd ) = alhf * gaf(kd) * qil(kd) + one
1681!
1682 st1 = qil(kd)
1683 st2 = c0i * st1
1684 if (c0ifac > 1.0e-6_kp) st2 = st2 * exp(c0ifac*min(tol(kd)-t0c,zero))
1685 tem = c0 * (one-st1)
1686 tem2 = st2*qi0 + tem*qw0
1687!
1688 DO l=kd,kb1
1689 lp1 = l + 1
1690 tx2 = akt(l) * eta(l)
1691 tx1 = tx2 * tem2
1692 q0u(l) = tx1
1693 fco(l) = fco(lp1) - fco(l) + tx1
1694 rnn(l) = rnn(lp1) - rnn(l) &
1695 & + eta(l)*(qol(l)+cll(l)+cil(l)) + tx1*zet(l)
1696 gmh(l) = gmh(lp1) - gmh(l) &
1697 & + gms(l)*(qol(l)+cll(l)+cil(l)) + tx1*xi(l)
1698!
1699 tem1 = (one-akt(l)) * eta(l)
1700
1701 akt(l) = qll(l) + (st2 + tem) * tx2
1702
1703 akc(l) = one / akt(l)
1704!
1705 st1 = half * (qil(l)+qil(lp1))
1706 st2 = c0i * st1
1707 if (c0ifac > 1.0e-6_kp) st2 = st2 * exp(c0ifac*min(tol(lp1)-t0c,zero))
1708 tem = c0 * (one-st1)
1709 tem2 = st2*qi0 + tem*qw0
1710!
1711 bkc(l) = qll(lp1) - (st2 + tem) * tem1
1712!
1713 tx1 = tem1*tem2
1714 q0d(l) = tx1
1715 fco(l) = fco(l) + tx1
1716 rnn(l) = rnn(l) + tx1*zet(lp1)
1717 gmh(l) = gmh(l) + tx1*xi(lp1)
1718 ENDDO
1719
1720 qw00 = qw0
1721 qi00 = qi0
1722 ii = 0
1723 777 continue
1724!
1725 ep_wfn = .false.
1726 rnn(kbl) = zero
1727 tx3 = bkc(kb1) * (qib + qlb)
1728 tx4 = zero
1729 tx5 = zero
1730 DO l=kb1,kd1,-1
1731 tem = bkc(l-1) * akc(l)
1732 tx3 = (tx3 + fco(l)) * tem
1733 tx4 = (tx4 + rnn(l)) * tem
1734 tx5 = (tx5 + gmh(l)) * tem
1735 ENDDO
1736 IF (kd < kb1) THEN
1737 hsd = hst(kd1) + ltl(kd1) * nu *(qol(kd1)-qst(kd1))
1738 ELSE
1739 hsd = hbl
1740 ENDIF
1741!
1742 tx3 = (tx3 + fco(kd)) * akc(kd)
1743 tx4 = (tx4 + rnn(kd)) * akc(kd)
1744 tx5 = (tx5 + gmh(kd)) * akc(kd)
1745 alm = alhf*qil(kd) - ltl(kd) * vtf(kd)
1746!
1747 hsu = hst(kd) + ltl(kd) * nu * (qol(kd)-qst(kd))
1748!
1749!===> VERTICAL INTEGRALS NEEDED TO COMPUTE THE ENTRAINMENT PARAMETER
1750!
1751 tx1 = alm * tx4
1752 tx2 = alm * tx5
1753
1754 DO l=kd,kb1
1755 tau = hol(l) - hsu
1756 tx1 = tx1 + tau * eta(l)
1757 tx2 = tx2 + tau * gms(l)
1758 ENDDO
1759!
1760! MODIFY HSU TO INCLUDE CLOUD LIQUID WATER AND ICE TERMS
1761!
1762 hsu = hsu - alm * tx3
1763!
1764 clp = zero
1765 alm = -100.0_kp
1766 hos = hol(kd)
1767 qos = qol(kd)
1768 qis = cil(kd)
1769 qls = cll(kd)
1770
1771 cnvflg = hbl > hsu .and. abs(tx1) > 1.0e-4_kp
1772
1773!***********************************************************************
1774
1775 st1 = half*(hsu + hsd)
1776
1777 IF (cnvflg) THEN
1778!
1779! STANDARD CASE:
1780! CLOUD CAN BE NEUTRALLY BOUYANT AT MIDDLE OF LEVEL KD W/ +VE LAMBDA.
1781! EPP < .25 IS REQUIRED TO HAVE REAL ROOTS.
1782!
1783 clp = one
1784 st2 = hbl - hsu
1785!
1786 if (tx2 == zero) then
1787 alm = - st2 / tx1
1788 if (alm > almax) alm = -100.0_kp
1789 else
1790 x00 = tx2 + tx2
1791 epp = tx1 * tx1 - (x00+x00)*st2
1792 if (epp > zero) then
1793 x00 = one / x00
1794 tem = sqrt(epp)
1795 tem1 = (-tx1-tem)*x00
1796 tem2 = (-tx1+tem)*x00
1797 if (tem1 > almax) tem1 = -100.0_kp
1798 if (tem2 > almax) tem2 = -100.0_kp
1799 alm = max(tem1,tem2)
1800
1801 endif
1802 endif
1803!
1804! CLIP CASE:
1805! NON-ENTRAINIG CLOUD DETRAINS IN LOWER HALF OF TOP LAYER.
1806! NO CLOUDS ARE ALLOWED TO DETRAIN BELOW THE TOP LAYER.
1807!
1808 ELSEIF (hbl <= hsu .AND. hbl > st1) THEN
1809 alm = zero
1810! CLP = (HBL-ST1) / (HSU-ST1) ! commented on Jan 16, 2010
1811 ENDIF
1812!
1813 cnvflg = .true.
1814 IF (almin1 > zero) THEN
1815 IF (alm >= almin1) cnvflg = .false.
1816 ELSE
1817 lowest = kd == kb1
1818 IF ( (alm > zero) .OR. &
1819 & (.NOT. lowest .AND. alm == zero) ) cnvflg = .false.
1820 ENDIF
1821!
1822!===> IF NO SOUNDING MEETS SECOND CONDITION, RETURN
1823!
1824 IF (cnvflg) THEN
1825 IF (ii > 0 .or. (qw00 == zero .and. qi00 == zero)) RETURN
1826 clp = one
1827 ep_wfn = .true.
1828 GO TO 888
1829 ENDIF
1830!
1831 st1s = one
1832 IF(clp > zero .AND. clp < one) THEN
1833 st1 = half*(one+clp)
1834 st2 = one - st1
1835 st1s = st1
1836 hstkd = hst(kd)
1837 qstkd = qst(kd)
1838 ltlkd = ltl(kd)
1839 q0ukd = q0u(kd)
1840 q0dkd = q0d(kd)
1841 dlbkd = dlb(kd)
1842 qrbkd = qrb(kd)
1843!
1844 hst(kd) = hst(kd)*st1 + hst(kd1)*st2
1845 hos = hol(kd)*st1 + hol(kd1)*st2
1846 qst(kd) = qst(kd)*st1 + qst(kd1)*st2
1847 qos = qol(kd)*st1 + qol(kd1)*st2
1848 qls = cll(kd)*st1 + cll(kd1)*st2
1849 qis = cil(kd)*st1 + cil(kd1)*st2
1850 ltl(kd) = ltl(kd)*st1 + ltl(kd1)*st2
1851!
1852 dlb(kd) = dlb(kd)*clp
1853 qrb(kd) = qrb(kd)*clp
1854 eta(kd) = eta(kd)*clp
1855 gms(kd) = gms(kd)*clp
1856 q0u(kd) = q0u(kd)*clp
1857 q0d(kd) = q0d(kd)*clp
1858 ENDIF
1859!
1860!
1861!***********************************************************************
1862!
1863! Critical workfunction is included in this version
1864!
1865 acr = zero
1866 tem = prl(kd1) - (prl(kd1)-prl(kd)) * clp * half
1867 tx1 = prl(kbl) - tem
1868 tx2 = min(900.0_kp, max(tx1,100.0_kp))
1869 tem1 = log(tx2*0.01_kp) * oneolog10
1870 tem2 = one - tem1
1871 if ( kdt == 1 ) then
1872! rel_fac = (dt * facdt) / (tem1*12.0_kp + tem2*3.0)
1873 rel_fac = (dt * facdt) / (tem1*6.0_kp + tem2*adjts_s)
1874 else
1875 rel_fac = (dt * facdt) / (tem1*adjts_d + tem2*adjts_s)
1876 endif
1877!
1878! rel_fac = max(zero, min(one,rel_fac))
1879 rel_fac = max(zero, min(half,rel_fac))
1880
1881 IF (crtfun) THEN
1882 iwk = tem*0.02_kp - 0.999999999_kp
1883 iwk = max(1, min(iwk, 16))
1884 acr = tx1 * (ac(iwk) + tem * ad(iwk)) * ccwf
1885 ENDIF
1886!
1887!===> NORMALIZED MASSFLUX
1888!
1889! ETA IS THE THICKNESS COMING IN AND normalized MASS FLUX GOING OUT.
1890! GMS IS THE THICKNESS SQUARE ; IT IS LATER REUSED FOR GAMMA_S
1891!
1892! ETA(K) = ONE
1893
1894 DO l=kb1,kd,-1
1895 eta(l) = eta(l+1) + alm * (eta(l) + alm * gms(l))
1896 etai(l) = one / eta(l)
1897 ENDDO
1898 etai(kbl) = one
1899!
1900!===> CLOUD WORKFUNCTION
1901!
1902 wfn = zero
1903 akm = zero
1904 det = zero
1905 hcc = hbl
1906 cnvflg = .false.
1907 qtl = qst(kb1) - gaf(kb1)*hst(kb1)
1908 tx1 = hbl
1909!
1910 qtv = qbl
1911 det = qlb + qib
1912!
1913 tx2 = zero
1914 dpneg = zero
1915!
1916 DO l=kb1,kd1,-1
1917 lm1 = l - 1
1918 lp1 = l + 1
1919 del_eta = eta(l) - eta(lp1)
1920 hccp = hcc + del_eta*hol(l)
1921!
1922 qtlp = qst(lm1) - gaf(lm1)*hst(lm1)
1923 qtvp = half * ((qtlp+qtl)*eta(l) &
1924 & + (gaf(l)+gaf(lm1))*hccp)
1925 st1 = eta(l)*q0u(l) + eta(lp1)*q0d(l)
1926 detp = (bkc(l)*det - (qtvp-qtv) &
1927 & + del_eta*(qol(l)+cll(l)+cil(l)) + st1) * akc(l)
1928
1929 tem1 = akt(l) - qll(l)
1930 tem2 = qll(lp1) - bkc(l)
1931 rns(l) = tem1*detp + tem2*det - st1
1932
1933 qtp = half * (qil(l)+qil(lm1))
1934 tem2 = min(qtp*(detp-eta(l)*qw00), &
1935 & (one-qtp)*(detp-eta(l)*qi00))
1936 st1 = min(tx2,tem2)
1937 tx2 = tem2
1938!
1939 IF (rns(l) < zero .or. st1 < zero) ep_wfn = .true.
1940 IF (detp <= zero) cnvflg = .true.
1941
1942 st1 = hst(l) - ltl(l)*nu*(qst(l)-qol(l))
1943
1944
1945 tem2 = hccp + detp * qtp * alhf
1946!
1947 st2 = ltl(l) * vtf(l)
1948 tem5 = cll(l) + cil(l)
1949 tem3 = (tx1 - eta(lp1)*st1 - st2*(det-tem5*eta(lp1))) * dlb(l)
1950 tem4 = (tem2 - eta(l )*st1 - st2*(detp-tem5*eta(l))) * dlt(l)
1951!
1952 st1 = tem3 + tem4
1953
1954 wfn = wfn + st1
1955 akm = akm - min(st1,zero)
1956
1957 if (st1 < zero .and. wfn < zero) then
1958 dpneg = dpneg + prl(lp1) - prl(l)
1959 endif
1960
1961 buy(l) = half * (tem3/(eta(lp1)*qrb(l)) + tem4/(eta(l)*qrt(l)))
1962!
1963 hcc = hccp
1964 det = detp
1965 qtl = qtlp
1966 qtv = qtvp
1967 tx1 = tem2
1968
1969 ENDDO
1970
1971 del_eta = eta(kd) - eta(kd1)
1972 hccp = hcc + del_eta*hos
1973!
1974 qtlp = qst(kd) - gaf(kd)*hst(kd)
1975 qtvp = qtlp*eta(kd) + gaf(kd)*hccp
1976 st1 = eta(kd)*q0u(kd) + eta(kd1)*q0d(kd)
1977 detp = (bkc(kd)*det - (qtvp-qtv) &
1978 & + del_eta*(qos+qls+qis) + st1) * akc(kd)
1979!
1980 tem1 = akt(kd) - qll(kd)
1981 tem2 = qll(kd1) - bkc(kd)
1982 rns(kd) = tem1*detp + tem2*det - st1
1983!
1984 IF (rns(kd) < zero) ep_wfn = .true.
1985 IF (detp <= zero) cnvflg = .true.
1986!
1987 888 continue
1988
1989 if (ep_wfn) then
1990 IF ((qw00 == zero .and. qi00 == zero)) RETURN
1991 if (ii == 0) then
1992 ii = 1
1993 if (clp > zero .and. clp < one) then
1994 hst(kd) = hstkd
1995 qst(kd) = qstkd
1996 ltl(kd) = ltlkd
1997 q0u(kd) = q0ukd
1998 q0d(kd) = q0dkd
1999 dlb(kd) = dlbkd
2000 qrb(kd) = qrbkd
2001 endif
2002 do l=kd,kb1
2003 lp1 = l + 1
2004 fco(l) = fco(l) - q0u(l) - q0d(l)
2005 rnn(l) = rnn(l) - q0u(l)*zet(l) - q0d(l)*zet(lp1)
2006 gmh(l) = gmh(l) - q0u(l)*xi(l) - q0d(l)*zet(lp1)
2007 eta(l) = zet(l) - zet(lp1)
2008 gms(l) = xi(l) - xi(lp1)
2009 q0u(l) = zero
2010 q0d(l) = zero
2011 ENDDO
2012 qw00 = zero
2013 qi00 = zero
2014
2015 go to 777
2016 else
2017 cnvflg = .true.
2018 endif
2019 endif
2020!
2021!
2022! ST1 = 0.5 * (HST(KD) - LTL(KD)*NU*(QST(KD)-QOS) &
2023! & + HST(KD1) - LTL(KD1)*NU*(QST(KD1)-QOL(KD1)))
2024!
2025 st1 = hst(kd) - ltl(kd)*nu*(qst(kd)-qos)
2026 st2 = ltl(kd) * vtf(kd)
2027 tem5 = (qls + qis) * eta(kd1)
2028 st1 = half * (tx1-eta(kd1)*st1-st2*(det-tem5))*dlb(kd)
2029!
2030 wfn = wfn + st1
2031 akm = akm - min(st1,zero) ! Commented on 08/26/02 - does not include top
2032!
2033
2034 buy(kd) = st1 / (eta(kd1)*qrb(kd))
2035!
2036 det = detp
2037 hcc = hccp
2038 akm = akm / wfn
2039
2040
2041!***********************************************************************
2042!
2043 IF (wrkfun) THEN ! If only to calculate workfunction save it and return
2044 IF (wfn >= zero) wfnc = wfn
2045 RETURN
2046 ELSEIF (.NOT. crtfun) THEN
2047 acr = wfnc
2048 ENDIF
2049!
2050!===> THIRD CHECK BASED ON CLOUD WORKFUNCTION
2051!
2052 calcup = .false.
2053
2054 tem = max(0.05_kp, min(cd*200.0_kp, max_neg_bouy))
2055 IF (.not. cnvflg .and. wfn > acr .and. &
2056 & dpneg < dpnegcr .and. akm <= tem) calcup = .true.
2057!
2058!===> IF NO SOUNDING MEETS THIRD CONDITION, RETURN
2059!
2060 IF (.NOT. calcup) RETURN
2061!
2062! This is for not LL - 20050601
2063! IF (ALMIN2 .NE. zero) THEN
2064! IF (ALMIN1 .NE. ALMIN2) ST1 = one / max(ONE_M10,(ALMIN2-ALMIN1))
2065! IF (ALM < ALMIN2) THEN
2066! CLP = CLP * max(zero, min(one,(0.3 + 0.7*(ALM-ALMIN1)*ST1)))
2067!! CLP = CLP * max(0.0, min(1.0,(0.2 + 0.8*(ALM-ALMIN1)*ST1)))
2068!! CLP = CLP * max(0.0, min(1.0,(0.1 + 0.9*(ALM-ALMIN1)*ST1)))
2069! ENDIF
2070! ENDIF
2071!
2072 clp = clp * rhc
2073 dlq = zero
2074 tem = one / (one + dlq_fac)
2075 do l=kd,kb1
2076 rnn(l) = rns(l) * tem
2077 dlq(l) = rns(l) * tem * dlq_fac
2078 enddo
2079 DO l=kbl,k
2080 rnn(l) = zero
2081 ENDDO
2082!
2083! If downdraft is to be invoked, do preliminary check to see
2084! if enough rain is available and then call DDRFT.
2085!
2086 ddft = .false.
2087 IF (dpd > zero) THEN
2088 train = zero
2089 IF (clp > zero) THEN
2090 DO l=kd,kb1
2091 train = train + rnn(l)
2092 ENDDO
2093 ENDIF
2094
2095 pl = (prl(kd1) + prl(kd))*half
2096 IF (train > 1.0e-4_kp .AND. pl <= dpd*prl(kp1)) ddft = .true.
2097 ENDIF
2098!
2099 IF (ddft) THEN ! Downdraft scheme based on (Cheng and Arakawa, 1997)
2100 CALL ddrft( &
2101 & k, kp1, kd &
2102 &, tla, alfind, wcbase &
2103 &, tol, qol, hol, prl, qst, hst, gam, gaf &
2104! &, TOL, QOL, HOL, PRL, QST, HST, GAM, GAF, HBL, QBL &
2105 &, qrb, qrt, buy, kbl, idh, eta, rnn, etai &
2106 &, alm, wfn, train, ddft &
2107 &, etd, hod, qod, evp, dof, cldfr, etz &
2108 &, gms, gsd, ghd, wvl)
2109
2110 ENDIF
2111!
2112! No Downdraft case (including case with no downdraft solution)
2113! ---------------------------------------------------------
2114!
2115 IF (.NOT. ddft) THEN
2116 DO l=kd,kp1
2117 etd(l) = zero
2118 hod(l) = zero
2119 qod(l) = zero
2120 wvl(l) = zero
2121 ENDDO
2122 DO l=kd,k
2123 evp(l) = zero
2124 etz(l) = zero
2125 ENDDO
2126
2127 ENDIF
2128
2129!
2130!===> CALCULATE GAMMAS i.e. TENDENCIES PER UNIT CLOUD BASE MASSFLUX
2131! Includes downdraft terms!
2132
2133 avh = zero
2134
2135!
2136! Fraction of detrained condensate evaporated
2137!
2138! tem1 = max(ZERO, min(HALF, (prl(kd)-FOUR_P2)*ONE_M2))
2139! tem1 = max(ZERO, min(HALF, (prl(kd)-300.0)*0.005))
2140 tem1 = zero
2141! tem1 = 1.0
2142! if (kd1 == kbl) tem1 = 0.0
2143!
2144 tem2 = one - tem1
2145 tem = det * qil(kd)
2146
2147
2148 st1 = (hcc+alhf*tem-eta(kd)*hst(kd)) / (one+gam(kd))
2149 ds = eta(kd1) * (hos- hol(kd)) - alhl*(qos - qol(kd))
2150 dh = eta(kd1) * (hos- hol(kd))
2151
2152
2153 gms(kd) = (ds + st1 - tem1*det*alhl-tem*alhf) * pri(kd)
2154 gmh(kd) = pri(kd) * (hcc-eta(kd)*hos + dh)
2155!
2156! TENDENCY FOR SUSPENDED ENVIRONMENTAL ICE AND/OR LIQUID WATER
2157!
2158 qll(kd) = (tem2*(det-tem) + eta(kd1)*(qls-cll(kd)) &
2159 & + (one-qil(kd))*dlq(kd) - eta(kd)*qls ) * pri(kd)
2160
2161 qil(kd) = (tem2*tem + eta(kd1)*(qis-cil(kd)) &
2162 & + qil(kd)*dlq(kd) - eta(kd)*qis ) * pri(kd)
2163!
2164 ghd(kd) = zero
2165 gsd(kd) = zero
2166!
2167 DO l=kd1,k
2168 lm1 = l - 1
2169 st1 = one - alfint(l,1)
2170 st2 = one - alfint(l,2)
2171 st3 = one - alfint(l,3)
2172 st4 = one - alfint(l,4)
2173 st5 = one - alfind(l)
2174 hb = alfint(l,1)*hol(lm1) + st1*hol(l)
2175 qb = alfint(l,2)*qol(lm1) + st2*qol(l)
2176
2177 tem = alfint(l,4)*cil(lm1) + st4*cil(l)
2178 tem2 = alfint(l,3)*cll(lm1) + st3*cll(l)
2179
2180 tem1 = eta(l) * (tem - cil(l))
2181 tem3 = eta(l) * (tem2 - cll(l))
2182
2183 hbd = alfind(l)*hol(lm1) + st5*hol(l)
2184 qbd = alfind(l)*qol(lm1) + st5*qol(l)
2185
2186 tem5 = etd(l) * (hod(l) - hbd)
2187 tem6 = etd(l) * (qod(l) - qbd)
2188!
2189 dh = eta(l) * (hb - hol(l)) + tem5
2190 ds = dh - alhl * (eta(l) * (qb - qol(l)) + tem6)
2191
2192 gmh(l) = dh * pri(l)
2193 gms(l) = ds * pri(l)
2194!
2195 ghd(l) = tem5 * pri(l)
2196 gsd(l) = (tem5 - alhl * tem6) * pri(l)
2197!
2198 qll(l) = (tem3 + (one-qil(l))*dlq(l)) * pri(l)
2199 qil(l) = (tem1 + qil(l)*dlq(l)) * pri(l)
2200
2201 tem1 = eta(l) * (cil(lm1) - tem)
2202 tem3 = eta(l) * (cll(lm1) - tem2)
2203
2204 dh = eta(l) * (hol(lm1) - hb) - tem5
2205 ds = dh - alhl * eta(l) * (qol(lm1) - qb) &
2206 & + alhl * (tem6 - evp(lm1))
2207
2208 gmh(lm1) = gmh(lm1) + dh * pri(lm1)
2209 gms(lm1) = gms(lm1) + ds * pri(lm1)
2210!
2211 ghd(lm1) = ghd(lm1) - tem5 * pri(lm1)
2212 gsd(lm1) = gsd(lm1) - (tem5-alhl*(tem6-evp(lm1))) * pri(lm1)
2213
2214 qil(lm1) = qil(lm1) + tem1 * pri(lm1)
2215 qll(lm1) = qll(lm1) + tem3 * pri(lm1)
2216!
2217 avh = avh + gmh(lm1)*(prs(l)-prs(lm1))
2218
2219 ENDDO
2220!
2221 hbd = hol(k)
2222 qbd = qol(k)
2223 tem5 = etd(kp1) * (hod(kp1) - hbd)
2224 tem6 = etd(kp1) * (qod(kp1) - qbd)
2225 dh = - tem5
2226 ds = dh + alhl * tem6
2227 tem1 = dh * pri(k)
2228 tem2 = (ds - alhl * evp(k)) * pri(k)
2229 gmh(k) = gmh(k) + tem1
2230 gms(k) = gms(k) + tem2
2231 ghd(k) = ghd(k) + tem1
2232 gsd(k) = gsd(k) + tem2
2233
2234!
2235 avh = avh + gmh(k)*(prs(kp1)-prs(k))
2236!
2237 tem4 = - gravfac * pris
2238 tx1 = dh * tem4
2239 tx2 = ds * tem4
2240!
2241 DO l=kbl,k
2242 gmh(l) = gmh(l) + tx1
2243 gms(l) = gms(l) + tx2
2244 ghd(l) = ghd(l) + tx1
2245 gsd(l) = gsd(l) + tx2
2246!
2247 avh = avh + tx1*(prs(l+1)-prs(l))
2248 ENDDO
2249!
2250!***********************************************************************
2251!***********************************************************************
2252
2253!===> KERNEL (AKM) CALCULATION BEGINS
2254
2255!===> MODIFY SOUNDING WITH UNIT MASS FLUX
2256!
2257 DO l=kd,k
2258
2259 tem1 = gmh(l)
2260 tem2 = gms(l)
2261 hol(l) = hol(l) + tem1*testmb
2262 qol(l) = qol(l) + (tem1-tem2) * testmboalhl
2263 hst(l) = hst(l) + tem2*(one+gam(l))*testmb
2264 qst(l) = qst(l) + tem2*gam(l) * testmboalhl
2265 cll(l) = cll(l) + qll(l) * testmb
2266 cil(l) = cil(l) + qil(l) * testmb
2267 ENDDO
2268!
2269 if (alm > zero) then
2270 hos = hos + gmh(kd) * testmb
2271 qos = qos + (gmh(kd)-gms(kd)) * testmboalhl
2272 qls = qls + qll(kd) * testmb
2273 qis = qis + qil(kd) * testmb
2274 else
2275 st2 = one - st1s
2276 hos = hos + (st1s*gmh(kd)+st2*gmh(kd1)) * testmb
2277 qos = qos + (st1s * (gmh(kd)-gms(kd)) &
2278 & + st2 * (gmh(kd1)-gms(kd1))) * testmboalhl
2279 hst(kd) = hst(kd) + (st1s*gms(kd)*(one+gam(kd)) &
2280 & + st2*gms(kd1)*(one+gam(kd1))) * testmb
2281 qst(kd) = qst(kd) + (st1s*gms(kd)*gam(kd) &
2282 & + st2*gms(kd1)*gam(kd1)) * testmboalhl
2283
2284 qls = qls + (st1s*qll(kd)+st2*qll(kd1)) * testmb
2285 qis = qis + (st1s*qil(kd)+st2*qil(kd1)) * testmb
2286 endif
2287
2288!
2289 tem = prl(kmaxp1) - prl(kmax)
2290 hbl = hol(kmax) * tem
2291 qbl = qol(kmax) * tem
2292 qlb = cll(kmax) * tem
2293 qib = cil(kmax) * tem
2294 DO l=kmaxm1,kbl,-1
2295 tem = prl(l+1) - prl(l)
2296 hbl = hbl + hol(l) * tem
2297 qbl = qbl + qol(l) * tem
2298 qlb = qlb + cll(l) * tem
2299 qib = qib + cil(l) * tem
2300 ENDDO
2301 hbl = hbl * prism
2302 qbl = qbl * prism
2303 qlb = qlb * prism
2304 qib = qib * prism
2305
2306! if (ctei .and. sgcs(kd) > 0.65) then
2307! hbl = hbl * hpert_fac
2308! qbl = qbl * hpert_fac
2309! endif
2310
2311!***********************************************************************
2312
2313!===> CLOUD WORKFUNCTION FOR MODIFIED SOUNDING, THEN KERNEL (AKM)
2314!
2315 akm = zero
2316 tx1 = zero
2317 qtl = qst(kb1) - gaf(kb1)*hst(kb1)
2318 qtv = qbl
2319 hcc = hbl
2320 tx2 = hcc
2321 tx4 = (alhf*half)*max(zero,min(one,(tcr-tcl-tol(kb1))*tcrf))
2322!
2323 qtv = qbl
2324 tx1 = qib + qlb
2325!
2326
2327 DO l=kb1,kd1,-1
2328 lm1 = l - 1
2329 lp1 = l + 1
2330 del_eta = eta(l) - eta(lp1)
2331 hccp = hcc + del_eta*hol(l)
2332!
2333 qtlp = qst(lm1) - gaf(lm1)*hst(lm1)
2334 qtvp = half * ((qtlp+qtl)*eta(l) + (gaf(l)+gaf(lm1))*hccp)
2335
2336 detp = (bkc(l)*tx1 - (qtvp-qtv) &
2337 & + del_eta*(qol(l)+cll(l)+cil(l)) &
2338 & + eta(l)*q0u(l) + eta(lp1)*q0d(l)) * akc(l)
2339 IF (detp <= zero) cnvflg = .true.
2340
2341 st1 = hst(l) - ltl(l)*nu*(qst(l)-qol(l))
2342
2343 tem2 = (alhf*half)*max(zero,min(one,(tcr-tcl-tol(lm1))*tcrf))
2344 tem1 = hccp + detp * (tem2+tx4)
2345
2346 st2 = ltl(l) * vtf(l)
2347 tem5 = cll(l) + cil(l)
2348 akm = akm + &
2349 & ( (tx2 -eta(lp1)*st1-st2*(tx1-tem5*eta(lp1))) * dlb(l) &
2350 & + (tem1 -eta(l )*st1-st2*(detp-tem5*eta(l))) * dlt(l) )
2351!
2352 hcc = hccp
2353 tx1 = detp
2354 tx2 = tem1
2355 qtl = qtlp
2356 qtv = qtvp
2357 tx4 = tem2
2358 ENDDO
2359!
2360 if (cnvflg) return
2361!
2362! Eventhough we ignore the change in lambda, we still assume
2363! that the cLoud-top contribution is zero; as though we still
2364! had non-bouyancy there.
2365!
2366!
2367 st1 = hst(kd) - ltl(kd)*nu*(qst(kd)-qos)
2368 st2 = ltl(kd) * vtf(kd)
2369 tem5 = (qls + qis) * eta(kd1)
2370 akm = akm + half * (tx2-eta(kd1)*st1-st2*(tx1-tem5)) * dlb(kd)
2371!
2372 akm = (akm - wfn) * testmbi
2373
2374
2375!***********************************************************************
2376
2377!===> MASS FLUX
2378!
2379 amb = - (wfn-acr) / akm
2380!
2381!===> RELAXATION AND CLIPPING FACTORS
2382!
2383 amb = amb * clp * rel_fac
2384
2385!!! if (DDFT) AMB = MIN(AMB, ONE/CLDFRD)
2386
2387!===> SUB-CLOUD LAYER DEPTH LIMIT ON MASS FLUX
2388
2389 ambmax = (prl(kmaxp1)-prl(kbl))*(fracbl*gravcon)
2390 amb = max(min(amb, ambmax),zero)
2391
2392!***********************************************************************
2393!*************************RESULTS***************************************
2394!***********************************************************************
2395
2396!===> PRECIPITATION AND CLW DETRAINMENT
2397!
2398 if (amb > zero) then
2399
2400!
2401! if (wvl(kd) > zero) then
2402! tx1 = one - amb * eta(kd) / (rho(kd)*wvl(kd))
2403! sigf(kd) = max(zero, min(one, tx1 * tx1))
2404! endif
2405 if (aw_scal) then
2406 tx1 = (0.2_kp / max(alm, 1.0e-5_kp))
2407 tx2 = one - min(one, pi * tx1 * tx1 / area)
2408
2409 tx2 = tx2 * tx2
2410
2411! comnet out the following for now - 07/23/18
2412! do l=kd1,kbl
2413! lp1 = min(K, l+1)
2414! if (wvl(l) > zero .and. wvl(lp1) > zero) then
2415! tx1 = one - amb * (eta(l)+eta(lp1))
2416! & / ((wvl(l)+wvl(lp1))*rho(l)*grav)
2417! sigf(l) = max(zero, min(one, tx1 * tx1))
2418! else
2419! sigf(l) = min(one,tx2)
2420! endif
2421! sigf(l) = max(sigf(l), tx2)
2422! enddo
2423! sigf(kd) = sigf(kd1)
2424! if (kbl < k) then
2425! sigf(kbl+1:k) = sigf(kbl)
2426! endif
2427 sigf(kd:k) = tx2
2428 else
2429 sigf(kd:k) = one
2430 endif
2431
2432 tx1 = max(1.0e-6_kp, abs(gms(kd) * onebcp * sigf(kd)))
2433 amb = min(tx1*amb, tfrac_max*toi(kd)) / tx1
2434
2435!
2436 avt = zero
2437 avq = zero
2438 avr = dof * sigf(kbl)
2439!
2440 dsfc = dsfc + amb * etd(k) * (one/dt) * sigf(kbl)
2441!
2442 DO l=k,kd,-1
2443 pcu(l) = pcu(l) + amb*rnn(l)*sigf(l) ! (A40)
2444 avr = avr + rnn(l) * sigf(l)
2445 ENDDO
2446 pcu(k) = pcu(k) + amb * dof * sigf(kbl)
2447!
2448!===> TEMPARATURE AND Q CHANGE AND CLOUD MASS FLUX DUE TO CLOUD TYPE KD
2449!
2450 tx1 = amb * onebcp
2451 tx2 = amb * oneoalhl
2452 DO l=kd,k
2453 delp = prs(l+1) - prs(l)
2454 tx3 = amb * sigf(l)
2455 st1 = gms(l) * tx1 * sigf(l)
2456 toi(l) = toi(l) + st1
2457 tcu(l) = tcu(l) + st1
2458 tcd(l) = tcd(l) + gsd(l) * tx1 * sigf(l)
2459!
2460 st1 = st1 - elocp * (qil(l) + qll(l)) * tx3
2461
2462 avt = avt + st1 * delp
2463
2464 flx(l) = flx(l) + eta(l) * tx3
2465 flxd(l) = flxd(l) + etd(l) * tx3
2466!
2467 qii(l) = qii(l) + qil(l) * tx3
2468 tem = zero
2469
2470 qli(l) = qli(l) + qll(l) * tx3 + tem
2471
2472 st1 = (gmh(l)-gms(l)) * tx2 * sigf(l)
2473
2474 qoi(l) = qoi(l) + st1
2475 qcu(l) = qcu(l) + st1
2476 qcd(l) = qcd(l) + (ghd(l)-gsd(l)) * tx2 * sigf(l)
2477!
2478 avq = avq + (st1 + (qll(l)+qil(l))*tx3) * delp
2479! avq = avq + st1 * (prs(l+1)-prs(l))
2480! avr = avr + (QLL(L) + QIL(L)*(1+alhf/alhl))
2481 avr = avr + (qll(l) + qil(l)) * delp * sigf(l) * gravcon
2482
2483! Correction for negative condensate!
2484 if (qii(l) < zero) then
2485 tem = qii(l) * elfocp
2486 qoi(l) = qoi(l) + qii(l)
2487 qcu(l) = qcu(l) + qii(l)
2488 toi(l) = toi(l) - tem
2489 tcu(l) = tcu(l) - tem
2490 qii(l) = zero
2491 endif
2492 if (qli(l) < zero) then
2493 tem = qli(l) * elocp
2494 qoi(l) = qoi(l) + qli(l)
2495 qcu(l) = qcu(l) + qli(l)
2496 toi(l) = toi(l) - tem
2497 tcu(l) = tcu(l) - tem
2498 qli(l) = zero
2499 endif
2500
2501 ENDDO
2502 avr = avr * amb
2503!
2504! Correction for negative condensate!
2505! if (advcld) then
2506! do l=kd,k
2507! if (qli(l) < zero) then
2508! qoi(l) = qoi(l) + qli(l)
2509! toi(l) = toi(l) - (alhl/cp) * qli(l)
2510! qli(l) = zero
2511! endif
2512! if (qii(l) < zero) then
2513! qoi(l) = qoi(l) + qii(l)
2514! toi(l) = toi(l) - ((alhl+alhf)/cp) * qii(l)
2515! qii(l) = zero
2516! endif
2517! enddo
2518! endif
2519
2520!
2521!
2522 tx1 = zero
2523 tx2 = zero
2524!
2525 IF (revap) THEN ! REEVAPORATION OF FALLING CONVECTIVE RAIN
2526!
2527 tem = zero
2528 do l=kd,kbl
2529 IF (l < idh .or. (.not. ddft)) THEN
2530 tem = tem + amb * rnn(l) * sigf(l)
2531 endif
2532 enddo
2533 tem = tem + amb * dof * sigf(kbl)
2534 tem = tem * (3600.0_kp/dt)
2535 tem1 = sqrt(max(one, min(100.0_kp,(6.25e10_kp/max(area,one))))) ! 20110530
2536
2537 clfrac = max(zero, min(half, rknob*clf(tem)*tem1))
2538 cldfrd = clfrac
2539!
2540 DO l=kd,kbl ! Testing on 20070926
2541! for L=KD,K
2542 IF (l >= idh .AND. ddft) THEN
2543 tem = amb * sigf(l)
2544 tx2 = tx2 + tem * rnn(l)
2545 cldfrd = min(tem*cldfr(l), clfrac)
2546 ELSE
2547 tx1 = tx1 + amb * rnn(l) * sigf(l)
2548 ENDIF
2549 tx4 = zfac * phil(l)
2550 tx4 = (one - tx4 * (one - half*tx4)) * afc
2551!
2552 IF (tx1 > zero .OR. tx2 > zero) THEN
2553 teq = toi(l)
2554 qeq = qoi(l)
2555 pl = half * (prl(l+1)+prl(l))
2556
2557 st1 = max(zero, min(one, (tcr-teq)*tcrf))
2558 st2 = st1*elfocp + (one-st1)*elocp
2559
2560 CALL qsatcn ( teq,pl,qsteq,dqdt)
2561!
2562 deltaq = half * (qsteq*rhc_ls(l)-qeq) / (one+st2*dqdt)
2563!
2564 qeq = qeq + deltaq
2565 teq = teq - deltaq*st2
2566!
2567 tem1 = max(zero, min(one, (tcr-teq)*tcrf))
2568 tem2 = tem1*elfocp + (one-tem1)*elocp
2569
2570 CALL qsatcn ( teq,pl,qsteq,dqdt)
2571!
2572 deltaq = (qsteq*rhc_ls(l)-qeq) / (one+tem2*dqdt)
2573!
2574 qeq = qeq + deltaq
2575 teq = teq - deltaq*tem2
2576
2577 IF (qeq > qoi(l)) THEN
2578 potevap = (qeq-qoi(l))*(prl(l+1)-prl(l))*gravcon
2579
2580 tem4 = zero
2581 if (tx1 > zero) &
2582 & tem4 = potevap * (one - exp( tx4*tx1**0.57777778_kp ))
2583 actevap = min(tx1, tem4*clfrac)
2584
2585 if (tx1 < rainmin*dt) actevap = min(tx1, potevap)
2586!
2587 tem4 = zero
2588 if (tx2 > zero) &
2589 & tem4 = potevap * (one - exp( tx4*tx2**0.57777778_kp ))
2590 tem4 = min(min(tx2, tem4*cldfrd), potevap-actevap)
2591 if (tx2 < rainmin*dt) tem4 = min(tx2, potevap-actevap)
2592!
2593 tx1 = tx1 - actevap
2594 tx2 = tx2 - tem4
2595 st1 = (actevap+tem4) * pri(l)
2596 qoi(l) = qoi(l) + st1
2597 qcu(l) = qcu(l) + st1
2598!
2599
2600 st1 = st1 * elocp
2601 toi(l) = toi(l) - st1
2602 tcu(l) = tcu(l) - st1
2603 ENDIF
2604 ENDIF
2605 ENDDO
2606!
2607 cup = cup + tx1 + tx2 + dof * amb * sigf(kbl)
2608 ELSE
2609 DO l=kd,k
2610 tx1 = tx1 + amb * rnn(l) * sigf(l)
2611 ENDDO
2612 cup = cup + tx1 + dof * amb * sigf(kbl)
2613 ENDIF
2614!
2615! Convective transport (mixing) of passive tracers
2616!
2617 if (ntrc > 0) then
2618 do l=kd,km1
2619 if (etz(l) /= zero) etzi(l) = one / etz(l)
2620 enddo
2621 DO n=1,ntrc ! Tracer loop ; first two are u and v
2622
2623 DO l=kd,k
2624 hol(l) = roi(l,n)
2625 ENDDO
2626!
2627 hcc = rbl(n)
2628 hod(kd) = hol(kd)
2629! Compute downdraft properties for the tracer
2630 DO l=kd1,k
2631 lm1 = l - 1
2632 st1 = one - alfind(l)
2633 hb = alfind(l) * hol(lm1) + st1 * hol(l)
2634 IF (etz(lm1) /= zero) THEN
2635 tem = etzi(lm1)
2636 IF (etd(l) > etd(lm1)) THEN
2637 hod(l) = (etd(lm1)*(hod(lm1)-hol(lm1)) &
2638 & + etd(l) *(hol(lm1)-hb) + etz(lm1)*hb) * tem
2639 ELSE
2640 hod(l) = (etd(lm1)*(hod(lm1)-hb) + etz(lm1)*hb) * tem
2641 ENDIF
2642 ELSE
2643 hod(l) = hb
2644 ENDIF
2645 ENDDO
2646
2647 DO l=kb1,kd,-1
2648 hcc = hcc + (eta(l)-eta(l+1))*hol(l)
2649 ENDDO
2650!
2651! Scavenging -- fscav - fraction scavenged [km-1]
2652! delz - distance from the entrainment to detrainment layer [km]
2653! fnoscav - the fraction not scavenged
2654! following Liu et al. [JGR,2001] Eq 1
2655
2656 if (fscav_(n) > zero) then
2657 delzkm = ( phil(kd) - phih(kd1) ) *(onebg*0.001_kp)
2658 fnoscav = exp(- fscav_(n) * delzkm)
2659 else
2660 fnoscav = one
2661 endif
2662
2663 gmh(kd) = pri(kd) * (hcc-eta(kd)*hol(kd)) * trcfac(kd,n) &
2664 & * fnoscav
2665 DO l=kd1,k
2666 if (fscav_(n) > zero) then
2667 delzkm = ( phil(kd) - phih(l+1) ) *(onebg*0.001_kp)
2668 fnoscav = exp(- fscav_(n) * delzkm)
2669 endif
2670 lm1 = l - 1
2671 st1 = one - alfint(l,n+4)
2672 st2 = one - alfind(l)
2673 hb = alfint(l,n+4) * hol(lm1) + st1 * hol(l)
2674 hbd = alfind(l) * hol(lm1) + st2 * hol(l)
2675 tem5 = etd(l) * (hod(l) - hbd)
2676 dh = eta(l) * (hb - hol(l)) * fnoscav + tem5
2677 gmh(l ) = dh * pri(l) * trcfac(l,n)
2678 dh = eta(l) * (hol(lm1) - hb) * fnoscav - tem5
2679 gmh(lm1) = gmh(lm1) + dh * pri(lm1) * trcfac(l,n)
2680 ENDDO
2681!
2682 st2 = zero
2683 DO l=kd,k
2684 st1 = gmh(l)*amb*sigf(l) + st2
2685 st3 = hol(l) + st1
2686 st2 = st3 - trcmin(n) ! if trcmin is defined limit change
2687 if (st2 < zero) then
2688 roi(l,n) = trcmin(n)
2689 rcu(l,n) = rcu(l,n) + st1
2690 if (l < k) &
2691 & st2 = st2 * (prl(l+1)-prl(l))*pri(l+1) * (cmb2pa/grav)
2692 else
2693 roi(l,n) = st3
2694 rcu(l,n) = rcu(l,n) + st1
2695 st2 = zero
2696 endif
2697
2698 ENDDO
2699 ENDDO ! Tracer loop NTRC
2700 endif
2701 endif ! amb > zero
2702
2703 RETURN
2704 end subroutine cloud
2705
2707 SUBROUTINE ddrft( &
2708 & K, KP1, KD &
2709 &, TLA, ALFIND, wcbase &
2710 &, TOL, QOL, HOL, PRL, QST, HST, GAM, GAF &
2711! &, TOL, QOL, HOL, PRL, QST, HST, GAM, GAF, HBL, QBL&
2712 &, QRB, QRT, BUY, KBL, IDH, ETA, RNN, ETAI &
2713 &, ALM, WFN, TRAIN, DDFT &
2714 &, ETD, HOD, QOD, EVP, DOF, CLDFRD, WCB &
2715 &, GMS, GSD, GHD, wvlu)
2716
2717!
2718!***********************************************************************
2719!******************** Cumulus Downdraft Subroutine *********************
2720!****************** Based on Cheng and Arakawa (1997) ****** **********
2721!************************ SUBROUTINE DDRFT ****************************
2722!************************* October 2004 ******************************
2723!***********************************************************************
2724!***********************************************************************
2725!************* Shrinivas.Moorthi@noaa.gov (301) 683-3718 ***************
2726!***********************************************************************
2727!***********************************************************************
2728!23456789012345678901234567890123456789012345678901234567890123456789012
2729!
2730!===> TOL(K) INPUT TEMPERATURE KELVIN
2731!===> QOL(K) INPUT SPECIFIC HUMIDITY NON-DIMENSIONAL
2732
2733!===> PRL(KP1) INPUT PRESSURE @ EDGES MB
2734
2735!===> K INPUT THE RISE & THE INDEX OF THE SUBCLOUD LAYER
2736!===> KD INPUT DETRAINMENT LEVEL ( 1<= KD < K )
2737!
2738 IMPLICIT NONE
2739!
2740! INPUT ARGUMENTS
2741!
2742 INTEGER K, KP1, KD, KBL
2743 real(kind=kind_phys) alfind(k), wcbase
2744
2745 real(kind=kind_phys), dimension(kd:k) :: hol, qol, hst, qst &
2746 &, tol, qrb, qrt, rnn &
2747 &, rns, etai
2748 real(kind=kind_phys), dimension(kd:kp1) :: gaf, buy, gam, eta &
2749 &, prl
2750!
2751! real(kind=kind_phys) HBL, QBL, PRIS &
2752! &, TRAIN, WFN, ALM
2753!
2754! TEMPORARY WORK SPACE
2755!
2756 real(kind=kind_phys), dimension(KD:K) :: rnf, wcb, evp, stlt &
2757 &, ghd, gsd, cldfrd &
2758 &, gqw, qrpi, qrps, bud
2759
2760 real(kind=kind_phys), dimension(KD:KP1) :: qrp, wvl, wvlu, etd &
2761 &, hod, qod, ror, gms
2762
2763 real(kind=kind_phys) tl, pl, ql, qs, dqs, st1 &
2764 &, qqq, del_eta, hb, qb, tb &
2765 &, tem, tem1, tem2, tem3, tem4, st2 &
2766 &, errmin, errmi2, errh, errw, erre, tem5 &
2767 &, tem6, hbd, qbd, tx1, tx2, tx3 &
2768 &, tx4, tx5, tx6, tx7, tx8, tx9 &
2769 &, wfn, alm, al2 &
2770 &, train, gmf, onpg, ctla, vtrm &
2771 &, rpart, qrmin, aa1, bb1, cc1, dd1 &
2772! &, WC2MIN, WCMIN, WCBASE, F2, F3, F5 &
2773 &, wc2min, wcmin, f2, f3, f5 &
2774 &, gmf1, gmf5, qraf, qrbf, del_tla &
2775 &, tla, stla, ctl2, ctl3 &
2776! &, TLA, STLA, CTL2, CTL3, ASIN &
2777! &, RNT, RNB, ERRQ, RNTP, QRPF, VTPF &
2778 &, rnt, rnb, errq, rntp &
2779 &, edz, ddz, ce, qhs, fac, facg &
2780 &, rsum1, rsum2, rsum3, cee, dof, dofw
2781! &, sialf
2782
2783 INTEGER I, L, N, IX, KD1, II, kb1, IP1, JJ, ntla &
2784 &, IT, KM1, KTEM, KK, KK1, LM1, LL, LP1 &
2785 &, IDW, IDH, IDN(K), idnm, itr
2786!
2787 parameter(errmin=0.0001_kp, errmi2=0.1_kp*errmin)
2788! parameter (ERRMIN=0.00001, ERRMI2=0.1*ERRMIN)
2789!
2790! real (kind=kind_phys), parameter :: PIINV=one/PI, pio2=half*pi
2791!
2792 parameter(onpg=one+half, gmf=one/onpg, rpart=zero)
2793! parameter (ONPG=1.0+0.5, GMF=1.0/ONPG, RPART=1.0)
2794! parameter (ONPG=1.0+0.5, GMF=1.0/ONPG, RPART=0.5)
2795! PARAMETER (AA1=1.0, BB1=1.5, CC1=1.1, DD1=0.85, F3=CC1, F5=2.5)
2796! PARAMETER (AA1=2.0, BB1=1.5, CC1=1.1, DD1=0.85, F3=CC1, F5=2.5)
2797 parameter(aa1=1.0_kp, bb1=1.0_kp, cc1=1.0_kp, dd1=1.0_kp, &
2798 & f3=cc1, f5=1.0_kp)
2799 parameter(qrmin=1.0e-6_kp, wc2min=0.01_kp, gmf1=gmf/aa1, gmf5=gmf/f5)
2800! parameter (QRMIN=1.0E-6, WC2MIN=1.00, GMF1=GMF/AA1, GMF5=GMF/F5)
2801 parameter(wcmin=sqrt(wc2min))
2802! parameter (sialf=0.5)
2803!
2804 integer, parameter :: itrmu=25, itrmd=25 &
2805 &, itrmin=15, itrmnd=12, numtla=2
2806
2807! uncentering for vvel in dd
2808 real(kind=kind_phys), parameter :: ddunc1=0.25_kp &
2809 &, ddunc2=one-ddunc1 &
2810! &, ddunc1=0.4, ddunc2=one-ddunc1 &
2811! &, ddunc1=0.3, ddunc2=one-ddunc1 &
2812 &, vtpexp=-0.3636_kp &
2813 &, vtp=36.34_kp*sqrt(1.2_kp)*(0.001_kp)**0.1364_kp
2814!
2815! real(kind=kind_phys) EM(K*K), ELM(K)
2816 real(kind=kind_phys) elm(k), aa(kd:k,kd:kp1), qw(kd:k,kd:k) &
2817 &, vt(2), vrw(2), trw(2), qa(3), wa(3)
2818
2819 LOGICAL SKPUP, cnvflg, DDFT, UPDRET, DDLGK
2820
2821!***********************************************************************
2822
2823
2824 kd1 = kd + 1
2825 km1 = k - 1
2826 kb1 = kbl - 1
2827!
2828! VTP = 36.34*SQRT(1.2)* (0.001)**0.1364
2829! VTPEXP = -0.3636
2830! PIINV = 1.0 / PI
2831! PICON = PIO2 * ONEBG
2832!
2833! Compute Rain Water Budget of the Updraft (Cheng and Arakawa, 1997)
2834!
2835 cldfrd = zero
2836 rntp = zero
2837 dof = zero
2838 errq = 10.0_kp
2839 rnb = zero
2840 rnt = zero
2841 tx2 = prl(kbl)
2842!
2843 tx1 = (prl(kd) + prl(kd1)) * half
2844 ror(kd) = cmpor*tx1 / (tol(kd)*(one+nu*qol(kd)))
2845! GMS(KD) = VTP * ROR(KD) ** VTPEXP
2846 gms(kd) = vtp * vtpf(ror(kd))
2847!
2848 qrp(kd) = qrmin
2849!
2850 tem = tol(k) * (one + nu * qol(k))
2851 ror(kp1) = half * cmpor * (prl(kp1)+prl(k)) / tem
2852 gms(kp1) = vtp * vtpf(ror(kp1))
2853 qrp(kp1) = qrmin
2854!
2855 kk = kbl
2856 DO l=kd1,k
2857 tem = half * (tol(l)+tol(l-1)) &
2858 & * (one + (half*nu) * (qol(l)+qol(l-1)))
2859 ror(l) = cmpor * prl(l) / tem
2860! GMS(L) = VTP * ROR(L) ** VTPEXP
2861 gms(l) = vtp * vtpf(ror(l))
2862 qrp(l) = qrmin
2863 if (buy(l) <= zero .and. kk == kbl) then
2864 kk = l
2865 endif
2866 ENDDO
2867 if (kk /= kbl) then
2868 do l=kk,kbl
2869 buy(l) = 0.9_kp * buy(l-1)
2870 enddo
2871 endif
2872!
2873 do l=kd,k
2874 qrpi(l) = buy(l)
2875 enddo
2876 do l=kd1,kb1
2877 buy(l) = 0.25_kp * (qrpi(l-1)+qrpi(l)+qrpi(l)+qrpi(l+1))
2878 enddo
2879
2880!
2881! CALL ANGRAD(TX1, ALM, STLA, CTL2, AL2, PI, TLA, TX2, WFN, TX3)
2882 tx1 = 1000.0_kp + tx1 - prl(kp1)
2883! CALL ANGRAD(TX1, ALM, AL2, TLA, TX2, WFN, TX3)
2884 CALL angrad(tx1, alm, al2, tla)
2885!
2886! Following Ucla approach for rain profile
2887!
2888 f2 = (bb1+bb1)*onebg/(pi*0.2_kp)
2889! WCMIN = SQRT(WC2MIN)
2890! WCBASE = WCMIN
2891!
2892! del_tla = TLA * 0.2
2893! del_tla = TLA * 0.25
2894 del_tla = tla * 0.3_kp
2895 tla = tla - del_tla
2896!
2897 DO l=kd,k
2898 rnf(l) = zero
2899 rns(l) = zero
2900 stlt(l) = zero
2901 gqw(l) = zero
2902 qrp(l) = qrmin
2903 DO n=kd,k
2904 qw(n,l) = zero
2905 ENDDO
2906 ENDDO
2907! DO L=KD,KP1
2908! WVL(L) = zero
2909! ENDDO
2910!
2911!-----QW(N,L) = D(W(N)*W(N))/DQR(L)
2912!
2913 kk = kbl
2914 qw(kd,kd) = -qrb(kd) * gmf1
2915 ghd(kd) = eta(kd) * eta(kd)
2916 gqw(kd) = qw(kd,kd) * ghd(kd)
2917 gsd(kd) = etai(kd) * etai(kd)
2918!
2919 gqw(kk) = - qrb(kk-1) * (gmf1+gmf1)
2920!
2921 wcb(kk) = wcbase * wcbase
2922
2923 tx1 = wcb(kk)
2924 gsd(kk) = one
2925 ghd(kk) = one
2926!
2927 tem = gmf1 + gmf1
2928 DO l=kb1,kd1,-1
2929 ghd(l) = eta(l) * eta(l)
2930 gsd(l) = etai(l) * etai(l)
2931 gqw(l) = - ghd(l) * (qrb(l-1)+qrt(l)) * tem
2932 qw(l,l) = - qrt(l) * tem
2933!
2934 st1 = half * (eta(l) + eta(l+1))
2935 tx1 = tx1 + buy(l) * tem * (qrb(l)+qrt(l)) * st1 * st1
2936 wcb(l) = tx1 * gsd(l)
2937 ENDDO
2938!
2939 tem1 = (qrb(kd) + qrt(kd1) + qrt(kd1)) * gmf1
2940 gqw(kd1) = - ghd(kd1) * tem1
2941 qw(kd1,kd1) = - qrt(kd1) * tem
2942 st1 = half * (eta(kd) + eta(kd1))
2943 wcb(kd) = (tx1 + buy(kd)*tem*qrb(kd)*st1*st1) * gsd(kd)
2944!
2945 DO l=kd1,kbl
2946 DO n=kd,l-1
2947 qw(n,l) = gqw(l) * gsd(n)
2948 ENDDO
2949 ENDDO
2950 qw(kbl,kbl) = zero
2951!
2952 do ntla=1,numtla ! numtla is the the maximimu number of tilting angle tries
2953 ! ------
2954! if (errq < 1.0 .or. tla > 45.0) cycle
2955 if (errq < 0.1_kp .or. tla > 45.0_kp) cycle
2956!
2957 tla = tla + del_tla
2958 stla = sin(tla*deg2rad) ! sine of tilting angle
2959 ctl2 = one - stla * stla ! cosine square of tilting angle
2960!
2961 stla = f2 * stla * al2
2962 ctl2 = dd1 * ctl2
2963 ctl3 = 0.1364_kp * ctl2
2964!
2965 DO l=kd,k
2966 rnf(l) = zero
2967 stlt(l) = zero
2968 qrp(l) = qrmin
2969 ENDDO
2970 DO l=kd,kp1
2971 wvl(l) = zero
2972 ENDDO
2973 wvl(kbl) = wcbase
2974 stlt(kbl) = one / wcbase
2975!
2976 DO l=kd,kp1
2977 DO n=kd,k
2978 aa(n,l) = zero
2979 ENDDO
2980 ENDDO
2981!
2982 skpup = .false.
2983!
2984 DO itr=1,itrmu ! Rain Profile Iteration starts!
2985 IF (.NOT. skpup) THEN
2986! wvlu = wvl
2987!
2988!-----CALCULATING THE VERTICAL VELOCITY
2989!
2990 tx1 = zero
2991 qrpi(kbl) = one / qrp(kbl)
2992 DO l=kb1,kd,-1
2993 tx1 = tx1 + qrp(l+1)*gqw(l+1)
2994 st1 = wcb(l) + qw(l,l)*qrp(l) + tx1*gsd(l)
2995! if (st1 > wc2min) then
2996 if (st1 > zero) then
2997 wvl(l) = max(ddunc1*sqrt(st1) + ddunc2*wvl(l), wcmin)
2998! WVL(L) = SQRT(ST1)
2999! WVL(L) = max(half * (SQRT(ST1) + WVL(L)), wcmin)
3000! qrp(l) = half*((wvl(l)*wvl(l)-wcb(l)-tx1*gsd(l))/qw(l,l)&
3001! & + qrp(l))
3002 else
3003
3004! wvl(l) = 0.5*(wcmin+wvl(l))
3005! wvl(l) = max(half*(wvl(l) + wvl(l+1)), wcmin)
3006 wvl(l) = max(wvl(l),wcmin)
3007 qrp(l) = (wvl(l)*wvl(l) - wcb(l) - tx1*gsd(l))/qw(l,l)
3008! qrp(l) = half*((wvl(l)*wvl(l)-wcb(l)-tx1*gsd(l))/qw(l,l)&
3009! & + qrp(l))
3010 endif
3011 qrp(l) = max(qrp(l), qrmin)
3012
3013 stlt(l) = one / wvl(l)
3014 qrpi(l) = one / qrp(l)
3015 ENDDO
3016!
3017!-----CALCULATING TRW, VRW AND OF
3018!
3019! VT(1) = GMS(KD) * QRP(KD)**0.1364
3020 vt(1) = gms(kd) * qrpf(qrp(kd))
3021 trw(1) = eta(kd) * qrp(kd) * stlt(kd)
3022 tx6 = trw(1) * vt(1)
3023 vrw(1) = f3*wvl(kd) - ctl2*vt(1)
3024 bud(kd) = stla * tx6 * qrb(kd) * half
3025 rnf(kd) = bud(kd)
3026 dof = 1.1364_kp * bud(kd) * qrpi(kd)
3027 dofw = -bud(kd) * stlt(kd)
3028!
3029 rnt = trw(1) * vrw(1)
3030 tx2 = zero
3031 tx4 = zero
3032 rnb = rnt
3033 tx1 = half
3034 tx8 = zero
3035!
3036 IF (rnt >= zero) THEN
3037 tx3 = (rnt-ctl3*tx6) * qrpi(kd)
3038 tx5 = ctl2 * tx6 * stlt(kd)
3039 ELSE
3040 tx3 = zero
3041 tx5 = zero
3042 rnt = zero
3043 rnb = zero
3044 ENDIF
3045!
3046 DO l=kd1,kb1
3047 ktem = max(l-2, kd)
3048 ll = l - 1
3049!
3050! VT(2) = GMS(L) * QRP(L)**0.1364
3051 vt(2) = gms(l) * qrpf(qrp(l))
3052 trw(2) = eta(l) * qrp(l) * stlt(l)
3053 vrw(2) = f3*wvl(l) - ctl2*vt(2)
3054 qqq = stla * trw(2) * vt(2)
3055 st1 = tx1 * qrb(ll)
3056 bud(l) = qqq * (st1 + qrt(l))
3057!
3058 qa(2) = dof
3059 wa(2) = dofw
3060 dof = 1.1364_kp * bud(l) * qrpi(l)
3061 dofw = -bud(l) * stlt(l)
3062!
3063 rnf(ll) = rnf(ll) + qqq * st1
3064 rnf(l) = qqq * qrt(l)
3065!
3066 tem3 = vrw(1) + vrw(2)
3067 tem4 = trw(1) + trw(2)
3068!
3069 tx6 = pt25 * tem3 * tem4
3070 tem4 = tem4 * ctl3
3071!
3072!-----BY QR ABOVE
3073!
3074! TEM1 = pt25*(TRW(1)*TEM3 - TEM4*VT(1))*TX7
3075 tem1 = pt25*(trw(1)*tem3 - tem4*vt(1))*qrpi(ll)
3076 st1 = pt25*(trw(1)*(ctl2*vt(1)-vrw(2)) &
3077 & * stlt(ll) + f3*trw(2))
3078!-----BY QR BELOW
3079 tem2 = pt25*(trw(2)*tem3 - tem4*vt(2))*qrpi(l)
3080 st2 = pt25*(trw(2)*(ctl2*vt(2)-vrw(1)) &
3081 & * stlt(l) + f3*trw(1))
3082!
3083! From top to the KBL-2 layer
3084!
3085 qa(1) = tx2
3086 qa(2) = qa(2) + tx3 - tem1
3087 qa(3) = -tem2
3088!
3089 wa(1) = tx4
3090 wa(2) = wa(2) + tx5 - st1
3091 wa(3) = -st2
3092!
3093 tx2 = tem1
3094 tx3 = tem2
3095 tx4 = st1
3096 tx5 = st2
3097!
3098 vt(1) = vt(2)
3099 trw(1) = trw(2)
3100 vrw(1) = vrw(2)
3101!
3102 IF (wvl(ktem) == wcmin) wa(1) = zero
3103 IF (wvl(ll) == wcmin) wa(2) = zero
3104 IF (wvl(l) == wcmin) wa(3) = zero
3105 DO n=ktem,kbl
3106 aa(ll,n) = (wa(1)*qw(ktem,n) * stlt(ktem) &
3107 & + wa(2)*qw(ll,n) * stlt(ll) &
3108 & + wa(3)*qw(l,n) * stlt(l) ) * half
3109 ENDDO
3110 aa(ll,ktem) = aa(ll,ktem) + qa(1)
3111 aa(ll,ll) = aa(ll,ll) + qa(2)
3112 aa(ll,l) = aa(ll,l) + qa(3)
3113 bud(ll) = (tx8 + rnn(ll)) * half &
3114 & - rnb + tx6 - bud(ll)
3115 aa(ll,kbl+1) = bud(ll)
3116 rnb = tx6
3117 tx1 = one
3118 tx8 = rnn(ll)
3119 ENDDO
3120 l = kbl
3121 ll = l - 1
3122! VT(2) = GMS(L) * QRP(L)**0.1364
3123 vt(2) = gms(l) * qrpf(qrp(l))
3124 trw(2) = eta(l) * qrp(l) * stlt(l)
3125 vrw(2) = f3*wvl(l) - ctl2*vt(2)
3126 st1 = stla * trw(2) * vt(2) * qrb(ll)
3127 bud(l) = st1
3128
3129 qa(2) = dof
3130 wa(2) = dofw
3131 dof = 1.1364_kp * bud(l) * qrpi(l)
3132 dofw = -bud(l) * stlt(l)
3133!
3134 rnf(ll) = rnf(ll) + st1
3135!
3136 tem3 = vrw(1) + vrw(2)
3137 tem4 = trw(1) + trw(2)
3138!
3139 tx6 = pt25 * tem3 * tem4
3140 tem4 = tem4 * ctl3
3141!
3142!-----BY QR ABOVE
3143!
3144 tem1 = pt25*(trw(1)*tem3 - tem4*vt(1))*qrpi(ll)
3145 st1 = pt25*(trw(1)*(ctl2*vt(1)-vrw(2)) &
3146 & * stlt(ll) + f3*trw(2))
3147!-----BY QR BELOW
3148 tem2 = pt25*(trw(2)*tem3 - tem4*vt(2))*qrpi(l)
3149 st2 = pt25*(trw(2)*(ctl2*vt(2)-vrw(1)) &
3150 & * stlt(l) + f3*trw(1))
3151!
3152! For the layer next to the top of the boundary layer
3153!
3154 qa(1) = tx2
3155 qa(2) = qa(2) + tx3 - tem1
3156 qa(3) = -tem2
3157!
3158 wa(1) = tx4
3159 wa(2) = wa(2) + tx5 - st1
3160 wa(3) = -st2
3161!
3162 tx2 = tem1
3163 tx3 = tem2
3164 tx4 = st1
3165 tx5 = st2
3166!
3167 idw = max(l-2, kd)
3168!
3169 IF (wvl(idw) == wcmin) wa(1) = zero
3170 IF (wvl(ll) == wcmin) wa(2) = zero
3171 IF (wvl(l) == wcmin) wa(3) = zero
3172!
3173 kk = idw
3174 DO n=kk,l
3175 aa(ll,n) = (wa(1)*qw(kk,n) * stlt(kk) &
3176 & + wa(2)*qw(ll,n) * stlt(ll) &
3177 & + wa(3)*qw(l,n) * stlt(l) ) * half
3178
3179 ENDDO
3180!
3181 aa(ll,idw) = aa(ll,idw) + qa(1)
3182 aa(ll,ll) = aa(ll,ll) + qa(2)
3183 aa(ll,l) = aa(ll,l) + qa(3)
3184 bud(ll) = (tx8+rnn(ll)) * half - rnb + tx6 - bud(ll)
3185!
3186 aa(ll,l+1) = bud(ll)
3187!
3188 rnb = trw(2) * vrw(2)
3189!
3190! For the top of the boundary layer
3191!
3192 IF (rnb < zero) THEN
3193 kk = kbl
3194 tem = vt(2) * trw(2)
3195 qa(2) = (rnb - ctl3*tem) * qrpi(kk)
3196 wa(2) = ctl2 * tem * stlt(kk)
3197 ELSE
3198 rnb = zero
3199 qa(2) = zero
3200 wa(2) = zero
3201 ENDIF
3202!
3203 qa(1) = tx2
3204 qa(2) = dof + tx3 - qa(2)
3205 qa(3) = zero
3206!
3207 wa(1) = tx4
3208 wa(2) = dofw + tx5 - wa(2)
3209 wa(3) = zero
3210!
3211 kk = kbl
3212 IF (wvl(kk-1) == wcmin) wa(1) = zero
3213 IF (wvl(kk) == wcmin) wa(2) = zero
3214!
3215 DO ii=1,2
3216 n = kk + ii - 2
3217 aa(kk,n) = (wa(1)*qw(kk-1,n) * stlt(kk-1) &
3218 & + wa(2)*qw(kk,n) * stlt(kk)) * half
3219 ENDDO
3220 fac = half
3221 ll = kbl
3222 l = ll + 1
3223 lm1 = ll - 1
3224 aa(ll,lm1) = aa(ll,lm1) + qa(1)
3225 aa(ll,ll) = aa(ll,ll) + qa(2)
3226 bud(ll) = half*rnn(lm1) - tx6 + rnb - bud(ll)
3227 aa(ll,ll+1) = bud(ll)
3228!
3229!-----SOLVING THE BUDGET EQUATIONS FOR DQR
3230!
3231 DO l=kd1,kbl
3232 lm1 = l - 1
3233 cnvflg = abs(aa(lm1,lm1)) < abs(aa(l,lm1))
3234 DO n=lm1,kbl+1
3235 IF (cnvflg) THEN
3236 tx1 = aa(lm1,n)
3237 aa(lm1,n) = aa(l,n)
3238 aa(l,n) = tx1
3239 ENDIF
3240 ENDDO
3241 tx1 = aa(l,lm1) / aa(lm1,lm1)
3242 DO n=l,kbl+1
3243 aa(l,n) = aa(l,n) - tx1 * aa(lm1,n)
3244 ENDDO
3245 ENDDO
3246!
3247!-----BACK SUBSTITUTION AND CHECK IF THE SOLUTION CONVERGES
3248!
3249 kk = kbl
3250 kk1 = kk + 1
3251 aa(kk,kk1) = aa(kk,kk1) / aa(kk,kk) ! Qr correction !
3252 tx2 = abs(aa(kk,kk1)) * qrpi(kk) ! Error Measure !
3253!
3254 kk = kbl + 1
3255 DO l=kb1,kd,-1
3256 lp1 = l + 1
3257 tx1 = zero
3258 DO n=lp1,kbl
3259 tx1 = tx1 + aa(l,n) * aa(n,kk)
3260 ENDDO
3261 aa(l,kk) = (aa(l,kk) - tx1) / aa(l,l) ! Qr correction !
3262 tx2 = max(tx2, abs(aa(l,kk))*qrpi(l)) ! Error Measure !
3263 ENDDO
3264!
3265! tem = 0.5
3266 if (tx2 > one .and. abs(errq-tx2) > 0.1_kp) then
3267 tem = half
3268!! elseif (tx2 < 0.1) then
3269!! tem = 1.2
3270 else
3271 tem = one
3272 endif
3273!
3274 DO l=kd,kbl
3275! QRP(L) = MAX(QRP(L)+AA(L,KBL+1), QRMIN)
3276 qrp(l) = max(qrp(l)+aa(l,kbl+1)*tem, qrmin)
3277 ENDDO
3278!
3279 IF (itr < itrmin) THEN
3280 tem = abs(errq-tx2)
3281 IF (tem >= errmi2 .AND. tx2 >= errmin) THEN
3282 errq = tx2 ! Further iteration !
3283 ELSE
3284 skpup = .true. ! Converges !
3285 errq = zero ! Rain profile exists!
3286 ENDIF
3287 ELSE
3288 tem = errq - tx2
3289! IF (TEM < ZERO .AND. ERRQ > 0.1_kp) THEN
3290 IF (tem < zero .AND. errq > 0.5_kp) THEN
3291! IF (TEM < ZERO .and. &
3292! & (ntla < numtla .or. ERRQ > 0.5_kp)) THEN
3293 skpup = .true. ! No convergence !
3294 errq = 10.0_kp ! No rain profile!
3295!!!! ELSEIF (ABS(TEM) < ERRMI2 .OR. TX2 < ERRMIN) THEN
3296 ELSEIF (tx2 < errmin) THEN
3297 skpup = .true. ! Converges !
3298 errq = zero ! Rain profile exists!
3299 elseif (tem < zero .and. errq < 0.1_kp) then
3300 skpup = .true.
3301! if (ntla == numtla .or. tem > -0.003) then
3302 errq = zero
3303! else
3304! errq = 10.0
3305! endif
3306 ELSE
3307 errq = tx2 ! Further iteration !
3308
3309! if (itr == itrmu .and. ERRQ > ERRMIN*10 &
3310! & .and. ntla == 1) ERRQ = 10.0
3311 ENDIF
3312 ENDIF
3313!
3314 ENDIF ! SKPUP ENDIF!
3315!
3316 ENDDO ! End of the ITR Loop!!
3317!
3318 IF (errq < 0.1_kp) THEN
3319 ddft = .true.
3320 rnb = - rnb
3321! do l=kd1,kb1-1
3322! if (wvl(l)-wcbase < 1.0E-9) ddft = .false.
3323! enddo
3324 ELSE
3325 ddft = .false.
3326 ENDIF
3327
3328 enddo ! End of ntla loop
3329!
3330! Caution !! Below is an adjustment to rain flux to maintain
3331! conservation of precip!
3332!
3333 IF (ddft) THEN
3334 tx1 = zero
3335 DO l=kd,kb1
3336 tx1 = tx1 + rnf(l)
3337 ENDDO
3338 tx1 = train / (tx1+rnt+rnb)
3339 IF (abs(tx1-one) < 0.2_kp) THEN
3340 rnt = max(rnt*tx1,zero)
3341 rnb = rnb * tx1
3342 DO l=kd,kb1
3343 rnf(l) = rnf(l) * tx1
3344 ENDDO
3345! rain flux adjustment is over
3346
3347 ELSE
3348 ddft = .false.
3349 errq = 10.0_kp
3350 ENDIF
3351 ENDIF
3352!
3353 dof = zero
3354 IF (.NOT. ddft) then
3355 wvlu(kd:kp1) = zero
3356 RETURN ! Rain profile did not converge!
3357 ! No down draft for this case - rerurn
3358 ! ------------------------------------
3359!
3360 else ! rain profile converged - do downdraft calculation
3361 ! ------------------------------------------------
3362
3363 wvlu(kd:kp1) = wvl(kd:kp1) ! save updraft vertical velocity for output
3364
3365!
3366! Downdraft calculation begins
3367! ----------------------------
3368!
3369 DO l=kd,k
3370 wcb(l) = zero
3371 ENDDO
3372!
3373 errq = 10.0_kp
3374! At this point stlt contains inverse of updraft vertical velocity 1/Wu.
3375
3376 kk = max(kb1,kd1)
3377 DO l=kk,k
3378 stlt(l) = stlt(l-1)
3379 ENDDO
3380 tem = stla / bb1 ! this is 2/(pi*radius*grav)
3381!
3382 DO l=kd,k
3383 IF (l <= kbl) THEN
3384 stlt(l) = eta(l) * stlt(l) * tem / ror(l)
3385 ELSE
3386 stlt(l) = zero
3387 ENDIF
3388 ENDDO
3389
3390 rsum1 = zero
3391 rsum2 = zero
3392!
3393 idn(:) = idnmax
3394 DO l=kd,kp1
3395 etd(l) = zero
3396 wvl(l) = zero
3397! QRP(L) = zero
3398 ENDDO
3399 DO l=kd,k
3400 evp(l) = zero
3401 buy(l) = zero
3402 qrp(l+1) = zero
3403 ENDDO
3404 hod(kd) = hol(kd)
3405 qod(kd) = qol(kd)
3406 tx1 = zero
3407!!! TX1 = STLT(KD)*QRB(KD)*ONE ! sigma at the top
3408! TX1 = MIN(STLT(KD)*QRB(KD)*ONE, ONE) ! sigma at the top
3409! TX1 = MIN(STLT(KD)*QRB(KD)*0.5, ONE) ! sigma at the top
3410 rntp = zero
3411 tx5 = tx1
3412 qa(1) = zero
3413!
3414! Here we assume RPART of detrained rain RNT goes to Pd
3415!
3416 IF (rnt > zero) THEN
3417 if (tx1 > zero) THEN
3418 qrp(kd) = (rpart*rnt / (ror(kd)*tx1*gms(kd))) &
3419 & ** (one/1.1364_kp)
3420 else
3421 tx1 = rpart*rnt / (ror(kd)*gms(kd)*qrp(kd)**1.1364_kp)
3422 endif
3423 rntp = (one - rpart) * rnt
3424 buy(kd) = - ror(kd) * tx1 * qrp(kd)
3425 ELSE
3426 qrp(kd) = zero
3427 ENDIF
3428!
3429! L-loop for the downdraft iteration from KD1 to KP1 (bottom surface)
3430!
3431! BUD(KD) = ROR(KD)
3432 idnm = 1
3433 DO l=kd1,kp1
3434
3435 qa(1) = zero
3436 ddlgk = idn(idnm) == idnmax
3437 if (.not. ddlgk) cycle
3438 IF (l <= k) THEN
3439 st1 = one - alfind(l)
3440 wa(1) = alfind(l)*hol(l-1) + st1*hol(l)
3441 wa(2) = alfind(l)*qol(l-1) + st1*qol(l)
3442 wa(3) = alfind(l)*tol(l-1) + st1*tol(l)
3443 qa(2) = alfind(l)*hst(l-1) + st1*hst(l)
3444 qa(3) = alfind(l)*qst(l-1) + st1*qst(l)
3445 ELSE
3446 wa(1) = hol(k)
3447 wa(2) = qol(k)
3448 wa(3) = tol(k)
3449 qa(2) = hst(k)
3450 qa(3) = qst(k)
3451 ENDIF
3452!
3453 fac = two
3454 IF (l == kd1) fac = one
3455
3456 facg = fac * half * gmf5 ! 12/17/97
3457!
3458! DDLGK = IDN(idnm) == 99
3459
3460 bud(kd) = ror(l)
3461
3462 tx1 = tx5
3463 wvl(l) = max(wvl(l-1),one_m1)
3464
3465 qrp(l) = max(qrp(l-1),qrp(l))
3466!
3467! VT(1) = GMS(L-1) * QRP(L-1) ** 0.1364
3468 vt(1) = gms(l-1) * qrpf(qrp(l-1))
3469 rnt = ror(l-1) * (wvl(l-1)+vt(1))*qrp(l-1)
3470!
3471! TEM = MAX(ALM, 2.5E-4) * MAX(ETA(L), 1.0)
3472 tem = max(alm,one_m6) * max(eta(l), one)
3473! TEM = MAX(ALM, 1.0E-5) * MAX(ETA(L), 1.0)
3474 trw(1) = picon*tem*(qrb(l-1)+qrt(l-1))
3475 trw(2) = one / trw(1)
3476!
3477 vrw(1) = half * (gam(l-1) + gam(l))
3478 vrw(2) = one / (vrw(1) + vrw(1))
3479!
3480 tx4 = (qrt(l-1)+qrb(l-1))*(onebg*fac*500.0_kp*eknob)
3481!
3482 dofw = one / (wa(3) * (one + nu*wa(2))) ! 1.0 / TVbar!
3483!
3484 etd(l) = etd(l-1)
3485 hod(l) = hod(l-1)
3486 qod(l) = qod(l-1)
3487!
3488 errq = 10.0_kp
3489
3490!
3491 IF (l <= kbl) THEN
3492 tx3 = stlt(l-1) * qrt(l-1) * (half*fac)
3493 tx8 = stlt(l) * qrb(l-1) * (half*fac)
3494 tx9 = tx8 + tx3
3495 ELSE
3496 tx3 = zero
3497 tx8 = zero
3498 tx9 = zero
3499 ENDIF
3500!
3501 tem = wvl(l-1) + vt(1)
3502 IF (tem > zero) THEN
3503 tem1 = one / (tem*ror(l-1))
3504 tx3 = vt(1) * tem1 * ror(l-1) * tx3
3505 tx6 = tx1 * tem1
3506 ELSE
3507 tx6 = one
3508 ENDIF
3509!
3510 IF (l == kd1) THEN
3511 IF (rnt > zero) THEN
3512 tem = max(qrp(l-1),qrp(l))
3513 wvl(l) = tx1 * tem * qrb(l-1)*(facg*5.0_kp)
3514 ENDIF
3515 wvl(l) = max(one_m2, wvl(l))
3516 trw(1) = trw(1) * half
3517 trw(2) = trw(2) + trw(2)
3518 ELSE
3519 IF (ddlgk) evp(l-1) = evp(l-2)
3520 ENDIF
3521!
3522! No downdraft above level IDH
3523!
3524
3525 IF (l < idh) THEN
3526
3527 etd(l) = zero
3528 hod(l) = wa(1)
3529 qod(l) = wa(2)
3530 evp(l-1) = zero
3531 wvl(l) = zero
3532 qrp(l) = zero
3533 buy(l) = zero
3534 tx5 = tx9
3535 errq = zero
3536 rntp = rntp + rnt * tx1
3537 rnt = zero
3538 wcb(l-1) = zero
3539
3540! ENDIF
3541! BUD(KD) = ROR(L)
3542!
3543! Iteration loop for a given level L begins
3544!
3545 else
3546 DO itr=1,itrmd
3547!
3548! cnvflg = DDLGK .AND. (ERRQ > ERRMIN)
3549 cnvflg = errq > errmin
3550 IF (cnvflg) THEN
3551!
3552! VT(1) = GMS(L) * QRP(L) ** 0.1364
3553 vt(1) = gms(l) * qrpf(qrp(l))
3554 tem = wvl(l) + vt(1)
3555!
3556 IF (tem > zero) THEN
3557 st1 = ror(l) * tem * qrp(l) + rnt
3558 IF (st1 /= zero) st1 = two * evp(l-1) / st1
3559 tem1 = one / (tem*ror(l))
3560 tem2 = vt(1) * tem1 * ror(l) * tx8
3561 ELSE
3562 tem1 = zero
3563 tem2 = tx8
3564 st1 = zero
3565 ENDIF
3566!
3567 st2 = tx5
3568 tem = ror(l)*wvl(l) - ror(l-1)*wvl(l-1)
3569 if (tem > zero) then
3570 tx5 = (tx1 - st1 + tem2 + tx3)/(one+tem*tem1)
3571 else
3572 tx5 = tx1 - tem*tx6 - st1 + tem2 + tx3
3573 endif
3574 tx5 = max(tx5,zero)
3575 tx5 = half * (tx5 + st2)
3576!
3577! qqq = 1.0 + tem * tem1 * (1.0 - sialf)
3578!
3579! if (qqq > 0.0) then
3580! TX5 = (TX1 - sialf*tem*tx6 - ST1 + TEM2 + TX3) / qqq
3581! else
3582! TX5 = (TX1 - tem*tx6 - ST1 + TEM2 + TX3)
3583! endif
3584!
3585 tem1 = etd(l)
3586 etd(l) = ror(l) * tx5 * max(wvl(l),zero)
3587!
3588 if (etd(l) > zero) etd(l) = half * (etd(l) + tem1)
3589!
3590
3591 del_eta = etd(l) - etd(l-1)
3592
3593! TEM = DEL_ETA * TRW(2)
3594! TEM2 = MAX(MIN(TEM, 1.0), -1.0)
3595! IF (ABS(TEM) > 1.0 .AND. ETD(L) > 0.0 ) THEN
3596! DEL_ETA = TEM2 * TRW(1)
3597! ETD(L) = ETD(L-1) + DEL_ETA
3598! ENDIF
3599! IF (WVL(L) > 0.0) TX5 = ETD(L) / (ROR(L)*WVL(L))
3600!
3601 erre = etd(l) - tem1
3602!
3603 tem = max(abs(del_eta), trw(1))
3604 tem2 = del_eta / tem
3605 tem1 = sqrt(max((tem+del_eta)*(tem-del_eta),zero))
3606! TEM1 = SQRT(MAX((TRW(1)+DEL_ETA)*(TRW(1)-DEL_ETA),0.0))
3607
3608 edz = (half + asin(tem2)*piinv)*del_eta + tem1*piinv
3609
3610 ddz = edz - del_eta
3611 wcb(l-1) = etd(l) + ddz
3612!
3613 tem1 = hod(l)
3614 IF (del_eta > zero) THEN
3615 qqq = one / (etd(l) + ddz)
3616 hod(l) = (etd(l-1)*hod(l-1) + del_eta*hol(l-1) &
3617 & + ddz*wa(1)) * qqq
3618 qod(l) = (etd(l-1)*qod(l-1) + del_eta*qol(l-1) &
3619 & + ddz*wa(2)) * qqq
3620 ELSEif((etd(l-1) + edz) > zero) then
3621 qqq = one / (etd(l-1) + edz)
3622 hod(l) = (etd(l-1)*hod(l-1) + edz*wa(1)) * qqq
3623 qod(l) = (etd(l-1)*qod(l-1) + edz*wa(2)) * qqq
3624 ENDIF
3625 errh = hod(l) - tem1
3626 errq = abs(errh/hod(l)) + abs(erre/max(etd(l),one_m5))
3627
3628 dof = ddz
3629 vt(2) = qqq
3630!
3631 ddz = dof
3632 tem4 = qod(l)
3633 tem1 = vrw(1)
3634!
3635 qhs = qa(3) + half * (gaf(l-1)+gaf(l)) * (hod(l)-qa(2))
3636!
3637! First iteration !
3638!
3639 st2 = prl(l) * (qhs + tem1 * (qhs-qod(l)))
3640 tem2 = ror(l) * qrp(l)
3641 CALL qrabf(tem2,qraf,qrbf)
3642 tem6 = tx5 * (1.6_kp + 124.9_kp * qraf) * qrbf * tx4
3643!
3644 ce = tem6 * st2 / ((5.4e5_kp*st2 + 2.55e6_kp)*(etd(l)+ddz))
3645!
3646 tem2 = - ((one+tem1)*(qhs+ce) + tem1*qod(l))
3647 tem3 = (one + tem1) * qhs * (qod(l)+ce)
3648 tem = max(tem2*tem2 - four*tem1*tem3,zero)
3649 qod(l) = max(tem4, (- tem2 - sqrt(tem)) * vrw(2))
3650!
3651!
3652! second iteration !
3653!
3654 st2 = prl(l) * (qhs + tem1 * (qhs-qod(l)))
3655 ce = tem6 * st2 / ((5.4e5_kp*st2 + 2.55e6_kp)*(etd(l)+ddz))
3656! CEE = CE * (ETD(L)+DDZ)
3657!
3658
3659
3660 tem2 = - ((one+tem1)*(qhs+ce) + tem1*tem4)
3661 tem3 = (one + tem1) * qhs * (tem4+ce)
3662 tem = max(tem2*tem2 - four*tem1*tem3,zero)
3663 qod(l) = max(tem4, (- tem2 - sqrt(tem)) * vrw(2))
3664! Evaporation in Layer L-1
3665!
3666 evp(l-1) = (qod(l)-tem4) * (etd(l)+ddz)
3667! Calculate Pd (L+1/2)
3668 qa(1) = tx1*rnt + rnf(l-1) - evp(l-1)
3669!
3670 if (qa(1) > zero) then
3671 IF (etd(l) > zero) THEN
3672 tem = qa(1) / (etd(l)+ror(l)*tx5*vt(1))
3673 qrp(l) = max(tem,zero)
3674 ELSEIF (tx5 > zero) THEN
3675 qrp(l) = (max(zero,qa(1)/(ror(l)*tx5*gms(l)))) &
3676 & ** (one/1.1364_kp)
3677 ELSE
3678 qrp(l) = zero
3679 ENDIF
3680 else
3681 qrp(l) = half * qrp(l)
3682 endif
3683! Compute Buoyancy
3684 tem1 = wa(3) + (hod(l)-wa(1)-alhl*(qod(l)-wa(2))) &
3685 & * onebcp
3686
3687 tem1 = tem1 * (one + nu*qod(l))
3688 ror(l) = cmpor * prl(l) / tem1
3689 tem1 = tem1 * dofw
3690!!! TEM1 = TEM1 * (1.0 + NU*QOD(L)) * DOFW
3691
3692 buy(l) = (tem1 - one - qrp(l)) * ror(l) * tx5
3693! Compute W (L+1/2)
3694
3695 tem1 = wvl(l)
3696 wvl(l) = vt(2) * (etd(l-1)*wvl(l-1) - facg &
3697 & * (buy(l-1)*qrt(l-1)+buy(l)*qrb(l-1)))
3698
3699!
3700 if (wvl(l) < zero) then
3701! WVL(L) = max(wvl(l), 0.1*tem1)
3702! WVL(L) = 0.5*tem1
3703! WVL(L) = 0.1*tem1
3704! WVL(L) = 0.0
3705 wvl(l) = 1.0e-10_kp
3706 else
3707 wvl(l) = half*(wvl(l)+tem1)
3708 endif
3709
3710!
3711! WVL(L) = max(0.5*(WVL(L)+TEM1), 0.0)
3712
3713 errw = wvl(l) - tem1
3714!
3715 errq = errq + abs(errw/max(wvl(l),one_m5))
3716
3717
3718! IF (ITR >= MIN(ITRMIN,ITRMD/2)) THEN
3719 IF (itr >= min(itrmnd,itrmd/2)) THEN
3720
3721 IF (etd(l-1) == zero .AND. errq > 0.2_kp) THEN
3722
3723 ror(l) = bud(kd)
3724 etd(l) = zero
3725 wvl(l) = zero
3726 errq = zero
3727 hod(l) = wa(1)
3728 qod(l) = wa(2)
3729! TX5 = TX1 + TX9
3730 if (l <= kbl) then
3731 tx5 = tx9
3732 else
3733 tx5 = (stlt(kb1) * qrt(kb1) &
3734 & + stlt(kbl) * qrb(kb1)) * (0.5_kp*fac)
3735 endif
3736
3737 evp(l-1) = zero
3738 tem = max(tx1*rnt+rnf(l-1),zero)
3739 qa(1) = tem - evp(l-1)
3740! IF (QA(1) > 0.0) THEN
3741
3742 qrp(l) = (qa(1) / (ror(l)*tx5*gms(l))) &
3743 & ** (one/1.1364_kp)
3744! endif
3745 buy(l) = - ror(l) * tx5 * qrp(l)
3746 wcb(l-1) = zero
3747 ENDIF
3748!
3749 del_eta = etd(l) - etd(l-1)
3750 IF(del_eta < zero .AND. errq > 0.1_kp) THEN
3751 ror(l) = bud(kd)
3752 etd(l) = zero
3753 wvl(l) = zero
3754!!!!! TX5 = TX1 + TX9
3755 cldfrd(l-1) = tx5
3756!
3757 del_eta = - etd(l-1)
3758 edz = zero
3759 ddz = -del_eta
3760 wcb(l-1) = ddz
3761!
3762 hod(l) = hod(l-1)
3763 qod(l) = qod(l-1)
3764!
3765 tem4 = qod(l)
3766 tem1 = vrw(1)
3767!
3768 qhs = qa(3) + half * (gaf(l-1)+gaf(l)) &
3769 & * (hod(l)-qa(2))
3770
3771!
3772! First iteration !
3773!
3774 st2 = prl(l) * (qhs + tem1 * (qhs-qod(l)))
3775 tem2 = ror(l) * qrp(l-1)
3776 CALL qrabf(tem2,qraf,qrbf)
3777 tem6 = tx5 * (1.6_kp + 124.9_kp * qraf) * qrbf * tx4
3778!
3779 ce = tem6*st2/((5.4e5_kp*st2 + 2.55e6_kp)*(etd(l)+ddz))
3780!
3781
3782 tem2 = - ((one+tem1)*(qhs+ce) + tem1*qod(l))
3783 tem3 = (one + tem1) * qhs * (qod(l)+ce)
3784 tem = max(tem2*tem2 -four*tem1*tem3,zero)
3785 qod(l) = max(tem4, (- tem2 - sqrt(tem)) * vrw(2))
3786!
3787! second iteration !
3788!
3789 st2 = prl(l) * (qhs + tem1 * (qhs-qod(l)))
3790 ce = tem6*st2/((5.4e5_kp*st2 + 2.55e6_kp)*(etd(l)+ddz))
3791! CEE = CE * (ETD(L)+DDZ)
3792!
3793
3794
3795 tem2 = - ((one+tem1)*(qhs+ce) + tem1*tem4)
3796 tem3 = (one + tem1) * qhs * (tem4+ce)
3797 tem = max(tem2*tem2 -four*tem1*tem3,zero)
3798 qod(l) = max(tem4, (- tem2 - sqrt(tem)) * vrw(2))
3799
3800! Evaporation in Layer L-1
3801!
3802 evp(l-1) = (qod(l)-tem4) * (etd(l)+ddz)
3803
3804! Calculate Pd (L+1/2)
3805! RNN(L-1) = TX1*RNT + RNF(L-1) - EVP(L-1)
3806
3807 qa(1) = tx1*rnt + rnf(l-1)
3808 evp(l-1) = min(evp(l-1), qa(1))
3809 qa(1) = qa(1) - evp(l-1)
3810 qrp(l) = zero
3811
3812!
3813! IF (QA(1) > 0.0) THEN
3814!! RNS(L-1) = QA(1)
3815!!! tx5 = tx9
3816! QRP(L) = (QA(1) / (ROR(L)*TX5*GMS(L))) &
3817! & ** (1.0/1.1364)
3818! endif
3819! ERRQ = 0.0
3820! Compute Buoyancy
3821! TEM1 = WA(3)+(HOD(L)-WA(1)-ALHL*(QOD(L)-WA(2))) &
3822! & * (1.0/CP)
3823! TEM1 = TEM1 * (1.0 + NU*QOD(L)) * DOFW
3824! BUY(L) = (TEM1 - 1.0 - QRP(L)) * ROR(L) * TX5
3825!
3826! IF (QA(1) > 0.0) RNS(L) = QA(1)
3827
3828 IF (l .LE. k) THEN
3829 rns(l) = qa(1)
3830 qa(1) = zero
3831 ENDIF
3832 tx5 = tx9
3833 errq = zero
3834 qrp(l) = zero
3835 buy(l) = zero
3836!
3837 ENDIF
3838 ENDIF
3839 ENDIF
3840!
3841 ENDDO ! End of the iteration loop for a given L!
3842 IF (l <= k) THEN
3843 IF (etd(l-1) == zero .AND. errq > 0.1_kp .and. l <= kbl) THEN
3844!!! & .AND. ERRQ > ERRMIN*10.0 .and. l <= kbl) THEN
3845! & .AND. ERRQ > ERRMIN*10.0) THEN
3846 ror(l) = bud(kd)
3847 hod(l) = wa(1)
3848 qod(l) = wa(2)
3849 tx5 = tx9 ! Does not make too much difference!
3850! TX5 = TX1 + TX9
3851 evp(l-1) = zero
3852! EVP(L-1) = CEE * (1.0 - qod(l)/qa(3))
3853 qa(1) = tx1*rnt + rnf(l-1)
3854 evp(l-1) = min(evp(l-1), qa(1))
3855 qa(1) = qa(1) - evp(l-1)
3856
3857! QRP(L) = 0.0
3858! if (tx5 == 0.0 .or. gms(l) == 0.0) then
3859! write(0,*)' Ctx5=',tx5,' gms=',gms(l),' ror=',ror(l) &
3860! &, ' L=',L,' QA=',QA(1),' tx1=',tx1,' tx9=',tx9 &
3861! &, ' kbl=',kbl,' etd1=',etd(l-1),' DEL_ETA=',DEL_ETA
3862! endif
3863! IF (QA(1) > 0.0) THEN
3864
3865 qrp(l) = (qa(1) / (ror(l)*tx5*gms(l))) &
3866 & ** (one/1.1364_kp)
3867! ENDIF
3868 etd(l) = zero
3869 wvl(l) = zero
3870 st1 = one - alfind(l)
3871
3872 errq = zero
3873 buy(l) = - ror(l) * tx5 * qrp(l)
3874 wcb(l-1) = zero
3875 ENDIF
3876 ENDIF
3877!
3878 ll = min(idn(idnm), kp1)
3879 IF (errq < one .AND. l <= ll) THEN
3880 IF (etd(l-1) > zero .AND. etd(l) == zero) THEN
3881 idn(idnm) = l
3882 wvl(l) = zero
3883 if (l < kbl .or. tx5 > zero) idnm = idnm + 1
3884 errq = zero
3885 ENDIF
3886 if (etd(l) == zero .and. l > kbl) then
3887 idn(idnm) = l
3888 if (tx5 > zero) idnm = idnm + 1
3889 endif
3890 ENDIF
3891!
3892! If downdraft properties are not obtainable, (i.e.solution does
3893! not converge) , no downdraft is assumed
3894!
3895! IF (ERRQ > ERRMIN*100.0 .AND. IDN(idnm) == 99) &
3896 IF (errq > 0.1_kp .AND. idn(idnm) == idnmax) ddft = .false.
3897!
3898 dof = zero
3899 IF (.NOT. ddft) RETURN
3900!
3901! if (ddlgk .or. l .le. idn(idnm)) then
3902! rsum2 = rsum2 + evp(l-1)
3903! write(0,*)' rsum1=',rsum1,' rsum2=',rsum2,' L=',L,' qa=',qa(1)&
3904! &, ' evp=',evp(l-1)
3905! else
3906! rsum1 = rsum1 + rnf(l-1)
3907! write(0,*)' rsum1=',rsum1,' rsum2=',rsum2,' L=',L,' rnf=', &
3908! & rnf(l-1)
3909! endif
3910
3911 endif ! if (l < idh)
3912 ENDDO ! End of the L Loop of downdraft !
3913
3914 tx1 = zero
3915
3916 dof = qa(1)
3917!
3918! write(0,*)' dof=',dof,' rntp=',rntp,' rnb=',rnb
3919! write(0,*)' total=',(rsum1+dof+rntp+rnb)
3920!
3921 dof = max(dof, zero)
3922 rnn(kd) = rntp
3923 tx1 = evp(kd)
3924 tx2 = rntp + rnb + dof
3925
3926 ii = idh
3927 IF (ii >= kd1+1) THEN
3928 rnn(kd) = rnn(kd) + rnf(kd)
3929 tx2 = tx2 + rnf(kd)
3930 rnn(ii-1) = zero
3931 tx1 = evp(ii-1)
3932 ENDIF
3933 DO l=kd,k
3934 ii = idh
3935
3936 IF (l > kd1 .AND. l < ii) THEN
3937 rnn(l-1) = rnf(l-1)
3938 tx2 = tx2 + rnn(l-1)
3939 ELSEIF (l >= ii .AND. l < idn(idnm)) THEN
3940 rnn(l) = rns(l)
3941 tx2 = tx2 + rnn(l)
3942 tx1 = tx1 + evp(l)
3943 ELSEIF (l >= idn(idnm)) THEN
3944 etd(l+1) = zero
3945 hod(l+1) = zero
3946 qod(l+1) = zero
3947 evp(l) = zero
3948 rnn(l) = rnf(l) + rns(l)
3949 tx2 = tx2 + rnn(l)
3950 ENDIF
3951 ENDDO
3952!
3953! For Downdraft case the rain is that falls thru the bottom
3954
3955 l = kbl
3956
3957 rnn(l) = rnn(l) + rnb
3958 cldfrd(l) = tx5
3959
3960!
3961! Caution !! Below is an adjustment to rain flux to maintain
3962! conservation of precip!
3963
3964!
3965 IF (tx1 > zero) THEN
3966 tx1 = (train - tx2) / tx1
3967 ELSE
3968 tx1 = zero
3969 ENDIF
3970
3971 DO l=kd,k
3972 evp(l) = evp(l) * tx1
3973 ENDDO
3974
3975 ENDIF ! if (.not. DDFT) loop endif
3976!
3977!***********************************************************************
3978!***********************************************************************
3979
3980 RETURN
3981 end subroutine ddrft
3982
3984 SUBROUTINE qsatcn(TT,P,Q,DQDT)
3985!
3986 USE funcphys , ONLY : fpvs
3987
3988 implicit none
3989!
3990 real(kind=kind_phys) tt, p, q, dqdt
3991!
3992! real(kind=kind_phys), parameter :: ZERO=0.0, ONE=1.0 &
3993! &, rvi=one/rv, facw=CVAP-CLIQ &
3994! &, faci=CVAP-CSOL, hsub=alhl+alhf &
3995! &, tmix=TTP-20.0 &
3996! &, DEN=one/(TTP-TMIX)
3997!
3998 real(kind=kind_phys) es, d, hlorv, w
3999!
4000! es = 10.0 * fpvs(tt) ! fpvs is in centibars!
4001 es = min(p, 0.01_kp * fpvs(tt)) ! fpvs is in Pascals!
4002! D = one / max(p+epsm1*es,ONE_M10)
4003 d = one / (p+epsm1*es)
4004!
4005 q = min(eps*es*d, one)
4006!
4007 w = max(zero, min(one, (tt - tmix)*den))
4008 hlorv = ( w * (alhl + facw * (tt-ttp)) &
4009 & + (one-w) * (hsub + faci * (tt-ttp)) ) * rvi
4010 dqdt = p * q * hlorv * d / (tt*tt)
4011!
4012 return
4013 end subroutine qsatcn
4014
4016 SUBROUTINE angrad(PRES, ALM, AL2, TLA)
4017 implicit none
4018
4019 real(kind=kind_phys) pres, alm, al2, tla, tem
4020!
4021 integer i
4022!
4023 IF (tla < 0.0_kp) THEN
4024 IF (pres <= plac(1)) THEN
4025 tla = tlac(1)
4026 ELSEIF (pres <= plac(2)) THEN
4027 tla = tlac(2) + (pres-plac(2))*tlbpl(1)
4028 ELSEIF (pres <= plac(3)) THEN
4029 tla = tlac(3) + (pres-plac(3))*tlbpl(2)
4030 ELSEIF (pres <= plac(4)) THEN
4031 tla = tlac(4) + (pres-plac(4))*tlbpl(3)
4032 ELSEIF (pres <= plac(5)) THEN
4033 tla = tlac(5) + (pres-plac(5))*tlbpl(4)
4034 ELSEIF (pres <= plac(6)) THEN
4035 tla = tlac(6) + (pres-plac(6))*tlbpl(5)
4036 ELSEIF (pres <= plac(7)) THEN
4037 tla = tlac(7) + (pres-plac(7))*tlbpl(6)
4038 ELSEIF (pres <= plac(8)) THEN
4039 tla = tlac(8) + (pres-plac(8))*tlbpl(7)
4040 ELSE
4041 tla = tlac(8)
4042 ENDIF
4043 ENDIF
4044 IF (pres >= refp(1)) THEN
4045 tem = refr(1)
4046 ELSEIF (pres >= refp(2)) THEN
4047 tem = refr(1) + (pres-refp(1)) * drdp(1)
4048 ELSEIF (pres >= refp(3)) THEN
4049 tem = refr(2) + (pres-refp(2)) * drdp(2)
4050 ELSEIF (pres >= refp(4)) THEN
4051 tem = refr(3) + (pres-refp(3)) * drdp(3)
4052 ELSEIF (pres >= refp(5)) THEN
4053 tem = refr(4) + (pres-refp(4)) * drdp(4)
4054 ELSEIF (pres >= refp(6)) THEN
4055 tem = refr(5) + (pres-refp(5)) * drdp(5)
4056 ELSE
4057 tem = refr(6)
4058 ENDIF
4059!
4060 tem = 2.0e-4_kp / tem
4061 al2 = min(4.0_kp*tem, max(alm, tem))
4062!
4063 RETURN
4064 end subroutine angrad
4065
4067 SUBROUTINE setqrp
4068 implicit none
4069
4070 real(kind=kind_phys) tem2,tem1,x,xinc,xmax,xmin
4071 integer jx
4072! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4073! XMIN = 1.0E-6
4074 xmin = 0.0_kp
4075 xmax = 5.0_kp
4076 xinc = (xmax-xmin)/(nqrp-1)
4077 c2xqrp = one / xinc
4078 c1xqrp = one - xmin*c2xqrp
4079 tem1 = 0.001_kp ** 0.2046_kp
4080 tem2 = 0.001_kp ** 0.525_kp
4081 DO jx=1,nqrp
4082 x = xmin + (jx-1)*xinc
4083 tbqrp(jx) = x ** 0.1364_kp
4084 tbqra(jx) = tem1 * x ** 0.2046_kp
4085 tbqrb(jx) = tem2 * x ** 0.525_kp
4086 ENDDO
4087! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4088 RETURN
4089 end subroutine setqrp
4090
4092 SUBROUTINE qrabf(QRP,QRAF,QRBF)
4093 implicit none
4094!
4095 real(kind=kind_phys) qrp, qraf, qrbf, xj, real_nqrp
4096 INTEGER JX
4097! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4098 real_nqrp = real(nqrp)
4099 xj = min(max(c1xqrp+c2xqrp*qrp,one),real_nqrp)
4100 jx = min(xj,nqrp-one)
4101 xj = xj - jx
4102 qraf = tbqra(jx) + xj * (tbqra(jx+1)-tbqra(jx))
4103 qrbf = tbqrb(jx) + xj * (tbqrb(jx+1)-tbqrb(jx))
4104! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4105 RETURN
4106 end subroutine qrabf
4107
4109 SUBROUTINE setvtp
4110 implicit none
4111
4112 real(kind=kind_phys), parameter :: vtpexp=-0.3636_kp
4113 real(kind=kind_phys) xinc,x,xmax,xmin
4114 integer jx
4115! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4116 xmin = 0.05_kp
4117 xmax = 1.5_kp
4118 xinc = (xmax-xmin)/(nvtp-1)
4119 c2xvtp = one / xinc
4120 c1xvtp = one - xmin*c2xvtp
4121 DO jx=1,nvtp
4122 x = xmin + (jx-1)*xinc
4123 tbvtp(jx) = x ** vtpexp
4124 ENDDO
4125! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4126 RETURN
4127 end subroutine setvtp
4128!
4130 real(kind=kind_phys) FUNCTION qrpf(QRP)
4131!
4132 implicit none
4133
4134 real(kind=kind_phys) qrp, xj, real_nqrp
4135 INTEGER jx
4136! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4137 real_nqrp = real(nqrp)
4138 xj = min(max(c1xqrp+c2xqrp*qrp,one),real_nqrp)
4139! XJ = MIN(MAX(C1XQRP+C2XQRP*QRP,ONE),FLOAT(NQRP))
4140 jx = min(xj,nqrp-one)
4141 qrpf = tbqrp(jx) + (xj-jx) * (tbqrp(jx+1)-tbqrp(jx))
4142! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4143 RETURN
4144 end function qrpf
4145
4147 real(kind=kind_phys) FUNCTION vtpf(ROR)
4148!
4149 implicit none
4150 real(kind=kind_phys) ror, xj, real_nvtp
4151 INTEGER jx
4152! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4153 real_nvtp = real(nvtp)
4154 xj = min(max(c1xvtp+c2xvtp*ror,one),real_nvtp)
4155 jx = min(xj,nvtp-one)
4156 vtpf = tbvtp(jx) + (xj-jx) * (tbvtp(jx+1)-tbvtp(jx))
4157! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4158 RETURN
4159 end function vtpf
4160
4162 real(kind=kind_phys) FUNCTION clf(PRATE)
4163!
4164 implicit none
4165 real(kind=kind_phys) prate
4166!
4167 real (kind=kind_phys), parameter :: ccf1=0.30_kp, ccf2=0.09_kp &
4168 &, ccf3=0.04_kp, ccf4=0.01_kp &
4169 &, pr1=1.0_kp, pr2=5.0_kp &
4170 &, pr3=20.0_kp
4171!
4172 if (prate < pr1) then
4173 clf = ccf1
4174 elseif (prate < pr2) then
4175 clf = ccf2
4176 elseif (prate < pr3) then
4177 clf = ccf3
4178 else
4179 clf = ccf4
4180 endif
4181!
4182 RETURN
4183 end function clf
4184 end module rascnv
subroutine, public rascnv_run(im, k, itc, ntc, ntr, dt, dtf, ccwf, area, dxmin, dxinv, psauras, prauras, wminras, dlqf, flipv, me, rannum, nrcm, mp_phys, mp_phys_mg, ntk, kdt, rhc, tin, qin, uin, vin, ccin, fscav, prsi, prsl, prsik, prslk, phil, phii, kpbl, cdrag, rainc, kbot, ktop, kcnv, ddvel, ud_mf, dd_mf, dt_mf, qlcn, qicn, w_upi, cf_upi, cnv_mfd, cnv_dqldt, clcn, cnv_fice, cnv_ndrop, cnv_nice, errmsg, errflg)
Definition rascnv.F90:296
subroutine setqrp
Definition rascnv.F90:4068
subroutine cloud(k, kp1, kd, ntrc, kblmx, kblmn, fracbl, max_neg_bouy, vsmooth, aw_scal, revap, wrkfun, calkbl, crtfun, dt, kdt, tla, dpd, alfint, rhfacl, rhfacs, area, ccwf, cd, trcfac, alfind, rhc_ls, phil, phih, prs, prsm, sgcs, toi, qoi, roi, qli, qii, kpbl, dsfc, tcu, qcu, rcu, pcu, flx, flxd, cup, wfnc, fscav_, trcmin, ntk, c0, qw0, c0i, qi0, dlq_fac)
Definition rascnv.F90:1090
real(kind=kind_phys) function qrpf(qrp)
Definition rascnv.F90:4131
real(kind=kind_phys) function clf(prate)
Definition rascnv.F90:4163
subroutine, public rascnv_init(me, dt, con_g, con_cp, con_rd, con_rv, con_hvap, con_hfus, con_fvirt, con_t0c, con_ttp, con_cvap, con_cliq, con_csol, con_eps, con_epsm1, errmsg, errflg)
The subroutine initializes rascnv.
Definition rascnv.F90:118
subroutine qrabf(qrp, qraf, qrbf)
Definition rascnv.F90:4093
subroutine setvtp
Definition rascnv.F90:4110
real(kind=kind_phys) function vtpf(ror)
Definition rascnv.F90:4148
subroutine angrad(pres, alm, al2, tla)
Definition rascnv.F90:4017
subroutine ddrft(k, kp1, kd, tla, alfind, wcbase, tol, qol, hol, prl, qst, hst, gam, gaf, qrb, qrt, buy, kbl, idh, eta, rnn, etai, alm, wfn, train, ddft, etd, hod, qod, evp, dof, cldfrd, wcb, gms, gsd, ghd, wvlu)
Definition rascnv.F90:2716
subroutine qsatcn(tt, p, q, dqdt)
Definition rascnv.F90:3985
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26