
!#######################################################################
!------- Initialize constants & lookup tables for microphysics ---------
!#######################################################################

    SUBROUTINE GSMCONST (DTPH)

!-------------------------------------------------------------------------------
!---  SUBPROGRAM DOCUMENTATION BLOCK
!   PRGRMMR: Ferrier         ORG: W/NP22     DATE: February 2001
!-------------------------------------------------------------------------------
! 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 GSMCONST 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

!-----------------------------------------------------------------------

!--- INCLUDE COMMON BLOCKS

    INCLUDE "parmeta.f90"
    INTEGER, PARAMETER :: LP1=LM+1
    INCLUDE "COMM_MAPOT.f90"

!-------------------------------------------------------------------------
!-------------- Parameters & arrays for lookup tables --------------------
!-------------------------------------------------------------------------

!--- Common block of constants used in column microphysics

    COMMON /CMICRO_CONS/ ABFR, CBFR, CIACW, CIACR, C_N0r0, &
    CN0r0, CN0r_DMRmin, CN0r_DMRmax, CRACW, CRAUT, ESW0, &
    QAUT0, RFmax, RHgrd, RQR_DR1, RQR_DR2, RQR_DR3, RQR_DRmin, &
    RQR_DRmax, RR_DRmin, RR_DR1, RR_DR2, RR_DR3, RR_DRmax

!--- Common block for lookup table used in calculating growth rates of
!    nucleated ice crystals growing in water saturated conditions

    INTEGER, PARAMETER :: MY_T1=1, MY_T2=35
    COMMON /CMY600/ MY_GROWTH(MY_T1:MY_T2)
    REAL :: MY_GROWTH

!--- Mean ice particle diameters vary from 50 microns to 1000 microns

    REAL, PARAMETER :: DMImin=.05e-3, DMImax=1.e-3, DelDMI=1.e-6, &
    XMImin=1.e6*DMImin, XMImax=1.e6*DMImax
    INTEGER, PARAMETER :: MDImin=XMImin, MDImax=XMImax

!--- Various ice lookup tables

    COMMON /IACCR_TABLES/ ACCRI(MDImin:MDImax)
    COMMON /IMASS_TABLES/ MASSI(MDImin:MDImax)
    COMMON /SDENS_TABLES/ SDENS(MDImin:MDImax)
    REAL :: MASSI
    COMMON /IRATE_TABLES/ VSNOWI(MDImin:MDImax)
    COMMON /IVENT_TABLES/ VENTI1(MDImin:MDImax), VENTI2(MDImin:MDImax)

!--- Common block for riming tables

    INTEGER, PARAMETER :: Nrime=40
    COMMON /IRIME_TABLES/ VEL_RF(2:9,0:Nrime)

!--- Mean rain drop diameters vary from 50 microns to 450 microns

    REAL, PARAMETER :: DMRmin=.05E-3, DMRmax=.45E-3, DelDMR=1.E-6, &
    XMRmin=1.E6*DMRmin, XMRmax=1.E6*DMRmax
    INTEGER, PARAMETER :: MDRmin=XMRmin, MDRmax=XMRmax

!--- Various rain lookup tables

    COMMON /RACCR_TABLES/ ACCRR(MDRmin:MDRmax)
    COMMON /RMASS_TABLES/ MASSR(MDRmin:MDRmax)
    REAL :: MASSR
    COMMON /RRATE_TABLES/ RRATE(MDRmin:MDRmax)
    COMMON /RVELR_TABLES/ VRAIN(MDRmin:MDRmax)
    COMMON /RVENT_TABLES/ VENTR1(MDRmin:MDRmax), VENTR2(MDRmin:MDRmax)

!--- Parameters & data statement for local calculations

    REAL, PARAMETER :: C1=1./3., DMR1=.1E-3, DMR2=.2E-3, DMR3=.32E-3, &
    N0r0=8.E6, N0s0=4.E6, RHOL=1000., RHOS=100., T0C=273.15, &
    XMR1=1.e6*DMR1, XMR2=1.e6*DMR2, XMR3=1.e6*DMR3
    INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3

!------------------------------------------------------------------------
!--------- Constants passed through /CMICRO_CONS/ common block ----------
!------------------------------------------------------------------------

!--- 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+0.08*[(100.-dX)/95.]**.5

    DX=111.*(DPHD**2+DLMD**2)**.5         ! Model resolution at equator (km)
    DX=MIN(100., MAX(5., DX) )
!    RHgrd=0.90+.08*((100.-DX)/95.)**.5
! OU      RHgrd=0.998
      RHgrd=1.05

!--- Create lookup tables for saturation vapor pressure w/r/t water & ice

    CALL GPVS

!--- Read in various lookup tables

    OPEN (UNIT=1,FILE="eta_micro_lookup.dat",FORM="UNFORMATTED")
    READ(1) VENTR1
    READ(1) VENTR2
    READ(1) ACCRR
    READ(1) MASSR
    READ(1) VRAIN
    READ(1) RRATE
    READ(1) VENTI1
    READ(1) VENTI2
    READ(1) ACCRI
    READ(1) MASSI
    READ(1) VSNOWI
! LG
! SM	VSNOWI=VSNOWI*1.2
! SM	VSNOWI=VSNOWI*1.1
    VSNOWI=VSNOWI*1.2
! LG
    READ(1) VEL_RF
!      read(1) my_growth    ! Applicable only for DTPH=180 s
    CLOSE (1)

!--- Calculates coefficients for growth rates of ice nucleated in water
!    saturated conditions, scaled by physics time step (lookup table)

    CALL MY_GROWTH_RATES (DTPH)

    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-42

!--- 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

    XNCW=200.E6                 ! 300 cm**-3 droplet concentration
    QAUT0=PI*RHOL*XNCW*(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

!-----------------------------------------------------------------------

    RETURN
    END SUBROUTINE GSMCONST

!#######################################################################
!--- Sets up lookup table for calculating initial ice crystal growth ---
!#######################################################################

    SUBROUTINE MY_GROWTH_RATES (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).

    integer, parameter :: MY_T1=1, MY_T2=35
    COMMON /CMY600/ MY_GROWTH(MY_T1:MY_T2)
    REAL :: MY_GROWTH, MY_600(MY_T1:MY_T2)

    DATA MY_600 / &
    &  5.5e-8, 1.4E-7, 2.8E-7, 6.E-7, 3.3E-6,      & !  -1 to  -5 deg C
    &  2.E-6, 9.E-7, 8.8E-7, 8.2E-7, 9.4e-7,       & !  -6 to -10 deg C
    &  1.2E-6, 1.85E-6, 5.5E-6, 1.5E-5, 1.7E-5,    & ! -11 to -15 deg C
    &  1.5E-5, 1.E-5, 3.4E-6, 1.85E-6, 1.35E-6,    & ! -16 to -20 deg C
    &  1.05E-6, 1.E-6, 9.5E-7, 9.0E-7, 9.5E-7,     & ! -21 to -25 deg C
    &  9.5E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7,         & ! -26 to -30 deg C
    &  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

!-----------------------------------------------------------------------

    RETURN
    END SUBROUTINE MY_GROWTH_RATES

!#######################################################################
!-- Lookup tables for the saturation vapor pressure w/r/t water & ice --
!#######################################################################

    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

! 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
!   MACHINE:  IBM SP

!$$$
!----------------------------------------------------------------------
    PARAMETER(NX=7501)
    REAL :: TBPVS(NX),TBPVS0(NX)
    COMMON/COMPVS0/ C1XPVS0,C2XPVS0,TBPVS0
    COMMON/COMPVS/ C1XPVS,C2XPVS,TBPVS
!----------------------------------------------------------------------
    XMIN=180.0
    XMAX=330.0
    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

    RETURN
    END SUBROUTINE GPVS
!-----------------------------------------------------------------------
!***********************************************************************
!-----------------------------------------------------------------------
    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

! USAGE:   PVS=FPVS(T)

!   INPUT ARGUMENT LIST:
!     T        - REAL TEMPERATURE IN KELVIN

!   OUTPUT ARGUMENT LIST:
!     FPVS     - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)

! COMMON BLOCKS:
!   COMPVS   - SCALING PARAMETERS AND TABLE COMPUTED IN GPVS.

! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE:  IBM SP

!$$$
!-----------------------------------------------------------------------
    PARAMETER(NX=7501)
    REAL :: TBPVS(NX)
    COMMON/COMPVS/ C1XPVS,C2XPVS,TBPVS
!-----------------------------------------------------------------------
    XJ=MIN(MAX(C1XPVS+C2XPVS*T,1.),FLOAT(NX))
    JX=MIN(XJ,NX-1.)
    FPVS=TBPVS(JX)+(XJ-JX)*(TBPVS(JX+1)-TBPVS(JX))

    RETURN
    END FUNCTION FPVS
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
    FUNCTION FPVS0(T)
!-----------------------------------------------------------------------
    PARAMETER(NX=7501)
    REAL :: TBPVS0(NX)
    COMMON/COMPVS0/ C1XPVS0,C2XPVS0,TBPVS0
!-----------------------------------------------------------------------
    XJ1=MIN(MAX(C1XPVS0+C2XPVS0*T,1.),FLOAT(NX))
    JX1=MIN(XJ1,NX-1.)
    FPVS0=TBPVS0(JX1)+(XJ1-JX1)*(TBPVS0(JX1+1)-TBPVS0(JX1))

    RETURN
    END FUNCTION FPVS0
!-----------------------------------------------------------------------
!***********************************************************************
!-----------------------------------------------------------------------
    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

! 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
!   MACHINE:  IBM SP

!$$$
!-----------------------------------------------------------------------
    PARAMETER(CP=1.0046E+3,RD=287.04,RV=4.6150E+2 &
    ,         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)
    PARAMETER(PSATK=PSAT*1.E-3)
    PARAMETER(DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP))
    PARAMETER(DLDTI=CVAP-CICE,XAI=-DLDTI/RV,XBI=XAI+HSUB/(RV*TTP))
!-----------------------------------------------------------------------
    TR=TTP/T

    IF(T >= TTP)THEN
        FPVSX=PSATK*(TR**XA)*EXP(XB*(1.-TR))
    ELSE
        FPVSX=PSATK*(TR**XAI)*EXP(XBI*(1.-TR))
    ENDIF

    RETURN
    END FUNCTION FPVSX
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
    FUNCTION FPVSX0(T)
!-----------------------------------------------------------------------
    PARAMETER(CP=1.0046E+3,RD=287.04,RV=4.6150E+2 &
    ,         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)
    PARAMETER(PSATK=PSAT*1.E-3)
    PARAMETER(DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP))
    PARAMETER(DLDTI=CVAP-CICE,XAI=-DLDT/RV,XBI=XA+HSUB/(RV*TTP))
!-----------------------------------------------------------------------
    TR=TTP/T
    FPVSX0=PSATK*(TR**XA)*EXP(XB*(1.-TR))

    RETURN
    END FUNCTION FPVSX0
