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