

MODULE module_sf_gfdl



CONTAINS


   SUBROUTINE SF_GFDL(U3D,V3D,T3D,QV3D,P3D, &
                     CP,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2, CPM, &
                     DT, SMOIS,num_soil_layers,ISLTYP,ZNT, &



                     UST,PSIM,PSIH, &
                     XLAND,HFX,QFX,TAUX,TAUY,LH,GSW,GLW,TSK,FLHC,FLQC, &
                     QGH,QSFC,U10,V10, &
                     GZ1OZ0,WSPD,BR,ISFFLX, &
                     EP1,EP2,KARMAN,NTSFLG,SFENTH, &
                     ids,ide, jds,jde, kds,kde, &
                     ims,ime, jms,jme, kms,kme, &
                     its,ite, jts,jte, kts,kte )

      USE MODULE_GFS_MACHINE, ONLY : kind_phys
      USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys,fpvs
      USE MODULE_GFS_PHYSCONS, grav => con_g

      IMPLICIT NONE
      character*255 :: message
      INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
                                        ims,ime, jms,jme, kms,kme, &
                                        its,ite, jts,jte, kts,kte, &
                                        ISFFLX,NUM_SOIL_LAYERS,NTSFLG
      REAL, INTENT(IN) :: &
                                        CP, &
                                        EP1, &
                                        EP2, &
                                        KARMAN, &
                                        R, &
                                        ROVCP, &
                                        DT, &
                                        SFENTH, &
                                        XLV
      REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
                                        P3D, &
                                        QV3D, &
                                        T3D, &
                                        U3D, &
                                        V3D
      INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: ISLTYP
      REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), INTENT(INOUT):: SMOIS
      REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
                                        PSFC, &
                                        GLW, &
                                        GSW, &
                                        XLAND
      REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
                                        TSK, &
                                        BR, &
                                        CHS, &
                                        CHS2, &
                                        CPM, &
                                        CQS2, &
                                        FLHC, &
                                        FLQC, &
                                        GZ1OZ0, &
                                        HFX, &
                                        LH, &
                                        PSIM, &
                                        PSIH, &
                                        QFX, &
                                        QGH, &
                                        QSFC, &
                                        UST, &
                                        ZNT, &
                                        WSPD, &
                                        TAUX, &
                                        TAUY
      REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
                                        U10, &
                                        V10
      REAL :: ESAT, &
                                        cpcgs, &
                                        smc, &
                                        smcdry, &
                                        smcmax
      REAL (kind=kind_phys) :: &
                                        RHOX
      REAL, DIMENSION(1:30) :: MAXSMC, &
                                        DRYSMC
      REAL (kind=kind_phys), DIMENSION(its:ite) :: &
                                        CH, &
                                        CM, &
                                        DDVEL, &
                                        DRAIN, &
                                        EP, &
                                        EVAP, &
                                        FH, &
                                        FH2, &
                                        FM, &
                                        HFLX, &
                                        PH, &
                                        PM, &
                                        PRSL1, &
                                        PRSLKI, &
                                        PS, &
                                        Q1, &
                                        Q2M, &
                                        QSS, &
                                        QSURF, &
                                        RB, &
                                        RCL, &
                                        RHO1, &
                                        SLIMSK, &
                                        STRESS, &
                                        T1, &
                                        T2M, &
                                        THGB, &
                                        THX, &
                                        TSKIN, &
                                        SHELEG, &
                                        U1, &
                                        U10M, &
                                        USTAR, &
                                        V1, &
                                        V10M, &
                                        WIND, &
                                        Z0RL, &
                                        Z1
      REAL, DIMENSION(kms:kme, ims:ime) :: &
                                        rpc, &
                                        tpc, &
                                        upc, &
                                        vpc
     REAL, DIMENSION(ims:ime) :: &
                                        pspc, &
                                        pkmax, &
                                        tstrc, &
                                        zoc, &
                                        mzoc, &
                                        wetc, &
                                        slwdc, &
                                        rib, &
                                        zkmax, &
                                        tkmax, &
                                        fxmx, &
                                        fxmy, &
                                        cdm, &
                                        fxh, &
                                        fxe, &
                                        xxfh, &
                                        xxfh2, &
                                        wind10, &
                                        tjloc
      INTEGER :: &
                                        I, &
                                        II, &
                                        IGPVS, &
                                        IM, &
                                        J, &
                                        K, &
                                        KM
        DATA MAXSMC/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, &
                    0.404, 0.464, 0.465, 0.406, 0.468, 0.468, &
                    0.439, 1.000, 0.200, 0.421, 0.000, 0.000, &
                    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
                    0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
        DATA DRYSMC/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, &
                    0.067, 0.120, 0.103, 0.100, 0.126, 0.138, &
                    0.066, 0.000, 0.006, 0.028, 0.000, 0.000, &
                    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
                    0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
      DATA IGPVS/0/
      save igpvs
   if(igpvs.eq.0) then
   endif
   igpvs=1
   IM=ITE-ITS+1
   KM=KTE-KTS+1
   WRITE(message,*)'WITHIN THE GFDL SCHEME, NTSFLG=1 FOR GFDL SLAB  2010 UPGRADS',NTSFLG
   call wrf_debug(1,message)
   DO J=jts,jte
      DO i=its,ite
        DDVEL(I)=0.
        RCL(i)=1.
        PRSL1(i)=P3D(i,kts,j)*.001
         wetc(i)=1.0
         if(xland(i,j).lt.1.99) then
         smc=smois(i,1,j)
         smcdry=drysmc(isltyp(i,j))
         smcmax=maxsmc(isltyp(i,j))
         wetc(i)=(smc-smcdry)/(smcmax-smcdry)
         wetc(i)=amin1(1.,amax1(wetc(i),0.))
         endif
        pspc(i)=PSFC(i,j)*10.
        pkmax(i)=P3D(i,kts,j)*10.
        PS(i)=PSFC(i,j)*.001
        Q1(I) = QV3D(i,kts,j)
        rpc(kts,i)=QV3D(i,kts,j)
        QSURF(I)=0.
        SHELEG(I)=0.
        SLIMSK(i)=ABS(XLAND(i,j)-2.)
        TSKIN(i)=TSK(i,j)
        tstrc(i)=TSK(i,j)
        T1(I) = T3D(i,kts,j)
        tpc(kts,i)=T3D(i,kts,j)
        U1(I) = U3D(i,kts,j)
        upc(kts,i)=U3D(i,kts,j) * 100.
        USTAR(I) = UST(i,j)
        V1(I) = V3D(i,kts,j)
        vpc(kts,i)=v3D(i,kts,j) * 100.
        Z0RL(I) = ZNT(i,j)*100.
        zoc(i)=ZNT(i,j)*100.
         if(XLAND(i,j).gt.1.99) zoc(i)=- zoc(i)
        slwdc(i)=gsw(i,j)+glw(i,j)
        slwdc(i)=0.239*60.*slwdc(i)*1.e-4
         tjloc(i)=float(j)
      ENDDO
      DO i=its,ite
         PRSLKI(i)=(PS(I)/PRSL1(I))**ROVCP
         THGB(I)=TSKIN(i)*(100./PS(I))**ROVCP
         THX(I)=T1(i)*(100./PRSL1(I))**ROVCP
         RHO1(I)=PRSL1(I)*1000./(R*T1(I)*(1.+EP1*Q1(I)))
         Q1(I)=Q1(I)/(1.+Q1(I))
      ENDDO
     CALL MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, &
                   pspc,pkmax,wetc,slwdc,tjloc, &
                   upc,vpc,tpc,rpc,dt,J,wind10,xxfh2,ntsflg,SFENTH, &
                   ids,ide, jds,jde, kds,kde, &
                   ims,ime, jms,jme, kms,kme, &
                   its,ite, jts,jte, kts,kte )
1010 format(2I4,9F11.6)
      DO i=its,ite
        IF(NTSFLG==1) then
        tsk(i,j) = tstrc(i)
        if(j.eq.jde) then
         tsk(i,j) = tsk(i,j-1)
         endif
        if(j.eq.jds) then
         tsk(i,j) = tsk(i,j+1)
         endif
        if(i.eq.ide) tsk(i,j) = tsk(i-1,j)
        if(i.eq.ids) tsk(i,j) = tsk(i+1,j)
        endif
        znt(i,j)= 0.01*abs(zoc(i))
        wspd(i,j) = SQRT(upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i))
        wspd(i,j) = amax1(wspd(i,j) ,100.)/100.
        u10m(i) = u1(i)*(wind10(i)/wspd(i,j))/100.
        v10m(i) = v1(i)*(wind10(i)/wspd(i,j))/100.
        zkmax(i) = -R*tpc(kts,i)*alog(pkmax(i)/pspc(i))/grav
        gz1oz0(i,j)=alog(zkmax(i)/znt(i,j))
        ustar (i)= 0.01*sqrt(cdm(i)* &
                   (upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i)))
        qfx (i,j)=-10.*fxe(i)
        hfx (i,j)= -10.*CP*fxh(i)
        taux (i,j)= fxmx(i)/10.
        tauy (i,j)= fxmy(i)/10.
        fm(i) = karman/sqrt(cdm(i))
        fh(i) = karman*xxfh(i)
        PSIM(i,j)=GZ1OZ0(i,j)-FM(i)
        PSIH(i,j)=GZ1OZ0(i,j)-FH(i)
        fh2(i) = karman*xxfh2(i)
        ch(i) = karman*karman/(fm(i) * fh(i))
        cm(i) = cdm(i)
        U10(i,j)=U10M(i)
        V10(i,j)=V10M(i)
        BR(i,j)=rib(i)
        CHS(I,J)=CH(I)*wspd (i,j)
        CHS2(I,J)=USTAR(I)*KARMAN/FH2(I)
        CPM(I,J)=CP*(1.+0.8*QV3D(i,kts,j))
        esat = fpvs(t1(i))
        QGH(I,J)=ep2*esat/(1000.*ps(i)-esat)
        esat = fpvs(tskin(i))
        qss(i) = ep2*esat/(1000.*ps(i)-esat)
        QSFC(I,J)=qss(i)
        UST(i,j)=ustar(i)
      ENDDO
      DO i=its,ite
        FLHC(i,j)=CPM(I,J)*RHO1(I)*CHS(I,J)
        FLQC(i,j)=RHO1(I)*CHS(I,J)
        CQS2(i,j)=CHS2(I,J)
      ENDDO
      IF (ISFFLX.EQ.0) THEN
        DO i=its,ite
          HFX(i,j)=0.
          LH(i,j)=0.
          QFX(i,j)=0.
        ENDDO
      ELSE
        DO i=its,ite
          IF(XLAND(I,J)-1.5.GT.0.)THEN
            hfx (i,j)= -10.*CP*fxh(i)
          ELSEIF(XLAND(I,J)-1.5.LT.0.)THEN
            hfx (i,j)= -10.*CP*fxh(i)
            HFX(I,J)=AMAX1(HFX(I,J),-250.)
          ENDIF
          qfx(i,j)=-10.*fxe(i)
          QFX(I,J)=AMAX1(QFX(I,J),0.)
          LH(I,J)=XLV*QFX(I,J)
        ENDDO
      ENDIF
    if(j.eq.jds.or.j.eq.jde) then
    write(message,*) "TSFC in gfdl sf mod,dt, its,ite,jts,jts", dt,its,ite,jts,jte,ids,ide,jds,jde
    call wrf_debug(1,message)
    write(message,*) "TSFC", (TSK(i,j),i=its,ite)
    call wrf_debug(1,message)
    endif
   ENDDO
   END SUBROUTINE SF_GFDL
       SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, &
                          pspc,pkmax,wetc,slwdc,tjloc, &
                          upc,vpc,tpc,rpc,dt,jfix,wind10,xxfh2,ntsflg,sfenth, &
                          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
      integer,intent(in) :: ims,ime, jms,jme, kms,kme
      integer,intent(in) :: its,ite, jts,jte, kts,kte
      integer,intent(in) :: jfix,ntsflg
      real, intent (out), dimension (ims :ime ) :: fxh
      real, intent (out), dimension (ims :ime ) :: fxe
      real, intent (out), dimension (ims :ime ) :: fxmx
      real, intent (out), dimension (ims :ime ) :: fxmy
      real, intent (out), dimension (ims :ime ) :: cdm
      real, intent (out), dimension (ims :ime ) :: rib
      real, intent (out), dimension (ims :ime ) :: xxfh
      real, intent (out), dimension (ims :ime ) :: xxfh2
      real, intent (out), dimension (ims :ime ) :: wind10
      real, intent ( inout), dimension (ims :ime ) :: zoc,mzoc
      real, intent ( inout), dimension (ims :ime ) :: tstrc
      real, intent ( in) :: dt
      real, intent ( in) :: sfenth
      real, intent ( in), dimension (ims :ime ) :: pspc
      real, intent ( in), dimension (ims :ime ) :: pkmax
      real, intent ( in), dimension (ims :ime ) :: wetc
      real, intent ( in), dimension (ims :ime ) :: slwdc
      real, intent ( in), dimension (ims :ime ) :: tjloc
      real, intent ( in), dimension (kms:kme, ims :ime ) :: upc
      real, intent ( in), dimension (kms:kme, ims :ime ) :: vpc
      real, intent ( in), dimension (kms:kme, ims :ime ) :: tpc
      real, intent ( in), dimension (kms:kme, ims :ime ) :: rpc
      integer, parameter :: icntx = 30
      integer, dimension(1 :ime) :: ifz
      integer, dimension(1 :ime) :: indx
      integer, dimension(1 :ime) :: istb
      integer, dimension(1 :ime) :: it
      integer, dimension(1 :ime) :: iutb
      real, dimension(1 :ime) :: aap
      real, dimension(1 :ime) :: bq1
      real, dimension(1 :ime) :: bq1p
      real, dimension(1 :ime) :: delsrad
      real, dimension(1 :ime) :: ecof
      real, dimension(1 :ime) :: ecofp
      real, dimension(1 :ime) :: estso
      real, dimension(1 :ime) :: estsop
      real, dimension(1 :ime) :: fmz1
      real, dimension(1 :ime) :: fmz10
      real, dimension(1 :ime) :: fmz2
      real, dimension(1 :ime) :: fmzo1
      real, dimension(1 :ime) :: foft
      real, dimension(1 :ime) :: foftm
      real, dimension(1 :ime) :: frac
      real, dimension(1 :ime) :: land
      real, dimension(1 :ime) :: pssp
      real, dimension(1 :ime) :: qf
      real, dimension(1 :ime) :: rdiff
      real, dimension(1 :ime) :: rho
      real, dimension(1 :ime) :: rkmaxp
      real, dimension(1 :ime) :: rstso
      real, dimension(1 :ime) :: rstsop
      real, dimension(1 :ime) :: sf10
      real, dimension(1 :ime) :: sf2
      real, dimension(1 :ime) :: sfm
      real, dimension(1 :ime) :: sfzo
      real, dimension(1 :ime) :: sgzm
      real, dimension(1 :ime) :: slwa
      real, dimension(1 :ime) :: szeta
      real, dimension(1 :ime) :: szetam
      real, dimension(1 :ime) :: t1
      real, dimension(1 :ime) :: t2
      real, dimension(1 :ime) :: tab1
      real, dimension(1 :ime) :: tab2
      real, dimension(1 :ime) :: tempa1
      real, dimension(1 :ime) :: tempa2
      real, dimension(1 :ime) :: theta
      real, dimension(1 :ime) :: thetap
      real, dimension(1 :ime) :: tsg
      real, dimension(1 :ime) :: tsm
      real, dimension(1 :ime) :: tsp
      real, dimension(1 :ime) :: tss
      real, dimension(1 :ime) :: ucom
      real, dimension(1 :ime) :: uf10
      real, dimension(1 :ime) :: uf2
      real, dimension(1 :ime) :: ufh
      real, dimension(1 :ime) :: ufm
      real, dimension(1 :ime) :: ufzo
      real, dimension(1 :ime) :: ugzm
      real, dimension(1 :ime) :: uzeta
      real, dimension(1 :ime) :: uzetam
      real, dimension(1 :ime) :: vcom
      real, dimension(1 :ime) :: vrtkx
      real, dimension(1 :ime) :: vrts
      real, dimension(1 :ime) :: wind
      real, dimension(1 :ime) :: windp
      real, dimension(1 :ime) :: xxfm
      real, dimension(1 :ime) :: xxsh
      real, dimension(1 :ime) :: z10
      real, dimension(1 :ime) :: z2
      real, dimension(1 :ime) :: zeta
      real, dimension(1 :ime) :: zkmax
      real, dimension(1 :ime) :: pss
      real, dimension(1 :ime) :: tstar
      real, dimension(1 :ime) :: ukmax
      real, dimension(1 :ime) :: vkmax
      real, dimension(1 :ime) :: tkmax
      real, dimension(1 :ime) :: rkmax
      real, dimension(1 :ime) :: zot
      real, dimension(1 :ime) :: fhzo1
      real, dimension(1 :ime) :: sfh
      real :: ux13, yo, y,xo,x,ux21,ugzzo,ux11,ux12,uzetao,xnum,alll
      real :: ux1,ugz,x10,uzo,uq,ux2,ux3,xtan,xden,y10,uzet1o,ugz10
      real :: szet2, zal2,ugz2
      real :: rovcp,boycon,cmo2,psps1,zog,enrca,rca,cmo1,amask,en,ca,a,c
      real :: sgz,zal10,szet10,fmz,szo,sq,fmzo,rzeta1,zal1g,szetao,rzeta2,zal2g
      real :: hcap,xks,pith,teps,diffot,delten,alevp,psps2,alfus,nstep
      real :: shfx,sigt4,reflect
      real :: cor1,cor2,szetho,zal2gh,cons_p000001,cons_7,vis,ustar,restar,rat
      real :: wndm,ckg
      real :: yz,y1,y2,y3,y4,windmks,znott,znotm
      integer:: i,j,ii,iq,nnest,icnt,ngd,ip
      real, dimension (223) :: tab
      real, dimension (223) :: table
      real, dimension (101) :: tab11
      real, dimension (41) :: table4
      real, dimension (42) :: tab3
      real, dimension (54) :: table2
      real, dimension (54) :: table3
      real, dimension (74) :: table1
      real, dimension (80) :: tab22
      equivalence (tab(1),tab11(1))
      equivalence (tab(102),tab22(1))
      equivalence (tab(182),tab3(1))
      equivalence (table(1),table1(1))
      equivalence (table(75),table2(1))
      equivalence (table(129),table3(1))
      equivalence (table(183),table4(1))
      data amask/ -98.0/
      data tab11/21*0.01403,0.01719,0.02101,0.02561,0.03117,0.03784, &
     &.04584,.05542,.06685,.08049,.09672,.1160,.1388,.1658,.1977,.2353, &
     &.2796,.3316,.3925,.4638,.5472,.6444,.7577,.8894,1.042,1.220,1.425, &
     &1.662,1.936,2.252,2.615,3.032,3.511,4.060,4.688,5.406,6.225,7.159, &
     &8.223,9.432,10.80,12.36,14.13,16.12,18.38,20.92,23.80,27.03,30.67, &
     &34.76,39.35,44.49,50.26,56.71,63.93,71.98,80.97,90.98,102.1,114.5, &
     &128.3,143.6,160.6,179.4,200.2,223.3,248.8,276.9,307.9,342.1,379.8, &
     &421.3,466.9,517.0,572.0,632.3,698.5,770.9,850.2,937.0,1032./
      data tab22/1146.6,1272.0,1408.1,1556.7,1716.9,1890.3,2077.6,2279.6 &
     &,2496.7,2729.8,2980.0,3247.8,3534.1,3839.8,4164.8,4510.5,4876.9, &
     &5265.1,5675.2,6107.8,6566.2,7054.7,7575.3,8129.4,8719.2,9346.5, &
     &10013.,10722.,11474.,12272.,13119.,14017.,14969.,15977.,17044., &
     &18173.,19367.,20630.,21964.,23373.,24861.,26430.,28086.,29831., &
     &31671.,33608.,35649.,37796.,40055.,42430.,44927.,47551.,50307., &
     &53200.,56236.,59422.,62762.,66264.,69934.,73777.,77802.,82015., &
     &86423.,91034.,95855.,100890.,106160.,111660.,117400.,123400., &
     &129650.,136170.,142980.,150070.,157460.,165160.,173180.,181530., &
     &190220.,199260./
      data tab3/208670.,218450.,228610.,239180.,250160.,261560.,273400., &
     &285700.,298450.,311690.,325420.,339650.,354410.,369710.,385560., &
     &401980.,418980.,436590.,454810.,473670.,493170.,513350.,534220., &
     &555800.,578090.,601130.,624940.,649530.,674920.,701130.,728190., &
     &756110.,784920.,814630.,845280.,876880.,909450.,943020.,977610., &
     &1013250.,1049940.,1087740./
      data table1/20*0.0,.3160e-02,.3820e-02,.4600e-02,.5560e-02,.6670e-02, &
     & .8000e-02,.9580e-02,.1143e-01,.1364e-01,.1623e-01,.1928e-01, &
     &.2280e-01,.2700e-01,.3190e-01,.3760e-01,.4430e-01,.5200e-01, &
     &.6090e-01,.7130e-01,.8340e-01,.9720e-01,.1133e+00,.1317e-00, &
     &.1526e-00,.1780e-00,.2050e-00,.2370e-00,.2740e-00,.3160e-00, &
     &.3630e-00,.4170e-00,.4790e-00,.5490e-00,.6280e-00,.7180e-00, &
     &.8190e-00,.9340e-00,.1064e+01,.1209e+01,.1368e+01,.1560e+01, &
     &.1770e+01,.1990e+01,.2260e+01,.2540e+01,.2880e+01,.3230e+01, &
     &.3640e+01,.4090e+01,.4590e+01,.5140e+01,.5770e+01,.6450e+01, &
     &.7220e+01/
      data table2/.8050e+01,.8990e+01,.1001e+02,.1112e+02,.1240e+02, &
     &.1380e+02,.1530e+02,.1700e+02,.1880e+02,.2080e+02,.2310e+02, &
     &.2550e+02,.2810e+02,.3100e+02,.3420e+02,.3770e+02,.4150e+02, &
     &.4560e+02,.5010e+02,.5500e+02,.6030e+02,.6620e+02,.7240e+02, &
     &.7930e+02,.8680e+02,.9500e+02,.1146e+03,.1254e+03,.1361e+03, &
     &.1486e+03,.1602e+03,.1734e+03,.1873e+03,.2020e+03,.2171e+03, &
     &.2331e+03,.2502e+03,.2678e+03,.2863e+03,.3057e+03,.3250e+03, &
     &.3457e+03,.3664e+03,.3882e+03,.4101e+03,.4326e+03,.4584e+03, &
     &.4885e+03,.5206e+03,.5541e+03,.5898e+03,.6273e+03,.6665e+03, &
     &.7090e+03/
      data table3/.7520e+03,.7980e+03,.8470e+03,.8980e+03,.9520e+03, &
     &.1008e+04,.1067e+04,.1129e+04,.1194e+04,.1263e+04,.1334e+04, &
     &.1409e+04,.1488e+04,.1569e+04,.1656e+04,.1745e+04,.1840e+04, &
     &.1937e+04,.2041e+04,.2147e+04,.2259e+04,.2375e+04,.2497e+04, &
     &.2624e+04,.2756e+04,.2893e+04,.3036e+04,.3186e+04,.3340e+04, &
     &.3502e+04,.3670e+04,.3843e+04,.4025e+04,.4213e+04,.4408e+04, &
     &.4611e+04,.4821e+04,.5035e+04,.5270e+04,.5500e+04,.5740e+04, &
     &.6000e+04,.6250e+04,.6520e+04,.6810e+04,.7090e+04,.7390e+04, &
     &.7700e+04,.8020e+04,.8350e+04,.8690e+04,.9040e+04,.9410e+04, &
     &.9780e+04/
      data table4/.1016e+05,.1057e+05,.1098e+05,.1140e+05,.1184e+05, &
     &.1230e+05,.1275e+05,.1324e+05,.1373e+05,.1423e+05,.1476e+05, &
     &.1530e+05,.1585e+05,.1642e+05,.1700e+05,.1761e+05,.1822e+05, &
     &.1886e+05,.1950e+05,.2018e+05,.2087e+05,.2158e+05,.2229e+05, &
     &.2304e+05,.2381e+05,.2459e+05,.2539e+05,.2621e+05,.2706e+05, &
     &.2792e+05,.2881e+05,.2971e+05,.3065e+05,.3160e+05,.3257e+05, &
     &.3357e+05,.3459e+05,.3564e+05,.3669e+05,.3780e+05,.0000e+00/
      real,parameter :: cp = 1.00464e7
      real,parameter :: g = 980.6
      real,parameter :: rgas = 2.87e6
      real,parameter :: og = 1./g
      character*255 :: message
      j = IFIX(tjloc(its))
       cor1 = .120
       cor2 = 720.
       yz= -0.0001344
       y1= 3.015e-05
       y2= 1.517e-06
       y3= -3.567e-08
       y4= 2.046e-10
      do i = its,ite
        z10(i) = 1000.
        z2 (i) = 200.
        pss(i) = pspc(i)
        tstar(i) = tstrc(i)
        ukmax(i) = upc(1,i)
        vkmax(i) = vpc(1,i)
        tkmax(i) = tpc(1,i)
        rkmax(i) = rpc(1,i)
      enddo
      do i = its,ite
        windp(i) = SQRT(ukmax(i)*ukmax(i) + vkmax(i)*vkmax(i))
        wind (i) = amax1(windp(i),100.)
        if (zoc(i) .LT. amask) zoc(i) = -0.0185*0.001*wind(i)*wind(i)*og
        if (zoc(i) .GT. 0.0) then
          ecof(i) = wetc(i)
          land(i) = 1.0
          zot (i) = zoc(i)
        else
          ecof(i) = wetc(i)
          land(i) = 0.0
          zot (i) = zoc(i)
          windmks=wind(i)*.01
          znott=0.2375*exp(-0.5250*windmks) + 0.0025*exp(-0.0211*windmks)
          znott=0.01*znott
          if(windmks.le. 7.) &
          znott = (0.0185/9.8*(7.59e-8*wind(i)**2+ &
                                  2.46e-4*wind(i))**2)
          zot (i) = 100.*znott
              zoc(i)=-(0.0185/9.8*(7.59e-8*wind(i)**2+ &
                                  2.46e-4*wind(i))**2)*100.
            if(wind(i).ge.1250.0) &
             zoc(i)=-(.000739793 * wind(i) -0.58)/10
        if(wind(i).ge.3000.) then
        windmks=wind(i)*.01
        znotm = yz +windmks*y1 +windmks**2*y2 +windmks**3*y3 +windmks**4*y4
        zoc(i) = 100.*znotm
        endif
        endif
      mzoc(i) = zoc(i)
      enddo
      en = 2.
      c = .76
      a = 5.
      ca = .4
      cmo1 = .5*a - 1.648
      cmo2 = 17.193 + .5*a - 10.*c
      boycon = .61
        rovcp=rgas/cp
      do i = its,ite
        theta(i) = tkmax(i)/((pkmax(i)/pspc(i))**rovcp)
        vrtkx(i) = 1.0 + boycon*rkmax(i)
        zkmax(i) = -rgas*tkmax(i)*alog(pkmax(i)/pspc(i))*og
      enddo
      do i = its,ite
        tsg (i) = tstar(i)
        tab1 (i) = tstar(i) - 153.16
        it (i) = IFIX(tab1(i))
        tab2 (i) = tab1(i) - FLOAT(it(i))
        t1 (i) = tab(it(i) + 1)
        t2 (i) = table(it(i) + 1)
        estso(i) = t1(i) + tab2(i)*t2(i)
         psps1 = (pss(i) - estso(i))
          if(psps1 .EQ. 0.0)then
           psps1 = .1
          endif
        rstso(i) = 0.622*estso(i)/psps1
        vrts (i) = 1. + boycon*ecof(i)*rstso(i)
      enddo
      do i = its,ite
        tempa1(i) = theta(i)*vrtkx(i) - tstar(i)*vrts(i)
        tempa2(i) = tempa1(i)*(theta(i) - tstar(i))
        if (tempa2(i) .LT. 0.) tempa1(i) = 1.0e-4
        tab1(i) = ABS(tempa1(i))
        if (tab1(i) .LT. 1.0e-4) tempa1(i) = 1.0e-4
        rib (i) = g*zkmax(i)*tempa1(i)/ &
                                    (tkmax(i)*vrtkx(i)*wind(i)*wind(i))
        tab2(i) = ABS(zoc(i))
        tab1(i) = 0.95/(c*(1. - tab2(i)/zkmax(i)))
        if (rib(i) .GT. tab1(i)) rib(i) = tab1(i)
      enddo
      do i = its,ite
        zeta(i) = ca*rib(i)/0.03
      enddo
      rca = 1./ca
      enrca = en*rca
      enrca = 0.0
      zog = .0185*og
      ip = 0
      do i = its,ite
        if (zeta(i) .GE. 0.0) then
          ip = ip + 1
          istb(ip) = i
        endif
      enddo
      if (ip .EQ. 0) go to 170
      do i = 1,ip
        szetam(i) = 1.0e+30
        sgzm(i) = 0.0e+00
        szeta(i) = zeta(istb(i))
        ifz(i) = 1
      enddo
      do icnt = 1,icntx
        do i = 1,ip
          if (ifz(i) .EQ. 0) go to 80
          zal1g = ALOG(szeta(i))
          if (szeta(i) .LE. 0.5) then
            fmz1(i) = (zal1g + a*szeta(i))*rca
          else if (szeta(i) .GT. 0.5 .AND. szeta(i) .LE. 10.) then
            rzeta1 = 1./szeta(i)
            fmz1(i) = (8.*zal1g + 4.25*rzeta1 - &
                                          0.5*rzeta1*rzeta1 + cmo1)*rca
          else if (szeta(i) .GT. 10.) then
            fmz1(i) = (c*szeta(i) + cmo2)*rca
          endif
          szetao = ABS(zoc(istb(i)))/zkmax(istb(i))*szeta(i)
          zal2g = ALOG(szetao)
          if (szetao .LE. 0.5) then
            fmzo1(i) = (zal2g + a*szetao)*rca
            sfzo (i) = 1. + a*szetao
          else if (szetao .GT. 0.5 .AND. szetao .LE. 10.) then
            rzeta2 = 1./szetao
            fmzo1(i) = (8.*zal2g + 4.25*rzeta2 - &
                                          0.5*rzeta2*rzeta2 + cmo1)*rca
            sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
          else if (szetao .GT. 10.) then
            fmzo1(i) = (c*szetao + cmo2)*rca
            sfzo (i) = c*szetao
          endif
          szetho = ABS(zot(istb(i)))/zkmax(istb(i))*szeta(i)
          zal2gh = ALOG(szetho)
          if (szetho .LE. 0.5) then
            fhzo1(i) = (zal2gh + a*szetho)*rca
            sfzo (i) = 1. + a*szetho
          else if (szetho .GT. 0.5 .AND. szetho .LE. 10.) then
            rzeta2 = 1./szetho
            fhzo1(i) = (8.*zal2gh + 4.25*rzeta2 - &
                                          0.5*rzeta2*rzeta2 + cmo1)*rca
            sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
          else if (szetho .GT. 10.) then
            fhzo1(i) = (c*szetho + cmo2)*rca
            sfzo (i) = c*szetho
          endif
            szet10 = ABS(z10(istb(i)))/zkmax(istb(i))*szeta(i)
            zal10 = ALOG(szet10)
            if (szet10 .LE. 0.5) then
              fmz10(i) = (zal10 + a*szet10)*rca
            else if (szet10 .GT. 0.5 .AND. szet10 .LE. 10.) then
              rzeta2 = 1./szet10
              fmz10(i) = (8.*zal10 + 4.25*rzeta2 - &
                                          0.5*rzeta2*rzeta2 + cmo1)*rca
            else if (szet10 .GT. 10.) then
              fmz10(i) = (c*szet10 + cmo2)*rca
            endif
            sf10(i) = fmz10(i) - fmzo1(i)
            szet2 = ABS(z2 (istb(i)))/zkmax(istb(i))*szeta(i)
            zal2 = ALOG(szet2 )
            if (szet2 .LE. 0.5) then
              fmz2 (i) = (zal2 + a*szet2 )*rca
            else if (szet2 .GT. 0.5 .AND. szet2 .LE. 2.) then
              rzeta2 = 1./szet2
              fmz2 (i) = (8.*zal2 + 4.25*rzeta2 - &
                                          0.5*rzeta2*rzeta2 + cmo1)*rca
            else if (szet2 .GT. 2.) then
              fmz2 (i) = (c*szet2 + cmo2)*rca
            endif
            sf2 (i) = fmz2 (i) - fmzo1(i)
          sfm(i) = fmz1(i) - fmzo1(i)
          sfh(i) = fmz1(i) - fhzo1(i)
          sgz = ca*rib(istb(i))*sfm(i)*sfm(i)/ &
                                               (sfh(i) + enrca*sfzo(i))
          fmz = (sgz - szeta(i))/szeta(i)
          fmzo = ABS(fmz)
          if (fmzo .GE. 5.0e-5) then
            sq = (sgz - sgzm(i))/(szeta(i) - szetam(i))
            if(sq .EQ. 1) then
             write(message,*)'NCO ERROR DIVIDE BY ZERO IN MFLUX2 (STABLE CASE)'
             call wrf_message(message)
             write(message,*)'sq is 1 ',fmzo,sgz,sgzm(i),szeta(i),szetam(i)
             call wrf_message(message)
            endif
            szetam(i) = szeta(i)
            szeta (i) = (sgz - szeta(i)*sq)/(1.0 - sq)
            sgzm (i) = sgz
          else
            ifz(i) = 0
          endif
80 continue
        enddo
      enddo
      do i = 1,ip
        if (ifz(i) .GE. 1) go to 110
      enddo
      go to 130
110 continue
      write(6,120)
120 format(2X, ' NON-CONVERGENCE FOR STABLE ZETA IN ROW ')
130 continue
      do i = 1,ip
        szo = zoc(istb(i))
        if (szo .LT. 0.0) then
        wndm=wind(istb(i))*0.01
          if(wndm.lt.15.0) then
          ckg=0.0185*og
          else
           ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og
          endif
       szo = - ckg*wind(istb(i))*wind(istb(i))/ &
                             (sfm(i)*sfm(i))
        cons_p000001 = .000001
        cons_7 = 7.
        vis = 1.4E-1
        ustar = sqrt( -szo / zog)
        restar = -ustar * szo / vis
        restar = max(restar,cons_p000001)
        rat = 2.67 * restar ** .25 - 2.57
        rat = min(rat ,cons_7)
 rat=0.
        zot(istb(i)) = szo * exp(-rat)
         else
         zot(istb(i)) = zoc(istb(i))
        endif
             zoc(istb(i)) = szo
      enddo
      do i = 1,ip
        xxfm(istb(i)) = sfm(i)
        xxfh(istb(i)) = sfh(i)
        xxfh2(istb(i)) = sf2 (i)
        xxsh(istb(i)) = sfzo(i)
      enddo
        do i = 1,ip
         wind10(istb(i)) = sf10(i)*wind(istb(i))/sfm(i)
         wind10(istb(i)) = wind10(istb(i)) * 1.944
           if(wind10(istb(i)) .GT. 6000.0) then
       wind10(istb(i))=wind10(istb(i))+wind10(istb(i))*cor1 &
                        - cor2
           endif
         wind10(istb(i)) = wind10(istb(i)) / 1.944
        enddo
170 continue
      iq = 0
      do i = its,ite
        if (zeta(i) .LT. 0.0) then
          iq = iq + 1
          iutb(iq) = i
        endif
      enddo
      if (iq .EQ. 0) go to 290
      do i = 1,iq
        uzeta (i) = zeta(iutb(i))
        ifz (i) = 1
        uzetam(i) = 1.0e+30
        ugzm (i) = 0.0e+00
      enddo
      do icnt = 1,icntx
        do i = 1,iq
          if (ifz(i) .EQ. 0) go to 200
          ugzzo = ALOG(zkmax(iutb(i))/ABS(zot(iutb(i))))
          uzetao = ABS(zot(iutb(i)))/zkmax(iutb(i))*uzeta(i)
          ux11 = 1. - 16.*uzeta(i)
          ux12 = 1. - 16.*uzetao
          y = SQRT(ux11)
          yo = SQRT(ux12)
          ufzo(i) = 1./yo
          ux13 = (1. + y)/(1. + yo)
          ux21 = ALOG(ux13)
          ufh(i) = (ugzzo - 2.*ux21)*rca
          ugzzo = ALOG(zkmax(iutb(i))/ABS(zoc(iutb(i))))
          uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
          ux11 = 1. - 16.*uzeta(i)
          ux12 = 1. - 16.*uzetao
          y = SQRT(ux11)
          yo = SQRT(ux12)
          ux13 = (1. + y)/(1. + yo)
          ux21 = ALOG(ux13)
          x = SQRT(y)
          xo = SQRT(yo)
          xnum = (x**2 + 1.)*((x + 1.)**2)
          xden = (xo**2 + 1.)*((xo + 1.)**2)
          xtan = ATAN(x) - ATAN(xo)
          ux3 = ALOG(xnum/xden)
          ufm(i) = (ugzzo - ux3 + 2.*xtan)*rca
            ugz10 = ALOG(z10(iutb(i))/ABS(zoc(iutb(i))))
            uzet1o = ABS(z10(iutb(i)))/zkmax(iutb(i))*uzeta(i)
            uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
            ux11 = 1. - 16.*uzet1o
            ux12 = 1. - 16.*uzetao
            y = SQRT(ux11)
            y10 = SQRT(ux12)
            ux13 = (1. + y)/(1. + y10)
            ux21 = ALOG(ux13)
            x = SQRT(y)
            x10 = SQRT(y10)
            xnum = (x**2 + 1.)*((x + 1.)**2)
            xden = (x10**2 + 1.)*((x10 + 1.)**2)
            xtan = ATAN(x) - ATAN(x10)
            ux3 = ALOG(xnum/xden)
            uf10(i) = (ugz10 - ux3 + 2.*xtan)*rca
          ugz2 = ALOG(z2 (iutb(i))/ABS(zoc(iutb(i))))
          uzet1o = ABS(z2 (iutb(i)))/zkmax(iutb(i))*uzeta(i)
          uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
          ux11 = 1. - 16.*uzet1o
          ux12 = 1. - 16.*uzetao
          y = SQRT(ux11)
          yo = SQRT(ux12)
          ux13 = (1. + y)/(1. + yo)
          ux21 = ALOG(ux13)
          uf2 (i) = (ugzzo - 2.*ux21)*rca
          ugz = ca*rib(iutb(i))*ufm(i)*ufm(i)/(ufh(i) + enrca*ufzo(i))
          ux1 = (ugz - uzeta(i))/uzeta(i)
          ux2 = ABS(ux1)
          if (ux2 .GE. 5.0e-5) then
            uq = (ugz - ugzm(i))/(uzeta(i) - uzetam(i))
            uzetam(i) = uzeta(i)
            if(uq .EQ. 1) then
             write(0,*)'NCO ERROR DIVIDE BY ZERO IN MFLUX2 (UNSTABLE CASE)'
             write(0,*)'uq is 1 ',ux2,ugz,ugzm(i),uzeta(i),uzetam(i)
            endif
            uzeta (i) = (ugz - uzeta(i)*uq)/(1.0 - uq)
            ugzm (i) = ugz
          else
            ifz(i) = 0
          endif
200 continue
        enddo
      enddo
      do i = 1,iq
        if (ifz(i) .GE. 1) go to 230
      enddo
      go to 250
230 continue
      write(6,240)
240 format(2X, ' NON-CONVERGENCE FOR UNSTABLE ZETA IN ROW ')
250 continue
      do i = 1,iq
        uzo = zoc(iutb(i))
        if (zoc(iutb(i)) .LT. 0.0) then
        wndm=wind(iutb(i))*0.01
         if(wndm.lt.15.0) then
         ckg=0.0185*og
         else
            ckg=(4*0.000308*wndm)*og
            ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og
         endif
        uzo =-ckg*wind(iutb(i))*wind(iutb(i))/(ufm(i)*ufm(i))
        cons_p000001 = .000001
        cons_7 = 7.
        vis = 1.4E-1
        ustar = sqrt( -uzo / zog)
        restar = -ustar * uzo / vis
        restar = max(restar,cons_p000001)
        rat = 2.67 * restar ** .25 - 2.57
        rat = min(rat ,cons_7)
 rat=0.0
        zot(iutb(i)) = uzo * exp(-rat)
         else
         zot(iutb(i)) = zoc(iutb(i))
        endif
            zoc(iutb(i)) = uzo
      enddo
        do i = 1,iq
          wind10(iutb(i)) = uf10(i)*wind(iutb(i))/ufm(i)
          wind10(iutb(i)) = wind10(iutb(i)) * 1.944
           if(wind10(iutb(i)) .GT. 6000.0) then
         wind10(iutb(i))=wind10(iutb(i))+wind10(iutb(i))*cor1 &
                         - cor2
           endif
          wind10(iutb(i)) = wind10(iutb(i)) / 1.944
        enddo
      do i = 1,iq
        xxfm(iutb(i)) = ufm(i)
        xxfh(iutb(i)) = ufh(i)
        xxfh2(iutb(i)) = uf2 (i)
        xxsh(iutb(i)) = ufzo(i)
      enddo
290 continue
      do i = its,ite
        ucom(i) = ukmax(i)
        vcom(i) = vkmax(i)
        if (windp(i) .EQ. 0.0) then
          windp(i) = 100.0
          ucom (i) = 100.0/SQRT(2.0)
          vcom (i) = 100.0/SQRT(2.0)
        endif
        rho(i) = pss(i)/(rgas*(tsg(i) + enrca*(theta(i) - &
                tsg(i))*xxsh(i)/(xxfh(i) + enrca*xxsh(i))))
        bq1(i) = wind(i)*rho(i)/(xxfm(i)*(xxfh(i) + enrca*xxsh(i)))
      enddo
      if (ntsflg .EQ. 0) go to 370
      alll = 600.
      xks = 0.01
      hcap = .5/2.39e-8
      pith = SQRT(4.*ATAN(1.0))
      alfus = alll/2.39e-8
      teps = 0.1
      ip = 0
      do i = its,ite
        if (land(i) .EQ. 1) then
          ip = ip + 1
          indx (ip) = i
          slwa (ip) = slwdc(i)/(2.39e-8*60.)
          tss (ip) = tstar(i)
          thetap (ip) = theta(i)
          rkmaxp (ip) = rkmax(i)
          aap (ip) = 5.673e-5
          pssp (ip) = pss(i)
          ecofp (ip) = ecof(i)
          estsop (ip) = estso(i)
          rstsop (ip) = rstso(i)
          bq1p (ip) = bq1(i)
          bq1p (ip) = amax1(bq1p(ip),0.1e-3)
          delsrad(ip) = dt *pith/(hcap*SQRT(3600.*24.*xks))
        endif
      enddo
      do i = 1,ip
        ifz (i) = 1
        tsm (i) = tss(i)
        rdiff(i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
300 format(2X, ' SURFACE EQUILIBRIUM CALCULATION ')
          foftm(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsm(i)**4 - &
           cp*bq1p(i)*(tsm(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* &
           rdiff(i))
        tsp(i) = foftm(i)
      enddo
      do icnt = 1,icntx
        do i = 1,ip
          if (ifz(i) .EQ. 0) go to 330
          tab1 (i) = tsp(i) - 153.16
          it (i) = IFIX(tab1(i))
          tab2 (i) = tab1(i) - FLOAT(it(i))
          t1 (i) = tab(it(i) + 1)
          t2 (i) = table(it(i) + 1)
          estsop(i) = t1(i) + tab2(i)*t2(i)
             psps2 = (pssp(i) - estsop(i))
             if(psps2 .EQ. 0.0)then
               psps2 = .1
             endif
          rstsop(i) = 0.622*estsop(i)/psps2
          rdiff (i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
            foft(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsp(i)**4 - &
             cp*bq1p(i)*(tsp(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* &
             rdiff(i))
          frac(i) = ABS((foft(i) - tsp(i))/tsp(i))
          if (frac(i) .GE. teps) then
            qf (i) = (foft(i) - foftm(i))/(tsp(i) - tsm(i))
            tsm (i) = tsp(i)
            tsp (i) = (foft(i) - tsp(i)*qf(i))/(1. - qf(i))
            foftm(i) = foft(i)
          else
            ifz(i) = 0
          endif
330 continue
        enddo
      enddo
      do i = 1,ip
        if (ifz(i) .EQ. 1) then
        write(6, 340) tsp(i), i, j
340 format(2X, ' NON-CONVERGENCE OF T* PREDICTED (T*,I,J) = ', E14.8, &
            2I5)
        write(6,345) indx(i), j, tstar(indx(i)), tsp(i), ip
345 format(2X, ' I, J, OLD T*, NEW T*, NPTS ', 2I5, 2E14.8, I5)
350 format(2X, ' REFLECT, SIGT4, SHFX, ALEVP, DELTEN, DIFFOT ', &
            6E14.8)
        endif
      enddo
      do i = 1,ip
        ii = indx(i)
        tstrc(ii) = tsp (i)
      enddo
370 continue
      do i = its,ite
        fxh(i) = bq1(i)*(theta(i) - tsg(i))
        fxe(i) = ecof(i)*bq1(i)*(rkmax(i) - rstso(i))
        if (fxe(i) .GT. 0.0) fxe(i) = 0.0
        fxmx(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*ucom(i)/ &
                 windp(i)
        fxmy(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*vcom(i)/ &
        windp(i)
        cdm(i) = 1./(xxfm(i)*xxfm(i))
      enddo
      return
      end subroutine MFLUX2
  SUBROUTINE hwrfsfcinit(isn,XICE,VEGFRA,SNOW,SNOWC,CANWAT,SMSTAV, &
                        SMSTOT, SFCRUNOFF,UDRUNOFF,GRDFLX,ACSNOW, &
                        ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,DZS,SFCEVP, &
                        TMN, &
                        num_soil_layers, &
                        allowed_to_read, &
                        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
   INTEGER, INTENT(IN) :: num_soil_layers
   REAL, DIMENSION( num_soil_layers), INTENT(IN) :: DZS
   REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
            INTENT(INOUT) :: SMOIS, &
                                                         TSLB
   REAL, DIMENSION( ims:ime, jms:jme ) , &
            INTENT(INOUT) :: SNOW, &
                                                         SNOWC, &
                                                        CANWAT, &
                                                        SMSTAV, &
                                                        SMSTOT, &
                                                     SFCRUNOFF, &
                                                      UDRUNOFF, &
                                                        SFCEVP, &
                                                        GRDFLX, &
                                                        ACSNOW, &
                                                          XICE, &
                                                        VEGFRA, &
                                                        TMN, &
                                                        ACSNOM
   INTEGER, DIMENSION( ims:ime, jms:jme ) , &
            INTENT(INOUT) :: IVGTYP, &
                                                        ISLTYP
  INTEGER, INTENT(IN) :: isn
  LOGICAL, INTENT(IN) :: allowed_to_read
  INTEGER :: iseason
  INTEGER :: icm,jcm,itf,jtf
  INTEGER :: I,J,L
   itf=min0(ite,ide-1)
   jtf=min0(jte,jde-1)
   icm = ide/2
   jcm = jde/2
   iseason=isn
   DO J=jts,jtf
       DO I=its,itf
       SNOWC(i,j)=0.
    ENDDO
   ENDDO
  END SUBROUTINE hwrfsfcinit
END MODULE module_sf_gfdl
