module_cu_gd.F
References to this file elsewhere.
1 !WRF:MODEL_LAYER:PHYSICS
2 !
3
4 MODULE module_cu_gd
5
6 CONTAINS
7
8 !-------------------------------------------------------------
9 SUBROUTINE GRELLDRV( &
10 DT,itimestep,DX &
11 ,rho,RAINCV &
12 ,U,V,t,W,q,p,pi &
13 ,dz8w,p8w,XLV,CP,G,r_v &
14 ,STEPCU,htop,hbot &
15 ,CU_ACT_FLAG,warm_rain &
16 ,APR_GR,APR_W,APR_MC,APR_ST,APR_AS &
17 ,APR_CAPMA,APR_CAPME,APR_CAPMI &
18 ,MASS_FLUX,XF_ENS,PR_ENS,HT,XLAND,gsw &
19 ,GDC,GDC2 &
20 ,ensdim,maxiens,maxens,maxens2,maxens3 &
21 ,ids,ide, jds,jde, kds,kde &
22 ,ims,ime, jms,jme, kms,kme &
23 ,its,ite, jts,jte, kts,kte &
24 ,RQVCUTEN,RQCCUTEN,RQICUTEN &
25 ,RQVFTEN,RQVBLTEN &
26 ,RTHFTEN,RTHCUTEN,RTHRATEN,RTHBLTEN &
27 ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
28 )
29 !-------------------------------------------------------------
30 IMPLICIT NONE
31 !-------------------------------------------------------------
32 INTEGER, INTENT(IN ) :: &
33 ids,ide, jds,jde, kds,kde, &
34 ims,ime, jms,jme, kms,kme, &
35 its,ite, jts,jte, kts,kte
36
37 integer, intent (in ) :: &
38 ensdim,maxiens,maxens,maxens2,maxens3
39
40 INTEGER, INTENT(IN ) :: STEPCU, ITIMESTEP
41 LOGICAL, INTENT(IN ) :: warm_rain
42
43 REAL, INTENT(IN ) :: XLV, R_v
44 REAL, INTENT(IN ) :: CP,G
45
46 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
47 INTENT(IN ) :: &
48 U, &
49 V, &
50 W, &
51 pi, &
52 t, &
53 q, &
54 p, &
55 dz8w, &
56 p8w, &
57 rho
58 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
59 OPTIONAL , &
60 INTENT(INOUT ) :: &
61 GDC,GDC2
62
63
64 REAL, DIMENSION( ims:ime , jms:jme ),INTENT(IN) :: GSW,HT,XLAND
65 !
66 REAL, INTENT(IN ) :: DT, DX
67 !
68
69 REAL, DIMENSION( ims:ime , jms:jme ), &
70 INTENT(INOUT) :: RAINCV, MASS_FLUX, &
71 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
72 APR_CAPMA,APR_CAPME,APR_CAPMI,htop,hbot
73 !+lxz
74 ! REAL, DIMENSION( ims:ime , jms:jme ) :: & !, INTENT(INOUT) :: &
75 ! HTOP, &! highest model layer penetrated by cumulus since last reset in radiation_driver
76 ! HBOT ! lowest model layer penetrated by cumulus since last reset in radiation_driver
77 ! ! HBOT>HTOP follow physics leveling convention
78
79 LOGICAL, DIMENSION( ims:ime , jms:jme ), &
80 INTENT(INOUT) :: CU_ACT_FLAG
81
82 !
83 ! Optionals
84 !
85 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
86 OPTIONAL, &
87 INTENT(INOUT) :: RTHFTEN, &
88 RQVFTEN
89
90 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
91 OPTIONAL, &
92 INTENT(IN ) :: &
93 RTHRATEN, &
94 RTHBLTEN, &
95 RQVBLTEN
96 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
97 OPTIONAL, &
98 INTENT(INOUT) :: &
99 RTHCUTEN, &
100 RQVCUTEN, &
101 RQCCUTEN, &
102 RQICUTEN
103 !
104 ! Flags relating to the optional tendency arrays declared above
105 ! Models that carry the optional tendencies will provdide the
106 ! optional arguments at compile time; these flags all the model
107 ! to determine at run-time whether a particular tracer is in
108 ! use or not.
109 !
110 LOGICAL, OPTIONAL :: &
111 F_QV &
112 ,F_QC &
113 ,F_QR &
114 ,F_QI &
115 ,F_QS
116
117
118
119 ! LOCAL VARS
120 real, dimension ( ims:ime , jms:jme , 1:ensdim) :: &
121 massfln,xf_ens,pr_ens
122 real, dimension (its:ite,kts:kte+1) :: &
123 OUTT,OUTQ,OUTQC,phh,cupclw
124 real, dimension (its:ite,kts:kte+1) :: phf
125 real, dimension (its:ite) :: &
126 pret, ter11, aa0, fp
127 !+lxz
128 integer, dimension (its:ite) :: &
129 kbcon, ktop
130 !.lxz
131 integer, dimension (its:ite,jts:jte) :: &
132 iact_old_gr
133 integer :: ichoice,iens,ibeg,iend,jbeg,jend
134
135 !
136 ! basic environmental input includes moisture convergence (mconv)
137 ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
138 ! convection for this call only and at that particular gridpoint
139 !
140 real, dimension (its:ite,kts:kte+1) :: &
141 T2d,TN,q2d,qo,PO,P2d,US,VS,omeg
142 real, dimension (its:ite) :: &
143 Z1,PSUR,AAEQ,direction,mconv,cuten,umean,vmean,pmean
144
145 INTEGER :: i,j,k,ICLDCK,ipr,jpr
146 REAL :: tcrit,dp,dq
147 INTEGER :: itf,jtf,ktf
148 REAL :: rkbcon,rktop !-lxz
149
150 ichoice=0
151 iens=1
152 ipr=0
153 jpr=0
154 jbeg=max(jts,jds+4)
155 jend=min(jte,jde-5)
156 ibeg=max(its,ids+4)
157 iend=min(ite,ide-5)
158 tcrit=258.
159
160 itf=MIN(ite,ide-1)
161 ktf=MIN(kte,kde-1)
162 jtf=MIN(jte,jde-1)
163 !
164 DO 100 J = jts,jtf
165 DO I= its,itf
166 cuten(i)=0.
167 iact_old_gr(i,j)=0
168 mass_flux(i,j)=0.
169 raincv(i,j)=0.
170 CU_ACT_FLAG(i,j) = .true.
171 ENDDO
172 DO k=1,ensdim
173 DO I= its,itf
174 massfln(i,j,k)=0.
175 ENDDO
176 ENDDO
177 #if ( EM_CORE == 1 )
178 DO k= kts,ktf
179 DO I= its,itf
180 RTHFTEN(i,k,j)=(RTHFTEN(i,k,j)+RTHRATEN(i,k,j)+RTHBLTEN(i,k,j))*pi(i,k,j)
181 RQVFTEN(i,k,j)=RQVFTEN(i,k,j)+RQVBLTEN(i,k,j)
182 ENDDO
183 ENDDO
184 ! hydrostatic pressure, first on full levels
185 DO I=ITS,ITF
186 phf(i,1) = p8w(i,1,j)
187 ENDDO
188 ! integrate up, dp = -rho * g * dz
189 DO K=kts+1,ktf+1
190 DO I=ITS,ITF
191 phf(i,k) = phf(i,k-1) - rho(i,k-1,j) * g * dz8w(i,k-1,j)
192 ENDDO
193 ENDDO
194 ! scale factor so that pressure is not zero after integration
195 DO I=ITS,ITF
196 fp(i) = (p8w(i,kts,j)-p8w(i,kte,j))/(phf(i,kts)-phf(i,kte))
197 ENDDO
198 ! re-integrate up, dp = -rho * g * dz * scale_factor
199 DO K=kts+1,ktf+1
200 DO I=ITS,ITF
201 phf(i,k) = phf(i,k-1) - rho(i,k-1,j) * g * dz8w(i,k-1,j) * fp(i)
202 ENDDO
203 ENDDO
204 ! put hydrostatic pressure on half levels
205 DO K=kts,ktf
206 DO I=ITS,ITF
207 phh(i,k) = (phf(i,k) + phf(i,k+1))*0.5
208 ENDDO
209 ENDDO
210 #endif
211 DO I=ITS,ITF
212 #if ( EM_CORE == 1 )
213 PSUR(I)=p8w(I,1,J)*.01
214 #endif
215 #if ( NMM_CORE == 1 )
216 PSUR(I)=p(I,1,J)*.01
217 #endif
218 TER11(I)=HT(i,j)
219 mconv(i)=0.
220 aaeq(i)=0.
221 direction(i)=0.
222 pret(i)=0.
223 umean(i)=0.
224 vmean(i)=0.
225 pmean(i)=0.
226 ENDDO
227 DO K=kts,ktf
228 DO I=ITS,ITF
229 omeg(i,k)=0.
230 ! cupclw(i,k)=0.
231 #if ( EM_CORE == 1 )
232 po(i,k)=phh(i,k)*.01
233 #endif
234
235 #if ( NMM_CORE == 1 )
236 po(i,k)=p(i,k,j)*.01
237 #endif
238 P2d(I,K)=PO(i,k)
239 US(I,K) =u(i,k,j)
240 VS(I,K) =v(i,k,j)
241 T2d(I,K)=t(i,k,j)
242 q2d(I,K)=q(i,k,j)
243 omeg(I,K)= -g*rho(i,k,j)*w(i,k,j)
244 TN(I,K)=t2d(i,k)+RTHFTEN(i,k,j)*dt
245 IF(TN(I,K).LT.200.)TN(I,K)=T2d(I,K)
246 QO(I,K)=q2d(i,k)+RQVFTEN(i,k,j)*dt
247 IF(Q2d(I,K).LT.1.E-08)Q2d(I,K)=1.E-08
248 IF(QO(I,K).LT.1.E-08)QO(I,K)=1.E-08
249 OUTT(I,K)=0.
250 OUTQ(I,K)=0.
251 OUTQC(I,K)=0.
252 ! RTHFTEN(i,k,j)=0.
253 ! RQVFTEN(i,k,j)=0.
254 ENDDO
255 ENDDO
256 do k= kts+1,ktf-1
257 DO I = its,itf
258 if((p2d(i,1)-p2d(i,k)).gt.150.and.p2d(i,k).gt.300)then
259 dp=-.5*(p2d(i,k+1)-p2d(i,k-1))
260 umean(i)=umean(i)+us(i,k)*dp
261 vmean(i)=vmean(i)+vs(i,k)*dp
262 pmean(i)=pmean(i)+dp
263 endif
264 enddo
265 enddo
266 DO I = its,itf
267 umean(i)=umean(i)/pmean(i)
268 vmean(i)=vmean(i)/pmean(i)
269 direction(i)=(atan2(umean(i),vmean(i))+3.1415926)*57.29578
270 if(direction(i).gt.360.)direction(i)=direction(i)-360.
271 ENDDO
272 DO K=kts,ktf-1
273 DO I = its,itf
274 dq=(q2d(i,k+1)-q2d(i,k))
275 mconv(i)=mconv(i)+omeg(i,k)*dq/g
276 ENDDO
277 ENDDO
278 DO I = its,itf
279 if(mconv(i).lt.0.)mconv(i)=0.
280 ENDDO
281 !
282 !---- CALL CUMULUS PARAMETERIZATION
283 !
284 CALL CUP_enss(outqc,j,AAEQ,T2d,Q2d,TER11,TN,QO,PO,PRET, &
285 P2d,OUTT,OUTQ,DT,PSUR,US,VS,tcrit,iens, &
286 mconv,massfln,iact_old_gr,omeg,direction,MASS_FLUX, &
287 maxiens,maxens,maxens2,maxens3,ensdim, &
288 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
289 APR_CAPMA,APR_CAPME,APR_CAPMI,kbcon,ktop, &
290 xf_ens,pr_ens,XLAND,gsw,cupclw, &
291 xlv,r_v,cp,g,ichoice,ipr,jpr, &
292 ids,ide, jds,jde, kds,kde, &
293 ims,ime, jms,jme, kms,kme, &
294 its,ite, jts,jte, kts,kte )
295
296 CALL neg_check(dt,q2d,outq,outt,outqc,pret,its,ite,kts,kte,itf,ktf)
297 if(j.ge.jbeg.and.j.le.jend)then
298
299 DO I=its,itf
300 ! cuten(i)=0.
301 if(i.ge.ibeg.and.i.le.iend)then
302 if(pret(i).gt.0.)then
303 raincv(i,j)=pret(i)*dt
304 cuten(i)=1.
305 rkbcon = kte+kts - kbcon(i)
306 rktop = kte+kts - ktop(i)
307 if (ktop(i) > HTOP(i,j)) HTOP(i,j) = ktop(i)+.001
308 if (kbcon(i) < HBOT(i,j)) HBOT(i,j) = kbcon(i)+.001
309 endif
310 else
311 pret(i)=0.
312 endif
313 ENDDO
314 DO K=kts,ktf
315 DO I=its,itf
316 RTHCUTEN(I,K,J)=outt(i,k)*cuten(i)/pi(i,k,j)
317 RQVCUTEN(I,K,J)=outq(i,k)*cuten(i)
318 ENDDO
319 ENDDO
320
321 IF(PRESENT(RQCCUTEN)) THEN
322 IF ( F_QC ) THEN
323 DO K=kts,ktf
324 DO I=its,itf
325 RQCCUTEN(I,K,J)=outqc(I,K)*cuten(i)
326 IF ( PRESENT( GDC ) ) GDC(I,K,J)=CUPCLW(I,K)*cuten(i)
327 IF ( PRESENT( GDC2 ) ) GDC2(I,K,J)=0.
328 ENDDO
329 ENDDO
330 ENDIF
331 ENDIF
332
333 !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
334
335 IF(PRESENT(RQICUTEN).AND.PRESENT(RQCCUTEN))THEN
336 IF (F_QI) THEN
337 DO K=kts,ktf
338 DO I=its,itf
339 if(t2d(i,k).lt.258.)then
340 RQICUTEN(I,K,J)=outqc(I,K)*cuten(i)
341 RQCCUTEN(I,K,J)=0.
342 IF ( PRESENT( GDC2 ) ) GDC2(I,K,J)=CUPCLW(I,K)*cuten(i)
343 else
344 RQICUTEN(I,K,J)=0.
345 RQCCUTEN(I,K,J)=outqc(I,K)*cuten(i)
346 IF ( PRESENT( GDC ) ) GDC(I,K,J)=CUPCLW(I,K)*cuten(i)
347 endif
348 ENDDO
349 ENDDO
350 ENDIF
351 ENDIF
352 endif !jbeg,jend
353
354 100 continue
355
356 END SUBROUTINE GRELLDRV
357
358
359 SUBROUTINE CUP_enss(OUTQC,J,AAEQ,T,Q,Z1, &
360 TN,QO,PO,PRE,P,OUTT,OUTQ,DTIME,PSUR,US,VS, &
361 TCRIT,iens,mconv,massfln,iact, &
362 omeg,direction,massflx,maxiens, &
363 maxens,maxens2,maxens3,ensdim, &
364 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
365 APR_CAPMA,APR_CAPME,APR_CAPMI,kbcon,ktop, & !-lxz
366 xf_ens,pr_ens,xland,gsw,cupclw, &
367 xl,rv,cp,g,ichoice,ipr,jpr, &
368 ids,ide, jds,jde, kds,kde, &
369 ims,ime, jms,jme, kms,kme, &
370 its,ite, jts,jte, kts,kte )
371
372 IMPLICIT NONE
373
374 integer &
375 ,intent (in ) :: &
376 ids,ide, jds,jde, kds,kde, &
377 ims,ime, jms,jme, kms,kme, &
378 its,ite, jts,jte, kts,kte,ipr,jpr
379 integer, intent (in ) :: &
380 j,ensdim,maxiens,maxens,maxens2,maxens3,ichoice,iens
381 !
382 !
383 !
384 real, dimension (ims:ime,jms:jme,1:ensdim) &
385 ,intent (inout) :: &
386 massfln,xf_ens,pr_ens
387 real, dimension (ims:ime,jms:jme) &
388 ,intent (inout ) :: &
389 APR_GR,APR_W,APR_MC,APR_ST,APR_AS,APR_CAPMA, &
390 APR_CAPME,APR_CAPMI,massflx
391 real, dimension (ims:ime,jms:jme) &
392 ,intent (in ) :: &
393 xland,gsw
394 integer, dimension (ims:ime,jms:jme) &
395 ,intent (in ) :: &
396 iact
397 ! outtem = output temp tendency (per s)
398 ! outq = output q tendency (per s)
399 ! outqc = output qc tendency (per s)
400 ! pre = output precip
401 real, dimension (its:ite,kts:kte) &
402 ,intent (out ) :: &
403 OUTT,OUTQ,OUTQC,CUPCLW
404 real, dimension (its:ite) &
405 ,intent (out ) :: &
406 pre
407 !+lxz
408 integer, dimension (its:ite) &
409 ,intent (out ) :: &
410 kbcon,ktop
411 !.lxz
412 !
413 ! basic environmental input includes moisture convergence (mconv)
414 ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
415 ! convection for this call only and at that particular gridpoint
416 !
417 real, dimension (its:ite,kts:kte) &
418 ,intent (in ) :: &
419 T,TN,PO,P,US,VS,omeg
420 real, dimension (its:ite,kts:kte) &
421 ,intent (inout) :: &
422 Q,QO
423 real, dimension (its:ite) &
424 ,intent (in ) :: &
425 Z1,PSUR,AAEQ,direction,mconv
426
427
428 real &
429 ,intent (in ) :: &
430 dtime,tcrit,xl,cp,rv,g
431
432
433 !
434 ! local ensemble dependent variables in this routine
435 !
436 real, dimension (its:ite,1:maxens) :: &
437 xaa0_ens
438 real, dimension (1:maxens) :: &
439 mbdt_ens
440 real, dimension (1:maxens2) :: &
441 edt_ens
442 real, dimension (its:ite,1:maxens2) :: &
443 edtc
444 real, dimension (its:ite,kts:kte,1:maxens2) :: &
445 dellat_ens,dellaqc_ens,dellaq_ens,pwo_ens
446 !
447 !
448 !
449 !***************** the following are your basic environmental
450 ! variables. They carry a "_cup" if they are
451 ! on model cloud levels (staggered). They carry
452 ! an "o"-ending (z becomes zo), if they are the forced
453 ! variables. They are preceded by x (z becomes xz)
454 ! to indicate modification by some typ of cloud
455 !
456 ! z = heights of model levels
457 ! q = environmental mixing ratio
458 ! qes = environmental saturation mixing ratio
459 ! t = environmental temp
460 ! p = environmental pressure
461 ! he = environmental moist static energy
462 ! hes = environmental saturation moist static energy
463 ! z_cup = heights of model cloud levels
464 ! q_cup = environmental q on model cloud levels
465 ! qes_cup = saturation q on model cloud levels
466 ! t_cup = temperature (Kelvin) on model cloud levels
467 ! p_cup = environmental pressure
468 ! he_cup = moist static energy on model cloud levels
469 ! hes_cup = saturation moist static energy on model cloud levels
470 ! gamma_cup = gamma on model cloud levels
471 !
472 !
473 ! hcd = moist static energy in downdraft
474 ! zd normalized downdraft mass flux
475 ! dby = buoancy term
476 ! entr = entrainment rate
477 ! zd = downdraft normalized mass flux
478 ! entr= entrainment rate
479 ! hcd = h in model cloud
480 ! bu = buoancy term
481 ! zd = normalized downdraft mass flux
482 ! gamma_cup = gamma on model cloud levels
483 ! mentr_rate = entrainment rate
484 ! qcd = cloud q (including liquid water) after entrainment
485 ! qrch = saturation q in cloud
486 ! pwd = evaporate at that level
487 ! pwev = total normalized integrated evaoprate (I2)
488 ! entr= entrainment rate
489 ! z1 = terrain elevation
490 ! entr = downdraft entrainment rate
491 ! jmin = downdraft originating level
492 ! kdet = level above ground where downdraft start detraining
493 ! psur = surface pressure
494 ! z1 = terrain elevation
495 ! pr_ens = precipitation ensemble
496 ! xf_ens = mass flux ensembles
497 ! massfln = downdraft mass flux ensembles used in next timestep
498 ! omeg = omega from large scale model
499 ! mconv = moisture convergence from large scale model
500 ! zd = downdraft normalized mass flux
501 ! zu = updraft normalized mass flux
502 ! dir = "storm motion"
503 ! mbdt = arbitrary numerical parameter
504 ! dtime = dt over which forcing is applied
505 ! iact_gr_old = flag to tell where convection was active
506 ! kbcon = LFC of parcel from k22
507 ! k22 = updraft originating level
508 ! icoic = flag if only want one closure (usually set to zero!)
509 ! dby = buoancy term
510 ! ktop = cloud top (output)
511 ! xmb = total base mass flux
512 ! hc = cloud moist static energy
513 ! hkb = moist static energy at originating level
514 ! mentr_rate = entrainment rate
515
516 real, dimension (its:ite,kts:kte) :: &
517 he,hes,qes,z, &
518 heo,heso,qeso,zo, &
519 xhe,xhes,xqes,xz,xt,xq, &
520
521 qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
522 qeso_cup,qo_cup,heo_cup,heso_cup,zo_cup,po_cup,gammao_cup, &
523 tn_cup, &
524 xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup,xp_cup,xgamma_cup, &
525 xt_cup, &
526
527 dby,qc,qrcd,pwd,pw,hcd,qcd,dbyd,hc,qrc,zu,zd,clw_all, &
528 dbyo,qco,qrcdo,pwdo,pwo,hcdo,qcdo,dbydo,hco,qrco,zuo,zdo, &
529 xdby,xqc,xqrcd,xpwd,xpw,xhcd,xqcd,xhc,xqrc,xzu,xzd, &
530
531 ! cd = detrainment function for updraft
532 ! cdd = detrainment function for downdraft
533 ! dellat = change of temperature per unit mass flux of cloud ensemble
534 ! dellaq = change of q per unit mass flux of cloud ensemble
535 ! dellaqc = change of qc per unit mass flux of cloud ensemble
536
537 cd,cdd,scr1,DELLAH,DELLAQ,DELLAT,DELLAQC
538
539 ! aa0 cloud work function for downdraft
540 ! edt = epsilon
541 ! aa0 = cloud work function without forcing effects
542 ! aa1 = cloud work function with forcing effects
543 ! xaa0 = cloud work function with cloud effects (ensemble dependent)
544 ! edt = epsilon
545 real, dimension (its:ite) :: &
546 edt,edto,edtx,AA1,AA0,XAA0,HKB,HKBO,aad,XHKB,QKB,QKBO, &
547 XMB,XPWAV,XPWEV,PWAV,PWEV,PWAVO,PWEVO,BU,BUO,cap_max,xland1, &
548 cap_max_increment,closure_n
549 integer, dimension (its:ite) :: &
550 kzdown,KDET,K22,KB,JMIN,kstabi,kstabm,K22x, & !-lxz
551 KBCONx,KBx,KTOPx,ierr,ierr2,ierr3,KBMAX
552
553 integer :: &
554 nall,iedt,nens,nens3,ki,I,K,KK,iresult
555 real :: &
556 day,dz,mbdt,entr_rate,radius,entrd_rate,mentr_rate,mentrd_rate, &
557 zcutdown,edtmax,edtmin,depth_min,zkbmax,z_detr,zktop, &
558 massfld,dh,cap_maxs
559
560 integer :: itf,jtf,ktf
561 integer :: jmini
562 logical :: keep_going
563
564 itf=MIN(ite,ide-1)
565 ktf=MIN(kte,kde-1)
566 jtf=MIN(jte,jde-1)
567
568 !sms$distribute end
569 day=86400.
570 do i=its,itf
571 closure_n(i)=16.
572 xland1(i)=1.
573 if(xland(i,j).gt.1.5)xland1(i)=0.
574 cap_max_increment(i)=25.
575 enddo
576 !
577 !--- specify entrainmentrate and detrainmentrate
578 !
579 if(iens.le.4)then
580 radius=14000.-float(iens)*2000.
581 else
582 radius=12000.
583 endif
584 !
585 !--- gross entrainment rate (these may be changed later on in the
586 !--- program, depending what your detrainment is!!)
587 !
588 entr_rate=.2/radius
589 !
590 !--- entrainment of mass
591 !
592 mentrd_rate=0.
593 mentr_rate=entr_rate
594 !
595 !--- initial detrainmentrates
596 !
597 do k=kts,ktf
598 do i=its,itf
599 cupclw(i,k)=0.
600 cd(i,k)=0.1*entr_rate
601 cdd(i,k)=0.
602 enddo
603 enddo
604 !
605 !--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
606 ! base mass flux
607 !
608 edtmax=.8
609 edtmin=.2
610 !
611 !--- minimum depth (m), clouds must have
612 !
613 depth_min=500.
614 !
615 !--- maximum depth (mb) of capping
616 !--- inversion (larger cap = no convection)
617 !
618 cap_maxs=75.
619 !sms$to_local(grid_dh: <1, mix :size>, <2, mjx :size>) begin
620 DO 7 i=its,itf
621 kbmax(i)=1
622 aa0(i)=0.
623 aa1(i)=0.
624 aad(i)=0.
625 edt(i)=0.
626 kstabm(i)=ktf-1
627 IERR(i)=0
628 IERR2(i)=0
629 IERR3(i)=0
630 if(aaeq(i).lt.-1.)then
631 ierr(i)=20
632 endif
633 7 CONTINUE
634 !
635 !--- first check for upstream convection
636 !
637 do i=its,itf
638 cap_max(i)=cap_maxs
639 ! if(tkmax(i,j).lt.2.)cap_max(i)=25.
640 if(gsw(i,j).lt.1.)cap_max(i)=25.
641
642 iresult=0
643 ! massfld=0.
644 ! call cup_direction2(i,j,direction,iact, &
645 ! cu_mfx,iresult,0,massfld, &
646 ! ids,ide, jds,jde, kds,kde, &
647 ! ims,ime, jms,jme, kms,kme, &
648 ! its,ite, jts,jte, kts,kte)
649 ! cap_max(i)=cap_maxs
650 if(iresult.eq.1)then
651 cap_max(i)=cap_maxs+20.
652 endif
653 ! endif
654 enddo
655 !
656 !--- max height(m) above ground where updraft air can originate
657 !
658 zkbmax=4000.
659 !
660 !--- height(m) above which no downdrafts are allowed to originate
661 !
662 zcutdown=3000.
663 !
664 !--- depth(m) over which downdraft detrains all its mass
665 !
666 z_detr=1250.
667 !
668 do nens=1,maxens
669 mbdt_ens(nens)=(float(nens)-3.)*dtime*1.e-3+dtime*5.E-03
670 enddo
671 do nens=1,maxens2
672 edt_ens(nens)=.95-float(nens)*.01
673 enddo
674 ! if(j.eq.jpr)then
675 ! print *,'radius ensemble ',iens,radius
676 ! print *,mbdt_ens
677 ! print *,edt_ens
678 ! endif
679 !
680 !--- environmental conditions, FIRST HEIGHTS
681 !
682 do i=its,itf
683 if(ierr(i).ne.20)then
684 do k=1,maxens*maxens2*maxens3
685 xf_ens(i,j,(iens-1)*maxens*maxens2*maxens3+k)=0.
686 pr_ens(i,j,(iens-1)*maxens*maxens2*maxens3+k)=0.
687 enddo
688 endif
689 enddo
690 !
691 !--- calculate moist static energy, heights, qes
692 !
693 call cup_env(z,qes,he,hes,t,q,p,z1, &
694 psur,ierr,tcrit,0,xl,cp, &
695 ids,ide, jds,jde, kds,kde, &
696 ims,ime, jms,jme, kms,kme, &
697 its,ite, jts,jte, kts,kte)
698 call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
699 psur,ierr,tcrit,0,xl,cp, &
700 ids,ide, jds,jde, kds,kde, &
701 ims,ime, jms,jme, kms,kme, &
702 its,ite, jts,jte, kts,kte)
703 !
704 !--- environmental values on cloud levels
705 !
706 call cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup,he_cup, &
707 hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
708 ierr,z1,xl,rv,cp, &
709 ids,ide, jds,jde, kds,kde, &
710 ims,ime, jms,jme, kms,kme, &
711 its,ite, jts,jte, kts,kte)
712 call cup_env_clev(tn,qeso,qo,heo,heso,zo,po,qeso_cup,qo_cup, &
713 heo_cup,heso_cup,zo_cup,po_cup,gammao_cup,tn_cup,psur, &
714 ierr,z1,xl,rv,cp, &
715 ids,ide, jds,jde, kds,kde, &
716 ims,ime, jms,jme, kms,kme, &
717 its,ite, jts,jte, kts,kte)
718 do i=its,itf
719 if(ierr(i).eq.0)then
720 !
721 do k=kts,ktf-2
722 if(zo_cup(i,k).gt.zkbmax+z1(i))then
723 kbmax(i)=k
724 go to 25
725 endif
726 enddo
727 25 continue
728 !
729 !
730 !--- level where detrainment for downdraft starts
731 !
732 do k=kts,ktf
733 if(zo_cup(i,k).gt.z_detr+z1(i))then
734 kdet(i)=k
735 go to 26
736 endif
737 enddo
738 26 continue
739 !
740 endif
741 enddo
742 !
743 !
744 !
745 !------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
746 !
747 CALL cup_MAXIMI(HEO_CUP,3,KBMAX,K22,ierr, &
748 ids,ide, jds,jde, kds,kde, &
749 ims,ime, jms,jme, kms,kme, &
750 its,ite, jts,jte, kts,kte)
751 DO 36 i=its,itf
752 IF(ierr(I).eq.0.)THEN
753 IF(K22(I).GE.KBMAX(i))ierr(i)=2
754 endif
755 36 CONTINUE
756 !
757 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
758 !
759 call cup_kbcon(cap_max_increment,1,k22,kbcon,heo_cup,heso_cup, &
760 ierr,kbmax,po_cup,cap_max, &
761 ids,ide, jds,jde, kds,kde, &
762 ims,ime, jms,jme, kms,kme, &
763 its,ite, jts,jte, kts,kte)
764 ! call cup_kbcon_cin(1,k22,kbcon,heo_cup,heso_cup,z,tn_cup, &
765 ! qeso_cup,ierr,kbmax,po_cup,cap_max,xl,cp,&
766 ! ids,ide, jds,jde, kds,kde, &
767 ! ims,ime, jms,jme, kms,kme, &
768 ! its,ite, jts,jte, kts,kte)
769 !
770 !--- increase detrainment in stable layers
771 !
772 CALL cup_minimi(HEso_cup,Kbcon,kstabm,kstabi,ierr, &
773 ids,ide, jds,jde, kds,kde, &
774 ims,ime, jms,jme, kms,kme, &
775 its,ite, jts,jte, kts,kte)
776 do i=its,itf
777 IF(ierr(I).eq.0.)THEN
778 if(kstabm(i)-1.gt.kstabi(i))then
779 do k=kstabi(i),kstabm(i)-1
780 cd(i,k)=cd(i,k-1)+1.5*entr_rate
781 if(cd(i,k).gt.10.0*entr_rate)cd(i,k)=10.0*entr_rate
782 enddo
783 ENDIF
784 ENDIF
785 ENDDO
786 !
787 !--- calculate incloud moist static energy
788 !
789 call cup_up_he(k22,hkb,z_cup,cd,mentr_rate,he_cup,hc, &
790 kbcon,ierr,dby,he,hes_cup, &
791 ids,ide, jds,jde, kds,kde, &
792 ims,ime, jms,jme, kms,kme, &
793 its,ite, jts,jte, kts,kte)
794 call cup_up_he(k22,hkbo,zo_cup,cd,mentr_rate,heo_cup,hco, &
795 kbcon,ierr,dbyo,heo,heso_cup, &
796 ids,ide, jds,jde, kds,kde, &
797 ims,ime, jms,jme, kms,kme, &
798 its,ite, jts,jte, kts,kte)
799
800 !--- DETERMINE CLOUD TOP - KTOP
801 !
802 call cup_ktop(1,dbyo,kbcon,ktop,ierr, &
803 ids,ide, jds,jde, kds,kde, &
804 ims,ime, jms,jme, kms,kme, &
805 its,ite, jts,jte, kts,kte)
806 DO 37 i=its,itf
807 kzdown(i)=0
808 if(ierr(i).eq.0)then
809 zktop=(zo_cup(i,ktop(i))-z1(i))*.6
810 zktop=min(zktop+z1(i),zcutdown+z1(i))
811 do k=kts,kte
812 if(zo_cup(i,k).gt.zktop)then
813 kzdown(i)=k
814 go to 37
815 endif
816 enddo
817 endif
818 37 CONTINUE
819 !
820 !--- DOWNDRAFT ORIGINATING LEVEL - JMIN
821 !
822 call cup_minimi(HEso_cup,K22,kzdown,JMIN,ierr, &
823 ids,ide, jds,jde, kds,kde, &
824 ims,ime, jms,jme, kms,kme, &
825 its,ite, jts,jte, kts,kte)
826 DO 100 i=its,ite
827 IF(ierr(I).eq.0.)THEN
828 !
829 !--- check whether it would have buoyancy, if there where
830 !--- no entrainment/detrainment
831 !
832 !jm begin 20061212: the following code replaces code that
833 ! was too complex and causing problem for optimization.
834 ! Done in consultation with G. Grell.
835 jmini = jmin(i)
836 keep_going = .TRUE.
837 DO WHILE ( keep_going )
838 keep_going = .FALSE.
839 if ( jmini - 1 .lt. kdet(i) ) kdet(i) = jmini-1
840 if ( jmini .ge. ktop(i)-1 ) jmini = ktop(i) - 2
841 ki = jmini
842 hcdo(i,ki)=heso_cup(i,ki)
843 DZ=Zo_cup(i,Ki+1)-Zo_cup(i,Ki)
844 dh=0.
845 DO k=ki-1,1,-1
846 hcdo(i,k)=heso_cup(i,jmini)
847 DZ=Zo_cup(i,K+1)-Zo_cup(i,K)
848 dh=dh+dz*(HCDo(i,K)-heso_cup(i,k))
849 IF(dh.gt.0.)THEN
850 jmini=jmini-1
851 IF ( jmini .gt. 3 ) THEN
852 keep_going = .TRUE.
853 ELSE
854 ierr(i) = 9
855 EXIT
856 ENDIF
857 ENDIF
858 ENDDO
859 ENDDO
860 jmin(i) = jmini
861 IF ( jmini .le. 3 ) THEN
862 ierr(i)=4
863 ENDIF
864 !jm end 20061212
865 ENDIF
866 100 CONTINUE
867 !
868 ! - Must have at least depth_min m between cloud convective base
869 ! and cloud top.
870 !
871 do i=its,itf
872 IF(ierr(I).eq.0.)THEN
873 IF(-zo_cup(I,KBCON(I))+zo_cup(I,KTOP(I)).LT.depth_min)then
874 ierr(i)=6
875 endif
876 endif
877 enddo
878
879 !
880 !c--- normalized updraft mass flux profile
881 !
882 call cup_up_nms(zu,z_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
883 ids,ide, jds,jde, kds,kde, &
884 ims,ime, jms,jme, kms,kme, &
885 its,ite, jts,jte, kts,kte)
886 call cup_up_nms(zuo,zo_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
887 ids,ide, jds,jde, kds,kde, &
888 ims,ime, jms,jme, kms,kme, &
889 its,ite, jts,jte, kts,kte)
890 !
891 !c--- normalized downdraft mass flux profile,also work on bottom detrainment
892 !--- in this routine
893 !
894 call cup_dd_nms(zd,z_cup,cdd,mentrd_rate,jmin,ierr, &
895 0,kdet,z1, &
896 ids,ide, jds,jde, kds,kde, &
897 ims,ime, jms,jme, kms,kme, &
898 its,ite, jts,jte, kts,kte)
899 call cup_dd_nms(zdo,zo_cup,cdd,mentrd_rate,jmin,ierr, &
900 1,kdet,z1, &
901 ids,ide, jds,jde, kds,kde, &
902 ims,ime, jms,jme, kms,kme, &
903 its,ite, jts,jte, kts,kte)
904 !
905 !--- downdraft moist static energy
906 !
907 call cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,mentrd_rate, &
908 jmin,ierr,he,dbyd,he_cup, &
909 ids,ide, jds,jde, kds,kde, &
910 ims,ime, jms,jme, kms,kme, &
911 its,ite, jts,jte, kts,kte)
912 call cup_dd_he(heso_cup,zdo,hcdo,zo_cup,cdd,mentrd_rate, &
913 jmin,ierr,heo,dbydo,he_cup,&
914 ids,ide, jds,jde, kds,kde, &
915 ims,ime, jms,jme, kms,kme, &
916 its,ite, jts,jte, kts,kte)
917 !
918 !--- calculate moisture properties of downdraft
919 !
920 call cup_dd_moisture(zd,hcd,hes_cup,qcd,qes_cup, &
921 pwd,q_cup,z_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
922 pwev,bu,qrcd,q,he,t_cup,2,xl, &
923 ids,ide, jds,jde, kds,kde, &
924 ims,ime, jms,jme, kms,kme, &
925 its,ite, jts,jte, kts,kte)
926 call cup_dd_moisture(zdo,hcdo,heso_cup,qcdo,qeso_cup, &
927 pwdo,qo_cup,zo_cup,cdd,mentrd_rate,jmin,ierr,gammao_cup, &
928 pwevo,bu,qrcdo,qo,heo,tn_cup,1,xl, &
929 ids,ide, jds,jde, kds,kde, &
930 ims,ime, jms,jme, kms,kme, &
931 its,ite, jts,jte, kts,kte)
932 !
933 !--- calculate moisture properties of updraft
934 !
935 call cup_up_moisture(ierr,z_cup,qc,qrc,pw,pwav, &
936 kbcon,ktop,cd,dby,mentr_rate,clw_all, &
937 q,GAMMA_cup,zu,qes_cup,k22,q_cup,xl, &
938 ids,ide, jds,jde, kds,kde, &
939 ims,ime, jms,jme, kms,kme, &
940 its,ite, jts,jte, kts,kte)
941 do k=kts,ktf
942 do i=its,itf
943 cupclw(i,k)=qrc(i,k)
944 enddo
945 enddo
946 call cup_up_moisture(ierr,zo_cup,qco,qrco,pwo,pwavo, &
947 kbcon,ktop,cd,dbyo,mentr_rate,clw_all, &
948 qo,GAMMAo_cup,zuo,qeso_cup,k22,qo_cup,xl,&
949 ids,ide, jds,jde, kds,kde, &
950 ims,ime, jms,jme, kms,kme, &
951 its,ite, jts,jte, kts,kte)
952 !
953 !--- calculate workfunctions for updrafts
954 !
955 call cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
956 kbcon,ktop,ierr, &
957 ids,ide, jds,jde, kds,kde, &
958 ims,ime, jms,jme, kms,kme, &
959 its,ite, jts,jte, kts,kte)
960 call cup_up_aa0(aa1,zo,zuo,dbyo,GAMMAo_CUP,tn_cup, &
961 kbcon,ktop,ierr, &
962 ids,ide, jds,jde, kds,kde, &
963 ims,ime, jms,jme, kms,kme, &
964 its,ite, jts,jte, kts,kte)
965 do i=its,itf
966 if(ierr(i).eq.0)then
967 if(aa1(i).eq.0.)then
968 ierr(i)=17
969 endif
970 endif
971 enddo
972 !
973 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
974 !
975 call cup_dd_edt(ierr,us,vs,zo,ktop,kbcon,edt,po,pwavo, &
976 pwevo,edtmax,edtmin,maxens2,edtc, &
977 ids,ide, jds,jde, kds,kde, &
978 ims,ime, jms,jme, kms,kme, &
979 its,ite, jts,jte, kts,kte)
980 do 250 iedt=1,maxens2
981 do i=its,itf
982 if(ierr(i).eq.0)then
983 edt(i)=edtc(i,iedt)
984 edto(i)=edtc(i,iedt)
985 edtx(i)=edtc(i,iedt)
986 endif
987 enddo
988 do k=kts,ktf
989 do i=its,itf
990 dellat_ens(i,k,iedt)=0.
991 dellaq_ens(i,k,iedt)=0.
992 dellaqc_ens(i,k,iedt)=0.
993 pwo_ens(i,k,iedt)=0.
994 enddo
995 enddo
996 !
997 do i=its,itf
998 aad(i)=0.
999 enddo
1000 ! do i=its,itf
1001 ! if(ierr(i).eq.0)then
1002 ! eddt(i,j)=edt(i)
1003 ! EDTX(I)=EDT(I)
1004 ! BU(I)=0.
1005 ! BUO(I)=0.
1006 ! endif
1007 ! enddo
1008 !
1009 !---downdraft workfunctions
1010 !
1011 ! call cup_dd_aa0(edt,ierr,aa0,jmin,gamma_cup,t_cup, &
1012 ! hcd,hes_cup,z,zd, &
1013 ! ids,ide, jds,jde, kds,kde, &
1014 ! ims,ime, jms,jme, kms,kme, &
1015 ! its,ite, jts,jte, kts,kte)
1016 ! call cup_dd_aa0(edto,ierr,aad,jmin,gammao_cup,tn_cup, &
1017 ! hcdo,heso_cup,zo,zdo, &
1018 ! ids,ide, jds,jde, kds,kde, &
1019 ! ims,ime, jms,jme, kms,kme, &
1020 ! its,ite, jts,jte, kts,kte)
1021 !
1022 !--- change per unit mass that a model cloud would modify the environment
1023 !
1024 !--- 1. in bottom layer
1025 !
1026 call cup_dellabot(ipr,jpr,heo_cup,ierr,zo_cup,po,hcdo,edto, &
1027 zdo,cdd,heo,dellah,j,mentrd_rate,zo,g, &
1028 ids,ide, jds,jde, kds,kde, &
1029 ims,ime, jms,jme, kms,kme, &
1030 its,ite, jts,jte, kts,kte)
1031 call cup_dellabot(ipr,jpr,qo_cup,ierr,zo_cup,po,qrcdo,edto, &
1032 zdo,cdd,qo,dellaq,j,mentrd_rate,zo,g,&
1033 ids,ide, jds,jde, kds,kde, &
1034 ims,ime, jms,jme, kms,kme, &
1035 its,ite, jts,jte, kts,kte)
1036 !
1037 !--- 2. everywhere else
1038 !
1039 call cup_dellas(ierr,zo_cup,po_cup,hcdo,edto,zdo,cdd, &
1040 heo,dellah,j,mentrd_rate,zuo,g, &
1041 cd,hco,ktop,k22,kbcon,mentr_rate,jmin,heo_cup,kdet, &
1042 k22,ipr,jpr,'deep', &
1043 ids,ide, jds,jde, kds,kde, &
1044 ims,ime, jms,jme, kms,kme, &
1045 its,ite, jts,jte, kts,kte)
1046 !
1047 !-- take out cloud liquid water for detrainment
1048 !
1049 !?? do k=kts,ktf
1050 do k=kts,ktf-1
1051 do i=its,itf
1052 scr1(i,k)=0.
1053 dellaqc(i,k)=0.
1054 if(ierr(i).eq.0)then
1055 ! print *,'in vupnewg, after della ',ierr(i),aa0(i),i,j
1056 scr1(i,k)=qco(i,k)-qrco(i,k)
1057 if(k.eq.ktop(i)-0)dellaqc(i,k)= &
1058 .01*zuo(i,ktop(i))*qrco(i,ktop(i))* &
1059 9.81/(po_cup(i,k)-po_cup(i,k+1))
1060 if(k.lt.ktop(i).and.k.gt.kbcon(i))then
1061 dz=zo_cup(i,k+1)-zo_cup(i,k)
1062 dellaqc(i,k)=.01*9.81*cd(i,k)*dz*zuo(i,k) &
1063 *.5*(qrco(i,k)+qrco(i,k+1))/ &
1064 (po_cup(i,k)-po_cup(i,k+1))
1065 endif
1066 endif
1067 enddo
1068 enddo
1069 call cup_dellas(ierr,zo_cup,po_cup,qrcdo,edto,zdo,cdd, &
1070 qo,dellaq,j,mentrd_rate,zuo,g, &
1071 cd,scr1,ktop,k22,kbcon,mentr_rate,jmin,qo_cup,kdet, &
1072 k22,ipr,jpr,'deep', &
1073 ids,ide, jds,jde, kds,kde, &
1074 ims,ime, jms,jme, kms,kme, &
1075 its,ite, jts,jte, kts,kte )
1076 !
1077 !--- using dellas, calculate changed environmental profiles
1078 !
1079 ! do 200 nens=1,maxens
1080 mbdt=mbdt_ens(2)
1081 do i=its,itf
1082 xaa0_ens(i,1)=0.
1083 xaa0_ens(i,2)=0.
1084 xaa0_ens(i,3)=0.
1085 enddo
1086
1087 ! mbdt=mbdt_ens(nens)
1088 ! do i=its,itf
1089 ! xaa0_ens(i,nens)=0.
1090 ! enddo
1091 do k=kts,ktf
1092 do i=its,itf
1093 dellat(i,k)=0.
1094 if(ierr(i).eq.0)then
1095 XHE(I,K)=DELLAH(I,K)*MBDT+HEO(I,K)
1096 XQ(I,K)=DELLAQ(I,K)*MBDT+QO(I,K)
1097 DELLAT(I,K)=(1./cp)*(DELLAH(I,K)-xl*DELLAQ(I,K))
1098 XT(I,K)= DELLAT(I,K)*MBDT+TN(I,K)
1099 IF(XQ(I,K).LE.0.)XQ(I,K)=1.E-08
1100 ! if(i.eq.ipr.and.j.eq.jpr)then
1101 ! print *,k,DELLAH(I,K),DELLAQ(I,K),DELLAT(I,K)
1102 ! endif
1103 ENDIF
1104 enddo
1105 enddo
1106 do i=its,itf
1107 if(ierr(i).eq.0)then
1108 XHE(I,ktf)=HEO(I,ktf)
1109 XQ(I,ktf)=QO(I,ktf)
1110 XT(I,ktf)=TN(I,ktf)
1111 IF(XQ(I,ktf).LE.0.)XQ(I,ktf)=1.E-08
1112 endif
1113 enddo
1114 !
1115 !--- calculate moist static energy, heights, qes
1116 !
1117 call cup_env(xz,xqes,xhe,xhes,xt,xq,po,z1, &
1118 psur,ierr,tcrit,2,xl,cp, &
1119 ids,ide, jds,jde, kds,kde, &
1120 ims,ime, jms,jme, kms,kme, &
1121 its,ite, jts,jte, kts,kte)
1122 !
1123 !--- environmental values on cloud levels
1124 !
1125 call cup_env_clev(xt,xqes,xq,xhe,xhes,xz,po,xqes_cup,xq_cup, &
1126 xhe_cup,xhes_cup,xz_cup,po_cup,gamma_cup,xt_cup,psur, &
1127 ierr,z1,xl,rv,cp, &
1128 ids,ide, jds,jde, kds,kde, &
1129 ims,ime, jms,jme, kms,kme, &
1130 its,ite, jts,jte, kts,kte)
1131 !
1132 !
1133 !**************************** static control
1134 !
1135 !--- moist static energy inside cloud
1136 !
1137 do i=its,itf
1138 if(ierr(i).eq.0)then
1139 xhkb(i)=xhe(i,k22(i))
1140 endif
1141 enddo
1142 call cup_up_he(k22,xhkb,xz_cup,cd,mentr_rate,xhe_cup,xhc, &
1143 kbcon,ierr,xdby,xhe,xhes_cup, &
1144 ids,ide, jds,jde, kds,kde, &
1145 ims,ime, jms,jme, kms,kme, &
1146 its,ite, jts,jte, kts,kte)
1147 !
1148 !c--- normalized mass flux profile
1149 !
1150 call cup_up_nms(xzu,xz_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
1151 ids,ide, jds,jde, kds,kde, &
1152 ims,ime, jms,jme, kms,kme, &
1153 its,ite, jts,jte, kts,kte)
1154 !
1155 !--- moisture downdraft
1156 !
1157 call cup_dd_nms(xzd,xz_cup,cdd,mentrd_rate,jmin,ierr, &
1158 1,kdet,z1, &
1159 ids,ide, jds,jde, kds,kde, &
1160 ims,ime, jms,jme, kms,kme, &
1161 its,ite, jts,jte, kts,kte)
1162 call cup_dd_he(xhes_cup,xzd,xhcd,xz_cup,cdd,mentrd_rate, &
1163 jmin,ierr,xhe,dbyd,xhe_cup,&
1164 ids,ide, jds,jde, kds,kde, &
1165 ims,ime, jms,jme, kms,kme, &
1166 its,ite, jts,jte, kts,kte)
1167 call cup_dd_moisture(xzd,xhcd,xhes_cup,xqcd,xqes_cup, &
1168 xpwd,xq_cup,xz_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
1169 xpwev,bu,xqrcd,xq,xhe,xt_cup,3,xl, &
1170 ids,ide, jds,jde, kds,kde, &
1171 ims,ime, jms,jme, kms,kme, &
1172 its,ite, jts,jte, kts,kte)
1173
1174 !
1175 !------- MOISTURE updraft
1176 !
1177 call cup_up_moisture(ierr,xz_cup,xqc,xqrc,xpw,xpwav, &
1178 kbcon,ktop,cd,xdby,mentr_rate,clw_all, &
1179 xq,GAMMA_cup,xzu,xqes_cup,k22,xq_cup,xl, &
1180 ids,ide, jds,jde, kds,kde, &
1181 ims,ime, jms,jme, kms,kme, &
1182 its,ite, jts,jte, kts,kte)
1183 !
1184 !--- workfunctions for updraft
1185 !
1186 call cup_up_aa0(xaa0,xz,xzu,xdby,GAMMA_CUP,xt_cup, &
1187 kbcon,ktop,ierr, &
1188 ids,ide, jds,jde, kds,kde, &
1189 ims,ime, jms,jme, kms,kme, &
1190 its,ite, jts,jte, kts,kte)
1191 !
1192 !--- workfunctions for downdraft
1193 !
1194 !
1195 ! call cup_dd_aa0(edtx,ierr,xaa0,jmin,gamma_cup,xt_cup, &
1196 ! xhcd,xhes_cup,xz,xzd, &
1197 ! ids,ide, jds,jde, kds,kde, &
1198 ! ims,ime, jms,jme, kms,kme, &
1199 ! its,ite, jts,jte, kts,kte)
1200 do 200 nens=1,maxens
1201 do i=its,itf
1202 if(ierr(i).eq.0)then
1203 xaa0_ens(i,nens)=xaa0(i)
1204 nall=(iens-1)*maxens3*maxens*maxens2 &
1205 +(iedt-1)*maxens*maxens3 &
1206 +(nens-1)*maxens3
1207 do k=kts,ktf
1208 if(k.le.ktop(i))then
1209 do nens3=1,maxens3
1210 if(nens3.eq.7)then
1211 !--- b=0
1212 pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
1213 pwo(i,k)
1214 ! +edto(i)*pwdo(i,k)
1215 !--- b=beta
1216 else if(nens3.eq.8)then
1217 pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
1218 pwo(i,k)
1219 !--- b=beta/2
1220 else if(nens3.eq.9)then
1221 pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
1222 pwo(i,k)
1223 ! +.5*edto(i)*pwdo(i,k)
1224 else
1225 pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
1226 pwo(i,k)+edto(i)*pwdo(i,k)
1227 endif
1228 enddo
1229 endif
1230 enddo
1231 if(pr_ens(i,j,nall+7).lt.1.e-6)then
1232 ierr(i)=18
1233 do nens3=1,maxens3
1234 pr_ens(i,j,nall+nens3)=0.
1235 enddo
1236 endif
1237 do nens3=1,maxens3
1238 if(pr_ens(i,j,nall+nens3).lt.1.e-4)then
1239 pr_ens(i,j,nall+nens3)=0.
1240 endif
1241 enddo
1242 endif
1243 enddo
1244 200 continue
1245 !
1246 !--- LARGE SCALE FORCING
1247 !
1248 !
1249 !------- CHECK wether aa0 should have been zero
1250 !
1251 !
1252 CALL cup_MAXIMI(HEO_CUP,3,KBMAX,K22x,ierr, &
1253 ids,ide, jds,jde, kds,kde, &
1254 ims,ime, jms,jme, kms,kme, &
1255 its,ite, jts,jte, kts,kte)
1256 do i=its,itf
1257 ierr2(i)=ierr(i)
1258 ierr3(i)=ierr(i)
1259 enddo
1260 call cup_kbcon(cap_max_increment,2,k22x,kbconx,heo_cup, &
1261 heso_cup,ierr2,kbmax,po_cup,cap_max, &
1262 ids,ide, jds,jde, kds,kde, &
1263 ims,ime, jms,jme, kms,kme, &
1264 its,ite, jts,jte, kts,kte)
1265 call cup_kbcon(cap_max_increment,3,k22x,kbconx,heo_cup, &
1266 heso_cup,ierr3,kbmax,po_cup,cap_max, &
1267 ids,ide, jds,jde, kds,kde, &
1268 ims,ime, jms,jme, kms,kme, &
1269 its,ite, jts,jte, kts,kte)
1270 !
1271 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
1272 !
1273 call cup_forcing_ens(closure_n,xland1,aa0,aa1,xaa0_ens,mbdt_ens,dtime, &
1274 ierr,ierr2,ierr3,xf_ens,j,'deeps', &
1275 maxens,iens,iedt,maxens2,maxens3,mconv, &
1276 po_cup,ktop,omeg,zdo,k22,zuo,pr_ens,edto,kbcon, &
1277 massflx,iact,direction,ensdim,massfln,ichoice, &
1278 ids,ide, jds,jde, kds,kde, &
1279 ims,ime, jms,jme, kms,kme, &
1280 its,ite, jts,jte, kts,kte)
1281 !
1282 do k=kts,ktf
1283 do i=its,itf
1284 if(ierr(i).eq.0)then
1285 dellat_ens(i,k,iedt)=dellat(i,k)
1286 dellaq_ens(i,k,iedt)=dellaq(i,k)
1287 dellaqc_ens(i,k,iedt)=dellaqc(i,k)
1288 pwo_ens(i,k,iedt)=pwo(i,k)+edt(i)*pwdo(i,k)
1289 else
1290 dellat_ens(i,k,iedt)=0.
1291 dellaq_ens(i,k,iedt)=0.
1292 dellaqc_ens(i,k,iedt)=0.
1293 pwo_ens(i,k,iedt)=0.
1294 endif
1295 ! if(i.eq.ipr.and.j.eq.jpr)then
1296 ! print *,iens,iedt,dellat(i,k),dellat_ens(i,k,iedt), &
1297 ! dellaq(i,k), dellaqc(i,k)
1298 ! endif
1299 enddo
1300 enddo
1301 250 continue
1302 !
1303 !--- FEEDBACK
1304 !
1305 call cup_output_ens(xf_ens,ierr,dellat_ens,dellaq_ens, &
1306 dellaqc_ens,outt,outq,outqc,pre,pwo_ens,xmb,ktop, &
1307 j,'deep',maxens2,maxens,iens,ierr2,ierr3, &
1308 pr_ens,maxens3,ensdim,massfln, &
1309 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
1310 APR_CAPMA,APR_CAPME,APR_CAPMI,closure_n,xland1, &
1311 ids,ide, jds,jde, kds,kde, &
1312 ims,ime, jms,jme, kms,kme, &
1313 its,ite, jts,jte, kts,kte)
1314 do i=its,itf
1315 PRE(I)=MAX(PRE(I),0.)
1316 enddo
1317 !
1318 !---------------------------done------------------------------
1319 !
1320
1321 END SUBROUTINE CUP_enss
1322
1323
1324 SUBROUTINE cup_dd_aa0(edt,ierr,aa0,jmin,gamma_cup,t_cup, &
1325 hcd,hes_cup,z,zd, &
1326 ids,ide, jds,jde, kds,kde, &
1327 ims,ime, jms,jme, kms,kme, &
1328 its,ite, jts,jte, kts,kte )
1329
1330 IMPLICIT NONE
1331 !
1332 ! on input
1333 !
1334
1335 ! only local wrf dimensions are need as of now in this routine
1336
1337 integer &
1338 ,intent (in ) :: &
1339 ids,ide, jds,jde, kds,kde, &
1340 ims,ime, jms,jme, kms,kme, &
1341 its,ite, jts,jte, kts,kte
1342 ! aa0 cloud work function for downdraft
1343 ! gamma_cup = gamma on model cloud levels
1344 ! t_cup = temperature (Kelvin) on model cloud levels
1345 ! hes_cup = saturation moist static energy on model cloud levels
1346 ! hcd = moist static energy in downdraft
1347 ! edt = epsilon
1348 ! zd normalized downdraft mass flux
1349 ! z = heights of model levels
1350 ! ierr error value, maybe modified in this routine
1351 !
1352 real, dimension (its:ite,kts:kte) &
1353 ,intent (in ) :: &
1354 z,zd,gamma_cup,t_cup,hes_cup,hcd
1355 real, dimension (its:ite) &
1356 ,intent (in ) :: &
1357 edt
1358 integer, dimension (its:ite) &
1359 ,intent (in ) :: &
1360 jmin
1361 !
1362 ! input and output
1363 !
1364
1365
1366 integer, dimension (its:ite) &
1367 ,intent (inout) :: &
1368 ierr
1369 real, dimension (its:ite) &
1370 ,intent (out ) :: &
1371 aa0
1372 !
1373 ! local variables in this routine
1374 !
1375
1376 integer :: &
1377 i,k,kk
1378 real :: &
1379 dz
1380 !
1381 integer :: itf, ktf
1382 !
1383 itf=MIN(ite,ide-1)
1384 ktf=MIN(kte,kde-1)
1385 !
1386 !?? DO k=kts,kte-1
1387 DO k=kts,ktf-1
1388 do i=its,itf
1389 IF(ierr(I).eq.0.and.k.lt.jmin(i))then
1390 KK=JMIN(I)-K
1391 !
1392 !--- ORIGINAL
1393 !
1394 DZ=(Z(I,KK)-Z(I,KK+1))
1395 AA0(I)=AA0(I)+zd(i,kk)*EDT(I)*DZ*(9.81/(1004.*T_cup(I,KK))) &
1396 *((hcd(i,kk)-hes_cup(i,kk))/(1.+GAMMA_cup(i,kk)))
1397 endif
1398 enddo
1399 enddo
1400
1401 END SUBROUTINE CUP_dd_aa0
1402
1403
1404 SUBROUTINE cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, &
1405 pwev,edtmax,edtmin,maxens2,edtc, &
1406 ids,ide, jds,jde, kds,kde, &
1407 ims,ime, jms,jme, kms,kme, &
1408 its,ite, jts,jte, kts,kte )
1409
1410 IMPLICIT NONE
1411
1412 integer &
1413 ,intent (in ) :: &
1414 ids,ide, jds,jde, kds,kde, &
1415 ims,ime, jms,jme, kms,kme, &
1416 its,ite, jts,jte, kts,kte
1417 integer, intent (in ) :: &
1418 maxens2
1419 !
1420 ! ierr error value, maybe modified in this routine
1421 !
1422 real, dimension (its:ite,kts:kte) &
1423 ,intent (in ) :: &
1424 us,vs,z,p
1425 real, dimension (its:ite,1:maxens2) &
1426 ,intent (out ) :: &
1427 edtc
1428 real, dimension (its:ite) &
1429 ,intent (out ) :: &
1430 edt
1431 real, dimension (its:ite) &
1432 ,intent (in ) :: &
1433 pwav,pwev
1434 real &
1435 ,intent (in ) :: &
1436 edtmax,edtmin
1437 integer, dimension (its:ite) &
1438 ,intent (in ) :: &
1439 ktop,kbcon
1440 integer, dimension (its:ite) &
1441 ,intent (inout) :: &
1442 ierr
1443 !
1444 ! local variables in this routine
1445 !
1446
1447 integer i,k,kk
1448 real einc,pef,pefb,prezk,zkbc
1449 real, dimension (its:ite) :: &
1450 vshear,sdp,vws
1451
1452 integer :: itf, ktf
1453
1454 itf=MIN(ite,ide-1)
1455 ktf=MIN(kte,kde-1)
1456 !
1457 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
1458 !
1459 ! */ calculate an average wind shear over the depth of the cloud
1460 !
1461 do i=its,itf
1462 edt(i)=0.
1463 vws(i)=0.
1464 sdp(i)=0.
1465 vshear(i)=0.
1466 enddo
1467 do kk = kts,ktf-1
1468 do 62 i=its,itf
1469 IF(ierr(i).ne.0)GO TO 62
1470 if (kk .le. min0(ktop(i),ktf) .and. kk .ge. kbcon(i)) then
1471 vws(i) = vws(i)+ &
1472 (abs((us(i,kk+1)-us(i,kk))/(z(i,kk+1)-z(i,kk))) &
1473 + abs((vs(i,kk+1)-vs(i,kk))/(z(i,kk+1)-z(i,kk)))) * &
1474 (p(i,kk) - p(i,kk+1))
1475 sdp(i) = sdp(i) + p(i,kk) - p(i,kk+1)
1476 endif
1477 if (kk .eq. ktf)vshear(i) = 1.e3 * vws(i) / sdp(i)
1478 62 continue
1479 end do
1480 do i=its,itf
1481 IF(ierr(i).eq.0)then
1482 pef=(1.591-.639*VSHEAR(I)+.0953*(VSHEAR(I)**2) &
1483 -.00496*(VSHEAR(I)**3))
1484 if(pef.gt.edtmax)pef=edtmax
1485 if(pef.lt.edtmin)pef=edtmin
1486 !
1487 !--- cloud base precip efficiency
1488 !
1489 zkbc=z(i,kbcon(i))*3.281e-3
1490 prezk=.02
1491 if(zkbc.gt.3.)then
1492 prezk=.96729352+zkbc*(-.70034167+zkbc*(.162179896+zkbc &
1493 *(- 1.2569798E-2+zkbc*(4.2772E-4-zkbc*5.44E-6))))
1494 endif
1495 if(zkbc.gt.25)then
1496 prezk=2.4
1497 endif
1498 pefb=1./(1.+prezk)
1499 if(pefb.gt.edtmax)pefb=edtmax
1500 if(pefb.lt.edtmin)pefb=edtmin
1501 EDT(I)=1.-.5*(pefb+pef)
1502 !--- edt here is 1-precipeff!
1503 ! einc=(1.-edt(i))/float(maxens2)
1504 ! einc=edt(i)/float(maxens2+1)
1505 !--- 20 percent
1506 einc=.2*edt(i)
1507 do k=1,maxens2
1508 edtc(i,k)=edt(i)+float(k-2)*einc
1509 enddo
1510 endif
1511 enddo
1512 do i=its,itf
1513 IF(ierr(i).eq.0)then
1514 do k=1,maxens2
1515 EDTC(I,K)=-EDTC(I,K)*PWAV(I)/PWEV(I)
1516 IF(EDTC(I,K).GT.edtmax)EDTC(I,K)=edtmax
1517 IF(EDTC(I,K).LT.edtmin)EDTC(I,K)=edtmin
1518 enddo
1519 endif
1520 enddo
1521
1522 END SUBROUTINE cup_dd_edt
1523
1524
1525 SUBROUTINE cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,entr, &
1526 jmin,ierr,he,dby,he_cup, &
1527 ids,ide, jds,jde, kds,kde, &
1528 ims,ime, jms,jme, kms,kme, &
1529 its,ite, jts,jte, kts,kte )
1530
1531 IMPLICIT NONE
1532 !
1533 ! on input
1534 !
1535
1536 ! only local wrf dimensions are need as of now in this routine
1537
1538 integer &
1539 ,intent (in ) :: &
1540 ids,ide, jds,jde, kds,kde, &
1541 ims,ime, jms,jme, kms,kme, &
1542 its,ite, jts,jte, kts,kte
1543 ! hcd = downdraft moist static energy
1544 ! he = moist static energy on model levels
1545 ! he_cup = moist static energy on model cloud levels
1546 ! hes_cup = saturation moist static energy on model cloud levels
1547 ! dby = buoancy term
1548 ! cdd= detrainment function
1549 ! z_cup = heights of model cloud levels
1550 ! entr = entrainment rate
1551 ! zd = downdraft normalized mass flux
1552 !
1553 real, dimension (its:ite,kts:kte) &
1554 ,intent (in ) :: &
1555 he,he_cup,hes_cup,z_cup,cdd,zd
1556 ! entr= entrainment rate
1557 real &
1558 ,intent (in ) :: &
1559 entr
1560 integer, dimension (its:ite) &
1561 ,intent (in ) :: &
1562 jmin
1563 !
1564 ! input and output
1565 !
1566
1567 ! ierr error value, maybe modified in this routine
1568
1569 integer, dimension (its:ite) &
1570 ,intent (inout) :: &
1571 ierr
1572
1573 real, dimension (its:ite,kts:kte) &
1574 ,intent (out ) :: &
1575 hcd,dby
1576 !
1577 ! local variables in this routine
1578 !
1579
1580 integer :: &
1581 i,k,ki
1582 real :: &
1583 dz
1584
1585 integer :: itf, ktf
1586
1587 itf=MIN(ite,ide-1)
1588 ktf=MIN(kte,kde-1)
1589
1590 do k=kts+1,ktf
1591 do i=its,itf
1592 dby(i,k)=0.
1593 IF(ierr(I).eq.0)then
1594 hcd(i,k)=hes_cup(i,k)
1595 endif
1596 enddo
1597 enddo
1598 !
1599 do 100 i=its,itf
1600 IF(ierr(I).eq.0)then
1601 k=jmin(i)
1602 hcd(i,k)=hes_cup(i,k)
1603 dby(i,k)=hcd(i,jmin(i))-hes_cup(i,k)
1604 !
1605 do ki=jmin(i)-1,1,-1
1606 DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
1607 HCD(i,Ki)=(HCD(i,Ki+1)*(1.-.5*CDD(i,Ki)*DZ) &
1608 +entr*DZ*HE(i,Ki) &
1609 )/(1.+entr*DZ-.5*CDD(i,Ki)*DZ)
1610 dby(i,ki)=HCD(i,Ki)-hes_cup(i,ki)
1611 enddo
1612 !
1613 endif
1614 !--- end loop over i
1615 100 continue
1616
1617
1618 END SUBROUTINE cup_dd_he
1619
1620
1621 SUBROUTINE cup_dd_moisture(zd,hcd,hes_cup,qcd,qes_cup, &
1622 pwd,q_cup,z_cup,cdd,entr,jmin,ierr, &
1623 gamma_cup,pwev,bu,qrcd, &
1624 q,he,t_cup,iloop,xl, &
1625 ids,ide, jds,jde, kds,kde, &
1626 ims,ime, jms,jme, kms,kme, &
1627 its,ite, jts,jte, kts,kte )
1628
1629 IMPLICIT NONE
1630
1631 integer &
1632 ,intent (in ) :: &
1633 ids,ide, jds,jde, kds,kde, &
1634 ims,ime, jms,jme, kms,kme, &
1635 its,ite, jts,jte, kts,kte
1636 ! cdd= detrainment function
1637 ! q = environmental q on model levels
1638 ! q_cup = environmental q on model cloud levels
1639 ! qes_cup = saturation q on model cloud levels
1640 ! hes_cup = saturation h on model cloud levels
1641 ! hcd = h in model cloud
1642 ! bu = buoancy term
1643 ! zd = normalized downdraft mass flux
1644 ! gamma_cup = gamma on model cloud levels
1645 ! mentr_rate = entrainment rate
1646 ! qcd = cloud q (including liquid water) after entrainment
1647 ! qrch = saturation q in cloud
1648 ! pwd = evaporate at that level
1649 ! pwev = total normalized integrated evaoprate (I2)
1650 ! entr= entrainment rate
1651 !
1652 real, dimension (its:ite,kts:kte) &
1653 ,intent (in ) :: &
1654 zd,t_cup,hes_cup,hcd,qes_cup,q_cup,z_cup,cdd,gamma_cup,q,he
1655 real &
1656 ,intent (in ) :: &
1657 entr,xl
1658 integer &
1659 ,intent (in ) :: &
1660 iloop
1661 integer, dimension (its:ite) &
1662 ,intent (in ) :: &
1663 jmin
1664 integer, dimension (its:ite) &
1665 ,intent (inout) :: &
1666 ierr
1667 real, dimension (its:ite,kts:kte) &
1668 ,intent (out ) :: &
1669 qcd,qrcd,pwd
1670 real, dimension (its:ite) &
1671 ,intent (out ) :: &
1672 pwev,bu
1673 !
1674 ! local variables in this routine
1675 !
1676
1677 integer :: &
1678 i,k,ki
1679 real :: &
1680 dh,dz,dqeva
1681
1682 integer :: itf, ktf
1683
1684 itf=MIN(ite,ide-1)
1685 ktf=MIN(kte,kde-1)
1686
1687 do i=its,itf
1688 bu(i)=0.
1689 pwev(i)=0.
1690 enddo
1691 do k=kts,ktf
1692 do i=its,itf
1693 qcd(i,k)=0.
1694 qrcd(i,k)=0.
1695 pwd(i,k)=0.
1696 enddo
1697 enddo
1698 !
1699 !
1700 !
1701 do 100 i=its,itf
1702 IF(ierr(I).eq.0)then
1703 k=jmin(i)
1704 DZ=Z_cup(i,K+1)-Z_cup(i,K)
1705 qcd(i,k)=q_cup(i,k)
1706 ! qcd(i,k)=.5*(qes_cup(i,k)+q_cup(i,k))
1707 qrcd(i,k)=qes_cup(i,k)
1708 pwd(i,jmin(i))=min(0.,qcd(i,k)-qrcd(i,k))
1709 pwev(i)=pwev(i)+pwd(i,jmin(i))
1710 qcd(i,k)=qes_cup(i,k)
1711 !
1712 DH=HCD(I,k)-HES_cup(I,K)
1713 bu(i)=dz*dh
1714 do ki=jmin(i)-1,1,-1
1715 DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
1716 QCD(i,Ki)=(qCD(i,Ki+1)*(1.-.5*CDD(i,Ki)*DZ) &
1717 +entr*DZ*q(i,Ki) &
1718 )/(1.+entr*DZ-.5*CDD(i,Ki)*DZ)
1719 !
1720 !--- to be negatively buoyant, hcd should be smaller than hes!
1721 !
1722 DH=HCD(I,ki)-HES_cup(I,Ki)
1723 bu(i)=bu(i)+dz*dh
1724 QRCD(I,Ki)=qes_cup(i,ki)+(1./XL)*(GAMMA_cup(i,ki) &
1725 /(1.+GAMMA_cup(i,ki)))*DH
1726 dqeva=qcd(i,ki)-qrcd(i,ki)
1727 if(dqeva.gt.0.)dqeva=0.
1728 pwd(i,ki)=zd(i,ki)*dqeva
1729 qcd(i,ki)=qrcd(i,ki)
1730 pwev(i)=pwev(i)+pwd(i,ki)
1731 ! if(iloop.eq.1.and.i.eq.102.and.j.eq.62)then
1732 ! print *,'in cup_dd_moi ', hcd(i,ki),HES_cup(I,Ki),dh,dqeva
1733 ! endif
1734 enddo
1735 !
1736 !--- end loop over i
1737 if(pwev(I).eq.0.and.iloop.eq.1)then
1738 ! print *,'problem with buoy in cup_dd_moisture',i
1739 ierr(i)=7
1740 endif
1741 if(BU(I).GE.0.and.iloop.eq.1)then
1742 ! print *,'problem with buoy in cup_dd_moisture',i
1743 ierr(i)=7
1744 endif
1745 endif
1746 100 continue
1747
1748 END SUBROUTINE cup_dd_moisture
1749
1750
1751 SUBROUTINE cup_dd_nms(zd,z_cup,cdd,entr,jmin,ierr, &
1752 itest,kdet,z1, &
1753 ids,ide, jds,jde, kds,kde, &
1754 ims,ime, jms,jme, kms,kme, &
1755 its,ite, jts,jte, kts,kte )
1756
1757 IMPLICIT NONE
1758 !
1759 ! on input
1760 !
1761
1762 ! only local wrf dimensions are need as of now in this routine
1763
1764 integer &
1765 ,intent (in ) :: &
1766 ids,ide, jds,jde, kds,kde, &
1767 ims,ime, jms,jme, kms,kme, &
1768 its,ite, jts,jte, kts,kte
1769 ! z_cup = height of cloud model level
1770 ! z1 = terrain elevation
1771 ! entr = downdraft entrainment rate
1772 ! jmin = downdraft originating level
1773 ! kdet = level above ground where downdraft start detraining
1774 ! itest = flag to whether to calculate cdd
1775
1776 real, dimension (its:ite,kts:kte) &
1777 ,intent (in ) :: &
1778 z_cup
1779 real, dimension (its:ite) &
1780 ,intent (in ) :: &
1781 z1
1782 real &
1783 ,intent (in ) :: &
1784 entr
1785 integer, dimension (its:ite) &
1786 ,intent (in ) :: &
1787 jmin,kdet
1788 integer &
1789 ,intent (in ) :: &
1790 itest
1791 !
1792 ! input and output
1793 !
1794
1795 ! ierr error value, maybe modified in this routine
1796
1797 integer, dimension (its:ite) &
1798 ,intent (inout) :: &
1799 ierr
1800 ! zd is the normalized downdraft mass flux
1801 ! cdd is the downdraft detrainmen function
1802
1803 real, dimension (its:ite,kts:kte) &
1804 ,intent (out ) :: &
1805 zd
1806 real, dimension (its:ite,kts:kte) &
1807 ,intent (inout) :: &
1808 cdd
1809 !
1810 ! local variables in this routine
1811 !
1812
1813 integer :: &
1814 i,k,ki
1815 real :: &
1816 a,perc,dz
1817
1818 integer :: itf, ktf
1819
1820 itf=MIN(ite,ide-1)
1821 ktf=MIN(kte,kde-1)
1822 !
1823 !--- perc is the percentage of mass left when hitting the ground
1824 !
1825 perc=.03
1826
1827 do k=kts,ktf
1828 do i=its,itf
1829 zd(i,k)=0.
1830 if(itest.eq.0)cdd(i,k)=0.
1831 enddo
1832 enddo
1833 a=1.-perc
1834 !
1835 !
1836 !
1837 do 100 i=its,itf
1838 IF(ierr(I).eq.0)then
1839 zd(i,jmin(i))=1.
1840 !
1841 !--- integrate downward, specify detrainment(cdd)!
1842 !
1843 do ki=jmin(i)-1,1,-1
1844 DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
1845 if(ki.le.kdet(i).and.itest.eq.0)then
1846 cdd(i,ki)=entr+(1.- (a*(z_cup(i,ki)-z1(i)) &
1847 +perc*(z_cup(i,kdet(i))-z1(i)) ) &
1848 /(a*(z_cup(i,ki+1)-z1(i)) &
1849 +perc*(z_cup(i,kdet(i))-z1(i))))/dz
1850 endif
1851 zd(i,ki)=zd(i,ki+1)*(1.+(entr-cdd(i,ki))*dz)
1852 enddo
1853 !
1854 endif
1855 !--- end loop over i
1856 100 continue
1857
1858 END SUBROUTINE cup_dd_nms
1859
1860
1861 SUBROUTINE cup_dellabot(ipr,jpr,he_cup,ierr,z_cup,p_cup, &
1862 hcd,edt,zd,cdd,he,della,j,mentrd_rate,z,g, &
1863 ids,ide, jds,jde, kds,kde, &
1864 ims,ime, jms,jme, kms,kme, &
1865 its,ite, jts,jte, kts,kte )
1866
1867 IMPLICIT NONE
1868
1869 integer &
1870 ,intent (in ) :: &
1871 ids,ide, jds,jde, kds,kde, &
1872 ims,ime, jms,jme, kms,kme, &
1873 its,ite, jts,jte, kts,kte
1874 integer, intent (in ) :: &
1875 j,ipr,jpr
1876 !
1877 ! ierr error value, maybe modified in this routine
1878 !
1879 real, dimension (its:ite,kts:kte) &
1880 ,intent (out ) :: &
1881 della
1882 real, dimension (its:ite,kts:kte) &
1883 ,intent (in ) :: &
1884 z_cup,p_cup,hcd,zd,cdd,he,z,he_cup
1885 real, dimension (its:ite) &
1886 ,intent (in ) :: &
1887 edt
1888 real &
1889 ,intent (in ) :: &
1890 g,mentrd_rate
1891 integer, dimension (its:ite) &
1892 ,intent (inout) :: &
1893 ierr
1894 !
1895 ! local variables in this routine
1896 !
1897
1898 integer i
1899 real detdo,detdo1,detdo2,entdo,dp,dz,subin, &
1900 totmas
1901 !
1902 integer :: itf, ktf
1903
1904 itf=MIN(ite,ide-1)
1905 ktf=MIN(kte,kde-1)
1906 !
1907 !
1908 ! if(j.eq.jpr)print *,'in cup dellabot '
1909 do 100 i=its,itf
1910 della(i,1)=0.
1911 if(ierr(i).ne.0)go to 100
1912 dz=z_cup(i,2)-z_cup(i,1)
1913 DP=100.*(p_cup(i,1)-P_cup(i,2))
1914 detdo1=edt(i)*zd(i,2)*CDD(i,1)*DZ
1915 detdo2=edt(i)*zd(i,1)
1916 entdo=edt(i)*zd(i,2)*mentrd_rate*dz
1917 subin=-EDT(I)*zd(i,2)
1918 detdo=detdo1+detdo2-entdo+subin
1919 DELLA(I,1)=(detdo1*.5*(HCD(i,1)+HCD(i,2)) &
1920 +detdo2*hcd(i,1) &
1921 +subin*he_cup(i,2) &
1922 -entdo*he(i,1))*g/dp
1923 100 CONTINUE
1924
1925 END SUBROUTINE cup_dellabot
1926
1927
1928 SUBROUTINE cup_dellas(ierr,z_cup,p_cup,hcd,edt,zd,cdd, &
1929 he,della,j,mentrd_rate,zu,g, &
1930 cd,hc,ktop,k22,kbcon,mentr_rate,jmin,he_cup,kdet,kpbl, &
1931 ipr,jpr,name, &
1932 ids,ide, jds,jde, kds,kde, &
1933 ims,ime, jms,jme, kms,kme, &
1934 its,ite, jts,jte, kts,kte )
1935
1936 IMPLICIT NONE
1937
1938 integer &
1939 ,intent (in ) :: &
1940 ids,ide, jds,jde, kds,kde, &
1941 ims,ime, jms,jme, kms,kme, &
1942 its,ite, jts,jte, kts,kte
1943 integer, intent (in ) :: &
1944 j,ipr,jpr
1945 !
1946 ! ierr error value, maybe modified in this routine
1947 !
1948 real, dimension (its:ite,kts:kte) &
1949 ,intent (out ) :: &
1950 della
1951 real, dimension (its:ite,kts:kte) &
1952 ,intent (in ) :: &
1953 z_cup,p_cup,hcd,zd,cdd,he,hc,cd,zu,he_cup
1954 real, dimension (its:ite) &
1955 ,intent (in ) :: &
1956 edt
1957 real &
1958 ,intent (in ) :: &
1959 g,mentrd_rate,mentr_rate
1960 integer, dimension (its:ite) &
1961 ,intent (in ) :: &
1962 kbcon,ktop,k22,jmin,kdet,kpbl
1963 integer, dimension (its:ite) &
1964 ,intent (inout) :: &
1965 ierr
1966 character *(*), intent (in) :: &
1967 name
1968 !
1969 ! local variables in this routine
1970 !
1971
1972 integer i,k
1973 real detdo1,detdo2,entdo,dp,dz,subin,detdo,entup, &
1974 detup,subdown,entdoj,entupk,detupk,totmas
1975 !
1976 integer :: itf, ktf
1977
1978 itf=MIN(ite,ide-1)
1979 ktf=MIN(kte,kde-1)
1980 !
1981 !
1982 i=ipr
1983 ! if(j.eq.jpr)then
1984 ! print *,'in dellas kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)'
1985 ! print *,kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)
1986 ! endif
1987 DO K=kts+1,ktf
1988 do i=its,itf
1989 della(i,k)=0.
1990 enddo
1991 enddo
1992 !
1993 DO 100 k=kts+1,ktf-1
1994 DO 100 i=its,ite
1995 IF(ierr(i).ne.0)GO TO 100
1996 IF(K.Gt.KTOP(I))GO TO 100
1997 !
1998 !--- SPECIFY DETRAINMENT OF DOWNDRAFT, HAS TO BE CONSISTENT
1999 !--- WITH ZD CALCULATIONS IN SOUNDD.
2000 !
2001 DZ=Z_cup(I,K+1)-Z_cup(I,K)
2002 detdo=edt(i)*CDD(i,K)*DZ*ZD(i,k+1)
2003 entdo=edt(i)*mentrd_rate*dz*zd(i,k+1)
2004 subin=zu(i,k+1)-zd(i,k+1)*edt(i)
2005 entup=0.
2006 detup=0.
2007 if(k.ge.kbcon(i).and.k.lt.ktop(i))then
2008 entup=mentr_rate*dz*zu(i,k)
2009 detup=CD(i,K+1)*DZ*ZU(i,k)
2010 endif
2011 subdown=(zu(i,k)-zd(i,k)*edt(i))
2012 entdoj=0.
2013 entupk=0.
2014 detupk=0.
2015 !
2016 if(k.eq.jmin(i))then
2017 entdoj=edt(i)*zd(i,k)
2018 endif
2019
2020 if(k.eq.k22(i)-1)then
2021 entupk=zu(i,kpbl(i))
2022 endif
2023
2024 if(k.gt.kdet(i))then
2025 detdo=0.
2026 endif
2027
2028 if(k.eq.ktop(i)-0)then
2029 detupk=zu(i,ktop(i))
2030 subin=0.
2031 endif
2032 if(k.lt.kbcon(i))then
2033 detup=0.
2034 endif
2035 !C
2036 !C--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
2037 !C
2038 totmas=subin-subdown+detup-entup-entdo+ &
2039 detdo-entupk-entdoj+detupk
2040 ! if(j.eq.jpr.and.i.eq.ipr)print *,'k,totmas,sui,sud = ',k,
2041 ! 1 totmas,subin,subdown
2042 ! if(j.eq.jpr.and.i.eq.ipr)print *,'updr stuff = ',detup,
2043 ! 1 entup,entupk,detupk
2044 ! if(j.eq.jpr.and.i.eq.ipr)print *,'dddr stuff = ',entdo,
2045 ! 1 detdo,entdoj
2046 if(abs(totmas).gt.1.e-6)then
2047 ! print *,'*********************',i,j,k,totmas,name
2048 ! print *,kpbl(i),k22(i),kbcon(i),ktop(i)
2049 !c print *,'updr stuff = ',subin,
2050 !c 1 subdown,detup,entup,entupk,detupk
2051 !c print *,'dddr stuff = ',entdo,
2052 !c 1 detdo,entdoj
2053 ! call wrf_error_fatal ( 'totmas .gt.1.e-6' )
2054 endif
2055 dp=100.*(p_cup(i,k-1)-p_cup(i,k))
2056 della(i,k)=(subin*he_cup(i,k+1) &
2057 -subdown*he_cup(i,k) &
2058 +detup*.5*(HC(i,K+1)+HC(i,K)) &
2059 +detdo*.5*(HCD(i,K+1)+HCD(i,K)) &
2060 -entup*he(i,k) &
2061 -entdo*he(i,k) &
2062 -entupk*he_cup(i,k22(i)) &
2063 -entdoj*he_cup(i,jmin(i)) &
2064 +detupk*hc(i,ktop(i)) &
2065 )*g/dp
2066 ! if(i.eq.ipr.and.j.eq.jpr)then
2067 ! print *,k,della(i,k),subin*he_cup(i,k+1),subdown*he_cup(i,k),
2068 ! 1 detdo*.5*(HCD(i,K+1)+HCD(i,K))
2069 ! print *,k,detup*.5*(HC(i,K+1)+HC(i,K)),detupk*hc(i,ktop(i)),
2070 ! 1 entup*he(i,k),entdo*he(i,k)
2071 ! print *,k,he_cup(i,k+1),he_cup(i,k),entupk*he_cup(i,k)
2072 ! endif
2073
2074 100 CONTINUE
2075
2076 END SUBROUTINE cup_dellas
2077
2078
2079 SUBROUTINE cup_direction2(i,j,dir,id,massflx, &
2080 iresult,imass,massfld, &
2081 ids,ide, jds,jde, kds,kde, &
2082 ims,ime, jms,jme, kms,kme, &
2083 its,ite, jts,jte, kts,kte )
2084
2085 IMPLICIT NONE
2086
2087 integer &
2088 ,intent (in ) :: &
2089 ids,ide, jds,jde, kds,kde, &
2090 ims,ime, jms,jme, kms,kme, &
2091 its,ite, jts,jte, kts,kte
2092 integer, intent (in ) :: &
2093 i,j,imass
2094 integer, intent (out ) :: &
2095 iresult
2096 !
2097 ! ierr error value, maybe modified in this routine
2098 !
2099 integer, dimension (ims:ime,jms:jme) &
2100 ,intent (in ) :: &
2101 id
2102 real, dimension (ims:ime,jms:jme) &
2103 ,intent (in ) :: &
2104 massflx
2105 real, dimension (its:ite) &
2106 ,intent (inout) :: &
2107 dir
2108 real &
2109 ,intent (out ) :: &
2110 massfld
2111 !
2112 ! local variables in this routine
2113 !
2114
2115 integer k,ia,ja,ib,jb
2116 real diff
2117 !
2118 !
2119 !
2120 if(imass.eq.1)then
2121 massfld=massflx(i,j)
2122 endif
2123 iresult=0
2124 ! return
2125 diff=22.5
2126 if(dir(i).lt.22.5)dir(i)=360.+dir(i)
2127 if(id(i,j).eq.1)iresult=1
2128 ! ja=max(2,j-1)
2129 ! ia=max(2,i-1)
2130 ! jb=min(mjx-1,j+1)
2131 ! ib=min(mix-1,i+1)
2132 ja=j-1
2133 ia=i-1
2134 jb=j+1
2135 ib=i+1
2136 if(dir(i).gt.90.-diff.and.dir(i).le.90.+diff)then
2137 !--- steering flow from the east
2138 if(id(ib,j).eq.1)then
2139 iresult=1
2140 if(imass.eq.1)then
2141 massfld=max(massflx(ib,j),massflx(i,j))
2142 endif
2143 return
2144 endif
2145 else if(dir(i).gt.135.-diff.and.dir(i).le.135.+diff)then
2146 !--- steering flow from the south-east
2147 if(id(ib,ja).eq.1)then
2148 iresult=1
2149 if(imass.eq.1)then
2150 massfld=max(massflx(ib,ja),massflx(i,j))
2151 endif
2152 return
2153 endif
2154 !--- steering flow from the south
2155 else if(dir(i).gt.180.-diff.and.dir(i).le.180.+diff)then
2156 if(id(i,ja).eq.1)then
2157 iresult=1
2158 if(imass.eq.1)then
2159 massfld=max(massflx(i,ja),massflx(i,j))
2160 endif
2161 return
2162 endif
2163 !--- steering flow from the south west
2164 else if(dir(i).gt.225.-diff.and.dir(i).le.225.+diff)then
2165 if(id(ia,ja).eq.1)then
2166 iresult=1
2167 if(imass.eq.1)then
2168 massfld=max(massflx(ia,ja),massflx(i,j))
2169 endif
2170 return
2171 endif
2172 !--- steering flow from the west
2173 else if(dir(i).gt.270.-diff.and.dir(i).le.270.+diff)then
2174 if(id(ia,j).eq.1)then
2175 iresult=1
2176 if(imass.eq.1)then
2177 massfld=max(massflx(ia,j),massflx(i,j))
2178 endif
2179 return
2180 endif
2181 !--- steering flow from the north-west
2182 else if(dir(i).gt.305.-diff.and.dir(i).le.305.+diff)then
2183 if(id(ia,jb).eq.1)then
2184 iresult=1
2185 if(imass.eq.1)then
2186 massfld=max(massflx(ia,jb),massflx(i,j))
2187 endif
2188 return
2189 endif
2190 !--- steering flow from the north
2191 else if(dir(i).gt.360.-diff.and.dir(i).le.360.+diff)then
2192 if(id(i,jb).eq.1)then
2193 iresult=1
2194 if(imass.eq.1)then
2195 massfld=max(massflx(i,jb),massflx(i,j))
2196 endif
2197 return
2198 endif
2199 !--- steering flow from the north-east
2200 else if(dir(i).gt.45.-diff.and.dir(i).le.45.+diff)then
2201 if(id(ib,jb).eq.1)then
2202 iresult=1
2203 if(imass.eq.1)then
2204 massfld=max(massflx(ib,jb),massflx(i,j))
2205 endif
2206 return
2207 endif
2208 endif
2209
2210 END SUBROUTINE cup_direction2
2211
2212
2213 SUBROUTINE cup_env(z,qes,he,hes,t,q,p,z1, &
2214 psur,ierr,tcrit,itest,xl,cp, &
2215 ids,ide, jds,jde, kds,kde, &
2216 ims,ime, jms,jme, kms,kme, &
2217 its,ite, jts,jte, kts,kte )
2218
2219 IMPLICIT NONE
2220
2221 integer &
2222 ,intent (in ) :: &
2223 ids,ide, jds,jde, kds,kde, &
2224 ims,ime, jms,jme, kms,kme, &
2225 its,ite, jts,jte, kts,kte
2226 !
2227 ! ierr error value, maybe modified in this routine
2228 ! q = environmental mixing ratio
2229 ! qes = environmental saturation mixing ratio
2230 ! t = environmental temp
2231 ! tv = environmental virtual temp
2232 ! p = environmental pressure
2233 ! z = environmental heights
2234 ! he = environmental moist static energy
2235 ! hes = environmental saturation moist static energy
2236 ! psur = surface pressure
2237 ! z1 = terrain elevation
2238 !
2239 !
2240 real, dimension (its:ite,kts:kte) &
2241 ,intent (in ) :: &
2242 p,t
2243 real, dimension (its:ite,kts:kte) &
2244 ,intent (out ) :: &
2245 he,hes,qes
2246 real, dimension (its:ite,kts:kte) &
2247 ,intent (inout) :: &
2248 z,q
2249 real, dimension (its:ite) &
2250 ,intent (in ) :: &
2251 psur,z1
2252 real &
2253 ,intent (in ) :: &
2254 xl,cp
2255 integer, dimension (its:ite) &
2256 ,intent (inout) :: &
2257 ierr
2258 integer &
2259 ,intent (in ) :: &
2260 itest
2261 !
2262 ! local variables in this routine
2263 !
2264
2265 integer :: &
2266 i,k,iph
2267 real, dimension (1:2) :: AE,BE,HT
2268 real, dimension (its:ite,kts:kte) :: tv
2269 real :: tcrit,e,tvbar
2270
2271 integer :: itf, ktf
2272
2273 itf=MIN(ite,ide-1)
2274 ktf=MIN(kte,kde-1)
2275
2276 HT(1)=XL/CP
2277 HT(2)=2.834E6/CP
2278 BE(1)=.622*HT(1)/.286
2279 AE(1)=BE(1)/273.+ALOG(610.71)
2280 BE(2)=.622*HT(2)/.286
2281 AE(2)=BE(2)/273.+ALOG(610.71)
2282 ! print *, 'TCRIT = ', tcrit,its,ite
2283 DO k=kts,ktf
2284 do i=its,itf
2285 if(ierr(i).eq.0)then
2286 !Csgb - IPH is for phase, dependent on TCRIT (water or ice)
2287 IPH=1
2288 IF(T(I,K).LE.TCRIT)IPH=2
2289 ! print *, 'AE(IPH),BE(IPH) = ',AE(IPH),BE(IPH),AE(IPH)-BE(IPH),T(i,k),i,k
2290 E=EXP(AE(IPH)-BE(IPH)/T(I,K))
2291 ! print *, 'P, E = ', P(I,K), E
2292 QES(I,K)=.622*E/(100.*P(I,K)-E)
2293 IF(QES(I,K).LE.1.E-08)QES(I,K)=1.E-08
2294 IF(Q(I,K).GT.QES(I,K))Q(I,K)=QES(I,K)
2295 TV(I,K)=T(I,K)+.608*Q(I,K)*T(I,K)
2296 endif
2297 enddo
2298 enddo
2299 !
2300 !--- z's are calculated with changed h's and q's and t's
2301 !--- if itest=2
2302 !
2303 if(itest.ne.2)then
2304 do i=its,itf
2305 if(ierr(i).eq.0)then
2306 Z(I,1)=max(0.,Z1(I))-(ALOG(P(I,1))- &
2307 ALOG(PSUR(I)))*287.*TV(I,1)/9.81
2308 endif
2309 enddo
2310
2311 ! --- calculate heights
2312 DO K=kts+1,ktf
2313 do i=its,itf
2314 if(ierr(i).eq.0)then
2315 TVBAR=.5*TV(I,K)+.5*TV(I,K-1)
2316 Z(I,K)=Z(I,K-1)-(ALOG(P(I,K))- &
2317 ALOG(P(I,K-1)))*287.*TVBAR/9.81
2318 endif
2319 enddo
2320 enddo
2321 else
2322 do k=kts,ktf
2323 do i=its,itf
2324 if(ierr(i).eq.0)then
2325 z(i,k)=(he(i,k)-1004.*t(i,k)-2.5e6*q(i,k))/9.81
2326 z(i,k)=max(1.e-3,z(i,k))
2327 endif
2328 enddo
2329 enddo
2330 endif
2331 !
2332 !--- calculate moist static energy - HE
2333 ! saturated moist static energy - HES
2334 !
2335 DO k=kts,ktf
2336 do i=its,itf
2337 if(ierr(i).eq.0)then
2338 if(itest.eq.0)HE(I,K)=9.81*Z(I,K)+1004.*T(I,K)+2.5E06*Q(I,K)
2339 HES(I,K)=9.81*Z(I,K)+1004.*T(I,K)+2.5E06*QES(I,K)
2340 IF(HE(I,K).GE.HES(I,K))HE(I,K)=HES(I,K)
2341 ! if(i.eq.2)then
2342 ! print *,k,z(i,k),t(i,k),p(i,k),he(i,k),hes(i,k)
2343 ! endif
2344 endif
2345 enddo
2346 enddo
2347
2348 END SUBROUTINE cup_env
2349
2350
2351 SUBROUTINE cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup, &
2352 he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
2353 ierr,z1,xl,rv,cp, &
2354 ids,ide, jds,jde, kds,kde, &
2355 ims,ime, jms,jme, kms,kme, &
2356 its,ite, jts,jte, kts,kte )
2357
2358 IMPLICIT NONE
2359
2360 integer &
2361 ,intent (in ) :: &
2362 ids,ide, jds,jde, kds,kde, &
2363 ims,ime, jms,jme, kms,kme, &
2364 its,ite, jts,jte, kts,kte
2365 !
2366 ! ierr error value, maybe modified in this routine
2367 ! q = environmental mixing ratio
2368 ! q_cup = environmental mixing ratio on cloud levels
2369 ! qes = environmental saturation mixing ratio
2370 ! qes_cup = environmental saturation mixing ratio on cloud levels
2371 ! t = environmental temp
2372 ! t_cup = environmental temp on cloud levels
2373 ! p = environmental pressure
2374 ! p_cup = environmental pressure on cloud levels
2375 ! z = environmental heights
2376 ! z_cup = environmental heights on cloud levels
2377 ! he = environmental moist static energy
2378 ! he_cup = environmental moist static energy on cloud levels
2379 ! hes = environmental saturation moist static energy
2380 ! hes_cup = environmental saturation moist static energy on cloud levels
2381 ! gamma_cup = gamma on cloud levels
2382 ! psur = surface pressure
2383 ! z1 = terrain elevation
2384 !
2385 !
2386 real, dimension (its:ite,kts:kte) &
2387 ,intent (in ) :: &
2388 qes,q,he,hes,z,p,t
2389 real, dimension (its:ite,kts:kte) &
2390 ,intent (out ) :: &
2391 qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup
2392 real, dimension (its:ite) &
2393 ,intent (in ) :: &
2394 psur,z1
2395 real &
2396 ,intent (in ) :: &
2397 xl,rv,cp
2398 integer, dimension (its:ite) &
2399 ,intent (inout) :: &
2400 ierr
2401 !
2402 ! local variables in this routine
2403 !
2404
2405 integer :: &
2406 i,k
2407
2408 integer :: itf, ktf
2409
2410 itf=MIN(ite,ide-1)
2411 ktf=MIN(kte,kde-1)
2412
2413 do k=kts+1,ktf
2414 do i=its,itf
2415 if(ierr(i).eq.0)then
2416 qes_cup(i,k)=.5*(qes(i,k-1)+qes(i,k))
2417 q_cup(i,k)=.5*(q(i,k-1)+q(i,k))
2418 hes_cup(i,k)=.5*(hes(i,k-1)+hes(i,k))
2419 he_cup(i,k)=.5*(he(i,k-1)+he(i,k))
2420 if(he_cup(i,k).gt.hes_cup(i,k))he_cup(i,k)=hes_cup(i,k)
2421 z_cup(i,k)=.5*(z(i,k-1)+z(i,k))
2422 p_cup(i,k)=.5*(p(i,k-1)+p(i,k))
2423 t_cup(i,k)=.5*(t(i,k-1)+t(i,k))
2424 gamma_cup(i,k)=(xl/cp)*(xl/(rv*t_cup(i,k) &
2425 *t_cup(i,k)))*qes_cup(i,k)
2426 endif
2427 enddo
2428 enddo
2429 do i=its,itf
2430 if(ierr(i).eq.0)then
2431 qes_cup(i,1)=qes(i,1)
2432 q_cup(i,1)=q(i,1)
2433 hes_cup(i,1)=hes(i,1)
2434 he_cup(i,1)=he(i,1)
2435 z_cup(i,1)=.5*(z(i,1)+z1(i))
2436 p_cup(i,1)=.5*(p(i,1)+psur(i))
2437 t_cup(i,1)=t(i,1)
2438 gamma_cup(i,1)=xl/cp*(xl/(rv*t_cup(i,1) &
2439 *t_cup(i,1)))*qes_cup(i,1)
2440 endif
2441 enddo
2442
2443 END SUBROUTINE cup_env_clev
2444
2445
2446 SUBROUTINE cup_forcing_ens(closure_n,xland,aa0,aa1,xaa0,mbdt,dtime,ierr,ierr2,ierr3,&
2447 xf_ens,j,name,maxens,iens,iedt,maxens2,maxens3,mconv, &
2448 p_cup,ktop,omeg,zd,k22,zu,pr_ens,edt,kbcon,massflx, &
2449 iact_old_gr,dir,ensdim,massfln,icoic, &
2450 ids,ide, jds,jde, kds,kde, &
2451 ims,ime, jms,jme, kms,kme, &
2452 its,ite, jts,jte, kts,kte )
2453
2454 IMPLICIT NONE
2455
2456 integer &
2457 ,intent (in ) :: &
2458 ids,ide, jds,jde, kds,kde, &
2459 ims,ime, jms,jme, kms,kme, &
2460 its,ite, jts,jte, kts,kte
2461 integer, intent (in ) :: &
2462 j,ensdim,maxens,iens,iedt,maxens2,maxens3
2463 !
2464 ! ierr error value, maybe modified in this routine
2465 ! pr_ens = precipitation ensemble
2466 ! xf_ens = mass flux ensembles
2467 ! massfln = downdraft mass flux ensembles used in next timestep
2468 ! omeg = omega from large scale model
2469 ! mconv = moisture convergence from large scale model
2470 ! zd = downdraft normalized mass flux
2471 ! zu = updraft normalized mass flux
2472 ! aa0 = cloud work function without forcing effects
2473 ! aa1 = cloud work function with forcing effects
2474 ! xaa0 = cloud work function with cloud effects (ensemble dependent)
2475 ! edt = epsilon
2476 ! dir = "storm motion"
2477 ! mbdt = arbitrary numerical parameter
2478 ! dtime = dt over which forcing is applied
2479 ! iact_gr_old = flag to tell where convection was active
2480 ! kbcon = LFC of parcel from k22
2481 ! k22 = updraft originating level
2482 ! icoic = flag if only want one closure (usually set to zero!)
2483 ! name = deep or shallow convection flag
2484 !
2485 real, dimension (ims:ime,jms:jme,1:ensdim) &
2486 ,intent (inout) :: &
2487 pr_ens
2488 real, dimension (ims:ime,jms:jme,1:ensdim) &
2489 ,intent (out ) :: &
2490 xf_ens,massfln
2491 real, dimension (ims:ime,jms:jme) &
2492 ,intent (in ) :: &
2493 massflx
2494 real, dimension (its:ite,kts:kte) &
2495 ,intent (in ) :: &
2496 omeg,zd,zu,p_cup
2497 real, dimension (its:ite,1:maxens) &
2498 ,intent (in ) :: &
2499 xaa0
2500 real, dimension (its:ite) &
2501 ,intent (in ) :: &
2502 aa1,edt,dir,mconv,xland
2503 real, dimension (its:ite) &
2504 ,intent (inout) :: &
2505 aa0,closure_n
2506 real, dimension (1:maxens) &
2507 ,intent (in ) :: &
2508 mbdt
2509 real &
2510 ,intent (in ) :: &
2511 dtime
2512 integer, dimension (ims:ime,jms:jme) &
2513 ,intent (in ) :: &
2514 iact_old_gr
2515 integer, dimension (its:ite) &
2516 ,intent (in ) :: &
2517 k22,kbcon,ktop
2518 integer, dimension (its:ite) &
2519 ,intent (inout) :: &
2520 ierr,ierr2,ierr3
2521 integer &
2522 ,intent (in ) :: &
2523 icoic
2524 character *(*), intent (in) :: &
2525 name
2526 !
2527 ! local variables in this routine
2528 !
2529
2530 real, dimension (1:maxens3) :: &
2531 xff_ens3
2532 real, dimension (1:maxens) :: &
2533 xk
2534 integer :: &
2535 i,k,nall,n,ne,nens,nens3,iresult,iresultd,iresulte,mkxcrt,kclim
2536 parameter (mkxcrt=15)
2537 real :: &
2538 a1,massfld,xff0,xomg,aclim1,aclim2,aclim3,aclim4
2539 real, dimension(1:mkxcrt) :: &
2540 pcrit,acrit,acritt
2541
2542 integer :: itf,nall2
2543
2544 itf=MIN(ite,ide-1)
2545
2546 DATA PCRIT/850.,800.,750.,700.,650.,600.,550.,500.,450.,400., &
2547 350.,300.,250.,200.,150./
2548 DATA ACRIT/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
2549 .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
2550 ! GDAS DERIVED ACRIT
2551 DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688, &
2552 .743,.813,.886,.947,1.138,1.377,1.896/
2553 !
2554 nens=0
2555
2556 !--- LARGE SCALE FORCING
2557 !
2558 DO 100 i=its,itf
2559 ! if(i.eq.ipr.and.j.eq.jpr)print *,'ierr = ',ierr(i)
2560 if(name.eq.'deeps'.and.ierr(i).gt.995)then
2561 ! print *,i,j,ierr(i),aa0(i)
2562 aa0(i)=0.
2563 ierr(i)=0
2564 endif
2565 IF(ierr(i).eq.0)then
2566 ! kclim=0
2567 do k=mkxcrt,1,-1
2568 if(p_cup(i,ktop(i)).lt.pcrit(k))then
2569 kclim=k
2570 go to 9
2571 endif
2572 enddo
2573 if(p_cup(i,ktop(i)).ge.pcrit(1))kclim=1
2574 9 continue
2575 kclim=max(kclim,1)
2576 k=max(kclim-1,1)
2577 aclim1=acrit(kclim)*1.e3
2578 aclim2=acrit(k)*1.e3
2579 aclim3=acritt(kclim)*1.e3
2580 aclim4=acritt(k)*1.e3
2581 ! print *,'p_cup(ktop(i)),kclim,pcrit(kclim)'
2582 ! print *,p_cup(i,ktop(i)),kclim,pcrit(kclim)
2583 ! print *,'aclim1,aclim2,aclim3,aclim4'
2584 ! print *,aclim1,aclim2,aclim3,aclim4
2585 ! print *,dtime,name,ierr(i),aa1(i),aa0(i)
2586 ! print *,dtime,name,ierr(i),aa1(i),aa0(i)
2587 !
2588 !--- treatment different for this closure
2589 !
2590 if(name.eq.'deeps')then
2591 !
2592 xff0= (AA1(I)-AA0(I))/DTIME
2593 xff_ens3(1)=(AA1(I)-AA0(I))/dtime
2594 xff_ens3(2)=.9*xff_ens3(1)
2595 xff_ens3(3)=1.1*xff_ens3(1)
2596 !
2597 !--- more original Arakawa-Schubert (climatologic value of aa0)
2598 !
2599 !
2600 !--- omeg is in bar/s, mconv done with omeg in Pa/s
2601 ! more like Brown (1979), or Frank-Cohen (199?)
2602 !
2603 xff_ens3(4)=-omeg(i,k22(i))/9.81
2604 xff_ens3(5)=-omeg(i,kbcon(i))/9.81
2605 xff_ens3(6)=-omeg(i,1)/9.81
2606 do k=2,kbcon(i)-1
2607 xomg=-omeg(i,k)/9.81
2608 if(xomg.gt.xff_ens3(6))xff_ens3(6)=xomg
2609 enddo
2610 !
2611 !--- more like Krishnamurti et al.
2612 !
2613 xff_ens3(7)=mconv(i)
2614 xff_ens3(8)=mconv(i)
2615 xff_ens3(9)=mconv(i)
2616 !
2617 !--- more like Fritsch Chappel or Kain Fritsch (plus triggers)
2618 !
2619 xff_ens3(10)=AA1(I)/(60.*20.)
2620 xff_ens3(11)=AA1(I)/(60.*30.)
2621 xff_ens3(12)=AA1(I)/(60.*40.)
2622 !
2623 !--- more original Arakawa-Schubert (climatologic value of aa0)
2624 !
2625 xff_ens3(13)=max(0.,(AA1(I)-aclim1)/dtime)
2626 xff_ens3(14)=max(0.,(AA1(I)-aclim2)/dtime)
2627 xff_ens3(15)=max(0.,(AA1(I)-aclim3)/dtime)
2628 xff_ens3(16)=max(0.,(AA1(I)-aclim4)/dtime)
2629 ! if(ierr2(i).gt.0.or.ierr3(i).gt.0)then
2630 ! xff_ens3(10)=0.
2631 ! xff_ens3(11)=0.
2632 ! xff_ens3(12)=0.
2633 ! xff_ens3(13)=0.
2634 ! xff_ens3(14)=0.
2635 ! xff_ens3(15)=0.
2636 ! xff_ens3(16)=0.
2637 ! endif
2638
2639 do nens=1,maxens
2640 XK(nens)=(XAA0(I,nens)-AA1(I))/MBDT(2)
2641 if(xk(nens).le.0.and.xk(nens).gt.-1.e-6) &
2642 xk(nens)=-1.e-6
2643 if(xk(nens).gt.0.and.xk(nens).lt.1.e-6) &
2644 xk(nens)=1.e-6
2645 enddo
2646 !
2647 !--- add up all ensembles
2648 !
2649 do 350 ne=1,maxens
2650 !
2651 !--- for every xk, we have maxens3 xffs
2652 !--- iens is from outermost ensemble (most expensive!
2653 !
2654 !--- iedt (maxens2 belongs to it)
2655 !--- is from second, next outermost, not so expensive
2656 !
2657 !--- so, for every outermost loop, we have maxens*maxens2*3
2658 !--- ensembles!!! nall would be 0, if everything is on first
2659 !--- loop index, then ne would start counting, then iedt, then iens....
2660 !
2661 iresult=0
2662 iresultd=0
2663 iresulte=0
2664 nall=(iens-1)*maxens3*maxens*maxens2 &
2665 +(iedt-1)*maxens*maxens3 &
2666 +(ne-1)*maxens3
2667 !
2668 ! over water, enfor!e small cap for some of the closures
2669 !
2670 if(xland(i).lt.0.1)then
2671 if(ierr2(i).gt.0.or.ierr3(i).gt.0)then
2672 ! - ierr2 - 75 mb cap thickness, ierr3 - 125 cap thickness
2673
2674 ! - for larger cap, set Grell closure to zero
2675 xff_ens3(1) =0.
2676 massfln(i,j,nall+1)=0.
2677 xff_ens3(2) =0.
2678 massfln(i,j,nall+2)=0.
2679 xff_ens3(3) =0.
2680 massfln(i,j,nall+3)=0.
2681 closure_n(i)=closure_n(i)-1.
2682
2683 xff_ens3(7) =0.
2684 massfln(i,j,nall+7)=0.
2685 xff_ens3(8) =0.
2686 massfln(i,j,nall+8)=0.
2687 xff_ens3(9) =0.
2688 ! massfln(i,j,nall+9)=0.
2689 closure_n(i)=closure_n(i)-1.
2690 endif
2691 !
2692 ! also take out some closures in general
2693 !
2694 xff_ens3(4) =0.
2695 massfln(i,j,nall+4)=0.
2696 xff_ens3(5) =0.
2697 massfln(i,j,nall+5)=0.
2698 xff_ens3(6) =0.
2699 massfln(i,j,nall+6)=0.
2700 closure_n(i)=closure_n(i)-3.
2701
2702 xff_ens3(10)=0.
2703 massfln(i,j,nall+10)=0.
2704 xff_ens3(11)=0.
2705 massfln(i,j,nall+11)=0.
2706 xff_ens3(12)=0.
2707 massfln(i,j,nall+12)=0.
2708 if(ne.eq.1)closure_n(i)=closure_n(i)-3
2709 xff_ens3(13)=0.
2710 massfln(i,j,nall+13)=0.
2711 xff_ens3(14)=0.
2712 massfln(i,j,nall+14)=0.
2713 xff_ens3(15)=0.
2714 massfln(i,j,nall+15)=0.
2715 massfln(i,j,nall+16)=0.
2716 if(ne.eq.1)closure_n(i)=closure_n(i)-4
2717
2718 endif
2719 !
2720 ! end water treatment
2721 !
2722 !--- check for upwind convection
2723 ! iresult=0
2724 massfld=0.
2725
2726 ! call cup_direction2(i,j,dir,iact_old_gr, &
2727 ! massflx,iresult,1, &
2728 ! massfld, &
2729 ! ids,ide, jds,jde, kds,kde, &
2730 ! ims,ime, jms,jme, kms,kme, &
2731 ! its,ite, jts,jte, kts,kte )
2732 ! if(i.eq.ipr.and.j.eq.jpr.and.iedt.eq.1.and.ne.eq.1)then
2733 ! if(iedt.eq.1.and.ne.eq.1)then
2734 ! print *,massfld,ne,iedt,iens
2735 ! print *,xk(ne),xff_ens3(1),xff_ens3(2),xff_ens3(3)
2736 ! endif
2737 ! print *,i,j,massfld,aa0(i),aa1(i)
2738 IF(XK(ne).lt.0.and.xff0.gt.0.)iresultd=1
2739 iresulte=max(iresult,iresultd)
2740 iresulte=1
2741 if(iresulte.eq.1)then
2742 !
2743 !--- special treatment for stability closures
2744 !
2745
2746 if(xff0.gt.0.)then
2747 xf_ens(i,j,nall+1)=max(0.,-xff_ens3(1)/xk(ne)) &
2748 +massfld
2749 xf_ens(i,j,nall+2)=max(0.,-xff_ens3(2)/xk(ne)) &
2750 +massfld
2751 xf_ens(i,j,nall+3)=max(0.,-xff_ens3(3)/xk(ne)) &
2752 +massfld
2753 xf_ens(i,j,nall+13)=max(0.,-xff_ens3(13)/xk(ne)) &
2754 +massfld
2755 xf_ens(i,j,nall+14)=max(0.,-xff_ens3(14)/xk(ne)) &
2756 +massfld
2757 xf_ens(i,j,nall+15)=max(0.,-xff_ens3(15)/xk(ne)) &
2758 +massfld
2759 xf_ens(i,j,nall+16)=max(0.,-xff_ens3(16)/xk(ne)) &
2760 +massfld
2761 else
2762 xf_ens(i,j,nall+1)=massfld
2763 xf_ens(i,j,nall+2)=massfld
2764 xf_ens(i,j,nall+3)=massfld
2765 xf_ens(i,j,nall+13)=massfld
2766 xf_ens(i,j,nall+14)=massfld
2767 xf_ens(i,j,nall+15)=massfld
2768 xf_ens(i,j,nall+16)=massfld
2769 endif
2770 !
2771 !--- if iresult.eq.1, following independent of xff0
2772 !
2773 xf_ens(i,j,nall+4)=max(0.,xff_ens3(4) &
2774 +massfld)
2775 xf_ens(i,j,nall+5)=max(0.,xff_ens3(5) &
2776 +massfld)
2777 xf_ens(i,j,nall+6)=max(0.,xff_ens3(6) &
2778 +massfld)
2779 a1=max(1.e-3,pr_ens(i,j,nall+7))
2780 xf_ens(i,j,nall+7)=max(0.,xff_ens3(7) &
2781 /a1)
2782 a1=max(1.e-3,pr_ens(i,j,nall+8))
2783 xf_ens(i,j,nall+8)=max(0.,xff_ens3(8) &
2784 /a1)
2785 a1=max(1.e-3,pr_ens(i,j,nall+9))
2786 xf_ens(i,j,nall+9)=max(0.,xff_ens3(9) &
2787 /a1)
2788 if(XK(ne).lt.0.)then
2789 xf_ens(i,j,nall+10)=max(0., &
2790 -xff_ens3(10)/xk(ne)) &
2791 +massfld
2792 xf_ens(i,j,nall+11)=max(0., &
2793 -xff_ens3(11)/xk(ne)) &
2794 +massfld
2795 xf_ens(i,j,nall+12)=max(0., &
2796 -xff_ens3(12)/xk(ne)) &
2797 +massfld
2798 else
2799 xf_ens(i,j,nall+10)=massfld
2800 xf_ens(i,j,nall+11)=massfld
2801 xf_ens(i,j,nall+12)=massfld
2802 endif
2803 if(icoic.ge.1)then
2804 closure_n(i)=0.
2805 xf_ens(i,j,nall+1)=xf_ens(i,j,nall+icoic)
2806 xf_ens(i,j,nall+2)=xf_ens(i,j,nall+icoic)
2807 xf_ens(i,j,nall+3)=xf_ens(i,j,nall+icoic)
2808 xf_ens(i,j,nall+4)=xf_ens(i,j,nall+icoic)
2809 xf_ens(i,j,nall+5)=xf_ens(i,j,nall+icoic)
2810 xf_ens(i,j,nall+6)=xf_ens(i,j,nall+icoic)
2811 xf_ens(i,j,nall+7)=xf_ens(i,j,nall+icoic)
2812 xf_ens(i,j,nall+8)=xf_ens(i,j,nall+icoic)
2813 xf_ens(i,j,nall+9)=xf_ens(i,j,nall+icoic)
2814 xf_ens(i,j,nall+10)=xf_ens(i,j,nall+icoic)
2815 xf_ens(i,j,nall+11)=xf_ens(i,j,nall+icoic)
2816 xf_ens(i,j,nall+12)=xf_ens(i,j,nall+icoic)
2817 xf_ens(i,j,nall+13)=xf_ens(i,j,nall+icoic)
2818 xf_ens(i,j,nall+14)=xf_ens(i,j,nall+icoic)
2819 xf_ens(i,j,nall+15)=xf_ens(i,j,nall+icoic)
2820 xf_ens(i,j,nall+16)=xf_ens(i,j,nall+icoic)
2821 endif
2822 !
2823 ! replace 13-16 for now with other stab closures
2824 ! (13 gave problems for mass model)
2825 !
2826 ! xf_ens(i,j,nall+13)=xf_ens(i,j,nall+1)
2827 if(icoic.eq.0)xf_ens(i,j,nall+14)=xf_ens(i,j,nall+13)
2828 ! xf_ens(i,j,nall+15)=xf_ens(i,j,nall+11)
2829 ! xf_ens(i,j,nall+16)=xf_ens(i,j,nall+11)
2830 ! xf_ens(i,j,nall+7)=xf_ens(i,j,nall+4)
2831 ! xf_ens(i,j,nall+8)=xf_ens(i,j,nall+5)
2832 ! xf_ens(i,j,nall+9)=xf_ens(i,j,nall+6)
2833 !
2834 !--- store new for next time step
2835 !
2836 do nens3=1,maxens3
2837 massfln(i,j,nall+nens3)=edt(i) &
2838 *xf_ens(i,j,nall+nens3)
2839 massfln(i,j,nall+nens3)=max(0., &
2840 massfln(i,j,nall+nens3))
2841 enddo
2842 !
2843 !
2844 !--- do some more on the caps!!! ne=1 for 175, ne=2 for 100,....
2845 !
2846 ! do not care for caps here for closure groups 1 and 5,
2847 ! they are fine, do not turn them off here
2848 !
2849 !
2850 if(ne.eq.2.and.ierr2(i).gt.0)then
2851 xf_ens(i,j,nall+1) =0.
2852 xf_ens(i,j,nall+2) =0.
2853 xf_ens(i,j,nall+3) =0.
2854 xf_ens(i,j,nall+4) =0.
2855 xf_ens(i,j,nall+5) =0.
2856 xf_ens(i,j,nall+6) =0.
2857 xf_ens(i,j,nall+7) =0.
2858 xf_ens(i,j,nall+8) =0.
2859 xf_ens(i,j,nall+9) =0.
2860 xf_ens(i,j,nall+10)=0.
2861 xf_ens(i,j,nall+11)=0.
2862 xf_ens(i,j,nall+12)=0.
2863 xf_ens(i,j,nall+13)=0.
2864 xf_ens(i,j,nall+14)=0.
2865 xf_ens(i,j,nall+15)=0.
2866 xf_ens(i,j,nall+16)=0.
2867 massfln(i,j,nall+1)=0.
2868 massfln(i,j,nall+2)=0.
2869 massfln(i,j,nall+3)=0.
2870 massfln(i,j,nall+4)=0.
2871 massfln(i,j,nall+5)=0.
2872 massfln(i,j,nall+6)=0.
2873 massfln(i,j,nall+7)=0.
2874 massfln(i,j,nall+8)=0.
2875 massfln(i,j,nall+9)=0.
2876 massfln(i,j,nall+10)=0.
2877 massfln(i,j,nall+11)=0.
2878 massfln(i,j,nall+12)=0.
2879 massfln(i,j,nall+13)=0.
2880 massfln(i,j,nall+14)=0.
2881 massfln(i,j,nall+15)=0.
2882 massfln(i,j,nall+16)=0.
2883 endif
2884 if(ne.eq.3.and.ierr3(i).gt.0)then
2885 xf_ens(i,j,nall+1) =0.
2886 xf_ens(i,j,nall+2) =0.
2887 xf_ens(i,j,nall+3) =0.
2888 xf_ens(i,j,nall+4) =0.
2889 xf_ens(i,j,nall+5) =0.
2890 xf_ens(i,j,nall+6) =0.
2891 xf_ens(i,j,nall+7) =0.
2892 xf_ens(i,j,nall+8) =0.
2893 xf_ens(i,j,nall+9) =0.
2894 xf_ens(i,j,nall+10)=0.
2895 xf_ens(i,j,nall+11)=0.
2896 xf_ens(i,j,nall+12)=0.
2897 xf_ens(i,j,nall+13)=0.
2898 xf_ens(i,j,nall+14)=0.
2899 xf_ens(i,j,nall+15)=0.
2900 xf_ens(i,j,nall+16)=0.
2901 massfln(i,j,nall+1)=0.
2902 massfln(i,j,nall+2)=0.
2903 massfln(i,j,nall+3)=0.
2904 massfln(i,j,nall+4)=0.
2905 massfln(i,j,nall+5)=0.
2906 massfln(i,j,nall+6)=0.
2907 massfln(i,j,nall+7)=0.
2908 massfln(i,j,nall+8)=0.
2909 massfln(i,j,nall+9)=0.
2910 massfln(i,j,nall+10)=0.
2911 massfln(i,j,nall+11)=0.
2912 massfln(i,j,nall+12)=0.
2913 massfln(i,j,nall+13)=0.
2914 massfln(i,j,nall+14)=0.
2915 massfln(i,j,nall+15)=0.
2916 massfln(i,j,nall+16)=0.
2917 endif
2918
2919 endif
2920 350 continue
2921 ! ne=1, cap=175
2922 !
2923 nall=(iens-1)*maxens3*maxens*maxens2 &
2924 +(iedt-1)*maxens*maxens3
2925 ! ne=2, cap=100
2926 !
2927 nall2=(iens-1)*maxens3*maxens*maxens2 &
2928 +(iedt-1)*maxens*maxens3 &
2929 +(2-1)*maxens3
2930 xf_ens(i,j,nall+4) = xf_ens(i,j,nall2+4)
2931 xf_ens(i,j,nall+5) =xf_ens(i,j,nall2+5)
2932 xf_ens(i,j,nall+6) =xf_ens(i,j,nall2+6)
2933 xf_ens(i,j,nall+7) =xf_ens(i,j,nall2+7)
2934 xf_ens(i,j,nall+8) =xf_ens(i,j,nall2+8)
2935 xf_ens(i,j,nall+9) =xf_ens(i,j,nall2+9)
2936 xf_ens(i,j,nall+10)=xf_ens(i,j,nall2+10)
2937 xf_ens(i,j,nall+11)=xf_ens(i,j,nall2+11)
2938 xf_ens(i,j,nall+12)=xf_ens(i,j,nall2+12)
2939 go to 100
2940 endif
2941 elseif(ierr(i).ne.20.and.ierr(i).ne.0)then
2942 do n=1,ensdim
2943 xf_ens(i,j,n)=0.
2944 massfln(i,j,n)=0.
2945 enddo
2946 endif
2947 100 continue
2948
2949 END SUBROUTINE cup_forcing_ens
2950
2951
2952 SUBROUTINE cup_kbcon(cap_inc,iloop,k22,kbcon,he_cup,hes_cup, &
2953 ierr,kbmax,p_cup,cap_max, &
2954 ids,ide, jds,jde, kds,kde, &
2955 ims,ime, jms,jme, kms,kme, &
2956 its,ite, jts,jte, kts,kte )
2957
2958 IMPLICIT NONE
2959 !
2960
2961 ! only local wrf dimensions are need as of now in this routine
2962
2963 integer &
2964 ,intent (in ) :: &
2965 ids,ide, jds,jde, kds,kde, &
2966 ims,ime, jms,jme, kms,kme, &
2967 its,ite, jts,jte, kts,kte
2968 !
2969 !
2970 !
2971 ! ierr error value, maybe modified in this routine
2972 !
2973 real, dimension (its:ite,kts:kte) &
2974 ,intent (in ) :: &
2975 he_cup,hes_cup,p_cup
2976 real, dimension (its:ite) &
2977 ,intent (in ) :: &
2978 cap_max,cap_inc
2979 integer, dimension (its:ite) &
2980 ,intent (in ) :: &
2981 kbmax
2982 integer, dimension (its:ite) &
2983 ,intent (inout) :: &
2984 kbcon,k22,ierr
2985 integer &
2986 ,intent (in ) :: &
2987 iloop
2988 !
2989 ! local variables in this routine
2990 !
2991
2992 integer :: &
2993 i
2994 real :: &
2995 pbcdif,plus
2996 integer :: itf
2997
2998 itf=MIN(ite,ide-1)
2999 !
3000 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
3001 !
3002 DO 27 i=its,itf
3003 kbcon(i)=1
3004 IF(ierr(I).ne.0)GO TO 27
3005 KBCON(I)=K22(I)
3006 GO TO 32
3007 31 CONTINUE
3008 KBCON(I)=KBCON(I)+1
3009 IF(KBCON(I).GT.KBMAX(i)+2)THEN
3010 if(iloop.lt.4)ierr(i)=3
3011 ! if(iloop.lt.4)ierr(i)=997
3012 GO TO 27
3013 ENDIF
3014 32 CONTINUE
3015 IF(HE_cup(I,K22(I)).LT.HES_cup(I,KBCON(I)))GO TO 31
3016
3017 ! cloud base pressure and max moist static energy pressure
3018 ! i.e., the depth (in mb) of the layer of negative buoyancy
3019 if(KBCON(I)-K22(I).eq.1)go to 27
3020 PBCDIF=-P_cup(I,KBCON(I))+P_cup(I,K22(I))
3021 plus=max(25.,cap_max(i)-float(iloop-1)*cap_inc(i))
3022 if(iloop.eq.4)plus=cap_max(i)
3023 IF(PBCDIF.GT.plus)THEN
3024 K22(I)=K22(I)+1
3025 KBCON(I)=K22(I)
3026 GO TO 32
3027 ENDIF
3028 27 CONTINUE
3029
3030 END SUBROUTINE cup_kbcon
3031
3032
3033 SUBROUTINE cup_kbcon_cin(iloop,k22,kbcon,he_cup,hes_cup, &
3034 z,tmean,qes,ierr,kbmax,p_cup,cap_max,xl,cp, &
3035 ids,ide, jds,jde, kds,kde, &
3036 ims,ime, jms,jme, kms,kme, &
3037 its,ite, jts,jte, kts,kte )
3038
3039 IMPLICIT NONE
3040 !
3041
3042 ! only local wrf dimensions are need as of now in this routine
3043
3044 integer &
3045 ,intent (in ) :: &
3046 ids,ide, jds,jde, kds,kde, &
3047 ims,ime, jms,jme, kms,kme, &
3048 its,ite, jts,jte, kts,kte
3049 !
3050 !
3051 !
3052 ! ierr error value, maybe modified in this routine
3053 !
3054 real, dimension (its:ite,kts:kte) &
3055 ,intent (in ) :: &
3056 he_cup,hes_cup,p_cup,z,tmean,qes
3057 real, dimension (its:ite) &
3058 ,intent (in ) :: &
3059 cap_max
3060 real &
3061 ,intent (in ) :: &
3062 xl,cp
3063 integer, dimension (its:ite) &
3064 ,intent (in ) :: &
3065 kbmax
3066 integer, dimension (its:ite) &
3067 ,intent (inout) :: &
3068 kbcon,k22,ierr
3069 integer &
3070 ,intent (in ) :: &
3071 iloop
3072 !
3073 ! local variables in this routine
3074 !
3075
3076 integer :: &
3077 i,k
3078 real :: &
3079 cin,cin_max,dh,tprim,gamma
3080 !
3081 integer :: itf
3082
3083 itf=MIN(ite,ide-1)
3084 !
3085 !
3086
3087 !
3088 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
3089 !
3090 DO 27 i=its,itf
3091 cin_max=-cap_max(i)
3092 kbcon(i)=1
3093 cin = 0.
3094 IF(ierr(I).ne.0)GO TO 27
3095 KBCON(I)=K22(I)
3096 GO TO 32
3097 31 CONTINUE
3098 KBCON(I)=KBCON(I)+1
3099 IF(KBCON(I).GT.KBMAX(i)+2)THEN
3100 if(iloop.eq.1)ierr(i)=3
3101 ! if(iloop.eq.2)ierr(i)=997
3102 GO TO 27
3103 ENDIF
3104 32 CONTINUE
3105 dh = HE_cup(I,K22(I)) - HES_cup(I,KBCON(I))
3106 if (dh.lt. 0.) then
3107 GAMMA=(xl/cp)*(xl/(461.525*(Tmean(I,K22(i))**2)))*QES(I,K22(i))
3108 tprim = dh/(cp*(1.+gamma))
3109
3110 cin = cin + 9.8066 * tprim &
3111 *(z(i,k22(i))-z(i,k22(i)-1)) / tmean(i,k22(i))
3112 go to 31
3113 end if
3114
3115
3116 ! If negative energy in negatively buoyant layer
3117 ! exceeds convective inhibition (CIN) threshold,
3118 ! then set K22 level one level up and see if that
3119 ! will work.
3120
3121 IF(cin.lT.cin_max)THEN
3122 ! print *,i,cin,cin_max,k22(i),kbcon(i)
3123 K22(I)=K22(I)+1
3124 KBCON(I)=K22(I)
3125 GO TO 32
3126 ENDIF
3127 27 CONTINUE
3128
3129 END SUBROUTINE cup_kbcon_cin
3130
3131
3132 SUBROUTINE cup_ktop(ilo,dby,kbcon,ktop,ierr, &
3133 ids,ide, jds,jde, kds,kde, &
3134 ims,ime, jms,jme, kms,kme, &
3135 its,ite, jts,jte, kts,kte )
3136
3137 IMPLICIT NONE
3138 !
3139 ! on input
3140 !
3141
3142 ! only local wrf dimensions are need as of now in this routine
3143
3144 integer &
3145 ,intent (in ) :: &
3146 ids,ide, jds,jde, kds,kde, &
3147 ims,ime, jms,jme, kms,kme, &
3148 its,ite, jts,jte, kts,kte
3149 ! dby = buoancy term
3150 ! ktop = cloud top (output)
3151 ! ilo = flag
3152 ! ierr error value, maybe modified in this routine
3153 !
3154 real, dimension (its:ite,kts:kte) &
3155 ,intent (inout) :: &
3156 dby
3157 integer, dimension (its:ite) &
3158 ,intent (in ) :: &
3159 kbcon
3160 integer &
3161 ,intent (in ) :: &
3162 ilo
3163 integer, dimension (its:ite) &
3164 ,intent (out ) :: &
3165 ktop
3166 integer, dimension (its:ite) &
3167 ,intent (inout) :: &
3168 ierr
3169 !
3170 ! local variables in this routine
3171 !
3172
3173 integer :: &
3174 i,k
3175 !
3176 integer :: itf, ktf
3177
3178 itf=MIN(ite,ide-1)
3179 ktf=MIN(kte,kde-1)
3180 !
3181 !
3182 DO 42 i=its,itf
3183 ktop(i)=1
3184 IF(ierr(I).EQ.0)then
3185 DO 40 K=KBCON(I)+1,ktf-1
3186 IF(DBY(I,K).LE.0.)THEN
3187 KTOP(I)=K-1
3188 GO TO 41
3189 ENDIF
3190 40 CONTINUE
3191 if(ilo.eq.1)ierr(i)=5
3192 ! if(ilo.eq.2)ierr(i)=998
3193 GO TO 42
3194 41 CONTINUE
3195 do k=ktop(i)+1,ktf
3196 dby(i,k)=0.
3197 enddo
3198 endif
3199 42 CONTINUE
3200
3201 END SUBROUTINE cup_ktop
3202
3203
3204 SUBROUTINE cup_MAXIMI(ARRAY,KS,KE,MAXX,ierr, &
3205 ids,ide, jds,jde, kds,kde, &
3206 ims,ime, jms,jme, kms,kme, &
3207 its,ite, jts,jte, kts,kte )
3208
3209 IMPLICIT NONE
3210 !
3211 ! on input
3212 !
3213
3214 ! only local wrf dimensions are need as of now in this routine
3215
3216 integer &
3217 ,intent (in ) :: &
3218 ids,ide, jds,jde, kds,kde, &
3219 ims,ime, jms,jme, kms,kme, &
3220 its,ite, jts,jte, kts,kte
3221 ! array input array
3222 ! x output array with return values
3223 ! kt output array of levels
3224 ! ks,kend check-range
3225 real, dimension (its:ite,kts:kte) &
3226 ,intent (in ) :: &
3227 array
3228 integer, dimension (its:ite) &
3229 ,intent (in ) :: &
3230 ierr,ke
3231 integer &
3232 ,intent (in ) :: &
3233 ks
3234 integer, dimension (its:ite) &
3235 ,intent (out ) :: &
3236 maxx
3237 real, dimension (its:ite) :: &
3238 x
3239 real :: &
3240 xar
3241 integer :: &
3242 i,k
3243 integer :: itf
3244
3245 itf=MIN(ite,ide-1)
3246
3247 DO 200 i=its,itf
3248 MAXX(I)=KS
3249 if(ierr(i).eq.0)then
3250 X(I)=ARRAY(I,KS)
3251 !
3252 DO 100 K=KS,KE(i)
3253 XAR=ARRAY(I,K)
3254 IF(XAR.GE.X(I)) THEN
3255 X(I)=XAR
3256 MAXX(I)=K
3257 ENDIF
3258 100 CONTINUE
3259 endif
3260 200 CONTINUE
3261
3262 END SUBROUTINE cup_MAXIMI
3263
3264
3265 SUBROUTINE cup_minimi(ARRAY,KS,KEND,KT,ierr, &
3266 ids,ide, jds,jde, kds,kde, &
3267 ims,ime, jms,jme, kms,kme, &
3268 its,ite, jts,jte, kts,kte )
3269
3270 IMPLICIT NONE
3271 !
3272 ! on input
3273 !
3274
3275 ! only local wrf dimensions are need as of now in this routine
3276
3277 integer &
3278 ,intent (in ) :: &
3279 ids,ide, jds,jde, kds,kde, &
3280 ims,ime, jms,jme, kms,kme, &
3281 its,ite, jts,jte, kts,kte
3282 ! array input array
3283 ! x output array with return values
3284 ! kt output array of levels
3285 ! ks,kend check-range
3286 real, dimension (its:ite,kts:kte) &
3287 ,intent (in ) :: &
3288 array
3289 integer, dimension (its:ite) &
3290 ,intent (in ) :: &
3291 ierr,ks,kend
3292 integer, dimension (its:ite) &
3293 ,intent (out ) :: &
3294 kt
3295 real, dimension (its:ite) :: &
3296 x
3297 integer :: &
3298 i,k,kstop
3299
3300 integer :: itf
3301
3302 itf=MIN(ite,ide-1)
3303
3304 DO 200 i=its,itf
3305 KT(I)=KS(I)
3306 if(ierr(i).eq.0)then
3307 X(I)=ARRAY(I,KS(I))
3308 KSTOP=MAX(KS(I)+1,KEND(I))
3309 !
3310 DO 100 K=KS(I)+1,KSTOP
3311 IF(ARRAY(I,K).LT.X(I)) THEN
3312 X(I)=ARRAY(I,K)
3313 KT(I)=K
3314 ENDIF
3315 100 CONTINUE
3316 endif
3317 200 CONTINUE
3318
3319 END SUBROUTINE cup_MINIMI
3320
3321
3322 SUBROUTINE cup_output_ens(xf_ens,ierr,dellat,dellaq,dellaqc, &
3323 outtem,outq,outqc,pre,pw,xmb,ktop, &
3324 j,name,nx,nx2,iens,ierr2,ierr3,pr_ens, &
3325 maxens3,ensdim,massfln, &
3326 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
3327 APR_CAPMA,APR_CAPME,APR_CAPMI,closure_n,xland1, &
3328 ids,ide, jds,jde, kds,kde, &
3329 ims,ime, jms,jme, kms,kme, &
3330 its,ite, jts,jte, kts,kte)
3331
3332 IMPLICIT NONE
3333 !
3334 ! on input
3335 !
3336
3337 ! only local wrf dimensions are need as of now in this routine
3338
3339 integer &
3340 ,intent (in ) :: &
3341 ids,ide, jds,jde, kds,kde, &
3342 ims,ime, jms,jme, kms,kme, &
3343 its,ite, jts,jte, kts,kte
3344 integer, intent (in ) :: &
3345 j,ensdim,nx,nx2,iens,maxens3
3346 ! xf_ens = ensemble mass fluxes
3347 ! pr_ens = precipitation ensembles
3348 ! dellat = change of temperature per unit mass flux of cloud ensemble
3349 ! dellaq = change of q per unit mass flux of cloud ensemble
3350 ! dellaqc = change of qc per unit mass flux of cloud ensemble
3351 ! outtem = output temp tendency (per s)
3352 ! outq = output q tendency (per s)
3353 ! outqc = output qc tendency (per s)
3354 ! pre = output precip
3355 ! xmb = total base mass flux
3356 ! xfac1 = correction factor
3357 ! pw = pw -epsilon*pd (ensemble dependent)
3358 ! ierr error value, maybe modified in this routine
3359 !
3360 real, dimension (ims:ime,jms:jme,1:ensdim) &
3361 ,intent (inout) :: &
3362 xf_ens,pr_ens,massfln
3363 real, dimension (ims:ime,jms:jme) &
3364 ,intent (inout) :: &
3365 APR_GR,APR_W,APR_MC,APR_ST,APR_AS,APR_CAPMA, &
3366 APR_CAPME,APR_CAPMI
3367
3368 real, dimension (its:ite,kts:kte) &
3369 ,intent (out ) :: &
3370 outtem,outq,outqc
3371 real, dimension (its:ite) &
3372 ,intent (out ) :: &
3373 pre,xmb
3374 real, dimension (its:ite) &
3375 ,intent (inout ) :: &
3376 closure_n,xland1
3377 real, dimension (its:ite,kts:kte,1:ensdim) &
3378 ,intent (in ) :: &
3379 dellat,dellaqc,dellaq,pw
3380 integer, dimension (its:ite) &
3381 ,intent (in ) :: &
3382 ktop
3383 integer, dimension (its:ite) &
3384 ,intent (inout) :: &
3385 ierr,ierr2,ierr3
3386 !
3387 ! local variables in this routine
3388 !
3389
3390 integer :: &
3391 i,k,n,ncount
3392 real :: &
3393 outtes,ddtes,dtt,dtq,dtqc,dtpw,tuning,prerate,clos_wei
3394 real, dimension (its:ite) :: &
3395 xfac1
3396 real, dimension (its:ite):: &
3397 xmb_ske,xmb_ave,xmb_std,xmb_cur,xmbweight
3398 real, dimension (its:ite):: &
3399 pr_ske,pr_ave,pr_std,pr_cur
3400 real, dimension (its:ite,jts:jte):: &
3401 pr_gr,pr_w,pr_mc,pr_st,pr_as,pr_capma, &
3402 pr_capme,pr_capmi
3403
3404 !
3405 character *(*), intent (in) :: &
3406 name
3407 !
3408 integer :: itf, ktf
3409
3410 itf=MIN(ite,ide-1)
3411 ktf=MIN(kte,kde-1)
3412 tuning=0.
3413 !
3414 !
3415 DO k=kts,ktf
3416 do i=its,itf
3417 outtem(i,k)=0.
3418 outq(i,k)=0.
3419 outqc(i,k)=0.
3420 enddo
3421 enddo
3422 do i=its,itf
3423 pre(i)=0.
3424 xmb(i)=0.
3425 xfac1(i)=1.
3426 xmbweight(i)=1.
3427 enddo
3428 do i=its,itf
3429 IF(ierr(i).eq.0)then
3430 do n=(iens-1)*nx*nx2*maxens3+1,iens*nx*nx2*maxens3
3431 if(pr_ens(i,j,n).le.0.)then
3432 xf_ens(i,j,n)=0.
3433 endif
3434 enddo
3435 endif
3436 enddo
3437 !
3438 !--- calculate ensemble average mass fluxes
3439 !
3440 call massflx_stats(xf_ens,ensdim,nx2,nx,maxens3, &
3441 xmb_ave,xmb_std,xmb_cur,xmb_ske,j,ierr,1, &
3442 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
3443 APR_CAPMA,APR_CAPME,APR_CAPMI, &
3444 pr_gr,pr_w,pr_mc,pr_st,pr_as, &
3445 pr_capma,pr_capme,pr_capmi, &
3446 ids,ide, jds,jde, kds,kde, &
3447 ims,ime, jms,jme, kms,kme, &
3448 its,ite, jts,jte, kts,kte )
3449 call massflx_stats(pr_ens,ensdim,nx2,nx,maxens3, &
3450 pr_ave,pr_std,pr_cur,pr_ske,j,ierr,2, &
3451 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
3452 APR_CAPMA,APR_CAPME,APR_CAPMI, &
3453 pr_gr,pr_w,pr_mc,pr_st,pr_as, &
3454 pr_capma,pr_capme,pr_capmi, &
3455 ids,ide, jds,jde, kds,kde, &
3456 ims,ime, jms,jme, kms,kme, &
3457 its,ite, jts,jte, kts,kte )
3458 !
3459 !-- now do feedback
3460 !
3461 ddtes=200.
3462 ! if(name.eq.'shal')ddtes=200.
3463 do i=its,itf
3464 if(ierr(i).eq.0)then
3465 if(xmb_ave(i).le.0.)then
3466 ierr(i)=13
3467 xmb_ave(i)=0.
3468 endif
3469 ! xmb(i)=max(0.,xmb_ave(i))
3470 xmb(i)=max(.1*xmb_ave(i),xmb_ave(i)-tuning*xmb_std(i))
3471 ! --- Now use proper count of how many closures were actually
3472 ! used in cup_forcing_ens (including screening of some
3473 ! closures over water) to properly normalize xmb
3474 clos_wei=16./max(1.,closure_n(i))
3475 if (xland1(i).lt.0.5)xmb(i)=xmb(i)*clos_wei
3476 if(xmb(i).eq.0.)then
3477 ierr(i)=19
3478 endif
3479 if(xmb(i).gt.100.)then
3480 ierr(i)=19
3481 endif
3482 xfac1(i)=xmb(i)
3483
3484 endif
3485 xfac1(i)=xmb_ave(i)
3486 ENDDO
3487 DO k=kts,ktf
3488 do i=its,itf
3489 dtt=0.
3490 dtq=0.
3491 dtqc=0.
3492 dtpw=0.
3493 IF(ierr(i).eq.0.and.k.le.ktop(i))then
3494 do n=1,nx
3495 dtt=dtt+dellat(i,k,n)
3496 dtq=dtq+dellaq(i,k,n)
3497 dtqc=dtqc+dellaqc(i,k,n)
3498 dtpw=dtpw+pw(i,k,n)
3499 enddo
3500 outtes=dtt*XMB(I)*86400./float(nx)
3501 IF((OUTTES.GT.2.*ddtes.and.k.gt.2))THEN
3502 XMB(I)= 2.*ddtes/outtes * xmb(i)
3503 outtes=1.*ddtes
3504 endif
3505 if (outtes .lt. -ddtes) then
3506 XMB(I)= -ddtes/outtes * xmb(i)
3507 outtes=-ddtes
3508 endif
3509 if (outtes .gt. .5*ddtes.and.k.le.2) then
3510 XMB(I)= ddtes/outtes * xmb(i)
3511 outtes=.5*ddtes
3512 endif
3513 OUTTEM(I,K)=XMB(I)*dtt/float(nx)
3514 OUTQ(I,K)=XMB(I)*dtq/float(nx)
3515 OUTQC(I,K)=XMB(I)*dtqc/float(nx)
3516 PRE(I)=PRE(I)+XMB(I)*dtpw/float(nx)
3517 endif
3518 enddo
3519 enddo
3520 do i=its,itf
3521 if(ierr(i).eq.0)then
3522 prerate=pre(i)*3600.
3523 if(prerate.lt.0.1)then
3524 if(ierr2(i).gt.0.or.ierr3(i).gt.0)then
3525 pre(i)=0.
3526 ierr(i)=221
3527 do k=kts,ktf
3528 outtem(i,k)=0.
3529 outq(i,k)=0.
3530 outqc(i,k)=0.
3531 enddo
3532 do k=(iens-1)*nx*nx2*maxens3+1,iens*nx*nx2*maxens3
3533 massfln(i,j,k)=0.
3534 xf_ens(i,j,k)=0.
3535 enddo
3536 endif
3537 endif
3538
3539 endif
3540 ENDDO
3541
3542 do i=its,itf
3543 if(ierr(i).eq.0)then
3544 xfac1(i)=xmb(i)/xfac1(i)
3545 do k=(iens-1)*nx*nx2*maxens3+1,iens*nx*nx2*maxens3
3546 massfln(i,j,k)=massfln(i,j,k)*xfac1(i)
3547 xf_ens(i,j,k)=xf_ens(i,j,k)*xfac1(i)
3548 enddo
3549 endif
3550 ENDDO
3551
3552 END SUBROUTINE cup_output_ens
3553
3554
3555 SUBROUTINE cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
3556 kbcon,ktop,ierr, &
3557 ids,ide, jds,jde, kds,kde, &
3558 ims,ime, jms,jme, kms,kme, &
3559 its,ite, jts,jte, kts,kte )
3560
3561 IMPLICIT NONE
3562 !
3563 ! on input
3564 !
3565
3566 ! only local wrf dimensions are need as of now in this routine
3567
3568 integer &
3569 ,intent (in ) :: &
3570 ids,ide, jds,jde, kds,kde, &
3571 ims,ime, jms,jme, kms,kme, &
3572 its,ite, jts,jte, kts,kte
3573 ! aa0 cloud work function
3574 ! gamma_cup = gamma on model cloud levels
3575 ! t_cup = temperature (Kelvin) on model cloud levels
3576 ! dby = buoancy term
3577 ! zu= normalized updraft mass flux
3578 ! z = heights of model levels
3579 ! ierr error value, maybe modified in this routine
3580 !
3581 real, dimension (its:ite,kts:kte) &
3582 ,intent (in ) :: &
3583 z,zu,gamma_cup,t_cup,dby
3584 integer, dimension (its:ite) &
3585 ,intent (in ) :: &
3586 kbcon,ktop
3587 !
3588 ! input and output
3589 !
3590
3591
3592 integer, dimension (its:ite) &
3593 ,intent (inout) :: &
3594 ierr
3595 real, dimension (its:ite) &
3596 ,intent (out ) :: &
3597 aa0
3598 !
3599 ! local variables in this routine
3600 !
3601
3602 integer :: &
3603 i,k
3604 real :: &
3605 dz,da
3606 !
3607 integer :: itf, ktf
3608
3609 itf = MIN(ite,ide-1)
3610 ktf = MIN(kte,kde-1)
3611
3612 do i=its,itf
3613 aa0(i)=0.
3614 enddo
3615 DO 100 k=kts+1,ktf
3616 DO 100 i=its,itf
3617 IF(ierr(i).ne.0)GO TO 100
3618 IF(K.LE.KBCON(I))GO TO 100
3619 IF(K.Gt.KTOP(I))GO TO 100
3620 DZ=Z(I,K)-Z(I,K-1)
3621 da=zu(i,k)*DZ*(9.81/(1004.*( &
3622 (T_cup(I,K)))))*DBY(I,K-1)/ &
3623 (1.+GAMMA_CUP(I,K))
3624 IF(K.eq.KTOP(I).and.da.le.0.)go to 100
3625 AA0(I)=AA0(I)+da
3626 if(aa0(i).lt.0.)aa0(i)=0.
3627 100 continue
3628
3629 END SUBROUTINE cup_up_aa0
3630
3631
3632 SUBROUTINE cup_up_he(k22,hkb,z_cup,cd,entr,he_cup,hc, &
3633 kbcon,ierr,dby,he,hes_cup, &
3634 ids,ide, jds,jde, kds,kde, &
3635 ims,ime, jms,jme, kms,kme, &
3636 its,ite, jts,jte, kts,kte )
3637
3638 IMPLICIT NONE
3639 !
3640 ! on input
3641 !
3642
3643 ! only local wrf dimensions are need as of now in this routine
3644
3645 integer &
3646 ,intent (in ) :: &
3647 ids,ide, jds,jde, kds,kde, &
3648 ims,ime, jms,jme, kms,kme, &
3649 its,ite, jts,jte, kts,kte
3650 ! hc = cloud moist static energy
3651 ! hkb = moist static energy at originating level
3652 ! he = moist static energy on model levels
3653 ! he_cup = moist static energy on model cloud levels
3654 ! hes_cup = saturation moist static energy on model cloud levels
3655 ! dby = buoancy term
3656 ! cd= detrainment function
3657 ! z_cup = heights of model cloud levels
3658 ! entr = entrainment rate
3659 !
3660 real, dimension (its:ite,kts:kte) &
3661 ,intent (in ) :: &
3662 he,he_cup,hes_cup,z_cup,cd
3663 ! entr= entrainment rate
3664 real &
3665 ,intent (in ) :: &
3666 entr
3667 integer, dimension (its:ite) &
3668 ,intent (in ) :: &
3669 kbcon,k22
3670 !
3671 ! input and output
3672 !
3673
3674 ! ierr error value, maybe modified in this routine
3675
3676 integer, dimension (its:ite) &
3677 ,intent (inout) :: &
3678 ierr
3679
3680 real, dimension (its:ite,kts:kte) &
3681 ,intent (out ) :: &
3682 hc,dby
3683 real, dimension (its:ite) &
3684 ,intent (out ) :: &
3685 hkb
3686 !
3687 ! local variables in this routine
3688 !
3689
3690 integer :: &
3691 i,k
3692 real :: &
3693 dz
3694 !
3695 integer :: itf, ktf
3696
3697 itf = MIN(ite,ide-1)
3698 ktf = MIN(kte,kde-1)
3699 !
3700 !--- moist static energy inside cloud
3701 !
3702 do i=its,itf
3703 if(ierr(i).eq.0.)then
3704 hkb(i)=he_cup(i,k22(i))
3705 do k=1,k22(i)
3706 hc(i,k)=he_cup(i,k)
3707 DBY(I,K)=0.
3708 enddo
3709 do k=k22(i),kbcon(i)-1
3710 hc(i,k)=hkb(i)
3711 DBY(I,K)=0.
3712 enddo
3713 k=kbcon(i)
3714 hc(i,k)=hkb(i)
3715 DBY(I,Kbcon(i))=Hkb(I)-HES_cup(I,K)
3716 endif
3717 enddo
3718 do k=kts+1,ktf
3719 do i=its,itf
3720 if(k.gt.kbcon(i).and.ierr(i).eq.0.)then
3721 DZ=Z_cup(i,K)-Z_cup(i,K-1)
3722 HC(i,K)=(HC(i,K-1)*(1.-.5*CD(i,K)*DZ)+entr* &
3723 DZ*HE(i,K-1))/(1.+entr*DZ-.5*cd(i,k)*dz)
3724 DBY(I,K)=HC(I,K)-HES_cup(I,K)
3725 endif
3726 enddo
3727
3728 enddo
3729
3730 END SUBROUTINE cup_up_he
3731
3732
3733 SUBROUTINE cup_up_moisture(ierr,z_cup,qc,qrc,pw,pwav, &
3734 kbcon,ktop,cd,dby,mentr_rate,clw_all, &
3735 q,GAMMA_cup,zu,qes_cup,k22,qe_cup,xl, &
3736 ids,ide, jds,jde, kds,kde, &
3737 ims,ime, jms,jme, kms,kme, &
3738 its,ite, jts,jte, kts,kte )
3739
3740 IMPLICIT NONE
3741 !
3742 ! on input
3743 !
3744
3745 ! only local wrf dimensions are need as of now in this routine
3746
3747 integer &
3748 ,intent (in ) :: &
3749 ids,ide, jds,jde, kds,kde, &
3750 ims,ime, jms,jme, kms,kme, &
3751 its,ite, jts,jte, kts,kte
3752 ! cd= detrainment function
3753 ! q = environmental q on model levels
3754 ! qe_cup = environmental q on model cloud levels
3755 ! qes_cup = saturation q on model cloud levels
3756 ! dby = buoancy term
3757 ! cd= detrainment function
3758 ! zu = normalized updraft mass flux
3759 ! gamma_cup = gamma on model cloud levels
3760 ! mentr_rate = entrainment rate
3761 !
3762 real, dimension (its:ite,kts:kte) &
3763 ,intent (in ) :: &
3764 q,zu,gamma_cup,qe_cup,dby,qes_cup,z_cup,cd
3765 ! entr= entrainment rate
3766 real &
3767 ,intent (in ) :: &
3768 mentr_rate,xl
3769 integer, dimension (its:ite) &
3770 ,intent (in ) :: &
3771 kbcon,ktop,k22
3772 !
3773 ! input and output
3774 !
3775
3776 ! ierr error value, maybe modified in this routine
3777
3778 integer, dimension (its:ite) &
3779 ,intent (inout) :: &
3780 ierr
3781 ! qc = cloud q (including liquid water) after entrainment
3782 ! qrch = saturation q in cloud
3783 ! qrc = liquid water content in cloud after rainout
3784 ! pw = condensate that will fall out at that level
3785 ! pwav = totan normalized integrated condensate (I1)
3786 ! c0 = conversion rate (cloud to rain)
3787
3788 real, dimension (its:ite,kts:kte) &
3789 ,intent (out ) :: &
3790 qc,qrc,pw,clw_all
3791 real, dimension (its:ite) &
3792 ,intent (out ) :: &
3793 pwav
3794 !
3795 ! local variables in this routine
3796 !
3797
3798 integer :: &
3799 iall,i,k
3800 real :: &
3801 dh,qrch,c0,dz,radius
3802 !
3803 integer :: itf, ktf
3804
3805 itf = MIN(ite,ide-1)
3806 ktf = MIN(kte,kde-1)
3807
3808 iall=0
3809 c0=.002
3810 !
3811 !--- no precip for small clouds
3812 !
3813 if(mentr_rate.gt.0.)then
3814 radius=.2/mentr_rate
3815 if(radius.lt.900.)c0=0.
3816 ! if(radius.lt.900.)iall=0
3817 endif
3818 do i=its,itf
3819 pwav(i)=0.
3820 enddo
3821 do k=kts,ktf
3822 do i=its,itf
3823 pw(i,k)=0.
3824 qc(i,k)=qes_cup(i,k)
3825 clw_all(i,k)=0.
3826 qrc(i,k)=0.
3827 enddo
3828 enddo
3829 do i=its,itf
3830 if(ierr(i).eq.0.)then
3831 do k=k22(i),kbcon(i)-1
3832 qc(i,k)=qe_cup(i,k22(i))
3833 enddo
3834 endif
3835 enddo
3836
3837 DO 100 k=kts+1,ktf
3838 DO 100 i=its,itf
3839 IF(ierr(i).ne.0)GO TO 100
3840 IF(K.Lt.KBCON(I))GO TO 100
3841 IF(K.Gt.KTOP(I))GO TO 100
3842 DZ=Z_cup(i,K)-Z_cup(i,K-1)
3843 !
3844 !------ 1. steady state plume equation, for what could
3845 !------ be in cloud without condensation
3846 !
3847 !
3848 QC(i,K)=(QC(i,K-1)*(1.-.5*CD(i,K)*DZ)+mentr_rate* &
3849 DZ*Q(i,K-1))/(1.+mentr_rate*DZ-.5*cd(i,k)*dz)
3850 !
3851 !--- saturation in cloud, this is what is allowed to be in it
3852 !
3853 QRCH=QES_cup(I,K)+(1./XL)*(GAMMA_cup(i,k) &
3854 /(1.+GAMMA_cup(i,k)))*DBY(I,K)
3855 !
3856 !------- LIQUID WATER CONTENT IN cloud after rainout
3857 !
3858 clw_all(i,k)=QC(I,K)-QRCH
3859 QRC(I,K)=(QC(I,K)-QRCH)/(1.+C0*DZ)
3860 if(qrc(i,k).lt.0.)then
3861 qrc(i,k)=0.
3862 endif
3863 !
3864 !------- 3.Condensation
3865 !
3866 PW(i,k)=c0*dz*QRC(I,K)*zu(i,k)
3867 if(iall.eq.1)then
3868 qrc(i,k)=0.
3869 pw(i,k)=(QC(I,K)-QRCH)*zu(i,k)
3870 if(pw(i,k).lt.0.)pw(i,k)=0.
3871 endif
3872 !
3873 !----- set next level
3874 !
3875 QC(I,K)=QRC(I,K)+qrch
3876 !
3877 !--- integrated normalized ondensate
3878 !
3879 PWAV(I)=PWAV(I)+PW(I,K)
3880 100 CONTINUE
3881
3882 END SUBROUTINE cup_up_moisture
3883
3884
3885 SUBROUTINE cup_up_nms(zu,z_cup,entr,cd,kbcon,ktop,ierr,k22, &
3886 ids,ide, jds,jde, kds,kde, &
3887 ims,ime, jms,jme, kms,kme, &
3888 its,ite, jts,jte, kts,kte )
3889
3890 IMPLICIT NONE
3891
3892 !
3893 ! on input
3894 !
3895
3896 ! only local wrf dimensions are need as of now in this routine
3897
3898 integer &
3899 ,intent (in ) :: &
3900 ids,ide, jds,jde, kds,kde, &
3901 ims,ime, jms,jme, kms,kme, &
3902 its,ite, jts,jte, kts,kte
3903 ! cd= detrainment function
3904 real, dimension (its:ite,kts:kte) &
3905 ,intent (in ) :: &
3906 z_cup,cd
3907 ! entr= entrainment rate
3908 real &
3909 ,intent (in ) :: &
3910 entr
3911 integer, dimension (its:ite) &
3912 ,intent (in ) :: &
3913 kbcon,ktop,k22
3914 !
3915 ! input and output
3916 !
3917
3918 ! ierr error value, maybe modified in this routine
3919
3920 integer, dimension (its:ite) &
3921 ,intent (inout) :: &
3922 ierr
3923 ! zu is the normalized mass flux
3924
3925 real, dimension (its:ite,kts:kte) &
3926 ,intent (out ) :: &
3927 zu
3928 !
3929 ! local variables in this routine
3930 !
3931
3932 integer :: &
3933 i,k
3934 real :: &
3935 dz
3936 integer :: itf, ktf
3937
3938 itf = MIN(ite,ide-1)
3939 ktf = MIN(kte,kde-1)
3940 !
3941 ! initialize for this go around
3942 !
3943 do k=kts,ktf
3944 do i=its,itf
3945 zu(i,k)=0.
3946 enddo
3947 enddo
3948 !
3949 ! do normalized mass budget
3950 !
3951 do i=its,itf
3952 IF(ierr(I).eq.0)then
3953 do k=k22(i),kbcon(i)
3954 zu(i,k)=1.
3955 enddo
3956 DO K=KBcon(i)+1,KTOP(i)
3957 DZ=Z_cup(i,K)-Z_cup(i,K-1)
3958 ZU(i,K)=ZU(i,K-1)*(1.+(entr-cd(i,k))*DZ)
3959 enddo
3960 endif
3961 enddo
3962
3963 END SUBROUTINE cup_up_nms
3964
3965 !====================================================================
3966 SUBROUTINE gdinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
3967 MASS_FLUX,cp,restart, &
3968 P_QC,P_QI,P_FIRST_SCALAR, &
3969 RTHFTEN, RQVFTEN, &
3970 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
3971 APR_CAPMA,APR_CAPME,APR_CAPMI, &
3972 allowed_to_read, &
3973 ids, ide, jds, jde, kds, kde, &
3974 ims, ime, jms, jme, kms, kme, &
3975 its, ite, jts, jte, kts, kte )
3976 !--------------------------------------------------------------------
3977 IMPLICIT NONE
3978 !--------------------------------------------------------------------
3979 LOGICAL , INTENT(IN) :: restart,allowed_to_read
3980 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
3981 ims, ime, jms, jme, kms, kme, &
3982 its, ite, jts, jte, kts, kte
3983 INTEGER , INTENT(IN) :: P_FIRST_SCALAR, P_QI, P_QC
3984 REAL, INTENT(IN) :: cp
3985
3986 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
3987 RTHCUTEN, &
3988 RQVCUTEN, &
3989 RQCCUTEN, &
3990 RQICUTEN
3991
3992 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
3993 RTHFTEN, &
3994 RQVFTEN
3995
3996 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(OUT) :: &
3997 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
3998 APR_CAPMA,APR_CAPME,APR_CAPMI, &
3999 MASS_FLUX
4000
4001 INTEGER :: i, j, k, itf, jtf, ktf
4002
4003 jtf=min0(jte,jde-1)
4004 ktf=min0(kte,kde-1)
4005 itf=min0(ite,ide-1)
4006
4007 IF(.not.restart)THEN
4008 DO j=jts,jtf
4009 DO k=kts,ktf
4010 DO i=its,itf
4011 RTHCUTEN(i,k,j)=0.
4012 RQVCUTEN(i,k,j)=0.
4013 ENDDO
4014 ENDDO
4015 ENDDO
4016
4017 DO j=jts,jtf
4018 DO k=kts,ktf
4019 DO i=its,itf
4020 RTHFTEN(i,k,j)=0.
4021 RQVFTEN(i,k,j)=0.
4022 ENDDO
4023 ENDDO
4024 ENDDO
4025
4026 IF (P_QC .ge. P_FIRST_SCALAR) THEN
4027 DO j=jts,jtf
4028 DO k=kts,ktf
4029 DO i=its,itf
4030 RQCCUTEN(i,k,j)=0.
4031 ENDDO
4032 ENDDO
4033 ENDDO
4034 ENDIF
4035
4036 IF (P_QI .ge. P_FIRST_SCALAR) THEN
4037 DO j=jts,jtf
4038 DO k=kts,ktf
4039 DO i=its,itf
4040 RQICUTEN(i,k,j)=0.
4041 ENDDO
4042 ENDDO
4043 ENDDO
4044 ENDIF
4045
4046 DO j=jts,jtf
4047 DO i=its,itf
4048 mass_flux(i,j)=0.
4049 ENDDO
4050 ENDDO
4051
4052 ENDIF
4053 DO j=jts,jtf
4054 DO i=its,itf
4055 APR_GR(i,j)=0.
4056 APR_ST(i,j)=0.
4057 APR_W(i,j)=0.
4058 APR_MC(i,j)=0.
4059 APR_AS(i,j)=0.
4060 APR_CAPMA(i,j)=0.
4061 APR_CAPME(i,j)=0.
4062 APR_CAPMI(i,j)=0.
4063 ENDDO
4064 ENDDO
4065
4066 END SUBROUTINE gdinit
4067
4068
4069 SUBROUTINE massflx_stats(xf_ens,ensdim,maxens,maxens2,maxens3, &
4070 xt_ave,xt_std,xt_cur,xt_ske,j,ierr,itest, &
4071 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
4072 APR_CAPMA,APR_CAPME,APR_CAPMI, &
4073 pr_gr,pr_w,pr_mc,pr_st,pr_as, &
4074 pr_capma,pr_capme,pr_capmi, &
4075 ids,ide, jds,jde, kds,kde, &
4076 ims,ime, jms,jme, kms,kme, &
4077 its,ite, jts,jte, kts,kte)
4078
4079 IMPLICIT NONE
4080
4081 integer, intent (in ) :: &
4082 j,ensdim,maxens3,maxens,maxens2,itest
4083 INTEGER, INTENT(IN ) :: &
4084 ids,ide, jds,jde, kds,kde, &
4085 ims,ime, jms,jme, kms,kme, &
4086 its,ite, jts,jte, kts,kte
4087
4088
4089 real, dimension (its:ite) &
4090 , intent(inout) :: &
4091 xt_ave,xt_cur,xt_std,xt_ske
4092 integer, dimension (its:ite), intent (in) :: &
4093 ierr
4094 real, dimension (ims:ime,jms:jme,1:ensdim) &
4095 , intent(in ) :: &
4096 xf_ens
4097 real, dimension (ims:ime,jms:jme) &
4098 , intent(inout) :: &
4099 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
4100 APR_CAPMA,APR_CAPME,APR_CAPMI
4101 real, dimension (its:ite,jts:jte) &
4102 , intent(inout) :: &
4103 pr_gr,pr_w,pr_mc,pr_st,pr_as, &
4104 pr_capma,pr_capme,pr_capmi
4105
4106 !
4107 ! local stuff
4108 !
4109 real, dimension (its:ite , 1:maxens3 ) :: &
4110 x_ave,x_cur,x_std,x_ske
4111 real, dimension (its:ite , 1:maxens ) :: &
4112 x_ave_cap
4113
4114
4115 integer, dimension (1:maxens3) :: nc1
4116 integer :: i,k
4117 integer :: num,kk,num2,iedt
4118 real :: a3,a4
4119
4120 num=ensdim/maxens3
4121 num2=ensdim/maxens
4122 if(itest.eq.1)then
4123 do i=its,ite
4124 pr_gr(i,j) = 0.
4125 pr_w(i,j) = 0.
4126 pr_mc(i,j) = 0.
4127 pr_st(i,j) = 0.
4128 pr_as(i,j) = 0.
4129 pr_capma(i,j) = 0.
4130 pr_capme(i,j) = 0.
4131 pr_capmi(i,j) = 0.
4132 enddo
4133 endif
4134
4135 do k=1,maxens
4136 do i=its,ite
4137 x_ave_cap(i,k)=0.
4138 enddo
4139 enddo
4140 do k=1,maxens3
4141 do i=its,ite
4142 x_ave(i,k)=0.
4143 x_std(i,k)=0.
4144 x_ske(i,k)=0.
4145 x_cur(i,k)=0.
4146 enddo
4147 enddo
4148 do i=its,ite
4149 xt_ave(i)=0.
4150 xt_std(i)=0.
4151 xt_ske(i)=0.
4152 xt_cur(i)=0.
4153 enddo
4154 do kk=1,num
4155 do k=1,maxens3
4156 do i=its,ite
4157 if(ierr(i).eq.0)then
4158 x_ave(i,k)=x_ave(i,k)+xf_ens(i,j,maxens3*(kk-1)+k)
4159 endif
4160 enddo
4161 enddo
4162 enddo
4163 do iedt=1,maxens2
4164 do k=1,maxens
4165 do kk=1,maxens3
4166 do i=its,ite
4167 if(ierr(i).eq.0)then
4168 x_ave_cap(i,k)=x_ave_cap(i,k) &
4169 +xf_ens(i,j,maxens3*(k-1)+(iedt-1)*maxens*maxens3+kk)
4170 endif
4171 enddo
4172 enddo
4173 enddo
4174 enddo
4175 do k=1,maxens
4176 do i=its,ite
4177 if(ierr(i).eq.0)then
4178 x_ave_cap(i,k)=x_ave_cap(i,k)/float(num2)
4179 endif
4180 enddo
4181 enddo
4182
4183 do k=1,maxens3
4184 do i=its,ite
4185 if(ierr(i).eq.0)then
4186 x_ave(i,k)=x_ave(i,k)/float(num)
4187 endif
4188 enddo
4189 enddo
4190 do k=1,maxens3
4191 do i=its,ite
4192 if(ierr(i).eq.0)then
4193 xt_ave(i)=xt_ave(i)+x_ave(i,k)
4194 endif
4195 enddo
4196 enddo
4197 do i=its,ite
4198 if(ierr(i).eq.0)then
4199 xt_ave(i)=xt_ave(i)/float(maxens3)
4200 endif
4201 enddo
4202 !
4203 !--- now do std, skewness,curtosis
4204 !
4205 do kk=1,num
4206 do k=1,maxens3
4207 do i=its,ite
4208 if(ierr(i).eq.0.and.x_ave(i,k).gt.0.)then
4209 ! print *,i,j,k,kk,x_std(i,k),xf_ens(i,j,maxens3*(kk-1)+k),x_ave(i,k)
4210 x_std(i,k)=x_std(i,k)+(xf_ens(i,j,maxens3*(kk-1)+k)-x_ave(i,k))**2
4211 x_ske(i,k)=x_ske(i,k)+(xf_ens(i,j,maxens3*(kk-1)+k)-x_ave(i,k))**3
4212 x_cur(i,k)=x_cur(i,k)+(xf_ens(i,j,maxens3*(kk-1)+k)-x_ave(i,k))**4
4213 endif
4214 enddo
4215 enddo
4216 enddo
4217 do k=1,maxens3
4218 do i=its,ite
4219 if(ierr(i).eq.0.and.xt_ave(i).gt.0.)then
4220 xt_std(i)=xt_std(i)+(x_ave(i,k)-xt_ave(i))**2
4221 xt_ske(i)=xt_ske(i)+(x_ave(i,k)-xt_ave(i))**3
4222 xt_cur(i)=xt_cur(i)+(x_ave(i,k)-xt_ave(i))**4
4223 endif
4224 enddo
4225 enddo
4226 do k=1,maxens3
4227 do i=its,ite
4228 if(ierr(i).eq.0.and.x_std(i,k).gt.0.)then
4229 x_std(i,k)=x_std(i,k)/float(num)
4230 a3=max(1.e-6,x_std(i,k))
4231 x_std(i,k)=sqrt(a3)
4232 a3=max(1.e-6,x_std(i,k)**3)
4233 a4=max(1.e-6,x_std(i,k)**4)
4234 x_ske(i,k)=x_ske(i,k)/float(num)/a3
4235 x_cur(i,k)=x_cur(i,k)/float(num)/a4
4236 endif
4237 ! print*,' '
4238 ! print*,'Some statistics at gridpoint i,j, ierr',i,j,ierr(i)
4239 ! print*,'statistics for closure number ',k
4240 ! print*,'Average= ',x_ave(i,k),' Std= ',x_std(i,k)
4241 ! print*,'Skewness= ',x_ske(i,k),' Curtosis= ',x_cur(i,k)
4242 ! print*,' '
4243
4244 enddo
4245 enddo
4246 do i=its,ite
4247 if(ierr(i).eq.0.and.xt_std(i).gt.0.)then
4248 xt_std(i)=xt_std(i)/float(maxens3)
4249 a3=max(1.e-6,xt_std(i))
4250 xt_std(i)=sqrt(a3)
4251 a3=max(1.e-6,xt_std(i)**3)
4252 a4=max(1.e-6,xt_std(i)**4)
4253 xt_ske(i)=xt_ske(i)/float(maxens3)/a3
4254 xt_cur(i)=xt_cur(i)/float(maxens3)/a4
4255 ! print*,' '
4256 ! print*,'Total ensemble independent statistics at i =',i
4257 ! print*,'Average= ',xt_ave(i),' Std= ',xt_std(i)
4258 ! print*,'Skewness= ',xt_ske(i),' Curtosis= ',xt_cur(i)
4259 ! print*,' '
4260 !
4261 ! first go around: store massflx for different closures/caps
4262 !
4263 if(itest.eq.1)then
4264 pr_gr(i,j) = .333*(x_ave(i,1)+x_ave(i,2)+x_ave(i,3))
4265 pr_w(i,j) = .333*(x_ave(i,4)+x_ave(i,5)+x_ave(i,6))
4266 pr_mc(i,j) = .333*(x_ave(i,7)+x_ave(i,8)+x_ave(i,9))
4267 pr_st(i,j) = .333*(x_ave(i,10)+x_ave(i,11)+x_ave(i,12))
4268 pr_as(i,j) = .25*(x_ave(i,13)+x_ave(i,14)+x_ave(i,15) &
4269 + x_ave(i,16))
4270 pr_capma(i,j) = x_ave_cap(i,1)
4271 pr_capme(i,j) = x_ave_cap(i,2)
4272 pr_capmi(i,j) = x_ave_cap(i,3)
4273 !
4274 ! second go around: store preciprates (mm/hour) for different closures/caps
4275 !
4276 else if (itest.eq.2)then
4277 APR_GR(i,j)=.333*(x_ave(i,1)+x_ave(i,2)+x_ave(i,3))* &
4278 3600.*pr_gr(i,j) +APR_GR(i,j)
4279 APR_W(i,j)=.333*(x_ave(i,4)+x_ave(i,5)+x_ave(i,6))* &
4280 3600.*pr_w(i,j) +APR_W(i,j)
4281 APR_MC(i,j)=.333*(x_ave(i,7)+x_ave(i,8)+x_ave(i,9))* &
4282 3600.*pr_mc(i,j) +APR_MC(i,j)
4283 APR_ST(i,j)=.333*(x_ave(i,10)+x_ave(i,11)+x_ave(i,12))* &
4284 3600.*pr_st(i,j) +APR_ST(i,j)
4285 APR_AS(i,j)=.25*(x_ave(i,13)+x_ave(i,14)+x_ave(i,15) &
4286 + x_ave(i,16))* &
4287 3600.*pr_as(i,j) +APR_AS(i,j)
4288 APR_CAPMA(i,j) = x_ave_cap(i,1)* &
4289 3600.*pr_capma(i,j) +APR_CAPMA(i,j)
4290 APR_CAPME(i,j) = x_ave_cap(i,2)* &
4291 3600.*pr_capme(i,j) +APR_CAPME(i,j)
4292 APR_CAPMI(i,j) = x_ave_cap(i,3)* &
4293 3600.*pr_capmi(i,j) +APR_CAPMI(i,j)
4294 endif
4295 endif
4296 enddo
4297
4298 END SUBROUTINE massflx_stats
4299
4300
4301 SUBROUTINE neg_check(dt,q,outq,outt,outqc,pret,its,ite,kts,kte,itf,ktf)
4302
4303 INTEGER, INTENT(IN ) :: its,ite,kts,kte,itf,ktf
4304
4305 real, dimension (its:ite,kts:kte ) , &
4306 intent(inout ) :: &
4307 q,outq,outt,outqc
4308 real, dimension (its:ite ) , &
4309 intent(inout ) :: &
4310 pret
4311 real &
4312 ,intent (in ) :: &
4313 dt
4314 real :: thresh,qmem,qmemf,qmem2,qtest,qmem1
4315 !
4316 ! first do check on vertical heating rate
4317 !
4318 thresh=200.01
4319 do i=its,itf
4320 qmemf=1.
4321 qmem=0.
4322 do k=kts,ktf
4323 qmem=outt(i,k)*86400.
4324 if(qmem.gt.2.*thresh)then
4325 qmem2=2.*thresh/qmem
4326 qmemf=min(qmemf,qmem2)
4327 !
4328 !
4329 ! print *,'1',' adjusted massflux by factor ',i,k,qmem,qmem2,qmemf
4330 endif
4331 if(qmem.lt.-thresh)then
4332 qmem2=-thresh/qmem
4333 qmemf=min(qmemf,qmem2)
4334 !
4335 !
4336 ! print *,'2',' adjusted massflux by factor ',i,k,qmem,qmem2,qmemf
4337 endif
4338 enddo
4339 do k=kts,ktf
4340 outq(i,k)=outq(i,k)*qmemf
4341 outt(i,k)=outt(i,k)*qmemf
4342 outqc(i,k)=outqc(i,k)*qmemf
4343 enddo
4344 pret(i)=pret(i)*qmemf
4345 enddo
4346 !
4347 ! check whether routine produces negative q's. This can happen, since
4348 ! tendencies are calculated based on forced q's. This should have no
4349 ! influence on conservation properties, it scales linear through all
4350 ! tendencies
4351 !
4352 thresh=1.e-10
4353 do i=its,itf
4354 qmemf=1.
4355 do k=kts,ktf
4356 qmem=outq(i,k)
4357 if(abs(qmem).gt.0.)then
4358 qtest=q(i,k)+outq(i,k)*dt
4359 if(qtest.lt.thresh)then
4360 !
4361 ! qmem2 would be the maximum allowable tendency
4362 !
4363 qmem1=outq(i,k)
4364 qmem2=(thresh-q(i,k))/dt
4365 qmemf=min(qmemf,qmem2/qmem1)
4366 ! qmemf=max(0.,qmemf)
4367 ! print *,'4 adjusted tendencies ',i,k,qmem,qmem2,qmemf
4368 endif
4369 endif
4370 enddo
4371 do k=kts,ktf
4372 outq(i,k)=outq(i,k)*qmemf
4373 outt(i,k)=outt(i,k)*qmemf
4374 outqc(i,k)=outqc(i,k)*qmemf
4375 enddo
4376 pret(i)=pret(i)*qmemf
4377 enddo
4378
4379 END SUBROUTINE neg_check
4380
4381
4382 !-------------------------------------------------------
4383 END MODULE module_cu_gd