| blob | 2481a8f6fb5cf8b29be68f3cf70c7eaf6f94943b |
1 !WRF:MODEL_LAYER:PHYSICS
2 !
4 MODULE module_cu_kf
6 USE module_wrf_error
8 REAL , PARAMETER :: RAD = 1500.
10 CONTAINS
12 !-------------------------------------------------------------
13 SUBROUTINE KFCPS( &
14 ids,ide, jds,jde, kds,kde &
15 ,ims,ime, jms,jme, kms,kme &
16 ,its,ite, jts,jte, kts,kte &
17 ,DT,KTAU,DX,CUDT,ADAPT_STEP_FLAG &
18 ,rho &
19 ,RAINCV,PRATEC,NCA &
20 ,U,V,TH,T,W,QV,dz8w,Pcps,pi &
21 ,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
22 ,EP2,SVP1,SVP2,SVP3,SVPT0 &
23 ,STEPCU,CU_ACT_FLAG,warm_rain &
24 ! optional arguments
25 ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
26 ,RQVCUTEN,RQCCUTEN,RQRCUTEN,RQICUTEN,RQSCUTEN &
27 ,RTHCUTEN &
28 !kf_edrates
29 ,UDR_KF,DDR_KF &
30 ,UER_KF,DER_KF &
31 ,TIMEC_KF,KF_EDRATES &
32 )
33 !-------------------------------------------------------------
34 IMPLICIT NONE
35 !-------------------------------------------------------------
36 INTEGER, INTENT(IN ) :: &
37 ids,ide, jds,jde, kds,kde, &
38 ims,ime, jms,jme, kms,kme, &
39 its,ite, jts,jte, kts,kte
41 INTEGER, INTENT(IN ) :: STEPCU
42 LOGICAL, INTENT(IN ) :: warm_rain
44 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1
45 REAL, INTENT(IN ) :: CP,R,G,EP1,EP2
46 REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
48 INTEGER, INTENT(IN ) :: KTAU
50 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
51 INTENT(IN ) :: &
52 U, &
53 V, &
54 W, &
55 TH, &
56 QV, &
57 T, &
58 dz8w, &
59 Pcps, &
60 rho, &
61 pi
62 !
63 REAL, INTENT(IN ) :: DT, DX
64 REAL, INTENT(IN ) :: CUDT
65 LOGICAL,OPTIONAL, INTENT(IN ) :: ADAPT_STEP_FLAG
67 REAL, DIMENSION( ims:ime , jms:jme ), &
68 INTENT(INOUT) :: &
69 RAINCV &
70 ,PRATEC &
71 , NCA
73 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
74 INTENT(INOUT) :: &
75 W0AVG
77 LOGICAL, DIMENSION( ims:ime , jms:jme ), &
78 INTENT(INOUT) :: CU_ACT_FLAG
80 !
81 ! Optional arguments
82 !
84 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
85 OPTIONAL, &
86 INTENT(INOUT) :: &
87 RTHCUTEN &
88 ,RQVCUTEN &
89 ,RQCCUTEN &
90 ,RQRCUTEN &
91 ,RQICUTEN &
92 ,RQSCUTEN
94 !
95 ! Flags relating to the optional tendency arrays declared above
96 ! Models that carry the optional tendencies will provdide the
97 ! optional arguments at compile time; these flags all the model
98 ! to determine at run-time whether a particular tracer is in
99 ! use or not.
100 !
102 LOGICAL, OPTIONAL :: &
103 F_QV &
104 ,F_QC &
105 ,F_QR &
106 ,F_QI &
107 ,F_QS
109 !kf_edrates
110 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
111 INTENT(INOUT) :: &
112 UDR_KF, &
113 DDR_KF, &
114 UER_KF, &
115 DER_KF
117 REAL, DIMENSION( ims:ime , jms:jme ) , &
118 INTENT(INOUT) :: &
119 TIMEC_KF
121 INTEGER, INTENT(IN) :: KF_EDRATES
123 ! LOCAL VARS
125 REAL, DIMENSION( kts:kte ) :: &
126 U1D, &
127 V1D, &
128 T1D, &
129 DZ1D, &
130 QV1D, &
131 P1D, &
132 RHO1D, &
133 W0AVG1D
135 REAL, DIMENSION( kts:kte ):: &
136 DQDT, &
137 DQIDT, &
138 DQCDT, &
139 DQRDT, &
140 DQSDT, &
141 DTDT
143 REAL :: TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp
145 INTEGER :: i,j,k,NTST,ICLDCK
147 LOGICAL :: qi_flag , qs_flag
148 ! adjustable time step changes
149 REAL :: lastdt = -1.0
150 REAL :: W0AVGfctr, W0fctr, W0den
152 !----------------------------------------------------------------------
154 !--- CALL CUMULUS PARAMETERIZATION
155 !
156 !...TST IS THE NUMBER OF TIME STEPS IN 10 MINUTES...W0AVG IS CLOSE TO A
157 !...RUNNING MEAN VERTICAL VELOCITY...NOTE THAT IF YOU CHANGE TST, IT WIL
158 !...CHANGE THE FREQUENCY OF THE CONVECTIVE INTITIATION CHECK (SEE BELOW)
159 !...NOTE THAT THE ORDERING OF VERTICAL LAYERS MUST BE REVERSED FOR W0AVG
160 !...BECAUSE THE ORDERING IS REVERSED IN KFPARA...
161 !
162 DXSQ=DX*DX
163 qi_flag = .FALSE.
164 qs_flag = .FALSE.
165 IF ( PRESENT( F_QI ) ) qi_flag = f_qi
166 IF ( PRESENT( F_QS ) ) qs_flag = f_qs
168 !----------------------
169 NTST=STEPCU
170 TST=float(NTST*2)
171 !----------------------
172 ! NTST=NINT(1200./(DT*2.))
173 ! TST=float(NTST)
174 ! NTST=NINT(0.5*TST)
175 ! NTST=MAX0(NTST,1)
176 !----------------------
177 ! ICLDCK=MOD(KTAU,NTST)
178 !----------------------
179 ! write(0,*) 'DT = ',DT,' KTAU = ',KTAU,' DX = ',DX
180 ! write(0,*) 'CUDT = ',CUDT
181 ! write(0,*) 'ADAPT_STEP_FLAG = ',ADAPT_STEP_FLAG,' IDS = ',IDS
182 ! write(0,*) 'STEPCU = ',STEPCU,' warm_rain = ',warm_rain
183 ! write(0,*) 'F_QV = ',F_QV,' F_QC = ',F_QV
184 ! write(0,*) 'F_QI = ',F_QI,' F_QS = ',F_QS
185 ! write(0,*) 'F_QR = ',F_QR
186 ! stop
187 if (lastdt < 0) then
188 lastdt = dt
189 endif
191 if (ADAPT_STEP_FLAG) then
192 W0AVGfctr = 2 * MAX(CUDT*60,dt) - dt
193 W0fctr = dt
194 W0den = 2 * MAX(CUDT*60,dt)
195 else
196 W0AVGfctr = (TST-1.)
197 W0fctr = 1.
198 W0den = TST
199 endif
201 DO J = jts,jte
202 DO K=kts,kte
203 DO I= its,ite
204 ! SCR1=-5.0E-4*G*rho(I,K,J)*(w(I,K,J)+w(I,K+1,J))
205 ! TV=T(I,K,J)*(1.+EP1*QV(I,K,J))
206 ! RHOE=Pcps(I,K,J)/(R*TV)
207 ! W0=-101.9368*SCR1/RHOE
208 W0=0.5*(w(I,K,J)+w(I,K+1,J))
210 ! Old:
211 !
212 ! W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST
213 ! New, to support adaptive time step:
214 !
215 W0AVG(I,K,J) = ( W0AVG(I,K,J) * W0AVGfctr + W0 * W0fctr ) / W0den
217 ENDDO
218 ENDDO
219 ENDDO
220 lastdt = dt
222 DO J = jts,jte
223 DO I= its,ite
224 CU_ACT_FLAG(i,j) = .true.
225 ENDDO
226 ENDDO
228 DO J = jts,jte
229 DO I=its,ite
230 ! if (i.eq. 110 .and. j .eq. 59 ) then
231 ! write(0,*) 'nca = ',nca(i,j),' CU_ACT_FLAG = ',CU_ACT_FLAG(i,j)
232 ! write(0,*) 'dt = ',dt,' ADAPT_STEP_FLAG = ',ADAPT_STEP_FLAG
233 ! endif
234 ! IF ( NINT(NCA(I,J)) .gt. 0 ) then
235 IF ( NCA(I,J) .gt. 0.5*DT ) then
236 CU_ACT_FLAG(i,j) = .false.
237 ELSE
239 DO k=kts,kte
240 DQDT(k)=0.
241 DQIDT(k)=0.
242 DQCDT(k)=0.
243 DQRDT(k)=0.
244 DQSDT(k)=0.
245 DTDT(k)=0.
246 ENDDO
247 !kf_edrates
248 IF (KF_EDRATES == 1) THEN
249 DO k=kts,kte
250 UDR_KF(I,k,J)=0.
251 DDR_KF(I,k,J)=0.
252 UER_KF(I,k,J)=0.
253 DER_KF(I,k,J)=0.
254 ENDDO
255 TIMEC_KF(I,J)=0.
256 ENDIF
257 RAINCV(I,J)=0.
258 PRATEC(I,J)=0.
259 !
260 ! assign vars from 3D to 1D
262 DO K=kts,kte
263 U1D(K) =U(I,K,J)
264 V1D(K) =V(I,K,J)
265 T1D(K) =T(I,K,J)
266 RHO1D(K) =rho(I,K,J)
267 QV1D(K)=QV(I,K,J)
268 P1D(K) =Pcps(I,K,J)
269 W0AVG1D(K) =W0AVG(I,K,J)
270 DZ1D(k)=dz8w(I,K,J)
271 ENDDO
273 !
274 CALL KFPARA(I, J, &
275 U1D,V1D,T1D,QV1D,P1D,DZ1D, &
276 W0AVG1D,DT,DX,DXSQ,RHO1D, &
277 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
278 EP2,SVP1,SVP2,SVP3,SVPT0, &
279 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
280 RAINCV,PRATEC,NCA, &
281 warm_rain,qi_flag,qs_flag, &
282 ids,ide, jds,jde, kds,kde, &
283 ims,ime, jms,jme, kms,kme, &
284 its,ite, jts,jte, kts,kte, &
285 !kf_edrates
286 UDR_KF,DDR_KF, &
287 UER_KF,DER_KF, &
288 TIMEC_KF,KF_EDRATES )
290 IF ( PRESENT( RTHCUTEN ) .AND. PRESENT( RQVCUTEN ) ) THEN
291 DO K=kts,kte
292 RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
293 RQVCUTEN(I,K,J)=DQDT(K)
294 ENDDO
295 ENDIF
297 IF( PRESENT(RQRCUTEN) .AND. PRESENT(RQCCUTEN) .AND. &
298 PRESENT(F_QR) ) THEN
299 IF ( F_QR ) THEN
300 DO K=kts,kte
301 RQRCUTEN(I,K,J)=DQRDT(K)
302 RQCCUTEN(I,K,J)=DQCDT(K)
303 ENDDO
304 ELSE
305 ! This is the case for Eta microphysics without 3d rain field
306 DO K=kts,kte
307 RQRCUTEN(I,K,J)=0.
308 RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
309 ENDDO
310 ENDIF
311 ENDIF
313 !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
315 IF( PRESENT( RQICUTEN ) .AND. qi_flag )THEN
316 DO K=kts,kte
317 RQICUTEN(I,K,J)=DQIDT(K)
318 ENDDO
319 ENDIF
321 IF( PRESENT ( RQSCUTEN ) .AND. qs_flag )THEN
322 DO K=kts,kte
323 RQSCUTEN(I,K,J)=DQSDT(K)
324 ENDDO
325 ENDIF
326 !
327 ENDIF
328 ENDDO
329 ENDDO
331 END SUBROUTINE KFCPS
333 !-----------------------------------------------------------
334 SUBROUTINE KFPARA (I, J, &
335 U0,V0,T0,QV0,P0,DZQ,W0AVG1D, &
336 DT,DX,DXSQ,rho, &
337 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
338 EP2,SVP1,SVP2,SVP3,SVPT0, &
339 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
340 RAINCV,PRATEC,NCA, &
341 warm_rain,qi_flag,qs_flag, &
342 ids,ide, jds,jde, kds,kde, &
343 ims,ime, jms,jme, kms,kme, &
344 its,ite, jts,jte, kts,kte, &
345 !kf_edrates
346 UDR_KF,DDR_KF, &
347 UER_KF,DER_KF, &
348 TIMEC_KF,KF_EDRATES )
349 !-----------------------------------------------------------
350 IMPLICIT NONE
351 !-----------------------------------------------------------
352 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
353 ims,ime, jms,jme, kms,kme, &
354 its,ite, jts,jte, kts,kte, &
355 I,J
356 LOGICAL, INTENT(IN ) :: warm_rain
357 LOGICAL :: qi_flag, qs_flag
359 !
360 REAL, DIMENSION( kts:kte ), &
361 INTENT(IN ) :: U0, &
362 V0, &
363 T0, &
364 QV0, &
365 P0, &
366 rho, &
367 DZQ, &
368 W0AVG1D
369 !
370 REAL, INTENT(IN ) :: DT,DX,DXSQ
371 !
373 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
374 REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
375 !
376 REAL, DIMENSION( kts:kte ), INTENT(INOUT) :: &
377 DQDT, &
378 DQIDT, &
379 DQCDT, &
380 DQRDT, &
381 DQSDT, &
382 DTDT
384 REAL, DIMENSION( ims:ime , jms:jme ), &
385 INTENT(INOUT) :: RAINCV, &
386 PRATEC, &
387 NCA
389 !kf_edrates
390 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
391 INTENT(INOUT) :: UDR_KF, &
392 DDR_KF, &
393 UER_KF, &
394 DER_KF
396 REAL, DIMENSION( ims:ime , jms:jme ), &
397 INTENT(INOUT) :: TIMEC_KF
399 INTEGER, INTENT(IN) :: KF_EDRATES
401 !
402 !...DEFINE LOCAL VARIABLES...
403 !
404 REAL, DIMENSION( kts:kte ) :: &
405 Q0,Z0,TV0,TU,TVU,QU,TZ,TVD, &
406 QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD, &
407 UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2, &
408 UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE, &
409 THTAU,THETEU,THTAD,THETED,QLIQ,QICE, &
410 QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC, &
411 DETLQ2,DETIC2,RATIO,RATIO2
413 REAL, DIMENSION( kts:kte ) :: &
414 DOMGDP,EXN,RHOE,TVQU,DP,RH,EQFRC,WSPD, &
415 QDT,FXM,THTAG,THTESG,THPA,THFXTOP, &
416 THFXBOT,QPA,QFXTOP,QFXBOT,QLPA,QLFXIN, &
417 QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA, &
418 QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT, &
419 QL0,QLG,QI0,QIG,QR0,QRG,QS0,QSG
421 REAL, DIMENSION( kts:kte+1 ) :: OMG
422 REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
424 ! LOCAL VARS
426 REAL :: P00,T00,CV,B61,RLF,RHIC,RHBC,PIE, &
427 TTFRZ,TBFRZ,C5,RATE
428 REAL :: GDRY,ROCP,ALIQ,BLIQ, &
429 CLIQ,DLIQ,AICE,BICE,CICE,DICE
430 REAL :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX, &
431 ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL, &
432 CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR, &
433 ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
434 TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD, &
435 UPNEW,ABE,WKLCL,THTUDL,TUDL,TTEMP,FRC1, &
436 QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
437 DZZ,WSQ,UDLBE,REI,EE2,UD2,TTMP,F1,F2, &
438 THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1, &
439 UD1,CLDHGT,DPTT,QNEWLQ,DUMFDP,EE,TSAT, &
440 THTA,P150,USR,VCONV,TIMEC,SHSIGN,VWS,PEF, &
441 CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN, &
442 DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1, &
443 DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF, &
444 UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF, &
445 DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
446 AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1, &
447 DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF, &
448 TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR, &
449 UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2, &
450 RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
451 DDFRC,TDC,DEFRC
453 INTEGER :: KX,K,KL
454 !
455 INTEGER :: ISTOP,ML,L5,L4,KMIX,LOW, &
456 LC,MXLAYR,LLFC,NLAYRS,NK, &
457 KPBL,KLCL,LCL,LET,IFLAG, &
458 KFRZ,NK1,LTOP,NJ,LTOP1, &
459 LTOPM1,LVF,KSTART,KMIN,LFS, &
460 ND,NIC,LDB,LDT,ND1,NDK, &
461 NM,LMAX,NCOUNT,NOITR, &
462 NSTEP,NTC
463 !
464 DATA P00,T00/1.E5,273.16/
465 DATA CV,B61,RLF/717.,0.608,3.339E5/
466 DATA RHIC,RHBC/1.,0.90/
467 DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
468 DATA RATE/0.01/
469 !-----------------------------------------------------------
470 GDRY=-G/CP
471 ROCP=R/CP
472 KL=kte
473 KX=kte
474 !
475 ! ALIQ = 613.3
476 ! BLIQ = 17.502
477 ! CLIQ = 4780.8
478 ! DLIQ = 32.19
479 ALIQ = SVP1*1000.
480 BLIQ = SVP2
481 CLIQ = SVP2*SVPT0
482 DLIQ = SVP3
483 AICE = 613.2
484 BICE = 22.452
485 CICE = 6133.0
486 DICE = 0.61
487 !
489 !...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER
490 !...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL)
491 !...FIELD. 'FBFRC' IS THE FRACTION OF AVAILABLE
492 !...PRECIPITATION TO BE FED BACK (0.0 - 1.0)...
493 !
494 FBFRC=0.0
495 !
496 !...SCHEME IS CALLED ONCE ON EACH NORTH-SOUTH SLICE, THE LOOP BELOW
497 !...CHECKS FOR THE POSSIBILITY OF INITIATING PARAMETERIZED
498 !...CONVECTION AT EACH POINT WITHIN THE SLICE
499 !
500 !...SEE IF IT IS NECESSARY TO CHECK FOR CONVECTIVE TRIGGERING AT THIS
501 !...GRID POINT. IF NCA>0, CONVECTION IS ALREADY ACTIVE AT THIS POINT,
502 !...JUST FEED BACK THE TENDENCIES SAVED FROM THE TIME WHEN CONVECTION
503 !...WAS INITIATED. IF NCA<0, CONVECTION IS NOT ACTIVE
504 !...AND YOU MAY WANT TO CHECK TO SEE IF IT CAN BE ACTIVATED FOR THE
505 !...CURRENT CONDITIONS. IN PREVIOUS APLICATIONS OF THIS SCHEME,
506 !...THE VARIABLE ICLDCK WAS USED BELOW TO SAVE TIME BY ONLY CHECKING
507 !...FOR THE POSSIBILITY OF CONVECTIVE INITIATION EVERY 5 OR 10
508 !...MINUTES...
509 !
511 ! 10 CONTINUE
512 !SUE P300=1000.*(PSB(I,J)*A(KL)+PTOP-30.)+PP3D(I,J,KL)
514 P300=P0(1)-30000.
515 !
516 !...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
517 !...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
518 !...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
519 !
520 !...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
521 !...FROM BOTTOM-UP IN THE KF SCHEME...
522 !
523 ML=0
524 !SUE tmprpsb=1./PSB(I,J)
525 !SUE CELL=PTOP*tmprpsb
527 DO 15 K=1,KX
528 !SUE P0(K)=1.E3*(A(NK)*PSB(I,J)+PTOP)+PP3D(I,J,NK)
529 !
530 !...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
531 !
532 ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
533 QES(K)=EP2*ES/(P0(K)-ES)
534 Q0(K)=AMIN1(QES(K),QV0(K))
535 Q0(K)=AMAX1(0.000001,Q0(K))
536 QL0(K)=0.
537 QI0(K)=0.
538 QR0(K)=0.
539 QS0(K)=0.
541 TV0(K)=T0(K)*(1.+B61*Q0(K))
542 RHOE(K)=P0(K)/(R*TV0(K))
544 DP(K)=rho(k)*g*DZQ(k)
545 !
546 !...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
547 ! DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
548 !
549 IF(P0(K).GE.500E2)L5=K
550 IF(P0(K).GE.400E2)L4=K
551 IF(P0(K).GE.P300)LLFC=K
552 IF(T0(K).GT.T00)ML=K
553 15 CONTINUE
555 Z0(1)=.5*DZQ(1)
556 DO 20 K=2,KL
557 Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
558 DZA(K-1)=Z0(K)-Z0(K-1)
559 20 CONTINUE
560 DZA(KL)=0.
561 KMIX=1
562 25 LOW=KMIX
564 IF(LOW.GT.LLFC)GOTO 325
566 LC=LOW
567 MXLAYR=0
568 !
569 !...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
570 !...UNSTABLE AIR 50 TO 100 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
571 !...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
572 !...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 60 mb..
573 !
574 NLAYRS=0
575 DPTHMX=0.
576 DO 63 NK=LC,KX
577 DPTHMX=DPTHMX+DP(NK)
578 NLAYRS=NLAYRS+1
579 63 IF(DPTHMX.GT.6.E3)GOTO 64
580 GOTO 325
581 64 KPBL=LC+NLAYRS-1
582 KMIX=LC+1
583 18 THMIX=0.
584 QMIX=0.
585 ZMIX=0.
586 PMIX=0.
587 DPTHMX=0.
588 !
589 !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
590 !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
591 !...LAYERS...
592 !
593 DO 17 NK=LC,KPBL
594 DPTHMX=DPTHMX+DP(NK)
595 ROCPQ=0.2854*(1.-0.28*Q0(NK))
596 THMIX=THMIX+DP(NK)*T0(NK)*(P00/P0(NK))**ROCPQ
597 QMIX=QMIX+DP(NK)*Q0(NK)
598 ZMIX=ZMIX+DP(NK)*Z0(NK)
599 17 PMIX=PMIX+DP(NK)*P0(NK)
600 THMIX=THMIX/DPTHMX
601 QMIX=QMIX/DPTHMX
602 ZMIX=ZMIX/DPTHMX
603 PMIX=PMIX/DPTHMX
604 ROCPQ=0.2854*(1.-0.28*QMIX)
605 TMIX=THMIX*(PMIX/P00)**ROCPQ
606 EMIX=QMIX*PMIX/(EP2+QMIX)
607 !
608 !...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL, PRESSURE
609 !...LEVEL OF LCL...
610 !
611 TLOG=ALOG(EMIX/ALIQ)
612 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
613 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX- &
614 TDPT)
615 TLCL=AMIN1(TLCL,TMIX)
616 TVLCL=TLCL*(1.+0.608*QMIX)
617 CPORQ=1./ROCPQ
618 PLCL=P00*(TLCL/THMIX)**CPORQ
619 DO 29 NK=LC,KL
620 KLCL=NK
621 IF(PLCL.GE.P0(NK))GOTO 35
622 29 CONTINUE
623 GOTO 325
624 35 K=KLCL-1
625 DLP=ALOG(PLCL/P0(K))/ALOG(P0(KLCL)/P0(K))
626 !
627 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
628 !
629 TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
630 QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
631 TVEN=TENV*(1.+0.608*QENV)
632 TVBAR=0.5*(TV0(K)+TVEN)
633 ! ZLCL=Z0(K)+R*TVBAR*ALOG(P0(K)/PLCL)/G
634 ZLCL=Z0(K)+(Z0(KLCL)-Z0(K))*DLP
635 !
636 !...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
637 !...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0AVG IS AN
638 !...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
639 !...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
640 !...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
641 !...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
642 !...SUCCESS AT GRID LENGTHS NEAR 25 km. FOR DIFFERENT GRID-LENGTHS,
643 !...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
644 !...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
645 !
646 WKLCL=0.02*ZLCL/2.5E3
647 WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)*DX/25.E3- &
648 WKLCL
649 WABS=ABS(WKL)+1.E-10
650 WSIGNE=WKL/WABS
651 DTLCL=4.64*WSIGNE*WABS**0.33
652 GDT=G*DTLCL*(ZLCL-Z0(LC))/(TV0(LC)+TVEN)
653 WLCL=1.+.5*WSIGNE*SQRT(ABS(GDT)+1.E-10)
654 IF(TLCL+DTLCL.GT.TENV)GOTO 45
655 IF(KPBL.GE.LLFC)GOTO 325
656 GOTO 25
657 !
658 !...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
659 !...EQUIVALENT POTENTIAL TEMPERATURE
660 !...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
661 !
662 45 THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
663 EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
664 ES=ALIQ*EXP((TENV*BLIQ-CLIQ)/(TENV-DLIQ))
665 TVAVG=0.5*(TV0(KLCL)+TENV*(1.+0.608*QENV))
666 PLCL=P0(KLCL)*EXP(G/(R*TVAVG)*(Z0(KLCL)-ZLCL))
667 QESE=EP2*ES/(PLCL-ES)
668 GDT=G*DTLCL*(ZLCL-Z0(LC))/(TV0(LC)+TVEN)
669 WLCL=1.+.5*WSIGNE*SQRT(ABS(GDT)+1.E-10)
670 THTES(K)=TENV*(1.E5/PLCL)**(0.2854*(1.-0.28*QESE))* &
671 EXP((3374.6525/TENV-2.5403)*QESE*(1.+0.81*QESE))
672 WTW=WLCL*WLCL
673 IF(WLCL.LT.0.)GOTO 25
674 TVLCL=TLCL*(1.+0.608*QMIX)
675 RHOLCL=PLCL/(R*TVLCL)
676 !
677 LCL=KLCL
678 LET=LCL
679 !
680 !*******************************************************************
681 ! *
682 ! COMPUTE UPDRAFT PROPERTIES *
683 ! *
684 !*******************************************************************
685 !
686 !
687 !...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
688 !
689 WU(K)=WLCL
690 AU0=PIE*RAD*RAD
691 UMF(K)=RHOLCL*AU0
692 VMFLCL=UMF(K)
693 UPOLD=VMFLCL
694 UPNEW=UPOLD
695 !
696 !...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
697 !...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE BUOYANT ENERGY,
698 ! TRPPT IS THE TOTAL RATE OF PRECIPITATION PRODUCTION...
699 !
700 RATIO2(K)=0.
701 UER(K)=0.
702 ABE=0.
703 TRPPT=0.
704 TU(K)=TLCL
705 TVU(K)=TVLCL
706 QU(K)=QMIX
707 EQFRC(K)=1.
708 QLIQ(K)=0.
709 QICE(K)=0.
710 QLQOUT(K)=0.
711 QICOUT(K)=0.
712 DETLQ(K)=0.
713 DETIC(K)=0.
714 PPTLIQ(K)=0.
715 PPTICE(K)=0.
716 IFLAG=0
717 KFRZ=LC
718 !
719 !...THE AMOUNT OF CONV AVAIL POT ENERGY (CAPE) IS CALCULATED WITH
720 ! RESPECT TO UNDILUTE PARCEL ASCENT; EQ POT TEMP OF UNDILUTE
721 ! PARCEL IS THTUDL, UNDILUTE TEMPERATURE IS GIVEN BY TUDL...
722 !
723 THTUDL=THETEU(K)
724 TUDL=TLCL
725 !
726 !...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
727 ! PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
728 ! FREEZING IS SPECIFIED TO BEGIN. WITHIN THE GLACIATION
729 ! INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
730 ! PREVIOUS MODEL LEVEL...
731 !
732 TTEMP=TTFRZ
733 !
734 !...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
735 ! MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
736 ! MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
737 !
738 DO 60 NK=K,KL-1
739 NK1=NK+1
740 RATIO2(NK1)=RATIO2(NK)
741 !
742 !...UPDATE UPDRAFT PROPERTIES AT THE NEXT MODEL LVL TO REFLECT
743 ! ENTRAINMENT OF ENVIRONMENTAL AIR...
744 !
745 FRC1=0.
746 TU(NK1)=T0(NK1)
747 THETEU(NK1)=THETEU(NK)
748 QU(NK1)=QU(NK)
749 QLIQ(NK1)=QLIQ(NK)
750 QICE(NK1)=QICE(NK)
752 CALL TPMIX(P0(NK1),THETEU(NK1),TU(NK1),QU(NK1),QLIQ(NK1), &
753 QICE(NK1),QNEWLQ,QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0, &
754 XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
755 TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
756 !
757 !...CHECK TO SEE IF UPDRAFT TEMP IS WITHIN THE FREEZING INTERVAL,
758 ! IF IT IS, CALCULATE THE FRACTIONAL CONVERSION TO GLACIATION
759 ! AND ADJUST QNEWLQ TO REFLECT THE GRADUAL CHANGE IN THETAU
760 ! SINCE THE LAST MODEL LEVEL...THE GLACIATION EFFECTS WILL BE
761 ! DETERMINED AFTER THE AMOUNT OF CONDENSATE AVAILABLE AFTER
762 ! PRECIP FALLOUT IS DETERMINED...TTFRZ IS THE TEMP AT WHICH
763 ! GLACIATION BEGINS, TBFRZ THE TEMP AT WHICH IT ENDS...
764 !
765 IF(TU(NK1).LE.TTFRZ.AND.IFLAG.LT.1)THEN
766 IF(TU(NK1).GT.TBFRZ)THEN
767 IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
768 FRC1=(TTEMP-TU(NK1))/(TTFRZ-TBFRZ)
769 R1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
770 ELSE
771 FRC1=(TTEMP-TBFRZ)/(TTFRZ-TBFRZ)
772 R1=1.
773 IFLAG=1
774 ENDIF
775 QNWFRZ=QNEWLQ
776 QNEWIC=QNEWIC+QNEWLQ*R1*0.5
777 QNEWLQ=QNEWLQ-QNEWLQ*R1*0.5
778 EFFQ=(TTFRZ-TBFRZ)/(TTEMP-TBFRZ)
779 TTEMP=TU(NK1)
780 ENDIF
781 !
782 ! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
783 !
784 IF(NK.EQ.K)THEN
785 BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
786 BOTERM=2.*(Z0(NK1)-ZLCL)*G*BE/1.5
787 ENTERM=0.
788 DZZ=Z0(NK1)-ZLCL
789 ELSE
790 BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
791 BOTERM=2.*DZA(NK)*G*BE/1.5
792 ENTERM=2.*UER(NK)*WTW/UPOLD
793 DZZ=DZA(NK)
794 ENDIF
795 WSQ=WTW
796 CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM,RATE, &
797 QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1), G)
799 !...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
800 ! IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
801 !
802 IF(WTW.LE.0.)GOTO 65
803 WABS=SQRT(ABS(WTW))
804 WU(NK1)=WTW/WABS
805 !
806 ! UPDATE THE ABE FOR UNDILUTE ASCENT...
807 !
808 THTES(NK1)=T0(NK1)*(1.E5/P0(NK1))**(0.2854*(1.-0.28*QES(NK1))) &
809 * &
810 EXP((3374.6525/T0(NK1)-2.5403)*QES(NK1)*(1.+0.81* &
811 QES(NK1)))
812 UDLBE=((2.*THTUDL)/(THTES(NK)+THTES(NK1))-1.)*DZZ
813 IF(UDLBE.GT.0.)ABE=ABE+UDLBE*G
814 !
815 ! DETERMINE THE EFFECTS OF CLOUD GLACIATION IF WITHIN THE SPECIFIED
816 ! TEMP INTERVAL...
817 !
818 IF(FRC1.GT.1.E-6)THEN
819 CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QLIQ(NK1), &
820 QICE(NK1),RATIO2(NK1),TTFRZ,TBFRZ,QNWFRZ,RL,FRC1,EFFQ, &
821 IFLAG,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE &
822 ,CICE,DICE)
823 ENDIF
824 !
825 ! CALL SUBROUTINE TO CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMP.
826 ! WITHIN GLACIATION INTERVAL, THETAE MUST BE CALCULATED WITH RESPECT TO
827 ! SAME DEGREE OF GLACIATION FOR ALL ENTRAINING AIR...
828 !
829 CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),RATIO2(NK1), &
830 RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
832 !...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
833 !
834 REI=VMFLCL*DP(NK1)*0.03/RAD
835 TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
836 !
837 !...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, NO
838 ! ENTRAINMENT IS ALLOWED AT THIS LEVEL...
839 !
840 IF(TVQU(NK1).LE.TV0(NK1))THEN
841 UER(NK1)=0.0
842 UDR(NK1)=REI
843 EE2=0.
844 UD2=1.
845 EQFRC(NK1)=0.
846 GOTO 55
847 ENDIF
848 LET=NK1
849 TTMP=TVQU(NK1)
850 !
851 !...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL
852 ! AIR FOR ESTIMATION OF ENTRAINMENT AND DETRAINMENT RATES...
853 !
854 F1=0.95
855 F2=1.-F1
856 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
857 QTMP=F1*Q0(NK1)+F2*QU(NK1)
858 TMPLIQ=F2*QLIQ(NK1)
859 TMPICE=F2*QICE(NK1)
860 CALL TPMIX(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,QNEWLQ, &
861 QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ, &
862 DLIQ,AICE,BICE,CICE,DICE)
863 TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
864 IF(TU95.GT.TV0(NK1))THEN
865 EE2=1.
866 UD2=0.
867 EQFRC(NK1)=1.0
868 GOTO 50
869 ENDIF
870 F1=0.10
871 F2=1.-F1
872 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
873 QTMP=F1*Q0(NK1)+F2*QU(NK1)
874 TMPLIQ=F2*QLIQ(NK1)
875 TMPICE=F2*QICE(NK1)
876 CALL TPMIX(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,QNEWLQ, &
877 QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ, &
878 DLIQ,AICE,BICE,CICE,DICE)
879 TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
880 IF(TU10.EQ.TVQU(NK1))THEN
881 EE2=1.
882 UD2=0.
883 EQFRC(NK1)=1.0
884 GOTO 50
885 ENDIF
886 EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
887 EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
888 EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
889 IF(EQFRC(NK1).EQ.1)THEN
890 EE2=1.
891 UD2=0.
892 GOTO 50
893 ELSEIF(EQFRC(NK1).EQ.0.)THEN
894 EE2=0.
895 UD2=1.
896 GOTO 50
897 ELSE
898 !
899 !...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
900 ! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
901 !
902 CALL PROF5(EQFRC(NK1),EE2,UD2)
903 ENDIF
904 !
905 50 IF(NK.EQ.K)THEN
906 EE1=1.
907 UD1=0.
908 ENDIF
909 !
910 !...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE
911 ! FRACTIONAL VALUES IN THE LAYER...
912 !
913 UER(NK1)=0.5*REI*(EE1+EE2)
914 UDR(NK1)=0.5*REI*(UD1+UD2)
915 !
916 !...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
917 ! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATION
918 !
919 55 IF(UMF(NK)-UDR(NK1).LT.10.)THEN
920 !
921 !...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL
922 ! UPDRAFT FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE
923 ! PREVIOUS MODEL
924 !
925 IF(UDLBE.GT.0.)ABE=ABE-UDLBE*G
926 LET=NK
927 ! WRITE(98,1015)P0(NK1)/100.
928 GOTO 65
929 ENDIF
930 EE1=EE2
931 UD1=UD2
932 UPOLD=UMF(NK)-UDR(NK1)
933 UPNEW=UPOLD+UER(NK1)
934 UMF(NK1)=UPNEW
935 !
936 !...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND ICE IN
937 ! THE DETRAINING UPDRAFT MASS...
938 !
939 DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
940 DETIC(NK1)=QICE(NK1)*UDR(NK1)
941 QDT(NK1)=QU(NK1)
942 QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
943 THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
944 QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
945 QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
946 !
947 !...KFRZ IS THE HIGHEST MODEL LEVEL AT WHICH LIQUID CONDENSATE IS
948 ! GENERATING PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF LIQUID
949 ! PRECIP AT A GIVING MODEL LVL, PPTICE THE SAME FOR ICE, TRPPT IS
950 ! THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE CURRENT MODEL LEVEL
951 !
952 IF(ABS(RATIO2(NK1)-1.).GT.1.E-6)KFRZ=NK1
953 PPTLIQ(NK1)=QLQOUT(NK1)*(UMF(NK)-UDR(NK1))
954 PPTICE(NK1)=QICOUT(NK1)*(UMF(NK)-UDR(NK1))
955 TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
956 IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
957 60 CONTINUE
958 !
959 !...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIU
960 ! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
961 ! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE
962 ! BETWEEN THE LET AND CLOUD TOP...
963 !
964 !...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL
965 ! VELOCITY FIRST BECOMES NEGATIVE...
966 !
967 65 LTOP=NK
968 CLDHGT=Z0(LTOP)-ZLCL
969 !
970 !...IF CLOUD TOP HGT IS LESS THAN SPECIFIED MINIMUM HEIGHT, GO BACK AND
971 ! THE NEXT HIGHEST 60MB LAYER TO SEE IF A BIGGER CLOUD CAN BE OBTAINED
972 ! THAT SOURCE AIR...
973 !
974 ! IF(CLDHGT.LT.4.E3.OR.ABE.LT.1.)THEN
975 IF(CLDHGT.LT.3.E3.OR.ABE.LT.1.)THEN
976 DO 70 NK=K,LTOP
977 UMF(NK)=0.
978 UDR(NK)=0.
979 UER(NK)=0.
980 DETLQ(NK)=0.
981 DETIC(NK)=0.
982 PPTLIQ(NK)=0.
983 70 PPTICE(NK)=0.
984 GOTO 25
985 ENDIF
986 !
987 !...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS
988 ! FLUX THIS LEVEL...
989 !
990 IF(LET.EQ.LTOP)THEN
991 UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
992 DETLQ(LTOP)=QLIQ(LTOP)*UDR(LTOP)*UPNEW/UPOLD
993 DETIC(LTOP)=QICE(LTOP)*UDR(LTOP)*UPNEW/UPOLD
994 TRPPT=TRPPT-(PPTLIQ(LTOP)+PPTICE(LTOP))
995 UER(LTOP)=0.
996 UMF(LTOP)=0.
997 GOTO 85
998 ENDIF
999 !
1000 ! BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
1001 !
1002 DPTT=0.
1003 DO 71 NJ=LET+1,LTOP
1004 71 DPTT=DPTT+DP(NJ)
1005 DUMFDP=UMF(LET)/DPTT
1006 !
1007 !...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
1008 ! RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
1009 ! PTOP
1010 !
1011 DO 75 NK=LET+1,LTOP
1012 UDR(NK)=DP(NK)*DUMFDP
1013 UMF(NK)=UMF(NK-1)-UDR(NK)
1014 DETLQ(NK)=QLIQ(NK)*UDR(NK)
1015 DETIC(NK)=QICE(NK)*UDR(NK)
1016 TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
1017 PPTLIQ(NK)=(UMF(NK-1)-UDR(NK))*QLQOUT(NK)
1018 PPTICE(NK)=(UMF(NK-1)-UDR(NK))*QICOUT(NK)
1019 TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
1020 75 CONTINUE
1021 !
1022 !...SEND UPDRAFT CHARACTERISTICS TO OUTPUT FILES...
1023 !
1024 85 CONTINUE
1025 !
1026 !...EXTEND THE UPDRAFT MASS FLUX PROFILE DOWN TO THE SOURCE LAYER FOR
1027 ! THE UPDRAFT AIR...ALSO, DEFINE THETAE FOR LEVELS BELOW THE LCL...
1028 !
1029 DO 90 NK=1,K
1030 IF(NK.GE.LC)THEN
1031 IF(NK.EQ.LC)THEN
1032 UMF(NK)=VMFLCL*DP(NK)/DPTHMX
1033 UER(NK)=VMFLCL*DP(NK)/DPTHMX
1034 ELSEIF(NK.LE.KPBL)THEN
1035 UER(NK)=VMFLCL*DP(NK)/DPTHMX
1036 UMF(NK)=UMF(NK-1)+UER(NK)
1037 ELSE
1038 UMF(NK)=VMFLCL
1039 UER(NK)=0.
1040 ENDIF
1041 TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
1042 QU(NK)=QMIX
1043 WU(NK)=WLCL
1044 ELSE
1045 TU(NK)=0.
1046 QU(NK)=0.
1047 UMF(NK)=0.
1048 WU(NK)=0.
1049 UER(NK)=0.
1050 ENDIF
1051 UDR(NK)=0.
1052 QDT(NK)=0.
1053 QLIQ(NK)=0.
1054 QICE(NK)=0.
1055 QLQOUT(NK)=0.
1056 QICOUT(NK)=0.
1057 PPTLIQ(NK)=0.
1058 PPTICE(NK)=0.
1059 DETLQ(NK)=0.
1060 DETIC(NK)=0.
1061 RATIO2(NK)=0.
1062 EE=Q0(NK)*P0(NK)/(EP2+Q0(NK))
1063 TLOG=ALOG(EE/ALIQ)
1064 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
1065 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T0(NK)-T00))*( &
1066 T0(NK)-TDPT)
1067 THTA=T0(NK)*(1.E5/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
1068 THETEE(NK)=THTA* &
1069 EXP((3374.6525/TSAT-2.5403)*Q0(NK)*(1.+0.81*Q0(NK)) &
1070 )
1071 THTES(NK)=THTA* &
1072 EXP((3374.6525/T0(NK)-2.5403)*QES(NK)*(1.+0.81* &
1073 QES(NK)))
1074 EQFRC(NK)=1.0
1075 90 CONTINUE
1076 !
1077 LTOP1=LTOP+1
1078 LTOPM1=LTOP-1
1079 !
1080 !...DEFINE VARIABLES ABOVE CLOUD TOP...
1081 !
1082 DO 95 NK=LTOP1,KX
1083 UMF(NK)=0.
1084 UDR(NK)=0.
1085 UER(NK)=0.
1086 QDT(NK)=0.
1087 QLIQ(NK)=0.
1088 QICE(NK)=0.
1089 QLQOUT(NK)=0.
1090 QICOUT(NK)=0.
1091 DETLQ(NK)=0.
1092 DETIC(NK)=0.
1093 PPTLIQ(NK)=0.
1094 PPTICE(NK)=0.
1095 IF(NK.GT.LTOP1)THEN
1096 TU(NK)=0.
1097 QU(NK)=0.
1098 WU(NK)=0.
1099 ENDIF
1100 THTA0(NK)=0.
1101 THTAU(NK)=0.
1102 EMS(NK)=DP(NK)*DXSQ/G
1103 EMSD(NK)=1./EMS(NK)
1104 TG(NK)=T0(NK)
1105 QG(NK)=Q0(NK)
1106 QLG(NK)=0.
1107 QIG(NK)=0.
1108 QRG(NK)=0.
1109 QSG(NK)=0.
1110 95 OMG(NK)=0.
1111 OMG(KL+1)=0.
1112 P150=P0(KLCL)-1.50E4
1113 DO 100 NK=1,LTOP
1114 THTAD(NK)=0.
1115 EMS(NK)=DP(NK)*DXSQ/G
1116 EMSD(NK)=1./EMS(NK)
1117 !
1118 !...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION
1119 ! SCHEME
1120 !
1121 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
1122 THTAU(NK)=TU(NK)*EXN(NK)
1123 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
1124 THTA0(NK)=T0(NK)*EXN(NK)
1125 !
1126 !...LVF IS THE LEVEL AT WHICH MOISTURE FLUX IS ESTIMATED AS THE BASIS
1127 !...FOR PRECIPITATION EFFICIENCY CALCULATIONS...
1128 !
1129 IF(P0(NK).GT.P150)LVF=NK
1130 100 OMG(NK)=0.
1131 LVF=MIN0(LVF,LET)
1132 USR=UMF(LVF+1)*(QU(LVF+1)+QLIQ(LVF+1)+QICE(LVF+1))
1133 USR=AMIN1(USR,TRPPT)
1134 IF(USR.LT.1.E-8)USR=TRPPT
1135 !
1136 ! WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
1137 ! * TMIX-T00,PMIX,QMIX,ABE
1138 ! WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
1139 ! * WLCL,CLDHGT
1140 !
1141 !...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
1142 !...AND MIDTROPOSPHERE IS USED.
1143 !
1144 WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
1145 WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
1146 WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
1147 VCONV=.5*(WSPD(KLCL)+WSPD(L5))
1148 if (VCONV .gt. 0.) then
1149 TIMEC=DX/VCONV
1150 else
1151 TIMEC=3600.
1152 endif
1153 ! TIMEC=DX/VCONV
1154 TADVEC=TIMEC
1155 TIMEC=AMAX1(1800.,TIMEC)
1156 TIMEC=AMIN1(3600.,TIMEC)
1157 NIC=NINT(TIMEC/DT)
1158 TIMEC=FLOAT(NIC)*DT
1159 !
1160 !...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
1161 !
1162 ! SHSIGN = CVMGT(1.,-1.,WSPD(LTOP).GT.WSPD(KLCL))
1163 IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
1164 SHSIGN=1.
1165 ELSE
1166 SHSIGN=-1.
1167 ENDIF
1168 VWS=(U0(LTOP)-U0(KLCL))*(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL))* &
1169 (V0(LTOP)-V0(KLCL))
1170 VWS=1.E3*SHSIGN*SQRT(VWS)/(Z0(LTOP)-Z0(LCL))
1171 PEF=1.591+VWS*(-.639+VWS*(9.53E-2-VWS*4.96E-3))
1172 PEF=AMAX1(PEF,.2)
1173 PEF=AMIN1(PEF,.9)
1174 !
1175 !...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
1176 !
1177 CBH=(ZLCL-Z0(1))*3.281E-3
1178 IF(CBH.LT.3.)THEN
1179 RCBH=.02
1180 ELSE
1181 RCBH=.96729352+CBH*(-.70034167+CBH*(.162179896+CBH*(- &
1182 1.2569798E-2+CBH*(4.2772E-4-CBH*5.44E-6))))
1183 ENDIF
1184 IF(CBH.GT.25)RCBH=2.4
1185 PEFCBH=1./(1.+RCBH)
1186 PEFCBH=AMIN1(PEFCBH,.9)
1187 !
1188 !... MEAN PEF. IS USED TO COMPUTE RAINFALL.
1189 !
1190 PEFF=.5*(PEF+PEFCBH)
1191 PEFF2=PEFF
1192 ! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1193 !
1194 !*****************************************************************
1195 ! *
1196 ! COMPUTE DOWNDRAFT PROPERTIES *
1197 ! *
1198 !*****************************************************************
1199 !
1200 !...LET DOWNDRAFT ORIGINATE AT THE LEVEL OF MINIMUM SATURATION EQUIVALEN
1201 !...POTENTIAL TEMPERATURE (SEQT) IN THE CLOUD LAYER, EXTEND DOWNWARD TO
1202 !...SURFACE, OR TO THE LAYER BELOW CLOUD BASE AT WHICH ENVIR SEQT IS LES
1203 !...THAN MIN SEQT IN THE CLOUD LAYER...LET DOWNDRAFT DETRAIN OVER A LAYE
1204 !...OF SPECIFIED PRESSURE-DEPTH (DPDD)...
1205 !
1206 TDER=0.
1207 KSTART=MAX0(KPBL,KLCL)
1208 THTMIN=THTES(KSTART+1)
1209 KMIN=KSTART+1
1210 DO 104 NK=KSTART+2,LTOP-1
1211 THTMIN=AMIN1(THTMIN,THTES(NK))
1212 IF(THTMIN.EQ.THTES(NK))KMIN=NK
1213 104 CONTINUE
1214 LFS=KMIN
1215 IF(RATIO2(LFS).GT.0.)CALL ENVIRTHT(P0(LFS),T0(LFS),Q0(LFS), &
1216 THETEE(LFS),0.,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1217 EQFRC(LFS)=(THTES(LFS)-THETEU(LFS))/(THETEE(LFS)-THETEU(LFS))
1218 EQFRC(LFS)=AMAX1(EQFRC(LFS),0.)
1219 EQFRC(LFS)=AMIN1(EQFRC(LFS),1.)
1220 THETED(LFS)=THTES(LFS)
1221 !
1222 !...ESTIMATE THE EFFECT OF MELTING PRECIPITATION IN THE DOWNDRAFT...
1223 !
1224 IF(ML.GT.0)THEN
1225 DTMLTD=0.5*(QU(KLCL)-QU(LTOP))*RLF/CP
1226 ELSE
1227 DTMLTD=0.
1228 ENDIF
1229 TZ(LFS)=T0(LFS)-DTMLTD
1230 ES=ALIQ*EXP((TZ(LFS)*BLIQ-CLIQ)/(TZ(LFS)-DLIQ))
1231 QS=EP2*ES/(P0(LFS)-ES)
1232 QD(LFS)=EQFRC(LFS)*Q0(LFS)+(1.-EQFRC(LFS))*QU(LFS)
1233 THTAD(LFS)=TZ(LFS)*(P00/P0(LFS))**(0.2854*(1.-0.28*QD(LFS)))
1234 IF(QD(LFS).GE.QS)THEN
1235 THETED(LFS)=THTAD(LFS)* &
1236 EXP((3374.6525/TZ(LFS)-2.5403)*QS*(1.+0.81*QS))
1237 ELSE
1238 CALL ENVIRTHT(P0(LFS),TZ(LFS),QD(LFS),THETED(LFS),0.,RL,EP2,ALIQ, &
1239 BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1240 ENDIF
1241 DO 107 NK=1,LFS
1242 ND=LFS-NK
1243 IF(THETED(LFS).GT.THTES(ND).OR.ND.EQ.1)THEN
1244 LDB=ND
1245 !
1246 !...IF DOWNDRAFT NEVER BECOMES NEGATIVELY BUOYANT OR IF IT
1247 !...IS SHALLOWER 50 mb, DON'T ALLOW IT TO OCCUR AT ALL...
1248 !
1249 IF(NK.EQ.1.OR.(P0(LDB)-P0(LFS)).LT.50.E2)GOTO 141
1250 ! testing ---- no downdraft
1251 ! GOTO 141
1252 GOTO 110
1253 ENDIF
1254 107 CONTINUE
1255 !
1256 !...ALLOW DOWNDRAFT TO DETRAIN IN A SINGLE LAYER, BUT WITH DOWNDRAFT AIR
1257 !...TYPICALLY FLUSHED UP INTO HIGHER LAYERS AS ALLOWED IN THE TOTAL
1258 !...VERTICAL ADVECTION CALCULATIONS FARTHER DOWN IN THE CODE...
1259 !
1260 110 DPDD=DP(LDB)
1261 LDT=LDB
1262 FRC=1.
1263 DPT=0.
1264 ! DO 115 NK=LDB,LFS
1265 ! DPT=DPT+DP(NK)
1266 ! IF(DPT.GT.DPDD)THEN
1267 ! LDT=NK
1268 ! FRC=(DPDD+DP(NK)-DPT)/DP(NK)
1269 ! GOTO 120
1270 ! ENDIF
1271 ! IF(NK.EQ.LFS-1)THEN
1272 ! LDT=NK
1273 ! FRC=1.
1274 ! DPDD=DPT
1275 ! GOTO 120
1276 ! ENDIF
1277 !115 CONTINUE
1278 120 CONTINUE
1279 !
1280 !...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX..
1281 !
1282 TVD(LFS)=T0(LFS)*(1.+0.608*QES(LFS))
1283 RDD=P0(LFS)/(R*TVD(LFS))
1284 A1=(1.-PEFF)*AU0
1285 DMF(LFS)=-A1*RDD
1286 DER(LFS)=EQFRC(LFS)*DMF(LFS)
1287 DDR(LFS)=0.
1288 DO 140 ND=LFS-1,LDB,-1
1289 ND1=ND+1
1290 IF(ND.LE.LDT)THEN
1291 DER(ND)=0.
1292 DDR(ND)=-DMF(LDT+1)*DP(ND)*FRC/DPDD
1293 DMF(ND)=DMF(ND1)+DDR(ND)
1294 FRC=1.
1295 THETED(ND)=THETED(ND1)
1296 QD(ND)=QD(ND1)
1297 ELSE
1298 DER(ND)=DMF(LFS)*0.03*DP(ND)/RAD
1299 DDR(ND)=0.
1300 DMF(ND)=DMF(ND1)+DER(ND)
1301 IF(RATIO2(ND).GT.0.)CALL ENVIRTHT(P0(ND),T0(ND),Q0(ND), &
1302 THETEE(ND),0.,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1303 THETED(ND)=(THETED(ND1)*DMF(ND1)+THETEE(ND)*DER(ND))/DMF(ND)
1304 QD(ND)=(QD(ND1)*DMF(ND1)+Q0(ND)*DER(ND))/DMF(ND)
1305 ENDIF
1306 140 CONTINUE
1307 TDER=0.
1308 !
1309 !...CALCULATION AN EVAPORATION RATE FOR GIVEN MASS FLUX...
1310 !
1311 DO 135 ND=LDB,LDT
1312 TZ(ND)= &
1313 TPDD(P0(ND),THETED(LDT),T0(ND),QS,QD(ND),1.0,XLV0,XLV1, &
1314 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1315 ES=ALIQ*EXP((TZ(ND)*BLIQ-CLIQ)/(TZ(ND)-DLIQ))
1316 QS=EP2*ES/(P0(ND)-ES)
1317 DSSDT=(CLIQ-BLIQ*DLIQ)/((TZ(ND)-DLIQ)*(TZ(ND)-DLIQ))
1318 RL=XLV0-XLV1*TZ(ND)
1319 DTMP=RL*QS*(1.-RHBC)/(CP+RL*RHBC*QS*DSSDT)
1320 T1RH=TZ(ND)+DTMP
1321 ES=RHBC*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
1322 QSRH=EP2*ES/(P0(ND)-ES)
1323 !
1324 !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
1325 !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
1326 !
1327 IF(QSRH.LT.QD(ND))THEN
1328 QSRH=QD(ND)
1329 ! T1RH=T1+(QS-QSRH)*RL/CP
1330 T1RH=TZ(ND)
1331 ENDIF
1332 TZ(ND)=T1RH
1333 QS=QSRH
1334 TDER=TDER+(QS-QD(ND))*DDR(ND)
1335 QD(ND)=QS
1336 135 THTAD(ND)=TZ(ND)*(P00/P0(ND))**(0.2854*(1.-0.28*QD(ND)))
1337 !
1338 !...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
1339 !...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
1340 !
1341 141 IF(TDER.LT.1.)THEN
1342 ! WRITE(98,3004)I,J
1343 3004 FORMAT(' ','I=',I3,2X,'J=',I3)
1344 PPTFLX=TRPPT
1345 CPR=TRPPT
1346 TDER=0.
1347 CNDTNF=0.
1348 UPDINC=1.
1349 LDB=LFS
1350 DO 117 NDK=1,LTOP
1351 DMF(NDK)=0.
1352 DER(NDK)=0.
1353 DDR(NDK)=0.
1354 THTAD(NDK)=0.
1355 WD(NDK)=0.
1356 TZ(NDK)=0.
1357 117 QD(NDK)=0.
1358 AINCM2=100.
1359 GOTO 165
1360 ENDIF
1361 !
1362 !...ADJUST DOWNDRAFT MASS FLUX SO THAT EVAPORATION RATE IN DOWNDRAFT IS
1363 !...CONSISTENT WITH PRECIPITATION EFFICIENCY RELATIONSHIP...
1364 !
1365 DEVDMF=TDER/DMF(LFS)
1366 PPR=0.
1367 PPTFLX=PEFF*USR
1368 RCED=TRPPT-PPTFLX
1369 !
1370 !...PPR IS THE TOTAL AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE
1371 !...UPDRAFT FROM CLOUD BASE TO THE LFS...UPDRAFT MASS FLUX WILL BE
1372 !...INCREASED UP TO THE LFS TO ACCOUNT FOR UPDRAFT AIR MIXING WITH
1373 !...ENVIRONMENTAL AIR TO THE UPDRAFT, SO PPR WILL INCREASE
1374 !...PROPORTIONATELY...
1375 !
1376 DO 132 NM=KLCL,LFS
1377 132 PPR=PPR+PPTLIQ(NM)+PPTICE(NM)
1378 IF(LFS.GE.KLCL)THEN
1379 DPPTDF=(1.-PEFF)*PPR*(1.-EQFRC(LFS))/UMF(LFS)
1380 ELSE
1381 DPPTDF=0.
1382 ENDIF
1383 !
1384 !...CNDTNF IS THE AMOUNT OF CONDENSATE TRANSFERRED ALONG WITH UPDRAFT
1385 !...MASS THE DOWNDRAFT AT THE LFS...
1386 !
1387 CNDTNF=(QLIQ(LFS)+QICE(LFS))*(1.-EQFRC(LFS))
1388 DMFLFS=RCED/(DEVDMF+DPPTDF+CNDTNF)
1389 IF(DMFLFS.GT.0.)THEN
1390 TDER=0.
1391 GOTO 141
1392 ENDIF
1393 !
1394 !...DDINC IS THE FACTOR BY WHICH TO INCREASE THE FIRST-GUESS DOWNDRAFT
1395 !...MASS FLUX TO SATISFY THE PRECIP EFFICIENCY RELATIONSHIP, UPDINC IS T
1396 !...WHICH TO INCREASE THE UPDRAFT MASS FLUX BELOW THE LFS TO ACCOUNT FOR
1397 !...TRANSFER OF MASS FROM UPDRAFT TO DOWNDRAFT...
1398 !
1399 ! DDINC=DMFLFS/DMF(LFS)
1400 IF(LFS.GE.KLCL)THEN
1401 UPDINC=(UMF(LFS)-(1.-EQFRC(LFS))*DMFLFS)/UMF(LFS)
1402 !
1403 !...LIMIT UPDINC TO LESS THAN OR EQUAL TO 1.5...
1404 !
1405 IF(UPDINC.GT.1.5)THEN
1406 UPDINC=1.5
1407 DMFLFS2=UMF(LFS)*(UPDINC-1.)/(EQFRC(LFS)-1.)
1408 RCED2=DMFLFS2*(DEVDMF+DPPTDF+CNDTNF)
1409 PPTFLX=PPTFLX+(RCED-RCED2)
1410 PEFF2=PPTFLX/USR
1411 RCED=RCED2
1412 DMFLFS=DMFLFS2
1413 ENDIF
1414 ELSE
1415 UPDINC=1.
1416 ENDIF
1417 DDINC=DMFLFS/DMF(LFS)
1418 DO 149 NK=LDB,LFS
1419 DMF(NK)=DMF(NK)*DDINC
1420 DER(NK)=DER(NK)*DDINC
1421 DDR(NK)=DDR(NK)*DDINC
1422 149 CONTINUE
1423 CPR=TRPPT+PPR*(UPDINC-1.)
1424 PPTFLX=PPTFLX+PEFF*PPR*(UPDINC-1.)
1425 PEFF=PEFF2
1426 TDER=TDER*DDINC
1427 !
1428 !...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
1429 ! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
1430 ! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
1431 !
1432 DO 155 NK=LC,LFS
1433 UMF(NK)=UMF(NK)*UPDINC
1434 UDR(NK)=UDR(NK)*UPDINC
1435 UER(NK)=UER(NK)*UPDINC
1436 PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
1437 PPTICE(NK)=PPTICE(NK)*UPDINC
1438 DETLQ(NK)=DETLQ(NK)*UPDINC
1439 155 DETIC(NK)=DETIC(NK)*UPDINC
1440 !
1441 !...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
1442 !...DOWNDRAFT...
1443 !
1444 IF(LDB.GT.1)THEN
1445 DO 156 NK=1,LDB-1
1446 DMF(NK)=0.
1447 DER(NK)=0.
1448 DDR(NK)=0.
1449 WD(NK)=0.
1450 TZ(NK)=0.
1451 QD(NK)=0.
1452 THTAD(NK)=0.
1453 156 CONTINUE
1454 ENDIF
1455 DO 157 NK=LFS+1,KX
1456 DMF(NK)=0.
1457 DER(NK)=0.
1458 DDR(NK)=0.
1459 WD(NK)=0.
1460 TZ(NK)=0.
1461 QD(NK)=0.
1462 THTAD(NK)=0.
1463 157 CONTINUE
1464 DO 158 NK=LDT+1,LFS-1
1465 TZ(NK)=0.
1466 QD(NK)=0.
1467 158 CONTINUE
1468 !
1469 !
1470 !...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE
1471 ! INFLOW INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN
1472 ! IS AVAILABLE IN THAT LAYER INITIALLY...
1473 !
1474 165 AINCMX=1000.
1475 LMAX=MAX0(KLCL,LFS)
1476 DO 166 NK=LC,LMAX
1477 IF((UER(NK)-DER(NK)).GT.0.)AINCM1=EMS(NK)/((UER(NK)-DER(NK))* &
1478 TIMEC)
1479 AINCMX=AMIN1(AINCMX,AINCM1)
1480 166 CONTINUE
1481 AINC=1.
1482 IF(AINCMX.LT.AINC)AINC=AINCMX
1483 !
1484 !...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRFT AND DOWNDRFT...THEY
1485 !...WILL ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE
1486 !...STABILIZATION CLOSURE...
1487 !
1488 NCOUNT=0
1489 TDER2=TDER
1490 PPTFL2=PPTFLX
1491 DO 170 NK=1,LTOP
1492 DETLQ2(NK)=DETLQ(NK)
1493 DETIC2(NK)=DETIC(NK)
1494 UDR2(NK)=UDR(NK)
1495 UER2(NK)=UER(NK)
1496 DDR2(NK)=DDR(NK)
1497 DER2(NK)=DER(NK)
1498 UMF2(NK)=UMF(NK)
1499 DMF2(NK)=DMF(NK)
1500 170 CONTINUE
1501 FABE=1.
1502 STAB=0.95
1503 NOITR=0
1504 IF(AINC/AINCMX.GT.0.999)THEN
1505 NCOUNT=0
1506 GOTO 255
1507 ENDIF
1508 ISTOP=0
1509 175 NCOUNT=NCOUNT+1
1510 !
1511 !*****************************************************************
1512 ! *
1513 ! COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE *
1514 ! *
1515 !*****************************************************************
1516 !
1517 !...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
1518 !...SATISFY MASS CONTINUITY...
1519 !
1520 185 CONTINUE
1521 DTT=TIMEC
1522 DO 200 NK=1,LTOP
1523 DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
1524 IF(NK.GT.1)THEN
1525 OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
1526 DTT1=0.75*DP(NK-1)/(ABS(OMG(NK))+1.E-10)
1527 DTT=AMIN1(DTT,DTT1)
1528 ENDIF
1529 200 CONTINUE
1530 DO 488 NK=1,LTOP
1531 THPA(NK)=THTA0(NK)
1532 QPA(NK)=Q0(NK)
1533 NSTEP=NINT(TIMEC/DTT+1)
1534 DTIME=TIMEC/FLOAT(NSTEP)
1535 FXM(NK)=OMG(NK)*DXSQ/G
1536 488 CONTINUE
1537 !
1538 !...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
1539 !
1540 DO 495 NTC=1,NSTEP
1541 !
1542 !...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED
1543 !...SIGN OF OMEGA...
1544 !
1545 DO 493 NK=1,LTOP
1546 THFXTOP(NK)=0.
1547 THFXBOT(NK)=0.
1548 QFXTOP(NK)=0.
1549 493 QFXBOT(NK)=0.
1550 DO 494 NK=2,LTOP
1551 IF(OMG(NK).LE.0.)THEN
1552 THFXBOT(NK)=-FXM(NK)*THPA(NK-1)
1553 QFXBOT(NK)=-FXM(NK)*QPA(NK-1)
1554 THFXTOP(NK-1)=THFXTOP(NK-1)-THFXBOT(NK)
1555 QFXTOP(NK-1)=QFXTOP(NK-1)-QFXBOT(NK)
1556 ELSE
1557 THFXBOT(NK)=-FXM(NK)*THPA(NK)
1558 QFXBOT(NK)=-FXM(NK)*QPA(NK)
1559 THFXTOP(NK-1)=THFXTOP(NK-1)-THFXBOT(NK)
1560 QFXTOP(NK-1)=QFXTOP(NK-1)-QFXBOT(NK)
1561 ENDIF
1562 494 CONTINUE
1563 !
1564 !...UPDATE THE THETA AND QV VALUES AT EACH LEVEL..
1565 !
1566 DO 492 NK=1,LTOP
1567 THPA(NK)=THPA(NK)+(THFXBOT(NK)+UDR(NK)*THTAU(NK)+DDR(NK)* &
1568 THTAD(NK)+THFXTOP(NK)-(UER(NK)-DER(NK))*THTA0(NK))* &
1569 DTIME*EMSD(NK)
1570 QPA(NK)=QPA(NK)+(QFXBOT(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK)+ &
1571 QFXTOP(NK)-(UER(NK)-DER(NK))*Q0(NK))*DTIME*EMSD(NK)
1573 492 CONTINUE
1574 495 CONTINUE
1575 DO 498 NK=1,LTOP
1576 THTAG(NK)=THPA(NK)
1577 QG(NK)=QPA(NK)
1578 498 CONTINUE
1579 !
1580 !...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO,
1581 !...BORROW MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO.
1582 !
1583 DO 499 NK=1,LTOP
1584 IF(QG(NK).LT.0.)THEN
1585 IF(NK.EQ.1)THEN
1586 CALL wrf_error_fatal ( 'module_cu_kf.F: problem with kf scheme: qg = 0 at the surface' )
1587 ENDIF
1588 NK1=NK+1
1589 IF(NK.EQ.LTOP)NK1=KLCL
1590 TMA=QG(NK1)*EMS(NK1)
1591 TMB=QG(NK-1)*EMS(NK-1)
1592 TMM=(QG(NK)-1.E-9)*EMS(NK)
1593 BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
1594 ACOEFF=BCOEFF*TMA/TMB
1595 TMB=TMB*(1.-BCOEFF)
1596 TMA=TMA*(1.-ACOEFF)
1597 IF(NK.EQ.LTOP)THEN
1598 QVDIFF=(QG(NK1)-TMA*EMSD(NK1))*100./QG(NK1)
1599 IF(ABS(QVDIFF).GT.1.)THEN
1600 PRINT *,'--WARNING-- CLOUD BASE WATER VAPOR CHANGES BY ', &
1601 QVDIFF, &
1602 ' PERCENT WHEN MOISTURE IS BORROWED TO PREVENT NEG VALUES', &
1603 ' IN KAIN-FRITSCH'
1604 ENDIF
1605 ENDIF
1606 QG(NK)=1.E-9
1607 QG(NK1)=TMA*EMSD(NK1)
1608 QG(NK-1)=TMB*EMSD(NK-1)
1609 ENDIF
1610 499 CONTINUE
1611 TOPOMG=(UDR(LTOP)-UER(LTOP))*DP(LTOP)*EMSD(LTOP)
1612 IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
1613 ! WRITE(98,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME;'
1614 ! * ,'TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1615 WRITE(6,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME;' &
1616 ,'TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1617 ISTOP=1
1618 GOTO 265
1619 ENDIF
1620 !
1621 !...CONVERT THETA TO T...
1622 !
1623 ! PAY ATTENTION ...
1624 !
1625 DO 230 NK=1,LTOP
1626 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QG(NK)))
1627 TG(NK)=THTAG(NK)/EXN(NK)
1628 TVG(NK)=TG(NK)*(1.+0.608*QG(NK))
1629 230 CONTINUE
1630 !
1631 !*******************************************************************
1632 ! *
1633 ! COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY. *
1634 ! *
1635 !*******************************************************************
1636 !
1637 !...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
1638 !
1639 THMIX=0.
1640 QMIX=0.
1641 PMIX=0.
1642 DO 217 NK=LC,KPBL
1643 ROCPQ=0.2854*(1.-0.28*QG(NK))
1644 THMIX=THMIX+DP(NK)*TG(NK)*(P00/P0(NK))**ROCPQ
1645 QMIX=QMIX+DP(NK)*QG(NK)
1646 217 PMIX=PMIX+DP(NK)*P0(NK)
1647 THMIX=THMIX/DPTHMX
1648 QMIX=QMIX/DPTHMX
1649 PMIX=PMIX/DPTHMX
1650 ROCPQ=0.2854*(1.-0.28*QMIX)
1651 TMIX=THMIX*(PMIX/P00)**ROCPQ
1652 ES=ALIQ*EXP((TMIX*BLIQ-CLIQ)/(TMIX-DLIQ))
1653 QS=EP2*ES/(PMIX-ES)
1654 !
1655 !...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
1656 !
1657 IF(QMIX.GT.QS)THEN
1658 RL=XLV0-XLV1*TMIX
1659 CPM=CP*(1.+0.887*QMIX)
1660 DSSDT=QS*(CLIQ-BLIQ*DLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
1661 DQ=(QMIX-QS)/(1.+RL*DSSDT/CPM)
1662 TMIX=TMIX+RL/CP*DQ
1663 QMIX=QMIX-DQ
1664 ROCPQ=0.2854*(1.-0.28*QMIX)
1665 THMIX=TMIX*(P00/PMIX)**ROCPQ
1666 TLCL=TMIX
1667 PLCL=PMIX
1668 ELSE
1669 QMIX=AMAX1(QMIX,0.)
1670 EMIX=QMIX*PMIX/(EP2+QMIX)
1671 TLOG=ALOG(EMIX/ALIQ)
1672 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
1673 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX- &
1674 TDPT)
1675 TLCL=AMIN1(TLCL,TMIX)
1676 CPORQ=1./ROCPQ
1677 PLCL=P00*(TLCL/THMIX)**CPORQ
1678 ENDIF
1679 TVLCL=TLCL*(1.+0.608*QMIX)
1680 DO 235 NK=LC,KL
1681 KLCL=NK
1682 235 IF(PLCL.GE.P0(NK))GOTO 240
1683 240 K=KLCL-1
1684 DLP=ALOG(PLCL/P0(K))/ALOG(P0(KLCL)/P0(K))
1685 !
1686 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
1687 !
1688 TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
1689 QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
1690 TVEN=TENV*(1.+0.608*QENV)
1691 TVBAR=0.5*(TVG(K)+TVEN)
1692 ! ZLCL=Z0(K)+R*TVBAR*ALOG(P0(K)/PLCL)/G
1693 ZLCL=Z0(K)+(Z0(KLCL)-Z0(K))*DLP
1694 TVAVG=0.5*(TVEN+TG(KLCL)*(1.+0.608*QG(KLCL)))
1695 PLCL=P0(KLCL)*EXP(G/(R*TVAVG)*(Z0(KLCL)-ZLCL))
1696 THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
1697 EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
1698 ES=ALIQ*EXP((TENV*BLIQ-CLIQ)/(TENV-DLIQ))
1699 QESE=EP2*ES/(PLCL-ES)
1700 THTESG(K)=TENV*(1.E5/PLCL)**(0.2854*(1.-0.28*QESE))* &
1701 EXP((3374.6525/TENV-2.5403)*QESE*(1.+0.81*QESE))
1702 !
1703 !...COMPUTE ADJUSTED ABE(ABEG).
1704 !
1705 ABEG=0.
1706 THTUDL=THETEU(K)
1707 DO 245 NK=K,LTOPM1
1708 NK1=NK+1
1709 ES=ALIQ*EXP((TG(NK1)*BLIQ-CLIQ)/(TG(NK1)-DLIQ))
1710 QESE=EP2*ES/(P0(NK1)-ES)
1711 THTESG(NK1)=TG(NK1)*(1.E5/P0(NK1))**(0.2854*(1.-0.28*QESE))* &
1712 EXP((3374.6525/TG(NK1)-2.5403)*QESE*(1.+0.81*QESE) &
1713 )
1714 ! DZZ=CVMGT(Z0(KLCL)-ZLCL,DZA(NK),NK.EQ.K)
1715 IF(NK.EQ.K)THEN
1716 DZZ=Z0(KLCL)-ZLCL
1717 ELSE
1718 DZZ=DZA(NK)
1719 ENDIF
1720 BE=((2.*THTUDL)/(THTESG(NK1)+THTESG(NK))-1.)*DZZ
1721 245 IF(BE.GT.0.)ABEG=ABEG+BE*G
1722 !
1723 !...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
1724 !...THE PERIOD TIMEC...
1725 !
1726 IF(NOITR.EQ.1)THEN
1727 ! WRITE(98,1060)FABE
1728 GOTO 265
1729 ENDIF
1730 DABE=AMAX1(ABE-ABEG,0.1*ABE)
1731 FABE=ABEG/(ABE+1.E-8)
1732 IF(FABE.GT.1.)THEN
1733 ! WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS '
1734 ! *,'GRID POINT; NO CONVECTION ALLOWED!'
1735 GOTO 325
1736 ENDIF
1737 IF(NCOUNT.NE.1)THEN
1738 DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
1739 IF(DFDA.GT.0.)THEN
1740 NOITR=1
1741 AINC=AINCOLD
1742 GOTO 255
1743 ENDIF
1744 ENDIF
1745 AINCOLD=AINC
1746 FABEOLD=FABE
1747 IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
1748 ! WRITE(98,1055)FABE
1749 GOTO 265
1750 ENDIF
1751 IF(FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB)GOTO 265
1752 IF(NCOUNT.GT.10)THEN
1753 ! WRITE(98,1060)FABE
1754 GOTO 265
1755 ENDIF
1756 !
1757 !...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE
1758 !...CONVECTIVE MASS FLUX BY THE FACTOR AINC:
1759 !
1760 IF(FABE.EQ.0.)THEN
1761 AINC=AINC*0.5
1762 ELSE
1763 AINC=AINC*STAB*ABE/(DABE+1.E-8)
1764 ENDIF
1765 255 AINC=AMIN1(AINCMX,AINC)
1766 !...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION
1767 !...WILL BE MINIMAL SO JUST IGNORE IT...
1768 IF(AINC.LT.0.05)GOTO 325
1769 ! AINC=AMAX1(AINC,0.05)
1770 TDER=TDER2*AINC
1771 PPTFLX=PPTFL2*AINC
1772 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABEOLD,AINCOLD
1773 DO 260 NK=1,LTOP
1774 UMF(NK)=UMF2(NK)*AINC
1775 DMF(NK)=DMF2(NK)*AINC
1776 DETLQ(NK)=DETLQ2(NK)*AINC
1777 DETIC(NK)=DETIC2(NK)*AINC
1778 UDR(NK)=UDR2(NK)*AINC
1779 UER(NK)=UER2(NK)*AINC
1780 DER(NK)=DER2(NK)*AINC
1781 DDR(NK)=DDR2(NK)*AINC
1782 260 CONTINUE
1783 !
1784 !...GO BACK UP FOR ANOTHER ITERATION...
1785 !
1786 GOTO 175
1787 265 CONTINUE
1789 !kf_edrates
1790 !Save up/down entrainment/detrainment rates as 3D variables
1791 IF (KF_EDRATES == 1) THEN
1792 DO NK=1,LTOP
1793 UDR_KF(I,NK,J)=UDR(NK)
1794 DDR_KF(I,NK,J)=DDR(NK)
1795 UER_KF(I,NK,J)=UER(NK)
1796 DER_KF(I,NK,J)=DER(NK)
1797 ENDDO
1798 ENDIF
1800 !
1801 !...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
1802 !...GRID POINT...
1803 !
1804 !...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
1805 !
1806 !...FRC2 IS THE FRACTION OF TOTAL CONDENSATE
1807 !...GENERATED THAT GOES INTO PRECIPITIATION
1808 FRC2=PPTFLX/(CPR*AINC)
1809 DO 270 NK=1,LTOP
1810 QLPA(NK)=QL0(NK)
1811 QIPA(NK)=QI0(NK)
1812 QRPA(NK)=QR0(NK)
1813 QSPA(NK)=QS0(NK)
1814 RAINFB(NK)=PPTLIQ(NK)*AINC*FBFRC*FRC2
1815 SNOWFB(NK)=PPTICE(NK)*AINC*FBFRC*FRC2
1816 270 CONTINUE
1817 DO 290 NTC=1,NSTEP
1818 !
1819 !...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH
1820 !...LAYER BASED ON THE SIGN OF OMEGA...
1821 !
1822 DO 275 NK=1,LTOP
1823 QLFXIN(NK)=0.
1824 QLFXOUT(NK)=0.
1825 QIFXIN(NK)=0.
1826 QIFXOUT(NK)=0.
1827 QRFXIN(NK)=0.
1828 QRFXOUT(NK)=0.
1829 QSFXIN(NK)=0.
1830 QSFXOUT(NK)=0.
1831 275 CONTINUE
1832 DO 280 NK=2,LTOP
1833 IF(OMG(NK).LE.0.)THEN
1834 QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
1835 QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
1836 QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
1837 QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
1838 QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
1839 QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
1840 QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
1841 QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
1842 ELSE
1843 QLFXOUT(NK)=FXM(NK)*QLPA(NK)
1844 QIFXOUT(NK)=FXM(NK)*QIPA(NK)
1845 QRFXOUT(NK)=FXM(NK)*QRPA(NK)
1846 QSFXOUT(NK)=FXM(NK)*QSPA(NK)
1847 QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
1848 QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
1849 QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
1850 QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
1851 ENDIF
1852 280 CONTINUE
1853 !
1854 !...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
1855 !
1856 DO 285 NK=1,LTOP
1857 QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))*DTIME* &
1858 EMSD(NK)
1859 QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))*DTIME* &
1860 EMSD(NK)
1861 QRPA(NK)=QRPA(NK)+(QRFXIN(NK)+QLQOUT(NK)*UDR(NK)-QRFXOUT(NK) &
1862 +RAINFB(NK))*DTIME*EMSD(NK)
1863 QSPA(NK)=QSPA(NK)+(QSFXIN(NK)+QICOUT(NK)*UDR(NK)-QSFXOUT(NK) &
1864 +SNOWFB(NK))*DTIME*EMSD(NK)
1865 285 CONTINUE
1866 290 CONTINUE
1867 DO 295 NK=1,LTOP
1868 QLG(NK)=QLPA(NK)
1869 QIG(NK)=QIPA(NK)
1870 QRG(NK)=QRPA(NK)
1871 QSG(NK)=QSPA(NK)
1872 295 CONTINUE
1874 !kf_edrates
1875 !Save convective timescale (TIMEC) as 2D variable
1876 IF (KF_EDRATES == 1) THEN
1877 TIMEC_KF(I,J)=TIMEC
1878 ENDIF
1880 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABE,AINC
1881 !
1882 !...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
1883 !
1884 IF(ISTOP.EQ.1)THEN
1885 WRITE(6,1070)' P ',' DP ',' DT K/D ',' DR K/D ',' OMG ', &
1886 ' DOMGDP ',' UMF ',' UER ',' UDR ',' DMF ',' DER ' &
1887 ,' DDR ',' EMS ',' W0 ',' DETLQ ',' DETIC '
1888 DO 300 K=LTOP,1,-1
1889 DTT=(TG(K)-T0(K))*86400./TIMEC
1890 RL=XLV0-XLV1*TG(K)
1891 DR=-(QG(K)-Q0(K))*RL*86400./(TIMEC*CP)
1892 UDFRC=UDR(K)*TIMEC*EMSD(K)
1893 UEFRC=UER(K)*TIMEC*EMSD(K)
1894 DDFRC=DDR(K)*TIMEC*EMSD(K)
1895 DEFRC=-DER(K)*TIMEC*EMSD(K)
1896 WRITE (6,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)* &
1897 1.E4,UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC &
1898 ,DDFRC,EMS(K)/1.E11,W0AVG1D(K)*1.E2,DETLQ(K) &
1899 *TIMEC*EMSD(K)*1.E3,DETIC(K)*TIMEC*EMSD(K)* &
1900 1.E3
1901 300 CONTINUE
1902 WRITE(6,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG', &
1903 'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
1904 DO 305 K=KX,1,-1
1905 DTT=TG(K)-T0(K)
1906 TUC=TU(K)-T00
1907 IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
1908 TDC=TZ(K)-T00
1909 IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
1910 ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
1911 QGS=ES*EP2/(P0(K)-ES)
1912 RH0=Q0(K)/QES(K)
1913 RHG=QG(K)/QGS
1914 WRITE (6,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC &
1915 ,TDC,Q0(K)*1000.,QG(K)*1000.,(QG(K)-Q0(K))* &
1916 1000.,QU(K)*1000.,QD(K)*1000.,QLG(K)*1000., &
1917 QIG(K)*1000.,QRG(K)*1000.,QSG(K)*1000.,RH0,RHG
1918 305 CONTINUE
1919 !
1920 !...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
1921 !...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
1922 !
1923 IF(ISTOP.EQ.1)THEN
1924 DO 310 K=1,KX
1925 WRITE ( wrf_err_message , 1115 ) &
1926 Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
1927 U0(K),V0(K),DP(K)/100.,W0AVG1D(K)
1928 CALL wrf_message ( TRIM( wrf_err_message ) )
1929 310 CONTINUE
1930 CALL wrf_error_fatal ( 'module_cu_kf.F: KAIN-FRITSCH' )
1931 ENDIF
1932 ENDIF
1933 CNDTNF=(1.-EQFRC(LFS))*(QLIQ(LFS)+QICE(LFS))*DMF(LFS)
1934 ! WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
1935 !
1936 ! EVALUATE MOISTURE BUDGET...
1937 !
1938 QINIT=0.
1939 QFNL=0.
1940 DPT=0.
1941 DO 315 NK=1,LTOP
1942 DPT=DPT+DP(NK)
1943 QINIT=QINIT+Q0(NK)*EMS(NK)
1944 QFNL=QFNL+QG(NK)*EMS(NK)
1945 QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
1946 315 CONTINUE
1947 QFNL=QFNL+PPTFLX*TIMEC*(1.-FBFRC)
1948 ERR2=(QFNL-QINIT)*100./QINIT
1949 ! WRITE(98,1110)QINIT,QFNL,ERR2
1950 ! IF(ABS(ERR2).GT.0.05)STOP 'QVERR'
1951 IF(ABS(ERR2).GT.0.05)CALL wrf_error_fatal( 'module_cu_kf.F: QVERR' )
1952 RELERR=ERR2*QINIT/(PPTFLX*TIMEC+1.E-10)
1953 ! WRITE(98,1120)RELERR
1954 ! WRITE(98,*)'TDER, CPR, USR, TRPPT =',
1955 ! *TDER,CPR*AINC,USR*AINC,TRPPT*AINC
1956 !
1957 !...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
1958 !
1959 !...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
1960 !...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
1961 !
1962 IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
1963 NCA(I,J)=FLOAT(NIC)*DT
1964 DO 320 K=1,KX
1965 ! IF(IMOIST.NE.2)THEN
1966 !
1967 !...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMAT
1968 !...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
1969 !...NOTE: THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
1970 !...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
1971 !...OF QG...
1972 !
1973 ! RLC=XLV0-XLV1*TG(K)
1974 ! RLS=XLS0-XLS1*TG(K)
1975 ! CPM=CP*(1.+0.887*QG(K))
1976 ! TG(K)=TG(K)-(RLC*(QLG(K)+QRG(K))+RLS*(QIG(K)+QSG(K)))/CPM
1977 ! QG(K)=QG(K)+(QLG(K)+QRG(K)+QIG(K)+QSG(K))
1978 ! DQCDT(K)=0.
1979 ! DQIDT(K)=0.
1980 ! DQRDT(K)=0.
1981 ! DQSDT(K)=0.
1982 ! ELSE
1983 IF(.NOT. qi_flag .and. warm_rain)THEN
1984 !
1985 !...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
1986 !
1987 CPM=CP*(1.+0.887*QG(K))
1988 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
1989 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
1990 DQIDT(K)=0.
1991 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
1992 DQSDT(K)=0.
1993 ELSEIF(.NOT. qi_flag .and. .not. warm_rain)THEN
1994 !
1995 !...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROME
1996 !...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
1997 !
1998 CPM=CP*(1.+0.887*QG(K))
1999 IF(K.LE.ML)THEN
2000 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
2001 ELSEIF(K.GT.ML)THEN
2002 TG(K)=TG(K)+(QLG(K)+QRG(K))*RLF/CPM
2003 ENDIF
2004 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
2005 DQIDT(K)=0.
2006 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
2007 DQSDT(K)=0.
2008 ELSEIF(qi_flag) THEN
2009 !
2010 !...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE
2011 !...TENDENCY OF HYDROMETEORS DIRECTLY...
2012 !
2013 DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
2014 DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
2015 DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
2016 IF (qs_flag ) THEN
2017 DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
2018 ELSE
2019 DQIDT(K)=DQIDT(K)+(QSG(K)-QS0(K))/TIMEC
2020 ENDIF
2021 ELSE
2022 CALL wrf_error_fatal ( 'module_cu_kf: THIS COMBINATION OF IMOIST, IICE NOT ALLOWED' )
2023 ENDIF
2024 ! ENDIF
2025 DTDT(K)=(TG(K)-T0(K))/TIMEC
2026 DQDT(K)=(QG(K)-Q0(K))/TIMEC
2027 320 CONTINUE
2029 ! RAINCV is in the unit of mm
2031 PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
2032 RAINCV(I,J)=DT*PRATEC(I,J)
2033 RNC=RAINCV(I,J)*NIC
2034 ! WRITE(98,909)RNC
2035 909 FORMAT(' CONVECTIVE RAINFALL =',F8.4,' CM')
2037 325 CONTINUE
2039 1000 FORMAT(' ',10A8)
2040 1005 FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
2041 1010 FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
2042 1015 FORMAT(' ','ALL REMAINING MASS DETRAINS BELOW ',F4.0,' MB')
2043 1025 FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M', &
2044 ' DTLCL=',F5.2,' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=', &
2045 I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1, &
2046 ' CAPE=',0PF7.1)
2047 1030 FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =', &
2048 E12.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =', &
2049 F8.1)
2050 1035 FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'WKL=' &
2051 ,F6.3,'VWS=',F5.2)
2052 1040 FORMAT(' ','PRECIP EFF = 100%, ENVIR CANNOT SUPPORT DOWND' &
2053 ,'RAFTS')
2054 !1045 FORMAT('NUMBER OF DOWNDRAFT ITERATIONS EXCEEDS 10...PPTFLX' &
2055 ! ' IS DIFFERENT FROM THAT GIVEN BY PRECIP EFF RELATION')
2056 ! FLIC HAS TROUBLE WITH THIS ONE.
2057 1045 FORMAT('NUMBER OF DOWNDRAFT ITERATIONS EXCEEDS 10')
2058 1050 FORMAT(' ','LCOUNT= ',I3,' PPTFLX/CPR, PEFF= ',F5.3,1X,F5.3, &
2059 'DMF(LFS)/UMF(LCL)= ',F5.3)
2060 1055 FORMAT(/'*** DEGREE OF STABILIZATION =',F5.3,', NO MORE MASS F' &
2061 ,'LUX IS ALLOWED')
2062 !1060 FORMAT(/' ITERATION DOES NOT CONVERGE TO GIVE THE SPECIFIED ' &
2063 ! 'DEGREE OF STABILIZATION! FABE= ',F6.4)
2064 1060 FORMAT(/' ITERATION DOES NOT CONVERGE. FABE= ',F6.4)
2065 1070 FORMAT (16A8)
2066 1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3)
2067 1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, NSTEP=',F5.0,I3, &
2068 'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
2069 1085 FORMAT (A3,16A7,2A8)
2070 1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3)
2071 1095 FORMAT(' ',' PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ', &
2072 F10.0)
2073 1105 FORMAT(' ','NET LATENT HEAT RELEASE =',E12.5,' ACTUAL HEATING =', &
2074 E12.5,' J/KG-S, DIFFERENCE = ',F9.3,'PERCENT')
2075 1110 FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5, &
2076 ' TOTAL WATER CHANGE =',F8.2,'PERCENT')
2077 1115 FORMAT (2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4 &
2078 )
2079 1120 FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3, &
2080 'PERCENT')
2082 END SUBROUTINE KFPARA
2084 !-----------------------------------------------------------------------
2085 SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
2086 QNEWIC,QLQOUT,QICOUT,G)
2087 !-----------------------------------------------------------------------
2088 IMPLICIT NONE
2089 !-----------------------------------------------------------------------
2090 ! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
2091 ! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
2092 ! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
2093 ! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
2094 ! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
2096 REAL, INTENT(IN ) :: G
2097 REAL, INTENT(IN ) :: DZ,BOTERM,ENTERM,RATE
2098 REAL, INTENT(INOUT) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
2100 REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
2102 QTOT=QLIQ+QICE
2103 QNEW=QNEWLQ+QNEWIC
2104 !
2105 ! ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY C
2106 ! BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
2107 ! LEVELS...
2108 !
2109 QEST=0.5*(QTOT+QNEW)
2110 G1=WTW+BOTERM-ENTERM-2.*G*DZ*QEST/1.5
2111 IF(G1.LT.0.0)G1=0.
2112 WAVG=(SQRT(WTW)+SQRT(G1))/2.
2113 CONV=RATE*DZ/WAVG
2114 !
2115 ! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
2116 ! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
2117 ! IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
2118 ! SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...
2119 !
2120 RATIO3=QNEWLQ/(QNEW+1.E-10)
2121 ! OLDQ=QTOT
2122 QTOT=QTOT+0.6*QNEW
2123 OLDQ=QTOT
2124 RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-10)
2125 QTOT=QTOT*EXP(-CONV)
2126 !
2127 ! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
2128 ! PARCEL AT THIS LEVEL...
2129 !
2130 DQ=OLDQ-QTOT
2131 QLQOUT=RATIO4*DQ
2132 QICOUT=(1.-RATIO4)*DQ
2133 !
2134 ! ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
2135 ! LATE VERTICAL VELOCITY
2136 !
2137 PPTDRG=0.5*(OLDQ+QTOT-0.2*QNEW)
2138 WTW=WTW+BOTERM-ENTERM-2.*G*DZ*PPTDRG/1.5
2139 !
2140 ! DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
2141 ! DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...
2142 !
2143 QLIQ=RATIO4*QTOT+RATIO3*0.4*QNEW
2144 QICE=(1.-RATIO4)*QTOT+(1.-RATIO3)*0.4*QNEW
2145 QNEWLQ=0.
2146 QNEWIC=0.
2148 END SUBROUTINE CONDLOAD
2150 !-----------------------------------------------------------------------
2151 SUBROUTINE DTFRZNEW(TU,P,THTEU,QVAP,QLIQ,QICE,RATIO2,TTFRZ,TBFRZ, &
2152 QNWFRZ,RL,FRC1,EFFQ,IFLAG,XLV0,XLV1,XLS0,XLS1, &
2153 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2154 !-----------------------------------------------------------------------
2155 IMPLICIT NONE
2156 !-----------------------------------------------------------------------
2157 REAL, INTENT(IN ) :: XLV0,XLV1
2158 REAL, INTENT(IN ) :: P,TTFRZ,TBFRZ,EFFQ,XLS0,XLS1,EP2,ALIQ, &
2159 BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE
2160 REAL, INTENT(INOUT) :: TU,THTEU,QVAP,QLIQ,QICE,RATIO2, &
2161 FRC1,RL,QNWFRZ
2162 INTEGER, INTENT(INOUT) :: IFLAG
2164 REAL :: CCP,RV,C5,QLQFRZ,QNEW,ESLIQ,ESICE,RLC,RLS,PI,ES,RLF,A, &
2165 B,C,DQVAP,DTFRZ,TU1,QVAP1
2166 !-----------------------------------------------------------------------
2167 !
2168 !...ALLOW GLACIATION OF THE UPDRAFT TO OCCUR AS AN APPROXIMATELY LINEAR
2169 ! FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE TTFRZ TO TBFRZ...
2170 !
2172 RV=461.5
2173 C5=1.0723E-3
2174 !
2175 !...ADJUST THE LIQUID WATER CONCENTRATIONS FROM FRESH CONDENSATE AND THA
2176 ! BROUGHT UP FROM LOWER LEVELS TO AN AMOUNT THAT WOULD BE PRESENT IF N
2177 ! LIQUID WATER HAD FROZEN THUS FAR...THIS IS NECESSARY BECAUSE THE
2178 ! EXPRESSION FOR TEMP CHANGE IS MULTIPLIED BY THE FRACTION EQUAL TO TH
2179 ! PARCEL TEMP DECREASE SINCE THE LAST MODEL LEVEL DIVIDED BY THE TOTAL
2180 ! GLACIATION INTERVAL, SO THAT EFFECTIVELY THIS APPROXIMATELY ALLOWS A
2181 ! AMOUNT OF LIQUID WATER TO FREEZE WHICH IS EQUAL TO THIS SAME FRACTIO
2182 ! OF THE LIQUID WATER THAT WAS PRESENT BEFORE THE GLACIATION PROCESS W
2183 ! INITIATED...ALSO, TO ALLOW THETAU TO CONVERT APPROXIMATELY LINEARLY
2184 ! ITS VALUE WITH RESPECT TO ICE, WE NEED TO ALLOW A PORTION OF THE FRE
2185 ! CONDENSATE TO CONTRIBUTE TO THE GLACIATION PROCESS; THE FRACTIONAL
2186 ! AMOUNT THAT APPLIES TO THIS PORTION IS 1/2 OF THE FRACTIONAL AMOUNT
2187 ! FROZEN OF THE "OLD" CONDENSATE BECAUSE THIS FRESH CONDENSATE IS ONLY
2188 ! PRODUCED GRADUALLY OVER THE LAYER...NOTE THAT IN TERMS OF THE DYNAMI
2189 ! OF THE PRECIPITATION PROCESS, IE. PRECIPITATION FALLOUT, THIS FRACTI
2190 ! AMNT OF FRESH CONDENSATE HAS ALREADY BEEN INCLUDED IN THE ICE CATEGO
2191 !
2192 QLQFRZ=QLIQ*EFFQ
2193 QNEW=QNWFRZ*EFFQ*0.5
2194 ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2195 ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
2196 RLC=2.5E6-2369.276*(TU-273.16)
2197 RLS=2833922.-259.532*(TU-273.16)
2198 RLF=RLS-RLC
2199 CCP=1005.7*(1.+0.89*QVAP)
2200 !
2201 ! A = D(ES)/DT IS THAT CALCULATED FROM BUCK`S (1981) EMPIRICAL FORMULAS
2202 ! FOR SATURATION VAPOR PRESSURE...
2203 !
2204 A=(CICE-BICE*DICE)/((TU-DICE)*(TU-DICE))
2205 B=RLS*EP2/P
2206 C=A*B*ESICE/CCP
2207 DQVAP=B*(ESLIQ-ESICE)/(RLS+RLS*C)-RLF*(QLQFRZ+QNEW)/(RLS+RLS/C)
2208 DTFRZ=(RLF*(QLQFRZ+QNEW)+B*(ESLIQ-ESICE))/(CCP+A*B*ESICE)
2209 TU1=TU
2210 QVAP1=QVAP
2211 TU=TU+FRC1*DTFRZ
2212 QVAP=QVAP-FRC1*DQVAP
2213 ES=QVAP*P/(EP2+QVAP)
2214 ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2215 ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
2216 RATIO2=(ESLIQ-ES)/(ESLIQ-ESICE)
2217 !
2218 ! TYPICALLY, RATIO2 IS VERY CLOSE TO (TTFRZ-TU)/(TTFRZ-TBFRZ), USUALLY
2219 ! WITHIN 1% (USING TU BEFORE GALCIATION EFFECTS ARE APPLIED); IF THE
2220 ! INITIAL UPDRAFT TEMP IS BELOW TBFRZ AND RATIO2 IS STILL LESS THAN 1,
2221 ! AN ADJUSTMENT TO FRC1 AND RATIO2 IS INTRODUCED SO THAT GLACIATION
2222 ! EFFECTS ARE NOT UNDERESTIMATED; CONVERSELY, IF RATIO2 IS GREATER THAN
2223 ! FRC1 IS ADJUSTED SO THAT GLACIATION EFFECTS ARE NOT OVERESTIMATED...
2224 !
2225 IF(IFLAG.GT.0.AND.RATIO2.LT.1)THEN
2226 FRC1=FRC1+(1.-RATIO2)
2227 TU=TU1+FRC1*DTFRZ
2228 QVAP=QVAP1-FRC1*DQVAP
2229 RATIO2=1.
2230 IFLAG=1
2231 GOTO 20
2232 ENDIF
2233 IF(RATIO2.GT.1.)THEN
2234 FRC1=FRC1-(RATIO2-1.)
2235 FRC1=AMAX1(0.0,FRC1)
2236 TU=TU1+FRC1*DTFRZ
2237 QVAP=QVAP1-FRC1*DQVAP
2238 RATIO2=1.
2239 IFLAG=1
2240 ENDIF
2241 !
2242 ! CALCULATE A HYBRID VALUE OF THETAU, ASSUMING THAT THE LATENT HEAT OF
2243 ! VAPORIZATION/SUBLIMATION CAN BE ESTIMATED USING THE SAME WEIGHTING
2244 ! FUNCTION AS THAT USED TO CALCULATE SATURATION VAPOR PRESSURE, CALCU-
2245 ! LATE NEW LIQUID WATER AND ICE CONCENTRATIONS...
2246 !
2247 20 RLC=XLV0-XLV1*TU
2248 RLS=XLS0-XLS1*TU
2249 RL=RATIO2*RLS+(1.-RATIO2)*RLC
2250 PI=(1.E5/P)**(0.2854*(1.-0.28*QVAP))
2251 THTEU=TU*PI*EXP(RL*QVAP*C5/TU*(1.+0.81*QVAP))
2252 IF(IFLAG.EQ.1)THEN
2253 QICE=QICE+FRC1*DQVAP+QLIQ
2254 QLIQ=0.
2255 ELSE
2256 QICE=QICE+FRC1*(DQVAP+QLQFRZ)
2257 QLIQ=QLIQ-FRC1*QLQFRZ
2258 ENDIF
2259 QNWFRZ=0.
2261 END SUBROUTINE DTFRZNEW
2263 !-----------------------------------------------------------------------
2264 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2265 ! THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN
2266 ! DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN F
2267 ! HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMA
2268 ! TABLES ED. BY ABRAMOWITZ AND STEGUN, NAT L BUREAU OF STANDARDS APPLI
2269 ! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
2270 ! JACK KAIN
2271 ! 7/6/89
2272 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2273 !***********************************************************************
2274 !***** GAUSSIAN TYPE MIXING PROFILE....******************************
2275 SUBROUTINE PROF5(EQ,EE,UD)
2276 !-----------------------------------------------------------------------
2277 IMPLICIT NONE
2278 !-----------------------------------------------------------------------
2279 REAL, INTENT(IN ) :: EQ
2280 REAL, INTENT(INOUT) :: EE,UD
2281 REAL :: SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
2283 DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676, &
2284 0.9372980,0.33267,0.166666667,0.202765151/
2285 X=(EQ-0.5)/SIGMA
2286 Y=6.*EQ-3.
2287 EY=EXP(Y*Y/(-2))
2288 E45=EXP(-4.5)
2289 T2=1./(1.+P*ABS(Y))
2290 T1=0.500498
2291 C1=A1*T1+A2*T1*T1+A3*T1*T1*T1
2292 C2=A1*T2+A2*T2*T2+A3*T2*T2*T2
2293 IF(Y.GE.0.)THEN
2294 EE=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
2295 UD=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*(0.5+EQ*EQ/2.- &
2296 EQ)
2297 ELSE
2298 EE=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
2299 UD=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*(0.5+EQ* &
2300 EQ/2.-EQ)
2301 ENDIF
2302 EE=EE/FE
2303 UD=UD/FE
2305 END SUBROUTINE PROF5
2307 !-----------------------------------------------------------------------
2308 SUBROUTINE TPMIX(P,THTU,TU,QU,QLIQ,QICE,QNEWLQ,QNEWIC,RATIO2,RL, &
2309 XLV0,XLV1,XLS0,XLS1, &
2310 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2311 !-----------------------------------------------------------------------
2312 IMPLICIT NONE
2313 !-----------------------------------------------------------------------
2314 REAL, INTENT(IN ) :: XLV0,XLV1
2315 REAL, INTENT(IN ) :: P,THTU,RATIO2,RL,XLS0, &
2316 XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,&
2317 CICE,DICE
2318 REAL, INTENT(INOUT) :: QU,QLIQ,QICE,TU,QNEWLQ,QNEWIC
2319 REAL :: ES,QS,PI,THTGS,F0,T1,T0,C5,RV,ESLIQ,ESICE,F1,DT,QNEW, &
2320 DQ, QTOT,DQICE,DQLIQ,RLL,CCP
2321 INTEGER :: ITCNT
2322 !-----------------------------------------------------------------------
2323 !
2324 !...THIS SUBROUTINE ITERATIVELY EXTRACTS WET-BULB TEMPERATURE FROM EQUIV
2325 ! POTENTIAL TEMPERATURE, THEN CHECKS TO SEE IF SUFFICIENT MOISTURE IS
2326 ! AVAILABLE TO ACHIEVE SATURATION...IF NOT, TEMPERATURE IS ADJUSTED
2327 ! ACCORDINGLY, IF SO, THE RESIDUAL LIQUID WATER/ICE CONCENTRATION IS
2328 ! DETERMINED...
2329 !
2330 C5=1.0723E-3
2331 RV=461.5
2332 !
2333 ! ITERATE TO FIND WET BULB TEMPERATURE AS A FUNCTION OF EQUIVALENT POT
2334 ! TEMP AND PRS, ASSUMING SATURATION VAPOR PRESSURE...RATIO2 IS THE DEG
2335 ! OF GLACIATION...
2336 !
2337 IF(RATIO2.LT.1.E-6)THEN
2338 ES=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2339 QS=EP2*ES/(P-ES)
2340 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2341 THTGS=TU*PI*EXP((3374.6525/TU-2.5403)*QS*(1.+0.81*QS))
2342 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
2343 ES=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
2344 QS=EP2*ES/(P-ES)
2345 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2346 THTGS=TU*PI*EXP((3114.834/TU-0.278296)*QS*(1.+0.81*QS))
2347 ELSE
2348 ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2349 ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
2350 ES=(1.-RATIO2)*ESLIQ+RATIO2*ESICE
2351 QS=EP2*ES/(P-ES)
2352 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2353 THTGS=TU*PI*EXP(RL*QS*C5/TU*(1.+0.81*QS))
2354 ENDIF
2355 F0=THTGS-THTU
2356 T1=TU-0.5*F0
2357 T0=TU
2358 ITCNT=0
2359 90 IF(RATIO2.LT.1.E-6)THEN
2360 ES=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
2361 QS=EP2*ES/(P-ES)
2362 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2363 THTGS=T1*PI*EXP((3374.6525/T1-2.5403)*QS*(1.+0.81*QS))
2364 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
2365 ES=AICE*EXP((BICE*T1-CICE)/(T1-DICE))
2366 QS=EP2*ES/(P-ES)
2367 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2368 THTGS=T1*PI*EXP((3114.834/T1-0.278296)*QS*(1.+0.81*QS))
2369 ELSE
2370 ESLIQ=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
2371 ESICE=AICE*EXP((BICE*T1-CICE)/(T1-DICE))
2372 ES=(1.-RATIO2)*ESLIQ+RATIO2*ESICE
2373 QS=EP2*ES/(P-ES)
2374 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2375 THTGS=T1*PI*EXP(RL*QS*C5/T1*(1.+0.81*QS))
2376 ENDIF
2377 F1=THTGS-THTU
2378 IF(ABS(F1).LT.0.01)GOTO 50
2379 ITCNT=ITCNT+1
2380 IF(ITCNT.GT.10)GOTO 50
2381 DT=F1*(T1-T0)/(F1-F0)
2382 T0=T1
2383 F0=F1
2384 T1=T1-DT
2385 GOTO 90
2386 !
2387 ! IF THE PARCEL IS SUPERSATURATED, CALCULATE CONCENTRATION OF FRESH
2388 ! CONDENSATE...
2389 !
2390 50 IF(QS.LE.QU)THEN
2391 QNEW=QU-QS
2392 QU=QS
2393 GOTO 96
2394 ENDIF
2395 !
2396 ! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
2397 ! ADJUSTED...IF LIQUID WATER OR ICE IS PRESENT, IT IS ALLOWED TO EVAPO
2398 ! SUBLIMATE.
2399 !
2400 QNEW=0.
2401 DQ=QS-QU
2402 QTOT=QLIQ+QICE
2403 !
2404 ! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS
2405 ! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MI
2406 ! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURA
2407 ! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPR
2408 ! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE
2409 !
2410 !...NOTE THAT THE LIQ AND ICE MAY BE PRESENT IN PROPORTIONS SLIGHTLY DIF
2411 ! THAN SUGGESTED BY THE VALUE OF RATIO2...CHECK TO MAKE SURE THAT LIQ
2412 ! ICE CONCENTRATIONS ARE NOT REDUCED TO BELOW ZERO WHEN EVAPORATION/
2413 ! SUBLIMATION OCCURS...
2414 !
2415 IF(QTOT.GE.DQ)THEN
2416 DQICE=0.0
2417 DQLIQ=0.0
2418 QLIQ=QLIQ-(1.-RATIO2)*DQ
2419 IF(QLIQ.LT.0.)THEN
2420 DQICE=0.0-QLIQ
2421 QLIQ=0.0
2422 ENDIF
2423 QICE=QICE-RATIO2*DQ+DQICE
2424 IF(QICE.LT.0.)THEN
2425 DQLIQ=0.0-QICE
2426 QICE=0.0
2427 ENDIF
2428 QLIQ=QLIQ+DQLIQ
2429 QU=QS
2430 GOTO 96
2431 ELSE
2432 IF(RATIO2.LT.1.E-6)THEN
2433 RLL=XLV0-XLV1*T1
2434 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
2435 RLL=XLS0-XLS1*T1
2436 ELSE
2437 RLL=RL
2438 ENDIF
2439 CCP=1005.7*(1.+0.89*QU)
2440 IF(QTOT.LT.1.E-10)THEN
2441 !
2442 !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
2443 T1=T1+RLL*(DQ/(1.+DQ))/CCP
2444 GOTO 96
2445 ELSE
2446 !
2447 !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURA
2448 ! THE TEMPERATURE IS GIVEN BY:
2449 T1=T1+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CCP
2450 QU=QU+QTOT
2451 QTOT=0.
2452 ENDIF
2453 QLIQ=0
2454 QICE=0.
2455 ENDIF
2456 96 TU=T1
2457 QNEWLQ=(1.-RATIO2)*QNEW
2458 QNEWIC=RATIO2*QNEW
2459 IF(ITCNT.GT.10)PRINT*,'***** NUMBER OF ITERATIONS IN TPMIX =', &
2460 ITCNT
2462 END SUBROUTINE TPMIX
2463 !-----------------------------------------------------------------------
2464 SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,R1,RL, &
2465 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2466 !-----------------------------------------------------------------------
2467 IMPLICIT NONE
2468 !-----------------------------------------------------------------------
2469 REAL, INTENT(IN ) :: P1,T1,Q1,R1,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,&
2470 BICE,CICE,DICE
2471 REAL, INTENT(INOUT) :: THT1
2472 REAL:: T00,P00,C1,C2,C3,C4,C5,EE,TLOG,TDPT,TSAT,THT,TFPT,TLOGIC, &
2473 TSATLQ,TSATIC
2475 DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834,&
2476 0.278296,1.0723E-3/
2477 !
2478 ! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
2479 !
2481 IF(R1.LT.1.E-6)THEN
2482 EE=Q1*P1/(EP2+Q1)
2483 TLOG=ALOG(EE/ALIQ)
2484 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2485 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
2486 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
2487 THT1=THT*EXP((C1/TSAT-C2)*Q1*(1.+0.81*Q1))
2488 ELSEIF(ABS(R1-1.).LT.1.E-6)THEN
2489 EE=Q1*P1/(EP2+Q1)
2490 TLOG=ALOG(EE/AICE)
2491 TFPT=(CICE-DICE*TLOG)/(BICE-TLOG)
2492 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
2493 TSAT=TFPT-(.182+1.13E-3*(TFPT-T00)-3.58E-4*(T1-T00))*(T1-TFPT)
2494 THT1=THT*EXP((C3/TSAT-C4)*Q1*(1.+0.81*Q1))
2495 ELSE
2496 EE=Q1*P1/(EP2+Q1)
2497 TLOG=ALOG(EE/ALIQ)
2498 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2499 TLOGIC=ALOG(EE/AICE)
2500 TFPT=(CICE-DICE*TLOGIC)/(BICE-TLOGIC)
2501 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
2502 TSATLQ=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
2503 TSATIC=TFPT-(.182+1.13E-3*(TFPT-T00)-3.58E-4*(T1-T00))*(T1-TFPT)
2504 TSAT=R1*TSATIC+(1.-R1)*TSATLQ
2505 THT1=THT*EXP(RL*Q1*C5/TSAT*(1.+0.81*Q1))
2506 ENDIF
2508 END SUBROUTINE ENVIRTHT
2510 !-----------------------------------------------------------------------
2511 !************************* TPDD.FOR ************************************
2512 ! THIS SUBROUTINE ITERATIVELY EXTRACTS TEMPERATURE FROM EQUIVALENT *
2513 ! POTENTIAL TEMP. IT IS DESIGNED FOR USE WITH DOWNDRAFT CALCULATIONS.
2514 ! IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, PARCEL *
2515 ! TEMP, SPECIFIC HUMIDITY, AND LIQUID WATER CONTENT ARE ITERATIVELY *
2516 ! CALCULATED. *
2517 !***********************************************************************
2518 FUNCTION TPDD(P,THTED,TGS,RS,RD,RH,XLV0,XLV1, &
2519 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2520 !-----------------------------------------------------------------------
2521 IMPLICIT NONE
2522 !-----------------------------------------------------------------------
2523 REAL, INTENT(IN ) :: XLV0,XLV1
2524 REAL, INTENT(IN ) :: P,THTED,TGS,RD,RH,EP2,ALIQ,BLIQ, &
2525 CLIQ,DLIQ,AICE,BICE,CICE,DICE
2526 REAL, INTENT(INOUT) :: RS
2527 REAL :: TPDD,ES,PI,THTGS,F0,T1,T0,CCP,F1,DT,RL,DSSDT,T1RH,RSRH
2528 INTEGER :: ITCNT
2529 !-----------------------------------------------------------------------
2530 ES=ALIQ*EXP((BLIQ*TGS-CLIQ)/(TGS-DLIQ))
2531 RS=EP2*ES/(P-ES)
2532 PI=(1.E5/P)**(0.2854*(1.-0.28*RS))
2533 THTGS=TGS*PI*EXP((3374.6525/TGS-2.5403)*RS*(1.+0.81*RS))
2534 F0=THTGS-THTED
2535 T1=TGS-0.5*F0
2536 T0=TGS
2537 CCP=1005.7
2538 !
2539 !...ITERATE TO FIND WET-BULB TEMPERATURE...
2540 !
2541 ITCNT=0
2542 90 ES=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
2543 RS=EP2*ES/(P-ES)
2544 PI=(1.E5/P)**(0.2854*(1.-0.28*RS))
2545 THTGS=T1*PI*EXP((3374.6525/T1-2.5403)*RS*(1.+0.81*RS))
2546 F1=THTGS-THTED
2547 IF(ABS(F1).LT.0.05)GOTO 50
2548 ITCNT=ITCNT+1
2549 IF(ITCNT.GT.10)GOTO 50
2550 DT=F1*(T1-T0)/(F1-F0)
2551 T0=T1
2552 F0=F1
2553 T1=T1-DT
2554 GOTO 90
2555 50 RL=XLV0-XLV1*T1
2556 !
2557 !...IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, ESTIMATE THE
2558 ! TEMPERATURE AND MIXING RATIO WHICH WILL YIELD THE APPROPRIATE VALUE.
2559 !
2560 IF(RH.EQ.1.)GOTO 110
2561 DSSDT=(CLIQ-BLIQ*DLIQ)/((T1-DLIQ)*(T1-DLIQ))
2562 DT=RL*RS*(1.-RH)/(CCP+RL*RH*RS*DSSDT)
2563 T1RH=T1+DT
2564 ES=RH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
2565 RSRH=EP2*ES/(P-ES)
2566 !
2567 !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
2568 !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
2569 !
2570 IF(RSRH.LT.RD)THEN
2571 RSRH=RD
2572 T1RH=T1+(RS-RSRH)*RL/CCP
2573 ENDIF
2574 T1=T1RH
2575 RS=RSRH
2576 110 TPDD=T1
2577 IF(ITCNT.GT.10)PRINT*,'***** NUMBER OF ITERATIONS IN TPDD = ', &
2578 ITCNT
2580 END FUNCTION TPDD
2582 !====================================================================
2583 SUBROUTINE kfinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, &
2584 RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS, &
2585 P_FIRST_SCALAR,restart,allowed_to_read, &
2586 ids, ide, jds, jde, kds, kde, &
2587 ims, ime, jms, jme, kms, kme, &
2588 its, ite, jts, jte, kts, kte )
2589 !--------------------------------------------------------------------
2590 IMPLICIT NONE
2591 !--------------------------------------------------------------------
2592 LOGICAL , INTENT(IN) :: restart, allowed_to_read
2593 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
2594 ims, ime, jms, jme, kms, kme, &
2595 its, ite, jts, jte, kts, kte
2596 INTEGER , INTENT(IN) :: P_QI,P_QS,P_FIRST_SCALAR
2598 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
2599 RTHCUTEN, &
2600 RQVCUTEN, &
2601 RQCCUTEN, &
2602 RQRCUTEN, &
2603 RQICUTEN, &
2604 RQSCUTEN
2606 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG
2608 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA
2610 INTEGER :: i, j, k, itf, jtf, ktf
2612 jtf=min0(jte,jde-1)
2613 ktf=min0(kte,kde-1)
2614 itf=min0(ite,ide-1)
2616 IF(.not.restart)THEN
2617 DO j=jts,jtf
2618 DO k=kts,ktf
2619 DO i=its,itf
2620 RTHCUTEN(i,k,j)=0.
2621 RQVCUTEN(i,k,j)=0.
2622 RQCCUTEN(i,k,j)=0.
2623 RQRCUTEN(i,k,j)=0.
2624 ENDDO
2625 ENDDO
2626 ENDDO
2628 IF (P_QI .ge. P_FIRST_SCALAR) THEN
2629 DO j=jts,jtf
2630 DO k=kts,ktf
2631 DO i=its,itf
2632 RQICUTEN(i,k,j)=0.
2633 ENDDO
2634 ENDDO
2635 ENDDO
2636 ENDIF
2638 IF (P_QS .ge. P_FIRST_SCALAR) THEN
2639 DO j=jts,jtf
2640 DO k=kts,ktf
2641 DO i=its,itf
2642 RQSCUTEN(i,k,j)=0.
2643 ENDDO
2644 ENDDO
2645 ENDDO
2646 ENDIF
2648 DO j=jts,jtf
2649 DO i=its,itf
2650 NCA(i,j)=-100.
2651 ENDDO
2652 ENDDO
2654 DO j=jts,jtf
2655 DO k=kts,ktf
2656 DO i=its,itf
2657 W0AVG(i,k,j)=0.
2658 ENDDO
2659 ENDDO
2660 ENDDO
2662 ENDIF
2664 END SUBROUTINE kfinit
2666 !-------------------------------------------------------
2668 END MODULE module_cu_kf