#include "options.inc"

c
c--------------------------------------------------------
c     preparation of the topography
c--------------------------------------------------------
c

#ifdef partial_cell
c I have used this one for E5 setup
c#define auto_find_digged_points

c this one for the DS setup
c#define preserve_mask_in_regions
#endif

      subroutine prep_topo
      use prep_module
      use island_module
      implicit none
      integer nx,ny
      real, allocatable :: rose(:,:),x(:),y(:),nrose(:,:)
      real, allocatable :: tmprose(:,:)
      real, allocatable :: tx(:)
      real, allocatable :: fmsk(:,:),srose(:,:)
      integer, allocatable :: map(:,:)
      integer i,j ,k
      integer :: indp
#ifdef partial_cell
      integer :: kmt_bak(imt,jmt),is,ie,js,je
      real  :: rose_digged(imt,jmt)
#endif
      real :: p_cell_min

      integer, parameter :: maxipp=40000,mnisle=500
      integer iperm(maxipp),jperm(maxipp)
      integer nippts(mnisle), iofs(mnisle)
      integer nisle
      logical error
      character(len=80) :: FLAME_ID ='no flame id '
      character(len=32+80+60) :: iotext='no stamp, no iotext '

      integer :: lenw,init
      real, allocatable :: work(:)

      print*,''
      print*,' Setting up the topography mask'
      print*,''

      allocate(map(imt,jmt) )
c
c     call templates to read topography file
c
      call read_rose_dims(nx,ny)
      print*,' data set dimensions are ',nx,ny
      allocate( rose( nx,ny), x(nx), y(ny), tx(nx) )
      print*,' now reading data '
      call read_rose(rose,x,y,nx,ny)
c
      print*,' interpolating data'
      allocate( nrose(imt,jmt), fmsk(imt,jmt), srose(imt,jmt) )
      if (.not.rotated_grid) then
c       use cheaper routine here
        lenw=10*max(nx,ny); allocate(work(lenw))
        init=1
        if (cyclic) then
         call tranlon (rose, nx, nx-2, ny, x, xt(2:imt-1),imt-2, tx)
         call ctf (rose,nx, ny, tx, y, 
     &                nrose,imt, jmt, 2, imt-1, 1, jmt,
     &                xt, yt, init, work, lenw,spval)
        else
         call tranlon (rose, nx, nx-2, ny, x, xt,imt, tx)
         call ctf (rose,nx, ny, tx, y, 
     &                nrose,imt, jmt, 1, imt, 1, jmt,
     &                xt, yt, init, work, lenw,spval)
        endif
        deallocate(work)
      else
c      there is one point at i=1 which is the same as i=nx
c      this is no good for interpolator
c       call rot_intrp_sclr (rose(1:nx-1,:), x, y, 
       allocate(tmprose(nx-1,ny))
       do i=1,nx-1;do j=1,ny
        tmprose(i,j)=rose(i,j)
       enddo;enddo
       print*,'ok'
       call rot_intrp_sclr (tmprose, x, y, 
     &              nx-1, ny, nrose, xt, yt,
     &              imt, jmt, psir, thetar, phir,spval,.true.)
       print*,'ok'
       deallocate(tmprose)
      endif
      print*,' interpolation done'
      deallocate( rose,x,y, tx)

      if (cyclic) then
        nrose(1,:)=nrose(imt-1,:); nrose(imt,:)=nrose(2,:)
      else
        nrose(1,:)=nrose(2,:); nrose(imt,:)=nrose(imt-1,:)
      endif
c
c-----------------------------------------------------------------------
c     convert depths from meters to cm and change sign. set min 
c     allowable depth
c     also set a mask 'fmsk' to be used in 'sfltr'.
c     -> since etopo5 changed its sign. the following procedure was
c        altered slightly; jk, may-05-1997
c-----------------------------------------------------------------------
c
      do j=1,jmt
        do i=1,imt
          fmsk(i,j)=1.0
          map(i,j)=0
          if (nrose(i,j) < 0. ) then
            nrose(i,j) = -nrose(i,j)*100.0
          else
            nrose(i,j) = 0.
            fmsk(i,j)=0.
            map(i,j)=1.
          endif
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     set topographic modifications on model grid here
c-----------------------------------------------------------------------
c
c     set artificial walls at northern & southern boundaries
c
      nrose(:,1 )  = 0; fmsk(:,1 )  = 0
      nrose(:,jmt) = 0; fmsk(:,jmt) = 0

      if (.not. cyclic) then
        nrose(1,: )  = 0; fmsk(1,: )  = 0
        nrose(imt,:) = 0; fmsk(imt,:) = 0
      endif
c
c     smooth, but preserve coastline
c
      srose=nrose
      if (topo_smooth_ntimes .ne. 0) then

        print*,' smoothing ',topo_smooth_ntimes,' times'

        call sfltr (srose, imt, imt, jmt, fmsk, 
c     &         2, imt-1, 2, jmt-1, topo_smooth_ntimes,cyclic)
     &         1, imt, 1, jmt, topo_smooth_ntimes,cyclic)
      endif

      allocate( kmt(imt,jmt) )

      if (read_grid_topo) then
c
c-----------------------------------------------------------------------
c       now read a previosly generated kmt mask 
c       and forget about the data
c-----------------------------------------------------------------------
c
        call read_bin_kmt()

      else
c
c-----------------------------------------------------------------------
c     discretize topography to  "kmt" (min = 2 levels)
c-----------------------------------------------------------------------
c
        do j=1,jmt
         do i=2,imt-1
          kmt(i,j) = 0
          if (srose(i,j) /= 0) then
           kmt(i,j) = max(kmt_min,indp (srose(i,j), zw, km))
c            kmt(i,j) = indp (srose(i,j), zw, km)
          endif
         enddo
        enddo
c
c-----------------------------------------------------------------------
c      apply boundary conditions
c-----------------------------------------------------------------------
c
       call set_kmt_bc()
c
c-----------------------------------------------------------------------
c     limit the minimum number of levels. kmt_min should be >= 2
c-----------------------------------------------------------------------
c
       if (kmt_min <2) then
        print*,' kmt_min must be greater than 2'
        stop
       endif
 
       do i=2,imt-1
        do j=jmt-1,2,-1
          if (kmt(i,j) /= 0 .and. kmt(i,j) < kmt_min) then
            kmt(i,j)=max(kmt(i,j), kmt_min, 2)
          endif
        enddo
       enddo

      endif ! read_grid_topo

#ifdef partial_cell
c
c-----------------------------------------------------------------------
c     construct partial cell from undiscretised kmt mask
c     htp = depth to ocean bottom at T-cells (cm)
c-----------------------------------------------------------------------
c
      allocate( htp(imt,jmt), hup(imt,jmt) )
      htp=0. ; hup=0.
      rose_digged=spval
      kmt_bak=kmt

      do j=1,jmt
       do i=1,imt
         if (kmt(i,j)==0) then
          htp(i,j)=0.  ! preserve coastline
         else
          if (srose(i,j) > 0.) then
            if (zw(kmt(i,j)) > srose(i,j) ) then
               ! fill with a partial cell
               htp(i,j)=max(zw(kmt_min),srose(i,j))
               if (zw( max(1,kmt(i,j)-1) ) > htp(i,j)) then
                  ! seems to be an artifically digged point
#ifdef auto_find_digged_points
                  if (zw( max(1,kmt(i,j)-1) ) 
     &                -dzt( max(1,kmt(i,j)))/5.> htp(i,j)) then
#endif
                     ! artifically digged point, preserve depth
                     htp(i,j)=zw(kmt(i,j))
c                     print*,' found an articially deep point at ',
c     &                     ' i=',i,' j=',j,' ,... preserving '
                     rose_digged(i,j)=srose(i,j)/100.
#ifdef auto_find_digged_points
                  else
c                     ! artifically digged point, remove that stuff
 30                   kmt(i,j) = kmt(i,j)-1
                      if (kmt(i,j)<kmt_min) then
                        ! no partial cell possible, too shallow
                        kmt(i,j)=kmt_min
                        htp(i,j)=zw(kmt(i,j))
                        rose_digged(i,j)=srose(i,j)/100.
                        goto 35
                      endif
                      if (zw( max(1,kmt(i,j)-1) ) > htp(i,j)) goto 30
 35                  continue

                  endif
#endif
               endif
            elseif (zw(kmt(i,j)) < srose(i,j) ) then
 20               if (kmt(i,j)==km) then
                    htp(i,j)=zw(kmt(i,j)) ! max depth reached
c                    print*,' max. depth reached for i=',i,' j=',j
                  else  
                   ! increase kmt(i,j) and fill with a partial cell beneath
                   kmt(i,j)=kmt(i,j)+1
                   htp(i,j)=srose(i,j)
                   if (zw(kmt(i,j)) < htp(i,j)) then
c                      print*,' digging',i,j
                      goto 20 ! still to shallow, go to the next deeper box
                   else
c                    print*,' depth increased at i=',i,' j=',j
                   endif
                   
                  endif
            else
              htp(i,j)=zw(kmt(i,j)) ! the best case
            endif
          else
           htp(i,j)=zw(kmt(i,j))  ! we have here an artificial water point
           rose_digged(i,j)=srose(i,j)/100.
          endif
         endif ! kmt==0
       enddo
      enddo

#ifdef preserve_mask_in_regions
c
c     now preserve to water depth from the old mask in some sensitive regions
c

c 1. Denmark Strait

      js = indp (65.6, yt, jmt);je = indp (66.7, yt, jmt)
      is = indp (360-28.0,xt,imt);ie=indp (360-26., xt, imt)
      do j=js,je
       do i=is,ie
         if (htp(i,j) > 0.) then
           if (zw(kmt(i,j)) > htp(i,j) ) then
              htp(i,j)=zw(kmt(i,j))
           endif
         endif
       enddo
      enddo

c 2. Faroer channel
c    okay

c 3. Florida Strait
c     okay

c
c    I do not care about the rest ...
c

#endif
c
c     a check
c
      do j=1,jmt
       do i=1,imt
        if (htp(i,j)>zw(max(1,kmt(i,j)))) then
         print*,' error at i=',i,' j=',j
         print*,' htp(i,j)=',htp(i,j)
         print*,' zw(kmt(i,j))=',zw(max(1, kmt(i,j) ))
         print*,' srose(i,j)=',srose(i,j)
         stop
        endif
        if (htp(i,j)<zw(max(1, kmt(i,j)-1 )).and.kmt(i,j)/=0) then
         print*,' error at i=',i,' j=',j
         print*,' htp(i,j)=',htp(i,j)
         print*,' zw(kmt(i,j)-1)=',zw(max(1, kmt(i,j)-1 ))
         print*,' srose(i,j)=',srose(i,j)
         print*,' kmt(i,j)=',kmt(i,j)
         stop
        endif
       enddo
      enddo
c
c---------------------------------------------------------------------
c     do not allow partial cells thickness less than min(dzt(k))
c     (arbitrarily 50m minimum in coarse models ... otherwise dzt(1))
c---------------------------------------------------------------------
c
      p_cell_min = min(50.0e2,zw(1))
      do k=2,km
        if (p_cell_min .gt. (zw(k)-zw(k-1))) p_cell_min = zw(k)-zw(k-1) 
      enddo
      do j=2,jmt-1
        do i=2,imt-1
	  k = kmt(i,j)
	  if (k .lt. 2) then
	    htp(i,j) = 0.0
	  else
	    if ((htp(i,j)-zw(k-1)) .lt. p_cell_min) then
c             print*,' WARNING: partial cell limit violated'
c             print*,' at i=',i,' j=',j,' htp(i,j)=',htp(i,j)
c             print*,' zw(k-1)=',zw(k-1)
c             print*,' zw(k)=',zw(k),' p_cell_min',p_cell_min
             htp(i,j) = zw(k-1) + p_cell_min
	    endif
	  endif
        enddo
      enddo 
#endif
c
c-----------------------------------------------------------------------
c     test for "island perimeter violations"
c     channels between different land masses must be at least 2 ocean
c     cells wide to prevent conflicting values of the stream function
c     on these island perimeter points.
c-----------------------------------------------------------------------
c
 10   continue
#ifdef partial_cell
      kmt_bak=kmt
#endif
       write (6,'(/a/)')
     &  'Searching for and correcting island PERIMETER VIOLATIONS...'
       call isleperim(kmt,map, iperm, jperm, iofs, 
     &        nippts, nisle
     &,       imt, jmt, mnisle, maxipp,0,cyclic,
     &        enable_obc_north,enable_obc_south,
     &        enable_obc_west,enable_obc_east,.true.)

       if (.not. read_grid_topo) then
        write (6,'(/a/)')
     &   'Searching for isolated bays, potholes, etc ...'
        call fill_bays_and_stuff(kmt,imt,jmt,kmt_min,error,
                            ! check for potholes, trenches/corners
     &                          .false., .false.)
       else
        print*,' searching for depth violations and isolated bays only'
        call fill_bays_and_stuff(kmt,imt,jmt,kmt_min,error,
                            ! leave potholes, trenches/corners
                            ! perform only necessary checks
     &                          .true., .true.)
       endif

#ifdef partial_cell
       do j=1,jmt
        do i=1,imt
          if (kmt_bak(i,j)/=kmt(i,j) ) then
           if (kmt(i,j)==0) then
            htp(i,j)=0.
           else
            htp(i,j)=zw(kmt(i,j))
           endif
          endif
        enddo
       enddo
#endif
      if (error) goto 10
c
c---------------------------------------------------------------------
c     topography for u grid,  it is needed afterwards
c---------------------------------------------------------------------
c
      allocate(kmu(imt,jmt))
      kmu=0
      do j=1,jmt-1
        do i=1,imt-1
           kmu(i,j) = min (kmt(i,j), kmt(i+1,j), 
     &                     kmt(i,j+1), kmt(i+1,j+1))
        enddo
      enddo
      if (enable_obc_north) kmu(:,jmt) = kmu(:,jmt-1)
      if (enable_obc_south) kmu(:,1) = kmu(:,2)
      if (enable_obc_west) kmu(1,:) = kmu(2,:)
      if (enable_obc_east) kmu(imt,:) = kmu(imt-1,:)
      if ((enable_obc_north .or. enable_obc_south .or.
     & enable_obc_west  .or. enable_obc_east) .and..not.cyclic ) then 
c      checks in case of o.b.c
       kmu(imt-1:imt,1)=0.
       kmu(1,1)=0.
       kmu(1,jmt-1:jmt)=0.
c      done
      endif
      if (cyclic) kmu(imt,:) = kmu(2,:)

#ifdef partial_cell

c
c     construct hup to show it in the netcdf file
c
      hup=0.
      do j=1,jmt-1
        do i=2,imt-1
          hup(i,j) = min( htp(i,j),htp(i+1,j),htp(i,j+1),htp(i+1,j+1))
	enddo
      enddo

      if (enable_obc_north) hup(:,jmt) = hup(:,jmt-1)
      if (enable_obc_south) hup(:,1)   = hup(:,2)
      if (enable_obc_west)  hup(1,:)   = hup(2,:)
      if (enable_obc_east)  hup(imt,:) = hup(imt-1,:)

      if ((enable_obc_north .or. enable_obc_south .or.
     &     enable_obc_west  .or. enable_obc_east).and..not.cyclic) then 
c      checks in case of o.b.c
       hup(imt-1:imt,1)=0.
       hup(1,1)=0.
       hup(1,jmt-1:jmt)=0.
c      done
      endif

      if (cyclic) hup(imt,:) = hup(2,:)
      if (cyclic) hup(1,:) = hup(imt-1,:)

#endif

c      call showmap(map,imt,jmt)
      deallocate(map) 
c
c---------------------------------------------------------------------
c     write the topography to netcdf file 
c---------------------------------------------------------------------
c
      call write_kmt_to_cdf(nrose,srose
#ifdef partial_cell
     & ,rose_digged
#endif
     &     )

      deallocate( nrose, fmsk, srose )
c
c---------------------------------------------------------------------
c     write the kmt to binary file new_kmt.dta
c---------------------------------------------------------------------
c
      call write_bin_kmt()

      print*,''
      print*,' Topography generation done'
      print*,''

      end subroutine

 
      subroutine sfltr (a, im, il, jl, f, i1, i2, j1, j2, 
     &                  nfil,cyclic)
c
c-----------------------------------------------------------------------
c     this 2 dimensional symmetric filter is based on a weighting of
c     1/4  1/2  1/4. 
c     f is the mask field = 0 over land & 1 in ocean
c     a is the data to be smooth. It must be zero over land !
c     nfil is the number of smoothing passes
c     i1 .. i2 are the points along the 1st dimension to smooth
c     j1 .. j2 are the points along the 2nd dimension to smooth
c     this filter behaves as if there is a no flux boundary condition
c     across coastlines. It therefore preserves coasts.
c     note: this filter does not take into account variable grids
c-----------------------------------------------------------------------
c           
      dimension a(im,jl), s(im,jl), f(im,jl)
      logical cyclic
c
      if (nfil .eq. 0) return
c
      do 20 j=1,jl
        do 10 i=1,il
          if (f(i,j) .eq. 0.0 .and. a(i,j) .ne. 0.0) then
            print *,' "a" is not zero but "f" is zero in "sfltr"'
            print *,' i=',i,' j=',j,' a(i,j)=',a(i,j)
            stop
          endif
10      continue
20    continue
c
      is = max (i1,2); ie = min (i2,il-1)
      js = max (j1,2); je = min (j2,jl-1)
      print *,' '
      print *,' ==> Applying smoothing filter ', nfil,' times over'
      print *,
     $'     domain defined by is=',is,', ie=',ie,', js=',js,', je=',je
      print *,' '
c
      do 100 n=1,nfil
        do 50 j=js,je
          do 40 i=2,il-1
            if (i .ge. is .and. i .le. ie) then
              x1p = 0.25*f(i,j)
              x2p = 1.0-0.25*(f(i-1,j+1)+f(i+1,j+1))
              x3p = 0.25*f(i,j)
              x1c = 0.25*f(i,j)
              x2c = 1.0-0.25*(f(i-1,j)+f(i+1,j))
              x3c= 0.25*f(i,j)
              x1m = 0.25*f(i,j)
              x2m = 1.0-0.25*(f(i-1,j-1)+f(i+1,j-1))
              x3m = 0.25*f(i,j)
              y1m = 0.25*f(i,j)
              y2m = 1.0-0.25*(f(i-1,j+1)+f(i-1,j-1))
              y3m = 0.25*f(i,j)
              y1c = 0.25*f(i,j)
              y2c = 1.0-0.25*(f(i,j+1)+f(i,j-1))
              y3c = 0.25*f(i,j)
              y1p = 0.25*f(i,j)
              y2p = 1.0-0.25*(f(i+1,j+1)+f(i+1,j-1))
              y3p = 0.25*f(i,j)
              s(i,j) = 
     $        x1p*y1m*a(i-1,j+1) + x2p*y1c*a(i,j+1) + x3p*y1p*a(i+1,j+1)
     $      + x1c*y2m*a(i-1,j)   + x2c*y2c*a(i,j)   + x3c*y2p*a(i+1,j)
     $      + x1m*y3m*a(i-1,j-1) + x2m*y3c*a(i,j-1) + x3m*y3p*a(i+1,j-1)
            else
              s(i,j) = a(i,j)
            endif
40        continue  
50      continue
        if (cyclic) then
         do  j=js,je
          s(1,j)  = s(il-1,j); s(il,j) = s(2,j)
         enddo
        else
         do j=js,je
          s(1,j)  = 0.0; s(il,j) = 0.0
         enddo
        endif
        do 80 j=js,je
          do 70 i=1,il
            a(i,j) = s(i,j)
70        continue  
80      continue
100   continue
      end subroutine


      subroutine fill_bays_and_stuff(kmt,imt,jmt,kmt_min,error,
     &   leave_potholes,leave_trenches)
      implicit none
      integer imt,jmt,kmt(imt,jmt),nkmt(imt,jmt),kmt_min
      integer kmu,i,j,k
      logical :: error , leave_potholes,leave_trenches

      kmu(i,j) = min(kmt(i,j), kmt(i+1,j), kmt(i,j+1), kmt(i+1,j+1))
      error=.false.
c
c-----------------------------------------------------------------------
c     limit the minimum number of levels. kmt_min should be >= 2
c-----------------------------------------------------------------------
c
      if (kmt_min <2) then
        print*,' kmt_min must be greater than 2'
        stop
      endif

      do i=2,imt-1
        do j=jmt-1,2,-1
          if (kmt(i,j) /= 0 .and. kmt(i,j) < kmt_min) then
            kmt(i,j)=max(kmt(i,j), kmt_min, 2)
            error=.true.
            print*,' minimal depth violation at i=',i,' j=',j
          endif
        enddo
      enddo
      call set_kmt_bc()
c
c-----------------------------------------------------------------------
c     detect isolated bays:
c     These are situations where an ocean point on the "t" grid is
c     surrounded by 4 land points on the "u" grid (the 4 verticies of
c     the "t" cell) which prevents lateral advection of tracers (T & S).
c     These points sometimes lead to numerical instabilities and can be
c     eliminated by changing the ocean "t" cell to a land value.
c-----------------------------------------------------------------------
c
      nkmt=kmt
      do i=2,imt-1
        do j=jmt-1,2,-1
          if (kmt(i,j) .ne. 0) then
            if (kmu(i  ,j) == 0 .and. kmu(i  ,j-1) == 0 .and.
     &          kmu(i-1,j) == 0 .and. kmu(i-1,j-1) == 0) then
c             option 1: fill cell to land
              nkmt(i,j)=0
c             option 2: interactive user changes
c             option 3: leave nonadvective bay (requires no kmt changes)
              error=.true.
              print*,' found isolated bay at i=',i,' j=',j
            end if
          endif
        enddo
      enddo
      kmt=nkmt
      call set_kmt_bc()


      if (.not. leave_potholes) then
c
c-----------------------------------------------------------------------
c     detect potholes:
c     These are situations where an ocean point on the "t" grid is
c     deeper than its 4 surrounding neighbor cells, which prevents
c     lateral and vertical advection of tracers at this depth.
c     These points usually cause no problems, but they also do not
c     contribute to ocean circulation.  They may be eliminated by
c     changing kmt at this "t" grid location to the maximum of kmt on
c     the 4 neighboring locations.
c-----------------------------------------------------------------------

      nkmt=kmt
      do i=2,imt-1
        do j=jmt-1,2,-1
          if (kmt(i+1,j) .lt. kmt(i,j) .and.
     &        kmt(i-1,j) .lt. kmt(i,j) .and.
     &        kmt(i,j+1) .lt. kmt(i,j) .and.
     &        kmt(i,j-1) .lt. kmt(i,j)) then
c
c           option 1: fill this cell to the max of its neighbors
c           note: this still is not enough to guarantee flow
c
            nkmt(i,j) = max(kmt(i+1,j), kmt(i-1,j),
     &                      kmt(i,j+1), kmt(i,j-1))
c
c           option 2: interactive user changes
c           option 3: leave pothole (requires no kmt changes)
c
              print*,' found pothole at i=',i,' j=',j
              error=.true.
          endif
        enddo
      enddo
      kmt=nkmt
      call set_kmt_bc()
 
      endif


      if (.not. leave_trenches) then
c
c     Search for isolated corners and trenches, i.e. "t" grid boxes 
c     at any depth that can`t be influenced by advection
c     -------------------------------------------------------------
c
c      Examples: (x=land, .=sea)
c
c      1) x...x   2) .....   3) xxx..   4) x....
c         xx.xx      .....      xx..x      xx...
c         xx.xx      xx.xx      x..xx      x..xx
c         x...x      xxxxx      ..xxx      x....
c
c      Note: Case 4 is not covered by the present routine! However, it
c            is probably arguable if a diffusive flux should be 
c            possible through a gap like that.
c
c      This procedure mimics a check that will later be repeated for
c      the surface level in MOM`s iperim.F.

c     Compute number of vertical levels on the "u" grid
c     -------------------------------------------------

      nkmt=kmt
      do j=2, jmt-1
        do i=2, imt-1
          k = kmt(i,j)
          if (k .ne. 0) then
            if (kmu(i,j) .lt. k .and. kmu(i-1,j) .lt. k .and.
     &        kmu(i,j-1) .lt. k .and. kmu(i-1,j-1) .lt. k) then
              k = max (kmu(i,j),   kmu(i-1,j), 
     &                 kmu(i,j-1), kmu(i-1,j-1))
c             guarantee min. number of levels:
              if (k.ne.0 .and. k.lt.kmt_min) k = kmt_min 
              nkmt(i,j)=k
              error=.true.
              print*,' found isolated trench/corner at i=',i,' j=',j
            endif
          endif
        end do
      end do
      kmt=nkmt
      call set_kmt_bc()

      endif

      end subroutine 


      subroutine set_kmt_bc()
      use prep_module
c-----------------------------------------------------------------------
c      apply boundary conditions
c-----------------------------------------------------------------------
       if (cyclic) then
         kmt(1,:)   = kmt(imt-1,:); kmt(imt,:) = kmt(2,:)
       else
         kmt(imt,:) = 0; kmt(1,:) = 0
         if (enable_obc_east) then
          kmt(imt,:) = kmt(imt-1,:)
          i=imt-1
          do j=1,jmt
           k=min(min(kmt(i,j),kmt(i-1,j)),min(kmt(i,j),kmt(i-2,j)))
           if (k /= kmt(i,j)) then
            print*,' Changing Topography at eastern OB. '
            print*,'  at i=',i,' j=',j
            print*,'  new kmt(i,j)=',k, 'old kmt(i,j)=',kmt(i,j)
            kmt(i-2:i+1,j)=k
           endif
          enddo
         endif
         if (enable_obc_west) then
          kmt(1,:) = kmt(2,:)
          i=2
          do j=1,jmt
           k=min(min(kmt(i,j),kmt(i+1,j)),min(kmt(i+2,j),kmt(i,j)))
           if (k /= kmt(i,j)) then
            print*,' Changing Topography at western OB. '
            print*,'  at i=',i,' j=',j
            print*,'  new kmt(i,j)=',k, 'old kmt(i,j)=',kmt(i,j)
            kmt(i-1:i+2,j)=k
           endif
          enddo
         endif
       endif
       kmt(:,jmt) = 0 ; kmt(:,1) = 0
       if (enable_obc_north) then
         kmt(:,jmt) = kmt(:,jmt-1)
         j=jmt-1
         do i=1,imt
          k=min(min(kmt(i,j),kmt(i,j-1)),min(kmt(i,j),kmt(i,j-2)))
          if (k /= kmt(i,j)) then
           print*,' Changing Topography at northern OB. '
           print*,'  at i=',i,' j=',j
           print*,'  new kmt(i,j)=',k,' old kmt(i,j)=',kmt(i,j)
           kmt(i,j-2:j+1)=k
          endif
         enddo
       endif
       if (enable_obc_south) then
         kmt(:,1) = kmt(:,2)
         j=2
         do i=1,imt
          k=min(min(kmt(i,j),kmt(i,j+1)),min(kmt(i,j),kmt(i,j+2)))
          if (k /= kmt(i,j)) then
           print*,' Changing Topography at southern OB. '
           print*,'  at i=',i,' j=',j
           print*,'  new kmt(i,j)=',k,' old kmt(i,j)=',kmt(i,j)
           kmt(i,j-1:j+2)=k
          endif
         enddo
       endif
      end subroutine




      subroutine read_bin_kmt()
      use prep_module
      implicit none
      integer  :: i,j,k,io
      character(len=80) :: FLAME_ID ='no flame id '
      character(len=32+80+60) :: iotext='no stamp, no iotext '

      print*,' reading topography mask from file ',
     &    old_kmt_file(1:len_trim(old_kmt_file))

      call getunit(io,old_kmt_file,'usr ieee')
      read (io) flame_id  ! no checks
c       skip header record and read dimensional info
      read (io) 
      read (io) i, j, k
c       check that grid sizes in argument list 
      if (i/=imt .or. j /= jmt .or. k/= km) then
c        check i,j,k
         print*,' Error: '
         print*,' read imt,jmt,km=',i,j,k
         print*,' previously defined as ',imt,jmt,km
         call halt_stop(' reading kmt file ')
      endif
      read (io) 
      read (io) kmt
      close (io)
      end subroutine



      subroutine write_bin_kmt()
      use prep_module
      implicit none
      integer  :: io
      character(len=80) :: FLAME_ID ='no flame id '
      character(len=32+80+60) :: iotext='no stamp, no iotext '

      print*,' writing topography mask to file ',
     &    new_kmt_file(1:len_trim(new_kmt_file))

      call getunit(io,new_kmt_file,'usr ieee')
      write (io) flame_id ! no checks
      write (io) iotext
      write (io) imt, jmt, km
      write (io) iotext
      write (io) kmt
      close (io)
#ifdef partial_cell
      print*,' writing partial cell mask to file ',
     &    new_htp_file(1:len_trim(new_htp_file))

      call getunit(io,new_htp_file,'usr ieee')
      write (io) flame_id ! no checks
      write (io)  iotext
      write (io) imt, jmt, km
      write (io)  iotext
      write (io)  htp
      close (io)
#endif
      end subroutine
