! source file: /Users/oschlies/UVIC/master/source/mom/topog.F
      subroutine topog (kmt, map, xt, yt, zt, xu, yu, zw, imt2,jmt2,km2)

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

!     A base kmt field is constructed using one of the following:

!      (A) Construct a kmt field using one of the following:
!        -Dscripps_kmt : interpolate Scripps 1 degree topography to
!                        "kmt" on MOM grid
!        -Dread_my_kmt : read in your own "kmt" from a file "kmt.dta"
!        -Didealized_kmt : generate an "idealized" not very realistic
!                          topography without reading the database
!        -Drectangular_box : sets "kmt" to a flat bottom rectangular box

!       Options:
!        -Dcyclic        : sets 1 conditons in zonal direction
!                          (otherwise solid_walls is assumed)
!        -Dflat_bottom   : sets "kmt" over ocean to maximum levels (km)

!     Problems may exist in the base kmt field. They are fixed by
!     iterating over the following sequence of operations.

!      (B) Fill isolated cells (optional)

!       Option:
!        -Dfill_isolated_cells : These are places where "T" cells are
!                                 deeper than nearest neighbors. Also
!                                 fills in "T" cells with land at all
!                                 four corners.

!      (C) Add USER CHANGES to kmt field (optional)

!      (D) Limit minimum levels based on parameter "kmt_min"
!          "kmt_min" is the minimum number of vertical levels
!           permitted at any ocean "T" cell.  kmt_min must be at >= 2.

!           Options:
!           -Dfill_shallow: makes all "t" cells for which
!                           kmt(i,jrow) < kmt_min into land.
!           -Ddeepen_shallow: deepens such cell to depth kmt_min.
!           -Dround_shallow: changes depth to the nearer of 0 and
!                            zw(kmt_min)

!      (E) Fix perimeter violations (cells in the ocean perimeter of
!                                    two different land masses)
!          Options:
!           -Dfill_perimeter_violations : build a land bridge between
!                                         these land masses (default)
!           -Dwiden_perimeter_violations: widen straits between two
!                                         different land masses to a
!                                         minimum of two "t" cells

!-----------------------------------------------------------------------
!     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|

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

      include "stdunits.h"
      include "size.h"
      include "isleperim.h"
      real tmpij(imt,jmt)
      integer ib(10), ic(10)

!     kmt_min is the minimum permitted number of levels at an ocean
!     "T" cell.  kmt_min must be at least 2.

      parameter (kmt_min = 2)

      dimension kmt(imt,jmt), kmu(imt,jmt), xt(imt), yt(jmt), zt(km)
      dimension xu(imt), yu(jmt), zw(km)

      dimension map(imt,jmt)
      dimension iperm(maxipp)
      dimension jperm(maxipp)
      dimension nippts(mnisle)
      dimension iofs(mnisle)

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

      character(120) :: fname, new_file_name
      logical exists

!     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 ("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
      ic(2) = jmt
      call openfile (fname, iou)
      call getvara ('kmt', iou, imt*jmt, ib, ic, tmpij, c1, c0)
      call closefile (iou)
      kmt(:,:) = tmpij(:,:)

      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.'

!-----------------------------------------------------------------------
!     check that only one method of generating topography was chosen
!-----------------------------------------------------------------------

      if (ncase .eq. 0) then
        write (stdout,'(/a/a/a/)')
     &  '=>Error: One of the following options must be enabled:'
     &, '      idealized_kmt, scripps_kmt, read_my_kmt, rectangular_box'
     &, '      etopo_kmt'
        stop '=>topog'
      elseif (ncase .gt. 1) then
        write (stdout,'(/a/a/)')
     &  '=>Error: Only one of the following options may be enabled:'
     &, '      idealized_kmt, scripps_kmt, read_my_kmt, rectangular_box'
        stop '=>topog'
      endif

!-----------------------------------------------------------------------
!     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 consistant
!-----------------------------------------------------------------------

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

!-----------------------------------------------------------------------
!     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 incompatable'
          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
!-----------------------------------------------------------------------

      print '(/a/)', 'Use -Dskip_island_map eliminate the following map'
      call showmap (map, imt, jmt, linewidth)

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

      print '(/a/)', 'Use -Dskip_kmt_map to eliminate the following map'
      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 =')

      return
      end

      function checksumi (kmt, im, jm, text)
      dimension kmt(im,jm)
      character(*) :: text
      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)
!-----------------------------------------------------------------------

      include "size.h"
      include "scalar.h"
      dimension kmt(imt, jmt), xu(imt), yu(jmt), zw(km)
      dimension tcella(km)
      dimension cst(jmt)
      dimension dxt(imt), 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 1, solid wall and symmetry conditions on "t" grid
!-----------------------------------------------------------------------

      dimension kmt(imt,jmt)

!     set 1 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 notthern edge of parallelogram (deg)
!     elon2 = ending longitude of notthern edge of parallelogram (deg)
!     num   = number of vertical levels

!     outputs:

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

!     based on code by: R. C. Pacanowski
!-----------------------------------------------------------------------

      dimension kmt(imt,jmt), 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

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

      dimension kmt(imt,jmt)
      dimension 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