! source file: /raid23/csomes/UVic/2.9/c_n_isotopes/last_glacial_experiments/dlgp15/updates/topog.F
      subroutine topog (kmt, kmu, map, xt, yt, zt, xu, yu, zw, imt2
     &,                 jmt2, km2, sg_bathy, sg_ocean_mask)

!-----------------------------------------------------------------------
!     construct the topographic mask "kmt" which determines the geometry
!     and topography by defining the number of vertical levels on model
!     "t" points.

!-----------------------------------------------------------------------
!     input:

!     xt  = longitudes of model grid "t" points (degrees)
!     yt  = latitudes of model grid "t" points (degrees)
!     zt  = depths of box "u/v" points on model grid (cm)
!     xu  = longitudes of model grid "u" points (degrees)
!     yu  = latitudes of model grid "u" points (degrees)
!     zw  = depths of box bottoms on model grid "t" boxes (cm)
!     imt = number of model longitudes
!     jmt = number of model latitudes
!     km  = maximum number of levels

!     output:

!     kmt = number of vertical levels on "t" points
!           0 is a land point.
!           n implies n boxes between the surface and bottom
!             where (n >= kmt_min) & (n <= km)
!     map = array indicating land masses and their perimeters
!           0 is a mid ocean point
!           n > 0 is a land mass
!           n < 0 is the ocean perimeter of land mass |n|

!-----------------------------------------------------------------------

      implicit none

      character(40) :: cktext
      character(120) :: iotext
      character(60) :: expnam
      character(32) :: stamp
      character(120) :: fname, new_file_name

!     kmt_min is the minimum permitted number of levels at an ocean
!     "T" cell.  kmt_min must be at least 2.
      integer kmt_min, imt2, jmt2, km2, ncase, iou, j, i, n
      integer nerror, jrow, npass, k, nisle, linewidth
      parameter (kmt_min = 2)

      real c0, c1, cksum1, checksumi, cksum2

      logical exists

      include "stdunits.h"
      include "size.h"
      include "isleperim.h"

      integer map(imt,jmt), iperm(maxipp), jperm(maxipp),nippts(mnisle)
      integer iofs(mnisle), kmt(imt,jmt), kmu(imt,jmt), ib(10), ic(10)

      real tmpij(imt-2,jmt-2), tmpijk(imt-2,jmt-2,km)
      real xt(imt), yt(jmt), zt(km), xu(imt), yu(jmt), zw(km)
      real sg_bathy(imt,jmt,km), sg_ocean_mask(imt,jmt)

!     the variable auto_kmt_changes tells whether any changes have been
!     made to the kmt field as a result of define options.
!     nchanges counts how many changes in each group.

      auto_kmt_changes = .false.
      nchanges = 0

!-----------------------------------------------------------------------
!     check that grid sizes in argument list and file "size.h" agree
!-----------------------------------------------------------------------

      call size_check (imt2, jmt2, km2, 'topog', 'stop')

      ncase = 0

      write (stdout,'(//,t37,a,/)')
     & 'T O P O G R A P H Y    G E N E R A T I O N'
      write (stdout,'(a,i4,/)') 'kmt_min = ', kmt_min

      ncase = ncase + 1
      c0 = 0.
      c1 = 1.
      fname = new_file_name ("G_kmt.nc")
      inquire (file=trim(fname), exist=exists)
      if (.not. exists) then
        print*, "Error => ", trim(fname), " does not exist."
        stop 'topog in topog.f'
      endif
      ib(:) = 1
      ic(:) = imt-2
      ic(2) = jmt-2
      call openfile (fname, iou)
      if (.not. exists) then
        print*, "==>  Warning: G_kmt data does not exist."
      else
        call getvara ('G_kmt', iou, ic(1)*ic(2), ib, ic, tmpij, c1, c0)
      endif
      kmt(2:imt-1,2:jmt-1) = tmpij(1:imt-2,1:jmt-2)

      nerror = 0
      do i=1,imt
        do j=1,jmt
          if (kmt(i,j) .lt. 0 .or. kmt(i,j) .gt. km) then
            nerror = nerror + 1
          endif
        enddo
      enddo
      if (nerror .gt. 0) then
        write (stdout,'(a,i4,a,a,i4,/)')
     &    '==>ERROR: There are ', nerror, ' depth values in the file '
     &,   '"kmt.dta" that are outside the range 0 <= kmt(i,j) <= km = '
     &,   km
        stop '==>topog'
      endif
      write (stdout,'(/a/)') '==> The kmt field has been imported.'

      fname = new_file_name ("G_subgrid_bathy.nc")
      inquire (file=trim(fname), exist=exists)
      if (.not. exists) then
         print*, "Error => ", trim(fname), " does not exist."
         stop 'topog in topog.f'
      endif
      ib(:) = 1
      ic(:) = imt-2
      ic(2) = jmt-2
      ic(3) = km
      call openfile (fname, iou)
      if (.not. exists) then
         print*, "==> Warning: G_subgrid_bathy data does not exist."
      else
         call getvara ('SUBGRID_BATHY', iou, ic(1)*ic(2)*ic(3), ib, ic
     &,             tmpijk, c1, c0)
      endif
      sg_bathy(2:imt-1,2:jmt-1,:) = tmpijk(1:imt-2,1:jmt-2,:)

      fname = new_file_name ("G_subgrid_ocean_mask.nc")
      inquire (file=trim(fname), exist=exists)
      if (.not. exists) then
         print*, "Error => ", trim(fname), " does not exist."
         stop 'topog in topog.f'
      endif
      ib(:) = 1
      ic(:) = imt-2
      ic(2) = jmt-2
      ic(3) = km
      call openfile (fname, iou)
      if (.not. exists) then
         print*, "==> Warning: G_subgrid_bathy data does not exist."
      else
         call getvara ('OCEAN_MASK_MG', iou, ic(1)*ic(2), ib, ic
     &,             tmpij, c1, c0)
      endif
      sg_ocean_mask(2:imt-1,2:jmt-1) = tmpij(1:imt-2,1:jmt-2)

!-----------------------------------------------------------------------
!     apply boundary conditions
!-----------------------------------------------------------------------

      call kmtbc (kmt, imt, jmt)
      call area_volume (kmt, xu, yu, zw)
      cksum1 = checksumi (kmt, imt, jmt, 'Original "kmt" checksum =')

!-----------------------------------------------------------------------
!     iterate on changes to "kmt" until all are consistent
!-----------------------------------------------------------------------

      npass = 0
      n_del_kmt = 0
1000  continue
      npass = npass + 1

!-----------------------------------------------------------------------
!     fill isolated cells or potholes at all levels in kmt
!     filled kmt array is just the maximum of the surrounding kmu levels
!     symmetry conditions are automatic, if set on kmu
!-----------------------------------------------------------------------

      do i=1,imt-1
        do j=1,jmt-1
          kmu(i,j) = min(kmt(i,j), kmt(i+1,j), kmt(i,j+1), kmt(i+1,j+1))
        enddo
      enddo
      n = 0
      write (*,'(/a)') 'Searching for isolated T cells...'
      do i=2,imt-1
        do j=2,jmt-1
          k = max (kmu(i,j), kmu(i-1,j), kmu(i,j-1), kmu(i-1,j-1))
          if (kmt(i,j) .ne. k) n = n + 1
          kmt(i,j) = k
        enddo
      enddo
      if (n .gt. 0) then
        write (*,*) '-> Found ',n,' and filled them in.'
      else
        write (*,*) '-> None Found.'
      endif

      call kmtbc (kmt, imt, jmt)

!-----------------------------------------------------------------------
!     USER INPUT
!     other changes to "kmt" may be made here.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!     ALL USER MODIFICATIONS MUST PRECEDE THESE FINAL CHECKS
!-----------------------------------------------------------------------

      call kmtbc (kmt, imt, jmt)

!-----------------------------------------------------------------------
!     limit the minimum number of levels. kmt_min should be >= 2
!-----------------------------------------------------------------------

      call min_depth (kmt, kmt_min, xu, yu, zw)
      call kmtbc (kmt, imt, jmt)

!-----------------------------------------------------------------------
!     test for "island perimeter violations"
!     channels between different land masses must be at least 2 ocean
!     cells wide to prevent conflicting values of the stream function
!     on these island perimeter points.
!-----------------------------------------------------------------------

      write (stdout,'(/a/)')
     & 'Searching for and correcting island PERIMETER VIOLATIONS...'
      call isleperim (kmt, map, iperm, jperm, iofs, nippts, nisle
     &,               imt, jmt, km, mnisle, maxipp, xu, yu, zw)
      call kmtbc (kmt, imt, jmt)
      call area_volume (kmt, xu, yu, zw)

      write (cktext,'(a,i2,a)') 'Pass # ',npass
     &,                         ' complete.  kmt checksum ='
      cksum2 = checksumi (kmt, imt, jmt, cktext)
      if (cksum1 .ne. cksum2) then
        cksum1 = cksum2
        if (npass .gt. 15) then
          print *,'=>Error: changes to "kmt" may be incompatible'
          print *,'         arbitrarily stopping after 15 passes'
          stop
        endif
        go to 1000
      endif
      print '(/a/)', 'All changes have now been made to kmt'

!-----------------------------------------------------------------------
!     end of iteration on "kmt" changes
!-----------------------------------------------------------------------

      call showmap (map, imt, jmt, linewidth)

!-----------------------------------------------------------------------
!     show the kmt field
!-----------------------------------------------------------------------

      write (stdout,'(/,20x,a/)') ' The "kmt" field follows:'
      call iplot (kmt, imt, imt, jmt)
      call area_volume (kmt, xu, yu, zw)

!---------------------------------------------------------------------
!     compute a topography checksum
!---------------------------------------------------------------------

      cksum1 = checksumi (kmt, imt, jmt, 'Final "kmt" checksum =')

!     construct depth arrays associated with "u" cells

      do j=1,jmt
        kmu(imt,j) = 0
      enddo
      do i=1,imt
        kmu(i,jmt) = 0
      enddo
      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
      do j=1,jmt
        kmu(imt,j) = kmu(2,j)
      enddo

      return
      end

      function checksumi (kmt, im, jm, text)

      implicit none

      integer i, im, jm, jrow, kmt(im,jm)
      character(*) :: text
      real checksumi, cksum

      cksum = 0.0
      do jrow=1,jm
        do i=1,im
          cksum = cksum + i*jrow*kmt(i,jrow)
        enddo
      enddo
      print *, text, cksum
      checksumi = cksum
      return
      end

      subroutine area_volume(kmt, xu, yu, zw)

!-----------------------------------------------------------------------
!     compute surface area and volume of ocean ("t" cells and "u" cells)
!     (note that areas are defined at each level)
!-----------------------------------------------------------------------

      implicit none

      integer k, i, jrow
      real ocnp, tcellv, degrad

      include "size.h"
      include "scalar.h"

      integer kmt(imt,jmt)

      real xu(imt), yu(jmt), zw(km), tcella(km), cst(jmt), dxt(imt)
      real dyt(jmt)

      do k=1,km
        tcella(k) = 0
      enddo
      ocnp   = 0
      tcellv = 0
      pi = 4.0 * atan(1.0)
      degrad = pi / 180.0
      do i=2,imt-1
        dxt(i) = radius * degrad * (xu(i) - xu(i-1))
      enddo
      do jrow = 2,jmt-1
        cst(jrow) = cos(0.5*(yu(jrow)+yu(jrow-1))*degrad)
        dyt(jrow) = radius * degrad * (yu(jrow) - yu(jrow-1))
      enddo

      do jrow=2,jmt-1
        do i=2,imt-1
          if (kmt(i,jrow) .gt. 0) then
            do k=1,kmt(i,jrow)
              tcella(k) = tcella(k) +
     &           cst(jrow)*dxt(i)*dyt(jrow)
            enddo
            tcellv = tcellv + cst(jrow)*dxt(i)*dyt(jrow)*zw(kmt(i,jrow))
            ocnp   = ocnp + float(kmt(i,jrow))
          endif
        enddo
      enddo

      print '(/)'
      print '(a,f8.0)',   'number of ocean T cells =', ocnp
      print '(a,e15.6,a)','surface area (T cells)  =',tcella(1),' cm**2'
      print '(a,e15.6,a)', 'ocean volume (T cells)  =', tcellv,' cm**3'
      print '(/)'

      return
      end

      subroutine kmtbc (kmt, imt, jmt)

!-----------------------------------------------------------------------
!     set lateral boundary conditions on kmt
!     set cyclic, solid wall and symmetry conditions on "t" grid
!-----------------------------------------------------------------------

      implicit none

      integer kmt(imt,jmt), imt, jmt, jrow, i

!     set cyclic conditions on eastern and western boundary

      do jrow=1,jmt
        kmt(1,jrow)   = kmt(imt-1,jrow)
        kmt(imt,jrow) = kmt(2,jrow)
      enddo

!     close the domain at the northernmost row

      do i=1,imt
        kmt(i,jmt) = 0
      enddo

!     close the domain at the southernmost row

      do i=1,imt
        kmt(i,1)   = 0
      enddo

      return
      end

      subroutine setkmt (kmt, xt, yt, imt, jmt, alat1, slon1, elon1
     &,                  alat2, slon2, elon2, num)

!-----------------------------------------------------------------------
!     set the topography mask "kmt(i,j)" = "num" within the area of
!     the parallelogram bounded by vertices:
!     (alat1,slon1), (alat1,elon1), (alat2,slon2), & (alat2,elon2)

!     inputs:

!     xt = longitudes of "t" points in degrees
!     yt = latitudes of "t" points in degrees
!     imt = number of model longitudes
!     jmt = number of model latitudes
!     alat1 = southern latitude of parallelogram (degrees)
!     slon1 = starting longitude of southern edge of parallelogram (deg)
!     elon1 = ending longitude of southern edge of parallelogram (deg)
!     alat2 = northern latitude of parallelogram (degrees)
!     slon2 = starting longitude of northern edge of parallelogram (deg)
!     elon2 = ending longitude of northern edge of parallelogram (deg)
!     num   = number of vertical levels

!     outputs:

!     kmt   = mask of vertical levels on model "t" points
!-----------------------------------------------------------------------

      implicit none

      integer imt, jmt, j1, indp, j2, js, je, i1, i2, is1, ie1
      integer is2, ie2, is, ie, jrow, i, num, kmt(imt,jmt)

      real alat1, alat2, slon1, elon1, slon2, elon2, c1, rdj
      real xt(imt), yt(jmt)

!     convert the four vertices into model indices
!     (js,is1), (js,ie1), (je,is2), (je,ie2)

      j1 = indp (alat1, yt, jmt)
      j2 = indp (alat2, yt, jmt)
      js = min (j1,j2)
      je = max (j1,j2)

      i1  = indp (slon1, xt, imt)
      i2  = indp (elon1, xt, imt)
      is1 = min (i1,i2)
      ie1 = max (i1,i2)

      i1  = indp (slon2, xt, imt)
      i2  = indp (elon2, xt, imt)
      is2 = min (i1,i2)
      ie2 = max (i1,i2)

      is = is1
      ie = ie1

!     fill in the area bounded by (js,is1), (js,ie1), (je,is2), (je,ie2)
!     the nudging of 1.e-5 is to insure that the test case resolution
!     generates the same topography and geometry on various computers.

      c1 = 1.0
      if (js .eq. je) then
        rdj = c1
      else
        rdj = c1/(je-js)
      endif
      do jrow=js,je
        do i=is,ie
          kmt(i,jrow) = num
        enddo
        is = nint(rdj*((jrow-js)*is2 + (je-jrow)*is1) + 1.0e-5)
        ie = nint(rdj*((jrow-js)*ie2 + (je-jrow)*ie1) + 1.0e-5)
      enddo
      return
      end

      subroutine min_depth (kmt, kmt_min, xu, yu, zw)

!     limit the minimum number of levels. kmt_min should be >= 2

      implicit none

      integer n, i, jrow, kmt_min, kmt_shallow

      include "size.h"
      include "stdunits.h"

      integer kmt(imt,jmt)
      real xu(imt), yu(jmt), zw(km)

      write (stdout,'(/a/)')
     & 'Searching for and correcting minimum level violations'
      n = 0
      do i=2,imt-1
        do jrow=jmt-1,2,-1
          if (kmt(i,jrow) .ne. 0 .and. kmt(i,jrow) .lt. kmt_min) then
            n = n + 1
            kmt_shallow = kmt(i,jrow)
            if (zw(kmt(i,jrow)) .lt. 0.5*zw(kmt_min)) then
              kmt(i,jrow) = 0
            else
              kmt(i,jrow) = kmt_min
            endif
          endif
        enddo
      enddo
      if (n .gt. 0) then
        write (stdout,'(a,i5,a/)')  '->Modified', n,' shallow cells'
      else
        write (stdout,'(a/)') '->No modifications needed'
      endif
      return
      end