c23456789012345678901234567890123456789012345678901234567890123456789012


      SUBROUTINE modgrid(iopt,opt)

c     ******************************************************************
c
c     This routine establishes the model grid
c
c     ******************************************************************

c     -------------
c     Declarations
c     -------------

c     ''''''''''''''''''''''
      implicit none
      include 'commonvar.f'      
c     ''''''''''''''''''''''

c--   Argument variables

      integer	iopt

      character	opt*5


c--   Internal variables

      integer	i,j,k,jump,i1,i2,j1,j2

      real*8	z0,zm

      logical	flag


c     ETOPO5 (NH)

      integer	idim_nh,jdim_nh,kco_nh

      parameter (idim_nh=1801,
     >		 jdim_nh= 961)

      real*8    z_nh  (idim_nh,jdim_nh),
     >		lon_nh(idim_nh        ),
     >	        lat_nh(        jdim_nh)


c     ETOPO5 (SH)

      integer   idim_sh,jdim_sh,kco_sh

      parameter (idim_sh=1801,
     >           jdim_sh= 961)

      real*8    z_sh  (idim_sh,jdim_sh),
     >          lon_sh(idim_sh        ),
     >          lat_sh(        jdim_sh)


c     ETOPO5 (roughness)

      integer	kcomp,kco

      real*8	z_mean




c     -------------------------------------
c     Atlantic domain (5-degree resolution) 
c     -------------------------------------

c---- Test

c     ++++++++++++++++++++++++
      if (opt.eq.'AT_05') then
c     ++++++++++++++++++++++++


c---- Zonal length of grid cells

      do i=1,itot
         dx(i) = 5.d0*rad
      enddo


c---- Meridional length of grid cells

      do j=1,jtot
         dy(j) = 5.d0*rad
      enddo


c---- Vertical length of grid cells

      do k=1,ktot
         dz(k) = 1000.d0
      enddo


c---- Longitudes and latitudes

      lon_w = -75.d0 

      lon(1) = lon_w*rad 
      do i=1+1,itot+1
         lon(i) = lon(i-1) + dx(i-1)
      enddo

c--

      lat_s = -45.d0

      lat(1) = lat_s*rad 
      do j=1+1,jtot+1
         lat(j) = lat(j-1) + dy(j-1)
      enddo


c---- End of test

c     ++++++
      endif
c     ++++++




c     -------------------------------------
c     Atlantic domain (3-degree resolution) 
c     -------------------------------------

c---- Test

c     ++++++++++++++++++++++++
      if (opt.eq.'AT_03') then
c     ++++++++++++++++++++++++


c---- Zonal length of grid cells

      do i=1,itot
         dx(i) = 3.d0*rad
      enddo


c---- Meridional length of grid cells

      do j=1,jtot
         dy(j) = 3.d0*rad
      enddo


c---- Vertical length of grid cells

      do k=1,ktot
         dz(k) = 1000.d0
      enddo


c---- Longitudes and latitudes

      lon_w = -75.d0 

      lon(1) = lon_w*rad 
      do i=1+1,itot+1
         lon(i) = lon(i-1) + dx(i-1)
      enddo

c--

      lat_s = -45.d0

      lat(1) = lat_s*rad 
      do j=1+1,jtot+1
         lat(j) = lat(j-1) + dy(j-1)
      enddo


c---- End of test

c     ++++++
      endif
c     ++++++




c     -------------------------------------
c     Atlantic domain (3-degree resolution)
c     -------------------------------------

c---- Test

c     ++++++++++++++++++++++++
      if (opt.eq.'AT003') then
c     ++++++++++++++++++++++++


c---- Zonal length of grid cells

      do i=1,itot
         dx(i) = 3.d0*rad
      enddo


c---- Meridional length of grid cells

      do j=1,jtot
         dy(j) = 3.d0*rad
      enddo


c---- Vertical length of grid cells

      dz(1) = 1000.d0   ! AABW
      dz(2) = 1250.d0   ! NADW (lower layer)
      dz(3) = 1500.d0   ! NADW (upper layer)


c---- Depth levels

      z_l(1) = z_b
      
      do k=1+1,ktot+1
         z_l(k) = z_l(k-1) - dz(k-1)
      enddo 


c---- Longitudes and latitudes

      lon_w = -75.d0

      lon(1) = lon_w*rad
      do i=1+1,itot+1
         lon(i) = lon(i-1) + dx(i-1)
      enddo

c--

      lat_s = -45.d0

      lat(1) = lat_s*rad
      do j=1+1,jtot+1
         lat(j) = lat(j-1) + dy(j-1)
      enddo


c---- End of test

c     ++++++
      endif
c     ++++++




c     -------------------------------------
c     Atlantic domain (2-degree resolution)
c     -------------------------------------

c---- Test

c     ++++++++++++++++++++++++
      if (opt.eq.'AT002') then
c     ++++++++++++++++++++++++


c---- Zonal length of grid cells

      do i=1,itot
         dx(i) = 2.d0*rad
      enddo


c---- Meridional length of grid cells

      do j=1,jtot
         dy(j) = 2.d0*rad
      enddo


c---- Vertical length of grid cells

      dz(1) = 1000.d0   ! AABW
      dz(2) = 1250.d0	! NADW (lower layer)
      dz(3) = 1500.d0	! NADW (upper layer)


c---- Depth levels

      z_l(1) = z_b
      
      do k=1+1,ktot+1
         z_l(k) = z_l(k-1) - dz(k-1)
      enddo 


c---- Longitudes and latitudes

      lon_w = -77.d0

      lon(1) = lon_w*rad
      do i=1+1,itot+1
         lon(i) = lon(i-1) + dx(i-1)
      enddo

c--

      lat_s = -46.d0

      lat(1) = lat_s*rad
      do j=1+1,jtot+1
         lat(j) = lat(j-1) + dy(j-1)
      enddo


c---- End of test

c     ++++++
      endif
c     ++++++




c     ----------------------
c     Read ETOPO5 bathymetry
c     ----------------------

c---- Read NH 

      if (iopt.eq.0) then
       open (20,file=
     > '~omarchal/Data/Bathymetry/ETOPO5/seltoAAAa003hU.d',status='old') 
      else
       open (20,file='seltoAAAa003hU.d',status='old')
      endif

      do jump=1,20
         read(20,*)
      enddo

      do j=jdim_nh,1,-1
         read(20,*) (z_nh(i,j),i=1,idim_nh)
      enddo

      close (20)


c---- Read SH

      if (iopt.eq.0) then
       open (20,file=
     > '~omarchal/Data/Bathymetry/ETOPO5/seltoAAAa003kH.d',status='old') 
      else
       open (20,file='seltoAAAa003kH.d',status='old')
      endif 

      do jump=1,20
         read(20,*)
      enddo

      do j=jdim_sh,1,-1
         read(20,*) (z_sh(i,j),i=1,idim_sh)
      enddo

      close (20)


c---- Make values below sea level positive

      do i=1,idim_nh
      do j=1,jdim_nh
         z_nh(i,j) = -z_nh(i,j)
      enddo
      enddo

c--

      do i=1,idim_sh
      do j=1,jdim_sh
         z_sh(i,j) = -z_sh(i,j)
      enddo
      enddo


c---- Coordinates of depth values

      do i=1,idim_nh
         lon_nh(i) = (-90.+(i-1)*5.d0/60.d0)*rad
      enddo

      do j=1,jdim_nh
         lat_nh(j) = (  0.+(j-1)*5.d0/60.d0)*rad
      enddo

c--

      do i=1,idim_sh
         lon_sh(i) = (-90.+(i-1)*5.d0/60.d0)*rad
      enddo

      do j=1,jdim_sh
         lat_sh(j) = (-80.+(j-1)*5.d0/60.d0)*rad
      enddo




c     ------------------------------------------------
c     Average of ETOPO5 depths in each model grid cell
c     ------------------------------------------------

c++
      do i1=1,itot
      do j1=1,jtot
c++

         kco_nh   = 0
         kco_sh   = 0
         z(i1,j1) = 0.d0

c--
         do i2=1,idim_nh

            if (lon_nh(i2).ge.lon(i1  ) .and.
     >          lon_nh(i2).lt.lon(i1+1)) then
              do j2=1,jdim_nh
                 if (lat_nh(j2).ge.lat(j1  ) .and.  
     >               lat_nh(j2).lt.lat(j1+1)) then
                     kco_nh   = kco_nh   + 1
                     z(i1,j1) = z(i1,j1) + z_nh(i2,j2)
                 endif
              enddo
             endif

          enddo

          if (kco_nh.ne.0) z(i1,j1) = z(i1,j1) / dfloat(kco_nh)
c--

          do i2=1,idim_sh

            if (lon_sh(i2).ge.lon(i1  ) .and.
     >          lon_sh(i2).lt.lon(i1+1)) then
              do j2=1,jdim_sh-1 
                 if (lat_sh(j2).ge.lat(j1  ) .and.
     >               lat_sh(j2).lt.lat(j1+1)) then
                     kco_sh   = kco_sh   + 1
                     z(i1,j1) = z(i1,j1) + z_sh(i2,j2)
                 endif
              enddo
             endif

          enddo
       
          if (kco_sh.ne.0) z(i1,j1) = z(i1,j1) / dfloat(kco_sh)
c--

c++
      enddo
      enddo
c++




c     ------------------------------------------
c     Round model grid cell depths to nearest dz
c     ------------------------------------------

      do i=1,itot
      do j=1,jtot

         flag = .false.
         do k=1,ktot
            zm = (z_l(k)+z_l(k+1)) / 2.
            if (z(i,j).lt.zm) then
                z0   = z_l(k+1)
                flag =  .true.
            endif
         enddo

         if (flag) then
           z(i,j) = z0
         else
           z(i,j) = z_b
         endif

      enddo
      enddo




c     ---------------------
c     End of the subroutine
c     ---------------------

      return

      end

c     ******************************************************************
