module_cu_kfeta.F
References to this file elsewhere.
1 MODULE module_cu_kfeta
2
3 USE module_wrf_error
4
5 !--------------------------------------------------------------------
6 ! Lookup table variables:
7 INTEGER, PARAMETER :: KFNT=250,KFNP=220
8 REAL, DIMENSION(KFNT,KFNP),PRIVATE, SAVE :: TTAB,QSTAB
9 REAL, DIMENSION(KFNP),PRIVATE, SAVE :: THE0K
10 REAL, DIMENSION(200),PRIVATE, SAVE :: ALU
11 REAL, PRIVATE, SAVE :: RDPR,RDTHK,PLUTOP
12 ! Note: KF Lookup table is used by subroutines KF_eta_PARA, TPMIX2,
13 ! TPMIX2DD, ENVIRTHT
14 ! End of Lookup table variables:
15
16 CONTAINS
17
18 SUBROUTINE KF_eta_CPS( &
19 ids,ide, jds,jde, kds,kde &
20 ,ims,ime, jms,jme, kms,kme &
21 ,its,ite, jts,jte, kts,kte &
22 ,DT,KTAU,DX &
23 ,rho,RAINCV,NCA &
24 ,U,V,TH,T,W,dz8w,Pcps,pi &
25 ,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
26 ,EP2,SVP1,SVP2,SVP3,SVPT0 &
27 ,STEPCU,CU_ACT_FLAG,warm_rain,CUTOP,CUBOT &
28 ,QV &
29 ! optionals
30 ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
31 ,RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN &
32 ,RQICUTEN,RQSCUTEN &
33 )
34 !
35 !-------------------------------------------------------------
36 IMPLICIT NONE
37 !-------------------------------------------------------------
38 INTEGER, INTENT(IN ) :: &
39 ids,ide, jds,jde, kds,kde, &
40 ims,ime, jms,jme, kms,kme, &
41 its,ite, jts,jte, kts,kte
42
43 INTEGER, INTENT(IN ) :: STEPCU
44 LOGICAL, INTENT(IN ) :: warm_rain
45
46 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1
47 REAL, INTENT(IN ) :: CP,R,G,EP1,EP2
48 REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
49
50 INTEGER, INTENT(IN ) :: KTAU
51
52 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
53 INTENT(IN ) :: &
54 U, &
55 V, &
56 W, &
57 TH, &
58 T, &
59 QV, &
60 dz8w, &
61 Pcps, &
62 rho, &
63 pi
64 !
65 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
66 INTENT(INOUT) :: &
67 W0AVG
68
69 REAL, INTENT(IN ) :: DT, DX
70 !
71 REAL, DIMENSION( ims:ime , jms:jme ), &
72 INTENT(INOUT) :: RAINCV
73
74 REAL, DIMENSION( ims:ime , jms:jme ), &
75 INTENT(INOUT) :: NCA
76
77 REAL, DIMENSION( ims:ime , jms:jme ), &
78 INTENT(OUT) :: CUBOT, &
79 CUTOP
80
81 LOGICAL, DIMENSION( ims:ime , jms:jme ), &
82 INTENT(INOUT) :: CU_ACT_FLAG
83
84 !
85 ! Optional arguments
86 !
87
88 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
89 OPTIONAL, &
90 INTENT(INOUT) :: &
91 RTHCUTEN, &
92 RQVCUTEN, &
93 RQCCUTEN, &
94 RQRCUTEN, &
95 RQICUTEN, &
96 RQSCUTEN
97
98 !
99 ! Flags relating to the optional tendency arrays declared above
100 ! Models that carry the optional tendencies will provdide the
101 ! optional arguments at compile time; these flags all the model
102 ! to determine at run-time whether a particular tracer is in
103 ! use or not.
104 !
105 LOGICAL, OPTIONAL :: &
106 F_QV &
107 ,F_QC &
108 ,F_QR &
109 ,F_QI &
110 ,F_QS
111
112
113 ! LOCAL VARS
114
115 LOGICAL :: flag_qr, flag_qi, flag_qs
116
117 REAL, DIMENSION( kts:kte ) :: &
118 U1D, &
119 V1D, &
120 T1D, &
121 DZ1D, &
122 QV1D, &
123 P1D, &
124 RHO1D, &
125 W0AVG1D
126
127 REAL, DIMENSION( kts:kte ):: &
128 DQDT, &
129 DQIDT, &
130 DQCDT, &
131 DQRDT, &
132 DQSDT, &
133 DTDT
134
135 REAL :: TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp,RTHCUMAX
136
137 INTEGER :: i,j,k,NTST,ICLDCK
138 !
139 DXSQ=DX*DX
140
141 !----------------------
142 NTST=STEPCU
143 TST=float(NTST*2)
144 flag_qr = .FALSE.
145 flag_qi = .FALSE.
146 flag_qs = .FALSE.
147 IF ( PRESENT(F_QR) ) flag_qr = F_QR
148 IF ( PRESENT(F_QI) ) flag_qi = F_QI
149 IF ( PRESENT(F_QS) ) flag_qs = F_QS
150 !
151 DO J = jts,jte
152 DO K=kts,kte
153 DO I= its,ite
154 ! SCR1=-5.0E-4*G*rho(I,K,J)*(w(I,K,J)+w(I,K+1,J))
155 ! TV=T(I,K,J)*(1.+EP1*QV(I,K,J))
156 ! RHOE=Pcps(I,K,J)/(R*TV)
157 ! W0=-101.9368*SCR1/RHOE
158 W0=0.5*(w(I,K,J)+w(I,K+1,J))
159 W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST
160 ENDDO
161 ENDDO
162 ENDDO
163 !
164 !...CHECK FOR CONVECTIVE INITIATION EVERY 5 MINUTES (OR NTST/2)...
165 !
166 !----------------------
167 ICLDCK=MOD(KTAU,NTST)
168 IF(ICLDCK.EQ.0 .or. KTAU .eq. 1) then
169 !
170 DO J = jts,jte
171 DO I= its,ite
172 CU_ACT_FLAG(i,j) = .true.
173 ENDDO
174 ENDDO
175
176 DO J = jts,jte
177 DO I=its,ite
178
179 IF ( NINT(NCA(I,J)) .gt. 0 ) then
180 CU_ACT_FLAG(i,j) = .false.
181 ELSE
182
183 DO k=kts,kte
184 DQDT(k)=0.
185 DQIDT(k)=0.
186 DQCDT(k)=0.
187 DQRDT(k)=0.
188 DQSDT(k)=0.
189 DTDT(k)=0.
190 ENDDO
191 RAINCV(I,J)=0.
192 CUTOP(I,J)=KTS
193 CUBOT(I,J)=KTE+1
194 !
195 ! assign vars from 3D to 1D
196
197 DO K=kts,kte
198 U1D(K) =U(I,K,J)
199 V1D(K) =V(I,K,J)
200 T1D(K) =T(I,K,J)
201 RHO1D(K) =rho(I,K,J)
202 QV1D(K)=QV(I,K,J)
203 P1D(K) =Pcps(I,K,J)
204 W0AVG1D(K) =W0AVG(I,K,J)
205 DZ1D(k)=dz8w(I,K,J)
206 ENDDO
207 CALL KF_eta_PARA(I, J, &
208 U1D,V1D,T1D,QV1D,P1D,DZ1D, &
209 W0AVG1D,DT,DX,DXSQ,RHO1D, &
210 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
211 EP2,SVP1,SVP2,SVP3,SVPT0, &
212 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
213 RAINCV,NCA,NTST, &
214 flag_QI,flag_QS,warm_rain, &
215 CUTOP,CUBOT, &
216 ids,ide, jds,jde, kds,kde, &
217 ims,ime, jms,jme, kms,kme, &
218 its,ite, jts,jte, kts,kte)
219 IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
220 DO K=kts,kte
221 RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
222 ! RTHCUMAX=max(abs(RTHCUTEN(I,K,J)),RTHCUMAX)
223 RQVCUTEN(I,K,J)=DQDT(K)
224 ENDDO
225 ENDIF
226
227 IF(PRESENT(rqrcuten).AND.PRESENT(rqccuten)) THEN
228 IF( F_QR )THEN
229 DO K=kts,kte
230 RQRCUTEN(I,K,J)=DQRDT(K)
231 RQCCUTEN(I,K,J)=DQCDT(K)
232 ENDDO
233 ELSE
234 ! This is the case for Eta microphysics without 3d rain field
235 DO K=kts,kte
236 RQRCUTEN(I,K,J)=0.
237 RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
238 ENDDO
239 ENDIF
240 ENDIF
241
242 !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
243
244
245 IF(PRESENT( rqicuten )) THEN
246 IF ( F_QI ) THEN
247 DO K=kts,kte
248 RQICUTEN(I,K,J)=DQIDT(K)
249 ENDDO
250 ENDIF
251 ENDIF
252
253 IF(PRESENT( rqscuten )) THEN
254 IF ( F_QS ) THEN
255 DO K=kts,kte
256 RQSCUTEN(I,K,J)=DQSDT(K)
257 ENDDO
258 ENDIF
259 ENDIF
260 !
261 ENDIF
262 ENDDO
263 ENDDO
264 ENDIF
265 !
266 END SUBROUTINE KF_eta_CPS
267 ! ****************************************************************************
268 !-----------------------------------------------------------
269 SUBROUTINE KF_eta_PARA (I, J, &
270 U0,V0,T0,QV0,P0,DZQ,W0AVG1D, &
271 DT,DX,DXSQ,rhoe, &
272 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
273 EP2,SVP1,SVP2,SVP3,SVPT0, &
274 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
275 RAINCV,NCA,NTST, &
276 F_QI,F_QS,warm_rain, &
277 CUTOP,CUBOT, &
278 ids,ide, jds,jde, kds,kde, &
279 ims,ime, jms,jme, kms,kme, &
280 its,ite, jts,jte, kts,kte)
281 !-----------------------------------------------------------
282 !***** The KF scheme that is currently used in experimental runs of EMCs
283 !***** Eta model....jsk 8/00
284 !
285 IMPLICIT NONE
286 !-----------------------------------------------------------
287 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
288 ims,ime, jms,jme, kms,kme, &
289 its,ite, jts,jte, kts,kte, &
290 I,J,NTST
291 ! ,P_QI,P_QS,P_FIRST_SCALAR
292
293 LOGICAL, INTENT(IN ) :: F_QI, F_QS
294
295 LOGICAL, INTENT(IN ) :: warm_rain
296 !
297 REAL, DIMENSION( kts:kte ), &
298 INTENT(IN ) :: U0, &
299 V0, &
300 T0, &
301 QV0, &
302 P0, &
303 rhoe, &
304 DZQ, &
305 W0AVG1D
306 !
307 REAL, INTENT(IN ) :: DT,DX,DXSQ
308 !
309
310 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
311 REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
312
313 !
314 REAL, DIMENSION( kts:kte ), INTENT(INOUT) :: &
315 DQDT, &
316 DQIDT, &
317 DQCDT, &
318 DQRDT, &
319 DQSDT, &
320 DTDT
321
322 REAL, DIMENSION( ims:ime , jms:jme ), &
323 INTENT(INOUT) :: NCA
324
325 REAL, DIMENSION( ims:ime , jms:jme ), &
326 INTENT(INOUT) :: RAINCV
327
328 REAL, DIMENSION( ims:ime , jms:jme ), &
329 INTENT(OUT) :: CUBOT, &
330 CUTOP
331 !
332 !...DEFINE LOCAL VARIABLES...
333 !
334 REAL, DIMENSION( kts:kte ) :: &
335 Q0,Z0,TV0,TU,TVU,QU,TZ,TVD, &
336 QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD, &
337 UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2, &
338 UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE, &
339 THTAU,THETEU,THTAD,THETED,QLIQ,QICE, &
340 QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC, &
341 DETLQ2,DETIC2,RATIO,RATIO2
342
343
344 REAL, DIMENSION( kts:kte ) :: &
345 DOMGDP,EXN,TVQU,DP,RH,EQFRC,WSPD, &
346 QDT,FXM,THTAG,THPA,THFXOUT, &
347 THFXIN,QPA,QFXOUT,QFXIN,QLPA,QLFXIN, &
348 QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA, &
349 QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT, &
350 QL0,QLG,QI0,QIG,QR0,QRG,QS0,QSG
351
352
353 REAL, DIMENSION( kts:kte+1 ) :: OMG
354 REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
355 REAL, DIMENSION( kts:kte ) :: &
356 CLDHGT,QSD,DILFRC,DDILFRC,TKE,TGU,QGU,THTEEG
357
358 ! LOCAL VARS
359
360 REAL :: P00,T00,RLF,RHIC,RHBC,PIE, &
361 TTFRZ,TBFRZ,C5,RATE
362 REAL :: GDRY,ROCP,ALIQ,BLIQ, &
363 CLIQ,DLIQ
364 REAL :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX, &
365 ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL, &
366 CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR, &
367 ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
368 TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD, &
369 UPNEW,ABE,WKLCL,TTEMP,FRC1, &
370 QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
371 DZZ,UDLBE,REI,EE2,UD2,TTMP,F1,F2, &
372 THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1, &
373 UD1,DPTT,QNEWLQ,DUMFDP,EE,TSAT, &
374 THTA,VCONV,TIMEC,SHSIGN,VWS,PEF, &
375 CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN, &
376 DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1, &
377 DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF, &
378 UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF, &
379 DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
380 AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1, &
381 DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF, &
382 TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR, &
383 UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2, &
384 RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
385 DDFRC,TDC,DEFRC,RHBAR,DMFFRC,DPMIN,DILBE
386 REAL :: ASTRT,TP,VALUE,AINTRP,TKEMAX,QFRZ,&
387 QSS,PPTMLT,DTMELT,RHH,EVAC,BINC
388 !
389 INTEGER :: INDLU,NU,NUCHM,NNN,KLFS
390 REAL :: CHMIN,PM15,CHMAX,DTRH,RAD,DPPP
391 REAL :: TVDIFF,DTTOT,ABSOMG,ABSOMGTC,FRDP
392
393 INTEGER :: KX,K,KL
394 !
395 INTEGER :: NCHECK
396 INTEGER, DIMENSION (kts:kte) :: KCHECK
397
398 INTEGER :: ISTOP,ML,L5,KMIX,LOW, &
399 LC,MXLAYR,LLFC,NLAYRS,NK, &
400 KPBL,KLCL,LCL,LET,IFLAG, &
401 NK1,LTOP,NJ,LTOP1, &
402 LTOPM1,LVF,KSTART,KMIN,LFS, &
403 ND,NIC,LDB,LDT,ND1,NDK, &
404 NM,LMAX,NCOUNT,NOITR, &
405 NSTEP,NTC,NCHM,ISHALL,NSHALL
406 LOGICAL :: IPRNT
407 CHARACTER*1024 message
408 !
409 DATA P00,T00/1.E5,273.16/
410 DATA RLF/3.339E5/
411 DATA RHIC,RHBC/1.,0.90/
412 DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
413 DATA RATE/0.03/
414 !-----------------------------------------------------------
415 IPRNT=.FALSE.
416 GDRY=-G/CP
417 ROCP=R/CP
418 NSHALL = 0
419 KL=kte
420 KX=kte
421 !
422 ! ALIQ = 613.3
423 ! BLIQ = 17.502
424 ! CLIQ = 4780.8
425 ! DLIQ = 32.19
426 ALIQ = SVP1*1000.
427 BLIQ = SVP2
428 CLIQ = SVP2*SVPT0
429 DLIQ = SVP3
430 !
431 !
432 !****************************************************************************
433 ! ! PPT FB MODS
434 !...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER ! PPT FB MODS
435 !...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL) ! PPT FB MODS
436 !...FIELD. "FBFRC" IS THE FRACTION OF AVAILABLE ! PPT FB MODS
437 !...PRECIPITATION TO BE FED BACK (0.0 - 1.0)... ! PPT FB MODS
438 FBFRC=0.0 ! PPT FB MODS
439 !...mods to allow shallow convection...
440 NCHM = 0
441 ISHALL = 0
442 DPMIN = 5.E3
443 !...
444 P300=P0(1)-30000.
445 !
446 !...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
447 !...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
448 !...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
449 !
450 !...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
451 !...FROM BOTTOM-UP IN THE KF SCHEME...
452 !
453 ML=0
454 !SUE tmprpsb=1./PSB(I,J)
455 !SUE CELL=PTOP*tmprpsb
456 !
457 DO K=1,KX
458 !
459 !...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
460 !
461 ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
462 QES(K)=0.622*ES/(P0(K)-ES)
463 Q0(K)=AMIN1(QES(K),QV0(K))
464 Q0(K)=AMAX1(0.000001,Q0(K))
465 QL0(K)=0.
466 QI0(K)=0.
467 QR0(K)=0.
468 QS0(K)=0.
469 RH(K) = Q0(K)/QES(K)
470 DILFRC(K) = 1.
471 TV0(K)=T0(K)*(1.+0.608*Q0(K))
472 ! RHOE(K)=P0(K)/(R*TV0(K))
473 ! DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
474 DP(K)=rhoe(k)*g*DZQ(k)
475 ! IF Turbulent Kinetic Energy (TKE) is available from turbulent mixing scheme
476 ! use it for shallow convection...For now, assume it is not available....
477 ! TKE(K) = Q2(I,J,NK)
478 TKE(K) = 0.
479 CLDHGT(K) = 0.
480 ! IF(P0(K).GE.500E2)L5=K
481 IF(P0(K).GE.0.5*P0(1))L5=K
482 IF(P0(K).GE.P300)LLFC=K
483 IF(T0(K).GT.T00)ML=K
484 ENDDO
485 !
486 !...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
487 Z0(1)=.5*DZQ(1)
488 !cdir novector
489 DO K=2,KL
490 Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
491 DZA(K-1)=Z0(K)-Z0(K-1)
492 ENDDO
493 DZA(KL)=0.
494 !
495 !
496 ! To save time, specify a pressure interval to move up in sequential
497 ! check of different ~50 mb deep groups of adjacent model layers in
498 ! the process of identifying updraft source layer (USL). Note that
499 ! this search is terminated as soon as a buoyant parcel is found and
500 ! this parcel can produce a cloud greater than specifed minimum depth
501 ! (CHMIN)...For now, set interval at 15 mb...
502 !
503 NCHECK = 1
504 KCHECK(NCHECK)=1
505 PM15 = P0(1)-15.E2
506 DO K=2,LLFC
507 IF(P0(K).LT.PM15)THEN
508 NCHECK = NCHECK+1
509 KCHECK(NCHECK) = K
510 PM15 = PM15-15.E2
511 ENDIF
512 ENDDO
513 !
514 NU=0
515 NUCHM=0
516 usl: DO
517 NU = NU+1
518 IF(NU.GT.NCHECK)THEN
519 IF(ISHALL.EQ.1)THEN
520 CHMAX = 0.
521 NCHM = 0
522 DO NK = 1,NCHECK
523 NNN=KCHECK(NK)
524 IF(CLDHGT(NNN).GT.CHMAX)THEN
525 NCHM = NNN
526 NUCHM = NK
527 CHMAX = CLDHGT(NNN)
528 ENDIF
529 ENDDO
530 NU = NUCHM-1
531 FBFRC=1.
532 CYCLE usl
533 ELSE
534 RETURN
535 ENDIF
536 ENDIF
537 KMIX = KCHECK(NU)
538 LOW=KMIX
539 !...
540 LC = LOW
541 !
542 !...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
543 !...UNSTABLE AIR AT LEAST 50 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
544 !...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
545 !...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 50 mb..
546 !
547 NLAYRS=0
548 DPTHMX=0.
549 NK=LC-1
550 IF ( NK+1 .LT. KTS ) THEN
551 WRITE(message,*)'WOULD GO OFF BOTTOM: KF_ETA_PARA I,J,NK',I,J,NK
552 CALL wrf_message (TRIM(message))
553 ELSE
554 DO
555 NK=NK+1
556 IF ( NK .GT. KTE ) THEN
557 WRITE(message,*)'WOULD GO OFF TOP: KF_ETA_PARA I,J,DPTHMX,DPMIN',I,J,DPTHMX,DPMIN
558 CALL wrf_message (TRIM(message))
559 EXIT
560 ENDIF
561 DPTHMX=DPTHMX+DP(NK)
562 NLAYRS=NLAYRS+1
563 IF(DPTHMX.GT.DPMIN)THEN
564 EXIT
565 ENDIF
566 END DO
567 ENDIF
568 IF(DPTHMX.LT.DPMIN)THEN
569 RETURN
570 ENDIF
571 KPBL=LC+NLAYRS-1
572 !
573 !...********************************************************
574 !...for computational simplicity without much loss in accuracy,
575 !...mix temperature instead of theta for evaluating convective
576 !...initiation (triggering) potential...
577 ! THMIX=0.
578 TMIX=0.
579 QMIX=0.
580 ZMIX=0.
581 PMIX=0.
582 !
583 !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
584 !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
585 !...LAYERS...
586 !
587 !cdir novector
588 DO NK=LC,KPBL
589 TMIX=TMIX+DP(NK)*T0(NK)
590 QMIX=QMIX+DP(NK)*Q0(NK)
591 ZMIX=ZMIX+DP(NK)*Z0(NK)
592 PMIX=PMIX+DP(NK)*P0(NK)
593 ENDDO
594 ! THMIX=THMIX/DPTHMX
595 TMIX=TMIX/DPTHMX
596 QMIX=QMIX/DPTHMX
597 ZMIX=ZMIX/DPTHMX
598 PMIX=PMIX/DPTHMX
599 EMIX=QMIX*PMIX/(0.622+QMIX)
600 !
601 !...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL...
602 !
603 ! TLOG=ALOG(EMIX/ALIQ)
604 ! ...calculate dewpoint using lookup table...
605 !
606 astrt=1.e-3
607 ainc=0.075
608 a1=emix/aliq
609 tp=(a1-astrt)/ainc
610 indlu=int(tp)+1
611 value=(indlu-1)*ainc+astrt
612 aintrp=(a1-value)/ainc
613 tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
614 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
615 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
616 TLCL=AMIN1(TLCL,TMIX)
617 TVLCL=TLCL*(1.+0.608*QMIX)
618 ZLCL = ZMIX+(TLCL-TMIX)/GDRY
619 NK = LC-1
620 DO
621 NK = NK+1
622 KLCL=NK
623 IF(ZLCL.LE.Z0(NK) .or. NK.GT.KL)THEN
624 EXIT
625 ENDIF
626 ENDDO
627 IF(NK.GT.KL)THEN
628 RETURN
629 ENDIF
630 K=KLCL-1
631 DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
632 !
633 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
634 !
635 TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
636 QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
637 TVEN=TENV*(1.+0.608*QENV)
638 !
639 !...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
640 !...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0 IS AN
641 !...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
642 !...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
643 !...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
644 !...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
645 !...SUCCESS AT GRID LENGTHS NEAR 25 km. FOR DIFFERENT GRID-LENGTHS,
646 !...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
647 !...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
648 !
649 IF(ZLCL.LT.2.E3)THEN
650 WKLCL=0.02*ZLCL/2.E3
651 ELSE
652 WKLCL=0.02
653 ENDIF
654 WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)*DX/25.E3-WKLCL
655 IF(WKL.LT.0.0001)THEN
656 DTLCL=0.
657 ELSE
658 DTLCL=4.64*WKL**0.33
659 ENDIF
660 !
661 !...for ETA model, give parcel an extra temperature perturbation based
662 !...the threshold RH for condensation (U00)...
663 !
664 !...for now, just assume U00=0.75...
665 !...!!!!!! for MM5, SET DTRH = 0. !!!!!!!!
666 ! U00 = 0.75
667 ! IF(U00.lt.1.)THEN
668 ! QSLCL=QES(K)+(QES(KLCL)-QES(K))*DLP
669 ! RHLCL = QENV/QSLCL
670 ! DQSSDT = QMIX*(CLIQ-BLIQ*DLIQ)/((TLCL-DLIQ)*(TLCL-DLIQ))
671 ! IF(RHLCL.ge.0.75 .and. RHLCL.le.0.95)then
672 ! DTRH = 0.25*(RHLCL-0.75)*QMIX/DQSSDT
673 ! ELSEIF(RHLCL.GT.0.95)THEN
674 ! DTRH = (1./RHLCL-1.)*QMIX/DQSSDT
675 ! ELSE
676 DTRH = 0.
677 ! ENDIF
678 ! ENDIF
679 ! IF(ISHALL.EQ.1)IPRNT=.TRUE.
680 ! IPRNT=.TRUE.
681 ! IF(TLCL+DTLCL.GT.TENV)GOTO 45
682 !
683 trigger: IF(TLCL+DTLCL+DTRH.LT.TENV)THEN
684 !
685 ! Parcel not buoyant, CYCLE back to start of trigger and evaluate next potential USL...
686 !
687 CYCLE usl
688 !
689 ELSE ! Parcel is buoyant, determine updraft
690 !
691 !...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
692 !...EQUIVALENT POTENTIAL TEMPERATURE
693 !...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
694 !
695 CALL ENVIRTHT(PMIX,TMIX,QMIX,THETEU(K),ALIQ,BLIQ,CLIQ,DLIQ)
696 !
697 !...modify calculation of initial parcel vertical velocity...jsk 11/26/97
698 !
699 DTTOT = DTLCL+DTRH
700 IF(DTTOT.GT.1.E-4)THEN
701 GDT=2.*G*DTTOT*500./TVEN
702 WLCL=1.+0.5*SQRT(GDT)
703 WLCL = AMIN1(WLCL,3.)
704 ELSE
705 WLCL=1.
706 ENDIF
707 PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
708 WTW=WLCL*WLCL
709 !
710 TVLCL=TLCL*(1.+0.608*QMIX)
711 RHOLCL=PLCL/(R*TVLCL)
712 !
713 LCL=KLCL
714 LET=LCL
715 ! make RAD a function of background vertical velocity...
716 IF(WKL.LT.0.)THEN
717 RAD = 1000.
718 ELSEIF(WKL.GT.0.1)THEN
719 RAD = 2000.
720 ELSE
721 RAD = 1000.+1000*WKL/0.1
722 ENDIF
723 !
724 !*******************************************************************
725 ! *
726 ! COMPUTE UPDRAFT PROPERTIES *
727 ! *
728 !*******************************************************************
729 !
730 !
731 !...
732 !...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
733 !
734 WU(K)=WLCL
735 AU0=0.01*DXSQ
736 UMF(K)=RHOLCL*AU0
737 VMFLCL=UMF(K)
738 UPOLD=VMFLCL
739 UPNEW=UPOLD
740 !
741 !...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
742 !...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE
743 !...BUOYANT ENERGY, TRPPT IS THE TOTAL RATE OF PRECIPITATION
744 !...PRODUCTION...
745 !
746 RATIO2(K)=0.
747 UER(K)=0.
748 ABE=0.
749 TRPPT=0.
750 TU(K)=TLCL
751 TVU(K)=TVLCL
752 QU(K)=QMIX
753 EQFRC(K)=1.
754 QLIQ(K)=0.
755 QICE(K)=0.
756 QLQOUT(K)=0.
757 QICOUT(K)=0.
758 DETLQ(K)=0.
759 DETIC(K)=0.
760 PPTLIQ(K)=0.
761 PPTICE(K)=0.
762 IFLAG=0
763 !
764 !...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
765 !...PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
766 !...FREEZING IS SPECIFIED TO BEGIN. WITHIN THE GLACIATION
767 !...INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
768 !...PREVIOUS MODEL LEVEL...
769 !
770 TTEMP=TTFRZ
771 !
772 !...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
773 !...MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
774 !...MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
775 !
776 !
777 EE1=1.
778 UD1=0.
779 REI = 0.
780 DILBE = 0.
781 updraft: DO NK=K,KL-1
782 NK1=NK+1
783 RATIO2(NK1)=RATIO2(NK)
784 FRC1=0.
785 TU(NK1)=T0(NK1)
786 THETEU(NK1)=THETEU(NK)
787 QU(NK1)=QU(NK)
788 QLIQ(NK1)=QLIQ(NK)
789 QICE(NK1)=QICE(NK)
790 call tpmix2(p0(nk1),theteu(nk1),tu(nk1),qu(nk1),qliq(nk1), &
791 qice(nk1),qnewlq,qnewic,XLV1,XLV0)
792 !
793 !
794 !...CHECK TO SEE IF UPDRAFT TEMP IS ABOVE THE TEMPERATURE AT WHICH
795 !...GLACIATION IS ASSUMED TO INITIATE; IF IT IS, CALCULATE THE
796 !...FRACTION OF REMAINING LIQUID WATER TO FREEZE...TTFRZ IS THE
797 !...TEMP AT WHICH FREEZING BEGINS, TBFRZ THE TEMP BELOW WHICH ALL
798 !...LIQUID WATER IS FROZEN AT EACH LEVEL...
799 !
800 IF(TU(NK1).LE.TTFRZ)THEN
801 IF(TU(NK1).GT.TBFRZ)THEN
802 IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
803 FRC1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
804 ELSE
805 FRC1=1.
806 IFLAG=1
807 ENDIF
808 TTEMP=TU(NK1)
809 !
810 ! DETERMINE THE EFFECTS OF LIQUID WATER FREEZING WHEN TEMPERATURE
811 !...IS BELOW TTFRZ...
812 !
813 QFRZ = (QLIQ(NK1)+QNEWLQ)*FRC1
814 QNEWIC=QNEWIC+QNEWLQ*FRC1
815 QNEWLQ=QNEWLQ-QNEWLQ*FRC1
816 QICE(NK1) = QICE(NK1)+QLIQ(NK1)*FRC1
817 QLIQ(NK1) = QLIQ(NK1)-QLIQ(NK1)*FRC1
818 CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QFRZ, &
819 QICE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
820 ENDIF
821 TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
822 !
823 ! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
824 !
825 IF(NK.EQ.K)THEN
826 BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
827 BOTERM=2.*(Z0(NK1)-ZLCL)*G*BE/1.5
828 DZZ=Z0(NK1)-ZLCL
829 ELSE
830 BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
831 BOTERM=2.*DZA(NK)*G*BE/1.5
832 DZZ=DZA(NK)
833 ENDIF
834 ENTERM=2.*REI*WTW/UPOLD
835
836 CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM, &
837 RATE,QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1),G)
838 !
839 !...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
840 !...IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
841 !
842 IF(WTW.LT.1.E-3)THEN
843 EXIT
844 ELSE
845 WU(NK1)=SQRT(WTW)
846 ENDIF
847 !...Calculate value of THETA-E in environment to entrain into updraft...
848 !
849 CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
850 !
851 !...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
852 !
853 REI=VMFLCL*DP(NK1)*0.03/RAD
854 TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
855 IF(NK.EQ.K)THEN
856 DILBE=((TVLCL+TVQU(NK1))/(TVEN+TV0(NK1))-1.)*DZZ
857 ELSE
858 DILBE=((TVQU(NK)+TVQU(NK1))/(TV0(NK)+TV0(NK1))-1.)*DZZ
859 ENDIF
860 IF(DILBE.GT.0.)ABE=ABE+DILBE*G
861 !
862 !...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, MINIMAL
863 !...ENTRAINMENT (0.5*REI) IS IMPOSED...
864 !
865 IF(TVQU(NK1).LE.TV0(NK1))THEN ! Entrain/Detrain IF BLOCK
866 EE2=0.5
867 UD2=1.
868 EQFRC(NK1)=0.
869 ELSE
870 LET=NK1
871 TTMP=TVQU(NK1)
872 !
873 !...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL AIR...
874 !
875 F1=0.95
876 F2=1.-F1
877 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
878 QTMP=F1*Q0(NK1)+F2*QU(NK1)
879 TMPLIQ=F2*QLIQ(NK1)
880 TMPICE=F2*QICE(NK1)
881 call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
882 qnewlq,qnewic,XLV1,XLV0)
883 TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
884 IF(TU95.GT.TV0(NK1))THEN
885 EE2=1.
886 UD2=0.
887 EQFRC(NK1)=1.0
888 ELSE
889 F1=0.10
890 F2=1.-F1
891 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
892 QTMP=F1*Q0(NK1)+F2*QU(NK1)
893 TMPLIQ=F2*QLIQ(NK1)
894 TMPICE=F2*QICE(NK1)
895 call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
896 qnewlq,qnewic,XLV1,XLV0)
897 TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
898 TVDIFF = ABS(TU10-TVQU(NK1))
899 IF(TVDIFF.LT.1.e-3)THEN
900 EE2=1.
901 UD2=0.
902 EQFRC(NK1)=1.0
903 ELSE
904 EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
905 EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
906 EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
907 IF(EQFRC(NK1).EQ.1)THEN
908 EE2=1.
909 UD2=0.
910 ELSEIF(EQFRC(NK1).EQ.0.)THEN
911 EE2=0.
912 UD2=1.
913 ELSE
914 !
915 !...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
916 ! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
917 !
918 CALL PROF5(EQFRC(NK1),EE2,UD2)
919 ENDIF
920 ENDIF
921 ENDIF
922 ENDIF ! End of Entrain/Detrain IF BLOCK
923 !
924 !
925 !...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE FRACTIONAL
926 ! VALUES IN THE LAYER...
927 !
928 EE2 = AMAX1(EE2,0.5)
929 UD2 = 1.5*UD2
930 UER(NK1)=0.5*REI*(EE1+EE2)
931 UDR(NK1)=0.5*REI*(UD1+UD2)
932 !
933 !...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
934 ! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATIONS...
935 !
936 IF(UMF(NK)-UDR(NK1).LT.10.)THEN
937 !
938 !...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL UPD MASS
939 ! FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE PREVIOUS MODEL LVL..
940 ! First, correct ABE calculation if needed...
941 !
942 IF(DILBE.GT.0.)THEN
943 ABE=ABE-DILBE*G
944 ENDIF
945 LET=NK
946 ! WRITE(98,1015)P0(NK1)/100.
947 EXIT
948 ELSE
949 EE1=EE2
950 UD1=UD2
951 UPOLD=UMF(NK)-UDR(NK1)
952 UPNEW=UPOLD+UER(NK1)
953 UMF(NK1)=UPNEW
954 DILFRC(NK1) = UPNEW/UPOLD
955 !
956 !...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND
957 !...ICE IN THE DETRAINING UPDRAFT MASS...
958 !
959 DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
960 DETIC(NK1)=QICE(NK1)*UDR(NK1)
961 QDT(NK1)=QU(NK1)
962 QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
963 THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
964 QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
965 QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
966 !
967 !...PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF
968 !...LIQUID PRECIP AT A GIVEN MODEL LVL, PPTICE THE SAME FOR ICE,
969 !...TRPPT IS THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE
970 !...CURRENT MODEL LEVEL...
971 !
972 PPTLIQ(NK1)=QLQOUT(NK1)*UMF(NK)
973 PPTICE(NK1)=QICOUT(NK1)*UMF(NK)
974 !
975 TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
976 IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
977 ENDIF
978 !
979 END DO updraft
980 !
981 !...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIU
982 ! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
983 ! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE BE
984 ! THE LET AND CLOUD TOP...
985 !
986 !...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL VELOC
987 ! FIRST BECOMES NEGATIVE...
988 !
989 LTOP=NK
990 CLDHGT(LC)=Z0(LTOP)-ZLCL
991 !
992 !...Instead of using the same minimum cloud height (for deep convection)
993 !...everywhere, try specifying minimum cloud depth as a function of TLCL...
994 !
995 !
996 !
997 IF(TLCL.GT.293.)THEN
998 CHMIN = 4.E3
999 ELSEIF(TLCL.LE.293. .and. TLCL.GE.273)THEN
1000 CHMIN = 2.E3 + 100.*(TLCL-273.)
1001 ELSEIF(TLCL.LT.273.)THEN
1002 CHMIN = 2.E3
1003 ENDIF
1004
1005 !
1006 !...If cloud top height is less than the specified minimum for deep
1007 !...convection, save value to consider this level as source for
1008 !...shallow convection, go back up to check next level...
1009 !
1010 !...Try specifying minimum cloud depth as a function of TLCL...
1011 !
1012 !
1013 !...DO NOT ALLOW ANY CLOUD FROM THIS LAYER IF:
1014 !
1015 !... 1.) if there is no CAPE, or
1016 !... 2.) cloud top is at model level just above LCL, or
1017 !... 3.) cloud top is within updraft source layer, or
1018 !... 4.) cloud-top detrainment layer begins within
1019 !... updraft source layer.
1020 !
1021 IF(LTOP.LE.KLCL .or. LTOP.LE.KPBL .or. LET+1.LE.KPBL)THEN ! No Convection Allowed
1022 CLDHGT(LC)=0.
1023 DO NK=K,LTOP
1024 UMF(NK)=0.
1025 UDR(NK)=0.
1026 UER(NK)=0.
1027 DETLQ(NK)=0.
1028 DETIC(NK)=0.
1029 PPTLIQ(NK)=0.
1030 PPTICE(NK)=0.
1031 ENDDO
1032 !
1033 ELSEIF(CLDHGT(LC).GT.CHMIN .and. ABE.GT.1)THEN ! Deep Convection allowed
1034 ISHALL=0
1035 EXIT usl
1036 ELSE
1037 !
1038 !...TO DISALLOW SHALLOW CONVECTION, COMMENT OUT NEXT LINE !!!!!!!!
1039 ISHALL = 1
1040 IF(NU.EQ.NUCHM)THEN
1041 EXIT usl ! Shallow Convection from this layer
1042 ELSE
1043 ! Remember this layer (by virtue of non-zero CLDHGT) as potential shallow-cloud layer
1044 DO NK=K,LTOP
1045 UMF(NK)=0.
1046 UDR(NK)=0.
1047 UER(NK)=0.
1048 DETLQ(NK)=0.
1049 DETIC(NK)=0.
1050 PPTLIQ(NK)=0.
1051 PPTICE(NK)=0.
1052 ENDDO
1053 ENDIF
1054 ENDIF
1055 ENDIF trigger
1056 END DO usl
1057 IF(ISHALL.EQ.1)THEN
1058 KSTART=MAX0(KPBL,KLCL)
1059 LET=KSTART
1060 endif
1061 !
1062 !...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS FL
1063 ! THIS LEVEL...
1064 !
1065 IF(LET.EQ.LTOP)THEN
1066 UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
1067 DETLQ(LTOP)=QLIQ(LTOP)*UDR(LTOP)*UPNEW/UPOLD
1068 DETIC(LTOP)=QICE(LTOP)*UDR(LTOP)*UPNEW/UPOLD
1069 UER(LTOP)=0.
1070 UMF(LTOP)=0.
1071 ELSE
1072 !
1073 ! BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
1074 !
1075 DPTT=0.
1076 DO NJ=LET+1,LTOP
1077 DPTT=DPTT+DP(NJ)
1078 ENDDO
1079 DUMFDP=UMF(LET)/DPTT
1080 !
1081 !...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
1082 ! RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
1083 !
1084 DO NK=LET+1,LTOP
1085 !
1086 !...entrainment is allowed at every level except for LTOP, so disallow
1087 !...entrainment at LTOP and adjust entrainment rates between LET and LTOP
1088 !...so the the dilution factor due to entyrianment is not changed but
1089 !...the actual entrainment rate will change due due forced total
1090 !...detrainment in this layer...
1091 !
1092 IF(NK.EQ.LTOP)THEN
1093 UDR(NK) = UMF(NK-1)
1094 UER(NK) = 0.
1095 DETLQ(NK) = UDR(NK)*QLIQ(NK)*DILFRC(NK)
1096 DETIC(NK) = UDR(NK)*QICE(NK)*DILFRC(NK)
1097 ELSE
1098 UMF(NK)=UMF(NK-1)-DP(NK)*DUMFDP
1099 UER(NK)=UMF(NK)*(1.-1./DILFRC(NK))
1100 UDR(NK)=UMF(NK-1)-UMF(NK)+UER(NK)
1101 DETLQ(NK)=UDR(NK)*QLIQ(NK)*DILFRC(NK)
1102 DETIC(NK)=UDR(NK)*QICE(NK)*DILFRC(NK)
1103 ENDIF
1104 IF(NK.GE.LET+2)THEN
1105 TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
1106 PPTLIQ(NK)=UMF(NK-1)*QLQOUT(NK)
1107 PPTICE(NK)=UMF(NK-1)*QICOUT(NK)
1108 TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
1109 ENDIF
1110 ENDDO
1111 ENDIF
1112 !
1113 ! Initialize some arrays below cloud base and above cloud top...
1114 !
1115 DO NK=1,K
1116 IF(NK.GE.LC)THEN
1117 IF(NK.EQ.LC)THEN
1118 UMF(NK)=VMFLCL*DP(NK)/DPTHMX
1119 UER(NK)=VMFLCL*DP(NK)/DPTHMX
1120 ELSEIF(NK.LE.KPBL)THEN
1121 UER(NK)=VMFLCL*DP(NK)/DPTHMX
1122 UMF(NK)=UMF(NK-1)+UER(NK)
1123 ELSE
1124 UMF(NK)=VMFLCL
1125 UER(NK)=0.
1126 ENDIF
1127 TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
1128 QU(NK)=QMIX
1129 WU(NK)=WLCL
1130 ELSE
1131 TU(NK)=0.
1132 QU(NK)=0.
1133 UMF(NK)=0.
1134 WU(NK)=0.
1135 UER(NK)=0.
1136 ENDIF
1137 UDR(NK)=0.
1138 QDT(NK)=0.
1139 QLIQ(NK)=0.
1140 QICE(NK)=0.
1141 QLQOUT(NK)=0.
1142 QICOUT(NK)=0.
1143 PPTLIQ(NK)=0.
1144 PPTICE(NK)=0.
1145 DETLQ(NK)=0.
1146 DETIC(NK)=0.
1147 RATIO2(NK)=0.
1148 CALL ENVIRTHT(P0(NK),T0(NK),Q0(NK),THETEE(NK),ALIQ,BLIQ,CLIQ,DLIQ)
1149 EQFRC(NK)=1.0
1150 ENDDO
1151 !
1152 LTOP1=LTOP+1
1153 LTOPM1=LTOP-1
1154 !
1155 !...DEFINE VARIABLES ABOVE CLOUD TOP...
1156 !
1157 DO NK=LTOP1,KX
1158 UMF(NK)=0.
1159 UDR(NK)=0.
1160 UER(NK)=0.
1161 QDT(NK)=0.
1162 QLIQ(NK)=0.
1163 QICE(NK)=0.
1164 QLQOUT(NK)=0.
1165 QICOUT(NK)=0.
1166 DETLQ(NK)=0.
1167 DETIC(NK)=0.
1168 PPTLIQ(NK)=0.
1169 PPTICE(NK)=0.
1170 IF(NK.GT.LTOP1)THEN
1171 TU(NK)=0.
1172 QU(NK)=0.
1173 WU(NK)=0.
1174 ENDIF
1175 THTA0(NK)=0.
1176 THTAU(NK)=0.
1177 EMS(NK)=0.
1178 EMSD(NK)=0.
1179 TG(NK)=T0(NK)
1180 QG(NK)=Q0(NK)
1181 QLG(NK)=0.
1182 QIG(NK)=0.
1183 QRG(NK)=0.
1184 QSG(NK)=0.
1185 OMG(NK)=0.
1186 ENDDO
1187 OMG(KX+1)=0.
1188 DO NK=1,LTOP
1189 EMS(NK)=DP(NK)*DXSQ/G
1190 EMSD(NK)=1./EMS(NK)
1191 !
1192 !...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION SCH
1193 !
1194 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
1195 THTAU(NK)=TU(NK)*EXN(NK)
1196 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
1197 THTA0(NK)=T0(NK)*EXN(NK)
1198 DDILFRC(NK) = 1./DILFRC(NK)
1199 OMG(NK)=0.
1200 ENDDO
1201 ! IF (XTIME.LT.10.)THEN
1202 ! WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
1203 ! * TMIX-T00,PMIX,QMIX,ABE
1204 ! WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
1205 ! * WLCL,CLDHGT
1206 ! ENDIF
1207 !
1208 !...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
1209 !...AND MIDTROPOSPHERE IS USED.
1210 !
1211 WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
1212 WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
1213 WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
1214 VCONV=.5*(WSPD(KLCL)+WSPD(L5))
1215 !...for ETA model, DX is a function of location...
1216 ! TIMEC=DX(I,J)/VCONV
1217 TIMEC=DX/VCONV
1218 TADVEC=TIMEC
1219 TIMEC=AMAX1(1800.,TIMEC)
1220 TIMEC=AMIN1(3600.,TIMEC)
1221 IF(ISHALL.EQ.1)TIMEC=2400.
1222 NIC=NINT(TIMEC/DT)
1223 TIMEC=FLOAT(NIC)*DT
1224 !
1225 !...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
1226 !
1227 IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
1228 SHSIGN=1.
1229 ELSE
1230 SHSIGN=-1.
1231 ENDIF
1232 VWS=(U0(LTOP)-U0(KLCL))*(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL))* &
1233 (V0(LTOP)-V0(KLCL))
1234 VWS=1.E3*SHSIGN*SQRT(VWS)/(Z0(LTOP)-Z0(LCL))
1235 PEF=1.591+VWS*(-.639+VWS*(9.53E-2-VWS*4.96E-3))
1236 PEF=AMAX1(PEF,.2)
1237 PEF=AMIN1(PEF,.9)
1238 !
1239 !...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
1240 !
1241 CBH=(ZLCL-Z0(1))*3.281E-3
1242 IF(CBH.LT.3.)THEN
1243 RCBH=.02
1244 ELSE
1245 RCBH=.96729352+CBH*(-.70034167+CBH*(.162179896+CBH*(- &
1246 1.2569798E-2+CBH*(4.2772E-4-CBH*5.44E-6))))
1247 ENDIF
1248 IF(CBH.GT.25)RCBH=2.4
1249 PEFCBH=1./(1.+RCBH)
1250 PEFCBH=AMIN1(PEFCBH,.9)
1251 !
1252 !... MEAN PEF. IS USED TO COMPUTE RAINFALL.
1253 !
1254 PEFF=.5*(PEF+PEFCBH)
1255 PEFF2 = PEFF ! JSK MODS
1256 IF(IPRNT)THEN
1257 WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1258 ! call flush(98)
1259 endif
1260 ! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1261 !*****************************************************************
1262 ! *
1263 ! COMPUTE DOWNDRAFT PROPERTIES *
1264 ! *
1265 !*****************************************************************
1266 !
1267 !
1268 TDER=0.
1269 devap:IF(ISHALL.EQ.1)THEN
1270 LFS = 1
1271 ELSE
1272 !
1273 !...start downdraft about 150 mb above cloud base...
1274 !
1275 ! KSTART=MAX0(KPBL,KLCL)
1276 ! KSTART=KPBL ! Changed 7/23/99
1277 KSTART=KPBL+1 ! Changed 7/23/99
1278 KLFS = LET-1
1279 DO NK = KSTART+1,KL
1280 DPPP = P0(KSTART)-P0(NK)
1281 ! IF(DPPP.GT.200.E2)THEN
1282 IF(DPPP.GT.150.E2)THEN
1283 KLFS = NK
1284 EXIT
1285 ENDIF
1286 ENDDO
1287 KLFS = MIN0(KLFS,LET-1)
1288 LFS = KLFS
1289 !
1290 !...if LFS is not at least 50 mb above cloud base (implying that the
1291 !...level of equil temp, LET, is just above cloud base) do not allow a
1292 !...downdraft...
1293 !
1294 IF((P0(KSTART)-P0(LFS)).GT.50.E2)THEN
1295 THETED(LFS) = THETEE(LFS)
1296 QD(LFS) = Q0(LFS)
1297 !
1298 !...call tpmix2dd to find wet-bulb temp, qv...
1299 !
1300 call tpmix2dd(p0(lfs),theted(lfs),tz(lfs),qss,i,j)
1301 THTAD(LFS)=TZ(LFS)*(P00/P0(LFS))**(0.2854*(1.-0.28*QSS))
1302 !
1303 !...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX...
1304 !
1305 TVD(LFS)=TZ(LFS)*(1.+0.608*QSS)
1306 RDD=P0(LFS)/(R*TVD(LFS))
1307 A1=(1.-PEFF)*AU0
1308 DMF(LFS)=-A1*RDD
1309 DER(LFS)=DMF(LFS)
1310 DDR(LFS)=0.
1311 RHBAR = RH(LFS)*DP(LFS)
1312 DPTT = DP(LFS)
1313 DO ND = LFS-1,KSTART,-1
1314 ND1 = ND+1
1315 DER(ND)=DER(LFS)*EMS(ND)/EMS(LFS)
1316 DDR(ND)=0.
1317 DMF(ND)=DMF(ND1)+DER(ND)
1318 THETED(ND)=(THETED(ND1)*DMF(ND1)+THETEE(ND)*DER(ND))/DMF(ND)
1319 QD(ND)=(QD(ND1)*DMF(ND1)+Q0(ND)*DER(ND))/DMF(ND)
1320 DPTT = DPTT+DP(ND)
1321 RHBAR = RHBAR+RH(ND)*DP(ND)
1322 ENDDO
1323 RHBAR = RHBAR/DPTT
1324 DMFFRC = 2.*(1.-RHBAR)
1325 DPDD = 0.
1326 !...Calculate melting effect
1327 !... first, compute total frozen precipitation generated...
1328 !
1329 pptmlt = 0.
1330 DO NK = KLCL,LTOP
1331 PPTMLT = PPTMLT+PPTICE(NK)
1332 ENDDO
1333 if(lc.lt.ml)then
1334 !...For now, calculate melting effect as if DMF = -UMF at KLCL, i.e., as
1335 !...if DMFFRC=1. Otherwise, for small DMFFRC, DTMELT gets too large!
1336 !...12/14/98 jsk...
1337 DTMELT = RLF*PPTMLT/(CP*UMF(KLCL))
1338 else
1339 DTMELT = 0.
1340 endif
1341 LDT = MIN0(LFS-1,KSTART-1)
1342 !
1343 call tpmix2dd(p0(kstart),theted(kstart),tz(kstart),qss,i,j)
1344 !
1345 tz(kstart) = tz(kstart)-dtmelt
1346 ES=ALIQ*EXP((BLIQ*TZ(KSTART)-CLIQ)/(TZ(KSTART)-DLIQ))
1347 QSS=0.622*ES/(P0(KSTART)-ES)
1348 THETED(KSTART)=TZ(KSTART)*(1.E5/P0(KSTART))**(0.2854*(1.-0.28*QSS))* &
1349 EXP((3374.6525/TZ(KSTART)-2.5403)*QSS*(1.+0.81*QSS))
1350 !....
1351 LDT = MIN0(LFS-1,KSTART-1)
1352 DO ND = LDT,1,-1
1353 DPDD = DPDD+DP(ND)
1354 THETED(ND) = THETED(KSTART)
1355 QD(ND) = QD(KSTART)
1356 !
1357 !...call tpmix2dd to find wet bulb temp, saturation mixing ratio...
1358 !
1359 call tpmix2dd(p0(nd),theted(nd),tz(nd),qss,i,j)
1360 qsd(nd) = qss
1361 !
1362 !...specify RH decrease of 20%/km in downdraft...
1363 !
1364 RHH = 1.-0.2/1000.*(Z0(KSTART)-Z0(ND))
1365 !
1366 !...adjust downdraft TEMP, Q to specified RH:
1367 !
1368 IF(RHH.LT.1.)THEN
1369 DSSDT=(CLIQ-BLIQ*DLIQ)/((TZ(ND)-DLIQ)*(TZ(ND)-DLIQ))
1370 RL=XLV0-XLV1*TZ(ND)
1371 DTMP=RL*QSS*(1.-RHH)/(CP+RL*RHH*QSS*DSSDT)
1372 T1RH=TZ(ND)+DTMP
1373 ES=RHH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
1374 QSRH=0.622*ES/(P0(ND)-ES)
1375 !
1376 !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
1377 !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
1378 !
1379 IF(QSRH.LT.QD(ND))THEN
1380 QSRH=QD(ND)
1381 T1RH=TZ(ND)+(QSS-QSRH)*RL/CP
1382 ENDIF
1383 TZ(ND)=T1RH
1384 QSS=QSRH
1385 QSD(ND) = QSS
1386 ENDIF
1387 TVD(nd) = tz(nd)*(1.+0.608*qsd(nd))
1388 IF(TVD(ND).GT.TV0(ND).OR.ND.EQ.1)THEN
1389 LDB=ND
1390 EXIT
1391 ENDIF
1392 ENDDO
1393 IF((P0(LDB)-P0(LFS)) .gt. 50.E2)THEN ! minimum Downdraft depth!
1394 DO ND=LDT,LDB,-1
1395 ND1 = ND+1
1396 DDR(ND) = -DMF(KSTART)*DP(ND)/DPDD
1397 DER(ND) = 0.
1398 DMF(ND) = DMF(ND1)+DDR(ND)
1399 TDER=TDER+(QSD(nd)-QD(ND))*DDR(ND)
1400 QD(ND)=QSD(nd)
1401 THTAD(ND)=TZ(ND)*(P00/P0(ND))**(0.2854*(1.-0.28*QD(ND)))
1402 ENDDO
1403 ENDIF
1404 ENDIF
1405 ENDIF devap
1406 !
1407 !...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
1408 !...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
1409 !
1410 d_mf: IF(TDER.LT.1.)THEN
1411 ! WRITE(98,3004)I,J
1412 !3004 FORMAT(' ','No Downdraft!; I=',I3,2X,'J=',I3,'ISHALL =',I2)
1413 PPTFLX=TRPPT
1414 CPR=TRPPT
1415 TDER=0.
1416 CNDTNF=0.
1417 UPDINC=1.
1418 LDB=LFS
1419 DO NDK=1,LTOP
1420 DMF(NDK)=0.
1421 DER(NDK)=0.
1422 DDR(NDK)=0.
1423 THTAD(NDK)=0.
1424 WD(NDK)=0.
1425 TZ(NDK)=0.
1426 QD(NDK)=0.
1427 ENDDO
1428 AINCM2=100.
1429 ELSE
1430 DDINC = -DMFFRC*UMF(KLCL)/DMF(KSTART)
1431 UPDINC=1.
1432 IF(TDER*DDINC.GT.TRPPT)THEN
1433 DDINC = TRPPT/TDER
1434 ENDIF
1435 TDER = TDER*DDINC
1436 DO NK=LDB,LFS
1437 DMF(NK)=DMF(NK)*DDINC
1438 DER(NK)=DER(NK)*DDINC
1439 DDR(NK)=DDR(NK)*DDINC
1440 ENDDO
1441 CPR=TRPPT
1442 PPTFLX = TRPPT-TDER
1443 PEFF=PPTFLX/TRPPT
1444 IF(IPRNT)THEN
1445 write(98,*)'PRECIP EFFICIENCY =',PEFF
1446 ! call flush(98)
1447 ENDIF
1448 !
1449 !
1450 !...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
1451 ! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
1452 ! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
1453 !
1454 ! DO NK=LC,LFS
1455 ! UMF(NK)=UMF(NK)*UPDINC
1456 ! UDR(NK)=UDR(NK)*UPDINC
1457 ! UER(NK)=UER(NK)*UPDINC
1458 ! PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
1459 ! PPTICE(NK)=PPTICE(NK)*UPDINC
1460 ! DETLQ(NK)=DETLQ(NK)*UPDINC
1461 ! DETIC(NK)=DETIC(NK)*UPDINC
1462 ! ENDDO
1463 !
1464 !...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
1465 !...DOWNDRAFT...
1466 !
1467 IF(LDB.GT.1)THEN
1468 DO NK=1,LDB-1
1469 DMF(NK)=0.
1470 DER(NK)=0.
1471 DDR(NK)=0.
1472 WD(NK)=0.
1473 TZ(NK)=0.
1474 QD(NK)=0.
1475 THTAD(NK)=0.
1476 ENDDO
1477 ENDIF
1478 DO NK=LFS+1,KX
1479 DMF(NK)=0.
1480 DER(NK)=0.
1481 DDR(NK)=0.
1482 WD(NK)=0.
1483 TZ(NK)=0.
1484 QD(NK)=0.
1485 THTAD(NK)=0.
1486 ENDDO
1487 DO NK=LDT+1,LFS-1
1488 TZ(NK)=0.
1489 QD(NK)=0.
1490 THTAD(NK)=0.
1491 ENDDO
1492 ENDIF d_mf
1493 !
1494 !...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE INFL
1495 ! INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN IS AVAILAB
1496 ! IN THAT LAYER INITIALLY...
1497 !
1498 AINCMX=1000.
1499 LMAX=MAX0(KLCL,LFS)
1500 DO NK=LC,LMAX
1501 IF((UER(NK)-DER(NK)).GT.1.e-3)THEN
1502 AINCM1=EMS(NK)/((UER(NK)-DER(NK))*TIMEC)
1503 AINCMX=AMIN1(AINCMX,AINCM1)
1504 ENDIF
1505 ENDDO
1506 AINC=1.
1507 IF(AINCMX.LT.AINC)AINC=AINCMX
1508 !
1509 !...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRAFT AND DOWNDRAFT...THEY WILL
1510 !...BE ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE STABILIZATION
1511 !...CLOSURE...
1512 !
1513 TDER2=TDER
1514 PPTFL2=PPTFLX
1515 DO NK=1,LTOP
1516 DETLQ2(NK)=DETLQ(NK)
1517 DETIC2(NK)=DETIC(NK)
1518 UDR2(NK)=UDR(NK)
1519 UER2(NK)=UER(NK)
1520 DDR2(NK)=DDR(NK)
1521 DER2(NK)=DER(NK)
1522 UMF2(NK)=UMF(NK)
1523 DMF2(NK)=DMF(NK)
1524 ENDDO
1525 FABE=1.
1526 STAB=0.95
1527 NOITR=0
1528 ISTOP=0
1529 !
1530 IF(ISHALL.EQ.1)THEN ! First for shallow convection
1531 !
1532 ! No iteration for shallow convection; if turbulent kinetic energy (TKE) is available
1533 ! from a turbulence parameterization, scale cloud-base updraft mass flux as a function
1534 ! of TKE, but for now, just specify shallow-cloud mass flux using TKEMAX = 5...
1535 !
1536 !...find the maximum TKE value between LC and KLCL...
1537 ! TKEMAX = 0.
1538 TKEMAX = 5.
1539 ! DO 173 K = LC,KLCL
1540 ! NK = KX-K+1
1541 ! TKEMAX = AMAX1(TKEMAX,Q2(I,J,NK))
1542 ! 173 CONTINUE
1543 ! TKEMAX = AMIN1(TKEMAX,10.)
1544 ! TKEMAX = AMAX1(TKEMAX,5.)
1545 !c TKEMAX = 10.
1546 !c...3_24_99...DPMIN was changed for shallow convection so that it is the
1547 !c... the same as for deep convection (5.E3). Since this doubles
1548 !c... (roughly) the value of DPTHMX, add a factor of 0.5 to calcu-
1549 !c... lation of EVAC...
1550 !c EVAC = TKEMAX*0.1
1551 EVAC = 0.5*TKEMAX*0.1
1552 ! AINC = 0.1*DPTHMX*DXIJ*DXIJ/(VMFLCL*G*TIMEC)
1553 ! AINC = EVAC*DPTHMX*DX(I,J)*DX(I,J)/(VMFLCL*G*TIMEC)
1554 AINC = EVAC*DPTHMX*DXSQ/(VMFLCL*G*TIMEC)
1555 TDER=TDER2*AINC
1556 PPTFLX=PPTFL2*AINC
1557 DO NK=1,LTOP
1558 UMF(NK)=UMF2(NK)*AINC
1559 DMF(NK)=DMF2(NK)*AINC
1560 DETLQ(NK)=DETLQ2(NK)*AINC
1561 DETIC(NK)=DETIC2(NK)*AINC
1562 UDR(NK)=UDR2(NK)*AINC
1563 UER(NK)=UER2(NK)*AINC
1564 DER(NK)=DER2(NK)*AINC
1565 DDR(NK)=DDR2(NK)*AINC
1566 ENDDO
1567 ENDIF ! Otherwise for deep convection
1568 ! use iterative procedure to find mass fluxes...
1569 iter: DO NCOUNT=1,10
1570 !
1571 !*****************************************************************
1572 ! *
1573 ! COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE *
1574 ! *
1575 !*****************************************************************
1576 !
1577 !...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
1578 !...SATISFY MASS CONTINUITY...
1579 !
1580 DTT=TIMEC
1581 DO NK=1,LTOP
1582 DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
1583 IF(NK.GT.1)THEN
1584 OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
1585 ABSOMG = ABS(OMG(NK))
1586 ABSOMGTC = ABSOMG*TIMEC
1587 FRDP = 0.75*DP(NK-1)
1588 IF(ABSOMGTC.GT.FRDP)THEN
1589 DTT1 = FRDP/ABSOMG
1590 DTT=AMIN1(DTT,DTT1)
1591 ENDIF
1592 ENDIF
1593 ENDDO
1594 DO NK=1,LTOP
1595 THPA(NK)=THTA0(NK)
1596 QPA(NK)=Q0(NK)
1597 NSTEP=NINT(TIMEC/DTT+1)
1598 DTIME=TIMEC/FLOAT(NSTEP)
1599 FXM(NK)=OMG(NK)*DXSQ/G
1600 ENDDO
1601 !
1602 !...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
1603 !
1604 DO NTC=1,NSTEP
1605 !
1606 !...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED
1607 !...SIGN OF OMEGA...
1608 !
1609 DO NK=1,LTOP
1610 THFXIN(NK)=0.
1611 THFXOUT(NK)=0.
1612 QFXIN(NK)=0.
1613 QFXOUT(NK)=0.
1614 ENDDO
1615 DO NK=2,LTOP
1616 IF(OMG(NK).LE.0.)THEN
1617 THFXIN(NK)=-FXM(NK)*THPA(NK-1)
1618 QFXIN(NK)=-FXM(NK)*QPA(NK-1)
1619 THFXOUT(NK-1)=THFXOUT(NK-1)+THFXIN(NK)
1620 QFXOUT(NK-1)=QFXOUT(NK-1)+QFXIN(NK)
1621 ELSE
1622 THFXOUT(NK)=FXM(NK)*THPA(NK)
1623 QFXOUT(NK)=FXM(NK)*QPA(NK)
1624 THFXIN(NK-1)=THFXIN(NK-1)+THFXOUT(NK)
1625 QFXIN(NK-1)=QFXIN(NK-1)+QFXOUT(NK)
1626 ENDIF
1627 ENDDO
1628 !
1629 !...UPDATE THE THETA AND QV VALUES AT EACH LEVEL...
1630 !
1631 DO NK=1,LTOP
1632 THPA(NK)=THPA(NK)+(THFXIN(NK)+UDR(NK)*THTAU(NK)+DDR(NK)* &
1633 THTAD(NK)-THFXOUT(NK)-(UER(NK)-DER(NK))*THTA0(NK))* &
1634 DTIME*EMSD(NK)
1635 QPA(NK)=QPA(NK)+(QFXIN(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK)- &
1636 QFXOUT(NK)-(UER(NK)-DER(NK))*Q0(NK))*DTIME*EMSD(NK)
1637 ENDDO
1638 ENDDO
1639 DO NK=1,LTOP
1640 THTAG(NK)=THPA(NK)
1641 QG(NK)=QPA(NK)
1642 ENDDO
1643 !
1644 !...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO, BORRO
1645 !...MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO...
1646 !
1647 DO NK=1,LTOP
1648 IF(QG(NK).LT.0.)THEN
1649 IF(NK.EQ.1)THEN ! JSK MODS
1650 ! PRINT *,' PROBLEM WITH KF SCHEME: ' ! JSK MODS
1651 ! PRINT *,'QG = 0 AT THE SURFACE!!!!!!!' ! JSK MODS
1652 CALL wrf_error_fatal ( 'QG, QG(NK).LT.0') ! JSK MODS
1653 ENDIF ! JSK MODS
1654 NK1=NK+1
1655 IF(NK.EQ.LTOP)THEN
1656 NK1=KLCL
1657 ENDIF
1658 TMA=QG(NK1)*EMS(NK1)
1659 TMB=QG(NK-1)*EMS(NK-1)
1660 TMM=(QG(NK)-1.E-9)*EMS(NK )
1661 BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
1662 ACOEFF=BCOEFF*TMA/TMB
1663 TMB=TMB*(1.-BCOEFF)
1664 TMA=TMA*(1.-ACOEFF)
1665 IF(NK.EQ.LTOP)THEN
1666 QVDIFF=(QG(NK1)-TMA*EMSD(NK1))*100./QG(NK1)
1667 ! IF(ABS(QVDIFF).GT.1.)THEN
1668 ! PRINT *,'!!!WARNING!!! CLOUD BASE WATER VAPOR CHANGES BY ', &
1669 ! QVDIFF, &
1670 ! '% WHEN MOISTURE IS BORROWED TO PREVENT NEGATIVE ', &
1671 ! 'VALUES IN KAIN-FRITSCH'
1672 ! ENDIF
1673 ENDIF
1674 QG(NK)=1.E-9
1675 QG(NK1)=TMA*EMSD(NK1)
1676 QG(NK-1)=TMB*EMSD(NK-1)
1677 ENDIF
1678 ENDDO
1679 TOPOMG=(UDR(LTOP)-UER(LTOP))*DP(LTOP)*EMSD(LTOP)
1680 IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
1681 ! WRITE(99,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME; &
1682 ! TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1683 ! TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1684 ISTOP=1
1685 IPRNT=.TRUE.
1686 EXIT iter
1687 ENDIF
1688 !
1689 !...CONVERT THETA TO T...
1690 !
1691 DO NK=1,LTOP
1692 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QG(NK)))
1693 TG(NK)=THTAG(NK)/EXN(NK)
1694 TVG(NK)=TG(NK)*(1.+0.608*QG(NK))
1695 ENDDO
1696 IF(ISHALL.EQ.1)THEN
1697 EXIT iter
1698 ENDIF
1699 !
1700 !*******************************************************************
1701 ! *
1702 ! COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY. *
1703 ! *
1704 !*******************************************************************
1705 !
1706 !...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
1707 !
1708 ! THMIX=0.
1709 TMIX=0.
1710 QMIX=0.
1711 !
1712 !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
1713 !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
1714 !...LAYERS...
1715 !
1716 DO NK=LC,KPBL
1717 TMIX=TMIX+DP(NK)*TG(NK)
1718 QMIX=QMIX+DP(NK)*QG(NK)
1719 ENDDO
1720 TMIX=TMIX/DPTHMX
1721 QMIX=QMIX/DPTHMX
1722 ES=ALIQ*EXP((TMIX*BLIQ-CLIQ)/(TMIX-DLIQ))
1723 QSS=0.622*ES/(PMIX-ES)
1724 !
1725 !...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
1726 !
1727 IF(QMIX.GT.QSS)THEN
1728 RL=XLV0-XLV1*TMIX
1729 CPM=CP*(1.+0.887*QMIX)
1730 DSSDT=QSS*(CLIQ-BLIQ*DLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
1731 DQ=(QMIX-QSS)/(1.+RL*DSSDT/CPM)
1732 TMIX=TMIX+RL/CP*DQ
1733 QMIX=QMIX-DQ
1734 TLCL=TMIX
1735 ELSE
1736 QMIX=AMAX1(QMIX,0.)
1737 EMIX=QMIX*PMIX/(0.622+QMIX)
1738 astrt=1.e-3
1739 binc=0.075
1740 a1=emix/aliq
1741 tp=(a1-astrt)/binc
1742 indlu=int(tp)+1
1743 value=(indlu-1)*binc+astrt
1744 aintrp=(a1-value)/binc
1745 tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
1746 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
1747 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
1748 TLCL=AMIN1(TLCL,TMIX)
1749 ENDIF
1750 TVLCL=TLCL*(1.+0.608*QMIX)
1751 ZLCL = ZMIX+(TLCL-TMIX)/GDRY
1752 DO NK = LC,KL
1753 KLCL=NK
1754 IF(ZLCL.LE.Z0(NK))THEN
1755 EXIT
1756 ENDIF
1757 ENDDO
1758 K=KLCL-1
1759 DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
1760 !
1761 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
1762 !
1763 TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
1764 QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
1765 TVEN=TENV*(1.+0.608*QENV)
1766 PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
1767 THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
1768 EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
1769 !
1770 !...COMPUTE ADJUSTED ABE(ABEG).
1771 !
1772 ABEG=0.
1773 DO NK=K,LTOPM1
1774 NK1=NK+1
1775 THETEU(NK1) = THETEU(NK)
1776 !
1777 call tpmix2dd(p0(nk1),theteu(nk1),tgu(nk1),qgu(nk1),i,j)
1778 !
1779 TVQU(NK1)=TGU(NK1)*(1.+0.608*QGU(NK1)-QLIQ(NK1)-QICE(NK1))
1780 IF(NK.EQ.K)THEN
1781 DZZ=Z0(KLCL)-ZLCL
1782 DILBE=((TVLCL+TVQU(NK1))/(TVEN+TVG(NK1))-1.)*DZZ
1783 ELSE
1784 DZZ=DZA(NK)
1785 DILBE=((TVQU(NK)+TVQU(NK1))/(TVG(NK)+TVG(NK1))-1.)*DZZ
1786 ENDIF
1787 IF(DILBE.GT.0.)ABEG=ABEG+DILBE*G
1788 !
1789 !...DILUTE BY ENTRAINMENT BY THE RATE AS ORIGINAL UPDRAFT...
1790 !
1791 CALL ENVIRTHT(P0(NK1),TG(NK1),QG(NK1),THTEEG(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
1792 THETEU(NK1)=THETEU(NK1)*DDILFRC(NK1)+THTEEG(NK1)*(1.-DDILFRC(NK1))
1793 ENDDO
1794 !
1795 !...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
1796 !...THE PERIOD TIMEC...
1797 !
1798 IF(NOITR.EQ.1)THEN
1799 ! write(98,*)' '
1800 ! write(98,*)'TAU, I, J, =',NTSD,I,J
1801 ! WRITE(98,1060)FABE
1802 ! GOTO 265
1803 EXIT iter
1804 ENDIF
1805 DABE=AMAX1(ABE-ABEG,0.1*ABE)
1806 FABE=ABEG/ABE
1807 IF(FABE.GT.1. .and. ISHALL.EQ.0)THEN
1808 ! WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS
1809 ! *GRID POINT; NO CONVECTION ALLOWED!'
1810 RETURN
1811 ENDIF
1812 IF(NCOUNT.NE.1)THEN
1813 IF(ABS(AINC-AINCOLD).LT.0.0001)THEN
1814 NOITR=1
1815 AINC=AINCOLD
1816 CYCLE iter
1817 ENDIF
1818 DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
1819 IF(DFDA.GT.0.)THEN
1820 NOITR=1
1821 AINC=AINCOLD
1822 CYCLE iter
1823 ENDIF
1824 ENDIF
1825 AINCOLD=AINC
1826 FABEOLD=FABE
1827 IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
1828 ! write(98,*)' '
1829 ! write(98,*)'TAU, I, J, =',NTSD,I,J
1830 ! WRITE(98,1055)FABE
1831 ! GOTO 265
1832 EXIT
1833 ENDIF
1834 IF((FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB) .or. NCOUNT.EQ.10)THEN
1835 EXIT iter
1836 ELSE
1837 IF(NCOUNT.GT.10)THEN
1838 ! write(98,*)' '
1839 ! write(98,*)'TAU, I, J, =',NTSD,I,J
1840 ! WRITE(98,1060)FABE
1841 ! GOTO 265
1842 EXIT
1843 ENDIF
1844 !
1845 !...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE CONVECTI
1846 !...MASS FLUX BY THE FACTOR AINC:
1847 !
1848 IF(FABE.EQ.0.)THEN
1849 AINC=AINC*0.5
1850 ELSE
1851 IF(DABE.LT.1.e-4)THEN
1852 NOITR=1
1853 AINC=AINCOLD
1854 CYCLE iter
1855 ELSE
1856 AINC=AINC*STAB*ABE/DABE
1857 ENDIF
1858 ENDIF
1859 ! AINC=AMIN1(AINCMX,AINC)
1860 AINC=AMIN1(AINCMX,AINC)
1861 !...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION ! JSK MODS
1862 !...WILL BE MINIMAL SO JUST IGNORE IT... ! JSK MODS
1863 IF(AINC.LT.0.05)then
1864 RETURN ! JSK MODS
1865 ENDIF
1866 ! AINC=AMAX1(AINC,0.05) ! JSK MODS
1867 TDER=TDER2*AINC
1868 PPTFLX=PPTFL2*AINC
1869 ! IF (XTIME.LT.10.)THEN
1870 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,
1871 ! * FABEOLD,AINCOLD
1872 ! ENDIF
1873 DO NK=1,LTOP
1874 UMF(NK)=UMF2(NK)*AINC
1875 DMF(NK)=DMF2(NK)*AINC
1876 DETLQ(NK)=DETLQ2(NK)*AINC
1877 DETIC(NK)=DETIC2(NK)*AINC
1878 UDR(NK)=UDR2(NK)*AINC
1879 UER(NK)=UER2(NK)*AINC
1880 DER(NK)=DER2(NK)*AINC
1881 DDR(NK)=DDR2(NK)*AINC
1882 ENDDO
1883 !
1884 !...GO BACK UP FOR ANOTHER ITERATION...
1885 !
1886 ENDIF
1887 ENDDO iter
1888 !
1889 !...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
1890 !
1891 !...FRC2 IS THE FRACTION OF TOTAL CONDENSATE ! PPT FB MODS
1892 !...GENERATED THAT GOES INTO PRECIPITIATION ! PPT FB MODS
1893 !
1894 ! Redistribute hydormeteors according to the final mass-flux values:
1895 !
1896 IF(CPR.GT.0.)THEN
1897 FRC2=PPTFLX/(CPR*AINC) ! PPT FB MODS
1898 ELSE
1899 FRC2=0.
1900 ENDIF
1901 DO NK=1,LTOP
1902 QLPA(NK)=QL0(NK)
1903 QIPA(NK)=QI0(NK)
1904 QRPA(NK)=QR0(NK)
1905 QSPA(NK)=QS0(NK)
1906 RAINFB(NK)=PPTLIQ(NK)*AINC*FBFRC*FRC2 ! PPT FB MODS
1907 SNOWFB(NK)=PPTICE(NK)*AINC*FBFRC*FRC2 ! PPT FB MODS
1908 ENDDO
1909 DO NTC=1,NSTEP
1910 !
1911 !...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH LAY
1912 !...BASED ON THE SIGN OF OMEGA...
1913 !
1914 DO NK=1,LTOP
1915 QLFXIN(NK)=0.
1916 QLFXOUT(NK)=0.
1917 QIFXIN(NK)=0.
1918 QIFXOUT(NK)=0.
1919 QRFXIN(NK)=0.
1920 QRFXOUT(NK)=0.
1921 QSFXIN(NK)=0.
1922 QSFXOUT(NK)=0.
1923 ENDDO
1924 DO NK=2,LTOP
1925 IF(OMG(NK).LE.0.)THEN
1926 QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
1927 QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
1928 QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
1929 QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
1930 QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
1931 QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
1932 QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
1933 QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
1934 ELSE
1935 QLFXOUT(NK)=FXM(NK)*QLPA(NK)
1936 QIFXOUT(NK)=FXM(NK)*QIPA(NK)
1937 QRFXOUT(NK)=FXM(NK)*QRPA(NK)
1938 QSFXOUT(NK)=FXM(NK)*QSPA(NK)
1939 QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
1940 QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
1941 QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
1942 QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
1943 ENDIF
1944 ENDDO
1945 !
1946 !...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
1947 !
1948 DO NK=1,LTOP
1949 QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))*DTIME*EMSD(NK)
1950 QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))*DTIME*EMSD(NK)
1951 QRPA(NK)=QRPA(NK)+(QRFXIN(NK)-QRFXOUT(NK)+RAINFB(NK))*DTIME*EMSD(NK) ! PPT FB MODS
1952 QSPA(NK)=QSPA(NK)+(QSFXIN(NK)-QSFXOUT(NK)+SNOWFB(NK))*DTIME*EMSD(NK) ! PPT FB MODS
1953 ENDDO
1954 ENDDO
1955 DO NK=1,LTOP
1956 QLG(NK)=QLPA(NK)
1957 QIG(NK)=QIPA(NK)
1958 QRG(NK)=QRPA(NK)
1959 QSG(NK)=QSPA(NK)
1960 ENDDO
1961 !
1962 !...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
1963 !...GRID POINT...
1964 !
1965 ! IF (XTIME.LT.10.)THEN
1966 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
1967 ! ENDIF
1968 IF(IPRNT)THEN
1969 WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
1970 ! call flush(98)
1971 endif
1972 !
1973 !...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
1974 !
1975 !297 IF(IPRNT)then
1976 IF(IPRNT)then
1977 ! if(I.eq.16 .and. J.eq.41)then
1978 ! IF(ISTOP.EQ.1)THEN
1979 write(98,*)
1980 ! write(98,*)'At t(h), I, J =',float(NTSD)*72./3600.,I,J
1981 write(98,*)'P(LC), DTP, WKL, WKLCL =',p0(LC)/100., &
1982 TLCL+DTLCL+dtrh-TENV,WKL,WKLCL
1983 write(98,*)'TLCL, DTLCL, DTRH, TENV =',TLCL,DTLCL, &
1984 DTRH,TENV
1985 WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG, &
1986 TMIX-T00,PMIX,QMIX,ABE
1987 WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100., &
1988 WLCL,CLDHGT(LC)
1989 WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1990 write(98,*)'PRECIP EFFICIENCY =',PEFF
1991 WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
1992 ! ENDIF
1993 !!!!! HERE !!!!!!!
1994 WRITE(98,1070)' P ',' DP ',' DT K/D ',' DR K/D ',' OMG ', &
1995 ' DOMGDP ',' UMF ',' UER ',' UDR ',' DMF ',' DER ' &
1996 ,' DDR ',' EMS ',' W0 ',' DETLQ ',' DETIC '
1997 write(98,*)'just before DO 300...'
1998 ! call flush(98)
1999 DO NK=1,LTOP
2000 K=LTOP-NK+1
2001 DTT=(TG(K)-T0(K))*86400./TIMEC
2002 RL=XLV0-XLV1*TG(K)
2003 DR=-(QG(K)-Q0(K))*RL*86400./(TIMEC*CP)
2004 UDFRC=UDR(K)*TIMEC*EMSD(K)
2005 UEFRC=UER(K)*TIMEC*EMSD(K)
2006 DDFRC=DDR(K)*TIMEC*EMSD(K)
2007 DEFRC=-DER(K)*TIMEC*EMSD(K)
2008 WRITE(98,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)*1.E4, &
2009 UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC,DDFRC,EMS(K)/1.E11, &
2010 W0AVG1D(K)*1.E2,DETLQ(K)*TIMEC*EMSD(K)*1.E3,DETIC(K)* &
2011 TIMEC*EMSD(K)*1.E3
2012 ENDDO
2013 WRITE(98,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG', &
2014 'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
2015 DO NK=1,KL
2016 K=KX-NK+1
2017 DTT=TG(K)-T0(K)
2018 TUC=TU(K)-T00
2019 IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
2020 TDC=TZ(K)-T00
2021 IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
2022 IF(T0(K).LT.T00)THEN
2023 ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
2024 ELSE
2025 ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
2026 ENDIF
2027 QGS=ES*0.622/(P0(K)-ES)
2028 RH0=Q0(K)/QES(K)
2029 RHG=QG(K)/QGS
2030 WRITE(98,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC, &
2031 TDC,Q0(K)*1000.,QG(K)*1000.,(QG(K)-Q0(K))*1000.,QU(K)* &
2032 1000.,QD(K)*1000.,QLG(K)*1000.,QIG(K)*1000.,QRG(K)*1000., &
2033 QSG(K)*1000.,RH0,RHG
2034 ENDDO
2035 !
2036 !...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
2037 !...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
2038 !
2039 ! IF(ISTOP.EQ.1 .or. ISHALL.EQ.1)THEN
2040
2041 ! IF(ISHALL.NE.1)THEN
2042 ! write(98,4421)i,j,iyr,imo,idy,ihr,imn
2043 ! write(98)i,j,iyr,imo,idy,ihr,imn,kl
2044 ! 4421 format(7i4)
2045 ! write(98,4422)kl
2046 ! 4422 format(i6)
2047 DO 310 NK = 1,KL
2048 k = kl - nk + 1
2049 write(98,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000., &
2050 u0(k),v0(k),W0AVG1D(K),dp(k),tke(k)
2051 ! write(98) p0,t0,q0,u0,v0,w0,dp,tke
2052 ! WRITE(98,1115)Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000.,
2053 ! * U0(K),V0(K),DP(K)/100.,W0AVG(I,J,K)
2054 310 CONTINUE
2055 IF(ISTOP.EQ.1)THEN
2056 CALL wrf_error_fatal ( 'KAIN-FRITSCH, istop=1, diags' )
2057 ENDIF
2058 ! ENDIF
2059 4455 format(8f11.3)
2060 ENDIF
2061 CNDTNF=(1.-EQFRC(LFS))*(QLIQ(LFS)+QICE(LFS))*DMF(LFS)
2062 RAINCV(I,J)=DT*PPTFLX*(1.-FBFRC)/DXSQ ! PPT FB MODS
2063 ! RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ ! PPT FB MODS
2064 ! RNC=0.1*TIMEC*PPTFLX/DXSQ
2065 RNC=RAINCV(I,J)*NIC
2066 IF(ISHALL.EQ.0.AND.IPRNT)write (98,909)I,J,RNC
2067
2068 ! WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
2069 !
2070 ! EVALUATE MOISTURE BUDGET...
2071 !
2072
2073 QINIT=0.
2074 QFNL=0.
2075 DPT=0.
2076 DO 315 NK=1,LTOP
2077 DPT=DPT+DP(NK)
2078 QINIT=QINIT+Q0(NK)*EMS(NK)
2079 QFNL=QFNL+QG(NK)*EMS(NK)
2080 QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
2081 315 CONTINUE
2082 QFNL=QFNL+PPTFLX*TIMEC*(1.-FBFRC) ! PPT FB MODS
2083 ! QFNL=QFNL+PPTFLX*TIMEC ! PPT FB MODS
2084 ERR2=(QFNL-QINIT)*100./QINIT
2085 IF(IPRNT)WRITE(98,1110)QINIT,QFNL,ERR2
2086 IF(ABS(ERR2).GT.0.05 .AND. ISTOP.EQ.0)THEN
2087 ! write(99,*)'!!!!!!!! MOISTURE BUDGET ERROR IN KFPARA !!!'
2088 ! WRITE(99,1110)QINIT,QFNL,ERR2
2089 IPRNT=.TRUE.
2090 ISTOP=1
2091 write(98,4422)kl
2092 4422 format(i6)
2093 DO 311 NK = 1,KL
2094 k = kl - nk + 1
2095 ! write(99,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000., &
2096 ! u0(k),v0(k),W0AVG1D(K),dp(k)
2097 ! write(98) p0,t0,q0,u0,v0,w0,dp,tke
2098 ! WRITE(98,1115)P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
2099 ! U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
2100 WRITE(98,4456)P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
2101 U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
2102 311 CONTINUE
2103 ! call flush(98)
2104
2105 ! GOTO 297
2106 ! STOP 'QVERR'
2107 ENDIF
2108 1115 FORMAT (2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
2109 4456 format(8f12.3)
2110 IF(PPTFLX.GT.0.)THEN
2111 RELERR=ERR2*QINIT/(PPTFLX*TIMEC)
2112 ELSE
2113 RELERR=0.
2114 ENDIF
2115 IF(IPRNT)THEN
2116 WRITE(98,1120)RELERR
2117 WRITE(98,*)'TDER, CPR, TRPPT =', &
2118 TDER,CPR*AINC,TRPPT*AINC
2119 ENDIF
2120 !
2121 !...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
2122 !
2123 !...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
2124 !...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
2125 !
2126 IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
2127 NCA(I,J)=FLOAT(NIC)
2128 IF(ISHALL.EQ.1)THEN
2129 TIMEC = 2400.
2130 NCA(I,J) = FLOAT(NTST)
2131 NSHALL = NSHALL+1
2132 ENDIF
2133 DO K=1,KX
2134 ! IF(IMOIST(INEST).NE.2)THEN
2135 !
2136 !...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMAT
2137 !...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
2138 !...NOTE: THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
2139 !...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
2140 !...OF QG...
2141 !
2142 ! RLC=XLV0-XLV1*TG(K)
2143 ! RLS=XLS0-XLS1*TG(K)
2144 ! CPM=CP*(1.+0.887*QG(K))
2145 ! TG(K)=TG(K)-(RLC*(QLG(K)+QRG(K))+RLS*(QIG(K)+QSG(K)))/CPM
2146 ! QG(K)=QG(K)+(QLG(K)+QRG(K)+QIG(K)+QSG(K))
2147 ! DQLDT(I,J,NK)=0.
2148 ! DQIDT(I,J,NK)=0.
2149 ! DQRDT(I,J,NK)=0.
2150 ! DQSDT(I,J,NK)=0.
2151 ! ELSE
2152 !
2153 !...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
2154 !
2155 IF(.NOT. F_QI .and. warm_rain)THEN
2156
2157 CPM=CP*(1.+0.887*QG(K))
2158 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
2159 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
2160 DQIDT(K)=0.
2161 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
2162 DQSDT(K)=0.
2163 ELSEIF(.NOT. F_QI .and. .not. warm_rain)THEN
2164 !
2165 !...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROME
2166 !...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
2167 !
2168 CPM=CP*(1.+0.887*QG(K))
2169 IF(K.LE.ML)THEN
2170 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
2171 ELSEIF(K.GT.ML)THEN
2172 TG(K)=TG(K)+(QLG(K)+QRG(K))*RLF/CPM
2173 ENDIF
2174 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
2175 DQIDT(K)=0.
2176 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
2177 DQSDT(K)=0.
2178 ELSEIF(F_QI) THEN
2179 !
2180 !...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE TENDEN
2181 !...OF HYDROMETEORS DIRECTLY...
2182 !
2183 DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
2184 DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
2185 DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
2186 IF (F_QS) THEN
2187 DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
2188 ELSE
2189 DQIDT(K)=DQIDT(K)+(QSG(K)-QS0(K))/TIMEC
2190 ENDIF
2191 ELSE
2192 ! PRINT *,'THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED!'
2193 CALL wrf_error_fatal ( 'KAIN-FRITSCH, THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED' )
2194 ENDIF
2195 DTDT(K)=(TG(K)-T0(K))/TIMEC
2196 DQDT(K)=(QG(K)-Q0(K))/TIMEC
2197 ENDDO
2198 RAINCV(I,J)=DT*PPTFLX*(1.-FBFRC)/DXSQ ! PPT FB MODS
2199 ! RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ ! PPT FB MODS
2200 ! RNC=0.1*TIMEC*PPTFLX/DXSQ
2201 RNC=RAINCV(I,J)*NIC
2202 909 FORMAT('AT I, J =',i3,1x,i3,' CONVECTIVE RAINFALL =',F8.4,' mm')
2203 ! write (98,909)I,J,RNC
2204 ! write (6,909)I,J,RNC
2205 ! WRITE(98,*)'at NTSD =',NTSD,',No. of KF points activated =',
2206 ! * NCCNT
2207 ! call flush(98)
2208 1000 FORMAT(' ',10A8)
2209 1005 FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
2210 1010 FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
2211 1015 FORMAT(' ','ALL REMAINING MASS DETRAINS BELOW ',F4.0,' MB')
2212 1025 FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M', &
2213 ' DTLCL=',F5.2,' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=', &
2214 I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1, &
2215 ' CAPE=',0PF7.1)
2216 1030 FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =', &
2217 E12.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =', &
2218 F8.1)
2219 1035 FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'WKL=' &
2220 ,F6.3,'VWS=',F5.2)
2221 !1055 FORMAT('*** DEGREE OF STABILIZATION =',F5.3, &
2222 ! ', NO MORE MASS FLUX IS ALLOWED!')
2223 !1060 FORMAT(' ITERATION DOES NOT CONVERGE TO GIVE THE SPECIFIED &
2224 ! &DEGREE OF STABILIZATION! FABE= ',F6.4)
2225 1070 FORMAT (16A8)
2226 1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3)
2227 1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, TADVEC, NSTEP=', &
2228 2(1X,F5.0),I3,'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
2229 1085 FORMAT (A3,16A7,2A8)
2230 1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3)
2231 1095 FORMAT(' ',' PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ',F10.0)
2232 1105 FORMAT(' ','NET LATENT HEAT RELEASE =',E12.5,' ACTUAL HEATING =',&
2233 E12.5,' J/KG-S, DIFFERENCE = ',F9.3,'%')
2234 1110 FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5, &
2235 ' TOTAL WATER CHANGE =',F8.2,'%')
2236 ! 1115 FORMAT (2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
2237 1120 FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3,'%')
2238 !
2239 !-----------------------------------------------------------------------
2240 !--------------SAVE CLOUD TOP AND BOTTOM FOR RADIATION------------------
2241 !-----------------------------------------------------------------------
2242 !
2243 CUTOP(I,J)=REAL(LTOP)
2244 CUBOT(I,J)=REAL(LCL)
2245 !
2246 !-----------------------------------------------------------------------
2247 END SUBROUTINE KF_eta_PARA
2248 !********************************************************************
2249 ! ***********************************************************************
2250 SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0)
2251 !
2252 ! Lookup table variables:
2253 ! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
2254 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
2255 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
2256 ! REAL, SAVE, DIMENSION(1:200) :: ALU
2257 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
2258 ! End of Lookup table variables:
2259 !-----------------------------------------------------------------------
2260 IMPLICIT NONE
2261 !-----------------------------------------------------------------------
2262 REAL, INTENT(IN ) :: P,THES,XLV1,XLV0
2263 REAL, INTENT(OUT ) :: QNEWLQ,QNEWIC
2264 REAL, INTENT(INOUT) :: TU,QU,QLIQ,QICE
2265 REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11, &
2266 TEMP,QS,QNEW,DQ,QTOT,RLL,CPP
2267 INTEGER :: IPTB,ITHTB
2268 !-----------------------------------------------------------------------
2269
2270 !c******** LOOKUP TABLE VARIABLES... ****************************
2271 ! parameter(kfnt=250,kfnp=220)
2272 !c
2273 ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),
2274 ! * alu(200),rdpr,rdthk,plutop
2275 !C***************************************************************
2276 !c
2277 !c***********************************************************************
2278 !c scaling pressure and tt table index
2279 !c***********************************************************************
2280 !c
2281 tp=(p-plutop)*rdpr
2282 qq=tp-aint(tp)
2283 iptb=int(tp)+1
2284
2285 !
2286 !***********************************************************************
2287 ! base and scaling factor for the
2288 !***********************************************************************
2289 !
2290 ! scaling the and tt table index
2291 bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
2292 tth=(thes-bth)*rdthk
2293 pp =tth-aint(tth)
2294 ithtb=int(tth)+1
2295 IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
2296 write(98,*)'**** OUT OF BOUNDS *********'
2297 ! call flush(98)
2298 ENDIF
2299 !
2300 t00=ttab(ithtb ,iptb )
2301 t10=ttab(ithtb+1,iptb )
2302 t01=ttab(ithtb ,iptb+1)
2303 t11=ttab(ithtb+1,iptb+1)
2304 !
2305 q00=qstab(ithtb ,iptb )
2306 q10=qstab(ithtb+1,iptb )
2307 q01=qstab(ithtb ,iptb+1)
2308 q11=qstab(ithtb+1,iptb+1)
2309 !
2310 !***********************************************************************
2311 ! parcel temperature
2312 !***********************************************************************
2313 !
2314 temp=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
2315 !
2316 qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
2317 !
2318 DQ=QS-QU
2319 IF(DQ.LE.0.)THEN
2320 QNEW=QU-QS
2321 QU=QS
2322 ELSE
2323 !
2324 ! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
2325 ! ADJUSTED...IF LIQUID WATER IS PRESENT, IT IS ALLOWED TO EVAPORATE
2326 !
2327 QNEW=0.
2328 QTOT=QLIQ+QICE
2329 !
2330 ! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS AT ITS
2331 ! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MIXING
2332 ! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURATION
2333 ! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPRIATE
2334 ! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE MADE.
2335 !
2336 !...subsaturated values only occur in calculations involving various mixtures of
2337 !...updraft and environmental air for estimation of entrainment and detrainment.
2338 !...For these purposes, assume that reasonable estimates can be given using
2339 !...liquid water saturation calculations only - i.e., ignore the effect of the
2340 !...ice phase in this process only...will not affect conservative properties...
2341 !
2342 IF(QTOT.GE.DQ)THEN
2343 qliq=qliq-dq*qliq/(qtot+1.e-10)
2344 qice=qice-dq*qice/(qtot+1.e-10)
2345 QU=QS
2346 ELSE
2347 RLL=XLV0-XLV1*TEMP
2348 CPP=1004.5*(1.+0.89*QU)
2349 IF(QTOT.LT.1.E-10)THEN
2350 !
2351 !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
2352 TEMP=TEMP+RLL*(DQ/(1.+DQ))/CPP
2353 ELSE
2354 !
2355 !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURATION,
2356 ! THE TEMPERATURE IS GIVEN BY:
2357 !
2358 TEMP=TEMP+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CPP
2359 QU=QU+QTOT
2360 QTOT=0.
2361 QLIQ=0.
2362 QICE=0.
2363 ENDIF
2364 ENDIF
2365 ENDIF
2366 TU=TEMP
2367 qnewlq=qnew
2368 qnewic=0.
2369 !
2370 END SUBROUTINE TPMIX2
2371 !******************************************************************************
2372 SUBROUTINE DTFRZNEW(TU,P,THTEU,QU,QFRZ,QICE,ALIQ,BLIQ,CLIQ,DLIQ)
2373 !-----------------------------------------------------------------------
2374 IMPLICIT NONE
2375 !-----------------------------------------------------------------------
2376 REAL, INTENT(IN ) :: P,QFRZ,ALIQ,BLIQ,CLIQ,DLIQ
2377 REAL, INTENT(INOUT) :: TU,THTEU,QU,QICE
2378 REAL :: RLC,RLS,RLF,CPP,A,DTFRZ,ES,QS,DQEVAP,PII
2379 !-----------------------------------------------------------------------
2380 !
2381 !...ALLOW THE FREEZING OF LIQUID WATER IN THE UPDRAFT TO PROCEED AS AN
2382 !...APPROXIMATELY LINEAR FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE
2383 !...TTFRZ TO TBFRZ...
2384 !...FOR COLDER TERMPERATURES, FREEZE ALL LIQUID WATER...
2385 !...THERMODYNAMIC PROPERTIES ARE STILL CALCULATED WITH RESPECT TO LIQUID WATER
2386 !...TO ALLOW THE USE OF LOOKUP TABLE TO EXTRACT TMP FROM THETAE...
2387 !
2388 RLC=2.5E6-2369.276*(TU-273.16)
2389 RLS=2833922.-259.532*(TU-273.16)
2390 RLF=RLS-RLC
2391 CPP=1004.5*(1.+0.89*QU)
2392 !
2393 ! A = D(es)/DT IS THAT CALCULATED FROM BUCK (1981) EMPERICAL FORMULAS
2394 ! FOR SATURATION VAPOR PRESSURE...
2395 !
2396 A=(CLIQ-BLIQ*DLIQ)/((TU-DLIQ)*(TU-DLIQ))
2397 DTFRZ = RLF*QFRZ/(CPP+RLS*QU*A)
2398 TU = TU+DTFRZ
2399
2400 ES = ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2401 QS = ES*0.622/(P-ES)
2402 !
2403 !...FREEZING WARMS THE AIR AND IT BECOMES UNSATURATED...ASSUME THAT SOME OF THE
2404 !...LIQUID WATER THAT IS AVAILABLE FOR FREEZING EVAPORATES TO MAINTAIN SATURA-
2405 !...TION...SINCE THIS WATER HAS ALREADY BEEN TRANSFERRED TO THE ICE CATEGORY,
2406 !...SUBTRACT IT FROM ICE CONCENTRATION, THEN SET UPDRAFT MIXING RATIO AT THE NEW
2407 !...TEMPERATURE TO THE SATURATION VALUE...
2408 !
2409 DQEVAP = QS-QU
2410 QICE = QICE-DQEVAP
2411 QU = QU+DQEVAP
2412 PII=(1.E5/P)**(0.2854*(1.-0.28*QU))
2413 THTEU=TU*PII*EXP((3374.6525/TU-2.5403)*QU*(1.+0.81*QU))
2414 !
2415 END SUBROUTINE DTFRZNEW
2416 ! --------------------------------------------------------------------------------
2417
2418 SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
2419 QNEWIC,QLQOUT,QICOUT,G)
2420
2421 !-----------------------------------------------------------------------
2422 IMPLICIT NONE
2423 !-----------------------------------------------------------------------
2424 ! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
2425 ! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
2426 ! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
2427 ! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
2428 ! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
2429
2430 REAL, INTENT(IN ) :: G
2431 REAL, INTENT(IN ) :: DZ,BOTERM,ENTERM,RATE
2432 REAL, INTENT(INOUT) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
2433 REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
2434
2435 !
2436 ! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
2437 ! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
2438 ! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
2439 ! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
2440 ! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
2441 QTOT=QLIQ+QICE
2442 QNEW=QNEWLQ+QNEWIC
2443 !
2444 ! ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY
2445 ! BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
2446 ! LEVELS...
2447 !
2448 QEST=0.5*(QTOT+QNEW)
2449 G1=WTW+BOTERM-ENTERM-2.*G*DZ*QEST/1.5
2450 IF(G1.LT.0.0)G1=0.
2451 WAVG=0.5*(SQRT(WTW)+SQRT(G1))
2452 CONV=RATE*DZ/WAVG
2453 !
2454 ! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
2455 ! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
2456 ! IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
2457 ! SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...
2458 !
2459 RATIO3=QNEWLQ/(QNEW+1.E-8)
2460 ! OLDQ=QTOT
2461 QTOT=QTOT+0.6*QNEW
2462 OLDQ=QTOT
2463 RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-8)
2464 QTOT=QTOT*EXP(-CONV)
2465 !
2466 ! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
2467 ! PARCEL AT THIS LEVEL...
2468 !
2469 DQ=OLDQ-QTOT
2470 QLQOUT=RATIO4*DQ
2471 QICOUT=(1.-RATIO4)*DQ
2472 !
2473 ! ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
2474 ! LATE VERTICAL VELOCITY
2475 !
2476 PPTDRG=0.5*(OLDQ+QTOT-0.2*QNEW)
2477 WTW=WTW+BOTERM-ENTERM-2.*G*DZ*PPTDRG/1.5
2478 IF(ABS(WTW).LT.1.E-4)WTW=1.E-4
2479 !
2480 ! DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
2481 ! DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...
2482 !
2483 QLIQ=RATIO4*QTOT+RATIO3*0.4*QNEW
2484 QICE=(1.-RATIO4)*QTOT+(1.-RATIO3)*0.4*QNEW
2485 QNEWLQ=0.
2486 QNEWIC=0.
2487
2488 END SUBROUTINE CONDLOAD
2489
2490 ! ----------------------------------------------------------------------
2491 SUBROUTINE PROF5(EQ,EE,UD)
2492 !
2493 !***********************************************************************
2494 !***** GAUSSIAN TYPE MIXING PROFILE....******************************
2495 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2496 ! THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN
2497 ! DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN FROM
2498 ! "HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICS TABLES"
2499 ! ED. BY ABRAMOWITZ AND STEGUN, NATL BUREAU OF STANDARDS APPLIED
2500 ! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
2501 ! JACK KAIN
2502 ! 7/6/89
2503 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2504 !-----------------------------------------------------------------------
2505 IMPLICIT NONE
2506 !-----------------------------------------------------------------------
2507 REAL, INTENT(IN ) :: EQ
2508 REAL, INTENT(INOUT) :: EE,UD
2509 REAL :: SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
2510
2511 DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676, &
2512 0.9372980,0.33267,0.166666667,0.202765151/
2513 X=(EQ-0.5)/SIGMA
2514 Y=6.*EQ-3.
2515 EY=EXP(Y*Y/(-2))
2516 E45=EXP(-4.5)
2517 T2=1./(1.+P*ABS(Y))
2518 T1=0.500498
2519 C1=A1*T1+A2*T1*T1+A3*T1*T1*T1
2520 C2=A1*T2+A2*T2*T2+A3*T2*T2*T2
2521 IF(Y.GE.0.)THEN
2522 EE=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
2523 UD=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*(0.5+EQ*EQ/2.- &
2524 EQ)
2525 ELSE
2526 EE=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
2527 UD=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*(0.5+EQ* &
2528 EQ/2.-EQ)
2529 ENDIF
2530 EE=EE/FE
2531 UD=UD/FE
2532
2533 END SUBROUTINE PROF5
2534
2535 ! ------------------------------------------------------------------------
2536 SUBROUTINE TPMIX2DD(p,thes,ts,qs,i,j)
2537 !
2538 ! Lookup table variables:
2539 ! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
2540 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
2541 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
2542 ! REAL, SAVE, DIMENSION(1:200) :: ALU
2543 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
2544 ! End of Lookup table variables:
2545 !-----------------------------------------------------------------------
2546 IMPLICIT NONE
2547 !-----------------------------------------------------------------------
2548 REAL, INTENT(IN ) :: P,THES
2549 REAL, INTENT(INOUT) :: TS,QS
2550 INTEGER, INTENT(IN ) :: i,j ! avail for debugging
2551 REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11
2552 INTEGER :: IPTB,ITHTB
2553 CHARACTER*256 :: MESS
2554 !-----------------------------------------------------------------------
2555
2556 !
2557 !******** LOOKUP TABLE VARIABLES (F77 format)... ****************************
2558 ! parameter(kfnt=250,kfnp=220)
2559 !
2560 ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp), &
2561 ! alu(200),rdpr,rdthk,plutop
2562 !***************************************************************
2563 !
2564 !***********************************************************************
2565 ! scaling pressure and tt table index
2566 !***********************************************************************
2567 !
2568 tp=(p-plutop)*rdpr
2569 qq=tp-aint(tp)
2570 iptb=int(tp)+1
2571 !
2572 !***********************************************************************
2573 ! base and scaling factor for the
2574 !***********************************************************************
2575 !
2576 ! scaling the and tt table index
2577 bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
2578 tth=(thes-bth)*rdthk
2579 pp =tth-aint(tth)
2580 ithtb=int(tth)+1
2581 !
2582 t00=ttab(ithtb ,iptb )
2583 t10=ttab(ithtb+1,iptb )
2584 t01=ttab(ithtb ,iptb+1)
2585 t11=ttab(ithtb+1,iptb+1)
2586 !
2587 q00=qstab(ithtb ,iptb )
2588 q10=qstab(ithtb+1,iptb )
2589 q01=qstab(ithtb ,iptb+1)
2590 q11=qstab(ithtb+1,iptb+1)
2591 !
2592 !***********************************************************************
2593 ! parcel temperature and saturation mixing ratio
2594 !***********************************************************************
2595 !
2596 ts=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
2597 !
2598 qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
2599 !
2600 END SUBROUTINE TPMIX2DD
2601
2602 ! -----------------------------------------------------------------------
2603 SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,ALIQ,BLIQ,CLIQ,DLIQ)
2604 !
2605 !-----------------------------------------------------------------------
2606 IMPLICIT NONE
2607 !-----------------------------------------------------------------------
2608 REAL, INTENT(IN ) :: P1,T1,Q1,ALIQ,BLIQ,CLIQ,DLIQ
2609 REAL, INTENT(INOUT) :: THT1
2610 REAL :: EE,TLOG,ASTRT,AINC,A1,TP,VALUE,AINTRP,TDPT,TSAT,THT, &
2611 T00,P00,C1,C2,C3,C4,C5
2612 INTEGER :: INDLU
2613 !-----------------------------------------------------------------------
2614 DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834, &
2615 0.278296,1.0723E-3/
2616 !
2617 ! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
2618 !
2619 ! NOTE: Calculations for mixed/ice phase no longer used...jsk 8/00
2620 !
2621 EE=Q1*P1/(0.622+Q1)
2622 ! TLOG=ALOG(EE/ALIQ)
2623 ! ...calculate LOG term using lookup table...
2624 !
2625 astrt=1.e-3
2626 ainc=0.075
2627 a1=ee/aliq
2628 tp=(a1-astrt)/ainc
2629 indlu=int(tp)+1
2630 value=(indlu-1)*ainc+astrt
2631 aintrp=(a1-value)/ainc
2632 tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
2633 !
2634 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2635 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
2636 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
2637 THT1=THT*EXP((C1/TSAT-C2)*Q1*(1.+0.81*Q1))
2638 !
2639 END SUBROUTINE ENVIRTHT
2640 ! ***********************************************************************
2641 !====================================================================
2642 SUBROUTINE kf_eta_init(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, &
2643 RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS, &
2644 SVP1,SVP2,SVP3,SVPT0, &
2645 P_FIRST_SCALAR,restart,allowed_to_read, &
2646 ids, ide, jds, jde, kds, kde, &
2647 ims, ime, jms, jme, kms, kme, &
2648 its, ite, jts, jte, kts, kte )
2649 !--------------------------------------------------------------------
2650 IMPLICIT NONE
2651 !--------------------------------------------------------------------
2652 LOGICAL , INTENT(IN) :: restart,allowed_to_read
2653 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
2654 ims, ime, jms, jme, kms, kme, &
2655 its, ite, jts, jte, kts, kte
2656 INTEGER , INTENT(IN) :: P_QI,P_QS,P_FIRST_SCALAR
2657
2658 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
2659 RTHCUTEN, &
2660 RQVCUTEN, &
2661 RQCCUTEN, &
2662 RQRCUTEN, &
2663 RQICUTEN, &
2664 RQSCUTEN
2665
2666 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG
2667
2668 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA
2669
2670 INTEGER :: i, j, k, itf, jtf, ktf
2671 REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
2672
2673 jtf=min0(jte,jde-1)
2674 ktf=min0(kte,kde-1)
2675 itf=min0(ite,ide-1)
2676
2677 IF(.not.restart)THEN
2678
2679 DO j=jts,jtf
2680 DO k=kts,ktf
2681 DO i=its,itf
2682 RTHCUTEN(i,k,j)=0.
2683 RQVCUTEN(i,k,j)=0.
2684 RQCCUTEN(i,k,j)=0.
2685 RQRCUTEN(i,k,j)=0.
2686 ENDDO
2687 ENDDO
2688 ENDDO
2689
2690 IF (P_QI .ge. P_FIRST_SCALAR) THEN
2691 DO j=jts,jtf
2692 DO k=kts,ktf
2693 DO i=its,itf
2694 RQICUTEN(i,k,j)=0.
2695 ENDDO
2696 ENDDO
2697 ENDDO
2698 ENDIF
2699
2700 IF (P_QS .ge. P_FIRST_SCALAR) THEN
2701 DO j=jts,jtf
2702 DO k=kts,ktf
2703 DO i=its,itf
2704 RQSCUTEN(i,k,j)=0.
2705 ENDDO
2706 ENDDO
2707 ENDDO
2708 ENDIF
2709
2710 DO j=jts,jtf
2711 DO i=its,itf
2712 NCA(i,j)=-100.
2713 ENDDO
2714 ENDDO
2715
2716 DO j=jts,jtf
2717 DO k=kts,ktf
2718 DO i=its,itf
2719 W0AVG(i,k,j)=0.
2720 ENDDO
2721 ENDDO
2722 ENDDO
2723
2724 endif
2725
2726 CALL KF_LUTAB(SVP1,SVP2,SVP3,SVPT0)
2727
2728 END SUBROUTINE kf_eta_init
2729
2730 !-------------------------------------------------------
2731
2732 subroutine kf_lutab(SVP1,SVP2,SVP3,SVPT0)
2733 !
2734 ! This subroutine is a lookup table.
2735 ! Given a series of series of saturation equivalent potential
2736 ! temperatures, the temperature is calculated.
2737 !
2738 !--------------------------------------------------------------------
2739 IMPLICIT NONE
2740 !--------------------------------------------------------------------
2741 ! Lookup table variables
2742 ! INTEGER, SAVE, PARAMETER :: KFNT=250,KFNP=220
2743 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
2744 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
2745 ! REAL, SAVE, DIMENSION(1:200) :: ALU
2746 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
2747 ! End of Lookup table variables
2748
2749 INTEGER :: KP,IT,ITCNT,I
2750 REAL :: DTH,TMIN,TOLER,PBOT,DPR, &
2751 TEMP,P,ES,QS,PI,THES,TGUES,THGUES,F0,T1,T0,THGS,F1,DT, &
2752 ASTRT,AINC,A1,THTGS
2753 ! REAL :: ALIQ,BLIQ,CLIQ,DLIQ,SVP1,SVP2,SVP3,SVPT0
2754 REAL :: ALIQ,BLIQ,CLIQ,DLIQ
2755 REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
2756 !
2757 ! equivalent potential temperature increment
2758 data dth/1./
2759 ! minimum starting temp
2760 data tmin/150./
2761 ! tolerance for accuracy of temperature
2762 data toler/0.001/
2763 ! top pressure (pascals)
2764 plutop=5000.0
2765 ! bottom pressure (pascals)
2766 pbot=110000.0
2767
2768 ALIQ = SVP1*1000.
2769 BLIQ = SVP2
2770 CLIQ = SVP2*SVPT0
2771 DLIQ = SVP3
2772
2773 !
2774 ! compute parameters
2775 !
2776 ! 1._over_(sat. equiv. theta increment)
2777 rdthk=1./dth
2778 ! pressure increment
2779 !
2780 DPR=(PBOT-PLUTOP)/REAL(KFNP-1)
2781 ! dpr=(pbot-plutop)/REAL(kfnp-1)
2782 ! 1._over_(pressure increment)
2783 rdpr=1./dpr
2784 ! compute the spread of thes
2785 ! thespd=dth*(kfnt-1)
2786 !
2787 ! calculate the starting sat. equiv. theta
2788 !
2789 temp=tmin
2790 p=plutop-dpr
2791 do kp=1,kfnp
2792 p=p+dpr
2793 es=aliq*exp((bliq*temp-cliq)/(temp-dliq))
2794 qs=0.622*es/(p-es)
2795 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
2796 the0k(kp)=temp*pi*exp((3374.6525/temp-2.5403)*qs* &
2797 (1.+0.81*qs))
2798 enddo
2799 !
2800 ! compute temperatures for each sat. equiv. potential temp.
2801 !
2802 p=plutop-dpr
2803 do kp=1,kfnp
2804 thes=the0k(kp)-dth
2805 p=p+dpr
2806 do it=1,kfnt
2807 ! define sat. equiv. pot. temp.
2808 thes=thes+dth
2809 ! iterate to find temperature
2810 ! find initial guess
2811 if(it.eq.1) then
2812 tgues=tmin
2813 else
2814 tgues=ttab(it-1,kp)
2815 endif
2816 es=aliq*exp((bliq*tgues-cliq)/(tgues-dliq))
2817 qs=0.622*es/(p-es)
2818 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
2819 thgues=tgues*pi*exp((3374.6525/tgues-2.5403)*qs* &
2820 (1.+0.81*qs))
2821 f0=thgues-thes
2822 t1=tgues-0.5*f0
2823 t0=tgues
2824 itcnt=0
2825 ! iteration loop
2826 do itcnt=1,11
2827 es=aliq*exp((bliq*t1-cliq)/(t1-dliq))
2828 qs=0.622*es/(p-es)
2829 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
2830 thtgs=t1*pi*exp((3374.6525/t1-2.5403)*qs*(1.+0.81*qs))
2831 f1=thtgs-thes
2832 if(abs(f1).lt.toler)then
2833 exit
2834 endif
2835 ! itcnt=itcnt+1
2836 dt=f1*(t1-t0)/(f1-f0)
2837 t0=t1
2838 f0=f1
2839 t1=t1-dt
2840 enddo
2841 ttab(it,kp)=t1
2842 qstab(it,kp)=qs
2843 enddo
2844 enddo
2845 !
2846 ! lookup table for tlog(emix/aliq)
2847 !
2848 ! set up intial values for lookup tables
2849 !
2850 astrt=1.e-3
2851 ainc=0.075
2852 !
2853 a1=astrt-ainc
2854 do i=1,200
2855 a1=a1+ainc
2856 alu(i)=alog(a1)
2857 enddo
2858 !
2859 END SUBROUTINE KF_LUTAB
2860
2861 END MODULE module_cu_kfeta