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