c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c     begin file topography.f
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 
      program topography_driver

      implicit none
      integer nlat,nlng,nlev

      print *, " "
      print *, "enter grid dimensions: nlng,nlat,nlev "
      read(*,*) nlng,nlat,nlev
      print *, "--->",nlng,nlat,nlev

      call topography(nlng,nlat,nlev)

      end program topography_driver

c---------------------------------------------------------------------

      subroutine topography(nlng,nlat,nlev)

      implicit none

      integer nlng,nlat,nlev,n,iocean,jocean,nsmooth
      double precision ulat(nlng,0:nlat),ulng(nlng,0:nlat)
      double precision h(nlng,0:nlat)
      integer kmt(nlng,nlat),lmask(nlng,0:nlat),smask(nlng,0:nlat)
      double precision dz(nlev),zt(nlev)

      character*80 filename

      print *, " "
      print *, "enter name of direct access file containing"
      print *, "double-precision fields:"
      print *, "ulat(1:nlng,0:nlat),ulng(1:nlng,0:nlat)"
      read(*,'(a80)') filename

      open (53, file=filename, form='unformatted', status='old',
     & access='direct',action='read',recl=(nlat+1)*nlng*8)
      read (53,rec=1) ulat(:,0:nlat)
      read (53,rec=2) ulng(:,0:nlat)
      close(53)
      print *, "--->",filename

      call zcalc(nlev,dz,zt)

      call add_topo(nlng,nlat,nlev,ulng,ulat,h,lmask)

      print *, " "
      print *, "enter: "
      print *, "  0 to continue "
      print *, "  1 to output depth field H(i,j) to direct access file"
      read(*,*) n
      print *, "--->",n
      if(n.eq.1) then
        print *, " "
        print *, "enter filename"
        read(*,'(a80)') filename
        open (54, file=filename, form='unformatted', status='unknown',
     &   access='direct',action='write',recl=nlat*nlng*8)
        write(54,rec=1) h(:,1:nlat)
        close(54)
        open (54, file='lmask.da', form='unformatted', status='unknown',
     &   access='direct',action='write',recl=nlat*nlng*8)
        write(54,rec=1) lmask(:,1:nlat)
        close(54)
      endif

      call min_h(nlng,nlat,nlev,dz,zt,h,kmt)

      print *, " "
      print *, "enter: "
      print *, "  0 to continue "
      print *, "  n to smooth depth field H(i,j) n times"
      read(*,*) nsmooth
      print *, "--->",nsmooth
      if(nsmooth.gt.0) then
        call depth_smooth(nlng,nlat,h,lmask,nsmooth)
      endif

      call compute_kmt(nlng,nlat,nlev,dz,zt,h,kmt)

      call get_smask(nlng,nlat,kmt,lmask,smask)
      open (54, file='smask0.da', form='unformatted', status='unknown',
     &   access='direct',action='write',recl=nlat*nlng*4)
      write(54,rec=1) smask(:,1:nlat)
      close(54)

      print *, " "
      print *, "enter: "
      print *, "  0 to continue "
      print *, "  n to smooth kmt array kmt(i,j) n times"
      read(*,*) nsmooth
      print *, "--->",nsmooth
      if(nsmooth.gt.0) then
        call kmt_smooth(nlng,nlat,kmt,smask,nsmooth)
        call get_smask(nlng,nlat,kmt,lmask,smask)
      endif

      print *, " "
      print *, "enter: "
      print *, "  0 to continue "
      print *, "  1 to remove isolated ridges at all levels from KMT"
      read(*,*) n
      print *, "--->",n
      if (n.eq.1)
     &    call remove_isolated_ridges(nlng,nlat,nlev,kmt,smask)

      call get_smask(nlng,nlat,kmt,lmask,smask)
      open (54, file='smask1.da', form='unformatted', status='unknown',
     &   access='direct',action='write',recl=nlat*nlng*4)
      write(54,rec=1) smask(:,1:nlat)
      close(54)

      print *, " "
      print *, "enter: "
      print *, "  0 to continue "
      print *, "  1 to remove isolated points at all levels from KMT"
      read(*,*) n
      print *, "--->",n
      if (n.eq.1) 
     &    call remove_isolated_points(nlng,nlat,nlev,kmt)

      call get_smask(nlng,nlat,kmt,lmask,smask)
      open (54, file='smask2.da', form='unformatted', status='unknown',
     &   access='direct',action='write',recl=nlat*nlng*4)
      write(54,rec=1) smask(:,1:nlat)
      close(54)

      print *, " "
      print *, "enter: "
      print *, "  0 to continue "
      print *, "  1 to remove disconnected bays that are not"
      print *, "    connected to a selected point"
      read(*,*) n
      print *, "--->",n
      if (n.eq.1) then
2000    print *, " "
        print *, " enter indices (iocean, jocean) "
        print *, " for a select point in the ocean"
        read(*,*) iocean,jocean
        print *, "--->",iocean,jocean
        if (kmt(iocean,jocean).eq.0) then
          print *, " "
          print *, " that is a land point, try again "
          go to 2000
        endif
        call remove_disconnected_bays(iocean,jocean,nlng,nlat,kmt)
      endif

      print *, " "
      print *, "enter: "
      print *, "  0,     to continue "
      print *, 
     &"  n > 0, to force KMT to have a minimum value of n at ocean pts"
      read(*,*) n
      print *, "--->",n
      if(n.gt.0) then
        where(kmt.gt.0) kmt = max(kmt,n)
      endif

      print *, " "
      print *, "enter: "
      print *, "  0 to continue "
      print *, "  1 to output KMT field to direct access file"
      read(*,*) n
      print *, "--->",n
      if(n.eq.1) then
        print *, " "
        print *, "enter filename"
        read(*,'(a80)') filename
        open (55, file=filename, form='unformatted', status='unknown',
     &    access='direct',action='write',recl=nlat*nlng*4)
        write(55,rec=1) kmt
        close(55)
        print *, "--->",filename
      endif

      end subroutine topography

c----------------------------------------------------------------------
c----------------------------------------------------------------------
 
      subroutine zcalc(nlev,dz,zt)
 
c----------------------------------------------------------------------
c     Reads depth profile for number of levels in model grid
c     calculates depths from the profile
c----------------------------------------------------------------------
c      input:
c        nlev = No. of depths in model
 
c      output:
c        zt = depth to center of each level
c        dz = depth profile
 
c      external requirement:
c        'in_depths.dat' = ascii file of dz for nlev levels
c----------------------------------------------------------------------
 
      implicit none

      integer nlev,k
      double precision dz(nlev),zt(nlev)
      character*80 filename

      integer k
c-------------------------------------------------

c     read dz depth profile from file
c
      print *, " "
      print *, " enter name of ascii grid file with depth profile in cm"
      read(*,'(a80)') filename

      open(54,file=filename,status='old',form='formatted')
      do k = 1,nlev
         read (54,*) dz(k)
      enddo
c     print *,k,dz
      close (54)
 
      print *, '--->',filename
c     calculate z from dz
 
      zt(0) = 0.
      zt(1)=.5*dz(1)
      do k=2,nlev
         zt(k)=zt(k-1) + 0.5*(dz(k)+dz(k-1))
      enddo
      print *, " "
      print *, "  k       dz(m)       zt(m)"
 11   format(x,i3,x,f12.3,x,f12.3)
      do k=0,nlev
         write(*,11) k,dz(k)*0.01,zt(k)*0.01
      enddo
 
      return
 
      end
 
c---------------------------------------------------------------------
c---------------------------------------------------------------------
c     subroutine linking the mesh and the topography
c---------------------------------------------------------------------
c---------------------------------------------------------------------
 
      subroutine add_topo(nlng,nlat,nlev,ulng,ulat,h,lmask)

      implicit none

      integer nlng,nlat,nlev
      double precision ulat(nlng,0:nlat),ulng(nlng,0:nlat)
      double precision h(nlng,0:nlat)
      integer ntlng,ntlat,lmask(nlng,0:nlat)
      parameter (ntlng=10800,ntlat=5400)
      real rtopo (ntlng,ntlat)

      double precision x(3,5),ullng(4),ullat(4),pi
      double precision real_lat,real_lng,hh,hland,hocean,xp,yp,zp
      double precision xv1,xv2,xpv,yv1,yv2,ypv,zv1,zv2,zpv
      integer i,j,k,i1,i2,im1,j1,l,kp1,ln,la
      integer lngmax,lngmin,latmax,latmin,nbland,nbpocean,nbpland
      logical incell, land
      character*80 fname

      pi = 3.141592653589793d0

 
c-----------------------------------------------------------------
c     reading topography...
c     The  integer values (depth at each point) are in meters, positive
c     above sealevel, neg. below.
c     The resolution is 2 min. x 2min.
c     The data is organized in full 180 deg. latitude circles
c     beginning at 90N (constant val.) and stepping south to 89 55.0S
c-----------------------------------------------------------------

      print *, " "
      print *, " reading topo data..."

c  dataset below archived at /project/oce/bryan/etopo_merge_bedmap.da
      fname = '/cgd/oce/yeager/obs/topo/etopo_merge_bedmap.da'
      open (55, file=fname, form='unformatted', status='old',
     & access='direct',action='read',recl=ntlat*ntlng*4)
      read (55,rec=1) rtopo
      close(55)
 
      print *, " "
      print *, "read rtopo. r(5000,2500) = ",rtopo(5000,2500) 
      print *, "calculating depth field..."
 
c----------------------------------------------------------
c     is there land at T-points?
c----------------------------------------------------------

      where(ulng.gt.pi) ulng=ulng-2.0*pi
 
      nbland = 0
      do j=1,nlat
         do i=1,nlng
            if (i.eq.1) then
               im1 = nlng
            else
               im1 = i-1
            endif
 
c     what are the coordonates of the 4 edges of the cell?
            do k=1,4
               if (k.eq.1.or.k.eq.2) then
                  j1 = j-1
               else
                  j1 = j
               endif
               if (k.eq.1.or.k.eq.4) then
                  i1 =  im1
               else
                  i1 = i
               endif
               x (1,k) = cos(ulng(i1,j1))*cos(ulat(i1,j1))
               x (2,k) = sin(ulng(i1,j1))*cos(ulat(i1,j1))
               x (3,k) = sin(ulat(i1,j1))
               ullng (k) = ulng(i1,j1)
               ullat (k) = ulat(i1,j1)
            enddo
 
 
c----------------------------------------------------------
c     what are the (long-lat) boundaries of the cell ?
c----------------------------------------------------------

            lngmax = (180.0*30+int(max(ulng(i,j),ulng(im1,j),
     &           ulng(i,j-1),ulng(im1,j-1))*180.*30./pi))+1
            lngmin = (180.0*30+int(min(ulng(i,j),ulng(im1,j),
     &           ulng(i,j-1),ulng(im1,j-1))*180.*30./pi))+1
            latmin = (90*30 - int(max(ulat(i,j),ulat(i,j-1),
     &           ulat(im1,j),ulat(im1,j-1))*30.*180./pi)) + 1
            latmax = (90*30 - int(min(ulat(i,j),ulat(i,j-1),
     &           ulat(im1,j),ulat(im1,j-1))*30.*180./pi)) + 1
            if (float(lngmax-lngmin)/float(ntlng).gt.0.25) then
               l = lngmin
               lngmin = lngmax
               lngmax = l+ntlng
            endif
 
c----------------------------------------------------------
c     calculation of depth
c----------------------------------------------------------
 
            hland = 0.
            hocean = 0.
            incell = .true.
            land = .false.
            nbpocean =0
            nbpland =0
            do j1=latmin,latmax
               do i2=lngmin,lngmax
 
                  i1 = mod(i2,ntlng)
                  incell = .true.
                  real_lng = float(i1-1)*2.*pi/float(ntlng)-pi
                  real_lat = pi/2.-float(j1-1)*pi/float(ntlat)
                  xp = cos(real_lng)*cos(real_lat)
                  yp = sin(real_lng)*cos(real_lat)
                  zp = sin(real_lat)
 
                  do k=1,4
                     kp1 = mod(k,4)+1
                     xv1 = x(1,k) - xp
                     yv1 = x(2,k) - yp
                     zv1 = x(3,k) - zp
                     xv2 = x(1,kp1) - xp
                     yv2 = x(2,kp1) - yp
                     zv2 = x(3,kp1) - zp
                     xpv = yv1*zv2 - zv1*yv2
                     ypv =-xv1*zv2 + zv1*xv2
                     zpv = xv1*yv2 - yv1*xv2
                     incell = incell.and.(xpv*xp+ypv*yp+zpv*zp.ge.0)
                  enddo
 
                  if (incell) goto 20
                  goto 10
 20               continue
                  if (rtopo(i1,j1).ge.0.0) then
                     nbpland  = nbpland +1
                  else
                     hocean = hocean + rtopo(i1,j1)
                     nbpocean  = nbpocean +1
                  endif
 
 
 10               continue
 
               enddo
            enddo
 
 
            if (nbpland.gt.(nbpocean)) then
               hh =  0.0
            else
               if (nbpocean+nbpland.eq.0) then
c                 print *,'nbpocean+nbpland.eq.0: ',i,j
                  hh=0.0
                  do k=1,4
                     ln = (-180.0*30-int(ullng(k)*180.*30./pi))+1
                     la = (90*12 - int(ullat(k)*12.*180./pi)) + 1
c                     print *,'     ln,la,ullng,ullat = ',ln, la,
c    &                        ullng(k),ullat(k)
                     hh = hh+ rtopo(ln,la)
c                    print *,'     hh = ',hh
                     nbpocean = 1
                  enddo
                  hh = hh/4.
c                 print *,i,j,tlng(i,j)*180./pi,tlat(i,j)*180./pi,hh
               else
C                 hh = (hocean+hland)/float(nbpland+nbpocean)
                  hh = (hocean)/float(nbpocean)
               endif
            endif
            if (hh.ge.0.0) then
               nbland = nbland +1
            endif
            h(i,j) = -hh*100.0
c	if there is any land in this grid cell, don't smooth later
            if (nbpland.gt.0) then
              lmask(i,j) = 1
            else
              lmask(i,j) = 0
            endif
 
	 enddo
c        print *,'nb de points :',nbland,' pour j ',j
      enddo
 
      end
 
c----------------------------------------------------------------------

      subroutine min_h(nlng,nlat,nlev,dz,zt,h,kmt)

c-----------------------------------------------------------------------
c     enforce lower bounds on h(i,j) depth field.
c-----------------------------------------------------------------------

      implicit none

      integer nlat,nlng,nlev
      double precision h(nlng,0:nlat)
      integer kmt(nlng,nlat)
      double precision dz(nlev),zt(nlev)
      double precision hmin
      integer i,j,k

      print *, " "
      print *, " modifying h field ..."

      print *, " "
      print *, "enter shallow water cuttoff (m)' "
      print *, "  (points with depths shallower than this cuttoff"
      print *, "   will be automatically set to land)"
      read(*,*) hmin
      print *,'--->',hmin
c     hmin=5.
      hmin = hmin*100.   ! convert to cm
 
      do j = 1,nlat
         do i = 1,nlng
            if(h(i,j) .lt. hmin) h(i,j) = 0.0
            if((h(i,j) .ge. hmin) .and. (h(i,j) .lt. zt(1)))
     &         h(i,j) = zt(1)
         enddo
      enddo
 
      end  ! min_h
 
c-----------------------------------------------------------------------

      subroutine compute_kmt(nlng,nlat,nlev,dz,zt,h,kmt)
 
c-----------------------------------------------------------------------
c     find KMT field given H field
c
c     input:
c       H = depth field (cm)
c       Z = model depths (cm)
c       nlng, nlat, nlev = No. of longitudes, latitudes, depths in model
c
c     output:
c       KMT field
c-----------------------------------------------------------------------

      implicit none

      integer nlat,nlng,nlev
      double precision h(nlng,0:nlat)
      integer kmt(nlng,nlat)
      double precision dz(nlev),zt(nlev)
      integer i,j,k
 
c***********************************************************************

      print *, " "
      print *, " computing kmt field ..."
 
      do j = 1,nlat
         do i = 1,nlng
            kmt(i,j) = 0
            if(h(i,j) .lt. zt(1))  kmt(i,j) = 0
            if(h(i,j) .ge. zt(nlev)) kmt(i,j) = nlev
         enddo
      enddo
 
      do k = 1,nlev-1
         do j = 1,nlat
            do i = 1,nlng
               if((h(i,j) .ge. zt(k)) .and. (h(i,j) .lt. zt(k+1)))
     &            kmt(i,j) = k
            enddo
         enddo
      enddo
 
      do i = 1,nlng
         kmt(i,1)    = 0
         kmt(i,nlat) = 0
      enddo
 
      end  ! compute_kmt

c---------------------------------------------------------------------

      subroutine depth_smooth(nlng,nlat,h,lmask,ns)

c---------------------------------------------------------------------
c     smooths H(i,j) depth field using a predefined weight array (wts).
c---------------------------------------------------------------------

      implicit none

      integer i,j,n,nlng,nlat,ns,wtsum,ii,jj
      integer ist,jst,ien,jen,iw,jw
      integer wts(3,3), nmiss, lmask(nlng,0:nlat)
      double precision h(nlng,0:nlat), hnew(nlng,0:nlat)
      double precision hsum

      data wts(1:3,1) /1,2,1/
      data wts(1:3,2) /2,8,2/
      data wts(1:3,3) /1,2,1/

      print *, " "
c     print *, "smoothing depth array h(i,j) ",ns," times"

      hnew = h

      do n=1,ns
      do j=1,nlat
        do i=1,nlng    
        if (h(i,j).ne.0.0.and.lmask(i,j).eq.0) then
          wtsum = 0
	  hsum = 0.0
	  if (j.eq.1) then 	
		jst=0 
	  else 
		jst=-1
	  end if
	  if (j.eq.nlat) then 
		jen=0 
	  else 
		jen=1
	  end if
	  do ii=-1,1
	  do jj=jst,jen
	    if (i.eq.1.and.ii.eq.-1) then 
		iw=nlng 
	    else if (i.eq.nlng.and.ii.eq.1) then 
		iw=1 
	    else 
		iw=i+ii 
	    end if
	    jw = j+jj
	    if (h(iw,jw).ne.0.0.and.lmask(iw,jw).eq.0) then
	      wtsum = wtsum+wts(ii+2,jj+2)
	      hsum = hsum+h(iw,jw)*wts(ii+2,jj+2)
	    end if
	  enddo
	  enddo
!  nmiss = sum(wts) - wtsum
!  if (nmiss.ne.0) then
!     hsum = hsum+h(i,j)*nmiss
       !  end if
!  hnew(i,j) = hsum/sum(wts)
	  hnew(i,j) = hsum/wtsum
        end if
        enddo    ! i
      enddo    ! j
      h = hnew
      enddo    ! n

      end subroutine depth_smooth

c---------------------------------------------------------------------

      subroutine kmt_smooth(nlng,nlat,kmt,smask,ns)

c---------------------------------------------------------------------
c     smooths KMT(i,j) depth field using a predefined weight array (wts).
c---------------------------------------------------------------------

      implicit none

      integer i,j,n,nlng,nlat,ns,wtsum,ii,jj
      integer ist,jst,ien,jen,iw,jw
      integer wts(3,3), nmiss, smask(nlng,0:nlat),kmt(nlng,nlat)
      integer kmtsum,kmtnew(nlng,nlat)

      data wts(1:3,1) /1,1,1/
      data wts(1:3,2) /1,1,1/
      data wts(1:3,3) /1,1,1/

      print *, " "
c     print *, "smoothing kmt array kmt(i,j) ",ns," times"

      kmtnew = kmt

      do n=1,ns
      do j=1,nlat
        do i=1,nlng    
        if (kmt(i,j).ne.0.and.smask(i,j).ne.0) then
          wtsum = 0
	  kmtsum = 0
	  if (j.eq.1) then 	
		jst=0 
	  else 
		jst=-1
	  end if
	  if (j.eq.nlat) then 
		jen=0 
	  else 
		jen=1
	  end if
	  do ii=-1,1
	  do jj=jst,jen
	    if (i.eq.1.and.ii.eq.-1) then 
		iw=nlng 
	    else if (i.eq.nlng.and.ii.eq.1) then 
		iw=1 
	    else 
		iw=i+ii 
	    end if
	    jw = j+jj
	    if (kmt(iw,jw).ne.0) then
	      wtsum = wtsum+wts(ii+2,jj+2)
	      kmtsum = kmtsum+kmt(iw,jw)*wts(ii+2,jj+2)
	    end if
	  enddo
	  enddo
	  kmtnew(i,j) = kmtsum/wtsum
        end if
        enddo    ! i
      enddo    ! j
      kmt = kmtnew
      enddo    ! n

      end subroutine kmt_smooth

c---------------------------------------------------------------------

      subroutine get_smask(nlng,nlat,kmt,lmask,smask)

c---------------------------------------------------------------------
c     flags (i,j) where kmt needs to be smoothed by inspecting
c	i,j and nearest neighbors to N,S,E,W:
c			(i,jn)
c		(iw,j)  (i,j)	(ie,j)
c			(i,js)
c  KEY:
c	0  = land OR monotonic kmt change in 5-pt stencil
c	4  = kmt saddle point
c	3  = local bathymetric maximum (kmt minimum) in both i and j
c      -3  = local bathymetric hole (kmt maximum) in both i and j
c	2  = local bathymetric maximum (kmt minimum) in  i-direction
c      -2  = local bathymetric hole (kmt maximum) in i-direction
c	1  = local bathymetric maximum (kmt minimum) in  j-direction
c      -1  = local bathymetric hole (kmt maximum) in j-direction
c---------------------------------------------------------------------

      implicit none

      integer i,j,n,nlng,nlat,tmp1,tmp2,tmp3,tmp4
      integer iw,js,ie,jn
      integer lmask(nlng,0:nlat),smask(nlng,0:nlat),kmt(nlng,nlat)

      do j=2,nlat-1
        do i=1,nlng    
        smask(i,j) = 0
        if (lmask(i,j).eq.0) then
          js = j-1
          jn = j+1
          iw = i-1
          ie = i+1
	  if (i.eq.1) iw=nlng 
	  if (i.eq.nlng) ie=1
          tmp1 = kmt(i,j)-kmt(iw,j)
          tmp2 = kmt(ie,j)-kmt(i,j)
          tmp3 = kmt(i,j)-kmt(i,js)
          tmp4 = kmt(i,jn)-kmt(i,j)
          if (tmp1*tmp2.lt.0.and.tmp3*tmp4.lt.0) then 
             if (tmp1*tmp3.lt.0) then
               smask(i,j) = 4
             else if (tmp1.lt.0) then
	       smask(i,j) = 3
             else
	       smask(i,j) = -3
             endif
          else if (tmp1*tmp2.lt.0.and.tmp3*tmp4.ge.0) then
             if (tmp1.lt.0) then
               smask(i,j) = 2
             else
               smask(i,j) = -2
             endif
          else if (tmp1*tmp2.ge.0.and.tmp3*tmp4.lt.0) then 
             if (tmp3.lt.0) then
               smask(i,j) = 1
             else
               smask(i,j) = -1
             endif
          else    
             smask(i,j) = 0
          endif
        endif
        enddo
      enddo

      end subroutine get_smask

      subroutine remove_isolated_ridges(nlng,nlat,nlev,kmt,smask)

c---------------------------------------------------------------------
c     removes one-grid cell isolated ridges at all levels.
c     (land points sandwiched between ocean on lateral boundaries:
c           land to N and S, or, land to E and W)
c---------------------------------------------------------------------

      implicit none

      integer i,j,k,nlng,nlat,nlev
      integer im1,ip1,jm1,jp1,nremoved
      integer kmt(nlng,nlat), kmtx(nlng,nlat), smask(nlng,0:nlat)
      logical land,lande,landw,landn,lands

      print *, " "
      print *, "removing isolated ridges... "

      kmtx = kmt   ! copy kmt into temp array

      do k = 1,nlev

1000    nremoved = 0

        do i=1,nlng

          if(i.eq.1) then
            im1 = nlng
          else
            im1 = i-1
          endif
          if(i.eq.nlng) then
            ip1 = 1
          else
            ip1 = i+1
          endif

        do j=2,nlat-1    ! j=1 and j=nlat are assumed to be land

          jm1 = j-1 
          jp1 = j+1 

c         ...note: the following must be modified for tripole grids

          land   = k.gt.kmtx(i,j)
          lande  = k.gt.kmtx(ip1,j)
          landw  = k.gt.kmtx(im1,j)
          landn  = k.gt.kmtx(i,jp1)
          lands  = k.gt.kmtx(i,jm1)


          if ( land  .and. (smask(i,j).ne.0) .and.
     &         ((.not. lande .and. .not. landw) .or. 
     &         (.not. landn .and. .not. lands))) then
            kmtx(i,j) = kmtx(i,j) + 1   ! remove cell at level k
            nremoved = nremoved + 1
          endif

        enddo  ! j
        enddo  ! i 

c       print *, 'points removed this pass: ', k, nremoved
        if (nremoved .ne. 0) go to 1000

      enddo    ! k

      kmt = kmtx    ! overwrite kmt 

      end subroutine remove_isolated_ridges

 
c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine remove_isolated_points(nlng,nlat,nlev,kmt)

c---------------------------------------------------------------------
c     removes isolated points at all levels.
c     (points sandwiched between land on lateral boundaries:
c           land to N and S, or, land to E and W)
c---------------------------------------------------------------------

      implicit none

      integer i,j,k,nlng,nlat,nlev
      integer im1,ip1,jm1,jp1,nremoved
      integer kmt(nlng,nlat), kmtx(nlng,nlat)
      logical land,lande,landw,landn,lands

      print *, " "
      print *, "removing isolated points... "

      kmtx = kmt   ! copy kmt into temp array

      do k = nlev,1,-1

1000    nremoved = 0

        do i=1,nlng

          if(i.eq.1) then
            im1 = nlng
          else
            im1 = i-1
          endif
          if(i.eq.nlng) then
            ip1 = 1
          else
            ip1 = i+1
          endif

        do j=2,nlat-1    ! j=1 and j=nlat are assumed to be land

          jm1 = j-1 
          jp1 = j+1 

c         ...note: the following must be modified for tripole grids

          land   = k.gt.kmtx(i,j)
          lande  = k.gt.kmtx(ip1,j)
          landw  = k.gt.kmtx(im1,j)
          landn  = k.gt.kmtx(i,jp1)
          lands  = k.gt.kmtx(i,jm1)


          if ( .not.land  .and. 
     &         ((lande .and. landw) .or. (landn .and. lands))
     &       ) then
            kmtx(i,j) = kmtx(i,j) - 1   ! remove cell at level k
            nremoved = nremoved + 1
          endif

        enddo  ! j
        enddo  ! i 

c       print *, 'points removed this pass: ', k, nremoved
        if (nremoved .ne. 0) go to 1000

      enddo    ! k

      kmt = kmtx    ! overwrite kmt 

      end subroutine remove_isolated_points

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine remove_disconnected_bays(iocean,jocean,imt,jmt,KMTX)

c-----------------------------------------------------------------------
c     remove isolated lakes, bays, single points from grid (KMT field)
c     note: this routine assumes land along northern and southern
c     boundaries: KMTX(:,1) = KMTX(:,jmt) = 0
c-----------------------------------------------------------------------

      implicit none

      integer i,j,iocean,jocean,imax,imt,jmt
      logical, dimension(imt,jmt) :: MASKC
      integer, dimension(imt,jmt) :: KMTX,KMUX,KINT
      integer, dimension(imt,jmt) :: CFN,CFS,CFE,CFW,I1,I2
      real*4, dimension(imt,jmt) :: CN,CS,CE,CW,C0

      print *, " "
      print *, " removing disconnected bays..."

      do i=1,imt
      do j=1,jmt
        I1(i,j) = i
        I2(i,j) = j
      enddo
      enddo
      KINT = I1 + (I2-1)*imt

      KMUX = min(KMTX,cshift(KMTX,+1,1))
      KMUX = min(KMUX,cshift(KMUX,+1,2))

      CFE = KMUX + cshift(KMUX,-1,2)
      CFW = cshift(CFE,-1,1)
      CFN = KMUX + cshift(KMUX,-1,1)
      CFS = cshift(CFN,-1,2)

      where (CFE.gt.0)
         CFE = 1
      elsewhere
         CFE = imt*jmt + 1
      endwhere

      where (CFW.gt.0)
         CFW = 1
      elsewhere
         CFW = imt*jmt + 1
      endwhere

      where (CFN.gt.0)
         CFN = 1
      elsewhere
         CFN = imt*jmt + 1
      endwhere

      where (CFS.gt.0)
         CFS = 1
      elsewhere
         CFS = imt*jmt + 1
      endwhere

      MASKC = .true.
      where (KMTX.eq.0)
         MASKC = .false.
         KINT = imt*jmt + 1
      endwhere

      imax = (imt+jmt)*3
      do i = 1,imax
         CE = cshift(KINT,+1,1)*(1.0*CFE)
         CW = cshift(KINT,-1,1)*(1.0*CFW)
         CN = cshift(KINT,+1,2)*(1.0*CFN)
         CS = cshift(KINT,-1,2)*(1.0*CFS)
         C0 = KINT*1.0
         where (MASKC)
            KINT = nint(min(CE,CW,CN,CS,C0))
         endwhere
      enddo

c     print *, 'KMTX(iocean,jocean) =', KMTX(iocean,jocean)
c     print *, 'KINT(iocean,jocean) =', KINT(iocean,jocean)

      where (KINT.ne.KINT(iocean,jocean))
         KINT = 0
      elsewhere
         KINT = KMTX
      endwhere

      print *, ' '
      print *, 'surface count =',count(KINT.ne.0)

      KMTX = KINT

      end subroutine remove_disconnected_bays

 
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c     end file topography.f
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||