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