!======================================================================
! 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"

      subroutine ana_initial (tile)
      implicit none
#include "param.h"
       integer tile, trd
#include "compute_tile_bounds.h"
      call ana_initial_tile (Istr,Iend,Jstr,Jend)
      return
      end
!
      subroutine ana_initial_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  Set initial conditions for momentum and tracer variables using
!  analytical expressions.
!----------------------------------------------------------------------
!
#ifdef AGRIF
      use Agrif_Util
#endif
#ifdef SUBSTANCE
      USE comsubstance, ONLY: cini_wat
#endif
      implicit none
#include "param.h"
#include "grid.h"
#include "ocean2d.h"
#include "ocean3d.h"
#include "scalars.h"
#ifdef SEDIMENT
# include "sediment.h"
#endif
#ifdef SST_SKIN
# include "forces.h"
#endif
      integer Istr,Iend,Jstr,Jend, i,j,k, itrc
#if defined SUBSTANCE || defined SEDIMENT
      integer :: iv
#endif
#ifdef SOLITON
      real fac, x,y, x0,y0, cff1,cff2,cff3
#endif
#ifdef BASIN
      real cff1,cff2
#endif
#ifdef INTERNAL
      real cff
#endif
#ifdef IGW
      real cff,sig,rhom,depth
#endif
#ifdef BIOLOGY
      real temp, SiO4
#endif
#ifdef THACKER
      real ETA0, omega, depth, Lt
#endif
#ifdef TANK
      real ETA0, depth, Lt,cff
#endif
#ifdef MOVING_BATHY
      real cff
#endif
#ifdef RIP
      integer iper
      real xper,yper,xx,yy,xs,db,eps
#endif
#ifdef KH_INST
      real rho1_exp,rho2_exp,h1_exp, Lt
      real d(0:N)
!$acc declare create(d)
      real zu0,zr0,zint,depth,dz,gam,del,hint,du,drho
      real hintu
      real zdye10,zdye20
      real hu0,hup0,htemp0,hsalt0,eps3D
#endif
#ifdef ISOLITON
      real zint,depth,dz,gam,del
      real rho1_exp,rho2_exp,h1_exp,amp_exp, Lt
      real d(0:N)
#endif
#ifdef TS_HADV_TEST
      real cff,cff1,cx,cy,cmax,xx,yy
# if defined SOLID_BODY_ROT || defined SOLID_BODY_PER
#  include "work.h"
# endif
#endif
#ifdef ESTUARY
      real cff
#endif
#ifdef DUNE
      real cff, cff2
# ifdef DUNE3D
      real Rij,cffi,cffj
# endif
#endif
#ifdef JET
      integer istr1,iend1,jstr1,jend1
#endif
!
#include "compute_auxiliary_bounds.h"
!
#if defined EW_PERIODIC && !defined MPI
# define IR_RANGE Istr,Iend
# define IU_RANGE Istr,Iend
#else
# define IR_RANGE IstrR,IendR
# define IU_RANGE Istr,IendR
#endif
#if defined NS_PERIODIC && !defined MPI
# define JR_RANGE Jstr,Jend
# define JV_RANGE Jstr,Jend
#else
# define JR_RANGE JstrR,JendR
# define JV_RANGE Jstr,JendR
#endif

#ifdef ANA_INITIAL

!$acc kernels if(compute_on_device) default(present)
!
!----------------------------------------------------------------------
!  Initial conditions for free surface and 2D momentum components.
!----------------------------------------------------------------------
!
      if (nrrec.eq.0) then
!
# if defined SOLITON
#  ifdef AGRIF
        if (Agrif_Root()) then
#  endif
          x0=2.*xl/3.
          y0=el/2.
#  ifdef AGRIF
        else
          x0=2.*Agrif_Parent_xl/3.
          y0=Agrif_Parent_el/2.
        endif
#  endif
        cff1=0.395
        cff2=0.771*(cff1*cff1)
        do j=JR_RANGE
          do i=IR_RANGE
            x=xr(i,j)-x0
            y=yr(i,j)-y0
            cff3=exp(-cff1*x)
            fac=cff2*(2.*cff3/(1.+cff3*cff3))**2
            zeta(i,j,1)=0.25*fac*(6.*y*y+3.)*exp(-0.5*y*y)
          enddo
        enddo
        do j=JR_RANGE
          do i=IU_RANGE
            x=0.5*(xr(i-1,j)+xr(i,j))-x0
            y=0.5*(yr(i-1,j)+yr(i,j))-y0
            cff3=exp(-cff1*x)
            fac=cff2 * (2.*cff3/(1.+cff3*cff3))**2
            ubar(i,j,1)=0.25*fac*(6.*y*y-9.)*exp(-0.5*y*y)
          enddo
        enddo
        do j=JV_RANGE
          do i=IR_RANGE
            x=0.5*(xr(i,j-1)+xr(i,j))-x0
            y=0.5*(yr(i,j-1)+yr(i,j))-y0
            cff3=exp(-cff1*x)
            fac=cff2 * (2.*cff3/(1.+cff3*cff3))**2
            vbar(i,j,1)=2.*fac*y*(-2.*cff1*tanh(cff1*x))
     &                                    *exp(-0.5*y*y)
          enddo
        enddo

# elif defined KH_INST
        do j=JR_RANGE
          do i=IR_RANGE
            zeta(i,j,1)=0.
            zeta(i,j,2)=0.
          enddo
        enddo

# elif defined THACKER
        ETA0=0.1
        depth=10.
        Lt=80.e3
        omega=0.5*f(1,1) + sqrt( 0.25*f(1,1)**2 + 2*g*depth/Lt**2 )
        do j=JR_RANGE
          do i=IR_RANGE
            zeta(i,j,1)=2*ETA0*depth*(xr(i,j)/Lt - .5*ETA0/Lt)
            zeta(i,j,2)=zeta(i,j,1)
            ubar(i,j,1)=0.
            vbar(i,j,1)=-ETA0*omega*Lt
            ubar(i,j,2)=0
            vbar(i,j,2)=vbar(i,j,1)
          enddo
        enddo

# elif defined TANK
        ETA0=0.001
        Lt=10.
        do j=JR_RANGE
          do i=IR_RANGE
#  ifndef TANKY
            zeta(i,j,1)=ETA0*cos(pi*xr(i,j)/Lt)
#  else
            zeta(i,j,1)=ETA0*cos(pi*yr(i,j)/Lt)
#  endif
            ubar(i,j,1)=0.
            vbar(i,j,1)=0.
            zeta(i,j,2)=zeta(i,j,1)
            ubar(i,j,2)=ubar(i,j,1)
            vbar(i,j,2)=vbar(i,j,1)
          enddo
        enddo

# elif defined DUNE
        cff=0.5*(1./10)**2
#  ifdef ANA_DUNE
        cff2=1.67
#  else
        cff2=0.58
#  endif
        do j=JR_RANGE
          do i=IR_RANGE
#  if defined DUNE3D
#   ifdef MPI
              cffi=(float(i +ii*Lm))
              cffj=(float(j +jj*Mm))
#   else
              cffi=(float(i))
              cffj=(float(j))
#   endif
              Rij = sqrt( (cffi-28.)**2.+
     &                    (cffj-25.)**2. )
              ubar(i,j,1)=cff2/(1-0.5*exp(-cff*(2.*(Rij))**2.))
#  else
              ubar(i,j,1) = cff2/(1-
     &               0.5*exp(-cff*(xp(i,j)-0.5*xl)**2))
#  endif
              ubar(i,j,:) = ubar(i,j,1)
              vbar(i,j,1) = 0.
              vbar(i,j,:) = vbar(i,j,1)
              zeta(i,j,1)=0.0464+xr(i,j)*(0.0394-0.0464)/xl
              zeta(i,j,:) = zeta(i,j,1)
          enddo
        enddo

# elif defined SED_TOY && ( defined SED_TOY_RESUSP || \
                            defined SED_TOY_CONSOLID )
        do j=JR_RANGE
          do i=IR_RANGE
            zeta(i,j,1)=0.
            zeta(i,j,2)=0.
            ubar(i,j,1)=0.
            vbar(i,j,1)=0.
            ubar(i,j,2)=ubar(i,j,1)
            vbar(i,j,2)=0.
          enddo
        enddo

# else
        do j=JR_RANGE
          do i=IR_RANGE
            zeta(i,j,1)=0.
            zeta(i,j,2)=0.
            ubar(i,j,1)=0.
            vbar(i,j,1)=0.
            ubar(i,j,2)=0.
            vbar(i,j,2)=0.
          enddo
        enddo
# endif

      else  !  nrrec.ne.0.

# if defined RIP && defined RIP_TOPO_2D
!
!  Add perturbation to restart field
!
#  ifdef GRANDPOPO
        xs=85        ! inner surfzone at GPP
        db=50        ! distance from xs to sand bar
#  else
        xs=150       ! inner surfzone at Duck
        db=80        ! distance from xs to sand bar
#  endif
        eps=0.01     ! perturbation magnitude
        do j=JstrR,JendR
          do i=IstrR,IendR
#  ifdef OBC_EAST
            xx=xr(i,j)
#  else
            xx=xl-xr(i,j)
#  endif
            yy=yr(i,j);
            xper=1.
            yper=0.
            do iper=1,12
              yper= yper+eps*cos(2*pi*float(iper)*yy/(xs+db)
     &                           +2*pi*tanh(float(iper)/10))
            enddo
            ubar(i,j,1)=(1+xper)*ubar(i,j,1);
            ubar(i,j,2)=ubar(i,j,1);
            vbar(i,j,1)=(1+yper)*vbar(i,j,1);
            vbar(i,j,2)=vbar(i,j,1);
          enddo
        enddo
# endif /* RIP */

      endif  !  nrrec.eq.0.

# ifdef SOLVE3D
!
!----------------------------------------------------------------------
!  Initial conditions for momentum components [m/s].
!----------------------------------------------------------------------
!
      if (nrrec.eq.0) then
!
#  if defined THACKER
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              u(i,j,k,1)=0.
              v(i,j,k,1)=-ETA0*omega*Lt
              u(i,j,k,2)=0.
              v(i,j,k,2)=v(i,j,k,1)
            enddo
          enddo
        enddo
#  elif defined TS_HADV_TEST
#   ifdef DIAGONAL_ADV
        do j=Jstr,Jend
          do k=1,N
            do i=IstrU,Iend
              u(i,j,k,1)  = 1.D0
              u(i,j,k,2)  = u(i,j,k,1)
            enddo
            do i=Istr,Iend
              if (j.ge.JstrV) then
                v(i,j,k,1) = 1.D0
                v(i,j,k,2) = v(i,j,k,1)
              end if
            enddo
          enddo
        enddo
#   endif
#   if defined SOLID_BODY_ROT || defined SOLID_BODY_PER

#    define psi work2d
        cff  = - xl / pi
        cff1 =   pi / xl
        do j=Jstr,Jend+1
          do i=Istr,Iend+1
            psi(i,j) =  cff * (sin(cff1*xp(i,j)))**2
     &                      * (sin(cff1*yp(i,j)))**2
          enddo
        enddo

#    if defined EW_PERIODIC || defined NS_PERIODIC || defined  MPI
        call exchange_p2d_tile(Istr,Iend,Jstr,Jend,psi)
#    endif

        cff  = pm(1,1)
        do k=1,N
          do j=Jstr,Jend
            do i=IstrU,Iend
              u(i,j,k,1)  = - cff * (psi(i,j+1)-psi(i,j))
              u(i,j,k,2)  =  u(i,j,k,1)
#    if defined SOLID_BODY_PER
              u(i,j,k,3)  =  u(i,j,k,1)
#    endif
            enddo
          enddo
        enddo
        do k=1,N
          do j=JstrV,Jend
            do i=Istr,Iend
              v(i,j,k,1) =    cff * (psi(i+1,j)-psi(i,j))
              v(i,j,k,2) = v(i,j,k,1)
#    if defined SOLID_BODY_PER
              v(i,j,k,3) = v(i,j,k,1)
#    endif
            enddo
          enddo
        enddo
#   endif
#   undef psi
       cmax = 0.D0
       do j=JstrV,Jend-1
         do i=IstrU,Iend-1
            cx = dt*( MAX(u(i+1,j,N,1),0.)-MIN(u(i,j,N,1),0.) )*pm(i,j)
            cy = dt*( MAX(v(i,j+1,N,1),0.)-MIN(v(i,j,N,1),0.) )*pn(i,j)
            cmax = MAX( cx+cy, cmax )
         enddo
       enddo
       print*,'CFL_max = ',cmax

#  elif defined SED_TOY && ( defined SED_TOY_RESUSP || \
                             defined SED_TOY_CONSOLID )
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              u(i,j,k,1)=0.
              u(i,j,k,2)=u(i,j,k,1)
              v(i,j,k,1)=0.
              v(i,j,k,2)=v(i,j,k,1)
            enddo
          enddo
        enddo

#  elif defined KH_INST
!
! Traditional KH instability:
!
!  Configuration from Pelletier et al.: linear shear within mixed layer
!
        du     = 1.0    ! Value to be replaced by the RunBuilder script.
        hu0    = 10.
        eps3D  = 0.   ! 3D perturbation
!
!$acc loop private (hint,hintu)
        do j=JR_RANGE
          do i=IU_RANGE
            depth=h(i,j)
            zu0=-0.5*depth   ! u & rho shears at same depth (zu0=zr0)
!           zu0=-0.4*depth   ! shifted profile for thermocline entrainment studies
            hup0=0.15*depth  ! perturbation
            hint = 0.
            hintu = 0.
            do k=1,N
              u(i,j,k,1)  = -du*tanh((z_r(1,1,k)-zu0)/hu0)
              u(i,j,k,1)  = u(i,j,k,1) + 0.01*du *               ! add perturbation
     &                       (-2.*tanh((z_r(1,1,k)-zu0)/hup0))*
     &                       (1. -tanh((z_r(1,1,k)-zu0)/hup0)**2.)
     &                                      *sin(2.*2.*pi/xl*xp(i,j))
     &                       *(1.+eps3D*(sin( 2.*4.*pi/el*yr(i,j))))

              u(i,j,k,2)  = u(i,j,k,1)
              hintu = hintu + 0.5*(Hz(i,j,k)+Hz(i-1,j,k))
     &                          * u(i,j,k,1)
              hint = hint   + 0.5*(Hz(i,j,k)+Hz(i-1,j,k))
            enddo
            hintu = hintu / hint
            do k=1,N
                u(i,j,k,1) = u(i,j,k,1) - hintu
                u(i,j,k,2) = u(i,j,k,1)
            enddo 
            ubar(i,j,:) = 0.
          enddo
        enddo

        do j=JR_RANGE
          do i=IR_RANGE
            depth=h(i,j)
            zu0=-0.5*depth   ! u & rho shears at same depth (zu0=zr0)
!           zu0=-0.4*depth   ! shifted profile for thermocline entrainment studies
            hup0=0.15*depth  ! perturbation
            do k=1,N
              v(i,j,k,1)  = 0.
              v(i,j,k,2)  = 0.
#   ifdef NBQ
              wz(i,j,k,1) = -0.01*du*hup0*
     &                              (1.-tanh((z_w(1,1,k)-zu0)/hup0)**2)
     &                                  *2.*2.*pi/xl*cos(2.*2.*pi/xl*xp(i,j))
     &                                  *(1.+eps3D*(sin( 2.*4.*pi/el*yp(i,j))))

              wz(i,j,k,2) = wz(i,j,k,1)
#   endif
            enddo
          enddo
        enddo
#   ifdef NBQ
        do j=JR_RANGE
          do i=IR_RANGE
              wz(i,j,0,:)=0.
          enddo
        enddo
#   endif

#  else
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              u(i,j,k,1)=0.
              v(i,j,k,1)=0.
              u(i,j,k,2)=0.
              v(i,j,k,2)=0.
#   ifdef NBQ
              wz(i,j,k,:)=0.
#   endif
            enddo
          enddo
        enddo
#   ifdef NBQ
        do j=JR_RANGE
          do i=IR_RANGE
             wz(i,j,0,1)=0.
             wz(i,j,0,2)=0.
          enddo
        enddo
#   endif

#  endif

      else   ! nrrec.ne.0.

#  if defined RIP && defined RIP_TOPO_2D
!
!  Add perturbation to restart field
!
        do k=1,N
          do j=JstrR,JendR
            do i=IstrR,IendR
#   ifdef OBC_EAST
              xx=xr(i,j)
#   else
              xx=xl-xr(i,j)
#   endif
              yy=yr(i,j);
              xper=1.
              yper=0.
              do iper=1,12
                yper= yper+eps*cos(2*pi*float(iper)*yy/(xs+db)
     &                            +2*pi*tanh(float(iper)/10))
              enddo
              u(i,j,k,1)=(1+xper)*u(i,j,k,1)
              u(i,j,k,2)=u(i,j,k,1)
              v(i,j,k,1)=(1+yper)*v(i,j,k,1)
              v(i,j,k,2)=v(i,j,k,1)
            enddo
          enddo
        enddo
#  endif /* RIP */
!
      endif   !nrrec.eq.0.
!
!----------------------------------------------------------------------
!  Set initial conditions for potential temperature [degC] and
!  salinity [PSU] and passive tracer
!----------------------------------------------------------------------
!
      if (nrrec.eq.0) then
!
#  ifdef TRACERS
#   ifdef BASIN
        cff1=(44.690/39.382)**2
        cff2=cff1*(rho0*800./g)*(5.0e-5/((42.689/44.690)**2))
#    ifdef TEMPERATURE
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              t(i,j,k,1,itemp)=cff2*exp(z_r(i,j,k)/800.)
     &                   *(0.6-0.4*tanh(z_r(i,j,k)/800.))
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
             enddo
          enddo
        enddo
#    endif /* TEMPERATURE */
#   elif defined SINGLE_COLUMN && defined WILLIS_DEARDORFF
#    ifdef TEMPERATURE
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              t(i,j,k,1,itemp)=0.1 * z_r(i,j,k) + 22.
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
             enddo
          enddo
        enddo
#    endif /* TEMPERATURE */
#   elif defined SINGLE_COLUMN && defined DIURNAL_CYCLE
#    ifdef TEMPERATURE
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              t(i,j,k,1,itemp)=0.01 * z_r(i,j,k) + 22.
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
             enddo
          enddo
        enddo
#    endif /* TEMPERATURE */
#   elif defined SINGLE_COLUMN && defined KATO_PHILIPS
#    ifdef TEMPERATURE
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              t(i,j,k,1,itemp)=16.0+(0.01)**2/(g*Tcoef/rho0)
     &             * (z_r(i,j,k))
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
             enddo
          enddo
        enddo
#    endif /* TEMPERATURE */
#   elif defined TS_HADV_TEST && defined DIAGONAL_ADV
        cff = 1./(50.*50.)
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=exp( -cff * ( (xr(i,j)-0.25*xl)**2
     &                                   +   (yr(i,j)-0.25*el)**2 ) )
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif /* TEMPERATURE */
#    ifdef SALINITY
              t(i,j,k,1,isalt)=1.0    !<-- test constancy preservation
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo
#   elif defined TS_HADV_TEST && ( defined SOLID_BODY_ROT || defined SOLID_BODY_PER )
        cff = 1./(50.*50.)
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=exp( -cff * ( (xr(i,j)-0.25*xl)**2
     &                                   +   (yr(i,j)-0.25*el)**2 ) )
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif /* TEMPERATURE */
#    ifdef SALINITY
              t(i,j,k,1,isalt)=1.0    !<-- test constancy preservation
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo
#   elif defined CANYON
#    ifdef CANYON_STRAT
#     ifdef TEMPERATURE
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              t(i,j,k,1,itemp)=3.488*exp(z_r(i,j,k)/800.)
     &          *(1.-0.666666666666*tanh(z_r(i,j,k)/800.))
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
            enddo
          enddo
        enddo
#     endif /* TEMPERATURE */
#    else
#     ifdef TEMPERATURE
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              t(i,j,k,1,itemp)=T0
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
            enddo
          enddo
        enddo
#     endif /* TEMPERATURE */
#    endif
#   elif defined SHOREFACE || defined RIP
        do k=1,N
          do j=jstrR,jendR
            do i=istrR,iendR
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=18.0
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif /* TEMPERATURE */
#    ifdef SALINITY
              t(i,j,k,1,isalt)=32.0
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo

#   elif defined SWASH
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=20.
!              t(i,j,k,1,itemp)=20.-(zeta(i,j,1)-z_r(i,j,k))
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif /* TEMPERATURE */
            enddo
          enddo
        enddo

#   elif defined THACKER
#    ifdef TEMPERATURE
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              t(i,j,k,1,itemp)=20.
!              t(i,j,k,1,itemp)=20.-(zeta(i,j,1)-z_r(i,j,k))
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
            enddo
          enddo
        enddo
#    endif /* TEMPERATURE */

#   elif defined TANK
#    ifdef TEMPERATURE
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              t(i,j,k,1,itemp)=20.
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
            enddo
          enddo
        enddo
#    endif /* TEMPERATURE */
#   elif defined MOVING_BATHY
        cff=((0.92461)**2)/g
#    ifdef TEMPERATURE
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              t(i,j,k,1,itemp)=rho0*(1-cff*z_r(i,j,k)) - 1000.
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
            enddo
          enddo
        enddo
#    endif /* TEMPERATURE */
#   elif defined EQUATOR
        do k=1,N
          do j=JstrR,JendR
            do i=IstrR,IendR
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=10.0
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif /* TEMPERATURE */
#    ifdef SALINITY
              t(i,j,k,1,isalt)=35.0
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo

#   elif defined KH_INST
!
! Traditional KH instability:
!
!  Configuration from Pelletier et al.: linear shear within mixed layer
!
#    ifdef TEMPERATURE
        htemp0 = 10.
#    endif /* TEMPERATURE */
#    ifdef SALINITY
        hsalt0 = 10.
#    endif
        drho   = 1.2    ! Value to be replaced by script (originally 1.2)
!
        do j=JR_RANGE
          do i=IR_RANGE
            depth=h(i,j)
            zr0=-0.5*depth
            hup0=0.15*depth

            do k=1,N
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp) = 13.+ 1.*drho * (
     &                            0.1*(z_r(1,1,k)-zr0)/depth +
     &                          dtanh((z_r(1,1,k)-zr0)/htemp0)
     &                                              )  / Tcoef
!    &            +0.001*
!!   &                       (tanh(abs(z_r(1,1,k)-zr0)/hup0))*
!    &                       (1. -tanh((z_r(1,1,k)-zr0)/hup0)**2.)
!    &                                      *sin(2.*2.*pi/xl*xr(i,j))

              t(i,j,k,2,itemp) = t(i,j,k,1,itemp)
#    endif
#    ifdef SALINITY
              t(i,j,k,1,isalt) = 35. +0.05*drho* ( 1. -
     &                             (dtanh((z_r(1,1,k)-zr0)/hsalt0))**2 )
              t(i,j,k,2,isalt) = t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo

#   elif defined ISOLITON
#    define HORN_PROFILE
#    ifdef HORN_PROFILE
!
!  Horn, D.A., J. Imberger, & G.N. Ivey, (2001).
!  The degeneration of large-scale interfacial gravity waves in lakes.
!  J. Fluid Mech., 434:181-207.
!  --> initialy setup with variable salinity in the tank experiment of
!      Horn et al (2001); the conversion to temperature variations
!      through linear EOS gives unrealistic range
!
        rho1_exp = R0 - 10.  ! * 0.7286 /2.
        rho2_exp = R0 + 10.  ! * 0.7286 /2.
        h1_exp = 0.7 * h(1,1)
        amp_exp = 0.0391
        Lt=6.
#     ifdef TEMPERATURE
        do j=JstrR,JendR
          do i=IstrR,IendR
            do k=1,N
              d(k)=0.5*(rho1_exp+rho2_exp)
     &           + 0.5*(rho2_exp-rho1_exp)
     &               *tanh(exp(1.)*(-h1_exp
     &                  -2.*amp_exp*(xr(i,j)-Lt/2.)/Lt
     &                    +abs(z_r(i,j,k)))/(0.015/2.))   ! 0.005=dz_q
              t(i,j,k,1,itemp) = 13. - ( d(k) - R0 )  / Tcoef
              t(i,j,k,2,itemp) = t(i,j,k,1,itemp)
            enddo
          enddo
        enddo
#     endif /* TEMPERATURE */
#    else
        depth=h(1,1)
        gam=0.9
        del=0.3
        dz=2*(z_r(1,1,2)-z_r(1,1,1));
#     ifdef TEMPERATURE
        do j=JstrR,JendR
          do i=IstrR,IendR
            zint=depth*(del-1)+gam*del*depth*(2*xr(i,j)/xl - 1)
            do k=1,N
              t(i,j,k,1,itemp)=0.5*(5/Tcoef+10/Tcoef)
     &                        -0.5*(5/Tcoef-10/Tcoef)
     &                            *tanh((z_r(i,j,k)-zint)/dz);
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
            enddo
          enddo
        enddo
#     endif /* TEMPERATURE */
#    endif /* HORN_PROFILE */

#   elif defined GRAV_ADJ
#    ifdef TEMPERATURE
        do k=1,N
          do j=JstrR,JendR
            do i=IstrR,IendR
#     ifdef MPI
              if (i+iminmpi-1 .LE. LLm/2) then
#     else
              if (i .LE. LLm/2) then
#     endif
                t(i,j,k,1,itemp)=5/Tcoef
                t(i,j,k,2,itemp)=5/Tcoef
              else
                t(i,j,k,1,itemp)=10/Tcoef
                t(i,j,k,2,itemp)=10/Tcoef
              endif
            enddo
          enddo
        enddo
#    endif /* TEMPERATURE */

#   elif defined ACOUSTIC
#    ifdef TEMPERATURE

        do j=JstrR,JendR
          do i=IstrR,IendR
            do k=1,N
              t(i,j,k,1,itemp) = 13./Tcoef
              t(i,j,k,2,itemp) = t(i,j,k,1,itemp)
            enddo
          enddo
        enddo
#    endif /* TEMPERATURE */

#   elif defined INNERSHELF
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef INNERSHELF_EKMAN
#     ifdef TEMPERATURE
              t(i,j,k,1,itemp)=15.
#     endif /* TEMPERATURE */
#    else
#     define Z0 (-80.)
#     define THKNSS 50.
#     define Z1 0.
#     define STRAT 1000.
#     ifdef TEMPERATURE
              t(i,j,k,1,itemp)=14.+5.*THKNSS*log(
     &                         cosh((z_w(i,j,k  )-Z0)/THKNSS)
     &                        /cosh((z_w(i,j,k-1)-Z0)/THKNSS)
     &                            )/(z_w(i,j,k)-z_w(i,j,k-1))
     &               +((z_w(i,j,k)+z_w(i,j,k-1))/2.-Z1)/STRAT
#     endif /* TEMPERATURE */
#     undef Z0
#     undef THKNSS
#     undef Z1
#     undef STRAT
#    endif
#    ifdef TEMPERATURE
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif /* TEMPERATURE */
#    ifdef SALINITY
              t(i,j,k,1,isalt)=35.
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo

#   elif defined INTERNAL
        cff=((2.e-3)**2)/g
#    ifdef TEMPERATURE
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              t(i,j,k,1,itemp)=rho0*(1-cff*z_r(i,j,k)) - 1000.
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
            enddo
          enddo
        enddo
#    endif /* TEMPERATURE */

#   elif defined IGW
        cff=((2.e-3)**2)/g
        rhom=1027.037885
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              depth = -z_r(i,j,k)
#    ifdef EXPERIMENT1
              if(depth.le.50.) then
                sig = 25.0
              elseif( depth.gt.50. .and. depth.le.60.) then
                sig = 26.2 + rhom*(depth-50.)*cff
                sig = sig + 1.2*(depth-60.)/10.
              else
                sig = 26.2 + rhom*(depth-50.)*cff
              endif
#     ifdef TEMPERATURE
              t(i,j,k,1,itemp)=(R0-sig)/Tcoef
#     endif /* TEMPERATURE */
#    elif defined EXPERIMENT2
              if (depth.le.50.) then
               sig = 25.8
              else
               sig = 27.0
              endif
#     ifdef TEMPERATURE
              t(i,j,k,1,itemp)=(R0-sig)/Tcoef
#     endif /* TEMPERATURE */
#    else /* default EXPERIMENT3 */
              sig = 26.2 + rhom*(depth-50.)*cff
#     ifdef TEMPERATURE
              t(i,j,k,1,itemp)=(R0-sig)/Tcoef
#     endif /* TEMPERATURE */
#    endif
#    ifdef TEMPERATURE
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif /* TEMPERATURE */
            enddo
          enddo
        enddo

#   elif defined OVERFLOW
#    ifdef TEMPERATURE
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              t(i,j,k,1,itemp)=T0*(0.5-0.5*tanh(yr(i,j)/1000.-25.))
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
            enddo
          enddo
        enddo
#    endif /* TEMPERATURE */

#   elif defined RIVER 
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=10.  !4.+10.*exp(z_r(i,j,k)/50.)
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif /* TEMPERATURE */
#    ifdef SALINITY
              t(i,j,k,1,isalt)=36.
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo

#   elif defined SEAGRASS
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=10.  
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif /* TEMPERATURE */
#    ifdef SALINITY
              t(i,j,k,1,isalt)=36.
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo

#   elif defined SEAMOUNT
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=5.*15.*exp(z_r(i,j,k)/100.)
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif /* TEMPERATURE */
#    ifdef SALINITY
              t(i,j,k,1,isalt)=35.
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo

#   elif defined SHELFRONT
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=T0+2.5*tanh((yr(i,j)-50000.0)/20000.0)
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif /* TEMPERATURE */
#    ifdef SALINITY
              t(i,j,k,1,isalt)=S0  +  tanh((yr(i,j)-50000.0)/20000.0)
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo

#   elif defined UPWELLING
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    define Z0 (-35.)
#    define THKNSS 6.5
#    define Z1 (-75.)
#    define STRAT 150.
#    ifdef TEMPERATURE
             t(i,j,k,1,itemp)=14.+4.*THKNSS*log(
     &                        cosh((z_w(i,j,k )-Z0)/THKNSS)
     &                       /cosh((z_w(i,j,k-1)-Z0)/THKNSS)
     &                           )/(z_w(i,j,k)-z_w(i,j,k-1))
     &             +((z_w(i,j,k)+z_w(i,j,k-1))/2.-Z1)/STRAT
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif /* TEMPERATURE */
#    undef Z0
#    undef THKNSS
#    undef Z1
#    undef STRAT
#    ifdef SALINITY
              t(i,j,k,1,isalt)=35.
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo

#   elif defined REGIONAL
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    define Z0 (-80.)
#    define THKNSS 50.
#    define Z1 0.
#    define STRAT 1000.
#    ifdef TEMPERATURE
             t(i,j,k,1,itemp)=14.+5.*THKNSS*log(
     &                        cosh((z_w(i,j,k )-Z0)/THKNSS)
     &                       /cosh((z_w(i,j,k-1)-Z0)/THKNSS)
     &                           )/(z_w(i,j,k)-z_w(i,j,k-1))
     &              +((z_w(i,j,k)+z_w(i,j,k-1))/2.-Z1)/STRAT
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif /* TEMPERATURE */
#    undef Z0
#    undef THKNSS
#    undef Z1
#    undef STRAT
#    ifdef SALINITY
              t(i,j,k,1,isalt)=35.
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo

#   elif defined JET
!
! Compute analytical t and geostrophically adjusted
! zeta,ub,vb,u,v fields with no zonal perturbation first
!  --> store in climatological arrays
!  ... This is done both during initialization
!      and restart procedures
!
      endif  ! nrrec.eq.0.
!$acc end kernels
!
! This is a quick fix to get the initialisation reproducible
! (and right...) with OpenMP or tiling. This initialisation is performed
! on the whole domain (or MPI subdomain).
! Not optimal at all but the routine ana_jet should be totally
! re-written to make it properly
!
#    ifdef MPI
#     define LOCALLM Lmmpi
#     define LOCALMM Mmmpi
#    else
#     define LOCALLM Lm
#     define LOCALMM Mm
#    endif
C$OMP MASTER
      Istr1=1
      Iend1=LOCALLM
      Jstr1=1
      Jend1=LOCALMM
      call ana_jet_tile (Istr1,Iend1,Jstr1,Jend1,0)
C$OMP END MASTER
!
! Compute same fields with zonal perturbation
! --> store in initial arrays
!
      if (nrrec.eq.0) then  ! only during initialization
!
C$OMP MASTER
        call ana_jet_tile (Istr1,Iend1,Jstr1,Jend1,1)
C$OMP END MASTER
C$OMP BARRIER

#   elif defined TIDAL_FLAT
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=10.
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif
#    ifdef SALINITY
              t(i,j,k,1,isalt)=35.5
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo
#   elif defined ESTUARY
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=10.
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif
#    ifdef SALINITY
#     ifdef MPI
        cff  = (float(i +ii*Lm) )
#     else
        cff  = (float(i       ) )
#     endif
              t(i,j,k,1,isalt)= max(0., 35.5 * (1. - cff / LLm))  ! linear gradient from offshore to river
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo

#   elif defined DUNE
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=10.
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif
#    ifdef SALINITY
              t(i,j,k,1,isalt)=35.5
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo

#   elif defined SED_TOY
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=10.
              t(i,j,k,2,itemp)=t(i,j,k,1,itemp)
#    endif
#    ifdef SALINITY
              t(i,j,k,1,isalt)=35.5
              t(i,j,k,2,isalt)=t(i,j,k,1,isalt)
#    endif
            enddo
          enddo
        enddo

#   elif defined KILPATRICK
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=(288.95-273.16)+1.5*TANH(xr(i,j)/100.E+3)
              t(i,j,k,2,itemp)=(288.95-273.16)+1.5*TANH(xr(i,j)/100.E+3)
#    endif /* TEMPERATURE */
#    ifdef SALINITY
              t(i,j,k,1,isalt)=35.
              t(i,j,k,2,isalt)=35.
#    endif /* SALINITY */
            enddo
          enddo
        enddo

#   else
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
#    ifdef TEMPERATURE
              t(i,j,k,1,itemp)=20.
              t(i,j,k,2,itemp)=20.
#    endif /* TEMPERATURE */
#    ifdef SALINITY
              t(i,j,k,1,isalt)=35.
              t(i,j,k,2,isalt)=35.
#    endif
            enddo
          enddo
        enddo
#   endif
!
#  endif /*TRACERS */
      endif   !  nrrec.eq.0.
!
# endif /*SOLVE3D */
# if !defined JET
!$acc end kernels
# endif
#endif /* ANA_INITIAL */

#ifdef SOLVE3D
!
!--------------------------------------------------------------------
!  Initial conditions for passive tracer
!--------------------------------------------------------------------
!
# ifdef PASSIVE_TRACER
      do k=1,N
        do j=JR_RANGE
          do i=IR_RANGE
            if (.not.got_tini(itpas)) then
#  if defined REGIONAL && !defined ANA_PSOURCE
              if (h(i,j).lt.300.) then
!             if (z_r(i,j,k).lt.-800 .and.
!     &           z_r(i,j,k).gt.-1200 ) then
                t(i,j,k,1,itpas)=1
              else
                t(i,j,k,1,itpas)=0
              endif
#  elif defined SHOREFACE || defined RIP
            if (h(i,j).lt.2) then
              t(i,j,k,1,itpas)=1.
            else
              t(i,j,k,1,itpas)=0.
            endif
            t(i,j,k,2,itpas)=t(i,j,k,1,itpas)
#  elif defined KH_INST
              zdye10 = -0.65*depth
              t(i,j,k,1,itpas) =
     &            + 35 +0.05* drho* (
     &              1. - (dtanh( (z_r(1,1,k)-zdye10)/hsalt0 ) )**2.
     &                     )
!
              t(i,j,k,2,itpas) = t(i,j,k,1,itpas)

              zdye20 = -0.35*depth
              t(i,j,k,1,itpas+1) =
     &            + 35 +0.05* drho* (
     &              1. - (dtanh( (z_r(1,1,k)-zdye20)/hsalt0 ) )**2.
     &                     )!
              t(i,j,k,2,itpas+1) = t(i,j,k,1,itpas+1)
#  elif defined SED_TOY
            if (k.gt.int(N/2)) then
              t(i,j,k,1,itpas)=1.
            else
              t(i,j,k,1,itpas)=0.
            endif
#  elif defined TIDAL_FLAT || defined ESTUARY
            if (k.gt.int(N/2)) then
              t(i,j,k,1,itpas)=1.
            else
              t(i,j,k,1,itpas)=0.
            endif
#  elif defined FLASH_RIP
            if (h(i,j).lt.0.5) then
              t(i,j,k,1,itpas)=1.
            else
              t(i,j,k,1,itpas)=0.
            endif
#  else
            if (h(i,j).lt.300.) then
              t(i,j,k,1,itpas)=1.
            else
              t(i,j,k,1,itpas)=0.
            endif
#  endif
            endif
          enddo
        enddo
      enddo
# endif /* PASSIVE_TRACER */
!
!
!--------------------------------------------------------------------
!  Initial conditions for tracer-Substance
!--------------------------------------------------------------------
!
# ifdef SUBSTANCE
      do k=1,N
        do j=JR_RANGE
          do i=IR_RANGE
            do iv=itsubs1,itsubs2
              if (.not.got_tini(iv)) then
#  if defined REGIONAL && !defined ANA_PSOURCE
                t(i,j,k,1,iv)=cini_wat(iv)
#  elif defined TIDAL_FLAT || defined ESTUARY
                if (h(i,j).gt.0) then
                  t(i,j,k,1,iv)=cini_wat(iv)
                else
                  t(i,j,k,1,iv)=0.
                endif
#  elif defined DUNE && defined MUSTANG
                t(i,j,k,1,iv)=cini_wat(iv)
#  elif defined SED_TOY
                t(i,j,k,1,iv)=cini_wat(iv)
#  else
                if (h(i,j).lt.300) then
                  t(i,j,k,1,iv)=cini_wat(iv)
                else
                  t(i,j,k,1,iv)=0.
                endif
#  endif
              endif
            enddo
          enddo
        enddo
      enddo
# endif /*SUBSTANCE */

!----------------------------------------------------------------------
!  Initial conditions for sediment tracer variables.
!----------------------------------------------------------------------
!
# ifdef SEDIMENT
      do itrc=1,NST
        iv=itrc_sed+itrc-1
        do k=1,N
          do j=JR_RANGE
            do i=IR_RANGE
              if (.not.got_tini(iv)) then
                t(i,j,k,1,iv)=Csed(itrc)
                t(i,j,k,2,iv)=t(i,j,k,1,iv)
              endif
            enddo
          enddo
        enddo
      enddo

#  if defined SED_TOY && ( defined SED_TOY_RESUSP || \
                           defined SED_TOY_CONSOLID )
      do iv=itrc_mud,itrc_mud+NMUD-1
        do k=1,5
          do j=JR_RANGE
            do i=IR_RANGE
              if (.not.got_tini(iv)) then
!                t(i,j,k,1,iv)=1./7.5
                 t(i,j,k,1,iv)=0
                 t(i,j,k,2,iv)=t(i,j,k,1,iv)
              endif
            enddo
          enddo
        enddo
      enddo

#  endif
# endif
!
!----------------------------------------------------------------------
!  Initial conditions for skin temperature
!----------------------------------------------------------------------
!
# ifdef SST_SKIN
#  ifdef TEMPERATURE
      do j=JR_RANGE
        do i=IR_RANGE
          sst_skin(i,j)=t(i,j,N,1,itemp)
          dT_skin(i,j)=0.
        enddo
      enddo
#  endif /* TEMPERATURE */
# endif
!
!--------------------------------------------------------------------
!  Initial conditions for biology tracer variables.
!--------------------------------------------------------------------
!
# ifdef BIOLOGY
      do k=1,N
        do j=JR_RANGE
          do i=IR_RANGE
            if (.not.got_tini(iNO3_)) then
#  ifdef TEMPERATURE
              temp=t(i,j,k,1,itemp)
#  else
              temp=25.
#  endif /* TEMPERATURE */
               if (temp.lt.8.) then
                 SiO4=30.
               elseif (temp.ge.8. .and. temp.le.11.) then
                 SiO4=30.-((temp-8.)*(20./3.))
               elseif (temp.gt.11. .and. temp.le.13.) then
                 SiO4=10.-((temp-11.)*(8./2.))
               elseif (temp.gt.13. .and. temp.le.16.) then
                 SiO4=2.-((temp-13.)*(2./3.))
               elseif (temp.gt.16.) then
                 SiO4=0.
               endif
               t(i,j,k,1,iNO3_)=1.67+0.5873*SiO4+0.0144*SiO4**2
     &                               +0.0003099*SiO4**3
               t(i,j,k,2,iNO3_)=t(i,j,k,1,iNO3_)
            endif
#  ifdef PISCES
            if (.not.got_tini(iDIC_))  then
              t(i,j,k,1,iDIC_)=2150.
              t(i,j,k,2,iDIC_)=t(i,j,k,1,iDIC_)
            endif
            if (.not.got_tini(iTAL_))  then
              t(i,j,k,1,iTAL_)=2350.
              t(i,j,k,2,iTAL_)=t(i,j,k,1,iTAL_)
            endif
            if (.not.got_tini(iOXY_))  then
              t(i,j,k,1,iOXY_)=200.
              t(i,j,k,2,iOXY_)=t(i,j,k,1,iOXY_)
            endif
            if (.not.got_tini(iPOC_))  then
              t(i,j,k,1,iPOC_)=1.e-2
              t(i,j,k,2,iPOC_)=t(i,j,k,1,iPOC_)
            endif
            if (.not.got_tini(iPHY_))  then
              t(i,j,k,1,iPHY_)=0.01
              t(i,j,k,2,iPHY_)=t(i,j,k,1,iPHY_)
            endif
            if (.not.got_tini(iZOO_))  then
              t(i,j,k,1,iZOO_)=0.01
              t(i,j,k,2,iZOO_)=t(i,j,k,1,iZOO_)
            endif
            if (.not.got_tini(iDOC_))  then
              t(i,j,k,1,iDOC_)=5.
              t(i,j,k,2,iDOC_)=t(i,j,k,1,iDOC_)
            endif
            if (.not.got_tini(iFER_))  then
              t(i,j,k,1,iFER_)=6.e-4
              t(i,j,k,2,iFER_)=t(i,j,k,1,iFER_)
            endif
#   if ! defined key_pisces_npzd
            if (.not.got_tini(iCAL_))  then
              t(i,j,k,1,iCAL_)=0.01
              t(i,j,k,2,iCAL_)=t(i,j,k,1,iCAL_)
            endif
            if (.not.got_tini(iPO4_))  then
#    ifdef TEMPERATURE
              temp=t(i,j,k,1,itemp)
#    else
              temp=25.
#    endif /* TEMPERATURE */
              if (temp.lt.8.) then
                SiO4=30.
              elseif (temp.ge.8. .and. temp.le.11.) then
                SiO4=30.-((temp-8.)*(20./3.))
              elseif (temp.gt.11. .and. temp.le.13.) then
                SiO4=10.-((temp-11.)*(8./2.))
              elseif (temp.gt.13. .and. temp.le.16.) then
                SiO4=2.-((temp-13.)*(2./3.))
              elseif (temp.gt.16.) then
                SiO4=0.
              endif
              t(i,j,k,1,iPO4_)=(1.67+0.5873*SiO4+0.0144*SiO4**2
     &                               +0.0003099*SiO4**3)/16.
              t(i,j,k,2,iPO4_)=t(i,j,k,1,iPO4_)
            endif
            if (.not.got_tini(iSIL_))  then
              t(i,j,k,1,iSIL_)=91.51
              t(i,j,k,2,iSIL_)=t(i,j,k,1,iSIL_)
            endif
            if (.not.got_tini(iDIA_))  then
              t(i,j,k,1,iDIA_)=0.01
              t(i,j,k,2,iDIA_)=t(i,j,k,1,iDIA_)
            endif
            if (.not.got_tini(iMES_))  then
              t(i,j,k,1,iMES_)=0.01
              t(i,j,k,2,iMES_)=t(i,j,k,1,iMES_)
            endif
            if (.not.got_tini(iDSI_))  then
              t(i,j,k,1,iDSI_)=1.5e-3
              t(i,j,k,2,iDSI_)=t(i,j,k,1,iDSI_)
            endif
            if (.not.got_tini(iBFE_))  then
              t(i,j,k,1,iBFE_)=1.E-2*5E-6
              t(i,j,k,2,iBFE_)=t(i,j,k,1,iBFE_)
            endif
            if (.not.got_tini(iGOC_))  then
              t(i,j,k,1,iGOC_)=1.e-2
              t(i,j,k,2,iGOC_)=t(i,j,k,1,iGOC_)
            endif
            if (.not.got_tini(iSFE_))  then
              t(i,j,k,1,iSFE_)=1.e-2*5.E-6
              t(i,j,k,2,iSFE_)=t(i,j,k,1,iSFE_)
            endif
            if (.not.got_tini(iDFE_))  then
              t(i,j,k,1,iDFE_)=1.e-2*5.E-6
              t(i,j,k,2,iDFE_)=t(i,j,k,1,iDFE_)
            endif
            if (.not.got_tini(iGSI_))  then
              t(i,j,k,1,iGSI_)=1.e-2*0.15
              t(i,j,k,2,iGSI_)=t(i,j,k,1,iGSI_)
            endif
            if (.not.got_tini(iNFE_))  then
              t(i,j,k,1,iNFE_)=1.e-2*5.E-6
              t(i,j,k,2,iNFE_)=t(i,j,k,1,iNFE_)
            endif
            if (.not.got_tini(iNCH_))  then
              t(i,j,k,1,iNCH_)=1.e-2*12./55.
              t(i,j,k,2,iNCH_)=t(i,j,k,1,iNCH_)
            endif
            if (.not.got_tini(iDCH_))  then
              t(i,j,k,1,iDCH_)=1.e-2*12./55.
              t(i,j,k,2,iDCH_)=t(i,j,k,1,iDCH_)
            endif
            if (.not.got_tini(iNH4_))  then
              t(i,j,k,1,iNH4_)=1.e-2
              t(i,j,k,2,iNH4_)=t(i,j,k,1,iNH4_)
            endif
#   endif
#   ifdef key_ligand
            if (.not.got_tini(iLGW_))  then
              t(i,j,k,1,iLGW_)=1.e-3
              t(i,j,k,2,iLGW_)=t(i,j,k,1,iLGW_)
            endif
#   endif
#   ifdef key_pisces_quota
            if (.not.got_tini(iDON_))  then
              t(i,j,k,1,iDON_)=5.
              t(i,j,k,2,iDON_)=t(i,j,k,1,iDON_)
            endif
            if (.not.got_tini(iDOP_))  then
              t(i,j,k,1,iDOP_)=5.
              t(i,j,k,2,iDOP_)=t(i,j,k,1,iDOP_)
            endif
            if (.not.got_tini(iPON_))  then
              t(i,j,k,1,iPON_)=1.e-2
              t(i,j,k,2,iPON_)=t(i,j,k,1,iPON_)
            endif
            if (.not.got_tini(iPOP_))  then
              t(i,j,k,1,iPOP_)=1.e-2
              t(i,j,k,2,iPOP_)=t(i,j,k,1,iPOP_)
            endif
            if (.not.got_tini(iNPH_))  then
              t(i,j,k,1,iNPH_)=1.e-2
              t(i,j,k,2,iNPH_)=t(i,j,k,1,iNPH_)
            endif
            if (.not.got_tini(iPPH_))  then
              t(i,j,k,1,iPPH_)=1.e-2
              t(i,j,k,2,iPPH_)=t(i,j,k,1,iPPH_)
            endif
            if (.not.got_tini(iNDI_))  then
              t(i,j,k,1,iNDI_)=1.e-2
              t(i,j,k,2,iNDI_)=t(i,j,k,1,iNDI_)
            endif
            if (.not.got_tini(iPDI_))  then
              t(i,j,k,1,iPDI_)=1.e-2
              t(i,j,k,2,iPDI_)=t(i,j,k,1,iPDI_)
            endif
            if (.not.got_tini(iPIC_))  then
              t(i,j,k,1,iPIC_)=1.e-2
              t(i,j,k,2,iPIC_)=t(i,j,k,1,iPIC_)
            endif
            if (.not.got_tini(iNPI_))  then
              t(i,j,k,1,iNPI_)=1.e-2
              t(i,j,k,2,iNPI_)=t(i,j,k,1,iNPI_)
            endif
            if (.not.got_tini(iPPI_))  then
              t(i,j,k,1,iPPI_)=1.e-2
              t(i,j,k,2,iPPI_)=t(i,j,k,1,iPPI_)
            endif
            if (.not.got_tini(iPFE_))  then
              t(i,j,k,1,iPFE_)=1.e-2*5.E-6
              t(i,j,k,2,iPFE_)=t(i,j,k,1,iPFE_)
            endif
            if (.not.got_tini(iPCH_))  then
              t(i,j,k,1,iPCH_)=1.e-2*12./55.
              t(i,j,k,2,iPCH_)=t(i,j,k,1,iPCH_)
            endif
            if (.not.got_tini(iGON_))  then
              t(i,j,k,1,iGON_)=1.e-2
              t(i,j,k,2,iGON_)=t(i,j,k,1,iGON_)
            endif
            if (.not.got_tini(iGOP_))  then
              t(i,j,k,1,iGOP_)=1.e-2
              t(i,j,k,2,iGOP_)=t(i,j,k,1,iGOP_)
            endif
#   endif
#  elif defined BIO_NChlPZD
            if (.not.got_tini(iChla)) then
              t(i,j,k,1,iChla)=0.08
              t(i,j,k,2,iChla)=t(i,j,k,1,iChla)
            endif
            if (.not.got_tini(iPhy1)) then
              t(i,j,k,1,iPhy1)=0.1
              t(i,j,k,2,iPhy1)=t(i,j,k,1,iPhy1)
            endif
            if (.not.got_tini(iZoo1)) then
              t(i,j,k,1,iZoo1)=0.06
              t(i,j,k,2,iZoo1)=t(i,j,k,1,iZoo1)
            endif
            if (.not.got_tini(iDet1)) then
              t(i,j,k,1,iDet1)=0.02
              t(i,j,k,2,iDet1)=t(i,j,k,1,iDet1)
            endif
#   ifdef OXYGEN
            if (.not.got_tini(iO2))  then
              t(i,j,k,1,iO2)=250.0
              t(i,j,k,2,iO2)=t(i,j,k,1,iO2)
            endif
#   endif
#  elif defined BIO_N2ChlPZD2
            if (.not.got_tini(iNH4_))  then
              t(i,j,k,1,iNH4_)=0.1
              t(i,j,k,2,iNH4_)=t(i,j,k,1,iNH4_)
            endif
            if (.not.got_tini(iChla)) then
              t(i,j,k,1,iChla)=0.08
              t(i,j,k,2,iChla)=t(i,j,k,1,iChla)
            endif
            if (.not.got_tini(iPhy1)) then
              t(i,j,k,1,iPhy1)=0.06
              t(i,j,k,2,iPhy1)=t(i,j,k,1,iPhy1)
            endif
            if (.not.got_tini(iZoo1)) then
              t(i,j,k,1,iZoo1)=0.04
              t(i,j,k,2,iZoo1)=t(i,j,k,1,iZoo1)
            endif
            if (.not.got_tini(iDet1)) then
              t(i,j,k,1,iDet1)=0.02
              t(i,j,k,2,iDet1)=t(i,j,k,1,iDet1)
            endif
            if (.not.got_tini(iDet2)) then
              t(i,j,k,1,iDet2)=0.02
              t(i,j,k,2,iDet2)=t(i,j,k,1,iDet2)
            endif
#  elif defined BIO_BioEBUS
            if (.not.got_tini(iNO2_))  then
              t(i,j,k,1,iNO2_)=0.05*exp(z_r(i,j,k)/100.)
              t(i,j,k,2,iNO2_)=t(i,j,k,1,iNO2_)
            endif
            if (.not.got_tini(iNH4_))  then
              t(i,j,k,1,iNH4_)=0.1*exp(z_r(i,j,k)/100.)
              t(i,j,k,2,iNH4_)=t(i,j,k,1,iNH4_)
            endif
            if (.not.got_tini(iPhy1)) then
              t(i,j,k,1,iPhy1)=0.04*exp(z_r(i,j,k)/50.)
              t(i,j,k,2,iPhy1)=t(i,j,k,1,iPhy1)
            endif
            if (.not.got_tini(iPhy2)) then
              t(i,j,k,1,iPhy2)=0.06*exp(z_r(i,j,k)/50.)
              t(i,j,k,2,iPhy2)=t(i,j,k,1,iPhy2)
            endif
            if (.not.got_tini(iZoo1)) then
              t(i,j,k,1,iZoo1)=0.04*exp(z_r(i,j,k)/100.)
              t(i,j,k,2,iZoo1)=t(i,j,k,1,iZoo1)
            endif
            if (.not.got_tini(iZoo2)) then
              t(i,j,k,1,iZoo2)=0.04*exp(z_r(i,j,k)/100.)
              t(i,j,k,2,iZoo2)=t(i,j,k,1,iZoo2)
            endif
            if (.not.got_tini(iDet1)) then
              t(i,j,k,1,iDet1)=0.02
              t(i,j,k,2,iDet1)=t(i,j,k,1,iDet1)
            endif
            if (.not.got_tini(iDet2)) then
              t(i,j,k,1,iDet2)=0.02
              t(i,j,k,2,iDet2)=t(i,j,k,1,iDet2)
            endif
            if (.not.got_tini(iDON)) then
              t(i,j,k,1,iDON)=0.5*exp(z_r(i,j,k)/100.)
              t(i,j,k,2,iDON)=t(i,j,k,1,iDON)
            endif
            if (.not.got_tini(iO2)) then
              t(i,j,k,1,iO2)=250.*exp(z_r(i,j,k)/100.)
              t(i,j,k,2,iO2)=t(i,j,k,1,iO2)
            endif
#   ifdef NITROUS_OXIDE
            if (.not.got_tini(iN2O)) then
              t(i,j,k,1,iN2O)=-(0.008*exp(z_r(i,j,k)/100.)-0.01)
              if (t(i,j,k,1,iN2O).lt.0.) then
                t(i,j,k,1,iN2O)=0.
              endif
              t(i,j,k,2,iN2O)=t(i,j,k,1,iN2O)
            endif
#   endif
#  endif /* PISCES */
          enddo
        enddo
      enddo
# endif /* BIOLOGY */
#endif /* SOLVE3D */
!
!--------------------------------------------------------------------
!  Exchange boundary information
!--------------------------------------------------------------------
!

#if defined EW_PERIODIC || defined NS_PERIODIC || defined  MPI
# ifdef ANA_INITIAL
#  ifndef BASIN      
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                        zeta(START_2D_ARRAY,1))
      call exchange_u2d_tile (Istr,Iend,Jstr,Jend,
     &                        ubar(START_2D_ARRAY,1))
      call exchange_v2d_tile (Istr,Iend,Jstr,Jend,
     &                        vbar(START_2D_ARRAY,1))
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                        zeta(START_2D_ARRAY,2))
      call exchange_u2d_tile (Istr,Iend,Jstr,Jend,
     &                        ubar(START_2D_ARRAY,2))
      call exchange_v2d_tile (Istr,Iend,Jstr,Jend,
     &                        vbar(START_2D_ARRAY,2))
#  endif      
#  ifdef SOLVE3D
#   ifdef THREE_GHOST_POINTS_UV
      call exchange_u3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             u(START_2D_ARRAY,1,1))
      call exchange_v3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             v(START_2D_ARRAY,1,1))
      call exchange_u3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             u(START_2D_ARRAY,1,2))
      call exchange_v3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             v(START_2D_ARRAY,1,2))
#    if defined SOLID_BODY_ROT || defined SOLID_BODY_PER
      call exchange_u3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             u(START_2D_ARRAY,1,3))
      call exchange_v3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             v(START_2D_ARRAY,1,3))
#    endif
#    ifdef NBQ
      call exchange_w3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             wz(START_2D_ARRAY,0,1))
      call exchange_w3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             wz(START_2D_ARRAY,0,2))
#    endif
#   else
#    ifndef BASIN      
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
     &                        u(START_2D_ARRAY,1,1))
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
     &                        v(START_2D_ARRAY,1,1))
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
     &                        u(START_2D_ARRAY,1,2))
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
     &                        v(START_2D_ARRAY,1,2))
#    endif      
#    if defined SOLID_BODY_ROT || defined SOLID_BODY_PER
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,
     &                        u(START_2D_ARRAY,1,3))
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
     &                        v(START_2D_ARRAY,1,3))
#    endif
#    ifdef NBQ
      call exchange_w3d_tile (Istr,Iend,Jstr,Jend,
     &                        wz(START_2D_ARRAY,0,1))
      call exchange_w3d_tile (Istr,Iend,Jstr,Jend,
     &                        wz(START_2D_ARRAY,0,2))
#    endif
#   endif
#  endif /* SOLVE3D */
# endif /* ANA_INITIAL */
# ifdef SOLVE3D
#  ifdef TRACERS
      do itrc=1,NT
        if (.not.got_tini(itrc)) then
#   ifdef THREE_GHOST_POINTS_TS
          call exchange_r3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                                 t(START_2D_ARRAY,1,1,itrc))
          call exchange_r3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                                 t(START_2D_ARRAY,1,2,itrc))
#   else
          call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                            t(START_2D_ARRAY,1,1,itrc))
          call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                            t(START_2D_ARRAY,1,2,itrc))
#   endif
        endif
      enddo
#  endif /* TRACERS */
# endif
#endif
#undef IR_RANGE
#undef IU_RANGE
#undef JR_RANGE
#undef JV_RANGE

      return
      end

#ifdef JET
!======================================================================
!
      subroutine ana_jet_tile (Istr,Iend,Jstr,Jend, jet_perturb)
!
!======================================================================
!
!  Set initial conditions for momentum and tracer variables using
!  analytical expressions for the JET test case
!
!----------------------------------------------------------------------
!
# ifdef AGRIF
      use Agrif_Util
# endif
      implicit none
# include "param.h"
# include "grid.h"
# include "ocean2d.h"
# include "ocean3d.h"
# include "scalars.h"
# include "mpi_cpl.h"
# include "private_scratch.h"
# ifdef CLIMATOLOGY
#  include "climat.h"
# endif
      integer Istr,Iend,Jstr,Jend, i,j,k, itrc,
     &        jet_perturb, tndx
      real cff,x,y,Fyz,
     &     rhomaxs,bgdrhos,zs1,dzs,drhos,drhosfs,
     &     rhomaxn,bgdrhon,zn1,dzn,drhon,drhosfn,
     &     drhosf,z0,z0p,dzs_a,dzn_a,zk,
     &     Ljet,cff_perturb,ierr,tanh
      real*QUAD mssh,area,cff1
      real rhoprof(0:N,2)
      real ru(PRIVATE_2D_SCRATCH_ARRAY,N),
     &     rv(PRIVATE_2D_SCRATCH_ARRAY,N),
     &      P(PRIVATE_2D_SCRATCH_ARRAY,N),
     &     FC(PRIVATE_2D_SCRATCH_ARRAY),
     &    dZx(PRIVATE_2D_SCRATCH_ARRAY),
     &     rx(PRIVATE_2D_SCRATCH_ARRAY),
     &    dRx(PRIVATE_2D_SCRATCH_ARRAY),
     &    dR(PRIVATE_2D_SCRATCH_ARRAY,0:N),
     &    dZ(PRIVATE_2D_SCRATCH_ARRAY,0:N)
!$acc declare create(rhoprof,ru,rv,P)
# ifdef MPI
      include 'mpif.h'
      real*QUAD allmssh(1,NNODES),allarea(1,NNODES)
# endif
      real asymth,zz,zasym,dzasym
      asymth(zz,zasym,dzasym) =
     &         (zz-zasym)
     &         *sqrt(1+0.5*((zz-zasym)+abs(zz-zasym))**2/dzasym**2)
     &         +zasym
!
# include "compute_auxiliary_bounds.h"

! ---------------------------------------------------------------------
!                     BAROCLINIC JET TEST CASE
! ---------------------------------------------------------------------
!
! --- Jet and density profile parameters ---
!
      Ljet        = 1600.e3;  ! Jet width in km
!
! southern profile
      rhomaxs=27.75     ! max density
      bgdrhos=9.8e-6    ! background stratification
      zs1=-1000
      dzs=700
      drhos=1.4         ! gives strength of Phillips mode (0 to 1.4)
! northern profile
      rhomaxn=27.7573   ! max density
      bgdrhon=9.8e-6    ! background stratification
      zn1=-400
      dzn=300
! surface Charney mode
      drhosfs=1.5       ! northern anomaly  (0 for Phillips mode only)
      drhosfn=0.0       ! southern anomaly  (0 for Phillips mode only)
      z0=-300.          ! vertical penetration (-600 for strong Charney mode)
      drhosf=0.00       ! weakens Charney mode (canceled for 1.05)
      z0p=-110.         ! only used for drhosf non zero
!
! zonal perturbation
!
      if (jet_perturb .eq. 0) then
        cff_perturb = 0.0    ! zonal perturbation coefficient
        tndx=2
      else
        cff_perturb = 0.02
        tndx=1
      endif
!
! Set flat free surface zeta and associated vertical grid
! --> zeta will be diagnozed from pressure after
!     from first call to PGF routine
!
      knew=tndx  ! --> set_depth

!$acc kernels if(compute_on_device) default(present)
      zeta(:,:,tndx)=0.
!$acc end kernels

      call set_depth_tile (Istr,Iend,Jstr,Jend)

!$acc kernels if(compute_on_device) default(present)

!
! --- Build northern and southern density profiles ---
!
! First compute background density profiles associated with
! a gentle stratification that does not depend
! on jet side. It is there to ensure static stability.
      do k=1,N
        zk = z_r(istr,jstrR,k)
        rhoprof(k,1)=rhomaxs-bgdrhos*(zk+h(istr,jstrR))
        rhoprof(k,2)=rhomaxn-bgdrhon*(zk+h(istr,jstrR))
      enddo
!
! Second, get main north/south contrast with a distorded
! tanh shape of variable amplitude.Distorsion of the tanh
! is done with the asym function that increases
! stratification in the upper ocean
      dzs_a=1.3*dzs
      do k=1,N
        zk=asymth(z_r(istr,jstrR,k),zs1,dzs_a)
        rhoprof(k,1)=rhoprof(k,1)-drhos*(0.5+0.5*tanh((zk-zs1)/dzs))
      enddo
      dzn_a=1.3*dzn
      drhon=-(rhoprof(N,1)-rhoprof(N,2))/(0.5+0.5*tanh((zk-zn1)/dzn))
      do k=1,N
        zk=asymth(z_r(istr,jstrR,k),zn1,dzn_a)
        rhoprof(k,2)=rhoprof(k,2)-drhon*(0.5+0.5*tanh((zk-zn1)/dzn))
      enddo
      do k=1,N
        zk = z_r(istr,jstrR,k)
        rhoprof(k,1)=rhoprof(k,1)
     &             -drhosf*(exp((zk-z0p)/abs(z0p)))/(exp(1.))
        rhoprof(k,1)=rhoprof(k,1)
     &             -drhosfs*0.5*(1+tanh((zk-z0)/abs(z0)))/tanh(1.)
        rhoprof(k,2)=rhoprof(k,2)
     &             -drhosf*(exp((zk-z0p)/abs(z0p)))/(exp(1.))
        rhoprof(k,2)=rhoprof(k,2)
     &             -drhosfn*0.5*(1+tanh((zk-z0)/abs(z0)))/tanh(1.)
      enddo
!
! --- Fit southern and northern profiles ---
!
      do k=1,N
        do j=JstrR,JendR
          do i=IstrR,IendR
# ifdef MPI
            y=j-0.5+jj*Mm
            x=i-0.5+ii*Lm
# else
            y=j-0.5
            x=i-0.5
# endif
            y=y/real(MMm)-0.5
            x=x/real(LLm)
!
! Add vertical attenuation and XI perturbation
!
            y = y + cff_perturb*exp(z_r(i,j,k)/1000.)
     &                         *exp(-(x-0.5)**2/0.05)
     &                         *( 0.5*sin(2.*pi*x)
     &                           +0.5*sin(6.*pi*x) )
!
! Set Jet structure and place its center at y=pi/2
!
            y = y*pi*el/Ljet + pi/2

            if (y.lt.0) then
              Fyz=1.0
            elseif (y.gt.pi) then
              Fyz=0.0
            else
              Fyz=1.0-(y-sin(y)*cos(y))/pi
            endif
            rho(i,j,k)=Fyz*rhoprof(k,1)+(1.-Fyz)*rhoprof(k,2)
          enddo
        enddo
      enddo
!
! Fill tracer array and make exchanges
!
      do k=1,N
        do j=JstrR,JendR
          do i=IstrR,IendR
            t(i,j,k,tndx,1)=rho(i,j,k)
          enddo
        enddo
      enddo
!$acc end kernels

# if defined EW_PERIODIC || defined NS_PERIODIC || defined  MPI
#  ifdef THREE_GHOST_POINTS_TS
      call exchange_r3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                             t(START_2D_ARRAY,1,tndx,1))
#  else
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                        t(START_2D_ARRAY,1,tndx,1))
#  endif
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                        rho(START_2D_ARRAY,1))
# endif

!
! --- Compute zeta ---
!
! Get pressure using pressure gradient routine
!
      call prsgrd_tile (Istr,Iend,Jstr,Jend, ru,rv, P,
     &                  dR,dZ, FC,dZx,rx,dRx)

!$acc kernels if(compute_on_device) default(present)

!
! Compute zeta so as to cancel bottom pressure perturbation
!
      do j=Jstr,Jend
        do i=Istr,Iend
          zeta(i,j,tndx)=-P(i,j,1)/g
        enddo
      enddo
!
! Remove mean sea level (mssh)  !!! MPI DEBUG problems here !!!
!
      mssh=QuadZero
      area=QuadZero
      do j=Jstr,Jend
        do i=Istr,Iend
           cff1=1./(pm(i,j)*pn(i,j))
           mssh=mssh+cff1*zeta(i,j,tndx)
           area=area+cff1
        enddo
      enddo
      mssh=mssh/area
!$acc end kernels
# ifdef MPI
#  ifdef OPENACC
      NOTREADY
#  endif
      call MPI_ALLGATHER(mssh,1,MPI_DOUBLE_PRECISION,
     &                allmssh,1,MPI_DOUBLE_PRECISION,
     &                          MPI_COMM_WORLD,ierr)
      call MPI_ALLGATHER(area,1,MPI_DOUBLE_PRECISION,
     &                allarea,1,MPI_DOUBLE_PRECISION,
     &                          MPI_COMM_WORLD,ierr)
      mssh=QuadZero
      area=QuadZero
      do i=1,NNODES
        mssh=mssh+allmssh(1,i)*allarea(1,i)
        area=area+allarea(1,i)
      enddo
      mssh=mssh/area
# endif
      mssh=nint(mssh*1.d6,kind=8)/1.d6 ! round to 6 digits

!$acc kernels if(compute_on_device) default(present)

      do j=Jstr,Jend
        do i=Istr,Iend
          zeta(i,j,tndx)=zeta(i,j,tndx)
     &                          - mssh
        enddo
      enddo
!
! Apply zeta boundary conditions
!
      if (SOUTHERN_EDGE) then
        do i=Istr,Iend                     ! Southern edge gradient BC
          zeta(i,Jstr-1,tndx)=zeta(i,Jstr,tndx)
          enddo
      endif
      if (NORTHERN_EDGE) then
        do i=istr,iend                     ! Northern edge gradient BC
          zeta(i,Jend+1,tndx)=zeta(i,Jend,tndx)
        enddo
      endif
!$acc end kernels

# if defined EW_PERIODIC || defined NS_PERIODIC || defined  MPI
      call exchange_r2d_tile (istr,iend,jstr,jend,
     &                        zeta(START_2D_ARRAY,tndx))
# endif

!
! --- Get gressure gradient (ru,rv) using updated zeta field ---
!
      call set_depth_tile (Istr,Iend,Jstr,Jend)
      call prsgrd_tile (Istr,Iend,Jstr,Jend, ru,rv, P,
     &                  dR,dZ, FC,dZx,rx,dRx)

!$acc kernels if(compute_on_device) default(present)
# define WORK rho1
      do k=1,N
        do j=JstrV,Jend
          do i=Istr,Iend
            WORK(i,j,k)=rv(i,j,k)
          enddo
        enddo
      enddo
      if (NORTHERN_EDGE) then
        do k=1,N
          do i=Istr,Iend
            WORK(i,Jend+1,k)=0.
          enddo
        enddo
      endif
      if (SOUTHERN_EDGE) then
        do k=1,N
          do i=Istr,Iend
            WORK(i,JstrV-1,k)=0.
          enddo
        enddo
      endif
!$acc end kernels

# if defined EW_PERIODIC || defined NS_PERIODIC || defined  MPI
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,
     &                        WORK(START_2D_ARRAY,1))
# endif

!$acc kernels if(compute_on_device) default(present)
      do k=1,N
        do j=JstrV-1,Jend+1
          do i=Istr-1,Iend+1
            rv(i,j,k)=WORK(i,j,k)
          enddo
        enddo
      enddo

      do k=1,N
        do j=Jstr,Jend
          do i=IstrU,Iend
            WORK(i,j,k)=ru(i,j,k)
          enddo
        enddo
      enddo
      if (NORTHERN_EDGE) then
        do k=1,N
          do i=IstrU,Iend
            WORK(i,Jend+1,k)=0.
          enddo
        enddo
      endif
      if (SOUTHERN_EDGE) then
        do k=1,N
          do i=Istr,Iend
            WORK(i,0,k)=0.
          enddo
        enddo
      endif
!$acc end kernels

# if defined EW_PERIODIC || defined NS_PERIODIC || defined  MPI
      call exchange_u3d_tile (istr,iend,jstr,jend,
     &                        WORK(START_2D_ARRAY,1))
# endif

!$acc kernels if(compute_on_device) default(present)
      do k=1,N
        do j=Jstr-1,Jend+1
          do i=Istr-1,Iend+1
            ru(i,j,k)=WORK(i,j,k)
          enddo
        enddo
      enddo
# undef WORK
!
! --- Compute geostrophic u,v velocities ---
!
      do k=1,N
        do j=Jstr,Jend
          do i=Istr,Iend
            cff=0.5*(pm(i-1,j)+pm(i,j))
     &             *(pn(i-1,j)+pn(i,j))
     &             /(Hz(i-1,j,k)+Hz(i,j,k))/f(i,j)
            u(i,j,k,tndx) = 0.25*cff*
     &            (rv(i,j,k)+rv(i,j+1,k)+rv(i-1,j,k)+rv(i-1,j+1,k))
          enddo
        enddo
      enddo
      do k=1,N
        do j=Jstr,Jend
          do i=Istr,Iend
            cff=0.5*(pm(i,j)+pm(i,j-1))
     &             *(pn(i,j)+pn(i,j-1))
     &             /(Hz(i,j,k)+Hz(i,j-1,k))/f(i,j)
            v(i,j,k,tndx) =-0.25*cff*
     &            (ru(i,j,k)+ru(i+1,j,k)+ru(i,j-1,k)+ru(i+1,j-1,k))
          enddo
        enddo
      enddo
!
! --- Compute barotropic u velocity ---
!
      ubar(:,:,tndx)=0.
      do k=1,N
        do j=Jstr,Jend
          do i=Istr,Iend
            ubar(i,j,tndx)=ubar(i,j,tndx)+u(i,j,k,tndx)
     &                             *0.5*(Hz(i-1,j,k)+Hz(i,j,k))
          enddo
        enddo
      enddo
      do j=Jstr,Jend
        do i=Istr,Iend
          ubar(i,j,tndx) = ubar(i,j,tndx) /
     &                  ( 0.5*(zeta(i,j,tndx)+zeta(i-1,j,tndx))
     &                  + 0.5*(h(i,j)+h(i-1,j)) )
        enddo
      enddo
!
! --- Compute barotropic v velocity ---
!
      vbar(:,:,tndx)=0.
      do k=1,N
        do j=Jstr,Jend
          do i=Istr,Iend
            vbar(i,j,tndx)=vbar(i,j,tndx)+v(i,j,k,tndx)
     &                           *0.5*(Hz(i,j-1,k)+Hz(i,j,k))
          enddo
        enddo
      enddo
      do j=Jstr,Jend
        do i=Istr,Iend
          vbar(i,j,tndx) = vbar(i,j,tndx) /
     &                ( 0.5*(zeta(i,j,tndx)+zeta(i,j-1,tndx))
     &                + 0.5*(h(i,j)+h(i,j-1)) )
        enddo
      enddo

# ifdef CLIMATOLOGY
      if (jet_perturb .eq. 0) then
!
! Fill climatological arrays. MPI exchanges are done in analytical
!
        do k=1,N
          do j=Jstr,Jend
            do i=Istr,Iend
              uclm(i,j,k) = u(i,j,k,tndx)
              vclm(i,j,k) = v(i,j,k,tndx)
              tclm(i,j,k,1) = t(i,j,k,tndx,1)
            enddo
          enddo
        enddo
        do j=Jstr,Jend
          do i=Istr,Iend
            ssh(i,j) = zeta(i,j,tndx)
          enddo
        enddo
        do j=Jstr,Jend
          do i=Istr,Iend
            ubclm(i,j) = ubar(i,j,tndx)
          enddo
        enddo
        do j=Jstr,Jend
          do i=Istr,Iend
            vbclm(i,j) = vbar(i,j,tndx)
          enddo
        enddo
      endif
# endif
!$acc end kernels
      return
      end
#endif /* JET */

