subroutine sfcEta(eta,aeta,dfl,im,jm,lm,nm,flat)

      implicit none

      integer, intent(in) :: im,jm,lm,nm
      logical, intent(in) :: flat
      real,dimension(lm+1), intent(in) ::eta,dfl
      real,dimension(lm), intent(in) :: aeta

      include 'const.h'

      integer::i,j,n,l,lmp1,lmnnw,lmnne,lmnsw,lmnse   &
              ,lnbmx,lrais,mmx,m,lij ,k
      real,dimension(0:im+1,0:jm+1,nm)::hgt,sm,fis,etas,wm &
                                       ,res,ref
      integer,dimension(0:im+1,0:jm+1,nm)::lhgt
      real,dimension(0:im+1,0:jm+1,nm,lm)::htm,vtm
      integer,dimension(0:im+1,0:jm+1,nm)::lmh,lmv
      integer,dimension(0:im+1,0:jm+1,nm)::lmin
      real,dimension(0:im+1,0:jm+1,nm)::dum
      real::prfs
      real,parameter::prf0=101325., t0=288.,g=9.80616,r=287.04,  &
             gamma=0.0065,gorg = g/(r*gamma)
      real,dimension(4)::seta
      real::hgtmax
      integer::imaxh,jmaxh,nmaxh
      integer::isCorner2,isCorner2v

      hgtmax=0.

      if(.not.flat)then

      lmp1=lm+1
      open(1,file='topo.dat', status='old',form='unformatted')
      read(1) hgt,sm
!      print*,hgt
      close(1)

      do l=1,lm+1
        print *,l,dfl(l)
      enddo

      do n=1,nm
      do j=1,jm    
      do i=1,im    


         prfs=prf0*((t0-gamma*hgt(i,j,n))/t0)**gorg
         etas(i,j,n)=(prfs-pt)/(prf0-pt)

!mp        through this point hgt represents the unchanged raw hgt value
!
         l=lm
         do while(etas(i,j,n).lt.eta(l).and.l.gt.0) 
     
!!!!!DRAGAN test      do while(etas(i,j,n).lt.eta(l))
           l=l-1
!!!!! DRAGAN test       if(l.eq.1) go to 450
	
	 enddo
!!!!!DRAGAN test    450       continue
           if (etas(i,j,n).lt.aeta(l)) then
               ref (i,j,n)= eta(l)
               hgt (i,j,n)=dfl(l)
               lhgt(i,j,n)=     l
            else
               ref (i,j,n)=eta(l+1)
               if (l+1 .eq. lmp1 .and. hgt(i,j,n) .gt. 0.) then
                  hgt(i,j,n)=0.004
               else
                  hgt(i,j,n)=dfl(l+1)
               endif
               lhgt(i,j,n)=    l+1
            endif
            if(hgt(i,j,n).gt.hgtmax)then
              hgtmax=hgt(i,j,n)
              imaxh=i
              jmaxh=j
              nmaxh=n
            endif
      enddo
      enddo
      enddo

      print *,"max hgt ...",imaxh,jmaxh,nmaxh,hgtmax

      call bocoh(etas,im,jm,nm)
      call bocoh(ref,im,jm,nm)
      call bocoh(hgt,im,jm,nm)
      call bocoh(lhgt,im,jm,nm)

!      find points having no neighboring wind point because of mountains
!c     on all four sides (sunken points), and
!c     *  if land, just raise to reach the lowest wind;
!c     *  if an isolated water point, or a water point with only one
!c          neighboring water point and at sea level: raise, and declare
!c          land (in order not to trigger convection over artificially
!c          elevated warm water);
!c     *  otherwise:
!c          remove (level off) land at one of the four neighboring wind
!c          points -- the one having the smallest three-point averaged
!c          elevation -- so as to restore one of the wind points
!c
!c

!c
!c     Save first the sea mask, sm, into a back-up mask ("water mask",
!c     wm) to use in testing the neighbors, since sm will be changed
!c     along the way; and initialize lmin, level beyond which points
!c     should not be raised because they have participated in leveling
!c     off at that level
!c
       do n=1,nm
       do j=0,jm+1
       do i=0,im+1
          wm  (i,j,n)=sm  (i,j,n)
          lmin(i,j,n)=2
       enddo
       enddo
       enddo


       do n=1,nm
       do j=1,jm
       do i=1,im

        
!
!c ************ Find the maximum elevation (minimum l) of four neighboring wind
!c              points -- to the nw, ne, sw and se
!c

        lmnnw=min(lhgt(i-1,j,n),lhgt(i-1,j+1,n),lhgt(i,j+1,n))
        lmnne=min(lhgt(i+1,j,n),lhgt(i+1,j+1,n),lhgt(i,j+1,n))
        lmnsw=min(lhgt(i-1,j,n),lhgt(i-1,j-1,n),lhgt(i,j-1,n))
        lmnse=min(lhgt(i,j-1,n),lhgt(i+1,j-1,n),lhgt(i+1,j,n))
!c
!c ************ Which one is the lowest of the four maximum elevations?
!c
        lnbmx=max(lmnnw,lmnne,lmnsw,lmnse)
!c
!c ************ Is the point's original elevation below that of the lowest
!c              available original wind, and, at the same time, in case it
!c              participated in leveling off is it still below the level of the
!c              lowest wind created by leveling off?
!c

           if (lhgt(i,j,n) .gt. lnbmx .and.    &
             ref(i,j,n) .gt. eta(lmin(i,j,n))) then
                  lrais=max0(lnbmx,lmin(i,j,n))
                  if (sm(i,j,n) .eq. 0.) then
                     ref (i,j,n)= eta(lrais)
                     hgt (i,j,n)=dfl(lrais)

                  else if(wm(i,j-1,n)+wm(i,j+1,n)        &
                    +wm(i-1,j,n)+wm(i+1,j,n).eq.0..or.   &
                     wm(i,j-1,n)+wm(i,j+1,n)             &
                    +wm(i-1,j,n)+wm(i+1,j,n).eq.1..and.  &
                       lhgt(i,j,n).eq.lm+1) then
                     ref (i,j,n)= eta(lrais)
                     hgt (i,j,n)=dfl(lrais)

                  else
!c *************** There is a neighboring water point, and if it is only a
!c                 single neighboring water point the water points are above sea
!c                 level.  denote wind points to the east, north, west
!c                 and south by 1, 2, 3, and 4, respectively, and check the
!c                 three-point volume at each of them.
!c

                   seta(1)=etas(i-1,j,n)+etas(i-1,j+1,n)+etas(i,j+1,n)
                   seta(2)=etas(i,j+1,n)+etas(i+1,j+1,n)+etas(i+1,j,n)
                   seta(3)=etas(i-1,j,n)+etas(i-1,j-1,n)+etas(i,j-1,n)
                   seta(4)=etas(i,j-1,n)+etas(i+1,j-1,n)+etas(i+1,j,n)


                     mmx=1
                     do m=2,4
                        if(seta(m) .gt. seta(mmx)) mmx=m
                     enddo
!c
!c ****************** n.b.: When removing now the obstructing land around the 
!c                    wind point vm, lhgt is not changed.  Thus, hopefully (?) 
!c                    resulting topography does not depend on the direction of 
!c                    the sweep within the current loop.
!c
                     lij=lhgt(i,j,n)

                     if (mmx .eq. 1 .or. mmx .eq. 3) then
                             if(ref(i-1,j,n).lt.eta(lij)) then
                                     ref(i-1,j,n)=eta(lij)
                                     hgt(i-1,j,n)=dfl(lij)

                             endif
                         lmin(i-1,j,n)=max(lmin(i-1,j,n),lij)
                     endif

                     if (mmx .eq. 1) then
                        if (ref(i-1,j-1,n) .lt. eta(lij)) then
                           ref(i-1,j-1,n)=   eta(lij)
                           hgt(i-1,j-1,n)=  dfl(lij)

                        endif
                        lmin(i-1,j-1,n)=max(lmin(i-1,j-1,n),lij)
                     endif

                      if (mmx .eq. 1 .or. mmx .eq. 2) then
                             if(ref(i,j+1,n).lt.eta(lij)) then
                                      ref(i,j+1,n)=eta(lij)
                                     hgt(i,j+1,n)=dfl(lij)

                      endif
                        lmin(i,j+1,n)=max(lmin(i,j+1,n),lij)
                      endif

                     if (mmx .eq. 2) then
                        if (ref(i+1,j+1,n) .lt. eta(lij)) then
                           ref(i+1  ,j+1 ,n )=   eta(lij)
                           hgt(i+1  ,j+1  ,n)=  dfl(lij)

                        endif
                        lmin(i+1,j+1,n)=max(lmin(i+1,j+1,n),lij)
                     endif

                     if (mmx .eq. 2 .or. mmx .eq. 4) then
                         if(ref(i+1,j,n).lt.eta(lij)) then
                               ref(i+1,j,n)=eta(lij)
                                hgt(i+1,j,n)=dfl(lij)

                          endif
                          lmin(i+1,j,n)=max(lmin(i+1,j,n),lij)
                     endif

                     if (mmx .eq. 4) then
                        if (ref(i+1,j-1,n) .lt. eta(lij)) then
                           ref(i+1,j-1,n)=   eta(lij)
                           hgt(i+1,j-1,n)=  dfl(lij)

                        endif
                        lmin(i+1,j-1,n)=max(lmin(i+1,j-1,n),lij)
                     endif

                     if (mmx .eq. 3 .or. mmx .eq. 4) then
                      if(ref(i,j-1,n).lt.eta(lij)) then
                       ref(i,j-1,n)=eta(lij)
                       hgt(i,j-1,n)=dfl(lij)

                        endif
                       lmin(i,j-1,n)=max(lmin(i,j-1,n),lij)
                     endif

                     if (mmx .eq. 3) then
                       if(ref(i-1,j-1,n).lt.eta(lij)) then
                        ref(i-1,j-1,n)=eta(lij)
                        hgt(i-1,j-1,n)=dfl(lij)

                       endif
                        lmin(i-1,j-1,n)=max(lmin(i-1,j-1,n),lij)
 
                      endif
                  endif
               endif

         enddo
         enddo
         enddo
   
         call bocoh(ref,im,jm,nm)
         call bocoh(hgt,im,jm,nm)
 
         do n=1,nm
         do i=0,im+1
         do j=0,jm+1
           if(ref(i,j,n).ne.0)then
             res(i,j,n)=1./ref(i,j,n)
           endif
         enddo
         enddo
         enddo

         htm=1.                                                       
         vtm=1.                                                       
     
         do l=1,lm
         do n=1,nm                              
         do j=0,jm+1                                          
         do i=0,im+1
            if (eta(l+1) .gt. ref(i,j,n))then
               htm(i,j,n,l)=0.
            endif
         enddo
         enddo
         enddo
         enddo

        do n=1,nm
        do j=1,jm
        do i=1,im
        do l=1,lm
           if(eta(l+1).gt.ref(i,j,n)) then
               vtm(i-1,j-1,n,l)=0.
               vtm(i-1,j,n,l)=0.
               vtm(i,j-1,n,l)=0.
               vtm(i,j,n,l)=0.
            endif
        enddo
        enddo
        enddo
        enddo

        do n=1,nm
        do i=0,im
        do j=0,jm
          if(isCorner2v(i,j,n,im,jm,nm).ne.0)then
            vtm(i,j,n,:)=0.
          endif
        enddo
        enddo
        enddo


        lmh=lm
        lmv=lm

        do n=1,nm
        DO J=0,JM+1
        DO I=0,IM+1
        DO L=1,LM
          IF(LMH(I,J,n).EQ.LM.AND.ETA(L+1).GE.REf(I,J,n)) LMH(I,J,n)=L
        ENDDO
          if(isCorner2(i,j,n,im,jm,nm).ne.0)then
            lmh(i,j,n)=0
          endif
        ENDDO
        ENDDO
        ENDDO

        DO L=1,LM
        do n=1,nm
        DO J=1,JM-1
        DO I=1,IM-1
           IF(LMV(I,J,n).EQ.LM.AND.VTM(I,J,n,L).LT.0.1) LMV(I,J,n)=L-1
        ENDDO
        ENDDO
        ENDDO
        ENDDO

        else           !flat surface
          fis=0.
          res=1.
          htm=1.
          vtm=1.
          lmh=lm
          lmv=lm
        endif


      fis=hgt*g

      if(sigma) then
         lmh=lm
         lmv=lm
         vtm=1.

         do n=1,nm
         do i=0,im+1
         do j=0,jm+1
           if(isCorner2(i,j,n,im,jm,nm).ne.0)then
             htm(i,j,n,:)=0
             lmh(i,j,n)=0
             res(i,j,n)=0.
           else
             htm(i,j,n,:)=1
             lmh(i,j,n)=lm
             res(i,j,n)=1.
           endif
         enddo
         enddo
         enddo
      endif
           

      do l=1,lm
        call bocov1(vtm(:,:,:,l),im,jm,nm)
      enddo
       dum=lhgt
      call bocoh(dum,im,jm,nm)
      lhgt=dum
     
      open(1,file='sfceta.dat',form='unformatted')
      write(1)fis,res,htm,vtm,lmh,lmv        
      close(1)

      open(1,file='forVt.dat',form='unformatted')
      write(1)lhgt,hgt,ref
      close(1)

      open(1,file='sm.dat',form='unformatted')
      write(1)sm
      close(1)

      END SUBROUTINE sfcEta