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