!>\file module_MP_FER_HIRES.F90 !! "Modified" fer_hires microphysics - 11 July 2016 version !! !! (1) Ice nucleation: Fletcher (1962) replaces Meyers et al. (1992) !! (2) Cloud ice is a simple function of the number concentration from (1), and it !! is no longer a fractional function of the large ice. Thus, the FLARGE & ! FSMALL parameters are no longer used. ! (3) T_ICE_init=-12 deg C provides a slight delay in the initial onset of ice. ! (4) NLImax is a function of rime factor (RF) and temperature. ! a) For RF>10, NLImax=1.e3. Mean ice diameters can exceed the 1 mm maximum ! size in the tables so that NLICE=NLImax=1.e3. ! b) Otherwise, NLImax is 10 L-1 at 0C and decreasing to 5 L-1 at <=-40C. ! NLICE>NLImax at the maximum ice diameter of 1 mm. ! (5) Can turn off ice processes by setting T_ICE & T_ICE_init to be < -100 deg C ! (6) Modified the homogeneous freezing of cloud water when TNLImax. ! (10) Ice deposition does not change the rime factor (RF) when RF>=10 & T>T_ICE. ! (11) Limit GAMMAS to <=1.5 (air resistance impact on ice fall speeds) ! (12) NSImax is maximum # conc of ice crsytals. At cold temperature NSImax is ! calculated based on assuming 10% of total ice content is due to cloud ice. ! !-- Further modifications starting on 23 July 2015 ! (13) RHgrd is passed in as an input argument so that it can vary for different ! domains (RHgrd=0.98 for 12-km parent, 1.0 for 3-km nests) ! (14) Use the old "PRAUT" cloud water autoconversion *threshold* (QAUT0) !-- Further modifications starting on 28 July 2015 ! (15) Added calculations for radar reflectivity and number concentrations of ! rain (Nrain) and precipitating ice (Nsnow). ! (16) Removed double counting of air resistance term for riming onto ice (PIACW) ! (17) The maximum rime factor (RFmx) is now a function of MASSI(INDEXS), accounting ! for the increase in unrimed ice particle densities as values of INDEXS ! decrease from the maximum upper limit of 1000 microns to the lower limit of ! 50 microns, coinciding with the assumed size of cloud ice; see lines 1128-1134. ! (18) A new closure is used for updating the rime factor, which is described in ! detail near lines 1643-1682. The revised code is near lines 1683-1718. ! (19) Restructured the two-pass algorithm to be more robust, removed the HAIL ! & LARGE_RF logical variables so that NLICE>NLImax can occur. ! (20) Increased nsimax (see !aug27 below) ! (21) Modified the rain sedimentation (see two !aug27 blocks below) ! (22) NInuclei is the lower of Fletcher (1962), Cooper (1986), or NSImax. ! (23) NLImax is no longer used or enforced. Instead, INDEXS=MDImax when RF>20, ! else INDEXS is a function of temperature. Look for !sep10 comment. ! (24) An override was inserted for (18), such that the rime density is not diluted ! diluted when RF>20. Look for !sep10 comment. ! (25) Radar reflectivity calculations were changes to reduce radar bright bands, ! limit enhanced, mixed-phase reflectivity to RF>=20. Look for !sep10 comments. ! (26) NLICE is not to exceed NSI_max (250 L^-1) when RF<20. Look for !sep16 comments. ! Commented out! (28) Increase hail fall speeds using Thompson et al. (2008). Look for !sep22 comments. ! (29) Modify NLImax, INDEXS for RF>=20. Look for !sep22 comments. ! (30) Check on NSmICE, Vci based on whether FLIMASS<1. Look for !sep22a comments. ! Revised in (34)! (31) Introduced RFlag logical, which if =T enforces a lower limit of drop sizes not ! to go below INDEXRmin and N0r is adjusted. Look for !nov25 comments (corrections, ! refinements to sep25 & nov18 versions, includes an additional fix in nov25-fix). ! Also set INDEXRmin=500 rather than 250 microns. !----------------------------------------------------------------------------- !--- The following changes now refer to dates when those were made in 2016. !----------------------------------------------------------------------------- ! (32) Convective (RF>=20, Ng~10 L^-1, RHOg~500 kg m^-3), transition (RF=10, Ng~25 L^-1, ! RHOg~300 kg m^-3), & stratiform (RF<2) profiles are blended based on RF. !mar08 ! (33) Fixed bug in Biggs' freezing, put back in collisional drop freezing. !mar03 ! (34) Changes in (31) are revised so that INDEXRmin at and below 0C level is ! based on a rain rate equal to the snowfall rate above the 0C level. !mar03 ! (35) Increase radar reflectivity when RF>10 and RQSnew > 2.5 g m^-3. !mar12 ! (36) !mar10 combines all elements of (32)-(35) together. ! (37) Bug fixes for the changes in (34) and the RFLAG variable !apr18 ! (38) Revised Schumann-Ludlam limit. !apr18 ! (39) Simplified PCOND (cloud cond/evap) calculation !apr21 ! (40) Slight change in calculating RF. !apr22 ! (41) Reduce RF values for calculating mean sizes of snow, graupel, sleet/hail !apr22a ! (42) Increase reflectivity from large, wet, high rime factor ice (graupel) by ! assuming |Kw|**2/|Ki|**2 = 0.224 (Smith, 1984, JCAM). ! (43) Major restructuring of code to allow N0r to vary from N0r0 !may11 ! (44) More major restructuring of code to use fixed XLS, XLV, XLF !may12 ! (45) Increased VEL_INC ~ VrimeF**2, put the enhanced graupel/hail fall speeds ! from Thompson into the code but only in limited circumstances, restructured ! and streamlined the INDEXS calculation, removed the upper limit for ! for the vapor mixing ratio is at water saturation when calculating ice ! deposition, and N0r is gradually increased for conditions supporting ! drizzle when rain contents decrease below 0.25 g/m**3. !may17 ! (46) The may11 code changes that increase N0r0 when rain contents exceed 1 g m^-3 ! have been removed, limit the number of iterations calculating final rain ! parameters, remove the revised N0r calculation for reflectivity. All of ! the changes following those made in the may10 code. !may20 ! (47) Reduce the assumed # concentration of hail/sleet when RF>10 from 5 L^-1 to ! 1 L^-1, and also reduce it for graupel when RF>5 from 10 L^-1 to 5 L^-1. ! This is being done to try and make greater use of the Thompson graupel/hail ! fallspeeds by having INDEXS==MDImax. ! (48) Increased NCW from 200e6 to 300e6 for a more delayed onset of drizzle, ! simplified drizzle algorithm to reduce/eliminate N0r bulls eyes and to allow ! for supercooled drizzle, and set limits for 8.e6 <= N0r (m^-4) <= 1.e9 !may31 ! (49) Further restructuring of code to better define STRAT, DRZL logicals, ! add these rain flags to mprates arrays !jun01 ! (50) Increase in reflectivity due to wet ice was commented out. ! (51) Fixed minor bug to update INDEXR2 in the "rain_pass: do" loop. !jun13 ! (52) Final changes to Nsnow for boosting reflectivities from ice for ! mass contents exceeding 5 g m^-3. !jun16 ! (53) Cosmetic changes only that do not affect the calculations. Removed old, unused ! diagnostic arrays. Updated comments. ! !----------------------------------------------------------------------------- ! MODULE MODULE_MP_FER_HIRES ! !----------------------------------------------------------------------------- #ifdef MPI USE mpi #endif USE machine !MZ !MZ USE MODULE_CONSTANTS,ONLY : PI, CP, EPSQ, GRAV=>G, RHOL=>RHOWATER, & !MZ RD=>R_D, RV=>R_V, T0C=>TIW, EPS=>EP_2, EPS1=>EP_1, CLIQ, CICE, & !MZ XLV !MZ !MZ temporary values copied from module_CONSTANTS; ideally they come from host model !side REAL, PARAMETER :: pi=3.141592653589793 ! ludolf number REAL, PARAMETER :: cp=1004.6 ! spec. heat for dry air at constant pressure REAL, PARAMETER :: epsq=1.e-12 ! floor value for specific humidity (kg/kg) REAL, PARAMETER :: grav= 9.8060226 ! gravity REAL, PARAMETER :: RHOL=1000. ! density of water (kg/m3) REAL, PARAMETER :: RD=287.04 ! gas constant for dry air REAL, PARAMETER :: RV=461.6 ! gas constant for water vapor REAL, PARAMETER :: T0C= 273.15 ! melting point REAL, PARAMETER :: EPS=RD/RV REAL, PARAMETER :: EPS1=RV/RD-1. REAL, PARAMETER :: CLIQ = 4190. ! MZ: inconsistent value below REAL, PARAMETER :: CICE = 2106. REAL, PARAMETER :: XLV = 2.5E6 !----------------------------------------------------------------------------- PUBLIC :: FERRIER_INIT_HR, GPVS_HR,FPVS,FPVS0,NX !----------------------------------------------------------------------------- REAL,PRIVATE,SAVE :: ABFR, CBFR, CIACW, CIACR, C_N0r0, C_NR, Crain, & !jul28 & CRACW, ARAUT, BRAUT, ESW0, RFmx1, ARcw, RH_NgC, RH_NgT, & !jul31 !mar08 & RR_DRmin, RR_DR1, RR_DR2, RR_DR3, RR_DR4, RR_DR5, RR_DRmax, & !may17 & BETA6, & & RQhail, AVhail, BVhail, QAUT0 !may17 ! INTEGER,PRIVATE,PARAMETER :: INDEXRstrmax=500 !mar03, stratiform maximum REAL,PUBLIC,SAVE :: CN0r0, CN0r_DMRmin, CN0r_DMRmax, & RFmax, RQR_DRmax, RQR_DRmin ! INTEGER, PRIVATE,PARAMETER :: MY_T1=1, MY_T2=35 REAL,PRIVATE,DIMENSION(MY_T1:MY_T2),SAVE :: MY_GROWTH_NMM ! REAL, PRIVATE,PARAMETER :: DMImin=.05e-3, DMImax=1.e-3, & & DelDMI=1.e-6,XMImin=1.e6*DMImin REAL, PUBLIC,PARAMETER :: XMImax=1.e6*DMImax, XMIexp=.0536 INTEGER, PUBLIC,PARAMETER :: MDImin=XMImin, MDImax=XMImax REAL, ALLOCATABLE, DIMENSION(:) :: & & ACCRI,VSNOWI,VENTI1,VENTI2 REAL, PUBLIC,DIMENSION(MDImin:MDImax) :: SDENS !-- For RRTM ! REAL, PRIVATE,PARAMETER :: DMRmin=.05e-3, DMRmax=1.0e-3, & & DelDMR=1.e-6, XMRmin=1.e6*DMRmin, XMRmax=1.e6*DMRmax INTEGER, PUBLIC,PARAMETER :: MDRmin=XMRmin, MDRmax=XMRmax ! REAL, ALLOCATABLE, DIMENSION(:):: & & ACCRR,MASSR,RRATE,VRAIN,VENTR1,VENTR2 ! INTEGER, PRIVATE,PARAMETER :: Nrime=40 REAL, ALLOCATABLE, DIMENSION(:,:) :: VEL_RF ! INTEGER,PARAMETER :: NX=7501 REAL, PARAMETER :: XMIN=180.0,XMAX=330.0 REAL, DIMENSION(NX),PUBLIC,SAVE :: TBPVS,TBPVS0 REAL, PUBLIC,SAVE :: C1XPVS0,C2XPVS0,C1XPVS,C2XPVS ! REAL,DIMENSION(MY_T2+8) :: MP_RESTART_STATE REAL,DIMENSION(nx) :: TBPVS_STATE,TBPVS0_STATE ! REAL, PRIVATE,PARAMETER :: CVAP=1846., XLF=3.3358e+5, XLS=XLV+XLF & & ,EPSQ1=1.001*EPSQ, RCP=1./CP, RCPRV=RCP/RV, RGRAV=1./GRAV & & ,RRHOL=1./RHOL, XLV1=XLV/CP, XLF1=XLF/CP, XLS1=XLS/CP & & ,XLV2=XLV*XLV/RV, XLS2=XLS*XLS/RV & & ,XLV3=XLV*XLV*RCPRV, XLS3=XLS*XLS*RCPRV & !--- 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, DMR4=0.45E-3 & & ,DMR5=0.67E-3 & & ,XMR1=1.e6*DMR1, XMR2=1.e6*DMR2, XMR3=1.e6*DMR3 & & ,XMR4=1.e6*DMR4, XMR5=1.e6*DMR5, RQRmix=0.05E-3, RQSmix=1.E-3 & !jul28 !apr27 & ,Cdry=1.634e13, Cwet=1./.224 !jul28 !apr27 INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3, MDR4=XMR4 & & , MDR5=XMR5 !-- Debug 20120111 LOGICAL, SAVE :: WARN1=.TRUE.,WARN2=.TRUE.,WARN3=.TRUE.,WARN5=.TRUE. REAL, SAVE :: Pwarn=75.E2, QTwarn=1.E-3 INTEGER, PARAMETER :: MAX_ITERATIONS=10 ! ! ====================================================================== !--- Important tunable parameters that are exported to other modules ! * 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 ! !-- To turn off ice processes, set T_ICE & T_ICE_init to <= -100. (i.e., -100 C) ! ! * NSImax - maximum number concentrations (m**-3) of small ice crystals ! * 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 1.0 mm ! * N0rmin - minimum intercept (m**-4) for rain drops ! * NCW - number concentrations of cloud droplets (m**-3) ! ====================================================================== REAL, PUBLIC,PARAMETER :: & & RHgrd_in=1. & &, P_RHgrd_out=850.E2 & & ,T_ICE=-40. & & ,T_ICEK=T0C+T_ICE & & ,T_ICE_init=-12. & & ,NSI_max=250.E3 & & ,NLImin=1.0E3 & & ,N0r0=8.E6 & & ,N0rmin=1.E4 & !! based on Aligo's email,NCW is changed to 250E6 & ,NCW=250.E6 !HWRF & ,NCW=300.E6 !- 100.e6 (maritime), 500.e6 (continental) !--- Other public variables passed to other routines: REAL, ALLOCATABLE ,DIMENSION(:) :: MASSI ! CONTAINS !----------------------------------------------------------------------- !----------------------------------------------------------------------- !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& !----------------------------------------------------------------------- !>\ingroup hafs_famp !! This is the driver scheme of Ferrier-Aligo microphysics scheme. !! NOTE: The only differences between FER_HIRES and FER_HIRES_ADVECT !! is that the QT, and F_* are all local variables in the advected !! version, and QRIMEF is only in the advected version. The innards !! are all the same. SUBROUTINE FER_HIRES (DT,RHgrd, & & prsi,p_phy,t_phy, & & q,qt, & & LOWLYR,SR,TRAIN_PHY, & & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, & & QC,QR,QS, & & RAINNC,RAINNCV, & & threads, & & ims,ime, lm, & & d_ss, & & refl_10cm,DX1 ) !----------------------------------------------------------------------- IMPLICIT NONE !----------------------------------------------------------------------- INTEGER,INTENT(IN) :: D_SS,IMS,IME,LM,DX1 REAL, INTENT(IN) :: DT,RHgrd INTEGER, INTENT(IN) :: THREADS REAL, INTENT(IN), DIMENSION(ims:ime, lm+1):: & & prsi REAL, INTENT(IN), DIMENSION(ims:ime, lm):: & & p_phy REAL, INTENT(INOUT), DIMENSION(ims:ime, lm):: & & q,qt,t_phy REAL, INTENT(INOUT), DIMENSION(ims:ime, lm ):: & !Aligo Oct 23,2019: dry mixing ratio for cloud species & qc,qr,qs REAL, INTENT(INOUT), DIMENSION(ims:ime, lm) :: & & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY REAL, INTENT(OUT), DIMENSION(ims:ime, lm) :: & & refl_10cm REAL, INTENT(INOUT), DIMENSION(ims:ime) :: & & RAINNC,RAINNCV REAL, INTENT(OUT), DIMENSION(ims:ime):: SR REAL, INTENT(OUT), DIMENSION( ims:ime, lm ) :: & & TRAIN_PHY ! INTEGER, DIMENSION( ims:ime ),INTENT(INOUT) :: LOWLYR !----------------------------------------------------------------------- ! LOCAL VARS !----------------------------------------------------------------------- ! 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):: APREC,PREC,ACPREC INTEGER :: I,K,KK REAL :: wc, RDIS, BETA6 !------------------------------------------------------------------------ ! For subroutine EGCP01COLUMN_hr !----------------------------------------------------------------------- INTEGER :: LSFC,I_index,J_index,L INTEGER,DIMENSION(ims:ime) :: LMH REAL :: TC,QI,QRdum,QW,Fice,Frain,DUM,ASNOW,ARAIN REAL,DIMENSION(lm) :: P_col,Q_col,T_col,WC_col, & RimeF_col,QI_col,QR_col,QW_col, THICK_col,DPCOL,pcond1d, & pidep1d,piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d, & pimlt1d,praut1d,pracw1d,prevp1d,pisub1d,pevap1d,DBZ_col, & NR_col,NS_col,vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d, & INDEXS1d,INDEXR1d,RFlag1d,RHC_col ! !----------------------------------------------------------------------- !********************************************************************** !----------------------------------------------------------------------- ! ! MZ: HWRF start !---------- !2015-03-30, recalculate some constants which may depend on phy time step CALL MY_GROWTH_RATES_NMM_hr (DT) !--- CIACW is used in calculating riming rates ! The assumed effective collection efficiency of cloud water rimed onto ! ice is =0.5 below: ! CIACW=DT*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*DT ! !--- CRACW is used in calculating collection of cloud water by rain (an ! assumed collection efficiency of 1.0) ! CRACW=DT*0.25*PI*1.0 ! !-- See comments in subroutine etanewhr_init starting with variable RDIS= ! !-- Relative dispersion == standard deviation of droplet spectrum / mean radius ! (see pp 1542-1543, Liu & Daum, JAS, 2004) RDIS=0.5 !-- relative dispersion of droplet spectrum BETA6=( (1.+3.*RDIS*RDIS)*(1.+4.*RDIS*RDIS)*(1.+5.*RDIS*RDIS)/ & & ((1.+RDIS*RDIS)*(1.+2.*RDIS*RDIS) ) ) BRAUT=DT*1.1E10*BETA6/NCW !! END OF adding, 2015-03-30 !----------- ! MZ: HWRF end ! DO i = ims,ime ACPREC(i)=0. APREC (i)=0. PREC (i)=0. SR (i)=0. ENDDO DO k = 1,lm DO i = ims,ime TRAIN_PHY (i,k)=0. ENDDO ENDDO !----------------------------------------------------------------------- !-- Start of original driver for EGCP01COLUMN_hr !----------------------------------------------------------------------- ! DO I=IMS,IME LSFC=LM-LOWLYR(I)+1 ! "L" of surface DO K=1,LM DPCOL(K)=prsi(I,K)-prsi(I,K+1) ENDDO ! !--- Initialize column data (1D arrays) ! L=LM !-- qt = CWM, total condensate IF (qt(I,L) .LE. EPSQ) qt(I,L)=EPSQ F_ice_phy(I,L)=1. F_rain_phy(I,L)=0. F_RimeF_phy(I,L)=1. do L=LM,1,-1 P_col(L)=P_phy(I,L) ! !--- Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc) ! THICK_col(L)=DPCOL(L)*RGRAV T_col(L)=T_phy(I,L) TC=T_col(L)-T0C Q_col(L)=max(EPSQ, q(I,L)) IF (qt(I,L) .LE. EPSQ1) THEN WC_col(L)=0. IF (TC .LT. T_ICE) THEN F_ice_phy(I,L)=1. ELSE F_ice_phy(I,L)=0. ENDIF F_rain_phy(I,L)=0. F_RimeF_phy(I,L)=1. ELSE WC_col(L)=qt(I,L) !-- Debug 20120111 ! TC==TC will fail if NaN, preventing unnecessary error messages IF (WC_col(L)>QTwarn .AND. P_col(L)1 g/kg condensate in stratosphere; I,L,TC,P,QT=', & I,L,TC,.01*P_col(L),1000.*WC_col(L) QTwarn=MAX(WC_col(L),10.*QTwarn) Pwarn=MIN(P_col(L),0.5*Pwarn) ENDIF !-- TC/=TC will pass if TC is NaN IF (WARN5 .AND. TC/=TC) THEN WRITE(0,*) 'WARN5: NaN temperature; I,L,P=',I,L,.01*P_col(L) WARN5=.FALSE. ENDIF ENDIF IF (T_ICE<=-100.) F_ice_phy(I,L)=0. ! ! ! !--- Determine composition of condensate in terms of ! ! cloud water, ice, & rain ! ! WC=WC_col(L) QI=0. QRdum=0. QW=0. Fice=F_ice_phy(I,L) Frain=F_rain_phy(I,L) ! 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 QRdum=QW QW=0. ELSE QRdum=Frain*QW QW=QW-QRdum ENDIF ENDIF IF (QI .LE. 0.) F_RimeF_phy(I,L)=1. RimeF_col(L)=F_RimeF_phy(I,L) ! (real) QI_col(L)=QI QR_col(L)=QRdum 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 !mz 05/06/2020 - 10km !------------------------------------------------------------ IF(DX1 .GE. 10000 .AND. P_col(L)0) associated with snow ! APREC(I)=(ARAIN+ASNOW)*RRHOL ! Accumulated surface precip (depth in m) !<--- Ying PREC(I)=PREC(I)+APREC(I) ACPREC(I)=ACPREC(I)+APREC(I) IF(APREC(I) .LT. 1.E-8) THEN SR(I)=0. ELSE SR(I)=RRHOL*ASNOW/APREC(I) ENDIF ! !####################################################################### !####################################################################### ! enddo ! End "I" loop ! !----------------------------------------------------------------------- !-- End of original driver for EGCP01COLUMN_hr !----------------------------------------------------------------------- ! do k = lm, 1, -1 DO i = ims,ime WC=qt(I,K) QS(I,K)=0. QR(I,K)=0. QC(I,K)=0. ! IF(F_ICE_PHY(I,K)>=1.)THEN QS(I,K)=WC ELSEIF(F_ICE_PHY(I,K)<=0.)THEN QC(I,K)=WC ELSE QS(I,K)=F_ICE_PHY(I,K)*WC QC(I,K)=WC-QS(I,K) ENDIF ! IF(QC(I,K)>0..AND.F_RAIN_PHY(I,K)>0.)THEN IF(F_RAIN_PHY(I,K).GE.1.)THEN QR(I,K)=QC(I,K) QC(I,K)=0. ELSE QR(I,K)=F_RAIN_PHY(I,K)*QC(I,K) QC(I,K)=QC(I,K)-QR(I,K) ENDIF ENDIF ENDDO !- i ENDDO !- k ! !- Update rain (convert from m to kg/m**2, which is also equivalent to mm depth) ! DO i=ims,ime RAINNC(i)=APREC(i)*1000.+RAINNC(i) RAINNCV(i)=APREC(i)*1000. ENDDO ! !----------------------------------------------------------------------- ! END SUBROUTINE FER_HIRES ! !----------------------------------------------------------------------- ! !############################################################################### ! ***** 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. ! (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. There are two modes of ice nucleation ! following Meyers et al. (JAM, 1992): ! a) Deposition & condensation freezing nucleation - eq. (2.4) when ! air is supersaturated w/r/t ice ! b) Contact freezing nucleation - eq. (2.6) in presence of cloud water ! (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). ! (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. !############################################################################### !############################################################################### ! !>\ingroup hafs_famp !! This is the grid-scale microphysical processes of Ferrier-Aligo microphysics !! scheme (i.e., condensation and precipitation). !!\param arain accumulated rainfall at the surface (kg) !!\param asnow accumulated snowfall at the surface (kg) !!\param dtph physics time step (s) !!\param rhc_col vertical column of threshold relative humidity for onset of !! condensation (ratio) !!\param i_index i index !!\param j_index j index !!\param lsfc Eta level of level above surface, ground !!\param p_col vertical column of model pressure (Pa) !!\param qi_col vertical column of model ice mixing ratio (kg/kg) !!\param qr_col vertical column of model rain ratio (kg/kg) !!\param q_col vertical column of model water vapor specific humidity (kg/kg) !!\param qw_col vertical column of model cloud water mixing ratio (kg/kg) !!\param rimef_col vertical column of rime factor for ice in model (ratio, defined below) !!\param t_col vertical column of model temperature (deg K) !!\param thick_col vertical column of model mass thickness (density*height increment) !!\param wc_col vertical column of model mixing ratio of total condensate (kg/kg) !!\param lm vertical dimension !!\param pcond1d net cloud water condensation (>0) or evaporation (<0) (kg/kg) !!\param pidep1d net ice deposition (>0) or sublimation (<0) (kg/kg) !!\param piacw1d cloud water collection by precipitation ice (kg/kg) !!\param piacwi1d cloud water riming onto precipitation ice at <0 (kg/kg) !!\param piacwr1d accreted cloud water shed to form rain at >0 (kg/kg) !!\param piacr1d freezing of supercooled rain to precipitation ice (kg/kg) !!\param picnd1d condensation onto wet, melting ice (kg/kg) !!\param pievp1d evaporation from wet, melting ice (kg/kg) !!\param pimlt1d melting of precipitation ice to form rain (kg/kg) !!\param praut1d droplet self_collection (autoconversion) to form rain (kg/kg) !!\param pracw1d cloud water collection (accretion) by rain (kg/kg) !!\param prevp1d rain evaporation (kg/kg) !!\param pisub1d !!\param pevap1d !!\param DBZ_col vertical column of radar reflectivity (dBZ) !!\param NR_col vertical column of rain number concentration (m^-3) !!\param NS_col vertical column of snow number concentration (m^-3) !!\param vsnow1d fall speed of rimed snow w/ air resistance correction !!\param vrain11d fall speed of rain into grid from above (m/s) !!\param vrain21d fall speed of rain out of grid box to the level below (m/s) !!\param vci1d Fall speed of 50-micron ice crystals w/ air resistance correction !!\param NSmICE1d number concentration of small ice crystals at current level !!\param INDEXS1d !!\param INDEXR1d !!\param RFlag1d !!\param DX1 SUBROUTINE EGCP01COLUMN_hr ( ARAIN, ASNOW, DTPH, RHC_col, & & I_index, J_index, LSFC, & & P_col, QI_col, QR_col, Q_col, QW_col, RimeF_col, T_col, & & THICK_col, WC_col ,LM,pcond1d,pidep1d, & & piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d,pimlt1d, & & praut1d,pracw1d,prevp1d,pisub1d,pevap1d, DBZ_col,NR_col,NS_col, & & vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d,INDEXS1d,INDEXR1d, & !jul28 & RFlag1d,DX1) !jun01 ! !############################################################################### !############################################################################### ! !------------------------------------------------------------------------------- !----- 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_hr FROM SUBROUTINE EGCP01DRV ! ! INPUT ARGUMENT LIST: ! DTPH - physics time step (s) ! RHgrd - threshold relative humidity (ratio) for onset of condensation ! 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) ! Q_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) ! Q_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) ! DBZ_col - vertical column of radar reflectivity (dBZ) ! NR_col - vertical column of rain number concentration (m^-3) ! NS_col - vertical column of snow number concentration (m^-3) ! ! OUTPUT FILES: ! NONE ! ! Subprograms & Functions called: ! * Real Function CONDENSE - cloud water condensation ! * Real Function DEPOSIT - ice deposition (not sublimation) ! * Integer Function GET_INDEXR - estimate the mean size of raindrops (microns) ! ! UNIQUE: NONE ! ! LIBRARY: NONE ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 90 ! MACHINE : IBM SP ! !------------------------------------------------------------------------- !--------------- Arrays & constants in argument list --------------------- !------------------------------------------------------------------------- ! IMPLICIT NONE ! INTEGER,INTENT(IN) :: LM,I_index, J_index, LSFC,DX1 REAL,INTENT(IN) :: DTPH REAL,INTENT(INOUT) :: ARAIN, ASNOW REAL,DIMENSION(LM),INTENT(INOUT) :: P_col, QI_col,QR_col & & ,Q_col ,QW_col, RimeF_col, T_col, THICK_col,WC_col,pcond1d & & ,pidep1d,piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d & & ,pimlt1d,praut1d,pracw1d,prevp1d,pisub1d,pevap1d,DBZ_col,NR_col & & ,NS_col,vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d,INDEXS1d & !jun01 & ,INDEXR1d,RFlag1d,RHC_col !jun01 ! !-------------------------------------------------------------------------------- !--- The following arrays are integral calculations based on the mean ! snow/graupel diameters, which vary from 50 microns to 1000 microns ! (1 mm) at 1-micron intervals and assume exponential size distributions. ! The values are normalized and require being multipled by the number ! concentration of large ice (NLICE). !--------------------------------------- ! - 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). !-------------------------------------------------------------------------------- !--- The following arrays are integral calculations based on the mean ! rain diameters, which vary from 50 microns to 1000 microns ! (1 mm) at 1-micron intervals and assume exponential size distributions. ! The values are normalized and require being multiplied by the rain intercept ! (N0r). !--------------------------------------- ! - VENTR1 - integrated quantity associated w/ ventilation effects ! (capacitance only) for calculating evaporation from rain ! - VENTR2 - integrated quantity associated w/ ventilation effects ! (with fall speed) for calculating evaporation from rain ! - ACCRR - integrated quantity associated w/ cloud water collection by rain ! - MASSR - integrated quantity associated w/ rain ! - VRAIN - mass-weighted fall speed of rain, used to calculate ! precipitation rates ! - RRATE - precipitation rates, which should also be equal to RHO*QR*VRAIN ! !------------------------------------------------------------------------- !------- 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, Zmin=0.01, DBZmin=-20. ! !--- 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. ! !--- This variable is for debugging purposes (if .true.) ! LOGICAL, PARAMETER :: PRINT_diag=.false. ! !----------------------------------------------------------------------- !--- Local variables !----------------------------------------------------------------------- ! REAL :: EMAIRI, N0r, NLICE, NSmICE, NInuclei, Nrain, Nsnow, Nmix REAL :: RHgrd LOGICAL :: CLEAR, ICE_logical, DBG_logical, RAIN_logical, & STRAT, DRZL INTEGER :: INDEX_MY,INDEXR,INDEXR1,INDEXR2,INDEXS,IPASS,ITDX,IXRF,& & IXS,LBEF,L,INDEXRmin,INDEXS0C,IDR !mar03 !may20 ! ! REAL :: ABI,ABW,AIEVP,ARAINnew,ASNOWnew,BLDTRH,BUDGET, & & CREVP,DELI,DELR,DELT,DELV,DELW,DENOMF, & & DENOMI,DENOMW,DENOMWI,DIDEP, & & DIEVP,DIFFUS,DLI,DTRHO,DUM,DUM1,DUM2,DUM3, & & DWV0,DWVI,DWVR,DYNVIS,ESI,ESW,FIR,FLIMASS, & & FWR,FWS,GAMMAR,GAMMAS, & & PCOND,PIACR,PIACW,PIACWI,PIACWR,PICND,PIDEP,PIDEP_max, & & PIEVP,PILOSS,PIMLT,PINIT,PP,PRACW,PRAUT,PREVP,PRLOSS, & & QI,QInew,QLICE,QR,QRnew,QSI,QSIgrd,QSInew,QSW,QSW0, & & QSWgrd,QSWnew,QT,QTICE,QTnew,QTRAIN,Q,QW,QWnew,Rcw, & & RFACTOR,RFmx,RFrime,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, & & VSNOW1,WC,WCnew,WSgrd,WS,WSnew,WV,WVnew, & & XLI,XLIMASS,XRF, & & NSImax,QRdum,QSmICE,QLgIce,RQLICE,VCI,TIMLT, & & RQSnew,RQRnew,Zrain,Zsnow,Ztot,RHOX0C,RFnew,PSDEP,DELS !mar03 !apr22 REAL, SAVE :: Revised_LICE=1.e-3 !-- kg/m**3 ! !####################################################################### !########################## Begin Execution ############################ !####################################################################### ! ! ARAIN=0. ! Accumulated rainfall into grid box from above (kg/m**2) VRAIN1=0. ! Rain fall speeds into grib box from above (m/s) VSNOW1=0. ! Ice fall speeds into grib box from above (m/s) ASNOW=0. ! Accumulated snowfall into grid box from above (kg/m**2) INDEXS0C=MDImin ! Mean snow/graupel diameter just above (<0C) freezing level (height) RHOX0C=22.5 ! Estimated ice density at 0C (kg m^-3) !mar03 TIMLT=0. ! Total ice melting in a layer (drizzle detection) STRAT=.FALSE. ! Stratiform rain DSD below melting level !may11 DRZL=.FALSE. ! Drizzle DSD below melting level !may23 ! !----------------------------------------------------------------------- !------------ Loop from top (L=1) to surface (L=LSFC) ------------------ !----------------------------------------------------------------------- ! big_loop: DO L=LM,1,-1 pcond1d(L)=0. pidep1d(L)=0. piacw1d(L)=0. piacwi1d(L)=0. piacwr1d(L)=0. piacr1d(L)=0. picnd1d(L)=0. pievp1d(L)=0. pimlt1d(L)=0. praut1d(L)=0. pracw1d(L)=0. prevp1d(L)=0. pisub1d(L)=0. pevap1d(L)=0. vsnow1d(L)=0. vrain11d(L)=0. vrain21d(L)=0. vci1d(L)=0. NSmICE1d(L)=0. DBZ_col(L)=DBZmin NR_col(L)=0. NS_col(L)=0. INDEXR1d(L)=0. INDEXS1d(L)=0. RFlag1d(L)=0. !jun01 ! !--- Skip this level and go to the next lower level if no condensate ! and very low specific humidities ! !--- Check if any rain is falling into layer from above ! IF (ARAIN .GT. CLIMIT) THEN CLEAR=.FALSE. VRAIN1=0. ELSE CLEAR=.TRUE. ARAIN=0. ENDIF ! !--- Check if any ice is falling into layer from above ! !--- NOTE that "SNOW" in variable names is often synonomous with ! large, precipitation ice particles ! IF (ASNOW .GT. CLIMIT) THEN CLEAR=.FALSE. VSNOW1=0. ELSE ASNOW=0. ENDIF ! !----------------------------------------------------------------------- !------------ Proceed with cloud microphysics calculations ------------- !----------------------------------------------------------------------- ! TK=T_col(L) ! Temperature (deg K) TC=TK-T0C ! Temperature (deg C) PP=P_col(L) ! Pressure (Pa) Q=Q_col(L) ! Specific humidity of water vapor (kg/kg) WV=Q/(1.-Q) ! 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 !----------------------------------------------------------------------- ! !--- This check is to determine grid-scale saturation when no condensate is present ! ESW=MIN(1000.*FPVS0(TK),0.99*PP) ! Saturation vapor pressure w/r/t water QSW=EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water WS=QSW ! General saturation mixing ratio (water/ice) QSI=QSW ! Saturation mixing ratio w/r/t ice IF (TC .LT. 0.) THEN ESI=MIN(1000.*FPVS(TK),0.99*PP) ! 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. ! !----------------------------------------------------------------------- !-- Loop to the end if in clear, subsaturated air free of condensate --- !----------------------------------------------------------------------- ! IF (CLEAR) THEN STRAT=.FALSE. !- Reset stratiform rain flag DRZL=.FALSE. !- Reset drizzle flag INDEXRmin=MDRmin !- Reset INDEXRmin TIMLT=0. !- Reset accumulated ice melting CYCLE big_loop ENDIF ! !----------------------------------------------------------------------- !--------- 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*Q)) ! 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) PINIT=0. ! Ice initiation (part of PIDEP calculation, 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) NSmICE=0. ! Cloud ice number concentration (m^-3) Nrain=0. ! Rain number concentration (m^-3) !jul28 begin Nsnow=0. ! "Snow" number concentration (m^-3) RQRnew=0. ! Final rain content (kg/m**3) RQSnew=0. ! Final "snow" content (kg/m**3) Zrain=0. ! Radar reflectivity from rain (mm**6 m-3) Zsnow=0. ! Radar reflectivity from snow (mm**6 m-3) Ztot=0. ! Radar reflectivity from rain+snow (mm**6 m-3) INDEXR=MDRmin ! Mean diameter of rain (microns) INDEXR1=INDEXR ! 1st updated mean diameter of rain (microns) INDEXR2=INDEXR ! 2nd updated mean diameter of rain (microns) N0r=0. ! 1st estimate for rain intercept (m^-4) DUM1=MIN(0.,TC) DUM=XMImax*EXP(XMIexp*DUM1) INDEXS=MIN(MDImax, MAX(MDImin, INT(DUM) ) ) ! 1st estimate for mean diameter of snow (microns) VCI=0. ! Cloud ice fall speeds (m/s) VSNOW=0. ! "Snow" (snow/graupel/sleet/hail) fall speeds (m/s) VRAIN2=0. ! Rain fall speeds out of bottom of grid box (m/s) RimeF1=1. ! Rime Factor (ratio, >=1, defined below) ! !--- 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(0,"(/2(a,i4),a,i2)") '{@ i=',I_index,' j=',J_index, ! & ' L=',L ! WRITE(0,"(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) ! TK2=1./(TK*TK) ! 1./TK**2 ! !--- 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=SQRT(TK*TK*TK)/(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=MIN(1.5, (1.E5/PP)**C1) !-- limited to 1.5x ! !--- 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 IF (T_ICE <= -100.) THEN 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. ! ice_test: IF (ICE_logical) THEN ! !--- IMPORTANT: Estimate time-averaged properties. ! !--- ! -> Small ice particles are assumed to have a mean diameter of 50 microns. ! * QSmICE - estimated mixing ratio for small cloud ice !--- ! * 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 ! * VCI - Fall speed of 50-micron ice crystals w/ air resistance correction ! * EMAIRI - equivalent mass of air associated layer and with fall of snow into layer ! * XLIMASS - used for debugging, associated with 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 ! * QSmICE - mixing ratio 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. ! NInuclei=0. NSmICE=0. QSmICE=0. Rcw=0. IF (TC<0.) THEN ! !--- Max # conc of small ice crystals based on 10% of total ice content ! or the parameter NSI_max ! NSImax=MAX(NSI_max, 0.1*RHO*QI/MASSI(MDImin) ) !aug27 ! !-- Specify Fletcher, Cooper, Meyers, etc. here for ice nuclei concentrations ! Cooper (1986): NInuclei=MIN(5.*EXP(-0.304*TC), NSImax) ! Fletcher (1962): NInuclei=MIN(0.01*EXP(-0.6*TC), NSImax) ! !aug28: The formulas below mean that Fletcher is used for >-21C and Cooper at colder ! temperatures. In areas of high ice contents near the tops of deep convection, ! the number concentrations will be determined by the lower value of the "FQi" ! contribution to NSImax or Cooper. ! NInuclei=MIN(0.01*EXP(-0.6*TC), NSImax) !aug28 - Fletcher (1962) NInuclei=MIN(5.*EXP(-0.304*TC), NInuclei) !aug28 - Cooper (1984) IF (QI>EPSQ) THEN DUM=RRHO*MASSI(MDImin) NSmICE=MIN(NInuclei, QI/DUM) QSmICE=NSmICE*DUM ENDIF ! End IF (QI>EPSQ) ENDIF ! End IF (TC<0.) init_ice: IF (QI<=EPSQ .AND. ASNOW<=CLIMIT) THEN TOT_ICE=0. PILOSS=0. RimeF1=1. VrimeF=1. VEL_INC=GAMMAS VSNOW=0. VSNOW1=0. VCI=0. EMAIRI=THICK XLIMASS=RimeF1*MASSI(INDEXS) FLIMASS=1. QLICE=0. RQLICE=0. QTICE=0. NLICE=0. ELSE init_ice ! !--- 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). ! ! !sep10 - Start of changes described in (23) at top of code. ! TOT_ICE=THICK*QI+BLEND*ASNOW PILOSS=-TOT_ICE/THICK QLgICE=MAX(0., QI-QSmICE) !-- 1st-guess estimate of large ice VCI=GAMMAS*VSNOWI(MDImin) ! !-- Need to save this original value before two_pass iteration ! LBEF=MAX(1,L-1) RimeF1=(RimeF_col(L)*THICK*QLgICE & & +RimeF_col(LBEF)*BLEND*ASNOW)/TOT_ICE ! !mar08 see (32); !apr22a see (41) start - Estimate mean diameter (INDEXS in microns) IF (RimeF1>2.) THEN DUM3=RH_NgC*(RHO*QLgICE)**C1 !- convective mean diameter DUM2=RH_NgT*(RHO*QLgICE)**C1 !- transition mean diameter IF (RimeF1>=10.) THEN DUM=DUM3 ELSE IF (RimeF1>=5.) THEN DUM=0.2*(RimeF1-5.) !- Blend at 5<=RF<10 DUM=DUM3*DUM+DUM2*(1.-DUM) ELSE DUM1=REAL(INDEXS) !- stratiform mean diameter DUM=0.33333*(RimeF1-2.) !- Blend at 2=5. .AND. INDEXS==MDImax .AND. RQLICE>RQhail) THEN !- Additional increase using Thompson graupel/hail fall speeds DUM=GAMMAS*AVhail*RQLICE**BVhail IF (DUM>VSNOW) THEN VSNOW=DUM VEL_INC=VSNOW/VSNOWI(INDEXS) ENDIF ENDIF XLIMASS=RimeF1*MASSI(INDEXS) NLICE=RQLICE/XLIMASS ! !sep16 - End of change described in (26) at top of code. ! !=========================================== IF (IPASS>=2 .OR. & (NLICE>=NLImin .AND. NLICE<=NSI_max)) EXIT two_pass !may17 - end !=========================================== ! !--- Force NLICE to be between NLImin and NSI_max when IPASS=1 ! ! IF (PRINT_diag .AND. RQLICE>Revised_LICE) DUM2=NLICE !-- For debugging (see DUM2 below) NLICE=MAX(NLImin, MIN(NSI_max, NLICE) ) !sep16 - End of changes ! XLI=RQLICE/(NLICE*RimeF1) !- Mean mass of unrimed ice new_size: IF (XLI<=MASSI(MDImin) ) THEN INDEXS=MDImin ELSE IF (XLI<=MASSI(450) ) THEN new_size DLI=9.5885E5*XLI**.42066 ! DLI in microns INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) ) ELSE IF (XLIRevised_LICE) THEN ! WRITE(0,"(5(a12,g11.4,1x))") '{$ RimeF1=',RimeF1, & ! & ' RHO*QLICE=',RQLICE,' TC=',TC,' NLICE=',NLICE, & ! & ' NLICEold=',DUM2 ! Revised_LICE=1.2*RQLICE ! ENDIF ENDIF new_size !=========================================== ENDDO two_pass !=========================================== ENDIF init_ice ENDIF ice_test ! !mar03 !may11 !jun01 - start; see (34) above IF(TC<=0.) THEN STRAT=.FALSE. INDEXRmin=MDRmin TIMLT=0. INDEXS0C=INDEXS RHOX0C=22.5*MAX(1.,MIN(RimeF1,40.)) !- Estimated ice density at 0C (kg m^-3) ELSE ! TC>0. IF(.NOT.RAIN_logical) THEN STRAT=.FALSE. !- Reset STRAT INDEXRmin=MDRmin !- Reset INDEXRmin IF(.NOT.ICE_logical) TIMLT=0. ELSE !- STRAT=T for stratiform rain IF(TIMLT>EPSQ .AND. RHOX0C<=225.) THEN STRAT=.TRUE. ELSE STRAT=.FALSE. !- Reset STRAT INDEXRmin=MDRmin !- Reset INDEXRmin ENDIF IF(STRAT .AND. INDEXRmin<=MDRmin) THEN INDEXRmin=INDEXS0C*(0.001*RHOX0C)**C1 INDEXRmin=MAX(MDRmin, MIN(INDEXRmin, INDEXRstrmax) ) ENDIF ENDIF ENDIF ! IF(STRAT .OR. TIMLT>EPSQ) THEN DRZL=.FALSE. ELSE !- DRZL=T for drizzle (no melted ice falling from above) DRZL=.TRUE. !mar30 ENDIF !jun01 - end ! !---------------------------------------------------------------------- !--------------- Calculate individual processes ----------------------- !---------------------------------------------------------------------- ! !--- Cloud water autoconversion to rain (PRAUT) and collection of cloud ! water by precipitation ice (PIACW) ! IF (QW.GT.EPSQ .AND. TC.GE.T_ICE) THEN !-- The old autoconversion threshold returns DUM2=RHO*QW IF (DUM2>QAUT0) THEN !-- July 2010 version follows Liu & Daum (JAS, 2004) and Liu et al. (JAS, 2006) DUM2=DUM2*DUM2 DUM=BRAUT*DUM2*QW DUM1=ARAUT*DUM2 PRAUT=MIN(QW, DUM*(1.-EXP(-DUM1*DUM1)) ) ENDIF 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) ) !jul28 (16) PIACW=FWS*QW IF (TC<0.) THEN PIACWI=PIACW !- Large ice riming Rcw=ARcw*DUM2**C1 !- Cloud droplet radius in microns ENDIF ENDIF ! End IF (QLICE .GT. EPSQ) ENDIF ! End IF (QW.GT.EPSQ .AND. TC.GE.T_ICE) ! !---------------------------------------------------------------------- !--- Calculate homogeneous freezing of cloud water (PIACW, PIACWI) and ! ice deposition (PIDEP), which also includes ice initiation (PINIT) ! ice_only: IF (TC.LT.T_ICE .AND. (WV.GT.QSWgrd .OR. QW.GT.EPSQ)) THEN ! !--- Adjust to ice saturation at T10): ! RH_NgC=PI*500.*5.E3 ! RHOg=500 kg m^-3, Ng=5.e3 m^-3 => "hail" content >7.85 g m^-3 for INDEXS=MDImax !- or - ! RH_NgC=PI*500.*1.E3 ! RHOg=500 kg m^-3, Ng=1.e3 m^-3 => "hail" content >1.57 g m^-3 for INDEXS=MDImax ! RH_NgC=PI*500.*1.E3 !- RHOg=500 kg m^-3, Ng=1.e3 m^-3 RQhail=RH_NgC*(1.E-3)**3 !- Threshold "hail" content ! becomes 1.57 g m^-3 Bvhail=0.82*C1 !- Exponent for Thompson graupel Avhail=1353.*(1./RH_NgC)**Bvhail !- 1353 ~ constant for Thompson graupel RH_NgC=1.E6*(1./RH_NgC)**C1 !mar08 (convection, convert to microns) ! !- An explanation for the following settings assumed for graupel (RF>5): ! RH_NgT=PI*300.*10.E3 ! RHOg=300 kg m^-3, Ng=10.e3 m^-3 => "graupel" content must exceed 9.43 g m^-3 for INDEXS=MDImax !- or - ! RH_NgT=PI*300.*5.E3 ! RHOg=300 kg m^-3, Ng=5.e3 m^-3 => "graupel" content must exceed 4.71 g m^-3 for INDEXS=MDImax ! RH_NgT=PI*300.*5.E3 !- RHOg=300 kg m^-3, Ng=5.e3 m^-3 RH_NgT=1.E6*(1./RH_NgT)**C1 !mar08 (transition, convert to microns) !may20 - end ! !--- 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. ! RETURN ! !----------------------------------------------------------------------- END SUBROUTINE FERRIER_INIT_hr ! !>\ingroup hafs_famp SUBROUTINE MY_GROWTH_RATES_NMM_hr (DTPH) ! !--- 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 ! !----------------------------------------------------------------------- !-- 20090714: These values are in g and need to be converted to kg below 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_NMM=DT_ICE*MY_600*1.E-3 !-- 20090714: Convert from g to kg ! !----------------------------------------------------------------------- ! END SUBROUTINE MY_GROWTH_RATES_NMM_hr ! !----------------------------------------------------------------------- !--------- Old GFS saturation vapor pressure lookup tables ----------- !----------------------------------------------------------------------- ! !>\ingroup hafs_famp SUBROUTINE GPVS_hr ! ****************************************************************** !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . ! SUBPROGRAM: GPVS_hr 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_hr ! ! SUBPROGRAMS CALLED: ! (FPVSX) - INLINABLE FUNCTION TO COMPUTE SATURATION VAPOR PRESSURE ! ! 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_hr !----------------------------------------------------------------------- !*********************************************************************** !----------------------------------------------------------------------- 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 !----------------------------------------------------------------------- IF (T>=XMIN .AND. T<=XMAX) THEN XJ=MIN(MAX(C1XPVS+C2XPVS*T,1.),FLOAT(NX)) JX=MIN(XJ,NX-1.) FPVS=TBPVS(JX)+(XJ-JX)*(TBPVS(JX+1)-TBPVS(JX)) ELSE IF (T>XMAX) THEN !-- Magnus Tetens formula for water saturation (Murray, 1967) ! (saturation vapor pressure in kPa) FPVS=0.61078*exp(17.2694*(T-273.16)/(T-35.86)) ELSE !-- Magnus Tetens formula for ice saturation(Murray, 1967) ! (saturation vapor pressure in kPa) FPVS=0.61078*exp(21.8746*(T-273.16)/(T-7.66)) ENDIF ! END FUNCTION FPVS !----------------------------------------------------------------------- !----------------------------------------------------------------------- REAL FUNCTION FPVS0(T) !----------------------------------------------------------------------- IMPLICIT NONE real,INTENT(IN) :: T real :: XJ1 integer :: JX1 !----------------------------------------------------------------------- IF (T>=XMIN .AND. T<=XMAX) THEN XJ1=MIN(MAX(C1XPVS0+C2XPVS0*T,1.),FLOAT(NX)) JX1=MIN(XJ1,NX-1.) FPVS0=TBPVS0(JX1)+(XJ1-JX1)*(TBPVS0(JX1+1)-TBPVS0(JX1)) ELSE !-- Magnus Tetens formula for water saturation (Murray, 1967) ! (saturation vapor pressure in kPa) FPVS0=0.61078*exp(17.2694*(T-273.16)/(T-35.86)) ENDIF ! 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 SUBROUTINE ferhires_finalize() IMPLICIT NONE if (ALLOCATED(ventr1)) DEALLOCATE(ventr1) if (ALLOCATED(ventr2)) DEALLOCATE(ventr2) if (ALLOCATED(accrr)) DEALLOCATE(accrr) if (ALLOCATED(massr)) DEALLOCATE(massr) if (ALLOCATED(vrain)) DEALLOCATE(vrain) if (ALLOCATED(rrate)) DEALLOCATE(rrate) if (ALLOCATED(venti1)) DEALLOCATE(venti1) if (ALLOCATED(venti2)) DEALLOCATE(venti2) if (ALLOCATED(accri)) DEALLOCATE(accri) if (ALLOCATED(massi)) DEALLOCATE(massi) if (ALLOCATED(vsnowi)) DEALLOCATE(vsnowi) if (ALLOCATED(vel_rf)) DEALLOCATE(vel_rf) END SUBROUTINE ferhires_finalize ! END MODULE module_mp_fer_hires