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