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