MODULE module_sf_noahmpdrv
  USE module_sf_noahmplsm
  USE module_sf_urban
  USE module_model_constants, ONLY : R_D, CP, XLF, XLV, RHOWATER, KARMAN
  USE module_sf_noahdrv, ONLY : SOIL_VEG_GEN_PARM
  USE module_sf_noah_seaice
  USE module_sf_noahlsm_glacial_only
  USE MODULE_RA_GFDLETA, ONLY: CAL_MON_DAY
CONTAINS
  SUBROUTINE noahmplsm(DZ8W,QV3D,P8W3D,T3D,TSK, &
       HFX,QFX,LH,GRDFLX, QGH,GSW,SWDOWN,GLW,SMSTAV,SMSTOT, &
       SFCRUNOFF, UDRUNOFF,IVGTYP,ISLTYP,VEGFRA, &
       ALBEDO,ALBBCK,ZNT,Z0,TMN,XLAND,XICE,XICE_THRESHOLD,ISICE,EMISS,EMBCK, &
       SNOWC,QSFC,RAINBL, &
       num_soil_layers,DT,DZS,ITIMESTEP, &
       SMOIS,TSLB,SNOW,CANWAT, &
       CHS,CHS2,CQS2,CPM,ROVCP,SR,chklowq,qz0, &
       myj,RIB,frpcpn, &
       SH2O,SNOWH, &
       U_PHY,V_PHY, &
       COSZ_URB2D, XLAT_URB2D, &
       SNOALB, &
       SNOTIME,ACSNOM,ACSNOW, &
       idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc ,iopt_frz ,iopt_inf , &
       iopt_rad ,iopt_alb ,iopt_snf ,iopt_tbot,iopt_stc , &
       isnowxy ,tvxy ,tgxy ,canicexy , &
       canliqxy ,eahxy ,tahxy ,cmxy ,chxy , &
       fwetxy ,sneqvoxy ,alboldxy ,qsnowxy ,wslakexy ,zwtxy ,waxy , &
       wtxy ,tsnoxy ,zsnsoxy ,snicexy ,snliqxy ,lfmassxy ,rtmassxy , &
       stmassxy ,woodxy ,stblcpxy ,fastcpxy ,xlaixy ,xsaixy , &
       tradxy ,tsxy ,neexy ,gppxy ,nppxy ,fvegxy ,qinxy , &
       runsfxy ,runsbxy ,ecanxy ,edirxy ,etranxy ,fsaxy ,firaxy , &
       aparxy ,psnxy ,savxy ,sagxy , &
       fsnoxy ,YR ,JULIAN , &
       potevp, &
       qcxy ,pblhxy ,isurban ,iz0tlnd ,dx , &
       chstarxy ,t2mvxy ,t2mbxy ,rssunxy ,rsshaxy, bgapxy ,wgapxy ,gapxy, &
       tgvxy ,tgbxy ,q2mvxy ,q2mbxy ,shdmaxxy, chvxy ,chbxy , &
       ids,ide, jds,jde, kds,kde, &
       ims,ime, jms,jme, kms,kme, &
       its,ite, jts,jte, kts,kte )
    IMPLICIT NONE
    INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
         & ims,ime, jms,jme, kms,kme, &
         & its,ite, jts,jte, kts,kte
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U_PHY
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: V_PHY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT_URB2D
    REAL, DIMENSION( ims:ime, jms:jme ) , &
         & INTENT(IN ) :: TMN, &
         & XLAND, &
         & XICE, &
         & VEGFRA, &
         & SNOALB, &
         & GSW, &
         & SWDOWN, &
         & GLW, &
         & Z0, &
         & ALBBCK, &
         & RAINBL, &
         & EMBCK, &
         & SR
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
         INTENT(IN ) :: QV3D, &
         p8w3D, &
         DZ8W, &
         T3D
    REAL, DIMENSION( ims:ime, jms:jme ) , &
         INTENT(INOUT ) :: QGH, &
         CHS, &
         CPM
    INTEGER, DIMENSION( ims:ime, jms:jme ) , &
         INTENT(IN ) :: IVGTYP, &
         ISLTYP
    INTEGER, INTENT(IN) :: num_soil_layers,ITIMESTEP
    REAL, INTENT(IN ) :: DT,ROVCP,XICE_THRESHOLD
    INTEGER, INTENT(IN ) :: ISICE
    REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
    REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
         & INTENT(INOUT) :: SMOIS, &
         & SH2O, &
         & TSLB
    REAL, DIMENSION( ims:ime, jms:jme ) , &
         & INTENT(INOUT) :: TSK, &
         & HFX, &
         & QFX, &
         & LH, &
         & GRDFLX, &
         & QSFC, &
         & CQS2, &
         & CHS2, &
         & SNOW, &
         & SNOWC, &
         & SNOWH, &
         & CANWAT, &
         & SMSTAV, &
         & SMSTOT, &
         & SFCRUNOFF, &
         & UDRUNOFF, &
         & ACSNOM, &
         & ACSNOW, &
         & EMISS, &
         & POTEVP, &
         & RIB, &
         & ALBEDO, &
         & ZNT
    REAL, DIMENSION( ims:ime, jms:jme ) , &
         INTENT(OUT) :: CHKLOWQ
    REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: QZ0
    INTEGER, INTENT(IN) :: idveg
    INTEGER, INTENT(IN) :: iopt_crs
    INTEGER, INTENT(IN) :: iopt_btr
    INTEGER, INTENT(IN) :: iopt_run
    INTEGER, INTENT(IN) :: iopt_sfc
    INTEGER, INTENT(IN) :: iopt_frz
    INTEGER, INTENT(IN) :: iopt_inf
    INTEGER, INTENT(IN) :: iopt_rad
    INTEGER, INTENT(IN) :: iopt_alb
    INTEGER, INTENT(IN) :: iopt_snf
    INTEGER, INTENT(IN) :: iopt_tbot
    INTEGER, INTENT(IN) :: iopt_stc
    INTEGER, INTENT(IN) :: YR
    REAL, INTENT(IN) :: JULIAN
    INTEGER, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: isnowxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tvxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tgxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canicexy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canliqxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: eahxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tahxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: cmxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: chxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fwetxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: sneqvoxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: alboldxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: qsnowxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wslakexy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: zwtxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: waxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wtxy
    REAL, DIMENSION(ims:ime,-2:num_soil_layers,jms:jme), INTENT(INOUT) :: zsnsoxy
    REAL, DIMENSION(ims:ime,-2: 0,jms:jme), INTENT(INOUT) :: tsnoxy
    REAL, DIMENSION(ims:ime,-2: 0,jms:jme), INTENT(INOUT) :: snicexy
    REAL, DIMENSION(ims:ime,-2: 0,jms:jme), INTENT(INOUT) :: snliqxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: lfmassxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: rtmassxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stmassxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: woodxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stblcpxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fastcpxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: xlaixy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: xsaixy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: SNOTIME
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: tradxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: tsxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: neexy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: gppxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: nppxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: fvegxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: qinxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: runsfxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: runsbxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: ecanxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: edirxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: etranxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: fsaxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: firaxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: aparxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: psnxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: savxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: sagxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: fsnoxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: chstarxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: t2mvxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: t2mbxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: q2mvxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: q2mbxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: qcxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: pblhxy
    INTEGER , INTENT(IN) :: isurban
    INTEGER , INTENT(IN) :: iz0tlnd
    REAL , INTENT(IN) :: dx
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: rssunxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: rsshaxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: bgapxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: wgapxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: gapxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: tgvxy
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: tgbxy
    REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: shdmaxxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: chvxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: chbxy
    INTEGER :: YEARLEN
    REAL :: ETP, SSOIL,EC, ESNOW, &
         FLX1,FLX2,FLX3,DEW,FDOWN,RC,PC,FFROZP
    INTEGER :: isnow
    REAL, DIMENSION(-2:num_soil_layers) :: stc
    REAL, DIMENSION( 1:num_soil_layers) :: smc
    REAL, DIMENSION( 1:num_soil_layers) :: smh2o
    REAL :: tv
    REAL :: tg
    REAL :: canice
    REAL :: canliq
    REAL :: snowd
    REAL :: swe
    REAL :: eah
    REAL :: tah
    REAL :: cm
    REAL :: ch
    REAL :: fwet
    REAL :: sneqvo
    REAL :: albold
    REAL :: qsnow
    REAL :: wslake
    REAL :: zwt
    REAL :: wa
    REAL :: wt
    REAL, DIMENSION(-2:num_soil_layers) :: zsnso
    REAL, DIMENSION(-2: 0) :: tsno
    REAL, DIMENSION(-2: 0) :: snice
    REAL, DIMENSION(-2: 0) :: snliq
    REAL :: lfmass
    REAL :: rtmass
    REAL :: stmass
    REAL :: wood
    REAL :: stblcp
    REAL :: fastcp
    REAL :: plai
    REAL :: psai
    REAL :: chstar2
    REAL :: cqstar2
    REAL :: chstar
    REAL :: tstar
    REAL :: t2mv
    REAL :: t2mb
    REAL :: q2mv
    REAL :: q2mb
    REAL :: qc
    REAL :: t2m
    REAL :: pblh
    REAL :: qsfc1d
    REAL, DIMENSION(ims:ime,jms:jme) :: tstarxy
    REAL, DIMENSION(ims:ime,jms:jme) :: chstar2xy
    REAL :: rssun
    REAL :: rssha
    REAL :: bgap
    REAL :: wgap
    REAL :: gap
    REAL :: tgv
    REAL :: tgb
    REAL :: snowhk
    REAL :: snotime1
    REAL :: qv1d
    REAL :: dz8w1d
    REAL :: shdmax
    REAL :: chv
    REAL :: chb
    REAL :: trad
    REAL :: ts
    REAL :: nee
    REAL :: gpp
    REAL :: npp
    REAL :: fveg
    REAL :: qin
    REAL :: runsf
    REAL :: runsb
    REAL :: ecan
    REAL :: esoil
    REAL :: etran
    REAL :: fsa
    REAL :: fira
    REAL :: fsh
    REAL :: flh
    REAL :: apar
    REAL :: psn
    REAL :: sav
    REAL :: sag
    REAL :: fsno
    REAL :: salb
    REAL :: errwat
    REAL :: qmelt
    REAL :: ponding
    REAL :: ponding1
    REAL :: ponding2
    real :: fsr
    real :: fcev
    real :: fgev
    real :: fctr
    real, dimension(-2: 0) :: ficeold
    INTEGER :: ILOC
    INTEGER :: JLOC
    INTEGER :: ISC
    INTEGER :: IST
    LOGICAL, INTENT(IN ) :: myj,frpcpn
    LOGICAL, PARAMETER :: LOCAL=.false.
    LOGICAL :: FRZGRA, SNOWNG
    LOGICAL :: IPRINT
    INTEGER :: I,J, ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
    INTEGER :: NROOT
    INTEGER :: KZ ,K
    INTEGER :: NS
    REAL :: DQSDT2, LWDN, PRCP, PSFC, UU, VV, CO2AIR, O2AIR, &
         & Q2SAT,Q2SATI,SFCPRS,SFCTMP,SHDFAC,SNOALB1, &
         & SOLDN,TBOT,ZLVL, Q2K,ALBBRD, ETA, ETA_KINEMATIC, &
         & EMBRD, FOLN, LAT, &
         & Z0K,RUNOFF1,RUNOFF2,SOLNET,E2SAT,SFCTSNO
    REAL :: RIBB
    REAL :: FDTW
    REAL :: EMISSI
    REAL :: SNCOVR,SNEQV,CHK,TH2
    REAL :: SMCMAX,SNOMLT,SOILM,SOILW,Q1,T1
    REAL :: Z0BRD
    REAL :: COSZ
    REAL, DIMENSION(1:num_soil_layers):: SLDPTH,SWC
    REAL, DIMENSION(1:num_soil_layers):: STCNEW
    REAL, DIMENSION(1:num_soil_layers) :: ZSOIL, RTDIS
    REAL, PARAMETER :: TRESH=.95E0, A2=17.67,A3=273.15,A4=29.65, &
         T0=273.16E0, ELWV=2.50E6, A23M4=A2*(A3-A4)
    REAL, PARAMETER :: CAPA=R_D/CP
    REAL :: APELM
    REAL :: APES
    REAL :: SFCTH2
    FDTW=DT/(XLV*RHOWATER)
    IPRINT=.false.
    SLOPETYP=1
    YEARLEN = 365
    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
    NSOIL=num_soil_layers
    DO NS=1,NSOIL
       SLDPTH(NS)=DZS(NS)
    ENDDO
    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 )
    ISC = 4
    ZSOIL(1) = -SLDPTH(1)
    DO KZ = 2, NSOIL
       ZSOIL(KZ) = -SLDPTH(KZ) + ZSOIL(KZ-1)
    END DO
    FOLN = 1.0
    DO J=jts,jte
       IF(ITIMESTEP.EQ.1)THEN
          DO I=its,ite
             IF((XLAND(I,J)-1.5).GE.0.)THEN
                IF(XICE(I,J).EQ.1..and.IPRINT)PRINT*,' sea-ice at water point, I=',I, &
                     'J=',J
                SMSTAV(I,J)=1.0
                SMSTOT(I,J)=1.0
                DO NS=1,NSOIL
                   SMOIS(I,NS,J)=1.0
                   TSLB(I,NS,J)=273.16
                ENDDO
             ELSE
                IF(XICE(I,J).EQ.1.)THEN
                   SMSTAV(I,J)=1.0
                   SMSTOT(I,J)=1.0
                   DO NS=1,NSOIL
                      SMOIS(I,NS,J)=1.0
                   ENDDO
                ENDIF
             ENDIF
          ENDDO
       ENDIF
       DO I=its,ite
          PSFC=P8w3D(i,1,j)
          SFCPRS=(P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
          Q2K=QV3D(i,1,j)/(1.0+QV3D(i,1,j))
          Q2SAT=QGH(I,J)/(1.0+QGH(I,J))
          IF((myj).AND.(Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
             CHKLOWQ(I,J)=0.
          ELSE
             CHKLOWQ(I,J)=1.
          ENDIF
          SFCTMP=T3D(i,1,j)
          ZLVL=0.5*DZ8W(i,1,j)
          APES=(1.E5/PSFC)**CAPA
          APELM=(1.E5/SFCPRS)**CAPA
          SFCTH2=SFCTMP*APELM
          TH2=SFCTH2/APES
          EMISSI = EMISS(I,J)
          LWDN=GLW(I,J)
          SOLDN=SWDOWN(I,J)
          SOLNET=SOLDN*(1.-ALBEDO(I,J))
          PRCP=RAINBL(i,j)/DT
          VEGTYP=IVGTYP(I,J)
          SOILTYP=ISLTYP(I,J)
          SHDFAC=VEGFRA(I,J)/100.
          T1=TSK(I,J)
          CHK=CHS(I,J)
          SNOALB1=SNOALB(I,J)
          IF(FRPCPN) THEN
             FFROZP=SR(I,J)
          ELSE
             IF (SFCTMP <= 273.15) THEN
                FFROZP = 1.0
             ELSE
                FFROZP = 0.0
             ENDIF
          ENDIF
          IF((XLAND(I,J)-1.5).GE.0.)THEN
          ELSE
             IF (XICE(I,J) .GT. 0.5) THEN
                ICE=1
             ELSE
                ICE=0
             ENDIF
             DQSDT2=Q2SAT*A23M4/(SFCTMP-A4)**2
             IF(SNOW(I,J).GT.0.0)THEN
                SFCTSNO=SFCTMP
                E2SAT=611.2*EXP(6174.*(1./273.15 - 1./SFCTSNO))
                Q2SATI=0.622*E2SAT/(SFCPRS-E2SAT)
                Q2SATI=Q2SATI/(1.0+Q2SATI)
                IF(T1 .GT. 273.15)THEN
                   Q2SAT=Q2SAT*(1.-SNOWC(I,J)) + Q2SATI*SNOWC(I,J)
                   DQSDT2=DQSDT2*(1.-SNOWC(I,J)) + Q2SATI*6174./(SFCTSNO**2)*SNOWC(I,J)
                ELSE
                   Q2SAT=Q2SATI
                   DQSDT2=Q2SATI*6174./(SFCTSNO**2)
                ENDIF
                IF(T1 .GT. 273. .AND. SNOWC(I,J) .GT. 0.)DQSDT2=DQSDT2*(1.-SNOWC(I,J))
             ENDIF
             IF(ICE.EQ.0)THEN
                TBOT=TMN(I,J)
             ELSE
                TBOT=271.16
             ENDIF
             IF(VEGTYP.EQ.25) SHDFAC=0.0000
             IF(VEGTYP.EQ.26) SHDFAC=0.0000
             IF(VEGTYP.EQ.27) SHDFAC=0.0000
             IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.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
             ALBBRD=ALBBCK(I,J)
             Z0BRD=Z0(I,J)
             EMBRD=EMBCK(I,J)
             RIBB=RIB(I,J)
             SNOTIME1 = SNOTIME(I,J)
             IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == 31 .or. &
                  IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
                VEGTYP = ISURBAN
             ENDIF
             IST = 1
             IF(VEGTYP == 16) IST = 2
             CALL REDPRM (VEGTYP,SOILTYP,SLOPETYP,SLDPTH,ZSOIL,NSOIL,ISURBAN)
             UU = U_PHY(I,1,J)
             VV = V_PHY(I,1,J)
             CO2AIR = 395.E-06 * SFCPRS
             O2AIR = 0.209 * SFCPRS
             COSZ = COSZ_URB2D(I,J)
             LAT = XLAT_URB2D(I,J)
             isnow = isnowxy (i,j)
             stc (isnow+1: 0) = tsnoxy (i,isnow+1: 0,j)
             stc ( 1:nsoil) = tslb (i, 1:nsoil,j)
             smc ( 1:nsoil) = smois (i, 1:nsoil,j)
             smh2o( 1:nsoil) = sh2o (i, 1:nsoil,j)
             tv = tvxy (i,j)
             tg = tgxy (i,j)
             canliq = canliqxy(i,j)
             canice = canicexy(i,j)
             snowd = snowh (i,j)
             swe = snow (i,j)
             eah = eahxy (i,j)
             tah = tahxy (i,j)
             cm = cmxy (i,j)
             ch = chxy (i,j)
             chstar = chs (i,j)
             chstar2 = chs2 (i,j)
             cqstar2 = cqs2 (i,j)
             tstar = T1
             qc = qcxy (i,j)
             pblh = pblhxy (i,j)
             qsfc1d = qsfc (i,j)
             t2mv = t2mvxy (i,j)
             t2mb = t2mbxy (i,j)
             q2mv = q2mvxy (i,j)
             q2mb = q2mbxy (i,j)
             qv1d = qv3d (i,1,j)
             dz8w1d = dz8w (i,1,j)
             shdmax = shdmaxxy (i,j)/100.
             fwet = fwetxy (i,j)
             sneqvo = sneqvoxy(i,j)
             albold = alboldxy(i,j)
             qsnow = qsnowxy (i,j)
             wslake = wslakexy(i,j)
             zwt = zwtxy (i,j)
             wa = waxy (i,j)
             wt = wtxy (i,j)
             zsnso(isnow+1:nsoil) = zsnsoxy (i,isnow+1:nsoil,j)
             snice(isnow+1: 0) = snicexy (i,isnow+1: 0,j)
             snliq(isnow+1: 0) = snliqxy (i,isnow+1: 0,j)
             lfmass = lfmassxy(i,j)
             rtmass = rtmassxy(i,j)
             stmass = stmassxy(i,j)
             wood = woodxy (i,j)
             stblcp = stblcpxy(i,j)
             fastcp = fastcpxy(i,j)
             plai = xlaixy (i,j)
             psai = xsaixy (i,j)
             ficeold(isnow+1:0) = snicexy(i,isnow+1:0,j) &
                  /(snicexy(i,isnow+1:0,j)+snliqxy(i,isnow+1:0,j))
    IF ( XICE(I,J) >= XICE_THRESHOLD ) THEN
       SH2O (i,1:nsoil,j) = 1.0
       XLAIXY(i,j) = 0.01
       cycle
    ELSE IF ( VEGTYP == ISICE ) THEN
       SNCOVR = SNOWC(I,J)
       swe = swe*0.001
       if ( (swe.ne.0..AND.snowd.eq.0.).or.(snowd.le.swe) )THEN
             snowd= 5.*swe
       endif
       CALL SFLX_GLACIAL(I,J,ISICE,FFROZP,DT,ZLVL,NSOIL,SLDPTH, &
            & LWDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K, &
            & TH2,Q2SAT,DQSDT2, &
            & ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, &
            & tstar,STC(1:NSOIL),snowd,swe,salb,chstar, &
            & ETA,fsh, ETA_KINEMATIC,FDOWN, &
            & ESNOW,DEW, &
            & ETP,SSOIL, &
            & FLX1,FLX2,FLX3, &
            & SNOMLT,SNCOVR, &
            & runsf, &
            & Q1, &
            & SNOTIME1, &
            & RIBB)
       tgb = sfctmp
       tgv = 0.0
       swe = swe*1000.
       plai = 0.01
       smc = 1.00
       smh2o = 1.00
       runsb = 0.00
       fgev = ETA
       fcev = 0.
       fctr = 0.
       soilm = 1.0
       SNOWC(I,J) = 1.0
       QFX(I,J) = eta_kinematic
       POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
       CHS2(I,J) = CQS2(I,J)
       IF ( Q1 .GT. QSFC(I,J) ) THEN
          CQS2(I,J) = CHS(I,J)
       ENDIF
    ELSE
             nee = -1.E36
             npp = -1.E36
             CALL NOAHMP_SFLX (&
                  I , J , LAT , YEARLEN , JULIAN , COSZ , &
                  DT , DX , DZ8W1D , NSOIL , ZSOIL , 3 , &
                  SHDFAC , SHDMAX , VEGTYP , ISURBAN , ICE , IST , &
                  ISC , &
                  IZ0TLND , &
                  SFCTMP , SFCPRS , PSFC , UU , VV , QV1D , &
                  QC , SOLDN , LWDN , PRCP , TBOT , CO2AIR , &
                  O2AIR , FOLN , FICEOLD , PBLH , &
                  ZLVL , ALBOLD , SNEQVO , &
                  STC , SMH2O , SMC , TAH , EAH , FWET , &
                  CANLIQ , CANICE , TV , TG , QSFC1D , QSNOW , &
                  ISNOW , ZSNSO , SNOWD , SWE , SNICE , SNLIQ , &
                  ZWT , WA , WT , WSLAKE , LFMASS , RTMASS , &
                  STMASS , WOOD , STBLCP , FASTCP , PLAI , PSAI , &
                  CM , CH , CHSTAR , &
                  FSA , FSR , FIRA , FSH , SSOIL , FCEV , &
                  FGEV , FCTR , ECAN , ETRAN , ESOIL , TRAD , &
                  TS , TGB , TGV , T2MV , T2MB , TSTAR , &
                  Q1 , Q2MV , Q2MB , RUNSF , RUNSB , APAR , &
                  PSN , SAV , SAG , FSNO , NEE , GPP , &
                  NPP , FVEG , SALB , QMELT , PONDING , PONDING1, &
                  PONDING2, RSSUN , RSSHA , BGAP , WGAP , GAP , &
                  ERRWAT , CHV , CHB , EMISSI)
             chs2 (i,j) = chstar2
             cqs2 (i,j) = cqstar2
             QFX (I,J) = ecan + esoil + etran
             SNOWC (I,J) = fsno
   ENDIF
             isnowxy (i,j) = isnow
             canliqxy (i,j) = canliq
             canicexy (i,j) = canice
             snowh (i,j) = snowd
             snow (i,j) = swe
             zsnsoxy (i,isnow+1:nsoil,j) = zsnso (isnow+1:nsoil)
             tslb (i, 1:nsoil,j) = stc ( 1:nsoil)
             tsnoxy (i,isnow+1: 0,j) = stc (isnow+1: 0)
             smois (i, 1:nsoil,j) = smc ( 1:nsoil)
             sh2o (i, 1:nsoil,j) = smh2o ( 1:nsoil)
             snicexy (i,isnow+1: 0,j) = snice (isnow+1: 0)
             snliqxy (i,isnow+1: 0,j) = snliq (isnow+1: 0)
             tvxy (i,j) = tv
             tgxy (i,j) = tg
             zwtxy (i,j) = zwt
             waxy (i,j) = wa
             wtxy (i,j) = wt
             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
             emiss (i,j) = emissi
             eahxy (i,j) = eah
             tahxy (i,j) = tah
             fwetxy (i,j) = fwet
             sneqvoxy (i,j) = sneqvo
             alboldxy (i,j) = albold
             qsnowxy (i,j) = qsnow
             wslakexy (i,j) = wslake
             cmxy (i,j) = cm
             chxy (i,j) = chstar
             rssunxy (i,j) = rssun
             rsshaxy (i,j) = rssha
             bgapxy (i,j) = bgap
             wgapxy (i,j) = wgap
             gapxy (i,j) = gap
             tgvxy (i,j) = tgv
             tgbxy (i,j) = tgb
             chvxy (i,j) = chv
             chbxy (i,j) = chb
             runsfxy (i,j) = runsf
             runsbxy (i,j) = runsb
             ecanxy (i,j) = ecan
             edirxy (i,j) = esoil
             etranxy (i,j) = etran
             aparxy (i,j) = apar
             psnxy (i,j) = psn
             savxy (i,j) = sav
             sagxy (i,j) = sag
             fsnoxy (i,j) = fsno
             fsaxy (i,j) = fsa
             firaxy (i,j) = fira
             hfx (i,j) = fsh
             lh (i,j) = fcev + fgev + fctr
             grdflx (i,j) = ssoil
             tradxy (i,j) = trad
             tsxy (i,j) = ts
             neexy (i,j) = nee
             gppxy (i,j) = gpp
             nppxy (i,j) = npp
             fvegxy (i,j) = fveg
             t2mvxy (i,j) = t2mv
             t2mbxy (i,j) = t2mb
             q2mvxy (i,j) = q2mv
             q2mbxy (i,j) = q2mb
             chstarxy (i,j) = chstar
             chs (i,j) = chstar
             tstarxy (i,j) = tstar
             CANWAT(I,J) = canliqxy (i,j) + canicexy (i,j)
             IF ( SALB > -999 ) THEN
                ALBEDO(I,J) = salb
             ENDIF
             TSK(I,J) = tradxy (i,j)
             QSFC(I,J) = Q1/(1.0-Q1)
             q2mvxy(i,j) = q2mvxy(i,j)/(1.0-q2mvxy(i,j))
             q2mbxy(i,j) = q2mbxy(i,j)/(1.0-q2mbxy(i,j))
             SNOTIME(I,J) = SNOTIME1
             SMSTAV(I,J)=SOILW
             SMSTOT(I,J)=SOILM*1000.
             SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+runsfxy(i,j)*DT*1000.0
             UDRUNOFF(I,J)=UDRUNOFF(I,J)+runsbxy(i,j)*DT*1000.0
             IF(FFROZP.GT.0.5)THEN
                ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
             ENDIF
             IF(SNOW(I,J).GT.0.)THEN
             ENDIF
          ENDIF
       ENDDO
    ENDDO
  END SUBROUTINE noahmplsm
  SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP , &
       TSLB , SMOIS , SH2O , DZS , FNDSOILW , FNDSNOWH , &
       TSK, isnowxy , tvxy ,tgxy ,canicexy , &
       canliqxy ,eahxy ,tahxy ,cmxy ,chxy , &
       fwetxy ,sneqvoxy ,alboldxy ,qsnowxy ,wslakexy ,zwtxy ,waxy , &
       wtxy ,tsnoxy ,zsnsoxy ,snicexy ,snliqxy ,lfmassxy ,rtmassxy , &
       stmassxy ,woodxy ,stblcpxy ,fastcpxy ,xsaixy , &
       t2mvxy ,t2mbxy ,chstarxy , &
       num_soil_layers, restart, &
       allowed_to_read , &
       ids,ide, jds,jde, kds,kde, &
       ims,ime, jms,jme, kms,kme, &
       its,ite, jts,jte, kts,kte )
    INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
         & ims,ime, jms,jme, kms,kme, &
         & its,ite, jts,jte, kts,kte
    INTEGER, INTENT(IN) :: num_soil_layers
    LOGICAL, INTENT(IN) :: restart, &
         & allowed_to_read
    REAL, DIMENSION( num_soil_layers), INTENT(IN) :: DZS
    REAL, DIMENSION( ims:ime, num_soil_layers, 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
    LOGICAL, INTENT(IN) :: FNDSOILW, &
         & FNDSNOWH
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: TSK
    INTEGER, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: isnowxy
    REAL, DIMENSION(ims:ime,-2:num_soil_layers,jms:jme), INTENT(INOUT) :: zsnsoxy
    REAL, DIMENSION(ims:ime,-2: 0,jms:jme), INTENT(INOUT) :: tsnoxy
    REAL, DIMENSION(ims:ime,-2: 0,jms:jme), INTENT(INOUT) :: snicexy
    REAL, DIMENSION(ims:ime,-2: 0,jms:jme), INTENT(INOUT) :: snliqxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tvxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tgxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canicexy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canliqxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: eahxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tahxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: cmxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: chxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fwetxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: sneqvoxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: alboldxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: qsnowxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wslakexy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: zwtxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: waxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wtxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: lfmassxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: rtmassxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stmassxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: woodxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stblcpxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fastcpxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: xsaixy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: t2mvxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: t2mbxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: chstarxy
    REAL, DIMENSION(1:num_soil_layers) :: ZSOIL
    REAL :: BX, SMCMAX, PSISAT, FREE
    REAL, PARAMETER :: BLIM = 5.5
    REAL, PARAMETER :: HLICE = 3.335E5
    REAL, PARAMETER :: GRAV = 9.81
    REAL, PARAMETER :: T0 = 273.15
    INTEGER :: errflag
    character(len=80) :: err_message
    character(len=4) :: MMINSL
    character(len=*), intent(in) :: MMINLU
    MMINSL='STAS'
    call read_mp_veg_parameters(trim(MMINLU))
    IF ( allowed_to_read ) THEN
       CALL wrf_message( 'INITIALIZE THREE Noah LSM RELATED TABLES' )
       CALL SOIL_VEG_GEN_PARM( MMINLU, MMINSL )
    ENDIF
    IF( .NOT. restart ) THEN
       itf=min0(ite,ide-1)
       jtf=min0(jte,jde-1)
       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_fatal3("<stdin>",1163,&
"module_sf_noahlsm.F: lsminit: out of range value "
               "of ISLTYP. Is this field in the input?" )
       ENDIF
       DO J = jts , jtf
          DO I = its , itf
             BX = BB(ISLTYP(I,J))
             SMCMAX = MAXSMC(ISLTYP(I,J))
             PSISAT = SATPSI(ISLTYP(I,J))
             IF ( ( bx > 0.0 ) .AND. ( smcmax > 0.0 ) .AND. ( psisat > 0.0 ) ) THEN
                DO NS=1, num_soil_layers
                   IF ( TSLB(I,NS,J) < 273.149 ) THEN
                      BX = BB(ISLTYP(I,J))
                      SMCMAX = MAXSMC(ISLTYP(I,J))
                      PSISAT = SATPSI(ISLTYP(I,J))
                      IF ( BX > BLIM ) BX = BLIM
                      FK=(( (HLICE/(GRAV*(-PSISAT))) * &
                           ((TSLB(I,NS,J)-T0)/TSLB(I,NS,J)) )**(-1/BX) )*SMCMAX
                      FK = MAX(FK, 0.02)
                      SH2O(I,NS,J) = MIN( FK, SMOIS(I,NS,J) )
                      CALL FRH2O ( FREE , TSLB(I,NS,J) , SMOIS(I,NS,J) , SH2O(I,NS,J) )
                      SH2O(I,NS,J) = FREE
                   ELSE
                      SH2O(I,NS,J)=SMOIS(I,NS,J)
                   ENDIF
                END DO
             ELSE
                DO NS=1, num_soil_layers
                   SH2O(I,NS,J)=SMOIS(I,NS,J)
                END DO
             ENDIF
          ENDDO
       ENDDO
       IF(.NOT.FNDSNOWH)THEN
          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
             ENDDO
          ENDDO
       ENDIF
       DO J = jts,jtf
          DO I = its,itf
             tvxy (I,J) = TSK(I,J)
             tgxy (I,J) = TSK(I,J)
             CANWAT (I,J) = 0.0
             canliqxy (I,J) = CANWAT(I,J)
             canicexy (I,J) = 0.
             eahxy (I,J) = 2000.
             tahxy (I,J) = 287.
             t2mvxy (I,J) = TSK(I,J)
             t2mbxy (I,J) = TSK(I,J)
             chstarxy (I,J) = 0.0
             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
             waxy (I,J) = 4900.
             wtxy (I,J) = waxy(i,j)
             zwtxy (I,J) = (25. + 2.0) - waxy(i,j)/1000/0.2
             lfmassxy (I,J) = 50.
             stmassxy (I,J) = 50.0
             rtmassxy (I,J) = 500.0
             woodxy (I,J) = 500.0
             stblcpxy (I,J) = 1000.0
             fastcpxy (I,J) = 1000.0
             xsaixy (I,J) = 0.1
          enddo
       enddo
       ZSOIL(1) = -DZS(1)
       DO NS=2, num_soil_layers
          ZSOIL(NS) = ZSOIL(NS-1) - DZS(NS)
       END DO
       CALL snow_init ( ims , ime , jms , jme , its , itf , jts , jtf , 3 , &
            & num_soil_layers , zsoil , snow , tgxy , snowh , &
            & zsnsoxy , tsnoxy , snicexy , snliqxy , isnowxy )
    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 )
    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
    REAL, INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:NSOIL,jms:jme) :: ZSNSOXY
    REAL, INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1: 0,jms:jme) :: TSNOXY
    REAL, INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1: 0,jms:jme) :: SNICEXY
    REAL, INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1: 0,jms:jme) :: SNLIQXY
    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.35 ) ) 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.35 ) THEN
                ISNOWXY(I,J) = -3
                DZSNO(-2) = 0.05
                DZSNO(-1) = 0.10
                DZSNO( 0) = SNODEP(I,J) - DZSNO(-1) - DZSNO(-2)
             ELSE
                CALL wrf_error_fatal3("<stdin>",1354,&
"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)
             SNLIQXY(I,IZ,J) = 0.00
             SNICEXY(I,IZ,J) = 1.00 * DZSNO(IZ) * (SWE(I,J)/SNODEP(I,J))
          END DO
          DO IZ = ISNOWXY(I,J)+1 , 0
             DZSNSO(IZ) = -DZSNO(IZ)
          END DO
          DZSNSO(1) = ZSOIL(1)
          DO IZ = 2 , NSOIL
             DZSNSO(IZ) = (ZSOIL(IZ) - ZSOIL(IZ-1))
          END DO
          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
END MODULE module_sf_noahmpdrv