CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
cu_c3_sh.F90
1
3
4module cu_c3_sh
5 use machine , only : kind_phys
6 use progsigma, only : progsigma_calc
7
8 !real(kind=kind_phys), parameter:: c1_shal=0.0015! .0005
9 real(kind=kind_phys), parameter:: g =9.81
10 real(kind=kind_phys), parameter:: cp =1004.
11 real(kind=kind_phys), parameter:: xlv=2.5e6
12 real(kind=kind_phys), parameter:: r_v=461.
13 real(kind=kind_phys) :: c0_shal=.004
14 real(kind=kind_phys) :: c1_shal=0. !0.005! .0005
15 real(kind=kind_phys), parameter:: fluxtune=1.5
16
17contains
18
66 subroutine cu_c3_sh_run ( &
67 us,vs,zo,t,q,z1,tn,qo,po,psur,dhdt,kpbl,rho, & ! input variables, must be supplied
68 hfx,qfx,xland,ichoice,tcrit,dtime, &
69 zuo,xmb_out,kbcon,ktop,k22,ierr,ierrc, &
70 flag_init, flag_restart,fv,r_d,delp,tmf,qmicro, &
71 forceqv_spechum,betascu,betamcu,betadcu,sigmain,&
72 sigmaout,progsigma,dx, &
73 outt,outq,outqc,outu,outv,cnvwt,pre,cupclw, & ! output tendencies
74 itf,ktf,its,ite, kts,kte,ipr,tropics) ! dimesnional variables
75!
76! this module needs some subroutines from gf_deep
77!
82
83 implicit none
84 integer &
85 ,intent (in ) :: &
86 itf,ktf, &
87 its,ite, kts,kte,ipr
88 logical, intent(in) :: flag_init, flag_restart, progsigma
89 logical :: make_calc_for_xk = .true.
90 integer, intent (in ) :: &
91 ichoice
92 !
93 !
94 !
95 ! outtem = output temp tendency (per s)
96 ! outq = output q tendency (per s)
97 ! outqc = output qc tendency (per s)
98 ! pre = output precip
99 real(kind=kind_phys), dimension (its:,kts:) &
100 ,intent (inout ) :: &
101 cnvwt,outt,outq,outqc,cupclw,zuo,outu,outv
102!$acc declare copy(cnvwt,outt,outq,outqc,cupclw,zuo,outu,outv)
103 real(kind=kind_phys), dimension (its:,kts:) &
104 ,intent (in ) :: &
105 tmf
106 real(kind=kind_phys), dimension (its:,kts:) &
107 ,intent (in ), optional :: &
108 qmicro, sigmain, forceqv_spechum
109
110 real(kind=kind_phys), dimension (its:) &
111 ,intent (out ) :: &
112 xmb_out
113 integer, dimension (its:) &
114 ,intent (inout ) :: &
115 ierr
116 integer, dimension (its:) &
117 ,intent (out ) :: &
118 kbcon,ktop,k22
119 integer, dimension (its:) &
120 ,intent (in ) :: &
121 kpbl,tropics
122!$acc declare copyout(xmb_out,kbcon,ktop,k22) copyin(kpbl,tropics) copy(ierr)
123 !
124 ! basic environmental input includes a flag (ierr) to turn off
125 ! convection for this call only and at that particular gridpoint
126 !
127 real(kind=kind_phys), dimension (its:,kts:) &
128 ,intent (in ) :: &
129 t,po,tn,dhdt,rho,us,vs,delp
130 real(kind=kind_phys), dimension (its:,kts:) &
131 ,intent (inout) :: &
132 q,qo
133 real(kind=kind_phys), dimension (its:) &
134 ,intent (in ) :: &
135 xland,z1,psur,hfx,qfx,dx
136
137 real(kind=kind_phys) &
138 ,intent (in ) :: &
139 dtime,tcrit,fv,r_d,betascu,betamcu,betadcu
140!$acc declare sigmaout
141 real(kind=kind_phys), dimension (its:,kts:) &
142 ,intent (out), optional :: &
143 sigmaout
144
145
146!$acc declare copyin(t,po,tn,dhdt,rho,us,vs) copy(q,qo) copyin(xland,z1,psur,hfx,qfx) copyin(dtime,tcrit)
147 !
148 !***************** the following are your basic environmental
149 ! variables. they carry a "_cup" if they are
150 ! on model cloud levels (staggered). they carry
151 ! an "o"-ending (z becomes zo), if they are the forced
152 ! variables.
153 !
154 ! z = heights of model levels
155 ! q = environmental mixing ratio
156 ! qes = environmental saturation mixing ratio
157 ! t = environmental temp
158 ! p = environmental pressure
159 ! he = environmental moist static energy
160 ! hes = environmental saturation moist static energy
161 ! z_cup = heights of model cloud levels
162 ! q_cup = environmental q on model cloud levels
163 ! qes_cup = saturation q on model cloud levels
164 ! t_cup = temperature (kelvin) on model cloud levels
165 ! p_cup = environmental pressure
166 ! he_cup = moist static energy on model cloud levels
167 ! hes_cup = saturation moist static energy on model cloud levels
168 ! gamma_cup = gamma on model cloud levels
169 ! dby = buoancy term
170 ! entr = entrainment rate
171 ! bu = buoancy term
172 ! gamma_cup = gamma on model cloud levels
173 ! qrch = saturation q in cloud
174 ! pwev = total normalized integrated evaoprate (i2)
175 ! z1 = terrain elevation
176 ! psur = surface pressure
177 ! zu = updraft normalized mass flux
178 ! kbcon = lfc of parcel from k22
179 ! k22 = updraft originating level
180 ! ichoice = flag if only want one closure (usually set to zero!)
181 ! dby = buoancy term
182 ! ktop = cloud top (output)
183 ! xmb = total base mass flux
184 ! hc = cloud moist static energy
185 ! hkb = moist static energy at originating level
186
187 real(kind=kind_phys), dimension (its:ite,kts:kte) :: &
188 entr_rate_2d,he,hes,qes,z, &
189 heo,heso,qeso,zo, &
190 xhe,xhes,xqes,xz,xt,xq, &
191 qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
192 qeso_cup,qo_cup,heo_cup,heso_cup,zo_cup,po_cup,gammao_cup, &
193 tn_cup, &
194 xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup, &
195 xt_cup,dby,hc,zu, &
196 dbyo,qco,pwo,hco,qrco, &
197 dbyt,xdby,xhc,xzu, &
198
199 ! cd = detrainment function for updraft
200 ! dellat = change of temperature per unit mass flux of cloud ensemble
201 ! dellaq = change of q per unit mass flux of cloud ensemble
202 ! dellaqc = change of qc per unit mass flux of cloud ensemble
203
204 cd,dellah,dellaq,dellat,dellaqc,uc,vc,dellu,dellv,u_cup,v_cup, &
205 wu2,omega_u,zeta,zdqca,del,clw_all
206
207!$acc declare create( &
208!$acc entr_rate_2d,he,hes,qes,z, &
209!$acc heo,heso,qeso,zo, &
210!$acc xhe,xhes,xqes,xz,xt,xq, &
211!$acc qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
212!$acc qeso_cup,qo_cup,heo_cup,heso_cup,zo_cup,po_cup,gammao_cup, &
213!$acc tn_cup, &
214!$acc xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup, &
215!$acc xt_cup,dby,hc,zu, &
216!$acc dbyo,qco,pwo,hco,qrco, &
217!$acc dbyt,xdby,xhc,xzu, &
218!$acc cd,dellah,dellaq,dellat,dellaqc,uc,vc,dellu,dellv,u_cup,v_cup)
219
220 ! aa0 cloud work function for downdraft
221 ! aa0 = cloud work function without forcing effects
222 ! aa1 = cloud work function with forcing effects
223 ! xaa0 = cloud work function with cloud effects (ensemble dependent)
224
225 real(kind=kind_phys), dimension (its:ite) :: &
226 zws,ztexec,zqexec,pre,aa1,aa0,xaa0,hkb, &
227 flux_tun,hkbo,xhkb, &
228 rand_vmas,xmbmax,xmb, &
229 cap_max,entr_rate, &
230 cap_max_increment,lambau,wc,omegac,sigmab, &
231 scaldfunc
232 integer, dimension (its:ite) :: &
233 kstabi,xland1,kbmax,ktopx
234!$acc declare create( &
235!$acc zws,ztexec,zqexec,pre,aa1,aa0,xaa0,hkb, &
236!$acc flux_tun,hkbo,xhkb, &
237!$acc rand_vmas,xmbmax,xmb, &
238!$acc cap_max,entr_rate, &
239!$acc cap_max_increment,lambau, &
240!$acc kstabi,xland1,kbmax,ktopx)
241
242 logical :: flag_shallow,flag_mid
243 logical, dimension(its:ite) :: cnvflg
244 integer :: &
245 kstart,i,k,ki
246 real(kind=kind_phys) :: &
247 dz,mbdt,zkbmax, &
248 cap_maxs,trash,trash2,frh,el2orc,gravinv
249
250 real(kind=kind_phys) buo_flux,pgeoh,dp,entup,detup,totmas
251 real(kind=kind_phys) :: &
252 sigmind,sigminm,sigmins
253 parameter(sigmind=0.005,sigmins=0.03,sigminm=0.01)
254
255 real(kind=kind_phys) xff_shal(3),blqe,xkshal
256 character*50 :: ierrc(its:)
257 real(kind=kind_phys), dimension (its:ite,kts:kte) :: &
258 up_massentr,up_massdetr,up_massentro,up_massdetro,up_massentru,up_massdetru
259!$acc declare create(up_massentr,up_massdetr,up_massentro,up_massdetro,up_massentru,up_massdetru)
260 real(kind=kind_phys) :: c_up,x_add,qaver,dts,fp,fpi
261 real(kind=kind_phys), dimension (its:ite,kts:kte) :: c1d,dtempdz
262 integer, dimension (its:ite,kts:kte) :: k_inv_layers
263 integer, dimension (its:ite) :: start_level, pmin_lev
264!$acc declare create(c1d,dtempdz,k_inv_layers,start_level, pmin_lev)
265
266 real(kind=kind_phys), parameter :: zero = 0
267
268!$acc kernels
269 start_level(:)=0
270 rand_vmas(:)=0.
271 flux_tun(:)=fluxtune
272 lambau(:)=2.
273 c1d(:,:)=0.
274 scaldfunc(:)=0.
275!$acc end kernels
276
277 el2orc=xlv*xlv/(r_v*cp)
278
279!$acc kernels
280 do i=its,itf
281 xland1(i)=int(xland(i)+.001) ! 1.
282 ktopx(i)=0
283 if(xland(i).gt.1.5 .or. xland(i).lt.0.5)then
284 xland1(i)=0
285 c0_shal=.001
286 c1_shal=.001
287! ierr(i)=100
288 endif
289 pre(i)=0.
290 xmb_out(i)=0.
291 cap_max_increment(i)=25.
292 entr_rate(i) = 1.e-3 !9.e-5 ! 1.75e-3 ! 1.2e-3 ! .2/50.
293 enddo
294!$acc end kernels
295
296 do i=its,itf
297 ierrc(i)=" "
298 enddo
299!
300!--- initial entrainment rate (these may be changed later on in the
301!--- program
302!
303
304!
306!
307!$acc kernels
308 do k=kts,ktf
309 do i=its,itf
310 up_massentro(i,k)=0.
311 up_massdetro(i,k)=0.
312 up_massentru(i,k)=0.
313 up_massdetru(i,k)=0.
314 z(i,k)=zo(i,k)
315 xz(i,k)=zo(i,k)
316 qrco(i,k)=0.
317 pwo(i,k)=0.
318 cd(i,k)=.75*entr_rate(i)
319 dellaqc(i,k)=0.
320 cupclw(i,k)=0.
321 enddo
322 enddo
323!$acc end kernels
324!
325!--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
326!
327!--- minimum depth (m), clouds must have
328!
329!
330!--- maximum depth (mb) of capping
331!--- inversion (larger cap = no convection)
332!
333!$acc kernels
334 cap_maxs=175.
335 do i=its,itf
336 kbmax(i)=1
337 aa0(i)=0.
338 aa1(i)=0.
339 enddo
340 do i=its,itf
341 cap_max(i)=cap_maxs
342 ztexec(i) = 0.
343 zqexec(i) = 0.
344 zws(i) = 0.
345 enddo
346 do i=its,itf
347 !- buoyancy flux (h+le)
348 buo_flux= (hfx(i)/cp+0.608*t(i,1)*qfx(i)/xlv)/rho(i,1)
349 pgeoh = zo(i,2)*g
350 !-convective-scale velocity w*
351 zws(i) = max(0.,flux_tun(i)*0.41*buo_flux*zo(i,2)*g/t(i,1))
352 if(zws(i) > tiny(pgeoh)) then
353 !-convective-scale velocity w*
354 zws(i) = 1.2*zws(i)**.3333
355 !- temperature excess
356 ztexec(i) = max(flux_tun(i)*hfx(i)/(rho(i,1)*zws(i)*cp),0.0)
357 !- moisture excess
358 zqexec(i) = max(flux_tun(i)*qfx(i)/xlv/(rho(i,1)*zws(i)),0.)
359 endif
361 !- height of the pbl
362 zws(i) = max(0.,flux_tun(i)*0.41*buo_flux*zo(i,kpbl(i))*g/t(i,kpbl(i)))
363 zws(i) = 1.2*zws(i)**.3333
364 zws(i) = zws(i)*rho(i,kpbl(i)) !check if zrho is correct
365
366 enddo
367!$acc end kernels
368!
370!
371 zkbmax=3000.
372!
374!
375 call cup_env(z,qes,he,hes,t,q,po,z1, &
376 psur,ierr,tcrit,-1, &
377 itf,ktf, &
378 its,ite, kts,kte)
379 call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
380 psur,ierr,tcrit,-1, &
381 itf,ktf, &
382 its,ite, kts,kte)
383
384!
386!
387 call cup_env_clev(t,qes,q,he,hes,z,po,qes_cup,q_cup,he_cup, &
388 hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
389 ierr,z1, &
390 itf,ktf, &
391 its,ite, kts,kte)
392 call cup_env_clev(tn,qeso,qo,heo,heso,zo,po,qeso_cup,qo_cup, &
393 heo_cup,heso_cup,zo_cup,po_cup,gammao_cup,tn_cup,psur, &
394 ierr,z1, &
395 itf,ktf, &
396 its,ite, kts,kte)
397
398!$acc kernels
399 do i=its,itf
400 if(ierr(i).eq.0)then
401 u_cup(i,kts)=us(i,kts)
402 v_cup(i,kts)=vs(i,kts)
403 do k=kts+1,ktf
404 u_cup(i,k)=.5*(us(i,k-1)+us(i,k))
405 v_cup(i,k)=.5*(vs(i,k-1)+vs(i,k))
406 enddo
407 endif
408 enddo
409
410 do i=its,itf
411 if(ierr(i).eq.0)then
412!
413!$acc loop seq
414 do k=kts,ktf
415 if(zo_cup(i,k).gt.zkbmax+z1(i))then
416 kbmax(i)=k
417 go to 25
418 endif
419 enddo
420 25 continue
421!
422 kbmax(i)=min(kbmax(i),ktf/2)
423 endif
424 enddo
425!$acc end kernels
426
427!
428!
429!
431!
432!$acc parallel loop
433 do 36 i=its,itf
434 if(kpbl(i).gt.3)cap_max(i)=po_cup(i,kpbl(i))
435 if(ierr(i) == 0)then
436 k22(i)=maxloc(heo_cup(i,2:kbmax(i)),1)
437 k22(i)=max(2,k22(i))
438 if(k22(i).gt.kbmax(i))then
439 ierr(i)=2
440#ifndef _OPENACC
441 ierrc(i)="could not find k22"
442#endif
443 ktop(i)=0
444 k22(i)=0
445 kbcon(i)=0
446 endif
447 endif
448 36 continue
449!$acc end parallel
450!
453!
454!$acc parallel loop private(x_add)
455 do i=its,itf
456 if(ierr(i).eq.0)then
457 x_add = xlv*zqexec(i)+cp*ztexec(i)
458 call get_cloud_bc(kte,he_cup(i,1:kte),hkb(i),k22(i),x_add)
459 call get_cloud_bc(kte,heo_cup(i,1:kte),hkbo(i),k22(i),x_add)
460 endif ! ierr
461 enddo
462!$acc end parallel
463
464!joe-georg and saulo's new idea:
465
466!$acc kernels
467 do i=its,itf
468 do k=kts,ktf
469 dbyo(i,k)= 0. !hkbo(i)-heso_cup(i,k)
470 clw_all(i,k)=0.
471 enddo
472 enddo
473!$acc end kernels
474
475
476 call cup_kbcon(ierrc,cap_max_increment,5,k22,kbcon,heo_cup,heso_cup, &
477 hkbo,ierr,kbmax,po_cup,cap_max, &
478 ztexec,zqexec, &
479 0,itf,ktf, &
480 its,ite, kts,kte, &
481 z_cup,entr_rate,heo,0)
482
484 call cup_minimi(heso_cup,kbcon,kbmax,kstabi,ierr, &
485 itf,ktf, &
486 its,ite, kts,kte)
487!
488 call get_inversion_layers(ierr,p_cup,t_cup,z_cup,q_cup,qes_cup,k_inv_layers,&
489 kbcon,kstabi,dtempdz,itf,ktf,its,ite, kts,kte)
490!
491!
492!$acc parallel loop private(frh,kstart,x_add)
493 do i=its,itf
494 entr_rate_2d(i,:)=entr_rate(i)
495 if(ierr(i) == 0)then
496 start_level(i)=k22(i)
497 x_add = xlv*zqexec(i)+cp*ztexec(i)
498 call get_cloud_bc(kte,he_cup(i,1:kte),hkb(i),k22(i),x_add)
499 if(kbcon(i).gt.ktf-4)then
500 ierr(i)=231
501 endif
502 do k=kts,ktf
503 frh = 2.*min(qo_cup(i,k)/qeso_cup(i,k),1.)
504 entr_rate_2d(i,k)=entr_rate(i) !*(2.3-frh)
505 cd(i,k)=.75*entr_rate_2d(i,k)
506 enddo
507!
508! first estimate for shallow convection
509!
510 ktop(i)=1
511 kstart=kpbl(i)
512 if(kpbl(i).lt.5)kstart=kbcon(i)
513! if(k_inv_layers(i,1).gt.0)then
514!! ktop(i)=min(k_inv_layers(i,1),k_inv_layers(i,2))
515 if(k_inv_layers(i,1).gt.0 .and. &
516 (po_cup(i,kstart)-po_cup(i,k_inv_layers(i,1))).lt.200.)then
517 ktop(i)=k_inv_layers(i,1)
518 else
519 do k=kbcon(i)+1,ktf
520 if((po_cup(i,kstart)-po_cup(i,k)).gt.200.)then
521 ktop(i)=k
522 exit
523 endif
524 enddo
525 endif
526 endif
527 enddo
528!$acc end parallel
530 call rates_up_pdf(rand_vmas,ipr,'shallow',ktop,ierr,po_cup,entr_rate_2d,hkbo,heo,heso_cup,zo_cup, &
531 xland1,kstabi,k22,kbcon,its,ite,itf,kts,kte,ktf,zuo,kpbl,ktopx,kbcon,pmin_lev)
532!$acc kernels
533 do i=its,itf
534 if(ierr(i).eq.0)then
535! do k=maxloc(zuo(i,:),1),1,-1 ! ktop(i)-1,1,-1
536! if(zuo(i,k).lt.1.e-6)then
537! k22(i)=k+1
538! start_level(i)=k22(i)
539! exit
540! endif
541! enddo
542 if(k22(i).gt.1)then
543!$acc loop independent
544 do k=1,k22(i)-1
545 zuo(i,k)=0.
546 zu(i,k)=0.
547 xzu(i,k)=0.
548 enddo
549 endif
550!$acc loop seq
551 do k=maxloc(zuo(i,:),1),ktop(i)
552 if(zuo(i,k).lt.1.e-6)then
553 ktop(i)=k-1
554 exit
555 endif
556 enddo
557!$acc loop independent
558 do k=k22(i),ktop(i)
559 xzu(i,k)= zuo(i,k)
560 zu(i,k)= zuo(i,k)
561 enddo
562!$acc loop independent
563 do k=ktop(i)+1,ktf
564 zuo(i,k)=0.
565 zu(i,k)=0.
566 xzu(i,k)=0.
567 enddo
568 k22(i)=max(2,k22(i))
569 endif
570 enddo
571!$acc end kernels
572!
574!
575 call get_lateral_massflux(itf,ktf, its,ite, kts,kte &
576 ,ierr,ktop,zo_cup,zuo,cd,entr_rate_2d &
577 ,up_massentro, up_massdetro ,up_massentr, up_massdetr &
578 ,2,kbcon,k22,up_massentru,up_massdetru,lambau)
579!$acc kernels
580 do k=kts,ktf
581 do i=its,itf
582 hc(i,k)=0.
583 qco(i,k)=0.
584 qrco(i,k)=0.
585 dby(i,k)=0.
586 hco(i,k)=0.
587 dbyo(i,k)=0.
588 enddo
589 enddo
590 do i=its,itf
591 if(ierr(i) /= 0) cycle
592 do k=1,start_level(i)
593 uc(i,k)=u_cup(i,k)
594 vc(i,k)=v_cup(i,k)
595 enddo
596 do k=1,start_level(i)-1
597 hc(i,k)=he_cup(i,k)
598 hco(i,k)=heo_cup(i,k)
599 enddo
600 k=start_level(i)
601 hc(i,k)=hkb(i)
602 hco(i,k)=hkbo(i)
603 enddo
604!$acc end kernels
605!
606!
607
608!$acc parallel loop private(ki,qaver,k,trash,trash2,dz,dp)
609 do 42 i=its,itf
610 dbyt(i,:)=0.
611 if(ierr(i) /= 0) cycle
612!$acc loop seq
613 do k=start_level(i)+1,ktop(i)
614 hc(i,k)=(hc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*hc(i,k-1)+ &
615 up_massentr(i,k-1)*he(i,k-1)) / &
616 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
617 uc(i,k)=(uc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*uc(i,k-1)+ &
618 up_massentr(i,k-1)*us(i,k-1)) / &
619 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
620 vc(i,k)=(vc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*vc(i,k-1)+ &
621 up_massentr(i,k-1)*vs(i,k-1))/ &
622 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
623 dby(i,k)=max(0.,hc(i,k)-hes_cup(i,k))
624 hco(i,k)=(hco(i,k-1)*zuo(i,k-1)-.5*up_massdetro(i,k-1)*hco(i,k-1)+ &
625 up_massentro(i,k-1)*heo(i,k-1)) / &
626 (zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
627 dbyo(i,k)=hco(i,k)-heso_cup(i,k)
628 dz=zo_cup(i,k+1)-zo_cup(i,k)
629 if(k.ge.kbcon(i))dbyt(i,k)=dbyt(i,k-1)+dbyo(i,k)*dz
630 enddo
631 ki=maxloc(dbyt(i,:),1)
632 if(ktop(i).gt.ki+1)then
633 ktop(i)=ki+1
634 zuo(i,ktop(i)+1:ktf)=0.
635 zu(i,ktop(i)+1:ktf)=0.
636 cd(i,ktop(i)+1:ktf)=0.
637 up_massdetro(i,ktop(i))=zuo(i,ktop(i))
638! up_massentro(i,ktop(i))=0.
639 up_massentro(i,ktop(i):ktf)=0.
640 up_massdetro(i,ktop(i)+1:ktf)=0.
641 entr_rate_2d(i,ktop(i)+1:ktf)=0.
642
643! ierr(i)=423
644 endif
645
646 if(ktop(i).lt.kbcon(i)+1)then
647 ierr(i)=5
648#ifndef _OPENACC
649 ierrc(i)='ktop is less than kbcon+1'
650#endif
651 go to 42
652 endif
653 if(ktop(i).gt.ktf-2)then
654 ierr(i)=5
655#ifndef _OPENACC
656 ierrc(i)="ktop is larger than ktf-2"
657#endif
658 go to 42
659 endif
660!
661 call get_cloud_bc(kte,qo_cup(i,1:kte),qaver,k22(i),zero)
662 qaver = qaver + zqexec(i)
663 do k=1,start_level(i)-1
664 qco(i,k)= qo_cup(i,k)
665 enddo
666 k=start_level(i)
667 qco(i,k)= qaver
668!
669!$acc loop seq
670 do k=start_level(i)+1,ktop(i)
671 trash=qeso_cup(i,k)+(1./xlv)*(gammao_cup(i,k) &
672 /(1.+gammao_cup(i,k)))*dbyo(i,k)
673 !- total water liq+vapour
674 trash2 = qco(i,k-1) ! +qrco(i,k-1)
675 qco(i,k)= (trash2* ( zuo(i,k-1)-0.5*up_massdetr(i,k-1)) + &
676 up_massentr(i,k-1)*qo(i,k-1)) / &
677 (zuo(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
678
679 if(qco(i,k)>=trash ) then
680 dz=z_cup(i,k)-z_cup(i,k-1)
681 ! cloud liquid water
682 c1d(i,k)=c1_shal! 0. !.02*up_massdetr(i,k-1)
683 clw_all(i,k)=max(0._kind_phys,qco(i,k)-trash)
684 qrco(i,k)= (qco(i,k)-trash)/(1.+(c0_shal+c1d(i,k))*dz)
685 if(qrco(i,k).lt.0.)then ! hli new test 02/12/19
686 qrco(i,k)=0.
687 !c1d(i,k)=0.
688 endif
689 pwo(i,k)=c0_shal*dz*qrco(i,k)*zuo(i,k)
690 ! cloud water vapor
691 qco(i,k)= trash+qrco(i,k)
692
693 else
694 qrco(i,k)= 0.0
695 endif
696 cupclw(i,k)=qrco(i,k)
697 enddo
698 trash=0.
699 trash2=0.
700!$acc loop independent
701 do k=k22(i)+1,ktop(i)
702 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
703 cnvwt(i,k)=zuo(i,k)*cupclw(i,k)*g/dp
704!$acc atomic
705 trash2=trash2+entr_rate_2d(i,k)
706!$acc atomic
707 qco(i,k)=qco(i,k)-qrco(i,k)
708 enddo
709!$acc loop independent
710 do k=k22(i)+1,max(kbcon(i),k22(i)+1)
711!$acc atomic
712 trash=trash+entr_rate_2d(i,k)
713 enddo
714!$acc loop independent
715 do k=ktop(i)+1,ktf-1
716 hc(i,k)=hes_cup(i,k)
717 hco(i,k)=heso_cup(i,k)
718 qco(i,k)=qeso_cup(i,k)
719 uc(i,k)=u_cup(i,k)
720 vc(i,k)=v_cup(i,k)
721 qrco(i,k)=0.
722 dby(i,k)=0.
723 dbyo(i,k)=0.
724 zu(i,k)=0.
725 xzu(i,k)=0.
726 zuo(i,k)=0.
727 enddo
728 42 continue
729!$acc end parallel
730!
731!--- calculate workfunctions for updrafts
732!
733 if(make_calc_for_xk) then
734 call cup_up_aa0(aa0,z,zu,dby,gamma_cup,t_cup, &
735 kbcon,ktop,ierr, &
736 itf,ktf, its,ite, kts,kte)
737 call cup_up_aa0(aa1,zo,zuo,dbyo,gammao_cup,tn_cup, &
738 kbcon,ktop,ierr, &
739 itf,ktf, its,ite, kts,kte)
740!$acc kernels
741 do i=its,itf
742 if(ierr(i) == 0)then
743 if(aa1(i) <= 0.)then
744 ierr(i)=17
745#ifndef _OPENACC
746 ierrc(i)="cloud work function zero"
747#endif
748 endif
749 endif
750 enddo
751!$acc end kernels
752 endif
753
754!LB: insert calls to updraft vertical veloicity and prognostic area fraction here:
755 call calculate_updraft_velocity(its,itf,ktf,ite,kts,kte,ierr,progsigma, &
756 k22,kbcon,ktop,zo,entr_rate_2d,cd,fv,r_d,el2orc,qeso,tn,qo,po,dbyo, &
757 clw_all,qrco,delp,zu,wu2,omega_u,zeta,wc,omegac,zdqca)
758
759
760!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
761!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
762
763!
764!--- change per unit mass that a model cloud would modify the environment
765!
766!--- 1. in bottom layer
767!
768!$acc kernels
769 do k=kts,kte
770 do i=its,itf
771 dellah(i,k)=0.
772 dellaq(i,k)=0.
773 dellaqc(i,k)=0.
774 dellu(i,k)=0.
775 dellv(i,k)=0.
776 enddo
777 enddo
778!$acc end kernels
779!
780!---------------------------------------------- cloud level ktop
781!
782!- - - - - - - - - - - - - - - - - - - - - - - - model level ktop-1
783! . . .
784! . . .
785! . . .
786! . . .
787! . . .
788! . . .
789!
790!---------------------------------------------- cloud level k+2
791!
792!- - - - - - - - - - - - - - - - - - - - - - - - model level k+1
793!
794!---------------------------------------------- cloud level k+1
795!
796!- - - - - - - - - - - - - - - - - - - - - - - - model level k
797!
798!---------------------------------------------- cloud level k
799!
800! . . .
801! . . .
802! . . .
803! . . .
804! . . .
805! . . .
806! . . .
807! . . .
808! . . .
809! . . .
810!
811!---------------------------------------------- cloud level 3
812!
813!- - - - - - - - - - - - - - - - - - - - - - - - model level 2
814!
815!---------------------------------------------- cloud level 2
816!
817!- - - - - - - - - - - - - - - - - - - - - - - - model level 1
818 trash2=0.
819!$acc kernels
820!$acc loop independent
821 do i=its,itf
822 if(ierr(i).eq.0)then
823 dp=100.*(po_cup(i,1)-po_cup(i,2))
824 dellu(i,1)= -zuo(i,2)*(uc(i,2)-u_cup(i,2)) *g/dp
825 dellv(i,1)= -zuo(i,2)*(vc(i,2)-v_cup(i,2)) *g/dp
826 dellah(i,1)=-zuo(i,2)*(hco(i,2)-heo_cup(i,2))*g/dp
827
828 dellaq(i,1)=-zuo(i,2)*(qco(i,2)-qo_cup(i,2))*g/dp
829
830 do k=k22(i),ktop(i)
831 ! entrainment/detrainment for updraft
832 entup=up_massentro(i,k)
833 detup=up_massdetro(i,k)
834 totmas=detup-entup+zuo(i,k+1)-zuo(i,k)
835#ifndef _OPENACC
836 if(abs(totmas).gt.1.e-6)then
837 write(0,*)'*********************',i,k,totmas
838 write(0,*)k22(i),kbcon(i),ktop(i)
839 endif
840#endif
841 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
842 dellah(i,k) =-(zuo(i,k+1)*(hco(i,k+1)-heo_cup(i,k+1) )- &
843 zuo(i,k )*(hco(i,k )-heo_cup(i,k ) ))*g/dp
844
845 !-- take out cloud liquid water for detrainment
846 dz=zo_cup(i,k+1)-zo_cup(i,k)
847 if(k.lt.ktop(i) .and. c1d(i,k) > 0)then
848 dellaqc(i,k)= zuo(i,k)*c1d(i,k)*qrco(i,k)*dz/dp*g ! detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
849 else
850 dellaqc(i,k)=detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
851! dellaqc(i,k)= detup*qrco(i,k) *g/dp
852 endif
853
854 !-- condensation source term = detrained + flux divergence of
855 !-- cloud liquid water (qrco)
856 c_up = dellaqc(i,k)+(zuo(i,k+1)* qrco(i,k+1) - &
857 zuo(i,k )* qrco(i,k ) )*g/dp
858! c_up = dellaqc(i,k)
859 !-- water vapor budget (flux divergence of q_up-q_env - condensation
860 !term)
861 dellaq(i,k) =-(zuo(i,k+1)*(qco(i,k+1)-qo_cup(i,k+1) ) - &
862 zuo(i,k )*(qco(i,k )-qo_cup(i,k ) ) )*g/dp &
863 - c_up - 0.5*(pwo(i,k)+pwo(i,k+1))*g/dp
864 dellu(i,k) =-(zuo(i,k+1)*(uc(i,k+1)-u_cup(i,k+1) ) - &
865 zuo(i,k )*(uc(i,k )-u_cup(i,k ) ) )*g/dp
866 dellv(i,k) =-(zuo(i,k+1)*(vc(i,k+1)-v_cup(i,k+1) ) - &
867 zuo(i,k )*(vc(i,k )-v_cup(i,k ) ) )*g/dp
868
869 enddo
870 endif
871 enddo
872!$acc end kernels
873
874!
875!--- using dellas, calculate changed environmental profiles
876!
877 mbdt=.5 !3.e-4
878!$acc kernels
879 do k=kts,ktf
880 do i=its,itf
881 dellat(i,k)=0.
882 if(ierr(i)/=0)cycle
883 xhe(i,k)=dellah(i,k)*mbdt+heo(i,k)
884 xq(i,k)=max(1.e-16,(dellaq(i,k)+dellaqc(i,k))*mbdt+qo(i,k))
885 dellat(i,k)=(1./cp)*(dellah(i,k)-xlv*(dellaq(i,k)))
886 xt(i,k)= (-dellaqc(i,k)*xlv/cp+dellat(i,k))*mbdt+tn(i,k)
887 xt(i,k)= max(190.,xt(i,k))
888
889 enddo
890 enddo
891 do i=its,itf
892 if(ierr(i).eq.0)then
893! xhkb(i)=hkbo(i)+(dellah(i,k22(i)))*mbdt
894 xhe(i,ktf)=heo(i,ktf)
895 xq(i,ktf)=qo(i,ktf)
896 xt(i,ktf)=tn(i,ktf)
897 endif
898 enddo
899!$acc end kernels
900!
901!
902 if(make_calc_for_xk) then
903!
904!--- calculate moist static energy, heights, qes
905!
906 call cup_env(xz,xqes,xhe,xhes,xt,xq,po,z1, &
907 psur,ierr,tcrit,-1, &
908 itf,ktf, &
909 its,ite, kts,kte)
910!
911!--- environmental values on cloud levels
912!
913 call cup_env_clev(xt,xqes,xq,xhe,xhes,xz,po,xqes_cup,xq_cup, &
914 xhe_cup,xhes_cup,xz_cup,po_cup,gamma_cup,xt_cup,psur, &
915 ierr,z1, &
916 itf,ktf, &
917 its,ite, kts,kte)
918!
919!
920!**************************** static control
921!$acc kernels
922 do k=kts,ktf
923 do i=its,itf
924 xhc(i,k)=0.
925 xdby(i,k)=0.
926 enddo
927 enddo
928!$acc end kernels
929
930!$acc parallel loop private(x_add)
931 do i=its,itf
932 if(ierr(i).eq.0)then
933 x_add = xlv*zqexec(i)+cp*ztexec(i)
934 call get_cloud_bc(kte,xhe_cup(i,1:kte),xhkb(i),k22(i),x_add)
935 do k=1,start_level(i)-1
936 xhc(i,k)=xhe_cup(i,k)
937 enddo
938 k=start_level(i)
939 xhc(i,k)=xhkb(i)
940 endif !ierr
941 enddo
942!$acc end parallel
943!
944!
945!$acc kernels
946 do i=its,itf
947 if(ierr(i).eq.0)then
948 xzu(i,1:ktf)=zuo(i,1:ktf)
949!$acc loop seq
950 do k=start_level(i)+1,ktop(i)
951 xhc(i,k)=(xhc(i,k-1)*xzu(i,k-1)-.5*up_massdetro(i,k-1)*xhc(i,k-1)+ &
952 up_massentro(i,k-1)*xhe(i,k-1)) / &
953 (xzu(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
954 xdby(i,k)=xhc(i,k)-xhes_cup(i,k)
955 enddo
956!$acc loop independent
957 do k=ktop(i)+1,ktf
958 xhc(i,k)=xhes_cup(i,k)
959 xdby(i,k)=0.
960 xzu(i,k)=0.
961 enddo
962 endif
963 enddo
964!$acc end kernels
965
966
968! equation 8, call progsigma_calc() to compute updraft area fraction based on a moisture budget
969 if(progsigma)then
970 flag_shallow = .true.
971 flag_mid = .false.
972 do k=kts,ktf
973 do i=its,itf
974 del(i,k) = delp(i,k)*0.001
975 enddo
976 enddo
977 do i=its,itf
978 cnvflg(i)=.false.
979 enddo
980 do i=its,itf
981 if(ierr(i)==0)then
982 cnvflg(i)=.true.
983 endif
984 enddo
985 call progsigma_calc(itf,ktf,flag_init,flag_restart,flag_shallow, &
986 flag_mid,del,tmf,qmicro,dbyo,zdqca,omega_u,zeta,xlv,dtime, &
987 forceqv_spechum,kbcon,ktop,cnvflg,betascu,betamcu,betadcu, &
988 sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
989
990 endif
991
992!--- workfunctions for updraft
993!
994 call cup_up_aa0(xaa0,xz,xzu,xdby,gamma_cup,xt_cup, &
995 kbcon,ktop,ierr, &
996 itf,ktf, &
997 its,ite, kts,kte)
998!
999 endif
1000!
1001!
1002! now for shallow forcing
1003!
1004!$acc kernels
1005!$acc loop private(xff_shal)
1006 do i=its,itf
1007 xmb(i)=0.
1008
1009 if(progsigma)then
1010 gravinv = 1./g
1011 if(ierr(i)==0)then
1012 xmb(i) = sigmab(i)*((-1.0*omegac(i))*gravinv)
1013 if (dx(i) < 10.e3) then
1014 scaldfunc(i)=(1.-sigmab(i))*(1.-sigmab(i))
1015 scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.)
1016 else
1017 scaldfunc(i) = 1.0
1018 endif
1019 xmb(i)=scaldfunc(i)*xmb(i)
1020 endif
1021
1022 else
1023
1024 xff_shal(1:3)=0.
1025 if(ierr(i).eq.0)then
1026 xmbmax(i)=1.0
1027! xmbmax(i)=100.*(p(i,kbcon(i))-p(i,kbcon(i)+1))/(g*dtime)
1028!
1029!-stabilization closure
1030 xkshal=(xaa0(i)-aa1(i))/mbdt
1031 if(xkshal.le.0.and.xkshal.gt.-.01*mbdt) &
1032 xkshal=-.01*mbdt
1033 if(xkshal.gt.0.and.xkshal.lt.1.e-2) &
1034 xkshal=1.e-2
1035
1036 xff_shal(1)=max(0.,-(aa1(i)-aa0(i))/(xkshal*dtime))
1037!
1038!- closure from grant (2001)
1039 xff_shal(2)=.03*zws(i)
1040!- boundary layer qe closure
1041 blqe=0.
1042 trash=0.
1043 do k=1,kbcon(i)
1044 blqe=blqe+100.*dhdt(i,k)*(po_cup(i,k)-po_cup(i,k+1))/g
1045 enddo
1046 trash=max((hc(i,kbcon(i))-he_cup(i,kbcon(i))),1.e1)
1047 xff_shal(3)=max(0.,blqe/trash)
1048 xff_shal(3)=min(xmbmax(i),xff_shal(3))
1049!- average
1050 xmb(i)=(xff_shal(1)+xff_shal(2)+xff_shal(3))/3.
1051 xmb(i)=min(xmbmax(i),xmb(i))
1052 if(ichoice > 0)xmb(i)=min(xmbmax(i),xff_shal(ichoice))
1053 if(xmb(i) <= 0.)then
1054 ierr(i)=21
1055#ifndef _OPENACC
1056 ierrc(i)="21"
1057#endif
1058 endif
1059 endif
1060
1061 endif !progsigma
1062
1063 if(ierr(i).ne.0)then
1064 k22(i)=0
1065 kbcon(i)=0
1066 ktop(i)=0
1067 xmb(i)=0.
1068 outt(i,:)=0.
1069 outu(i,:)=0.
1070 outv(i,:)=0.
1071 outq(i,:)=0.
1072 outqc(i,:)=0.
1073 else if(ierr(i).eq.0)then
1074 xmb_out(i)=xmb(i)
1075!
1076! final tendencies
1077!
1078 pre(i)=0.
1079!$acc loop independent
1080 do k=2,ktop(i)
1081 outt(i,k)= dellat(i,k)*xmb(i)
1082 outq(i,k)= dellaq(i,k)*xmb(i)
1083 outqc(i,k)= dellaqc(i,k)*xmb(i)
1084!$acc atomic
1085 pre(i) = pre(i)+pwo(i,k)*xmb(i)
1086 enddo
1087 outt(i,1)= dellat(i,1)*xmb(i)
1088 outq(i,1)= dellaq(i,1)*xmb(i)
1089 outu(i,1)=dellu(i,1)*xmb(i)
1090 outv(i,1)=dellv(i,1)*xmb(i)
1091 do k=kts+1,ktop(i)
1092 outu(i,k)=.25*(dellu(i,k-1)+2.*dellu(i,k)+dellu(i,k+1))*xmb(i)
1093 outv(i,k)=.25*(dellv(i,k-1)+2.*dellv(i,k)+dellv(i,k+1))*xmb(i)
1094 enddo
1095
1096 endif
1097
1098 enddo
1099!
1100! since kinetic energy is being dissipated, add heating accordingly (from ecmwf)
1101!
1102 do i=its,itf
1103 if(ierr(i).eq.0) then
1104 dts=0.
1105 fpi=0.
1106 do k=kts,ktop(i)
1107 dp=(po_cup(i,k)-po_cup(i,k+1))*100.
1108!total ke dissiptaion estimate
1109 dts= dts -(outu(i,k)*us(i,k)+outv(i,k)*vs(i,k))*dp/g
1110! fpi needed for calcualtion of conversion to pot. energyintegrated
1111 fpi = fpi +sqrt(outu(i,k)*outu(i,k) + outv(i,k)*outv(i,k))*dp
1112 enddo
1113 if(fpi.gt.0.)then
1114 do k=kts,ktop(i)
1115 fp= sqrt((outu(i,k)*outu(i,k)+outv(i,k)*outv(i,k)))/fpi
1116 outt(i,k)=outt(i,k)+fp*dts*g/cp
1117 enddo
1118 endif
1119 endif
1120 enddo
1121!$acc end kernels
1122!
1123! done shallow
1124!--------------------------done------------------------------
1125!
1126! do k=1,30
1127! print*,'hlisq',qco(1,k),qrco(1,k),pwo(1,k)
1128! enddo
1129
1130 end subroutine cu_c3_sh_run
1132end module cu_c3_sh
subroutine rates_up_pdf(rand_vmas, ipr, name, ktop, ierr, p_cup, entr_rate_2d, hkbo, heo, heso_cup, z_cup, xland, kstabi, k22, kbcon, its, ite, itf, kts, kte, ktf, zuo, kpbl, ktopdby, csum, pmin_lev)
Driver for the normalized mass-flux routine.
subroutine calculate_updraft_velocity(its, itf, ktf, ite, kts, kte, ierr, progsigma, k22, kbcon, ktcon, zo, entr_rate_2d, cd, fv, rd, el2orc, qeso, to, qo, po, dbyo, clw_all, qlk, delp, zu, wu2, omega_u, zeta, wc, omegac, zdqca)
subroutine get_lateral_massflux(itf, ktf, its, ite, kts, kte, ierr, ktop, zo_cup, zuo, cd, entr_rate_2d, up_massentro, up_massdetro, up_massentr, up_massdetr, draft, kbcon, k22, up_massentru, up_massdetru, lambau)
Calculates mass entranment and detrainment rates.
subroutine get_cloud_bc(mzp, array, x_aver, k22, add)
Calculates the average value of a variable at the updraft originating level.
subroutine cup_minimi(array, ks, kend, kt, ierr, itf, ktf, its, ite, kts, kte)
Calculates the level at which the minimum value in an array occurs.
subroutine cup_env(z, qes, he, hes, t, q, p, z1, psur, ierr, tcrit, itest, itf, ktf, its, ite, kts, kte)
Calculates environmental moist static energy, saturation moist static energy, heights,...
subroutine get_inversion_layers(ierr, p_cup, t_cup, z_cup, qo_cup, qeso_cup, k_inv_layers, kstart, kend, dtempdz, itf, ktf, its, ite, kts, kte)
Finds temperature inversions using the first and second derivative of temperature.
subroutine cup_env_clev(t, qes, q, he, hes, z, p, qes_cup, q_cup, he_cup, hes_cup, z_cup, p_cup, gamma_cup, t_cup, psur, ierr, z1, itf, ktf, its, ite, kts, kte)
Calculates environmental values on cloud levels.
subroutine cup_up_aa0(aa0, z, zu, dby, gamma_cup, t_cup, kbcon, ktop, ierr, itf, ktf, its, ite, kts, kte)
Calculates the cloud work functions for updrafts.
subroutine cup_kbcon(ierrc, cap_inc, iloop_in, k22, kbcon, he_cup, hes_cup, hkb, ierr, kbmax, p_cup, cap_max, ztexec, zqexec, jprnt, itf, ktf, its, ite, kts, kte, z_cup, entr_rate, heo, imid)
Calculates the level of convective cloud base.
subroutine cu_c3_sh_run(us, vs, zo, t, q, z1, tn, qo, po, psur, dhdt, kpbl, rho, hfx, qfx, xland, ichoice, tcrit, dtime, zuo, xmb_out, kbcon, ktop, k22, ierr, ierrc, flag_init, flag_restart, fv, r_d, delp, tmf, qmicro, forceqv_spechum, betascu, betamcu, betadcu, sigmain, sigmaout, progsigma, dx, outt, outq, outqc, outu, outv, cnvwt, pre, cupclw, itf, ktf, its, ite, kts, kte, ipr, tropics)
Definition cu_c3_sh.F90:75