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,sice_2,sst2    &
                ,vegfrc,sm,hlam,hphi,sno,albedo,albase,mxsnal
      integer::i,j,n,ndays, nmonths, t ,INDAY,IOUTNUM,IIOUT,IHRST
      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  !!!!!,msstt(im+1,jm+1,nm)
      real,parameter::rtd=180.0/3.1415926
      character(len=100)::fname1
!      character(len=3)::cn!!!!!!!!!!!
      character(len=4)::cn  !dragan
     
!!!!!!!D july 2015      logical lusaf

      lvsst=.true. !DRAGAN, may 2015
!!!      lvsst=.false.
! print*,'before ssthires'

!GSM
      NAMELIST /CNSTDATA/INDAY,IHRST,IDAT,IOUTNUM,IIOUT

      open(19,file='cnstdata',status='old')
      READ(19,CNSTDATA)
!GSM

      if(sstc.eq.1)then

        write(0,*)' in sfcEta2 ', sstc,' before ssthires'
        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).eq.9.9990003E+20.or.sst(i,j,n).ge.350)then
       print*,'SST>350',sst(i,j,n), i,j,n
         endif
	 
	 
	  if(sst(i,j,n).gt.500.)then
	    sst(i,j,n)=sst2(i,j,n)
          endif
        enddo
        enddo
        enddo
      endif
!print*,'after ssthires'

!!!july2015      call snohires(si,sm,hphi,hlam,6,lusaf,im,jm,nm)
      call snohires(si,sm,hphi,hlam,6,.true.,im,jm,nm)       
!      print*,'SI=',si(326,112,1:38)
! print*,'after snohires'     
      do n=1,nm
      do i=1,im
      do j=1,jm
!!!        si(i,j,n)=sfcgridhi(i,j,n,14)*0.2
!!!        sice_2(i,j,n)=sfcgridhi(i,j,n,13)
      enddo
      enddo
      enddo
 
        call vfrac(hphi,hlam,sm,si,vegfrc,im,jm,nm,idat)

        sno=0.

        do n=1,nm
        do i=1,im
        do j=1,jm
          
!%%%	  do i=0,im+1
!          do j=0,jm+1
	  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
            
            
                  if(hphip.gt.(-60).and.hphip.lt.0)then !!!!!dragan,june 2012
                    
                    sm(i,j,n)=1.
                    sno(i,j,n)=sfcgridhi(i,j,n,13)/200*0.2
                    si(i,j,n)=sno(i,j,n)/0.2
                 !   si(i,j,n)=0.
                    sice(i,j,n)=0.
            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,idat)

       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(2,FILE='SICE_sfcEta2.txt')
       WRITE(2,*) SICE
       CLOSE(2)
!
       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)
!!      write(0,*),'ivgtyp=',ivgtyp
       end subroutine sfceta2