#include "cppdefs.h"
      MODULE mod_tides
#if defined SSH_TIDES || defined UV_TIDES || defined POT_TIDES
!
!svn $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2020 The ROMS/TOMS Group                         !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.txt                                              !
!=======================================================================
!                                                                      !
!  Tidal Components:                                                   !
!                                                                      !
!  Each of the following arrays has a dimension in tidal components    !
!  classified by period:                                               !
!                                                                      !
!    semi-diurnal:  M2, S2, N2, K2  (12.42, 12.00, 12.66, 11.97h)      !
!         diurnal:  K1, O1, P1, Q1  (23.93, 25.82, 24.07, 26.87h)      !
!                                                                      !
!  and other longer periods. The order of these tidal components is    !
!  irrelevant here.  The number of components to USE is depends on     !
!  the regional application.                                           !
!                                                                      !
!  CosOmega     Cosine tidal harmonics for current omega(t).           !
!  SinOmega     Sine tidal harmonics for current omega(t).             !
!  SSH_Tamp     Tidal elevation amplitude (m) at RHO-points.           !
!  SSH_Tphase   Tidal elevation phase (degrees/360) at RHO-points.     !
!  Tperiod      Tidal period (s).                                      !
!  UV_Tangle    Tidal current angle (radians; counterclockwise         !
!                 from EAST and rotated to curvilinear grid) at        !
!                 RHO-points.                                          !
!  UV_Tmajor    Maximum tidal current: tidal ellipse major axis        !
!                 (m/s) at RHO-points.                                 !
!  UV_Tminor    Minimum tidal current: tidal ellipse minor axis        !
!                 (m/s) at RHO-points.                                 !
!  UV_Tphase    Tidal current phase (degrees/360) at RHO-points.       !
!                                                                      !
# if defined AVERAGES && defined AVERAGES_DETIDE
!                                                                      !
!  Detided time-averaged fields via least-squares fitting. Notice that !
!  the harmonics for the state variable have an extra dimension of     !
!  size (0:2*NTC) to store several terms:                              !
!                                                                      !
!               index 0               mean term (accumulated sum)      !
!                     1:NTC           accumulated sine terms           !
!                     NTC+1:2*NTC     accumulated cosine terms         !
!                                                                      !
!  CosW_avg     Current time-average window COS(omega(k)*t).           !
!  CosW_sum     Time-accumulated COS(omega(k)*t).                      !
!  SinW_avg     Current time-average window SIN(omega(k)*t).           !
!  SinW_sum     Time-accumulated SIN(omega(k)*t).                      !
!  CosWCosW     Time-accumulated COS(omega(k)*t)*COS(omega(l)*t).      !
!  SinWSinW     Time-accumulated SIN(omega(k)*t)*SIN(omega(l)*t).      !
!  SinWCosW     Time-accumulated SIN(omega(k)*t)*COS(omega(l)*t).      !
!                                                                      !
!  ubar_detided Time-averaged and detided 2D u-momentum.               !
!  ubar_tide    Time-accumulated 2D u-momentum tide harmonics.         !
!  vbar_detided Time-averaged and detided 2D v-momentum.               !
!  vbar_tide    Time-accumulated 2D v-momentum tide harmonics.         !
!  zeta_detided Time-averaged and detided free-surface.                !
!  zeta_tide    Time-accumulated free-surface tide harmonics.          !
#  ifdef SOLVE3D
!  t_detided    Time-averaged and detided tracers (T,S).               !
!  t_tide       Time-accumulated 3D tracers (T,S) tide harmonics.      !
!  u_detided    Time-averaged and detided 3D u-momentum.               !
!  u_tide       Time-accumulated 3D u-momentum tide harmonics.         !
!  v_detided    Time-averaged and detided 3D v-momentum.               !
!  v_tide       Time-accumulated 3D v-momentum tide harmonics.         !
#  endif
!                                                                      !
# endif
# ifdef UV_WAVEDRAG
!  We need to accumulate a detided estimate of the bottom velocities   !
!  which is the subtracted from the bottom velocity before the         !
!  wave drag due to conversion of tidal energy into internal tides     !
!  is applied.                                                         !
!  sum_ubot      Time-accumulated bottom u-momentum.                   !
!  sum_vbot      Time-accumulated bottom v-momentum.                   !
!  filt_ubot     Time-mean bottom u-momentum.                          !
!  filt_vbot     Time-mean bottom v-momentum.                          !
# endif
!=======================================================================
!
        USE mod_kinds
        USE mod_param
        USE mod_stepping

        implicit none

        TYPE T_TIDES

          real(r8), pointer :: Tperiod(:)
# if defined AVERAGES && defined AVERAGES_DETIDE
          real(r8), pointer :: CosOmega(:)
          real(r8), pointer :: SinOmega(:)
          real(r8), pointer :: CosW_avg(:)
          real(r8), pointer :: CosW_sum(:)
          real(r8), pointer :: SinW_avg(:)
          real(r8), pointer :: SinW_sum(:)
          real(r8), pointer :: CosWCosW(:,:)
          real(r8), pointer :: SinWSinW(:,:)
          real(r8), pointer :: SinWCosW(:,:)
# endif
# if defined SSH_TIDES
          real(r8), pointer :: SSH_Tamp(:,:,:)
          real(r8), pointer :: SSH_Tphase(:,:,:)
# endif
# if defined UV_TIDES
          real(r8), pointer :: UV_Tangle(:,:,:)
          real(r8), pointer :: UV_Tmajor(:,:,:)
          real(r8), pointer :: UV_Tminor(:,:,:)
          real(r8), pointer :: UV_Tphase(:,:,:)
# endif
# if defined POT_TIDES
          real(r8), pointer :: POT_Tamp(:,:,:)
          real(r8), pointer :: POT_Tphase(:,:,:)
          real(r8), pointer :: Ptide(:,:)
# endif
# if defined TIDES_ASTRO
          real(r8), pointer :: Vu_sat(:,:,:)
          real(r8), pointer :: f_sat(:,:,:)
# endif
# if defined AVERAGES && defined AVERAGES_DETIDE
          real(r8), pointer :: ubar_detided(:,:)
          real(r8), pointer :: ubar_tide(:,:,:)

          real(r8), pointer :: vbar_detided(:,:)
          real(r8), pointer :: vbar_tide(:,:,:)

          real(r8), pointer :: zeta_detided(:,:)
          real(r8), pointer :: zeta_tide(:,:,:)
#  ifdef SOLVE3D
          real(r8), pointer :: t_detided(:,:,:,:)
          real(r8), pointer :: t_tide(:,:,:,:,:)

          real(r8), pointer :: u_detided(:,:,:)
          real(r8), pointer :: u_tide(:,:,:,:)

          real(r8), pointer :: v_detided(:,:,:)
          real(r8), pointer :: v_tide(:,:,:,:)
#  endif
# endif
# ifdef UV_WAVEDRAG
          real(r8), pointer :: sum_ubot(:,:)
          real(r8), pointer :: sum_vbot(:,:)
          real(r8), pointer :: filt_ubot(:,:)
          real(r8), pointer :: filt_vbot(:,:)
# endif

        END TYPE T_TIDES

        TYPE (T_TIDES), allocatable :: TIDES(:)

# if defined TIDES_ASTRO
        integer, parameter :: MAX_SAT = 10
        integer, allocatable  :: doodson(:,:)
        real(r8), allocatable :: phase_corr(:)
        integer, allocatable  :: num_sat(:)
        integer, allocatable  :: sat_doodson(:,:,:)
        real(r8), allocatable :: sat_phase_corr(:,:)
        real(r8), allocatable :: sat_amp(:,:)
        real(r8), allocatable :: sat_flag(:,:)
# endif
# ifdef UV_WAVEDRAG
        integer :: nstep_hour
        integer :: hour_steps = 0
        real(r8), parameter :: drag_scale = 0.375
# endif

      CONTAINS

      SUBROUTINE allocate_tides (ng, LBi, UBi, LBj, UBj)
!
!=======================================================================
!                                                                      !
!  This routine allocates all variables in the module for all nested   !
!  grids.                                                              !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_iounits
      USE mod_ncparam
      USE mod_netcdf
      USE mod_scalars
      USE mod_stepping
!
      USE strings_mod, ONLY : FoundError
!
! Inported variable declarations.
!
      integer, intent(in) :: ng, LBi, UBi, LBj, UBj
!
!  Local variable declarations.
!
      logical :: foundit

      integer :: Nfiles, Vid, i, ifile, mg, nvatt, nvdim

      real(r8) :: size2d
!
!-----------------------------------------------------------------------
!  Allocate module variables.
!-----------------------------------------------------------------------
!
!  Inquire about the maximum number of tidal components. Notice that
!  currently we only support nested applications where the tidal
!  forcing is applied to the main coarser grid (RefineScale(ng)=0)
!  and the other grids get the tidal forcing from the contact areas.
!
      IF (LprocessTides(ng)) THEN
        MTC=0
        foundit=.FALSE.
        CALL netcdf_inq_var (ng, iNLM, TIDE(ng)%name,                   &
     &                       MyVarName = TRIM(Vname(1,idTper)),         &
     &                       SearchVar = foundit,                       &
     &                       VarID = Vid,                               &
     &                       nVardim = nvdim,                           &
     &                       nVarAtt = nvatt)
        IF (FoundError(exit_flag, NoError, __LINE__,                    &
     &                 __FILE__//", allocate_tides")) RETURN
!
!  Set maximum number of tidal components.  Allocate and initialize
!  TIDE I/O structure. Notice that in nested applications, all the
!  nested grids need to have the same number of tidal component.
!
        IF (foundit) THEN
          MTC=MAX(MTC,var_Dsize(1))            ! first dimension
          DO mg=1,Ngrids
            NTC(mg)=var_Dsize(1)
          END DO
        END IF
      END IF
!
!  Allocate structure.
!
      IF (ng.eq.1) allocate ( TIDES(Ngrids) )
!
!  Set horizontal array size.
!
      size2d=REAL((UBi-LBi+1)*(UBj-LBj+1),r8)
!
!  Allocate tidal forcing variables.
!
      allocate ( TIDES(ng) % Tperiod(MTC)  )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)

# if defined SSH_TIDES
      allocate ( TIDES(ng) % SSH_Tamp(LBi:UBi,LBj:UBj,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)*size2d

      allocate ( TIDES(ng) % SSH_Tphase(LBi:UBi,LBj:UBj,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)*size2d
# endif

# if defined UV_TIDES
      allocate ( TIDES(ng) % UV_Tangle(LBi:UBi,LBj:UBj,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)*size2d

      allocate ( TIDES(ng) % UV_Tmajor(LBi:UBi,LBj:UBj,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)*size2d

      allocate ( TIDES(ng) % UV_Tminor(LBi:UBi,LBj:UBj,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)*size2d

      allocate ( TIDES(ng) % UV_Tphase(LBi:UBi,LBj:UBj,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)*size2d
# endif

# if defined POT_TIDES
      allocate ( TIDES(ng) % POT_Tamp(LBi:UBi,LBj:UBj,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)*size2d
      allocate ( TIDES(ng) % POT_Tphase(LBi:UBi,LBj:UBj,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)*size2d
      allocate ( TIDES(ng) % Ptide(LBi:UBi,LBj:UBj) )
      Dmem(ng)=Dmem(ng)+size2d
# endif

# ifdef TIDES_ASTRO
      allocate( TIDES(ng) % Vu_sat(LBi:UBi,LBj:UBj,MTC))
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)*size2d
      allocate( TIDES(ng) % f_sat(LBi:UBi,LBj:UBj,MTC))
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)*size2d
# endif

# if defined TIDES_ASTRO
      allocate ( doodson(6,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC*6,r8)
      allocate ( phase_corr(MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)
      allocate ( num_sat(MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)
      allocate ( sat_doodson(3,MAX_SAT,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(3*MAX_SAT*MTC,r8)
      allocate ( sat_phase_corr(MAX_SAT,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MAX_SAT*MTC,r8)
      allocate ( sat_amp(MAX_SAT,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MAX_SAT*MTC,r8)
      allocate ( sat_flag(MAX_SAT,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MAX_SAT*MTC,r8)
# endif

# if defined AVERAGES && defined AVERAGES_DETIDE
!
!  Allocate variables used for the least-squares detiding of
!  time-averaged fields.
!
      allocate ( TIDES(ng) % CosOmega(MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)

      allocate ( TIDES(ng) % SinOmega(MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)

      allocate ( TIDES(ng) % CosW_avg(MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)

      allocate ( TIDES(ng) % CosW_sum(MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)

      allocate ( TIDES(ng) % SinW_avg(MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)

      allocate ( TIDES(ng) % SinW_sum(MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC,r8)

      allocate ( TIDES(ng) % CosWCosW(MTC,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC*MTC,r8)

      allocate ( TIDES(ng) % SinWSinW(MTC,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC*MTC,r8)

      allocate ( TIDES(ng) % SinWCosW(MTC,MTC) )
      Dmem(ng)=Dmem(ng)+REAL(MTC*MTC,r8)

      IF (Aout(idFsuD,ng)) THEN
        allocate ( TIDES(ng) % zeta_detided(LBi:UBi,LBj:UBj) )
        Dmem(ng)=Dmem(ng)+size2d

        allocate ( TIDES(ng) % zeta_tide(LBi:UBi,LBj:UBj,0:2*MTC) )
        Dmem(ng)=Dmem(ng)+REAL(2*MTC+1,r8)*size2d
      END IF

      IF (Aout(idu2dD,ng)) THEN
        allocate ( TIDES(ng) % ubar_detided(LBi:UBi,LBj:UBj) )

        allocate ( TIDES(ng) % ubar_tide(LBi:UBi,LBj:UBj,0:2*MTC) )
        Dmem(ng)=Dmem(ng)+REAL(2*MTC+1,r8)*size2d
      END IF

      IF (Aout(idv2dD,ng)) THEN
        allocate ( TIDES(ng) % vbar_detided(LBi:UBi,LBj:UBj) )

        allocate ( TIDES(ng) % vbar_tide(LBi:UBi,LBj:UBj,0:2*MTC) )
        Dmem(ng)=Dmem(ng)+REAL(2*MTC+1,r8)*size2d
      END IF

#  ifdef SOLVE3D
      IF (Aout(idu3dD,ng)) THEN
        allocate ( TIDES(ng) % u_detided(LBi:UBi,LBj:UBj,N(ng)) )
        Dmem(ng)=Dmem(ng)+REAL(N(ng),r8)*size2d

        allocate ( TIDES(ng) % u_tide(LBi:UBi,LBj:UBj,N(ng),0:2*MTC) )
        Dmem(ng)=Dmem(ng)+REAL(N(ng)*(2*MTC+1),r8)*size2d
      END IF

      IF (Aout(idv3dD,ng)) THEN
        allocate ( TIDES(ng) % v_detided(LBi:UBi,LBj:UBj,N(ng)) )
        Dmem(ng)=Dmem(ng)+REAL(N(ng),r8)*size2d

        allocate ( TIDES(ng) % v_tide(LBi:UBi,LBj:UBj,N(ng),0:2*MTC) )
        Dmem(ng)=Dmem(ng)+REAL(N(ng)*(2*MTC+1),r8)*size2d
      END IF

      IF (ANY(Aout(idTrcD(:),ng))) THEN
        allocate ( TIDES(ng) % t_detided(LBi:UBi,LBj:UBj,N(ng),NAT) )
        Dmem(ng)=Dmem(ng)+REAL(N(ng)*NAT,r8)*size2d

        allocate ( TIDES(ng) % t_tide(LBi:UBi,LBj:UBj,N(ng),            &
     &                                0:2*MTC,NAT) )
        Dmem(ng)=Dmem(ng)+REAL(N(ng)*(2*MTC+1)*NAT,r8)*size2d
      END IF
#  endif
# endif
# ifdef UV_WAVEDRAG
      allocate ( TIDES(ng) % sum_ubot(LBi:UBi,LBj:UBj) )
      Dmem(ng)=Dmem(ng)+size2d
      allocate ( TIDES(ng) % sum_vbot(LBi:UBi,LBj:UBj) )
      Dmem(ng)=Dmem(ng)+size2d
      allocate ( TIDES(ng) % filt_ubot(LBi:UBi,LBj:UBj) )
      Dmem(ng)=Dmem(ng)+size2d
      allocate ( TIDES(ng) % filt_vbot(LBi:UBi,LBj:UBj) )
      Dmem(ng)=Dmem(ng)+size2d
# endif

      RETURN
      END SUBROUTINE allocate_tides

      SUBROUTINE initialize_tides (ng, tile)
!
!=======================================================================
!                                                                      !
!  This routine initialize all variables in the module using first     !
!  touch distribution policy. In shared-memory configuration, this     !
!  operation actually performs propagation of the  "shared arrays"     !
!  across the cluster, unless another policy is specified to           !
!  override the default.                                               !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_ncparam
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile
!
!  Local variable declarations.
!
      integer :: Imin, Imax, Jmin, Jmax
      integer :: i, itide, itrc, j, jtide, k

      real(r8), parameter :: IniVal = 0.0_r8

# include "set_bounds.h"
!
!  Set array initialization range.
!
# ifdef DISTRIBUTE
      Imin=BOUNDS(ng)%LBi(tile)
      Imax=BOUNDS(ng)%UBi(tile)
      Jmin=BOUNDS(ng)%LBj(tile)
      Jmax=BOUNDS(ng)%UBj(tile)
# else
      IF (DOMAIN(ng)%Western_Edge(tile)) THEN
        Imin=BOUNDS(ng)%LBi(tile)
      ELSE
        Imin=Istr
      END IF
      IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
        Imax=BOUNDS(ng)%UBi(tile)
      ELSE
        Imax=Iend
      END IF
      IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
        Jmin=BOUNDS(ng)%LBj(tile)
      ELSE
        Jmin=Jstr
      END IF
      IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
        Jmax=BOUNDS(ng)%UBj(tile)
      ELSE
        Jmax=Jend
      END IF
# endif
!
!-----------------------------------------------------------------------
!  Initialize module variables.
!-----------------------------------------------------------------------
!
!  Initialize tidal forcing variables.
!
      IF (DOMAIN(ng)%SouthWest_Test(tile)) THEN
        DO itide=1,MTC
          TIDES(ng) % Tperiod(itide) = IniVal
        END DO
      END IF

      DO itide=1,MTC
# if defined SSH_TIDES
        DO j=Jmin,Jmax
          DO i=Imin,Imax
            TIDES(ng) % SSH_Tamp(i,j,itide) = IniVal
            TIDES(ng) % SSH_Tphase(i,j,itide) = IniVal
          END DO
        END DO
# endif
# if defined UV_TIDES
        DO j=Jmin,Jmax
          DO i=Imin,Imax
            TIDES(ng) % UV_Tangle(i,j,itide) = IniVal
            TIDES(ng) % UV_Tmajor(i,j,itide) = IniVal
            TIDES(ng) % UV_Tminor(i,j,itide) = IniVal
            TIDES(ng) % UV_Tphase(i,j,itide) = IniVal
          END DO
        END DO
# endif
      END DO

# if defined AVERAGES && defined AVERAGES_DETIDE
!
!  Initialize cariables used for the least-squares detiding of
!  time-averaged fields.
!
      IF (DOMAIN(ng)%SouthWest_Test(tile)) THEN
        DO jtide=1,MTC
          TIDES(ng) % CosOmega(jtide) = IniVal
          TIDES(ng) % SinOmega(jtide) = IniVal
          TIDES(ng) % CosW_avg(jtide) = IniVal
          TIDES(ng) % CosW_sum(jtide) = IniVal
          TIDES(ng) % SinW_avg(jtide) = IniVal
          TIDES(ng) % SinW_sum(jtide) = IniVal
          DO itide=1,MTC
            TIDES(ng) % CosWCosW(itide,jtide) = IniVal
            TIDES(ng) % SinWSinW(itide,jtide) = IniVal
            TIDES(ng) % SinWCosW(itide,jtide) = IniVal
          END DO
        END DO
      END IF

      IF (Aout(idFsuD,ng)) THEN
        DO j=Jmin,Jmax
          DO i=Imin,Imax
            TIDES(ng) % zeta_detided(i,j) = IniVal
          END DO
        END DO
        DO itide=0,2*MTC
          DO j=Jmin,Jmax
            DO i=Imin,Imax
              TIDES(ng) % zeta_tide(i,j,itide) = IniVal
            END DO
          END DO
        END DO
      END IF

      IF (Aout(idu2dD,ng)) THEN
        DO j=Jmin,Jmax
          DO i=Imin,Imax
            TIDES(ng) % ubar_detided(i,j) = IniVal
          END DO
        END DO
        DO itide=0,2*MTC
          DO j=Jmin,Jmax
            DO i=Imin,Imax
              TIDES(ng) % ubar_tide(i,j,itide) = IniVal
            END DO
          END DO
        END DO
      END IF
# endif

# if defined POT_TIDES
      DO itide=1,MTC
        DO j=Jmin,Jmax
          DO i=Imin,Imax
            TIDES(ng) % POT_Tamp(i,j,itide) = IniVal
            TIDES(ng) % POT_Tphase(i,j,itide) = IniVal
          END DO
        END DO
      END DO
      DO j=Jmin,Jmax
        DO i=Imin,Imax
          TIDES(ng) % Ptide(i,j) = IniVal
        END DO
      END DO
# endif

# ifdef TIDES_ASTRO
      DO itide=1,MTC
        DO j=Jmin,Jmax
          DO i=Imin,Imax
            TIDES(ng) % Vu_sat(i,j,itide) = 0.0_r8
            TIDES(ng) % f_sat(i,j,itide) = 1.0_r8
          END DO
        END DO
      END DO
# endif

# if defined TIDES_ASTRO
! Assumes tides are in the order Q1, O1, P1, K1, N2, M2, S2, K2

      doodson = reshape( (/ 1, -2, 0, 1, 0, 0,                          &
     &                      1, -1, 0, 0, 0, 0,                          &
     &                      1, 1, -2, 0, 0, 0,                          &
     &                      1, 1,  0, 0, 0, 0,                          &
     &                      2, -1, 0, 1, 0, 0,                          &
     &                      2, 0,  0, 0, 0, 0,                          &
     &                      2, 2, -2, 0, 0, 0,                          &
     &                      2, 2,  0, 0, 0, 0/), (/ 6,MTC /) )

      phase_corr = (/ -0.25, -0.25, -0.25, -0.75,                       &
     &                             0.0, 0.0, 0.0, 0.0 /)

      num_sat = (/ 10, 8, 6, 10, 4, 9, 3, 5 /)

      sat_doodson(:,:,1) = reshape( (/ -2, -3, 0, -2, -2, 0,            &
     &                                 -1, -2, 0, -1, -1, 0,            &
     &                                 -1,  0, 0,  0, -2, 0,            &
     &                                 -1,  0, 1,  0, -1, 0,            &
     &                                  1,  0, 0,  2,  0, 0 /),         &
     &                               (/ 3,MAX_SAT /) )
      sat_doodson(:,:,2) = reshape( (/ -1, 0, 0, 0, -2, 0,              &
     &                                 0, -1, 0, 1, -1, 0,              &
     &                                 1,  0, 0, 1,  1, 0,              &
     &                                 2,  0, 0, 2,  1, 0,              &
     &                                 -99, -99, -99, -99, -99, -99 /), &
     &                                (/ 3,MAX_SAT /) )
      sat_doodson(:,:,3) = reshape( (/ 0, -2, 0, 0, -1, 0,              &
     &                                 0,  0, 2, 1,  0, 0,              &
     &                                 2,  0, 0, 2,  1, 0,              &
     &                                 -99, -99, -99, -99, -99, -99,    &
     &                                 -99, -99, -99, -99, -99, -99 /), &
     &                                (/ 3,MAX_SAT /) )
      sat_doodson(:,:,4) = reshape( (/ -2, -1, 0, -1, -1, 0,            &
     &                                 -1,  0, 0, -1,  1, 0,            &
     &                                 0, -2, 0,  0, -1, 0,             &
     &                                 0,  1, 0,  0,  2, 0,             &
     &                                 1,  0, 0,  1,  1, 0 /),          &
     &                                (/ 3,MAX_SAT /) )
      sat_doodson(:,:,5) = reshape( (/ -2, -2, 0, -1,  0, 1,            &
     &                                  0, -2, 0,  0, -1, 0,            &
     &                                 -99, -99, -99, -99, -99, -99,    &
     &                                 -99, -99, -99, -99, -99, -99,    &
     &                                 -99, -99, -99, -99, -99, -99 /), &
     &                                (/ 3,MAX_SAT /) )
      sat_doodson(:,:,6) = reshape( (/ -1, -1, 0, -1,  0, 0,            &
     &                                  0, -2, 0,  0, -1, 0,            &
     &                                  1, -1, 0,  1,  0, 0,            &
     &                                  1,  1, 0,  2,  0, 0,            &
     &                                  2,  1, 0, -99, -99, -99 /),     &
     &                                (/ 3,MAX_SAT /) )
      sat_doodson(:,:,7) = reshape( (/ 0, -1, 0, 1, 0, 0,               &
     &                                 2,  0, 0, -99, -99, -99,         &
     &                                 -99, -99, -99, -99, -99, -99,    &
     &                                 -99, -99, -99, -99, -99, -99,    &
     &                                 -99, -99, -99, -99, -99, -99 /), &
     &                                (/ 3,MAX_SAT /) )
      sat_doodson(:,:,8) = reshape( (/ -1, 0, 0, -1, 1, 0,              &
     &                                 0, -1, 0,  0, 1, 0,              &
     &                                 0,  2, 0, -99, -99, -99,         &
     &                                 -99, -99, -99, -99, -99, -99,    &
     &                                 -99, -99, -99, -99, -99, -99 /), &
     &                                (/ 3,MAX_SAT /) )

      sat_phase_corr = reshape(                                         &
     & (/ 0.5, 0.5, 0.75, 0.75, 0.75, 0.5, 0.0, 0.0, 0.75, 0.5,         &
     &    0.25, 0.5, 0.0, 0.25, 0.75, 0.25, 0.5, 0.5, -99., -99.,       &
     &    0.0, 0.5, 0.5, 0.75, 0.5, 0.5, -99., -99., -99., -99.,        &
     &    0.0, 0.75, 0.25, 0.75, 0.0, 0.5, 0.0, 0.5, 0.25, 0.25,        &
     &    0.5, 0.0, 0.0, 0.5, -99., -99., -99., -99., -99., -99.,       &
     &    0.75, 0.75, 0.0, 0.5, 0.25, 0.75, 0.75, 0.0, 0.0, -99.,       &
     &    0.0, 0.75, 0.0, -99., -99., -99., -99., -99., -99., -99.,     &
     &    0.75, 0.75, 0.5, 0.0, 0.0, -99., -99., -99., -99., -99. /),   &
     &      (/ MAX_SAT,MTC /) )

      sat_amp = reshape(                                                &
     &      (/ 0.0007, 0.0038, 0.0010, 0.0115, 0.0292,                  &
     &         0.0057, 0.0008, 0.1884, 0.0018, 0.0028,                  &
     &         0.0003, 0.0058, 0.1885, 0.0004, 0.0029,                  &
     &         0.0004, 0.0064, 0.0010, 9999.0, 9999.0,                  &
     &         0.0008, 0.0112, 0.0004, 0.0004, 0.0015,                  &
     &         0.0003, 9999.0, 9999.0, 9999.0, 9999.0,                  &
     &         0.0002, 0.0001, 0.0007, 0.0001, 0.0001,                  &
     &         0.0198, 0.1356, 0.0029, 0.0002, 0.0001,                  &
     &         0.0039, 0.0008, 0.0005, 0.0373, 9999.0,                  &
     &         9999.0, 9999.0, 9999.0, 9999.0, 9999.0,                  &
     &         0.0001, 0.0004, 0.0005, 0.0373, 0.0001,                  &
     &         0.0009, 0.0002, 0.0006, 0.0002, 9999.0,                  &
     &         0.0022, 0.0001, 0.0001, 9999.0, 9999.0,                  &
     &         9999.0, 9999.0, 9999.0, 9999.0, 9999.0,                  &
     &         0.0024, 0.0004, 0.0128, 0.2980, 0.0324,                  &
     &         9999.0, 9999.0, 9999.0, 9999.0, 9999.0 /),               &
     &      (/ MAX_SAT,MTC /) )

      sat_flag = reshape(                                               &
     &       (/ 0, 0, 1, 1, 1, 0, 0, 0, 1, 0,                           &
     &          1, 0, 0, 1, 1, 1, 0, 0, -99, -99,                       &
     &          0, 0, 0, 1, 0, 0, -99, -99, -99, -99,                   &
     &          0, 1, 1, 1, 0, 0, 0, 0, 1, 1,                           &
     &          0, 0, 0, 0, -99, -99, -99, -99, -99, -99,               &
     &          2, 2, 0, 0, 2, 2, 2, 0, 0, -99,                         &
     &          0, 2, 0, -99, -99, -99, -99, -99, -99, -99,             &
     &          2, 2, 0, 0, 0, -99, -99, -99, -99, -99 /),              &
     &      (/ MAX_SAT,MTC /) )

# endif
# if defined AVERAGES_DETIDE && defined AVERAGES
      IF (Aout(idv2dD,ng)) THEN
        DO j=Jmin,Jmax
          DO i=Imin,Imax
            TIDES(ng) % vbar_detided(i,j) = IniVal
          END DO
        END DO
        DO itide=0,2*MTC
          DO j=Jmin,Jmax
            DO i=Imin,Imax
              TIDES(ng) % vbar_tide(i,j,itide) = IniVal
            END DO
          END DO
        END DO
      END IF

#  ifdef SOLVE3D
      IF (Aout(idu3dD,ng)) THEN
        DO k=1,N(ng)
          DO j=Jmin,Jmax
            DO i=Imin,Imax
              TIDES(ng) % u_detided(i,j,k) = IniVal
            END DO
          END DO
        END DO
        DO itide=0,2*MTC
          DO k=1,N(ng)
            DO j=Jmin,Jmax
              DO i=Imin,Imax
                TIDES(ng) % u_tide(i,j,k,itide) = IniVal
              END DO
            END DO
          END DO
        END DO
      END IF

      IF (Aout(idv3dD,ng)) THEN
        DO k=1,N(ng)
          DO j=Jmin,Jmax
            DO i=Imin,Imax
              TIDES(ng) % v_detided(i,j,k) = IniVal
            END DO
          END DO
        END DO
        DO itide=0,2*MTC
          DO k=1,N(ng)
            DO j=Jmin,Jmax
              DO i=Imin,Imax
                TIDES(ng) % v_tide(i,j,k,itide) = IniVal
              END DO
            END DO
          END DO
        END DO
      END IF

      IF (ANY(Aout(idTrcD(:),ng))) THEN
        DO itrc=1,NAT
          DO k=1,N(ng)
            DO j=Jmin,Jmax
              DO i=Imin,Imax
                TIDES(ng) % t_detided(i,j,k,itrc) = IniVal
              END DO
            END DO
          END DO
        END DO
        DO itrc=1,NAT
          DO itide=0,2*MTC
            DO k=1,N(ng)
              DO j=Jmin,Jmax
                DO i=Imin,Imax
                  TIDES(ng) % t_tide(i,j,k,itide,itrc) = IniVal
                END DO
              END DO
            END DO
          END DO
        END DO
      END IF
#  endif
# endif
# ifdef UV_WAVEDRAG
      DO j=Jmin,Jmax
        DO i=Imin,Imax
          TIDES(ng) % sum_ubot(i,j) = IniVal
          TIDES(ng) % sum_vbot(i,j) = IniVal
          TIDES(ng) % filt_ubot(i,j) = IniVal
          TIDES(ng) % filt_vbot(i,j) = IniVal
        END DO
      END DO
# endif

      RETURN
      END SUBROUTINE initialize_tides
#endif

#ifdef TIDES_ASTRO
      SUBROUTINE tide_astro(ttime,vux,fx,xlat,ng,tile,                  &
     &                      LBi, UBi, LBj, UBj)
      USE mod_scalars
      implicit none

      integer, intent(in)  :: ng, tile
      integer, intent(in) :: LBi, UBi, LBj, UBj
      real(r8), intent(in)  :: ttime
# ifdef ASSUMED_SHAPE
      real(r8), intent(in)  :: xlat(LBi:,LBj:)
      real(r8), intent(out) :: vux(LBi:,LBj:,:)
      real(r8), intent(out) :: fx(LBi:,LBj:,:)
# else
      real(r8), intent(in)  :: xlat(LBi:UBi,LBj:UBj)
      real(r8), intent(out) :: vux(LBi:UBi,LBj:UBj,MTC)
      real(r8), intent(out) :: fx(LBi:UBi,LBj:UBj,MTC)
# endif

! Local variables
      real(r8) :: slat, hh, d1, tau, dtau_2, freq(MTC)
      real(r8) :: h, pp, s, p, enp, dh, dpp, ds, dp2, dnp
      real(r8) :: uu, twopi, sumc, sums, v, rr, vdbl, uudbl
      integer  :: i, iv, j, k, iuu, itide, intdays
      integer  :: Imin, Imax, Jmin, Jmax
      real(r8), parameter :: eps = 1.e-20_r8
!
!  Set array initialization range.
!
#include "set_bounds.h"

#ifdef _OPENMP
      IF (DOMAIN(ng)%Western_Edge(tile)) THEN
        Imin=BOUNDS(ng)%LBi(tile)
      ELSE
        Imin=Istr
      END IF
      IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
        Imax=BOUNDS(ng)%UBi(tile)
      ELSE
        Imax=Iend
      END IF
      IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
        Jmin=BOUNDS(ng)%LBj(tile)
      ELSE
        Jmin=Jstr
      END IF
      IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
        Jmax=BOUNDS(ng)%UBj(tile)
      ELSE
        Jmax=Jend
      END IF
#else
      Imin=BOUNDS(ng)%LBi(tile)
      Imax=BOUNDS(ng)%UBi(tile)
      Jmin=BOUNDS(ng)%LBj(tile)
      Jmax=BOUNDS(ng)%UBj(tile)
#endif
!
!***********************************************************************
!*  THIS SUBROUTINE CALCULATES V (ASTRONOMICAL PHASE ARGUMENT), U AND F
!*  (NODAL MODULATION PHASE AND AMPLITUDE CORRECTIONS) FOR ALL CONSTITU-
!*  ENTS.
!***********************************************************************
!*  NTIDAL IS THE NUMBER OF MAIN CONSTITUENTS
!*  NTOTAL IS THE NUMBER OF CONSTITUENTS (MAIN + SHALLOW WATER)
!*  FOR  THE GIVEN TIME KH, THE TABLE OF F AND V+U VALUES IS
!*  CALCULATED FOR ALL THE CONSTITUENTS.
!*     F IS THE NODAL MODULATION ADJUSTMENT FACTOR FOR AMPLITUDE
!*     U IS THE NODAL MODULATION ADJUSTMENT FACTOR FOR PHASE
!*     V IS THE ASTRONOMICAL ARGUMENT ADJUSTMENT FOR PHASE.
!
!
!***********************************************************************
!  The astromical arguments are computed at ttime
!
!     d1 is days measured from 1200 UT December 31, 1899
!
      d1=ttime*sec2day
      twopi = 2.0_r8*pi
      call astr(d1,h,pp,s,p,enp,dh,dpp,ds,dp2,dnp)
      intdays=int(ttime*sec2day)
      hh=real(ttime-intdays*day2sec,r8)/3600._r8
      tau = hh/24._r8 + h - s
      dtau_2 = 365.0_r8 + dh - ds
!
!***********************************************************************
!*  ONLY THE FRACTIONAL PART OF A SOLAR DAY NEED BE RETAINED FOR COMPU-
!*  TING THE LUNAR TIME TAU.
!
      DO itide=1,MTC
        freq(itide)=(doodson(1,itide)*dtau_2+doodson(2,itide)*ds+       &
     &               doodson(3,itide)*dh+doodson(4,itide)*dp2+          &
     &               doodson(5,itide)*dnp+doodson(6,itide)*dpp)         &
     &               /(365._r8*24._r8)
        vdbl=doodson(1,itide)*tau+doodson(2,itide)*s+                   &
     &               doodson(3,itide)*h+doodson(4,itide)*p+             &
     &               doodson(5,itide)*enp+doodson(6,itide)*pp+          &
     &               phase_corr(itide)
        iv=vdbl
        iv=(iv/2)*2
        DO j=Jmin,Jmax
          DO i=Imin,Imax
            v=vdbl-iv
            sumc=1.
            sums=0.
            slat=SIN(deg2rad*xlat(i,j))
            DO k=1,num_sat(itide)
!
!***********************************************************************
!*  HERE THE SATELLITE AMPLITUDE RATIO ADJUSTMENT FOR LATITUDE IS MADE
!
              rr=sat_amp(k,itide)
              IF (sat_flag(k,itide) == 1) THEN
                rr=sat_amp(k,itide)*0.36309*(1.-5.*slat*slat)/(eps+slat)
              ELSE IF (sat_flag(k,itide) == 2) THEN
                rr=sat_amp(k,itide)*2.59808*slat
              END IF
              uudbl=sat_doodson(1,k,itide)*p+sat_doodson(2,k,itide)*enp+&
     &              sat_doodson(3,k,itide)*pp+sat_phase_corr(k,itide)
              iuu=uudbl
              uu=uudbl-iuu
              sumc=sumc+rr*COS(uu*twopi)
              sums=sums+rr*SIN(uu*twopi)
            END DO
            fx(i,j,itide)=SQRT(sumc*sumc+sums*sums)
            vux(i,j,itide)=v+ATAN2(sums,sumc)/twopi
            v=v-int(v)
            if (v.lt.0.) v=v+1.
            vux(i,j,itide)=vux(i,j,itide)-int(vux(i,j,itide))
            if (vux(i,j,itide).lt.0.) vux(i,j,itide)=vux(i,j,itide)+1.
!       write(6,998) kon(k),v*360.,f(k),360.*(vux(k)-v)
!       write(7,998) kon(k),v*360.,f(k),360.*(vux(k)-v)
!998    format(' ',a5,5x,3f12.5)
          END DO
        END DO
      END DO
!
      RETURN
      END SUBROUTINE tide_astro
      SUBROUTINE astr(d1,h,pp,s,p,np,dh,dpp,ds,dp2,dnp)
!        this subroutine calculates the following five ephermides
!        of the sun and moon
!        h = mean longitude of the sum
!        pp = mean longitude of the solar perigee
!        s = mean longitude of the moon
!        p = mean longitude of the lunar perigee
!        np = negative of the longitude of the mean ascending node
!        and their rates of change.
!        Units for the ephermides are cycles and for their derivatives
!        are cycles/365 days
!        The formulae for calculating this ephermides were taken from
!        pages 98 and 107 of the Explanatory Supplement to the
!        Astronomical Ephermeris and the American Ephermis and
!        Nautical Almanac (1961)
!
        implicit none
      real(r8), intent(in)  :: d1
      real(r8), intent(out) :: h, pp, s, p, np, dh, dpp, ds, dp2, dnp
      real(r8)  :: d2, f, f2

        d2=d1*1.d-4
        f=360.d0
        f2=f/365.d0
        h=279.696678d0+.9856473354d0*d1+.00002267d0*d2*d2
        pp=281.220833d0+.0000470684d0*d1+.0000339d0*d2*d2+              &
     &  .00000007d0*d2**3
        s=270.434164d0+13.1763965268d0*d1-.000085d0*d2*d2+              &
     &  .000000039d0*d2**3
        p=334.329556d0+.1114040803d0*d1-.0007739d0*d2*d2-               &
     &  .00000026d0*d2**3
        np=-259.183275d0+.0529539222d0*d1-.0001557d0*d2*d2-             &
     &  .00000005d0*d2**3
        h=h/f
        pp=pp/f
        s=s/f
        p=p/f
        np=np/f
        h=h-int(h)
        pp=pp-int(pp)
        s=s-int(s)
        p=p-int(p)
        np=np-int(np)
        dh=.9856473354d0+2.d-8*.00002267d0*d1
        dpp=.0000470684d0+2.d-8*.0000339d0*d1                           &
     &  +3.d-12*.00000007d0*d1**2
        ds=13.1763965268d0-2.d-8*.000085d0*d1+                          &
     &  3.d-12*.000000039d0*d1**2
        dp2=.1114040803d0-2.d-8*.0007739d0*d1-                          &
     &  3.d-12*.00000026d0*d1**2
        dnp=+.0529539222d0-2.d-8*.0001557d0*d1-                         &
     &  3.d-12*.00000005d0*d1**2
        dh=dh/f2
        dpp=dpp/f2
        ds=ds/f2
        dp2=dp2/f2
        dnp=dnp/f2
        RETURN
      END SUBROUTINE astr
#endif
      END MODULE mod_tides
