!  Program Name:
!  Author(s)/Contact(s):
!  Abstract:
!  History Log:
!
!  Usage:
!  Parameters: <Specify typical arguments passed>
!  Input Files:
!        <list file names and briefly describe the data they include>
!  Output Files:
!        <list file names and briefly describe the information they include>
!
!  Condition codes:
!        <list exit condition or error codes returned >
!        If appropriate, descriptive troubleshooting instructions or
!        likely causes for failures could be mentioned here with the
!        appropriate error code
!
!  User controllable options: <if applicable>

MODULE module_sf_noahmpdrv

!-------------------------------
#if ( WRF_CHEM == 1 )
  USE module_data_gocart_dust
#endif
!-------------------------------

!
CONTAINS
!
  SUBROUTINE noahmplsm(ITIMESTEP,        YR,   JULIAN,   COSZIN,   XLATIN,  & ! IN : Time/Space-related
                  DZ8W,       DT,       DZS,    NSOIL,       DX,            & ! IN : Model configuration
	        IVGTYP,   ISLTYP,    VEGFRA,   VEGMAX,      TMN,            & ! IN : Vegetation/Soil characteristics
		 XLAND,     XICE,XICE_THRES,                                & ! IN : Vegetation/Soil characteristics
                 IDVEG, IOPT_CRS,  IOPT_BTR, IOPT_RUN, IOPT_SFC, IOPT_FRZ,  & ! IN : User options
              IOPT_INF, IOPT_RAD,  IOPT_ALB, IOPT_SNF,IOPT_TBOT, IOPT_STC,  & ! IN : User options
              IOPT_GLA, IOPT_RSF, IOPT_SOIL,IOPT_PEDO,IOPT_CROP,            & ! IN : User options
              IOPT_IMPERV, IZ0TLND,                                         & ! IN : User options
                   T3D,     QV3D,     U_PHY,    V_PHY,   SWDOWN,      GLW,  & ! IN : Forcing
		 P8W3D,PRECIP_IN,        SR,                                & ! IN : Forcing
                   TSK,      HFX,      QFX,        LH,   GRDFLX,    SMSTAV, & ! IN/OUT LSM eqv
                SMSTOT,SFCRUNOFF, UDRUNOFF,    ALBEDO,    SNOWC,     SMOIS, & ! IN/OUT LSM eqv
		  SH2O,     TSLB,     SNOW,     SNOWH,   CANWAT,    ACSNOM, & ! IN/OUT LSM eqv
		ACSNOW,    EMISS,     QSFC,                                 & ! IN/OUT LSM eqv
 		    Z0,      ZNT,                                           & ! IN/OUT LSM eqv
               ISNOWXY,     TVXY,     TGXY,  CANICEXY, CANLIQXY,     EAHXY, & ! IN/OUT Noah MP only
	         TAHXY,     CMXY,     CHXY,    FWETXY, SNEQVOXY,  ALBOLDXY, & ! IN/OUT Noah MP only
               QSNOWXY, QRAINXY, WSLAKEXY, ZWTXY,  WAXY,   WTXY,    TSNOXY, & ! IN/OUT Noah MP only
	       ZSNSOXY,  SNICEXY,  SNLIQXY,  LFMASSXY, RTMASSXY,  STMASSXY, & ! IN/OUT Noah MP only
	        WOODXY, STBLCPXY, FASTCPXY,    XLAIXY,   XSAIXY,   TAUSSXY, & ! IN/OUT Noah MP only
	       SMOISEQ, SMCWTDXY,DEEPRECHXY,   RECHXY,                      & ! IN/OUT Noah MP only
	        T2MVXY,   T2MBXY,    Q2MVXY,   Q2MBXY,                      & ! OUT Noah MP only
	        TRADXY,    NEEXY,    GPPXY,     NPPXY,   FVEGXY,   RUNSFXY, & ! OUT Noah MP only
	       RUNSBXY,   ECANXY,   EDIRXY,   ETRANXY,    FSAXY,    FIRAXY, & ! OUT Noah MP only
	        APARXY,    PSNXY,    SAVXY,     SAGXY,  RSSUNXY,   RSSHAXY, & ! OUT Noah MP only
              ALBSNDXY, ALBSNIXY,                                           & ! OUT Noah MP only
		BGAPXY,   WGAPXY,    TGVXY,     TGBXY,    CHVXY,     CHBXY, & ! OUT Noah MP only
		 SHGXY,    SHCXY,    SHBXY,     EVGXY,    EVBXY,     GHVXY, & ! OUT Noah MP only
		 GHBXY,    IRGXY,    IRCXY,     IRBXY,     TRXY,     EVCXY, & ! OUT Noah MP only
              CHLEAFXY,   CHUCXY,   CHV2XY,    CHB2XY,                      & ! OUT Noah MP only
              PSNOWLIQXY, PSNOWTEMPXY, PSNOWDZXY,                           & ! OUT    for crocus
              PSNOWHEATXY, PSNOWRHOXY, PSNOWSWEXY,                          & ! IN/OUT for crocus
              PSNOWGRAN1XY, PSNOWGRAN2XY, PSNOWHISTXY,                      & ! IN/OUT for crocus
              PSNOWALBXY, PSNOWAGEXY, act_level,                            & ! IN/OUT for crocus
              PSNOWHEIGHTXY, PSNOWTOTSWEXY,  PSNOWTHRUFALXY,                & ! OUT    for crocus
              GLACINFO ,GLACT,                                              & ! IN/OUT for crocus
              FLOW_ICE, FLOW_SNOW,                                          & ! OUT    for crocus
              crocus_opt,                                                   & ! IN     for crocus
#ifdef SPATIAL_SOIL
                 BEXP_3D,SMCDRY_3D,SMCWLT_3D,SMCREF_3D,SMCMAX_3D,             &
		 DKSAT_3D,DWSAT_3D,PSISAT_3D,QUARTZ_3D,                       &
		 REFDK_2D,REFKDT_2D,SLOPE_2D,                                 &
		 CWPVT_2D,VCMX25_2D,MP_2D,HVT_2D,MFSNO_2D,RSURFEXP_2D,        &
                 AXAJ_2D,BXAJ_2D,XXAJ_2D,                                     &
                 IMPERV_2D,                                                   &
                 SSI_2D,SNOWRETFAC_2D,TAU0_2D,RSURFSNOW_2D,SCAMAX_2D,         &
#endif
#ifdef WRF_HYDRO
               sfcheadrt,INFXSRT,soldrain,                                  &
#endif
               ids,ide,  jds,jde,  kds,kde,                    &
               ims,ime,  jms,jme,  kms,kme,                    &
               its,ite,  jts,jte,  kts,kte,                    &
               MP_RAINC, MP_RAINNC, MP_SHCV, MP_SNOW, MP_GRAUP, MP_HAIL, &
               VIS_ICEALB  &
#ifdef WRF_HYDRO
               , ACCPRCP,  ACCECAN, ACCETRAN,   ACCEDIR        &         ! NEW output accumulator variables
               , SOILSAT_TOP, SOILSAT, SOILICE, SNOWT_AVG      &         ! NEW soil saturation and snow temp
#endif
               )
!----------------------------------------------------------------
    USE MODULE_SF_NOAHMPLSM
    USE module_sf_noahmp_glacier
    USE NOAHMP_TABLES, ONLY: ISICE_TABLE, CO2_TABLE, O2_TABLE, HVB_TABLE, HVT_TABLE
    USE module_snowcro
    USE MODD_SNOW_PAR, ONLY: XZ0SN, XZ0HSN

!----------------------------------------------------------------
    IMPLICIT NONE
!----------------------------------------------------------------

! IN only

    INTEGER,                                         INTENT(IN   ) ::  ITIMESTEP ! timestep number
    INTEGER,                                         INTENT(IN   ) ::  YR        ! 4-digit year
    REAL,                                            INTENT(IN   ) ::  JULIAN    ! Julian day
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) ::  COSZIN    ! cosine zenith angle
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) ::  XLATIN    ! latitude [rad]
    REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) ::  DZ8W      ! thickness of atmo layers [m]
    REAL,                                            INTENT(IN   ) ::  DT        ! timestep [s]
    REAL,    DIMENSION(1:nsoil),                     INTENT(IN   ) ::  DZS       ! thickness of soil layers [m]
    INTEGER,                                         INTENT(IN   ) ::  NSOIL     ! number of soil layers
    REAL,                                            INTENT(IN   ) ::  DX        ! horizontal grid spacing [m]
    INTEGER, DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) ::  IVGTYP    ! vegetation type
    INTEGER, DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) ::  ISLTYP    ! soil type
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) ::  VEGFRA    ! vegetation fraction []
    REAL,    DIMENSION( ims:ime ,         jms:jme ), INTENT(IN   ) ::  VEGMAX    ! annual max vegetation fraction []
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) ::  TMN       ! deep soil temperature [K]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) ::  XLAND     ! =2 ocean; =1 land/seaice
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) ::  XICE      ! fraction of grid that is seaice
    REAL,                                            INTENT(IN   ) ::  XICE_THRES! fraction of grid determining seaice
    INTEGER,                                         INTENT(IN   ) ::  IDVEG     ! dynamic vegetation (1 -> off ; 2 -> on) with opt_crs = 1
    INTEGER,                                         INTENT(IN   ) ::  IOPT_CRS  ! canopy stomatal resistance (1-> Ball-Berry; 2->Jarvis)
    INTEGER,                                         INTENT(IN   ) ::  IOPT_BTR  ! soil moisture factor for stomatal resistance (1-> Noah; 2-> CLM; 3-> SSiB)
    INTEGER,                                         INTENT(IN   ) ::  IOPT_RUN  ! runoff and groundwater (1->SIMGM; 2->SIMTOP; 3->Schaake96; 4->BATS)
    INTEGER,                                         INTENT(IN   ) ::  IOPT_SFC  ! surface layer drag coeff (CH & CM) (1->M-O; 2->Chen97)
    INTEGER,                                         INTENT(IN   ) ::  IOPT_FRZ  ! supercooled liquid water (1-> NY06; 2->Koren99)
    INTEGER,                                         INTENT(IN   ) ::  IOPT_INF  ! frozen soil permeability (1-> NY06; 2->Koren99)
    INTEGER,                                         INTENT(IN   ) ::  IOPT_RAD  ! radiation transfer (1->gap=F(3D,cosz); 2->gap=0; 3->gap=1-Fveg)
    INTEGER,                                         INTENT(IN   ) ::  IOPT_ALB  ! snow surface albedo (1->BATS; 2->CLASS)
    INTEGER,                                         INTENT(IN   ) ::  IOPT_SNF  ! rainfall & snowfall (1-Jordan91; 2->BATS; 3->Noah)
    INTEGER,                                         INTENT(IN   ) ::  IOPT_TBOT ! lower boundary of soil temperature (1->zero-flux; 2->Noah)
    INTEGER,                                         INTENT(IN   ) ::  IOPT_STC  ! snow/soil temperature time scheme
    INTEGER,                                         INTENT(IN   ) ::  IOPT_GLA  ! glacier option (1->phase change; 2->simple)
    INTEGER,                                         INTENT(IN   ) ::  IOPT_RSF  ! surface resistance (1->Sakaguchi/Zeng; 2->Seller; 3->mod Sellers; 4->1+snow)
    INTEGER,                                         INTENT(IN   ) ::  IOPT_SOIL ! soil configuration option
    INTEGER,                                         INTENT(IN   ) ::  IOPT_PEDO ! soil pedotransfer function option
    INTEGER,                                         INTENT(IN   ) ::  IOPT_CROP ! crop model option (0->none; 1->Liu et al.)
    INTEGER,                                         INTENT(IN   ) ::  IOPT_IMPERV !imperviousness infiltration adjustment (0->none; 1->total;
                                                                                        !9->old)
    INTEGER,                                         INTENT(IN   ) ::  IZ0TLND   ! option of Chen adjustment of Czil (not used)
    REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) ::  T3D       ! 3D atmospheric temperature valid at mid-levels [K]
    REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) ::  QV3D      ! 3D water vapor mixing ratio [kg/kg_dry]
    REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) ::  U_PHY     ! 3D U wind component [m/s]
    REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) ::  V_PHY     ! 3D V wind component [m/s]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) ::  SWDOWN    ! solar down at surface [W m-2]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) ::  GLW       ! longwave down at surface [W m-2]
    REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) ::  P8W3D     ! 3D pressure, valid at interface [Pa]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) ::  PRECIP_IN ! total input precipitation [mm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) ::  SR        ! frozen precipitation ratio [-]

!Optional Detailed Precipitation Partitioning Inputs
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ), OPTIONAL ::  MP_RAINC  ! convective precipitation entering land model [mm] ! MB/AN : v3.7
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ), OPTIONAL ::  MP_RAINNC ! large-scale precipitation entering land model [mm]! MB/AN : v3.7
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ), OPTIONAL ::  MP_SHCV   ! shallow conv precip entering land model [mm]      ! MB/AN : v3.7
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ), OPTIONAL ::  MP_SNOW   ! snow precipitation entering land model [mm]       ! MB/AN : v3.7
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ), OPTIONAL ::  MP_GRAUP  ! graupel precipitation entering land model [mm]    ! MB/AN : v3.7
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ), OPTIONAL ::  MP_HAIL   ! hail precipitation entering land model [mm]       ! MB/AN : v3.7
#ifdef WRF_HYDRO
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  sfcheadrt,INFXSRT,soldrain   ! for WRF-Hydro
#endif
#ifdef SPATIAL_SOIL
    REAL,    DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(IN) ::  BEXP_3D   ! C-H B exponent
    REAL,    DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(IN) ::  SMCDRY_3D ! Soil Moisture Limit: Dry
    REAL,    DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(IN) ::  SMCWLT_3D ! Soil Moisture Limit: Wilt
    REAL,    DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(IN) ::  SMCREF_3D ! Soil Moisture Limit: Reference
    REAL,    DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(IN) ::  SMCMAX_3D ! Soil Moisture Limit: Max
    REAL,    DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(IN) ::  DKSAT_3D  ! Saturated Soil Conductivity
    REAL,    DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(IN) ::  DWSAT_3D  ! Saturated Soil Diffusivity
    REAL,    DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(IN) ::  PSISAT_3D ! Saturated Matric Potential
    REAL,    DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(IN) ::  QUARTZ_3D ! Soil quartz content
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  REFDK_2D  ! Reference Soil Conductivity
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  REFKDT_2D ! Soil Infiltration Parameter
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  SLOPE_2D  ! Soil Drainage Parameter
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  CWPVT_2D  ! Canopy wind parameter
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  VCMX25_2D ! VCmax at 25C
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  MP_2D     ! Slope of Ball-Berry rs-P relationship
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  HVT_2D    ! Canopy Height
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  MFSNO_2D  ! Snow cover m parameter
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  RSURFEXP_2D ! exponent in the shape parameter for soil resistance option 1
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  AXAJ_2D   ! Xinanjiang: Tension water distribution inflection parameter [-]
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  BXAJ_2D   ! Xinanjiang: Tension water distribution shape parameter [-]
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  XXAJ_2D   ! Xinanjiang: Free water distribution shape parameter [-]
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  IMPERV_2D ! impervious fraction
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  SSI_2D         ! liquid water holding capacity for snowpack (m3/m3)
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  SNOWRETFAC_2D  ! snowpack water release timescale factor (1/s)
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  TAU0_2D        ! tau0 from Yang97 eqn. 10a
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  RSURFSNOW_2D   ! surface resistence for snow [s/m]
    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN)          ::  SCAMAX_2D      ! max fractional snow covered area (0.0-1.0)
#endif

! INOUT (with generic LSM equivalent)

    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  TSK       ! surface radiative temperature [K]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  HFX       ! sensible heat flux [W m-2]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  QFX       ! latent heat flux [kg s-1 m-2]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  LH        ! latent heat flux [W m-2]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  GRDFLX    ! ground/snow heat flux [W m-2]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  SMSTAV    ! soil moisture avail. [not used]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  SMSTOT    ! total soil water [mm][not used]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  SFCRUNOFF ! accumulated surface runoff [m]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  UDRUNOFF  ! accumulated sub-surface runoff [m]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  ALBEDO    ! total grid albedo []
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  SNOWC     ! snow cover fraction []
    REAL,    DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(INOUT) ::  SMOIS     ! volumetric soil moisture [m3/m3]
    REAL,    DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(INOUT) ::  SH2O      ! volumetric liquid soil moisture [m3/m3]
    REAL,    DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(INOUT) ::  TSLB      ! soil temperature [K]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  SNOW      ! snow water equivalent [mm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  SNOWH     ! physical snow depth [m]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  CANWAT    ! total canopy water + ice [mm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  ACSNOM    ! accumulated snow melt leaving pack
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  ACSNOW    ! accumulated snow on grid
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  EMISS     ! surface bulk emissivity
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  QSFC      ! bulk surface specific humidity
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  Z0        ! combined z0 sent to coupled model
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  ZNT       ! combined z0 sent to coupled model

! INOUT (with no Noah LSM equivalent)

    INTEGER, DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  ISNOWXY   ! actual no. of snow layers
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  TVXY      ! vegetation leaf temperature
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  TGXY      ! bulk ground surface temperature
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  CANICEXY  ! canopy-intercepted ice (mm)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  CANLIQXY  ! canopy-intercepted liquid water (mm)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  EAHXY     ! canopy air vapor pressure (pa)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  TAHXY     ! canopy air temperature (k)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  CMXY      ! bulk momentum drag coefficient
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  CHXY      ! bulk sensible heat exchange coefficient
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  FWETXY    ! wetted or snowed fraction of the canopy (-)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  SNEQVOXY  ! snow mass at last time step(mm h2o)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  ALBOLDXY  ! snow albedo at last time step (-)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  QSNOWXY   ! snowfall on the ground [mm/s]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  QRAINXY   ! rainfall on the ground [mm/s]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  WSLAKEXY  ! lake water storage [mm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  ZWTXY     ! water table depth [m]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  WAXY      ! water in the "aquifer" [mm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  WTXY      ! groundwater storage [mm]
    REAL,    DIMENSION( ims:ime,-2:0,     jms:jme ), INTENT(INOUT) ::  TSNOXY    ! snow temperature [K]
    REAL,    DIMENSION( ims:ime,-2:NSOIL, jms:jme ), INTENT(INOUT) ::  ZSNSOXY   ! snow layer depth [m]
    REAL,    DIMENSION( ims:ime,-2:0,     jms:jme ), INTENT(INOUT) ::  SNICEXY   ! snow layer ice [mm]
    REAL,    DIMENSION( ims:ime,-2:0,     jms:jme ), INTENT(INOUT) ::  SNLIQXY   ! snow layer liquid water [mm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  LFMASSXY  ! leaf mass [g/m2]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  RTMASSXY  ! mass of fine roots [g/m2]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  STMASSXY  ! stem mass [g/m2]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  WOODXY    ! mass of wood (incl. woody roots) [g/m2]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  STBLCPXY  ! stable carbon in deep soil [g/m2]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  FASTCPXY  ! short-lived carbon, shallow soil [g/m2]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  XLAIXY    ! leaf area index
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  XSAIXY    ! stem area index
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  TAUSSXY   ! snow age factor
    REAL,    DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(INOUT) ::  SMOISEQ   ! eq volumetric soil moisture [m3/m3]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  SMCWTDXY  ! soil moisture content in the layer to the water table when deep
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  DEEPRECHXY ! recharge to the water table when deep
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(INOUT) ::  RECHXY    ! recharge to the water table (diagnostic)

! OUT (with no Noah LSM equivalent)

    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  T2MVXY    ! 2m temperature of vegetation part
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  T2MBXY    ! 2m temperature of bare ground part
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  Q2MVXY    ! 2m mixing ratio of vegetation part
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  Q2MBXY    ! 2m mixing ratio of bare ground part
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  TRADXY    ! surface radiative temperature (k)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  NEEXY     ! net ecosys exchange (g/m2/s CO2)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  GPPXY     ! gross primary assimilation [g/m2/s C]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  NPPXY     ! net primary productivity [g/m2/s C]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  FVEGXY    ! Noah-MP vegetation fraction [-]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  RUNSFXY   ! surface runoff [mm/s]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  RUNSBXY   ! subsurface runoff [mm/s]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  ECANXY    ! evaporation of intercepted water (mm/s)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  EDIRXY    ! soil surface evaporation rate (mm/s]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  ETRANXY   ! transpiration rate (mm/s)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  FSAXY     ! total absorbed solar radiation (w/m2)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  FIRAXY    ! total net longwave rad (w/m2) [+ to atm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  APARXY    ! photosyn active energy by canopy (w/m2)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  PSNXY     ! total photosynthesis (umol co2/m2/s) [+]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  SAVXY     ! solar rad absorbed by veg. (w/m2)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  SAGXY     ! solar rad absorbed by ground (w/m2)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  RSSUNXY   ! sunlit leaf stomatal resistance (s/m)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  RSSHAXY   ! shaded leaf stomatal resistance (s/m)
    REAL,    DIMENSION( ims:ime,1:2,      jms:jme ), INTENT(OUT  ) ::  ALBSNDXY  ! snow albedo (direct)
    REAL,    DIMENSION( ims:ime,1:2,      jms:jme ), INTENT(OUT  ) ::  ALBSNIXY  ! snow albedo (diffuse)
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  BGAPXY    ! between gap fraction
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  WGAPXY    ! within gap fraction
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  TGVXY     ! under canopy ground temperature[K]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  TGBXY     ! bare ground temperature [K]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  CHVXY     ! sensible heat exchange coefficient vegetated
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  CHBXY     ! sensible heat exchange coefficient bare-ground
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  SHGXY     ! veg ground sen. heat [w/m2]   [+ to atm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  SHCXY     ! canopy sen. heat [w/m2]   [+ to atm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  SHBXY     ! bare sensible heat [w/m2]     [+ to atm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  EVGXY     ! veg ground evap. heat [w/m2]  [+ to atm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  EVBXY     ! bare soil evaporation [w/m2]  [+ to atm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  GHVXY     ! veg ground heat flux [w/m2]  [+ to soil]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  GHBXY     ! bare ground heat flux [w/m2] [+ to soil]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  IRGXY     ! veg ground net LW rad. [w/m2] [+ to atm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  IRCXY     ! canopy net LW rad. [w/m2] [+ to atm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  IRBXY     ! bare net longwave rad. [w/m2] [+ to atm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  TRXY      ! transpiration [w/m2]  [+ to atm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  EVCXY     ! canopy evaporation heat [w/m2]  [+ to atm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  CHLEAFXY  ! leaf exchange coefficient
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  CHUCXY    ! under canopy exchange coefficient
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  CHV2XY    ! veg 2m exchange coefficient
    REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(OUT  ) ::  CHB2XY    ! bare 2m exchange coefficient
    INTEGER,  INTENT(IN   )   ::     ids,ide, jds,jde, kds,kde,  &  ! d -> domain
         &                           ims,ime, jms,jme, kms,kme,  &  ! m -> memory
         &                           its,ite, jts,jte, kts,kte      ! t -> tile

#ifdef WRF_HYDRO
    REAL,    DIMENSION( ims:ime,          jms:jme ), OPTIONAL, INTENT(INOUT) ::  ACCPRCP     ! accumulated precip [mm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), OPTIONAL, INTENT(INOUT) ::  ACCECAN     ! accumulated canopy evap [mm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), OPTIONAL, INTENT(INOUT) ::  ACCETRAN    ! accumulated transpiration [mm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), OPTIONAL, INTENT(INOUT) ::  ACCEDIR     ! accumulated direct soil evap [mm]
    REAL,    DIMENSION( ims:ime,          jms:jme ), OPTIONAL, INTENT(INOUT) ::  SOILSAT_TOP ! soil saturation in the top layer [fraction]
    REAL,    DIMENSION( ims:ime,          jms:jme ), OPTIONAL, INTENT(INOUT) ::  SOILSAT     ! soil saturation column integrated [fraction]
    REAL,    DIMENSION( ims:ime,          jms:jme ), OPTIONAL, INTENT(INOUT) ::  SOILICE     ! fraction of soil moisture that is ice [fraction]
    REAL,    DIMENSION( ims:ime,          jms:jme ), OPTIONAL, INTENT(INOUT) ::  SNOWT_AVG   ! average snowpack temperature (by layer mass) [K]
#endif

! 1D equivalent of 2D/3D fields

! IN only

    REAL                                :: COSZ         ! cosine zenith angle
    REAL                                :: LAT          ! latitude [rad]
    REAL                                :: Z_ML         ! model height [m]
    INTEGER                             :: VEGTYP       ! vegetation type
    INTEGER                             :: SOILTYP      ! soil type
    INTEGER                             :: CROPTYPE     ! crop type
    REAL                                :: FVEG         ! vegetation fraction [-]
    REAL                                :: FVGMAX       ! annual max vegetation fraction []
    REAL                                :: TBOT         ! deep soil temperature [K]
    REAL                                :: T_ML         ! temperature valid at mid-levels [K]
    REAL                                :: Q_ML         ! water vapor mixing ratio [kg/kg_dry]
    REAL                                :: U_ML         ! U wind component [m/s]
    REAL                                :: V_ML         ! V wind component [m/s]
    REAL                                :: SWDN         ! solar down at surface [W m-2]
    REAL                                :: LWDN         ! longwave down at surface [W m-2]
    REAL                                :: P_ML         ! pressure, valid at interface [Pa]
    REAL                                :: PSFC         ! surface pressure [Pa]
    REAL                                :: PRCP         ! total precipitation entering  [mm]         ! MB/AN : v3.7
    REAL                                :: PRCPCONV     ! convective precipitation entering  [mm]    ! MB/AN : v3.7
    REAL                                :: PRCPNONC     ! non-convective precipitation entering [mm] ! MB/AN : v3.7
    REAL                                :: PRCPSHCV     ! shallow convective precip entering  [mm]   ! MB/AN : v3.7
    REAL                                :: PRCPSNOW     ! snow entering land model [mm]              ! MB/AN : v3.7
    REAL                                :: PRCPGRPL     ! graupel entering land model [mm]           ! MB/AN : v3.7
    REAL                                :: PRCPHAIL     ! hail entering land model [mm]              ! MB/AN : v3.7
    REAL                                :: PRCPOTHR     ! other precip, e.g. fog [mm]                ! MB/AN : v3.7

! INOUT (with generic LSM equivalent)

    REAL                                :: FSH          ! total sensible heat (w/m2) [+ to atm]
    REAL                                :: SSOIL        ! soil heat heat (w/m2)
    REAL                                :: SALB         ! surface albedo (-)
    REAL                                :: FSNO         ! snow cover fraction (-)
    REAL,   DIMENSION( 1:NSOIL)         :: SMCEQ        ! eq vol. soil moisture (m3/m3)
    REAL,   DIMENSION( 1:NSOIL)         :: SMC          ! vol. soil moisture (m3/m3)
    REAL,   DIMENSION( 1:NSOIL)         :: SMH2O        ! vol. soil liquid water (m3/m3)
    REAL,   DIMENSION(-2:NSOIL)         :: STC          ! snow/soil tmperatures
    REAL                                :: SWE          ! snow water equivalent (mm)
    REAL                                :: SNDPTH       ! snow depth (m)
    REAL                                :: EMISSI       ! net surface emissivity
    REAL                                :: QSFC1D       ! bulk surface specific humidity

! INOUT (with no Noah LSM equivalent)

    INTEGER                             :: ISNOW        ! actual no. of snow layers
    REAL                                :: TV           ! vegetation canopy temperature
    REAL                                :: TG           ! ground surface temperature
    REAL                                :: CANICE       ! canopy-intercepted ice (mm)
    REAL                                :: CANLIQ       ! canopy-intercepted liquid water (mm)
    REAL                                :: EAH          ! canopy air vapor pressure (pa)
    REAL                                :: TAH          ! canopy air temperature (k)
    REAL                                :: CM           ! momentum drag coefficient
    REAL                                :: CH           ! sensible heat exchange coefficient
    REAL                                :: FWET         ! wetted or snowed fraction of the canopy (-)
    REAL                                :: SNEQVO       ! snow mass at last time step(mm h2o)
    REAL                                :: ALBOLD       ! snow albedo at last time step (-)
    REAL                                :: QSNOW        ! snowfall on the ground [mm/s]
    REAL                                :: QRAIN        ! rainfall on the ground [mm/s]
    REAL                                :: WSLAKE       ! lake water storage [mm]
    REAL                                :: ZWT          ! water table depth [m]
    REAL                                :: WA           ! water in the "aquifer" [mm]
    REAL                                :: WT           ! groundwater storage [mm]
    REAL                                :: SMCWTD       ! soil moisture content in the layer to the water table when deep
    REAL                                :: DEEPRECH     ! recharge to the water table when deep
    REAL                                :: RECH         ! recharge to the water table (diagnostic)
    REAL, DIMENSION(-2:NSOIL)           :: ZSNSO        ! snow layer depth [m]
    REAL, DIMENSION(-2:              0) :: SNICE        ! snow layer ice [mm]
    REAL, DIMENSION(-2:              0) :: SNLIQ        ! snow layer liquid water [mm]
    REAL                                :: LFMASS       ! leaf mass [g/m2]
    REAL                                :: RTMASS       ! mass of fine roots [g/m2]
    REAL                                :: STMASS       ! stem mass [g/m2]
    REAL                                :: WOOD         ! mass of wood (incl. woody roots) [g/m2]
    REAL                                :: GRAIN        ! mass of grain XING [g/m2]
    REAL                                :: GDD          ! mass of grain XING[g/m2]
    INTEGER                             :: PGS          !stem respiration [g/m2/s]
    REAL                                :: STBLCP       ! stable carbon in deep soil [g/m2]
    REAL                                :: FASTCP       ! short-lived carbon, shallow soil [g/m2]
    REAL                                :: PLAI         ! leaf area index
    REAL                                :: PSAI         ! stem area index
    REAL                                :: TAUSS        ! non-dimensional snow age

! OUT (with no Noah LSM equivalent)

    REAL                                :: Z0WRF        ! combined z0 sent to coupled model
    REAL                                :: T2MV         ! 2m temperature of vegetation part
    REAL                                :: T2MB         ! 2m temperature of bare ground part
    REAL                                :: Q2MV         ! 2m mixing ratio of vegetation part
    REAL                                :: Q2MB         ! 2m mixing ratio of bare ground part
    REAL                                :: TRAD         ! surface radiative temperature (k)
    REAL                                :: NEE          ! net ecosys exchange (g/m2/s CO2)
    REAL                                :: GPP          ! gross primary assimilation [g/m2/s C]
    REAL                                :: NPP          ! net primary productivity [g/m2/s C]
    REAL                                :: FVEGMP       ! greenness vegetation fraction [-]
    REAL                                :: RUNSF        ! surface runoff [mm/s]
    REAL                                :: RUNSB        ! subsurface runoff [mm/s]
    REAL                                :: ECAN         ! evaporation of intercepted water (mm/s)
    REAL                                :: ETRAN        ! transpiration rate (mm/s)
    REAL                                :: ESOIL        ! soil surface evaporation rate (mm/s]
    REAL                                :: FSA          ! total absorbed solar radiation (w/m2)
    REAL                                :: FIRA         ! total net longwave rad (w/m2) [+ to atm]
    REAL                                :: APAR         ! photosyn active energy by canopy (w/m2)
    REAL                                :: PSN          ! total photosynthesis (umol co2/m2/s) [+]
    REAL                                :: SAV          ! solar rad absorbed by veg. (w/m2)
    REAL                                :: SAG          ! solar rad absorbed by ground (w/m2)
    REAL                                :: RSSUN        ! sunlit leaf stomatal resistance (s/m)
    REAL                                :: RSSHA        ! shaded leaf stomatal resistance (s/m)
    REAL, DIMENSION(1:2)                :: ALBSND       ! snow albedo (direct)
    REAL, DIMENSION(1:2)                :: ALBSNI       ! snow albedo (diffuse)
    REAL                                :: RB           ! leaf boundary layer resistance (s/m)
    REAL                                :: LAISUN       ! sunlit leaf area index (m2/m2)
    REAL                                :: LAISHA       ! shaded leaf area index (m2/m2)
    REAL                                :: BGAP         ! between gap fraction
    REAL                                :: WGAP         ! within gap fraction
    REAL                                :: TGV          ! under canopy ground temperature[K]
    REAL                                :: TGB          ! bare ground temperature [K]
    REAL                                :: CHV          ! sensible heat exchange coefficient vegetated
    REAL                                :: CHB          ! sensible heat exchange coefficient bare-ground
    REAL                                :: IRC          ! canopy net LW rad. [w/m2] [+ to atm]
    REAL                                :: IRG          ! veg ground net LW rad. [w/m2] [+ to atm]
    REAL                                :: SHC          ! canopy sen. heat [w/m2]   [+ to atm]
    REAL                                :: SHG          ! veg ground sen. heat [w/m2]   [+ to atm]
    REAL                                :: EVG          ! veg ground evap. heat [w/m2]  [+ to atm]
    REAL                                :: GHV          ! veg ground heat flux [w/m2]  [+ to soil]
    REAL                                :: IRB          ! bare net longwave rad. [w/m2] [+ to atm]
    REAL                                :: SHB          ! bare sensible heat [w/m2]     [+ to atm]
    REAL                                :: EVB          ! bare evaporation heat [w/m2]  [+ to atm]
    REAL                                :: GHB          ! bare ground heat flux [w/m2] [+ to soil]
    REAL                                :: TR           ! transpiration [w/m2]  [+ to atm]
    REAL                                :: EVC          ! canopy evaporation heat [w/m2]  [+ to atm]
    REAL                                :: CHLEAF       ! leaf exchange coefficient
    REAL                                :: CHUC         ! under canopy exchange coefficient
    REAL                                :: CHV2         ! veg 2m exchange coefficient
    REAL                                :: CHB2         ! bare 2m exchange coefficient
    REAL                                :: PAHV         ! precipitation advected heat - vegetation net (W/m2)
    REAL                                :: PAHG         ! precipitation advected heat - under canopy net (W/m2)
    REAL                                :: PAHB         ! precipitation advected heat - bare ground net (W/m2)
    REAL                                :: PAH          ! precipitation advected heat - total (W/m2)

! Intermediate terms

    REAL                                :: FPICE        ! snow fraction of precip
    REAL                                :: FCEV         ! canopy evaporation heat (w/m2) [+ to atm]
    REAL                                :: FGEV         ! ground evaporation heat (w/m2) [+ to atm]
    REAL                                :: FCTR         ! transpiration heat flux (w/m2) [+ to atm]
    REAL                                :: QSNBOT       ! snowmelt out bottom of pack [mm/s]
    REAL                                :: PONDING      ! snowmelt with no pack [mm]
    REAL                                :: PONDING1     ! snowmelt with no pack [mm]
    REAL                                :: PONDING2     ! snowmelt with no pack [mm]

! Local terms

    REAL                                :: FSR          ! total reflected solar radiation (w/m2)
    REAL, DIMENSION(-2:0)               :: FICEOLD      ! snow layer ice fraction []
    REAL                                :: CO2PP        ! CO2 partial pressure [Pa]
    REAL                                :: O2PP         ! O2 partial pressure [Pa]
    REAL, DIMENSION(1:NSOIL)            :: ZSOIL        ! depth to soil interfaces [m]
    REAL                                :: FOLN         ! nitrogen saturation [%]

    REAL                                :: QC           ! cloud specific humidity for MYJ [not used]
    REAL                                :: PBLH         ! PBL height for MYJ [not used]
    REAL                                :: DZ8W1D       ! model level heights for MYJ [not used]

    INTEGER                             :: I
    INTEGER                             :: J
    INTEGER                             :: K
    INTEGER                             :: ICE
    INTEGER                             :: SLOPETYP
    LOGICAL                             :: IPRINT

    INTEGER                             :: SOILCOLOR    ! soil color index
    INTEGER                             :: IST          ! surface type 1-soil; 2-lake
    INTEGER                             :: YEARLEN

    INTEGER, PARAMETER                  :: NSNOW = 3    ! number of snow layers fixed to 3
    REAL, PARAMETER                     :: undefined_value = -1.E36

    type(noahmp_parameters) :: parameters


    INTEGER,                INTENT(IN)                 :: act_level
    REAL, DIMENSION(:,:),   INTENT(INOUT), allocatable :: PSNOWALBXY    ! Prognostic surface snow albedo
    REAL, DIMENSION(:,:,:), INTENT(INOUT), allocatable :: PSNOWHEATXY
    REAL, DIMENSION(:,:,:), INTENT(INOUT), allocatable :: PSNOWRHOXY
    REAL, DIMENSION(:,:,:), INTENT(INOUT), allocatable :: PSNOWSWEXY
    REAL, DIMENSION(:,:,:), INTENT(INOUT), allocatable :: PSNOWGRAN1XY   ! Snow layer(s) grain parameter 1
    REAL, DIMENSION(:,:,:), INTENT(INOUT), allocatable :: PSNOWGRAN2XY   ! Snow layer(s) grain parameter 2
    REAL, DIMENSION(:,:,:), INTENT(INOUT), allocatable :: PSNOWHISTXY       ! Snow layer(s) grain historical parameter
    REAL, DIMENSION(:,:,:), INTENT(INOUT), allocatable :: PSNOWAGEXY ! Snow grain age

    REAL, DIMENSION(:,:,:), INTENT(INOUT), allocatable   :: PSNOWLIQXY
    REAL, DIMENSION(:,:,:), INTENT(INOUT), allocatable   :: PSNOWTEMPXY
    REAL, DIMENSION(:,:,:), INTENT(INOUT), allocatable   :: PSNOWDZXY

    REAL, DIMENSION(:,:),   INTENT(INOUT), allocatable  :: PSNOWTHRUFALXY ! rate that liquid water leaves snow pack [kg/(m2 s)]   (should only be out)
    REAL, DIMENSION(:,:),   INTENT(INOUT), allocatable  :: PSNOWHEIGHTXY   ! integrated crocus snowheight PSNOWDZXY
    REAL, DIMENSION(:,:),   INTENT(INOUT), allocatable  :: PSNOWTOTSWEXY   ! crocus snowheight PSNOWSWE
    REAL, DIMENSION(:,:),   INTENT(INOUT), allocatable  ::   FLOW_ICE
    REAL, DIMENSION(:,:),   INTENT(INOUT), allocatable  ::   FLOW_SNOW
    INTEGER,DIMENSION(:,:), INTENT(INOUT), allocatable  :: GLACINFO
    REAL, DIMENSION(:,:),   INTENT(INOUT), allocatable  :: GLACT
    INTEGER,                INTENT(IN)                  :: crocus_opt
    REAL, DIMENSION(:,:),   INTENT(INOUT), allocatable  :: VIS_ICEALB

    INTEGER                                             :: GLACINFOH,GLACR ! glacier info
    CHARACTER(LEN=11)                                   :: HSNOWRES
    CHARACTER(LEN=3)                                    :: HIMPLICIT_WIND
    LOGICAL                                             :: OGLACIER

    REAL, DIMENSION(1)                :: ZP_RRSNOW ,ZP_PSN3L, ZP_QA
    REAL, DIMENSION(1)                :: ZP_PEW_A_COEF, ZP_PEW_B_COEF
    REAL, DIMENSION(1)                :: ZP_PET_A_COEF, ZP_PEQ_A_COEF
    REAL, DIMENSION(1)                :: ZP_PET_B_COEF, ZP_PEQ_B_COEF
    REAL, PARAMETER                   :: ZTFRZ = 273.16 !freezing/melting point (k)
    REAL                              :: ZFPICE ! snow fraction of precip
    REAL, DIMENSION(1)                :: FOO1, FOO2
    REAL, DIMENSION(:,:), allocatable :: ZP_SNOWSWE
    REAL, DIMENSION(:,:), allocatable :: ZP_SNOWDZ
    REAL, DIMENSION(:,:), allocatable :: ZP_SNOWRHO
    REAL, DIMENSION(:,:), allocatable :: ZP_SNOWHEAT
    REAL, DIMENSION(:,:), allocatable :: ZP_SNOWTEMP
    REAL, DIMENSION(:,:), allocatable :: ZP_SNOWLIQ
    REAL, DIMENSION(:,:), allocatable :: ZP_SNOWGRAN1
    REAL, DIMENSION(:,:), allocatable :: ZP_SNOWGRAN2
    REAL, DIMENSION(:,:), allocatable :: ZP_SNOWHIST
    REAL, DIMENSION(:,:), allocatable :: ZP_SNOWAGE

    REAL, DIMENSION(1) :: ZP_VMOD,ZP_RHOA, ZP_UREF
    REAL, DIMENSION(1) :: ZP_EXNS, ZP_EXNA, ZP_DIRCOSZW
    REAL, DIMENSION(1) :: ZP_Z0NAT
    REAL, DIMENSION(1) :: ZP_Z0HNAT, ZP_ALB, ZP_SOILCOND
    REAL, DIMENSION(1) :: ZP_D_G,ZP_SNOWALB, ZP_VIS_ICEALB
    REAL, DIMENSION(1) :: ZFLOW_ICE, ZFLOW_SNOW
    REAL, DIMENSION(1) :: ZP_THRUFAL, ZP_GRNDFLUX, ZP_EVAPCOR
    REAL, DIMENSION(1) :: ZP_RNSNOW, ZP_GFLUXSNOW, ZP_HPSNOW
    REAL, DIMENSION(1) :: ZP_LEL3L, ZP_SNDRIFT, ZP_RI
    REAL, DIMENSION(1) :: ZP_SNOWHMASS
    REAL, DIMENSION(1) :: ZP_USTARSNOW
    REAL, DIMENSION(1) :: ZP_VEGTYPE, ZP_ZENITH
    REAL, DIMENSION(1) :: ZP_PS        ! PSFC
    REAL, DIMENSION(1) :: ZP_SRSNOW    ! QSNOW
    REAL, DIMENSION(1) :: ZP_TA        ! T_ML
    REAL, DIMENSION(1) :: ZP_TG        ! TG?
    REAL, DIMENSION(1) :: ZP_SW_RAD    ! SWDN
    REAL, DIMENSION(1) :: ZP_LW_RAD    ! LWDN
    REAL, DIMENSION(1) :: ZP_ZREF      ! Z_ML
    REAL, DIMENSION(1) :: ZP_HSNOW     ! FSH
    REAL, DIMENSION(1) :: ZP_LES3L     ! FGEV
    REAL, DIMENSION(1) :: ZP_EVAP      ! ESOIL
    REAL, DIMENSION(1) :: ZP_EMISNOW   ! EMISSI
    REAL, DIMENSION(1) :: ZP_CDSNOW    ! CM
    REAL, DIMENSION(1) :: ZP_CHSNOW    ! CH
    REAL, DIMENSION(1) :: ZP_QS        ! QSFC1D
    REAL, DIMENSION(1) :: ZP_Z0EFF     ! effective roughness

    !     Added some temporary fields for using noahmp-glacier to calculate 2m variables
    REAL, DIMENSION(-2:NSOIL) :: ZSNSOT ! snow layer depth [m]
    REAL, DIMENSION(-2:0)     :: SNICET ! snow layer ice [mm]
    REAL, DIMENSION(-2:0)     :: SNLIQT ! snow layer liquid water [mm]

    REAL, DIMENSION(-2:NSOIL) :: ZSNSOTH ! snow layer depth [m]
    REAL, DIMENSION(-2:0)     :: SNICEH  ! snow layer ice [mm]
    REAL, DIMENSION(-2:NSOIL) :: ZSNSOH  ! snow layer depth [m]

    !     Added  fluxes:
    REAL,DIMENSION(1) :: ZP_FSA_CROCUS
    REAL,DIMENSION(1) :: ZP_FSR_CROCUS
    REAL,DIMENSION(1) :: ZP_FIRA_CROCUS

    LOGICAL      :: OSNOWDRIFT, OSNOWDRIFT_SUBLIM ! activate snowdrift, sublimation during drift
    LOGICAL      :: OSNOW_ABS_ZENITH ! activate parametrization of solar absorption for polar regions
    REAL         :: SWEH
    REAL         :: SNDPTHH
    REAL         :: SNEQVOH
    CHARACTER(3) :: HSNOWMETAMO, HSNOWRAD

    REAL,   DIMENSION( 1:NSOIL)    :: liq_frac

    ! crocus, assign mock local values, and initialize. These values should be in namelist.
    if (crocus_opt /= 0) then
       allocate(ZP_SNOWSWE  (1,1:act_level))
       allocate(ZP_SNOWDZ   (1,1:act_level))
       allocate(ZP_SNOWRHO  (1,1:act_level))
       allocate(ZP_SNOWHEAT (1,1:act_level))
       allocate(ZP_SNOWTEMP (1,1:act_level))
       allocate(ZP_SNOWLIQ  (1,1:act_level))
       allocate(ZP_SNOWGRAN1(1,1:act_level))
       allocate(ZP_SNOWGRAN2(1,1:act_level))
       allocate(ZP_SNOWHIST (1,1:act_level))
       allocate(ZP_SNOWAGE  (1,1:act_level))

       HSNOWRES = 'ISBA-SNOW3L'  !   'DEF'
       OGLACIER = .false.
       !      OGLACIER = .true.
       HIMPLICIT_WIND = 'NEW'
       !      HIMPLICIT_WIND = 'OLD'

       ZP_VEGTYPE = 1.0
       OSNOWDRIFT = .true.
       OSNOWDRIFT_SUBLIM = .true.
       OSNOW_ABS_ZENITH = .false.
       HSNOWMETAMO = 'B92'
       !      HSNOWMETAMO = 'C13'
       HSNOWRAD         = 'B92'
       !     HSNOWRAD         = 'TAR2'
    else
       allocate(ZP_SNOWDZ   (1,1)) ! This is allocated to avoid SEGFAULT when crocus_opt=0 and
       ZP_SNOWDZ = 0.0             ! optimization is off
    end if ! crocus_opt /= 0

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

    CALL NOAHMP_OPTIONS(IDVEG  ,IOPT_CRS  ,IOPT_BTR  ,IOPT_RUN  ,IOPT_SFC  ,IOPT_FRZ , &
                     IOPT_INF  ,IOPT_RAD  ,IOPT_ALB  ,IOPT_SNF  ,IOPT_TBOT, IOPT_STC , &
		     IOPT_RSF  ,IOPT_SOIL ,IOPT_PEDO ,IOPT_CROP, IOPT_IMPERV )

    IPRINT    =  .false.                     ! debug printout

    YEARLEN = 365                            ! find length of year for phenology (also S Hemisphere)
    if (mod(YR,4) == 0) then
       YEARLEN = 366
       if (mod(YR,100) == 0) then
          YEARLEN = 365
          if (mod(YR,400) == 0) then
             YEARLEN = 366
          endif
       endif
    endif

    ZSOIL(1) = -DZS(1)                    ! depth to soil interfaces (<0) [m]
    DO K = 2, NSOIL
       ZSOIL(K) = -DZS(K) + ZSOIL(K-1)
    END DO

    JLOOP : DO J=jts,jte
       IF(ITIMESTEP == 1)THEN
          DO I=its,ite
             IF((XLAND(I,J)-1.5) >= 0.) THEN    ! Open water case
                IF(XICE(I,J) == 1. .AND. IPRINT) PRINT *,' sea-ice at water point, I=',I,'J=',J
                SMSTAV(I,J) = 1.0
                SMSTOT(I,J) = 1.0
                DO K = 1, NSOIL
                   SMOIS(I,K,J) = 1.0
                    TSLB(I,K,J) = 273.16
                ENDDO
             ELSE
                IF(XICE(I,J) == 1.) THEN        ! Sea-ice case
                   SMSTAV(I,J) = 1.0
                   SMSTOT(I,J) = 1.0
                   DO K = 1, NSOIL
                      SMOIS(I,K,J) = 1.0
                   ENDDO
                ENDIF
             ENDIF
          ENDDO
       ENDIF                                                               ! end of initialization over ocean


!-----------------------------------------------------------------------
   ILOOP : DO I = its, ite
      ! Reset GLACR and ICE at the start of every iteration
      GLACR = 0
      ICE   = 0

    IF (XICE(I,J) >= XICE_THRES) THEN
       ICE = 1                            ! Sea-ice point

       SH2O  (i,1:NSOIL,j) = 1.0
       XLAIXY(i,j)         = 0.01

       CYCLE ILOOP ! Skip any processing at sea-ice points

    ELSE

       IF((XLAND(I,J)-1.5) >= 0.) CYCLE ILOOP   ! Open water case

!     2D to 1D

! IN only

       COSZ   = COSZIN  (I,J)                         ! cos zenith angle []
       LAT    = XLATIN  (I,J)                         ! latitude [rad]
       Z_ML   = 0.5*DZ8W(I,1,J)                       ! DZ8W: thickness of full levels; ZLVL forcing height [m]
       VEGTYP = IVGTYP(I,J)                           ! vegetation type
       SOILTYP= ISLTYP(I,J)                           ! soil type
       FVEG   = VEGFRA(I,J)/100.                      ! vegetation fraction [0-1]
       FVGMAX = VEGMAX (I,J)/100.                     ! Vegetation fraction annual max [0-1]
       TBOT = TMN(I,J)                                ! Fixed deep soil temperature for land
       T_ML   = T3D(I,1,J)                            ! temperature defined at intermediate level [K]
       Q_ML   = QV3D(I,1,J)/(1.0+QV3D(I,1,J))         ! convert from mixing ratio to specific humidity [kg/kg]
       U_ML   = U_PHY(I,1,J)                          ! u-wind at interface [m/s]
       V_ML   = V_PHY(I,1,J)                          ! v-wind at interface [m/s]
       SWDN   = SWDOWN(I,J)                           ! shortwave down from SW scheme [W/m2]
       LWDN   = GLW(I,J)                              ! total longwave down from LW scheme [W/m2]
       P_ML   =(P8W3D(I,KTS+1,J)+P8W3D(I,KTS,J))*0.5  ! surface pressure defined at intermediate level [Pa]
	                                              !    consistent with temperature, mixing ratio
       PSFC   = P8W3D(I,1,J)                          ! surface pressure defined a full levels [Pa]
       PRCP   = PRECIP_IN (I,J) / DT                  ! timestep total precip rate (glacier) [mm/s]! MB: v3.7

       CROPTYPE = 0
       if (crocus_opt /= 0) then
          ZFLOW_ICE = 0
          ZFLOW_SNOW = 0
       end if

       IF (PRESENT(MP_RAINC) .AND. PRESENT(MP_RAINNC) .AND. PRESENT(MP_SHCV) .AND. &
           PRESENT(MP_SNOW)  .AND. PRESENT(MP_GRAUP)  .AND. PRESENT(MP_HAIL)   ) THEN

         PRCPCONV  = MP_RAINC (I,J)/DT                ! timestep convective precip rate [mm/s]     ! MB: v3.7
         PRCPNONC  = MP_RAINNC(I,J)/DT                ! timestep non-convective precip rate [mm/s] ! MB: v3.7
         PRCPSHCV  = MP_SHCV(I,J)  /DT                ! timestep shallow conv precip rate [mm/s]   ! MB: v3.7
         PRCPSNOW  = MP_SNOW(I,J)  /DT                ! timestep snow precip rate [mm/s]           ! MB: v3.7
         PRCPGRPL  = MP_GRAUP(I,J) /DT                ! timestep graupel precip rate [mm/s]        ! MB: v3.7
         PRCPHAIL  = MP_HAIL(I,J)  /DT                ! timestep hail precip rate [mm/s]           ! MB: v3.7

         PRCPOTHR  = PRCP - PRCPCONV - PRCPNONC - PRCPSHCV ! take care of other (fog) contained in rainbl
	 PRCPOTHR  = MAX(0.0,PRCPOTHR)
	 PRCPNONC  = PRCPNONC + PRCPOTHR
         PRCPSNOW  = PRCPSNOW + SR(I,J)  * PRCPOTHR
       ELSE
         PRCPCONV  = 0.
         PRCPNONC  = PRCP
         PRCPSHCV  = 0.
         PRCPSNOW  = SR(I,J) * PRCP
         PRCPGRPL  = 0.
         PRCPHAIL  = 0.
       ENDIF

! IN/OUT fields

       ISNOW                 = ISNOWXY (I,J)                ! snow layers []
       SMC  (      1:NSOIL)  = SMOIS   (I,      1:NSOIL,J)  ! soil total moisture [m3/m3]
       SMH2O(      1:NSOIL)  = SH2O    (I,      1:NSOIL,J)  ! soil liquid moisture [m3/m3]
       STC  (-NSNOW+1:    0) = TSNOXY  (I,-NSNOW+1:    0,J) ! snow temperatures [K]
       STC  (      1:NSOIL)  = TSLB    (I,      1:NSOIL,J)  ! soil temperatures [K]
       SWE                   = SNOW    (I,J)                ! snow water equivalent [mm]
       SNDPTH                = SNOWH   (I,J)                ! snow depth [m]
       QSFC1D                = QSFC    (I,J)

! INOUT (with no Noah LSM equivalent)

       TV                    = TVXY    (I,J)                ! leaf temperature [K]
       TG                    = TGXY    (I,J)                ! ground temperature [K]
       CANLIQ                = CANLIQXY(I,J)                ! canopy liquid water [mm]
       CANICE                = CANICEXY(I,J)                ! canopy frozen water [mm]
       EAH                   = EAHXY   (I,J)                ! canopy vapor pressure [Pa]
       TAH                   = TAHXY   (I,J)                ! canopy temperature [K]
       CM                    = CMXY    (I,J)                ! avg. momentum exchange (MP only) [m/s]
       CH                    = CHXY    (I,J)                ! avg. heat exchange (MP only) [m/s]
       FWET                  = FWETXY  (I,J)                ! canopy fraction wet or snow
       SNEQVO                = SNEQVOXY(I,J)                ! SWE previous timestep
       ALBOLD                = ALBOLDXY(I,J)                ! albedo previous timestep, for snow aging
       QSNOW                 = QSNOWXY (I,J)                ! snow falling on ground
       QRAIN                 = QRAINXY (I,J)                ! rain falling on ground
       WSLAKE                = WSLAKEXY(I,J)                ! lake water storage (can be neg.) (mm)
       ZWT                   = ZWTXY   (I,J)                ! depth to water table [m]
       WA                    = WAXY    (I,J)                ! water storage in aquifer [mm]
       WT                    = WTXY    (I,J)                ! water in aquifer&saturated soil [mm]
       ZSNSO(-NSNOW+1:NSOIL) = ZSNSOXY (I,-NSNOW+1:NSOIL,J) ! depth to layer interface
       SNICE(-NSNOW+1:    0) = SNICEXY (I,-NSNOW+1:    0,J) ! snow layer ice content
       SNLIQ(-NSNOW+1:    0) = SNLIQXY (I,-NSNOW+1:    0,J) ! snow layer water content
       LFMASS                = LFMASSXY(I,J)                ! leaf mass
       RTMASS                = RTMASSXY(I,J)                ! root mass
       STMASS                = STMASSXY(I,J)                ! stem mass
       WOOD                  = WOODXY  (I,J)                ! mass of wood (incl. woody roots) [g/m2]
       STBLCP                = STBLCPXY(I,J)                ! stable carbon pool
       FASTCP                = FASTCPXY(I,J)                ! fast carbon pool
       PLAI                  = XLAIXY  (I,J)                ! leaf area index [-] (no snow effects)
       PSAI                  = XSAIXY  (I,J)                ! stem area index [-] (no snow effects)
       TAUSS                 = TAUSSXY (I,J)                ! non-dimensional snow age
       SMCEQ(       1:NSOIL) = SMOISEQ (I,       1:NSOIL,J)
       SMCWTD                = SMCWTDXY(I,J)
       RECH                  = 0.
       DEEPRECH              = 0.

       SLOPETYP     = 1                               ! set underground runoff slope term
       IST          = 1                               ! MP surface type: 1 = land; 2 = lake
       SOILCOLOR    = 4                               ! soil color: assuming a middle color category ?????????

       IF(SOILTYP == 14 .AND. XICE(I,J) == 0.) THEN
          IF(IPRINT) PRINT *, ' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
          IF(IPRINT) PRINT *, i,j,'RESET SOIL in surfce.F'
          SOILTYP = 7
       ENDIF

#ifdef SPATIAL_SOIL
       parameters%bexp   = BEXP_3D  (I,1:NSOIL,J) ! C-H B exponent
       parameters%smcdry = SMCDRY_3D(I,1:NSOIL,J) ! Soil Moisture Limit: Dry
       parameters%smcwlt = SMCWLT_3D(I,1:NSOIL,J) ! Soil Moisture Limit: Wilt
       parameters%smcref = SMCREF_3D(I,1:NSOIL,J) ! Soil Moisture Limit: Reference
       parameters%smcmax = SMCMAX_3D(I,1:NSOIL,J) ! Soil Moisture Limit: Max
       parameters%dksat  = DKSAT_3D (I,1:NSOIL,J) ! Saturated Soil Conductivity
       parameters%dwsat  = DWSAT_3D (I,1:NSOIL,J) ! Saturated Soil Diffusivity
       parameters%psisat = PSISAT_3D(I,1:NSOIL,J) ! Saturated Matric Potential
       parameters%quartz = QUARTZ_3D(I,1:NSOIL,J) ! Soil quartz content
       parameters%refdk  = REFDK_2D (I,J)         ! Reference Soil Conductivity
       parameters%refkdt = REFKDT_2D(I,J)         ! Soil Infiltration Parameter
       parameters%slope  = SLOPE_2D (I,J)         ! Soil Drainage Parameter
       parameters%cwpvt  = CWPVT_2D(I,J)          ! Canopy wind parameter
       parameters%vcmx25 = VCMX25_2D(I,J)         ! VCmax at 25C
       parameters%mp     = MP_2D(I,J)             ! Slope of Ball-Berry rs-P relationship
       parameters%hvt    = HVT_2D(I,J)            ! Canopy Height
       parameters%mfsno  = MFSNO_2D(I,J)          ! Snow cover m parameter
       parameters%rsurf_exp = RSURFEXP_2D(I,J)    ! exponent in the shape parameter for soil resistance option 1
       parameters%axaj   = AXAJ_2D(I,J)           ! Xinanjiang: Tension water distribution inflection parameter [-]
       parameters%bxaj   = BXAJ_2D(I,J)           ! Xinanjiang: Tension water distribution shape parameter [-]
       parameters%xxaj   = XXAJ_2D(I,J)           ! Xinanjiang: Free water distribution shape parameter [-]
       parameters%imperv = IMPERV_2D(I,J)         ! impervious fraction
       parameters%ssi    = SSI_2D(I,J)            ! liquid water holding capacity for snowpack (m3/m3)
       parameters%snow_ret_fac = SNOWRETFAC_2D(I,J)   ! snowpack water release timescale factor (1/s)
       parameters%tau0   = TAU0_2D(I,J)           ! tau0 from Yang97 eqn. 10a
       parameters%rsurf_snow = RSURFSNOW_2D(I,J)  ! surface resistence for snow [s/m]
       parameters%scamax = SCAMAX_2D(I,J)         ! maximum fractional snow covered area (0.0-1.0)
#endif
       CALL TRANSFER_MP_PARAMETERS(VEGTYP,SOILTYP,SLOPETYP,SOILCOLOR,CROPTYPE,parameters,iopt_imperv)

       GRAIN = 0.0                ! mass of grain [g/m2]
       GDD   = 0.0                ! growing degree days
       PGS   = 0                  ! crop growth stage

       if (HVT_TABLE(VEGTYP) > 0.0) then
         parameters%hvb = parameters%hvt * HVB_TABLE(VEGTYP) / HVT_TABLE(VEGTYP)
       else
         parameters%hvb = 0.0
       endif

       !ADCHANGE: Add some sanity checks in case calibration knocks the order of these out of sequence.
       !The min diffs were pulled from the existing SOILPARM.TBL defaults.
       !Currently water is 0, so enforcing 0 as the absolute min.
       parameters%smcmax = min(parameters%smcmax, 1.0)
       parameters%smcref = max(min(parameters%smcref, parameters%smcmax - 0.01), 0.0)
       parameters%smcwlt = max(min(parameters%smcwlt, parameters%smcref - 0.01), 0.0)
       parameters%smcdry = max(min(parameters%smcdry, parameters%smcref - 0.01), 0.0)

! Initialized local

       FICEOLD = 0.0
       FICEOLD(ISNOW+1:0) = SNICEXY(I,ISNOW+1:0,J) &  ! snow ice fraction
           /(SNICEXY(I,ISNOW+1:0,J)+SNLIQXY(I,ISNOW+1:0,J))
       CO2PP  = CO2_TABLE * P_ML                      ! partial pressure co2 [Pa]
       O2PP   = O2_TABLE  * P_ML                      ! partial pressure  o2 [Pa]
       FOLN   = 1.0                                   ! for now, set to nitrogen saturation
       QC     = undefined_value                       ! test dummy value
       PBLH   = undefined_value                       ! test dummy value ! PBL height
       DZ8W1D = DZ8W (I,1,J)                          ! thickness of atmospheric layers

! Initialize crocus, !  Do this before NOAHMP, so that outputs (such as TG) from Noah-mp glacier affect Crocus this step
       if (crocus_opt /= 0) then
          GLACINFOH = GLACINFO(I,J)   ! vegetation type

          ZP_SNOWSWE(1,1:act_level) = PSNOWSWEXY(I,1:act_level,J)
          ZP_SNOWRHO(1,1:act_level) = PSNOWRHOXY(I,1:act_level,J)
          ZP_SNOWHEAT(1,1:act_level) = PSNOWHEATXY(I,1:act_level,J)
          ZP_SNOWGRAN1(1,1:act_level) = PSNOWGRAN1XY(I,1:act_level,J)
          ZP_SNOWGRAN2(1,1:act_level) = PSNOWGRAN2XY(I,1:act_level,J)
          ZP_SNOWHIST(1,1:act_level) = PSNOWHISTXY(I,1:act_level,J)
          ZP_SNOWAGE(1,1:act_level)= PSNOWAGEXY(I,1:act_level,J)
          ZP_SNOWLIQ(1,1:act_level) = PSNOWLIQXY(I,1:act_level,J)
          ZP_SNOWTEMP(1,1:act_level) = PSNOWTEMPXY(I,1:act_level,J)
          ZP_SNOWDZ(1,1:act_level) = PSNOWDZXY(I,1:act_level,J)
          ZP_SNOWALB(1) =  PSNOWALBXY(I,J)
          ZP_EVAP(1)    =  0 !PSNOWEVAPXY(I,J)
          ZP_vis_icealb(1) = VIS_ICEALB(I,J)

          ZP_ZENITH    = ACOS(COSZIN(I,J))

          ! Currently the surfex naming convention is being used
          ZP_PS(1)  =   PSFC
          ZP_SRSNOW(1) = QSNOW
          ZP_TA(1) = T_ML
          ZP_TG(1)  =  TSLB(i,1,j)
          IF (ZP_SNOWTEMP(1,3) .ne. 0 ) then
             ZP_TG(1)  =  (TSLB(i,1,j)+ZP_SNOWTEMP(1,3))/2.
          end if
          ZP_RHOA(1) = P_ML/(T_ML*287.058)
          ZP_SW_RAD(1)  = SWDN
          ZP_LW_RAD(1)  = LWDN
          ZP_ZREF(1)       = Z_ML
          ZP_HSNOW(1)   = FSH
          ZP_LES3L(1)      = FGEV
          ZP_EVAP(1)       = ESOIL
          ZP_EMISNOW(1) = EMISSI
          ZP_CDSNOW(1)  = CM
          ZP_CHSNOW(1)  = CH
          ZP_QS(1)             = QSFC1D
          ZP_QA(1)             =Q_ML
          ZP_ALB(1)            = ALBEDO(I,J) !soil/vegetation albedo. Is SALB combined snow/soil albedo?
          ZP_EXNA(1)             = (P_ML/1.e5)** (287.04/1004.6) !(P/XP00)**(XRD/XCPD)
          ZP_EXNS(1)             = (PSFC/1.e5)** (287.04/1004.6) !(P/XP00)**(XRD/XCPD)
          ZP_DIRCOSZW(1)   = 1.0    !(for now we use a flat surface)
          ZP_UREF(1)            = Z_ML

          ZP_Z0EFF(1)          =  XZ0SN
          ZP_Z0NAT(1)         =  XZ0SN
          ZP_Z0HNAT(1)      = XZ0HSN

          !     in offline mode, A_COEF = 0, and B_COEF = variable of interest at atmospheric level  (Boone et al, 2017; The interaction between soil-biasphere-------)
          ZP_PSN3L(1)    = 1.0
          ZP_VMOD(1)     = 1.0      !(assume terrain is flat)
          ZP_D_G(1)         = 0.1   ! should read from namelits (soil_thick_input)

          !     Fluxes
          ZP_FSA_CROCUS(1)  = 0
          ZP_FSR_CROCUS(1) = 0
          ZP_FIRA_CROCUS(1) = 0

          ZP_SOILCOND(1) = 0.32333+(0.10073*(-ZSNSO(1)))

          ZP_PEW_A_COEF(1)  = 0.0
          ZP_PEW_B_COEF(1)  = sqrt(U_ML**2+V_ML**2)
          ZP_PET_A_COEF(1)    =0.0
          ZP_PEQ_A_COEF(1)   = 0.0
          ZP_PET_B_COEF(1)   = T_ML*(1.e5/P_ML)** (287.04/1004.6)
          ZP_PEQ_B_COEF(1)   = Q_ML

          ZP_VMOD(1)  =      ZP_PEW_B_COEF(1)


          ! Added from module_sf_noahmp_glacier:
          ! partition precipitation into rain and snow (from CANWATER)
          ! Jordan (1991)
          IF(IOPT_SNF == 1 .OR. IOPT_SNF == 4) THEN
             IF(T_ML > ZTFRZ+2.5)THEN
                ZFPICE = 0.
             ELSE
                IF(T_ML <= ZTFRZ+0.5)THEN
                   ZFPICE = 1.0
                ELSE IF(T_ML <= ZTFRZ+2.)THEN
                   ZFPICE = 1.-(-54.632 + 0.2*T_ML)
                ELSE
                   ZFPICE = 0.6
                ENDIF
             ENDIF
          ENDIF

          IF(IOPT_SNF == 2) THEN
             IF(T_ML >= ZTFRZ+2.2) THEN
                ZFPICE = 0.
             ELSE
                ZFPICE = 1.0
             ENDIF
          ENDIF

          IF(IOPT_SNF == 3) THEN
             IF(T_ML >= ZTFRZ) THEN
                ZFPICE = 0.
             ELSE
                ZFPICE = 1.0
             ENDIF
          ENDIF

          !     Hedstrom NR and JW Pomeroy (1998), Hydrol. Processes, 12, 1611-1625
          !     fresh snow density

          ZP_RRSNOW(1)   = PRCP* (1.-ZFPICE)
          ZP_SRSNOW(1)   = PRCP*ZFPICE

          ZP_THRUFAL(:) = 0.0
      end if ! crocus_opt /= 0

       IF(VEGTYP == 25) FVEG = 0.0                  ! Set playa, lava, sand to bare
       IF(VEGTYP == 25) PLAI = 0.0
       IF(VEGTYP == 26) FVEG = 0.0                  ! hard coded for USGS
       IF(VEGTYP == 26) PLAI = 0.0
       IF(VEGTYP == 27) FVEG = 0.0
       IF(VEGTYP == 27) PLAI = 0.0

       ! If running crocus, test if glaicer has has comletely melted, then RUN NoahMP
       IF ( (crocus_opt == 0) .or. (GLACINFOH == 1 .and. ZP_SNOWDZ(1,1) .eq. 0.0)) THEN
          IF ( VEGTYP == ISICE_TABLE ) THEN
             ICE = -1                           ! Land-ice point
             if (crocus_opt /= 0) &
                  GLACR = 1
             CALL NOAHMP_OPTIONS_GLACIER(IOPT_ALB  ,IOPT_SNF  ,IOPT_TBOT, IOPT_STC, IOPT_GLA )

             TBOT = MIN(TBOT,263.15)                      ! set deep temp to at most -10C
             CALL NOAHMP_GLACIER(I,       J,    COSZ,   NSNOW,   NSOIL,      DT, & ! IN : Time/Space/Model-related
                              T_ML,    P_ML,    U_ML,    V_ML,    Q_ML,    SWDN, & ! IN : Forcing
                              PRCP,    LWDN,    TBOT,    Z_ML, FICEOLD,   ZSOIL, & ! IN : Forcing
                              parameters%swe_limit,                              & ! IN : Forcing
                             QSNOW, QRAIN, SNEQVO,  ALBOLD,   CM,   CH,   ISNOW, & ! IN/OUT :
                               SWE,     SMC,   ZSNSO,  SNDPTH,   SNICE,   SNLIQ, & ! IN/OUT :
                                TG,     STC,   SMH2O,   TAUSS,  QSFC1D,          & ! IN/OUT :
                               FSA,     FSR,    FIRA,     FSH,    FGEV,   SSOIL, & ! OUT :
                              TRAD,   ESOIL,   RUNSF,   RUNSB,     SAG,    SALB, & ! OUT :
                             QSNBOT,PONDING,PONDING1,PONDING2,    T2MB,    Q2MB, & ! OUT :
                             EMISSI,  FPICE,    CHB2 &                             ! OUT :
#ifdef WRF_HYDRO
                             , sfcheadrt(i,j)                                      &
#endif
                             )
          ENDIF ! VEGTYP == ISICE_TABLE
       ELSE IF (crocus_opt /= 0 .and. GLACINFOH == 1 .and. ZP_SNOWDZ(1,1) .gt. 0.0) THEN
          ICE = -1               ! Land-ice point
          GLACR= 1
          CALL NOAHMP_OPTIONS_GLACIER(IOPT_ALB  ,IOPT_SNF  ,IOPT_TBOT, IOPT_STC, IOPT_GLA )
          TBOT = MIN(TBOT,263.15) ! set deep temp to at most -10C
                                  ! send snow info from crocus to Noahmp_glacier to calc 2m temps etc
                                  ! move NOAHMP_GLACIER to after CROCUS
                                  ! for temperatre glacier, set top soil temp to 272
          ZP_TG(1) = 272.0
          ! Why is this if block here? Not sure top condition would ever trigger  
          ! given we are already in an if block requiring ZP_SNOWDZ(1,1) .gt. 0.0. (AD)
          IF (ZP_SNOWDZ(1,1) .eq. 0.0 .and. ZP_SRSNOW(1) .eq. 0.0) THEN
          ELSE
            CALL SNOWCRO(HSNOWRES, OGLACIER, HIMPLICIT_WIND,                   & !(IN)
                  ZP_PEW_A_COEF, ZP_PEW_B_COEF,                                & !(IN)
                  ZP_PET_A_COEF, ZP_PEQ_A_COEF, ZP_PET_B_COEF, ZP_PEQ_B_COEF,  & !(IN)
                  ZP_SNOWSWE,ZP_SNOWRHO, ZP_SNOWHEAT, ZP_SNOWALB,              & ! (INOUT) !several layers
                  ZP_SNOWGRAN1, ZP_SNOWGRAN2, ZP_SNOWHIST, ZP_SNOWAGE, DT,     & !(INOUT) PTSTEP: (IN)
                  ZP_PS, ZP_SRSNOW, ZP_RRSNOW ,ZP_PSN3L, ZP_TA, ZP_TG,         & !(IN)
                  ZP_SW_RAD, ZP_QA, ZP_VMOD, ZP_LW_RAD, ZP_RHOA, ZP_UREF,      & !(IN)
                  ZP_EXNS, ZP_EXNA, ZP_DIRCOSZW, ZP_ZREF, ZP_Z0NAT, ZP_Z0EFF,  & !(IN)
                  ZP_Z0HNAT, ZP_ALB, ZP_SOILCOND, ZP_D_G,                      & !(IN)
                  ZP_SNOWLIQ,                                                  & !(OUT)
                  ZP_SNOWTEMP, ZP_SNOWDZ, ZP_THRUFAL, ZP_GRNDFLUX, ZP_EVAPCOR, & !(OUT)
                  ZP_RNSNOW, ZP_HSNOW, ZP_GFLUXSNOW, ZP_HPSNOW, ZP_LES3L,      & !(OUT)
                  ZP_LEL3L, ZP_EVAP, ZP_SNDRIFT, ZP_RI,                        & !(OUT)
                  ZP_EMISNOW, ZP_CDSNOW, ZP_USTARSNOW,                         & !(OUT)
                  ZP_CHSNOW, ZP_SNOWHMASS, ZP_QS,                              & !(OUT)
                  ZP_VEGTYPE, ZP_ZENITH,                                       & !(IN)
                  OSNOWDRIFT,OSNOWDRIFT_SUBLIM,                                & !(IN)
                  OSNOW_ABS_ZENITH, HSNOWMETAMO,HSNOWRAD,                      &
                  act_level, ZP_VIS_ICEALB,                                    &
                  ZP_FSA_CROCUS, ZP_FSR_CROCUS, ZP_FIRA_CROCUS,                &
                  ZFLOW_ICE, ZFLOW_SNOW,                                       &
                  I,J) !(OUT)

             !  Call to NOAHM_GLACIER to after CROCUS to calculte 2m temp and other 2m variables.
             !  But use snow info from CROCUS.

             ZSNSO(-NSNOW:0) = 0.0
             ALBOLD= ZP_SNOWALB(1)
             SNDPTHH = sum(ZP_SNOWDZ)
             SNEQVOH= sum(PSNOWSWEXY(I,:,J))
             SWEH = sum(ZP_SNOWDZ) ! duplicate of SNDPTHH?
             ZSNSOTH(-1)=-SNDPTHH/2
             zsnsoTH(0)=-SNDPTHH
             ZSNSOTH(1)=-SNDPTHH-0.1
             zSNSOTH(2)=ZSNSO(1)-0.3
             zSNSOTH(3)=ZSNSO(2)-0.6
             zSNSOTH(4)=ZSNSO(3)-1.0
             STC(-NSNOW+3)=273.15
             STC(-NSNOW+2)=(ZP_SNOWTEMP(1,1))
             ISNOW=-NSNOW+1 ! (for glacier, always assume a snow layer)
             SNICEH=sum(PSNOWDZXY(I,:,J))*1000./(NSNOW)
             if (SNDPTHH .eq. 0) then
                if(VEGTYP .eq. 24) then
                   ISNOW=1
                else
                   ISNOW=0
                endif
             endif

             CALL NOAHMP_GLACIER(     I,       J,    COSZ,   NSNOW,   NSOIL,      DT, & ! IN : Time/Space/Model-related
                  T_ML,    P_ML,    U_ML,    V_ML,    Q_ML,    SWDN,   & ! IN : Forcing
                  PRCP,    LWDN,    TBOT,    Z_ML, FICEOLD,   ZSOIL,   & ! IN : Forcing
                  parameters%swe_limit                               , & ! IN : Forcing
                  QSNOW, QRAIN, SNEQVOH,  ALBOLD,   CM,   CH,   ISNOW, & ! IN/OUT :
!                 SWE,     SMC,   ZSNSO,  SNDPTH,   SNICE,   SNLIQ,    & ! IN/OUT :
                  SWEH,     SMC,   ZSNSOTH,  SNDPTHH,   SNICEH,   SNLIQ, & ! IN/OUT :
                  TG,     STC,   SMH2O,   TAUSS,  QSFC1D,              & ! IN/OUT :
                  FSA,     FSR,    FIRA,     FSH,    FGEV,   SSOIL,    & ! OUT :
                  TRAD,   ESOIL,   RUNSF,   RUNSB,     SAG,    SALB,   & ! OUT :
                  QSNBOT,PONDING,PONDING1,PONDING2,    T2MB,    Q2MB,  & ! OUT :
                  EMISSI,  FPICE,    CHB2                              & ! OUT :
#ifdef WRF_HYDRO
                  , sfcheadrt(i,j)                                     &
#endif
             )


             RUNSF=ZP_THRUFAL(1)
             RUNSB=0.            ! assume that there are no subsurface runoff under the glacier.
#ifdef WRF_HYDRO
             RUNSF=ZP_THRUFAL(1)+sfcheadrt(i,j)/DT
#endif

             SNLIQ = 0
             FOO1 = MAX((ZP_LEL3L+ZP_LES3L)/2.8440E06,0.)
             FOO2 = ABS(MIN((ZP_LEL3L+ZP_LES3L)/2.8440E06,0.))

             FSA = ZP_FSA_CROCUS(1)
             FSR = ZP_FSR_CROCUS(1)
             FIRA = ZP_FIRA_CROCUS(1)
             FSH = ZP_HSNOW(1)
             FGEV = ZP_LEL3L(1)+ZP_LES3L(1)
             SSOIL=ZP_GRNDFLUX(1)
             ESOIL = FOO1(1)-FOO2(1)

             TG=ZP_SNOWTEMP(1,1)
             STC   (      1:NSOIL)=273.16

             ! If a crocus glacier cell has melted out, need to reset so NoahMP can take over.
             ! Note that NOAHMP_GLACIER call above updates SNDPTHH, so we are checking ZP_SNOWDZ instead.
             ! Also note that Crocus can generate VERY small negative ZP_SNOWDZ, so we use a <= test.
             ! The negative outputs should probably be caught and fixed in Crocus.
             if (ZP_SNOWDZ(1,1) .le. 0.0) then
               print *, "Crocus glacier cell has melted out. Switching to NoahMP.", i, j
               GLACR = 0
               ! Resetting soil moisture to max possible for a soil (smc=1 while glacier).
               ! This will not conserve mass, but avoids unphysical condition when we transition to NoahMP.
               ! Scaling liquid content based on original ratio (likely small since frozen).
               if (all(smc .gt. 0.)) then
                 liq_frac = smh2o/smc
                 smc = min(smc, parameters%smcmax)
                 smh2o = smc*liq_frac
               endif
             end if

             !----------------
          ENDIF ! If not (ZP_SNOWDZ(1,1) .eq. 0.0 .and. ZP_SRSNOW(1) .eq. 0.0)
                ! If crocus is not running, use the original NOAH-MP glacier
       ELSE IF (crocus_opt /= 0 .and. &
            GLACINFOH .ne. 1 .and. &
            VEGTYP == ISICE_TABLE) THEN ! all other glacier grid points
          ICE = -1
          GLACR= 1
          CALL NOAHMP_OPTIONS_GLACIER(IOPT_ALB  ,IOPT_SNF  ,IOPT_TBOT, IOPT_STC, IOPT_GLA )
          TBOT = MIN(TBOT,263.15) ! set deep temp to at most -10C

          CALL NOAHMP_GLACIER(     I,       J,    COSZ,   NSNOW,   NSOIL,      DT, & ! IN : Time/Space/Model-related
               T_ML,    P_ML,    U_ML,    V_ML,    Q_ML,    SWDN,  & ! IN : Forcing
               PRCP,    LWDN,    TBOT,    Z_ML, FICEOLD,   ZSOIL,  & ! IN : Forcing
               parameters%swe_limit,                               &
               QSNOW, QRAIN, SNEQVO,  ALBOLD,   CM,   CH,   ISNOW, & ! IN/OUT :
               SWE,     SMC,   ZSNSO,  SNDPTH,   SNICE,   SNLIQ,   & ! IN/OUT :
               TG,     STC,   SMH2O,   TAUSS,  QSFC1D,             & ! IN/OUT :
               FSA,     FSR,    FIRA,     FSH,    FGEV,   SSOIL,   & ! OUT :
               TRAD,   ESOIL,   RUNSF,   RUNSB,     SAG,    SALB,  & ! OUT :
               QSNBOT,PONDING,PONDING1,PONDING2,    T2MB,    Q2MB, & ! OUT :
               EMISSI,  FPICE,    CHB2                             & ! OUT :
#ifdef WRF_HYDRO
               , sfcheadrt(i,j)                                    &
#endif
               )
       ENDIF

       IF ((crocus_opt == 0 .and. VEGTYP == ISICE_TABLE) .or. GLACR == 1) THEN
         FSNO   = 1.0
         TV     = undefined_value     ! Output from standard Noah-MP undefined for glacier points
         TGB    = TG
         CANICE = undefined_value
         CANLIQ = undefined_value
         EAH    = undefined_value
         TAH    = undefined_value
         FWET   = undefined_value
         WSLAKE = undefined_value
         ZWT    = undefined_value
         WA     = undefined_value
         WT     = undefined_value
         LFMASS = undefined_value
         RTMASS = undefined_value
         STMASS = undefined_value
         WOOD   = undefined_value
         STBLCP = undefined_value
         FASTCP = undefined_value
         PLAI   = undefined_value
         PSAI   = undefined_value
         T2MV   = undefined_value
         Q2MV   = undefined_value
         NEE    = undefined_value
         GPP    = undefined_value
         NPP    = undefined_value
         FVEGMP = 0.0
         ECAN   = undefined_value
         ETRAN  = undefined_value
         APAR   = undefined_value
         PSN    = undefined_value
         SAV    = undefined_value
         RSSUN  = undefined_value
         RSSHA  = undefined_value
         BGAP   = undefined_value
         WGAP   = undefined_value
         TGV    = undefined_value
         CHV    = undefined_value
         CHB    = CH
         IRC    = undefined_value
         IRG    = undefined_value
         SHC    = undefined_value
         SHG    = undefined_value
         EVG    = undefined_value
         GHV    = undefined_value
         IRB    = FIRA
         SHB    = FSH
         EVB    = FGEV
         GHB    = SSOIL
         TR     = undefined_value
         EVC    = undefined_value
         CHLEAF = undefined_value
         CHUC   = undefined_value
         CHV2   = undefined_value
         FCEV   = undefined_value
         FCTR   = undefined_value
         Z0WRF  = 0.002
         QFX(I,J) = ESOIL
         LH (I,J) = FGEV
         ALBSND = undefined_value
         ALBSNI = undefined_value

      ELSE
         ICE=0 ! Neither sea ice or land ice.

         if (crocus_opt /= 0 .and. GLACINFOH == 1) then
            IF (PLAI .lt. -1000) THEN
               PLAI = 0.5            ! trude test. Glacier points does not have a LAI value (-1e+36). Give LAI a value.
            endif
            IF (TAH .lt. -1000) THEN
               TAH = 273.2           ! trude test. Glacier points does not have a TAH value (-1e+36). Give TAH a value after glaicer has melted.
            endif
            IF (EAH .lt. -1000) THEN
               EAH = 750             ! trude test. Glacier points does not have a EAH value (-1e+36). Give EAH a value  after glaicer has melted.
            endif
            IF (FWET .lt. -1000) THEN
               FWET = 0.05           ! trude test. Glacier points does not have a FWET value (-1e+36). Give FWET a value  after glaicer has melted.
            endif
            IF (CANLIQ .lt. -1000) THEN
               CANLIQ = 0.05         ! trude test. Glacier points does not have a CANLIQ value (-1e+36). Give CANLIQ a value  after glaicer has melted.
            endif
            IF (CANICE .lt. -1000) THEN
               CANICE = 0.00         ! trude test. Glacier points does not have a CANICE value (-1e+36). Give CANICE a value  after glaicer has melted.
            endif
            IF (TV .lt. -1000) THEN
               TV = T_ML +0.5        ! trude test. Glacier points does not have a TV value (-1e+36). Give TV a value  after glaicer has melted.
            endif
            IF (TG .le. -1000) THEN
               TG = T_ML +1.         ! trude test. Glacier points does not have a TG value (-1e+36). Give TG a value  after glaicer has melted.
            endif
            ! IF (ZSNSO(1) .le. -1000) THEN
            !    ZSNSO(:) = 0.       ! trude test. Glacier points does not have a ZSNSO(1) value (-1e+36). Give ZSNSO(1) a value  after glaicer has melted.
            ! endif
            ! next varables could probably be defined from input somewhere
            IF (ZWT .lt. -1000) THEN
               ZWT = 2.5             ! trude test. Glacier points does not have a LAI value (-1e+36). Give LAI a value.
            endif
            IF (WA .lt. -1000) THEN
               WA = 4900.            ! trude test. Glacier points does not have a LAI value (-1e+36). Give LAI a value.
            endif
            IF (WT .lt. -1000) THEN
               WT = 4900.            ! trude test. Glacier points does not have a LAI value (-1e+36). Give LAI a value.
            endif
            IF (WSLAKE .lt. -1000) THEN
               WSLAKE = 0.           ! trude test. Glacier points does not have a LAI value (-1e+36). Give LAI a value.
            endif
            IF (LFMASS .lt. -1000) THEN
               LFMASS = 6.68219      ! trude test. Glacier points does not have a LAI value (-1e+36). Give LAI a value.
            endif
            IF (RTMASS .lt. -1000) THEN
               RTMASS = 500.         ! trude test. Glacier points does not have a LAI value (-1e+36). Give LAI a value.
            endif
            IF (STMASS .lt. -1000) THEN
               STMASS = 17.82        ! trude test. Glacier points does not have a LAI value (-1e+36). Give LAI a value.
            endif
            IF (WOOD .lt. -1000) THEN
               WOOD = 500.           ! trude test. Glacier points does not have a LAI value (-1e+36). Give LAI a value.
            endif
            IF (STBLCP .lt. -1000) THEN
               STBLCP = 1000.        ! trude test. Glacier points does not have a LAI value (-1e+36). Give LAI a value.
            endif
            IF (FASTCP .lt. -1000) THEN
               FASTCP = 1000.        ! trude test. Glacier points does not have a LAI value (-1e+36). Give LAI a value.
            endif
            IF (PSAI .lt. -1000) THEN
               PSAI = 0.49           ! trude test. Glacier points does not have a LAI value (-1e+36). Give LAI a value.
            endif
         end if ! crocus_opt /= 0


         CALL NOAHMP_SFLX (parameters, &
            I       , J       , LAT     , YEARLEN , JULIAN  , COSZ    , & ! IN : Time/Space-related
            DT      , DX      , DZ8W1D  , NSOIL   , ZSOIL   , NSNOW   , & ! IN : Model configuration
            FVEG    , FVGMAX  , VEGTYP  , ICE     , IST     , CROPTYPE, & ! IN : Vegetation/Soil characteristics
            SMCEQ   ,                                                   & ! IN : Vegetation/Soil characteristics
            T_ML    , P_ML    , PSFC    , U_ML    , V_ML    , Q_ML    , & ! IN : Forcing
            QC      , SWDN    , LWDN    ,                               & ! IN : Forcing
	    PRCPCONV, PRCPNONC, PRCPSHCV, PRCPSNOW, PRCPGRPL, PRCPHAIL, & ! IN : Forcing
            TBOT    , CO2PP   , O2PP    , FOLN    , FICEOLD , Z_ML    , & ! IN : Forcing
            ALBOLD  , SNEQVO  ,                                         & ! IN/OUT :
            STC     , SMH2O   , SMC     , TAH     , EAH     , FWET    , & ! IN/OUT :
            CANLIQ  , CANICE  , TV      , TG    , QSFC1D, QSNOW, QRAIN, & ! IN/OUT :
            ISNOW   , ZSNSO   , SNDPTH  , SWE     , SNICE   , SNLIQ   , & ! IN/OUT :
            ZWT     , WA      , WT      , WSLAKE  , LFMASS  , RTMASS  , & ! IN/OUT :
            STMASS  , WOOD    , STBLCP  , FASTCP  , PLAI    , PSAI    , & ! IN/OUT :
            CM      , CH      , TAUSS   ,                               & ! IN/OUT :
            GRAIN   , GDD     , PGS     ,                               & ! IN/OUT
            SMCWTD  ,DEEPRECH , RECH    ,                               & ! IN/OUT :
            Z0WRF   ,                                                   &
            FSA     , FSR     , FIRA    , FSH     , SSOIL   , FCEV    , & ! OUT :
            FGEV    , FCTR    , ECAN    , ETRAN   , ESOIL   , TRAD    , & ! OUT :
            TGB     , TGV     , T2MV    , T2MB    , Q2MV    , Q2MB    , & ! OUT :
            RUNSF   , RUNSB   , APAR    , PSN     , SAV     , SAG     , & ! OUT :
            FSNO    , NEE     , GPP     , NPP     , FVEGMP  , SALB    , & ! OUT :
            QSNBOT  , PONDING , PONDING1, PONDING2, RSSUN   , RSSHA   , & ! OUT :
            ALBSND  , ALBSNI  ,                                         & ! OUT :
            BGAP    , WGAP    , CHV     , CHB     , EMISSI  ,           & ! OUT :
            SHG     , SHC     , SHB     , EVG     , EVB     , GHV     , & ! OUT :
	    GHB     , IRG     , IRC     , IRB     , TR      , EVC     , & ! OUT :
	    CHLEAF  , CHUC    , CHV2    , CHB2    , FPICE   , PAHV    , &
            PAHG    , PAHB    , PAH     , LAISUN  , LAISHA  , RB        &
#ifdef WRF_HYDRO
            , sfcheadrt(i,j)                               &
#endif
            )            ! OUT :

            QFX(I,J) = ECAN + ESOIL + ETRAN
            LH       (I,J)                = FCEV + FGEV + FCTR

   ENDIF ! glacial split ends

#ifdef WRF_HYDRO

!---LPR Attention: lpr added this part 2013-09-04 below to avoid restart NaN values----
!----------------- which may cause the RAPID crash--------------------------
!yw    if(isnan(RUNSF)) then
     if((RUNSF+1.0) .eq. RUNSF) then
         RUNSF   = 0
         RUNSB   = 0
     endif

!----lpr add end------------------------------------------------------------------------

!AD_CHANGE: Glacier cells can produce small negative subsurface runoff for mass balance.
!           This will crash channel routing, so only pass along positive runoff.

     if((RUNSF+1.0) .eq. RUNSF) then
         RUNSF   = 0
         RUNSB   = 0
#ifdef HYDRO_D
        print*, "Warning: NAN found from RUNSF and set to 0."
#endif
     endif

            soldrain(i,j) = max(RUNSB*dt, 0.)        !mm , underground runoff
            INFXSRT(i,j) = RUNSF*dt        !mm , surface runoff

#endif


! INPUT/OUTPUT

             TSK      (I,J)                = TRAD
             HFX      (I,J)                = FSH
             GRDFLX   (I,J)                = SSOIL
	     SMSTAV   (I,J)                = 0.0  ! [maintained as Noah consistency]
             SMSTOT   (I,J)                = 0.0  ! [maintained as Noah consistency]
             SFCRUNOFF(I,J)                = SFCRUNOFF(I,J) + RUNSF * DT
             UDRUNOFF (I,J)                = UDRUNOFF(I,J)  + RUNSB * DT
#ifdef WRF_HYDRO
             if(present(ACCPRCP)) ACCPRCP  (I,J)                = ACCPRCP(I,J) + PRCP * DT
             if(present(ACCECAN)) ACCECAN  (I,J)                = ACCECAN(I,J) + ECAN * DT
             if(present(ACCETRAN)) ACCETRAN (I,J)                = ACCETRAN(I,J) + ETRAN * DT
             if(present(ACCEDIR)) ACCEDIR  (I,J)                = ACCEDIR(I,J) + ESOIL * DT

	     ! Top 2 layer soil saturation
             if(present(SOILSAT_TOP)) SOILSAT_TOP (I,J)           = SUM(SMC(1:2) * DZS(1:2)) / SUM(parameters%smcmax(1:2) * DZS(1:2))
	     ! Soil column integrated soil saturation
             if(present(SOILSAT)) SOILSAT (I,J)  = SUM(SMC(1:NSOIL) * DZS(1:NSOIL)) / SUM(parameters%smcmax(1:NSOIL) * DZS(1:NSOIL))
	     ! Soil column integrated soil ice fraction
	     if(present(SOILICE)) SOILICE (I,J)  = 1. - ( SUM(SMH2O(1:NSOIL) * DZS(1:NSOIL)) / SUM(SMC(1:NSOIL) * DZS(1:NSOIL)) )
             ! Snowpack average temperature (weighted by layer mass)
             if (present(SNOWT_AVG)) then
		if ( SUM(SNICE(-NSNOW+1:0) + SNLIQ(-NSNOW+1:0)) .gt. 0 ) then
		  SNOWT_AVG (I,J)  = SUM(STC(-NSNOW+1:0) * (SNICE(-NSNOW+1:0) + SNLIQ(-NSNOW+1:0))) / SUM(SNICE(-NSNOW+1:0) + SNLIQ(-NSNOW+1:0))
		else
		  SNOWT_AVG (I,J) = undefined_value
		endif
	     endif
#endif
             IF ( SALB > -999 ) THEN
                ALBEDO(I,J)                = SALB
             ENDIF
             SNOWC    (I,J)                = FSNO
             SMOIS    (I,      1:NSOIL,J)  = SMC   (      1:NSOIL)
             SH2O     (I,      1:NSOIL,J)  = SMH2O (      1:NSOIL)
             TSLB     (I,      1:NSOIL,J)  = STC   (      1:NSOIL)
             SNOW     (I,J)                = SWE
             SNOWH    (I,J)                = SNDPTH
             CANWAT   (I,J)                = CANLIQ + CANICE
             ACSNOW   (I,J)                = ACSNOW(I,J) + PRCP * FPICE
             ACSNOM   (I,J)                = ACSNOM(I,J) + QSNBOT*DT + PONDING + PONDING1 + PONDING2
             EMISS    (I,J)                = EMISSI
             QSFC     (I,J)                = QSFC1D

             ISNOWXY  (I,J)                = ISNOW
             TVXY     (I,J)                = TV
             TGXY     (I,J)                = TG
             CANLIQXY (I,J)                = CANLIQ
             CANICEXY (I,J)                = CANICE
             EAHXY    (I,J)                = EAH
             TAHXY    (I,J)                = TAH
             CMXY     (I,J)                = CM
             CHXY     (I,J)                = CH
             FWETXY   (I,J)                = FWET
             SNEQVOXY (I,J)                = SNEQVO
             ALBOLDXY (I,J)                = ALBOLD
             QSNOWXY  (I,J)                = QSNOW
             QRAINXY  (I,J)                = QRAIN
             WSLAKEXY (I,J)                = WSLAKE
             ZWTXY    (I,J)                = ZWT
             WAXY     (I,J)                = WA
             WTXY     (I,J)                = WT
             TSNOXY   (I,-NSNOW+1:    0,J) = STC   (-NSNOW+1:    0)
             ZSNSOXY  (I,-NSNOW+1:NSOIL,J) = ZSNSO (-NSNOW+1:NSOIL)
             SNICEXY  (I,-NSNOW+1:    0,J) = SNICE (-NSNOW+1:    0)
             SNLIQXY  (I,-NSNOW+1:    0,J) = SNLIQ (-NSNOW+1:    0)
             LFMASSXY (I,J)                = LFMASS
             RTMASSXY (I,J)                = RTMASS
             STMASSXY (I,J)                = STMASS
             WOODXY   (I,J)                = WOOD
             STBLCPXY (I,J)                = STBLCP
             FASTCPXY (I,J)                = FASTCP
             XLAIXY   (I,J)                = PLAI
             XSAIXY   (I,J)                = PSAI
             TAUSSXY  (I,J)                = TAUSS

! OUTPUT

             Z0       (I,J)                = Z0WRF
             ZNT      (I,J)                = Z0WRF
             T2MVXY   (I,J)                = T2MV
             T2MBXY   (I,J)                = T2MB
             Q2MVXY   (I,J)                = Q2MV/(1.0 - Q2MV)  ! specific humidity to mixing ratio
             Q2MBXY   (I,J)                = Q2MB/(1.0 - Q2MB)  ! consistent with registry def of Q2
             TRADXY   (I,J)                = TRAD
             NEEXY    (I,J)                = NEE
             GPPXY    (I,J)                = GPP
             NPPXY    (I,J)                = NPP
             FVEGXY   (I,J)                = FVEGMP
             RUNSFXY  (I,J)                = RUNSF
             RUNSBXY  (I,J)                = RUNSB
             ECANXY   (I,J)                = ECAN
             EDIRXY   (I,J)                = ESOIL
             ETRANXY  (I,J)                = ETRAN
             FSAXY    (I,J)                = FSA
             FIRAXY   (I,J)                = FIRA
             APARXY   (I,J)                = APAR
             PSNXY    (I,J)                = PSN
             SAVXY    (I,J)                = SAV
             SAGXY    (I,J)                = SAG
             RSSUNXY  (I,J)                = RSSUN
             RSSHAXY  (I,J)                = RSSHA
             IF(SWE > 0.0) THEN
               ALBSNDXY (I,:,J)            = ALBSND   !snow albedo (direct)
               ALBSNIXY (I,:,J)            = ALBSNI   !snow albedo (diffuse)
             ELSE
               ALBSNDXY (I,:,J)            = 0.0   !snow albedo (direct)
               ALBSNIXY (I,:,J)            = 0.0   !snow albedo (diffuse)
             ENDIF
             BGAPXY   (I,J)                = BGAP
             WGAPXY   (I,J)                = WGAP
             TGVXY    (I,J)                = TGV
             TGBXY    (I,J)                = TGB
             CHVXY    (I,J)                = CHV
             CHBXY    (I,J)                = CHB
             IRCXY    (I,J)                = IRC
             IRGXY    (I,J)                = IRG
             SHCXY    (I,J)                = SHC
             SHGXY    (I,J)                = SHG
             EVGXY    (I,J)                = EVG
             GHVXY    (I,J)                = GHV
             IRBXY    (I,J)                = IRB
             SHBXY    (I,J)                = SHB
             EVBXY    (I,J)                = EVB
             GHBXY    (I,J)                = GHB
             TRXY     (I,J)                = TR
             EVCXY    (I,J)                = EVC
             CHLEAFXY (I,J)                = CHLEAF
             CHUCXY   (I,J)                = CHUC
             CHV2XY   (I,J)                = CHV2
             CHB2XY   (I,J)                = CHB2
             RECHXY   (I,J)                = RECHXY(I,J) + RECH*1.E3 !RECHARGE TO THE WATER TABLE
             DEEPRECHXY(I,J)               = DEEPRECHXY(I,J) + DEEPRECH
             SMCWTDXY(I,J)                 = SMCWTD

             if (crocus_opt /= 0) then
                PSNOWSWEXY(I,1:act_level,J)   = ZP_SNOWSWE(1,1:act_level)
                PSNOWRHOXY(I,1:act_level,J)   = ZP_SNOWRHO(1,1:act_level)
                PSNOWHEATXY(I,1:act_level,J)  = ZP_SNOWHEAT(1,1:act_level)
                PSNOWGRAN1XY(I,1:act_level,J) = ZP_SNOWGRAN1(1,1:act_level)
                PSNOWGRAN2XY(I,1:act_level,J) = ZP_SNOWGRAN2(1,1:act_level)
                PSNOWHISTXY(I,1:act_level,J)  = ZP_SNOWHIST(1,1:act_level)
                PSNOWAGEXY(I,1:act_level,J)   = ZP_SNOWAGE(1,1:act_level)
                PSNOWLIQXY(I,1:act_level,J)   = ZP_SNOWLIQ(1,1:act_level)
                PSNOWTEMPXY(I,1:act_level,J)  = ZP_SNOWTEMP(1,1:act_level)
                PSNOWDZXY(I,1:act_level,J)    = ZP_SNOWDZ(1,1:act_level)

                PSNOWALBXY(I,J)     = ZP_SNOWALB(1)
                PSNOWTHRUFALXY(I,J) = ZP_THRUFAL(1)*DT
                FLOW_ICE(I,J)       = FLOW_ICE(I,J)+ZFLOW_ICE(1)*DT
                FLOW_SNOW(I,J)      = FLOW_SNOW(I,J)+ZFLOW_SNOW(1)*DT
                PSNOWHEIGHTXY(I,J)  = sum(PSNOWDZXY(I,:,J))
                PSNOWTOTSWEXY(I,J)  = sum(PSNOWSWEXY(I,:,J))

                if (GLACINFO(I,J).eq.1 .and. PSNOWHEIGHTXY(I,J).gt.0) then
                   ALBEDO   (I,J) = ZP_SNOWALB(1)
                   SNOW     (I,J) = PSNOWTOTSWEXY(I,J)
                   SNOWH    (I,J) = PSNOWHEIGHTXY(I,J)
                   ACSNOM   (I,J) = FLOW_ICE(I,J)+FLOW_SNOW(I,J)
                   CANWAT   (I,J) = 0. 
#ifdef WRF_HYDRO
                   if(present(ACCETRAN)) ACCETRAN(I,J) = 0.
                   if(present(ACCECAN))  ACCECAN(I,J)  = 0.
                   if(present(ACCEDIR))  ACCEDIR(I,J)  = ACCEDIR(I,J) + (zp_evap(1)*DT) - (ESOIL*DT)
                   if(present(SOILSAT_TOP)) SOILSAT_TOP (I,J) = undefined_value
                   if(present(SOILSAT))   SOILSAT(I,J) = undefined_value
#endif
                endif
             end if ! crocus_opt /= 0
          ENDIF                                                         ! endif of land-sea test
      ENDDO ILOOP                                                       ! of I loop
   ENDDO JLOOP                                                          ! of J loop

!------------------------------------------------------
  END SUBROUTINE noahmplsm
!------------------------------------------------------

SUBROUTINE TRANSFER_MP_PARAMETERS(VEGTYPE,SOILTYPE,SLOPETYPE,SOILCOLOR,CROPTYPE,parameters,iopt_imperv)

  USE NOAHMP_TABLES
  USE MODULE_SF_NOAHMPLSM

  implicit none

  INTEGER, INTENT(IN)    :: VEGTYPE
  INTEGER, INTENT(IN)    :: SOILTYPE
  INTEGER, INTENT(IN)    :: SLOPETYPE
  INTEGER, INTENT(IN)    :: SOILCOLOR
  INTEGER, INTENT(IN)    :: CROPTYPE
  INTEGER, INTENT(IN)    :: iopt_imperv

  type (noahmp_parameters), intent(inout) :: parameters

  REAL    :: REFDK
  REAL    :: REFKDT
  REAL    :: FRZK
  REAL    :: FRZFACT

  parameters%ISWATER   =   ISWATER_TABLE
  parameters%ISBARREN  =  ISBARREN_TABLE
  parameters%ISICE     =     ISICE_TABLE
  parameters%ISCROP    =    ISCROP_TABLE
  parameters%EBLFOREST = EBLFOREST_TABLE

  parameters%URBAN_FLAG = .FALSE.
  IF( VEGTYPE == ISURBAN_TABLE                  .or. VEGTYPE == LOW_DENSITY_RESIDENTIAL_TABLE  .or. &
      VEGTYPE == HIGH_DENSITY_RESIDENTIAL_TABLE .or. VEGTYPE == HIGH_INTENSITY_INDUSTRIAL_TABLE ) THEN
     parameters%URBAN_FLAG = .TRUE.
  ENDIF

!------------------------------------------------------------------------------------------!
! Transfer veg parameters
!------------------------------------------------------------------------------------------!

  parameters%CH2OP  =  CH2OP_TABLE(VEGTYPE)       !maximum intercepted h2o per unit lai+sai (mm)
  parameters%DLEAF  =  DLEAF_TABLE(VEGTYPE)       !characteristic leaf dimension (m)
  parameters%Z0MVT  =  Z0MVT_TABLE(VEGTYPE)       !momentum roughness length (m)
!  parameters%HVT    =    HVT_TABLE(VEGTYPE)       !top of canopy (m)
  parameters%HVB    =    HVB_TABLE(VEGTYPE)       !bottom of canopy (m)
  parameters%DEN    =    DEN_TABLE(VEGTYPE)       !tree density (no. of trunks per m2)
  parameters%RC     =     RC_TABLE(VEGTYPE)       !tree crown radius (m)
!  parameters%MFSNO  =  MFSNO_TABLE(VEGTYPE)       !snowmelt m parameter ()
  parameters%SAIM   =   SAIM_TABLE(VEGTYPE,:)     !monthly stem area index, one-sided
  parameters%LAIM   =   LAIM_TABLE(VEGTYPE,:)     !monthly leaf area index, one-sided
  parameters%SLA    =    SLA_TABLE(VEGTYPE)       !single-side leaf area per Kg [m2/kg]
  parameters%DILEFC = DILEFC_TABLE(VEGTYPE)       !coeficient for leaf stress death [1/s]
  parameters%DILEFW = DILEFW_TABLE(VEGTYPE)       !coeficient for leaf stress death [1/s]
  parameters%FRAGR  =  FRAGR_TABLE(VEGTYPE)       !fraction of growth respiration  !original was 0.3
  parameters%LTOVRC = LTOVRC_TABLE(VEGTYPE)       !leaf turnover [1/s]

  parameters%C3PSN  =  C3PSN_TABLE(VEGTYPE)       !photosynthetic pathway: 0. = c4, 1. = c3
  parameters%KC25   =   KC25_TABLE(VEGTYPE)       !co2 michaelis-menten constant at 25c (pa)
  parameters%AKC    =    AKC_TABLE(VEGTYPE)       !q10 for kc25
  parameters%KO25   =   KO25_TABLE(VEGTYPE)       !o2 michaelis-menten constant at 25c (pa)
  parameters%AKO    =    AKO_TABLE(VEGTYPE)       !q10 for ko25
!  parameters%VCMX25 = VCMX25_TABLE(VEGTYPE)       !maximum rate of carboxylation at 25c (umol co2/m**2/s)
  parameters%AVCMX  =  AVCMX_TABLE(VEGTYPE)       !q10 for vcmx25
  parameters%BP     =     BP_TABLE(VEGTYPE)       !minimum leaf conductance (umol/m**2/s)
!  parameters%MP     =     MP_TABLE(VEGTYPE)       !slope of conductance-to-photosynthesis relationship
  parameters%QE25   =   QE25_TABLE(VEGTYPE)       !quantum efficiency at 25c (umol co2 / umol photon)
  parameters%AQE    =    AQE_TABLE(VEGTYPE)       !q10 for qe25
  parameters%RMF25  =  RMF25_TABLE(VEGTYPE)       !leaf maintenance respiration at 25c (umol co2/m**2/s)
  parameters%RMS25  =  RMS25_TABLE(VEGTYPE)       !stem maintenance respiration at 25c (umol co2/kg bio/s)
  parameters%RMR25  =  RMR25_TABLE(VEGTYPE)       !root maintenance respiration at 25c (umol co2/kg bio/s)
  parameters%ARM    =    ARM_TABLE(VEGTYPE)       !q10 for maintenance respiration
  parameters%FOLNMX = FOLNMX_TABLE(VEGTYPE)       !foliage nitrogen concentration when f(n)=1 (%)
  parameters%TMIN   =   TMIN_TABLE(VEGTYPE)       !minimum temperature for photosynthesis (k)

  parameters%XL     =     XL_TABLE(VEGTYPE)       !leaf/stem orientation index
  parameters%RHOL   =   RHOL_TABLE(VEGTYPE,:)     !leaf reflectance: 1=vis, 2=nir
  parameters%RHOS   =   RHOS_TABLE(VEGTYPE,:)     !stem reflectance: 1=vis, 2=nir
  parameters%TAUL   =   TAUL_TABLE(VEGTYPE,:)     !leaf transmittance: 1=vis, 2=nir
  parameters%TAUS   =   TAUS_TABLE(VEGTYPE,:)     !stem transmittance: 1=vis, 2=nir

  parameters%MRP    =    MRP_TABLE(VEGTYPE)       !microbial respiration parameter (umol co2 /kg c/ s)
!  parameters%CWPVT  =  CWPVT_TABLE(VEGTYPE)       !empirical canopy wind parameter

  parameters%WRRAT  =  WRRAT_TABLE(VEGTYPE)       !wood to non-wood ratio
  parameters%WDPOOL = WDPOOL_TABLE(VEGTYPE)       !wood pool (switch 1 or 0) depending on woody or not [-]
  parameters%TDLEF  =  TDLEF_TABLE(VEGTYPE)       !characteristic T for leaf freezing [K]

  parameters%NROOT  =  NROOT_TABLE(VEGTYPE)       !number of soil layers with root present
  parameters%RGL    =    RGL_TABLE(VEGTYPE)       !Parameter used in radiation stress function
  parameters%RSMIN  =     RS_TABLE(VEGTYPE)       !Minimum stomatal resistance [s m-1]
  parameters%HS     =     HS_TABLE(VEGTYPE)       !Parameter used in vapor pressure deficit function
  parameters%TOPT   =   TOPT_TABLE(VEGTYPE)       !Optimum transpiration air temperature [K]
  parameters%RSMAX  =  RSMAX_TABLE(VEGTYPE)       !Maximal stomatal resistance [s m-1]

!------------------------------------------------------------------------------------------!
! Transfer rad parameters
!------------------------------------------------------------------------------------------!

   parameters%ALBSAT    = ALBSAT_TABLE(SOILCOLOR,:)
   parameters%ALBDRY    = ALBDRY_TABLE(SOILCOLOR,:)
   parameters%ALBICE    = ALBICE_TABLE
   parameters%ALBLAK    = ALBLAK_TABLE
   parameters%OMEGAS    = OMEGAS_TABLE
   parameters%BETADS    = BETADS_TABLE
   parameters%BETAIS    = BETAIS_TABLE
   parameters%EG        = EG_TABLE

!------------------------------------------------------------------------------------------!
! Transfer crop parameters
!------------------------------------------------------------------------------------------!

  IF(CROPTYPE > 0) THEN
   parameters%PLTDAY    =    PLTDAY_TABLE(CROPTYPE)    ! Planting date
   parameters%HSDAY     =     HSDAY_TABLE(CROPTYPE)    ! Harvest date
   parameters%PLANTPOP  =  PLANTPOP_TABLE(CROPTYPE)    ! Plant density [per ha] - used?
   parameters%IRRI      =      IRRI_TABLE(CROPTYPE)    ! Irrigation strategy 0= non-irrigation 1=irrigation (no water-stress)
   parameters%GDDTBASE  =  GDDTBASE_TABLE(CROPTYPE)    ! Base temperature for GDD accumulation [C]
   parameters%GDDTCUT   =   GDDTCUT_TABLE(CROPTYPE)    ! Upper temperature for GDD accumulation [C]
   parameters%GDDS1     =     GDDS1_TABLE(CROPTYPE)    ! GDD from seeding to emergence
   parameters%GDDS2     =     GDDS2_TABLE(CROPTYPE)    ! GDD from seeding to initial vegetative
   parameters%GDDS3     =     GDDS3_TABLE(CROPTYPE)    ! GDD from seeding to post vegetative
   parameters%GDDS4     =     GDDS4_TABLE(CROPTYPE)    ! GDD from seeding to intial reproductive
   parameters%GDDS5     =     GDDS5_TABLE(CROPTYPE)    ! GDD from seeding to pysical maturity
   parameters%C3C4      =      C3C4_TABLE(CROPTYPE)    ! photosynthetic pathway:  1. = c3 2. = c4
   parameters%AREF      =      AREF_TABLE(CROPTYPE)    ! reference maximum CO2 assimulation rate
   parameters%PSNRF     =     PSNRF_TABLE(CROPTYPE)    ! CO2 assimulation reduction factor(0-1) (caused by non-modeling part,e.g.pest,weeds)
   parameters%I2PAR     =     I2PAR_TABLE(CROPTYPE)    ! Fraction of incoming solar radiation to photosynthetically active radiation
   parameters%TASSIM0   =   TASSIM0_TABLE(CROPTYPE)    ! Minimum temperature for CO2 assimulation [C]
   parameters%TASSIM1   =   TASSIM1_TABLE(CROPTYPE)    ! CO2 assimulation linearly increasing until temperature reaches T1 [C]
   parameters%TASSIM2   =   TASSIM2_TABLE(CROPTYPE)    ! CO2 assmilation rate remain at Aref until temperature reaches T2 [C]
   parameters%K         =         K_TABLE(CROPTYPE)    ! light extinction coefficient
   parameters%EPSI      =      EPSI_TABLE(CROPTYPE)    ! initial light use efficiency
   parameters%Q10MR     =     Q10MR_TABLE(CROPTYPE)    ! q10 for maintainance respiration
   parameters%FOLN_MX   =   FOLN_MX_TABLE(CROPTYPE)    ! foliage nitrogen concentration when f(n)=1 (%)
   parameters%LEFREEZ   =   LEFREEZ_TABLE(CROPTYPE)    ! characteristic T for leaf freezing [K]
   parameters%DILE_FC   =   DILE_FC_TABLE(CROPTYPE,:)  ! coeficient for temperature leaf stress death [1/s]
   parameters%DILE_FW   =   DILE_FW_TABLE(CROPTYPE,:)  ! coeficient for water leaf stress death [1/s]
   parameters%FRA_GR    =    FRA_GR_TABLE(CROPTYPE)    ! fraction of growth respiration
   parameters%LF_OVRC   =   LF_OVRC_TABLE(CROPTYPE,:)  ! fraction of leaf turnover  [1/s]
   parameters%ST_OVRC   =   ST_OVRC_TABLE(CROPTYPE,:)  ! fraction of stem turnover  [1/s]
   parameters%RT_OVRC   =   RT_OVRC_TABLE(CROPTYPE,:)  ! fraction of root tunrover  [1/s]
   parameters%LFMR25    =    LFMR25_TABLE(CROPTYPE)    ! leaf maintenance respiration at 25C [umol CO2/m**2  /s]
   parameters%STMR25    =    STMR25_TABLE(CROPTYPE)    ! stem maintenance respiration at 25C [umol CO2/kg bio/s]
   parameters%RTMR25    =    RTMR25_TABLE(CROPTYPE)    ! root maintenance respiration at 25C [umol CO2/kg bio/s]
   parameters%GRAINMR25 = GRAINMR25_TABLE(CROPTYPE)    ! grain maintenance respiration at 25C [umol CO2/kg bio/s]
   parameters%LFPT      =      LFPT_TABLE(CROPTYPE,:)  ! fraction of carbohydrate flux to leaf
   parameters%STPT      =      STPT_TABLE(CROPTYPE,:)  ! fraction of carbohydrate flux to stem
   parameters%RTPT      =      RTPT_TABLE(CROPTYPE,:)  ! fraction of carbohydrate flux to root
   parameters%GRAINPT   =   GRAINPT_TABLE(CROPTYPE,:)  ! fraction of carbohydrate flux to grain
   parameters%BIO2LAI   =   BIO2LAI_TABLE(CROPTYPE)    ! leaf are per living leaf biomass [m^2/kg]
  END IF

!------------------------------------------------------------------------------------------!
! Transfer global parameters
!------------------------------------------------------------------------------------------!

   parameters%CO2        =         CO2_TABLE
   parameters%O2         =          O2_TABLE
   parameters%TIMEAN     =      TIMEAN_TABLE
   parameters%FSATMX     =      FSATMX_TABLE
   parameters%Z0SNO      =       Z0SNO_TABLE
   parameters%SSI        =         SSI_TABLE
   parameters%SNOW_RET_FAC =  SNOW_RET_FAC_TABLE
   parameters%SWEMX      =       SWEMX_TABLE
   parameters%GRAIN_GROWTH = GRAIN_GROWTH_TABLE
   parameters%EXTRA_GROWTH = EXTRA_GROWTH_TABLE
   parameters%DIRT_SOOT    =    DIRT_SOOT_TABLE
   parameters%BATS_COSZ    =    BATS_COSZ_TABLE
   parameters%BATS_VIS_NEW = BATS_VIS_NEW_TABLE
   parameters%BATS_NIR_NEW = BATS_NIR_NEW_TABLE
   parameters%BATS_VIS_AGE = BATS_VIS_AGE_TABLE
   parameters%BATS_NIR_AGE = BATS_NIR_AGE_TABLE
   parameters%BATS_VIS_DIR = BATS_VIS_DIR_TABLE
   parameters%BATS_NIR_DIR = BATS_NIR_DIR_TABLE
   parameters%RSURF_SNOW =  RSURF_SNOW_TABLE
   parameters%SWE_LIMIT  =  SWE_LIMIT_TABLE

! ----------------------------------------------------------------------
!  Transfer soil parameters
! ----------------------------------------------------------------------

#ifndef SPATIAL_SOIL
    parameters%BEXP   = BEXP_TABLE   (SOILTYPE)
    parameters%DKSAT  = DKSAT_TABLE  (SOILTYPE)
    parameters%DWSAT  = DWSAT_TABLE  (SOILTYPE)
    parameters%F1     = F1_TABLE     (SOILTYPE)
    parameters%PSISAT = PSISAT_TABLE (SOILTYPE)
    parameters%QUARTZ = QUARTZ_TABLE (SOILTYPE)
    parameters%SMCDRY = SMCDRY_TABLE (SOILTYPE)
    parameters%SMCMAX = SMCMAX_TABLE (SOILTYPE)
    parameters%SMCREF = SMCREF_TABLE (SOILTYPE)
    parameters%SMCWLT = SMCWLT_TABLE (SOILTYPE)
    parameters%REFDK  = REFDK_TABLE
    parameters%REFKDT = REFKDT_TABLE
    parameters%SLOPE  = SLOPE_TABLE(SLOPETYPE)
    parameters%CWPVT  =  CWPVT_TABLE(VEGTYPE)       !empirical canopy wind parameter
    parameters%VCMX25 = VCMX25_TABLE(VEGTYPE)       !maximum rate of carboxylation at 25c (umol co2/m**2/s)
    parameters%MFSNO  =  MFSNO_TABLE(VEGTYPE)       !snowmelt m parameter ()
    parameters%HVT    =    HVT_TABLE(VEGTYPE)       !top of canopy (m)
    parameters%MP     =     MP_TABLE(VEGTYPE)       !slope of conductance-to-photosynthesis relationship
    parameters%RSURF_EXP  =   RSURF_EXP_TABLE
    parameters%AXAJ = AXAJ_TABLE(SOILTYPE)
    parameters%BXAJ = BXAJ_TABLE(SOILTYPE)
    parameters%XXAJ = XXAJ_TABLE(SOILTYPE)
    IF (parameters%URBAN_FLAG) THEN
      parameters%IMPERV = IMPERV_URBAN_TABLE
    ELSE
      parameters%IMPERV = 0.0
    ENDIF
    parameters%SSI    = SSI_TABLE
    parameters%SNOW_RET_FAC = SNOW_RET_FAC_TABLE
    parameters%TAU0   = TAU0_TABLE
    parameters%RSURF_SNOW = RSURF_SNOW_TABLE
    parameters%SCAMAX = SCAMAX_TABLE
#endif
! ----------------------------------------------------------------------
! Transfer GENPARM parameters
! ----------------------------------------------------------------------
    parameters%CSOIL  = CSOIL_TABLE
    parameters%ZBOT   = ZBOT_TABLE
    parameters%CZIL   = CZIL_TABLE

    FRZK   = FRZK_TABLE
    parameters%KDT    = parameters%REFKDT * parameters%DKSAT(1) / parameters%REFDK

    IF(parameters%URBAN_FLAG)THEN  ! Hardcoding some urban parameters for soil
       if (iopt_imperv .eq. 9) then
         parameters%SMCMAX = 0.45
         parameters%SMCREF = 0.42
         parameters%SMCWLT = 0.40
         parameters%SMCDRY = 0.40
       endif
       parameters%CSOIL  = 3.E6
    ENDIF

! adjust FRZK parameter to actual soil type: FRZK * FRZFACT

    IF(SOILTYPE /= 14) then
      FRZFACT = (parameters%SMCMAX(1) / parameters%SMCREF(1)) * (0.412 / 0.468)
      parameters%FRZX = FRZK * FRZFACT
    END IF

 END SUBROUTINE TRANSFER_MP_PARAMETERS

SUBROUTINE PEDOTRANSFER_SR2006(nsoil,sand,clay,orgm,parameters)

  use module_sf_noahmplsm
  use noahmp_tables

  implicit none

  integer,                    intent(in   ) :: nsoil     ! number of soil layers
  real, dimension( 1:nsoil ), intent(inout) :: sand
  real, dimension( 1:nsoil ), intent(inout) :: clay
  real, dimension( 1:nsoil ), intent(inout) :: orgm

  real, dimension( 1:nsoil ) :: theta_1500t
  real, dimension( 1:nsoil ) :: theta_1500
  real, dimension( 1:nsoil ) :: theta_33t
  real, dimension( 1:nsoil ) :: theta_33
  real, dimension( 1:nsoil ) :: theta_s33t
  real, dimension( 1:nsoil ) :: theta_s33
  real, dimension( 1:nsoil ) :: psi_et
  real, dimension( 1:nsoil ) :: psi_e

  type(noahmp_parameters), intent(inout) :: parameters
  integer :: k

  do k = 1,4
    if(sand(k) <= 0 .or. clay(k) <= 0) then
      sand(k) = 0.41
      clay(k) = 0.18
    end if
    if(orgm(k) <= 0 ) orgm(k) = 0.0
  end do

  theta_1500t =   sr2006_theta_1500t_a*sand       &
                + sr2006_theta_1500t_b*clay       &
                + sr2006_theta_1500t_c*orgm       &
                + sr2006_theta_1500t_d*sand*orgm  &
                + sr2006_theta_1500t_e*clay*orgm  &
                + sr2006_theta_1500t_f*sand*clay  &
                + sr2006_theta_1500t_g

  theta_1500  =   theta_1500t                      &
                + sr2006_theta_1500_a*theta_1500t  &
                + sr2006_theta_1500_b

  theta_33t   =   sr2006_theta_33t_a*sand       &
                + sr2006_theta_33t_b*clay       &
                + sr2006_theta_33t_c*orgm       &
                + sr2006_theta_33t_d*sand*orgm  &
                + sr2006_theta_33t_e*clay*orgm  &
                + sr2006_theta_33t_f*sand*clay  &
                + sr2006_theta_33t_g

  theta_33    =   theta_33t                              &
                + sr2006_theta_33_a*theta_33t*theta_33t  &
                + sr2006_theta_33_b*theta_33t            &
                + sr2006_theta_33_c

  theta_s33t  =   sr2006_theta_s33t_a*sand      &
                + sr2006_theta_s33t_b*clay      &
                + sr2006_theta_s33t_c*orgm      &
                + sr2006_theta_s33t_d*sand*orgm &
                + sr2006_theta_s33t_e*clay*orgm &
                + sr2006_theta_s33t_f*sand*clay &
                + sr2006_theta_s33t_g

  theta_s33   = theta_s33t                       &
                + sr2006_theta_s33_a*theta_s33t  &
                + sr2006_theta_s33_b

  psi_et      =   sr2006_psi_et_a*sand           &
                + sr2006_psi_et_b*clay           &
                + sr2006_psi_et_c*theta_s33      &
                + sr2006_psi_et_d*sand*theta_s33 &
                + sr2006_psi_et_e*clay*theta_s33 &
                + sr2006_psi_et_f*sand*clay      &
                + sr2006_psi_et_g

  psi_e       =   psi_et                        &
                + sr2006_psi_e_a*psi_et*psi_et  &
                + sr2006_psi_e_b*psi_et         &
                + sr2006_psi_e_c

  parameters%smcwlt = theta_1500
  parameters%smcref = theta_33
  parameters%smcmax =   theta_33    &
                      + theta_s33            &
                      + sr2006_smcmax_a*sand &
                      + sr2006_smcmax_b

  parameters%bexp   = 3.816712826 / (log(theta_33) - log(theta_1500) )
  parameters%psisat = psi_e
  parameters%dksat  = 1930.0 * (parameters%smcmax - theta_33) ** (3.0 - 1.0/parameters%bexp)
  parameters%quartz = sand

! Units conversion

  parameters%psisat = max(0.1,parameters%psisat)     ! arbitrarily impose a limit of 0.1kpa
  parameters%psisat = 0.101997 * parameters%psisat   ! convert kpa to m
  parameters%dksat  = parameters%dksat / 3600000.0   ! convert mm/h to m/s
  parameters%dwsat  = parameters%dksat * parameters%psisat *parameters%bexp / parameters%smcmax  ! units should be m*m/s
  parameters%smcdry = parameters%smcwlt

! Introducing somewhat arbitrary limits (based on SOILPARM) to prevent bad things

  parameters%smcmax = max(0.32 ,min(parameters%smcmax,             0.50 ))
  parameters%smcref = max(0.17 ,min(parameters%smcref,parameters%smcmax ))
  parameters%smcwlt = max(0.01 ,min(parameters%smcwlt,parameters%smcref ))
  parameters%smcdry = max(0.01 ,min(parameters%smcdry,parameters%smcref ))
  parameters%bexp   = max(2.50 ,min(parameters%bexp,               12.0 ))
  parameters%psisat = max(0.03 ,min(parameters%psisat,             1.00 ))
  parameters%dksat  = max(5.e-7,min(parameters%dksat,              1.e-5))
  parameters%dwsat  = max(1.e-6,min(parameters%dwsat,              3.e-5))
  parameters%quartz = max(0.05 ,min(parameters%quartz,             0.95 ))

 END SUBROUTINE PEDOTRANSFER_SR2006

  SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP ,   IVGTYP, &
       TSLB , SMOIS , SH2O , DZS , FNDSOILW , FNDSNOWH ,             &
       TSK, isnowxy , tvxy     ,tgxy     ,canicexy ,         TMN,     XICE,   &
       canliqxy ,eahxy    ,tahxy    ,cmxy     ,chxy     ,                     &
       fwetxy   ,sneqvoxy ,alboldxy ,qsnowxy  ,wslakexy ,zwtxy    ,waxy     , &
       wtxy     ,tsnoxy   ,zsnsoxy  ,snicexy  ,snliqxy  ,lfmassxy ,rtmassxy , &
       stmassxy ,woodxy   ,stblcpxy ,fastcpxy ,xsaixy   ,lai      ,           &
!jref:start
       t2mvxy   ,t2mbxy   ,chstarxy,            &
!jref:end
       NSOIL, restart,                 &
       allowed_to_read , iopt_run, iopt_crop, iopt_pedo,   &
       ids,ide, jds,jde, kds,kde,                &
       ims,ime, jms,jme, kms,kme,                &
       its,ite, jts,jte, kts,kte,                &
       smoiseq  ,smcwtdxy ,rechxy   ,deeprechxy, areaxy, dx, dy, msftx, msfty,&     ! Optional groundwater
       wtddt    ,stepwtd  ,dt       ,qrfsxy     ,qspringsxy  , qslatxy    ,  &      ! Optional groundwater
       fdepthxy ,ht     ,riverbedxy ,eqzwt     ,rivercondxy ,pexpxy            )    ! Optional groundwater

  USE NOAHMP_TABLES

  IMPLICIT NONE

! Initializing Canopy air temperature to 287 K seems dangerous to me [KWM].

    INTEGER, INTENT(IN   )    ::     ids,ide, jds,jde, kds,kde,  &
         &                           ims,ime, jms,jme, kms,kme,  &
         &                           its,ite, jts,jte, kts,kte

    INTEGER, INTENT(IN)       ::     NSOIL, iopt_run, iopt_crop, iopt_pedo

    LOGICAL, INTENT(IN)       ::     restart,                    &
         &                           allowed_to_read

    REAL,    DIMENSION( NSOIL), INTENT(IN)    ::     DZS  ! Thickness of the soil layers [m]
    REAL,    INTENT(IN) , OPTIONAL ::     DX, DY
    REAL,    DIMENSION( ims:ime, jms:jme ) ,  INTENT(IN) , OPTIONAL :: MSFTX,MSFTY

    REAL,    DIMENSION( ims:ime, NSOIL, jms:jme ) ,    &
         &   INTENT(INOUT)    ::     SMOIS,                      &
         &                           SH2O,                       &
         &                           TSLB

    REAL,    DIMENSION( ims:ime, jms:jme ) ,                     &
         &   INTENT(INOUT)    ::     SNOW,                       &
         &                           SNOWH,                      &
         &                           CANWAT

    INTEGER, DIMENSION( ims:ime, jms:jme ),                      &
         &   INTENT(IN)       ::     ISLTYP,  &
                                     IVGTYP

    LOGICAL, INTENT(IN)       ::     FNDSOILW,                   &
         &                           FNDSNOWH

    REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: TSK         !skin temperature (k)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: TMN         !deep soil temperature (k)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: XICE         !sea ice fraction
    INTEGER, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: isnowxy     !actual no. of snow layers
    REAL, DIMENSION(ims:ime,-2:NSOIL,jms:jme), INTENT(INOUT) :: zsnsoxy  !snow layer depth [m]
    REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: tsnoxy   !snow temperature [K]
    REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: snicexy  !snow layer ice [mm]
    REAL, DIMENSION(ims:ime,-2:              0,jms:jme), INTENT(INOUT) :: snliqxy  !snow layer liquid water [mm]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tvxy        !vegetation canopy temperature
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tgxy        !ground surface temperature
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canicexy    !canopy-intercepted ice (mm)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canliqxy    !canopy-intercepted liquid water (mm)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: eahxy       !canopy air vapor pressure (pa)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tahxy       !canopy air temperature (k)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: cmxy        !momentum drag coefficient
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: chxy        !sensible heat exchange coefficient
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fwetxy      !wetted or snowed fraction of the canopy (-)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: sneqvoxy    !snow mass at last time step(mm h2o)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: alboldxy    !snow albedo at last time step (-)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: qsnowxy     !snowfall on the ground [mm/s]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wslakexy    !lake water storage [mm]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: zwtxy       !water table depth [m]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: waxy        !water in the "aquifer" [mm]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wtxy        !groundwater storage [mm]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: lfmassxy    !leaf mass [g/m2]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: rtmassxy    !mass of fine roots [g/m2]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stmassxy    !stem mass [g/m2]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: woodxy      !mass of wood (incl. woody roots) [g/m2]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stblcpxy    !stable carbon in deep soil [g/m2]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fastcpxy    !short-lived carbon, shallow soil [g/m2]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: xsaixy      !stem area index
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: lai         !leaf area index

! IOPT_RUN = 5 option

    REAL, DIMENSION(ims:ime,1:nsoil,jms:jme), INTENT(INOUT) , OPTIONAL :: smoiseq !equilibrium soil moisture content [m3m-3]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) , OPTIONAL :: smcwtdxy    !deep soil moisture content [m3m-3]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) , OPTIONAL :: deeprechxy  !deep recharge [m]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) , OPTIONAL :: rechxy      !accumulated recharge [mm]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) , OPTIONAL :: qrfsxy      !accumulated flux from groundwater to rivers [mm]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) , OPTIONAL :: qspringsxy  !accumulated seeping water [mm]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) , OPTIONAL :: qslatxy     !accumulated lateral flow [mm]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) , OPTIONAL :: areaxy      !grid cell area [m2]
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) , OPTIONAL :: FDEPTHXY    !efolding depth for transmissivity (m)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) , OPTIONAL :: HT          !terrain height (m)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) , OPTIONAL :: RIVERBEDXY  !riverbed depth (m)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) , OPTIONAL :: EQZWT       !equilibrium water table depth (m)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) , OPTIONAL :: RIVERCONDXY !river conductance
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) , OPTIONAL :: PEXPXY      !factor for river conductance

    INTEGER,  INTENT(INOUT) , OPTIONAL :: STEPWTD
    REAL, INTENT(IN) , OPTIONAL :: DT, WTDDT

!jref:start
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: t2mvxy        !2m temperature vegetation part (k)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: t2mbxy        !2m temperature bare ground part (k)
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: chstarxy        !dummy
!jref:end


    REAL, DIMENSION(1:NSOIL)  :: ZSOIL      ! Depth of the soil layer bottom (m) from
    !                                                   the surface (negative)

    REAL                      :: BEXP, SMCMAX, PSISAT
    REAL                      :: FK, masslai,masssai

    REAL, PARAMETER           :: BLIM  = 5.5
    REAL, PARAMETER           :: HLICE = 3.335E5
    REAL, PARAMETER           :: GRAV = 9.81
    REAL, PARAMETER           :: T0 = 273.15

    INTEGER                   :: errflag, i,j,itf,jtf,ns

    character(len=80) :: err_message
    character(len=4)  :: MMINSL
    character(len=*), intent(in) :: MMINLU
    MMINSL='STAS'

    call read_mp_veg_parameters(trim(MMINLU))
    call read_mp_soil_parameters()
    call read_mp_rad_parameters()
    call read_mp_global_parameters()
    if(iopt_crop == 1) call read_mp_crop_parameters()
    if(iopt_pedo == 1) call read_mp_optional_parameters()

    IF( .NOT. restart ) THEN

       itf=min0(ite,ide-1)
       jtf=min0(jte,jde-1)

       !
       ! initialize physical snow height SNOWH
       !
       IF(.NOT.FNDSNOWH)THEN
          ! If no SNOWH do the following
          CALL wrf_message( 'SNOW HEIGHT NOT FOUND - VALUE DEFINED IN LSMINIT' )
          DO J = jts,jtf
             DO I = its,itf
                SNOWH(I,J)=SNOW(I,J)*0.005               ! SNOW in mm and SNOWH in m
             ENDDO
          ENDDO
       ENDIF

       errflag = 0
       DO j = jts,jtf
          DO i = its,itf
             IF ( ISLTYP( i,j ) .LT. 1 ) THEN
                errflag = 1
                WRITE(err_message,*)"module_sf_noahlsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
                CALL wrf_message(err_message)
             ENDIF
          ENDDO
       ENDDO
       IF ( errflag .EQ. 1 ) THEN
          CALL wrf_error_fatal( "module_sf_noahlsm.F: lsminit: out of range value "// &
               "of ISLTYP. Is this field in the input?" )
       ENDIF
! GAC-->LATERALFLOW
! 20130219 - No longer need this - see module_data_gocart_dust
!#if ( WRF_CHEM == 1 )
!       !
!       ! need this parameter for dust parameterization in wrf/chem
!       !
!       do I=1,NSLTYPE
!          porosity(i)=maxsmc(i)
!       enddo
!#endif
! <--GAC

! initialize soil liquid water content SH2O

       DO J = jts , jtf
          DO I = its , itf
	    IF(IVGTYP(I,J)==ISICE_TABLE .AND. XICE(I,J) <= 0.0) THEN
              DO NS=1, NSOIL
	        SMOIS(I,NS,J) = 1.0                     ! glacier starts all frozen
	        SH2O(I,NS,J) = 0.0
	        TSLB(I,NS,J) = MIN(TSLB(I,NS,J),263.15) ! set glacier temp to at most -10C
              END DO
	        !TMN(I,J) = MIN(TMN(I,J),263.15)         ! set deep temp to at most -10C
		SNOW(I,J) = MAX(SNOW(I,J), 10.0)        ! set SWE to at least 10mm
                SNOWH(I,J)=SNOW(I,J)*0.01               ! SNOW in mm and SNOWH in m
	    ELSE

              BEXP   =   BEXP_TABLE(ISLTYP(I,J))
              SMCMAX = SMCMAX_TABLE(ISLTYP(I,J))
              PSISAT = PSISAT_TABLE(ISLTYP(I,J))

              DO NS=1, NSOIL
	        IF ( SMOIS(I,NS,J) > SMCMAX )  SMOIS(I,NS,J) = SMCMAX
              END DO
              IF ( ( BEXP > 0.0 ) .AND. ( SMCMAX > 0.0 ) .AND. ( PSISAT > 0.0 ) ) THEN
                DO NS=1, NSOIL
                   IF ( TSLB(I,NS,J) < 273.149 ) THEN    ! Use explicit as initial soil ice
                      FK=(( (HLICE/(GRAV*(-PSISAT))) *                              &
                           ((TSLB(I,NS,J)-T0)/TSLB(I,NS,J)) )**(-1/BEXP) )*SMCMAX
                      FK = MAX(FK, 0.02)
                      SH2O(I,NS,J) = MIN( FK, SMOIS(I,NS,J) )
                   ELSE
                      SH2O(I,NS,J)=SMOIS(I,NS,J)
                   ENDIF
                END DO
              ELSE
                DO NS=1, NSOIL
                   SH2O(I,NS,J)=SMOIS(I,NS,J)
                END DO
              ENDIF
            ENDIF
          ENDDO
       ENDDO
!  ENDIF


       DO J = jts,jtf
          DO I = its,itf
             tvxy       (I,J) = TSK(I,J)
	       if(snow(i,j) > 0.0 .and. tsk(i,j) > 273.15) tvxy(I,J) = 273.15
             tgxy       (I,J) = TSK(I,J)
	       if(snow(i,j) > 0.0 .and. tsk(i,j) > 273.15) tgxy(I,J) = 273.15
             CANWAT     (I,J) = 0.0
             canliqxy   (I,J) = CANWAT(I,J)
             canicexy   (I,J) = 0.
             eahxy      (I,J) = 2000.
             tahxy      (I,J) = TSK(I,J)
	       if(snow(i,j) > 0.0 .and. tsk(i,j) > 273.15) tahxy(I,J) = 273.15
!             tahxy      (I,J) = 287.
!jref:start
             t2mvxy     (I,J) = TSK(I,J)
	       if(snow(i,j) > 0.0 .and. tsk(i,j) > 273.15) t2mvxy(I,J) = 273.15
             t2mbxy     (I,J) = TSK(I,J)
	       if(snow(i,j) > 0.0 .and. tsk(i,j) > 273.15) t2mbxy(I,J) = 273.15
             chstarxy     (I,J) = 0.1
!jref:end

             cmxy       (I,J) = 0.0
             chxy       (I,J) = 0.0
             fwetxy     (I,J) = 0.0
             sneqvoxy   (I,J) = 0.0
             alboldxy   (I,J) = 0.65
             qsnowxy    (I,J) = 0.0
             wslakexy   (I,J) = 0.0

             if(iopt_run.ne.5) then
                   waxy       (I,J) = 4900.                                       !???
                   wtxy       (I,J) = waxy(i,j)                                   !???
                   zwtxy      (I,J) = (25. + 2.0) - waxy(i,j)/1000/0.2            !???
             else
                   waxy       (I,J) = 0.
                   wtxy       (I,J) = 0.
                   areaxy     (I,J) = (DX * DY) / ( MSFTX(I,J) * MSFTY(I,J) )
             endif

           IF(IVGTYP(I,J) == ISBARREN_TABLE .OR. IVGTYP(I,J) == ISICE_TABLE .OR. &
	      IVGTYP(I,J) == ISURBAN_TABLE  .OR. IVGTYP(I,J) == ISWATER_TABLE ) THEN

	     lai        (I,J) = 0.0
             xsaixy     (I,J) = 0.0
             lfmassxy   (I,J) = 0.0
             stmassxy   (I,J) = 0.0
             rtmassxy   (I,J) = 0.0
             woodxy     (I,J) = 0.0
             stblcpxy   (I,J) = 0.0
             fastcpxy   (I,J) = 0.0

	   ELSE

	     lai        (I,J) = max(lai(i,j),0.05)             ! at least start with 0.05 for arbitrary initialization (v3.7)
             xsaixy     (I,J) = max(0.1*lai(I,J),0.05)         ! MB: arbitrarily initialize SAI using input LAI (v3.7)
             masslai = 1000. / max(SLA_TABLE(IVGTYP(I,J)),1.0) ! conversion from lai to mass  (v3.7)
             lfmassxy   (I,J) = lai(i,j)*masslai               ! use LAI to initialize (v3.7)
             masssai = 1000. / 3.0                             ! conversion from lai to mass (v3.7)
             stmassxy   (I,J) = xsaixy(i,j)*masssai            ! use SAI to initialize (v3.7)
             rtmassxy   (I,J) = 500.0                          ! these are all arbitrary and probably should be
             woodxy     (I,J) = 500.0                          ! in the table or read from initialization
             stblcpxy   (I,J) = 1000.0                         !
             fastcpxy   (I,J) = 1000.0                         !

	   END IF

          enddo
       enddo

       ! Given the soil layer thicknesses (in DZS), initialize the soil layer
       ! depths from the surface.
       ZSOIL(1)         = -DZS(1)          ! negative
       DO NS=2, NSOIL
          ZSOIL(NS)       = ZSOIL(NS-1) - DZS(NS)
       END DO

       ! Initialize snow/soil layer arrays ZSNSOXY, TSNOXY, SNICEXY, SNLIQXY,
       ! and ISNOWXY
       CALL snow_init ( ims , ime , jms , jme , its , itf , jts , jtf , 3 , &
            &           NSOIL , zsoil , snow , tgxy , snowh ,     &
            &           zsnsoxy , tsnoxy , snicexy , snliqxy , isnowxy )

       !initialize arrays for groundwater dynamics iopt_run=5

       if(iopt_run.eq.5) then
          IF ( PRESENT(smoiseq)     .AND. &
            PRESENT(smcwtdxy)    .AND. &
            PRESENT(rechxy)      .AND. &
            PRESENT(deeprechxy)  .AND. &
            PRESENT(areaxy)      .AND. &
            PRESENT(dx)          .AND. &
            PRESENT(dy)          .AND. &
            PRESENT(msftx)       .AND. &
            PRESENT(msfty)       .AND. &
            PRESENT(wtddt)       .AND. &
            PRESENT(stepwtd)     .AND. &
            PRESENT(dt)          .AND. &
            PRESENT(qrfsxy)      .AND. &
            PRESENT(qspringsxy)  .AND. &
            PRESENT(qslatxy)     .AND. &
            PRESENT(fdepthxy)    .AND. &
            PRESENT(ht)          .AND. &
            PRESENT(riverbedxy)  .AND. &
            PRESENT(eqzwt)       .AND. &
            PRESENT(rivercondxy) .AND. &
            PRESENT(pexpxy)            ) THEN

             STEPWTD = nint(WTDDT*60./DT)
             STEPWTD = max(STEPWTD,1)

              CALL groundwater_init ( &
      &       nsoil, zsoil , dzs  ,isltyp, ivgtyp,wtddt , &
      &       fdepthxy, ht, riverbedxy, eqzwt, rivercondxy, pexpxy , areaxy, zwtxy,   &
      &       smois,sh2o, smoiseq, smcwtdxy, deeprechxy, rechxy, qslatxy, qrfsxy, qspringsxy, &
      &       ids,ide, jds,jde, kds,kde,                    &
      &       ims,ime, jms,jme, kms,kme,                    &
      &       its,ite, jts,jte, kts,kte                     )

          ELSE
             CALL wrf_error_fatal ('Not enough fields to use groundwater option in Noah-MP')
          END IF
       endif

    ENDIF
  END SUBROUTINE NOAHMP_INIT

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

  SUBROUTINE SNOW_INIT ( ims , ime , jms , jme , its , itf , jts , jtf ,                  &
       &                 NSNOW , NSOIL , ZSOIL , SWE , TGXY , SNODEP ,                    &
       &                 ZSNSOXY , TSNOXY , SNICEXY ,SNLIQXY , ISNOWXY )
!------------------------------------------------------------------------------------------
!   Initialize snow arrays for Noah-MP LSM, based in input SNOWDEP, NSNOW
!   ISNOWXY is an index array, indicating the index of the top snow layer.  Valid indices
!           for snow layers range from 0 (no snow) and -1 (shallow snow) to (-NSNOW)+1 (deep snow).
!   TSNOXY holds the temperature of the snow layer.  Snow layers are initialized with
!          temperature = ground temperature [?].  Snow-free levels in the array have value 0.0
!   SNICEXY is the frozen content of a snow layer.  Initial estimate based on SNODEP and SWE
!   SNLIQXY is the liquid content of a snow layer.  Initialized to 0.0
!   ZNSNOXY is the layer depth from the surface.
!------------------------------------------------------------------------------------------
    IMPLICIT NONE
!------------------------------------------------------------------------------------------
    INTEGER, INTENT(IN)                              :: ims, ime, jms, jme
    INTEGER, INTENT(IN)                              :: its, itf, jts, jtf
    INTEGER, INTENT(IN)                              :: NSNOW
    INTEGER, INTENT(IN)                              :: NSOIL
    REAL,    INTENT(IN), DIMENSION(ims:ime, jms:jme) :: SWE
    REAL,    INTENT(IN), DIMENSION(ims:ime, jms:jme) :: SNODEP
    REAL,    INTENT(IN), DIMENSION(ims:ime, jms:jme) :: TGXY
    REAL,    INTENT(IN), DIMENSION(1:NSOIL)          :: ZSOIL

    INTEGER, INTENT(OUT), DIMENSION(ims:ime, jms:jme)                :: ISNOWXY ! Top snow layer index
    REAL,    INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:NSOIL,jms:jme) :: ZSNSOXY ! Snow/soil layer depth from surface [m]
    REAL,    INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:    0,jms:jme) :: TSNOXY  ! Snow layer temperature [K]
    REAL,    INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:    0,jms:jme) :: SNICEXY ! Snow layer ice content [mm]
    REAL,    INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:    0,jms:jme) :: SNLIQXY ! snow layer liquid content [mm]

! Local variables:
!   DZSNO   holds the thicknesses of the various snow layers.
!   DZSNOSO holds the thicknesses of the various soil/snow layers.
    INTEGER                           :: I,J,IZ
    REAL,   DIMENSION(-NSNOW+1:    0) :: DZSNO
    REAL,   DIMENSION(-NSNOW+1:NSOIL) :: DZSNSO

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

    DO J = jts , jtf
       DO I = its , itf
          IF ( SNODEP(I,J) < 0.025 ) THEN
             ISNOWXY(I,J) = 0
             DZSNO(-NSNOW+1:0) = 0.
          ELSE
             IF ( ( SNODEP(I,J) >= 0.025 ) .AND. ( SNODEP(I,J) <= 0.05 ) ) THEN
                ISNOWXY(I,J)    = -1
                DZSNO(0)  = SNODEP(I,J)
             ELSE IF ( ( SNODEP(I,J) > 0.05 ) .AND. ( SNODEP(I,J) <= 0.10 ) ) THEN
                ISNOWXY(I,J)    = -2
                DZSNO(-1) = SNODEP(I,J)/2.
                DZSNO( 0) = SNODEP(I,J)/2.
             ELSE IF ( (SNODEP(I,J) > 0.10 ) .AND. ( SNODEP(I,J) <= 0.25 ) ) THEN
                ISNOWXY(I,J)    = -2
                DZSNO(-1) = 0.05
                DZSNO( 0) = SNODEP(I,J) - DZSNO(-1)
             ELSE IF ( ( SNODEP(I,J) > 0.25 ) .AND. ( SNODEP(I,J) <= 0.45 ) ) THEN
                ISNOWXY(I,J)    = -3
                DZSNO(-2) = 0.05
                DZSNO(-1) = 0.5*(SNODEP(I,J)-DZSNO(-2))
                DZSNO( 0) = 0.5*(SNODEP(I,J)-DZSNO(-2))
             ELSE IF ( SNODEP(I,J) > 0.45 ) THEN
                ISNOWXY(I,J)     = -3
                DZSNO(-2) = 0.05
                DZSNO(-1) = 0.20
                DZSNO( 0) = SNODEP(I,J) - DZSNO(-1) - DZSNO(-2)
             ELSE
                CALL wrf_error_fatal("Problem with the logic assigning snow layers.")
             END IF
          END IF

          TSNOXY (I,-NSNOW+1:0,J) = 0.
          SNICEXY(I,-NSNOW+1:0,J) = 0.
          SNLIQXY(I,-NSNOW+1:0,J) = 0.
          DO IZ = ISNOWXY(I,J)+1 , 0
             TSNOXY(I,IZ,J)  = TGXY(I,J)  ! [k]
             SNLIQXY(I,IZ,J) = 0.00
             SNICEXY(I,IZ,J) = 1.00 * DZSNO(IZ) * (SWE(I,J)/SNODEP(I,J))  ! [kg/m3]
          END DO

          ! Assign local variable DZSNSO, the soil/snow layer thicknesses, for snow layers
          DO IZ = ISNOWXY(I,J)+1 , 0
             DZSNSO(IZ) = -DZSNO(IZ)
          END DO

          ! Assign local variable DZSNSO, the soil/snow layer thicknesses, for soil layers
          DZSNSO(1) = ZSOIL(1)
          DO IZ = 2 , NSOIL
             DZSNSO(IZ) = (ZSOIL(IZ) - ZSOIL(IZ-1))
          END DO

          ! Assign ZSNSOXY, the layer depths, for soil and snow layers
          ZSNSOXY(I,ISNOWXY(I,J)+1,J) = DZSNSO(ISNOWXY(I,J)+1)
          DO IZ = ISNOWXY(I,J)+2 , NSOIL
             ZSNSOXY(I,IZ,J) = ZSNSOXY(I,IZ-1,J) + DZSNSO(IZ)
          ENDDO

       END DO
    END DO

  END SUBROUTINE SNOW_INIT
! ==================================================================================================
! ----------------------------------------------------------------------
    SUBROUTINE GROUNDWATER_INIT (   &
            &            NSOIL , ZSOIL , DZS, ISLTYP, IVGTYP, WTDDT , &
            &            FDEPTH, TOPO, RIVERBED, EQWTD, RIVERCOND, PEXP , AREA ,WTD ,  &
            &            SMOIS,SH2O, SMOISEQ, SMCWTDXY, DEEPRECHXY, RECHXY ,  &
            &            QSLATXY, QRFSXY, QSPRINGSXY,                  &
            &            ids,ide, jds,jde, kds,kde,                    &
            &            ims,ime, jms,jme, kms,kme,                    &
            &            its,ite, jts,jte, kts,kte                     )


  USE NOAHMP_TABLES, ONLY : BEXP_TABLE,SMCMAX_TABLE,PSISAT_TABLE,SMCWLT_TABLE,DWSAT_TABLE,DKSAT_TABLE, &
                                ISURBAN_TABLE, ISICE_TABLE ,ISWATER_TABLE
  USE module_sf_noahmp_groundwater, ONLY : LATERALFLOW

! ----------------------------------------------------------------------
  IMPLICIT NONE
! ----------------------------------------------------------------------

    INTEGER,  INTENT(IN   )   ::     ids,ide, jds,jde, kds,kde,  &
         &                           ims,ime, jms,jme, kms,kme,  &
         &                           its,ite, jts,jte, kts,kte
    INTEGER, INTENT(IN)                              :: NSOIL
    REAL,   INTENT(IN)                               ::     WTDDT
    REAL,    INTENT(IN), DIMENSION(1:NSOIL)          :: ZSOIL,DZS
    INTEGER, INTENT(IN), DIMENSION(ims:ime, jms:jme) :: ISLTYP, IVGTYP
    REAL,    INTENT(IN), DIMENSION(ims:ime, jms:jme) :: FDEPTH, TOPO, RIVERBED, EQWTD, RIVERCOND, PEXP , AREA
    REAL,    INTENT(INOUT), DIMENSION(ims:ime, jms:jme) :: WTD
    REAL,     DIMENSION( ims:ime , 1:nsoil, jms:jme ), &
         &    INTENT(INOUT)   ::                          SMOIS, &
         &                                                 SH2O, &
         &                                                 SMOISEQ
    REAL,    INTENT(INOUT), DIMENSION(ims:ime, jms:jme) ::  &
                                                           SMCWTDXY, &
                                                           DEEPRECHXY, &
                                                           RECHXY, &
                                                           QSLATXY, &
                                                           QRFSXY, &
                                                           QSPRINGSXY
! local
    INTEGER  :: I,J,K,ITER,itf,jtf
    REAL :: BEXP,SMCMAX,PSISAT,SMCWLT,DWSAT,DKSAT
    REAL :: FRLIQ,SMCEQDEEP
    REAL :: DELTAT,RCOND
    REAL :: AA,BBB,CC,DD,DX,FUNC,DFUNC,DDZ,EXPON,SMC,FLUX
    REAL, DIMENSION(1:NSOIL) :: SMCEQ
    REAL,      DIMENSION( ims:ime, jms:jme )    :: QLAT, QRF
    INTEGER,   DIMENSION( ims:ime, jms:jme )    :: LANDMASK !-1 for water (ice or no ice) and glacial areas, 1 for land where the LSM does its soil moisture calculations


       itf=min0(ite,ide-1)
       jtf=min0(jte,jde-1)

!first compute lateral flow and flow to rivers to initialize deep soil moisture

    DELTAT = WTDDT * 60. !timestep in seconds for this calculation

    WHERE(IVGTYP.NE.ISWATER_TABLE.AND.IVGTYP.NE.ISICE_TABLE)
         LANDMASK=1
    ELSEWHERE
         LANDMASK=-1
    ENDWHERE

!Calculate lateral flow

    QLAT = 0.
CALL LATERALFLOW(ISLTYP,WTD,QLAT,FDEPTH,TOPO,LANDMASK,DELTAT,AREA       &
                        ,ids,ide,jds,jde,kds,kde                      &
                        ,ims,ime,jms,jme,kms,kme                      &
                        ,its,ite,jts,jte,kts,kte                      )


!compute flux from grounwater to rivers in the cell

    DO J=jts,jtf
       DO I=its,itf
          IF(LANDMASK(I,J).GT.0)THEN
             IF(WTD(I,J) .GT. RIVERBED(I,J) .AND.  EQWTD(I,J) .GT. RIVERBED(I,J)) THEN
               RCOND = RIVERCOND(I,J) * EXP(PEXP(I,J)*(WTD(I,J)-EQWTD(I,J)))
             ELSE
               RCOND = RIVERCOND(I,J)
             ENDIF
             QRF(I,J) = RCOND * (WTD(I,J)-RIVERBED(I,J)) * DELTAT/AREA(I,J)
!for now, dont allow it to go from river to groundwater
             QRF(I,J) = MAX(QRF(I,J),0.)
          ELSE
             QRF(I,J) = 0.
          ENDIF
       ENDDO
    ENDDO

!now compute eq. soil moisture, change soil moisture to be compatible with the water table and compute deep soil moisture

       DO J = jts,jtf
          DO I = its,itf
             BEXP   =   BEXP_TABLE(ISLTYP(I,J))
             SMCMAX = SMCMAX_TABLE(ISLTYP(I,J))
             SMCWLT = SMCWLT_TABLE(ISLTYP(I,J))
             IF(IVGTYP(I,J)==ISURBAN_TABLE)THEN
                 SMCMAX = 0.45
                 SMCWLT = 0.40
             ENDIF
             DWSAT  =   DWSAT_TABLE(ISLTYP(I,J))
             DKSAT  =   DKSAT_TABLE(ISLTYP(I,J))
             PSISAT = -PSISAT_TABLE(ISLTYP(I,J))
           IF ( ( BEXP > 0.0 ) .AND. ( smcmax > 0.0 ) .AND. ( -psisat > 0.0 ) ) THEN
             !initialize equilibrium soil moisture for water table diagnostic
                    CALL EQSMOISTURE(NSOIL ,  ZSOIL , SMCMAX , SMCWLT ,DWSAT, DKSAT  ,BEXP  , & !in
                                     SMCEQ                          )  !out

             SMOISEQ (I,1:NSOIL,J) = SMCEQ (1:NSOIL)


              !make sure that below the water table the layers are saturated and initialize the deep soil moisture
             IF(WTD(I,J) < ZSOIL(NSOIL)-DZS(NSOIL)) THEN

!initialize deep soil moisture so that the flux compensates qlat+qrf
!use Newton-Raphson method to find soil moisture

                         EXPON = 2. * BEXP + 3.
                         DDZ = ZSOIL(NSOIL) - WTD(I,J)
                         CC = PSISAT/DDZ
                         FLUX = (QLAT(I,J)-QRF(I,J))/DELTAT

                         SMC = 0.5 * SMCMAX

                         DO ITER = 1, 100
                           DD = (SMC+SMCMAX)/(2.*SMCMAX)
                           AA = -DKSAT * DD  ** EXPON
                           BBB = CC * ( (SMCMAX/SMC)**BEXP - 1. ) + 1.
                           FUNC =  AA * BBB - FLUX
                           DFUNC = -DKSAT * (EXPON/(2.*SMCMAX)) * DD ** (EXPON - 1.) * BBB &
                                   + AA * CC * (-BEXP) * SMCMAX ** BEXP * SMC ** (-BEXP-1.)

                           DX = FUNC/DFUNC
                           SMC = SMC - DX
                           IF ( ABS (DX) < 1.E-6)EXIT
                         ENDDO

                  SMCWTDXY(I,J) = MAX(SMC,1.E-4)

             ELSEIF(WTD(I,J) < ZSOIL(NSOIL))THEN
                  SMCEQDEEP = SMCMAX * ( PSISAT / ( PSISAT - DZS(NSOIL) ) ) ** (1./BEXP)
!                  SMCEQDEEP = MAX(SMCEQDEEP,SMCWLT)
                  SMCEQDEEP = MAX(SMCEQDEEP,1.E-4)
                  SMCWTDXY(I,J) = SMCMAX * ( WTD(I,J) -  (ZSOIL(NSOIL)-DZS(NSOIL))) + &
                                  SMCEQDEEP * (ZSOIL(NSOIL) - WTD(I,J))

             ELSE !water table within the resolved layers
                  SMCWTDXY(I,J) = SMCMAX
                  DO K=NSOIL,2,-1
                     IF(WTD(I,J) .GE. ZSOIL(K-1))THEN
                          FRLIQ = SH2O(I,K,J) / SMOIS(I,K,J)
                          SMOIS(I,K,J) = SMCMAX
                          SH2O(I,K,J) = SMCMAX * FRLIQ
                     ELSE
                          IF(SMOIS(I,K,J).LT.SMCEQ(K))THEN
                              WTD(I,J) = ZSOIL(K)
                          ELSE
                              WTD(I,J) = ( SMOIS(I,K,J)*DZS(K) - SMCEQ(K)*ZSOIL(K-1) + SMCMAX*ZSOIL(K) ) / &
                                         (SMCMAX - SMCEQ(K))
                          ENDIF
                          EXIT
                     ENDIF
                  ENDDO
             ENDIF
            ELSE
              SMOISEQ (I,1:NSOIL,J) = SMCMAX
              SMCWTDXY(I,J) = SMCMAX
              WTD(I,J) = 0.
            ENDIF

!zero out some arrays

             DEEPRECHXY(I,J) = 0.
             RECHXY(I,J) = 0.
             QSLATXY(I,J) = 0.
             QRFSXY(I,J) = 0.
             QSPRINGSXY(I,J) = 0.

          ENDDO
       ENDDO




    END  SUBROUTINE GROUNDWATER_INIT
! ==================================================================================================
! ----------------------------------------------------------------------
  SUBROUTINE EQSMOISTURE(NSOIL  ,  ZSOIL , SMCMAX , SMCWLT, DWSAT , DKSAT ,BEXP , & !in
                         SMCEQ                          )  !out
! ----------------------------------------------------------------------
  IMPLICIT NONE
! ----------------------------------------------------------------------
! input
  INTEGER,                         INTENT(IN) :: NSOIL !no. of soil layers
  REAL, DIMENSION(       1:NSOIL), INTENT(IN) :: ZSOIL !depth of soil layer-bottom [m]
  REAL,                            INTENT(IN) :: SMCMAX , SMCWLT, BEXP , DWSAT, DKSAT
!output
  REAL,  DIMENSION(      1:NSOIL), INTENT(OUT) :: SMCEQ  !equilibrium soil water  content [m3/m3]
!local
  INTEGER                                     :: K , ITER
  REAL                                        :: DDZ , SMC, FUNC, DFUNC , AA, BB , EXPON, DX

!gmmcompute equilibrium soil moisture content for the layer when wtd=zsoil(k)


   DO K=1,NSOIL

            IF ( K == 1 )THEN
                DDZ = -ZSOIL(K+1) * 0.5
            ELSEIF ( K < NSOIL ) THEN
                DDZ = ( ZSOIL(K-1) - ZSOIL(K+1) ) * 0.5
            ELSE
                DDZ = ZSOIL(K-1) - ZSOIL(K)
            ENDIF

!use Newton-Raphson method to find eq soil moisture

            EXPON = BEXP +1.
            AA = DWSAT/DDZ
            BB = DKSAT / SMCMAX ** EXPON

            SMC = 0.5 * SMCMAX

         DO ITER = 1, 100
            FUNC = (SMC - SMCMAX) * AA +  BB * SMC ** EXPON
            DFUNC = AA + BB * EXPON * SMC ** BEXP

            DX = FUNC/DFUNC
            SMC = SMC - DX
            IF ( ABS (DX) < 1.E-6)EXIT
         ENDDO

!             SMCEQ(K) = MIN(MAX(SMC,SMCWLT),SMCMAX*0.99)
             SMCEQ(K) = MIN(MAX(SMC,1.E-4),SMCMAX*0.99)
   ENDDO

END  SUBROUTINE EQSMOISTURE
!
!------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------
!
END MODULE module_sf_noahmpdrv