/wrfv2_fire/phys/module_cu_kf.F
FORTRAN Legacy | 2669 lines | 1750 code | 115 blank | 804 comment | 3 complexity | 6487caca6f6ab8133f78922b148dcf06 MD5 | raw file
Possible License(s): AGPL-1.0
- !WRF:MODEL_LAYER:PHYSICS
- !
- MODULE module_cu_kf
- USE module_wrf_error
- REAL , PARAMETER :: RAD = 1500.
- CONTAINS
- !-------------------------------------------------------------
- SUBROUTINE KFCPS( &
- ids,ide, jds,jde, kds,kde &
- ,ims,ime, jms,jme, kms,kme &
- ,its,ite, jts,jte, kts,kte &
- ,DT,KTAU,DX,CUDT,CURR_SECS,ADAPT_STEP_FLAG &
- ,CUDTACTTIME &
- ,rho &
- ,RAINCV,PRATEC,NCA &
- ,U,V,TH,T,W,QV,dz8w,Pcps,pi &
- ,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
- ,EP2,SVP1,SVP2,SVP3,SVPT0 &
- ,STEPCU,CU_ACT_FLAG,warm_rain &
- ! optional arguments
- ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
- ,RQVCUTEN,RQCCUTEN,RQRCUTEN,RQICUTEN,RQSCUTEN &
- ,RTHCUTEN &
- )
- !-------------------------------------------------------------
- IMPLICIT NONE
- !-------------------------------------------------------------
- INTEGER, INTENT(IN ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- INTEGER, INTENT(IN ) :: STEPCU
- LOGICAL, INTENT(IN ) :: warm_rain
- REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1
- REAL, INTENT(IN ) :: CP,R,G,EP1,EP2
- REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
- INTEGER, INTENT(IN ) :: KTAU
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
- INTENT(IN ) :: &
- U, &
- V, &
- W, &
- TH, &
- QV, &
- T, &
- dz8w, &
- Pcps, &
- rho, &
- pi
- !
- REAL, INTENT(IN ) :: DT, DX
- REAL, INTENT(IN ) :: CUDT
- REAL, INTENT(IN ) :: CURR_SECS
- LOGICAL,OPTIONAL, INTENT(IN ) :: ADAPT_STEP_FLAG
- REAL, INTENT (INOUT) :: CUDTACTTIME
- REAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: &
- RAINCV &
- ,PRATEC &
- , NCA
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
- INTENT(INOUT) :: &
- W0AVG
- LOGICAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: CU_ACT_FLAG
- !
- ! Optional arguments
- !
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
- OPTIONAL, &
- INTENT(INOUT) :: &
- RTHCUTEN &
- ,RQVCUTEN &
- ,RQCCUTEN &
- ,RQRCUTEN &
- ,RQICUTEN &
- ,RQSCUTEN
- !
- ! Flags relating to the optional tendency arrays declared above
- ! Models that carry the optional tendencies will provdide the
- ! optional arguments at compile time; these flags all the model
- ! to determine at run-time whether a particular tracer is in
- ! use or not.
- !
- LOGICAL, OPTIONAL :: &
- F_QV &
- ,F_QC &
- ,F_QR &
- ,F_QI &
- ,F_QS
- ! LOCAL VARS
- REAL, DIMENSION( kts:kte ) :: &
- U1D, &
- V1D, &
- T1D, &
- DZ1D, &
- QV1D, &
- P1D, &
- RHO1D, &
- W0AVG1D
- REAL, DIMENSION( kts:kte ):: &
- DQDT, &
- DQIDT, &
- DQCDT, &
- DQRDT, &
- DQSDT, &
- DTDT
- REAL :: TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp
- INTEGER :: i,j,k,NTST,ICLDCK
- LOGICAL :: qi_flag , qs_flag
- ! adjustable time step changes
- REAL :: lastdt = -1.0
- REAL :: W0AVGfctr, W0fctr, W0den
- LOGICAL :: run_param , doing_adapt_dt , decided
- !----------------------------------------------------------------------
- !--- CALL CUMULUS PARAMETERIZATION
- !
- !...TST IS THE NUMBER OF TIME STEPS IN 10 MINUTES...W0AVG IS CLOSE TO A
- !...RUNNING MEAN VERTICAL VELOCITY...NOTE THAT IF YOU CHANGE TST, IT WIL
- !...CHANGE THE FREQUENCY OF THE CONVECTIVE INTITIATION CHECK (SEE BELOW)
- !...NOTE THAT THE ORDERING OF VERTICAL LAYERS MUST BE REVERSED FOR W0AVG
- !...BECAUSE THE ORDERING IS REVERSED IN KFPARA...
- !
- DXSQ=DX*DX
- qi_flag = .FALSE.
- qs_flag = .FALSE.
- IF ( PRESENT( F_QI ) ) qi_flag = f_qi
- IF ( PRESENT( F_QS ) ) qs_flag = f_qs
- !----------------------
- NTST=STEPCU
- TST=float(NTST*2)
- !----------------------
- ! NTST=NINT(1200./(DT*2.))
- ! TST=float(NTST)
- ! NTST=NINT(0.5*TST)
- ! NTST=MAX0(NTST,1)
- !----------------------
- ! ICLDCK=MOD(KTAU,NTST)
- !----------------------
- ! write(0,*) 'DT = ',DT,' KTAU = ',KTAU,' DX = ',DX
- ! write(0,*) 'CUDT = ',CUDT,' CURR_SECS = ',CURR_SECS
- ! write(0,*) 'ADAPT_STEP_FLAG = ',ADAPT_STEP_FLAG,' IDS = ',IDS
- ! write(0,*) 'STEPCU = ',STEPCU,' warm_rain = ',warm_rain
- ! write(0,*) 'F_QV = ',F_QV,' F_QC = ',F_QV
- ! write(0,*) 'F_QI = ',F_QI,' F_QS = ',F_QS
- ! write(0,*) 'F_QR = ',F_QR
- ! stop
- if (lastdt < 0) then
- lastdt = dt
- endif
- if (ADAPT_STEP_FLAG) then
- W0AVGfctr = 2 * MAX(CUDT*60,dt) - dt
- W0fctr = dt
- W0den = 2 * MAX(CUDT*60,dt)
- else
- W0AVGfctr = (TST-1.)
- W0fctr = 1.
- W0den = TST
- endif
- DO J = jts,jte
- DO K=kts,kte
- DO I= its,ite
- ! SCR1=-5.0E-4*G*rho(I,K,J)*(w(I,K,J)+w(I,K+1,J))
- ! TV=T(I,K,J)*(1.+EP1*QV(I,K,J))
- ! RHOE=Pcps(I,K,J)/(R*TV)
- ! W0=-101.9368*SCR1/RHOE
- W0=0.5*(w(I,K,J)+w(I,K+1,J))
- ! Old:
- !
- ! W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST
- ! New, to support adaptive time step:
- !
- W0AVG(I,K,J) = ( W0AVG(I,K,J) * W0AVGfctr + W0 * W0fctr ) / W0den
- ENDDO
- ENDDO
- ENDDO
- lastdt = dt
- ! Initialization for adaptive time step.
- doing_adapt_dt = .FALSE.
- IF ( PRESENT(adapt_step_flag) ) THEN
- IF ( adapt_step_flag ) THEN
- doing_adapt_dt = .TRUE.
- IF ( cudtacttime .EQ. 0. ) THEN
- cudtacttime = curr_secs + cudt*60.
- END IF
- END IF
- END IF
- ! Do we run through this scheme or not?
- ! Test 1: If this is the initial model time, then yes.
- ! KTAU=1
- ! Test 2: If the user asked for the cumulus to be run every time step, then yes.
- ! CUDT=0 or STEPCU=1
- ! Test 3: If not adaptive dt, and this is on the requested cumulus frequency, then yes.
- ! MOD(KTAU,NST)=0
- ! Test 4: If using adaptive dt and the current time is past the last requested activate cumulus time, then yes.
- ! CURR_SECS >= CUDTACTTIME
- ! If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag
- ! to TRUE. The decided flag says that one of these tests was able to say "yes", run the scheme.
- ! We only proceed to other tests if the previous tests all have left decided as FALSE.
- ! If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next
- ! cumulus run.
- decided = .FALSE.
- run_param = .FALSE.
- IF ( ( .NOT. decided ) .AND. &
- ( ktau .EQ. 1 ) ) THEN
- run_param = .TRUE.
- decided = .TRUE.
- END IF
- IF ( ( .NOT. decided ) .AND. &
- ( ( cudt .EQ. 0. ) .OR. ( stepcu .EQ. 1 ) ) ) THEN
- run_param = .TRUE.
- decided = .TRUE.
- END IF
- IF ( ( .NOT. decided ) .AND. &
- ( .NOT. doing_adapt_dt ) .AND. &
- ( MOD(ktau,ntst) .EQ. 0 ) ) THEN
- run_param = .TRUE.
- decided = .TRUE.
- END IF
- IF ( ( .NOT. decided ) .AND. &
- ( doing_adapt_dt ) .AND. &
- ( curr_secs .GE. cudtacttime ) ) THEN
- run_param = .TRUE.
- decided = .TRUE.
- cudtacttime = curr_secs + cudt*60
- END IF
- IF (run_param) then
- DO J = jts,jte
- DO I= its,ite
- CU_ACT_FLAG(i,j) = .true.
- ENDDO
- ENDDO
- DO J = jts,jte
- DO I=its,ite
- ! if (i.eq. 110 .and. j .eq. 59 ) then
- ! write(0,*) 'nca = ',nca(i,j),' CU_ACT_FLAG = ',CU_ACT_FLAG(i,j)
- ! write(0,*) 'dt = ',dt,' ADAPT_STEP_FLAG = ',ADAPT_STEP_FLAG
- ! endif
- ! IF ( NINT(NCA(I,J)) .gt. 0 ) then
- IF ( NCA(I,J) .gt. 0.5*DT ) then
- CU_ACT_FLAG(i,j) = .false.
- ELSE
- DO k=kts,kte
- DQDT(k)=0.
- DQIDT(k)=0.
- DQCDT(k)=0.
- DQRDT(k)=0.
- DQSDT(k)=0.
- DTDT(k)=0.
- ENDDO
- RAINCV(I,J)=0.
- PRATEC(I,J)=0.
- !
- ! assign vars from 3D to 1D
- DO K=kts,kte
- U1D(K) =U(I,K,J)
- V1D(K) =V(I,K,J)
- T1D(K) =T(I,K,J)
- RHO1D(K) =rho(I,K,J)
- QV1D(K)=QV(I,K,J)
- P1D(K) =Pcps(I,K,J)
- W0AVG1D(K) =W0AVG(I,K,J)
- DZ1D(k)=dz8w(I,K,J)
- ENDDO
- !
- CALL KFPARA(I, J, &
- U1D,V1D,T1D,QV1D,P1D,DZ1D, &
- W0AVG1D,DT,DX,DXSQ,RHO1D, &
- XLV0,XLV1,XLS0,XLS1,CP,R,G, &
- EP2,SVP1,SVP2,SVP3,SVPT0, &
- DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
- RAINCV,PRATEC,NCA, &
- warm_rain,qi_flag,qs_flag, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
-
- IF ( PRESENT( RTHCUTEN ) .AND. PRESENT( RQVCUTEN ) ) THEN
- DO K=kts,kte
- RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
- RQVCUTEN(I,K,J)=DQDT(K)
- ENDDO
- ENDIF
- IF( PRESENT(RQRCUTEN) .AND. PRESENT(RQCCUTEN) .AND. &
- PRESENT(F_QR) ) THEN
- IF ( F_QR ) THEN
- DO K=kts,kte
- RQRCUTEN(I,K,J)=DQRDT(K)
- RQCCUTEN(I,K,J)=DQCDT(K)
- ENDDO
- ELSE
- ! This is the case for Eta microphysics without 3d rain field
- DO K=kts,kte
- RQRCUTEN(I,K,J)=0.
- RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
- ENDDO
- ENDIF
- ENDIF
- !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
- IF( PRESENT( RQICUTEN ) .AND. qi_flag )THEN
- DO K=kts,kte
- RQICUTEN(I,K,J)=DQIDT(K)
- ENDDO
- ENDIF
- IF( PRESENT ( RQSCUTEN ) .AND. qs_flag )THEN
- DO K=kts,kte
- RQSCUTEN(I,K,J)=DQSDT(K)
- ENDDO
- ENDIF
- !
- ENDIF
- ENDDO
- ENDDO
- ENDIF
- END SUBROUTINE KFCPS
-
- !-----------------------------------------------------------
- SUBROUTINE KFPARA (I, J, &
- U0,V0,T0,QV0,P0,DZQ,W0AVG1D, &
- DT,DX,DXSQ,rho, &
- XLV0,XLV1,XLS0,XLS1,CP,R,G, &
- EP2,SVP1,SVP2,SVP3,SVPT0, &
- DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
- RAINCV,PRATEC,NCA, &
- warm_rain,qi_flag,qs_flag, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- !-----------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------
- INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte, &
- I,J
- LOGICAL, INTENT(IN ) :: warm_rain
- LOGICAL :: qi_flag, qs_flag
- !
- REAL, DIMENSION( kts:kte ), &
- INTENT(IN ) :: U0, &
- V0, &
- T0, &
- QV0, &
- P0, &
- rho, &
- DZQ, &
- W0AVG1D
- !
- REAL, INTENT(IN ) :: DT,DX,DXSQ
- !
- REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
- REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
- !
- REAL, DIMENSION( kts:kte ), INTENT(INOUT) :: &
- DQDT, &
- DQIDT, &
- DQCDT, &
- DQRDT, &
- DQSDT, &
- DTDT
- REAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: RAINCV, &
- PRATEC, &
- NCA
- !
- !...DEFINE LOCAL VARIABLES...
- !
- REAL, DIMENSION( kts:kte ) :: &
- Q0,Z0,TV0,TU,TVU,QU,TZ,TVD, &
- QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD, &
- UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2, &
- UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE, &
- THTAU,THETEU,THTAD,THETED,QLIQ,QICE, &
- QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC, &
- DETLQ2,DETIC2,RATIO,RATIO2
- REAL, DIMENSION( kts:kte ) :: &
- DOMGDP,EXN,RHOE,TVQU,DP,RH,EQFRC,WSPD, &
- QDT,FXM,THTAG,THTESG,THPA,THFXTOP, &
- THFXBOT,QPA,QFXTOP,QFXBOT,QLPA,QLFXIN, &
- QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA, &
- QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT, &
- QL0,QLG,QI0,QIG,QR0,QRG,QS0,QSG
-
- REAL, DIMENSION( kts:kte+1 ) :: OMG
- REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
- ! LOCAL VARS
- REAL :: P00,T00,CV,B61,RLF,RHIC,RHBC,PIE, &
- TTFRZ,TBFRZ,C5,RATE
- REAL :: GDRY,ROCP,ALIQ,BLIQ, &
- CLIQ,DLIQ,AICE,BICE,CICE,DICE
- REAL :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX, &
- ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL, &
- CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR, &
- ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
- TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD, &
- UPNEW,ABE,WKLCL,THTUDL,TUDL,TTEMP,FRC1, &
- QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
- DZZ,WSQ,UDLBE,REI,EE2,UD2,TTMP,F1,F2, &
- THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1, &
- UD1,CLDHGT,DPTT,QNEWLQ,DUMFDP,EE,TSAT, &
- THTA,P150,USR,VCONV,TIMEC,SHSIGN,VWS,PEF, &
- CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN, &
- DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1, &
- DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF, &
- UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF, &
- DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
- AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1, &
- DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF, &
- TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR, &
- UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2, &
- RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
- DDFRC,TDC,DEFRC
- INTEGER :: KX,K,KL
- !
- INTEGER :: ISTOP,ML,L5,L4,KMIX,LOW, &
- LC,MXLAYR,LLFC,NLAYRS,NK, &
- KPBL,KLCL,LCL,LET,IFLAG, &
- KFRZ,NK1,LTOP,NJ,LTOP1, &
- LTOPM1,LVF,KSTART,KMIN,LFS, &
- ND,NIC,LDB,LDT,ND1,NDK, &
- NM,LMAX,NCOUNT,NOITR, &
- NSTEP,NTC
- !
- DATA P00,T00/1.E5,273.16/
- DATA CV,B61,RLF/717.,0.608,3.339E5/
- DATA RHIC,RHBC/1.,0.90/
- DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
- DATA RATE/0.01/
- !-----------------------------------------------------------
- GDRY=-G/CP
- ROCP=R/CP
- KL=kte
- KX=kte
- !
- ! ALIQ = 613.3
- ! BLIQ = 17.502
- ! CLIQ = 4780.8
- ! DLIQ = 32.19
- ALIQ = SVP1*1000.
- BLIQ = SVP2
- CLIQ = SVP2*SVPT0
- DLIQ = SVP3
- AICE = 613.2
- BICE = 22.452
- CICE = 6133.0
- DICE = 0.61
- !
- !...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER
- !...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL)
- !...FIELD. 'FBFRC' IS THE FRACTION OF AVAILABLE
- !...PRECIPITATION TO BE FED BACK (0.0 - 1.0)...
- !
- FBFRC=0.0
- !
- !...SCHEME IS CALLED ONCE ON EACH NORTH-SOUTH SLICE, THE LOOP BELOW
- !...CHECKS FOR THE POSSIBILITY OF INITIATING PARAMETERIZED
- !...CONVECTION AT EACH POINT WITHIN THE SLICE
- !
- !...SEE IF IT IS NECESSARY TO CHECK FOR CONVECTIVE TRIGGERING AT THIS
- !...GRID POINT. IF NCA>0, CONVECTION IS ALREADY ACTIVE AT THIS POINT,
- !...JUST FEED BACK THE TENDENCIES SAVED FROM THE TIME WHEN CONVECTION
- !...WAS INITIATED. IF NCA<0, CONVECTION IS NOT ACTIVE
- !...AND YOU MAY WANT TO CHECK TO SEE IF IT CAN BE ACTIVATED FOR THE
- !...CURRENT CONDITIONS. IN PREVIOUS APLICATIONS OF THIS SCHEME,
- !...THE VARIABLE ICLDCK WAS USED BELOW TO SAVE TIME BY ONLY CHECKING
- !...FOR THE POSSIBILITY OF CONVECTIVE INITIATION EVERY 5 OR 10
- !...MINUTES...
- !
- ! 10 CONTINUE
- !SUE P300=1000.*(PSB(I,J)*A(KL)+PTOP-30.)+PP3D(I,J,KL)
- P300=P0(1)-30000.
- !
- !...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
- !...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
- !...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
- !
- !...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
- !...FROM BOTTOM-UP IN THE KF SCHEME...
- !
- ML=0
- !SUE tmprpsb=1./PSB(I,J)
- !SUE CELL=PTOP*tmprpsb
- DO 15 K=1,KX
- !SUE P0(K)=1.E3*(A(NK)*PSB(I,J)+PTOP)+PP3D(I,J,NK)
- !
- !...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
- !
- ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
- QES(K)=EP2*ES/(P0(K)-ES)
- Q0(K)=AMIN1(QES(K),QV0(K))
- Q0(K)=AMAX1(0.000001,Q0(K))
- QL0(K)=0.
- QI0(K)=0.
- QR0(K)=0.
- QS0(K)=0.
- TV0(K)=T0(K)*(1.+B61*Q0(K))
- RHOE(K)=P0(K)/(R*TV0(K))
- DP(K)=rho(k)*g*DZQ(k)
- !
- !...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
- ! DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
- !
- IF(P0(K).GE.500E2)L5=K
- IF(P0(K).GE.400E2)L4=K
- IF(P0(K).GE.P300)LLFC=K
- IF(T0(K).GT.T00)ML=K
- 15 CONTINUE
- Z0(1)=.5*DZQ(1)
- DO 20 K=2,KL
- Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
- DZA(K-1)=Z0(K)-Z0(K-1)
- 20 CONTINUE
- DZA(KL)=0.
- KMIX=1
- 25 LOW=KMIX
- IF(LOW.GT.LLFC)GOTO 325
- LC=LOW
- MXLAYR=0
- !
- !...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
- !...UNSTABLE AIR 50 TO 100 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
- !...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
- !...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 60 mb..
- !
- NLAYRS=0
- DPTHMX=0.
- DO 63 NK=LC,KX
- DPTHMX=DPTHMX+DP(NK)
- NLAYRS=NLAYRS+1
- 63 IF(DPTHMX.GT.6.E3)GOTO 64
- GOTO 325
- 64 KPBL=LC+NLAYRS-1
- KMIX=LC+1
- 18 THMIX=0.
- QMIX=0.
- ZMIX=0.
- PMIX=0.
- DPTHMX=0.
- !
- !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
- !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
- !...LAYERS...
- !
- DO 17 NK=LC,KPBL
- DPTHMX=DPTHMX+DP(NK)
- ROCPQ=0.2854*(1.-0.28*Q0(NK))
- THMIX=THMIX+DP(NK)*T0(NK)*(P00/P0(NK))**ROCPQ
- QMIX=QMIX+DP(NK)*Q0(NK)
- ZMIX=ZMIX+DP(NK)*Z0(NK)
- 17 PMIX=PMIX+DP(NK)*P0(NK)
- THMIX=THMIX/DPTHMX
- QMIX=QMIX/DPTHMX
- ZMIX=ZMIX/DPTHMX
- PMIX=PMIX/DPTHMX
- ROCPQ=0.2854*(1.-0.28*QMIX)
- TMIX=THMIX*(PMIX/P00)**ROCPQ
- EMIX=QMIX*PMIX/(EP2+QMIX)
- !
- !...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL, PRESSURE
- !...LEVEL OF LCL...
- !
- TLOG=ALOG(EMIX/ALIQ)
- TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
- TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX- &
- TDPT)
- TLCL=AMIN1(TLCL,TMIX)
- TVLCL=TLCL*(1.+0.608*QMIX)
- CPORQ=1./ROCPQ
- PLCL=P00*(TLCL/THMIX)**CPORQ
- DO 29 NK=LC,KL
- KLCL=NK
- IF(PLCL.GE.P0(NK))GOTO 35
- 29 CONTINUE
- GOTO 325
- 35 K=KLCL-1
- DLP=ALOG(PLCL/P0(K))/ALOG(P0(KLCL)/P0(K))
- !
- !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
- !
- TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
- QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
- TVEN=TENV*(1.+0.608*QENV)
- TVBAR=0.5*(TV0(K)+TVEN)
- ! ZLCL=Z0(K)+R*TVBAR*ALOG(P0(K)/PLCL)/G
- ZLCL=Z0(K)+(Z0(KLCL)-Z0(K))*DLP
- !
- !...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
- !...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0AVG IS AN
- !...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
- !...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
- !...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
- !...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
- !...SUCCESS AT GRID LENGTHS NEAR 25 km. FOR DIFFERENT GRID-LENGTHS,
- !...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
- !...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
- !
- WKLCL=0.02*ZLCL/2.5E3
- WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)*DX/25.E3- &
- WKLCL
- WABS=ABS(WKL)+1.E-10
- WSIGNE=WKL/WABS
- DTLCL=4.64*WSIGNE*WABS**0.33
- GDT=G*DTLCL*(ZLCL-Z0(LC))/(TV0(LC)+TVEN)
- WLCL=1.+.5*WSIGNE*SQRT(ABS(GDT)+1.E-10)
- IF(TLCL+DTLCL.GT.TENV)GOTO 45
- IF(KPBL.GE.LLFC)GOTO 325
- GOTO 25
- !
- !...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
- !...EQUIVALENT POTENTIAL TEMPERATURE
- !...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
- !
- 45 THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
- EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
- ES=ALIQ*EXP((TENV*BLIQ-CLIQ)/(TENV-DLIQ))
- TVAVG=0.5*(TV0(KLCL)+TENV*(1.+0.608*QENV))
- PLCL=P0(KLCL)*EXP(G/(R*TVAVG)*(Z0(KLCL)-ZLCL))
- QESE=EP2*ES/(PLCL-ES)
- GDT=G*DTLCL*(ZLCL-Z0(LC))/(TV0(LC)+TVEN)
- WLCL=1.+.5*WSIGNE*SQRT(ABS(GDT)+1.E-10)
- THTES(K)=TENV*(1.E5/PLCL)**(0.2854*(1.-0.28*QESE))* &
- EXP((3374.6525/TENV-2.5403)*QESE*(1.+0.81*QESE))
- WTW=WLCL*WLCL
- IF(WLCL.LT.0.)GOTO 25
- TVLCL=TLCL*(1.+0.608*QMIX)
- RHOLCL=PLCL/(R*TVLCL)
- !
- LCL=KLCL
- LET=LCL
- !
- !*******************************************************************
- ! *
- ! COMPUTE UPDRAFT PROPERTIES *
- ! *
- !*******************************************************************
- !
- !
- !...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
- !
- WU(K)=WLCL
- AU0=PIE*RAD*RAD
- UMF(K)=RHOLCL*AU0
- VMFLCL=UMF(K)
- UPOLD=VMFLCL
- UPNEW=UPOLD
- !
- !...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
- !...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE BUOYANT ENERGY,
- ! TRPPT IS THE TOTAL RATE OF PRECIPITATION PRODUCTION...
- !
- RATIO2(K)=0.
- UER(K)=0.
- ABE=0.
- TRPPT=0.
- TU(K)=TLCL
- TVU(K)=TVLCL
- QU(K)=QMIX
- EQFRC(K)=1.
- QLIQ(K)=0.
- QICE(K)=0.
- QLQOUT(K)=0.
- QICOUT(K)=0.
- DETLQ(K)=0.
- DETIC(K)=0.
- PPTLIQ(K)=0.
- PPTICE(K)=0.
- IFLAG=0
- KFRZ=LC
- !
- !...THE AMOUNT OF CONV AVAIL POT ENERGY (CAPE) IS CALCULATED WITH
- ! RESPECT TO UNDILUTE PARCEL ASCENT; EQ POT TEMP OF UNDILUTE
- ! PARCEL IS THTUDL, UNDILUTE TEMPERATURE IS GIVEN BY TUDL...
- !
- THTUDL=THETEU(K)
- TUDL=TLCL
- !
- !...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
- ! PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
- ! FREEZING IS SPECIFIED TO BEGIN. WITHIN THE GLACIATION
- ! INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
- ! PREVIOUS MODEL LEVEL...
- !
- TTEMP=TTFRZ
- !
- !...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
- ! MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
- ! MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
- !
- DO 60 NK=K,KL-1
- NK1=NK+1
- RATIO2(NK1)=RATIO2(NK)
- !
- !...UPDATE UPDRAFT PROPERTIES AT THE NEXT MODEL LVL TO REFLECT
- ! ENTRAINMENT OF ENVIRONMENTAL AIR...
- !
- FRC1=0.
- TU(NK1)=T0(NK1)
- THETEU(NK1)=THETEU(NK)
- QU(NK1)=QU(NK)
- QLIQ(NK1)=QLIQ(NK)
- QICE(NK1)=QICE(NK)
- CALL TPMIX(P0(NK1),THETEU(NK1),TU(NK1),QU(NK1),QLIQ(NK1), &
- QICE(NK1),QNEWLQ,QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0, &
- XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
- TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
- !
- !...CHECK TO SEE IF UPDRAFT TEMP IS WITHIN THE FREEZING INTERVAL,
- ! IF IT IS, CALCULATE THE FRACTIONAL CONVERSION TO GLACIATION
- ! AND ADJUST QNEWLQ TO REFLECT THE GRADUAL CHANGE IN THETAU
- ! SINCE THE LAST MODEL LEVEL...THE GLACIATION EFFECTS WILL BE
- ! DETERMINED AFTER THE AMOUNT OF CONDENSATE AVAILABLE AFTER
- ! PRECIP FALLOUT IS DETERMINED...TTFRZ IS THE TEMP AT WHICH
- ! GLACIATION BEGINS, TBFRZ THE TEMP AT WHICH IT ENDS...
- !
- IF(TU(NK1).LE.TTFRZ.AND.IFLAG.LT.1)THEN
- IF(TU(NK1).GT.TBFRZ)THEN
- IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
- FRC1=(TTEMP-TU(NK1))/(TTFRZ-TBFRZ)
- R1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
- ELSE
- FRC1=(TTEMP-TBFRZ)/(TTFRZ-TBFRZ)
- R1=1.
- IFLAG=1
- ENDIF
- QNWFRZ=QNEWLQ
- QNEWIC=QNEWIC+QNEWLQ*R1*0.5
- QNEWLQ=QNEWLQ-QNEWLQ*R1*0.5
- EFFQ=(TTFRZ-TBFRZ)/(TTEMP-TBFRZ)
- TTEMP=TU(NK1)
- ENDIF
- !
- ! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
- !
- IF(NK.EQ.K)THEN
- BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
- BOTERM=2.*(Z0(NK1)-ZLCL)*G*BE/1.5
- ENTERM=0.
- DZZ=Z0(NK1)-ZLCL
- ELSE
- BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
- BOTERM=2.*DZA(NK)*G*BE/1.5
- ENTERM=2.*UER(NK)*WTW/UPOLD
- DZZ=DZA(NK)
- ENDIF
- WSQ=WTW
- CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM,RATE, &
- QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1), G)
-
- !...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
- ! IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
- !
- IF(WTW.LE.0.)GOTO 65
- WABS=SQRT(ABS(WTW))
- WU(NK1)=WTW/WABS
- !
- ! UPDATE THE ABE FOR UNDILUTE ASCENT...
- !
- THTES(NK1)=T0(NK1)*(1.E5/P0(NK1))**(0.2854*(1.-0.28*QES(NK1))) &
- * &
- EXP((3374.6525/T0(NK1)-2.5403)*QES(NK1)*(1.+0.81* &
- QES(NK1)))
- UDLBE=((2.*THTUDL)/(THTES(NK)+THTES(NK1))-1.)*DZZ
- IF(UDLBE.GT.0.)ABE=ABE+UDLBE*G
- !
- ! DETERMINE THE EFFECTS OF CLOUD GLACIATION IF WITHIN THE SPECIFIED
- ! TEMP INTERVAL...
- !
- IF(FRC1.GT.1.E-6)THEN
- CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QLIQ(NK1), &
- QICE(NK1),RATIO2(NK1),TTFRZ,TBFRZ,QNWFRZ,RL,FRC1,EFFQ, &
- IFLAG,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE &
- ,CICE,DICE)
- ENDIF
- !
- ! CALL SUBROUTINE TO CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMP.
- ! WITHIN GLACIATION INTERVAL, THETAE MUST BE CALCULATED WITH RESPECT TO
- ! SAME DEGREE OF GLACIATION FOR ALL ENTRAINING AIR...
- !
- CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),RATIO2(NK1), &
- RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
-
- !...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
- !
- REI=VMFLCL*DP(NK1)*0.03/RAD
- TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
- !
- !...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, NO
- ! ENTRAINMENT IS ALLOWED AT THIS LEVEL...
- !
- IF(TVQU(NK1).LE.TV0(NK1))THEN
- UER(NK1)=0.0
- UDR(NK1)=REI
- EE2=0.
- UD2=1.
- EQFRC(NK1)=0.
- GOTO 55
- ENDIF
- LET=NK1
- TTMP=TVQU(NK1)
- !
- !...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL
- ! AIR FOR ESTIMATION OF ENTRAINMENT AND DETRAINMENT RATES...
- !
- F1=0.95
- F2=1.-F1
- THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
- QTMP=F1*Q0(NK1)+F2*QU(NK1)
- TMPLIQ=F2*QLIQ(NK1)
- TMPICE=F2*QICE(NK1)
- CALL TPMIX(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,QNEWLQ, &
- QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ, &
- DLIQ,AICE,BICE,CICE,DICE)
- TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
- IF(TU95.GT.TV0(NK1))THEN
- EE2=1.
- UD2=0.
- EQFRC(NK1)=1.0
- GOTO 50
- ENDIF
- F1=0.10
- F2=1.-F1
- THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
- QTMP=F1*Q0(NK1)+F2*QU(NK1)
- TMPLIQ=F2*QLIQ(NK1)
- TMPICE=F2*QICE(NK1)
- CALL TPMIX(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,QNEWLQ, &
- QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ, &
- DLIQ,AICE,BICE,CICE,DICE)
- TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
- IF(TU10.EQ.TVQU(NK1))THEN
- EE2=1.
- UD2=0.
- EQFRC(NK1)=1.0
- GOTO 50
- ENDIF
- EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
- EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
- EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
- IF(EQFRC(NK1).EQ.1)THEN
- EE2=1.
- UD2=0.
- GOTO 50
- ELSEIF(EQFRC(NK1).EQ.0.)THEN
- EE2=0.
- UD2=1.
- GOTO 50
- ELSE
- !
- !...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
- ! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
- !
- CALL PROF5(EQFRC(NK1),EE2,UD2)
- ENDIF
- !
- 50 IF(NK.EQ.K)THEN
- EE1=1.
- UD1=0.
- ENDIF
- !
- !...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE
- ! FRACTIONAL VALUES IN THE LAYER...
- !
- UER(NK1)=0.5*REI*(EE1+EE2)
- UDR(NK1)=0.5*REI*(UD1+UD2)
- !
- !...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
- ! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATION
- !
- 55 IF(UMF(NK)-UDR(NK1).LT.10.)THEN
- !
- !...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL
- ! UPDRAFT FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE
- ! PREVIOUS MODEL
- !
- IF(UDLBE.GT.0.)ABE=ABE-UDLBE*G
- LET=NK
- ! WRITE(98,1015)P0(NK1)/100.
- GOTO 65
- ENDIF
- EE1=EE2
- UD1=UD2
- UPOLD=UMF(NK)-UDR(NK1)
- UPNEW=UPOLD+UER(NK1)
- UMF(NK1)=UPNEW
- !
- !...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND ICE IN
- ! THE DETRAINING UPDRAFT MASS...
- !
- DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
- DETIC(NK1)=QICE(NK1)*UDR(NK1)
- QDT(NK1)=QU(NK1)
- QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
- THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
- QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
- QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
- !
- !...KFRZ IS THE HIGHEST MODEL LEVEL AT WHICH LIQUID CONDENSATE IS
- ! GENERATING PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF LIQUID
- ! PRECIP AT A GIVING MODEL LVL, PPTICE THE SAME FOR ICE, TRPPT IS
- ! THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE CURRENT MODEL LEVEL
- !
- IF(ABS(RATIO2(NK1)-1.).GT.1.E-6)KFRZ=NK1
- PPTLIQ(NK1)=QLQOUT(NK1)*(UMF(NK)-UDR(NK1))
- PPTICE(NK1)=QICOUT(NK1)*(UMF(NK)-UDR(NK1))
- TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
- IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
- 60 CONTINUE
- !
- !...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIU
- ! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
- ! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE
- ! BETWEEN THE LET AND CLOUD TOP...
- !
- !...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL
- ! VELOCITY FIRST BECOMES NEGATIVE...
- !
- 65 LTOP=NK
- CLDHGT=Z0(LTOP)-ZLCL
- !
- !...IF CLOUD TOP HGT IS LESS THAN SPECIFIED MINIMUM HEIGHT, GO BACK AND
- ! THE NEXT HIGHEST 60MB LAYER TO SEE IF A BIGGER CLOUD CAN BE OBTAINED
- ! THAT SOURCE AIR...
- !
- ! IF(CLDHGT.LT.4.E3.OR.ABE.LT.1.)THEN
- IF(CLDHGT.LT.3.E3.OR.ABE.LT.1.)THEN
- DO 70 NK=K,LTOP
- UMF(NK)=0.
- UDR(NK)=0.
- UER(NK)=0.
- DETLQ(NK)=0.
- DETIC(NK)=0.
- PPTLIQ(NK)=0.
- 70 PPTICE(NK)=0.
- GOTO 25
- ENDIF
- !
- !...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS
- ! FLUX THIS LEVEL...
- !
- IF(LET.EQ.LTOP)THEN
- UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
- DETLQ(LTOP)=QLIQ(LTOP)*UDR(LTOP)*UPNEW/UPOLD
- DETIC(LTOP)=QICE(LTOP)*UDR(LTOP)*UPNEW/UPOLD
- TRPPT=TRPPT-(PPTLIQ(LTOP)+PPTICE(LTOP))
- UER(LTOP)=0.
- UMF(LTOP)=0.
- GOTO 85
- ENDIF
- !
- ! BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
- !
- DPTT=0.
- DO 71 NJ=LET+1,LTOP
- 71 DPTT=DPTT+DP(NJ)
- DUMFDP=UMF(LET)/DPTT
- !
- !...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
- ! RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
- ! PTOP
- !
- DO 75 NK=LET+1,LTOP
- UDR(NK)=DP(NK)*DUMFDP
- UMF(NK)=UMF(NK-1)-UDR(NK)
- DETLQ(NK)=QLIQ(NK)*UDR(NK)
- DETIC(NK)=QICE(NK)*UDR(NK)
- TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
- PPTLIQ(NK)=(UMF(NK-1)-UDR(NK))*QLQOUT(NK)
- PPTICE(NK)=(UMF(NK-1)-UDR(NK))*QICOUT(NK)
- TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
- 75 CONTINUE
- !
- !...SEND UPDRAFT CHARACTERISTICS TO OUTPUT FILES...
- !
- 85 CONTINUE
- !
- !...EXTEND THE UPDRAFT MASS FLUX PROFILE DOWN TO THE SOURCE LAYER FOR
- ! THE UPDRAFT AIR...ALSO, DEFINE THETAE FOR LEVELS BELOW THE LCL...
- !
- DO 90 NK=1,K
- IF(NK.GE.LC)THEN
- IF(NK.EQ.LC)THEN
- UMF(NK)=VMFLCL*DP(NK)/DPTHMX
- UER(NK)=VMFLCL*DP(NK)/DPTHMX
- ELSEIF(NK.LE.KPBL)THEN
- UER(NK)=VMFLCL*DP(NK)/DPTHMX
- UMF(NK)=UMF(NK-1)+UER(NK)
- ELSE
- UMF(NK)=VMFLCL
- UER(NK)=0.
- ENDIF
- TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
- QU(NK)=QMIX
- WU(NK)=WLCL
- ELSE
- TU(NK)=0.
- QU(NK)=0.
- UMF(NK)=0.
- WU(NK)=0.
- UER(NK)=0.
- ENDIF
- UDR(NK)=0.
- QDT(NK)=0.
- QLIQ(NK)=0.
- QICE(NK)=0.
- QLQOUT(NK)=0.
- QICOUT(NK)=0.
- PPTLIQ(NK)=0.
- PPTICE(NK)=0.
- DETLQ(NK)=0.
- DETIC(NK)=0.
- RATIO2(NK)=0.
- EE=Q0(NK)*P0(NK)/(EP2+Q0(NK))
- TLOG=ALOG(EE/ALIQ)
- TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
- TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T0(NK)-T00))*( &
- T0(NK)-TDPT)
- THTA=T0(NK)*(1.E5/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
- THETEE(NK)=THTA* &
- EXP((3374.6525/TSAT-2.5403)*Q0(NK)*(1.+0.81*Q0(NK)) &
- )
- THTES(NK)=THTA* &
- EXP((3374.6525/T0(NK)-2.5403)*QES(NK)*(1.+0.81* &
- QES(NK)))
- EQFRC(NK)=1.0
- 90 CONTINUE
- !
- LTOP1=LTOP+1
- LTOPM1=LTOP-1
- !
- !...DEFINE VARIABLES ABOVE CLOUD TOP...
- !
- DO 95 NK=LTOP1,KX
- UMF(NK)=0.
- UDR(NK)=0.
- UER(NK)=0.
- QDT(NK)=0.
- QLIQ(NK)=0.
- QICE(NK)=0.
- QLQOUT(NK)=0.
- QICOUT(NK)=0.
- DETLQ(NK)=0.
- DETIC(NK)=0.
- PPTLIQ(NK)=0.
- PPTICE(NK)=0.
- IF(NK.GT.LTOP1)THEN
- TU(NK)=0.
- QU(NK)=0.
- WU(NK)=0.
- ENDIF
- THTA0(NK)=0.
- THTAU(NK)=0.
- EMS(NK)=DP(NK)*DXSQ/G
- EMSD(NK)=1./EMS(NK)
- TG(NK)=T0(NK)
- QG(NK)=Q0(NK)
- QLG(NK)=0.
- QIG(NK)=0.
- QRG(NK)=0.
- QSG(NK)=0.
- 95 OMG(NK)=0.
- OMG(KL+1)=0.
- P150=P0(KLCL)-1.50E4
- DO 100 NK=1,LTOP
- THTAD(NK)=0.
- EMS(NK)=DP(NK)*DXSQ/G
- EMSD(NK)=1./EMS(NK)
- !
- !...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION
- ! SCHEME
- !
- EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
- THTAU(NK)=TU(NK)*EXN(NK)
- EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
- THTA0(NK)=T0(NK)*EXN(NK)
- !
- !...LVF IS THE LEVEL AT WHICH MOISTURE FLUX IS ESTIMATED AS THE BASIS
- !...FOR PRECIPITATION EFFICIENCY CALCULATIONS...
- !
- IF(P0(NK).GT.P150)LVF=NK
- 100 OMG(NK)=0.
- LVF=MIN0(LVF,LET)
- USR=UMF(LVF+1)*(QU(LVF+1)+QLIQ(LVF+1)+QICE(LVF+1))
- USR=AMIN1(USR,TRPPT)
- IF(USR.LT.1.E-8)USR=TRPPT
- !
- ! WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
- ! * TMIX-T00,PMIX,QMIX,ABE
- ! WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
- ! * WLCL,CLDHGT
- !
- !...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
- !...AND MIDTROPOSPHERE IS USED.
- !
- WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
- WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
- WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
- VCONV=.5*(WSPD(KLCL)+WSPD(L5))
- if (VCONV .gt. 0.) then
- TIMEC=DX/VCONV
- else
- TIMEC=3600.
- endif
- ! TIMEC=DX/VCONV
- TADVEC=TIMEC
- TIMEC=AMAX1(1800.,TIMEC)
- TIMEC=AMIN1(3600.,TIMEC)
- NIC=NINT(TIMEC/DT)
- TIMEC=FLOAT(NIC)*DT
- !
- !...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
- !
- ! SHSIGN = CVMGT(1.,-1.,WSPD(LTOP).GT.WSPD(KLCL))
- IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
- SHSIGN=1.
- ELSE
- SHSIGN=-1.
- ENDIF
- VWS=(U0(LTOP)-U0(KLCL))*(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL))* &
- (V0(LTOP)-V0(KLCL))
- VWS=1.E3*SHSIGN*SQRT(VWS)/(Z0(LTOP)-Z0(LCL))
- PEF=1.591+VWS*(-.639+VWS*(9.53E-2-VWS*4.96E-3))
- PEF=AMAX1(PEF,.2)
- PEF=AMIN1(PEF,.9)
- !
- !...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
- !
- CBH=(ZLCL-Z0(1))*3.281E-3
- IF(CBH.LT.3.)THEN
- RCBH=.02
- ELSE
- RCBH=.96729352+CBH*(-.70034167+CBH*(.162179896+CBH*(- &
- 1.2569798E-2+CBH*(4.2772E-4-CBH*5.44E-6))))
- ENDIF
- IF(CBH.GT.25)RCBH=2.4
- PEFCBH=1./(1.+RCBH)
- PEFCBH=AMIN1(PEFCBH,.9)
- !
- !... MEAN PEF. IS USED TO COMPUTE RAINFALL.
- !
- PEFF=.5*(PEF+PEFCBH)
- PEFF2=PEFF
- ! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
- !
- !*****************************************************************
- ! *
- ! COMPUTE DOWNDRAFT PROPERTIES *
- ! *
- !*****************************************************************
- !
- !...LET DOWNDRAFT ORIGINATE AT THE LEVEL OF MINIMUM SATURATION EQUIVALEN
- !...POTENTIAL TEMPERATURE (SEQT) IN THE CLOUD LAYER, EXTEND DOWNWARD TO
- !...SURFACE, OR TO THE LAYER BELOW CLOUD BASE AT WHICH ENVIR SEQT IS LES
- !...THAN MIN SEQT IN THE CLOUD LAYER...LET DOWNDRAFT DETRAIN OVER A LAYE
- !...OF SPECIFIED PRESSURE-DEPTH (DPDD)...
- !
- TDER=0.
- KSTART=MAX0(KPBL,KLCL)
- THTMIN=THTES(KSTART+1)
- KMIN=KSTART+1
- DO 104 NK=KSTART+2,LTOP-1
- THTMIN=AMIN1(THTMIN,THTES(NK))
- IF(THTMIN.EQ.THTES(NK))KMIN=NK
- 104 CONTINUE
- LFS=KMIN
- IF(RATIO2(LFS).GT.0.)CALL ENVIRTHT(P0(LFS),T0(LFS),Q0(LFS), &
- THETEE(LFS),0.,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
- EQFRC(LFS)=(THTES(LFS)-THETEU(LFS))/(THETEE(LFS)-THETEU(LFS))
- EQFRC(LFS)=AMAX1(EQFRC(LFS),0.)
- EQFRC(LFS)=AMIN1(EQFRC(LFS),1.)
- THETED(LFS)=THTES(LFS)
- !
- !...ESTIMATE THE EFFECT OF MELTING PRECIPITATION IN THE DOWNDRAFT...
- !
- IF(ML.GT.0)THEN
- DTMLTD=0.5*(QU(KLCL)-QU(LTOP))*RLF/CP
- ELSE
- DTMLTD=0.
- ENDIF
- TZ(LFS)=T0(LFS)-DTMLTD
- ES=ALIQ*EXP((TZ(LFS)*BLIQ-CLIQ)/(TZ(LFS)-DLIQ))
- QS=EP2*ES/(P0(LFS)-ES)
- QD(LFS)=EQFRC(LFS)*Q0(LFS)+(1.-EQFRC(LFS))*QU(LFS)
- THTAD(LFS)=TZ(LFS)*(P00/P0(LFS))**(0.2854*(1.-0.28*QD(LFS)))
- IF(QD(LFS).GE.QS)THEN
- THETED(LFS)=THTAD(LFS)* &
- EXP((3374.6525/TZ(LFS)-2.5403)*QS*(1.+0.81*QS))
- ELSE
- CALL ENVIRTHT(P0(LFS),TZ(LFS),QD(LFS),THETED(LFS),0.,RL,EP2,ALIQ, &
- BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
- ENDIF
- DO 107 NK=1,LFS
- ND=LFS-NK
- IF(THETED(LFS).GT.THTES(ND).OR.ND.EQ.1)THEN
- LDB=ND
- !
- !...IF DOWNDRAFT NEVER BECOMES NEGATIVELY BUOYANT OR IF IT
- !...IS SHALLOWER 50 mb, DON'T ALLOW IT TO OCCUR AT ALL...
- !
- IF(NK.EQ.1.OR.(P0(LDB)-P0(LFS)).LT.50.E2)GOTO 141
- ! testing ---- no downdraft
- ! GOTO 141
- GOTO 110
- ENDIF
- 107 CONTINUE
- !
- !...ALLOW DOWNDRAFT TO DETRAIN IN A SINGLE LAYER, BUT WITH DOWNDRAFT AIR
- !...TYPICALLY FLUSHED UP INTO HIGHER LAYERS AS ALLOWED IN THE TOTAL
- !...VERTICAL ADVECTION CALCULATIONS FARTHER DOWN IN THE CODE...
- !
- 110 DPDD=DP(LDB)
- LDT=LDB
- FRC=1.
- DPT=0.
- ! DO 115 NK=LDB,LFS
- ! DPT=DPT+DP(NK)
- ! IF(DPT.GT.DPDD)THEN
- ! LDT=NK
- ! FRC=(DPDD+DP(NK)-DPT)/DP(NK)
- ! GOTO 120
- ! ENDIF
- ! IF(NK.EQ.LFS-1)THEN
- ! LDT=NK
- ! FRC=1.
- ! DPDD=DPT
- ! GOTO 120
- ! ENDIF
- !115 CONTINUE
- 120 CONTINUE
- !
- !...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX..
- !
- TVD(LFS)=T0(LFS)*(1.+0.608*QES(LFS))
- RDD=P0(LFS)/(R*TVD(LFS))
- A1=(1.-PEFF)*AU0
- DMF(LFS)=-A1*RDD
- DER(LFS)=EQFRC(LFS)*DMF(LFS)
- DDR(LFS)=0.
- DO 140 ND=LFS-1,LDB,-1
- ND1=ND+1
- IF(ND.LE.LDT)THEN
- DER(ND)=0.
- DDR(ND)=-DMF(LDT+1)*DP(ND)*FRC/DPDD
- DMF(ND)=DMF(ND1)+DDR(ND)
- FRC=1.
- THETED(ND)=THETED(ND1)
- QD(ND)=QD(ND1)
- ELSE
- DER(ND)=DMF(LFS)*0.03*DP(ND)/RAD
- DDR(ND)=0.
- DMF(ND)=DMF(ND1)+DER(ND)
- IF(RATIO2(ND).GT.0.)CALL ENVIRTHT(P0(ND),T0(ND),Q0(ND), &
- THETEE(ND),0.,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
- THETED(ND)=(THETED(ND1)*DMF(ND1)+THETEE(ND)*DER(ND))/DMF(ND)
- QD(ND)=(QD(ND1)*DMF(ND1)+Q0(ND)*DER(ND))/DMF(ND)
- ENDIF
- 140 CONTINUE
- TDER=0.
- !
- !...CALCULATION AN EVAPORATION RATE FOR GIVEN MASS FLUX...
- !
- DO 135 ND=LDB,LDT
- TZ(ND)= &
- TPDD(P0(ND),THETED(LDT),T0(ND),QS,QD(ND),1.0,XLV0,XLV1, &
- EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
- ES=ALIQ*EXP((TZ(ND)*BLIQ-CLIQ)/(TZ(ND)-DLIQ))
- QS=EP2*ES/(P0(ND)-ES)
- DSSDT=(CLIQ-BLIQ*DLIQ)/((TZ(ND)-DLIQ)*(TZ(ND)-DLIQ))
- RL=XLV0-XLV1*TZ(ND)
- DTMP=RL*QS*(1.-RHBC)/(CP+RL*RHBC*QS*DSSDT)
- T1RH=TZ(ND)+DTMP
- ES=RHBC*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
- QSRH=EP2*ES/(P0(ND)-ES)
- !
- !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
- !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
- !
- IF(QSRH.LT.QD(ND))THEN
- QSRH=QD(ND)
- ! T1RH=T1+(QS-QSRH)*RL/CP
- T1RH=TZ(ND)
- ENDIF
- TZ(ND)=T1RH
- QS=QSRH
- TDER=TDER+(QS-QD(ND))*DDR(ND)
- QD(ND)=QS
- 135 THTAD(ND)=TZ(ND)*(P00/P0(ND))**(0.2854*(1.-0.28*QD(ND)))
- !
- !...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
- !...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
- !
- 141 IF(TDER.LT.1.)THEN
- ! WRITE(98,3004)I,J
- 3004 FORMAT(' ','I=',I3,2X,'J=',I3)
- PPTFLX=TRPPT
- CPR=TRPPT
- TDER=0.
- CNDTNF=0.
- UPDINC=1.
- LDB=LFS
- DO 117 NDK=1,LTOP
- DMF(NDK)=0.
- DER(NDK)=0.
- DDR(NDK)=0.
- THTAD(NDK)=0.
- WD(NDK)=0.
- TZ(NDK)=0.
- 117 QD(NDK)=0.
- AINCM2=100.
- GOTO 165
- ENDIF
- !
- !...ADJUST DOWNDRAFT MASS FLUX SO THAT EVAPORATION RATE IN DOWNDRAFT IS
- !...CONSISTENT WITH PRECIPITATION EFFICIENCY RELATIONSHIP...
- !
- DEVDMF=TDER/DMF(LFS)
- PPR=0.
- PPTFLX=PEFF*USR
- RCED=TRPPT-PPTFLX
- !
- !...PPR IS THE TOTAL AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE
- !...UPDRAFT FROM CLOUD BASE TO THE LFS...UPDRAFT MASS FLUX WILL BE
- !...INCREASED UP TO THE LFS TO ACCOUNT FOR UPDRAFT AIR MIXING WITH
- !...ENVIRONMENTAL AIR TO THE UPDRAFT, SO PPR WILL INCREASE
- !...PROPORTIONATELY...
- !
- DO 132 NM=KLCL,LFS
- 132 PPR=PPR+PPTLIQ(NM)+PPTICE(NM)
- IF(LFS.GE.KLCL)THEN
- DPPTDF=(1.-PEFF)*PPR*(1.-EQFRC(LFS))/UMF(LFS)
- ELSE
- DPPTDF=0.
- ENDIF
- !
- !...CNDTNF IS THE AMOUNT OF CONDENSATE TRANSFERRED ALONG WITH UPDRAFT
- !...MASS THE DOWNDRAFT AT THE LFS...
- !
- CNDTNF=(QLIQ(LFS)+QICE(LFS))*(1.-EQFRC(LFS))
- DMFLFS=RCED/(DEVDMF+DPPTDF+CNDTNF)
- IF(DMFLFS.GT.0.)THEN
- TDER=0.
- GOTO 141
- ENDIF
- !
- !...DDINC IS THE FACTOR BY WHICH TO INCREASE THE FIRST-GUESS DOWNDRAFT
- !...MASS FLUX TO SATISFY THE PRECIP EFFICIENCY RELATIONSHIP, UPDINC IS T
- !...WHICH TO INCREASE THE UPDRAFT MASS FLUX BELOW THE LFS TO ACCOUNT FOR
- !...TRANSFER OF MASS FROM UPDRAFT TO DOWNDRAFT...
- !
- ! DDINC=DMFLFS/DMF(LFS)
- IF(LFS.GE.KLCL)THEN
- UPDINC=(UMF(LFS)-(1.-EQFRC(LFS))*DMFLFS)/UMF(LFS)
- !
- !...LIMIT UPDINC TO LESS THAN OR EQUAL TO 1.5...
- !
- IF(UPDINC.GT.1.5)THEN
- UPDINC=1.5
- DMFLFS2=UMF(LFS)*(UPDINC-1.)/(EQFRC(LFS)-1.)
- RCED2=DMFLFS2*(DEVDMF+DPPTDF+CNDTNF)
- PPTFLX=PPTFLX+(RCED-RCED2)
- PEFF2=PPTFLX/USR
- RCED=RCED2
- DMFLFS=DMFLFS2
- ENDIF
- ELSE
- UPDINC=1.
- ENDIF
- DDINC=DMFLFS/DMF(LFS)
- DO 149 NK=LDB,LFS
- DMF(NK)=DMF(NK)*DDINC
- DER(NK)=DER(NK)*DDINC
- DDR(NK)=DDR(NK)*DDINC
- 149 CONTINUE
- CPR=TRPPT+PPR*(UPDINC-1.)
- PPTFLX=PPTFLX+PEFF*PPR*(UPDINC-1.)
- PEFF=PEFF2
- TDER=TDER*DDINC
- !
- !...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
- ! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
- ! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
- !
- DO 155 NK=LC,LFS
- UMF(NK)=UMF(NK)*UPDINC
- UDR(NK)=UDR(NK)*UPDINC
- UER(NK)=UER(NK)*UPDINC
- PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
- PPTICE(NK)=PPTICE(NK)*UPDINC
- DETLQ(NK)=DETLQ(NK)*UPDINC
- 155 DETIC(NK)=DETIC(NK)*UPDINC
- !
- !...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
- !...DOWNDRAFT...
- !
- IF(LDB.GT.1)THEN
- DO 156 NK=1,LDB-1
- DMF(NK)=0.
- DER(NK)=0.
- DDR(NK)=0.
- WD(NK)=0.
- TZ(NK)=0.
- QD(NK)=0.
- THTAD(NK)=0.
- 156 CONTINUE
- ENDIF
- DO 157 NK=LFS+1,KX
- DMF(NK)=0.
- DER(NK)=0.
- DDR(NK)=0.
- WD(NK)=0.
- TZ(NK)=0.
- QD(NK)=0.
- THTAD(NK)=0.
- 157 CONTINUE
- DO 158 NK=LDT+1,LFS-1
- TZ(NK)=0.
- QD(NK)=0.
- 158 CONTINUE
- !
- !
- !...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE
- ! INFLOW INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN
- ! IS AVAILABLE IN THAT LAYER INITIALLY...
- !
- 165 AINCMX=1000.
- LMAX=MAX0(KLCL,LFS)
- DO 166 NK=LC,LMAX
- IF((UER(NK)-DER(NK)).GT.0.)AINCM1=EMS(NK)/((UER(NK)-DER(NK))* &
- TIMEC)
- AINCMX=AMIN1(AINCMX,AINCM1)
- 166 CONTINUE
- AINC=1.
- IF(AINCMX.LT.AINC)AINC=AINCMX
- !
- !...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRFT AND DOWNDRFT...THEY
- !...WILL ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE
- !...STABILIZATION CLOSURE...
- !
- NCOUNT=0
- TDER2=TDER
- PPTFL2=PPTFLX
- DO 170 NK=1,LTOP
- DETLQ2(NK)=DETLQ(NK)
- DETIC2(NK)=DETIC(NK)
- UDR2(NK)=UDR(NK)
- UER2(NK)=UER(NK)
- DDR2(NK)=DDR(NK)
- DER2(NK)=DER(NK)
- UMF2(NK)=UMF(NK)
- DMF2(NK)=DMF(NK)
- 170 CONTINUE
- FABE=1.
- STAB=0.95
- NOITR=0
- IF(AINC/AINCMX.GT.0.999)THEN
- NCOUNT=0
- GOTO 255
- ENDIF
- ISTOP=0
- 175 NCOUNT=NCOUNT+1
- !
- !*****************************************************************
- ! *
- ! COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE *
- ! *
- !*****************************************************************
- !
- !...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
- !...SATISFY MASS CONTINUITY...
- !
- 185 CONTINUE
- DTT=TIMEC
- DO 200 NK=1,LTOP
- DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
- IF(NK.GT.1)THEN
- OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
- DTT1=0.75*DP(NK-1)/(ABS(OMG(NK))+1.E-10)
- DTT=AMIN1(DTT,DTT1)
- ENDIF
- 200 CONTINUE
- DO 488 NK=1,LTOP
- THPA(NK)=THTA0(NK)
- QPA(NK)=Q0(NK)
- NSTEP=NINT(TIMEC/DTT+1)
- DTIME=TIMEC/FLOAT(NSTEP)
- FXM(NK)=OMG(NK)*DXSQ/G
- 488 CONTINUE
- !
- !...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
- !
- DO 495 NTC=1,NSTEP
- !
- !...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED
- !...SIGN OF OMEGA...
- !
- DO 493 NK=1,LTOP
- THFXTOP(NK)=0.
- THFXBOT(NK)=0.
- QFXTOP(NK)=0.
- 493 QFXBOT(NK)=0.
- DO 494 NK=2,LTOP
- IF(OMG(NK).LE.0.)THEN
- THFXBOT(NK)=-FXM(NK)*THPA(NK-1)
- QFXBOT(NK)=-FXM(NK)*QPA(NK-1)
- THFXTOP(NK-1)=THFXTOP(NK-1)-THFXBOT(NK)
- QFXTOP(NK-1)=QFXTOP(NK-1)-QFXBOT(NK)
- ELSE
- THFXBOT(NK)=-FXM(NK)*THPA(NK)
- QFXBOT(NK)=-FXM(NK)*QPA(NK)
- THFXTOP(NK-1)=THFXTOP(NK-1)-THFXBOT(NK)
- QFXTOP(NK-1)=QFXTOP(NK-1)-QFXBOT(NK)
- ENDIF
- 494 CONTINUE
- !
- !...UPDATE THE THETA AND QV VALUES AT EACH LEVEL..
- !
- DO 492 NK=1,LTOP
- THPA(NK)=THPA(NK)+(THFXBOT(NK)+UDR(NK)*THTAU(NK)+DDR(NK)* &
- THTAD(NK)+THFXTOP(NK)-(UER(NK)-DER(NK))*THTA0(NK))* &
- DTIME*EMSD(NK)
- QPA(NK)=QPA(NK)+(QFXBOT(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK)+ &
- QFXTOP(NK)-(UER(NK)-DER(NK))*Q0(NK))*DTIME*EMSD(NK)
- 492 CONTINUE
- 495 CONTINUE
- DO 498 NK=1,LTOP
- THTAG(NK)=THPA(NK)
- QG(NK)=QPA(NK)
- 498 CONTINUE
- !
- !...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO,
- !...BORROW MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO.
- !
- DO 499 NK=1,LTOP
- IF(QG(NK).LT.0.)THEN
- IF(NK.EQ.1)THEN
- CALL wrf_error_fatal ( 'module_cu_kf.F: problem with kf scheme: qg = 0 at the surface' )
- ENDIF
- NK1=NK+1
- IF(NK.EQ.LTOP)NK1=KLCL
- TMA=QG(NK1)*EMS(NK1)
- TMB=QG(NK-1)*EMS(NK-1)
- TMM=(QG(NK)-1.E-9)*EMS(NK)
- BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
- ACOEFF=BCOEFF*TMA/TMB
- TMB=TMB*(1.-BCOEFF)
- TMA=TMA*(1.-ACOEFF)
- IF(NK.EQ.LTOP)THEN
- QVDIFF=(QG(NK1)-TMA*EMSD(NK1))*100./QG(NK1)
- IF(ABS(QVDIFF).GT.1.)THEN
- PRINT *,'--WARNING-- CLOUD BASE WATER VAPOR CHANGES BY ', &
- QVDIFF, &
- ' PERCENT WHEN MOISTURE IS BORROWED TO PREVENT NEG VALUES', &
- ' IN KAIN-FRITSCH'
- ENDIF
- ENDIF
- QG(NK)=1.E-9
- QG(NK1)=TMA*EMSD(NK1)
- QG(NK-1)=TMB*EMSD(NK-1)
- ENDIF
- 499 CONTINUE
- TOPOMG=(UDR(LTOP)-UER(LTOP))*DP(LTOP)*EMSD(LTOP)
- IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
- ! WRITE(98,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME;'
- ! * ,'TOPOMG, OMG =',TOPOMG,OMG(LTOP)
- WRITE(6,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME;' &
- ,'TOPOMG, OMG =',TOPOMG,OMG(LTOP)
- ISTOP=1
- GOTO 265
- ENDIF
- !
- !...CONVERT THETA TO T...
- !
- ! PAY ATTENTION ...
- !
- DO 230 NK=1,LTOP
- EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QG(NK)))
- TG(NK)=THTAG(NK)/EXN(NK)
- TVG(NK)=TG(NK)*(1.+0.608*QG(NK))
- 230 CONTINUE
- !
- !*******************************************************************
- ! *
- ! COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY. *
- ! *
- !*******************************************************************
- !
- !...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
- !
- THMIX=0.
- QMIX=0.
- PMIX=0.
- DO 217 NK=LC,KPBL
- ROCPQ=0.2854*(1.-0.28*QG(NK))
- THMIX=THMIX+DP(NK)*TG(NK)*(P00/P0(NK))**ROCPQ
- QMIX=QMIX+DP(NK)*QG(NK)
- 217 PMIX=PMIX+DP(NK)*P0(NK)
- THMIX=THMIX/DPTHMX
- QMIX=QMIX/DPTHMX
- PMIX=PMIX/DPTHMX
- ROCPQ=0.2854*(1.-0.28*QMIX)
- TMIX=THMIX*(PMIX/P00)**ROCPQ
- ES=ALIQ*EXP((TMIX*BLIQ-CLIQ)/(TMIX-DLIQ))
- QS=EP2*ES/(PMIX-ES)
- !
- !...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
- !
- IF(QMIX.GT.QS)THEN
- RL=XLV0-XLV1*TMIX
- CPM=CP*(1.+0.887*QMIX)
- DSSDT=QS*(CLIQ-BLIQ*DLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
- DQ=(QMIX-QS)/(1.+RL*DSSDT/CPM)
- TMIX=TMIX+RL/CP*DQ
- QMIX=QMIX-DQ
- ROCPQ=0.2854*(1.-0.28*QMIX)
- THMIX=TMIX*(P00/PMIX)**ROCPQ
- TLCL=TMIX
- PLCL=PMIX
- ELSE
- QMIX=AMAX1(QMIX,0.)
- EMIX=QMIX*PMIX/(EP2+QMIX)
- TLOG=ALOG(EMIX/ALIQ)
- TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
- TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX- &
- TDPT)
- TLCL=AMIN1(TLCL,TMIX)
- CPORQ=1./ROCPQ
- PLCL=P00*(TLCL/THMIX)**CPORQ
- ENDIF
- TVLCL=TLCL*(1.+0.608*QMIX)
- DO 235 NK=LC,KL
- KLCL=NK
- 235 IF(PLCL.GE.P0(NK))GOTO 240
- 240 K=KLCL-1
- DLP=ALOG(PLCL/P0(K))/ALOG(P0(KLCL)/P0(K))
- !
- !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
- !
- TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
- QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
- TVEN=TENV*(1.+0.608*QENV)
- TVBAR=0.5*(TVG(K)+TVEN)
- ! ZLCL=Z0(K)+R*TVBAR*ALOG(P0(K)/PLCL)/G
- ZLCL=Z0(K)+(Z0(KLCL)-Z0(K))*DLP
- TVAVG=0.5*(TVEN+TG(KLCL)*(1.+0.608*QG(KLCL)))
- PLCL=P0(KLCL)*EXP(G/(R*TVAVG)*(Z0(KLCL)-ZLCL))
- THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
- EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
- ES=ALIQ*EXP((TENV*BLIQ-CLIQ)/(TENV-DLIQ))
- QESE=EP2*ES/(PLCL-ES)
- THTESG(K)=TENV*(1.E5/PLCL)**(0.2854*(1.-0.28*QESE))* &
- EXP((3374.6525/TENV-2.5403)*QESE*(1.+0.81*QESE))
- !
- !...COMPUTE ADJUSTED ABE(ABEG).
- !
- ABEG=0.
- THTUDL=THETEU(K)
- DO 245 NK=K,LTOPM1
- NK1=NK+1
- ES=ALIQ*EXP((TG(NK1)*BLIQ-CLIQ)/(TG(NK1)-DLIQ))
- QESE=EP2*ES/(P0(NK1)-ES)
- THTESG(NK1)=TG(NK1)*(1.E5/P0(NK1))**(0.2854*(1.-0.28*QESE))* &
- EXP((3374.6525/TG(NK1)-2.5403)*QESE*(1.+0.81*QESE) &
- )
- ! DZZ=CVMGT(Z0(KLCL)-ZLCL,DZA(NK),NK.EQ.K)
- IF(NK.EQ.K)THEN
- DZZ=Z0(KLCL)-ZLCL
- ELSE
- DZZ=DZA(NK)
- ENDIF
- BE=((2.*THTUDL)/(THTESG(NK1)+THTESG(NK))-1.)*DZZ
- 245 IF(BE.GT.0.)ABEG=ABEG+BE*G
- !
- !...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
- !...THE PERIOD TIMEC...
- !
- IF(NOITR.EQ.1)THEN
- ! WRITE(98,1060)FABE
- GOTO 265
- ENDIF
- DABE=AMAX1(ABE-ABEG,0.1*ABE)
- FABE=ABEG/(ABE+1.E-8)
- IF(FABE.GT.1.)THEN
- ! WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS '
- ! *,'GRID POINT; NO CONVECTION ALLOWED!'
- GOTO 325
- ENDIF
- IF(NCOUNT.NE.1)THEN
- DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
- IF(DFDA.GT.0.)THEN
- NOITR=1
- AINC=AINCOLD
- GOTO 255
- ENDIF
- ENDIF
- AINCOLD=AINC
- FABEOLD=FABE
- IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
- ! WRITE(98,1055)FABE
- GOTO 265
- ENDIF
- IF(FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB)GOTO 265
- IF(NCOUNT.GT.10)THEN
- ! WRITE(98,1060)FABE
- GOTO 265
- ENDIF
- !
- !...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE
- !...CONVECTIVE MASS FLUX BY THE FACTOR AINC:
- !
- IF(FABE.EQ.0.)THEN
- AINC=AINC*0.5
- ELSE
- AINC=AINC*STAB*ABE/(DABE+1.E-8)
- ENDIF
- 255 AINC=AMIN1(AINCMX,AINC)
- !...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION
- !...WILL BE MINIMAL SO JUST IGNORE IT...
- IF(AINC.LT.0.05)GOTO 325
- ! AINC=AMAX1(AINC,0.05)
- TDER=TDER2*AINC
- PPTFLX=PPTFL2*AINC
- ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABEOLD,AINCOLD
- DO 260 NK=1,LTOP
- UMF(NK)=UMF2(NK)*AINC
- DMF(NK)=DMF2(NK)*AINC
- DETLQ(NK)=DETLQ2(NK)*AINC
- DETIC(NK)=DETIC2(NK)*AINC
- UDR(NK)=UDR2(NK)*AINC
- UER(NK)=UER2(NK)*AINC
- DER(NK)=DER2(NK)*AINC
- DDR(NK)=DDR2(NK)*AINC
- 260 CONTINUE
- !
- !...GO BACK UP FOR ANOTHER ITERATION...
- !
- GOTO 175
- 265 CONTINUE
- !
- !...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
- !...GRID POINT...
- !
- !...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
- !
- !...FRC2 IS THE FRACTION OF TOTAL CONDENSATE
- !...GENERATED THAT GOES INTO PRECIPITIATION
- FRC2=PPTFLX/(CPR*AINC)
- DO 270 NK=1,LTOP
- QLPA(NK)=QL0(NK)
- QIPA(NK)=QI0(NK)
- QRPA(NK)=QR0(NK)
- QSPA(NK)=QS0(NK)
- RAINFB(NK)=PPTLIQ(NK)*AINC*FBFRC*FRC2
- SNOWFB(NK)=PPTICE(NK)*AINC*FBFRC*FRC2
- 270 CONTINUE
- DO 290 NTC=1,NSTEP
- !
- !...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH
- !...LAYER BASED ON THE SIGN OF OMEGA...
- !
- DO 275 NK=1,LTOP
- QLFXIN(NK)=0.
- QLFXOUT(NK)=0.
- QIFXIN(NK)=0.
- QIFXOUT(NK)=0.
- QRFXIN(NK)=0.
- QRFXOUT(NK)=0.
- QSFXIN(NK)=0.
- QSFXOUT(NK)=0.
- 275 CONTINUE
- DO 280 NK=2,LTOP
- IF(OMG(NK).LE.0.)THEN
- QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
- QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
- QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
- QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
- QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
- QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
- QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
- QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
- ELSE
- QLFXOUT(NK)=FXM(NK)*QLPA(NK)
- QIFXOUT(NK)=FXM(NK)*QIPA(NK)
- QRFXOUT(NK)=FXM(NK)*QRPA(NK)
- QSFXOUT(NK)=FXM(NK)*QSPA(NK)
- QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
- QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
- QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
- QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
- ENDIF
- 280 CONTINUE
- !
- !...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
- !
- DO 285 NK=1,LTOP
- QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))*DTIME* &
- EMSD(NK)
- QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))*DTIME* &
- EMSD(NK)
- QRPA(NK)=QRPA(NK)+(QRFXIN(NK)+QLQOUT(NK)*UDR(NK)-QRFXOUT(NK) &
- +RAINFB(NK))*DTIME*EMSD(NK)
- QSPA(NK)=QSPA(NK)+(QSFXIN(NK)+QICOUT(NK)*UDR(NK)-QSFXOUT(NK) &
- +SNOWFB(NK))*DTIME*EMSD(NK)
- 285 CONTINUE
- 290 CONTINUE
- DO 295 NK=1,LTOP
- QLG(NK)=QLPA(NK)
- QIG(NK)=QIPA(NK)
- QRG(NK)=QRPA(NK)
- QSG(NK)=QSPA(NK)
- 295 CONTINUE
- ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABE,AINC
- !
- !...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
- !
- IF(ISTOP.EQ.1)THEN
- WRITE(6,1070)' P ',' DP ',' DT K/D ',' DR K/D ',' OMG ', &
- ' DOMGDP ',' UMF ',' UER ',' UDR ',' DMF ',' DER ' &
- ,' DDR ',' EMS ',' W0 ',' DETLQ ',' DETIC '
- DO 300 K=LTOP,1,-1
- DTT=(TG(K)-T0(K))*86400./TIMEC
- RL=XLV0-XLV1*TG(K)
- DR=-(QG(K)-Q0(K))*RL*86400./(TIMEC*CP)
- UDFRC=UDR(K)*TIMEC*EMSD(K)
- UEFRC=UER(K)*TIMEC*EMSD(K)
- DDFRC=DDR(K)*TIMEC*EMSD(K)
- DEFRC=-DER(K)*TIMEC*EMSD(K)
- WRITE (6,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)* &
- 1.E4,UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC &
- ,DDFRC,EMS(K)/1.E11,W0AVG1D(K)*1.E2,DETLQ(K) &
- *TIMEC*EMSD(K)*1.E3,DETIC(K)*TIMEC*EMSD(K)* &
- 1.E3
- 300 CONTINUE
- WRITE(6,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG', &
- 'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
- DO 305 K=KX,1,-1
- DTT=TG(K)-T0(K)
- TUC=TU(K)-T00
- IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
- TDC=TZ(K)-T00
- IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
- ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
- QGS=ES*EP2/(P0(K)-ES)
- RH0=Q0(K)/QES(K)
- RHG=QG(K)/QGS
- WRITE (6,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC &
- ,TDC,Q0(K)*1000.,QG(K)*1000.,(QG(K)-Q0(K))* &
- 1000.,QU(K)*1000.,QD(K)*1000.,QLG(K)*1000., &
- QIG(K)*1000.,QRG(K)*1000.,QSG(K)*1000.,RH0,RHG
- 305 CONTINUE
- !
- !...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
- !...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
- !
- IF(ISTOP.EQ.1)THEN
- DO 310 K=1,KX
- WRITE ( wrf_err_message , 1115 ) &
- Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
- U0(K),V0(K),DP(K)/100.,W0AVG1D(K)
- CALL wrf_message ( TRIM( wrf_err_message ) )
- 310 CONTINUE
- CALL wrf_error_fatal ( 'module_cu_kf.F: KAIN-FRITSCH' )
- ENDIF
- ENDIF
- CNDTNF=(1.-EQFRC(LFS))*(QLIQ(LFS)+QICE(LFS))*DMF(LFS)
- ! WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
- !
- ! EVALUATE MOISTURE BUDGET...
- !
- QINIT=0.
- QFNL=0.
- DPT=0.
- DO 315 NK=1,LTOP
- DPT=DPT+DP(NK)
- QINIT=QINIT+Q0(NK)*EMS(NK)
- QFNL=QFNL+QG(NK)*EMS(NK)
- QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
- 315 CONTINUE
- QFNL=QFNL+PPTFLX*TIMEC*(1.-FBFRC)
- ERR2=(QFNL-QINIT)*100./QINIT
- ! WRITE(98,1110)QINIT,QFNL,ERR2
- ! IF(ABS(ERR2).GT.0.05)STOP 'QVERR'
- IF(ABS(ERR2).GT.0.05)CALL wrf_error_fatal( 'module_cu_kf.F: QVERR' )
- RELERR=ERR2*QINIT/(PPTFLX*TIMEC+1.E-10)
- ! WRITE(98,1120)RELERR
- ! WRITE(98,*)'TDER, CPR, USR, TRPPT =',
- ! *TDER,CPR*AINC,USR*AINC,TRPPT*AINC
- !
- !...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
- !
- !...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
- !...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
- !
- IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
- NCA(I,J)=FLOAT(NIC)*DT
- DO 320 K=1,KX
- ! IF(IMOIST.NE.2)THEN
- !
- !...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMAT
- !...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
- !...NOTE: THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
- !...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
- !...OF QG...
- !
- ! RLC=XLV0-XLV1*TG(K)
- ! RLS=XLS0-XLS1*TG(K)
- ! CPM=CP*(1.+0.887*QG(K))
- ! TG(K)=TG(K)-(RLC*(QLG(K)+QRG(K))+RLS*(QIG(K)+QSG(K)))/CPM
- ! QG(K)=QG(K)+(QLG(K)+QRG(K)+QIG(K)+QSG(K))
- ! DQCDT(K)=0.
- ! DQIDT(K)=0.
- ! DQRDT(K)=0.
- ! DQSDT(K)=0.
- ! ELSE
- IF(.NOT. qi_flag .and. warm_rain)THEN
- !
- !...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
- !
- CPM=CP*(1.+0.887*QG(K))
- TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
- DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
- DQIDT(K)=0.
- DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
- DQSDT(K)=0.
- ELSEIF(.NOT. qi_flag .and. .not. warm_rain)THEN
- !
- !...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROME
- !...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
- !
- CPM=CP*(1.+0.887*QG(K))
- IF(K.LE.ML)THEN
- TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
- ELSEIF(K.GT.ML)THEN
- TG(K)=TG(K)+(QLG(K)+QRG(K))*RLF/CPM
- ENDIF
- DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
- DQIDT(K)=0.
- DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
- DQSDT(K)=0.
- ELSEIF(qi_flag) THEN
- !
- !...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE
- !...TENDENCY OF HYDROMETEORS DIRECTLY...
- !
- DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
- DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
- DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
- IF (qs_flag ) THEN
- DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
- ELSE
- DQIDT(K)=DQIDT(K)+(QSG(K)-QS0(K))/TIMEC
- ENDIF
- ELSE
- CALL wrf_error_fatal ( 'module_cu_kf: THIS COMBINATION OF IMOIST, IICE NOT ALLOWED' )
- ENDIF
- ! ENDIF
- DTDT(K)=(TG(K)-T0(K))/TIMEC
- DQDT(K)=(QG(K)-Q0(K))/TIMEC
- 320 CONTINUE
- ! RAINCV is in the unit of mm
- PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
- RAINCV(I,J)=DT*PRATEC(I,J)
- RNC=RAINCV(I,J)*NIC
- ! WRITE(98,909)RNC
- 909 FORMAT(' CONVECTIVE RAINFALL =',F8.4,' CM')
- 325 CONTINUE
- 1000 FORMAT(' ',10A8)
- 1005 FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
- 1010 FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
- 1015 FORMAT(' ','ALL REMAINING MASS DETRAINS BELOW ',F4.0,' MB')
- 1025 FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M', &
- ' DTLCL=',F5.2,' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=', &
- I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1, &
- ' CAPE=',0PF7.1)
- 1030 FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =', &
- E12.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =', &
- F8.1)
- 1035 FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'WKL=' &
- ,F6.3,'VWS=',F5.2)
- 1040 FORMAT(' ','PRECIP EFF = 100%, ENVIR CANNOT SUPPORT DOWND' &
- ,'RAFTS')
- !1045 FORMAT('NUMBER OF DOWNDRAFT ITERATIONS EXCEEDS 10...PPTFLX' &
- ! ' IS DIFFERENT FROM THAT GIVEN BY PRECIP EFF RELATION')
- ! FLIC HAS TROUBLE WITH THIS ONE.
- 1045 FORMAT('NUMBER OF DOWNDRAFT ITERATIONS EXCEEDS 10')
- 1050 FORMAT(' ','LCOUNT= ',I3,' PPTFLX/CPR, PEFF= ',F5.3,1X,F5.3, &
- 'DMF(LFS)/UMF(LCL)= ',F5.3)
- 1055 FORMAT(/'*** DEGREE OF STABILIZATION =',F5.3,', NO MORE MASS F' &
- ,'LUX IS ALLOWED')
- !1060 FORMAT(/' ITERATION DOES NOT CONVERGE TO GIVE THE SPECIFIED ' &
- ! 'DEGREE OF STABILIZATION! FABE= ',F6.4)
- 1060 FORMAT(/' ITERATION DOES NOT CONVERGE. FABE= ',F6.4)
- 1070 FORMAT (16A8)
- 1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3)
- 1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, NSTEP=',F5.0,I3, &
- 'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
- 1085 FORMAT (A3,16A7,2A8)
- 1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3)
- 1095 FORMAT(' ',' PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ', &
- F10.0)
- 1105 FORMAT(' ','NET LATENT HEAT RELEASE =',E12.5,' ACTUAL HEATING =', &
- E12.5,' J/KG-S, DIFFERENCE = ',F9.3,'PERCENT')
- 1110 FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5, &
- ' TOTAL WATER CHANGE =',F8.2,'PERCENT')
- 1115 FORMAT (2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4 &
- )
- 1120 FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3, &
- 'PERCENT')
- END SUBROUTINE KFPARA
- !-----------------------------------------------------------------------
- SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
- QNEWIC,QLQOUT,QICOUT,G)
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- ! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
- ! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
- ! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
- ! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
- ! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
- REAL, INTENT(IN ) :: G
- REAL, INTENT(IN ) :: DZ,BOTERM,ENTERM,RATE
- REAL, INTENT(INOUT) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
- REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
- QTOT=QLIQ+QICE
- QNEW=QNEWLQ+QNEWIC
- !
- ! ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY C
- ! BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
- ! LEVELS...
- !
- QEST=0.5*(QTOT+QNEW)
- G1=WTW+BOTERM-ENTERM-2.*G*DZ*QEST/1.5
- IF(G1.LT.0.0)G1=0.
- WAVG=(SQRT(WTW)+SQRT(G1))/2.
- CONV=RATE*DZ/WAVG
- !
- ! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
- ! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
- ! IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
- ! SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...
- !
- RATIO3=QNEWLQ/(QNEW+1.E-10)
- ! OLDQ=QTOT
- QTOT=QTOT+0.6*QNEW
- OLDQ=QTOT
- RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-10)
- QTOT=QTOT*EXP(-CONV)
- !
- ! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
- ! PARCEL AT THIS LEVEL...
- !
- DQ=OLDQ-QTOT
- QLQOUT=RATIO4*DQ
- QICOUT=(1.-RATIO4)*DQ
- !
- ! ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
- ! LATE VERTICAL VELOCITY
- !
- PPTDRG=0.5*(OLDQ+QTOT-0.2*QNEW)
- WTW=WTW+BOTERM-ENTERM-2.*G*DZ*PPTDRG/1.5
- !
- ! DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
- ! DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...
- !
- QLIQ=RATIO4*QTOT+RATIO3*0.4*QNEW
- QICE=(1.-RATIO4)*QTOT+(1.-RATIO3)*0.4*QNEW
- QNEWLQ=0.
- QNEWIC=0.
- END SUBROUTINE CONDLOAD
- !-----------------------------------------------------------------------
- SUBROUTINE DTFRZNEW(TU,P,THTEU,QVAP,QLIQ,QICE,RATIO2,TTFRZ,TBFRZ, &
- QNWFRZ,RL,FRC1,EFFQ,IFLAG,XLV0,XLV1,XLS0,XLS1, &
- EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- REAL, INTENT(IN ) :: XLV0,XLV1
- REAL, INTENT(IN ) :: P,TTFRZ,TBFRZ,EFFQ,XLS0,XLS1,EP2,ALIQ, &
- BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE
- REAL, INTENT(INOUT) :: TU,THTEU,QVAP,QLIQ,QICE,RATIO2, &
- FRC1,RL,QNWFRZ
- INTEGER, INTENT(INOUT) :: IFLAG
- REAL :: CCP,RV,C5,QLQFRZ,QNEW,ESLIQ,ESICE,RLC,RLS,PI,ES,RLF,A, &
- B,C,DQVAP,DTFRZ,TU1,QVAP1
- !-----------------------------------------------------------------------
- !
- !...ALLOW GLACIATION OF THE UPDRAFT TO OCCUR AS AN APPROXIMATELY LINEAR
- ! FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE TTFRZ TO TBFRZ...
- !
- RV=461.5
- C5=1.0723E-3
- !
- !...ADJUST THE LIQUID WATER CONCENTRATIONS FROM FRESH CONDENSATE AND THA
- ! BROUGHT UP FROM LOWER LEVELS TO AN AMOUNT THAT WOULD BE PRESENT IF N
- ! LIQUID WATER HAD FROZEN THUS FAR...THIS IS NECESSARY BECAUSE THE
- ! EXPRESSION FOR TEMP CHANGE IS MULTIPLIED BY THE FRACTION EQUAL TO TH
- ! PARCEL TEMP DECREASE SINCE THE LAST MODEL LEVEL DIVIDED BY THE TOTAL
- ! GLACIATION INTERVAL, SO THAT EFFECTIVELY THIS APPROXIMATELY ALLOWS A
- ! AMOUNT OF LIQUID WATER TO FREEZE WHICH IS EQUAL TO THIS SAME FRACTIO
- ! OF THE LIQUID WATER THAT WAS PRESENT BEFORE THE GLACIATION PROCESS W
- ! INITIATED...ALSO, TO ALLOW THETAU TO CONVERT APPROXIMATELY LINEARLY
- ! ITS VALUE WITH RESPECT TO ICE, WE NEED TO ALLOW A PORTION OF THE FRE
- ! CONDENSATE TO CONTRIBUTE TO THE GLACIATION PROCESS; THE FRACTIONAL
- ! AMOUNT THAT APPLIES TO THIS PORTION IS 1/2 OF THE FRACTIONAL AMOUNT
- ! FROZEN OF THE "OLD" CONDENSATE BECAUSE THIS FRESH CONDENSATE IS ONLY
- ! PRODUCED GRADUALLY OVER THE LAYER...NOTE THAT IN TERMS OF THE DYNAMI
- ! OF THE PRECIPITATION PROCESS, IE. PRECIPITATION FALLOUT, THIS FRACTI
- ! AMNT OF FRESH CONDENSATE HAS ALREADY BEEN INCLUDED IN THE ICE CATEGO
- !
- QLQFRZ=QLIQ*EFFQ
- QNEW=QNWFRZ*EFFQ*0.5
- ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
- ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
- RLC=2.5E6-2369.276*(TU-273.16)
- RLS=2833922.-259.532*(TU-273.16)
- RLF=RLS-RLC
- CCP=1005.7*(1.+0.89*QVAP)
- !
- ! A = D(ES)/DT IS THAT CALCULATED FROM BUCK`S (1981) EMPIRICAL FORMULAS
- ! FOR SATURATION VAPOR PRESSURE...
- !
- A=(CICE-BICE*DICE)/((TU-DICE)*(TU-DICE))
- B=RLS*EP2/P
- C=A*B*ESICE/CCP
- DQVAP=B*(ESLIQ-ESICE)/(RLS+RLS*C)-RLF*(QLQFRZ+QNEW)/(RLS+RLS/C)
- DTFRZ=(RLF*(QLQFRZ+QNEW)+B*(ESLIQ-ESICE))/(CCP+A*B*ESICE)
- TU1=TU
- QVAP1=QVAP
- TU=TU+FRC1*DTFRZ
- QVAP=QVAP-FRC1*DQVAP
- ES=QVAP*P/(EP2+QVAP)
- ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
- ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
- RATIO2=(ESLIQ-ES)/(ESLIQ-ESICE)
- !
- ! TYPICALLY, RATIO2 IS VERY CLOSE TO (TTFRZ-TU)/(TTFRZ-TBFRZ), USUALLY
- ! WITHIN 1% (USING TU BEFORE GALCIATION EFFECTS ARE APPLIED); IF THE
- ! INITIAL UPDRAFT TEMP IS BELOW TBFRZ AND RATIO2 IS STILL LESS THAN 1,
- ! AN ADJUSTMENT TO FRC1 AND RATIO2 IS INTRODUCED SO THAT GLACIATION
- ! EFFECTS ARE NOT UNDERESTIMATED; CONVERSELY, IF RATIO2 IS GREATER THAN
- ! FRC1 IS ADJUSTED SO THAT GLACIATION EFFECTS ARE NOT OVERESTIMATED...
- !
- IF(IFLAG.GT.0.AND.RATIO2.LT.1)THEN
- FRC1=FRC1+(1.-RATIO2)
- TU=TU1+FRC1*DTFRZ
- QVAP=QVAP1-FRC1*DQVAP
- RATIO2=1.
- IFLAG=1
- GOTO 20
- ENDIF
- IF(RATIO2.GT.1.)THEN
- FRC1=FRC1-(RATIO2-1.)
- FRC1=AMAX1(0.0,FRC1)
- TU=TU1+FRC1*DTFRZ
- QVAP=QVAP1-FRC1*DQVAP
- RATIO2=1.
- IFLAG=1
- ENDIF
- !
- ! CALCULATE A HYBRID VALUE OF THETAU, ASSUMING THAT THE LATENT HEAT OF
- ! VAPORIZATION/SUBLIMATION CAN BE ESTIMATED USING THE SAME WEIGHTING
- ! FUNCTION AS THAT USED TO CALCULATE SATURATION VAPOR PRESSURE, CALCU-
- ! LATE NEW LIQUID WATER AND ICE CONCENTRATIONS...
- !
- 20 RLC=XLV0-XLV1*TU
- RLS=XLS0-XLS1*TU
- RL=RATIO2*RLS+(1.-RATIO2)*RLC
- PI=(1.E5/P)**(0.2854*(1.-0.28*QVAP))
- THTEU=TU*PI*EXP(RL*QVAP*C5/TU*(1.+0.81*QVAP))
- IF(IFLAG.EQ.1)THEN
- QICE=QICE+FRC1*DQVAP+QLIQ
- QLIQ=0.
- ELSE
- QICE=QICE+FRC1*(DQVAP+QLQFRZ)
- QLIQ=QLIQ-FRC1*QLQFRZ
- ENDIF
- QNWFRZ=0.
-
- END SUBROUTINE DTFRZNEW
- !-----------------------------------------------------------------------
- !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
- ! THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN
- ! DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN F
- ! HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMA
- ! TABLES ED. BY ABRAMOWITZ AND STEGUN, NAT L BUREAU OF STANDARDS APPLI
- ! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
- ! JACK KAIN
- ! 7/6/89
- !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
- !***********************************************************************
- !***** GAUSSIAN TYPE MIXING PROFILE....******************************
- SUBROUTINE PROF5(EQ,EE,UD)
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- REAL, INTENT(IN ) :: EQ
- REAL, INTENT(INOUT) :: EE,UD
- REAL :: SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
- DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676, &
- 0.9372980,0.33267,0.166666667,0.202765151/
- X=(EQ-0.5)/SIGMA
- Y=6.*EQ-3.
- EY=EXP(Y*Y/(-2))
- E45=EXP(-4.5)
- T2=1./(1.+P*ABS(Y))
- T1=0.500498
- C1=A1*T1+A2*T1*T1+A3*T1*T1*T1
- C2=A1*T2+A2*T2*T2+A3*T2*T2*T2
- IF(Y.GE.0.)THEN
- EE=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
- UD=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*(0.5+EQ*EQ/2.- &
- EQ)
- ELSE
- EE=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
- UD=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*(0.5+EQ* &
- EQ/2.-EQ)
- ENDIF
- EE=EE/FE
- UD=UD/FE
- END SUBROUTINE PROF5
- !-----------------------------------------------------------------------
- SUBROUTINE TPMIX(P,THTU,TU,QU,QLIQ,QICE,QNEWLQ,QNEWIC,RATIO2,RL, &
- XLV0,XLV1,XLS0,XLS1, &
- EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- REAL, INTENT(IN ) :: XLV0,XLV1
- REAL, INTENT(IN ) :: P,THTU,RATIO2,RL,XLS0, &
- XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,&
- CICE,DICE
- REAL, INTENT(INOUT) :: QU,QLIQ,QICE,TU,QNEWLQ,QNEWIC
- REAL :: ES,QS,PI,THTGS,F0,T1,T0,C5,RV,ESLIQ,ESICE,F1,DT,QNEW, &
- DQ, QTOT,DQICE,DQLIQ,RLL,CCP
- INTEGER :: ITCNT
- !-----------------------------------------------------------------------
- !
- !...THIS SUBROUTINE ITERATIVELY EXTRACTS WET-BULB TEMPERATURE FROM EQUIV
- ! POTENTIAL TEMPERATURE, THEN CHECKS TO SEE IF SUFFICIENT MOISTURE IS
- ! AVAILABLE TO ACHIEVE SATURATION...IF NOT, TEMPERATURE IS ADJUSTED
- ! ACCORDINGLY, IF SO, THE RESIDUAL LIQUID WATER/ICE CONCENTRATION IS
- ! DETERMINED...
- !
- C5=1.0723E-3
- RV=461.5
- !
- ! ITERATE TO FIND WET BULB TEMPERATURE AS A FUNCTION OF EQUIVALENT POT
- ! TEMP AND PRS, ASSUMING SATURATION VAPOR PRESSURE...RATIO2 IS THE DEG
- ! OF GLACIATION...
- !
- IF(RATIO2.LT.1.E-6)THEN
- ES=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
- QS=EP2*ES/(P-ES)
- PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
- THTGS=TU*PI*EXP((3374.6525/TU-2.5403)*QS*(1.+0.81*QS))
- ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
- ES=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
- QS=EP2*ES/(P-ES)
- PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
- THTGS=TU*PI*EXP((3114.834/TU-0.278296)*QS*(1.+0.81*QS))
- ELSE
- ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
- ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
- ES=(1.-RATIO2)*ESLIQ+RATIO2*ESICE
- QS=EP2*ES/(P-ES)
- PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
- THTGS=TU*PI*EXP(RL*QS*C5/TU*(1.+0.81*QS))
- ENDIF
- F0=THTGS-THTU
- T1=TU-0.5*F0
- T0=TU
- ITCNT=0
- 90 IF(RATIO2.LT.1.E-6)THEN
- ES=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
- QS=EP2*ES/(P-ES)
- PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
- THTGS=T1*PI*EXP((3374.6525/T1-2.5403)*QS*(1.+0.81*QS))
- ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
- ES=AICE*EXP((BICE*T1-CICE)/(T1-DICE))
- QS=EP2*ES/(P-ES)
- PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
- THTGS=T1*PI*EXP((3114.834/T1-0.278296)*QS*(1.+0.81*QS))
- ELSE
- ESLIQ=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
- ESICE=AICE*EXP((BICE*T1-CICE)/(T1-DICE))
- ES=(1.-RATIO2)*ESLIQ+RATIO2*ESICE
- QS=EP2*ES/(P-ES)
- PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
- THTGS=T1*PI*EXP(RL*QS*C5/T1*(1.+0.81*QS))
- ENDIF
- F1=THTGS-THTU
- IF(ABS(F1).LT.0.01)GOTO 50
- ITCNT=ITCNT+1
- IF(ITCNT.GT.10)GOTO 50
- DT=F1*(T1-T0)/(F1-F0)
- T0=T1
- F0=F1
- T1=T1-DT
- GOTO 90
- !
- ! IF THE PARCEL IS SUPERSATURATED, CALCULATE CONCENTRATION OF FRESH
- ! CONDENSATE...
- !
- 50 IF(QS.LE.QU)THEN
- QNEW=QU-QS
- QU=QS
- GOTO 96
- ENDIF
- !
- ! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
- ! ADJUSTED...IF LIQUID WATER OR ICE IS PRESENT, IT IS ALLOWED TO EVAPO
- ! SUBLIMATE.
- !
- QNEW=0.
- DQ=QS-QU
- QTOT=QLIQ+QICE
- !
- ! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS
- ! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MI
- ! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURA
- ! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPR
- ! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE
- !
- !...NOTE THAT THE LIQ AND ICE MAY BE PRESENT IN PROPORTIONS SLIGHTLY DIF
- ! THAN SUGGESTED BY THE VALUE OF RATIO2...CHECK TO MAKE SURE THAT LIQ
- ! ICE CONCENTRATIONS ARE NOT REDUCED TO BELOW ZERO WHEN EVAPORATION/
- ! SUBLIMATION OCCURS...
- !
- IF(QTOT.GE.DQ)THEN
- DQICE=0.0
- DQLIQ=0.0
- QLIQ=QLIQ-(1.-RATIO2)*DQ
- IF(QLIQ.LT.0.)THEN
- DQICE=0.0-QLIQ
- QLIQ=0.0
- ENDIF
- QICE=QICE-RATIO2*DQ+DQICE
- IF(QICE.LT.0.)THEN
- DQLIQ=0.0-QICE
- QICE=0.0
- ENDIF
- QLIQ=QLIQ+DQLIQ
- QU=QS
- GOTO 96
- ELSE
- IF(RATIO2.LT.1.E-6)THEN
- RLL=XLV0-XLV1*T1
- ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
- RLL=XLS0-XLS1*T1
- ELSE
- RLL=RL
- ENDIF
- CCP=1005.7*(1.+0.89*QU)
- IF(QTOT.LT.1.E-10)THEN
- !
- !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
- T1=T1+RLL*(DQ/(1.+DQ))/CCP
- GOTO 96
- ELSE
- !
- !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURA
- ! THE TEMPERATURE IS GIVEN BY:
- T1=T1+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CCP
- QU=QU+QTOT
- QTOT=0.
- ENDIF
- QLIQ=0
- QICE=0.
- ENDIF
- 96 TU=T1
- QNEWLQ=(1.-RATIO2)*QNEW
- QNEWIC=RATIO2*QNEW
- IF(ITCNT.GT.10)PRINT*,'***** NUMBER OF ITERATIONS IN TPMIX =', &
- ITCNT
- END SUBROUTINE TPMIX
- !-----------------------------------------------------------------------
- SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,R1,RL, &
- EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- REAL, INTENT(IN ) :: P1,T1,Q1,R1,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,&
- BICE,CICE,DICE
- REAL, INTENT(INOUT) :: THT1
- REAL:: T00,P00,C1,C2,C3,C4,C5,EE,TLOG,TDPT,TSAT,THT,TFPT,TLOGIC, &
- TSATLQ,TSATIC
- DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834,&
- 0.278296,1.0723E-3/
- !
- ! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
- !
- IF(R1.LT.1.E-6)THEN
- EE=Q1*P1/(EP2+Q1)
- TLOG=ALOG(EE/ALIQ)
- TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
- TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
- THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
- THT1=THT*EXP((C1/TSAT-C2)*Q1*(1.+0.81*Q1))
- ELSEIF(ABS(R1-1.).LT.1.E-6)THEN
- EE=Q1*P1/(EP2+Q1)
- TLOG=ALOG(EE/AICE)
- TFPT=(CICE-DICE*TLOG)/(BICE-TLOG)
- THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
- TSAT=TFPT-(.182+1.13E-3*(TFPT-T00)-3.58E-4*(T1-T00))*(T1-TFPT)
- THT1=THT*EXP((C3/TSAT-C4)*Q1*(1.+0.81*Q1))
- ELSE
- EE=Q1*P1/(EP2+Q1)
- TLOG=ALOG(EE/ALIQ)
- TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
- TLOGIC=ALOG(EE/AICE)
- TFPT=(CICE-DICE*TLOGIC)/(BICE-TLOGIC)
- THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
- TSATLQ=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
- TSATIC=TFPT-(.182+1.13E-3*(TFPT-T00)-3.58E-4*(T1-T00))*(T1-TFPT)
- TSAT=R1*TSATIC+(1.-R1)*TSATLQ
- THT1=THT*EXP(RL*Q1*C5/TSAT*(1.+0.81*Q1))
- ENDIF
- END SUBROUTINE ENVIRTHT
- !-----------------------------------------------------------------------
- !************************* TPDD.FOR ************************************
- ! THIS SUBROUTINE ITERATIVELY EXTRACTS TEMPERATURE FROM EQUIVALENT *
- ! POTENTIAL TEMP. IT IS DESIGNED FOR USE WITH DOWNDRAFT CALCULATIONS.
- ! IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, PARCEL *
- ! TEMP, SPECIFIC HUMIDITY, AND LIQUID WATER CONTENT ARE ITERATIVELY *
- ! CALCULATED. *
- !***********************************************************************
- FUNCTION TPDD(P,THTED,TGS,RS,RD,RH,XLV0,XLV1, &
- EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- REAL, INTENT(IN ) :: XLV0,XLV1
- REAL, INTENT(IN ) :: P,THTED,TGS,RD,RH,EP2,ALIQ,BLIQ, &
- CLIQ,DLIQ,AICE,BICE,CICE,DICE
- REAL, INTENT(INOUT) :: RS
- REAL :: TPDD,ES,PI,THTGS,F0,T1,T0,CCP,F1,DT,RL,DSSDT,T1RH,RSRH
- INTEGER :: ITCNT
- !-----------------------------------------------------------------------
- ES=ALIQ*EXP((BLIQ*TGS-CLIQ)/(TGS-DLIQ))
- RS=EP2*ES/(P-ES)
- PI=(1.E5/P)**(0.2854*(1.-0.28*RS))
- THTGS=TGS*PI*EXP((3374.6525/TGS-2.5403)*RS*(1.+0.81*RS))
- F0=THTGS-THTED
- T1=TGS-0.5*F0
- T0=TGS
- CCP=1005.7
- !
- !...ITERATE TO FIND WET-BULB TEMPERATURE...
- !
- ITCNT=0
- 90 ES=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
- RS=EP2*ES/(P-ES)
- PI=(1.E5/P)**(0.2854*(1.-0.28*RS))
- THTGS=T1*PI*EXP((3374.6525/T1-2.5403)*RS*(1.+0.81*RS))
- F1=THTGS-THTED
- IF(ABS(F1).LT.0.05)GOTO 50
- ITCNT=ITCNT+1
- IF(ITCNT.GT.10)GOTO 50
- DT=F1*(T1-T0)/(F1-F0)
- T0=T1
- F0=F1
- T1=T1-DT
- GOTO 90
- 50 RL=XLV0-XLV1*T1
- !
- !...IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, ESTIMATE THE
- ! TEMPERATURE AND MIXING RATIO WHICH WILL YIELD THE APPROPRIATE VALUE.
- !
- IF(RH.EQ.1.)GOTO 110
- DSSDT=(CLIQ-BLIQ*DLIQ)/((T1-DLIQ)*(T1-DLIQ))
- DT=RL*RS*(1.-RH)/(CCP+RL*RH*RS*DSSDT)
- T1RH=T1+DT
- ES=RH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
- RSRH=EP2*ES/(P-ES)
- !
- !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
- !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
- !
- IF(RSRH.LT.RD)THEN
- RSRH=RD
- T1RH=T1+(RS-RSRH)*RL/CCP
- ENDIF
- T1=T1RH
- RS=RSRH
- 110 TPDD=T1
- IF(ITCNT.GT.10)PRINT*,'***** NUMBER OF ITERATIONS IN TPDD = ', &
- ITCNT
- END FUNCTION TPDD
- !====================================================================
- SUBROUTINE kfinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, &
- RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS, &
- P_FIRST_SCALAR,restart,allowed_to_read, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte )
- !--------------------------------------------------------------------
- IMPLICIT NONE
- !--------------------------------------------------------------------
- LOGICAL , INTENT(IN) :: restart, allowed_to_read
- INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte
- INTEGER , INTENT(IN) :: P_QI,P_QS,P_FIRST_SCALAR
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
- RTHCUTEN, &
- RQVCUTEN, &
- RQCCUTEN, &
- RQRCUTEN, &
- RQICUTEN, &
- RQSCUTEN
- REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG
- REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA
- INTEGER :: i, j, k, itf, jtf, ktf
-
- jtf=min0(jte,jde-1)
- ktf=min0(kte,kde-1)
- itf=min0(ite,ide-1)
-
- IF(.not.restart)THEN
- DO j=jts,jtf
- DO k=kts,ktf
- DO i=its,itf
- RTHCUTEN(i,k,j)=0.
- RQVCUTEN(i,k,j)=0.
- RQCCUTEN(i,k,j)=0.
- RQRCUTEN(i,k,j)=0.
- ENDDO
- ENDDO
- ENDDO
- IF (P_QI .ge. P_FIRST_SCALAR) THEN
- DO j=jts,jtf
- DO k=kts,ktf
- DO i=its,itf
- RQICUTEN(i,k,j)=0.
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- IF (P_QS .ge. P_FIRST_SCALAR) THEN
- DO j=jts,jtf
- DO k=kts,ktf
- DO i=its,itf
- RQSCUTEN(i,k,j)=0.
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- DO j=jts,jtf
- DO i=its,itf
- NCA(i,j)=-100.
- ENDDO
- ENDDO
- DO j=jts,jtf
- DO k=kts,ktf
- DO i=its,itf
- W0AVG(i,k,j)=0.
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- END SUBROUTINE kfinit
- !-------------------------------------------------------
- END MODULE module_cu_kf