6 ! V3.3: A new trigger function is added based Ma and Tan (2009):
7 ! Ma, L.-M. and Z.-M. Tan, 2009: Improving the behavior of
8 ! the cumulus parameterization for tropical cyclone prediction:
9 ! Convection trigger. Atmospheric Research, 92, 190 - 211.
11 !--------------------------------------------------------------------
12 ! Lookup table variables:
13 INTEGER, PARAMETER :: KFNT=250,KFNP=220
14 REAL, DIMENSION(KFNT,KFNP),PRIVATE, SAVE :: TTAB,QSTAB
15 REAL, DIMENSION(KFNP),PRIVATE, SAVE :: THE0K
16 REAL, DIMENSION(200),PRIVATE, SAVE :: ALU
17 REAL, PRIVATE, SAVE :: RDPR,RDTHK,PLUTOP
18 ! Note: KF Lookup table is used by subroutines KF_eta_PARA, TPMIX2,
20 ! End of Lookup table variables:
24 SUBROUTINE KF_eta_CPS( &
25 ids,ide, jds,jde, kds,kde &
26 ,ims,ime, jms,jme, kms,kme &
27 ,its,ite, jts,jte, kts,kte &
29 ,DT,KTAU,DX,CUDT,CURR_SECS,ADAPT_STEP_FLAG &
31 ,rho,RAINCV,PRATEC,NCA &
32 ,U,V,TH,T,W,dz8w,Pcps,pi &
33 ,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
34 ,EP2,SVP1,SVP2,SVP3,SVPT0 &
35 ,STEPCU,CU_ACT_FLAG,warm_rain,CUTOP,CUBOT &
38 ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
39 ,RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN &
40 ,RQICUTEN,RQSCUTEN, RQVFTEN &
43 !-------------------------------------------------------------
45 !-------------------------------------------------------------
46 INTEGER, INTENT(IN ) :: &
47 ids,ide, jds,jde, kds,kde, &
48 ims,ime, jms,jme, kms,kme, &
49 its,ite, jts,jte, kts,kte
51 INTEGER, INTENT(IN ) :: trigger
52 INTEGER, INTENT(IN ) :: STEPCU
53 LOGICAL, INTENT(IN ) :: warm_rain
55 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1
56 REAL, INTENT(IN ) :: CP,R,G,EP1,EP2
57 REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
59 INTEGER, INTENT(IN ) :: KTAU
61 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
74 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
78 REAL, INTENT(IN ) :: DT, DX
79 REAL, INTENT(IN ) :: CUDT
80 REAL, INTENT(IN ) :: CURR_SECS
81 LOGICAL,OPTIONAL,INTENT(IN ) :: ADAPT_STEP_FLAG
82 REAL, INTENT (INOUT) :: CUDTACTTIME
84 REAL, DIMENSION( ims:ime , jms:jme ), &
85 INTENT(INOUT) :: RAINCV
87 REAL, DIMENSION( ims:ime , jms:jme ), &
88 INTENT(INOUT) :: PRATEC
90 REAL, DIMENSION( ims:ime , jms:jme ), &
93 REAL, DIMENSION( ims:ime , jms:jme ), &
94 INTENT(OUT) :: CUBOT, &
97 LOGICAL, DIMENSION( ims:ime , jms:jme ), &
98 INTENT(INOUT) :: CU_ACT_FLAG
104 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
115 ! Flags relating to the optional tendency arrays declared above
116 ! Models that carry the optional tendencies will provdide the
117 ! optional arguments at compile time; these flags all the model
118 ! to determine at run-time whether a particular tracer is in
121 LOGICAL, OPTIONAL :: &
131 LOGICAL :: flag_qr, flag_qi, flag_qs
133 REAL, DIMENSION( kts:kte ) :: &
145 REAL, DIMENSION( kts:kte ):: &
153 REAL, DIMENSION (its-1:ite+1,kts:kte,jts-1:jte+1) :: aveh_t, aveh_q
154 REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: aveh_qmax, aveh_qmin
155 REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: avev_t, avev_q
156 REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: avev_qmax, avev_qmin
157 REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: coef_v, coef_h, tpart_h, tpart_v
161 REAL, DIMENSION (kts:kte) :: z0
163 REAL :: TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp
164 integer :: ibegh,iendh,jbegh,jendh
165 integer :: istart,iend,jstart,jend
166 INTEGER :: i,j,k,NTST
167 REAL :: lastdt = -1.0
168 REAL :: W0AVGfctr, W0fctr, W0den
169 LOGICAL :: run_param , doing_adapt_dt , decided
174 !----------------------
180 IF ( PRESENT(F_QR) ) flag_qr = F_QR
181 IF ( PRESENT(F_QI) ) flag_qi = F_QI
182 IF ( PRESENT(F_QS) ) flag_qs = F_QS
188 if (ADAPT_STEP_FLAG) then
189 W0AVGfctr = 2 * MAX(CUDT*60,dt) - dt
191 W0den = 2 * MAX(CUDT*60,dt)
201 ! SCR1=-5.0E-4*G*rho(I,K,J)*(w(I,K,J)+w(I,K+1,J))
202 ! TV=T(I,K,J)*(1.+EP1*QV(I,K,J))
203 ! RHOE=Pcps(I,K,J)/(R*TV)
204 ! W0=-101.9368*SCR1/RHOE
205 W0=0.5*(w(I,K,J)+w(I,K+1,J))
209 ! W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST
211 ! New, to support adaptive time step:
213 W0AVG(I,K,J) = ( W0AVG(I,K,J) * W0AVGfctr + W0 * W0fctr ) / W0den
219 ! Initialization for adaptive time step.
221 doing_adapt_dt = .FALSE.
222 IF ( PRESENT(adapt_step_flag) ) THEN
223 IF ( adapt_step_flag ) THEN
224 doing_adapt_dt = .TRUE.
225 IF ( cudtacttime .EQ. 0. ) THEN
226 cudtacttime = curr_secs + cudt*60.
231 ! Do we run through this scheme or not?
233 ! Test 1: If this is the initial model time, then yes.
235 ! Test 2: If the user asked for the cumulus to be run every time step, then yes.
237 ! Test 3: If not adaptive dt, and this is on the requested cumulus frequency, then yes.
239 ! Test 4: If using adaptive dt and the current time is past the last requested activate cumulus time, then yes.
240 ! CURR_SECS >= CUDTACTTIME
242 ! If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag
243 ! to TRUE. The decided flag says that one of these tests was able to say "yes", run the scheme.
244 ! We only proceed to other tests if the previous tests all have left decided as FALSE.
246 ! If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next
251 IF ( ( .NOT. decided ) .AND. &
252 ( ktau .EQ. 1 ) ) THEN
257 IF ( ( .NOT. decided ) .AND. &
258 ( ( cudt .EQ. 0. ) .OR. ( stepcu .EQ. 1 ) ) ) THEN
263 IF ( ( .NOT. decided ) .AND. &
264 ( .NOT. doing_adapt_dt ) .AND. &
265 ( MOD(ktau,ntst) .EQ. 0 ) ) THEN
270 IF ( ( .NOT. decided ) .AND. &
271 ( doing_adapt_dt ) .AND. &
272 ( curr_secs .GE. cudtacttime ) ) THEN
275 cudtacttime = curr_secs + cudt*60
280 ! New trigger function
281 IF (trigger.eq.2) THEN
283 ! calculate 9-point average of moisture advection and temperature using halo (Horizontal)
285 aveh_t=-999 ! horizontal 9-point ave
287 avev_t=0 ! vertical 3-level ave
297 ibegh=max(its-1, ids+1) ! start from 2
298 jbegh=max(jts-1, jds+1)
299 iendh=min(ite+1, ide-2) ! end at ide-2
300 jendh=min(jte+1, jde-2)
304 aveh_t(i,k,j)=(T(i-1,k,j-1)+T(i-1,k,j) +T(i-1,k,j+1)+ &
305 T(i,k,j-1) +T(i,k,j) +T(i,k,j+1)+ &
306 T(i+1,k,j-1) +T(i+1,k,j) +T(i+1,k,j+1))/9.
307 aveh_q(i,k,j)=(rqvften(i-1,k,j-1)+rqvften(i-1,k,j) +rqvften(i-1,k,j+1)+ &
308 rqvften(i,k,j-1) +rqvften(i,k,j) +rqvften(i,k,j+1)+ &
309 rqvften(i+1,k,j-1) +rqvften(i+1,k,j) +rqvften(i+1,k,j+1))/9.
313 ! boundary value ( all processors will do the following? Or just those processsors handling sub-area including boundary)
319 aveh_t(i,k,j)=aveh_t(i+1,k,j)
320 aveh_q(i,k,j)=aveh_q(i+1,k,j)
321 elseif(i.eq.ide-1) then
322 aveh_t(i,k,j)=aveh_t(i-1,k,j)
323 aveh_q(i,k,j)=aveh_q(i-1,k,j)
327 aveh_t(i,k,j)=aveh_t(i,k,j+1)
328 aveh_q(i,k,j)=aveh_q(i,k,j+1)
329 elseif(j.eq.jde-1) then
330 aveh_t(i,k,j)=aveh_t(i,k,j-1)
331 aveh_q(i,k,j)=aveh_q(i,k,j-1)
334 if(j.eq.jds.and.i.eq.ids) then
335 aveh_q(i,k,j)=aveh_q(i+1,k,j+1)
336 aveh_t(i,k,j)=aveh_t(i+1,k,j+1)
339 if(j.eq.jde-1.and.i.eq.ids) then
340 aveh_q(i,k,j)=aveh_q(i+1,k,j-1)
341 aveh_t(i,k,j)=aveh_t(i+1,k,j-1)
344 if(j.eq.jde-1.and.i.eq.ide-1) then
345 aveh_q(i,k,j)=aveh_q(i-1,k,j-1)
346 aveh_t(i,k,j)=aveh_t(i-1,k,j-1)
349 if(j.eq.jds.and.i.eq.ide-1) then
350 aveh_q(i,k,j)=aveh_q(i-1,k,j+1)
351 aveh_t(i,k,j)=aveh_t(i-1,k,j+1)
357 ! search for max/min moisture advection in 9-point range, calculate horizontal T-perturbation (tpart_h)
358 istart=max(its, ids+1) ! start from 2
359 jstart=max(jts, jds+1)
360 iend=min(ite, ide-2) ! end at ide-2
365 aveh_qmax(i,k,j)=aveh_q(i,k,j)
366 aveh_qmin(i,k,j)=aveh_q(i,k,j)
369 if(aveh_q(i+II,k,j+JJ).gt.aveh_qmax(i,k,j)) aveh_qmax(i,k,j)=aveh_q(i+II,k,j+JJ)
370 if(aveh_q(i+II,k,j+JJ).lt.aveh_qmin(i,k,j)) aveh_qmin(i,k,j)=aveh_q(i+II,k,j+JJ)
373 if(aveh_qmax(i,k,j).gt.aveh_qmin(i,k,j))then
374 coef_h(i,k,j)=(aveh_q(i,k,j)-aveh_qmin(i,k,j))/(aveh_qmax(i,k,j)-aveh_qmin(i,k,j))
378 coef_h(i,k,j)=amin1(coef_h(i,k,j),1.0)
379 coef_h(i,k,j)=amax1(coef_h(i,k,j),0.0)
380 tpart_h(i,k,j)=coef_h(i,k,j)*(T(i,k,j)-aveh_t(i,k,j))
385 ! vertical 3-layer calculation
388 z0(1) = 0.5 * dz8w(i,1,j)
390 Z0(K) = Z0(K-1) + .5 * (DZ8W(i,K,j) + DZ8W(i,K-1,j))
393 ttop = t(i,k,j) + ((t(i,k,j) - t(i,k+1,j)) / (z0(k) - z0(k+1))) * (z0(k)-z0(k-1))
394 avev_t(i,k,j)=(T(i,k-1,j) + T(i,k,j) + ttop)/3.
395 ! avev_t(i,k,j)=(T(i,k-1,j)+T(i,k,j) + T(i,k+1,j))/3.
396 avev_q(i,k,j)=(rqvften(i,k-1,j)+rqvften(i,k,j) + rqvften(i,k+1,j))/3.
398 avev_t(i,kts,j)=avev_t(i,kts+1,j) ! lowest level value, is it the same as avev_t(i,kds,j)=avev_t(i,kds+1,j)?
399 avev_q(i,kts,j)=avev_q(i,kts+1,j)
400 avev_t(i,kte,j)=avev_t(i,kte-1,j) ! highest level value
401 avev_q(i,kte,j)=avev_q(i,kte-1,j)
408 avev_qmax(i,k,j)=avev_q(i,k,j)
409 avev_qmin(i,k,j)=avev_q(i,k,j)
411 if(avev_q(i,k+kk,j).gt.avev_qmax(i,k,j)) avev_qmax(i,k,j)=avev_q(i,k+kk,j)
412 if(avev_q(i,k+kk,j).lt.avev_qmin(i,k,j)) avev_qmin(i,k,j)=avev_q(i,k+kk,j)
414 if(avev_qmax(i,k,j).gt.avev_qmin(i,k,j)) then
415 coef_v(i,k,j)=(avev_q(i,k,j)-avev_qmin(i,k,j))/(avev_qmax(i,k,j)-avev_qmin(i,k,j))
419 tpart_v(i,k,j)=coef_v(i,k,j)*(T(i,k,j)-avev_t(i,k,j))
421 tpart_v(i,kts,j)= tpart_v(i,kts+1,j) ! lowest level
422 tpart_v(i,kte,j)= tpart_v(i,kte-1,j) ! highest level
425 ENDIF ! endif (trigger.eq.2)
429 CU_ACT_FLAG(i,j) = .true.
437 IF ( NCA(I,J) .ge. 0.5*DT ) then
438 CU_ACT_FLAG(i,j) = .false.
454 ! assign vars from 3D to 1D
463 W0AVG1D(K) =W0AVG(I,K,J)
465 IF (trigger.eq.2) THEN
466 tpart_h1D(K) =tpart_h(I,K,J)
467 tpart_v1D(K) =tpart_v(I,K,J)
473 CALL KF_eta_PARA(I, J, &
474 U1D,V1D,T1D,QV1D,P1D,DZ1D,W0AVG1D, &
475 tpart_h1D,tpart_v1D, &
478 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
479 EP2,SVP1,SVP2,SVP3,SVPT0, &
480 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
482 flag_QI,flag_QS,warm_rain, &
484 ids,ide, jds,jde, kds,kde, &
485 ims,ime, jms,jme, kms,kme, &
486 its,ite, jts,jte, kts,kte)
487 IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
489 RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
490 RQVCUTEN(I,K,J)=DQDT(K)
494 IF(PRESENT(rqrcuten).AND.PRESENT(rqccuten)) THEN
497 RQRCUTEN(I,K,J)=DQRDT(K)
498 RQCCUTEN(I,K,J)=DQCDT(K)
501 ! This is the case for Eta microphysics without 3d rain field
504 RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
509 !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
512 IF(PRESENT( rqicuten )) THEN
515 RQICUTEN(I,K,J)=DQIDT(K)
520 IF(PRESENT( rqscuten )) THEN
523 RQSCUTEN(I,K,J)=DQSDT(K)
533 END SUBROUTINE KF_eta_CPS
534 ! ****************************************************************************
535 !-----------------------------------------------------------
536 SUBROUTINE KF_eta_PARA (I, J, &
537 U0,V0,T0,QV0,P0,DZQ,W0AVG1D, &
541 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
542 EP2,SVP1,SVP2,SVP3,SVPT0, &
543 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
545 F_QI,F_QS,warm_rain, &
547 ids,ide, jds,jde, kds,kde, &
548 ims,ime, jms,jme, kms,kme, &
549 its,ite, jts,jte, kts,kte)
550 !-----------------------------------------------------------
551 !***** The KF scheme that is currently used in experimental runs of EMCs
552 !***** Eta model....jsk 8/00
555 !-----------------------------------------------------------
556 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
557 ims,ime, jms,jme, kms,kme, &
558 its,ite, jts,jte, kts,kte, &
560 ! ,P_QI,P_QS,P_FIRST_SCALAR
561 INTEGER, INTENT(IN ) :: trigger
563 LOGICAL, INTENT(IN ) :: F_QI, F_QS
565 LOGICAL, INTENT(IN ) :: warm_rain
567 REAL, DIMENSION( kts:kte ), &
579 REAL, INTENT(IN ) :: DT,DX,DXSQ
582 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
583 REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
586 REAL, DIMENSION( kts:kte ), INTENT(INOUT) :: &
594 REAL, DIMENSION( ims:ime , jms:jme ), &
597 REAL, DIMENSION( ims:ime , jms:jme ), &
598 INTENT(INOUT) :: RAINCV
600 REAL, DIMENSION( ims:ime , jms:jme ), &
601 INTENT(INOUT) :: PRATEC
603 REAL, DIMENSION( ims:ime , jms:jme ), &
604 INTENT(OUT) :: CUBOT, &
606 REAL, INTENT(IN ) :: CUDT
608 !...DEFINE LOCAL VARIABLES...
610 REAL, DIMENSION( kts:kte ) :: &
611 Q0,Z0,TV0,TU,TVU,QU,TZ,TVD, &
612 QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD, &
613 UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2, &
614 UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE, &
615 THTAU,THETEU,THTAD,THETED,QLIQ,QICE, &
616 QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC, &
617 DETLQ2,DETIC2,RATIO,RATIO2
620 REAL, DIMENSION( kts:kte ) :: &
621 DOMGDP,EXN,TVQU,DP,RH,EQFRC,WSPD, &
622 QDT,FXM,THTAG,THPA,THFXOUT, &
623 THFXIN,QPA,QFXOUT,QFXIN,QLPA,QLFXIN, &
624 QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA, &
625 QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT, &
626 QL0,QLG,QI0,QIG,QR0,QRG,QS0,QSG
629 REAL, DIMENSION( kts:kte+1 ) :: OMG
630 REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
631 REAL, DIMENSION( kts:kte ) :: &
632 CLDHGT,QSD,DILFRC,DDILFRC,TKE,TGU,QGU,THTEEG
636 REAL :: P00,T00,RLF,RHIC,RHBC,PIE, &
638 REAL :: GDRY,ROCP,ALIQ,BLIQ, &
640 REAL :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX, &
641 ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL, &
642 CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR, &
643 ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
644 TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD, &
645 UPNEW,ABE,WKLCL,TTEMP,FRC1, &
646 QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
647 DZZ,UDLBE,REI,EE2,UD2,TTMP,F1,F2, &
648 THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1, &
649 UD1,DPTT,QNEWLQ,DUMFDP,EE,TSAT, &
650 THTA,VCONV,TIMEC,SHSIGN,VWS,PEF, &
651 CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN, &
652 DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1, &
653 DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF, &
654 UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF, &
655 DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
656 AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1, &
657 DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF, &
658 TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR, &
659 UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2, &
660 RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
661 DDFRC,TDC,DEFRC,RHBAR,DMFFRC,DPMIN,DILBE
662 REAL :: ASTRT,TP,VALUE,AINTRP,TKEMAX,QFRZ,&
663 QSS,PPTMLT,DTMELT,RHH,EVAC,BINC
665 INTEGER :: INDLU,NU,NUCHM,NNN,KLFS
666 REAL :: CHMIN,PM15,CHMAX,DTRH,RAD,DPPP
667 REAL :: TVDIFF,DTTOT,ABSOMG,ABSOMGTC,FRDP
672 INTEGER, DIMENSION (kts:kte) :: KCHECK
674 INTEGER :: ISTOP,ML,L5,KMIX,LOW, &
675 LC,MXLAYR,LLFC,NLAYRS,NK, &
676 KPBL,KLCL,LCL,LET,IFLAG, &
678 LTOPM1,LVF,KSTART,KMIN,LFS, &
679 ND,NIC,LDB,LDT,ND1,NDK, &
680 NM,LMAX,NCOUNT,NOITR, &
681 NSTEP,NTC,NCHM,ISHALL,NSHALL
683 REAL :: u00,qslcl,rhlcl,dqssdt !jfb
684 CHARACTER*1024 message
686 DATA P00,T00/1.E5,273.16/
688 DATA RHIC,RHBC/1.,0.90/
689 DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
690 DATA RATE/0.03/ ! wrf default
691 ! DATA RATE/0.01/ ! value used in NRCM
692 ! DATA RATE/0.001/ ! effectively turn off autoconversion
693 !-----------------------------------------------------------
711 !****************************************************************************
713 !...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER ! PPT FB MODS
714 !...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL) ! PPT FB MODS
715 !...FIELD. "FBFRC" IS THE FRACTION OF AVAILABLE ! PPT FB MODS
716 !...PRECIPITATION TO BE FED BACK (0.0 - 1.0)... ! PPT FB MODS
717 FBFRC=0.0 ! PPT FB MODS
718 !...mods to allow shallow convection...
725 !...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
726 !...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
727 !...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
729 !...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
730 !...FROM BOTTOM-UP IN THE KF SCHEME...
733 !SUE tmprpsb=1./PSB(I,J)
734 !SUE CELL=PTOP*tmprpsb
738 ! Saturation vapor pressure (ES) is calculated following Buck (1981)
739 !...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
741 ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
742 QES(K)=0.622*ES/(P0(K)-ES)
743 Q0(K)=AMIN1(QES(K),QV0(K))
744 Q0(K)=AMAX1(0.000001,Q0(K))
751 TV0(K)=T0(K)*(1.+0.608*Q0(K))
752 ! RHOE(K)=P0(K)/(R*TV0(K))
753 ! DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
754 DP(K)=rhoe(k)*g*DZQ(k)
755 ! IF Turbulent Kinetic Energy (TKE) is available from turbulent mixing scheme
756 ! use it for shallow convection...For now, assume it is not available....
757 ! TKE(K) = Q2(I,J,NK)
760 ! IF(P0(K).GE.500E2)L5=K
761 IF(P0(K).GE.0.5*P0(1))L5=K
762 IF(P0(K).GE.P300)LLFC=K
765 !...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
769 Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
770 DZA(K-1)=Z0(K)-Z0(K-1)
775 ! To save time, specify a pressure interval to move up in sequential
776 ! check of different ~50 mb deep groups of adjacent model layers in
777 ! the process of identifying updraft source layer (USL). Note that
778 ! this search is terminated as soon as a buoyant parcel is found and
779 ! this parcel can produce a cloud greater than specifed minimum depth
780 ! (CHMIN)...For now, set interval at 15 mb...
786 IF(P0(K).LT.PM15)THEN
803 IF(CLDHGT(NNN).GT.CHMAX)THEN
821 !...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
822 !...UNSTABLE AIR AT LEAST 50 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
823 !...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
824 !...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 50 mb..
829 IF ( NK+1 .LT. KTS ) THEN
830 WRITE(message,*)'WOULD GO OFF BOTTOM: KF_ETA_PARA I,J,NK',I,J,NK
831 CALL wrf_message (TRIM(message))
835 IF ( NK .GT. KTE ) THEN
836 WRITE(message,*)'WOULD GO OFF TOP: KF_ETA_PARA I,J,DPTHMX,DPMIN',I,J,DPTHMX,DPMIN
837 CALL wrf_message (TRIM(message))
842 IF(DPTHMX.GT.DPMIN)THEN
847 IF(DPTHMX.LT.DPMIN)THEN
852 !...********************************************************
853 !...for computational simplicity without much loss in accuracy,
854 !...mix temperature instead of theta for evaluating convective
855 !...initiation (triggering) potential...
862 !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
863 !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
868 TMIX=TMIX+DP(NK)*T0(NK)
869 QMIX=QMIX+DP(NK)*Q0(NK)
870 ZMIX=ZMIX+DP(NK)*Z0(NK)
871 PMIX=PMIX+DP(NK)*P0(NK)
878 EMIX=QMIX*PMIX/(0.622+QMIX)
880 !...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL...
882 ! TLOG=ALOG(EMIX/ALIQ)
883 ! ...calculate dewpoint using lookup table...
890 value=(indlu-1)*ainc+astrt
891 aintrp=(a1-value)/ainc
892 tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
893 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
894 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
895 TLCL=AMIN1(TLCL,TMIX)
896 TVLCL=TLCL*(1.+0.608*QMIX)
897 ZLCL = ZMIX+(TLCL-TMIX)/GDRY
902 IF(ZLCL.LE.Z0(NK) .or. NK.GT.KL)THEN
910 ! calculate DLP using Z instead of log(P)
911 DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
913 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
915 TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
916 QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
917 TVEN=TENV*(1.+0.608*QENV)
919 !...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
920 !...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0 IS AN
921 !...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
922 !...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
923 !...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
924 !...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
925 !...SUCCESS AT GRID LENGTHS NEAR 25 km. FOR DIFFERENT GRID-LENGTHS,
926 !...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
927 !...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
929 IF(ZLCL.LT.2.E3)THEN ! Kain (2004) Eq. 2
932 WKLCL=0.02 ! units of m/s
934 WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)*DX/25.E3-WKLCL
935 IF(WKL.LT.0.0001)THEN
938 DTLCL=4.64*WKL**0.33 ! Kain (2004) Eq. 1
941 ! New trigger function
942 IF(trigger.eq.2) then
943 DTLCL=amax1(TPART_H0(KLCL)+TPART_V0(KLCL),0.0)
945 ! end for trigger function
948 if (trigger .eq. 3) then
949 !...for ETA model, give parcel an extra temperature perturbation based
950 !...the threshold RH for condensation (U00)...
951 ! as described in Narita and Ohmori (2007, 12th Mesoscale Conf.
953 !...for now, just assume U00=0.75...
954 !...!!!!!! for MM5, SET DTRH = 0. !!!!!!!!
957 QSLCL=QES(K)+(QES(KLCL)-QES(K))*DLP
959 DQSSDT = QMIX*(CLIQ-BLIQ*DLIQ)/((TLCL-DLIQ)*(TLCL-DLIQ))
960 IF(RHLCL.ge.0.75 .and. RHLCL.le.0.95)then
961 DTRH = 0.25*(RHLCL-0.75)*QMIX/DQSSDT
962 ELSEIF(RHLCL.GT.0.95)THEN
963 DTRH = (1./RHLCL-1.)*QMIX/DQSSDT
969 ! IF(ISHALL.EQ.1)IPRNT=.TRUE.
971 ! IF(TLCL+DTLCL.GT.TENV)GOTO 45
973 trigger2: IF(TLCL+DTLCL+DTRH.LT.TENV)THEN
975 ! Parcel not buoyant, CYCLE back to start of trigger and evaluate next potential USL...
979 ELSE ! Parcel is buoyant, determine updraft
981 !...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
982 !...EQUIVALENT POTENTIAL TEMPERATURE
983 !...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
985 CALL ENVIRTHT(PMIX,TMIX,QMIX,THETEU(K),ALIQ,BLIQ,CLIQ,DLIQ)
987 !...modify calculation of initial parcel vertical velocity...jsk 11/26/97
990 IF(DTTOT.GT.1.E-4)THEN
991 GDT=2.*G*DTTOT*500./TVEN ! Kain (2004) Eq. 3 (sort of)
992 WLCL=1.+0.5*SQRT(GDT)
993 WLCL = AMIN1(WLCL,3.)
997 PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
1000 TVLCL=TLCL*(1.+0.608*QMIX)
1001 RHOLCL=PLCL/(R*TVLCL)
1005 ! make RAD a function of background vertical velocity... (Kain (2004) Eq. 6)
1008 ELSEIF(WKL.GT.0.1)THEN
1011 RAD = 1000.+1000*WKL/0.1
1014 !*******************************************************************
1016 ! COMPUTE UPDRAFT PROPERTIES *
1018 !*******************************************************************
1022 !...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
1031 !...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
1032 !...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE
1033 !...BUOYANT ENERGY, TRPPT IS THE TOTAL RATE OF PRECIPITATION
1054 !...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
1055 !...PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
1056 !...FREEZING IS SPECIFIED TO BEGIN. WITHIN THE GLACIATION
1057 !...INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
1058 !...PREVIOUS MODEL LEVEL...
1062 !...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
1063 !...MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
1064 !...MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
1066 ! **1 variables indicate the bottom of a model layer
1067 ! **2 variables indicate the top of a model layer
1073 updraft: DO NK=K,KL-1
1075 RATIO2(NK1)=RATIO2(NK)
1078 THETEU(NK1)=THETEU(NK)
1082 call tpmix2(p0(nk1),theteu(nk1),tu(nk1),qu(nk1),qliq(nk1), &
1083 qice(nk1),qnewlq,qnewic,XLV1,XLV0)
1086 !...CHECK TO SEE IF UPDRAFT TEMP IS ABOVE THE TEMPERATURE AT WHICH
1087 !...GLACIATION IS ASSUMED TO INITIATE; IF IT IS, CALCULATE THE
1088 !...FRACTION OF REMAINING LIQUID WATER TO FREEZE...TTFRZ IS THE
1089 !...TEMP AT WHICH FREEZING BEGINS, TBFRZ THE TEMP BELOW WHICH ALL
1090 !...LIQUID WATER IS FROZEN AT EACH LEVEL...
1092 IF(TU(NK1).LE.TTFRZ)THEN
1093 IF(TU(NK1).GT.TBFRZ)THEN
1094 IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
1095 FRC1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
1102 ! DETERMINE THE EFFECTS OF LIQUID WATER FREEZING WHEN TEMPERATURE
1103 !...IS BELOW TTFRZ...
1105 QFRZ = (QLIQ(NK1)+QNEWLQ)*FRC1
1106 QNEWIC=QNEWIC+QNEWLQ*FRC1
1107 QNEWLQ=QNEWLQ-QNEWLQ*FRC1
1108 QICE(NK1) = QICE(NK1)+QLIQ(NK1)*FRC1
1109 QLIQ(NK1) = QLIQ(NK1)-QLIQ(NK1)*FRC1
1110 CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QFRZ, &
1111 QICE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
1113 TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
1115 ! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
1118 BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
1119 BOTERM=2.*(Z0(NK1)-ZLCL)*G*BE/1.5
1122 BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
1123 BOTERM=2.*DZA(NK)*G*BE/1.5
1126 ENTERM=2.*REI*WTW/UPOLD
1128 CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM, &
1129 RATE,QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1),G)
1131 !...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
1132 !...IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
1134 IF(WTW.LT.1.E-3)THEN
1139 !...Calculate value of THETA-E in environment to entrain into updraft...
1141 CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
1143 !...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
1145 REI=VMFLCL*DP(NK1)*0.03/RAD ! KF (1990) Eq. 1; Kain (2004) Eq. 5
1146 TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
1148 DILBE=((TVLCL+TVQU(NK1))/(TVEN+TV0(NK1))-1.)*DZZ
1150 DILBE=((TVQU(NK)+TVQU(NK1))/(TV0(NK)+TV0(NK1))-1.)*DZZ
1152 IF(DILBE.GT.0.)ABE=ABE+DILBE*G
1154 !...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, MINIMAL
1155 !...ENTRAINMENT (0.5*REI) IS IMPOSED...
1157 IF(TVQU(NK1).LE.TV0(NK1))THEN ! Entrain/Detrain IF BLOCK
1158 EE2=0.5 ! Kain (2004) Eq. 4
1165 !...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL AIR...
1169 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
1170 QTMP=F1*Q0(NK1)+F2*QU(NK1)
1173 call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
1174 qnewlq,qnewic,XLV1,XLV0)
1175 TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
1176 IF(TU95.GT.TV0(NK1))THEN
1183 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
1184 QTMP=F1*Q0(NK1)+F2*QU(NK1)
1187 call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
1188 qnewlq,qnewic,XLV1,XLV0)
1189 TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
1190 TVDIFF = ABS(TU10-TVQU(NK1))
1191 IF(TVDIFF.LT.1.e-3)THEN
1196 EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
1197 EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
1198 EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
1199 IF(EQFRC(NK1).EQ.1)THEN
1202 ELSEIF(EQFRC(NK1).EQ.0.)THEN
1207 !...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
1208 ! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
1210 CALL PROF5(EQFRC(NK1),EE2,UD2)
1214 ENDIF ! End of Entrain/Detrain IF BLOCK
1217 !...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE FRACTIONAL
1218 ! VALUES IN THE LAYER...
1220 EE2 = AMAX1(EE2,0.5)
1222 UER(NK1)=0.5*REI*(EE1+EE2)
1223 UDR(NK1)=0.5*REI*(UD1+UD2)
1225 !...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
1226 ! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATIONS...
1228 IF(UMF(NK)-UDR(NK1).LT.10.)THEN
1230 !...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL UPD MASS
1231 ! FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE PREVIOUS MODEL LVL..
1232 ! First, correct ABE calculation if needed...
1238 ! WRITE(98,1015)P0(NK1)/100.
1243 UPOLD=UMF(NK)-UDR(NK1)
1244 UPNEW=UPOLD+UER(NK1)
1246 DILFRC(NK1) = UPNEW/UPOLD
1248 !...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND
1249 !...ICE IN THE DETRAINING UPDRAFT MASS...
1251 DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
1252 DETIC(NK1)=QICE(NK1)*UDR(NK1)
1254 QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
1255 THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
1256 QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
1257 QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
1259 !...PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF
1260 !...LIQUID PRECIP AT A GIVEN MODEL LVL, PPTICE THE SAME FOR ICE,
1261 !...TRPPT IS THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE
1262 !...CURRENT MODEL LEVEL...
1264 PPTLIQ(NK1)=QLQOUT(NK1)*UMF(NK)
1265 PPTICE(NK1)=QICOUT(NK1)*UMF(NK)
1267 TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
1268 IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
1273 !...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIUM
1274 ! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
1275 ! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE BETWEEN
1276 ! THE LET AND CLOUD TOP...
1278 !...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL VELOCITY
1279 ! FIRST BECOMES NEGATIVE...
1282 CLDHGT(LC)=Z0(LTOP)-ZLCL
1284 !...Instead of using the same minimum cloud height (for deep convection)
1285 !...everywhere, try specifying minimum cloud depth as a function of TLCL...
1289 IF(TLCL.GT.293.)THEN
1291 ELSEIF(TLCL.LE.293. .and. TLCL.GE.273)THEN
1292 CHMIN = 2.E3 + 100.*(TLCL-273.)
1293 ELSEIF(TLCL.LT.273.)THEN
1298 !...If cloud top height is less than the specified minimum for deep
1299 !...convection, save value to consider this level as source for
1300 !...shallow convection, go back up to check next level...
1302 !...Try specifying minimum cloud depth as a function of TLCL...
1305 !...DO NOT ALLOW ANY CLOUD FROM THIS LAYER IF:
1307 !... 1.) if there is no CAPE, or
1308 !... 2.) cloud top is at model level just above LCL, or
1309 !... 3.) cloud top is within updraft source layer, or
1310 !... 4.) cloud-top detrainment layer begins within
1311 !... updraft source layer.
1313 IF(LTOP.LE.KLCL .or. LTOP.LE.KPBL .or. LET+1.LE.KPBL)THEN ! No Convection Allowed
1325 ELSEIF(CLDHGT(LC).GT.CHMIN .and. ABE.GT.1)THEN ! Deep Convection allowed
1330 !...TO DISALLOW SHALLOW CONVECTION, COMMENT OUT NEXT LINE !!!!!!!!
1333 EXIT usl ! Shallow Convection from this layer
1335 ! Remember this layer (by virtue of non-zero CLDHGT) as potential shallow-cloud layer
1350 KSTART=MAX0(KPBL,KLCL)
1354 !...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS FL
1358 UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
1359 DETLQ(LTOP)=QLIQ(LTOP)*UDR(LTOP)*UPNEW/UPOLD
1360 DETIC(LTOP)=QICE(LTOP)*UDR(LTOP)*UPNEW/UPOLD
1365 ! BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
1371 DUMFDP=UMF(LET)/DPTT
1373 !...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
1374 ! RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
1378 !...entrainment is allowed at every level except for LTOP, so disallow
1379 !...entrainment at LTOP and adjust entrainment rates between LET and LTOP
1380 !...so the the dilution factor due to entyrianment is not changed but
1381 !...the actual entrainment rate will change due due forced total
1382 !...detrainment in this layer...
1387 DETLQ(NK) = UDR(NK)*QLIQ(NK)*DILFRC(NK)
1388 DETIC(NK) = UDR(NK)*QICE(NK)*DILFRC(NK)
1390 UMF(NK)=UMF(NK-1)-DP(NK)*DUMFDP
1391 UER(NK)=UMF(NK)*(1.-1./DILFRC(NK))
1392 UDR(NK)=UMF(NK-1)-UMF(NK)+UER(NK)
1393 DETLQ(NK)=UDR(NK)*QLIQ(NK)*DILFRC(NK)
1394 DETIC(NK)=UDR(NK)*QICE(NK)*DILFRC(NK)
1397 TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
1398 PPTLIQ(NK)=UMF(NK-1)*QLQOUT(NK)
1399 PPTICE(NK)=UMF(NK-1)*QICOUT(NK)
1400 TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
1405 ! Initialize some arrays below cloud base and above cloud top...
1408 IF(T0(NK).GT.T00)ML=NK
1413 UMF(NK)=VMFLCL*DP(NK)/DPTHMX
1414 UER(NK)=VMFLCL*DP(NK)/DPTHMX
1415 ELSEIF(NK.LE.KPBL)THEN
1416 UER(NK)=VMFLCL*DP(NK)/DPTHMX
1417 UMF(NK)=UMF(NK-1)+UER(NK)
1422 TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
1443 CALL ENVIRTHT(P0(NK),T0(NK),Q0(NK),THETEE(NK),ALIQ,BLIQ,CLIQ,DLIQ)
1450 !...DEFINE VARIABLES ABOVE CLOUD TOP...
1484 EMS(NK)=DP(NK)*DXSQ/G
1487 !...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION SCHEME
1489 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
1490 THTAU(NK)=TU(NK)*EXN(NK)
1491 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
1492 THTA0(NK)=T0(NK)*EXN(NK)
1493 DDILFRC(NK) = 1./DILFRC(NK)
1496 ! IF (XTIME.LT.10.)THEN
1497 ! WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
1498 ! * TMIX-T00,PMIX,QMIX,ABE
1499 ! WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
1503 !...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
1504 !...AND MIDTROPOSPHERE IS USED.
1506 WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
1507 WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
1508 WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
1509 VCONV=.5*(WSPD(KLCL)+WSPD(L5))
1510 !...for ETA model, DX is a function of location...
1511 ! TIMEC=DX(I,J)/VCONV
1514 TIMEC=AMAX1(1800.,TIMEC) ! 30 minutes >= TIMEC <= 60 minutes
1515 TIMEC=AMIN1(3600.,TIMEC)
1516 IF(ISHALL.EQ.1)TIMEC=2400. ! shallow convection TIMEC = 40 minutes
1520 !...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
1522 IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
1527 VWS=(U0(LTOP)-U0(KLCL))*(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL))* &
1529 VWS=1.E3*SHSIGN*SQRT(VWS)/(Z0(LTOP)-Z0(LCL))
1530 PEF=1.591+VWS*(-.639+VWS*(9.53E-2-VWS*4.96E-3))
1534 !...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
1536 CBH=(ZLCL-Z0(1))*3.281E-3
1540 RCBH=.96729352+CBH*(-.70034167+CBH*(.162179896+CBH*(- &
1541 1.2569798E-2+CBH*(4.2772E-4-CBH*5.44E-6))))
1543 IF(CBH.GT.25)RCBH=2.4
1545 PEFCBH=AMIN1(PEFCBH,.9)
1547 !... MEAN PEF. IS USED TO COMPUTE RAINFALL.
1549 PEFF=.5*(PEF+PEFCBH)
1550 PEFF2 = PEFF ! JSK MODS
1552 ! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1553 WRITE(message,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1554 CALL wrf_message( message )
1557 ! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1558 !*****************************************************************
1560 ! COMPUTE DOWNDRAFT PROPERTIES *
1562 !*****************************************************************
1566 devap:IF(ISHALL.EQ.1)THEN
1570 !...start downdraft about 150 mb above cloud base...
1572 ! KSTART=MAX0(KPBL,KLCL)
1573 ! KSTART=KPBL ! Changed 7/23/99
1574 KSTART=KPBL+1 ! Changed 7/23/99
1577 DPPP = P0(KSTART)-P0(NK)
1578 ! IF(DPPP.GT.200.E2)THEN
1579 IF(DPPP.GT.150.E2)THEN
1584 KLFS = MIN0(KLFS,LET-1)
1587 !...if LFS is not at least 50 mb above cloud base (implying that the
1588 !...level of equil temp, LET, is just above cloud base) do not allow a
1591 IF((P0(KSTART)-P0(LFS)).GT.50.E2)THEN
1592 THETED(LFS) = THETEE(LFS)
1595 !...call tpmix2dd to find wet-bulb temp, qv...
1597 call tpmix2dd(p0(lfs),theted(lfs),tz(lfs),qss,i,j)
1598 THTAD(LFS)=TZ(LFS)*(P00/P0(LFS))**(0.2854*(1.-0.28*QSS))
1600 !...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX...
1602 TVD(LFS)=TZ(LFS)*(1.+0.608*QSS)
1603 RDD=P0(LFS)/(R*TVD(LFS))
1608 RHBAR = RH(LFS)*DP(LFS)
1610 DO ND = LFS-1,KSTART,-1
1612 DER(ND)=DER(LFS)*EMS(ND)/EMS(LFS)
1614 DMF(ND)=DMF(ND1)+DER(ND)
1615 THETED(ND)=(THETED(ND1)*DMF(ND1)+THETEE(ND)*DER(ND))/DMF(ND)
1616 QD(ND)=(QD(ND1)*DMF(ND1)+Q0(ND)*DER(ND))/DMF(ND)
1618 RHBAR = RHBAR+RH(ND)*DP(ND)
1621 DMFFRC = 2.*(1.-RHBAR) ! Kain (2004) eq. 11
1623 !...Calculate melting effect
1624 !... first, compute total frozen precipitation generated...
1628 PPTMLT = PPTMLT+PPTICE(NK)
1631 !...For now, calculate melting effect as if DMF = -UMF at KLCL, i.e., as
1632 !...if DMFFRC=1. Otherwise, for small DMFFRC, DTMELT gets too large!
1634 DTMELT = RLF*PPTMLT/(CP*UMF(KLCL))
1638 LDT = MIN0(LFS-1,KSTART-1)
1640 call tpmix2dd(p0(kstart),theted(kstart),tz(kstart),qss,i,j)
1642 tz(kstart) = tz(kstart)-dtmelt
1643 ES=ALIQ*EXP((BLIQ*TZ(KSTART)-CLIQ)/(TZ(KSTART)-DLIQ))
1644 QSS=0.622*ES/(P0(KSTART)-ES)
1645 THETED(KSTART)=TZ(KSTART)*(1.E5/P0(KSTART))**(0.2854*(1.-0.28*QSS))* &
1646 EXP((3374.6525/TZ(KSTART)-2.5403)*QSS*(1.+0.81*QSS))
1648 LDT = MIN0(LFS-1,KSTART-1)
1651 THETED(ND) = THETED(KSTART)
1654 !...call tpmix2dd to find wet bulb temp, saturation mixing ratio...
1656 call tpmix2dd(p0(nd),theted(nd),tz(nd),qss,i,j)
1659 !...specify RH decrease of 20%/km in downdraft...
1661 RHH = 1.-0.2/1000.*(Z0(KSTART)-Z0(ND))
1663 !...adjust downdraft TEMP, Q to specified RH:
1666 DSSDT=(CLIQ-BLIQ*DLIQ)/((TZ(ND)-DLIQ)*(TZ(ND)-DLIQ))
1668 DTMP=RL*QSS*(1.-RHH)/(CP+RL*RHH*QSS*DSSDT)
1670 ES=RHH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
1671 QSRH=0.622*ES/(P0(ND)-ES)
1673 !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
1674 !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
1676 IF(QSRH.LT.QD(ND))THEN
1678 T1RH=TZ(ND)+(QSS-QSRH)*RL/CP
1684 TVD(nd) = tz(nd)*(1.+0.608*qsd(nd))
1685 IF(TVD(ND).GT.TV0(ND).OR.ND.EQ.1)THEN
1690 IF((P0(LDB)-P0(LFS)) .gt. 50.E2)THEN ! minimum Downdraft depth!
1693 DDR(ND) = -DMF(KSTART)*DP(ND)/DPDD
1695 DMF(ND) = DMF(ND1)+DDR(ND)
1696 TDER=TDER+(QSD(nd)-QD(ND))*DDR(ND)
1698 THTAD(ND)=TZ(ND)*(P00/P0(ND))**(0.2854*(1.-0.28*QD(ND)))
1704 !...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
1705 !...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
1707 d_mf: IF(TDER.LT.1.)THEN
1709 !3004 FORMAT(' ','No Downdraft!; I=',I3,2X,'J=',I3,'ISHALL =',I2)
1727 DDINC = -DMFFRC*UMF(KLCL)/DMF(KSTART)
1729 IF(TDER*DDINC.GT.TRPPT)THEN
1734 DMF(NK)=DMF(NK)*DDINC
1735 DER(NK)=DER(NK)*DDINC
1736 DDR(NK)=DDR(NK)*DDINC
1742 ! write(98,*)'PRECIP EFFICIENCY =',PEFF
1743 write(message,*)'PRECIP EFFICIENCY =',PEFF
1744 CALL wrf_message(message)
1749 !...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
1750 ! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
1751 ! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
1754 ! UMF(NK)=UMF(NK)*UPDINC
1755 ! UDR(NK)=UDR(NK)*UPDINC
1756 ! UER(NK)=UER(NK)*UPDINC
1757 ! PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
1758 ! PPTICE(NK)=PPTICE(NK)*UPDINC
1759 ! DETLQ(NK)=DETLQ(NK)*UPDINC
1760 ! DETIC(NK)=DETIC(NK)*UPDINC
1763 !...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
1793 !...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE INFLOW
1794 ! INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN IS AVAILABLE
1795 ! IN THAT LAYER INITIALLY...
1800 IF((UER(NK)-DER(NK)).GT.1.e-3)THEN
1801 AINCM1=EMS(NK)/((UER(NK)-DER(NK))*TIMEC)
1802 AINCMX=AMIN1(AINCMX,AINCM1)
1806 IF(AINCMX.LT.AINC)AINC=AINCMX
1808 !...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRAFT AND DOWNDRAFT...THEY WILL
1809 !...BE ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE STABILIZATION
1815 DETLQ2(NK)=DETLQ(NK)
1816 DETIC2(NK)=DETIC(NK)
1829 IF(ISHALL.EQ.1)THEN ! First for shallow convection
1831 ! No iteration for shallow convection; if turbulent kinetic energy (TKE) is available
1832 ! from a turbulence parameterization, scale cloud-base updraft mass flux as a function
1833 ! of TKE, but for now, just specify shallow-cloud mass flux using TKEMAX = 5...
1835 !...find the maximum TKE value between LC and KLCL...
1838 ! DO 173 K = LC,KLCL
1840 ! TKEMAX = AMAX1(TKEMAX,Q2(I,J,NK))
1842 ! TKEMAX = AMIN1(TKEMAX,10.)
1843 ! TKEMAX = AMAX1(TKEMAX,5.)
1845 !c...3_24_99...DPMIN was changed for shallow convection so that it is the
1846 !c... the same as for deep convection (5.E3). Since this doubles
1847 !c... (roughly) the value of DPTHMX, add a factor of 0.5 to calcu-
1848 !c... lation of EVAC...
1849 !c EVAC = TKEMAX*0.1
1850 EVAC = 0.5*TKEMAX*0.1
1851 ! AINC = 0.1*DPTHMX*DXIJ*DXIJ/(VMFLCL*G*TIMEC)
1852 ! AINC = EVAC*DPTHMX*DX(I,J)*DX(I,J)/(VMFLCL*G*TIMEC)
1853 AINC = EVAC*DPTHMX*DXSQ/(VMFLCL*G*TIMEC)
1857 UMF(NK)=UMF2(NK)*AINC
1858 DMF(NK)=DMF2(NK)*AINC
1859 DETLQ(NK)=DETLQ2(NK)*AINC
1860 DETIC(NK)=DETIC2(NK)*AINC
1861 UDR(NK)=UDR2(NK)*AINC
1862 UER(NK)=UER2(NK)*AINC
1863 DER(NK)=DER2(NK)*AINC
1864 DDR(NK)=DDR2(NK)*AINC
1866 ENDIF ! Otherwise for deep convection
1867 ! use iterative procedure to find mass fluxes...
1868 iter: DO NCOUNT=1,10
1870 !*****************************************************************
1872 ! COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE *
1874 !*****************************************************************
1876 !...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
1877 !...SATISFY MASS CONTINUITY...
1881 DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
1883 OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
1884 ABSOMG = ABS(OMG(NK))
1885 ABSOMGTC = ABSOMG*TIMEC
1886 FRDP = 0.75*DP(NK-1)
1887 IF(ABSOMGTC.GT.FRDP)THEN
1896 NSTEP=NINT(TIMEC/DTT+1)
1897 DTIME=TIMEC/FLOAT(NSTEP)
1898 FXM(NK)=OMG(NK)*DXSQ/G
1901 !...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
1905 !...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED ON
1906 !...SIGN OF OMEGA...
1915 IF(OMG(NK).LE.0.)THEN
1916 THFXIN(NK)=-FXM(NK)*THPA(NK-1)
1917 QFXIN(NK)=-FXM(NK)*QPA(NK-1)
1918 THFXOUT(NK-1)=THFXOUT(NK-1)+THFXIN(NK)
1919 QFXOUT(NK-1)=QFXOUT(NK-1)+QFXIN(NK)
1921 THFXOUT(NK)=FXM(NK)*THPA(NK)
1922 QFXOUT(NK)=FXM(NK)*QPA(NK)
1923 THFXIN(NK-1)=THFXIN(NK-1)+THFXOUT(NK)
1924 QFXIN(NK-1)=QFXIN(NK-1)+QFXOUT(NK)
1928 !...UPDATE THE THETA AND QV VALUES AT EACH LEVEL...
1931 THPA(NK)=THPA(NK)+(THFXIN(NK)+UDR(NK)*THTAU(NK)+DDR(NK)* &
1932 THTAD(NK)-THFXOUT(NK)-(UER(NK)-DER(NK))*THTA0(NK))* &
1934 QPA(NK)=QPA(NK)+(QFXIN(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK)- &
1935 QFXOUT(NK)-(UER(NK)-DER(NK))*Q0(NK))*DTIME*EMSD(NK)
1943 !...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO, BORROW
1944 !...MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO...
1947 IF(QG(NK).LT.0.)THEN
1948 IF(NK.EQ.1)THEN ! JSK MODS
1949 ! PRINT *,' PROBLEM WITH KF SCHEME: ' ! JSK MODS
1950 ! PRINT *,'QG = 0 AT THE SURFACE!!!!!!!' ! JSK MODS
1951 CALL wrf_error_fatal ( 'QG, QG(NK).LT.0') ! JSK MODS
1957 TMA=QG(NK1)*EMS(NK1)
1958 TMB=QG(NK-1)*EMS(NK-1)
1959 TMM=(QG(NK)-1.E-9)*EMS(NK )
1960 BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
1961 ACOEFF=BCOEFF*TMA/TMB
1965 QVDIFF=(QG(NK1)-TMA*EMSD(NK1))*100./QG(NK1)
1966 ! IF(ABS(QVDIFF).GT.1.)THEN
1967 ! PRINT *,'!!!WARNING!!! CLOUD BASE WATER VAPOR CHANGES BY ', &
1969 ! '% WHEN MOISTURE IS BORROWED TO PREVENT NEGATIVE ', &
1970 ! 'VALUES IN KAIN-FRITSCH'
1974 QG(NK1)=TMA*EMSD(NK1)
1975 QG(NK-1)=TMB*EMSD(NK-1)
1978 TOPOMG=(UDR(LTOP)-UER(LTOP))*DP(LTOP)*EMSD(LTOP)
1979 IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
1980 ! WRITE(99,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME; &
1981 ! TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1982 ! TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1988 !...CONVERT THETA TO T...
1991 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QG(NK)))
1992 TG(NK)=THTAG(NK)/EXN(NK)
1993 TVG(NK)=TG(NK)*(1.+0.608*QG(NK))
1999 !*******************************************************************
2001 ! COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY. *
2003 !*******************************************************************
2005 !...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
2011 !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
2012 !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
2016 TMIX=TMIX+DP(NK)*TG(NK)
2017 QMIX=QMIX+DP(NK)*QG(NK)
2021 ES=ALIQ*EXP((TMIX*BLIQ-CLIQ)/(TMIX-DLIQ))
2022 QSS=0.622*ES/(PMIX-ES)
2024 !...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
2028 CPM=CP*(1.+0.887*QMIX)
2029 DSSDT=QSS*(CLIQ-BLIQ*DLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
2030 DQ=(QMIX-QSS)/(1.+RL*DSSDT/CPM)
2036 EMIX=QMIX*PMIX/(0.622+QMIX)
2042 value=(indlu-1)*binc+astrt
2043 aintrp=(a1-value)/binc
2044 tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
2045 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2046 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
2047 TLCL=AMIN1(TLCL,TMIX)
2049 TVLCL=TLCL*(1.+0.608*QMIX)
2050 ZLCL = ZMIX+(TLCL-TMIX)/GDRY
2053 IF(ZLCL.LE.Z0(NK))THEN
2058 DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
2060 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
2062 TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
2063 QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
2064 TVEN=TENV*(1.+0.608*QENV)
2065 PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
2066 THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
2067 EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
2069 !...COMPUTE ADJUSTED ABE(ABEG).
2074 THETEU(NK1) = THETEU(NK)
2076 call tpmix2dd(p0(nk1),theteu(nk1),tgu(nk1),qgu(nk1),i,j)
2078 TVQU(NK1)=TGU(NK1)*(1.+0.608*QGU(NK1)-QLIQ(NK1)-QICE(NK1))
2081 DILBE=((TVLCL+TVQU(NK1))/(TVEN+TVG(NK1))-1.)*DZZ
2084 DILBE=((TVQU(NK)+TVQU(NK1))/(TVG(NK)+TVG(NK1))-1.)*DZZ
2086 IF(DILBE.GT.0.)ABEG=ABEG+DILBE*G
2088 !...DILUTE BY ENTRAINMENT BY THE RATE AS ORIGINAL UPDRAFT...
2090 CALL ENVIRTHT(P0(NK1),TG(NK1),QG(NK1),THTEEG(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
2091 THETEU(NK1)=THETEU(NK1)*DDILFRC(NK1)+THTEEG(NK1)*(1.-DDILFRC(NK1))
2094 !...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
2095 !...THE PERIOD TIMEC...
2099 ! write(98,*)'TAU, I, J, =',NTSD,I,J
2100 ! WRITE(98,1060)FABE
2104 DABE=AMAX1(ABE-ABEG,0.1*ABE)
2106 IF(FABE.GT.1. .and. ISHALL.EQ.0)THEN
2107 ! WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS
2108 ! *GRID POINT; NO CONVECTION ALLOWED!'
2112 IF(ABS(AINC-AINCOLD).LT.0.0001)THEN
2117 DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
2126 IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
2128 ! write(98,*)'TAU, I, J, =',NTSD,I,J
2129 ! WRITE(98,1055)FABE
2133 IF((FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB) .or. NCOUNT.EQ.10)THEN
2136 IF(NCOUNT.GT.10)THEN
2138 ! write(98,*)'TAU, I, J, =',NTSD,I,J
2139 ! WRITE(98,1060)FABE
2144 !...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE CONVECTIVE
2145 !...MASS FLUX BY THE FACTOR AINC:
2150 IF(DABE.LT.1.e-4)THEN
2155 AINC=AINC*STAB*ABE/DABE
2158 ! AINC=AMIN1(AINCMX,AINC)
2159 AINC=AMIN1(AINCMX,AINC)
2160 !...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION ! JSK MODS
2161 !...WILL BE MINIMAL SO JUST IGNORE IT... ! JSK MODS
2162 IF(AINC.LT.0.05)then
2165 ! AINC=AMAX1(AINC,0.05) ! JSK MODS
2168 ! IF (XTIME.LT.10.)THEN
2169 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,
2173 UMF(NK)=UMF2(NK)*AINC
2174 DMF(NK)=DMF2(NK)*AINC
2175 DETLQ(NK)=DETLQ2(NK)*AINC
2176 DETIC(NK)=DETIC2(NK)*AINC
2177 UDR(NK)=UDR2(NK)*AINC
2178 UER(NK)=UER2(NK)*AINC
2179 DER(NK)=DER2(NK)*AINC
2180 DDR(NK)=DDR2(NK)*AINC
2183 !...GO BACK UP FOR ANOTHER ITERATION...
2188 !...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
2190 !...FRC2 IS THE FRACTION OF TOTAL CONDENSATE ! PPT FB MODS
2191 !...GENERATED THAT GOES INTO PRECIPITIATION ! PPT FB MODS
2193 ! Redistribute hydormeteors according to the final mass-flux values:
2196 FRC2=PPTFLX/(CPR*AINC) ! PPT FB MODS
2205 RAINFB(NK)=PPTLIQ(NK)*AINC*FBFRC*FRC2 ! PPT FB MODS
2206 SNOWFB(NK)=PPTICE(NK)*AINC*FBFRC*FRC2 ! PPT FB MODS
2210 !...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH LAYER
2211 !...BASED ON THE SIGN OF OMEGA...
2224 IF(OMG(NK).LE.0.)THEN
2225 QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
2226 QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
2227 QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
2228 QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
2229 QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
2230 QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
2231 QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
2232 QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
2234 QLFXOUT(NK)=FXM(NK)*QLPA(NK)
2235 QIFXOUT(NK)=FXM(NK)*QIPA(NK)
2236 QRFXOUT(NK)=FXM(NK)*QRPA(NK)
2237 QSFXOUT(NK)=FXM(NK)*QSPA(NK)
2238 QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
2239 QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
2240 QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
2241 QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
2245 !...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
2248 QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))*DTIME*EMSD(NK)
2249 QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))*DTIME*EMSD(NK)
2250 QRPA(NK)=QRPA(NK)+(QRFXIN(NK)-QRFXOUT(NK)+RAINFB(NK))*DTIME*EMSD(NK) ! PPT FB MODS
2251 QSPA(NK)=QSPA(NK)+(QSFXIN(NK)-QSFXOUT(NK)+SNOWFB(NK))*DTIME*EMSD(NK) ! PPT FB MODS
2261 !...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
2264 ! IF (XTIME.LT.10.)THEN
2265 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
2268 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
2269 WRITE(message,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
2270 CALL wrf_message(message)
2274 !...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
2278 ! if(I.eq.16 .and. J.eq.41)then
2279 ! IF(ISTOP.EQ.1)THEN
2281 ! write(98,*)'At t(h), I, J =',float(NTSD)*72./3600.,I,J
2282 write(message,*)'P(LC), DTP, WKL, WKLCL =',p0(LC)/100., &
2283 TLCL+DTLCL+dtrh-TENV,WKL,WKLCL
2284 call wrf_message(message)
2285 write(message,*)'TLCL, DTLCL, DTRH, TENV =',TLCL,DTLCL, &
2287 call wrf_message(message)
2288 WRITE(message,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG, &
2289 TMIX-T00,PMIX,QMIX,ABE
2290 call wrf_message(message)
2291 WRITE(message,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100., &
2293 call wrf_message(message)
2294 WRITE(message,1035)PEF,PEFCBH,LC,LET,WKL,VWS
2295 call wrf_message(message)
2296 write(message,*)'PRECIP EFFICIENCY =',PEFF
2297 call wrf_message(message)
2298 WRITE(message,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
2299 call wrf_message(message)
2302 WRITE(message,1070)' P ',' DP ',' DT K/D ',' DR K/D ',' OMG ', &
2303 ' DOMGDP ',' UMF ',' UER ',' UDR ',' DMF ',' DER ' &
2304 ,' DDR ',' EMS ',' W0 ',' DETLQ ',' DETIC '
2305 call wrf_message(message)
2306 write(message,*)'just before DO 300...'
2307 call wrf_message(message)
2311 DTT=(TG(K)-T0(K))*86400./TIMEC
2313 DR=-(QG(K)-Q0(K))*RL*86400./(TIMEC*CP)
2314 UDFRC=UDR(K)*TIMEC*EMSD(K)
2315 UEFRC=UER(K)*TIMEC*EMSD(K)
2316 DDFRC=DDR(K)*TIMEC*EMSD(K)
2317 DEFRC=-DER(K)*TIMEC*EMSD(K)
2318 WRITE(message,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)*1.E4, &
2319 UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC,DDFRC,EMS(K)/1.E11, &
2320 W0AVG1D(K)*1.E2,DETLQ(K)*TIMEC*EMSD(K)*1.E3,DETIC(K)* &
2322 call wrf_message(message)
2324 WRITE(message,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG', &
2325 'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
2326 call wrf_message(message)
2331 IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
2333 IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
2334 IF(T0(K).LT.T00)THEN
2335 ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
2337 ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
2339 QGS=ES*0.622/(P0(K)-ES)
2342 WRITE(message,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC, &
2343 TDC,Q0(K)*1000.,QG(K)*1000.,(QG(K)-Q0(K))*1000.,QU(K)* &
2344 1000.,QD(K)*1000.,QLG(K)*1000.,QIG(K)*1000.,QRG(K)*1000., &
2345 QSG(K)*1000.,RH0,RHG
2346 call wrf_message(message)
2349 !...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
2350 !...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
2352 ! IF(ISTOP.EQ.1 .or. ISHALL.EQ.1)THEN
2354 ! IF(ISHALL.NE.1)THEN
2355 ! write(98,4421)i,j,iyr,imo,idy,ihr,imn
2356 ! write(98)i,j,iyr,imo,idy,ihr,imn,kl
2362 write(98,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000., &
2363 u0(k),v0(k),W0AVG1D(K),dp(k),tke(k)
2364 ! write(98) p0,t0,q0,u0,v0,w0,dp,tke
2365 ! WRITE(98,1115)Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000.,
2366 ! * U0(K),V0(K),DP(K)/100.,W0AVG(I,J,K)
2369 CALL wrf_error_fatal ( 'KAIN-FRITSCH, istop=1, diags' )
2374 CNDTNF=(1.-EQFRC(LFS))*(QLIQ(LFS)+QICE(LFS))*DMF(LFS)
2375 PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
2376 RAINCV(I,J)=DT*PRATEC(I,J) ! PPT FB MODS
2377 ! RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ ! PPT FB MODS
2378 ! RNC=0.1*TIMEC*PPTFLX/DXSQ
2380 IF(ISHALL.EQ.0.AND.IPRNT)write (98,909)I,J,RNC
2382 ! WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
2384 ! EVALUATE MOISTURE BUDGET...
2392 QINIT=QINIT+Q0(NK)*EMS(NK)
2393 QFNL=QFNL+QG(NK)*EMS(NK)
2394 QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
2396 QFNL=QFNL+PPTFLX*TIMEC*(1.-FBFRC) ! PPT FB MODS
2397 ! QFNL=QFNL+PPTFLX*TIMEC ! PPT FB MODS
2398 ERR2=(QFNL-QINIT)*100./QINIT
2399 IF(IPRNT)WRITE(98,1110)QINIT,QFNL,ERR2
2400 IF(ABS(ERR2).GT.0.05 .AND. ISTOP.EQ.0)THEN
2401 ! write(99,*)'!!!!!!!! MOISTURE BUDGET ERROR IN KFPARA !!!'
2402 ! WRITE(99,1110)QINIT,QFNL,ERR2
2409 ! write(99,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000., &
2410 ! u0(k),v0(k),W0AVG1D(K),dp(k)
2411 ! write(98) p0,t0,q0,u0,v0,w0,dp,tke
2412 ! WRITE(98,1115)P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
2413 ! U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
2414 WRITE(98,4456)P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
2415 U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
2422 1115 FORMAT (2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
2424 IF(PPTFLX.GT.0.)THEN
2425 RELERR=ERR2*QINIT/(PPTFLX*TIMEC)
2430 WRITE(98,1120)RELERR
2431 WRITE(98,*)'TDER, CPR, TRPPT =', &
2432 TDER,CPR*AINC,TRPPT*AINC
2435 !...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
2437 !...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
2438 !...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
2440 IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
2441 NCA(I,J) = REAL(NIC)*DT
2449 ! IF(IMOIST(INEST).NE.2)THEN
2451 !...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMATED
2452 !...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
2453 !...NOTE: THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
2454 !...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
2457 ! RLC=XLV0-XLV1*TG(K)
2458 ! RLS=XLS0-XLS1*TG(K)
2459 ! CPM=CP*(1.+0.887*QG(K))
2460 ! TG(K)=TG(K)-(RLC*(QLG(K)+QRG(K))+RLS*(QIG(K)+QSG(K)))/CPM
2461 ! QG(K)=QG(K)+(QLG(K)+QRG(K)+QIG(K)+QSG(K))
2468 !...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
2470 IF(.NOT. F_QI .and. warm_rain)THEN
2472 CPM=CP*(1.+0.887*QG(K))
2473 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
2474 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
2476 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
2478 ELSEIF(.NOT. F_QI .and. .not. warm_rain)THEN
2480 !...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROMETEORS
2481 !...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
2483 CPM=CP*(1.+0.887*QG(K))
2485 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
2487 TG(K)=TG(K)+(QLG(K)+QRG(K))*RLF/CPM
2489 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
2491 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
2495 !...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE TENDENCIES
2496 !...OF HYDROMETEORS DIRECTLY...
2498 DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
2499 DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
2500 DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
2502 DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
2504 DQIDT(K)=DQIDT(K)+(QSG(K)-QS0(K))/TIMEC
2507 ! PRINT *,'THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED!'
2508 CALL wrf_error_fatal ( 'KAIN-FRITSCH, THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED' )
2510 DTDT(K)=(TG(K)-T0(K))/TIMEC
2511 DQDT(K)=(QG(K)-Q0(K))/TIMEC
2513 PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
2514 RAINCV(I,J)=DT*PRATEC(I,J)
2515 ! RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ ! PPT FB MODS
2516 ! RNC=0.1*TIMEC*PPTFLX/DXSQ
2518 909 FORMAT('AT I, J =',i3,1x,i3,' CONVECTIVE RAINFALL =',F8.4,' mm')
2519 ! write (98,909)I,J,RNC
2520 ! write (6,909)I,J,RNC
2521 ! WRITE(98,*)'at NTSD =',NTSD,',No. of KF points activated =',
2524 1000 FORMAT(' ',10A8)
2525 1005 FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
2526 1010 FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
2527 1015 FORMAT(' ','ALL REMAINING MASS DETRAINS BELOW ',F4.0,' MB')
2528 1025 FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M', &
2529 ' DTLCL=',F5.2,' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=', &
2530 I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1, &
2532 1030 FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =', &
2533 E12.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =', &
2535 1035 FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'WKL=' &
2537 !1055 FORMAT('*** DEGREE OF STABILIZATION =',F5.3, &
2538 ! ', NO MORE MASS FLUX IS ALLOWED!')
2539 !1060 FORMAT(' ITERATION DOES NOT CONVERGE TO GIVE THE SPECIFIED &
2540 ! &DEGREE OF STABILIZATION! FABE= ',F6.4)
2542 1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3)
2543 1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, TADVEC, NSTEP=', &
2544 2(1X,F5.0),I3,'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
2545 1085 FORMAT (A3,16A7,2A8)
2546 1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3)
2547 1095 FORMAT(' ',' PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ',F10.0)
2548 1105 FORMAT(' ','NET LATENT HEAT RELEASE =',E12.5,' ACTUAL HEATING =',&
2549 E12.5,' J/KG-S, DIFFERENCE = ',F9.3,'%')
2550 1110 FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5, &
2551 ' TOTAL WATER CHANGE =',F8.2,'%')
2552 ! 1115 FORMAT (2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
2553 1120 FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3,'%')
2555 !-----------------------------------------------------------------------
2556 !--------------SAVE CLOUD TOP AND BOTTOM FOR RADIATION------------------
2557 !-----------------------------------------------------------------------
2559 CUTOP(I,J)=REAL(LTOP)
2560 CUBOT(I,J)=REAL(LCL)
2562 !-----------------------------------------------------------------------
2563 END SUBROUTINE KF_eta_PARA
2564 !********************************************************************
2565 ! ***********************************************************************
2566 SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0)
2568 ! Lookup table variables:
2569 ! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
2570 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
2571 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
2572 ! REAL, SAVE, DIMENSION(1:200) :: ALU
2573 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
2574 ! End of Lookup table variables:
2575 !-----------------------------------------------------------------------
2577 !-----------------------------------------------------------------------
2578 REAL, INTENT(IN ) :: P,THES,XLV1,XLV0
2579 REAL, INTENT(OUT ) :: QNEWLQ,QNEWIC
2580 REAL, INTENT(INOUT) :: TU,QU,QLIQ,QICE
2581 REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11, &
2582 TEMP,QS,QNEW,DQ,QTOT,RLL,CPP
2583 INTEGER :: IPTB,ITHTB
2584 !-----------------------------------------------------------------------
2586 !c******** LOOKUP TABLE VARIABLES... ****************************
2587 ! parameter(kfnt=250,kfnp=220)
2589 ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),
2590 ! * alu(200),rdpr,rdthk,plutop
2591 !C***************************************************************
2593 !c***********************************************************************
2594 !c scaling pressure and tt table index
2595 !c***********************************************************************
2602 !***********************************************************************
2603 ! base and scaling factor for the
2604 !***********************************************************************
2606 ! scaling the and tt table index
2607 bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
2608 tth=(thes-bth)*rdthk
2611 IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
2612 write(98,*)'**** OUT OF BOUNDS *********'
2616 t00=ttab(ithtb ,iptb )
2617 t10=ttab(ithtb+1,iptb )
2618 t01=ttab(ithtb ,iptb+1)
2619 t11=ttab(ithtb+1,iptb+1)
2621 q00=qstab(ithtb ,iptb )
2622 q10=qstab(ithtb+1,iptb )
2623 q01=qstab(ithtb ,iptb+1)
2624 q11=qstab(ithtb+1,iptb+1)
2626 !***********************************************************************
2627 ! parcel temperature
2628 !***********************************************************************
2630 temp=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
2632 qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
2640 ! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
2641 ! ADJUSTED...IF LIQUID WATER IS PRESENT, IT IS ALLOWED TO EVAPORATE
2646 ! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS AT ITS
2647 ! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MIXING
2648 ! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURATION
2649 ! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPRIATE
2650 ! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE MADE.
2652 !...subsaturated values only occur in calculations involving various mixtures of
2653 !...updraft and environmental air for estimation of entrainment and detrainment.
2654 !...For these purposes, assume that reasonable estimates can be given using
2655 !...liquid water saturation calculations only - i.e., ignore the effect of the
2656 !...ice phase in this process only...will not affect conservative properties...
2659 qliq=qliq-dq*qliq/(qtot+1.e-10)
2660 qice=qice-dq*qice/(qtot+1.e-10)
2664 CPP=1004.5*(1.+0.89*QU)
2665 IF(QTOT.LT.1.E-10)THEN
2667 !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
2668 TEMP=TEMP+RLL*(DQ/(1.+DQ))/CPP
2671 !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURATION,
2672 ! THE TEMPERATURE IS GIVEN BY:
2674 TEMP=TEMP+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CPP
2686 END SUBROUTINE TPMIX2
2687 !******************************************************************************
2688 SUBROUTINE DTFRZNEW(TU,P,THTEU,QU,QFRZ,QICE,ALIQ,BLIQ,CLIQ,DLIQ)
2689 !-----------------------------------------------------------------------
2691 !-----------------------------------------------------------------------
2692 REAL, INTENT(IN ) :: P,QFRZ,ALIQ,BLIQ,CLIQ,DLIQ
2693 REAL, INTENT(INOUT) :: TU,THTEU,QU,QICE
2694 REAL :: RLC,RLS,RLF,CPP,A,DTFRZ,ES,QS,DQEVAP,PII
2695 !-----------------------------------------------------------------------
2697 !...ALLOW THE FREEZING OF LIQUID WATER IN THE UPDRAFT TO PROCEED AS AN
2698 !...APPROXIMATELY LINEAR FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE
2699 !...TTFRZ TO TBFRZ...
2700 !...FOR COLDER TEMPERATURES, FREEZE ALL LIQUID WATER...
2701 !...THERMODYNAMIC PROPERTIES ARE STILL CALCULATED WITH RESPECT TO LIQUID WATER
2702 !...TO ALLOW THE USE OF LOOKUP TABLE TO EXTRACT TMP FROM THETAE...
2704 RLC=2.5E6-2369.276*(TU-273.16)
2705 RLS=2833922.-259.532*(TU-273.16)
2707 CPP=1004.5*(1.+0.89*QU)
2709 ! A = D(es)/DT IS THAT CALCULATED FROM BUCK (1981) EMPERICAL FORMULAS
2710 ! FOR SATURATION VAPOR PRESSURE...
2712 A=(CLIQ-BLIQ*DLIQ)/((TU-DLIQ)*(TU-DLIQ))
2713 DTFRZ = RLF*QFRZ/(CPP+RLS*QU*A)
2716 ES = ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2717 QS = ES*0.622/(P-ES)
2719 !...FREEZING WARMS THE AIR AND IT BECOMES UNSATURATED...ASSUME THAT SOME OF THE
2720 !...LIQUID WATER THAT IS AVAILABLE FOR FREEZING EVAPORATES TO MAINTAIN SATURA-
2721 !...TION...SINCE THIS WATER HAS ALREADY BEEN TRANSFERRED TO THE ICE CATEGORY,
2722 !...SUBTRACT IT FROM ICE CONCENTRATION, THEN SET UPDRAFT MIXING RATIO AT THE NEW
2723 !...TEMPERATURE TO THE SATURATION VALUE...
2728 PII=(1.E5/P)**(0.2854*(1.-0.28*QU))
2729 THTEU=TU*PII*EXP((3374.6525/TU-2.5403)*QU*(1.+0.81*QU))
2731 END SUBROUTINE DTFRZNEW
2732 ! --------------------------------------------------------------------------------
2734 SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
2735 QNEWIC,QLQOUT,QICOUT,G)
2737 !-----------------------------------------------------------------------
2739 !-----------------------------------------------------------------------
2740 ! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
2741 ! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
2742 ! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
2743 ! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
2744 ! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
2746 REAL, INTENT(IN ) :: G
2747 REAL, INTENT(IN ) :: DZ,BOTERM,ENTERM,RATE
2748 REAL, INTENT(INOUT) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
2749 REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
2754 ! ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY
2755 ! BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
2758 QEST=0.5*(QTOT+QNEW)
2759 G1=WTW+BOTERM-ENTERM-2.*G*DZ*QEST/1.5
2761 WAVG=0.5*(SQRT(WTW)+SQRT(G1))
2762 CONV=RATE*DZ/WAVG ! KF90 Eq. 9
2764 ! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
2765 ! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
2766 ! IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
2767 ! SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...
2769 RATIO3=QNEWLQ/(QNEW+1.E-8)
2773 RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-8)
2774 QTOT=QTOT*EXP(-CONV) ! KF90 Eq. 9
2776 ! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
2777 ! PARCEL AT THIS LEVEL...
2781 QICOUT=(1.-RATIO4)*DQ
2783 ! ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
2784 ! LATE VERTICAL VELOCITY
2786 PPTDRG=0.5*(OLDQ+QTOT-0.2*QNEW)
2787 WTW=WTW+BOTERM-ENTERM-2.*G*DZ*PPTDRG/1.5
2788 IF(ABS(WTW).LT.1.E-4)WTW=1.E-4
2790 ! DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
2791 ! DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...
2793 QLIQ=RATIO4*QTOT+RATIO3*0.4*QNEW
2794 QICE=(1.-RATIO4)*QTOT+(1.-RATIO3)*0.4*QNEW
2798 END SUBROUTINE CONDLOAD
2800 ! ----------------------------------------------------------------------
2801 SUBROUTINE PROF5(EQ,EE,UD)
2803 !***********************************************************************
2804 !***** GAUSSIAN TYPE MIXING PROFILE....******************************
2805 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2806 ! THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN
2807 ! DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN FROM
2808 ! "HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICS TABLES"
2809 ! ED. BY ABRAMOWITZ AND STEGUN, NATL BUREAU OF STANDARDS APPLIED
2810 ! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
2813 ! Solves for KF90 Eq. 2
2815 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2816 !-----------------------------------------------------------------------
2818 !-----------------------------------------------------------------------
2819 REAL, INTENT(IN ) :: EQ
2820 REAL, INTENT(INOUT) :: EE,UD
2821 REAL :: SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
2823 DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676, &
2824 0.9372980,0.33267,0.166666667,0.202765151/
2831 C1=A1*T1+A2*T1*T1+A3*T1*T1*T1
2832 C2=A1*T2+A2*T2*T2+A3*T2*T2*T2
2834 EE=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
2835 UD=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*(0.5+EQ*EQ/2.- &
2838 EE=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
2839 UD=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*(0.5+EQ* &
2845 END SUBROUTINE PROF5
2847 ! ------------------------------------------------------------------------
2848 SUBROUTINE TPMIX2DD(p,thes,ts,qs,i,j)
2850 ! Lookup table variables:
2851 ! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
2852 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
2853 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
2854 ! REAL, SAVE, DIMENSION(1:200) :: ALU
2855 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
2856 ! End of Lookup table variables:
2857 !-----------------------------------------------------------------------
2859 !-----------------------------------------------------------------------
2860 REAL, INTENT(IN ) :: P,THES
2861 REAL, INTENT(INOUT) :: TS,QS
2862 INTEGER, INTENT(IN ) :: i,j ! avail for debugging
2863 REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11
2864 INTEGER :: IPTB,ITHTB
2865 CHARACTER*256 :: MESS
2866 !-----------------------------------------------------------------------
2869 !******** LOOKUP TABLE VARIABLES (F77 format)... ****************************
2870 ! parameter(kfnt=250,kfnp=220)
2872 ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp), &
2873 ! alu(200),rdpr,rdthk,plutop
2874 !***************************************************************
2876 !***********************************************************************
2877 ! scaling pressure and tt table index
2878 !***********************************************************************
2884 !***********************************************************************
2885 ! base and scaling factor for the
2886 !***********************************************************************
2888 ! scaling the and tt table index
2889 bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
2890 tth=(thes-bth)*rdthk
2894 t00=ttab(ithtb ,iptb )
2895 t10=ttab(ithtb+1,iptb )
2896 t01=ttab(ithtb ,iptb+1)
2897 t11=ttab(ithtb+1,iptb+1)
2899 q00=qstab(ithtb ,iptb )
2900 q10=qstab(ithtb+1,iptb )
2901 q01=qstab(ithtb ,iptb+1)
2902 q11=qstab(ithtb+1,iptb+1)
2904 !***********************************************************************
2905 ! parcel temperature and saturation mixing ratio
2906 !***********************************************************************
2908 ts=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
2910 qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
2912 END SUBROUTINE TPMIX2DD
2914 ! -----------------------------------------------------------------------
2915 SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,ALIQ,BLIQ,CLIQ,DLIQ)
2917 !-----------------------------------------------------------------------
2919 !-----------------------------------------------------------------------
2920 REAL, INTENT(IN ) :: P1,T1,Q1,ALIQ,BLIQ,CLIQ,DLIQ
2921 REAL, INTENT(INOUT) :: THT1
2922 REAL :: EE,TLOG,ASTRT,AINC,A1,TP,VALUE,AINTRP,TDPT,TSAT,THT, &
2923 T00,P00,C1,C2,C3,C4,C5
2925 !-----------------------------------------------------------------------
2926 DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834, &
2929 ! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
2931 ! NOTE: Calculations for mixed/ice phase no longer used...jsk 8/00
2932 ! For example, KF90 Eq. 10 no longer used
2935 ! TLOG=ALOG(EE/ALIQ)
2936 ! ...calculate LOG term using lookup table...
2943 value=(indlu-1)*ainc+astrt
2944 aintrp=(a1-value)/ainc
2945 tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
2947 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2948 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
2949 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
2950 THT1=THT*EXP((C1/TSAT-C2)*Q1*(1.+0.81*Q1))
2952 END SUBROUTINE ENVIRTHT
2953 ! ***********************************************************************
2954 !====================================================================
2955 SUBROUTINE kf_eta_init(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, &
2956 RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS, &
2957 SVP1,SVP2,SVP3,SVPT0, &
2958 P_FIRST_SCALAR,restart,allowed_to_read, &
2959 ids, ide, jds, jde, kds, kde, &
2960 ims, ime, jms, jme, kms, kme, &
2961 its, ite, jts, jte, kts, kte )
2962 !--------------------------------------------------------------------
2964 !--------------------------------------------------------------------
2965 LOGICAL , INTENT(IN) :: restart,allowed_to_read
2966 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
2967 ims, ime, jms, jme, kms, kme, &
2968 its, ite, jts, jte, kts, kte
2969 INTEGER , INTENT(IN) :: P_QI,P_QS,P_FIRST_SCALAR
2971 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
2979 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG
2981 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA
2983 INTEGER :: i, j, k, itf, jtf, ktf
2984 REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
2990 IF(.not.restart)THEN
3003 IF (P_QI .ge. P_FIRST_SCALAR) THEN
3013 IF (P_QS .ge. P_FIRST_SCALAR) THEN
3039 CALL KF_LUTAB(SVP1,SVP2,SVP3,SVPT0)
3041 END SUBROUTINE kf_eta_init
3043 !-------------------------------------------------------
3045 subroutine kf_lutab(SVP1,SVP2,SVP3,SVPT0)
3047 ! This subroutine is a lookup table.
3048 ! Given a series of series of saturation equivalent potential
3049 ! temperatures, the temperature is calculated.
3051 !--------------------------------------------------------------------
3053 !--------------------------------------------------------------------
3054 ! Lookup table variables
3055 ! INTEGER, SAVE, PARAMETER :: KFNT=250,KFNP=220
3056 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
3057 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
3058 ! REAL, SAVE, DIMENSION(1:200) :: ALU
3059 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
3060 ! End of Lookup table variables
3062 INTEGER :: KP,IT,ITCNT,I
3063 REAL :: DTH,TMIN,TOLER,PBOT,DPR, &
3064 TEMP,P,ES,QS,PI,THES,TGUES,THGUES,F0,T1,T0,THGS,F1,DT, &
3066 ! REAL :: ALIQ,BLIQ,CLIQ,DLIQ,SVP1,SVP2,SVP3,SVPT0
3067 REAL :: ALIQ,BLIQ,CLIQ,DLIQ
3068 REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
3070 ! equivalent potential temperature increment
3072 ! minimum starting temp
3074 ! tolerance for accuracy of temperature
3076 ! top pressure (pascals)
3078 ! bottom pressure (pascals)
3087 ! compute parameters
3089 ! 1._over_(sat. equiv. theta increment)
3091 ! pressure increment
3093 DPR=(PBOT-PLUTOP)/REAL(KFNP-1)
3094 ! dpr=(pbot-plutop)/REAL(kfnp-1)
3095 ! 1._over_(pressure increment)
3097 ! compute the spread of thes
3098 ! thespd=dth*(kfnt-1)
3100 ! calculate the starting sat. equiv. theta
3106 es=aliq*exp((bliq*temp-cliq)/(temp-dliq))
3108 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
3109 the0k(kp)=temp*pi*exp((3374.6525/temp-2.5403)*qs* &
3113 ! compute temperatures for each sat. equiv. potential temp.
3120 ! define sat. equiv. pot. temp.
3122 ! iterate to find temperature
3123 ! find initial guess
3129 es=aliq*exp((bliq*tgues-cliq)/(tgues-dliq))
3131 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
3132 thgues=tgues*pi*exp((3374.6525/tgues-2.5403)*qs* &
3140 es=aliq*exp((bliq*t1-cliq)/(t1-dliq))
3142 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
3143 thtgs=t1*pi*exp((3374.6525/t1-2.5403)*qs*(1.+0.81*qs))
3145 if(abs(f1).lt.toler)then
3149 dt=f1*(t1-t0)/(f1-f0)
3159 ! lookup table for tlog(emix/aliq)
3161 ! set up intial values for lookup tables
3172 END SUBROUTINE KF_LUTAB
3174 END MODULE module_cu_kfeta