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