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