/wrfv2_fire/phys/module_mp_HWRF.F
FORTRAN Legacy | 2653 lines | 1168 code | 55 blank | 1430 comment | 0 complexity | d9d0c6dbc1e2fdaab2f34f4034d6d2c2 MD5 | raw file
Possible License(s): AGPL-1.0
- !WRF:MODEL_MP:PHYSICS
- !
- MODULE module_mp_HWRF
- !-----------------------------------------------------------------------
- !-- The following changes were made on 24 July 2006.
- ! (1) All known version 2.1 dependencies were removed from the
- ! operational WRF NMM model code (search for "!HWRF")
- ! (2) Incorporated code changes from the GFDL model (search for "!GFDL")
- !-----------------------------------------------------------------------
- REAL,PRIVATE,SAVE :: ABFR, CBFR, CIACW, CIACR, C_N0r0, &
- & CN0r0, CN0r_DMRmin, CN0r_DMRmax, CRACW, CRAUT, ESW0, &
- & RFmax, RQR_DR1, RQR_DR2, RQR_DR3, RQR_DRmin, &
- & RQR_DRmax, RR_DRmin, RR_DR1, RR_DR2, RR_DR3, RR_DRmax
- !
- INTEGER, PRIVATE,PARAMETER :: MY_T1=1, MY_T2=35
- REAL,PRIVATE,DIMENSION(MY_T1:MY_T2),SAVE :: MY_GROWTH
- !
- REAL, PRIVATE,PARAMETER :: DMImin=.05e-3, DMImax=1.e-3, &
- & DelDMI=1.e-6,XMImin=1.e6*DMImin
- INTEGER, PUBLIC,PARAMETER :: XMImax=1.e6*DMImax, XMIexp=.0536, &
- & MDImin=XMImin, MDImax=XMImax
- REAL, PRIVATE,DIMENSION(MDImin:MDImax) :: &
- & ACCRI,SDENS,VSNOWI,VENTI1,VENTI2
- !
- REAL, PRIVATE,PARAMETER :: DMRmin=.05e-3, DMRmax=.45e-3, &
- & DelDMR=1.e-6,XMRmin=1.e6*DMRmin, XMRmax=1.e6*DMRmax
- INTEGER, PRIVATE,PARAMETER :: MDRmin=XMRmin, MDRmax=XMRmax
- REAL, PRIVATE,DIMENSION(MDRmin:MDRmax):: &
- & ACCRR,MASSR,RRATE,VRAIN,VENTR1,VENTR2
- !
- INTEGER, PRIVATE,PARAMETER :: Nrime=40
- REAL, DIMENSION(2:9,0:Nrime),PRIVATE,SAVE :: VEL_RF
- !
- INTEGER,PARAMETER :: NX=7501
- REAL, PARAMETER :: XMIN=180.0,XMAX=330.0
- REAL, DIMENSION(NX),PRIVATE,SAVE :: TBPVS,TBPVS0
- REAL, PRIVATE,SAVE :: C1XPVS0,C2XPVS0,C1XPVS,C2XPVS
- !
- REAL, PRIVATE,PARAMETER :: &
- !--- Physical constants follow:
- & CP=1004.6, EPSQ=1.E-12, GRAV=9.806, RHOL=1000., RD=287.04 &
- & ,RV=461.5, T0C=273.15, XLS=2.834E6 &
- !--- Derived physical constants follow:
- & ,EPS=RD/RV, EPS1=RV/RD-1., EPSQ1=1.001*EPSQ &
- & ,RCP=1./CP, RCPRV=RCP/RV, RGRAV=1./GRAV, RRHOL=1./RHOL &
- & ,XLS1=XLS*RCP, XLS2=XLS*XLS*RCPRV, XLS3=XLS*XLS/RV &
- !--- Constants specific to the parameterization follow:
- !--- CLIMIT/CLIMIT1 are lower limits for treating accumulated precipitation
- & ,CLIMIT=10.*EPSQ, CLIMIT1=-CLIMIT &
- & ,C1=1./3. &
- & ,DMR1=.1E-3, DMR2=.2E-3, DMR3=.32E-3 &
- & ,XMR1=1.e6*DMR1, XMR2=1.e6*DMR2, XMR3=1.e6*DMR3
- INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3
- REAL:: WC
- !
- ! ======================================================================
- !--- Important tunable parameters that are exported to other modules
- !GFDL * RHgrd - generic reference to the threshold relative humidity for
- !GFDL the onset of condensation
- !GFDL (new) * RHgrd_in - "RHgrd" for the inner domain
- !GFDL (new) * RHgrd_out - "RHgrd" for the outer domain
- !HWRF 6/11/2010 mod - use lower RHgrd_out for p < 850 hPa
- ! * T_ICE - temperature (C) threshold at which all remaining liquid water
- ! is glaciated to ice
- ! * T_ICE_init - maximum temperature (C) at which ice nucleation occurs
- ! * NLImax - maximum number concentrations (m**-3) of large ice (snow/graupel/sleet)
- ! * NLImin - minimum number concentrations (m**-3) of large ice (snow/graupel/sleet)
- ! * N0r0 - assumed intercept (m**-4) of rain drops if drop diameters are between 0.2 and 0.45 mm
- ! * N0rmin - minimum intercept (m**-4) for rain drops
- ! * NCW - number concentrations of cloud droplets (m**-3)
- ! * FLARGE1, FLARGE2 - number fraction of large ice to total (large+snow) ice
- ! at T>0C and in presence of sublimation (FLARGE1), otherwise in
- ! presence of ice saturated/supersaturated conditions
- ! * PRINT_diag - for extended model diagnostics (code currently commented out)
- ! * PRINT_err - for run-time prints when water budgets are not conserved (for debugging)
- ! ======================================================================
- REAL, PUBLIC,PARAMETER :: &
- ! & RHgrd=1. &
- & RHgrd_in=1. & !GFDL
- & ,RHgrd_out=0.975 & !GFDL
- & ,P_RHgrd_out=850.E2 & !HWRF 6/11/2010
- & ,T_ICE=-40. & !GFDL
- !GFDL & ,T_ICE=-30. &
- & ,T_ICEK=T0C+T_ICE &
- & ,T_ICE_init=-5. &
- & ,NLImax=5.E3 &
- & ,NLImin=1.E3 &
- & ,N0r0=8.E6 &
- & ,N0rmin=1.E4 &
- & ,NCW=60.E6 & !GFDL
- !HWRF & ,NCW=100.E6 &
- & ,FLARGE1=1. &
- & ,FLARGE2=.2
- ! LOGICAL, PARAMETER :: PRINT_diag=.FALSE. !GFDL
- LOGICAL, PARAMETER :: PRINT_err=.TRUE. !GFDL** => eventually set this to .false.
- !--- Other public variables passed to other routines:
- REAL,PUBLIC,SAVE :: QAUT0
- REAL, PUBLIC,DIMENSION(MDImin:MDImax) :: MASSI
- !
- !
- CONTAINS
- !-----------------------------------------------------------------------
- !-----------------------------------------------------------------------
- SUBROUTINE ETAMP_NEW_HWRF (itimestep,DT,DX,DY,GID,RAINNC,RAINNCV, & !GID
- & dz8w,rho_phy,p_phy,pi_phy,th_phy,qv,qt, & !gopal's doing
- & LOWLYR,SR, &
- & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
- & QC,QR,QI, &
- & ids,ide, jds,jde, kds,kde, &
- & ims,ime, jms,jme, kms,kme, &
- & its,ite, jts,jte, kts,kte )
- !HWRF SUBROUTINE ETAMP_NEW (itimestep,DT,DX,DY, &
- !HWRF & dz8w,rho_phy,p_phy,pi_phy,th_phy,qv,qc, &
- !HWRF & LOWLYR,SR, &
- !HWRF & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
- !HWRF & mp_restart_state,tbpvs_state,tbpvs0_state, &
- !HWRF & RAINNC,RAINNCV, &
- !HWRF & ids,ide, jds,jde, kds,kde, &
- !HWRF & ims,ime, jms,jme, kms,kme, &
- !HWRF & its,ite, jts,jte, kts,kte )
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- INTEGER, PARAMETER :: ITLO=-60, ITHI=40
- INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
- & ,IMS,IME,JMS,JME,KMS,KME &
- & ,ITS,ITE,JTS,JTE,KTS,KTE &
- & ,ITIMESTEP,GID ! GID gopal's doing
- REAL, INTENT(IN) :: DT,DX,DY
- REAL, INTENT(IN), DIMENSION(ims:ime, kms:kme, jms:jme):: &
- & dz8w,p_phy,pi_phy,rho_phy
- REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme):: &
- & th_phy,qv,qt,qc,qr,qi
- REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme ) :: &
- & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
- REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: &
- & RAINNC,RAINNCV
- REAL, INTENT(OUT), DIMENSION(ims:ime,jms:jme):: SR
- !
- !HWRF REAL,DIMENSION(*),INTENT(INOUT) :: MP_RESTART_STATE
- !
- !HWRF REAL,DIMENSION(nx),INTENT(INOUT) :: TBPVS_STATE,TBPVS0_STATE
- !
- INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR
- !-----------------------------------------------------------------------
- ! LOCAL VARS
- !-----------------------------------------------------------------------
- ! NSTATS,QMAX,QTOT are diagnostic vars
- INTEGER,DIMENSION(ITLO:ITHI,4) :: NSTATS
- REAL, DIMENSION(ITLO:ITHI,5) :: QMAX
- REAL, DIMENSION(ITLO:ITHI,22):: QTOT
- ! SOME VARS WILL BE USED FOR DATA ASSIMILATION (DON'T NEED THEM NOW).
- ! THEY ARE TREATED AS LOCAL VARS, BUT WILL BECOME STATE VARS IN THE
- ! FUTURE. SO, WE DECLARED THEM AS MEMORY SIZES FOR THE FUTURE USE
- ! TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related
- ! the microphysics scheme. Instead, they will be used by Eta precip
- ! assimilation.
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
- & TLATGS_PHY,TRAIN_PHY
- REAL, DIMENSION(ims:ime,jms:jme):: APREC,PREC,ACPREC
- REAL, DIMENSION(its:ite, kts:kte, jts:jte):: t_phy
- INTEGER :: I,J,K,KFLIP
- !
- !-----------------------------------------------------------------------
- !**********************************************************************
- !-----------------------------------------------------------------------
- !
- !HWRF MY_GROWTH(MY_T1:MY_T2)=MP_RESTART_STATE(MY_T1:MY_T2)
- !HWRF!
- !HWRF C1XPVS0=MP_RESTART_STATE(MY_T2+1)
- !HWRF C2XPVS0=MP_RESTART_STATE(MY_T2+2)
- !HWRF C1XPVS =MP_RESTART_STATE(MY_T2+3)
- !HWRF C2XPVS =MP_RESTART_STATE(MY_T2+4)
- !HWRF CIACW =MP_RESTART_STATE(MY_T2+5)
- !HWRF CIACR =MP_RESTART_STATE(MY_T2+6)
- !HWRF CRACW =MP_RESTART_STATE(MY_T2+7)
- !HWRF CRAUT =MP_RESTART_STATE(MY_T2+8)
- !HWRF!
- !HWRF TBPVS(1:NX) =TBPVS_STATE(1:NX)
- !HWRF TBPVS0(1:NX)=TBPVS0_STATE(1:NX)
- !
- DO j = jts,jte
- DO k = kts,kte
- DO i = its,ite
- t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
- qv(i,k,j)=qv(i,k,j)/(1.+qv(i,k,j)) !Convert to specific humidity
- ENDDO
- ENDDO
- ENDDO
- ! initial diagnostic variables and data assimilation vars
- ! (will need to delete this part in the future)
- DO k = 1,4
- DO i = ITLO,ITHI
- NSTATS(i,k)=0.
- ENDDO
- ENDDO
- DO k = 1,5
- DO i = ITLO,ITHI
- QMAX(i,k)=0.
- ENDDO
- ENDDO
- DO k = 1,22
- DO i = ITLO,ITHI
- QTOT(i,k)=0.
- ENDDO
- ENDDO
- ! initial data assimilation vars (will need to delete this part in the future)
- DO j = jts,jte
- DO k = kts,kte
- DO i = its,ite
- TLATGS_PHY (i,k,j)=0.
- TRAIN_PHY (i,k,j)=0.
- ENDDO
- ENDDO
- ENDDO
- DO j = jts,jte
- DO i = its,ite
- ACPREC(i,j)=0.
- APREC (i,j)=0.
- PREC (i,j)=0.
- SR (i,j)=0.
- ENDDO
- ENDDO
- !-- 6/11/2010: Update QT, F_ice, F_rain arrays
- DO j = jts,jte
- DO k = kts,kte
- DO i = its,ite
- QT(I,K,J)=QC(I,K,J)+QR(I,K,J)+QI(I,K,J)
- IF (QI(I,K,J) <= EPSQ) THEN
- F_ICE_PHY(I,K,J)=0.
- IF (T_PHY(I,K,J) < T_ICEK) F_ICE_PHY(I,K,J)=1.
- ELSE
- F_ICE_PHY(I,K,J)=MAX( 0., MIN(1., QI(I,K,J)/QT(I,K,J) ) )
- ENDIF
- IF (QR(I,K,J) <= EPSQ) THEN
- F_RAIN_PHY(I,K,J)=0.
- ELSE
- F_RAIN_PHY(I,K,J)=QR(I,K,J)/(QR(I,K,J)+QC(I,K,J))
- ENDIF
- ENDDO
- ENDDO
- ENDDO
- !-----------------------------------------------------------------------
- CALL EGCP01DRV(GID,DT,LOWLYR, &
- & APREC,PREC,ACPREC,SR,NSTATS,QMAX,QTOT, &
- & dz8w,rho_phy,qt,t_phy,qv,F_ICE_PHY,P_PHY, &
- & F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY,TRAIN_PHY, &
- & ids,ide, jds,jde, kds,kde, &
- & ims,ime, jms,jme, kms,kme, &
- & its,ite, jts,jte, kts,kte )
- !-----------------------------------------------------------------------
- DO j = jts,jte
- DO k = kts,kte
- DO i = its,ite
- th_phy(i,k,j) = t_phy(i,k,j)/pi_phy(i,k,j)
- qv(i,k,j)=qv(i,k,j)/(1.-qv(i,k,j)) !Convert to mixing ratio
- WC=qt(I,K,J)
- QI(I,K,J)=0.
- QR(I,K,J)=0.
- QC(I,K,J)=0.
- IF(F_ICE_PHY(I,K,J)>=1.)THEN
- QI(I,K,J)=WC
- ELSEIF(F_ICE_PHY(I,K,J)<=0.)THEN
- QC(I,K,J)=WC
- ELSE
- QI(I,K,J)=F_ICE_PHY(I,K,J)*WC
- QC(I,K,J)=WC-QI(I,K,J)
- ENDIF
- !
- IF(QC(I,K,J)>0..AND.F_RAIN_PHY(I,K,J)>0.)THEN
- IF(F_RAIN_PHY(I,K,J).GE.1.)THEN
- QR(I,K,J)=QC(I,K,J)
- QC(I,K,J)=0.
- ELSE
- QR(I,K,J)=F_RAIN_PHY(I,K,J)*QC(I,K,J)
- QC(I,K,J)=QC(I,K,J)-QR(I,K,J)
- ENDIF
- endif
- ENDDO
- ENDDO
- ENDDO
- !
- ! update rain (from m to mm)
- DO j=jts,jte
- DO i=its,ite
- RAINNC(i,j)=APREC(i,j)*1000.+RAINNC(i,j)
- RAINNCV(i,j)=APREC(i,j)*1000.
- ENDDO
- ENDDO
- !
- !HWRF MP_RESTART_STATE(MY_T1:MY_T2)=MY_GROWTH(MY_T1:MY_T2)
- !HWRF MP_RESTART_STATE(MY_T2+1)=C1XPVS0
- !HWRF MP_RESTART_STATE(MY_T2+2)=C2XPVS0
- !HWRF MP_RESTART_STATE(MY_T2+3)=C1XPVS
- !HWRF MP_RESTART_STATE(MY_T2+4)=C2XPVS
- !HWRF MP_RESTART_STATE(MY_T2+5)=CIACW
- !HWRF MP_RESTART_STATE(MY_T2+6)=CIACR
- !HWRF MP_RESTART_STATE(MY_T2+7)=CRACW
- !HWRF MP_RESTART_STATE(MY_T2+8)=CRAUT
- !HWRF!
- !HWRF TBPVS_STATE(1:NX) =TBPVS(1:NX)
- !HWRF TBPVS0_STATE(1:NX)=TBPVS0(1:NX)
- !-----------------------------------------------------------------------
- END SUBROUTINE ETAMP_NEW_HWRF
- !-----------------------------------------------------------------------
- SUBROUTINE EGCP01DRV(GID, & !GID gopal's doing
- & DTPH,LOWLYR,APREC,PREC,ACPREC,SR, &
- & NSTATS,QMAX,QTOT, &
- & dz8w,RHO_PHY,CWM_PHY,T_PHY,Q_PHY,F_ICE_PHY,P_PHY, &
- & F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY,TRAIN_PHY, &
- & ids,ide, jds,jde, kds,kde, &
- & ims,ime, jms,jme, kms,kme, &
- & its,ite, jts,jte, kts,kte)
- !-----------------------------------------------------------------------
- ! DTPH Physics time step (s)
- ! CWM_PHY (qt) Mixing ratio of the total condensate. kg/kg
- ! Q_PHY Mixing ratio of water vapor. kg/kg
- ! F_RAIN_PHY Fraction of rain.
- ! F_ICE_PHY Fraction of ice.
- ! F_RIMEF_PHY Mass ratio of rimed ice (rime factor).
- !
- !TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related the
- !micrphysics sechme. Instead, they will be used by Eta precip assimilation.
- !
- !NSTATS,QMAX,QTOT are used for diagnosis purposes.
- !
- !-----------------------------------------------------------------------
- !--- Variables APREC,PREC,ACPREC,SR are calculated for precip assimilation
- ! and/or ZHAO's scheme in Eta and are not required by this microphysics
- ! scheme itself.
- !--- NSTATS,QMAX,QTOT are used for diagnosis purposes only. They will be
- ! printed out when PRINT_diag is true.
- !
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- !
- INTEGER, PARAMETER :: ITLO=-60, ITHI=40
- ! VARIABLES PASSED IN/OUT
- INTEGER,INTENT(IN ) :: ids,ide, jds,jde, kds,kde &
- & ,ims,ime, jms,jme, kms,kme &
- & ,its,ite, jts,jte, kts,kte
- INTEGER,INTENT(IN ) :: GID ! grid%id gopal's doing
- REAL,INTENT(IN) :: DTPH
- INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR
- INTEGER,DIMENSION(ITLO:ITHI,4),INTENT(INOUT) :: NSTATS
- REAL,DIMENSION(ITLO:ITHI,5),INTENT(INOUT) :: QMAX
- REAL,DIMENSION(ITLO:ITHI,22),INTENT(INOUT) :: QTOT
- REAL,DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: &
- & APREC,PREC,ACPREC,SR
- REAL,DIMENSION( its:ite, kts:kte, jts:jte ),INTENT(INOUT) :: t_phy
- REAL,DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: &
- & dz8w,P_PHY,RHO_PHY
- REAL,DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(INOUT) :: &
- & CWM_PHY, F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY &
- & ,Q_PHY,TRAIN_PHY
- !
- !-----------------------------------------------------------------------
- !LOCAL VARIABLES
- !-----------------------------------------------------------------------
- !
- !HWRF - Below are directives in the operational code that have been removed,
- ! where "TEMP_DEX" has been replaced with "I,J,L" and "TEMP_DIMS" has
- ! been replaced with "its:ite,jts:jte,kts:kte"
- !HWRF#define CACHE_FRIENDLY_MP_ETANEW
- !HWRF#ifdef CACHE_FRIENDLY_MP_ETANEW
- !HWRF# define TEMP_DIMS kts:kte,its:ite,jts:jte
- !HWRF# define TEMP_DEX L,I,J
- !HWRF#else
- !HWRF# define TEMP_DIMS its:ite,jts:jte,kts:kte
- !HWRF# define TEMP_DEX I,J,L
- !HWRF#endif
- !HWRF!
- INTEGER :: LSFC,I,J,I_index,J_index,L,K,KFLIP
- !HWRF REAL,DIMENSION(TEMP_DIMS) :: CWM,T,Q,TRAIN,TLATGS,P
- REAL,DIMENSION(its:ite,jts:jte,kts:kte) :: &
- & CWM,T,Q,TRAIN,TLATGS,P
- REAL,DIMENSION(kts:kte,its:ite,jts:jte) :: F_ice,F_rain,F_RimeF
- INTEGER,DIMENSION(its:ite,jts:jte) :: LMH
- REAL :: TC,WC,QI,QR,QW,Fice,Frain,DUM,ASNOW,ARAIN
- REAL,DIMENSION(kts:kte) :: P_col,Q_col,T_col,QV_col,WC_col, &
- RimeF_col,QI_col,QR_col,QW_col, THICK_col, RHC_col, DPCOL !GFDL
- REAL,DIMENSION(2) :: PRECtot,PRECmax
- !-----------------------------------------------------------------------
- !
- DO J=JTS,JTE
- DO I=ITS,ITE
- LMH(I,J) = KTE-LOWLYR(I,J)+1
- ENDDO
- ENDDO
- DO 98 J=JTS,JTE
- DO 98 I=ITS,ITE
- DO L=KTS,KTE
- KFLIP=KTE+1-L
- CWM(I,J,L)=CWM_PHY(I,KFLIP,J)
- T(I,J,L)=T_PHY(I,KFLIP,J)
- Q(I,J,L)=Q_PHY(I,KFLIP,J)
- P(I,J,L)=P_PHY(I,KFLIP,J)
- TLATGS(I,J,L)=TLATGS_PHY(I,KFLIP,J)
- TRAIN(I,J,L)=TRAIN_PHY(I,KFLIP,J)
- F_ice(L,I,J)=F_ice_PHY(I,KFLIP,J)
- F_rain(L,I,J)=F_rain_PHY(I,KFLIP,J)
- F_RimeF(L,I,J)=F_RimeF_PHY(I,KFLIP,J)
- ENDDO
- 98 CONTINUE
-
- DO 100 J=JTS,JTE
- DO 100 I=ITS,ITE
- LSFC=LMH(I,J) ! "L" of surface
- !
- DO K=KTS,KTE
- KFLIP=KTE+1-K
- DPCOL(K)=RHO_PHY(I,KFLIP,J)*GRAV*dz8w(I,KFLIP,J)
- ENDDO
- !
- !
- !--- Initialize column data (1D arrays)
- !
- IF (CWM(I,J,1) .LE. EPSQ) CWM(I,J,1)=EPSQ
- F_ice(1,I,J)=1.
- F_rain(1,I,J)=0.
- F_RimeF(1,I,J)=1.
- DO L=1,LSFC
- !
- !--- Pressure (Pa) = (Psfc-Ptop)*(ETA/ETA_sfc)+Ptop
- !
- P_col(L)=P(I,J,L)
- !
- !--- Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
- !
- THICK_col(L)=DPCOL(L)*RGRAV
- T_col(L)=T(I,J,L)
- TC=T_col(L)-T0C
- QV_col(L)=max(EPSQ, Q(I,J,L))
- IF (CWM(I,J,L) .LE. EPSQ1) THEN
- WC_col(L)=0.
- IF (TC .LT. T_ICE) THEN
- F_ice(L,I,J)=1.
- ELSE
- F_ice(L,I,J)=0.
- ENDIF
- F_rain(L,I,J)=0.
- F_RimeF(L,I,J)=1.
- ELSE
- WC_col(L)=CWM(I,J,L)
- ENDIF
- !
- !--- Determine composition of condensate in terms of
- ! cloud water, ice, & rain
- !
- WC=WC_col(L)
- QI=0.
- QR=0.
- QW=0.
- Fice=F_ice(L,I,J)
- Frain=F_rain(L,I,J)
- IF (Fice .GE. 1.) THEN
- QI=WC
- ELSE IF (Fice .LE. 0.) THEN
- QW=WC
- ELSE
- QI=Fice*WC
- QW=WC-QI
- ENDIF
- IF (QW.GT.0. .AND. Frain.GT.0.) THEN
- IF (Frain .GE. 1.) THEN
- QR=QW
- QW=0.
- ELSE
- QR=Frain*QW
- QW=QW-QR
- ENDIF
- ENDIF
- RimeF_col(L)=F_RimeF(L,I,J)
- QI_col(L)=QI
- QR_col(L)=QR
- QW_col(L)=QW
- !GFDL => New. Added RHC_col to allow for height- and grid-dependent values for
- !GFDL the relative humidity threshold for condensation ("RHgrd")
- !6/11/2010 mod - Use lower RHgrd_out threshold for < 850 hPa
- !------------------------------------------------------------
- IF(GID .EQ. 1 .AND. P_col(L)<P_RHgrd_out) THEN ! gopal's doing based on GFDL
- RHC_col(L)=RHgrd_out
- ELSE
- RHC_col(L)=RHgrd_in
- ENDIF
- !------------------------------------------------------------
- ENDDO
- !
- !#######################################################################
- !
- !--- Perform the microphysical calculations in this column
- !
- I_index=I
- J_index=J
- CALL EGCP01COLUMN ( ARAIN, ASNOW, DTPH, I_index, J_index, LSFC, &
- & P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col, &
- & THICK_col, WC_col, RHC_col, KTS,KTE,NSTATS,QMAX,QTOT ) !GFDL
- !
- !#######################################################################
- !
- !
- !--- Update storage arrays
- !
- DO L=1,LSFC
- TRAIN(I,J,L)=(T_col(L)-T(I,J,L))/DTPH
- TLATGS(I,J,L)=T_col(L)-T(I,J,L)
- T(I,J,L)=T_col(L)
- Q(I,J,L)=QV_col(L)
- CWM(I,J,L)=WC_col(L)
- !
- !--- REAL*4 array storage
- !
- F_RimeF(L,I,J)=MAX(1., RimeF_col(L))
- IF (QI_col(L) .LE. EPSQ) THEN
- F_ice(L,I,J)=0.
- IF (T_col(L) .LT. T_ICEK) F_ice(L,I,J)=1.
- ELSE
- F_ice(L,I,J)=MAX( 0., MIN(1., QI_col(L)/WC_col(L)) )
- ENDIF
- IF (QR_col(L) .LE. EPSQ) THEN
- DUM=0
- ELSE
- DUM=QR_col(L)/(QR_col(L)+QW_col(L))
- ENDIF
- F_rain(L,I,J)=DUM
- !
- ENDDO
- !
- !--- Update accumulated precipitation statistics
- !
- !--- Surface precipitation statistics; SR is fraction of surface
- ! precipitation (if >0) associated with snow
- !
- APREC(I,J)=(ARAIN+ASNOW)*RRHOL ! Accumulated surface precip (depth in m) !<--- Ying
- PREC(I,J)=PREC(I,J)+APREC(I,J)
- ACPREC(I,J)=ACPREC(I,J)+APREC(I,J)
- IF(APREC(I,J) .LT. 1.E-8) THEN
- SR(I,J)=0.
- ELSE
- SR(I,J)=RRHOL*ASNOW/APREC(I,J)
- ENDIF
- ! !
- ! !--- Debug statistics
- ! !
- ! IF (PRINT_diag) THEN
- ! PRECtot(1)=PRECtot(1)+ARAIN
- ! PRECtot(2)=PRECtot(2)+ASNOW
- ! PRECmax(1)=MAX(PRECmax(1), ARAIN)
- ! PRECmax(2)=MAX(PRECmax(2), ASNOW)
- ! ENDIF
- !#######################################################################
- !#######################################################################
- !
- 100 CONTINUE ! End "I" & "J" loops
- DO 101 J=JTS,JTE
- DO 101 I=ITS,ITE
- DO L=KTS,KTE
- KFLIP=KTE+1-L
- CWM_PHY(I,KFLIP,J)=CWM(I,J,L)
- T_PHY(I,KFLIP,J)=T(I,J,L)
- Q_PHY(I,KFLIP,J)=Q(I,J,L)
- TLATGS_PHY(I,KFLIP,J)=TLATGS(I,J,L)
- TRAIN_PHY(I,KFLIP,J)=TRAIN(I,J,L)
- F_ice_PHY(I,KFLIP,J)=F_ice(L,I,J)
- F_rain_PHY(I,KFLIP,J)=F_rain(L,I,J)
- F_RimeF_PHY(I,KFLIP,J)=F_RimeF(L,I,J)
- ENDDO
- 101 CONTINUE
- END SUBROUTINE EGCP01DRV
- !
- !
- !###############################################################################
- ! ***** VERSION OF MICROPHYSICS DESIGNED FOR HIGHER RESOLUTION MESO ETA MODEL
- ! (1) Represents sedimentation by preserving a portion of the precipitation
- ! through top-down integration from cloud-top. Modified procedure to
- ! Zhao and Carr (1997).
- ! (2) Microphysical equations are modified to be less sensitive to time
- ! steps by use of Clausius-Clapeyron equation to account for changes in
- ! saturation mixing ratios in response to latent heating/cooling.
- ! (3) Prevent spurious temperature oscillations across 0C due to
- ! microphysics.
- ! (4) Uses lookup tables for: calculating two different ventilation
- ! coefficients in condensation and deposition processes; accretion of
- ! cloud water by precipitation; precipitation mass; precipitation rate
- ! (and mass-weighted precipitation fall speeds).
- ! (5) Assumes temperature-dependent variation in mean diameter of large ice
- ! (Houze et al., 1979; Ryan et al., 1996).
- ! -> 8/22/01: This relationship has been extended to colder temperatures
- ! to parameterize smaller large-ice particles down to mean sizes of MDImin,
- ! which is 50 microns reached at -55.9C.
- ! (6) Attempts to differentiate growth of large and small ice, mainly for
- ! improved transition from thin cirrus to thick, precipitating ice
- ! anvils.
- ! -> 8/22/01: This feature has been diminished by effectively adjusting to
- ! ice saturation during depositional growth at temperatures colder than
- ! -10C. Ice sublimation is calculated more explicitly. The logic is
- ! that sources of are either poorly understood (e.g., nucleation for NWP)
- ! or are not represented in the Eta model (e.g., detrainment of ice from
- ! convection). Otherwise the model is too wet compared to the radiosonde
- ! observations based on 1 Feb - 18 March 2001 retrospective runs.
- ! (7) Top-down integration also attempts to treat mixed-phase processes,
- ! allowing a mixture of ice and water. Based on numerous observational
- ! studies, ice growth is based on nucleation at cloud top &
- ! subsequent growth by vapor deposition and riming as the ice particles
- ! fall through the cloud. Effective nucleation rates are a function
- ! of ice supersaturation following Meyers et al. (JAM, 1992).
- ! -> 8/22/01: The simulated relative humidities were far too moist compared
- ! to the rawinsonde observations. This feature has been substantially
- ! diminished, limited to a much narrower temperature range of 0 to -10C.
- ! (8) Depositional growth of newly nucleated ice is calculated for large time
- ! steps using Fig. 8 of Miller and Young (JAS, 1979), at 1 deg intervals
- ! using their ice crystal masses calculated after 600 s of growth in water
- ! saturated conditions. The growth rates are normalized by time step
- ! assuming 3D growth with time**1.5 following eq. (6.3) in Young (1993).
- ! -> 8/22/01: This feature has been effectively limited to 0 to -10C.
- ! (9) Ice precipitation rates can increase due to increase in response to
- ! cloud water riming due to (a) increased density & mass of the rimed
- ! ice, and (b) increased fall speeds of rimed ice.
- ! -> 8/22/01: This feature has been effectively limited to 0 to -10C.
- !###############################################################################
- !###############################################################################
- !
- SUBROUTINE EGCP01COLUMN ( ARAIN, ASNOW, DTPH, I_index, J_index, &
- & LSFC, P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col, &
- & THICK_col, WC_col, RHC_col, KTS,KTE,NSTATS,QMAX,QTOT) !GFDL
- !
- !###############################################################################
- !###############################################################################
- !
- !-------------------------------------------------------------------------------
- !----- NOTE: Code is currently set up w/o threading!
- !-------------------------------------------------------------------------------
- !$$$ SUBPROGRAM DOCUMENTATION BLOCK
- ! . . .
- ! SUBPROGRAM: Grid-scale microphysical processes - condensation & precipitation
- ! PRGRMMR: Ferrier ORG: W/NP22 DATE: 08-2001
- ! PRGRMMR: Jin (Modification for WRF structure)
- !-------------------------------------------------------------------------------
- ! ABSTRACT:
- ! * Merges original GSCOND & PRECPD subroutines.
- ! * Code has been substantially streamlined and restructured.
- ! * Exchange between water vapor & small cloud condensate is calculated using
- ! the original Asai (1965, J. Japan) algorithm. See also references to
- ! Yau and Austin (1979, JAS), Rutledge and Hobbs (1983, JAS), and Tao et al.
- ! (1989, MWR). This algorithm replaces the Sundqvist et al. (1989, MWR)
- ! parameterization.
- !-------------------------------------------------------------------------------
- !
- ! USAGE:
- ! * CALL EGCP01COLUMN FROM SUBROUTINE EGCP01DRV
- !
- ! INPUT ARGUMENT LIST:
- ! DTPH - physics time step (s)
- ! I_index - I index
- ! J_index - J index
- ! LSFC - Eta level of level above surface, ground
- ! P_col - vertical column of model pressure (Pa)
- ! QI_col - vertical column of model ice mixing ratio (kg/kg)
- ! QR_col - vertical column of model rain ratio (kg/kg)
- ! QV_col - vertical column of model water vapor specific humidity (kg/kg)
- ! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
- ! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
- ! T_col - vertical column of model temperature (deg K)
- ! THICK_col - vertical column of model mass thickness (density*height increment)
- ! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
- ! RHC_col - vertical column of threshold relative humidity for onset of condensation (ratio) !GFDL
- !
- !
- ! OUTPUT ARGUMENT LIST:
- ! ARAIN - accumulated rainfall at the surface (kg)
- ! ASNOW - accumulated snowfall at the surface (kg)
- ! QV_col - vertical column of model water vapor specific humidity (kg/kg)
- ! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
- ! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
- ! QI_col - vertical column of model ice mixing ratio (kg/kg)
- ! QR_col - vertical column of model rain ratio (kg/kg)
- ! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
- ! T_col - vertical column of model temperature (deg K)
- !
- ! OUTPUT FILES:
- ! NONE
- !
- ! Subprograms & Functions called:
- ! * Real Function CONDENSE - cloud water condensation
- ! * Real Function DEPOSIT - ice deposition (not sublimation)
- !
- ! UNIQUE: NONE
- !
- ! LIBRARY: NONE
- !
- ! COMMON BLOCKS:
- ! CMICRO_CONS - key constants initialized in GSMCONST
- ! CMICRO_STATS - accumulated and maximum statistics
- ! CMY_GROWTH - lookup table for growth of ice crystals in
- ! water saturated conditions (Miller & Young, 1979)
- ! IVENT_TABLES - lookup tables for ventilation effects of ice
- ! IACCR_TABLES - lookup tables for accretion rates of ice
- ! IMASS_TABLES - lookup tables for mass content of ice
- ! IRATE_TABLES - lookup tables for precipitation rates of ice
- ! IRIME_TABLES - lookup tables for increase in fall speed of rimed ice
- ! RVENT_TABLES - lookup tables for ventilation effects of rain
- ! RACCR_TABLES - lookup tables for accretion rates of rain
- ! RMASS_TABLES - lookup tables for mass content of rain
- ! RVELR_TABLES - lookup tables for fall speeds of rain
- ! RRATE_TABLES - lookup tables for precipitation rates of rain
- !
- ! ATTRIBUTES:
- ! LANGUAGE: FORTRAN 90
- ! MACHINE : IBM SP
- !
- !
- !-------------------------------------------------------------------------
- !--------------- Arrays & constants in argument list ---------------------
- !-------------------------------------------------------------------------
- !
- IMPLICIT NONE
- !
- INTEGER,INTENT(IN) :: KTS,KTE,I_index, J_index, LSFC
- REAL,INTENT(INOUT) :: ARAIN, ASNOW
- REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: P_col, QI_col,QR_col &
- & ,QV_col ,QW_col, RimeF_col, T_col, THICK_col, WC_col, RHC_col !GFDL
- !
- !-------------------------------------------------------------------------
- !-------------- Common blocks for microphysical statistics ---------------
- !-------------------------------------------------------------------------
- !
- !-------------------------------------------------------------------------
- !--------- Common blocks for constants initialized in GSMCONST ----------
- !
- INTEGER, PARAMETER :: ITLO=-60, ITHI=40
- INTEGER,INTENT(INOUT) :: NSTATS(ITLO:ITHI,4)
- REAL,INTENT(INOUT) :: QMAX(ITLO:ITHI,5),QTOT(ITLO:ITHI,22)
- !
- !-------------------------------------------------------------------------
- !--------------- Common blocks for various lookup tables -----------------
- !
- !--- Discretized growth rates of small ice crystals after their nucleation
- ! at 1 C intervals from -1 C to -35 C, based on calculations by Miller
- ! and Young (1979, JAS) after 600 s of growth. Resultant growth rates
- ! are multiplied by physics time step in GSMCONST.
- !
- !-------------------------------------------------------------------------
- !
- !--- Mean ice-particle diameters varying from 50 microns to 1000 microns
- ! (1 mm), assuming an exponential size distribution.
- !
- !---- Meaning of the following arrays:
- ! - mdiam - mean diameter (m)
- ! - VENTI1 - integrated quantity associated w/ ventilation effects
- ! (capacitance only) for calculating vapor deposition onto ice
- ! - VENTI2 - integrated quantity associated w/ ventilation effects
- ! (with fall speed) for calculating vapor deposition onto ice
- ! - ACCRI - integrated quantity associated w/ cloud water collection by ice
- ! - MASSI - integrated quantity associated w/ ice mass
- ! - VSNOWI - mass-weighted fall speed of snow (large ice), used to calculate
- ! precipitation rates
- !
- !
- !-------------------------------------------------------------------------
- !
- !--- VEL_RF - velocity increase of rimed particles as functions of crude
- ! particle size categories (at 0.1 mm intervals of mean ice particle
- ! sizes) and rime factor (different values of Rime Factor of 1.1**N,
- ! where N=0 to Nrime).
- !
- !-------------------------------------------------------------------------
- !
- !--- Mean rain drop diameters varying from 50 microns (0.05 mm) to 450 microns
- ! (0.45 mm), assuming an exponential size distribution.
- !
- !-------------------------------------------------------------------------
- !------- Key parameters, local variables, & important comments ---------
- !-----------------------------------------------------------------------
- !
- !--- TOLER => Tolerance or precision for accumulated precipitation
- !
- REAL, PARAMETER :: TOLER=5.E-7, C2=1./6., RHO0=1.194, Xratio=.025
- !
- !--- If BLEND=1:
- ! precipitation (large) ice amounts are estimated at each level as a
- ! blend of ice falling from the grid point above and the precip ice
- ! present at the start of the time step (see TOT_ICE below).
- !--- If BLEND=0:
- ! precipitation (large) ice amounts are estimated to be the precip
- ! ice present at the start of the time step.
- !
- !--- Extended to include sedimentation of rain on 2/5/01
- !
- REAL, PARAMETER :: BLEND=1.
- !
- !-----------------------------------------------------------------------
- !--- Local variables
- !-----------------------------------------------------------------------
- !
- REAL EMAIRI, N0r, NLICE, NSmICE, RHgrd !GFDL
- LOGICAL CLEAR, ICE_logical, DBG_logical, RAIN_logical
- INTEGER :: IDR,INDEX_MY,INDEXR,INDEXR1,INDEXS,IPASS,ITDX,IXRF, &
- & IXS,LBEF,L
- !
- REAL :: ABI,ABW,AIEVP,ARAINnew,ASNOWnew,BLDTRH,BUDGET, &
- & CREVP,DELI,DELR,DELT,DELV,DELW,DENOMF, &
- & DENOMI,DENOMW,DENOMWI,DIDEP, &
- & DIEVP,DIFFUS,DLI,DTPH,DTRHO,DUM,DUM1, &
- & DUM2,DWV0,DWVI,DWVR,DYNVIS,ESI,ESW,FIR,FLARGE,FLIMASS, &
- & FSMALL,FWR,FWS,GAMMAR,GAMMAS, &
- & PCOND,PIACR,PIACW,PIACWI,PIACWR,PICND,PIDEP,PIDEP_max, &
- & PIEVP,PILOSS,PIMLT,PP,PRACW,PRAUT,PREVP,PRLOSS, &
- & QI,QInew,QLICE,QR,QRnew,QSI,QSIgrd,QSInew,QSW,QSW0, &
- & QSWgrd,QSWnew,QT,QTICE,QTnew,QTRAIN,QV,QW,QW0,QWnew, &
- & RFACTOR,RHO,RIMEF,RIMEF1,RQR,RR,RRHO,SFACTOR, &
- & TC,TCC,TFACTOR,THERM_COND,THICK,TK,TK2,TNEW, &
- & TOT_ICE,TOT_ICEnew,TOT_RAIN,TOT_RAINnew, &
- & VEL_INC,VENTR,VENTIL,VENTIS,VRAIN1,VRAIN2,VRIMEF,VSNOW, &
- & WC,WCnew,WSgrd,WS,WSnew,WV,WVnew,WVQW, &
- & XLF,XLF1,XLI,XLV,XLV1,XLV2,XLIMASS,XRF,XSIMASS
- !
- !#######################################################################
- !########################## Begin Execution ############################
- !#######################################################################
- !
- !
- ARAIN=0. ! Accumulated rainfall into grid box from above (kg/m**2)
- ASNOW=0. ! Accumulated snowfall into grid box from above (kg/m**2)
- !
- !-----------------------------------------------------------------------
- !------------ Loop from top (L=1) to surface (L=LSFC) ------------------
- !-----------------------------------------------------------------------
- !
- DO 10 L=1,LSFC
- !--- Skip this level and go to the next lower level if no condensate
- ! and very low specific humidities
- !
- IF (QV_col(L).LE.EPSQ .AND. WC_col(L).LE.EPSQ) GO TO 10
- !
- !-----------------------------------------------------------------------
- !------------ Proceed with cloud microphysics calculations -------------
- !-----------------------------------------------------------------------
- !
- TK=T_col(L) ! Temperature (deg K)
- TC=TK-T0C ! Temperature (deg C)
- PP=P_col(L) ! Pressure (Pa)
- QV=QV_col(L) ! Specific humidity of water vapor (kg/kg)
- WV=QV/(1.-QV) ! Water vapor mixing ratio (kg/kg)
- WC=WC_col(L) ! Grid-scale mixing ratio of total condensate (water or ice; kg/kg)
- RHgrd=RHC_col(L) ! Threshold relative humidity for the onset of condensation
- !
- !-----------------------------------------------------------------------
- !--- Moisture variables below are mixing ratios & not specifc humidities
- !-----------------------------------------------------------------------
- !
- CLEAR=.TRUE.
- !
- !--- This check is to determine grid-scale saturation when no condensate is present
- !
- ESW=1000.*FPVS0(TK) ! Saturation vapor pressure w/r/t water
- QSW=EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water
- QSI = QSW ! Initialize variable
- WS=QSW ! General saturation mixing ratio (water/ice)
- IF (TC .LT. 0.) THEN
- ESI=1000.*FPVS(TK) ! Saturation vapor pressure w/r/t ice
- QSI=EPS*ESI/(PP-ESI) ! Saturation mixing ratio w/r/t water
- WS=QSI ! General saturation mixing ratio (water/ice)
- ENDIF
- !
- !--- Effective grid-scale Saturation mixing ratios
- !
- QSWgrd=RHgrd*QSW
- QSIgrd=RHgrd*QSI
- WSgrd=RHgrd*WS
- !
- !--- Check if air is subsaturated and w/o condensate
- !
- IF (WV.GT.WSgrd .OR. WC.GT.EPSQ) CLEAR=.FALSE.
- !
- !--- Check if any rain is falling into layer from above
- !
- IF (ARAIN .GT. CLIMIT) THEN
- CLEAR=.FALSE.
- ELSE
- ARAIN=0.
- ENDIF
- !
- !--- Check if any ice is falling into layer from above
- !
- !--- NOTE that "SNOW" in variable names is synonomous with
- ! large, precipitation ice particles
- !
- IF (ASNOW .GT. CLIMIT) THEN
- CLEAR=.FALSE.
- ELSE
- ASNOW=0.
- ENDIF
- !
- !-----------------------------------------------------------------------
- !-- Loop to the end if in clear, subsaturated air free of condensate ---
- !-----------------------------------------------------------------------
- !
- IF (CLEAR) GO TO 10
- !
- !-----------------------------------------------------------------------
- !--------- Initialize RHO, THICK & microphysical processes -------------
- !-----------------------------------------------------------------------
- !
- !
- !--- Virtual temperature, TV=T*(1./EPS-1)*Q, Q is specific humidity;
- ! (see pp. 63-65 in Fleagle & Businger, 1963)
- !
- RHO=PP/(RD*TK*(1.+EPS1*QV)) ! Air density (kg/m**3)
- RRHO=1./RHO ! Reciprocal of air density
- DTRHO=DTPH*RHO ! Time step * air density
- BLDTRH=BLEND*DTRHO ! Blend parameter * time step * air density
- THICK=THICK_col(L) ! Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
- !
- ARAINnew=0. ! Updated accumulated rainfall
- ASNOWnew=0. ! Updated accumulated snowfall
- QI=QI_col(L) ! Ice mixing ratio
- QInew=0. ! Updated ice mixing ratio
- QR=QR_col(L) ! Rain mixing ratio
- QRnew=0. ! Updated rain ratio
- QW=QW_col(L) ! Cloud water mixing ratio
- QWnew=0. ! Updated cloud water ratio
- !
- PCOND=0. ! Condensation (>0) or evaporation (<0) of cloud water (kg/kg)
- PIDEP=0. ! Deposition (>0) or sublimation (<0) of ice crystals (kg/kg)
- PIACW=0. ! Cloud water collection (riming) by precipitation ice (kg/kg; >0)
- PIACWI=0. ! Growth of precip ice by riming (kg/kg; >0)
- PIACWR=0. ! Shedding of accreted cloud water to form rain (kg/kg; >0)
- PIACR=0. ! Freezing of rain onto large ice at supercooled temps (kg/kg; >0)
- PICND=0. ! Condensation (>0) onto wet, melting ice (kg/kg)
- PIEVP=0. ! Evaporation (<0) from wet, melting ice (kg/kg)
- PIMLT=0. ! Melting ice (kg/kg; >0)
- PRAUT=0. ! Cloud water autoconversion to rain (kg/kg; >0)
- PRACW=0. ! Cloud water collection (accretion) by rain (kg/kg; >0)
- PREVP=0. ! Rain evaporation (kg/kg; <0)
- !
- !--- Double check input hydrometeor mixing ratios
- !
- ! DUM=WC-(QI+QW+QR)
- ! DUM1=ABS(DUM)
- ! DUM2=TOLER*MIN(WC, QI+QW+QR)
- ! IF (DUM1 .GT. DUM2) THEN
- ! WRITE(6,"(/2(a,i4),a,i2)") '{@ i=',I_index,' j=',J_index,
- ! & ' L=',L
- ! WRITE(6,"(4(a12,g11.4,1x))")
- ! & '{@ TCold=',TC,'P=',.01*PP,'DIFF=',DUM,'WCold=',WC,
- ! & '{@ QIold=',QI,'QWold=',QW,'QRold=',QR
- ! ENDIF
- !
- !***********************************************************************
- !*********** MAIN MICROPHYSICS CALCULATIONS NOW FOLLOW! ****************
- !***********************************************************************
- !
- !--- Calculate a few variables, which are used more than once below
- !
- !--- Latent heat of vaporization as a function of temperature from
- ! Bolton (1980, JAS)
- !
- XLV=3.148E6-2370*TK ! Latent heat of vaporization (Lv)
- XLF=XLS-XLV ! Latent heat of fusion (Lf)
- XLV1=XLV*RCP ! Lv/Cp
- XLF1=XLF*RCP ! Lf/Cp
- TK2=1./(TK*TK) ! 1./TK**2
- !GFDL XLV2=XLV*XLV*QSW*TK2/RV ! Lv**2*Qsw/(Rv*TK**2)
- XLV2=XLV*XLV*QSWgrd*TK2/RV ! Lv**2*QSWgrd/(Rv*TK**2) !GFDL
- DENOMW=1.+XLV2*RCP ! Denominator term, Clausius-Clapeyron correction
- !
- !--- Basic thermodynamic quantities
- ! * DYNVIS - dynamic viscosity [ kg/(m*s) ]
- ! * THERM_COND - thermal conductivity [ J/(m*s*K) ]
- ! * DIFFUS - diffusivity of water vapor [ m**2/s ]
- !
- TFACTOR=TK**1.5/(TK+120.)
- DYNVIS=1.496E-6*TFACTOR
- THERM_COND=2.116E-3*TFACTOR
- DIFFUS=8.794E-5*TK**1.81/PP
- !
- !--- Air resistance term for the fall speed of ice following the
- ! basic research by Heymsfield, Kajikawa, others
- !
- GAMMAS=(1.E5/PP)**C1
- !
- !--- Air resistance for rain fall speed (Beard, 1985, JAS, p.470)
- !
- GAMMAR=(RHO0/RHO)**.4
- !
- !----------------------------------------------------------------------
- !------------- IMPORTANT MICROPHYSICS DECISION TREE -----------------
- !----------------------------------------------------------------------
- !
- !--- Determine if conditions supporting ice are present
- !
- IF (TC.LT.0. .OR. QI.GT.EPSQ .OR. ASNOW.GT.CLIMIT) THEN
- ICE_logical=.TRUE.
- ELSE
- ICE_logical=.FALSE.
- QLICE=0.
- QTICE=0.
- ENDIF
- !
- !--- Determine if rain is present
- !
- RAIN_logical=.FALSE.
- IF (ARAIN.GT.CLIMIT .OR. QR.GT.EPSQ) RAIN_logical=.TRUE.
- !
- IF (ICE_logical) THEN
- !
- !--- IMPORTANT: Estimate time-averaged properties.
- !
- !---
- ! * FLARGE - ratio of number of large ice to total (large & small) ice
- ! * FSMALL - ratio of number of small ice crystals to large ice particles
- ! -> Small ice particles are assumed to have a mean diameter of 50 microns.
- ! * XSIMASS - used for calculating small ice mixing ratio
- !---
- ! * TOT_ICE - total mass (small & large) ice before microphysics,
- ! which is the sum of the total mass of large ice in the
- ! current layer and the input flux of ice from above
- ! * PILOSS - greatest loss (<0) of total (small & large) ice by
- ! sublimation, removing all of the ice falling from above
- ! and the ice within the layer
- ! * RimeF1 - Rime Factor, which is the mass ratio of total (unrimed & rimed)
- ! ice mass to the unrimed ice mass (>=1)
- ! * VrimeF - the velocity increase due to rime factor or melting (ratio, >=1)
- ! * VSNOW - Fall speed of rimed snow w/ air resistance correction
- ! * EMAIRI - equivalent mass of air associated layer and with fall of snow into layer
- ! * XLIMASS - used for calculating large ice mixing ratio
- ! * FLIMASS - mass fraction of large ice
- ! * QTICE - time-averaged mixing ratio of total ice
- ! * QLICE - time-averaged mixing ratio of large ice
- ! * NLICE - time-averaged number concentration of large ice
- ! * NSmICE - number concentration of small ice crystals at current level
- !---
- !--- Assumed number fraction of large ice particles to total (large & small)
- ! ice particles, which is based on a general impression of the literature.
- !
- WVQW=WV+QW ! Water vapor & cloud water
- !
- IF (TC.GE.0. .OR. WVQW.LT.QSIgrd) THEN
- !
- !--- Eliminate small ice particle contributions for melting & sublimation
- !
- FLARGE=FLARGE1
- ELSE
- !
- !--- Enhanced number of small ice particles during depositional growth
- ! (effective only when 0C > T >= T_ice [-10C] )
- !
- FLARGE=FLARGE2
- !
- !--- Larger number of small ice particles due to rime splintering
- !
- IF (TC.GE.-8. .AND. TC.LE.-3.) FLARGE=.5*FLARGE
- !
- ENDIF ! End IF (TC.GE.0. .OR. WVQW.LT.QSIgrd)
- !GFDL => turned on in GFDL code, but not here => FLARGE=1.0
- FSMALL=(1.-FLARGE)/FLARGE
- XSIMASS=RRHO*MASSI(MDImin)*FSMALL
- IF (QI.LE.EPSQ .AND. ASNOW.LE.CLIMIT) THEN
- INDEXS=MDImin
- TOT_ICE=0.
- PILOSS=0.
- RimeF1=1.
- VrimeF=1.
- VEL_INC=GAMMAS
- VSNOW=0.
- EMAIRI=THICK
- XLIMASS=RRHO*RimeF1*MASSI(INDEXS)
- FLIMASS=XLIMASS/(XLIMASS+XSIMASS)
- QLICE=0.
- QTICE=0.
- NLICE=0.
- NSmICE=0.
- ELSE
- !
- !--- For T<0C mean particle size follows Houze et al. (JAS, 1979, p. 160),
- ! converted from Fig. 5 plot of LAMDAs. Similar set of relationships
- ! also shown in Fig. 8 of Ryan (BAMS, 1996, p. 66).
- !
- DUM=XMImax*EXP(.0536*TC)
- INDEXS=MIN(MDImax, MAX(MDImin, INT(DUM) ) )
- TOT_ICE=THICK*QI+BLEND*ASNOW
- PILOSS=-TOT_ICE/THICK
- LBEF=MAX(1,L-1)
- DUM1=RimeF_col(LBEF)
- DUM2=RimeF_col(L)
- RimeF1=(DUM2*THICK*QI+DUM1*BLEND*ASNOW)/TOT_ICE
- RimeF1=MIN(RimeF1, RFmax)
- DO IPASS=0,1
- IF (RimeF1 .LE. 1.) THEN
- RimeF1=1.
- VrimeF=1.
- ELSE
- IXS=MAX(2, MIN(INDEXS/100, 9))
- XRF=10.492*ALOG(RimeF1)
- IXRF=MAX(0, MIN(INT(XRF), Nrime))
- IF (IXRF .GE. Nrime) THEN
- VrimeF=VEL_RF(IXS,Nrime)
- ELSE
- VrimeF=VEL_RF(IXS,IXRF)+(XRF-FLOAT(IXRF))* &
- & (VEL_RF(IXS,IXRF+1)-VEL_RF(IXS,IXRF))
- ENDIF
- ENDIF ! End IF (RimeF1 .LE. 1.)
- VEL_INC=GAMMAS*VrimeF
- VSNOW=VEL_INC*VSNOWI(INDEXS)
- EMAIRI=THICK+BLDTRH*VSNOW
- XLIMASS=RRHO*RimeF1*MASSI(INDEXS)
- FLIMASS=XLIMASS/(XLIMASS+XSIMASS)
- QTICE=TOT_ICE/EMAIRI
- QLICE=FLIMASS*QTICE
- NLICE=QLICE/XLIMASS
- NSmICE=Fsmall*NLICE
- !
- IF ( (NLICE.GE.NLImin .AND. NLICE.LE.NLImax) &
- & .OR. IPASS.EQ.1) THEN
- EXIT
- ELSE
- !
- !--- Reduce excessive accumulation of ice at upper levels
- ! associated with strong grid-resolved ascent
- !
- !--- Force NLICE to be between NLImin and NLImax
- !
- DUM=MAX(NLImin, MIN(NLImax, NLICE) )
- XLI=RHO*(QTICE/DUM-XSIMASS)/RimeF1
- IF (XLI .LE. MASSI(MDImin) ) THEN
- INDEXS=MDImin
- ELSE IF (XLI .LE. MASSI(450) ) THEN
- DLI=9.5885E5*XLI**.42066 ! DLI in microns
- INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
- ELSE IF (XLI .LE. MASSI(MDImax) ) THEN
- DLI=3.9751E6*XLI**.49870 ! DLI in microns
- INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
- ELSE
- INDEXS=MDImax
- !
- !--- 8/22/01: Increase density of large ice if maximum limits
- ! are reached for number concentration (NLImax) and mean size
- ! (MDImax). Done to increase fall out of ice.
- !
- IF (DUM .GE. NLImax) &
- & RimeF1=RHO*(QTICE/NLImax-XSIMASS)/MASSI(INDEXS)
- ENDIF ! End IF (XLI .LE. MASSI(MDImin) )
- ! WRITE(6,"(4(a12,g11.4,1x))")
- ! & '{$ TC=',TC,'P=',.01*PP,'NLICE=',NLICE,'DUM=',DUM,
- ! & '{$ XLI=',XLI,'INDEXS=',FLOAT(INDEXS),'RHO=',RHO,'QTICE=',QTICE,
- ! & '{$ XSIMASS=',XSIMASS,'RimeF1=',RimeF1
- ENDIF ! End IF ( (NLICE.GE.NLImin .AND. NLICE.LE.NLImax) ...
- ENDDO ! End DO IPASS=0,1
- ENDIF ! End IF (QI.LE.EPSQ .AND. ASNOW.LE.CLIMIT)
- ENDIF ! End IF (ICE_logical)
- !
- !----------------------------------------------------------------------
- !--------------- Calculate individual processes -----------------------
- !----------------------------------------------------------------------
- !
- !--- Cloud water autoconversion to rain and collection by rain
- !
- IF (QW.GT.EPSQ .AND. TC.GE.T_ICE) THEN
- !
- !--- QW0 could be modified based on land/sea properties,
- ! presence of convection, etc. This is why QAUT0 and CRAUT
- ! are passed into the subroutine as externally determined
- ! parameters. Can be changed in the future if desired.
- !
- QW0=QAUT0*RRHO
- PRAUT=MAX(0., QW-QW0)*CRAUT
- IF (QLICE .GT. EPSQ) THEN
- !
- !--- Collection of cloud water by large ice particles ("snow")
- ! PIACWI=PIACW for riming, PIACWI=0 for shedding
- !
- FWS=MIN(1., CIACW*VEL_INC*NLICE*ACCRI(INDEXS)/PP**C1)
- PIACW=FWS*QW
- IF (TC .LT. 0.) PIACWI=PIACW ! Large ice riming
- ENDIF ! End IF (QLICE .GT. EPSQ)
- ENDIF ! End IF (QW.GT.EPSQ .AND. TC.GE.T_ICE)
- !
- !----------------------------------------------------------------------
- !--- Loop around some of the ice-phase processes if no ice should be present
- !----------------------------------------------------------------------
- !
- IF (ICE_logical .EQV. .FALSE.) GO TO 20
- !
- !--- Now the pretzel logic of calculating ice deposition
- !
- IF (TC.LT.T_ICE .AND. (WV.GT.QSIgrd .OR. QW.GT.EPSQ)) THEN
- !
- !--- Adjust to ice saturation at T<T_ICE (-10C) if supersaturated.
- ! Sources of ice due to nucleation and convective detrainment are
- ! either poorly understood, poorly resolved at typical NWP
- ! resolutions, or are not represented (e.g., no detrained
- ! condensate in BMJ Cu scheme).
- !
- PCOND=-QW
- DUM1=TK+XLV1*PCOND ! Updated (dummy) temperature (deg K)
- DUM2=WV+QW ! Updated (dummy) water vapor mixing ratio
- DUM=1000.*FPVS(DUM1) ! Updated (dummy) saturation vapor pressure w/r/t ice
- DUM=RHgrd*EPS*DUM/(PP-DUM) ! Updated (dummy) saturation mixing ratio w/r/t ice
- IF (DUM2 .GT. DUM) PIDEP=DEPOSIT (PP, DUM1, DUM2, RHgrd) !GFDL
- DWVi=0. ! Used only for debugging
- !
- ELSE IF (TC .LT. 0.) THEN
- !
- !--- These quantities are handy for ice deposition/sublimation
- ! PIDEP_max - max deposition or minimum sublimation to ice saturation
- !
- !GFDL DENOMI=1.+XLS2*QSI*TK2
- !GFDL DWVi=MIN(WVQW,QSW)-QSI
- DENOMI=1.+XLS2*QSIgrd*TK2 !GFDL
- DWVi=MIN(WVQW,QSWgrd)-QSIgrd !GFDL
- PIDEP_max=MAX(PILOSS, DWVi/DENOMI)
- IF (QTICE .GT. 0.) THEN
- !
- !--- Calculate ice deposition/sublimation
- ! * SFACTOR - [VEL_INC**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5],
- ! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
- ! * Units: SFACTOR - s**.5/m ; ABI - m**2/s ; NLICE - m**-3 ;
- ! VENTIL, VENTIS - m**-2 ; VENTI1 - m ;
- ! VENTI2 - m**2/s**.5 ; DIDEP - unitless
- !
- SFACTOR=SQRT(VEL_INC)*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2 !GFDL
- ABI=1./(RHO*XLS3*QSI*TK2/THERM_COND+1./DIFFUS)
- !
- !--- VENTIL - Number concentration * ventilation factors for large ice
- !--- VENTIS - Number concentration * ventilation factors for small ice
- !
- !--- Variation in the number concentration of ice with time is not
- ! accounted for in these calculations (could be in the future).
- !
- VENTIL=(VENTI1(INDEXS)+SFACTOR*VENTI2(INDEXS))*NLICE
- VENTIS=(VENTI1(MDImin)+SFACTOR*VENTI2(MDImin))*NSmICE
- DIDEP=ABI*(VENTIL+VENTIS)*DTPH
- !
- !--- Account for change in water vapor supply w/ time
- !
- IF (DIDEP .GE. Xratio)then
- DIDEP=(1.-EXP(-DIDEP*DENOMI))/DENOMI
- endif
- IF (DWVi .GT. 0.) THEN
- PIDEP=MIN(DWVi*DIDEP, PIDEP_max)
- ELSE IF (DWVi .LT. 0.) THEN
- PIDEP=MAX(DWVi*DIDEP, PIDEP_max)
- ENDIF
- !
- !GFDL ELSE IF (WVQW.GT.QSI .AND. TC.LE.T_ICE_init) THEN
- ELSE IF (WVQW.GT.QSIgrd .AND. TC.LE.T_ICE_init) THEN !GFDL
- !
- !--- Ice nucleation in near water-saturated conditions. Ice crystal
- ! growth during time step calculated using Miller & Young (1979, JAS).
- !--- These deposition rates could drive conditions below water saturation,
- ! which is the basis of these calculations. Intended to approximate
- ! more complex & computationally intensive calculations.
- !
- INDEX_MY=MAX(MY_T1, MIN( INT(.5-TC), MY_T2 ) )
- !
- !--- DUM1 is the supersaturation w/r/t ice at water-saturated conditions
- !
- !--- DUM2 is the number of ice crystals nucleated at water-saturated
- ! conditions based on Meyers et al. (JAM, 1992).
- !
- !--- Prevent unrealistically large ice initiation (limited by PIDEP_max)
- ! if DUM2 values are increased in future experiments
- !
- DUM1=QSW/QSI-1.
- DUM2=1.E3*EXP(12.96*DUM1-.639)
- PIDEP=MIN(PIDEP_max, DUM2*MY_GROWTH(INDEX_MY)*RRHO)
- !
- ENDIF ! End IF (QTICE .GT. 0.)
- !
- ENDIF ! End IF (TC.LT.T_ICE .AND. (WV.GT.QSIgrd .OR. QW.GT.EPSQ))
- !
- !------------------------------------------------------------------------
- !
- 20 CONTINUE ! Jump here if conditions for ice are not present
- !
- !------------------------------------------------------------------------
- !
- !--- Cloud water condensation
- !
- IF (TC.GE.T_ICE .AND. (QW.GT.EPSQ .OR. WV.GT.QSWgrd)) THEN
- IF (PIACWI.EQ.0. .AND. PIDEP.EQ.0.) THEN
- PCOND=CONDENSE (PP, QW, TK, WV, RHgrd) !GFDL
- ELSE
- !
- !--- Modify cloud condensation in response to ice processes
- !
- DUM=XLV*QSWgrd*RCPRV*TK2
- DENOMWI=1.+XLS*DUM
- DENOMF=XLF*DUM
- DUM=MAX(0., PIDEP)
- PCOND=(WV-QSWgrd-DENOMWI*DUM-DENOMF*PIACWI)/DENOMW
- DUM1=-QW
- DUM2=PCOND-PIACW
- IF (DUM2 .LT. DUM1) THEN
- !
- !--- Limit cloud water sinks
- !
- DUM=DUM1/DUM2
- PCOND=DUM*PCOND
- PIACW=DUM*PIACW
- PIACWI=DUM*PIACWI
- ENDIF ! End IF (DUM2 .LT. DUM1)
- ENDIF ! End IF (PIACWI.EQ.0. .AND. PIDEP.EQ.0.)
- ENDIF ! End IF (TC.GE.T_ICE .AND. (QW.GT.EPSQ .OR. WV.GT.QSWgrd))
- !
- !--- Limit freezing of accreted rime to prevent temperature oscillations,
- ! a crude Schumann-Ludlam limit (p. 209 of Young, 1993).
- !
- TCC=TC+XLV1*PCOND+XLS1*PIDEP+XLF1*PIACWI
- IF (TCC .GT. 0.) THEN
- PIACWI=0.
- TCC=TC+XLV1*PCOND+XLS1*PIDEP
- ENDIF
- IF (TC.GT.0. .AND. TCC.GT.0. .AND. ICE_logical) THEN
- !
- !--- Calculate melting and evaporation/condensation
- ! * Units: SFACTOR - s**.5/m ; ABI - m**2/s ; NLICE - m**-3 ;
- ! VENTIL - m**-2 ; VENTI1 - m ;
- ! VENTI2 - m**2/s**.5 ; CIEVP - /s
- !
- SFACTOR=SQRT(VEL_INC)*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2 !GFDL
- VENTIL=NLICE*(VENTI1(INDEXS)+SFACTOR*VENTI2(INDEXS))
- AIEVP=VENTIL*DIFFUS*DTPH
- IF (AIEVP .LT. Xratio) THEN
- DIEVP=AIEVP
- ELSE
- DIEVP=1.-EXP(-AIEVP)
- ENDIF
- QSW0=EPS*ESW0/(PP-ESW0)
- !GFDL DWV0=MIN(WV,QSW)-QSW0
- DWV0=MIN(WV,QSWgrd)-QSW0*RHgrd !GFDL
- DUM=QW+PCOND
- !GFDL IF (WV.LT.QSW .AND. DUM.LE.EPSQ) THEN
- IF (WV.LT.QSWgrd .AND. DUM.LE.EPSQ) THEN !GFDL
- !
- !--- Evaporation from melting snow (sink of snow) or shedding
- ! of water condensed onto melting snow (source of rain)
- !
- DUM=DWV0*DIEVP
- PIEVP=MAX( MIN(0., DUM), PILOSS)
- PICND=MAX(0., DUM)
- ENDIF ! End IF (WV.LT.QSW .AND. DUM.LE.EPSQ)
- PIMLT=THERM_COND*TCC*VENTIL*RRHO*DTPH/XLF
- !
- !--- Limit melting to prevent temperature oscillations across 0C
- !
- DUM1=MAX( 0., (TCC+XLV1*PIEVP)/XLF1 )
- PIMLT=MIN(PIMLT, DUM1)
- !
- !--- Limit loss of snow by melting (>0) and evaporation
- !
- DUM=PIEVP-PIMLT
- IF (DUM .LT. PILOSS) THEN
- DUM1=PILOSS/DUM
- PIMLT=PIMLT*DUM1
- PIEVP=PIEVP*DUM1
- ENDIF ! End IF (DUM .GT. QTICE)
- ENDIF ! End IF (TC.GT.0. .AND. TCC.GT.0. .AND. ICE_logical)
- !
- !--- IMPORTANT: Estimate time-averaged properties.
- !
- ! * TOT_RAIN - total mass of rain before microphysics, which is the sum of
- ! the total mass of rain in the current layer and the input
- ! flux of rain from above
- ! * VRAIN1 - fall speed of rain into grid from above (with air resistance correction)
- ! * QTRAIN - time-averaged mixing ratio of rain (kg/kg)
- ! * PRLOSS - greatest loss (<0) of rain, removing all rain falling from
- ! above and the rain within the layer
- ! * RQR - rain content (kg/m**3)
- ! * INDEXR - mean size of rain drops to the nearest 1 micron in size
- ! * N0r - intercept of rain size distribution (typically 10**6 m**-4)
- !
- TOT_RAIN=0.
- VRAIN1=0.
- QTRAIN=0.
- PRLOSS=0.
- RQR=0.
- N0r=0.
- INDEXR=MDRmin
- INDEXR1=INDEXR ! For debugging only
- IF (RAIN_logical) THEN
- IF (ARAIN .LE. 0.) THEN
- INDEXR=MDRmin
- VRAIN1=0.
- ELSE
- !
- !--- INDEXR (related to mean diameter) & N0r could be modified
- ! by land/sea properties, presence of convection, etc.
- !
- !--- Rain rate normalized to a density of 1.194 kg/m**3
- !
- RR=ARAIN/(DTPH*GAMMAR)
- !
- IF (RR .LE. RR_DRmin) THEN
- !
- !--- Assume fixed mean diameter of rain (0.2 mm) for low rain rates,
- ! instead vary N0r with rain rate
- !
- INDEXR=MDRmin
- ELSE IF (RR .LE. RR_DR1) THEN
- !
- !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
- ! for mean diameters (Dr) between 0.05 and 0.10 mm:
- ! V(Dr)=5.6023e4*Dr**1.136, V in m/s and Dr in m
- ! RR = PI*1000.*N0r0*5.6023e4*Dr**(4+1.136) = 1.408e15*Dr**5.136,
- ! RR in kg/(m**2*s)
- ! Dr (m) = 1.123e-3*RR**.1947 -> Dr (microns) = 1.123e3*RR**.1947
- !
- INDEXR=INT( 1.123E3*RR**.1947 + .5 )
- INDEXR=MAX( MDRmin, MIN(INDEXR, MDR1) )
- ELSE IF (RR .LE. RR_DR2) THEN
- !
- !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
- ! for mean diameters (Dr) between 0.10 and 0.20 mm:
- ! V(Dr)=1.0867e4*Dr**.958, V in m/s and Dr in m
- ! RR = PI*1000.*N0r0*1.0867e4*Dr**(4+.958) = 2.731e14*Dr**4.958,
- ! RR in kg/(m**2*s)
- ! Dr (m) = 1.225e-3*RR**.2017 -> Dr (microns) = 1.225e3*RR**.2017
- !
- INDEXR=INT( 1.225E3*RR**.2017 + .5 )
- INDEXR=MAX( MDR1, MIN(INDEXR, MDR2) )
- ELSE IF (RR .LE. RR_DR3) THEN
- !
- !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
- ! for mean diameters (Dr) between 0.20 and 0.32 mm:
- ! V(Dr)=2831.*Dr**.80, V in m/s and Dr in m
- ! RR = PI*1000.*N0r0*2831.*Dr**(4+.80) = 7.115e13*Dr**4.80,
- ! RR in kg/(m**2*s)
- ! Dr (m) = 1.3006e-3*RR**.2083 -> Dr (microns) = 1.3006e3*RR**.2083
- !
- INDEXR=INT( 1.3006E3*RR**.2083 + .5 )
- INDEXR=MAX( MDR2, MIN(INDEXR, MDR3) )
- ELSE IF (RR .LE. RR_DRmax) THEN
- !
- !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
- ! for mean diameters (Dr) between 0.32 and 0.45 mm:
- ! V(Dr)=944.8*Dr**.6636, V in m/s and Dr in m
- ! RR = PI*1000.*N0r0*944.8*Dr**(4+.6636) = 2.3745e13*Dr**4.6636,
- ! RR in kg/(m**2*s)
- ! Dr (m) = 1.355e-3*RR**.2144 -> Dr (microns) = 1.355e3*RR**.2144
- !
- INDEXR=INT( 1.355E3*RR**.2144 + .5 )
- INDEXR=MAX( MDR3, MIN(INDEXR, MDRmax) )
- ELSE
- !
- !--- Assume fixed mean diameter of rain (0.45 mm) for high rain rates,
- ! instead vary N0r with rain rate
- !
- INDEXR=MDRmax
- ENDIF ! End IF (RR .LE. RR_DRmin) etc.
- VRAIN1=GAMMAR*VRAIN(INDEXR)
- ENDIF ! End IF (ARAIN .LE. 0.)
- INDEXR1=INDEXR ! For debugging only
- TOT_RAIN=THICK*QR+BLEND*ARAIN
- QTRAIN=TOT_RAIN/(THICK+BLDTRH*VRAIN1)
- PRLOSS=-TOT_RAIN/THICK
- RQR=RHO*QTRAIN
- !
- !--- RQR - time-averaged rain content (kg/m**3)
- !
- IF (RQR .LE. RQR_DRmin) THEN
- N0r=MAX(N0rmin, CN0r_DMRmin*RQR)
- INDEXR=MDRmin
- ELSE IF (RQR .GE. RQR_DRmax) THEN
- N0r=CN0r_DMRmax*RQR
- INDEXR=MDRmax
- ELSE
- N0r=N0r0
- INDEXR=MAX( XMRmin, MIN(CN0r0*RQR**.25, XMRmax) )
- ENDIF
- !
- IF (TC .LT. T_ICE) THEN
- PIACR=-PRLOSS
- ELSE
- !GFDL DWVr=WV-PCOND-QSW
- DWVr=WV-PCOND-QSWgrd !GFDL
- DUM=QW+PCOND
- IF (DWVr.LT.0. .AND. DUM.LE.EPSQ) THEN
- !
- !--- Rain evaporation
- !
- ! * RFACTOR - [GAMMAR**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5],
- ! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
- !
- ! * Units: RFACTOR - s**.5/m ; ABW - m**2/s ; VENTR - m**-2 ;
- ! N0r - m**-4 ; VENTR1 - m**2 ; VENTR2 - m**3/s**.5 ;
- ! CREVP - unitless
- !
- RFACTOR=SQRT(GAMMAR)*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2 !GFDL
- ABW=1./(RHO*XLV2/THERM_COND+1./DIFFUS)
- !
- !--- Note that VENTR1, VENTR2 lookup tables do not include the
- ! 1/Davg multiplier as in the ice tables
- !
- VENTR=N0r*(VENTR1(INDEXR)+RFACTOR*VENTR2(INDEXR))
- CREVP=ABW*VENTR*DTPH
- IF (CREVP .LT. Xratio) THEN
- DUM=DWVr*CREVP
- ELSE
- DUM=DWVr*(1.-EXP(-CREVP*DENOMW))/DENOMW
- ENDIF
- PREVP=MAX(DUM, PRLOSS)
- ELSE IF (QW .GT. EPSQ) THEN
- FWR=CRACW*GAMMAR*N0r*ACCRR(INDEXR)
- PRACW=MIN(1.,FWR)*QW
- ENDIF ! End IF (DWVr.LT.0. .AND. DUM.LE.EPSQ)
- !
- IF (TC.LT.0. .AND. TCC.LT.0.) THEN
- !
- !--- Biggs (1953) heteorogeneous freezing (e.g., Lin et al., 1983)
- ! - Rescaled mean drop diameter from microns (INDEXR) to mm (DUM) to prevent underflow
- !
- DUM=.001*FLOAT(INDEXR)
- DUM=(EXP(ABFR*TC)-1.)*DUM*DUM*DUM*DUM*DUM*DUM*DUM
- PIACR=MIN(CBFR*N0r*RRHO*DUM, QTRAIN)
- IF (QLICE .GT. EPSQ) THEN
- !
- !--- Freezing of rain by collisions w/ large ice
- !
- DUM=GAMMAR*VRAIN(INDEXR)
- DUM1=DUM-VSNOW
- !
- !--- DUM2 - Difference in spectral fall speeds of rain and
- ! large ice, parameterized following eq. (48) on p. 112 of
- ! Murakami (J. Meteor. Soc. Japan, 1990)
- !
- DUM2=SQRT(DUM1*DUM1+.04*DUM*VSNOW) !GFDL
- DUM1=5.E-12*INDEXR*INDEXR+2.E-12*INDEXR*INDEXS &
- & +.5E-12*INDEXS*INDEXS
- FIR=MIN(1., CIACR*NLICE*DUM1*DUM2)
- !
- !--- Future? Should COLLECTION BY SMALL ICE SHOULD BE INCLUDED???
- !
- PIACR=MIN(PIACR+FIR*QTRAIN, QTRAIN)
- ENDIF ! End IF (QLICE .GT. EPSQ)
- DUM=PREVP-PIACR
- If (DUM .LT. PRLOSS) THEN
- DUM1=PRLOSS/DUM
- PREVP=DUM1*PREVP
- PIACR=DUM1*PIACR
- ENDIF ! End If (DUM .LT. PRLOSS)
- ENDIF ! End IF (TC.LT.0. .AND. TCC.LT.0.)
- ENDIF ! End IF (TC .LT. T_ICE)
- ENDIF ! End IF (RAIN_logical)
- !
- !----------------------------------------------------------------------
- !---------------------- Main Budget Equations -------------------------
- !----------------------------------------------------------------------
- !
- !
- !-----------------------------------------------------------------------
- !--- Update fields, determine characteristics for next lower layer ----
- !-----------------------------------------------------------------------
- !
- !--- Carefully limit sinks of cloud water
- !
- DUM1=PIACW+PRAUT+PRACW-MIN(0.,PCOND)
- IF (DUM1 .GT. QW) THEN
- DUM=QW/DUM1
- PIACW=DUM*PIACW
- PIACWI=DUM*PIACWI
- PRAUT=DUM*PRAUT
- PRACW=DUM*PRACW
- IF (PCOND .LT. 0.) PCOND=DUM*PCOND
- ENDIF
- PIACWR=PIACW-PIACWI ! TC >= 0C
- !
- !--- QWnew - updated cloud water mixing ratio
- !
- DELW=PCOND-PIACW-PRAUT-PRACW
- QWnew=QW+DELW
- IF (QWnew .LE. EPSQ) QWnew=0.
- IF (QW.GT.0. .AND. QWnew.NE.0.) THEN
- DUM=QWnew/QW
- IF (DUM .LT. TOLER) QWnew=0.
- ENDIF
- !
- !--- Update temperature and water vapor mixing ratios
- !
- DELT= XLV1*(PCOND+PIEVP+PICND+PREVP) &
- & +XLS1*PIDEP+XLF1*(PIACWI+PIACR-PIMLT)
- Tnew=TK+DELT
- !
- DELV=-PCOND-PIDEP-PIEVP-PICND-PREVP
- WVnew=WV+DELV
- !
- !--- Update ice mixing ratios
- !
- !---
- ! * TOT_ICEnew - total mass (small & large) ice after microphysics,
- ! which is the sum of the total mass of large ice in the
- ! current layer and the flux of ice out of the grid box below
- ! * RimeF - Rime Factor, which is the mass ratio of total (unrimed &
- ! rimed) ice mass to the unrimed ice mass (>=1)
- ! * QInew - updated mixing ratio of total (large & small) ice in layer
- ! -> TOT_ICEnew=QInew*THICK+BLDTRH*QLICEnew*VSNOW
- ! -> But QLICEnew=QInew*FLIMASS, so
- ! -> TOT_ICEnew=QInew*(THICK+BLDTRH*FLIMASS*VSNOW)
- ! * ASNOWnew - updated accumulation of snow at bottom of grid cell
- !---
- !
- DELI=0.
- RimeF=1.
- IF (ICE_logical) THEN
- DELI=PIDEP+PIEVP+PIACWI+PIACR-PIMLT
- TOT_ICEnew=TOT_ICE+THICK*DELI
- IF (TOT_ICE.GT.0. .AND. TOT_ICEnew.NE.0.) THEN
- DUM=TOT_ICEnew/TOT_ICE
- IF (DUM .LT. TOLER) TOT_ICEnew=0.
- ENDIF
- IF (TOT_ICEnew .LE. CLIMIT) THEN
- TOT_ICEnew=0.
- RimeF=1.
- QInew=0.
- ASNOWnew=0.
- ELSE
- !
- !--- Update rime factor if appropriate
- !
- DUM=PIACWI+PIACR
- IF (DUM.LE.EPSQ .AND. PIDEP.LE.EPSQ) THEN
- RimeF=RimeF1
- ELSE
- !
- !--- Rime Factor, RimeF = (Total ice mass)/(Total unrimed ice mass)
- ! DUM1 - Total ice mass, rimed & unrimed
- ! DUM2 - Estimated mass of *unrimed* ice
- !
- DUM1=TOT_ICE+THICK*(PIDEP+DUM)
- DUM2=TOT_ICE/RimeF1+THICK*PIDEP
- IF (DUM2 .LE. 0.) THEN
- RimeF=RFmax
- ELSE
- RimeF=MIN(RFmax, MAX(1., DUM1/DUM2) )
- ENDIF
- ENDIF ! End IF (DUM.LE.EPSQ .AND. PIDEP.LE.EPSQ)
- QInew=TOT_ICEnew/(THICK+BLDTRH*FLIMASS*VSNOW)
- IF (QInew .LE. EPSQ) QInew=0.
- IF (QI.GT.0. .AND. QInew.NE.0.) THEN
- DUM=QInew/QI
- IF (DUM .LT. TOLER) QInew=0.
- ENDIF
- ASNOWnew=BLDTRH*FLIMASS*VSNOW*QInew
- IF (ASNOW.GT.0. .AND. ASNOWnew.NE.0.) THEN
- DUM=ASNOWnew/ASNOW
- IF (DUM .LT. TOLER) ASNOWnew=0.
- ENDIF
- ENDIF ! End IF (TOT_ICEnew .LE. CLIMIT)
- ENDIF ! End IF (ICE_logical)
- !
- !--- Update rain mixing ratios
- !
- !---
- ! * TOT_RAINnew - total mass of rain after microphysics
- ! current layer and the input flux of ice from above
- ! * VRAIN2 - time-averaged fall speed of rain in grid and below
- ! (with air resistance correction)
- ! * QRnew - updated rain mixing ratio in layer
- ! -> TOT_RAINnew=QRnew*(THICK+BLDTRH*VRAIN2)
- ! * ARAINnew - updated accumulation of rain at bottom of grid cell
- !---
- !
- DELR=PRAUT+PRACW+PIACWR-PIACR+PIMLT+PREVP+PICND
- TOT_RAINnew=TOT_RAIN+THICK*DELR
- IF (TOT_RAIN.GT.0. .AND. TOT_RAINnew.NE.0.) THEN
- DUM=TOT_RAINnew/TOT_RAIN
- IF (DUM .LT. TOLER) TOT_RAINnew=0.
- ENDIF
- IF (TOT_RAINnew .LE. CLIMIT) THEN
- TOT_RAINnew=0.
- VRAIN2=0.
- QRnew=0.
- ARAINnew=0.
- ELSE
- !
- !--- 1st guess time-averaged rain rate at bottom of grid box
- !
- RR=TOT_RAINnew/(DTPH*GAMMAR)
- !
- !--- Use same algorithm as above for calculating mean drop diameter
- ! (IDR, in microns), which is used to estimate the time-averaged
- ! fall speed of rain drops at the bottom of the grid layer. This
- ! isn't perfect, but the alternative is solving a transcendental
- ! equation that is numerically inefficient and nasty to program
- ! (coded in earlier versions of GSMCOLUMN prior to 8-22-01).
- !
- IF (RR .LE. RR_DRmin) THEN
- IDR=MDRmin
- ELSE IF (RR .LE. RR_DR1) THEN
- IDR=INT( 1.123E3*RR**.1947 + .5 )
- IDR=MAX( MDRmin, MIN(IDR, MDR1) )
- ELSE IF (RR .LE. RR_DR2) THEN
- IDR=INT( 1.225E3*RR**.2017 + .5 )
- IDR=MAX( MDR1, MIN(IDR, MDR2) )
- ELSE IF (RR .LE. RR_DR3) THEN
- IDR=INT( 1.3006E3*RR**.2083 + .5 )
- IDR=MAX( MDR2, MIN(IDR, MDR3) )
- ELSE IF (RR .LE. RR_DRmax) THEN
- IDR=INT( 1.355E3*RR**.2144 + .5 )
- IDR=MAX( MDR3, MIN(IDR, MDRmax) )
- ELSE
- IDR=MDRmax
- ENDIF ! End IF (RR .LE. RR_DRmin)
- VRAIN2=GAMMAR*VRAIN(IDR)
- QRnew=TOT_RAINnew/(THICK+BLDTRH*VRAIN2)
- IF (QRnew .LE. EPSQ) QRnew=0.
- IF (QR.GT.0. .AND. QRnew.NE.0.) THEN
- DUM=QRnew/QR
- IF (DUM .LT. TOLER) QRnew=0.
- ENDIF
- ARAINnew=BLDTRH*VRAIN2*QRnew
- IF (ARAIN.GT.0. .AND. ARAINnew.NE.0.) THEN
- DUM=ARAINnew/ARAIN
- IF (DUM .LT. TOLER) ARAINnew=0.
- ENDIF
- ENDIF
- !
- WCnew=QWnew+QRnew+QInew
- !
- !----------------------------------------------------------------------
- !-------------- Begin debugging & verification ------------------------
- !----------------------------------------------------------------------
- !
- !--- QT, QTnew - total water (vapor & condensate) before & after microphysics, resp.
- !
- QT=THICK*(WV+WC)+ARAIN+ASNOW
- QTnew=THICK*(WVnew+WCnew)+ARAINnew+ASNOWnew
- BUDGET=QT-QTnew
- !
- !--- Additional check on budget preservation, accounting for truncation effects
- !
- IF (PRINT_err) THEN
- DBG_logical=.FALSE.
- DUM=ABS(BUDGET)
- IF (DUM .GT. TOLER) THEN
- DUM=DUM/MIN(QT, QTnew)
- IF (DUM .GT. TOLER) DBG_logical=.TRUE.
- ENDIF
- !
- ! DUM=(RHgrd+.001)*QSInew
- ! IF ( (QWnew.GT.EPSQ) .OR. QRnew.GT.EPSQ .OR. WVnew.GT.DUM) &
- ! & .AND. TC.LT.T_ICE ) DBG_logical=.TRUE.
- !
- ! IF (TC.GT.5. .AND. QInew.GT.EPSQ) DBG_logical=.TRUE.
- !
- IF (WVnew.LT.EPSQ .OR. DBG_logical) THEN
- !
- WRITE(6,"(/2(a,i4),2(a,i2))") '{} i=',I_index,' j=',J_index, &
- & ' L=',L,' LSFC=',LSFC
- !
- ESW=1000.*FPVS0(Tnew)
- QSWnew=EPS*ESW/(PP-ESW)
- IF (TC.LT.0. .OR. Tnew .LT. 0.) THEN
- ESI=1000.*FPVS(Tnew)
- QSInew=EPS*ESI/(PP-ESI)
- ELSE
- QSI=QSW
- QSInew=QSWnew
- ENDIF
- WSnew=QSInew
- WRITE(6,"(4(a12,g11.4,1x))") &
- & '{} TCold=',TC,'TCnew=',Tnew-T0C,'P=',.01*PP,'RHO=',RHO, &
- & '{} THICK=',THICK,'RHold=',WV/WS,'RHnew=',WVnew/WSnew, &
- & 'RHgrd=',RHgrd, &
- & '{} RHWold=',WV/QSW,'RHWnew=',WVnew/QSWnew,'RHIold=',WV/QSI, &
- & 'RHInew=',WVnew/QSInew, &
- & '{} QSWold=',QSW,'QSWnew=',QSWnew,'QSIold=',QSI,'QSInew=',QSInew, &
- & '{} WSold=',WS,'WSnew=',WSnew,'WVold=',WV,'WVnew=',WVnew, &
- & '{} WCold=',WC,'WCnew=',WCnew,'QWold=',QW,'QWnew=',QWnew, &
- & '{} QIold=',QI,'QInew=',QInew,'QRold=',QR,'QRnew=',QRnew, &
- & '{} ARAINold=',ARAIN,'ARAINnew=',ARAINnew,'ASNOWold=',ASNOW, &
- & 'ASNOWnew=',ASNOWnew, &
- & '{} TOT_RAIN=',TOT_RAIN,'TOT_RAINnew=',TOT_RAINnew, &
- & 'TOT_ICE=',TOT_ICE,'TOT_ICEnew=',TOT_ICEnew, &
- & '{} BUDGET=',BUDGET,'QTold=',QT,'QTnew=',QTnew
- !
- WRITE(6,"(4(a12,g11.4,1x))") &
- & '{} DELT=',DELT,'DELV=',DELV,'DELW=',DELW,'DELI=',DELI, &
- & '{} DELR=',DELR,'PCOND=',PCOND,'PIDEP=',PIDEP,'PIEVP=',PIEVP, &
- & '{} PICND=',PICND,'PREVP=',PREVP,'PRAUT=',PRAUT,'PRACW=',PRACW, &
- & '{} PIACW=',PIACW,'PIACWI=',PIACWI,'PIACWR=',PIACWR,'PIMLT=', &
- & PIMLT, &
- & '{} PIACR=',PIACR
- !
- IF (ICE_logical) WRITE(6,"(4(a12,g11.4,1x))") &
- & '{} RimeF1=',RimeF1,'GAMMAS=',GAMMAS,'VrimeF=',VrimeF, &
- & 'VSNOW=',VSNOW, &
- & '{} INDEXS=',FLOAT(INDEXS),'FLARGE=',FLARGE,'FSMALL=',FSMALL, &
- & 'FLIMASS=',FLIMASS, &
- & '{} XSIMASS=',XSIMASS,'XLIMASS=',XLIMASS,'QLICE=',QLICE, &
- & 'QTICE=',QTICE, &
- & '{} NLICE=',NLICE,'NSmICE=',NSmICE,'PILOSS=',PILOSS, &
- & 'EMAIRI=',EMAIRI, &
- & '{} RimeF=',RimeF
- !
- IF (TOT_RAIN.GT.0. .OR. TOT_RAINnew.GT.0.) &
- & WRITE(6,"(4(a12,g11.4,1x))") &
- & '{} INDEXR1=',FLOAT(INDEXR1),'INDEXR=',FLOAT(INDEXR), &
- & 'GAMMAR=',GAMMAR,'N0r=',N0r, &
- & '{} VRAIN1=',VRAIN1,'VRAIN2=',VRAIN2,'QTRAIN=',QTRAIN,'RQR=',RQR, &
- & '{} PRLOSS=',PRLOSS,'VOLR1=',THICK+BLDTRH*VRAIN1, &
- & 'VOLR2=',THICK+BLDTRH*VRAIN2
- !
- IF (PRAUT .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} QW0=',QW0
- !
- IF (PRACW .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} FWR=',FWR
- !
- IF (PIACR .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} FIR=',FIR
- !
- DUM=PIMLT+PICND-PREVP-PIEVP
- IF (DUM.GT.0. .or. DWVi.NE.0.) &
- & WRITE(6,"(4(a12,g11.4,1x))") &
- & '{} TFACTOR=',TFACTOR,'DYNVIS=',DYNVIS, &
- & 'THERM_CON=',THERM_COND,'DIFFUS=',DIFFUS
- !
- IF (PREVP .LT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
- & '{} RFACTOR=',RFACTOR,'ABW=',ABW,'VENTR=',VENTR,'CREVP=',CREVP, &
- & '{} DWVr=',DWVr,'DENOMW=',DENOMW
- !
- IF (PIDEP.NE.0. .AND. DWVi.NE.0.) &
- & WRITE(6,"(4(a12,g11.4,1x))") &
- & '{} DWVi=',DWVi,'DENOMI=',DENOMI,'PIDEP_max=',PIDEP_max, &
- & 'SFACTOR=',SFACTOR, &
- & '{} ABI=',ABI,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
- & 'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
- & '{} VENTIS=',VENTIS,'DIDEP=',DIDEP
- !
- ! IF (PIDEP.GT.0. .AND. PCOND.NE.0.) &
- ! & WRITE(6,"(4(a12,g11.4,1x))") &
- ! & '{} DENOMW=',DENOMW,'DENOMWI=',DENOMWI,'DENOMF=',DENOMF, &
- ! & 'DUM2=',PCOND-PIACW
- !
- ! IF (FWS .GT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
- ! & '{} FWS=',FWS
- !
- DUM=PIMLT+PICND-PIEVP
- IF (DUM.GT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
- & '{} SFACTOR=',SFACTOR,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
- & 'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
- & '{} AIEVP=',AIEVP,'DIEVP=',DIEVP,'QSW0=',QSW0,'DWV0=',DWV0
- !
- ENDIF !-- IF (WVnew.LT.EPSQ .OR. DBG_logical) THEN
- ENDIF !-- IF (PRINT_err) THEN
- !
- !-----------------------------------------------------------------------
- !--------------- Water budget statistics & maximum values --------------
- !-----------------------------------------------------------------------
- !
- ! IF (PRINT_diag) THEN
- ! ITdx=MAX( ITLO, MIN( INT(Tnew-T0C), ITHI ) )
- ! IF (QInew .GT. EPSQ) NSTATS(ITdx,1)=NSTATS(ITdx,1)+1
- ! IF (QInew.GT.EPSQ .AND. QRnew+QWnew.GT.EPSQ) &
- ! & NSTATS(ITdx,2)=NSTATS(ITdx,2)+1
- ! IF (QWnew .GT. EPSQ) NSTATS(ITdx,3)=NSTATS(ITdx,3)+1
- ! IF (QRnew .GT. EPSQ) NSTATS(ITdx,4)=NSTATS(ITdx,4)+1
- ! !
- ! QMAX(ITdx,1)=MAX(QMAX(ITdx,1), QInew)
- ! QMAX(ITdx,2)=MAX(QMAX(ITdx,2), QWnew)
- ! QMAX(ITdx,3)=MAX(QMAX(ITdx,3), QRnew)
- ! QMAX(ITdx,4)=MAX(QMAX(ITdx,4), ASNOWnew)
- ! QMAX(ITdx,5)=MAX(QMAX(ITdx,5), ARAINnew)
- ! QTOT(ITdx,1)=QTOT(ITdx,1)+QInew*THICK
- ! QTOT(ITdx,2)=QTOT(ITdx,2)+QWnew*THICK
- ! QTOT(ITdx,3)=QTOT(ITdx,3)+QRnew*THICK
- ! !
- ! QTOT(ITdx,4)=QTOT(ITdx,4)+PCOND*THICK
- ! QTOT(ITdx,5)=QTOT(ITdx,5)+PICND*THICK
- ! QTOT(ITdx,6)=QTOT(ITdx,6)+PIEVP*THICK
- ! QTOT(ITdx,7)=QTOT(ITdx,7)+PIDEP*THICK
- ! QTOT(ITdx,8)=QTOT(ITdx,8)+PREVP*THICK
- ! QTOT(ITdx,9)=QTOT(ITdx,9)+PRAUT*THICK
- ! QTOT(ITdx,10)=QTOT(ITdx,10)+PRACW*THICK
- ! QTOT(ITdx,11)=QTOT(ITdx,11)+PIMLT*THICK
- ! QTOT(ITdx,12)=QTOT(ITdx,12)+PIACW*THICK
- ! QTOT(ITdx,13)=QTOT(ITdx,13)+PIACWI*THICK
- ! QTOT(ITdx,14)=QTOT(ITdx,14)+PIACWR*THICK
- ! QTOT(ITdx,15)=QTOT(ITdx,15)+PIACR*THICK
- ! !
- ! QTOT(ITdx,16)=QTOT(ITdx,16)+(WVnew-WV)*THICK
- ! QTOT(ITdx,17)=QTOT(ITdx,17)+(QWnew-QW)*THICK
- ! QTOT(ITdx,18)=QTOT(ITdx,18)+(QInew-QI)*THICK
- ! QTOT(ITdx,19)=QTOT(ITdx,19)+(QRnew-QR)*THICK
- ! QTOT(ITdx,20)=QTOT(ITdx,20)+(ARAINnew-ARAIN)
- ! QTOT(ITdx,21)=QTOT(ITdx,21)+(ASNOWnew-ASNOW)
- ! IF (QInew .GT. 0.) &
- ! & QTOT(ITdx,22)=QTOT(ITdx,22)+QInew*THICK/RimeF
- ! !
- ! ENDIF
- !
- !----------------------------------------------------------------------
- !------------------------- Update arrays ------------------------------
- !----------------------------------------------------------------------
- !
- T_col(L)=Tnew ! Updated temperature
- !
- QV_col(L)=max(EPSQ, WVnew/(1.+WVnew)) ! Updated specific humidity
- WC_col(L)=max(EPSQ, WCnew) ! Updated total condensate mixing ratio
- QI_col(L)=max(EPSQ, QInew) ! Updated ice mixing ratio
- QR_col(L)=max(EPSQ, QRnew) ! Updated rain mixing ratio
- QW_col(L)=max(EPSQ, QWnew) ! Updated cloud water mixing ratio
- RimeF_col(L)=RimeF ! Updated rime factor
- ASNOW=ASNOWnew ! Updated accumulated snow
- ARAIN=ARAINnew ! Updated accumulated rain
- !
- !#######################################################################
- !
- 10 CONTINUE ! ##### End "L" loop through model levels #####
- !
- !#######################################################################
- !
- !-----------------------------------------------------------------------
- !--------------------------- Return to GSMDRIVE -----------------------
- !-----------------------------------------------------------------------
- !
- CONTAINS
- !#######################################################################
- !--------- Produces accurate calculation of cloud condensation ---------
- !#######################################################################
- !
- REAL FUNCTION CONDENSE (PP, QW, TK, WV, RHgrd) !GFDL
- !
- !---------------------------------------------------------------------------------
- !------ The Asai (1965) algorithm takes into consideration the release of ------
- !------ latent heat in increasing the temperature & in increasing the ------
- !------ saturation mixing ratio (following the Clausius-Clapeyron eqn.). ------
- !---------------------------------------------------------------------------------
- !
- IMPLICIT NONE
- !
- INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
- REAL (KIND=HIGH_PRES), PARAMETER :: &
- & RHLIMIT=.001, RHLIMIT1=-RHLIMIT
- REAL (KIND=HIGH_PRES) :: COND, SSAT, WCdum
- !
- REAL,INTENT(IN) :: QW,PP,WV,TK,RHgrd !GFDL
- REAL WVdum,Tdum,XLV2,DWV,WS,ESW,XLV1,XLV
- integer nsteps
- !
- !-----------------------------------------------------------------------
- !
- !--- LV (T) is from Bolton (JAS, 1980)
- !
- XLV=3.148E6-2370.*TK
- XLV1=XLV*RCP
- XLV2=XLV*XLV*RCPRV
- Tdum=TK
- WVdum=WV
- WCdum=QW
- ESW=1000.*FPVS0(Tdum) ! Saturation vapor press w/r/t water
- WS=RHgrd*EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water
- DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
- SSAT=DWV/WS ! Supersaturation ratio
- CONDENSE=0.
- nsteps = 0
- DO WHILE ((SSAT.LT.RHLIMIT1 .AND. WCdum.GT.EPSQ) &
- & .OR. SSAT.GT.RHLIMIT)
- nsteps = nsteps + 1
- COND=DWV/(1.+XLV2*WS/(Tdum*Tdum)) ! Asai (1965, J. Japan)
- COND=MAX(COND, -WCdum) ! Limit cloud water evaporation
- Tdum=Tdum+XLV1*COND ! Updated temperature
- WVdum=WVdum-COND ! Updated water vapor mixing ratio
- WCdum=WCdum+COND ! Updated cloud water mixing ratio
- CONDENSE=CONDENSE+COND ! Total cloud water condensation
- ESW=1000.*FPVS0(Tdum) ! Updated saturation vapor press w/r/t water
- WS=RHgrd*EPS*ESW/(PP-ESW) ! Updated saturation mixing ratio w/r/t water
- DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
- SSAT=DWV/WS ! Grid-scale supersaturation ratio
- ENDDO
- !
- END FUNCTION CONDENSE
- !
- !#######################################################################
- !---------------- Calculate ice deposition at T<T_ICE ------------------
- !#######################################################################
- !
- REAL FUNCTION DEPOSIT (PP, Tdum, WVdum, RHgrd) !GFDL
- !
- !--- Also uses the Asai (1965) algorithm, but uses a different target
- ! vapor pressure for the adjustment
- !
- IMPLICIT NONE
- !
- INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
- REAL (KIND=HIGH_PRES), PARAMETER :: RHLIMIT=.001, &
- & RHLIMIT1=-RHLIMIT
- REAL (KIND=HIGH_PRES) :: DEP, SSAT
- !
- real,INTENT(IN) :: PP,RHgrd !GFDL
- real,INTENT(INOUT) :: WVdum,Tdum
- real ESI,WS,DWV
- !
- !-----------------------------------------------------------------------
- !
- ESI=1000.*FPVS(Tdum) ! Saturation vapor press w/r/t ice
- WS=RHgrd*EPS*ESI/(PP-ESI) ! Saturation mixing ratio
- DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
- SSAT=DWV/WS ! Supersaturation ratio
- DEPOSIT=0.
- DO WHILE (SSAT.GT.RHLIMIT .OR. SSAT.LT.RHLIMIT1)
- !
- !--- Note that XLVS2=LS*LV/(CP*RV)=LV*WS/(RV*T*T)*(LS/CP*DEP1),
- ! where WS is the saturation mixing ratio following Clausius-
- ! Clapeyron (see Asai,1965; Young,1993,p.405)
- !
- DEP=DWV/(1.+XLS2*WS/(Tdum*Tdum)) ! Asai (1965, J. Japan)
- Tdum=Tdum+XLS1*DEP ! Updated temperature
- WVdum=WVdum-DEP ! Updated ice mixing ratio
- DEPOSIT=DEPOSIT+DEP ! Total ice deposition
- ESI=1000.*FPVS(Tdum) ! Updated saturation vapor press w/r/t ice
- WS=RHgrd*EPS*ESI/(PP-ESI) ! Updated saturation mixing ratio w/r/t ice
- DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
- SSAT=DWV/WS ! Grid-scale supersaturation ratio
- ENDDO
- !
- END FUNCTION DEPOSIT
- !
- END SUBROUTINE EGCP01COLUMN
- !#######################################################################
- !------- Initialize constants & lookup tables for microphysics ---------
- !#######################################################################
- !
- ! SH 0211/2002
- !-----------------------------------------------------------------------
- SUBROUTINE etanewinit_HWRF (GSMDT,DT,DELX,DELY,LOWLYR,restart, &
- & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
- !HWRF & MP_RESTART_STATE,TBPVS_STATE,TBPVS0_STATE, &
- & ALLOWED_TO_READ, &
- & IDS,IDE,JDS,JDE,KDS,KDE, &
- & IMS,IME,JMS,JME,KMS,KME, &
- & ITS,ITE,JTS,JTE,KTS,KTE )
- !-----------------------------------------------------------------------
- !-------------------------------------------------------------------------------
- !--- SUBPROGRAM DOCUMENTATION BLOCK
- ! PRGRMMR: Ferrier ORG: W/NP22 DATE: February 2001
- ! Jin ORG: W/NP22 DATE: January 2002
- ! (Modification for WRF structure)
- !
- !-------------------------------------------------------------------------------
- ! ABSTRACT:
- ! * Reads various microphysical lookup tables used in COLUMN_MICRO
- ! * Lookup tables were created "offline" and are read in during execution
- ! * Creates lookup tables for saturation vapor pressure w/r/t water & ice
- !-------------------------------------------------------------------------------
- !
- ! USAGE: CALL etanewinit FROM SUBROUTINE GSMDRIVE AT MODEL START TIME
- !
- ! INPUT ARGUMENT LIST:
- ! DTPH - physics time step (s)
- !
- ! OUTPUT ARGUMENT LIST:
- ! NONE
- !
- ! OUTPUT FILES:
- ! NONE
- !
- ! SUBROUTINES:
- ! MY_GROWTH_RATES - lookup table for growth of nucleated ice
- ! GPVS - lookup tables for saturation vapor pressure (water, ice)
- !
- ! UNIQUE: NONE
- !
- ! LIBRARY: NONE
- !
- ! COMMON BLOCKS:
- ! CMICRO_CONS - constants used in GSMCOLUMN
- ! CMY600 - lookup table for growth of ice crystals in
- ! water saturated conditions (Miller & Young, 1979)
- ! IVENT_TABLES - lookup tables for ventilation effects of ice
- ! IACCR_TABLES - lookup tables for accretion rates of ice
- ! IMASS_TABLES - lookup tables for mass content of ice
- ! IRATE_TABLES - lookup tables for precipitation rates of ice
- ! IRIME_TABLES - lookup tables for increase in fall speed of rimed ice
- ! MAPOT - Need lat/lon grid resolution
- ! RVENT_TABLES - lookup tables for ventilation effects of rain
- ! RACCR_TABLES - lookup tables for accretion rates of rain
- ! RMASS_TABLES - lookup tables for mass content of rain
- ! RVELR_TABLES - lookup tables for fall speeds of rain
- ! RRATE_TABLES - lookup tables for precipitation rates of rain
- !
- ! ATTRIBUTES:
- ! LANGUAGE: FORTRAN 90
- ! MACHINE : IBM SP
- !
- !-----------------------------------------------------------------------
- !
- !
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- !-------------------------------------------------------------------------
- !-------------- Parameters & arrays for lookup tables --------------------
- !-------------------------------------------------------------------------
- !
- !--- Common block of constants used in column microphysics
- !
- !WRF
- ! real DLMD,DPHD
- !WRF
- !
- !-----------------------------------------------------------------------
- !--- Parameters & data statement for local calculations
- !-----------------------------------------------------------------------
- !
- INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3
- !
- ! VARIABLES PASSED IN
- integer,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
- & ,IMS,IME,JMS,JME,KMS,KME &
- & ,ITS,ITE,JTS,JTE,KTS,KTE
- !WRF
- INTEGER, DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: LOWLYR
- !
- real, INTENT(IN) :: DELX,DELY
- !HWRF real,DIMENSION(*), INTENT(INOUT) :: MP_RESTART_STATE
- !HWRF real,DIMENSION(NX), INTENT(INOUT) :: TBPVS_STATE,TBPVS0_STATE
- real,DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(OUT) :: &
- & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
- INTEGER, PARAMETER :: ITLO=-60, ITHI=40
- ! integer,DIMENSION(ITLO:ITHI,4),INTENT(INOUT) :: NSTATS
- ! real,DIMENSION(ITLO:ITHI,5),INTENT(INOUT) :: QMAX
- ! real,DIMENSION(ITLO:ITHI,22),INTENT(INOUT) :: QTOT
- ! real,INTENT(INOUT) :: PRECtot(2),PRECmax(2)
- real,INTENT(IN) :: DT,GSMDT
- LOGICAL,INTENT(IN) :: allowed_to_read,restart
- !
- !-----------------------------------------------------------------------
- ! LOCAL VARIABLES
- !-----------------------------------------------------------------------
- REAL :: BBFR,DTPH,PI,DX,Thour_print
- INTEGER :: I,IM,J,L,K,JTF,KTF,ITF
- INTEGER :: etampnew_unit1
- LOGICAL :: opened
- LOGICAL , EXTERNAL :: wrf_dm_on_monitor
- CHARACTER*80 errmess
- !
- !-----------------------------------------------------------------------
- !
- JTF=MIN0(JTE,JDE-1)
- KTF=MIN0(KTE,KDE-1)
- ITF=MIN0(ITE,IDE-1)
- !
- DO J=JTS,JTF
- DO I=ITS,ITF
- LOWLYR(I,J)=1
- ENDDO
- ENDDO
- !
- IF(.NOT.RESTART .AND. ALLOWED_TO_READ) THEN !HWRF
- CALL wrf_debug(1,'WARNING: F_ICE_PHY,F_RAIN_PHY AND F_RIMEF_PHY IS REINITIALIZED') !HWRF
- DO J = jts,jte
- DO K = kts,kte
- DO I= its,ite
- F_ICE_PHY(i,k,j)=0.
- F_RAIN_PHY(i,k,j)=0.
- F_RIMEF_PHY(i,k,j)=1.
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- !
- !-----------------------------------------------------------------------
- IF(ALLOWED_TO_READ)THEN
- !-----------------------------------------------------------------------
- !
- DX=SQRT((DELX)**2+(DELY)**2)/1000. ! Model resolution at equator (km) !GFDL
- DX=MIN(100., MAX(5., DX) )
- !
- !-- Relative humidity threshold for the onset of grid-scale condensation
- !!-- 9/1/01: Assume the following functional dependence for 5 - 100 km resolution:
- !! RHgrd=0.90 for dx=100 km, 0.98 for dx=5 km, where
- ! RHgrd=0.90+.08*((100.-DX)/95.)**.5
- !
- DTPH=MAX(GSMDT*60.,DT)
- DTPH=NINT(DTPH/DT)*DT
- !
- !--- Create lookup tables for saturation vapor pressure w/r/t water & ice
- !
- CALL GPVS
- !
- !--- Read in various lookup tables
- !
- IF ( wrf_dm_on_monitor() ) THEN
- DO i = 31,99
- INQUIRE ( i , OPENED = opened )
- IF ( .NOT. opened ) THEN
- etampnew_unit1 = i
- GOTO 2061
- ENDIF
- ENDDO
- etampnew_unit1 = -1
- 2061 CONTINUE
- ENDIF
- !
- CALL wrf_dm_bcast_bytes ( etampnew_unit1 , IWORDSIZE )
- !
- IF ( etampnew_unit1 < 0 ) THEN
- CALL wrf_error_fatal ( 'module_mp_hwrf: etanewinit: Can not find unused fortran unit to read in lookup table.' )
- ENDIF
- !
- IF ( wrf_dm_on_monitor() ) THEN
- !!was OPEN (UNIT=1,FILE="eta_micro_lookup.dat",FORM="UNFORMATTED")
- OPEN(UNIT=etampnew_unit1,FILE="ETAMPNEW_DATA", &
- & FORM="UNFORMATTED",STATUS="OLD",ERR=9061)
- !
- READ(etampnew_unit1) VENTR1
- READ(etampnew_unit1) VENTR2
- READ(etampnew_unit1) ACCRR
- READ(etampnew_unit1) MASSR
- READ(etampnew_unit1) VRAIN
- READ(etampnew_unit1) RRATE
- READ(etampnew_unit1) VENTI1
- READ(etampnew_unit1) VENTI2
- READ(etampnew_unit1) ACCRI
- READ(etampnew_unit1) MASSI
- READ(etampnew_unit1) VSNOWI
- READ(etampnew_unit1) VEL_RF
- ! read(etampnew_unit1) my_growth ! Applicable only for DTPH=180 s
- CLOSE (etampnew_unit1)
- ENDIF
- !
- CALL wrf_dm_bcast_bytes ( VENTR1 , size ( VENTR1 ) * RWORDSIZE )
- CALL wrf_dm_bcast_bytes ( VENTR2 , size ( VENTR2 ) * RWORDSIZE )
- CALL wrf_dm_bcast_bytes ( ACCRR , size ( ACCRR ) * RWORDSIZE )
- CALL wrf_dm_bcast_bytes ( MASSR , size ( MASSR ) * RWORDSIZE )
- CALL wrf_dm_bcast_bytes ( VRAIN , size ( VRAIN ) * RWORDSIZE )
- CALL wrf_dm_bcast_bytes ( RRATE , size ( RRATE ) * RWORDSIZE )
- CALL wrf_dm_bcast_bytes ( VENTI1 , size ( VENTI1 ) * RWORDSIZE )
- CALL wrf_dm_bcast_bytes ( VENTI2 , size ( VENTI2 ) * RWORDSIZE )
- CALL wrf_dm_bcast_bytes ( ACCRI , size ( ACCRI ) * RWORDSIZE )
- CALL wrf_dm_bcast_bytes ( MASSI , size ( MASSI ) * RWORDSIZE )
- CALL wrf_dm_bcast_bytes ( VSNOWI , size ( VSNOWI ) * RWORDSIZE )
- CALL wrf_dm_bcast_bytes ( VEL_RF , size ( VEL_RF ) * RWORDSIZE )
- !
- !--- Calculates coefficients for growth rates of ice nucleated in water
- ! saturated conditions, scaled by physics time step (lookup table)
- !
- CALL MY_GROWTH_RATES (DTPH)
- ! CALL MY_GROWTH_RATES (DTPH,MY_GROWTH)
- !
- PI=ACOS(-1.)
- !
- !--- Constants associated with Biggs (1953) freezing of rain, as parameterized
- ! following Lin et al. (JCAM, 1983) & Reisner et al. (1998, QJRMS).
- !
- ABFR=-0.66
- BBFR=100.
- CBFR=20.*PI*PI*BBFR*RHOL*1.E-21
- !
- !--- CIACW is used in calculating riming rates
- ! The assumed effective collection efficiency of cloud water rimed onto
- ! ice is =0.5 below:
- !
- CIACW=DTPH*0.25*PI*0.5*(1.E5)**C1
- !
- !--- CIACR is used in calculating freezing of rain colliding with large ice
- ! The assumed collection efficiency is 1.0
- !
- CIACR=PI*DTPH
- !
- !--- Based on rain lookup tables for mean diameters from 0.05 to 0.45 mm
- ! * Four different functional relationships of mean drop diameter as
- ! a function of rain rate (RR), derived based on simple fits to
- ! mass-weighted fall speeds of rain as functions of mean diameter
- ! from the lookup tables.
- !
- RR_DRmin=N0r0*RRATE(MDRmin) ! RR for mean drop diameter of .05 mm
- RR_DR1=N0r0*RRATE(MDR1) ! RR for mean drop diameter of .10 mm
- RR_DR2=N0r0*RRATE(MDR2) ! RR for mean drop diameter of .20 mm
- RR_DR3=N0r0*RRATE(MDR3) ! RR for mean drop diameter of .32 mm
- RR_DRmax=N0r0*RRATE(MDRmax) ! RR for mean drop diameter of .45 mm
- !
- RQR_DRmin=N0r0*MASSR(MDRmin) ! Rain content for mean drop diameter of .05 mm
- RQR_DR1=N0r0*MASSR(MDR1) ! Rain content for mean drop diameter of .10 mm
- RQR_DR2=N0r0*MASSR(MDR2) ! Rain content for mean drop diameter of .20 mm
- RQR_DR3=N0r0*MASSR(MDR3) ! Rain content for mean drop diameter of .32 mm
- RQR_DRmax=N0r0*MASSR(MDRmax) ! Rain content for mean drop diameter of .45 mm
- C_N0r0=PI*RHOL*N0r0
- CN0r0=1.E6/C_N0r0**.25
- CN0r_DMRmin=1./(PI*RHOL*DMRmin**4)
- CN0r_DMRmax=1./(PI*RHOL*DMRmax**4)
- !
- !--- CRACW is used in calculating collection of cloud water by rain (an
- ! assumed collection efficiency of 1.0)
- !
- CRACW=DTPH*0.25*PI*1.0
- !
- ESW0=1000.*FPVS0(T0C) ! Saturation vapor pressure at 0C
- RFmax=1.1**Nrime ! Maximum rime factor allowed
- !
- !------------------------------------------------------------------------
- !--------------- Constants passed through argument list -----------------
- !------------------------------------------------------------------------
- !
- !--- Important parameters for self collection (autoconversion) of
- ! cloud water to rain.
- !
- !--- CRAUT is proportional to the rate that cloud water is converted by
- ! self collection to rain (autoconversion rate)
- !
- CRAUT=1.-EXP(-1.E-3*DTPH)
- !
- !--- QAUT0 is the threshold cloud content for autoconversion to rain
- ! needed for droplets to reach a diameter of 20 microns (following
- ! Manton and Cotton, 1977; Banta and Hanson, 1987, JCAM)
- !--- QAUT0=1.2567, 0.8378, or 0.4189 g/m**3 for droplet number concentrations
- ! of 300, 200, and 100 cm**-3, respectively
- !
- QAUT0=PI*RHOL*NCW*(20.E-6)**3/6.
- !
- !--- For calculating snow optical depths by considering bulk density of
- ! snow based on emails from Q. Fu (6/27-28/01), where optical
- ! depth (T) = 1.5*SWP/(Reff*DENS), SWP is snow water path, Reff
- ! is effective radius, and DENS is the bulk density of snow.
- !
- ! SWP (kg/m**2)=(1.E-3 kg/g)*SWPrad, SWPrad in g/m**2 used in radiation
- ! T = 1.5*1.E3*SWPrad/(Reff*DENS)
- !
- ! See derivation for MASSI(INDEXS), note equal to RHO*QSNOW/NSNOW
- !
- ! SDENS=1.5e3/DENS, DENS=MASSI(INDEXS)/[PI*(1.E-6*INDEXS)**3]
- !
- DO I=MDImin,MDImax
- SDENS(I)=PI*1.5E-15*FLOAT(I*I*I)/MASSI(I)
- ENDDO
- !
- Thour_print=-DTPH/3600.
- ! SH 0211/2002
- ! IF (PRINT_diag) THEN
- !
- ! !-------- Total and maximum quantities
- ! !
- ! NSTATS=0 ! Microphysical statistics dealing w/ grid-point counts
- ! QMAX=0. ! Microphysical statistics dealing w/ hydrometeor mass
- ! QTOT=0. ! Microphysical statistics dealing w/ hydrometeor mass
- ! PRECmax=0. ! Maximum precip rates (rain, snow) at surface (mm/h)
- ! PRECtot=0. ! Total precipitation (rain, snow) accumulation at surface
- ! ENDIF
- !wrf
- !HWRF IF(.NOT.RESTART)THEN
- !HWRF MP_RESTART_STATE(MY_T1:MY_T2)=MY_GROWTH(MY_T1:MY_T2)
- !HWRF MP_RESTART_STATE(MY_T2+1)=C1XPVS0
- !HWRF MP_RESTART_STATE(MY_T2+2)=C2XPVS0
- !HWRF MP_RESTART_STATE(MY_T2+3)=C1XPVS
- !HWRF MP_RESTART_STATE(MY_T2+4)=C2XPVS
- !HWRF MP_RESTART_STATE(MY_T2+5)=CIACW
- !HWRF MP_RESTART_STATE(MY_T2+6)=CIACR
- !HWRF MP_RESTART_STATE(MY_T2+7)=CRACW
- !HWRF MP_RESTART_STATE(MY_T2+8)=CRAUT
- !HWRF TBPVS_STATE(1:NX) =TBPVS(1:NX)
- !HWRF TBPVS0_STATE(1:NX)=TBPVS0(1:NX)
- !HWRF ENDIF
- ENDIF ! Allowed_to_read
- RETURN
- !
- !-----------------------------------------------------------------------
- !
- 9061 CONTINUE
- WRITE( errmess , '(A,I4)' ) &
- 'module_mp_hwrf: error opening ETAMPNEW_DATA on unit ' &
- &, etampnew_unit1
- CALL wrf_error_fatal(errmess)
- !
- !-----------------------------------------------------------------------
- END SUBROUTINE etanewinit_HWRF
- !
- SUBROUTINE MY_GROWTH_RATES (DTPH)
- ! SUBROUTINE MY_GROWTH_RATES (DTPH,MY_GROWTH)
- !
- !--- Below are tabulated values for the predicted mass of ice crystals
- ! after 600 s of growth in water saturated conditions, based on
- ! calculations from Miller and Young (JAS, 1979). These values are
- ! crudely estimated from tabulated curves at 600 s from Fig. 6.9 of
- ! Young (1993). Values at temperatures colder than -27C were
- ! assumed to be invariant with temperature.
- !
- !--- Used to normalize Miller & Young (1979) calculations of ice growth
- ! over large time steps using their tabulated values at 600 s.
- ! Assumes 3D growth with time**1.5 following eq. (6.3) in Young (1993).
- !
- IMPLICIT NONE
- !
- REAL,INTENT(IN) :: DTPH
- !
- REAL DT_ICE
- REAL,DIMENSION(35) :: MY_600
- !WRF
- !
- !-----------------------------------------------------------------------
- DATA MY_600 / &
- & 5.5e-8, 1.4E-7, 2.8E-7, 6.E-7, 3.3E-6, &
- & 2.E-6, 9.E-7, 8.8E-7, 8.2E-7, 9.4e-7, &
- & 1.2E-6, 1.85E-6, 5.5E-6, 1.5E-5, 1.7E-5, &
- & 1.5E-5, 1.E-5, 3.4E-6, 1.85E-6, 1.35E-6, &
- & 1.05E-6, 1.E-6, 9.5E-7, 9.0E-7, 9.5E-7, &
- & 9.5E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7, &
- & 9.E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7 / ! -31 to -35 deg C
- !
- !-----------------------------------------------------------------------
- !
- DT_ICE=(DTPH/600.)**1.5
- MY_GROWTH=DT_ICE*MY_600
- !
- !-----------------------------------------------------------------------
- !
- END SUBROUTINE MY_GROWTH_RATES
- !
- !-----------------------------------------------------------------------
- !--------- Old GFS saturation vapor pressure lookup tables -----------
- !-----------------------------------------------------------------------
- !
- SUBROUTINE GPVS
- ! ******************************************************************
- !$$$ SUBPROGRAM DOCUMENTATION BLOCK
- ! . . .
- ! SUBPROGRAM: GPVS COMPUTE SATURATION VAPOR PRESSURE TABLE
- ! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
- !
- ! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE TABLE AS A FUNCTION OF
- ! TEMPERATURE FOR THE TABLE LOOKUP FUNCTION FPVS.
- ! EXACT SATURATION VAPOR PRESSURES ARE CALCULATED IN SUBPROGRAM FPVSX.
- ! THE CURRENT IMPLEMENTATION COMPUTES A TABLE WITH A LENGTH
- ! OF 7501 FOR TEMPERATURES RANGING FROM 180.0 TO 330.0 KELVIN.
- !
- ! PROGRAM HISTORY LOG:
- ! 91-05-07 IREDELL
- ! 94-12-30 IREDELL EXPAND TABLE
- ! 96-02-19 HONG ICE EFFECT
- ! 01-11-29 JIN MODIFIED FOR WRF
- !
- ! USAGE: CALL GPVS
- !
- ! SUBPROGRAMS CALLED:
- ! (FPVSX) - INLINABLE FUNCTION TO COMPUTE SATURATION VAPOR PRESSURE
- !
- ! COMMON BLOCKS:
- ! COMPVS - SCALING PARAMETERS AND TABLE FOR FUNCTION FPVS.
- !
- ! ATTRIBUTES:
- ! LANGUAGE: FORTRAN 90
- !
- !$$$
- IMPLICIT NONE
- real :: X,XINC,T
- integer :: JX
- !----------------------------------------------------------------------
- XINC=(XMAX-XMIN)/(NX-1)
- C1XPVS=1.-XMIN/XINC
- C2XPVS=1./XINC
- C1XPVS0=1.-XMIN/XINC
- C2XPVS0=1./XINC
- !
- DO JX=1,NX
- X=XMIN+(JX-1)*XINC
- T=X
- TBPVS(JX)=FPVSX(T)
- TBPVS0(JX)=FPVSX0(T)
- ENDDO
- !
- END SUBROUTINE GPVS
- !-----------------------------------------------------------------------
- !***********************************************************************
- !-----------------------------------------------------------------------
- REAL FUNCTION FPVS(T)
- !-----------------------------------------------------------------------
- !$$$ SUBPROGRAM DOCUMENTATION BLOCK
- ! . . .
- ! SUBPROGRAM: FPVS COMPUTE SATURATION VAPOR PRESSURE
- ! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
- !
- ! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE FROM THE TEMPERATURE.
- ! A LINEAR INTERPOLATION IS DONE BETWEEN VALUES IN A LOOKUP TABLE
- ! COMPUTED IN GPVS. SEE DOCUMENTATION FOR FPVSX FOR DETAILS.
- ! INPUT VALUES OUTSIDE TABLE RANGE ARE RESET TO TABLE EXTREMA.
- ! THE INTERPOLATION ACCURACY IS ALMOST 6 DECIMAL PLACES.
- ! ON THE CRAY, FPVS IS ABOUT 4 TIMES FASTER THAN EXACT CALCULATION.
- ! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
- !
- ! PROGRAM HISTORY LOG:
- ! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
- ! 94-12-30 IREDELL EXPAND TABLE
- ! 96-02-19 HONG ICE EFFECT
- ! 01-11-29 JIN MODIFIED FOR WRF
- !
- ! USAGE: PVS=FPVS(T)
- !
- ! INPUT ARGUMENT LIST:
- ! T - REAL TEMPERATURE IN KELVIN
- !
- ! OUTPUT ARGUMENT LIST:
- ! FPVS - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
- !
- ! ATTRIBUTES:
- ! LANGUAGE: FORTRAN 90
- !$$$
- IMPLICIT NONE
- real,INTENT(IN) :: T
- real XJ
- integer :: JX
- !-----------------------------------------------------------------------
- XJ=MIN(MAX(C1XPVS+C2XPVS*T,1.),FLOAT(NX))
- JX=MIN(XJ,NX-1.)
- FPVS=TBPVS(JX)+(XJ-JX)*(TBPVS(JX+1)-TBPVS(JX))
- !
- END FUNCTION FPVS
- !-----------------------------------------------------------------------
- !-----------------------------------------------------------------------
- REAL FUNCTION FPVS0(T)
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- real,INTENT(IN) :: T
- real :: XJ1
- integer :: JX1
- !-----------------------------------------------------------------------
- XJ1=MIN(MAX(C1XPVS0+C2XPVS0*T,1.),FLOAT(NX))
- JX1=MIN(XJ1,NX-1.)
- FPVS0=TBPVS0(JX1)+(XJ1-JX1)*(TBPVS0(JX1+1)-TBPVS0(JX1))
- !
- END FUNCTION FPVS0
- !-----------------------------------------------------------------------
- !***********************************************************************
- !-----------------------------------------------------------------------
- REAL FUNCTION FPVSX(T)
- !-----------------------------------------------------------------------
- !$$$ SUBPROGRAM DOCUMENTATION BLOCK
- ! . . .
- ! SUBPROGRAM: FPVSX COMPUTE SATURATION VAPOR PRESSURE
- ! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
- !
- ! ABSTRACT: EXACTLY COMPUTE SATURATION VAPOR PRESSURE FROM TEMPERATURE.
- ! THE WATER MODEL ASSUMES A PERFECT GAS, CONSTANT SPECIFIC HEATS
- ! FOR GAS AND LIQUID, AND NEGLECTS THE VOLUME OF THE LIQUID.
- ! THE MODEL DOES ACCOUNT FOR THE VARIATION OF THE LATENT HEAT
- ! OF CONDENSATION WITH TEMPERATURE. THE ICE OPTION IS NOT INCLUDED.
- ! THE CLAUSIUS-CLAPEYRON EQUATION IS INTEGRATED FROM THE TRIPLE POINT
- ! TO GET THE FORMULA
- ! PVS=PSATK*(TR**XA)*EXP(XB*(1.-TR))
- ! WHERE TR IS TTP/T AND OTHER VALUES ARE PHYSICAL CONSTANTS
- ! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
- !
- ! PROGRAM HISTORY LOG:
- ! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
- ! 94-12-30 IREDELL EXACT COMPUTATION
- ! 96-02-19 HONG ICE EFFECT
- ! 01-11-29 JIN MODIFIED FOR WRF
- !
- ! USAGE: PVS=FPVSX(T)
- ! REFERENCE: EMANUEL(1994),116-117
- !
- ! INPUT ARGUMENT LIST:
- ! T - REAL TEMPERATURE IN KELVIN
- !
- ! OUTPUT ARGUMENT LIST:
- ! FPVSX - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
- !
- ! ATTRIBUTES:
- ! LANGUAGE: FORTRAN 90
- !$$$
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- real, parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 &
- , CLIQ=4.1855E+3,CVAP= 1.8460E+3 &
- , CICE=2.1060E+3,HSUB=2.8340E+6
- !
- real, parameter :: PSATK=PSAT*1.E-3
- real, parameter :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP)
- real, parameter :: DLDTI=CVAP-CICE &
- , XAI=-DLDTI/RV,XBI=XAI+HSUB/(RV*TTP)
- real T,TR
- !-----------------------------------------------------------------------
- TR=TTP/T
- !
- IF(T.GE.TTP)THEN
- FPVSX=PSATK*(TR**XA)*EXP(XB*(1.-TR))
- ELSE
- FPVSX=PSATK*(TR**XAI)*EXP(XBI*(1.-TR))
- ENDIF
- !
- END FUNCTION FPVSX
- !-----------------------------------------------------------------------
- !-----------------------------------------------------------------------
- REAL FUNCTION FPVSX0(T)
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- real,parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 &
- , CLIQ=4.1855E+3,CVAP=1.8460E+3 &
- , CICE=2.1060E+3,HSUB=2.8340E+6
- real,PARAMETER :: PSATK=PSAT*1.E-3
- real,PARAMETER :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP)
- real,PARAMETER :: DLDTI=CVAP-CICE &
- , XAI=-DLDT/RV,XBI=XA+HSUB/(RV*TTP)
- real :: T,TR
- !-----------------------------------------------------------------------
- TR=TTP/T
- FPVSX0=PSATK*(TR**XA)*EXP(XB*(1.-TR))
- !
- END FUNCTION FPVSX0
- !
- END MODULE module_mp_HWRF