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