#include "options.inc"


      subroutine prep_grid
      use prep_module
      implicit none

      integer,parameter :: maxlen = 2000
      real :: dxtdeg_t(maxlen),dxudeg_t(maxlen)
      real :: dytdeg_t(maxlen),dyudeg_t(maxlen)
      real :: dzt_t(maxlen),dzw_t(0:maxlen)
      real :: yt_t(maxlen),yu_t(0:maxlen)
c
c     FLAME vertical grid
      integer, parameter :: km_flame = 45
      real dz_flame(km_flame)
      data dz_flame/ 
     &  1000.0,   1044.051, 1102.236, 1169.269, 1246.635, 1336.169,
     &  1440.147, 1561.417, 1703.569, 1871.169, 2070.068, 2307.849,
     &  2594.425, 2992.892, 3518.793, 4100.175, 4863.114, 5680.496,
     &  6513.772, 7038.668, 7397.827, 7758.212, 8130.866, 8726.499,
     &  9576.454,10908.060,12732.562,15187.976,18426.63,
     &  22000.,  24182.902,24817.098, 13*250.00e2 /
  
c   POP 1/10 deg vertical grid
      integer, parameter :: km_pop = 40
      real dz_pop(km_pop)
      data dz_pop/ 
     &   10.01244e+2, 10.11258e+2, 10.31682e+2, 10.6333e+2, 11.0751e+2, 
     &   11.6615e+2, 12.4193e+2, 13.3861e+2, 14.6140e+2, 16.1756e+2,
     &   18.1737e+2, 20.7556e+2, 24.1368e+2, 28.6382e+2, 34.7464e+2,
     &   43.2086e+2, 55.1681e+2, 72.3046e+2, 96.7490e+2, 130.039e+2,
     &   170.049e+2, 207.993e+2, 233.569e+2, 245.272e+2, 248.980e+2,
     &   249.832e+2, 249.979e+2, 249.998e+2, 250.000e+2, 250.000e+2,
     &   250.000e+2, 250.000e+2, 250.000e+2, 250.000e+2, 250.000e+2,
     &   250.000e+2, 250.000e+2, 250.000e+2, 250.000e+2, 250.000e+2  /

      integer i,j,n,m,k
      integer :: nbpts =2
      real :: offset=0.
      real :: stretch_x=1.0, stretch_y=1.0,stretch_z=1.0
      real :: num=0.,tolr=1.e-5,dxubar,dyubar,dzwbar,zz
      real sumt,sumu

      print*,''
      print*,' Setting up the grid'
      print*,''
      if (centered_t) then 
         print*,' T grid points are centered in their boxes'
         print*,' while U grid points are uncentered'
         print*,' Note that this is non-standard in MOM2'
      else
         print*,' U grid points are centered in their boxes'
         print*,' while T grid points are uncentered'
         print*,' This is standard setting in MOM2'
      endif

      if (read_grid_topo) then

       call  read_grid_bin()

      else
c
c-----------------------------------------------------------------------
c     Calculate resolution in longitude. Add one boundary cell at the
c     start and end of the domain so that calculations are meaningful
c     for grid cells i=2,imt-1
c-----------------------------------------------------------------------
c
      nbpts = 2
      write (6,'(/a,i1,a)')
     & ' Generating the longitudinal resolution ( ', nbpts
     &, ' extra boundary points will be used).'
c
c     convert negative longitudes (degrees west of Greenwich) to
c     positive longitudes
c
      do n=1,nxlon
        if (x_lon(n) < 0.0) then
          offset = 360.0
        endif
      enddo
      if (offset > 0.0) then
          write (6,'(/,a/,a/,a)')
     &    '=>Warning: adding 360.0 degrees to all longitudes to convert'
     &,   '           negative values (degrees west of Greenwich) to'
     &,   '           positive values (degrees east of Greenwich)'
      endif
      x_lon = x_lon + offset
c
c     if a region contains the greenwich meridian, compensate by adding
c     360.0 degrees to all remaining regions
c
      do n=2,nxlon
        if (x_lon(n-1) >  x_lon(n)) then
          x_lon(n) = x_lon(n) + 360.0
          write (6,'(/,a,i1,a,a/)')
     &    ' Warning: adding 360.0 degrees to x_lon(',n,')'
     &,   ' to insure region boundaries increase monotonically'
        endif
      enddo
c
c     if width of domain exceeds 360.0 deg... limit width to 360.0
c
      if (x_lon(nxlon) - x_lon(1) .gt. 360.0) then
        write (6,'(/a,a,g14.7/)')
     &  '=>Warning: Domain width exceeds 360 deg. Restricting last'
     &, ' x_lon to ',x_lon(1)+360.0
          do n=1,nxlon
            write (6,'(i3,f10.5)') n, x_lon(n)
          end do
      endif
      if (x_lon(nxlon) - x_lon(1) .eq. 360.0) then
        if (dx_lon(nxlon) .ne. dx_lon(1)) then
          write (6,'(/a,a)')
     &    '=>Error: dx_lon(1) must equal dx_lon(last) when domain'
     &,   ' width = 360.0 degrees'
          stop
        endif
      endif

      if (centered_t) then
       call gcell_centered_t (maxlen, nxlon, x_lon,  dx_lon, nbpts,
     &           imt, dxtdeg_t, dxudeg_t, stretch_x,0,tol_gcell )
      else
       call gcell (maxlen, nxlon, x_lon,  dx_lon, nbpts,
     &           imt, dxtdeg_t, dxudeg_t, stretch_x,0,tol_gcell )
      endif

      allocate( dxtdeg(imt), dxudeg(imt), xt(imt), xu(imt))
      dxtdeg=dxtdeg_t(1:imt)
      dxudeg=dxudeg_t(1:imt)
c
c-----------------------------------------------------------------------
c     Build the longitudinal grid points for a "B" grid
c     (account for an extra boundary point at the start)
c-----------------------------------------------------------------------
c
      xt(1) = x_lon(1) - 0.5*dx_lon(1)
      xu(1) = x_lon(1)
c
      do i=2,imt
        xu(i) = xu(i-1) +  dxtdeg(i)
        xt(i) = xt(i-1) +  dxudeg(i-1)
      enddo
c
c-----------------------------------------------------------------------
c     Calculate resolution in latitude. Add one boundary cell at the
c     start and end of the domain so that calculations are meaningful
c     for grid cells jrow=2,jmt-1
c-----------------------------------------------------------------------
c
      nbpts = 2
      write (6,'(/a,i1,a)')
     & ' Generating the latitudinal resolution ( ', nbpts
     &, ' extra boundary points will be used).'
c
c     if width of domain exceeds 180.0 deg... limit width to 180.0
c
      if (y_lat(nylat) - y_lat(1) .gt. 180.0) then
        write (6,'(/a,a,g14.7/)')
     &  '=>Warning: Latitudinal domain width exceeds 180 deg.'
     &, ' Restricting last y_lat to ',y_lat(1)+180.0
      endif

      if (isotropic) then

        print*,' constructing a Mercator grid '
c
c       this subroutine from MOM 3 seem to be buggy
c
c       call iso_grid (maxlen, nylat, y_lat, dy_lat, jmt, 
c     &                dytdeg_t, dyudeg_t, yt_t, yu_t)
c       allocate( dytdeg(jmt), dyudeg(jmt), yt(jmt), yu(jmt))
c       dytdeg=dytdeg_t(1:jmt)
c       dyudeg=dyudeg_t(1:jmt)
c       yt=yt_t(1:jmt)
c       yu=yu_t(1:jmt)
c      bugfix
c       if (yu(jmt) == 0) yu(jmt)=yu(jmt-1)+dyudeg(jmt)

c
c      use the old MOM2 routine here
c
       if (centered_t) then
        call isotro_centered_t(maxlen, nylat, y_lat, dy_lat, 
     &                 jmt, dytdeg_t, dyudeg_t)
       else
        call sub_isotro(maxlen, nylat, y_lat, dy_lat, 
     &                  jmt, dytdeg_t, dyudeg_t)
       endif
       allocate( dytdeg(jmt), dyudeg(jmt), yt(jmt), yu(jmt))
       dytdeg=dytdeg_t(1:jmt); dyudeg=dyudeg_t(1:jmt)
       yu(1) = y_lat(1) + dytdeg(1)
       yt(1) = y_lat(1) + dy_lat(1)
       do j=2,jmt
        yu(j) = yu(j-1) +  dytdeg(j)
        yt(j) = yt(j-1) +  dyudeg(j-1)
        if (abs(yu(j)).lt.1.0e-10)  yu(j) = 0.0
       enddo

      else

       do n=2,nylat
        if (y_lat(n-1) .gt. y_lat(n)) then
          write (6,'(/,a,/a/)')
     &    ' =>Error: latitude boundaries must increase monotonically'
     &,   '          check the specifications in the USER INPUT section'
          do m=1,nylat
            write (6,'(i3,f10.5)') m, y_lat(m)
          end do
          stop
        endif
       enddo

      if (centered_t) then
       call gcell_centered_t (maxlen, nylat, y_lat, dy_lat, nbpts,
     &            jmt, dytdeg_t, dyudeg_t, stretch_y,0,tol_gcell)
      else
       call gcell (maxlen, nylat, y_lat, dy_lat, nbpts,
     &            jmt, dytdeg_t, dyudeg_t, stretch_y,0,tol_gcell)
       endif

       allocate( dytdeg(jmt), dyudeg(jmt), yt(jmt), yu(jmt))
       dytdeg=dytdeg_t(1:jmt)
       dyudeg=dyudeg_t(1:jmt)
c
c-----------------------------------------------------------------------
c     Build the latitudinal grid points on a "B" grid
c     (account for an extra boundary point at the start)
c-----------------------------------------------------------------------
c
       yt(1) = y_lat(1) - 0.5*dy_lat(1)
       yu(1) = y_lat(1)
       do j=2,jmt
        yu(j) = yu(j-1) +  dytdeg(j)
        yt(j) = yt(j-1) +  dyudeg(j-1)
       enddo

      endif 

c
c-----------------------------------------------------------------------
c     Calculate resolution in depth. No boundary cells added here so
c     calculations are meaningful for k=1,km. Allow the bottom region 
c     to be stretched further if desired.
c-----------------------------------------------------------------------
c
      if (flame_vert_grid) then
       km=km_flame
       allocate( dzt(km), dzw(0:km), zt(km), zw(km))
       dzt(1:km)=dz_flame
       dzw(1) = dzt(1)
       sumt = dzt(1)
       sumu = 0.5*dzt(1) + dzw(1)
       do i=2,km
        sumt = sumt + dzt(i)
        dzw(i) = 2*(sumt-sumu)
        sumu = sumu + dzw(i)
       enddo
       zt(1) = 0.5*dzt(1)
       zw(1) = dzt(1)
       do i=2,km
        zw(i) = zw(i-1) + dzt(i)
        zt(i) = zt(i-1) + dzw(i-1)
       enddo
       dzw(0)  = zt(1)
       dzw(km) = zw(km) - zt(km)
      elseif (pop_vert_grid) then
       km=km_pop ! POP is t centered 
       allocate( dzt(km), dzw(0:km), zt(km), zw(km))
       do k=1,km
        dzt(k)=dz_pop(k)
       enddo
       sumt = 0.0
       do k=1,km-1
        dzw(k) = 0.5*(dzt(k+1) + dzt(k))
        sumt = sumt + dzt(k)
       enddo
       dzw(km) = .5*(dzt(km))
       sumt = sumt + dzt(km)
       zt(1) = .5*dzt(1)
       zw(1) = dzt(1)
       do k=2,km
        zw(k) = zw(k-1) + dzt(k)
        zt(k) = zt(k-1) + dzw(k-1)
       enddo
       dzw(0)  = zt(1)
       dzw(km) = zw(km) - zt(km)
      else
       nbpts = 0
       write (6,'(/a,i1,a)')
     & ' Generating the vertical resolution ( ', nbpts
     &, ' extra boundary points will be used).'
       do n=2,nzdepth
        if (z_depth(n-1) .gt. z_depth(n)) then
          write (6,'(/,a,/a/)')
     &    ' =>Error: depth boundaries must increase monotonically'
     &,   '          check the specifications in the USER INPUT section'
          stop
        endif
       enddo

      if (centered_t) then
       call gcell_centered_t(maxlen, nzdepth, z_depth, dz_depth, nbpts,
     &            km, dzt_t, dzw_t(1), stretch_z,0,tol_gcell)
      else
       call gcell(maxlen, nzdepth, z_depth, dz_depth, nbpts,
     &            km, dzt_t, dzw_t(1), stretch_z,0,tol_gcell)
      endif

       allocate( dzt(km), dzw(0:km), zt(km), zw(km))
       dzt=dzt_t(1:km)
       dzw=dzw_t(0:km)
c
c-----------------------------------------------------------------------
c     Build the vertical grid points on a "B" grid. The T and U
c     cells are staggered in the horizontal but at the same level in
c     the vertical. However, the W cells here refer to the vertical
c     advection velocities at the bottoms of the U and T cells.
c     (no extra boundary point at the start)
c-----------------------------------------------------------------------
c
       zt(1) = z_depth(1) + 0.5*dz_depth(1)
       zw(1) = z_depth(1) + dzt(1)
       do k=2,km
        zw(k) = zw(k-1) + dzt(k)
        zt(k) = zt(k-1) + dzw(k-1)
       enddo
c
c     set W cell thickness at surface and bottom to ocean part of cell
c
       dzw(0)  = zt(1)
       dzw(km) = zw(km) - zt(km)
      endif
c
c-----------------------------------------------------------------------
c     Check if the T grid resolution is an average of the
c     U cell resolution. This insures more accurate advection of
c     tracers within a stretched grid.
c-----------------------------------------------------------------------
c
      tolr=1.e-5
      write (6,'(/)')
      do i=2,imt-1
        dxubar = 0.5*(dxudeg(i) + dxudeg(i-1))
        if (abs(dxubar-dxtdeg(i)) .gt. tolr) then
          num = num + 1
          write (6,'(a,i5,a)')
     &    '=>Warning: T cell delta x at i=',i
     &,   ' is not an average of adjacent U cell delta x`s'     
        endif     
      enddo
c
      do j=2,jmt-1
        dyubar = .5*(dyudeg(j) + dyudeg(j-1))
        if (abs(dyubar-dytdeg(j)) .gt. tolr) then
          num = num + 1
          write (6,'(a,i5,a)')
     &    '=>Warning: T cell delta y at jrow=',j
     &,   ' is not an average of adjacent U cell delta y`s' 
        endif     
      enddo
c
      tolr = 1.e0
      do k=2,km-1
        dzwbar = 0.5*(dzw(k) + dzw(k-1))
        if (abs(dzwbar-dzt(k)) .gt. tolr) then
          num = num + 1
          write (6,'(a,i5,a)')
     &    '=>Warning: T cell delta z at k=',k
     &,   ' is not an average of adjacent W cell delta z`s'     
        endif     
      enddo

      endif ! read_grid_topo 
c
c-----------------------------------------------------------------------
c     write the old mom output file
c-----------------------------------------------------------------------
c
      call write_grid_bin()
c
c-----------------------------------------------------------------------
c     Print all grid coordinates
c-----------------------------------------------------------------------
c
      zz = 0.0
      write (6
     &,'(//,40x,a,//,a,g14.7,a,/a/,a,g14.7,a/a,/,a,g14.7,a/)') 
     &  ' Grid Point Coordinate details: '
     &, ' The western edge of the T-cell at i=2 is at:'
     &,  xu(1),' (deg) longitude',' (the 1st T-cell is a boundary cell)'
     &, ' The southern edge of the T-cell at jrow=2 is at:'
     &,  yu(1),' (deg) latitude',' (the 1st T-cell is a boundary cell)'
     &,' The top edge of the T-cell at k=1 is at z =',zz,' (cm)'
      write (6,'(a,i4,a,g14.7,a/a,i4,a)') 
     & ' The eastern edge of the T-cell at i=',imt-1,' is at:'
     &,  xu(imt-1),' (deg) longitude'
     &,' (the last T-cell (imt=',imt,') is a boundary cell)'
     &, ' The northern edge of the T-cell at jrow=',jmt-1,' is at:'
     &,  yu(jmt-1),' (deg) latitude'
     &,' (the last T-cell (jmt=',jmt,') is a boundary cell)'
     &,' The bottom edge of T-cell at k=',km,' is at z ='
     &,  zw(km),' (cm)'
      write (6,'(/,a,g14.7,a/a/,a,g14.7,a/a/,a,g14.7,a/)') 
     &  ' The western edge of the U-cell at i=2 is at:', xt(1)
     &, ' (deg) longitude',' (the 1st U-cell is a boundary point)'
     &, ' The southern edge of the U cell at jrow=1 is at:', yt(1)
     &, ' (deg) latitude',' (the 1st U cell is a boundary point)'
     &, ' The top edge of the W-cell at k=1 is at z =',zt(1),' (cm)'
      write (6,'(a,i4,a,g14.7,a/a,i4,a)') 
     & ' The eastern edge of the U-cell at i=',imt-1,' is at:'
     &,  xt(imt),' (deg) longitude'
     &,' (the last U-cell (imt=',imt,') is a boundary cell)'
     &, ' The northern edge of the U-cell at jrow=',jmt-1,' is at:'
     &,  yt(jmt),' (deg) latitude'
     &,' (the last U-cell (jmt=',jmt,') is a boundary cell)'
     &,' The bottom edge of W-cell at k=',km,' is at z ='
     &,  zw(km),' (cm)'
      write (6,9103) km
      write (6,9002) (zt(k),k=1,km)
      write (6     ,9104) km
      write (6     ,9002) (zw(k),k=1,km)
      write (6     ,9105) jmt
      write (6     ,9001) (yt(j),j=1,jmt)
      write (6     ,9106) jmt
      write (6     ,9001) (yu(j),j=1,jmt)
      write (6     ,9107) imt
      write (6     ,9001) (xt(i),i=1,imt)
      write (6     ,9108) imt
      write (6     ,9001) (xu(i),i=1,imt)

      print*,' '
      print*,' grid generation done ...'
      print*,' '

9001  format (1x,10f10.4)
9002  format (1x,10f10.2)
9101  format (/,  a,i4,' in units of ',a,' as follows:')
9103  format (/,' Depth to T & U grid points (cm): zt(k) k=1,',i3)
9104  format (/,' Depth to W grid points (cm): zw(k) k=1,',i3)
9105  format (/,' Latitude of T points (deg): yt(j) j=1,',i4)
9106  format (/,' Latitude of U points (deg): yu(j) j=1,',i4)
9107  format (/,' Longitude of T points (deg): xt(i) i=1,',i4)
9108  format (/,' Longitude of U points (deg): xu(i) i=1,',i4)


      end 
      




      subroutine gcell (maxlen, n_bounds, bounds, d_bounds, nbpts
     &,                 num, deltat, deltau, stretch,my_pe,tol_in)
c
c=======================================================================
c
c              G R I D   C E L L   C O N S T R U C T I O N
c
c     A domain is composed of one or more regions:
c     Build "num" T  cells with resolution "deltat(n) n=1,num" 
c     within the domain composed of regions bounded by "bounds".
c     Also construct "num" U  cells of resolution "deltau(n) n=1,num"
c     with the relation between T and U cells given by:
c     deltat(n) = 0.5*(deltau(n-1) + deltau(n))
c     Resolution may be constant or smoothly varying within each
c     region AND there must be an integral number of grid cells within 
c     each region. The domain is the sum of all regions.
c
c     inputs:
c
c     maxlen   = maximum length of "deltat" and "deltau"
c     n_bounds = number of bounds needed to define the regions
c     bounds   = latitude, longitude, or depth at each bound
c     d_bounds = delta (resolution) at each of the "bounds"
c     nbpts    = number of extra boundary cells to add to the domain.
c                (usually one at the beginning and end)
c     stretch  = stretching factor for last region (should only be used
c                in the vertical to provide increased stretching of grid
c                points. "stretch" = 1.0 gives no increased stretching.
c                "stretch" = 1.2 gives increased stretching...etc
c             
c     outputs:
c
c     num    = total number of grid cells within the domain
c     deltau = resolution of U grid cells: n=1,num
c     deltat = resolution of T grid cells: n=1,num
c
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
c=======================================================================
c
      implicit none
      integer stdout,maxlen,n_bounds,num,nbpts,my_pe
      parameter( stdout=6)
      real deltat(maxlen), deltau(maxlen)
      real d_bounds(n_bounds), bounds(n_bounds), stretch,tol_in

      logical keep_going
      double precision p5,pi ,avg_res, chg_res, tol, wid, an
      double precision del,d_new,sum
      integer l,n,i,ierr,m
c
c     Set some constants
c
      p5 = 0.5d0
      pi = 4.0d0*atan(1.0d0)
c
c     Do all regions, one at a time, to construct the domain
c
      num  = 1
      do l=1,n_bounds-1
c
       if (my_pe == 0) 
     &  write (stdout,'(2x,a,i2,a,g14.7,a,g14.7,a,g14.7,a,g14.7,a)')
     & ' region # ',l,'  going from ',bounds(l),' (res=',d_bounds(l)
     &,') to ',  bounds(l+1),' (res=',d_bounds(l+1),')' 
c
c       avg_res = average resolution of T cells within region
c       chg_res = change in resolution across the region
c       wid     = width of region
c       tol     = tolerance for fitting T cels within region width
c
c       provide for stretching last region if needed
c
        if (l .eq. n_bounds-1) then
          avg_res = p5*(d_bounds(l) + stretch*d_bounds(l+1))
          chg_res = (stretch*d_bounds(l+1) - d_bounds(l))
        else
          avg_res = p5*(d_bounds(l) + d_bounds(l+1))
          chg_res = (d_bounds(l+1) - d_bounds(l))
        endif
c
c
c       had to change this ...
c        tol = 1.d-3
c      all the time
        tol=tol_in
c
        wid = abs(bounds(l+1) - bounds(l))
        an  = wid/avg_res
        m   = nint(an)
c
c       Calculate resolution of U cells: "deltau"
c       U grid points will be centered in these cells
c       n = number of T cells fitting within the region boundaries
c       note: "sum" initially discounts half of the U cells widths
c       at the boundaries
c
        sum = p5*d_bounds(l) - p5*d_bounds(l+1)
        n   = 0
        i = 0
        keep_going = .true.
        do while (i .le. 100000 .and. keep_going)
          i = i + 1
          del = avg_res - p5*chg_res*cos((pi/m)*i)
          if (sum + del .le. wid*(1.0d0 + tol)) then
            sum = sum + del
            if (num+i-1 .gt. maxlen) then
              write (stdout,*) "=>Error: maxlen exceeded in gcell. "
     &,                        " ...increase size of maxlen"
             stop
            endif
            deltau(num+i-1) = del
            n = n + 1
          else
            keep_going = .false.
          endif
        enddo 

        if (l .eq. n_bounds-1 .and. stretch .ne. 1.0) then
       if (my_pe == 0) 
     &    write (stdout,'(a,i3,a,f5.2)')
     & '    constructed ',n,' cells with a stretch factor of ', stretch
       if (my_pe == 0) 
     &    write (stdout,'(/2(a,g14.7),/2(a,g14.7),/a,a/)')
     &    'Note: you specified the ocean bottom at ',bounds(l+1)
     &,   ' cm with a bottom cell thickness of ',d_bounds(l+1)
     &,   '      The stretch factor puts the bottom at ',bounds(l)+sum
     &,   ' cm with a bottom cell thickness of '
     &,   p5*(deltau(num+n-1) + deltau(num+n-2))
     &,   '      Adjust "stretch_z" in subroutine "gcoord" to get'
     &,   ' closer to the desired specifications if needed.'

        else
       if (my_pe == 0) 
     &    write (stdout,'(a,g14.7,a)')
     &   '    constructed ',an,' grid cells for this region'
       if (my_pe == 0) print*,' abs(an-n) = ',abs(an-n)
          if (abs(an-n) .gt. 0.01d0) then
       if (my_pe == 0) 
     &      write (stdout, '(/,a,i2,/,a,g14.7/,a,g14.7,a//a/a)')
     & '==>Error: non integral number of cells in region #',l
     &,'          average resolution within region =',avg_res
     &,'          this implies ',an,' grid cells'
     &,'          Change grid specifications within USER INPUT section'
     &,'          Here is some help...'
            d_new = (2.0*wid)/(n-1) - d_bounds(l)
       if (my_pe == 0) 
     &      write (stdout,'(/a,i4,a,i2,a,1pe14.7,a,1pe14.7/)')
     &      ' Note: to get ',n-1,' grid cells within region ',l
     &,     ', change resolution from ', d_bounds(l+1), ' to ', d_new
            d_new = (2.0*wid)/n - d_bounds(l)
       if (my_pe == 0) 
     &      write (stdout,'(/a,i4,a,i2,a,1pe14.7,a,1pe14.7/)')
     &      ' Note: to get ',n,' grid cells within region ',l
     &,     ', change resolution from ', d_bounds(l+1), ' to ', d_new
            d_new = (2.0*wid)/(n+1) - d_bounds(l)
       if (my_pe == 0) 
     &      write (stdout,'(/a,i4,a,i2,a,1pe14.7,a,1pe14.7/)')
     &      ' Note: to get ',n+1,' grid cells within region ',l
     &,     ', change resolution from ', d_bounds(l+1), ' to ', d_new
             stop
          endif
        endif
        num = num + n
      enddo
c
c     adjust "num" to reflect the total number of cells contained in
c     all regions
c
      num = num - 1
c
      do i=1,num
c
c       build resolution for T cells: "deltat". Note that
c       variable resolution (stretched grid) implies T points are
c       off center
c
        if (i .eq. 1) then
          deltat(i) = p5*(d_bounds(1) + deltau(i))
        else
          deltat(i) = p5*(deltau(i) + deltau(i-1))
        endif
      enddo     
c
c     add boundary points if needed
c
      if (nbpts .ne. 0) then
        do i=num,1,-1
          deltat(i+1) = deltat(i) 
          deltau(i+1) = deltau(i)
        enddo
        deltat(1)     = deltat(2)
        deltau(1)     = d_bounds(1)
        deltat(num+2) = deltat(num+1) 
        deltau(num+2) = deltau(num+1)
        num           = num + 2 
      endif
      end subroutine gcell





      subroutine gcell_centered_t 
     &                 (maxlen, n_bounds, bounds, d_bounds, nbpts
     &,                 num, deltat, deltau, stretch,my_pe,tol_in)
c
c=======================================================================
c
c              G R I D   C E L L   C O N S T R U C T I O N
c
c     A domain is composed of one or more regions:
c     Build "num" T  cells with resolution "deltat(n) n=1,num" 
c     within the domain composed of regions bounded by "bounds".
c     Also construct "num" U  cells of resolution "deltau(n) n=1,num"
c     with the relation between T and U cells given by:
c     deltau(n) = 0.5*(deltat(n+1) + deltat(n))
c     Resolution may be constant or smoothly varying within each
c     region AND there must be an integral number of grid cells within 
c     each region. The domain is the sum of all regions.
c
c     inputs:
c
c     maxlen   = maximum length of "deltat" and "deltau"
c     n_bounds = number of bounds needed to define the regions
c     bounds   = latitude, longitude, or depth at each bound
c     d_bounds = delta (resolution) at each of the "bounds"
c     nbpts    = number of extra boundary cells to add to the domain.
c                (usually one at the beginning and end)
c     stretch  = stretching factor for last region (should only be used
c                in the vertical to provide increased stretching of grid
c                points. "stretch" = 1.0 gives no increased stretching.
c                "stretch" = 1.2 gives increased stretching...etc
c             
c     outputs:
c
c     num    = total number of grid cells within the domain
c     deltau = resolution of U grid cells: n=1,num
c     deltat = resolution of T grid cells: n=1,num
c
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
c=======================================================================
c
      implicit none
      integer stdout,maxlen,n_bounds,num,nbpts,my_pe
      parameter( stdout=6)
      real deltat(maxlen), deltau(maxlen)
      real d_bounds(n_bounds), bounds(n_bounds), stretch,tol_in

      logical keep_going
      double precision p5,pi ,avg_res, chg_res, tol, wid, an
      double precision del,d_new,sum
      integer l,n,i,ierr,m
c
c     Set some constants
c
      p5 = 0.5d0
      pi = 4.0d0*atan(1.0d0)
c
c     Do all regions, one at a time, to construct the domain
c
      num  = 1
      do l=1,n_bounds-1
c
       if (my_pe == 0) 
     &  write (stdout,'(2x,a,i2,a,g14.7,a,g14.7,a,g14.7,a,g14.7,a)')
     & ' region # ',l,'  going from ',bounds(l),' (res=',d_bounds(l)
     &,') to ',  bounds(l+1),' (res=',d_bounds(l+1),')' 
c
c       avg_res = average resolution of T cells within region
c       chg_res = change in resolution across the region
c       wid     = width of region
c       tol     = tolerance for fitting T cels within region width
c
c       provide for stretching last region if needed
c
        if (l .eq. n_bounds-1) then
          avg_res = p5*(d_bounds(l) + stretch*d_bounds(l+1))
          chg_res = (stretch*d_bounds(l+1) - d_bounds(l))
        else
          avg_res = p5*(d_bounds(l) + d_bounds(l+1))
          chg_res = (d_bounds(l+1) - d_bounds(l))
        endif
c
c
c       had to change this ...
c        tol = 1.d-3
c      all the time
        tol=tol_in
c
        wid = abs(bounds(l+1) - bounds(l))
        an  = wid/avg_res
        m   = nint(an)


c
c       Calculate resolution of "T" cells: "deltat". Note that
c       "T" grid points are centered in these cells (as in MOM 1)
c       n = number of "t" cells fitting within the region boundaries
c
        sum = 0.0
	n   = 0
        do i = 1,100000
	  if (m.le.1) m=2
	  del = avg_res - p5*chg_res*cos((pi/(m-1))*(i-1))
	  if (sum + del .le. wid*(1.0 + tol)) then
	    sum = sum + del
	    if (num+i-1 .gt. maxlen) then
	      write (stdout,*) "=>Error: maxlen exceeded in gcell. "
     &,                        " ...increase size of maxlen"
	      stop
	    endif
	    deltat(num+i-1) = del
	    n = n + 1
	  else
	    go to 100
	  endif
	enddo 
 100    continue


        if (l .eq. n_bounds-1 .and. stretch .ne. 1.0) then
       if (my_pe == 0) 
     &    write (stdout,'(a,i3,a,f5.2)')
     & '    constructed ',n,' cells with a stretch factor of ', stretch
       if (my_pe == 0) 
     &    write (stdout,'(/2(a,g14.7),/2(a,g14.7),/a,a/)')
     &    'Note: you specified the ocean bottom at ',bounds(l+1)
     &,   ' cm with a bottom cell thickness of ',d_bounds(l+1)
     &,   '      The stretch factor puts the bottom at ',bounds(l)+sum
     &,   ' cm with a bottom cell thickness of ',deltat(num+n-1)
     &,   '      Adjust "stretch_z" in subroutine "gcoord" to get'
     &,   ' closer to the desired specifications if needed.'

        else
       if (my_pe == 0) 
     &    write (stdout,'(a,g14.7,a)')
     &   '    constructed ',an,' grid cells for this region'
       if (my_pe == 0) print*,' abs(an-n) = ',abs(an-n)
          if (abs(an-n) .gt. 0.01d0) then
       if (my_pe == 0) 
     &      write (stdout, '(/,a,i2,/,a,g14.7/,a,g14.7,a//a/a)')
     & '==>Error: non integral number of cells in region #',l
     &,'          average resolution within region =',avg_res
     &,'          this implies ',an,' grid cells'
     &,'          Change grid specifications within USER INPUT section'
     &,'          Here is some help...'
            d_new = (2.0*wid)/(n-1) - d_bounds(l)
       if (my_pe == 0) 
     &      write (stdout,'(/a,i4,a,i2,a,1pe14.7,a,1pe14.7/)')
     &      ' Note: to get ',n-1,' grid cells within region ',l
     &,     ', change resolution from ', d_bounds(l+1), ' to ', d_new
            d_new = (2.0*wid)/n - d_bounds(l)
       if (my_pe == 0) 
     &      write (stdout,'(/a,i4,a,i2,a,1pe14.7,a,1pe14.7/)')
     &      ' Note: to get ',n,' grid cells within region ',l
     &,     ', change resolution from ', d_bounds(l+1), ' to ', d_new
            d_new = (2.0*wid)/(n+1) - d_bounds(l)
       if (my_pe == 0) 
     &      write (stdout,'(/a,i4,a,i2,a,1pe14.7,a,1pe14.7/)')
     &      ' Note: to get ',n+1,' grid cells within region ',l
     &,     ', change resolution from ', d_bounds(l+1), ' to ', d_new
             stop
          endif
        endif
        num = num + n
      enddo
c
c     adjust "num" to reflect the total number of cells contained in
c     all regions
c
      num = num - 1
c
      do i=1,num
c
c       build resolution for "u" cells: "deltau". Note that
c       variable resolution (stretched grid) implies "u" points are
c       off center as in MOM 1
c
        if (i .eq. num) then
          deltau(i) = d_bounds(n_bounds)
        else
	  deltau(i) = p5*(deltat(i+1) + deltat(i))
        endif
      enddo     
c
c     add boundary points if needed
c
      if (nbpts .ne. 0) then
        do i=num,1,-1
          deltat(i+1) = deltat(i) 
          deltau(i+1) = deltau(i)
        enddo
        deltat(1)     = deltat(2)
        deltau(1)     = d_bounds(1)
        deltat(num+2) = deltat(num+1) 
        deltau(num+2) = deltau(num+1)
        num           = num + 2 
      endif
      end subroutine gcell_centered_t





      subroutine iso_grid (maxlen, nylats, y_lat, dy_lat, jmt
     &,                    dytdeg, dyudeg, yt, yu)
c
c=======================================================================
c     compute latitude resolution of grid cells to match the convergence
c     of meridians.
c     
ce    Note: this routine was taken from MOM 3 and seems to be buggy
c     
c     inputs:
c
c     maxlen   = maximum length of "dytdeg" and "dyudeg"
c     nylats   = should equal 2 (defines one region)
c     dy_lat   = latitudinal resolution of grid cell on equator
c     y_lat(1) = southern boundary of the domain (it will be adjusted
c                to fit an integral number of cells)
c     y_lat(2) = northern boundary of the domain (it will be adjusted
c                to fit an integral number of cells)
c
c     outputs:
c
c     jmt  = number of grid cells
c     dytdeg = resolution of T cells (degrees)
c     dyudeg = resolution of U cells (degrees)
c     yt     = latitude of point within T cell (degrees)
c     yu     - latitude of point within U cell (degrees)
c
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
c=======================================================================
c
      implicit none

      integer maxlen,nylats,jmt
      real y_lat(nylats), dy_lat(nylats)
      real dytdeg(maxlen), dyudeg(maxlen)
      real yt(maxlen), yu(maxlen)

      integer, parameter :: jmaxlat = 1000, lenjmax = 2*jmaxlat+1
      real dusq(-jmaxlat-1:jmaxlat+1), dtsq(-jmaxlat-1:jmaxlat+1)
      real usq(-jmaxlat-1:jmaxlat+1), tsq(-jmaxlat-1:jmaxlat+1)
      real y_bound(2), dy_bound(2)
c
      real pi, rrad
      integer n,n1,n2,j

      if (nylats .gt. 2) then
        write (6,'(/a,a/)')
     &  '=>Error: Must only have nylats=2 when option isotropic_grid'
     &,' is enabled'
        call abort()
      endif
      if (dy_lat(1) .ne. dy_lat(2)) then
        write (6,'(/a,a/)')
     &  '=>Error: dy_lat(1) must equal dy_lat(2) when option '
     &,' isotropic_grid is enabled'
        call abort()
      endif
      pi   = 4.0*atan(1.0)
      rrad = pi/180.0
c
c     build an isotropic grid
c
      usq(0)  = 0.0
      dusq(0) = dy_lat(1)
      tsq(1)  = 0.5*dusq(0)
      tsq(0)  = -tsq(1)
      dtsq(1) = dusq(0)*cos(tsq(1)*rrad)
      dtsq(0) = dtsq(1)
      do n=1,jmaxlat
        dusq(n)   = 2.0*dtsq(n) - dusq(n-1)
	usq(n)    = tsq(n) + 0.5*dusq(n)
	if (tsq(n) .lt. 90) then
	  tsq(n+1)  = tsq(n) + dusq(n)
	  dtsq(n+1) = dusq(0)*cos(tsq(n+1)*rrad)
	else
	  tsq(n+1)  = tsq(n)
	  dtsq(n+1) = 0.0
	endif
      enddo
      do n=-1,-jmaxlat,-1
        dusq(n) = dusq(-n)
	usq(n)  = -usq(-n)
	dtsq(n) = dtsq(-(n-1)) 
	tsq(n)  = -tsq(-(n-1))
      enddo
c
c     pick out cells between bounding latitudes
c
      n1 = -jmaxlat
      n2 = jmaxlat
      do n=-jmaxlat,jmaxlat
        if (usq(n) .lt. y_lat(1)) n1 = n
	if (usq(n) .lt. y_lat(2)) n2 = n+1
      enddo
c
c     re-define bounding latitudes to match isotropic boundaries
c
      y_lat(1) = usq(n1)
      y_lat(2) = usq(n2) 
c
      if (n1 .eq. -jmaxlat .or. n2 .eq. jmaxlat) then
        write (6,'(/a,a/)')
     &  '=>Error: Need to increase jmaxlat to reach the max latitude'
     &,' of the domain specified for option isotropic_grid.'
        call abort()
      endif

      jmt = n2 - n1 + 1
      do j=1,jmt
        dytdeg(j) = dtsq(j+n1-1)
	dyudeg(j) = dusq(j+n1-1)
        yu(j) = usq(j+n1-1)
        yt(j) = tsq(j+n1-1)
      enddo
      return
      end



      subroutine read_grid_bin()
      use prep_module
      implicit none
      integer  :: io,i,j,k
      character(len=80) :: FLAME_ID ='no flame id '

      real (kind=4), allocatable  :: xt4(:),xu4(:), yt4(:),yu4(:)
      real (kind=4), allocatable  :: zt4(:),zw4(:),dzt4(:),dzw4(:)
      real (kind=4), allocatable  :: dxtdeg4(:),dxudeg4(:)
      real (kind=4), allocatable  :: dytdeg4(:),dyudeg4(:)

      print*,' reading grid infos from file ',
     &       old_grid_file(1:len_trim(old_grid_file))

      call getunit(io,old_grid_file,'usr ieee')
      read (io) flame_id   ! this the the flame ID check
      read (io) imt, jmt, km

      allocate( dxtdeg(imt), dxudeg(imt), xt(imt), xu(imt))
      allocate( dytdeg(jmt), dyudeg(jmt), yt(jmt), yu(jmt))
      allocate( dzt(km), dzw(0:km), zt(km), zw(km))

      allocate( dxtdeg4(imt), dxudeg4(imt), xt4(imt), xu4(imt))
      allocate( dytdeg4(jmt), dyudeg4(jmt), yt4(jmt), yu4(jmt))
      allocate( dzt4(km), dzw4(0:km), zt4(km), zw4(km))

      read (io)(dxtdeg4(i),i=1,imt)
     &,        (dytdeg4(j),j=1,jmt)
     &,        (dxudeg4(i),i=1,imt)
     &,        (dyudeg4(j),j=1,jmt)
     &,        (dzt4(k),k=1,km)
     &,        (dzw4(k),k=0,km)
     &,        (xt4(i),i=1,imt)
     &,        (xu4(i),i=1,imt)
     &,        (yt4(j),j=1,jmt)
     &,        (yu4(j),j=1,jmt)
     &,        (zt4(k),k=1,km)
     &,        (zw4(k),k=1,km)
       close (io)

      dxtdeg=dxtdeg4; dxudeg=dxudeg4; xt= xt4; xu= xu4
      dytdeg=dytdeg4; dyudeg=dyudeg4; yt= yt4; yu= yu4
      dzt=dzt4; dzw= dzw4;zt= zt4; zw=zw4

      deallocate( dxtdeg4, dxudeg4, xt4, xu4)
      deallocate( dytdeg4, dyudeg4, yt4, yu4)
      deallocate( dzt4, dzw4, zt4, zw4)
      end subroutine


      subroutine write_grid_bin()
      use prep_module
      implicit none

      integer :: i,j,k,io
      character(len=80) :: FLAME_ID ='no flame id '
      real (kind=4) :: xt4(imt),xu4(imt), yt4(jmt),yu4(jmt)
      real (kind=4) :: zt4(km),zw4(km),dzt4(km),dzw4(0:km)
      real (kind=4) :: dxtdeg4(imt),dxudeg4(imt)
      real (kind=4) :: dytdeg4(jmt),dyudeg4(jmt)

      print*,' writing grid infos to file ',
     &       new_grid_file(1:len_trim(new_grid_file))

      dxtdeg4=dxtdeg; dxudeg4=dxudeg; xt4= xt; xu4= xu
      dytdeg4=dytdeg; dyudeg4=dyudeg; yt4= yt; yu4= yu
      dzt4=dzt; dzw4= dzw;zt4= zt; zw4=zw

      call getunit(io,new_grid_file,'usr ieee')
      write (io) flame_id   ! this the the flame ID check
      write (io) imt, jmt, km
      write (io) (dxtdeg4(i),i=1,imt)
     &,         (dytdeg4(j),j=1,jmt)
     &,         (dxudeg4(i),i=1,imt)
     &,         (dyudeg4(j),j=1,jmt)
     &,         (dzt4(k),k=1,km)
     &,         (dzw4(k),k=0,km)
     &,         (xt4(i),i=1,imt)
     &,         (xu4(i),i=1,imt)
     &,         (yt4(j),j=1,jmt)
     &,         (yu4(j),j=1,jmt)
     &,         (zt4(k),k=1,km)
     &,         (zw4(k),k=1,km)
       close (io)

       end subroutine



      subroutine sub_isotro(maxlen, n_bounds, bounds, d_bounds, 
     &                      num, deltat, deltau)
c
c=======================================================================
c
c Subroutine to determine the fine resolution DYNAMO-2 grid
c The ARAKAWA B grid is a mesh of velocity grid cells, with tracer
c points in the corners of a cell. The zonal width of a velocity grid
c cell dxu is given by the distance of the two surrounding meridians,
c while the meridional extent is given by the distance of the two 
c westernmost tracer points:
c
c      +->  T(i,j+1) -------------- T(i+1,j+1)
c      |    |                                |
c   dyu(j)  |              U(i,j)            |
c      |    |                                |
c      +->  T(i,j) ------------------ T(i+1,j)
c
c           ^                                ^
c           |                                |
c           +------------ dxu(i) ------------+
c
c
c There arises no problem from the zonal grid, because it is defined to
c be a regular third (or four third in the coarse resolution case) 
c degree grid. To yield an isotropic grid, the meridional extent of a 
c grid cell has to be a function of latitude. Let phi_s denote
c the southern edge of the gridcell, and phi_n the northern one. We are 
c now looking for the latitude phi_c of the center of the gridcell,
c which is defined by 
c
c  phi_s + 0.5*dyu = phi_c   ; dyu=dxu*cos(phi_c)      and thus
c  phi_s + 0.5*dxu*cos(phi_c) - phi_c = 0
c
c This equation will be solved by an iteration. We will start at the 
c Equator, proceeding north until bounds(2) is reached. Then, all until
c bounds(1) will be mirrored at the Equator.
c
c     The relation between "t" and "u" cells is given by:
c     deltat(n) = 0.5*(deltau(n-1) + deltau(n))
c
c     There must be an integral number of grid cells within 
c     the region to provide an easy switch between coarse and fine 
c     grid-resolution. 
c
c     inputs:
c
c     maxlen   = maximum length of "deltat" and "deltau"
c     n_bounds = number of bounds needed to define the regions
c     bounds   = latitude, longitude, or depth at each bound
c     d_bounds = delta (resolution) at the equator: d_bounds(1)
c     nbpts    = number of extra boundary cells to add to the domain.
c
c     outputs:
c
c     num      = total number of grid cells within the domain
c     deltau   = resolution of "u" grid cells: n=1,num
c     deltat   = resolution of "t" grid cells: n=1,num
c     bounds   = new latitude at each bound
c     d_bounds = resolution between yu(0) and yt(1): d_bounds(1)
c
c     author:      j.kroeger    e-mail=> jkroeger@ifm.uni-kiel.de
c=======================================================================
c
c
      implicit none
      integer :: nbpts = 0, stdout = 6, maxlen, n_bounds, num
      integer :: i,j,k,jmt,nrc,itc,nrs
      real :: hdtone,endlat,soubou,stlat,dummy,p5
      real deltat(maxlen), deltau(maxlen)
      real d_bounds(n_bounds), bounds(n_bounds)
      double precision dely(maxlen)
      double precision hdy0,pi,pifac,phi_n,phi_s,phi_c,abort,
     &                 hdcos,residual,dxt,dynew
c
      dummy=-9999.0
      dxt=d_bounds(1)
      hdy0=0.5*dxt         ! half dy cosine amplitude
      p5 = 0.5
      pi = 4.0*atan(1.0)
      pifac=pi/180.0
      nrc=0                ! gridcell counter
      abort=1.0d-12        ! should be small enough
      phi_n=hdy0           ! let`s start north of the Equator
  100 phi_s=phi_n          ! define new southern edge
      phi_c=phi_s+hdy0*dcos(phi_s*pifac)
      itc=0                ! iteration counter
  200 hdcos=hdy0*dcos(phi_c*pifac)
      residual=phi_c-hdcos-phi_s
      itc=itc+1            ! increase iteration counter 
      if(dabs(residual).lt.abort)then
       goto 300           ! fine, we did convergence
      else
       phi_c=phi_c-0.5*residual ! new estimate
       goto 200                 ! next iterative step
      endif
  300 phi_n=phi_s+2.0*hdcos     ! new northern edge
      nrc=nrc+1                 ! nr of grid cells
      dely(nrc)=2.0*hdcos
ce      if ((phi_c.ge.bounds(2)).and.(mod(nrc,4).eq.0)) goto 500
      if ((phi_c.ge.bounds(2))) goto 500
      goto 100
  500 continue
      phi_s=hdy0      ! hit northern edge, now look for southern
      nrs=0
      soubou=abs(bounds(1))
  600 nrs=nrs+1       ! number of tracer cells south of the Equator
      phi_c=phi_s+p5*dely(nrs)
      phi_s=phi_s+dely(nrs)
c      if((phi_c.ge.soubou).and.(mod(nrs,4).eq.0))goto 700
      if((phi_c.ge.soubou))goto 700
      goto 600
  700 continue
      stlat=-phi_c    ! hit southern edge
      do i=nrc,1,-1   ! shift northern hemisphere 
       dely(i+nrs)=dely(i)
      end do
      dely(nrs)=2*hdy0
      do i=1,nrs      ! mirror southern region at Equator
       dely(i)=dely(2*nrs-i)
      end do
      jmt=nrs+nrc
c      build resolution for "t" cells: "deltat".
      do j=1,jmt
       deltau(j)=dely(j)
      enddo
      deltat(1) = p5*(dely(2*nrs)+deltau(1))
      do i=2,jmt
       deltat(i) = p5*(deltau(i-1) + deltau(i))
      enddo
      hdtone=p5*deltau(2*nrs)  ! halfdeltatone
      write (stdout,'(a,g14.7,a)')
     &   '    constructed ',jmt,' meridional grid cells '
      endlat=stlat
      do j=1,jmt
        endlat=endlat+deltat(j)
      end do
      write (stdout,'(a,g14.7,a,g14.7)')
     &   ' startlatitude =    ',stlat,
     &   ' endlatitude   =    ',endlat
      bounds(1)=stlat
      bounds(2)=endlat
      d_bounds(1)=hdtone
      num=jmt
      end subroutine
c


      subroutine isotro_centered_t
     &   (maxlen, n_bounds, bounds, d_bounds,num, deltat, deltau)
c
c=======================================================================
c
c Subroutine to determine the fine resolution DYNAMO-2 grid
c The ARAKAWA B grid is a mesh of tracer grid cells, with velocity 
c points in the corners of a cell. The zonal width of a tracer grid
c cell dxt is given by the distance of the two surrounding meridians,
c while the meridional extent is given by the distance of the two 
c easternmost velocity points:
c
c      +->  U(i-1,j  ) ------------ U(i  ,j  )
c      |    |                                |
c   dyt(j)  |              T(i,j)            |
c      |    |                                |
c      +->  U(i-1,j-1) ------------ U(i  ,j-1)
c
c           ^                                ^
c           |                                |
c           +------------ dxt(i) ------------+
c
c There arises no problem from the zonal grid, because it is defined to
c be a regular third (or four third in the coarse resolution case) 
c degree grid. To yield an isotropic grid, the meridional extent of a 
c grid cell has to be a function of latitude. Let phi_s denote
c the southern edge of the gridcell, and phi_n the northern one. We are 
c now looking for the latitude phi_c of the center of the gridcell,
c which is defined by 
c
c  phi_s + 0.5*dyt = phi_c   ; dyt=dxt*cos(phi_c)      and thus
c  phi_s + 0.5*dxt*cos(phi_c) - phi_c = 0
c
c This equation will be solved by an iteration. We will start at the 
c Equator, proceeding north until bounds(2) is reached. Then, all until
c bounds(1) will be mirrored at the Equator.
c
c     The relation between "t" and "u" cells is given by:
c     deltau(n) = 0.5*(deltat(n+1) + deltat(n))
c
c     There must be an integral number of grid cells within 
c     the region to provide an easy switch between coarse and fine 
c     grid-resolution. 
c
c     inputs:
c
c     maxlen   = maximum length of "deltat" and "deltau"
c     n_bounds = number of bounds needed to define the regions
c     bounds   = latitude, longitude, or depth at each bound
c     d_bounds = delta (resolution) at the equator: d_bounds(1)
c     nbpts    = number of extra boundary cells to add to the domain.
c                (usually one at the beginning and end)
c     outputs:
c
c     num      = total number of grid cells within the domain
c     deltau   = resolution of "u" grid cells: n=1,num
c     deltat   = resolution of "t" grid cells: n=1,num
c     bounds   = new latitude at each bound
c     d_bounds = resolution between yu(0) and yt(1): d_bounds(1)
c
c     author:      j.kroeger    e-mail=> jkroeger@ifm.uni-kiel.de
c=======================================================================
c
      implicit none
      integer :: nbpts = 0, stdout = 6, maxlen, n_bounds, num
      integer :: i,j,k,jmt,nrc,itc,nrs
      real :: hdtone,endlat,soubou,stlat,dummy,p5
      real deltat(maxlen), deltau(maxlen)
      real d_bounds(n_bounds), bounds(n_bounds)
      double precision dely(maxlen)
      double precision hdy0,pi,pifac,phi_n,phi_s,phi_c,abort,
     &                 hdcos,residual,dxt,dynew
c
      dummy=-9999.0
      dxt=d_bounds(1)
      hdy0=0.5*dxt         ! half dy cosine amplitude
      p5 = 0.5
      pi = 4.0*atan(1.0)
      pifac=pi/180.0
      nrc=0                ! gridcell counter
      abort=1.0d-12        ! should be small enough
      phi_n=0.0            ! let`s start at the Equator
  100 phi_s=phi_n          ! define new southern edge
      phi_c=phi_s+hdy0*dcos(phi_s*pifac)
      itc=0                ! iteration counter
  200 hdcos=hdy0*dcos(phi_c*pifac)
      residual=phi_c-hdcos-phi_s
      itc=itc+1            ! increase iteration counter 
      if(dabs(residual).lt.abort)then
       goto 300           ! fine, we did convergence
      else
       phi_c=phi_c-0.5*residual ! new estimate
       goto 200                 ! next iterative step
      endif
  300 phi_n=phi_s+2.0*hdcos     ! new northern edge
      nrc=nrc+1                 ! nr of grid cells
      dely(nrc)=2.0*hdcos
      if ((phi_n.ge.bounds(2)).and.(mod(nrc,4).eq.0)) goto 500
      goto 100
  500 continue
      phi_s=0.0       ! hit northern edge, now look for southern
      nrs=0
      soubou=abs(bounds(1))
  600 nrs=nrs+1       ! number of tracer cells south of the Equator
      phi_c=phi_s+p5*dely(nrs)
      phi_s=phi_s+dely(nrs)
      if((phi_s.ge.soubou).and.(mod(nrs,4).eq.0))goto 700
      goto 600
  700 continue
      stlat=-phi_s    ! hit southern edge
      do i=nrc,1,-1   ! shift northern hemisphere 
       dely(i+nrs)=dely(i)
      end do
      dely(nrs)=2*hdy0
      do i=1,nrs      ! mirror southern region at Equator
       dely(i)=dely(2*nrs-i+1)
      end do
c
      jmt=nrs+nrc
c
c       build resolution for "u" cells: "deltau". Note that
c       variable resolution (stretched grid) implies "u" points are
c       off center as in MOM 1
c
      do j=1,jmt
       deltat(j)=dely(j)
      enddo
c
      hdtone=p5*deltat(1)  ! halfdeltatone
c
      do i=1,jmt-1
       deltau(i) = p5*(deltat(i+1) + deltat(i))
      enddo
      deltau(jmt) = deltat(jmt)
c
c
          write (stdout,'(a,g14.7,a)')
     &   '    constructed ',jmt,' meridional grid cells '
c
      endlat=stlat
      do j=1,jmt
        endlat=endlat+deltat(j)
      end do
c
          write (stdout,'(a,g14.7,a,g14.7)')
     &   ' startlatitude =    ',stlat,
     &   ' endlatitude   =    ',endlat
c
      bounds(1)=stlat
      bounds(2)=endlat
      d_bounds(1)=hdtone
      num=jmt
      end subroutine
c


      subroutine gcell_old (maxlen, n_bounds, bounds, d_bounds, nbpts
     &,                 num, deltat, deltau, stretch,my_pe)
c
c=======================================================================
c
c              G R I D   C E L L   C O N S T R U C T I O N
c
c     A domain is composed of one or more regions:
c     Build "num" T  cells with resolution "deltat(n) n=1,num" 
c     within the domain composed of regions bounded by "bounds".
c     Also construct "num" U  cells of resolution "deltau(n) n=1,num"
c     with the relation between T and U cells given by:
c     deltat(n) = 0.5*(deltau(n-1) + deltau(n))
c     Resolution may be constant or smoothly varying within each
c     region AND there must be an integral number of grid cells within 
c     each region. The domain is the sum of all regions.
c
c     inputs:
c
c     maxlen   = maximum length of "deltat" and "deltau"
c     n_bounds = number of bounds needed to define the regions
c     bounds   = latitude, longitude, or depth at each bound
c     d_bounds = delta (resolution) at each of the "bounds"
c     nbpts    = number of extra boundary cells to add to the domain.
c                (usually one at the beginning and end)
c     stretch  = stretching factor for last region (should only be used
c                in the vertical to provide increased stretching of grid
c                points. "stretch" = 1.0 gives no increased stretching.
c                "stretch" = 1.2 gives increased stretching...etc
c             
c     outputs:
c
c     num    = total number of grid cells within the domain
c     deltau = resolution of U grid cells: n=1,num
c     deltat = resolution of T grid cells: n=1,num
c
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
c=======================================================================
c
      implicit none
      integer stdout,maxlen,n_bounds,num,nbpts,my_pe
      parameter( stdout=6)
      real deltat(maxlen), deltau(maxlen)
      real d_bounds(n_bounds), bounds(n_bounds), stretch

      logical keep_going
      double precision p5,pi ,avg_res, chg_res, tol, wid, an
      double precision del,d_new,sum
      integer l,n,i,ierr,m
c
c     Set some constants
c
      p5 = 0.5d0
      pi = 4.0d0*atan(1.0d0)
c
c     Do all regions, one at a time, to construct the domain
c
      num  = 1
      do l=1,n_bounds-1
c
       if (my_pe == 0) 
     &  write (stdout,'(2x,a,i2,a,g14.7,a,g14.7,a,g14.7,a,g14.7,a)')
     & ' region # ',l,'  going from ',bounds(l),' (res=',d_bounds(l)
     &,') to ',  bounds(l+1),' (res=',d_bounds(l+1),')' 
c
c       avg_res = average resolution of T cells within region
c       chg_res = change in resolution across the region
c       wid     = width of region
c       tol     = tolerance for fitting T cels within region width
c
c       provide for stretching last region if needed
c
        if (l .eq. n_bounds-1) then
          avg_res = p5*(d_bounds(l) + stretch*d_bounds(l+1))
          chg_res = (stretch*d_bounds(l+1) - d_bounds(l))
        else
          avg_res = p5*(d_bounds(l) + d_bounds(l+1))
          chg_res = (d_bounds(l+1) - d_bounds(l))
        endif
c
c
c       had to change this ...
c        tol = 1.d-5
        tol = 1.d-3
c
        wid = abs(bounds(l+1) - bounds(l))
        an  = wid/avg_res
        m   = nint(an)
c
c       Calculate resolution of U cells: "deltau"
c       U grid points will be centered in these cells
c       n = number of T cells fitting within the region boundaries
c       note: "sum" initially discounts half of the U cells widths
c       at the boundaries
c
        sum = p5*d_bounds(l) - p5*d_bounds(l+1)
        n   = 0
        i = 0
        keep_going = .true.
        do while (i .le. 100000 .and. keep_going)
          i = i + 1
          del = avg_res - p5*chg_res*cos((pi/m)*i)
          if (sum + del .le. wid*(1.0d0 + tol)) then
            sum = sum + del
            if (num+i-1 .gt. maxlen) then
              write (stdout,*) "=>Error: maxlen exceeded in gcell. "
     &,                        " ...increase size of maxlen"
              call halt_stop(' in gcell')
            endif
            deltau(num+i-1) = del
            n = n + 1
          else
            keep_going = .false.
          endif
        enddo 

        if (l .eq. n_bounds-1 .and. stretch .ne. 1.0) then
       if (my_pe == 0) 
     &    write (stdout,'(a,i3,a,f5.2)')
     & '    constructed ',n,' cells with a stretch factor of ', stretch
       if (my_pe == 0) 
     &    write (stdout,'(/2(a,g14.7),/2(a,g14.7),/a,a/)')
     &    'Note: you specified the ocean bottom at ',bounds(l+1)
     &,   ' cm with a bottom cell thickness of ',d_bounds(l+1)
     &,   '      The stretch factor puts the bottom at ',bounds(l)+sum
     &,   ' cm with a bottom cell thickness of '
     &,   p5*(deltau(num+n-1) + deltau(num+n-2))
     &,   '      Adjust "stretch_z" in subroutine "gcoord" to get'
     &,   ' closer to the desired specifications if needed.'

        else
       if (my_pe == 0) 
     &    write (stdout,'(a,g14.7,a)')
     &   '    constructed ',an,' grid cells for this region'
       if (my_pe == 0) print*,' abs(an-n) = ',abs(an-n)
          if (abs(an-n) .gt. 0.01) then
       if (my_pe == 0) 
     &      write (stdout, '(/,a,i2,/,a,g14.7/,a,g14.7,a//a/a)')
     & '==>Error: non integral number of cells in region #',l
     &,'          average resolution within region =',avg_res
     &,'          this implies ',an,' grid cells'
     &,'          Change grid specifications within USER INPUT section'
     &,'          Here is some help...'
            d_new = (2.0*wid)/(n-1) - d_bounds(l)
       if (my_pe == 0) 
     &      write (stdout,'(/a,i4,a,i2,a,1pe14.7,a,1pe14.7/)')
     &      ' Note: to get ',n-1,' grid cells within region ',l
     &,     ', change resolution from ', d_bounds(l+1), ' to ', d_new
            d_new = (2.0*wid)/n - d_bounds(l)
       if (my_pe == 0) 
     &      write (stdout,'(/a,i4,a,i2,a,1pe14.7,a,1pe14.7/)')
     &      ' Note: to get ',n,' grid cells within region ',l
     &,     ', change resolution from ', d_bounds(l+1), ' to ', d_new
            d_new = (2.0*wid)/(n+1) - d_bounds(l)
       if (my_pe == 0) 
     &      write (stdout,'(/a,i4,a,i2,a,1pe14.7,a,1pe14.7/)')
     &      ' Note: to get ',n+1,' grid cells within region ',l
     &,     ', change resolution from ', d_bounds(l+1), ' to ', d_new
            call sub_flush(6)
            call halt_stop(' in gcell')
          endif
        endif
        num = num + n
      enddo
c
c     adjust "num" to reflect the total number of cells contained in
c     all regions
c
      num = num - 1
c
      do i=1,num
c
c       build resolution for T cells: "deltat". Note that
c       variable resolution (stretched grid) implies T points are
c       off center
c
        if (i .eq. 1) then
          deltat(i) = p5*(d_bounds(1) + deltau(i))
        else
          deltat(i) = p5*(deltau(i) + deltau(i-1))
        endif
      enddo     
c
c     add boundary points if needed
c
      if (nbpts .ne. 0) then
        do i=num,1,-1
          deltat(i+1) = deltat(i) 
          deltau(i+1) = deltau(i)
        enddo
        deltat(1)     = deltat(2)
        deltau(1)     = d_bounds(1)
        deltat(num+2) = deltat(num+1) 
        deltau(num+2) = deltau(num+1)
        num           = num + 2 
      endif
      end subroutine gcell_old

