source: NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE/tide_mod.F90 @ 11864

MODULE tide_mod
   !!======================================================================
   !!                       ***  MODULE  tide_mod  ***
   !! Compute nodal modulations corrections and pulsations
   !!======================================================================
   !! History :  1.0  !  2007  (O. Le Galloudec)  Original code
   !!----------------------------------------------------------------------
   USE oce, ONLY : sshn       ! sea-surface height
   USE par_oce                ! ocean parameters
   USE phycst, ONLY : rpi     ! pi
   USE daymod, ONLY : ndt05   ! half-length of time step
   USE in_out_manager         ! I/O units
   USE iom                    ! xIOs server

   IMPLICIT NONE
   PRIVATE

   PUBLIC   tide_init
   PUBLIC   tide_update          ! called by stp
   PUBLIC   tide_init_harmonics  ! called internally and by module diaharm
   PUBLIC   upd_tide         ! called in dynspg_... modules

   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 64   !: maximum number of harmonic components

   TYPE ::    tide
      CHARACTER(LEN=4) ::   cname_tide = ''
      REAL(wp)         ::   equitide
      INTEGER          ::   nt, ns, nh, np, np1, shift
      INTEGER          ::   nksi, nnu0, nnu1, nnu2, R
      INTEGER          ::   nformula
   END TYPE tide

   TYPE(tide), DIMENSION(:), POINTER ::   tide_components !: Array of selected tidal component parameters

   TYPE, PUBLIC ::   tide_harmonic       !:   Oscillation parameters of harmonic tidal components
      CHARACTER(LEN=4) ::   cname_tide   !    Name of component
      REAL(wp)         ::   equitide     !    Amplitude of equilibrium tide
      REAL(wp)         ::   f            !    Node factor
      REAL(wp)         ::   omega        !    Angular velocity
      REAL(wp)         ::   v0           !    Initial phase at prime meridian
      REAL(wp)         ::   u            !    Phase correction
   END type tide_harmonic

   TYPE(tide_harmonic), PUBLIC, DIMENSION(:), POINTER ::   tide_harmonics !: Oscillation parameters of selected tidal components

   LOGICAL , PUBLIC ::   ln_tide         !:
   LOGICAL , PUBLIC ::   ln_tide_pot     !:
   INTEGER          ::   nn_tide_var     !  Variant of tidal parameter set and tide-potential computation
   LOGICAL          ::   ln_tide_dia     !  Enable tidal diagnostic output
   LOGICAL          ::   ln_read_load    !:
   LOGICAL , PUBLIC ::   ln_scal_load    !:
   LOGICAL , PUBLIC ::   ln_tide_ramp    !:
   INTEGER , PUBLIC ::   nb_harmo        !: Number of active tidal components
   REAL(wp), PUBLIC ::   rn_tide_ramp_dt     !:
   REAL(wp), PUBLIC ::   rn_scal_load    !:
   CHARACTER(lc), PUBLIC ::   cn_tide_load   !: 
   REAL(wp)         ::   rn_tide_gamma   ! Tidal tilt factor

   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   pot_astro !: tidal potential
   REAL(wp),         ALLOCATABLE, DIMENSION(:,:)   ::   pot_astro_comp   !  tidal-potential component
   REAL(wp),         ALLOCATABLE, DIMENSION(:,:,:) ::   amp_pot, phi_pot
   REAL(wp),         ALLOCATABLE, DIMENSION(:,:,:) ::   amp_load, phi_load

   REAL(wp) ::   rn_tide_ramp_t !   Elapsed time in seconds

   REAL(wp) ::   sh_T, sh_s, sh_h, sh_p, sh_p1             ! astronomic angles
   REAL(wp) ::   sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R   !
   REAL(wp) ::   sh_I, sh_x1ra, sh_N                       !

   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id$ 
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE tide_init
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE tide_init  ***
      !!----------------------------------------------------------------------      
      INTEGER  :: ji, jk
      CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: sn_tide_cnames ! Names of selected tidal components
      INTEGER  ::   ios                 ! Local integer output status for namelist read
      !
      NAMELIST/nam_tide/ln_tide, nn_tide_var, ln_tide_dia, ln_tide_pot, rn_tide_gamma, &
         &              ln_scal_load, ln_read_load, cn_tide_load,         &
         &              ln_tide_ramp, rn_scal_load, rn_tide_ramp_dt,      &
         &              sn_tide_cnames
      !!----------------------------------------------------------------------
      !
      ! Initialise all array elements of sn_tide_cnames, as some of them
      ! typically do not appear in namelist_ref or namelist_cfg
      sn_tide_cnames(:) = ''
      ! Read Namelist nam_tide
      REWIND( numnam_ref )              ! Namelist nam_tide in reference namelist : Tides
      READ  ( numnam_ref, nam_tide, IOSTAT = ios, ERR = 901)
901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_tide in reference namelist', lwp )
      !
      REWIND( numnam_cfg )              ! Namelist nam_tide in configuration namelist : Tides
      READ  ( numnam_cfg, nam_tide, IOSTAT = ios, ERR = 902 )
902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_tide in configuration namelist', lwp )
      IF(lwm) WRITE ( numond, nam_tide )
      !
      IF( ln_tide ) THEN
         IF (lwp) THEN
            WRITE(numout,*)
            WRITE(numout,*) 'tide_init : Initialization of the tidal components'
            WRITE(numout,*) '~~~~~~~~~ '
            WRITE(numout,*) '   Namelist nam_tide'
            WRITE(numout,*) '      Use tidal components                       ln_tide         = ', ln_tide
            WRITE(numout,*) '         Variant (1: default; 0: legacy option)  nn_tide_var     = ', nn_tide_var
            WRITE(numout,*) '         Tidal diagnostic output                 ln_tide_dia     = ', ln_tide_dia
            WRITE(numout,*) '         Apply astronomical potential            ln_tide_pot     = ', ln_tide_pot
            WRITE(numout,*) '         Tidal tilt factor                       rn_tide_gamma   = ', rn_tide_gamma
            WRITE(numout,*) '         Use scalar approx. for load potential   ln_scal_load    = ', ln_scal_load
            WRITE(numout,*) '         Read load potential from file           ln_read_load    = ', ln_read_load
            WRITE(numout,*) '         Apply ramp on tides at startup          ln_tide_ramp    = ', ln_tide_ramp
            WRITE(numout,*) '         Fraction of SSH used in scal. approx.   rn_scal_load    = ', rn_scal_load
            WRITE(numout,*) '         Duration (days) of ramp                 rn_tide_ramp_dt = ', rn_tide_ramp_dt
         ENDIF
      ELSE
         rn_scal_load = 0._wp 

         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'tide_init : tidal components not used (ln_tide = F)'
         IF(lwp) WRITE(numout,*) '~~~~~~~~~ '
         RETURN
      ENDIF
      !
      IF( ln_read_load.AND.(.NOT.ln_tide_pot) ) &
          &   CALL ctl_stop('ln_read_load requires ln_tide_pot')
      IF( ln_scal_load.AND.(.NOT.ln_tide_pot) ) &
          &   CALL ctl_stop('ln_scal_load requires ln_tide_pot')
      IF( ln_scal_load.AND.ln_read_load ) &
          &   CALL ctl_stop('Choose between ln_scal_load and ln_read_load')
      IF( ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rn_tide_ramp_dt) )   &
         &   CALL ctl_stop('rn_tide_ramp_dt must be lower than run duration')
      IF( ln_tide_ramp.AND.(rn_tide_ramp_dt<0.) ) &
         &   CALL ctl_stop('rn_tide_ramp_dt must be positive')
      !
      ! Initialise array used to store tidal oscillation parameters (frequency,
      ! amplitude, phase); also retrieve and store array of information about
      ! selected tidal components
      CALL tide_init_harmonics(sn_tide_cnames, tide_harmonics, tide_components)
      !
      ! Number of active tidal components
      nb_harmo = size(tide_components)
      !       
      ! Ensure that tidal components have been set in namelist_cfg
      IF( nb_harmo == 0 )   CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' )
      !
      IF (.NOT.ln_scal_load ) rn_scal_load = 0._wp
      !
      ALLOCATE( amp_pot(jpi,jpj,nb_harmo),                      &
           &      phi_pot(jpi,jpj,nb_harmo), pot_astro(jpi,jpj)   )
      IF( ln_tide_dia ) ALLOCATE( pot_astro_comp(jpi,jpj) )
      IF( ln_read_load ) THEN
         ALLOCATE( amp_load(jpi,jpj,nb_harmo), phi_load(jpi,jpj,nb_harmo) )
         CALL tide_init_load
         amp_pot(:,:,:) = amp_load(:,:,:)
         phi_pot(:,:,:) = phi_load(:,:,:)
      ELSE     
         amp_pot(:,:,:) = 0._wp
         phi_pot(:,:,:) = 0._wp   
      ENDIF
      !
   END SUBROUTINE tide_init


   SUBROUTINE tide_init_components(pcnames, ptide_comp)
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE tide_init_components  ***
      !!
      !! Returns pointer to array of variables of type 'tide' that contain
      !! information about the selected tidal components
      !! ----------------------------------------------------------------------
      CHARACTER(LEN=4),              DIMENSION(jpmax_harmo), INTENT(in)  ::   pcnames             ! Names of selected components
      TYPE(tide),       POINTER,     DIMENSION(:),           INTENT(out) ::   ptide_comp          ! Selected components
      INTEGER,          ALLOCATABLE, DIMENSION(:)                        ::   icomppos            ! Indices of selected components
      INTEGER                                                            ::   icomp, jk, jj, ji   ! Miscellaneous integers
      LOGICAL                                                            ::   llmatch             ! Local variables used for
      INTEGER                                                            ::   ic1, ic2            !    string comparison
      TYPE(tide),       POINTER,     DIMENSION(:)                        ::   tide_components     ! All available components
      
      ! Populate local array with information about all available tidal
      ! components
      !
      ! Note, here 'tide_components' locally overrides the global module
      ! variable of the same name to enable the use of the global name in the
      ! include file that contains the initialisation of elements of array
      ! 'tide_components'
      ALLOCATE(tide_components(jpmax_harmo), icomppos(jpmax_harmo))
      ! Initialise array of indices of the selected componenents
      icomppos(:) = 0
      ! Include tidal component parameters for all available components
      IF (nn_tide_var < 1) THEN
#define TIDE_VAR_0
#include "tide.h90"
#undef TIDE_VAR_0
      ELSE
#include "tide.h90"
      END IF
      ! Identify the selected components that are availble
      icomp = 0
      DO jk = 1, jpmax_harmo
         IF (TRIM(pcnames(jk)) /= '') THEN
            DO jj = 1, jpmax_harmo
               ! Find matches between selected and available constituents
               ! (ignore capitalisation unless legacy variant has been selected)
               IF (nn_tide_var < 1) THEN
                  llmatch = (TRIM(pcnames(jk)) == TRIM(tide_components(jj)%cname_tide))
               ELSE
                  llmatch = .TRUE.
                  ji = MAX(LEN_TRIM(pcnames(jk)), LEN_TRIM(tide_components(jj)%cname_tide))
                  DO WHILE (llmatch.AND.(ji > 0))
                     ic1 = IACHAR(pcnames(jk)(ji:ji))
                     IF ((ic1 >= 97).AND.(ic1 <= 122)) ic1 = ic1 - 32
                     ic2 = IACHAR(tide_components(jj)%cname_tide(ji:ji))
                     IF ((ic2 >= 97).AND.(ic2 <= 122)) ic2 = ic2 - 32
                     llmatch = (ic1 == ic2)
                     ji = ji - 1
                  END DO
               END IF
               IF (llmatch) THEN
                  ! Count and record the match
                  icomp = icomp + 1
                  icomppos(icomp) = jj
                  ! Set the capitalisation of the tidal constituent identifier
                  ! as specified in the namelist
                  tide_components(jj)%cname_tide = pcnames(jk)
                  IF (lwp) WRITE(numout, '(10X,"Tidal component #",I2.2,36X,"= ",A4)') icomp, tide_components(jj)%cname_tide
                  EXIT
               END IF
            END DO
            IF ((lwp).AND.(jj > jpmax_harmo)) WRITE(numout, '(10X,"Tidal component ",A4," is not available!")') pcnames(jk)
         END IF
      END DO
      
      ! Allocate and populate reduced list of components
      ALLOCATE(ptide_comp(icomp))
      DO jk = 1, icomp
         ptide_comp(jk) = tide_components(icomppos(jk))
      END DO
      
      ! Release local array of available components and list of selected
      ! components
      DEALLOCATE(tide_components, icomppos)
      
   END SUBROUTINE tide_init_components


   SUBROUTINE tide_init_harmonics(pcnames, ptide_harmo, ptide_comp)
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE tide_init_harmonics  ***
      !!
      !! Returns pointer to array of variables of type 'tide_harmonics' that
      !! contain oscillation parameters of the selected harmonic tidal
      !! components
      !! ----------------------------------------------------------------------
      CHARACTER(LEN=4),             DIMENSION(jpmax_harmo), INTENT(in) ::   pcnames     ! Names of selected components
      TYPE(tide_harmonic), POINTER, DIMENSION(:)                       ::   ptide_harmo ! Oscillation parameters of tidal components
      TYPE(tide),          POINTER, DIMENSION(:), OPTIONAL             ::   ptide_comp  ! Selected components
      TYPE(tide),          POINTER, DIMENSION(:)                       ::   ztcomp      ! Selected components

      ! Retrieve information about selected tidal components
      ! If requested, prepare tidal component array for returning
      IF (PRESENT(ptide_comp)) THEN
         CALL tide_init_components(pcnames, ptide_comp)
         ztcomp => ptide_comp
      ELSE
         CALL tide_init_components(pcnames, ztcomp)
      END IF

      ! Allocate and populate array of oscillation parameters
      ALLOCATE(ptide_harmo(size(ztcomp)))
      ptide_harmo(:)%cname_tide = ztcomp(:)%cname_tide
      ptide_harmo(:)%equitide = ztcomp(:)%equitide
      CALL tide_harmo(ztcomp, ptide_harmo)

   END SUBROUTINE tide_init_harmonics


   SUBROUTINE tide_init_potential
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE tide_init_potential  ***
      !!
      !! *** Reference:
      !!      CT71) Cartwright, D. E. and Tayler, R. J. (1971): New computations of
      !!            the Tide-generating Potential. Geophys. J. R. astr. Soc. 23,
      !!            pp. 45-74. DOI: 10.1111/j.1365-246X.1971.tb01803.x
      !!
      !!----------------------------------------------------------------------
      INTEGER  ::   ji, jj, jk   ! dummy loop indices
      REAL(wp) ::   zcons, ztmp1, ztmp2, zlat, zlon, ztmp, zamp, zcs   ! local scalar
      !!----------------------------------------------------------------------

      IF( ln_read_load ) THEN
    amp_pot(:,:,:) = amp_load(:,:,:)
    phi_pot(:,:,:) = phi_load(:,:,:)
      ELSE     
         amp_pot(:,:,:) = 0._wp
         phi_pot(:,:,:) = 0._wp   
      ENDIF
      DO jk = 1, nb_harmo
         zcons = rn_tide_gamma * tide_components(jk)%equitide * tide_harmonics(jk)%f
         DO ji = 1, jpi
            DO jj = 1, jpj
               ztmp1 =  tide_harmonics(jk)%f * amp_pot(ji,jj,jk) * COS( phi_pot(ji,jj,jk) &
                  &                                                   + tide_harmonics(jk)%v0 + tide_harmonics(jk)%u )
               ztmp2 = -tide_harmonics(jk)%f * amp_pot(ji,jj,jk) * SIN( phi_pot(ji,jj,jk) &
                  &                                                   + tide_harmonics(jk)%v0 + tide_harmonics(jk)%u )
               zlat = gphit(ji,jj)*rad !! latitude en radian
               zlon = glamt(ji,jj)*rad !! longitude en radian
               ztmp = tide_harmonics(jk)%v0 + tide_harmonics(jk)%u + tide_components(jk)%nt * zlon
               ! le potentiel est composé des effets des astres:
               SELECT CASE( tide_components(jk)%nt )
               CASE( 0 )                                                  !   long-periodic tidal constituents (included unless
                  zcs = zcons * ( 0.5_wp - 1.5_wp * SIN( zlat )**2 )      !   compatibility with original formulation is requested)
                  IF ( nn_tide_var < 1 ) zcs = 0.0_wp
               CASE( 1 )                                                  !   diurnal tidal constituents
                  zcs = zcons * SIN( 2.0_wp*zlat )
               CASE( 2 )                                                  !   semi-diurnal tidal constituents
                  zcs = zcons * COS( zlat )**2
               CASE( 3 )                                                  !   Terdiurnal tidal constituents; the colatitude-dependent
                  zcs = zcons * COS( zlat )**3                            !   factor is sin(theta)^3 (Table 2 of CT71)
               CASE DEFAULT                                               !   constituents of higher frequency are not included
                  zcs = 0.0_wp
               END SELECT
               ztmp1 = ztmp1 + zcs * COS( ztmp )
               ztmp2 = ztmp2 - zcs * SIN( ztmp )
               zamp = SQRT( ztmp1*ztmp1 + ztmp2*ztmp2 )
               amp_pot(ji,jj,jk) = zamp
               phi_pot(ji,jj,jk) = ATAN2( -ztmp2 / MAX( 1.e-10_wp , zamp ) ,   &
                  &                        ztmp1 / MAX( 1.e-10_wp,  zamp )   )
            END DO
         END DO
      END DO
      !
   END SUBROUTINE tide_init_potential


   SUBROUTINE tide_init_load
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE tide_init_load  ***
      !!----------------------------------------------------------------------
      INTEGER :: inum                 ! Logical unit of input file
      INTEGER :: ji, jj, itide        ! dummy loop indices
      REAL(wp), DIMENSION(jpi,jpj) ::   ztr, zti   !: workspace to read in tidal harmonics data 
      !!----------------------------------------------------------------------
      IF(lwp) THEN
         WRITE(numout,*)
         WRITE(numout,*) 'tide_init_load : Initialization of load potential from file'
         WRITE(numout,*) '~~~~~~~~~~~~~~ '
      ENDIF
      !
      CALL iom_open ( cn_tide_load , inum )
      !
      DO itide = 1, nb_harmo
         CALL iom_get  ( inum, jpdom_data,TRIM(tide_components(itide)%cname_tide)//'_z1', ztr(:,:) )
         CALL iom_get  ( inum, jpdom_data,TRIM(tide_components(itide)%cname_tide)//'_z2', zti(:,:) )
         !
         DO ji=1,jpi
            DO jj=1,jpj
               amp_load(ji,jj,itide) =  SQRT( ztr(ji,jj)**2. + zti(ji,jj)**2. )
               phi_load(ji,jj,itide) = ATAN2(-zti(ji,jj), ztr(ji,jj) )
            END DO
         END DO
         !
      END DO
      CALL iom_close( inum )
      !
   END SUBROUTINE tide_init_load


   SUBROUTINE tide_update( kt )
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE tide_update  ***
      !!----------------------------------------------------------------------      
      INTEGER, INTENT( in ) ::   kt     ! ocean time-step
      INTEGER               ::   jk     ! dummy loop index
      !!----------------------------------------------------------------------
      
      IF( nsec_day == NINT(0.5_wp * rdt) .OR. kt == nit000 ) THEN      ! start a new day
         !
         CALL tide_harmo(tide_components, tide_harmonics, ndt05) ! Update oscillation parameters of tidal components for start of current day
         !
         !
         IF(lwp) THEN
            WRITE(numout,*)
            WRITE(numout,*) 'tide_update : Update of the components and (re)Init. the potential at kt=', kt
            WRITE(numout,*) '~~~~~~~~~~~ '
            DO jk = 1, nb_harmo
               WRITE(numout,*) tide_harmonics(jk)%cname_tide, tide_harmonics(jk)%u, &
                  &            tide_harmonics(jk)%f,tide_harmonics(jk)%v0, tide_harmonics(jk)%omega
            END DO
         ENDIF
         !
         IF( ln_tide_pot )   CALL tide_init_potential
         !
         rn_tide_ramp_t = (kt - nit000)*rdt !   Elapsed time in seconds
      ENDIF
      !
   END SUBROUTINE tide_update


   SUBROUTINE tide_harmo( ptide_comp, ptide_harmo, psec_day )
      !
      TYPE(tide),          DIMENSION(:), POINTER ::   ptide_comp   ! Array of selected tidal component parameters
      TYPE(tide_harmonic), DIMENSION(:), POINTER ::   ptide_harmo  ! Oscillation parameters of selected tidal components
      INTEGER, OPTIONAL ::   psec_day                              ! Number of seconds since the start of the current day
      !
      IF (PRESENT(psec_day)) THEN 
         CALL astronomic_angle(psec_day)
      ELSE
         CALL astronomic_angle(nsec_day)
      END IF
      CALL tide_pulse( ptide_comp, ptide_harmo )
      CALL tide_vuf(   ptide_comp, ptide_harmo )
      !
   END SUBROUTINE tide_harmo


   SUBROUTINE astronomic_angle(psec_day)
      !!----------------------------------------------------------------------
      !!  tj is time elapsed since 1st January 1900, 0 hour, counted in julian
      !!  century (e.g. time in days divide by 36525)
      !!----------------------------------------------------------------------
      INTEGER  ::   psec_day !   Number of seconds from midnight
      REAL(wp) ::   cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2
      REAL(wp) ::   zqy , zsy, zday, zdj, zhfrac
      !!----------------------------------------------------------------------
      !
      zqy = AINT( (nyear-1901.)/4. )
      zsy = nyear - 1900.
      !
      zdj  = dayjul( nyear, nmonth, nday )
      zday = zdj + zqy - 1.
      !
      zhfrac = psec_day / 3600.
      !
      !----------------------------------------------------------------------
      !  Sh_n Longitude of ascending lunar node
      !----------------------------------------------------------------------
      sh_N=(259.1560564-19.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*rad
      !----------------------------------------------------------------------
      ! T mean solar angle (Greenwhich time)
      !----------------------------------------------------------------------
      sh_T=(180.+zhfrac*(360./24.))*rad
      !----------------------------------------------------------------------
      ! h mean solar Longitude
      !----------------------------------------------------------------------
      sh_h=(280.1895014-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad
      !----------------------------------------------------------------------
      ! s mean lunar Longitude
      !----------------------------------------------------------------------
      sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad
      !----------------------------------------------------------------------
      ! p1 Longitude of solar perigee
      !----------------------------------------------------------------------
      sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad
      !----------------------------------------------------------------------
      ! p Longitude of lunar perigee
      !----------------------------------------------------------------------
      sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad

      sh_N = MOD( sh_N ,2*rpi )
      sh_s = MOD( sh_s ,2*rpi )
      sh_h = MOD( sh_h, 2*rpi )
      sh_p = MOD( sh_p, 2*rpi )
      sh_p1= MOD( sh_p1,2*rpi )

      cosI = 0.913694997 -0.035692561 *cos(sh_N)

      sh_I = ACOS( cosI )

      sin2I   = sin(sh_I)
      sh_tgn2 = tan(sh_N/2.0)

      at1=atan(1.01883*sh_tgn2)
      at2=atan(0.64412*sh_tgn2)

      sh_xi=-at1-at2+sh_N

      IF( sh_N > rpi )   sh_xi=sh_xi-2.0*rpi

      sh_nu = at1 - at2

      !----------------------------------------------------------------------
      ! For constituents l2 k1 k2
      !----------------------------------------------------------------------

      tgI2 = tan(sh_I/2.0)
      P1   = sh_p-sh_xi

      t2 = tgI2*tgI2
      t4 = t2*t2
      sh_x1ra = sqrt( 1.0-12.0*t2*cos(2.0*P1)+36.0*t4 )

      p = sin(2.0*P1)
      q = 1.0/(6.0*t2)-cos(2.0*P1)
      sh_R = atan(p/q)

      p = sin(2.0*sh_I)*sin(sh_nu)
      q = sin(2.0*sh_I)*cos(sh_nu)+0.3347
      sh_nuprim = atan(p/q)

      s2 = sin(sh_I)*sin(sh_I)
      p  = s2*sin(2.0*sh_nu)
      q  = s2*cos(2.0*sh_nu)+0.0727
      sh_nusec = 0.5*atan(p/q)
      !
   END SUBROUTINE astronomic_angle


   SUBROUTINE tide_pulse( ptide_comp, ptide_harmo )
      !!----------------------------------------------------------------------
      !!                     ***  ROUTINE tide_pulse  ***
      !!                      
      !! ** Purpose : Compute tidal frequencies
      !!----------------------------------------------------------------------
      TYPE(tide),          DIMENSION(:), POINTER ::   ptide_comp   ! Array of selected tidal component parameters
      TYPE(tide_harmonic), DIMENSION(:), POINTER ::   ptide_harmo  ! Oscillation parameters of selected tidal components
      !
      INTEGER  ::   jh
      REAL(wp) ::   zscale
      REAL(wp) ::   zomega_T =  13149000.0_wp
      REAL(wp) ::   zomega_s =    481267.892_wp
      REAL(wp) ::   zomega_h =     36000.76892_wp
      REAL(wp) ::   zomega_p =      4069.0322056_wp
      REAL(wp) ::   zomega_n =      1934.1423972_wp
      REAL(wp) ::   zomega_p1=         1.719175_wp
      !!----------------------------------------------------------------------
      !
      zscale =  rad / ( 36525._wp * 86400._wp ) 
      !
      DO jh = 1, size(ptide_harmo)
         ptide_harmo(jh)%omega = (  zomega_T * ptide_comp( jh )%nT   &
            &                         + zomega_s * ptide_comp( jh )%ns   &
            &                         + zomega_h * ptide_comp( jh )%nh   &
            &                         + zomega_p * ptide_comp( jh )%np   &
            &                         + zomega_p1* ptide_comp( jh )%np1  ) * zscale
      END DO
      !
   END SUBROUTINE tide_pulse


   SUBROUTINE tide_vuf( ptide_comp, ptide_harmo )
      !!----------------------------------------------------------------------
      !!                     ***  ROUTINE tide_vuf  ***
      !!                      
      !! ** Purpose : Compute nodal modulation corrections
      !!
      !! ** Outputs : vt: Phase of tidal potential relative to Greenwich (radians)
      !!              ut: Phase correction u due to nodal motion (radians)
      !!              ft: Nodal correction factor
      !!----------------------------------------------------------------------
      TYPE(tide),          DIMENSION(:), POINTER ::   ptide_comp   ! Array of selected tidal component parameters
      TYPE(tide_harmonic), DIMENSION(:), POINTER ::   ptide_harmo  ! Oscillation parameters of selected tidal components
      !
      INTEGER ::   jh   ! dummy loop index
      !!----------------------------------------------------------------------
      !
      DO jh = 1, size(ptide_harmo)
         !  Phase of the tidal potential relative to the Greenwhich 
         !  meridian (e.g. the position of the fictuous celestial body). Units are radian:
         ptide_harmo(jh)%v0 = sh_T * ptide_comp( jh )%nT    &
            &                   + sh_s * ptide_comp( jh )%ns    &
            &                   + sh_h * ptide_comp( jh )%nh    &
            &                   + sh_p * ptide_comp( jh )%np    &
            &                   + sh_p1* ptide_comp( jh )%np1   &
            &                   +        ptide_comp( jh )%shift * rad
         !
         !  Phase correction u due to nodal motion. Units are radian:
         ptide_harmo(jh)%u = sh_xi     * ptide_comp( jh )%nksi   &
            &                  + sh_nu     * ptide_comp( jh )%nnu0   &
            &                  + sh_nuprim * ptide_comp( jh )%nnu1   &
            &                  + sh_nusec  * ptide_comp( jh )%nnu2   &
            &                  + sh_R      * ptide_comp( jh )%R

         !  Nodal correction factor:
         ptide_harmo(jh)%f = nodal_factort( ptide_comp( jh )%nformula )
      END DO
      !
   END SUBROUTINE tide_vuf


   RECURSIVE FUNCTION nodal_factort( kformula ) RESULT( zf )
      !!----------------------------------------------------------------------
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) :: kformula
      !
      REAL(wp)            :: zf
      REAL(wp)            :: zs, zf1, zf2
      CHARACTER(LEN=3)    :: clformula
      !!----------------------------------------------------------------------
      !
      SELECT CASE( kformula )
      !
      CASE( 0 )                  !==  formule 0, solar waves
         zf = 1.0
         !
      CASE( 1 )                  !==  formule 1, compound waves (78 x 78)
         zf=nodal_factort(78)
         zf = zf * zf
         !
      CASE ( 4 )                 !==  formule 4,  compound waves (78 x 235) 
         zf1 = nodal_factort( 78)
         zf  = nodal_factort(235)
         zf  = zf1 * zf
         !
      CASE( 18 )                 !==  formule 18,  compound waves (78 x 78 x 78 )
         zf1 = nodal_factort(78)
         zf  = zf1 * zf1 * zf1
         !
      CASE( 20 )                 !==  formula 20, compound waves ( 78 x 78 x 78 x 78 )
         zf1 = nodal_factort(78)
         zf  = zf1 * zf1 * zf1 * zf1
         !
      CASE( 73 )                 !==  formule 73
         zs = sin(sh_I)
         zf = (2./3.-zs*zs)/0.5021
         !
      CASE( 74 )                 !==  formule 74
         zs = sin(sh_I)
         zf = zs * zs / 0.1578
         !
      CASE( 75 )                 !==  formule 75
         zs = cos(sh_I/2)
         zf = sin(sh_I) * zs * zs / 0.3800
         !
      CASE( 76 )                 !==  formule 76
         zf = sin(2*sh_I) / 0.7214
         !
      CASE( 78 )                 !==  formule 78
         zs = cos(sh_I/2)
         zf = zs * zs * zs * zs / 0.9154
         !
      CASE( 149 )                !==  formule 149
         zs = cos(sh_I/2)
         zf = zs*zs*zs*zs*zs*zs / 0.8758
         !
      CASE( 215 )                !==  formule 215
         zs = cos(sh_I/2)
         zf = zs*zs*zs*zs / 0.9154 * sh_x1ra
         !
      CASE( 227 )                !==  formule 227 
         zs = sin(2*sh_I)
         zf = sqrt( 0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006 )
         !
      CASE ( 235 )               !==  formule 235 
         zs = sin(sh_I)
         zf = sqrt( 19.0444*zs*zs*zs*zs + 2.7702*zs*zs*cos(2*sh_nu) + .0981 )
         !
      CASE DEFAULT
         WRITE( clformula, '(I3)' ) kformula
         CALL ctl_stop('nodal_factort: formula ' // clformula // ' is not available')
      END SELECT
      !
   END FUNCTION nodal_factort


   FUNCTION dayjul( kyr, kmonth, kday )
      !!----------------------------------------------------------------------
      !!  *** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE)
      !!----------------------------------------------------------------------
      INTEGER,INTENT(in) ::   kyr, kmonth, kday
      !
      INTEGER,DIMENSION(12) ::  idayt, idays
      INTEGER  ::   inc, ji
      REAL(wp) ::   dayjul, zyq
      !
      DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./
      !!----------------------------------------------------------------------
      !
      idays(1) = 0.
      idays(2) = 31.
      inc = 0.
      zyq = MOD( kyr-1900. , 4. )
      IF( zyq == 0.)   inc = 1.
      DO ji = 3, 12
         idays(ji)=idayt(ji)+inc
      END DO
      dayjul = idays(kmonth) + kday
      !
   END FUNCTION dayjul


   SUBROUTINE upd_tide(pdelta)
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE upd_tide  ***
      !!
      !! ** Purpose :   provide at each time step the astronomical potential
      !!
      !! ** Method  :   computed from pulsation and amplitude of all tide components
      !!
      !! ** Action  :   pot_astro   actronomical potential
      !!----------------------------------------------------------------------      
      REAL, INTENT(in)              ::   pdelta      ! Temporal offset in seconds
      INTEGER                       ::   jk          ! Dummy loop index
      REAL(wp)                      ::   zt, zramp   ! Local scalars
      REAL(wp), DIMENSION(nb_harmo) ::   zwt         ! Temporary array
      !!----------------------------------------------------------------------      
      !
      zwt(:) = tide_harmonics(:)%omega * pdelta
      !
      IF( ln_tide_ramp ) THEN         ! linear increase if asked
         zt = rn_tide_ramp_t + pdelta
         zramp = MIN(  MAX( zt / (rn_tide_ramp_dt*rday) , 0._wp ) , 1._wp  )
      ENDIF
      !
      pot_astro(:,:) = 0._wp          ! update tidal potential (sum of all harmonics)
      DO jk = 1, nb_harmo
         IF ( .NOT. ln_tide_dia ) THEN
            pot_astro(:,:) = pot_astro(:,:) + amp_pot(:,:,jk) * COS( zwt(jk) + phi_pot(:,:,jk) )
         ELSE
            pot_astro_comp(:,:) = amp_pot(:,:,jk) * COS( zwt(jk) + phi_pot(:,:,jk) )
            pot_astro(:,:) = pot_astro(:,:) + pot_astro_comp(:,:)
            IF ( iom_use( "tide_pot_" // TRIM( tide_harmonics(jk)%cname_tide ) ) ) THEN   ! Output tidal potential (incl. load potential)
               IF ( ln_tide_ramp ) pot_astro_comp(:,:) = zramp * pot_astro_comp(:,:)
               CALL iom_put( "tide_pot_" // TRIM( tide_harmonics(jk)%cname_tide ), pot_astro_comp(:,:) )
            END IF
         END IF
      END DO
      !
      IF ( ln_tide_ramp ) pot_astro(:,:) = zramp * pot_astro(:,:)
      !
      IF( ln_tide_dia ) THEN          ! Output total tidal potential (incl. load potential)
         IF ( iom_use( "tide_pot" ) ) CALL iom_put( "tide_pot", pot_astro(:,:) + rn_scal_load * sshn(:,:) )
      END IF
      !
   END SUBROUTINE upd_tide

   !!======================================================================
END MODULE tide_mod