#include "options.inc"

      subroutine grids
c
c=======================================================================
c     set up a staggered "B" grid for MOM and compute grid related
c     variables
c     original author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
c     SPFLAME version : c.eden
c=======================================================================
c
      use spflame_module
      implicit none
      real y,dfdy,amix
      integer is,ie,js,je,i,j,k,jm1,jp1,im1,ip1,io
      integer :: ii,jj
      real :: scal , scal_length, fxa , fxb, degtcm, distance, dstnc
      real, allocatable:: glat(:,:)

      degtcm   = radius/180*pi

      if (mother<0) then
c
c-----------------------------------------------------------------------
c      read in or set up the grid for the mother domain
c-----------------------------------------------------------------------
c
       if (enable_simple_grid)  then
        call grid_template()
       else
        call read_grid
       endif

      endif
c
c-----------------------------------------------------------------------
c     exchange grid over sub domains
c-----------------------------------------------------------------------
c
      call exchg_sub_grid
c
c-----------------------------------------------------------------------
c     convert grid resolution to cm
c-----------------------------------------------------------------------
c
      do j=1,jmt
        dyt(j) = dytdeg(j)*degtcm
        dyu(j) = dyudeg(j)*degtcm
      enddo
      do i=1,imt
        dxt(i) = dxtdeg(i)*degtcm
        dxu(i) = dxudeg(i)*degtcm
      enddo
      if (cyclic) then
       dxt(1)   = dxt(imt-1)
       dxt(imt) = dxt(2)
       dxu(1)   = dxu(imt-1)
       dxu(imt) = dxu(2)
      endif
c
c-----------------------------------------------------------------------
c     compute all quantities derived from the grid spacings
c-----------------------------------------------------------------------
c
      do k=1,km
        c2dzt(k) = 2.*dzt(k)
        dzt2r(k) = 1./c2dzt(k)
      enddo
      dzwr(km)  = 1./dzw(km)
      dzw2r(km) = 0.5/dzw(km)
      do k=1,km
        dzwr(k-1)    = 1./dzw(k-1)
        dzw2r(k-1)   = 0.5/dzw(k-1)
      enddo
      do k=1,km
        dztr(k)  = 1./dzt(k)
        rho0dztr(k) = 1./(rho0*dzt(k))
        dztur(k) = 1./(dzw(k-1)*dzt(k))
        dztlr(k) = 1./(dzw(k)*dzt(k))
      enddo
c
      if (enable_beta_plane .or. enable_f_plane) then
       do j=1,jmt
        phi(j)   = yu(beta_plane_j0)/180*pi
	phit(j)  = yt(beta_plane_j0)/180*pi
       enddo
      else
       do j=1,jmt
        phi(j)   = yu(j)/180*pi
	phit(j)  = yt(j)/180*pi
       enddo
      endif

      do j=1,jmt
        dytr(j)  = 1./dyt(j)
        dyt2r(j) = 0.5/dyt(j)
        dyt4r(j) = 0.25/dyt(j)
        dyur(j)  = 1./dyu(j)
        dyu2r(j) = 0.5/dyu(j)
        dyu4r(j) = 0.25/dyu(j)
        cst(j)   = cos(phit(j))
        csu(j)   = cos(phi (j))
        sine(j)  = sin(phi(j))
        cstr(j)     = 1./cst(j)
        csur(j)     = 1./csu(j)
        tng(j)      = sine(j)/csu(j)
	cstdytr(j)  = 1./(cst(j)*dyt(j))
	cstdyt2r(j) = cstdytr(j)*0.5
        rho0csudyur(j)  = 1./(rho0*csu(j)*dyu(j))
        csudyur(j) = 1./(csu(j)*dyu(j))
        csudyu2r(j) = 0.5/(csu(j)*dyu(j))
        cst_dytr(j) = cst(j)/dyt(j)
        csu_dyur(j) = csu(j)/dyu(j)
      enddo

      do i=1,imt
        dxtr(i)  = 1./dxt(i)
        dxt2r(i) = 0.5/dxt(i)
        dxt4r(i) = 0.25/dxt(i)
        dxur(i)  = 1./dxu(i)
        dxu2r(i) = 0.5/dxu(i)
        dxu4r(i) = 0.25/dxu(i)
      enddo
      do i=2,imt-1
        dxmetr(i) = 1./(dxt(i) + dxt(i+1))
      enddo
      do i=1,imt
        duw(i) = (xu(i) - xt(i))*degtcm
      enddo
      do i=1,imt-1
	due(i) = (xt(i+1) - xu(i))*degtcm
      enddo
      if (cyclic) then
       due(imt) = due(2)
      else
       due(imt) = due(imt-1)
      endif
      do j=1,jmt
        dus(j) = (yu(j) - yt(j))*degtcm
      enddo
      do j=1,jmt-1
	dun(j) = (yt(j+1) - yu(j))*degtcm
      enddo
      dun(jmt) = dun(jmt-1)
c
c-----------------------------------------------------------------------
c     build coefficients to minimize advection and diffusion computation 
c-----------------------------------------------------------------------
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      do j=js-1,je+1
        do i=is-1,ie+1
          ip1=min(i+1,imt)
	  rho0csudxur(i,j) = rho0r*csur(j)*dxur(i)
	  csudxu2r(i,j)    = csur(j)*dxur(i)*0.5
	  am_csudxtr(i,j)  = am*csur(j)*dxtr(ip1)
	enddo
      enddo

      do j=js-1,je+1
        do i=is-1,ie+1
	  cstdxtr(i,j)    = cstr(j)*dxtr(i)
	  cstdxt2r(i,j)   = cstr(j)*dxtr(i)*.5
	  ah_cstdxur(i,j) = ah*cstr(j)*dxur(i)
	enddo
      enddo
c
c-----------------------------------------------------------------------
c     coriolis factors
c-----------------------------------------------------------------------
c
      if (enable_beta_plane .or. enable_f_plane) then
c
c      beta plane with f = f0 + beta*y where f0 is at yu(1)
c      if f_plane then beta = 0
c
       if (enable_f_plane) then
        dfdy   = 0.
        if (my_pe==0) then
         print*,''
         print*,' using a f-plane approximation '
         print*,' referenced to latitude ',yu(beta_plane_j0)
         print*,' f0=',2*omega*sine(beta_plane_j0),' 1/s'
         print*,''
        endif
       else
        dfdy   = 2.*omega*csu(beta_plane_j0)/radius
        if (my_pe==0) then
         print*,''
         print*,' using a beta-plane approximation '
         print*,' referenced to latitude ',yu(beta_plane_j0)
         print*,' f0=',2*omega*sine(beta_plane_j0),' 1/s'
         print*,' beta=',dfdy,' 1/s/cm'
         print*,''
        endif
       endif
       do j=js_pe,je_pe
        y            = (yu(j)-yu(beta_plane_j0))*degtcm
        cori(:,j,1) = 2.*omega*sine(beta_plane_j0) + y*dfdy
        cori(:,j,2) = -(2.*omega*sine(beta_plane_j0) + y*dfdy)
       enddo
        
      elseif (enable_rotated_grid) then

       if (my_pe==0) then
         print*,''
         print*,' using a rotated grid '
         print*,''
       endif

       if (mother>=0) then
        print*,' rotated grid for sub domains not yet'
        print*,' implemented' ! but possible,
!       exchange glat over sub domains here, but interpolation
c       of glat is no good ...
        call halt_stop('in grids (glat)')
       endif

       allocate(glat(imt,jmt))
       if (my_pe==0) print*,
     &  ' =>Reading a "glat" field from binary file glat.dta'
       call getunit (io, 'glat.dta', 'usr ieee')
       read (io) ! flame id
c       skip header record and read dimensional info
       read (io) ! iotext
       read (io) i, j, k
       if (i/=imt .or. j /= jmt .or. k/= km) then
c      check i,j,k
        if (my_pe==0) then
         print*,' Error: '
         print*,' read imt,jmt,km=',i,j,k
         print*,' but input from namelist was ',imt,jmt,km
         call sub_flush(6)
        endif
        call halt_stop(' reading glat file ')
       endif
       read (io) ! iotext
       read (io) glat
       close(io)
       if (my_pe==0) print*,' done'

       do j=js_pe,je_pe
	do i=is_pe,ie_pe
         cori(i,j,1) = 2.0*omega*sin(glat(i,j)/180*pi)
	 cori(i,j,2) = -cori(i,j,1)
	enddo
       enddo
       deallocate(glat)

      else

       do j=js_pe,je_pe
        cori(:,j,1) = 2.*omega*sine(j)
        cori(:,j,2) = -2.*omega*sine(j)
       enddo

      endif
c
c-----------------------------------------------------------------------
c     metric diffusion factors and metric advection factors
c-----------------------------------------------------------------------
c
      amix = sqrt(abs(ambi))
      if (enable_beta_plane .or. enable_f_plane) then
       do j=1,jmt
c
c       set metric factors = 0 on beta and f plane
c       (strictly, they should be non zero... but they are small)
c
	am3_biha(j)   = 0.; am4_biha(j,1) = 0.; am4_biha(j,2) = 0.
	am3_ha(j)   = 0.; am4_ha(j,1) = 0.; am4_ha(j,2) = 0.
        advmet(j,1) = 0.; advmet(j,2) = 0.
       enddo
      else
       do j=1,jmt
	am3_biha(j)   = amix*(1.-tng(j)*tng(j))/(rho0*radius**2)
	am4_biha(j,1) = -amix*2.*sine(j)/(rho0*radius*csu(j)*csu(j))
	am4_biha(j,2) = -am4_biha(j,1)
	am3_ha(j)   = am*(1.-tng(j)*tng(j))/(rho0*radius**2)
	am4_ha(j,1) = -am*2.*sine(j)/(rho0*radius*csu(j)*csu(j))
	am4_ha(j,2) = -am4_ha(j,1)
        advmet(j,1) = tng(j)/radius
        advmet(j,2) = -advmet(j,1)
       enddo
      endif
c
c     for momentum
c
      visc_cnu = sqrt(abs(ambi)); visc_ceu = sqrt(abs(ambi))
c
      do j=1,jmt
	jm1 = max(1,j-1); jp1 = min(jmt,j+1)
	amc_north_ha(j) = am/rho0*cst(jp1)*dytr(jp1)*csur(j)*dyur(j)
	amc_south_ha(j) = am/rho0*cst(j)*dytr(j)*csur(j)*dyur(j)
      enddo
c
c     for tracers
c
      diff_cet = sqrt(abs(ahbi)); diff_cnt = sqrt(abs(ahbi))
c
      do j=1,jmt
          jm1 = max(1,j-1); jp1 = min(jmt,j+1)
          ahc_north_ha(j) = ah*csu(j)*dyur(j)*cstr(j)*dytr(j)
          ahc_south_ha(j) = ah*csu(jm1)*dyur(jm1)*cstr(j)*dytr(j)
      enddo

      if (enable_obc_south) then 
       amc_south_ha(2) = 0.; ahc_south_ha(2) = 0.
      endif
      if (enable_obc_north) then 
       amc_north_ha(jmt-2) = 0.; amc_north_ha(jmt-1) = 0.
c      ahc_north_ha(jmt-2) = 0.; ahc_north_ha(jmt-1) = 0.
      endif
c
c-----------------------------------------------------------------------
c     scaling of am/ah/ambi/ahbi
c-----------------------------------------------------------------------
c
      if (mother<0 .and. (enable_obc_south_sponge .or. 
     &                      enable_obc_south_sponge_harm)) then
c
c      define a frictional (bi)-harmonic sponge layer for obc south
c      set scaling for mother domain only
c
       if (my_pe==0) then
        print*,''
        if (enable_obc_south_sponge) then
         print*,' defining a frictional biharm. sponge layer for',
     &    ' southern open boundary'
        endif
        if (enable_obc_south_sponge_harm) then
         print*,' defining a frictional harm. sponge layer for',
     &    ' southern open boundary'
        endif
        if (enable_obc_south_sponge_diff) then
          print*,' scaling also (z-level) diffusivity'
        endif
        print*,' width in zonal slabs : ',obc_south_sponge_width
        print*,' e folding scale, same unit :',obc_south_sponge_scale
        print*,' factor for increased harm. frict. :',
     &              obc_south_sponge_fac 

        if (.not.enable_friction_biharmonic) then
         print*,' However, the biharm. sponge layer will not work'
         print*,' if you do not enable biharmonic friction '
        endif
        if (.not.enable_friction_harmonic) then
         print*,' However, the harm. sponge layer will not work'
         print*,' if you do not enable harmonic friction '
        endif
        print*,''
       endif

       do j=js_pe-2,je_pe+2
          fxa = 1.0
          if (j < obc_south_sponge_width) then   
            fxb   = -float(j-1)/obc_south_sponge_scale
            fxa = 1.0 + (obc_south_sponge_fac - 1.0) * exp(fxb) 
          end if
          if (enable_obc_south_sponge) then
            ambi_scale(:,j) = sqrt( abs(fxa) )
            if (enable_obc_south_sponge_diff) then
             ahbi_scale(:,j) = sqrt( abs(fxa) )
            endif
          endif
          if (enable_obc_south_sponge_harm) then
            am_scale(:,j) = abs(fxa) 
            if (enable_obc_south_sponge_diff) then
             ah_scale(:,j) = abs(fxa)
            endif
          endif
       end do

      endif

      if (mother<0 .and. (enable_obc_north_sponge.or.
     &                      enable_obc_north_sponge_harm) ) then
c
c      define a frictional (bi)harmonic sponge layer for obc north
c      set scaling for mother domain only
c
       if (my_pe==0) then
        print*,''
        if (enable_obc_north_sponge) then
         print*,' defining a frictional biharm. sponge layer for',
     &    ' northern open boundary'
        endif
        if (enable_obc_north_sponge_harm) then
         print*,' defining a frictional harm. sponge layer for',
     &    ' northern open boundary'
        endif
        if (enable_obc_north_sponge_diff) then
          print*,' scaling also (z-level) diffusivity'
        endif
        print*,' width in zonal slabs : ',obc_north_sponge_width
        print*,' e folding scale, same unit :',obc_north_sponge_scale
        print*,' factor for increased harm. frict. :',
     &              obc_north_sponge_fac 

        if (.not.enable_friction_biharmonic) then
         print*,' However, the biharm. sponge layer will not work'
         print*,' if you do not enable biharmonic friction '
        endif
        if (.not.enable_friction_harmonic) then
         print*,' However, the harm. sponge layer will not work'
         print*,' if you do not enable harmonic friction '
        endif
        print*,''
       endif

       do j=js_pe-2,je_pe+2
          fxa = 1.0
          if (j > jmt-obc_north_sponge_width) then   
            fxb   = -float(jmt-j)/obc_north_sponge_scale
            fxa = 1.0 + (obc_north_sponge_fac - 1.0) * exp(fxb) 
          end if
          if (enable_obc_north_sponge) then
           ambi_scale(:,j) = ambi_scale(:,j)*sqrt( abs(fxa) )
           if (enable_obc_north_sponge_diff) then
            ahbi_scale(:,j) = ahbi_scale(:,j)*sqrt( abs(fxa) )
           endif
          endif
          if (enable_obc_north_sponge_harm) then
           am_scale(:,j) = am_scale(:,j)*abs(fxa) 
           if (enable_obc_north_sponge_diff) then
            ah_scale(:,j) = ah_scale(:,j)*abs(fxa) 
           endif
          endif
       end do

      endif


      if (mother<0 .and. (enable_obc_east_sponge.or.
     &                      enable_obc_east_sponge_harm) ) then
c
c      define a frictional (bi)harmonic sponge layer for obc east
c      set scaling for mother domain only
c
       if (my_pe==0) then
        print*,''
        if (enable_obc_east_sponge) then
         print*,' defining a frictional biharm. sponge layer for',
     &    ' eastern open boundary'
        endif
        if (enable_obc_east_sponge_harm) then
         print*,' defining a frictional harm. sponge layer for',
     &    ' eastern open boundary'
        endif
        if (enable_obc_east_sponge_diff) then
          print*,' scaling also (z-level) diffusivity'
        endif
        print*,' width in meridional slabs : ',obc_east_sponge_width
        print*,' e folding scale, same unit :',obc_east_sponge_scale
        print*,' factor for increased harm. frict. :',
     &              obc_east_sponge_fac 

        if (.not.enable_friction_biharmonic) then
         print*,' However, the biharm. sponge layer will not work'
         print*,' if you do not enable biharmonic friction '
        endif
        if (.not.enable_friction_harmonic) then
         print*,' However, the harm. sponge layer will not work'
         print*,' if you do not enable harmonic friction '
        endif
        print*,''
       endif

       do i=is_pe-2,ie_pe+2
c
c       add here amplification factor 
c
          fxa = 0.0
          if (i > imt-obc_east_sponge_width) then   
            fxb   = -float(imt-i)/obc_east_sponge_scale
            fxa = obc_east_sponge_fac * exp(fxb) 
          end if
          if (enable_obc_east_sponge) then
           ambi_scale(i,:) = ambi_scale(i,:)+sqrt( abs(fxa) )
           if (enable_obc_east_sponge_diff) then
            ahbi_scale(i,:) = ahbi_scale(i,:)+sqrt( abs(fxa) )
           endif
          endif
          if (enable_obc_east_sponge_harm) then
           am_scale(i,:) = am_scale(i,:)+abs(fxa) 
           if (enable_obc_east_sponge_diff) then
            ah_scale(i,:) = ah_scale(i,:)+abs(fxa) 
           endif
          endif
       end do

      endif

      if (mother<0 .and. (enable_obc_west_sponge.or.
     &                      enable_obc_west_sponge_harm) ) then
c
c      define a frictional (bi)harmonic sponge layer for obc east
c      set scaling for mother domain only
c
       if (my_pe==0) then
        print*,''
        if (enable_obc_west_sponge) then
         print*,' defining a frictional biharm. sponge layer for',
     &    ' western open boundary'
        endif
        if (enable_obc_west_sponge_harm) then
         print*,' defining a frictional harm. sponge layer for',
     &    ' western open boundary'
        endif
        if (enable_obc_west_sponge_diff) then
          print*,' scaling also (z-level) diffusivity'
        endif
        print*,' width in meridional slabs : ',obc_west_sponge_width
        print*,' e folding scale, same unit :',obc_west_sponge_scale
        print*,' factor for increased harm. frict. :',
     &              obc_west_sponge_fac 

        if (.not.enable_friction_biharmonic) then
         print*,' However, the biharm. sponge layer will not work'
         print*,' if you do not enable biharmonic friction '
        endif
        if (.not.enable_friction_harmonic) then
         print*,' However, the harm. sponge layer will not work'
         print*,' if you do not enable harmonic friction '
        endif
        print*,''
       endif

       do i=is_pe-2,ie_pe+2
c
c       add here amplification factor 
c
          fxa = 0.0
          if (i < obc_west_sponge_width) then   
            fxb   = -float(i-1)/obc_west_sponge_scale
            fxa = obc_west_sponge_fac * exp(fxb) 
          end if
          if (enable_obc_west_sponge) then
           ambi_scale(i,:) = ambi_scale(i,:)+sqrt( abs(fxa) )
           if (enable_obc_west_sponge_diff) then
            ahbi_scale(i,:) = ahbi_scale(i,:)+sqrt( abs(fxa) )
           endif
          endif
          if (enable_obc_west_sponge_harm) then
           am_scale(i,:) = am_scale(i,:)+abs(fxa) 
           if (enable_obc_west_sponge_diff) then
            ah_scale(i,:) = ah_scale(i,:)+abs(fxa) 
           endif
          endif
       end do

      endif


      if (enable_friction_cosine_scal)  then

       if (my_pe==0) then
         print*,''
         print*,
     &    ' scaling hor. friction parameter with cosine of latitude'
         print*,''
       endif

       do j=js_pe-2,je_pe+2
        i=min(max(1,j),jmt)
        am_scale(:,j)  =am_scale(:,j)*csu(i)
        ambi_scale(:,j)=ambi_scale(:,j)*csu(i)**1.5
       enddo
      endif

      if (enable_diffusion_cosine_scal)  then

       if (my_pe==0) then
         print*,''
         print*,
     &    ' scaling hor. diffusion parameter with cosine of latitude'
         print*,''
       endif

       do j=js_pe-2,je_pe+2
        i=min(max(1,j),jmt)
        ah_scale(:,j)  =ah_scale(:,j)*cst(i)
        ahbi_scale(:,j)=ahbi_scale(:,j)*cst(i)**1.5
       enddo
      endif

      if (enable_no_isopyc_medwater) then

       if (my_pe==0) then
         print*,''
         print*,
     &    ' increasing hor. harm. diffusivity in medwater outflow-area'
         print*,''
       endif

       do j=js_pe-2,je_pe+2 
        do i=is_pe-2,ie_pe+2
c     position Gibraltar: for HS_4 i,j = 112,682 (folded)      
c        distance = dstnc(yt(j),xt(i),yt(682),xt(112))
         distance = dstnc(yt(j),xt(i),35.8585,269.2083)
         if (distance < 500.0) then
          if (distance <= 300.0) then
           ah_scale(i,j) = ah_scale(i,j)*14.
          else
           ah_scale(i,j) = ah_scale(i,j)*(0.07*(500.0-distance))
          endif
         endif 
        enddo
       enddo
      endif  
c
      end subroutine grids



      function dstnc(phi1,rla1,phi2,rla2)
c returns distance (in km) between two geographical
c points
c        written by   p.herrmann
c      phi1 : first latitude
c      rla1 : first longitude
c      phi2 : second latitude
c      rla2 : second longitude

       real p1,p2,rl1,rl2,x,xx
       data s/0.0174533/
       p1=phi1*s; p2=phi2*s; rl1=rla1*s; rl2=rla2*s
       x=sin(p1)*sin(p2)+cos(p1)*cos(p2)*cos(rl2-rl1)
       xx=abs(x)
       if(xx.gt.1.0) x=1.0
       d=atan(sqrt((1-x)/(1+x)))*222.24/s; dstnc=d
       end function dstnc



      subroutine read_grid ()
c
c-----------------------------------------------------------------------
c     read in the grid from a binary file
c-----------------------------------------------------------------------
c
      use spflame_module
      implicit none
      integer i,j,k,imt2,jmt2,km2,io

      if (my_pe==0) print*,
     & ' =>Reading the grid definition from file grid.dta'
      call getunit (io, 'grid.dta', 'u s r ieee')
      read (io) ! no checks for awi id
      read (io) imt2, jmt2, km2
      if (imt2/=imt .or. jmt2 /= jmt .or. km2 /= km) then
c      check imt2,jmt2,km2
       if (my_pe==0) then
        print*,' Error: '
        print*,' read imt,jmt,km=',imt2,jmt2,km2
        print*,' but input from namelist was ',imt,jmt,km
        call sub_flush(6)
       endif
       call halt_stop(' reading grid file ')
      endif
      read (io) (dxtdeg(i),i=1,imt)
     &,         (dytdeg(j),j=1,jmt)
     &,         (dxudeg(i),i=1,imt)
     &,         (dyudeg(j),j=1,jmt)
     &,         (dzt(k),k=1,km)
     &,         (dzw(k),k=0,km)
     &,         (xt(i),i=1,imt)
     &,         (xu(i),i=1,imt)
     &,         (yt(j),j=1,jmt)
     &,         (yu(j),j=1,jmt)
     &,         (zt(k),k=1,km)
     &,         (zw(k),k=1,km)
      close (io)
      if (my_pe==0) print*,' done'
      end subroutine read_grid




      subroutine topog 
c
c-----------------------------------------------------------------------
c     construct the topographic mask "kmt" which determines the geometry
c     and topography by defining the number of vertical levels on model
c     T points. Construct also other topography related variables
c     SPFLAME version : c.eden
c---------------------------------------------------------------------
c
      use spflame_module
      implicit none
      integer i,j,k,io,kk,jj,ii,is,ie,js,je,kb
      integer kmu_big(imt,jmt),indp
      real p_cell_min,hy,hcm
#ifdef partial_cell
      real dhue,dhun
      dhue(i,k,j)   = min(dhu(i+1,k,j),dhu(i,k,j))
      dhun(i,k,j)   = min(dhu(i,k,j+1),dhu(i,k,j))
#endif

      if (mother<0 ) then
c
c---------------------------------------------------------------------
c     read or construct the kmt-mask
c---------------------------------------------------------------------
c
       if (enable_simple_topo) then

         call kmt_template()

       else

        if (my_pe==0) print*,
     &   ' =>Reading a "kmt" field from binary file kmt.dta'
        call getunit (io, 'kmt.dta', 'usr ieee')
        read (io) ! no checks
c       skip header record and read dimensional info
        read (io) 
        read (io) i, j, k
        if (i/=imt .or. j /= jmt .or. k/= km) then
c        check i,j,k
        if (my_pe==0) then
         print*,' Error: '
         print*,' read imt,jmt,km=',i,j,k
         print*,' but input from namelist was ',imt,jmt,km
         call sub_flush(6)
        endif
        call halt_stop(' reading kmt file ')
       endif
       read (io) 
       read (io) kmt_big
       close(io)
       if (my_pe==0) print*,' done'
      endif
     
#ifdef partial_cell
c
c---------------------------------------------------------------------
c     read or contruct  htp = depth to ocean bottom at T-cells (cm)
c---------------------------------------------------------------------
c
      if (enable_simple_topo) then

        call htp_template()

      else

       if (my_pe==0) print*,
     & ' =>Reading a "htp" field from binary file htp.dta'
       call getunit (io, 'htp.dta', 'usr ieee')
       read (io) ! no checks
c       skip header record and read dimensional info
       read (io) 
       read (io) i, j, k
       if (i/=imt .or. j /= jmt .or. k/= km) then
c       check i,j,k
        if (my_pe==0) then
         print*,' Error: '
         print*,' read imt,jmt,km=',i,j,k
         print*,' but input from namelist was ',imt,jmt,km
         call sub_flush(6)
        endif
        call halt_stop(' reading htp file ')
       endif
       read (io) 
       read (io) htp
       close(io)
       if (my_pe==0) print*,' done'
      endif
#endif

      endif  ! mother<0
c
c---------------------------------------------------------------------
c     exchange kmt over sub domains
c---------------------------------------------------------------------
c
      call exchg_sub_topo   ! exchange also ht 
c
c---------------------------------------------------------------------
c     check kmt at the boundaries
c---------------------------------------------------------------------
c
      if (enable_obc_south) then
       j=2
       do i=1,imt
        k=min(min(kmt_big(i,j),kmt_big(i,j+1)),
     &        min(kmt_big(i,j),kmt_big(i,j+2)))
        if (k /= kmt_big(i,j)) then
         if (my_pe==0) then
          write(6,*)' WARNING: Changing Topography at southern OB. '
          print*,'  at i=',i,' j=',j
          print*,'  new kmt(i,j)=',k
         endif
         kmt_big(i,j)=k
         kmt_big(i,j-1)=kmt_big(i,j)
         kmt_big(i,j+1)=kmt_big(i,j)
         kmt_big(i,j+2)=kmt_big(i,j+1)
#ifdef partial_cell
         do jj=j-1,j+2; htp(i,jj) = zw(k); enddo
#endif
        endif
       enddo
      else
       kmt_big(:,1)=0.
#ifdef partial_cell
       htp(:,1)=0.
#endif
      endif

      if (enable_obc_north) then
       j=jmt-1
       do i=1,imt
        k=min(min(kmt_big(i,j),kmt_big(i,j-1)),
     &        min(kmt_big(i,j),kmt_big(i,j-2)))
        if (k /= kmt_big(i,j)) then
         if (my_pe==0) then
          write(6,*)' WARNING: Changing Topography at northern OB. '
          print*,'  at i=',i,' j=',j
          print*,'  new kmt(i,j)=',k
         endif
         kmt_big(i,j)=k
         kmt_big(i,j+1)=kmt_big(i,j)
         kmt_big(i,j-1)=kmt_big(i,j)
         kmt_big(i,j-2)=kmt_big(i,j-1)
#ifdef partial_cell
         do jj=j-2,j+1; htp(i,jj) = zw(k); enddo
#endif
        endif
       enddo
      else
       kmt_big(:,jmt)=0
#ifdef partial_cell
       htp(:,jmt)=0.
#endif
      endif

      if (enable_obc_west) then
       i=2
       do j=1,jmt
        k=min(min(kmt_big(i,j),kmt_big(i+1,j)),
     &        min(kmt_big(i,j),kmt_big(i+2,j)))
        if (k /= kmt_big(i,j)) then
         if (my_pe==0) then
          write(6,*)' WARNING: Changing Topography at western OB. '
          print*,'  at i=',i,' j=',j
          print*,'  new kmt(i,j)=',k
         endif
         kmt_big(i,j)=k
         kmt_big(i-1,j)=kmt_big(i,j)
         kmt_big(i+1,j)=kmt_big(i,j)
         kmt_big(i+2,j)=kmt_big(i+1,j)
#ifdef partial_cell
         do ii=i-1,i+2; htp(ii,j) = zw(k); enddo
#endif
        endif
       enddo
      elseif (cyclic) then
        kmt_big(1,:)=kmt_big(imt-1,:)
#ifdef partial_cell
        htp(1,:)=htp(imt-1,:)
#endif
      else
        kmt_big(1,:)=0
#ifdef partial_cell
        htp(1,:)=0.
#endif
      endif

      if (enable_obc_east) then
       i=imt-1
       do j=1,jmt
        k=min(min(kmt_big(i,j),kmt_big(i-1,j)),
     &        min(kmt_big(i,j),kmt_big(i-2,j)))
        if (k /= kmt_big(i,j)) then
         if (my_pe==0) then
          write(6,*)' WARNING: Changing Topography at eastern OB. '
          print*,'  at i=',i,' j=',j
          print*,'  new kmt(i,j)=',k
         endif
         kmt_big(i,j)=k
         kmt_big(i+1,j)=kmt_big(i,j)
         kmt_big(i-1,j)=kmt_big(i,j)
         kmt_big(i-2,j)=kmt_big(i-1,j)
#ifdef partial_cell
         do ii=i-2,i+1; htp(ii,j) = zw(k); enddo
#endif
        endif
       enddo
      elseif (cyclic) then
       kmt_big(imt,:)=kmt_big(2,:)
#ifdef partial_cell
       htp(imt,:)=htp(imt,:)
#endif
      else
       kmt_big(imt,:)=0
#ifdef partial_cell
       htp(imt,:)=0.
#endif
      endif

      if ((enable_obc_north .or. enable_obc_south .or.
     &  enable_obc_west .or. enable_obc_east) .and. .not.cyclic ) then 
c      last checks in case of o.b.c
       kmt_big(imt,1)=0
       kmt_big(imt,jmt)=0
       kmt_big(1,jmt)=0
       kmt_big(1,1)=0
#ifdef partial_cell
       htp(imt,1)=0.
       htp(imt,jmt)=0.
       htp(1,jmt)=0.
       htp(1,1)=0.
#endif
c      done
      endif
c
      kmt=kmt_big(is_pe:ie_pe,js_pe:je_pe)
c
c---------------------------------------------------------------------
c     construct topography for u grid
c---------------------------------------------------------------------
c
      kmu_big=0
      do j=1,jmt-1
        do i=1,imt-1
           kmu_big(i,j) = min (kmt_big(i,j), kmt_big(i+1,j), 
     &                         kmt_big(i,j+1), kmt_big(i+1,j+1))
        enddo
      enddo

      if (enable_obc_north) kmu_big(:,jmt) = kmu_big(:,jmt-1)
      if (enable_obc_south) kmu_big(:,1) = kmu_big(:,2)
      if (enable_obc_west) kmu_big(1,:) = kmu_big(2,:)
      if (enable_obc_east) kmu_big(imt,:) = kmu_big(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_big(imt-1:imt,1)=0.
       kmu_big(1,1)=0.
       kmu_big(1,jmt-1:jmt)=0.
c      done
      endif

      if (cyclic) kmu_big(imt,:) = kmu_big(2,:)

      kmu=kmu_big(is_pe:ie_pe,js_pe:je_pe)

#ifdef partial_cell
c
c---------------------------------------------------------------------
c     a check
c---------------------------------------------------------------------
c
      if (my_pe==0) print*,' checking consitency of htp/kmt masks'
      do j=1,jmt
       do i=1,imt
        if (kmt_big(i,j)/=0) then
         if (htp(i,j)>zw(max(1,kmt_big(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_big(i,j) ))
          stop
         endif
         if (htp(i,j)<=zw(max(1, kmt_big(i,j)-1 ))) 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_big(i,j)-1 ))
          print*,' kmt(i,j)=',kmt_big(i,j)
          stop
         endif
        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(10.0e2,zw(1)) ! schiet wat
      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_big(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
              if (my_pe==0) then
               print*,' WARNING: partial cell limit violated'
               print*,' at i=',i,' j=',j,' htp(i,j)=',htp(i,j)
               print*,' zw(k-1)=',zw(k-1)
               print*,' zw(k)=',zw(k),' p_cell_min',p_cell_min
              endif
	      htp(i,j) = zw(k-1) + p_cell_min
	    endif
	  endif
        enddo
      enddo 
      if (my_pe==0) print*,' done'
c
c-----------------------------------------------------------------------
c     construct hu as the depth to bottom on U cells
c-----------------------------------------------------------------------
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,:)
c
c     fracdz(k,0) = fractional distance between grid point and cell top
c     fracdz(k,1) = fractional distance between grid point and cell bot
c
      do k=1,km
        if (k .gt. 1) then
          fracdz(k,0) = (zt(k) - zw(k-1))/dzt(k)
        else
          fracdz(k,0) = (zt(k) - 0.)/dzt(k)
        endif
        fracdz(k,1) = (zw(k) - zt(k))/dzt(k)
      enddo
c
c-----------------------------------------------------------------------
c     construct thickness of T cells and U cells
c     (an aberration of the canonical indexing is used here)
c-----------------------------------------------------------------------
c
      do k=1,km
       dht(:,k,:)=dzt(k)
       dhu(:,k,:)=dzt(k)
       ztp(:,k,:)=zt(k)
      enddo
      do k=0,km
       dhwt(:,k,:)=dzw(k)
      enddo

      do j=js_pe,je_pe
	do i=is_pe,ie_pe
	  kb = kmt_big(i,j)
	  if (kb .gt. 1) then
	    dht(i,kb,j) = htp(i,j) - zw(kb-1)
	  endif
	enddo
      enddo
      call border_exchg(dht,km,2)
      if (cyclic) call set_cyclic(dht,km,2)
c
      do j=js_pe,je_pe
	do i=is_pe,ie_pe
	  kb = kmt_big(i,j)
	  if (kb .gt. 1) then
	    ztp(i,kb,j) = zw(kb-1) + fracdz(kb,0)*dht(i,kb,j)
            dhwt(i,kb-1,j) = ztp(i,kb,j) - zt(kb-1)
	  endif
	enddo
      enddo

      call border_exchg(dhwt,km+1,2)
      call border_exchg(ztp,km,1)
      if (cyclic) call set_cyclic(dhwt,km+1,2)
      if (cyclic) call set_cyclic(ztp,km,1)

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      do j=js,je
        do k=1,km
          do i=is,ie
	    dhu(i,k,j) = min(dht(i,k,j), dht(i+1,k,j), dht(i,k,j+1)
     &,                      dht(i+1,k,j+1))
	  enddo
          if (my_blk_i==1)       dhu(1,k,j)   = dht(1,k,j)
          if (my_blk_i==n_pes_i) dhu(imt,k,j) = dht(imt,k,j)
	enddo
      enddo

      call border_exchg(dhu,km,2)
      if (cyclic) call set_cyclic(dhu,km,2)
      if (my_blk_j==1) dhu(:,:,1)=dhu(:,:,2)
      if (my_blk_j==n_pes_j) dhu(:,:,jmt)=dhu(:,:,jmt-1)
c
c       compute momentum sink by partial-cells
c
c from here:   dhu(is_pe-2:ie_pe+2,km,js_pe-2:je_pe+2)
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      if (enable_friction_biharmonic) then
        do j=js,je
           do k=1,km
              do i=is,ie
                pc_sink_biha(i,k,j) =
     &           - visc_ceu*         
     &          (csur(j)*dxur(i)/dhu(i,k,j))*
     &          (csur(j)*dxtr(i+1)*(dhu(i,k,j)-dhue(i,k,j))
     &         + csur(j)*dxtr(i)*(dhu(i,k,j)-dhue(i-1,k,j)))
     &         - visc_cnu*
     &          (csudyur(j)/dhu(i,k,j))*
     &          (cst_dytr(j+1)*(dhu(i,k,j)-dhun(i,k,j))
     &         + cst_dytr(j)*(dhu(i,k,j)-dhun(i,k,j-1)))
              enddo
           enddo
        enddo
        call border_exchg(pc_sink_biha,km,1)
        if (cyclic) call set_cyclic(pc_sink_biha,km,1)
      endif

      if (enable_friction_harmonic) then
        do j=js,je
         do k=1,km
              do i=is,ie
                pc_sink_ha(i,k,j) =
     &           - am*         
     &          (csur(j)*dxur(i)/dhu(i,k,j))*
     &          (csur(j)*dxtr(i+1)*(dhu(i,k,j)-dhue(i,k,j))
     &         + csur(j)*dxtr(i)*(dhu(i,k,j)-dhue(i-1,k,j)))
     &         - am*
     &          (csudyur(j)/dhu(i,k,j))*
     &          (cst_dytr(j+1)*(dhu(i,k,j)-dhun(i,k,j))
     &         + cst_dytr(j)*(dhu(i,k,j)-dhun(i,k,j-1)))
              enddo
         enddo
        enddo
      endif
#endif

c---------------------------------------------------------------------
c     compute depths and reciprocal depths over "u" cells
c---------------------------------------------------------------------

      hr=0.
      do j=max(1,js_pe-1),min(je_pe+1,jmt)
        do i=max(is_pe-1,1),min(ie_pe+1,imt)
          if (kmu_big(i,j) .ne. 0) then
#ifdef partial_cell
            hr(i,j) = 1./hup(i,j)
#else
            hr(i,j) = 1./zw(kmu_big(i,j))
#endif
          endif
        enddo
      enddo
c
c---------------------------------------------------------------------
c     contruct masks for tracer and velocities
c---------------------------------------------------------------------
c
      umask=0.; tmask=0.
      if ( mother<0 ) then
        kk = 2
      else 
        kk = max(2,zoom_fac/2)
      endif
      do j=max(1,js_pe-kk),min(je_pe+kk,jmt)
        do k=1,km
           do i=max(is_pe-kk,1),min(ie_pe+kk,imt)
              if (kmt_big(i,j) .ge. k) then
                tmask(i,k,j) = 1.
              else
                tmask(i,k,j) = 0.
              endif
              if (kmu_big(i,j) .ge. k) then
                umask(i,k,j) = 1.
              else
                umask(i,k,j) = 0.
              endif
           enddo
        enddo
      enddo
c
c     we are using partially a 4th order scheme
c
      if (enable_quicker_advection) 
     &      call set_cyclic_4th_order(tmask,km,kk)

c-----------------------------------------------------------------------
c     compute surface area and volume of ocean ("t" cells and "u" cells)
c     (note that areas are defined at each level)
c-----------------------------------------------------------------------
c
      tcella=0.; ucella=0.; tcellv = 0.; ucellv = 0.
      do j=js_pe,je_pe
        do i=is_pe,ie_pe
          if (kmt(i,j) .gt. 0) then
	    do k=1,kmt(i,j)
              tcella(k) = tcella(k) + cst(j)*dxt(i)*dyt(j)
	    enddo
            tcellv = tcellv +cst(j)*dxt(i)*dyt(j)
#ifdef partial_cell
     &              *htp(i,j)
#else
     &              *zw(kmt(i,j))
#endif
          endif
          if (kmu(i,j) .gt. 0) then
	    do k=1,kmu(i,j)
              ucella(k) = ucella(k) + csu(j)*dxu(i)*dyu(j)
	    enddo
            ucellv = ucellv +csu(j)*dxu(i)*dyu(j)
#ifdef partial_cell
     &              *hup(i,j)
#else
     &              *zw(kmu(i,j))
#endif
          endif
        enddo
      enddo
      call global_sum(tcellv)
      call global_sum(ucellv)
      do k=1,km
       call global_sum(tcella(k))
       call global_sum(ucella(k))
      enddo
      if (my_pe==0) write (6,9341) tcella(1), tcellv, ucella(1), ucellv
9341  format (//,'  Global ocean statistics:'
     &/,'  the total ocean surface area (t cells) =',1pe15.8,'cm**2'
     &/,'  the total ocean volume (t cells)       =',1pe15.8,'cm**3'
     &/,'  the total ocean surface area (u cells) =',1pe15.8,'cm**2'
     &/,'  the total ocean volume (u cells)       =',1pe15.8,'cm**3')
      end subroutine topog


