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