| blob | 762d87cb8d74899d059c340cd7527b2799312692 |
1 !WRF:MODEL_LAYER:PHYSICS
2 !
4 MODULE module_cu_g3
6 CONTAINS
8 !-------------------------------------------------------------
9 SUBROUTINE G3DRV( &
10 DT,itimestep,DX &
11 ,rho,RAINCV,PRATEC &
12 ,U,V,t,W,q,p,pi &
13 ,dz8w,p8w,XLV,CP,G,r_v &
14 ,STEPCU,htop,hbot &
15 ,CU_ACT_FLAG,warm_rain &
16 ,APR_GR,APR_W,APR_MC,APR_ST,APR_AS &
17 ,APR_CAPMA,APR_CAPME,APR_CAPMI &
18 ,MASS_FLUX,XF_ENS,PR_ENS,HT,XLAND,gsw,edt_out &
19 ,GDC,GDC2 ,kpbl,k22_shallow,kbcon_shallow &
20 ,ktop_shallow,xmb_shallow &
21 ,cugd_tten,cugd_qvten ,cugd_qcten &
22 ,cugd_ttens,cugd_qvtens,cugd_avedx,imomentum &
23 ,ensdim,maxiens,maxens,maxens2,maxens3,ichoice &
24 ,ishallow_g3,ids,ide, jds,jde, kds,kde &
25 ,ims,ime, jms,jme, kms,kme &
26 ,ips,ipe, jps,jpe, kps,kpe &
27 ,its,ite, jts,jte, kts,kte &
28 ,periodic_x,periodic_y &
29 ,RQVCUTEN,RQCCUTEN,RQICUTEN &
30 ,RQVFTEN,RTHFTEN,RTHCUTEN &
31 ,rqvblten,rthblten &
32 ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
33 #if ( WRF_DFI_RADAR == 1 )
34 ! Optional CAP suppress option
35 ,do_capsuppress,cap_suppress_loc &
36 #endif
37 )
38 !-------------------------------------------------------------
39 IMPLICIT NONE
40 !-------------------------------------------------------------
41 INTEGER, INTENT(IN ) :: &
42 ids,ide, jds,jde, kds,kde, &
43 ims,ime, jms,jme, kms,kme, &
44 ips,ipe, jps,jpe, kps,kpe, &
45 its,ite, jts,jte, kts,kte
46 LOGICAL periodic_x,periodic_y
47 integer, parameter :: ens4_spread = 3 ! max(3,cugd_avedx)
48 integer, parameter :: ens4=ens4_spread*ens4_spread
50 integer, intent (in ) :: &
51 ensdim,maxiens,maxens,maxens2,maxens3,ichoice
53 INTEGER, INTENT(IN ) :: STEPCU, ITIMESTEP,cugd_avedx, &
54 ishallow_g3,imomentum
55 LOGICAL, INTENT(IN ) :: warm_rain
57 REAL, INTENT(IN ) :: XLV, R_v
58 REAL, INTENT(IN ) :: CP,G
60 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
61 INTENT(IN ) :: &
62 U, &
63 V, &
64 W, &
65 pi, &
66 t, &
67 q, &
68 p, &
69 dz8w, &
70 p8w, &
71 rho
72 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
73 OPTIONAL , &
74 INTENT(INOUT ) :: &
75 GDC,GDC2
77 REAL, DIMENSION( ims:ime , jms:jme ),INTENT(IN) :: GSW,HT,XLAND
78 INTEGER, DIMENSION( ims:ime , jms:jme ),INTENT(IN) :: KPBL
79 INTEGER, DIMENSION( ims:ime , jms:jme ),INTENT(INOUT) :: k22_shallow, &
80 kbcon_shallow,ktop_shallow
81 !
82 REAL, INTENT(IN ) :: DT, DX
83 !
85 REAL, DIMENSION( ims:ime , jms:jme ), &
86 INTENT(INOUT) :: pratec,RAINCV, MASS_FLUX, &
87 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
88 edt_out,APR_CAPMA,APR_CAPME,APR_CAPMI, &
89 htop,hbot,xmb_shallow
90 !+lxz
91 ! REAL, DIMENSION( ims:ime , jms:jme ) :: & !, INTENT(INOUT) :: &
92 ! HTOP, &! highest model layer penetrated by cumulus since last reset in radiation_driver
93 ! HBOT ! lowest model layer penetrated by cumulus since last reset in radiation_driver
94 ! ! HBOT>HTOP follow physics leveling convention
96 LOGICAL, DIMENSION( ims:ime , jms:jme ), &
97 INTENT(INOUT) :: CU_ACT_FLAG
99 !
100 ! Optionals
101 !
102 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
103 OPTIONAL, &
104 INTENT(INOUT) :: RTHFTEN, &
105 cugd_tten,cugd_qvten,cugd_qcten, &
106 cugd_ttens,cugd_qvtens, &
107 RQVFTEN
109 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
110 OPTIONAL, &
111 INTENT(INOUT) :: &
112 RTHCUTEN, &
113 RQVCUTEN, &
114 RQVBLTEN, &
115 RTHBLTEN, &
116 RQCCUTEN, &
117 RQICUTEN
118 !
119 ! Flags relating to the optional tendency arrays declared above
120 ! Models that carry the optional tendencies will provdide the
121 ! optional arguments at compile time; these flags all the model
122 ! to determine at run-time whether a particular tracer is in
123 ! use or not.
124 !
125 LOGICAL, OPTIONAL :: &
126 F_QV &
127 ,F_QC &
128 ,F_QR &
129 ,F_QI &
130 ,F_QS
133 #if ( WRF_DFI_RADAR == 1 )
134 !
135 ! option of cap suppress:
136 ! do_capsuppress = 1 do
137 ! do_capsuppress = other don't
138 !
139 !
140 INTEGER, INTENT(IN ) ,OPTIONAL :: do_capsuppress
141 REAL, DIMENSION( ims:ime, jms:jme ),INTENT(IN ),OPTIONAL :: cap_suppress_loc
142 REAL, DIMENSION( its:ite ) :: cap_suppress_j
143 #endif
146 ! LOCAL VARS
147 real, dimension(ims:ime,jms:jme,1:ensdim),intent(inout) :: &
148 xf_ens,pr_ens
149 real, dimension ( its:ite , jts:jte , 1:ensdim) :: &
150 massflni,xfi_ens,pri_ens
151 REAL, DIMENSION( its:ite , jts:jte ) :: MASSI_FLX, &
152 APRi_GR,APRi_W,APRi_MC,APRi_ST,APRi_AS, &
153 edti_out,APRi_CAPMA,APRi_CAPME,APRi_CAPMI,gswi
154 real, dimension (its:ite,kts:kte) :: &
155 SUBT,SUBQ,OUTT,OUTQ,OUTQC,phh,subm,cupclw,dhdt, &
156 outts,outqs
157 real, dimension (its:ite) :: &
158 pret, ter11, aa0, fp,xlandi
159 !+lxz
160 integer, dimension (its:ite) :: &
161 kbcon, ktop,kpbli,k22s,kbcons,ktops
162 !.lxz
163 integer, dimension (its:ite,jts:jte) :: &
164 iact_old_gr
165 integer :: iens,ibeg,iend,jbeg,jend,n,nn,ens4n
166 integer :: ibegh,iendh,jbegh,jendh
167 integer :: ibegc,iendc,jbegc,jendc
169 !
170 ! basic environmental input includes moisture convergence (mconv)
171 ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
172 ! convection for this call only and at that particular gridpoint
173 !
174 real, dimension (its:ite,kts:kte) :: &
175 T2d,q2d,PO,P2d,US,VS,tn,qo,tshall,qshall
176 real, dimension (ips-2:ipe+2,kps:kpe,jps-2:jpe+2) :: &
177 ave_f_t,ave_f_q
178 real, dimension (its:ite,kts:kte,1:ens4) :: &
179 omeg,tx,qx
180 real, dimension (its:ite) :: &
181 Z1,PSUR,AAEQ,direction,cuten,umean,vmean,pmean,xmbs
182 real, dimension (its:ite,1:ens4) :: &
183 mconv
185 INTEGER :: i,j,k,ICLDCK,ipr,jpr
186 REAL :: tcrit,tscl_KF,dp,dq,sub_spread,subcenter
187 INTEGER :: itf,jtf,ktf,iss,jss,nbegin,nend
188 INTEGER :: high_resolution
189 REAL :: rkbcon,rktop !-lxz
190 ! ruc variable
191 real, dimension (its:ite) :: tkm
193 ! A. Betts for shallow convection: suggestion for the KF timescale < DELTAX / 25 m/s
194 tscl_kf=dx/25.
195 !
196 ! write(0,*)'ishallow = ',ishallow_g3
197 high_resolution=0
198 if(cugd_avedx.gt.1) high_resolution=1
199 subcenter=0.
200 ! subcenter=1./float(cugd_avedx)
201 sub_spread=max(1.,float(cugd_avedx*cugd_avedx-1))
202 sub_spread=(1.-subcenter)/sub_spread
203 iens=1
204 ipr=37
205 jpr=1
206 ipr=0
207 jpr=0
208 ! if(itimestep.eq.8)then
209 ! ipr=37
210 ! jpr=16
211 ! endif
212 IF ( periodic_x ) THEN
213 ibeg=max(its,ids)
214 iend=min(ite,ide-1)
215 ibegc=max(its,ids)
216 iendc=min(ite,ide-1)
217 ELSE
218 ibeg=max(its,ids)
219 iend=min(ite,ide-1)
220 ibegc=max(its,ids+4)
221 iendc=min(ite,ide-5)
222 END IF
223 IF ( periodic_y ) THEN
224 jbeg=max(jts,jds)
225 jend=min(jte,jde-1)
226 jbegc=max(jts,jds)
227 jendc=min(jte,jde-1)
228 ELSE
229 jbeg=max(jts,jds)
230 jend=min(jte,jde-1)
231 jbegc=max(jts,jds+4)
232 jendc=min(jte,jde-5)
233 END IF
234 do j=jts,jte
235 do i=its,ite
236 ktop_shallow(i,j)=0
237 xmb_shallow(i,j)=0
238 enddo
239 enddo
240 tcrit=258.
241 ave_f_t=0.
242 ave_f_q=0.
244 itf=MIN(ite,ide-1)
245 ktf=MIN(kte,kde-1)
246 jtf=MIN(jte,jde-1)
247 !
248 #if ( EM_CORE == 1 )
249 if(high_resolution.eq.1)then
250 !
251 ! calculate these on the halo...the incominh tendencies have been exchanged on a 24pt halo
252 ! only neede for high resolution run
253 !
254 ibegh=its
255 jbegh=jts
256 iendh=ite
257 jendh=jte
258 if(its.eq.ips)ibegh=max(its-1,ids)
259 if(jts.eq.jps)jbegh=max(jts-1,jds)
260 if(jte.eq.jpe)jendh=min(jte+1,jde-1)
261 if(ite.eq.ipe)iendh=min(ite+1,ide-1)
262 DO J = jbegh,jendh
263 DO k= kts,ktf
264 DO I= ibegh,iendh
265 ave_f_t(i,k,j)=(rthften(i-1,k,j-1)+rthften(i-1,k,j) + rthften(i-1,k,j+1)+ &
266 rthften(i,k,j-1) +rthften(i,k,j) +rthften(i,k,j+1)+ &
267 rthften(i+1,k,j-1) +rthften(i+1,k,j) +rthften(i+1,k,j+1))/9.
268 ave_f_q(i,k,j)=(rqvften(i-1,k,j-1)+rqvften(i-1,k,j) + rqvften(i-1,k,j+1)+ &
269 rqvften(i,k,j-1) +rqvften(i,k,j) +rqvften(i,k,j+1)+ &
270 rqvften(i+1,k,j-1) +rqvften(i+1,k,j) +rqvften(i+1,k,j+1))/9.
271 ! ave_f_t(i,k,j)=rthften(i,k,j)
272 ! ave_f_q(i,k,j)=rqvften(i,k,j)
273 ENDDO
274 ENDDO
275 ENDDO
276 endif
277 #endif
278 DO 100 J = jts,jtf
279 DO n= 1,ensdim
280 DO I= its,itf
281 xfi_ens(i,j,n)=0.
282 pri_ens(i,j,n)=0.
283 ! xfi_ens(i,j,n)=xf_ens(i,j,n)
284 ! pri_ens(i,j,n)=pr_ens(i,j,n)
285 ENDDO
286 ENDDO
287 DO I= its,itf
288 kbcon(i)=0
289 ktop(i)=0
290 tkm(i)=0.
291 HBOT(I,J) =REAL(KTE)
292 HTOP(I,J) =REAL(KTS)
293 iact_old_gr(i,j)=0
294 mass_flux(i,j)=0.
295 massi_flx(i,j)=0.
296 raincv(i,j)=0.
297 pratec (i,j)=0.
298 edt_out(i,j)=0.
299 edti_out(i,j)=0.
300 gswi(i,j)=gsw(i,j)
301 xlandi(i)=xland(i,j)
302 APRi_GR(i,j)=apr_gr(i,j)
303 APRi_w(i,j)=apr_w(i,j)
304 APRi_mc(i,j)=apr_mc(i,j)
305 APRi_st(i,j)=apr_st(i,j)
306 APRi_as(i,j)=apr_as(i,j)
307 APRi_capma(i,j)=apr_capma(i,j)
308 APRi_capme(i,j)=apr_capme(i,j)
309 APRi_capmi(i,j)=apr_capmi(i,j)
310 CU_ACT_FLAG(i,j) = .true.
311 ENDDO
312 do k=kts,kte
313 DO I= its,itf
314 cugd_tten(i,k,j)=0.
315 cugd_ttens(i,k,j)=0.
316 cugd_qvten(i,k,j)=0.
317 cugd_qvtens(i,k,j)=0.
318 cugd_qcten(i,k,j)=0.
319 ENDDO
320 ENDDO
321 DO n=1,ens4
322 DO I= its,itf
323 mconv(i,n)=0.
324 ENDDO
325 do k=kts,kte
326 DO I= its,itf
327 omeg(i,k,n)=0.
328 tx(i,k,n)=0.
329 qx(i,k,n)=0.
330 ENDDO
331 ENDDO
332 ENDDO
333 DO k=1,ensdim
334 DO I= its,itf
335 massflni(i,j,k)=0.
336 ENDDO
337 ENDDO
338 ! put hydrostatic pressure on half levels
339 DO K=kts,ktf
340 DO I=ITS,ITF
341 phh(i,k) = p(i,k,j)
342 ENDDO
343 ENDDO
345 DO I=ITS,ITF
346 PSUR(I)=p8w(I,1,J)*.01
347 ! PSUR(I)=p(I,1,J)*.01
348 TER11(I)=HT(i,j)
349 aaeq(i)=0.
350 direction(i)=0.
351 pret(i)=0.
352 umean(i)=0.
353 vmean(i)=0.
354 pmean(i)=0.
355 kpbli(i)=kpbl(i,j)
356 ENDDO
357 DO K=kts,ktf
358 DO I=ITS,ITF
359 po(i,k)=phh(i,k)*.01
360 subm(i,k)=0.
361 P2d(I,K)=PO(i,k)
362 US(I,K) =u(i,k,j)
363 VS(I,K) =v(i,k,j)
364 T2d(I,K)=t(i,k,j)
365 q2d(I,K)=q(i,k,j)
366 IF(Q2d(I,K).LT.1.E-08)Q2d(I,K)=1.E-08
367 SUBT(I,K)=0.
368 SUBQ(I,K)=0.
369 OUTT(I,K)=0.
370 OUTQ(I,K)=0.
371 OUTQC(I,K)=0.
372 OUTTS(I,K)=0.
373 OUTQS(I,K)=0.
374 TN(I,K)=t2d(i,k)+RTHFTEN(i,k,j)*dt
375 QO(I,K)=q2d(i,k)+RQVFTEN(i,k,j)*dt
376 TSHALL(I,K)=t2d(i,k)+RTHBLTEN(i,k,j)*pi(i,k,j)*dt
377 DHDT(I,K)=cp*RTHBLTEN(i,k,j)*pi(i,k,j)+ XLV*RQVBLTEN(i,k,j)
378 QSHALL(I,K)=q2d(i,k)+RQVBLTEN(i,k,j)*dt
379 if(high_resolution.eq.1)then
380 TN(I,K)=t2d(i,k)+ave_f_t(i,k,j)*dt
381 QO(I,K)=q2d(i,k)+ave_f_q(i,k,j)*dt
382 endif
383 IF(TN(I,K).LT.200.)TN(I,K)=T2d(I,K)
384 IF(QO(I,K).LT.1.E-08)QO(I,K)=1.E-08
385 ENDDO
386 ENDDO
387 ens4n=0
388 nbegin=0
389 nend=0
390 if(ens4_spread.gt.1)then
391 nbegin=-ens4_spread/2
392 nend=ens4_spread/2
393 endif
394 do nn=nbegin,nend,1
395 jss=max(j+nn,jds+0)
396 jss=min(jss,jde-1)
397 do n=nbegin,nend,1
398 ens4n=ens4n+1
399 DO K=kts,ktf
400 DO I=ITS,ITF
401 iss=max(i+n,ids+0)
402 iss=min(iss,ide-1)
403 omeg(I,K,ens4n)= -g*rho(i,k,j)*w(iss,k,jss)
404 ! omeg(I,K,ens4n)= -g*rho(i,k,j)*w(i,k,j)
405 Tx(I,K,ens4n)=t2d(i,k)+RTHFTEN(iss,k,jss)*dt
406 ! Tx(I,K,ens4n)=t2d(i,k)+RTHFTEN(i,k,j)*dt
407 if(high_resolution.eq.1)Tx(I,K,ens4n)=t2d(i,k)+ave_f_t(iss,k,jss)*dt
408 IF(Tx(I,K,ens4n).LT.200.)Tx(I,K,ens4n)=T2d(I,K)
409 Qx(I,K,ens4n)=q2d(i,k)+RQVFTEN(iss,k,jss)*dt
410 Qx(I,K,ens4n)=q2d(i,k)+RQVFTEN(i,k,j)*dt
411 if(high_resolution.eq.1)qx(I,K,ens4n)=q2d(i,k)+ave_f_q(iss,k,jss)*dt
412 IF(Qx(I,K,ens4n).LT.1.E-08)Qx(I,K,ens4n)=1.E-08
413 enddo
414 enddo
415 enddo !n
416 enddo !nn
417 do k= kts+1,ktf-1
418 DO I = its,itf
419 if((p2d(i,1)-p2d(i,k)).gt.150.and.p2d(i,k).gt.300)then
420 dp=-.5*(p2d(i,k+1)-p2d(i,k-1))
421 umean(i)=umean(i)+us(i,k)*dp
422 vmean(i)=vmean(i)+vs(i,k)*dp
423 pmean(i)=pmean(i)+dp
424 endif
425 enddo
426 enddo
427 DO I = its,itf
428 umean(i)=umean(i)/pmean(i)
429 vmean(i)=vmean(i)/pmean(i)
430 direction(i)=(atan2(umean(i),vmean(i))+3.1415926)*57.29578
431 if(direction(i).gt.360.)direction(i)=direction(i)-360.
432 ENDDO
433 do n=1,ens4
434 DO K=kts,ktf-1
435 DO I = its,itf
436 dq=(q2d(i,k+1)-q2d(i,k))
437 mconv(i,n)=mconv(i,n)+omeg(i,k,n)*dq/g
438 enddo
439 ENDDO
440 ENDDO
441 do n=1,ens4
442 DO I = its,itf
443 if(mconv(i,n).lt.0.)mconv(i,n)=0.
444 ENDDO
445 ENDDO
446 !
447 !---- CALL CUMULUS PARAMETERIZATION
448 !
449 #if ( WRF_DFI_RADAR == 1 )
450 if(do_capsuppress == 1 ) then
451 DO I= its,itf
452 cap_suppress_j(i)=cap_suppress_loc(i,j)
453 ENDDO
454 endif
455 #endif
456 CALL CUP_enss_3d(outqc,j,AAEQ,T2d,Q2d,TER11,subm,TN,QO,PO,PRET,&
457 P2d,OUTT,OUTQ,DT,itimestep,tkm,PSUR,US,VS,tcrit,iens,tx,qx, &
458 tshall,qshall,kpbli,DHDT,outts,outqs,tscl_kf, &
459 k22s,kbcons,ktops,xmbs, &
460 mconv,massflni,iact_old_gr,omeg,direction,MASSi_FLX, &
461 maxiens,maxens,maxens2,maxens3,ensdim, &
462 APRi_GR,APRi_W,APRi_MC,APRi_ST,APRi_AS, &
463 APRi_CAPMA,APRi_CAPME,APRi_CAPMI,kbcon,ktop,cupclw, &
464 xfi_ens,pri_ens,XLANDi,gswi,edti_out,subt,subq, &
465 ! ruc lv_p,rv_p,cpd_p,g0_p,ichoice,ipr,jpr, &
466 xlv,r_v,cp,g,ichoice,ipr,jpr,ens4,high_resolution, &
467 ishallow_g3,itf,jtf,ktf, &
468 its,ite, jts,jte, kts,kte &
469 #if ( WRF_DFI_RADAR == 1 )
470 ,do_capsuppress,cap_suppress_j &
471 #endif
472 )
475 if(j.lt.jbegc.or.j.gt.jendc)go to 100
476 DO I=ibegc,iendc
477 xmb_shallow(i,j)=xmbs(i)
478 k22_shallow(i,j)=k22s(i)
479 kbcon_shallow(i,j)=kbcons(i)
480 ktop_shallow(i,j)=ktops(i)
481 cuten(i)=0.
482 if(pret(i).gt.0.)then
483 cuten(i)=1.
484 ! raincv(i,j)=pret(i)*dt
485 endif
486 ENDDO
487 DO I=ibegc,iendc
488 DO K=kts,ktf
489 cugd_ttens(I,K,J)=subt(i,k)*cuten(i)*sub_spread
490 cugd_qvtens(I,K,J)=subq(i,k)*cuten(i)*sub_spread
491 ! cugd_tten(I,K,J)=outt(i,k)*cuten(i)
492 ! cugd_qvten(I,K,J)=outq(i,k)*cuten(i)
493 cugd_tten(I,K,J)=outts(i,k)+outt(i,k)*cuten(i)
494 cugd_qvten(I,K,J)=outqs(i,k)+outq(i,k)*cuten(i)
495 cugd_qcten(I,K,J)=outqc(i,k)*cuten(i)
496 ENDDO
497 ENDDO
498 DO I=ibegc,iendc
499 if(pret(i).gt.0.)then
500 raincv(i,j)=pret(i)*dt
501 pratec(i,j)=pret(i)
502 rkbcon = kte+kts - kbcon(i)
503 rktop = kte+kts - ktop(i)
504 if (ktop(i) > HTOP(i,j)) HTOP(i,j) = ktop(i)+.001
505 if (kbcon(i) < HBOT(i,j)) HBOT(i,j) = kbcon(i)+.001
506 endif
507 ENDDO
508 DO n= 1,ensdim
509 DO I= ibegc,iendc
510 xf_ens(i,j,n)=xfi_ens(i,j,n)
511 pr_ens(i,j,n)=pri_ens(i,j,n)
512 ENDDO
513 ENDDO
514 DO I= ibegc,iendc
515 APR_GR(i,j)=apri_gr(i,j)
516 APR_w(i,j)=apri_w(i,j)
517 APR_mc(i,j)=apri_mc(i,j)
518 APR_st(i,j)=apri_st(i,j)
519 APR_as(i,j)=apri_as(i,j)
520 APR_capma(i,j)=apri_capma(i,j)
521 APR_capme(i,j)=apri_capme(i,j)
522 APR_capmi(i,j)=apri_capmi(i,j)
523 mass_flux(i,j)=massi_flx(i,j)
524 edt_out(i,j)=edti_out(i,j)
525 ENDDO
526 IF(PRESENT(RQCCUTEN)) THEN
527 IF ( F_QC ) THEN
528 DO K=kts,ktf
529 DO I=ibegc,iendc
530 RQCCUTEN(I,K,J)=outqc(I,K)*cuten(i)
531 IF ( PRESENT( GDC ) ) GDC(I,K,J)=CUPCLW(I,K)*cuten(i)
532 IF ( PRESENT( GDC2 ) ) GDC2(I,K,J)=0.
533 ENDDO
534 ENDDO
535 ENDIF
536 ENDIF
538 !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
540 IF(PRESENT(RQICUTEN).AND.PRESENT(RQCCUTEN))THEN
541 IF (F_QI) THEN
542 DO K=kts,ktf
543 DO I=ibegc,iendc
544 if(t2d(i,k).lt.258.)then
545 RQICUTEN(I,K,J)=outqc(I,K)*cuten(i)
546 cugd_qcten(i,k,j)=0.
547 RQCCUTEN(I,K,J)=0.
548 IF ( PRESENT( GDC2 ) ) GDC2(I,K,J)=CUPCLW(I,K)*cuten(i)
549 else
550 RQICUTEN(I,K,J)=0.
551 RQCCUTEN(I,K,J)=outqc(I,K)*cuten(i)
552 IF ( PRESENT( GDC ) ) GDC(I,K,J)=CUPCLW(I,K)*cuten(i)
553 endif
554 ENDDO
555 ENDDO
556 ENDIF
557 ENDIF
559 100 continue
561 END SUBROUTINE G3DRV
563 SUBROUTINE CUP_enss_3d(OUTQC,J,AAEQ,T,Q,Z1,sub_mas, &
564 TN,QO,PO,PRE,P,OUTT,OUTQ,DTIME,ktau,tkmax,PSUR,US,VS, &
565 TCRIT,iens,tx,qx, &
566 tshall,qshall,kpbl,dhdt,outts,outqs,tscl_kf, &
567 k23,kbcon3,ktop3,xmb3, &
568 mconv,massfln,iact, &
569 omeg,direction,massflx,maxiens, &
570 maxens,maxens2,maxens3,ensdim, &
571 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
572 APR_CAPMA,APR_CAPME,APR_CAPMI,kbcon,ktop,cupclw, & !-lxz
573 xf_ens,pr_ens,xland,gsw,edt_out,subt,subq, &
574 xl,rv,cp,g,ichoice,ipr,jpr,ens4,high_resolution, &
575 ishallow_g3,itf,jtf,ktf, &
576 its,ite, jts,jte, kts,kte &
577 #if ( WRF_DFI_RADAR == 1 )
578 ! Optional CAP suppress option
579 ,do_capsuppress,cap_suppress_j &
580 #endif
581 )
583 IMPLICIT NONE
585 integer &
586 ,intent (in ) :: &
587 itf,jtf,ktf,ktau, &
588 its,ite, jts,jte, kts,kte,ipr,jpr,ens4,high_resolution
589 integer, intent (in ) :: &
590 j,ensdim,maxiens,ishallow_g3,maxens,maxens2,maxens3,ichoice,iens
591 !
592 !
593 !
594 real, dimension (its:ite,jts:jte,1:ensdim) &
595 ,intent (inout) :: &
596 massfln,xf_ens,pr_ens
597 real, dimension (its:ite,jts:jte) &
598 ,intent (inout ) :: &
599 APR_GR,APR_W,APR_MC,APR_ST,APR_AS,APR_CAPMA, &
600 APR_CAPME,APR_CAPMI,massflx,edt_out
601 real, dimension (its:ite,jts:jte) &
602 ,intent (in ) :: &
603 gsw
604 integer, dimension (its:ite,jts:jte) &
605 ,intent (in ) :: &
606 iact
607 ! outtem = output temp tendency (per s)
608 ! outq = output q tendency (per s)
609 ! outqc = output qc tendency (per s)
610 ! pre = output precip
611 real, dimension (its:ite,kts:kte) &
612 ,intent (inout ) :: &
613 DHDT,OUTT,OUTQ,OUTQC,subt,subq,sub_mas,cupclw,outts,outqs
614 real, dimension (its:ite) &
615 ,intent (out ) :: &
616 pre,xmb3
617 integer, dimension (its:ite) &
618 ,intent (out ) :: &
619 kbcon,ktop,k23,kbcon3,ktop3
620 integer, dimension (its:ite) &
621 ,intent (in ) :: &
622 kpbl
623 !
624 ! basic environmental input includes moisture convergence (mconv)
625 ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
626 ! convection for this call only and at that particular gridpoint
627 !
628 real, dimension (its:ite,kts:kte) &
629 ,intent (in ) :: &
630 T,PO,P,US,VS,tn,tshall,qshall
631 real, dimension (its:ite,kts:kte,1:ens4) &
632 ,intent (inout ) :: &
633 omeg,tx,qx
634 real, dimension (its:ite,kts:kte) &
635 ,intent (inout) :: &
636 Q,QO
637 real, dimension (its:ite) &
638 ,intent (in ) :: &
639 Z1,PSUR,AAEQ,direction,tkmax,xland
640 real, dimension (its:ite,1:ens4) &
641 ,intent (in ) :: &
642 mconv
645 real &
646 ,intent (in ) :: &
647 dtime,tcrit,xl,cp,rv,g,tscl_kf
649 #if ( WRF_DFI_RADAR == 1 )
650 !
651 ! option of cap suppress:
652 ! do_capsuppress = 1 do
653 ! do_capsuppress = other don't
654 !
655 !
656 INTEGER, INTENT(IN ) ,OPTIONAL :: do_capsuppress
657 REAL, DIMENSION( its:ite ),INTENT(IN ) ,OPTIONAL :: cap_suppress_j
658 #endif
660 !
661 ! local ensemble dependent variables in this routine
662 !
663 real, dimension (its:ite,1:maxens) :: &
664 xaa0_ens
665 real, dimension (1:maxens) :: &
666 mbdt_ens
667 real, dimension (1:maxens2) :: &
668 edt_ens
669 real, dimension (its:ite,1:maxens2) :: &
670 edtc
671 real, dimension (its:ite,kts:kte,1:maxens2) :: &
672 dellat_ens,dellaqc_ens,dellaq_ens,pwo_ens,subt_ens,subq_ens
673 !
674 !
675 !
676 !***************** the following are your basic environmental
677 ! variables. They carry a "_cup" if they are
678 ! on model cloud levels (staggered). They carry
679 ! an "o"-ending (z becomes zo), if they are the forced
680 ! variables. They are preceded by x (z becomes xz)
681 ! to indicate modification by some typ of cloud
682 !
683 ! z = heights of model levels
684 ! q = environmental mixing ratio
685 ! qes = environmental saturation mixing ratio
686 ! t = environmental temp
687 ! p = environmental pressure
688 ! he = environmental moist static energy
689 ! hes = environmental saturation moist static energy
690 ! z_cup = heights of model cloud levels
691 ! q_cup = environmental q on model cloud levels
692 ! qes_cup = saturation q on model cloud levels
693 ! t_cup = temperature (Kelvin) on model cloud levels
694 ! p_cup = environmental pressure
695 ! he_cup = moist static energy on model cloud levels
696 ! hes_cup = saturation moist static energy on model cloud levels
697 ! gamma_cup = gamma on model cloud levels
698 !
699 !
700 ! hcd = moist static energy in downdraft
701 ! zd normalized downdraft mass flux
702 ! dby = buoancy term
703 ! entr = entrainment rate
704 ! zd = downdraft normalized mass flux
705 ! entr= entrainment rate
706 ! hcd = h in model cloud
707 ! bu = buoancy term
708 ! zd = normalized downdraft mass flux
709 ! gamma_cup = gamma on model cloud levels
710 ! mentr_rate = entrainment rate
711 ! qcd = cloud q (including liquid water) after entrainment
712 ! qrch = saturation q in cloud
713 ! pwd = evaporate at that level
714 ! pwev = total normalized integrated evaoprate (I2)
715 ! entr= entrainment rate
716 ! z1 = terrain elevation
717 ! entr = downdraft entrainment rate
718 ! jmin = downdraft originating level
719 ! kdet = level above ground where downdraft start detraining
720 ! psur = surface pressure
721 ! z1 = terrain elevation
722 ! pr_ens = precipitation ensemble
723 ! xf_ens = mass flux ensembles
724 ! massfln = downdraft mass flux ensembles used in next timestep
725 ! omeg = omega from large scale model
726 ! mconv = moisture convergence from large scale model
727 ! zd = downdraft normalized mass flux
728 ! zu = updraft normalized mass flux
729 ! dir = "storm motion"
730 ! mbdt = arbitrary numerical parameter
731 ! dtime = dt over which forcing is applied
732 ! iact_gr_old = flag to tell where convection was active
733 ! kbcon = LFC of parcel from k22
734 ! k22 = updraft originating level
735 ! icoic = flag if only want one closure (usually set to zero!)
736 ! dby = buoancy term
737 ! ktop = cloud top (output)
738 ! xmb = total base mass flux
739 ! hc = cloud moist static energy
740 ! hkb = moist static energy at originating level
741 ! mentr_rate = entrainment rate
742 real, dimension (its:ite,kts:kte) :: &
743 he3,hes3,qes3,z3,zdo3, &
744 qes3_cup,q3_cup,he3_cup,hes3_cup,z3_cup,gamma3_cup,t3_cup, &
745 xhe3,xhes3,xqes3,xz3,xt3,xq3, &
746 xqes3_cup,xq3_cup,xhe3_cup,xhes3_cup,xz3_cup,xgamma3_cup, &
747 xt3_cup, &
748 xdby3,xqc3,xhc3,xqrc3,xzu3, &
749 dby3,qc3,pw3,hc3,qrc3,zu3,cd3,DELLAH3,DELLAQ3, &
750 dsubt3,dsubq3,DELLAT3,DELLAQC3
752 real, dimension (its:ite,kts:kte) :: &
753 he,hes,qes,z, &
754 heo,heso,qeso,zo, &
755 xhe,xhes,xqes,xz,xt,xq, &
757 qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
758 qeso_cup,qo_cup,heo_cup,heso_cup,zo_cup,po_cup,gammao_cup, &
759 tn_cup, &
760 xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup,xp_cup,xgamma_cup, &
761 xt_cup, &
763 dby,qc,qrcd,pwd,pw,hcd,qcd,dbyd,hc,qrc,zu,zd,clw_all, &
764 dbyo,qco,qrcdo,pwdo,pwo,hcdo,qcdo,dbydo,hco,qrco,zuo,zdo, &
765 xdby,xqc,xqrcd,xpwd,xpw,xhcd,xqcd,xhc,xqrc,xzu,xzd, &
767 ! cd = detrainment function for updraft
768 ! cdd = detrainment function for downdraft
769 ! dellat = change of temperature per unit mass flux of cloud ensemble
770 ! dellaq = change of q per unit mass flux of cloud ensemble
771 ! dellaqc = change of qc per unit mass flux of cloud ensemble
773 cd,cdd,scr1,DELLAH,DELLAQ,DELLAT,DELLAQC,dsubt,dsubq
775 ! aa0 cloud work function for downdraft
776 ! edt = epsilon
777 ! aa0 = cloud work function without forcing effects
778 ! aa1 = cloud work function with forcing effects
779 ! xaa0 = cloud work function with cloud effects (ensemble dependent)
780 ! edt = epsilon
782 real, dimension (its:ite) :: &
783 aa3,hkb3,qkb3,pwav3,bu3,xaa3,xhkb3, &
784 edt,edto,edtx,AA1,AA0,XAA0,HKB, &
785 HKBO,aad,XHKB,QKB,QKBO,edt3, &
786 XMB,XPWAV,XPWEV,PWAV,PWEV,PWAVO, &
787 PWEVO,BU,BUO,cap_max,xland1, &
788 cap_max_increment,closure_n,cap_max3
789 real, dimension (its:ite,1:ens4) :: &
790 axx
791 integer, dimension (its:ite) :: &
792 kzdown,KDET,K22,KB,JMIN,kstabi,kstabm,K22x,jmin3,kdet3, & !-lxz
793 KBCONx,KBx,KTOPx,ierr,ierr2,ierr3,KBMAX,ierr5
795 integer :: &
796 nall,iedt,nens,nens3,ki,I,K,KK,iresult
797 real :: &
798 day,dz,mbdt,entr_rate,radius,entrd_rate,mentr_rate,mentrd_rate, &
799 zcutdown,edtmax,edtmin,depth_min,zkbmax,z_detr,zktop, &
800 massfld,dh,cap_maxs,trash,entr_rate3,mentr_rate3
802 integer :: jmini
803 logical :: keep_going
804 real xff_shal(9),blqe,xkshal
808 day=86400.
809 do i=its,itf
810 xmb3(i)=0.
811 closure_n(i)=16.
812 xland1(i)=1.
813 if(xland(i).gt.1.5)xland1(i)=0.
814 ! cap_max_increment(i)=50.
815 cap_max_increment(i)=25.
816 enddo
817 !
818 !--- specify entrainmentrate and detrainmentrate
819 !
820 if(iens.le.4)then
821 radius=14000.-float(iens)*2000.
822 else
823 radius=12000.
824 endif
825 !
826 !--- gross entrainment rate (these may be changed later on in the
827 !--- program, depending what your detrainment is!!)
828 !
829 entr_rate =.2/radius
830 entr_rate3=.2/500.
831 !
832 !--- entrainment of mass
833 !
834 mentrd_rate=0.
835 mentr_rate=entr_rate
836 mentr_rate3=entr_rate3
837 !
838 !--- initial detrainmentrates
839 !
840 do k=kts,ktf
841 do i=its,itf
842 cupclw(i,k)=0.
843 cd(i,k)=0.01*entr_rate
844 cd3(i,k)=entr_rate3
845 cdd(i,k)=0.
846 zdo3(i,k)=0.
847 enddo
848 enddo
849 !
850 !--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
851 ! base mass flux
852 !
853 edtmax=1.
854 edtmin=.2
855 !
856 !--- minimum depth (m), clouds must have
857 !
858 depth_min=500.
859 !
860 !--- maximum depth (mb) of capping
861 !--- inversion (larger cap = no convection)
862 !
863 ! cap_maxs=125.
864 cap_maxs=75.
865 DO i=its,itf
866 kbmax(i)=1
867 jmin3(i)=0
868 kdet3(i)=0
869 aa0(i)=0.
870 aa1(i)=0.
871 aa3(i)=0.
872 aad(i)=0.
873 edt(i)=0.
874 kstabm(i)=ktf-1
875 IERR(i)=0
876 IERR2(i)=0
877 IERR3(i)=0
878 IERR5(i)=0
879 if(aaeq(i).lt.-0.1)then
880 ierr(i)=20
881 endif
882 enddo
883 !
884 !--- first check for upstream convection
885 !
886 #if ( WRF_DFI_RADAR == 1 )
887 if(do_capsuppress == 1) then
888 do i=its,itf
889 cap_max(i)=cap_maxs
890 cap_max3(i)=25.
891 if(gsw(i,j).lt.1.or.high_resolution.eq.1)cap_max(i)=25.
892 if (abs(cap_suppress_j(i) - 1.0 ) < 0.1 ) then
893 cap_max(i)=cap_maxs+75.
894 elseif (abs(cap_suppress_j(i) - 0.0 ) < 0.1 ) then
895 cap_max(i)=10.0
896 endif
897 iresult=0
898 enddo
899 else
900 do i=its,itf
901 cap_max(i)=cap_maxs
902 if(gsw(i,j).lt.1.or.high_resolution.eq.1)cap_max(i)=25.
903 iresult=0
904 enddo
905 endif
907 do i=its,itf
908 edt_out(i,j)=cap_max(i)
909 enddo
910 #else
911 do i=its,itf
912 cap_max(i)=cap_maxs
913 if(gsw(i,j).lt.1.or.high_resolution.eq.1)cap_max(i)=25.
914 iresult=0
916 enddo
917 #endif
918 !
919 !--- max height(m) above ground where updraft air can originate
920 !
921 zkbmax=4000.
922 !
923 !--- height(m) above which no downdrafts are allowed to originate
924 !
925 zcutdown=3000.
926 !
927 !--- depth(m) over which downdraft detrains all its mass
928 !
929 z_detr=1250.
930 !
931 do nens=1,maxens
932 mbdt_ens(nens)=(float(nens)-3.)*dtime*1.e-3+dtime*5.E-03
933 enddo
934 do nens=1,maxens2
935 edt_ens(nens)=.95-float(nens)*.01
936 enddo
937 !
938 !--- environmental conditions, FIRST HEIGHTS
939 !
940 do i=its,itf
941 if(ierr(i).ne.20)then
942 do k=1,maxens*maxens2*maxens3
943 xf_ens(i,j,(iens-1)*maxens*maxens2*maxens3+k)=0.
944 pr_ens(i,j,(iens-1)*maxens*maxens2*maxens3+k)=0.
945 enddo
946 endif
947 enddo
948 !
949 !--- calculate moist static energy, heights, qes
950 !
951 call cup_env(z,qes,he,hes,t,q,p,z1, &
952 psur,ierr,tcrit,0,xl,cp, &
953 itf,jtf,ktf, &
954 its,ite, jts,jte, kts,kte)
955 call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
956 psur,ierr,tcrit,0,xl,cp, &
957 itf,jtf,ktf, &
958 its,ite, jts,jte, kts,kte)
959 !
960 !--- environmental values on cloud levels
961 !
962 call cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup,he_cup, &
963 hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
964 ierr,z1,xl,rv,cp, &
965 itf,jtf,ktf, &
966 its,ite, jts,jte, kts,kte)
967 call cup_env_clev(tn,qeso,qo,heo,heso,zo,po,qeso_cup,qo_cup, &
968 heo_cup,heso_cup,zo_cup,po_cup,gammao_cup,tn_cup,psur, &
969 ierr,z1,xl,rv,cp, &
970 itf,jtf,ktf, &
971 its,ite, jts,jte, kts,kte)
972 do i=its,itf
973 if(ierr(i).eq.0)then
974 !
975 do k=kts,ktf
976 if(zo_cup(i,k).gt.zkbmax+z1(i))then
977 kbmax(i)=k
978 go to 25
979 endif
980 enddo
981 25 continue
982 !
983 !--- level where detrainment for downdraft starts
984 !
985 do k=kts,ktf
986 if(zo_cup(i,k).gt.z_detr+z1(i))then
987 kdet(i)=k
988 go to 26
989 endif
990 enddo
991 26 continue
992 !
993 endif
994 enddo
995 !
996 !
997 !
998 !------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
999 !
1000 CALL cup_MAXIMI(HEO_CUP,3,KBMAX,K22,ierr, &
1001 itf,jtf,ktf, &
1002 its,ite, jts,jte, kts,kte)
1003 DO 36 i=its,itf
1004 IF(ierr(I).eq.0.)THEN
1005 IF(K22(I).GE.KBMAX(i))ierr(i)=2
1006 endif
1007 36 CONTINUE
1008 !
1009 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
1010 !
1011 call cup_kbcon(cap_max_increment,1,k22,kbcon,heo_cup,heso_cup, &
1012 ierr,kbmax,po_cup,cap_max, &
1013 itf,jtf,ktf, &
1014 its,ite, jts,jte, kts,kte)
1015 !
1016 !--- increase detrainment in stable layers
1017 !
1018 CALL cup_minimi(HEso_cup,Kbcon,kstabm,kstabi,ierr, &
1019 itf,jtf,ktf, &
1020 its,ite, jts,jte, kts,kte)
1021 do i=its,itf
1022 IF(ierr(I).eq.0.)THEN
1023 if(kstabm(i)-1.gt.kstabi(i))then
1024 do k=kstabi(i),kstabm(i)-1
1025 cd(i,k)=cd(i,k-1)+.15*entr_rate
1026 if(cd(i,k).gt.1.0*entr_rate)cd(i,k)=1.0*entr_rate
1027 enddo
1028 ENDIF
1029 ENDIF
1030 ENDDO
1031 !
1032 !--- calculate incloud moist static energy
1033 !
1034 call cup_up_he(k22,hkb,z_cup,cd,mentr_rate,he_cup,hc, &
1035 kbcon,ierr,dby,he,hes_cup, &
1036 itf,jtf,ktf, &
1037 its,ite, jts,jte, kts,kte)
1038 call cup_up_he(k22,hkbo,zo_cup,cd,mentr_rate,heo_cup,hco, &
1039 kbcon,ierr,dbyo,heo,heso_cup, &
1040 itf,jtf,ktf, &
1041 its,ite, jts,jte, kts,kte)
1043 !--- DETERMINE CLOUD TOP - KTOP
1044 !
1045 call cup_ktop(1,dbyo,kbcon,ktop,ierr, &
1046 itf,jtf,ktf, &
1047 its,ite, jts,jte, kts,kte)
1048 DO 37 i=its,itf
1049 kzdown(i)=0
1050 if(ierr(i).eq.0)then
1051 zktop=(zo_cup(i,ktop(i))-z1(i))*.6
1052 zktop=min(zktop+z1(i),zcutdown+z1(i))
1053 do k=kts,kte
1054 if(zo_cup(i,k).gt.zktop)then
1055 kzdown(i)=k
1056 go to 37
1057 endif
1058 enddo
1059 endif
1060 37 CONTINUE
1061 !
1062 !--- DOWNDRAFT ORIGINATING LEVEL - JMIN
1063 !
1064 call cup_minimi(HEso_cup,K22,kzdown,JMIN,ierr, &
1065 itf,jtf,ktf, &
1066 its,ite, jts,jte, kts,kte)
1067 DO 100 i=its,ite
1068 IF(ierr(I).eq.0.)THEN
1069 !
1070 !--- check whether it would have buoyancy, if there where
1071 !--- no entrainment/detrainment
1072 !
1073 jmini = jmin(i)
1074 keep_going = .TRUE.
1075 do while ( keep_going )
1076 keep_going = .FALSE.
1077 if ( jmini - 1 .lt. kdet(i) ) kdet(i) = jmini-1
1078 if ( jmini .ge. ktop(i)-1 ) jmini = ktop(i) - 2
1079 ki = jmini
1080 hcdo(i,ki)=heso_cup(i,ki)
1081 DZ=Zo_cup(i,Ki+1)-Zo_cup(i,Ki)
1082 dh=0.
1083 do k=ki-1,1,-1
1084 hcdo(i,k)=heso_cup(i,jmini)
1085 DZ=Zo_cup(i,K+1)-Zo_cup(i,K)
1086 dh=dh+dz*(HCDo(i,K)-heso_cup(i,k))
1087 if(dh.gt.0.)then
1088 jmini=jmini-1
1089 if ( jmini .gt. 3 ) then
1090 keep_going = .TRUE.
1091 else
1092 ierr(i) = 9
1093 exit
1094 endif
1095 endif
1096 enddo
1097 enddo
1098 jmin(i) = jmini
1099 if ( jmini .le. 3 ) then
1100 ierr(i)=4
1101 endif
1102 ENDIF
1103 100 continue
1104 !
1105 ! - Must have at least depth_min m between cloud convective base
1106 ! and cloud top.
1107 !
1108 do i=its,itf
1109 IF(ierr(I).eq.0.)THEN
1110 IF(-zo_cup(I,KBCON(I))+zo_cup(I,KTOP(I)).LT.depth_min)then
1111 ierr(i)=6
1112 endif
1113 endif
1114 enddo
1116 !
1117 !c--- normalized updraft mass flux profile
1118 !
1119 call cup_up_nms(zu,z_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
1120 itf,jtf,ktf, &
1121 its,ite, jts,jte, kts,kte)
1122 call cup_up_nms(zuo,zo_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
1123 itf,jtf,ktf, &
1124 its,ite, jts,jte, kts,kte)
1125 !
1126 !c--- normalized downdraft mass flux profile,also work on bottom detrainment
1127 !--- in this routine
1128 !
1129 call cup_dd_nms(zd,z_cup,cdd,mentrd_rate,jmin,ierr, &
1130 0,kdet,z1, &
1131 itf,jtf,ktf, &
1132 its,ite, jts,jte, kts,kte)
1133 call cup_dd_nms(zdo,zo_cup,cdd,mentrd_rate,jmin,ierr, &
1134 1,kdet,z1, &
1135 itf,jtf,ktf, &
1136 its,ite, jts,jte, kts,kte)
1137 !
1138 !--- downdraft moist static energy
1139 !
1140 call cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,mentrd_rate, &
1141 jmin,ierr,he,dbyd,he_cup, &
1142 itf,jtf,ktf, &
1143 its,ite, jts,jte, kts,kte)
1144 call cup_dd_he(heso_cup,zdo,hcdo,zo_cup,cdd,mentrd_rate, &
1145 jmin,ierr,heo,dbydo,he_cup,&
1146 itf,jtf,ktf, &
1147 its,ite, jts,jte, kts,kte)
1148 !
1149 !--- calculate moisture properties of downdraft
1150 !
1151 call cup_dd_moisture_3d(zd,hcd,hes_cup,qcd,qes_cup, &
1152 pwd,q_cup,z_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
1153 pwev,bu,qrcd,q,he,t_cup,2,xl,high_resolution, &
1154 itf,jtf,ktf, &
1155 its,ite, jts,jte, kts,kte)
1156 call cup_dd_moisture_3d(zdo,hcdo,heso_cup,qcdo,qeso_cup, &
1157 pwdo,qo_cup,zo_cup,cdd,mentrd_rate,jmin,ierr,gammao_cup, &
1158 pwevo,bu,qrcdo,qo,heo,tn_cup,1,xl,high_resolution, &
1159 itf,jtf,ktf, &
1160 its,ite, jts,jte, kts,kte)
1161 !
1162 !--- calculate moisture properties of updraft
1163 !
1164 call cup_up_moisture('deep',ierr,z_cup,qc,qrc,pw,pwav, &
1165 kbcon,ktop,cd,dby,mentr_rate,clw_all, &
1166 q,GAMMA_cup,zu,qes_cup,k22,q_cup,xl, &
1167 itf,jtf,ktf, &
1168 its,ite, jts,jte, kts,kte)
1169 do k=kts,ktf
1170 do i=its,itf
1171 cupclw(i,k)=qrc(i,k)
1172 enddo
1173 enddo
1174 call cup_up_moisture('deep',ierr,zo_cup,qco,qrco,pwo,pwavo, &
1175 kbcon,ktop,cd,dbyo,mentr_rate,clw_all, &
1176 qo,GAMMAo_cup,zuo,qeso_cup,k22,qo_cup,xl,&
1177 itf,jtf,ktf, &
1178 its,ite, jts,jte, kts,kte)
1179 !
1180 !--- calculate workfunctions for updrafts
1181 !
1182 call cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
1183 kbcon,ktop,ierr, &
1184 itf,jtf,ktf, &
1185 its,ite, jts,jte, kts,kte)
1186 call cup_up_aa0(aa1,zo,zuo,dbyo,GAMMAo_CUP,tn_cup, &
1187 kbcon,ktop,ierr, &
1188 itf,jtf,ktf, &
1189 its,ite, jts,jte, kts,kte)
1190 do i=its,itf
1191 if(ierr(i).eq.0)then
1192 if(aa1(i).eq.0.)then
1193 ierr(i)=17
1194 endif
1195 endif
1196 enddo
1197 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1198 ! NEXT section for shallow convection
1199 !
1200 if(ishallow_g3.eq.1)then
1201 ! write(0,*)'now do shallow for j = ',j
1202 call cup_env(z3,qes3,he3,hes3,tshall,qshall,po,z1, &
1203 psur,ierr5,tcrit,0,xl,cp, &
1204 itf,jtf,ktf, &
1205 its,ite, jts,jte, kts,kte)
1206 call cup_env_clev(tshall,qes3,qshall,he3,hes3,z3,po,qes3_cup,q3_cup, &
1207 he3_cup,hes3_cup,z3_cup,po_cup,gamma3_cup,t3_cup,psur, &
1208 ierr5,z1,xl,rv,cp, &
1209 itf,jtf,ktf, &
1210 its,ite, jts,jte, kts,kte)
1211 CALL cup_MAXIMI(HE3_CUP,1,kbmax,K23,ierr5, &
1212 itf,jtf,ktf, &
1213 its,ite, jts,jte, kts,kte)
1214 DO i=its,itf
1215 if(kpbl(i).gt.1)cap_max3(i)=po_cup(i,kpbl(i))
1216 IF(ierr5(I).eq.0.)THEN
1217 IF(K23(I).Gt.Kbmax(i))ierr5(i)=2
1218 ! if(kpbl(i).gt.5)k23(i)=kpbl(i)
1219 endif
1220 ENDDO
1221 call cup_kbcon(cap_max_increment,5,k23,kbcon3,he3_cup,hes3_cup, &
1222 ierr5,kbmax,po_cup,cap_max3, &
1223 ! ierr5,kpbl,po_cup,cap_max3, &
1224 itf,jtf,ktf, &
1225 its,ite, jts,jte, kts,kte)
1226 call cup_up_he(k23,hkb3,z3_cup,cd3,mentr_rate3,he3_cup,hc3, &
1227 kbcon3,ierr5,dby3,he3,hes3_cup, &
1228 itf,jtf,ktf, &
1229 its,ite, jts,jte, kts,kte)
1230 call cup_ktop(1,dby3,kbcon3,ktop3,ierr5, &
1231 itf,jtf,ktf, &
1232 its,ite, jts,jte, kts,kte)
1233 call cup_up_nms(zu3,z3_cup,mentr_rate3,cd3,kbcon3,ktop3, &
1234 ierr5,k23, &
1235 itf,jtf,ktf, &
1236 its,ite, jts,jte, kts,kte)
1237 call cup_up_moisture('shallow',ierr5,z3_cup,qc3,qrc3,pw3,pwav3, &
1238 kbcon3,ktop3,cd3,dby3,mentr_rate3,clw_all, &
1239 qshall,GAMMA3_cup,zu3,qes3_cup,k23,q3_cup,xl,&
1240 itf,jtf,ktf, &
1241 its,ite, jts,jte, kts,kte)
1242 call cup_up_aa0(aa3,z3,zu3,dby3,GAMMA3_CUP,t3_cup, &
1243 kbcon3,ktop3,ierr5, &
1244 itf,jtf,ktf, &
1245 its,ite, jts,jte, kts,kte)
1246 do i=its,itf
1247 if(ierr5(i).eq.0)then
1248 if(aa3(i).eq.0.)then
1249 ierr5(i)=17
1250 endif
1251 endif
1252 enddo
1253 ! call cup_dellabot('shallow',ipr,jpr,q3_cup,ierr5,z3_cup,po,qrcdo,edto, &
1254 ! zdo,cdd,q3,dellaq3,dsubq,j,mentrd_rate,z3,g,&
1255 ! itf,jtf,ktf, &
1256 ! its,ite, jts,jte, kts,kte)
1257 call cup_dellas_3d(ierr5,z3_cup,po_cup,hcdo,edt3,zdo3,cdd, &
1258 he3,dellah3,dsubt3,j,mentrd_rate,zu3,g, &
1259 cd3,hc3,ktop3,k23,kbcon3,mentr_rate3,jmin,he3_cup,kdet, &
1260 k23,ipr,jpr,'shallow',0, &
1261 itf,jtf,ktf, &
1262 its,ite, jts,jte, kts,kte)
1263 call cup_dellas_3d(ierr5,z3_cup,po_cup,qrcdo,edt3,zdo3,cdd, &
1264 qshall,dellaq3,dsubq3,j,mentrd_rate,zu3,g, &
1265 cd3,qc3,ktop3,k23,kbcon3,mentr_rate3,jmin,q3_cup,kdet, &
1266 k23,ipr,jpr,'shallow',0, &
1267 itf,jtf,ktf, &
1268 its,ite, jts,jte, kts,kte )
1269 mbdt=mbdt_ens(1)
1270 do k=kts,ktf
1271 do i=its,itf
1272 dellat3(i,k)=0.
1273 if(ierr5(i).eq.0)then
1274 trash=dsubt3(i,k)
1275 XHE3(I,K)=(dsubt3(i,k)+DELLAH3(I,K))*MBDT+HE3(I,K)
1276 XQ3(I,K)=(dsubq3(i,k)+DELLAQ3(I,K))*MBDT+QSHALL(I,K)
1277 DELLAT3(I,K)=(1./cp)*(DELLAH3(I,K)-xl*DELLAQ3(I,K))
1278 dSUBT3(I,K)=(1./cp)*(dsubt3(i,k)-xl*dsubq3(i,k))
1279 XT3(I,K)= (DELLAT3(I,K)+dsubt3(i,k))*MBDT+TSHALL(I,K)
1280 IF(XQ3(I,K).LE.0.)XQ3(I,K)=1.E-08
1281 if(i.eq.ipr.and.j.eq.jpr)then
1282 write(0,*)k,trash,DELLAQ3(I,K),dsubq3(I,K),dsubt3(i,k)
1283 endif
1284 ENDIF
1285 enddo
1286 enddo
1287 do i=its,itf
1288 if(ierr5(i).eq.0)then
1289 XHE3(I,ktf)=HE3(I,ktf)
1290 XQ3(I,ktf)=QSHALL(I,ktf)
1291 XT3(I,ktf)=TSHALL(I,ktf)
1292 IF(XQ3(I,ktf).LE.0.)XQ3(I,ktf)=1.E-08
1293 endif
1294 enddo
1295 !
1296 !--- calculate moist static energy, heights, qes
1297 !
1298 call cup_env(xz3,xqes3,xhe3,xhes3,xt3,xq3,po,z1, &
1299 psur,ierr5,tcrit,2,xl,cp, &
1300 itf,jtf,ktf, &
1301 its,ite, jts,jte, kts,kte)
1302 !
1303 !--- environmental values on cloud levels
1304 !
1305 call cup_env_clev(xt3,xqes3,xq3,xhe3,xhes3,xz3,po,xqes3_cup,xq3_cup, &
1306 xhe3_cup,xhes3_cup,xz3_cup,po_cup,gamma3_cup,xt3_cup,psur, &
1307 ierr5,z1,xl,rv,cp, &
1308 itf,jtf,ktf, &
1309 its,ite, jts,jte, kts,kte)
1310 !
1311 !
1312 !**************************** static control
1313 !
1314 !--- moist static energy inside cloud
1315 !
1316 do i=its,itf
1317 if(ierr5(i).eq.0)then
1318 xhkb3(i)=xhe3(i,k23(i))
1319 endif
1320 enddo
1321 call cup_up_he(k23,xhkb3,xz3_cup,cd3,mentr_rate3,xhe3_cup,xhc3, &
1322 kbcon3,ierr5,xdby3,xhe3,xhes3_cup, &
1323 itf,jtf,ktf, &
1324 its,ite, jts,jte, kts,kte)
1325 !
1326 !c--- normalized mass flux profile and CWF
1327 !
1328 call cup_up_nms(xzu3,xz3_cup,mentr_rate3,cd3,kbcon3,ktop3,ierr5,k23, &
1329 itf,jtf,ktf, &
1330 its,ite, jts,jte, kts,kte)
1331 call cup_up_aa0(xaa3,xz3,xzu3,xdby3,GAMMA3_CUP,xt3_cup, &
1332 kbcon3,ktop3,ierr5, &
1333 itf,jtf,ktf, &
1334 its,ite, jts,jte, kts,kte)
1335 !
1336 ! now for shallow forcing
1337 !
1338 do i=its,itf
1339 xmb3(i)=0.
1340 xff_shal(1:9)=0.
1341 if(ierr5(i).eq.0)then
1342 xkshal=(xaa3(i)-aa0(i))/mbdt
1343 if(xkshal.le.0.and.xkshal.gt.-1.e-6)xkshal=-1.e-6
1344 if(xkshal.gt.0.and.xkshal.lt.+1.e-6)xkshal=+1.e-6
1345 xff_shal(1)=max(0.,-(aa3(i)-aa0(i))/(xkshal*dtime))
1346 xff_shal(2)=max(0.,-(aa3(i)-aa0(i))/(xkshal*dtime))
1347 xff_shal(3)=max(0.,-(aa3(i)-aa0(i))/(xkshal*dtime))
1348 xff_shal(4)=max(0.,-aa3(i)/(xkshal*tscl_KF))
1349 xff_shal(5)=max(0.,-aa3(i)/(xkshal*tscl_KF))
1350 xff_shal(6)=max(0.,-aa3(i)/(xkshal*tscl_KF))
1351 ! boundary layer QE (from Saulo Freitas)
1352 blqe=0.
1353 if(k23(i).lt.kpbl(i)+1)then
1354 do k=1,kbcon3(i)-1
1355 blqe=blqe+100.*dhdt(i,k)*(p_cup(i,k)-p_cup(i,k+1))/g
1356 enddo
1357 trash=max((hc3(i,kbcon3(i))-he_cup(i,kbcon3(i))),1.e-4)
1358 xff_shal(7)=max(0.,blqe/trash)
1359 xff_shal(7)=min(0.5,xff_shal(7))
1360 else
1361 xff_shal(7)=0.
1362 endif
1363 xff_shal(8)= xff_shal(7)
1364 xff_shal(9)= xff_shal(7)
1365 do k=1,9
1366 xmb3(i)=xmb3(i)+xff_shal(k)
1367 enddo
1368 xmb3(i)=min(.1,xmb3(i)/9.)
1369 if(xmb3(i).eq.0.)ierr5(i)=22
1370 if(xmb3(i).lt.0.)then
1371 ierr5(i)=21
1372 ! write(0,*)'neg xmb,i,j,xmb3 for shallow = ',i,j,k23(i),ktop3(i),kbcon3(i),kpbl(i)
1373 endif
1374 endif
1375 ! if(ierr5(i).eq.0)write(0,*)'i,j,xmb3 for shallow = ',i,j,xmb3(i),k23(i),ktop3(i)
1376 ! if(ierr5(i).eq.0.and.i.eq.12.and.j.eq.25)write(0,*)'i,j,xmb3 for shallow = ',k23(i),ktop3(i),kbcon3(i),kpbl(i)
1377 ! if(ierr5(i).eq.0)write(0,*)'i,j,xmb3 for shallow = ',i,j,k23(i),ktop3(i),kbcon3(i),kpbl(i)
1378 if(ierr5(i).eq.0)then
1379 !
1380 ! got the mass flux, now do feedback
1381 trash=0.
1382 !
1383 do k=2,ktop3(i)
1384 trash=max(trash,86400.*(dsubt3(i,k)+dellat3(i,k))*xmb3(i))
1385 enddo
1386 if(trash.gt.150.)xmb3(i)=xmb3(i)*150./trash
1387 do k=2,ktop3(i)
1388 outts(i,k)=(dsubt3(i,k)+dellat3(i,k))*xmb3(i)
1389 outqs(i,k)=(dsubq3(i,k)+dellaq3(i,k))*xmb3(i)
1390 enddo
1391 endif
1392 enddo
1393 if(j.eq.-25)then
1394 write(0,*)'!!!!!!!! j = ',j,' !!!!!!!!!!!!!!!!!!!!'
1395 i=12
1396 write(0,*)k23(i),kbcon3(i),ktop3(i)
1397 write(0,*)kpbl(i),ierr5(i),ierr(i)
1398 write(0,*)xmb3(i),xff_shal(1:9)
1399 write(0,*)xaa3(i),aa1(i),aa0(i),aa3(i)
1400 do k=1,ktf
1401 write(0,*)po(i,k),he3(i,k),hes3(i,k),dellah3(i,k)
1402 enddo
1403 do k=1,ktf
1404 write(0,*)zu3(i,k),hc3(i,k),dsubt3(i,k),dellat3(i,k)
1405 enddo
1406 do k=1,ktop3(i)+1
1407 blqe=cp*outts(i,k)+xl*outqs(i,k)
1408 write(0,*)outts(i,k),outqs(i,k),blqe
1409 enddo
1410 endif
1411 !
1412 ! done shallow
1413 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1414 ENDIF
1415 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1418 call cup_axx(tcrit,kbmax,z1,p,psur,xl,rv,cp,tx,qx,axx,ierr, &
1419 cap_max,cap_max_increment,entr_rate,mentr_rate,&
1420 j,itf,jtf,ktf, &
1421 its,ite, jts,jte, kts,kte,ens4)
1423 !
1424 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
1425 !
1426 call cup_dd_edt(ierr,us,vs,zo,ktop,kbcon,edt,po,pwavo, &
1427 pwevo,edtmax,edtmin,maxens2,edtc, &
1428 itf,jtf,ktf, &
1429 its,ite, jts,jte, kts,kte)
1430 do 250 iedt=1,maxens2
1431 do i=its,itf
1432 if(ierr(i).eq.0)then
1433 edt(i)=edtc(i,iedt)
1434 edto(i)=edtc(i,iedt)
1435 edtx(i)=edtc(i,iedt)
1436 edt_out(i,j)=edtc(i,2)
1437 if(high_resolution.eq.1)then
1438 edt(i)=edtc(i,3)
1439 edto(i)=edtc(i,3)
1440 edtx(i)=edtc(i,3)
1441 edt_out(i,j)=edtc(i,3)
1442 endif
1443 endif
1444 enddo
1445 do k=kts,ktf
1446 do i=its,itf
1447 subt_ens(i,k,iedt)=0.
1448 subq_ens(i,k,iedt)=0.
1449 dellat_ens(i,k,iedt)=0.
1450 dellaq_ens(i,k,iedt)=0.
1451 dellaqc_ens(i,k,iedt)=0.
1452 pwo_ens(i,k,iedt)=0.
1453 enddo
1454 enddo
1455 !
1456 if(j.eq.jpr.and.iedt.eq.1.and.ipr.gt.its.and.ipr.lt.ite)then
1457 ! if(j.eq.jpr)then
1458 i=ipr
1459 ! write(0,*)'in 250 loop ',iedt,edt(ipr),ierr(ipr)
1460 ! if(ierr(i).eq.0.or.ierr(i).eq.3)then
1461 write(0,*)'250',k22(I),kbcon(i),ktop(i),jmin(i)
1462 write(0,*)edt(i),aa0(i),aa1(i)
1463 do k=kts,ktf
1464 write(0,*)k,z(i,k),he(i,k),hes(i,k)
1465 enddo
1466 write(0,*)'end 250 loop ',iedt,edt(ipr),ierr(ipr)
1467 do k=1,ktop(i)+1
1468 write(0,*)zu(i,k),zd(i,k),pw(i,k),pwd(i,k)
1469 enddo
1470 ! endif
1471 endif
1472 do i=its,itf
1473 aad(i)=0.
1474 enddo
1475 !
1476 !--- change per unit mass that a model cloud would modify the environment
1477 !
1478 !--- 1. in bottom layer
1479 !
1480 call cup_dellabot('deep',ipr,jpr,heo_cup,ierr,zo_cup,po,hcdo,edto, &
1481 zdo,cdd,heo,dellah,dsubt,j,mentrd_rate,zo,g, &
1482 itf,jtf,ktf, &
1483 its,ite, jts,jte, kts,kte)
1484 call cup_dellabot('deep',ipr,jpr,qo_cup,ierr,zo_cup,po,qrcdo,edto, &
1485 zdo,cdd,qo,dellaq,dsubq,j,mentrd_rate,zo,g,&
1486 itf,jtf,ktf, &
1487 its,ite, jts,jte, kts,kte)
1488 !
1489 !--- 2. everywhere else
1490 !
1491 call cup_dellas_3d(ierr,zo_cup,po_cup,hcdo,edto,zdo,cdd, &
1492 heo,dellah,dsubt,j,mentrd_rate,zuo,g, &
1493 cd,hco,ktop,k22,kbcon,mentr_rate,jmin,heo_cup,kdet, &
1494 k22,ipr,jpr,'deep',high_resolution, &
1495 itf,jtf,ktf, &
1496 its,ite, jts,jte, kts,kte)
1497 !
1498 !-- take out cloud liquid water for detrainment
1499 !
1500 !?? do k=kts,ktf
1501 do k=kts,ktf-1
1502 do i=its,itf
1503 scr1(i,k)=0.
1504 dellaqc(i,k)=0.
1505 if(ierr(i).eq.0)then
1506 scr1(i,k)=qco(i,k)-qrco(i,k)
1507 if(k.eq.ktop(i)-0)dellaqc(i,k)= &
1508 .01*zuo(i,ktop(i))*qrco(i,ktop(i))* &
1509 9.81/(po_cup(i,k)-po_cup(i,k+1))
1510 if(k.lt.ktop(i).and.k.gt.kbcon(i))then
1511 dz=zo_cup(i,k+1)-zo_cup(i,k)
1512 dellaqc(i,k)=.01*9.81*cd(i,k)*dz*zuo(i,k) &
1513 *.5*(qrco(i,k)+qrco(i,k+1))/ &
1514 (po_cup(i,k)-po_cup(i,k+1))
1515 endif
1516 endif
1517 enddo
1518 enddo
1519 call cup_dellas_3d(ierr,zo_cup,po_cup,qrcdo,edto,zdo,cdd, &
1520 qo,dellaq,dsubq,j,mentrd_rate,zuo,g, &
1521 cd,qco,ktop,k22,kbcon,mentr_rate,jmin,qo_cup,kdet, &
1522 k22,ipr,jpr,'deep',high_resolution, &
1523 itf,jtf,ktf, &
1524 its,ite, jts,jte, kts,kte )
1525 !
1526 !--- using dellas, calculate changed environmental profiles
1527 !
1528 ! do 200 nens=1,maxens
1529 mbdt=mbdt_ens(2)
1530 do i=its,itf
1531 xaa0_ens(i,1)=0.
1532 xaa0_ens(i,2)=0.
1533 xaa0_ens(i,3)=0.
1534 enddo
1536 if(j.eq.jpr)then
1537 write(0,*)'xt',xl,'DELLAH(I,K),DELLAQ(I,K),dsubq(I,K),dsubt(i,k)'
1538 endif
1539 do k=kts,ktf
1540 do i=its,itf
1541 dellat(i,k)=0.
1542 if(ierr(i).eq.0)then
1543 trash=dsubt(i,k)
1544 XHE(I,K)=(dsubt(i,k)+DELLAH(I,K))*MBDT+HEO(I,K)
1545 XQ(I,K)=(dsubq(i,k)+DELLAQ(I,K))*MBDT+QO(I,K)
1546 DELLAT(I,K)=(1./cp)*(DELLAH(I,K)-xl*DELLAQ(I,K))
1547 dSUBT(I,K)=(1./cp)*(dsubt(i,k)-xl*dsubq(i,k))
1548 XT(I,K)= (DELLAT(I,K)+dsubt(i,k))*MBDT+TN(I,K)
1549 IF(XQ(I,K).LE.0.)XQ(I,K)=1.E-08
1550 if(i.eq.ipr.and.j.eq.jpr)then
1551 write(0,*)k,trash,DELLAQ(I,K),dsubq(I,K),dsubt(i,k)
1552 endif
1553 ENDIF
1554 enddo
1555 enddo
1556 do i=its,itf
1557 if(ierr(i).eq.0)then
1558 XHE(I,ktf)=HEO(I,ktf)
1559 XQ(I,ktf)=QO(I,ktf)
1560 XT(I,ktf)=TN(I,ktf)
1561 IF(XQ(I,ktf).LE.0.)XQ(I,ktf)=1.E-08
1562 endif
1563 enddo
1564 !
1565 !--- calculate moist static energy, heights, qes
1566 !
1567 call cup_env(xz,xqes,xhe,xhes,xt,xq,po,z1, &
1568 psur,ierr,tcrit,2,xl,cp, &
1569 itf,jtf,ktf, &
1570 its,ite, jts,jte, kts,kte)
1571 !
1572 !--- environmental values on cloud levels
1573 !
1574 call cup_env_clev(xt,xqes,xq,xhe,xhes,xz,po,xqes_cup,xq_cup, &
1575 xhe_cup,xhes_cup,xz_cup,po_cup,gamma_cup,xt_cup,psur, &
1576 ierr,z1,xl,rv,cp, &
1577 itf,jtf,ktf, &
1578 its,ite, jts,jte, kts,kte)
1579 !
1580 !
1581 !**************************** static control
1582 !
1583 !--- moist static energy inside cloud
1584 !
1585 do i=its,itf
1586 if(ierr(i).eq.0)then
1587 xhkb(i)=xhe(i,k22(i))
1588 endif
1589 enddo
1590 call cup_up_he(k22,xhkb,xz_cup,cd,mentr_rate,xhe_cup,xhc, &
1591 kbcon,ierr,xdby,xhe,xhes_cup, &
1592 itf,jtf,ktf, &
1593 its,ite, jts,jte, kts,kte)
1594 !
1595 !c--- normalized mass flux profile
1596 !
1597 call cup_up_nms(xzu,xz_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
1598 itf,jtf,ktf, &
1599 its,ite, jts,jte, kts,kte)
1600 !
1601 !--- moisture downdraft
1602 !
1603 call cup_dd_nms(xzd,xz_cup,cdd,mentrd_rate,jmin,ierr, &
1604 1,kdet,z1, &
1605 itf,jtf,ktf, &
1606 its,ite, jts,jte, kts,kte)
1607 call cup_dd_he(xhes_cup,xzd,xhcd,xz_cup,cdd,mentrd_rate, &
1608 jmin,ierr,xhe,dbyd,xhe_cup,&
1609 itf,jtf,ktf, &
1610 its,ite, jts,jte, kts,kte)
1611 call cup_dd_moisture_3d(xzd,xhcd,xhes_cup,xqcd,xqes_cup, &
1612 xpwd,xq_cup,xz_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
1613 xpwev,bu,xqrcd,xq,xhe,xt_cup,3,xl,high_resolution, &
1614 itf,jtf,ktf, &
1615 its,ite, jts,jte, kts,kte)
1617 !
1618 !------- MOISTURE updraft
1619 !
1620 call cup_up_moisture('deep',ierr,xz_cup,xqc,xqrc,xpw,xpwav, &
1621 kbcon,ktop,cd,xdby,mentr_rate,clw_all, &
1622 xq,GAMMA_cup,xzu,xqes_cup,k22,xq_cup,xl, &
1623 itf,jtf,ktf, &
1624 its,ite, jts,jte, kts,kte)
1625 !
1626 !--- workfunctions for updraft
1627 !
1628 call cup_up_aa0(xaa0,xz,xzu,xdby,GAMMA_CUP,xt_cup, &
1629 kbcon,ktop,ierr, &
1630 itf,jtf,ktf, &
1631 its,ite, jts,jte, kts,kte)
1632 do 200 nens=1,maxens
1633 do i=its,itf
1634 if(ierr(i).eq.0)then
1635 xaa0_ens(i,nens)=xaa0(i)
1636 nall=(iens-1)*maxens3*maxens*maxens2 &
1637 +(iedt-1)*maxens*maxens3 &
1638 +(nens-1)*maxens3
1639 do k=kts,ktf
1640 if(k.le.ktop(i))then
1641 do nens3=1,maxens3
1642 if(nens3.eq.7)then
1643 !--- b=0
1644 pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3) &
1645 +edto(i)*pwdo(i,k) &
1646 +pwo(i,k)
1647 !--- b=beta
1648 else if(nens3.eq.8)then
1649 pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
1650 pwo(i,k)
1651 !--- b=beta/2
1652 else if(nens3.eq.9)then
1653 pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3) &
1654 +.5*edto(i)*pwdo(i,k) &
1655 + pwo(i,k)
1656 else
1657 pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
1658 pwo(i,k)+edto(i)*pwdo(i,k)
1659 endif
1660 enddo
1661 endif
1662 enddo
1663 if(pr_ens(i,j,nall+7).lt.1.e-6)then
1664 ierr(i)=18
1665 do nens3=1,maxens3
1666 pr_ens(i,j,nall+nens3)=0.
1667 enddo
1668 endif
1669 do nens3=1,maxens3
1670 if(pr_ens(i,j,nall+nens3).lt.1.e-4)then
1671 pr_ens(i,j,nall+nens3)=0.
1672 endif
1673 enddo
1674 endif
1675 enddo
1676 200 continue
1677 !
1678 !--- LARGE SCALE FORCING
1679 !
1680 !
1681 !------- CHECK wether aa0 should have been zero
1682 !
1683 !
1684 CALL cup_MAXIMI(HEO_CUP,3,KBMAX,K22x,ierr, &
1685 itf,jtf,ktf, &
1686 its,ite, jts,jte, kts,kte)
1687 do i=its,itf
1688 ierr2(i)=ierr(i)
1689 ierr3(i)=ierr(i)
1690 enddo
1691 call cup_kbcon(cap_max_increment,2,k22x,kbconx,heo_cup, &
1692 heso_cup,ierr2,kbmax,po_cup,cap_max, &
1693 itf,jtf,ktf, &
1694 its,ite, jts,jte, kts,kte)
1695 call cup_kbcon(cap_max_increment,3,k22x,kbconx,heo_cup, &
1696 heso_cup,ierr3,kbmax,po_cup,cap_max, &
1697 itf,jtf,ktf, &
1698 its,ite, jts,jte, kts,kte)
1699 !
1700 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
1701 !
1703 call cup_forcing_ens_3d(closure_n,xland1,aa0,aa1,xaa0_ens,mbdt_ens,dtime, &
1704 ierr,ierr2,ierr3,xf_ens,j,'deeps',axx, &
1705 maxens,iens,iedt,maxens2,maxens3,mconv, &
1706 po_cup,ktop,omeg,zdo,k22,zuo,pr_ens,edto,kbcon, &
1707 massflx,iact,direction,ensdim,massfln,ichoice,edt_out, &
1708 high_resolution,itf,jtf,ktf, &
1709 its,ite, jts,jte, kts,kte,ens4,ktau)
1710 !
1711 do k=kts,ktf
1712 do i=its,itf
1713 if(ierr(i).eq.0)then
1714 subt_ens(i,k,iedt)=dsubt(i,k)
1715 subq_ens(i,k,iedt)=dsubq(i,k)
1716 dellat_ens(i,k,iedt)=dellat(i,k)
1717 dellaq_ens(i,k,iedt)=dellaq(i,k)
1718 dellaqc_ens(i,k,iedt)=dellaqc(i,k)
1719 pwo_ens(i,k,iedt)=pwo(i,k)+edt(i)*pwdo(i,k)
1720 else
1721 subt_ens(i,k,iedt)=0.
1722 subq_ens(i,k,iedt)=0.
1723 dellat_ens(i,k,iedt)=0.
1724 dellaq_ens(i,k,iedt)=0.
1725 dellaqc_ens(i,k,iedt)=0.
1726 pwo_ens(i,k,iedt)=0.
1727 endif
1728 if(i.eq.ipr.and.j.eq.jpr)then
1729 write(0,*)'1',iens,iedt,dellat(i,k),dellat_ens(i,k,iedt), &
1730 dellaq(i,k), dellaqc(i,k)
1731 write(0,*)'2',k,subt_ens(i,k,iedt),subq_ens(i,k,iedt)
1732 endif
1733 enddo
1734 enddo
1735 250 continue
1736 !
1737 !--- FEEDBACK
1738 !
1739 call cup_output_ens_3d(xf_ens,ierr,dellat_ens,dellaq_ens, &
1740 dellaqc_ens,subt_ens,subq_ens,subt,subq,outt, &
1741 outq,outqc,zuo,sub_mas,pre,pwo_ens,xmb,ktop, &
1742 j,'deep',maxens2,maxens,iens,ierr2,ierr3, &
1743 pr_ens,maxens3,ensdim,massfln, &
1744 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
1745 APR_CAPMA,APR_CAPME,APR_CAPMI,closure_n,xland1, &
1746 itf,jtf,ktf, &
1747 its,ite, jts,jte, kts,kte)
1748 k=1
1749 do i=its,itf
1750 if(ierr(i).eq.0.and.ierr5(i).eq.0.and.kbcon(i).lt.ktop3(i)+1)then
1751 ! write(0,*)'both ier and ier5=0 at i,j=',i,j,kbcon(i),ktop3(i)
1752 if(high_resolution.eq.1)then
1753 outts(i,kts:kte)=0.
1754 outqs(i,kts:kte)=0.
1755 endif
1756 elseif (ierr5(i).eq.0)then
1757 ! write(0,*)'ier5=0 at i,j=',i,j,k23(i),ktop3(i)
1758 endif
1760 PRE(I)=MAX(PRE(I),0.)
1761 if(i.eq.ipr.and.j.eq.jpr)then
1762 write(0,*)'i,j,pre(i),aa0(i),aa1(i)'
1763 write(0,*)i,j,pre(i),aa0(i)
1764 endif
1765 enddo
1766 !
1767 !---------------------------done------------------------------
1768 !
1769 do i=its,itf
1770 if(ierr(i).eq.0)then
1771 if(i.eq.ipr.and.j.eq.jpr)then
1772 write(0,*)'on output, pre =',pre(i),its,itf,kts,ktf
1773 do k=kts,ktf
1774 write(0,*)z(i,k),outt(i,k)*86400.,subt(i,k)*86400.
1775 enddo
1776 write(0,*)i,j,(axx(i,k),k=1,ens4)
1777 endif
1778 endif
1779 enddo
1780 ! print *,'ierr(i) = ',ierr(i),pre(i)
1782 END SUBROUTINE CUP_enss_3d
1785 SUBROUTINE cup_dd_aa0(edt,ierr,aa0,jmin,gamma_cup,t_cup, &
1786 hcd,hes_cup,z,zd, &
1787 itf,jtf,ktf, &
1788 its,ite, jts,jte, kts,kte )
1790 IMPLICIT NONE
1791 !
1792 ! on input
1793 !
1795 ! only local wrf dimensions are need as of now in this routine
1797 integer &
1798 ,intent (in ) :: &
1799 itf,jtf,ktf, &
1800 its,ite, jts,jte, kts,kte
1801 ! aa0 cloud work function for downdraft
1802 ! gamma_cup = gamma on model cloud levels
1803 ! t_cup = temperature (Kelvin) on model cloud levels
1804 ! hes_cup = saturation moist static energy on model cloud levels
1805 ! hcd = moist static energy in downdraft
1806 ! edt = epsilon
1807 ! zd normalized downdraft mass flux
1808 ! z = heights of model levels
1809 ! ierr error value, maybe modified in this routine
1810 !
1811 real, dimension (its:ite,kts:kte) &
1812 ,intent (in ) :: &
1813 z,zd,gamma_cup,t_cup,hes_cup,hcd
1814 real, dimension (its:ite) &
1815 ,intent (in ) :: &
1816 edt
1817 integer, dimension (its:ite) &
1818 ,intent (in ) :: &
1819 jmin
1820 !
1821 ! input and output
1822 !
1825 integer, dimension (its:ite) &
1826 ,intent (inout) :: &
1827 ierr
1828 real, dimension (its:ite) &
1829 ,intent (out ) :: &
1830 aa0
1831 !
1832 ! local variables in this routine
1833 !
1835 integer :: &
1836 i,k,kk
1837 real :: &
1838 dz
1839 !
1840 do i=its,itf
1841 aa0(i)=0.
1842 enddo
1843 !
1844 !?? DO k=kts,kte-1
1845 DO k=kts,ktf-1
1846 do i=its,itf
1847 IF(ierr(I).eq.0.and.k.lt.jmin(i))then
1848 KK=JMIN(I)-K
1849 !
1850 !--- ORIGINAL
1851 !
1852 DZ=(Z(I,KK)-Z(I,KK+1))
1853 AA0(I)=AA0(I)+zd(i,kk)*EDT(I)*DZ*(9.81/(1004.*T_cup(I,KK))) &
1854 *((hcd(i,kk)-hes_cup(i,kk))/(1.+GAMMA_cup(i,kk)))
1855 endif
1856 enddo
1857 enddo
1859 END SUBROUTINE CUP_dd_aa0
1862 SUBROUTINE cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, &
1863 pwev,edtmax,edtmin,maxens2,edtc, &
1864 itf,jtf,ktf, &
1865 its,ite, jts,jte, kts,kte )
1867 IMPLICIT NONE
1869 integer &
1870 ,intent (in ) :: &
1871 itf,jtf,ktf, &
1872 its,ite, jts,jte, kts,kte
1873 integer, intent (in ) :: &
1874 maxens2
1875 !
1876 ! ierr error value, maybe modified in this routine
1877 !
1878 real, dimension (its:ite,kts:kte) &
1879 ,intent (in ) :: &
1880 us,vs,z,p
1881 real, dimension (its:ite,1:maxens2) &
1882 ,intent (out ) :: &
1883 edtc
1884 real, dimension (its:ite) &
1885 ,intent (out ) :: &
1886 edt
1887 real, dimension (its:ite) &
1888 ,intent (in ) :: &
1889 pwav,pwev
1890 real &
1891 ,intent (in ) :: &
1892 edtmax,edtmin
1893 integer, dimension (its:ite) &
1894 ,intent (in ) :: &
1895 ktop,kbcon
1896 integer, dimension (its:ite) &
1897 ,intent (inout) :: &
1898 ierr
1899 !
1900 ! local variables in this routine
1901 !
1903 integer i,k,kk
1904 real einc,pef,pefb,prezk,zkbc
1905 real, dimension (its:ite) :: &
1906 vshear,sdp,vws
1908 !
1909 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
1910 !
1911 ! */ calculate an average wind shear over the depth of the cloud
1912 !
1913 do i=its,itf
1914 edt(i)=0.
1915 vws(i)=0.
1916 sdp(i)=0.
1917 vshear(i)=0.
1918 enddo
1919 do k=1,maxens2
1920 do i=its,itf
1921 edtc(i,k)=0.
1922 enddo
1923 enddo
1924 do kk = kts,ktf-1
1925 do 62 i=its,itf
1926 IF(ierr(i).ne.0)GO TO 62
1927 if (kk .le. min0(ktop(i),ktf) .and. kk .ge. kbcon(i)) then
1928 vws(i) = vws(i)+ &
1929 (abs((us(i,kk+1)-us(i,kk))/(z(i,kk+1)-z(i,kk))) &
1930 + abs((vs(i,kk+1)-vs(i,kk))/(z(i,kk+1)-z(i,kk)))) * &
1931 (p(i,kk) - p(i,kk+1))
1932 sdp(i) = sdp(i) + p(i,kk) - p(i,kk+1)
1933 endif
1934 if (kk .eq. ktf)vshear(i) = 1.e3 * vws(i) / sdp(i)
1935 62 continue
1936 end do
1937 do i=its,itf
1938 IF(ierr(i).eq.0)then
1939 pef=(1.591-.639*VSHEAR(I)+.0953*(VSHEAR(I)**2) &
1940 -.00496*(VSHEAR(I)**3))
1941 if(pef.gt.1.)pef=1.
1942 if(pef.lt.0.)pef=0.
1943 !
1944 !--- cloud base precip efficiency
1945 !
1946 zkbc=z(i,kbcon(i))*3.281e-3
1947 prezk=.02
1948 if(zkbc.gt.3.)then
1949 prezk=.96729352+zkbc*(-.70034167+zkbc*(.162179896+zkbc &
1950 *(- 1.2569798E-2+zkbc*(4.2772E-4-zkbc*5.44E-6))))
1951 endif
1952 if(zkbc.gt.25)then
1953 prezk=2.4
1954 endif
1955 pefb=1./(1.+prezk)
1956 if(pefb.gt.1.)pefb=1.
1957 if(pefb.lt.0.)pefb=0.
1958 EDT(I)=1.-.5*(pefb+pef)
1959 !--- edt here is 1-precipeff!
1960 einc=.2*edt(i)
1961 do k=1,maxens2
1962 edtc(i,k)=edt(i)+float(k-2)*einc
1963 enddo
1964 endif
1965 enddo
1966 do i=its,itf
1967 IF(ierr(i).eq.0)then
1968 do k=1,maxens2
1969 EDTC(I,K)=-EDTC(I,K)*PWAV(I)/PWEV(I)
1970 IF(EDTC(I,K).GT.edtmax)EDTC(I,K)=edtmax
1971 IF(EDTC(I,K).LT.edtmin)EDTC(I,K)=edtmin
1972 enddo
1973 endif
1974 enddo
1976 END SUBROUTINE cup_dd_edt
1979 SUBROUTINE cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,entr, &
1980 jmin,ierr,he,dby,he_cup, &
1981 itf,jtf,ktf, &
1982 its,ite, jts,jte, kts,kte )
1984 IMPLICIT NONE
1985 !
1986 ! on input
1987 !
1989 ! only local wrf dimensions are need as of now in this routine
1991 integer &
1992 ,intent (in ) :: &
1993 itf,jtf,ktf, &
1994 its,ite, jts,jte, kts,kte
1995 ! hcd = downdraft moist static energy
1996 ! he = moist static energy on model levels
1997 ! he_cup = moist static energy on model cloud levels
1998 ! hes_cup = saturation moist static energy on model cloud levels
1999 ! dby = buoancy term
2000 ! cdd= detrainment function
2001 ! z_cup = heights of model cloud levels
2002 ! entr = entrainment rate
2003 ! zd = downdraft normalized mass flux
2004 !
2005 real, dimension (its:ite,kts:kte) &
2006 ,intent (in ) :: &
2007 he,he_cup,hes_cup,z_cup,cdd,zd
2008 ! entr= entrainment rate
2009 real &
2010 ,intent (in ) :: &
2011 entr
2012 integer, dimension (its:ite) &
2013 ,intent (in ) :: &
2014 jmin
2015 !
2016 ! input and output
2017 !
2019 ! ierr error value, maybe modified in this routine
2021 integer, dimension (its:ite) &
2022 ,intent (inout) :: &
2023 ierr
2025 real, dimension (its:ite,kts:kte) &
2026 ,intent (out ) :: &
2027 hcd,dby
2028 !
2029 ! local variables in this routine
2030 !
2032 integer :: &
2033 i,k,ki
2034 real :: &
2035 dz
2038 do k=kts+1,ktf
2039 do i=its,itf
2040 dby(i,k)=0.
2041 IF(ierr(I).eq.0)then
2042 hcd(i,k)=hes_cup(i,k)
2043 endif
2044 enddo
2045 enddo
2046 !
2047 do 100 i=its,itf
2048 IF(ierr(I).eq.0)then
2049 k=jmin(i)
2050 hcd(i,k)=hes_cup(i,k)
2051 dby(i,k)=hcd(i,jmin(i))-hes_cup(i,k)
2052 !
2053 do ki=jmin(i)-1,1,-1
2054 DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
2055 HCD(i,Ki)=(HCD(i,Ki+1)*(1.-.5*CDD(i,Ki)*DZ) &
2056 +entr*DZ*HE(i,Ki) &
2057 )/(1.+entr*DZ-.5*CDD(i,Ki)*DZ)
2058 dby(i,ki)=HCD(i,Ki)-hes_cup(i,ki)
2059 enddo
2060 !
2061 endif
2062 !--- end loop over i
2063 100 continue
2066 END SUBROUTINE cup_dd_he
2069 SUBROUTINE cup_dd_moisture_3d(zd,hcd,hes_cup,qcd,qes_cup, &
2070 pwd,q_cup,z_cup,cdd,entr,jmin,ierr, &
2071 gamma_cup,pwev,bu,qrcd, &
2072 q,he,t_cup,iloop,xl,high_resolution, &
2073 itf,jtf,ktf, &
2074 its,ite, jts,jte, kts,kte )
2076 IMPLICIT NONE
2078 integer &
2079 ,intent (in ) :: &
2080 itf,jtf,ktf, &
2081 its,ite, jts,jte, kts,kte,high_resolution
2082 ! cdd= detrainment function
2083 ! q = environmental q on model levels
2084 ! q_cup = environmental q on model cloud levels
2085 ! qes_cup = saturation q on model cloud levels
2086 ! hes_cup = saturation h on model cloud levels
2087 ! hcd = h in model cloud
2088 ! bu = buoancy term
2089 ! zd = normalized downdraft mass flux
2090 ! gamma_cup = gamma on model cloud levels
2091 ! mentr_rate = entrainment rate
2092 ! qcd = cloud q (including liquid water) after entrainment
2093 ! qrch = saturation q in cloud
2094 ! pwd = evaporate at that level
2095 ! pwev = total normalized integrated evaoprate (I2)
2096 ! entr= entrainment rate
2097 !
2098 real, dimension (its:ite,kts:kte) &
2099 ,intent (in ) :: &
2100 zd,t_cup,hes_cup,hcd,qes_cup,q_cup,z_cup,cdd,gamma_cup,q,he
2101 real &
2102 ,intent (in ) :: &
2103 entr,xl
2104 integer &
2105 ,intent (in ) :: &
2106 iloop
2107 integer, dimension (its:ite) &
2108 ,intent (in ) :: &
2109 jmin
2110 integer, dimension (its:ite) &
2111 ,intent (inout) :: &
2112 ierr
2113 real, dimension (its:ite,kts:kte) &
2114 ,intent (out ) :: &
2115 qcd,qrcd,pwd
2116 real, dimension (its:ite) &
2117 ,intent (out ) :: &
2118 pwev,bu
2119 !
2120 ! local variables in this routine
2121 !
2123 integer :: &
2124 i,k,ki
2125 real :: &
2126 dh,dz,dqeva
2128 do i=its,itf
2129 bu(i)=0.
2130 pwev(i)=0.
2131 enddo
2132 do k=kts,ktf
2133 do i=its,itf
2134 qcd(i,k)=0.
2135 qrcd(i,k)=0.
2136 pwd(i,k)=0.
2137 enddo
2138 enddo
2139 !
2140 !
2141 !
2142 do 100 i=its,itf
2143 IF(ierr(I).eq.0)then
2144 k=jmin(i)
2145 DZ=Z_cup(i,K+1)-Z_cup(i,K)
2146 qcd(i,k)=q_cup(i,k)
2147 if(high_resolution.eq.1)qcd(i,k)=.5*(qes_cup(i,k)+q_cup(i,k))
2148 qrcd(i,k)=qes_cup(i,k)
2149 pwd(i,jmin(i))=min(0.,qcd(i,k)-qrcd(i,k))
2150 pwev(i)=pwev(i)+pwd(i,jmin(i))
2151 qcd(i,k)=qes_cup(i,k)
2152 !
2153 DH=HCD(I,k)-HES_cup(I,K)
2154 bu(i)=dz*dh
2155 do ki=jmin(i)-1,1,-1
2156 DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
2157 QCD(i,Ki)=(qCD(i,Ki+1)*(1.-.5*CDD(i,Ki)*DZ) &
2158 +entr*DZ*q(i,Ki) &
2159 )/(1.+entr*DZ-.5*CDD(i,Ki)*DZ)
2160 !
2161 !--- to be negatively buoyant, hcd should be smaller than hes!
2162 !
2163 DH=HCD(I,ki)-HES_cup(I,Ki)
2164 bu(i)=bu(i)+dz*dh
2165 QRCD(I,Ki)=qes_cup(i,ki)+(1./XL)*(GAMMA_cup(i,ki) &
2166 /(1.+GAMMA_cup(i,ki)))*DH
2167 dqeva=qcd(i,ki)-qrcd(i,ki)
2168 if(dqeva.gt.0.)dqeva=0.
2169 pwd(i,ki)=zd(i,ki)*dqeva
2170 qcd(i,ki)=qrcd(i,ki)
2171 pwev(i)=pwev(i)+pwd(i,ki)
2172 ! if(iloop.eq.1.and.i.eq.102.and.j.eq.62)then
2173 ! print *,'in cup_dd_moi ', hcd(i,ki),HES_cup(I,Ki),dh,dqeva
2174 ! endif
2175 enddo
2176 !
2177 !--- end loop over i
2178 if(pwev(I).eq.0.and.iloop.eq.1)then
2179 ! print *,'problem with buoy in cup_dd_moisture',i
2180 ierr(i)=7
2181 endif
2182 if(BU(I).GE.0.and.iloop.eq.1)then
2183 ! print *,'problem with buoy in cup_dd_moisture',i
2184 ierr(i)=7
2185 endif
2186 endif
2187 100 continue
2189 END SUBROUTINE cup_dd_moisture_3d
2192 SUBROUTINE cup_dd_nms(zd,z_cup,cdd,entr,jmin,ierr, &
2193 itest,kdet,z1, &
2194 itf,jtf,ktf, &
2195 its,ite, jts,jte, kts,kte )
2197 IMPLICIT NONE
2198 !
2199 ! on input
2200 !
2202 ! only local wrf dimensions are need as of now in this routine
2204 integer &
2205 ,intent (in ) :: &
2206 itf,jtf,ktf, &
2207 its,ite, jts,jte, kts,kte
2208 ! z_cup = height of cloud model level
2209 ! z1 = terrain elevation
2210 ! entr = downdraft entrainment rate
2211 ! jmin = downdraft originating level
2212 ! kdet = level above ground where downdraft start detraining
2213 ! itest = flag to whether to calculate cdd
2215 real, dimension (its:ite,kts:kte) &
2216 ,intent (in ) :: &
2217 z_cup
2218 real, dimension (its:ite) &
2219 ,intent (in ) :: &
2220 z1
2221 real &
2222 ,intent (in ) :: &
2223 entr
2224 integer, dimension (its:ite) &
2225 ,intent (in ) :: &
2226 jmin,kdet
2227 integer &
2228 ,intent (in ) :: &
2229 itest
2230 !
2231 ! input and output
2232 !
2234 ! ierr error value, maybe modified in this routine
2236 integer, dimension (its:ite) &
2237 ,intent (inout) :: &
2238 ierr
2239 ! zd is the normalized downdraft mass flux
2240 ! cdd is the downdraft detrainmen function
2242 real, dimension (its:ite,kts:kte) &
2243 ,intent (out ) :: &
2244 zd
2245 real, dimension (its:ite,kts:kte) &
2246 ,intent (inout) :: &
2247 cdd
2248 !
2249 ! local variables in this routine
2250 !
2252 integer :: &
2253 i,k,ki
2254 real :: &
2255 a,perc,dz
2257 !
2258 !--- perc is the percentage of mass left when hitting the ground
2259 !
2260 perc=.03
2262 do k=kts,ktf
2263 do i=its,itf
2264 zd(i,k)=0.
2265 if(itest.eq.0)cdd(i,k)=0.
2266 enddo
2267 enddo
2268 a=1.-perc
2269 !
2270 !
2271 !
2272 do 100 i=its,itf
2273 IF(ierr(I).eq.0)then
2274 zd(i,jmin(i))=1.
2275 !
2276 !--- integrate downward, specify detrainment(cdd)!
2277 !
2278 do ki=jmin(i)-1,1,-1
2279 DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
2280 if(ki.le.kdet(i).and.itest.eq.0)then
2281 cdd(i,ki)=entr+(1.- (a*(z_cup(i,ki)-z1(i)) &
2282 +perc*(z_cup(i,kdet(i))-z1(i)) ) &
2283 /(a*(z_cup(i,ki+1)-z1(i)) &
2284 +perc*(z_cup(i,kdet(i))-z1(i))))/dz
2285 endif
2286 zd(i,ki)=zd(i,ki+1)*(1.+(entr-cdd(i,ki))*dz)
2287 enddo
2288 !
2289 endif
2290 !--- end loop over i
2291 100 continue
2293 END SUBROUTINE cup_dd_nms
2296 SUBROUTINE cup_dellabot(name,ipr,jpr,he_cup,ierr,z_cup,p_cup, &
2297 hcd,edt,zd,cdd,he,della,subs,j,mentrd_rate,z,g, &
2298 itf,jtf,ktf, &
2299 its,ite, jts,jte, kts,kte )
2301 IMPLICIT NONE
2303 integer &
2304 ,intent (in ) :: &
2305 itf,jtf,ktf, &
2306 its,ite, jts,jte, kts,kte
2307 integer, intent (in ) :: &
2308 j,ipr,jpr
2309 character *(*), intent (in) :: &
2310 name
2311 !
2312 ! ierr error value, maybe modified in this routine
2313 !
2314 real, dimension (its:ite,kts:kte) &
2315 ,intent (out ) :: &
2316 della,subs
2317 real, dimension (its:ite,kts:kte) &
2318 ,intent (in ) :: &
2319 z_cup,p_cup,hcd,zd,cdd,he,z,he_cup
2320 real, dimension (its:ite) &
2321 ,intent (in ) :: &
2322 edt
2323 real &
2324 ,intent (in ) :: &
2325 g,mentrd_rate
2326 integer, dimension (its:ite) &
2327 ,intent (inout) :: &
2328 ierr
2329 !
2330 ! local variables in this routine
2331 !
2333 integer i
2334 real detdo,detdo1,detdo2,entdo,dp,dz,subin, &
2335 totmas
2336 !
2337 !
2338 ! if(name.eq.'shallow')then
2339 ! edt(:)=0.
2340 ! cdd(:,:)=0.
2341 ! endif
2342 do 100 i=its,itf
2343 della(i,1)=0.
2344 subs(i,1)=0.
2345 if(ierr(i).ne.0)go to 100
2346 dz=z_cup(i,2)-z_cup(i,1)
2347 DP=100.*(p_cup(i,1)-P_cup(i,2))
2348 detdo1=edt(i)*zd(i,2)*CDD(i,1)*DZ
2349 detdo2=edt(i)*zd(i,1)
2350 entdo=edt(i)*zd(i,2)*mentrd_rate*dz
2351 subin=-EDT(I)*zd(i,2)
2352 detdo=detdo1+detdo2-entdo+subin
2353 DELLA(I,1)=(detdo1*.5*(HCD(i,1)+HCD(i,2)) &
2354 +detdo2*hcd(i,1) &
2355 +subin*he_cup(i,2) &
2356 -entdo*he(i,1))*g/dp
2357 SUBS(I,1)=0.
2358 if(i.eq.ipr.and.j.eq.jpr)then
2359 write(0,*)'db1',della(i,1),subs(i,1),subin,entdo
2360 write(0,*)'db2',detdo1,detdo2,detdo1+detdo2-entdo+subin
2361 endif
2362 100 CONTINUE
2364 END SUBROUTINE cup_dellabot
2367 SUBROUTINE cup_dellas_3d(ierr,z_cup,p_cup,hcd,edt,zd,cdd, &
2368 he,della,subs,j,mentrd_rate,zu,g, &
2369 cd,hc,ktop,k22,kbcon,mentr_rate,jmin,he_cup,kdet,kpbl, &
2370 ipr,jpr,name,high_res, &
2371 itf,jtf,ktf, &
2372 its,ite, jts,jte, kts,kte )
2374 IMPLICIT NONE
2376 integer &
2377 ,intent (in ) :: &
2378 itf,jtf,ktf, &
2379 its,ite, jts,jte, kts,kte
2380 integer, intent (in ) :: &
2381 j,ipr,jpr,high_res
2382 !
2383 ! ierr error value, maybe modified in this routine
2384 !
2385 real, dimension (its:ite,kts:kte) &
2386 ,intent (out ) :: &
2387 della,subs
2388 real, dimension (its:ite,kts:kte) &
2389 ,intent (in ) :: &
2390 z_cup,p_cup,hcd,zd,cdd,he,hc,cd,zu,he_cup
2391 real, dimension (its:ite) &
2392 ,intent (in ) :: &
2393 edt
2394 real &
2395 ,intent (in ) :: &
2396 g,mentrd_rate,mentr_rate
2397 integer, dimension (its:ite) &
2398 ,intent (in ) :: &
2399 kbcon,ktop,k22,jmin,kdet,kpbl
2400 integer, dimension (its:ite) &
2401 ,intent (inout) :: &
2402 ierr
2403 character *(*), intent (in) :: &
2404 name
2405 !
2406 ! local variables in this routine
2407 !
2409 integer i,k,kstart
2410 real detdo1,detdo2,entdo,dp,dz,subin,detdo,entup, &
2411 detup,subdown,entdoj,entupk,detupk,totmas
2412 !
2413 i=ipr
2414 kstart=kts+1
2415 if(name.eq.'shallow')kstart=kts
2416 DO K=kstart,ktf
2417 do i=its,itf
2418 della(i,k)=0.
2419 subs(i,k)=0.
2420 enddo
2421 enddo
2422 !
2423 ! no downdrafts for shallow convection
2424 !
2425 DO 100 k=kts+1,ktf-1
2426 DO 100 i=its,ite
2427 IF(ierr(i).ne.0)GO TO 100
2428 IF(K.Gt.KTOP(I))GO TO 100
2429 if(k.lt.k22(i)-1.and.name.eq.'shallow')GO TO 100
2430 !
2431 !--- SPECIFY DETRAINMENT OF DOWNDRAFT, HAS TO BE CONSISTENT
2432 !--- WITH ZD CALCULATIONS IN SOUNDD.
2433 !
2434 DZ=Z_cup(I,K+1)-Z_cup(I,K)
2435 detdo=edt(i)*CDD(i,K)*DZ*ZD(i,k+1)
2436 entdo=edt(i)*mentrd_rate*dz*zd(i,k+1)
2437 !3d subin=zu(i,k+1)-zd(i,k+1)*edt(i)
2438 subin=-zd(i,k+1)*edt(i)
2439 entup=0.
2440 detup=0.
2441 if(k.ge.kbcon(i).and.k.lt.ktop(i))then
2442 entup=mentr_rate*dz*zu(i,k)
2443 detup=CD(i,K+1)*DZ*ZU(i,k)
2444 endif
2445 !3d subdown=(zu(i,k)-zd(i,k)*edt(i))
2446 subdown=-zd(i,k)*edt(i)
2447 entdoj=0.
2448 entupk=0.
2449 detupk=0.
2450 !
2451 if(k.eq.jmin(i))then
2452 entdoj=edt(i)*zd(i,k)
2453 endif
2455 if(k.eq.k22(i)-1)then
2456 entupk=zu(i,kpbl(i))
2457 subin=zu(i,k+1)-zd(i,k+1)*edt(i)
2458 if(high_res.eq.1)subin=-zd(i,k+1)*edt(i)
2459 ! subin=-zd(i,k+1)*edt(i)
2460 endif
2462 if(k.gt.kdet(i))then
2463 detdo=0.
2464 endif
2466 if(k.eq.ktop(i)-0)then
2467 detupk=zu(i,ktop(i))
2468 subin=0.
2469 !
2470 ! this subsidene for ktop now in subs term!
2471 ! subdown=zu(i,k)
2472 subdown=0.
2473 endif
2474 if(k.lt.kbcon(i))then
2475 detup=0.
2476 endif
2477 !C
2478 !C--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
2479 !C
2480 totmas=subin-subdown+detup-entup-entdo+ &
2481 detdo-entupk-entdoj+detupk
2482 ! if(j.eq.jpr.and.i.eq.ipr)print *,'k,totmas,sui,sud = ',k,
2483 ! 1 totmas,subin,subdown
2484 ! if(j.eq.jpr.and.i.eq.ipr)print *,'updr stuff = ',detup,
2485 ! 1 entup,entupk,detupk
2486 ! if(j.eq.jpr.and.i.eq.ipr)print *,'dddr stuff = ',entdo,
2487 ! 1 detdo,entdoj
2488 if(abs(totmas).gt.1.e-6)then
2489 ! print *,'*********************',i,j,k,totmas,name
2490 ! print *,kpbl(i),k22(i),kbcon(i),ktop(i)
2491 !c print *,'updr stuff = ',subin,
2492 !c 1 subdown,detup,entup,entupk,detupk
2493 !c print *,'dddr stuff = ',entdo,
2494 !c 1 detdo,entdoj
2495 ! call wrf_error_fatal ( 'totmas .gt.1.e-6' )
2496 endif
2497 dp=100.*(p_cup(i,k-1)-p_cup(i,k))
2498 della(i,k)=(detup*.5*(HC(i,K+1)+HC(i,K)) &
2499 +detdo*.5*(HCD(i,K+1)+HCD(i,K)) &
2500 -entup*he(i,k) &
2501 -entdo*he(i,k) &
2502 +subin*he_cup(i,k+1) &
2503 -subdown*he_cup(i,k) &
2504 +detupk*(hc(i,ktop(i))-he_cup(i,ktop(i))) &
2505 -entupk*he_cup(i,k22(i)) &
2506 -entdoj*he_cup(i,jmin(i)) &
2507 )*g/dp
2508 if(high_res.eq.1)then
2509 ! the first term includes entr and detr into/from updraft as well as (entup-detup)*he(i,k) from
2510 ! neighbouring point, to make things mass consistent....
2511 ! if(k.ge.k22(i))then
2512 della(i,k)=( &
2513 detup*.5*(HC(i,K+1)+HC(i,K))-entup*he(i,k)+(entup-detup)*he(i,k) &
2514 +detdo*.5*(HCD(i,K+1)+HCD(i,K)) &
2515 -entdo*he(i,k) &
2516 +subin*he_cup(i,k+1) &
2517 -subdown*he_cup(i,k) &
2518 +detupk*(hc(i,ktop(i))-he(i,ktop(i))) &
2519 -entdoj*he_cup(i,jmin(i)) &
2520 -entupk*he_cup(i,k22(i))+entupk*he(i,k) &
2521 )*g/dp
2522 ! else if(k.eq.k22(i)-1)then
2523 ! della(i,k)=(-entupk*he_cup(i,k22(i))+entupk*he(i,k))*g/dp
2524 endif
2525 !3d subin=zu(i,k+1)-zd(i,k+1)*edt(i)
2526 !
2527 ! updraft subsidence only
2528 !
2529 if(k.ge.k22(i).and.k.lt.ktop(i))then
2530 subs(i,k)=(zu(i,k+1)*he_cup(i,k+1) &
2531 -zu(i,k)*he_cup(i,k))*g/dp
2532 ! else if(k.eq.ktop(i))then
2533 ! subs(i,k)=-detupk*he_cup(i,k)*g/dp
2534 endif
2535 !
2536 ! in igh res case, subsidence terms are for meighbouring points only. This has to be
2537 ! done mass consistent with the della term
2538 if(high_res.eq.1)then
2539 if(k.ge.k22(i).and.k.lt.ktop(i))then
2540 subs(i,k)=(zu(i,k+1)*he_cup(i,k+1)-zu(i,k)*he_cup(i,k)-(entup-detup)*he(i,k))*g/dp
2541 else if(k.eq.ktop(i))then
2542 subs(i,k)=detupk*(he(i,ktop(i))-he_cup(i,ktop(i)))*g/dp
2543 else if(k.eq.k22(i)-1)then
2544 subs(i,k)=(entupk*he(i,k)-entupk*he_cup(i,k))*g/dp
2545 endif
2546 endif
2547 if(i.eq.ipr.and.j.eq.jpr)then
2548 write(0,*)'d',k,della(i,k),subs(i,k),subin,subdown
2549 ! write(0,*)'d',detup,entup,entdo,entupk,entdoj
2550 ! print *,k,della(i,k),subin*he_cup(i,k+1),subdown*he_cup(i,k),
2551 ! 1 detdo*.5*(HCD(i,K+1)+HCD(i,K))
2552 ! print *,k,detup*.5*(HC(i,K+1)+HC(i,K)),detupk*hc(i,ktop(i)),
2553 ! 1 entup*he(i,k),entdo*he(i,k)
2554 ! print *,k,he_cup(i,k+1),he_cup(i,k),entupk*he_cup(i,k)
2555 endif
2557 100 CONTINUE
2559 END SUBROUTINE cup_dellas_3d
2562 SUBROUTINE cup_direction2(i,j,dir,id,massflx, &
2563 iresult,imass,massfld, &
2564 itf,jtf,ktf, &
2565 its,ite, jts,jte, kts,kte )
2567 IMPLICIT NONE
2569 integer &
2570 ,intent (in ) :: &
2571 itf,jtf,ktf, &
2572 its,ite, jts,jte, kts,kte
2573 integer, intent (in ) :: &
2574 i,j,imass
2575 integer, intent (out ) :: &
2576 iresult
2577 !
2578 ! ierr error value, maybe modified in this routine
2579 !
2580 integer, dimension (its:ite,jts:jte) &
2581 ,intent (in ) :: &
2582 id
2583 real, dimension (its:ite,jts:jte) &
2584 ,intent (in ) :: &
2585 massflx
2586 real, dimension (its:ite) &
2587 ,intent (inout) :: &
2588 dir
2589 real &
2590 ,intent (out ) :: &
2591 massfld
2592 !
2593 ! local variables in this routine
2594 !
2596 integer k,ia,ja,ib,jb
2597 real diff
2598 !
2599 !
2600 !
2601 if(imass.eq.1)then
2602 massfld=massflx(i,j)
2603 endif
2604 iresult=0
2605 ! return
2606 diff=22.5
2607 if(dir(i).lt.22.5)dir(i)=360.+dir(i)
2608 if(id(i,j).eq.1)iresult=1
2609 ! ja=max(2,j-1)
2610 ! ia=max(2,i-1)
2611 ! jb=min(mjx-1,j+1)
2612 ! ib=min(mix-1,i+1)
2613 ja=j-1
2614 ia=i-1
2615 jb=j+1
2616 ib=i+1
2617 if(dir(i).gt.90.-diff.and.dir(i).le.90.+diff)then
2618 !--- steering flow from the east
2619 if(id(ib,j).eq.1)then
2620 iresult=1
2621 if(imass.eq.1)then
2622 massfld=max(massflx(ib,j),massflx(i,j))
2623 endif
2624 return
2625 endif
2626 else if(dir(i).gt.135.-diff.and.dir(i).le.135.+diff)then
2627 !--- steering flow from the south-east
2628 if(id(ib,ja).eq.1)then
2629 iresult=1
2630 if(imass.eq.1)then
2631 massfld=max(massflx(ib,ja),massflx(i,j))
2632 endif
2633 return
2634 endif
2635 !--- steering flow from the south
2636 else if(dir(i).gt.180.-diff.and.dir(i).le.180.+diff)then
2637 if(id(i,ja).eq.1)then
2638 iresult=1
2639 if(imass.eq.1)then
2640 massfld=max(massflx(i,ja),massflx(i,j))
2641 endif
2642 return
2643 endif
2644 !--- steering flow from the south west
2645 else if(dir(i).gt.225.-diff.and.dir(i).le.225.+diff)then
2646 if(id(ia,ja).eq.1)then
2647 iresult=1
2648 if(imass.eq.1)then
2649 massfld=max(massflx(ia,ja),massflx(i,j))
2650 endif
2651 return
2652 endif
2653 !--- steering flow from the west
2654 else if(dir(i).gt.270.-diff.and.dir(i).le.270.+diff)then
2655 if(id(ia,j).eq.1)then
2656 iresult=1
2657 if(imass.eq.1)then
2658 massfld=max(massflx(ia,j),massflx(i,j))
2659 endif
2660 return
2661 endif
2662 !--- steering flow from the north-west
2663 else if(dir(i).gt.305.-diff.and.dir(i).le.305.+diff)then
2664 if(id(ia,jb).eq.1)then
2665 iresult=1
2666 if(imass.eq.1)then
2667 massfld=max(massflx(ia,jb),massflx(i,j))
2668 endif
2669 return
2670 endif
2671 !--- steering flow from the north
2672 else if(dir(i).gt.360.-diff.and.dir(i).le.360.+diff)then
2673 if(id(i,jb).eq.1)then
2674 iresult=1
2675 if(imass.eq.1)then
2676 massfld=max(massflx(i,jb),massflx(i,j))
2677 endif
2678 return
2679 endif
2680 !--- steering flow from the north-east
2681 else if(dir(i).gt.45.-diff.and.dir(i).le.45.+diff)then
2682 if(id(ib,jb).eq.1)then
2683 iresult=1
2684 if(imass.eq.1)then
2685 massfld=max(massflx(ib,jb),massflx(i,j))
2686 endif
2687 return
2688 endif
2689 endif
2691 END SUBROUTINE cup_direction2
2694 SUBROUTINE cup_env(z,qes,he,hes,t,q,p,z1, &
2695 psur,ierr,tcrit,itest,xl,cp, &
2696 itf,jtf,ktf, &
2697 its,ite, jts,jte, kts,kte )
2699 IMPLICIT NONE
2701 integer &
2702 ,intent (in ) :: &
2703 itf,jtf,ktf, &
2704 its,ite, jts,jte, kts,kte
2705 !
2706 ! ierr error value, maybe modified in this routine
2707 ! q = environmental mixing ratio
2708 ! qes = environmental saturation mixing ratio
2709 ! t = environmental temp
2710 ! tv = environmental virtual temp
2711 ! p = environmental pressure
2712 ! z = environmental heights
2713 ! he = environmental moist static energy
2714 ! hes = environmental saturation moist static energy
2715 ! psur = surface pressure
2716 ! z1 = terrain elevation
2717 !
2718 !
2719 real, dimension (its:ite,kts:kte) &
2720 ,intent (in ) :: &
2721 p,t,q
2722 real, dimension (its:ite,kts:kte) &
2723 ,intent (out ) :: &
2724 he,hes,qes
2725 real, dimension (its:ite,kts:kte) &
2726 ,intent (inout) :: &
2727 z
2728 real, dimension (its:ite) &
2729 ,intent (in ) :: &
2730 psur,z1
2731 real &
2732 ,intent (in ) :: &
2733 xl,cp
2734 integer, dimension (its:ite) &
2735 ,intent (inout) :: &
2736 ierr
2737 integer &
2738 ,intent (in ) :: &
2739 itest
2740 !
2741 ! local variables in this routine
2742 !
2744 integer :: &
2745 i,k,iph
2746 real, dimension (1:2) :: AE,BE,HT
2747 real, dimension (its:ite,kts:kte) :: tv
2748 real :: tcrit,e,tvbar
2751 HT(1)=XL/CP
2752 HT(2)=2.834E6/CP
2753 BE(1)=.622*HT(1)/.286
2754 AE(1)=BE(1)/273.+ALOG(610.71)
2755 BE(2)=.622*HT(2)/.286
2756 AE(2)=BE(2)/273.+ALOG(610.71)
2757 ! print *, 'TCRIT = ', tcrit,its,ite
2758 DO k=kts,ktf
2759 do i=its,itf
2760 if(ierr(i).eq.0)then
2761 !Csgb - IPH is for phase, dependent on TCRIT (water or ice)
2762 IPH=1
2763 IF(T(I,K).LE.TCRIT)IPH=2
2764 ! print *, 'AE(IPH),BE(IPH) = ',AE(IPH),BE(IPH),AE(IPH)-BE(IPH),T(i,k),i,k
2765 E=EXP(AE(IPH)-BE(IPH)/T(I,K))
2766 ! print *, 'P, E = ', P(I,K), E
2767 QES(I,K)=.622*E/(100.*P(I,K)-E)
2768 IF(QES(I,K).LE.1.E-08)QES(I,K)=1.E-08
2769 IF(QES(I,K).LT.Q(I,K))QES(I,K)=Q(I,K)
2770 ! IF(Q(I,K).GT.QES(I,K))Q(I,K)=QES(I,K)
2771 TV(I,K)=T(I,K)+.608*Q(I,K)*T(I,K)
2772 endif
2773 enddo
2774 enddo
2775 !
2776 !--- z's are calculated with changed h's and q's and t's
2777 !--- if itest=2
2778 !
2779 if(itest.ne.2)then
2780 do i=its,itf
2781 if(ierr(i).eq.0)then
2782 Z(I,1)=max(0.,Z1(I))-(ALOG(P(I,1))- &
2783 ALOG(PSUR(I)))*287.*TV(I,1)/9.81
2784 endif
2785 enddo
2787 ! --- calculate heights
2788 DO K=kts+1,ktf
2789 do i=its,itf
2790 if(ierr(i).eq.0)then
2791 TVBAR=.5*TV(I,K)+.5*TV(I,K-1)
2792 Z(I,K)=Z(I,K-1)-(ALOG(P(I,K))- &
2793 ALOG(P(I,K-1)))*287.*TVBAR/9.81
2794 endif
2795 enddo
2796 enddo
2797 else
2798 do k=kts,ktf
2799 do i=its,itf
2800 if(ierr(i).eq.0)then
2801 z(i,k)=(he(i,k)-1004.*t(i,k)-2.5e6*q(i,k))/9.81
2802 z(i,k)=max(1.e-3,z(i,k))
2803 endif
2804 enddo
2805 enddo
2806 endif
2807 !
2808 !--- calculate moist static energy - HE
2809 ! saturated moist static energy - HES
2810 !
2811 DO k=kts,ktf
2812 do i=its,itf
2813 if(ierr(i).eq.0)then
2814 if(itest.eq.0)HE(I,K)=9.81*Z(I,K)+1004.*T(I,K)+2.5E06*Q(I,K)
2815 HES(I,K)=9.81*Z(I,K)+1004.*T(I,K)+2.5E06*QES(I,K)
2816 IF(HE(I,K).GE.HES(I,K))HE(I,K)=HES(I,K)
2817 endif
2818 enddo
2819 enddo
2821 END SUBROUTINE cup_env
2824 SUBROUTINE cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup, &
2825 he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
2826 ierr,z1,xl,rv,cp, &
2827 itf,jtf,ktf, &
2828 its,ite, jts,jte, kts,kte )
2830 IMPLICIT NONE
2832 integer &
2833 ,intent (in ) :: &
2834 itf,jtf,ktf, &
2835 its,ite, jts,jte, kts,kte
2836 !
2837 ! ierr error value, maybe modified in this routine
2838 ! q = environmental mixing ratio
2839 ! q_cup = environmental mixing ratio on cloud levels
2840 ! qes = environmental saturation mixing ratio
2841 ! qes_cup = environmental saturation mixing ratio on cloud levels
2842 ! t = environmental temp
2843 ! t_cup = environmental temp on cloud levels
2844 ! p = environmental pressure
2845 ! p_cup = environmental pressure on cloud levels
2846 ! z = environmental heights
2847 ! z_cup = environmental heights on cloud levels
2848 ! he = environmental moist static energy
2849 ! he_cup = environmental moist static energy on cloud levels
2850 ! hes = environmental saturation moist static energy
2851 ! hes_cup = environmental saturation moist static energy on cloud levels
2852 ! gamma_cup = gamma on cloud levels
2853 ! psur = surface pressure
2854 ! z1 = terrain elevation
2855 !
2856 !
2857 real, dimension (its:ite,kts:kte) &
2858 ,intent (in ) :: &
2859 qes,q,he,hes,z,p,t
2860 real, dimension (its:ite,kts:kte) &
2861 ,intent (out ) :: &
2862 qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup
2863 real, dimension (its:ite) &
2864 ,intent (in ) :: &
2865 psur,z1
2866 real &
2867 ,intent (in ) :: &
2868 xl,rv,cp
2869 integer, dimension (its:ite) &
2870 ,intent (inout) :: &
2871 ierr
2872 !
2873 ! local variables in this routine
2874 !
2876 integer :: &
2877 i,k
2880 do k=kts+1,ktf
2881 do i=its,itf
2882 if(ierr(i).eq.0)then
2883 qes_cup(i,k)=.5*(qes(i,k-1)+qes(i,k))
2884 q_cup(i,k)=.5*(q(i,k-1)+q(i,k))
2885 hes_cup(i,k)=.5*(hes(i,k-1)+hes(i,k))
2886 he_cup(i,k)=.5*(he(i,k-1)+he(i,k))
2887 if(he_cup(i,k).gt.hes_cup(i,k))he_cup(i,k)=hes_cup(i,k)
2888 z_cup(i,k)=.5*(z(i,k-1)+z(i,k))
2889 p_cup(i,k)=.5*(p(i,k-1)+p(i,k))
2890 t_cup(i,k)=.5*(t(i,k-1)+t(i,k))
2891 gamma_cup(i,k)=(xl/cp)*(xl/(rv*t_cup(i,k) &
2892 *t_cup(i,k)))*qes_cup(i,k)
2893 endif
2894 enddo
2895 enddo
2896 do i=its,itf
2897 if(ierr(i).eq.0)then
2898 qes_cup(i,1)=qes(i,1)
2899 q_cup(i,1)=q(i,1)
2900 hes_cup(i,1)=hes(i,1)
2901 he_cup(i,1)=he(i,1)
2902 z_cup(i,1)=.5*(z(i,1)+z1(i))
2903 p_cup(i,1)=.5*(p(i,1)+psur(i))
2904 t_cup(i,1)=t(i,1)
2905 gamma_cup(i,1)=xl/cp*(xl/(rv*t_cup(i,1) &
2906 *t_cup(i,1)))*qes_cup(i,1)
2907 endif
2908 enddo
2910 END SUBROUTINE cup_env_clev
2913 SUBROUTINE cup_forcing_ens_3d(closure_n,xland,aa0,aa1,xaa0,mbdt,dtime,ierr,ierr2,ierr3,&
2914 xf_ens,j,name,axx,maxens,iens,iedt,maxens2,maxens3,mconv, &
2915 p_cup,ktop,omeg,zd,k22,zu,pr_ens,edt,kbcon,massflx, &
2916 iact_old_gr,dir,ensdim,massfln,icoic,edt_out, &
2917 high_resolution,itf,jtf,ktf, &
2918 its,ite, jts,jte, kts,kte,ens4,ktau )
2920 IMPLICIT NONE
2922 integer &
2923 ,intent (in ) :: &
2924 itf,jtf,ktf, &
2925 its,ite, jts,jte, kts,kte,ens4,high_resolution,ktau
2926 integer, intent (in ) :: &
2927 j,ensdim,maxens,iens,iedt,maxens2,maxens3
2928 !
2929 ! ierr error value, maybe modified in this routine
2930 ! pr_ens = precipitation ensemble
2931 ! xf_ens = mass flux ensembles
2932 ! massfln = downdraft mass flux ensembles used in next timestep
2933 ! omeg = omega from large scale model
2934 ! mconv = moisture convergence from large scale model
2935 ! zd = downdraft normalized mass flux
2936 ! zu = updraft normalized mass flux
2937 ! aa0 = cloud work function without forcing effects
2938 ! aa1 = cloud work function with forcing effects
2939 ! xaa0 = cloud work function with cloud effects (ensemble dependent)
2940 ! edt = epsilon
2941 ! dir = "storm motion"
2942 ! mbdt = arbitrary numerical parameter
2943 ! dtime = dt over which forcing is applied
2944 ! iact_gr_old = flag to tell where convection was active
2945 ! kbcon = LFC of parcel from k22
2946 ! k22 = updraft originating level
2947 ! icoic = flag if only want one closure (usually set to zero!)
2948 ! name = deep or shallow convection flag
2949 !
2950 real, dimension (its:ite,jts:jte,1:ensdim) &
2951 ,intent (inout) :: &
2952 pr_ens
2953 real, dimension (its:ite,jts:jte,1:ensdim) &
2954 ,intent (out ) :: &
2955 xf_ens,massfln
2956 real, dimension (its:ite,jts:jte) &
2957 ,intent (inout ) :: &
2958 edt_out
2959 real, dimension (its:ite,jts:jte) &
2960 ,intent (in ) :: &
2961 massflx
2962 real, dimension (its:ite,kts:kte) &
2963 ,intent (in ) :: &
2964 zd,zu,p_cup
2965 real, dimension (its:ite,kts:kte,1:ens4) &
2966 ,intent (in ) :: &
2967 omeg
2968 real, dimension (its:ite,1:maxens) &
2969 ,intent (in ) :: &
2970 xaa0
2971 real, dimension (its:ite) &
2972 ,intent (in ) :: &
2973 aa1,edt,dir,xland
2974 real, dimension (its:ite,1:ens4) &
2975 ,intent (in ) :: &
2976 mconv,axx
2977 real, dimension (its:ite) &
2978 ,intent (inout) :: &
2979 aa0,closure_n
2980 real, dimension (1:maxens) &
2981 ,intent (in ) :: &
2982 mbdt
2983 real &
2984 ,intent (in ) :: &
2985 dtime
2986 integer, dimension (its:ite,jts:jte) &
2987 ,intent (in ) :: &
2988 iact_old_gr
2989 integer, dimension (its:ite) &
2990 ,intent (in ) :: &
2991 k22,kbcon,ktop
2992 integer, dimension (its:ite) &
2993 ,intent (inout) :: &
2994 ierr,ierr2,ierr3
2995 integer &
2996 ,intent (in ) :: &
2997 icoic
2998 character *(*), intent (in) :: &
2999 name
3000 !
3001 ! local variables in this routine
3002 !
3004 real, dimension (1:maxens3) :: &
3005 xff_ens3
3006 real, dimension (1:maxens) :: &
3007 xk
3008 integer :: &
3009 i,k,nall,n,ne,nens,nens3,iresult,iresultd,iresulte,mkxcrt,kclim
3010 parameter (mkxcrt=15)
3011 real :: &
3012 fens4,a1,massfld,a_ave,xff0,xff00,xxx,xomg,aclim1,aclim2,aclim3,aclim4
3013 real, dimension(1:mkxcrt) :: &
3014 pcrit,acrit,acritt
3016 integer :: nall2,ixxx,irandom
3017 integer, dimension (8) :: seed
3020 DATA PCRIT/850.,800.,750.,700.,650.,600.,550.,500.,450.,400., &
3021 350.,300.,250.,200.,150./
3022 DATA ACRIT/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
3023 .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
3024 ! GDAS DERIVED ACRIT
3025 DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688, &
3026 .743,.813,.886,.947,1.138,1.377,1.896/
3027 !
3028 seed=0
3029 seed(2)=j
3030 seed(3)=ktau
3031 nens=0
3032 irandom=1
3033 if(high_resolution.eq.1)irandom=0
3034 irandom=0
3035 fens4=float(ens4)
3037 !--- LARGE SCALE FORCING
3038 !
3039 DO 100 i=its,itf
3040 if(name.eq.'deeps'.and.ierr(i).gt.995)then
3041 aa0(i)=0.
3042 ierr(i)=0
3043 endif
3044 IF(ierr(i).eq.0)then
3045 !
3046 !---
3047 !
3048 if(name.eq.'deeps')then
3049 !
3050 a_ave=0.
3051 do ne=1,ens4
3052 a_ave=a_ave+axx(i,ne)
3053 enddo
3054 a_ave=max(0.,a_ave/fens4)
3055 a_ave=min(a_ave,aa1(i))
3056 a_ave=max(0.,a_ave)
3057 do ne=1,16
3058 xff_ens3(ne)=0.
3059 enddo
3060 xff0= (AA1(I)-AA0(I))/DTIME
3061 if(high_resolution.eq.1)xff0= (a_ave-AA0(I))/DTIME
3062 xff_ens3(1)=(AA1(I)-AA0(I))/dtime
3063 xff_ens3(2)=(a_ave-AA0(I))/dtime
3064 if(irandom.eq.1)then
3065 seed(1)=i
3066 call random_seed (PUT=seed)
3067 call random_number (xxx)
3068 ixxx=min(ens4,max(1,int(fens4*xxx+1.e-8)))
3069 xff_ens3(3)=(axx(i,ixxx)-AA0(I))/dtime
3070 call random_number (xxx)
3071 ixxx=min(ens4,max(1,int(fens4*xxx+1.e-8)))
3072 xff_ens3(13)=(axx(i,ixxx)-AA0(I))/dtime
3073 else
3074 xff_ens3(3)=(AA1(I)-AA0(I))/dtime
3075 xff_ens3(13)=(AA1(I)-AA0(I))/dtime
3076 endif
3077 if(high_resolution.eq.1)then
3078 xff_ens3(1)=(a_ave-AA0(I))/dtime
3079 xff_ens3(2)=(a_ave-AA0(I))/dtime
3080 xff_ens3(3)=(a_ave-AA0(I))/dtime
3081 xff_ens3(13)=(a_ave-AA0(I))/dtime
3082 endif
3083 !
3084 !--- more original Arakawa-Schubert (climatologic value of aa0)
3085 !
3086 !
3087 !--- omeg is in bar/s, mconv done with omeg in Pa/s
3088 ! more like Brown (1979), or Frank-Cohen (199?)
3089 !
3090 xff_ens3(14)=0.
3091 do ne=1,ens4
3092 xff_ens3(14)=xff_ens3(14)-omeg(i,k22(i),ne)/(fens4*9.81)
3093 enddo
3094 if(xff_ens3(14).lt.0.)xff_ens3(14)=0.
3095 xff_ens3(5)=0.
3096 do ne=1,ens4
3097 xff_ens3(5)=xff_ens3(5)-omeg(i,kbcon(i),ne)/(fens4*9.81)
3098 enddo
3099 if(xff_ens3(5).lt.0.)xff_ens3(5)=0.
3100 !
3101 ! minimum below kbcon
3102 !
3103 if(high_resolution.eq.0)then
3104 xff_ens3(4)=-omeg(i,2,1)/9.81
3105 do k=2,kbcon(i)-1
3106 do ne=1,ens4
3107 xomg=-omeg(i,k,ne)/9.81
3108 if(xomg.lt.xff_ens3(4))xff_ens3(4)=xomg
3109 enddo
3110 enddo
3111 if(xff_ens3(4).lt.0.)xff_ens3(4)=0.
3112 !
3113 ! max below kbcon
3114 xff_ens3(6)=-omeg(i,2,1)/9.81
3115 do k=2,kbcon(i)-1
3116 do ne=1,ens4
3117 xomg=-omeg(i,k,ne)/9.81
3118 if(xomg.gt.xff_ens3(6))xff_ens3(6)=xomg
3119 enddo
3120 enddo
3121 if(xff_ens3(6).lt.0.)xff_ens3(6)=0.
3122 endif
3123 if(high_resolution.eq.1)then
3124 xff_ens3(5)=min(xff_ens3(5),xff_ens3(14))
3125 xff_ens3(4)=xff_ens3(5)
3126 xff_ens3(6)=xff_ens3(5)
3127 endif
3128 !
3129 !--- more like Krishnamurti et al.; pick max and average values
3130 !
3131 xff_ens3(7)=mconv(i,1)
3132 xff_ens3(8)=mconv(i,1)
3133 xff_ens3(9)=mconv(i,1)
3134 if(ens4.gt.1)then
3135 do ne=2,ens4
3136 if (mconv(i,ne).gt.xff_ens3(7))xff_ens3(7)=mconv(i,ne)
3137 enddo
3138 do ne=2,ens4
3139 if (mconv(i,ne).lt.xff_ens3(8))xff_ens3(8)=mconv(i,ne)
3140 enddo
3141 do ne=2,ens4
3142 xff_ens3(9)=xff_ens3(9)+mconv(i,ne)
3143 enddo
3144 xff_ens3(9)=xff_ens3(9)/fens4
3145 endif
3146 if(high_resolution.eq.1)then
3147 xff_ens3(7)=xff_ens3(9)
3148 xff_ens3(8)=xff_ens3(9)
3149 xff_ens3(15)=xff_ens3(9)
3150 endif
3151 !
3152 if(high_resolution.eq.0)then
3153 if(irandom.eq.1)then
3154 seed(1)=i
3155 call random_seed (PUT=seed)
3156 call random_number (xxx)
3157 ixxx=min(ens4,max(1,int(fens4*xxx+1.e-8)))
3158 xff_ens3(15)=mconv(i,ixxx)
3159 else
3160 xff_ens3(15)=mconv(i,1)
3161 endif
3162 endif
3163 !
3164 !--- more like Fritsch Chappel or Kain Fritsch (plus triggers)
3165 !
3166 xff_ens3(10)=A_AVE/(60.*40.)
3167 xff_ens3(11)=AA1(I)/(60.*40.)
3168 if(irandom.eq.1)then
3169 seed(1)=i
3170 call random_seed (PUT=seed)
3171 call random_number (xxx)
3172 ixxx=min(ens4,max(1,int(fens4*xxx+1.e-8)))
3173 xff_ens3(12)=AXX(I,ixxx)/(60.*40.)
3174 else
3175 xff_ens3(12)=AA1(I)/(60.*40.)
3176 endif
3177 if(high_resolution.eq.1)then
3178 xff_ens3(11)=xff_ens3(10)
3179 xff_ens3(12)=xff_ens3(10)
3180 endif
3181 !
3182 !--- more original Arakawa-Schubert (climatologic value of aa0)
3183 !
3184 ! edt_out(i,j)=xff0
3185 if(icoic.eq.0)then
3186 if(xff0.lt.0.)then
3187 xff_ens3(1)=0.
3188 xff_ens3(2)=0.
3189 xff_ens3(3)=0.
3190 xff_ens3(13)=0.
3191 xff_ens3(10)=0.
3192 xff_ens3(11)=0.
3193 xff_ens3(12)=0.
3194 endif
3195 endif
3199 do nens=1,maxens
3200 XK(nens)=(XAA0(I,nens)-AA1(I))/MBDT(2)
3201 if(xk(nens).le.0.and.xk(nens).gt.-1.e-6) &
3202 xk(nens)=-1.e-6
3203 if(xk(nens).gt.0.and.xk(nens).lt.1.e-6) &
3204 xk(nens)=1.e-6
3205 enddo
3206 !
3207 !--- add up all ensembles
3208 !
3209 do 350 ne=1,maxens
3210 !
3211 !--- for every xk, we have maxens3 xffs
3212 !--- iens is from outermost ensemble (most expensive!
3213 !
3214 !--- iedt (maxens2 belongs to it)
3215 !--- is from second, next outermost, not so expensive
3216 !
3217 !--- so, for every outermost loop, we have maxens*maxens2*3
3218 !--- ensembles!!! nall would be 0, if everything is on first
3219 !--- loop index, then ne would start counting, then iedt, then iens....
3220 !
3221 iresult=0
3222 iresultd=0
3223 iresulte=0
3224 nall=(iens-1)*maxens3*maxens*maxens2 &
3225 +(iedt-1)*maxens*maxens3 &
3226 +(ne-1)*maxens3
3227 !
3228 ! over water, enfor!e small cap for some of the closures
3229 !
3230 if(xland(i).lt.0.1)then
3231 if(ierr2(i).gt.0.or.ierr3(i).gt.0)then
3232 xff_ens3(1) =0.
3233 massfln(i,j,nall+1)=0.
3234 xff_ens3(2) =0.
3235 massfln(i,j,nall+2)=0.
3236 xff_ens3(3) =0.
3237 massfln(i,j,nall+3)=0.
3238 xff_ens3(10) =0.
3239 massfln(i,j,nall+10)=0.
3240 xff_ens3(11) =0.
3241 massfln(i,j,nall+11)=0.
3242 xff_ens3(12) =0.
3243 massfln(i,j,nall+12)=0.
3244 xff_ens3(7) =0.
3245 massfln(i,j,nall+7)=0.
3246 xff_ens3(8) =0.
3247 massfln(i,j,nall+8)=0.
3248 xff_ens3(9) =0.
3249 massfln(i,j,nall+9)=0.
3250 closure_n(i)=closure_n(i)-1.
3251 xff_ens3(13) =0.
3252 massfln(i,j,nall+13)=0.
3253 xff_ens3(15) =0.
3254 massfln(i,j,nall+15)=0.
3255 endif
3256 endif
3257 !
3258 ! end water treatment
3259 !
3260 !
3261 !--- check for upwind convection
3262 ! iresult=0
3263 massfld=0.
3265 ! call cup_direction2(i,j,dir,iact_old_gr, &
3266 ! massflx,iresult,1, &
3267 ! massfld, &
3268 ! itf,jtf,ktf, &
3269 ! ims,ime, jms,jme, kms,kme, &
3270 ! its,ite, jts,jte, kts,kte )
3271 ! if(i.eq.ipr.and.j.eq.jpr.and.iedt.eq.1.and.ne.eq.1)then
3272 ! if(iedt.eq.1.and.ne.eq.1)then
3273 ! print *,massfld,ne,iedt,iens
3274 ! print *,xk(ne),xff_ens3(1),xff_ens3(2),xff_ens3(3)
3275 ! endif
3276 ! print *,i,j,massfld,aa0(i),aa1(i)
3277 IF(XK(ne).lt.0.and.xff0.gt.0.)iresultd=1
3278 iresulte=max(iresult,iresultd)
3279 iresulte=1
3280 if(iresulte.eq.1)then
3281 !
3282 !--- special treatment for stability closures
3283 !
3285 if(xff0.ge.0.)then
3286 xf_ens(i,j,nall+1)=massfld
3287 xf_ens(i,j,nall+2)=massfld
3288 xf_ens(i,j,nall+3)=massfld
3289 xf_ens(i,j,nall+13)=massfld
3290 if(xff_ens3(1).gt.0)xf_ens(i,j,nall+1)=max(0.,-xff_ens3(1)/xk(ne)) &
3291 +massfld
3292 if(xff_ens3(2).gt.0)xf_ens(i,j,nall+2)=max(0.,-xff_ens3(2)/xk(ne)) &
3293 +massfld
3294 if(xff_ens3(3).gt.0)xf_ens(i,j,nall+3)=max(0.,-xff_ens3(3)/xk(ne)) &
3295 +massfld
3296 if(xff_ens3(13).gt.0)xf_ens(i,j,nall+13)=max(0.,-xff_ens3(13)/xk(ne)) &
3297 +massfld
3298 ! endif
3299 else
3300 xf_ens(i,j,nall+1)=massfld
3301 xf_ens(i,j,nall+2)=massfld
3302 xf_ens(i,j,nall+3)=massfld
3303 xf_ens(i,j,nall+13)=massfld
3304 endif
3305 !
3306 !--- if iresult.eq.1, following independent of xff0
3307 !
3308 xf_ens(i,j,nall+4)=max(0.,xff_ens3(4) &
3309 +massfld)
3310 xf_ens(i,j,nall+5)=max(0.,xff_ens3(5) &
3311 +massfld)
3312 xf_ens(i,j,nall+6)=max(0.,xff_ens3(6) &
3313 +massfld)
3314 xf_ens(i,j,nall+14)=max(0.,xff_ens3(14) &
3315 +massfld)
3316 a1=max(1.e-3,pr_ens(i,j,nall+7))
3317 xf_ens(i,j,nall+7)=max(0.,xff_ens3(7) &
3318 /a1)
3319 a1=max(1.e-3,pr_ens(i,j,nall+8))
3320 xf_ens(i,j,nall+8)=max(0.,xff_ens3(8) &
3321 /a1)
3322 a1=max(1.e-3,pr_ens(i,j,nall+9))
3323 xf_ens(i,j,nall+9)=max(0.,xff_ens3(9) &
3324 /a1)
3325 a1=max(1.e-3,pr_ens(i,j,nall+15))
3326 xf_ens(i,j,nall+15)=max(0.,xff_ens3(15) &
3327 /a1)
3328 if(XK(ne).lt.0.)then
3329 xf_ens(i,j,nall+10)=max(0., &
3330 -xff_ens3(10)/xk(ne)) &
3331 +massfld
3332 xf_ens(i,j,nall+11)=max(0., &
3333 -xff_ens3(11)/xk(ne)) &
3334 +massfld
3335 xf_ens(i,j,nall+12)=max(0., &
3336 -xff_ens3(12)/xk(ne)) &
3337 +massfld
3338 else
3339 xf_ens(i,j,nall+10)=massfld
3340 xf_ens(i,j,nall+11)=massfld
3341 xf_ens(i,j,nall+12)=massfld
3342 endif
3343 if(icoic.ge.1)then
3344 closure_n(i)=0.
3345 xf_ens(i,j,nall+1)=xf_ens(i,j,nall+icoic)
3346 xf_ens(i,j,nall+2)=xf_ens(i,j,nall+icoic)
3347 xf_ens(i,j,nall+3)=xf_ens(i,j,nall+icoic)
3348 xf_ens(i,j,nall+4)=xf_ens(i,j,nall+icoic)
3349 xf_ens(i,j,nall+5)=xf_ens(i,j,nall+icoic)
3350 xf_ens(i,j,nall+6)=xf_ens(i,j,nall+icoic)
3351 xf_ens(i,j,nall+7)=xf_ens(i,j,nall+icoic)
3352 xf_ens(i,j,nall+8)=xf_ens(i,j,nall+icoic)
3353 xf_ens(i,j,nall+9)=xf_ens(i,j,nall+icoic)
3354 xf_ens(i,j,nall+10)=xf_ens(i,j,nall+icoic)
3355 xf_ens(i,j,nall+11)=xf_ens(i,j,nall+icoic)
3356 xf_ens(i,j,nall+12)=xf_ens(i,j,nall+icoic)
3357 xf_ens(i,j,nall+13)=xf_ens(i,j,nall+icoic)
3358 xf_ens(i,j,nall+14)=xf_ens(i,j,nall+icoic)
3359 xf_ens(i,j,nall+15)=xf_ens(i,j,nall+icoic)
3360 xf_ens(i,j,nall+16)=xf_ens(i,j,nall+icoic)
3361 endif
3362 !
3363 ! 16 is a randon pick from the oher 15
3364 !
3365 if(irandom.eq.1)then
3366 call random_number (xxx)
3367 ixxx=min(15,max(1,int(15.*xxx+1.e-8)))
3368 xf_ens(i,j,nall+16)=xf_ens(i,j,nall+ixxx)
3369 else
3370 xf_ens(i,j,nall+16)=xf_ens(i,j,nall+1)
3371 endif
3372 !
3373 !
3374 !--- store new for next time step
3375 !
3376 do nens3=1,maxens3
3377 massfln(i,j,nall+nens3)=edt(i) &
3378 *xf_ens(i,j,nall+nens3)
3379 massfln(i,j,nall+nens3)=max(0., &
3380 massfln(i,j,nall+nens3))
3381 enddo
3382 !
3383 !
3384 !--- do some more on the caps!!! ne=1 for 175, ne=2 for 100,....
3385 !
3386 ! do not care for caps here for closure groups 1 and 5,
3387 ! they are fine, do not turn them off here
3388 !
3389 !
3390 if(ne.eq.2.and.ierr2(i).gt.0)then
3391 xf_ens(i,j,nall+1) =0.
3392 xf_ens(i,j,nall+2) =0.
3393 xf_ens(i,j,nall+3) =0.
3394 xf_ens(i,j,nall+4) =0.
3395 xf_ens(i,j,nall+5) =0.
3396 xf_ens(i,j,nall+6) =0.
3397 xf_ens(i,j,nall+7) =0.
3398 xf_ens(i,j,nall+8) =0.
3399 xf_ens(i,j,nall+9) =0.
3400 xf_ens(i,j,nall+10)=0.
3401 xf_ens(i,j,nall+11)=0.
3402 xf_ens(i,j,nall+12)=0.
3403 xf_ens(i,j,nall+13)=0.
3404 xf_ens(i,j,nall+14)=0.
3405 xf_ens(i,j,nall+15)=0.
3406 xf_ens(i,j,nall+16)=0.
3407 massfln(i,j,nall+1)=0.
3408 massfln(i,j,nall+2)=0.
3409 massfln(i,j,nall+3)=0.
3410 massfln(i,j,nall+4)=0.
3411 massfln(i,j,nall+5)=0.
3412 massfln(i,j,nall+6)=0.
3413 massfln(i,j,nall+7)=0.
3414 massfln(i,j,nall+8)=0.
3415 massfln(i,j,nall+9)=0.
3416 massfln(i,j,nall+10)=0.
3417 massfln(i,j,nall+11)=0.
3418 massfln(i,j,nall+12)=0.
3419 massfln(i,j,nall+13)=0.
3420 massfln(i,j,nall+14)=0.
3421 massfln(i,j,nall+15)=0.
3422 massfln(i,j,nall+16)=0.
3423 endif
3424 if(ne.eq.3.and.ierr3(i).gt.0)then
3425 xf_ens(i,j,nall+1) =0.
3426 xf_ens(i,j,nall+2) =0.
3427 xf_ens(i,j,nall+3) =0.
3428 xf_ens(i,j,nall+4) =0.
3429 xf_ens(i,j,nall+5) =0.
3430 xf_ens(i,j,nall+6) =0.
3431 xf_ens(i,j,nall+7) =0.
3432 xf_ens(i,j,nall+8) =0.
3433 xf_ens(i,j,nall+9) =0.
3434 xf_ens(i,j,nall+10)=0.
3435 xf_ens(i,j,nall+11)=0.
3436 xf_ens(i,j,nall+12)=0.
3437 xf_ens(i,j,nall+13)=0.
3438 xf_ens(i,j,nall+14)=0.
3439 xf_ens(i,j,nall+15)=0.
3440 xf_ens(i,j,nall+16)=0.
3441 massfln(i,j,nall+1)=0.
3442 massfln(i,j,nall+2)=0.
3443 massfln(i,j,nall+3)=0.
3444 massfln(i,j,nall+4)=0.
3445 massfln(i,j,nall+5)=0.
3446 massfln(i,j,nall+6)=0.
3447 massfln(i,j,nall+7)=0.
3448 massfln(i,j,nall+8)=0.
3449 massfln(i,j,nall+9)=0.
3450 massfln(i,j,nall+10)=0.
3451 massfln(i,j,nall+11)=0.
3452 massfln(i,j,nall+12)=0.
3453 massfln(i,j,nall+13)=0.
3454 massfln(i,j,nall+14)=0.
3455 massfln(i,j,nall+15)=0.
3456 massfln(i,j,nall+16)=0.
3457 endif
3459 endif
3460 350 continue
3461 ! ne=1, cap=175
3462 !
3463 nall=(iens-1)*maxens3*maxens*maxens2 &
3464 +(iedt-1)*maxens*maxens3
3465 ! ne=2, cap=100
3466 !
3467 nall2=(iens-1)*maxens3*maxens*maxens2 &
3468 +(iedt-1)*maxens*maxens3 &
3469 +(2-1)*maxens3
3470 xf_ens(i,j,nall+4) = xf_ens(i,j,nall2+4)
3471 xf_ens(i,j,nall+5) =xf_ens(i,j,nall2+5)
3472 xf_ens(i,j,nall+6) =xf_ens(i,j,nall2+6)
3473 xf_ens(i,j,nall+14) =xf_ens(i,j,nall2+14)
3474 xf_ens(i,j,nall+7) =xf_ens(i,j,nall2+7)
3475 xf_ens(i,j,nall+8) =xf_ens(i,j,nall2+8)
3476 xf_ens(i,j,nall+9) =xf_ens(i,j,nall2+9)
3477 xf_ens(i,j,nall+15) =xf_ens(i,j,nall2+15)
3478 xf_ens(i,j,nall+10)=xf_ens(i,j,nall2+10)
3479 xf_ens(i,j,nall+11)=xf_ens(i,j,nall2+11)
3480 xf_ens(i,j,nall+12)=xf_ens(i,j,nall2+12)
3481 go to 100
3482 endif
3483 elseif(ierr(i).ne.20.and.ierr(i).ne.0)then
3484 do n=1,ensdim
3485 xf_ens(i,j,n)=0.
3486 massfln(i,j,n)=0.
3487 enddo
3488 endif
3489 100 continue
3491 END SUBROUTINE cup_forcing_ens_3d
3494 SUBROUTINE cup_kbcon(cap_inc,iloop,k22,kbcon,he_cup,hes_cup, &
3495 ierr,kbmax,p_cup,cap_max, &
3496 itf,jtf,ktf, &
3497 its,ite, jts,jte, kts,kte )
3499 IMPLICIT NONE
3500 !
3502 ! only local wrf dimensions are need as of now in this routine
3504 integer &
3505 ,intent (in ) :: &
3506 itf,jtf,ktf, &
3507 its,ite, jts,jte, kts,kte
3508 !
3509 !
3510 !
3511 ! ierr error value, maybe modified in this routine
3512 !
3513 real, dimension (its:ite,kts:kte) &
3514 ,intent (in ) :: &
3515 he_cup,hes_cup,p_cup
3516 real, dimension (its:ite) &
3517 ,intent (in ) :: &
3518 cap_max,cap_inc
3519 integer, dimension (its:ite) &
3520 ,intent (in ) :: &
3521 kbmax
3522 integer, dimension (its:ite) &
3523 ,intent (inout) :: &
3524 kbcon,k22,ierr
3525 integer &
3526 ,intent (in ) :: &
3527 iloop
3528 !
3529 ! local variables in this routine
3530 !
3532 integer :: &
3533 i
3534 real :: &
3535 pbcdif,plus
3536 !
3537 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
3538 !
3539 DO 27 i=its,itf
3540 kbcon(i)=1
3541 IF(ierr(I).ne.0)GO TO 27
3542 KBCON(I)=K22(I)
3543 GO TO 32
3544 31 CONTINUE
3545 KBCON(I)=KBCON(I)+1
3546 IF(KBCON(I).GT.KBMAX(i)+2)THEN
3547 if(iloop.ne.4)ierr(i)=3
3548 ! if(iloop.lt.4)ierr(i)=997
3549 GO TO 27
3550 ENDIF
3551 32 CONTINUE
3552 IF(HE_cup(I,K22(I)).LT.HES_cup(I,KBCON(I)))GO TO 31
3554 ! cloud base pressure and max moist static energy pressure
3555 ! i.e., the depth (in mb) of the layer of negative buoyancy
3556 if(KBCON(I)-K22(I).eq.1)go to 27
3557 PBCDIF=-P_cup(I,KBCON(I))+P_cup(I,K22(I))
3558 plus=max(25.,cap_max(i)-float(iloop-1)*cap_inc(i))
3559 if(iloop.eq.4)plus=cap_max(i)
3560 !
3561 ! for shallow convection, if cap_max is greater than 25, it is the pressure at pbltop
3562 if(iloop.eq.5)plus=25.
3563 if(iloop.eq.5.and.cap_max(i).gt.25)pbcdif=-P_cup(I,KBCON(I))+cap_max(i)
3564 IF(PBCDIF.GT.plus)THEN
3565 K22(I)=K22(I)+1
3566 KBCON(I)=K22(I)
3567 GO TO 32
3568 ENDIF
3569 27 CONTINUE
3571 END SUBROUTINE cup_kbcon
3574 SUBROUTINE cup_ktop(ilo,dby,kbcon,ktop,ierr, &
3575 itf,jtf,ktf, &
3576 its,ite, jts,jte, kts,kte )
3578 IMPLICIT NONE
3579 !
3580 ! on input
3581 !
3583 ! only local wrf dimensions are need as of now in this routine
3585 integer &
3586 ,intent (in ) :: &
3587 itf,jtf,ktf, &
3588 its,ite, jts,jte, kts,kte
3589 ! dby = buoancy term
3590 ! ktop = cloud top (output)
3591 ! ilo = flag
3592 ! ierr error value, maybe modified in this routine
3593 !
3594 real, dimension (its:ite,kts:kte) &
3595 ,intent (inout) :: &
3596 dby
3597 integer, dimension (its:ite) &
3598 ,intent (in ) :: &
3599 kbcon
3600 integer &
3601 ,intent (in ) :: &
3602 ilo
3603 integer, dimension (its:ite) &
3604 ,intent (out ) :: &
3605 ktop
3606 integer, dimension (its:ite) &
3607 ,intent (inout) :: &
3608 ierr
3609 !
3610 ! local variables in this routine
3611 !
3613 integer :: &
3614 i,k
3615 !
3616 DO 42 i=its,itf
3617 ktop(i)=1
3618 IF(ierr(I).EQ.0)then
3619 DO 40 K=KBCON(I)+1,ktf-1
3620 IF(DBY(I,K).LE.0.)THEN
3621 KTOP(I)=K-1
3622 GO TO 41
3623 ENDIF
3624 40 CONTINUE
3625 if(ilo.eq.1)ierr(i)=5
3626 ! if(ilo.eq.2)ierr(i)=998
3627 GO TO 42
3628 41 CONTINUE
3629 do k=ktop(i)+1,ktf
3630 dby(i,k)=0.
3631 enddo
3632 endif
3633 42 CONTINUE
3635 END SUBROUTINE cup_ktop
3638 SUBROUTINE cup_MAXIMI(ARRAY,KS,KE,MAXX,ierr, &
3639 itf,jtf,ktf, &
3640 its,ite, jts,jte, kts,kte )
3642 IMPLICIT NONE
3643 !
3644 ! on input
3645 !
3647 ! only local wrf dimensions are need as of now in this routine
3649 integer &
3650 ,intent (in ) :: &
3651 itf,jtf,ktf, &
3652 its,ite, jts,jte, kts,kte
3653 ! array input array
3654 ! x output array with return values
3655 ! kt output array of levels
3656 ! ks,kend check-range
3657 real, dimension (its:ite,kts:kte) &
3658 ,intent (in ) :: &
3659 array
3660 integer, dimension (its:ite) &
3661 ,intent (in ) :: &
3662 ierr,ke
3663 integer &
3664 ,intent (in ) :: &
3665 ks
3666 integer, dimension (its:ite) &
3667 ,intent (out ) :: &
3668 maxx
3669 real, dimension (its:ite) :: &
3670 x
3671 real :: &
3672 xar
3673 integer :: &
3674 i,k
3676 DO 200 i=its,itf
3677 MAXX(I)=KS
3678 if(ierr(i).eq.0)then
3679 X(I)=ARRAY(I,KS)
3680 !
3681 DO 100 K=KS,KE(i)
3682 XAR=ARRAY(I,K)
3683 IF(XAR.GE.X(I)) THEN
3684 X(I)=XAR
3685 MAXX(I)=K
3686 ENDIF
3687 100 CONTINUE
3688 endif
3689 200 CONTINUE
3691 END SUBROUTINE cup_MAXIMI
3694 SUBROUTINE cup_minimi(ARRAY,KS,KEND,KT,ierr, &
3695 itf,jtf,ktf, &
3696 its,ite, jts,jte, kts,kte )
3698 IMPLICIT NONE
3699 !
3700 ! on input
3701 !
3703 ! only local wrf dimensions are need as of now in this routine
3705 integer &
3706 ,intent (in ) :: &
3707 itf,jtf,ktf, &
3708 its,ite, jts,jte, kts,kte
3709 ! array input array
3710 ! x output array with return values
3711 ! kt output array of levels
3712 ! ks,kend check-range
3713 real, dimension (its:ite,kts:kte) &
3714 ,intent (in ) :: &
3715 array
3716 integer, dimension (its:ite) &
3717 ,intent (in ) :: &
3718 ierr,ks,kend
3719 integer, dimension (its:ite) &
3720 ,intent (out ) :: &
3721 kt
3722 real, dimension (its:ite) :: &
3723 x
3724 integer :: &
3725 i,k,kstop
3727 DO 200 i=its,itf
3728 KT(I)=KS(I)
3729 if(ierr(i).eq.0)then
3730 X(I)=ARRAY(I,KS(I))
3731 KSTOP=MAX(KS(I)+1,KEND(I))
3732 !
3733 DO 100 K=KS(I)+1,KSTOP
3734 IF(ARRAY(I,K).LT.X(I)) THEN
3735 X(I)=ARRAY(I,K)
3736 KT(I)=K
3737 ENDIF
3738 100 CONTINUE
3739 endif
3740 200 CONTINUE
3742 END SUBROUTINE cup_MINIMI
3745 SUBROUTINE cup_output_ens_3d(xf_ens,ierr,dellat,dellaq,dellaqc, &
3746 subt_ens,subq_ens,subt,subq,outtem,outq,outqc, &
3747 zu,sub_mas,pre,pw,xmb,ktop, &
3748 j,name,nx,nx2,iens,ierr2,ierr3,pr_ens, &
3749 maxens3,ensdim,massfln, &
3750 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
3751 APR_CAPMA,APR_CAPME,APR_CAPMI,closure_n,xland1, &
3752 itf,jtf,ktf, &
3753 its,ite, jts,jte, kts,kte)
3755 IMPLICIT NONE
3756 !
3757 ! on input
3758 !
3760 ! only local wrf dimensions are need as of now in this routine
3762 integer &
3763 ,intent (in ) :: &
3764 itf,jtf,ktf, &
3765 its,ite, jts,jte, kts,kte
3766 integer, intent (in ) :: &
3767 j,ensdim,nx,nx2,iens,maxens3
3768 ! xf_ens = ensemble mass fluxes
3769 ! pr_ens = precipitation ensembles
3770 ! dellat = change of temperature per unit mass flux of cloud ensemble
3771 ! dellaq = change of q per unit mass flux of cloud ensemble
3772 ! dellaqc = change of qc per unit mass flux of cloud ensemble
3773 ! outtem = output temp tendency (per s)
3774 ! outq = output q tendency (per s)
3775 ! outqc = output qc tendency (per s)
3776 ! pre = output precip
3777 ! xmb = total base mass flux
3778 ! xfac1 = correction factor
3779 ! pw = pw -epsilon*pd (ensemble dependent)
3780 ! ierr error value, maybe modified in this routine
3781 !
3782 real, dimension (its:ite,jts:jte,1:ensdim) &
3783 ,intent (inout) :: &
3784 xf_ens,pr_ens,massfln
3785 real, dimension (its:ite,jts:jte) &
3786 ,intent (inout) :: &
3787 APR_GR,APR_W,APR_MC,APR_ST,APR_AS,APR_CAPMA, &
3788 APR_CAPME,APR_CAPMI
3790 real, dimension (its:ite,kts:kte) &
3791 ,intent (out ) :: &
3792 outtem,outq,outqc,subt,subq,sub_mas
3793 real, dimension (its:ite,kts:kte) &
3794 ,intent (in ) :: &
3795 zu
3796 real, dimension (its:ite) &
3797 ,intent (out ) :: &
3798 pre,xmb
3799 real, dimension (its:ite) &
3800 ,intent (inout ) :: &
3801 closure_n,xland1
3802 real, dimension (its:ite,kts:kte,1:nx) &
3803 ,intent (in ) :: &
3804 subt_ens,subq_ens,dellat,dellaqc,dellaq,pw
3805 integer, dimension (its:ite) &
3806 ,intent (in ) :: &
3807 ktop
3808 integer, dimension (its:ite) &
3809 ,intent (inout) :: &
3810 ierr,ierr2,ierr3
3811 !
3812 ! local variables in this routine
3813 !
3815 integer :: &
3816 i,k,n,ncount
3817 real :: &
3818 outtes,ddtes,dtt,dtq,dtqc,dtpw,tuning,prerate,clos_wei,xmbhelp
3819 real :: &
3820 dtts,dtqs
3821 real, dimension (its:ite) :: &
3822 xfac1,xfac2
3823 real, dimension (its:ite):: &
3824 xmb_ske,xmb_ave,xmb_std,xmb_cur,xmbweight
3825 real, dimension (its:ite):: &
3826 pr_ske,pr_ave,pr_std,pr_cur
3827 real, dimension (its:ite,jts:jte):: &
3828 pr_gr,pr_w,pr_mc,pr_st,pr_as,pr_capma, &
3829 pr_capme,pr_capmi
3830 real, dimension (5) :: weight,wm,wm1,wm2,wm3
3831 real, dimension (its:ite,5) :: xmb_w
3833 !
3834 character *(*), intent (in) :: &
3835 name
3836 !
3837 weight(1) = -999. !this will turn off weights
3838 wm(1)=-999.
3840 tuning=0.
3841 !
3842 !
3843 DO k=kts,ktf
3844 do i=its,itf
3845 outtem(i,k)=0.
3846 outq(i,k)=0.
3847 outqc(i,k)=0.
3848 subt(i,k)=0.
3849 subq(i,k)=0.
3850 sub_mas(i,k)=0.
3851 enddo
3852 enddo
3853 do i=its,itf
3854 pre(i)=0.
3855 xmb(i)=0.
3856 xfac1(i)=0.
3857 xfac2(i)=0.
3858 xmbweight(i)=1.
3859 enddo
3860 do i=its,itf
3861 IF(ierr(i).eq.0)then
3862 do n=(iens-1)*nx*nx2*maxens3+1,iens*nx*nx2*maxens3
3863 if(pr_ens(i,j,n).le.0.)then
3864 xf_ens(i,j,n)=0.
3865 endif
3866 enddo
3867 endif
3868 enddo
3869 !
3870 !--- calculate ensemble average mass fluxes
3871 !
3872 call massflx_stats(xf_ens,ensdim,nx2,nx,maxens3, &
3873 xmb_ave,xmb_std,xmb_cur,xmb_ske,j,ierr,1, &
3874 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
3875 APR_CAPMA,APR_CAPME,APR_CAPMI, &
3876 pr_gr,pr_w,pr_mc,pr_st,pr_as, &
3877 pr_capma,pr_capme,pr_capmi, &
3878 itf,jtf,ktf, &
3879 its,ite, jts,jte, kts,kte )
3880 xmb_w=0.
3881 call massflx_stats(pr_ens,ensdim,nx2,nx,maxens3, &
3882 pr_ave,pr_std,pr_cur,pr_ske,j,ierr,2, &
3883 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
3884 APR_CAPMA,APR_CAPME,APR_CAPMI, &
3885 pr_gr,pr_w,pr_mc,pr_st,pr_as, &
3886 pr_capma,pr_capme,pr_capmi, &
3887 itf,jtf,ktf, &
3888 its,ite, jts,jte, kts,kte )
3889 !
3890 !-- now do feedback
3891 !
3892 ddtes=100.
3893 do i=its,itf
3894 if(ierr(i).eq.0)then
3895 if(xmb_ave(i).le.0.)then
3896 ierr(i)=13
3897 xmb_ave(i)=0.
3898 endif
3899 xmb(i)=max(.1*xmb_ave(i),xmb_ave(i)-tuning*xmb_std(i))
3900 ! --- Now use proper count of how many closures were actually
3901 ! used in cup_forcing_ens (including screening of some
3902 ! closures over water) to properly normalize xmb
3903 clos_wei=16./max(1.,closure_n(i))
3904 if (xland1(i).lt.0.5)xmb(i)=xmb(i)*clos_wei
3905 if(xmb(i).eq.0.)then
3906 ierr(i)=19
3907 endif
3908 if(xmb(i).gt.100.)then
3909 ierr(i)=19
3910 endif
3911 xfac1(i)=xmb(i)
3912 xfac2(i)=xmb(i)
3914 endif
3915 ! if(weight(1).lt.-100.)xfac1(i)=xmb_ave(i)
3916 ! if(weight(1).lt.-100.)xfac2(i)=xmb_ave(i)
3917 ENDDO
3918 DO k=kts,ktf
3919 do i=its,itf
3920 dtt=0.
3921 dtts=0.
3922 dtq=0.
3923 dtqs=0.
3924 dtqc=0.
3925 dtpw=0.
3926 IF(ierr(i).eq.0.and.k.le.ktop(i))then
3927 do n=1,nx
3928 dtt=dtt+dellat(i,k,n)
3929 dtts=dtts+subt_ens(i,k,n)
3930 dtq=dtq+dellaq(i,k,n)
3931 dtqs=dtqs+subq_ens(i,k,n)
3932 dtqc=dtqc+dellaqc(i,k,n)
3933 dtpw=dtpw+pw(i,k,n)
3934 enddo
3935 OUTTEM(I,K)=XMB(I)*dtt/float(nx)
3936 SUBT(I,K)=XMB(I)*dtts/float(nx)
3937 OUTQ(I,K)=XMB(I)*dtq/float(nx)
3938 SUBQ(I,K)=XMB(I)*dtqs/float(nx)
3939 OUTQC(I,K)=XMB(I)*dtqc/float(nx)
3940 PRE(I)=PRE(I)+XMB(I)*dtpw/float(nx)
3941 sub_mas(i,k)=zu(i,k)*xmb(i)
3942 endif
3943 enddo
3944 enddo
3946 do i=its,itf
3947 if(ierr(i).eq.0)then
3948 do k=(iens-1)*nx*nx2*maxens3+1,iens*nx*nx2*maxens3
3949 massfln(i,j,k)=massfln(i,j,k)*xfac1(i)
3950 xf_ens(i,j,k)=xf_ens(i,j,k)*xfac1(i)
3951 enddo
3952 endif
3953 ENDDO
3955 END SUBROUTINE cup_output_ens_3d
3958 SUBROUTINE cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
3959 kbcon,ktop,ierr, &
3960 itf,jtf,ktf, &
3961 its,ite, jts,jte, kts,kte )
3963 IMPLICIT NONE
3964 !
3965 ! on input
3966 !
3968 ! only local wrf dimensions are need as of now in this routine
3970 integer &
3971 ,intent (in ) :: &
3972 itf,jtf,ktf, &
3973 its,ite, jts,jte, kts,kte
3974 ! aa0 cloud work function
3975 ! gamma_cup = gamma on model cloud levels
3976 ! t_cup = temperature (Kelvin) on model cloud levels
3977 ! dby = buoancy term
3978 ! zu= normalized updraft mass flux
3979 ! z = heights of model levels
3980 ! ierr error value, maybe modified in this routine
3981 !
3982 real, dimension (its:ite,kts:kte) &
3983 ,intent (in ) :: &
3984 z,zu,gamma_cup,t_cup,dby
3985 integer, dimension (its:ite) &
3986 ,intent (in ) :: &
3987 kbcon,ktop
3988 !
3989 ! input and output
3990 !
3993 integer, dimension (its:ite) &
3994 ,intent (inout) :: &
3995 ierr
3996 real, dimension (its:ite) &
3997 ,intent (out ) :: &
3998 aa0
3999 !
4000 ! local variables in this routine
4001 !
4003 integer :: &
4004 i,k
4005 real :: &
4006 dz,da
4007 !
4008 do i=its,itf
4009 aa0(i)=0.
4010 enddo
4011 DO 100 k=kts+1,ktf
4012 DO 100 i=its,itf
4013 IF(ierr(i).ne.0)GO TO 100
4014 IF(K.LE.KBCON(I))GO TO 100
4015 IF(K.Gt.KTOP(I))GO TO 100
4016 DZ=Z(I,K)-Z(I,K-1)
4017 da=zu(i,k)*DZ*(9.81/(1004.*( &
4018 (T_cup(I,K)))))*DBY(I,K-1)/ &
4019 (1.+GAMMA_CUP(I,K))
4020 IF(K.eq.KTOP(I).and.da.le.0.)go to 100
4021 AA0(I)=AA0(I)+da
4022 if(aa0(i).lt.0.)aa0(i)=0.
4023 100 continue
4025 END SUBROUTINE cup_up_aa0
4028 SUBROUTINE cup_up_he(k22,hkb,z_cup,cd,entr,he_cup,hc, &
4029 kbcon,ierr,dby,he,hes_cup, &
4030 itf,jtf,ktf, &
4031 its,ite, jts,jte, kts,kte )
4033 IMPLICIT NONE
4034 !
4035 ! on input
4036 !
4038 ! only local wrf dimensions are need as of now in this routine
4040 integer &
4041 ,intent (in ) :: &
4042 itf,jtf,ktf, &
4043 its,ite, jts,jte, kts,kte
4044 ! hc = cloud moist static energy
4045 ! hkb = moist static energy at originating level
4046 ! he = moist static energy on model levels
4047 ! he_cup = moist static energy on model cloud levels
4048 ! hes_cup = saturation moist static energy on model cloud levels
4049 ! dby = buoancy term
4050 ! cd= detrainment function
4051 ! z_cup = heights of model cloud levels
4052 ! entr = entrainment rate
4053 !
4054 real, dimension (its:ite,kts:kte) &
4055 ,intent (in ) :: &
4056 he,he_cup,hes_cup,z_cup,cd
4057 ! entr= entrainment rate
4058 real &
4059 ,intent (in ) :: &
4060 entr
4061 integer, dimension (its:ite) &
4062 ,intent (in ) :: &
4063 kbcon,k22
4064 !
4065 ! input and output
4066 !
4068 ! ierr error value, maybe modified in this routine
4070 integer, dimension (its:ite) &
4071 ,intent (inout) :: &
4072 ierr
4074 real, dimension (its:ite,kts:kte) &
4075 ,intent (out ) :: &
4076 hc,dby
4077 real, dimension (its:ite) &
4078 ,intent (out ) :: &
4079 hkb
4080 !
4081 ! local variables in this routine
4082 !
4084 integer :: &
4085 i,k
4086 real :: &
4087 dz
4088 !
4089 !--- moist static energy inside cloud
4090 !
4091 do k=kts,ktf
4092 do i=its,itf
4093 hc(i,k)=0.
4094 DBY(I,K)=0.
4095 enddo
4096 enddo
4097 do i=its,itf
4098 hkb(i)=0.
4099 enddo
4100 do i=its,itf
4101 if(ierr(i).eq.0.)then
4102 hkb(i)=he_cup(i,k22(i))
4103 do k=1,k22(i)
4104 hc(i,k)=he_cup(i,k)
4105 ! DBY(I,K)=0.
4106 enddo
4107 do k=k22(i),kbcon(i)-1
4108 hc(i,k)=hkb(i)
4109 ! DBY(I,K)=0.
4110 enddo
4111 k=kbcon(i)
4112 hc(i,k)=hkb(i)
4113 DBY(I,Kbcon(i))=Hkb(I)-HES_cup(I,K)
4114 endif
4115 enddo
4116 do k=kts+1,ktf
4117 do i=its,itf
4118 if(k.gt.kbcon(i).and.ierr(i).eq.0.)then
4119 DZ=Z_cup(i,K)-Z_cup(i,K-1)
4120 HC(i,K)=(HC(i,K-1)*(1.-.5*CD(i,K)*DZ)+entr* &
4121 DZ*HE(i,K-1))/(1.+entr*DZ-.5*cd(i,k)*dz)
4122 DBY(I,K)=HC(I,K)-HES_cup(I,K)
4123 endif
4124 enddo
4126 enddo
4128 END SUBROUTINE cup_up_he
4131 SUBROUTINE cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, &
4132 kbcon,ktop,cd,dby,mentr_rate,clw_all, &
4133 q,GAMMA_cup,zu,qes_cup,k22,qe_cup,xl, &
4134 itf,jtf,ktf, &
4135 its,ite, jts,jte, kts,kte )
4137 IMPLICIT NONE
4138 !
4139 ! on input
4140 !
4142 ! only local wrf dimensions are need as of now in this routine
4144 integer &
4145 ,intent (in ) :: &
4146 itf,jtf,ktf, &
4147 its,ite, jts,jte, kts,kte
4148 ! cd= detrainment function
4149 ! q = environmental q on model levels
4150 ! qe_cup = environmental q on model cloud levels
4151 ! qes_cup = saturation q on model cloud levels
4152 ! dby = buoancy term
4153 ! cd= detrainment function
4154 ! zu = normalized updraft mass flux
4155 ! gamma_cup = gamma on model cloud levels
4156 ! mentr_rate = entrainment rate
4157 !
4158 real, dimension (its:ite,kts:kte) &
4159 ,intent (in ) :: &
4160 q,zu,gamma_cup,qe_cup,dby,qes_cup,z_cup,cd
4161 ! entr= entrainment rate
4162 real &
4163 ,intent (in ) :: &
4164 mentr_rate,xl
4165 integer, dimension (its:ite) &
4166 ,intent (in ) :: &
4167 kbcon,ktop,k22
4168 !
4169 ! input and output
4170 !
4172 ! ierr error value, maybe modified in this routine
4174 integer, dimension (its:ite) &
4175 ,intent (inout) :: &
4176 ierr
4177 character *(*), intent (in) :: &
4178 name
4179 ! qc = cloud q (including liquid water) after entrainment
4180 ! qrch = saturation q in cloud
4181 ! qrc = liquid water content in cloud after rainout
4182 ! pw = condensate that will fall out at that level
4183 ! pwav = totan normalized integrated condensate (I1)
4184 ! c0 = conversion rate (cloud to rain)
4186 real, dimension (its:ite,kts:kte) &
4187 ,intent (out ) :: &
4188 qc,qrc,pw,clw_all
4189 real, dimension (its:ite) &
4190 ,intent (out ) :: &
4191 pwav
4192 !
4193 ! local variables in this routine
4194 !
4196 integer :: &
4197 iall,i,k
4198 real :: &
4199 dh,qrch,c0,dz,radius
4200 !
4201 iall=0
4202 c0=.002
4203 !
4204 !--- no precip for small clouds
4205 !
4206 if(name.eq.'shallow')c0=0.
4207 do i=its,itf
4208 pwav(i)=0.
4209 enddo
4210 do k=kts,ktf
4211 do i=its,itf
4212 pw(i,k)=0.
4213 qc(i,k)=0.
4214 if(ierr(i).eq.0)qc(i,k)=qes_cup(i,k)
4215 clw_all(i,k)=0.
4216 qrc(i,k)=0.
4217 enddo
4218 enddo
4219 do i=its,itf
4220 if(ierr(i).eq.0.)then
4221 do k=k22(i),kbcon(i)-1
4222 qc(i,k)=qe_cup(i,k22(i))
4223 enddo
4224 endif
4225 enddo
4227 DO 100 k=kts+1,ktf
4228 DO 100 i=its,itf
4229 IF(ierr(i).ne.0)GO TO 100
4230 IF(K.Lt.KBCON(I))GO TO 100
4231 IF(K.Gt.KTOP(I))GO TO 100
4232 DZ=Z_cup(i,K)-Z_cup(i,K-1)
4233 !
4234 !------ 1. steady state plume equation, for what could
4235 !------ be in cloud without condensation
4236 !
4237 !
4238 QC(i,K)=(QC(i,K-1)*(1.-.5*CD(i,K)*DZ)+mentr_rate* &
4239 DZ*Q(i,K-1))/(1.+mentr_rate*DZ-.5*cd(i,k)*dz)
4240 !
4241 !--- saturation in cloud, this is what is allowed to be in it
4242 !
4243 QRCH=QES_cup(I,K)+(1./XL)*(GAMMA_cup(i,k) &
4244 /(1.+GAMMA_cup(i,k)))*DBY(I,K)
4245 !
4246 !------- LIQUID WATER CONTENT IN cloud after rainout
4247 !
4248 clw_all(i,k)=QC(I,K)-QRCH
4249 QRC(I,K)=(QC(I,K)-QRCH)/(1.+C0*DZ)
4250 if(qrc(i,k).lt.0.)then
4251 qrc(i,k)=0.
4252 endif
4253 !
4254 !------- 3.Condensation
4255 !
4256 PW(i,k)=c0*dz*QRC(I,K)*zu(i,k)
4257 if(iall.eq.1)then
4258 qrc(i,k)=0.
4259 pw(i,k)=(QC(I,K)-QRCH)*zu(i,k)
4260 if(pw(i,k).lt.0.)pw(i,k)=0.
4261 endif
4262 !
4263 !----- set next level
4264 !
4265 QC(I,K)=QRC(I,K)+qrch
4266 !
4267 !--- integrated normalized ondensate
4268 !
4269 PWAV(I)=PWAV(I)+PW(I,K)
4270 100 CONTINUE
4272 END SUBROUTINE cup_up_moisture
4275 SUBROUTINE cup_up_nms(zu,z_cup,entr,cd,kbcon,ktop,ierr,k22, &
4276 itf,jtf,ktf, &
4277 its,ite, jts,jte, kts,kte )
4279 IMPLICIT NONE
4281 !
4282 ! on input
4283 !
4285 ! only local wrf dimensions are need as of now in this routine
4287 integer &
4288 ,intent (in ) :: &
4289 itf,jtf,ktf, &
4290 its,ite, jts,jte, kts,kte
4291 ! cd= detrainment function
4292 real, dimension (its:ite,kts:kte) &
4293 ,intent (in ) :: &
4294 z_cup,cd
4295 ! entr= entrainment rate
4296 real &
4297 ,intent (in ) :: &
4298 entr
4299 integer, dimension (its:ite) &
4300 ,intent (in ) :: &
4301 kbcon,ktop,k22
4302 !
4303 ! input and output
4304 !
4306 ! ierr error value, maybe modified in this routine
4308 integer, dimension (its:ite) &
4309 ,intent (inout) :: &
4310 ierr
4311 ! zu is the normalized mass flux
4313 real, dimension (its:ite,kts:kte) &
4314 ,intent (out ) :: &
4315 zu
4316 !
4317 ! local variables in this routine
4318 !
4320 integer :: &
4321 i,k
4322 real :: &
4323 dz
4324 !
4325 ! initialize for this go around
4326 !
4327 do k=kts,ktf
4328 do i=its,itf
4329 zu(i,k)=0.
4330 enddo
4331 enddo
4332 !
4333 ! do normalized mass budget
4334 !
4335 do i=its,itf
4336 IF(ierr(I).eq.0)then
4337 do k=k22(i),kbcon(i)
4338 zu(i,k)=1.
4339 enddo
4340 DO K=KBcon(i)+1,KTOP(i)
4341 DZ=Z_cup(i,K)-Z_cup(i,K-1)
4342 ZU(i,K)=ZU(i,K-1)*(1.+(entr-cd(i,k))*DZ)
4343 enddo
4344 endif
4345 enddo
4347 END SUBROUTINE cup_up_nms
4349 !====================================================================
4350 SUBROUTINE g3init(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
4351 MASS_FLUX,cp,restart, &
4352 P_QC,P_QI,P_FIRST_SCALAR, &
4353 RTHFTEN, RQVFTEN, &
4354 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
4355 APR_CAPMA,APR_CAPME,APR_CAPMI, &
4356 cugd_tten,cugd_ttens,cugd_qvten, &
4357 cugd_qvtens,cugd_qcten, &
4358 allowed_to_read, &
4359 ids, ide, jds, jde, kds, kde, &
4360 ims, ime, jms, jme, kms, kme, &
4361 its, ite, jts, jte, kts, kte )
4362 !--------------------------------------------------------------------
4363 IMPLICIT NONE
4364 !--------------------------------------------------------------------
4365 LOGICAL , INTENT(IN) :: restart,allowed_to_read
4366 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
4367 ims, ime, jms, jme, kms, kme, &
4368 its, ite, jts, jte, kts, kte
4369 INTEGER , INTENT(IN) :: P_FIRST_SCALAR, P_QI, P_QC
4370 REAL, INTENT(IN) :: cp
4372 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
4373 CUGD_TTEN, &
4374 CUGD_TTENS, &
4375 CUGD_QVTEN, &
4376 CUGD_QVTENS, &
4377 CUGD_QCTEN
4378 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
4379 RTHCUTEN, &
4380 RQVCUTEN, &
4381 RQCCUTEN, &
4382 RQICUTEN
4384 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
4385 RTHFTEN, &
4386 RQVFTEN
4388 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(OUT) :: &
4389 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
4390 APR_CAPMA,APR_CAPME,APR_CAPMI, &
4391 MASS_FLUX
4393 INTEGER :: i, j, k, itf, jtf, ktf
4395 jtf=min0(jte,jde-1)
4396 ktf=min0(kte,kde-1)
4397 itf=min0(ite,ide-1)
4399 IF(.not.restart)THEN
4400 DO j=jts,jte
4401 DO k=kts,kte
4402 DO i=its,ite
4403 RTHCUTEN(i,k,j)=0.
4404 RQVCUTEN(i,k,j)=0.
4405 ENDDO
4406 ENDDO
4407 ENDDO
4408 DO j=jts,jte
4409 DO k=kts,kte
4410 DO i=its,ite
4411 cugd_tten(i,k,j)=0.
4412 cugd_ttens(i,k,j)=0.
4413 cugd_qvten(i,k,j)=0.
4414 cugd_qvtens(i,k,j)=0.
4415 ENDDO
4416 ENDDO
4417 ENDDO
4419 DO j=jts,jtf
4420 DO k=kts,ktf
4421 DO i=its,itf
4422 RTHFTEN(i,k,j)=0.
4423 RQVFTEN(i,k,j)=0.
4424 ENDDO
4425 ENDDO
4426 ENDDO
4428 IF (P_QC .ge. P_FIRST_SCALAR) THEN
4429 DO j=jts,jtf
4430 DO k=kts,ktf
4431 DO i=its,itf
4432 RQCCUTEN(i,k,j)=0.
4433 cugd_qcten(i,k,j)=0.
4434 ENDDO
4435 ENDDO
4436 ENDDO
4437 ENDIF
4439 IF (P_QI .ge. P_FIRST_SCALAR) THEN
4440 DO j=jts,jtf
4441 DO k=kts,ktf
4442 DO i=its,itf
4443 RQICUTEN(i,k,j)=0.
4444 ENDDO
4445 ENDDO
4446 ENDDO
4447 ENDIF
4449 DO j=jts,jtf
4450 DO i=its,itf
4451 mass_flux(i,j)=0.
4452 ENDDO
4453 ENDDO
4455 ENDIF
4456 DO j=jts,jtf
4457 DO i=its,itf
4458 APR_GR(i,j)=0.
4459 APR_ST(i,j)=0.
4460 APR_W(i,j)=0.
4461 APR_MC(i,j)=0.
4462 APR_AS(i,j)=0.
4463 APR_CAPMA(i,j)=0.
4464 APR_CAPME(i,j)=0.
4465 APR_CAPMI(i,j)=0.
4466 ENDDO
4467 ENDDO
4469 END SUBROUTINE g3init
4472 SUBROUTINE massflx_stats(xf_ens,ensdim,maxens,maxens2,maxens3, &
4473 xt_ave,xt_std,xt_cur,xt_ske,j,ierr,itest, &
4474 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
4475 APR_CAPMA,APR_CAPME,APR_CAPMI, &
4476 pr_gr,pr_w,pr_mc,pr_st,pr_as, &
4477 pr_capma,pr_capme,pr_capmi, &
4478 itf,jtf,ktf, &
4479 its,ite, jts,jte, kts,kte)
4481 IMPLICIT NONE
4483 integer, intent (in ) :: &
4484 j,ensdim,maxens3,maxens,maxens2,itest
4485 INTEGER, INTENT(IN ) :: &
4486 itf,jtf,ktf, &
4487 its,ite, jts,jte, kts,kte
4490 real, dimension (its:ite) &
4491 , intent(inout) :: &
4492 xt_ave,xt_cur,xt_std,xt_ske
4493 integer, dimension (its:ite), intent (in) :: &
4494 ierr
4495 real, dimension (its:ite,jts:jte,1:ensdim) &
4496 , intent(in ) :: &
4497 xf_ens
4498 real, dimension (its:ite,jts:jte) &
4499 , intent(inout) :: &
4500 APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
4501 APR_CAPMA,APR_CAPME,APR_CAPMI
4502 real, dimension (its:ite,jts:jte) &
4503 , intent(inout) :: &
4504 pr_gr,pr_w,pr_mc,pr_st,pr_as, &
4505 pr_capma,pr_capme,pr_capmi
4507 !
4508 ! local stuff
4509 !
4510 real, dimension (its:ite , 1:maxens3 ) :: &
4511 x_ave,x_cur,x_std,x_ske
4512 real, dimension (its:ite , 1:maxens ) :: &
4513 x_ave_cap
4516 integer, dimension (1:maxens3) :: nc1
4517 integer :: i,k
4518 integer :: num,kk,num2,iedt
4519 real :: a3,a4
4521 num=ensdim/maxens3
4522 num2=ensdim/maxens
4523 if(itest.eq.1)then
4524 do i=its,ite
4525 pr_gr(i,j) = 0.
4526 pr_w(i,j) = 0.
4527 pr_mc(i,j) = 0.
4528 pr_st(i,j) = 0.
4529 pr_as(i,j) = 0.
4530 pr_capma(i,j) = 0.
4531 pr_capme(i,j) = 0.
4532 pr_capmi(i,j) = 0.
4533 enddo
4534 endif
4536 do k=1,maxens
4537 do i=its,ite
4538 x_ave_cap(i,k)=0.
4539 enddo
4540 enddo
4541 do k=1,maxens3
4542 do i=its,ite
4543 x_ave(i,k)=0.
4544 x_std(i,k)=0.
4545 x_ske(i,k)=0.
4546 x_cur(i,k)=0.
4547 enddo
4548 enddo
4549 do i=its,ite
4550 xt_ave(i)=0.
4551 xt_std(i)=0.
4552 xt_ske(i)=0.
4553 xt_cur(i)=0.
4554 enddo
4555 do kk=1,num
4556 do k=1,maxens3
4557 do i=its,ite
4558 if(ierr(i).eq.0)then
4559 x_ave(i,k)=x_ave(i,k)+xf_ens(i,j,maxens3*(kk-1)+k)
4560 endif
4561 enddo
4562 enddo
4563 enddo
4564 do iedt=1,maxens2
4565 do k=1,maxens
4566 do kk=1,maxens3
4567 do i=its,ite
4568 if(ierr(i).eq.0)then
4569 x_ave_cap(i,k)=x_ave_cap(i,k) &
4570 +xf_ens(i,j,maxens3*(k-1)+(iedt-1)*maxens*maxens3+kk)
4571 endif
4572 enddo
4573 enddo
4574 enddo
4575 enddo
4576 do k=1,maxens
4577 do i=its,ite
4578 if(ierr(i).eq.0)then
4579 x_ave_cap(i,k)=x_ave_cap(i,k)/float(num2)
4580 endif
4581 enddo
4582 enddo
4584 do k=1,maxens3
4585 do i=its,ite
4586 if(ierr(i).eq.0)then
4587 x_ave(i,k)=x_ave(i,k)/float(num)
4588 endif
4589 enddo
4590 enddo
4591 do k=1,maxens3
4592 do i=its,ite
4593 if(ierr(i).eq.0)then
4594 xt_ave(i)=xt_ave(i)+x_ave(i,k)
4595 endif
4596 enddo
4597 enddo
4598 do i=its,ite
4599 if(ierr(i).eq.0)then
4600 xt_ave(i)=xt_ave(i)/float(maxens3)
4601 endif
4602 enddo
4603 !
4604 !--- now do std, skewness,curtosis
4605 !
4606 do kk=1,num
4607 do k=1,maxens3
4608 do i=its,ite
4609 if(ierr(i).eq.0.and.x_ave(i,k).gt.0.)then
4610 ! print *,i,j,k,kk,x_std(i,k),xf_ens(i,j,maxens3*(kk-1)+k),x_ave(i,k)
4611 x_std(i,k)=x_std(i,k)+(xf_ens(i,j,maxens3*(kk-1)+k)-x_ave(i,k))**2
4612 x_ske(i,k)=x_ske(i,k)+(xf_ens(i,j,maxens3*(kk-1)+k)-x_ave(i,k))**3
4613 x_cur(i,k)=x_cur(i,k)+(xf_ens(i,j,maxens3*(kk-1)+k)-x_ave(i,k))**4
4614 endif
4615 enddo
4616 enddo
4617 enddo
4618 do k=1,maxens3
4619 do i=its,ite
4620 if(ierr(i).eq.0.and.xt_ave(i).gt.0.)then
4621 xt_std(i)=xt_std(i)+(x_ave(i,k)-xt_ave(i))**2
4622 xt_ske(i)=xt_ske(i)+(x_ave(i,k)-xt_ave(i))**3
4623 xt_cur(i)=xt_cur(i)+(x_ave(i,k)-xt_ave(i))**4
4624 endif
4625 enddo
4626 enddo
4627 do k=1,maxens3
4628 do i=its,ite
4629 if(ierr(i).eq.0.and.x_std(i,k).gt.0.)then
4630 x_std(i,k)=x_std(i,k)/float(num)
4631 a3=max(1.e-6,x_std(i,k))
4632 x_std(i,k)=sqrt(a3)
4633 a3=max(1.e-6,x_std(i,k)**3)
4634 a4=max(1.e-6,x_std(i,k)**4)
4635 x_ske(i,k)=x_ske(i,k)/float(num)/a3
4636 x_cur(i,k)=x_cur(i,k)/float(num)/a4
4637 endif
4638 ! print*,' '
4639 ! print*,'Some statistics at gridpoint i,j, ierr',i,j,ierr(i)
4640 ! print*,'statistics for closure number ',k
4641 ! print*,'Average= ',x_ave(i,k),' Std= ',x_std(i,k)
4642 ! print*,'Skewness= ',x_ske(i,k),' Curtosis= ',x_cur(i,k)
4643 ! print*,' '
4645 enddo
4646 enddo
4647 do i=its,ite
4648 if(ierr(i).eq.0.and.xt_std(i).gt.0.)then
4649 xt_std(i)=xt_std(i)/float(maxens3)
4650 a3=max(1.e-6,xt_std(i))
4651 xt_std(i)=sqrt(a3)
4652 a3=max(1.e-6,xt_std(i)**3)
4653 a4=max(1.e-6,xt_std(i)**4)
4654 xt_ske(i)=xt_ske(i)/float(maxens3)/a3
4655 xt_cur(i)=xt_cur(i)/float(maxens3)/a4
4656 ! print*,' '
4657 ! print*,'Total ensemble independent statistics at i =',i
4658 ! print*,'Average= ',xt_ave(i),' Std= ',xt_std(i)
4659 ! print*,'Skewness= ',xt_ske(i),' Curtosis= ',xt_cur(i)
4660 ! print*,' '
4661 !
4662 ! first go around: store massflx for different closures/caps
4663 !
4664 if(itest.eq.1)then
4665 pr_gr(i,j) = .25*(x_ave(i,1)+x_ave(i,2)+x_ave(i,3)+x_ave(i,13))
4666 pr_w(i,j) = .25*(x_ave(i,4)+x_ave(i,5)+x_ave(i,6)+x_ave(i,14))
4667 pr_mc(i,j) = .25*(x_ave(i,7)+x_ave(i,8)+x_ave(i,9)+x_ave(i,15))
4668 pr_st(i,j) = .333*(x_ave(i,10)+x_ave(i,11)+x_ave(i,12))
4669 pr_as(i,j) = x_ave(i,16)
4670 pr_capma(i,j) = x_ave_cap(i,1)
4671 pr_capme(i,j) = x_ave_cap(i,2)
4672 pr_capmi(i,j) = x_ave_cap(i,3)
4673 !
4674 ! second go around: store preciprates (mm/hour) for different closures/caps
4675 !
4676 else if (itest.eq.2)then
4677 APR_GR(i,j)=.25*(x_ave(i,1)+x_ave(i,2)+x_ave(i,3)+x_ave(i,13))* &
4678 3600.*pr_gr(i,j) +APR_GR(i,j)
4679 APR_W(i,j)=.25*(x_ave(i,4)+x_ave(i,5)+x_ave(i,6)+x_ave(i,14))* &
4680 3600.*pr_w(i,j) +APR_W(i,j)
4681 APR_MC(i,j)=.25*(x_ave(i,7)+x_ave(i,8)+x_ave(i,9)+x_ave(i,15))* &
4682 3600.*pr_mc(i,j) +APR_MC(i,j)
4683 APR_ST(i,j)=.333*(x_ave(i,10)+x_ave(i,11)+x_ave(i,12))* &
4684 3600.*pr_st(i,j) +APR_ST(i,j)
4685 APR_AS(i,j)=x_ave(i,16)* &
4686 3600.*pr_as(i,j) +APR_AS(i,j)
4687 APR_CAPMA(i,j) = x_ave_cap(i,1)* &
4688 3600.*pr_capma(i,j) +APR_CAPMA(i,j)
4689 APR_CAPME(i,j) = x_ave_cap(i,2)* &
4690 3600.*pr_capme(i,j) +APR_CAPME(i,j)
4691 APR_CAPMI(i,j) = x_ave_cap(i,3)* &
4692 3600.*pr_capmi(i,j) +APR_CAPMI(i,j)
4693 endif
4694 endif
4695 enddo
4697 END SUBROUTINE massflx_stats
4699 SUBROUTINE cup_axx(tcrit,kbmax,z1,p,psur,xl,rv,cp,tx,qx,axx,ierr, &
4700 cap_max,cap_max_increment,entr_rate,mentr_rate,&
4701 j,itf,jtf,ktf, &
4702 its,ite, jts,jte, kts,kte,ens4)
4703 IMPLICIT NONE
4704 INTEGER, INTENT(IN ) :: &
4705 j,itf,jtf,ktf, &
4706 its,ite, jts,jte, kts,kte,ens4
4707 real, dimension (its:ite,kts:kte,1:ens4) &
4708 , intent(inout) :: &
4709 tx,qx
4710 real, dimension (its:ite,kts:kte) &
4711 , intent(in) :: &
4712 p
4713 real, dimension (its:ite) &
4714 , intent(in) :: &
4715 z1,psur,cap_max,cap_max_increment
4716 real, intent(in) :: &
4717 tcrit,xl,rv,cp,mentr_rate,entr_rate
4718 real, dimension (its:ite,1:ens4) &
4719 , intent(out) :: &
4720 axx
4721 integer, dimension (its:ite), intent (in) :: &
4722 ierr,kbmax
4723 integer, dimension (its:ite) :: &
4724 ierrxx,k22xx,kbconxx,ktopxx,kstabm,kstabi
4725 real, dimension (1:2) :: AE,BE,HT
4726 real, dimension (its:ite,kts:kte) :: tv
4727 real :: e,tvbar
4728 integer n,i,k,iph
4729 real, dimension (its:ite,kts:kte) :: &
4730 he,hes,qes,z, &
4731 qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
4732 tn_cup, &
4733 dby,qc,qrcd,pwd,pw,hcd,qcd,dbyd,hc,qrc,zu,zd,cd
4735 real, dimension (its:ite) :: &
4736 AA0,HKB,QKB, &
4737 PWAV,BU
4738 do n=1,ens4
4739 do i=its,ite
4740 axx(i,n)=0.
4741 enddo
4742 enddo
4743 HT(1)=XL/CP
4744 HT(2)=2.834E6/CP
4745 BE(1)=.622*HT(1)/.286
4746 AE(1)=BE(1)/273.+ALOG(610.71)
4747 BE(2)=.622*HT(2)/.286
4748 AE(2)=BE(2)/273.+ALOG(610.71)
4749 !
4750 !
4751 do 100 n=1,ens4
4753 do k=kts,ktf
4754 do i=its,itf
4755 cd(i,k)=0.1*entr_rate
4756 enddo
4757 enddo
4760 do i=its,itf
4761 ierrxx(i)=ierr(i)
4762 k22xx(i)=1
4763 kbconxx(i)=1
4764 ktopxx(i)=1
4765 kstabm(i)=ktf-1
4766 enddo
4767 DO k=kts,ktf
4768 do i=its,itf
4769 if(ierrxx(i).eq.0)then
4770 IPH=1
4771 IF(Tx(I,K,n).LE.TCRIT)IPH=2
4772 E=EXP(AE(IPH)-BE(IPH)/TX(I,K,N))
4773 QES(I,K)=.622*E/(100.*P(I,K)-E)
4774 IF(QES(I,K).LE.1.E-08)QES(I,K)=1.E-08
4775 IF(Qx(I,K,N).GT.QES(I,K))Qx(I,K,N)=QES(I,K)
4776 TV(I,K)=Tx(I,K,N)+.608*Qx(I,K,N)*Tx(I,K,N)
4777 endif
4778 enddo
4779 enddo
4780 !
4781 do i=its,itf
4782 if(ierrxx(i).eq.0)then
4783 Z(I,KTS)=max(0.,Z1(I))-(ALOG(P(I,KTS))- &
4784 ALOG(PSUR(I)))*287.*TV(I,KTS)/9.81
4785 endif
4786 enddo
4788 ! --- calculate heights
4789 DO K=kts+1,ktf
4790 do i=its,itf
4791 if(ierrxx(i).eq.0)then
4792 TVBAR=.5*TV(I,K)+.5*TV(I,K-1)
4793 Z(I,K)=Z(I,K-1)-(ALOG(P(I,K))- &
4794 ALOG(P(I,K-1)))*287.*TVBAR/9.81
4795 endif
4796 enddo
4797 enddo
4798 !
4799 !--- calculate moist static energy - HE
4800 ! saturated moist static energy - HES
4801 !
4802 DO k=kts,ktf
4803 do i=its,itf
4804 if(ierrxx(i).eq.0)then
4805 HE(I,K)=9.81*Z(I,K)+1004.*Tx(I,K,n)+2.5E06*Qx(I,K,n)
4806 HES(I,K)=9.81*Z(I,K)+1004.*Tx(I,K,n)+2.5E06*QES(I,K)
4807 IF(HE(I,K).GE.HES(I,K))HE(I,K)=HES(I,K)
4808 endif
4809 enddo
4810 enddo
4812 ! cup levels
4813 !
4814 do k=kts+1,ktf
4815 do i=its,itf
4816 if(ierrxx(i).eq.0)then
4817 qes_cup(i,k)=.5*(qes(i,k-1)+qes(i,k))
4818 q_cup(i,k)=.5*(qx(i,k-1,n)+qx(i,k,n))
4819 hes_cup(i,k)=.5*(hes(i,k-1)+hes(i,k))
4820 he_cup(i,k)=.5*(he(i,k-1)+he(i,k))
4821 if(he_cup(i,k).gt.hes_cup(i,k))he_cup(i,k)=hes_cup(i,k)
4822 z_cup(i,k)=.5*(z(i,k-1)+z(i,k))
4823 p_cup(i,k)=.5*(p(i,k-1)+p(i,k))
4824 t_cup(i,k)=.5*(tx(i,k-1,n)+tx(i,k,n))
4825 gamma_cup(i,k)=(xl/cp)*(xl/(rv*t_cup(i,k) &
4826 *t_cup(i,k)))*qes_cup(i,k)
4827 endif
4828 enddo
4829 enddo
4830 do i=its,itf
4831 if(ierrxx(i).eq.0)then
4832 qes_cup(i,1)=qes(i,1)
4833 q_cup(i,1)=qx(i,1,n)
4834 hes_cup(i,1)=hes(i,1)
4835 he_cup(i,1)=he(i,1)
4836 z_cup(i,1)=.5*(z(i,1)+z1(i))
4837 p_cup(i,1)=.5*(p(i,1)+psur(i))
4838 t_cup(i,1)=tx(i,1,n)
4839 gamma_cup(i,1)=xl/cp*(xl/(rv*t_cup(i,1) &
4840 *t_cup(i,1)))*qes_cup(i,1)
4841 endif
4842 enddo
4843 !
4844 !
4845 !------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
4846 !
4847 CALL cup_MAXIMI(HE_CUP,3,KBMAX,K22XX,ierrxx, &
4848 itf,jtf,ktf, &
4849 its,ite, jts,jte, kts,kte)
4850 DO 36 i=its,itf
4851 IF(ierrxx(I).eq.0.)THEN
4852 IF(K22xx(I).GE.KBMAX(i))ierrxx(i)=2
4853 endif
4854 36 CONTINUE
4855 !
4856 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
4857 !
4858 call cup_kbcon(cap_max_increment,1,k22xx,kbconxx,he_cup,hes_cup, &
4859 ierrxx,kbmax,p_cup,cap_max, &
4860 itf,jtf,ktf, &
4861 its,ite, jts,jte, kts,kte)
4862 !
4863 !--- increase detrainment in stable layers
4864 !
4865 CALL cup_minimi(HEs_cup,Kbconxx,kstabm,kstabi,ierrxx, &
4866 itf,jtf,ktf, &
4867 its,ite, jts,jte, kts,kte)
4868 do i=its,itf
4869 IF(ierrxx(I).eq.0.)THEN
4870 if(kstabm(i)-1.gt.kstabi(i))then
4871 do k=kstabi(i),kstabm(i)-1
4872 cd(i,k)=cd(i,k-1)+1.5*entr_rate
4873 if(cd(i,k).gt.10.0*entr_rate)cd(i,k)=10.0*entr_rate
4874 enddo
4875 ENDIF
4876 ENDIF
4877 ENDDO
4878 !
4879 !--- calculate incloud moist static energy
4880 !
4881 call cup_up_he(k22xx,hkb,z_cup,cd,mentr_rate,he_cup,hc, &
4882 kbconxx,ierrxx,dby,he,hes_cup, &
4883 itf,jtf,ktf, &
4884 its,ite, jts,jte, kts,kte)
4886 !--- DETERMINE CLOUD TOP - KTOP
4887 !
4888 call cup_ktop(1,dby,kbconxx,ktopxx,ierrxx, &
4889 itf,jtf,ktf, &
4890 its,ite, jts,jte, kts,kte)
4891 !
4892 !c--- normalized updraft mass flux profile
4893 !
4894 call cup_up_nms(zu,z_cup,mentr_rate,cd,kbconxx,ktopxx,ierrxx,k22xx, &
4895 itf,jtf,ktf, &
4896 its,ite, jts,jte, kts,kte)
4897 !
4898 !--- calculate workfunctions for updrafts
4899 !
4900 call cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
4901 kbconxx,ktopxx,ierrxx, &
4902 itf,jtf,ktf, &
4903 its,ite, jts,jte, kts,kte)
4904 do i=its,itf
4905 if(ierrxx(i).eq.0)axx(i,n)=aa0(i)
4906 enddo
4907 100 continue
4908 END SUBROUTINE cup_axx
4910 SUBROUTINE conv_grell_spread3d(rthcuten,rqvcuten,rqccuten,raincv, &
4911 & cugd_avedx,cugd_tten,cugd_qvten,rqicuten,cugd_ttens, &
4912 & cugd_qvtens,cugd_qcten,pi_phy,moist_qv,pratec,dt,num_tiles,&
4913 & imomentum,F_QV ,F_QC ,F_QR ,F_QI ,F_QS, &
4914 & ids, ide, jds, jde, kds, kde, &
4915 & ips, ipe, jps, jpe, kps, kpe, &
4916 & ims, ime, jms, jme, kms, kme, &
4917 & its, ite, jts, jte, kts, kte )
4919 !
4921 INTEGER, INTENT(IN ) :: num_tiles,imomentum
4922 INTEGER, INTENT(IN ) :: ids, ide, jds, jde, kds, kde,&
4923 ims,ime, jms,jme, kms,kme, &
4924 ips,ipe, jps,jpe, kps,kpe, &
4925 its,ite, jts,jte, kts,kte, &
4926 cugd_avedx
4927 REAL, DIMENSION (ims:ime,kms:kme,jms:jme), optional,INTENT (INOUT) :: &
4928 & rthcuten,rqvcuten,rqccuten,rqicuten
4929 REAL, DIMENSION (ims:ime,kms:kme,jms:jme), optional,INTENT (IN ) :: &
4930 & cugd_tten,cugd_qvten,cugd_ttens,cugd_qvtens,cugd_qcten
4931 REAL, DIMENSION (ims:ime,kms:kme,jms:jme),INTENT (IN) :: &
4932 moist_qv
4933 REAL, DIMENSION (ims:ime,kms:kme,jms:jme), INTENT (IN) :: &
4934 PI_PHY
4935 REAL, DIMENSION (ims:ime,jms:jme), INTENT (INOUT) :: &
4936 raincv,pratec
4937 REAL, INTENT(IN) :: dt
4938 INTEGER :: ikk1,ikk2,ikk11,i,j,k,kk,nn,smoothh,smoothv
4939 INTEGER :: ifs,ife,jfs,jfe,ido,jdo,cugd_spread
4940 LOGICAL :: new
4941 !
4942 ! Flags relating to the optional tendency arrays declared above
4943 ! Models that carry the optional tendencies will provdide the
4944 ! optional arguments at compile time; these flags all the model
4945 ! to determine at run-time whether a particular tracer is in
4946 ! use or not.
4947 !
4948 LOGICAL, OPTIONAL :: &
4949 F_QV &
4950 ,F_QC &
4951 ,F_QR &
4952 ,F_QI &
4953 ,F_QS
4954 REAL, DIMENSION (its-2:ite+2,kts:kte,jts-2:jte+2) :: &
4955 RTHcutent,RQVcutent
4956 real, dimension (its-2:ite+2,jts-2:jte+2) :: Qmem
4957 real, dimension (its-1:ite+1,jts-1:jte+1) :: smTT,smTQ
4958 real, dimension (kts:kte) :: conv_TRASHT,conv_TRASHQ
4959 REAL :: Qmem1,Qmem2,Qmemf,Thresh
4961 smoothh=1
4962 smoothv=1
4963 cugd_spread=cugd_avedx/2
4965 ifs=max(its,ids)
4966 jfs=max(jts,jds)
4967 ife=min(ite,ide-1)
4968 jfe=min(jte,jde-1)
4970 do j=jfs-2,jfe+2
4971 do i=ifs-2,ife+2
4972 Qmem(i,j)=1.
4973 enddo
4974 enddo
4975 do j=jfs-1,jfe+1
4976 do i=ifs-1,ife+1
4977 smTT(i,j)=0.
4978 smTQ(i,j)=0.
4979 enddo
4980 enddo
4981 do j=jfs,jfe
4982 do k=kts,kte
4983 do i=ifs,ife
4984 rthcuten(i,k,j)=0.
4985 rqvcuten(i,k,j)=0.
4986 enddo
4987 enddo
4988 enddo
4989 do j=jfs-2,jfe+2
4990 do k=kts,kte
4991 do i=ifs-2,ife+2
4992 RTHcutent(i,k,j)=0.
4993 RQVcutent(i,k,j)=0.
4994 enddo
4995 enddo
4996 enddo
4997 !
4998 !
4999 !
5000 !
5001 ! prelims finished, now go real for every grid point
5002 !
5003 if(cugd_spread.gt.0.or.smoothh.eq.1)then
5004 !if(its.eq.ips)ifs=max(its-1,ids)
5005 !if(ite.eq.ipe)ife=min(ite+1,ide-1)
5006 !if(jts.eq.jps)jfs=max(jts-1,jds)
5007 !if(jte.eq.jpe)jfe=min(jte+1,jde-1)
5008 ifs=max(its-1,ids)
5009 ife=min(ite+1,ide-1)
5010 jfs=max(jts-1,jds)
5011 jfe=min(jte+1,jde-1)
5012 endif
5014 ! *** jm note -- for smoothing this goes out one row/column beyond tile in i and j
5015 do j=jfs,jfe
5016 do i=ifs,ife
5017 !
5018 do k=kts,kte
5019 RTHcutent(i,k,j)=cugd_tten(i,k,j)
5020 RQVcutent(i,k,j)=cugd_qvten(i,k,j)
5021 enddo
5022 !
5023 ! for high res run, spread the subsidence
5024 ! this is tricky......only consider grid points where there was no rain,
5025 ! so cugd_tten and such are zero!
5026 !
5027 if(cugd_spread.gt.0)then
5028 do k=kts,kte
5029 do nn=-1,1,1
5030 jdo=max(j+nn,jds)
5031 jdo=min(jdo,jde-1)
5032 do kk=-1,1,1
5033 ido=max(i+kk,ids)
5034 ido=min(ido,ide-1)
5035 RTHcutent(i,k,j)=RTHcutent(i,k,j) &
5036 +Qmem(ido,jdo)*cugd_ttens(ido,k,jdo)
5037 RQVcutent(i,k,j)=RQVcutent(i,k,j) &
5038 +Qmem(ido,jdo)*cugd_qvtens(ido,k,jdo)
5039 enddo
5040 enddo
5041 enddo
5042 endif
5043 !
5044 ! end spreading
5046 if(cugd_spread.eq.0)then
5047 do k=kts,kte
5048 RTHcutent(i,k,j)=RTHcutent(i,k,j)+cugd_ttens(i,k,j)
5049 RQVcutent(i,k,j)=RQVcutent(i,k,j)+cugd_qvtens(i,k,j)
5050 enddo
5051 endif
5052 enddo ! end j
5053 enddo ! end i
5055 ! smooth
5056 do k=kts,kte
5057 if(smoothh.eq.0)then
5058 ifs=max(its,ids+4)
5059 ife=min(ite,ide-5)
5060 jfs=max(jts,jds+4)
5061 jfe=min(jte,jde-5)
5062 do i=ifs,ife
5063 do j=jfs,jfe
5064 rthcuten(i,k,j)=RTHcutent(i,k,j)
5065 rqvcuten(i,k,j)=RQVcutent(i,k,j)
5066 enddo ! end j
5067 enddo ! end j
5068 else if(smoothh.eq.1)then ! smooth
5069 ifs=max(its,ids)
5070 ife=min(ite,ide-1)
5071 jfs=max(jts,jds)
5072 jfe=min(jte,jde-1)
5073 ! we need an extra row for j (halo comp)
5074 !if(jts.eq.jps)jfs=max(jts-1,jds)
5075 !if(jte.eq.jpe)jfe=min(jte+1,jde-1)
5076 jfs=max(jts-1,jds)
5077 jfe=min(jte+1,jde-1)
5078 do i=ifs,ife
5079 do j=jfs,jfe
5080 smTT(i,j)=.25*(RTHcutent(i-1,k,j)+2.*RTHcutent(i,k,j)+RTHcutent(i+1,k,j))
5081 smTQ(i,j)=.25*(RQVcutent(i-1,k,j)+2.*RQVcutent(i,k,j)+RQVcutent(i+1,k,j))
5082 enddo ! end j
5083 enddo ! end j
5084 ifs=max(its,ids+4)
5085 ife=min(ite,ide-5)
5086 jfs=max(jts,jds+4)
5087 jfe=min(jte,jde-5)
5088 do i=ifs,ife
5089 do j=jfs,jfe
5090 rthcuten(i,k,j)=.25*(smTT(i,j-1)+2.*smTT(i,j)+smTT(i,j+1))
5091 rqvcuten(i,k,j)=.25*(smTQ(i,j-1)+2.*smTQ(i,j)+smTQ(i,j+1))
5092 enddo ! end j
5093 enddo ! end i
5094 endif ! smoothh
5095 enddo ! end k
5096 !
5097 ! check moistening rates
5098 !
5099 ifs=max(its,ids+4)
5100 ife=min(ite,ide-5)
5101 jfs=max(jts,jds+4)
5102 jfe=min(jte,jde-5)
5103 do j=jfs,jfe
5104 do i=ifs,ife
5105 Qmemf=1.
5106 Thresh=1.e-20
5107 do k=kts,kte
5108 if(rqvcuten(i,k,j).lt.0.)then
5109 Qmem1=moist_qv(i,k,j)+rqvcuten(i,k,j)*dt
5110 if(Qmem1.lt.Thresh)then
5111 Qmem1=rqvcuten(i,k,j)
5112 Qmem2=(Thresh-moist_qv(i,k,j))/dt
5113 Qmemf=min(Qmemf,Qmem2/Qmem1)
5114 Qmemf=max(0.,Qmemf)
5115 Qmemf=min(1.,Qmemf)
5116 endif
5117 endif
5118 enddo
5119 do k=kts,kte
5120 rqvcuten(i,k,j)=rqvcuten(i,k,j)*Qmemf
5121 rthcuten(i,k,j)=rthcuten(i,k,j)*Qmemf
5122 enddo
5123 if(present(rqccuten))then
5124 if(f_qc) then
5125 do k=kts,kte
5126 rqccuten(i,k,j)=rqccuten(i,k,j)*Qmemf
5127 enddo
5128 endif
5129 endif
5130 if(present(rqicuten))then
5131 if(f_qi) then
5132 do k=kts,kte
5133 rqicuten(i,k,j)=rqicuten(i,k,j)*Qmemf
5134 enddo
5135 endif
5136 endif
5137 raincv(I,J)=raincv(I,J)*Qmemf
5138 pratec(I,J)=pratec(I,J)*Qmemf
5139 !
5140 ! check heating rates
5142 !
5143 Thresh=200.
5144 Qmemf=1.
5145 Qmem1=0.
5146 do k=kts,kte
5147 Qmem1=abs(rthcuten(i,k,j))*86400.
5149 if(Qmem1.gt.Thresh)then
5150 Qmem2=Thresh/Qmem1
5151 Qmemf=min(Qmemf,Qmem2)
5152 Qmemf=max(0.,Qmemf)
5153 endif
5155 enddo
5156 raincv(i,j)=raincv(i,j)*Qmemf
5157 pratec(i,j)=pratec(i,j)*Qmemf
5158 do k=kts,kte
5159 rqvcuten(i,k,j)=rqvcuten(i,k,j)*Qmemf
5160 rthcuten(i,k,j)=rthcuten(i,k,j)*Qmemf
5161 enddo
5162 if(present(rqccuten))then
5163 if(f_qc) then
5164 do k=kts,kte
5165 rqccuten(i,k,j)=rqccuten(i,k,j)*Qmemf
5166 enddo
5167 endif
5168 endif
5169 if(present(rqicuten))then
5170 if(f_qi) then
5171 do k=kts,kte
5172 rqicuten(i,k,j)=rqicuten(i,k,j)*Qmemf
5173 enddo
5174 endif
5175 endif
5176 if(smoothv.eq.1)then
5177 !
5178 ! smooth for now
5179 !
5180 do k=kts+2,kte-2
5181 conv_TRASHT(k)= .25*(rthcuten(i,k-1,j)+2.*rthcuten(i,k,j)+rthcuten(i,k+1,j))
5182 conv_TRASHQ(k)= .25*(rqvcuten(i,k-1,j)+2.*rqvcuten(i,k,j)+rqvcuten(i,k+1,j))
5183 enddo
5184 do k=kts+2,kte-2
5185 rthcuten(i,k,j)=conv_TRASHT(k)
5186 rqvcuten(i,k,j)=conv_TRASHQ(k)
5187 enddo
5188 endif
5189 do k=kts,kte
5190 rthcuten(i,k,j)=rthcuten(i,k,j)/pi_phy(i,k,j)
5191 enddo
5192 enddo ! end j
5193 enddo ! end i
5195 END SUBROUTINE CONV_GRELL_SPREAD3D
5196 !-------------------------------------------------------
5197 END MODULE module_cu_g3