CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
cs_conv.F90
1
3
4module cs_conv
5!---------------------------------------------------------------------------------
6! Purpose:
7!
9! Purpose:
10!
14!---------------------------------------------------------------------------------
15!
16 use machine , only : kind_phys
17 use physcons, only : cp => con_cp, grav => con_g, &
18 & rair => con_rd, rvap => con_rv, &
19 & cliq => con_cliq, cvap => con_cvap, &
20 & epsv => con_eps, epsvm1 => con_epsm1, &
21 & epsvt => con_fvirt, &
22 & el => con_hvap, emelt => con_hfus, t0c => con_t0c
23 use funcphys, only : fpvs ! this is saturation vapor pressure in funcphys.f
24
25
26 implicit none
27
28 private ! Make default type private to the module
29
30 real(kind_phys), parameter :: zero=0.0d0, one=1.0d0, half=0.5d0
31 real(kind_phys), parameter :: cpoel=cp/el, cpoesub=cp/(el+emelt), esubocp=1.0/cpoesub, &
32 elocp=el/cp, oneocp=one/cp, gocp=grav/cp, gravi=one/grav,&
33 emeltocp=emelt/cp, cpoemelt=cp/emelt, epsln=1.e-10_kind_phys
34
35 real(kind_phys), parameter :: fact1=(cvap-cliq)/rvap, fact2=el/rvap-fact1*t0c
36
37 logical, parameter :: adjustp=.true.
38! logical, parameter :: adjustp=.false.
39
40! Tuning parameters set from namelist
41!
42! real(kind_phys), parameter, public :: CLMD = 0.60, & !< entrainment efficiency (now thru argument)
43 real(kind_phys), parameter, public :: &
44 pa=0.15, & !< factor for buoyancy to affect updraft velocity
45 cpres = 0.55, &
46 alp0 = 5.0e7, &
47! ALP0 = 8.0e7, & !< alpha parameter in prognostic closure
48! CLMP = (one-CLMD)*(PA+PA), &
49! CLMDPA = CLMD*PA, &
50 spblmin=0.05, &
51 spblmax=0.30, &
52! spblcrit=0.03, & !< minimum cloudbase height in p/ps
53! spblcrit=0.035,& !< minimum cloudbase height in p/ps
54! spblcrit=0.025,& !< minimum cloudbase height in p/ps
55 cincrit= -150.0
56! cincrit= -120.0
57! cincrit= -100.0
58
59!DD precz0 and preczh control partitioning of water between detrainment
60!DD and precipitation. Decrease for more precip
61
62 real(kind_phys), public :: precz0, preczh, clmd, clmp, clmdpa
63 real(kind_phys), public, parameter :: c0t=0.002, d0t=0.002
64!
65! Private data
66!
67 real(kind_phys), parameter :: unset_kind_phys = -999._kind_phys ! missing value
68!
69 integer :: iulog
70 !DD Note - see if I can find corresponding variable in a GFS module
71!
72! Shared variables
73!
74 integer, parameter :: iti = 2, itl = 3
75
76 integer, save, dimension(50) :: imfxr
86
87! PUBLIC: interfaces
88!
89 public cs_conv_run ! CS scheme main driver
90
91 contains
92
154 subroutine cs_conv_run( IJSDIM , KMAX , ntracp1 , NN, &
155 NTR , nctp , & !DD dimensions
156 otspt , lat , kdt , &
157 t , q , rain1 , clw , &
158 zm , zi , pap , paph , &
159 delta , delti , ud_mf , dd_mf , dt_mf, &
160 u , v , fscav , fswtr, &
161 cbmfx , mype , wcbmaxm , precz0in, preczhin, &
162 clmdin , sigma , do_aw , do_awdd , flx_form, &
163 lprnt , ipr, kcnv, &
164 QLCN, QICN, w_upi, cf_upi, CNV_MFD, & ! for coupling to MG microphysics
165 CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE, &
166 mp_phys,errmsg,errflg)
167
168
169 implicit none
170!
171! input arguments
172!
173 INTEGER, INTENT(IN) :: ijsdim, kmax, ntracp1, nn, ntr, mype, nctp, mp_phys, kdt, lat !! DD, for GFS, pass in
174 logical, intent(in) :: otspt(:,:) ! otspt(:,1) - on/off switch for tracer transport by updraft and
175 ! downdraft. should not include subgrid PDF and turbulence
176 ! otspt(:,2) - on/off switch for tracer transport by subsidence
177 ! should include subgrid PDF and turbulence
178
179 real(kind_phys), intent(inout) :: t(:,:) ! temperature at mid-layer (K)
180 real(kind_phys), intent(inout) :: q(:,:) ! water vapor array including moisture (kg/kg)
181 real(kind_phys), intent(inout) :: clw(:,:,:) ! tracer array including cloud condensate (kg/kg)
182 real(kind_phys), intent(in) :: pap(:,:) ! pressure at mid-layer (Pa)
183 real(kind_phys), intent(in) :: paph(:,:) ! pressure at boundaries (Pa)
184 real(kind_phys), intent(in) :: zm(:,:) ! geopotential at mid-layer (m)
185 real(kind_phys), intent(in) :: zi(:,:) ! geopotential at boundaries (m)
186 real(kind_phys), intent(in) :: fscav(:), fswtr(:), wcbmaxm(:)
187 real(kind_phys), intent(in) :: precz0in, preczhin, clmdin
188! added for cs_convr
189 real(kind_phys), intent(inout) :: u(:,:) ! zonal wind at mid-layer (m/s)
190 real(kind_phys), intent(inout) :: v(:,:) ! meridional wind at mid-layer (m/s)
191
192 real(kind_phys), intent(in) :: delta ! physics time step
193 real(kind_phys), intent(in) :: delti ! dynamics time step (model time increment in seconds)
194 logical, intent(in) :: do_aw, do_awdd, flx_form
195!
196! modified arguments
197!
198 real(kind_phys), intent(inout), optional :: cbmfx(:,:) ! cloud base mass flux (kg/m2/s)
199!
200! output arguments
201!
202! updraft, downdraft, and detrainment mass flux (kg/m2/s)
203 real(kind_phys), intent(inout), dimension(:,:), optional :: ud_mf
204 real(kind_phys), intent(inout), dimension(:,:) :: dd_mf, dt_mf
205
206 real(kind_phys), intent(out) :: rain1(:) ! lwe thickness of deep convective precipitation amount (m)
207! GJF* These variables are conditionally allocated depending on whether the
208! Morrison-Gettelman microphysics is used, so they must be declared
209! using assumed shape.
210 real(kind_phys), intent(out), dimension(:,:), optional :: qlcn, qicn, w_upi,cnv_mfd, &
211 cnv_dqldt, clcn, cnv_fice, &
212 cnv_ndrop, cnv_nice, cf_upi
213! *GJF
214 logical, intent(in) :: lprnt
215 integer, intent(in) :: ipr
216 integer, intent(inout) :: kcnv(:) ! zero if no deep convection and 1 otherwise
217 character(len=*), intent(out) :: errmsg
218 integer, intent(out) :: errflg
219
220!DDsigma - output added for AW sigma diagnostics
221! interface sigma and vertical velocity by cloud type (1=sfc)
222! real(kind_phys), intent(out), dimension(:,:,:) :: sigmai, vverti
223 real(kind_phys), intent(out), dimension(:,:) :: sigma ! sigma sigma totaled over cloud type - on interfaces (1=sfc)
224! sigma terms in eq 91 and 92
225! real(kind_phys), dimension(IJSDIM,KMAX) :: sfluxterm, qvfluxterm, condterm
226!DDsigma
227!
228! output arguments of CS_CUMLUS
229!
230 real(kind_phys), dimension(IJSDIM,KMAX+1,nctp) :: vverti, sigmai
231
232 real(kind_phys) gtt(ijsdim,kmax)
233 real(kind_phys) gtq(ijsdim,kmax,ntr)
234 real(kind_phys) gtu(ijsdim,kmax)
235 real(kind_phys) gtv(ijsdim,kmax)
236 real(kind_phys) cmdet(ijsdim,kmax)
237 real(kind_phys) gtprp(ijsdim,kmax+1)
238 real(kind_phys) gsnwp(ijsdim,kmax+1)
239 real(kind_phys) gmfx0(ijsdim,kmax+1)
240 real(kind_phys) gmfx1(ijsdim,kmax+1)
241 integer kt(ijsdim,nctp)
242
243 real(kind_phys) :: cape(ijsdim)
244 real(kind_phys) :: prec(ijsdim)
245 real(kind_phys) :: snow(ijsdim)
246!
247! input arguments of CS_CUMLUS
248!
249 real(kind_phys) gdt(ijsdim,kmax)
250 real(kind_phys) gdq(ijsdim,kmax,ntr)
251 real(kind_phys) gdu(ijsdim,kmax)
252 real(kind_phys) gdv(ijsdim,kmax)
253 real(kind_phys) gdtm(ijsdim,kmax+1)
254 real(kind_phys) gdp(ijsdim,kmax)
255 real(kind_phys) gdpm(ijsdim,kmax+1)
256 real(kind_phys) gdz(ijsdim,kmax)
257 real(kind_phys) gdzm(ijsdim,kmax+1)
258 real(kind_phys) delp(ijsdim,kmax)
259 real(kind_phys) delpi(ijsdim,kmax)
260!
261! local variables
262!
263!DD real(kind_phys) :: zs(IJSDIM) !< surface height [m]
264
265 integer ktmax(ijsdim)
266 real(kind_phys) :: ftintm, wrk, wrk1, tem
267 integer i, k, n, ists, iens, kp1
268
269!DD borrowed from RAS to go form total condensate to ice/water separately
270! parameter (tf=130.16, tcr=160.16, tcrf=1.0/(tcr-tf),tcl=2.0)
271! parameter (tf=230.16, tcr=260.16, tcrf=1.0/(tcr-tf))
272 real(kind_phys), parameter :: tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf), tcl=2.0
273 logical, save :: first=.true.
274
275 ! Initialize CCPP error handling variables
276 errmsg = ''
277 errflg = 0
278
279! lprnt = kdt == 1 .and. mype == 38
280! ipr = 43
281
282 precz0 = precz0in
283 preczh = preczhin
284 clmd = clmdin
285 clmp = (one-clmd)*(pa+pa)
286 clmdpa = clmd*pa
287!
288 if (first) then
289 do i=1,ntr
290 imfxr(i) = 0
291 enddo
292! IMFXR(1) = 1
293! IMFXR(ITL) = 1
294! IMFXR(ITI) = 1
295 first = .false.
296 endif
297!
298 ists = 1
299 iens = ijsdim
300
301 do k=1,kmax+1
302 do i=1,ijsdim
303 gdzm(i,k) = zi(i,k) * gravi
304 gdpm(i,k) = paph(i,k)
305 enddo
306 enddo
307
308 do k=1,kmax
309 do i=1,ijsdim
310 gdt(i,k) = t(i,k)
311 gdu(i,k) = u(i,k)
312 gdv(i,k) = v(i,k)
313 gdz(i,k) = zm(i,k) * gravi
314 gdp(i,k) = pap(i,k)
315 gdq(i,k,1) = q(i,k)
316 delp(i,k) = paph(i,k) - paph(i,k+1)
317 delpi(i,k) = grav / delp(i,k)
318 enddo
319 enddo
320
321!DD following adapted from ras
332 if (clw(1,1,2) <= -999.0) then ! input ice/water are together
333 do k=1,kmax
334 do i=1,ijsdim
335 tem = clw(i,k,1) * max(zero, min(one, (tcr-t(i,k))*tcrf))
336 clw(i,k,2) = clw(i,k,1) - tem
337 clw(i,k,1) = tem
338 enddo
339 enddo
340 endif
341!DD end ras adaptation
342 do k=1,kmax
343 do i=1,ijsdim
344 tem = min(clw(i,k,1), 0.0)
345 wrk = min(clw(i,k,2), 0.0)
346 clw(i,k,1) = clw(i,k,1) - tem
347 clw(i,k,2) = clw(i,k,2) - wrk
348 gdq(i,k,1) = gdq(i,k,1) + tem + wrk
349 enddo
350 enddo
351! if (lprnt) write(0,*)'in cs clw1b=',clw(ipr,:,1),' kdt=',kdt
352! if (lprnt) write(0,*)'in cs clw2b=',clw(ipr,:,2),' kdt=',kdt
353
354 do n=2,ntr
355 do k=1,kmax
356 do i=1,ijsdim
357 gdq(i,k,n) = clw(i,k,n-1)
358 enddo
359 enddo
360 enddo
361! if (lprnt) write(0,*)' incs tke=',gdq(ipr,1:25,ntr)
362!
363!***************************************************************************************
364!
366!
367
368 DO k=2,kmax
369 DO i=ists,iens
370 wrk = one / gdp(i,k)
371 wrk1 = one / log(gdp(i,k-1)*wrk)
372 ftintm = wrk1 * log(gdpm(i,k)*wrk)
373 gdtm(i,k) = ftintm*gdt(i,k-1) + (one-ftintm)*gdt(i,k)
374 ENDDO
375 ENDDO
376
377 DO i=ists,iens
378 gdtm(i,kmax+1) = gdt(i,kmax)
379 gdtm(i,1) = gdt(i,1) ! Is this a good approximation ? - Moorthi
380 ENDDO
381
383 do n=1,nctp
384 do k=1,kmax+1
385 do i=ists,iens
386 vverti(i,k,n) = zero
387 sigmai(i,k,n) = zero
388 enddo
389 enddo
390 enddo
391 do k=1,kmax+1
392 do i=ists,iens
393 sigma(i,k) = zero
394 enddo
395 enddo
396!
398 call cs_cumlus (ijsdim, ijsdim, kmax , ntr , & !DD dimensions
399 otspt(1:ntr,1), otspt(1:ntr,2), &
400 lprnt , ipr , &
401 gtt , gtq , gtu , gtv , & ! output
402 cmdet , & ! output
403 gtprp , gsnwp , gmfx0 , & ! output
404 gmfx1 , cape , kt , & ! output
405 cbmfx , & ! modified
406 gdt , gdq , gdu , gdv , & ! input
407 gdtm , & ! input
408 gdp , gdpm , gdz , gdzm , & ! input
409 delp , delpi , &
410 delta , delti , ists , iens, mype,& ! input
411 fscav, fswtr, wcbmaxm, nctp, &
412 sigmai, sigma, vverti, & ! input/output !DDsigma
413 do_aw, do_awdd, flx_form)
414!
415!
416!DD detrainment has to be added in for GFS
417!
418! if (lprnt) write(0,*)' aft cs_cum gtqi=',gtq(ipr,:,2)
419! if (lprnt) write(0,*)' aft cs_cum gtql=',gtq(ipr,:,3)
420
421
422 do n=2,ntr
423 do k=1,kmax
424 do i=1,ijsdim
425 clw(i,k,n-1) = max(zero, gdq(i,k,n) + gtq(i,k,n) * delta)
426 enddo
427 enddo
428 enddo
429! if (lprnt) write(0,*)' aftcs_cum tkein=',gdq(ipr,1:25,ntr),' delta=',delta
430! if (lprnt) write(0,*)' aftcs_cum tke=',clw(ipr,1:25,ntr-1)
431! if (lprnt) write(0,*)'in cs clw1a=',clw(ipr,:,1),' kdt=',kdt
432! if (lprnt) write(0,*)'in cs clw2a=',clw(ipr,:,2),' kdt=',kdt
433!
434 do k=1,kmax
435 do i=1,ijsdim
436 q(i,k) = max(zero, gdq(i,k,1) + gtq(i,k,1) * delta)
437 t(i,k) = gdt(i,k) + gtt(i,k) * delta
438 u(i,k) = gdu(i,k) + gtu(i,k) * delta
439 v(i,k) = gdv(i,k) + gtv(i,k) * delta
440! Set the mass fluxes.
441 ud_mf(i,k) = gmfx0(i,k)
442 dd_mf(i,k) = gmfx1(i,k)
443 dt_mf(i,k) = cmdet(i,k)
444 enddo
445 enddo
446
447
448 if (mp_phys == 10) then ! for 2M microphysics, always output these variables
449 if (do_aw) then
450 do k=1,kmax
451 kp1 = min(k+1,kmax)
452 do i=1,ijsdim
453 qicn(i,k) = max(0.0, clw(i,k,1)-gdq(i,k,2))
454 qlcn(i,k) = max(0.0, clw(i,k,2)-gdq(i,k,3))
455
456
457 wrk = qicn(i,k) + qlcn(i,k)
458 if (wrk > 1.0e-12) then
459 cnv_fice(i,k) = qicn(i,k) / wrk
460 else
461 cnv_fice(i,k) = 0.0
462 endif
463!
464! CNV_MFD(i,k) = dt_mf(i,k) * (1.0/delta)
465 cnv_mfd(i,k) = dt_mf(i,k)
466 cnv_dqldt(i,k) = wrk / delta
467! CNV_PRC3(i,k) = 0.0
468 cnv_ndrop(i,k) = 0.0
469 cnv_nice(i,k) = 0.0
470 cf_upi(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)),0.1))
471! CLCN(i,k) = cf_upi(i,k) !downdraft is below updraft
472!! clcn(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)/delta),0.25))
473
474 w_upi(i,k) = 0.0
475!! w_upi(i,k) = ud_mf(i,k)*(t(i,k)+epsvt*gdq(i,k,1)) * rair &
476!! / (delta*max(cf_upi(i,k),1.e-12)*gdp(i,k))
477 enddo
478 enddo
479 do k=1,kmax
480 do i=1,ijsdim
481 do n=1,nctp
482 w_upi(i,k) = w_upi(i,k) + vverti(i,k,n)
483 enddo
484 if (sigma(i,k) > 1.0e-10) then
485 w_upi(i,k) = w_upi(i,k) / sigma(i,k)
486 else
487 w_upi(i,k) = 0.0
488 endif
489 enddo
490 enddo
491 else
492 do k=1,kmax
493 do i=1,ijsdim
494 qicn(i,k) = max(0.0, clw(i,k,1)-gdq(i,k,2))
495 qlcn(i,k) = max(0.0, clw(i,k,2)-gdq(i,k,3))
496 cnv_fice(i,k) = qicn(i,k) / max(1.0e-10,qicn(i,k)+qlcn(i,k))
497!
498! CNV_MFD(i,k) = dt_mf(i,k) * (1/delta)
499 cnv_mfd(i,k) = dt_mf(i,k)
500 cnv_dqldt(i,k) = (qicn(i,k)+qlcn(i,k)) / delta
501! CNV_PRC3(i,k) = 0.0
502 cnv_ndrop(i,k) = 0.0
503 cnv_nice(i,k) = 0.0
504 cf_upi(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)),0.1))
505! & 500*ud_mf(i,k)),0.60))
506! CLCN(i,k) = cf_upi(i,k) !downdraft is below updraft
507
508 w_upi(i,k) = ud_mf(i,k)*(t(i,k)+epsvt*gdq(i,k,1)) * rair &
509 / (max(cf_upi(i,k),1.e-12)*gdp(i,k))
510 enddo
511 enddo
512 endif
513 endif
514
515!****************************************************************************
516
517 ktmax = 1
518 do n=1,nctp
519 do i=1,ijsdim
520 ktmax(i) = max(ktmax(i), kt(i,n))
521 enddo
522 enddo
523!
524 do i=1,ijsdim
525 prec(i) = gtprp(i,1)
526 snow(i) = gsnwp(i,1)
527 if (prec(i)+snow(i) > 0.0) then
528 kcnv(i) = 1
529 else
530 kcnv(i) = 0
531 endif
532 enddo
533
535
536 do k=1,kmax
537 do i=1,ijsdim
538 ud_mf(i,k) = ud_mf(i,k) * delta
539 dd_mf(i,k) = dd_mf(i,k) * delta
540 dt_mf(i,k) = dt_mf(i,k) * delta
541 enddo
542 enddo
543
544! rain1(:) = prec(:) * (delta*0.001) ! Convert prec flux kg/m2/sec to rain1 in m
545 do i= 1, ijsdim
546 rain1(i) = prec(i) * (delta*0.001)
547 enddo
548
549! if (lprnt) then
550! write(0,*)' aft cs_cum prec=',prec(ipr),'GTPRP=',GTPRP(ipr,1)
551! endif
552
553
554! if (do_aw) then
555! call moist_bud(ijsdim,ijsdim,im,kmax,mype,kdt,grav,delta,delp,prec &
556! , gdq(1,1,1), gdq(1,1,2), gdq(1,1,3) &
557! , q,clw(1,1,1),clw(1,1,2),'cs_conv_aw')
558! endif
559
560 end subroutine cs_conv_run
562
563
564!************************************************************************
588 SUBROUTINE cs_cumlus (im , IJSDIM, KMAX , NTR , & !DD dimensions
589 otspt1, otspt2, lprnt , ipr , &
590 GTT , GTQ , GTU , GTV , & ! output
591 CMDET , & ! output
592 GTPRP , GSNWP , GMFX0 , & ! output
593 GMFX1 , CAPE , KT , & ! output
594 CBMFX , & ! modified
595 GDT , GDQ , GDU , GDV , & ! input
596 GDTM , & ! input
597 GDP , GDPM , GDZ , GDZM , & ! input
598 delp , delpinv , &
599 DELTA , DELTI , ISTS , IENS, mype,& ! input
600 fscav, fswtr, wcbmaxm, nctp, & !
601 sigmai, sigma, vverti, & ! input/output !DDsigma
602 do_aw, do_awdd, flx_form)
603!
604 IMPLICIT NONE
605
606 Integer, parameter :: ntrq=4 ! starting index for tracers
607 INTEGER, INTENT(IN) :: im, IJSDIM, KMAX, NTR, mype, nctp, ipr !! DD, for GFS, pass in
608 logical, intent(in) :: do_aw, do_awdd, flx_form ! switch to apply Arakawa-Wu to the tendencies
609 logical, intent(in) :: otspt1(ntr), otspt2(ntr), lprnt
610 REAL(kind_phys),intent(in) :: DELP (IJSDIM, KMAX)
611 REAL(kind_phys),intent(in) :: DELPINV (IJSDIM, KMAX)
612!
613! [OUTPUT]
614 REAL(kind_phys), INTENT(OUT) :: GTT (IJSDIM, KMAX ) ! heating rate
615 REAL(kind_phys), INTENT(OUT) :: GTQ (IJSDIM, KMAX, NTR) ! change in q
616 REAL(kind_phys), INTENT(OUT) :: GTU (IJSDIM, KMAX ) ! tendency of u
617 REAL(kind_phys), INTENT(OUT) :: GTV (IJSDIM, KMAX ) ! tendency of v
618 REAL(kind_phys), INTENT(OUT) :: CMDET (IJSDIM, KMAX ) ! detrainment mass flux
619 REAL(kind_phys) :: GTLDET( IJSDIM, KMAX ) ! cloud liquid tendency by detrainment
620 REAL(kind_phys) :: GTIDET( IJSDIM, KMAX ) ! cloud ice tendency by detrainment
621! assuming there is no flux at the top of the atmospherea - Moorthi
622 REAL(kind_phys), INTENT(OUT) :: GTPRP (IJSDIM, KMAX+1 ) ! rain+snow flux
623 REAL(kind_phys), INTENT(OUT) :: GSNWP (IJSDIM, KMAX+1 ) ! snowfall flux
624 REAL(kind_phys), INTENT(OUT) :: GMFX0 (IJSDIM, KMAX+1 ) ! updraft mass flux
625 REAL(kind_phys), INTENT(OUT) :: GMFX1 (IJSDIM, KMAX+1 ) ! downdraft mass flux
626
627 REAL(kind_phys), INTENT(OUT) :: CAPE (IJSDIM )
628 INTEGER , INTENT(OUT) :: KT (IJSDIM, NCTP ) ! cloud top
629!
630! [MODIFIED]
631 REAL(kind_phys), INTENT(INOUT) :: CBMFX ( IM, NCTP ) !! cloud base mass flux
632
633 !DDsigma - output added for AW sigma diagnostics
634 real(kind_phys), intent(out) :: sigmai(IM,KMAX+1,nctp) !DDsigma sigma by cloud type - on interfaces (1=sfc)
635 real(kind_phys), intent(out) :: vverti(IM,KMAX+1,nctp) !DDsigma vert. vel. by cloud type - on interfaces (1=sfc)
636 real(kind_phys), intent(out) :: sigma(IM,KMAX+1) !DDsigma sigma totaled over cloud type - on interfaces (1=sfc)
637
638! for computing AW flux form of tendencies
639! real(kind_phys), dimension(IM,KMAX) :: & !DDsigmadiag
640! sfluxterm, qvfluxterm
641! real(kind_phys), dimension(IM,KMAX) :: & !DDsigmadiag
642! qlfluxterm, qifluxterm
643! real(kind_phys), dimension(ijsdim,kmax,ntrq:ntr) :: trfluxterm ! tendencies of tracers due to eddy mass flux
644 real(kind_phys), dimension(IM,KMAX) :: & !DDsigmadiag
645 condtermt, condtermq, frzterm, prectermq, prectermfrz
646 !DDsigma
647
648!
649! [INPUT]
650 REAL(kind_phys), INTENT(IN) :: GDT (IJSDIM, KMAX ) ! temperature T
651 REAL(kind_phys), INTENT(IN) :: GDQ (IJSDIM, KMAX, NTR) ! humidity, tracer !DDsigmadiag
652 REAL(kind_phys), INTENT(IN) :: GDU (IJSDIM, KMAX ) ! westerly u
653 REAL(kind_phys), INTENT(IN) :: GDV (IJSDIM, KMAX ) ! southern wind v
654 REAL(kind_phys), INTENT(IN) :: GDTM (IJSDIM, KMAX+1 ) ! temperature T
655 REAL(kind_phys), INTENT(IN) :: GDP (IJSDIM, KMAX ) ! pressure P
656 REAL(kind_phys), INTENT(IN) :: GDPM (IJSDIM, KMAX+1 ) ! pressure (half lev)
657 REAL(kind_phys), INTENT(IN) :: GDZ (IJSDIM, KMAX ) ! altitude
658 REAL(kind_phys), INTENT(IN) :: GDZM (IJSDIM, KMAX+1 ) ! altitude
659 REAL(kind_phys), INTENT(IN) :: DELTA ! delta(t) (dynamics)
660 REAL(kind_phys), INTENT(IN) :: DELTI ! delta(t) (internal variable)
661 INTEGER, INTENT(IN) :: ISTS, IENS ! array range
662
663 real(kind_phys), intent(in) :: fscav(ntr), fswtr(ntr), wcbmaxm(ijsdim)
664!
665! [INTERNAL WORK]
666 REAL(kind_phys), allocatable :: GPRCC (:, :) ! rainfall
667 REAL(kind_phys) GSNWC ( IJSDIM ) ! snowfall
668 REAL(kind_phys) CUMCLW( IJSDIM, KMAX ) ! cloud water in cumulus
669 REAL(kind_phys) CUMFRC( IJSDIM ) ! cumulus cloud fraction
670!COSP
671 REAL(kind_phys) QLIQC ( IJSDIM, KMAX ) ! cumulus cloud liquid water [kg/kg]
672 REAL(kind_phys) QICEC ( IJSDIM, KMAX ) ! cumulus cloud ice [kg/kg]
673 REAL(kind_phys) GPRCPF( IJSDIM, KMAX ) ! rainfall flux at full level
674 REAL(kind_phys) GSNWPF( IJSDIM, KMAX ) ! snowfall flux at full level
675!
676 REAL(kind_phys) GTCFRC( IJSDIM, KMAX ) ! change in cloud fraction
677 REAL(kind_phys) FLIQC ( IJSDIM, KMAX ) ! liquid ratio in cumulus
678!
679!#ifdef OPT_CHASER
680! REAL(kind_phys) RFXC ( IJSDIM, KMAX+1 ) ! precipi. flx [kg/m2/s]
681! REAL(kind_phys) SFXC ( IJSDIM, KMAX+1 ) ! ice/snow flx [kg/m2/s]
682! INTEGER LEVCUM( IJSDIM, KMAX ) ! flag for cum. cloud top
683! REAL(kind_phys) LNFRC ( IJSDIM, KMAX ) ! areal rates of clouds
684! REAL(kind_phys) REVC ( IJSDIM, KMAX ) ! evaporation rates
685!#endif
686!
687 REAL(kind_phys) GDCFRC( IJSDIM, KMAX ) ! cloud fraction
688!
689! REAL(kind_phys) GTQL ( IJSDIM, KMAX ) ! tendency of cloud liquid
690!
691 REAL(kind_phys) GDW ( IJSDIM, KMAX ) ! total water
692 REAL(kind_phys) GDQS ( IJSDIM, KMAX ) ! saturate moisture
693 REAL(kind_phys) FDQS ( IJSDIM, KMAX )
694 REAL(kind_phys) GAM ( IJSDIM, KMAX )
695 REAL(kind_phys) GDS ( IJSDIM, KMAX ) ! dry static energy
696 REAL(kind_phys) GDH ( IJSDIM, KMAX ) ! moist static energy
697 REAL(kind_phys) GDHS ( IJSDIM, KMAX ) ! saturate MSE
698!
699 REAL(kind_phys) GCYM ( IJSDIM, KMAX, NCTP ) ! norm. mass flux (half lev)
700 REAL(kind_phys) GCHB ( IJSDIM ) ! cloud base MSE-Li*Qi
701 REAL(kind_phys) GCWB ( IJSDIM ) ! cloud base total water
702 REAL(kind_phys) GCtrB ( IJSDIM, ntrq:ntr ) ! cloud base water vapor tracer
703 REAL(kind_phys) GCUB ( IJSDIM ) ! cloud base U
704 REAL(kind_phys) GCVB ( IJSDIM ) ! cloud base V
705 REAL(kind_phys) GCIB ( IJSDIM ) ! cloud base ice
706 REAL(kind_phys) ELAM ( IJSDIM, KMAX, NCTP ) ! entrainment (rate*massflux)
707 REAL(kind_phys) GCYT ( IJSDIM, NCTP ) ! norm. mass flux @top
708 REAL(kind_phys) GCHT ( IJSDIM, NCTP ) ! cloud top MSE
709 REAL(kind_phys) GCQT ( IJSDIM, NCTP ) ! cloud top q
710 REAL(kind_phys) GCwT ( IJSDIM ) ! cloud top total water
711 REAL(kind_phys) GCUT ( IJSDIM, NCTP ) ! cloud top U
712 REAL(kind_phys) GCVT ( IJSDIM, NCTP ) ! cloud top V
713 REAL(kind_phys) GCLT ( IJSDIM, NCTP ) ! cloud top cloud water
714 REAL(kind_phys) GCIT ( IJSDIM, NCTP ) ! cloud top cloud ice
715 REAL(kind_phys) GCtrT (IJSDIM, ntrq:ntr, NCTP) ! cloud top tracer
716 REAL(kind_phys) GTPRT ( IJSDIM, NCTP ) ! precipitation/M
717 REAL(kind_phys) GCLZ ( IJSDIM, KMAX ) ! cloud liquid for each CTP
718 REAL(kind_phys) GCIZ ( IJSDIM, KMAX ) ! cloud ice for each CTP
719
720 REAL(kind_phys) ACWF ( IJSDIM ) ! cloud work function
721 REAL(kind_phys) GPRCIZ( IJSDIM, KMAX+1, NCTP ) ! precipitation
722 REAL(kind_phys) GSNWIZ( IJSDIM, KMAX+1, NCTP ) ! snowfall
723 REAL(kind_phys) GTPRC0( IJSDIM ) ! precip. before evap.
724
725 REAL(kind_phys) GMFLX ( IJSDIM, KMAX+1 ) ! mass flux (updraft+downdraft)
726 REAL(kind_phys) QLIQ ( IJSDIM, KMAX ) ! total cloud liquid
727 REAL(kind_phys) QICE ( IJSDIM, KMAX ) ! total cloud ice
728 REAL(kind_phys) GPRCI ( IJSDIM, KMAX ) ! rainfall generation
729 REAL(kind_phys) GSNWI ( IJSDIM, KMAX ) ! snowfall generation
730
731 REAL(kind_phys) GPRCP ( IJSDIM, KMAX+1 ) ! rainfall flux
732!
733 REAL(kind_phys) GTEVP ( IJSDIM, KMAX ) ! evaporation+sublimation
734 REAL(kind_phys) GMDD ( IJSDIM, KMAX+1 ) ! downdraft mass flux
735
736 REAL(kind_phys) CUMHGT( IJSDIM, NCTP ) ! cloud top height
737 REAL(kind_phys) CTOPP ( IJSDIM ) ! cloud top pressure
738
739 REAL(kind_phys) GDZTR ( IJSDIM ) ! tropopause height
740 REAL(kind_phys) FLIQOU( IJSDIM, KMAX ) ! liquid ratio in cumulus
741!#ifdef OPT_CHASER
742! REAL(kind_phys) TOPFLX( IJSDIM, NCTP ) !! flux at each cloud top
743!#endif
744 INTEGER KB ( IJSDIM )
745 INTEGER KSTRT ( IJSDIM ) ! tropopause level
746 REAL(kind_phys) GAMX
747 REAL(kind_phys) CIN ( IJSDIM )
748 INTEGER JBUOY ( IJSDIM )
749 REAL(kind_phys) DELZ, BUOY, DELWC, DELER
750!M REAL(kind_phys) WCB ( NCTP ) ! updraft velocity**2 @base
751!M SAVE WCB
752 REAL(kind_phys) WCBX (IJSDIM)
753! REAL(kind_phys) ERMR ( NCTP ) ! entrainment rate (ASMODE)
754! SAVE ERMR
755 INTEGER KTMX ( NCTP ) ! max of cloud top
756 INTEGER KTMXT ! max of cloud top
757 REAL(kind_phys) TIMED
758 REAL(kind_phys) GDCLDX, GDMU2X, GDMU3X
759!
760 LOGICAL OOUT1, OOUT2
761 INTEGER KBMX, I, K, CTP, ierr, n, kp1, l, l1, kk, kbi, kmi, km1
762 real(kind_phys) tem1, tem2, tem3, cbmfl, mflx_e, teme, tems
763
764 REAL(kind_phys) HBGT ( IJSDIM ) ! imbalance in column heat
765 REAL(kind_phys) WBGT ( IJSDIM ) ! imbalance in column water
766
767 !DDsigma begin local work variables - all on model interfaces (sfc=1)
768 REAL(kind_phys) lamdai( IJSDIM, KMAX+1, nctp ) ! lamda for cloud type ctp
769 REAL(kind_phys) lamdaprod( IJSDIM, KMAX+1 ) ! product of (1+lamda) through cloud type ctp
770 REAL(kind_phys) gdrhom ! density
771 REAL(kind_phys) gdtvm ! virtual temperature
772 REAL(kind_phys) gdqm, gdwm,gdlm, gdim ! water vaper
773 REAL(kind_phys) gdtrm(ntrq:ntr) ! tracer
774 character(len=4) :: cproc !DDsigmadiag
775
776 ! the following are new arguments to cumup to get them out
777 REAL(kind_phys) wcv( IJSDIM, KMAX+1, nctp) ! in-cloud vertical velocity
778 REAL(kind_phys) GCTM ( IJSDIM, KMAX+1 ) ! cloud T (half lev) !DDsigmadiag make output
779 REAL(kind_phys) GCQM ( IJSDIM, KMAX+1, nctp ) ! cloud q (half lev) !DDsigmadiag make output
780 REAL(kind_phys) GCwM ( IJSDIM, KMAX+1, nctp ) ! cloud q (half lev) !DDsigmadiag make output
781 REAL(kind_phys) GCiM ( IJSDIM, KMAX+1 ) ! cloud q (half lev) !DDsigmadiag make output
782 REAL(kind_phys) GClM ( IJSDIM, KMAX+1 ) ! cloud q (half lev) !DDsigmadiag make output
783 REAL(kind_phys) GChM ( IJSDIM, KMAX+1, nctp ) ! cloud q (half lev) !DDsigmadiag make output
784 REAL(kind_phys) GCtrM (IJSDIM, KMAX, ntrq:ntr) ! cloud tracer (half lev) !DDsigmadiag make output
785
786! these are the fluxes at the interfaces - AW will operate on them
787 REAL(kind_phys), dimension(ijsdim,Kmax+1,nctp) :: sfluxtem, qvfluxtem, qlfluxtem, qifluxtem
788 REAL(kind_phys), dimension(ijsdim,Kmax+1,ntrq:ntr,nctp) :: trfluxtem ! tracer
789
790 REAL(kind_phys), dimension(ijsdim,Kmax+1) :: dtcondtem, dqcondtem, dtfrztem, dqprectem,dfrzprectem
791 REAL(kind_phys), dimension(ijsdim,Kmax) :: dtevap, dqevap, dtmelt, dtsubl
792 REAL(kind_phys), dimension(ijsdim) :: moistening_aw
793 real(kind_phys) rhs_q, rhs_h, sftem, qftem, qlftem, qiftem
794 real(kind_phys), dimension(ijsdim,kmax+1) :: gctbl, gcqbl,gcwbl, gcqlbl, gcqibl !DDsigmadiag updraft profiles below cloud Base
795 real(kind_phys), dimension(ijsdim,kmax,ntrq:ntr) :: gctrbl !DDsigmadiag tracer updraft profiles below cloud Base
796 real(kind_phys), dimension(ijsdim,kmax+1) :: sigmad
797 real(kind_phys) :: fsigma( IJSDIM, KMAX+1 ) ! factor to reduce mass flux terms (1-sigma**2) for AW
798 real(kind_phys) :: lamdamax ! for sorting lamda values
799 integer loclamdamax
800 real(kind_phys) :: pr_tot, pr_ice, pr_liq
801!DDsigma end local work variables
802!
803! [INTERNAL PARM]
804 REAL(kind_phys) :: WCBMIN = 0._kind_phys ! min. of updraft velocity at cloud base
805!M REAL(kind_phys) :: WCBMAX = 1.4_kind_phys ! max. of updraft velocity at cloud base
806!M wcbas commented by Moorthi since it is not used
807!M REAL(kind_phys) :: WCBAS = 2._kind_phys ! updraft velocity**2 at cloud base (ASMODE)
808!M REAL(kind_phys) :: ERAMIN = 1.e-5_kind_phys ! min. of entrainment rate
809 ! used only in OPT_ASMODE
810!M REAL(kind_phys) :: ERAMAX = 2.e-3_kind_phys ! max. of entrainment rate
811 ! used only in OPT_ASMODE
812! downdraft mass flux terms now slot nctp+1 in the *fluxterm arrays
813 REAL(kind_phys) dtdwn ( IJSDIM, KMAX ) ! t tendency downdraft detrainment
814 REAL(kind_phys) dqvdwn ( IJSDIM, KMAX ) ! qv tendency downdraft detrainment
815 REAL(kind_phys) dqldwn ( IJSDIM, KMAX ) ! ql tendency downdraft detrainment
816 REAL(kind_phys) dqidwn ( IJSDIM, KMAX ) ! qi tendency downdraft detrainment
817 REAL(kind_phys), dimension(ijsdim,kmax,ntrq:ntr) :: dtrdwn ! tracer tendency downdraft detrainment
818
819 LOGICAL :: OINICB = .false. ! set 0.d0 to CBMFX
820
821 REAL(kind_phys) :: VARMIN = 1.e-13_kind_phys ! minimum of PDF variance
822 REAL(kind_phys) :: VARMAX = 5.e-7_kind_phys ! maximum of PDF variance
823 REAL(kind_phys) :: SKWMAX = 0.566_kind_phys ! maximum of PDF skewness
824
825 REAL(kind_phys) :: PSTRMX = 400.e2_kind_phys ! max P of tropopause
826 REAL(kind_phys) :: PSTRMN = 50.e2_kind_phys ! min P of tropopause
827 REAL(kind_phys) :: GCRSTR = 1.e-4_kind_phys ! crit. dT/dz tropopause
828
829 ! 0: mass fixer is not applied
830 ! tracers which may become negative values
831 ! e.g. subgrid-PDFs
832 ! 1: mass fixer is applied, total mass may change through cumulus scheme
833 ! e.g. moisture, liquid cloud, ice cloud, aerosols
834 ! 2: mass fixer is applied, total mass never change through cumulus scheme
835 ! e.g. CO2
836 real(kind=kind_phys), parameter :: zero=0.0, one=1.0
837 real(kind=kind_phys) :: tem, esat
838!
839 LOGICAL, SAVE :: OFIRST = .true. ! called first time?
840!
841 IF (ofirst) THEN
842 ofirst = .false.
843 IF (oinicb) THEN
844 cbmfx = zero
845 ENDIF
846 ENDIF
847!
848
849 kp1 = kmax + 1
850 do n=1,ntr
851 do k=1,kmax
852 do i=1,ijsdim
853 gtq(i,k,n) = zero
854 enddo
855 enddo
856 enddo
857 do k=1,kmax+1
858 do i=1,ijsdim
859 gmflx(i,k) = zero
860 gmfx0(i,k) = zero
861 enddo
862 enddo
863 do k=1,kmax
864 do i=1,ijsdim
865 gtt(i,k) = zero
866 gtu(i,k) = zero
867 gtv(i,k) = zero
868 gprci(i,k) = zero
869 gsnwi(i,k) = zero
870 qliq(i,k) = zero
871 qice(i,k) = zero
872! gtcfrc(i,k) = zero
873! cumclw(i,k) = zero
874! fliqc(i,k) = zero
875 fliqou(i,k) = zero
876 gprcpf(i,k) = zero
877 gsnwpf(i,k) = zero
878 cmdet(i,k) = zero
879 enddo
880 enddo
881 if (flx_form) then
882 do ctp = 1,nctp
883 do k=1,kp1
884 do i=1,ijsdim
885 sfluxtem(i,k,ctp) = zero
886 qvfluxtem(i,k,ctp) = zero
887 qlfluxtem(i,k,ctp) = zero
888 qifluxtem(i,k,ctp) = zero
889 enddo
890 enddo
891 do n = ntrq,ntr
892 do k=1,kp1
893 do i=1,ijsdim
894 trfluxtem(i,k,n,ctp) = zero
895 enddo
896 enddo
897 enddo
898 enddo
899 do k=1,kmax
900 do i=1,ijsdim
901 condtermt(i,k) = zero
902 condtermq(i,k) = zero
903 frzterm(i,k) = zero
904 prectermq(i,k) = zero
905 prectermfrz(i,k) = zero
906 enddo
907 enddo
908 do k=1,kmax
909 do i=1,ijsdim
910 dtdwn(i,k) = zero
911 dqvdwn(i,k) = zero
912 dqldwn(i,k) = zero
913 dqidwn(i,k) = zero
914 enddo
915 enddo
916 do n = ntrq,ntr
917 do k=1,kmax
918 do i=1,ijsdim
919 dtrdwn(i,k,n) = zero
920 enddo
921 enddo
922 enddo
923 endif
924 do i=1,ijsdim
925! gprcc(i,:) = zero
926! gmflx(i,kp1) = zero
927 gmfx0(i,kp1) = zero
928 gtprc0(i) = zero
929! hbgt(i) = zero
930! wbgt(i) = zero
931 gdztr(i) = zero
932 kstrt(i) = kmax
933 enddo
934
935 do k=1,kmax
936 do i=1,ijsdim
937 gdw(i,k) = gdq(i,k,1) + gdq(i,k,itl) + gdq(i,k,iti)
938 enddo
939 enddo
943 DO k=1,kmax
944 DO i=ists,iens
945 esat = min(gdp(i,k), fpvs(gdt(i,k)))
946 gdqs(i,k) = min(epsv*esat/max(gdp(i,k)+epsvm1*esat, 1.0e-10), 0.1)
947 tem = one / gdt(i,k)
948 fdqs(i,k) = gdqs(i,k) * tem * (fact1 + fact2*tem) ! calculate d(qs)/dT
949 gam(i,k) = elocp*fdqs(i,k)
950 gds(i,k) = cp*gdt(i,k) + grav*gdz(i,k) ! layer dry static energy
951 gdh(i,k) = gds(i,k) + el*gdq(i,k,1) ! layer moist static energy
952 gdhs(i,k) = gds(i,k) + el*gdqs(i,k) ! layer sat. moist static energy
953 ENDDO
954 ENDDO
955!
956! < tropopause >
957!
959 DO k=1,kmax
960 DO i=ists,iens
961 gamx = (gdtm(i,k+1)-gdtm(i,k)) / (gdzm(i,k+1)-gdzm(i,k))
962 IF ((gdp(i,k) < pstrmx .AND. gamx > gcrstr) .OR. gdp(i,k) < pstrmn) THEN
963 kstrt(i) = min(k, kstrt(i))
964 ENDIF
965 ENDDO
966 ENDDO
967 DO i=ists,iens
968 k = kstrt(i)
969 gdztr(i) = gdzm(i,k)
970 ENDDO
971!
972!DDsigma - arguments added to get subcloud profiles in updraft
973! so AW eddy flux tendencies can be computed
974
976 CALL cumbas(ijsdim, kmax , & !DD dimensions
977 kb , gcym(:,:,1) , kbmx , & ! output
978 ntr , ntrq , &
979 gchb , gcwb , gcub , gcvb , & ! output
980 gcib , gctrb, & ! output
981 gdh , gdw , gdhs , gdqs , & ! input
982 gdq(:,:,iti) , gdu , gdv , gdzm , & ! input
983 gdpm , fdqs , gam , & ! input
984 lprnt, ipr, &
985 ists , iens , & !) ! input
986 gctbl, gcqbl,gdq,gcwbl, gcqlbl, gcqibl, gctrbl) ! sub cloud tendencies
987!
989!
990 DO i=ists,iens
991 cape(i) = zero
992 cin(i) = zero
993 jbuoy(i) = 0
994 enddo
995 DO k=2,kmax
996 DO i=ists,iens
997 if (kb(i) > 0) then
998 IF (k >= kb(i)) THEN
999 buoy = (gdh(i,1)-gdhs(i,k)) / ((one+elocp*fdqs(i,k)) * cp*gdt(i,k))
1000 ELSE
1001 buoy = (gds(i,1)-gds(i,k)) / (cp*gdt(i,k))
1002 END IF
1003 IF (buoy > zero .AND. jbuoy(i) >= -1) THEN
1004 cape(i) = cape(i) + buoy * grav * (gdzm(i,k+1) - gdzm(i,k))
1005 jbuoy(i) = 2
1006 ELSEIF (buoy < zero .AND. jbuoy(i) /= 2) THEN
1007 cin(i) = cin(i) + buoy * grav * (gdzm(i,k+1) - gdzm(i,k))
1008 jbuoy(i) = -1
1009 ENDIF
1010 endif
1011 ENDDO
1012 ENDDO
1013 DO i=ists,iens
1014 IF (jbuoy(i) /= 2) cin(i) = -999.d0
1015 if (cin(i) < cincrit) kb(i) = -1
1016 ENDDO
1017
1018!DDsigma some initialization before summing over cloud type
1020 if(flx_form) then
1021 do k=1,kp1 ! Moorthi
1022 do i=1,ijsdim
1023 lamdaprod(i,k) = one
1024 sigma(i,k) = 0.0
1025 enddo
1026 enddo
1027
1028 do ctp=1,nctp
1029 do k=1,kp1
1030 do i=1,ijsdim
1031 lamdai(i,k,ctp) = zero
1032 sigmai(i,k,ctp) = zero
1033 vverti(i,k,ctp) = zero
1034 enddo
1035 enddo
1036 enddo
1037 endif
1038
1039 do ctp=2,nctp
1040 do k=1,kmax
1041 do i=1,ijsdim
1042 gcym(i,k,ctp) = gcym(i,k,1)
1043 enddo
1044 enddo
1045 enddo
1046
1047 DO ctp=1,nctp ! loop over cloud types
1048
1049 tem = ctp / dble(nctp)
1050 do i=1,ijsdim
1051 delwc = tem * (wcbmaxm(i) - wcbmin)
1052 wcbx(i) = delwc * delwc
1053 enddo
1054
1055! getting more incloud profiles of variables to compute eddy flux tendencies
1056! and condensation rates
1057
1058
1059! DH* GNU crashes - check all arguments to CUMUP for their dimensions
1060! before and after CUMUP (i.e. here), and inside the routine, in
1061! particular: gctm, gcqm, gcwm, gchm, gcwt, gclm, gcim,gctrm
1062! also, inside, check that no reads/writes out of bounds occur *DH
1064 CALL cumup(ijsdim, kmax, ntr, ntrq, & !DD dimensions
1065 acwf , & ! output
1066 gclz , gciz , gprciz(:,:,ctp), gsnwiz(:,:,ctp), & ! output
1067 gcyt(:,ctp) , gcht(:,ctp) , gcqt(:,ctp), & ! output
1068 gclt(:,ctp) , gcit(:,ctp) , gtprt(:,ctp), & ! output
1069 gcut(:,ctp) , gcvt(:,ctp) , gctrt(:,ntrq:ntr,ctp), & ! output
1070 kt(:,ctp) , ktmx(ctp) , & ! output
1071 gcym(:,:,ctp) , & ! modified
1072 wcv(:,:,ctp) , & ! !DD-sigma new output
1073 gchb , gcwb , gcub , gcvb , & ! input !DDsigmadiag
1074 gcib , gctrb , & ! input
1075 gdu , gdv , gdh , gdw , & ! input
1076 gdhs , gdqs , gdt , gdtm , & ! input
1077 gdq , gdq(:,:,iti) , gdz , gdzm , & ! input
1078 gdpm , fdqs , gam , gdztr , & ! input
1079 cpres , wcbx , & ! input
1080 kb , ctp , ists , iens , & ! input
1081 gctm , gcqm(:,:,ctp), gcwm(:,:,ctp), gchm(:,:,ctp),&
1082 gcwt, gclm, gcim, gctrm, & ! additional incloud profiles and cloud top total water
1083 lprnt , ipr )
1084!
1086 CALL cumbmx(ijsdim, kmax, & !DD dimensions
1087 cbmfx(:,ctp), & ! modified
1088 acwf , gcyt(:,ctp), gdzm , & ! input
1089 gdw , gdqs , delp , & ! input
1090 kt(:,ctp), ktmx(ctp) , kb , & ! input
1091 delti , ists , iens )
1092
1093!DDsigma - begin sigma computation
1094! At this point cbmfx is updated and we have everything we need to compute sigma
1095
1096 if (flx_form) then
1097 do k=1,kmax + 1 ! Moorthi
1098 do i=1,ijsdim
1099 dqcondtem(i,k) = zero
1100 dqprectem(i,k) = zero
1101 dfrzprectem(i,k) = zero
1102 dtfrztem(i,k) = zero
1103 dtcondtem(i,k) = zero
1104 enddo
1105 enddo
1106
1107 do i=ists,iens
1108 cbmfl = cbmfx(i,ctp)
1109 kk = kt(i,ctp) ! cloud top index
1110
1111 if(cbmfl > zero) then ! this should avoid zero wcv in the denominator
1112 kbi = kb(i) ! cloud base index
1113 do k=kbi,kk ! loop from cloud base to cloud top
1114 km1 = k - 1
1115! get environment variables interpolated to layer interface
1116 gdqm = half * (gdq(i,k,1) + gdq(i,km1,1)) ! as computed in cumup
1117! GDwM = half * (GDw(I,K) + GDw(I,KM1 ))
1118 gdlm = half * (gdq(i,k,itl) + gdq(i,km1,itl))
1119 gdim = half * (gdq(i,k,iti) + gdq(i,km1,iti))
1120 do n = ntrq,ntr
1121 gdtrm(n) = half * (gdq(i,k,n) + gdq(i,km1,n)) ! as computed in cumup
1122 enddo
1123 mflx_e = gcym(i,k,ctp) * cbmfl ! mass flux at level k for cloud ctp
1124
1125
1129
1130 lamdai(i,k,ctp) = mflx_e * rair * gdtm(i,k)*(one+epsvt*gdqm) &
1131 / (gdpm(i,k)*wcv(i,k,ctp))
1132
1133! just compute lamdai here, we will compute sigma, sigmai, and vverti outside
1134! the cloud type loop after we can sort lamdai
1135! lamdaprod(i,k) = lamdaprod(i,k) * (one+lamdai(i,k,ctp))
1136!
1137!! vverti(i,k,ctp) = wcv(i,k)
1138!! sigmai(i,k,ctp) = lamdai / lamdaprod(i,k)
1139!! sigma(i,k) = max(zero, min(one, sigma(i,k) + sigmai(i,k,ctp)))
1140!
1141! sigmai(i,k,ctp) = lamdai(i,k,ctp) / lamdaprod(i,k)
1142! sigma(i,k) = max(zero, min(one, sigma(i,k) + sigmai(i,k,ctp)))
1143! vverti(i,k,ctp) = sigmai(i,k,ctp) * wcv(i,k,ctp)
1144
1145
1146! sigma effect won't be applied until later, when lamda is sorted
1147! fsigma = 1.0 ! no aw effect, comment following lines to undo AW
1148! fsigma = one - sigma(i,k)
1149
1151! fsigma is the AW reduction of flux tendencies
1152
1153 if(k == kbi) then
1154 do l=2,kbi ! compute eddy fluxes below cloud base
1155! tem = - fsigma * gcym(i,l,ctp) * cbmfl
1156 tem = - gcym(i,l,ctp) * cbmfl
1157
1158! first get environment variables at layer interface
1159 l1 = l - 1
1160 gdqm = half * (gdq(i,l,1) + gdq(i,l1,1))
1161 gdlm = half * (gdq(i,l,itl) + gdq(i,l1,itl))
1162 gdim = half * (gdq(i,l,iti) + gdq(i,l1,iti))
1163!! GDwM = half * (GDw(I,l) + GDw(I,l1))
1164 do n = ntrq,ntr
1165 gdtrm(n) = half * (gdq(i,l,n) + gdq(i,l1,n)) ! as computed in cumup
1166 enddo
1167
1168! flux = mass flux * (updraft variable minus environment variable)
1169!centered differences
1170 sfluxtem(i,l,ctp) = tem * (gdtm(i,l)-gctbl(i,l))
1171 qvfluxtem(i,l,ctp) = tem * (gdqm-gcqbl(i,l))
1172 qlfluxtem(i,l,ctp) = tem * (gdlm-gcqlbl(i,l))
1173 qifluxtem(i,l,ctp) = tem * (gdim-gcqibl(i,l))
1174 do n = ntrq,ntr
1175 trfluxtem(i,l,n,ctp) = tem * (gdtrm(n)-gctrbl(i,l,n))
1176 enddo
1177! if(lprnt .and. i == ipr) write(0,*)' l=',l,' kbi=',kbi,' tem =', tem,' trfluxtem=',trfluxtem(l,ntr),&
1178! ' gdtrm=',gdtrm(ntr),' gctrbl=',gctrbl(i,l,ntr),' gq=',GDQ(I,l,ntr),GDQ(I,l1,ntr),' l1=',l1,' ctp=',ctp,&
1179! ' fsigma=',fsigma,' gcym=',gcym(i,l,ctp),' cbmfl=',cbmfl,' sigma=',sigma(i,k)
1180
1181! The following commented out by Moorthi on April 13, 2018 because tke below
1182! cloud base becomes too large otherwise when shoc is used
1183
1184!upstream - This better matches what the original CS tendencies do
1185! sfluxtem(l) = tem * (gdt(i,l)+gocp*(gdz(i,l)-gdzm(i,l))-gctbl(i,l))
1186! qvfluxtem(l) = tem * (gdq(i,l,1)-gcqbl(i,l))
1187! qlfluxtem(l) = tem * (gdq(i,l,3)-gcqlbl(i,l))
1188! qifluxtem(l) = tem * (gdq(i,l,2)-gcqibl(i,l))
1189! do n = ntrq,NTR
1190! trfluxtem(l,n) = tem * (gdq(i,l,n)-gctrbl(i,l,n))
1191! enddo
1192
1193 enddo
1194 else
1195! flux = mass flux * (updraft variable minus environment variable)
1196
1197! tem = - fsigma * mflx_e
1198 tem = - mflx_e
1199!centered
1200 sfluxtem(i,k,ctp) = tem * (gdtm(i,k)+gocp*gdzm(i,k)-gctm(i,k))
1201 qvfluxtem(i,k,ctp) = tem * (gdqm-gcqm(i,k,ctp))
1202 qlfluxtem(i,k,ctp) = tem * (gdlm-gclm(i,k))
1203 qifluxtem(i,k,ctp) = tem * (gdim-gcim(i,k))
1204 do n = ntrq,ntr
1205 trfluxtem(i,k,n,ctp) = tem * (gdtrm(n)-gctrm(i,k,n))
1206 enddo
1207
1208!upstream - This better matches what the original CS tendencies do
1209! if(k < kk) then
1210! sfluxtem(k) = tem * (gdt(i,k)+gocp*gdz(i,k)-gctm(i,k))
1211! qvfluxtem(k) = tem * (gdq(i,k,1)-gcqm(i,k))
1212! qlfluxtem(k) = tem * (gdq(i,k,3)-gclm(i,k))
1213! qifluxtem(k) = tem * (gdq(i,k,2)-gcim(i,k))
1214! do n = ntrq,NTR
1215! trfluxtem(k,n) = tem * (gdq(i,k,n)-gctrm(i,k,n))
1216! enddo
1217! if(lprnt .and. i == ipr) write(0,*)' k=',k,' kbi=',kbi,' tem =', tem,' kk=',kk,&
1218! ' gctrm=',gctrm(i,k,ntr),' gdq=',gdq(I,k,ntr),' gctrm=',gctrm(I,k,ntr),' ctp=',ctp,&
1219! ' fsigma=',fsigma,' mflx_e=',mflx_e,' trfluxtemk=',trfluxtem(k,ntr),' sigma=',sigma(i,k)
1220
1221! else
1222! centered at top of cloud
1223! sfluxtem(k) = tem * (gdtm(i,k)+gocp*gdzm(i,k)-gctm(i,k))
1224! qvfluxtem(k) = tem * (gdqm-gcqm(i,k))
1225! qlfluxtem(k) = tem * (gdlm-gclm(i,k))
1226! qifluxtem(k) = tem * (gdim-gcim(i,k))
1227! do n = ntrq,NTR
1228! trfluxtem(k,n) = tem * (gdtrm(n)-gctrm(i,k,n))
1229! enddo
1230! endif
1231
1232! if(lprnt .and. i == ipr) write(0,*)' k=',k,' kbi=',kbi,' tem =', tem,' kk=',kk,&
1233! ' gctrm=',gctrm(i,k,ntr),' gdtrm=',gdtrm(ntr),' gctrm=',gctrm(I,k,ntr),' ctp=',ctp,&
1234! ' fsigma=',fsigma,' mflx_e=',mflx_e,' trfluxtemk=',trfluxtem(k,ntr),' sigma=',sigma(i,k)
1235
1236
1237 endif ! if(k > kbi) then
1238 enddo ! end of k=kbi,kk loop
1239
1240 endif ! end of if(cbmfl > zero)
1241
1242
1243
1244 enddo ! end of i loop
1245 endif ! if (flx_form)
1246!
1247! we don't reduce these values in AW, just the tendencies due to fluxes
1248! do i=ists,iens
1249! if (cbmfx(i,ctp) > zero) then
1250! tem = one - sigma(i,kt(i,ctp))
1251! gcyt(i,ctp) = tem * gcyt(i,ctp)
1252! gtprt(i,ctp) = tem * gtprt(i,ctp)
1253! gclt(i,ctp) = tem * gclt(i,ctp)
1254! gcht(i,ctp) = tem * gcht(i,ctp)
1255! gcqt(i,ctp) = tem * gcqt(i,ctp)
1256! gcit(i,ctp) = tem * gcit(i,ctp)
1257! do n = ntrq,ntr
1258! gctrt(i,n,ctp) = tem * gctrt(i,n,ctp)
1259! enddo
1260! gcut(i,ctp) = tem * gcut(i,ctp)
1261! gcvt(i,ctp) = tem * gcvt(i,ctp)
1262! do k=1,kmax
1263! kk = kb(i)
1264! if (k < kk) then
1265! tem = one - sigma(i,kk)
1266! tem1 = tem
1267! else
1268! tem = one - sigma(i,k)
1269! tem1 = one - 0.5*(sigma(i,k)+sigma(i,k-1))
1270! endif
1271! gcym(i,k,ctp) = tem * gcym(i,k,ctp)
1272! gprciz(i,k) = tem1 * gprciz(i,k)
1273! gsnwiz(i,k) = tem1 * gsnwiz(i,k)
1274! gclz(i,k) = tem1 * gclz(i,k)
1275! gciz(i,k) = tem1 * gciz(i,k)
1276! enddo
1277! endif
1278! enddo
1279
1280!
1282 CALL cumflx(im , ijsdim, kmax , & !DD dimensions
1283 gmfx0 , gprci , gsnwi , cmdet, & ! output
1284 qliq , qice , gtprc0, & ! output
1285 cbmfx(:,ctp) , gcym(:,:,ctp), gprciz(:,:,ctp), gsnwiz(:,:,ctp) , & ! input
1286 gtprt(:,ctp) , gclz , gciz , gcyt(:,ctp),& ! input
1287 kb , kt(:,ctp) , ktmx(ctp) , & ! input
1288 ists , iens ) ! input
1289
1290 ENDDO ! end of cloud type ctp loop
1291
1293 do k=1,kmax
1294 do i=ists,iens
1295 gmflx(i,k) = gmfx0(i,k) ! contains net updraft mass flux for all clouds
1296 enddo
1297 enddo
1298 ktmxt = 3
1299 DO ctp=1,nctp
1300 IF (ktmx(ctp) > ktmxt) ktmxt = ktmx(ctp)
1301 ENDDO
1302
1303! DO K=1,KTMXT
1304! DO I=ISTS,IENS
1305! CUMCLW(I,K) = QLIQ(I,K) + QICE(I,K)
1306! IF (CUMCLW(I,K) > zero) THEN
1307! FLIQC(I,K) = QLIQ(I,K) / CUMCLW(I,K)
1308! ENDIF
1309! ENDDO
1310! ENDDO
1311!
1312! Cumulus Cloudiness
1313! CALL CUMCLD(IJSDIM, KMAX , & !DD dimensions
1314! CUMCLW, QLIQ , QICE , FLIQC , & ! modified
1315! CUMFRC, & ! output
1316! GMFLX , KTMXT , ISTS , IENS ) ! input
1317!
1318! - Call cumdet() to compute cloud detrainment heating.
1319 if (.not. flx_form) then
1320 CALL cumdet(im , ijsdim, kmax , ntr , ntrq , & !DD dimensions
1321 gtt , gtq , gtu , gtv , & ! modified
1322 gdh , gdq , gdu , gdv , & ! input
1323! GTT , GTQ , GTCFRC, GTU , GTV , & ! modified
1324! GDH , GDQ , GDCFRC, GDU , GDV , & ! input
1325 cbmfx , gcyt , delpinv , gcht , gcqt , & ! input
1326 gclt , gcit , gcut , gcvt , gdq(:,:,iti),& ! input
1327 gctrt , &
1328 kt , ists , iens, nctp ) ! input
1329 endif
1330
1331!for now area fraction of the downdraft is zero, it will be computed
1332! within cumdwn and applied there. So we will get the total sigma now before calling it,
1333! and apply to the diabatic terms from the updrafts.
1334
1335! if (do_aw.and.flx_form) then
1336 if (flx_form) then
1337 do k=1,kp1
1338 do i=ists,iens
1339 lamdamax = maxval(lamdai(i,k,:))
1340 do while (lamdamax > zero)
1341 loclamdamax = maxloc(lamdai(i,k,:),dim=1)
1342 lamdaprod(i,k) = lamdaprod(i,k) * (one+lamdai(i,k,loclamdamax))
1343 sigmai(i,k,loclamdamax) = lamdai(i,k,loclamdamax) / lamdaprod(i,k)
1344 sigma(i,k) = max(zero, min(one, sigma(i,k) + sigmai(i,k,loclamdamax)))
1345 vverti(i,k,loclamdamax) = sigmai(i,k,loclamdamax) * wcv(i,k,loclamdamax)
1346
1347 ! make this lamdai negative so it won't be counted again
1348 lamdai(i,k,loclamdamax) = -lamdai(i,k,loclamdamax)
1349 ! get new lamdamax
1350 lamdamax = maxval(lamdai(i,k,:))
1351 enddo
1352 ! restore original values of lamdai
1353 lamdai(i,k,:) = abs(lamdai(i,k,:))
1354! write(6,'(i2,14f7.4)') k,sigmai(i,k,:)
1355 enddo
1356 enddo
1357 endif
1358
1359! the condensation terms - these come from the MSE and condensed water budgets for
1360! an entraining updraft
1361 if(flx_form) then
1362 DO ctp=1,nctp ! loop over cloud types
1363 dtcondtem(:,:) = zero
1364 dqcondtem(:,:) = zero
1365 dqprectem(:,:) = zero
1366 dfrzprectem(:,:) = zero
1367 dtfrztem(:,:) = zero
1368 do i=ists,iens
1369 cbmfl = cbmfx(i,ctp)
1370 kk = kt(i,ctp) ! cloud top index
1371 if(cbmfl > zero) then ! this should avoid zero wcv in the denominator
1372 kbi = kb(i) ! cloud base index
1373 do k=kbi,kk ! loop from cloud base to cloud top
1374 km1 = k - 1
1375 rhs_h = zero
1376 rhs_q = zero
1377 if(k > kbi) then
1378! tem = cbmfl * (one - sigma(i,k))
1379 tem = cbmfl * (one - 0.5*(sigma(i,k)+sigma(i,km1)))
1380 tem1 = gcym(i,k,ctp) * (one - sigma(i,k))
1381 tem2 = gcym(i,km1,ctp) * (one - sigma(i,km1))
1382 rhs_h = cbmfl * (tem1*gchm(i,k,ctp) - (tem2*gchm(i,km1,ctp) &
1383 + gdh(i,km1)*(tem1-tem2)) )
1384 rhs_q = cbmfl * (tem1*(gcwm(i,k,ctp)-gcqm(i,k,ctp)) &
1385 - (tem2*(gcwm(i,km1,ctp)-gcqm(i,km1,ctp)) &
1386 + (gdw(i,km1)-gdq(i,km1,1))*(tem1-tem2)) )
1387!
1388 dqcondtem(i,km1) = -rhs_q ! condensation
1389 dqprectem(i,km1) = tem * (gprciz(i,k,ctp) + gsnwiz(i,k,ctp)) ! total precip production
1390 dfrzprectem(i,km1) = tem * gsnwiz(i,k,ctp) ! production of frozen precip
1391 dtfrztem(i,km1) = rhs_h*oneocp ! heating due to freezing
1392! total temperature tendency due to in cloud microphysics
1393 dtcondtem(i,km1) = - elocp * dqcondtem(i,km1) + dtfrztem(i,km1)
1394
1395 endif ! if(k > kbi) then
1396 enddo ! end of k=kbi,kk loop
1397
1398 endif ! end of if(cbmfl > zero)
1399
1400
1401! get tendencies by difference of fluxes, sum over cloud type
1402
1403 do k = 1,kk
1404! sum single cloud microphysical tendencies over all cloud types
1405 condtermt(i,k) = condtermt(i,k) + dtcondtem(i,k) * delpinv(i,k)
1406 condtermq(i,k) = condtermq(i,k) + dqcondtem(i,k) * delpinv(i,k)
1407 prectermq(i,k) = prectermq(i,k) + dqprectem(i,k) * delpinv(i,k)
1408 prectermfrz(i,k) = prectermfrz(i,k) + dfrzprectem(i,k) * delpinv(i,k)
1409 frzterm(i,k) = frzterm(i,k) + dtfrztem(i,k) * delpinv(i,k)
1410
1411! if (lprnt .and. i == ipr) write(0,*)' k=',k,' trfluxtem=',trfluxtem(k+1,ntr),trfluxtem(k,ntr),&
1412! ' ctp=',ctp,' trfluxterm=',trfluxterm(i,k,ntr)
1413 enddo
1414
1415 enddo ! end of i loop
1416 enddo ! end of nctp loop
1417 endif
1418!downdraft sigma and mass-flux tendency terms are now put into
1419! the nctp+1 slot of the cloud-type dimensiond variables
1420
1421 do k=1,kmax
1422 do i=ists,iens
1423 sigmad(i,k) = zero
1424 enddo
1425 enddo
1426
1429 CALL cumdwn(im, ijsdim, kmax, ntr, ntrq, nctp, & ! DD dimensions
1430 gtt , gtq , gtu , gtv , & ! modified
1431 gmflx , & ! modified updraft+downdraft flux
1432 gprcp , gsnwp , gtevp , gmdd , & ! output
1433 gprci , gsnwi , & ! input
1434 gdh , gdw , gdq , gdq(:,:,iti) , & ! input
1435 gdqs , gds , gdhs , gdt , & ! input
1436 gdu , gdv , gdz , & ! input
1437 gdzm , fdqs , delp , delpinv , & ! input
1438 sigmad, do_aw , do_awdd, flx_form, & ! DDsigma input
1439 dtmelt, dtevap, dtsubl, & ! DDsigma input
1440 dtdwn , dqvdwn, dqldwn, dqidwn, & ! DDsigma input
1441 dtrdwn, &
1442 kb , ktmxt , ists , iens ) ! input
1443
1444
1445! sigma = sigma + sigmad
1446
1448 if (.not. flx_form) then
1449! Cloud Subsidence Heating
1450! -----------------------=
1451 CALL cumsbh(im , ijsdim, kmax , ntr , ntrq , & !DD dimensions
1452 gtt , gtq , & ! modified
1453 gtu , gtv , & ! modified
1454 gdh , gdq , gdq(:,:,iti) , & ! input
1455 gdu , gdv , & ! input
1456 delpinv , gmflx , gmfx0 , & ! input
1457 ktmxt , cpres , kb, ists , iens ) ! input
1458 else
1459 CALL cumsbw(im , ijsdim, kmax , & !DD dimensions
1460 gtu , gtv , & ! modified
1461 gdu , gdv , & ! input
1462 delpinv , gmflx , gmfx0 , & ! input
1463 ktmxt , cpres , kb, ists , iens ) ! input
1464
1465 endif
1466!
1467! for now the following routines appear to be of no consequence - DD
1468!
1469 if (.not. flx_form) then
1470! Tracer Updraft properties
1471! -------------
1472 allocate (gprcc(ijsdim,ntr))
1473 do n=1,ntr
1474 do i=1,ijsdim
1475 gprcc(i,n) = zero
1476 enddo
1477 enddo
1478 CALL cumupr(im , ijsdim, kmax , ntr , & !DD dimensions
1479 gtq , gprcc , & ! modified
1480 gdq , cbmfx , & ! input
1481 gcym , gcyt , gcqt , gclt , gcit , & ! input
1482 gtprt , gtevp , gtprc0, & ! input
1483 kb , kbmx , kt , ktmx , ktmxt , & ! input
1484 delpinv , otspt1, ists , iens, & ! input
1485 fscav , fswtr, nctp)
1486!
1487! Tracer Change due to Downdraft
1488! ---------------
1489 CALL cumdnr(im ,ijsdim , kmax , ntr , & !DD dimensions
1490 gtq , & ! modified
1491 gdq , gmdd , delpinv , & ! input
1492 ktmxt , otspt1, ists , iens ) ! input
1493!!
1494!! Tracer change due to Subsidence
1495!! ---------------
1496!! This will be done by cumsbh, now DD 20170907
1497! CALL CUMSBR(im , IJSDIM, KMAX , NTR ,NCTP, & !DD dimensions
1498! GTQ , & ! modified
1499! GDQ , DELPI , & ! input
1500! GMFLX , KTMXT , OTSPT2, & ! input
1501! ISTS , IENS ) ! input
1502
1503 endif
1504
1505! if this tracer not advected zero it out
1506 DO n = ntrq,ntr
1507 if (.not. otspt2(n)) then
1508 DO k=1,kmax
1509 DO i=ists,iens
1510 gtq(i,k,n) = 0.0
1511 ENDDO
1512 ENDDO
1513 endif
1514 ENDDO
1515
1516! if(do_aw .and. flx_form) then ! compute AW tendencies
1518 if(flx_form) then ! compute AW tendencies
1519 ! AW lump all heating together, compute qv term
1520
1521! sigma interpolated to the layer for condensation, etc. terms, precipitation
1522 if(do_aw) then
1523 do k=1,kmax
1524 kp1 = k+1
1525 do i=1,ijsdim
1526 fsigma(i,k) = one - half*(sigma(i,k)+sigma(i,kp1))
1527 enddo
1528 enddo
1529 else
1530 do k=1,kmax+1
1531 do i=1,ijsdim
1532 fsigma(i,k) = one
1533 enddo
1534 enddo
1535 endif
1536
1537! AW adjustment of precip fluxes from downdraft model
1538 if(do_aw) then
1539 kp1 = kmax+1
1540 DO i=ists,iens
1541 gsnwp( i,kp1 ) = zero
1542 gprcp( i,kp1 ) = zero
1543 ENDDO
1544 tem1 = cpoemelt/grav
1545 tem2 = cpoel/grav
1546 tem3 = cpoesub/grav
1547 DO k=kmax,1,-1
1548 kp1 = k+1
1549 DO i=ists,iens
1550 tem = -dtmelt(i,k) * delp(i,k) * tem1
1551 teme = -dtevap(i,k) * delp(i,k) * tem2
1552 tems = -dtsubl(i,k) * delp(i,k) * tem3
1553 gsnwp(i,k) = gsnwp(i,kp1) + fsigma(i,k) * (gsnwi(i,k) - tem - tems)
1554 gprcp(i,k) = gprcp(i,kp1) + fsigma(i,k) * (gprci(i,k) + tem - teme)
1555 ENDDO
1556 ENDDO
1557 endif
1558
1559
1560! some of the above routines have set the tendencies and they need to be
1561! reinitialized, gtt not needed, but gtq needed Anning 5/25/2020
1562 do n=1,ntr
1563 do k=1,kmax
1564 do i=1,ijsdim
1565 gtq(i,k,n) = zero
1566 enddo
1567 enddo
1568 enddo
1569! do k=1,kmax
1570! do i=1,ijsdim
1571! gtt(i,k) = zero
1572! enddo
1573! enddo
1574 do k=1,kmax
1575 do i=ists,iens
1576 dqevap(i,k) = - dtevap(i,k)*cpoel - dtsubl(i,k)*cpoesub
1577 dtevap(i,k) = dtevap(i,k) + dtsubl(i,k)
1578 dtsubl(i,k) = zero
1579 enddo
1580 enddo
1581
1582
1583! diabatic terms from updraft and downdraft models
1584 DO k=1,kmax
1585 DO i=ists,iens
1586 tem = frzterm(i,k)*cpoemelt - prectermfrz(i,k)
1587! gtt(i,k) = gtt(i,k) + fsigma(i,k)*(dtmelt(i,k) + dtevap(i,k)) + condtermt(i,k)
1588! gtq(i,k,1) = gtq(i,k,1) + fsigma(i,k)*dqevap(i,k) + condtermq(i,k)
1589! gtq(i,k,itl) = gtq(i,k,itl) - (condtermq(i,k) + prectermq(i,k) + tem)
1590! gtq(i,k,iti) = gtq(i,k,iti) + tem
1591 gtt(i,k) = dtdwn(i,k) + condtermt(i,k) &
1592 + fsigma(i,k)*(dtmelt(i,k) + dtevap(i,k))
1593 gtq(i,k,1) = dqvdwn(i,k) + condtermq(i,k) &
1594 + fsigma(i,k) * dqevap(i,k)
1595 gtq(i,k,itl) = dqldwn(i,k) - condtermq(i,k) &
1596 - prectermq(i,k) - tem
1597 gtq(i,k,iti) = dqidwn(i,k) + tem
1598
1599
1600! detrainment terms get zeroed
1601! gtldet(i,k) = zero
1602! gtidet(i,k) = zero
1603 ENDDO
1604 ENDDO
1605!! flux tendencies - compute the vertical flux divergence
1606 DO ctp =1,nctp
1607 DO i=ists,iens
1608 cbmfl = cbmfx(i,ctp)
1609 kk = kt(i,ctp) ! cloud top index
1610 if(cbmfl > zero) then ! this should avoid zero wcv in the denominator
1611 DO k=1,kk
1612 kp1 = k+1
1613 gtt(i,k) = gtt(i,k) - (fsigma(i,kp1)*sfluxtem(i,kp1,ctp) &
1614 - fsigma(i,k)*sfluxtem(i,k,ctp)) * delpinv(i,k)
1615 gtq(i,k,1) = gtq(i,k,1) - (fsigma(i,kp1)*qvfluxtem(i,kp1,ctp) &
1616 - fsigma(i,k)*qvfluxtem(i,k,ctp)) * delpinv(i,k)
1617 gtq(i,k,itl) = gtq(i,k,itl) - (fsigma(i,kp1)*qlfluxtem(i,kp1,ctp) &
1618 - fsigma(i,k)*qlfluxtem(i,k,ctp)) * delpinv(i,k)
1619 gtq(i,k,iti) = gtq(i,k,iti) - (fsigma(i,kp1)*qifluxtem(i,kp1,ctp) &
1620 - fsigma(i,k)*qifluxtem(i,k,ctp)) * delpinv(i,k)
1621 ENDDO
1622! replace tracer tendency only if to be advected.
1623 DO n = ntrq,ntr
1624 if (otspt2(n)) then
1625 DO k=1,kk
1626 kp1 = k+1
1627 gtq(i,k,n) = - (fsigma(i,kp1)*trfluxtem(i,kp1,n,ctp) &
1628 - fsigma(i,k)*trfluxtem(i,k,n,ctp)) * delpinv(i,k)
1629 ENDDO
1630 endif
1631 ENDDO
1632 end if
1633 ENDDO
1634 ENDDO
1635
1636! if(kdt>4) stop 1000
1637 DO i=ists,iens
1638 moistening_aw(i) = zero
1639 enddo
1640
1641! adjust tendencies that will lead to negative water mixing ratios
1642 tem2 = one / delta
1643 DO k=1,kmax
1644 DO i=ists,iens
1645 tem1 = - gdq(i,k,itl)*tem2
1646 if (gtq(i,k,itl) < tem1) then
1647 tem3 = gtq(i,k,itl) - tem1
1648 gtq(i,k,1) = gtq(i,k,1) + tem3
1649 gtq(i,k,itl) = tem1
1650 gtt(i,k) = gtt(i,k) - elocp*tem3
1651 endif
1652 tem1 = - gdq(i,k,iti)*tem2
1653 if (gtq(i,k,iti) < tem1) then
1654 tem3 = gtq(i,k,iti) - tem1
1655 gtq(i,k,1) = gtq(i,k,1) + tem3
1656 gtq(i,k,iti) = tem1
1657 gtt(i,k) = gtt(i,k) - esubocp*tem3
1658 endif
1659 tem1 = - gdq(i,k,1)*tem2
1660 if (gtq(i,k,1) < tem1) then
1661 gtt(i,k) = gtt(i,k) + elocp*(gtq(i,k,1)-tem1)
1662 gtq(i,k,1) = tem1
1663 endif
1664
1665! column-integrated total water tendency - to be used to impose water conservation
1666 moistening_aw(i) = moistening_aw(i) &
1667 + (gtq(i,k,1)+gtq(i,k,itl)+gtq(i,k,iti)) * delp(i,k) * gravi
1668 ENDDO
1669 ENDDO
1670
1671! replace tracer tendency only if to be advected.
1672 DO n = ntrq,ntr
1673 if (otspt2(n)) then
1674 DO k=1,kmax
1675 DO i=ists,iens
1676 gtq(i,k,n) = gtq(i,k,n) + dtrdwn(i,k,n)
1677 ENDDO
1678 ENDDO
1679 endif
1680 ENDDO
1681! if (lprnt) write(0,*)' endcs_cum gtq=',gtq(ipr,1:25,ntr)
1682! if (lprnt) write(0,*)' endcs_cum trfluxterm=',trfluxterm(ipr,1:25,ntr)
1683
1684 endif ! if (flx_form)
1685
1686!!!! this section may need adjustment for cloud ice and water with flux_form
1687!
1688! do k=1,kmax
1689! do i=ISTS,IENS
1690! GTQ(I,k,ITI) = GTQI(I,k)
1691! enddo
1692! enddo
1693!
1695 CALL cumfxr(im , ijsdim, kmax , ntr , & !DD dimensions
1696 gtq , & ! modified
1697 gdq , delp , delta , ktmxt , imfxr, & ! input
1698 ists , iens ) ! input
1699
1700!
1701! do k=1,kmax
1702! do i=ISTS,IENS
1703! GTQL(I,k) = GTQ(I,k,ITL) + GTLDET(I,k) + GTIDET(I,k)
1704! enddo
1705! enddo
1706!
1707! Tracer mass fixer with detrainment
1708! CALL CUMFXR1(IM , IJSDIM, KMAX , & !DD dimensions
1709! GTQL , & ! modified
1710! GDQ(1,1,ITL), DELP, DELTA, KTMXT, IMFXR(ITL), & ! input
1711! ISTS , IENS ) ! input
1712!
1713!!!!! end fixer section
1714
1715! DO K=1,KMAX
1716! DO I=ISTS, IENS
1717! GTLDET(I,k) = GTQL(I,k) - GTQ(I,k,ITL) - GTIDET(I,k)
1718
1719! tendencies of subgrid PDF (turned off)
1720! GDCLDX = GDCFRC( I,K ) + GTCFRC( I,K )*DELTA
1721! GDCLDX = MIN( MAX( GDCLDX, 0.D0 ), one )
1722! GTCFRC( I,K ) = ( GDCLDX - GDCFRC( I,K ) )/DELTA
1723!
1724! GDMU2X = GDQ( I,K,IMU2 ) + GTQ( I,K,IMU2 )*DELTA
1725! GDMU2X = MIN( MAX( GDMU2X,VARMIN ),VARMAX )
1726! GDMU3X = GDQ( I,K,IMU3 ) + GTQ( I,K,IMU3 )*DELTA
1727! GDMU3X = MIN( MAX( GDMU3X,-SKWMAX ),SKWMAX )
1728! GTQ( I,K,IMU2 ) = ( GDMU2X - GDQ( I,K,IMU2 ))/DELTA
1729! GTQ( I,K,IMU3 ) = ( GDMU3X - GDQ( I,K,IMU3 ))/DELTA
1730!
1731! tem = DELP(I,K)*GRAVI
1732! HBGT(I) = HBGT(I) + (CP*GTT(I,K) + EL*GTQ(I,K,1) &
1733! - EMELT*GTQ(I,K,ITI)) * tem
1734! - EMELT*(GTQ(I,K,ITI)+GTIDET(I,K))) * tem
1735! WBGT(I) = WBGT(I) + (GTQ(I,K,1) + GTQ(I,K,ITL) + GTQ(I,K,ITI)) * tem
1736! + GTLDET(I,K) + GTIDET(I,K)) * tem
1737! ENDDO
1738! ENDDO
1739
1740
1741!
1742! DO I=ISTS,IENS
1743! HBGT(I) = HBGT(I) - EMELT*GSNWC(I)
1744! WBGT(I) = WBGT(I) + GPRCC(I,1) + GSNWC(I)
1745! CTOPP(I) = 1.D6
1746! ENDDO
1747!
1748! The following commented out because they are unused
1749! DO CTP=1,NCTP
1750! DO I=ISTS, IENS
1751! kk = kt(i,ctp)
1752! IF (KK > KB(I) ) THEN
1753! CUMHGT(I,CTP) = GDZ(I,KK)
1754! CTOPP(I) = MIN(CTOPP(I), GDP(I,KK))
1755! ELSE
1756! CUMHGT (I,CTP) = -999.D0
1757! ENDIF
1758! ENDDO
1759! ENDDO
1760! DO I=ISTS,IENS
1761! IF(CTOPP(I) >= 1.D6) THEN
1762! CTOPP(I) = -999.D0
1763! ENDIF
1764! ENDDO
1765!
1767!In fact, no adjustment of the precip
1768! is occuring now which is a good sign! DD
1769 if(flx_form) then
1770 DO i = ists, iens
1771 if(gprcp(i,1)+gsnwp(i,1) > 1.e-12_kind_phys) then
1772 moistening_aw(i) = -moistening_aw(i) / (gprcp(i,1)+gsnwp(i,1))
1773! print*,'moistening_aw',moistening_aw(i)
1774 gprcp(i,:) = gprcp(i,:) * moistening_aw(i)
1775 gsnwp(i,:) = gsnwp(i,:) * moistening_aw(i)
1776 endif
1777 END DO
1778 endif
1779
1780! second method of determining sfc precip only
1781! if(flx_form) then
1782! DO I = ISTS, IENS
1783! pr_tot = zero
1784! pr_liq = zero
1785! pr_ice = zero
1786! do k = 1,kmax
1787! pr_tot = pr_tot - (gtq(i,k,1)+gtq(i,k,itl)+gtq(i,k,iti)) * delp(i,k) * gravi
1788! pr_ice = pr_ice + ( cp*gtt(i,k) + el*gtq(i,k,1) - emelt*gtq(i,k,iti) ) &
1789! * delp(i,k)*gravi
1790! enddo
1791 !pr_ice = max( min(pr_tot, pr_ice / (emelt)),zero)
1792! pr_ice = pr_ice / emelt
1793! pr_liq = pr_tot - pr_ice
1794! END DO
1795! print *,'precip1',pr_tot*86400.,gprcp(1,1)*86400.,gsnwp(1,1)*86400.
1796! print *,'precip2',pr_tot*86400.,pr_liq*86400.,pr_ice*86400.
1797! endif
1798
1799 DO k = 1, kmax
1800 DO i = ists, iens
1801 gprcpf( i,k ) = 0.5*( gprcp( i,k )+gprcp( i,k+1 ) )
1802 gsnwpf( i,k ) = 0.5*( gsnwp( i,k )+gsnwp( i,k+1 ) )
1803 END DO
1804 END DO
1805
1806!
1807! do i=ISTS,IENS
1808! GPRCC( I,1 ) = GPRCP( I,1 )
1809! GSNWC( I ) = GSNWP( I,1 )
1810! enddo
1811
1812! adjust sfc precip consistently with vertically integrated
1813! temperature and moisture tendencies
1814
1815 do k=1,kmax+1
1816 do i=ists,iens
1817 gtprp(i,k) = gprcp(i,k) + gsnwp(i,k)
1818 enddo
1819 enddo
1820!
1821!DD provide GFS with a separate downdraft mass flux
1822 if(do_aw) then
1823 DO k = 1, kmax+1
1824 DO i = ists, iens
1825 fsigma(i,k) = one - sigma(i,k)
1826 gmfx0( i,k ) = gmfx0( i,k ) * fsigma(i,k)
1827 gmflx( i,k ) = gmflx( i,k ) * fsigma(i,k)
1828 END DO
1829 END DO
1830 endif
1831 DO k = 1, kmax+1
1832 DO i = ists, iens
1833 gmfx1( i,k ) = gmfx0( i,k ) - gmflx( i,k )
1834 END DO
1835 END DO
1836
1837 if (allocated(gprcc)) deallocate(gprcc)
1838
1839!
1840 END SUBROUTINE cs_cumlus
1842!***********************************************************************
1845 SUBROUTINE cumbas & ! cloud base
1846 ( ijsdim, kmax , & !DD dimensions
1847 kb , gcym , kbmx , & ! output
1848 ntr , ntrq , &
1849 gchb , gcwb , gcub , gcvb , & ! output
1850 gcib , gctrb, & ! output
1851 gdh , gdw , gdhs , gdqs , & ! input
1852 gdqi , gdu , gdv , gdzm , & ! input
1853 gdpm , fdqs , gam , & ! input
1854 lprnt, ipr, &
1855 ists , iens , gctbl, gcqbl ,gdq, &
1856 gcwbl, gcqlbl, gcqibl, gctrbl ) ! input !DDsigmadiag add updraft profiles below cloud base
1857!
1858!
1859 IMPLICIT NONE
1860! integer, parameter :: crtrh=0.80
1861 integer, parameter :: crtrh=0.70
1862 INTEGER, INTENT(IN) :: IJSDIM, KMAX , ntr, ntrq ! DD, for GFS, pass in
1863 integer ipr
1864 logical lprnt
1865!
1866! [OUTPUT]
1867 INTEGER KB (IJSDIM) ! cloud base
1868 REAL(kind_phys) GCYM (IJSDIM, KMAX) ! norm. mass flux (half lev)
1869 INTEGER KBMX
1870 REAL(kind_phys) GCHB (IJSDIM) ! cloud base MSE
1871 REAL(kind_phys) GCWB (IJSDIM) ! cloud base total water
1872 REAL(kind_phys) GCUB (IJSDIM) ! cloud base U
1873 REAL(kind_phys) GCVB (IJSDIM) ! cloud base V
1874 REAL(kind_phys) GCIB (IJSDIM) ! cloud base ice
1875 REAL(kind_phys) GCtrB (IJSDIM,ntrq:ntr) ! cloud base tracer
1876
1877!DDsigma added to arglist for AW, subcloud updraft profiles: temperature, water vapor
1878! total water, cloud water, and cloud ice respectively
1879 REAL(kind_phys), dimension(ijsdim,kmax) :: gctbl, gcqbl, gcwbl, gcqlbl, gcqibl
1880 REAL(kind_phys), dimension(ijsdim,kmax,ntrq:ntr) :: gctrbl !DDsigmadiag
1881!
1882! [INPUT]
1883 REAL(kind_phys) GDH (IJSDIM, KMAX) ! moist static energy
1884 REAL(kind_phys) GDW (IJSDIM, KMAX) ! total water
1885 REAL(kind_phys) GDq (IJSDIM, KMAX, ntr) ! water vapor and tracer
1886 REAL(kind_phys) GDHS (IJSDIM, KMAX) ! saturate MSE
1887 REAL(kind_phys) GDQS (IJSDIM, KMAX) ! saturate humidity
1888 REAL(kind_phys) GDQI (IJSDIM, KMAX) ! cloud ice
1889 REAL(kind_phys) GDU (IJSDIM, KMAX) ! u-velocity
1890 REAL(kind_phys) GDV (IJSDIM, KMAX) ! v-velocity
1891 REAL(kind_phys) GDZM (IJSDIM, KMAX+1) ! Altitude (half lev)
1892 REAL(kind_phys) GDPM (IJSDIM, KMAX+1) ! pressure (half lev)
1893 REAL(kind_phys) FDQS (IJSDIM, KMAX)
1894 REAL(kind_phys) GAM (IJSDIM, KMAX)
1895 INTEGER ISTS, IENS
1896!
1897! [INTERNAL WORK]
1898 REAL(kind_phys) CBASE (IJSDIM) ! one over cloud base height
1899! REAL(kind_phys) CBASEP(IJSDIM) ! cloud base pressure
1900 REAL(kind_phys) DELZ, GAMX, wrk
1901! REAL(kind_phys) DELZ, QSL, GAMX, wrk
1902! REAL(kind_phys), dimension(ijsdim,kmax) :: gchbl !DDsigmadiag
1903 real(kind_phys), dimension(ijsdim) :: gcqb, tx1, spbl, qsl
1904 INTEGER I, K, kp1, n
1905!
1906! [INTERNAL PARM]
1907 INTEGER :: KMAXM1
1908 INTEGER :: KLCLB !! LCL base level
1909 INTEGER :: KCB !! fix cloud bottom
1910 INTEGER :: KBMAX !! cloud base max
1911 INTEGER :: KBOFS !! cloud base offset
1912
1913 kmaxm1 = kmax-1
1914 klclb = 1 ! LCL base level
1915 kcb = 0 ! fix cloud bottom
1916! KBMAX = KMAXM1 ! cloud base max
1917 kbmax = kmax/2 ! cloud base max
1918 kbofs = 0 ! cloud base offset
1919!
1920 do k=1,kmax
1921 do i=ists,iens
1922 gcym(i,k) = zero
1923 enddo
1924 enddo
1925!
1926 IF (kcb > 0) THEN
1927 DO i=ists,iens
1928 kb(i) = kcb
1929 ENDDO
1930 ELSE
1931 DO i=ists,iens
1932! KB(I) = KBMAX
1933 kb(i) = -1
1934 tx1(i) = one / gdpm(i,1)
1935 ENDDO
1936 DO k=klclb+1,kbmax-1
1937 DO i=ists,iens
1938 gamx = fdqs(i,k) / (one+gam(i,k)) * oneocp
1939 qsl(i) = gdqs(i,k) + gamx * (gdh(i,klclb)-gdhs(i,k))
1940 spbl(i) = one - gdpm(i,k) * tx1(i)
1941 IF (gdw(i,klclb) >= qsl(i) .and. kb(i) < 0 &
1942 .and. spbl(i) >= spblmin) THEN
1943! .and. spbl(i) >= spblcrit .and. spbl(i) < spblcrit*10.0) THEN
1944 kb(i) = k + kbofs
1945 ENDIF
1946 ENDDO
1947 ENDDO
1948 ENDIF
1949!
1950 do i=ists,iens
1951 tx1(i) = zero
1952 qsl(i) = zero
1953 enddo
1954 do k=1,kbmax
1955 do i=ists,iens
1956 if (k < kb(i)) then
1957 tx1(i) = tx1(i) + gdw(i,k) * (gdzm(i,k+1)-gdzm(i,k))
1958 qsl(i) = qsl(i) + gdqs(i,k) * (gdzm(i,k+1)-gdzm(i,k))
1959 endif
1960 enddo
1961 enddo
1962 do i=ists,iens
1963 if (qsl(i) > zero) tx1(i) = tx1(i) / qsl(i)
1964 if (tx1(i) < crtrh) kb(i) = -1
1965 enddo
1966
1967!
1968 kbmx = 1
1969 DO i=ists,iens
1970 kbmx = max(kbmx, kb(i))
1971 if (kb(i) > 0) then
1972 cbase(i) = one / (gdzm(i,kb(i)) - gdzm(i,1))
1973! CBASEP(I) = GDPM(I,KB(I))
1974 endif
1975 ENDDO
1976!
1977 DO k=2,kbmx
1978 DO i=ists,iens
1979 IF (k <= kb(i)) THEN
1980! GCYM(I,K) = sqrt((GDZM(I,K)-GDZM(I,1))*CBASE(i))
1981 gcym(i,k) = (gdzm(i,k)-gdzm(i,1))*cbase(i)
1982 ENDIF
1983 ENDDO
1984 ENDDO
1985!
1986 DO i=ists,iens
1987 gchb(i) = zero
1988 gcwb(i) = zero
1989 gcub(i) = zero
1990 gcvb(i) = zero
1991 gcib(i) = zero
1992 gcqb(i) = zero
1993 ENDDO
1994 do n = ntrq,ntr
1995 DO i=ists,iens
1996 gctrb(i,n) = zero
1997 enddo
1998 enddo
1999 do k=1,kmax
2000 do i=ists,iens
2001! GChbl(i,k) = zero
2002 gcqbl(i,k) = zero
2003 gcqlbl(i,k) = zero
2004 gcqibl(i,k) = zero
2005 gctbl(i,k) = zero
2006 gcwbl(i,k) = zero
2007 enddo
2008 enddo
2009! do n=ntrq,ntr
2010! do k=1,kmax
2011! do i=ists,iens
2012! gtrqbl(i,k,n) = zero
2013! enddo
2014! enddo
2015! enddo
2016!
2017 DO k=1,kbmx
2018 kp1 = min(k+1, kmax)
2019 DO i=ists,iens
2020 IF (k < kb(i)) THEN
2021 delz = gcym(i,kp1) - gcym(i,k)
2022 gchb(i) = gchb(i) + delz * gdh(i,k)
2023 gcwb(i) = gcwb(i) + delz * gdw(i,k)
2024 gcub(i) = gcub(i) + delz * gdu(i,k)
2025 gcvb(i) = gcvb(i) + delz * gdv(i,k)
2026 gcib(i) = gcib(i) + delz * gdqi(i,k)
2027 gcqb(i) = gcqb(i) + delz * gdq(i,k,1)
2028! do n = ntrq,ntr
2029! GCtrB(I,n) = GCtrB(I,n) + DELZ * GDQ (I,K,n)
2030! enddo
2031! get subcloud profiles to pass out and do AW eddy flux tendencies
2032! removing the normalized mass flux weighting
2033 wrk = one / gcym(i,kp1)
2034! gchbl(i,kp1) = gchb(i) * wrk
2035 gcqbl(i,kp1) = gcqb(i) * wrk
2036 gcqibl(i,kp1) = gcib(i) * wrk
2037 gcwbl(i,kp1) = gcwb(i) * wrk
2038 gcqlbl(i,kp1) = gcwbl(i,kp1) - (gcqibl(i,kp1)+gcqbl(i,kp1))
2039! gctbl(i,kp1) = (gchbl(i,kp1) - grav*gdzm(i,kp1) - el*gcqbl(i,kp1)) * oneocp
2040 gctbl(i,kp1) = (gchb(i)*wrk - grav*gdzm(i,kp1) - el*gcqbl(i,kp1)) * oneocp
2041! tracers
2042 do n=ntrq,ntr
2043 gctrb(i,n) = gctrb(i,n) + delz * gdq(i,k,n)
2044 gctrbl(i,kp1,n) = gctrb(i,n) * wrk
2045 enddo
2046 ENDIF
2047 ENDDO
2048 ENDDO
2049!
2050 END SUBROUTINE cumbas
2051!***********************************************************************
2054 SUBROUTINE cumup & !! in-cloud properties
2055 ( ijsdim, kmax , ntr , ntrq , & !DD dimensions
2056 acwf , & ! output
2057 gclz , gciz , gprciz, gsnwiz, & ! output
2058 gcyt , gcht , gcqt , & ! output
2059 gclt , gcit , gtprt , & ! output
2060 gcut , gcvt , gctrt , & ! output
2061 kt , ktmx , & ! output
2062 gcym , & ! modified
2063 wcv , & ! !DDsigma new output
2064 gchb , gcwb , gcub , gcvb , & ! input !DDsigmadiag
2065 gcib , gctrb , & ! input
2066 gdu , gdv , gdh , gdw , & ! input
2067 gdhs , gdqs , gdt , gdtm , & ! input
2068 gdq , gdqi , gdz , gdzm , & ! input
2069 gdpm , fdqs , gam , gdztr , & ! input
2070 cpres , wcb , & ! input
2071! CPRES , WCB , ERMR , & ! input
2072 kb , ctp , ists , iens, & ! input
2073 gctm , gcqm , gcwm , gchm, gcwt,&
2074 gclm, gcim , gctrm , lprnt, ipr )
2075!
2076!DD AW the above line of arguments were previously local, and often scalars.
2077! Dimensions were added to them to save profiles for each grid point.
2078!
2079 IMPLICIT NONE
2080
2081 INTEGER, INTENT(IN) :: IJSDIM, KMAX, NTR, ipr , ntrq ! DD, for GFS, pass in
2082 logical :: lprnt
2083!
2084! [OUTPUT]
2085 REAL(kind_phys) ACWF (IJSDIM)
2086 REAL(kind_phys) GCLZ (IJSDIM, KMAX)
2087 REAL(kind_phys) GCIZ (IJSDIM, KMAX)
2088 REAL(kind_phys) GPRCIZ(IJSDIM, KMAX+1)
2089 REAL(kind_phys) GSNWIZ(IJSDIM, KMAX+1)
2090 REAL(kind_phys) GCYT (IJSDIM)
2091 REAL(kind_phys) GCHT (IJSDIM)
2092 REAL(kind_phys) GCQT (IJSDIM)
2093 REAL(kind_phys) GCLT (IJSDIM)
2094 REAL(kind_phys) GCIT (IJSDIM)
2095 REAL(kind_phys) GCtrT (IJSDIM, ntrq:ntr)
2096 REAL(kind_phys) GTPRT (IJSDIM)
2097 REAL(kind_phys) GCUT (IJSDIM)
2098 REAL(kind_phys) GCVT (IJSDIM)
2099 REAL(kind_phys) GCwT (IJSDIM)
2100 INTEGER KT (IJSDIM)
2101 INTEGER KTMX
2102 REAL(kind_phys) WCV (IJSDIM, KMAX+1)
2103!
2104! [MODIFIED]
2105 REAL(kind_phys) GCYM (IJSDIM, KMAX)
2106!
2107! [INPUT]
2108 REAL(kind_phys) GCHB (IJSDIM)
2109 REAL(kind_phys) GCWB (IJSDIM)
2110 REAL(kind_phys) GCUB (IJSDIM)
2111 REAL(kind_phys) GCVB (IJSDIM)
2112 REAL(kind_phys) GCIB (IJSDIM)
2113 REAL(kind_phys) GCtrB (IJSDIM,ntrq:ntr)
2114 REAL(kind_phys) GDU (IJSDIM, KMAX)
2115 REAL(kind_phys) GDV (IJSDIM, KMAX)
2116 REAL(kind_phys) GDH (IJSDIM, KMAX)
2117 REAL(kind_phys) GDW (IJSDIM, KMAX)
2118 REAL(kind_phys) GDHS (IJSDIM, KMAX)
2119 REAL(kind_phys) GDQS (IJSDIM, KMAX)
2120 REAL(kind_phys) GDT (IJSDIM, KMAX)
2121 REAL(kind_phys) GDTM (IJSDIM, KMAX+1)
2122 REAL(kind_phys) GDQ (IJSDIM, KMAX, NTR)
2123 REAL(kind_phys) GDQI (IJSDIM, KMAX)
2124 REAL(kind_phys) GDZ (IJSDIM, KMAX)
2125 REAL(kind_phys) GDZM (IJSDIM, KMAX+1)
2126 REAL(kind_phys) GDPM (IJSDIM, KMAX+1)
2127 REAL(kind_phys) FDQS (IJSDIM, KMAX)
2128 REAL(kind_phys) GAM (IJSDIM, KMAX)
2129 REAL(kind_phys) GDZTR (IJSDIM)
2130 REAL(kind_phys) CPRES
2131 REAL(kind_phys) WCB(ijsdim)
2132! REAL(kind_phys) ERMR !< entrainment rate (ASMODE)
2133 INTEGER KB (IJSDIM)
2134 INTEGER CTP, ISTS, IENS
2135!
2136! [INTERNAL WORK]
2137 REAL(kind_phys) ACWFK (IJSDIM,KMAX)
2138 REAL(kind_phys) ACWFN (IJSDIM,KMAX)
2139 REAL(kind_phys) myGCHt
2140 REAL(kind_phys) GCHMZ (IJSDIM, KMAX)
2141 REAL(kind_phys) GCWMZ (IJSDIM, KMAX)
2142 REAL(kind_phys) GCUMZ (IJSDIM, KMAX)
2143 REAL(kind_phys) GCVMZ (IJSDIM, KMAX)
2144 REAL(kind_phys) GCqMZ (IJSDIM )
2145 REAL(kind_phys) GCIMZ (IJSDIM, KMAX)
2146 REAL(kind_phys) GCtrMZ(IJSDIM, KMAX,ntrq:ntr)
2147 REAL(kind_phys) GTPRMZ(IJSDIM, KMAX)
2148!
2149 REAL(kind_phys) BUOY (IJSDIM, KMAX)
2150 REAL(kind_phys) BUOYM (IJSDIM, KMAX)
2151 REAL(kind_phys) WCM (IJSDIM )
2152! REAL(kind_phys) WCM (IJSDIM, KMAX) !< updraft velocity**2 (half lev)
2153!DD sigma make output REAL(kind_phys) WCV ( IJSDIM, KMAX+1 ) !! updraft velocity (half lev)
2154 REAL(kind_phys) GCY (IJSDIM, KMAX)
2155! REAL(kind_phys) ELAR (IJSDIM, KMAX) !< entrainment rate
2156 REAL(kind_phys) ELAR
2157!
2158 REAL(kind_phys) GCHM (IJSDIM, KMAX+1)
2159 REAL(kind_phys) GCWM (IJSDIM, KMAX+1)
2160 REAL(kind_phys) GCTM (IJSDIM, KMAX+1)
2161 REAL(kind_phys) GCQM (IJSDIM, KMAX+1)
2162 REAL(kind_phys) GCLM (IJSDIM, KMAX+1)
2163 REAL(kind_phys) GCIM (IJSDIM, KMAX+1)
2164 REAL(kind_phys) GCUM (IJSDIM, KMAX)
2165 REAL(kind_phys) GCVM (IJSDIM, KMAX)
2166 REAL(kind_phys) GCtrM (IJSDIM, KMAX,ntrq:ntr)
2167!
2168 REAL(kind_phys), dimension(IJSDIM) :: WCM_, ELARM1, GDZMKB
2169 REAL(kind_phys) GDQSM, GDHSM, GDQM, GDSM, GDCM, FDQSM, GCCM, &
2170 DELZ, ELADZ, DCTM , CPGMI, DELC, FICE, ELARM2,GCCMZ, &
2171 PRECR, GTPRIZ, DELZL, GCCT, DCT, WCVX, PRCZH, wrk
2172 INTEGER K, I, kk, km1, kp1, n
2173
2174! CHARACTER CTNUM*2
2175!
2176!DD#ifdef OPT_CUMBGT
2177!DD REAL(kind_phys) HBGT (IJSDIM) ! heat budget
2178!DD REAL(kind_phys) WBGT (IJSDIM) ! water budget
2179!DD REAL(kind_phys) PBGT (IJSDIM) ! precipitation budget
2180!DD REAL(kind_phys) MBGT (IJSDIM) ! mass budget
2181!DD REAL(kind_phys) GTPRX (IJSDIM) ! (rain+snow)*eta at top
2182!DD REAL(kind_phys) GSNWT (IJSDIM) ! cloud top snow*eta
2183!DD REAL(kind_phys) HBMX, WBMX, PBMX, MBMX
2184!DD SAVE HBMX, WBMX, PBMX, MBMX
2185!DD#endif
2186!
2187! [INTERNAL PARAM]
2188
2189 REAL(kind_phys), parameter :: ZTREF = 1._kind_phys, ztrefi = one/ztref, &
2190 elamin = zero, elamax = 4.e-3 ! min and max entrainment rate
2191 REAL(kind_phys) :: PB = 1.0_kind_phys
2192!m REAL(kind_phys) :: TAUZ = 5.0e3_kind_phys
2193 REAL(kind_phys) :: TAUZ = 1.0e4_kind_phys
2194!m REAL(kind_phys) :: ELMD = 2.4e-3 ! for Neggers and Siebesma (2002)
2195!m REAL(kind_phys) :: ELAMAX = 5.e-3 ! max. of entrainment rate
2196! REAL(kind_phys) :: WCCRT = zero
2197!m REAL(kind_phys) :: WCCRT = 0.01
2198 REAL(kind_phys) :: WCCRT = 1.0e-6_kind_phys, wvcrt=1.0e-3_kind_phys
2199 REAL(kind_phys) :: TSICE = 273.15_kind_phys ! compatible with macrop_driver
2200 REAL(kind_phys) :: TWICE = 233.15_kind_phys ! compatible with macrop_driver
2201 REAL(kind_phys) :: c1t
2202
2203! REAL(kind_phys) :: wfn_neg = 0.1
2204 REAL(kind_phys) :: wfn_neg = 0.15
2205! REAL(kind_phys) :: wfn_neg = 0.25
2206! REAL(kind_phys) :: wfn_neg = 0.30
2207! REAL(kind_phys) :: wfn_neg = 0.35
2208
2209 REAL(kind_phys) :: esat, tem
2210! REAL(kind_phys) :: esat, tem, rhs_h, rhs_q
2211!
2212! [INTERNAL FUNC]
2213 REAL(kind_phys) FPREC ! precipitation ratio in condensate
2214 REAL(kind_phys) FRICE ! ice ratio in cloud water
2215 REAL(kind_phys) Z ! altitude
2216 REAL(kind_phys) ZH ! scale height
2217 REAL(kind_phys) T ! temperature
2218!
2219 fprec(z,zh) = min(max(one-exp(-(z-precz0)/zh), zero), one)
2220 frice(t) = min(max((tsice-t)/(tsice-twice), zero), one)
2221!
2222! Note: iteration is not made to diagnose cloud ice for simplicity
2223!
2224 do i=ists,iens
2225 acwf(i) = zero
2226 gcyt(i) = zero
2227 gcht(i) = zero
2228 gcqt(i) = zero
2229 gclt(i) = zero
2230 gcit(i) = zero
2231 gtprt(i) = zero
2232 gcut(i) = zero
2233 gcvt(i) = zero
2234 gcwt(i) = zero
2235 enddo
2236 do k=1,kmax+1
2237 do i=ists,iens
2238 gprciz(i,k) = zero
2239 gsnwiz(i,k) = zero
2240 enddo
2241 enddo
2242 do k=1,kmax
2243 do i=ists,iens
2244 wcv(i,k) = unset_kind_phys
2245 gclm(i,k) = unset_kind_phys
2246 gcim(i,k) = unset_kind_phys
2247 enddo
2248 enddo
2249 do k=1,kmax
2250 do i=ists,iens
2251 acwfk(i,k) = unset_kind_phys
2252 acwfn(i,k) = unset_kind_phys
2253 gclz(i,k) = zero
2254 gciz(i,k) = zero
2255!
2256 gchmz(i,k) = zero
2257 gcwmz(i,k) = zero
2258 gcimz(i,k) = zero
2259 gcumz(i,k) = zero
2260 gcvmz(i,k) = zero
2261 gtprmz(i,k) = zero
2262!
2263 buoy(i,k) = unset_kind_phys
2264 buoym(i,k) = unset_kind_phys
2265 gcy(i,k) = unset_kind_phys
2266!
2267 gchm(i,k) = unset_kind_phys
2268 gcwm(i,k) = unset_kind_phys
2269 gctm(i,k) = unset_kind_phys
2270 gcqm(i,k) = unset_kind_phys
2271 gcum(i,k) = unset_kind_phys
2272 gcvm(i,k) = unset_kind_phys
2273 enddo
2274 enddo
2275 do i=ists,iens
2276 gcqmz(i) = zero
2277 wcm(i) = unset_kind_phys
2278 wcm_(i) = zero
2279 enddo
2280! tracers
2281 do n=ntrq,ntr
2282 do i=ists,iens
2283 gctrt(i,n) = zero
2284 enddo
2285 do k=1,kmax
2286 do i=ists,iens
2287 gctrm(i,k,n) = unset_kind_phys
2288 enddo
2289 enddo
2290 enddo
2291
2292! DO I=ISTS,IENS
2293! if (kb(i) > 0) then
2294! GDZMKB(I) = GDZM(I,KB(I)) ! cloud base height
2295! endif
2296! ENDDO
2297!
2298! < cloud base properties >
2299!
2300 DO i=ists,iens
2301 k = kb(i)
2302 if (k > 0) then
2303 gdzmkb(i) = gdzm(i,k) ! cloud base height
2304 gchm(i,k) = gchb(i)
2305 gcwm(i,k) = gcwb(i)
2306 wcm_(i) = wcb(i)
2307 gcum(i,k) = gcub(i)
2308 gcvm(i,k) = gcvb(i)
2309 do n = ntrq,ntr
2310 gctrm(i,k,n) = gctrb(i,n)
2311 enddo
2312!
2313 esat = min(gdpm(i,k), fpvs(gdtm(i,k)))
2314 gdqsm = min(epsv*esat/max(gdpm(i,k)+epsvm1*esat, 1.0e-10), 0.1)
2315 gdsm = cp*gdtm(i,k) + grav*gdzmkb(i) ! dse
2316 gdhsm = gdsm + el*gdqsm ! saturated mse
2317! FDQSM = FDQSAT(GDTM(I,K), GDQSM)
2318 tem = one / gdtm(i,k)
2319 fdqsm = gdqsm * tem * (fact1 + fact2*tem) ! calculate d(qs)/dT
2320!
2321 tem = one / (cp+el*fdqsm)
2322 dctm = (gchb(i) - gdhsm) * tem
2323 gcqm(i,k) = min(gdqsm + fdqsm*dctm, gcwm(i,k))
2324 gccm = max(gcwm(i,k)-gcqm(i,k), zero)
2325! GCTM(I,K) = GDT(I,K) + DCTM ! old
2326! GCTM(I,K) = (GCHB(I) - gdsm - el*gcqm(i,k)) * oneocp + dctm ! new
2327 gctm(i,k) = (gchb(i) - grav*gdzm(i,k) - el*gcqm(i,k)) * oneocp + dctm ! new
2328!
2329 gcim(i,k) = frice(gctm(i,k)) * gccm ! cloud base ice
2330 gclm(i,k) = max(gccm-gcim(i,k), zero) ! cloud base liquid
2331 gchm(i,k) = gchm(i,k) + emelt * (gcim(i,k)-gcib(i))
2332 dctm = (gchm(i,k) - gdhsm) * tem
2333! GCTM(I,K) = dctm + GDT(I,K) + gocp*gdzm(i,k) !DD old AW convert to DSE
2334 gctm(i,k) = dctm + (gchb(i) - el*gcqm(i,k)) * oneocp ! new, make dse
2335!
2336 gdqm = half * (gdq(i,k,1) + gdq(i,k-1,1))
2337 gdcm = half * (gdq(i,k,itl) + gdqi(i,k) &
2338 + gdq(i,k-1,itl) + gdqi(i,k-1))
2339
2340!
2341 buoym(i,k) = (dctm*tem + epsvt*(gcqm(i,k)-gdqm) - gccm + gdcm )*grav
2342!
2343 acwfk(i,k) = zero
2344 acwfn(i,k) = zero
2345!
2346!DD#ifdef OPT_ASMODE
2347!DD ELARM1(I) = ERMR
2348!DD#elif defined OPT_NS02
2349!DD ELARM1(I) = ELMD / SQRT(WCM(I,K))
2350!DD#else
2351! ELARM1(I) = CLMD*PA*BUOYM(I,K)/WCM(I,K)
2352! ELARM1(I) = min(max(CLMD*PA*BUOYM(I,K)/WCM_(I), ELAMIN), ELAMAX)
2353!DD#endif
2354! ELARM1(I) = MIN(MAX(ELARM1(I), ELAMIN), ELAMAX)
2355!
2356 gchmz(i,k) = gchm(i,k)
2357 gcwmz(i,k) = gcwm(i,k)
2358 gcumz(i,k) = gcum(i,k)
2359 gcvmz(i,k) = gcvm(i,k)
2360 gcqmz(i) = gcqm(i,k)
2361 gcimz(i,k) = gcim(i,k)
2362 do n = ntrq,ntr
2363 gctrmz(i,k,n) = gctrm(i,k,n)
2364 enddo
2365 endif
2366 ENDDO
2367!
2368! < in-cloud properties >
2369!
2370 DO k=3,kmax
2371 km1 = k - 1
2372 DO i=ists,iens
2373 IF (kb(i) > 0 .and. k > kb(i) .AND. wcm_(i) > wccrt) THEN
2374 wcv(i,km1) = sqrt(max(wcm_(i), zero))
2375 delz = gdzm(i,k) - gdzm(i,km1)
2376 elarm1(i) = min(max(clmdpa*buoym(i,km1)/wcm_(i), elamin), elamax)
2377 gcym(i,k) = gcym(i,km1) * exp(elarm1(i)*delz)
2378 eladz = gcym(i,k) - gcym(i,km1)
2379!
2380 gchmz(i,k) = gchmz(i,km1) + gdh(i,km1)*eladz
2381 gcwmz(i,k) = gcwmz(i,km1) + gdw(i,km1)*eladz
2382!
2383 esat = min(gdpm(i,k), fpvs(gdtm(i,k)))
2384 gdqsm = min(epsv*esat/max(gdpm(i,k)+epsvm1*esat, 1.0e-10), 0.1)
2385 gdhsm = cp*gdtm(i,k ) + grav*gdzm(i,k) + el*gdqsm
2386! FDQSM = FDQSAT(GDTM(I,K), GDQSM)
2387 tem = one / gdtm(i,k)
2388 fdqsm = gdqsm * tem * (fact1 + fact2*tem) ! calculate d(qs)/dT
2389 cpgmi = one / (cp + el*fdqsm)
2390
2391!
2392 wrk = one / gcym(i,k)
2393 dctm = (gchmz(i,k)*wrk - gdhsm) * cpgmi
2394 gcqmz(i) = min((gdqsm+fdqsm*dctm)*gcym(i,k), gcwmz(i,k))
2395 if(preczh > zero) then
2396 prczh = preczh * min(gdztr(i)*ztrefi, one)
2397 precr = fprec(gdzm(i,k)-gdzmkb(i), prczh )
2398 gtprmz(i,k) = precr * (gcwmz(i,k)-gcqmz(i))
2399 else
2400 delc=gdz(i,k)-gdz(i,km1)
2401 if(gdtm(i,k)>tsice) then
2402 c1t=c0t*delc
2403 else
2404 c1t=c0t*exp(d0t*(gdtm(i,k)-tsice))*delc
2405 end if
2406 c1t=min(c1t, one)
2407 gtprmz(i,k) = c1t * (gcwmz(i,k)-gcqmz(i))
2408 end if
2409 gtprmz(i,k) = max(gtprmz(i,k), gtprmz(i,km1))
2410 gccmz = gcwmz(i,k) - gcqmz(i) - gtprmz(i,k )
2411 delc = min(gccmz, zero)
2412 gccmz = gccmz - delc
2413 gcqmz(i) = gcqmz(i) + delc
2414!
2415 fice = frice(gdtm(i,k)+dctm )
2416 gcimz(i,k) = fice * gccmz
2417 gsnwiz(i,km1) = fice * (gtprmz(i,k)-gtprmz(i,km1))
2418 gchmz(i,k) = gchmz(i,k) + emelt * (gcimz(i,k) + gsnwiz(i,km1) &
2419 - gcimz(i,km1) - gdqi(i,km1)*eladz)
2420 dctm = (gchmz(i,k)*wrk - gdhsm) * cpgmi
2421!
2422 gdqm = half * (gdq(i,k,1) + gdq(i,km1,1))
2423 gdcm = half * (gdq(i,k,itl) + gdqi(i,k) &
2424 + gdq(i,km1,itl) + gdqi(i,km1))
2425 gcqm(i,k) = wrk * gcqmz(i)
2426 gccm = wrk * gccmz
2427!
2428 buoym(i,k) = (dctm*tem + epsvt*(gcqm(i,k)-gdqm)-gccm+gdcm) * grav
2429 buoy(i,km1) = half * (buoym(i,k)+buoym(i,km1))
2430!
2431!DD#ifdef OPT_ASMODE
2432!DD WCM(I,K ) &
2433!DD = (WCM_(I) + 2.D0*PA*DELZ*BUOY(I,KM1) ) &
2434!DD / (one + 2.D0*PB*DELZ*ERMR)
2435!DD#elif OPT_NS02
2436!DD WCM(I,K ) = WCM_(I ) &
2437!DD + 2.D0*DELZ*(PA*BUOYM(I,KM1)-ELMD*WCV(I,KM1))
2438!DD WCM(I,K ) = MAX(WCM(I,K ), zero )
2439!DD WCVX = SQRT(half*(WCM(I,K )+WCM_(I)))
2440!DD WCM(I,K) = WCM_(I) + 2.D0*DELZ*(PA*BUOY(I,KM1)-ELMD*WCVX)
2441!DD#else
2442 IF (buoy(i,km1) > zero) THEN
2443 wcm(i) = (wcm_(i) + clmp*delz*buoy(i,km1)) / (one + delz/tauz)
2444 ELSE
2445 wcm(i) = (wcm_(i) + pa*(delz+delz)*buoy(i,km1) ) &
2446 / (one + delz/tauz + (delz+delz)*elamin )
2447 ENDIF
2448!DD#endif
2449!
2450!DD#ifdef OPT_ASMODE
2451!DD ELARM2 = ERMR
2452!DD#elif OPT_NS02
2453!DD ELARM2 = ELMD/SQRT(MAX(WCM(I,K), EPSln))
2454!DD#else
2455! ELARM2 = CLMD*PA*BUOYM(I,K) / MAX(WCM(I), EPSln)
2456!DD#endif
2457 if (wcm(i) > zero) then
2458 elarm2 = min(max(clmdpa*buoym(i,k)/wcm(i),elamin), elamax)
2459 else
2460 elarm2 = zero
2461 endif
2462 elar = half * (elarm1(i) + elarm2)
2463 gcym(i,k) = gcym(i,km1) * exp(elar*delz)
2464 eladz = gcym(i,k) - gcym(i,km1)
2465!
2466 gchmz(i,k) = gchmz(i,km1) + gdh(i,km1)*eladz
2467 gcwmz(i,k) = gcwmz(i,km1) + gdw(i,km1)*eladz
2468 gcumz(i,k) = gcumz(i,km1) + gdu(i,km1)*eladz
2469 gcvmz(i,k) = gcvmz(i,km1) + gdv(i,km1)*eladz
2470 do n = ntrq,ntr
2471 gctrmz(i,k,n) = gctrmz(i,km1,n) + gdq(i,km1,n)*eladz
2472 enddo
2473!
2474 wrk = one / gcym(i,k)
2475 dctm = (gchmz(i,k)*wrk - gdhsm) * cpgmi
2476 gcqmz(i) = min((gdqsm+fdqsm*dctm)*gcym(i,k), gcwmz(i,k))
2477 if(preczh > zero) then
2478 gtprmz(i,k) = precr * (gcwmz(i,k)-gcqmz(i))
2479 else
2480 gtprmz(i,k) = c1t * (gcwmz(i,k)-gcqmz(i))
2481 end if
2482 gtprmz(i,k) = max(gtprmz(i,k), gtprmz(i,km1))
2483 gccmz = gcwmz(i,k) - gcqmz(i) - gtprmz(i,k)
2484 delc = min(gccmz, zero)
2485 gccmz = gccmz - delc
2486 gcqmz(i) = gcqmz(i) + delc
2487 gccm = wrk * gccmz
2488 gcqm(i,k) = wrk * gcqmz(i)
2489!
2490 fice = frice(gdtm(i,k)+dctm )
2491 gcimz(i,k) = fice*gccmz
2492 gcim(i,k) = gcimz(i,k)*wrk
2493 gclm(i,k) = max(gccm-gcim(i,k), zero)
2494 gtpriz = gtprmz(i,k) - gtprmz(i,km1)
2495 gsnwiz(i,km1) = fice*gtpriz
2496
2497 gprciz(i,km1) = (one-fice )*gtpriz
2498 gchmz(i,k) = gchmz(i,k) + emelt*(gcimz(i,k) + gsnwiz(i,km1) &
2499 - gcimz(i,km1) - gdqi(i,km1)*eladz )
2500 gchm(i,k) = gchmz(i,k)*wrk
2501 dctm = (gchm(i,k)-gdhsm) * cpgmi
2502! GCTM(I,K) = dctm + GDTM(I,K) + gocp*gdzm(i,k) ! old, make dse
2503 gctm(i,k) = dctm + (gchm(i,k) - el*gcqm(i,k)) * oneocp ! new, make dse
2504!
2505 gcwm(i,k) = gcwmz(i,k) * wrk
2506 gcum(i,k) = gcumz(i,k) * wrk
2507 gcvm(i,k) = gcvmz(i,k) * wrk
2508 do n = ntrq,ntr
2509 gctrm(i,k,n) = gctrmz(i,k,n) * wrk
2510 enddo
2511 delzl = gdz(i,km1)-gdzm(i,km1)
2512 gcy(i,km1) = gcym(i,km1) * exp(elar*delzl)
2513 gclz(i,km1) = half * (gclm(i,k) + gclm(i,km1)) * gcy(i,km1)
2514 gciz(i,km1) = half * (gcim(i,k) + gcim(i,km1)) * gcy(i,km1)
2515
2516!
2517 buoym(i,k) = (dctm*tem + epsvt*(gcqm(i,k)-gdqm)-gccm+gdcm) * grav
2518 buoy(i,km1) = half * (buoym(i,k)+buoym(i,km1))
2519!
2520 IF (buoy(i,km1) > zero) THEN
2521 wcm(i) = (wcm_(i) + clmp*delz*buoy(i,km1)) / (one + delz/tauz)
2522 ELSE
2523 wcm(i) = (wcm_(i) + pa*(delz+delz)*buoy(i,km1) ) &
2524 / (one + delz/tauz + (delz+delz)*elamin )
2525 ENDIF
2526
2527!
2528! IF (BUOY(I,KM1) > zero) THEN
2529! ACWF(I) = ACWF(I) + BUOY(I,KM1)*GCY(I,KM1)*DELZ
2530! ENDIF
2531! ACWF(I) = ACWF(I) + BUOY(I,KM1)*GCY(I,KM1)*DELZ
2532!!! wrk = BUOY(I,KM1)*GCY(I,KM1)*DELZ
2533!!! ACWFK(I,K) = ACWFK(I,KM1) + wrk
2534!!! ACWFN(I,K) = ACWFN(I,KM1) - min(wrk,0.0)
2535! ACWFN(I,K) = ACWFN(I,KM1) + abs(min(wrk,0.0))
2536!
2537
2538 wrk = buoy(i,km1)*gcy(i,km1)*delz
2539 acwfk(i,k) = acwfk(i,km1) + wrk
2540 acwfn(i,k) = acwfn(i,km1) - min(wrk,0.0)
2541
2542 wcm_(i) = wcm(i)
2543
2544! if (lprnt .and. i == ipr) write(0,*) ' in cumup k=',k,' km1=',km1,' WCM_=',WCM_(I),' gcy=',gcy(i,km1),' buoym=',buoym(i,km1)
2545
2546 ENDIF ! IF (K > KB(I) .AND. WCM_(I) > WCCRT) THEN
2547 ENDDO
2548 ENDDO
2549!
2550! < find cloud top >
2551!
2552 DO i=ists,iens
2553 kt(i) = -1
2554 enddo
2555 DO k=kmax,2,-1
2556 DO i=ists,iens
2557 if (kb(i) > 0 .and. k > kb(i) .and. acwfk(i,k) > 1.0e-10) then
2558 wrk = acwfn(i,k) / acwfk(i,k)
2559 IF (kt(i) == -1 .and. wrk < wfn_neg .AND. wcv(i,k) > wvcrt) THEN
2560 kt(i) = k
2561 acwf(i) = acwfk(i,k)
2562 ENDIF
2563 endif
2564 ENDDO
2565 ENDDO
2566! if (lprnt .and. kt(ipr) > 0) write(0,*) ' in cumup kt=',kt(ipr),' gcy=',gcy(ipr,kt(ipr))
2567!
2568 ktmx = 2
2569 DO i=ists,iens
2570 kt(i) = min(kt(i), kmax-1)
2571 ktmx = max(ktmx, kt(i))
2572 ENDDO
2573!
2574 DO i=ists,iens
2575 kk = max(1, kt(i)+1)
2576 do k=kk,kmax
2577 gcym(i,k) = zero
2578 gclz(i,k) = zero
2579 gciz(i,k) = zero
2580 gprciz(i,k) = zero
2581 gsnwiz(i,k) = zero
2582 enddo
2583 ENDDO
2584! if (lprnt .and. kt(ipr) > 0) write(0,*) ' in cumup2 kt=',kt(ipr),' gcy=',gcy(ipr,kt(ipr))
2585!
2586! < cloud top properties >
2587!
2588 DO i=ists,iens
2589 IF (kb(i) > 0 .and. kt(i) > kb(i)) THEN
2590 k = kt(i)
2591 kp1 = k + 1
2592 gcyt(i) = gcy(i,k)
2593 eladz = gcyt(i) - gcym(i,k)
2594!
2595 gcht(i) = gchmz(i,k) + gdh(i,k)*eladz
2596 gcwt(i) = gcwmz(i,k) + gdw(i,k)*eladz
2597 gcut(i) = gcumz(i,k) + gdu(i,k)*eladz
2598 gcvt(i) = gcvmz(i,k) + gdv(i,k)*eladz
2599 do n = ntrq,ntr
2600 gctrt(i,n) = gctrmz(i,k,n) + gdq(i,k,n)*eladz
2601 enddo
2602!
2603 wrk = one / gcyt(i)
2604 dct = (gcht(i)*wrk - gdhs(i,k)) / (cp*(one + gam(i,k)))
2605 gcqt(i) = min((gdqs(i,k) + fdqs(i,k)*dct) * gcyt(i), gcwt(i))
2606 if(preczh > zero) then
2607 prczh = preczh * min(gdztr(i)*ztrefi, one)
2608 gtprt(i) = fprec(gdz(i,k)-gdzmkb(i), prczh) * (gcwt(i)-gcqt(i))
2609 else
2610 delc=gdz(i,k)-gdz(i,k-1)
2611 if(gdtm(i,k)>tsice) then
2612 c1t=c0t*delc
2613 else
2614 c1t=c0t*exp(d0t*(gdtm(i,k)-tsice))*delc
2615 end if
2616 c1t=min(c1t, one)
2617 gtprt(i) = c1t * (gcwt(i)-gcqt(i))
2618 end if
2619 gtprt(i) = max(gtprt(i), gtprmz(i,k))
2620 gcct = gcwt(i) - gcqt(i) - gtprt(i)
2621 delc = min(gcct, zero)
2622 gcct = gcct - delc
2623 gcqt(i) = gcqt(i) + delc
2624!
2625 fice = frice(gdt(i,k)+dct)
2626 gcit(i) = fice*gcct
2627 gclt(i) = (one-fice) * gcct
2628 gtpriz = gtprt(i) - gtprmz(i,k)
2629 gprciz(i,k) = (one-fice) * gtpriz
2630 gsnwiz(i,k) = fice * gtpriz
2631 gcht(i) = gcht(i) &
2632 + emelt * (gcit(i) + gsnwiz(i,k) - gcimz(i,k) - gdqi(i,k)*eladz)
2633!
2634 gcut(i) = gcut(i)*(one-cpres) + gcy(i,k)*gdu(i,k)*cpres
2635 gcvt(i) = gcvt(i)*(one-cpres) + gcy(i,k)*gdv(i,k)*cpres
2636 do n = ntrq,ntr
2637! GCtrT(I,n) = GCtrT(I,n)*(one-CPRES) + GCY(I,K)*GDq(I,K,n)*CPRES
2638 gctrt(i,n) = gctrt(i,n) + gcy(i,k)*gdq(i,k,n)
2639 enddo
2640 gclz(i,k) = gclt(i)
2641 gciz(i,k) = gcit(i)
2642
2643!DD AW get the cloud top values denormalized and put into profile
2644 mygcht = gcht(i) - el*(gcwt(i) - gcqt(i))
2645
2646 gctm(i,kp1) = wrk * (mygcht - el*gcqt(i)) * oneocp
2647!Moorthi gcqm(i,kp1) = gcqt(i)
2648 gcqm(i,kp1) = gcqt(i)*wrk ! check this - oct17 2016
2649 gcim(i,kp1) = gcit(i)*wrk
2650 gclm(i,kp1) = gclt(i)*wrk
2651 do n = ntrq,ntr
2652 gctrm(i,kp1,n) = gctrt(i,n)*wrk
2653 enddo
2654!
2655 ENDIF
2656 ENDDO
2657!
2658!DD#ifdef OPT_CUMBGT /* budget check */
2659!DD HBGT ( ISTS:IENS ) = 0.D0
2660!DD WBGT ( ISTS:IENS ) = 0.D0
2661!DD PBGT ( ISTS:IENS ) = 0.D0
2662!DD MBGT ( ISTS:IENS ) = 0.D0
2663!DD GTPRX( ISTS:IENS ) = 0.D0
2664!DD GSNWT( ISTS:IENS ) = 0.D0
2665!DD!
2666!DD IF ( CTP .EQ. 1 ) THEN
2667!DD HBMX = 0.D0
2668!DD WBMX = 0.D0
2669!DD PBMX = 0.D0
2670!DD MBMX = 0.D0
2671!DD END IF
2672!DD!
2673!DD DO K = 2, KMAX
2674!DD DO I = ISTS, IENS
2675!DD IF ( K .GE. KB( I ) .AND. K .LT. KT( I ) ) THEN
2676!DD ELADZ = GCYM( I,K+1 ) - GCYM( I,K )
2677!DD DELZ = GDZM( I,K+1 ) - GDZM( I,K )
2678!DD HBGT( I ) = HBGT( I ) + ( GDH( I,K )-EMELT*GDQI( I,K ) )*ELADZ
2679!DD WBGT( I ) = WBGT( I ) + GDW( I,K )*ELADZ
2680!DD MBGT( I ) = MBGT( I ) + ELAM( I,K )*DELZ
2681!DD GTPRX( I ) = GTPRX( I ) + GPRCIZ( I,K ) + GSNWIZ( I,K )
2682!DD GSNWT( I ) = GSNWT( I ) + GSNWIZ( I,K )
2683!DD END IF
2684!DD END DO
2685!DD END DO
2686!DD!
2687!DD DO I = ISTS, IENS
2688!DD IF ( KT( I ) .GT. KB( I ) ) THEN
2689!DD ELADZ = GCYT( I ) - GCYM( I,KT(I) )
2690!DD DELZ = GDZ( I,KT(I) )-GDZM( I,KT(I) )
2691!DD GTPRX( I ) = GTPRX( I ) + GPRCIZ( I,KT(I) ) + GSNWIZ( I,KT(I) )
2692!DD GSNWT( I ) = GSNWT( I ) + GSNWIZ( I,KT(I) )
2693!DD HBGT( I ) = HBGT( I ) + GCHB( I ) - EMELT*GCIB( I ) &
2694!DD + ( GDH( I,KT(I) )-EMELT*GDQI( I,KT(I) ) ) *ELADZ &
2695!DD - ( GCHT(I)-EMELT*( GCIT(I)+GSNWT(I) ) )
2696!DD WBGT( I ) = WBGT( I ) &
2697!DD + GCWB( I ) + GDW( I,KT(I) )*ELADZ &
2698!DD - GCQT( I ) - GCLT( I ) - GCIT( I ) &
2699!DD - GTPRT( I )
2700!DD MBGT( I ) = MBGT( I ) + one + ELAM( I,KT(I) )*DELZ &
2701!DD - GCYT( I )
2702!DD PBGT( I ) = GTPRT( I ) - GTPRX( I )
2703!DD!
2704!DD IF ( ABS( HBGT(I) ) .GT. ABS( HBMX ) ) HBMX = HBGT(I)
2705!DD IF ( ABS( WBGT(I) ) .GT. ABS( WBMX ) ) WBMX = WBGT(I)
2706!DD IF ( ABS( PBGT(I) ) .GT. ABS( PBMX ) ) PBMX = PBGT(I)
2707!DD IF ( ABS( MBGT(I) ) .GT. ABS( MBMX ) ) MBMX = MBGT(I)
2708!DD END IF
2709!DD END DO
2710!DD!
2711!DD IF ( CTP .EQ. NCTP ) THEN
2712!DD WRITE( iulog,* ) '### CUMUP(rank=',irank,'): energy imbalance =', HBMX
2713!DD WRITE( iulog,* ) '### CUMUP(rank=',irank,'): water imbalance =', WBMX
2714!DD WRITE( iulog,* ) '### CUMUP(rank=',irank,'): precipitation imbalance =', PBMX
2715!DD WRITE( iulog,* ) '### CUMUP(rank=',irank,'): mass imbalance =', MBMX
2716!DD END IF
2717!DD#endif
2718!
2719! WRITE( CTNUM, '(I2.2)' ) CTP
2720!
2721 END SUBROUTINE cumup
2722!***********************************************************************
2725 SUBROUTINE cumbmx & !! cloud base mass flux
2726 ( ijsdim, kmax , & !DD dimensions
2727 cbmfx , & ! modified
2728 acwf , gcyt , gdzm , & ! input
2729 gdw , gdqs , delp , & ! input
2730 kt , ktmx , kb , & ! input
2731 delt , ists , iens ) ! input
2732!
2733!
2734 IMPLICIT NONE
2735
2736 INTEGER, INTENT(IN) :: IJSDIM, KMAX ! DD, for GFS, pass in
2737!
2738! [MODIFY]
2739 REAL(kind_phys) CBMFX (IJSDIM)
2740!
2741! [INPUT]
2742 REAL(kind_phys) ACWF (IJSDIM)
2743 REAL(kind_phys) GCYT (IJSDIM)
2744 REAL(kind_phys) GDZM (IJSDIM, KMAX+1)
2745 REAL(kind_phys) GDW (IJSDIM, KMAX)
2746 REAL(kind_phys) GDQS (IJSDIM, KMAX)
2747 REAL(kind_phys) DELP (IJSDIM, KMAX)
2748 INTEGER KT (IJSDIM)
2749 INTEGER KTMX
2750 INTEGER KB (IJSDIM)
2751 REAL(kind_phys) DELT
2752 INTEGER ISTS, IENS
2753!
2754! [INTERNAL WORK]
2755 REAL(kind_phys), dimension(ijsdim) :: QX, QSX, RHM
2756 INTEGER I, K
2757 REAL(kind_phys) ALP, FMAX1, wrk
2758!
2759! [INTERNAL PARAM]
2760 REAL(kind_phys) :: FMAX = 1.5e-2_kind_phys ! maximum flux
2761! REAL(kind_phys) :: RHMCRT = zero ! critical val. of cloud mean RH
2762 REAL(kind_phys) :: RHMCRT = 0.25_kind_phys ! critical val. of cloud mean RH
2763! REAL(kind_phys) :: RHMCRT = 0.50_kind_phys ! critical val. of cloud mean RH
2764 REAL(kind_phys) :: ALP1 = zero
2765 REAL(kind_phys) :: TAUD = 1.e3_kind_phys
2766! REAL(kind_phys) :: TAUD = 6.e2_kind_phys
2767 REAL(kind_phys) :: ZFMAX = 3.5e3_kind_phys
2768 REAL(kind_phys) :: ZDFMAX = 5.e2_kind_phys
2769! REAL(kind_phys) :: FMAXP = 2._kind_phys
2770!
2771 do i=ists,iens
2772 qx(i) = zero
2773 qsx(i) = zero
2774 enddo
2775!
2776 DO k=1,ktmx
2777 DO i=ists,iens
2778 IF (kb(i) > 0 .and. k >= kb(i) .AND. k <= kt(i)) THEN
2779 qx(i) = qx(i) + gdw(i,k) * delp(i,k)
2780 qsx(i) = qsx(i) + gdqs(i,k) * delp(i,k)
2781 ENDIF
2782 ENDDO
2783 ENDDO
2784 DO i=ists,iens
2785 rhm(i) = min(one, max(zero, qx(i)/max(qsx(i),epsln)))
2786 ENDDO
2787!
2788 wrk = one + delt/(taud+taud)
2789 DO i=ists,iens
2790 k = kb(i)
2791 IF (k > 0 .and. kt(i) > k .AND. rhm(i) >= rhmcrt) THEN
2792 cbmfx(i) = max(cbmfx(i), zero)
2793 alp = alp0 + alp1 * (gdzm(i,kt(i))-gdzm(i,k))
2794 fmax1 = (one - tanh((gdzm(i,1)-zfmax)/zdfmax)) * half
2795! FMAX1 = FMAX * FMAX1**FMAXP
2796 fmax1 = fmax * fmax1 * fmax1
2797! CBMFX(I) = CBMFX(I) + MAX(ACWF(I), zero)/(ALP+ALP)*DELT
2798! CBMFX(I) = CBMFX(I) / (one + DELT/(TAUD+TAUD))
2799 cbmfx(i) = (cbmfx(i) + acwf(i)*(delt/(alp+alp))) * wrk
2800 cbmfx(i) = min(max(cbmfx(i), zero), fmax1/gcyt(i))
2801 ELSE
2802 cbmfx(i) = zero
2803 ENDIF
2804 ENDDO
2805!
2806 END SUBROUTINE cumbmx
2807!***********************************************************************
2810 SUBROUTINE cumflx & !! cloud mass flux
2811 ( im , ijsdim, kmax , & !DD dimensions
2812 gmflx , gprci , gsnwi , cmdet, & ! output
2813 qliq , qice , gtprc0, & ! output
2814 cbmfx , gcym , gprciz, gsnwiz, & ! input
2815 gtprt , gclz , gciz , gcyt, & ! input
2816 kb , kt , ktmx , & ! input
2817 ists , iens ) ! input
2818!
2819 IMPLICIT NONE
2820
2821 INTEGER, INTENT(IN) :: IJSDIM, KMAX, IM !! DD, for GFS, pass in
2822!
2823! [OUTPUT]
2824 REAL(kind_phys) GMFLX (IJSDIM, KMAX+1)
2825 REAL(kind_phys) CMDET (IJSDIM, KMAX)
2826 REAL(kind_phys) GPRCI (IJSDIM, KMAX)
2827 REAL(kind_phys) GSNWI (IJSDIM, KMAX)
2828 REAL(kind_phys) QLIQ (IJSDIM, KMAX)
2829 REAL(kind_phys) QICE (IJSDIM, KMAX)
2830 REAL(kind_phys) GTPRC0(IJSDIM)
2831!
2832! [INPUT]
2833 REAL(kind_phys) CBMFX (IJSDIM)
2834 REAL(kind_phys) GCYM (IJSDIM, KMAX)
2835 REAL(kind_phys) GCYT (IJSDIM)
2836 REAL(kind_phys) GPRCIZ(IJSDIM, KMAX+1)
2837 REAL(kind_phys) GSNWIZ(IJSDIM, KMAX+1)
2838 REAL(kind_phys) GTPRT (IJSDIM)
2839 REAL(kind_phys) GCLZ (IJSDIM, KMAX)
2840 REAL(kind_phys) GCIZ (IJSDIM, KMAX)
2841 INTEGER KB (IJSDIM)
2842 INTEGER KT (IJSDIM)
2843 INTEGER KTMX
2844 INTEGER ISTS, IENS, I, K
2845!
2846 DO k=1,ktmx
2847 DO i=ists,iens
2848 if (kb(i) > 0) then
2849 gmflx(i,k) = gmflx(i,k) + cbmfx(i) * gcym(i,k)
2850 gprci(i,k) = gprci(i,k) + cbmfx(i) * gprciz(i,k)
2851 gsnwi(i,k) = gsnwi(i,k) + cbmfx(i) * gsnwiz(i,k)
2852 qliq(i,k) = qliq(i,k) + cbmfx(i) * gclz(i,k)
2853 qice(i,k) = qice(i,k) + cbmfx(i) * gciz(i,k)
2854 endif
2855 ENDDO
2856 ENDDO
2857!
2858 DO i= ists,iens
2859 k = kt(i)
2860 if (kb(i) > 0 .and. k > kb(i)) then
2861 gtprc0(i) = gtprc0(i) + cbmfx(i) * gtprt(i)
2862 cmdet(i,k) = cmdet(i,k) + cbmfx(i) * gcyt(i)
2863 endif
2864 ENDDO
2865!
2866 END SUBROUTINE cumflx
2867!***********************************************************************
2870 SUBROUTINE cumdet & !! detrainment
2871 ( im , ijsdim, kmax , ntr , ntrq , & !DD dimensions
2872 gtt , gtq , gtu , gtv , & ! modified
2873 gdh , gdq , gdu , gdv , & ! input
2874! GTT , GTQ , GTCFRC, GTU , GTV , & ! modified
2875! GDH , GDQ , GDCFRC, GDU , GDV , & ! input
2876 cbmfx , gcyt , delpi , gcht , gcqt , & ! input
2877 gclt , gcit , gcut , gcvt , gdqi , & ! input
2878 gctrt, &
2879 kt , ists , iens , nctp ) ! input
2880!
2881 IMPLICIT NONE
2882
2883 INTEGER, INTENT(IN) :: im, IJSDIM, KMAX, NTR, nctp, ntrq !! DD, for GFS, pass in
2884!
2885! [MODIFY]
2886 REAL(kind_phys) GTT (IJSDIM, KMAX)
2887 REAL(kind_phys) GTQ (IJSDIM, KMAX, NTR)
2888! REAL(kind_phys) GTCFRC(IJSDIM, KMAX) !< cloud fraction tendency
2889 REAL(kind_phys) GTU (IJSDIM, KMAX)
2890 REAL(kind_phys) GTV (IJSDIM, KMAX)
2891!
2892! [INPUT]
2893 REAL(kind_phys) GDH (IJSDIM, KMAX)
2894 REAL(kind_phys) GDQ (IJSDIM, KMAX, NTR)
2895! REAL(kind_phys) GDCFRC(IJSDIM, KMAX) !< cloud fraction
2896 REAL(kind_phys) GDU (IJSDIM, KMAX)
2897 REAL(kind_phys) GDV (IJSDIM, KMAX)
2898 REAL(kind_phys) DELPI (IJSDIM, KMAX)
2899 REAL(kind_phys) CBMFX (IM, NCTP)
2900 REAL(kind_phys) GCYT (IJSDIM, NCTP)
2901 REAL(kind_phys) GCHT (IJSDIM, NCTP)
2902 REAL(kind_phys) GCQT (IJSDIM, NCTP)
2903 REAL(kind_phys) GCLT (IJSDIM, NCTP)
2904 REAL(kind_phys) GCIT (IJSDIM, NCTP)
2905 REAL(kind_phys) GCtrT (IJSDIM, ntrq:ntr, NCTP)
2906 REAL(kind_phys) GCUT (IJSDIM, NCTP)
2907 REAL(kind_phys) GCVT (IJSDIM, NCTP)
2908 REAL(kind_phys) GDQI (IJSDIM, KMAX)
2909 INTEGER KT (IJSDIM, NCTP)
2910 INTEGER ISTS, IENS
2911!
2912! [INTERNAL WORK]
2913 REAL(kind_phys) GTHCI, GTQVCI, GTXCI
2914 integer I, K, CTP, kk,n
2915!
2916
2917 DO ctp=1,nctp
2918 DO i=ists,iens
2919 k = kt(i,ctp)
2920 IF (k > 0) THEN
2921 gtxci = delpi(i,k)*cbmfx(i,ctp)
2922
2923 gthci = gtxci * (gcht(i,ctp) - gcyt(i,ctp)*gdh(i,k))
2924 gtqvci = gtxci * (gcqt(i,ctp) - gcyt(i,ctp)*gdq(i,k,1))
2925!
2926 gtq(i,k,1) = gtq(i,k,1) + gtqvci
2927 gtt(i,k) = gtt(i,k) + (gthci - el*gtqvci) * oneocp
2928! ql tendency by detrainment is treated by stratiform scheme
2929 gtq(i,k,itl) = gtq(i,k,itl) + gtxci * (gclt(i,ctp) - gcyt(i,ctp)*gdq(i,k,itl))
2930! qi tendency by detrainment is treated by stratiform scheme
2931 gtq(i,k,iti) = gtq(i,k,iti) + gtxci * (gcit(i,ctp) - gcyt(i,ctp)*gdqi(i,k))
2932 do n = ntrq,ntr
2933 gtq(i,k,n) = gtq(i,k,n) + gtxci * (gctrt(i,n,ctp) - gcyt(i,ctp)*gdq(i,k,n))
2934 enddo
2935
2936! GTCFRC(I,K) = GTCFRC(I,K) + GTXCI * (GCYT(I,CTP) - GCYT(I,CTP)*GDCFRC(I,K))
2937 gtu(i,k) = gtu(i,k) + gtxci * (gcut(i,ctp) - gcyt(i,ctp)*gdu(i,k))
2938 gtv(i,k) = gtv(i,k) + gtxci * (gcvt(i,ctp) - gcyt(i,ctp)*gdv(i,k))
2939 ENDIF
2940 ENDDO
2941 ENDDO
2942!
2943 END SUBROUTINE cumdet
2944!***********************************************************************
2946 SUBROUTINE cumsbh & !! adiabat. descent
2947 ( im , ijsdim, kmax , ntr, ntrq, & !DD dimensions
2948 gtt , gtq , & ! modified
2949 gtu , gtv , & ! modified
2950 gdh , gdq , gdqi , & ! input
2951 gdu , gdv , & ! input
2952 delpi , gmflx , gmfx0 , & ! input
2953 ktmx , cpres , kb, ists , iens ) ! input
2954!
2955!
2956 IMPLICIT NONE
2957
2958 INTEGER, INTENT(IN) :: IJSDIM, IM, KMAX, NTR, ntrq !! DD, for GFS, pass in
2959!
2960! [MODIFY]
2961 REAL(kind_phys) GTT (IJSDIM, KMAX)
2962 REAL(kind_phys) GTQ (IJSDIM, KMAX, NTR)
2963 REAL(kind_phys) GTU (IJSDIM, KMAX)
2964 REAL(kind_phys) GTV (IJSDIM, KMAX)
2965!
2966! [INPUT]
2967 REAL(kind_phys) GDH (IJSDIM, KMAX)
2968 REAL(kind_phys) GDQ (IJSDIM, KMAX, NTR)
2969 REAL(kind_phys) GDQI (IJSDIM, KMAX)
2970 REAL(kind_phys) GDU (IJSDIM, KMAX)
2971 REAL(kind_phys) GDV (IJSDIM, KMAX)
2972 REAL(kind_phys) DELPI (IJSDIM, KMAX)
2973 REAL(kind_phys) GMFLX (IJSDIM, KMAX+1)
2974 REAL(kind_phys) GMFX0 (IJSDIM, KMAX+1)
2975 INTEGER KB(IJSDIM)
2976 INTEGER KTMX
2977 REAL(kind_phys) CPRES
2978 INTEGER ISTS, IENS
2979!
2980! [INTERNAL WORK]
2981 REAL(kind_phys) SBH0, SBQ0, SBL0, SBI0, SBC0, SBS0, &
2982 SBH1, SBQ1, SBL1, SBI1, SBC1, SBS1, FX1, &
2983 SBU0, SBV0, SBU1, SBV1, GTHCI, GTQVCI, &
2984 GTQLCI, GTQICI, GTM2CI, GTM3CI, wrk, wrk1
2985 REAL(kind_phys) FX(ISTS:IENS)
2986
2987 REAL(kind_phys), dimension(IJSDIM, KMAX) :: GTLSBH, GTISBH
2988 integer :: I, K, KM, KP, n
2989!
2990!
2991 fx = zero
2992 do k=1,kmax
2993 do i=ists,iens
2994 gtlsbh(i,k) = zero
2995 gtisbh(i,k) = zero
2996 enddo
2997 enddo
2998!
2999 DO k=ktmx,1,-1
3000 km = max(k-1, 1)
3001 kp = min(k+1, kmax)
3002 DO i=ists,iens
3003 if (kb(i) > 0) then
3004 sbh0 = gmflx(i,kp) * (gdh(i,kp)-gdh(i,k))
3005 sbq0 = gmflx(i,kp) * (gdq(i,kp,1)-gdq(i,k,1))
3006 sbl0 = gmflx(i,kp) * (gdq(i,kp,itl )-gdq(i,k,itl))
3007 sbi0 = gmflx(i,kp) * (gdqi(i,kp)-gdqi(i,k))
3008 sbu0 = gmflx(i,kp) * (gdu(i,kp)-gdu(i,k)) &
3009 - gmfx0(i,kp) * (gdu(i,kp)-gdu(i,k))*cpres
3010 sbv0 = gmflx(i,kp) * (gdv(i,kp)-gdv(i,k)) &
3011 - gmfx0(i,kp) * (gdv(i,kp)-gdv(i,k))*cpres
3012!
3013 sbh1 = gmflx(i,k) * (gdh(i,k)-gdh(i,km))
3014 sbq1 = gmflx(i,k) * (gdq(i,k,1)-gdq(i,km,1))
3015 sbl1 = gmflx(i,k) * (gdq(i,k,itl)-gdq(i,km,itl))
3016 sbi1 = gmflx(i,k) * (gdqi(i,k)-gdqi(i,km))
3017 sbu1 = gmflx(i,k) * (gdu(i,k)-gdu(i,km)) &
3018 - gmfx0(i,k) * (gdu(i,k)-gdu(i,km))*cpres
3019 sbv1 = gmflx(i,k) * (gdv(i,k)-gdv(i,km)) &
3020 - gmfx0(i,k) * (gdv(i,k)-gdv(i,km))*cpres
3021!
3022 IF (gmflx(i,k) > gmflx(i,kp)) THEN
3023 fx1 = half
3024 ELSE
3025 fx1 = zero
3026 ENDIF
3027!
3028 wrk = delpi(i,k)
3029 wrk1 = one - fx(i)
3030 gthci = wrk * (wrk1*sbh0 + fx1 *sbh1)
3031 gtqvci = wrk * (wrk1*sbq0 + fx1 *sbq1)
3032 gtqlci = wrk * (wrk1*sbl0 + fx1 *sbl1)
3033 gtqici = wrk * (wrk1*sbi0 + fx1 *sbi1)
3034!
3035 gtt(i,k ) = gtt(i,k) + (gthci-el*gtqvci)*oneocp
3036 gtq(i,k,1 ) = gtq(i,k,1) + gtqvci
3037 gtq(i,k,itl) = gtq(i,k,itl) + gtqlci
3038 gtq(i,k,iti) = gtq(i,k,iti) + gtqici
3039 gtu(i,k) = gtu(i,k) + wrk * (wrk1*sbu0 + fx1*sbu1)
3040 gtv(i,k) = gtv(i,k) + wrk * (wrk1*sbv0 + fx1*sbv1)
3041 DO n = ntrq, ntr
3042 gtq(i,k,n) = gtq(i,k,n) + wrk &
3043 * ( wrk1 * (gmflx(i,kp) * (gdq(i,kp,n)-gdq(i,k ,n))) &
3044 + fx1 * (gmflx(i,k ) * (gdq(i,k ,n)-gdq(i,km,n))) )
3045 ENDDO
3046
3047 gtlsbh(i,k) = gtqlci
3048 gtisbh(i,k) = gtqici
3049!
3050! SBC0 = GMFLX(I,K+1) * (GDQ(I,KP,IMU2)-GDQ(I,K,IMU2))
3051! SBS0 = GMFLX(I,K+1) * (GDQ(I,KP,IMU3)-GDQ(I,K,IMU3))
3052! SBC1 = GMFLX(I,K) * (GDQ(I,K,IMU2)-GDQ(I,KM,IMU2))
3053! SBS1 = GMFLX(I,K) * (GDQ(I,K,IMU3)-GDQ(I,KM,IMU3))
3054! GTM2CI = DELPI(I,K) * (( one-FX(I))*SBC0 + FX1 *SBC1)
3055! GTM3CI = DELPI(I,K) * ((one-FX(I))*SBS0 + FX1 *SBS1)
3056! GTQ(I,K,IMU2) = GTQ(I,K,IMU2) + GTM2CI
3057! GTQ(I,K,IMU3) = GTQ(I,K,IMU3) + GTM3CI
3058!
3059 fx(i) = fx1
3060 endif
3061 enddo
3062 enddo
3063!
3064 END SUBROUTINE cumsbh
3065!***********************************************************************
3066!
3067!***********************************************************************
3070 SUBROUTINE cumsbw & !! adiabat. descent
3071 ( im , ijsdim, kmax , & !DD dimensions
3072 gtu , gtv , & ! modified
3073 gdu , gdv , & ! input
3074 delpi , gmflx , gmfx0 , & ! input
3075 ktmx , cpres , kb, ists , iens ) ! input
3076!
3077!
3078 IMPLICIT NONE
3079
3080 INTEGER, INTENT(IN) :: IJSDIM, IM, KMAX!! DD, for GFS, pass in
3081!
3082! [MODIFY]
3083 REAL(kind_phys) GTU (IJSDIM, KMAX)
3084 REAL(kind_phys) GTV (IJSDIM, KMAX)
3085!
3086! [INPUT]
3087 REAL(kind_phys) GDU (IJSDIM, KMAX)
3088 REAL(kind_phys) GDV (IJSDIM, KMAX)
3089 REAL(kind_phys) DELPI (IJSDIM, KMAX)
3090 REAL(kind_phys) GMFLX (IJSDIM, KMAX+1)
3091 REAL(kind_phys) GMFX0 (IJSDIM, KMAX+1)
3092 INTEGER KB(IJSDIM)
3093 INTEGER KTMX, ISTS, IENS
3094 REAL(kind_phys) CPRES
3095!
3096! [INTERNAL WORK]
3097 REAL(kind_phys) FX1, SBU0, SBV0, SBU1, SBV1, wrk, wrk1
3098 REAL(kind_phys) FX(ISTS:IENS)
3099
3100 integer :: I, K, KM, KP
3101!
3102!
3103 fx = zero
3104!
3105 DO k=ktmx,1,-1
3106 km = max(k-1, 1)
3107 kp = min(k+1, kmax)
3108 DO i=ists,iens
3109 if (kb(i) > 0) then
3110 sbu0 = gmflx(i,kp) * (gdu(i,kp)-gdu(i,k)) &
3111 - gmfx0(i,kp) * (gdu(i,kp)-gdu(i,k))*cpres
3112 sbv0 = gmflx(i,kp) * (gdv(i,kp)-gdv(i,k)) &
3113 - gmfx0(i,kp) * (gdv(i,kp)-gdv(i,k))*cpres
3114!
3115 sbu1 = gmflx(i,k) * (gdu(i,k)-gdu(i,km)) &
3116 - gmfx0(i,k) * (gdu(i,k)-gdu(i,km))*cpres
3117 sbv1 = gmflx(i,k) * (gdv(i,k)-gdv(i,km)) &
3118 - gmfx0(i,k) * (gdv(i,k)-gdv(i,km))*cpres
3119!
3120 IF (gmflx(i,k) > gmflx(i,kp)) THEN
3121 fx1 = half
3122 ELSE
3123 fx1 = zero
3124 ENDIF
3125!
3126 wrk = delpi(i,k)
3127 wrk1 = one - fx(i)
3128!
3129 gtu(i,k) = gtu(i,k) + wrk * (wrk1*sbu0 + fx1*sbu1)
3130 gtv(i,k) = gtv(i,k) + wrk * (wrk1*sbv0 + fx1*sbv1)
3131!
3132 fx(i) = fx1
3133 endif
3134 enddo
3135 enddo
3136!
3137 END SUBROUTINE cumsbw
3138!***********************************************************************
3141 SUBROUTINE cumdwn & ! Freeze & Melt & Evaporation
3142 ( im , ijsdim, kmax , ntr,ntrq,nctp, & !DD dimensions
3143 gtt , gtq , gtu , gtv , & ! modified
3144 gmflx , & ! modified
3145 gprcp , gsnwp , gtevp , gmdd , & ! output
3146 gprci , gsnwi , & ! input
3147 gdh , gdw , gdq , gdqi , & ! input
3148 gdqs , gds , gdhs , gdt , & ! input
3149 gdu , gdv , gdz , & ! input
3150 gdzm , fdqs , delp , & ! input
3151 delpi , &
3152 sigmad, do_aw , do_awdd, flx_form, & !DDsigma input
3153 gtmelt, gtevap, gtsubl, & !DDsigma input
3154 dtdwn , dqvdwn, dqldwn, dqidwn, & !DDsigma input
3155 dtrdwn, &
3156 kb , ktmx , ists , iens ) ! input
3157!
3158! DD AW : modify to get eddy fluxes and microphysical tendencies for AW
3159!
3160 IMPLICIT NONE
3161
3162 INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR , ntrq, nctp !! DD, for GFS, pass in
3163 logical, intent(in) :: do_aw, do_awdd, flx_form
3164!
3165! [MODIFY]
3166 REAL(kind_phys) GTT (IJSDIM, KMAX)
3167 REAL(kind_phys) GTQ (IJSDIM, KMAX, NTR)
3168 REAL(kind_phys) GTU (IJSDIM, KMAX)
3169 REAL(kind_phys) GTV (IJSDIM, KMAX)
3170 REAL(kind_phys) GMFLX (IJSDIM, KMAX+1)
3171!
3172! [OUTPUT]
3173 REAL(kind_phys) GPRCP (IJSDIM, KMAX+1)
3174 REAL(kind_phys) GSNWP (IJSDIM, KMAX+1)
3175 REAL(kind_phys) GTEVP (IJSDIM, KMAX)
3176 REAL(kind_phys) GMDD (IJSDIM, KMAX+1)
3177
3178!AW microphysical tendencies
3179 REAL(kind_phys) gtmelt (IJSDIM, KMAX)
3180 REAL(kind_phys) gtevap (IJSDIM, KMAX)
3181 REAL(kind_phys) gtsubl (IJSDIM, KMAX)
3182!AW eddy flux tendencies
3183 REAL(kind_phys) dtdwn (IJSDIM, KMAX)
3184 REAL(kind_phys) dqvdwn (IJSDIM, KMAX)
3185 REAL(kind_phys) dqldwn (IJSDIM, KMAX)
3186 REAL(kind_phys) dqidwn (IJSDIM, KMAX)
3187 REAL(kind_phys) dtrdwn (IJSDIM, KMAX, ntrq:ntr)
3188
3189! [INPUT]
3190 REAL(kind_phys) GPRCI (IJSDIM, KMAX)
3191 REAL(kind_phys) GSNWI (IJSDIM, KMAX)
3192 REAL(kind_phys) GDH (IJSDIM, KMAX)
3193 REAL(kind_phys) GDW (IJSDIM, KMAX)
3194 REAL(kind_phys) GDQ (IJSDIM, KMAX, NTR)
3195 REAL(kind_phys) GDQI (IJSDIM, KMAX)
3196 REAL(kind_phys) GDQS (IJSDIM, KMAX)
3197 REAL(kind_phys) GDS (IJSDIM, KMAX)
3198 REAL(kind_phys) GDHS (IJSDIM, KMAX)
3199 REAL(kind_phys) GDT (IJSDIM, KMAX)
3200 REAL(kind_phys) GDU (IJSDIM, KMAX)
3201 REAL(kind_phys) GDV (IJSDIM, KMAX)
3202 REAL(kind_phys) GDZ (IJSDIM, KMAX)
3203 REAL(kind_phys) GDZM (IJSDIM, KMAX+1)
3204 REAL(kind_phys) FDQS (IJSDIM, KMAX)
3205 REAL(kind_phys) DELP (IJSDIM, KMAX)
3206 REAL(kind_phys) DELPI (IJSDIM, KMAX)
3207 INTEGER KB (IJSDIM)
3208 INTEGER KTMX, ISTS, IENS
3209 REAL(kind_phys) sigmad (IM,KMAX+1)
3210
3211!
3212! [INTERNAL WORK]
3213! Note: Some variables have 3-dimensions for the purpose of budget check.
3214 REAL(kind_phys) EVAPD (IJSDIM, KMAX)
3215 REAL(kind_phys) SUBLD (IJSDIM, KMAX)
3216 REAL(kind_phys) EVAPE (IJSDIM, KMAX)
3217 REAL(kind_phys) SUBLE (IJSDIM, KMAX)
3218 REAL(kind_phys) EVAPX (IJSDIM, KMAX)
3219 REAL(kind_phys) SUBLX (IJSDIM, KMAX)
3220 REAL(kind_phys) GMDDE (IJSDIM, KMAX)
3221 REAL(kind_phys) SNMLT (IJSDIM, KMAX)
3222 REAL(kind_phys) GCHDD (IJSDIM, KMAX)
3223 REAL(kind_phys) GCWDD (IJSDIM, KMAX)
3224 REAL(kind_phys) GTTEV (IJSDIM, KMAX)
3225 REAL(kind_phys) GTQEV (IJSDIM, KMAX)
3226 REAL(kind_phys) GCHD (ISTS:IENS)
3227 REAL(kind_phys) GCWD (ISTS:IENS)
3228! profiles of downdraft variables for AW flux tendencies
3229 REAL(kind_phys) GCdseD(ISTS:IENS, KMAX)
3230 REAL(kind_phys) GCqvD (ISTS:IENS, KMAX)
3231 REAL(kind_phys) GCqlD (ISTS:IENS, KMAX)
3232 REAL(kind_phys) GCqiD (ISTS:IENS, KMAX)
3233 REAL(kind_phys) GCtrD (ISTS:IENS, ntrq:ntr)
3234
3235 REAL(kind_phys) GCUD (ISTS:IENS)
3236 REAL(kind_phys) GCVD (ISTS:IENS)
3237 REAL(kind_phys) FSNOW (ISTS:IENS)
3238 REAL(kind_phys) GMDDD (ISTS:IENS)
3239 INTEGER I, K
3240 REAL(kind_phys) GDTW
3241 REAL(kind_phys) GCHX, GCTX, GCQSX, GTPRP, EVSU, GTEVE, LVIC
3242 REAL(kind_phys) DQW, DTW, GDQW, DZ, GCSD, FDET, GDHI
3243 REAL(kind_phys) GMDDX, GMDDMX
3244 REAL(kind_phys) GCHDX, GCWDX
3245 REAL(kind_phys) GCUDD, GCVDD
3246 REAL(kind_phys) GTHCI, GTQVCI, GTQLCI, GTQICI
3247!M REAL(kind_phys) GTHCI, GTQVCI, GTQLCI, GTQICI, GTUCI, GTVCI
3248 real(kind_phys) wrk, fmelt, fevp, gctrdd(ntrq:ntr)
3249!DD#ifdef OPT_CUMBGT
3250 REAL(kind_phys) WBGT ( ISTS:IENS ) !! water budget
3251 REAL(kind_phys) HBGT ( ISTS:IENS ) !! energy budget
3252 REAL(kind_phys) DDWBGT( ISTS:IENS ) !! downdraft water budget
3253 REAL(kind_phys) DDHBGT( ISTS:IENS ) !! downdraft energy budget
3254 REAL(kind_phys) WMX, HMX, DDWMX, DDHMX, tx1, wrk1, wrk2, wrk3, wrk4, wrkn
3255 REAL(kind_phys) dp_above, dp_below
3256 real(kind_phys) fsigma
3257 integer ij, n, kp1
3258!DD#endif
3259!
3260! [INTERNAL PARM]
3261 REAL(kind_phys), parameter :: TWSNOW = 273.15_kind_phys
3262 REAL(kind_phys), parameter :: FTMLT = 4._kind_phys
3263 REAL(kind_phys), parameter :: GMFLXC = 5.e-2_kind_phys
3264 REAL(kind_phys), parameter :: VTERMS = 2._kind_phys
3265! REAL(kind_phys), parameter :: MELTAU = 10._kind_phys !< melting timescale
3266 REAL(kind_phys), parameter :: MELTAU = 20._kind_phys
3267!
3268! REAL(kind_phys), parameter :: EVAPR = 0.4_kind_phys !< evaporation factor ! Moorthi June 28, 2017
3269 REAL(kind_phys), parameter :: EVAPR = 0.3_kind_phys
3270! REAL(kind_phys), parameter :: EVAPR = 0._kind_phys !< evaporation factor
3271 REAL(kind_phys), parameter :: REVPDD = 1._kind_phys
3272 REAL(kind_phys), parameter :: RDDR = 5.e-4_kind_phys
3273! REAL(kind_phys), parameter :: RDDR = 0._kind_phys !< DD rate (T0 R0 W0)^-1
3274 REAL(kind_phys), parameter :: RDDMX = 0.5_kind_phys
3275 REAL(kind_phys), parameter :: VTERM = 5._kind_phys
3276! REAL(kind_phys), parameter :: VTERM = 4._kind_phys !< term. vel. of precip. ! Moorthi June 28, 2017
3277 REAL(kind_phys), parameter :: EVATAU = 2._kind_phys
3278 REAL(kind_phys), parameter :: ZDMIN = 5.e2_kind_phys
3279 real(kind_phys), parameter :: evapovtrm=evapr/vterm
3280
3281!NOTE
3282! downdraft area ffraction still needs to be computed for AW, assumed zero for now,
3283! as passed in.
3284
3285!
3286! Note: It is assumed that condensate evaporates in downdraft air.
3287!
3288 do k=1,kmax
3289 do i=ists,iens
3290 gprcp(i,k) = zero
3291 gsnwp(i,k) = zero
3292 gmdd(i,k) = zero
3293 gtevp(i,k) = zero
3294 evapd(i,k) = zero
3295 subld(i,k) = zero
3296 evape(i,k) = zero
3297 suble(i,k) = zero
3298 evapx(i,k) = zero
3299 sublx(i,k) = zero
3300 gmdde(i,k) = zero
3301 snmlt(i,k) = zero
3302 gchdd(i,k) = zero
3303 gcwdd(i,k) = zero
3304 gttev(i,k) = zero
3305 gtqev(i,k) = zero
3306 gcdsed(i,k) = zero
3307 gcqvd(i,k) = zero
3308! GCqlD (I,k) = zero
3309! GCqiD (I,k) = zero
3310 gtevap(i,k) = zero
3311 gtmelt(i,k) = zero
3312 gtsubl(i,k) = zero
3313 enddo
3314 enddo
3315
3316! These are zeroed by the calling routine, cs_cumlus
3317! do k=1,kmax
3318! do i=ists,iens
3319! dtdwn (I,k) = zero
3320! dqvdwn(I,k) = zero
3321! dqldwn(I,k) = zero
3322! dqidwn(I,k) = zero
3323! enddo
3324! enddo
3325! do n=ntrq,ntr
3326! do k=1,kmax
3327! do i=ists,iens
3328! dtrdwn(i,k,n) = zero
3329! enddo
3330! enddo
3331! enddo
3332!
3333 do i=ists,iens
3334 gchd(i) = zero
3335 gcwd(i) = zero
3336 gcud(i) = zero
3337 gcvd(i) = zero
3338 enddo
3339 do n=ntrq,ntr
3340 do i=ists,iens
3341 gctrd(i,n) = zero
3342 enddo
3343 enddo
3344!
3345 DO k=ktmx,1,-1 ! loop A
3346 kp1 = min(k+1,kmax)
3347!
3348! < precipitation melt & freeze >
3349!
3350 DO i=ists,iens
3351 if (kb(i) > 0) then
3352 gtprp = gprcp(i,kp1) + gsnwp(i,kp1)
3353 IF (gtprp > zero) THEN
3354 fsnow(i) = gsnwp(i,kp1) / gtprp
3355 ELSE
3356 fsnow(i) = zero
3357 ENDIF
3358 lvic = elocp + emeltocp*fsnow(i)
3359 gdtw = gdt(i,k) - lvic*(gdqs(i,k) - gdq(i,k,1)) &
3360 / (one + lvic*fdqs(i,k))
3361
3362 dz = gdzm(i,kp1) - gdzm(i,k)
3363 fmelt = (one + ftmlt*(gdtw - twsnow)) &
3364 * (one - tanh(gmflx(i,kp1)/gmflxc)) &
3365 * (one - tanh(vterms*meltau/dz))
3366 IF (gdtw < twsnow) THEN
3367 snmlt(i,k) = gprcp(i,kp1)*min(max(fmelt, one), zero)
3368 ELSE
3369 snmlt(i,k) = gsnwp(i,kp1)*max(min(fmelt, one), zero)
3370 ENDIF
3371 gsnwp(i,k) = gsnwp(i,kp1)+gsnwi(i,k) - snmlt(i,k)
3372 gprcp(i,k) = gprcp(i,kp1)+gprci(i,k) + snmlt(i,k)
3373 gttev(i,k) = -emeltocp * snmlt(i,k) * delpi(i,k)
3374!DD heating rate due to precip melting for AW
3375 gtmelt(i,k) = gtmelt(i,k) + gttev(i,k)
3376 endif
3377 ENDDO
3378!
3379! < downdraft >
3380!
3381 DO i=ists,iens ! loop B
3382 if (kb(i) > 0) then
3383 wrk = delpi(i,k)
3384 wrk1 = oneocp * wrk
3385 dz = gdzm(i,kp1) - gdzm(i,k)
3386 fevp = (one - tanh(evatau*vterm/dz))
3387 IF (gmdd(i,kp1) > zero) THEN
3388 gchx = gchd(i) / gmdd(i,kp1)
3389 gctx = gdt(i,k) + (gchx-gdhs(i,k)) / (cp+el*fdqs(i,k))
3390 gcqsx = gdqs(i,k) + fdqs(i,k) * (gctx - gdt(i,k))
3391 gcqsx = gcqsx*gmdd(i,kp1)
3392 evsu = max(gcqsx-gcwd(i), zero) * fevp
3393 gtprp = gprcp(i,k) + gsnwp(i,k)
3394 IF (gtprp > zero) THEN
3395 fsnow(i) = gsnwp(i,k) / gtprp
3396 ELSE
3397 fsnow(i) = zero
3398 ENDIF
3399 evapd(i,k) = min(evsu*(one-fsnow(i)), gprcp(i,k))
3400 subld(i,k) = min(evsu*fsnow(i), gsnwp(i,k))
3401 gprcp(i,k) = gprcp(i,k) - evapd(i,k)
3402 gsnwp(i,k) = gsnwp(i,k) - subld(i,k)
3403! temperature tendencies due to evaporation and sublimation of precip
3404! This is within downdraft
3405 gtevap(i,k) = gtevap(i,k) - elocp * evapd(i,k) * wrk
3406 gtsubl(i,k) = gtsubl(i,k) - esubocp * subld(i,k) * wrk
3407 gcwd(i) = gcwd(i) + evapd(i,k) + subld(i,k)
3408 gchd(i) = gchd(i) - emelt*subld(i,k)
3409 ENDIF
3410
3411 gmdd(i,k) = gmdd(i,kp1)
3412!
3413 lvic = elocp + emeltocp*fsnow(i)
3414 dqw = (gdqs(i,k) - gdw(i,k)) / (one + lvic*fdqs(i,k))
3415 dqw = max(dqw, zero)
3416 dtw = lvic*dqw
3417 gdqw = gdw(i,k) + dqw*fevp
3418!
3419 evsu = min(one, evapovtrm*dqw*dz*fevp)
3420 evape(i,k) = evsu*gprcp(i,k)
3421 suble(i,k) = evsu*gsnwp(i,k)
3422 gtevp(i,k) = evapd(i,k) + subld(i,k) + evape(i,k) + suble(i,k)
3423!
3424 gtprp = gprcp(i,k) + gsnwp(i,k)
3425 gprcp(i,k) = gprcp(i,k) - evape(i,k)
3426 gsnwp(i,k) = gsnwp(i,k) - suble(i,k)
3427! additional temperature tendencies due to evaporation and sublimation of precip
3428! This is outside of downdraft
3429 gtevap(i,k) = gtevap(i,k) - el*evape(i,k) * wrk1
3430 gtsubl(i,k) = gtsubl(i,k) - (el+emelt)*suble(i,k) * wrk1
3431!
3432 gmddd(i) = zero
3433 IF (gdz(i,k)-gdzm(i,1) > zdmin) THEN
3434 gteve = evape(i,k) + suble(i,k)
3435 gmddmx = revpdd*gteve/max(dqw, 1.d-10)
3436 gmdde(i,k) = rddr * (dtw*gtprp*delp(i,k))
3437 gmdde(i,k) = max(min(gmdde(i,k), gmddmx), zero)
3438 gmddx = gmdd(i,kp1) + gmdde(i,k)
3439 evsu = gmdde(i,k)*dqw*fevp
3440 IF (gteve > zero) THEN
3441 fsnow(i) = suble(i,k) / gteve
3442 ELSE
3443 fsnow(i) = zero
3444 END IF
3445 evapx(i,k) = (one-fsnow(i)) * evsu
3446 sublx(i,k) = fsnow(i) * evsu
3447!
3448 IF (gmddx > zero) THEN
3449 gdhi = gdh(i,k) - emelt*gdqi(i,k)
3450 gchdx = gchd(i) + gdhi*gmdde(i,k) - emelt*sublx(i,k)
3451 gcwdx = gcwd(i) + gdqw*gmdde(i,k)
3452 gcsd = (gchdx - el*gcwdx) / gmddx
3453 IF (gcsd < gds(i,k)) THEN
3454 gchd(i) = gchdx
3455 gcwd(i) = gcwdx
3456 gcud(i) = gcud(i) + gdu(i,k)*gmdde(i,k)
3457 gcvd(i) = gcvd(i) + gdv(i,k)*gmdde(i,k)
3458 do n = ntrq,ntr
3459 gctrd(i,n) = gctrd(i,n) + gdq(i,k,n)*gmdde(i,k)
3460 enddo
3461 gmdd(i,k) = gmddx
3462 evape(i,k) = evape(i,k) - evapx(i,k)
3463 suble(i,k) = suble(i,k) - sublx(i,k)
3464 evapd(i,k) = evapd(i,k) + evapx(i,k)
3465 subld(i,k) = subld(i,k) + sublx(i,k)
3466 gmddd(i) = zero
3467 ELSE
3468 gmdde(i,k) = zero
3469 gmddd(i) = gmdd(i,kp1)
3470 ENDIF
3471 ENDIF
3472 ELSE
3473 gmddd(i) = dz / (gdzm(i,kp1)-gdzm(i,1)) * gmdd(i,kp1)
3474 ENDIF
3475!
3476 gmddd(i) = max(gmddd(i), gmdd(i,k)-rddmx*gmflx(i,k))
3477!
3478 IF (gmddd(i) > zero) THEN
3479 fdet = gmddd(i)/gmdd(i,k)
3480 gchdd(i,k) = fdet*gchd(i)
3481 gcwdd(i,k) = fdet*gcwd(i)
3482 gcudd = fdet*gcud(i)
3483 gcvdd = fdet*gcvd(i)
3484 do n = ntrq,ntr
3485 gctrdd(n) = fdet*gctrd(i,n)
3486 enddo
3487!
3488 gthci = wrk * (gchdd(i,k) - gmddd(i)*gdh(i,k))
3489 gtqvci = wrk * (gcwdd(i,k) - gmddd(i)*gdq(i,k,1))
3490!
3491 gtt(i,k) = gtt(i,k) + (gthci - el*gtqvci)*oneocp
3492 gtq(i,k,1) = gtq(i,k,1) + gtqvci
3493 gtq(i,k,itl) = gtq(i,k,itl) - wrk * gmddd(i)*gdq(i,k,itl)
3494 gtq(i,k,iti) = gtq(i,k,iti) - wrk * gmddd(i)*gdqi(i,k)
3495
3496 do n = ntrq,ntr
3497 gtq(i,k,n) = gtq(i,k,n) + wrk * (gctrdd(n) - gmddd(i)*gdq(i,k,n))
3498 gctrd(i,n) = gctrd(i,n) - gctrdd(n)
3499 enddo
3500
3501 gtu(i,k) = gtu(i,k) + wrk * (gcudd - gmddd(i)*gdu(i,k))
3502 gtv(i,k) = gtv(i,k) + wrk * (gcvdd - gmddd(i)*gdv(i,k))
3503!
3504 gchd(i) = gchd(i) - gchdd(i,k)
3505 gcwd(i) = gcwd(i) - gcwdd(i,k)
3506 gcud(i) = gcud(i) - gcudd
3507 gcvd(i) = gcvd(i) - gcvdd
3508 gmdd(i,k) = gmdd(i,k) - gmddd(i)
3509 ENDIF
3510 gcdsed(i,k) = gchd(i) - el*gcwd(i)
3511 gcqvd(i,k) = gcwd(i)
3512 endif
3513 ENDDO ! loop B
3514!
3515 ENDDO ! loop A
3516!
3517 DO k=1,ktmx
3518 kp1 = min(k+1,kmax)
3519 DO i=ists,iens
3520 if (kb(i) > 0) then
3521 wrk = delpi(i,k)
3522 tx1 = delpi(i,kp1)
3523
3524 gttev(i,k) = gttev(i,k) - wrk &
3525 * (elocp*evape(i,k)+(elocp+emeltocp)*suble(i,k))
3526 gtt(i,k) = gtt(i,k) + gttev(i,k)
3527!
3528 gtqev(i,k) = gtqev(i,k) + (evape(i,k)+suble(i,k)) * wrk
3529 gtq(i,k,1) = gtq(i,k,1) + gtqev(i,k)
3530!
3531 gmflx(i,k) = gmflx(i,k) - gmdd(i,k)
3532 endif
3533 ENDDO ! end of i loop
3534 ENDDO ! end of k loop
3535
3536! AW tendencies due to vertical divergence of eddy fluxes
3537 DO k=2,ktmx
3538 kp1 = min(k+1,kmax)
3539 DO i=ists,iens
3540 if (kb(i) > 0) then
3541 if (k > 1 .and. flx_form) then
3542 fsigma = one - sigmad(i,kp1)
3543 dp_below = wrk * (one - sigmad(i,k))
3544 dp_above = tx1 * (one - sigmad(i,kp1))
3545
3546 wrk1 = gmdd(i,kp1) * (gdt(i,k)+gocp*gdz(i,k)) - gcdsed(i,kp1)*oneocp
3547 wrk2 = gmdd(i,kp1) * gdq(i,k,1) - gcqvd(i,kp1)
3548 wrk3 = gmdd(i,kp1) * gdq(i,k,itl)
3549 wrk4 = gmdd(i,kp1) * gdqi(i,k)
3550
3551 dtdwn(i,k) = dtdwn(i,k) + dp_below * wrk1
3552 dqvdwn(i,k) = dqvdwn(i,k) + dp_below * wrk2
3553 dqldwn(i,k) = dqldwn(i,k) + dp_below * wrk3 ! gcqld=0 - gcqld(i,k))
3554 dqidwn(i,k) = dqidwn(i,k) + dp_below * wrk4 ! gcqid=0 - gcqid(i,k))
3555
3556 dtdwn(i,kp1) = dtdwn(i,kp1) - dp_above * wrk1
3557 dqvdwn(i,kp1) = dqvdwn(i,kp1) - dp_above * wrk2
3558 dqldwn(i,kp1) = dqldwn(i,kp1) - dp_above * wrk3 ! gcqld=0 - gcqld(i,k))
3559 dqidwn(i,kp1) = dqidwn(i,kp1) - dp_above * wrk4 ! gcqid=0 - gcqid(i,k))
3560 do n = ntrq,ntr
3561 wrkn = gmdd(i,kp1) * gdq(i,k,n)
3562 dtrdwn(i,k,n) = dtrdwn(i,k,n) + dp_below * wrkn
3563 dtrdwn(i,kp1,n) = dtrdwn(i,kp1,n) - dp_above * wrkn
3564 enddo
3565 endif
3566
3567 endif
3568 ENDDO ! end of i loop
3569 ENDDO ! end of k loop
3570!
3571 END SUBROUTINE cumdwn
3572!***********************************************************************
3575 SUBROUTINE cumcld & !! cloudiness
3576 ( ijsdim, kmax , & !DD dimensions
3577 cumclw, qliq , qice , fliqc , & ! modified
3578 cumfrc, & ! output
3579 gmflx , ktmx , ists, iens ) ! input
3580!
3581 IMPLICIT NONE
3582
3583 INTEGER, INTENT(IN) :: IJSDIM, KMAX ! DD, for GFS, pass in
3584!
3585! [OUTPUT]
3586 REAL(kind_phys) CUMFRC(IJSDIM)
3587!
3588! [MODIFY]
3589 REAL(kind_phys) CUMCLW(IJSDIM, KMAX)
3590 REAL(kind_phys) QLIQ (IJSDIM, KMAX)
3591 REAL(kind_phys) QICE (IJSDIM, KMAX)
3592 REAL(kind_phys) FLIQC (IJSDIM, KMAX)
3593!
3594! [INPUT]
3595 REAL(kind_phys) GMFLX (IJSDIM, KMAX+1) ! cumulus mass flux
3596 INTEGER KTMX
3597 INTEGER ISTS, IENS
3598!
3599! [WORK]
3600 INTEGER I, K
3601 REAL(kind_phys) CUMF, QC, wrk
3602 LOGICAL, SAVE :: OFIRST = .true.
3603!
3604! [INTERNAL PARAM]
3605 REAL(kind_phys) :: FACLW = 0.1_kind_phys
3606 REAL(kind_phys) :: CMFMIN = 2.e-3_kind_phys
3607 REAL(kind_phys) :: CMFMAX = 3.e-1_kind_phys
3608 REAL(kind_phys) :: CLMIN = 1.e-3_kind_phys
3609 REAL(kind_phys) :: CLMAX = 0.1_kind_phys
3610 REAL(kind_phys), SAVE :: FACLF
3611!
3612 IF ( ofirst ) THEN
3613 faclf = (clmax-clmin)/log(cmfmax/cmfmin)
3614 ofirst = .false.
3615 END IF
3616
3617 cumfrc(ists:iens) = zero
3618 DO k=1,ktmx
3619 DO i=ists,iens
3620 cumfrc(i) = max(cumfrc(i), gmflx(i,k))
3621 ENDDO
3622 ENDDO
3623 DO i=ists,iens
3624 IF (cumfrc(i) > zero) THEN
3625 cumf = log(max(cumfrc(i), cmfmin)/cmfmin)
3626 cumfrc(i) = min(faclf*cumf+clmin, clmax)
3627 ENDIF
3628 ENDDO
3629!
3630 DO k=1,ktmx
3631 DO i=ists,iens
3632 IF (gmflx(i,k) > zero) THEN
3633 wrk = faclw / gmflx(i,k) * cumfrc(i)
3634 qliq(i,k) = wrk * qliq(i,k)
3635 qice(i,k) = wrk * qice(i,k)
3636 cumclw(i,k) = wrk * cumclw(i,k)
3637 qc = qliq(i,k) + qice(i,k)
3638 IF (qc > zero) THEN
3639 fliqc(i,k) = qliq(i,k) / qc
3640 ENDIF
3641 ENDIF
3642 ENDDO
3643 ENDDO
3644!
3645 END SUBROUTINE cumcld
3646!***********************************************************************
3649 SUBROUTINE cumupr & !! Tracer Updraft
3650 ( im , ijsdim, kmax , ntr , & !DD dimensions
3651 gtr , gprcc , & ! modified
3652 gdr , cbmfx , & ! input
3653 gcym , gcyt , gcqt , gclt , gcit , & ! input
3654 gtprt , gtevp , gtprc0, & ! input
3655 kb , kbmx , kt , ktmx , ktmxt , & ! input
3656 delpi , otspt , ists , iens, & ! input
3657 fscav, fswtr, nctp)
3658!
3659 IMPLICIT NONE
3660
3661 INTEGER, INTENT(IN) :: im, IJSDIM, KMAX, NTR, nctp !! DD, for GFS, pass in
3662!
3663! [MODIFY]
3664 REAL(kind_phys) GTR (IJSDIM, KMAX, NTR)
3665 REAL(kind_phys) GPRCC (IJSDIM, NTR)
3666!
3667! [INPUT]
3668 REAL(kind_phys) GDR (IJSDIM, KMAX, NTR)
3669 REAL(kind_phys) CBMFX (IM, NCTP)
3670 REAL(kind_phys) GCYM (IJSDIM, KMAX, nctp)
3671 REAL(kind_phys) GCYT (IJSDIM, NCTP)
3672 REAL(kind_phys) GCQT (IJSDIM, NCTP)
3673 REAL(kind_phys) GCLT (IJSDIM, NCTP)
3674 REAL(kind_phys) GCIT (IJSDIM, NCTP)
3675 REAL(kind_phys) GTPRT (IJSDIM, NCTP)
3676 REAL(kind_phys) GTEVP (IJSDIM, KMAX)
3677 REAL(kind_phys) GTPRC0(IJSDIM) !! precip. before evap.
3678 real(kind_phys) fscav(ntr), fswtr(ntr)
3679 INTEGER KB (IJSDIM )
3680 INTEGER KBMX
3681 INTEGER KT (IJSDIM, NCTP)
3682 INTEGER KTMX (NCTP)
3683 INTEGER KTMXT
3684 REAL(kind_phys) DELPI (IJSDIM, KMAX)
3685 LOGICAL OTSPT (NTR) !! transport with this routine?
3686 INTEGER ISTS, IENS
3687!
3688! [INTERNAL WORK]
3689 INTEGER I, K, LT, TP, CTP
3690 REAL(kind_phys) :: GCRTD, SCAV, GCWT, GPRCR, evpf, cbmfxl
3691 REAL(kind_phys), dimension(ists:iens) :: GCRB, GCRT, DR, gtprc0i
3692! REAL(kind_phys), dimension(ists:iens,kmax) :: DGCB, DZ, RDZM, EVPF
3693! REAL(kind_phys), dimension(ists:iens,nctp) :: DZT, RGCWT, MASK1, MASK2
3694 REAL(kind_phys), dimension(ists:iens,nctp) :: RGCWT, MASK1
3695!
3696 do i=ists,iens
3697 if (gtprc0(i) > zero) then
3698 gtprc0i(i) = one / gtprc0(i)
3699 else
3700 gtprc0i(i) = zero
3701 endif
3702 enddo
3703 DO ctp=1,nctp
3704 DO i=ists,iens
3705 k = kt(i,ctp)
3706!
3707 gcwt = gcqt(i,ctp) + gclt(i,ctp) + gcit(i,ctp)
3708 rgcwt(i,ctp) = zero
3709 IF (gcwt > zero) THEN
3710 rgcwt(i,ctp) = one / gcwt
3711 ENDIF
3712!
3713 mask1(i,ctp) = zero
3714 IF (kb(i) > 0 .and. k > kb(i)) THEN
3715 mask1(i,ctp) = one
3716 ENDIF
3717! MASK2(I,CTP) = zero
3718! IF (CBMFX(I,CTP) > zero) then
3719! MASK2(I,CTP) = one
3720! ENDIF
3721 ENDDO
3722 ENDDO
3723!
3724 DO lt=1,ntr ! outermost tracer LT loop
3725!
3726 IF (otspt(lt)) THEN
3727 DO ctp=1,nctp
3728 DO i=ists,iens
3729 gcrb(i) = zero
3730 dr(i) = zero
3731 enddo
3732 DO k=1,kbmx
3733 DO i=ists,iens
3734 IF (kb(i) > 0 .and. k < kb(i)) THEN
3735 gcrb(i) = gcrb(i) + (gcym(i,k+1,ctp)-gcym(i,k,ctp))* gdr(i,k,lt)
3736 ENDIF
3737 ENDDO
3738 ENDDO
3739!
3740 DO k=2,ktmx(ctp)
3741 DO i=ists,iens
3742 IF (kb(i) > 0 .and. k >= kb(i) .AND. k < kt(i,ctp)) THEN
3743 dr(i) = dr(i) + (gcym(i,k+1,ctp)-gcym(i,k,ctp)) * gdr(i,k,lt)
3744 ENDIF
3745 ENDDO
3746 ENDDO
3747!
3748 DO i=ists,iens
3749 k = kt(i,ctp)
3750 if (kb(i) > 0 .and. k > kb(i)) then
3751 dr(i) = dr(i) + (gcyt(i,ctp) - gcym(i,k,ctp)) * gdr(i,k,lt) &
3752 * mask1(i,ctp)
3753 gcrt(i) = (gcrb(i) + dr(i)) * mask1(i,ctp)
3754!
3755 scav = fscav(lt)*gtprt(i,ctp) + fswtr(lt)*gtprt(i,ctp)*rgcwt(i,ctp)
3756 scav = min(scav, one)
3757 gcrtd = gcrt(i) * (one - scav)
3758 cbmfxl = max(zero, cbmfx(i,ctp))
3759 gprcr = scav * gcrt(i) * cbmfxl
3760
3761 gtr(i,k,lt) = gtr(i,k,lt) + delpi(i,k) * cbmfxl &
3762 * (gcrtd - gcyt(i,ctp) * gdr(i,k,lt))
3763 gprcc(i,lt) = gprcc(i,lt) + gprcr
3764
3765! GPRCR = SCAV * GCRT(I) * CBMFX(I,CTP)
3766! GTR(I,K,LT) = GTR(I,K,LT) + DELPI(I,K) * CBMFX(I,CTP) &
3767! * (GCRTD - GCYT(I,CTP) * GDR(I,K,LT)) * MASK2(I,CTP)
3768! GPRCC(I,LT) = GPRCC(I,LT) + GPRCR * MASK2(I,CTP)
3769 endif
3770 ENDDO
3771 ENDDO
3772!
3773 DO k=ktmxt,1,-1
3774 DO i=ists,iens
3775 evpf = gtevp(i,k) * gtprc0i(i)
3776 gtr(i,k,lt) = gtr(i,k,lt) + delpi(i,k) * gprcc(i,lt) * evpf
3777 gprcc(i,lt) = gprcc(i,lt) * (one - evpf)
3778! GTR(I,K,LT) = GTR(I,K,LT) + DELPI(I,K) * GPRCC(I,LT) * EVPF(I,K)
3779! GPRCC(I,LT) = GPRCC(I,LT) * (one - EVPF(I,K))
3780 ENDDO
3781 ENDDO
3782!
3783 ENDIF
3784!
3785 ENDDO ! outermost tracer LT loop
3786!
3787 END SUBROUTINE cumupr
3788!***********************************************************************
3790 SUBROUTINE cumdnr & !! Tracer Downdraft
3791 ( im , ijsdim, kmax , ntr , & !DD dimensions
3792 gtr , & ! modified
3793 gdr , gmdd , delpi , & ! input
3794 ktmx , otspt , ists , iens ) ! input
3795!
3796 IMPLICIT NONE
3797
3798 INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR !! DD, for GFS, pass in
3799!
3800! [MODIFY]
3801 REAL(kind_phys) GTR (IJSDIM, KMAX, NTR) ! Temperature tendency
3802!
3803! [INPUT]
3804 REAL(kind_phys) GDR (IJSDIM, KMAX, NTR)
3805 REAL(kind_phys) GMDD (IJSDIM, KMAX) ! downdraft mass flux
3806 REAL(kind_phys) DELPI (IJSDIM, KMAX )
3807 LOGICAL OTSPT (NTR)
3808 INTEGER KTMX, ISTS, IENS
3809!
3810! [INTERNAL WORK]
3811 REAL(kind_phys) GCRD (ISTS:IENS) ! downdraft q
3812 REAL(kind_phys) GMDDE, GMDDD, GCRDD
3813 INTEGER I, K, LT, kp1
3814!
3815!
3816 DO lt=1,ntr
3817 IF (otspt(lt)) THEN
3818 gcrd = zero
3819 DO k=ktmx,1,-1
3820 kp1 = min(k+1,kmax)
3821 DO i=ists,iens
3822 gmdde = gmdd(i,k) - gmdd(i,kp1)
3823 IF (gmdde >= zero) THEN
3824 gcrd(i) = gcrd(i) + gdr(i,k,lt)*gmdde
3825 ELSEIF (gmdd(i,kp1) > zero) THEN
3826 gmddd = - gmdde
3827 gcrdd = gmddd/gmdd(i,kp1) * gcrd(i)
3828 gtr(i,k,lt) = gtr(i,k,lt) + delpi(i,k) &
3829 * (gcrdd - gmddd*gdr(i,k,lt))
3830 gcrd(i) = gcrd(i) - gcrdd
3831 ENDIF
3832 ENDDO
3833 ENDDO
3834 ENDIF
3835 ENDDO
3836!
3837 END SUBROUTINE cumdnr
3838!***********************************************************************
3840 SUBROUTINE cumsbr & !! Tracer Subsidence
3841 ( im , ijsdim, kmax, ntr, nctp, & !DD dimensions
3842 gtr , & ! modified
3843 gdr , delp , & ! input
3844 gmflx , ktmx , otspt , & ! input
3845 sigmai , sigma , & !DDsigma input
3846 ists, iens ) ! input
3847!
3848 IMPLICIT NONE
3849
3850 INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR, nctp !! DD, for GFS, pass in
3851!
3852! [MODIFY]
3853 REAL(kind_phys) GTR (IJSDIM, KMAX, NTR) !! tracer tendency
3854!
3855! [INPUT]
3856 REAL(kind_phys) GDR (IJSDIM, KMAX, NTR) !! tracer
3857 REAL(kind_phys) DELP (IJSDIM, KMAX)
3858 REAL(kind_phys) GMFLX (IJSDIM, KMAX+1) !! mass flux
3859 INTEGER KTMX
3860 LOGICAL OTSPT (NTR) !! tracer transport on/off
3861 INTEGER ISTS, IENS
3862 REAL(kind_phys) sigmai (IM,KMAX+1,NCTP), sigma(IM,KMAX+1) !!DDsigma cloud updraft fraction
3863!
3864! [INTERNAL WORK]
3865 INTEGER I, K, KM, KP, LT
3866 REAL(kind_phys) SBR0, SBR1, FX1
3867 REAL(kind_phys) FX(ISTS:IENS)
3868!
3869 DO lt=1,ntr
3870 IF (otspt(lt)) THEN
3871 DO i=ists,iens
3872 fx(i) = zero
3873 enddo
3874 DO k=ktmx,1,-1
3875 km = max(k-1, 1)
3876 kp = min(k+1, kmax)
3877 DO i=ists,iens
3878 sbr0 = gmflx(i,k+1) * (gdr(i,kp,lt) - gdr(i,k,lt))
3879 sbr1 = gmflx(i,k) * (gdr(i,k,lt) - gdr(i,km,lt))
3880 IF (gmflx(i,k) > gmflx(i,k+1)) THEN
3881 fx1 = half
3882 ELSE
3883 fx1 = zero
3884 END IF
3885 gtr(i,k,lt) = gtr(i,k,lt) + grav/delp(i,k) &
3886 * ((one-fx(i))*sbr0 + fx1*sbr1)
3887 fx(i) = fx1
3888 ENDDO
3889 ENDDO
3890 ENDIF
3891 ENDDO
3892!
3893 END SUBROUTINE cumsbr
3894!*********************************************************************
3897 SUBROUTINE cumfxr & ! Tracer mass fixer
3898 ( im , ijsdim, kmax , ntr , & !DD dimensions
3899 gtr , & ! modified
3900 gdr , delp , delta , ktmx , imfxr , & ! input
3901 ists , iens ) ! input
3902!
3903 IMPLICIT NONE
3904
3905 INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR !! DD, for GFS, pass in
3906!
3907! [MODIFY]
3908 REAL(kind_phys) GTR (IJSDIM, KMAX, NTR) ! tracer tendency
3909!
3910! [INPUT]
3911 REAL(kind_phys) GDR (IJSDIM, KMAX, NTR) ! tracer
3912 REAL(kind_phys) DELP (IJSDIM, KMAX)
3913 REAL(kind_phys) DELTA ! time step
3914 INTEGER KTMX
3915 INTEGER IMFXR (NTR)
3916 ! 0: mass fixer is not applied
3917 ! tracers which may become negative values
3918 ! e.g. subgrid-PDFs
3919 ! 1: mass fixer is applied, total mass may change through cumulus scheme
3920 ! e.g. moisture, liquid cloud, ice cloud, aerosols
3921 ! 2: mass fixer is applied, total mass never change through cumulus scheme
3922 ! e.g. CO2
3923 !DD add new CASE
3924 ! 3: just fill holes, no attempt to conserve
3925 INTEGER ISTS, IENS
3926!
3927! [INTERNAL WORK]
3928 REAL(kind_phys) GDR1
3929 REAL(kind_phys) GDR2 (ISTS:IENS, KMAX)
3930 REAL(kind_phys), dimension(ISTS:IENS) :: TOT0, TOT1, TRAT
3931 REAL(kind_phys) FWAT
3932 INTEGER I, K, LT
3933!
3934! Attention: tracers are forced to be positive unless IMFXR=0.
3935!
3936 DO lt=1,ntr
3937 SELECT CASE (imfxr(lt))
3938 CASE (0)
3939 cycle
3940 CASE (1)
3941 fwat = one
3942 CASE (2)
3943 fwat = zero
3944 CASE (3)
3945 CASE DEFAULT
3946 EXIT
3947 END SELECT
3948!
3949 DO i=ists,iens
3950 tot0(i) = zero
3951 tot1(i) = zero
3952 trat(i) = one
3953 enddo
3954!
3955 DO k=ktmx,1,-1
3956 DO i=ists,iens
3957 IF (gtr(i,k,lt) /= zero) THEN
3958 gdr1 = gdr(i,k,lt) + delta*gtr(i,k,lt)
3959 gdr2(i,k) = max(gdr1, zero)
3960 gdr1 = gdr1 * fwat + gdr(i,k,lt)*(one - fwat)
3961 tot0(i) = tot0(i) + gdr1 *(delp(i,k)*gravi)
3962 tot1(i) = tot1(i) + gdr2(i,k)*(delp(i,k)*gravi)
3963 ENDIF
3964 ENDDO
3965 ENDDO
3966!
3967 if(imfxr(lt) .ne. 3) then
3968 DO i=ists,iens
3969 IF (tot1(i) > zero ) THEN
3970 trat(i) = max(tot0(i), zero) / tot1(i)
3971 ENDIF
3972 ENDDO
3973 endif
3974!
3975 DO k=ktmx,1,-1
3976 DO i=ists,iens
3977 IF (gtr(i,k,lt) /= zero ) THEN
3978 gdr2(i,k ) = gdr2(i,k)*trat(i)
3979 gtr(i,k,lt) = (gdr2(i,k)-gdr(i,k,lt)) / delta
3980 ENDIF
3981 ENDDO
3982 ENDDO
3983!
3984 ENDDO ! LT-loop
3985!
3986 END SUBROUTINE cumfxr
3987!*********************************************************************
3989 SUBROUTINE cumfxr1 & ! Tracer mass fixer
3990 ( im , ijsdim, kmax ,nctp, & !DD dimensions
3991 gtr , & ! modified
3992 gdr , delp , delta , ktmx , imfxr , & ! input
3993 ists , iens ) ! input
3994!
3995 IMPLICIT NONE
3996
3997 INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, nctp !! DD, for GFS, pass in
3998!
3999! [MODIFY]
4000 REAL(kind_phys) GTR (IJSDIM, KMAX) ! tracer tendency
4001!
4002! [INPUT]
4003 REAL(kind_phys) GDR (IJSDIM, KMAX) ! tracer
4004 REAL(kind_phys) DELP (IJSDIM, KMAX)
4005 REAL(kind_phys) DELTA ! time step
4006 INTEGER KTMX
4007 INTEGER IMFXR
4008 ! 0: mass fixer is not applied
4009 ! tracers which may become negative values
4010 ! e.g. subgrid-PDFs
4011 ! 1: mass fixer is applied, total mass may change through cumulus scheme
4012 ! e.g. moisture, liquid cloud, ice cloud, aerosols
4013 ! 2: mass fixer is applied, total mass never change through cumulus scheme
4014 ! e.g. CO2
4015 INTEGER ISTS, IENS
4016!
4017! [INTERNAL WORK]
4018 REAL(kind_phys) GDR1
4019 REAL(kind_phys) GDR2 (ISTS:IENS, KMAX)
4020 REAL(kind_phys), dimension(ISTS:IENS) :: TOT0, TOT1, TRAT
4021 REAL(kind_phys) FWAT
4022 INTEGER I, K
4023!
4024! Attention: tracers are forced to be positive unless IMFXR=0.
4025!
4026 SELECT CASE (imfxr)
4027 CASE (0)
4028 RETURN
4029 CASE (1)
4030 fwat = one
4031 CASE (2)
4032 fwat = zero
4033 CASE DEFAULT
4034 RETURN
4035 END SELECT
4036!
4037 DO i=ists,iens
4038 tot0(i) = zero
4039 tot1(i) = zero
4040 enddo
4041!
4042 DO k=ktmx,1,-1
4043 DO i=ists,iens
4044 IF (gtr(i,k) /= zero) THEN
4045 gdr1 = gdr(i,k) + delta*gtr(i,k)
4046 gdr2(i,k) = max(gdr1, zero)
4047 gdr1 = gdr1*fwat + gdr(i,k)*(one - fwat)
4048 tot0(i) = tot0(i) + gdr1 *(delp(i,k)*gravi)
4049 tot1(i) = tot1(i) + gdr2(i,k)*(delp(i,k)*gravi)
4050 ENDIF
4051 ENDDO
4052 ENDDO
4053!
4054 DO i=ists,iens
4055 IF (tot1(i) > zero) THEN
4056 trat(i) = max(tot0(i), zero) / tot1(i)
4057 ELSE
4058 trat(i) = one
4059 ENDIF
4060 ENDDO
4061!
4062 DO k=ktmx,1,-1
4063 DO i=ists,iens
4064 IF (gtr(i,k) /= zero) THEN
4065 gdr2(i,k) = gdr2(i,k)*trat(i)
4066 gtr(i,k) = (gdr2(i,k)-gdr(i,k)) / delta
4067 ENDIF
4068 ENDDO
4069 ENDDO
4070!
4071 END SUBROUTINE cumfxr1
4072!*********************************************************************
4074 SUBROUTINE cumchk & ! check range of output values
4075 ( ijsdim, kmax , ntr , & !DD dimensions
4076 gtt , gtq , gtu , gtv , & ! input
4077 gtcfrc, gprcc , gsnwc , cumclw, & ! input
4078 cumfrc, fliqc , gtprp , & ! input
4079 ists , iens ) ! input
4080!
4081 IMPLICIT NONE
4082
4083 INTEGER, INTENT(IN) :: IJSDIM, KMAX, NTR ! DD, for GFS, pass in
4084!
4085! [INPUT]
4086 REAL(kind_phys) GTT (IJSDIM, KMAX) ! heating rate
4087 REAL(kind_phys) GTQ (IJSDIM, KMAX, NTR) ! change in q
4088 REAL(kind_phys) GTU (IJSDIM, KMAX) ! tendency of u
4089 REAL(kind_phys) GTV (IJSDIM, KMAX) ! tendency of v
4090 REAL(kind_phys) GPRCC (IJSDIM, NTR ) ! rainfall
4091 REAL(kind_phys) GSNWC (IJSDIM) ! snowfall
4092 REAL(kind_phys) CUMCLW(IJSDIM, KMAX) ! cloud water in cumulus
4093 REAL(kind_phys) CUMFRC(IJSDIM) ! cumulus cloud fraction
4094 REAL(kind_phys) GTCFRC(IJSDIM, KMAX) ! change in cloud fraction
4095 REAL(kind_phys) FLIQC (IJSDIM, KMAX) ! liquid ratio in cumulus
4096 REAL(kind_phys) GTPRP (IJSDIM, KMAX) ! rain+snow flux
4097!
4098 INTEGER ISTS, IENS
4099!
4100! [INTERNAL WORK]
4101 INTEGER I, K
4102!
4103! [INTERNAL PARM]
4104 REAL(kind_phys) :: GTTMAX = 1.e-2_kind_phys
4105 REAL(kind_phys) :: GTQVMAX = 1.e-4_kind_phys
4106 REAL(kind_phys) :: GTQLMAX = 1.e-5_kind_phys
4107 REAL(kind_phys) :: GTUMAX = 1.e-2_kind_phys
4108 REAL(kind_phys) :: GTVMAX = 1.e-2_kind_phys
4109 REAL(kind_phys) :: GTCFMAX = 1.e-3_kind_phys
4110 REAL(kind_phys) :: PRCCMAX = 1.e-2_kind_phys
4111 REAL(kind_phys) :: SNWCMAX = 1.e-2_kind_phys
4112 REAL(kind_phys) :: CLWMAX = 1.e-3_kind_phys
4113 REAL(kind_phys) :: TPRPMAX = 1.e-2_kind_phys
4114 REAL(kind_phys) :: GTQIMAX = 1.e-5_kind_phys
4115 !REAL(kind_phys) :: GTM2MAX = 1._kind_phys
4116 !REAL(kind_phys) :: GTM3MAX = 1._kind_phys
4117!
4118 DO k=1,kmax
4119 DO i=ists, iens
4120 IF (abs(gtt(i,k)) > gttmax) THEN
4121 WRITE(iulog,*) '### CUMCHK: GTT(',i,',',k,')=',gtt(i,k)
4122 ENDIF
4123 IF (abs(gtq(i,k,1) ) > gtqvmax) THEN
4124 WRITE(iulog,*) '### CUMCHK: GTQ(',i,',',k,',1 )=', gtq(i,k,1)
4125 ENDIF
4126 IF (abs(gtq(i,k,itl)) > gtqlmax) THEN
4127 WRITE(iulog,*) '### CUMCHK: GTQ(',i,',',k,',ITL )=', gtq(i,k,itl)
4128 ENDIF
4129 IF (abs(gtu(i,k)) > gtumax) THEN
4130 WRITE(iulog,*) '### CUMCHK: GTU(',i,',',k,')=',gtu(i,k)
4131 END IF
4132 IF (abs(gtv(i,k)) > gtvmax) THEN
4133 WRITE(iulog,*) '### CUMCHK: GTV(',i,',',k,')=',gtv(i,k)
4134 ENDIF
4135 IF (abs(gtcfrc(i,k)) > gtcfmax) THEN
4136 WRITE(iulog,*) '### CUMCHK: GTCFRC(',i,',',k,')=', gtcfrc(i,k)
4137 ENDIF
4138 IF (cumclw(i,k) > clwmax .OR. cumclw(i,k) < zero) THEN
4139 WRITE(iulog,*) '### CUMCHK: CUMCLW(',i,',',k,')=', cumclw(i,k)
4140 ENDIF
4141 IF (fliqc(i,k) > one .OR. fliqc(i,k) < zero) THEN
4142 WRITE(iulog,*) '### CUMCHK: FLIQC(',i,',',k,')=', fliqc(i,k)
4143 ENDIF
4144 IF (gtprp(i,k) > tprpmax .OR. gtprp(i,k) < zero) THEN
4145 WRITE(iulog,*) '### CUMCHK: GTPRP(',i,',',k,')=', gtprp(i,k)
4146 ENDIF
4147 IF (abs(gtq(i,k,iti)) > gtqimax) THEN
4148 WRITE(iulog,*) '### CUMCHK: GTQ(',i,',',k,',ITI )=', gtq(i,k,iti)
4149 ENDIF
4150! IF (ABS(GTQ(I,K,IMU2) ) > GTM2MAX) THEN
4151! WRITE(iulog,*) '### CUMCHK: GTQ(',I,',',K,',IMU2 )=', GTQ(I,K,IMU2)
4152! ENDIF
4153! IF (ABS(GTQ(I,K,IMU3)) > GTM3MAX) THEN
4154! WRITE(iulog,*) '### CUMCHK: GTQ(',I,',',K,',IMU3 )=', GTQ(I,K,IMU3)
4155! ENDIF
4156 ENDDO
4157 ENDDO
4158!
4159 DO i=ists,iens
4160 IF (gprcc(i,1) > prccmax .OR. gprcc(i,1) < zero) THEN
4161 WRITE(iulog,*) '### CUMCHK: GPRCC(',i,')=',gprcc(i,1)
4162 END IF
4163 IF (gsnwc(i) > snwcmax .OR. gsnwc(i) < zero) THEN
4164 WRITE(iulog,*) '### CUMCHK: GSNWC(',i,')=',gsnwc(i)
4165 END IF
4166 IF (cumfrc(i) > one .OR. cumfrc(i) < zero) THEN
4167 WRITE(iulog,*) '### CUMCHK: CUMFRC(',i,')=',cumfrc(i)
4168 ENDIF
4169 ENDDO
4170!
4171 END SUBROUTINE cumchk
4172
4173!***********************************************************************
4174
4175end module cs_conv
subroutine esat(t, esw, esi, desw, desi)
use polynomials to calculate saturation vapor pressure and derivative with respect to temperature: ov...
subroutine cumsbr(im, ijsdim, kmax, ntr, nctp, gtr, gdr, delp, gmflx, ktmx, otspt, sigmai, sigma, ists, iens)
Definition cs_conv.F90:3847
subroutine, public cs_conv_run(ijsdim, kmax, ntracp1, nn, ntr, nctp, otspt, lat, kdt, t, q, rain1, clw, zm, zi, pap, paph, delta, delti, ud_mf, dd_mf, dt_mf, u, v, fscav, fswtr, cbmfx, mype, wcbmaxm, precz0in, preczhin, clmdin, sigma, do_aw, do_awdd, flx_form, lprnt, ipr, kcnv, qlcn, qicn, w_upi, cf_upi, cnv_mfd, cnv_dqldt, clcn, cnv_fice, cnv_ndrop, cnv_nice, mp_phys, errmsg, errflg)
Definition cs_conv.F90:167
subroutine cumchk(ijsdim, kmax, ntr, gtt, gtq, gtu, gtv, gtcfrc, gprcc, gsnwc, cumclw, cumfrc, fliqc, gtprp, ists, iens)
Definition cs_conv.F90:4080
subroutine cumsbw(im, ijsdim, kmax, gtu, gtv, gdu, gdv, delpi, gmflx, gmfx0, ktmx, cpres, kb, ists, iens)
This subroutine calculate cloud subsidence heating.
Definition cs_conv.F90:3076
subroutine cumup(ijsdim, kmax, ntr, ntrq, acwf, gclz, gciz, gprciz, gsnwiz, gcyt, gcht, gcqt, gclt, gcit, gtprt, gcut, gcvt, gctrt, kt, ktmx, gcym, wcv, gchb, gcwb, gcub, gcvb, gcib, gctrb, gdu, gdv, gdh, gdw, gdhs, gdqs, gdt, gdtm, gdq, gdqi, gdz, gdzm, gdpm, fdqs, gam, gdztr, cpres, wcb, kb, ctp, ists, iens, gctm, gcqm, gcwm, gchm, gcwt, gclm, gcim, gctrm, lprnt, ipr)
This subroutine calculates in-cloud properties.
Definition cs_conv.F90:2075
subroutine cumcld(ijsdim, kmax, cumclw, qliq, qice, fliqc, cumfrc, gmflx, ktmx, ists, iens)
This subroutine computes cumulus cloudiness.
Definition cs_conv.F90:3580
subroutine cumsbh(im, ijsdim, kmax, ntr, ntrq, gtt, gtq, gtu, gtv, gdh, gdq, gdqi, gdu, gdv, delpi, gmflx, gmfx0, ktmx, cpres, kb, ists, iens)
Definition cs_conv.F90:2954
subroutine cumbmx(ijsdim, kmax, cbmfx, acwf, gcyt, gdzm, gdw, gdqs, delp, kt, ktmx, kb, delt, ists, iens)
This subroutine computes cloud base mass flux.
Definition cs_conv.F90:2732
subroutine cumfxr(im, ijsdim, kmax, ntr, gtr, gdr, delp, delta, ktmx, imfxr, ists, iens)
This subroutine calculates tracer mass fixer without detrainment.
Definition cs_conv.F90:3902
subroutine cumdnr(im, ijsdim, kmax, ntr, gtr, gdr, gmdd, delpi, ktmx, otspt, ists, iens)
Definition cs_conv.F90:3795
subroutine cumfxr1(im, ijsdim, kmax, nctp, gtr, gdr, delp, delta, ktmx, imfxr, ists, iens)
Definition cs_conv.F90:3994
subroutine cumdwn(im, ijsdim, kmax, ntr, ntrq, nctp, gtt, gtq, gtu, gtv, gmflx, gprcp, gsnwp, gtevp, gmdd, gprci, gsnwi, gdh, gdw, gdq, gdqi, gdqs, gds, gdhs, gdt, gdu, gdv, gdz, gdzm, fdqs, delp, delpi, sigmad, do_aw, do_awdd, flx_form, gtmelt, gtevap, gtsubl, dtdwn, dqvdwn, dqldwn, dqidwn, dtrdwn, kb, ktmx, ists, iens)
This subroution calculates freeze, melt and evaporation in cumulus downdraft.
Definition cs_conv.F90:3157
subroutine cumupr(im, ijsdim, kmax, ntr, gtr, gprcc, gdr, cbmfx, gcym, gcyt, gcqt, gclt, gcit, gtprt, gtevp, gtprc0, kb, kbmx, kt, ktmx, ktmxt, delpi, otspt, ists, iens, fscav, fswtr, nctp)
This subroutine calculates.
Definition cs_conv.F90:3658
subroutine cumflx(im, ijsdim, kmax, gmflx, gprci, gsnwi, cmdet, qliq, qice, gtprc0, cbmfx, gcym, gprciz, gsnwiz, gtprt, gclz, gciz, gcyt, kb, kt, ktmx, ists, iens)
This subroutine computes cloud mass flux & precip.
Definition cs_conv.F90:2818
subroutine cumbas(ijsdim, kmax, kb, gcym, kbmx, ntr, ntrq, gchb, gcwb, gcub, gcvb, gcib, gctrb, gdh, gdw, gdhs, gdqs, gdqi, gdu, gdv, gdzm, gdpm, fdqs, gam, lprnt, ipr, ists, iens, gctbl, gcqbl, gdq, gcwbl, gcqlbl, gcqibl, gctrbl)
This subroutine calculates cloud base properties.
Definition cs_conv.F90:1857
subroutine cumdet(im, ijsdim, kmax, ntr, ntrq, gtt, gtq, gtu, gtv, gdh, gdq, gdu, gdv, cbmfx, gcyt, delpi, gcht, gcqt, gclt, gcit, gcut, gcvt, gdqi, gctrt, kt, ists, iens, nctp)
This subroutine calculates cloud detrainment heating.
Definition cs_conv.F90:2880
subroutine cs_cumlus(im, ijsdim, kmax, ntr, otspt1, otspt2, lprnt, ipr, gtt, gtq, gtu, gtv, cmdet, gtprp, gsnwp, gmfx0, gmfx1, cape, kt, cbmfx, gdt, gdq, gdu, gdv, gdtm, gdp, gdpm, gdz, gdzm, delp, delpinv, delta, delti, ists, iens, mype, fscav, fswtr, wcbmaxm, nctp, sigmai, sigma, vverti, do_aw, do_awdd, flx_form)
Main subroutine for the cumulus parameterization with state-dependent entrainment rate developed by M...
Definition cs_conv.F90:603
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26
This module contains some of the most frequently used math and physics constants for GCM models.
Definition physcons.F90:36