CCPP SciDoc v7.0.0  v7.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 go to 25
386 endif
387 enddo
388 25 continue
389!
390 kbmax(i)=min(kbmax(i),ktf/2)
391 endif
392 enddo
393!$acc end kernels
394
395!
396!
397!
399!
400!$acc parallel loop
401 do 36 i=its,itf
402 if(kpbl(i).gt.3)cap_max(i)=po_cup(i,kpbl(i))
403 if(ierr(i) == 0)then
404 k22(i)=maxloc(heo_cup(i,2:kbmax(i)),1)
405 k22(i)=max(2,k22(i))
406 if(k22(i).gt.kbmax(i))then
407 ierr(i)=2
408#ifndef _OPENACC
409 ierrc(i)="could not find k22"
410#endif
411 ktop(i)=0
412 k22(i)=0
413 kbcon(i)=0
414 endif
415 endif
416 36 continue
417!$acc end parallel
418!
421!
422!$acc parallel loop private(x_add)
423 do i=its,itf
424 if(ierr(i).eq.0)then
425 x_add = xlv*zqexec(i)+cp*ztexec(i)
426 call get_cloud_bc(kte,he_cup(i,1:kte),hkb(i),k22(i),x_add)
427 call get_cloud_bc(kte,heo_cup(i,1:kte),hkbo(i),k22(i),x_add)
428 endif ! ierr
429 enddo
430!$acc end parallel
431
432!joe-georg and saulo's new idea:
433
434!$acc kernels
435 do i=its,itf
436 do k=kts,ktf
437 dbyo(i,k)= 0. !hkbo(i)-heso_cup(i,k)
438 enddo
439 enddo
440!$acc end kernels
441
442
443 call cup_kbcon(ierrc,cap_max_increment,5,k22,kbcon,heo_cup,heso_cup, &
444 hkbo,ierr,kbmax,po_cup,cap_max, &
445 ztexec,zqexec, &
446 0,itf,ktf, &
447 its,ite, kts,kte, &
448 z_cup,entr_rate,heo,0)
449
451 call cup_minimi(heso_cup,kbcon,kbmax,kstabi,ierr, &
452 itf,ktf, &
453 its,ite, kts,kte)
454!
455 call get_inversion_layers(ierr,p_cup,t_cup,z_cup,q_cup,qes_cup,k_inv_layers,&
456 kbcon,kstabi,dtempdz,itf,ktf,its,ite, kts,kte)
457!
458!
459!$acc parallel loop private(frh,kstart,x_add)
460 do i=its,itf
461 entr_rate_2d(i,:)=entr_rate(i)
462 if(ierr(i) == 0)then
463 start_level(i)=k22(i)
464 x_add = xlv*zqexec(i)+cp*ztexec(i)
465 call get_cloud_bc(kte,he_cup(i,1:kte),hkb(i),k22(i),x_add)
466 if(kbcon(i).gt.ktf-4)then
467 ierr(i)=231
468 endif
469 do k=kts,ktf
470 frh = 2.*min(qo_cup(i,k)/qeso_cup(i,k),1.)
471 entr_rate_2d(i,k)=entr_rate(i) !*(2.3-frh)
472 cd(i,k)=.75*entr_rate_2d(i,k)
473 enddo
474!
475! first estimate for shallow convection
476!
477 ktop(i)=1
478 kstart=kpbl(i)
479 if(kpbl(i).lt.5)kstart=kbcon(i)
480! if(k_inv_layers(i,1).gt.0)then
481!! ktop(i)=min(k_inv_layers(i,1),k_inv_layers(i,2))
482 if(k_inv_layers(i,1).gt.0 .and. &
483 (po_cup(i,kstart)-po_cup(i,k_inv_layers(i,1))).lt.200.)then
484 ktop(i)=k_inv_layers(i,1)
485 else
486 do k=kbcon(i)+1,ktf
487 if((po_cup(i,kstart)-po_cup(i,k)).gt.200.)then
488 ktop(i)=k
489 exit
490 endif
491 enddo
492 endif
493 endif
494 enddo
495!$acc end parallel
497 call rates_up_pdf(rand_vmas,ipr,'shallow',ktop,ierr,po_cup,entr_rate_2d,hkbo,heo,heso_cup,zo_cup, &
498 xland1,kstabi,k22,kbcon,its,ite,itf,kts,kte,ktf,zuo,kpbl,ktopx,kbcon,pmin_lev)
499!$acc kernels
500 do i=its,itf
501 if(ierr(i).eq.0)then
502! do k=maxloc(zuo(i,:),1),1,-1 ! ktop(i)-1,1,-1
503! if(zuo(i,k).lt.1.e-6)then
504! k22(i)=k+1
505! start_level(i)=k22(i)
506! exit
507! endif
508! enddo
509 if(k22(i).gt.1)then
510!$acc loop independent
511 do k=1,k22(i)-1
512 zuo(i,k)=0.
513 zu(i,k)=0.
514 xzu(i,k)=0.
515 enddo
516 endif
517!$acc loop seq
518 do k=maxloc(zuo(i,:),1),ktop(i)
519 if(zuo(i,k).lt.1.e-6)then
520 ktop(i)=k-1
521 exit
522 endif
523 enddo
524!$acc loop independent
525 do k=k22(i),ktop(i)
526 xzu(i,k)= zuo(i,k)
527 zu(i,k)= zuo(i,k)
528 enddo
529!$acc loop independent
530 do k=ktop(i)+1,ktf
531 zuo(i,k)=0.
532 zu(i,k)=0.
533 xzu(i,k)=0.
534 enddo
535 k22(i)=max(2,k22(i))
536 endif
537 enddo
538!$acc end kernels
539!
541!
542 call get_lateral_massflux(itf,ktf, its,ite, kts,kte &
543 ,ierr,ktop,zo_cup,zuo,cd,entr_rate_2d &
544 ,up_massentro, up_massdetro ,up_massentr, up_massdetr &
545 ,2,kbcon,k22,up_massentru,up_massdetru,lambau)
546!$acc kernels
547 do k=kts,ktf
548 do i=its,itf
549 hc(i,k)=0.
550 qco(i,k)=0.
551 qrco(i,k)=0.
552 dby(i,k)=0.
553 hco(i,k)=0.
554 dbyo(i,k)=0.
555 enddo
556 enddo
557 do i=its,itf
558 if(ierr(i) /= 0) cycle
559 do k=1,start_level(i)
560 uc(i,k)=u_cup(i,k)
561 vc(i,k)=v_cup(i,k)
562 enddo
563 do k=1,start_level(i)-1
564 hc(i,k)=he_cup(i,k)
565 hco(i,k)=heo_cup(i,k)
566 enddo
567 k=start_level(i)
568 hc(i,k)=hkb(i)
569 hco(i,k)=hkbo(i)
570 enddo
571!$acc end kernels
572!
573!
574
575!$acc parallel loop private(ki,qaver,k,trash,trash2,dz,dp)
576 do 42 i=its,itf
577 dbyt(i,:)=0.
578 if(ierr(i) /= 0) cycle
579!$acc loop seq
580 do k=start_level(i)+1,ktop(i)
581 hc(i,k)=(hc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*hc(i,k-1)+ &
582 up_massentr(i,k-1)*he(i,k-1)) / &
583 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
584 uc(i,k)=(uc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*uc(i,k-1)+ &
585 up_massentr(i,k-1)*us(i,k-1)) / &
586 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
587 vc(i,k)=(vc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*vc(i,k-1)+ &
588 up_massentr(i,k-1)*vs(i,k-1))/ &
589 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
590 dby(i,k)=max(0.,hc(i,k)-hes_cup(i,k))
591 hco(i,k)=(hco(i,k-1)*zuo(i,k-1)-.5*up_massdetro(i,k-1)*hco(i,k-1)+ &
592 up_massentro(i,k-1)*heo(i,k-1)) / &
593 (zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
594 dbyo(i,k)=hco(i,k)-heso_cup(i,k)
595 dz=zo_cup(i,k+1)-zo_cup(i,k)
596 if(k.ge.kbcon(i))dbyt(i,k)=dbyt(i,k-1)+dbyo(i,k)*dz
597 enddo
598 ki=maxloc(dbyt(i,:),1)
599 if(ktop(i).gt.ki+1)then
600 ktop(i)=ki+1
601 zuo(i,ktop(i)+1:ktf)=0.
602 zu(i,ktop(i)+1:ktf)=0.
603 cd(i,ktop(i)+1:ktf)=0.
604 up_massdetro(i,ktop(i))=zuo(i,ktop(i))
605! up_massentro(i,ktop(i))=0.
606 up_massentro(i,ktop(i):ktf)=0.
607 up_massdetro(i,ktop(i)+1:ktf)=0.
608 entr_rate_2d(i,ktop(i)+1:ktf)=0.
609
610! ierr(i)=423
611 endif
612
613 if(ktop(i).lt.kbcon(i)+1)then
614 ierr(i)=5
615#ifndef _OPENACC
616 ierrc(i)='ktop is less than kbcon+1'
617#endif
618 go to 42
619 endif
620 if(ktop(i).gt.ktf-2)then
621 ierr(i)=5
622#ifndef _OPENACC
623 ierrc(i)="ktop is larger than ktf-2"
624#endif
625 go to 42
626 endif
627!
628 call get_cloud_bc(kte,qo_cup(i,1:kte),qaver,k22(i),zero)
629 qaver = qaver + zqexec(i)
630 do k=1,start_level(i)-1
631 qco(i,k)= qo_cup(i,k)
632 enddo
633 k=start_level(i)
634 qco(i,k)= qaver
635!
636!$acc loop seq
637 do k=start_level(i)+1,ktop(i)
638 trash=qeso_cup(i,k)+(1./xlv)*(gammao_cup(i,k) &
639 /(1.+gammao_cup(i,k)))*dbyo(i,k)
640 !- total water liq+vapour
641 trash2 = qco(i,k-1) ! +qrco(i,k-1)
642 qco(i,k)= (trash2* ( zuo(i,k-1)-0.5*up_massdetr(i,k-1)) + &
643 up_massentr(i,k-1)*qo(i,k-1)) / &
644 (zuo(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
645
646 if(qco(i,k)>=trash ) then
647 dz=z_cup(i,k)-z_cup(i,k-1)
648 ! cloud liquid water
649 c1d(i,k)=.02*up_massdetr(i,k-1)
650 qrco(i,k)= (qco(i,k)-trash)/(1.+(c0_shal+c1d(i,k))*dz)
651 if(qrco(i,k).lt.0.)then ! hli new test 02/12/19
652 qrco(i,k)=0.
653 c1d(i,k)=0.
654 endif
655 pwo(i,k)=c0_shal*dz*qrco(i,k)*zuo(i,k)
656 ! cloud water vapor
657 qco(i,k)= trash+qrco(i,k)
658
659 else
660 qrco(i,k)= 0.0
661 endif
662 cupclw(i,k)=qrco(i,k)
663 enddo
664 trash=0.
665 trash2=0.
666!$acc loop independent
667 do k=k22(i)+1,ktop(i)
668 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
669 cnvwt(i,k)=zuo(i,k)*cupclw(i,k)*g/dp
670!$acc atomic
671 trash2=trash2+entr_rate_2d(i,k)
672!$acc atomic
673 qco(i,k)=qco(i,k)-qrco(i,k)
674 enddo
675!$acc loop independent
676 do k=k22(i)+1,max(kbcon(i),k22(i)+1)
677!$acc atomic
678 trash=trash+entr_rate_2d(i,k)
679 enddo
680!$acc loop independent
681 do k=ktop(i)+1,ktf-1
682 hc(i,k)=hes_cup(i,k)
683 hco(i,k)=heso_cup(i,k)
684 qco(i,k)=qeso_cup(i,k)
685 uc(i,k)=u_cup(i,k)
686 vc(i,k)=v_cup(i,k)
687 qrco(i,k)=0.
688 dby(i,k)=0.
689 dbyo(i,k)=0.
690 zu(i,k)=0.
691 xzu(i,k)=0.
692 zuo(i,k)=0.
693 enddo
694 42 continue
695!$acc end parallel
696!
697!--- calculate workfunctions for updrafts
698!
699 if(make_calc_for_xk) then
700 call cup_up_aa0(aa0,z,zu,dby,gamma_cup,t_cup, &
701 kbcon,ktop,ierr, &
702 itf,ktf, its,ite, kts,kte)
703 call cup_up_aa0(aa1,zo,zuo,dbyo,gammao_cup,tn_cup, &
704 kbcon,ktop,ierr, &
705 itf,ktf, its,ite, kts,kte)
706!$acc kernels
707 do i=its,itf
708 if(ierr(i) == 0)then
709 if(aa1(i) <= 0.)then
710 ierr(i)=17
711#ifndef _OPENACC
712 ierrc(i)="cloud work function zero"
713#endif
714 endif
715 endif
716 enddo
717!$acc end kernels
718 endif
719!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
720!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
721
722!
723!--- change per unit mass that a model cloud would modify the environment
724!
725!--- 1. in bottom layer
726!
727!$acc kernels
728 do k=kts,kte
729 do i=its,itf
730 dellah(i,k)=0.
731 dellaq(i,k)=0.
732 dellaqc(i,k)=0.
733 dellu(i,k)=0.
734 dellv(i,k)=0.
735 enddo
736 enddo
737!$acc end kernels
738!
739!---------------------------------------------- cloud level ktop
740!
741!- - - - - - - - - - - - - - - - - - - - - - - - model level ktop-1
742! . . .
743! . . .
744! . . .
745! . . .
746! . . .
747! . . .
748!
749!---------------------------------------------- cloud level k+2
750!
751!- - - - - - - - - - - - - - - - - - - - - - - - model level k+1
752!
753!---------------------------------------------- cloud level k+1
754!
755!- - - - - - - - - - - - - - - - - - - - - - - - model level k
756!
757!---------------------------------------------- cloud level k
758!
759! . . .
760! . . .
761! . . .
762! . . .
763! . . .
764! . . .
765! . . .
766! . . .
767! . . .
768! . . .
769!
770!---------------------------------------------- cloud level 3
771!
772!- - - - - - - - - - - - - - - - - - - - - - - - model level 2
773!
774!---------------------------------------------- cloud level 2
775!
776!- - - - - - - - - - - - - - - - - - - - - - - - model level 1
777 trash2=0.
778!$acc kernels
779!$acc loop independent
780 do i=its,itf
781 if(ierr(i).eq.0)then
782 dp=100.*(po_cup(i,1)-po_cup(i,2))
783 dellu(i,1)= -zuo(i,2)*(uc(i,2)-u_cup(i,2)) *g/dp
784 dellv(i,1)= -zuo(i,2)*(vc(i,2)-v_cup(i,2)) *g/dp
785 dellah(i,1)=-zuo(i,2)*(hco(i,2)-heo_cup(i,2))*g/dp
786
787 dellaq(i,1)=-zuo(i,2)*(qco(i,2)-qo_cup(i,2))*g/dp
788
789 do k=k22(i),ktop(i)
790 ! entrainment/detrainment for updraft
791 entup=up_massentro(i,k)
792 detup=up_massdetro(i,k)
793 totmas=detup-entup+zuo(i,k+1)-zuo(i,k)
794#ifndef _OPENACC
795 if(abs(totmas).gt.1.e-6)then
796 write(0,*)'*********************',i,k,totmas
797 write(0,*)k22(i),kbcon(i),ktop(i)
798 endif
799#endif
800 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
801 dellah(i,k) =-(zuo(i,k+1)*(hco(i,k+1)-heo_cup(i,k+1) )- &
802 zuo(i,k )*(hco(i,k )-heo_cup(i,k ) ))*g/dp
803
804 !-- take out cloud liquid water for detrainment
805 dz=zo_cup(i,k+1)-zo_cup(i,k)
806 if(k.lt.ktop(i) .and. c1d(i,k) > 0)then
807 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
808 else
809 dellaqc(i,k)=detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
810! dellaqc(i,k)= detup*qrco(i,k) *g/dp
811 endif
812
813 !-- condensation source term = detrained + flux divergence of
814 !-- cloud liquid water (qrco)
815 c_up = dellaqc(i,k)+(zuo(i,k+1)* qrco(i,k+1) - &
816 zuo(i,k )* qrco(i,k ) )*g/dp
817! c_up = dellaqc(i,k)
818 !-- water vapor budget (flux divergence of q_up-q_env - condensation
819 !term)
820 dellaq(i,k) =-(zuo(i,k+1)*(qco(i,k+1)-qo_cup(i,k+1) ) - &
821 zuo(i,k )*(qco(i,k )-qo_cup(i,k ) ) )*g/dp &
822 - c_up - 0.5*(pwo(i,k)+pwo(i,k+1))*g/dp
823 dellu(i,k) =-(zuo(i,k+1)*(uc(i,k+1)-u_cup(i,k+1) ) - &
824 zuo(i,k )*(uc(i,k )-u_cup(i,k ) ) )*g/dp
825 dellv(i,k) =-(zuo(i,k+1)*(vc(i,k+1)-v_cup(i,k+1) ) - &
826 zuo(i,k )*(vc(i,k )-v_cup(i,k ) ) )*g/dp
827
828 enddo
829 endif
830 enddo
831!$acc end kernels
832
833!
834!--- using dellas, calculate changed environmental profiles
835!
836 mbdt=.5 !3.e-4
837!$acc kernels
838 do k=kts,ktf
839 do i=its,itf
840 dellat(i,k)=0.
841 if(ierr(i)/=0)cycle
842 xhe(i,k)=dellah(i,k)*mbdt+heo(i,k)
843 xq(i,k)=max(1.e-16,(dellaq(i,k)+dellaqc(i,k))*mbdt+qo(i,k))
844 dellat(i,k)=(1./cp)*(dellah(i,k)-xlv*(dellaq(i,k)))
845 xt(i,k)= (-dellaqc(i,k)*xlv/cp+dellat(i,k))*mbdt+tn(i,k)
846 xt(i,k)= max(190.,xt(i,k))
847
848 enddo
849 enddo
850 do i=its,itf
851 if(ierr(i).eq.0)then
852! xhkb(i)=hkbo(i)+(dellah(i,k22(i)))*mbdt
853 xhe(i,ktf)=heo(i,ktf)
854 xq(i,ktf)=qo(i,ktf)
855 xt(i,ktf)=tn(i,ktf)
856 endif
857 enddo
858!$acc end kernels
859!
860!
861 if(make_calc_for_xk) then
862!
863!--- calculate moist static energy, heights, qes
864!
865 call cup_env(xz,xqes,xhe,xhes,xt,xq,po,z1, &
866 psur,ierr,tcrit,-1, &
867 itf,ktf, &
868 its,ite, kts,kte)
869!
870!--- environmental values on cloud levels
871!
872 call cup_env_clev(xt,xqes,xq,xhe,xhes,xz,po,xqes_cup,xq_cup, &
873 xhe_cup,xhes_cup,xz_cup,po_cup,gamma_cup,xt_cup,psur, &
874 ierr,z1, &
875 itf,ktf, &
876 its,ite, kts,kte)
877!
878!
879!**************************** static control
880!$acc kernels
881 do k=kts,ktf
882 do i=its,itf
883 xhc(i,k)=0.
884 xdby(i,k)=0.
885 enddo
886 enddo
887!$acc end kernels
888
889!$acc parallel loop private(x_add)
890 do i=its,itf
891 if(ierr(i).eq.0)then
892 x_add = xlv*zqexec(i)+cp*ztexec(i)
893 call get_cloud_bc(kte,xhe_cup(i,1:kte),xhkb(i),k22(i),x_add)
894 do k=1,start_level(i)-1
895 xhc(i,k)=xhe_cup(i,k)
896 enddo
897 k=start_level(i)
898 xhc(i,k)=xhkb(i)
899 endif !ierr
900 enddo
901!$acc end parallel
902!
903!
904!$acc kernels
905 do i=its,itf
906 if(ierr(i).eq.0)then
907 xzu(i,1:ktf)=zuo(i,1:ktf)
908!$acc loop seq
909 do k=start_level(i)+1,ktop(i)
910 xhc(i,k)=(xhc(i,k-1)*xzu(i,k-1)-.5*up_massdetro(i,k-1)*xhc(i,k-1)+ &
911 up_massentro(i,k-1)*xhe(i,k-1)) / &
912 (xzu(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
913 xdby(i,k)=xhc(i,k)-xhes_cup(i,k)
914 enddo
915!$acc loop independent
916 do k=ktop(i)+1,ktf
917 xhc(i,k)=xhes_cup(i,k)
918 xdby(i,k)=0.
919 xzu(i,k)=0.
920 enddo
921 endif
922 enddo
923!$acc end kernels
924
925!
926!--- workfunctions for updraft
927!
928 call cup_up_aa0(xaa0,xz,xzu,xdby,gamma_cup,xt_cup, &
929 kbcon,ktop,ierr, &
930 itf,ktf, &
931 its,ite, kts,kte)
932!
933 endif
934!
935!
936! now for shallow forcing
937!
938!$acc kernels
939!$acc loop private(xff_shal)
940 do i=its,itf
941 xmb(i)=0.
942 xff_shal(1:3)=0.
943 if(ierr(i).eq.0)then
944 xmbmax(i)=1.0
945! xmbmax(i)=100.*(p(i,kbcon(i))-p(i,kbcon(i)+1))/(g*dtime)
946!
947!-stabilization closure
948 xkshal=(xaa0(i)-aa1(i))/mbdt
949 if(xkshal.le.0.and.xkshal.gt.-.01*mbdt) &
950 xkshal=-.01*mbdt
951 if(xkshal.gt.0.and.xkshal.lt.1.e-2) &
952 xkshal=1.e-2
953
954 xff_shal(1)=max(0.,-(aa1(i)-aa0(i))/(xkshal*dtime))
955!
956!- closure from grant (2001)
957 xff_shal(2)=.03*zws(i)
958!- boundary layer qe closure
959 blqe=0.
960 trash=0.
961 do k=1,kbcon(i)
962 blqe=blqe+100.*dhdt(i,k)*(po_cup(i,k)-po_cup(i,k+1))/g
963 enddo
964 trash=max((hc(i,kbcon(i))-he_cup(i,kbcon(i))),1.e1)
965 xff_shal(3)=max(0.,blqe/trash)
966 xff_shal(3)=min(xmbmax(i),xff_shal(3))
967!- average
968 xmb(i)=(xff_shal(1)+xff_shal(2)+xff_shal(3))/3.
969 xmb(i)=min(xmbmax(i),xmb(i))
970 if(ichoice > 0)xmb(i)=min(xmbmax(i),xff_shal(ichoice))
971 if(xmb(i) <= 0.)then
972 ierr(i)=21
973#ifndef _OPENACC
974 ierrc(i)="21"
975#endif
976 endif
977 endif
978 if(ierr(i).ne.0)then
979 k22(i)=0
980 kbcon(i)=0
981 ktop(i)=0
982 xmb(i)=0.
983 outt(i,:)=0.
984 outu(i,:)=0.
985 outv(i,:)=0.
986 outq(i,:)=0.
987 outqc(i,:)=0.
988 else if(ierr(i).eq.0)then
989 xmb_out(i)=xmb(i)
990!
991! final tendencies
992!
993 pre(i)=0.
994!$acc loop independent
995 do k=2,ktop(i)
996 outt(i,k)= dellat(i,k)*xmb(i)
997 outq(i,k)= dellaq(i,k)*xmb(i)
998 outqc(i,k)= dellaqc(i,k)*xmb(i)
999!$acc atomic
1000 pre(i) = pre(i)+pwo(i,k)*xmb(i)
1001 enddo
1002 outt(i,1)= dellat(i,1)*xmb(i)
1003 outq(i,1)= dellaq(i,1)*xmb(i)
1004 outu(i,1)=dellu(i,1)*xmb(i)
1005 outv(i,1)=dellv(i,1)*xmb(i)
1006 do k=kts+1,ktop(i)
1007 outu(i,k)=.25*(dellu(i,k-1)+2.*dellu(i,k)+dellu(i,k+1))*xmb(i)
1008 outv(i,k)=.25*(dellv(i,k-1)+2.*dellv(i,k)+dellv(i,k+1))*xmb(i)
1009 enddo
1010
1011 endif
1012 enddo
1013!
1014! since kinetic energy is being dissipated, add heating accordingly (from ecmwf)
1015!
1016 do i=its,itf
1017 if(ierr(i).eq.0) then
1018 dts=0.
1019 fpi=0.
1020 do k=kts,ktop(i)
1021 dp=(po_cup(i,k)-po_cup(i,k+1))*100.
1022!total ke dissiptaion estimate
1023 dts= dts -(outu(i,k)*us(i,k)+outv(i,k)*vs(i,k))*dp/g
1024! fpi needed for calcualtion of conversion to pot. energyintegrated
1025 fpi = fpi +sqrt(outu(i,k)*outu(i,k) + outv(i,k)*outv(i,k))*dp
1026 enddo
1027 if(fpi.gt.0.)then
1028 do k=kts,ktop(i)
1029 fp= sqrt((outu(i,k)*outu(i,k)+outv(i,k)*outv(i,k)))/fpi
1030 outt(i,k)=outt(i,k)+fp*dts*g/cp
1031 enddo
1032 endif
1033 endif
1034 enddo
1035!$acc end kernels
1036!
1037! done shallow
1038!--------------------------done------------------------------
1039!
1040! do k=1,30
1041! print*,'hlisq',qco(1,k),qrco(1,k),pwo(1,k)
1042! enddo
1043
1044 end subroutine cu_gf_sh_run
1046end 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:5203
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:4715
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:5403
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:3864
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:4743
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:3805
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