CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
cu_gf_deep.F90
1
3
5module cu_gf_deep
6 use machine , only : kind_phys
7 use physcons, only : qamin
8 real(kind=kind_phys), parameter::g=9.81
9 real(kind=kind_phys), parameter:: cp=1004.
10 real(kind=kind_phys), parameter:: xlv=2.5e6
11 real(kind=kind_phys), parameter::r_v=461.
12 real(kind=kind_phys), parameter :: tcrit=258.
14 real(kind=kind_phys), parameter:: c1= 0.003 !.002 ! .0005
16 integer, parameter :: irainevap=1
18 real(kind=kind_phys), parameter::frh_thresh = .9
20 real(kind=kind_phys), parameter::rh_thresh = .97
22 real(kind=kind_phys), parameter::betajb=1.2
24 integer, parameter:: use_excess=0
25 real(kind=kind_phys), parameter :: fluxtune=1.5
27 real(kind=kind_phys), parameter :: pgcd = 0.1
28!
30 integer, parameter :: autoconv=1 !2
31 integer, parameter :: aeroevap=1 !3
32 real(kind=kind_phys), parameter :: scav_factor = 0.5
33
34 real(kind=kind_phys), parameter :: dx_thresh = 6500.
36 integer, parameter:: maxens3=16
37
38!---meltglac-------------------------------------------------
39 logical, parameter :: melt_glac = .true.
40 real(kind=kind_phys), parameter :: &
41 t_0 = 273.16, &
42 t_ice = 250.16, &
43 xlf = 0.333e6
44!---meltglac-------------------------------------------------
45!-----srf-08aug2017-----begin
46 real(kind=kind_phys), parameter :: qrc_crit= 2.e-4
47!-----srf-08aug2017-----end
48
49contains
50
55 integer function my_maxloc1d(A,N)
56!$acc routine vector
57 implicit none
58 real(kind_phys), intent(in) :: a(:)
59 integer, intent(in) :: n
60
61 real(kind_phys) :: imaxval
62 integer :: i
63
64 imaxval = maxval(a)
65 my_maxloc1d = 1
66!$acc loop
67 do i = 1, n
68 if ( a(i) == imaxval ) then
69 my_maxloc1d = i
70 return
71 endif
72 end do
73 return
74 end function my_maxloc1d
75
78 subroutine cu_gf_deep_run( &
79 itf,ktf,its,ite, kts,kte &
80 ,dicycle & ! diurnal cycle flag
81 ,ichoice & ! choice of closure, use "0" for ensemble average
82 ,ipr & ! this flag can be used for debugging prints
83 ,ccn & ! not well tested yet
84 ,ccnclean &
85 ,dtime & ! dt over which forcing is applied
86 ,imid & ! flag to turn on mid level convection
87 ,kpbl & ! level of boundary layer height
88 ,dhdt & ! boundary layer forcing (one closure for shallow)
89 ,xland & ! land mask
90 ,zo & ! heights above surface
91 ,forcing & ! only diagnostic
92 ,t & ! t before forcing
93 ,q & ! q before forcing
94 ,z1 & ! terrain
95 ,tn & ! t including forcing
96 ,qo & ! q including forcing
97 ,po & ! pressure (mb)
98 ,psur & ! surface pressure (mb)
99 ,us & ! u on mass points
100 ,vs & ! v on mass points
101 ,rho & ! density
102 ,hfx & ! w/m2, positive upward
103 ,qfx & ! w/m2, positive upward
104 ,dx & ! dx is grid point dependent here
105 ,mconv & ! integrated vertical advection of moisture
106 ,omeg & ! omega (pa/s)
107 ,csum & ! used to implement memory, set to zero if not avail
108 ,cnvwt & ! gfs needs this
109 ,zuo & ! nomalized updraft mass flux
110 ,zdo & ! nomalized downdraft mass flux
111 ,zdm & ! nomalized downdraft mass flux from mid scheme
112 ,edto & !
113 ,edtm & !
114 ,xmb_out & ! the xmb's may be needed for dicycle
115 ,xmbm_in & !
116 ,xmbs_in & !
117 ,pre & !
118 ,outu & ! momentum tendencies at mass points
119 ,outv & !
120 ,outt & ! temperature tendencies
121 ,outq & ! q tendencies
122 ,outqc & ! ql/qice tendencies
123 ,kbcon & ! lfc of parcel from k22
124 ,ktop & ! cloud top
125 ,cupclw & ! used for direct coupling to radiation, but with tuning factors
126 ,frh_out & ! fractional coverage
127 ,ierr & ! ierr flags are error flags, used for debugging
128 ,ierrc & ! the following should be set to zero if not available
129 ,nchem &
130 ,fscav &
131 ,chem3d &
132 ,wetdpc_deep &
133 ,do_smoke_transport &
134 ,rand_mom & ! for stochastics mom, if temporal and spatial patterns exist
135 ,rand_vmas & ! for stochastics vertmass, if temporal and spatial patterns exist
136 ,rand_clos & ! for stochastics closures, if temporal and spatial patterns exist
137 ,nranflag & ! flag to what you want perturbed
138 !! 1 = momentum transport
139 !! 2 = normalized vertical mass flux profile
140 !! 3 = closures
141 !! more is possible, talk to developer or
142 !! implement yourself. pattern is expected to be
143 !! betwee -1 and +1
144 ,do_capsuppress,cap_suppress_j & !
145 ,k22 & !
146 ,jmin,kdt,tropics) !
147
148 implicit none
149
150 integer &
151 ,intent (in ) :: &
152 nranflag,itf,ktf,its,ite, kts,kte,ipr,imid,kdt
153 integer, intent (in ) :: &
154 ichoice,nchem
155 real(kind=kind_phys), dimension (its:ite,4) &
156 ,intent (in ) :: rand_clos
157 real(kind=kind_phys), dimension (its:ite) &
158 ,intent (in ) :: rand_mom,rand_vmas
159!$acc declare copyin(rand_clos,rand_mom,rand_vmas)
160
161 integer, intent(in) :: do_capsuppress
162 real(kind=kind_phys), intent(in), dimension(:), optional :: cap_suppress_j
163!$acc declare create(cap_suppress_j)
164 !
165 !
166 !
167 real(kind=kind_phys), dimension (its:ite,1:maxens3) :: xf_ens,pr_ens
168!$acc declare create(xf_ens,pr_ens)
169 ! outtem = output temp tendency (per s)
170 ! outq = output q tendency (per s)
171 ! outqc = output qc tendency (per s)
172 ! pre = output precip
173 real(kind=kind_phys), dimension (its:ite,kts:kte) &
174 ,intent (inout ) :: &
175 cnvwt,outu,outv,outt,outq,outqc,cupclw
176 real(kind=kind_phys), dimension (its:ite) &
177 ,intent (out ) :: &
178 frh_out
179 real(kind=kind_phys), dimension (its:ite) &
180 ,intent (inout ) :: &
181 pre,xmb_out
182!$acc declare copy(cnvwt,outu,outv,outt,outq,outqc,cupclw,frh_out,pre,xmb_out)
183 real(kind=kind_phys), dimension (its:ite) &
184 ,intent (in ) :: &
185 hfx,qfx,xmbm_in,xmbs_in
186!$acc declare copyin(hfx,qfx,xmbm_in,xmbs_in)
187 integer, dimension (its:ite) &
188 ,intent (inout ) :: &
189 kbcon,ktop
190!$acc declare copy(kbcon,ktop)
191 integer, dimension (its:ite) &
192 ,intent (in ) :: &
193 kpbl,tropics
194!$acc declare copyin(kpbl,tropics)
195 !
196 ! basic environmental input includes moisture convergence (mconv)
197 ! omega (omeg), windspeed (us,vs), and a flag (ierr) to turn off
198 ! convection for this call only and at that particular gridpoint
199 !
200 real(kind=kind_phys), dimension (its:ite,kts:kte) &
201 ,intent (in ) :: &
202 dhdt,rho,t,po,us,vs,tn
203!$acc declare copyin(dhdt,rho,t,po,us,vs,tn)
204 real(kind=kind_phys), dimension (its:ite,kts:kte) &
205 ,intent (inout ) :: &
206 omeg
207!$acc declare copy(omeg)
208 real(kind=kind_phys), dimension (its:ite,kts:kte) &
209 ,intent (inout) :: &
210 q,qo,zuo,zdo,zdm
211!$acc declare copy(q,qo,zuo,zdo,zdm)
212 real(kind=kind_phys), dimension (its:ite) &
213 ,intent (in ) :: &
214 dx,z1,psur,xland
215!$acc declare copyin(dx,z1,psur,xland)
216 real(kind=kind_phys), dimension (its:ite) &
217 ,intent (inout ) :: &
218 mconv,ccn
219!$acc declare copy(mconv,ccn)
220 real(kind=kind_phys), dimension (:,:,:) &
221 ,intent (inout), optional :: &
222 chem3d
223 logical, intent (in) :: do_smoke_transport
224 real(kind=kind_phys), dimension (:,:) &
225 , intent (out), optional :: wetdpc_deep
226 real(kind=kind_phys), intent (in) :: fscav(:)
227!$acc declare copy(chem3d) copyout(wetdpc_deep) copyin(fscav)
228
229 real(kind=kind_phys) &
230 ,intent (in ) :: &
231 dtime,ccnclean
232
233
234!
235! local ensemble dependent variables in this routine
236!
237 real(kind=kind_phys), dimension (its:ite,1) :: &
238 xaa0_ens
239 real(kind=kind_phys), dimension (its:ite,1) :: &
240 edtc
241 real(kind=kind_phys), dimension (its:ite,kts:kte,1) :: &
242 dellat_ens,dellaqc_ens,dellaq_ens,pwo_ens
243!$acc declare create(xaa0_ens,edtc,dellat_ens,dellaqc_ens,dellaq_ens,pwo_ens)
244!
245!
246!
247!***************** the following are your basic environmental
248! variables. they carry a "_cup" if they are
249! on model cloud levels (staggered). they carry
250! an "o"-ending (z becomes zo), if they are the forced
251! variables. they are preceded by x (z becomes xz)
252! to indicate modification by some typ of cloud
253!
254 ! z = heights of model levels
255 ! q = environmental mixing ratio
256 ! qes = environmental saturation mixing ratio
257 ! t = environmental temp
258 ! p = environmental pressure
259 ! he = environmental moist static energy
260 ! hes = environmental saturation moist static energy
261 ! z_cup = heights of model cloud levels
262 ! q_cup = environmental q on model cloud levels
263 ! qes_cup = saturation q on model cloud levels
264 ! t_cup = temperature (kelvin) on model cloud levels
265 ! p_cup = environmental pressure
266 ! he_cup = moist static energy on model cloud levels
267 ! hes_cup = saturation moist static energy on model cloud levels
268 ! gamma_cup = gamma on model cloud levels
269!
270!
271 ! hcd = moist static energy in downdraft
272 ! zd normalized downdraft mass flux
273 ! dby = buoancy term
274 ! entr = entrainment rate
275 ! zd = downdraft normalized mass flux
276 ! entr= entrainment rate
277 ! hcd = h in model cloud
278 ! bu = buoancy term
279 ! zd = normalized downdraft mass flux
280 ! gamma_cup = gamma on model cloud levels
281 ! qcd = cloud q (including liquid water) after entrainment
282 ! qrch = saturation q in cloud
283 ! pwd = evaporate at that level
284 ! pwev = total normalized integrated evaoprate (i2)
285 ! entr= entrainment rate
286 ! z1 = terrain elevation
287 ! entr = downdraft entrainment rate
288 ! jmin = downdraft originating level
289 ! kdet = level above ground where downdraft start detraining
290 ! psur = surface pressure
291 ! z1 = terrain elevation
292 ! pr_ens = precipitation ensemble
293 ! xf_ens = mass flux ensembles
294 ! omeg = omega from large scale model
295 ! mconv = moisture convergence from large scale model
296 ! zd = downdraft normalized mass flux
297 ! zu = updraft normalized mass flux
298 ! dir = "storm motion"
299 ! mbdt = arbitrary numerical parameter
300 ! dtime = dt over which forcing is applied
301 ! kbcon = lfc of parcel from k22
302 ! k22 = updraft originating level
303 ! ichoice = flag if only want one closure (usually set to zero!)
304 ! dby = buoancy term
305 ! ktop = cloud top (output)
306 ! xmb = total base mass flux
307 ! hc = cloud moist static energy
308 ! hkb = moist static energy at originating level
309 real(kind=kind_phys), dimension (its:ite,kts:kte,nchem) :: &
310 chem
311 real(kind=kind_phys), dimension (its:ite,kts:kte,nchem) :: &
312 chem_cup,chem_up,chem_down,dellac,dellac2,chem_c,chem_pw,chem_pwd
313 real(kind=kind_phys), dimension (its:ite,nchem) :: &
314 chem_pwav,chem_psum
315 real(kind=kind_phys):: dtime_max,sum1,sum2
316 real(kind=kind_phys), dimension (kts:kte) :: trac,trcflx_in,trcflx_out,trc,trco
317 real(kind=kind_phys), dimension (its:ite,kts:kte) :: pwdper, massflx
318 integer :: nv
319!$acc declare create(chem,chem_cup,chem_up,chem_down,dellac,dellac2,chem_c,chem_pw,chem_pwd, &
320!$acc chem_pwav,chem_psum,pwdper,massflx)
321
322 real(kind=kind_phys), dimension (its:ite,kts:kte) :: &
323 entr_rate_2d,mentrd_rate_2d,he,hes,qes,z, heo,heso,qeso,zo, &
324 xhe,xhes,xqes,xz,xt,xq,qes_cup,q_cup,he_cup,hes_cup,z_cup, &
325 p_cup,gamma_cup,t_cup, qeso_cup,qo_cup,heo_cup,heso_cup, &
326 zo_cup,po_cup,gammao_cup,tn_cup, &
327 xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup, &
328 xt_cup, dby,hc,zu,clw_all, &
329 dbyo,qco,qrcdo,pwdo,pwo,hcdo,qcdo,dbydo,hco,qrco, &
330 dbyt,xdby,xhc,xzu, &
331
332 ! cd = detrainment function for updraft
333 ! cdd = detrainment function for downdraft
334 ! dellat = change of temperature per unit mass flux of cloud ensemble
335 ! dellaq = change of q per unit mass flux of cloud ensemble
336 ! dellaqc = change of qc per unit mass flux of cloud ensemble
337
338 cd,cdd,dellah,dellaq,dellat,dellaqc, &
339 u_cup,v_cup,uc,vc,ucd,vcd,dellu,dellv
340!$acc declare create( &
341!$acc entr_rate_2d,mentrd_rate_2d,he,hes,qes,z, heo,heso,qeso,zo, &
342!$acc xhe,xhes,xqes,xz,xt,xq,qes_cup,q_cup,he_cup,hes_cup,z_cup, &
343!$acc p_cup,gamma_cup,t_cup, qeso_cup,qo_cup,heo_cup,heso_cup, &
344!$acc zo_cup,po_cup,gammao_cup,tn_cup, &
345!$acc xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup, &
346!$acc xt_cup, dby,hc,zu,clw_all, &
347!$acc dbyo,qco,qrcdo,pwdo,pwo,hcdo,qcdo,dbydo,hco,qrco, &
348!$acc dbyt,xdby,xhc,xzu, &
349!$acc cd,cdd,dellah,dellaq,dellat,dellaqc, &
350!$acc u_cup,v_cup,uc,vc,ucd,vcd,dellu,dellv)
351
352 ! aa0 cloud work function for downdraft
353 ! edt = epsilon
354 ! aa0 = cloud work function without forcing effects
355 ! aa1 = cloud work function with forcing effects
356 ! xaa0 = cloud work function with cloud effects (ensemble dependent)
357 ! edt = epsilon
358
359 real(kind=kind_phys), dimension (its:ite) :: &
360 edt,edto,edtm,aa1,aa0,xaa0,hkb, &
361 hkbo,xhkb, &
362 xmb,pwavo,ccnloss, &
363 pwevo,bu,bud,cap_max, &
364 cap_max_increment,closure_n,psum,psumh,sig,sigd
365 real(kind=kind_phys), dimension (its:ite) :: &
366 axx,edtmax,edtmin,entr_rate
367 integer, dimension (its:ite) :: &
368 kzdown,kdet,k22,jmin,kstabi,kstabm,k22x,xland1, &
369 ktopdby,kbconx,ierr2,ierr3,kbmax
370!$acc declare create(edt,edto,edtm,aa1,aa0,xaa0,hkb, &
371!$acc hkbo,xhkb, &
372!$acc xmb,pwavo,ccnloss, &
373!$acc pwevo,bu,bud,cap_max, &
374!$acc cap_max_increment,closure_n,psum,psumh,sig,sigd, &
375!$acc axx,edtmax,edtmin,entr_rate, &
376!$acc kzdown,kdet,k22,jmin,kstabi,kstabm,k22x,xland1, &
377!$acc ktopdby,kbconx,ierr2,ierr3,kbmax)
378
379 integer, dimension (its:ite), intent(inout) :: ierr
380 integer, dimension (its:ite), intent(in), optional :: csum
381!$acc declare copy(ierr) copyin(csum)
382 integer :: &
383 iloop,nens3,ki,kk,i,k
384 real(kind=kind_phys) :: &
385 dz,dzo,mbdt,radius, &
386 zcutdown,depth_min,zkbmax,z_detr,zktop, &
387 dh,cap_maxs,trash,trash2,frh,sig_thresh
388 real(kind=kind_phys), dimension (its:ite) :: pefc
389 real(kind=kind_phys) entdo,dp,subin,detdo,entup, &
390 detup,subdown,entdoj,entupk,detupk,totmas
391
392 real(kind=kind_phys), dimension (its:ite) :: lambau,flux_tun,zws,ztexec,zqexec
393!$acc declare create(lambau,flux_tun,zws,ztexec,zqexec)
394
395 integer :: jprnt,jmini,start_k22
396 logical :: keep_going,flg(its:ite)
397!$acc declare create(flg)
398
399 character*50 :: ierrc(its:ite)
400 character*4 :: cumulus
401 real(kind=kind_phys), dimension (its:ite,kts:kte) :: &
402 up_massentr,up_massdetr,c1d &
403 ,up_massentro,up_massdetro,dd_massentro,dd_massdetro
404 real(kind=kind_phys), dimension (its:ite,kts:kte) :: &
405 up_massentru,up_massdetru,dd_massentru,dd_massdetru
406!$acc declare create(up_massentr,up_massdetr,c1d,up_massentro,up_massdetro,dd_massentro,dd_massdetro, &
407!$acc up_massentru,up_massdetru,dd_massentru,dd_massdetru)
408 real(kind=kind_phys) c1_max,buo_flux,pgcon,pgc,blqe
409
410 real(kind=kind_phys) :: xff_mid(its:ite,2)
411!$acc declare create(xff_mid)
412 integer :: iversion=1
413 real(kind=kind_phys) :: denom,h_entr,umean,t_star,dq
414 integer, intent(in) :: dicycle
415 real(kind=kind_phys), dimension (its:ite) :: aa1_bl,hkbo_bl,tau_bl,tau_ecmwf,wmean
416 real(kind=kind_phys), dimension (its:ite,kts:kte) :: tn_bl, qo_bl, qeso_bl, heo_bl, heso_bl &
417 ,qeso_cup_bl,qo_cup_bl, heo_cup_bl,heso_cup_bl &
418 ,gammao_cup_bl,tn_cup_bl,hco_bl,dbyo_bl
419 real(kind=kind_phys), dimension(its:ite) :: xf_dicycle
420!$acc declare create(aa1_bl,hkbo_bl,tau_bl,tau_ecmwf,wmean, &
421!$acc tn_bl, qo_bl, qeso_bl, heo_bl, heso_bl, &
422!$acc qeso_cup_bl,qo_cup_bl, heo_cup_bl,heso_cup_bl, &
423!$acc gammao_cup_bl,tn_cup_bl,hco_bl,dbyo_bl,xf_dicycle)
424 real(kind=kind_phys), intent(inout), dimension(its:ite,10) :: forcing
425!$acc declare copy(forcing)
426 integer :: turn,pmin_lev(its:ite),start_level(its:ite),ktopkeep(its:ite)
427 real(kind=kind_phys), dimension (its:ite,kts:kte) :: dtempdz
428 integer, dimension (its:ite,kts:kte) :: k_inv_layers
429 real(kind=kind_phys), dimension (its:ite) :: c0 ! HCB
430 real(kind=kind_phys), dimension (its:ite,kts:kte) :: c0t3d ! hli for smoke/dust wet scavenging
431!$acc declare create(pmin_lev,start_level,ktopkeep,dtempdz,k_inv_layers,c0,c0t3d)
432
433! rainevap from sas
434 real(kind=kind_phys) zuh2(40)
435 real(kind=kind_phys), dimension (its:ite) :: rntot,delqev,delq2,qevap,rn,qcond
436!$acc declare create(zuh2,rntot,delqev,delq2,qevap,rn,qcond)
437 real(kind=kind_phys) :: rain,t1,q1,elocp,evef,el2orc,evfact,evfactl,g_rain,e_dn,c_up
438 real(kind=kind_phys) :: pgeoh,dts,fp,fpi,pmin,x_add,beta,beta_u
439 real(kind=kind_phys) :: cbeg,cmid,cend,const_a,const_b,const_c
440!---meltglac-------------------------------------------------
441
442 real(kind=kind_phys), dimension (its:ite,kts:kte) :: p_liq_ice,melting_layer,melting
443!$acc declare create(p_liq_ice,melting_layer,melting)
444
445 integer :: itemp
446
447!---meltglac-------------------------------------------------
448!$acc kernels
449 melting_layer(:,:)=0.
450 melting(:,:)=0.
451 flux_tun(:)=fluxtune
452!$acc end kernels
453! if(imid.eq.1)flux_tun(:)=fluxtune+.5
454 cumulus='deep'
455 if(imid.eq.1)cumulus='mid'
456 pmin=150.
457 if(imid.eq.1)pmin=75.
458!$acc kernels
459 ktopdby(:)=0
460!$acc end kernels
461 c1_max=c1
462 elocp=xlv/cp
463 el2orc=xlv*xlv/(r_v*cp)
464 evfact=0.25 ! .4
465 evfactl=0.25 ! .2
466
467
468!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
469!
470!proportionality constant to estimate pressure gradient of updraft (zhang and wu, 2003, jas
471!
472! ecmwf
473 pgcon=0.
474!$acc kernels
475 lambau(:)=2.0
476 if(imid.eq.1)lambau(:)=2.0
477! here random must be between -1 and 1
478 if(nranflag == 1)then
479 lambau(:)=1.5+rand_mom(:)
480 endif
481!$acc end kernels
482! sas
483! lambau=0.
484! pgcon=-.55
485!
486!---------------------------------------------------- ! HCB
487! Set cloud water to rain water conversion rate (c0)
488!$acc kernels
489 c0(:)=0.004
490 do i=its,itf
491 xland1(i)=int(xland(i)+.0001) ! 1.
492 if(xland(i).gt.1.5 .or. xland(i).lt.0.5)then
493 xland1(i)=0
494 endif
495 if(xland1(i).eq.1)c0(i)=0.002
496 if(imid.eq.1)then
497 c0(i)=0.002
498 endif
499 enddo
500!$acc end kernels
501
502!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
503!$acc kernels
504 ztexec(:) = 0.
505 zqexec(:) = 0.
506 zws(:) = 0.
507
508 do i=its,itf
509 !- buoyancy flux (h+le)
510 buo_flux= (hfx(i)/cp+0.608*t(i,1)*qfx(i)/xlv)/rho(i,1)
511 pgeoh = zo(i,2)*g
512 !-convective-scale velocity w*
513 zws(i) = max(0.,flux_tun(i)*0.41*buo_flux*zo(i,2)*g/t(i,1))
514 if(zws(i) > tiny(pgeoh)) then
515 !-convective-scale velocity w*
516 zws(i) = 1.2*zws(i)**.3333
517 !- temperature excess
518 ztexec(i) = max(flux_tun(i)*hfx(i)/(rho(i,1)*zws(i)*cp),0.0)
519 !- moisture excess
520 zqexec(i) = max(flux_tun(i)*qfx(i)/xlv/(rho(i,1)*zws(i)),0.)
521 endif
522 !- zws for shallow convection closure (grant 2001)
523 !- height of the pbl
524 zws(i) = max(0.,.001-flux_tun(i)*0.41*buo_flux*zo(i,kpbl(i))*g/t(i,kpbl(i)))
525 zws(i) = 1.2*zws(i)**.3333
526 zws(i) = zws(i)*rho(i,kpbl(i)) !check if zrho is correct
527 enddo
528!$acc end kernels
529 cap_maxs=75. ! 150.
530!$acc kernels
531 do i=its,itf
532 edto(i)=0.
533 closure_n(i)=16.
534 xmb_out(i)=0.
535 cap_max(i)=cap_maxs
536 cap_max_increment(i)=20.
537!
538! for water or ice
539!
540 if (xland1(i)==0) then
541 cap_max_increment(i)=20.
542 else
543 if(ztexec(i).gt.0.)cap_max(i)=cap_max(i)+25.
544 if(ztexec(i).lt.0.)cap_max(i)=cap_max(i)-25.
545 endif
546#ifndef _OPENACC
547 ierrc(i)=" "
548#endif
549 enddo
550!$acc end kernels
551 if(use_excess == 0 )then
552!$acc kernels
553 ztexec(:)=0
554 zqexec(:)=0
555!$acc end kernels
556 endif
557 if(do_capsuppress == 1) then
558!$acc kernels
559 do i=its,itf
560 cap_max(i)=cap_maxs
561 if (abs(cap_suppress_j(i) - 1.0 ) < 0.1 ) then
562 cap_max(i)=cap_maxs+75.
563 elseif (abs(cap_suppress_j(i) - 0.0 ) < 0.1 ) then
564 cap_max(i)=10.0
565 endif
566 enddo
567!$acc end kernels
568 endif
569!
570!--- initial entrainment rate (these may be changed later on in the
571!--- program
572!
573!$acc kernels
574 start_level(:)=kte
575!$acc end kernels
576
577!$acc kernels
578!$acc loop private(radius,frh)
579 do i=its,ite
580 c1d(i,:)= 0. !c1 ! 0. ! c1 ! max(.003,c1+float(csum(i))*.0001)
581 entr_rate(i)=7.e-5 - min(20.,float(csum(i))) * 3.e-6
582 if(xland1(i) == 0)entr_rate(i)=7.e-5
583 if(dx(i)<dx_thresh) entr_rate(i)=2.e-4
584 if(imid.eq.1)entr_rate(i)=3.e-4
585 radius=.2/entr_rate(i)
586 frh=min(1.,3.14*radius*radius/dx(i)/dx(i))
587 if(frh > frh_thresh)then
588 frh=frh_thresh
589 radius=sqrt(frh*dx(i)*dx(i)/3.14)
590 entr_rate(i)=.2/radius
591 endif
592 sig(i)=(1.-frh)**2
593 !frh_out(i) = frh
594 if(forcing(i,7).eq.0.)sig(i)=1.
595 if(kdt.le.(3600./dtime))sig(i)=1.
596 frh_out(i) = frh*sig(i)
597 enddo
598!$acc end kernels
599 sig_thresh = (1.-frh_thresh)**2
600
601
602!
603!--- entrainment of mass
604!
605!
606!--- initial detrainmentrates
607!
608!$acc kernels
609 do k=kts,ktf
610 do i=its,itf
611 cnvwt(i,k)=0.
612 zuo(i,k)=0.
613 zdo(i,k)=0.
614 z(i,k)=zo(i,k)
615 xz(i,k)=zo(i,k)
616 cupclw(i,k)=0.
617 cd(i,k)=.1*entr_rate(i) !1.e-9 ! 1.*entr_rate
618 if(imid.eq.1)cd(i,k)=.5*entr_rate(i)
619 cdd(i,k)=1.e-9
620 hcdo(i,k)=0.
621 qrcdo(i,k)=0.
622 dellaqc(i,k)=0.
623 enddo
624 enddo
625!$acc end kernels
626!
627!--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
628! base mass flux
629!
630!$acc kernels
631 edtmax(:)=1.
632! if(imid.eq.1)edtmax(:)=.15
633 edtmin(:)=.1
634! if(imid.eq.1)edtmin(:)=.05
635!$acc end kernels
636!
637!--- minimum depth (m), clouds must have
638!
639 depth_min=3000.
640!--- for RRFS allow only very deep convection
641 if(dx(its)<dx_thresh)depth_min=5000.
642 if(imid.eq.1)depth_min=2500.
643!
644!--- maximum depth (mb) of capping
645!--- inversion (larger cap = no convection)
646!
647!$acc kernels
648 do i=its,itf
649 kbmax(i)=1
650 aa0(i)=0.
651 aa1(i)=0.
652 edt(i)=0.
653 kstabm(i)=ktf-1
654 ierr2(i)=0
655 ierr3(i)=0
656 enddo
657!$acc end kernels
658 x_add=0.
659!
660!--- max height(m) above ground where updraft air can originate
661!
662 zkbmax=4000.
663 if(imid.eq.1)zkbmax=2000.
664!
665!--- height(m) above which no downdrafts are allowed to originate
666!
667 zcutdown=4000.
668!
669!--- depth(m) over which downdraft detrains all its mass
670!
671 z_detr=500.
672!
673
674!
675!--- environmental conditions, first heights
676!
677!$acc kernels
678 do i=its,itf
679 do k=1,maxens3
680 xf_ens(i,k)=0.
681 pr_ens(i,k)=0.
682 enddo
683 enddo
684!$acc end kernels
685!
687!
688 call cup_env(z,qes,he,hes,t,q,po,z1, &
689 psur,ierr,tcrit,-1, &
690 itf,ktf, &
691 its,ite, kts,kte)
692 call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
693 psur,ierr,tcrit,-1, &
694 itf,ktf, &
695 its,ite, kts,kte)
696
697!
699!
700 call cup_env_clev(t,qes,q,he,hes,z,po,qes_cup,q_cup,he_cup, &
701 hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
702 ierr,z1, &
703 itf,ktf, &
704 its,ite, kts,kte)
705 call cup_env_clev(tn,qeso,qo,heo,heso,zo,po,qeso_cup,qo_cup, &
706 heo_cup,heso_cup,zo_cup,po_cup,gammao_cup,tn_cup,psur, &
707 ierr,z1, &
708 itf,ktf, &
709 its,ite, kts,kte)
710!---meltglac-------------------------------------------------
712 call get_partition_liq_ice(ierr,tn,po_cup,p_liq_ice,melting_layer,&
713 itf,ktf,its,ite,kts,kte,cumulus)
714!---meltglac-------------------------------------------------
715!$acc kernels
716 do i=its,itf
717 if(ierr(i).eq.0)then
718 if(kpbl(i).gt.5 .and. imid.eq.1)cap_max(i)=po_cup(i,kpbl(i))
719 u_cup(i,kts)=us(i,kts)
720 v_cup(i,kts)=vs(i,kts)
721 do k=kts+1,ktf
722 u_cup(i,k)=.5*(us(i,k-1)+us(i,k))
723 v_cup(i,k)=.5*(vs(i,k-1)+vs(i,k))
724 enddo
725 endif
726 enddo
727 do i=its,itf
728 if(ierr(i).eq.0)then
729 do k=kts,ktf
730 if(zo_cup(i,k).gt.zkbmax+z1(i))then
731 kbmax(i)=k
732 go to 25
733 endif
734 enddo
735 25 continue
736!
738!
739 do k=kts,ktf
740 if(zo_cup(i,k).gt.z_detr+z1(i))then
741 kdet(i)=k
742 go to 26
743 endif
744 enddo
745 26 continue
746!
747 endif
748 enddo
749!$acc end kernels
750
751!
752!
753!
755!
756 start_k22=2
757!$acc parallel loop
758 do 36 i=its,itf
759 if(ierr(i).eq.0)then
760 k22(i)=maxloc(heo_cup(i,start_k22:kbmax(i)+2),1)+start_k22-1
761 if(k22(i).ge.kbmax(i))then
762 ierr(i)=2
763#ifndef _OPENACC
764 ierrc(i)="could not find k22"
765#endif
766 ktop(i)=0
767 k22(i)=0
768 kbcon(i)=0
769 endif
770 endif
771 36 continue
772!$acc end parallel
773
774!
777!
778!$acc parallel loop private(x_add)
779 do i=its,itf
780 if(ierr(i).eq.0)then
781 x_add = xlv*zqexec(i)+cp*ztexec(i)
782 call get_cloud_bc(kte,he_cup(i,1:kte),hkb(i),k22(i),x_add)
783 call get_cloud_bc(kte,heo_cup(i,1:kte),hkbo(i),k22(i),x_add)
784 endif ! ierr
785 enddo
786!$acc end parallel
787
788 jprnt=0
789 iloop=1
790 if(imid.eq.1)iloop=5
791 call cup_kbcon(ierrc,cap_max_increment,iloop,k22,kbcon,heo_cup,heso_cup, &
792 hkbo,ierr,kbmax,po_cup,cap_max, &
793 ztexec,zqexec, &
794 jprnt,itf,ktf, &
795 its,ite, kts,kte, &
796 z_cup,entr_rate,heo,imid)
797!
799!
800 call cup_minimi(heso_cup,kbcon,kstabm,kstabi,ierr, &
801 itf,ktf, &
802 its,ite, kts,kte)
803!$acc parallel loop private(frh,x_add)
804 do i=its,itf
805 if(ierr(i) == 0)then
806 frh = min(qo_cup(i,kbcon(i))/qeso_cup(i,kbcon(i)),1.)
807 if(frh >= rh_thresh .and. sig(i) <= sig_thresh )then
808 ierr(i)=231
809 cycle
810 endif
811!
812! never go too low...
813!
814 x_add=0.
815!$acc loop seq
816 do k=kbcon(i)+1,ktf
817 if(po(i,kbcon(i))-po(i,k) > pmin+x_add)then
818 pmin_lev(i)=k
819 exit
820 endif
821 enddo
822!
824!
825 start_level(i)=k22(i)
826 x_add = xlv*zqexec(i)+cp*ztexec(i)
827 call get_cloud_bc(kte,he_cup(i,1:kte),hkb(i),k22(i),x_add)
828 endif
829 enddo
830!$acc end parallel
831
832!
834!
835 if(imid.eq.1)then
836 call get_inversion_layers(ierr,p_cup,t_cup,z_cup,q_cup,qes_cup,k_inv_layers, &
837 kbcon,kstabi,dtempdz,itf,ktf,its,ite, kts,kte)
838 endif
839!$acc kernels
840 do i=its,itf
841 if(kstabi(i).lt.kbcon(i))then
842 kbcon(i)=1
843 ierr(i)=42
844 endif
845 do k=kts,ktf
846 entr_rate_2d(i,k)=entr_rate(i)
847 enddo
848 if(ierr(i).eq.0)then
849 kbcon(i)=max(2,kbcon(i))
850 do k=kts+1,ktf
851 frh = min(qo_cup(i,k)/qeso_cup(i,k),1.)
852 entr_rate_2d(i,k)=entr_rate(i) *(1.3-frh)
853 enddo
854 if(imid.eq.1)then
855 if(k_inv_layers(i,2).gt.0 .and. &
856 (po_cup(i,k22(i))-po_cup(i,k_inv_layers(i,2))).lt.500.)then
857
858 ktop(i)=min(kstabi(i),k_inv_layers(i,2))
859 ktopdby(i)=ktop(i)
860 else
861!$acc loop seq
862 do k=kbcon(i)+1,ktf
863 if((po_cup(i,k22(i))-po_cup(i,k)).gt.500.)then
864 ktop(i)=k
865 ktopdby(i)=ktop(i)
866 exit
867 endif
868 enddo
869 endif ! k_inv_lay
870 endif
871
872 endif
873 enddo
874!$acc end kernels
875
876!
878!
879 i=0
880 !- for mid level clouds we do not allow clouds taller than where stability
881 !- changes
882 if(imid.eq.1)then
883 call rates_up_pdf(rand_vmas,ipr,'mid',ktop,ierr,po_cup,entr_rate_2d,hkbo,heo,heso_cup,zo_cup, &
884 xland1,kstabi,k22,kbcon,its,ite,itf,kts,kte,ktf,zuo,kpbl,ktopdby,csum,pmin_lev)
885 else
886 call rates_up_pdf(rand_vmas,ipr,'deep',ktop,ierr,po_cup,entr_rate_2d,hkbo,heo,heso_cup,zo_cup, &
887 xland1,kstabi,k22,kbcon,its,ite,itf,kts,kte,ktf,zuo,kbcon,ktopdby,csum,pmin_lev)
888 endif
889!
890!
891!
892!$acc kernels
893 do i=its,itf
894 if(ierr(i).eq.0)then
895
896 if(k22(i).gt.1)then
897!$acc loop independent
898 do k=1,k22(i) -1
899 zuo(i,k)=0.
900 zu(i,k)=0.
901 xzu(i,k)=0.
902 enddo
903 endif
904!$acc loop independent
905 do k=k22(i),ktop(i)
906 xzu(i,k)= zuo(i,k)
907 zu(i,k)= zuo(i,k)
908 enddo
909!$acc loop independent
910 do k=ktop(i)+1,kte
911 zuo(i,k)=0.
912 zu(i,k)=0.
913 xzu(i,k)=0.
914 enddo
915 endif
916 enddo
917!$acc end kernels
918!
920!
921 if(imid.eq.1)then
922 call get_lateral_massflux(itf,ktf, its,ite, kts,kte &
923 ,ierr,ktop,zo_cup,zuo,cd,entr_rate_2d &
924 ,up_massentro, up_massdetro ,up_massentr, up_massdetr &
925 ,3,kbcon,k22,up_massentru,up_massdetru,lambau)
926 else
927 call get_lateral_massflux(itf,ktf, its,ite, kts,kte &
928 ,ierr,ktop,zo_cup,zuo,cd,entr_rate_2d &
929 ,up_massentro, up_massdetro ,up_massentr, up_massdetr &
930 ,1,kbcon,k22,up_massentru,up_massdetru,lambau)
931 endif
932
933
934!
935! note: ktop here already includes overshooting, ktopdby is without
936! overshooting
937!
938!$acc kernels
939 do k=kts,ktf
940 do i=its,itf
941 uc(i,k)=0.
942 vc(i,k)=0.
943 hc(i,k)=0.
944 dby(i,k)=0.
945 hco(i,k)=0.
946 dbyo(i,k)=0.
947 enddo
948 enddo
949 do i=its,itf
950 if(ierr(i).eq.0)then
951 do k=1,start_level(i)
952 uc(i,k)=u_cup(i,k)
953 vc(i,k)=v_cup(i,k)
954 enddo
955 do k=1,start_level(i)-1
956 hc(i,k)=he_cup(i,k)
957 hco(i,k)=heo_cup(i,k)
958 enddo
959 k=start_level(i)
960 hc(i,k)=hkb(i)
961 hco(i,k)=hkbo(i)
962 endif
963 enddo
964!$acc end kernels
965!
966!---meltglac-------------------------------------------------
967 !
968 !--- 1st guess for moist static energy and dbyo (not including ice phase)
969 !
970!$acc parallel loop private(denom,kk,ki)
971 do i=its,itf
972 ktopkeep(i)=0
973 dbyt(i,:)=0.
974 if(ierr(i) /= 0) cycle
975 ktopkeep(i)=ktop(i)
976!$acc loop seq
977 do k=start_level(i) +1,ktop(i) !mass cons option
978
979 denom=zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1)
980 if(denom.lt.1.e-8)then
981 ierr(i)=51
982 exit
983 endif
984 hco(i,k)=(hco(i,k-1)*zuo(i,k-1)-.5*up_massdetro(i,k-1)*hco(i,k-1)+ &
985 up_massentro(i,k-1)*heo(i,k-1)) / &
986 (zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
987 dbyo(i,k)=hco(i,k)-heso_cup(i,k)
988 enddo
989 ! for now no overshooting (only very little)
990!$acc loop seq
991 do k=ktop(i)-1,kbcon(i),-1
992 if(dbyo(i,k).gt.0.)then
993 ktopkeep(i)=k+1
994 exit
995 endif
996 enddo
997 enddo
998!$acc end parallel
999
1000!$acc kernels
1001 do 37 i=its,itf
1002 kzdown(i)=0
1003 if(ierr(i).eq.0)then
1004 zktop=(zo_cup(i,ktop(i))-z1(i))*.6
1005 if(imid.eq.1)zktop=(zo_cup(i,ktop(i))-z1(i))*.4
1006 zktop=min(zktop+z1(i),zcutdown+z1(i))
1007!$acc loop seq
1008 do k=kts,ktf
1009 if(zo_cup(i,k).gt.zktop)then
1010 kzdown(i)=k
1011 kzdown(i)=min(kzdown(i),kstabi(i)-1) !
1012 go to 37
1013 endif
1014 enddo
1015 endif
1016 37 continue
1017!$acc end kernels
1018
1019!
1021!
1022 call cup_minimi(heso_cup,k22,kzdown,jmin,ierr, &
1023 itf,ktf, &
1024 its,ite, kts,kte)
1025!$acc kernels
1026 do 100 i=its,itf
1027 if(ierr(i).eq.0)then
1028!
1029!-----srf-08aug2017-----begin
1030! if(imid .ne. 1 .and. melt_glac) jmin(i)=max(jmin(i),maxloc(melting_layer(i,:),1))
1031!-----srf-08aug2017-----end
1032
1033!--- check whether it would have buoyancy, if there where
1034!--- no entrainment/detrainment
1035!
1036 jmini = jmin(i)
1037 keep_going = .true.
1038 do while ( keep_going )
1039 keep_going = .false.
1040 if ( jmini - 1 .lt. kdet(i) ) kdet(i) = jmini-1
1041 if ( jmini .ge. ktop(i)-1 ) jmini = ktop(i) - 2
1042 ki = jmini
1043 hcdo(i,ki)=heso_cup(i,ki)
1044 dz=zo_cup(i,ki+1)-zo_cup(i,ki)
1045 dh=0.
1046!$acc loop seq
1047 do k=ki-1,1,-1
1048 hcdo(i,k)=heso_cup(i,jmini)
1049 dz=zo_cup(i,k+1)-zo_cup(i,k)
1050 dh=dh+dz*(hcdo(i,k)-heso_cup(i,k))
1051 if(dh.gt.0.)then
1052 jmini=jmini-1
1053 if ( jmini .gt. 5 ) then
1054 keep_going = .true.
1055 else
1056 ierr(i) = 9
1057#ifndef _OPENACC
1058 ierrc(i) = "could not find jmini9"
1059#endif
1060 exit
1061 endif
1062 endif
1063 enddo
1064 enddo
1065 jmin(i) = jmini
1066 if ( jmini .le. 5 ) then
1067 ierr(i)=4
1068#ifndef _OPENACC
1069 ierrc(i) = "could not find jmini4"
1070#endif
1071 endif
1072 endif
1073100 continue
1074 do i=its,itf
1075 if(ierr(i) /= 0) cycle
1076!$acc loop independent
1077 do k=ktop(i)+1,ktf
1078 hco(i,k)=heso_cup(i,k)
1079 dbyo(i,k)=0.
1080 enddo
1081 enddo
1082!$acc end kernels
1083 !
1085 !
1086 if(imid.eq.1)then
1087 call cup_up_moisture('mid',ierr,zo_cup,qco,qrco,pwo,pwavo, &
1088 p_cup,kbcon,ktop,dbyo,clw_all,xland1, &
1089 qo,gammao_cup,zuo,qeso_cup,k22,qo_cup,c0,c0t3d, &
1090 zqexec,ccn,ccnclean,rho,c1d,tn_cup,autoconv,up_massentr,up_massdetr,psum,psumh, &
1091 1,itf,ktf, &
1092 its,ite, kts,kte)
1093 else
1094 call cup_up_moisture('deep',ierr,zo_cup,qco,qrco,pwo,pwavo, &
1095 p_cup,kbcon,ktop,dbyo,clw_all,xland1, &
1096 qo,gammao_cup,zuo,qeso_cup,k22,qo_cup,c0,c0t3d, &
1097 zqexec,ccn,ccnclean,rho,c1d,tn_cup,autoconv,up_massentr,up_massdetr,psum,psumh, &
1098 1,itf,ktf, &
1099 its,ite, kts,kte)
1100 endif
1101!---meltglac-------------------------------------------------
1102
1103!$acc kernels
1104 do i=its,itf
1105
1106 ktopkeep(i)=0
1107 dbyt(i,:)=0.
1108 if(ierr(i) /= 0) cycle
1109 ktopkeep(i)=ktop(i)
1110!$acc loop seq
1111 do k=start_level(i) +1,ktop(i) !mass cons option
1112
1113 denom=zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1)
1114 if(denom.lt.1.e-8)then
1115 ierr(i)=51
1116 exit
1117 endif
1118
1119 hc(i,k)=(hc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*hc(i,k-1)+ &
1120 up_massentr(i,k-1)*he(i,k-1)) / &
1121 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
1122 uc(i,k)=(uc(i,k-1)*zu(i,k-1)-.5*up_massdetru(i,k-1)*uc(i,k-1)+ &
1123 up_massentru(i,k-1)*us(i,k-1) &
1124 -pgcon*.5*(zu(i,k)+zu(i,k-1))*(u_cup(i,k)-u_cup(i,k-1))) / &
1125 (zu(i,k-1)-.5*up_massdetru(i,k-1)+up_massentru(i,k-1))
1126 vc(i,k)=(vc(i,k-1)*zu(i,k-1)-.5*up_massdetru(i,k-1)*vc(i,k-1)+ &
1127 up_massentru(i,k-1)*vs(i,k-1) &
1128 -pgcon*.5*(zu(i,k)+zu(i,k-1))*(v_cup(i,k)-v_cup(i,k-1))) / &
1129 (zu(i,k-1)-.5*up_massdetru(i,k-1)+up_massentru(i,k-1))
1130 dby(i,k)=hc(i,k)-hes_cup(i,k)
1131 hco(i,k)=(hco(i,k-1)*zuo(i,k-1)-.5*up_massdetro(i,k-1)*hco(i,k-1)+ &
1132 up_massentro(i,k-1)*heo(i,k-1)) / &
1133 (zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
1134!---meltglac-------------------------------------------------
1135 !
1136 !- include glaciation effects on hc,hco
1137 ! ------ ice content --------
1138 hc(i,k)= hc(i,k)+(1.-p_liq_ice(i,k))*qrco(i,k)*xlf
1139 hco(i,k)= hco(i,k)+(1.-p_liq_ice(i,k))*qrco(i,k)*xlf
1140
1141 dby(i,k)=hc(i,k)-hes_cup(i,k)
1142!---meltglac-------------------------------------------------
1143 dbyo(i,k)=hco(i,k)-heso_cup(i,k)
1144 dz=zo_cup(i,k+1)-zo_cup(i,k)
1145 dbyt(i,k)=dbyt(i,k-1)+dbyo(i,k)*dz
1146
1147 enddo
1148! for now no overshooting (only very little)
1149 kk=maxloc(dbyt(i,:),1)
1150 ki=maxloc(zuo(i,:),1)
1151!
1152!$acc loop seq
1153 do k=ktop(i)-1,kbcon(i),-1
1154 if(dbyo(i,k).gt.0.)then
1155 ktopkeep(i)=k+1
1156 exit
1157 endif
1158 enddo
1159 enddo
1160!$acc end kernels
1161
116241 continue
1163!$acc kernels
1164 do i=its,itf
1165 if(ierr(i) /= 0) cycle
1166 do k=ktop(i)+1,ktf
1167 hc(i,k)=hes_cup(i,k)
1168 uc(i,k)=u_cup(i,k)
1169 vc(i,k)=v_cup(i,k)
1170 hco(i,k)=heso_cup(i,k)
1171 dby(i,k)=0.
1172 dbyo(i,k)=0.
1173 zu(i,k)=0.
1174 zuo(i,k)=0.
1175 cd(i,k)=0.
1176 entr_rate_2d(i,k)=0.
1177 up_massentr(i,k)=0.
1178 up_massdetr(i,k)=0.
1179 up_massentro(i,k)=0.
1180 up_massdetro(i,k)=0.
1181 enddo
1182 enddo
1183!
1184 do i=its,itf
1185 if(ierr(i)/=0)cycle
1186 if(ktop(i).lt.kbcon(i)+2)then
1187 ierr(i)=5
1188#ifndef _OPENACC
1189 ierrc(i)='ktop too small deep'
1190#endif
1191 ktop(i)=0
1192 endif
1193 enddo
1194!$acc end kernels
1195!
1196! - must have at least depth_min m between cloud convective base
1197! and cloud top.
1198!
1199!$acc kernels
1200 do i=its,itf
1201 if(ierr(i).eq.0)then
1202 if ( jmin(i) - 1 .lt. kdet(i) ) kdet(i) = jmin(i)-1
1203 if(-zo_cup(i,kbcon(i))+zo_cup(i,ktop(i)).lt.depth_min)then
1204 ierr(i)=6
1205#ifndef _OPENACC
1206 ierrc(i)="cloud depth very shallow"
1207#endif
1208 endif
1209 endif
1210 enddo
1211!$acc end kernels
1212
1213!
1214!--- normalized downdraft mass flux profile,also work on bottom detrainment
1215!--- in this routine
1216!
1217!$acc kernels
1218 do k=kts,ktf
1219 do i=its,itf
1220 zdo(i,k)=0.
1221 cdd(i,k)=0.
1222 dd_massentro(i,k)=0.
1223 dd_massdetro(i,k)=0.
1224 dd_massentru(i,k)=0.
1225 dd_massdetru(i,k)=0.
1226 hcdo(i,k)=heso_cup(i,k)
1227 ucd(i,k)=u_cup(i,k)
1228 vcd(i,k)=v_cup(i,k)
1229 dbydo(i,k)=0.
1230 mentrd_rate_2d(i,k)=entr_rate(i)
1231 enddo
1232 enddo
1233!$acc end kernels
1234
1235!$acc parallel loop private(beta,itemp,dzo,h_entr)
1236 do i=its,itf
1237 if(ierr(i)/=0)cycle
1238 beta=max(.025,.055-float(csum(i))*.0015) !.02
1239 if(imid.eq.1)beta=.025
1240 bud(i)=0.
1241 cdd(i,1:jmin(i))=.1*entr_rate(i)
1242 cdd(i,jmin(i))=0.
1243 dd_massdetro(i,:)=0.
1244 dd_massentro(i,:)=0.
1245 call get_zu_zd_pdf_fim(0,po_cup(i,:),rand_vmas(i),0.0_kind_phys,ipr,xland1(i),zuh2,4, &
1246 ierr(i),kdet(i),jmin(i)+1,zdo(i,:),kts,kte,ktf,beta,kpbl(i),csum(i),pmin_lev(i))
1247 if(zdo(i,jmin(i)) .lt.1.e-8)then
1248 zdo(i,jmin(i))=0.
1249 jmin(i)=jmin(i)-1
1250 cdd(i,jmin(i):ktf)=0.
1251 zdo(i,jmin(i)+1:ktf)=0.
1252 if(zdo(i,jmin(i)) .lt.1.e-8)then
1253 ierr(i)=876
1254 cycle
1255 endif
1256 endif
1257
1258 itemp = maxloc(zdo(i,:),1)
1259 do ki=jmin(i) , itemp,-1
1260 !=> from jmin to maximum value zd -> change entrainment
1261 dzo=zo_cup(i,ki+1)-zo_cup(i,ki)
1262 dd_massdetro(i,ki)=cdd(i,ki)*dzo*zdo(i,ki+1)
1263 dd_massentro(i,ki)=zdo(i,ki)-zdo(i,ki+1)+dd_massdetro(i,ki)
1264 if(dd_massentro(i,ki).lt.0.)then
1265 dd_massentro(i,ki)=0.
1266 dd_massdetro(i,ki)=zdo(i,ki+1)-zdo(i,ki)
1267 if(zdo(i,ki+1).gt.0.)cdd(i,ki)=dd_massdetro(i,ki)/(dzo*zdo(i,ki+1))
1268 endif
1269 if(zdo(i,ki+1).gt.0.)mentrd_rate_2d(i,ki)=dd_massentro(i,ki)/(dzo*zdo(i,ki+1))
1270 enddo
1271 mentrd_rate_2d(i,1)=0.
1272 do ki=itemp-1,1,-1
1273 !=> from maximum value zd to surface -> change detrainment
1274 dzo=zo_cup(i,ki+1)-zo_cup(i,ki)
1275 dd_massentro(i,ki)=mentrd_rate_2d(i,ki)*dzo*zdo(i,ki+1)
1276 dd_massdetro(i,ki) = zdo(i,ki+1)+dd_massentro(i,ki)-zdo(i,ki)
1277 if(dd_massdetro(i,ki).lt.0.)then
1278 dd_massdetro(i,ki)=0.
1279 dd_massentro(i,ki)=zdo(i,ki)-zdo(i,ki+1)
1280 if(zdo(i,ki+1).gt.0.)mentrd_rate_2d(i,ki)=dd_massentro(i,ki)/(dzo*zdo(i,ki+1))
1281 endif
1282 if(zdo(i,ki+1).gt.0.)cdd(i,ki)= dd_massdetro(i,ki)/(dzo*zdo(i,ki+1))
1283 enddo
1284!
1286 do k=2,jmin(i)+1
1287 dd_massentru(i,k-1)=dd_massentro(i,k-1)+lambau(i)*dd_massdetro(i,k-1)
1288 dd_massdetru(i,k-1)=dd_massdetro(i,k-1)+lambau(i)*dd_massdetro(i,k-1)
1289 enddo
1290 dbydo(i,jmin(i))=hcdo(i,jmin(i))-heso_cup(i,jmin(i))
1291 bud(i)=dbydo(i,jmin(i))*(zo_cup(i,jmin(i)+1)-zo_cup(i,jmin(i)))
1292 ucd(i,jmin(i)+1)=.5*(uc(i,jmin(i)+1)+u_cup(i,jmin(i)+1))
1293!$acc loop seq
1294 do ki=jmin(i) ,1,-1
1295 dzo=zo_cup(i,ki+1)-zo_cup(i,ki)
1296 h_entr=.5*(heo(i,ki)+.5*(hco(i,ki)+hco(i,ki+1)))
1297 ucd(i,ki)=(ucd(i,ki+1)*zdo(i,ki+1) &
1298 -.5*dd_massdetru(i,ki)*ucd(i,ki+1)+ &
1299 dd_massentru(i,ki)*us(i,ki) &
1300 -pgcon*zdo(i,ki+1)*(us(i,ki+1)-us(i,ki))) / &
1301 (zdo(i,ki+1)-.5*dd_massdetru(i,ki)+dd_massentru(i,ki))
1302 vcd(i,ki)=(vcd(i,ki+1)*zdo(i,ki+1) &
1303 -.5*dd_massdetru(i,ki)*vcd(i,ki+1)+ &
1304 dd_massentru(i,ki)*vs(i,ki) &
1305 -pgcon*zdo(i,ki+1)*(vs(i,ki+1)-vs(i,ki))) / &
1306 (zdo(i,ki+1)-.5*dd_massdetru(i,ki)+dd_massentru(i,ki))
1307 hcdo(i,ki)=(hcdo(i,ki+1)*zdo(i,ki+1) &
1308 -.5*dd_massdetro(i,ki)*hcdo(i,ki+1)+ &
1309 dd_massentro(i,ki)*h_entr) / &
1310 (zdo(i,ki+1)-.5*dd_massdetro(i,ki)+dd_massentro(i,ki))
1311 dbydo(i,ki)=hcdo(i,ki)-heso_cup(i,ki)
1312 bud(i)=bud(i)+dbydo(i,ki)*dzo
1313 enddo
1314
1315 if(bud(i).gt.0)then
1316 ierr(i)=7
1317#ifndef _OPENACC
1318 ierrc(i)='downdraft is not negatively buoyant '
1319#endif
1320 endif
1321 enddo
1322!$acc end parallel
1323
1324!
1326!
1327 call cup_dd_moisture(ierrc,zdo,hcdo,heso_cup,qcdo,qeso_cup, &
1328 pwdo,qo_cup,zo_cup,dd_massentro,dd_massdetro,jmin,ierr,gammao_cup, &
1329 pwevo,bu,qrcdo,po_cup,qo,heo,1, &
1330 itf,ktf, &
1331 its,ite, kts,kte)
1332!
1333!$acc kernels
1334 do i=its,itf
1335 if(ierr(i)/=0)cycle
1336 do k=kts+1,ktop(i)
1337 dp=100.*(po_cup(i,1)-po_cup(i,2))
1338 cupclw(i,k)=qrco(i,k) ! my mod
1339 cnvwt(i,k)=zuo(i,k)*cupclw(i,k)*g/dp
1340 enddo
1341 enddo
1342!$acc end kernels
1343!
1345!
1346 call cup_up_aa0(aa0,z,zu,dby,gamma_cup,t_cup, &
1347 kbcon,ktop,ierr, &
1348 itf,ktf, &
1349 its,ite, kts,kte)
1350 call cup_up_aa0(aa1,zo,zuo,dbyo,gammao_cup,tn_cup, &
1351 kbcon,ktop,ierr, &
1352 itf,ktf, &
1353 its,ite, kts,kte)
1354
1355!$acc kernels
1356 do i=its,itf
1357 if(ierr(i)/=0)cycle
1358 if(aa1(i).eq.0.)then
1359 ierr(i)=17
1360#ifndef _OPENACC
1361 ierrc(i)="cloud work function zero"
1362#endif
1363 endif
1364 enddo
1365!$acc end kernels
1366
1367!
1368!--- diurnal cycle closure
1369!
1370 !--- aa1 from boundary layer (bl) processes only
1371!$acc kernels
1372 aa1_bl(:) = 0.0
1373 xf_dicycle(:) = 0.0
1374 tau_ecmwf(:) = 0.
1375!$acc end kernels
1376 !- way to calculate the fraction of cape consumed by shallow convection
1377 !iversion=1 ! ecmwf
1378 iversion=0 ! orig
1379 !
1380 ! betchold et al 2008 time-scale of cape removal
1381!
1382! wmean is of no meaning over land....
1383! still working on replacing it over water
1384!
1385!$acc kernels
1386 do i=its,itf
1387 if(ierr(i).eq.0)then
1388 !- mean vertical velocity
1389 wmean(i) = 3.0 ! m/s ! in the future change for wmean == integral( w dz) / cloud_depth
1390 if(imid.eq.1)wmean(i) = 3.0
1391 !- time-scale cape removal from betchold et al. 2008
1392 tau_ecmwf(i)=( zo_cup(i,ktop(i))- zo_cup(i,kbcon(i)) ) / wmean(i)
1393 tau_ecmwf(i)=max(tau_ecmwf(i),720.)
1394 tau_ecmwf(i)= tau_ecmwf(i) * (1.0061 + 1.23e-2 * (dx(i)/1000.))! dx(i) must be in meters
1395 endif
1396 enddo
1397 tau_bl(:) = 0.
1398!$acc end kernels
1399
1400 !
1401 if(dicycle == 1) then
1402!$acc kernels
1403 do i=its,itf
1404
1405 if(ierr(i).eq.0)then
1406 if(xland1(i) == 0 ) then
1407 !- over water
1408 umean= 2.0+sqrt(0.5*(us(i,1)**2+vs(i,1)**2+us(i,kbcon(i))**2+vs(i,kbcon(i))**2))
1409 tau_bl(i) = (zo_cup(i,kbcon(i))- z1(i)) /umean
1410 else
1411 !- over land
1412 tau_bl(i) =( zo_cup(i,ktopdby(i))- zo_cup(i,kbcon(i)) ) / wmean(i)
1413 endif
1414
1415 endif
1416 enddo
1417!$acc end kernels
1418!$acc kernels
1419 !-get the profiles modified only by bl tendencies
1420 do i=its,itf
1421 tn_bl(i,:)=0.;qo_bl(i,:)=0.
1422 if ( ierr(i) == 0 )then
1423 !below kbcon -> modify profiles
1424 tn_bl(i,1:kbcon(i)) = tn(i,1:kbcon(i))
1425 qo_bl(i,1:kbcon(i)) = qo(i,1:kbcon(i))
1426 !above kbcon -> keep environment profiles
1427 tn_bl(i,kbcon(i)+1:ktf) = t(i,kbcon(i)+1:ktf)
1428 qo_bl(i,kbcon(i)+1:ktf) = q(i,kbcon(i)+1:ktf)
1429 endif
1430 enddo
1431!$acc end kernels
1433 call cup_env(zo,qeso_bl,heo_bl,heso_bl,tn_bl,qo_bl,po,z1, &
1434 psur,ierr,tcrit,-1, &
1435 itf,ktf, its,ite, kts,kte)
1437 call cup_env_clev(tn_bl,qeso_bl,qo_bl,heo_bl,heso_bl,zo,po,qeso_cup_bl,qo_cup_bl, &
1438 heo_cup_bl,heso_cup_bl,zo_cup,po_cup,gammao_cup_bl,tn_cup_bl,psur, &
1439 ierr,z1, &
1440 itf,ktf,its,ite, kts,kte)
1441
1442 if(iversion == 1) then
1443 !-- version ecmwf
1444 t_star=1.
1445
1446 !-- calculate pcape from bl forcing only
1448 call cup_up_aa1bl(aa1_bl,t,tn,q,qo,dtime, &
1449 zo_cup,zuo,dbyo_bl,gammao_cup_bl,tn_cup_bl, &
1450 kbcon,ktop,ierr, &
1451 itf,ktf,its,ite, kts,kte)
1452!$acc kernels
1453 do i=its,itf
1454
1455 if(ierr(i).eq.0)then
1456
1457 !- only for convection rooting in the pbl
1458 !if(zo_cup(i,kbcon(i))-z1(i) > zo(i,kpbl(i)+1)) then
1459 ! aa1_bl(i) = 0.0
1460 !else
1461 !- multiply aa1_bl the " time-scale" - tau_bl
1462 ! aa1_bl(i) = max(0.,aa1_bl(i)/t_star* tau_bl(i))
1463 aa1_bl(i) = ( aa1_bl(i)/t_star)* tau_bl(i)
1464 !endif
1465 endif
1466 enddo
1467!$acc end kernels
1468
1469 else
1470
1471 !- version for real cloud-work function
1472
1473!$acc kernels
1474 do i=its,itf
1475 if(ierr(i).eq.0)then
1476 hkbo_bl(i)=heo_cup_bl(i,k22(i))
1477 endif ! ierr
1478 enddo
1479 do k=kts,ktf
1480 do i=its,itf
1481 hco_bl(i,k)=0.
1482 dbyo_bl(i,k)=0.
1483 enddo
1484 enddo
1485 do i=its,itf
1486 if(ierr(i).eq.0)then
1487 do k=1,kbcon(i)-1
1488 hco_bl(i,k)=hkbo_bl(i)
1489 enddo
1490 k=kbcon(i)
1491 hco_bl(i,k)=hkbo_bl(i)
1492 dbyo_bl(i,k)=hkbo_bl(i) - heso_cup_bl(i,k)
1493 endif
1494 enddo
1495!
1496!
1497 do i=its,itf
1498 if(ierr(i).eq.0)then
1499 do k=kbcon(i)+1,ktop(i)
1500 hco_bl(i,k)=(hco_bl(i,k-1)*zuo(i,k-1)-.5*up_massdetro(i,k-1)*hco_bl(i,k-1)+ &
1501 up_massentro(i,k-1)*heo_bl(i,k-1)) / &
1502 (zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
1503 dbyo_bl(i,k)=hco_bl(i,k)-heso_cup_bl(i,k)
1504 enddo
1505 do k=ktop(i)+1,ktf
1506 hco_bl(i,k)=heso_cup_bl(i,k)
1507 dbyo_bl(i,k)=0.0
1508 enddo
1509 endif
1510 enddo
1511!$acc end kernels
1513 call cup_up_aa0(aa1_bl,zo,zuo,dbyo_bl,gammao_cup_bl,tn_cup_bl, &
1514 kbcon,ktop,ierr, &
1515 itf,ktf,its,ite, kts,kte)
1516!$acc kernels
1517 do i=its,itf
1518
1519 if(ierr(i).eq.0)then
1520 !- get the increment on aa0 due the bl processes
1521 aa1_bl(i) = aa1_bl(i) - aa0(i)
1522 !- only for convection rooting in the pbl
1523 !if(zo_cup(i,kbcon(i))-z1(i) > 500.0) then !- instead 500 -> zo_cup(kpbl(i))
1524 ! aa1_bl(i) = 0.0
1525 !else
1526 ! !- multiply aa1_bl the "normalized time-scale" - tau_bl/ model_timestep
1527 aa1_bl(i) = aa1_bl(i)* tau_bl(i)/ dtime
1528 !endif
1529#ifndef _OPENACC
1530! print*,'aa0,aa1bl=',aa0(i),aa1_bl(i),aa0(i)-aa1_bl(i),tau_bl(i)!,dtime,xland(i)
1531#endif
1532 endif
1533 enddo
1534!$acc end kernels
1535 endif
1536 endif ! version of implementation
1537
1538!$acc kernels
1539 axx(:)=aa1(:)
1540!$acc end kernels
1541
1542!
1544!
1545 call cup_dd_edt(ierr,us,vs,zo,ktop,kbcon,edt,po,pwavo, &
1546 pwo,ccn,ccnclean,pwevo,edtmax,edtmin,edtc,psum,psumh, &
1547 rho,aeroevap,pefc,xland1,itf,ktf, &
1548 its,ite, kts,kte)
1549 do i=its,itf
1550 if(ierr(i)/=0)cycle
1551 edto(i)=edtc(i,1)
1552 enddo
1553
1555 call get_melting_profile(ierr,tn_cup,po_cup, p_liq_ice,melting_layer,qrco &
1556 ,pwo,edto,pwdo,melting &
1557 ,itf,ktf,its,ite, kts,kte, cumulus )
1558!$acc kernels
1559 do k=kts,ktf
1560 do i=its,itf
1561 dellat_ens(i,k,1)=0.
1562 dellaq_ens(i,k,1)=0.
1563 dellaqc_ens(i,k,1)=0.
1564 pwo_ens(i,k,1)=0.
1565 enddo
1566 enddo
1567!
1568!--- change per unit mass that a model cloud would modify the environment
1569!
1570!--- 1. in bottom layer
1571!
1572 do k=kts,kte
1573 do i=its,itf
1574 dellu(i,k)=0.
1575 dellv(i,k)=0.
1576 dellah(i,k)=0.
1577 dellat(i,k)=0.
1578 dellaq(i,k)=0.
1579 dellaqc(i,k)=0.
1580 enddo
1581 enddo
1582!$acc end kernels
1583!
1584!---------------------------------------------- cloud level ktop
1585!
1586!- - - - - - - - - - - - - - - - - - - - - - - - model level ktop-1
1587! . . .
1588! . . .
1589! . . .
1590! . . .
1591! . . .
1592! . . .
1593!
1594!---------------------------------------------- cloud level k+2
1595!
1596!- - - - - - - - - - - - - - - - - - - - - - - - model level k+1
1597!
1598!---------------------------------------------- cloud level k+1
1599!
1600!- - - - - - - - - - - - - - - - - - - - - - - - model level k
1601!
1602!---------------------------------------------- cloud level k
1603!
1604! . . .
1605! . . .
1606! . . .
1607! . . .
1608! . . .
1609! . . .
1610! . . .
1611! . . .
1612! . . .
1613! . . .
1614!
1615!---------------------------------------------- cloud level 3
1616!
1617!- - - - - - - - - - - - - - - - - - - - - - - - model level 2
1618!
1619!---------------------------------------------- cloud level 2
1620!
1621!- - - - - - - - - - - - - - - - - - - - - - - - model level 1
1622!$acc kernels
1623 do i=its,itf
1624 if(ierr(i)/=0)cycle
1625 dp=100.*(po_cup(i,1)-po_cup(i,2))
1626 dellu(i,1)=pgcd*(edto(i)*zdo(i,2)*ucd(i,2) &
1627 -edto(i)*zdo(i,2)*u_cup(i,2))*g/dp &
1628 -zuo(i,2)*(uc(i,2)-u_cup(i,2)) *g/dp
1629 dellv(i,1)=pgcd*(edto(i)*zdo(i,2)*vcd(i,2) &
1630 -edto(i)*zdo(i,2)*v_cup(i,2))*g/dp &
1631 -zuo(i,2)*(vc(i,2)-v_cup(i,2)) *g/dp
1632
1633 do k=kts+1,ktop(i)
1634 ! these three are only used at or near mass detrainment and/or entrainment levels
1635 pgc=pgcon
1636 entupk=0.
1637 if(k == k22(i)-1) entupk=zuo(i,k+1)
1638 detupk=0.
1639 entdoj=0.
1640 ! detrainment and entrainment for fowndrafts
1641 detdo=edto(i)*dd_massdetro(i,k)
1642 entdo=edto(i)*dd_massentro(i,k)
1643 ! entrainment/detrainment for updraft
1644 entup=up_massentro(i,k)
1645 detup=up_massdetro(i,k)
1646 ! subsidence by downdrafts only
1647 subin=-zdo(i,k+1)*edto(i)
1648 subdown=-zdo(i,k)*edto(i)
1649 ! special levels
1650 if(k.eq.ktop(i))then
1651 detupk=zuo(i,ktop(i))
1652 subin=0.
1653 subdown=0.
1654 detdo=0.
1655 entdo=0.
1656 entup=0.
1657 detup=0.
1658 endif
1659 totmas=subin-subdown+detup-entup-entdo+ &
1660 detdo-entupk-entdoj+detupk+zuo(i,k+1)-zuo(i,k)
1661 if(abs(totmas).gt.1.e-6)then
1662#ifndef _OPENACC
1663 write(0,123)'totmas=',k22(i),kbcon(i),k,entup,detup,edto(i),zdo(i,k+1),dd_massdetro(i,k),dd_massentro(i,k)
1664123 format(a7,1x,3i3,2e12.4,1(1x,f5.2),3e12.4)
1665#endif
1666 endif
1667 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
1668 pgc=pgcon
1669 if(k.ge.ktop(i))pgc=0.
1670
1671 dellu(i,k) =-(zuo(i,k+1)*(uc(i,k+1)-u_cup(i,k+1) ) - &
1672 zuo(i,k )*(uc(i,k )-u_cup(i,k ) ) )*g/dp &
1673 +(zdo(i,k+1)*(ucd(i,k+1)-u_cup(i,k+1) ) - &
1674 zdo(i,k )*(ucd(i,k )-u_cup(i,k ) ) )*g/dp*edto(i)*pgcd
1675 dellv(i,k) =-(zuo(i,k+1)*(vc(i,k+1)-v_cup(i,k+1) ) - &
1676 zuo(i,k )*(vc(i,k )-v_cup(i,k ) ) )*g/dp &
1677 +(zdo(i,k+1)*(vcd(i,k+1)-v_cup(i,k+1) ) - &
1678 zdo(i,k )*(vcd(i,k )-v_cup(i,k ) ) )*g/dp*edto(i)*pgcd
1679
1680 enddo ! k
1681
1682 enddo
1683
1684
1685 do i=its,itf
1686 !trash = 0.0
1687 !trash2 = 0.0
1688 if(ierr(i).eq.0)then
1689
1690 dp=100.*(po_cup(i,1)-po_cup(i,2))
1691
1692 dellah(i,1)=(edto(i)*zdo(i,2)*hcdo(i,2) &
1693 -edto(i)*zdo(i,2)*heo_cup(i,2))*g/dp &
1694 -zuo(i,2)*(hco(i,2)-heo_cup(i,2))*g/dp
1695
1696 dellaq(i,1)=(edto(i)*zdo(i,2)*qcdo(i,2) &
1697 -edto(i)*zdo(i,2)*qo_cup(i,2))*g/dp &
1698 -zuo(i,2)*(qco(i,2)-qo_cup(i,2))*g/dp
1699
1700 g_rain= 0.5*(pwo(i,1)+pwo(i,2))*g/dp
1701 e_dn = -0.5*(pwdo(i,1)+pwdo(i,2))*g/dp*edto(i) ! pwdo < 0 and e_dn must > 0
1702 dellaq(i,1) = dellaq(i,1)+ e_dn-g_rain
1703
1704 !--- conservation check
1705 !- water mass balance
1706 !trash = trash + (dellaq(i,1)+dellaqc(i,1)+g_rain-e_dn)*dp/g
1707 !- h budget
1708 !trash2 = trash2+ (dellah(i,1))*dp/g
1709
1710
1711 do k=kts+1,ktop(i)
1712 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
1713 ! these three are only used at or near mass detrainment and/or entrainment levels
1714
1715 dellah(i,k) =-(zuo(i,k+1)*(hco(i,k+1)-heo_cup(i,k+1) ) - &
1716 zuo(i,k )*(hco(i,k )-heo_cup(i,k ) ) )*g/dp &
1717 +(zdo(i,k+1)*(hcdo(i,k+1)-heo_cup(i,k+1) ) - &
1718 zdo(i,k )*(hcdo(i,k )-heo_cup(i,k ) ) )*g/dp*edto(i)
1719
1720!---meltglac-------------------------------------------------
1721
1722 dellah(i,k) = dellah(i,k) + xlf*((1.-p_liq_ice(i,k))*0.5*(qrco(i,k+1)+qrco(i,k)) &
1723 - melting(i,k))*g/dp
1724
1725!---meltglac-------------------------------------------------
1726
1727 !- check h conservation
1728 ! trash2 = trash2+ (dellah(i,k))*dp/g
1729
1730
1731 !-- take out cloud liquid water for detrainment
1732 detup=up_massdetro(i,k)
1733 dz=zo_cup(i,k)-zo_cup(i,k-1)
1734 if(k.lt.ktop(i)) then
1735 dellaqc(i,k) = zuo(i,k)*c1d(i,k)*qrco(i,k)*dz/dp*g
1736 else
1737 dellaqc(i,k)= detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
1738 endif
1739 !---
1740 g_rain= 0.5*(pwo(i,k)+pwo(i,k+1))*g/dp
1741 e_dn = -0.5*(pwdo(i,k)+pwdo(i,k+1))*g/dp*edto(i) ! pwdo < 0 and e_dn must > 0
1742 !-- condensation source term = detrained + flux divergence of
1743 !-- cloud liquid water (qrco) + converted to rain
1744
1745 c_up = dellaqc(i,k)+(zuo(i,k+1)* qrco(i,k+1) - &
1746 zuo(i,k )* qrco(i,k ) )*g/dp + g_rain
1747! c_up = dellaqc(i,k)+ g_rain
1748 !-- water vapor budget
1749 !-- = flux divergence z*(q_c - q_env)_up_and_down &
1750 !-- - condensation term + evaporation
1751 dellaq(i,k) =-(zuo(i,k+1)*(qco(i,k+1)-qo_cup(i,k+1) ) - &
1752 zuo(i,k )*(qco(i,k )-qo_cup(i,k ) ) )*g/dp &
1753 +(zdo(i,k+1)*(qcdo(i,k+1)-qo_cup(i,k+1) ) - &
1754 zdo(i,k )*(qcdo(i,k )-qo_cup(i,k ) ) )*g/dp*edto(i) &
1755 - c_up + e_dn
1756 !- check water conservation liq+condensed (including rainfall)
1757 ! trash= trash+ (dellaq(i,k)+dellaqc(i,k)+ g_rain-e_dn)*dp/g
1758
1759 enddo ! k
1760 endif
1761
1762 enddo
1763!$acc end kernels
1764
1765444 format(1x,i2,1x,7e12.4) !,1x,f7.2,2x,e13.5)
1766!
1767!--- using dellas, calculate changed environmental profiles
1768!
1769 mbdt=.1
1770!$acc kernels
1771 do i=its,itf
1772 xaa0_ens(i,1)=0.
1773 enddo
1774
1775 do i=its,itf
1776 if(ierr(i).eq.0)then
1777 do k=kts,ktf
1778 xhe(i,k)=dellah(i,k)*mbdt+heo(i,k)
1779! xq(i,k)=max(1.e-16,(dellaqc(i,k)+dellaq(i,k))*mbdt+qo(i,k))
1780 xq(i,k)=max(1.e-16,dellaq(i,k)*mbdt+qo(i,k))
1781 dellat(i,k)=(1./cp)*(dellah(i,k)-xlv*dellaq(i,k))
1782! xt(i,k)= (dellat(i,k)-xlv/cp*dellaqc(i,k))*mbdt+tn(i,k)
1783 xt(i,k)= dellat(i,k)*mbdt+tn(i,k)
1784 xt(i,k)=max(190.,xt(i,k))
1785 enddo
1786
1787 ! Smooth dellas (HCB)
1788 do k=kts+1,ktf
1789 xt(i,k)=tn(i,k)+0.25*(dellat(i,k-1) + 2.*dellat(i,k) + dellat(i,k+1)) * mbdt
1790 xt(i,k)=max(190.,xt(i,k))
1791 xq(i,k)=max(1.e-16, qo(i,k)+0.25*(dellaq(i,k-1) + 2.*dellaq(i,k) + dellaq(i,k+1)) * mbdt)
1792 xhe(i,k)=heo(i,k)+0.25*(dellah(i,k-1) + 2.*dellah(i,k) + dellah(i,k+1)) * mbdt
1793 enddo
1794 endif
1795 enddo
1796 do i=its,itf
1797 if(ierr(i).eq.0)then
1798 xhe(i,ktf)=heo(i,ktf)
1799 xq(i,ktf)=qo(i,ktf)
1800 xt(i,ktf)=tn(i,ktf)
1801 endif
1802 enddo
1803!$acc end kernels
1804!
1806!
1807 call cup_env(xz,xqes,xhe,xhes,xt,xq,po,z1, &
1808 psur,ierr,tcrit,-1, &
1809 itf,ktf, &
1810 its,ite, kts,kte)
1811!
1813!
1814 call cup_env_clev(xt,xqes,xq,xhe,xhes,xz,po,xqes_cup,xq_cup, &
1815 xhe_cup,xhes_cup,xz_cup,po_cup,gamma_cup,xt_cup,psur, &
1816 ierr,z1, &
1817 itf,ktf, &
1818 its,ite, kts,kte)
1819!
1820!
1821!**************************** static control
1822!
1823!--- moist static energy inside cloud
1824!
1825!$acc kernels
1826 do k=kts,ktf
1827 do i=its,itf
1828 xhc(i,k)=0.
1829 xdby(i,k)=0.
1830 enddo
1831 enddo
1832!$acc end kernels
1833!$acc parallel loop private(x_add,k)
1834 do i=its,itf
1835 if(ierr(i).eq.0)then
1836 x_add = xlv*zqexec(i)+cp*ztexec(i)
1837 call get_cloud_bc(kte,xhe_cup(i,1:kte),xhkb(i),k22(i),x_add)
1838 do k=1,start_level(i)-1
1839 xhc(i,k)=xhe_cup(i,k)
1840 enddo
1841 k=start_level(i)
1842 xhc(i,k)=xhkb(i)
1843 endif !ierr
1844 enddo
1845!$acc end parallel
1846!
1847!
1848!$acc kernels
1849 do i=its,itf
1850 if(ierr(i).eq.0)then
1851!$acc loop seq
1852 do k=start_level(i) +1,ktop(i)
1853 xhc(i,k)=(xhc(i,k-1)*xzu(i,k-1)-.5*up_massdetro(i,k-1)*xhc(i,k-1) + &
1854 up_massentro(i,k-1)*xhe(i,k-1)) / &
1855 (xzu(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
1856
1857
1858!---meltglac-------------------------------------------------
1859 !
1860 !- include glaciation effects on xhc
1861 ! ------ ice content --------
1862 xhc(i,k)= xhc(i,k)+ xlf*(1.-p_liq_ice(i,k))*qrco(i,k)
1863!---meltglac-------------------------------------------------
1864
1865 xdby(i,k)=xhc(i,k)-xhes_cup(i,k)
1866 enddo
1867!$acc loop independent
1868 do k=ktop(i)+1,ktf
1869 xhc(i,k)=xhes_cup(i,k)
1870 xdby(i,k)=0.
1871 enddo
1872 endif
1873 enddo
1874!$acc end kernels
1875!
1877!
1878 call cup_up_aa0(xaa0,xz,xzu,xdby,gamma_cup,xt_cup, &
1879 kbcon,ktop,ierr, &
1880 itf,ktf, &
1881 its,ite, kts,kte)
1882!$acc parallel loop
1883 do i=its,itf
1884 if(ierr(i).eq.0)then
1885 xaa0_ens(i,1)=xaa0(i)
1886!$acc loop seq
1887 do k=kts,ktop(i)
1888!$acc loop independent
1889 do nens3=1,maxens3
1890 if(nens3.eq.7)then
1891!--- b=0
1892 pr_ens(i,nens3)=pr_ens(i,nens3) &
1893 +pwo(i,k)+edto(i)*pwdo(i,k)
1894!--- b=beta
1895 else if(nens3.eq.8)then
1896 pr_ens(i,nens3)=pr_ens(i,nens3)+ &
1897 pwo(i,k)+edto(i)*pwdo(i,k)
1898!--- b=beta/2
1899 else if(nens3.eq.9)then
1900 pr_ens(i,nens3)=pr_ens(i,nens3) &
1901 + pwo(i,k)+edto(i)*pwdo(i,k)
1902 else
1903 pr_ens(i,nens3)=pr_ens(i,nens3)+ &
1904 pwo(i,k) +edto(i)*pwdo(i,k)
1905 endif
1906 enddo
1907 enddo
1908 if(pr_ens(i,7).lt.1.e-6)then
1909 ierr(i)=18
1910#ifndef _OPENACC
1911 ierrc(i)="total normalized condensate too small"
1912#endif
1913 do nens3=1,maxens3
1914 pr_ens(i,nens3)=0.
1915 enddo
1916 endif
1917 do nens3=1,maxens3
1918 if(pr_ens(i,nens3).lt.1.e-5)then
1919 pr_ens(i,nens3)=0.
1920 endif
1921 enddo
1922 endif
1923 enddo
1924!$acc end parallel
1925 200 continue
1926!
1927!--- large scale forcing
1928!
1929!
1930!------- check wether aa0 should have been zero, assuming this
1931! ensemble is chosen
1932!
1933!
1934!$acc kernels
1935 do i=its,itf
1936 ierr2(i)=ierr(i)
1937 ierr3(i)=ierr(i)
1938 k22x(i)=k22(i)
1939 enddo
1940!$acc end kernels
1941 call cup_maximi(heo_cup,2,kbmax,k22x,ierr, &
1942 itf,ktf, &
1943 its,ite, kts,kte)
1944 iloop=2
1945 call cup_kbcon(ierrc,cap_max_increment,iloop,k22x,kbconx,heo_cup, &
1946 heso_cup,hkbo,ierr2,kbmax,po_cup,cap_max, &
1947 ztexec,zqexec, &
1948 0,itf,ktf, &
1949 its,ite, kts,kte, &
1950 z_cup,entr_rate,heo,imid)
1951 iloop=3
1952 call cup_kbcon(ierrc,cap_max_increment,iloop,k22x,kbconx,heo_cup, &
1953 heso_cup,hkbo,ierr3,kbmax,po_cup,cap_max, &
1954 ztexec,zqexec, &
1955 0,itf,ktf, &
1956 its,ite, kts,kte, &
1957 z_cup,entr_rate,heo,imid)
1958!
1960!
1961!$acc kernels
1962 do i = its,itf
1963 mconv(i) = 0
1964 if(ierr(i)/=0)cycle
1965!$acc loop independent
1966 do k=1,ktop(i)
1967 dq=(qo_cup(i,k+1)-qo_cup(i,k))
1968!$acc atomic update
1969 mconv(i)=mconv(i)+omeg(i,k)*dq/g
1970 enddo
1971 enddo
1972!$acc end kernels
1973 call cup_forcing_ens_3d(closure_n,xland1,aa0,aa1,xaa0_ens,mbdt,dtime, &
1974 ierr,ierr2,ierr3,xf_ens,axx,forcing, &
1975 maxens3,mconv,rand_clos, &
1976 po_cup,ktop,omeg,zdo,zdm,k22,zuo,pr_ens,edto,edtm,kbcon, &
1977 ichoice, &
1978 imid,ipr,itf,ktf, &
1979 its,ite, kts,kte, &
1980 dicycle,tau_ecmwf,aa1_bl,xf_dicycle)
1981!
1982!$acc kernels
1983 do k=kts,ktf
1984 do i=its,itf
1985 if(ierr(i).eq.0)then
1986 dellat_ens(i,k,1)=dellat(i,k)
1987 dellaq_ens(i,k,1)=dellaq(i,k)
1988 dellaqc_ens(i,k,1)=dellaqc(i,k)
1989 pwo_ens(i,k,1)=pwo(i,k) +edto(i)*pwdo(i,k)
1990 else
1991 dellat_ens(i,k,1)=0.
1992 dellaq_ens(i,k,1)=0.
1993 dellaqc_ens(i,k,1)=0.
1994 pwo_ens(i,k,1)=0.
1995 endif
1996 enddo
1997 enddo
1998!$acc end kernels
1999
2000 250 continue
2001!
2002!--- feedback
2003!
2004 if(imid.eq.1 .and. ichoice .le.2)then
2005!$acc kernels
2006 do i=its,itf
2007 !-boundary layer qe
2008 xff_mid(i,1)=0.
2009 xff_mid(i,2)=0.
2010 if(ierr(i).eq.0)then
2011 blqe=0.
2012 trash=0.
2013 if(k22(i).lt.kpbl(i)+1)then
2014 do k=1,kpbl(i)
2015 blqe=blqe+100.*dhdt(i,k)*(po_cup(i,k)-po_cup(i,k+1))/g
2016 enddo
2017 trash=max((hco(i,kbcon(i))-heo_cup(i,kbcon(i))),1.e1)
2018 xff_mid(i,1)=max(0.,blqe/trash)
2019 xff_mid(i,1)=min(0.1,xff_mid(i,1))
2020 endif
2021 xff_mid(i,2)=min(0.1,.03*zws(i))
2022 forcing(i,1)=xff_mid(i,1)
2023 forcing(i,2)=xff_mid(i,2)
2024 endif
2025 enddo
2026!$acc end kernels
2027 endif
2028 call cup_output_ens_3d(xff_mid,xf_ens,ierr,dellat_ens,dellaq_ens, &
2029 dellaqc_ens,outt,outq,outqc,dx, &
2030 zuo,pre,pwo_ens,xmb,ktop, &
2031 edto,pwdo,'deep',ierr2,ierr3, &
2032 po_cup,pr_ens,maxens3, &
2033 sig,closure_n,xland1,xmbm_in,xmbs_in, &
2034 ichoice,imid,ipr,itf,ktf, &
2035 its,ite, kts,kte, &
2036 dicycle,xf_dicycle )
2037
2039
2040 call rain_evap_below_cloudbase(itf,ktf,its,ite, &
2041 kts,kte,ierr,kbcon,xmb,psur,xland,qo_cup, &
2042 po_cup,qes_cup,pwavo,edto,pwevo,pre,outt,outq) !,outbuoy)
2043
2044!
2045!
2047!
2049 if (do_smoke_transport .and. nchem > 0) then
2050!
2051! initialize tracers if they exist
2052!
2053 chem(:,:,:) = 0.
2054!$acc kernels
2055 do nv = 1,nchem
2056 do k = 1, ktf
2057 do i = 1, itf
2058 chem(i,k,nv) = max(qamin, chem3d(i,k,nv))
2059 enddo
2060 enddo
2061 enddo
2062
2063 wetdpc_deep = 0.
2064
2065 chem_pwav(:,:) = 0.
2066 chem_psum(:,:) = 0.
2067 chem_pw(:,:,:) = 0.
2068 chem_pwd(:,:,:) = 0.
2069 pwdper(:,:) = 0.
2070 chem_down(:,:,:) = 0.
2071 chem_up(:,:,:) = 0.
2072 chem_c(:,:,:) = 0.
2073 chem_cup(:,:,:) = 0.
2074
2075 do i=its,itf
2076 if(ierr(i).eq.0)then
2077 do k=kts,jmin(i)
2078 if(pwavo(i).ne.0.) pwdper(i,k)=-edtc(i,1)*pwdo(i,k)/pwavo(i)
2079 enddo
2080 pwdper(i,:)=0.
2081 do nv=1,nchem
2082 do k=kts+1,ktf
2083 chem_cup(i,k,nv)=.5*(chem(i,k-1,nv)+chem(i,k,nv))
2084 enddo
2085 chem_cup(i,kts,nv)=chem(i,kts,nv)
2086!
2087! in updraft
2088!
2089 do k=1,k22(i)
2090 chem_up(i,k,nv)=chem_cup(i,k,nv)
2091 enddo
2092 do k=k22(i)+1,ktop(i)
2093 chem_up(i,k,nv)=(chem_up(i,k-1,nv)*zuo(i,k-1) &
2094 -.5*up_massdetr(i,k-1)*chem_up(i,k-1,nv)+ &
2095 up_massentr(i,k-1)*chem(i,k-1,nv)) / &
2096 (zuo(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
2097 chem_c(i,k,nv)=fscav(nv)*chem_up(i,k,nv)
2098 dz=zo_cup(i,k)-zo_cup(i,k-1)
2099 trash2=chem_up(i,k,nv)-chem_c(i,k,nv)
2100 trash=chem_c(i,k,nv)/(1.+c0t3d(i,k)*dz)
2101 chem_pw=c0t3d(i,k)*dz*trash*zuo(i,k)
2102 chem_up(i,k,nv)=trash2+trash
2103 chem_pwav(i,nv)=chem_pwav(i,nv)+chem_pw(i,k,nv)! *g/dp
2104 enddo
2105 do k=ktop(i)+1,ktf
2106 chem_up(i,k,nv)=chem_cup(i,k,nv)
2107 enddo
2108!
2109! in downdraft
2110!
2111 chem_down(i,jmin(i)+1,nv)=chem_cup(i,jmin(i)+1,nv)
2112 chem_psum(i,nv)=0.
2113 do ki=jmin(i),2,-1
2114 dp=100.*(po_cup(i,ki)-po_cup(i,ki+1))
2115 chem_down(i,ki,nv)=(chem_down(i,ki+1,nv)*zdo(i,ki+1) &
2116 -.5_kind_phys*dd_massdetro(i,ki)*chem_down(i,ki+1,nv)+ &
2117 dd_massentro(i,ki)*chem(i,ki,nv)) / &
2118 (zdo(i,ki+1)-.5_kind_phys*dd_massdetro(i,ki)+dd_massentro(i,ki))
2119 chem_down(i,ki,nv)=chem_down(i,ki,nv)+pwdper(i,ki)*chem_pwav(i,nv)
2120 chem_pwd(i,ki,nv)=max(0._kind_phys,pwdper(i,ki)*chem_pwav(i,nv))
2121 enddo
2122! total wet deposition
2123 do k=1,ktf-1
2124 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
2125 chem_psum(i,nv)=chem_psum(i,nv)+chem_pw(i,k,nv)*g !/dp
2126 enddo
2127 chem_psum(i,nv)=chem_psum(i,nv)*xmb(i)*dtime
2128!
2129 enddo ! nchem
2130 endif ! ierr=0
2131 enddo ! i
2132
2133 dellac(:,:,:)=0.
2134
2135 do nv=1,nchem
2136 do i=its,itf
2137 if(ierr(i).eq.0)then
2138 dp=100.*(po_cup(i,1)-po_cup(i,2))
2139 dellac(i,1,nv)=dellac(i,1,nv)+(edto(i)*zdo(i,2)*chem_down(i,2,nv))*g/dp*xmb(i)
2140 if(k22(i).eq.2)then
2141 entupk=zuo(i,2)
2142 dellac(i,1,nv)=dellac(i,1,nv)-entupk*chem_cup(i,2,nv)*g/dp*xmb(i)
2143 endif
2144 do k=kts+1,ktop(i)-1
2145 detup=0.
2146 detdo=0.
2147 entup=0.
2148 entdo=0.
2149 entdoj=0.
2150 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
2151 ! entrainment/detrainment for updraft
2152 entdo=edto(i)*dd_massentro(i,k)*chem(i,k,nv)
2153 detdo=edto(i)*dd_massdetro(i,k)*.5*(chem_down(i,k+1,nv)+chem_down(i,k,nv))
2154 entup=up_massentro(i,k)*chem(i,k,nv)
2155 detup=up_massdetro(i,k)*.5*(chem_up(i,k+1,nv)+chem_up(i,k,nv))
2156 ! special levels
2157 if(k == k22(i)-1) then
2158 entup=zuo(i,k+1)*chem_cup(i,k+1,nv)
2159 detup=0.
2160 endif
2161 if(k.eq.jmin(i))entdoj=edto(i)*zdo(i,k)*chem_cup(i,k,nv)
2162! mass budget
2163 dellac(i,k,nv) =dellac(i,k,nv) + (detup+detdo-entdo-entup-entdoj)*g/dp*xmb(i)
2164 enddo
2165 dellac(i,ktop(i),nv)=zuo(i,ktop(i))*chem_up(i,ktop(i),nv)*g/dp*xmb(i)
2166 endif ! ierr
2167 enddo ! i
2168 enddo ! nchem loop
2169
2170! fct for subsidence
2171 dellac2(:,:,:)=0.
2172 massflx(:,:)=0.
2173 do nv=1,nchem
2174!$acc loop private(trcflx_in)
2175 do i=its,itf
2176 if(ierr(i).eq.0)then
2177 trcflx_in(:)=0.
2178 dtime_max=dtime
2179
2180! initialize fct routine
2181 do k=kts,ktop(i)
2182 dp=100._kind_phys*(po_cup(i,k)-po_cup(i,k+1))
2183 dtime_max=min(dtime_max,.5_kind_phys*dp)
2184 massflx(i,k)=-xmb(i)*(zuo(i,k)-edto(i)*zdo(i,k))
2185 trcflx_in(k)=massflx(i,k)*chem_cup(i,k,nv)
2186 enddo
2187 trcflx_in(1)=0.
2188 massflx(i,1)=0.
2189 call fct1d3(ktop(i),kte,dtime_max,po_cup(i,:),chem(i,:,nv),massflx(i,:), &
2190 trcflx_in,dellac2(i,:,nv),g)
2191 do k=kts,ktop(i)
2192 trash=chem(i,k,nv)
2193 chem(i,k,nv)=chem(i,k,nv) + (dellac(i,k,nv)+dellac2(i,k,nv))*dtime
2194 if(chem(i,k,nv).lt.qamin)then
2195 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
2196 wetdpc_deep(i,nv)=wetdpc_deep(i,nv)+(qamin-chem(i,k,nv))*dp/g/dtime
2197 chem(i,k,nv)=qamin
2198 endif
2199 enddo
2200 endif
2201
2202 enddo ! i
2203 enddo ! nchem loop
2204
2206 do nv = 1, nchem
2207 do i = 1, itf
2208 do k = 1, ktf
2209 if(ierr(i).eq.0) then
2210 if (k <= ktop(i)) then
2211 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
2212 wetdpc_deep(i,nv)=wetdpc_deep(i,nv) + ((chem3d(i,k,nv)-chem(i,k,nv))*dp/(g*dtime))
2213 chem3d(i,k,nv) = chem(i,k,nv)
2214 endif
2215 endif
2216 enddo
2217 wetdpc_deep(i,nv)=max(wetdpc_deep(i,nv),qamin)
2218 enddo
2219 enddo
2220!$acc end kernels
2221
2222 endif ! nchem > 0
2223
2224 k=1
2225!$acc kernels
2226 do i=its,itf
2227 if(ierr(i).eq.0 .and.pre(i).gt.0.) then
2228 forcing(i,6)=sig(i)
2229 pre(i)=max(pre(i),0.)
2230 xmb_out(i)=xmb(i)
2231 outu(i,1)=dellu(i,1)*xmb(i)
2232 outv(i,1)=dellv(i,1)*xmb(i)
2233 do k=kts+1,ktop(i)
2234 outu(i,k)=.25*(dellu(i,k-1)+2.*dellu(i,k)+dellu(i,k+1))*xmb(i)
2235 outv(i,k)=.25*(dellv(i,k-1)+2.*dellv(i,k)+dellv(i,k+1))*xmb(i)
2236 enddo
2237 elseif(ierr(i).ne.0 .or. pre(i).eq.0.)then
2238 ktop(i)=0
2239 do k=kts,kte
2240 outt(i,k)=0.
2241 outq(i,k)=0.
2242 outqc(i,k)=0.
2243 outu(i,k)=0.
2244 outv(i,k)=0.
2245 enddo
2246 endif
2247 enddo
2248!$acc end kernels
2249! rain evaporation as in sas
2250!
2251 if(irainevap.eq.1)then
2252!$acc kernels
2253 do i = its,itf
2254 rntot(i) = 0.
2255 delqev(i) = 0.
2256 delq2(i) = 0.
2257 rn(i) = 0.
2258 rntot(i) = 0.
2259 rain=0.
2260 if(ierr(i).eq.0)then
2261!$acc loop independent
2262 do k = ktop(i), 1, -1
2263 rain = pwo(i,k) + edto(i) * pwdo(i,k)
2264!$acc atomic
2265 rntot(i) = rntot(i) + rain * xmb(i)* .001 * dtime
2266 enddo
2267 endif
2268 enddo
2269 do i = its,itf
2270 qevap(i) = 0.
2271 flg(i) = .true.
2272 if(ierr(i).eq.0)then
2273 evef = edt(i) * evfact * sig(i)**2
2274 if(xland(i).gt.0.5 .and. xland(i).lt.1.5) evef = edt(i) * evfactl * sig(i)**2
2275!$acc loop seq
2276 do k = ktop(i), 1, -1
2277 rain = pwo(i,k) + edto(i) * pwdo(i,k)
2278 rn(i) = rn(i) + rain * xmb(i) * .001 * dtime
2279 !if(po(i,k).gt.400.)then
2280 if(flg(i))then
2281 q1=qo(i,k)+(outq(i,k))*dtime
2282 t1=tn(i,k)+(outt(i,k))*dtime
2283 qcond(i) = evef * (q1 - qeso(i,k)) &
2284 & / (1. + el2orc * qeso(i,k) / t1**2)
2285 dp = -100.*(p_cup(i,k+1)-p_cup(i,k))
2286 if(rn(i).gt.0. .and. qcond(i).lt.0.) then
2287 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dtime*rn(i))))
2288 qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
2289 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
2290 endif
2291 if(rn(i).gt.0..and.qcond(i).lt.0..and. &
2292 & delq2(i).gt.rntot(i)) then
2293 qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
2294 flg(i) = .false.
2295 endif
2296 if(rn(i).gt.0..and.qevap(i).gt.0.) then
2297 outq(i,k) = outq(i,k) + qevap(i)/dtime
2298 outt(i,k) = outt(i,k) - elocp * qevap(i)/dtime
2299 rn(i) = max(0.,rn(i) - .001 * qevap(i) * dp / g)
2300 pre(i) = pre(i) - qevap(i) * dp /g/dtime
2301 pre(i)=max(pre(i),0.)
2302 delqev(i) = delqev(i) + .001*dp*qevap(i)/g
2303 endif
2304 !endif ! 400mb
2305 endif
2306 enddo
2307! pre(i)=1000.*rn(i)/dtime
2308 endif
2309 enddo
2310!$acc end kernels
2311 endif
2312
2313!$acc kernels
2314 do i=its,itf
2315 if(ierr(i).eq.0) then
2316 if(aeroevap.gt.1)then
2317 ! aerosol scavagening
2318 ccnloss(i)=ccn(i)*pefc(i)*xmb(i) ! HCB
2319 ccn(i) = ccn(i) - ccnloss(i)*scav_factor
2320 endif
2321 endif
2322 enddo
2323!$acc end kernels
2324
2325!
2327!
2328!$acc kernels
2329 do i=its,itf
2330 if(ierr(i).eq.0) then
2331 dts=0.
2332 fpi=0.
2333 do k=kts,ktop(i)
2334 dp=(po_cup(i,k)-po_cup(i,k+1))*100.
2335!total ke dissiptaion estimate
2336 dts= dts -(outu(i,k)*us(i,k)+outv(i,k)*vs(i,k))*dp/g
2337! fpi needed for calcualtion of conversion to pot. energyintegrated
2338 fpi = fpi +sqrt(outu(i,k)*outu(i,k) + outv(i,k)*outv(i,k))*dp
2339 enddo
2340 if(fpi.gt.0.)then
2341 do k=kts,ktop(i)
2342 fp= sqrt((outu(i,k)*outu(i,k)+outv(i,k)*outv(i,k)))/fpi
2343 outt(i,k)=outt(i,k)+fp*dts*g/cp
2344 enddo
2345 endif
2346 endif
2347 enddo
2348!$acc end kernels
2349
2350!
2351!---------------------------done------------------------------
2352!
2353
2354 end subroutine cu_gf_deep_run
2355
2356
2359 subroutine fct1d3 (ktop,n,dt,z,tracr,massflx,trflx_in,dellac,g)
2360!$acc routine vector
2361! --- modify a 1-D array of tracer fluxes for the purpose of maintaining
2362! --- monotonicity (including positive-definiteness) in the tracer field
2363! --- during tracer transport.
2364
2365! --- the underlying transport equation is (d tracr/dt) = - (d trflx/dz)
2366! --- where dz = |z(k+1)-z(k)| (k=1,...,n) and trflx = massflx * tracr
2367! --- physical dimensions of tracr,trflx,dz are arbitrary to some extent
2368! --- but are subject to the constraint dim[trflx] = dim[tracr*(dz/dt)].
2369
2370! --- note: tracr is carried in grid cells while z and fluxes are carried on
2371! --- interfaces. interface variables at index k are at grid location k-1/2.
2372! --- sign convention: mass fluxes are considered positive in +k direction.
2373
2374! --- massflx and trflx_in must be provided independently to allow the
2375! --- algorithm to generate an auxiliary low-order (diffusive) tracer flux
2376! --- as a stepping stone toward the final product trflx_out.
2377
2378 implicit none
2379 integer,intent(in) :: n,ktop ! number of grid cells
2380 real(kind=kind_phys) ,intent(in) :: dt,g ! transport time step
2381 real(kind=kind_phys) ,intent(in) :: z(n+0) ! location of cell interfaces
2382 real(kind=kind_phys) ,intent(in) :: tracr(n) ! the transported variable
2383 real(kind=kind_phys) ,intent(in) :: massflx(n+0) ! mass flux across interfaces
2384 real(kind=kind_phys) ,intent(in) :: trflx_in(n+0) ! original tracer flux
2385 real(kind=kind_phys) ,intent(out):: dellac(n+0) ! modified tracr flux
2386 real(kind=kind_phys) :: trflx_out(n+0) ! modified tracr flux
2387 integer k,km1,kp1
2388 logical :: NaN, error=.false., vrbos=.true.
2389 real(kind=kind_phys) dtovdz(n),trmax(n),trmin(n),flx_lo(n+0),antifx(n+0),clipped(n+0), &
2390 soln_hi(n),totlin(n),totlout(n),soln_lo(n),clipin(n),clipout(n),arg
2391 real(kind=kind_phys),parameter :: epsil=1.e-22 ! prevent division by zero
2392 real(kind=kind_phys),parameter :: damp=1. ! damper of antidff flux (1=no damping)
2393 nan(arg) = .not. (arg.ge.0. .or. arg.le.0.) ! NaN detector
2394 dtovdz(:)=0.
2395 soln_lo(:)=0.
2396 antifx(:)=0.
2397 clipin(:)=0.
2398 totlin(:)=0.
2399 totlout(:)=0.
2400 clipout(:)=0.
2401 flx_lo(:)=0.
2402 trmin(:)=0.
2403 trmax(:)=0.
2404 clipped(:)=0.
2405 trflx_out(:)=0.
2406 do k=1,ktop
2407 dtovdz(k)=.01*dt/abs(z(k+1)-z(k))*g ! time step / grid spacing
2408 if (z(k).eq.z(k+1)) error=.true.
2409 end do
2410! if (vrbos .or. error) print '(a/(8es10.3))','(fct1d) dtovdz =',dtovdz
2411
2412 do k=2,ktop
2413 if (massflx(k).ge.0.) then
2414 flx_lo(k)=massflx(k)*tracr(k-1) ! low-order flux, upstream
2415 else
2416 flx_lo(k)=massflx(k)*tracr(k) ! low-order flux, upstream
2417 end if
2418 antifx(k)=trflx_in(k)-flx_lo(k) ! antidiffusive flux
2419 end do
2420 flx_lo( 1)=trflx_in( 1)
2421 flx_lo(ktop+1)=trflx_in(ktop+1)
2422 antifx( 1)=0.
2423 antifx(ktop+1)=0.
2424! --- clip low-ord fluxes to make sure they don't violate positive-definiteness
2425 do k=1,ktop
2426 totlout(k)=max(0.,flx_lo(k+1))-min(0.,flx_lo(k )) ! total flux out
2427 clipout(k)=min(1.,tracr(k)/max(epsil,totlout(k))/ (1.0001*dtovdz(k)))
2428 end do
2429
2430 do k=2,ktop
2431 if (massflx(k).ge.0.) then
2432 flx_lo(k)=flx_lo(k)*clipout(k-1)
2433 else
2434 flx_lo(k)=flx_lo(k)*clipout(k)
2435 end if
2436 end do
2437 if (massflx( 1).lt.0.) flx_lo( 1)=flx_lo( 1)*clipout(1)
2438 if (massflx(ktop+1).gt.0.)flx_lo(ktop+1)=flx_lo(ktop+1)*clipout(ktop)
2439
2440! --- a positive-definite low-order (diffusive) solution can now be constructed
2441
2442 do k=1,ktop
2443 soln_lo(k)=tracr(k)-(flx_lo(k+1)-flx_lo(k))*dtovdz(k) ! low-ord solutn
2444 dellac(k)=-(flx_lo(k+1)-flx_lo(k))*dtovdz(k)/dt
2445 !dellac(k)=soln_lo(k)
2446 end do
2447 return
2448 do k=1,ktop
2449 km1=max(1,k-1)
2450 kp1=min(ktop,k+1)
2451 trmax(k)= max(soln_lo(km1),soln_lo(k),soln_lo(kp1), &
2452 tracr(km1),tracr(k),tracr(kp1)) ! upper bound
2453 trmin(k)=max(0.,min(soln_lo(km1),soln_lo(k),soln_lo(kp1), &
2454 tracr(km1),tracr(k),tracr(kp1))) ! lower bound
2455 end do
2456
2457 do k=1,ktop
2458 totlin(k)=max(0.,antifx(k ))-min(0.,antifx(k+1)) ! total flux in
2459 totlout(k)=max(0.,antifx(k+1))-min(0.,antifx(k )) ! total flux out
2460
2461 clipin(k)=min(damp,(trmax(k)-soln_lo(k))/max(epsil,totlin(k)) &
2462 / (1.0001*dtovdz(k)))
2463 clipout(k)=min(damp,(soln_lo(k)-trmin(k))/max(epsil,totlout(k)) &
2464 / (1.0001*dtovdz(k)))
2465#ifndef _OPENACC
2466 if (nan(clipin(k))) print *,'(fct1d) error: clipin is NaN, k=',k
2467 if (nan(clipout(k))) print *,'(fct1d) error: clipout is NaN, k=',k
2468#endif
2469
2470 if (clipin(k).lt.0.) then
2471! print 100,'(fct1d) error: clipin < 0 at k =',k, &
2472! 'clipin',clipin(k),'trmax',trmax(k),'soln_lo',soln_lo(k), &
2473! 'totlin',totlin(k),'dt/dz',dtovdz(k)
2474 error=.true.
2475 end if
2476 if (clipout(k).lt.0.) then
2477! print 100,'(fct1d) error: clipout < 0 at k =',k, &
2478! 'clipout',clipout(k),'trmin',trmin(k),'soln_lo',soln_lo(k), &
2479! 'totlout',totlout(k),'dt/dz',dtovdz(k)
2480 error=.true.
2481 end if
2482! 100 format (a,i3/(4(a10,"=",es9.2)))
2483 end do
2484
2485 do k=2,ktop
2486 if (antifx(k).gt.0.) then
2487 clipped(k)=antifx(k)*min(clipout(k-1),clipin(k))
2488 else
2489 clipped(k)=antifx(k)*min(clipout(k),clipin(k-1))
2490 end if
2491 trflx_out(k)=flx_lo(k)+clipped(k)
2492 if (nan(trflx_out(k))) then
2493#ifndef _OPENACC
2494 print *,'(fct1d) error: trflx_out is NaN, k=',k
2495#endif
2496 error=.true.
2497 end if
2498 end do
2499 trflx_out( 1)=trflx_in( 1)
2500 trflx_out(ktop+1)=trflx_in(ktop+1)
2501 do k=1,ktop
2502 soln_hi(k)=tracr(k)-(trflx_out(k+1)-trflx_out(k))*dtovdz(k)
2503 dellac(k)=-g*(trflx_out(k+1)-trflx_out(k))*dtovdz(k)/dt
2504 !dellac(k)=soln_hi(k)
2505 end do
2506
2507#ifndef _OPENACC
2508 if (vrbos .or. error) then
2509! do k=2,ktop
2510! write(32,99)k, &
2511! 'tracr(k)', tracr(k), &
2512! 'flx_in(k)', trflx_in(k), &
2513! 'flx_in(k+1)', trflx_in(k+1), &
2514! 'flx_lo(k)', flx_lo(k), &
2515! 'flx_lo(k+1)', flx_lo(k+1), &
2516! 'soln_lo(k)', soln_lo(k), &
2517! 'trmin(k)', trmin(k), &
2518! 'trmax(k)', trmax(k), &
2519! 'totlin(k)', totlin(k), &
2520! 'totlout(k)', totlout(k), &
2521! 'clipin(k-1)', clipin(k-1), &
2522! 'clipin(k)', clipin(k), &
2523! 'clipout(k-1)', clipout(k-1), &
2524! 'clipout(k)', clipout(k), &
2525! 'antifx(k)', antifx(k), &
2526! 'antifx(k+1)', antifx(k+1), &
2527! 'clipped(k)', clipped(k), &
2528! 'clipped(k+1)', clipped(k+1), &
2529! 'flx_out(k)', trflx_out(k), &
2530! 'flx_out(k+1)', trflx_out(k+1), &
2531! 'dt/dz(k)', dtovdz(k), &
2532! 'final', tracr(k)-(trflx_out(k+1)-trflx_out(k))*dtovdz(k)
2533! 99 format ('(trc1d) k =',i4/(3(a13,'=',es13.6)))
2534! end do
2535 if (error) stop '(fct1d error)'
2536 end if
2537#endif
2538
2539 return
2540 end subroutine fct1d3
2541
2543 subroutine rain_evap_below_cloudbase(itf,ktf, its,ite, kts,kte,ierr, &
2544 kbcon,xmb,psur,xland,qo_cup, &
2545 po_cup,qes_cup,pwavo,edto,pwevo,pre,outt,outq) !,outbuoy)
2546
2547 implicit none
2548 real(kind=kind_phys), parameter :: alp1=5.44e-4 & !1/sec
2549 ,alp2=5.09e-3 & !unitless
2550 ,alp3=0.5777 & !unitless
2551 ,c_conv=0.05 !conv fraction area, unitless
2552
2553
2554 integer ,intent(in) :: itf,ktf, its,ite, kts,kte
2555 integer, dimension(its:ite) ,intent(in) :: ierr,kbcon
2556 real(kind=kind_phys), dimension(its:ite) ,intent(in) ::psur,xland,pwavo,edto,pwevo,xmb
2557 real(kind=kind_phys), dimension(its:ite,kts:kte),intent(in) :: po_cup,qo_cup,qes_cup
2558 real(kind=kind_phys), dimension(its:ite) ,intent(inout) :: pre
2559 real(kind=kind_phys), dimension(its:ite,kts:kte),intent(inout) :: outt,outq !,outbuoy
2560!$acc declare copyin(ierr,kbcon,psur,xland,pwavo,edto,pwevo,xmb,po_cup,qo_cup,qes_cup)
2561!$acc declare copy(pre,outt,outq)
2562
2563 !real, dimension(its:ite) ,intent(out) :: tot_evap_bcb
2564 !real, dimension(its:ite,kts:kte),intent(out) :: evap_bcb,net_prec_bcb
2565
2566 !-- locals
2567 integer :: i,k
2568 real(kind=kind_phys) :: rh_cr , del_t,del_q,dp,q_deficit
2569 real(kind=kind_phys), dimension(its:ite,kts:kte) :: evap_bcb,net_prec_bcb
2570 real(kind=kind_phys), dimension(its:ite) :: tot_evap_bcb
2571!$acc declare create(evap_bcb,net_prec_bcb,tot_evap_bcb)
2572
2573!$acc kernels
2574 do i=its,itf
2575 evap_bcb(i,:)= 0.0
2576 net_prec_bcb(i,:)= 0.0
2577 tot_evap_bcb(i) = 0.0
2578 if(ierr(i) /= 0) cycle
2579
2580 !-- critical rel humidity
2581 rh_cr=0.9*xland(i)+0.7*(1-xland(i))
2582 !RH_cr=1.
2583
2584 !-- net precipitation (after downdraft evap) at cloud base, available to
2585 !evap
2586 k=kbcon(i)
2587 !net_prec_bcb(i,k) = xmb(i)*(pwavo(i)+edto(i)*pwevo(i)) !-- pwevo<0.
2588 net_prec_bcb(i,k) = pre(i)
2589
2590!$acc loop seq
2591 do k=kbcon(i)-1, kts, -1
2592
2593 q_deficit = max(0.,(rh_cr*qes_cup(i,k) -qo_cup(i,k)))
2594
2595 if(q_deficit < 1.e-6) then
2596 net_prec_bcb(i,k)= net_prec_bcb(i,k+1)
2597 cycle
2598 endif
2599
2600 dp = 100.*(po_cup(i,k)-po_cup(i,k+1))
2601
2602 !--units here: kg[water]/kg[air}/sec
2603 evap_bcb(i,k) = c_conv * alp1 * q_deficit * &
2604 ( sqrt(po_cup(i,k)/psur(i))/alp2 *net_prec_bcb(i,k+1)/c_conv )**alp3
2605
2606 !--units here: kg[water]/kg[air}/sec * kg[air]/m3 * m = kg[water]/m2/sec
2607 evap_bcb(i,k)= evap_bcb(i,k)*dp/g
2608
2609 if((net_prec_bcb(i,k+1) - evap_bcb(i,k)).lt.0.) cycle
2610 if((pre(i) - evap_bcb(i,k)).lt.0.) cycle
2611 net_prec_bcb(i,k)= net_prec_bcb(i,k+1) - evap_bcb(i,k)
2612
2613 tot_evap_bcb(i) = tot_evap_bcb(i)+evap_bcb(i,k)
2614
2615 !-- feedback
2616 del_q = evap_bcb(i,k)*g/dp ! > 0., units: kg[water]/kg[air}/sec
2617 del_t = -evap_bcb(i,k)*g/dp*(xlv/cp) ! < 0., units: K/sec
2618
2619! print*,"ebcb2",k,del_q*86400,del_t*86400
2620
2621 outq(i,k) = outq(i,k) + del_q
2622 outt(i,k) = outt(i,k) + del_t
2623 !outbuoy(i,k) = outbuoy(i,k) + cp*del_t+xlv*del_q
2624
2625 pre(i) = pre(i) - evap_bcb(i,k)
2626 enddo
2627 enddo
2628!$acc end kernels
2629
2630 end subroutine rain_evap_below_cloudbase
2631
2634 subroutine cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, &
2635 pw,ccn,ccnclean,pwev,edtmax,edtmin,edtc,psum2,psumh, &
2636 rho,aeroevap,pefc,xland1,itf,ktf, &
2637 its,ite, kts,kte )
2638
2639 implicit none
2640
2641 integer &
2642 ,intent (in ) :: &
2643 aeroevap,itf,ktf, &
2644 its,ite, kts,kte
2645 !
2646 ! ierr error value, maybe modified in this routine
2647 !
2648 real(kind=kind_phys), dimension (its:ite,kts:kte) &
2649 ,intent (in ) :: &
2650 rho,us,vs,z,p,pw
2651 real(kind=kind_phys), dimension (its:ite,1) &
2652 ,intent (out ) :: &
2653 edtc
2654 real(kind=kind_phys), dimension (its:ite) &
2655 ,intent (out ) :: &
2656 pefc
2657 real(kind=kind_phys), dimension (its:ite) &
2658 ,intent (out ) :: &
2659 edt
2660 real(kind=kind_phys), dimension (its:ite) &
2661 ,intent (in ) :: &
2662 pwav,pwev,psum2,psumh,edtmax,edtmin
2663 integer, dimension (its:ite) &
2664 ,intent (in ) :: &
2665 ktop,kbcon,xland1
2666 real(kind=kind_phys), intent (in ) :: & !HCB
2667 ccnclean
2668 real(kind=kind_phys), dimension (its:ite) &
2669 ,intent (inout ) :: &
2670 ccn
2671 integer, dimension (its:ite) &
2672 ,intent (inout) :: &
2673 ierr
2674!$acc declare copyin(rho,us,vs,z,p,pw,pwav,pwev,psum2,psumh,edtmax,edtmin,ktop,kbcon)
2675!$acc declare copyout(edtc,edt) copy(ccn,ierr)
2676!
2677! local variables in this routine
2678!
2679
2680 integer i,k,kk
2681 real(kind=kind_phys) einc,pef,pefb,prezk,zkbc
2682 real(kind=kind_phys), dimension (its:ite) :: &
2683 vshear,sdp,vws
2684!$acc declare create(vshear,sdp,vws)
2685 real(kind=kind_phys) :: prop_c,aeroadd,alpha3,beta3
2686 prop_c=0. !10.386
2687 alpha3 = 0.75
2688 beta3 = -0.15
2689 pefc(:)=0.
2690 pefb=0.
2691 pef=0.
2692
2693!
2694!--- determine downdraft strength in terms of windshear
2695!
2696! */ calculate an average wind shear over the depth of the cloud
2697!
2698!$acc kernels
2699 do i=its,itf
2700 edt(i)=0.
2701 vws(i)=0.
2702 sdp(i)=0.
2703 vshear(i)=0.
2704 enddo
2705 do i=its,itf
2706 edtc(i,1)=0.
2707 enddo
2708 do kk = kts,ktf-1
2709 do 62 i=its,itf
2710 if(ierr(i).ne.0)go to 62
2711 if (kk .le. min0(ktop(i),ktf) .and. kk .ge. kbcon(i)) then
2712 vws(i) = vws(i)+ &
2713 (abs((us(i,kk+1)-us(i,kk))/(z(i,kk+1)-z(i,kk))) &
2714 + abs((vs(i,kk+1)-vs(i,kk))/(z(i,kk+1)-z(i,kk)))) * &
2715 (p(i,kk) - p(i,kk+1))
2716 sdp(i) = sdp(i) + p(i,kk) - p(i,kk+1)
2717 endif
2718 if (kk .eq. ktf-1)vshear(i) = 1.e3 * vws(i) / sdp(i)
2719 62 continue
2720 end do
2721 do i=its,itf
2722 if(ierr(i).eq.0)then
2723 pef=(1.591-.639*vshear(i)+.0953*(vshear(i)**2) &
2724 -.00496*(vshear(i)**3))
2725 if(pef.gt.0.9)pef=0.9
2726 if(pef.lt.0.1)pef=0.1
2727!
2728!--- cloud base precip efficiency
2729!
2730 zkbc=z(i,kbcon(i))*3.281e-3
2731 prezk=.02
2732 if(zkbc.gt.3.)then
2733 prezk=.96729352+zkbc*(-.70034167+zkbc*(.162179896+zkbc &
2734 *(- 1.2569798e-2+zkbc*(4.2772e-4-zkbc*5.44e-6))))
2735 endif
2736 if(zkbc.gt.25)then
2737 prezk=2.4
2738 endif
2739 pefb=1./(1.+prezk)
2740 if(pefb.gt.0.9)pefb=0.9
2741 if(pefb.lt.0.1)pefb=0.1
2742 pefb=pef
2743
2744 edt(i)=1.-.5*(pefb+pef)
2745 if(aeroevap.gt.1)then
2746 pefb=.5
2747 if(xland1(i) == 1)pefb=.3
2748 aeroadd=0.
2749 if((psumh(i)>0.).and.(psum2(i)>0.))then
2750 aeroadd=((ccnclean)**beta3)*(psumh(i)**(alpha3-1))
2751 prop_c=pefb/aeroadd
2752 aeroadd=((ccn(i))**beta3)*(psum2(i)**(alpha3-1))
2753 aeroadd=prop_c*aeroadd
2754 pefc(i)=aeroadd
2755
2756 if(pefc(i).gt.0.9)pefc(i)=0.9
2757 if(pefc(i).lt.0.1)pefc(i)=0.1
2758 edt(i)=1.-pefc(i)
2759 endif
2760 endif
2761
2762
2763!--- edt here is 1-precipeff!
2764 edtc(i,1)=edt(i)
2765 endif
2766 enddo
2767 do i=its,itf
2768 if(ierr(i).eq.0)then
2769 edtc(i,1)=-edtc(i,1)*psum2(i)/pwev(i)
2770 if(edtc(i,1).gt.edtmax(i))edtc(i,1)=edtmax(i)
2771 if(edtc(i,1).lt.edtmin(i))edtc(i,1)=edtmin(i)
2772 endif
2773 enddo
2774!$acc end kernels
2775
2776 end subroutine cup_dd_edt
2777
2779 subroutine cup_dd_moisture(ierrc,zd,hcd,hes_cup,qcd,qes_cup, &
2780 pwd,q_cup,z_cup,dd_massentr,dd_massdetr,jmin,ierr, &
2781 gamma_cup,pwev,bu,qrcd,p_cup, &
2782 q,he,iloop, &
2783 itf,ktf, &
2784 its,ite, kts,kte )
2785
2786 implicit none
2787
2788 integer &
2789 ,intent (in ) :: &
2790 itf,ktf, &
2791 its,ite, kts,kte
2792 ! cdd= detrainment function
2793 ! q = environmental q on model levels
2794 ! q_cup = environmental q on model cloud levels
2795 ! qes_cup = saturation q on model cloud levels
2796 ! hes_cup = saturation h on model cloud levels
2797 ! hcd = h in model cloud
2798 ! bu = buoancy term
2799 ! zd = normalized downdraft mass flux
2800 ! gamma_cup = gamma on model cloud levels
2801 ! mentr_rate = entrainment rate
2802 ! qcd = cloud q (including liquid water) after entrainment
2803 ! qrch = saturation q in cloud
2804 ! pwd = evaporate at that level
2805 ! pwev = total normalized integrated evaoprate (i2)
2806 ! entr= entrainment rate
2807 !
2808 real(kind=kind_phys), dimension (its:ite,kts:kte) &
2809 ,intent (in ) :: &
2810 zd,hes_cup,hcd,qes_cup,q_cup,z_cup, &
2811 dd_massentr,dd_massdetr,gamma_cup,q,he,p_cup
2812!$acc declare copyin(zd,hes_cup,hcd,qes_cup,q_cup,z_cup,dd_massentr,dd_massdetr,gamma_cup,q,he)
2813 integer &
2814 ,intent (in ) :: &
2815 iloop
2816 integer, dimension (its:ite) &
2817 ,intent (in ) :: &
2818 jmin
2819!$acc declare copyin(jmin)
2820 integer, dimension (its:ite) &
2821 ,intent (inout) :: &
2822 ierr
2823!$acc declare copy(ierr)
2824 real(kind=kind_phys), dimension (its:ite,kts:kte)&
2825 ,intent (out ) :: &
2826 qcd,qrcd,pwd
2827 real(kind=kind_phys), dimension (its:ite)&
2828 ,intent (out ) :: &
2829 pwev,bu
2830!$acc declare copyout(qcd,qrcd,pwd,pwev,bu)
2831 character*50 :: ierrc(its:ite)
2832!
2833! local variables in this routine
2834!
2835
2836 integer :: &
2837 i,k,ki
2838 real(kind=kind_phys) :: &
2839 denom,dp,dh,dz,dqeva
2840
2841!$acc kernels
2842 do i=its,itf
2843 bu(i)=0.
2844 pwev(i)=0.
2845 enddo
2846 do k=kts,ktf
2847 do i=its,itf
2848 qcd(i,k)=0.
2849 qrcd(i,k)=0.
2850 pwd(i,k)=0.
2851 enddo
2852 enddo
2853!
2854!
2855!
2856 do 100 i=its,itf
2857 if(ierr(i).eq.0)then
2858 k=jmin(i)
2859 dz=z_cup(i,k+1)-z_cup(i,k)
2860 dp=-100.*(p_cup(i,k+1)-p_cup(i,k))
2861 qcd(i,k)=q_cup(i,k)
2862 dh=hcd(i,k)-hes_cup(i,k)
2863 if(dh.lt.0)then
2864 qrcd(i,k)=(qes_cup(i,k)+(1./xlv)*(gamma_cup(i,k) &
2865 /(1.+gamma_cup(i,k)))*dh)
2866 else
2867 qrcd(i,k)=qes_cup(i,k)
2868 endif
2869 pwd(i,jmin(i))=zd(i,jmin(i))*min(0.,qcd(i,k)-qrcd(i,k))
2870 qcd(i,k)=qrcd(i,k)
2871 pwev(i)=pwev(i)+pwd(i,jmin(i))*g/dp ! *dz
2872!
2873 bu(i)=dz*dh
2874!$acc loop seq
2875 do ki=jmin(i)-1,1,-1
2876 dz=z_cup(i,ki+1)-z_cup(i,ki)
2877 dp=-100.*(p_cup(i,ki+1)-p_cup(i,ki))
2878! qcd(i,ki)=(qcd(i,ki+1)*(1.-.5*cdd(i,ki+1)*dz) &
2879! +entr*dz*q(i,ki) &
2880! )/(1.+entr*dz-.5*cdd(i,ki+1)*dz)
2881! dz=qcd(i,ki)
2882!print*,"i=",i," k=",ki," qcd(i,ki+1)=",qcd(i,ki+1)
2883!print*,"zd=",zd(i,ki+1)," dd_ma=",dd_massdetr(i,ki)," q=",q(i,ki)
2884!joe-added check for non-zero denominator:
2885 denom=zd(i,ki+1)-.5*dd_massdetr(i,ki)+dd_massentr(i,ki)
2886 if(denom.lt.1.e-16)then
2887 ierr(i)=51
2888 exit
2889 endif
2890 qcd(i,ki)=(qcd(i,ki+1)*zd(i,ki+1) &
2891 -.5*dd_massdetr(i,ki)*qcd(i,ki+1)+ &
2892 dd_massentr(i,ki)*q(i,ki)) / &
2893 (zd(i,ki+1)-.5*dd_massdetr(i,ki)+dd_massentr(i,ki))
2894!
2895!--- to be negatively buoyant, hcd should be smaller than hes!
2896!--- ideally, dh should be negative till dd hits ground, but that is not always
2897!--- the case
2898!
2899 dh=hcd(i,ki)-hes_cup(i,ki)
2900 bu(i)=bu(i)+dz*dh
2901 qrcd(i,ki)=qes_cup(i,ki)+(1./xlv)*(gamma_cup(i,ki) &
2902 /(1.+gamma_cup(i,ki)))*dh
2903 dqeva=qcd(i,ki)-qrcd(i,ki)
2904 if(dqeva.gt.0.)then
2905 dqeva=0.
2906 qrcd(i,ki)=qcd(i,ki)
2907 endif
2908 pwd(i,ki)=zd(i,ki)*dqeva
2909 qcd(i,ki)=qrcd(i,ki)
2910 pwev(i)=pwev(i)+pwd(i,ki)*g/dp
2911 enddo
2912!
2913!--- end loop over i
2914 if( (pwev(i).eq.0.) .and. (iloop.eq.1))then
2915! print *,'problem with buoy in cup_dd_moisture',i
2916 ierr(i)=7
2917#ifndef _OPENACC
2918 ierrc(i)="problem with buoy in cup_dd_moisture"
2919#endif
2920 endif
2921 if(bu(i).ge.0.and.iloop.eq.1)then
2922! print *,'problem with buoy in cup_dd_moisture',i
2923 ierr(i)=7
2924#ifndef _OPENACC
2925 ierrc(i)="problem2 with buoy in cup_dd_moisture"
2926#endif
2927 endif
2928 endif
2929100 continue
2930!$acc end kernels
2931
2932 end subroutine cup_dd_moisture
2933
2936 subroutine cup_env(z,qes,he,hes,t,q,p,z1, &
2937 psur,ierr,tcrit,itest, &
2938 itf,ktf, &
2939 its,ite, kts,kte )
2940
2941 implicit none
2942
2943 integer &
2944 ,intent (in ) :: &
2945 itf,ktf, &
2946 its,ite, kts,kte
2947 !
2948 !
2949 real(kind=kind_phys), dimension (its:ite,kts:kte) &
2950 ,intent (in ) :: &
2951 p,t,q
2952!$acc declare copyin(p,t,q)
2953 real(kind=kind_phys), dimension (its:ite,kts:kte) &
2954 ,intent (out ) :: &
2955 hes,qes
2956!$acc declare copyout(hes,qes)
2957 real(kind=kind_phys), dimension (its:ite,kts:kte) &
2958 ,intent (inout) :: &
2959 he,z
2960!$acc declare copy(he,z)
2961 real(kind=kind_phys), dimension (its:ite) &
2962 ,intent (in ) :: &
2963 psur,z1
2964!$acc declare copyin(psur,z1)
2965 integer, dimension (its:ite) &
2966 ,intent (inout) :: &
2967 ierr
2968!$acc declare copy(ierr)
2969 integer &
2970 ,intent (in ) :: &
2971 itest
2972!
2973! local variables in this routine
2974!
2975
2976 integer :: &
2977 i,k
2978! real(kind=kind_phys), dimension (1:2) :: ae,be,ht
2979 real(kind=kind_phys), dimension (its:ite,kts:kte) :: tv
2980!$acc declare create(tv)
2981 real(kind=kind_phys) :: tcrit,e,tvbar
2982! real(kind=kind_phys), external :: satvap
2983! real(kind=kind_phys) :: satvap
2984
2985
2986! ht(1)=xlv/cp
2987! ht(2)=2.834e6/cp
2988! be(1)=.622*ht(1)/.286
2989! ae(1)=be(1)/273.+alog(610.71)
2990! be(2)=.622*ht(2)/.286
2991! ae(2)=be(2)/273.+alog(610.71)
2992!$acc parallel loop collapse(2) private(e)
2993 do k=kts,ktf
2994 do i=its,itf
2995 if(ierr(i).eq.0)then
2996!csgb - iph is for phase, dependent on tcrit (water or ice)
2997! iph=1
2998! if(t(i,k).le.tcrit)iph=2
2999! print *, 'ae(iph),be(iph) = ',ae(iph),be(iph),ae(iph)-be(iph),t(i,k),i,k
3000! e=exp(ae(iph)-be(iph)/t(i,k))
3001! print *, 'p, e = ', p(i,k), e
3002! qes(i,k)=.622*e/(100.*p(i,k)-e)
3003 e=satvap(t(i,k))
3004 qes(i,k)=0.622*e/max(1.e-8,(p(i,k)-e))
3005 if(qes(i,k).le.1.e-16)qes(i,k)=1.e-16
3006 if(qes(i,k).lt.q(i,k))qes(i,k)=q(i,k)
3007! if(q(i,k).gt.qes(i,k))q(i,k)=qes(i,k)
3008 tv(i,k)=t(i,k)+.608*q(i,k)*t(i,k)
3009 endif
3010 enddo
3011 enddo
3012!$acc end parallel
3013!
3014!--- z's are calculated with changed h's and q's and t's
3015!--- if itest=2
3016!
3017 if(itest.eq.1 .or. itest.eq.0)then
3018!$acc kernels
3019 do i=its,itf
3020 if(ierr(i).eq.0)then
3021 z(i,1)=max(0.,z1(i))-(log(p(i,1))- &
3022 log(psur(i)))*287.*tv(i,1)/9.81
3023 endif
3024 enddo
3025
3026! --- calculate heights
3027!$acc loop seq
3028 do k=kts+1,ktf
3029!$acc loop private(tvbar)
3030 do i=its,itf
3031 if(ierr(i).eq.0)then
3032 tvbar=.5*tv(i,k)+.5*tv(i,k-1)
3033 z(i,k)=z(i,k-1)-(log(p(i,k))- &
3034 log(p(i,k-1)))*287.*tvbar/9.81
3035 endif
3036 enddo
3037 enddo
3038!$acc end kernels
3039 else if(itest.eq.2)then
3040!$acc kernels
3041 do k=kts,ktf
3042 do i=its,itf
3043 if(ierr(i).eq.0)then
3044 z(i,k)=(he(i,k)-1004.*t(i,k)-2.5e6*q(i,k))/9.81
3045 z(i,k)=max(1.e-3,z(i,k))
3046 endif
3047 enddo
3048 enddo
3049!$acc end kernels
3050 else if(itest.eq.-1)then
3051 endif
3052!
3053!--- calculate moist static energy - he
3054! saturated moist static energy - hes
3055!
3056!$acc kernels
3057 do k=kts,ktf
3058 do i=its,itf
3059 if(ierr(i).eq.0)then
3060 if(itest.le.0)he(i,k)=9.81*z(i,k)+1004.*t(i,k)+2.5e06*q(i,k)
3061 hes(i,k)=9.81*z(i,k)+1004.*t(i,k)+2.5e06*qes(i,k)
3062 if(he(i,k).ge.hes(i,k))he(i,k)=hes(i,k)
3063 endif
3064 enddo
3065 enddo
3066!$acc end kernels
3067
3068 end subroutine cup_env
3069
3090 subroutine cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup, &
3091 he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
3092 ierr,z1, &
3093 itf,ktf, &
3094 its,ite, kts,kte )
3095
3096 implicit none
3097
3098 integer &
3099 ,intent (in ) :: &
3100 itf,ktf, &
3101 its,ite, kts,kte
3102 !
3103 real(kind=kind_phys), dimension (its:ite,kts:kte) &
3104 ,intent (in ) :: &
3105 qes,q,he,hes,z,p,t
3106!$acc declare copyin(qes,q,he,hes,z,p,t)
3107 real(kind=kind_phys), dimension (its:ite,kts:kte) &
3108 ,intent (out ) :: &
3109 qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup
3110!$acc declare copyout(qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup)
3111 real(kind=kind_phys), dimension (its:ite) &
3112 ,intent (in ) :: &
3113 psur,z1
3114!$acc declare copyin(psur,z1)
3115 integer, dimension (its:ite) &
3116 ,intent (inout) :: &
3117 ierr
3118!$acc declare copy(ierr)
3119!
3120! local variables in this routine
3121!
3122
3123 integer :: &
3124 i,k
3125
3126!$acc kernels
3127 do k=kts,ktf
3128 do i=its,itf
3129 qes_cup(i,k)=0.
3130 q_cup(i,k)=0.
3131 hes_cup(i,k)=0.
3132 he_cup(i,k)=0.
3133 z_cup(i,k)=0.
3134 p_cup(i,k)=0.
3135 t_cup(i,k)=0.
3136 gamma_cup(i,k)=0.
3137 enddo
3138 enddo
3139 do k=kts+1,ktf
3140 do i=its,itf
3141 if(ierr(i).eq.0)then
3142 qes_cup(i,k)=.5*(qes(i,k-1)+qes(i,k))
3143 q_cup(i,k)=.5*(q(i,k-1)+q(i,k))
3144 hes_cup(i,k)=.5*(hes(i,k-1)+hes(i,k))
3145 he_cup(i,k)=.5*(he(i,k-1)+he(i,k))
3146 if(he_cup(i,k).gt.hes_cup(i,k))he_cup(i,k)=hes_cup(i,k)
3147 z_cup(i,k)=.5*(z(i,k-1)+z(i,k))
3148 p_cup(i,k)=.5*(p(i,k-1)+p(i,k))
3149 t_cup(i,k)=.5*(t(i,k-1)+t(i,k))
3150 gamma_cup(i,k)=(xlv/cp)*(xlv/(r_v*t_cup(i,k) &
3151 *t_cup(i,k)))*qes_cup(i,k)
3152 endif
3153 enddo
3154 enddo
3155 do i=its,itf
3156 if(ierr(i).eq.0)then
3157 qes_cup(i,1)=qes(i,1)
3158 q_cup(i,1)=q(i,1)
3159! hes_cup(i,1)=hes(i,1)
3160! he_cup(i,1)=he(i,1)
3161 hes_cup(i,1)=9.81*z1(i)+1004.*t(i,1)+2.5e6*qes(i,1)
3162 he_cup(i,1)=9.81*z1(i)+1004.*t(i,1)+2.5e6*q(i,1)
3163 z_cup(i,1)=.5*(z(i,1)+z1(i))
3164 p_cup(i,1)=.5*(p(i,1)+psur(i))
3165 z_cup(i,1)=z1(i)
3166 p_cup(i,1)=psur(i)
3167 t_cup(i,1)=t(i,1)
3168 gamma_cup(i,1)=xlv/cp*(xlv/(r_v*t_cup(i,1) &
3169 *t_cup(i,1)))*qes_cup(i,1)
3170 endif
3171 enddo
3172!$acc end kernels
3173 end subroutine cup_env_clev
3174
3177 subroutine cup_forcing_ens_3d(closure_n,xland,aa0,aa1,xaa0,mbdt,dtime,ierr,ierr2,ierr3,&
3178 xf_ens,axx,forcing,maxens3,mconv,rand_clos, &
3179 p_cup,ktop,omeg,zd,zdm,k22,zu,pr_ens,edt,edtm,kbcon, &
3180 ichoice, &
3181 imid,ipr,itf,ktf, &
3182 its,ite, kts,kte, &
3183 dicycle,tau_ecmwf,aa1_bl,xf_dicycle )
3184
3185 implicit none
3186
3187 integer &
3188 ,intent (in ) :: &
3189 imid,ipr,itf,ktf, &
3190 its,ite, kts,kte
3191 integer, intent (in ) :: &
3192 maxens3
3193 !
3194 ! ierr error value, maybe modified in this routine
3195 ! pr_ens = precipitation ensemble
3196 ! xf_ens = mass flux ensembles
3197 ! massfln = downdraft mass flux ensembles used in next timestep
3198 ! omeg = omega from large scale model
3199 ! mconv = moisture convergence from large scale model
3200 ! zd = downdraft normalized mass flux
3201 ! zu = updraft normalized mass flux
3202 ! aa0 = cloud work function without forcing effects
3203 ! aa1 = cloud work function with forcing effects
3204 ! xaa0 = cloud work function with cloud effects (ensemble dependent)
3205 ! edt = epsilon
3206 ! dir = "storm motion"
3207 ! mbdt = arbitrary numerical parameter
3208 ! dtime = dt over which forcing is applied
3209 ! iact_gr_old = flag to tell where convection was active
3210 ! kbcon = lfc of parcel from k22
3211 ! k22 = updraft originating level
3212 ! ichoice = flag if only want one closure (usually set to zero!)
3213 !
3214 real(kind=kind_phys), dimension (its:ite,1:maxens3) &
3215 ,intent (inout) :: &
3216 pr_ens
3217 real(kind=kind_phys), dimension (its:ite,1:maxens3) &
3218 ,intent (inout ) :: &
3219 xf_ens
3220!$acc declare copy(pr_ens,xf_ens)
3221 real(kind=kind_phys), dimension (its:ite,kts:kte) &
3222 ,intent (in ) :: &
3223 zd,zu,p_cup,zdm
3224 real(kind=kind_phys), dimension (its:ite,kts:kte) &
3225 ,intent (in ) :: &
3226 omeg
3227 real(kind=kind_phys), dimension (its:ite,1) &
3228 ,intent (in ) :: &
3229 xaa0
3230 real(kind=kind_phys), dimension (its:ite,4) &
3231 ,intent (in ) :: &
3232 rand_clos
3233 real(kind=kind_phys), dimension (its:ite) &
3234 ,intent (in ) :: &
3235 aa1,edt,edtm
3236 real(kind=kind_phys), dimension (its:ite) &
3237 ,intent (in ) :: &
3238 mconv,axx
3239!$acc declare copyin(zd,zu,p_cup,zdm,omeg,xaa0,rand_clos,aa1,edt,edtm,mconv,axx)
3240 real(kind=kind_phys), dimension (its:ite) &
3241 ,intent (inout) :: &
3242 aa0,closure_n
3243!$acc declare copy(aa0,closure_n)
3244 real(kind=kind_phys) &
3245 ,intent (in ) :: &
3246 mbdt
3247 real(kind=kind_phys) &
3248 ,intent (in ) :: &
3249 dtime
3250 integer, dimension (its:ite) &
3251 ,intent (inout ) :: &
3252 k22,kbcon,ktop
3253 integer, dimension (its:ite) &
3254 ,intent (in ) :: &
3255 xland
3256 integer, dimension (its:ite) &
3257 ,intent (inout) :: &
3258 ierr,ierr2,ierr3
3259!$acc declare copy(k22,kbcon,ktop,ierr,ierr2,ierr3) copyin(xland)
3260 integer &
3261 ,intent (in ) :: &
3262 ichoice
3263 integer, intent(in) :: dicycle
3264 real(kind=kind_phys), intent(in) , dimension (its:ite) :: aa1_bl,tau_ecmwf
3265 real(kind=kind_phys), intent(inout), dimension (its:ite) :: xf_dicycle
3266 real(kind=kind_phys), intent(inout), dimension (its:ite,10) :: forcing
3267!$acc declare copyin(aa1_bl,tau_ecmwf) copy(xf_dicycle,forcing)
3268 !- local var
3269 real(kind=kind_phys) :: xff_dicycle
3270!
3271! local variables in this routine
3272!
3273
3274 real(kind=kind_phys), dimension (1:maxens3) :: &
3275 xff_ens3
3276 real(kind=kind_phys), dimension (1) :: &
3277 xk
3278 integer :: &
3279 kk,i,k,n,ne
3280! integer, parameter :: mkxcrt=15
3281! real(kind=kind_phys), dimension(1:mkxcrt) :: &
3282! pcrit,acrit,acritt
3283 integer, dimension (its:ite) :: kloc
3284 real(kind=kind_phys) :: &
3285 a1,a_ave,xff0,xomg!,aclim1,aclim2,aclim3,aclim4
3286
3287 real(kind=kind_phys), dimension (its:ite) :: ens_adj
3288!$acc declare create(kloc,ens_adj)
3289
3290
3291
3292!
3293!$acc kernels
3294 ens_adj(:)=1.
3295!$acc end kernels
3296 xff_dicycle = 0.
3297
3298!--- large scale forcing
3299!
3300!$acc kernels
3301!$acc loop private(xff_ens3,xk)
3302 do 100 i=its,itf
3303 kloc(i)=1
3304 if(ierr(i).eq.0)then
3305! kloc(i)=maxloc(zu(i,:),1)
3306 kloc(i)=kbcon(i)
3307 ens_adj(i)=1.
3308!ss --- comment out adjustment over ocean
3309!ss if(ierr2(i).gt.0.and.ierr3(i).eq.0)ens_adj(i)=0.666 ! 2./3.
3310!ss if(ierr2(i).gt.0.and.ierr3(i).gt.0)ens_adj(i)=0.333
3311!
3312 a_ave=0.
3313 a_ave=axx(i)
3314 a_ave=max(0.,a_ave)
3315 a_ave=min(a_ave,aa1(i))
3316 a_ave=max(0.,a_ave)
3317 xff_ens3(:)=0.
3318 xff0= (aa1(i)-aa0(i))/dtime
3319 xff_ens3(1)=max(0.,(aa1(i)-aa0(i))/dtime)
3320 xff_ens3(2)=max(0.,(aa1(i)-aa0(i))/dtime)
3321 xff_ens3(3)=max(0.,(aa1(i)-aa0(i))/dtime)
3322 xff_ens3(16)=max(0.,(aa1(i)-aa0(i))/dtime)
3323 forcing(i,1)=xff_ens3(2)
3324!
3325!--- omeg is in bar/s, mconv done with omeg in pa/s
3326! more like brown (1979), or frank-cohen (199?)
3327!
3328! average aaround kbcon
3329!
3330 xomg=0.
3331 kk=0
3332 xff_ens3(4)=0.
3333 xff_ens3(5)=0.
3334 xff_ens3(6)=0.
3335 do k=kbcon(i)-1,kbcon(i)+1
3336 if(zu(i,k).gt.0.)then
3337 xomg=xomg-omeg(i,k)/9.81/max(0.3,(1.-(edt(i)*zd(i,k)-edtm(i)*zdm(i,k))/zu(i,k)))
3338 kk=kk+1
3339 endif
3340 enddo
3341 if(kk.gt.0)xff_ens3(4)=xomg/float(kk)
3342
3343!
3344! max below kbcon
3345! xff_ens3(6)=-omeg(i,k22(i))/9.81
3346! do k=k22(i),kbcon(i)
3347! xomg=-omeg(i,k)/9.81
3348! if(xomg.gt.xff_ens3(6))xff_ens3(6)=xomg
3349! enddo
3350!
3351! if(zu(i,kbcon(i)) > 0)xff_ens3(6)=betajb*xff_ens3(6)/zu(i,kbcon(i))
3352 xff_ens3(4)=betajb*xff_ens3(4)
3353 xff_ens3(5)=xff_ens3(4)
3354 xff_ens3(6)=xff_ens3(4)
3355 forcing(i,2)=xff_ens3(4)
3356 if(xff_ens3(4).lt.0.)xff_ens3(4)=0.
3357 if(xff_ens3(5).lt.0.)xff_ens3(5)=0.
3358 if(xff_ens3(6).lt.0.)xff_ens3(6)=0.
3359 xff_ens3(14)=xff_ens3(4)
3360!
3361!--- more like krishnamurti et al.; pick max and average values
3362!
3363 xff_ens3(7)= mconv(i) !/max(0.5,(1.-edt(i)*zd(i,kbcon(i))/zu(i,kloc(i))))
3364 xff_ens3(8)= mconv(i) !/max(0.5,(1.-edt(i)*zd(i,kbcon(i))/zu(i,kloc(i))))
3365 xff_ens3(9)= mconv(i) !/max(0.5,(1.-edt(i)*zd(i,kbcon(i))/zu(i,kloc(i))))
3366 xff_ens3(15)=mconv(i) !/max(0.5,(1.-edt(i)*zd(i,kbcon(i))/zu(i,kloc(i))))
3367 forcing(i,3)=xff_ens3(8)
3368!
3369!--- more like fritsch chappel or kain fritsch (plus triggers)
3370!
3371 xff_ens3(10)=aa1(i)/tau_ecmwf(i)
3372 xff_ens3(11)=aa1(i)/tau_ecmwf(i)
3373 xff_ens3(12)=aa1(i)/tau_ecmwf(i)
3374 xff_ens3(13)=(aa1(i))/tau_ecmwf(i) !(60.*15.) !tau_ecmwf(i)
3375 forcing(i,4)=xff_ens3(10)
3376! forcing(i,5)= aa1_bl(i)/tau_ecmwf(i)
3377
3378!!- more like bechtold et al. (jas 2014)
3379!! if(dicycle == 1) xff_dicycle = max(0.,aa1_bl(i)/tau_ecmwf(i)) !(60.*30.) !tau_ecmwf(i)
3380!gtest
3381 if(ichoice.eq.0)then
3382 if(xff0.lt.0.)then
3383 xff_ens3(1)=0.
3384 xff_ens3(2)=0.
3385 xff_ens3(3)=0.
3386 xff_ens3(10)=0.
3387 xff_ens3(11)=0.
3388 xff_ens3(12)=0.
3389 xff_ens3(13)= 0.
3390 xff_ens3(16)= 0.
3391! closure_n(i)=12.
3392! xff_dicycle = 0.
3393 endif !xff0
3394 endif ! ichoice
3395
3396 xk(1)=(xaa0(i,1)-aa1(i))/mbdt
3397 forcing(i,8)=mbdt*xk(1)/aa1(i)
3398! if(forcing(i,1).lt.0. .or. forcing(i,8).gt.-4.)ierr(i)=333
3399! if(forcing(i,2).lt.-0.05)ierr(i)=333
3400! forcing(i,4)=aa0(i)
3401! forcing(i,5)=aa1(i)
3402! forcing(i,6)=xaa0(i,1)
3403! forcing(i,7)=xk(1)
3404 if(xk(1).lt.0.and.xk(1).gt.-.01*mbdt) &
3405 xk(1)=-.01*mbdt
3406 if(xk(1).ge.0.and.xk(1).lt.1.e-2) &
3407 xk(1)=1.e-2
3408 ! enddo
3409!
3410!--- add up all ensembles
3411!
3412!
3413! over water, enfor!e small cap for some of the closures
3414!
3415 if(xland(i).lt.0.1)then
3416 if(ierr2(i).gt.0.or.ierr3(i).gt.0)then
3417 xff_ens3(1) =ens_adj(i)*xff_ens3(1)
3418 xff_ens3(2) =ens_adj(i)*xff_ens3(2)
3419 xff_ens3(3) =ens_adj(i)*xff_ens3(3)
3420 xff_ens3(4) =ens_adj(i)*xff_ens3(4)
3421 xff_ens3(5) =ens_adj(i)*xff_ens3(5)
3422 xff_ens3(6) =ens_adj(i)*xff_ens3(6)
3423 xff_ens3(7) =ens_adj(i)*xff_ens3(7)
3424 xff_ens3(8) =ens_adj(i)*xff_ens3(8)
3425 xff_ens3(9) =ens_adj(i)*xff_ens3(9)
3426 xff_ens3(10) =ens_adj(i)*xff_ens3(10)
3427 xff_ens3(11) =ens_adj(i)*xff_ens3(11)
3428 xff_ens3(12) =ens_adj(i)*xff_ens3(12)
3429 xff_ens3(13) =ens_adj(i)*xff_ens3(13)
3430 xff_ens3(14) =ens_adj(i)*xff_ens3(14)
3431 xff_ens3(15) =ens_adj(i)*xff_ens3(15)
3432 xff_ens3(16) =ens_adj(i)*xff_ens3(16)
3433!! !srf
3434!! xff_dicycle = ens_adj(i)*xff_dicycle
3435!! !srf end
3436! xff_ens3(7) =0.
3437! xff_ens3(8) =0.
3438! xff_ens3(9) =0.
3439 endif ! ierr2
3440 endif ! xland
3441!
3442! end water treatment
3443!
3444!
3445
3446!
3447!--- special treatment for stability closures
3448!
3449 if(xk(1).lt.0.)then
3450 if(xff_ens3(1).gt.0)xf_ens(i,1)=max(0.,-xff_ens3(1)/xk(1))
3451 if(xff_ens3(2).gt.0)xf_ens(i,2)=max(0.,-xff_ens3(2)/xk(1))
3452 if(xff_ens3(3).gt.0)xf_ens(i,3)=max(0.,-xff_ens3(3)/xk(1))
3453 if(xff_ens3(16).gt.0)xf_ens(i,16)=max(0.,-xff_ens3(16)/xk(1))
3454 xf_ens(i,1)= xf_ens(i,1)+xf_ens(i,1)*rand_clos(i,1)
3455 xf_ens(i,2)= xf_ens(i,2)+xf_ens(i,2)*rand_clos(i,1)
3456 xf_ens(i,3)= xf_ens(i,3)+xf_ens(i,3)*rand_clos(i,1)
3457 xf_ens(i,16)=xf_ens(i,16)+xf_ens(i,16)*rand_clos(i,1)
3458 else
3459 xff_ens3(1)= 0
3460 xff_ens3(2)= 0
3461 xff_ens3(3)= 0
3462 xff_ens3(16)=0
3463 endif
3464!
3465!--- if iresult.eq.1, following independent of xff0
3466!
3467 xf_ens(i,4)=max(0.,xff_ens3(4))
3468 xf_ens(i,5)=max(0.,xff_ens3(5))
3469 xf_ens(i,6)=max(0.,xff_ens3(6))
3470 xf_ens(i,14)=max(0.,xff_ens3(14))
3471 a1=max(1.e-3,pr_ens(i,7))
3472 xf_ens(i,7)=max(0.,xff_ens3(7)/a1)
3473 a1=max(1.e-3,pr_ens(i,8))
3474 xf_ens(i,8)=max(0.,xff_ens3(8)/a1)
3475! forcing(i,7)=xf_ens(i,8)
3476 a1=max(1.e-3,pr_ens(i,9))
3477 xf_ens(i,9)=max(0.,xff_ens3(9)/a1)
3478 a1=max(1.e-3,pr_ens(i,15))
3479 xf_ens(i,15)=max(0.,xff_ens3(15)/a1)
3480 xf_ens(i,4)=xf_ens(i,4)+xf_ens(i,4)*rand_clos(i,2)
3481 xf_ens(i,5)=xf_ens(i,5)+xf_ens(i,5)*rand_clos(i,2)
3482 xf_ens(i,6)=xf_ens(i,6)+xf_ens(i,6)*rand_clos(i,2)
3483 xf_ens(i,14)=xf_ens(i,14)+xf_ens(i,14)*rand_clos(i,2)
3484 xf_ens(i,7)=xf_ens(i,7)+xf_ens(i,7)*rand_clos(i,3)
3485 xf_ens(i,8)=xf_ens(i,8)+xf_ens(i,8)*rand_clos(i,3)
3486 xf_ens(i,9)=xf_ens(i,9)+xf_ens(i,9)*rand_clos(i,3)
3487 xf_ens(i,15)=xf_ens(i,15)+xf_ens(i,15)*rand_clos(i,3)
3488 if(xk(1).lt.0.)then
3489 xf_ens(i,10)=max(0.,-xff_ens3(10)/xk(1))
3490 xf_ens(i,11)=max(0.,-xff_ens3(11)/xk(1))
3491 xf_ens(i,12)=max(0.,-xff_ens3(12)/xk(1))
3492 xf_ens(i,13)=max(0.,-xff_ens3(13)/xk(1))
3493 xf_ens(i,10)=xf_ens(i,10)+xf_ens(i,10)*rand_clos(i,4)
3494 xf_ens(i,11)=xf_ens(i,11)+xf_ens(i,11)*rand_clos(i,4)
3495 xf_ens(i,12)=xf_ens(i,12)+xf_ens(i,12)*rand_clos(i,4)
3496 xf_ens(i,13)=xf_ens(i,13)+xf_ens(i,13)*rand_clos(i,4)
3497! forcing(i,8)=xf_ens(i,11)
3498 else
3499 xf_ens(i,10)=0.
3500 xf_ens(i,11)=0.
3501 xf_ens(i,12)=0.
3502 xf_ens(i,13)=0.
3503 !forcing(i,8)=0.
3504 endif
3505!srf-begin
3506!! if(xk(1).lt.0.)then
3507!! xf_dicycle(i) = max(0.,-xff_dicycle /xk(1))
3508!! forcing(i,9)=xf_dicycle(i)
3509!! else
3510!! xf_dicycle(i) = 0.
3511!! endif
3512!srf-end
3513 if(ichoice.ge.1)then
3514! closure_n(i)=0.
3515 xf_ens(i,1)=xf_ens(i,ichoice)
3516 xf_ens(i,2)=xf_ens(i,ichoice)
3517 xf_ens(i,3)=xf_ens(i,ichoice)
3518 xf_ens(i,4)=xf_ens(i,ichoice)
3519 xf_ens(i,5)=xf_ens(i,ichoice)
3520 xf_ens(i,6)=xf_ens(i,ichoice)
3521 xf_ens(i,7)=xf_ens(i,ichoice)
3522 xf_ens(i,8)=xf_ens(i,ichoice)
3523 xf_ens(i,9)=xf_ens(i,ichoice)
3524 xf_ens(i,10)=xf_ens(i,ichoice)
3525 xf_ens(i,11)=xf_ens(i,ichoice)
3526 xf_ens(i,12)=xf_ens(i,ichoice)
3527 xf_ens(i,13)=xf_ens(i,ichoice)
3528 xf_ens(i,14)=xf_ens(i,ichoice)
3529 xf_ens(i,15)=xf_ens(i,ichoice)
3530 xf_ens(i,16)=xf_ens(i,ichoice)
3531 endif
3532 elseif(ierr(i).ne.20.and.ierr(i).ne.0)then
3533 do n=1,maxens3
3534 xf_ens(i,n)=0.
3535!!
3536!! xf_dicycle(i) = 0.
3537!!
3538 enddo
3539 endif ! ierror
3540 100 continue
3541 !$acc end kernels
3542
3543
3544!-
3545!- diurnal cycle mass flux
3546!-
3547if(dicycle == 1 )then
3548!$acc kernels
3549!$acc loop private(xk)
3550 do i=its,itf
3551 xf_dicycle(i) = 0.
3552 if(ierr(i) /= 0)cycle
3553
3554 xk(1)=(xaa0(i,1)-aa1(i))/mbdt
3555! forcing(i,8)=xk(1)
3556 if(xk(1).lt.0.and.xk(1).gt.-.01*mbdt) xk(1)=-.01*mbdt
3557 if(xk(1).ge.0.and.xk(1).lt.1.e-2) xk(1)=1.e-2
3558
3559 xff_dicycle = (aa1(i)-aa1_bl(i))/tau_ecmwf(i)
3560! forcing(i,8)=xff_dicycle
3561 if(xk(1).lt.0) xf_dicycle(i)= max(0.,-xff_dicycle/xk(1))
3562
3563 xf_dicycle(i)= xf_ens(i,10)-xf_dicycle(i)
3564! forcing(i,6)=xf_dicycle(i)
3565 enddo
3566!$acc end kernels
3567else
3568!$acc kernels
3569 xf_dicycle(:) = 0.
3570!$acc end kernels
3571endif
3572!---------
3573
3574
3575
3576 end subroutine cup_forcing_ens_3d
3577
3579 subroutine cup_kbcon(ierrc,cap_inc,iloop_in,k22,kbcon,he_cup,hes_cup, &
3580 hkb,ierr,kbmax,p_cup,cap_max, &
3581 ztexec,zqexec, &
3582 jprnt,itf,ktf, &
3583 its,ite, kts,kte, &
3584 z_cup,entr_rate,heo,imid )
3585
3586 implicit none
3587!
3588
3589 ! only local wrf dimensions are need as of now in this routine
3590
3591 integer &
3592 ,intent (in ) :: &
3593 jprnt,itf,ktf,imid, &
3594 its,ite, kts,kte
3595 !
3596 !
3597 !
3598 ! ierr error value, maybe modified in this routine
3599 !
3600 real(kind=kind_phys), dimension (its:ite,kts:kte) &
3601 ,intent (in ) :: &
3602 he_cup,hes_cup,p_cup
3603!$acc declare copyin(he_cup,hes_cup,p_cup)
3604 real(kind=kind_phys), dimension (its:ite) &
3605 ,intent (in ) :: &
3606 entr_rate,ztexec,zqexec,cap_inc,cap_max
3607!$acc declare copyin(entr_rate,ztexec,zqexec,cap_inc,cap_max)
3608 real(kind=kind_phys), dimension (its:ite) &
3609 ,intent (inout ) :: &
3610 hkb !,cap_max
3611!$acc declare copy(hkb)
3612 integer, dimension (its:ite) &
3613 ,intent (in ) :: &
3614 kbmax
3615!$acc declare copyin(kbmax)
3616 integer, dimension (its:ite) &
3617 ,intent (inout) :: &
3618 kbcon,k22,ierr
3619!$acc declare copy(kbcon,k22,ierr)
3620 integer &
3621 ,intent (in ) :: &
3622 iloop_in
3623 character*50 :: ierrc(its:ite)
3624 real(kind=kind_phys), dimension (its:ite,kts:kte),intent (in) :: z_cup,heo
3625!$acc declare copyin(z_cup,heo)
3626 integer, dimension (its:ite) :: iloop,start_level
3627!$acc declare create(iloop,start_level)
3628!
3629! local variables in this routine
3630!
3631
3632 integer :: &
3633 i,k
3634 real(kind=kind_phys) :: &
3635 x_add,pbcdif,plus,hetest,dz
3636 real(kind=kind_phys), dimension (its:ite,kts:kte) ::hcot
3637!$acc declare create(hcot)
3638
3639!
3640!--- determine the level of convective cloud base - kbcon
3641!
3642!$acc kernels
3643 iloop(:)=iloop_in
3644!$acc end kernels
3645
3646!$acc parallel loop
3647 do 27 i=its,itf
3648 kbcon(i)=1
3649!
3650! reset iloop for mid level convection
3651 if(cap_max(i).gt.200 .and. imid.eq.1)iloop(i)=5
3652!
3653 if(ierr(i).ne.0)go to 27
3654 start_level(i)=k22(i)
3655 kbcon(i)=k22(i)+1
3656 if(iloop(i).eq.5)kbcon(i)=k22(i)
3657! if(iloop_in.eq.5)start_level(i)=kbcon(i)
3658 !== including entrainment for hetest
3659 hcot(i,1:start_level(i)) = hkb(i)
3660!$acc loop seq
3661 do k=start_level(i)+1,kbmax(i)+3
3662 dz=z_cup(i,k)-z_cup(i,k-1)
3663 hcot(i,k)= ( (1.-0.5*entr_rate(i)*dz)*hcot(i,k-1) &
3664 + entr_rate(i)*dz*heo(i,k-1) )/ &
3665 (1.+0.5*entr_rate(i)*dz)
3666 enddo
3667 !==
3668
3669 go to 32
3670 31 continue
3671 kbcon(i)=kbcon(i)+1
3672 if(kbcon(i).gt.kbmax(i)+2)then
3673 if(iloop(i).ne.4)then
3674 ierr(i)=3
3675#ifndef _OPENACC
3676 ierrc(i)="could not find reasonable kbcon in cup_kbcon"
3677#endif
3678 endif
3679 go to 27
3680 endif
3681 32 continue
3682 hetest=hcot(i,kbcon(i)) !hkb(i) ! he_cup(i,k22(i))
3683 if(hetest.lt.hes_cup(i,kbcon(i)))then
3684 go to 31
3685 endif
3686
3687! cloud base pressure and max moist static energy pressure
3688! i.e., the depth (in mb) of the layer of negative buoyancy
3689 if(kbcon(i)-k22(i).eq.1)go to 27
3690 if(iloop(i).eq.5 .and. (kbcon(i)-k22(i)).le.2)go to 27
3691 pbcdif=-p_cup(i,kbcon(i))+p_cup(i,k22(i))
3692 plus=max(25.,cap_max(i)-float(iloop(i)-1)*cap_inc(i))
3693 if(iloop(i).eq.4)plus=cap_max(i)
3694!
3695! for shallow convection, if cap_max is greater than 25, it is the pressure at pbltop
3696 if(iloop(i).eq.5)plus=150.
3697 if(iloop(i).eq.5.and.cap_max(i).gt.200)pbcdif=-p_cup(i,kbcon(i))+cap_max(i)
3698 if(pbcdif.le.plus)then
3699 go to 27
3700 elseif(pbcdif.gt.plus)then
3701 k22(i)=k22(i)+1
3702 kbcon(i)=k22(i)+1
3703!== since k22 has be changed, hkb has to be re-calculated
3704 x_add = xlv*zqexec(i)+cp*ztexec(i)
3705 call get_cloud_bc(kte,he_cup(i,1:kte),hkb(i),k22(i),x_add)
3706
3707 start_level(i)=k22(i)
3708! if(iloop_in.eq.5)start_level(i)=kbcon(i)
3709 hcot(i,1:start_level(i)) = hkb(i)
3710!$acc loop seq
3711 do k=start_level(i)+1,kbmax(i)+3
3712 dz=z_cup(i,k)-z_cup(i,k-1)
3713
3714 hcot(i,k)= ( (1.-0.5*entr_rate(i)*dz)*hcot(i,k-1) &
3715 + entr_rate(i)*dz*heo(i,k-1) )/ &
3716 (1.+0.5*entr_rate(i)*dz)
3717 enddo
3718 !==
3719
3720 if(iloop(i).eq.5)kbcon(i)=k22(i)
3721 if(kbcon(i).gt.kbmax(i)+2)then
3722 if(iloop(i).ne.4)then
3723 ierr(i)=3
3724#ifndef _OPENACC
3725 ierrc(i)="could not find reasonable kbcon in cup_kbcon"
3726#endif
3727 endif
3728 go to 27
3729 endif
3730 go to 32
3731 endif
3732 27 continue
3733 !$acc end parallel
3734
3735 end subroutine cup_kbcon
3736
3739 subroutine cup_maximi(array,ks,ke,maxx,ierr, &
3740 itf,ktf, &
3741 its,ite, kts,kte )
3742
3743 implicit none
3744!
3745! on input
3746!
3747
3748 ! only local wrf dimensions are need as of now in this routine
3749
3750 integer &
3751 ,intent (in ) :: &
3752 itf,ktf, &
3753 its,ite, kts,kte
3754 ! array input array
3755 ! x output array with return values
3756 ! kt output array of levels
3757 ! ks,kend check-range
3758 real(kind=kind_phys), dimension (its:ite,kts:kte) &
3759 ,intent (in ) :: &
3760 array
3761!$acc declare copyin(array)
3762 integer, dimension (its:ite) &
3763 ,intent (in ) :: &
3764 ierr,ke
3765!$acc declare copyin(ierr,ke)
3766 integer &
3767 ,intent (in ) :: &
3768 ks
3769 integer, dimension (its:ite) &
3770 ,intent (out ) :: &
3771 maxx
3772!$acc declare copyout(maxx)
3773 real(kind=kind_phys), dimension (its:ite) :: &
3774 x
3775!$acc declare create(x)
3776 real(kind=kind_phys) :: &
3777 xar
3778 integer :: &
3779 i,k
3780
3781!$acc kernels
3782 do 200 i=its,itf
3783 maxx(i)=ks
3784 if(ierr(i).eq.0)then
3785 x(i)=array(i,ks)
3786!
3787!$acc loop seq
3788 do 100 k=ks,ke(i)
3789 xar=array(i,k)
3790 if(xar.ge.x(i)) then
3791 x(i)=xar
3792 maxx(i)=k
3793 endif
3794 100 continue
3795 endif
3796 200 continue
3797 !$acc end kernels
3798
3799 end subroutine cup_maximi
3800
3802 subroutine cup_minimi(array,ks,kend,kt,ierr, &
3803 itf,ktf, &
3804 its,ite, kts,kte )
3805
3806 implicit none
3807!
3808! on input
3809!
3810
3811 ! only local wrf dimensions are need as of now in this routine
3812
3813 integer &
3814 ,intent (in ) :: &
3815 itf,ktf, &
3816 its,ite, kts,kte
3817 ! array input array
3818 ! x output array with return values
3819 ! kt output array of levels
3820 ! ks,kend check-range
3821 real(kind=kind_phys), dimension (its:ite,kts:kte) &
3822 ,intent (in ) :: &
3823 array
3824!$acc declare copyin(array)
3825 integer, dimension (its:ite) &
3826 ,intent (in ) :: &
3827 ierr,ks,kend
3828!$acc declare copyin(ierr,ks,kend)
3829 integer, dimension (its:ite) &
3830 ,intent (out ) :: &
3831 kt
3832!$acc declare copyout(kt)
3833 real(kind=kind_phys), dimension (its:ite) :: &
3834 x
3835!$acc declare create(x)
3836 integer :: &
3837 i,k,kstop
3838
3839!$acc kernels
3840 do 200 i=its,itf
3841 kt(i)=ks(i)
3842 if(ierr(i).eq.0)then
3843 x(i)=array(i,ks(i))
3844 kstop=max(ks(i)+1,kend(i))
3845!
3846!$acc loop seq
3847 do 100 k=ks(i)+1,kstop
3848 if(array(i,k).lt.x(i)) then
3849 x(i)=array(i,k)
3850 kt(i)=k
3851 endif
3852 100 continue
3853 endif
3854 200 continue
3855 !$acc end kernels
3856
3857 end subroutine cup_minimi
3858
3860 subroutine cup_up_aa0(aa0,z,zu,dby,gamma_cup,t_cup, &
3861 kbcon,ktop,ierr, &
3862 itf,ktf, &
3863 its,ite, kts,kte )
3864
3865 implicit none
3866!
3867! on input
3868!
3869
3870 ! only local wrf dimensions are need as of now in this routine
3871
3872 integer &
3873 ,intent (in ) :: &
3874 itf,ktf, &
3875 its,ite, kts,kte
3876 ! aa0 cloud work function
3877 ! gamma_cup = gamma on model cloud levels
3878 ! t_cup = temperature (kelvin) on model cloud levels
3879 ! dby = buoancy term
3880 ! zu= normalized updraft mass flux
3881 ! z = heights of model levels
3882 ! ierr error value, maybe modified in this routine
3883 !
3884 real(kind=kind_phys), dimension (its:ite,kts:kte) &
3885 ,intent (in ) :: &
3886 z,zu,gamma_cup,t_cup,dby
3887 integer, dimension (its:ite) &
3888 ,intent (in ) :: &
3889 kbcon,ktop
3890!$acc declare copyin(z,zu,gamma_cup,t_cup,dby,kbcon,ktop)
3891!
3892! input and output
3893!
3894
3895
3896 integer, dimension (its:ite) &
3897 ,intent (inout) :: &
3898 ierr
3899!$acc declare copy(ierr)
3900 real(kind=kind_phys), dimension (its:ite) &
3901 ,intent (out ) :: &
3902 aa0
3903!$acc declare copyout(aa0)
3904!
3905! local variables in this routine
3906!
3907
3908 integer :: &
3909 i,k
3910 real(kind=kind_phys) :: &
3911 dz,da
3912!
3913!$acc kernels
3914 do i=its,itf
3915 aa0(i)=0.
3916 enddo
3917 do k=kts+1,ktf
3918 do i=its,itf
3919 if(ierr(i).ne.0) cycle
3920 if(k.lt.kbcon(i)) cycle
3921 if(k.gt.ktop(i)) cycle
3922 dz=z(i,k)-z(i,k-1)
3923 da=zu(i,k)*dz*(9.81/(1004.*( &
3924 (t_cup(i,k)))))*dby(i,k-1)/ &
3925 (1.+gamma_cup(i,k))
3926 ! if(k.eq.ktop(i).and.da.le.0.)go to 100
3927 aa0(i)=aa0(i)+max(0.,da)
3928 if(aa0(i).lt.0.)aa0(i)=0.
3929 enddo
3930 enddo
3931!$acc end kernels
3932
3933 end subroutine cup_up_aa0
3934
3935!====================================================================
3936
3939 subroutine neg_check(name,j,dt,q,outq,outt,outu,outv, &
3940 outqc,pret,its,ite,kts,kte,itf,ktf,ktop)
3941
3942 integer, intent(in ) :: j,its,ite,kts,kte,itf,ktf
3943 integer, dimension (its:ite ), intent(in ) :: ktop
3944
3945 real(kind=kind_phys), dimension (its:ite,kts:kte ) , &
3946 intent(inout ) :: &
3947 outq,outt,outqc,outu,outv
3948 real(kind=kind_phys), dimension (its:ite,kts:kte ) , &
3949 intent(inout ) :: &
3950 q
3951 real(kind=kind_phys), dimension (its:ite ) , &
3952 intent(inout ) :: &
3953 pret
3954!$acc declare copy(outq,outt,outqc,outu,outv,q,pret)
3955 character *(*), intent (in) :: &
3956 name
3957 real(kind=kind_phys) &
3958 ,intent (in ) :: &
3959 dt
3960 real(kind=kind_phys) :: names,scalef,thresh,qmem,qmemf,qmem2,qtest,qmem1
3961 integer :: icheck
3962!
3963! first do check on vertical heating rate
3964!
3965 thresh=300.01
3966! thresh=200.01 !ss
3967! thresh=250.01
3968 names=1.
3969 if(name == 'shallow' .or. name == 'mid')then
3970 thresh=148.01
3971 names=1.
3972 endif
3973 scalef=86400.
3974!$acc kernels
3975!$acc loop private(qmemf,qmem,icheck)
3976 do i=its,itf
3977 if(ktop(i) <= 2)cycle
3978 icheck=0
3979 qmemf=1.
3980 qmem=0.
3981!$acc loop reduction(min:qmemf)
3982 do k=kts,ktop(i)
3983 qmem=(outt(i,k))*86400.
3984 if(qmem.gt.thresh)then
3985 qmem2=thresh/qmem
3986 qmemf=min(qmemf,qmem2)
3987 icheck=1
3988!
3989!
3990! print *,'1',' adjusted massflux by factor ',i,j,k,qmem,qmem2,qmemf,dt
3991 endif
3992 if(qmem.lt.-.5*thresh*names)then
3993 qmem2=-.5*names*thresh/qmem
3994 qmemf=min(qmemf,qmem2)
3995 icheck=2
3996!
3997!
3998 endif
3999 enddo
4000 do k=kts,ktop(i)
4001 outq(i,k)=outq(i,k)*qmemf
4002 outt(i,k)=outt(i,k)*qmemf
4003 outu(i,k)=outu(i,k)*qmemf
4004 outv(i,k)=outv(i,k)*qmemf
4005 outqc(i,k)=outqc(i,k)*qmemf
4006 enddo
4007 pret(i)=pret(i)*qmemf
4008 enddo
4009!$acc end kernels
4010! return
4011!
4012! check whether routine produces negative q's. this can happen, since
4013! tendencies are calculated based on forced q's. this should have no
4014! influence on conservation properties, it scales linear through all
4015! tendencies
4016!
4017! return
4018! write(14,*)'return'
4019 thresh=1.e-32
4020!$acc kernels
4021!$acc loop private(qmemf,qmem,icheck)
4022 do i=its,itf
4023 if(ktop(i) <= 2)cycle
4024 qmemf=1.
4025!$acc loop reduction(min:qmemf)
4026 do k=kts,ktop(i)
4027 qmem=outq(i,k)
4028 if(abs(qmem).gt.0. .and. q(i,k).gt.1.e-6)then
4029 qtest=q(i,k)+(outq(i,k))*dt
4030 if(qtest.lt.thresh)then
4031!
4032! qmem2 would be the maximum allowable tendency
4033!
4034 qmem1=abs(outq(i,k))
4035 qmem2=abs((thresh-q(i,k))/dt)
4036 qmemf=min(qmemf,qmem2/qmem1)
4037 qmemf=max(0.,qmemf)
4038 endif
4039 endif
4040 enddo
4041 do k=kts,ktop(i)
4042 outq(i,k)=outq(i,k)*qmemf
4043 outt(i,k)=outt(i,k)*qmemf
4044 outu(i,k)=outu(i,k)*qmemf
4045 outv(i,k)=outv(i,k)*qmemf
4046 outqc(i,k)=outqc(i,k)*qmemf
4047 enddo
4048 pret(i)=pret(i)*qmemf
4049 enddo
4050!$acc end kernels
4051 end subroutine neg_check
4052
4055 subroutine cup_output_ens_3d(xff_mid,xf_ens,ierr,dellat,dellaq,dellaqc, &
4056 outtem,outq,outqc,dx, &
4057 zu,pre,pw,xmb,ktop, &
4058 edt,pwd,name,ierr2,ierr3,p_cup,pr_ens, &
4059 maxens3, &
4060 sig,closure_n,xland1,xmbm_in,xmbs_in, &
4061 ichoice,imid,ipr,itf,ktf, &
4062 its,ite, kts,kte, &
4063 dicycle,xf_dicycle )
4064
4065 implicit none
4066!
4067! on input
4068!
4069 ! only local wrf dimensions are need as of now in this routine
4070
4071 integer &
4072 ,intent (in ) :: &
4073 ichoice,imid,ipr,itf,ktf, &
4074 its,ite, kts,kte
4075 integer, intent (in ) :: &
4076 maxens3
4077 ! xf_ens = ensemble mass fluxes
4078 ! pr_ens = precipitation ensembles
4079 ! dellat = change of temperature per unit mass flux of cloud ensemble
4080 ! dellaq = change of q per unit mass flux of cloud ensemble
4081 ! dellaqc = change of qc per unit mass flux of cloud ensemble
4082 ! outtem = output temp tendency (per s)
4083 ! outq = output q tendency (per s)
4084 ! outqc = output qc tendency (per s)
4085 ! pre = output precip
4086 ! xmb = total base mass flux
4087 ! xfac1 = correction factor
4088 ! pw = pw -epsilon*pd (ensemble dependent)
4089 ! ierr error value, maybe modified in this routine
4090 !
4091 real(kind=kind_phys), dimension (its:ite,1:maxens3) &
4092 ,intent (inout) :: &
4093 xf_ens,pr_ens
4094 real(kind=kind_phys), dimension (its:ite,kts:kte) &
4095 ,intent (inout ) :: &
4096 outtem,outq,outqc
4097 real(kind=kind_phys), dimension (its:ite,kts:kte) &
4098 ,intent (in ) :: &
4099 zu,pwd,p_cup
4100 real(kind=kind_phys), dimension (its:ite) &
4101 ,intent (in ) :: &
4102 sig,xmbm_in,xmbs_in,edt,dx
4103 real(kind=kind_phys), dimension (its:ite,2) &
4104 ,intent (in ) :: &
4105 xff_mid
4106 real(kind=kind_phys), dimension (its:ite) &
4107 ,intent (inout ) :: &
4108 pre,xmb
4109 real(kind=kind_phys), dimension (its:ite) &
4110 ,intent (inout ) :: &
4111 closure_n
4112 real(kind=kind_phys), dimension (its:ite,kts:kte,1) &
4113 ,intent (in ) :: &
4114 dellat,dellaqc,dellaq,pw
4115 integer, dimension (its:ite) &
4116 ,intent (in ) :: &
4117 ktop,xland1
4118 integer, dimension (its:ite) &
4119 ,intent (inout) :: &
4120 ierr,ierr2,ierr3
4121 integer, intent(in) :: dicycle
4122 real(kind=kind_phys), intent(in), dimension (its:ite) :: xf_dicycle
4123!$acc declare copyin(zu,pwd,p_cup,sig,xmbm_in,xmbs_in,edt,xff_mid,dellat,dellaqc,dellaq,pw,ktop,xland1,xf_dicycle)
4124!$acc declare copy(xf_ens,pr_ens,outtem,outq,outqc,pre,xmb,closure_n,ierr,ierr2,ierr3)
4125!
4126! local variables in this routine
4127!
4128
4129 integer :: &
4130 i,k,n
4131 real(kind=kind_phys) :: &
4132 clos_wei,dtt,dp,dtq,dtqc,dtpw,dtpwd
4133 real(kind=kind_phys), dimension (its:ite) :: &
4134 pre2,xmb_ave,pwtot
4135!$acc declare create(pre2,xmb_ave,pwtot)
4136!
4137 character *(*), intent (in) :: &
4138 name
4139
4140!
4141!$acc kernels
4142 do k=kts,kte
4143 do i=its,ite
4144 outtem(i,k)=0.
4145 outq(i,k)=0.
4146 outqc(i,k)=0.
4147 enddo
4148 enddo
4149 do i=its,itf
4150 pre(i)=0.
4151 xmb(i)=0.
4152 enddo
4153 do i=its,itf
4154 if(ierr(i).eq.0)then
4155 do n=1,maxens3
4156 if(pr_ens(i,n).le.0.)then
4157 xf_ens(i,n)=0.
4158 endif
4159 enddo
4160 endif
4161 enddo
4162!$acc end kernels
4163!
4164!--- calculate ensemble average mass fluxes
4165!
4166
4167!
4168!-- now do feedback
4169!
4170!!!!! deep convection !!!!!!!!!!
4171 if(imid.eq.0)then
4172!$acc kernels
4173 do i=its,itf
4174 if(ierr(i).eq.0)then
4175 k=0
4176 xmb_ave(i)=0.
4177!$acc loop seq
4178 do n=1,maxens3
4179 k=k+1
4180 xmb_ave(i)=xmb_ave(i)+xf_ens(i,n)
4181
4182 enddo
4183 !print *,'xf_ens',xf_ens
4184 xmb_ave(i)=xmb_ave(i)/float(k)
4185 !print *,'k,xmb_ave',k,xmb_ave
4186 !srf begin
4187 if(dicycle == 2 )then
4188 xmb_ave(i)=xmb_ave(i)-max(0.,xmbs_in(i))
4189 xmb_ave(i)=max(0.,xmb_ave(i))
4190 else if (dicycle == 1) then
4191! xmb_ave(i)=min(xmb_ave(i),xmb_ave(i) - xf_dicycle(i))
4192 xmb_ave(i)=xmb_ave(i) - xf_dicycle(i)
4193 xmb_ave(i)=max(0.,xmb_ave(i))
4194 endif
4195 !print *,"2 xmb_ave,xf_dicycle",xmb_ave,xf_dicycle
4196! --- now use proper count of how many closures were actually
4197! used in cup_forcing_ens (including screening of some
4198! closures over water) to properly normalize xmb
4199 clos_wei=16./max(1.,closure_n(i))
4200 xmb_ave(i)=min(xmb_ave(i),100.)
4201 xmb(i)=clos_wei*sig(i)*xmb_ave(i)
4202
4203 if(xmb(i) < 1.e-16)then
4204 ierr(i)=19
4205 endif
4206! xfac1(i)=xmb(i)
4207! xfac2(i)=xmb(i)
4208
4209 endif
4210 enddo
4211!$acc end kernels
4212!!!!! not so deep convection !!!!!!!!!!
4213 else ! imid == 1
4214!$acc kernels
4215 do i=its,itf
4216 xmb_ave(i)=0.
4217 if(ierr(i).eq.0)then
4218! ! first get xmb_ves, depend on ichoicee
4219!
4220 if(ichoice.eq.1 .or. ichoice.eq.2)then
4221 xmb_ave(i)=sig(i)*xff_mid(i,ichoice)
4222 else if(ichoice.gt.2)then
4223 k=0
4224!$acc loop seq
4225 do n=1,maxens3
4226 k=k+1
4227 xmb_ave(i)=xmb_ave(i)+xf_ens(i,n)
4228 enddo
4229 xmb_ave(i)=xmb_ave(i)/float(k)
4230 else if(ichoice == 0)then
4231 xmb_ave(i)=.5*sig(i)*(xff_mid(i,1)+xff_mid(i,2))
4232 endif ! ichoice gt.2
4233! which dicycle method
4234 if(dicycle == 2 )then
4235 xmb(i)=max(0.,xmb_ave(i)-xmbs_in(i))
4236 else if (dicycle == 1) then
4237! xmb(i)=min(xmb_ave(i),xmb_ave(i) - xf_dicycle(i))
4238 xmb(i)=xmb_ave(i) - xf_dicycle(i)
4239 xmb(i)=max(0.,xmb_ave(i))
4240 else if (dicycle == 0) then
4241 xmb(i)=max(0.,xmb_ave(i))
4242 endif ! dicycle=1,2
4243 endif ! ierr >0
4244 enddo ! i
4245!$acc end kernels
4246 endif ! imid=1
4247
4248!$acc kernels
4249 do i=its,itf
4250 if(ierr(i).eq.0)then
4251 dtpw=0.
4252 do k=kts,ktop(i)
4253 dtpw=dtpw+pw(i,k,1)
4254 outtem(i,k)= xmb(i)* dellat(i,k,1)
4255 outq(i,k)= xmb(i)* dellaq(i,k,1)
4256 outqc(i,k)= xmb(i)* dellaqc(i,k,1)
4257 enddo
4258 pre(i)=pre(i)+xmb(i)*dtpw
4259 endif
4260 enddo
4261!$acc end kernels
4262 return
4263
4264!$acc kernels
4265 do i=its,itf
4266 pwtot(i)=0.
4267 pre2(i)=0.
4268 if(ierr(i).eq.0)then
4269 do k=kts,ktop(i)
4270 pwtot(i)=pwtot(i)+pw(i,k,1)
4271 enddo
4272 do k=kts,ktop(i)
4273 dp=100.*(p_cup(i,k)-p_cup(i,k+1))/g
4274 dtt =dellat(i,k,1)
4275 dtq =dellaq(i,k,1)
4276! necessary to drive downdraft
4277 dtpwd=-pwd(i,k)*edt(i)
4278! take from dellaqc first
4279 dtqc=dellaqc(i,k,1)*dp - dtpwd
4280! if this is negative, use dellaqc first, rest needs to come from rain
4281 if(dtqc < 0.)then
4282 dtpwd=dtpwd-dellaqc(i,k,1)*dp
4283 dtqc=0.
4284! if this is positive, can come from clw detrainment
4285 else
4286 dtqc=dtqc/dp
4287 dtpwd=0.
4288 endif
4289 outtem(i,k)= xmb(i)* dtt
4290 outq(i,k)= xmb(i)* dtq
4291 outqc(i,k)= xmb(i)* dtqc
4292 xf_ens(i,:)=sig(i)*xf_ens(i,:)
4293! what is evaporated
4294 pre(i)=pre(i)-xmb(i)*dtpwd
4295 pre2(i)=pre2(i)+xmb(i)*(pw(i,k,1)+edt(i)*pwd(i,k))
4296! write(15,124)k,dellaqc(i,k,1),dtqc,-pwd(i,k)*edt(i),dtpwd
4297 enddo
4298 pre(i)=-pre(i)+xmb(i)*pwtot(i)
4299 endif
4300#ifndef _OPENACC
4301124 format(1x,i3,4e13.4)
4302125 format(1x,2e13.4)
4303#endif
4304 enddo
4305!$acc end kernels
4306
4307 end subroutine cup_output_ens_3d
4308!-------------------------------------------------------
4310 subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, &
4311 p_cup,kbcon,ktop,dby,clw_all,xland1, &
4312 q,gamma_cup,zu,qes_cup,k22,qe_cup,c0,c0t3d, &
4313 zqexec,ccn,ccnclean,rho,c1d,t,autoconv, &
4314 up_massentr,up_massdetr,psum,psumh, &
4315 itest,itf,ktf, &
4316 its,ite, kts,kte )
4317
4318 implicit none
4319 real(kind=kind_phys), parameter :: bdispm = 0.366
4320 real(kind=kind_phys), parameter :: bdispc = 0.146
4321!
4322! on input
4323!
4324
4325 ! only local wrf dimensions are need as of now in this routine
4326
4327 integer &
4328 ,intent (in ) :: &
4329 autoconv, &
4330 itest,itf,ktf, &
4331 its,ite, kts,kte
4332 ! cd= detrainment function
4333 ! q = environmental q on model levels
4334 ! qe_cup = environmental q on model cloud levels
4335 ! qes_cup = saturation q on model cloud levels
4336 ! dby = buoancy term
4337 ! cd= detrainment function
4338 ! zu = normalized updraft mass flux
4339 ! gamma_cup = gamma on model cloud levels
4340 !
4341 real(kind=kind_phys), dimension (its:ite,kts:kte) &
4342 ,intent (in ) :: &
4343 p_cup,rho,q,zu,gamma_cup,qe_cup, &
4344 up_massentr,up_massdetr,dby,qes_cup,z_cup
4345 real(kind=kind_phys), dimension (its:ite) &
4346 ,intent (in ) :: &
4347 zqexec,c0
4348 real(kind=kind_phys), dimension (its:ite,kts:kte), intent (out) :: c0t3d
4349 ! entr= entrainment rate
4350 integer, dimension (its:ite) &
4351 ,intent (in ) :: &
4352 kbcon,ktop,k22,xland1
4353!$acc declare copyin(p_cup,rho,q,zu,gamma_cup,qe_cup,up_massentr,up_massdetr,dby,qes_cup,z_cup,zqexec,c0,kbcon,ktop,k22,xland1)
4354!$acc declare copy(c0t3d)
4355 real(kind=kind_phys), intent (in ) :: & ! HCB
4356 ccnclean
4357!
4358! input and output
4359!
4360
4361 ! ierr error value, maybe modified in this routine
4362
4363 integer, dimension (its:ite) &
4364 ,intent (inout) :: &
4365 ierr
4366!$acc declare copy(ierr)
4367 character *(*), intent (in) :: &
4368 name
4369 ! qc = cloud q (including liquid water) after entrainment
4370 ! qrch = saturation q in cloud
4371 ! qrc = liquid water content in cloud after rainout
4372 ! pw = condensate that will fall out at that level
4373 ! pwav = totan normalized integrated condensate (i1)
4374 ! c0 = conversion rate (cloud to rain)
4375
4376 real(kind=kind_phys), dimension (its:ite,kts:kte) &
4377 ,intent (out ) :: &
4378 qc,qrc,pw,clw_all
4379!$acc declare copy(qc,qrc,pw,clw_all)
4380 real(kind=kind_phys), dimension (its:ite,kts:kte) &
4381 ,intent (inout) :: &
4382 c1d
4383!$acc declare copy(c1d)
4384 real(kind=kind_phys), dimension (its:ite,kts:kte) :: &
4385 qch,qrcb,pwh,clw_allh,c1d_b,t
4386!$acc declare create(qch,qrcb,pwh,clw_allh,c1d_b,t)
4387 real(kind=kind_phys), dimension (its:ite) :: &
4388 pwavh
4389!$acc declare create(pwavh)
4390 real(kind=kind_phys), dimension (its:ite) &
4391 ,intent (out ) :: &
4392 pwav,psum,psumh
4393!$acc declare copyout(pwav,psum,psumh)
4394 real(kind=kind_phys), dimension (its:ite) &
4395 ,intent (in ) :: &
4396 ccn
4397!$acc declare copyin(ccn)
4398!
4399! local variables in this routine
4400!
4401
4402 integer :: &
4403 iprop,iall,i,k
4404 integer :: start_level(its:ite),kklev(its:ite)
4405!$acc declare create(start_level,kklev)
4406 real(kind=kind_phys) :: &
4407 prop_ave,qrcb_h,dp,rhoc,qrch,qaver,clwdet, &
4408 dz,berryc0,q1,berryc
4409 real(kind=kind_phys) :: &
4410 denom, c0t, c0_iceconv
4411 real(kind=kind_phys), dimension (kts:kte) :: &
4412 prop_b
4413 real(kind=kind_phys), dimension (its:ite) :: &
4414 bdsp
4415!$acc declare create(prop_b)
4416!
4417 real(kind=kind_phys), parameter:: zero = 0
4418 logical :: is_mid, is_deep
4419
4420 is_mid = (name == 'mid')
4421 is_deep = (name == 'deep')
4422
4423!$acc kernels
4424 prop_b(kts:kte)=0
4425!$acc end kernels
4426 iall=0
4427 clwdet=0.1 !0.02
4428 c0_iceconv=0.01
4429 c1d_b=c1d
4430 bdsp(:)=bdispm
4431!$acc kernels
4432 c0t3d = 0.
4433!$acc end kernels
4434
4435!
4436!--- no precip for small clouds
4437!
4438! if(name.eq.'shallow')then
4439! c0=0.002
4440! endif
4441!$acc kernels
4442 do i=its,itf
4443 pwav(i)=0.
4444 pwavh(i)=0.
4445 psum(i)=0.
4446 psumh(i)=0.
4447 if (xland1(i) .eq. 0) then
4448 bdsp(i)=bdispm
4449 else
4450 bdsp(i)=bdispc
4451 endif
4452 enddo
4453 do k=kts,ktf
4454 do i=its,itf
4455 pw(i,k)=0.
4456 pwh(i,k)=0.
4457 qc(i,k)=0.
4458 if(ierr(i).eq.0)qc(i,k)=qe_cup(i,k)
4459 if(ierr(i).eq.0)qch(i,k)=qe_cup(i,k)
4460 clw_all(i,k)=0.
4461 clw_allh(i,k)=0.
4462 qrc(i,k)=0.
4463 qrcb(i,k)=0.
4464 enddo
4465 enddo
4466!$acc end kernels
4467
4468!$acc parallel loop private(start_level,qaver,k)
4469 do i=its,itf
4470 if(ierr(i).eq.0)then
4471 start_level=k22(i)
4472 call get_cloud_bc(kte,qe_cup(i,1:kte),qaver,k22(i),zero)
4473 qaver = qaver
4474 k=start_level(i)
4475 qc(i,k)= qaver
4476 qch(i,k)= qaver
4477 do k=1,start_level(i)-1
4478 qc(i,k)= qe_cup(i,k)
4479 qch(i,k)= qe_cup(i,k)
4480 enddo
4481!
4482! initialize below originating air
4483!
4484 endif
4485 enddo
4486!$acc end parallel
4487
4488
4489!$acc kernels
4490 do 100 i=its,itf
4491 !c0=.004 HCB tuning
4492 if(ierr(i).eq.0)then
4493
4494! below lfc, but maybe above lcl
4495!
4496! if(name == "deep" )then
4497!$acc loop seq
4498 do k=k22(i)+1,kbcon(i)
4499 if(t(i,k) > 273.16) then
4500 c0t = c0(i)
4501 else
4502 c0t = c0(i) * exp(c0_iceconv * (t(i,k) - 273.16))
4503 endif
4504 c0t3d(i,k)=c0t
4505 qc(i,k)= (qc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)* qc(i,k-1)+ &
4506 up_massentr(i,k-1)*q(i,k-1)) / &
4507 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
4508! qrch=qes_cup(i,k)
4509 qrch=qes_cup(i,k)+(1./xlv)*(gamma_cup(i,k) &
4510 /(1.+gamma_cup(i,k)))*dby(i,k)
4511 if(k.lt.kbcon(i))qrch=qc(i,k)
4512 if(qc(i,k).gt.qrch)then
4513 dz=z_cup(i,k)-z_cup(i,k-1)
4514 qrc(i,k)=(qc(i,k)-qrch)/(1.+c0t*dz)
4515 pw(i,k)=c0t*dz*qrc(i,k)*zu(i,k)
4516 qc(i,k)=qrch+qrc(i,k)
4517 clw_all(i,k)=qrc(i,k)
4518 endif
4519 clw_allh(i,k)=clw_all(i,k)
4520 qrcb(i,k)=qrc(i,k)
4521 pwh(i,k)=pw(i,k)
4522 qch(i,k)=qc(i,k)
4523 enddo
4524 ! endif
4525!
4526!now do the rest
4527!
4528 kklev(i)=maxloc(zu(i,:),1)
4529!$acc loop seq
4530 do k=kbcon(i)+1,ktop(i)
4531 if(t(i,k) > 273.16) then
4532 c0t = c0(i)
4533 else
4534 c0t = c0(i) * exp(c0_iceconv * (t(i,k) - 273.16))
4535 endif
4536 if(is_mid)c0t=0.004
4537 c0t3d(i,k)=c0t
4538
4539 if(autoconv .gt.1) c0t=c0(i)
4540 denom=zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1)
4541 if(denom.lt.1.e-16)then
4542 ierr(i)=51
4543 exit
4544 endif
4545
4546
4547 rhoc=.5*(rho(i,k)+rho(i,k-1))
4548 dz=z_cup(i,k)-z_cup(i,k-1)
4549 dp=-100.*(p_cup(i,k)-p_cup(i,k-1))
4550!
4551!--- saturation in cloud, this is what is allowed to be in it
4552!
4553 qrch=qes_cup(i,k)+(1./xlv)*(gamma_cup(i,k) &
4554 /(1.+gamma_cup(i,k)))*dby(i,k)
4555!
4556!------ 1. steady state plume equation, for what could
4557!------ be in cloud without condensation
4558!
4559!
4560 qc(i,k)= (qc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)* qc(i,k-1)+ &
4561 up_massentr(i,k-1)*q(i,k-1)) / &
4562 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
4563 qch(i,k)= (qch(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*qch(i,k-1)+ &
4564 up_massentr(i,k-1)*q(i,k-1)) / &
4565 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
4566
4567 if(qc(i,k).le.qrch)then
4568 qc(i,k)=qrch+1e-8
4569 endif
4570 if(qch(i,k).le.qrch)then
4571 qch(i,k)=qrch+1e-8
4572 endif
4573!
4574!------- total condensed water before rainout
4575!
4576 clw_all(i,k)=max(0.,qc(i,k)-qrch)
4577 qrc(i,k)=max(0.,(qc(i,k)-qrch)) ! /(1.+c0(i)*dz*zu(i,k))
4578 clw_allh(i,k)=max(0.,qch(i,k)-qrch)
4579 qrcb(i,k)=max(0.,(qch(i,k)-qrch)) ! /(1.+c0(i)*dz*zu(i,k))
4580 if(is_deep)then
4581 clwdet=0.1 !0.02 ! 05/11/2021
4582 !if(k.lt.kklev(i)) clwdet=0. ! 05/05/2021
4583 else
4584 clwdet=0.1 !0.02 ! 05/05/2021
4585 !if(k.lt.kklev(i)) clwdet=0. ! 05/25/2021
4586 endif
4587 if(k.gt.kbcon(i)+1)c1d(i,k)=clwdet*up_massdetr(i,k-1)
4588 if(k.gt.kbcon(i)+1)c1d_b(i,k)=clwdet*up_massdetr(i,k-1)
4589
4590 if(autoconv.eq.2) then
4591!
4592! normalized berry
4593!
4594! first calculate for average conditions, used in cup_dd_edt!
4595! this will also determine proportionality constant prop_b, which, if applied,
4596! would give the same results as c0 under these conditions
4597!
4598! Berry conversion for clean atmosphere
4599!
4600 q1=1.e3*rhoc*clw_allh(i,k)
4601! pwh units are kg/kg, but normalized by mass flux. So with massflux kg/m^2/s
4602 pwh(i,k)=c0t*dz*zu(i,k)*clw_allh(i,k)
4603 qrcb_h=(qch(i,k)-qrch)/(1.+(c1d_b(i,k)+c0t)*dz)
4604 qrcb(i,k)=0.
4605! unit (B) = g/m^3/s
4606 berryc0=(q1*q1/(60.0*(5.0 + 0.0366*ccnclean*1.e1/ &
4607 ( q1 * bdsp(i)) ) ))
4608! normalize Berry: berryc0=berryc0*g/dp*dz*zu = pwh, unts become kg/kg
4609! set 1:
4610 berryc0=1.e-3*berryc0*g/dp*dz
4611 prop_b(k)=pwh(i,k)/berryc0
4612 qrcb(i,k)=qrcb_h
4613 if(qrcb(i,k).le.0.)then
4614 pwh(i,k)=0.
4615 endif
4616 qch(i,k)=qrcb(i,k)+qrch
4617 pwavh(i)=pwavh(i)+pwh(i,k)
4618 psumh(i)=psumh(i)+pwh(i,k)*g/dp !dz !dp/g !*dp ! HCB
4619! then the real berry
4620!
4621 q1=1.e3*rhoc*clw_all(i,k)
4622 berryc=(q1*q1/(60.0*(5.0 + 0.0366*ccn(i)*1.e1/ &
4623 ( q1 * bdsp(i)) ) ))
4624 berryc=1.e-3*berryc*g/dp*dz
4625 pw(i,k)=prop_b(k)*berryc !*dz/zu(i,k)
4626! use berryc now as new c0 for this level
4627 berryc=pw(i,k)/(dz*zu(i,k)*clw_all(i,k))
4628 if(qrc(i,k).le.0.)then
4629 berryc=0.
4630 endif
4631 qrc(i,k)=(max(0.,(qc(i,k)-qrch))/(1+(c1d(i,k)+berryc)*dz))
4632 if(qrc(i,k).lt.0.)then
4633 qrc(i,k)=0.
4634 pw(i,k)=0.
4635 endif
4636 qc(i,k)=qrc(i,k)+qrch
4637
4638! if not running with berry at all, do the following
4639!
4640 else
4641! create clw detrainment profile that depends on mass detrainment and
4642! in-cloud clw/ice
4643!
4644 qrc(i,k)=(qc(i,k)-qrch)/(1.+(c1d(i,k)+c0t)*dz)
4645 if(qrc(i,k).lt.0.)then ! hli new test 02/12/19
4646 qrc(i,k)=0.
4647 endif
4648 pw(i,k)=c0t*dz*qrc(i,k)*zu(i,k)
4649
4650!-----srf-08aug2017-----begin
4651! pw(i,k)=(c1d(i,k)+c0)*dz*max(0.,qrc(i,k) -qrc_crit)! units kg[rain]/kg[air]
4652!-----srf-08aug2017-----end
4653 if(qrc(i,k).lt.0)then
4654 qrc(i,k)=0.
4655 pw(i,k)=0.
4656 endif
4657 qc(i,k)=qrc(i,k)+qrch
4658 endif !autoconv
4659 pwav(i)=pwav(i)+pw(i,k)
4660 psum(i)=psum(i)+pw(i,k)*g/dp ! HCB
4661 enddo ! k=kbcon,ktop
4662! do not include liquid/ice in qc
4663!$acc loop independent
4664 do k=k22(i)+1,ktop(i)
4665!$acc atomic
4666 qc(i,k)=qc(i,k)-qrc(i,k)
4667 enddo
4668 endif ! ierr
4669!
4670!--- integrated normalized ondensate
4671!
4672 100 continue
4673!$acc end kernels
4674 prop_ave=0.
4675 iprop=0
4676!$acc parallel loop reduction(+:prop_ave,iprop)
4677 do k=kts,kte
4678 prop_ave=prop_ave+prop_b(k)
4679 if(prop_b(k).gt.0)iprop=iprop+1
4680 enddo
4681!$acc end parallel
4682 iprop=max(iprop,1)
4683
4684 end subroutine cup_up_moisture
4685
4686!--------------------------------------------------------------------
4688 real function satvap(temp2)
4689!$acc routine seq
4690 implicit none
4691 real(kind=kind_phys) :: temp2, temp, toot, toto, eilog, tsot, &
4692 & ewlog, ewlog2, ewlog3, ewlog4
4693 temp = temp2-273.155
4694 if (temp.lt.-20.) then !!!! ice saturation
4695 toot = 273.16 / temp2
4696 toto = 1 / toot
4697 eilog = -9.09718 * (toot - 1) - 3.56654 * (log(toot) / &
4698 & log(10.)) + .876793 * (1 - toto) + (log(6.1071) / log(10.))
4699 satvap = 10 ** eilog
4700 else
4701 tsot = 373.16 / temp2
4702 ewlog = -7.90298 * (tsot - 1) + 5.02808 * &
4703 & (log(tsot) / log(10.))
4704 ewlog2 = ewlog - 1.3816e-07 * &
4705 & (10 ** (11.344 * (1 - (1 / tsot))) - 1)
4706 ewlog3 = ewlog2 + .0081328 * &
4707 & (10 ** (-3.49149 * (tsot - 1)) - 1)
4708 ewlog4 = ewlog3 + (log(1013.246) / log(10.))
4709 satvap = 10 ** ewlog4
4710 end if
4711 end function
4712!--------------------------------------------------------------------
4714 subroutine get_cloud_bc(mzp,array,x_aver,k22,add)
4715!$acc routine seq
4716 implicit none
4717 integer, intent(in) :: mzp,k22
4718 real(kind=kind_phys) , dimension(:), intent(in) :: array
4719 real(kind=kind_phys) , intent(in) :: add
4720 real(kind=kind_phys) , intent(out) :: x_aver
4721 integer :: i,local_order_aver,order_aver
4722
4723 !-- dimension of the average
4724 !-- a) to pick the value at k22 level, instead of a average between
4725 !-- k22-order_aver, ..., k22-1, k22 set order_aver=1)
4726 !-- b) to average between 1 and k22 => set order_aver = k22
4727 order_aver = 3 !=> average between k22, k22-1 and k22-2
4728
4729 local_order_aver=min(k22,order_aver)
4730
4731 x_aver=0.
4732 do i = 1,local_order_aver
4733 x_aver = x_aver + array(k22-i+1)
4734 enddo
4735 x_aver = x_aver/float(local_order_aver)
4736 x_aver = x_aver + add
4737
4738 end subroutine get_cloud_bc
4739 !========================================================================================
4741 subroutine rates_up_pdf(rand_vmas,ipr,name,ktop,ierr,p_cup,entr_rate_2d,hkbo,heo,heso_cup,z_cup, &
4742 xland,kstabi,k22,kbcon,its,ite,itf,kts,kte,ktf,zuo,kpbl,ktopdby,csum,pmin_lev)
4743 implicit none
4744 character *(*), intent (in) :: name
4745 integer, intent(in) :: ipr,its,ite,itf,kts,kte,ktf
4746 real(kind=kind_phys), dimension (its:ite,kts:kte),intent (inout) :: entr_rate_2d,zuo
4747 real(kind=kind_phys), dimension (its:ite,kts:kte),intent (in) ::p_cup, heo,heso_cup,z_cup
4748 real(kind=kind_phys), dimension (its:ite),intent (in) :: hkbo,rand_vmas
4749 integer, dimension (its:ite),intent (in) :: kstabi,k22,kpbl,csum,xland,pmin_lev
4750 integer, dimension (its:ite),intent (inout) :: kbcon,ierr,ktop,ktopdby
4751!$acc declare copy(entr_rate_2d,zuo,kbcon,ierr,ktop,ktopdby) &
4752!$acc copyin(p_cup, heo,heso_cup,z_cup,hkbo,rand_vmas,kstabi,k22,kpbl,csum,xland,pmin_lev)
4753
4754 !-local vars
4755 real(kind=kind_phys), dimension (its:ite,kts:kte) :: hcot
4756!$acc declare create(hcot)
4757 real(kind=kind_phys) :: entr_init,beta_u,dz,dbythresh,dzh2,zustart,zubeg,massent,massdetr
4758 real(kind=kind_phys) :: dby(kts:kte),dbm(kts:kte),zux(kts:kte)
4759 real(kind=kind_phys) zuh2(40),zh2(40)
4760 integer :: kklev,i,kk,kbegin,k,kfinalzu
4761 integer, dimension (its:ite) :: start_level
4762!$acc declare create(start_level)
4763 logical :: is_deep, is_mid, is_shallow
4764 !
4765 zustart=.1
4766 dbythresh= 0.8 !.0.95 ! 0.85, 0.6
4767 if(name == 'shallow' .or. name == 'mid') dbythresh=1.
4768
4769 !dby(:)=0.
4770
4771 is_deep = (name .eq. 'deep')
4772 is_mid = (name .eq. 'mid')
4773 is_shallow = (name .eq. 'shallow')
4774
4775!$acc parallel loop private(beta_u,entr_init,dz,massent,massdetr,zubeg,kklev,kfinalzu,dby,dbm,zux,zuh2,zh2)
4776 do i=its,itf
4777 if(ierr(i) > 0 )cycle
4778 zux(:)=0.
4779 beta_u=max(.1,.2-float(csum(i))*.01)
4780 zuo(i,:)=0.
4781 dby(:)=0.
4782 dbm(:)=0.
4783 kbcon(i)=max(kbcon(i),2)
4784 start_level(i)=k22(i)
4785 zuo(i,start_level(i))=zustart
4786 zux(start_level(i))=zustart
4787 entr_init=entr_rate_2d(i,kts)
4788!$acc loop seq
4789 do k=start_level(i)+1,kbcon(i)
4790 dz=z_cup(i,k)-z_cup(i,k-1)
4791 massent=dz*entr_rate_2d(i,k-1)*zuo(i,k-1)
4792! massdetr=dz*1.e-9*zuo(i,k-1)
4793 massdetr=dz*.1*entr_init*zuo(i,k-1)
4794 zuo(i,k)=zuo(i,k-1)+massent-massdetr
4795 zux(k)=zuo(i,k)
4796 enddo
4797 zubeg=zustart !zuo(i,kbcon(i))
4798 if(is_deep)then
4799 ktop(i)=0
4800 hcot(i,start_level(i))=hkbo(i)
4801 dz=z_cup(i,start_level(i))-z_cup(i,start_level(i)-1)
4802!$acc loop seq
4803 do k=start_level(i)+1,ktf-2
4804 dz=z_cup(i,k)-z_cup(i,k-1)
4805
4806 hcot(i,k)=( (1.-0.5*entr_rate_2d(i,k-1)*dz)*hcot(i,k-1) &
4807 + entr_rate_2d(i,k-1)*dz*heo(i,k-1))/ &
4808 (1.+0.5*entr_rate_2d(i,k-1)*dz)
4809 if(k >= kbcon(i)) dby(k)=dby(k-1)+(hcot(i,k)-heso_cup(i,k))*dz
4810 if(k >= kbcon(i)) dbm(k)=hcot(i,k)-heso_cup(i,k)
4811 enddo
4812 ktopdby(i)=maxloc(dby(:),1)
4813 kklev=maxloc(dbm(:),1)
4814!$acc loop seq
4815 do k=maxloc(dby(:),1)+1,ktf-2
4816 if(dby(k).lt.dbythresh*maxval(dby))then
4817 kfinalzu=k - 1
4818 ktop(i)=kfinalzu
4819 go to 412
4820 endif
4821 enddo
4822 kfinalzu=ktf-2
4823 ktop(i)=kfinalzu
4824412 continue
4825 ktop(i)=ktopdby(i) ! HCB
4826 kklev=min(kklev+3,ktop(i)-2)
4827!
4828! at least overshoot by one level
4829!
4830! kfinalzu=min(max(kfinalzu,ktopdby(i)+1),ktopdby(i)+2)
4831! ktop(i)=kfinalzu
4832 if(kfinalzu.le.kbcon(i)+2)then
4833 ierr(i)=41
4834 ktop(i)= 0
4835 else
4836 call get_zu_zd_pdf_fim(kklev,p_cup(i,:),rand_vmas(i),zubeg,ipr,xland(i),zuh2,1,ierr(i),k22(i), &
4837 kfinalzu+1,zuo(i,kts:kte),kts,kte,ktf,beta_u,kbcon(i),csum(i),pmin_lev(i))
4838 endif
4839 endif ! end deep
4840 if ( is_mid ) then
4841 if(ktop(i) <= kbcon(i)+2)then
4842 ierr(i)=41
4843 ktop(i)= 0
4844 else
4845 kfinalzu=ktop(i)
4846 ktopdby(i)=ktop(i)+1
4847 call get_zu_zd_pdf_fim(kklev,p_cup(i,:),rand_vmas(i),zubeg,ipr,xland(i),zuh2,3, &
4848 ierr(i),k22(i),ktopdby(i)+1,zuo(i,kts:kte),kts,kte,ktf,beta_u,kbcon(i),csum(i),pmin_lev(i))
4849 endif
4850 endif ! mid
4851 if ( is_shallow ) then
4852 if(ktop(i) <= kbcon(i)+2)then
4853 ierr(i)=41
4854 ktop(i)= 0
4855 else
4856 kfinalzu=ktop(i)
4857 ktopdby(i)=ktop(i)+1
4858 call get_zu_zd_pdf_fim(kbcon(i),p_cup(i,:),rand_vmas(i),zubeg,ipr,xland(i),zuh2,2,ierr(i),k22(i), &
4859 ktopdby(i)+1,zuo(i,kts:kte),kts,kte,ktf,beta_u,kbcon(i),csum(i),pmin_lev(i))
4860
4861 endif
4862 endif ! shal
4863 enddo
4864!$acc end parallel loop
4865
4866 end subroutine rates_up_pdf
4867!-------------------------------------------------------------------------
4869 subroutine get_zu_zd_pdf_fim(kklev,p,rand_vmas,zubeg,ipr,xland,zuh2,draft,ierr,kb,kt,zu,kts,kte,ktf,max_mass,kpbli,csum,pmin_lev)
4870!$acc routine vector
4871
4872 implicit none
4873! real(kind=kind_phys), parameter :: beta_deep=1.3,g_beta_deep=0.8974707
4874! real(kind=kind_phys), parameter :: beta_deep=1.2,g_beta_deep=0.8974707
4875! real(kind=kind_phys), parameter :: beta_sh=2.5,g_beta_sh=1.329340
4876 real(kind=kind_phys), parameter :: beta_sh=2.2,g_beta_sh=0.8974707
4877 real(kind=kind_phys), parameter :: beta_mid=1.3,g_beta_mid=0.8974707
4878! real(kind=kind_phys), parameter :: beta_mid=1.8,g_beta_mid=0.8974707
4879 real(kind=kind_phys), parameter :: beta_dd=4.0,g_beta_dd=6.
4880 integer, intent(in) ::ipr,xland,kb,kklev,kt,kts,kte,ktf,kpbli,csum,pmin_lev
4881 real(kind=kind_phys), intent(in) ::max_mass,zubeg
4882 real(kind=kind_phys), intent(inout) :: zu(kts:kte)
4883 real(kind=kind_phys), intent(in) :: p(kts:kte)
4884 real(kind=kind_phys) :: trash,beta_deep,zuh(kts:kte),zuh2(1:40)
4885 integer, intent(inout) :: ierr
4886 integer, intent(in) ::draft
4887
4888 !- local var
4889 integer :: k1,kk,k,kb_adj,kpbli_adj,kmax
4890 real(kind=kind_phys) :: maxlim,krmax,kratio,tunning,fzu,rand_vmas,lev_start
4891 real(kind=kind_phys) :: a,b,x1,y1,g_a,g_b,alpha2,g_alpha2
4892!
4893! very simple lookup tables
4894!
4895 real(kind=kind_phys), dimension(30) :: alpha,g_alpha
4896 data (alpha(k),k=1,30)/3.699999,3.699999,3.699999,3.699999,&
4897 3.024999,2.559999,2.249999,2.028571,1.862500, &
4898 1.733333,1.630000,1.545454,1.475000,1.415385, &
4899 1.364286,1.320000,1.281250,1.247059,1.216667, &
4900 1.189474,1.165000,1.142857,1.122727,1.104348, &
4901 1.087500,1.075000,1.075000,1.075000,1.075000,1.075000/
4902 data (g_alpha(k),k=1,30)/4.170645,4.170645,4.170645,4.170645, &
4903 2.046925 , 1.387837, 1.133003, 1.012418,0.9494680, &
4904 0.9153771,0.8972442,0.8885444,0.8856795,0.8865333, &
4905 0.8897996,0.8946404,0.9005030,0.9070138,0.9139161, &
4906 0.9210315,0.9282347,0.9354376,0.9425780,0.9496124, &
4907 0.9565111,0.9619183,0.9619183,0.9619183,0.9619183,0.9619183/
4908
4909 !- kb cannot be at 1st level
4910
4911 !-- fill zu with zeros
4912 zu(:)=0.0
4913 zuh(:)=0.0
4914 kb_adj=max(kb,2)
4915
4916! Dan: replaced draft string with integer
4917! up = 1
4918! sh2 = 2
4919! mid = 3
4920! down = 4
4921! downm = 5
4922
4923 if(draft == 1) then
4924 lev_start=min(.9,.1+csum*.013)
4925 kb_adj=max(kb,2)
4926! trash is the depth of the cloud
4927 trash=-p(kt)+p(kb_adj)
4928 tunning=p(kklev)
4929 if(rand_vmas.ne.0.) tunning=p(kklev-1)+.1*rand_vmas*trash
4930 beta_deep=1.3 +(1.-trash/1200.)
4931 tunning =min(0.95, (tunning-p(kb_adj))/(p(kt)-p(kb_adj))) !=.6
4932 tunning =max(0.02, tunning)
4933 alpha2= (tunning*(beta_deep -2.)+1.)/(1.-tunning)
4934 do k=27,3,-1
4935 if(alpha(k) >= alpha2)exit
4936 enddo
4937 k1=k+1
4938 if(alpha(k1) .ne. alpha(k1-1))then
4939! write(0,*)'k1 = ',k1
4940 a=alpha(k1)-alpha(k1-1)
4941 b=alpha(k1-1)*(k1) -(k1-1)*alpha(k1)
4942 x1= (alpha2-b)/a
4943 y1=a*x1+b
4944! write(0,*)'x1,y1,a,b ',x1,y1,a,b
4945 g_a=g_alpha(k1)-g_alpha(k1-1)
4946 g_b=g_alpha(k1-1)*k1 - (k1-1)*g_alpha(k1)
4947 g_alpha2=g_a*x1+g_b
4948! write(0,*)'g_a,g_b,g_alpha2 ',g_a,g_b,g_alpha2
4949 else
4950 g_alpha2=g_alpha(k1)
4951 endif
4952
4953! fzu = gamma(alpha2 + beta_deep)/(g_alpha2*g_beta_deep)
4954 fzu = gamma(alpha2 + beta_deep)/(gamma(alpha2)*gamma(beta_deep))
4955 zu(kb_adj)=zubeg
4956 do k=kb_adj+1,min(kte,kt-1)
4957 kratio= (p(k)-p(kb_adj))/(p(kt)-p(kb_adj)) !float(k)/float(kt+1)
4958 zu(k) = zubeg+fzu*kratio**(alpha2-1.0) * (1.0-kratio)**(beta_deep-1.0)
4959 enddo
4960
4961 if(zu(kpbli).gt.0.) &
4962 zu(kts:min(ktf,kt-1))= zu(kts:min(ktf,kt-1))/zu(kpbli)
4963 do k=my_maxloc1d(zu(:),kte),1,-1
4964 if(zu(k).lt.1.e-6)then
4965 kb_adj=k+1
4966 exit
4967 endif
4968 enddo
4969 kb_adj=max(2,kb_adj)
4970 do k=kts,kb_adj-1
4971 zu(k)=0.
4972 enddo
4973 maxlim=1.2
4974 a=maxval(zu)-zu(kb_adj)
4975 do k=kb_adj,kt
4976 trash=zu(k)
4977 if(a.gt.maxlim)then
4978 zu(k)=(zu(k)-zu(kb_adj))*maxlim/a+zu(kb_adj)
4979! if(p(kt).gt.400.)write(32,122)k,p(k),zu(k),trash
4980 endif
4981 enddo
4982#ifndef _OPENACC
4983122 format(1x,i4,1x,f8.1,1x,f6.2,1x,f6.2)
4984#endif
4985 elseif(draft == 2) then
4986 k=kklev
4987 if(kpbli.gt.5)k=kpbli
4988!new nov18
4989 tunning=p(kklev) !p(kpbli+1) !p(kklev) !p(kt)+(p(kpbli)-p(kt))*lev_start
4990!end new
4991 tunning =min(0.95, (tunning-p(kb_adj))/(p(kt)-p(kb_adj))) !=.6
4992 tunning =max(0.02, tunning)
4993 alpha2= (tunning*(beta_sh -2.)+1.)/(1.-tunning)
4994 do k=27,3,-1
4995 if(alpha(k) >= alpha2)exit
4996 enddo
4997 k1=k+1
4998 if(alpha(k1) .ne. alpha(k1-1))then
4999 a=alpha(k1)-alpha(k1-1)
5000 b=alpha(k1-1)*(k1) -(k1-1)*alpha(k1)
5001 x1= (alpha2-b)/a
5002 y1=a*x1+b
5003 g_a=g_alpha(k1)-g_alpha(k1-1)
5004 g_b=g_alpha(k1-1)*k1 - (k1-1)*g_alpha(k1)
5005 g_alpha2=g_a*x1+g_b
5006 else
5007 g_alpha2=g_alpha(k1)
5008 endif
5009
5010 fzu = gamma(alpha2 + beta_sh)/(g_alpha2*g_beta_sh)
5011 zu(kb_adj) = zubeg
5012 do k=kb_adj+1,min(kte,kt-1)
5013 kratio= (p(k)-p(kb_adj))/(p(kt)-p(kb_adj)) !float(k)/float(kt+1)
5014 zu(k) = zubeg+fzu*kratio**(alpha2-1.0) * (1.0-kratio)**(beta_sh-1.0)
5015 enddo
5016
5017! beta = 2.5 !2.5 ! max(2.5,2./tunning)
5018! if(maxval(zu(kts:min(ktf,kt+1))).gt.0.) &
5019! zu(kts:min(ktf,kt+1))= zu(kts:min(ktf,kt+1))/maxval(zu(kts:min(ktf,kt+1)))
5020 if(zu(kpbli).gt.0.) &
5021 zu(kts:min(ktf,kt-1))= zu(kts:min(ktf,kt-1))/zu(kpbli)
5022 do k=my_maxloc1d(zu(:),kte),1,-1
5023 if(zu(k).lt.1.e-6)then
5024 kb_adj=k+1
5025 exit
5026 endif
5027 enddo
5028 maxlim=1.
5029 a=maxval(zu)-zu(kb_adj)
5030 do k=kts,kt
5031 if(a.gt.maxlim)zu(k)=(zu(k)-zu(kb_adj))*maxlim/a+zu(kb_adj)
5032! write(32,122)k,p(k),zu(k)
5033 enddo
5034
5035 elseif(draft == 3) then
5036 kb_adj=max(kb,2)
5037 tunning=.5*(p(kt)+p(kpbli)) !p(kt)+(p(kb_adj)-p(kt))*.9 !*.33
5038!new nov18
5039! tunning=p(kpbli) !p(kpbli+1) !p(kklev) !p(kt)+(p(kpbli)-p(kt))*lev_start
5040!end new
5041 tunning =min(0.95, (tunning-p(kb_adj))/(p(kt)-p(kb_adj))) !=.6
5042 tunning =max(0.02, tunning)
5043 alpha2= (tunning*(beta_mid -2.)+1.)/(1.-tunning)
5044 do k=27,3,-1
5045 if(alpha(k) >= alpha2)exit
5046 enddo
5047 k1=k+1
5048 if(alpha(k1) .ne. alpha(k1-1))then
5049 a=alpha(k1)-alpha(k1-1)
5050 b=alpha(k1-1)*(k1) -(k1-1)*alpha(k1)
5051 x1= (alpha2-b)/a
5052 y1=a*x1+b
5053 g_a=g_alpha(k1)-g_alpha(k1-1)
5054 g_b=g_alpha(k1-1)*k1 - (k1-1)*g_alpha(k1)
5055 g_alpha2=g_a*x1+g_b
5056 else
5057 g_alpha2=g_alpha(k1)
5058 endif
5059
5060! fzu = gamma(alpha2 + beta_deep)/(g_alpha2*g_beta_deep)
5061 fzu = gamma(alpha2 + beta_mid)/(gamma(alpha2)*gamma(beta_mid))
5062! fzu = gamma(alpha2 + beta_mid)/(g_alpha2*g_beta_mid)
5063 zu(kb_adj) = zubeg
5064 do k=kb_adj+1,min(kte,kt-1)
5065 kratio= (p(k)-p(kb_adj))/(p(kt)-p(kb_adj)) !float(k)/float(kt+1)
5066 zu(k) = zubeg+fzu*kratio**(alpha2-1.0) * (1.0-kratio)**(beta_mid-1.0)
5067 enddo
5068
5069 if(zu(kpbli).gt.0.) &
5070 zu(kts:min(ktf,kt-1))= zu(kts:min(ktf,kt-1))/zu(kpbli)
5071 do k=my_maxloc1d(zu(:),kte),1,-1
5072 if(zu(k).lt.1.e-6)then
5073 kb_adj=k+1
5074 exit
5075 endif
5076 enddo
5077 kb_adj=max(2,kb_adj)
5078 do k=kts,kb_adj-1
5079 zu(k)=0.
5080 enddo
5081 maxlim=1.5
5082 a=maxval(zu)-zu(kb_adj)
5083 do k=kts,kt
5084 if(a.gt.maxlim)zu(k)=(zu(k)-zu(kb_adj))*maxlim/a+zu(kb_adj)
5085! write(33,122)k,p(k),zu(k)
5086 enddo
5087
5088 elseif(draft == 4 .or. draft == 5) then
5089
5090 tunning=p(kb)
5091 tunning =min(0.95, (tunning-p(1))/(p(kt)-p(1))) !=.6
5092 tunning =max(0.02, tunning)
5093 alpha2= (tunning*(beta_dd -2.)+1.)/(1.-tunning)
5094 do k=27,3,-1
5095 if(alpha(k) >= alpha2)exit
5096 enddo
5097 k1=k+1
5098 if(alpha(k1) .ne. alpha(k1-1))then
5099 a=alpha(k1)-alpha(k1-1)
5100 b=alpha(k1-1)*(k1) -(k1-1)*alpha(k1)
5101 x1= (alpha2-b)/a
5102 y1=a*x1+b
5103 g_a=g_alpha(k1)-g_alpha(k1-1)
5104 g_b=g_alpha(k1-1)*k1 - (k1-1)*g_alpha(k1)
5105 g_alpha2=g_a*x1+g_b
5106 else
5107 g_alpha2=g_alpha(k1)
5108 endif
5109
5110 fzu = gamma(alpha2 + beta_dd)/(g_alpha2*g_beta_dd)
5111! fzu = gamma(alpha2 + beta_dd)/(gamma(alpha2)*gamma(beta_dd))
5112 zu(:)=0.
5113 do k=2,min(kte,kt-1)
5114 kratio= (p(k)-p(1))/(p(kt)-p(1)) !float(k)/float(kt+1)
5115 zu(k) = fzu*kratio**(alpha2-1.0) * (1.0-kratio)**(beta_dd-1.0)
5116 enddo
5117
5118 fzu=maxval(zu(kts:min(ktf,kt-1)))
5119 if(fzu.gt.0.) &
5120 zu(kts:min(ktf,kt-1))= zu(kts:min(ktf,kt-1))/fzu
5121 zu(1)=0.
5122 do k=1,kb-2 !kb,2,-1
5123 zu(kb-k)=zu(kb-k+1)-zu(kb)*(p(kb-k)-p(kb-k+1))/(p(1)-p(kb))
5124 enddo
5125 zu(1)=0.
5126 endif
5127 end subroutine get_zu_zd_pdf_fim
5128
5129!-------------------------------------------------------------------------
5131 subroutine cup_up_aa1bl(aa0,t,tn,q,qo,dtime, &
5132 z_cup,zu,dby,gamma_cup,t_cup, &
5133 kbcon,ktop,ierr, &
5134 itf,ktf, &
5135 its,ite, kts,kte )
5136
5137 implicit none
5138!
5139! on input
5140!
5141
5142 ! only local wrf dimensions are need as of now in this routine
5143
5144 integer &
5145 ,intent (in ) :: &
5146 itf,ktf, &
5147 its,ite, kts,kte
5148 ! aa0 cloud work function
5149 ! gamma_cup = gamma on model cloud levels
5150 ! t_cup = temperature (kelvin) on model cloud levels
5151 ! dby = buoancy term
5152 ! zu= normalized updraft mass flux
5153 ! z = heights of model levels
5154 ! ierr error value, maybe modified in this routine
5155 !
5156 real(kind=kind_phys), dimension (its:ite,kts:kte) &
5157 ,intent (in ) :: &
5158 z_cup,zu,gamma_cup,t_cup,dby,t,tn,q,qo
5159 integer, dimension (its:ite) &
5160 ,intent (in ) :: &
5161 kbcon,ktop
5162 real(kind=kind_phys), intent(in) :: dtime
5163!
5164! input and output
5165!
5166 integer, dimension (its:ite) &
5167 ,intent (inout) :: &
5168 ierr
5169 real(kind=kind_phys), dimension (its:ite) &
5170 ,intent (out ) :: &
5171 aa0
5172!
5173! local variables in this routine
5174!
5175 integer :: &
5176 i,k
5177 real(kind=kind_phys) :: &
5178 dz,da
5179!
5180!$acc kernels
5181 do i=its,itf
5182 aa0(i)=0.
5183 enddo
5184 do i=its,itf
5185!$acc loop independent
5186 do k=kts,kbcon(i)
5187 if(ierr(i).ne.0 ) cycle
5188! if(k.gt.kbcon(i)) cycle
5189
5190 dz = (z_cup(i,k+1)-z_cup(i,k))*g
5191 da = dz*(tn(i,k)*(1.+0.608*qo(i,k))-t(i,k)*(1.+0.608*q(i,k)))/dtime
5192!$acc atomic
5193 aa0(i)=aa0(i)+da
5194 enddo
5195 enddo
5196!$acc end kernels
5197
5198 end subroutine cup_up_aa1bl
5199!----------------------------------------------------------------------
5201 subroutine get_inversion_layers(ierr,p_cup,t_cup,z_cup,qo_cup,qeso_cup,k_inv_layers,&
5202 kstart,kend,dtempdz,itf,ktf,its,ite, kts,kte)
5203
5204 implicit none
5205 integer ,intent (in ) :: itf,ktf,its,ite,kts,kte
5206 integer, dimension (its:ite) ,intent (in ) :: ierr,kstart,kend
5207!$acc declare copyin(ierr,kstart,kend)
5208 integer, dimension (its:ite) :: kend_p3
5209!$acc declare create(kend_p3)
5210
5211 real(kind=kind_phys), dimension (its:ite,kts:kte), intent (in ) :: p_cup,t_cup,z_cup,qo_cup,qeso_cup
5212 real(kind=kind_phys), dimension (its:ite,kts:kte), intent (out) :: dtempdz
5213 integer, dimension (its:ite,kts:kte), intent (out) :: k_inv_layers
5214!$acc declare copyin(p_cup,t_cup,z_cup,qo_cup,qeso_cup)
5215!$acc declare copyout(dtempdz,k_inv_layers)
5216 !-local vars
5217 real(kind=kind_phys) :: dp,l_mid,l_shal,first_deriv(kts:kte),sec_deriv(kts:kte)
5218 integer:: ken,kadd,kj,i,k,ilev,kk,ix,k800,k550,mid,shal
5219 !
5220 !-initialize k_inv_layers as undef
5221 l_mid=300.
5222 l_shal=100.
5223!$acc kernels
5224 k_inv_layers(:,:) = 1
5225!$acc end kernels
5226!$acc parallel loop private(first_deriv,sec_deriv,ilev,ix,k,kadd,ken)
5227 do i = its,itf
5228 if(ierr(i) == 0)then
5229 sec_deriv(:)=0.
5230 kend_p3(i)=kend(i)+3
5231 do k = kts+1,kend_p3(i)+4
5232 !- get the 1st der
5233 first_deriv(k)= (t_cup(i,k+1)-t_cup(i,k-1))/(z_cup(i,k+1)-z_cup(i,k-1))
5234 dtempdz(i,k)=first_deriv(k)
5235 enddo
5236 do k = kts+2,kend_p3(i)+3
5237 ! get the 2nd der
5238 sec_deriv(k)= (first_deriv(k+1)-first_deriv(k-1))/(z_cup(i,k+1)-z_cup(i,k-1))
5239 sec_deriv(k)= abs(sec_deriv(k))
5240 enddo
5241
5242 ilev=max(kts+3,kstart(i)+1)
5243 ix=1
5244 k=ilev
5245 do while (ilev < kend_p3(i)) !(z_cup(i,ilev)<15000.)
5246!$acc loop seq
5247 do kk=k,kend_p3(i)+2 !k,ktf-2
5248
5249 if(sec_deriv(kk) < sec_deriv(kk+1) .and. &
5250 sec_deriv(kk) < sec_deriv(kk-1) ) then
5251 k_inv_layers(i,ix)=kk
5252 ix=min(5,ix+1)
5253 ilev=kk+1
5254 exit
5255 endif
5256 ilev=kk+1
5257 enddo
5258 k=ilev
5259 enddo
5260 !- 2nd criteria
5261 kadd=0
5262 ken=maxloc(k_inv_layers(i,:),1)
5263!$acc loop seq
5264 do k=1,ken
5265 kk=k_inv_layers(i,k+kadd)
5266 if(kk.eq.1)exit
5267
5268 if( dtempdz(i,kk) < dtempdz(i,kk-1) .and. &
5269 dtempdz(i,kk) < dtempdz(i,kk+1) ) then ! the layer is not a local maximum
5270 kadd=kadd+1
5271 do kj = k,ken
5272 if(k_inv_layers(i,kj+kadd).gt.1)k_inv_layers(i,kj) = k_inv_layers(i,kj+kadd)
5273 if(k_inv_layers(i,kj+kadd).eq.1)k_inv_layers(i,kj) = 1
5274 enddo
5275 endif
5276 enddo
5277 endif
5278 enddo
5279!$acc end parallel
5280100 format(1x,16i3)
5281 !- find the locations of inversions around 800 and 550 hpa
5282!$acc parallel loop private(sec_deriv,shal,mid)
5283 do i = its,itf
5284 if(ierr(i) /= 0) cycle
5285
5286 !- now find the closest layers of 800 and 550 hpa.
5287 sec_deriv(:)=1.e9
5288 do k=1,maxloc(k_inv_layers(i,:),1) !kts,kte !kstart(i),kend(i) !kts,kte
5289 dp=p_cup(i,k_inv_layers(i,k))-p_cup(i,kstart(i))
5290 sec_deriv(k)=abs(dp)-l_shal
5291 enddo
5292 k800=minloc(abs(sec_deriv),1)
5293 sec_deriv(:)=1.e9
5294
5295 do k=1,maxloc(k_inv_layers(i,:),1) !kts,kte !kstart(i),kend(i) !kts,kte
5296 dp=p_cup(i,k_inv_layers(i,k))-p_cup(i,kstart(i))
5297 sec_deriv(k)=abs(dp)-l_mid
5298 enddo
5299 k550=minloc(abs(sec_deriv),1)
5300 !-save k800 and k550 in k_inv_layers array
5301 shal=1
5302 mid=2
5303 k_inv_layers(i,shal)=k_inv_layers(i,k800) ! this is for shallow convection
5304 k_inv_layers(i,mid )=k_inv_layers(i,k550) ! this is for mid/congestus convection
5305 k_inv_layers(i,mid+1:kte)=-1
5306 enddo
5307!$acc end parallel
5308
5309 end subroutine get_inversion_layers
5310!-----------------------------------------------------------------------------------
5311! DH* 20220604 - this isn't used at all
5312!!!!>\ingroup cu_gf_deep_group
5313!!!!> This function calcualtes
5314!!! function deriv3(xx, xi, yi, ni, m)
5315!!!!$acc routine vector
5316!!! !============================================================================*/
5317!!! ! evaluate first- or second-order derivatives
5318!!! ! using three-point lagrange interpolation
5319!!! ! written by: alex godunov (october 2009)
5320!!! ! input ...
5321!!! ! xx - the abscissa at which the interpolation is to be evaluated
5322!!! ! xi() - the arrays of data abscissas
5323!!! ! yi() - the arrays of data ordinates
5324!!! ! ni - size of the arrays xi() and yi()
5325!!! ! m - order of a derivative (1 or 2)
5326!!! ! output ...
5327!!! ! deriv3 - interpolated value
5328!!! !============================================================================*/
5329!!!
5330!!! implicit none
5331!!! integer, parameter :: n=3
5332!!! integer ni, m,i, j, k, ix
5333!!! real(kind=kind_phys):: deriv3, xx
5334!!! real(kind=kind_phys):: xi(ni), yi(ni), x(n), f(n)
5335!!!
5336!!! ! exit if too high-order derivative was needed,
5337!!! if (m > 2) then
5338!!! deriv3 = 0.0
5339!!! return
5340!!! end if
5341!!!
5342!!! ! if x is ouside the xi(1)-xi(ni) interval set deriv3=0.0
5343!!! if (xx < xi(1) .or. xx > xi(ni)) then
5344!!! deriv3 = 0.0
5345!!!#ifndef _OPENACC
5346!!! stop "problems with finding the 2nd derivative"
5347!!!#else
5348!!! return
5349!!!#endif
5350!!! end if
5351!!!
5352!!! ! a binary (bisectional) search to find i so that xi(i-1) < x < xi(i)
5353!!! i = 1
5354!!! j = ni
5355!!! do while (j > i+1)
5356!!! k = (i+j)/2
5357!!! if (xx < xi(k)) then
5358!!! j = k
5359!!! else
5360!!! i = k
5361!!! end if
5362!!! end do
5363!!!
5364!!! ! shift i that will correspond to n-th order of interpolation
5365!!! ! the search point will be in the middle in x_i, x_i+1, x_i+2 ...
5366!!! i = i + 1 - n/2
5367!!!
5368!!! ! check boundaries: if i is ouside of the range [1, ... n] -> shift i
5369!!! if (i < 1) i=1
5370!!! if (i + n > ni) i=ni-n+1
5371!!!
5372!!! ! old output to test i
5373!!! ! write(*,100) xx, i
5374!!! ! 100 format (f10.5, i5)
5375!!!
5376!!! ! just wanted to use index i
5377!!! ix = i
5378!!! ! initialization of f(n) and x(n)
5379!!! do i=1,n
5380!!! f(i) = yi(ix+i-1)
5381!!! x(i) = xi(ix+i-1)
5382!!! end do
5383!!!
5384!!! ! calculate the first-order derivative using lagrange interpolation
5385!!! if (m == 1) then
5386!!! deriv3 = (2.0*xx - (x(2)+x(3)))*f(1)/((x(1)-x(2))*(x(1)-x(3)))
5387!!! deriv3 = deriv3 + (2.0*xx - (x(1)+x(3)))*f(2)/((x(2)-x(1))*(x(2)-x(3)))
5388!!! deriv3 = deriv3 + (2.0*xx - (x(1)+x(2)))*f(3)/((x(3)-x(1))*(x(3)-x(2)))
5389!!! ! calculate the second-order derivative using lagrange interpolation
5390!!! else
5391!!! deriv3 = 2.0*f(1)/((x(1)-x(2))*(x(1)-x(3)))
5392!!! deriv3 = deriv3 + 2.0*f(2)/((x(2)-x(1))*(x(2)-x(3)))
5393!!! deriv3 = deriv3 + 2.0*f(3)/((x(3)-x(1))*(x(3)-x(2)))
5394!!! end if
5395!!! end function deriv3
5396! *DH 20220604
5397!=============================================================================================
5399 subroutine get_lateral_massflux(itf,ktf, its,ite, kts,kte &
5400 ,ierr,ktop,zo_cup,zuo,cd,entr_rate_2d &
5401 ,up_massentro, up_massdetro ,up_massentr, up_massdetr &
5402 ,draft,kbcon,k22,up_massentru,up_massdetru,lambau)
5403
5404 implicit none
5405 integer, intent (in) :: draft
5406 integer, intent(in):: itf,ktf, its,ite, kts,kte
5407 integer, intent(in) , dimension(its:ite) :: ierr,ktop,kbcon,k22
5408!$acc declare copyin(ierr,ktop,kbcon,k22)
5409 !real(kind=kind_phys), intent(in), optional , dimension(its:ite):: lambau
5410 real(kind=kind_phys), intent(inout), optional , dimension(its:ite):: lambau
5411 real(kind=kind_phys), intent(in) , dimension(its:ite,kts:kte) :: zo_cup,zuo
5412 real(kind=kind_phys), intent(inout), dimension(its:ite,kts:kte) :: cd,entr_rate_2d
5413 real(kind=kind_phys), intent( out), dimension(its:ite,kts:kte) :: up_massentro, up_massdetro &
5414 ,up_massentr, up_massdetr
5415 real(kind=kind_phys), intent( out), dimension(its:ite,kts:kte), optional :: &
5416 up_massentru,up_massdetru
5417!$acc declare copy(lambau,cd,entr_rate_2d) copyin(zo_cup,zuo) copyout(up_massentro, up_massdetro,up_massentr, up_massdetr)
5418!$acc declare copyout(up_massentro, up_massdetro,up_massentr, up_massdetr, up_massentru,up_massdetru)
5419 !-- local vars
5420 integer :: i,k, incr1,incr2,turn
5421 real(kind=kind_phys) :: dz,trash,trash2
5422
5423!$acc kernels
5424 do k=kts,kte
5425 do i=its,ite
5426 up_massentro(i,k)=0.
5427 up_massdetro(i,k)=0.
5428 up_massentr(i,k)=0.
5429 up_massdetr(i,k)=0.
5430 enddo
5431 enddo
5432!$acc end kernels
5433 if(present(up_massentru) .and. present(up_massdetru))then
5434!$acc kernels
5435 do k=kts,kte
5436 do i=its,ite
5437 up_massentru(i,k)=0.
5438 up_massdetru(i,k)=0.
5439 enddo
5440 enddo
5441!$acc end kernels
5442 endif
5443!$acc parallel loop
5444 do i=its,itf
5445 if(ierr(i).eq.0)then
5446
5447!$acc loop private(dz)
5448 do k=max(2,k22(i)+1),maxloc(zuo(i,:),1)
5449 !=> below maximum value zu -> change entrainment
5450 dz=zo_cup(i,k)-zo_cup(i,k-1)
5451
5452 up_massdetro(i,k-1)=cd(i,k-1)*dz*zuo(i,k-1)
5453 up_massentro(i,k-1)=zuo(i,k)-zuo(i,k-1)+up_massdetro(i,k-1)
5454 if(up_massentro(i,k-1).lt.0.)then
5455 up_massentro(i,k-1)=0.
5456 up_massdetro(i,k-1)=zuo(i,k-1)-zuo(i,k)
5457 if(zuo(i,k-1).gt.0.)cd(i,k-1)=up_massdetro(i,k-1)/(dz*zuo(i,k-1))
5458 endif
5459 if(zuo(i,k-1).gt.0.)entr_rate_2d(i,k-1)=(up_massentro(i,k-1))/(dz*zuo(i,k-1))
5460 enddo
5461!$acc loop private(dz)
5462 do k=maxloc(zuo(i,:),1)+1,ktop(i)
5463 !=> above maximum value zu -> change detrainment
5464 dz=zo_cup(i,k)-zo_cup(i,k-1)
5465 up_massentro(i,k-1)=entr_rate_2d(i,k-1)*dz*zuo(i,k-1)
5466 up_massdetro(i,k-1)=zuo(i,k-1)+up_massentro(i,k-1)-zuo(i,k)
5467 if(up_massdetro(i,k-1).lt.0.)then
5468 up_massdetro(i,k-1)=0.
5469 up_massentro(i,k-1)=zuo(i,k)-zuo(i,k-1)
5470 if(zuo(i,k-1).gt.0.)entr_rate_2d(i,k-1)=(up_massentro(i,k-1))/(dz*zuo(i,k-1))
5471 endif
5472
5473 if(zuo(i,k-1).gt.0.)cd(i,k-1)=up_massdetro(i,k-1)/(dz*zuo(i,k-1))
5474 enddo
5475 up_massdetro(i,ktop(i))=zuo(i,ktop(i))
5476 up_massentro(i,ktop(i))=0.
5477 do k=ktop(i)+1,ktf
5478 cd(i,k)=0.
5479 entr_rate_2d(i,k)=0.
5480 up_massentro(i,k)=0.
5481 up_massdetro(i,k)=0.
5482 enddo
5483 do k=2,ktf-1
5484 up_massentr(i,k-1)=up_massentro(i,k-1)
5485 up_massdetr(i,k-1)=up_massdetro(i,k-1)
5486 enddo
5487! Dan: draft
5488! deep = 1
5489! shallow = 2
5490! mid = 3
5491 if(present(up_massentru) .and. present(up_massdetru) .and. draft == 1)then
5492 !turn=maxloc(zuo(i,:),1)
5493 !do k=2,turn
5494 ! up_massentru(i,k-1)=up_massentro(i,k-1)+.1*lambau(i)*up_massentro(i,k-1)
5495 ! up_massdetru(i,k-1)=up_massdetro(i,k-1)+.1*lambau(i)*up_massentro(i,k-1)
5496 !enddo
5497 !do k=turn+1,ktf-1
5498 do k=2,ktf-1
5499 up_massentru(i,k-1)=up_massentro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
5500 up_massdetru(i,k-1)=up_massdetro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
5501 enddo
5502 else if(present(up_massentru) .and. present(up_massdetru) .and. draft == 2)then
5503 do k=2,ktf-1
5504 up_massentru(i,k-1)=up_massentro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
5505 up_massdetru(i,k-1)=up_massdetro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
5506 enddo
5507 else if(present(up_massentru) .and. present(up_massdetru) .and. draft == 3)then
5508 lambau(i)=0.
5509 do k=2,ktf-1
5510 up_massentru(i,k-1)=up_massentro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
5511 up_massdetru(i,k-1)=up_massdetro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
5512 enddo
5513 endif
5514
5515 trash=0.
5516 trash2=0.
5517 do k=k22(i)+1,ktop(i)
5518 trash2=trash2+entr_rate_2d(i,k)
5519 enddo
5520 do k=k22(i)+1,kbcon(i)
5521 trash=trash+entr_rate_2d(i,k)
5522 enddo
5523
5524 endif
5525 enddo
5526!$acc end parallel
5527 end subroutine get_lateral_massflux
5528!---meltglac-------------------------------------------------
5529!------------------------------------------------------------------------------------
5531 subroutine get_partition_liq_ice(ierr,tn,po_cup, p_liq_ice,melting_layer &
5532 ,itf,ktf,its,ite, kts,kte, cumulus )
5533 implicit none
5534 character *(*), intent (in) :: cumulus
5535 integer ,intent (in ) :: itf,ktf, its,ite, kts,kte
5536 real(kind=kind_phys), intent (in ), dimension(its:ite,kts:kte) :: tn,po_cup
5537 real(kind=kind_phys), intent (inout), dimension(its:ite,kts:kte) :: p_liq_ice,melting_layer
5538!$acc declare copyin(tn,po_cup) copy(p_liq_ice,melting_layer)
5539 integer , intent (in ), dimension(its:ite) :: ierr
5540!$acc declare copyin(ierr)
5541 integer :: i,k
5542 real(kind=kind_phys) :: dp
5543 real(kind=kind_phys), dimension(its:ite) :: norm
5544!$acc declare create(norm)
5545 real(kind=kind_phys), parameter :: t1=276.16
5546
5547 ! hli initialize at the very beginning
5548!$acc kernels
5549 p_liq_ice(:,:) = 1.
5550 melting_layer(:,:) = 0.
5551!$acc end kernels
5552 !-- get function of t for partition of total condensate into liq and ice phases.
5553 if(melt_glac .and. cumulus == 'deep') then
5554!$acc kernels
5555 do i=its,itf
5556 if(ierr(i).eq.0)then
5557 do k=kts,ktf
5558
5559 if (tn(i,k) <= t_ice) then
5560
5561 p_liq_ice(i,k) = 0.
5562 elseif( tn(i,k) > t_ice .and. tn(i,k) < t_0) then
5563
5564 p_liq_ice(i,k) = ((tn(i,k)-t_ice)/(t_0-t_ice))**2
5565 else
5566 p_liq_ice(i,k) = 1.
5567 endif
5568
5569 !melting_layer(i,k) = p_liq_ice(i,k) * (1.-p_liq_ice(i,k))
5570 enddo
5571 endif
5572 enddo
5573 !go to 655
5574 !-- define the melting layer (the layer will be between t_0+1 < temp < t_1
5575 do i=its,itf
5576 if(ierr(i).eq.0)then
5577 do k=kts,ktf
5578 if (tn(i,k) <= t_0+1) then
5579 melting_layer(i,k) = 0.
5580
5581 elseif( tn(i,k) > t_0+1 .and. tn(i,k) < t1) then
5582 melting_layer(i,k) = ((tn(i,k)-t_0+1)/(t1-t_0+1))**2
5583
5584 else
5585 melting_layer(i,k) = 1.
5586 endif
5587 melting_layer(i,k) = melting_layer(i,k)*(1-melting_layer(i,k))
5588 enddo
5589 endif
5590 enddo
5591 655 continue
5592 !-normalize vertical integral of melting_layer to 1
5593 norm(:)=0.
5594 !do k=kts,ktf
5595 do i=its,itf
5596 if(ierr(i).eq.0)then
5597!$acc loop independent
5598 do k=kts,ktf-1
5599 dp = 100.*(po_cup(i,k)-po_cup(i,k+1))
5600!$acc atomic update
5601 norm(i) = norm(i) + melting_layer(i,k)*dp/g
5602 enddo
5603 endif
5604 enddo
5605 do i=its,itf
5606 if(ierr(i).eq.0)then
5607 !print*,"i1=",i,maxval(melting_layer(i,:)),minval(melting_layer(i,:)),norm(i)
5608 melting_layer(i,:)=melting_layer(i,:)/(norm(i)+1.e-6)*(100*(po_cup(i,kts)-po_cup(i,ktf))/g)
5609 endif
5610 !print*,"i2=",i,maxval(melting_layer(i,:)),minval(melting_layer(i,:)),norm(i)
5611 enddo
5612 !--check
5613! norm(:)=0.
5614! do k=kts,ktf-1
5615! do i=its,itf
5616! dp = 100.*(po_cup(i,k)-po_cup(i,k+1))
5617! norm(i) = norm(i) + melting_layer(i,k)*dp/g/(100*(po_cup(i,kts)-po_cup(i,ktf))/g)
5618! !print*,"n=",i,k,norm(i)
5619! enddo
5620! enddo
5621!$acc end kernels
5622 else
5623!$acc kernels
5624 p_liq_ice(:,:) = 1.
5625 melting_layer(:,:) = 0.
5626!$acc end kernels
5627 endif
5628 end subroutine get_partition_liq_ice
5629
5630!------------------------------------------------------------------------------------
5632 subroutine get_melting_profile(ierr,tn_cup,po_cup, p_liq_ice,melting_layer,qrco &
5633 ,pwo,edto,pwdo,melting &
5634 ,itf,ktf,its,ite, kts,kte, cumulus )
5635 implicit none
5636 character *(*), intent (in) :: cumulus
5637 integer ,intent (in ) :: itf,ktf, its,ite, kts,kte
5638 integer ,intent (in ), dimension(its:ite) :: ierr
5639 real(kind=kind_phys) ,intent (in ), dimension(its:ite) :: edto
5640 real(kind=kind_phys) ,intent (in ), dimension(its:ite,kts:kte) :: tn_cup,po_cup,qrco,pwo &
5641 ,pwdo,p_liq_ice,melting_layer
5642 real(kind=kind_phys) ,intent (inout), dimension(its:ite,kts:kte) :: melting
5643!$acc declare copyin(ierr,edto,tn_cup,po_cup,qrco,pwo,pwdo,p_liq_ice,melting_layer,melting)
5644 integer :: i,k
5645 real(kind=kind_phys) :: dp
5646 real(kind=kind_phys), dimension(its:ite) :: norm,total_pwo_solid_phase
5647 real(kind=kind_phys), dimension(its:ite,kts:kte) :: pwo_solid_phase,pwo_eff
5648!$acc declare create(norm,total_pwo_solid_phase,pwo_solid_phase,pwo_eff)
5649
5650 if(melt_glac .and. cumulus == 'deep') then
5651!$acc kernels
5652 !-- set melting mixing ratio to zero for columns that do not have deep convection
5653 do i=its,itf
5654 if(ierr(i) > 0) melting(i,:) = 0.
5655 enddo
5656
5657 !-- now, get it for columns where deep convection is activated
5658 total_pwo_solid_phase(:)=0.
5659
5660 !do k=kts,ktf
5661 do k=kts,ktf-1
5662 do i=its,itf
5663 if(ierr(i) /= 0) cycle
5664 dp = 100.*(po_cup(i,k)-po_cup(i,k+1))
5665
5666 !-- effective precip (after evaporation by downdraft)
5667 pwo_eff(i,k) = 0.5*(pwo(i,k)+pwo(i,k+1) + edto(i)*(pwdo(i,k)+pwdo(i,k+1)))
5668
5669 !-- precipitation at solid phase(ice/snow)
5670 pwo_solid_phase(i,k) = (1.-p_liq_ice(i,k))*pwo_eff(i,k)
5671
5672 !-- integrated precip at solid phase(ice/snow)
5673 total_pwo_solid_phase(i) = total_pwo_solid_phase(i)+pwo_solid_phase(i,k)*dp/g
5674 enddo
5675 enddo
5676
5677 do k=kts,ktf
5678 do i=its,itf
5679 if(ierr(i) /= 0) cycle
5680 !-- melting profile (kg/kg)
5681 melting(i,k) = melting_layer(i,k)*(total_pwo_solid_phase(i)/(100*(po_cup(i,kts)-po_cup(i,ktf))/g))
5682 !print*,"mel=",k,melting(i,k),pwo_solid_phase(i,k),po_cup(i,k)
5683 enddo
5684 enddo
5685
5686!-- check conservation of total solid phase precip
5687! norm(:)=0.
5688! do k=kts,ktf-1
5689! do i=its,itf
5690! dp = 100.*(po_cup(i,k)-po_cup(i,k+1))
5691! norm(i) = norm(i) + melting(i,k)*dp/g
5692! enddo
5693! enddo
5694!
5695! do i=its,itf
5696! print*,"cons=",i,norm(i),total_pwo_solid_phase(i)
5697! enddo
5698!--
5699!$acc end kernels
5700 else
5701!$acc kernels
5702 !-- no melting allowed in this run
5703 melting(:,:) = 0.
5704!$acc end kernels
5705 endif
5706 end subroutine get_melting_profile
5707!---meltglac-------------------------------------------------
5708!-----srf-08aug2017-----begin
5710 subroutine get_cloud_top(name,ktop,ierr,p_cup,entr_rate_2d,hkbo,heo,heso_cup,z_cup, &
5711 kstabi,k22,kbcon,its,ite,itf,kts,kte,ktf,zuo,kpbl,klcl,hcot)
5712 implicit none
5713 integer, intent(in) :: its,ite,itf,kts,kte,ktf
5714 real(kind=kind_phys), dimension (its:ite,kts:kte),intent (inout) :: entr_rate_2d,zuo
5715 real(kind=kind_phys), dimension (its:ite,kts:kte),intent (in) ::p_cup, heo,heso_cup,z_cup
5716 real(kind=kind_phys), dimension (its:ite),intent (in) :: hkbo
5717 integer, dimension (its:ite),intent (in) :: kstabi,k22,kbcon,kpbl,klcl
5718 integer, dimension (its:ite),intent (inout) :: ierr,ktop
5719!$acc declare copy(entr_rate_2d,zuo,ierr,ktop) copyin(p_cup, heo,heso_cup,z_cup,hkbo,kstabi,k22,kbcon,kpbl,klcl)
5720 real(kind=kind_phys), dimension (its:ite,kts:kte) :: hcot
5721!$acc declare create(hcot)
5722 character *(*), intent (in) :: name
5723 real(kind=kind_phys) :: dz,dh, dbythresh
5724 real(kind=kind_phys) :: dby(kts:kte)
5725 integer :: i,k,ipr,kdefi,kstart,kbegzu,kfinalzu
5726 integer, dimension (its:ite) :: start_level
5727!$acc declare create(start_level)
5728 integer,parameter :: find_ktop_option = 1 !0=original, 1=new
5729
5730 dbythresh=0.8 !0.95 ! the range of this parameter is 0-1, higher => lower
5731 ! overshoting (cheque aa0 calculation)
5732 ! rainfall is too sensible this parameter
5733 ! for now, keep =1.
5734 if(name == 'shallow'.or. name == 'mid')then
5735 dbythresh=1.0
5736 endif
5737 ! print*,"================================cumulus=",name; call flush(6)
5738!$acc parallel loop private(dby,kfinalzu,dz)
5739 do i=its,itf
5740 kfinalzu=ktf-2
5741 ktop(i)=kfinalzu
5742 if(ierr(i).eq.0)then
5743 dby(kts:kte)=0.0
5744
5745 start_level(i)=kbcon(i)
5746 !-- hcot below kbcon
5747 hcot(i,kts:start_level(i))=hkbo(i)
5748
5749 dz=z_cup(i,start_level(i))-z_cup(i,start_level(i)-1)
5750 dby(start_level(i))=(hcot(i,start_level(i))-heso_cup(i,start_level(i)))*dz
5751
5752 !print*,'hco1=',start_level(i),kbcon(i),hcot(i,start_level(i))/heso_cup(i,start_level(i))
5753!$acc loop seq
5754 do k=start_level(i)+1,ktf-2
5755 dz=z_cup(i,k)-z_cup(i,k-1)
5756
5757 hcot(i,k)=( (1.-0.5*entr_rate_2d(i,k-1)*dz)*hcot(i,k-1) &
5758 +entr_rate_2d(i,k-1)*dz *heo(i,k-1) )/ &
5759 (1.+0.5*entr_rate_2d(i,k-1)*dz)
5760 dby(k)=dby(k-1)+(hcot(i,k)-heso_cup(i,k))*dz
5761 !print*,'hco2=',k,hcot(i,k)/heso_cup(i,k),dby(k),entr_rate_2d(i,k-1)
5762
5763 enddo
5764 if(find_ktop_option==0) then
5765 do k=maxloc(dby(:),1),ktf-2
5766 !~ print*,'hco30=',k,dby(k),dbythresh*maxval(dby)
5767
5768 if(dby(k).lt.dbythresh*maxval(dby))then
5769 kfinalzu = k - 1
5770 ktop(i) = kfinalzu
5771 !print*,'hco4=',k,kfinalzu,ktop(i),kbcon(i)+1;call flush(6)
5772 go to 412
5773 endif
5774 enddo
5775 412 continue
5776 else
5777 do k=start_level(i)+1,ktf-2
5778 !~ print*,'hco31=',k,dby(k),dbythresh*maxval(dby)
5779
5780 if(hcot(i,k) < heso_cup(i,k) )then
5781 kfinalzu = k - 1
5782 ktop(i) = kfinalzu
5783 !print*,'hco40=',k,kfinalzu,ktop(i),kbcon(i)+1;call flush(6)
5784 exit
5785 endif
5786 enddo
5787 endif
5788 if(kfinalzu.le.kbcon(i)+1) ierr(i)=41
5789 !~ print*,'hco5=',k,kfinalzu,ktop(i),kbcon(i)+1,ierr(i);call flush(6)
5790 !
5791 endif
5792 enddo
5793!$acc end parallel
5794 end subroutine get_cloud_top
5795!------------------------------------------------------------------------------------
5797end module cu_gf_deep
cu_gf_deep::cup_dd_moisture
subroutine cup_dd_moisture(ierrc, zd, hcd, hes_cup, qcd, qes_cup, pwd, q_cup, z_cup, dd_massentr, dd_massdetr, jmin, ierr, gamma_cup, pwev, bu, qrcd, p_cup, q, he, iloop, itf, ktf, its, ite, kts, kte)
Calcultes moisture properties of downdrafts.
Definition cu_gf_deep.F90:2785
cu_gf_deep::cup_output_ens_3d
subroutine cup_output_ens_3d(xff_mid, xf_ens, ierr, dellat, dellaq, dellaqc, outtem, outq, outqc, dx, zu, pre, pw, xmb, ktop, edt, pwd, name, ierr2, ierr3, p_cup, pr_ens, maxens3, sig, closure_n, xland1, xmbm_in, xmbs_in, ichoice, imid, ipr, itf, ktf, its, ite, kts, kte, dicycle, xf_dicycle)
This subroutine calculates final output fields including physical tendencies, precipitation,...
Definition cu_gf_deep.F90:4064
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::fct1d3
subroutine fct1d3(ktop, n, dt, z, tracr, massflx, trflx_in, dellac, g)
Calculates tracer fluxes due to subsidence, only up-stream differencing is currently used but flux co...
Definition cu_gf_deep.F90:2360
cu_gf_deep::get_zu_zd_pdf_fim
subroutine get_zu_zd_pdf_fim(kklev, p, rand_vmas, zubeg, ipr, xland, zuh2, draft, ierr, kb, kt, zu, kts, kte, ktf, max_mass, kpbli, csum, pmin_lev)
Calculates a normalized mass-flux profile for updrafts and downdrafts using the beta function.
Definition cu_gf_deep.F90:4870
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::cu_gf_deep_run
subroutine cu_gf_deep_run(itf, ktf, its, ite, kts, kte, dicycle, ichoice, ipr, ccn, ccnclean, dtime, imid, kpbl, dhdt, xland, zo, forcing, t, q, z1, tn, qo, po, psur, us, vs, rho, hfx, qfx, dx, mconv, omeg, csum, cnvwt, zuo, zdo, zdm, edto, edtm, xmb_out, xmbm_in, xmbs_in, pre, outu, outv, outt, outq, outqc, kbcon, ktop, cupclw, frh_out, ierr, ierrc, nchem, fscav, chem3d, wetdpc_deep, do_smoke_transport, rand_mom, rand_vmas, rand_clos, nranflag, do_capsuppress, cap_suppress_j, k22, jmin, kdt, tropics)
Driver for the deep or congestus GF routine.
Definition cu_gf_deep.F90:147
cu_gf_deep::get_cloud_top
subroutine get_cloud_top(name, ktop, ierr, p_cup, entr_rate_2d, hkbo, heo, heso_cup, z_cup, kstabi, k22, kbcon, its, ite, itf, kts, kte, ktf, zuo, kpbl, klcl, hcot)
Calculates the cloud top height.
Definition cu_gf_deep.F90:5712
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::get_partition_liq_ice
subroutine get_partition_liq_ice(ierr, tn, po_cup, p_liq_ice, melting_layer, itf, ktf, its, ite, kts, kte, cumulus)
Calculates the partition between cloud water and cloud ice.
Definition cu_gf_deep.F90:5533
cu_gf_deep::cup_up_aa1bl
subroutine cup_up_aa1bl(aa0, t, tn, q, qo, dtime, z_cup, zu, dby, gamma_cup, t_cup, kbcon, ktop, ierr, itf, ktf, its, ite, kts, kte)
Calculates the cloud work function based on boundary layer forcing.
Definition cu_gf_deep.F90:5136
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::satvap
real function satvap(temp2)
Calculates saturation vapor pressure.
Definition cu_gf_deep.F90:4689
cu_gf_deep::cup_maximi
subroutine cup_maximi(array, ks, ke, maxx, ierr, itf, ktf, its, ite, kts, kte)
Calculates the level at which the maximum value in an array occurs.
Definition cu_gf_deep.F90:3742
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::my_maxloc1d
integer function my_maxloc1d(a, n)
Definition cu_gf_deep.F90:56
cu_gf_deep::cup_up_moisture
subroutine cup_up_moisture(name, ierr, z_cup, qc, qrc, pw, pwav, p_cup, kbcon, ktop, dby, clw_all, xland1, q, gamma_cup, zu, qes_cup, k22, qe_cup, c0, c0t3d, zqexec, ccn, ccnclean, rho, c1d, t, autoconv, up_massentr, up_massdetr, psum, psumh, itest, itf, ktf, its, ite, kts, kte)
Calculates moisture properties of the updraft.
Definition cu_gf_deep.F90:4317
cu_gf_deep::neg_check
subroutine neg_check(name, j, dt, q, outq, outt, outu, outv, outqc, pret, its, ite, kts, kte, itf, ktf, ktop)
Checks for negative or excessive tendencies and corrects in a mass conversing way by adjusting the cl...
Definition cu_gf_deep.F90:3941
cu_gf_deep::cup_forcing_ens_3d
subroutine cup_forcing_ens_3d(closure_n, xland, aa0, aa1, xaa0, mbdt, dtime, ierr, ierr2, ierr3, xf_ens, axx, forcing, maxens3, mconv, rand_clos, p_cup, ktop, omeg, zd, zdm, k22, zu, pr_ens, edt, edtm, kbcon, ichoice, imid, ipr, itf, ktf, its, ite, kts, kte, dicycle, tau_ecmwf, aa1_bl, xf_dicycle)
Calculates an ensemble of closures and the resulting ensemble average to determine cloud base mass-fl...
Definition cu_gf_deep.F90:3184
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::get_melting_profile
subroutine get_melting_profile(ierr, tn_cup, po_cup, p_liq_ice, melting_layer, qrco, pwo, edto, pwdo, melting, itf, ktf, its, ite, kts, kte, cumulus)
Calculates the melting profile.
Definition cu_gf_deep.F90:5635
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_deep::rain_evap_below_cloudbase
subroutine rain_evap_below_cloudbase(itf, ktf, its, ite, kts, kte, ierr, kbcon, xmb, psur, xland, qo_cup, po_cup, qes_cup, pwavo, edto, pwevo, pre, outt, outq)
Calculates rain evaporation below cloud base.
Definition cu_gf_deep.F90:2546
cu_gf_deep::cup_dd_edt
subroutine cup_dd_edt(ierr, us, vs, z, ktop, kbcon, edt, p, pwav, pw, ccn, ccnclean, pwev, edtmax, edtmin, edtc, psum2, psumh, rho, aeroevap, pefc, xland1, itf, ktf, its, ite, kts, kte)
Calculates strength of downdraft based on windshear and/or aerosol content.
Definition cu_gf_deep.F90:2638
cu_gf_deep
This module contains the Grell_Freitas deep convection scheme.
Definition cu_gf_deep.F90:5