subroutine sfceta2(sm,hlam,hphi,sfcgridhi,smc,stc,im,jm,nm)

      implicit none

      integer, intent(in) :: im,jm,nm

      include 'const.h'
      real,dimension(0:im+1,0:jm+1,nm)::sst,epsr,si,sice,sst2    &
                ,vegfrc,sm,hlam,hphi,sno,albedo,albase,mxsnal
      integer::i,j,n,ndays
      integer,dimension(0:im+1,0:jm+1,nm)::ivgtyp,isltyp,islope
      real,dimension(0:im+1,0:jm+1,nm)::dum
      integer*4,dimension(360,180)::i1d
      real,dimension(0:im+1,0:jm+1,nm,14)::sfcgridhi
      real,dimension(0:im+1,0:jm+1,nm,4)::smc,stc,sh2o
      real::hphip
      real,parameter::rtd=180.0/3.1415926
      character(len=100)::fname1
      character(len=3)::cn

      lvsst=.false.   !dragan

      if(sstc.eq.1)then
        call ssthires(sst,sm,hphi,hlam,idat,im,jm,nm,0)
      else
        call ssthires(sst2,sm,hphi,hlam,idat,im,jm,nm,0)
        call ssthires2(sst,sm,hphi,hlam,idat,im,jm,nm,0)

        do n=1,nm
        do j=1,jm
        do i=1,im
          if(sst(i,j,n).gt.500.)then
            sst(i,j,n)=sst2(i,j,n)
          endif
        enddo
        enddo
        enddo
      endif
        call snohires(si,sm,hphi,hlam,6,.true.,im,jm,nm)
        call vfrac(hphi,hlam,sm,si,vegfrc,im,jm,nm)

        sno=0.

        do n=1,nm
        do i=1,im
        do j=1,jm
          hphip=hphi(i,j,n)*rtd

          if(hphip.gt.20)then
            if(sm(i,j,n).gt.0.9)then
              epsr(i,j,n)=0.97
              if(si(i,j,n).gt.0)then
                sm(i,j,n)=0.
                si(i,j,n)=0.
                sice(i,j,n)=1.
              endif
            else
              epsr(i,j,n)=1.0
              sice(i,j,n)=0.
              sno(i,j,n)=si(i,j,n)*0.2
            endif
          else                                     ! use low res analysis data
            if(sm(i,j,n).gt.0.9)then
              epsr(i,j,n)=0.97
              if(sfcgridhi(i,j,n,14).gt.0)then
                sm(i,j,n)=0.
                si(i,j,n)=0.
                sice(i,j,n)=1.
              endif
            else
              epsr(i,j,n)=1.0
              sice(i,j,n)=0.
              sno(i,j,n)=sfcgridhi(i,j,n,13)/200*0.2  ! assume snow density is 200 kg/m^3
              si(i,j,n)=sno(i,j,n)/0.2

              if(sno(i,j,n).lt.0)then
                sno(i,j,n)=0.
                si(i,j,n)=0.
              endif
            endif
          endif
    
        enddo
        enddo
        enddo


       rewind 30
       read (30) i1d
       ivgtyp=7
       call putem(hphi,hlam,i1d,sm,si,ivgtyp,im,jm,nm)

       rewind 31
       read (31) i1d
       isltyp=2
       call putem(hphi,hlam,i1d,sm,si,isltyp,im,jm,nm)

       do n=1,nm
       do i=1,im
       do j=1,jm
         if(isltyp(i,j,n).eq.13) isltyp(i,j,n)=9
         if(isltyp(i,j,n).eq.0) then
            isltyp(i,j,n)=9
         end if
       enddo
       enddo
       enddo

       call sfch2o(sm,stc,smc,isltyp,sh2o)

       rewind 32
       read (32) i1d
       isltyp=1
       call putem(hphi,hlam,i1d,sm,si,islope,im,jm,nm)

       call albsfc(sm,sice,sno,vegfrc,hphi,hlam,albedo,albase  &
                   ,mxsnal)

       call bocoh(sm,im,jm,nm)
       call bocoh(epsr,im,jm,nm)
       call bocoh(sice,im,jm,nm)
       call bocoh(vegfrc,im,jm,nm)
       dum=ivgtyp
       call bocoh(dum,im,jm,nm)
       ivgtyp=dum
       dum=isltyp
       call bocoh(dum,im,jm,nm)
       isltyp=dum
       dum=islope
       call bocoh(dum,im,jm,nm)
       islope=dum
       call bocoh(albedo,im,jm,nm)
       call bocoh(mxsnal,im,jm,nm)

       open(1,file='sfceta2.dat',form='unformatted')
       write(1)sst,epsr,sm,sice,si,sno,vegfrc,ivgtyp  &
               ,isltyp,islope,sh2o,albedo,albase,mxsnal
       close(1)

       if(lvsst)then
       ndays=ntstm*dt/(3600*24)

       do n=1,ndays
        call ssthires(sst,sm,hphi,hlam,idat,im,jm,nm,n)
        write(cn,'(i3.3)') n
        fname1='sst'//cn//'.dat'
        open(1,file=fname1,form='unformatted')
        write(1)sst
        close(1)
       enddo
       endif

       end subroutine sfceta2