!======================================================================
! 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"
!
!  ANALYTICAL PACKAGE:
!----------------------------------------------------------------------
!
!  This package is used to provide various analytical fields to the
!  model when appropriate.
!
!  Routines:
!
!  ana_bmflux     Analytical kinematic bottom momentum flux.
!  ana_btflux     Analytical kinematic bottom flux of tracer
!                            type variables.
!  ana_morphodyn  Analytical dynamic bathymetry
!  ana_bsedim     Analytical bottom sediment grain size
!                            and density.
!  ana_smflux     Analytical kinematic surface momentum flux
!                            (wind stress).
!  ana_srflux     Analytical kinematic surface shortwave
!                            radiation.
!  ana_ssh        Analytical sea surface height climatology.
!  ana_sst        Analytical sea surface temperature and dQdSST
!                            which are used during heat flux correction.
!  ana_sss        Analytical sea surface salinity which is used
!                            during salt flux correction.
!  ana_stflux     Analytical kinematic surface flux of tracer type
!                            variables.
!  ana_tclima     Analytical tracer climatology fields.
!  ana_uclima     Analytical tracer climatology fields.
!  ana_nbq_clima  Analytical NBQ variable climatology fields
!  ana_wwave      Analytical wind induced wave amplitude,
!                            direction and period.
!  ana_sediment   Analytical sediment
!  ana_psource    Analytical point source
!  ana_bry        Analytical boundary forcing.
!  ana_nbq_bry    Analytical boundary forcing for NBQ variables
!  ana_bry_bio    Analytical boundary forcing for BIO variables
!
!----------------------------------------------------------------------
!
#if !defined OPENMP
      integer function omp_get_thread_num()
      omp_get_thread_num=0
      return
      end
      integer function omp_get_num_threads()
      omp_get_num_threads=1
      return
      end
#endif
!
!======================================================================
!                   subroutine ana_bmflux
!======================================================================
!
#ifdef ANA_BMFLUX
      subroutine ana_bmflux_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  This routine sets kinematic bottom momentum flux (bottom stress)
! "bustr" and "bvstr" [m^2/s^2] using an analytical expression.
!----------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "grid.h"
# include "forces.h"
# include "scalars.h"
      integer Istr,Iend,Jstr,Jend, i,j
!
# include "compute_auxiliary_bounds.h"
!
      do j=JstrR,JendR
        do i=Istr,IendR
          bustr(i,j)=0.
        enddo
      enddo
      do j=Jstr,JendR
        do i=IstrR,IendR
          bvstr(i,j)=0.
        enddo
      enddo
      return
      end
#endif /* ANA_BMFLUX */
!
!======================================================================
!                   subroutine ana_btflux
!======================================================================
!
#if defined SOLVE3D && defined TRACERS
# if defined ANA_BTFLUX || defined ANA_BSFLUX || defined ANA_BPFLUX
      subroutine ana_btflux_tile (Istr,Iend,Jstr,Jend, itrc)
!
!----------------------------------------------------------------------
!  This routine sets kinematic bottom flux of tracer type variables
!  [tracer units m/s].
!
!  On Input:
!     itrc      Tracer type array index.
!----------------------------------------------------------------------
!
      implicit none
#  include "param.h"
#  include "grid.h"
#  include "forces.h"
#  include "scalars.h"
      integer itrc, Istr,Iend,Jstr,Jend, i,j
!
#  include "compute_auxiliary_bounds.h"
!
!$acc kernels if(compute_on_device) default(present)
!
!----------------------------------------------------------------------
! Set kinematic bottom heat flux [degC m/s] at horizontal RHO-points.
!----------------------------------------------------------------------
!
#  ifdef TEMPERATURE
      if (itrc.eq.itemp) then
        do j=JstrR,JendR
          do i=IstrR,IendR
#   if defined SINGLE_COLUMN
#    ifdef KATO_PHILIPS
            btflx(i,j,itemp)=-1.0e-5*(0.01)**2/(g*Tcoef/rho0)
#    elif defined WILLIS_DEARDORFF || defined DIURNAL_CYCLE
            btflx(i,j,itemp)=-1.0e-5*0.01
#    endif
#   else
            btflx(i,j,itemp)=0.
#   endif
          enddo
        enddo
      endif
#  endif /* TEMPERATURE */
!
!----------------------------------------------------------------------
!  Set kinematic bottom salt flux (m/s) at horizontal RHO-points,
!  scaling by bottom salinity is done in STEP3D.
!----------------------------------------------------------------------
!
#  ifdef SALINITY
      if (itrc.eq.isalt) then
        do j=JstrR,JendR
          do i=IstrR,IendR
            btflx(i,j,isalt)=0.
          enddo
        enddo
      endif
#  endif /* SALINITY */
!
!----------------------------------------------------------------------
!  Set kinematic surface flux of additional tracers,
! for example sediments, bio..., to zero
!----------------------------------------------------------------------
!
#  if defined SALINITY && defined TEMPERATURE
      if (itrc.ne.isalt .and. itrc.ne.itemp) then
#  endif /* SALINITY && TEMPERATURE*/
#  if !defined SALINITY && defined TEMPERATURE
      if (itrc.ne.itemp) then
#  endif /* !SALINITY && TEMPERATURE*/
#  if defined SALINITY && !defined TEMPERATURE
      if (itrc.ne.isalt) then
#  endif /* SALINITY && !TEMPERATURE*/
        do j=JstrR,JendR
          do i=IstrR,IendR
            btflx(i,j,itrc)=0.
          enddo
        enddo
#  if defined SALINITY || defined TEMPERATURE
      endif
#  endif /* SALINITY || TEMPERATURE*/     
!$acc end kernels 

      return
      end
# endif /* ANA_BTFLUX */
#endif /* SOLVE3D && TRACERS */
!
!======================================================================
!                   subroutine ana_moving_bathy
!======================================================================
!
#ifdef ANA_MORPHODYN
      subroutine ana_morphodyn (tile)
      implicit none
# include "param.h"
      integer tile
# ifdef  ALLOW_SINGLE_BLOCK_MODE
C$    integer  trd, omp_get_thread_num
# endif
# include "compute_tile_bounds.h"
      call ana_morphodyn_tile (Istr,Iend,Jstr,Jend)
      return
      end

      subroutine ana_morphodyn_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  This routine sets bathymetry increment
!  dh using an analytical expression:
!  example for SEAMOUNT test case
!----------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "grid.h"
# include "forces.h"
# include "scalars.h"
      integer Istr,Iend,Jstr,Jend, i,j
      real cff,cff1
!
# include "compute_auxiliary_bounds.h"
!
# ifdef SEAMOUNT
      cff=(1./50.e3)**2
      cff1=0.5*xl+time*1.
      ! 0.1m/s translation
      ! starting for mid-basin
      ! for example purpoe
      ! has to be coherent with initial bathymetry
      ! definition in ana_grid.F
      do j=JstrR,JendR
        do i=IstrR,IendR
          dh(i,j)=4500*(1-0.6*exp(-cff*((xr(i,j)-cff1  )**2+
     &                                  (yr(i,j)-0.5*el)**2)))
     &                                                 -h(i,j)
        enddo
      enddo
# elif defined MOVING_BATHY
      cff=1./(0.03686)**2
      cff1=1.e-3*sin(2.*pi/10.05*time)
      do j=JstrR,JendR
        do i=IstrR,IendR
          dh(i,j)=0.394 - 0.1*exp(-cff*(xr(i,j)-cff1)**2)
     &                                           - h(i,j)
        enddo
      enddo
# else
      dh(:,:) =???
# endif
      return
      end
#endif /* ANA_MORPHODYN */
!
!======================================================================
!                   subroutine ana_bsedim
!======================================================================
!
#if defined ANA_BSEDIM && defined BBL
      subroutine ana_bsedim (tile)
      implicit none
# include "param.h"
      integer tile
# ifdef ALLOW_SINGLE_BLOCK_MODE
C$    integer  trd, omp_get_thread_num
# endif
# include "compute_tile_bounds.h"
      call ana_bsedim_tile   (Istr,Iend,Jstr,Jend)
      return
      end
      subroutine ana_bsedim_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  This routine sets initial bottom sediment grain diameter size [m]
!  and density used in the bottom boundary formulation [kg/m^3].
!----------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "bbl.h"
# include "grid.h"
# include "scalars.h"
      integer Istr,Iend,Jstr,Jend, i,j
!
# include "compute_extended_bounds.h"
!
# ifdef REGIONAL
!
! taucb=critical threshold stress for initiation of motion
! (=bedload for coarse grains).
! critical suspension stress: ustar_crit=0.8*w_set
!
! determine taucb from Shields curve, fit provided by
! Soulsby & Whitehouse 1997, Threshold of sediment motion
! in coastal environments, Proc. Pacific Coasts and Ports
! '97 Conf., pp 149--154, Univ Canterbury, Nw Zealand.
!
! visk=1.3e-3/rhow; (kinem. visc., nu=mu/rhow)
! D=d50*(g*(Sdens/rhow-1)/(visk^2))^0.33333
! thetcr=0.3./(1+1.2*D) + 0.055*(1-exp(-0.02*D))
! taucb=thetcr.*(g*(sdens-rhow).*d50);
!
! Souslby's (1997) estimate of settling velocity
!   w_set = visk*(sqrt(10.36^2+1.049*D^3)-10.36)/d50 [m/s]
! with D as above
!
      do j=JstrR,JendR
        do i=IstrR,IendR
          Ssize(i,j)=1.5e-4     ! d50 [m]
          Sdens(i,j)=2650.0     ! rho sediment [kg/m^3]
          taucb(i,j)=0.16/rho0  ! critical bedload stress [m^2/s^2]
          w_set(i,j)=0.013      ! analytical settling velocity [m/s]
          Hripple(i,j)=0.01     ! analytical initial ripple height [m]
          Lripple(i,j)=0.10     ! analytical initial ripple length [m]
        enddo
      enddo
# else
      !ANA_BSEDIM: no values provided for SSIZE and SDENS.
      do j=JstrR,JendR
        do i=IstrR,IendR
          Ssize(i,j)=1.5e-4     ! d50 [m]
          Sdens(i,j)=2650.0     ! rho sediment [kg/m^3]
          taucb(i,j)=0.16/rho0  ! critical bedload stress [m^2/s^2]
          w_set(i,j)=0.013      ! analytical settling velocity [m/s]
          Hripple(i,j)=0.01     ! analytical initial ripple height [m]
          Lripple(i,j)=0.10     ! analytical initial ripple length [m]
        enddo
      enddo
# endif
      return
      end
#endif /* ANA_BSEDIM && BBL */
!
!======================================================================
!                   subroutine ana_smflux
!======================================================================
!
#ifdef ANA_SMFLUX
      subroutine ana_smflux_tile (Istr,Iend,Jstr,Jend)
!
!  Sets kinematic surface momentum flux (wind stress) "sustr" and "svstr"
!  [m^2/s^2] using an analytical expression.
!
# ifdef AGRIF
      use Agrif_UTIL
# endif
      implicit none
# include "param.h"
# include "ocean2d.h"
# include "grid.h"
# include "forces.h"
# include "scalars.h"
      integer Istr,Iend,Jstr,Jend, i,j
      real Ewind, Nwind, dircoef, windamp, wcos,wsin
      real cff1, cff2, cff3
!     data windamp /0./
      data Ewind, Nwind, dircoef /0., 0., 0./
!     save windamp
!
# include "compute_extended_bounds.h"
!
!  Set kinematic surface momentum flux (wind stress) component in the
!  XI-direction (m^2/s^2) at horizontal U-points.
!

      windamp = 0.

!$acc kernels if(compute_on_device) default(present)
# ifdef BASIN
      cff1=0.0001 * 0.5*(1.+tanh((time-6.*86400.)/(3.*86400.)))
      cff2=2.*pi/el
      do j=JstrR,JendR
        do i=IstrR,IendR
          sustr(i,j)=-cff1*cos(cff2*yr(i,j))
        enddo
      enddo
# elif defined SINGLE_COLUMN && defined KATO_PHILIPS
      do j=JstrR,JendR
        do i=IstrR,IendR
          sustr(i,j)=(0.01)**2
        enddo
      enddo
# elif defined SINGLE_COLUMN && defined FORCED_DBLEEK
      do j=JstrR,JendR
        do i=IstrR,IendR
          sustr(i,j)=-(0.0025)**2
        enddo
      enddo
# elif defined CANYON
      do j=JstrR,JendR
        do i=IstrR,IendR
          sustr(i,j)=0.0001*0.5*sin(2.*pi*tdays/10.)*
     &               (1.-tanh((yr(i,j)-0.5*el)/10000.))
        enddo
      enddo
# elif defined EQUATOR
      do j=JstrR,JendR
        do i=IstrR,IendR
          sustr(i,j)=-0.05/rho0
        enddo
      enddo
# elif defined UPWELLING
      if (tdays.le.2.) then
        windamp=-0.1*sin(pi*tdays/4.)/rho0
      else
        windamp=-0.1/rho0
      endif
      do j=JstrR,JendR
        do i=IstrR,IendR
          sustr(i,j)=windamp
        enddo
      enddo
# elif defined SHOREFACE
#  ifndef MRL_WCI
      windamp=tanh(24.0*dt*sec2day*float(iic-ntstart))
      wcos = cos(10.0/180.0*pi)
      do j=jstr-1,jend+1
        do i=istr-1,iend+1
          cff2  = h(i,j)+zeta(i,j,kstp)
          cff1  = 0.05/sqrt(9.8*cff2)  ! 0.6 is epsilon_b/rho
          sustr(i,j)=cff1*wcos*windamp
        enddo
      enddo
#  else
      do j=JstrR,JendR
        do i=IstrR,IendR
          sustr(i,j)=0.
        enddo
      enddo
#  endif

# elif defined SED_TOY
#  ifdef SED_TOY_ROUSE
      cff1=0.001875
#  elif defined SED_TOY_FLOC_0D
      cff1=0.0
#  elif defined SED_TOY_FLOC_1D
      cff1=0.001875 * sin(2. * pi / (12.0*3600) * time)
#  else
#   if defined SED_TOY_RESUSP || defined SED_TOY_FLOC
      if (tdays.lt.1.8) then
        cff1=0.001
      else if ((tdays.ge.3).and.(tdays.le.4.)) then
        cff1=0.0005
#   elif defined SED_TOY_CONSOLID
      if ((tdays.ge.0.7).and.((tdays.le.37.))) then
        cff1=0.001
#   endif
      else
        cff1=0.
      endif
#  endif

      do j=JstrR,JendR
        do i=IstrR,IendR
          sustr(i,j)=cff1
        enddo
      enddo

# else
      do j=JstrR,JendR
        do i=IstrR,IendR
          sustr(i,j)=0.
        enddo
      enddo
# endif
!
!  Set kinematic surface momentum flux (wind stress) component in the
!  ETA-direction (m^2/s^2) at horizontal V-points.
!
# ifdef INNERSHELF
      do j=JstrR,JendR
        do i=IstrR,IendR
          svstr(i,j)=0.07/rho0
        enddo
      enddo
# elif defined SHOREFACE
#  ifndef MRL_WCI
      windamp=tanh(24.0*dt*sec2day*float(iic-ntstart))
      wsin = sin(10.0/180.0*pi)
      do j=jstr-1,jend+1
        do i=istr-1,iend+1
          cff2  = h(i,j)+zeta(i,j,kstp)
          cff1  = 0.05/sqrt(9.8*cff2)  ! 0.6 is \epsilon_b/\rho
          svstr(i,j)=cff1*wsin*windamp
        enddo
      enddo
#  else
      do j=JstrR,JendR
        do i=IstrR,IendR
          svstr(i,j)=0.
        enddo
      enddo
#  endif
# else
      do j=JstrR,JendR
        do i=IstrR,IendR
          svstr(i,j)=0.
        enddo
      enddo
# endif
!$acc end kernels
      return
      end
#endif /* ANA_SMFLUX */
!
!======================================================================
!                   subroutine ana_srflux
!======================================================================
!
#ifdef ANA_SRFLUX
      subroutine ana_srflux_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  This subroutine sets kinematic surface solar shortwave radiation
!  flux "srflx" (degC m/s) using an analytical expression.
!----------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "grid.h"
# include "forces.h"
# include "scalars.h"
      integer Istr,Iend,Jstr,Jend, i,j
!
# include "compute_auxiliary_bounds.h"
!
!  Set kinematic surface solar shortwave radiation [degC m/s] at
!  horizontal RHO-points.
!
!$acc kernels if(compute_on_device) default(present)
# if defined INNERSHELF && !defined INNERSHELF_EKMAN
      do j=JstrR,JendR
        do i=IstrR,IendR
          srflx(i,j)=264.91/(rho0*Cp)
#  ifdef DIURNAL_INPUT_SRFLX
          srflxbio(i,j)=srflx(i,j)
#  endif
#  ifdef DIURNAL_CYCLE
          srflx(i,j)=srflx(i,j)*max(cos(2.*pi*(time/86400.-0.5)),0.)
#  endif
        enddo
      enddo
# else
      do j=JstrR,JendR
        do i=IstrR,IendR
          srflx(i,j)=0.   ! default is 0.
#  ifdef DIURNAL_INPUT_SRFLX
          srflxbio(i,j)=srflx(i,j)
#  endif
#  ifdef DIURNAL_CYCLE
          srflx(i,j)=264.91/(rho0*Cp)
     &                     *max(cos(2.*pi*(time/86400.-0.5)),0.)
#  endif
        enddo
      enddo
# endif /* INNERSHELF */
!$acc end kernels
      return
      end
#endif /* ANA_SRFLUX */
!
!======================================================================
!                   subroutine ana_ssh
!======================================================================
!
#if defined ANA_SSH && defined ZCLIMATOLOGY
      subroutine ana_ssh (tile)
      implicit none
# include "param.h"
      integer tile
# ifdef ALLOW_SINGLE_BLOCK_MODE
C$    integer  trd, omp_get_thread_num
# endif
# include "compute_tile_bounds.h"
      call ana_ssh_tile (Istr,Iend,Jstr,Jend)
      return
      end
!
      subroutine ana_ssh_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  This routine sets analytical sea surface height climatology [m].
!----------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "grid.h"
# include "ocean2d.h"
# include "climat.h"
# include "scalars.h"
# include "coupling.h"
      integer Istr,Iend,Jstr,Jend, i,j
# if defined INTERNAL || defined ANA_TIDES
      real U0,omega,kwave,ETA0
# endif
# if defined RIP && defined ANA_TIDES
      real surelev,omega1,omega2
      real tide_ampli_moy,tide_ampli_var
      real*QUAD mssh,area,cff1
#  ifdef MPI
#   include "mpi_cpl.h"
      include 'mpif.h'
      integer ierr
      real*QUAD allmssh(1,NNODES),allarea(1,NNODES)
#  endif
# endif
!
# include "compute_auxiliary_bounds.h"
!
!  Set sea surface height (meters).
!
# ifdef JET
!
!    --> climatology computed in ana_initial
!

# elif defined RIP && defined ANA_TIDES
!
! Remove mean sea level (mssh)
!
      mssh=QuadZero
      area=QuadZero
      do j=Jstr,Jend
        do i=Istr,Iend
          cff1=1./(pm(i,j)*pn(i,j))
          mssh=mssh+cff1*Zt_avg1(i,j)
          area=area+cff1
        enddo
      enddo
      mssh=mssh/area
#  ifdef MPI
      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
      omega1=2.*pi/(14.*86400.)
      omega2=2.*pi/(12.*3600. )
      tide_ampli_moy=1.
      tide_ampli_var=30./70.
      surelev=1.
      do j=JstrR,JendR
        do i=IstrR,IendR
          ssh(i,j)=zeta(i,j,knew)-mssh
     &                +tide_ampli_moy*(1.+tide_ampli_var
     &                *sin(time*omega1))
     &                *sin(time*omega2)
     &                    +time*surelev/31536000.
        enddo
      enddo

# elif (defined TIDAL_FLAT  || defined ESTUARY ) && defined ANA_TIDES
      omega=2.*pi/(12.0*3600)
      ETA0=2.0
      do j=JstrR,JendR
        do i=IstrR,IendR
          ssh(i,j)=ETA0*sin(omega*time)
        enddo
      enddo
#  if defined ESTUARY_TIDE
      ! adding another wave to modulate tide
      omega=2.*pi/((24.0+50./60.)/2.*3600)
      ETA0=2.0
      do j=JstrR,JendR
        do i=IstrR,IendR
          ssh(i,j) = ssh(i,j) + ETA0*sin(omega*time)
        enddo
      enddo
#  endif

# elif defined INTERNAL || defined ANA_TIDES
      U0=0.02
      omega=2.*pi/(12.4*3600)
      kwave=sqrt(((omega*omega)-(f(1,1)*f(1,1)))/(g*h(1,1)))
      ETA0=kwave*h(1,1)*U0/omega
      do j=JstrR,JendR
        do i=IstrR,IendR
          ssh(i,j)=ETA0*sin(omega*time-kwave*(xr(i,j)))
        enddo
      enddo

# elif defined DUNE
      do j=JstrR,JendR
        do i=IstrR,IendR
          ssh(i,j)=0.0464+xr(i,j)*(0.0394-0.0464)/xl
        enddo
      enddo
# else
      do j=JstrR,JendR
        do i=IstrR,IendR
          ssh(i,j)=0.
        enddo
      enddo
# endif

# ifdef WET_DRY
!
!  Check that min depth > Dcrit
!
      do j=JstrR,JendR
        do i=IstrR,IendR
          if (ssh(i,j) .lt. Dcrit(i,j)-h(i,j)) then
            ssh(i,j)=Dcrit(i,j)-h(i,j)
          endif
        enddo
      enddo
# endif

# if defined EW_PERIODIC || defined NS_PERIODIC || defined  MPI
      call exchange_r2d_tile (istr,iend,jstr,jend, ssh(START_2D_ARRAY))
# endif
      return
      end
#endif /* ANA_SSH && ZCLIMATOLOGY */
!
!======================================================================
!                   subroutine ana_sst
!======================================================================
!
#if defined ANA_SST && defined QCORRECTION
      subroutine ana_sst_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  This routine sets sea surface temperature SST[Celsius] and surface
!  net heat flux sensitivity dQdSTT to sea surface temperature using
!  analytical expressions. dQdSTT is usually computed in units of
!  [Watts/m^2/degC]. It needs to be scaled to [m/s] by dividing by
!  rho0*Cp.  These forcing fields are used when the heat flux
!  correction is activated:
!
!       Q_model ~ Q + dQdSST * (T_model - SST)
!----------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "grid.h"
# include "forces.h"
# include "scalars.h"
      integer Istr,Iend,Jstr,Jend, i,j
# ifdef EQUATOR
      real y1,y2,sst1,sst2
# endif
!
# include "compute_auxiliary_bounds.h"
!
# ifdef EQUATOR
! SST = 25C and lineraly decreases to 10C 1200km from the Equator.
      y1=1200.E+3
      y2=1500.E+3
      sst1=10.
      sst2=25.
!$acc kernels if(compute_on_device) default(present)
      do j=JstrR,JendR
        do i=IstrR,IendR
          sst(i,j)=sst1
          if ((yr(i,j).gt.-y1).and.(yr(i,j).lt.y1)) then
            sst(i,j)=sst2
          else
            if ((yr(i,j).gt.-y2).and.(yr(i,j).lt.-y1)) then
              sst(i,j)=((sst2-sst1)*yr(i,j)-sst1*y1+y2*sst2)/(y2-y1)
            endif
            if ((yr(i,j).gt.y1).and.(yr(i,j).lt.y2)) then
              sst(i,j)=((sst2-sst1)*yr(i,j)+sst1*y1-y2*sst2)/(y1-y2)
            endif
          endif
          dqdt(i,j)=-50.0/(rho0*Cp)
        enddo
      enddo
!$acc end kernels      
# elif defined INNERSHELF
!$acc kernels if(compute_on_device) default(present)
      do j=JstrR,JendR
        do i=IstrR,IendR
          sst(i,j)=24.
          dqdt(i,j)=-50.0/(rho0*Cp)
        enddo
      enddo
!$acc end kernels
# else
!$acc kernels if(compute_on_device) default(present)
      do j=JstrR,JendR
        do i=IstrR,IendR
          sst(i,j)=25
          dqdt(i,j)=0.
        enddo
      enddo
!$acc end kernels
# endif
      return
      end
#endif /* ANA_SST && QCORRECTION */

!
!======================================================================
!                   subroutine ana_sss
!======================================================================
!
#if defined SALINITY && defined SFLX_CORR && defined ANA_SSS
      subroutine ana_sss_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  This routine sets sea surface salinity SSS[PSU] using
!  analytical expressions. This forcing field is used when the
!  salt flux correction is activated:
!
!    SSSFLX_model ~ SSS*(E-P) + CST * (SSS_model - SSS)
!
!  we use DQDSST for CST....
!
!----------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "grid.h"
# include "forces.h"
# include "scalars.h"
      integer Istr,Iend,Jstr,Jend, i,j
!
# include "compute_auxiliary_bounds.h"
!
      do j=JstrR,JendR
        do i=IstrR,IendR
          sss(i,j)=???
        enddo
      enddo
      return
      end
#endif /* SALINITY && SFLX_CORR && ANA_SSS */
!
!======================================================================
!                   subroutine ana_stflux
!======================================================================
!
#if defined ANA_STFLUX || defined ANA_SSFLUX
      subroutine ana_stflux_tile (Istr,Iend,Jstr,Jend, itrc)
!
!----------------------------------------------------------------------
!  This routine sets kinematic surface flux of tracer type variables
!  "stflx" (tracer units m/s) using analytical expressions.
!
!  On Input:
!     itrc      Tracer type array index.
!----------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "grid.h"
# include "forces.h"
# include "scalars.h"
      integer itrc, Istr,Iend,Jstr,Jend, i,j
c
# include "compute_auxiliary_bounds.h"
c
# ifdef TRACERS
!$acc kernels if(compute_on_device) default(present)
#  ifdef TEMPERATURE
      if (itrc.eq.itemp) then
!
!  Set kinematic surface heat flux [degC m/s] at horizontal
!  RHO-points.
!
        do j=JstrR,JendR
          do i=IstrR,IendR
#   if defined SINGLE_COLUMN && defined WILLIS_DEARDORFF
            stflx(i,j,itemp)= - 100./(rho0*Cp)
#   elif defined SINGLE_COLUMN && defined DIURNAL_CYCLE
            stflx(i,j,itemp)= srflx(i,j) - 75./(rho0*Cp)
#   else
            stflx(i,j,itemp)=0.
#   endif
          enddo
        enddo
      endif
!
#  endif /* TEMPERATURE */
#  ifdef SALINITY
      if (itrc.eq.isalt) then
!
!  Set kinematic surface freshwater flux (m/s) at horizontal
!  RHO-points, scaling by surface salinity is done in STEP3D.
!
        do j=JstrR,JendR
          do i=IstrR,IendR
            stflx(i,j,isalt)=0.
          enddo
        enddo
#   ifdef RAIN_FLUX  
        do j=JstrR,JendR
          do i=IstrR,IendR
            !set a value in cm/day , converson is done in m/s
            EmP(i,j)=-0.1*0.01/86400.
          enddo
        enddo
#   endif
      endif
#  endif /* SALINITY */
#  if defined SALINITY && defined TEMPERATURE
      if (itrc.ne.isalt.and.itrc.ne.itemp) then
#  endif /* SALINITY && TEMPERATURE */
#  if !defined SALINITY && defined TEMPERATURE
      if (itrc.ne.itemp) then
#  endif /* !SALINITY && TEMPERATURE */
#  if defined SALINITY && !defined TEMPERATURE
      if (itrc.ne.isalt) then
#  endif /* SALINITY && !TEMPERATURE */

!
!  Set kinematic surface flux of additional tracers, if any.
!
        do j=JstrR,JendR
          do i=IstrR,IendR
            stflx(i,j,itrc)=0.
          enddo
        enddo
#  if defined SALINITY || defined TEMPERATURE
      endif
#  endif /* SALINITY || TEMPERATURE*/
!$acc end kernels 
# endif /* TRACERS */
      return
      end
#endif /* ANA_STFLUX || ANA_SSFLUX */
!
!======================================================================
!                   subroutine ana_tclima
!======================================================================
!
#ifdef TCLIMATOLOGY
      subroutine ana_tclima (tile)
      implicit none
# include"param.h"
      integer tile
# ifdef ALLOW_SINGLE_BLOCK_MODE
C$    integer  trd, omp_get_thread_num
# endif
# include "compute_tile_bounds.h"
      call ana_tclima_tile   (Istr,Iend,Jstr,Jend)
      return
      end
      subroutine ana_tclima_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  This routine sets analytical ACTIVE (T&S) tracer climatology fields.
!----------------------------------------------------------------------
!
# ifdef SUBSTANCE
      USE comsubstance, ONLY: cobc_wat
# endif
      implicit none
# include "param.h"
# include "grid.h"
# include "climat.h"
# include "ocean3d.h"
# include "scalars.h"
# include "sediment.h"
      integer Istr,Iend,Jstr,Jend, i,j,k, itrc
      real    cff,cff1
# ifdef IGW
      real    sig,rhom,depth
# endif
# if defined SUBSTANCE || defined SEDIMENT
      integer iv
# endif
!
# include "compute_auxiliary_bounds.h"
!
# ifdef ANA_TCLIMA
#  ifdef JET
!
!        --> climatology computed in ana_initial
!
#  elif defined IGW
      cff=((2.e-3)**2)/g
      rhom=1027.037885
      do k=1,N
        do j=JstrR,JendR
          do i=IstrR,IendR
            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
            tclm(i,j,k,itemp)=(R0-sig)/Tcoef
#   elif defined EXPERIMENT2
            if (depth.le.50.) then
             sig = 25.8
            else
             sig = 27.0
            endif
            tclm(i,j,k,itemp)=(R0-sig)/Tcoef
#   else /* default EXPERIMENT3 */
            sig = 26.2 + rhom*(depth-50.)*cff
            tclm(i,j,k,itemp)=(R0-sig)/Tcoef
#   endif
          enddo
        enddo
      enddo
#  else
      do k=1,N
        do j=JstrR,JendR
          do i=IstrR,IendR
#   ifdef TEMPERATURE
            tclm(i,j,k,itemp)=t(i,j,k,1,itemp)
#   endif /* TEMPERATURE */
#   ifdef SALINITY
            tclm(i,j,k,isalt)=t(i,j,k,1,isalt)
#   endif /* SALINITY */
          enddo
        enddo
      enddo
#  endif /* JET */
# endif

# ifdef BIOLOGY
#  define temp cff
#  define SiO4 cff1
      do k=1,N
        do j=JstrR,JendR
          do i=IstrR,IendR
#  ifdef ANA_TCLIMA
#   ifdef TEMPERATURE
            temp=t(i,j,k,1,itemp)
#   else
            temp=25.
#   endif
            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
            tclm(i,j,k,iNO3_)=1.67+0.5873*SiO4+0.0144*SiO4**2
     &                            +0.0003099*SiO4**3
#   ifdef PISCES
            tclm(i,j,k,iDIC_)=2150.
            tclm(i,j,k,iTAL_)=2350.
            tclm(i,j,k,iOXY_)=200.
            tclm(i,j,k,iPOC_)=0.01
            tclm(i,j,k,iPHY_)=0.01
            tclm(i,j,k,iZOO_)=0.01
            tclm(i,j,k,iDOC_)=5.
            tclm(i,j,k,iFER_)=6.e-4
#    if ! defined key_pisces_npzd
            tclm(i,j,k,iCAL_)=0.01
            tclm(i,j,k,iPO4_)=tclm(i,j,k,iNO3_)/16.
            tclm(i,j,k,iSIL_)=91.51
            tclm(i,j,k,iDIA_)=0.01
            tclm(i,j,k,iMES_)=0.01
            tclm(i,j,k,iDSI_)=1.5e-3
            tclm(i,j,k,iBFE_)=1.e-2*5.e-6
            tclm(i,j,k,iGOC_)=0.01
            tclm(i,j,k,iSFE_)=1.e-2*5.e-6
            tclm(i,j,k,iDFE_)=1.e-2*5.e-6
            tclm(i,j,k,iGSI_)=1.e-2*0.15
            tclm(i,j,k,iNFE_)=1.e-2*5.e-6
            tclm(i,j,k,iNCH_)=1.e-2*12./55.
            tclm(i,j,k,iDCH_)=1.e-2*12./55.
            tclm(i,j,k,iNH4_)=1.e-2
#    endif
#    ifdef key_ligand
            tclm(i,j,k,iLGW_)=1.e-3
#    endif
#    ifdef key_pisces_quota
            tclm(i,j,k,iDON_)=5.
            tclm(i,j,k,iDOP_)=5.
            tclm(i,j,k,iPON_)=1.e-2
            tclm(i,j,k,iPOP_)=1.e-2
            tclm(i,j,k,iNPH_)=1.e-2
            tclm(i,j,k,iPPH_)=1.e-2
            tclm(i,j,k,iNDI_)=1.e-2
            tclm(i,j,k,iPDI_)=1.e-2
            tclm(i,j,k,iPIC_)=1.e-2
            tclm(i,j,k,iNPI_)=1.e-2
            tclm(i,j,k,iPPI_)=1.e-2
            tclm(i,j,k,iPFE_)=1.e-2*5.E-6
            tclm(i,j,k,iPCH_)=1.e-2*12./55.
            tclm(i,j,k,iGON_)=1.e-2
            tclm(i,j,k,iGOP_)=1.e-2
#    endif
#   elif defined BIO_NChlPZD
            tclm(i,j,k,iChla)=0.08
            tclm(i,j,k,iPhy1)=0.1
            tclm(i,j,k,iZoo1)=0.06
            tclm(i,j,k,iDet1)=0.02
#   elif defined BIO_N2ChlPZD2
            tclm(i,j,k,iNH4_)=0.1
            tclm(i,j,k,iChla)=0.08
            tclm(i,j,k,iPhy1)=0.06
            tclm(i,j,k,iZoo1)=0.04
            tclm(i,j,k,iDet1)=0.02
            tclm(i,j,k,iDet2)=0.02
#   elif defined BIO_BioEBUS
            tclm(i,j,k,iNO2_)=0.05*exp(z_r(i,j,k)/100.)
            tclm(i,j,k,iNH4_)=0.1*exp(z_r(i,j,k)/100.)
            tclm(i,j,k,iPhy1)=0.04*exp(z_r(i,j,k)/50.)
            tclm(i,j,k,iPhy2)=0.06*exp(z_r(i,j,k)/50.)
            tclm(i,j,k,iZoo1)=0.04*exp(z_r(i,j,k)/100.)
            tclm(i,j,k,iZoo2)=0.04*exp(z_r(i,j,k)/100.)
            tclm(i,j,k,iDet1)=0.02
            tclm(i,j,k,iDet2)=0.02
            tclm(i,j,k,iDON) =0.5*exp(z_r(i,j,k)/100.)
#    ifdef NITROUS_OXIDE
            tclm(i,j,k,iN2O)=-(0.008*exp(z_r(i,j,k)/100.)-0.01)
            if (tclm(i,j,k,iN2O).lt.0.) then
              tclm(i,j,k,iN2O)=0.
            endif
#    endif
#   endif
#  else
            if (.not.got_tclm(iNO3_)) then
              temp=t(i,j,k,1,itemp)
              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
              tclm(i,j,k,iNO3_)=1.67+0.5873*SiO4+0.0144*SiO4**2
     &                              +0.0003099*SiO4**3
            endif
#   ifdef PISCES
            if (.not.got_tclm(iDIC_)) tclm(i,j,k,iDIC_)=2150.
            if (.not.got_tclm(iTAL_)) tclm(i,j,k,iTAL_)=2350.
            if (.not.got_tclm(iOXY_)) tclm(i,j,k,iOXY_)=200.
            if (.not.got_tclm(iPOC_)) tclm(i,j,k,iPOC_)=0.01
            if (.not.got_tclm(iPHY_)) tclm(i,j,k,iPHY_)=0.01
            if (.not.got_tclm(iZOO_)) tclm(i,j,k,iZOO_)=0.01
            if (.not.got_tclm(iDOC_)) tclm(i,j,k,iDOC_)=5.
            if (.not.got_tclm(iFER_)) tclm(i,j,k,iFER_)=6.e-4
#    if ! defined key_pisces_npzd
            if (.not.got_tclm(iCAL_)) tclm(i,j,k,iCAL_)=0.01
            if (.not.got_tclm(iPO4_)) then
#     ifdef TEMPERATURE
              temp=t(i,j,k,1,itemp)
#     else
            temp=25.
#     endif
              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
              tclm(i,j,k,iPO4_)=(1.67+0.5873*SiO4+0.0144*SiO4**2
     &                               +0.0003099*SiO4**3)/16.
            endif
            if (.not.got_tclm(iSIL_)) tclm(i,j,k,iSIL_)=91.51
            if (.not.got_tclm(iDIA_)) tclm(i,j,k,iDIA_)=0.01
            if (.not.got_tclm(iMES_)) tclm(i,j,k,iMES_)=0.01
            if (.not.got_tclm(iDSI_)) tclm(i,j,k,iDSI_)=0.0015
            if (.not.got_tclm(iBFE_)) tclm(i,j,k,iBFE_)=5.e-8
            if (.not.got_tclm(iGOC_)) tclm(i,j,k,iGOC_)=0.01
            if (.not.got_tclm(iSFE_)) tclm(i,j,k,iSFE_)=5.e-8
            if (.not.got_tclm(iDFE_)) tclm(i,j,k,iDFE_)=5.e-8
            if (.not.got_tclm(iGSI_)) tclm(i,j,k,iGSI_)=0.0015
            if (.not.got_tclm(iNFE_)) tclm(i,j,k,iNFE_)=5.e-8
            if (.not.got_tclm(iNCH_)) tclm(i,j,k,iNCH_)=1.e-2*12./55.
            if (.not.got_tclm(iDCH_)) tclm(i,j,k,iDCH_)=1.e-2*12./55.
            if (.not.got_tclm(iNH4_)) tclm(i,j,k,iNH4_)=0.01
#    endif
#    ifdef key_ligand
            if (.not.got_tclm(iLGW_)) tclm(i,j,k,iLGW_)=1.e-3
#    endif
#    ifdef key_pisces_quota
            if (.not.got_tclm(iDON_)) tclm(i,j,k,iDON_)=5.
            if (.not.got_tclm(iDOP_)) tclm(i,j,k,iDOP_)=5.
            if (.not.got_tclm(iPON_)) tclm(i,j,k,iPON_)=0.01
            if (.not.got_tclm(iPOP_)) tclm(i,j,k,iPOP_)=0.01
            if (.not.got_tclm(iNPH_)) tclm(i,j,k,iNPH_)=0.01
            if (.not.got_tclm(iPPH_)) tclm(i,j,k,iPPH_)=0.01
            if (.not.got_tclm(iNDI_)) tclm(i,j,k,iNDI_)=0.01
            if (.not.got_tclm(iPDI_)) tclm(i,j,k,iPDI_)=0.01
            if (.not.got_tclm(iPIC_)) tclm(i,j,k,iPIC_)=0.01
            if (.not.got_tclm(iNPI_)) tclm(i,j,k,iNPI_)=0.01
            if (.not.got_tclm(iPPI_)) tclm(i,j,k,iPPI_)=0.01
            if (.not.got_tclm(iPFE_)) tclm(i,j,k,iPFE_)=0.01*5.e-6
            if (.not.got_tclm(iPCH_)) tclm(i,j,k,iPCH_)=0.01*12./55.
            if (.not.got_tclm(iGON_)) tclm(i,j,k,iGON_)=0.01
            if (.not.got_tclm(iGOP_)) tclm(i,j,k,iGOP_)=0.01
#    endif
#   elif defined BIO_NChlPZD
            if (.not.got_tclm(iChla)) tclm(i,j,k,iChla)=0.08
            if (.not.got_tclm(iPhy1)) tclm(i,j,k,iPhy1)=0.1
            if (.not.got_tclm(iZoo1)) tclm(i,j,k,iZoo1)=0.06
            if (.not.got_tclm(iDet1)) tclm(i,j,k,iDet1)=0.02
#   elif defined BIO_N2ChlPZD2
            if (.not.got_tclm(iNH4_)) tclm(i,j,k,iNH4_)=0.1
            if (.not.got_tclm(iChla)) tclm(i,j,k,iChla)=0.08
            if (.not.got_tclm(iPhy1)) tclm(i,j,k,iPhy1)=0.06
            if (.not.got_tclm(iZoo1)) tclm(i,j,k,iZoo1)=0.04
            if (.not.got_tclm(iDet1)) tclm(i,j,k,iDet1)=0.02
            if (.not.got_tclm(iDet2)) tclm(i,j,k,iDet2)=0.02
#   elif defined BIO_BioEBUS
            if (.not.got_tclm(iNO2_)) tclm(i,j,k,iNO2_)=0.05*
!            if (.not.got_tclm(iNO2_)) tclm(i,j,k,iNO2_)=0.6*
     &                                  exp(z_r(i,j,k)/100.)
            if (.not.got_tclm(iNH4_)) tclm(i,j,k,iNH4_)=0.1*
     &                                  exp(z_r(i,j,k)/100.)
            if (.not.got_tclm(iPhy1)) tclm(i,j,k,iPhy1)=0.04*
     &                                  exp(z_r(i,j,k)/50.)
            if (.not.got_tclm(iPhy2)) tclm(i,j,k,iPhy2)=0.06*
     &                                  exp(z_r(i,j,k)/50.)
            if (.not.got_tclm(iZoo1)) tclm(i,j,k,iZoo1)=0.04*
     &                                  exp(z_r(i,j,k)/100.)
            if (.not.got_tclm(iZoo2)) tclm(i,j,k,iZoo2)=0.04*
     &                                  exp(z_r(i,j,k)/100.)
            if (.not.got_tclm(iDet1)) tclm(i,j,k,iDet1)=0.02
            if (.not.got_tclm(iDet2)) tclm(i,j,k,iDet2)=0.02
            if (.not.got_tclm(iDON))  tclm(i,j,k,iDON)=0.5*
!            if (.not.got_tclm(iDON))  tclm(i,j,k,iDON)=2.*
     &                                  exp(z_r(i,j,k)/100.)
            if (.not.got_tclm(iO2)) tclm(i,j,k,iO2)=250.*
     &                                  exp(z_r(i,j,k)/100.)
#    ifdef CARBON
            if (.not.got_tclm(iTCO2)) tclm(i,j,k,iTCO2)=2100.
            if (.not.got_tclm(iTALK)) tclm(i,j,k,iTALK)=2300.
#    endif
#    ifdef NITROUS_OXIDE
            if (.not.got_tclm(iN2O)) then
              tclm(i,j,k,iN2O)=-(0.008*exp(z_r(i,j,k)/100.)-0.01)
!	       tclm(i,j,k,iN2O)=-(0.01*exp(z_r(i,j,k)/100.)-0.01)
              if (tclm(i,j,k,iN2O).lt.0.) then
                tclm(i,j,k,iN2O)=0.
              endif
            endif
#    endif
#    ifdef HYDROGEN_SULFIDE
            if (.not.got_tclm(iH2S)) tclm(i,j,k,iH2S)=0.
#    endif
#   endif
#  endif /* ANA_TCLIMA */
          enddo
        enddo
      enddo
#  undef SiO4
#  undef temp
# endif /* BIOLOGY */

# ifdef SEDIMENT
      do itrc=1,NST
        iv=itrc_sed+itrc-1
        do k=1,N
          do j=JstrR,JendR
            do i=IstrR,IendR
#  ifdef ANA_TCLIMA
              tclm(i,j,k,iv)=0.   !Csed(itrc)
#  else
              if (.not.got_tclm(iv)) then
                tclm(i,j,k,iv)=0. !Csed(itrc)
              endif
#  endif
            enddo
          enddo
        enddo
      enddo
# endif /* SEDIMENT */

# ifdef PASSIVE_TRACER
      do k=1,N
        do j=JstrR,JendR
          do i=IstrR,IendR
#  ifdef ANA_TCLIMA
            tclm(i,j,k,itpas)=0.0
#  else
            if (.not.got_tclm(itpas)) then
              tclm(i,j,k,itpas)=0.0
            endif
#  endif
          enddo
        enddo
      enddo
# endif

# ifdef SUBSTANCE
      do k=1,N
        do j=JstrR,JendR
          do i=IstrR,IendR
            do iv=itsubs1,itsubs2
#  ifdef ANA_TCLIMA
              tclm(i,j,k,iv)=cobc_wat(iv)
#  else
              if (.not.got_tclm(iv)) then
                tclm(i,j,k,iv)=cobc_wat(iv)
              endif
#  endif
            enddo
          enddo
        enddo
      enddo
# endif

# if defined BD_EXCHANGE
      do itrc=1,NT
#  ifndef ANA_TCLIMA
        if (.not.got_tclm(itrc)) then
#  endif
#  ifdef THREE_GHOST_POINTS_TS
          call exchange_r3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &                                 tclm(START_2D_ARRAY,1,itrc))
#  else
          call exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &                            tclm(START_2D_ARRAY,1,itrc))
#  endif
#  ifndef ANA_TCLIMA
        endif
#  endif
      enddo
# endif

      return
      end
#endif /* TCLIMATOLOGY */
!
!======================================================================
!                   subroutine ana_uclima
!======================================================================
!
#if defined ANA_M2CLIMA && defined M2CLIMATOLOGY ||\
   (defined ANA_M3CLIMA && defined M3CLIMATOLOGY)
      subroutine ana_uclima (tile)
      implicit none
# include "param.h"
      integer tile
# ifdef ALLOW_SINGLE_BLOCK_MODE
C$    integer  trd, omp_get_thread_num
# endif
# include "compute_tile_bounds.h"
      call ana_uclima_tile   (Istr,Iend,Jstr,Jend)
      return
      end
!
      subroutine ana_uclima_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  This routine sets analytical momentum climatology fields.
!----------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "grid.h"
# include "ocean2d.h"
# include "ocean3d.h"
# include "climat.h"
# include "scalars.h"
      integer Istr,Iend,Jstr,Jend, i,j,k
# if defined INTERNAL || defined ANA_TIDES
      real U0,omega,kwave,V0
# endif
      real ramp
!
# include "compute_auxiliary_bounds.h"
!
# ifdef EW_PERIODIC
#  define IU_RANGE Istr,Iend
#  define IV_RANGE Istr,Iend
# else
#  define IU_RANGE Istr,IendR
#  define IV_RANGE IstrR,IendR
# endif

# ifdef NS_PERIODIC
#  define JU_RANGE Jstr,Jend
#  define JV_RANGE Jstr,Jend
# else
#  define JU_RANGE JstrR,JendR
#  define JV_RANGE Jstr,JendR
# endif
!
# if defined ANA_M2CLIMA && defined M2CLIMATOLOGY
#  if defined INTERNAL || (defined ANA_TIDES && !defined RIP  && !defined TIDAL_FLAT && !defined ESTUARY)
      U0=0.02
      omega=2.*pi/(12.4*3600)
      kwave=sqrt(((omega*omega)-(f(1,1)*f(1,1)))/(g*h(1,1)))
      V0=f(1,1)*U0/omega
      do j=JU_RANGE
        do i=IU_RANGE
          ubclm(i,j)=U0*sin(omega*time-kwave*0.5*(xr(i,j)+xr(i-1,j)))
        enddo
      enddo
      do j=JV_RANGE
        do i=IV_RANGE
          vbclm(i,j)=V0*cos(omega*time-kwave*0.5*(xr(i,j)+xr(i,j-1)))
        enddo
      enddo
#  elif defined JET
!
!       --> climatology computed in ana_initial
!
#  elif defined DUNE
      do j=JstrR,JendR
        do i=IstrR,IendR
#   ifdef ANA_DUNE
          ubclm(i,j)=1.67
#   else
          ubclm(i,j)=0.58
#   endif
          vbclm(i,j)=0.
        enddo
      enddo
#  else
      do j=JstrR,JendR
        do i=IstrR,IendR
          ubclm(i,j)=0.
          vbclm(i,j)=0.
        enddo
      enddo
#  endif
#  if defined BD_EXCHANGE
      call exchange_u2d_tile (Istr,Iend,Jstr,Jend,  ubclm)
      call exchange_v2d_tile (Istr,Iend,Jstr,Jend,  vbclm)
#  endif
# endif

# if defined ANA_M3CLIMA && defined M3CLIMATOLOGY && defined SOLVE3D
#  ifdef JET
!
!       --> climatology computed in ana_initial
!
#  else
      do k=1,N
        do j=JstrR,JendR
          do i=IstrR,IendR
            uclm(i,j,k)=0.
            vclm(i,j,k)=0.
          enddo
        enddo
      enddo
#  endif
#  if defined BD_EXCHANGE
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,  uclm)
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,  vclm)
#  endif
# endif
# undef IU_RANGE
# undef JU_RANGE
# undef IV_RANGE
# undef JV_RANGE
      return
      end
#endif /* ANA_M2CLIMA && M2CLIMATOLOGY || (ANA_M3CLIMA && M3CLIMATOLOGY) */
!
!======================================================================
!                   subroutine ana_nbq_clima
!======================================================================
!
#ifdef NBQCLIMATOLOGY
      subroutine ana_nbq_clima (tile)
      implicit none
# include "param.h"
      integer tile
# ifdef ALLOW_SINGLE_BLOCK_MODE
C$    integer  trd, omp_get_thread_num
# endif
# include "compute_tile_bounds.h"
      call ana_nbq_clima_tile (Istr,Iend,Jstr,Jend)
      return
      end
!
      subroutine ana_nbq_clima_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  This routine sets analytical NBQ climatology fields.
!----------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "grid.h"
# include "ocean3d.h"
# include "climat.h"
# include "scalars.h"
      integer Istr,Iend,Jstr,Jend, i,j,k
!
# include "compute_auxiliary_bounds.h"
!
# ifdef EW_PERIODIC
#  define IU_RANGE Istr,Iend
#  define IV_RANGE Istr,Iend
# else
#  define IU_RANGE Istr,IendR
#  define IV_RANGE IstrR,IendR
# endif

# ifdef NS_PERIODIC
#  define JU_RANGE Jstr,Jend
#  define JV_RANGE Jstr,Jend
# else
#  define JU_RANGE JstrR,JendR
#  define JV_RANGE Jstr,JendR
# endif
!

# if defined NBQ && (defined NBQCLIMATOLOGY || \
                    (defined AGRIF && !defined NBQ_FRC_BRY)) \
                 && (defined M3CLIMATOLOGY || \
                    (defined AGRIF && !defined M3_FRC_BRY))
      do k=1,N
        do j=JU_RANGE
          do i=IU_RANGE
            unbqclm(i,j,k)=0.5*(Hz(i,j,k)+Hz(i-1,j,k))*u(i,j,k,nrhs)
          enddo
        enddo
        do j=JV_RANGE
          do i=IV_RANGE
            vnbqclm(i,j,k)=0.5*(Hz(i,j,k)+Hz(i,j-1,k))*v(i,j,k,nrhs)
          enddo
        enddo
#  ifdef NBQ
        do j=JstrR,JendR
          do i=IstrR,IendR
            wnbqclm(i,j,k)=0.
            rnbqclm(i,j,k)=0.
          enddo
        enddo
#  endif
      enddo
# else
      do k=1,N
        do j=JU_RANGE
          do i=IU_RANGE
            unbqclm(i,j,k)=0.
          enddo
        enddo
        do j=JV_RANGE
          do i=IV_RANGE
            vnbqclm(i,j,k)=0.
          enddo
        enddo
#  ifdef NBQ
        do j=JstrR,JendR
          do i=IstrR,IendR
            wnbqclm(i,j,k)=0.
            rnbqclm(i,j,k)=0.
          enddo
        enddo
#  endif
      enddo
# endif

# if defined BD_EXCHANGE
      call exchange_u3d_tile (Istr,Iend,Jstr,Jend,  unbqclm)
      call exchange_v3d_tile (Istr,Iend,Jstr,Jend,  vnbqclm)
#  ifdef NBQ
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,  wnbqclm)
      call exchange_r3d_tile (Istr,Iend,Jstr,Jend,  rnbqclm)
#  endif
# endif
# undef IU_RANGE
# undef JU_RANGE
# undef IV_RANGE
# undef JV_RANGE
      return
      end
#endif /* NBQCLIMATOLOGY */
!
!======================================================================
!                   subroutine ana_wwave
!======================================================================
!
#ifdef ANA_WWAVE
      subroutine ana_wwave (tile)
      implicit none
# include "param.h"
      integer tile
# ifdef ALLOW_SINGLE_BLOCK_MODE
C$    integer  trd, omp_get_thread_num
# endif
# include "compute_tile_bounds.h"
      call ana_wwave_tile   (Istr,Iend,Jstr,Jend)
      return
      end
!
      subroutine ana_wwave_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  This routine sets wind induced wave amplitude, direction
!  and period used in the bottom boundary layer formulation.
!----------------------------------------------------------------------
!
      implicit none
      real wamp,wdir,wprd,
     &     ho,khd,kh,kw,kr,ks,wday,ramp,cff,co,gamw,
     &     cff1,cff2,Hrms,cfrq,cw,ch,cr,sbc,dd,cgo,
     &     dsup,Btg,cosw,sinw,eps,nw,wh
      parameter (eps=1.0e-10)
# include "param.h"
# include "grid.h"
# include "forces.h"
# include "ocean3d.h"
# include "scalars.h"
# ifdef WKB_WWAVE
#  include "wkb_wwave.h"
# endif
      integer Istr,Iend,Jstr,Jend, i,j,k
!
# include "compute_extended_bounds.h"
!
!  Set wind induced wave amplitude (m), direction (radians) and
!  period (s) at RHO-points.
!
# ifdef SHOREFACE
      do j=JstrR,JendR
        do i=IstrR,IendR
          Awave(i,j)=0.707           ! wave amplitude (m)
          Dwave(i,j)=-10.*deg2rad    ! wave direction (rad)
          Pwave(i,j)=10.             ! wave period (s)
        enddo
      enddo
      wamp = Awave(i,j)
      wdir = Dwave(i,j)
      wprd = Pwave(i,j)
      cfrq = 2.0*pi/wprd    ! peak wave frequency (rad/s)
      ho=h(1,1)             ! offshore depth
#  ifdef OBC_EAST
#   ifdef MPI
        ho=h(LLm/NP_XI,1)
        call MPI_Bcast(ho, 1, MPI_DOUBLE_PRECISION,
     &                     NP_XI-1, MPI_COMM_WORLD, ierr)
#   else
        ho=h(LLm,1)
#   endif
#  endif
      khd = ho*cfrq*cfrq/g
      kh  = sqrt(    khd*khd + khd/(1.0 + khd*(0.6666666666
     &               +khd*(0.3555555555 + khd*(0.1608465608
     &               +khd*(0.0632098765 + khd*(0.0217540484
     &                               +khd*0.0065407983)))))) )
      co=sqrt(g/kh*ho*tanh(kh))
      cgo=co*0.5*(1.+2.*kh/sinh(2.*kh))
      do j=jstrR,jendR
        do i=istrR,iendR
          dsup=0.
          do k=1,4
            dd = h(i,j) + z_w(i,j,N) + dsup
            khd = dd*cfrq*cfrq/g
            kh = sqrt(    khd*khd + khd/(1.0 + khd*(0.6666666666
     &                    +khd*(0.3555555555 + khd*(0.1608465608
     &                    +khd*(0.0632098765 + khd*(0.0217540484
     &                                    +khd*0.0065407983)))))) )
            kw=kh/max(dd,eps)
            cw=sqrt(g/kw*tanh(kh))
            nw=0.5*(1.+2.*kh/sinh(2.*kh))  ! n=Cg/C
            ks=sqrt(cgo/(2.*nw*cw))        ! shoaling coefficient
            kr=sqrt(cos(wdir)/cos(wdir))   ! refraction coefficient
            cosw=cos(wdir)
            sinw=sin(wdir)
            gamw=0.55
            cff1=gamw*dd
            cff2=2.0*wamp*ks*kr
            wh=min(cff1,cff2)
            dsup=-0.125*(wh**2)*kw/sinh(2.*kw*dd) ! wave set-up
          enddo
!
!  Fill global arrays
!
          whrm(i,j)=wh
          wfrq(i,j)=sqrt(g*kw*tanh(kw*dd))
          wdrx(i,j)=cosw
          wdre(i,j)=sinw
#  ifdef WET_DRY
          if (h(i,j)+z_w(i,j,N).le.Dcrit(i,j)+eps) then
            whrm(i,j)=0.
            wfrq(i,j)=0.
            wdrx(i,j)=0.
            wdre(i,j)=0.
          endif
#  endif
        enddo
      enddo

# else
      do j=JstrR,JendR
        do i=IstrR,IendR
          Awave(i,j)=0.707            ! wave amplitude (m)
          Dwave(i,j)=-10.*deg2rad     ! wave direction (rad)
          Pwave(i,j)=10.              ! wave period (s)
        enddo
      enddo
#  ifdef MRL_WCI
      Hrms=2.*Awave(1,1)         ! offshore RMS wave height (m)
      cfrq=2.*pi/Pwave(1,1)      ! peak wave frequency (rad/s)
      wdir=Dwave(1,1)            ! wave direction (rad)
!
!  Fill global arrays
!
      do j=jstr-1,jend+1
        do i=istr-1,iend+1
          whrm(i,j)=Hrms
          wfrq(i,j)=cfrq
          wdrx(i,j)=cos(wdir)
          wdre(i,j)=sin(wdir)
!  limit Hrms in surfzone
          gamw=0.55
          dd=h(i,j)+z_w(i,j,N)
          whrm(i,j)=min(whrm(i,j),gamw*dd)
#   ifdef WET_DRY
          if (h(i,j)+z_w(i,j,N) .lt. Dcrit(i,j)+eps) then
            whrm(i,j)=0.
            wfrq(i,j)=0.
            wdrx(i,j)=0.
            wdre(i,j)=0.
          endif
#   endif
        enddo
      enddo
#  endif /* MRL_WCI */
# endif /* SHOREFACE or else */
      return
      end
#endif /* ANA_WWAVE */
!
!======================================================================
!                   subroutine ana_sediment
!======================================================================
!
#ifdef SEDIMENT
      subroutine ana_sediment (tile)
      implicit none
# include "param.h"
      integer tile
# ifdef ALLOW_SINGLE_BLOCK_MODE
C$    integer  trd, omp_get_thread_num
# endif
# include "compute_tile_bounds.h"
      call ana_sediment_tile (Istr,Iend,Jstr,Jend)
      return
      end
!
      subroutine ana_sediment_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  This routine sets sediment ripple and bed parameters
!  from values found in sediment.in:
!
!    Hrip    initial ripple height [m]
!    Lrip    initial ripple length [m]
!    Bthk    initial thicknesses of bed layers [m]
!    Bpor    initial porosity of bed layers [m]
!    Bfrac   volume fraction of each size class in each bed layer
!
!----------------------------------------------------------------------
!
      implicit none
# include "param.h"
# include "grid.h"
# include "scalars.h"
# include "sediment.h"
# include "bbl.h"
      real    cff1,cff2,cff3,cff4,cff5
      integer Istr,Iend,Jstr,Jend, i,j, k,ilay, ised
      real bmz(NLAY)
      real alpha, bzactv, frt, tcb_top, tcb_bot
      real tcr(NLAY)
!
# undef DEBUG
# ifdef DEBUG
      integer ick,jck
      parameter (ick=60, jck=35)
# endif

# include "compute_extended_bounds.h"



      do j=JstrR,JendR
        do i=IstrR,IendR
# ifdef BBL
          if(.not.got_inibed(1)) Hripple(i,j)=Hrip
          if(.not.got_inibed(2)) Lripple(i,j)=Lrip
# endif
          do ilay=1,NLAY
            if(.not.got_inised(1)) bed_thick(i,j,ilay)=Bthk(ilay)
            if(.not.got_inised(2)) bed_poros(i,j,ilay)=Bpor(ilay)
            if(.not.got_inised(3)) then
              do ised=1,NST
                bed_frac(i,j,ilay,ised)=Bfr(ilay,ised)
              enddo
            endif
            bed_age(i,j,ilay)=time
!            bot_thick(i,j) = 0.003
            bot_thick(i,j) = 0.00

            do ised=1,NST
              bed_mass(i,j,ilay,1,ised)=
     &                     bed_thick(i,j,ilay)*
     &                     Srho(ised)*
     &                     (1.0-bed_poros(i,j,ilay))*
     &                     bed_frac(i,j,ilay,ised)
# ifdef MASKING
     &                                 *rmask(i,j)
# endif
              bed_mass(i,j,ilay,2,ised)=bed_mass(i,j,ilay,1,ised)
            enddo
          enddo
          cff1=1.0
          cff2=1.0
          cff3=1.0
          cff4=1.
          do ised=1,NST
            cff1=cff1*tau_ce(ised)**bed_frac(i,j,1,ised)
            cff2=cff2*Sd(ised)**bed_frac(i,j,1,ised)
            cff3=cff3*wsed(ised)**bed_frac(i,j,1,ised)
            cff4=cff4*Srho(ised)**bed_frac(i,j,1,ised)
          enddo
          taucb(i,j)=cff1                ! [m2/s2]
          Ssize(i,j)=cff2
          w_set(i,j)=cff3                ! [m/s]
          Sdens(i,j)=MAX(cff4,1050.)
!          cff5=0
!          do ilay=1,NLAY
!            do ised=1,2
!              cff5 = cff5 + bed_mass(i,j,ilay,1,ised)
!            enddo
!            bed_bctr(i,j,ilay) = exp ( (log(cff5) - tcr_off ) /tcr_slp   )
!          enddo
# ifdef DEBUG
          if(j.eq.jck.and.i.eq.ick) then
            write(6,*) '********** ANA_SEDIMENT ***********'
            do ised=1,NST
              write(6,*) 'Sd(ised)',Sd(ised)
              do ilay=1,NLAY
                write(6,*) 'i,j,ilay,ised',i,j,ilay,ised
                write(6,*) 'bed_frac',bed_frac(i,j,ilay,ised)
              enddo
            enddo
            do ilay=1,NLAY
              write(6,*) 'bed_thick',bed_thick(i,j,ilay)
              write(6,*) 'bed_por',bed_poros(i,j,ilay)
            enddo
          endif
# endif
# ifdef MORPHODYN
          bed_thick_tot(i,j,2)=0.
          do ilay=1,NLAY
            bed_thick_tot(i,j,2)=bed_thick_tot(i,j,2)+
     &                           bed_thick(i,j,ilay)
          enddo
          bed_thick_tot(i,j,1)=bed_thick_tot(i,j,2)
# endif
        enddo
      enddo

# if defined COHESIVE_BED || defined MIXED_BED

      do j=JstrR,JendR
        do i=IstrR,IendR
          bed_bctr(i,j,1) =  0.1
      !
!  Key to this algorithm: "mass available depth" (mad) [kg/m2]
!  present tau crit profile (tcp)
!  representative tau crit profile (tcr)
!  Compute depth and mass depth of sediment column
!  Compute cumulative depth
!  Compute mass depth
          bmz(1) = 0.0
          do ised=1,NST
            bmz(1) = bmz(1)+bed_mass(i,j,1,1,ised)
          enddo
          do k=2,NLAY
            bmz(k) = bmz(k-1)
            do ised=1,NST
              bmz(k)=bmz(k)+bed_mass(i,j,k,1,ised)
            enddo
          enddo

!  Calculate representative critical shear stress profile
!  Note that the values are for the TOP of the layer...we
!  assume the bottom of the bottom layer has tcr = tcr_max
          tcr(1) =  0.1
          do k=2,NLAY
            tcr(k) = tcr_min
            if (bmz(k-1).GT.1.e-14) then
!               tcr(k) = exp((log(bmz(k-1))-                            &
!     &               bottom(i,j,idoff))/                                &
!     &               bottom(i,j,idslp))
              tcr(k) = exp((log(bmz(k-1))-
     &                      tcr_off)/
     &                      tcr_slp)
            endif
          !   tcr(k) = MIN( MAX( tcr(k), tcr_min), tcr_max )
          enddo
!          ks=Ksed-2
!          do k=NLAY,NLAY-ks,-1
!             taucb(i,j,k)=tcr(k)
!          enddo
!  Relax tau crit profile taucb(i,j,k)
!  towards representative profile tcr...100 x slower for "swelling"
!  than consolidation
!  TODO - make the factor an input parameter
!           if( bed_bctr(i,j,1).LE.tcr(1)) then
!             alpha = MIN(dt/bottom(i,j,idtim),1.0)
!             alpha = MIN(dt/tcr_tim,1.0)
!           else
!             alpha = MIN(dt/(100.0*bottom(i,j,idtim)),1.0)
!             alpha = MIN(dt/(100.0*tcr_tim),1.0)
!           endif
          do k=1,NLAY
!            bed_bctr(i,j,k)=bed_bctr(i,j,k)+
! &            alpha*(tcr(k)-bed_bctr(i,j,k))
            bed_bctr(i,j,k)=
     &          MIN( MAX( tcr(k), tcr_min), tcr_max )
          enddo
!ALA    Enforce last layer = tcr_max
!ALA          bed_bctr(i,j,NLAY)= tcr_max
        enddo
      enddo

#  if defined SED_TOY_RESUSP || defined SED_TOY_CONSOLID
      bed_bctr(:,:,:) = 5.2
      bed_bctr(:,:,1) = 0.03
      bed_bctr(:,:,2) = 0.0911460784453711
      bed_bctr(:,:,3) = 0.918694902789539
      bed_bctr(:,:,4) = 3.54929559677213
#  endif

# endif

# if defined BD_EXCHANGE
      do ised=1,NST
        do ilay=1,NLAY
          call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_frac(START_2D_ARRAY,ilay,ised))
          call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_mass(START_2D_ARRAY,ilay,1,ised))
          call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_mass(START_2D_ARRAY,ilay,2,ised))
         enddo
      enddo
      do ilay=1,NLAY
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_poros(START_2D_ARRAY,ilay))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_thick(START_2D_ARRAY,ilay))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_age(START_2D_ARRAY,ilay))
#  if defined COHESIVE_BED || defined MIXED_BED
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_bctr(START_2D_ARRAY,ilay))
#  endif
      enddo
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bot_thick(START_2D_ARRAY))
#  ifdef MORPHODYN
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_thick_tot(START_2D_ARRAY,1))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_thick_tot(START_2D_ARRAY,2))
#  endif
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 taucb(START_2D_ARRAY))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 Ssize(START_2D_ARRAY))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 w_set(START_2D_ARRAY))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 Sdens(START_2D_ARRAY))
# endif
      return
      end
#endif /* SEDIMENT */
!
!======================================================================
!                   subroutine ana_psource
!======================================================================
!
#if (defined PSOURCE || defined PSOURCE_MASS) && defined ANA_PSOURCE
!
!----------------------------------------------------------------------
!  Set analytical tracer and mass point sources and sinks
!----------------------------------------------------------------------
!
      subroutine ana_psource_tile (Istr,Iend,Jstr,Jend)
      implicit none
# include "param.h"
# include "scalars.h"
# include "sources.h"
# include "ocean3d.h"
# include "grid.h"
!
      integer is, k, Istr,Iend,Jstr,Jend, i,j,itrc,ised
      real    cff, cff1, cff2, ramp
      real    xno3,temp,SiO4,zsrc
# include "compute_auxiliary_bounds.h"

!$acc kernels if(compute_on_device) default(present)
      if (iic.eq.ntstart) then
!
! Set-up nondimensional shape Qshape, must add to unity
!
# ifdef RIVER

#  ifdef PSOURCE_MASS
     !if the source is applied as volume vertical influx
     !the distribution must be either uniform or located at the bottom
     !if you want a bottom influx you should use GLS_MIXING instead of LMD
#   define CST_SHAPE
!#    define BOTTOM_SHAPE
#  else
#   define EXP_SHAPE
!#   define CST_SHAPE
#  endif

#  ifdef CST_SHAPE
        cff=1./float(N)
        do k=1,N                       ! Uniform vertical
          do is=1,Nsrc                 ! distribution
            Qshape(is,k)=cff
          enddo
        enddo
#  elif defined BOTTOM_SHAPE
        do is=1,Nsrc
         Qshape(is,1)=1
        enddo
        do k=2,N                       ! All flux at bottom
          do is=1,Nsrc
            Qshape(is,k)=0.
          enddo
        enddo
#  elif defined EXP_SHAPE
        do is=1,Nsrc
#   ifdef MPI
          i=Isrc_mpi(is,mynode)
          j=Jsrc_mpi(is,mynode)
#   else
          i=Isrc(is)
          j=Jsrc(is)
#   endif
#   ifdef MPI
          if (IstrR.le.i .and. i.le.IendR .and.
     &       JstrR.le.j .and. j.le.JendR) then
#   endif
             call get_qshape_exp(Qshape(is,:), h(i,j), hinv(i,j))
#   ifdef MPI 
          endif
#   endif            
        enddo
#  elif defined AL_SHAPE
        do is=1,Nsrc                   ! Set-up nondimensional shape
          do k=1,10
            Qshape(is,k)=0.0
          enddo
          do k=11,14
            Qshape(is,k)=0.05
          enddo
          Qshape(is,15)=0.1           ! These most add to unity!
          Qshape(is,16)=0.1           ! These most add to unity!
          Qshape(is,17)=0.1           ! These most add to unity!
          Qshape(is,18)=0.1           ! These most add to unity!
          Qshape(is,19)=0.2
          Qshape(is,20)=0.2
        enddo
#  endif

# elif defined SEAGRASS

        cff=1./float(N)
        do k=1,N                       ! Uniform vertical
          do is=1,Nsrc                 ! distribution
            Qshape(is,k)=cff
          enddo
        enddo

# elif defined REGIONAL
#  define EXP_SHAPE
#  ifdef CST_SHAPE
        cff=1./float(N)
        do k=1,N                       ! Uniform vertical
          do is=1,Nsrc                 ! distribution
            Qshape(is,k)=cff
          enddo
        enddo
#  elif defined EXP_SHAPE
        do is=1,Nsrc
#   ifdef MPI
          i=Isrc_mpi(is,mynode)
          j=Jsrc_mpi(is,mynode)
#   else
          i=Isrc(is)
          j=Jsrc(is)
#   endif
#   ifdef MPI
          if (IstrR.le.i .and. i.le.IendR .and.
     &       JstrR.le.j .and. j.le.JendR) then
#   endif
             call get_qshape_exp(Qshape(is,:), h(i,j), hinv(i,j))
#   ifdef MPI
          endif
#   endif 
        enddo
#  endif

# else
#  define EXP_SHAPE
#  ifdef CST_SHAPE
        cff=1./float(N)
        do k=1,N                       ! Uniform vertical
          do is=1,Nsrc                 ! distribution
            Qshape(is,k)=cff
          enddo
        enddo
#  elif defined EXP_SHAPE
        do is=1,Nsrc
#   ifdef MPI
          i=Isrc_mpi(is,mynode)
          j=Jsrc_mpi(is,mynode)
#   else
          i=Isrc(is)
          j=Jsrc(is)
#   endif
#   ifdef MPI
          if (IstrR.le.i .and. i.le.IendR .and.
     &        JstrR.le.j .and. j.le.JendR) then
#   endif
            call get_qshape_exp(Qshape(is,:), h(i,j), hinv(i,j))
#   ifdef MPI
          endif
#   endif
        enddo
#  endif
# endif /* RIVER || SEAGRASS || REGIONAL */

        do is=1,Nsrc
          Qbar0(is)=Qbar(is)
        enddo

      endif   ! iic.eq.ntstart

!
!  If Qbar is read in nc file, update Qbar0 at every time step
!
# ifdef PSOURCE_NCFILE
      do is=1,Nsrc
        Qbar0(is)=Qbar(is)
      enddo
# endif
!
! Set-up river ramp to avoid strong initial shocks
!
# ifdef RIVER_RAMP
      ramp=TANH(dt*sec2day*float(iic-ntstart)) ! 1-day ramp
# else
      ramp=1.0                                 ! no ramp (by default)
# endif
      do is=1,Nsrc
        Qbar(is)=ramp*Qbar0(is)
      enddo
!
! Set-up vertically integrated mass transport [m3/s] of point
! sources (these may be time-dependent; positive in the positive U-
! or V-direction and vice-versa) and vertically distribute them
! according to mass transport profile chosen above.
!
      do is=1,Nsrc
        do k=1,N
          Qsrc(is,k)=Qbar(is)*Qshape(is,k)
        enddo
      enddo
!
!  Set-up tracer (tracer units) point Sources/Sinks.
!
# ifdef SOLVE3D
#  if defined RIVER || defined SEAGRASS
      do k=1,N
        do is=1,Nsrc
#   ifdef TEMPERATURE
          Tsrc(is,k,itemp)=Tsrc0(is,itemp)
!          Tsrc(is,k,itemp)=4.+10.*exp(z_r(Isrc(is),Jsrc(is),k)/50.)
#   endif
#   ifdef SALINITY
          Tsrc(is,k,isalt)=Tsrc0(is,isalt)
#   endif
#   ifdef PASSIVE_TRACER
          Tsrc(is,k,itpas)=Tsrc0(is,itpas)
#   endif
#   ifdef SEDIMENT
          do itrc=1,NST
            ised=itrc_sed+itrc-1
            Tsrc(is,k,ised)=Tsrc0(is,ised)
          enddo
#   endif
#   ifdef SUBSTANCE
          do itrc=itsubs1,itsubs2
            Tsrc(is,k,itrc)=Tsrc0(is,itrc)
          enddo
#   endif
        enddo
      enddo
#  elif defined REGIONAL
      do k=1,N
         do is=1,Nsrc
#   ifdef PSOURCE_NCFILE_TS
#    ifdef TEMPERATURE
           if (.not.got_tsrc(itemp)) then
              Tsrc0(is,itemp)=5.
           endif
#    endif
#    ifdef SALINITY
           if (.not.got_tsrc(isalt)) then
              Tsrc0(is,isalt)=2.
           endif
#    endif
#   endif

#   ifdef TRACERS
           do itrc=1,NT
             Tsrc(is,k,itrc)=Tsrc0(is,itrc)
           enddo
#   endif
#   ifdef PASSIVE_TRACER
#    ifdef PSOURCE_NCFILE_TS
           if (.not.got_tsrc(itpas)) then
              Tsrc0(is,itpas)=1.
           endif
#    endif
           Tsrc(is,k,itpas)=Tsrc0(is,itpas)
#   endif
#   ifdef SEDIMENT
          do itrc=1,NST
            ised=itrc_sed+itrc-1
            Tsrc(is,k,ised)=Tsrc0(is,ised)
          enddo
#   endif
         enddo
       enddo
#  else
      do k=1,N
        do is=1,Nsrc
          do itrc=1,NT
            Tsrc(is,k,itrc)=Tsrc0(is,itrc)
          enddo
        enddo
      enddo
#  endif /* RIVER || SEAGRASS || REGIONAL */
!
!  Set-up biology tracer (tracer units) point Sources/Sinks.
!
#  ifdef BIOLOGY
      do k=1,N
        do is=1,Nsrc
!     Fill with analytical value  if not find in netcdf file
#   ifdef TEMPERATURE
          temp=Tsrc0(is,itemp)
#   else
          temp=25.
#   endif
          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
          xno3=1.67+0.5873*SiO4+0.0144*SiO4**2
     &         +0.0003099*SiO4**3
#   ifdef PSOURCE_NCFILE_TS
          if(.not.got_tsrc(iNO3_)) then
            Tsrc0(is,iNO3_)=xno3
          endif
#   endif
          Tsrc(is,k,iNO3_)=Tsrc0(is,iNO3_)

#   ifdef PISCES
#    ifdef PSOURCE_NCFILE_TS
          if (.not.got_tsrc(iPOC_)) then
            Tsrc0(is,iPOC_)=0.01
          endif
          if (.not.got_tsrc(iPHY_)) then
            Tsrc0(is,iPHY_)=0.01
          endif
          if (.not.got_tsrc(iZOO_)) then
            Tsrc0(is,iZOO_)=0.01
          endif
#     if ! defined key_pisces_npzd
          if (.not.got_tsrc(iCAL_)) then
            Tsrc0(is,iCAL_)=0.01
          endif
          if (.not.got_tsrc(iDIA_)) then
            Tsrc0(is,iDIA_)=0.01
          endif
          if (.not.got_tsrc(iDSI_)) then
            Tsrc0(is,iDSI_)=1.5e-3
          endif
          if (.not.got_tsrc(iBFE_)) then
            Tsrc0(is,iBFE_)=1.e-2*5.e-6
          endif
          if (.not.got_tsrc(iGOC_)) then
            Tsrc0(is,iGOC_)=0.01
          endif
          if (.not.got_tsrc(iSFE_)) then
            Tsrc0(is,iSFE_)=1.e-2*5.e-6
          endif
          if (.not.got_tsrc(iDFE_)) then
            Tsrc0(is,iDFE_)=1.e-2*5.e-6
          endif
          if (.not.got_tsrc(iGSI_)) then
            Tsrc0(is,iGSI_)=1.e-2*0.15
          endif
          if (.not.got_tsrc(iNFE_)) then
            Tsrc0(is,iNFE_)=1.e-2*5e-6
          endif
          if (.not.got_tsrc(iNCH_)) then
            Tsrc0(is,iNCH_)=1.e-2*12./55.
          endif
          if (.not.got_tsrc(iDCH_)) then
            Tsrc0(is,iDCH_)=1.e-2*12./55.
          endif
          if (.not.got_tsrc(iNH4_)) then
            Tsrc0(is,iNH4_)=1.e-2
          endif
#     endif
#     ifdef key_ligand
          if (.not.got_tsrc(iLGW_)) then
            Tsrc0(is,iLGW_)=1.e-3
          endif
#     endif
#     ifdef key_pisces_quota
!          if (.not.got_tsrc(iDON_)) then
!            Tsrc0(is,iDON_)=5.
!          endif
!          if (.not.got_tsrc(iDOP_)) then
!            Tsrc0(is,iDOP_)=5.
!          endif
          if (.not.got_tsrc(iPON_)) then
            Tsrc0(is,iPON_)=1.e-2
          endif
          if (.not.got_tsrc(iPOP_)) then
            Tsrc0(is,iPOP_)=1.e-2
          endif
          if (.not.got_tsrc(iNPH_)) then
            Tsrc0(is,iNPH_)=1.e-2
          endif
          if (.not.got_tsrc(iPPH_)) then
            Tsrc0(is,iPPH_)=1.e-2
          endif
          if (.not.got_tsrc(iNDI_)) then
            Tsrc0(is,iNDI_)=1.e-2
          endif
          if (.not.got_tsrc(iPDI_)) then
            Tsrc0(is,iPDI_)=1.e-2
          endif
          if (.not.got_tsrc(iPIC_)) then
            Tsrc0(is,iPIC_)=1.e-2
          endif
          if (.not.got_tsrc(iNPI_)) then
            Tsrc0(is,iNPI_)=1.e-2
          endif
          if (.not.got_tsrc(iPPI_)) then
            Tsrc0(is,iPPI_)=1.e-2
          endif
          if (.not.got_tsrc(iPFE_)) then
            Tsrc0(is,iPFE_)=1.e-2*5.e-6
          endif
          if (.not.got_tsrc(iPCH_)) then
            Tsrc0(is,iPCH_)=1.e-2*12./55.
          endif
          if (.not.got_tsrc(iGON_)) then
            Tsrc0(is,iGON_)=1.e-2
          endif
          if (.not.got_tsrc(iGOP_)) then
            Tsrc0(is,iGOP_)=1.e-2
          endif
#     endif
#    endif   /* PSOURCE_NCFILE_TS */
          Tsrc(is,k,iPOC_)=Tsrc0(is,iPOC_)
          Tsrc(is,k,iPHY_)=Tsrc0(is,iPHY_)
          Tsrc(is,k,iZOO_)=Tsrc0(is,iZOO_)
#    if ! defined key_pisces_npzd
          Tsrc(is,k,iCAL_)=Tsrc0(is,iCAL_)
          Tsrc(is,k,iDIA_)=Tsrc0(is,iDIA_)
          Tsrc(is,k,iDSI_)=Tsrc0(is,iDSI_)
          Tsrc(is,k,iBFE_)=Tsrc0(is,iBFE_)
          Tsrc(is,k,iGOC_)=Tsrc0(is,iGOC_)
          Tsrc(is,k,iSFE_)=Tsrc0(is,iSFE_)
          Tsrc(is,k,iDFE_)=Tsrc0(is,iDFE_)
          Tsrc(is,k,iGSI_)=Tsrc0(is,iGSI_)
          Tsrc(is,k,iNFE_)=Tsrc0(is,iNFE_)
          Tsrc(is,k,iNCH_)=Tsrc0(is,iNCH_)
          Tsrc(is,k,iDCH_)=Tsrc0(is,iDCH_)
          Tsrc(is,k,iNH4_)=Tsrc0(is,iNH4_)
#    endif
#    ifdef key_ligand
          Tsrc(is,k,iLGW_)=Tsrc0(is,iLGW_)
#    endif
#    ifdef key_pisces_quota
!          Tsrc(is,k,iDON_)=Tsrc0(is,iDON_)
!          Tsrc(is,k,iDOP_)=Tsrc0(is,iDOP_)
          Tsrc(is,k,iPON_)=Tsrc0(is,iPON_)
          Tsrc(is,k,iPOP_)=Tsrc0(is,iPOP_)
          Tsrc(is,k,iNPH_)=Tsrc0(is,iNPH_)
          Tsrc(is,k,iPPH_)=Tsrc0(is,iPPH_)
          Tsrc(is,k,iNDI_)=Tsrc0(is,iNDI_)
          Tsrc(is,k,iPDI_)=Tsrc0(is,iPDI_)
          Tsrc(is,k,iPIC_)=Tsrc0(is,iPIC_)
          Tsrc(is,k,iNPI_)=Tsrc0(is,iNPI_)
          Tsrc(is,k,iPPI_)=Tsrc0(is,iPPI_)
          Tsrc(is,k,iPFE_)=Tsrc0(is,iPFE_)
          Tsrc(is,k,iPCH_)=Tsrc0(is,iPCH_)
          Tsrc(is,k,iGON_)=Tsrc0(is,iGON_)
          Tsrc(is,k,iGOP_)=Tsrc0(is,iGOP_)
#    endif
#   elif defined BIO_NChlPZD
#    ifdef PSOURCE_NCFILE_TS
          if (.not.got_tsrc(iChla)) then
            Tsrc0(is,iChla)=0.08
          endif
          if (.not.got_tsrc(iPhy1)) then
            Tsrc0(is,iPhy1)=0.1
          endif
          if (.not.got_tsrc(iZoo1)) then
            Tsrc0(is,iZoo1)=0.06
          endif
          if (.not.got_tsrc(iDet1)) then
            Tsrc0(is,iDet1)=0.02
          endif
#    endif   /* PSOURCE_NCFILE_TS    */
          Tsrc(is,k,iChla)=Tsrc0(is,iChla)
          Tsrc(is,k,iPhy1)=Tsrc0(is,iPhy1)
          Tsrc(is,k,iZoo1)=Tsrc0(is,iZoo1)
          Tsrc(is,k,iDet1)=Tsrc0(is,iDet1)

#   elif defined BIO_N2ChlPZD2
#    ifdef PSOURCE_NCFILE_TS
          if (.not.got_tsrc(iNH4_)) then
            Tsrc0(is,iNH4_)=0.1
          endif
          if (.not.got_tsrc(iChla)) then
            Tsrc0(is,iChla)=0.08
          endif
          if (.not.got_tsrc(iPhy1)) then
            Tsrc0(is,iPhy1)=0.06

          endif
          if (.not.got_tsrc(iZoo1)) then
            Tsrc0(is,iZoo1)=0.04

          endif
          if (.not.got_tsrc(iDet1)) then
            Tsrc0(is,iDet1)=0.02

          endif
          if (.not.got_tsrc(iDet2)) then
            Tsrc0(is,iDet2)=0.02
          endif
#    endif
          Tsrc(is,k,iNH4_)= Tsrc0(is,iNH4_)
          Tsrc(is,k,iChla)=Tsrc0(is,iChla)
          Tsrc(is,k,iPhy1)=Tsrc0(is,iPhy1)
          Tsrc(is,k,iZoo1)=Tsrc0(is,iZoo1)
          Tsrc(is,k,iDet1)=Tsrc0(is,iDet1)
          Tsrc(is,k,iDet2)=Tsrc0(is,iDet2)
!----------------------------------------------------------------------
#   elif defined BIO_BioEBUS
          zsrc=z_r(Isrc(is),Jsrc(is),k)
#    ifdef PSOURCE_NCFILE_TS
          if (.not.got_tsrc(iNO2_)) then
            Tsrc0(is,iNO2_)=  0.05
          endif
          if (.not.got_tsrc(iNH4_)) then
            Tsrc0(is,iNH4_)= 0.1
          endif
          if (.not.got_tsrc(iPhy1)) then
            Tsrc0(is,iPhy1)= 0.04
          endif
          if (.not.got_tsrc(iPhy2)) then
            Tsrc0(is,iPhy2)= 0.06
          endif
          if (.not.got_tsrc(iZoo1)) then
            Tsrc0(is,iZoo1)= 0.04
          endif
          if (.not.got_tsrc(iZoo2)) then
            Tsrc0(is,iZoo2)= 0.
          endif
          if (.not.got_tsrc(iDet1)) then
            Tsrc0(is,iDet1)=0.02
          endif
          if (.not.got_tsrc(iDet2)) then
            Tsrc0(is,iDet2)=0.02
          endif
          if (.not.got_tsrc(iDON)) then
            Tsrc0(is,iDON)= 0.5
          endif
#     ifdef NITROUS_OXIDE
          if (.not.got_tsrc(iN2O)) then
            Tsrc0(is,iN2O)= 0.002 ! not used ...
          endif
#     endif
#    endif  /* PSOURCE_NCFILE_TS  */
          Tsrc(is,k,iNO2_)= Tsrc0(is,iNO2_)*exp(zsrc/100.)
          Tsrc(is,k,iNH4_)= Tsrc0(is,iNH4_)*exp(zsrc/100.)
          Tsrc(is,k,iPhy1)=Tsrc0(is,iPhy1)*exp(zsrc/50.)
          Tsrc(is,k,iPhy2)=Tsrc0(is,iPhy2)*exp(zsrc/50.)
          Tsrc(is,k,iZoo1)=Tsrc0(is,iZoo1)*exp(zsrc/100.)
          Tsrc(is,k,iZoo2)=Tsrc0(is,iZoo2)*exp(zsrc/100.)
          Tsrc(is,k,iDet1)=Tsrc0(is,iDet1)
          Tsrc(is,k,iDet2)=Tsrc0(is,iDet2)
          Tsrc(is,k,iDON)=Tsrc0(is,iDON)*exp(zsrc/100.)
          Tsrc(is,k,iN2O)=-(0.008*exp(zsrc/100.)-0.01)
          if (Tsrc(is,k,iN2O).lt.0.) Tsrc(is,k,iN2O)=0.
#   endif /* PISCES or BIO_NChlPZD or BIO_N2ChlPZD2 or BIO_BioEBUS*/
!----------------------------------------------------------------------
        enddo
      enddo
#  endif /* BIOLOGY */
# endif /* SOLVE3D */
!$acc end kernels 
      return
      end

      subroutine get_qshape_exp(Qshapeij, hij, hinvij)
        implicit none
# include "param.h"
# include "scalars.h"


        real, dimension(N), intent(inout) :: Qshapeij
        real, intent(in) :: hij, hinvij

        real :: cff_w, cff_r, cff1_w, cff1_r, Hs, cff, cff1, cff2
        real, dimension(0:N) :: z_w0
        real, dimension(1:N) :: z_r0
        integer :: k

        !! Compute z_w0 and z_r0 for a zeta=0. case to initialize Qshape
        !! without time dependancy
        z_w0(0)=-hij
        do k=1,N,+1
#  ifdef NEW_S_COORD
          cff_w =hc*sc_w(k)
          cff_r =hc*sc_r(k)
          cff1_w=Cs_w(k)
          cff1_r=Cs_r(k)
          z_w0(k)=(cff_w+cff1_w*abs(hij)) *hij*hinvij
          z_r0(k)=(cff_r+cff1_r*abs(hij)) *hij*hinvij
#  else
          cff_w =hc*(sc_w(k)-Cs_w(k))
          cff_r =hc*(sc_r(k)-Cs_r(k))
          cff1_w=Cs_w(k)
          cff1_r=Cs_r(k)
          z_w0(k)=cff_w+cff1_w*hij
          z_r0(k)=cff_r+cff1_r*hij
#  endif
        enddo

        Hs=hij                    ! Exponential vertical distribution
        cff=5.                    ! Hs/z0  (z0 surface layer depth)
        cff1=cff/(1-exp(-cff))
        cff2=0.
        do k=1,N
          Qshapeij(k)=cff1*exp(z_r0(k)*cff/Hs)*(z_w0(k)-z_w0(k-1))/Hs
          cff2=cff2+Qshapeij(k)
        enddo
        do k=1,N
          Qshapeij(k)=Qshapeij(k)/cff2
        enddo

      return
      end

#endif /* PSOURCE && ANA_PSOURCE */
!
!======================================================================
!                   subroutine ana_bry
!======================================================================
!
#ifdef ANA_BRY
      subroutine ana_bry (tile)
      implicit none
# include "param.h"
      integer tile
# ifdef ALLOW_SINGLE_BLOCK_MODE
C$    integer  trd, omp_get_thread_num
# endif
# include "compute_tile_bounds.h"
      call ana_bry_tile   (Istr,Iend,Jstr,Jend)
      return
      end
!
      subroutine ana_bry_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  Set analytical boundary forcing
!----------------------------------------------------------------------
!
# ifdef SUBSTANCE
      USE comsubstance, ONLY: cobc_wat
# endif
      implicit none
      integer Istr,Iend,Jstr,Jend, i,j,k, itrc
      real    eps
      parameter (eps=1.e-10)
# if defined WAVE_MAKER || defined WAVE_MAKER_INTERNAL
      integer iw,jw, jwtt0,jwtt
      real    Du,Dv,Zu,Zv,Zr,Dr,xu,yu,xv,yv,x0,y0,time0,
     &        h0,khd,kh,wa,wk,wp,wf,wd,wds,ramp,
     &        sigma,theta,gamma,Hs,cff_spread, sumspec,
     &        cff,cff0,cff1,cff2,cff3,cff4,cff5,
     &        fmin,fmax,df,dmin,dmax,dd, 
     &        mindiff,wkdom,wky,
     &        wa1,wk1,wp1,wf1,wa2,wk2,wp2,wf2
      real    K1, K2, K3, K4, K5, K6
      parameter (K1=0.6666666666, K2=0.3555555555,
     &           K3=0.1608465608, K4=0.0632098765,
     &           K5=0.0217540484, K6=0.0065407983)
      INTEGER nseed
      INTEGER, ALLOCATABLE :: seed(:)
# endif

# include "param.h"
# include "grid.h"
# include "scalars.h"
# include "boundary.h"
# include "ocean3d.h"
# include "ocean2d.h"
# include "coupling.h"
# include "forces.h"
# include "mixing.h"
# ifdef WAVE_MAKER
      real wd_bry_tmp(Nfrq)
# endif
!
# include "compute_auxiliary_bounds.h"

# ifdef OBC_WEST
!
!================== WESTERN BOUNDARY ======================
!
      if (WESTERN_EDGE) then

#  if defined WAVE_MAKER && !defined WAVE_MAKER_EAST
!
!  --- Surface wave maker on western boundary ---
!
#   include "wave_maker.h"

#  elif defined WAVE_MAKER_INTERNAL
!
!  --- Internal wave maker on western boundary ---
!
#   include "wave_maker_internal.h"

#  else

        if (FIRST_RST_TIME_STEP) then
#   ifdef Z_FRC_BRY
          do j=JstrR,JendR
            zetabry_west(j)=zeta(0,j,1)
#    ifdef OBC_PATM
            zetabry_west(j)=zetabry_west(j)
     &                     -(patm2d(Istr,j)-paref)/(g*rho0)
#    endif
#    if defined TIDAL_FLAT || defined ESTUARY
            zetabry_west(j)=2.*sin(2.*pi*time/(12.0*3600.0))
#     if defined ESTUARY_TIDE
            ! adding another wave to modulate tide
            zetabry_west(j) = zetabry_west(j) +
     &          2.*sin(2.*pi*time/((24.0+50./60.)/2.*3600.0))
#     endif
#    endif
          enddo
#   endif
#   ifdef M2_FRC_BRY
          do j=JstrR,JendR
#    if defined TIDAL_FLAT || defined ESTUARY
            ubarbry_west(j)=0.
#    else
            ubarbry_west(j)=ubar(1,j,1)
#    endif
            vbarbry_west(j)=vbar(0,j,1)
          enddo
#   endif
#   if defined SOLVE3D && \
     (defined M3_FRC_BRY || defined W_FRC_BRY || defined T_FRC_BRY)
          do k=1,N
            do j=JstrR,JendR
#    ifdef M3_FRC_BRY
#     if defined TIDAL_FLAT || defined ESTUARY
              ubry_west(j,k)=0.
#     else
              ubry_west(j,k)=u(1,j,k,1)
#     endif
              vbry_west(j,k)=v(0,j,k,1)
#    endif
#    ifdef W_FRC_BRY
              wbry_west(j,k)=wz(0,j,k,1)
#    endif
#    if defined T_FRC_BRY && defined TRACERS
              do itrc=1,NT
                tbry_west(j,k,itrc)=t(0,j,k,1,itrc)
              enddo
#    endif
#    ifdef SUBSTANCE
#     if defined TIDAL_FLAT || defined ESTUARY
            tbry_west(j,k,1)=10.
            tbry_west(j,k,2)=35.5
#     else
            do itrc=1,itsubs1-1
              tbry_west(j,k,itrc)=0.
            enddo
#     endif
            do itrc=itsubs1,NT
              tbry_west(j,k,itrc)=cobc_wat(itrc)
            enddo
#    else
#     if defined TIDAL_FLAT
            tbry_west(j,k,1)=10.
            tbry_west(j,k,2)=35.5
            tbry_west(j,k,3)=0.   ! tracer Coarse Sand
            tbry_west(j,k,4)=0.   ! tracer Fine Sand
            tbry_west(j,k,5)=0.02 ! tracer Mud
#     endif
#     if defined ESTUARY
            tbry_west(j,k,1)=10.
            tbry_west(j,k,2)=35.5
            tbry_west(j,k,3)=0.0   ! tracer Fine Sand
            tbry_west(j,k,4)=0.0   ! tracer Mud
#     endif
#    endif
            enddo
          enddo
#   endif    /* SOLVE3D && (M3_FRC_BRY || W_FRC_BRY || T_FRC_BRY)*/
        else
#   ifdef Z_FRC_BRY
          do j=JstrR,JendR
            zetabry_west(j)=0.0
#    ifdef OBC_PATM
            zetabry_west(j)=-(patm2d(Istr,j)-paref)/(g*rho0)
#    endif
#    if defined TIDAL_FLAT || defined ESTUARY
            zetabry_west(j)=2.*sin(2.*pi*time/(12.0*3600.0))
#     if defined ESTUARY_TIDE
            ! adding another wave to modulate tide
            zetabry_west(j) = zetabry_west(j) +
     &          2.*sin(2.*pi*time/((24.0+50./60.)/2.*3600.0))
#     endif
#    endif
          enddo
#   endif
#   ifdef M2_FRC_BRY
          do j=JstrR,JendR
            ubarbry_west(j)=0.0
            vbarbry_west(j)=0.0
          enddo
#   endif
#   if defined SOLVE3D && \
     (defined M3_FRC_BRY || defined W_FRC_BRY || defined T_FRC_BRY)
          do k=1,N
            do j=JstrR,JendR
#    ifdef M3_FRC_BRY
              ubry_west(j,k)=0.0
              vbry_west(j,k)=0.0
#    endif
#    ifdef W_FRC_BRY
              wbry_west(j,k)=wz(0,j,k,1)
#    endif
#    ifdef T_FRC_BRY
#     ifdef SUBSTANCE
#      if defined TIDAL_FLAT || defined ESTUARY
             tbry_west(j,k,1)=10.
             tbry_west(j,k,2)=35.5
#      else
             do itrc=1,itsubs1-1
               tbry_west(j,k,itrc)=0.
             enddo
#      endif
             do itrc=itsubs1,NT
               tbry_west(j,k,itrc)=cobc_wat(itrc)
             enddo
#     else
#      if defined TIDAL_FLAT
             tbry_west(j,k,1)=10.
             tbry_west(j,k,2)=35.
             tbry_west(j,k,3)=0. ! tracer Coarse Sand
             tbry_west(j,k,4)=0. ! tracer Fine Sand
             tbry_west(j,k,5)=0.02 ! tracer Mud
#      elif defined ESTUARY
             tbry_west(j,k,1)=10.
             tbry_west(j,k,2)=35.5
             tbry_west(j,k,3)=0.0   ! tracer Fine Sand
             tbry_west(j,k,4)=0.0   ! tracer Mud
#      else
             do itrc=1,NT
               tbry_west(j,k,itrc)=t(0,j,k,1,itrc)
             enddo
#      endif
#     endif
#    endif
            enddo
          enddo
#   endif    /* SOLVE3D && (M3_FRC_BRY || W_FRC_BRY || T_FRC_BRY)*/

        endif
#  endif /* WAVE_MAKER */
      endif
# endif /* OBC_WEST */
!
!================== EASTERN BOUNDARY ======================
!
# ifdef OBC_EAST
      if (EASTERN_EDGE) then

#  if defined WAVE_MAKER && defined WAVE_MAKER_EAST
!
!  --- Surface wave maker on eastern boundary ---
!
#   include "wave_maker.h"

!#  elif defined WAVE_MAKER_INTERNAL_EAST
!
!  --- Internal wave maker on eastern boundary ---
!
!#   include "wave_maker_internal.h"

#  else

        if (FIRST_RST_TIME_STEP) then
#   ifdef Z_FRC_BRY
          do j=JstrR,JendR
            zetabry_east(j)=zeta(Iend+1,j,1)
#    ifdef OBC_PATM
            zetabry_east(j)=zetabry_east(j)
     &                    -(patm2d(IendR,j)-paref)/(g*rho0)
#    endif
          enddo
#   endif
#   ifdef M2_FRC_BRY
          do j=JstrR,JendR
            ubarbry_east(j)=ubar(Iend+1,j,1)
            vbarbry_east(j)=vbar(Iend+1,j,1)
          enddo
#   endif
#   if defined SOLVE3D && \
     (defined M3_FRC_BRY || defined W_FRC_BRY || defined T_FRC_BRY)
          do k=1,N
            do j=JstrR,JendR
#    ifdef M3_FRC_BRY
              ubry_east(j,k)=u(Iend+1,j,k,1)
              vbry_east(j,k)=v(Iend+1,j,k,1)
#    endif
#    ifdef W_FRC_BRY
              wbry_east(j,k)=wz(Iend+1,j,k,1)
#    endif
#    if defined T_FRC_BRY && defined TRACERS
#     ifdef SUBSTANCE
              do itrc=1,itsubs1-1
                tbry_east(j,k,itrc)=0.
              enddo
              do itrc=itsubs1,NT
                tbry_east(j,k,itrc)=cobc_wat(itrc)
              enddo
#     else
              do itrc=1,NT
                tbry_east(j,k,itrc)=t(Iend+1,j,k,1,itrc)
              enddo
#     endif
#    endif
            enddo
          enddo
#   endif    /* SOLVE3D && (M3_FRC_BRY || W_FRC_BRY || T_FRC_BRY)*/
        else
#   ifdef Z_FRC_BRY
          do j=JstrR,JendR
            zetabry_east(j)=0.0
#    ifdef OBC_PATM
            zetabry_east(j)=-(patm2d(IendR,j)-paref)/(g*rho0)
#    endif
          enddo
#   endif
#   ifdef M2_FRC_BRY
          do j=JstrR,JendR
            ubarbry_east(j)=0.
            vbarbry_east(j)=0.
          enddo
#   endif
#   if defined SOLVE3D && \
     (defined M3_FRC_BRY || defined W_FRC_BRY || defined T_FRC_BRY)
          do k=1,N
            do j=JstrR,JendR
#    ifdef M3_FRC_BRY
              ubry_east(j,k)=0.0
              vbry_east(j,k)=0.0
#    endif
#    ifdef W_FRC_BRY
              wbry_east(j,k)=wz(Iend+1,j,k,1)
#    endif
#    ifdef T_FRC_BRY
#     ifdef SUBSTANCE
              do itrc=1,itsubs1-1
                tbry_east(j,k,itrc)=0.
              enddo
              do itrc=itsubs1,NT
                tbry_east(j,k,itrc)=cobc_wat(itrc)
              enddo
#     else
              do itrc=1,NT
                tbry_east(j,k,itrc)=t(Iend+1,j,k,1,itrc)
              enddo
#     endif
#    endif
            enddo
          enddo
#   endif    /* SOLVE3D && (M3_FRC_BRY || W_FRC_BRY || T_FRC_BRY)*/
        endif
#  endif /* WAVE_MAKER */
      endif
# endif /* OBC_EAST */
!
!================== SOUTHERN BOUNDARY ======================
!
# ifdef OBC_SOUTH
      if (SOUTHERN_EDGE) then
        if (FIRST_RST_TIME_STEP) then
#  ifdef Z_FRC_BRY
          do i=IstrR,IendR
            zetabry_south(i)=zeta(i,0,1)
#   ifdef OBC_PATM
            zetabry_south(i)=zetabry_south(i)
     &                    -(patm2d(i,JstrR)-paref)/(g*rho0)
#   endif
          enddo
#  endif
#  ifdef M2_FRC_BRY
          do i=IstrR,IendR
            ubarbry_south(i)=ubar(i,0,1)
            vbarbry_south(i)=vbar(i,1,1)
          enddo
#  endif
#  if defined SOLVE3D && \
     (defined M3_FRC_BRY || defined W_FRC_BRY || defined T_FRC_BRY)
          do k=1,N
            do i=IstrR,IendR
#   ifdef M3_FRC_BRY
              ubry_south(i,k)=u(i,0,k,1)
              vbry_south(i,k)=v(i,1,k,1)
#   endif
#   ifdef W_FRC_BRY
              wbry_south(i,k)=wz(i,0,k,1)
#   endif
#   if defined T_FRC_BRY && defined TRACERS
#    ifdef SUBSTANCE
              do itrc=1,itsubs1-1
                tbry_south(i,k,itrc)=0.
              enddo
              do itrc=itsubs1,NT
                tbry_south(i,k,itrc)=cobc_wat(itrc)
              enddo
#    else
              do itrc=1,NT
                tbry_south(i,k,itrc)= t(i,0,k,1,itrc)  ! =0.??
              enddo
#    endif
#   endif
            enddo
          enddo
#  endif    /* SOLVE3D && (M3_FRC_BRY || W_FRC_BRY || T_FRC_BRY)*/
        else
#  ifdef Z_FRC_BRY
          do i=IstrR,IendR
            zetabry_south(i)=0.0
#   ifdef OBC_PATM
            zetabry_south(i)=-(patm2d(i,JstrR)-paref)/(g*rho0)
#   endif
          enddo
#  endif
#  ifdef M2_FRC_BRY
          do i=IstrR,IendR
            ubarbry_south(i)=0.0
            vbarbry_south(i)=0.0
          enddo
#  endif
        endif
      endif
# endif /* OBC_SOUTH */
!
!================== NORTHERN BOUNDARY ======================
!
# ifdef OBC_NORTH
      if (NORTHERN_EDGE) then
        if (FIRST_RST_TIME_STEP) then
#  ifdef Z_FRC_BRY
          do i=IstrR,IendR
            zetabry_north(i)=zeta(i,Jend+1,1)
#   ifdef OBC_PATM
            zetabry_north(i)=zetabry_north(i)
     &                     -(patm2d(i,JendR)-paref)/(g*rho0)
#   endif
          enddo
#  endif
#  ifdef M2_FRC_BRY
          do i=IstrR,IendR
            ubarbry_north(i)=ubar(i,Jend+1,1)
            vbarbry_north(i)=vbar(i,Jend+1,1)
          enddo
#  endif
#  if defined SOLVE3D && \
     (defined M3_FRC_BRY || defined W_FRC_BRY || defined T_FRC_BRY)
          do k=1,N
            do i=IstrR,IendR
#   ifdef M3_FRC_BRY
              ubry_north(i,k)=u(i,Jend+1,k,1)
              vbry_north(i,k)=v(i,Jend+1,k,1)
#   endif
#   ifdef W_FRC_BRY
              wbry_north(i,k)=wz(i,Jend+1,k,1)
#   endif
#   if defined T_FRC_BRY && defined TRACERS
#    ifdef SUBSTANCE
              do itrc=1,itsubs1-1
                tbry_north(i,k,itrc)=0.
              enddo
              do itrc=itsubs1,NT
                tbry_north(i,k,itrc)=cobc_wat(itrc)
              enddo
#    else
              do itrc=1,NT
                tbry_north(i,k,itrc)=t(i,Jend+1,k,1,itrc) ! =0.??
              enddo
#    endif
#   endif
            enddo
          enddo
#  endif    /* SOLVE3D && (M3_FRC_BRY || W_FRC_BRY || T_FRC_BRY)*/
        else
#  ifdef Z_FRC_BRY
          do i=IstrR,IendR
            zetabry_north(i)=0.0
#   ifdef OBC_PATM
            zetabry_north(i)=-(patm2d(i,JendR)-paref)/(g*rho0)
#   endif
          enddo
#  endif
#  ifdef M2_FRC_BRY
          do i=IstrR,IendR
            ubarbry_north(i)=0.0
            vbarbry_north(i)=0.0
          enddo
#  endif
        endif
      endif
# endif /* OBC_NORTH */
!
      return
      end
#endif /* defined ANA_BRY */
!
!======================================================================
!                   subroutine ana_nbq_bry
!======================================================================
!
#ifdef NBQ_FRC_BRY
      subroutine ana_nbq_bry (tile)
      implicit none
# include "param.h"
      integer tile
# ifdef ALLOW_SINGLE_BLOCK_MODE
C$    integer  trd, omp_get_thread_num
# endif
# include "compute_tile_bounds.h"
      call ana_nbq_bry_tile (Istr,Iend,Jstr,Jend)
      return
      end
!
      subroutine ana_nbq_bry_tile (Istr,Iend,Jstr,Jend)
!
!----------------------------------------------------------------------
!  Set analytical boundary forcing
!----------------------------------------------------------------------
!
      implicit none
      integer Istr,Iend,Jstr,Jend, i,j,k, itrc
# include "param.h"
# include "grid.h"
# include "scalars.h"
# include "boundary.h"
# include "climat.h"
# include "ocean3d.h"
!
# include "compute_auxiliary_bounds.h"
!
# ifdef EW_PERIODIC
#  define IU_RANGE Istr,Iend
#  define IV_RANGE Istr,Iend
# else
#  define IU_RANGE Istr,IendR
#  define IV_RANGE IstrR,IendR
# endif

# ifdef NS_PERIODIC
#  define JU_RANGE Jstr,Jend
#  define JV_RANGE Jstr,Jend
# else
#  define JU_RANGE JstrR,JendR
#  define JV_RANGE Jstr,JendR
# endif

# ifdef OBC_WEST
      if (WESTERN_EDGE) then
        do k=1,N
          do j=JU_RANGE
            unbqbry_west(j,k)=0.5*(Hz(0,j,k)+Hz(1,j,k))
     &                                   *u(1,j,k,nrhs)
          enddo
          do j=JV_RANGE
            vnbqbry_west(j,k)=0.5*(Hz(0,j,k)+Hz(0,j-1,k))
     &                                   *v(0,j,k,nrhs)
          enddo
#  ifdef NBQ
          do j=IstrR,IendR
            wnbqbry_west(j,k)=0.
            rnbqbry_west(j,k)=0.
          enddo
#  endif
        enddo
      endif
# endif /* OBC_WEST */
!
# ifdef OBC_EAST
      if (EASTERN_EDGE) then
        do k=1,N
          do j=JU_RANGE
            unbqbry_east(j,k)=0.5*(Hz(IendR,j,k)+Hz(IendR-1,j,k))
     &                                         *u(IendR,j,k,nrhs)
#  ifdef NBQ
            wnbqbry_east(j,k)=0.
            rnbqbry_east(j,k)=0.
#  endif
          enddo
          do j=JV_RANGE
            vnbqbry_east(j,k)=0.5*(Hz(IendR,j,k)+Hz(IendR,j-1,k))
     &                                         *v(IendR,j,k,nrhs)
          enddo
        enddo
      endif
# endif /* OBC_EAST */
!
# ifdef OBC_SOUTH
      if (SOUTHERN_EDGE) then
        do k=1,N
          do i=IV_RANGE
            vnbqbry_south(i,k)=0.5*(Hz(i,0,k)+Hz(i,1,k))
     &                                    *v(i,1,k,nrhs)
#  ifdef NBQ
            wnbqbry_south(i,k)=0.
            rnbqbry_south(i,k)=0.
#  endif
          enddo
          do i=IU_RANGE
            unbqbry_south(i,k)=0.5*(Hz(i,0,k)+Hz(i-1,0,k))
     &                                    *u(i,0,k,nrhs)
          enddo
        enddo
      endif
# endif /* OBC_SOUTH */
!
# ifdef OBC_NORTH
      if (NORTHERN_EDGE) then
        do k=1,N
          do i=IV_RANGE
            vnbqbry_north(i,k)=0.5*(Hz(i,JendR,k)+Hz(i,JendR-1,k))
     &                                          *v(i,JendR,k,nrhs)
#  ifdef NBQ
            wnbqbry_north(i,k)=0.
            rnbqbry_north(i,k)=0.
#  endif
          enddo
          do i=IU_RANGE
            unbqbry_north(i,k)=0.5*(Hz(i,JendR,k)+Hz(i-1,JendR,k))
     &                                          *u(i,JendR,k,nrhs)
          enddo
        enddo
      endif
# endif /* OBC_NORTH */

      return
      end
#endif /* defined NBQ_FRC_BRY */
!
!======================================================================
!                   subroutine ana_bry_bio
!======================================================================
!
#if defined BIOLOGY && defined T_FRC_BRY
!
!----------------------------------------------------------------------
!  Set analytical boundary forcing for biological tracers
!----------------------------------------------------------------------
!
      subroutine ana_bry_bio (tile)
      implicit none
      integer tile
# include "param.h"
# ifdef ALLOW_SINGLE_BLOCK_MODE
C$    integer  trd, omp_get_thread_num
# endif
# include "compute_tile_bounds.h"
      call ana_bry_bio_tile (Istr,Iend,Jstr,Jend)
      return
      end
!
      subroutine ana_bry_bio_tile (Istr,Iend,Jstr,Jend)
      implicit none
      integer Istr,Iend,Jstr,Jend, i,j,k, itrc
      real    xno3,temp,SiO4,zbry
# include "param.h"
# include "boundary.h"
# include "scalars.h"
# include "ocean3d.h"
!
# include "compute_auxiliary_bounds.h"
!
# ifdef ANA_BRY
      do i=1,itrc
        got_tbry(i)=.false.
      enddo
# endif

# ifdef OBC_WEST
      if (WESTERN_EDGE) then
        do k=1,N
          do j=JstrR,JendR
#  ifdef TEMPERATURE
            temp=tbry_west(j,k,itemp)
#  else
            temp=25.
#  endif
            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
            xno3=1.67+0.5873*SiO4+0.0144*SiO4**2
     &           +0.0003099*SiO4**3
            if (.not.got_tbry(iNO3_)) tbry_west(j,k,iNO3_)=xno3
#  ifdef PISCES
            if (.not.got_tbry(iDIC_)) tbry_west(j,k,iDIC_)=2150.
            if (.not.got_tbry(iTAL_)) tbry_west(j,k,iTAL_)=2350.
            if (.not.got_tbry(iOXY_)) tbry_west(j,k,iOXY_)=200.
            if (.not.got_tbry(iPOC_)) tbry_west(j,k,iPOC_)=0.01
            if (.not.got_tbry(iPHY_)) tbry_west(j,k,iPHY_)=0.01
            if (.not.got_tbry(iZOO_)) tbry_west(j,k,iZOO_)=0.01
            if (.not.got_tbry(iDOC_)) tbry_west(j,k,iDOC_)=5.
            if (.not.got_tbry(iFER_)) tbry_west(j,k,iFER_)=6.e-4
#   if ! defined key_pisces_npzd
            if (.not.got_tbry(iCAL_)) tbry_west(j,k,iCAL_)=0.01
            if (.not.got_tbry(iPO4_)) tbry_west(j,k,iPO4_)=xno3/16.
            if (.not.got_tbry(iSIL_)) tbry_west(j,k,iSIL_)=91.51
            if (.not.got_tbry(iDIA_)) tbry_west(j,k,iDIA_)=0.01
            if (.not.got_tbry(iMES_)) tbry_west(j,k,iMES_)=0.01
            if (.not.got_tbry(iDSI_)) tbry_west(j,k,iDSI_)=1.5e-3
            if (.not.got_tbry(iBFE_)) tbry_west(j,k,iBFE_)=1.e-2*5.e-6
            if (.not.got_tbry(iGOC_)) tbry_west(j,k,iGOC_)=0.01
            if (.not.got_tbry(iSFE_)) tbry_west(j,k,iSFE_)=1.e-2*5.e-6
            if (.not.got_tbry(iDFE_)) tbry_west(j,k,iDFE_)=1.e-2*5.e-6
            if (.not.got_tbry(iGSI_)) tbry_west(j,k,iGSI_)=1.e-2*0.15
            if (.not.got_tbry(iNFE_)) tbry_west(j,k,iNFE_)=1.e-2*5e-6
            if (.not.got_tbry(iNCH_)) tbry_west(j,k,iNCH_)=1.e-2*12./55.
            if (.not.got_tbry(iDCH_)) tbry_west(j,k,iDCH_)=1.e-2*12./55.
            if (.not.got_tbry(iNH4_)) tbry_west(j,k,iNH4_)=1.e-2
#   endif
#   ifdef key_ligand
            if (.not.got_tbry(iLGW_)) tbry_west(j,k,iLGW_)=1.e-3
#   endif
#   ifdef key_pisces_quota
            if (.not.got_tbry(iDON_)) tbry_west(j,k,iDON_)=5.
            if (.not.got_tbry(iDOP_)) tbry_west(j,k,iDOP_)=5.
            if (.not.got_tbry(iPON_)) tbry_west(j,k,iPON_)=1.e-2
            if (.not.got_tbry(iPOP_)) tbry_west(j,k,iPOP_)=1.e-2
            if (.not.got_tbry(iNPH_)) tbry_west(j,k,iNPH_)=1.e-2
            if (.not.got_tbry(iPPH_)) tbry_west(j,k,iPPH_)=1.e-2
            if (.not.got_tbry(iNDI_)) tbry_west(j,k,iNDI_)=1.e-2
            if (.not.got_tbry(iPDI_)) tbry_west(j,k,iPDI_)=1.e-2
            if (.not.got_tbry(iPIC_)) tbry_west(j,k,iPIC_)=1.e-2
            if (.not.got_tbry(iNPI_)) tbry_west(j,k,iNPI_)=1.e-2
            if (.not.got_tbry(iPPI_)) tbry_west(j,k,iPPI_)=1.e-2
            if (.not.got_tbry(iPFE_)) tbry_west(j,k,iPFE_)=1.e-2*5.E-6
            if (.not.got_tbry(iPCH_)) tbry_west(j,k,iPCH_)=1.e-2*12./55.
            if (.not.got_tbry(iGON_)) tbry_west(j,k,iGON_)=1.e-2
            if (.not.got_tbry(iGOP_)) tbry_west(j,k,iGOP_)=1.e-2
#   endif
#  elif defined BIO_NChlPZD
            if (.not.got_tbry(iChla)) tbry_west(j,k,iChla)=0.08
            if (.not.got_tbry(iPhy1)) tbry_west(j,k,iPhy1)=0.1
            if (.not.got_tbry(iZoo1)) tbry_west(j,k,iZoo1)=0.06
            if (.not.got_tbry(iDet1)) tbry_west(j,k,iDet1)=0.02
#  elif defined BIO_N2ChlPZD2
            if (.not.got_tbry(iNH4_)) tbry_west(j,k,iNH4_)=0.1
            if (.not.got_tbry(iChla)) tbry_west(j,k,iChla)=0.08
            if (.not.got_tbry(iPhy1)) tbry_west(j,k,iPhy1)=0.06
            if (.not.got_tbry(iZoo1)) tbry_west(j,k,iZoo1)=0.04
            if (.not.got_tbry(iDet1)) tbry_west(j,k,iDet1)=0.02
            if (.not.got_tbry(iDet2)) tbry_west(j,k,iDet2)=0.02
#  elif defined BIO_BioEBUS
            zbry=z_r(Istr,j,k)
            if (.not.got_tbry(iNO2_)) tbry_west(j,k,iNO2_)=0.05*
     &                               exp(zbry/100.)
            if (.not.got_tbry(iNH4_)) tbry_west(j,k,iNH4_)=0.1*
     &                               exp(zbry/100.)
            if (.not.got_tbry(iPhy1)) tbry_west(j,k,iPhy1)=0.04*
     &                               exp(zbry/50.)
            if (.not.got_tbry(iPhy2)) tbry_west(j,k,iPhy2)=0.06*
     &                               exp(zbry/50.)
            if (.not.got_tbry(iZoo1)) tbry_west(j,k,iZoo1)=0.04*
     &                               exp(zbry/100.)
            if (.not.got_tbry(iZoo2)) tbry_west(j,k,iZoo2)=0.04*
     &                               exp(zbry/100.)
            if (.not.got_tbry(iDet1)) tbry_west(j,k,iDet1)=0.02
            if (.not.got_tbry(iDet2)) tbry_west(j,k,iDet2)=0.02
            if (.not.got_tbry(iDON))  tbry_west(j,k,iDON) =0.5*
     &                               exp(zbry/100.)
#   ifdef NITROUS_OXIDE
            if (.not.got_tbry(iN2O)) tbry_west(j,k,iN2O)=-(0.008*
     &                              exp(zbry/100.)-0.01)
            if (tbry_west(j,k,iN2O).lt.0.) tbry_west(j,k,iN2O)=0.
#   endif
#  endif /* PISCES or BIO_NChlPZD or BIO_N2ChlPZD2 or BIO_BioEBUS */
          enddo
        enddo
      endif
# endif /* OBC_WEST */
!
# ifdef OBC_EAST
      if (EASTERN_EDGE) then
        do k=1,N
          do j=JstrR,JendR
#  ifdef TEMPERATURE
            temp=tbry_east(j,k,itemp)
#  else
            temp=25.
#  endif
            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
            xno3=1.67+0.5873*SiO4+0.0144*SiO4**2
     &           +0.0003099*SiO4**3
            if (.not.got_tbry(iNO3_)) tbry_east(j,k,iNO3_)=xno3
#  ifdef PISCES
            if (.not.got_tbry(iDIC_)) tbry_east(j,k,iDIC_)=2150.
            if (.not.got_tbry(iTAL_)) tbry_east(j,k,iTAL_)=2350.
            if (.not.got_tbry(iOXY_)) tbry_east(j,k,iOXY_)=200.
            if (.not.got_tbry(iPOC_)) tbry_east(j,k,iPOC_)=0.01
            if (.not.got_tbry(iPHY_)) tbry_east(j,k,iPHY_)=0.01
            if (.not.got_tbry(iZOO_)) tbry_east(j,k,iZOO_)=0.01
            if (.not.got_tbry(iDOC_)) tbry_east(j,k,iDOC_)=5.
            if (.not.got_tbry(iFER_)) tbry_east(j,k,iFER_)=6.e-4    
#   if ! defined key_pisces_npzd
            if (.not.got_tbry(iCAL_)) tbry_east(j,k,iCAL_)=0.01
            if (.not.got_tbry(iPO4_)) tbry_east(j,k,iPO4_)=xno3/16.
            if (.not.got_tbry(iSIL_)) tbry_east(j,k,iSIL_)=91.51
            if (.not.got_tbry(iDIA_)) tbry_east(j,k,iDIA_)=0.01
            if (.not.got_tbry(iMES_)) tbry_east(j,k,iMES_)=0.01
            if (.not.got_tbry(iDSI_)) tbry_east(j,k,iDSI_)=1.5e-3
            if (.not.got_tbry(iBFE_)) tbry_east(j,k,iBFE_)=1.e-2*5.e-6
            if (.not.got_tbry(iGOC_)) tbry_east(j,k,iGOC_)=0.01
            if (.not.got_tbry(iSFE_)) tbry_east(j,k,iSFE_)=1.e-2*5.e-6
            if (.not.got_tbry(iDFE_)) tbry_east(j,k,iDFE_)=1.e-2*5.e-6
            if (.not.got_tbry(iGSI_)) tbry_east(j,k,iGSI_)=1.e-2*0.15
            if (.not.got_tbry(iNFE_)) tbry_east(j,k,iNFE_)=1.e-2*5e-6
            if (.not.got_tbry(iNCH_)) tbry_east(j,k,iNCH_)=1.e-2*12./55.
            if (.not.got_tbry(iDCH_)) tbry_east(j,k,iDCH_)=1.e-2*12./55.
            if (.not.got_tbry(iNH4_)) tbry_east(j,k,iNH4_)=1.e-2
#   endif
#   ifdef key_ligand
            if (.not.got_tbry(iLGW_)) tbry_east(j,k,iLGW_)=1.e-3
#   endif
#   ifdef key_pisces_quota
            if (.not.got_tbry(iDON_)) tbry_east(j,k,iDON_)=5.
            if (.not.got_tbry(iDOP_)) tbry_east(j,k,iDOP_)=5.
            if (.not.got_tbry(iPON_)) tbry_east(j,k,iPON_)=1.e-2
            if (.not.got_tbry(iPOP_)) tbry_east(j,k,iPOP_)=1.e-2
            if (.not.got_tbry(iNPH_)) tbry_east(j,k,iNPH_)=1.e-2
            if (.not.got_tbry(iPPH_)) tbry_east(j,k,iPPH_)=1.e-2
            if (.not.got_tbry(iNDI_)) tbry_east(j,k,iNDI_)=1.e-2
            if (.not.got_tbry(iPDI_)) tbry_east(j,k,iPDI_)=1.e-2
            if (.not.got_tbry(iPIC_)) tbry_east(j,k,iPIC_)=1.e-2
            if (.not.got_tbry(iNPI_)) tbry_east(j,k,iNPI_)=1.e-2
            if (.not.got_tbry(iPPI_)) tbry_east(j,k,iPPI_)=1.e-2
            if (.not.got_tbry(iPFE_)) tbry_east(j,k,iPFE_)=1.e-2*5.E-6
            if (.not.got_tbry(iPCH_)) tbry_east(j,k,iPCH_)=1.e-2*12./55.
            if (.not.got_tbry(iGON_)) tbry_east(j,k,iGON_)=1.e-2
            if (.not.got_tbry(iGOP_)) tbry_east(j,k,iGOP_)=1.e-2
#   endif
#  elif defined BIO_NChlPZD
            if (.not.got_tbry(iChla)) tbry_east(j,k,iChla)=0.08
            if (.not.got_tbry(iPhy1)) tbry_east(j,k,iPhy1)=0.1
            if (.not.got_tbry(iZoo1)) tbry_east(j,k,iZoo1)=0.06
            if (.not.got_tbry(iDet1)) tbry_east(j,k,iDet1)=0.02
#  elif defined BIO_N2ChlPZD2
            if (.not.got_tbry(iNH4_)) tbry_east(j,k,iNH4_)=0.1
            if (.not.got_tbry(iChla)) tbry_east(j,k,iChla)=0.08
            if (.not.got_tbry(iPhy1)) tbry_east(j,k,iPhy1)=0.06
            if (.not.got_tbry(iZoo1)) tbry_east(j,k,iZoo1)=0.04
            if (.not.got_tbry(iDet1)) tbry_east(j,k,iDet1)=0.02
            if (.not.got_tbry(iDet2)) tbry_east(j,k,iDet2)=0.02
#  elif defined BIO_BioEBUS
            zbry=z_r(IendR,j,k)
            if (.not.got_tbry(iNO2_)) tbry_east(j,k,iNO2_)=0.05*
     &                               exp(zbry/100.)
            if (.not.got_tbry(iNH4_)) tbry_east(j,k,iNH4_)=0.1*
     &                               exp(zbry/100.)
            if (.not.got_tbry(iPhy1)) tbry_east(j,k,iPhy1)=0.04*
     &                               exp(zbry/50.)
            if (.not.got_tbry(iPhy2)) tbry_east(j,k,iPhy2)=0.06*
     &                               exp(zbry/50.)
            if (.not.got_tbry(iZoo1)) tbry_east(j,k,iZoo1)=0.04*
     &                               exp(zbry/100.)
            if (.not.got_tbry(iZoo2)) tbry_east(j,k,iZoo2)=0.04*
     &                               exp(zbry/100.)
            if (.not.got_tbry(iDet1)) tbry_east(j,k,iDet1)=0.02
            if (.not.got_tbry(iDet2)) tbry_east(j,k,iDet2)=0.02
            if (.not.got_tbry(iDON))  tbry_east(j,k,iDON) =0.5*
     &                               exp(zbry/100.)
#   ifdef NITROUS_OXIDE
            if (.not.got_tbry(iN2O)) tbry_east(j,k,iN2O)=-(0.008*
     &                               exp(zbry/100.)-0.01)
             if (tbry_east(j,k,iN2O).lt.0.) tbry_east(j,k,iN2O)=0.
#   endif
#  endif /* PISCES or BIO_NChlPZD or BIO_N2ChlPZD2 or BIO_BioEBUS */
           enddo
         enddo
       endif
# endif /* OBC_EAST */
!
# ifdef OBC_NORTH
       if (NORTHERN_EDGE) then
         do k=1,N
           do j=IstrR,IendR
#  ifdef TEMPERATURE
             temp=tbry_north(j,k,itemp)
#  else
             temp=25.
#  endif
             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
             xno3=1.67+0.5873*SiO4+0.0144*SiO4**2
     &            +0.0003099*SiO4**3
             if (.not.got_tbry(iNO3_)) tbry_north(j,k,iNO3_)=xno3
#  ifdef PISCES
            if (.not.got_tbry(iDIC_)) tbry_north(j,k,iDIC_)=2150.
            if (.not.got_tbry(iTAL_)) tbry_north(j,k,iTAL_)=2350.
            if (.not.got_tbry(iOXY_)) tbry_north(j,k,iOXY_)=200.
            if (.not.got_tbry(iPOC_)) tbry_north(j,k,iPOC_)=0.01
            if (.not.got_tbry(iPHY_)) tbry_north(j,k,iPHY_)=0.01
            if (.not.got_tbry(iZOO_)) tbry_north(j,k,iZOO_)=0.01
            if (.not.got_tbry(iDOC_)) tbry_north(j,k,iDOC_)=5.
            if (.not.got_tbry(iFER_)) tbry_north(j,k,iFER_)=6.e-4
#   if ! defined key_pisces_npzd
            if (.not.got_tbry(iCAL_)) tbry_north(j,k,iCAL_)=0.01
            if (.not.got_tbry(iPO4_)) tbry_north(j,k,iPO4_)=xno3/16.
            if (.not.got_tbry(iSIL_)) tbry_north(j,k,iSIL_)=91.51
            if (.not.got_tbry(iDIA_)) tbry_north(j,k,iDIA_)=0.01
            if (.not.got_tbry(iMES_)) tbry_north(j,k,iMES_)=0.01
            if (.not.got_tbry(iDSI_)) tbry_north(j,k,iDSI_)=1.5e-3
            if (.not.got_tbry(iBFE_)) tbry_north(j,k,iBFE_)=1.e-2*5.e-6
            if (.not.got_tbry(iGOC_)) tbry_north(j,k,iGOC_)=0.01
            if (.not.got_tbry(iSFE_)) tbry_north(j,k,iSFE_)=1.e-2*5.e-6
            if (.not.got_tbry(iDFE_)) tbry_north(j,k,iDFE_)=1.e-2*5.e-6
            if (.not.got_tbry(iGSI_)) tbry_north(j,k,iGSI_)=1.e-2*0.15
            if (.not.got_tbry(iNFE_)) tbry_north(j,k,iNFE_)=1.e-2*5e-6
            if (.not.got_tbry(iNCH_)) tbry_north(j,k,iNCH_)=1.e-2*12./55.
            if (.not.got_tbry(iDCH_)) tbry_north(j,k,iDCH_)=1.e-2*12./55.
            if (.not.got_tbry(iNH4_)) tbry_north(j,k,iNH4_)=1.e-2
#   endif             
#   ifdef key_ligand
             if (.not.got_tbry(iLGW_)) tbry_north(j,k,iLGW_)=1.e-3
#   endif
#   ifdef key_pisces_quota
             if (.not.got_tbry(iDON_)) tbry_north(j,k,iDON_)=5.
             if (.not.got_tbry(iDOP_)) tbry_north(j,k,iDOP_)=5.
             if (.not.got_tbry(iPON_)) tbry_north(j,k,iPON_)=1.e-2
             if (.not.got_tbry(iPOP_)) tbry_north(j,k,iPOP_)=1.e-2
             if (.not.got_tbry(iNPH_)) tbry_north(j,k,iNPH_)=1.e-2
             if (.not.got_tbry(iPPH_)) tbry_north(j,k,iPPH_)=1.e-2
             if (.not.got_tbry(iNDI_)) tbry_north(j,k,iNDI_)=1.e-2
             if (.not.got_tbry(iPDI_)) tbry_north(j,k,iPDI_)=1.e-2
             if (.not.got_tbry(iPIC_)) tbry_north(j,k,iPIC_)=1.e-2
             if (.not.got_tbry(iNPI_)) tbry_north(j,k,iNPI_)=1.e-2
             if (.not.got_tbry(iPPI_)) tbry_north(j,k,iPPI_)=1.e-2
             if (.not.got_tbry(iPFE_)) tbry_north(j,k,iPFE_)=1.e-2*5.E-6
             if (.not.got_tbry(iPCH_)) tbry_north(j,k,iPCH_)=1.e-2*12./55.
             if (.not.got_tbry(iGON_)) tbry_north(j,k,iGON_)=1.e-2
             if (.not.got_tbry(iGOP_)) tbry_north(j,k,iGOP_)=1.e-2
#   endif
#  elif defined BIO_NChlPZD
             if (.not.got_tbry(iChla)) tbry_north(j,k,iChla)=0.08
             if (.not.got_tbry(iPhy1)) tbry_north(j,k,iPhy1)=0.1
             if (.not.got_tbry(iZoo1)) tbry_north(j,k,iZoo1)=0.06
             if (.not.got_tbry(iDet1)) tbry_north(j,k,iDet1)=0.02
#  elif defined BIO_N2ChlPZD2
             if (.not.got_tbry(iNH4_)) tbry_north(j,k,iNH4_)=0.1
             if (.not.got_tbry(iChla)) tbry_north(j,k,iChla)=0.08
             if (.not.got_tbry(iPhy1)) tbry_north(j,k,iPhy1)=0.06
             if (.not.got_tbry(iZoo1)) tbry_north(j,k,iZoo1)=0.04
             if (.not.got_tbry(iDet1)) tbry_north(j,k,iDet1)=0.02
             if (.not.got_tbry(iDet2)) tbry_north(j,k,iDet2)=0.02
#  elif defined BIO_BioEBUS
             zbry=z_r(j,JendR,k)
             if (.not.got_tbry(iNO2_)) tbry_north(j,k,iNO2_)=0.05*
     &                                             exp(zbry/100.)
             if (.not.got_tbry(iNH4_)) tbry_north(j,k,iNH4_)=0.1*
     &                                             exp(zbry/100.)
             if (.not.got_tbry(iPhy1)) tbry_north(j,k,iPhy1)=0.04*
     &                                              exp(zbry/50.)
             if (.not.got_tbry(iPhy2)) tbry_north(j,k,iPhy2)=0.06*
     &                                              exp(zbry/50.)
             if (.not.got_tbry(iZoo1)) tbry_north(j,k,iZoo1)=0.04*
     &                                             exp(zbry/100.)
             if (.not.got_tbry(iZoo2)) tbry_north(j,k,iZoo2)=0.04*
     &                                             exp(zbry/100.)
             if (.not.got_tbry(iDet1)) tbry_north(j,k,iDet1)=0.02
             if (.not.got_tbry(iDet2)) tbry_north(j,k,iDet2)=0.02
             if (.not.got_tbry(iDON))  tbry_north(j,k,iDON) =0.5*
     &                               exp(zbry/100.)
#   ifdef NITROUS_OXIDE
            if (.not.got_tbry(iN2O)) tbry_north(j,k,iN2O)=-(0.008*
     &                                 exp(zbry/100.)-0.01)
            if (tbry_north(j,k,iN2O).lt.0.) tbry_north(j,k,iN2O)=0.
#   endif
#  endif /* PISCES or BIO_NChlPZD or BIO_N2ChlPZD2 or BIO_BioEBUS */
           enddo
         enddo
       endif
# endif /* OBC_NORTH */

# ifdef OBC_SOUTH
       if (SOUTHERN_EDGE) then
         do k=1,N
           do j=IstrR,IendR
#  ifdef TEMPERATURE
             temp=tbry_south(j,k,itemp)
#  else
             temp=25.
#  endif
             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
             xno3=1.67+0.5873*SiO4+0.0144*SiO4**2
     &            +0.0003099*SiO4**3
             if (.not.got_tbry(iNO3_)) tbry_south(j,k,iNO3_)=xno3
#  ifdef PISCES
            if (.not.got_tbry(iDIC_)) tbry_south(j,k,iDIC_)=2150.
            if (.not.got_tbry(iTAL_)) tbry_south(j,k,iTAL_)=2350.
            if (.not.got_tbry(iOXY_)) tbry_south(j,k,iOXY_)=200.
            if (.not.got_tbry(iPOC_)) tbry_south(j,k,iPOC_)=0.01
            if (.not.got_tbry(iPHY_)) tbry_south(j,k,iPHY_)=0.01
            if (.not.got_tbry(iZOO_)) tbry_south(j,k,iZOO_)=0.01
            if (.not.got_tbry(iDOC_)) tbry_south(j,k,iDOC_)=5.
            if (.not.got_tbry(iFER_)) tbry_south(j,k,iFER_)=6.e-4
#   if ! defined key_pisces_npzd
            if (.not.got_tbry(iCAL_)) tbry_south(j,k,iCAL_)=0.01
            if (.not.got_tbry(iPO4_)) tbry_south(j,k,iPO4_)=xno3/16.
            if (.not.got_tbry(iSIL_)) tbry_south(j,k,iSIL_)=91.51
            if (.not.got_tbry(iDIA_)) tbry_south(j,k,iDIA_)=0.01
            if (.not.got_tbry(iMES_)) tbry_south(j,k,iMES_)=0.01
            if (.not.got_tbry(iDSI_)) tbry_south(j,k,iDSI_)=1.5e-3
            if (.not.got_tbry(iBFE_)) tbry_south(j,k,iBFE_)=1.e-2*5.e-6
            if (.not.got_tbry(iGOC_)) tbry_south(j,k,iGOC_)=0.01
            if (.not.got_tbry(iSFE_)) tbry_south(j,k,iSFE_)=1.e-2*5.e-6
            if (.not.got_tbry(iDFE_)) tbry_south(j,k,iDFE_)=1.e-2*5.e-6
            if (.not.got_tbry(iGSI_)) tbry_south(j,k,iGSI_)=1.e-2*0.15
            if (.not.got_tbry(iNFE_)) tbry_south(j,k,iNFE_)=1.e-2*5e-6
            if (.not.got_tbry(iNCH_)) tbry_south(j,k,iNCH_)=1.e-2*12./55.
            if (.not.got_tbry(iDCH_)) tbry_south(j,k,iDCH_)=1.e-2*12./55.
            if (.not.got_tbry(iNH4_)) tbry_south(j,k,iNH4_)=1.e-2
#   endif
#   ifdef key_ligand
             if (.not.got_tbry(iLGW_)) tbry_south(j,k,iLGW_)=1.e-3
#   endif
#   ifdef key_pisces_quota
             if (.not.got_tbry(iDON_)) tbry_south(j,k,iDON_)=5.
             if (.not.got_tbry(iDOP_)) tbry_south(j,k,iDOP_)=5.
             if (.not.got_tbry(iPON_)) tbry_south(j,k,iPON_)=1.e-2
             if (.not.got_tbry(iPOP_)) tbry_south(j,k,iPOP_)=1.e-2
             if (.not.got_tbry(iNPH_)) tbry_south(j,k,iNPH_)=1.e-2
             if (.not.got_tbry(iPPH_)) tbry_south(j,k,iPPH_)=1.e-2
             if (.not.got_tbry(iNDI_)) tbry_south(j,k,iNDI_)=1.e-2
             if (.not.got_tbry(iPDI_)) tbry_south(j,k,iPDI_)=1.e-2
             if (.not.got_tbry(iPIC_)) tbry_south(j,k,iPIC_)=1.e-2
             if (.not.got_tbry(iNPI_)) tbry_south(j,k,iNPI_)=1.e-2
             if (.not.got_tbry(iPPI_)) tbry_south(j,k,iPPI_)=1.e-2
             if (.not.got_tbry(iPFE_)) tbry_south(j,k,iPFE_)=1.e-2*5.E-6
             if (.not.got_tbry(iPCH_)) tbry_south(j,k,iPCH_)=1.e-2*12./55.
             if (.not.got_tbry(iGON_)) tbry_south(j,k,iGON_)=1.e-2
             if (.not.got_tbry(iGOP_)) tbry_south(j,k,iGOP_)=1.e-2
#   endif
#  elif defined BIO_NChlPZD
             if (.not.got_tbry(iChla)) tbry_south(j,k,iChla)=0.08
             if (.not.got_tbry(iPhy1)) tbry_south(j,k,iPhy1)=0.1
             if (.not.got_tbry(iZoo1)) tbry_south(j,k,iZoo1)=0.06
             if (.not.got_tbry(iDet1)) tbry_south(j,k,iDet1)=0.02
#  elif defined BIO_N2ChlPZD2
             if (.not.got_tbry(iNH4_)) tbry_south(j,k,iNH4_)=0.1
             if (.not.got_tbry(iChla)) tbry_south(j,k,iChla)=0.08
             if (.not.got_tbry(iPhy1)) tbry_south(j,k,iPhy1)=0.06
             if (.not.got_tbry(iZoo1)) tbry_south(j,k,iZoo1)=0.04
             if (.not.got_tbry(iDet1)) tbry_south(j,k,iDet1)=0.02
             if (.not.got_tbry(iDet2)) tbry_south(j,k,iDet2)=0.02
#  elif defined BIO_BioEBUS
             zbry=z_r(j,JstrR,k)
            if (.not.got_tbry(iNO2_)) tbry_south(j,k,iNO2_)=0.05*
     &                                     exp(zbry/100.)
            if (.not.got_tbry(iNH4_)) tbry_south(j,k,iNH4_)=0.1*
     &                                     exp(zbry/100.)
            if (.not.got_tbry(iPhy1)) tbry_south(j,k,iPhy1)=0.04*
     &                                     exp(zbry/50.)
            if (.not.got_tbry(iPhy2)) tbry_south(j,k,iPhy2)=0.06*
     &                                     exp(zbry/50.)
            if (.not.got_tbry(iZoo1)) tbry_south(j,k,iZoo1)=0.04*
     &                                     exp(zbry/100.)
            if (.not.got_tbry(iZoo2)) tbry_south(j,k,iZoo2)=0.04*
     &                                     exp(zbry/100.)
            if (.not.got_tbry(iDet1)) tbry_south(j,k,iDet1)=0.02
            if (.not.got_tbry(iDet2)) tbry_south(j,k,iDet2)=0.02
            if (.not.got_tbry(iDON))  tbry_south(j,k,iDON) =0.5*
     &                                     exp(zbry/100.)
#   ifdef NITROUS_OXIDE
            if (.not.got_tbry(iN2O)) tbry_south(j,k,iN2O)=-(0.008*
     &                              exp(zbry/100.)-0.01)
              if (tbry_south(j,k,iN2O).lt.0.) then
                tbry_south(j,k,iN2O)=0.
              endif
#   endif
#  endif /* PISCES or BIO_NChlPZD or BIO_N2ChlPZD2 or BIO_BioEBUS */
           enddo
         enddo
       endif
# endif /* OBC_SOUTH */
      return
      end
#else
      subroutine empty_analytical
      return
      end
#endif /* BIOLOGY && defined T_FRC_BRY */


