CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
cu_ntiedtke.F90
1
7!###########################################################
8
10
11!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12 use machine , only : kind_phys
13 ! DH* TODO - replace with arguments to subroutine calls,
14 ! this also requires redefining derived constants in the
15 ! parameter section below
16 use physcons, only:rd=>con_rd, rv=>con_rv, g=>con_g, &
17 & cpd=>con_cp, alv=>con_hvap, alf=>con_hfus
18
19 implicit none
20 real(kind=kind_phys),private :: rcpd,vtmpc1,tmelt,als,t13, &
21 c1es,c2es,c3les,c3ies,c4les,c4ies,c5les,c5ies,zrg
22
23 real(kind=kind_phys),private :: rovcp,r5alvcp,r5alscp,ralvdcp,ralsdcp,ralfdcp,rtwat,rtber,rtice
24 real(kind=kind_phys),private :: entrdd,cmfcmax,cmfcmin,cmfdeps,zdnoprc,cprcon
25 integer,private :: momtrans,p650
26
27 parameter( &
28 t13 = 0.333333333,&
29 rcpd=1.0/cpd, &
30 tmelt=273.16, &
31 zrg=1.0/g, &
32 c1es=610.78, &
33 c2es=c1es*rd/rv, &
34 c3les=17.2693882, &
35 c3ies=21.875, &
36 c4les=35.86, &
37 c4ies=7.66, &
38 als = alv+alf, &
39 c5les=c3les*(tmelt-c4les), &
40 c5ies=c3ies*(tmelt-c4ies), &
41 r5alvcp=c5les*alv*rcpd, &
42 r5alscp=c5ies*als*rcpd, &
43 ralvdcp=alv*rcpd, &
44 ralsdcp=als*rcpd, &
45 ralfdcp=alf*rcpd, &
46 rtwat=tmelt, &
47 rtber=tmelt-5., &
48 rtice=tmelt-23., &
49 vtmpc1=rv/rd-1.0, &
50 rovcp = rd*rcpd )
51!
52! entrdd: average entrainment & detrainment rate for downdrafts
53! ------
54!
55 parameter(entrdd = 2.0e-4)
56!
57! cmfcmax: maximum massflux value allowed for updrafts etc
58! -------
59!
60 parameter(cmfcmax = 1.0)
61!
62! cmfcmin: minimum massflux value (for safety)
63! -------
64!
65 parameter(cmfcmin = 1.e-10)
66!
67! cmfdeps: fractional massflux for downdrafts at lfs
68! -------
69!
70 parameter(cmfdeps = 0.30)
71
72! zdnoprc: deep cloud is thicker than this height (Unit:Pa)
73!
74 parameter(zdnoprc = 2.0e4)
75! -------
76!
77! cprcon: coefficient from cloud water to rain water
78!
79 parameter(cprcon = 1.4e-3)
80! -------
81!
82! momtrans: momentum transport method
83! ( 1 = IFS40r1 method; 2 = new method )
84!
85 parameter(momtrans = 2 )
86! -------
87!
88 logical :: isequil
89! isequil: representing equilibrium and nonequilibrium convection
90! ( .false. [default]; .true. [experimental]. Ref. Bechtold et al. 2014 JAS )
91!
92 parameter(isequil = .false. )
93!
94!--------------------
95! switches for deep, mid, shallow convections, downdraft, and momemtum transport
96! ------------------
97 logical :: lmfpen,lmfmid,lmfscv,lmfdd,lmfdudv
98 parameter(lmfpen=.true.,lmfmid=.true.,lmfscv=.true.,lmfdd=.true.,lmfdudv=.true.)
99!--------------------
100!#################### end of variables definition##########################
101!-----------------------------------------------------------------------
102!
103contains
109 subroutine cu_ntiedtke_init(imfshalcnv, imfshalcnv_ntiedtke, imfdeepcnv, &
110 imfdeepcnv_ntiedtke,mpirank, mpiroot, errmsg, errflg)
111
112 implicit none
113
114 integer, intent(in) :: imfshalcnv, imfshalcnv_ntiedtke
115 integer, intent(in) :: imfdeepcnv, imfdeepcnv_ntiedtke
116 integer, intent(in) :: mpirank
117 integer, intent(in) :: mpiroot
118 character(len=*), intent( out) :: errmsg
119 integer, intent( out) :: errflg
120
121 ! initialize ccpp error handling variables
122 errmsg = ''
123 errflg = 0
124
125 ! DH* temporary
126 if (mpirank==mpiroot) then
127 write(0,*) ' -----------------------------------------------------------------------------------------------------------------------------'
128 write(0,*) ' --- WARNING --- the CCPP New Tiedtke convection scheme is currently under development, use at your own risk --- WARNING ---'
129 write(0,*) ' -----------------------------------------------------------------------------------------------------------------------------'
130 end if
131 ! *DH temporary
132
133 ! Consistency checks
134 if (imfshalcnv/=imfshalcnv_ntiedtke) then
135 write(errmsg,'(*(a))') 'Logic error: namelist choice of', &
136 & ' shallow convection is different from new Tiedtke scheme'
137 errflg = 1
138 return
139 end if
140
141 if (imfdeepcnv/=imfdeepcnv_ntiedtke) then
142 write(errmsg,'(*(a))') 'Logic error: namelist choice of', &
143 & ' deep convection is different from new Tiedtke scheme'
144 errflg = 1
145 return
146 end if
147
148 end subroutine cu_ntiedtke_init
149
150! Tiedtke cumulus scheme from WRF with small modifications
151! This scheme includes both deep and shallow convections
152!===================
153!
157!-----------------------------------------------------------------------
158! level 1 subroutine 'tiecnvn'
159!-----------------------------------------------------------------
160 subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, &
161 evap,hfx,zprecc,lmask,lq,km,dt,dx,kbot,ktop,kcnv, &
162 ktrac,ud_mf,dd_mf,dt_mf,cnvw,cnvc,errmsg,errflg)
163!-----------------------------------------------------------------
164! this is the interface between the model and the mass
165! flux convection module
166!-----------------------------------------------------------------
167 implicit none
168! in&out variables
169 integer, intent(in) :: lq, km, ktrac
170 real(kind=kind_phys), intent(in ) :: dt
171 integer, dimension( : ), intent(in) :: lmask
172 real(kind=kind_phys), dimension( : ), intent(in ) :: evap, hfx, dx
173 real(kind=kind_phys), dimension( :, : ), intent(inout) :: pu, pv, pt, pqv
174 real(kind=kind_phys), dimension( :, :), intent(in ) :: tdi, qvdi, poz, prsl, pomg
175 real(kind=kind_phys), dimension( :, :), intent(in ), optional :: pqvf, ptf
176 real(kind=kind_phys), dimension( :, : ), intent(in ) :: pzz, prsi
177 real(kind=kind_phys), dimension( :, :, : ), intent(inout ) :: clw
178
179 integer, dimension( : ), intent(out) :: kbot, ktop, kcnv
180 real(kind=kind_phys), dimension( : ), intent(out) :: zprecc
181 real(kind=kind_phys), dimension (:, :), intent(out), optional :: ud_mf
182 real(kind=kind_phys), dimension (:, :), intent(out) :: dd_mf, dt_mf, cnvw, cnvc
183
184! error messages
185 character(len=*), intent(out) :: errmsg
186 integer, intent(out) :: errflg
187
188! local variables
189 real(kind=kind_phys) pum1(lq,km), pvm1(lq,km), ztt(lq,km), &
190 & ptte(lq,km), pqte(lq,km), pvom(lq,km), pvol(lq,km), &
191 & pverv(lq,km), pgeo(lq,km), pap(lq,km), paph(lq,km+1)
192 real(kind=kind_phys) pqhfl(lq), zqq(lq,km), &
193 & prsfc(lq), pssfc(lq), pcte(lq,km), &
194 & phhfl(lq), pgeoh(lq,km+1)
195 real(kind=kind_phys) ztp1(lq,km), zqp1(lq,km), ztu(lq,km), zqu(lq,km),&
196 & zlu(lq,km), zlude(lq,km), zmfu(lq,km), zmfd(lq,km), zmfude_rate(lq,km),&
197 & zqsat(lq,km), zrain(lq)
198 real(kind=kind_phys),allocatable :: pcen(:,:,:),ptenc(:,:,:)
199
200 integer icbot(lq), ictop(lq), ktype(lq), lndj(lq)
201 logical locum(lq)
202!
203 real(kind=kind_phys) ztmst,fliq,fice,ztc,zalf,tt
204 integer i,j,k,k1,n,km1,ktracer
205 real(kind=kind_phys) ztpp1
206 real(kind=kind_phys) zew,zqs,zcor
207!
208! Initialize CCPP error handling variables
209 errmsg = ''
210 errflg = 0
211
212 km1 = km + 1
213 ztmst=dt
214!
215! masv flux diagnostics.
216!
217 do j=1,lq
218 zrain(j)=0.0
219 locum(j)=.false.
220 prsfc(j)=0.0
221 pssfc(j)=0.0
222 pqhfl(j)=evap(j)
223 phhfl(j)=hfx(j)
224 pgeoh(j,km1)=pzz(j,1)
225 paph(j,km1)=prsi(j,1)
226 if(lmask(j).eq.1) then
227 lndj(j)=1
228 else
229 lndj(j)=0
230 end if
231 end do
232!
233! convert model variables for mflux scheme
234!
235 do k=1,km
236 k1=km-k+1
237 do j=1,lq
238 pcte(j,k1)=0.0
239 pvom(j,k1)=0.0
240 pvol(j,k1)=0.0
241 ztp1(j,k1)=pt(j,k)
242 zqp1(j,k1)=pqv(j,k)
243 pum1(j,k1)=pu(j,k)
244 pvm1(j,k1)=pv(j,k)
245 pverv(j,k1)=pomg(j,k)
246 pgeo(j,k1)=poz(j,k)
247 pgeoh(j,k1)=pzz(j,k+1)
248 pap(j,k1)=prsl(j,k)
249 paph(j,k1)=prsi(j,k+1)
250 tt=ztp1(j,k1)
251 zew = foeewm(tt)
252 zqs = zew/pap(j,k1)
253 zqs = min(0.5,zqs)
254 zcor = 1./(1.-vtmpc1*zqs)
255 zqsat(j,k1)=zqs*zcor
256 pqte(j,k1)=pqvf(j,k)+(pqv(j,k)-qvdi(j,k))/ztmst
257 zqq(j,k1) =pqte(j,k1)
258 ptte(j,k1)=ptf(j,k)+(pt(j,k)-tdi(j,k))/ztmst
259 ztt(j,k1) =ptte(j,k1)
260 ud_mf(j,k1)=0.
261 dd_mf(j,k1)=0.
262 dt_mf(j,k1)=0.
263 cnvw(j,k1)=0.
264 cnvc(j,k1)=0.
265 end do
266 end do
267
268 if(ktrac > 2) then
269 ktracer = ktrac - 2
270 allocate(pcen(lq,km,ktracer))
271 allocate(ptenc(lq,km,ktracer))
272 do n=1,ktracer
273 do k=1,km
274 k1=km-k+1
275 do j=1,lq
276 pcen(j,k1,n) = clw(j,k,n+2)
277 ptenc(j,k1,n)= 0.
278 end do
279 end do
280 end do
281 else
282 ktracer = 2
283 allocate(pcen(lq,km,ktracer))
284 allocate(ptenc(lq,km,ktracer))
285 do n=1,ktracer
286 do k=1,km
287 do j=1,lq
288 pcen(j,k,n) = 0.
289 ptenc(j,k,n)= 0.
290 end do
291 end do
292 end do
293 end if
294
295! print *, "pgeo=",pgeo(1,:)
296! print *, "pgeoh=",pgeoh(1,:)
297! print *, "pap=",pap(1,:)
298! print *, "paph=",paph(1,:)
299! print *, "ztp1=",ztp1(1,:)
300! print *, "zqp1=",zqp1(1,:)
301! print *, "pum1=",pum1(1,:)
302! print *, "pvm1=",pvm1(1,:)
303! print *, "pverv=",pverv(1,:)
304! print *, "pqte=",pqte(1,:)
305! print *, "ptte=",ptte(1,:)
306! print *, "hfx=", pqhfl(1),phhfl(1),dx(1)
307!
308!-----------------------------------------------------------------------
309!* 2. call 'cumastrn'(master-routine for cumulus parameterization)
310!
311 call cumastrn &
312 & (lq, km, km1, km-1, ztp1, &
313 & zqp1, pum1, pvm1, pverv, zqsat,&
314 & pqhfl, ztmst, pap, paph, pgeo, &
315 & ptte, pqte, pvom, pvol, prsfc,&
316 & pssfc, locum, ktracer, pcen, ptenc,&
317 & ktype, icbot, ictop, ztu, zqu, &
318 & zlu, zlude, zmfu, zmfd, zrain,&
319 & pcte, phhfl, lndj, pgeoh, zmfude_rate, dx)
320!
321! to include the cloud water and cloud ice detrained from convection
322!
323 do k=1,km
324 k1=km-k+1
325 do j=1,lq
326 if(pcte(j,k1).gt.0.) then
327 fliq=foealfa(ztp1(j,k1))
328 fice=1.0-fliq
329 clw(j,k,2)=clw(j,k,2)+fliq*pcte(j,k1)*ztmst
330 clw(j,k,1)=clw(j,k,1)+fice*pcte(j,k1)*ztmst
331 endif
332 end do
333 end do
334!
335 do k=1,km
336 k1 = km-k+1
337 do j=1,lq
338 pt(j,k) = ztp1(j,k1)+(ptte(j,k1)-ztt(j,k1))*ztmst
339 pqv(j,k)= zqp1(j,k1)+(pqte(j,k1)-zqq(j,k1))*ztmst
340 ud_mf(j,k)= zmfu(j,k1)*ztmst
341 dd_mf(j,k)= -zmfd(j,k1)*ztmst
342 dt_mf(j,k)= zmfude_rate(j,k1)*ztmst
343 cnvw(j,k) = zlude(j,k1)*ztmst*g/(prsi(j,k)-prsi(j,k+1))
344 cnvc(j,k) = 0.04 * log(1. + 675. * ud_mf(j,k))
345 cnvc(j,k) = min(cnvc(j,k), 0.6)
346 cnvc(j,k) = max(cnvc(j,k), 0.0)
347 end do
348 end do
349
350 do j=1,lq
351 zprecc(j)=amax1(0.0,(prsfc(j)+pssfc(j))*ztmst*0.001)
352 kbot(j) = km-icbot(j)+1
353 ktop(j) = km-ictop(j)+1
354 if(ktype(j).eq.1 .or. ktype(j).eq.3) then
355 kcnv(j)=1
356 else
357 kcnv(j)=0
358 end if
359 end do
360
361 if (lmfdudv) then
362 do k=1,km
363 k1=km-k+1
364 do j=1,lq
365 pu(j,k)=pu(j,k)+pvom(j,k1)*ztmst
366 pv(j,k)=pv(j,k)+pvol(j,k1)*ztmst
367 end do
368 end do
369 endif
370
371!
372! Currently, vertical mixing of tracers are turned off
373! if(ktrac > 2) then
374! do n=1,ktrac-2
375! do k=1,km
376! k1=km-k+1
377! do j=1,lq
378! clw(j,k,n+2)=pcen(j,k,n)+ptenc(j,k1,n)*ztmst
379! end do
380! end do
381! end do
382! end if
383 deallocate(pcen)
384 deallocate(ptenc)
385!
386 return
387 end subroutine cu_ntiedtke_run
388
389!#############################################################
390!
391! level 2 subroutines
392!
393!#############################################################
394!***********************************************************
395! subroutine cumastrn
396!***********************************************************
397 subroutine cumastrn &
398 & (klon, klev, klevp1, klevm1, pten, &
399 & pqen, puen, pven, pverv, pqsen,&
400 & pqhfl, ztmst, pap, paph, pgeo, &
401 & ptte, pqte, pvom, pvol, prsfc,&
402 & pssfc, ldcum, ktrac, pcen, ptenc,&
403 & ktype, kcbot, kctop, ptu, pqu,&
404 & plu, plude, pmfu, pmfd, prain,&
405 & pcte, phhfl, lndj, zgeoh, pmfude_rate, dx)
406 implicit none
407!
408!***cumastrn* master routine for cumulus massflux-scheme
409! m.tiedtke e.c.m.w.f. 1986/1987/1989
410! modifications
411! y.wang i.p.r.c 2001
412! c.zhang 2012
413!***purpose
414! -------
415! this routine computes the physical tendencies of the
416! prognostic variables t,q,u and v due to convective processes.
417! processes considered are: convective fluxes, formation of
418! precipitation, evaporation of falling rain below cloud base,
419! saturated cumulus downdrafts.
420!***method
421! ------
422! parameterization is done using a massflux-scheme.
423! (1) define constants and parameters
424! (2) specify values (t,q,qs...) at half levels and
425! initialize updraft- and downdraft-values in 'cuinin'
426! (3) calculate cloud base in 'cutypen', calculate cloud types in cutypen,
427! and specify cloud base massflux
428! (4) do cloud ascent in 'cuascn' in absence of downdrafts
429! (5) do downdraft calculations:
430! (a) determine values at lfs in 'cudlfsn'
431! (b) determine moist descent in 'cuddrafn'
432! (c) recalculate cloud base massflux considering the
433! effect of cu-downdrafts
434! (6) do final adjusments to convective fluxes in 'cuflxn',
435! do evaporation in subcloud layer
436! (7) calculate increments of t and q in 'cudtdqn'
437! (8) calculate increments of u and v in 'cududvn'
438!***externals.
439! ----------
440! cuinin: initializes values at vertical grid used in cu-parametr.
441! cutypen: cloud bypes, 1: deep cumulus 2: shallow cumulus
442! cuascn: cloud ascent for entraining plume
443! cudlfsn: determines values at lfs for downdrafts
444! cuddrafn:does moist descent for cumulus downdrafts
445! cuflxn: final adjustments to convective fluxes (also in pbl)
446! cudqdtn: updates tendencies for t and q
447! cududvn: updates tendencies for u and v
448!***switches.
449! --------
450! lmfmid=.t. midlevel convection is switched on
451! lmfdd=.t. cumulus downdrafts switched on
452! lmfdudv=.t. cumulus friction switched on
453!***
454! model parameters (defined in subroutine cuparam)
455! ------------------------------------------------
456! entrdd entrainment rate for cumulus downdrafts
457! cmfcmax maximum massflux value allowed for
458! cmfcmin minimum massflux value (for safety)
459! cmfdeps fractional massflux for downdrafts at lfs
460! cprcon coefficient for conversion from cloud water to rain
461!***reference.
462! ----------
463! paper on massflux scheme (tiedtke,1989)
464!-----------------------------------------------------------------
465 integer klev,klon,ktrac,klevp1,klevm1
466 real(kind=kind_phys) pten(klon,klev), pqen(klon,klev),&
467 & puen(klon,klev), pven(klon,klev),&
468 & ptte(klon,klev), pqte(klon,klev),&
469 & pvom(klon,klev), pvol(klon,klev),&
470 & pqsen(klon,klev), pgeo(klon,klev),&
471 & pap(klon,klev), paph(klon,klevp1),&
472 & pverv(klon,klev), pqhfl(klon),&
473 & phhfl(klon)
474 real(kind=kind_phys) ptu(klon,klev), pqu(klon,klev),&
475 & plu(klon,klev), plude(klon,klev),&
476 & pmfu(klon,klev), pmfd(klon,klev),&
477 & prain(klon),&
478 & prsfc(klon), pssfc(klon)
479 real(kind=kind_phys) ztenh(klon,klev), zqenh(klon,klev),&
480 & zgeoh(klon,klevp1), zqsenh(klon,klev),&
481 & ztd(klon,klev), zqd(klon,klev),&
482 & zmfus(klon,klev), zmfds(klon,klev),&
483 & zmfuq(klon,klev), zmfdq(klon,klev),&
484 & zdmfup(klon,klev), zdmfdp(klon,klev),&
485 & zmful(klon,klev), zrfl(klon),&
486 & zuu(klon,klev), zvu(klon,klev),&
487 & zud(klon,klev), zvd(klon,klev),&
488 & zlglac(klon,klev)
489 real(kind=kind_phys) pmflxr(klon,klevp1), pmflxs(klon,klevp1)
490 real(kind=kind_phys) zhcbase(klon),&
491 & zmfub(klon), zmfub1(klon),&
492 & zdhpbl(klon)
493 real(kind=kind_phys) zsfl(klon), zdpmel(klon,klev),&
494 & pcte(klon,klev), zcape(klon),&
495 & zcape1(klon), zcape2(klon),&
496 & ztauc(klon), ztaubl(klon),&
497 & zheat(klon)
498 real(kind=kind_phys) pcen(klon,klev,ktrac), ptenc(klon,klev,ktrac)
499 real(kind=kind_phys) wup(klon), zdqcv(klon)
500 real(kind=kind_phys) wbase(klon), zmfuub(klon)
501 real(kind=kind_phys) upbl(klon)
502 real(kind=kind_phys) dx(klon)
503 real(kind=kind_phys) pmfude_rate(klon,klev), pmfdde_rate(klon,klev)
504 real(kind=kind_phys) zmfuus(klon,klev), zmfdus(klon,klev)
505 real(kind=kind_phys) zmfudr(klon,klev), zmfddr(klon,klev)
506 real(kind=kind_phys) zuv2(klon,klev),ztenu(klon,klev),ztenv(klon,klev)
507 real(kind=kind_phys) zmfuvb(klon),zsum12(klon),zsum22(klon)
508 integer ilab(klon,klev), idtop(klon),&
509 & ictop0(klon), ilwmin(klon)
510 integer kdpl(klon)
511 integer kcbot(klon), kctop(klon),&
512 & ktype(klon), lndj(klon)
513 logical ldcum(klon), lldcum(klon)
514 logical loddraf(klon), llddraf3(klon), llo1, llo2(klon)
515
516! local varaiables
517 real(kind=kind_phys) zcons,zcons2,zqumqe,zdqmin,zdh,zmfmax
518 real(kind=kind_phys) zalfaw,zalv,zqalv,zc5ldcp,zc4les,zhsat,zgam,zzz,zhhat
519 real(kind=kind_phys) zpbmpt,zro,zdz,zdp,zeps,zfac,wspeed
520 integer jl,jk,ik
521 integer ikb,ikt,icum,itopm2
522 real(kind=kind_phys) ztmst,ztau,zerate,zderate,zmfa
523 real(kind=kind_phys) zmfs(klon),pmean(klev),zlon
524 real(kind=kind_phys) zduten,zdvten,ztdis,pgf_u,pgf_v
525!-------------------------------------------
526! 1. specify constants and parameters
527!-------------------------------------------
528 zcons=1./(g*ztmst)
529 zcons2=3./(g*ztmst)
530
531 zlon = real(klon)
532 do jk = klev , 1 , -1
533 pmean(jk) = sum(pap(:,jk))/zlon
534 end do
535 p650 = klev-2
536 do jk = klev , 3 , -1
537 if ( pmean(jk)/pmean(klev)*1.013250e5 > 650.e2 ) p650 = jk
538 end do
539
540!--------------------------------------------------------------
541!* 2. initialize values at vertical grid points in 'cuini'
542!--------------------------------------------------------------
543 call cuinin &
544 & (klon, klev, klevp1, klevm1, pten, &
545 & pqen, pqsen, puen, pven, pverv,&
546 & pgeo, paph, zgeoh, ztenh, zqenh,&
547 & zqsenh, ilwmin, ptu, pqu, ztd, &
548 & zqd, zuu, zvu, zud, zvd, &
549 & pmfu, pmfd, zmfus, zmfds, zmfuq,&
550 & zmfdq, zdmfup, zdmfdp, zdpmel, plu, &
551 & plude, ilab)
552
553!----------------------------------
554!* 3.0 cloud base calculations
555!----------------------------------
556!* (a) determine cloud base values in 'cutypen',
557! and the cumulus type 1 or 2
558! -------------------------------------------
559 call cutypen &
560 & ( klon, klev, klevp1, klevm1, pqen,&
561 & ztenh, zqenh, zqsenh, zgeoh, paph,&
562 & phhfl, pqhfl, pgeo, pqsen, pap,&
563 & pten, lndj, ptu, pqu, ilab,&
564 & ldcum, kcbot, ictop0, ktype, wbase, plu, kdpl)
565
566!* (b) assign the first guess mass flux at cloud base
567! ------------------------------------------
568 do jl=1,klon
569 zdhpbl(jl)=0.0
570 upbl(jl) = 0.0
571 idtop(jl)=0
572 end do
573
574 do jk=2,klev
575 do jl=1,klon
576 if(jk.ge.kcbot(jl) .and. ldcum(jl)) then
577 zdhpbl(jl)=zdhpbl(jl)+(alv*pqte(jl,jk)+cpd*ptte(jl,jk))&
578 & *(paph(jl,jk+1)-paph(jl,jk))
579 if(lndj(jl) .eq. 0) then
580 wspeed = sqrt(puen(jl,jk)**2 + pven(jl,jk)**2)
581 upbl(jl) = upbl(jl) + wspeed*(paph(jl,jk+1)-paph(jl,jk))
582 end if
583 end if
584 end do
585 end do
586
587 do jl=1,klon
588 if(ldcum(jl)) then
589 ikb=kcbot(jl)
590 zmfmax = (paph(jl,ikb)-paph(jl,ikb-1))*zcons2
591 if(ktype(jl) == 1) then
592 zmfub(jl)= 0.1*zmfmax
593 else if ( ktype(jl) == 2 ) then
594 zqumqe = pqu(jl,ikb) + plu(jl,ikb) - zqenh(jl,ikb)
595 zdqmin = max(0.01*zqenh(jl,ikb),1.e-10)
596 zdh = cpd*(ptu(jl,ikb)-ztenh(jl,ikb)) + alv*zqumqe
597 zdh = g*max(zdh,1.e5*zdqmin)
598 if ( zdhpbl(jl) > 0. ) then
599 zmfub(jl) = zdhpbl(jl)/zdh
600 zmfub(jl) = min(zmfub(jl),zmfmax)
601 else
602 zmfub(jl) = 0.1*zmfmax
603 ldcum(jl) = .false.
604 end if
605 end if
606 else
607 zmfub(jl) = 0.
608 end if
609 end do
610!------------------------------------------------------
611!* 4.0 determine cloud ascent for entraining plume
612!------------------------------------------------------
613!* (a) do ascent in 'cuasc'in absence of downdrafts
614!----------------------------------------------------------
615 call cuascn &
616 & (klon, klev, klevp1, klevm1, ztenh,&
617 & zqenh, puen, pven, pten, pqen,&
618 & pqsen, pgeo, zgeoh, pap, paph,&
619 & pqte, pverv, ilwmin, ldcum, zhcbase,&
620 & ktype, ilab, ptu, pqu, plu,&
621 & zuu, zvu, pmfu, zmfub,&
622 & zmfus, zmfuq, zmful, plude, zdmfup,&
623 & kcbot, kctop, ictop0, icum, ztmst,&
624 & zqsenh, zlglac, lndj, wup, wbase, kdpl, pmfude_rate )
625
626!* (b) check cloud depth and change entrainment rate accordingly
627! calculate precipitation rate (for downdraft calculation)
628!------------------------------------------------------------------
629 do jl=1,klon
630 if ( ldcum(jl) ) then
631 ikb = kcbot(jl)
632 itopm2 = kctop(jl)
633 zpbmpt = paph(jl,ikb) - paph(jl,itopm2)
634 if ( ktype(jl) == 1 .and. zpbmpt < zdnoprc ) ktype(jl) = 2
635 if ( ktype(jl) == 2 .and. zpbmpt >= zdnoprc ) ktype(jl) = 1
636 ictop0(jl) = kctop(jl)
637 end if
638 zrfl(jl)=zdmfup(jl,1)
639 end do
640
641 do jk=2,klev
642 do jl=1,klon
643 zrfl(jl)=zrfl(jl)+zdmfup(jl,jk)
644 end do
645 end do
646
647 do jk = 1,klev
648 do jl = 1,klon
649 pmfd(jl,jk) = 0.
650 zmfds(jl,jk) = 0.
651 zmfdq(jl,jk) = 0.
652 zdmfdp(jl,jk) = 0.
653 zdpmel(jl,jk) = 0.
654 end do
655 end do
656
657!-----------------------------------------
658!* 6.0 cumulus downdraft calculations
659!-----------------------------------------
660 if(lmfdd) then
661!* (a) determine lfs in 'cudlfsn'
662!--------------------------------------
663 call cudlfsn &
664 & (klon, klev,&
665 & kcbot, kctop, lndj, ldcum, &
666 & ztenh, zqenh, puen, pven, &
667 & pten, pqsen, pgeo, &
668 & zgeoh, paph, ptu, pqu, plu, &
669 & zuu, zvu, zmfub, zrfl, &
670 & ztd, zqd, zud, zvd, &
671 & pmfd, zmfds, zmfdq, zdmfdp, &
672 & idtop, loddraf)
673!* (b) determine downdraft t,q and fluxes in 'cuddrafn'
674!------------------------------------------------------------
675 call cuddrafn &
676 & ( klon, klev, loddraf, &
677 & ztenh, zqenh, puen, pven, &
678 & pgeo, zgeoh, paph, zrfl, &
679 & ztd, zqd, zud, zvd, pmfu, &
680 & pmfd, zmfds, zmfdq, zdmfdp, pmfdde_rate )
681!-----------------------------------------------------------
682 end if
683!
684!-----------------------------------------------------------------------
685!* 6.0 closure and clean work
686! ------
687!-- 6.1 recalculate cloud base massflux from a cape closure
688! for deep convection (ktype=1)
689!
690 do jl=1,klon
691 if(ldcum(jl) .and. ktype(jl) .eq. 1) then
692 ikb = kcbot(jl)
693 ikt = kctop(jl)
694 zheat(jl)=0.0
695 zcape(jl)=0.0
696 zcape1(jl)=0.0
697 zcape2(jl)=0.0
698 zmfub1(jl)=zmfub(jl)
699
700 ztauc(jl) = (zgeoh(jl,ikt)-zgeoh(jl,ikb)) / &
701 ((2.+ min(15.0,wup(jl)))*g)
702 if(lndj(jl) .eq. 0) then
703 upbl(jl) = 2.+ upbl(jl)/(paph(jl,klev+1)-paph(jl,ikb))
704 ztaubl(jl) = (zgeoh(jl,ikb)-zgeoh(jl,klev+1))/(g*upbl(jl))
705 ztaubl(jl) = min(300., ztaubl(jl))
706 else
707 ztaubl(jl) = ztauc(jl)
708 end if
709 end if
710 end do
711!
712 do jk = 1 , klev
713 do jl = 1 , klon
714 llo1 = ldcum(jl) .and. ktype(jl) .eq. 1
715 if ( llo1 .and. jk <= kcbot(jl) .and. jk > kctop(jl) ) then
716 ikb = kcbot(jl)
717 zdz = pgeo(jl,jk-1)-pgeo(jl,jk)
718 zdp = pap(jl,jk)-pap(jl,jk-1)
719 zheat(jl) = zheat(jl) + ((pten(jl,jk-1)-pten(jl,jk)+zdz*rcpd) / &
720 ztenh(jl,jk)+vtmpc1*(pqen(jl,jk-1)-pqen(jl,jk))) * &
721 (g*(pmfu(jl,jk)+pmfd(jl,jk)))
722 zcape1(jl) = zcape1(jl) + ((ptu(jl,jk)-ztenh(jl,jk))/ztenh(jl,jk) + &
723 vtmpc1*(pqu(jl,jk)-zqenh(jl,jk))-plu(jl,jk))*zdp
724 end if
725
726 if ( llo1 .and. jk >= kcbot(jl) ) then
727 if((paph(jl,klev+1)-paph(jl,kdpl(jl)))<50.e2) then
728 zdp = paph(jl,jk+1)-paph(jl,jk)
729 zcape2(jl) = zcape2(jl) + ztaubl(jl)* &
730 ((1.+vtmpc1*pqen(jl,jk))*ptte(jl,jk)+vtmpc1*pten(jl,jk)*pqte(jl,jk))*zdp
731 end if
732 end if
733 end do
734 end do
735
736 do jl=1,klon
737 if(ldcum(jl).and.ktype(jl).eq.1) then
738 ikb = kcbot(jl)
739 ikt = kctop(jl)
740 ztau = ztauc(jl) * (1.+1.33e-5*dx(jl))
741 ztau = max(ztmst,ztau)
742 ztau = max(720.,ztau)
743 ztau = min(10800.,ztau)
744 if(isequil) then
745 zcape2(jl)= max(0.,zcape2(jl))
746 zcape(jl) = max(0.,min(zcape1(jl)-zcape2(jl),5000.))
747 else
748 zcape(jl) = max(0.,min(zcape1(jl),5000.))
749 end if
750 zheat(jl) = max(1.e-4,zheat(jl))
751 zmfub1(jl) = (zcape(jl)*zmfub(jl))/(zheat(jl)*ztau)
752 zmfub1(jl) = max(zmfub1(jl),0.001)
753 zmfmax=(paph(jl,ikb)-paph(jl,ikb-1))*zcons2
754 zmfub1(jl)=min(zmfub1(jl),zmfmax)
755 end if
756 end do
757!
758!* 6.2 recalculate convective fluxes due to effect of
759! downdrafts on boundary layer moist static energy budget (ktype=2)
760!--------------------------------------------------------
761 do jl=1,klon
762 if(ldcum(jl) .and. ktype(jl) .eq. 2) then
763 ikb=kcbot(jl)
764 if(pmfd(jl,ikb).lt.0.0 .and. loddraf(jl)) then
765 zeps=-pmfd(jl,ikb)/max(zmfub(jl),cmfcmin)
766 else
767 zeps=0.
768 endif
769 zqumqe=pqu(jl,ikb)+plu(jl,ikb)- &
770 & zeps*zqd(jl,ikb)-(1.-zeps)*zqenh(jl,ikb)
771 zdqmin=max(0.01*zqenh(jl,ikb),cmfcmin)
772 zmfmax=(paph(jl,ikb)-paph(jl,ikb-1))*zcons2
773! using moist static engergy closure instead of moisture closure
774 zdh=cpd*(ptu(jl,ikb)-zeps*ztd(jl,ikb)- &
775 & (1.-zeps)*ztenh(jl,ikb))+alv*zqumqe
776 zdh=g*max(zdh,1.e5*zdqmin)
777 if(zdhpbl(jl).gt.0.)then
778 zmfub1(jl)=zdhpbl(jl)/zdh
779 else
780 zmfub1(jl) = zmfub(jl)
781 end if
782 zmfub1(jl) = min(zmfub1(jl),zmfmax)
783 end if
784
785!* 6.3 mid-level convection - nothing special
786!---------------------------------------------------------
787 if(ldcum(jl) .and. ktype(jl) .eq. 3 ) then
788 zmfub1(jl) = zmfub(jl)
789 end if
790
791 end do
792
793!* 6.4 scaling the downdraft mass flux
794!---------------------------------------------------------
795 do jk=1,klev
796 do jl=1,klon
797 if( ldcum(jl) ) then
798 zfac=zmfub1(jl)/max(zmfub(jl),cmfcmin)
799 pmfd(jl,jk)=pmfd(jl,jk)*zfac
800 zmfds(jl,jk)=zmfds(jl,jk)*zfac
801 zmfdq(jl,jk)=zmfdq(jl,jk)*zfac
802 zdmfdp(jl,jk)=zdmfdp(jl,jk)*zfac
803 pmfdde_rate(jl,jk) = pmfdde_rate(jl,jk)*zfac
804 end if
805 end do
806 end do
807
808!* 6.5 scaling the updraft mass flux
809! --------------------------------------------------------
810 do jl = 1,klon
811 if ( ldcum(jl) ) zmfs(jl) = zmfub1(jl)/max(cmfcmin,zmfub(jl))
812 end do
813 do jk = 2 , klev
814 do jl = 1,klon
815 if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
816 ikb = kcbot(jl)
817 if ( jk>ikb ) then
818 zdz = ((paph(jl,klev+1)-paph(jl,jk))/(paph(jl,klev+1)-paph(jl,ikb)))
819 pmfu(jl,jk) = pmfu(jl,ikb)*zdz
820 end if
821 zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons2
822 if ( pmfu(jl,jk)*zmfs(jl) > zmfmax ) then
823 zmfs(jl) = min(zmfs(jl),zmfmax/pmfu(jl,jk))
824 end if
825 end if
826 end do
827 end do
828 do jk = 2 , klev
829 do jl = 1,klon
830 if ( ldcum(jl) .and. jk <= kcbot(jl) .and. jk >= kctop(jl)-1 ) then
831 pmfu(jl,jk) = pmfu(jl,jk)*zmfs(jl)
832 zmfus(jl,jk) = zmfus(jl,jk)*zmfs(jl)
833 zmfuq(jl,jk) = zmfuq(jl,jk)*zmfs(jl)
834 zmful(jl,jk) = zmful(jl,jk)*zmfs(jl)
835 zdmfup(jl,jk) = zdmfup(jl,jk)*zmfs(jl)
836 plude(jl,jk) = plude(jl,jk)*zmfs(jl)
837 pmfude_rate(jl,jk) = pmfude_rate(jl,jk)*zmfs(jl)
838 end if
839 end do
840 end do
841
842!* 6.6 if ktype = 2, kcbot=kctop is not allowed
843! ---------------------------------------------------
844 do jl = 1,klon
845 if ( ktype(jl) == 2 .and. &
846 kcbot(jl) == kctop(jl) .and. kcbot(jl) >= klev-1 ) then
847 ldcum(jl) = .false.
848 ktype(jl) = 0
849 end if
850 end do
851
852 if ( .not. lmfscv .or. .not. lmfpen ) then
853 do jl = 1,klon
854 llo2(jl) = .false.
855 if ( (.not. lmfscv .and. ktype(jl) == 2) .or. &
856 (.not. lmfpen .and. ktype(jl) == 1) ) then
857 llo2(jl) = .true.
858 ldcum(jl) = .false.
859 end if
860 end do
861 end if
862
863!* 6.7 set downdraft mass fluxes to zero above cloud top
864!----------------------------------------------------
865 do jl = 1,klon
866 if ( loddraf(jl) .and. idtop(jl) <= kctop(jl) ) then
867 idtop(jl) = kctop(jl) + 1
868 end if
869 end do
870 do jk = 2 , klev
871 do jl = 1,klon
872 if ( loddraf(jl) ) then
873 if ( jk < idtop(jl) ) then
874 pmfd(jl,jk) = 0.
875 zmfds(jl,jk) = 0.
876 zmfdq(jl,jk) = 0.
877 pmfdde_rate(jl,jk) = 0.
878 zdmfdp(jl,jk) = 0.
879 else if ( jk == idtop(jl) ) then
880 pmfdde_rate(jl,jk) = 0.
881 end if
882 end if
883 end do
884 end do
885
886 itopm2 = 2
887!----------------------------------------------------------
888!* 7.0 determine final convective fluxes in 'cuflx'
889!----------------------------------------------------------
890 call cuflxn &
891 & ( klon, klev, ztmst &
892 & , pten, pqen, pqsen, ztenh, zqenh &
893 & , paph, pap, zgeoh, lndj, ldcum &
894 & , kcbot, kctop, idtop, itopm2 &
895 & , ktype, loddraf &
896 & , pmfu, pmfd, zmfus, zmfds &
897 & , zmfuq, zmfdq, zmful, plude &
898 & , zdmfup, zdmfdp, zdpmel, zlglac &
899 & , prain, pmfdde_rate, pmflxr, pmflxs )
900
901! some adjustments needed
902 do jl=1,klon
903 zmfs(jl) = 1.
904 zmfuub(jl)=0.
905 end do
906 do jk = 2 , klev
907 do jl = 1,klon
908 if ( loddraf(jl) .and. jk >= idtop(jl)-1 ) then
909 zmfmax = pmfu(jl,jk)*0.98
910 if ( pmfd(jl,jk)+zmfmax+1.e-15 < 0. ) then
911 zmfs(jl) = min(zmfs(jl),-zmfmax/pmfd(jl,jk))
912 end if
913 end if
914 end do
915 end do
916
917 do jk = 2 , klev
918 do jl = 1 , klon
919 if ( zmfs(jl) < 1. .and. jk >= idtop(jl)-1 ) then
920 pmfd(jl,jk) = pmfd(jl,jk)*zmfs(jl)
921 zmfds(jl,jk) = zmfds(jl,jk)*zmfs(jl)
922 zmfdq(jl,jk) = zmfdq(jl,jk)*zmfs(jl)
923 pmfdde_rate(jl,jk) = pmfdde_rate(jl,jk)*zmfs(jl)
924 zmfuub(jl) = zmfuub(jl) - (1.-zmfs(jl))*zdmfdp(jl,jk)
925 pmflxr(jl,jk+1) = pmflxr(jl,jk+1) + zmfuub(jl)
926 zdmfdp(jl,jk) = zdmfdp(jl,jk)*zmfs(jl)
927 end if
928 end do
929 end do
930
931 do jk = 2 , klev - 1
932 do jl = 1, klon
933 if ( loddraf(jl) .and. jk >= idtop(jl)-1 ) then
934 zerate = -pmfd(jl,jk) + pmfd(jl,jk-1) + pmfdde_rate(jl,jk)
935 if ( zerate < 0. ) then
936 pmfdde_rate(jl,jk) = pmfdde_rate(jl,jk) - zerate
937 end if
938 end if
939 if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
940 zerate = pmfu(jl,jk) - pmfu(jl,jk+1) + pmfude_rate(jl,jk)
941 if ( zerate < 0. ) then
942 pmfude_rate(jl,jk) = pmfude_rate(jl,jk) - zerate
943 end if
944 zdmfup(jl,jk) = pmflxr(jl,jk+1) + pmflxs(jl,jk+1) - &
945 pmflxr(jl,jk) - pmflxs(jl,jk)
946 zdmfdp(jl,jk) = 0.
947 end if
948 end do
949 end do
950
951! avoid negative humidities at ddraught top
952 do jl = 1,klon
953 if ( loddraf(jl) ) then
954 jk = idtop(jl)
955 ik = min(jk+1,klev)
956 if ( zmfdq(jl,jk) < 0.3*zmfdq(jl,ik) ) then
957 zmfdq(jl,jk) = 0.3*zmfdq(jl,ik)
958 end if
959 end if
960 end do
961
962! avoid negative humidities near cloud top because gradient of precip flux
963! and detrainment / liquid water flux are too large
964 do jk = 2 , klev
965 do jl = 1, klon
966 if ( ldcum(jl) .and. jk >= kctop(jl)-1 .and. jk < kcbot(jl) ) then
967 zdz = ztmst*g/(paph(jl,jk+1)-paph(jl,jk))
968 zmfa = zmfuq(jl,jk+1) + zmfdq(jl,jk+1) - &
969 zmfuq(jl,jk) - zmfdq(jl,jk) + &
970 zmful(jl,jk+1) - zmful(jl,jk) + zdmfup(jl,jk)
971 zmfa = (zmfa-plude(jl,jk))*zdz
972 if ( pqen(jl,jk)+zmfa < 0. ) then
973 plude(jl,jk) = plude(jl,jk) + 2.*(pqen(jl,jk)+zmfa)/zdz
974 end if
975 if ( plude(jl,jk) < 0. ) plude(jl,jk) = 0.
976 end if
977 if ( .not. ldcum(jl) ) pmfude_rate(jl,jk) = 0.
978 if ( abs(pmfd(jl,jk-1)) < 1.0e-20 ) pmfdde_rate(jl,jk) = 0.
979 end do
980 end do
981
982 do jl=1,klon
983 prsfc(jl) = pmflxr(jl,klev+1)
984 pssfc(jl) = pmflxs(jl,klev+1)
985 end do
986
987!----------------------------------------------------------------
988!* 8.0 update tendencies for t and q in subroutine cudtdq
989!----------------------------------------------------------------
990 call cudtdqn(klon,klev,itopm2,kctop,idtop,ldcum,loddraf, &
991 ztmst,paph,zgeoh,pgeo,pten,ztenh,pqen,zqenh,pqsen, &
992 zlglac,plude,pmfu,pmfd,zmfus,zmfds,zmfuq,zmfdq,zmful, &
993 zdmfup,zdmfdp,zdpmel,ptte,pqte,pcte)
994!----------------------------------------------------------------
995!* 9.0 update tendencies for u and u in subroutine cududv
996!----------------------------------------------------------------
997 if(lmfdudv) then
998 do jk = klev-1 , 2 , -1
999 ik = jk + 1
1000 do jl = 1,klon
1001 if ( ldcum(jl) ) then
1002 if ( jk == kcbot(jl) .and. ktype(jl) < 3 ) then
1003 ikb = kdpl(jl)
1004 zuu(jl,jk) = puen(jl,ikb-1)
1005 zvu(jl,jk) = pven(jl,ikb-1)
1006 else if ( jk == kcbot(jl) .and. ktype(jl) == 3 ) then
1007 zuu(jl,jk) = puen(jl,jk-1)
1008 zvu(jl,jk) = pven(jl,jk-1)
1009 end if
1010 if ( jk < kcbot(jl) .and. jk >= kctop(jl) ) then
1011 if(momtrans .eq. 1)then
1012 zfac = 0.
1013 if ( ktype(jl) == 1 .or. ktype(jl) == 3 ) zfac = 2.
1014 if ( ktype(jl) == 1 .and. jk <= kctop(jl)+2 ) zfac = 3.
1015 zerate = pmfu(jl,jk) - pmfu(jl,ik) + &
1016 (1.+zfac)*pmfude_rate(jl,jk)
1017 zderate = (1.+zfac)*pmfude_rate(jl,jk)
1018 zmfa = 1./max(cmfcmin,pmfu(jl,jk))
1019 zuu(jl,jk) = (zuu(jl,ik)*pmfu(jl,ik) + &
1020 zerate*puen(jl,jk)-zderate*zuu(jl,ik))*zmfa
1021 zvu(jl,jk) = (zvu(jl,ik)*pmfu(jl,ik) + &
1022 zerate*pven(jl,jk)-zderate*zvu(jl,ik))*zmfa
1023 else
1024 if(ktype(jl) == 1 .or. ktype(jl) == 3) then
1025 pgf_u = -0.7*0.5*(pmfu(jl,ik)*(puen(jl,ik)-puen(jl,jk))+&
1026 pmfu(jl,jk)*(puen(jl,jk)-puen(jl,jk-1)))
1027 pgf_v = -0.7*0.5*(pmfu(jl,ik)*(pven(jl,ik)-pven(jl,jk))+&
1028 pmfu(jl,jk)*(pven(jl,jk)-pven(jl,jk-1)))
1029 else
1030 pgf_u = 0.
1031 pgf_v = 0.
1032 end if
1033 zerate = pmfu(jl,jk) - pmfu(jl,ik) + pmfude_rate(jl,jk)
1034 zderate = pmfude_rate(jl,jk)
1035 zmfa = 1./max(cmfcmin,pmfu(jl,jk))
1036 zuu(jl,jk) = (zuu(jl,ik)*pmfu(jl,ik) + &
1037 zerate*puen(jl,jk)-zderate*zuu(jl,ik)+pgf_u)*zmfa
1038 zvu(jl,jk) = (zvu(jl,ik)*pmfu(jl,ik) + &
1039 zerate*pven(jl,jk)-zderate*zvu(jl,ik)+pgf_v)*zmfa
1040 end if
1041 end if
1042 end if
1043 end do
1044 end do
1045
1046 if(lmfdd) then
1047 do jk = 3 , klev
1048 ik = jk - 1
1049 do jl = 1,klon
1050 if ( ldcum(jl) ) then
1051 if ( jk == idtop(jl) ) then
1052 zud(jl,jk) = 0.5*(zuu(jl,jk)+puen(jl,ik))
1053 zvd(jl,jk) = 0.5*(zvu(jl,jk)+pven(jl,ik))
1054 else if ( jk > idtop(jl) ) then
1055 zerate = -pmfd(jl,jk) + pmfd(jl,ik) + pmfdde_rate(jl,jk)
1056 zmfa = 1./min(-cmfcmin,pmfd(jl,jk))
1057 zud(jl,jk) = (zud(jl,ik)*pmfd(jl,ik) - &
1058 zerate*puen(jl,ik)+pmfdde_rate(jl,jk)*zud(jl,ik))*zmfa
1059 zvd(jl,jk) = (zvd(jl,ik)*pmfd(jl,ik) - &
1060 zerate*pven(jl,ik)+pmfdde_rate(jl,jk)*zvd(jl,ik))*zmfa
1061 end if
1062 end if
1063 end do
1064 end do
1065 end if
1066! --------------------------------------------------
1067! rescale massfluxes for stability in Momentum
1068!------------------------------------------------------------------------
1069 zmfs(:) = 1.
1070 do jk = 2 , klev
1071 do jl = 1, klon
1072 if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
1073 zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons
1074 if ( pmfu(jl,jk) > zmfmax .and. jk >= kctop(jl) ) then
1075 zmfs(jl) = min(zmfs(jl),zmfmax/pmfu(jl,jk))
1076 end if
1077 end if
1078 end do
1079 end do
1080 do jk = 1 , klev
1081 do jl = 1, klon
1082 zmfuus(jl,jk) = pmfu(jl,jk)
1083 zmfdus(jl,jk) = pmfd(jl,jk)
1084 if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
1085 zmfuus(jl,jk) = pmfu(jl,jk)*zmfs(jl)
1086 zmfdus(jl,jk) = pmfd(jl,jk)*zmfs(jl)
1087 end if
1088 end do
1089 end do
1090!* 9.1 update u and v in subroutine cududvn
1091!-------------------------------------------------------------------
1092 do jk = 1 , klev
1093 do jl = 1, klon
1094 ztenu(jl,jk) = pvom(jl,jk)
1095 ztenv(jl,jk) = pvol(jl,jk)
1096 end do
1097 end do
1098
1099 call cududvn(klon,klev,itopm2,ktype,kcbot,kctop, &
1100 ldcum,ztmst,paph,puen,pven,zmfuus,zmfdus,zuu, &
1101 zud,zvu,zvd,pvom,pvol)
1102
1103! calculate KE dissipation
1104 do jl = 1, klon
1105 zsum12(jl) = 0.
1106 zsum22(jl) = 0.
1107 end do
1108 do jk = 1 , klev
1109 do jl = 1, klon
1110 zuv2(jl,jk) = 0.
1111 if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
1112 zdz = (paph(jl,jk+1)-paph(jl,jk))
1113 zduten = pvom(jl,jk) - ztenu(jl,jk)
1114 zdvten = pvol(jl,jk) - ztenv(jl,jk)
1115 zuv2(jl,jk) = sqrt(zduten**2+zdvten**2)
1116 zsum22(jl) = zsum22(jl) + zuv2(jl,jk)*zdz
1117 zsum12(jl) = zsum12(jl) - &
1118 (puen(jl,jk)*zduten+pven(jl,jk)*zdvten)*zdz
1119 end if
1120 end do
1121 end do
1122 do jk = 1 , klev
1123 do jl = 1, klon
1124 if ( ldcum(jl) .and. jk>=kctop(jl)-1 ) then
1125 ztdis = rcpd*zsum12(jl)*zuv2(jl,jk)/max(1.e-15,zsum22(jl))
1126 ptte(jl,jk) = ptte(jl,jk) + ztdis
1127 end if
1128 end do
1129 end do
1130
1131 end if
1132
1133!----------------------------------------------------------------------
1134!* 10. IN CASE THAT EITHER DEEP OR SHALLOW IS SWITCHED OFF
1135! NEED TO SET SOME VARIABLES A POSTERIORI TO ZERO
1136! ---------------------------------------------------
1137 if ( .not. lmfscv .or. .not. lmfpen ) then
1138 do jk = 2 , klev
1139 do jl = 1, klon
1140 if ( llo2(jl) .and. jk >= kctop(jl)-1 ) then
1141 ptu(jl,jk) = pten(jl,jk)
1142 pqu(jl,jk) = pqen(jl,jk)
1143 plu(jl,jk) = 0.
1144 pmfude_rate(jl,jk) = 0.
1145 pmfdde_rate(jl,jk) = 0.
1146 end if
1147 end do
1148 end do
1149 do jl = 1, klon
1150 if ( llo2(jl) ) then
1151 kctop(jl) = klev - 1
1152 kcbot(jl) = klev - 1
1153 end if
1154 end do
1155 end if
1156
1157 !----------------------------------------------------------------------
1158 !* 11.0 CHEMICAL TRACER TRANSPORT
1159 ! -------------------------
1160
1161 if ( ktrac > 0 ) then
1162 ! transport switched off for mid-level convection
1163 do jl = 1, klon
1164 if ( ldcum(jl) .and. ktype(jl) /= 3 .and. &
1165 kcbot(jl)-kctop(jl) >= 1 ) then
1166 lldcum(jl) = .true.
1167 llddraf3(jl) = loddraf(jl)
1168 else
1169 lldcum(jl) = .false.
1170 llddraf3(jl) = .false.
1171 end if
1172 end do
1173 ! check and correct mass fluxes for CFL criterium
1174 zmfs(:) = 1.
1175 do jk = 2 , klev
1176 do jl = 1, klon
1177 if ( lldcum(jl) .and. jk >= kctop(jl) ) then
1178 zmfmax = (paph(jl,jk)-paph(jl,jk-1))*0.8*zcons
1179 if ( pmfu(jl,jk) > zmfmax ) then
1180 zmfs(jl) = min(zmfs(jl),zmfmax/pmfu(jl,jk))
1181 end if
1182 end if
1183 end do
1184 end do
1185
1186 do jk = 1, klev
1187 do jl = 1, klon
1188 if ( lldcum(jl) .and. jk >= kctop(jl)-1 ) then
1189 zmfuus(jl,jk) = pmfu(jl,jk)*zmfs(jl)
1190 zmfudr(jl,jk) = pmfude_rate(jl,jk)*zmfs(jl)
1191 else
1192 zmfuus(jl,jk) = 0.
1193 zmfudr(jl,jk) = 0.
1194 end if
1195 if ( llddraf3(jl) .and. jk >= idtop(jl)-1 ) then
1196 zmfdus(jl,jk) = pmfd(jl,jk)*zmfs(jl)
1197 zmfddr(jl,jk) = pmfdde_rate(jl,jk)*zmfs(jl)
1198 else
1199 zmfdus(jl,jk) = 0.
1200 zmfddr(jl,jk) = 0.
1201 end if
1202 end do
1203 end do
1204
1205 call cuctracer(klon,klev,ktrac,kctop,idtop, &
1206 lldcum,llddraf3,ztmst,paph,zmfuus,zmfdus, &
1207 zmfudr,zmfddr,pcen,ptenc)
1208 end if
1209
1210 return
1211 end subroutine cumastrn
1212
1213!**********************************************
1214! level 3 subroutine cuinin
1215!**********************************************
1216!
1217 subroutine cuinin &
1218 & (klon, klev, klevp1, klevm1, pten,&
1219 & pqen, pqsen, puen, pven, pverv,&
1220 & pgeo, paph, pgeoh, ptenh, pqenh,&
1221 & pqsenh, klwmin, ptu, pqu, ptd,&
1222 & pqd, puu, pvu, pud, pvd,&
1223 & pmfu, pmfd, pmfus, pmfds, pmfuq,&
1224 & pmfdq, pdmfup, pdmfdp, pdpmel, plu,&
1225 & plude, klab)
1226 implicit none
1227! m.tiedtke e.c.m.w.f. 12/89
1228!***purpose
1229! -------
1230! this routine interpolates large-scale fields of t,q etc.
1231! to half levels (i.e. grid for massflux scheme),
1232! and initializes values for updrafts and downdrafts
1233!***interface
1234! ---------
1235! this routine is called from *cumastr*.
1236!***method.
1237! --------
1238! for extrapolation to half levels see tiedtke(1989)
1239!***externals
1240! ---------
1241! *cuadjtq* to specify qs at half levels
1242! ----------------------------------------------------------------
1243 integer klon,klev,klevp1,klevm1
1244 real(kind=kind_phys) pten(klon,klev), pqen(klon,klev),&
1245 & puen(klon,klev), pven(klon,klev),&
1246 & pqsen(klon,klev), pverv(klon,klev),&
1247 & pgeo(klon,klev), pgeoh(klon,klevp1),&
1248 & paph(klon,klevp1), ptenh(klon,klev),&
1249 & pqenh(klon,klev), pqsenh(klon,klev)
1250 real(kind=kind_phys) ptu(klon,klev), pqu(klon,klev),&
1251 & ptd(klon,klev), pqd(klon,klev),&
1252 & puu(klon,klev), pud(klon,klev),&
1253 & pvu(klon,klev), pvd(klon,klev),&
1254 & pmfu(klon,klev), pmfd(klon,klev),&
1255 & pmfus(klon,klev), pmfds(klon,klev),&
1256 & pmfuq(klon,klev), pmfdq(klon,klev),&
1257 & pdmfup(klon,klev), pdmfdp(klon,klev),&
1258 & plu(klon,klev), plude(klon,klev)
1259 real(kind=kind_phys) zwmax(klon), zph(klon), &
1260 & pdpmel(klon,klev)
1261 integer klab(klon,klev), klwmin(klon)
1262 logical loflag(klon)
1263! local variables
1264 integer jl,jk
1265 integer icall,ik
1266 real(kind=kind_phys) zzs
1267!------------------------------------------------------------
1268!* 1. specify large scale parameters at half levels
1269!* adjust temperature fields if staticly unstable
1270!* find level of maximum vertical velocity
1271! -----------------------------------------------------------
1272 do jk=2,klev
1273 do jl=1,klon
1274 ptenh(jl,jk)=(max(cpd*pten(jl,jk-1)+pgeo(jl,jk-1), &
1275 & cpd*pten(jl,jk)+pgeo(jl,jk))-pgeoh(jl,jk))*rcpd
1276 pqenh(jl,jk) = pqen(jl,jk-1)
1277 pqsenh(jl,jk)= pqsen(jl,jk-1)
1278 zph(jl)=paph(jl,jk)
1279 loflag(jl)=.true.
1280 end do
1281
1282 if ( jk >= klev-1 .or. jk < 2 ) cycle
1283 ik=jk
1284 icall=0
1285 call cuadjtqn(klon,klev,ik,zph,ptenh,pqsenh,loflag,icall)
1286 do jl=1,klon
1287 pqenh(jl,jk)=min(pqen(jl,jk-1),pqsen(jl,jk-1)) &
1288 & +(pqsenh(jl,jk)-pqsen(jl,jk-1))
1289 pqenh(jl,jk)=max(pqenh(jl,jk),0.)
1290 end do
1291 end do
1292
1293 do jl=1,klon
1294 ptenh(jl,klev)=(cpd*pten(jl,klev)+pgeo(jl,klev)- &
1295 & pgeoh(jl,klev))*rcpd
1296 pqenh(jl,klev)=pqen(jl,klev)
1297 ptenh(jl,1)=pten(jl,1)
1298 pqenh(jl,1)=pqen(jl,1)
1299 klwmin(jl)=klev
1300 zwmax(jl)=0.
1301 end do
1302
1303 do jk=klevm1,2,-1
1304 do jl=1,klon
1305 zzs=max(cpd*ptenh(jl,jk)+pgeoh(jl,jk), &
1306 & cpd*ptenh(jl,jk+1)+pgeoh(jl,jk+1))
1307 ptenh(jl,jk)=(zzs-pgeoh(jl,jk))*rcpd
1308 end do
1309 end do
1310
1311 do jk=klev,3,-1
1312 do jl=1,klon
1313 if(pverv(jl,jk).lt.zwmax(jl)) then
1314 zwmax(jl)=pverv(jl,jk)
1315 klwmin(jl)=jk
1316 end if
1317 end do
1318 end do
1319!-----------------------------------------------------------
1320!* 2.0 initialize values for updrafts and downdrafts
1321!-----------------------------------------------------------
1322 do jk=1,klev
1323 ik=jk-1
1324 if(jk.eq.1) ik=1
1325 do jl=1,klon
1326 ptu(jl,jk)=ptenh(jl,jk)
1327 ptd(jl,jk)=ptenh(jl,jk)
1328 pqu(jl,jk)=pqenh(jl,jk)
1329 pqd(jl,jk)=pqenh(jl,jk)
1330 plu(jl,jk)=0.
1331 puu(jl,jk)=puen(jl,ik)
1332 pud(jl,jk)=puen(jl,ik)
1333 pvu(jl,jk)=pven(jl,ik)
1334 pvd(jl,jk)=pven(jl,ik)
1335 klab(jl,jk)=0
1336 end do
1337 end do
1338 return
1339 end subroutine cuinin
1340
1341!---------------------------------------------------------
1342! level 3 souroutines
1343!--------------------------------------------------------
1344 subroutine cutypen &
1345 & ( klon, klev, klevp1, klevm1, pqen,&
1346 & ptenh, pqenh, pqsenh, pgeoh, paph,&
1347 & hfx, qfx, pgeo, pqsen, pap,&
1348 & pten, lndj, cutu, cuqu, culab,&
1349 & ldcum, cubot, cutop, ktype, wbase, culu, kdpl )
1350! zhang & wang iprc 2011-2013
1351!***purpose.
1352! --------
1353! to produce first guess updraught for cu-parameterizations
1354! calculates condensation level, and sets updraught base variables and
1355! first guess cloud type
1356!***interface
1357! ---------
1358! this routine is called from *cumastr*.
1359! input are environm. values of t,q,p,phi at half levels.
1360! it returns cloud types as follows;
1361! ktype=1 for deep cumulus
1362! ktype=2 for shallow cumulus
1363!***method.
1364! --------
1365! based on a simplified updraught equation
1366! partial(hup)/partial(z)=eta(h - hup)
1367! eta is the entrainment rate for test parcel
1368! h stands for dry static energy or the total water specific humidity
1369! references: christian jakob, 2003: a new subcloud model for
1370! mass-flux convection schemes
1371! influence on triggering, updraft properties, and model
1372! climate, mon.wea.rev.
1373! 131, 2765-2778
1374! and
1375! ifs documentation - cy36r1,cy38r1
1376!***input variables:
1377! ptenh [ztenh] - environment temperature on half levels. (cuini)
1378! pqenh [zqenh] - env. specific humidity on half levels. (cuini)
1379! pgeoh [zgeoh] - geopotential on half levels, (mssflx)
1380! paph - pressure of half levels. (mssflx)
1381! rho - density of the lowest model level
1382! qfx - net upward moisture flux at the surface (kg/m^2/s)
1383! hfx - net upward heat flux at the surface (w/m^2)
1384!***variables output by cutype:
1385! ktype - convection type - 1: penetrative (cumastr)
1386! 2: stratocumulus (cumastr)
1387! 3: mid-level (cuasc)
1388! information for updraft parcel (ptu,pqu,plu,kcbot,klab,kdpl...)
1389! ----------------------------------------------------------------
1390!-------------------------------------------------------------------
1391 implicit none
1392!-------------------------------------------------------------------
1393 integer klon, klev, klevp1, klevm1
1394 real(kind=kind_phys) ptenh(klon,klev), pqenh(klon,klev),&
1395 & pqsen(klon,klev), pqsenh(klon,klev),&
1396 & pgeoh(klon,klevp1), paph(klon,klevp1),&
1397 & pap(klon,klev), pqen(klon,klev)
1398 real(kind=kind_phys) pten(klon,klev)
1399 real(kind=kind_phys) ptu(klon,klev),pqu(klon,klev),plu(klon,klev)
1400 real(kind=kind_phys) pgeo(klon,klev)
1401 integer klab(klon,klev)
1402 integer kctop(klon),kcbot(klon)
1403
1404 real(kind=kind_phys) qfx(klon),hfx(klon)
1405 real(kind=kind_phys) zph(klon)
1406 integer lndj(klon)
1407 logical loflag(klon), deepflag(klon), resetflag(klon)
1408
1409! output variables
1410 real(kind=kind_phys) cutu(klon,klev), cuqu(klon,klev), culu(klon,klev)
1411 integer culab(klon,klev)
1412 real(kind=kind_phys) wbase(klon)
1413 integer ktype(klon),cubot(klon),cutop(klon),kdpl(klon)
1414 logical ldcum(klon)
1415
1416! local variables
1417 real(kind=kind_phys) zqold(klon)
1418 real(kind=kind_phys) rho, part1, part2, root, conw, deltt, deltq
1419 real(kind=kind_phys) eta(klon),dz(klon),coef(klon)
1420 real(kind=kind_phys) dhen(klon,klev), dh(klon,klev)
1421 real(kind=kind_phys) plude(klon,klev)
1422 real(kind=kind_phys) kup(klon,klev)
1423 real(kind=kind_phys) vptu(klon,klev),vten(klon,klev)
1424 real(kind=kind_phys) zbuo(klon,klev),abuoy(klon,klev)
1425
1426 real(kind=kind_phys) zz,zdken,zdq
1427 real(kind=kind_phys) fscale,crirh1,pp
1428 real(kind=kind_phys) atop1,atop2,abot
1429 real(kind=kind_phys) tmix,zmix,qmix,pmix
1430 real(kind=kind_phys) zlglac,dp
1431 integer nk,is,ikb,ikt
1432
1433 real(kind=kind_phys) zqsu,zcor,zdp,zesdp,zalfaw,zfacw,zfaci,zfac,zdsdp,zdqsdt,zdtdp
1434 real(kind=kind_phys) zpdifftop, zpdiffbot
1435 integer zcbase(klon), itoppacel(klon)
1436 integer jl,jk,ik,icall,levels
1437 logical needreset, lldcum(klon)
1438!--------------------------------------------------------------
1439 do jl=1,klon
1440 kcbot(jl)=klev
1441 kctop(jl)=klev
1442 kdpl(jl) =klev
1443 ktype(jl)=0
1444 wbase(jl)=0.
1445 ldcum(jl)=.false.
1446 end do
1447
1448!-----------------------------------------------------------
1449! let's do test,and check the shallow convection first
1450! the first level is klev
1451! define deltat and deltaq
1452!-----------------------------------------------------------
1453 do jk=1,klev
1454 do jl=1,klon
1455 plu(jl,jk)=culu(jl,jk) ! parcel liquid water
1456 ptu(jl,jk)=cutu(jl,jk) ! parcel temperature
1457 pqu(jl,jk)=cuqu(jl,jk) ! parcel specific humidity
1458 klab(jl,jk)=culab(jl,jk)
1459 dh(jl,jk)=0.0 ! parcel dry static energy
1460 dhen(jl,jk)=0.0 ! environment dry static energy
1461 kup(jl,jk)=0.0 ! updraught kinetic energy for parcel
1462 vptu(jl,jk)=0.0 ! parcel virtual temperature considering water-loading
1463 vten(jl,jk)=0.0 ! environment virtual temperature
1464 zbuo(jl,jk)=0.0 ! parcel buoyancy
1465 abuoy(jl,jk)=0.0
1466 end do
1467 end do
1468
1469 do jl=1,klon
1470 zqold(jl) = 0.
1471 lldcum(jl) = .false.
1472 loflag(jl) = .true.
1473 end do
1474
1475! check the levels from lowest level to second top level
1476 do jk=klevm1,2,-1
1477
1478! define the variables at the first level
1479 if(jk .eq. klevm1) then
1480 do jl=1,klon
1481 rho=pap(jl,klev)/ &
1482 & (rd*(pten(jl,klev)*(1.+vtmpc1*pqen(jl,klev))))
1483 hfx(jl) = hfx(jl)*rho*cpd
1484 qfx(jl) = qfx(jl)*rho
1485 part1 = 1.5*0.4*pgeo(jl,klev)/ &
1486 & (rho*pten(jl,klev))
1487 part2 = -hfx(jl)*rcpd-vtmpc1*pten(jl,klev)*qfx(jl)
1488 root = 0.001-part1*part2
1489 if(part2 .lt. 0.) then
1490 conw = 1.2*(root)**t13
1491 deltt = max(1.5*hfx(jl)/(rho*cpd*conw),0.)
1492 deltq = max(1.5*qfx(jl)/(rho*conw),0.)
1493 kup(jl,klev) = 0.5*(conw**2)
1494 pqu(jl,klev)= pqenh(jl,klev) + deltq
1495 dhen(jl,klev)= pgeoh(jl,klev) + ptenh(jl,klev)*cpd
1496 dh(jl,klev) = dhen(jl,klev) + deltt*cpd
1497 ptu(jl,klev) = (dh(jl,klev)-pgeoh(jl,klev))*rcpd
1498 vptu(jl,klev)=ptu(jl,klev)*(1.+vtmpc1*pqu(jl,klev))
1499 vten(jl,klev)=ptenh(jl,klev)*(1.+vtmpc1*pqenh(jl,klev))
1500 zbuo(jl,klev)=(vptu(jl,klev)-vten(jl,klev))/vten(jl,klev)
1501 klab(jl,klev) = 1
1502 else
1503 loflag(jl) = .false.
1504 end if
1505 end do
1506 end if
1507
1508 is=0
1509 do jl=1,klon
1510 if(loflag(jl))then
1511 is=is+1
1512 endif
1513 enddo
1514 if(is.eq.0) exit
1515
1516! the next levels, we use the variables at the first level as initial values
1517 do jl=1,klon
1518 if(loflag(jl)) then
1519 eta(jl) = 0.55/(pgeo(jl,jk)*zrg)+1.e-4
1520 dz(jl) = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
1521 coef(jl)= 0.5*eta(jl)*dz(jl)
1522 dhen(jl,jk) = pgeoh(jl,jk) + cpd*ptenh(jl,jk)
1523 dh(jl,jk) = (coef(jl)*(dhen(jl,jk+1)+dhen(jl,jk))&
1524 & +(1.-coef(jl))*dh(jl,jk+1))/(1.+coef(jl))
1525 pqu(jl,jk) =(coef(jl)*(pqenh(jl,jk+1)+pqenh(jl,jk))&
1526 & +(1.-coef(jl))*pqu(jl,jk+1))/(1.+coef(jl))
1527 ptu(jl,jk) = (dh(jl,jk)-pgeoh(jl,jk))*rcpd
1528 zqold(jl) = pqu(jl,jk)
1529 zph(jl)=paph(jl,jk)
1530 end if
1531 end do
1532! check if the parcel is saturated
1533 ik=jk
1534 icall=1
1535 call cuadjtqn(klon,klev,ik,zph,ptu,pqu,loflag,icall)
1536 do jl=1,klon
1537 if( loflag(jl) ) then
1538 zdq = max((zqold(jl) - pqu(jl,jk)),0.)
1539 plu(jl,jk) = plu(jl,jk+1) + zdq
1540 zlglac=zdq*((1.-foealfa(ptu(jl,jk))) - &
1541 (1.-foealfa(ptu(jl,jk+1))))
1542 plu(jl,jk) = min(plu(jl,jk),5.e-3)
1543 dh(jl,jk) = pgeoh(jl,jk) + cpd*(ptu(jl,jk)+ralfdcp*zlglac)
1544! compute the updraft speed
1545 vptu(jl,jk) = ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk))+&
1546 ralfdcp*zlglac
1547 vten(jl,jk) = ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
1548 zbuo(jl,jk) = (vptu(jl,jk) - vten(jl,jk))/vten(jl,jk)
1549 abuoy(jl,jk)=(zbuo(jl,jk)+zbuo(jl,jk+1))*0.5*g
1550 atop1 = 1.0 - 2.*coef(jl)
1551 atop2 = 2.0*dz(jl)*abuoy(jl,jk)
1552 abot = 1.0 + 2.*coef(jl)
1553 kup(jl,jk) = (atop1*kup(jl,jk+1) + atop2) / abot
1554
1555! let's find the exact cloud base
1556 if ( plu(jl,jk) > 0. .and. klab(jl,jk+1) == 1 ) then
1557 ik = jk + 1
1558 zqsu = foeewm(ptu(jl,ik))/paph(jl,ik)
1559 zqsu = min(0.5,zqsu)
1560 zcor = 1./(1.-vtmpc1*zqsu)
1561 zqsu = zqsu*zcor
1562 zdq = min(0.,pqu(jl,ik)-zqsu)
1563 zalfaw = foealfa(ptu(jl,ik))
1564 zfacw = c5les/((ptu(jl,ik)-c4les)**2)
1565 zfaci = c5ies/((ptu(jl,ik)-c4ies)**2)
1566 zfac = zalfaw*zfacw + (1.-zalfaw)*zfaci
1567 zesdp = foeewm(ptu(jl,ik))/paph(jl,ik)
1568 zcor = 1./(1.-vtmpc1*zesdp)
1569 zdqsdt = zfac*zcor*zqsu
1570 zdtdp = rd*ptu(jl,ik)/(cpd*paph(jl,ik))
1571 zdp = zdq/(zdqsdt*zdtdp)
1572 zcbase(jl) = paph(jl,ik) + zdp
1573! chose nearest half level as cloud base (jk or jk+1)
1574 zpdifftop = zcbase(jl) - paph(jl,jk)
1575 zpdiffbot = paph(jl,jk+1) - zcbase(jl)
1576 if ( zpdifftop > zpdiffbot .and. kup(jl,jk+1) > 0. ) then
1577 ikb = min(klev-1,jk+1)
1578 klab(jl,ikb) = 2
1579 klab(jl,jk) = 2
1580 kcbot(jl) = ikb
1581 plu(jl,jk+1) = 1.0e-8
1582 else if ( zpdifftop <= zpdiffbot .and.kup(jl,jk) > 0. ) then
1583 klab(jl,jk) = 2
1584 kcbot(jl) = jk
1585 end if
1586 end if
1587
1588 if(kup(jl,jk) .lt. 0.)then
1589 loflag(jl) = .false.
1590 if(plu(jl,jk+1) .gt. 0.) then
1591 kctop(jl) = jk
1592 lldcum(jl) = .true.
1593 else
1594 lldcum(jl) = .false.
1595 end if
1596 else
1597 if(plu(jl,jk) .gt. 0.)then
1598 klab(jl,jk)=2
1599 else
1600 klab(jl,jk)=1
1601 end if
1602 end if
1603 end if
1604 end do
1605
1606 end do ! end all the levels
1607
1608 do jl=1,klon
1609 ikb = kcbot(jl)
1610 ikt = kctop(jl)
1611 if(paph(jl,ikb) - paph(jl,ikt) > zdnoprc) lldcum(jl) = .false.
1612 if(lldcum(jl)) then
1613 ktype(jl) = 2
1614 ldcum(jl) = .true.
1615 wbase(jl) = sqrt(max(2.*kup(jl,ikb),0.))
1616 cubot(jl) = ikb
1617 cutop(jl) = ikt
1618 kdpl(jl) = klev
1619 else
1620 cutop(jl) = -1
1621 cubot(jl) = -1
1622 kdpl(jl) = klev - 1
1623 ldcum(jl) = .false.
1624 wbase(jl) = 0.
1625 end if
1626 end do
1627
1628 do jk=klev,1,-1
1629 do jl=1,klon
1630 ikt = kctop(jl)
1631 if(jk .ge. ikt)then
1632 culab(jl,jk) = klab(jl,jk)
1633 cutu(jl,jk) = ptu(jl,jk)
1634 cuqu(jl,jk) = pqu(jl,jk)
1635 culu(jl,jk) = plu(jl,jk)
1636 end if
1637 end do
1638 end do
1639
1640!-----------------------------------------------------------
1641! next, let's check the deep convection
1642! the first level is klevm1-1
1643! define deltat and deltaq
1644!----------------------------------------------------------
1645! we check the parcel starting level by level
1646! assume the mix-layer is 60hPa
1647 deltt = 0.2
1648 deltq = 1.0e-4
1649 do jl=1,klon
1650 deepflag(jl) = .false.
1651 end do
1652
1653 do jk=klev,1,-1
1654 do jl=1,klon
1655 if((paph(jl,klev+1)-paph(jl,jk)) .lt. 350.e2) itoppacel(jl) = jk
1656 end do
1657 end do
1658
1659 do levels=klevm1-1,klevm1-20,-1 ! loop starts
1660 do jk=1,klev
1661 do jl=1,klon
1662 plu(jl,jk)=0.0 ! parcel liquid water
1663 ptu(jl,jk)=0.0 ! parcel temperature
1664 pqu(jl,jk)=0.0 ! parcel specific humidity
1665 dh(jl,jk)=0.0 ! parcel dry static energy
1666 dhen(jl,jk)=0.0 ! environment dry static energy
1667 kup(jl,jk)=0.0 ! updraught kinetic energy for parcel
1668 vptu(jl,jk)=0.0 ! parcel virtual temperature consideringwater-loading
1669 vten(jl,jk)=0.0 ! environment virtual temperature
1670 abuoy(jl,jk)=0.0
1671 zbuo(jl,jk)=0.0
1672 klab(jl,jk)=0
1673 end do
1674 end do
1675
1676 do jl=1,klon
1677 kcbot(jl) = levels
1678 kctop(jl) = levels
1679 zqold(jl) = 0.
1680 lldcum(jl) = .false.
1681 resetflag(jl)= .false.
1682 loflag(jl) = (.not. deepflag(jl)) .and. (levels.ge.itoppacel(jl))
1683 end do
1684
1685! start the inner loop to search the deep convection points
1686 do jk=levels,2,-1
1687 is=0
1688 do jl=1,klon
1689 if(loflag(jl))then
1690 is=is+1
1691 endif
1692 enddo
1693 if(is.eq.0) exit
1694
1695! define the variables at the departure level
1696 if(jk .eq. levels) then
1697 do jl=1,klon
1698 if(loflag(jl)) then
1699 if((paph(jl,klev+1)-paph(jl,jk)) < 60.e2) then
1700 tmix=0.
1701 qmix=0.
1702 zmix=0.
1703 pmix=0.
1704 do nk=jk+2,jk,-1
1705 if(pmix < 50.e2) then
1706 dp = paph(jl,nk) - paph(jl,nk-1)
1707 tmix=tmix+dp*ptenh(jl,nk)
1708 qmix=qmix+dp*pqenh(jl,nk)
1709 zmix=zmix+dp*pgeoh(jl,nk)
1710 pmix=pmix+dp
1711 end if
1712 end do
1713 tmix=tmix/pmix
1714 qmix=qmix/pmix
1715 zmix=zmix/pmix
1716 else
1717 tmix=ptenh(jl,jk+1)
1718 qmix=pqenh(jl,jk+1)
1719 zmix=pgeoh(jl,jk+1)
1720 end if
1721
1722 pqu(jl,jk+1) = qmix + deltq
1723 dhen(jl,jk+1)= zmix + tmix*cpd
1724 dh(jl,jk+1) = dhen(jl,jk+1) + deltt*cpd
1725 ptu(jl,jk+1) = (dh(jl,jk+1)-pgeoh(jl,jk+1))*rcpd
1726 kup(jl,jk+1) = 0.5
1727 klab(jl,jk+1)= 1
1728 vptu(jl,jk+1)=ptu(jl,jk+1)*(1.+vtmpc1*pqu(jl,jk+1))
1729 vten(jl,jk+1)=ptenh(jl,jk+1)*(1.+vtmpc1*pqenh(jl,jk+1))
1730 zbuo(jl,jk+1)=(vptu(jl,jk+1)-vten(jl,jk+1))/vten(jl,jk+1)
1731 end if
1732 end do
1733 end if
1734
1735! the next levels, we use the variables at the first level as initial values
1736 do jl=1,klon
1737 if(loflag(jl)) then
1738! define the fscale
1739 fscale = min(1.,(pqsen(jl,jk)/pqsen(jl,levels))**3)
1740 eta(jl) = 1.75e-3*fscale
1741 dz(jl) = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
1742 coef(jl)= 0.5*eta(jl)*dz(jl)
1743 dhen(jl,jk) = pgeoh(jl,jk) + cpd*ptenh(jl,jk)
1744 dh(jl,jk) = (coef(jl)*(dhen(jl,jk+1)+dhen(jl,jk))&
1745 & +(1.-coef(jl))*dh(jl,jk+1))/(1.+coef(jl))
1746 pqu(jl,jk) =(coef(jl)*(pqenh(jl,jk+1)+pqenh(jl,jk))&
1747 & +(1.-coef(jl))*pqu(jl,jk+1))/(1.+coef(jl))
1748 ptu(jl,jk) = (dh(jl,jk)-pgeoh(jl,jk))*rcpd
1749 zqold(jl) = pqu(jl,jk)
1750 zph(jl)=paph(jl,jk)
1751 end if
1752 end do
1753! check if the parcel is saturated
1754 ik=jk
1755 icall=1
1756 call cuadjtqn(klon,klev,ik,zph,ptu,pqu,loflag,icall)
1757
1758 do jl=1,klon
1759 if( loflag(jl) ) then
1760 zdq = max((zqold(jl) - pqu(jl,jk)),0.)
1761 plu(jl,jk) = plu(jl,jk+1) + zdq
1762 zlglac=zdq*((1.-foealfa(ptu(jl,jk))) - &
1763 (1.-foealfa(ptu(jl,jk+1))))
1764 plu(jl,jk) = 0.5*plu(jl,jk)
1765 dh(jl,jk) = pgeoh(jl,jk) + cpd*(ptu(jl,jk)+ralfdcp*zlglac)
1766! compute the updraft speed
1767 vptu(jl,jk) = ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk))+&
1768 ralfdcp*zlglac
1769 vten(jl,jk) = ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
1770 zbuo(jl,jk) = (vptu(jl,jk) - vten(jl,jk))/vten(jl,jk)
1771 abuoy(jl,jk)=(zbuo(jl,jk)+zbuo(jl,jk+1))*0.5*g
1772 atop1 = 1.0 - 2.*coef(jl)
1773 atop2 = 2.0*dz(jl)*abuoy(jl,jk)
1774 abot = 1.0 + 2.*coef(jl)
1775 kup(jl,jk) = (atop1*kup(jl,jk+1) + atop2) / abot
1776! let's find the exact cloud base
1777 if ( plu(jl,jk) > 0. .and. klab(jl,jk+1) == 1 ) then
1778 ik = jk + 1
1779 zqsu = foeewm(ptu(jl,ik))/paph(jl,ik)
1780 zqsu = min(0.5,zqsu)
1781 zcor = 1./(1.-vtmpc1*zqsu)
1782 zqsu = zqsu*zcor
1783 zdq = min(0.,pqu(jl,ik)-zqsu)
1784 zalfaw = foealfa(ptu(jl,ik))
1785 zfacw = c5les/((ptu(jl,ik)-c4les)**2)
1786 zfaci = c5ies/((ptu(jl,ik)-c4ies)**2)
1787 zfac = zalfaw*zfacw + (1.-zalfaw)*zfaci
1788 zesdp = foeewm(ptu(jl,ik))/paph(jl,ik)
1789 zcor = 1./(1.-vtmpc1*zesdp)
1790 zdqsdt = zfac*zcor*zqsu
1791 zdtdp = rd*ptu(jl,ik)/(cpd*paph(jl,ik))
1792 zdp = zdq/(zdqsdt*zdtdp)
1793 zcbase(jl) = paph(jl,ik) + zdp
1794! chose nearest half level as cloud base (jk or jk+1)
1795 zpdifftop = zcbase(jl) - paph(jl,jk)
1796 zpdiffbot = paph(jl,jk+1) - zcbase(jl)
1797 if ( zpdifftop > zpdiffbot .and. kup(jl,jk+1) > 0. ) then
1798 ikb = min(klev-1,jk+1)
1799 klab(jl,ikb) = 2
1800 klab(jl,jk) = 2
1801 kcbot(jl) = ikb
1802 plu(jl,jk+1) = 1.0e-8
1803 else if ( zpdifftop <= zpdiffbot .and.kup(jl,jk) > 0. ) then
1804 klab(jl,jk) = 2
1805 kcbot(jl) = jk
1806 end if
1807 end if
1808
1809 if(kup(jl,jk) .lt. 0.)then
1810 loflag(jl) = .false.
1811 if(plu(jl,jk+1) .gt. 0.) then
1812 kctop(jl) = jk
1813 lldcum(jl) = .true.
1814 else
1815 lldcum(jl) = .false.
1816 end if
1817 else
1818 if(plu(jl,jk) .gt. 0.)then
1819 klab(jl,jk)=2
1820 else
1821 klab(jl,jk)=1
1822 end if
1823 end if
1824 end if
1825 end do
1826
1827 end do ! end all the levels
1828
1829 needreset = .false.
1830 do jl=1,klon
1831 ikb = kcbot(jl)
1832 ikt = kctop(jl)
1833 if(paph(jl,ikb) - paph(jl,ikt) < zdnoprc) lldcum(jl) = .false.
1834 if(lldcum(jl)) then
1835 ktype(jl) = 1
1836 ldcum(jl) = .true.
1837 deepflag(jl) = .true.
1838 wbase(jl) = sqrt(max(2.*kup(jl,ikb),0.))
1839 cubot(jl) = ikb
1840 cutop(jl) = ikt
1841 kdpl(jl) = levels+1
1842 needreset = .true.
1843 resetflag(jl)= .true.
1844 end if
1845 end do
1846
1847 if(needreset) then
1848 do jk=klev,1,-1
1849 do jl=1,klon
1850 if(resetflag(jl)) then
1851 ikt = kctop(jl)
1852 ikb = kdpl(jl)
1853 if(jk .le. ikb .and. jk .ge. ikt )then
1854 culab(jl,jk) = klab(jl,jk)
1855 cutu(jl,jk) = ptu(jl,jk)
1856 cuqu(jl,jk) = pqu(jl,jk)
1857 culu(jl,jk) = plu(jl,jk)
1858 else
1859 culab(jl,jk) = 1
1860 cutu(jl,jk) = ptenh(jl,jk)
1861 cuqu(jl,jk) = pqenh(jl,jk)
1862 culu(jl,jk) = 0.
1863 end if
1864 if ( jk .lt. ikt ) culab(jl,jk) = 0
1865 end if
1866 end do
1867 end do
1868 end if
1869
1870 end do ! end all cycles
1871
1872 return
1873 end subroutine cutypen
1874
1875!-----------------------------------------------------------------
1876! level 3 subroutines 'cuascn'
1877!-----------------------------------------------------------------
1878 subroutine cuascn &
1879 & (klon, klev, klevp1, klevm1, ptenh,&
1880 & pqenh, puen, pven, pten, pqen,&
1881 & pqsen, pgeo, pgeoh, pap, paph,&
1882 & pqte, pverv, klwmin, ldcum, phcbase,&
1883 & ktype, klab, ptu, pqu, plu,&
1884 & puu, pvu, pmfu, pmfub, &
1885 & pmfus, pmfuq, pmful, plude, pdmfup,&
1886 & kcbot, kctop, kctop0, kcum, ztmst,&
1887 & pqsenh, plglac, lndj, wup, wbase, kdpl, pmfude_rate)
1888 implicit none
1889! this routine does the calculations for cloud ascents
1890! for cumulus parameterization
1891! m.tiedtke e.c.m.w.f. 7/86 modif. 12/89
1892! y.wang iprc 11/01 modif.
1893! c.zhang iprc 05/12 modif.
1894!***purpose.
1895! --------
1896! to produce cloud ascents for cu-parametrization
1897! (vertical profiles of t,q,l,u and v and corresponding
1898! fluxes as well as precipitation rates)
1899!***interface
1900! ---------
1901! this routine is called from *cumastr*.
1902!***method.
1903! --------
1904! lift surface air dry-adiabatically to cloud base
1905! and then calculate moist ascent for
1906! entraining/detraining plume.
1907! entrainment and detrainment rates differ for
1908! shallow and deep cumulus convection.
1909! in case there is no penetrative or shallow convection
1910! check for possibility of mid level convection
1911! (cloud base values calculated in *cubasmc*)
1912!***externals
1913! ---------
1914! *cuadjtqn* adjust t and q due to condensation in ascent
1915! *cuentrn* calculate entrainment/detrainment rates
1916! *cubasmcn* calculate cloud base values for midlevel convection
1917!***reference
1918! ---------
1919! (tiedtke,1989)
1920!***input variables:
1921! ptenh [ztenh] - environ temperature on half levels. (cuini)
1922! pqenh [zqenh] - env. specific humidity on half levels. (cuini)
1923! puen - environment wind u-component. (mssflx)
1924! pven - environment wind v-component. (mssflx)
1925! pten - environment temperature. (mssflx)
1926! pqen - environment specific humidity. (mssflx)
1927! pqsen - environment saturation specific humidity. (mssflx)
1928! pgeo - geopotential. (mssflx)
1929! pgeoh [zgeoh] - geopotential on half levels, (mssflx)
1930! pap - pressure in pa. (mssflx)
1931! paph - pressure of half levels. (mssflx)
1932! pqte - moisture convergence (delta q/delta t). (mssflx)
1933! pverv - large scale vertical velocity (omega). (mssflx)
1934! klwmin [ilwmin] - level of minimum omega. (cuini)
1935! klab [ilab] - level label - 1: sub-cloud layer.
1936! 2: condensation level (cloud base)
1937! pmfub [zmfub] - updraft mass flux at cloud base. (cumastr)
1938!***variables modified by cuasc:
1939! ldcum - logical denoting profiles. (cubase)
1940! ktype - convection type - 1: penetrative (cumastr)
1941! 2: stratocumulus (cumastr)
1942! 3: mid-level (cuasc)
1943! ptu - cloud temperature.
1944! pqu - cloud specific humidity.
1945! plu - cloud liquid water (moisture condensed out)
1946! puu [zuu] - cloud momentum u-component.
1947! pvu [zvu] - cloud momentum v-component.
1948! pmfu - updraft mass flux.
1949! pmfus [zmfus] - updraft flux of dry static energy. (cubasmc)
1950! pmfuq [zmfuq] - updraft flux of specific humidity.
1951! pmful [zmful] - updraft flux of cloud liquid water.
1952! plude - liquid water returned to environment by detrainment.
1953! pdmfup [zmfup] -
1954! kcbot - cloud base level. (cubase)
1955! kctop - cloud top level
1956! kctop0 [ictop0] - estimate of cloud top. (cumastr)
1957! kcum [icum] - flag to control the call
1958
1959 integer klev,klon,klevp1,klevm1
1960 real(kind=kind_phys) ptenh(klon,klev), pqenh(klon,klev), &
1961 & puen(klon,klev), pven(klon,klev),&
1962 & pten(klon,klev), pqen(klon,klev),&
1963 & pgeo(klon,klev), pgeoh(klon,klevp1),&
1964 & pap(klon,klev), paph(klon,klevp1),&
1965 & pqsen(klon,klev), pqte(klon,klev),&
1966 & pverv(klon,klev), pqsenh(klon,klev)
1967 real(kind=kind_phys) ptu(klon,klev), pqu(klon,klev),&
1968 & puu(klon,klev), pvu(klon,klev),&
1969 & pmfu(klon,klev), zph(klon),&
1970 & pmfub(klon), &
1971 & pmfus(klon,klev), pmfuq(klon,klev),&
1972 & plu(klon,klev), plude(klon,klev),&
1973 & pmful(klon,klev), pdmfup(klon,klev)
1974 real(kind=kind_phys) zdmfen(klon), zdmfde(klon),&
1975 & zmfuu(klon), zmfuv(klon),&
1976 & zpbase(klon), zqold(klon)
1977 real(kind=kind_phys) phcbase(klon), zluold(klon)
1978 real(kind=kind_phys) zprecip(klon), zlrain(klon,klev)
1979 real(kind=kind_phys) zbuo(klon,klev), kup(klon,klev)
1980 real(kind=kind_phys) wup(klon)
1981 real(kind=kind_phys) wbase(klon), zodetr(klon,klev)
1982 real(kind=kind_phys) plglac(klon,klev)
1983
1984 real(kind=kind_phys) eta(klon),dz(klon)
1985
1986 integer klwmin(klon), ktype(klon),&
1987 & klab(klon,klev), kcbot(klon),&
1988 & kctop(klon), kctop0(klon)
1989 integer lndj(klon)
1990 logical ldcum(klon), loflag(klon)
1991 logical llo2,llo3, llo1(klon)
1992
1993 integer kdpl(klon)
1994 real(kind=kind_phys) zoentr(klon), zdpmean(klon)
1995 real(kind=kind_phys) pdmfen(klon,klev), pmfude_rate(klon,klev)
1996! local variables
1997 integer jl,jk
1998 integer ikb,icum,itopm2,ik,icall,is,kcum,jlm,jll
1999 integer jlx(klon)
2000 real(kind=kind_phys) ztmst,zcons2,zfacbuo,zprcdgw,z_cwdrag,z_cldmax,z_cwifrac,z_cprc2
2001 real(kind=kind_phys) zmftest,zmfmax,zqeen,zseen,zscde,zqude
2002 real(kind=kind_phys) zmfusk,zmfuqk,zmfulk
2003 real(kind=kind_phys) zbc,zbe,zkedke,zmfun,zwu,zprcon,zdt,zcbf,zzco
2004 real(kind=kind_phys) zlcrit,zdfi,zc,zd,zint,zlnew,zvw,zvi,zalfaw,zrold
2005 real(kind=kind_phys) zrnew,zz,zdmfeu,zdmfdu,dp
2006 real(kind=kind_phys) zfac,zbuoc,zdkbuo,zdken,zvv,zarg,zchange,zxe,zxs,zdshrd
2007 real(kind=kind_phys) atop1,atop2,abot
2008!--------------------------------
2009!* 1. specify parameters
2010!--------------------------------
2011 zcons2=3./(g*ztmst)
2012 zfacbuo = 0.5/(1.+0.5)
2013 zprcdgw = cprcon*zrg
2014 z_cldmax = 5.e-3
2015 z_cwifrac = 0.5
2016 z_cprc2 = 0.5
2017 z_cwdrag = (3.0/8.0)*0.506/0.2
2018!---------------------------------
2019! 2. set default values
2020!---------------------------------
2021 llo3 = .false.
2022 do jl=1,klon
2023 zluold(jl)=0.
2024 wup(jl)=0.
2025 zdpmean(jl)=0.
2026 zoentr(jl)=0.
2027 if(.not.ldcum(jl)) then
2028 ktype(jl)=0
2029 kcbot(jl) = -1
2030 pmfub(jl) = 0.
2031 pqu(jl,klev) = 0.
2032 end if
2033 end do
2034
2035 ! initialize variout quantities
2036 do jk=1,klev
2037 do jl=1,klon
2038 if(jk.ne.kcbot(jl)) plu(jl,jk)=0.
2039 pmfu(jl,jk)=0.
2040 pmfus(jl,jk)=0.
2041 pmfuq(jl,jk)=0.
2042 pmful(jl,jk)=0.
2043 plude(jl,jk)=0.
2044 plglac(jl,jk)=0.
2045 pdmfup(jl,jk)=0.
2046 zlrain(jl,jk)=0.
2047 zbuo(jl,jk)=0.
2048 kup(jl,jk)=0.
2049 pdmfen(jl,jk) = 0.
2050 pmfude_rate(jl,jk) = 0.
2051 if(.not.ldcum(jl).or.ktype(jl).eq.3) klab(jl,jk)=0
2052 if(.not.ldcum(jl).and.paph(jl,jk).lt.4.e4) kctop0(jl)=jk
2053 end do
2054 end do
2055
2056 do jl = 1,klon
2057 if ( ktype(jl) == 3 ) ldcum(jl) = .false.
2058 end do
2059!------------------------------------------------
2060! 3.0 initialize values at cloud base level
2061!------------------------------------------------
2062 do jl=1,klon
2063 kctop(jl)=kcbot(jl)
2064 if(ldcum(jl)) then
2065 ikb = kcbot(jl)
2066 kup(jl,ikb) = 0.5*wbase(jl)**2
2067 pmfu(jl,ikb) = pmfub(jl)
2068 pmfus(jl,ikb) = pmfub(jl)*(cpd*ptu(jl,ikb)+pgeoh(jl,ikb))
2069 pmfuq(jl,ikb) = pmfub(jl)*pqu(jl,ikb)
2070 pmful(jl,ikb) = pmfub(jl)*plu(jl,ikb)
2071 end if
2072 end do
2073!
2074!-----------------------------------------------------------------
2075! 4. do ascent: subcloud layer (klab=1) ,clouds (klab=2)
2076! by doing first dry-adiabatic ascent and then
2077! by adjusting t,q and l accordingly in *cuadjtqn*,
2078! then check for buoyancy and set flags accordingly
2079!-----------------------------------------------------------------
2080!
2081 do jk=klevm1,3,-1
2082! specify cloud base values for midlevel convection
2083! in *cubasmc* in case there is not already convection
2084! ---------------------------------------------------------------------
2085 ik=jk
2086 call cubasmcn&
2087 & (klon, klev, klevm1, ik, pten,&
2088 & pqen, pqsen, puen, pven, pverv,&
2089 & pgeo, pgeoh, ldcum, ktype, klab, zlrain,&
2090 & pmfu, pmfub, kcbot, ptu,&
2091 & pqu, plu, puu, pvu, pmfus,&
2092 & pmfuq, pmful, pdmfup)
2093 is = 0
2094 jlm = 0
2095 do jl = 1,klon
2096 loflag(jl) = .false.
2097 zprecip(jl) = 0.
2098 llo1(jl) = .false.
2099 is = is + klab(jl,jk+1)
2100 if ( klab(jl,jk+1) == 0 ) klab(jl,jk) = 0
2101 if ( (ldcum(jl) .and. klab(jl,jk+1) == 2) .or. &
2102 (ktype(jl) == 3 .and. klab(jl,jk+1) == 1) ) then
2103 loflag(jl) = .true.
2104 jlm = jlm + 1
2105 jlx(jlm) = jl
2106 end if
2107 zph(jl) = paph(jl,jk)
2108 if ( ktype(jl) == 3 .and. jk == kcbot(jl) ) then
2109 zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons2
2110 if ( pmfub(jl) > zmfmax ) then
2111 zfac = zmfmax/pmfub(jl)
2112 pmfu(jl,jk+1) = pmfu(jl,jk+1)*zfac
2113 pmfus(jl,jk+1) = pmfus(jl,jk+1)*zfac
2114 pmfuq(jl,jk+1) = pmfuq(jl,jk+1)*zfac
2115 pmfub(jl) = zmfmax
2116 end if
2117 pmfub(jl)=min(pmfub(jl),zmfmax)
2118 end if
2119 end do
2120
2121 if(is.gt.0) llo3 = .true.
2122!
2123!* specify entrainment rates in *cuentr*
2124! -------------------------------------
2125 ik=jk
2126 call cuentrn(klon,klev,ik,kcbot,ldcum,llo3, &
2127 pgeoh,pmfu,zdmfen,zdmfde)
2128!
2129! do adiabatic ascent for entraining/detraining plume
2130 if(llo3) then
2131! -------------------------------------------------------
2132!
2133 do jl = 1,klon
2134 zqold(jl) = 0.
2135 end do
2136 do jll = 1 , jlm
2137 jl = jlx(jll)
2138 zdmfde(jl) = min(zdmfde(jl),0.75*pmfu(jl,jk+1))
2139 if ( jk == kcbot(jl) ) then
2140 zoentr(jl) = -1.75e-3*(min(1.,pqen(jl,jk)/pqsen(jl,jk)) - &
2141 1.)*(pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
2142 zoentr(jl) = min(0.4,zoentr(jl))*pmfu(jl,jk+1)
2143 end if
2144 if ( jk < kcbot(jl) ) then
2145 zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons2
2146 zxs = max(pmfu(jl,jk+1)-zmfmax,0.)
2147 wup(jl) = wup(jl) + kup(jl,jk+1)*(pap(jl,jk+1)-pap(jl,jk))
2148 zdpmean(jl) = zdpmean(jl) + pap(jl,jk+1) - pap(jl,jk)
2149 zdmfen(jl) = zoentr(jl)
2150 if ( ktype(jl) >= 2 ) then
2151 zdmfen(jl) = 2.0*zdmfen(jl)
2152 zdmfde(jl) = zdmfen(jl)
2153 end if
2154 zdmfde(jl) = zdmfde(jl) * &
2155 (1.6-min(1.,pqen(jl,jk)/pqsen(jl,jk)))
2156 zmftest = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl)
2157 zchange = max(zmftest-zmfmax,0.)
2158 zxe = max(zchange-zxs,0.)
2159 zdmfen(jl) = zdmfen(jl) - zxe
2160 zchange = zchange - zxe
2161 zdmfde(jl) = zdmfde(jl) + zchange
2162 end if
2163 pdmfen(jl,jk) = zdmfen(jl) - zdmfde(jl)
2164 pmfu(jl,jk) = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl)
2165 zqeen = pqenh(jl,jk+1)*zdmfen(jl)
2166 zseen = (cpd*ptenh(jl,jk+1)+pgeoh(jl,jk+1))*zdmfen(jl)
2167 zscde = (cpd*ptu(jl,jk+1)+pgeoh(jl,jk+1))*zdmfde(jl)
2168 zqude = pqu(jl,jk+1)*zdmfde(jl)
2169 plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl)
2170 zmfusk = pmfus(jl,jk+1) + zseen - zscde
2171 zmfuqk = pmfuq(jl,jk+1) + zqeen - zqude
2172 zmfulk = pmful(jl,jk+1) - plude(jl,jk)
2173 plu(jl,jk) = zmfulk*(1./max(cmfcmin,pmfu(jl,jk)))
2174 pqu(jl,jk) = zmfuqk*(1./max(cmfcmin,pmfu(jl,jk)))
2175 ptu(jl,jk) = (zmfusk * &
2176 (1./max(cmfcmin,pmfu(jl,jk)))-pgeoh(jl,jk))*rcpd
2177 ptu(jl,jk) = max(100.,ptu(jl,jk))
2178 ptu(jl,jk) = min(400.,ptu(jl,jk))
2179 zqold(jl) = pqu(jl,jk)
2180 zlrain(jl,jk) = zlrain(jl,jk+1)*(pmfu(jl,jk+1)-zdmfde(jl)) * &
2181 (1./max(cmfcmin,pmfu(jl,jk)))
2182 zluold(jl) = plu(jl,jk)
2183 end do
2184! reset to environmental values if below departure level
2185 do jl = 1,klon
2186 if ( jk > kdpl(jl) ) then
2187 ptu(jl,jk) = ptenh(jl,jk)
2188 pqu(jl,jk) = pqenh(jl,jk)
2189 plu(jl,jk) = 0.
2190 zluold(jl) = plu(jl,jk)
2191 end if
2192 end do
2193!* do corrections for moist ascent
2194!* by adjusting t,q and l in *cuadjtq*
2195!------------------------------------------------
2196 ik=jk
2197 icall=1
2198!
2199 if ( jlm > 0 ) then
2200 call cuadjtqn(klon,klev,ik,zph,ptu,pqu,loflag,icall)
2201 end if
2202! compute the upfraft speed in cloud layer
2203 do jll = 1 , jlm
2204 jl = jlx(jll)
2205 if ( pqu(jl,jk) /= zqold(jl) ) then
2206 plglac(jl,jk) = plu(jl,jk) * &
2207 ((1.-foealfa(ptu(jl,jk)))- &
2208 (1.-foealfa(ptu(jl,jk+1))))
2209 ptu(jl,jk) = ptu(jl,jk) + ralfdcp*plglac(jl,jk)
2210 end if
2211 end do
2212 do jll = 1 , jlm
2213 jl = jlx(jll)
2214 if ( pqu(jl,jk) /= zqold(jl) ) then
2215 klab(jl,jk) = 2
2216 plu(jl,jk) = plu(jl,jk) + zqold(jl) - pqu(jl,jk)
2217 zbc = ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk+1) - &
2218 zlrain(jl,jk+1))
2219 zbe = ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
2220 zbuo(jl,jk) = zbc - zbe
2221! set flags for the case of midlevel convection
2222 if ( ktype(jl) == 3 .and. klab(jl,jk+1) == 1 ) then
2223 if ( zbuo(jl,jk) > -0.5 ) then
2224 ldcum(jl) = .true.
2225 kctop(jl) = jk
2226 kup(jl,jk) = 0.5
2227 else
2228 klab(jl,jk) = 0
2229 pmfu(jl,jk) = 0.
2230 plude(jl,jk) = 0.
2231 plu(jl,jk) = 0.
2232 end if
2233 end if
2234 if ( klab(jl,jk+1) == 2 ) then
2235 if ( zbuo(jl,jk) < 0. ) then
2236 ptenh(jl,jk) = 0.5*(pten(jl,jk)+pten(jl,jk-1))
2237 pqenh(jl,jk) = 0.5*(pqen(jl,jk)+pqen(jl,jk-1))
2238 zbuo(jl,jk) = zbc - ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
2239 end if
2240 zbuoc = (zbuo(jl,jk) / &
2241 (ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk)))+zbuo(jl,jk+1) / &
2242 (ptenh(jl,jk+1)*(1.+vtmpc1*pqenh(jl,jk+1))))*0.5
2243 zdkbuo = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zfacbuo*zbuoc
2244! mixing and "pressure" gradient term in upper troposphere
2245 if ( zdmfen(jl) > 0. ) then
2246 zdken = min(1.,(1.+z_cwdrag)*zdmfen(jl) / &
2247 max(cmfcmin,pmfu(jl,jk+1)))
2248 else
2249 zdken = min(1.,(1.+z_cwdrag)*zdmfde(jl) / &
2250 max(cmfcmin,pmfu(jl,jk+1)))
2251 end if
2252 kup(jl,jk) = (kup(jl,jk+1)*(1.-zdken)+zdkbuo) / &
2253 (1.+zdken)
2254 if ( zbuo(jl,jk) < 0. ) then
2255 zkedke = kup(jl,jk)/max(1.e-10,kup(jl,jk+1))
2256 zkedke = max(0.,min(1.,zkedke))
2257 zmfun = sqrt(zkedke)*pmfu(jl,jk+1) !* (1.6-min(1.,pqen(jl,jk) / &
2258 ! pqsen(jl,jk)))
2259 zdmfde(jl) = max(zdmfde(jl),pmfu(jl,jk+1)-zmfun)
2260 plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl)
2261 pmfu(jl,jk) = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl)
2262 end if
2263 if ( zbuo(jl,jk) > -0.2 ) then
2264 ikb = kcbot(jl)
2265 zoentr(jl) = 1.75e-3*(0.3-(min(1.,pqen(jl,jk-1) / &
2266 pqsen(jl,jk-1))-1.))*(pgeoh(jl,jk-1)-pgeoh(jl,jk)) * &
2267 zrg*min(1.,pqsen(jl,jk)/pqsen(jl,ikb))**3
2268 zoentr(jl) = min(0.4,zoentr(jl))*pmfu(jl,jk)
2269 else
2270 zoentr(jl) = 0.
2271 end if
2272! erase values if below departure level
2273 if ( jk > kdpl(jl) ) then
2274 pmfu(jl,jk) = pmfu(jl,jk+1)
2275 kup(jl,jk) = 0.5
2276 end if
2277 if ( kup(jl,jk) > 0. .and. pmfu(jl,jk) > 0. ) then
2278 kctop(jl) = jk
2279 llo1(jl) = .true.
2280 else
2281 klab(jl,jk) = 0
2282 pmfu(jl,jk) = 0.
2283 kup(jl,jk) = 0.
2284 zdmfde(jl) = pmfu(jl,jk+1)
2285 plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl)
2286 end if
2287! save detrainment rates for updraught
2288 if ( pmfu(jl,jk+1) > 0. ) pmfude_rate(jl,jk) = zdmfde(jl)
2289 end if
2290 else if ( ktype(jl) == 2 .and. pqu(jl,jk) == zqold(jl) ) then
2291 klab(jl,jk) = 0
2292 pmfu(jl,jk) = 0.
2293 kup(jl,jk) = 0.
2294 zdmfde(jl) = pmfu(jl,jk+1)
2295 plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl)
2296 pmfude_rate(jl,jk) = zdmfde(jl)
2297 end if
2298 end do
2299
2300 do jl = 1,klon
2301 if ( llo1(jl) ) then
2302! conversions only proceeds if plu is greater than a threshold liquid water
2303! content of 0.3 g/kg over water and 0.5 g/kg over land to prevent precipitation
2304! generation from small water contents.
2305 if ( lndj(jl).eq.1 ) then
2306 zdshrd = 5.e-4
2307 else
2308 zdshrd = 3.e-4
2309 end if
2310 ikb=kcbot(jl)
2311 if ( plu(jl,jk) > zdshrd )then
2312! if ((paph(jl,ikb)-paph(jl,jk))>zdnoprc) then
2313 zwu = min(15.0,sqrt(2.*max(0.1,kup(jl,jk+1))))
2314 zprcon = zprcdgw/(0.75*zwu)
2315! PARAMETERS FOR BERGERON-FINDEISEN PROCESS (T < -5C)
2316 zdt = min(rtber-rtice,max(rtber-ptu(jl,jk),0.))
2317 zcbf = 1. + z_cprc2*sqrt(zdt)
2318 zzco = zprcon*zcbf
2319 zlcrit = zdshrd/zcbf
2320 zdfi = pgeoh(jl,jk) - pgeoh(jl,jk+1)
2321 zc = (plu(jl,jk)-zluold(jl))
2322 zarg = (plu(jl,jk)/zlcrit)**2
2323 if ( zarg < 25.0 ) then
2324 zd = zzco*(1.-exp(-zarg))*zdfi
2325 else
2326 zd = zzco*zdfi
2327 end if
2328 zint = exp(-zd)
2329 zlnew = zluold(jl)*zint + zc/zd*(1.-zint)
2330 zlnew = max(0.,min(plu(jl,jk),zlnew))
2331 zlnew = min(z_cldmax,zlnew)
2332 zprecip(jl) = max(0.,zluold(jl)+zc-zlnew)
2333 pdmfup(jl,jk) = zprecip(jl)*pmfu(jl,jk)
2334 zlrain(jl,jk) = zlrain(jl,jk) + zprecip(jl)
2335 plu(jl,jk) = zlnew
2336 end if
2337 end if
2338 end do
2339 do jl = 1, klon
2340 if ( llo1(jl) ) then
2341 if ( zlrain(jl,jk) > 0. ) then
2342 zvw = 21.18*zlrain(jl,jk)**0.2
2343 zvi = z_cwifrac*zvw
2344 zalfaw = foealfa(ptu(jl,jk))
2345 zvv = zalfaw*zvw + (1.-zalfaw)*zvi
2346 zrold = zlrain(jl,jk) - zprecip(jl)
2347 zc = zprecip(jl)
2348 zwu = min(15.0,sqrt(2.*max(0.1,kup(jl,jk))))
2349 zd = zvv/zwu
2350 zint = exp(-zd)
2351 zrnew = zrold*zint + zc/zd*(1.-zint)
2352 zrnew = max(0.,min(zlrain(jl,jk),zrnew))
2353 zlrain(jl,jk) = zrnew
2354 end if
2355 end if
2356 end do
2357 do jll = 1 , jlm
2358 jl = jlx(jll)
2359 pmful(jl,jk) = plu(jl,jk)*pmfu(jl,jk)
2360 pmfus(jl,jk) = (cpd*ptu(jl,jk)+pgeoh(jl,jk))*pmfu(jl,jk)
2361 pmfuq(jl,jk) = pqu(jl,jk)*pmfu(jl,jk)
2362 end do
2363 end if
2364 end do
2365!----------------------------------------------------------------------
2366! 5. final calculations
2367! ------------------
2368 do jl = 1,klon
2369 if ( kctop(jl) == -1 ) ldcum(jl) = .false.
2370 kcbot(jl) = max(kcbot(jl),kctop(jl))
2371 if ( ldcum(jl) ) then
2372 wup(jl) = max(1.e-2,wup(jl)/max(1.,zdpmean(jl)))
2373 wup(jl) = sqrt(2.*wup(jl))
2374 end if
2375 end do
2376
2377 return
2378 end subroutine cuascn
2379!---------------------------------------------------------
2380! level 3 souroutines
2381!--------------------------------------------------------
2382 subroutine cudlfsn &
2383 & (klon, klev, &
2384 & kcbot, kctop, lndj, ldcum, &
2385 & ptenh, pqenh, puen, pven, &
2386 & pten, pqsen, pgeo, &
2387 & pgeoh, paph, ptu, pqu, plu,&
2388 & puu, pvu, pmfub, prfl, &
2389 & ptd, pqd, pud, pvd, &
2390 & pmfd, pmfds, pmfdq, pdmfdp, &
2391 & kdtop, lddraf)
2392
2393! this routine calculates level of free sinking for
2394! cumulus downdrafts and specifies t,q,u and v values
2395
2396! m.tiedtke e.c.m.w.f. 12/86 modif. 12/89
2397
2398! purpose.
2399! --------
2400! to produce lfs-values for cumulus downdrafts
2401! for massflux cumulus parameterization
2402
2403! interface
2404! ---------
2405! this routine is called from *cumastr*.
2406! input are environmental values of t,q,u,v,p,phi
2407! and updraft values t,q,u and v and also
2408! cloud base massflux and cu-precipitation rate.
2409! it returns t,q,u and v values and massflux at lfs.
2410
2411! method.
2412
2413! check for negative buoyancy of air of equal parts of
2414! moist environmental air and cloud air.
2415
2416! parameter description units
2417! --------- ----------- -----
2418! input parameters (integer):
2419
2420! *klon* number of grid points per packet
2421! *klev* number of levels
2422! *kcbot* cloud base level
2423! *kctop* cloud top level
2424
2425! input parameters (logical):
2426
2427! *lndj* land sea mask (1 for land)
2428! *ldcum* flag: .true. for convective points
2429
2430! input parameters (real(kind=kind_phys)):
2431
2432! *ptenh* env. temperature (t+1) on half levels k
2433! *pqenh* env. spec. humidity (t+1) on half levels kg/kg
2434! *puen* provisional environment u-velocity (t+1) m/s
2435! *pven* provisional environment v-velocity (t+1) m/s
2436! *pten* provisional environment temperature (t+1) k
2437! *pqsen* environment spec. saturation humidity (t+1) kg/kg
2438! *pgeo* geopotential m2/s2
2439! *pgeoh* geopotential on half levels m2/s2
2440! *paph* provisional pressure on half levels pa
2441! *ptu* temperature in updrafts k
2442! *pqu* spec. humidity in updrafts kg/kg
2443! *plu* liquid water content in updrafts kg/kg
2444! *puu* u-velocity in updrafts m/s
2445! *pvu* v-velocity in updrafts m/s
2446! *pmfub* massflux in updrafts at cloud base kg/(m2*s)
2447
2448! updated parameters (real(kind=kind_phys)):
2449
2450! *prfl* precipitation rate kg/(m2*s)
2451
2452! output parameters (real(kind=kind_phys)):
2453
2454! *ptd* temperature in downdrafts k
2455! *pqd* spec. humidity in downdrafts kg/kg
2456! *pud* u-velocity in downdrafts m/s
2457! *pvd* v-velocity in downdrafts m/s
2458! *pmfd* massflux in downdrafts kg/(m2*s)
2459! *pmfds* flux of dry static energy in downdrafts j/(m2*s)
2460! *pmfdq* flux of spec. humidity in downdrafts kg/(m2*s)
2461! *pdmfdp* flux difference of precip. in downdrafts kg/(m2*s)
2462
2463! output parameters (integer):
2464
2465! *kdtop* top level of downdrafts
2466
2467! output parameters (logical):
2468
2469! *lddraf* .true. if downdrafts exist
2470
2471! externals
2472! ---------
2473! *cuadjtq* for calculating wet bulb t and q at lfs
2474!----------------------------------------------------------------------
2475 implicit none
2476
2477 integer klev,klon
2478 real(kind=kind_phys) ptenh(klon,klev), pqenh(klon,klev), &
2479 & puen(klon,klev), pven(klon,klev), &
2480 & pten(klon,klev), pqsen(klon,klev), &
2481 & pgeo(klon,klev), &
2482 & pgeoh(klon,klev+1), paph(klon,klev+1),&
2483 & ptu(klon,klev), pqu(klon,klev), &
2484 & puu(klon,klev), pvu(klon,klev), &
2485 & plu(klon,klev), &
2486 & pmfub(klon), prfl(klon)
2487
2488 real(kind=kind_phys) ptd(klon,klev), pqd(klon,klev), &
2489 & pud(klon,klev), pvd(klon,klev), &
2490 & pmfd(klon,klev), pmfds(klon,klev), &
2491 & pmfdq(klon,klev), pdmfdp(klon,klev)
2492 integer kcbot(klon), kctop(klon), &
2493 & kdtop(klon), ikhsmin(klon)
2494 logical ldcum(klon), &
2495 & lddraf(klon)
2496 integer lndj(klon)
2497
2498 real(kind=kind_phys) ztenwb(klon,klev), zqenwb(klon,klev), &
2499 & zcond(klon), zph(klon), &
2500 & zhsmin(klon)
2501 logical llo2(klon)
2502! local variables
2503 integer jl,jk
2504 integer is,ik,icall,ike
2505 real(kind=kind_phys) zhsk,zttest,zqtest,zbuo,zmftop
2506
2507!----------------------------------------------------------------------
2508
2509! 1. set default values for downdrafts
2510! ---------------------------------
2511 do jl=1,klon
2512 lddraf(jl)=.false.
2513 kdtop(jl)=klev+1
2514 ikhsmin(jl)=klev+1
2515 zhsmin(jl)=1.e8
2516 enddo
2517!----------------------------------------------------------------------
2518
2519! 2. determine level of free sinking:
2520! downdrafts shall start at model level of minimum
2521! of saturation moist static energy or below
2522! respectively
2523
2524! for every point and proceed as follows:
2525
2526! (1) determine level of minimum of hs
2527! (2) determine wet bulb environmental t and q
2528! (3) do mixing with cumulus cloud air
2529! (4) check for negative buoyancy
2530! (5) if buoyancy>0 repeat (2) to (4) for next
2531! level below
2532
2533! the assumption is that air of downdrafts is mixture
2534! of 50% cloud air + 50% environmental air at wet bulb
2535! temperature (i.e. which became saturated due to
2536! evaporation of rain and cloud water)
2537! ----------------------------------------------------
2538 do jk=3,klev-2
2539 do jl=1,klon
2540 zhsk=cpd*pten(jl,jk)+pgeo(jl,jk) + &
2541 & foelhm(pten(jl,jk))*pqsen(jl,jk)
2542 if(zhsk .lt. zhsmin(jl)) then
2543 zhsmin(jl) = zhsk
2544 ikhsmin(jl)= jk
2545 end if
2546 end do
2547 end do
2548
2549
2550 ike=klev-3
2551 do jk=3,ike
2552
2553! 2.1 calculate wet-bulb temperature and moisture
2554! for environmental air in *cuadjtq*
2555! -------------------------------------------
2556 is=0
2557 do jl=1,klon
2558 ztenwb(jl,jk)=ptenh(jl,jk)
2559 zqenwb(jl,jk)=pqenh(jl,jk)
2560 zph(jl)=paph(jl,jk)
2561 llo2(jl)=ldcum(jl).and.prfl(jl).gt.0..and..not.lddraf(jl).and. &
2562 & (jk.lt.kcbot(jl).and.jk.gt.kctop(jl)).and. jk.ge.ikhsmin(jl)
2563 if(llo2(jl))then
2564 is=is+1
2565 endif
2566 enddo
2567 if(is.eq.0) cycle
2568
2569 ik=jk
2570 icall=2
2571 call cuadjtqn &
2572 & ( klon, klev, ik, zph, ztenwb, zqenwb, llo2, icall)
2573
2574! 2.2 do mixing of cumulus and environmental air
2575! and check for negative buoyancy.
2576! then set values for downdraft at lfs.
2577! ----------------------------------------
2578 do jl=1,klon
2579 if(llo2(jl)) then
2580 zttest=0.5*(ptu(jl,jk)+ztenwb(jl,jk))
2581 zqtest=0.5*(pqu(jl,jk)+zqenwb(jl,jk))
2582 zbuo=zttest*(1.+vtmpc1 *zqtest)- &
2583 & ptenh(jl,jk)*(1.+vtmpc1 *pqenh(jl,jk))
2584 zcond(jl)=pqenh(jl,jk)-zqenwb(jl,jk)
2585 zmftop=-cmfdeps*pmfub(jl)
2586 if(zbuo.lt.0..and.prfl(jl).gt.10.*zmftop*zcond(jl)) then
2587 kdtop(jl)=jk
2588 lddraf(jl)=.true.
2589 ptd(jl,jk)=zttest
2590 pqd(jl,jk)=zqtest
2591 pmfd(jl,jk)=zmftop
2592 pmfds(jl,jk)=pmfd(jl,jk)*(cpd*ptd(jl,jk)+pgeoh(jl,jk))
2593 pmfdq(jl,jk)=pmfd(jl,jk)*pqd(jl,jk)
2594 pdmfdp(jl,jk-1)=-0.5*pmfd(jl,jk)*zcond(jl)
2595 prfl(jl)=prfl(jl)+pdmfdp(jl,jk-1)
2596 endif
2597 endif
2598 enddo
2599
2600 enddo
2601
2602 return
2603 end subroutine cudlfsn
2604
2605!---------------------------------------------------------
2606! level 3 souroutines
2607!--------------------------------------------------------
2608!**********************************************
2609! subroutine cuddrafn
2610!**********************************************
2611 subroutine cuddrafn &
2612 & ( klon, klev, lddraf &
2613 & , ptenh, pqenh, puen, pven &
2614 & , pgeo, pgeoh, paph, prfl &
2615 & , ptd, pqd, pud, pvd, pmfu &
2616 & , pmfd, pmfds, pmfdq, pdmfdp, pmfdde_rate )
2617
2618! this routine calculates cumulus downdraft descent
2619
2620! m.tiedtke e.c.m.w.f. 12/86 modif. 12/89
2621
2622! purpose.
2623! --------
2624! to produce the vertical profiles for cumulus downdrafts
2625! (i.e. t,q,u and v and fluxes)
2626
2627! interface
2628! ---------
2629
2630! this routine is called from *cumastr*.
2631! input is t,q,p,phi,u,v at half levels.
2632! it returns fluxes of s,q and evaporation rate
2633! and u,v at levels where downdraft occurs
2634
2635! method.
2636! --------
2637! calculate moist descent for entraining/detraining plume by
2638! a) moving air dry-adiabatically to next level below and
2639! b) correcting for evaporation to obtain saturated state.
2640
2641! parameter description units
2642! --------- ----------- -----
2643! input parameters (integer):
2644
2645! *klon* number of grid points per packet
2646! *klev* number of levels
2647
2648! input parameters (logical):
2649
2650! *lddraf* .true. if downdrafts exist
2651
2652! input parameters (real(kind=kind_phys)):
2653
2654! *ptenh* env. temperature (t+1) on half levels k
2655! *pqenh* env. spec. humidity (t+1) on half levels kg/kg
2656! *puen* provisional environment u-velocity (t+1) m/s
2657! *pven* provisional environment v-velocity (t+1) m/s
2658! *pgeo* geopotential m2/s2
2659! *pgeoh* geopotential on half levels m2/s2
2660! *paph* provisional pressure on half levels pa
2661! *pmfu* massflux updrafts kg/(m2*s)
2662
2663! updated parameters (real(kind=kind_phys)):
2664
2665! *prfl* precipitation rate kg/(m2*s)
2666
2667! output parameters (real(kind=kind_phys)):
2668
2669! *ptd* temperature in downdrafts k
2670! *pqd* spec. humidity in downdrafts kg/kg
2671! *pud* u-velocity in downdrafts m/s
2672! *pvd* v-velocity in downdrafts m/s
2673! *pmfd* massflux in downdrafts kg/(m2*s)
2674! *pmfds* flux of dry static energy in downdrafts j/(m2*s)
2675! *pmfdq* flux of spec. humidity in downdrafts kg/(m2*s)
2676! *pdmfdp* flux difference of precip. in downdrafts kg/(m2*s)
2677
2678! externals
2679! ---------
2680! *cuadjtq* for adjusting t and q due to evaporation in
2681! saturated descent
2682!----------------------------------------------------------------------
2683 implicit none
2684
2685 integer klev,klon
2686 real(kind=kind_phys) ptenh(klon,klev), pqenh(klon,klev), &
2687 & puen(klon,klev), pven(klon,klev), &
2688 & pgeoh(klon,klev+1), paph(klon,klev+1), &
2689 & pgeo(klon,klev), pmfu(klon,klev)
2690
2691 real(kind=kind_phys) ptd(klon,klev), pqd(klon,klev), &
2692 & pud(klon,klev), pvd(klon,klev), &
2693 & pmfd(klon,klev), pmfds(klon,klev), &
2694 & pmfdq(klon,klev), pdmfdp(klon,klev), &
2695 & prfl(klon)
2696 real(kind=kind_phys) pmfdde_rate(klon,klev)
2697 logical lddraf(klon)
2698
2699 real(kind=kind_phys) zdmfen(klon), zdmfde(klon), &
2700 & zcond(klon), zoentr(klon), &
2701 & zbuoy(klon)
2702 real(kind=kind_phys) zph(klon)
2703 logical llo2(klon)
2704 logical llo1
2705! local variables
2706 integer jl,jk
2707 integer is,ik,icall,ike, itopde(klon)
2708 real(kind=kind_phys) zentr,zdz,zzentr,zseen,zqeen,zsdde,zqdde,zdmfdp
2709 real(kind=kind_phys) zmfdsk,zmfdqk,zbuo,zrain,zbuoyz,zmfduk,zmfdvk
2710
2711!----------------------------------------------------------------------
2712! 1. calculate moist descent for cumulus downdraft by
2713! (a) calculating entrainment/detrainment rates,
2714! including organized entrainment dependent on
2715! negative buoyancy and assuming
2716! linear decrease of massflux in pbl
2717! (b) doing moist descent - evaporative cooling
2718! and moistening is calculated in *cuadjtq*
2719! (c) checking for negative buoyancy and
2720! specifying final t,q,u,v and downward fluxes
2721! -------------------------------------------------
2722 do jl=1,klon
2723 zoentr(jl)=0.
2724 zbuoy(jl)=0.
2725 zdmfen(jl)=0.
2726 zdmfde(jl)=0.
2727 enddo
2728
2729 do jk=klev,1,-1
2730 do jl=1,klon
2731 pmfdde_rate(jl,jk) = 0.
2732 if((paph(jl,klev+1)-paph(jl,jk)).lt. 60.e2) itopde(jl) = jk
2733 end do
2734 end do
2735
2736 do jk=3,klev
2737 is=0
2738 do jl=1,klon
2739 zph(jl)=paph(jl,jk)
2740 llo2(jl)=lddraf(jl).and.pmfd(jl,jk-1).lt.0.
2741 if(llo2(jl)) then
2742 is=is+1
2743 endif
2744 enddo
2745 if(is.eq.0) cycle
2746
2747 do jl=1,klon
2748 if(llo2(jl)) then
2749 zentr = entrdd*pmfd(jl,jk-1)*(pgeoh(jl,jk-1)-pgeoh(jl,jk))*zrg
2750 zdmfen(jl)=zentr
2751 zdmfde(jl)=zentr
2752 endif
2753 enddo
2754
2755 do jl=1,klon
2756 if(llo2(jl)) then
2757 if(jk.gt.itopde(jl)) then
2758 zdmfen(jl)=0.
2759 zdmfde(jl)=pmfd(jl,itopde(jl))* &
2760 & (paph(jl,jk)-paph(jl,jk-1))/ &
2761 & (paph(jl,klev+1)-paph(jl,itopde(jl)))
2762 endif
2763 endif
2764 enddo
2765
2766 do jl=1,klon
2767 if(llo2(jl)) then
2768 if(jk.le.itopde(jl)) then
2769 zdz=-(pgeoh(jl,jk-1)-pgeoh(jl,jk))*zrg
2770 zzentr=zoentr(jl)*zdz*pmfd(jl,jk-1)
2771 zdmfen(jl)=zdmfen(jl)+zzentr
2772 zdmfen(jl)=max(zdmfen(jl),0.3*pmfd(jl,jk-1))
2773 zdmfen(jl)=max(zdmfen(jl),-0.75*pmfu(jl,jk)- &
2774 & (pmfd(jl,jk-1)-zdmfde(jl)))
2775 zdmfen(jl)=min(zdmfen(jl),0.)
2776 endif
2777 endif
2778 enddo
2779
2780 do jl=1,klon
2781 if(llo2(jl)) then
2782 pmfd(jl,jk)=pmfd(jl,jk-1)+zdmfen(jl)-zdmfde(jl)
2783 zseen=(cpd*ptenh(jl,jk-1)+pgeoh(jl,jk-1))*zdmfen(jl)
2784 zqeen=pqenh(jl,jk-1)*zdmfen(jl)
2785 zsdde=(cpd*ptd(jl,jk-1)+pgeoh(jl,jk-1))*zdmfde(jl)
2786 zqdde=pqd(jl,jk-1)*zdmfde(jl)
2787 zmfdsk=pmfds(jl,jk-1)+zseen-zsdde
2788 zmfdqk=pmfdq(jl,jk-1)+zqeen-zqdde
2789 pqd(jl,jk)=zmfdqk*(1./min(-cmfcmin,pmfd(jl,jk)))
2790 ptd(jl,jk)=(zmfdsk*(1./min(-cmfcmin,pmfd(jl,jk)))-&
2791 & pgeoh(jl,jk))*rcpd
2792 ptd(jl,jk)=min(400.,ptd(jl,jk))
2793 ptd(jl,jk)=max(100.,ptd(jl,jk))
2794 zcond(jl)=pqd(jl,jk)
2795 endif
2796 enddo
2797
2798 ik=jk
2799 icall=2
2800 call cuadjtqn(klon, klev, ik, zph, ptd, pqd, llo2, icall )
2801
2802 do jl=1,klon
2803 if(llo2(jl)) then
2804 zcond(jl)=zcond(jl)-pqd(jl,jk)
2805 zbuo=ptd(jl,jk)*(1.+vtmpc1 *pqd(jl,jk))- &
2806 & ptenh(jl,jk)*(1.+vtmpc1 *pqenh(jl,jk))
2807 if(prfl(jl).gt.0..and.pmfu(jl,jk).gt.0.) then
2808 zrain=prfl(jl)/pmfu(jl,jk)
2809 zbuo=zbuo-ptd(jl,jk)*zrain
2810 endif
2811 if(zbuo.ge.0 .or. prfl(jl).le.(pmfd(jl,jk)*zcond(jl))) then
2812 pmfd(jl,jk)=0.
2813 zbuo=0.
2814 endif
2815 pmfds(jl,jk)=(cpd*ptd(jl,jk)+pgeoh(jl,jk))*pmfd(jl,jk)
2816 pmfdq(jl,jk)=pqd(jl,jk)*pmfd(jl,jk)
2817 zdmfdp=-pmfd(jl,jk)*zcond(jl)
2818 pdmfdp(jl,jk-1)=zdmfdp
2819 prfl(jl)=prfl(jl)+zdmfdp
2820
2821! compute organized entrainment for use at next level
2822 zbuoyz=zbuo/ptenh(jl,jk)
2823 zbuoyz=min(zbuoyz,0.0)
2824 zdz=-(pgeo(jl,jk-1)-pgeo(jl,jk))
2825 zbuoy(jl)=zbuoy(jl)+zbuoyz*zdz
2826 zoentr(jl)=g*zbuoyz*0.5/(1.+zbuoy(jl))
2827 pmfdde_rate(jl,jk) = -zdmfde(jl)
2828 endif
2829 enddo
2830
2831 enddo
2832
2833 return
2834 end subroutine cuddrafn
2835!---------------------------------------------------------
2836! level 3 souroutines
2837!--------------------------------------------------------
2838 subroutine cuflxn &
2839 & ( klon, klev, ztmst &
2840 & , pten, pqen, pqsen, ptenh, pqenh &
2841 & , paph, pap, pgeoh, lndj, ldcum &
2842 & , kcbot, kctop, kdtop, ktopm2 &
2843 & , ktype, lddraf &
2844 & , pmfu, pmfd, pmfus, pmfds &
2845 & , pmfuq, pmfdq, pmful, plude &
2846 & , pdmfup, pdmfdp, pdpmel, plglac &
2847 & , prain, pmfdde_rate, pmflxr, pmflxs )
2848
2849! m.tiedtke e.c.m.w.f. 7/86 modif. 12/89
2850
2851! purpose
2852! -------
2853
2854! this routine does the final calculation of convective
2855! fluxes in the cloud layer and in the subcloud layer
2856
2857! interface
2858! ---------
2859! this routine is called from *cumastr*.
2860
2861
2862! parameter description units
2863! --------- ----------- -----
2864! input parameters (integer):
2865
2866! *klon* number of grid points per packet
2867! *klev* number of levels
2868! *kcbot* cloud base level
2869! *kctop* cloud top level
2870! *kdtop* top level of downdrafts
2871
2872! input parameters (logical):
2873
2874! *lndj* land sea mask (1 for land)
2875! *ldcum* flag: .true. for convective points
2876
2877! input parameters (real(kind=kind_phys)):
2878
2879! *ztmst* time step for the physics s
2880! *pten* provisional environment temperature (t+1) k
2881! *pqen* provisional environment spec. humidity (t+1) kg/kg
2882! *pqsen* environment spec. saturation humidity (t+1) kg/kg
2883! *ptenh* env. temperature (t+1) on half levels k
2884! *pqenh* env. spec. humidity (t+1) on half levels kg/kg
2885! *paph* provisional pressure on half levels pa
2886! *pap* provisional pressure on full levels pa
2887! *pgeoh* geopotential on half levels m2/s2
2888
2889! updated parameters (integer):
2890
2891! *ktype* set to zero if ldcum=.false.
2892
2893! updated parameters (logical):
2894
2895! *lddraf* set to .false. if ldcum=.false. or kdtop<kctop
2896
2897! updated parameters (real(kind=kind_phys)):
2898
2899! *pmfu* massflux in updrafts kg/(m2*s)
2900! *pmfd* massflux in downdrafts kg/(m2*s)
2901! *pmfus* flux of dry static energy in updrafts j/(m2*s)
2902! *pmfds* flux of dry static energy in downdrafts j/(m2*s)
2903! *pmfuq* flux of spec. humidity in updrafts kg/(m2*s)
2904! *pmfdq* flux of spec. humidity in downdrafts kg/(m2*s)
2905! *pmful* flux of liquid water in updrafts kg/(m2*s)
2906! *plude* detrained liquid water kg/(m3*s)
2907! *pdmfup* flux difference of precip. in updrafts kg/(m2*s)
2908! *pdmfdp* flux difference of precip. in downdrafts kg/(m2*s)
2909
2910! output parameters (real(kind=kind_phys)):
2911
2912! *pdpmel* change in precip.-fluxes due to melting kg/(m2*s)
2913! *plglac* flux of frozen cloud water in updrafts kg/(m2*s)
2914! *pmflxr* convective rain flux kg/(m2*s)
2915! *pmflxs* convective snow flux kg/(m2*s)
2916! *prain* total precip. produced in conv. updrafts kg/(m2*s)
2917! (no evaporation in downdrafts)
2918
2919! externals
2920! ---------
2921! none
2922!----------------------------------------------------------------------
2923 implicit none
2924
2925 integer klev,klon,ktopm2
2926 real(kind=kind_phys) pten(klon,klev), ptenh(klon,klev), &
2927 & pqen(klon,klev), pqsen(klon,klev), &
2928 & pqenh(klon,klev), pap(klon,klev), &
2929 & paph(klon,klev+1), pgeoh(klon,klev+1), &
2930 & plglac(klon,klev)
2931
2932 real(kind=kind_phys) pmfu(klon,klev), pmfd(klon,klev), &
2933 & pmfus(klon,klev), pmfds(klon,klev), &
2934 & pmfuq(klon,klev), pmfdq(klon,klev), &
2935 & pdmfup(klon,klev), pdmfdp(klon,klev), &
2936 & pdpmel(klon,klev), prain(klon), &
2937 & pmful(klon,klev), plude(klon,klev), &
2938 & pmflxr(klon,klev+1), pmflxs(klon,klev+1)
2939 real(kind=kind_phys) pmfdde_rate(klon,klev)
2940 integer kcbot(klon), kctop(klon), &
2941 & kdtop(klon), ktype(klon)
2942 logical lddraf(klon), &
2943 & ldcum(klon)
2944 integer lndj(klon)
2945! local variables
2946 integer jl,jk
2947 integer is,ik,icall,ike,ikb
2948 real(kind=kind_phys) ztmst,ztaumel,zcons1a,zcons1,zcons2,zcucov,zcpecons
2949 real(kind=kind_phys) zalfaw,zrfl,zdrfl1,zrnew,zrmin,zrfln,zdrfl,zdenom
2950 real(kind=kind_phys) zpdr,zpds,zzp,zfac,zsnmlt
2951 real(kind=kind_phys) rhevap(klon)
2952 integer idbas(klon)
2953 logical llddraf
2954!--------------------------------------------------------------------
2955!* specify constants
2956
2957 ztaumel=18000.
2958 zcons1a=cpd/(alf*g*ztaumel)
2959 zcons2=3./(g*ztmst)
2960 zcucov=0.05
2961 zcpecons=5.44e-4/g
2962
2963!* 1.0 determine final convective fluxes
2964! ---------------------------------
2965 do jl=1,klon
2966 prain(jl)=0.
2967 if(.not.ldcum(jl).or.kdtop(jl).lt.kctop(jl)) lddraf(jl)=.false.
2968 if(.not.ldcum(jl)) ktype(jl)=0
2969 idbas(jl) = klev
2970 if(lndj(jl) .eq. 1) then
2971 rhevap(jl) = 0.7
2972 else
2973 rhevap(jl) = 0.9
2974 end if
2975 enddo
2976
2977 do jk=ktopm2,klev
2978 ikb = min(jk+1,klev)
2979 do jl=1,klon
2980 pmflxr(jl,jk) = 0.
2981 pmflxs(jl,jk) = 0.
2982 pdpmel(jl,jk) = 0.
2983 if(ldcum(jl).and.jk.ge.kctop(jl)) then
2984 pmfus(jl,jk)=pmfus(jl,jk)-pmfu(jl,jk)* &
2985 & (cpd*ptenh(jl,jk)+pgeoh(jl,jk))
2986 pmfuq(jl,jk)=pmfuq(jl,jk)-pmfu(jl,jk)*pqenh(jl,jk)
2987 plglac(jl,jk)=pmfu(jl,jk)*plglac(jl,jk)
2988 llddraf = lddraf(jl) .and. jk >= kdtop(jl)
2989 if ( llddraf .and.jk.ge.kdtop(jl)) then
2990 pmfds(jl,jk) = pmfds(jl,jk)-pmfd(jl,jk) * &
2991 (cpd*ptenh(jl,jk)+pgeoh(jl,jk))
2992 pmfdq(jl,jk) = pmfdq(jl,jk)-pmfd(jl,jk)*pqenh(jl,jk)
2993 else
2994 pmfd(jl,jk) = 0.
2995 pmfds(jl,jk) = 0.
2996 pmfdq(jl,jk) = 0.
2997 pdmfdp(jl,jk-1) = 0.
2998 end if
2999 if ( llddraf .and. pmfd(jl,jk) < 0. .and. &
3000 abs(pmfd(jl,ikb)) < 1.e-20 ) then
3001 idbas(jl) = jk
3002 end if
3003 else
3004 pmfu(jl,jk)=0.
3005 pmfd(jl,jk)=0.
3006 pmfus(jl,jk)=0.
3007 pmfds(jl,jk)=0.
3008 pmfuq(jl,jk)=0.
3009 pmfdq(jl,jk)=0.
3010 pmful(jl,jk)=0.
3011 plglac(jl,jk)=0.
3012 pdmfup(jl,jk-1)=0.
3013 pdmfdp(jl,jk-1)=0.
3014 plude(jl,jk-1)=0.
3015 endif
3016 enddo
3017 enddo
3018
3019 do jl=1,klon
3020 pmflxr(jl,klev+1) = 0.
3021 pmflxs(jl,klev+1) = 0.
3022 end do
3023 do jl=1,klon
3024 if(ldcum(jl)) then
3025 ikb=kcbot(jl)
3026 ik=ikb+1
3027 zzp=((paph(jl,klev+1)-paph(jl,ik))/ &
3028 & (paph(jl,klev+1)-paph(jl,ikb)))
3029 if(ktype(jl).eq.3) then
3030 zzp=zzp**2
3031 endif
3032 pmfu(jl,ik)=pmfu(jl,ikb)*zzp
3033 pmfus(jl,ik)=(pmfus(jl,ikb)- &
3034 & foelhm(ptenh(jl,ikb))*pmful(jl,ikb))*zzp
3035 pmfuq(jl,ik)=(pmfuq(jl,ikb)+pmful(jl,ikb))*zzp
3036 pmful(jl,ik)=0.
3037 endif
3038 enddo
3039
3040 do jk=ktopm2,klev
3041 do jl=1,klon
3042 if(ldcum(jl).and.jk.gt.kcbot(jl)+1) then
3043 ikb=kcbot(jl)+1
3044 zzp=((paph(jl,klev+1)-paph(jl,jk))/ &
3045 & (paph(jl,klev+1)-paph(jl,ikb)))
3046 if(ktype(jl).eq.3) then
3047 zzp=zzp**2
3048 endif
3049 pmfu(jl,jk)=pmfu(jl,ikb)*zzp
3050 pmfus(jl,jk)=pmfus(jl,ikb)*zzp
3051 pmfuq(jl,jk)=pmfuq(jl,ikb)*zzp
3052 pmful(jl,jk)=0.
3053 endif
3054 ik = idbas(jl)
3055 llddraf = lddraf(jl) .and. jk > ik .and. ik < klev
3056 if ( llddraf .and. ik == kcbot(jl)+1 ) then
3057 zzp = ((paph(jl,klev+1)-paph(jl,jk))/(paph(jl,klev+1)-paph(jl,ik)))
3058 if ( ktype(jl) == 3 ) zzp = zzp*zzp
3059 pmfd(jl,jk) = pmfd(jl,ik)*zzp
3060 pmfds(jl,jk) = pmfds(jl,ik)*zzp
3061 pmfdq(jl,jk) = pmfdq(jl,ik)*zzp
3062 pmfdde_rate(jl,jk) = -(pmfd(jl,jk-1)-pmfd(jl,jk))
3063 else if ( llddraf .and. ik /= kcbot(jl)+1 .and. jk == ik+1 ) then
3064 pmfdde_rate(jl,jk) = -(pmfd(jl,jk-1)-pmfd(jl,jk))
3065 end if
3066 enddo
3067 enddo
3068!* 2. calculate rain/snow fall rates
3069!* calculate melting of snow
3070!* calculate evaporation of precip
3071! -------------------------------
3072
3073 do jk=ktopm2,klev
3074 do jl=1,klon
3075 if(ldcum(jl) .and. jk >=kctop(jl)-1 ) then
3076 prain(jl)=prain(jl)+pdmfup(jl,jk)
3077 if(pmflxs(jl,jk).gt.0..and.pten(jl,jk).gt.tmelt) then
3078 zcons1=zcons1a*(1.+0.5*(pten(jl,jk)-tmelt))
3079 zfac=zcons1*(paph(jl,jk+1)-paph(jl,jk))
3080 zsnmlt=min(pmflxs(jl,jk),zfac*(pten(jl,jk)-tmelt))
3081 pdpmel(jl,jk)=zsnmlt
3082 pqsen(jl,jk)=foeewm(pten(jl,jk)-zsnmlt/zfac)/pap(jl,jk)
3083 endif
3084 zalfaw=foealfa(pten(jl,jk))
3085 !
3086 ! No liquid precipitation above melting level
3087 !
3088 if ( pten(jl,jk) < tmelt .and. zalfaw > 0. ) then
3089 plglac(jl,jk) = plglac(jl,jk)+zalfaw*(pdmfup(jl,jk)+pdmfdp(jl,jk))
3090 zalfaw = 0.
3091 end if
3092 pmflxr(jl,jk+1)=pmflxr(jl,jk)+zalfaw* &
3093 & (pdmfup(jl,jk)+pdmfdp(jl,jk))+pdpmel(jl,jk)
3094 pmflxs(jl,jk+1)=pmflxs(jl,jk)+(1.-zalfaw)* &
3095 & (pdmfup(jl,jk)+pdmfdp(jl,jk))-pdpmel(jl,jk)
3096 if(pmflxr(jl,jk+1)+pmflxs(jl,jk+1).lt.0.0) then
3097 pdmfdp(jl,jk)=-(pmflxr(jl,jk)+pmflxs(jl,jk)+pdmfup(jl,jk))
3098 pmflxr(jl,jk+1)=0.0
3099 pmflxs(jl,jk+1)=0.0
3100 pdpmel(jl,jk) =0.0
3101 else if ( pmflxr(jl,jk+1) < 0. ) then
3102 pmflxs(jl,jk+1) = pmflxs(jl,jk+1)+pmflxr(jl,jk+1)
3103 pmflxr(jl,jk+1) = 0.
3104 else if ( pmflxs(jl,jk+1) < 0. ) then
3105 pmflxr(jl,jk+1) = pmflxr(jl,jk+1)+pmflxs(jl,jk+1)
3106 pmflxs(jl,jk+1) = 0.
3107 end if
3108 endif
3109 enddo
3110 enddo
3111 do jk=ktopm2,klev
3112 do jl=1,klon
3113 if(ldcum(jl).and.jk.ge.kcbot(jl)) then
3114 zrfl=pmflxr(jl,jk)+pmflxs(jl,jk)
3115 if(zrfl.gt.1.e-20) then
3116 zdrfl1=zcpecons*max(0.,pqsen(jl,jk)-pqen(jl,jk))*zcucov* &
3117 & (sqrt(paph(jl,jk)/paph(jl,klev+1))/5.09e-3* &
3118 & zrfl/zcucov)**0.5777* &
3119 & (paph(jl,jk+1)-paph(jl,jk))
3120 zrnew=zrfl-zdrfl1
3121 zrmin=zrfl-zcucov*max(0.,rhevap(jl)*pqsen(jl,jk) &
3122 & -pqen(jl,jk)) *zcons2*(paph(jl,jk+1)-paph(jl,jk))
3123 zrnew=max(zrnew,zrmin)
3124 zrfln=max(zrnew,0.)
3125 zdrfl=min(0.,zrfln-zrfl)
3126 zdenom=1./max(1.e-20,pmflxr(jl,jk)+pmflxs(jl,jk))
3127 zalfaw=foealfa(pten(jl,jk))
3128 if ( pten(jl,jk) < tmelt ) zalfaw = 0.
3129 zpdr=zalfaw*pdmfdp(jl,jk)
3130 zpds=(1.-zalfaw)*pdmfdp(jl,jk)
3131 pmflxr(jl,jk+1)=pmflxr(jl,jk)+zpdr &
3132 & +pdpmel(jl,jk)+zdrfl*pmflxr(jl,jk)*zdenom
3133 pmflxs(jl,jk+1)=pmflxs(jl,jk)+zpds &
3134 & -pdpmel(jl,jk)+zdrfl*pmflxs(jl,jk)*zdenom
3135 pdmfup(jl,jk)=pdmfup(jl,jk)+zdrfl
3136 if ( pmflxr(jl,jk+1)+pmflxs(jl,jk+1) < 0. ) then
3137 pdmfup(jl,jk) = pdmfup(jl,jk)-(pmflxr(jl,jk+1)+pmflxs(jl,jk+1))
3138 pmflxr(jl,jk+1) = 0.
3139 pmflxs(jl,jk+1) = 0.
3140 pdpmel(jl,jk) = 0.
3141 else if ( pmflxr(jl,jk+1) < 0. ) then
3142 pmflxs(jl,jk+1) = pmflxs(jl,jk+1)+pmflxr(jl,jk+1)
3143 pmflxr(jl,jk+1) = 0.
3144 else if ( pmflxs(jl,jk+1) < 0. ) then
3145 pmflxr(jl,jk+1) = pmflxr(jl,jk+1)+pmflxs(jl,jk+1)
3146 pmflxs(jl,jk+1) = 0.
3147 end if
3148 else
3149 pmflxr(jl,jk+1)=0.0
3150 pmflxs(jl,jk+1)=0.0
3151 pdmfdp(jl,jk)=0.0
3152 pdpmel(jl,jk)=0.0
3153 endif
3154 endif
3155 enddo
3156 enddo
3157
3158 return
3159 end subroutine cuflxn
3160!---------------------------------------------------------
3161! level 3 souroutines
3162!--------------------------------------------------------
3163 subroutine cudtdqn(klon,klev,ktopm2,kctop,kdtop,ldcum, &
3164 lddraf,ztmst,paph,pgeoh,pgeo,pten,ptenh,pqen, &
3165 pqenh,pqsen,plglac,plude,pmfu,pmfd,pmfus,pmfds, &
3166 pmfuq,pmfdq,pmful,pdmfup,pdmfdp,pdpmel,ptent,ptenq,pcte)
3167 implicit none
3168 integer klon,klev,ktopm2
3169 integer kctop(klon), kdtop(klon)
3170 logical ldcum(klon), lddraf(klon)
3171 real(kind=kind_phys) ztmst
3172 real(kind=kind_phys) paph(klon,klev+1), pgeoh(klon,klev+1)
3173 real(kind=kind_phys) pgeo(klon,klev), pten(klon,klev), &
3174 pqen(klon,klev), ptenh(klon,klev),&
3175 pqenh(klon,klev), pqsen(klon,klev),&
3176 plglac(klon,klev), plude(klon,klev)
3177 real(kind=kind_phys) pmfu(klon,klev), pmfd(klon,klev),&
3178 pmfus(klon,klev), pmfds(klon,klev),&
3179 pmfuq(klon,klev), pmfdq(klon,klev),&
3180 pmful(klon,klev), pdmfup(klon,klev),&
3181 pdpmel(klon,klev), pdmfdp(klon,klev)
3182 real(kind=kind_phys) ptent(klon,klev), ptenq(klon,klev)
3183 real(kind=kind_phys) pcte(klon,klev)
3184
3185! local variables
3186 integer jk , ik , jl
3187 real(kind=kind_phys) zalv , zzp
3188 real(kind=kind_phys) zdtdt(klon,klev) , zdqdt(klon,klev) , zdp(klon,klev)
3189 !* 1.0 SETUP AND INITIALIZATIONS
3190 ! -------------------------
3191 do jk = 1 , klev
3192 do jl = 1, klon
3193 if ( ldcum(jl) ) then
3194 zdp(jl,jk) = g/(paph(jl,jk+1)-paph(jl,jk))
3195 end if
3196 end do
3197 end do
3198
3199 !-----------------------------------------------------------------------
3200 !* 2.0 COMPUTE TENDENCIES
3201 ! ------------------
3202 do jk = ktopm2 , klev
3203 if ( jk < klev ) then
3204 do jl = 1,klon
3205 if ( ldcum(jl) ) then
3206 zalv = foelhm(pten(jl,jk))
3207 zdtdt(jl,jk) = zdp(jl,jk)*rcpd * &
3208 (pmfus(jl,jk+1)-pmfus(jl,jk)+pmfds(jl,jk+1) - &
3209 pmfds(jl,jk)+alf*plglac(jl,jk)-alf*pdpmel(jl,jk) - &
3210 zalv*(pmful(jl,jk+1)-pmful(jl,jk)-plude(jl,jk)-pdmfup(jl,jk)-pdmfdp(jl,jk)))
3211 zdqdt(jl,jk) = zdp(jl,jk)*(pmfuq(jl,jk+1) - &
3212 pmfuq(jl,jk)+pmfdq(jl,jk+1)-pmfdq(jl,jk)+pmful(jl,jk+1) - &
3213 pmful(jl,jk)-plude(jl,jk)-pdmfup(jl,jk)-pdmfdp(jl,jk))
3214 end if
3215 end do
3216 else
3217 do jl = 1,klon
3218 if ( ldcum(jl) ) then
3219 zalv = foelhm(pten(jl,jk))
3220 zdtdt(jl,jk) = -zdp(jl,jk)*rcpd * &
3221 (pmfus(jl,jk)+pmfds(jl,jk)+alf*pdpmel(jl,jk) - &
3222 zalv*(pmful(jl,jk)+pdmfup(jl,jk)+pdmfdp(jl,jk)+plude(jl,jk)))
3223 zdqdt(jl,jk) = -zdp(jl,jk)*(pmfuq(jl,jk) + plude(jl,jk) + &
3224 pmfdq(jl,jk)+(pmful(jl,jk)+pdmfup(jl,jk)+pdmfdp(jl,jk)))
3225 end if
3226 end do
3227 end if
3228 end do
3229 !---------------------------------------------------------------
3230 !* 3.0 UPDATE TENDENCIES
3231 ! -----------------
3232 do jk = ktopm2 , klev
3233 do jl = 1, klon
3234 if ( ldcum(jl) ) then
3235 ptent(jl,jk) = ptent(jl,jk) + zdtdt(jl,jk)
3236 ptenq(jl,jk) = ptenq(jl,jk) + zdqdt(jl,jk)
3237 pcte(jl,jk) = zdp(jl,jk)*plude(jl,jk)
3238 end if
3239 end do
3240 end do
3241
3242 return
3243 end subroutine cudtdqn
3244!---------------------------------------------------------
3245! level 3 souroutines
3246!--------------------------------------------------------
3247 subroutine cududvn(klon,klev,ktopm2,ktype,kcbot,kctop,ldcum, &
3248 ztmst,paph,puen,pven,pmfu,pmfd,puu,pud,pvu,pvd,ptenu, &
3249 ptenv)
3250 implicit none
3251 integer klon,klev,ktopm2
3252 integer ktype(klon), kcbot(klon), kctop(klon)
3253 logical ldcum(klon)
3254 real(kind=kind_phys) ztmst
3255 real(kind=kind_phys) paph(klon,klev+1)
3256 real(kind=kind_phys) puen(klon,klev), pven(klon,klev),&
3257 pmfu(klon,klev), pmfd(klon,klev),&
3258 puu(klon,klev), pud(klon,klev),&
3259 pvu(klon,klev), pvd(klon,klev)
3260 real(kind=kind_phys) ptenu(klon,klev), ptenv(klon,klev)
3261
3262!local variables
3263 real(kind=kind_phys) zuen(klon,klev) , zven(klon,klev) , zmfuu(klon,klev), &
3264 zmfdu(klon,klev), zmfuv(klon,klev), zmfdv(klon,klev)
3265
3266 integer ik , ikb , jk , jl
3267 real(kind=kind_phys) zzp, zdtdt
3268
3269 real(kind=kind_phys) zdudt(klon,klev), zdvdt(klon,klev), zdp(klon,klev)
3270!
3271 do jk = 1 , klev
3272 do jl = 1, klon
3273 if ( ldcum(jl) ) then
3274 zuen(jl,jk) = puen(jl,jk)
3275 zven(jl,jk) = pven(jl,jk)
3276 zdp(jl,jk) = g/(paph(jl,jk+1)-paph(jl,jk))
3277 end if
3278 end do
3279 end do
3280!----------------------------------------------------------------------
3281!* 1.0 CALCULATE FLUXES AND UPDATE U AND V TENDENCIES
3282! ----------------------------------------------
3283 do jk = ktopm2 , klev
3284 ik = jk - 1
3285 do jl = 1,klon
3286 if ( ldcum(jl) ) then
3287 zmfuu(jl,jk) = pmfu(jl,jk)*(puu(jl,jk)-zuen(jl,ik))
3288 zmfuv(jl,jk) = pmfu(jl,jk)*(pvu(jl,jk)-zven(jl,ik))
3289 zmfdu(jl,jk) = pmfd(jl,jk)*(pud(jl,jk)-zuen(jl,ik))
3290 zmfdv(jl,jk) = pmfd(jl,jk)*(pvd(jl,jk)-zven(jl,ik))
3291 end if
3292 end do
3293 end do
3294 ! linear fluxes below cloud
3295 do jk = ktopm2 , klev
3296 do jl = 1, klon
3297 if ( ldcum(jl) .and. jk > kcbot(jl) ) then
3298 ikb = kcbot(jl)
3299 zzp = ((paph(jl,klev+1)-paph(jl,jk))/(paph(jl,klev+1)-paph(jl,ikb)))
3300 if ( ktype(jl) == 3 ) zzp = zzp*zzp
3301 zmfuu(jl,jk) = zmfuu(jl,ikb)*zzp
3302 zmfuv(jl,jk) = zmfuv(jl,ikb)*zzp
3303 zmfdu(jl,jk) = zmfdu(jl,ikb)*zzp
3304 zmfdv(jl,jk) = zmfdv(jl,ikb)*zzp
3305 end if
3306 end do
3307 end do
3308!----------------------------------------------------------------------
3309!* 2.0 COMPUTE TENDENCIES
3310! ------------------
3311 do jk = ktopm2 , klev
3312 if ( jk < klev ) then
3313 ik = jk + 1
3314 do jl = 1,klon
3315 if ( ldcum(jl) ) then
3316 zdudt(jl,jk) = zdp(jl,jk) * &
3317 (zmfuu(jl,ik)-zmfuu(jl,jk)+zmfdu(jl,ik)-zmfdu(jl,jk))
3318 zdvdt(jl,jk) = zdp(jl,jk) * &
3319 (zmfuv(jl,ik)-zmfuv(jl,jk)+zmfdv(jl,ik)-zmfdv(jl,jk))
3320 end if
3321 end do
3322 else
3323 do jl = 1,klon
3324 if ( ldcum(jl) ) then
3325 zdudt(jl,jk) = -zdp(jl,jk)*(zmfuu(jl,jk)+zmfdu(jl,jk))
3326 zdvdt(jl,jk) = -zdp(jl,jk)*(zmfuv(jl,jk)+zmfdv(jl,jk))
3327 end if
3328 end do
3329 end if
3330 end do
3331!---------------------------------------------------------------------
3332!* 3.0 UPDATE TENDENCIES
3333! -----------------
3334 do jk = ktopm2 , klev
3335 do jl = 1, klon
3336 if ( ldcum(jl) ) then
3337 ptenu(jl,jk) = ptenu(jl,jk) + zdudt(jl,jk)
3338 ptenv(jl,jk) = ptenv(jl,jk) + zdvdt(jl,jk)
3339 end if
3340 end do
3341 end do
3342!----------------------------------------------------------------------
3343 return
3344 end subroutine cududvn
3345
3346!---------------------------------------------------------
3347! level 3 souroutines
3348!--------------------------------------------------------
3349 subroutine cuctracer(klon,klev,ktrac,kctop,kdtop, &
3350 ldcum,lddraf,ztmst,paph,pmfu,pmfd, &
3351 pudrate,pddrate,pcen,ptenc)
3352 implicit none
3353 integer klon,klev,ktrac
3354 integer kctop(klon), kdtop(klon)
3355 logical ldcum(klon), lddraf(klon)
3356 real(kind=kind_phys) ztmst
3357 real(kind=kind_phys) paph(klon,klev+1)
3358 real(kind=kind_phys) pmfu(klon,klev)
3359 real(kind=kind_phys) pmfd(klon,klev)
3360 real(kind=kind_phys) pudrate(klon,klev)
3361 real(kind=kind_phys) pddrate(klon,klev)
3362 real(kind=kind_phys) pcen(klon,klev,ktrac)
3363 real(kind=kind_phys) ptenc(klon,klev,ktrac)
3364 !----------------------------------------------------------------------
3365 integer ik , jk , jl , jn
3366 real(kind=kind_phys) zzp , zmfa , zerate , zposi
3367 ! ALLOCATABLE ARAYS
3368 real(kind=kind_phys) , dimension(:,:,:) , allocatable :: zcen , zcu , zcd , &
3369 ztenc , zmfc
3370 real(kind=kind_phys) , dimension(:,:) , allocatable :: zdp
3371 logical , dimension(:,:) , allocatable :: llcumask , llcumbas
3372 !----------------------------------------------------------------------
3373 allocate (zcen(klon,klev,ktrac)) ! Half-level environmental values
3374 allocate (zcu(klon,klev,ktrac)) ! Updraft values
3375 allocate (zcd(klon,klev,ktrac)) ! Downdraft values
3376 allocate (ztenc(klon,klev,ktrac)) ! Tendency
3377 allocate (zmfc(klon,klev,ktrac)) ! Fluxes
3378 allocate (zdp(klon,klev)) ! Pressure difference
3379 allocate (llcumask(klon,klev)) ! Mask for convection
3380 ! Initialize Cumulus mask + some setups
3381 do jk = 2, klev
3382 do jl = 1, klon
3383 llcumask(jl,jk) = .false.
3384 if ( ldcum(jl) ) then
3385 zdp(jl,jk) = g/(paph(jl,jk+1)-paph(jl,jk))
3386 if ( jk >= kctop(jl)-1 ) llcumask(jl,jk) = .true.
3387 end if
3388 end do
3389 end do
3390 !----------------------------------------------------------------------
3391 do jn = 1 , ktrac
3392 !* 1.0 DEFINE TRACERS AT HALF LEVELS
3393 ! -----------------------------
3394 do jk = 2 , klev
3395 ik = jk - 1
3396 do jl = 1, klon
3397 zcen(jl,jk,jn) = pcen(jl,jk,jn)
3398 zcd(jl,jk,jn) = pcen(jl,ik,jn)
3399 zcu(jl,jk,jn) = pcen(jl,ik,jn)
3400 zmfc(jl,jk,jn) = 0.
3401 ztenc(jl,jk,jn)= 0.
3402 end do
3403 end do
3404
3405 do jl = 1, klon
3406 zcu(jl,klev,jn) = pcen(jl,klev,jn)
3407 end do
3408 !* 2.0 COMPUTE UPDRAFT VALUES
3409 ! ----------------------
3410 do jk = klev - 1 , 3 , -1
3411 ik = jk + 1
3412 do jl = 1, klon
3413 if ( llcumask(jl,jk) ) then
3414 zerate = pmfu(jl,jk) - pmfu(jl,ik) + pudrate(jl,jk)
3415 zmfa = 1./max(cmfcmin,pmfu(jl,jk))
3416 if ( jk >= kctop(jl) ) then
3417 zcu(jl,jk,jn) = (pmfu(jl,ik)*zcu(jl,ik,jn) + &
3418 zerate*pcen(jl,jk,jn)-pudrate(jl,jk)*zcu(jl,ik,jn))*zmfa
3419 end if
3420 end if
3421 end do
3422 end do
3423 !* 3.0 COMPUTE DOWNDRAFT VALUES
3424 ! ------------------------
3425 do jk = 3 , klev
3426 ik = jk - 1
3427 do jl = 1, klon
3428 if ( lddraf(jl) .and. jk == kdtop(jl) ) then
3429 ! Nota: in order to avoid final negative Tracer values at LFS
3430 ! the allowed value of ZCD depends on the jump in mass flux
3431 ! at the LFS
3432 zcd(jl,jk,jn) = 0.1*zcu(jl,jk,jn) + 0.9*pcen(jl,ik,jn)
3433 else if ( lddraf(jl).and.jk>kdtop(jl) ) then
3434 zerate = -pmfd(jl,jk) + pmfd(jl,ik) + pddrate(jl,jk)
3435 zmfa = 1./min(-cmfcmin,pmfd(jl,jk))
3436 zcd(jl,jk,jn) = (pmfd(jl,ik)*zcd(jl,ik,jn) - &
3437 zerate*pcen(jl,ik,jn)+pddrate(jl,jk)*zcd(jl,ik,jn))*zmfa
3438 end if
3439 end do
3440 end do
3441 ! In order to avoid negative Tracer at KLEV adjust ZCD
3442 jk = klev
3443 ik = jk - 1
3444 do jl = 1, klon
3445 if ( lddraf(jl) ) then
3446 zposi = -zdp(jl,jk) *(pmfu(jl,jk)*zcu(jl,jk,jn) + &
3447 pmfd(jl,jk)*zcd(jl,jk,jn)-(pmfu(jl,jk)+pmfd(jl,jk))*pcen(jl,ik,jn))
3448 if ( pcen(jl,jk,jn)+zposi*ztmst < 0. ) then
3449 zmfa = 1./min(-cmfcmin,pmfd(jl,jk))
3450 zcd(jl,jk,jn) = ((pmfu(jl,jk)+pmfd(jl,jk))*pcen(jl,ik,jn) - &
3451 pmfu(jl,jk)*zcu(jl,jk,jn)+pcen(jl,jk,jn) / &
3452 (ztmst*zdp(jl,jk)))*zmfa
3453 end if
3454 end if
3455 end do
3456 end do
3457 !----------------------------------------------------------------------
3458 do jn = 1 , ktrac
3459 !* 4.0 COMPUTE FLUXES
3460 ! --------------
3461 do jk = 2 , klev
3462 ik = jk - 1
3463 do jl = 1, klon
3464 if ( llcumask(jl,jk) ) then
3465 zmfa = pmfu(jl,jk) + pmfd(jl,jk)
3466 zmfc(jl,jk,jn) = pmfu(jl,jk)*zcu(jl,jk,jn) + &
3467 pmfd(jl,jk)*zcd(jl,jk,jn) - zmfa*zcen(jl,ik,jn)
3468 end if
3469 end do
3470 end do
3471 !* 5.0 COMPUTE TENDENCIES = RHS
3472 ! ------------------------
3473 do jk = 2 , klev - 1
3474 ik = jk + 1
3475 do jl = 1, klon
3476 if ( llcumask(jl,jk) ) then
3477 ztenc(jl,jk,jn) = zdp(jl,jk)*(zmfc(jl,ik,jn)-zmfc(jl,jk,jn))
3478 end if
3479 end do
3480 end do
3481 jk = klev
3482 do jl = 1, klon
3483 if ( ldcum(jl) ) ztenc(jl,jk,jn) = -zdp(jl,jk)*zmfc(jl,jk,jn)
3484 end do
3485 end do
3486 !* 6.0 UPDATE TENDENCIES
3487 ! -----------------
3488 do jn = 1, ktrac
3489 do jk = 2, klev
3490 do jl = 1, klon
3491 if ( llcumask(jl,jk) ) then
3492 ptenc(jl,jk,jn) = ptenc(jl,jk,jn)+ztenc(jl,jk,jn)
3493 end if
3494 end do
3495 end do
3496 end do
3497 !---------------------------------------------------------------------------
3498 deallocate (llcumask)
3499 deallocate (zdp)
3500 deallocate (zmfc)
3501 deallocate (ztenc)
3502 deallocate (zcd)
3503 deallocate (zcu)
3504 deallocate (zcen)
3505 end subroutine cuctracer
3506
3507!---------------------------------------------------------
3508! level 4 souroutines
3509!--------------------------------------------------------
3510 subroutine cuadjtqn &
3511 & (klon, klev, kk, psp, pt, pq, ldflag, kcall)
3512! m.tiedtke e.c.m.w.f. 12/89
3513! purpose.
3514! --------
3515! to produce t,q and l values for cloud ascent
3516
3517! interface
3518! ---------
3519! this routine is called from subroutines:
3520! *cond* (t and q at condensation level)
3521! *cubase* (t and q at condensation level)
3522! *cuasc* (t and q at cloud levels)
3523! *cuini* (environmental t and qs values at half levels)
3524! input are unadjusted t and q values,
3525! it returns adjusted values of t and q
3526
3527! parameter description units
3528! --------- ----------- -----
3529! input parameters (integer):
3530
3531! *klon* number of grid points per packet
3532! *klev* number of levels
3533! *kk* level
3534! *kcall* defines calculation as
3535! kcall=0 env. t and qs in*cuini*
3536! kcall=1 condensation in updrafts (e.g. cubase, cuasc)
3537! kcall=2 evaporation in downdrafts (e.g. cudlfs,cuddraf)
3538! input parameters (real(kind=kind_phys)):
3539
3540! *psp* pressure pa
3541
3542! updated parameters (real(kind=kind_phys)):
3543
3544! *pt* temperature k
3545! *pq* specific humidity kg/kg
3546! externals
3547! ---------
3548! for condensation calculations.
3549! the tables are initialised in *suphec*.
3550
3551!----------------------------------------------------------------------
3552
3553 implicit none
3554
3555 integer klev,klon
3556 real(kind=kind_phys) pt(klon,klev), pq(klon,klev), &
3557 & psp(klon)
3558 logical ldflag(klon)
3559! local variables
3560 integer jl,jk
3561 integer isum,kcall,kk
3562 real(kind=kind_phys) zqmax,zqsat,zcor,zqp,zcond,zcond1,zl,zi,zf
3563!----------------------------------------------------------------------
3564! 1. define constants
3565! ----------------
3566 zqmax=0.5
3567
3568! 2. calculate condensation and adjust t and q accordingly
3569! -----------------------------------------------------
3570
3571 if ( kcall == 1 ) then
3572 do jl = 1,klon
3573 if ( ldflag(jl) ) then
3574 zqp = 1./psp(jl)
3575 zl = 1./(pt(jl,kk)-c4les)
3576 zi = 1./(pt(jl,kk)-c4ies)
3577 zqsat = c2es*(foealfa(pt(jl,kk))*exp(c3les*(pt(jl,kk)-tmelt)*zl) + &
3578 (1.-foealfa(pt(jl,kk)))*exp(c3ies*(pt(jl,kk)-tmelt)*zi))
3579 zqsat = zqsat*zqp
3580 zqsat = min(0.5,zqsat)
3581 zcor = 1. - vtmpc1*zqsat
3582 zf = foealfa(pt(jl,kk))*r5alvcp*zl**2 + &
3583 (1.-foealfa(pt(jl,kk)))*r5alscp*zi**2
3584 zcond = (pq(jl,kk)*zcor**2-zqsat*zcor)/(zcor**2+zqsat*zf)
3585 if ( zcond > 0. ) then
3586 pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond
3587 pq(jl,kk) = pq(jl,kk) - zcond
3588 zl = 1./(pt(jl,kk)-c4les)
3589 zi = 1./(pt(jl,kk)-c4ies)
3590 zqsat = c2es*(foealfa(pt(jl,kk)) * &
3591 exp(c3les*(pt(jl,kk)-tmelt)*zl)+(1.-foealfa(pt(jl,kk))) * &
3592 exp(c3ies*(pt(jl,kk)-tmelt)*zi))
3593 zqsat = zqsat*zqp
3594 zqsat = min(0.5,zqsat)
3595 zcor = 1. - vtmpc1*zqsat
3596 zf = foealfa(pt(jl,kk))*r5alvcp*zl**2 + &
3597 (1.-foealfa(pt(jl,kk)))*r5alscp*zi**2
3598 zcond1 = (pq(jl,kk)*zcor**2-zqsat*zcor)/(zcor**2+zqsat*zf)
3599 if ( abs(zcond) < 1.e-20 ) zcond1 = 0.
3600 pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1
3601 pq(jl,kk) = pq(jl,kk) - zcond1
3602 end if
3603 end if
3604 end do
3605 elseif ( kcall == 2 ) then
3606 do jl = 1,klon
3607 if ( ldflag(jl) ) then
3608 zqp = 1./psp(jl)
3609 zqsat = foeewm(pt(jl,kk))*zqp
3610 zqsat = min(0.5,zqsat)
3611 zcor = 1./(1.-vtmpc1*zqsat)
3612 zqsat = zqsat*zcor
3613 zcond = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk)))
3614 zcond = min(zcond,0.)
3615 pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond
3616 pq(jl,kk) = pq(jl,kk) - zcond
3617 zqsat = foeewm(pt(jl,kk))*zqp
3618 zqsat = min(0.5,zqsat)
3619 zcor = 1./(1.-vtmpc1*zqsat)
3620 zqsat = zqsat*zcor
3621 zcond1 = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk)))
3622 if ( abs(zcond) < 1.e-20 ) zcond1 = min(zcond1,0.)
3623 pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1
3624 pq(jl,kk) = pq(jl,kk) - zcond1
3625 end if
3626 end do
3627 else if ( kcall == 0 ) then
3628 do jl = 1,klon
3629 zqp = 1./psp(jl)
3630 zqsat = foeewm(pt(jl,kk))*zqp
3631 zqsat = min(0.5,zqsat)
3632 zcor = 1./(1.-vtmpc1*zqsat)
3633 zqsat = zqsat*zcor
3634 zcond1 = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk)))
3635 pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1
3636 pq(jl,kk) = pq(jl,kk) - zcond1
3637 zqsat = foeewm(pt(jl,kk))*zqp
3638 zqsat = min(0.5,zqsat)
3639 zcor = 1./(1.-vtmpc1*zqsat)
3640 zqsat = zqsat*zcor
3641 zcond1 = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk)))
3642 pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1
3643 pq(jl,kk) = pq(jl,kk) - zcond1
3644 end do
3645 end if
3646
3647 return
3648 end subroutine cuadjtqn
3649!---------------------------------------------------------
3650! level 4 souroutines
3651!--------------------------------------------------------
3652 subroutine cubasmcn &
3653 & (klon, klev, klevm1, kk, pten,&
3654 & pqen, pqsen, puen, pven, pverv,&
3655 & pgeo, pgeoh, ldcum, ktype, klab, plrain,&
3656 & pmfu, pmfub, kcbot, ptu,&
3657 & pqu, plu, puu, pvu, pmfus,&
3658 & pmfuq, pmful, pdmfup)
3659 implicit none
3660! m.tiedtke e.c.m.w.f. 12/89
3661! c.zhang iprc 05/2012
3662!***purpose.
3663! --------
3664! this routine calculates cloud base values
3665! for midlevel convection
3666!***interface
3667! ---------
3668! this routine is called from *cuasc*.
3669! input are environmental values t,q etc
3670! it returns cloudbase values for midlevel convection
3671!***method.
3672! -------
3673! s. tiedtke (1989)
3674!***externals
3675! ---------
3676! none
3677! ----------------------------------------------------------------
3678 real(kind=kind_phys) pten(klon,klev), pqen(klon,klev),&
3679 & puen(klon,klev), pven(klon,klev),&
3680 & pqsen(klon,klev), pverv(klon,klev),&
3681 & pgeo(klon,klev), pgeoh(klon,klev+1)
3682 real(kind=kind_phys) ptu(klon,klev), pqu(klon,klev),&
3683 & puu(klon,klev), pvu(klon,klev),&
3684 & plu(klon,klev), pmfu(klon,klev),&
3685 & pmfub(klon), &
3686 & pmfus(klon,klev), pmfuq(klon,klev),&
3687 & pmful(klon,klev), pdmfup(klon,klev),&
3688 & plrain(klon,klev)
3689 integer ktype(klon), kcbot(klon),&
3690 & klab(klon,klev)
3691 logical ldcum(klon)
3692! local variabels
3693 integer jl,kk,klev,klon,klevp1,klevm1
3694 real(kind=kind_phys) zzzmb
3695!--------------------------------------------------------
3696!* 1. calculate entrainment and detrainment rates
3697! -------------------------------------------------------
3698 do jl=1,klon
3699 if(.not.ldcum(jl) .and. klab(jl,kk+1).eq.0) then
3700 if(lmfmid .and. pqen(jl,kk) .gt. 0.80*pqsen(jl,kk).and. &
3701 pgeo(jl,kk)*zrg .gt. 5.0e2 .and. &
3702 & pgeo(jl,kk)*zrg .lt. 1.0e4 ) then
3703 ptu(jl,kk+1)=(cpd*pten(jl,kk)+pgeo(jl,kk)-pgeoh(jl,kk+1))&
3704 & *rcpd
3705 pqu(jl,kk+1)=pqen(jl,kk)
3706 plu(jl,kk+1)=0.
3707 zzzmb=max(cmfcmin,-pverv(jl,kk)*zrg)
3708 zzzmb=min(zzzmb,cmfcmax)
3709 pmfub(jl)=zzzmb
3710 pmfu(jl,kk+1)=pmfub(jl)
3711 pmfus(jl,kk+1)=pmfub(jl)*(cpd*ptu(jl,kk+1)+pgeoh(jl,kk+1))
3712 pmfuq(jl,kk+1)=pmfub(jl)*pqu(jl,kk+1)
3713 pmful(jl,kk+1)=0.
3714 pdmfup(jl,kk+1)=0.
3715 kcbot(jl)=kk
3716 klab(jl,kk+1)=1
3717 plrain(jl,kk+1)=0.0
3718 ktype(jl)=3
3719 end if
3720 end if
3721 end do
3722 return
3723 end subroutine cubasmcn
3724!
3725!---------------------------------------------------------
3726! level 4 souroutines
3727!---------------------------------------------------------
3728 subroutine cuentrn(klon,klev,kk,kcbot,ldcum,ldwork, &
3729 pgeoh,pmfu,pdmfen,pdmfde)
3730 implicit none
3731 integer klon,klev,kk
3732 integer kcbot(klon)
3733 logical ldcum(klon)
3734 logical ldwork
3735 real(kind=kind_phys) pgeoh(klon,klev+1)
3736 real(kind=kind_phys) pmfu(klon,klev)
3737 real(kind=kind_phys) pdmfen(klon)
3738 real(kind=kind_phys) pdmfde(klon)
3739 logical llo1
3740 integer jl
3741 real(kind=kind_phys) zdz , zmf
3742 real(kind=kind_phys) zentr(klon)
3743 !
3744 !* 1. CALCULATE ENTRAINMENT AND DETRAINMENT RATES
3745 ! -------------------------------------------
3746 if ( ldwork ) then
3747 do jl = 1,klon
3748 pdmfen(jl) = 0.
3749 pdmfde(jl) = 0.
3750 zentr(jl) = 0.
3751 end do
3752 !
3753 !* 1.1 SPECIFY ENTRAINMENT RATES
3754 ! -------------------------
3755 do jl = 1, klon
3756 if ( ldcum(jl) ) then
3757 zdz = (pgeoh(jl,kk)-pgeoh(jl,kk+1))*zrg
3758 zmf = pmfu(jl,kk+1)*zdz
3759 llo1 = kk < kcbot(jl)
3760 if ( llo1 ) then
3761 pdmfen(jl) = zentr(jl)*zmf
3762 pdmfde(jl) = 0.75e-4*zmf
3763 end if
3764 end if
3765 end do
3766 end if
3767 end subroutine cuentrn
3768!
3769!--------------------------------------------------------
3770! external functions
3771!------------------------------------------------------
3772 real(kind=kind_phys) function foealfa(tt)
3773! foealfa is calculated to distinguish the three cases:
3774!
3775! foealfa=1 water phase
3776! foealfa=0 ice phase
3777! 0 < foealfa < 1 mixed phase
3778!
3779! input : tt = temperature
3780!
3781 implicit none
3782 real(kind=kind_phys) tt
3783 foealfa = min(1.,((max(rtice,min(rtwat,tt))-rtice) &
3784 & /(rtwat-rtice))**2)
3785
3786 return
3787 end function foealfa
3788
3789 real(kind=kind_phys) function foelhm(tt)
3790 implicit none
3791 real(kind=kind_phys) tt
3792 foelhm = foealfa(tt)*alv + (1.-foealfa(tt))*als
3793 return
3794 end function foelhm
3795
3796 real(kind=kind_phys) function foeewm(tt)
3797 implicit none
3798 real(kind=kind_phys) tt
3799 foeewm = c2es * &
3800 & (foealfa(tt)*exp(c3les*(tt-tmelt)/(tt-c4les))+ &
3801 & (1.-foealfa(tt))*exp(c3ies*(tt-tmelt)/(tt-c4ies)))
3802 return
3803 end function foeewm
3804
3805 real(kind=kind_phys) function foedem(tt)
3806 implicit none
3807 real(kind=kind_phys) tt
3808 foedem = foealfa(tt)*r5alvcp*(1./(tt-c4les)**2)+ &
3809 & (1.-foealfa(tt))*r5alscp*(1./(tt-c4ies)**2)
3810 return
3811 end function foedem
3812
3813 real(kind=kind_phys) function foeldcpm(tt)
3814 implicit none
3815 real(kind=kind_phys) tt
3816 foeldcpm = foealfa(tt)*ralvdcp+ &
3817 & (1.-foealfa(tt))*ralsdcp
3818 return
3819 end function foeldcpm
3820
3821 real(kind=kind_phys) function foeldcp(tt)
3822 implicit none
3823 real(kind=kind_phys) tt
3824 foeldcp = foedelta(tt)*ralvdcp + (1.-foedelta(tt))*ralsdcp
3825 end function foeldcp
3826
3827 real(kind=kind_phys) function foedelta(tt)
3828 implicit none
3829 real(kind=kind_phys) tt
3830 foedelta = max(0.,sign(1.,tt-tmelt))
3831 end function foedelta
3832
3833end module cu_ntiedtke
3834
This module contains some of the most frequently used math and physics constants for GCM models.
Definition physcons.F90:36