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