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