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