1 MODULE module_cu_gf_deep
4 USE module_cu_gf_ctrans,only: ctrans_gf
6 real, parameter::g=9.81
7 real, parameter:: cp=1004.
8 real, parameter:: xlv=2.5e6
9 real, parameter::r_v=461.
10 real, parameter :: tcrit=258.
11 ! tuning constant for cloudwater/ice detrainment
12 real, parameter:: c1=.001 ! .0005
13 ! parameter to turn on or off evaporation of rainwater as done in SAS
14 integer, parameter :: irainevap=0
15 ! max allowed fractional coverage (frh_thresh)
16 real, parameter::frh_thresh = .9
17 ! rh threshold. if fractional coverage ~ frh_thres, do not use cupa any further
18 real, parameter::rh_thresh = .97
19 ! tuning constant for J. Brown closure (Ichoice = 4,5,6)
20 real, parameter::betajb=1.5
21 ! tuning for shallow and mid convection. EC uses 1.5
22 integer, parameter:: use_excess=1
23 real, parameter :: fluxtune=1.5
24 ! flag to turn off or modify mom transport by downdrafts
25 real, parameter :: pgcd = 1.
27 ! aerosol awareness, do not user yet!
29 integer, parameter :: autoconv=1
30 integer, parameter :: aeroevap=1
31 real, parameter :: ccnclean=250.
32 ! still 16 ensembles for clousres
33 integer, parameter:: maxens3=16
40 itf,ktf,its,ite, kts,kte &
42 ,dicycle & ! diurnal cycle flag
43 ,ichoice & ! choice of closure, use "0" for ensemble average
44 ,ipr & ! this flag can be used for debugging prints
45 ,ccn & ! not well tested yet
47 ,imid & ! flag to turn on mid level convection
49 ,kpbl & ! level of boundary layer height
50 ,dhdt & ! boundary layer forcing (one closure for shallow)
53 ,zo & ! heights above surface
54 ,forcing & ! only diagnostic
55 ,T & ! T before forcing
56 ,Q & ! Q before forcing
58 ,Tn & ! T including forcing
59 ,QO & ! Q including forcing
61 ,PSUR & ! surface pressure (mb)
62 ,US & ! u on mass points
63 ,VS & ! v on mass points
65 ,hfx & ! W/M2, positive upward
66 ,qfx & ! W/M2, positive upward
67 ,dx & ! dx is grid point dependent here
68 ,mconv & ! integrated vertical advection of moisture
69 ,omeg & ! omega (Pa/s)
71 ,csum & ! used to implement memory, set to zero if not avail
72 ,cnvwt & ! GFS needs this
73 ,zuo & ! nomalized updraft mass flux
74 ,zdo & ! nomalized downdraft mass flux
76 ,xmb_out & !the xmb's may be needed for dicycle
80 ,outu & ! momentum tendencies at mass points
82 ,outt & ! temperature tendencies
83 ,outq & ! q tendencies
84 ,outqc & ! ql/qice tendencies
87 ,cupclw & ! used for direct coupling to radiation, but with tuning factors
88 ,ierr & ! ierr flags are error flags, used for debugging
90 ! the following should be set to zero if not available
91 ,rand_mom & ! for stochastics mom, if temporal and spatial patterns exist
92 ,rand_vmas & ! for stochastics vertmass, if temporal and spatial patterns exist
93 ,rand_clos & ! for stochastics closures, if temporal and spatial patterns exist
94 ,nranflag & ! flag to what you want perturbed
95 ! 1 = momentum transport
96 ! 2 = normalized vertical mass flux profile
98 ! more is possible, talk to developer or
99 ! implement yourself. pattern is expected to be
101 #if ( WRF_DFI_RADAR == 1 )
102 ,do_capsuppress,cap_suppress_j &
104 #if ( WRF_CHEM == 1 )
105 ,num_chem,chem2d,outchemt &
106 ,num_tracer,tracer2d,outtracert &
107 ,numgas,chemopt,traceropt &
108 ,conv_tr_wetscav,conv_tr_aqchem &
118 nranflag,itf,ktf,its,ite, kts,kte,ipr,imid
119 integer, intent (in ) :: &
121 real, dimension (its:ite,4) &
122 ,intent (in ) :: rand_clos
123 real, dimension (its:ite) &
124 ,intent (in ) :: rand_mom,rand_vmas
126 #if ( WRF_DFI_RADAR == 1 )
128 ! option of cap suppress:
129 ! do_capsuppress = 1 do
130 ! do_capsuppress = other don't
133 INTEGER, INTENT(IN ) ,OPTIONAL :: do_capsuppress
134 REAL, DIMENSION( its:ite ) :: cap_suppress_j
139 real, dimension (its:ite,1:maxens3) :: xf_ens,pr_ens
140 ! outtem = output temp tendency (per s)
141 ! outq = output q tendency (per s)
142 ! outqc = output qc tendency (per s)
143 ! pre = output precip
144 real, dimension (its:ite,kts:kte) &
145 ,intent (inout ) :: &
146 cnvwt,outu,outv,OUTT,OUTQ,OUTQC,cupclw
147 real, dimension (its:ite) &
148 ,intent (inout ) :: &
150 real, dimension (its:ite) &
152 hfx,qfx,xmbm_in,xmbs_in
153 integer, dimension (its:ite) &
154 ,intent (inout ) :: &
156 integer, dimension (its:ite) &
160 ! basic environmental input includes moisture convergence (mconv)
161 ! omega (omeg), windspeed (us,vs), and a flag (ierr) to turn off
162 ! convection for this call only and at that particular gridpoint
164 real, dimension (its:ite,kts:kte) &
166 dhdt,rho,T,PO,US,VS,tn
167 real, dimension (its:ite,kts:kte) &
168 ,intent (inout ) :: &
170 real, dimension (its:ite,kts:kte) &
173 real, dimension (its:ite) &
176 real, dimension (its:ite) &
177 ,intent (inout ) :: &
185 #if ( WRF_CHEM == 1 )
186 INTEGER,INTENT(IN ) :: &
187 num_chem,num_tracer,numgas,chemopt,traceropt, &
188 conv_tr_wetscav,conv_tr_aqchem,chem_conv_tr
189 REAL,DIMENSION(its:ite , kts:kte , num_chem),INTENT(IN):: &
191 REAL,DIMENSION(its:ite , kts:kte , num_chem),INTENT(IN):: &
193 REAL,DIMENSION(its:ite , kts:kte , num_tracer),INTENT(INOUT):: &
195 REAL,DIMENSION(its:ite , kts:kte , num_tracer),INTENT(INOUT):: &
198 real,dimension(its:ite,kts:kte) :: tempco
203 ! local ensemble dependent variables in this routine
205 real, dimension (its:ite,1) :: &
207 real, dimension (its:ite,1) :: &
209 real, dimension (its:ite,kts:kte,1) :: &
210 dellat_ens,dellaqc_ens,dellaq_ens,pwo_ens
214 !***************** the following are your basic environmental
215 ! variables. They carry a "_cup" if they are
216 ! on model cloud levels (staggered). They carry
217 ! an "o"-ending (z becomes zo), if they are the forced
218 ! variables. They are preceded by x (z becomes xz)
219 ! to indicate modification by some typ of cloud
221 ! z = heights of model levels
222 ! q = environmental mixing ratio
223 ! qes = environmental saturation mixing ratio
224 ! t = environmental temp
225 ! p = environmental pressure
226 ! he = environmental moist static energy
227 ! hes = environmental saturation moist static energy
228 ! z_cup = heights of model cloud levels
229 ! q_cup = environmental q on model cloud levels
230 ! qes_cup = saturation q on model cloud levels
231 ! t_cup = temperature (Kelvin) on model cloud levels
232 ! p_cup = environmental pressure
233 ! he_cup = moist static energy on model cloud levels
234 ! hes_cup = saturation moist static energy on model cloud levels
235 ! gamma_cup = gamma on model cloud levels
238 ! hcd = moist static energy in downdraft
239 ! zd normalized downdraft mass flux
241 ! entr = entrainment rate
242 ! zd = downdraft normalized mass flux
243 ! entr= entrainment rate
244 ! hcd = h in model cloud
246 ! zd = normalized downdraft mass flux
247 ! gamma_cup = gamma on model cloud levels
248 ! qcd = cloud q (including liquid water) after entrainment
249 ! qrch = saturation q in cloud
250 ! pwd = evaporate at that level
251 ! pwev = total normalized integrated evaoprate (I2)
252 ! entr= entrainment rate
253 ! z1 = terrain elevation
254 ! entr = downdraft entrainment rate
255 ! jmin = downdraft originating level
256 ! kdet = level above ground where downdraft start detraining
257 ! psur = surface pressure
258 ! z1 = terrain elevation
259 ! pr_ens = precipitation ensemble
260 ! xf_ens = mass flux ensembles
261 ! omeg = omega from large scale model
262 ! mconv = moisture convergence from large scale model
263 ! zd = downdraft normalized mass flux
264 ! zu = updraft normalized mass flux
265 ! dir = "storm motion"
266 ! mbdt = arbitrary numerical parameter
267 ! dtime = dt over which forcing is applied
268 ! kbcon = LFC of parcel from k22
269 ! k22 = updraft originating level
270 ! ichoice = flag if only want one closure (usually set to zero!)
272 ! ktop = cloud top (output)
273 ! xmb = total base mass flux
274 ! hc = cloud moist static energy
275 ! hkb = moist static energy at originating level
277 real, dimension (its:ite,kts:kte) :: &
278 entr_rate_2d,mentrd_rate_2d,he,hes,qes,z, heo,heso,qeso,zo, &
279 xhe,xhes,xqes,xz,xt,xq,qes_cup,q_cup,he_cup,hes_cup,z_cup, &
280 p_cup,gamma_cup,t_cup, qeso_cup,qo_cup,heo_cup,heso_cup, &
281 zo_cup,po_cup,gammao_cup,tn_cup, &
282 xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup, &
283 xt_cup, dby,hc,zu,clw_all, &
284 dbyo,qco,qrcdo,pwdo,pwo,hcdo,qcdo,dbydo,hco,qrco, &
287 ! cd = detrainment function for updraft
288 ! cdd = detrainment function for downdraft
289 ! dellat = change of temperature per unit mass flux of cloud ensemble
290 ! dellaq = change of q per unit mass flux of cloud ensemble
291 ! dellaqc = change of qc per unit mass flux of cloud ensemble
293 cd,cdd,DELLAH,DELLAQ,DELLAT,DELLAQC, &
294 u_cup,v_cup,uc,vc,ucd,vcd,dellu,dellv
296 ! aa0 cloud work function for downdraft
298 ! aa0 = cloud work function without forcing effects
299 ! aa1 = cloud work function with forcing effects
300 ! xaa0 = cloud work function with cloud effects (ensemble dependent)
303 real, dimension (its:ite) :: &
304 edt,edto,AA1,AA0,XAA0,HKB, &
307 PWEVO,BU,BUD,cap_max, &
308 cap_max_increment,closure_n,psum,psumh,sig,sigd
309 real, dimension (its:ite) :: &
310 axx,edtmax,edtmin,entr_rate
311 integer, dimension (its:ite) :: &
312 kzdown,KDET,K22,JMIN,kstabi,kstabm,K22x,xland1, &
313 ktopdby,KBCONx,ierr2,ierr3,KBMAX
315 integer, dimension (its:ite), intent(inout) :: ierr
316 integer, dimension (its:ite), intent(in) :: csum
318 iloop,nens3,ki,kk,I,K
320 dz,dzo,mbdt,radius, &
321 zcutdown,depth_min,zkbmax,z_detr,zktop, &
322 dh,cap_maxs,trash,trash2,frh,sig_thresh
323 real entdo,dp,subin,detdo,entup, &
324 detup,subdown,entdoj,entupk,detupk,totmas
326 real, dimension (its:ite) :: lambau,flux_tun,zws,ztexec,zqexec
328 integer :: jprnt,jmini,start_k22
329 logical :: keep_going,flg(its:ite)
331 character*50 :: ierrc(its:ite)
332 real, dimension (its:ite,kts:kte) :: &
333 up_massentr,up_massdetr,c1d &
334 ,up_massentro,up_massdetro,dd_massentro,dd_massdetro
335 real, dimension (its:ite,kts:kte) :: &
336 up_massentru,up_massdetru,dd_massentru,dd_massdetru
337 real buo_flux,pgcon,pgc,blqe
339 real :: xff_mid(its:ite,2)
340 integer :: iversion=1
341 real :: denom,h_entr,umean,t_star,dq
342 integer, intent(IN) :: DICYCLE
343 real, dimension (its:ite) :: aa1_bl,hkbo_bl,tau_bl,tau_ecmwf,wmean
344 real, dimension (its:ite,kts:kte) :: tn_bl, qo_bl, qeso_bl, heo_bl, heso_bl &
345 ,qeso_cup_bl,qo_cup_bl, heo_cup_bl,heso_cup_bl &
346 ,gammao_cup_bl,tn_cup_bl,hco_bl,DBYo_bl
347 real, dimension(its:ite) :: xf_dicycle
348 real, intent(inout), dimension(its:ite,10) :: forcing
349 integer :: pmin_lev(its:ite),start_level(its:ite),ktopkeep(its:ite)
350 real, dimension (its:ite,kts:kte) :: dtempdz
351 integer, dimension (its:ite,kts:kte) :: k_inv_layers
355 real, dimension (its:ite) :: rntot,delqev,delq2,qevap,rn,qcond
356 real :: rain,t1,q1,elocp,evef,el2orc,evfact,evfactl,g_rain,e_dn,c_up
357 real :: pgeoh,dts,fp,fpi,pmin,x_add,beta,beta_u
358 real :: cbeg,cmid,cend,const_a,const_b,const_c
360 ! if(imid.eq.1)flux_tun(:)=fluxtune+.5
362 if(imid.eq.1)pmin=75.
365 el2orc=xlv*xlv/(r_v*cp)
368 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
370 !proportionality constant to estimate pressure gradient of updraft (Zhang and Wu, 2003, JAS
375 ! here random must be between -1 and 1
376 if(nranflag == 1)then
377 lambau(:)=1.5+rand_mom(:)
382 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
388 !- buoyancy flux (H+LE)
389 buo_flux= (hfx(i)/cp+0.608*t(i,1)*qfx(i)/xlv)/rho(i,1)
391 !-convective-scale velocity w*
392 zws(i) = max(0.,flux_tun(i)*0.41*buo_flux*zo(i,2)*g/t(i,1))
393 if(zws(i) > TINY(pgeoh)) then
394 !-convective-scale velocity w*
395 zws(i) = 1.2*zws(i)**.3333
396 !- temperature excess
397 ztexec(i) = MAX(flux_tun(i)*hfx(i)/(rho(i,1)*zws(i)*cp),0.0)
399 zqexec(i) = MAX(flux_tun(i)*qfx(i)/xlv/(rho(i,1)*zws(i)),0.)
401 !- zws for shallow convection closure (Grant 2001)
403 zws(i) = max(0.,flux_tun(i)*0.41*buo_flux*zo(i,kpbl(i))*g/t(i,kpbl(i)))
404 zws(i) = 1.2*zws(i)**.3333
405 zws(i) = zws(i)*rho(i,kpbl(i)) !check if zrho is correct
408 ! if(imid.eq.1)cap_maxs=150.
410 ! if(imid.eq.1)cap_maxs=100.
416 cap_max_increment(i)=20.
417 if(imid.eq.1)cap_max_increment(i)=10.
421 xland1(i)=int(xland(i)+.0001) ! 1.
422 if(xland(i).gt.1.5 .or. xland(i).lt.0.5)then
424 ! if(imid.eq.0)cap_max(i)=cap_maxs-25.
425 ! if(imid.eq.1)cap_max(i)=cap_maxs-50.
426 cap_max_increment(i)=20.
428 if(ztexec(i).gt.0.)cap_max(i)=cap_max(i)+25.
429 if(ztexec(i).lt.0.)cap_max(i)=cap_max(i)-25.
432 ! cap_max_increment(i)=1.
434 if(use_excess == 0 )then
438 #if ( WRF_DFI_RADAR == 1 )
439 if(do_capsuppress == 1) then
442 if (abs(cap_suppress_j(i) - 1.0 ) < 0.1 ) then
443 cap_max(i)=cap_maxs+75.
444 elseif (abs(cap_suppress_j(i) - 0.0 ) < 0.1 ) then
452 !--- initial entrainment rate (these may be changed later on in the
457 c1d(i,:)= 0. !c1 ! 0. ! c1 ! max(.003,c1+float(csum(i))*.0001)
458 entr_rate(i)=7.e-5 - min(20.,float(csum(i))) * 3.e-6
459 if(xland1(i) == 0)entr_rate(i)=7.e-5
460 if(imid.eq.1)entr_rate(i)=1.e-4
461 ! if(imid.eq.1)c1d(i,:)=c1
462 radius=.2/entr_rate(i)
463 frh=min(1.,3.14*radius*radius/dx(i)/dx(i))
464 if(frh > frh_thresh)then
466 radius=sqrt(frh*dx(i)*dx(i)/3.14)
467 entr_rate(i)=.2/radius
471 sig_thresh = (1.-frh_thresh)**2
475 !--- entrainment of mass
478 !--- initial detrainmentrates
488 cd(i,k)=1.e-9 ! 1.*entr_rate
489 ! if(imid.eq.1)cd(i,k)=entr_rate(i)
497 !--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
501 if(imid.eq.1)edtmax(:)=.15
503 if(imid.eq.1)edtmin(:)=.05
505 !--- minimum depth (m), clouds must have
508 if(imid.eq.1)depth_min=500.
510 !--- maximum depth (mb) of capping
511 !--- inversion (larger cap = no convection)
515 ! edtmax(i)=max(0.5,.8-float(csum(i))*.015) !.3)
516 ! if(xland1(i) == 1 )edtmax(i)=max(0.7,1.-float(csum(i))*.015) !.3)
528 ! cap_max(i)=cap_maxs
533 !--- max height(m) above ground where updraft air can originate
536 if(imid.eq.1)zkbmax=2000.
538 !--- height(m) above which no downdrafts are allowed to originate
542 !--- depth(m) over which downdraft detrains all its mass
545 ! if(imid.eq.1)z_detr=800.
549 !--- environmental conditions, FIRST HEIGHTS
559 !--- calculate moist static energy, heights, qes
561 call cup_env(z,qes,he,hes,t,q,po,z1, &
562 psur,ierr,tcrit,-1, &
565 call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
566 psur,ierr,tcrit,-1, &
571 !--- environmental values on cloud levels
573 call cup_env_clev(t,qes,q,he,hes,z,po,qes_cup,q_cup,he_cup, &
574 hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
578 call cup_env_clev(tn,qeso,qo,heo,heso,zo,po,qeso_cup,qo_cup, &
579 heo_cup,heso_cup,zo_cup,po_cup,gammao_cup,tn_cup,psur, &
585 if(kpbl(i).gt.5 .and. imid.eq.1)cap_max(i)=po_cup(i,kpbl(i))
586 u_cup(i,kts)=us(i,kts)
587 v_cup(i,kts)=vs(i,kts)
589 u_cup(i,k)=.5*(us(i,k-1)+us(i,k))
590 v_cup(i,k)=.5*(vs(i,k-1)+vs(i,k))
597 if(zo_cup(i,k).gt.zkbmax+z1(i))then
604 !--- level where detrainment for downdraft starts
607 if(zo_cup(i,k).gt.z_detr+z1(i))then
619 !------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
624 k22(i)=maxloc(HEO_CUP(i,start_k22:kbmax(i)+2),1)+start_k22-1
625 if(K22(I).GE.KBMAX(i))then
627 ierrc(i)="could not find k22"
635 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
640 x_add = xlv*zqexec(i)+cp*ztexec(i)
641 call get_cloud_bc(kte,he_cup (i,1:kte),hkb (i),k22(i),x_add)
642 call get_cloud_bc(kte,heo_cup (i,1:kte),hkbo (i),k22(i),x_add)
648 call cup_kbcon(ierrc,cap_max_increment,iloop,k22,kbcon,heo_cup,heso_cup, &
649 hkbo,ierr,kbmax,po_cup,cap_max, &
653 z_cup,entr_rate,heo,imid)
655 !--- increase detrainment in stable layers
657 CALL cup_minimi(HEso_cup,Kbcon,kstabm,kstabi,ierr, &
662 frh = min(qo_cup(i,kbcon(i))/qeso_cup(i,kbcon(i)),1.)
663 if(frh >= rh_thresh .and. sig(i) <= sig_thresh )then
668 ! never go too low...
670 ! if(imid.eq.0 .and. xland1(i).eq.0)x_add=150.
673 if(po(i,kbcon(i))-po(i,k) > pmin+x_add)then
679 ! initial conditions for updraft
681 start_level(i)=k22(i)
682 x_add = xlv*zqexec(i)+cp*ztexec(i)
683 call get_cloud_bc(kte,he_cup (i,1:kte),hkb (i),k22(i),x_add)
687 !--- get inversion layers for mid level cloud tops
690 call get_inversion_layers(ierr,p_cup,t_cup,z_cup,q_cup,qes_cup,k_inv_layers, &
691 kbcon,kstabi,dtempdz,itf,ktf,its,ite, kts,kte)
694 if(kstabi(i).lt.kbcon(i))then
699 entr_rate_2d(i,k)=entr_rate(i)
702 ! if(imid.eq.0 .and. pmin_lev(i).lt.kbcon(i)+3)pmin_lev(i)=kbcon(i)+3
703 kbcon(i)=max(2,kbcon(i))
705 frh = min(qo_cup(i,k)/qeso_cup(i,k),1.)
706 entr_rate_2d(i,k)=entr_rate(i) *(1.3-frh)
709 if(k_inv_layers(i,2).gt.0 .and. &
710 (po_cup(i,k22(i))-po_cup(i,k_inv_layers(i,2))).lt.500.)then
712 ktop(i)=min(kstabi(i),k_inv_layers(i,2))
716 if((po_cup(i,k22(i))-po_cup(i,k)).gt.500.)then
728 !-- get normalized mass flux, entrainment and detrainmentrates for updraft
731 !- for mid level clouds we do not allow clouds taller than where stability
734 call rates_up_pdf(rand_vmas,ipr,'mid',ktop,ierr,po_cup,entr_rate_2d,hkbo,heo,heso_cup,zo_cup, &
735 xland1,kstabi,k22,kbcon,its,ite,itf,kts,kte,ktf,zuo,kpbl,ktopdby,csum,pmin_lev)
737 call rates_up_pdf(rand_vmas,ipr,'deep',ktop,ierr,po_cup,entr_rate_2d,hkbo,heo,heso_cup,zo_cup, &
738 xland1,kstabi,k22,kbcon,its,ite,itf,kts,kte,ktf,zuo,kbcon,ktopdby,csum,pmin_lev)
765 ! calculate mass entrainment and detrainment
767 CALL get_lateral_massflux(itf,ktf, its,ite, kts,kte &
768 ,ierr,ktop,zo_cup,zuo,cd,entr_rate_2d &
769 ,up_massentro, up_massdetro ,up_massentr, up_massdetr &
770 ,'deep',kbcon,k22,up_massentru,up_massdetru,lambau)
774 ! NOTE: Ktop here already includes overshooting, ktopdby is without
789 do k=1,start_level(i)
793 do k=1,start_level(i)-1
795 hco(i,k)=heo_cup(i,k)
807 if(ierr(i) /= 0) cycle
809 DO k=start_level(i) +1,ktop(i) !mass cons option
811 denom=zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1)
812 if(denom.lt.1.e-8)then
817 hc(i,k)=(hc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*hc(i,k-1)+ &
818 up_massentr(i,k-1)*he(i,k-1)) / &
819 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
820 uc(i,k)=(uc(i,k-1)*zu(i,k-1)-.5*up_massdetru(i,k-1)*uc(i,k-1)+ &
821 up_massentru(i,k-1)*us(i,k-1) &
822 -pgcon*.5*(zu(i,k)+zu(i,k-1))*(u_cup(i,k)-u_cup(i,k-1))) / &
823 (zu(i,k-1)-.5*up_massdetru(i,k-1)+up_massentru(i,k-1))
824 vc(i,k)=(vc(i,k-1)*zu(i,k-1)-.5*up_massdetru(i,k-1)*vc(i,k-1)+ &
825 up_massentru(i,k-1)*vs(i,k-1) &
826 -pgcon*.5*(zu(i,k)+zu(i,k-1))*(v_cup(i,k)-v_cup(i,k-1))) / &
827 (zu(i,k-1)-.5*up_massdetru(i,k-1)+up_massentru(i,k-1))
828 dby(i,k)=hc(i,k)-hes_cup(i,k)
829 hco(i,k)=(hco(i,k-1)*zuo(i,k-1)-.5*up_massdetro(i,k-1)*hco(i,k-1)+ &
830 up_massentro(i,k-1)*heo(i,k-1)) / &
831 (zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
832 dbyo(i,k)=hco(i,k)-heso_cup(i,k)
833 DZ=Zo_cup(i,K+1)-Zo_cup(i,K)
834 dbyt(i,k)=dbyt(i,k-1)+dbyo(i,k)*dz
836 ! for now no overshooting (only very little)
837 kk=maxloc(dbyt(i,:),1)
838 ki=maxloc(zuo(i,:),1)
839 ! if(ipr .eq.1)write(16,*)'cupgf2',kk,ki
844 do k=ktop(i)-1,kbcon(i),-1
845 if(dbyo(i,k).gt.0.)then
851 if(ierr(i).eq.0)ktop(I)=ktopkeep(i)
855 if(ierr(i) /= 0) cycle
860 HCo(i,K)=heso_cup(i,k)
876 if(ktop(i).lt.kbcon(i)+2)then
878 ierrc(i)='ktop too small deep'
885 zktop=(zo_cup(i,ktop(i))-z1(i))*.6
886 if(imid.eq.1)zktop=(zo_cup(i,ktop(i))-z1(i))*.4
887 zktop=min(zktop+z1(i),zcutdown+z1(i))
889 if(zo_cup(i,k).gt.zktop)then
891 kzdown(i)=min(kzdown(i),kstabi(i)-1) !
898 !--- DOWNDRAFT ORIGINATING LEVEL - JMIN
900 call cup_minimi(HEso_cup,K22,kzdown,JMIN,ierr, &
906 !--- check whether it would have buoyancy, if there where
907 !--- no entrainment/detrainment
911 do while ( keep_going )
913 if ( jmini - 1 .lt. kdet(i) ) kdet(i) = jmini-1
914 if ( jmini .ge. ktop(i)-1 ) jmini = ktop(i) - 2
916 hcdo(i,ki)=heso_cup(i,ki)
917 DZ=Zo_cup(i,Ki+1)-Zo_cup(i,Ki)
920 hcdo(i,k)=heso_cup(i,jmini)
921 DZ=Zo_cup(i,K+1)-Zo_cup(i,K)
922 dh=dh+dz*(HCDo(i,K)-heso_cup(i,k))
925 if ( jmini .gt. 5 ) then
929 ierrc(i) = "could not find jmini9"
936 if ( jmini .le. 5 ) then
938 ierrc(i) = "could not find jmini4"
943 ! - Must have at least depth_min m between cloud convective base
948 if ( jmin(i) - 1 .lt. kdet(i) ) kdet(i) = jmin(i)-1
949 IF(-zo_cup(I,KBCON(I))+zo_cup(I,KTOP(I)).LT.depth_min)then
951 ierrc(i)="cloud depth very shallow"
957 !--- normalized downdraft mass flux profile,also work on bottom detrainment
968 hcdo(i,k)=heso_cup(i,k)
972 mentrd_rate_2d(i,k)=entr_rate(i)
976 beta=max(.02,.05-float(csum(i))*.0015) !.02
977 ! beta=max(.05,.08-float(csum(i))*.0015) !.02
978 if(imid.eq.0 .and. xland1(i) == 0)then
980 edtmax(i)=max(0.1,.4-float(csum(i))*.015) !.3)
982 if(imid.eq.1)beta=.02
985 cdd(i,1:jmin(i))=1.e-9
989 call get_zu_zd_pdf_fim(0,po_cup(i,:),rand_vmas(i),0.,ipr,xland1(i),zuh2,"DOWN",ierr(i),kdet(i),jmin(i),zdo(i,:),kts,kte,ktf,beta,kpbl(i),csum(i),pmin_lev(i))
990 if(zdo(i,jmin(i)) .lt.1.e-8)then
993 if(zdo(i,jmin(i)) .lt.1.e-8)then
999 do ki=jmin(i) ,maxloc(zdo(i,:),1),-1
1000 !=> from jmin to maximum value zd -> change entrainment
1001 dzo=zo_cup(i,ki+1)-zo_cup(i,ki)
1002 dd_massdetro(i,ki)=cdd(i,ki)*dzo*zdo(i,ki+1)
1003 dd_massentro(i,ki)=zdo(i,ki)-zdo(i,ki+1)+dd_massdetro(i,ki)
1004 if(dd_massentro(i,ki).lt.0.)then
1005 dd_massentro(i,ki)=0.
1006 dd_massdetro(i,ki)=zdo(i,ki+1)-zdo(i,ki)
1007 if(zdo(i,ki+1).gt.0.)cdd(i,ki)=dd_massdetro(i,ki)/(dzo*zdo(i,ki+1))
1009 if(zdo(i,ki+1).gt.0.)mentrd_rate_2d(i,ki)=dd_massentro(i,ki)/(dzo*zdo(i,ki+1))
1011 mentrd_rate_2d(i,1)=0.
1012 do ki=maxloc(zdo(i,:),1)-1,1,-1
1013 !=> from maximum value zd to surface -> change detrainment
1014 dzo=zo_cup(i,ki+1)-zo_cup(i,ki)
1015 dd_massentro(i,ki)=mentrd_rate_2d(i,ki)*dzo*zdo(i,ki+1)
1016 dd_massdetro(i,ki) = zdo(i,ki+1)+dd_massentro(i,ki)-zdo(i,ki)
1017 if(dd_massdetro(i,ki).lt.0.)then
1018 dd_massdetro(i,ki)=0.
1019 dd_massentro(i,ki)=zdo(i,ki)-zdo(i,ki+1)
1020 if(zdo(i,ki+1).gt.0.)mentrd_rate_2d(i,ki)=dd_massentro(i,ki)/(dzo*zdo(i,ki+1))
1022 if(zdo(i,ki+1).gt.0.)cdd(i,ki)= dd_massdetro(i,ki)/(dzo*zdo(i,ki+1))
1024 cbeg=po_cup(i,kbcon(i)) !850.
1025 cend=min(po_cup(i,ktop(i)),400.)
1026 cmid=.5*(cbeg+cend) !600.
1027 const_b=c1/((cmid*cmid-cbeg*cbeg)*(cbeg-cend)/(cend*cend-cbeg*cbeg)+cmid-cbeg)
1028 const_a=const_b*(cbeg-cend)/(cend*cend-cbeg*cbeg)
1029 const_c=-const_a*cbeg*cbeg-const_b*cbeg
1030 do k=kbcon(i)+1,ktop(i)-1
1031 c1d(i,k)=const_a*po_cup(i,k)*po_cup(i,k)+const_b*po_cup(i,k)+const_c
1032 c1d(i,k)=max(0.,c1d(i,k))
1035 if(imid.eq.1)c1d(i,:)=0.
1039 ! c1d(i,jmin(i)-2)=c1/40.
1040 ! if(imid.eq.1)c1d(i,jmin(i)-2)=c1/20.
1041 ! do k=jmin(i)-1,ktop(i)
1042 ! dz=zo_cup(i,ktop(i))-zo_cup(i,jmin(i))
1043 ! c1d(i,k)=c1d(i,k-1)+c1*(zo_cup(i,k+1)-zo_cup(i,k))/dz
1044 ! c1d(i,k)=max(0.,c1d(i,k))
1045 ! c1d(i,k)=min(.002,c1d(i,k))
1049 ! downdraft moist static energy + moisture budget
1051 dd_massentru(i,k-1)=dd_massentro(i,k-1)+lambau(i)*dd_massdetro(i,k-1)
1052 dd_massdetru(i,k-1)=dd_massdetro(i,k-1)+lambau(i)*dd_massdetro(i,k-1)
1054 dbydo(i,jmin(i))=hcdo(i,jmin(i))-heso_cup(i,jmin(i))
1055 bud(i)=dbydo(i,jmin(i))*(zo_cup(i,jmin(i)+1)-zo_cup(i,jmin(i)))
1057 dzo=zo_cup(i,ki+1)-zo_cup(i,ki)
1058 h_entr=.5*(heo(i,ki)+.5*(hco(i,ki)+hco(i,ki+1)))
1059 ucd(i,ki)=(ucd(i,ki+1)*zdo(i,ki+1) &
1060 -.5*dd_massdetru(i,ki)*ucd(i,ki+1)+ &
1061 dd_massentru(i,ki)*us(i,ki) &
1062 -pgcon*zdo(i,ki+1)*(us(i,ki+1)-us(i,ki))) / &
1063 (zdo(i,ki+1)-.5*dd_massdetru(i,ki)+dd_massentru(i,ki))
1064 vcd(i,ki)=(vcd(i,ki+1)*zdo(i,ki+1) &
1065 -.5*dd_massdetru(i,ki)*vcd(i,ki+1)+ &
1066 dd_massentru(i,ki)*vs(i,ki) &
1067 -pgcon*zdo(i,ki+1)*(vs(i,ki+1)-vs(i,ki))) / &
1068 (zdo(i,ki+1)-.5*dd_massdetru(i,ki)+dd_massentru(i,ki))
1069 hcdo(i,ki)=(hcdo(i,ki+1)*zdo(i,ki+1) &
1070 -.5*dd_massdetro(i,ki)*hcdo(i,ki+1)+ &
1071 dd_massentro(i,ki)*h_entr) / &
1072 (zdo(i,ki+1)-.5*dd_massdetro(i,ki)+dd_massentro(i,ki))
1073 dbydo(i,ki)=hcdo(i,ki)-heso_cup(i,ki)
1074 bud(i)=bud(i)+dbydo(i,ki)*dzo
1080 ierrc(i)='downdraft is not negatively buoyant '
1084 !--- calculate moisture properties of downdraft
1086 call cup_dd_moisture(ierrc,zdo,hcdo,heso_cup,qcdo,qeso_cup, &
1087 pwdo,qo_cup,zo_cup,dd_massentro,dd_massdetro,jmin,ierr,gammao_cup, &
1088 pwevo,bu,qrcdo,qo,heo,1, &
1092 !--- calculate moisture properties of updraft
1095 call cup_up_moisture('mid',ierr,zo_cup,qco,qrco,pwo,pwavo, &
1096 p_cup,kbcon,ktop,dbyo,clw_all,xland1, &
1097 qo,GAMMAo_cup,zuo,qeso_cup,k22,qo_cup, &
1098 ZQEXEC,ccn,rho,c1d,tn_cup,up_massentr,up_massdetr,psum,psumh, &
1102 call cup_up_moisture('deep',ierr,zo_cup,qco,qrco,pwo,pwavo, &
1103 p_cup,kbcon,ktop,dbyo,clw_all,xland1, &
1104 qo,GAMMAo_cup,zuo,qeso_cup,k22,qo_cup, &
1105 ZQEXEC,ccn,rho,c1d,tn_cup,up_massentr,up_massdetr,psum,psumh, &
1110 if(ierr(i).eq.0)then
1112 dp=100.*(po_cup(i,1)-po_cup(i,2))
1113 cupclw(i,k)=qrco(i,k) ! my mod
1114 cnvwt(i,k)=zuo(i,k)*cupclw(i,k)*g/dp
1119 !--- calculate workfunctions for updrafts
1121 call cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
1125 call cup_up_aa0(aa1,zo,zuo,dbyo,GAMMAo_CUP,tn_cup, &
1130 if(ierr(i).eq.0)then
1131 if(aa1(i).eq.0.)then
1133 ierrc(i)="cloud work function zero"
1138 !--- diurnal cycle closure
1140 !--- AA1 from boundary layer (bl) processes only
1142 xf_dicycle (:) = 0.0
1144 !- way to calculate the fraction of cape consumed by shallow convection
1148 ! Betchold et al 2008 time-scale of cape removal
1150 ! wmean is of no meaning over land....
1151 ! still working on replacing it over water
1154 if(ierr(i).eq.0)then
1155 !- mean vertical velocity
1156 wmean(i) = 7.0 ! m/s ! in the future change for Wmean == integral( W dz) / cloud_depth
1157 if(imid.eq.1)wmean(i) = 3.0
1158 !- time-scale cape removal from Betchold et al. 2008
1159 tau_ecmwf(i)=( zo_cup(i,ktopdby(i))- zo_cup(i,kbcon(i)) ) / wmean(i)
1160 tau_ecmwf(i)= tau_ecmwf(i) * (1.0061 + 1.23E-2 * (dx(i)/1000.))! dx(i) must be in meters
1165 IF(dicycle == 1) then
1168 if(ierr(i).eq.0)then
1169 if(xland1(i) == 0 ) then
1171 umean= 2.0+sqrt(2.0*(US(i,1)**2+VS(i,1)**2+US(i,kbcon(i))**2+VS(i,kbcon(i))**2))
1172 tau_bl(i) = (zo_cup(i,kbcon(i))- z1(i)) /umean
1175 tau_bl(i) =( zo_cup(i,ktopdby(i))- zo_cup(i,kbcon(i)) ) / wmean(i)
1181 if(iversion == 1) then
1183 t_star=4. !original =1
1185 !-- calculate pcape from BL forcing only
1186 call cup_up_aa1bl(aa1_bl,t,tn,q,qo,dtime, &
1187 zo_cup,zuo,dbyo_bl,GAMMAo_CUP_bl,tn_cup_bl, &
1189 itf,ktf,its,ite, kts,kte)
1193 if(ierr(i).eq.0)then
1195 !- only for convection rooting in the PBL
1196 if(zo_cup(i,kbcon(i))-z1(i) > zo(i,min(kte,kpbl(i)+1))) then
1199 !- multiply aa1_bl the " time-scale" - tau_bl
1200 aa1_bl(i) = max(0.,aa1_bl(i)/t_star* tau_bl(i))
1207 !- version for real cloud-work function
1209 !-get the profiles modified only by bl tendencies
1211 tn_bl(i,:)=0.;qo_bl(i,:)=0.
1212 if ( ierr(i) == 0 )then
1213 !below kbcon -> modify profiles
1214 tn_bl(i,1:kbcon(i)) = tn(i,1:kbcon(i))
1215 qo_bl(i,1:kbcon(i)) = qo(i,1:kbcon(i))
1216 !above kbcon -> keep environment profiles
1217 tn_bl(i,kbcon(i)+1:ktf) = t(i,kbcon(i)+1:ktf)
1218 qo_bl(i,kbcon(i)+1:ktf) = q(i,kbcon(i)+1:ktf)
1221 !--- calculate moist static energy, heights, qes, ... only by bl tendencies
1222 call cup_env(zo,qeso_bl,heo_bl,heso_bl,tn_bl,qo_bl,po,z1, &
1223 psur,ierr,tcrit,-1, &
1224 itf,ktf, its,ite, kts,kte)
1225 !--- environmental values on cloud levels only by bl tendencies
1226 call cup_env_clev(tn_bl,qeso_bl,qo_bl,heo_bl,heso_bl,zo,po,qeso_cup_bl,qo_cup_bl, &
1227 heo_cup_bl,heso_cup_bl,zo_cup,po_cup,gammao_cup_bl,tn_cup_bl,psur, &
1229 itf,ktf,its,ite, kts,kte)
1231 IF(ierr(I).eq.0)THEN
1232 hkbo_bl(i)=heo_cup_bl(i,k22(i))
1242 IF(ierr(I).eq.0)THEN
1244 hco_bl(i,k)=hkbo_bl(i)
1247 hco_bl (i,k)=hkbo_bl(i)
1248 DBYo_bl(i,k)=Hkbo_bl(i) - HESo_cup_bl(i,k)
1254 if(ierr(i).eq.0)then
1255 do k=kbcon(i)+1,ktop(i)
1256 hco_bl(i,k)=(hco_bl(i,k-1)*zuo(i,k-1)-.5*up_massdetro(i,k-1)*hco_bl(i,k-1)+ &
1257 up_massentro(i,k-1)*heo_bl(i,k-1)) / &
1258 (zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
1259 dbyo_bl(i,k)=hco_bl(i,k)-heso_cup_bl(i,k)
1262 hco_bl (i,k)=heso_cup_bl(i,k)
1268 !--- calculate workfunctions for updrafts
1269 call cup_up_aa0(aa1_bl,zo,zuo,dbyo_bl,GAMMAo_CUP_bl,tn_cup_bl, &
1271 itf,ktf,its,ite, kts,kte)
1275 if(ierr(i).eq.0)then
1276 !- get the increment on AA0 due the BL processes
1277 aa1_bl(i) = aa1_bl(i) - aa0(i)
1278 !- only for convection rooting in the PBL
1279 !if(zo_cup(i,kbcon(i))-z1(i) > 500.0) then !- instead 500 -> zo_cup(kpbl(i))
1282 ! !- multiply aa1_bl the "normalized time-scale" - tau_bl/ model_timestep
1283 aa1_bl(i) = aa1_bl(i)* tau_bl(i)/ dtime
1285 print*,'aa0,aa1bl=',aa0(i),aa1_bl(i),aa0(i)-aa1_bl(i),tau_bl(i)!,dtime,xland(i)
1289 ENDIF ! version of implementation
1295 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
1297 call cup_dd_edt(ierr,us,vs,zo,ktop,kbcon,edt,po,pwavo, &
1298 pwo,ccn,pwevo,edtmax,edtmin,edtc,psum,psumh, &
1299 rho,aeroevap,itf,ktf, &
1302 if(ierr(i).eq.0)then
1308 dellat_ens (i,k,1)=0.
1309 dellaq_ens (i,k,1)=0.
1310 dellaqc_ens(i,k,1)=0.
1315 !--- change per unit mass that a model cloud would modify the environment
1317 !--- 1. in bottom layer
1330 !---------------------------------------------- cloud level ktop
1332 !- - - - - - - - - - - - - - - - - - - - - - - - model level ktop-1
1340 !---------------------------------------------- cloud level k+2
1342 !- - - - - - - - - - - - - - - - - - - - - - - - model level k+1
1344 !---------------------------------------------- cloud level k+1
1346 !- - - - - - - - - - - - - - - - - - - - - - - - model level k
1348 !---------------------------------------------- cloud level k
1361 !---------------------------------------------- cloud level 3
1363 !- - - - - - - - - - - - - - - - - - - - - - - - model level 2
1365 !---------------------------------------------- cloud level 2
1367 !- - - - - - - - - - - - - - - - - - - - - - - - model level 1
1370 if(ierr(i).eq.0)then
1371 dp=100.*(po_cup(i,1)-po_cup(i,2))
1372 dellu(i,1)=pgcd*(edto(i)*zdo(i,2)*ucd(i,2) &
1373 -edto(i)*zdo(i,2)*u_cup(i,2))*g/dp
1374 dellv(i,1)=pgcd*(edto(i)*zdo(i,2)*vcd(i,2) &
1375 -edto(i)*zdo(i,2)*v_cup(i,2))*g/dp
1378 ! these three are only used at or near mass detrainment and/or entrainment levels
1381 if(k == k22(i)-1) entupk=zuo(i,k+1)
1384 ! detrainment and entrainment for fowndrafts
1385 detdo=edto(i)*dd_massdetro(i,k)
1386 entdo=edto(i)*dd_massentro(i,k)
1387 ! entrainment/detrainment for updraft
1388 entup=up_massentro(i,k)
1389 detup=up_massdetro(i,k)
1390 ! subsidence by downdrafts only
1391 subin=-zdo(i,k+1)*edto(i)
1392 subdown=-zdo(i,k)*edto(i)
1394 if(k.eq.ktop(i))then
1395 detupk=zuo(i,ktop(i))
1403 totmas=subin-subdown+detup-entup-entdo+ &
1404 detdo-entupk-entdoj+detupk+zuo(i,k+1)-zuo(i,k)
1405 if(abs(totmas).gt.1.e-6)then
1406 write(0,123)'totmas=',k22(i),kbcon(i),k,entup,detup,zuo(i,k+1),zuo(i,k),detdo,entdo
1407 123 formAT(a7,1X,3i3,2E12.4,2(1x,f5.2),2e12.4)
1409 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
1411 if(k.ge.ktop(i))pgc=0.
1413 dellu(i,k) =-(zuo(i,k+1)*(uc (i,k+1)-u_cup(i,k+1) ) - &
1414 zuo(i,k )*(uc (i,k )-u_cup(i,k ) ) )*g/dp &
1415 +(zdo(i,k+1)*(ucd(i,k+1)-u_cup(i,k+1) ) - &
1416 zdo(i,k )*(ucd(i,k )-u_cup(i,k ) ) )*g/dp*edto(i)*pgcd
1417 dellv(i,k) =-(zuo(i,k+1)*(vc (i,k+1)-v_cup(i,k+1) ) - &
1418 zuo(i,k )*(vc (i,k )-v_cup(i,k ) ) )*g/dp &
1419 +(zdo(i,k+1)*(vcd(i,k+1)-v_cup(i,k+1) ) - &
1420 zdo(i,k )*(vcd(i,k )-v_cup(i,k ) ) )*g/dp*edto(i)*pgcd
1431 if(ierr(i).eq.0)then
1433 dp=100.*(po_cup(i,1)-po_cup(i,2))
1435 dellah(i,1)=(edto(i)*zdo(i,2)*hcdo(i,2) &
1436 -edto(i)*zdo(i,2)*heo_cup(i,2))*g/dp
1438 dellaq (i,1)=(edto(i)*zdo(i,2)*qcdo(i,2) &
1439 -edto(i)*zdo(i,2)*qo_cup(i,2))*g/dp
1441 G_rain= 0.5*(pwo (i,1)+pwo (i,2))*g/dp
1442 E_dn = -0.5*(pwdo(i,1)+pwdo(i,2))*g/dp*edto(i) ! pwdo < 0 and E_dn must > 0
1443 dellaq(i,1) = dellaq(i,1)+ E_dn-G_rain
1445 !--- conservation check
1446 !- water mass balance
1447 !trash = trash + (dellaq(i,1)+dellaqc(i,1)+G_rain-E_dn)*dp/g
1449 !trash2 = trash2+ (dellah(i,1))*dp/g
1453 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
1454 ! these three are only used at or near mass detrainment and/or entrainment levels
1456 dellah(i,k) =-(zuo(i,k+1)*(hco (i,k+1)-heo_cup(i,k+1) ) - &
1457 zuo(i,k )*(hco (i,k )-heo_cup(i,k ) ) )*g/dp &
1458 +(zdo(i,k+1)*(hcdo(i,k+1)-heo_cup(i,k+1) ) - &
1459 zdo(i,k )*(hcdo(i,k )-heo_cup(i,k ) ) )*g/dp*edto(i)
1462 !- check H conservation
1463 ! trash2 = trash2+ (dellah(i,k))*dp/g
1466 !-- take out cloud liquid water for detrainment
1467 detup=up_massdetro(i,k)
1468 dz=zo_cup(i,k)-zo_cup(i,k-1)
1469 if(k.lt.ktop(i)) dellaqc(i,k) = zuo(i,k)*c1d(i,k)*qrco(i,k)*dz/dp*g
1470 ! dellaqc(i,k)= detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
1471 if(k.eq.ktop(i))dellaqc(i,k)= detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
1473 G_rain= 0.5*(pwo (i,k)+pwo (i,k+1))*g/dp
1474 E_dn = -0.5*(pwdo(i,k)+pwdo(i,k+1))*g/dp*edto(i) ! pwdo < 0 and E_dn must > 0
1475 !-- condensation source term = detrained + flux divergence of
1476 !-- cloud liquid water (qrco) + converted to rain
1478 C_up = dellaqc(i,k)+(zuo(i,k+1)* qrco(i,k+1) - &
1479 zuo(i,k )* qrco(i,k ) )*g/dp + G_rain
1480 ! C_up = dellaqc(i,k)+ G_rain
1481 !-- water vapor budget
1482 !-- = flux divergence z*(Q_c - Q_env)_up_and_down &
1483 !-- - condensation term + evaporation
1484 dellaq(i,k) =-(zuo(i,k+1)*(qco (i,k+1)-qo_cup(i,k+1) ) - &
1485 zuo(i,k )*(qco (i,k )-qo_cup(i,k ) ) )*g/dp &
1486 +(zdo(i,k+1)*(qcdo(i,k+1)-qo_cup(i,k+1) ) - &
1487 zdo(i,k )*(qcdo(i,k )-qo_cup(i,k ) ) )*g/dp*edto(i) &
1489 !- check water conservation liq+condensed (including rainfall)
1490 ! trash= trash+ (dellaq(i,k)+dellaqc(i,k)+ G_rain-E_dn)*dp/g
1496 444 format(1x,i2,1x,7e12.4) !,1x,f7.2,2x,e13.5)
1498 !--- using dellas, calculate changed environmental profiles
1506 if(ierr(i).eq.0)then
1508 XHE(I,K)=DELLAH(I,K)*MBDT+HEO(I,K)
1509 ! XQ(I,K)=max(1.e-16,(dellaqc(i,k)+DELLAQ(I,K))*MBDT+QO(I,K))
1510 XQ(I,K)=max(1.e-16,DELLAQ(I,K)*MBDT+QO(I,K))
1511 DELLAT(I,K)=(1./cp)*(DELLAH(I,K)-xlv*DELLAQ(I,K))
1512 ! XT(I,K)= (DELLAT(I,K)-xlv/cp*dellaqc(i,k))*MBDT+TN(I,K)
1513 XT(I,K)= DELLAT(I,K)*MBDT+TN(I,K)
1514 xt(i,k)=max(190.,xt(i,k))
1519 if(ierr(i).eq.0)then
1520 XHE(I,ktf)=HEO(I,ktf)
1526 !--- calculate moist static energy, heights, qes
1528 call cup_env(xz,xqes,xhe,xhes,xt,xq,po,z1, &
1529 psur,ierr,tcrit,-1, &
1533 !--- environmental values on cloud levels
1535 call cup_env_clev(xt,xqes,xq,xhe,xhes,xz,po,xqes_cup,xq_cup, &
1536 xhe_cup,xhes_cup,xz_cup,po_cup,gamma_cup,xt_cup,psur, &
1542 !**************************** static control
1544 !--- moist static energy inside cloud
1553 if(ierr(i).eq.0)then
1554 x_add = xlv*zqexec(i)+cp*ztexec(i)
1555 call get_cloud_bc(kte,xhe_cup (i,1:kte),xhkb (i),k22(i),x_add)
1556 do k=1,start_level(i)-1
1557 xhc(i,k)=xhe_cup(i,k)
1566 if(ierr(i).eq.0)then
1567 do k=start_level(i) +1,ktop(i)
1568 xhc(i,k)=(xhc(i,k-1)*xzu(i,k-1)-.5*up_massdetro(i,k-1)*xhc(i,k-1) + &
1569 up_massentro(i,k-1)*xhe(i,k-1)) / &
1570 (xzu(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
1571 xdby(i,k)=xhc(i,k)-xhes_cup(i,k)
1574 xHC (i,K)=xhes_cup(i,k)
1581 !--- workfunctions for updraft
1583 call cup_up_aa0(xaa0,xz,xzu,xdby,GAMMA_CUP,xt_cup, &
1588 if(ierr(i).eq.0)then
1589 xaa0_ens(i,1)=xaa0(i)
1594 pr_ens(i,nens3)=pr_ens(i,nens3) &
1595 +pwo(i,k)+edto(i)*pwdo(i,k)
1597 else if(nens3.eq.8)then
1598 pr_ens(i,nens3)=pr_ens(i,nens3)+ &
1599 pwo(i,k)+edto(i)*pwdo(i,k)
1601 else if(nens3.eq.9)then
1602 pr_ens(i,nens3)=pr_ens(i,nens3) &
1603 + pwo(i,k)+edto(i)*pwdo(i,k)
1605 pr_ens(i,nens3)=pr_ens(i,nens3)+ &
1606 pwo(i,k) +edto(i)*pwdo(i,k)
1610 if(pr_ens(i,7).lt.1.e-6)then
1612 ierrc(i)="total normalized condensate too small"
1618 if(pr_ens(i,nens3).lt.1.e-5)then
1626 !--- LARGE SCALE FORCING
1629 !------- CHECK wether aa0 should have been zero, assuming this
1630 ! ensemble is chosen
1638 CALL cup_MAXIMI(HEO_CUP,2,KBMAX,K22x,ierr, &
1642 call cup_kbcon(ierrc,cap_max_increment,iloop,k22x,kbconx,heo_cup, &
1643 heso_cup,hkbo,ierr2,kbmax,po_cup,cap_max, &
1647 z_cup,entr_rate,heo,imid)
1649 call cup_kbcon(ierrc,cap_max_increment,iloop,k22x,kbconx,heo_cup, &
1650 heso_cup,hkbo,ierr3,kbmax,po_cup,cap_max, &
1654 z_cup,entr_rate,heo,imid)
1656 !--- calculate cloud base mass flux
1663 dq=(qo_cup(i,k+1)-qo_cup(i,k))
1664 mconv(i)=mconv(i)+omeg(i,k)*dq/g
1667 call cup_forcing_ens_3d(closure_n,xland1,aa0,aa1,xaa0_ens,mbdt,dtime, &
1668 ierr,ierr2,ierr3,xf_ens,axx,forcing, &
1669 maxens3,mconv,rand_clos, &
1670 po_cup,ktop,omeg,zdo,k22,zuo,pr_ens,edto,kbcon, &
1674 dicycle,tau_ecmwf,aa1_bl,xf_dicycle)
1678 if(ierr(i).eq.0)then
1679 dellat_ens (i,k,1)=dellat(i,k)
1680 dellaq_ens (i,k,1)=dellaq(i,k)
1681 dellaqc_ens(i,k,1)=dellaqc(i,k)
1682 pwo_ens (i,k,1)=pwo(i,k) !+edto(i)*pwdo(i,k)
1684 dellat_ens (i,k,1)=0.
1685 dellaq_ens (i,k,1)=0.
1686 dellaqc_ens(i,k,1)=0.
1695 if(imid.eq.1 .and. ichoice .le.2)then
1700 if(ierr(i).eq.0)then
1703 if(k22(i).lt.kpbl(i)+1)then
1705 blqe=blqe+100.*dhdt(i,k)*(po_cup(i,k)-po_cup(i,k+1))/g
1707 trash=max((hco(i,kbcon(i))-heo_cup(i,kbcon(i))),1.e1)
1708 xff_mid(i,1)=max(0.,blqe/trash)
1709 xff_mid(i,1)=min(0.1,xff_mid(i,1))
1711 xff_mid(i,2)=min(0.1,.03*zws(i))
1715 call cup_output_ens_3d(xff_mid,xf_ens,ierr,dellat_ens,dellaq_ens, &
1717 outq,outqc,zuo,pre,pwo_ens,xmb,ktop, &
1718 edto,pwdo,'deep',ierr2,ierr3, &
1719 po_cup,pr_ens,maxens3, &
1720 sig,closure_n,xland1,xmbm_in,xmbs_in, &
1721 ichoice,imid,ipr,itf,ktf, &
1723 dicycle,xf_dicycle )
1726 if(ierr(i).eq.0 .and.pre(i).gt.0.) then
1727 PRE(I)=MAX(PRE(I),0.)
1730 outu(i,k)=dellu(i,k)*xmb(i)
1731 outv(i,k)=dellv(i,k)*xmb(i)
1733 elseif(ierr(i).ne.0 .or. pre(i).eq.0.)then
1744 ! rain evaporation as in SAS
1746 if(irainevap.eq.1)then
1754 if(ierr(i).eq.0)then
1755 do k = ktop(i), 1, -1
1756 rain = pwo(i,k) + edto(i) * pwdo(i,k)
1757 rntot(i) = rntot(i) + rain * xmb(i)* .001 * dtime
1764 if(ierr(i).eq.0)then
1765 evef = edt(i) * evfact
1766 if(xland(i).gt.0.5 .and. xland(i).lt.1.5) evef=edt(i) * evfactl
1767 do k = ktop(i), 1, -1
1768 rain = pwo(i,k) + edto(i) * pwdo(i,k)
1769 rn(i) = rn(i) + rain * xmb(i) * .001 * dtime
1771 q1=qo(i,k)+(outq(i,k))*dtime
1772 t1=tn(i,k)+(outt(i,k))*dtime
1773 qcond(i) = evef * (q1 - qeso(i,k)) &
1774 & / (1. + el2orc * qeso(i,k) / t1**2)
1775 dp = -100.*(p_cup(i,k+1)-p_cup(i,k))
1776 if(rn(i).gt.0. .and. qcond(i).lt.0.) then
1777 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dtime*rn(i))))
1778 qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
1779 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
1781 if(rn(i).gt.0..and.qcond(i).lt.0..and. &
1782 & delq2(i).gt.rntot(i)) then
1783 qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
1786 if(rn(i).gt.0..and.qevap(i).gt.0.) then
1787 outq(i,k) = outq(i,k) + qevap(i)/dtime
1788 outt(i,k) = outt(i,k) - elocp * qevap(i)/dtime
1789 rn(i) = max(0.,rn(i) - .001 * qevap(i) * dp / g)
1790 pre(i) = pre(i) - qevap(i) * dp /g/dtime
1791 PRE(I)=MAX(PRE(I),0.)
1792 delqev(i) = delqev(i) + .001*dp*qevap(i)/g
1796 ! pre(i)=1000.*rn(i)/dtime
1801 ! since kinetic energy is being dissipated, add heating accordingly (from ECMWF)
1804 if(ierr(i).eq.0) then
1808 dp=(po_cup(i,k)-po_cup(i,k+1))*100.
1809 !total KE dissiptaion estimate
1810 dts= dts -(outu(i,k)*us(i,k)+outv(i,k)*vs(i,k))*dp/g
1811 ! fpi needed for calcualtion of conversion to pot. energyintegrated
1812 fpi = fpi +sqrt(outu(i,k)*outu(i,k) + outv(i,k)*outv(i,k))*dp
1816 fp= sqrt((outu(i,k)*outu(i,k)+outv(i,k)*outv(i,k)))/fpi
1817 outt(i,k)=outt(i,k)+fp*dts*g/cp
1823 #if ( WRF_CHEM == 1 )
1824 !--- calculate in-cloud/updraft air temperature
1826 if (ierr(i)==0) then
1828 tempco(i,k)=(1./cp)*(hco(i,k)-g*zo_cup(i,k)-xlv*qco(i,k))
1832 tempco(i,k)=tn_cup(i,k)
1836 if ((chem_conv_tr>0).and.(chemopt>0)) then
1837 call ctrans_gf(numgas,num_chem,chem2d,chemopt,0 &
1838 ,outchemt,conv_tr_wetscav,conv_tr_aqchem &
1840 ,zuo,zdo,pwo,pwdo,pwevo,pwavo &
1841 ,up_massentro,up_massdetro &
1842 ,dd_massentro,dd_massdetro &
1844 ,ktop,k22,kbcon,jmin &
1846 ,itf,ktf,its,ite,kts,kte &
1849 if ((chem_conv_tr>0).and.(traceropt>0)) then
1850 call ctrans_gf(0,num_tracer,tracer2d,0,traceropt &
1853 ,zuo,zdo,pwo,pwdo,pwevo,pwavo &
1854 ,up_massentro,up_massdetro &
1855 ,dd_massentro,dd_massdetro &
1857 ,ktop,k22,kbcon,jmin &
1859 ,itf,ktf,its,ite,kts,kte &
1865 !---------------------------done------------------------------
1868 END SUBROUTINE CUP_gf
1871 SUBROUTINE cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, &
1872 pw,ccn,pwev,edtmax,edtmin,edtc,psum2,psumh, &
1873 rho,aeroevap,itf,ktf, &
1883 ! ierr error value, maybe modified in this routine
1885 real, dimension (its:ite,kts:kte) &
1888 real, dimension (its:ite,1) &
1891 real, dimension (its:ite) &
1894 real, dimension (its:ite) &
1896 pwav,pwev,ccn,psum2,psumh,edtmax,edtmin
1897 integer, dimension (its:ite) &
1900 integer, dimension (its:ite) &
1901 ,intent (inout) :: &
1904 ! local variables in this routine
1908 real einc,pef,pefb,prezk,zkbc
1909 real, dimension (its:ite) :: &
1911 real :: prop_c,pefc,aeroadd,alpha3,beta3
1918 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
1920 ! */ calculate an average wind shear over the depth of the cloud
1933 IF(ierr(i).ne.0)GO TO 62
1934 if (kk .le. min0(ktop(i),ktf) .and. kk .ge. kbcon(i)) then
1936 (abs((us(i,kk+1)-us(i,kk))/(z(i,kk+1)-z(i,kk))) &
1937 + abs((vs(i,kk+1)-vs(i,kk))/(z(i,kk+1)-z(i,kk)))) * &
1938 (p(i,kk) - p(i,kk+1))
1939 sdp(i) = sdp(i) + p(i,kk) - p(i,kk+1)
1941 if (kk .eq. ktf-1)vshear(i) = 1.e3 * vws(i) / sdp(i)
1945 IF(ierr(i).eq.0)then
1946 pef=(1.591-.639*VSHEAR(I)+.0953*(VSHEAR(I)**2) &
1947 -.00496*(VSHEAR(I)**3))
1948 if(pef.gt.0.9)pef=0.9
1949 if(pef.lt.0.1)pef=0.1
1951 !--- cloud base precip efficiency
1953 zkbc=z(i,kbcon(i))*3.281e-3
1956 prezk=.96729352+zkbc*(-.70034167+zkbc*(.162179896+zkbc &
1957 *(- 1.2569798E-2+zkbc*(4.2772E-4-zkbc*5.44E-6))))
1963 if(pefb.gt.0.9)pefb=0.9
1964 if(pefb.lt.0.1)pefb=0.1
1965 EDT(I)=1.-.5*(pefb+pef)
1966 if(aeroevap.gt.1)then
1967 aeroadd=(ccnclean**beta3)*((psumh(i))**(alpha3-1)) !*1.e6
1969 prop_c=.5*(pefb+pef)/aeroadd
1970 aeroadd=(ccn(i)**beta3)*((psum2(i))**(alpha3-1)) !*1.e6
1971 aeroadd=prop_c*aeroadd
1973 if(pefc.gt.0.9)pefc=0.9
1974 if(pefc.lt.0.1)pefc=0.1
1976 if(aeroevap.eq.2)EDT(I)=1.-.25*(pefb+pef+2.*pefc)
1980 !--- edt here is 1-precipeff!
1982 edtc(i,1)=edt(i)-einc
1986 IF(ierr(i).eq.0)then
1987 EDTC(I,1)=-EDTC(I,1)*pwav(I)/PWEV(I)
1988 IF(EDTC(I,1).GT.edtmax(i))EDTC(I,1)=edtmax(i)
1989 IF(EDTC(I,1).LT.edtmin(i))EDTC(I,1)=edtmin(i)
1993 END SUBROUTINE cup_dd_edt
1996 SUBROUTINE cup_dd_moisture(ierrc,zd,hcd,hes_cup,qcd,qes_cup, &
1997 pwd,q_cup,z_cup,dd_massentr,dd_massdetr,jmin,ierr, &
1998 gamma_cup,pwev,bu,qrcd, &
2009 ! cdd= detrainment function
2010 ! q = environmental q on model levels
2011 ! q_cup = environmental q on model cloud levels
2012 ! qes_cup = saturation q on model cloud levels
2013 ! hes_cup = saturation h on model cloud levels
2014 ! hcd = h in model cloud
2016 ! zd = normalized downdraft mass flux
2017 ! gamma_cup = gamma on model cloud levels
2018 ! mentr_rate = entrainment rate
2019 ! qcd = cloud q (including liquid water) after entrainment
2020 ! qrch = saturation q in cloud
2021 ! pwd = evaporate at that level
2022 ! pwev = total normalized integrated evaoprate (I2)
2023 ! entr= entrainment rate
2025 real, dimension (its:ite,kts:kte) &
2027 zd,hes_cup,hcd,qes_cup,q_cup,z_cup, &
2028 dd_massentr,dd_massdetr,gamma_cup,q,he
2032 integer, dimension (its:ite) &
2035 integer, dimension (its:ite) &
2036 ,intent (inout) :: &
2038 real, dimension (its:ite,kts:kte)&
2041 real, dimension (its:ite)&
2044 character*50 :: ierrc(its:ite)
2046 ! local variables in this routine
2069 IF(ierr(I).eq.0)then
2071 DZ=Z_cup(i,K+1)-Z_cup(i,K)
2073 DH=HCD(I,k)-HES_cup(I,K)
2075 QRCD(I,K)=(qes_cup(i,k)+(1./XLV)*(GAMMA_cup(i,k) &
2076 /(1.+GAMMA_cup(i,k)))*DH)
2078 qrcd(i,k)=qes_cup(i,k)
2080 pwd(i,jmin(i))=zd(i,jmin(i))*min(0.,qcd(i,k)-qrcd(i,k))
2082 pwev(i)=pwev(i)+pwd(i,jmin(i)) ! *dz
2085 do ki=jmin(i)-1,1,-1
2086 DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
2087 ! QCD(i,Ki)=(qCD(i,Ki+1)*(1.-.5*CDD(i,Ki+1)*DZ) &
2088 ! +entr*DZ*q(i,Ki) &
2089 ! )/(1.+entr*DZ-.5*CDD(i,Ki+1)*DZ)
2091 !print*,"i=",i," k=",ki," qcd(i,ki+1)=",qcd(i,ki+1)
2092 !print*,"zd=",zd(i,ki+1)," dd_ma=",dd_massdetr(i,ki)," q=",q(i,ki)
2093 !JOE-added check for non-zero denominator:
2094 denom=zd(i,ki+1)-.5*dd_massdetr(i,ki)+dd_massentr(i,ki)
2095 if(denom.lt.1.e-8)then
2099 qcd(i,ki)=(qcd(i,ki+1)*zd(i,ki+1) &
2100 -.5*dd_massdetr(i,ki)*qcd(i,ki+1)+ &
2101 dd_massentr(i,ki)*q(i,ki)) / &
2102 (zd(i,ki+1)-.5*dd_massdetr(i,ki)+dd_massentr(i,ki))
2104 !--- to be negatively buoyant, hcd should be smaller than hes!
2105 !--- ideally, dh should be negative till dd hits ground, but that is not always
2108 DH=HCD(I,ki)-HES_cup(I,Ki)
2110 QRCD(I,Ki)=qes_cup(i,ki)+(1./XLV)*(GAMMA_cup(i,ki) &
2111 /(1.+GAMMA_cup(i,ki)))*DH
2112 dqeva=qcd(i,ki)-qrcd(i,ki)
2115 qrcd(i,ki)=qcd(i,ki)
2117 pwd(i,ki)=zd(i,ki)*dqeva
2118 qcd(i,ki)=qrcd(i,ki)
2119 pwev(i)=pwev(i)+pwd(i,ki) ! *dz
2120 ! if(iloop.eq.1.and.i.eq.102.and.j.eq.62)then
2121 ! print *,'in cup_dd_moi ', hcd(i,ki),HES_cup(I,Ki),dh,dqeva
2125 !--- end loop over i
2126 if( (pwev(i).eq.0.) .and. (iloop.eq.1))then
2127 ! print *,'problem with buoy in cup_dd_moisture',i
2129 ierrc(i)="problem with buoy in cup_dd_moisture"
2131 if(BU(I).GE.0.and.iloop.eq.1)then
2132 ! print *,'problem with buoy in cup_dd_moisture',i
2134 ierrc(i)="problem2 with buoy in cup_dd_moisture"
2139 END SUBROUTINE cup_dd_moisture
2141 SUBROUTINE cup_env(z,qes,he,hes,t,q,p,z1, &
2142 psur,ierr,tcrit,itest, &
2153 ! ierr error value, maybe modified in this routine
2154 ! q = environmental mixing ratio
2155 ! qes = environmental saturation mixing ratio
2156 ! t = environmental temp
2157 ! tv = environmental virtual temp
2158 ! p = environmental pressure
2159 ! z = environmental heights
2160 ! he = environmental moist static energy
2161 ! hes = environmental saturation moist static energy
2162 ! psur = surface pressure
2163 ! z1 = terrain elevation
2166 real, dimension (its:ite,kts:kte) &
2169 real, dimension (its:ite,kts:kte) &
2172 real, dimension (its:ite,kts:kte) &
2173 ,intent (inout) :: &
2175 real, dimension (its:ite) &
2178 integer, dimension (its:ite) &
2179 ,intent (inout) :: &
2185 ! local variables in this routine
2190 ! real, dimension (1:2) :: AE,BE,HT
2191 real, dimension (its:ite,kts:kte) :: tv
2192 real :: tcrit,e,tvbar
2193 ! real, external :: satvap
2199 ! BE(1)=.622*HT(1)/.286
2200 ! AE(1)=BE(1)/273.+ALOG(610.71)
2201 ! BE(2)=.622*HT(2)/.286
2202 ! AE(2)=BE(2)/273.+ALOG(610.71)
2205 if(ierr(i).eq.0)then
2206 !Csgb - IPH is for phase, dependent on TCRIT (water or ice)
2208 ! IF(T(I,K).LE.TCRIT)IPH=2
2209 ! print *, 'AE(IPH),BE(IPH) = ',AE(IPH),BE(IPH),AE(IPH)-BE(IPH),T(i,k),i,k
2210 ! E=EXP(AE(IPH)-BE(IPH)/T(I,K))
2211 ! print *, 'P, E = ', P(I,K), E
2212 ! QES(I,K)=.622*E/(100.*P(I,K)-E)
2214 qes(i,k)=0.622*e/max(1.e-8,(p(i,k)-e))
2215 IF(QES(I,K).LE.1.E-16)QES(I,K)=1.E-16
2216 IF(QES(I,K).LT.Q(I,K))QES(I,K)=Q(I,K)
2217 ! IF(Q(I,K).GT.QES(I,K))Q(I,K)=QES(I,K)
2218 TV(I,K)=T(I,K)+.608*Q(I,K)*T(I,K)
2223 !--- z's are calculated with changed h's and q's and t's
2226 if(itest.eq.1 .or. itest.eq.0)then
2228 if(ierr(i).eq.0)then
2229 Z(I,1)=max(0.,Z1(I))-(ALOG(P(I,1))- &
2230 ALOG(PSUR(I)))*287.*TV(I,1)/9.81
2234 ! --- calculate heights
2237 if(ierr(i).eq.0)then
2238 TVBAR=.5*TV(I,K)+.5*TV(I,K-1)
2239 Z(I,K)=Z(I,K-1)-(ALOG(P(I,K))- &
2240 ALOG(P(I,K-1)))*287.*TVBAR/9.81
2244 else if(itest.eq.2)then
2247 if(ierr(i).eq.0)then
2248 z(i,k)=(he(i,k)-1004.*t(i,k)-2.5e6*q(i,k))/9.81
2249 z(i,k)=max(1.e-3,z(i,k))
2253 else if(itest.eq.-1)then
2256 !--- calculate moist static energy - HE
2257 ! saturated moist static energy - HES
2261 if(ierr(i).eq.0)then
2262 if(itest.le.0)HE(I,K)=9.81*Z(I,K)+1004.*T(I,K)+2.5E06*Q(I,K)
2263 HES(I,K)=9.81*Z(I,K)+1004.*T(I,K)+2.5E06*QES(I,K)
2264 IF(HE(I,K).GE.HES(I,K))HE(I,K)=HES(I,K)
2269 END SUBROUTINE cup_env
2272 SUBROUTINE cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup, &
2273 he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
2285 ! ierr error value, maybe modified in this routine
2286 ! q = environmental mixing ratio
2287 ! q_cup = environmental mixing ratio on cloud levels
2288 ! qes = environmental saturation mixing ratio
2289 ! qes_cup = environmental saturation mixing ratio on cloud levels
2290 ! t = environmental temp
2291 ! t_cup = environmental temp on cloud levels
2292 ! p = environmental pressure
2293 ! p_cup = environmental pressure on cloud levels
2294 ! z = environmental heights
2295 ! z_cup = environmental heights on cloud levels
2296 ! he = environmental moist static energy
2297 ! he_cup = environmental moist static energy on cloud levels
2298 ! hes = environmental saturation moist static energy
2299 ! hes_cup = environmental saturation moist static energy on cloud levels
2300 ! gamma_cup = gamma on cloud levels
2301 ! psur = surface pressure
2302 ! z1 = terrain elevation
2305 real, dimension (its:ite,kts:kte) &
2308 real, dimension (its:ite,kts:kte) &
2310 qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup
2311 real, dimension (its:ite) &
2314 integer, dimension (its:ite) &
2315 ,intent (inout) :: &
2318 ! local variables in this routine
2339 if(ierr(i).eq.0)then
2340 qes_cup(i,k)=.5*(qes(i,k-1)+qes(i,k))
2341 q_cup(i,k)=.5*(q(i,k-1)+q(i,k))
2342 hes_cup(i,k)=.5*(hes(i,k-1)+hes(i,k))
2343 he_cup(i,k)=.5*(he(i,k-1)+he(i,k))
2344 if(he_cup(i,k).gt.hes_cup(i,k))he_cup(i,k)=hes_cup(i,k)
2345 z_cup(i,k)=.5*(z(i,k-1)+z(i,k))
2346 p_cup(i,k)=.5*(p(i,k-1)+p(i,k))
2347 t_cup(i,k)=.5*(t(i,k-1)+t(i,k))
2348 gamma_cup(i,k)=(xlv/cp)*(xlv/(r_v*t_cup(i,k) &
2349 *t_cup(i,k)))*qes_cup(i,k)
2354 if(ierr(i).eq.0)then
2355 qes_cup(i,1)=qes(i,1)
2357 ! hes_cup(i,1)=hes(i,1)
2358 ! he_cup(i,1)=he(i,1)
2359 hes_cup(i,1)=9.81*z1(i)+1004.*t(i,1)+2.5e6*qes(i,1)
2360 he_cup(i,1)=9.81*z1(i)+1004.*t(i,1)+2.5e6*q(i,1)
2361 z_cup(i,1)=.5*(z(i,1)+z1(i))
2362 p_cup(i,1)=.5*(p(i,1)+psur(i))
2366 gamma_cup(i,1)=xlv/cp*(xlv/(r_v*t_cup(i,1) &
2367 *t_cup(i,1)))*qes_cup(i,1)
2371 END SUBROUTINE cup_env_clev
2373 SUBROUTINE cup_forcing_ens_3d(closure_n,xland,aa0,aa1,xaa0,mbdt,dtime,ierr,ierr2,ierr3,&
2374 xf_ens,axx,forcing,maxens3,mconv,rand_clos, &
2375 p_cup,ktop,omeg,zd,k22,zu,pr_ens,edt,kbcon, &
2379 dicycle,tau_ecmwf,aa1_bl,xf_dicycle )
2387 integer, intent (in ) :: &
2390 ! ierr error value, maybe modified in this routine
2391 ! pr_ens = precipitation ensemble
2392 ! xf_ens = mass flux ensembles
2393 ! massfln = downdraft mass flux ensembles used in next timestep
2394 ! omeg = omega from large scale model
2395 ! mconv = moisture convergence from large scale model
2396 ! zd = downdraft normalized mass flux
2397 ! zu = updraft normalized mass flux
2398 ! aa0 = cloud work function without forcing effects
2399 ! aa1 = cloud work function with forcing effects
2400 ! xaa0 = cloud work function with cloud effects (ensemble dependent)
2402 ! dir = "storm motion"
2403 ! mbdt = arbitrary numerical parameter
2404 ! dtime = dt over which forcing is applied
2405 ! iact_gr_old = flag to tell where convection was active
2406 ! kbcon = LFC of parcel from k22
2407 ! k22 = updraft originating level
2408 ! ichoice = flag if only want one closure (usually set to zero!)
2410 real, dimension (its:ite,1:maxens3) &
2411 ,intent (inout) :: &
2413 real, dimension (its:ite,1:maxens3) &
2414 ,intent (inout ) :: &
2416 real, dimension (its:ite,kts:kte) &
2419 real, dimension (its:ite,kts:kte) &
2422 real, dimension (its:ite,1) &
2425 real, dimension (its:ite,4) &
2428 real, dimension (its:ite) &
2431 real, dimension (its:ite) &
2434 real, dimension (its:ite) &
2435 ,intent (inout) :: &
2443 integer, dimension (its:ite) &
2444 ,intent (inout ) :: &
2446 integer, dimension (its:ite) &
2449 integer, dimension (its:ite) &
2450 ,intent (inout) :: &
2455 integer, intent(IN) :: DICYCLE
2456 real, intent(IN) , dimension (its:ite) :: aa1_bl,tau_ecmwf
2457 real, intent(INOUT), dimension (its:ite) :: xf_dicycle
2458 real, intent(INOUT), dimension (its:ite,10) :: forcing
2462 ! local variables in this routine
2465 real, dimension (1:maxens3) :: &
2467 real, dimension (1) :: &
2471 ! integer, parameter :: mkxcrt=15
2472 ! real, dimension(1:mkxcrt) :: &
2473 ! pcrit,acrit,acritt
2474 integer, dimension (its:ite) :: kloc
2476 a1,a_ave,xff0,xomg!,aclim1,aclim2,aclim3,aclim4
2478 real, dimension (its:ite) :: ens_adj
2486 !--- LARGE SCALE FORCING
2490 IF(ierr(i).eq.0)then
2491 kloc(i)=maxloc(zu(i,:),1)
2493 !ss --- comment out adjustment over ocean
2494 !ss if(ierr2(i).gt.0.and.ierr3(i).eq.0)ens_adj(i)=0.666 ! 2./3.
2495 !ss if(ierr2(i).gt.0.and.ierr3(i).gt.0)ens_adj(i)=0.333
2500 a_ave=min(a_ave,aa1(i))
2503 xff0= (AA1(I)-AA0(I))/DTIME
2504 xff_ens3(1)=max(0.,(AA1(I)-AA0(I))/dtime)
2505 xff_ens3(2)=max(0.,(AA1(I)-AA0(I))/dtime)
2506 xff_ens3(3)=max(0.,(AA1(I)-AA0(I))/dtime)
2507 xff_ens3(16)=max(0.,(AA1(I)-AA0(I))/dtime)
2508 forcing(i,1)=xff_ens3(2)
2510 !--- omeg is in bar/s, mconv done with omeg in Pa/s
2511 ! more like Brown (1979), or Frank-Cohen (199?)
2513 ! average aaround kbcon
2520 do k=kbcon(i)-1,kbcon(i)+1
2521 if(zu(i,k).gt.0.)then
2522 xomg=xomg-omeg(i,k)/9.81/max(0.5,(1.-edt(i)*zd(i,k)/zu(i,k)))
2526 if(kk.gt.0)xff_ens3(4)=xomg/float(kk)
2530 ! xff_ens3(6)=-omeg(i,k22(i))/9.81
2531 ! do k=k22(i),kbcon(i)
2532 ! xomg=-omeg(i,k)/9.81
2533 ! if(xomg.gt.xff_ens3(6))xff_ens3(6)=xomg
2536 ! if(zu(i,kbcon(i)) > 0)xff_ens3(6)=betajb*xff_ens3(6)/zu(i,kbcon(i))
2537 xff_ens3(4)=betajb*xff_ens3(4)
2538 xff_ens3(5)=xff_ens3(4)
2539 xff_ens3(6)=xff_ens3(4)
2540 if(xff_ens3(4).lt.0.)xff_ens3(4)=0.
2541 if(xff_ens3(5).lt.0.)xff_ens3(5)=0.
2542 if(xff_ens3(6).lt.0.)xff_ens3(6)=0.
2543 xff_ens3(14)=betajb*xff_ens3(4)
2544 forcing(i,2)=xff_ens3(4)
2546 !--- more like Krishnamurti et al.; pick max and average values
2548 xff_ens3(7)= mconv(i)/max(0.5,(1.-edt(i)*zd(i,kbcon(i))/zu(i,kloc(i))))
2549 xff_ens3(8)= mconv(i)/max(0.5,(1.-edt(i)*zd(i,kbcon(i))/zu(i,kloc(i))))
2550 xff_ens3(9)= mconv(i)/max(0.5,(1.-edt(i)*zd(i,kbcon(i))/zu(i,kloc(i))))
2551 xff_ens3(15)=mconv(i)/max(0.5,(1.-edt(i)*zd(i,kbcon(i))/zu(i,kloc(i))))
2552 forcing(i,3)=xff_ens3(8)
2554 !--- more like Fritsch Chappel or Kain Fritsch (plus triggers)
2556 xff_ens3(10)=AA1(i)/tau_ecmwf(i)
2557 xff_ens3(11)=AA1(I)/tau_ecmwf(i)
2558 xff_ens3(12)=AA1(I)/tau_ecmwf(i)
2559 xff_ens3(13)=(AA1(i))/tau_ecmwf(i) !(60.*15.) !tau_ecmwf(i)
2560 ! forcing(i,4)=xff_ens3(10)
2561 !- more like Bechtold et al. (JAS 2014)
2562 if(dicycle == 1) xff_dicycle = max(0.,AA1_BL(i)/tau_ecmwf(i)) !(60.*30.) !tau_ecmwf(i)
2564 if(ichoice.eq.0)then
2579 XK(1)=(XAA0(I,1)-AA1(I))/MBDT
2582 forcing(i,6)=xaa0(i,1)
2584 if(xk(1).le.0.and.xk(1).gt.-.01*mbdt) &
2586 if(xk(1).gt.0.and.xk(1).lt.1.e-2) &
2590 !--- add up all ensembles
2593 ! over water, enfor!e small cap for some of the closures
2595 if(xland(i).lt.0.1)then
2596 if(ierr2(i).gt.0.or.ierr3(i).gt.0)then
2597 xff_ens3(1) =ens_adj(i)*xff_ens3(1)
2598 xff_ens3(2) =ens_adj(i)*xff_ens3(2)
2599 xff_ens3(3) =ens_adj(i)*xff_ens3(3)
2600 xff_ens3(4) =ens_adj(i)*xff_ens3(4)
2601 xff_ens3(5) =ens_adj(i)*xff_ens3(5)
2602 xff_ens3(6) =ens_adj(i)*xff_ens3(6)
2603 xff_ens3(7) =ens_adj(i)*xff_ens3(7)
2604 xff_ens3(8) =ens_adj(i)*xff_ens3(8)
2605 xff_ens3(9) =ens_adj(i)*xff_ens3(9)
2606 xff_ens3(10) =ens_adj(i)*xff_ens3(10)
2607 xff_ens3(11) =ens_adj(i)*xff_ens3(11)
2608 xff_ens3(12) =ens_adj(i)*xff_ens3(12)
2609 xff_ens3(13) =ens_adj(i)*xff_ens3(13)
2610 xff_ens3(14) =ens_adj(i)*xff_ens3(14)
2611 xff_ens3(15) =ens_adj(i)*xff_ens3(15)
2612 xff_ens3(16) =ens_adj(i)*xff_ens3(16)
2614 xff_dicycle = ens_adj(i)*xff_dicycle
2622 ! end water treatment
2627 !--- special treatment for stability closures
2630 if(xff_ens3(1).gt.0)xf_ens(i,1)=max(0.,-xff_ens3(1)/xk(1))
2631 if(xff_ens3(2).gt.0)xf_ens(i,2)=max(0.,-xff_ens3(2)/xk(1))
2632 if(xff_ens3(3).gt.0)xf_ens(i,3)=max(0.,-xff_ens3(3)/xk(1))
2633 if(xff_ens3(16).gt.0)xf_ens(i,16)=max(0.,-xff_ens3(16)/xk(1))
2634 xf_ens(i,1)= xf_ens(i,1)+xf_ens(i,1)*rand_clos(i,1)
2635 xf_ens(i,2)= xf_ens(i,2)+xf_ens(i,2)*rand_clos(i,1)
2636 xf_ens(i,3)= xf_ens(i,3)+xf_ens(i,3)*rand_clos(i,1)
2637 xf_ens(i,16)=xf_ens(i,16)+xf_ens(i,16)*rand_clos(i,1)
2645 !--- if iresult.eq.1, following independent of xff0
2647 xf_ens(i,4)=max(0.,xff_ens3(4))
2648 xf_ens(i,5)=max(0.,xff_ens3(5))
2649 xf_ens(i,6)=max(0.,xff_ens3(6))
2650 xf_ens(i,14)=max(0.,xff_ens3(14))
2651 a1=max(1.e-5,pr_ens(i,7))
2652 xf_ens(i,7)=max(0.,xff_ens3(7)/a1)
2653 a1=max(1.e-5,pr_ens(i,8))
2654 xf_ens(i,8)=max(0.,xff_ens3(8)/a1)
2655 ! forcing(i,7)=xf_ens(i,8)
2656 a1=max(1.e-5,pr_ens(i,9))
2657 xf_ens(i,9)=max(0.,xff_ens3(9)/a1)
2658 a1=max(1.e-3,pr_ens(i,15))
2659 xf_ens(i,15)=max(0.,xff_ens3(15)/a1)
2660 xf_ens(i,4)=xf_ens(i,4)+xf_ens(i,4)*rand_clos(i,2)
2661 xf_ens(i,5)=xf_ens(i,5)+xf_ens(i,5)*rand_clos(i,2)
2662 xf_ens(i,6)=xf_ens(i,6)+xf_ens(i,6)*rand_clos(i,2)
2663 xf_ens(i,14)=xf_ens(i,14)+xf_ens(i,14)*rand_clos(i,2)
2664 xf_ens(i,7)=xf_ens(i,7)+xf_ens(i,7)*rand_clos(i,3)
2665 xf_ens(i,8)=xf_ens(i,8)+xf_ens(i,8)*rand_clos(i,3)
2666 xf_ens(i,9)=xf_ens(i,9)+xf_ens(i,9)*rand_clos(i,3)
2667 xf_ens(i,15)=xf_ens(i,15)+xf_ens(i,15)*rand_clos(i,3)
2669 xf_ens(i,10)=max(0.,-xff_ens3(10)/xk(1))
2670 xf_ens(i,11)=max(0.,-xff_ens3(11)/xk(1))
2671 xf_ens(i,12)=max(0.,-xff_ens3(12)/xk(1))
2672 xf_ens(i,13)=max(0.,-xff_ens3(13)/xk(1))
2673 xf_ens(i,10)=xf_ens(i,10)+xf_ens(i,10)*rand_clos(i,4)
2674 xf_ens(i,11)=xf_ens(i,11)+xf_ens(i,11)*rand_clos(i,4)
2675 xf_ens(i,12)=xf_ens(i,12)+xf_ens(i,12)*rand_clos(i,4)
2676 xf_ens(i,13)=xf_ens(i,13)+xf_ens(i,13)*rand_clos(i,4)
2677 forcing(i,8)=xf_ens(i,11)
2687 xf_dicycle(i) = max(0.,-xff_dicycle /xk(1))
2688 ! forcing(i,9)=xf_dicycle(i)
2693 if(ichoice.ge.1)then
2695 xf_ens(i,1)=xf_ens(i,ichoice)
2696 xf_ens(i,2)=xf_ens(i,ichoice)
2697 xf_ens(i,3)=xf_ens(i,ichoice)
2698 xf_ens(i,4)=xf_ens(i,ichoice)
2699 xf_ens(i,5)=xf_ens(i,ichoice)
2700 xf_ens(i,6)=xf_ens(i,ichoice)
2701 xf_ens(i,7)=xf_ens(i,ichoice)
2702 xf_ens(i,8)=xf_ens(i,ichoice)
2703 xf_ens(i,9)=xf_ens(i,ichoice)
2704 xf_ens(i,10)=xf_ens(i,ichoice)
2705 xf_ens(i,11)=xf_ens(i,ichoice)
2706 xf_ens(i,12)=xf_ens(i,ichoice)
2707 xf_ens(i,13)=xf_ens(i,ichoice)
2708 xf_ens(i,14)=xf_ens(i,ichoice)
2709 xf_ens(i,15)=xf_ens(i,ichoice)
2710 xf_ens(i,16)=xf_ens(i,ichoice)
2712 elseif(ierr(i).ne.20.and.ierr(i).ne.0)then
2720 END SUBROUTINE cup_forcing_ens_3d
2722 SUBROUTINE cup_kbcon(ierrc,cap_inc,iloop_in,k22,kbcon,he_cup,hes_cup, &
2723 hkb,ierr,kbmax,p_cup,cap_max, &
2727 z_cup,entr_rate,heo,imid )
2732 ! only local wrf dimensions are need as of now in this routine
2736 jprnt,itf,ktf,imid, &
2741 ! ierr error value, maybe modified in this routine
2743 real, dimension (its:ite,kts:kte) &
2745 he_cup,hes_cup,p_cup
2746 real, dimension (its:ite) &
2748 entr_rate,ztexec,zqexec,cap_inc,cap_max
2749 real, dimension (its:ite) &
2750 ,intent (inout ) :: &
2752 integer, dimension (its:ite) &
2755 integer, dimension (its:ite) &
2756 ,intent (inout) :: &
2761 character*50 :: ierrc(its:ite)
2762 real, dimension (its:ite,kts:kte),intent (in) :: z_cup,heo
2763 integer, dimension (its:ite) :: iloop,start_level
2765 ! local variables in this routine
2771 x_add,pbcdif,plus,hetest,dz
2772 real, dimension (its:ite,kts:kte) ::hcot
2774 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
2780 ! reset iloop for mid level convection
2781 if(cap_max(i).gt.200 .and. imid.eq.1)iloop(i)=5
2783 IF(ierr(I).ne.0)GO TO 27
2784 start_level(i)=k22(i)
2786 if(iloop(i).eq.5)KBCON(I)=K22(I)
2787 ! if(iloop_in.eq.5)start_level(i)=kbcon(i)
2788 !== including entrainment for hetest
2789 hcot(i,1:start_level(i)) = HKB(I)
2790 do k=start_level(i)+1,KBMAX(i)+3
2791 dz=z_cup(i,k)-z_cup(i,k-1)
2792 hcot(i,k)= ( (1.-0.5*entr_rate(i)*dz)*hcot(i,k-1) &
2793 + entr_rate(i)*dz*heo(i,k-1) )/ &
2794 (1.+0.5*entr_rate(i)*dz)
2801 IF(KBCON(I).GT.KBMAX(i)+2)THEN
2802 if(iloop(i).ne.4)then
2804 ierrc(i)="could not find reasonable kbcon in cup_kbcon"
2809 hetest=hcot(i,kbcon(i)) !hkb(i) ! HE_cup(I,K22(I))
2810 IF(HETEST.LT.HES_cup(I,KBCON(I)))then
2814 ! cloud base pressure and max moist static energy pressure
2815 ! i.e., the depth (in mb) of the layer of negative buoyancy
2816 if(KBCON(I)-K22(I).eq.1)go to 27
2817 if(iloop(i).eq.5 .and. (KBCON(I)-K22(I)).le.2)go to 27
2818 PBCDIF=-P_cup(I,KBCON(I))+P_cup(I,K22(I))
2819 plus=max(25.,cap_max(i)-float(iloop(i)-1)*cap_inc(i))
2820 if(iloop(i).eq.4)plus=cap_max(i)
2822 ! for shallow convection, if cap_max is greater than 25, it is the pressure at pbltop
2823 if(iloop(i).eq.5)plus=150.
2824 if(iloop(i).eq.5.and.cap_max(i).gt.200)pbcdif=-P_cup(I,KBCON(I))+cap_max(i)
2825 IF(PBCDIF.le.plus)THEN
2827 ElseIF(PBCDIF.GT.plus)THEN
2830 !== since k22 has be changed, HKB has to be re-calculated
2831 x_add = xlv*zqexec(i)+cp*ztexec(i)
2832 call get_cloud_bc(kte,he_cup (i,1:kte),hkb (i),k22(i),x_add)
2834 start_level(i)=k22(i)
2835 ! if(iloop_in.eq.5)start_level(i)=kbcon(i)
2836 hcot(i,1:start_level(i)) = hkb(I)
2837 do k=start_level(i)+1,KBMAX(i)+3
2838 dz=z_cup(i,k)-z_cup(i,k-1)
2840 hcot(i,k)= ( (1.-0.5*entr_rate(i)*dz)*hcot(i,k-1) &
2841 + entr_rate(i)*dz*heo(i,k-1) )/ &
2842 (1.+0.5*entr_rate(i)*dz)
2846 if(iloop(i).eq.5)KBCON(I)=K22(I)
2847 IF(KBCON(I).GT.KBMAX(i)+2)THEN
2848 if(iloop(i).ne.4)then
2850 ierrc(i)="could not find reasonable kbcon in cup_kbcon"
2858 END SUBROUTINE cup_kbcon
2861 SUBROUTINE cup_MAXIMI(ARRAY,KS,KE,MAXX,ierr, &
2870 ! only local wrf dimensions are need as of now in this routine
2877 ! x output array with return values
2878 ! kt output array of levels
2879 ! ks,kend check-range
2880 real, dimension (its:ite,kts:kte) &
2883 integer, dimension (its:ite) &
2889 integer, dimension (its:ite) &
2892 real, dimension (its:ite) :: &
2901 if(ierr(i).eq.0)then
2906 IF(XAR.GE.X(I)) THEN
2914 END SUBROUTINE cup_MAXIMI
2917 SUBROUTINE cup_minimi(ARRAY,KS,KEND,KT,ierr, &
2926 ! only local wrf dimensions are need as of now in this routine
2933 ! x output array with return values
2934 ! kt output array of levels
2935 ! ks,kend check-range
2936 real, dimension (its:ite,kts:kte) &
2939 integer, dimension (its:ite) &
2942 integer, dimension (its:ite) &
2945 real, dimension (its:ite) :: &
2952 if(ierr(i).eq.0)then
2954 KSTOP=MAX(KS(I)+1,KEND(I))
2956 DO 100 K=KS(I)+1,KSTOP
2957 IF(ARRAY(I,K).LT.X(I)) THEN
2965 END SUBROUTINE cup_MINIMI
2968 SUBROUTINE cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
2978 ! only local wrf dimensions are need as of now in this routine
2984 ! aa0 cloud work function
2985 ! gamma_cup = gamma on model cloud levels
2986 ! t_cup = temperature (Kelvin) on model cloud levels
2987 ! dby = buoancy term
2988 ! zu= normalized updraft mass flux
2989 ! z = heights of model levels
2990 ! ierr error value, maybe modified in this routine
2992 real, dimension (its:ite,kts:kte) &
2994 z,zu,gamma_cup,t_cup,dby
2995 integer, dimension (its:ite) &
3003 integer, dimension (its:ite) &
3004 ,intent (inout) :: &
3006 real, dimension (its:ite) &
3010 ! local variables in this routine
3023 IF(ierr(i).ne.0)GO TO 100
3024 IF(K.LT.KBCON(I))GO TO 100
3025 IF(K.Gt.KTOP(I))GO TO 100
3027 da=zu(i,k)*DZ*(9.81/(1004.*( &
3028 (T_cup(I,K)))))*DBY(I,K-1)/ &
3030 ! IF(K.eq.KTOP(I).and.da.le.0.)go to 100
3031 AA0(I)=AA0(I)+max(0.,da)
3032 if(aa0(i).lt.0.)aa0(i)=0.
3035 END SUBROUTINE cup_up_aa0
3037 !====================================================================
3038 SUBROUTINE neg_check(name,j,dt,q,outq,outt,outu,outv, &
3039 outqc,pret,its,ite,kts,kte,itf,ktf)
3041 INTEGER, INTENT(IN ) :: j,its,ite,kts,kte,itf,ktf
3043 real, dimension (its:ite,kts:kte ) , &
3045 outq,outt,outqc,outu,outv
3046 real, dimension (its:ite,kts:kte ) , &
3049 real, dimension (its:ite ) , &
3052 character *(*), intent (in) :: &
3057 real :: names,scalef,thresh,qmem,qmemf,qmem2,qtest,qmem1
3060 ! first do check on vertical heating rate
3066 if(name == 'shallow')then
3076 qmem=(outt(i,k))*86400.
3077 if(qmem.gt.thresh)then
3079 qmemf=min(qmemf,qmem2)
3083 ! print *,'1',' adjusted massflux by factor ',i,j,k,qmem,qmem2,qmemf,dt
3085 if(qmem.lt.-.5*thresh*names)then
3086 qmem2=-.5*names*thresh/qmem
3087 qmemf=min(qmemf,qmem2)
3094 outq(i,k)=outq(i,k)*qmemf
3095 outt(i,k)=outt(i,k)*qmemf
3096 outu(i,k)=outu(i,k)*qmemf
3097 outv(i,k)=outv(i,k)*qmemf
3098 outqc(i,k)=outqc(i,k)*qmemf
3100 pret(i)=pret(i)*qmemf
3104 ! check whether routine produces negative q's. This can happen, since
3105 ! tendencies are calculated based on forced q's. This should have no
3106 ! influence on conservation properties, it scales linear through all
3110 ! write(14,*)'return'
3116 if(abs(qmem).gt.0. .and. q(i,k).gt.1.e-6)then
3117 qtest=q(i,k)+(outq(i,k))*dt
3118 if(qtest.lt.thresh)then
3120 ! qmem2 would be the maximum allowable tendency
3122 qmem1=abs(outq(i,k))
3123 qmem2=abs((thresh-q(i,k))/dt)
3124 qmemf=min(qmemf,qmem2/qmem1)
3130 outq(i,k)=outq(i,k)*qmemf
3131 outt(i,k)=outt(i,k)*qmemf
3132 outu(i,k)=outu(i,k)*qmemf
3133 outv(i,k)=outv(i,k)*qmemf
3134 outqc(i,k)=outqc(i,k)*qmemf
3136 pret(i)=pret(i)*qmemf
3139 END SUBROUTINE neg_check
3142 SUBROUTINE cup_output_ens_3d(xff_mid,xf_ens,ierr,dellat,dellaq,dellaqc, &
3143 outtem,outq,outqc, &
3144 zu,pre,pw,xmb,ktop, &
3145 edt,pwd,name,ierr2,ierr3,p_cup,pr_ens, &
3147 sig,closure_n,xland1,xmbm_in,xmbs_in, &
3148 ichoice,imid,ipr,itf,ktf, &
3150 dicycle,xf_dicycle )
3156 ! only local wrf dimensions are need as of now in this routine
3160 ichoice,imid,ipr,itf,ktf, &
3162 integer, intent (in ) :: &
3164 ! xf_ens = ensemble mass fluxes
3165 ! pr_ens = precipitation ensembles
3166 ! dellat = change of temperature per unit mass flux of cloud ensemble
3167 ! dellaq = change of q per unit mass flux of cloud ensemble
3168 ! dellaqc = change of qc per unit mass flux of cloud ensemble
3169 ! outtem = output temp tendency (per s)
3170 ! outq = output q tendency (per s)
3171 ! outqc = output qc tendency (per s)
3172 ! pre = output precip
3173 ! xmb = total base mass flux
3174 ! xfac1 = correction factor
3175 ! pw = pw -epsilon*pd (ensemble dependent)
3176 ! ierr error value, maybe modified in this routine
3178 real, dimension (its:ite,1:maxens3) &
3179 ,intent (inout) :: &
3181 real, dimension (its:ite,kts:kte) &
3182 ,intent (inout ) :: &
3184 real, dimension (its:ite,kts:kte) &
3187 real, dimension (its:ite) &
3189 sig,xmbm_in,xmbs_in,edt
3190 real, dimension (its:ite,2) &
3193 real, dimension (its:ite) &
3194 ,intent (inout ) :: &
3196 real, dimension (its:ite) &
3197 ,intent (inout ) :: &
3199 real, dimension (its:ite,kts:kte,1) &
3201 dellat,dellaqc,dellaq,pw
3202 integer, dimension (its:ite) &
3205 integer, dimension (its:ite) &
3206 ,intent (inout) :: &
3208 integer, intent(IN) :: DICYCLE
3209 real, intent(IN), dimension (its:ite) :: xf_dicycle
3211 ! local variables in this routine
3217 clos_wei,dtt,dp,dtq,dtqc,dtpw,dtpwd
3218 real, dimension (its:ite) :: &
3221 character *(*), intent (in) :: &
3237 IF(ierr(i).eq.0)then
3239 if(pr_ens(i,n).le.0.)then
3246 !--- calculate ensemble average mass fluxes
3252 !!!!! DEEP Convection !!!!!!!!!!
3255 if(ierr(i).eq.0)then
3260 xmb_ave(i)=xmb_ave(i)+xf_ens(i,n)
3262 xmb_ave(i)=xmb_ave(i)/float(k)
3264 if(dicycle == 2 )then
3265 xmb_ave(i)=xmb_ave(i)-max(0.,xmbm_in(i),xmbs_in(i))
3266 xmb_ave(i)=max(0.,xmb_ave(i))
3267 else if (dicycle == 1) then
3268 xmb_ave(i)=min(xmb_ave(i),xmb_ave(i) - xf_dicycle(i))
3269 xmb_ave(i)=max(0.,xmb_ave(i))
3271 ! --- Now use proper count of how many closures were actually
3272 ! used in cup_forcing_ens (including screening of some
3273 ! closures over water) to properly normalize xmb
3274 clos_wei=16./max(1.,closure_n(i))
3275 xmb_ave(i)=min(xmb_ave(i),100.)
3276 xmb(i)=clos_wei*sig(i)*xmb_ave(i)
3278 if(xmb(i) < 1.e-16)then
3286 !!!!! NOT SO DEEP Convection !!!!!!!!!!
3290 IF(ierr(i).eq.0)then
3291 ! ! first get xmb_ves, depend on ichoicee
3293 if(ichoice.eq.1 .or. ichoice.eq.2)then
3294 xmb_ave(i)=sig(i)*xff_mid(i,ichoice)
3295 else if(ichoice.gt.2)then
3299 xmb_ave(i)=xmb_ave(i)+xf_ens(i,n)
3301 xmb_ave(i)=xmb_ave(i)/float(k)
3302 else if(ichoice == 0)then
3303 xmb_ave(i)=.5*sig(i)*(xff_mid(i,1)+xff_mid(i,2))
3304 endif ! ichoice gt.2
3305 ! which dicycle method
3306 if(dicycle == 2 )then
3307 xmb(i)=max(0.,xmb_ave(i)-xmbs_in(i))
3308 else if (dicycle == 1) then
3309 xmb(i)=min(xmb_ave(i),xmb_ave(i) - xf_dicycle(i))
3310 xmb(i)=max(0.,xmb_ave(i))
3311 else if (dicycle == 0) then
3312 xmb(i)=max(0.,xmb_ave(i))
3320 IF(ierr(i).eq.0)then
3322 pwtot(i)=pwtot(i)+pw(i,k,1)
3325 dp=100.*(p_cup(i,k)-p_cup(i,k+1))/g
3328 ! necessary to drive downdraft
3329 dtpwd=-pwd(i,k)*edt(i)
3330 ! take from dellaqc first
3331 dtqc=dellaqc (i,k,1)*dp - dtpwd
3332 ! if this is negative, it needs to come from rain
3334 dtpwd=dtpwd-dellaqc(i,k,1)*dp
3336 ! if this is positive, can come from clw detrainment
3341 OUTTEM(I,K)= XMB(I)* dtt
3342 OUTQ (I,K)= XMB(I)* dtq
3343 OUTQC (I,K)= XMB(I)* dtqc
3344 xf_ens(i,:)=sig(i)*xf_ens(i,:)
3345 ! what is evaporated
3346 PRE(I)=PRE(I)-XMB(I)*dtpwd
3348 PRE(I)=-PRE(I)+XMB(I)*pwtot(i)
3353 END SUBROUTINE cup_output_ens_3d
3354 !-------------------------------------------------------
3355 SUBROUTINE cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, &
3356 p_cup,kbcon,ktop,dby,clw_all,xland1, &
3357 q,GAMMA_cup,zu,qes_cup,k22,qe_cup, &
3358 ZQEXEC,ccn,rho,c1d,t, &
3359 up_massentr,up_massdetr,psum,psumh, &
3364 real, parameter :: BDISPM = 0.366 !Berry--size dispersion (martime)
3365 REAL, PARAMETER :: BDISPC = 0.146 !Berry--size dispersion (continental)
3370 ! only local wrf dimensions are need as of now in this routine
3376 ! cd= detrainment function
3377 ! q = environmental q on model levels
3378 ! qe_cup = environmental q on model cloud levels
3379 ! qes_cup = saturation q on model cloud levels
3380 ! dby = buoancy term
3381 ! cd= detrainment function
3382 ! zu = normalized updraft mass flux
3383 ! gamma_cup = gamma on model cloud levels
3385 real, dimension (its:ite,kts:kte) &
3387 p_cup,rho,q,zu,gamma_cup,qe_cup, &
3388 up_massentr,up_massdetr,dby,qes_cup,z_cup
3389 real, dimension (its:ite) &
3392 ! entr= entrainment rate
3393 integer, dimension (its:ite) &
3395 kbcon,ktop,k22,xland1
3400 ! ierr error value, maybe modified in this routine
3402 integer, dimension (its:ite) &
3403 ,intent (inout) :: &
3405 character *(*), intent (in) :: &
3407 ! qc = cloud q (including liquid water) after entrainment
3408 ! qrch = saturation q in cloud
3409 ! qrc = liquid water content in cloud after rainout
3410 ! pw = condensate that will fall out at that level
3411 ! pwav = totan normalized integrated condensate (I1)
3412 ! c0 = conversion rate (cloud to rain)
3414 real, dimension (its:ite,kts:kte) &
3417 real, dimension (its:ite,kts:kte) :: &
3418 qch,qrcb,pwh,clw_allh,c1d,t
3419 real, dimension (its:ite) :: &
3421 real, dimension (its:ite) &
3424 real, dimension (its:ite) &
3428 ! local variables in this routine
3433 integer :: start_level(its:ite)
3435 prop_ave,qrcb_h,bdsp,dp,rhoc,qrch,qaver, &
3436 c0,dz,berryc0,q1,berryc
3439 real, dimension (kts:kte) :: &
3444 c0=.004 ! Han et al. (2016); Li et al., (2019), updated from 0.002
3447 !--- no precip for small clouds
3449 ! if(name.eq.'shallow')then
3463 if(ierr(i).eq.0)qc(i,k)=qe_cup(i,k)
3464 if(ierr(i).eq.0)qch(i,k)=qe_cup(i,k)
3472 if(ierr(i).eq.0)then
3474 call get_cloud_bc(kte,qe_cup (i,1:kte),qaver,k22(i))
3479 do k=1,start_level(i)-1
3480 qc (i,k)= qe_cup(i,k)
3481 qch (i,k)= qe_cup(i,k)
3484 ! initialize below originating air
3491 IF(ierr(i).eq.0)then
3493 ! below LFC, but maybe above LCL
3495 ! if(name == "deep" )then
3496 DO k=k22(i)+1,kbcon(i)
3497 if (t(i,k).lt.273.15) c0=c0*exp(0.07*(t(i,k)-273.15)) ! Han et al. (2016); Li et al. (2019)
3498 qc(i,k)= (qc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)* qc(i,k-1)+ &
3499 up_massentr(i,k-1)*q(i,k-1)) / &
3500 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
3502 QRCH=QES_cup(I,K)+(1./XLV)*(GAMMA_cup(i,k) &
3503 /(1.+GAMMA_cup(i,k)))*DBY(I,K)
3504 if(k.lt.kbcon(i))qrch=qc(i,k)
3505 if(qc(i,k).gt.qrch)then
3506 DZ=Z_cup(i,K)-Z_cup(i,K-1)
3507 QRC(I,K)=(QC(I,K)-QRCH)/(1.+c0*DZ)
3508 PW(i,k)=c0*dz*QRC(I,K)*zu(i,k)
3509 qc(i,k)=qrch+qrc(i,k)
3510 clw_all(i,k)=qrc(i,k)
3517 DO k=kbcon(i)+1,ktop(i)
3519 if (t(i,k).lt.273.15) c0=c0*exp(0.07*(t(i,k)-273.15)) ! Han et al. (2016); Li et al. (2019)
3520 denom=zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1)
3521 if(denom.lt.1.e-8)then
3527 rhoc=.5*(rho(i,k)+rho(i,k-1))
3528 DZ=Z_cup(i,K)-Z_cup(i,K-1)
3529 DP=p_cup(i,K)-p_cup(i,K-1)
3531 !--- saturation in cloud, this is what is allowed to be in it
3533 QRCH=QES_cup(I,K)+(1./XLV)*(GAMMA_cup(i,k) &
3534 /(1.+GAMMA_cup(i,k)))*DBY(I,K)
3536 !------ 1. steady state plume equation, for what could
3537 !------ be in cloud without condensation
3540 qc(i,k)= (qc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)* qc(i,k-1)+ &
3541 up_massentr(i,k-1)*q(i,k-1)) / &
3542 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
3543 qch(i,k)= (qch(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*qch(i,k-1)+ &
3544 up_massentr(i,k-1)*q(i,k-1)) / &
3545 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
3547 if(qc(i,k).le.qrch)then
3550 if(qch(i,k).le.qrch)then
3554 !------- Total condensed water before rainout
3556 clw_all(i,k)=max(0.,QC(I,K)-QRCH)
3557 QRC(I,K)=max(0.,(QC(I,K)-QRCH)) ! /(1.+C0*DZ*zu(i,k))
3558 clw_allh(i,k)=max(0.,QCH(I,K)-QRCH)
3559 QRCB(I,K)=max(0.,(QCH(I,K)-QRCH)) ! /(1.+C0*DZ*zu(i,k))
3560 IF(autoconv.eq.2) then
3566 ! first calculate for average conditions, used in cup_dd_edt!
3567 ! this will also determine proportionality constant prop_b, which, if applied,
3568 ! would give the same results as c0 under these conditions
3570 q1=1.e3*rhoc*qrcb(i,k) ! g/m^3 ! g[h2o]/cm^3
3571 berryc0=q1*q1/(60.0*(5.0 + 0.0366*CCNclean/ &
3572 ( q1 * BDSP) ) ) !/(
3573 qrcb_h=((QCH(I,K)-QRCH)*zu(i,k)-qrcb(i,k-1)*(.5*up_massdetr(i,k-1)))/ &
3574 (zu(i,k)+.5*up_massdetr(i,k-1)+c0*dz*zu(i,k))
3575 prop_b(k)=c0*qrcb_h*zu(i,k)/(1.e-3*berryc0)
3576 pwh(i,k)=zu(i,k)*1.e-3*berryc0*dz*prop_b(k) ! 2.
3578 qrcb(i,k)=((QCh(I,K)-QRCH)*zu(i,k)-pwh(i,k)-qrcb(i,k-1)*(.5*up_massdetr(i,k-1)))/ &
3579 (zu(i,k)+.5*up_massdetr(i,k-1))
3580 if(qrcb(i,k).lt.0.)then
3581 berryc0=(qrcb(i,k-1)*(.5*up_massdetr(i,k-1))-(QCh(I,K)-QRCH)*zu(i,k))/zu(i,k)*1.e-3*dz*prop_b(k)
3582 pwh(i,k)=zu(i,k)*1.e-3*berryc0*dz*prop_b(k)
3585 QCh(I,K)=QRCb(I,K)+qrch
3586 PWAVH(I)=PWAVH(I)+pwh(I,K)
3587 Psumh(I)=Psumh(I)+clw_allh(I,K)*zu(i,k) *dz
3589 ! then the real berry
3591 q1=1.e3*rhoc*qrc(i,k) ! g/m^3 ! g[h2o]/cm^3
3592 berryc0=q1*q1/(60.0*(5.0 + 0.0366*CCN(i)/ &
3593 ( q1 * BDSP) ) ) !/(
3594 berryc0=1.e-3*berryc0*dz*prop_b(k) ! 2.
3596 qrc(i,k)=((QC(I,K)-QRCH)*zu(i,k)-zu(i,k)*berryc0-qrc(i,k-1)*(.5*up_massdetr(i,k-1)))/ &
3597 (zu(i,k)+.5*up_massdetr(i,k-1))
3598 if(qrc(i,k).lt.0.)then
3599 berryc0=((QC(I,K)-QRCH)*zu(i,k)-qrc(i,k-1)*(.5*up_massdetr(i,k-1)))/zu(i,k)
3602 pw(i,k)=berryc0*zu(i,k)
3603 QC(I,K)=QRC(I,K)+qrch
3605 ! if not running with berry at all, do the following
3610 pw(i,k)=(QC(I,K)-QRCH)*zu(i,k)
3611 if(pw(i,k).lt.0.)pw(i,k)=0.
3613 QRC(I,K)=(QC(I,K)-QRCH)/(1.+(c1d(i,k)+C0)*DZ)
3614 PW(i,k)=c0*dz*QRC(I,K)*zu(i,k)
3615 if(qrc(i,k).lt.0)then
3620 QC(I,K)=QRC(I,K)+qrch
3622 PWAV(I)=PWAV(I)+PW(I,K)
3623 Psum(I)=Psum(I)+clw_all(I,K)*zu(i,k) *dz
3624 enddo ! k=kbcon,ktop
3625 ! do not include liquid/ice in qc
3626 do k=k22(i)+1,ktop(i)
3627 qc(i,k)=qc(i,k)-qrc(i,k)
3631 !--- integrated normalized ondensate
3637 prop_ave=prop_ave+prop_b(k)
3638 if(prop_b(k).gt.0)iprop=iprop+1
3642 END SUBROUTINE cup_up_moisture
3644 !--------------------------------------------------------------------
3646 REAL FUNCTION satvap(temp2)
3648 real :: temp2, temp, toot, toto, eilog, tsot, &
3649 & ewlog, ewlog2, ewlog3, ewlog4
3650 temp = temp2-273.155
3651 if (temp.lt.-20.) then !!!! ice saturation
3652 toot = 273.16 / temp2
3654 eilog = -9.09718 * (toot - 1) - 3.56654 * (log(toot) / &
3655 & log(10.)) + .876793 * (1 - toto) + (log(6.1071) / log(10.))
3656 satvap = 10 ** eilog
3658 tsot = 373.16 / temp2
3659 ewlog = -7.90298 * (tsot - 1) + 5.02808 * &
3660 & (log(tsot) / log(10.))
3661 ewlog2 = ewlog - 1.3816e-07 * &
3662 & (10 ** (11.344 * (1 - (1 / tsot))) - 1)
3663 ewlog3 = ewlog2 + .0081328 * &
3664 & (10 ** (-3.49149 * (tsot - 1)) - 1)
3665 ewlog4 = ewlog3 + (log(1013.246) / log(10.))
3666 satvap = 10 ** ewlog4
3669 !--------------------------------------------------------------------
3670 SUBROUTINE get_cloud_bc(mzp,array,x_aver,k22,add)
3672 integer, intent(in) :: mzp,k22
3673 real , intent(in) :: array(mzp)
3674 real , optional , intent(in) :: add
3675 real , intent(out) :: x_aver
3676 integer :: i,local_order_aver,order_aver
3678 !-- dimension of the average
3679 !-- a) to pick the value at k22 level, instead of a average between
3680 !-- k22-order_aver, ..., k22-1, k22 set order_aver=1)
3681 !-- b) to average between 1 and k22 => set order_aver = k22
3682 order_aver = 3 !=> average between k22, k22-1 and k22-2
3684 local_order_aver=min(k22,order_aver)
3687 do i = 1,local_order_aver
3688 x_aver = x_aver + array(k22-i+1)
3690 x_aver = x_aver/float(local_order_aver)
3691 if(present(add)) x_aver = x_aver + add
3693 end SUBROUTINE get_cloud_bc
3694 !========================================================================================
3697 SUBROUTINE rates_up_pdf(rand_vmas,ipr,name,ktop,ierr,p_cup,entr_rate_2d,hkbo,heo,heso_cup,z_cup, &
3698 xland,kstabi,k22,kbcon,its,ite,itf,kts,kte,ktf,zuo,kpbl,ktopdby,csum,pmin_lev)
3700 character *(*), intent (in) :: name
3701 integer, intent(in) :: ipr,its,ite,itf,kts,kte,ktf
3702 real, dimension (its:ite,kts:kte),intent (inout) :: entr_rate_2d,zuo
3703 real, dimension (its:ite,kts:kte),intent (in) ::p_cup, heo,heso_cup,z_cup
3704 real, dimension (its:ite),intent (in) :: hkbo,rand_vmas
3705 integer, dimension (its:ite),intent (in) :: kstabi,k22,kpbl,csum,xland,pmin_lev
3706 integer, dimension (its:ite),intent (inout) :: kbcon,ierr,ktop,ktopdby
3708 real, dimension (its:ite,kts:kte) :: hcot
3709 real :: beta_u,dz,dbythresh,dzh2,zustart,zubeg,massent,massdetr
3710 real :: dby(kts:kte),dbm(kts:kte),zux(kts:kte)
3711 real zuh2(40),zh2(40)
3712 integer :: kklev,i,kk,kbegin,k,kfinalzu
3713 integer, dimension (its:ite) :: start_level
3716 dbythresh= 1. !.0.95 ! 0.85, 0.6
3717 if(name == 'shallow' .or. name == 'mid') dbythresh=1.
3722 beta_u=max(.1,.2-float(csum(i))*.01)
3726 kbcon(i)=max(kbcon(i),2)
3727 if(ierr(i).eq.0)then
3728 start_level(i)=k22(i)
3729 zuo(i,start_level(i))=zustart
3730 zux(start_level(i))=zustart
3731 do k=start_level(i)+1,kbcon(i)
3732 dz=z_cup(i,k)-z_cup(i,k-1)
3733 massent=dz*entr_rate_2d(i,k-1)*zuo(i,k-1)
3734 massdetr=dz*1.e-9*zuo(i,k-1)
3735 zuo(i,k)=zuo(i,k-1)+massent-massdetr
3738 zubeg=zustart !zuo(i,kbcon(i))
3739 if(name .eq. 'deep')then
3741 hcot(i,start_level(i))=hkbo(i)
3742 dz=z_cup(i,start_level(i))-z_cup(i,start_level(i)-1)
3743 do k=start_level(i)+1,ktf-2
3744 dz=z_cup(i,k)-z_cup(i,k-1)
3746 hcot(i,k)=( (1.-0.5*entr_rate_2d(i,k-1)*dz)*hcot(i,k-1) &
3747 + entr_rate_2d(i,k-1)*dz*heo(i,k-1))/ &
3748 (1.+0.5*entr_rate_2d(i,k-1)*dz)
3749 if(k >= kbcon(i)) dby(k)=dby(k-1)+(hcot(i,k)-heso_cup(i,k))*dz
3750 if(k >= kbcon(i)) dbm(k)=hcot(i,k)-heso_cup(i,k)
3752 ktopdby(i)=maxloc(dby(:),1)
3753 kklev=maxloc(dbm(:),1)
3754 do k=maxloc(dby(:),1)+1,ktf-2
3755 if(dby(k).lt.dbythresh*maxval(dby))then
3765 ! at least overshoot by one level
3767 ! kfinalzu=min(max(kfinalzu,ktopdby(i)+1),ktopdby(i)+2)
3769 if(kfinalzu.le.kbcon(i)+2)then
3773 ! call get_zu_zd_pdf_fim(ipr,xland(i),zuh2,"UP",ierr(i),start_level(i), &
3774 ! call get_zu_zd_pdf_fim(rand_vmas(i),zubeg,ipr,xland(i),zuh2,"UP",ierr(i),kbcon(i), &
3775 ! kfinalzu,zuo(i,kts:kte),kts,kte,ktf,beta_u,kpbl(i),csum(i),pmin_lev(i))
3776 call get_zu_zd_pdf_fim(kklev,p_cup(i,:),rand_vmas(i),zubeg,ipr,xland(i),zuh2,"UP",ierr(i),k22(i), &
3777 kfinalzu,zuo(i,kts:kte),kts,kte,ktf,beta_u,kstabi(i),csum(i),pmin_lev(i))
3780 if ( name == 'mid' ) then
3781 if(ktop(i) <= kbcon(i)+2)then
3786 ktopdby(i)=ktop(i)+1
3787 call get_zu_zd_pdf_fim(kklev,p_cup(i,:),rand_vmas(i),zubeg,ipr,xland(i),zuh2,"MID",ierr(i),k22(i),kfinalzu,zuo(i,kts:kte),kts,kte,ktf,beta_u,kbcon(i),csum(i),pmin_lev(i))
3789 ! dzh2=(z_cup(i,ktop(i))-z_cup(i,k22(i)))/40.
3790 ! zh2(1)=z_cup(i,k22(i))
3791 ! if(zuh2(1).gt.0.1 .and. kbegin.eq.0)kbegin=1
3793 ! zh2(k)=zh2(k-1)+dzh2
3794 ! if(zuh2(k).gt.0.1 .and. kbegin.eq.0)kbegin=k
3796 ! zuo(i,k22(i))=zuh2(kbegin)
3797 ! do k=k22(i)+1,kfinalzu+1
3799 ! if(z_cup(i,k).gt.zh2(kk) .and. z_cup(i,k).le.zh2(kk+1)) then
3805 ! if(zuo(i,ktop(i)).lt.1.e-4)ktop(i)=ktop(i)-1
3808 if ( name == 'shallow' ) then
3809 if(ktop(i) <= kbcon(i)+2)then
3815 call get_zu_zd_pdf_fim(kklev,p_cup(i,:),rand_vmas(i),zubeg,ipr,xland(i),zuh2,"SH2",ierr(i),k22(i), &
3816 kfinalzu,zuo(i,kts:kte),kts,kte,ktf,beta_u,kpbl(i),csum(i),pmin_lev(i))
3823 END SUBROUTINE rates_up_pdf
3824 !-------------------------------------------------------------------------
3825 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)
3828 integer, intent(in) ::ipr,xland,kb,kklev,kt,kts,kte,ktf,kpbli,csum,pmin_lev
3829 real, intent(in) ::max_mass,zubeg
3830 real, intent(inout) :: zu(kts:kte)
3831 real, intent(in) :: p(kts:kte)
3832 real :: zuh(kts:kte),zuh2(1:40)
3833 integer, intent(inout) :: ierr
3834 character*(*), intent(in) ::draft
3837 integer :: kk,k,kb_adj,kpbli_adj
3838 real :: krmax,beta, alpha,kratio,tunning,FZU,rand_vmas,lev_start
3839 !- kb cannot be at 1st level
3841 !-- fill zu with zeros
3845 IF(draft == "UP") then
3846 lev_start=min(.9,.4+csum*.013)
3848 tunning=p(kt)+(p(kpbli)-p(kt))*lev_start
3849 tunning =min(0.9, (tunning-p(kb_adj))/(p(kt)-p(kb_adj))) !=.6
3850 tunning =max(0.2, tunning)
3851 beta = 1.3 !2.5 ! max(2.5,2./tunning)
3852 alpha= (tunning*(beta -2.)+1.)/(1.-tunning)
3853 #if ( ! defined NO_GAMMA_SUPPORT )
3854 fzu = gamma(alpha + beta)/(gamma(alpha)*gamma(beta))
3856 call wrf_error_fatal ('compiler does not support 2008 gamma intrinsic')
3858 do k=kb_adj,min(kte,kt)
3859 kratio= (p(k)-p(kb_adj))/(p(kt)-p(kb_adj)) !float(k)/float(kt+1)
3860 zu(k) = zubeg+FZU*kratio**(alpha-1.0) * (1.0-kratio)**(beta-1.0)
3863 if(maxval(zu(kts:min(ktf,kt+1))).gt.0.) &
3864 zu(kts:min(ktf,kt+1))= zu(kts:min(ktf,kt+1))/maxval(zu(kts:min(ktf,kt+1)))
3865 do k=maxloc(zu(:),1),1,-1
3866 if(zu(k).lt.1.e-6)then
3871 kb_adj=max(2,kb_adj)
3876 ELSEIF(draft == "SH2") then
3877 tunning =min(0.8, (p(kpbli)-p(kb_adj))/(p(kt)-p(kb_adj))) !=.6
3878 tunning =max(0.2, tunning)
3879 beta = 2.5 !2.5 ! max(2.5,2./tunning)
3880 alpha= (tunning*(beta -2.)+1.)/(1.-tunning)
3881 #if ( ! defined NO_GAMMA_SUPPORT )
3882 fzu = gamma(alpha + beta)/(gamma(alpha)*gamma(beta))
3884 do k=kb_adj,min(kte,kt)
3885 kratio= (p(k)-p(kb_adj))/(p(kt)-p(kb_adj)) !float(k)/float(kt+1)
3886 zu(k) = zubeg+FZU*kratio**(alpha-1.0) * (1.0-kratio)**(beta-1.0)
3888 if(maxval(zu(kts:min(ktf,kt+1))).gt.0.) &
3889 zu(kts:min(ktf,kt+1))= zu(kts:min(ktf,kt+1))/maxval(zu(kts:min(ktf,kt+1)))
3890 do k=maxloc(zu(:),1),1,-1
3891 if(zu(k).lt.1.e-6)then
3897 ELSEIF(draft == "SH3") then
3900 alpha = tunning*beta
3901 beta = 3.5 ! max(2.5,2./tunning)
3902 alpha = beta -2. ! +1 !max(1.1,tunning*beta-abs(1.5-tunning)*5.)
3905 kratio= float(k)/float(40)
3906 zuh2(k) = zubeg+FZU*kratio**(alpha-1.0) * (1.0-kratio)**(beta-1.0)
3908 if(maxval(zuh2(1:40)).gt.0.) &
3909 zuh2(:)= zuh2(:)/ maxval(zuh2(1:40))
3910 ELSEIF(draft == "MID") then
3912 tunning=p(kt)+(p(kb_adj)-p(kt))*.9 !*.33
3913 tunning =min(0.9, (tunning-p(kb_adj))/(p(kt)-p(kb_adj))) !=.6
3914 tunning =max(0.2, tunning)
3915 beta = 1.3 !2.5 ! max(2.5,2./tunning)
3916 alpha= (tunning*(beta -2.)+1.)/(1.-tunning)
3917 #if ( ! defined NO_GAMMA_SUPPORT )
3918 fzu = gamma(alpha + beta)/(gamma(alpha)*gamma(beta))
3920 do k=kb_adj,min(kte,kt)
3921 kratio= (p(k)-p(kb_adj))/(p(kt)-p(kb_adj)) !float(k)/float(kt+1)
3922 zu(k) = zubeg+FZU*kratio**(alpha-1.0) * (1.0-kratio)**(beta-1.0)
3924 if(maxval(zu(kts:min(ktf,kt+1))).gt.0.) &
3925 zu(kts:min(ktf,kt+1))= zu(kts:min(ktf,kt+1))/maxval(zu(kts:min(ktf,kt+1)))
3926 do k=maxloc(zu(:),1),1,-1
3927 if(zu(k).lt.1.e-6)then
3932 kb_adj=max(2,kb_adj)
3937 ELSEIF(draft == "DOWN" .or. draft == "DOWNM") then
3940 ! beta = 3.0/tunning
3943 ! alpha = tunning*beta
3947 tunning =min(0.9, (tunning-p(1))/(p(kt)-p(1))) !=.6
3948 tunning =max(0.2, tunning)
3949 beta = 4. !2.5 ! max(2.5,2./tunning)
3950 alpha= (tunning*(beta -2.)+1.)/(1.-tunning)
3951 #if ( ! defined NO_GAMMA_SUPPORT )
3952 fzu = gamma(alpha + beta)/(gamma(alpha)*gamma(beta))
3956 kratio= (p(k)-p(1))/(p(kt)-p(1))
3957 zu(k) = FZU*kratio**(alpha-1.0) * (1.0-kratio)**(beta-1.0)
3959 ! if(maxloc(zuh(:),1).ge.kpbli)then
3960 ! do k=maxloc(zuh(:),1),1,-1
3961 ! kk=kpbli+k-maxloc(zuh(:),1)
3962 ! if(kk.gt.1)zu(kk)=zuh(k)
3964 ! do k=maxloc(zuh(:),1)+1,kt
3965 ! kk=kpbli+k-maxloc(zuh(:),1)
3966 ! if(kk.le.kt)zu(kk)=zuh(k)
3969 ! do k=2,kt ! maxloc(zuh(:),1)
3973 fzu=maxval(zu(kts:min(ktf,kt+1)))
3975 zu(kts:min(ktf,kt+1))= zu(kts:min(ktf,kt+1))/fzu
3976 ! if(zu(2).gt.max_mass)fzu=max_mass/zu(2) ! max(0.,zu(2)-max_mass)
3985 ! if(maxval(zu(kts:min(ktf,kt+1))).gt.0.) &
3986 ! zu(kts:min(ktf,kt+1))= zu(kts:min(ktf,kt+1))/ maxval(zu(kts:min(ktf,kt+1)))
3987 END SUBROUTINE get_zu_zd_pdf_fim
3989 !-------------------------------------------------------------------------
3990 SUBROUTINE cup_up_aa1bl(aa0,t,tn,q,qo,dtime, &
3991 z,zu,dby,gamma_cup,t_cup, &
4001 ! only local wrf dimensions are need as of now in this routine
4007 ! aa0 cloud work function
4008 ! gamma_cup = gamma on model cloud levels
4009 ! t_cup = temperature (Kelvin) on model cloud levels
4010 ! dby = buoancy term
4011 ! zu= normalized updraft mass flux
4012 ! z = heights of model levels
4013 ! ierr error value, maybe modified in this routine
4015 real, dimension (its:ite,kts:kte) &
4017 z,zu,gamma_cup,t_cup,dby,t,tn,q,qo
4018 integer, dimension (its:ite) &
4021 real, intent(in) :: dtime
4027 integer, dimension (its:ite) &
4028 ,intent (inout) :: &
4030 real, dimension (its:ite) &
4034 ! local variables in this routine
4047 IF(ierr(i).ne.0 )GO TO 100
4048 IF(k.gt.KBCON(i))GO TO 100
4051 !print*,"dz=",i,k,z(i,k),Z(I,K-1),dz
4052 !da=zu(i,k)*DZ*(9.81/(1004.*( &
4053 ! (T_cup(I,K)))))*DBY(I,K-1)/ &
4054 ! (1.+GAMMA_CUP(I,K))
4055 ! IF(K.eq.KTOP(I).and.da.le.0.)go to 100
4057 dA= DZ*9.81*( tn(i,k)-t(i,k) + 0.608*(qo(i,k)-q(i,k)))/dtime
4061 END SUBROUTINE cup_up_aa1bl
4062 !----------------------------------------------------------------------
4063 SUBROUTINE get_inversion_layers(ierr,p_cup,t_cup,z_cup,qo_cup,qeso_cup,k_inv_layers,&
4064 kstart,kend,dtempdz,itf,ktf,its,ite, kts,kte)
4067 integer ,intent (in ) :: itf,ktf,its,ite,kts,kte
4068 integer, dimension (its:ite) ,intent (in ) :: ierr,kstart,kend
4069 integer, dimension (its:ite) :: kend_p3
4071 real, dimension (its:ite,kts:kte), intent (in ) :: p_cup,t_cup,z_cup,qo_cup,qeso_cup
4072 real, dimension (its:ite,kts:kte), intent (out) :: dtempdz
4073 integer, dimension (its:ite,kts:kte), intent (out) :: k_inv_layers
4075 real :: dp,l_mid,l_shal,first_deriv(kts:kte),sec_deriv(kts:kte)
4076 integer:: ken,kadd,kj,i,k,ilev,kk,ix,k800,k550,mid,shal
4078 !-initialize k_inv_layers as undef
4081 k_inv_layers(:,:) = 1
4083 if(ierr(i) == 0)then
4084 kend_p3(i)=kend(i)+3
4085 DO k = kts+1,kend_p3(i)+4
4087 first_deriv(k)= (t_cup(i,k+1)-t_cup(i,k-1))/(z_cup(i,k+1)-z_cup(i,k-1))
4088 dtempdz(i,k)=first_deriv(k)
4090 DO k = kts+2,kend_p3(i)+3
4092 sec_deriv(k)= (first_deriv(k+1)-first_deriv(k-1))/(z_cup(i,k+1)-z_cup(i,k-1))
4093 sec_deriv(k)= abs(sec_deriv(k))
4096 ilev=max(kts+2,kstart(i)+1)
4099 DO WHILE (ilev < kend_p3(i)) !(z_cup(i,ilev)<15000.)
4100 do kk=k,kend_p3(i)+2 !k,ktf-2
4102 if(sec_deriv(kk) < sec_deriv(kk+1) .and. &
4103 sec_deriv(kk) < sec_deriv(kk-1) ) then
4104 k_inv_layers(i,ix)=kk
4115 ken=maxloc(k_inv_layers(i,:),1)
4117 kk=k_inv_layers(i,k+kadd)
4120 if( dtempdz(i,kk) < dtempdz(i,kk-1) .and. &
4121 dtempdz(i,kk) < dtempdz(i,kk+1) ) then ! the layer is not a local maximum
4124 if(k_inv_layers(i,kj+kadd).gt.1)k_inv_layers(i,kj) = k_inv_layers(i,kj+kadd)
4125 if(k_inv_layers(i,kj+kadd).eq.1)k_inv_layers(i,kj) = 1
4132 !- find the locations of inversions around 800 and 550 hPa
4135 if(ierr(i) /= 0) cycle
4137 !- now find the closest layers of 800 and 550 hPa.
4138 do k=1,maxloc(k_inv_layers(i,:),1) !kts,kte !kstart(i),kend(i) !kts,kte
4139 dp=p_cup(i,k_inv_layers(i,k))-p_cup(i,kstart(i))
4140 sec_deriv(k)=abs(dp)-l_shal
4142 k800=minloc(abs(sec_deriv),1)
4145 do k=1,maxloc(k_inv_layers(i,:),1) !kts,kte !kstart(i),kend(i) !kts,kte
4146 dp=p_cup(i,k_inv_layers(i,k))-p_cup(i,kstart(i))
4147 sec_deriv(k)=abs(dp)-l_mid
4149 k550=minloc(abs(sec_deriv),1)
4150 !-save k800 and k550 in k_inv_layers array
4153 k_inv_layers(i,shal)=k_inv_layers(i,k800) ! this is for shallow convection
4154 k_inv_layers(i,mid )=k_inv_layers(i,k550) ! this is for mid/congestus convection
4155 k_inv_layers(i,mid+1:kte)=-1
4159 END SUBROUTINE get_inversion_layers
4160 !-----------------------------------------------------------------------------------
4161 FUNCTION DERIV3(xx, xi, yi, ni, m)
4162 !============================================================================*/
4163 ! Evaluate first- or second-order derivatives
4164 ! using three-point Lagrange interpolation
4165 ! written by: Alex Godunov (October 2009)
4167 ! xx - the abscissa at which the interpolation is to be evaluated
4168 ! xi() - the arrays of data abscissas
4169 ! yi() - the arrays of data ordinates
4170 ! ni - size of the arrays xi() and yi()
4171 ! m - order of a derivative (1 or 2)
4173 ! deriv3 - interpolated value
4174 !============================================================================*/
4177 integer, parameter :: n=3
4178 integer ni, m,i, j, k, ix
4180 real:: xi(ni), yi(ni), x(n), f(n)
4182 ! exit if too high-order derivative was needed,
4188 ! if x is ouside the xi(1)-xi(ni) interval set deriv3=0.0
4189 if (xx < xi(1) .or. xx > xi(ni)) then
4191 stop "problems with finding the 2nd derivative"
4194 ! a binary (bisectional) search to find i so that xi(i-1) < x < xi(i)
4199 if (xx < xi(k)) then
4206 ! shift i that will correspond to n-th order of interpolation
4207 ! the search point will be in the middle in x_i, x_i+1, x_i+2 ...
4210 ! check boundaries: if i is ouside of the range [1, ... n] -> shift i
4212 if (i + n > ni) i=ni-n+1
4214 ! old output to test i
4215 ! write(*,100) xx, i
4216 ! 100 format (f10.5, I5)
4218 ! just wanted to use index i
4220 ! initialization of f(n) and x(n)
4226 ! calculate the first-order derivative using Lagrange interpolation
4228 deriv3 = (2.0*xx - (x(2)+x(3)))*f(1)/((x(1)-x(2))*(x(1)-x(3)))
4229 deriv3 = deriv3 + (2.0*xx - (x(1)+x(3)))*f(2)/((x(2)-x(1))*(x(2)-x(3)))
4230 deriv3 = deriv3 + (2.0*xx - (x(1)+x(2)))*f(3)/((x(3)-x(1))*(x(3)-x(2)))
4231 ! calculate the second-order derivative using Lagrange interpolation
4233 deriv3 = 2.0*f(1)/((x(1)-x(2))*(x(1)-x(3)))
4234 deriv3 = deriv3 + 2.0*f(2)/((x(2)-x(1))*(x(2)-x(3)))
4235 deriv3 = deriv3 + 2.0*f(3)/((x(3)-x(1))*(x(3)-x(2)))
4238 !=============================================================================================
4239 SUBROUTINE get_lateral_massflux(itf,ktf, its,ite, kts,kte &
4240 ,ierr,ktop,zo_cup,zuo,cd,entr_rate_2d &
4241 ,up_massentro, up_massdetro ,up_massentr, up_massdetr &
4242 ,draft,kbcon,k22,up_massentru,up_massdetru,lambau)
4245 character *(*), intent (in) :: draft
4246 integer, intent(in):: itf,ktf, its,ite, kts,kte
4247 integer, intent(in) , dimension(its:ite) :: ierr,ktop,kbcon,k22
4248 real, intent(in), OPTIONAL , dimension(its:ite):: lambau
4249 real, intent(in) , dimension(its:ite,kts:kte) :: zo_cup,zuo
4250 real, intent(inout), dimension(its:ite,kts:kte) :: cd,entr_rate_2d
4251 real, intent( out), dimension(its:ite,kts:kte) :: up_massentro, up_massdetro &
4252 ,up_massentr, up_massdetr
4253 real, intent( out), dimension(its:ite,kts:kte), OPTIONAL :: &
4254 up_massentru,up_massdetru
4256 Integer :: i,k, incr1,incr2
4257 REAL :: dz,trash,trash2
4261 up_massentro(i,k)=0.
4262 up_massdetro(i,k)=0.
4263 up_massentr (i,k)=0.
4264 up_massdetr (i,k)=0.
4267 if(present(up_massentru) .and. present(up_massdetru))then
4270 up_massentru(i,k)=0.
4271 up_massdetru(i,k)=0.
4276 if(ierr(i).eq.0)then
4278 do k=max(2,k22(i)+1),maxloc(zuo(i,:),1)
4279 !=> below maximum value zu -> change entrainment
4280 dz=zo_cup(i,k)-zo_cup(i,k-1)
4282 up_massdetro(i,k-1)=cd(i,k-1)*dz*zuo(i,k-1)
4283 up_massentro(i,k-1)=zuo(i,k)-zuo(i,k-1)+up_massdetro(i,k-1)
4284 if(up_massentro(i,k-1).lt.0.)then
4285 up_massentro(i,k-1)=0.
4286 up_massdetro(i,k-1)=zuo(i,k-1)-zuo(i,k)
4287 if(zuo(i,k-1).gt.0.)cd(i,k-1)=up_massdetro(i,k-1)/(dz*zuo(i,k-1))
4289 if(zuo(i,k-1).gt.0.)entr_rate_2d(i,k-1)=(up_massentro(i,k-1))/(dz*zuo(i,k-1))
4291 do k=maxloc(zuo(i,:),1)+1,ktop(i)
4292 !=> above maximum value zu -> change detrainment
4293 dz=zo_cup(i,k)-zo_cup(i,k-1)
4294 up_massentro(i,k-1)=entr_rate_2d(i,k-1)*dz*zuo(i,k-1)
4295 up_massdetro(i,k-1)=zuo(i,k-1)+up_massentro(i,k-1)-zuo(i,k)
4296 if(up_massdetro(i,k-1).lt.0.)then
4297 up_massdetro(i,k-1)=0.
4298 up_massentro(i,k-1)=zuo(i,k)-zuo(i,k-1)
4299 if(zuo(i,k-1).gt.0.)entr_rate_2d(i,k-1)=(up_massentro(i,k-1))/(dz*zuo(i,k-1))
4302 if(zuo(i,k-1).gt.0.)cd(i,k-1)=up_massdetro(i,k-1)/(dz*zuo(i,k-1))
4304 up_massdetro(i,ktop(i))=zuo(i,ktop(i))
4305 up_massentro(i,ktop(i))=0.
4308 entr_rate_2d(i,k)=0.
4309 up_massentro(i,k)=0.
4310 up_massdetro(i,k)=0.
4313 up_massentr (i,k-1)=up_massentro(i,k-1)
4314 up_massdetr (i,k-1)=up_massdetro(i,k-1)
4316 if(present(up_massentru) .and. present(up_massdetru))then
4318 up_massentru(i,k-1)=up_massentro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
4319 up_massdetru(i,k-1)=up_massdetro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
4325 do k=k22(i)+1,ktop(i)
4326 trash2=trash2+entr_rate_2d(i,k)
4328 do k=k22(i)+1,kbcon(i)
4329 trash=trash+entr_rate_2d(i,k)
4334 END SUBROUTINE get_lateral_massflux
4335 !==============================================================================
4336 !----------------------------------------------------------------------
4337 SUBROUTINE gfinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
4340 P_QC,P_QI,P_FIRST_SCALAR, &
4343 ids, ide, jds, jde, kds, kde, &
4344 ims, ime, jms, jme, kms, kme, &
4345 its, ite, jts, jte, kts, kte )
4346 !--------------------------------------------------------------------
4348 !--------------------------------------------------------------------
4349 LOGICAL , INTENT(IN) :: restart,allowed_to_read
4350 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
4351 ims, ime, jms, jme, kms, kme, &
4352 its, ite, jts, jte, kts, kte
4353 INTEGER , INTENT(IN) :: P_FIRST_SCALAR, P_QI, P_QC
4355 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
4361 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
4365 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
4369 INTEGER :: i, j, k, itf, jtf, ktf
4374 IF(.not.restart)THEN
4396 IF (P_QC .ge. P_FIRST_SCALAR) THEN
4406 IF (P_QI .ge. P_FIRST_SCALAR) THEN
4418 END SUBROUTINE gfinit
4419 END MODULE module_cu_gf_deep