| blob | a02f5e6bc18a62833a64e5151a212f16b161e0e3 |
1 MODULE module_cu_gf_deep
3 #if ( WRF_CHEM == 1)
4 USE module_cu_gf_ctrans,only: ctrans_gf
5 #endif
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.
26 !
27 ! aerosol awareness, do not user yet!
28 !
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
36 contains
39 SUBROUTINE CUP_gf( &
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
46 ,DTIME &
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)
51 ,xland & ! land mask
53 ,zo & ! heights above surface
54 ,forcing & ! only diagnostic
55 ,T & ! T before forcing
56 ,Q & ! Q before forcing
57 ,Z1 & ! terrain
58 ,Tn & ! T including forcing
59 ,QO & ! Q including forcing
60 ,PO & ! pressure (mb)
61 ,PSUR & ! surface pressure (mb)
62 ,US & ! u on mass points
63 ,VS & ! v on mass points
64 ,rho & ! density
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
75 ,edto &
76 ,xmb_out & !the xmb's may be needed for dicycle
77 ,xmbm_in &
78 ,xmbs_in &
79 ,pre &
80 ,outu & ! momentum tendencies at mass points
81 ,outv &
82 ,outt & ! temperature tendencies
83 ,outq & ! q tendencies
84 ,outqc & ! ql/qice tendencies
85 ,kbcon &
86 ,ktop &
87 ,cupclw & ! used for direct coupling to radiation, but with tuning factors
88 ,ierr & ! ierr flags are error flags, used for debugging
89 ,ierrc &
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
97 ! 3 = closures
98 ! more is possible, talk to developer or
99 ! implement yourself. pattern is expected to be
100 ! betwee -1 and +1
101 #if ( WRF_DFI_RADAR == 1 )
102 ,do_capsuppress,cap_suppress_j &
103 #endif
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 &
109 ,chem_conv_tr &
110 #endif
111 ,k22 &
112 ,jmin)
114 IMPLICIT NONE
116 integer &
117 ,intent (in ) :: &
118 nranflag,itf,ktf,its,ite, kts,kte,ipr,imid
119 integer, intent (in ) :: &
120 ichoice
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 )
127 !
128 ! option of cap suppress:
129 ! do_capsuppress = 1 do
130 ! do_capsuppress = other don't
131 !
132 !
133 INTEGER, INTENT(IN ) ,OPTIONAL :: do_capsuppress
134 REAL, DIMENSION( its:ite ) :: cap_suppress_j
135 #endif
136 !
137 !
138 !
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 ) :: &
149 pre,xmb_out
150 real, dimension (its:ite) &
151 ,intent (in ) :: &
152 hfx,qfx,xmbm_in,xmbs_in
153 integer, dimension (its:ite) &
154 ,intent (inout ) :: &
155 kbcon,ktop
156 integer, dimension (its:ite) &
157 ,intent (in ) :: &
158 kpbl
159 !
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
163 !
164 real, dimension (its:ite,kts:kte) &
165 ,intent (in ) :: &
166 dhdt,rho,T,PO,US,VS,tn
167 real, dimension (its:ite,kts:kte) &
168 ,intent (inout ) :: &
169 omeg
170 real, dimension (its:ite,kts:kte) &
171 ,intent (inout) :: &
172 Q,QO,zuo,zdo
173 real, dimension (its:ite) &
174 ,intent (in ) :: &
175 dx,ccn,Z1,PSUR,xland
176 real, dimension (its:ite) &
177 ,intent (inout ) :: &
178 mconv
181 real &
182 ,intent (in ) :: &
183 dtime
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):: &
190 tracer2d
191 REAL,DIMENSION(its:ite , kts:kte , num_chem),INTENT(IN):: &
192 chem2d
193 REAL,DIMENSION(its:ite , kts:kte , num_tracer),INTENT(INOUT):: &
194 outtracert
195 REAL,DIMENSION(its:ite , kts:kte , num_tracer),INTENT(INOUT):: &
196 outchemt
197 INTEGER :: nv
198 real,dimension(its:ite,kts:kte) :: tempco
200 #endif
202 !
203 ! local ensemble dependent variables in this routine
204 !
205 real, dimension (its:ite,1) :: &
206 xaa0_ens
207 real, dimension (its:ite,1) :: &
208 edtc
209 real, dimension (its:ite,kts:kte,1) :: &
210 dellat_ens,dellaqc_ens,dellaq_ens,pwo_ens
211 !
212 !
213 !
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
220 !
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
236 !
237 !
238 ! hcd = moist static energy in downdraft
239 ! zd normalized downdraft mass flux
240 ! dby = buoancy term
241 ! entr = entrainment rate
242 ! zd = downdraft normalized mass flux
243 ! entr= entrainment rate
244 ! hcd = h in model cloud
245 ! bu = buoancy term
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!)
271 ! dby = buoancy term
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, &
285 dbyt,xdby,xhc,xzu, &
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
297 ! edt = epsilon
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)
301 ! edt = epsilon
303 real, dimension (its:ite) :: &
304 edt,edto,AA1,AA0,XAA0,HKB, &
305 HKBO,XHKB, &
306 XMB,PWAVO, &
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
317 integer :: &
318 iloop,nens3,ki,kk,I,K
319 real :: &
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
353 ! rainevap from sas
354 real zuh2(40)
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
359 flux_tun(:)=fluxtune
360 ! if(imid.eq.1)flux_tun(:)=fluxtune+.5
361 pmin=150.
362 if(imid.eq.1)pmin=75.
363 ktopdby(:)=0
364 elocp=xlv/cp
365 el2orc=xlv*xlv/(r_v*cp)
366 evfact=.3
367 evfactl=.3
368 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
369 !
370 !proportionality constant to estimate pressure gradient of updraft (Zhang and Wu, 2003, JAS
371 !
372 ! ECMWF
373 pgcon=0.
374 lambau(:)=2.
375 ! here random must be between -1 and 1
376 if(nranflag == 1)then
377 lambau(:)=1.5+rand_mom(:)
378 endif
379 ! SAS
380 ! lambau=0.
381 ! pgcon=-.55
382 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
383 ztexec(:) = 0.
384 zqexec(:) = 0.
385 zws(:) = 0.
387 do i=its,itf
388 !- buoyancy flux (H+LE)
389 buo_flux= (hfx(i)/cp+0.608*t(i,1)*qfx(i)/xlv)/rho(i,1)
390 pgeoh = zo(i,2)*g
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)
398 !- moisture excess
399 zqexec(i) = MAX(flux_tun(i)*qfx(i)/xlv/(rho(i,1)*zws(i)),0.)
400 endif
401 !- zws for shallow convection closure (Grant 2001)
402 !- height of the pbl
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
406 enddo
407 ! cap_maxs=225.
408 ! if(imid.eq.1)cap_maxs=150.
409 cap_maxs=75. ! 150.
410 ! if(imid.eq.1)cap_maxs=100.
411 do i=its,itf
412 edto(i)=0.
413 closure_n(i)=16.
414 xmb_out(i)=0.
415 cap_max(i)=cap_maxs
416 cap_max_increment(i)=20.
417 if(imid.eq.1)cap_max_increment(i)=10.
418 !
419 ! for water or ice
420 !
421 xland1(i)=int(xland(i)+.0001) ! 1.
422 if(xland(i).gt.1.5 .or. xland(i).lt.0.5)then
423 xland1(i)=0
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.
427 else
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.
430 endif
431 ierrc(i)=" "
432 ! cap_max_increment(i)=1.
433 enddo
434 if(use_excess == 0 )then
435 ztexec(:)=0
436 zqexec(:)=0
437 endif
438 #if ( WRF_DFI_RADAR == 1 )
439 if(do_capsuppress == 1) then
440 do i=its,itf
441 cap_max(i)=cap_maxs
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
445 cap_max(i)=10.0
446 endif
447 enddo
448 endif
449 #endif
451 !
452 !--- initial entrainment rate (these may be changed later on in the
453 !--- program
454 !
455 start_level(:)=kte
456 do i=its,ite
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
465 frh=frh_thresh
466 radius=sqrt(frh*dx(i)*dx(i)/3.14)
467 entr_rate(i)=.2/radius
468 endif
469 sig(i)=(1.-frh)**2
470 enddo
471 sig_thresh = (1.-frh_thresh)**2
474 !
475 !--- entrainment of mass
476 !
477 !
478 !--- initial detrainmentrates
479 !
480 do k=kts,ktf
481 do i=its,itf
482 cnvwt(i,k)=0.
483 zuo(i,k)=0.
484 zdo(i,k)=0.
485 z(i,k)=zo(i,k)
486 xz(i,k)=zo(i,k)
487 cupclw(i,k)=0.
488 cd(i,k)=1.e-9 ! 1.*entr_rate
489 ! if(imid.eq.1)cd(i,k)=entr_rate(i)
490 cdd(i,k)=1.e-9
491 hcdo(i,k)=0.
492 qrcdo(i,k)=0.
493 dellaqc(i,k)=0.
494 enddo
495 enddo
496 !
497 !--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
498 ! base mass flux
499 !
500 edtmax(:)=1.
501 if(imid.eq.1)edtmax(:)=.15
502 edtmin(:)=.1
503 if(imid.eq.1)edtmin(:)=.05
504 !
505 !--- minimum depth (m), clouds must have
506 !
507 depth_min=1000.
508 if(imid.eq.1)depth_min=500.
509 !
510 !--- maximum depth (mb) of capping
511 !--- inversion (larger cap = no convection)
512 !
513 DO i=its,itf
514 ! if(imid.eq.0)then
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)
517 ! endif
518 kbmax(i)=1
519 aa0(i)=0.
520 aa1(i)=0.
521 edt(i)=0.
522 kstabm(i)=ktf-1
523 IERR2(i)=0
524 IERR3(i)=0
525 x_add=0.
526 enddo
527 ! do i=its,itf
528 ! cap_max(i)=cap_maxs
529 ! cap_max3(i)=25.
531 ! enddo
532 !
533 !--- max height(m) above ground where updraft air can originate
534 !
535 zkbmax=4000.
536 if(imid.eq.1)zkbmax=2000.
537 !
538 !--- height(m) above which no downdrafts are allowed to originate
539 !
540 zcutdown=4000.
541 !
542 !--- depth(m) over which downdraft detrains all its mass
543 !
544 z_detr=1000.
545 ! if(imid.eq.1)z_detr=800.
546 !
548 !
549 !--- environmental conditions, FIRST HEIGHTS
550 !
551 do i=its,itf
552 do k=1,maxens3
553 xf_ens(i,k)=0.
554 pr_ens(i,k)=0.
555 enddo
556 enddo
558 !
559 !--- calculate moist static energy, heights, qes
560 !
561 call cup_env(z,qes,he,hes,t,q,po,z1, &
562 psur,ierr,tcrit,-1, &
563 itf,ktf, &
564 its,ite, kts,kte)
565 call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
566 psur,ierr,tcrit,-1, &
567 itf,ktf, &
568 its,ite, kts,kte)
570 !
571 !--- environmental values on cloud levels
572 !
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, &
575 ierr,z1, &
576 itf,ktf, &
577 its,ite, kts,kte)
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, &
580 ierr,z1, &
581 itf,ktf, &
582 its,ite, kts,kte)
583 do i=its,itf
584 if(ierr(i).eq.0)then
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)
588 do k=kts+1,ktf
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))
591 enddo
592 endif
593 enddo
594 do i=its,itf
595 if(ierr(i).eq.0)then
596 do k=kts,ktf
597 if(zo_cup(i,k).gt.zkbmax+z1(i))then
598 kbmax(i)=k
599 go to 25
600 endif
601 enddo
602 25 continue
603 !
604 !--- level where detrainment for downdraft starts
605 !
606 do k=kts,ktf
607 if(zo_cup(i,k).gt.z_detr+z1(i))then
608 kdet(i)=k
609 go to 26
610 endif
611 enddo
612 26 continue
613 !
614 endif
615 enddo
616 !
617 !
618 !
619 !------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
620 !
621 start_k22=2
622 DO 36 i=its,itf
623 IF(ierr(I).eq.0)THEN
624 k22(i)=maxloc(HEO_CUP(i,start_k22:kbmax(i)+2),1)+start_k22-1
625 if(K22(I).GE.KBMAX(i))then
626 ierr(i)=2
627 ierrc(i)="could not find k22"
628 ktop(i)=0
629 k22(i)=0
630 kbcon(i)=0
631 endif
632 endif
633 36 CONTINUE
634 !
635 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
636 !
638 do i=its,itf
639 IF(ierr(I).eq.0)THEN
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)
643 endif ! ierr
644 enddo
645 jprnt=0
646 iloop=1
647 if(imid.eq.1)iloop=5
648 call cup_kbcon(ierrc,cap_max_increment,iloop,k22,kbcon,heo_cup,heso_cup, &
649 hkbo,ierr,kbmax,po_cup,cap_max, &
650 ztexec,zqexec, &
651 jprnt,itf,ktf, &
652 its,ite, kts,kte, &
653 z_cup,entr_rate,heo,imid)
654 !
655 !--- increase detrainment in stable layers
656 !
657 CALL cup_minimi(HEso_cup,Kbcon,kstabm,kstabi,ierr, &
658 itf,ktf, &
659 its,ite, kts,kte)
660 DO i=its,itf
661 IF(ierr(I) == 0)THEN
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
664 ierr(i)=231
665 cycle
666 endif
667 !
668 ! never go too low...
669 !
670 ! if(imid.eq.0 .and. xland1(i).eq.0)x_add=150.
671 x_add=0.
672 do k=kbcon(i)+1,ktf
673 if(po(i,kbcon(i))-po(i,k) > pmin+x_add)then
674 pmin_lev(i)=k
675 exit
676 endif
677 enddo
678 !
679 ! initial conditions for updraft
680 !
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)
684 endif
685 enddo
686 !
687 !--- get inversion layers for mid level cloud tops
688 !
689 if(imid.eq.1)then
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)
692 endif
693 DO i=its,itf
694 if(kstabi(i).lt.kbcon(i))then
695 kbcon(i)=1
696 ierr(i)=42
697 endif
698 do k=kts,ktf
699 entr_rate_2d(i,k)=entr_rate(i)
700 enddo
701 IF(ierr(I).eq.0)THEN
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))
704 do k=kts,ktf
705 frh = min(qo_cup(i,k)/qeso_cup(i,k),1.)
706 entr_rate_2d(i,k)=entr_rate(i) *(1.3-frh)
707 enddo
708 if(imid.eq.1)then
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))
713 ktopdby(i)=ktop(i)
714 else
715 do k=kbcon(i)+1,ktf
716 if((po_cup(i,k22(i))-po_cup(i,k)).gt.500.)then
717 ktop(i)=k
718 ktopdby(i)=ktop(i)
719 exit
720 endif
721 enddo
722 endif ! k_inv_lay
723 endif
725 endif
726 ENDDO
727 !
728 !-- get normalized mass flux, entrainment and detrainmentrates for updraft
729 !
730 i=0
731 !- for mid level clouds we do not allow clouds taller than where stability
732 !- changes
733 if(imid.eq.1)then
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)
736 else
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)
739 endif
740 !
741 !
742 !
743 do i=its,itf
744 if(ierr(i).eq.0)then
746 if(k22(i).gt.1)then
747 do k=1,k22(i) -1
748 zuo(i,k)=0.
749 zu (i,k)=0.
750 xzu(i,k)=0.
751 enddo
752 endif
753 do k=k22(i),ktop(i)
754 xzu(i,k)= zuo(i,k)
755 zu (i,k)= zuo(i,k)
756 enddo
757 do k=ktop(i)+1,kte
758 zuo(i,k)=0.
759 zu (i,k)=0.
760 xzu(i,k)=0.
761 enddo
762 endif
763 enddo
764 !
765 ! calculate mass entrainment and detrainment
766 !
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)
773 !
774 ! NOTE: Ktop here already includes overshooting, ktopdby is without
775 ! overshooting
776 !
777 do k=kts,ktf
778 do i=its,itf
779 uc (i,k)=0.
780 vc (i,k)=0.
781 hc (i,k)=0.
782 dby (i,k)=0.
783 hco (i,k)=0.
784 dbyo(i,k)=0.
785 enddo
786 enddo
787 do i=its,itf
788 IF(ierr(I).eq.0)THEN
789 do k=1,start_level(i)
790 uc(i,k)=u_cup(i,k)
791 vc(i,k)=v_cup(i,k)
792 enddo
793 do k=1,start_level(i)-1
794 hc (i,k)=he_cup(i,k)
795 hco(i,k)=heo_cup(i,k)
796 enddo
797 k=start_level(i)
798 hc (i,k)=hkb(i)
799 hco(i,k)=hkbo(i)
800 ENDIF
801 enddo
803 DO i=its,itf
805 ktopkeep(i)=0
806 dbyt(i,:)=0.
807 if(ierr(i) /= 0) cycle
808 ktopkeep(i)=ktop(i)
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
813 ierr(i)=51
814 exit
815 endif
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
835 ENDDO
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
840 ! if(kk.lt.ki+3)then
841 ! ierr(i)=423
842 ! endif
843 !
844 do k=ktop(i)-1,kbcon(i),-1
845 if(dbyo(i,k).gt.0.)then
846 ktopkeep(i)=k+1
847 exit
848 endif
849 enddo
850 ktop(I)=ktopkeep(i)
851 if(ierr(i).eq.0)ktop(I)=ktopkeep(i)
852 ENDDO
853 41 continue
854 DO i=its,itf
855 if(ierr(i) /= 0) cycle
856 do k=ktop(i)+1,ktf
857 HC(i,K)=hes_cup(i,k)
858 UC(i,K)=u_cup(i,k)
859 VC(i,K)=v_cup(i,k)
860 HCo(i,K)=heso_cup(i,k)
861 DBY(I,K)=0.
862 DBYo(I,K)=0.
863 zu(i,k)=0.
864 zuo(i,k)=0.
865 cd(i,k)=0.
866 entr_rate_2d(i,k)=0.
867 up_massentr(i,k)=0.
868 up_massdetr(i,k)=0.
869 up_massentro(i,k)=0.
870 up_massdetro(i,k)=0.
871 enddo
872 ENDDO
873 !
874 DO i=its,itf
875 if(ierr(i)/=0)cycle
876 if(ktop(i).lt.kbcon(i)+2)then
877 ierr(i)=5
878 ierrc(i)='ktop too small deep'
879 ktop(i)=0
880 endif
881 ENDDO
882 DO 37 i=its,itf
883 kzdown(i)=0
884 if(ierr(i).eq.0)then
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))
888 do k=kts,ktf
889 if(zo_cup(i,k).gt.zktop)then
890 kzdown(i)=k
891 kzdown(i)=min(kzdown(i),kstabi(i)-1) !
892 go to 37
893 endif
894 enddo
895 endif
896 37 CONTINUE
897 !
898 !--- DOWNDRAFT ORIGINATING LEVEL - JMIN
899 !
900 call cup_minimi(HEso_cup,K22,kzdown,JMIN,ierr, &
901 itf,ktf, &
902 its,ite, kts,kte)
903 DO 100 i=its,itf
904 IF(ierr(I).eq.0)THEN
905 !
906 !--- check whether it would have buoyancy, if there where
907 !--- no entrainment/detrainment
908 !
909 jmini = jmin(i)
910 keep_going = .TRUE.
911 do while ( keep_going )
912 keep_going = .FALSE.
913 if ( jmini - 1 .lt. kdet(i) ) kdet(i) = jmini-1
914 if ( jmini .ge. ktop(i)-1 ) jmini = ktop(i) - 2
915 ki = jmini
916 hcdo(i,ki)=heso_cup(i,ki)
917 DZ=Zo_cup(i,Ki+1)-Zo_cup(i,Ki)
918 dh=0.
919 do k=ki-1,1,-1
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))
923 if(dh.gt.0.)then
924 jmini=jmini-1
925 if ( jmini .gt. 5 ) then
926 keep_going = .TRUE.
927 else
928 ierr(i) = 9
929 ierrc(i) = "could not find jmini9"
930 exit
931 endif
932 endif
933 enddo
934 enddo
935 jmin(i) = jmini
936 if ( jmini .le. 5 ) then
937 ierr(i)=4
938 ierrc(i) = "could not find jmini4"
939 endif
940 ENDIF
941 100 continue
942 !
943 ! - Must have at least depth_min m between cloud convective base
944 ! and cloud top.
945 !
946 do i=its,itf
947 IF(ierr(I).eq.0)THEN
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
950 ierr(i)=6
951 ierrc(i)="cloud depth very shallow"
952 endif
953 endif
954 enddo
956 !
957 !--- normalized downdraft mass flux profile,also work on bottom detrainment
958 !--- in this routine
959 !
960 do k=kts,ktf
961 do i=its,itf
962 zdo(i,k)=0.
963 cdd(i,k)=0.
964 dd_massentro(i,k)=0.
965 dd_massdetro(i,k)=0.
966 dd_massentru(i,k)=0.
967 dd_massdetru(i,k)=0.
968 hcdo(i,k)=heso_cup(i,k)
969 ucd(i,k)=u_cup(i,k)
970 vcd(i,k)=v_cup(i,k)
971 dbydo(i,k)=0.
972 mentrd_rate_2d(i,k)=entr_rate(i)
973 enddo
974 enddo
975 do i=its,itf
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
979 ! beta=.01
980 edtmax(i)=max(0.1,.4-float(csum(i))*.015) !.3)
981 endif
982 if(imid.eq.1)beta=.02
983 bud(i)=0.
984 IF(ierr(I).eq.0)then
985 cdd(i,1:jmin(i))=1.e-9
986 cdd(i,jmin(i))=0.
987 dd_massdetro(i,:)=0.
988 dd_massentro(i,:)=0.
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
991 zdo(i,jmin(i))=0.
992 jmin(i)=jmin(i)-1
993 if(zdo(i,jmin(i)) .lt.1.e-8)then
994 ierr(i)=876
995 cycle
996 endif
997 endif
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))
1008 endif
1009 if(zdo(i,ki+1).gt.0.)mentrd_rate_2d(i,ki)=dd_massentro(i,ki)/(dzo*zdo(i,ki+1))
1010 enddo
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))
1021 endif
1022 if(zdo(i,ki+1).gt.0.)cdd(i,ki)= dd_massdetro(i,ki)/(dzo*zdo(i,ki+1))
1023 enddo
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))
1033 c1d(i,k)=c1
1034 enddo
1035 if(imid.eq.1)c1d(i,:)=0.
1036 ! do k=1,jmin(i)
1037 ! c1d(i,k)=0.
1038 ! enddo
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))
1046 ! enddo
1049 ! downdraft moist static energy + moisture budget
1050 do k=2,jmin(i)+1
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)
1053 enddo
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)))
1056 do ki=jmin(i) ,1,-1
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
1075 enddo
1076 endif
1078 if(bud(i).gt.0)then
1079 ierr(i)=7
1080 ierrc(i)='downdraft is not negatively buoyant '
1081 endif
1082 enddo
1083 !
1084 !--- calculate moisture properties of downdraft
1085 !
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, &
1089 itf,ktf, &
1090 its,ite, kts,kte)
1091 !
1092 !--- calculate moisture properties of updraft
1093 !
1094 if(imid.eq.1)then
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, &
1099 1,itf,ktf, &
1100 its,ite, kts,kte)
1101 else
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, &
1106 1,itf,ktf, &
1107 its,ite, kts,kte)
1108 endif
1109 do i=its,itf
1110 if(ierr(i).eq.0)then
1111 do k=kts+1,ktop(i)
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
1115 enddo
1116 endif
1117 enddo
1118 !
1119 !--- calculate workfunctions for updrafts
1120 !
1121 call cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
1122 kbcon,ktop,ierr, &
1123 itf,ktf, &
1124 its,ite, kts,kte)
1125 call cup_up_aa0(aa1,zo,zuo,dbyo,GAMMAo_CUP,tn_cup, &
1126 kbcon,ktop,ierr, &
1127 itf,ktf, &
1128 its,ite, kts,kte)
1129 do i=its,itf
1130 if(ierr(i).eq.0)then
1131 if(aa1(i).eq.0.)then
1132 ierr(i)=17
1133 ierrc(i)="cloud work function zero"
1134 endif
1135 endif
1136 enddo
1137 !
1138 !--- diurnal cycle closure
1139 !
1140 !--- AA1 from boundary layer (bl) processes only
1141 aa1_bl (:) = 0.0
1142 xf_dicycle (:) = 0.0
1143 tau_ecmwf (:) = 0.
1144 !- way to calculate the fraction of cape consumed by shallow convection
1145 iversion=1 ! ecmwf
1146 !iversion=0 ! orig
1147 !
1148 ! Betchold et al 2008 time-scale of cape removal
1149 !
1150 ! wmean is of no meaning over land....
1151 ! still working on replacing it over water
1152 !
1153 DO i=its,itf
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
1161 endif
1162 enddo
1163 tau_bl(:) = 0.
1164 !
1165 IF(dicycle == 1) then
1166 DO i=its,itf
1168 if(ierr(i).eq.0)then
1169 if(xland1(i) == 0 ) then
1170 !- over water
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
1173 else
1174 !- over land
1175 tau_bl(i) =( zo_cup(i,ktopdby(i))- zo_cup(i,kbcon(i)) ) / wmean(i)
1176 endif
1178 endif
1179 ENDDO
1181 if(iversion == 1) then
1182 !-- version ecmwf
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, &
1188 kbcon,ktop,ierr, &
1189 itf,ktf,its,ite, kts,kte)
1191 DO i=its,itf
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
1197 aa1_bl(i) = 0.0
1198 else
1199 !- multiply aa1_bl the " time-scale" - tau_bl
1200 aa1_bl(i) = max(0.,aa1_bl(i)/t_star* tau_bl(i))
1201 endif
1202 endif
1203 ENDDO
1205 else
1207 !- version for real cloud-work function
1209 !-get the profiles modified only by bl tendencies
1210 DO i=its,itf
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)
1219 endif
1220 ENDDO
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, &
1228 ierr,z1, &
1229 itf,ktf,its,ite, kts,kte)
1230 DO i=its,itf
1231 IF(ierr(I).eq.0)THEN
1232 hkbo_bl(i)=heo_cup_bl(i,k22(i))
1233 endif ! ierr
1234 ENDDO
1235 DO k=kts,ktf
1236 do i=its,itf
1237 hco_bl (i,k)=0.
1238 DBYo_bl(i,k)=0.
1239 enddo
1240 ENDDO
1241 DO i=its,itf
1242 IF(ierr(I).eq.0)THEN
1243 do k=1,kbcon(i)-1
1244 hco_bl(i,k)=hkbo_bl(i)
1245 enddo
1246 k=kbcon(i)
1247 hco_bl (i,k)=hkbo_bl(i)
1248 DBYo_bl(i,k)=Hkbo_bl(i) - HESo_cup_bl(i,k)
1249 ENDIF
1250 ENDDO
1251 !
1252 !
1253 DO i=its,itf
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)
1260 enddo
1261 do k=ktop(i)+1,ktf
1262 hco_bl (i,k)=heso_cup_bl(i,k)
1263 dbyo_bl(i,k)=0.0
1264 enddo
1265 endif
1266 ENDDO
1268 !--- calculate workfunctions for updrafts
1269 call cup_up_aa0(aa1_bl,zo,zuo,dbyo_bl,GAMMAo_CUP_bl,tn_cup_bl, &
1270 kbcon,ktop,ierr, &
1271 itf,ktf,its,ite, kts,kte)
1273 DO i=its,itf
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))
1280 ! aa1_bl(i) = 0.0
1281 !else
1282 ! !- multiply aa1_bl the "normalized time-scale" - tau_bl/ model_timestep
1283 aa1_bl(i) = aa1_bl(i)* tau_bl(i)/ dtime
1284 !endif
1285 print*,'aa0,aa1bl=',aa0(i),aa1_bl(i),aa0(i)-aa1_bl(i),tau_bl(i)!,dtime,xland(i)
1286 endif
1287 ENDDO
1288 ENDIF
1289 ENDIF ! version of implementation
1292 axx(:)=aa1(:)
1294 !
1295 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
1296 !
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, &
1300 its,ite, kts,kte)
1301 do i=its,itf
1302 if(ierr(i).eq.0)then
1303 edto(i)=edtc(i,1)
1304 endif
1305 enddo
1306 do k=kts,ktf
1307 do i=its,itf
1308 dellat_ens (i,k,1)=0.
1309 dellaq_ens (i,k,1)=0.
1310 dellaqc_ens(i,k,1)=0.
1311 pwo_ens (i,k,1)=0.
1312 enddo
1313 enddo
1314 !
1315 !--- change per unit mass that a model cloud would modify the environment
1316 !
1317 !--- 1. in bottom layer
1318 !
1319 do k=kts,kte
1320 do i=its,itf
1321 dellu (i,k)=0.
1322 dellv (i,k)=0.
1323 dellah (i,k)=0.
1324 dellat (i,k)=0.
1325 dellaq (i,k)=0.
1326 dellaqc(i,k)=0.
1327 enddo
1328 enddo
1329 !
1330 !---------------------------------------------- cloud level ktop
1331 !
1332 !- - - - - - - - - - - - - - - - - - - - - - - - model level ktop-1
1333 ! . . .
1334 ! . . .
1335 ! . . .
1336 ! . . .
1337 ! . . .
1338 ! . . .
1339 !
1340 !---------------------------------------------- cloud level k+2
1341 !
1342 !- - - - - - - - - - - - - - - - - - - - - - - - model level k+1
1343 !
1344 !---------------------------------------------- cloud level k+1
1345 !
1346 !- - - - - - - - - - - - - - - - - - - - - - - - model level k
1347 !
1348 !---------------------------------------------- cloud level k
1349 !
1350 ! . . .
1351 ! . . .
1352 ! . . .
1353 ! . . .
1354 ! . . .
1355 ! . . .
1356 ! . . .
1357 ! . . .
1358 ! . . .
1359 ! . . .
1360 !
1361 !---------------------------------------------- cloud level 3
1362 !
1363 !- - - - - - - - - - - - - - - - - - - - - - - - model level 2
1364 !
1365 !---------------------------------------------- cloud level 2
1366 !
1367 !- - - - - - - - - - - - - - - - - - - - - - - - model level 1
1369 do i=its,itf
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
1377 do k=kts+1,ktop(i)
1378 ! these three are only used at or near mass detrainment and/or entrainment levels
1379 pgc=pgcon
1380 entupk=0.
1381 if(k == k22(i)-1) entupk=zuo(i,k+1)
1382 detupk=0.
1383 entdoj=0.
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)
1393 ! SPECIAL LEVELS
1394 if(k.eq.ktop(i))then
1395 detupk=zuo(i,ktop(i))
1396 subin=0.
1397 subdown=0.
1398 detdo=0.
1399 entdo=0.
1400 entup=0.
1401 detup=0.
1402 endif
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)
1408 endif
1409 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
1410 pgc=pgcon
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
1422 enddo ! k
1424 endif
1425 enddo
1428 do i=its,itf
1429 !trash = 0.0
1430 !trash2 = 0.0
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
1448 !- H budget
1449 !trash2 = trash2+ (dellah(i,1))*dp/g
1452 do k=kts+1,ktop(i)
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
1472 !---
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) &
1488 - C_up + E_dn
1489 !- check water conservation liq+condensed (including rainfall)
1490 ! trash= trash+ (dellaq(i,k)+dellaqc(i,k)+ G_rain-E_dn)*dp/g
1492 enddo ! k
1493 endif
1495 enddo
1496 444 format(1x,i2,1x,7e12.4) !,1x,f7.2,2x,e13.5)
1497 !
1498 !--- using dellas, calculate changed environmental profiles
1499 !
1500 mbdt=.1
1501 do i=its,itf
1502 xaa0_ens(i,1)=0.
1503 enddo
1505 do i=its,itf
1506 if(ierr(i).eq.0)then
1507 do k=kts,ktf
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))
1515 enddo
1516 ENDIF
1517 enddo
1518 do i=its,itf
1519 if(ierr(i).eq.0)then
1520 XHE(I,ktf)=HEO(I,ktf)
1521 XQ(I,ktf)=QO(I,ktf)
1522 XT(I,ktf)=TN(I,ktf)
1523 endif
1524 enddo
1525 !
1526 !--- calculate moist static energy, heights, qes
1527 !
1528 call cup_env(xz,xqes,xhe,xhes,xt,xq,po,z1, &
1529 psur,ierr,tcrit,-1, &
1530 itf,ktf, &
1531 its,ite, kts,kte)
1532 !
1533 !--- environmental values on cloud levels
1534 !
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, &
1537 ierr,z1, &
1538 itf,ktf, &
1539 its,ite, kts,kte)
1540 !
1541 !
1542 !**************************** static control
1543 !
1544 !--- moist static energy inside cloud
1545 !
1546 do k=kts,ktf
1547 do i=its,itf
1548 xhc(i,k)=0.
1549 xDBY(I,K)=0.
1550 enddo
1551 enddo
1552 do i=its,itf
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)
1558 enddo
1559 k=start_level(i)
1560 xhc(i,k)=xhkb(i)
1561 endif !ierr
1562 enddo
1563 !
1564 !
1565 do i=its,itf
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)
1572 enddo
1573 do k=ktop(i)+1,ktf
1574 xHC (i,K)=xhes_cup(i,k)
1575 xDBY(I,K)=0.
1576 enddo
1577 endif
1578 enddo
1580 !
1581 !--- workfunctions for updraft
1582 !
1583 call cup_up_aa0(xaa0,xz,xzu,xdby,GAMMA_CUP,xt_cup, &
1584 kbcon,ktop,ierr, &
1585 itf,ktf, &
1586 its,ite, kts,kte)
1587 do i=its,itf
1588 if(ierr(i).eq.0)then
1589 xaa0_ens(i,1)=xaa0(i)
1590 do k=kts,ktop(i)
1591 do nens3=1,maxens3
1592 if(nens3.eq.7)then
1593 !--- b=0
1594 pr_ens(i,nens3)=pr_ens(i,nens3) &
1595 +pwo(i,k)+edto(i)*pwdo(i,k)
1596 !--- b=beta
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)
1600 !--- b=beta/2
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)
1604 else
1605 pr_ens(i,nens3)=pr_ens(i,nens3)+ &
1606 pwo(i,k) +edto(i)*pwdo(i,k)
1607 endif
1608 enddo
1609 enddo
1610 if(pr_ens(i,7).lt.1.e-6)then
1611 ierr(i)=18
1612 ierrc(i)="total normalized condensate too small"
1613 do nens3=1,maxens3
1614 pr_ens(i,nens3)=0.
1615 enddo
1616 endif
1617 do nens3=1,maxens3
1618 if(pr_ens(i,nens3).lt.1.e-5)then
1619 pr_ens(i,nens3)=0.
1620 endif
1621 enddo
1622 endif
1623 enddo
1624 200 continue
1625 !
1626 !--- LARGE SCALE FORCING
1627 !
1628 !
1629 !------- CHECK wether aa0 should have been zero, assuming this
1630 ! ensemble is chosen
1631 !
1632 !
1633 do i=its,itf
1634 ierr2(i)=ierr(i)
1635 ierr3(i)=ierr(i)
1636 k22x(i)=k22(i)
1637 enddo
1638 CALL cup_MAXIMI(HEO_CUP,2,KBMAX,K22x,ierr, &
1639 itf,ktf, &
1640 its,ite, kts,kte)
1641 iloop=2
1642 call cup_kbcon(ierrc,cap_max_increment,iloop,k22x,kbconx,heo_cup, &
1643 heso_cup,hkbo,ierr2,kbmax,po_cup,cap_max, &
1644 ztexec,zqexec, &
1645 0,itf,ktf, &
1646 its,ite, kts,kte, &
1647 z_cup,entr_rate,heo,imid)
1648 iloop=3
1649 call cup_kbcon(ierrc,cap_max_increment,iloop,k22x,kbconx,heo_cup, &
1650 heso_cup,hkbo,ierr3,kbmax,po_cup,cap_max, &
1651 ztexec,zqexec, &
1652 0,itf,ktf, &
1653 its,ite, kts,kte, &
1654 z_cup,entr_rate,heo,imid)
1655 !
1656 !--- calculate cloud base mass flux
1657 !
1659 DO I = its,itf
1660 mconv(i) = 0
1661 if(ierr(i)/=0)cycle
1662 DO K=1,ktop(i)
1663 dq=(qo_cup(i,k+1)-qo_cup(i,k))
1664 mconv(i)=mconv(i)+omeg(i,k)*dq/g
1665 ENDDO
1666 ENDDO
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, &
1671 ichoice, &
1672 imid,ipr,itf,ktf, &
1673 its,ite, kts,kte, &
1674 dicycle,tau_ecmwf,aa1_bl,xf_dicycle)
1675 !
1676 do k=kts,ktf
1677 do i=its,itf
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)
1683 else
1684 dellat_ens (i,k,1)=0.
1685 dellaq_ens (i,k,1)=0.
1686 dellaqc_ens(i,k,1)=0.
1687 pwo_ens (i,k,1)=0.
1688 endif
1689 enddo
1690 enddo
1691 250 continue
1692 !
1693 !--- FEEDBACK
1694 !
1695 if(imid.eq.1 .and. ichoice .le.2)then
1696 do i=its,itf
1697 !-boundary layer QE
1698 xff_mid(i,1)=0.
1699 xff_mid(i,2)=0.
1700 if(ierr(i).eq.0)then
1701 blqe=0.
1702 trash=0.
1703 if(k22(i).lt.kpbl(i)+1)then
1704 do k=1,kpbl(i)
1705 blqe=blqe+100.*dhdt(i,k)*(po_cup(i,k)-po_cup(i,k+1))/g
1706 enddo
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))
1710 endif
1711 xff_mid(i,2)=min(0.1,.03*zws(i))
1712 endif
1713 enddo
1714 endif
1715 call cup_output_ens_3d(xff_mid,xf_ens,ierr,dellat_ens,dellaq_ens, &
1716 dellaqc_ens,outt, &
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, &
1722 its,ite, kts,kte, &
1723 dicycle,xf_dicycle )
1724 k=1
1725 do i=its,itf
1726 if(ierr(i).eq.0 .and.pre(i).gt.0.) then
1727 PRE(I)=MAX(PRE(I),0.)
1728 xmb_out(i)=xmb(i)
1729 do k=kts,ktop(i)
1730 outu(i,k)=dellu(i,k)*xmb(i)
1731 outv(i,k)=dellv(i,k)*xmb(i)
1732 enddo
1733 elseif(ierr(i).ne.0 .or. pre(i).eq.0.)then
1734 ktop(i)=0
1735 do k=kts,kte
1736 outt(i,k)=0.
1737 outq(i,k)=0.
1738 outqc(i,k)=0.
1739 outu(i,k)=0.
1740 outv(i,k)=0.
1741 enddo
1742 endif
1743 enddo
1744 ! rain evaporation as in SAS
1745 !
1746 if(irainevap.eq.1)then
1747 do i = its,itf
1748 rntot(i) = 0.
1749 delqev(i) = 0.
1750 delq2(i) = 0.
1751 rn(i) = 0.
1752 rntot(i) = 0.
1753 rain=0.
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
1758 enddo
1759 endif
1760 enddo
1761 do i = its,itf
1762 qevap(i) = 0.
1763 flg(i) = .true.
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
1770 if(flg(i))then
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
1780 endif
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
1784 flg(i) = .false.
1785 endif
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
1793 endif
1794 endif
1795 enddo
1796 ! pre(i)=1000.*rn(i)/dtime
1797 endif
1798 enddo
1799 endif
1800 !
1801 ! since kinetic energy is being dissipated, add heating accordingly (from ECMWF)
1802 !
1803 do i=its,itf
1804 if(ierr(i).eq.0) then
1805 dts=0.
1806 fpi=0.
1807 do k=kts,ktop(i)
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
1813 enddo
1814 if(fpi.gt.0.)then
1815 do k=kts,ktop(i)
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
1818 enddo
1819 endif
1820 endif
1821 enddo
1823 #if ( WRF_CHEM == 1 )
1824 !--- calculate in-cloud/updraft air temperature
1825 do i=its,itf
1826 if (ierr(i)==0) then
1827 do k=kts,ktf
1828 tempco(i,k)=(1./cp)*(hco(i,k)-g*zo_cup(i,k)-xlv*qco(i,k))
1829 enddo
1830 else
1831 do k=kts,ktf
1832 tempco(i,k)=tn_cup(i,k)
1833 enddo
1834 endif
1835 enddo
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 &
1839 ,po,po_cup,zo_cup &
1840 ,zuo,zdo,pwo,pwdo,pwevo,pwavo &
1841 ,up_massentro,up_massdetro &
1842 ,dd_massentro,dd_massdetro &
1843 ,tempco,clw_all &
1844 ,ktop,k22,kbcon,jmin &
1845 ,xmb,ierr,edto &
1846 ,itf,ktf,its,ite,kts,kte &
1847 ,0)
1848 endif
1849 if ((chem_conv_tr>0).and.(traceropt>0)) then
1850 call ctrans_gf(0,num_tracer,tracer2d,0,traceropt &
1851 ,outtracert,0,0 &
1852 ,po,po_cup,zo_cup &
1853 ,zuo,zdo,pwo,pwdo,pwevo,pwavo &
1854 ,up_massentro,up_massdetro &
1855 ,dd_massentro,dd_massdetro &
1856 ,tempco,clw_all &
1857 ,ktop,k22,kbcon,jmin &
1858 ,xmb,ierr,edto &
1859 ,itf,ktf,its,ite,kts,kte &
1860 ,0)
1861 endif
1862 #endif
1864 !
1865 !---------------------------done------------------------------
1866 !
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, &
1874 its,ite, kts,kte )
1876 IMPLICIT NONE
1878 integer &
1879 ,intent (in ) :: &
1880 aeroevap,itf,ktf, &
1881 its,ite, kts,kte
1882 !
1883 ! ierr error value, maybe modified in this routine
1884 !
1885 real, dimension (its:ite,kts:kte) &
1886 ,intent (in ) :: &
1887 rho,us,vs,z,p,pw
1888 real, dimension (its:ite,1) &
1889 ,intent (out ) :: &
1890 edtc
1891 real, dimension (its:ite) &
1892 ,intent (out ) :: &
1893 edt
1894 real, dimension (its:ite) &
1895 ,intent (in ) :: &
1896 pwav,pwev,ccn,psum2,psumh,edtmax,edtmin
1897 integer, dimension (its:ite) &
1898 ,intent (in ) :: &
1899 ktop,kbcon
1900 integer, dimension (its:ite) &
1901 ,intent (inout) :: &
1902 ierr
1903 !
1904 ! local variables in this routine
1905 !
1907 integer i,k,kk
1908 real einc,pef,pefb,prezk,zkbc
1909 real, dimension (its:ite) :: &
1910 vshear,sdp,vws
1911 real :: prop_c,pefc,aeroadd,alpha3,beta3
1912 prop_c=8. !10.386
1913 alpha3 = 1.9
1914 beta3 = -1.13
1915 pefc=0.
1917 !
1918 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
1919 !
1920 ! */ calculate an average wind shear over the depth of the cloud
1921 !
1922 do i=its,itf
1923 edt(i)=0.
1924 vws(i)=0.
1925 sdp(i)=0.
1926 vshear(i)=0.
1927 enddo
1928 do i=its,itf
1929 edtc(i,1)=0.
1930 enddo
1931 do kk = kts,ktf-1
1932 do 62 i=its,itf
1933 IF(ierr(i).ne.0)GO TO 62
1934 if (kk .le. min0(ktop(i),ktf) .and. kk .ge. kbcon(i)) then
1935 vws(i) = vws(i)+ &
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)
1940 endif
1941 if (kk .eq. ktf-1)vshear(i) = 1.e3 * vws(i) / sdp(i)
1942 62 continue
1943 end do
1944 do i=its,itf
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
1950 !
1951 !--- cloud base precip efficiency
1952 !
1953 zkbc=z(i,kbcon(i))*3.281e-3
1954 prezk=.02
1955 if(zkbc.gt.3.)then
1956 prezk=.96729352+zkbc*(-.70034167+zkbc*(.162179896+zkbc &
1957 *(- 1.2569798E-2+zkbc*(4.2772E-4-zkbc*5.44E-6))))
1958 endif
1959 if(zkbc.gt.25)then
1960 prezk=2.4
1961 endif
1962 pefb=1./(1.+prezk)
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
1968 ! prop_c=.9/aeroadd
1969 prop_c=.5*(pefb+pef)/aeroadd
1970 aeroadd=(ccn(i)**beta3)*((psum2(i))**(alpha3-1)) !*1.e6
1971 aeroadd=prop_c*aeroadd
1972 pefc=aeroadd
1973 if(pefc.gt.0.9)pefc=0.9
1974 if(pefc.lt.0.1)pefc=0.1
1975 EDT(I)=1.-pefc
1976 if(aeroevap.eq.2)EDT(I)=1.-.25*(pefb+pef+2.*pefc)
1977 endif
1980 !--- edt here is 1-precipeff!
1981 einc=.2*edt(i)
1982 edtc(i,1)=edt(i)-einc
1983 endif
1984 enddo
1985 do i=its,itf
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)
1990 endif
1991 enddo
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, &
1999 q,he,iloop, &
2000 itf,ktf, &
2001 its,ite, kts,kte )
2003 IMPLICIT NONE
2005 integer &
2006 ,intent (in ) :: &
2007 itf,ktf, &
2008 its,ite, kts,kte
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
2015 ! bu = buoancy term
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
2024 !
2025 real, dimension (its:ite,kts:kte) &
2026 ,intent (in ) :: &
2027 zd,hes_cup,hcd,qes_cup,q_cup,z_cup, &
2028 dd_massentr,dd_massdetr,gamma_cup,q,he
2029 integer &
2030 ,intent (in ) :: &
2031 iloop
2032 integer, dimension (its:ite) &
2033 ,intent (in ) :: &
2034 jmin
2035 integer, dimension (its:ite) &
2036 ,intent (inout) :: &
2037 ierr
2038 real, dimension (its:ite,kts:kte)&
2039 ,intent (out ) :: &
2040 qcd,qrcd,pwd
2041 real, dimension (its:ite)&
2042 ,intent (out ) :: &
2043 pwev,bu
2044 character*50 :: ierrc(its:ite)
2045 !
2046 ! local variables in this routine
2047 !
2049 integer :: &
2050 i,k,ki
2051 real :: &
2052 denom,dh,dz,dqeva
2054 do i=its,itf
2055 bu(i)=0.
2056 pwev(i)=0.
2057 enddo
2058 do k=kts,ktf
2059 do i=its,itf
2060 qcd(i,k)=0.
2061 qrcd(i,k)=0.
2062 pwd(i,k)=0.
2063 enddo
2064 enddo
2065 !
2066 !
2067 !
2068 do 100 i=its,itf
2069 IF(ierr(I).eq.0)then
2070 k=jmin(i)
2071 DZ=Z_cup(i,K+1)-Z_cup(i,K)
2072 qcd(i,k)=q_cup(i,k)
2073 DH=HCD(I,k)-HES_cup(I,K)
2074 if(dh.lt.0)then
2075 QRCD(I,K)=(qes_cup(i,k)+(1./XLV)*(GAMMA_cup(i,k) &
2076 /(1.+GAMMA_cup(i,k)))*DH)
2077 else
2078 qrcd(i,k)=qes_cup(i,k)
2079 endif
2080 pwd(i,jmin(i))=zd(i,jmin(i))*min(0.,qcd(i,k)-qrcd(i,k))
2081 qcd(i,k)=qrcd(i,k)
2082 pwev(i)=pwev(i)+pwd(i,jmin(i)) ! *dz
2083 !
2084 bu(i)=dz*dh
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)
2090 ! dz=qcd(i,ki)
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
2096 ierr(i)=51
2097 exit
2098 endif
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))
2103 !
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
2106 !--- the case
2107 !
2108 DH=HCD(I,ki)-HES_cup(I,Ki)
2109 bu(i)=bu(i)+dz*dh
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)
2113 if(dqeva.gt.0.)then
2114 dqeva=0.
2115 qrcd(i,ki)=qcd(i,ki)
2116 endif
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
2122 ! endif
2123 enddo
2124 !
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
2128 ierr(i)=7
2129 ierrc(i)="problem with buoy in cup_dd_moisture"
2130 endif
2131 if(BU(I).GE.0.and.iloop.eq.1)then
2132 ! print *,'problem with buoy in cup_dd_moisture',i
2133 ierr(i)=7
2134 ierrc(i)="problem2 with buoy in cup_dd_moisture"
2135 endif
2136 endif
2137 100 continue
2139 END SUBROUTINE cup_dd_moisture
2141 SUBROUTINE cup_env(z,qes,he,hes,t,q,p,z1, &
2142 psur,ierr,tcrit,itest, &
2143 itf,ktf, &
2144 its,ite, kts,kte )
2146 IMPLICIT NONE
2148 integer &
2149 ,intent (in ) :: &
2150 itf,ktf, &
2151 its,ite, kts,kte
2152 !
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
2164 !
2165 !
2166 real, dimension (its:ite,kts:kte) &
2167 ,intent (in ) :: &
2168 p,t,q
2169 real, dimension (its:ite,kts:kte) &
2170 ,intent (out ) :: &
2171 he,hes,qes
2172 real, dimension (its:ite,kts:kte) &
2173 ,intent (inout) :: &
2174 z
2175 real, dimension (its:ite) &
2176 ,intent (in ) :: &
2177 psur,z1
2178 integer, dimension (its:ite) &
2179 ,intent (inout) :: &
2180 ierr
2181 integer &
2182 ,intent (in ) :: &
2183 itest
2184 !
2185 ! local variables in this routine
2186 !
2188 integer :: &
2189 i,k
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
2194 ! real :: satvap
2197 ! HT(1)=XLV/CP
2198 ! HT(2)=2.834E6/CP
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)
2203 do k=kts,ktf
2204 do i=its,itf
2205 if(ierr(i).eq.0)then
2206 !Csgb - IPH is for phase, dependent on TCRIT (water or ice)
2207 ! IPH=1
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)
2213 e=satvap(t(i,k))
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)
2219 endif
2220 enddo
2221 enddo
2222 !
2223 !--- z's are calculated with changed h's and q's and t's
2224 !--- if itest=2
2225 !
2226 if(itest.eq.1 .or. itest.eq.0)then
2227 do i=its,itf
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
2231 endif
2232 enddo
2234 ! --- calculate heights
2235 DO K=kts+1,ktf
2236 do i=its,itf
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
2241 endif
2242 enddo
2243 enddo
2244 else if(itest.eq.2)then
2245 do k=kts,ktf
2246 do i=its,itf
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))
2250 endif
2251 enddo
2252 enddo
2253 else if(itest.eq.-1)then
2254 endif
2255 !
2256 !--- calculate moist static energy - HE
2257 ! saturated moist static energy - HES
2258 !
2259 DO k=kts,ktf
2260 do i=its,itf
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)
2265 endif
2266 enddo
2267 enddo
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, &
2274 ierr,z1, &
2275 itf,ktf, &
2276 its,ite, kts,kte )
2278 IMPLICIT NONE
2280 integer &
2281 ,intent (in ) :: &
2282 itf,ktf, &
2283 its,ite, kts,kte
2284 !
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
2303 !
2304 !
2305 real, dimension (its:ite,kts:kte) &
2306 ,intent (in ) :: &
2307 qes,q,he,hes,z,p,t
2308 real, dimension (its:ite,kts:kte) &
2309 ,intent (out ) :: &
2310 qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup
2311 real, dimension (its:ite) &
2312 ,intent (in ) :: &
2313 psur,z1
2314 integer, dimension (its:ite) &
2315 ,intent (inout) :: &
2316 ierr
2317 !
2318 ! local variables in this routine
2319 !
2321 integer :: &
2322 i,k
2325 do k=kts,ktf
2326 do i=its,itf
2327 qes_cup(i,k)=0.
2328 q_cup(i,k)=0.
2329 hes_cup(i,k)=0.
2330 he_cup(i,k)=0.
2331 z_cup(i,k)=0.
2332 p_cup(i,k)=0.
2333 t_cup(i,k)=0.
2334 gamma_cup(i,k)=0.
2335 enddo
2336 enddo
2337 do k=kts+1,ktf
2338 do i=its,itf
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)
2350 endif
2351 enddo
2352 enddo
2353 do i=its,itf
2354 if(ierr(i).eq.0)then
2355 qes_cup(i,1)=qes(i,1)
2356 q_cup(i,1)=q(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))
2363 z_cup(i,1)=z1(i)
2364 p_cup(i,1)=psur(i)
2365 t_cup(i,1)=t(i,1)
2366 gamma_cup(i,1)=xlv/cp*(xlv/(r_v*t_cup(i,1) &
2367 *t_cup(i,1)))*qes_cup(i,1)
2368 endif
2369 enddo
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, &
2376 ichoice, &
2377 imid,ipr,itf,ktf, &
2378 its,ite, kts,kte, &
2379 dicycle,tau_ecmwf,aa1_bl,xf_dicycle )
2381 IMPLICIT NONE
2383 integer &
2384 ,intent (in ) :: &
2385 imid,ipr,itf,ktf, &
2386 its,ite, kts,kte
2387 integer, intent (in ) :: &
2388 maxens3
2389 !
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)
2401 ! edt = epsilon
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!)
2409 !
2410 real, dimension (its:ite,1:maxens3) &
2411 ,intent (inout) :: &
2412 pr_ens
2413 real, dimension (its:ite,1:maxens3) &
2414 ,intent (inout ) :: &
2415 xf_ens
2416 real, dimension (its:ite,kts:kte) &
2417 ,intent (in ) :: &
2418 zd,zu,p_cup
2419 real, dimension (its:ite,kts:kte) &
2420 ,intent (in ) :: &
2421 omeg
2422 real, dimension (its:ite,1) &
2423 ,intent (in ) :: &
2424 xaa0
2425 real, dimension (its:ite,4) &
2426 ,intent (in ) :: &
2427 rand_clos
2428 real, dimension (its:ite) &
2429 ,intent (in ) :: &
2430 aa1,edt
2431 real, dimension (its:ite) &
2432 ,intent (in ) :: &
2433 mconv,axx
2434 real, dimension (its:ite) &
2435 ,intent (inout) :: &
2436 aa0,closure_n
2437 real &
2438 ,intent (in ) :: &
2439 mbdt
2440 real &
2441 ,intent (in ) :: &
2442 dtime
2443 integer, dimension (its:ite) &
2444 ,intent (inout ) :: &
2445 k22,kbcon,ktop
2446 integer, dimension (its:ite) &
2447 ,intent (in ) :: &
2448 xland
2449 integer, dimension (its:ite) &
2450 ,intent (inout) :: &
2451 ierr,ierr2,ierr3
2452 integer &
2453 ,intent (in ) :: &
2454 ichoice
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
2459 !- local var
2460 real :: xff_dicycle
2461 !
2462 ! local variables in this routine
2463 !
2465 real, dimension (1:maxens3) :: &
2466 xff_ens3
2467 real, dimension (1) :: &
2468 xk
2469 integer :: &
2470 kk,i,k,n,ne
2471 ! integer, parameter :: mkxcrt=15
2472 ! real, dimension(1:mkxcrt) :: &
2473 ! pcrit,acrit,acritt
2474 integer, dimension (its:ite) :: kloc
2475 real :: &
2476 a1,a_ave,xff0,xomg!,aclim1,aclim2,aclim3,aclim4
2478 real, dimension (its:ite) :: ens_adj
2482 !
2483 ens_adj(:)=1.
2484 xff_dicycle = 0.
2486 !--- LARGE SCALE FORCING
2487 !
2488 DO 100 i=its,itf
2489 kloc(i)=1
2490 IF(ierr(i).eq.0)then
2491 kloc(i)=maxloc(zu(i,:),1)
2492 ens_adj(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
2496 !
2497 a_ave=0.
2498 a_ave=axx(i)
2499 a_ave=max(0.,a_ave)
2500 a_ave=min(a_ave,aa1(i))
2501 a_ave=max(0.,a_ave)
2502 xff_ens3(:)=0.
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)
2509 !
2510 !--- omeg is in bar/s, mconv done with omeg in Pa/s
2511 ! more like Brown (1979), or Frank-Cohen (199?)
2512 !
2513 ! average aaround kbcon
2514 !
2515 xomg=0.
2516 kk=0
2517 xff_ens3(4)=0.
2518 xff_ens3(5)=0.
2519 xff_ens3(6)=0.
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)))
2523 kk=kk+1
2524 endif
2525 enddo
2526 if(kk.gt.0)xff_ens3(4)=xomg/float(kk)
2528 !
2529 ! max below kbcon
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
2534 ! enddo
2535 !
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)
2545 !
2546 !--- more like Krishnamurti et al.; pick max and average values
2547 !
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)
2553 !
2554 !--- more like Fritsch Chappel or Kain Fritsch (plus triggers)
2555 !
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)
2563 !gtest
2564 if(ichoice.eq.0)then
2565 if(xff0.lt.0.)then
2566 xff_ens3(1)=0.
2567 xff_ens3(2)=0.
2568 xff_ens3(3)=0.
2569 xff_ens3(10)=0.
2570 xff_ens3(11)=0.
2571 xff_ens3(12)=0.
2572 xff_ens3(13)= 0.
2573 xff_ens3(16)= 0.
2574 closure_n(i)=12.
2575 ! xff_dicycle = 0.
2576 endif !xff0
2577 endif ! ichoice
2579 XK(1)=(XAA0(I,1)-AA1(I))/MBDT
2580 forcing(i,4)=aa0(i)
2581 forcing(i,5)=aa1(i)
2582 forcing(i,6)=xaa0(i,1)
2583 forcing(i,7)=xk(1)
2584 if(xk(1).le.0.and.xk(1).gt.-.01*mbdt) &
2585 xk(1)=-.01*mbdt
2586 if(xk(1).gt.0.and.xk(1).lt.1.e-2) &
2587 xk(1)=1.e-2
2588 ! enddo
2589 !
2590 !--- add up all ensembles
2591 !
2592 !
2593 ! over water, enfor!e small cap for some of the closures
2594 !
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)
2613 !srf
2614 xff_dicycle = ens_adj(i)*xff_dicycle
2615 !srf end
2616 ! xff_ens3(7) =0.
2617 ! xff_ens3(8) =0.
2618 ! xff_ens3(9) =0.
2619 endif ! ierr2
2620 endif ! xland
2621 !
2622 ! end water treatment
2623 !
2624 !
2626 !
2627 !--- special treatment for stability closures
2628 !
2629 if(XK(1).lt.0.)then
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)
2638 else
2639 xff_ens3(1)= 0
2640 xff_ens3(2)= 0
2641 xff_ens3(3)= 0
2642 xff_ens3(16)=0
2643 endif
2644 !
2645 !--- if iresult.eq.1, following independent of xff0
2646 !
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)
2668 if(XK(1).lt.0.)then
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)
2678 else
2679 xf_ens(i,10)=0.
2680 xf_ens(i,11)=0.
2681 xf_ens(i,12)=0.
2682 xf_ens(i,13)=0.
2683 forcing(i,8)=0.
2684 endif
2685 !srf-begin
2686 if(XK(1).lt.0.)then
2687 xf_dicycle(i) = max(0.,-xff_dicycle /xk(1))
2688 ! forcing(i,9)=xf_dicycle(i)
2689 else
2690 xf_dicycle(i) = 0.
2691 endif
2692 !srf-end
2693 if(ichoice.ge.1)then
2694 ! closure_n(i)=0.
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)
2711 endif
2712 elseif(ierr(i).ne.20.and.ierr(i).ne.0)then
2713 do n=1,maxens3
2714 xf_ens(i,n)=0.
2715 xf_dicycle(i) = 0.
2716 enddo
2717 endif ! ierror
2718 100 continue
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, &
2724 ztexec,zqexec, &
2725 jprnt,itf,ktf, &
2726 its,ite, kts,kte, &
2727 z_cup,entr_rate,heo,imid )
2729 IMPLICIT NONE
2730 !
2732 ! only local wrf dimensions are need as of now in this routine
2734 integer &
2735 ,intent (in ) :: &
2736 jprnt,itf,ktf,imid, &
2737 its,ite, kts,kte
2738 !
2739 !
2740 !
2741 ! ierr error value, maybe modified in this routine
2742 !
2743 real, dimension (its:ite,kts:kte) &
2744 ,intent (in ) :: &
2745 he_cup,hes_cup,p_cup
2746 real, dimension (its:ite) &
2747 ,intent (in ) :: &
2748 entr_rate,ztexec,zqexec,cap_inc,cap_max
2749 real, dimension (its:ite) &
2750 ,intent (inout ) :: &
2751 hkb !,cap_max
2752 integer, dimension (its:ite) &
2753 ,intent (in ) :: &
2754 kbmax
2755 integer, dimension (its:ite) &
2756 ,intent (inout) :: &
2757 kbcon,k22,ierr
2758 integer &
2759 ,intent (in ) :: &
2760 iloop_in
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
2764 !
2765 ! local variables in this routine
2766 !
2768 integer :: &
2769 i,k
2770 real :: &
2771 x_add,pbcdif,plus,hetest,dz
2772 real, dimension (its:ite,kts:kte) ::hcot
2773 !
2774 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
2775 !
2776 iloop(:)=iloop_in
2777 DO 27 i=its,itf
2778 kbcon(i)=1
2779 !
2780 ! reset iloop for mid level convection
2781 if(cap_max(i).gt.200 .and. imid.eq.1)iloop(i)=5
2782 !
2783 IF(ierr(I).ne.0)GO TO 27
2784 start_level(i)=k22(i)
2785 KBCON(I)=K22(I)+1
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)
2795 enddo
2796 !==
2798 GO TO 32
2799 31 CONTINUE
2800 KBCON(I)=KBCON(I)+1
2801 IF(KBCON(I).GT.KBMAX(i)+2)THEN
2802 if(iloop(i).ne.4)then
2803 ierr(i)=3
2804 ierrc(i)="could not find reasonable kbcon in cup_kbcon"
2805 endif
2806 GO TO 27
2807 ENDIF
2808 32 CONTINUE
2809 hetest=hcot(i,kbcon(i)) !hkb(i) ! HE_cup(I,K22(I))
2810 IF(HETEST.LT.HES_cup(I,KBCON(I)))then
2811 GO TO 31
2812 endif
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)
2821 !
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
2826 Go To 27
2827 ElseIF(PBCDIF.GT.plus)THEN
2828 K22(I)=K22(I)+1
2829 KBCON(I)=K22(I)+1
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)
2843 enddo
2844 !==
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
2849 ierr(i)=3
2850 ierrc(i)="could not find reasonable kbcon in cup_kbcon"
2851 endif
2852 GO TO 27
2853 ENDIF
2854 GO TO 32
2855 ENDIF
2856 27 CONTINUE
2858 END SUBROUTINE cup_kbcon
2861 SUBROUTINE cup_MAXIMI(ARRAY,KS,KE,MAXX,ierr, &
2862 itf,ktf, &
2863 its,ite, kts,kte )
2865 IMPLICIT NONE
2866 !
2867 ! on input
2868 !
2870 ! only local wrf dimensions are need as of now in this routine
2872 integer &
2873 ,intent (in ) :: &
2874 itf,ktf, &
2875 its,ite, kts,kte
2876 ! array input array
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) &
2881 ,intent (in ) :: &
2882 array
2883 integer, dimension (its:ite) &
2884 ,intent (in ) :: &
2885 ierr,ke
2886 integer &
2887 ,intent (in ) :: &
2888 ks
2889 integer, dimension (its:ite) &
2890 ,intent (out ) :: &
2891 maxx
2892 real, dimension (its:ite) :: &
2893 x
2894 real :: &
2895 xar
2896 integer :: &
2897 i,k
2899 DO 200 i=its,itf
2900 MAXX(I)=KS
2901 if(ierr(i).eq.0)then
2902 X(I)=ARRAY(I,KS)
2903 !
2904 DO 100 K=KS,KE(i)
2905 XAR=ARRAY(I,K)
2906 IF(XAR.GE.X(I)) THEN
2907 X(I)=XAR
2908 MAXX(I)=K
2909 ENDIF
2910 100 CONTINUE
2911 endif
2912 200 CONTINUE
2914 END SUBROUTINE cup_MAXIMI
2917 SUBROUTINE cup_minimi(ARRAY,KS,KEND,KT,ierr, &
2918 itf,ktf, &
2919 its,ite, kts,kte )
2921 IMPLICIT NONE
2922 !
2923 ! on input
2924 !
2926 ! only local wrf dimensions are need as of now in this routine
2928 integer &
2929 ,intent (in ) :: &
2930 itf,ktf, &
2931 its,ite, kts,kte
2932 ! array input array
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) &
2937 ,intent (in ) :: &
2938 array
2939 integer, dimension (its:ite) &
2940 ,intent (in ) :: &
2941 ierr,ks,kend
2942 integer, dimension (its:ite) &
2943 ,intent (out ) :: &
2944 kt
2945 real, dimension (its:ite) :: &
2946 x
2947 integer :: &
2948 i,k,kstop
2950 DO 200 i=its,itf
2951 KT(I)=KS(I)
2952 if(ierr(i).eq.0)then
2953 X(I)=ARRAY(I,KS(I))
2954 KSTOP=MAX(KS(I)+1,KEND(I))
2955 !
2956 DO 100 K=KS(I)+1,KSTOP
2957 IF(ARRAY(I,K).LT.X(I)) THEN
2958 X(I)=ARRAY(I,K)
2959 KT(I)=K
2960 ENDIF
2961 100 CONTINUE
2962 endif
2963 200 CONTINUE
2965 END SUBROUTINE cup_MINIMI
2968 SUBROUTINE cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
2969 kbcon,ktop,ierr, &
2970 itf,ktf, &
2971 its,ite, kts,kte )
2973 IMPLICIT NONE
2974 !
2975 ! on input
2976 !
2978 ! only local wrf dimensions are need as of now in this routine
2980 integer &
2981 ,intent (in ) :: &
2982 itf,ktf, &
2983 its,ite, kts,kte
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
2991 !
2992 real, dimension (its:ite,kts:kte) &
2993 ,intent (in ) :: &
2994 z,zu,gamma_cup,t_cup,dby
2995 integer, dimension (its:ite) &
2996 ,intent (in ) :: &
2997 kbcon,ktop
2998 !
2999 ! input and output
3000 !
3003 integer, dimension (its:ite) &
3004 ,intent (inout) :: &
3005 ierr
3006 real, dimension (its:ite) &
3007 ,intent (out ) :: &
3008 aa0
3009 !
3010 ! local variables in this routine
3011 !
3013 integer :: &
3014 i,k
3015 real :: &
3016 dz,da
3017 !
3018 do i=its,itf
3019 aa0(i)=0.
3020 enddo
3021 DO 100 k=kts+1,ktf
3022 DO 100 i=its,itf
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
3026 DZ=Z(I,K)-Z(I,K-1)
3027 da=zu(i,k)*DZ*(9.81/(1004.*( &
3028 (T_cup(I,K)))))*DBY(I,K-1)/ &
3029 (1.+GAMMA_CUP(I,K))
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.
3033 100 continue
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 ) , &
3044 intent(inout ) :: &
3045 outq,outt,outqc,outu,outv
3046 real, dimension (its:ite,kts:kte ) , &
3047 intent(inout ) :: &
3048 q
3049 real, dimension (its:ite ) , &
3050 intent(inout ) :: &
3051 pret
3052 character *(*), intent (in) :: &
3053 name
3054 real &
3055 ,intent (in ) :: &
3056 dt
3057 real :: names,scalef,thresh,qmem,qmemf,qmem2,qtest,qmem1
3058 integer :: icheck
3059 !
3060 ! first do check on vertical heating rate
3061 !
3062 thresh=300.01
3063 ! thresh=200.01 !ss
3064 ! thresh=250.01
3065 names=1.
3066 if(name == 'shallow')then
3067 thresh=148.01
3068 names=2.
3069 endif
3070 scalef=86400.
3071 do i=its,itf
3072 icheck=0
3073 qmemf=1.
3074 qmem=0.
3075 do k=kts,ktf
3076 qmem=(outt(i,k))*86400.
3077 if(qmem.gt.thresh)then
3078 qmem2=thresh/qmem
3079 qmemf=min(qmemf,qmem2)
3080 icheck=1
3081 !
3082 !
3083 ! print *,'1',' adjusted massflux by factor ',i,j,k,qmem,qmem2,qmemf,dt
3084 endif
3085 if(qmem.lt.-.5*thresh*names)then
3086 qmem2=-.5*names*thresh/qmem
3087 qmemf=min(qmemf,qmem2)
3088 icheck=2
3089 !
3090 !
3091 endif
3092 enddo
3093 do k=kts,ktf
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
3099 enddo
3100 pret(i)=pret(i)*qmemf
3101 enddo
3102 return
3103 !
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
3107 ! tendencies
3108 !
3109 ! return
3110 ! write(14,*)'return'
3111 thresh=1.e-16
3112 do i=its,itf
3113 qmemf=1.
3114 do k=kts,ktf-1
3115 qmem=outq(i,k)
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
3119 !
3120 ! qmem2 would be the maximum allowable tendency
3121 !
3122 qmem1=abs(outq(i,k))
3123 qmem2=abs((thresh-q(i,k))/dt)
3124 qmemf=min(qmemf,qmem2/qmem1)
3125 qmemf=max(0.,qmemf)
3126 endif
3127 endif
3128 enddo
3129 do k=kts,ktf
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
3135 enddo
3136 pret(i)=pret(i)*qmemf
3137 enddo
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, &
3146 maxens3, &
3147 sig,closure_n,xland1,xmbm_in,xmbs_in, &
3148 ichoice,imid,ipr,itf,ktf, &
3149 its,ite, kts,kte, &
3150 dicycle,xf_dicycle )
3152 IMPLICIT NONE
3153 !
3154 ! on input
3155 !
3156 ! only local wrf dimensions are need as of now in this routine
3158 integer &
3159 ,intent (in ) :: &
3160 ichoice,imid,ipr,itf,ktf, &
3161 its,ite, kts,kte
3162 integer, intent (in ) :: &
3163 maxens3
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
3177 !
3178 real, dimension (its:ite,1:maxens3) &
3179 ,intent (inout) :: &
3180 xf_ens,pr_ens
3181 real, dimension (its:ite,kts:kte) &
3182 ,intent (inout ) :: &
3183 outtem,outq,outqc
3184 real, dimension (its:ite,kts:kte) &
3185 ,intent (in ) :: &
3186 zu,pwd,p_cup
3187 real, dimension (its:ite) &
3188 ,intent (in ) :: &
3189 sig,xmbm_in,xmbs_in,edt
3190 real, dimension (its:ite,2) &
3191 ,intent (in ) :: &
3192 xff_mid
3193 real, dimension (its:ite) &
3194 ,intent (inout ) :: &
3195 pre,xmb
3196 real, dimension (its:ite) &
3197 ,intent (inout ) :: &
3198 closure_n
3199 real, dimension (its:ite,kts:kte,1) &
3200 ,intent (in ) :: &
3201 dellat,dellaqc,dellaq,pw
3202 integer, dimension (its:ite) &
3203 ,intent (in ) :: &
3204 ktop,xland1
3205 integer, dimension (its:ite) &
3206 ,intent (inout) :: &
3207 ierr,ierr2,ierr3
3208 integer, intent(IN) :: DICYCLE
3209 real, intent(IN), dimension (its:ite) :: xf_dicycle
3210 !
3211 ! local variables in this routine
3212 !
3214 integer :: &
3215 i,k,n
3216 real :: &
3217 clos_wei,dtt,dp,dtq,dtqc,dtpw,dtpwd
3218 real, dimension (its:ite) :: &
3219 xmb_ave,pwtot
3220 !
3221 character *(*), intent (in) :: &
3222 name
3224 !
3225 DO k=kts,kte
3226 do i=its,ite
3227 outtem (i,k)=0.
3228 outq (i,k)=0.
3229 outqc (i,k)=0.
3230 enddo
3231 enddo
3232 do i=its,itf
3233 pre(i)=0.
3234 xmb(i)=0.
3235 enddo
3236 do i=its,itf
3237 IF(ierr(i).eq.0)then
3238 do n=1,maxens3
3239 if(pr_ens(i,n).le.0.)then
3240 xf_ens(i,n)=0.
3241 endif
3242 enddo
3243 endif
3244 enddo
3245 !
3246 !--- calculate ensemble average mass fluxes
3247 !
3249 !
3250 !-- now do feedback
3251 !
3252 !!!!! DEEP Convection !!!!!!!!!!
3253 if(imid.eq.0)then
3254 do i=its,itf
3255 if(ierr(i).eq.0)then
3256 k=0
3257 xmb_ave(i)=0.
3258 do n=1,maxens3
3259 k=k+1
3260 xmb_ave(i)=xmb_ave(i)+xf_ens(i,n)
3261 enddo
3262 xmb_ave(i)=xmb_ave(i)/float(k)
3263 !srf begin
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))
3270 endif
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
3279 ierr(i)=19
3280 endif
3281 ! xfac1(i)=xmb(i)
3282 ! xfac2(i)=xmb(i)
3284 endif
3285 ENDDO
3286 !!!!! NOT SO DEEP Convection !!!!!!!!!!
3287 else ! imid == 1
3288 do i=its,itf
3289 xmb_ave(i)=0.
3290 IF(ierr(i).eq.0)then
3291 ! ! first get xmb_ves, depend on ichoicee
3292 !
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
3296 k=0
3297 do n=1,maxens3
3298 k=k+1
3299 xmb_ave(i)=xmb_ave(i)+xf_ens(i,n)
3300 enddo
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))
3313 endif ! dicycle=1,2
3314 endif ! ierr >0
3315 enddo ! i
3316 endif ! imid=1
3318 do i=its,itf
3319 pwtot(i)=0.
3320 IF(ierr(i).eq.0)then
3321 DO k=kts,ktop(i)
3322 pwtot(i)=pwtot(i)+pw(i,k,1)
3323 enddo
3324 DO k=kts,ktop(i)
3325 dp=100.*(p_cup(i,k)-p_cup(i,k+1))/g
3326 dtt =dellat (i,k,1)
3327 dtq =dellaq (i,k,1)
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
3333 if(dtqc < 0.)then
3334 dtpwd=dtpwd-dellaqc(i,k,1)*dp
3335 dtqc=0.
3336 ! if this is positive, can come from clw detrainment
3337 else
3338 dtpwd=0.
3339 dtqc=dtqc/dp
3340 endif
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
3347 enddo
3348 PRE(I)=-PRE(I)+XMB(I)*pwtot(i)
3349 endif
3350 enddo
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, &
3360 itest,itf,ktf, &
3361 its,ite, kts,kte )
3363 IMPLICIT NONE
3364 real, parameter :: BDISPM = 0.366 !Berry--size dispersion (martime)
3365 REAL, PARAMETER :: BDISPC = 0.146 !Berry--size dispersion (continental)
3366 !
3367 ! on input
3368 !
3370 ! only local wrf dimensions are need as of now in this routine
3372 integer &
3373 ,intent (in ) :: &
3374 itest,itf,ktf, &
3375 its,ite, kts,kte
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
3384 !
3385 real, dimension (its:ite,kts:kte) &
3386 ,intent (in ) :: &
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) &
3390 ,intent (in ) :: &
3391 zqexec
3392 ! entr= entrainment rate
3393 integer, dimension (its:ite) &
3394 ,intent (in ) :: &
3395 kbcon,ktop,k22,xland1
3396 !
3397 ! input and output
3398 !
3400 ! ierr error value, maybe modified in this routine
3402 integer, dimension (its:ite) &
3403 ,intent (inout) :: &
3404 ierr
3405 character *(*), intent (in) :: &
3406 name
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) &
3415 ,intent (out ) :: &
3416 qc,qrc,pw,clw_all
3417 real, dimension (its:ite,kts:kte) :: &
3418 qch,qrcb,pwh,clw_allh,c1d,t
3419 real, dimension (its:ite) :: &
3420 pwavh
3421 real, dimension (its:ite) &
3422 ,intent (out ) :: &
3423 pwav,psum,psumh
3424 real, dimension (its:ite) &
3425 ,intent (in ) :: &
3426 ccn
3427 !
3428 ! local variables in this routine
3429 !
3431 integer :: &
3432 iprop,iall,i,k
3433 integer :: start_level(its:ite)
3434 real :: &
3435 prop_ave,qrcb_h,bdsp,dp,rhoc,qrch,qaver, &
3436 c0,dz,berryc0,q1,berryc
3437 real :: &
3438 denom
3439 real, dimension (kts:kte) :: &
3440 prop_b
3441 !
3442 prop_b(kts:kte)=0
3443 iall=0
3444 c0=.004 ! Han et al. (2016); Li et al., (2019), updated from 0.002
3445 bdsp=BDISPM
3446 !
3447 !--- no precip for small clouds
3448 !
3449 ! if(name.eq.'shallow')then
3450 ! c0=0.002
3451 ! endif
3452 do i=its,itf
3453 pwav(i)=0.
3454 pwavh(i)=0.
3455 psum(i)=0.
3456 psumh(i)=0.
3457 enddo
3458 do k=kts,ktf
3459 do i=its,itf
3460 pw(i,k)=0.
3461 pwh(i,k)=0.
3462 qc(i,k)=0.
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)
3465 clw_all(i,k)=0.
3466 clw_allh(i,k)=0.
3467 qrc(i,k)=0.
3468 qrcb(i,k)=0.
3469 enddo
3470 enddo
3471 do i=its,itf
3472 if(ierr(i).eq.0)then
3473 start_level=k22(i)
3474 call get_cloud_bc(kte,qe_cup (i,1:kte),qaver,k22(i))
3475 qaver = qaver
3476 k=start_level(i)
3477 qc (i,k)= qaver
3478 qch (i,k)= qaver
3479 do k=1,start_level(i)-1
3480 qc (i,k)= qe_cup(i,k)
3481 qch (i,k)= qe_cup(i,k)
3482 enddo
3483 !
3484 ! initialize below originating air
3485 !
3486 endif
3487 enddo
3489 DO 100 i=its,itf
3490 c0=.004
3491 IF(ierr(i).eq.0)then
3493 ! below LFC, but maybe above LCL
3494 !
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))
3501 ! QRCH=QES_cup(I,K)
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)
3511 endif
3512 enddo
3513 ! endif
3514 !
3515 !now do the rest
3516 !
3517 DO k=kbcon(i)+1,ktop(i)
3518 c0=.004
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
3522 ierr(i)=51
3523 exit
3524 endif
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)
3530 !
3531 !--- saturation in cloud, this is what is allowed to be in it
3532 !
3533 QRCH=QES_cup(I,K)+(1./XLV)*(GAMMA_cup(i,k) &
3534 /(1.+GAMMA_cup(i,k)))*DBY(I,K)
3535 !
3536 !------ 1. steady state plume equation, for what could
3537 !------ be in cloud without condensation
3538 !
3539 !
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
3548 qc(i,k)=qrch
3549 endif
3550 if(qch(i,k).le.qrch)then
3551 qch(i,k)=qrch
3552 endif
3553 !
3554 !------- Total condensed water before rainout
3555 !
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
3563 !
3564 ! normalized berry
3565 !
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
3569 !
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.
3577 berryc=qrcb(i,k)
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)
3583 qrcb(i,k)=0.
3584 endif
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
3588 !
3589 ! then the real berry
3590 !
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.
3595 berryc=qrc(i,k)
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)
3600 qrc(i,k)=0.
3601 endif
3602 pw(i,k)=berryc0*zu(i,k)
3603 QC(I,K)=QRC(I,K)+qrch
3604 !
3605 ! if not running with berry at all, do the following
3606 !
3607 ELSE !c0=.002
3608 if(iall.eq.1)then
3609 qrc(i,k)=0.
3610 pw(i,k)=(QC(I,K)-QRCH)*zu(i,k)
3611 if(pw(i,k).lt.0.)pw(i,k)=0.
3612 else
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
3616 qrc(i,k)=0.
3617 pw(i,k)=0.
3618 endif
3619 endif
3620 QC(I,K)=QRC(I,K)+qrch
3621 endif !autoconv
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)
3628 enddo
3629 endif ! ierr
3630 !
3631 !--- integrated normalized ondensate
3632 !
3633 100 CONTINUE
3634 prop_ave=0.
3635 iprop=0
3636 do k=kts,kte
3637 prop_ave=prop_ave+prop_b(k)
3638 if(prop_b(k).gt.0)iprop=iprop+1
3639 enddo
3640 iprop=max(iprop,1)
3642 END SUBROUTINE cup_up_moisture
3644 !--------------------------------------------------------------------
3646 REAL FUNCTION satvap(temp2)
3647 implicit none
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
3653 toto = 1 / toot
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
3657 else
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
3667 end if
3668 END FUNCTION
3669 !--------------------------------------------------------------------
3670 SUBROUTINE get_cloud_bc(mzp,array,x_aver,k22,add)
3671 implicit none
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)
3686 x_aver=0.
3687 do i = 1,local_order_aver
3688 x_aver = x_aver + array(k22-i+1)
3689 enddo
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)
3699 implicit none
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
3707 !-local vars
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
3714 !
3715 zustart=.1
3716 dbythresh= 1. !.0.95 ! 0.85, 0.6
3717 if(name == 'shallow' .or. name == 'mid') dbythresh=1.
3718 dby(:)=0.
3720 DO i=its,itf
3721 zux(:)=0.
3722 beta_u=max(.1,.2-float(csum(i))*.01)
3723 zuo(i,:)=0.
3724 dby(:)=0.
3725 dbm(:)=0.
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
3736 zux(k)=zuo(i,k)
3737 enddo
3738 zubeg=zustart !zuo(i,kbcon(i))
3739 if(name .eq. 'deep')then
3740 ktop(i)=0
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)
3751 enddo
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
3756 kfinalzu=k - 1
3757 ktop(i)=kfinalzu
3758 go to 412
3759 endif
3760 enddo
3761 kfinalzu=ktf-2
3762 ktop(i)=kfinalzu
3763 412 continue
3764 !
3765 ! at least overshoot by one level
3766 !
3767 ! kfinalzu=min(max(kfinalzu,ktopdby(i)+1),ktopdby(i)+2)
3768 ! ktop(i)=kfinalzu
3769 if(kfinalzu.le.kbcon(i)+2)then
3770 ierr(i)=41
3771 ktop(i)= 0
3772 else
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))
3778 endif
3779 endif ! end deep
3780 if ( name == 'mid' ) then
3781 if(ktop(i) <= kbcon(i)+2)then
3782 ierr(i)=41
3783 ktop(i)= 0
3784 else
3785 kfinalzu=ktop(i)
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))
3788 ! kbegin=0
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
3792 ! do k=2,40
3793 ! zh2(k)=zh2(k-1)+dzh2
3794 ! if(zuh2(k).gt.0.1 .and. kbegin.eq.0)kbegin=k
3795 ! enddo
3796 ! zuo(i,k22(i))=zuh2(kbegin)
3797 ! do k=k22(i)+1,kfinalzu+1
3798 ! do kk=kbegin,39
3799 ! if(z_cup(i,k).gt.zh2(kk) .and. z_cup(i,k).le.zh2(kk+1)) then
3800 ! zuo(i,k)=zuh2(kk)
3801 ! exit
3802 ! endif
3803 ! enddo
3804 ! enddo
3805 ! if(zuo(i,ktop(i)).lt.1.e-4)ktop(i)=ktop(i)-1
3806 endif
3807 endif ! mid
3808 if ( name == 'shallow' ) then
3809 if(ktop(i) <= kbcon(i)+2)then
3810 ierr(i)=41
3811 ktop(i)= 0
3812 else
3813 kfinalzu=ktop(i)
3814 ktopdby(i)=ktop(i)
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))
3818 endif
3819 endif ! shal
3820 ENDIF ! ierr
3821 ENDDO
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)
3827 implicit none
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
3836 !- local var
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
3842 zu=0.0
3843 zuh=0.0
3844 kb_adj=max(kb,2)
3845 IF(draft == "UP") then
3846 lev_start=min(.9,.4+csum*.013)
3847 kb_adj=max(kb,2)
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))
3855 #else
3856 call wrf_error_fatal ('compiler does not support 2008 gamma intrinsic')
3857 #endif
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)
3861 enddo
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
3867 kb_adj=k+1
3868 exit
3869 endif
3870 enddo
3871 kb_adj=max(2,kb_adj)
3872 do k=kts,kb_adj-1
3873 zu(k)=0.
3874 enddo
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))
3883 #endif
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)
3887 enddo
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
3892 kb_adj=k+1
3893 exit
3894 endif
3895 enddo
3897 ELSEIF(draft == "SH3") then
3898 tunning = 0.6
3899 beta =2.2/tunning
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.)
3903 fzu=1.
3904 do k=1,40
3905 kratio= float(k)/float(40)
3906 zuh2(k) = zubeg+FZU*kratio**(alpha-1.0) * (1.0-kratio)**(beta-1.0)
3907 enddo
3908 if(maxval(zuh2(1:40)).gt.0.) &
3909 zuh2(:)= zuh2(:)/ maxval(zuh2(1:40))
3910 ELSEIF(draft == "MID") then
3911 kb_adj=max(kb,2)
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))
3919 #endif
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)
3923 enddo
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
3928 kb_adj=k+1
3929 exit
3930 endif
3931 enddo
3932 kb_adj=max(2,kb_adj)
3933 do k=kts,kb_adj-1
3934 zu(k)=0.
3935 enddo
3937 ELSEIF(draft == "DOWN" .or. draft == "DOWNM") then
3939 ! tunning = 0.8
3940 ! beta = 3.0/tunning
3941 ! tunning = 0.8
3942 ! beta =2.0/tunning
3943 ! alpha = tunning*beta
3944 ! fzu=1.
3945 ! zuh(:)=0.
3946 tunning=p(kb)
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))
3953 #endif
3954 zu(:)=0.
3955 do k=2,min(kt,ktf)
3956 kratio= (p(k)-p(1))/(p(kt)-p(1))
3957 zu(k) = FZU*kratio**(alpha-1.0) * (1.0-kratio)**(beta-1.0)
3958 enddo
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)
3963 ! enddo
3964 ! do k=maxloc(zuh(:),1)+1,kt
3965 ! kk=kpbli+k-maxloc(zuh(:),1)
3966 ! if(kk.le.kt)zu(kk)=zuh(k)
3967 ! enddo
3968 ! else
3969 ! do k=2,kt ! maxloc(zuh(:),1)
3970 ! zu(k)=zuh(k-1)
3971 ! enddo
3972 ! endif
3973 fzu=maxval(zu(kts:min(ktf,kt+1)))
3974 if(fzu.gt.0.) &
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)
3977 ! do k=2,kt+1
3978 ! zu(k)=fzu*zu(k)
3979 ! enddo
3980 zu(1)=0.
3983 ENDIF
3984 !- normalize ZU
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, &
3992 kbcon,ktop,ierr, &
3993 itf,ktf, &
3994 its,ite, kts,kte )
3996 IMPLICIT NONE
3997 !
3998 ! on input
3999 !
4001 ! only local wrf dimensions are need as of now in this routine
4003 integer &
4004 ,intent (in ) :: &
4005 itf,ktf, &
4006 its,ite, kts,kte
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
4014 !
4015 real, dimension (its:ite,kts:kte) &
4016 ,intent (in ) :: &
4017 z,zu,gamma_cup,t_cup,dby,t,tn,q,qo
4018 integer, dimension (its:ite) &
4019 ,intent (in ) :: &
4020 kbcon,ktop
4021 real, intent(in) :: dtime
4022 !
4023 ! input and output
4024 !
4027 integer, dimension (its:ite) &
4028 ,intent (inout) :: &
4029 ierr
4030 real, dimension (its:ite) &
4031 ,intent (out ) :: &
4032 aa0
4033 !
4034 ! local variables in this routine
4035 !
4037 integer :: &
4038 i,k
4039 real :: &
4040 dz,dA
4041 !
4042 DO i=its,itf
4043 AA0(I)=0.
4044 ENDDO
4045 DO 100 k=kts+1,ktf
4046 DO 100 i=its,itf
4047 IF(ierr(i).ne.0 )GO TO 100
4048 IF(k.gt.KBCON(i))GO TO 100
4050 DZ=Z(I,K)-Z(I,K-1)
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
4058 AA0(I)=AA0(I)+dA
4059 100 CONTINUE
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)
4066 IMPLICIT NONE
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
4074 !-local vars
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
4077 !
4078 !-initialize k_inv_layers as undef
4079 l_mid=300.
4080 l_shal=100.
4081 k_inv_layers(:,:) = 1
4082 do i = its,itf
4083 if(ierr(i) == 0)then
4084 kend_p3(i)=kend(i)+3
4085 DO k = kts+1,kend_p3(i)+4
4086 !- get the 1st der
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)
4089 enddo
4090 DO k = kts+2,kend_p3(i)+3
4091 ! get the 2nd der
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))
4094 enddo
4096 ilev=max(kts+2,kstart(i)+1)
4097 ix=1
4098 k=ilev
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
4105 ix=min(5,ix+1)
4106 ilev=kk+1
4107 exit
4108 endif
4109 ilev=kk+1
4110 enddo
4111 k=ilev
4112 ENDDO
4113 !- 2nd criteria
4114 kadd=0
4115 ken=maxloc(k_inv_layers(i,:),1)
4116 do k=1,ken
4117 kk=k_inv_layers(i,k+kadd)
4118 if(kk.eq.1)exit
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
4122 kadd=kadd+1
4123 do kj = k,ken
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
4126 enddo
4127 endif
4128 ENDDO
4129 endif
4130 ENDDO
4131 100 format(1x,16i3)
4132 !- find the locations of inversions around 800 and 550 hPa
4133 sec_deriv(:)=1.e9
4134 do i = its,itf
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
4141 enddo
4142 k800=minloc(abs(sec_deriv),1)
4143 sec_deriv(:)=1.e9
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
4148 enddo
4149 k550=minloc(abs(sec_deriv),1)
4150 !-save k800 and k550 in k_inv_layers array
4151 shal=1
4152 mid=2
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
4156 ENDDO
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)
4166 ! input ...
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)
4172 ! output ...
4173 ! deriv3 - interpolated value
4174 !============================================================================*/
4176 implicit none
4177 integer, parameter :: n=3
4178 integer ni, m,i, j, k, ix
4179 real:: deriv3, xx
4180 real:: xi(ni), yi(ni), x(n), f(n)
4182 ! exit if too high-order derivative was needed,
4183 if (m > 2) then
4184 deriv3 = 0.0
4185 return
4186 end if
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
4190 deriv3 = 0.0
4191 stop "problems with finding the 2nd derivative"
4192 end if
4194 ! a binary (bisectional) search to find i so that xi(i-1) < x < xi(i)
4195 i = 1
4196 j = ni
4197 do while (j > i+1)
4198 k = (i+j)/2
4199 if (xx < xi(k)) then
4200 j = k
4201 else
4202 i = k
4203 end if
4204 end do
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 ...
4208 i = i + 1 - n/2
4210 ! check boundaries: if i is ouside of the range [1, ... n] -> shift i
4211 if (i < 1) i=1
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
4219 ix = i
4220 ! initialization of f(n) and x(n)
4221 do i=1,n
4222 f(i) = yi(ix+i-1)
4223 x(i) = xi(ix+i-1)
4224 end do
4226 ! calculate the first-order derivative using Lagrange interpolation
4227 if (m == 1) then
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
4232 else
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)))
4236 end if
4237 END FUNCTION DERIV3
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)
4244 Implicit none
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
4255 !-- local vars
4256 Integer :: i,k, incr1,incr2
4257 REAL :: dz,trash,trash2
4259 do k=kts,kte
4260 do i=its,ite
4261 up_massentro(i,k)=0.
4262 up_massdetro(i,k)=0.
4263 up_massentr (i,k)=0.
4264 up_massdetr (i,k)=0.
4265 enddo
4266 enddo
4267 if(present(up_massentru) .and. present(up_massdetru))then
4268 do k=kts,kte
4269 do i=its,ite
4270 up_massentru(i,k)=0.
4271 up_massdetru(i,k)=0.
4272 enddo
4273 enddo
4274 endif
4275 DO i=its,itf
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))
4288 endif
4289 if(zuo(i,k-1).gt.0.)entr_rate_2d(i,k-1)=(up_massentro(i,k-1))/(dz*zuo(i,k-1))
4290 enddo
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))
4300 endif
4302 if(zuo(i,k-1).gt.0.)cd(i,k-1)=up_massdetro(i,k-1)/(dz*zuo(i,k-1))
4303 enddo
4304 up_massdetro(i,ktop(i))=zuo(i,ktop(i))
4305 up_massentro(i,ktop(i))=0.
4306 do k=ktop(i)+1,ktf
4307 cd(i,k)=0.
4308 entr_rate_2d(i,k)=0.
4309 up_massentro(i,k)=0.
4310 up_massdetro(i,k)=0.
4311 enddo
4312 do k=2,ktf-1
4313 up_massentr (i,k-1)=up_massentro(i,k-1)
4314 up_massdetr (i,k-1)=up_massdetro(i,k-1)
4315 enddo
4316 if(present(up_massentru) .and. present(up_massdetru))then
4317 do k=2,ktf-1
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)
4320 enddo
4321 endif
4323 trash=0.
4324 trash2=0.
4325 do k=k22(i)+1,ktop(i)
4326 trash2=trash2+entr_rate_2d(i,k)
4327 enddo
4328 do k=k22(i)+1,kbcon(i)
4329 trash=trash+entr_rate_2d(i,k)
4330 enddo
4332 endif
4333 ENDDO
4334 END SUBROUTINE get_lateral_massflux
4335 !==============================================================================
4336 !----------------------------------------------------------------------
4337 SUBROUTINE gfinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
4338 RUCUTEN,RVCUTEN, &
4339 restart, &
4340 P_QC,P_QI,P_FIRST_SCALAR, &
4341 RTHFTEN, RQVFTEN, &
4342 allowed_to_read, &
4343 ids, ide, jds, jde, kds, kde, &
4344 ims, ime, jms, jme, kms, kme, &
4345 its, ite, jts, jte, kts, kte )
4346 !--------------------------------------------------------------------
4347 IMPLICIT NONE
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) :: &
4356 RTHCUTEN, &
4357 RQVCUTEN, &
4358 RQCCUTEN, &
4359 RQICUTEN
4361 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
4362 RUCUTEN, &
4363 RVCUTEN
4365 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
4366 RTHFTEN, &
4367 RQVFTEN
4369 INTEGER :: i, j, k, itf, jtf, ktf
4370 jtf=min0(jte,jde-1)
4371 ktf=min0(kte,kde-1)
4372 itf=min0(ite,ide-1)
4374 IF(.not.restart)THEN
4375 DO j=jts,jte
4376 DO k=kts,kte
4377 DO i=its,ite
4378 RTHCUTEN(i,k,j)=0.
4379 RQVCUTEN(i,k,j)=0.
4380 RUCUTEN(i,k,j)=0.
4381 RVCUTEN(i,k,j)=0.
4382 ENDDO
4383 ENDDO
4384 ENDDO
4387 DO j=jts,jtf
4388 DO k=kts,ktf
4389 DO i=its,itf
4390 RTHFTEN(i,k,j)=0.
4391 RQVFTEN(i,k,j)=0.
4392 ENDDO
4393 ENDDO
4394 ENDDO
4396 IF (P_QC .ge. P_FIRST_SCALAR) THEN
4397 DO j=jts,jtf
4398 DO k=kts,ktf
4399 DO i=its,itf
4400 RQCCUTEN(i,k,j)=0.
4401 ENDDO
4402 ENDDO
4403 ENDDO
4404 ENDIF
4406 IF (P_QI .ge. P_FIRST_SCALAR) THEN
4407 DO j=jts,jtf
4408 DO k=kts,ktf
4409 DO i=its,itf
4410 RQICUTEN(i,k,j)=0.
4411 ENDDO
4412 ENDDO
4413 ENDDO
4414 ENDIF
4416 ENDIF
4418 END SUBROUTINE gfinit
4419 END MODULE module_cu_gf_deep