!======================================================================
! CROCO is a branch of ROMS developped at IRD, INRIA, 
! Ifremer, CNRS and Univ. Toulouse III  in France
! The two other branches from UCLA (Shchepetkin et al)
! and Rutgers University (Arango et al) are under MIT/X style license.
! CROCO specific routines (nesting) are under CeCILL-C license.
!
! CROCO website : http://www.croco-ocean.org
!======================================================================
!
#include "cppdefs.h"
#ifdef ANA_GRID
      subroutine ana_grid (tile)
      implicit none
# include "param.h"
      integer tile, trd
C$    integer omp_get_thread_num
# include "compute_tile_bounds.h"
      call ana_grid_tile (Istr,Iend,Jstr,Jend)
      return
      end

      subroutine ana_grid_tile (Istr,Iend,Jstr,Jend)
!
!  Set up model grid using analytical expressions:
!-------------------------------------------------------------------
!
!  INPUT: Grid configuration parameters
!  ------------------------------------
! Length_XI      are the physical dimensions of the computational
! Length_ETA     domain [usually measured in  meters];
!
! depth          is the maximum depth [meters, positive];
! f0,beta        are coriolis parameters which set up a beta-plane
!                                                  [1/s, 1/(m*s)].
!
!  OUTPUT: stored in common blocks, see files "scalars.h" "grid.h"
!  ---------------------------------------------------------------
!  xl,el      Physical dimensions of the computational domain
!                                    [usually measured in  meters];
!  h          Model bathymetry [meters, positive] at RHO-points.
!  hmin,hmax  Minimum and maximum values of depth of bathymetry [m].
!  f          Coriolis parameter (1/seconds) at RHO-points.
!  pm,pn      Coordinate transformation metric "m" [1/meters]
!                    associated with the differential distances in
!                    XI- and ETA-directions, both are at RHO-points.
!  xp,xr      XI-coordinates [m] at PSI- and RHO-points.
!  yp,yr      ETA-coordinates [m] at PSI- and RHO-points.
!-------------------------------------------------------------------
!
# ifdef AGRIF
      use Agrif_Util
# endif
      implicit none
      integer Istr,Iend,Jstr,Jend, i,j
# include "param.h"
# include "grid.h"
# include "scalars.h"
!
      real Length_XI, Length_ETA, depth, f0, beta,
     &     cff, y, x0, y0, dx, dy, NSUB
!
# ifdef DUNE
#  ifdef DUNE3D
      real Rij,cffi,cffj,cff1
#  elif defined ANA_DUNE
      real cff2
#  endif
# endif
# ifdef TIDAL_FLAT
      real depth_shallow, depth_deep
# endif
# ifdef ESTUARY
      real cff2, l2n, l2ml1m, l3ml1m
      integer jmiddle, imouth, itmp
      real hl1imouth, hl1mimouth, shtimouth, l1imouth, l2imouth, l3imouth
      real m, estCWM, estCWM_para, estHupstream, estHmouth, estHoff
      real estMS, estMS_para, estMS_para2, estCS, estCS2, estGamma2
      real xboti, l1i, l2i, l3i
      real hl1i, hl1mi, shti
      real ybotj
# endif
# ifdef INTERNAL
      real cff1, ridge_width, ridge_height
# endif
# ifdef RIP
      real h0,per,xx,yy,xs,db,alpha,lambda,eps
# endif
# ifdef SANDBAR
      real p1,p2,p3,p4,p5,p6,p7,xx
# endif
# if defined RIVER || defined SEAGRASS
      integer is,i0,j0
#  include "sources.h"
# endif
# ifdef MPI
#  define LOCALLM Lmmpi
#  define LOCALMM Mmmpi
# else
#  define LOCALLM Lm
#  define LOCALMM Mm
# endif
!
# include "compute_extended_bounds.h"


!
!----------------------------------------------------------------------
! Configuration parameters (h,f,beta)
!----------------------------------------------------------------------
!
# if defined BASIN
                     depth=5000.
                     f0=1.E-4
                     beta=2.E-11
# elif defined SINGLE_COLUMN
#  if defined KATO_PHILIPS || defined WILLIS_DEARDORFF
                     depth=50.
                     f0=0.
                     beta=0.
#  elif defined DIURNAL_CYCLE
                     depth=150.
                     f0=0.
                     beta=0.
#  elif defined FORCED_EKBBL
                     depth=1500.
                     f0=1.E-4
                     beta=0.
#  elif defined FORCED_DBLEEK
                     depth=25.0
                     f0=1.E-4
                     beta=0.
#  elif defined FORCED_NONROTBBL || defined FORCED_OSCNONROTBBL
                     depth=5.
                     f0=0.0
                     beta=0.
#  endif
# elif defined CANYON
                     depth=4000.
                     f0=1.E-4
                     beta=0.
# elif defined EQUATOR
                     depth=5000.
                     f0=0.
                     beta=2.2829e-11
# elif defined KH_INST
                     depth=256.
                     f0=0.
                     beta=0.
# elif defined GRAV_ADJ
#  ifdef NBQ
                     depth=0.3
#  else
                     depth=20.
#  endif
                     f0=0.
                     beta=0.
# elif defined ISOLITON
                     depth=0.29 ! Horn et al. (2001)
                     f0=0.
                     beta=0.
# elif defined ACOUSTIC
                     depth=128.
                     f0=0.
                     beta=0.
# elif defined INNERSHELF
                     depth=0.
                     f0=4*pi/86400*sin(-21*pi/180)
                     beta=0.
# elif defined INTERNAL
                     depth=2000.
                     f0=1.e-4
                     beta=0.
                     ridge_width=30.e3
                     ridge_height=1600.
# elif defined TS_HADV_TEST
                     depth=10.
                     f0=0.
                     beta=0.
# elif defined OVERFLOW
                     depth=40.
                     f0=0.
                     beta=0.
# elif defined SEAMOUNT
                     depth=4500.
                     beta=0.
                     f0=1.E-4
# elif defined SHELFRONT
                     depth=1660.
                     f0=1.E-4
                     beta=0.
# elif defined SOLITON
                     depth=1.
                     f0=0.
                     beta=1.
# elif defined UPWELLING
                     depth=150.
                     f0=-8.26E-5
                     beta=0.
# elif defined RIVER
                     depth=150.
                     f0=8.26E-5
                     beta=0.
# elif defined SEAGRASS
                     depth=0.215
                     f0=0.
                     beta=0.
# elif defined JET
                     depth=4000.
                     f0=1.E-4
                     beta=1.6E-11
# elif defined SHOREFACE
                     depth=15.
                     f0=0.
                     beta=0.
# elif defined SANDBAR
                     depth=4.1
                     f0=0.
                     beta=0.
# elif defined SWASH
                     depth=1.
                     f0=0.
                     beta=0.
# elif defined RIP
                     depth=15.
                     f0=0.
                     beta=0.
# elif defined THACKER
                     depth=10.
                     beta=0.
#  ifdef THACKER_2DV
                     f0=0.
#  else
                     f0=1.E-4
#  endif
# elif defined TANK
                     depth=10.
                     f0=0.
                     beta=0.
# elif defined MOVING_BATHY
                     depth=0.394
                     f0=0.
                     beta=0.
# elif defined DUNE
#  ifdef ANA_DUNE
                     depth=6.
#  else
                     depth=4.
#  endif
                     f0=0.
                     beta=0.
# elif defined SED_TOY
#  ifdef SED_TOY_ROUSE
                     depth=5.
#  elif defined SED_TOY_FLOC_1D
                     depth=5.
#  else
                     depth=20.
#  endif
                     f0=0.0
                     beta=0.
# elif defined TIDAL_FLAT
                     depth_shallow=-4.
                     depth_deep=16.
                     f0=0.
                     beta=0.
# elif defined ESTUARY
                     f0=0.
                     beta=0.
# elif defined KILPATRICK
                     f0=1.E-4
                     beta=0.
# else
                     depth=???
                     f0=???
                     beta=???
# endif
!
!----------------------------------------------------------------------
! Grid dimensions (Length_XI, Length_ETA)
!----------------------------------------------------------------------
!
# define Length_XI_ISO_GRID Length_ETA*float(LLm0)/float(MMm0)
# define Length_ETA_ISO_GRID Length_XI*float(MMm0)/float(LLm0)
!
# ifdef AGRIF
      if (Agrif_Root()) then
# endif
# if defined BASIN
                     Length_XI =3600.0e+3
                     Length_ETA=2800.0e+3
# elif defined SINGLE_COLUMN
                     Length_XI = 10.0e+3
                     Length_ETA= 10.0e+3
# elif defined CANYON
                     Length_XI =128.0e+3
                     Length_ETA= 96.0e+3
# elif defined EQUATOR
                     Length_ETA=3200.0e+3
                     Length_XI =4000.0e+3
# elif defined KH_INST
                     Length_XI =256.0
                     Length_ETA=Length_ETA_ISO_GRID
# elif defined ACOUSTIC
                     Length_XI =128.0
                     Length_ETA=Length_ETA_ISO_GRID
# elif defined GRAV_ADJ
#  ifdef NBQ
                     Length_XI =3.0
#  else
                     Length_XI =64.0e+3
#  endif
                     Length_ETA=Length_ETA_ISO_GRID
# elif defined ISOLITON
                     Length_XI =6.0
                     Length_ETA=Length_ETA_ISO_GRID
# elif defined INNERSHELF
                     Length_XI =200.0e+3
                     Length_ETA=Length_ETA_ISO_GRID
# elif defined INTERNAL
                     Length_XI =1200.0e+3
                     Length_ETA=Length_ETA_ISO_GRID
# elif defined TS_HADV_TEST
                     Length_XI = 1.0e+3
                     Length_ETA= 1.0e+3
# elif defined OVERFLOW
                     Length_ETA=64.0E+3
                     Length_XI =Length_XI_ISO_GRID
# elif defined SEAMOUNT
                     Length_XI =512.0e+3
                     Length_ETA=512.0e+3
# elif defined SHELFRONT
                     Length_ETA=200.0e+3
                     Length_XI =Length_XI_ISO_GRID
# elif defined SOLITON
                     Length_XI =48.
                     Length_ETA=16.
# elif defined UPWELLING
                     Length_ETA=80.e+3
                     Length_XI =Length_XI_ISO_GRID
# elif defined RIVER
                     Length_ETA=80.e+3
                     Length_XI =Length_XI_ISO_GRID
# elif defined SEAGRASS
                     Length_ETA=0.25
                     Length_XI =36.*0.05
# elif defined JET
                     Length_XI = 500.e+3
                     Length_ETA=2000.e+3
# elif defined SHOREFACE
                     Length_XI =1180.0
                     Length_ETA=Length_ETA_ISO_GRID
# elif defined SANDBAR
                     Length_XI =200.
                     Length_ETA=Length_ETA_ISO_GRID
# elif defined SWASH
                     Length_XI =110.0
                     Length_ETA=Length_ETA_ISO_GRID
# elif defined RIP
                     Length_XI =768.0
                     Length_ETA=768.0
# elif defined THACKER
                     Length_XI =200.e+3
                     Length_ETA=Length_ETA_ISO_GRID
# elif defined TANK
#  ifndef TANKY
                     Length_XI =10.0
                     Length_ETA=Length_ETA_ISO_GRID
#  else
                     Length_ETA=10.0
                     Length_XI =Length_XI_ISO_GRID
#  endif
# elif defined MOVING_BATHY
                     Length_XI =4.0
                     Length_ETA=Length_ETA_ISO_GRID
# elif defined DUNE
#  ifdef ANA_DUNE
                     Length_XI = 300.0
#  else
                     Length_XI = 100.0
#  endif
#  ifdef MUSTANG
                     Length_ETA=Length_XI
#  else
                     Length_ETA=Length_ETA_ISO_GRID
#  endif
# elif defined SED_TOY
#  ifdef SED_TOY_ROUSE
                     Length_XI = 10.e+3
                     Length_ETA= 10.e+3
#  else
                     Length_XI =40.0
                     Length_ETA=30.0
#  endif
# elif defined TIDAL_FLAT
                     Length_ETA= 1500.0
                     Length_XI = 100000.0
# elif defined ESTUARY
                     Length_ETA= 9000.0   ! dy =100m
                     Length_XI = 120000.0  ! dx = 600m or 100m if ESTUARY_SQUARE
# elif defined KILPATRICK
                     Length_ETA= 18.E+3
                     Length_XI = 3600.E+3
# else
                     Length_XI =???
                     Length_ETA=???
# endif

# undef Length_XI_ISO_GRID
# undef Length_ETA_ISO_GRID

!----------------------------------------------------------------------
! Copy physical dimensions of the grid into globally vizible variables
!----------------------------------------------------------------------
!
        xl=Length_XI
        el=Length_ETA

# ifdef AGRIF
      else
        xl=Agrif_Parent(xl)/Agrif_Parent(LLm)*LLm/Agrif_Rhox()
        el=Agrif_Parent(el)/Agrif_Parent(MMm)*MMm/Agrif_Rhoy()
        Length_XI=xl
        Length_ETA=el
      endif
# endif /* AGRIF */
!
!
!----------------------------------------------------------------------
! Set grid spacings for rectangular grids
!----------------------------------------------------------------------
!
      dx=Length_XI/float(LLm)
      dy=Length_ETA/float(MMm)
!
!----------------------------------------------------------------------
! Set reference point location
!----------------------------------------------------------------------
!
      x0=0.
      y0=0.
# ifdef EQUATOR
      y0=y0-0.5*Length_ETA
# endif
# if defined INTERNAL || defined THACKER || defined MOVING_BATHY
      x0=x0-0.5*Length_XI
      y0=y0-0.5*Length_ETA
# endif
# ifdef KILPATRICK
      x0=-1800.E+3
      y0=-9.E+3
# endif

CR      write(*,'(4(A3,pe15.9,1x),I3)') 'dx=',dx, 'dy=',dy,
CR     &                                'x0=',x0, 'y0=',y0, mynode
!
!----------------------------------------------------------------------
! Setup rectangular grid:  coordinates (XI,ETA) at PSI- and RHO-points.
!----------------------------------------------------------------------
!
!$acc kernels if(compute_on_device) default(present)
      do j=JstrR,JendR
        do i=IstrR,IendR
         xp(i,j) = x0 + dx*float(i+iminmpi-1-1)
         xr(i,j) = x0 + dx*(float(i+iminmpi-1-1)+0.5)
         yp(i,j) = y0 + dy*float(j+jminmpi-1-1)
         yr(i,j) = y0 + dy*(float(j+jminmpi-1-1)+0.5)
        enddo
      enddo
!$acc end kernels      
!$acc update host( xr, yr )      
!
!----------------------------------------------------------------------
! Compute coordinate transformation metrics at RHO-points "pm" and
! "pn"  (1/m) associated with the differential distances in XI and
! ETA, respectively.
!----------------------------------------------------------------------
!
!$acc kernels if(compute_on_device) default(present)
      do j=JstrR,JendR
        do i=IstrR,IendR
          pm(i,j)=1./dx
          pn(i,j)=1./dy
        enddo
      enddo
!$acc end kernels      
!
!----------------------------------------------------------------------
!  Set Coriolis parameter [1/s] at RHO-points.
!----------------------------------------------------------------------
!
!$acc kernels if(compute_on_device) default(present)
      do j=JstrR,JendR
        do i=IstrR,IendR
# ifdef EQUATOR
          f(i,j)=f0+beta*yr(i,j)
# else
          f(i,j)=f0+beta*(yr(i,j)-(0.5*el))
# endif
        enddo
      enddo
!$acc end kernels      

# if defined UV_COR_NT || defined CROCO_QH
!
!  Horizontal (non-traditional) Coriolis parameter
!
!$acc kernels if(compute_on_device) default(present)
      do j=JstrR,JendR
        do i=IstrR,IendR
          if (beta.eq.0) then     ! f-plan
            e(i,j)=2.*Erotation*cos(asin(f0/(2.*Erotation)))
          else
            e(i,j)=beta*Eradius   ! beta plan
          endif
        enddo
      enddo
!$acc end kernels      
# endif
!
!----------------------------------------------------------------------
!  Set bathymetry [meters; positive] at RHO-points.
!----------------------------------------------------------------------
!
# if defined BASIN || defined EQUATOR || defined GRAV_ADJ \
                   || defined SOLITON || defined ISOLITON \
                   || defined KH_INST || defined TS_HADV_TEST \
                   || defined ACOUSTIC || defined SINGLE_COLUMN \
                   || defined JET
!$acc kernels if(compute_on_device) default(present)
      do j=JstrR,JendR
        do i=IstrR,IendR
          h(i,j)=depth
        enddo
      enddo

# elif defined CANYON
      do j=JstrR,JendR
        do i=IstrR,IendR
          cff=32000.-16000.*(sin(pi*xr(i,j)/Length_XI))**24
          h(i,j)=20.+0.5*(depth-20.)*(1.+tanh((yr(i,j)-cff)/10000.))
        enddo
      enddo

# elif defined INNERSHELF
      do j=JstrR,JendR
        do i=IstrR,IendR
!         --- constant slope ---
          h(i,j)=1.e-3*(xl+0.5*dx-xr(i,j))+4.
!         --- shelf/slope tanh ---
!          h(i,j)=30.+500.*(1+sinh((150.e3-xr(i,j))/40.e3)
!     &                      /cosh((150.e3-xr(i,j))/40.e3))
        enddo
      enddo

# elif defined INTERNAL
      cff1=1./(ridge_width*ridge_width)
      do j=JstrR,JendR
        do i=IstrR,IendR
          if ((xr(i,j).gt.ridge_width).or.
     &       (xr(i,j).lt.(-ridge_width))) then
            cff=0.
          else
            cff=1-(xr(i,j)*xr(i,j)*cff1);
          endif
          h(i,j)=depth-ridge_height*cff*cff;
#  ifdef INTERNALSHELF
          if (xr(i,j).ge.0.) then
            h(i,j)=depth-ridge_height
          endif
#  endif
        enddo
      enddo

# elif defined OVERFLOW
      do j=JstrR,JendR
        do i=IstrR,IendR
          h(i,j)=20.+0.5*(depth-20.)*(1.+
     &                         tanh((yr(i,j)-40000.)/5000.))
        enddo
      enddo

# elif defined SEAMOUNT
      cff=(1./50.e3)**2
      do j=JstrR,JendR
        do i=IstrR,IendR
          h(i,j)=depth*(1-0.6*exp(-cff*((xr(i,j)-0.5*xl)**2+
     &                                  (yr(i,j)-0.5*el)**2)))
        enddo
      enddo

# elif defined SHELFRONT
      do j=JstrR,JendR
        do i=IstrR,IendR
          cff=yr(i,j)/1000.
          if (cff.lt.50.) then
            h(i,j)=50.+2.*cff
          elseif (cff.lt.60.) then
            h(i,j)=160.+1.5*(cff-50.)**2-0.1*(cff-60.0)**2
          elseif (cff.lt.100.) then
            h(i,j)=310.+30.*(cff-60.)
          elseif (cff.lt.110.) then
            h(i,j)=1660.-1.5*(cff-110.)**2
          else
            h(i,j)=1660.
          endif
        enddo
      enddo

# elif defined UPWELLING
      do j=JstrR,JendR
#  ifdef MPI
        y=dy*float(j+jj*Mm)
#  else
        y=dy*float(j)
#  endif
        if (y.gt.Length_ETA/2.) y=Length_ETA-y+dy
        cff=min(depth,84.5+66.526*tanh(0.00015*(y-0.125*Length_ETA)))
        do i=IstrR,IendR
          h(i,j)=cff
        enddo
      enddo

# elif defined RIVER
      do i=IstrR,IendR
#  ifdef MPI
        cff=(float(i +ii*Lm)-0.5)/float(LLm)
#  else
        cff=(float(i       )-0.5)/float(LLm)
#  endif
        if (cff.lt.0.05) then
          h(i,JstrR)=15.
        elseif (cff.lt.0.15) then
          h(i,JstrR)=15.+843.75*(cff-0.05)**2
        elseif (cff.lt.0.85) then
          h(i,JstrR)=15.+168.75*(cff-0.1)
        elseif (cff.lt.0.95) then
          h(i,JstrR)=150.-843.75*(cff-0.95)**2
        else
          h(i,JstrR)=150.
        endif
      enddo
      do j=JstrR+1,JendR
        do i=IstrR,IendR
          h(i,j)=h(i,JstrR)
        enddo
      enddo

# elif defined SEAGRASS
        do j=JstrR,JendR
            do i=IstrR,IendR
            h(i,j)=depth
            enddo
        enddo


# elif defined SHOREFACE
      do j=JstrR,JendR
        do i=IstrR,IendR
#  ifdef MPI
          cff=(float(i +ii*Lm) )
#  else
          cff=(float(i       ) )
#  endif
#  ifdef WET_DRY
          h(i,j)=12.0-0.0125*Length_XI/float(LLm+1)*cff ! original fonction
#  else
          h(i,j)=12.0-0.0125*Length_XI/float(LLm+1)*cff + 3.
#  endif
#  ifdef MASKING
          rmask(i,j)=1.0
          if (h(i,j).lt.0.01) then
            h(i,j)=0.01
            rmask(i,j)=0.0
            rmask(48,j)=1.0
          endif
#  endif
        enddo
      enddo

# elif defined SANDBAR
      do j=JstrR,JendR
        do i=IstrR,IendR
#  ifdef MPI
          cff=(float(i +ii*Lm) )
#  else
          cff=(float(i       ) )
#  endif
#  ifdef SANDBAR_OFFSHORE
! Roelvink & Reniers 1995: LIP-1B
          xx=xr(i,j)
          if (xx.lt.19) then
            h(i,j)=-depth
          elseif (xx.lt.52) then
            p1=+1.520229705740388e-08
            p2=-3.423212989640106e-06
            p3=+3.159005532947043e-04
            p4=-1.526398548035649e-02
            p5=+4.066341504316825e-01
            p6=-5.608106604351117e+00
            p7=+2.702732753680607e+01
            h(i,j) = p1*xx**6 + p2*xx**5 + p3*xx**4 +
     &               p4*xx**3 + p5*xx**2 + p6*xx    + p7
          elseif (xx.lt.120) then
            p1=-4.642289413679605e-11
            p2=+2.870629635648417e-08
            p3=-7.130431368316412e-06
            p4=+9.126154058771732e-04
            p5=-6.357269253701968e-02
            p6=+2.303201222063785e+00
            p7=-3.657085297281925e+01
            h(i,j) = p1*xx**6 + p2*xx**5 + p3*xx**4 +
     &               p4*xx**3 + p5*xx**2 + p6*xx    + p7
          elseif (xx.lt.140) then
            p1=+3.099533433018596e-08
            p2=-2.425187846848144e-05
            p3=+7.897142871067360e-03
            p4=-1.369896736259833e+00
            p5=+1.335157006848632e+02
            p6=-6.932492254089954e+03
            p7=+1.498142545396247e+05
            h(i,j) = p1*xx**6 + p2*xx**5 + p3*xx**4 +
     &               p4*xx**3 + p5*xx**2 + p6*xx    + p7
          elseif (xx.lt.160) then
            p1=-2.524608914887271e-08
            p2=+2.387297325503284e-05
            p3=-9.398023259277126e-03
            p4=+1.971403061349242e+00
            p5=-2.323952920877027e+02
            p6=+1.459649630832179e+04
            p7=-3.816025205750951e+05
            h(i,j) = p1*xx**6 + p2*xx**5 + p3*xx**4 +
     &               p4*xx**3 + p5*xx**2 + p6*xx    + p7
          else
            p1=+3.693458909496914e-09
            p2=-4.051059061933883e-06
            p3=+1.847329725370577e-03
            p4=-4.483034673182216e-01
            p5=+6.106360892732523e+01
            p6=-4.426507173590056e+03
            p7=+1.334155898180324e+05
            h(i,j) = p1*xx**6 + p2*xx**5 + p3*xx**4 +
     &               p4*xx**3 + p5*xx**2 + p6*xx    + p7
          endif
#  elif defined SANDBAR_ONSHORE
! Roelvink & Reniers 1995: LIP-1C
          xx=xr(i,j)
          if (xx.lt.19) then
            h(i,j)=-depth
          elseif (xx.lt.52) then
            p1=-1.550500196449383e-09
            p2=+3.095246505454406e-07
            p3=-2.524513106364847e-05
            p4=+1.072358697579965e-03
            p5=-2.475915084067423e-02
            p6=+3.341039769921187e-01
            p7=-6.243997999730277e+00
            h(i,j) = p1*xx**6 + p2*xx**5 + p3*xx**4 +
     &               p4*xx**3 + p5*xx**2 + p6*xx    + p7
          elseif (xx.lt.120) then
            p1=+2.879040803660182e-11
            p2=-3.112995372255505e-09
            p3=-2.128994193906542e-06
            p4=+5.806818410041643e-04
            p5=-5.889240225358890e-02
            p6=+2.719179108579184e+00
            p7=-5.021582880671649e+01
            h(i,j) = p1*xx**6 + p2*xx**5 + p3*xx**4 +
     &               p4*xx**3 + p5*xx**2 + p6*xx    + p7
          elseif (xx.lt.140) then
            p1=+6.722443006700593e-07
            p2=-5.191501807411085e-04
            p3=+1.669215462536683e-01
            p4=-2.860227732993812e+01
            p5=+2.754780791103959e+03
            p6=-1.414013613750512e+05
            p7=+3.022008907301721e+06
            h(i,j) = p1*xx**6 + p2*xx**5 + p3*xx**4 +
     &               p4*xx**3 + p5*xx**2 + p6*xx    + p7
          elseif (xx.lt.160) then
            p1=-3.300743083271716e-07
            p2=+2.941386493481447e-04
            p3=-1.091628962526408e-01
            p4=+2.159706871907915e+01
            p5=-2.402365333342808e+03
            p6=+1.424582143575532e+05
            p7=-3.518304916201898e+06
            h(i,j) = p1*xx**6 + p2*xx**5 + p3*xx**4 +
     &               p4*xx**3 + p5*xx**2 + p6*xx    + p7
          else
            p1=+3.341378689344706e-09
            p2=-3.683115456141188e-06
            p3=+1.686167358885289e-03
            p4=-4.104224500606568e-01
            p5=+5.602413257412504e+01
            p6=-4.066824640330583e+03
            p7=+1.226609757884486e+05
            h(i,j) = p1*xx**6 + p2*xx**5 + p3*xx**4 +
     &               p4*xx**3 + p5*xx**2 + p6*xx    + p7
          endif
#  endif
          h(i,j) = min(-h(i,j),depth)
        enddo
      enddo

# elif defined SWASH
      do j=JstrR,JendR
        do i=IstrR,IendR
#  ifdef MPI
          cff=(float(i +ii*Lm) )
#  else
          cff=(float(i       ) )
#  endif
          h(i,j)=1.0571-0.0125*Length_XI/float(LLm+1)*cff  ! GLOBEX profil
#  ifndef WET_DRY
     &                                              + 0.6
#  endif
          if (xr(i,j).lt.16.57) h(i,j)=0.85               ! GLOBEX
        enddo
      enddo

# elif defined RIP
#  ifdef GRANDPOPO
!
!  Idealized Grand Popo Beach in Benin
!
      xs=85;        ! inner surf zone
      db=50;        ! distance from xs to sand bar
      alpha=0.01;
      do j=JstrR,JendR
        do i=IstrR,IendR
#   ifdef OBC_EAST
          xx=xr(i,j)
#   else
          xx=Length_XI-xr(i,j)
#   endif
          h(i,j)=-4.5
     &       -1.7*exp(-6*(((xx-xs-db)/db)**2))
     &       +3.1*(1+tanh(0.025*(xx-xs)))
     &       +0.014*(xx+log(cosh(alpha*(xx-xs))/cosh(alpha*xs))/alpha)
        enddo
      enddo
#  else
!
! Weir et al. 2011 after Lippmann et al. 1999 (idealization of Duck beach)
!
      xs=150      ! inner surf zone
      db=80       ! distance from xs to sand bar
      alpha=0.02
      lambda=256  ! sand bar perturbation wavelength
#   ifdef RIP_TOPO_2D
      eps=0.0     ! no sand bar perturbation
#   else
      eps=0.1     ! sand bar perturbation magnitude
#   endif
      do j=JstrR,JendR
        do i=IstrR,IendR
#   ifdef OBC_EAST
          xx=xr(i,j)
#   else
          xx=Length_XI-xr(i,j)
#   endif
          yy=yr(i,j);
          h0=-1.5*exp(-5*(((xx-xs-db)/db)**2))
     &       +1.35*(1+tanh(0.0253*(xx-xs)))
     &       +0.0032*(xx+log(cosh(alpha*(xx-xs))/cosh(alpha*xs))/alpha)
          per=eps*cos(2*pi*yy/lambda)*exp(-5*(((xx-xs-db)/db)**2))
          h(i,j)=(1+per)*h0!-0.8;
        enddo
      enddo
#  endif

# elif defined THACKER
      cff=80.e3
      do j=JstrR,JendR
        do i=IstrR,IendR
#  ifdef THACKER_2DV
          h(i,j)=depth*(1-(xr(i,j)/cff)**2)
#  else
          h(i,j)=depth*(1-(xr(i,j)/cff)**2 -
     &                    (yr(i,j)/cff)**2 )
#  endif
        enddo
      enddo

# elif defined TANK
      do j=JstrR,JendR
        do i=IstrR,IendR
          h(i,j)=depth
        enddo
      enddo
# elif defined MOVING_BATHY
      cff=1./(0.03686)**2
      do j=JstrR,JendR
        do i=IstrR,IendR
          h(i,j)=depth-0.1*exp(-cff*xr(i,j)**2)
        enddo
      enddo

# elif defined DUNE
#  ifdef ANA_DUNE
      cff2=0.01
#  elif defined DUNE3D
      cff1=0.5*(1./10)**2
#  endif
      do j=JstrR,JendR
        do i=IstrR,IendR
#  ifdef MPI
          cff=(float(i +ii*Lm) )
#  else
          cff=(float(i       ) )
#  endif
#  ifdef ANA_DUNE
          h(i,j)=depth-2.*exp(-cff2*(xr(i,j)-0.5*xl)**2)
#  elif defined DUNE3D
#   ifdef MPI
          cffi=(float(i +ii*Lm) )
          cffj=(float(j +jj*Mm) )
#   else
          cffi=(float(i       ) )
          cffj=(float(j       ) )
#   endif
         if (xr(i,j).lt.4.) then
          h(i,j)=depth
         else
          Rij = sqrt( (cffi-28.)**2.+
     &                (cffj-25.)**2. )
          h(i,j)=depth*(1.-0.5*exp(-cff1*(2.*(Rij))**2.0))
         endif
#  else
         h(i,j)=depth*(1.-0.5*exp(-5.e-3*(2.*(cff-25.))**2))
#  endif
        enddo
      enddo

# elif defined SED_TOY
      do j=JstrR,JendR
        do i=IstrR,IendR
          h(i,j)=depth
!          rmask(i,j)=1.0
        enddo
      enddo

# elif defined TIDAL_FLAT
      do j=JstrR,JendR
        do i=IstrR,IendR
#  ifdef MPI
          cff=(float(i +ii*Lm) )
#  else
          cff=(float(i       ) )
#  endif
          if (cff.lt.2) then
            h(i,j)=depth_deep
          else
            h(i,j)=depth_deep-(cff-1)/float(LLm+1)
     &                        *(depth_deep-depth_shallow)
          endif
        enddo
      enddo

# elif defined ESTUARY

#  if defined ESTUARY_SQUARE
      imouth = 240
#  else
      imouth = 40
#  endif
      jmiddle = int(MMm/2)
      m = 2.0 ! half tide height (m)
      estCWM = 200.0 ! estuary Channel Widht at Mouth (m)
      estCWM_para = 1.0/36000.0  ! estuary Channel Widht at Mouth exponentiel parameter
      estHupstream = 2.5 ! estuary height upstream (m)
      estHmouth = 10.0 ! estuary height at mouth (m)
      estHoff = 20.0 ! estuary height offshore (m)
      estMS = 10000.0 ! estuary Mouth Section width (m)
      estMS_para = 1.0/24000.0 ! estuary Mouth Section width exponentiel parameter
      estMS_para2 = 1.0/3000.0 ! estuary Off Section width exponentiel parameter
      estCS = (estHupstream - estHmouth) / (dx * float(LLm - imouth)) ! estuary channel slope in estuary
      estCS2 = (estHmouth - estHoff) / (dx * float(imouth)) ! estuary channel slope offshore
      estGamma2 = 1.1**2

      do i = IstrR, IendR
#  ifdef MPI
        cff  = (float(i + iminmpi - 1) )
#  else
        cff  = (float(i       ) )
#  endif

        xboti = dx * (cff - imouth)

        if (cff .ge. imouth) then
          l1i = estCWM * exp( - estCWM_para * xboti)
          hl1i = estHmouth + estCS * xboti
          hl1mi = hl1i - m
          shti = estMS * exp( - estMS_para * xboti)! section half tide
          l2i = 100.0
          do itmp = 1, 1000
            l2n = (shti - l1i * hl1mi) / (hl1mi + 1.5 * m +
     &             0.5 * m * (1.0 + 4.0 * (m / hl1mi) * l2i /
     &            (l1i + l2i))**2 / estGamma2)
            l2i = l2n
          enddo
          l2i = l2n
          if (l2i .ge. l1i) then
            l3i = l2i * (1.0 + 4.0 * (m / hl1mi) * l2i /
     &         (l1i + l2i ))**2 / estGamma2
          else
            l2i = l1i
            l3i = (shti - hl1mi * (l1i + l2i) - 2.0 *
     &              m * l2i) * 2.0 / m - l2i
          endif
        else  ! cff < imouth

          l1imouth = estCWM
          hl1imouth = estHmouth
          hl1mimouth = hl1imouth - m
          shtimouth = estMS ! section half tide at imouth
          l2i = 100.0
          do itmp = 1, 1000
            l2n = (shtimouth - l1imouth * hl1mimouth) / (hl1mimouth +
     &            1.5 * m + 0.5 * m * (1.0 + 4.0 * (m / hl1mimouth) *
     &            l2i / (l1imouth + l2i) )**2 / estGamma2)
            l2i = l2n
          enddo
          l2imouth = l2n
          if (l2imouth .ge. l1imouth) then
            l3imouth = l2imouth * (1.0 + 4.0 * (m / hl1mimouth) *
     &            l2imouth / (l1imouth + l2imouth ))**2 / estGamma2
          else
            l2imouth = l1imouth
            l3imouth = (shtimouth - hl1mimouth * (l1imouth + l2imouth) -
     &             2.0 * m * l2imouth) * 2.0 / m - l2imouth
          endif

          l2ml1m = l2imouth - l1imouth
          l3ml1m = l3imouth - l1imouth

          hl1i = estHmouth + estCS2 * xboti
          hl1mi = hl1i - m
          shti = estMS * exp( - estMS_para2 * xboti)! section half tide
          l1i = (shti - l2ml1m * (hl1mi + 1.5 * m ) -
     &            0.5 * l3ml1m * m) / (2.0 * hl1mi + 2.0 * m)
          l2i = l1i + l2ml1m
          l3i = l1i + l3ml1m
        endif

        do j = JstrR, JendR
#  ifdef MPI
          cff2 = (float(j + jminmpi - 1) )
#  else
          cff2 = (float(j       ) )
#  endif
          ybotj = dy * abs(cff2 - jmiddle)

          if (ybotj .le. l1i) then
            h(i,j) = hl1i
          else if (ybotj .le. l2i) then
            h(i,j) = hl1i + (m - hl1i) / (l2i - l1i)
     &             * (ybotj - l1i)
          else
            h(i,j) = amax1(-7.0, m - 2.0 * m / (l3i - l2i)
     &             * (ybotj-l2i))
          endif

          ! walls at limit East
          if (cff .ge. LLm) then
            h(i, j) = -7.0
            if (cff2 .eq. jmiddle .and. cff .eq. LLm) then
                h(i, j) = estHupstream
            endif
          endif

        enddo  ! j


      enddo ! i

# elif defined KILPATRICK
      do j=JstrR,JendR
        do i=IstrR,IendR
          h(i,j)=500.0
        enddo
      enddo
# else
      do j=JstrR,JendR
        do i=IstrR,IendR
          h(i,j)=???
        enddo
      enddo
# endif

!----------------------------------------------------------------------
! initialisation of z0b in case of analytical cases
!----------------------------------------------------------------------
      do j=JstrR,JendR
        do i=IstrR,IendR
          zob(i,j)=zobt
        end do
      end do

# ifdef SEAGRASS
      do j=JstrR,JendR
         do i=IstrR,IendR
#   ifdef MPI
        cff  = (float(i +ii*Lm) )
#   else
        cff  = (float(i       ) )
#   endif
             if ((cff .lt. 14) .or. (cff .gt. 31)) then
                 zob(i,j) = 6.E-8 ! out of patch
             else
                 zob(i,j) = 1.E-6 ! in patch
             end if
         end do
       end do
# endif
!
!----------------------------------------------------------------------
!  Set masking at RHO-points.
!----------------------------------------------------------------------
!
# ifdef MASKING
      do j=JstrR,JendR             ! Set mask to all-water status
        do i=IstrR,IendR
          rmask(i,j)=1.
        enddo
      enddo
#  ifdef ESTUARY
      do i = IstrR, IendR
#   ifdef MPI
        cff  = (float(i + iminmpi - 1) )
#   else
        cff  = (float(i       ) )
#   endif
        ! limit East
        if (cff .ge. LLm) then
            rmask(i,:) = 0.
        endif
      end do

#  endif /* ESTUARY */
#  ifdef RIVER
!
! Mask out 8-point wide strip of land on the west, and
! carve two 1-point wide channels through that strip.

      if (WESTERN_EDGE) then
        do j=JstrR,JendR
          do i=IstrR,8
            rmask(i,j)=0.            ! <-- strip of land
          enddo
        enddo
        do is=1,Nsrc
#   ifdef MPI
          i0=Isrc_mpi(is,mynode)
          j0=Jsrc_mpi(is,mynode)
#   else
          i0=Isrc(is)
          j0=Jsrc(is)
#   endif
          if (IstrR.le.i0 .and. i0.le.IendR .and.
     &        JstrR.le.j0 .and. j0.le.JendR) then
            if (is.eq.1) then
              do j=j0,j0+6
                rmask(3,j)=1.      ! <-- upper channel with corner
              enddo
              do i=IstrR+3,8
                rmask(i,j0+6)=1.
              enddo
            else
              do i=IstrR+3,8
                rmask(i,j0)=1.     ! <-- lower channel along xi
              enddo
            endif
          endif
        enddo
      endif ! WESTERN_EDGE
#  endif /* RIVER */

#  ifdef SEAGRASS
      do j=JstrR,JendR
         do i=IstrR,IendR
         rmask(i,j)=0.      
         enddo
      enddo

      if (WESTERN_EDGE) then
         do i=IstrR+3,IendR
            rmask(i,3)=1.
         enddo
      else
         do i=IstrR,IendR
            rmask(i,3)=1.   
         enddo
      endif ! WESTERN_EDGE

#  endif /* SEAGRASS */

# endif /* MASKING */
!$acc end kernels
!
#else
      subroutine ana_grid_empty
#endif /* ANA_GRID */
      return
      end

