!
!#######################################################################
!------- 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"
      INTEGER, PARAMETER :: LP1=LM+1
      INCLUDE "MAPOT.comm"
!
!-------------------------------------------------------------------------
!-------------- 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
!
!--- 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
      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    !IBM
      CBFR=20.*PI*PI*BBFR*RHOL*1.E-33     !OTHERS
!
!--- 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
!
!#######################################################################
!--- 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
!
!#######################################################################
!-- Lookup tables for the saturation vapor pressure w/r/t water & ice --
!#######################################################################
!
      SUBROUTINE GPVS
C ******************************************************************
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . .
C SUBPROGRAM: GPVS COMPUTE SATURATION VAPOR PRESSURE TABLE
C AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
C
C ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE TABLE AS A FUNCTION OF
C TEMPERATURE FOR THE TABLE LOOKUP FUNCTION FPVS.
C EXACT SATURATION VAPOR PRESSURES ARE CALCULATED IN SUBPROGRAM FPVSX.
C THE CURRENT IMPLEMENTATION COMPUTES A TABLE WITH A LENGTH
C OF 7501 FOR TEMPERATURES RANGING FROM 180.0 TO 330.0 KELVIN.
C
C PROGRAM HISTORY LOG:
C 91-05-07 IREDELL
C 94-12-30 IREDELL EXPAND TABLE
C 96-02-19 HONG ICE EFFECT
C
C USAGE: CALL GPVS
C
C SUBPROGRAMS CALLED:
C (FPVSX) - INLINABLE FUNCTION TO COMPUTE SATURATION VAPOR PRESSURE
C
C COMMON BLOCKS:
C COMPVS - SCALING PARAMETERS AND TABLE FOR FUNCTION FPVS.
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 90
C MACHINE: IBM SP
C
C$$$
C----------------------------------------------------------------------
      PARAMETER(NX=7501)
      REAL TBPVS(NX),TBPVS0(NX)
      COMMON/COMPVS0/ C1XPVS0,C2XPVS0,TBPVS0
      COMMON/COMPVS/ C1XPVS,C2XPVS,TBPVS
C----------------------------------------------------------------------
      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
C
      DO JX=1,NX
        X=XMIN+(JX-1)*XINC
        T=X
        TBPVS(JX)=FPVSX(T)
        TBPVS0(JX)=FPVSX0(T)
      ENDDO
C
      RETURN
      END
C-----------------------------------------------------------------------
C***********************************************************************
C-----------------------------------------------------------------------
                           FUNCTION FPVS(T)
C-----------------------------------------------------------------------
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . .
C SUBPROGRAM: FPVS COMPUTE SATURATION VAPOR PRESSURE
C AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
C
C ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE FROM THE TEMPERATURE.
C A LINEAR INTERPOLATION IS DONE BETWEEN VALUES IN A LOOKUP TABLE
C COMPUTED IN GPVS. SEE DOCUMENTATION FOR FPVSX FOR DETAILS.
C INPUT VALUES OUTSIDE TABLE RANGE ARE RESET TO TABLE EXTREMA.
C THE INTERPOLATION ACCURACY IS ALMOST 6 DECIMAL PLACES.
C ON THE CRAY, FPVS IS ABOUT 4 TIMES FASTER THAN EXACT CALCULATION.
C THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
C
C PROGRAM HISTORY LOG:
C 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
C 94-12-30 IREDELL EXPAND TABLE
C 96-02-19 HONG ICE EFFECT
C
C USAGE: PVS=FPVS(T)
C
C INPUT ARGUMENT LIST:
C T - REAL TEMPERATURE IN KELVIN
C
C OUTPUT ARGUMENT LIST:
C FPVS - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
C
C COMMON BLOCKS:
C COMPVS - SCALING PARAMETERS AND TABLE COMPUTED IN GPVS.
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 90
C MACHINE: IBM SP
C
C$$$
C-----------------------------------------------------------------------
      PARAMETER(NX=7501)
      REAL TBPVS(NX)
      COMMON/COMPVS/ C1XPVS,C2XPVS,TBPVS
C-----------------------------------------------------------------------
      XJ=MIN(MAX(C1XPVS+C2XPVS*T,1.),FLOAT(NX))
      JX=MIN(XJ,NX-1.)
      FPVS=TBPVS(JX)+(XJ-JX)*(TBPVS(JX+1)-TBPVS(JX))
C
      RETURN
      END
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
                         FUNCTION FPVS0(T)
C-----------------------------------------------------------------------
      PARAMETER(NX=7501)
      REAL TBPVS0(NX)
      COMMON/COMPVS0/ C1XPVS0,C2XPVS0,TBPVS0
C-----------------------------------------------------------------------
      XJ1=MIN(MAX(C1XPVS0+C2XPVS0*T,1.),FLOAT(NX))
      JX1=MIN(XJ1,NX-1.)
      FPVS0=TBPVS0(JX1)+(XJ1-JX1)*(TBPVS0(JX1+1)-TBPVS0(JX1))
C
      RETURN
      END
C-----------------------------------------------------------------------
C***********************************************************************
C-----------------------------------------------------------------------
                         FUNCTION FPVSX(T)
C-----------------------------------------------------------------------
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . .
C SUBPROGRAM: FPVSX COMPUTE SATURATION VAPOR PRESSURE
C AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
C
C ABSTRACT: EXACTLY COMPUTE SATURATION VAPOR PRESSURE FROM TEMPERATURE.
C THE WATER MODEL ASSUMES A PERFECT GAS, CONSTANT SPECIFIC HEATS
C FOR GAS AND LIQUID, AND NEGLECTS THE VOLUME OF THE LIQUID.
C THE MODEL DOES ACCOUNT FOR THE VARIATION OF THE LATENT HEAT
C OF CONDENSATION WITH TEMPERATURE. THE ICE OPTION IS NOT INCLUDED.
C THE CLAUSIUS-CLAPEYRON EQUATION IS INTEGRATED FROM THE TRIPLE POINT
C TO GET THE FORMULA
C PVS=PSATK*(TR**XA)*EXP(XB*(1.-TR))
C WHERE TR IS TTP/T AND OTHER VALUES ARE PHYSICAL CONSTANTS
C THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
C
C PROGRAM HISTORY LOG:
C 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
C 94-12-30 IREDELL EXACT COMPUTATION
C 96-02-19 HONG ICE EFFECT
C
C USAGE: PVS=FPVSX(T)
C REFERENCE: EMANUEL(1994),116-117
C
C INPUT ARGUMENT LIST:
C T - REAL TEMPERATURE IN KELVIN
C
C OUTPUT ARGUMENT LIST:
C FPVSX - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 90
C MACHINE: IBM SP
C
C$$$
C-----------------------------------------------------------------------
      PARAMETER(CP=1.0046E+3,RD=287.04,RV=4.6150E+2
     1, TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2
     2, CLIQ=4.1855E+3,CVAP= 1.8460E+3
     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))
C-----------------------------------------------------------------------
      TR=TTP/T
C
      IF(T.GE.TTP)THEN
        FPVSX=PSATK*(TR**XA)*EXP(XB*(1.-TR))
      ELSE
        FPVSX=PSATK*(TR**XAI)*EXP(XBI*(1.-TR))
      ENDIF
C
      RETURN
      END
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
                        FUNCTION FPVSX0(T)
C-----------------------------------------------------------------------
      PARAMETER(CP=1.0046E+3,RD=287.04,RV=4.6150E+2
     1, TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2
     2, CLIQ=4.1855E+3,CVAP=1.8460E+3
     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))
C-----------------------------------------------------------------------
      TR=TTP/T
      FPVSX0=PSATK*(TR**XA)*EXP(XB*(1.-TR))
C
      RETURN
      END