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