1 !WRF:MODEL_LAYER:PHYSICS
8 REAL , PARAMETER :: RAD = 1500.
12 !-------------------------------------------------------------
14 ids,ide, jds,jde, kds,kde &
15 ,ims,ime, jms,jme, kms,kme &
16 ,its,ite, jts,jte, kts,kte &
17 ,DT,KTAU,DX,CUDT,CURR_SECS,ADAPT_STEP_FLAG &
21 ,U,V,TH,T,W,QV,dz8w,Pcps,pi &
22 ,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
23 ,EP2,SVP1,SVP2,SVP3,SVPT0 &
24 ,STEPCU,CU_ACT_FLAG,warm_rain &
26 ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
27 ,RQVCUTEN,RQCCUTEN,RQRCUTEN,RQICUTEN,RQSCUTEN &
30 !-------------------------------------------------------------
32 !-------------------------------------------------------------
33 INTEGER, INTENT(IN ) :: &
34 ids,ide, jds,jde, kds,kde, &
35 ims,ime, jms,jme, kms,kme, &
36 its,ite, jts,jte, kts,kte
38 INTEGER, INTENT(IN ) :: STEPCU
39 LOGICAL, INTENT(IN ) :: warm_rain
41 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1
42 REAL, INTENT(IN ) :: CP,R,G,EP1,EP2
43 REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
45 INTEGER, INTENT(IN ) :: KTAU
47 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
60 REAL, INTENT(IN ) :: DT, DX
61 REAL, INTENT(IN ) :: CUDT
62 REAL, INTENT(IN ) :: CURR_SECS
63 LOGICAL,OPTIONAL, INTENT(IN ) :: ADAPT_STEP_FLAG
64 REAL, INTENT (INOUT) :: CUDTACTTIME
66 REAL, DIMENSION( ims:ime , jms:jme ), &
72 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
76 LOGICAL, DIMENSION( ims:ime , jms:jme ), &
77 INTENT(INOUT) :: CU_ACT_FLAG
83 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
94 ! Flags relating to the optional tendency arrays declared above
95 ! Models that carry the optional tendencies will provdide the
96 ! optional arguments at compile time; these flags all the model
97 ! to determine at run-time whether a particular tracer is in
101 LOGICAL, OPTIONAL :: &
112 REAL, DIMENSION( kts:kte ) :: &
122 REAL, DIMENSION( kts:kte ):: &
130 REAL :: TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp
132 INTEGER :: i,j,k,NTST,ICLDCK
134 LOGICAL :: qi_flag , qs_flag
135 ! adjustable time step changes
136 REAL :: lastdt = -1.0
137 REAL :: W0AVGfctr, W0fctr, W0den
138 LOGICAL :: run_param , doing_adapt_dt , decided
140 !----------------------------------------------------------------------
142 !--- CALL CUMULUS PARAMETERIZATION
144 !...TST IS THE NUMBER OF TIME STEPS IN 10 MINUTES...W0AVG IS CLOSE TO A
145 !...RUNNING MEAN VERTICAL VELOCITY...NOTE THAT IF YOU CHANGE TST, IT WIL
146 !...CHANGE THE FREQUENCY OF THE CONVECTIVE INTITIATION CHECK (SEE BELOW)
147 !...NOTE THAT THE ORDERING OF VERTICAL LAYERS MUST BE REVERSED FOR W0AVG
148 !...BECAUSE THE ORDERING IS REVERSED IN KFPARA...
153 IF ( PRESENT( F_QI ) ) qi_flag = f_qi
154 IF ( PRESENT( F_QS ) ) qs_flag = f_qs
156 !----------------------
159 !----------------------
160 ! NTST=NINT(1200./(DT*2.))
164 !----------------------
165 ! ICLDCK=MOD(KTAU,NTST)
166 !----------------------
167 ! write(0,*) 'DT = ',DT,' KTAU = ',KTAU,' DX = ',DX
168 ! write(0,*) 'CUDT = ',CUDT,' CURR_SECS = ',CURR_SECS
169 ! write(0,*) 'ADAPT_STEP_FLAG = ',ADAPT_STEP_FLAG,' IDS = ',IDS
170 ! write(0,*) 'STEPCU = ',STEPCU,' warm_rain = ',warm_rain
171 ! write(0,*) 'F_QV = ',F_QV,' F_QC = ',F_QV
172 ! write(0,*) 'F_QI = ',F_QI,' F_QS = ',F_QS
173 ! write(0,*) 'F_QR = ',F_QR
179 if (ADAPT_STEP_FLAG) then
180 W0AVGfctr = 2 * MAX(CUDT*60,dt) - dt
182 W0den = 2 * MAX(CUDT*60,dt)
192 ! SCR1=-5.0E-4*G*rho(I,K,J)*(w(I,K,J)+w(I,K+1,J))
193 ! TV=T(I,K,J)*(1.+EP1*QV(I,K,J))
194 ! RHOE=Pcps(I,K,J)/(R*TV)
195 ! W0=-101.9368*SCR1/RHOE
196 W0=0.5*(w(I,K,J)+w(I,K+1,J))
200 ! W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST
201 ! New, to support adaptive time step:
203 W0AVG(I,K,J) = ( W0AVG(I,K,J) * W0AVGfctr + W0 * W0fctr ) / W0den
210 ! Initialization for adaptive time step.
212 doing_adapt_dt = .FALSE.
213 IF ( PRESENT(adapt_step_flag) ) THEN
214 IF ( adapt_step_flag ) THEN
215 doing_adapt_dt = .TRUE.
216 IF ( cudtacttime .EQ. 0. ) THEN
217 cudtacttime = curr_secs + cudt*60.
222 ! Do we run through this scheme or not?
224 ! Test 1: If this is the initial model time, then yes.
226 ! Test 2: If the user asked for the cumulus to be run every time step, then yes.
228 ! Test 3: If not adaptive dt, and this is on the requested cumulus frequency, then yes.
230 ! Test 4: If using adaptive dt and the current time is past the last requested activate cumulus time, then yes.
231 ! CURR_SECS >= CUDTACTTIME
233 ! If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag
234 ! to TRUE. The decided flag says that one of these tests was able to say "yes", run the scheme.
235 ! We only proceed to other tests if the previous tests all have left decided as FALSE.
237 ! If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next
242 IF ( ( .NOT. decided ) .AND. &
243 ( ktau .EQ. 1 ) ) THEN
248 IF ( ( .NOT. decided ) .AND. &
249 ( ( cudt .EQ. 0. ) .OR. ( stepcu .EQ. 1 ) ) ) THEN
254 IF ( ( .NOT. decided ) .AND. &
255 ( .NOT. doing_adapt_dt ) .AND. &
256 ( MOD(ktau,ntst) .EQ. 0 ) ) THEN
261 IF ( ( .NOT. decided ) .AND. &
262 ( doing_adapt_dt ) .AND. &
263 ( curr_secs .GE. cudtacttime ) ) THEN
266 cudtacttime = curr_secs + cudt*60
272 CU_ACT_FLAG(i,j) = .true.
278 ! if (i.eq. 110 .and. j .eq. 59 ) then
279 ! write(0,*) 'nca = ',nca(i,j),' CU_ACT_FLAG = ',CU_ACT_FLAG(i,j)
280 ! write(0,*) 'dt = ',dt,' ADAPT_STEP_FLAG = ',ADAPT_STEP_FLAG
282 ! IF ( NINT(NCA(I,J)) .gt. 0 ) then
283 IF ( NCA(I,J) .gt. 0.5*DT ) then
284 CU_ACT_FLAG(i,j) = .false.
298 ! assign vars from 3D to 1D
307 W0AVG1D(K) =W0AVG(I,K,J)
313 U1D,V1D,T1D,QV1D,P1D,DZ1D, &
314 W0AVG1D,DT,DX,DXSQ,RHO1D, &
315 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
316 EP2,SVP1,SVP2,SVP3,SVPT0, &
317 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
319 warm_rain,qi_flag,qs_flag, &
320 ids,ide, jds,jde, kds,kde, &
321 ims,ime, jms,jme, kms,kme, &
322 its,ite, jts,jte, kts,kte )
324 IF ( PRESENT( RTHCUTEN ) .AND. PRESENT( RQVCUTEN ) ) THEN
326 RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
327 RQVCUTEN(I,K,J)=DQDT(K)
331 IF( PRESENT(RQRCUTEN) .AND. PRESENT(RQCCUTEN) .AND. &
335 RQRCUTEN(I,K,J)=DQRDT(K)
336 RQCCUTEN(I,K,J)=DQCDT(K)
339 ! This is the case for Eta microphysics without 3d rain field
342 RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
347 !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
349 IF( PRESENT( RQICUTEN ) .AND. qi_flag )THEN
351 RQICUTEN(I,K,J)=DQIDT(K)
355 IF( PRESENT ( RQSCUTEN ) .AND. qs_flag )THEN
357 RQSCUTEN(I,K,J)=DQSDT(K)
369 !-----------------------------------------------------------
370 SUBROUTINE KFPARA (I, J, &
371 U0,V0,T0,QV0,P0,DZQ,W0AVG1D, &
373 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
374 EP2,SVP1,SVP2,SVP3,SVPT0, &
375 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
377 warm_rain,qi_flag,qs_flag, &
378 ids,ide, jds,jde, kds,kde, &
379 ims,ime, jms,jme, kms,kme, &
380 its,ite, jts,jte, kts,kte )
381 !-----------------------------------------------------------
383 !-----------------------------------------------------------
384 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
385 ims,ime, jms,jme, kms,kme, &
386 its,ite, jts,jte, kts,kte, &
388 LOGICAL, INTENT(IN ) :: warm_rain
389 LOGICAL :: qi_flag, qs_flag
392 REAL, DIMENSION( kts:kte ), &
402 REAL, INTENT(IN ) :: DT,DX,DXSQ
405 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
406 REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
408 REAL, DIMENSION( kts:kte ), INTENT(INOUT) :: &
416 REAL, DIMENSION( ims:ime , jms:jme ), &
417 INTENT(INOUT) :: RAINCV, &
421 !...DEFINE LOCAL VARIABLES...
423 REAL, DIMENSION( kts:kte ) :: &
424 Q0,Z0,TV0,TU,TVU,QU,TZ,TVD, &
425 QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD, &
426 UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2, &
427 UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE, &
428 THTAU,THETEU,THTAD,THETED,QLIQ,QICE, &
429 QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC, &
430 DETLQ2,DETIC2,RATIO,RATIO2
432 REAL, DIMENSION( kts:kte ) :: &
433 DOMGDP,EXN,RHOE,TVQU,DP,RH,EQFRC,WSPD, &
434 QDT,FXM,THTAG,THTESG,THPA,THFXTOP, &
435 THFXBOT,QPA,QFXTOP,QFXBOT,QLPA,QLFXIN, &
436 QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA, &
437 QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT, &
438 QL0,QLG,QI0,QIG,QR0,QRG,QS0,QSG
440 REAL, DIMENSION( kts:kte+1 ) :: OMG
441 REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
445 REAL :: P00,T00,CV,B61,RLF,RHIC,RHBC,PIE, &
447 REAL :: GDRY,ROCP,ALIQ,BLIQ, &
448 CLIQ,DLIQ,AICE,BICE,CICE,DICE
449 REAL :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX, &
450 ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL, &
451 CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR, &
452 ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
453 TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD, &
454 UPNEW,ABE,WKLCL,THTUDL,TUDL,TTEMP,FRC1, &
455 QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
456 DZZ,WSQ,UDLBE,REI,EE2,UD2,TTMP,F1,F2, &
457 THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1, &
458 UD1,CLDHGT,DPTT,QNEWLQ,DUMFDP,EE,TSAT, &
459 THTA,P150,USR,VCONV,TIMEC,SHSIGN,VWS,PEF, &
460 CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN, &
461 DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1, &
462 DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF, &
463 UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF, &
464 DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
465 AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1, &
466 DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF, &
467 TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR, &
468 UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2, &
469 RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
474 INTEGER :: ISTOP,ML,L5,L4,KMIX,LOW, &
475 LC,MXLAYR,LLFC,NLAYRS,NK, &
476 KPBL,KLCL,LCL,LET,IFLAG, &
477 KFRZ,NK1,LTOP,NJ,LTOP1, &
478 LTOPM1,LVF,KSTART,KMIN,LFS, &
479 ND,NIC,LDB,LDT,ND1,NDK, &
480 NM,LMAX,NCOUNT,NOITR, &
483 DATA P00,T00/1.E5,273.16/
484 DATA CV,B61,RLF/717.,0.608,3.339E5/
485 DATA RHIC,RHBC/1.,0.90/
486 DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
488 !-----------------------------------------------------------
508 !...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER
509 !...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL)
510 !...FIELD. 'FBFRC' IS THE FRACTION OF AVAILABLE
511 !...PRECIPITATION TO BE FED BACK (0.0 - 1.0)...
515 !...SCHEME IS CALLED ONCE ON EACH NORTH-SOUTH SLICE, THE LOOP BELOW
516 !...CHECKS FOR THE POSSIBILITY OF INITIATING PARAMETERIZED
517 !...CONVECTION AT EACH POINT WITHIN THE SLICE
519 !...SEE IF IT IS NECESSARY TO CHECK FOR CONVECTIVE TRIGGERING AT THIS
520 !...GRID POINT. IF NCA>0, CONVECTION IS ALREADY ACTIVE AT THIS POINT,
521 !...JUST FEED BACK THE TENDENCIES SAVED FROM THE TIME WHEN CONVECTION
522 !...WAS INITIATED. IF NCA<0, CONVECTION IS NOT ACTIVE
523 !...AND YOU MAY WANT TO CHECK TO SEE IF IT CAN BE ACTIVATED FOR THE
524 !...CURRENT CONDITIONS. IN PREVIOUS APLICATIONS OF THIS SCHEME,
525 !...THE VARIABLE ICLDCK WAS USED BELOW TO SAVE TIME BY ONLY CHECKING
526 !...FOR THE POSSIBILITY OF CONVECTIVE INITIATION EVERY 5 OR 10
531 !SUE P300=1000.*(PSB(I,J)*A(KL)+PTOP-30.)+PP3D(I,J,KL)
535 !...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
536 !...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
537 !...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
539 !...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
540 !...FROM BOTTOM-UP IN THE KF SCHEME...
543 !SUE tmprpsb=1./PSB(I,J)
544 !SUE CELL=PTOP*tmprpsb
547 !SUE P0(K)=1.E3*(A(NK)*PSB(I,J)+PTOP)+PP3D(I,J,NK)
549 !...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
551 ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
552 QES(K)=EP2*ES/(P0(K)-ES)
553 Q0(K)=AMIN1(QES(K),QV0(K))
554 Q0(K)=AMAX1(0.000001,Q0(K))
560 TV0(K)=T0(K)*(1.+B61*Q0(K))
561 RHOE(K)=P0(K)/(R*TV0(K))
563 DP(K)=rho(k)*g*DZQ(k)
565 !...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
566 ! DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
568 IF(P0(K).GE.500E2)L5=K
569 IF(P0(K).GE.400E2)L4=K
570 IF(P0(K).GE.P300)LLFC=K
576 Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
577 DZA(K-1)=Z0(K)-Z0(K-1)
583 IF(LOW.GT.LLFC)GOTO 325
588 !...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
589 !...UNSTABLE AIR 50 TO 100 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
590 !...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
591 !...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 60 mb..
598 63 IF(DPTHMX.GT.6.E3)GOTO 64
608 !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
609 !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
614 ROCPQ=0.2854*(1.-0.28*Q0(NK))
615 THMIX=THMIX+DP(NK)*T0(NK)*(P00/P0(NK))**ROCPQ
616 QMIX=QMIX+DP(NK)*Q0(NK)
617 ZMIX=ZMIX+DP(NK)*Z0(NK)
618 17 PMIX=PMIX+DP(NK)*P0(NK)
623 ROCPQ=0.2854*(1.-0.28*QMIX)
624 TMIX=THMIX*(PMIX/P00)**ROCPQ
625 EMIX=QMIX*PMIX/(EP2+QMIX)
627 !...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL, PRESSURE
631 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
632 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX- &
634 TLCL=AMIN1(TLCL,TMIX)
635 TVLCL=TLCL*(1.+0.608*QMIX)
637 PLCL=P00*(TLCL/THMIX)**CPORQ
640 IF(PLCL.GE.P0(NK))GOTO 35
644 DLP=ALOG(PLCL/P0(K))/ALOG(P0(KLCL)/P0(K))
646 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
648 TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
649 QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
650 TVEN=TENV*(1.+0.608*QENV)
651 TVBAR=0.5*(TV0(K)+TVEN)
652 ! ZLCL=Z0(K)+R*TVBAR*ALOG(P0(K)/PLCL)/G
653 ZLCL=Z0(K)+(Z0(KLCL)-Z0(K))*DLP
655 !...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
656 !...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0AVG IS AN
657 !...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
658 !...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
659 !...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
660 !...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
661 !...SUCCESS AT GRID LENGTHS NEAR 25 km. FOR DIFFERENT GRID-LENGTHS,
662 !...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
663 !...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
665 WKLCL=0.02*ZLCL/2.5E3
666 WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)*DX/25.E3- &
670 DTLCL=4.64*WSIGNE*WABS**0.33
671 GDT=G*DTLCL*(ZLCL-Z0(LC))/(TV0(LC)+TVEN)
672 WLCL=1.+.5*WSIGNE*SQRT(ABS(GDT)+1.E-10)
673 IF(TLCL+DTLCL.GT.TENV)GOTO 45
674 IF(KPBL.GE.LLFC)GOTO 325
677 !...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
678 !...EQUIVALENT POTENTIAL TEMPERATURE
679 !...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
681 45 THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
682 EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
683 ES=ALIQ*EXP((TENV*BLIQ-CLIQ)/(TENV-DLIQ))
684 TVAVG=0.5*(TV0(KLCL)+TENV*(1.+0.608*QENV))
685 PLCL=P0(KLCL)*EXP(G/(R*TVAVG)*(Z0(KLCL)-ZLCL))
686 QESE=EP2*ES/(PLCL-ES)
687 GDT=G*DTLCL*(ZLCL-Z0(LC))/(TV0(LC)+TVEN)
688 WLCL=1.+.5*WSIGNE*SQRT(ABS(GDT)+1.E-10)
689 THTES(K)=TENV*(1.E5/PLCL)**(0.2854*(1.-0.28*QESE))* &
690 EXP((3374.6525/TENV-2.5403)*QESE*(1.+0.81*QESE))
692 IF(WLCL.LT.0.)GOTO 25
693 TVLCL=TLCL*(1.+0.608*QMIX)
694 RHOLCL=PLCL/(R*TVLCL)
699 !*******************************************************************
701 ! COMPUTE UPDRAFT PROPERTIES *
703 !*******************************************************************
706 !...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
715 !...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
716 !...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE BUOYANT ENERGY,
717 ! TRPPT IS THE TOTAL RATE OF PRECIPITATION PRODUCTION...
738 !...THE AMOUNT OF CONV AVAIL POT ENERGY (CAPE) IS CALCULATED WITH
739 ! RESPECT TO UNDILUTE PARCEL ASCENT; EQ POT TEMP OF UNDILUTE
740 ! PARCEL IS THTUDL, UNDILUTE TEMPERATURE IS GIVEN BY TUDL...
745 !...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
746 ! PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
747 ! FREEZING IS SPECIFIED TO BEGIN. WITHIN THE GLACIATION
748 ! INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
749 ! PREVIOUS MODEL LEVEL...
753 !...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
754 ! MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
755 ! MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
759 RATIO2(NK1)=RATIO2(NK)
761 !...UPDATE UPDRAFT PROPERTIES AT THE NEXT MODEL LVL TO REFLECT
762 ! ENTRAINMENT OF ENVIRONMENTAL AIR...
766 THETEU(NK1)=THETEU(NK)
771 CALL TPMIX(P0(NK1),THETEU(NK1),TU(NK1),QU(NK1),QLIQ(NK1), &
772 QICE(NK1),QNEWLQ,QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0, &
773 XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
774 TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
776 !...CHECK TO SEE IF UPDRAFT TEMP IS WITHIN THE FREEZING INTERVAL,
777 ! IF IT IS, CALCULATE THE FRACTIONAL CONVERSION TO GLACIATION
778 ! AND ADJUST QNEWLQ TO REFLECT THE GRADUAL CHANGE IN THETAU
779 ! SINCE THE LAST MODEL LEVEL...THE GLACIATION EFFECTS WILL BE
780 ! DETERMINED AFTER THE AMOUNT OF CONDENSATE AVAILABLE AFTER
781 ! PRECIP FALLOUT IS DETERMINED...TTFRZ IS THE TEMP AT WHICH
782 ! GLACIATION BEGINS, TBFRZ THE TEMP AT WHICH IT ENDS...
784 IF(TU(NK1).LE.TTFRZ.AND.IFLAG.LT.1)THEN
785 IF(TU(NK1).GT.TBFRZ)THEN
786 IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
787 FRC1=(TTEMP-TU(NK1))/(TTFRZ-TBFRZ)
788 R1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
790 FRC1=(TTEMP-TBFRZ)/(TTFRZ-TBFRZ)
795 QNEWIC=QNEWIC+QNEWLQ*R1*0.5
796 QNEWLQ=QNEWLQ-QNEWLQ*R1*0.5
797 EFFQ=(TTFRZ-TBFRZ)/(TTEMP-TBFRZ)
801 ! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
804 BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
805 BOTERM=2.*(Z0(NK1)-ZLCL)*G*BE/1.5
809 BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
810 BOTERM=2.*DZA(NK)*G*BE/1.5
811 ENTERM=2.*UER(NK)*WTW/UPOLD
815 CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM,RATE, &
816 QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1), G)
818 !...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
819 ! IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
825 ! UPDATE THE ABE FOR UNDILUTE ASCENT...
827 THTES(NK1)=T0(NK1)*(1.E5/P0(NK1))**(0.2854*(1.-0.28*QES(NK1))) &
829 EXP((3374.6525/T0(NK1)-2.5403)*QES(NK1)*(1.+0.81* &
831 UDLBE=((2.*THTUDL)/(THTES(NK)+THTES(NK1))-1.)*DZZ
832 IF(UDLBE.GT.0.)ABE=ABE+UDLBE*G
834 ! DETERMINE THE EFFECTS OF CLOUD GLACIATION IF WITHIN THE SPECIFIED
837 IF(FRC1.GT.1.E-6)THEN
838 CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QLIQ(NK1), &
839 QICE(NK1),RATIO2(NK1),TTFRZ,TBFRZ,QNWFRZ,RL,FRC1,EFFQ, &
840 IFLAG,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE &
844 ! CALL SUBROUTINE TO CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMP.
845 ! WITHIN GLACIATION INTERVAL, THETAE MUST BE CALCULATED WITH RESPECT TO
846 ! SAME DEGREE OF GLACIATION FOR ALL ENTRAINING AIR...
848 CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),RATIO2(NK1), &
849 RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
851 !...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
853 REI=VMFLCL*DP(NK1)*0.03/RAD
854 TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
856 !...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, NO
857 ! ENTRAINMENT IS ALLOWED AT THIS LEVEL...
859 IF(TVQU(NK1).LE.TV0(NK1))THEN
870 !...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL
871 ! AIR FOR ESTIMATION OF ENTRAINMENT AND DETRAINMENT RATES...
875 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
876 QTMP=F1*Q0(NK1)+F2*QU(NK1)
879 CALL TPMIX(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,QNEWLQ, &
880 QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ, &
881 DLIQ,AICE,BICE,CICE,DICE)
882 TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
883 IF(TU95.GT.TV0(NK1))THEN
891 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
892 QTMP=F1*Q0(NK1)+F2*QU(NK1)
895 CALL TPMIX(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,QNEWLQ, &
896 QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ, &
897 DLIQ,AICE,BICE,CICE,DICE)
898 TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
899 IF(TU10.EQ.TVQU(NK1))THEN
905 EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
906 EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
907 EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
908 IF(EQFRC(NK1).EQ.1)THEN
912 ELSEIF(EQFRC(NK1).EQ.0.)THEN
918 !...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
919 ! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
921 CALL PROF5(EQFRC(NK1),EE2,UD2)
929 !...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE
930 ! FRACTIONAL VALUES IN THE LAYER...
932 UER(NK1)=0.5*REI*(EE1+EE2)
933 UDR(NK1)=0.5*REI*(UD1+UD2)
935 !...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
936 ! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATION
938 55 IF(UMF(NK)-UDR(NK1).LT.10.)THEN
940 !...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL
941 ! UPDRAFT FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE
944 IF(UDLBE.GT.0.)ABE=ABE-UDLBE*G
946 ! WRITE(98,1015)P0(NK1)/100.
951 UPOLD=UMF(NK)-UDR(NK1)
955 !...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND ICE IN
956 ! THE DETRAINING UPDRAFT MASS...
958 DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
959 DETIC(NK1)=QICE(NK1)*UDR(NK1)
961 QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
962 THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
963 QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
964 QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
966 !...KFRZ IS THE HIGHEST MODEL LEVEL AT WHICH LIQUID CONDENSATE IS
967 ! GENERATING PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF LIQUID
968 ! PRECIP AT A GIVING MODEL LVL, PPTICE THE SAME FOR ICE, TRPPT IS
969 ! THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE CURRENT MODEL LEVEL
971 IF(ABS(RATIO2(NK1)-1.).GT.1.E-6)KFRZ=NK1
972 PPTLIQ(NK1)=QLQOUT(NK1)*(UMF(NK)-UDR(NK1))
973 PPTICE(NK1)=QICOUT(NK1)*(UMF(NK)-UDR(NK1))
974 TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
975 IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
978 !...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIU
979 ! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
980 ! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE
981 ! BETWEEN THE LET AND CLOUD TOP...
983 !...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL
984 ! VELOCITY FIRST BECOMES NEGATIVE...
989 !...IF CLOUD TOP HGT IS LESS THAN SPECIFIED MINIMUM HEIGHT, GO BACK AND
990 ! THE NEXT HIGHEST 60MB LAYER TO SEE IF A BIGGER CLOUD CAN BE OBTAINED
993 ! IF(CLDHGT.LT.4.E3.OR.ABE.LT.1.)THEN
994 IF(CLDHGT.LT.3.E3.OR.ABE.LT.1.)THEN
1006 !...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS
1007 ! FLUX THIS LEVEL...
1010 UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
1011 DETLQ(LTOP)=QLIQ(LTOP)*UDR(LTOP)*UPNEW/UPOLD
1012 DETIC(LTOP)=QICE(LTOP)*UDR(LTOP)*UPNEW/UPOLD
1013 TRPPT=TRPPT-(PPTLIQ(LTOP)+PPTICE(LTOP))
1019 ! BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
1024 DUMFDP=UMF(LET)/DPTT
1026 !...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
1027 ! RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
1031 UDR(NK)=DP(NK)*DUMFDP
1032 UMF(NK)=UMF(NK-1)-UDR(NK)
1033 DETLQ(NK)=QLIQ(NK)*UDR(NK)
1034 DETIC(NK)=QICE(NK)*UDR(NK)
1035 TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
1036 PPTLIQ(NK)=(UMF(NK-1)-UDR(NK))*QLQOUT(NK)
1037 PPTICE(NK)=(UMF(NK-1)-UDR(NK))*QICOUT(NK)
1038 TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
1041 !...SEND UPDRAFT CHARACTERISTICS TO OUTPUT FILES...
1045 !...EXTEND THE UPDRAFT MASS FLUX PROFILE DOWN TO THE SOURCE LAYER FOR
1046 ! THE UPDRAFT AIR...ALSO, DEFINE THETAE FOR LEVELS BELOW THE LCL...
1051 UMF(NK)=VMFLCL*DP(NK)/DPTHMX
1052 UER(NK)=VMFLCL*DP(NK)/DPTHMX
1053 ELSEIF(NK.LE.KPBL)THEN
1054 UER(NK)=VMFLCL*DP(NK)/DPTHMX
1055 UMF(NK)=UMF(NK-1)+UER(NK)
1060 TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
1081 EE=Q0(NK)*P0(NK)/(EP2+Q0(NK))
1083 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
1084 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T0(NK)-T00))*( &
1086 THTA=T0(NK)*(1.E5/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
1088 EXP((3374.6525/TSAT-2.5403)*Q0(NK)*(1.+0.81*Q0(NK)) &
1091 EXP((3374.6525/T0(NK)-2.5403)*QES(NK)*(1.+0.81* &
1099 !...DEFINE VARIABLES ABOVE CLOUD TOP...
1121 EMS(NK)=DP(NK)*DXSQ/G
1131 P150=P0(KLCL)-1.50E4
1134 EMS(NK)=DP(NK)*DXSQ/G
1137 !...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION
1140 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
1141 THTAU(NK)=TU(NK)*EXN(NK)
1142 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
1143 THTA0(NK)=T0(NK)*EXN(NK)
1145 !...LVF IS THE LEVEL AT WHICH MOISTURE FLUX IS ESTIMATED AS THE BASIS
1146 !...FOR PRECIPITATION EFFICIENCY CALCULATIONS...
1148 IF(P0(NK).GT.P150)LVF=NK
1151 USR=UMF(LVF+1)*(QU(LVF+1)+QLIQ(LVF+1)+QICE(LVF+1))
1152 USR=AMIN1(USR,TRPPT)
1153 IF(USR.LT.1.E-8)USR=TRPPT
1155 ! WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
1156 ! * TMIX-T00,PMIX,QMIX,ABE
1157 ! WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
1160 !...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
1161 !...AND MIDTROPOSPHERE IS USED.
1163 WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
1164 WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
1165 WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
1166 VCONV=.5*(WSPD(KLCL)+WSPD(L5))
1167 if (VCONV .gt. 0.) then
1174 TIMEC=AMAX1(1800.,TIMEC)
1175 TIMEC=AMIN1(3600.,TIMEC)
1179 !...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
1181 ! SHSIGN = CVMGT(1.,-1.,WSPD(LTOP).GT.WSPD(KLCL))
1182 IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
1187 VWS=(U0(LTOP)-U0(KLCL))*(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL))* &
1189 VWS=1.E3*SHSIGN*SQRT(VWS)/(Z0(LTOP)-Z0(LCL))
1190 PEF=1.591+VWS*(-.639+VWS*(9.53E-2-VWS*4.96E-3))
1194 !...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
1196 CBH=(ZLCL-Z0(1))*3.281E-3
1200 RCBH=.96729352+CBH*(-.70034167+CBH*(.162179896+CBH*(- &
1201 1.2569798E-2+CBH*(4.2772E-4-CBH*5.44E-6))))
1203 IF(CBH.GT.25)RCBH=2.4
1205 PEFCBH=AMIN1(PEFCBH,.9)
1207 !... MEAN PEF. IS USED TO COMPUTE RAINFALL.
1209 PEFF=.5*(PEF+PEFCBH)
1211 ! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1213 !*****************************************************************
1215 ! COMPUTE DOWNDRAFT PROPERTIES *
1217 !*****************************************************************
1219 !...LET DOWNDRAFT ORIGINATE AT THE LEVEL OF MINIMUM SATURATION EQUIVALEN
1220 !...POTENTIAL TEMPERATURE (SEQT) IN THE CLOUD LAYER, EXTEND DOWNWARD TO
1221 !...SURFACE, OR TO THE LAYER BELOW CLOUD BASE AT WHICH ENVIR SEQT IS LES
1222 !...THAN MIN SEQT IN THE CLOUD LAYER...LET DOWNDRAFT DETRAIN OVER A LAYE
1223 !...OF SPECIFIED PRESSURE-DEPTH (DPDD)...
1226 KSTART=MAX0(KPBL,KLCL)
1227 THTMIN=THTES(KSTART+1)
1229 DO 104 NK=KSTART+2,LTOP-1
1230 THTMIN=AMIN1(THTMIN,THTES(NK))
1231 IF(THTMIN.EQ.THTES(NK))KMIN=NK
1234 IF(RATIO2(LFS).GT.0.)CALL ENVIRTHT(P0(LFS),T0(LFS),Q0(LFS), &
1235 THETEE(LFS),0.,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1236 EQFRC(LFS)=(THTES(LFS)-THETEU(LFS))/(THETEE(LFS)-THETEU(LFS))
1237 EQFRC(LFS)=AMAX1(EQFRC(LFS),0.)
1238 EQFRC(LFS)=AMIN1(EQFRC(LFS),1.)
1239 THETED(LFS)=THTES(LFS)
1241 !...ESTIMATE THE EFFECT OF MELTING PRECIPITATION IN THE DOWNDRAFT...
1244 DTMLTD=0.5*(QU(KLCL)-QU(LTOP))*RLF/CP
1248 TZ(LFS)=T0(LFS)-DTMLTD
1249 ES=ALIQ*EXP((TZ(LFS)*BLIQ-CLIQ)/(TZ(LFS)-DLIQ))
1250 QS=EP2*ES/(P0(LFS)-ES)
1251 QD(LFS)=EQFRC(LFS)*Q0(LFS)+(1.-EQFRC(LFS))*QU(LFS)
1252 THTAD(LFS)=TZ(LFS)*(P00/P0(LFS))**(0.2854*(1.-0.28*QD(LFS)))
1253 IF(QD(LFS).GE.QS)THEN
1254 THETED(LFS)=THTAD(LFS)* &
1255 EXP((3374.6525/TZ(LFS)-2.5403)*QS*(1.+0.81*QS))
1257 CALL ENVIRTHT(P0(LFS),TZ(LFS),QD(LFS),THETED(LFS),0.,RL,EP2,ALIQ, &
1258 BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1262 IF(THETED(LFS).GT.THTES(ND).OR.ND.EQ.1)THEN
1265 !...IF DOWNDRAFT NEVER BECOMES NEGATIVELY BUOYANT OR IF IT
1266 !...IS SHALLOWER 50 mb, DON'T ALLOW IT TO OCCUR AT ALL...
1268 IF(NK.EQ.1.OR.(P0(LDB)-P0(LFS)).LT.50.E2)GOTO 141
1269 ! testing ---- no downdraft
1275 !...ALLOW DOWNDRAFT TO DETRAIN IN A SINGLE LAYER, BUT WITH DOWNDRAFT AIR
1276 !...TYPICALLY FLUSHED UP INTO HIGHER LAYERS AS ALLOWED IN THE TOTAL
1277 !...VERTICAL ADVECTION CALCULATIONS FARTHER DOWN IN THE CODE...
1285 ! IF(DPT.GT.DPDD)THEN
1287 ! FRC=(DPDD+DP(NK)-DPT)/DP(NK)
1290 ! IF(NK.EQ.LFS-1)THEN
1299 !...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX..
1301 TVD(LFS)=T0(LFS)*(1.+0.608*QES(LFS))
1302 RDD=P0(LFS)/(R*TVD(LFS))
1305 DER(LFS)=EQFRC(LFS)*DMF(LFS)
1307 DO 140 ND=LFS-1,LDB,-1
1311 DDR(ND)=-DMF(LDT+1)*DP(ND)*FRC/DPDD
1312 DMF(ND)=DMF(ND1)+DDR(ND)
1314 THETED(ND)=THETED(ND1)
1317 DER(ND)=DMF(LFS)*0.03*DP(ND)/RAD
1319 DMF(ND)=DMF(ND1)+DER(ND)
1320 IF(RATIO2(ND).GT.0.)CALL ENVIRTHT(P0(ND),T0(ND),Q0(ND), &
1321 THETEE(ND),0.,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1322 THETED(ND)=(THETED(ND1)*DMF(ND1)+THETEE(ND)*DER(ND))/DMF(ND)
1323 QD(ND)=(QD(ND1)*DMF(ND1)+Q0(ND)*DER(ND))/DMF(ND)
1328 !...CALCULATION AN EVAPORATION RATE FOR GIVEN MASS FLUX...
1332 TPDD(P0(ND),THETED(LDT),T0(ND),QS,QD(ND),1.0,XLV0,XLV1, &
1333 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1334 ES=ALIQ*EXP((TZ(ND)*BLIQ-CLIQ)/(TZ(ND)-DLIQ))
1335 QS=EP2*ES/(P0(ND)-ES)
1336 DSSDT=(CLIQ-BLIQ*DLIQ)/((TZ(ND)-DLIQ)*(TZ(ND)-DLIQ))
1338 DTMP=RL*QS*(1.-RHBC)/(CP+RL*RHBC*QS*DSSDT)
1340 ES=RHBC*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
1341 QSRH=EP2*ES/(P0(ND)-ES)
1343 !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
1344 !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
1346 IF(QSRH.LT.QD(ND))THEN
1348 ! T1RH=T1+(QS-QSRH)*RL/CP
1353 TDER=TDER+(QS-QD(ND))*DDR(ND)
1355 135 THTAD(ND)=TZ(ND)*(P00/P0(ND))**(0.2854*(1.-0.28*QD(ND)))
1357 !...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
1358 !...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
1360 141 IF(TDER.LT.1.)THEN
1362 3004 FORMAT(' ','I=',I3,2X,'J=',I3)
1381 !...ADJUST DOWNDRAFT MASS FLUX SO THAT EVAPORATION RATE IN DOWNDRAFT IS
1382 !...CONSISTENT WITH PRECIPITATION EFFICIENCY RELATIONSHIP...
1384 DEVDMF=TDER/DMF(LFS)
1389 !...PPR IS THE TOTAL AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE
1390 !...UPDRAFT FROM CLOUD BASE TO THE LFS...UPDRAFT MASS FLUX WILL BE
1391 !...INCREASED UP TO THE LFS TO ACCOUNT FOR UPDRAFT AIR MIXING WITH
1392 !...ENVIRONMENTAL AIR TO THE UPDRAFT, SO PPR WILL INCREASE
1393 !...PROPORTIONATELY...
1396 132 PPR=PPR+PPTLIQ(NM)+PPTICE(NM)
1398 DPPTDF=(1.-PEFF)*PPR*(1.-EQFRC(LFS))/UMF(LFS)
1403 !...CNDTNF IS THE AMOUNT OF CONDENSATE TRANSFERRED ALONG WITH UPDRAFT
1404 !...MASS THE DOWNDRAFT AT THE LFS...
1406 CNDTNF=(QLIQ(LFS)+QICE(LFS))*(1.-EQFRC(LFS))
1407 DMFLFS=RCED/(DEVDMF+DPPTDF+CNDTNF)
1408 IF(DMFLFS.GT.0.)THEN
1413 !...DDINC IS THE FACTOR BY WHICH TO INCREASE THE FIRST-GUESS DOWNDRAFT
1414 !...MASS FLUX TO SATISFY THE PRECIP EFFICIENCY RELATIONSHIP, UPDINC IS T
1415 !...WHICH TO INCREASE THE UPDRAFT MASS FLUX BELOW THE LFS TO ACCOUNT FOR
1416 !...TRANSFER OF MASS FROM UPDRAFT TO DOWNDRAFT...
1418 ! DDINC=DMFLFS/DMF(LFS)
1420 UPDINC=(UMF(LFS)-(1.-EQFRC(LFS))*DMFLFS)/UMF(LFS)
1422 !...LIMIT UPDINC TO LESS THAN OR EQUAL TO 1.5...
1424 IF(UPDINC.GT.1.5)THEN
1426 DMFLFS2=UMF(LFS)*(UPDINC-1.)/(EQFRC(LFS)-1.)
1427 RCED2=DMFLFS2*(DEVDMF+DPPTDF+CNDTNF)
1428 PPTFLX=PPTFLX+(RCED-RCED2)
1436 DDINC=DMFLFS/DMF(LFS)
1438 DMF(NK)=DMF(NK)*DDINC
1439 DER(NK)=DER(NK)*DDINC
1440 DDR(NK)=DDR(NK)*DDINC
1442 CPR=TRPPT+PPR*(UPDINC-1.)
1443 PPTFLX=PPTFLX+PEFF*PPR*(UPDINC-1.)
1447 !...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
1448 ! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
1449 ! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
1452 UMF(NK)=UMF(NK)*UPDINC
1453 UDR(NK)=UDR(NK)*UPDINC
1454 UER(NK)=UER(NK)*UPDINC
1455 PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
1456 PPTICE(NK)=PPTICE(NK)*UPDINC
1457 DETLQ(NK)=DETLQ(NK)*UPDINC
1458 155 DETIC(NK)=DETIC(NK)*UPDINC
1460 !...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
1483 DO 158 NK=LDT+1,LFS-1
1489 !...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE
1490 ! INFLOW INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN
1491 ! IS AVAILABLE IN THAT LAYER INITIALLY...
1496 IF((UER(NK)-DER(NK)).GT.0.)AINCM1=EMS(NK)/((UER(NK)-DER(NK))* &
1498 AINCMX=AMIN1(AINCMX,AINCM1)
1501 IF(AINCMX.LT.AINC)AINC=AINCMX
1503 !...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRFT AND DOWNDRFT...THEY
1504 !...WILL ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE
1505 !...STABILIZATION CLOSURE...
1511 DETLQ2(NK)=DETLQ(NK)
1512 DETIC2(NK)=DETIC(NK)
1523 IF(AINC/AINCMX.GT.0.999)THEN
1530 !*****************************************************************
1532 ! COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE *
1534 !*****************************************************************
1536 !...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
1537 !...SATISFY MASS CONTINUITY...
1542 DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
1544 OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
1545 DTT1=0.75*DP(NK-1)/(ABS(OMG(NK))+1.E-10)
1552 NSTEP=NINT(TIMEC/DTT+1)
1553 DTIME=TIMEC/FLOAT(NSTEP)
1554 FXM(NK)=OMG(NK)*DXSQ/G
1557 !...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
1561 !...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED
1562 !...SIGN OF OMEGA...
1570 IF(OMG(NK).LE.0.)THEN
1571 THFXBOT(NK)=-FXM(NK)*THPA(NK-1)
1572 QFXBOT(NK)=-FXM(NK)*QPA(NK-1)
1573 THFXTOP(NK-1)=THFXTOP(NK-1)-THFXBOT(NK)
1574 QFXTOP(NK-1)=QFXTOP(NK-1)-QFXBOT(NK)
1576 THFXBOT(NK)=-FXM(NK)*THPA(NK)
1577 QFXBOT(NK)=-FXM(NK)*QPA(NK)
1578 THFXTOP(NK-1)=THFXTOP(NK-1)-THFXBOT(NK)
1579 QFXTOP(NK-1)=QFXTOP(NK-1)-QFXBOT(NK)
1583 !...UPDATE THE THETA AND QV VALUES AT EACH LEVEL..
1586 THPA(NK)=THPA(NK)+(THFXBOT(NK)+UDR(NK)*THTAU(NK)+DDR(NK)* &
1587 THTAD(NK)+THFXTOP(NK)-(UER(NK)-DER(NK))*THTA0(NK))* &
1589 QPA(NK)=QPA(NK)+(QFXBOT(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK)+ &
1590 QFXTOP(NK)-(UER(NK)-DER(NK))*Q0(NK))*DTIME*EMSD(NK)
1599 !...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO,
1600 !...BORROW MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO.
1603 IF(QG(NK).LT.0.)THEN
1605 CALL wrf_error_fatal ( 'module_cu_kf.F: problem with kf scheme: qg = 0 at the surface' )
1608 IF(NK.EQ.LTOP)NK1=KLCL
1609 TMA=QG(NK1)*EMS(NK1)
1610 TMB=QG(NK-1)*EMS(NK-1)
1611 TMM=(QG(NK)-1.E-9)*EMS(NK)
1612 BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
1613 ACOEFF=BCOEFF*TMA/TMB
1617 QVDIFF=(QG(NK1)-TMA*EMSD(NK1))*100./QG(NK1)
1618 IF(ABS(QVDIFF).GT.1.)THEN
1619 PRINT *,'--WARNING-- CLOUD BASE WATER VAPOR CHANGES BY ', &
1621 ' PERCENT WHEN MOISTURE IS BORROWED TO PREVENT NEG VALUES', &
1626 QG(NK1)=TMA*EMSD(NK1)
1627 QG(NK-1)=TMB*EMSD(NK-1)
1630 TOPOMG=(UDR(LTOP)-UER(LTOP))*DP(LTOP)*EMSD(LTOP)
1631 IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
1632 ! WRITE(98,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME;'
1633 ! * ,'TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1634 WRITE(6,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME;' &
1635 ,'TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1640 !...CONVERT THETA TO T...
1645 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QG(NK)))
1646 TG(NK)=THTAG(NK)/EXN(NK)
1647 TVG(NK)=TG(NK)*(1.+0.608*QG(NK))
1650 !*******************************************************************
1652 ! COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY. *
1654 !*******************************************************************
1656 !...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
1662 ROCPQ=0.2854*(1.-0.28*QG(NK))
1663 THMIX=THMIX+DP(NK)*TG(NK)*(P00/P0(NK))**ROCPQ
1664 QMIX=QMIX+DP(NK)*QG(NK)
1665 217 PMIX=PMIX+DP(NK)*P0(NK)
1669 ROCPQ=0.2854*(1.-0.28*QMIX)
1670 TMIX=THMIX*(PMIX/P00)**ROCPQ
1671 ES=ALIQ*EXP((TMIX*BLIQ-CLIQ)/(TMIX-DLIQ))
1674 !...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
1678 CPM=CP*(1.+0.887*QMIX)
1679 DSSDT=QS*(CLIQ-BLIQ*DLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
1680 DQ=(QMIX-QS)/(1.+RL*DSSDT/CPM)
1683 ROCPQ=0.2854*(1.-0.28*QMIX)
1684 THMIX=TMIX*(P00/PMIX)**ROCPQ
1689 EMIX=QMIX*PMIX/(EP2+QMIX)
1690 TLOG=ALOG(EMIX/ALIQ)
1691 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
1692 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX- &
1694 TLCL=AMIN1(TLCL,TMIX)
1696 PLCL=P00*(TLCL/THMIX)**CPORQ
1698 TVLCL=TLCL*(1.+0.608*QMIX)
1701 235 IF(PLCL.GE.P0(NK))GOTO 240
1703 DLP=ALOG(PLCL/P0(K))/ALOG(P0(KLCL)/P0(K))
1705 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
1707 TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
1708 QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
1709 TVEN=TENV*(1.+0.608*QENV)
1710 TVBAR=0.5*(TVG(K)+TVEN)
1711 ! ZLCL=Z0(K)+R*TVBAR*ALOG(P0(K)/PLCL)/G
1712 ZLCL=Z0(K)+(Z0(KLCL)-Z0(K))*DLP
1713 TVAVG=0.5*(TVEN+TG(KLCL)*(1.+0.608*QG(KLCL)))
1714 PLCL=P0(KLCL)*EXP(G/(R*TVAVG)*(Z0(KLCL)-ZLCL))
1715 THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
1716 EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
1717 ES=ALIQ*EXP((TENV*BLIQ-CLIQ)/(TENV-DLIQ))
1718 QESE=EP2*ES/(PLCL-ES)
1719 THTESG(K)=TENV*(1.E5/PLCL)**(0.2854*(1.-0.28*QESE))* &
1720 EXP((3374.6525/TENV-2.5403)*QESE*(1.+0.81*QESE))
1722 !...COMPUTE ADJUSTED ABE(ABEG).
1728 ES=ALIQ*EXP((TG(NK1)*BLIQ-CLIQ)/(TG(NK1)-DLIQ))
1729 QESE=EP2*ES/(P0(NK1)-ES)
1730 THTESG(NK1)=TG(NK1)*(1.E5/P0(NK1))**(0.2854*(1.-0.28*QESE))* &
1731 EXP((3374.6525/TG(NK1)-2.5403)*QESE*(1.+0.81*QESE) &
1733 ! DZZ=CVMGT(Z0(KLCL)-ZLCL,DZA(NK),NK.EQ.K)
1739 BE=((2.*THTUDL)/(THTESG(NK1)+THTESG(NK))-1.)*DZZ
1740 245 IF(BE.GT.0.)ABEG=ABEG+BE*G
1742 !...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
1743 !...THE PERIOD TIMEC...
1746 ! WRITE(98,1060)FABE
1749 DABE=AMAX1(ABE-ABEG,0.1*ABE)
1750 FABE=ABEG/(ABE+1.E-8)
1752 ! WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS '
1753 ! *,'GRID POINT; NO CONVECTION ALLOWED!'
1757 DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
1766 IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
1767 ! WRITE(98,1055)FABE
1770 IF(FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB)GOTO 265
1771 IF(NCOUNT.GT.10)THEN
1772 ! WRITE(98,1060)FABE
1776 !...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE
1777 !...CONVECTIVE MASS FLUX BY THE FACTOR AINC:
1782 AINC=AINC*STAB*ABE/(DABE+1.E-8)
1784 255 AINC=AMIN1(AINCMX,AINC)
1785 !...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION
1786 !...WILL BE MINIMAL SO JUST IGNORE IT...
1787 IF(AINC.LT.0.05)GOTO 325
1788 ! AINC=AMAX1(AINC,0.05)
1791 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABEOLD,AINCOLD
1793 UMF(NK)=UMF2(NK)*AINC
1794 DMF(NK)=DMF2(NK)*AINC
1795 DETLQ(NK)=DETLQ2(NK)*AINC
1796 DETIC(NK)=DETIC2(NK)*AINC
1797 UDR(NK)=UDR2(NK)*AINC
1798 UER(NK)=UER2(NK)*AINC
1799 DER(NK)=DER2(NK)*AINC
1800 DDR(NK)=DDR2(NK)*AINC
1803 !...GO BACK UP FOR ANOTHER ITERATION...
1808 !...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
1811 !...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
1813 !...FRC2 IS THE FRACTION OF TOTAL CONDENSATE
1814 !...GENERATED THAT GOES INTO PRECIPITIATION
1815 FRC2=PPTFLX/(CPR*AINC)
1821 RAINFB(NK)=PPTLIQ(NK)*AINC*FBFRC*FRC2
1822 SNOWFB(NK)=PPTICE(NK)*AINC*FBFRC*FRC2
1826 !...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH
1827 !...LAYER BASED ON THE SIGN OF OMEGA...
1840 IF(OMG(NK).LE.0.)THEN
1841 QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
1842 QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
1843 QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
1844 QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
1845 QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
1846 QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
1847 QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
1848 QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
1850 QLFXOUT(NK)=FXM(NK)*QLPA(NK)
1851 QIFXOUT(NK)=FXM(NK)*QIPA(NK)
1852 QRFXOUT(NK)=FXM(NK)*QRPA(NK)
1853 QSFXOUT(NK)=FXM(NK)*QSPA(NK)
1854 QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
1855 QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
1856 QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
1857 QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
1861 !...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
1864 QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))*DTIME* &
1866 QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))*DTIME* &
1868 QRPA(NK)=QRPA(NK)+(QRFXIN(NK)+QLQOUT(NK)*UDR(NK)-QRFXOUT(NK) &
1869 +RAINFB(NK))*DTIME*EMSD(NK)
1870 QSPA(NK)=QSPA(NK)+(QSFXIN(NK)+QICOUT(NK)*UDR(NK)-QSFXOUT(NK) &
1871 +SNOWFB(NK))*DTIME*EMSD(NK)
1880 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABE,AINC
1882 !...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
1885 WRITE(6,1070)' P ',' DP ',' DT K/D ',' DR K/D ',' OMG ', &
1886 ' DOMGDP ',' UMF ',' UER ',' UDR ',' DMF ',' DER ' &
1887 ,' DDR ',' EMS ',' W0 ',' DETLQ ',' DETIC '
1889 DTT=(TG(K)-T0(K))*86400./TIMEC
1891 DR=-(QG(K)-Q0(K))*RL*86400./(TIMEC*CP)
1892 UDFRC=UDR(K)*TIMEC*EMSD(K)
1893 UEFRC=UER(K)*TIMEC*EMSD(K)
1894 DDFRC=DDR(K)*TIMEC*EMSD(K)
1895 DEFRC=-DER(K)*TIMEC*EMSD(K)
1896 WRITE (6,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)* &
1897 1.E4,UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC &
1898 ,DDFRC,EMS(K)/1.E11,W0AVG1D(K)*1.E2,DETLQ(K) &
1899 *TIMEC*EMSD(K)*1.E3,DETIC(K)*TIMEC*EMSD(K)* &
1902 WRITE(6,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG', &
1903 'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
1907 IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
1909 IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
1910 ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
1911 QGS=ES*EP2/(P0(K)-ES)
1914 WRITE (6,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC &
1915 ,TDC,Q0(K)*1000.,QG(K)*1000.,(QG(K)-Q0(K))* &
1916 1000.,QU(K)*1000.,QD(K)*1000.,QLG(K)*1000., &
1917 QIG(K)*1000.,QRG(K)*1000.,QSG(K)*1000.,RH0,RHG
1920 !...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
1921 !...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
1925 WRITE ( wrf_err_message , 1115 ) &
1926 Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
1927 U0(K),V0(K),DP(K)/100.,W0AVG1D(K)
1928 CALL wrf_message ( TRIM( wrf_err_message ) )
1930 CALL wrf_error_fatal ( 'module_cu_kf.F: KAIN-FRITSCH' )
1933 CNDTNF=(1.-EQFRC(LFS))*(QLIQ(LFS)+QICE(LFS))*DMF(LFS)
1934 ! WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
1936 ! EVALUATE MOISTURE BUDGET...
1943 QINIT=QINIT+Q0(NK)*EMS(NK)
1944 QFNL=QFNL+QG(NK)*EMS(NK)
1945 QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
1947 QFNL=QFNL+PPTFLX*TIMEC*(1.-FBFRC)
1948 ERR2=(QFNL-QINIT)*100./QINIT
1949 ! WRITE(98,1110)QINIT,QFNL,ERR2
1950 ! IF(ABS(ERR2).GT.0.05)STOP 'QVERR'
1951 IF(ABS(ERR2).GT.0.05)CALL wrf_error_fatal( 'module_cu_kf.F: QVERR' )
1952 RELERR=ERR2*QINIT/(PPTFLX*TIMEC+1.E-10)
1953 ! WRITE(98,1120)RELERR
1954 ! WRITE(98,*)'TDER, CPR, USR, TRPPT =',
1955 ! *TDER,CPR*AINC,USR*AINC,TRPPT*AINC
1957 !...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
1959 !...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
1960 !...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
1962 IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
1963 NCA(I,J)=FLOAT(NIC)*DT
1965 ! IF(IMOIST.NE.2)THEN
1967 !...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMAT
1968 !...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
1969 !...NOTE: THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
1970 !...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
1973 ! RLC=XLV0-XLV1*TG(K)
1974 ! RLS=XLS0-XLS1*TG(K)
1975 ! CPM=CP*(1.+0.887*QG(K))
1976 ! TG(K)=TG(K)-(RLC*(QLG(K)+QRG(K))+RLS*(QIG(K)+QSG(K)))/CPM
1977 ! QG(K)=QG(K)+(QLG(K)+QRG(K)+QIG(K)+QSG(K))
1983 IF(.NOT. qi_flag .and. warm_rain)THEN
1985 !...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
1987 CPM=CP*(1.+0.887*QG(K))
1988 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
1989 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
1991 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
1993 ELSEIF(.NOT. qi_flag .and. .not. warm_rain)THEN
1995 !...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROME
1996 !...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
1998 CPM=CP*(1.+0.887*QG(K))
2000 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
2002 TG(K)=TG(K)+(QLG(K)+QRG(K))*RLF/CPM
2004 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
2006 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
2008 ELSEIF(qi_flag) THEN
2010 !...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE
2011 !...TENDENCY OF HYDROMETEORS DIRECTLY...
2013 DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
2014 DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
2015 DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
2017 DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
2019 DQIDT(K)=DQIDT(K)+(QSG(K)-QS0(K))/TIMEC
2022 CALL wrf_error_fatal ( 'module_cu_kf: THIS COMBINATION OF IMOIST, IICE NOT ALLOWED' )
2025 DTDT(K)=(TG(K)-T0(K))/TIMEC
2026 DQDT(K)=(QG(K)-Q0(K))/TIMEC
2029 ! RAINCV is in the unit of mm
2031 PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
2032 RAINCV(I,J)=DT*PRATEC(I,J)
2035 909 FORMAT(' CONVECTIVE RAINFALL =',F8.4,' CM')
2039 1000 FORMAT(' ',10A8)
2040 1005 FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
2041 1010 FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
2042 1015 FORMAT(' ','ALL REMAINING MASS DETRAINS BELOW ',F4.0,' MB')
2043 1025 FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M', &
2044 ' DTLCL=',F5.2,' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=', &
2045 I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1, &
2047 1030 FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =', &
2048 E12.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =', &
2050 1035 FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'WKL=' &
2052 1040 FORMAT(' ','PRECIP EFF = 100%, ENVIR CANNOT SUPPORT DOWND' &
2054 !1045 FORMAT('NUMBER OF DOWNDRAFT ITERATIONS EXCEEDS 10...PPTFLX' &
2055 ! ' IS DIFFERENT FROM THAT GIVEN BY PRECIP EFF RELATION')
2056 ! FLIC HAS TROUBLE WITH THIS ONE.
2057 1045 FORMAT('NUMBER OF DOWNDRAFT ITERATIONS EXCEEDS 10')
2058 1050 FORMAT(' ','LCOUNT= ',I3,' PPTFLX/CPR, PEFF= ',F5.3,1X,F5.3, &
2059 'DMF(LFS)/UMF(LCL)= ',F5.3)
2060 1055 FORMAT(/'*** DEGREE OF STABILIZATION =',F5.3,', NO MORE MASS F' &
2062 !1060 FORMAT(/' ITERATION DOES NOT CONVERGE TO GIVE THE SPECIFIED ' &
2063 ! 'DEGREE OF STABILIZATION! FABE= ',F6.4)
2064 1060 FORMAT(/' ITERATION DOES NOT CONVERGE. FABE= ',F6.4)
2066 1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3)
2067 1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, NSTEP=',F5.0,I3, &
2068 'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
2069 1085 FORMAT (A3,16A7,2A8)
2070 1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3)
2071 1095 FORMAT(' ',' PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ', &
2073 1105 FORMAT(' ','NET LATENT HEAT RELEASE =',E12.5,' ACTUAL HEATING =', &
2074 E12.5,' J/KG-S, DIFFERENCE = ',F9.3,'PERCENT')
2075 1110 FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5, &
2076 ' TOTAL WATER CHANGE =',F8.2,'PERCENT')
2077 1115 FORMAT (2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4 &
2079 1120 FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3, &
2082 END SUBROUTINE KFPARA
2084 !-----------------------------------------------------------------------
2085 SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
2086 QNEWIC,QLQOUT,QICOUT,G)
2087 !-----------------------------------------------------------------------
2089 !-----------------------------------------------------------------------
2090 ! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
2091 ! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
2092 ! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
2093 ! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
2094 ! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
2096 REAL, INTENT(IN ) :: G
2097 REAL, INTENT(IN ) :: DZ,BOTERM,ENTERM,RATE
2098 REAL, INTENT(INOUT) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
2100 REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
2105 ! ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY C
2106 ! BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
2109 QEST=0.5*(QTOT+QNEW)
2110 G1=WTW+BOTERM-ENTERM-2.*G*DZ*QEST/1.5
2112 WAVG=(SQRT(WTW)+SQRT(G1))/2.
2115 ! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
2116 ! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
2117 ! IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
2118 ! SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...
2120 RATIO3=QNEWLQ/(QNEW+1.E-10)
2124 RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-10)
2125 QTOT=QTOT*EXP(-CONV)
2127 ! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
2128 ! PARCEL AT THIS LEVEL...
2132 QICOUT=(1.-RATIO4)*DQ
2134 ! ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
2135 ! LATE VERTICAL VELOCITY
2137 PPTDRG=0.5*(OLDQ+QTOT-0.2*QNEW)
2138 WTW=WTW+BOTERM-ENTERM-2.*G*DZ*PPTDRG/1.5
2140 ! DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
2141 ! DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...
2143 QLIQ=RATIO4*QTOT+RATIO3*0.4*QNEW
2144 QICE=(1.-RATIO4)*QTOT+(1.-RATIO3)*0.4*QNEW
2148 END SUBROUTINE CONDLOAD
2150 !-----------------------------------------------------------------------
2151 SUBROUTINE DTFRZNEW(TU,P,THTEU,QVAP,QLIQ,QICE,RATIO2,TTFRZ,TBFRZ, &
2152 QNWFRZ,RL,FRC1,EFFQ,IFLAG,XLV0,XLV1,XLS0,XLS1, &
2153 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2154 !-----------------------------------------------------------------------
2156 !-----------------------------------------------------------------------
2157 REAL, INTENT(IN ) :: XLV0,XLV1
2158 REAL, INTENT(IN ) :: P,TTFRZ,TBFRZ,EFFQ,XLS0,XLS1,EP2,ALIQ, &
2159 BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE
2160 REAL, INTENT(INOUT) :: TU,THTEU,QVAP,QLIQ,QICE,RATIO2, &
2162 INTEGER, INTENT(INOUT) :: IFLAG
2164 REAL :: CCP,RV,C5,QLQFRZ,QNEW,ESLIQ,ESICE,RLC,RLS,PI,ES,RLF,A, &
2165 B,C,DQVAP,DTFRZ,TU1,QVAP1
2166 !-----------------------------------------------------------------------
2168 !...ALLOW GLACIATION OF THE UPDRAFT TO OCCUR AS AN APPROXIMATELY LINEAR
2169 ! FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE TTFRZ TO TBFRZ...
2175 !...ADJUST THE LIQUID WATER CONCENTRATIONS FROM FRESH CONDENSATE AND THA
2176 ! BROUGHT UP FROM LOWER LEVELS TO AN AMOUNT THAT WOULD BE PRESENT IF N
2177 ! LIQUID WATER HAD FROZEN THUS FAR...THIS IS NECESSARY BECAUSE THE
2178 ! EXPRESSION FOR TEMP CHANGE IS MULTIPLIED BY THE FRACTION EQUAL TO TH
2179 ! PARCEL TEMP DECREASE SINCE THE LAST MODEL LEVEL DIVIDED BY THE TOTAL
2180 ! GLACIATION INTERVAL, SO THAT EFFECTIVELY THIS APPROXIMATELY ALLOWS A
2181 ! AMOUNT OF LIQUID WATER TO FREEZE WHICH IS EQUAL TO THIS SAME FRACTIO
2182 ! OF THE LIQUID WATER THAT WAS PRESENT BEFORE THE GLACIATION PROCESS W
2183 ! INITIATED...ALSO, TO ALLOW THETAU TO CONVERT APPROXIMATELY LINEARLY
2184 ! ITS VALUE WITH RESPECT TO ICE, WE NEED TO ALLOW A PORTION OF THE FRE
2185 ! CONDENSATE TO CONTRIBUTE TO THE GLACIATION PROCESS; THE FRACTIONAL
2186 ! AMOUNT THAT APPLIES TO THIS PORTION IS 1/2 OF THE FRACTIONAL AMOUNT
2187 ! FROZEN OF THE "OLD" CONDENSATE BECAUSE THIS FRESH CONDENSATE IS ONLY
2188 ! PRODUCED GRADUALLY OVER THE LAYER...NOTE THAT IN TERMS OF THE DYNAMI
2189 ! OF THE PRECIPITATION PROCESS, IE. PRECIPITATION FALLOUT, THIS FRACTI
2190 ! AMNT OF FRESH CONDENSATE HAS ALREADY BEEN INCLUDED IN THE ICE CATEGO
2193 QNEW=QNWFRZ*EFFQ*0.5
2194 ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2195 ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
2196 RLC=2.5E6-2369.276*(TU-273.16)
2197 RLS=2833922.-259.532*(TU-273.16)
2199 CCP=1005.7*(1.+0.89*QVAP)
2201 ! A = D(ES)/DT IS THAT CALCULATED FROM BUCK`S (1981) EMPIRICAL FORMULAS
2202 ! FOR SATURATION VAPOR PRESSURE...
2204 A=(CICE-BICE*DICE)/((TU-DICE)*(TU-DICE))
2207 DQVAP=B*(ESLIQ-ESICE)/(RLS+RLS*C)-RLF*(QLQFRZ+QNEW)/(RLS+RLS/C)
2208 DTFRZ=(RLF*(QLQFRZ+QNEW)+B*(ESLIQ-ESICE))/(CCP+A*B*ESICE)
2212 QVAP=QVAP-FRC1*DQVAP
2213 ES=QVAP*P/(EP2+QVAP)
2214 ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2215 ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
2216 RATIO2=(ESLIQ-ES)/(ESLIQ-ESICE)
2218 ! TYPICALLY, RATIO2 IS VERY CLOSE TO (TTFRZ-TU)/(TTFRZ-TBFRZ), USUALLY
2219 ! WITHIN 1% (USING TU BEFORE GALCIATION EFFECTS ARE APPLIED); IF THE
2220 ! INITIAL UPDRAFT TEMP IS BELOW TBFRZ AND RATIO2 IS STILL LESS THAN 1,
2221 ! AN ADJUSTMENT TO FRC1 AND RATIO2 IS INTRODUCED SO THAT GLACIATION
2222 ! EFFECTS ARE NOT UNDERESTIMATED; CONVERSELY, IF RATIO2 IS GREATER THAN
2223 ! FRC1 IS ADJUSTED SO THAT GLACIATION EFFECTS ARE NOT OVERESTIMATED...
2225 IF(IFLAG.GT.0.AND.RATIO2.LT.1)THEN
2226 FRC1=FRC1+(1.-RATIO2)
2228 QVAP=QVAP1-FRC1*DQVAP
2233 IF(RATIO2.GT.1.)THEN
2234 FRC1=FRC1-(RATIO2-1.)
2235 FRC1=AMAX1(0.0,FRC1)
2237 QVAP=QVAP1-FRC1*DQVAP
2242 ! CALCULATE A HYBRID VALUE OF THETAU, ASSUMING THAT THE LATENT HEAT OF
2243 ! VAPORIZATION/SUBLIMATION CAN BE ESTIMATED USING THE SAME WEIGHTING
2244 ! FUNCTION AS THAT USED TO CALCULATE SATURATION VAPOR PRESSURE, CALCU-
2245 ! LATE NEW LIQUID WATER AND ICE CONCENTRATIONS...
2249 RL=RATIO2*RLS+(1.-RATIO2)*RLC
2250 PI=(1.E5/P)**(0.2854*(1.-0.28*QVAP))
2251 THTEU=TU*PI*EXP(RL*QVAP*C5/TU*(1.+0.81*QVAP))
2253 QICE=QICE+FRC1*DQVAP+QLIQ
2256 QICE=QICE+FRC1*(DQVAP+QLQFRZ)
2257 QLIQ=QLIQ-FRC1*QLQFRZ
2261 END SUBROUTINE DTFRZNEW
2263 !-----------------------------------------------------------------------
2264 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2265 ! THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN
2266 ! DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN F
2267 ! HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMA
2268 ! TABLES ED. BY ABRAMOWITZ AND STEGUN, NAT L BUREAU OF STANDARDS APPLI
2269 ! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
2272 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2273 !***********************************************************************
2274 !***** GAUSSIAN TYPE MIXING PROFILE....******************************
2275 SUBROUTINE PROF5(EQ,EE,UD)
2276 !-----------------------------------------------------------------------
2278 !-----------------------------------------------------------------------
2279 REAL, INTENT(IN ) :: EQ
2280 REAL, INTENT(INOUT) :: EE,UD
2281 REAL :: SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
2283 DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676, &
2284 0.9372980,0.33267,0.166666667,0.202765151/
2291 C1=A1*T1+A2*T1*T1+A3*T1*T1*T1
2292 C2=A1*T2+A2*T2*T2+A3*T2*T2*T2
2294 EE=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
2295 UD=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*(0.5+EQ*EQ/2.- &
2298 EE=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
2299 UD=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*(0.5+EQ* &
2305 END SUBROUTINE PROF5
2307 !-----------------------------------------------------------------------
2308 SUBROUTINE TPMIX(P,THTU,TU,QU,QLIQ,QICE,QNEWLQ,QNEWIC,RATIO2,RL, &
2309 XLV0,XLV1,XLS0,XLS1, &
2310 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2311 !-----------------------------------------------------------------------
2313 !-----------------------------------------------------------------------
2314 REAL, INTENT(IN ) :: XLV0,XLV1
2315 REAL, INTENT(IN ) :: P,THTU,RATIO2,RL,XLS0, &
2316 XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,&
2318 REAL, INTENT(INOUT) :: QU,QLIQ,QICE,TU,QNEWLQ,QNEWIC
2319 REAL :: ES,QS,PI,THTGS,F0,T1,T0,C5,RV,ESLIQ,ESICE,F1,DT,QNEW, &
2320 DQ, QTOT,DQICE,DQLIQ,RLL,CCP
2322 !-----------------------------------------------------------------------
2324 !...THIS SUBROUTINE ITERATIVELY EXTRACTS WET-BULB TEMPERATURE FROM EQUIV
2325 ! POTENTIAL TEMPERATURE, THEN CHECKS TO SEE IF SUFFICIENT MOISTURE IS
2326 ! AVAILABLE TO ACHIEVE SATURATION...IF NOT, TEMPERATURE IS ADJUSTED
2327 ! ACCORDINGLY, IF SO, THE RESIDUAL LIQUID WATER/ICE CONCENTRATION IS
2333 ! ITERATE TO FIND WET BULB TEMPERATURE AS A FUNCTION OF EQUIVALENT POT
2334 ! TEMP AND PRS, ASSUMING SATURATION VAPOR PRESSURE...RATIO2 IS THE DEG
2337 IF(RATIO2.LT.1.E-6)THEN
2338 ES=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2340 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2341 THTGS=TU*PI*EXP((3374.6525/TU-2.5403)*QS*(1.+0.81*QS))
2342 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
2343 ES=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
2345 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2346 THTGS=TU*PI*EXP((3114.834/TU-0.278296)*QS*(1.+0.81*QS))
2348 ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2349 ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
2350 ES=(1.-RATIO2)*ESLIQ+RATIO2*ESICE
2352 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2353 THTGS=TU*PI*EXP(RL*QS*C5/TU*(1.+0.81*QS))
2359 90 IF(RATIO2.LT.1.E-6)THEN
2360 ES=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
2362 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2363 THTGS=T1*PI*EXP((3374.6525/T1-2.5403)*QS*(1.+0.81*QS))
2364 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
2365 ES=AICE*EXP((BICE*T1-CICE)/(T1-DICE))
2367 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2368 THTGS=T1*PI*EXP((3114.834/T1-0.278296)*QS*(1.+0.81*QS))
2370 ESLIQ=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
2371 ESICE=AICE*EXP((BICE*T1-CICE)/(T1-DICE))
2372 ES=(1.-RATIO2)*ESLIQ+RATIO2*ESICE
2374 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2375 THTGS=T1*PI*EXP(RL*QS*C5/T1*(1.+0.81*QS))
2378 IF(ABS(F1).LT.0.01)GOTO 50
2380 IF(ITCNT.GT.10)GOTO 50
2381 DT=F1*(T1-T0)/(F1-F0)
2387 ! IF THE PARCEL IS SUPERSATURATED, CALCULATE CONCENTRATION OF FRESH
2396 ! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
2397 ! ADJUSTED...IF LIQUID WATER OR ICE IS PRESENT, IT IS ALLOWED TO EVAPO
2404 ! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS
2405 ! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MI
2406 ! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURA
2407 ! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPR
2408 ! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE
2410 !...NOTE THAT THE LIQ AND ICE MAY BE PRESENT IN PROPORTIONS SLIGHTLY DIF
2411 ! THAN SUGGESTED BY THE VALUE OF RATIO2...CHECK TO MAKE SURE THAT LIQ
2412 ! ICE CONCENTRATIONS ARE NOT REDUCED TO BELOW ZERO WHEN EVAPORATION/
2413 ! SUBLIMATION OCCURS...
2418 QLIQ=QLIQ-(1.-RATIO2)*DQ
2423 QICE=QICE-RATIO2*DQ+DQICE
2432 IF(RATIO2.LT.1.E-6)THEN
2434 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
2439 CCP=1005.7*(1.+0.89*QU)
2440 IF(QTOT.LT.1.E-10)THEN
2442 !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
2443 T1=T1+RLL*(DQ/(1.+DQ))/CCP
2447 !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURA
2448 ! THE TEMPERATURE IS GIVEN BY:
2449 T1=T1+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CCP
2457 QNEWLQ=(1.-RATIO2)*QNEW
2459 IF(ITCNT.GT.10)PRINT*,'***** NUMBER OF ITERATIONS IN TPMIX =', &
2462 END SUBROUTINE TPMIX
2463 !-----------------------------------------------------------------------
2464 SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,R1,RL, &
2465 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2466 !-----------------------------------------------------------------------
2468 !-----------------------------------------------------------------------
2469 REAL, INTENT(IN ) :: P1,T1,Q1,R1,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,&
2471 REAL, INTENT(INOUT) :: THT1
2472 REAL:: T00,P00,C1,C2,C3,C4,C5,EE,TLOG,TDPT,TSAT,THT,TFPT,TLOGIC, &
2475 DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834,&
2478 ! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
2484 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2485 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
2486 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
2487 THT1=THT*EXP((C1/TSAT-C2)*Q1*(1.+0.81*Q1))
2488 ELSEIF(ABS(R1-1.).LT.1.E-6)THEN
2491 TFPT=(CICE-DICE*TLOG)/(BICE-TLOG)
2492 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
2493 TSAT=TFPT-(.182+1.13E-3*(TFPT-T00)-3.58E-4*(T1-T00))*(T1-TFPT)
2494 THT1=THT*EXP((C3/TSAT-C4)*Q1*(1.+0.81*Q1))
2498 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2499 TLOGIC=ALOG(EE/AICE)
2500 TFPT=(CICE-DICE*TLOGIC)/(BICE-TLOGIC)
2501 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
2502 TSATLQ=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
2503 TSATIC=TFPT-(.182+1.13E-3*(TFPT-T00)-3.58E-4*(T1-T00))*(T1-TFPT)
2504 TSAT=R1*TSATIC+(1.-R1)*TSATLQ
2505 THT1=THT*EXP(RL*Q1*C5/TSAT*(1.+0.81*Q1))
2508 END SUBROUTINE ENVIRTHT
2510 !-----------------------------------------------------------------------
2511 !************************* TPDD.FOR ************************************
2512 ! THIS SUBROUTINE ITERATIVELY EXTRACTS TEMPERATURE FROM EQUIVALENT *
2513 ! POTENTIAL TEMP. IT IS DESIGNED FOR USE WITH DOWNDRAFT CALCULATIONS.
2514 ! IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, PARCEL *
2515 ! TEMP, SPECIFIC HUMIDITY, AND LIQUID WATER CONTENT ARE ITERATIVELY *
2517 !***********************************************************************
2518 FUNCTION TPDD(P,THTED,TGS,RS,RD,RH,XLV0,XLV1, &
2519 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2520 !-----------------------------------------------------------------------
2522 !-----------------------------------------------------------------------
2523 REAL, INTENT(IN ) :: XLV0,XLV1
2524 REAL, INTENT(IN ) :: P,THTED,TGS,RD,RH,EP2,ALIQ,BLIQ, &
2525 CLIQ,DLIQ,AICE,BICE,CICE,DICE
2526 REAL, INTENT(INOUT) :: RS
2527 REAL :: TPDD,ES,PI,THTGS,F0,T1,T0,CCP,F1,DT,RL,DSSDT,T1RH,RSRH
2529 !-----------------------------------------------------------------------
2530 ES=ALIQ*EXP((BLIQ*TGS-CLIQ)/(TGS-DLIQ))
2532 PI=(1.E5/P)**(0.2854*(1.-0.28*RS))
2533 THTGS=TGS*PI*EXP((3374.6525/TGS-2.5403)*RS*(1.+0.81*RS))
2539 !...ITERATE TO FIND WET-BULB TEMPERATURE...
2542 90 ES=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
2544 PI=(1.E5/P)**(0.2854*(1.-0.28*RS))
2545 THTGS=T1*PI*EXP((3374.6525/T1-2.5403)*RS*(1.+0.81*RS))
2547 IF(ABS(F1).LT.0.05)GOTO 50
2549 IF(ITCNT.GT.10)GOTO 50
2550 DT=F1*(T1-T0)/(F1-F0)
2557 !...IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, ESTIMATE THE
2558 ! TEMPERATURE AND MIXING RATIO WHICH WILL YIELD THE APPROPRIATE VALUE.
2560 IF(RH.EQ.1.)GOTO 110
2561 DSSDT=(CLIQ-BLIQ*DLIQ)/((T1-DLIQ)*(T1-DLIQ))
2562 DT=RL*RS*(1.-RH)/(CCP+RL*RH*RS*DSSDT)
2564 ES=RH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
2567 !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
2568 !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
2572 T1RH=T1+(RS-RSRH)*RL/CCP
2577 IF(ITCNT.GT.10)PRINT*,'***** NUMBER OF ITERATIONS IN TPDD = ', &
2582 !====================================================================
2583 SUBROUTINE kfinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, &
2584 RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS, &
2585 P_FIRST_SCALAR,restart,allowed_to_read, &
2586 ids, ide, jds, jde, kds, kde, &
2587 ims, ime, jms, jme, kms, kme, &
2588 its, ite, jts, jte, kts, kte )
2589 !--------------------------------------------------------------------
2591 !--------------------------------------------------------------------
2592 LOGICAL , INTENT(IN) :: restart, allowed_to_read
2593 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
2594 ims, ime, jms, jme, kms, kme, &
2595 its, ite, jts, jte, kts, kte
2596 INTEGER , INTENT(IN) :: P_QI,P_QS,P_FIRST_SCALAR
2598 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
2606 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG
2608 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA
2610 INTEGER :: i, j, k, itf, jtf, ktf
2616 IF(.not.restart)THEN
2628 IF (P_QI .ge. P_FIRST_SCALAR) THEN
2638 IF (P_QS .ge. P_FIRST_SCALAR) THEN
2664 END SUBROUTINE kfinit
2666 !-------------------------------------------------------
2668 END MODULE module_cu_kf