Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/cfgs/AMM12/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/cfgs/AMM12/EXPREF/namelist_cfg	(revision 12116)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/cfgs/AMM12/EXPREF/namelist_cfg	(revision 12117)
@@ -154,22 +154,21 @@
 &nam_tide      !   tide parameters                                      (default: OFF)
 !-----------------------------------------------------------------------
-   ln_tide     = .true.    ! Activate tides
-      ln_tide_pot   = .true.  !  use tidal potential forcing
-
-   clname(1)     =   'Q1'   !  name of constituent
-   clname(2)     =   'O1'
-   clname(3)     =   'P1'
-   clname(4)     =   'S1'
-   clname(5)     =   'K1'
-   clname(6)     =   '2N2'
-   clname(7)     =   'MU2'
-   clname(8)     =   'N2'
-   clname(9)     =   'NU2'
-   clname(10)    =   'M2'
-   clname(11)    =   'L2'
-   clname(12)    =   'T2'
-   clname(13)    =   'S2'
-   clname(14)    =   'K2'
-   clname(15)    =   'M4'
+   ln_tide     = .true.        ! Activate tides
+      ln_tide_pot     = .true. !  use tidal potential forcing
+   sn_tide_cnames(1)  = 'Q1'   !  name of constituent
+   sn_tide_cnames(2)  = 'O1'
+   sn_tide_cnames(3)  = 'P1'
+   sn_tide_cnames(4)  = 'S1'
+   sn_tide_cnames(5)  = 'K1'
+   sn_tide_cnames(6)  = '2N2'
+   sn_tide_cnames(7)  = 'MU2'
+   sn_tide_cnames(8)  = 'N2'
+   sn_tide_cnames(9)  = 'NU2'
+   sn_tide_cnames(10) = 'M2'
+   sn_tide_cnames(11) = 'L2'
+   sn_tide_cnames(12) = 'T2'
+   sn_tide_cnames(13) = 'S2'
+   sn_tide_cnames(14) = 'K2'
+   sn_tide_cnames(15) = 'M4'
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12116)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12117)
@@ -206,4 +206,30 @@
 
       </field_group> <!-- grid_T -->
+
+      <!-- Tides -->
+
+      <field_group id="Tides_T" grid_ref="grid_T_2D" operation="once" >
+        <!-- Tidal potential -->
+        <field id="tide_pot"      long_name="Total tidal potential" unit="m" />
+        <field id="tide_pot_M2"   long_name="M2 tidal potential"    unit="m" />
+        <field id="tide_pot_N2"   long_name="N2 tidal potential"    unit="m" />
+        <field id="tide_pot_2N2"  long_name="2N2 tidal potential"   unit="m" />
+        <field id="tide_pot_S2"   long_name="S2 tidal potential"    unit="m" />
+        <field id="tide_pot_K2"   long_name="K2 tidal potential"    unit="m" />
+        <field id="tide_pot_K1"   long_name="K1 tidal potential"    unit="m" />
+        <field id="tide_pot_O1"   long_name="O1 tidal potential"    unit="m" />
+        <field id="tide_pot_Q1"   long_name="Q1 tidal potential"    unit="m" />
+        <field id="tide_pot_P1"   long_name="P1 tidal potential"    unit="m" />
+        <field id="tide_pot_M4"   long_name="M4 tidal potential"    unit="m" />
+        <field id="tide_pot_Mf"   long_name="Mf tidal potential"    unit="m" />
+        <field id="tide_pot_Mm"   long_name="Mm tidal potential"    unit="m" />
+        <field id="tide_pot_Msqm" long_name="Msqm tidal potential"  unit="m" />
+        <field id="tide_pot_Mtm"  long_name="Mtm tidal potential"   unit="m" />
+        <field id="tide_pot_S1"   long_name="S1 tidal potential"    unit="m" />
+        <field id="tide_pot_MU2"  long_name="MU2 tidal potential"   unit="m" />
+        <field id="tide_pot_NU2"  long_name="NU2 tidal potential"   unit="m" />
+        <field id="tide_pot_L2"   long_name="L2 tidal potential"    unit="m" />
+        <field id="tide_pot_T2"   long_name="T2 tidal potential"    unit="m" />
+      </field_group>
 
       <!-- OSMOSIS variables (available with ln_zdfosm=.true.) -->
Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/cfgs/SHARED/namelist_ref	(revision 12116)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/cfgs/SHARED/namelist_ref	(revision 12117)
@@ -568,5 +568,9 @@
 !-----------------------------------------------------------------------
    ln_tide     = .false.      ! Activate tides
+      nn_tide_var   = 1          !  Variant of tidal parameter set and tide-potential computation
+      !                          !     (1: default; 0: compatibility with previous versions)
+      ln_tide_dia   = .false.    !  Enable tidal diagnostic output
       ln_tide_pot   = .true.                !  use tidal potential forcing
+         rn_tide_gamma = 0.7                   ! Tidal tilt factor
          ln_scal_load  = .false.               ! Use scalar approximation for
             rn_scal_load = 0.094               !     load potential
@@ -575,6 +579,6 @@
             !      
       ln_tide_ramp  = .false.               !  Use linear ramp for tides at startup
-         rdttideramp   =    0.                 !  ramp duration in days
-      clname(1)     = 'DUMMY'               !  name of constituent - all tidal components must be set in namelist_cfg
+         rn_tide_ramp_dt = 0.               !  ramp duration in days
+      sn_tide_cnames(1) = 'DUMMY'               !  name of constituent - all tidal components must be set in namelist_cfg
 /
 !-----------------------------------------------------------------------
@@ -652,5 +656,4 @@
    filtide          = 'bdydta/amm12_bdytide_'   !  file name root of tidal forcing files
    ln_bdytide_2ddta = .false.                   !
-   ln_bdytide_conj  = .false.                   ! 
 /
 
Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/cfgs/SPITZ12/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/cfgs/SPITZ12/EXPREF/namelist_cfg	(revision 12116)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/cfgs/SPITZ12/EXPREF/namelist_cfg	(revision 12117)
@@ -139,14 +139,14 @@
 &nam_tide      !   tide parameters                                      (default: OFF)
 !-----------------------------------------------------------------------
-   ln_tide     = .true.      ! Activate tides
-      ln_tide_pot   = .false.                !  use tidal potential forcing
-      clname(1)     =   'M2'    !  name of constituent
-      clname(2)     =   'S2'
-      clname(3)     =   'N2'
-      clname(4)     =   'K2'
-      clname(5)     =   'K1'
-      clname(6)     =   'O1'
-      clname(7)     =   'P1'
-      clname(8)     =   'Q1'
+   ln_tide     = .true.           ! Activate tides
+      ln_tide_pot       = .false. !  use tidal potential forcing
+      sn_tide_cnames(1) =   'M2'  !  name of constituent
+      sn_tide_cnames(2) =   'S2'
+      sn_tide_cnames(3) =   'N2'
+      sn_tide_cnames(4) =   'K2'
+      sn_tide_cnames(5) =   'K1'
+      sn_tide_cnames(6) =   'O1'
+      sn_tide_cnames(7) =   'P1'
+      sn_tide_cnames(8) =   'Q1'
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/BDY/bdydta.F90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/BDY/bdydta.F90	(revision 12116)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/BDY/bdydta.F90	(revision 12117)
@@ -23,5 +23,5 @@
    USE phycst         ! physical constants
    USE sbcapr         ! atmospheric pressure forcing
-   USE sbctide        ! Tidal forcing or not
+   USE tide_mod, ONLY: ln_tide ! tidal forcing
    USE bdy_oce        ! ocean open boundary conditions  
    USE bdytides       ! tidal forcing at boundaries
@@ -75,5 +75,5 @@
 CONTAINS
 
-   SUBROUTINE bdy_dta( kt, kit, kt_offset )
+   SUBROUTINE bdy_dta( kt, kt_offset )
       !!----------------------------------------------------------------------
       !!                   ***  SUBROUTINE bdy_dta  ***
@@ -85,7 +85,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in)           ::   kt           ! ocean time-step index 
-      INTEGER, INTENT(in), OPTIONAL ::   kit          ! subcycle time-step index (for timesplitting option)
-      INTEGER, INTENT(in), OPTIONAL ::   kt_offset    ! time offset in units of timesteps. NB. if kit
-      !                                               ! is present then units = subcycle timesteps.
+      INTEGER, INTENT(in), OPTIONAL ::   kt_offset    ! time offset in units of timesteps
       !                                               ! kt_offset = 0 => get data at "now" time level
       !                                               ! kt_offset = -1 => get data at "before" time level
@@ -105,5 +103,5 @@
       ! Initialise data arrays once for all from initial conditions where required
       !---------------------------------------------------------------------------
-      IF( kt == nit000 .AND. .NOT.PRESENT(kit) ) THEN
+      IF( kt == nit000 ) THEN
 
          ! Calculate depth-mean currents
@@ -216,5 +214,5 @@
          ! read/update all bdy data
          ! ------------------------
-         CALL fld_read( kt, 1, bf_alias, kit = kit, kt_offset = kt_offset )
+         CALL fld_read( kt, 1, bf_alias, kt_offset = kt_offset )
 
          ! apply some corrections in some specific cases...
@@ -275,10 +273,4 @@
             END DO
          ENDIF   ! ltotvel
-
-         ! update tidal harmonic forcing
-         IF( PRESENT(kit) .AND. nn_dyn2d_dta(jbdy) .GE. 2 ) THEN
-            CALL bdytide_update( kt = kt, idx = idx_bdy(jbdy), dta = dta_alias, td = tides(jbdy),   & 
-               &                 kit = kit, kt_offset = kt_offset )
-         ENDIF
 
          !  atm surface pressure : add inverted barometer effect to ssh if it was read
Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/BDY/bdyini.F90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/BDY/bdyini.F90	(revision 12116)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/BDY/bdyini.F90	(revision 12117)
@@ -22,5 +22,5 @@
    USE bdydta         ! open boundary cond. setting   (bdy_dta_init routine)
    USE bdytides       ! open boundary cond. setting   (bdytide_init routine)
-   USE sbctide        ! Tidal forcing or not
+   USE tide_mod, ONLY: ln_tide ! tidal forcing
    USE phycst   , ONLY: rday
    !
Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/BDY/bdytides.F90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/BDY/bdytides.F90	(revision 12116)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/BDY/bdytides.F90	(revision 12117)
@@ -18,5 +18,5 @@
    USE phycst         ! physical constants
    USE bdy_oce        ! ocean open boundary conditions
-   USE tideini        ! 
+   USE tide_mod       ! 
    USE daymod         ! calendar
    !
@@ -30,5 +30,4 @@
 
    PUBLIC   bdytide_init     ! routine called in bdy_init
-   PUBLIC   bdytide_update   ! routine called in bdy_dta
    PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts
 
@@ -45,4 +44,6 @@
    TYPE(OBC_DATA)  , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component)
 
+   INTEGER ::   kt_tide
+
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -64,5 +65,4 @@
       CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
       LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data
-      LOGICAL                                   ::   ln_bdytide_conj     !: If true, assume complex conjugate tidal data
       !!
       INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices
@@ -77,5 +77,5 @@
       TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
       !!
-      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
+      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta
       !!----------------------------------------------------------------------
       !
@@ -105,10 +105,9 @@
             IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
             IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta
-            IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj
             IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo
             IF(lwp) THEN 
                     WRITE(numout,*) '             Tidal components: ' 
                DO itide = 1, nb_harmo
-                  WRITE(numout,*)  '                 ', Wave(ntide(itide))%cname_tide 
+                  WRITE(numout,*)  '                 ', tide_harmonics(itide)%cname_tide 
                END DO
             ENDIF 
@@ -151,6 +150,6 @@
                igrd = 1                       ! Everything is at T-points here
                DO itide = 1, nb_harmo
-                  CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
-                  CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 
+                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) )
+                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) ) 
                   DO ib = 1, ilen0(igrd)
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -168,6 +167,6 @@
                igrd = 2                       ! Everything is at U-points here
                DO itide = 1, nb_harmo
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
+                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) )
+                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:) )
                   DO ib = 1, ilen0(igrd)
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -185,6 +184,6 @@
                igrd = 3                       ! Everything is at V-points here
                DO itide = 1, nb_harmo
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
+                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) )
+                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) )
                   DO ib = 1, ilen0(igrd)
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -210,5 +209,5 @@
                DO itide = 1, nb_harmo
                   !                                                              ! SSH fields
-                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
+                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc'
                   CALL iom_open( clfile, inum )
                   CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
@@ -218,5 +217,5 @@
                   CALL iom_close( inum )
                   !                                                              ! U fields
-                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
+                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc'
                   CALL iom_open( clfile, inum )
                   CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
@@ -226,5 +225,5 @@
                   CALL iom_close( inum )
                   !                                                              ! V fields
-                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
+                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc'
                   CALL iom_open( clfile, inum )
                   CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
@@ -240,10 +239,4 @@
             ENDIF ! ln_bdytide_2ddta=.true.
             !
-            IF( ln_bdytide_conj ) THEN    ! assume complex conjugate in data files
-               td%ssh0(:,:,2) = - td%ssh0(:,:,2)
-               td%u0  (:,:,2) = - td%u0  (:,:,2)
-               td%v0  (:,:,2) = - td%v0  (:,:,2)
-            ENDIF
-            !
             ! Allocate slow varying data in the case of time splitting:
             ! Do it anyway because at this stage knowledge of free surface scheme is unknown
@@ -260,93 +253,4 @@
       !
    END SUBROUTINE bdytide_init
-
-
-   SUBROUTINE bdytide_update( kt, idx, dta, td, kit, kt_offset )
-      !!----------------------------------------------------------------------
-      !!                 ***  SUBROUTINE bdytide_update  ***
-      !!                
-      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays. 
-      !!                
-      !!----------------------------------------------------------------------
-      INTEGER          , INTENT(in   ) ::   kt          ! Main timestep counter
-      TYPE(OBC_INDEX)  , INTENT(in   ) ::   idx         ! OBC indices
-      TYPE(OBC_DATA)   , INTENT(inout) ::   dta         ! OBC external data
-      TYPE(TIDES_DATA) , INTENT(inout) ::   td          ! tidal harmonics data
-      INTEGER, OPTIONAL, INTENT(in   ) ::   kit         ! Barotropic timestep counter (for timesplitting option)
-      INTEGER, OPTIONAL, INTENT(in   ) ::   kt_offset   ! time offset in units of timesteps. NB. if kit
-      !                                                 ! is present then units = subcycle timesteps.
-      !                                                 ! kt_offset = 0  => get data at "now"    time level
-      !                                                 ! kt_offset = -1 => get data at "before" time level
-      !                                                 ! kt_offset = +1 => get data at "after"  time level
-      !                                                 ! etc.
-      !
-      INTEGER  ::   itide, igrd, ib       ! dummy loop indices
-      INTEGER  ::   time_add              ! time offset in units of timesteps
-      INTEGER, DIMENSION(3) ::   ilen0    ! length of boundary data (from OBC arrays)
-      REAL(wp) ::   z_arg, z_sarg, zflag, zramp   ! local scalars    
-      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
-      !!----------------------------------------------------------------------
-      !
-      ilen0(1) =  SIZE(td%ssh(:,1,1))
-      ilen0(2) =  SIZE(td%u(:,1,1))
-      ilen0(3) =  SIZE(td%v(:,1,1))
-
-      zflag=1
-      IF ( PRESENT(kit) ) THEN
-        IF ( kit /= 1 ) zflag=0
-      ENDIF
-
-      IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN
-        !
-        kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
-        !
-        IF(lwp) THEN
-           WRITE(numout,*)
-           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt
-           WRITE(numout,*) '~~~~~~~~~~~~~~ '
-        ENDIF
-        !
-        CALL tide_init_elevation ( idx, td )
-        CALL tide_init_velocities( idx, td )
-        !
-      ENDIF 
-
-      time_add = 0
-      IF( PRESENT(kt_offset) ) THEN
-         time_add = kt_offset
-      ENDIF
-         
-      IF( PRESENT(kit) ) THEN  
-         z_arg = ((kt-kt_tide) * rdt + (kit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) )
-      ELSE                              
-         z_arg = ((kt-kt_tide)+time_add) * rdt
-      ENDIF
-
-      ! Linear ramp on tidal component at open boundaries 
-      zramp = 1._wp
-      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp)
-
-      DO itide = 1, nb_harmo
-         z_sarg = z_arg * omega_tide(itide)
-         z_cost(itide) = COS( z_sarg )
-         z_sist(itide) = SIN( z_sarg )
-      END DO
-
-      DO itide = 1, nb_harmo
-         igrd=1                              ! SSH on tracer grid
-         DO ib = 1, ilen0(igrd)
-            dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide))
-         END DO
-         igrd=2                              ! U grid
-         DO ib = 1, ilen0(igrd)
-            dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide))
-         END DO
-         igrd=3                              ! V grid
-         DO ib = 1, ilen0(igrd) 
-            dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide))
-         END DO
-      END DO
-      !
-   END SUBROUTINE bdytide_update
 
 
@@ -392,5 +296,5 @@
       ! Linear ramp on tidal component at open boundaries 
       zramp = 1.
-      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.)
+      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rn_tide_ramp_dt*rday),0.),1.)
 
       DO ib_bdy = 1,nb_bdy
@@ -433,5 +337,5 @@
             DO itide = 1, nb_harmo
                !
-               z_sarg = (z_arg + zoff) * omega_tide(itide)
+               z_sarg = (z_arg + zoff) * tide_harmonics(itide)%omega
                z_cost = zramp * COS( z_sarg )
                z_sist = zramp * SIN( z_sarg )
@@ -491,6 +395,6 @@
          END DO
          DO ib = 1 , ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*ftide(itide)
-            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
+            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u
          ENDDO
          DO ib = 1 , ilen0(igrd)
@@ -530,6 +434,6 @@
          END DO
          DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*ftide(itide)
-            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
+            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
          ENDDO
          DO ib = 1, ilen0(igrd)
@@ -551,6 +455,6 @@
          END DO
          DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*ftide(itide)
-            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
+            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
          ENDDO
          DO ib = 1, ilen0(igrd)
Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/DYN/dynspg.F90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/DYN/dynspg.F90	(revision 12116)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/DYN/dynspg.F90	(revision 12117)
@@ -21,6 +21,5 @@
    USE dynspg_exp     ! surface pressure gradient     (dyn_spg_exp routine)
    USE dynspg_ts      ! surface pressure gradient     (dyn_spg_ts  routine)
-   USE sbctide        ! 
-   USE updtide        ! 
+   USE tide_mod       !
    USE trd_oce        ! trends: ocean variables
    USE trddyn         ! trend manager: dynamics
@@ -43,4 +42,6 @@
    INTEGER, PARAMETER ::   np_EXP = 0   !       explicit time stepping
    INTEGER, PARAMETER ::   np_NO  =-1   ! no surface pressure gradient, no scheme
+   !
+   REAL(wp) ::   zt0step !   Time of day at the beginning of the time step
 
    !! * Substitutions
@@ -113,5 +114,7 @@
          IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. ln_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case
             !
-            CALL upd_tide( kt )                      ! update tide potential
+            ! Update tide potential at the beginning of current time step
+            zt0step = REAL(nsec_day, wp)-0.5_wp*rdt
+            CALL upd_tide(zt0step)
             !
             DO jj = 2, jpjm1                         ! add tide potential forcing
Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/DYN/dynspg_ts.F90	(revision 12116)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/DYN/dynspg_ts.F90	(revision 12117)
@@ -44,6 +44,5 @@
    USE bdytides        ! open boundary condition data
    USE bdydyn2d        ! open boundary conditions on barotropic variables
-   USE sbctide         ! tides
-   USE updtide         ! tide potential
+   USE tide_mod        !
    USE sbcwave         ! surface wave
    USE diatmb          ! Top,middle,bottom output
@@ -170,4 +169,5 @@
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2    ! averages over the sub-steps of zuwdmask and zvwdmask
+      REAL(wp) ::   zt0substep !   Time of day at the beginning of the time substep
       !!----------------------------------------------------------------------
       !
@@ -424,5 +424,9 @@
          !
          IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, kt_offset= noffset+1 )
-         IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, kt_offset= noffset   )
+         ! Update tide potential at the beginning of current time substep
+         IF( ln_tide_pot .AND. ln_tide ) THEN
+            zt0substep = REAL(nsec_day, wp) - 0.5_wp*rdt + (jn + noffset - 1) * rdt / REAL(nn_baro, wp)
+            CALL upd_tide(zt0substep)
+         END IF
          !
          !                    !==  extrapolation at mid-step  ==!   (jn+1/2)
Index: MO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/SBC/sbctide.F90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/SBC/sbctide.F90	(revision 12116)
+++ 	(revision )
@@ -1,168 +1,0 @@
-MODULE sbctide
-   !!======================================================================
-   !!                       ***  MODULE  sbctide  ***
-   !! Initialization of tidal forcing
-   !!======================================================================
-   !! History :  9.0  !  2007  (O. Le Galloudec)  Original code
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and tracers variables
-   USE dom_oce        ! ocean space and time domain
-   USE phycst         ! physical constant
-   USE daymod         ! calandar
-   USE tideini        ! 
-   !
-   USE in_out_manager ! I/O units
-   USE iom            ! xIOs server
-   USE ioipsl         ! NetCDF IPSL library
-   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-
-   IMPLICIT NONE
-   PUBLIC
-
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   pot_astro   !
-
-   !!----------------------------------------------------------------------
-   !!   tidal potential
-   !!----------------------------------------------------------------------
-   !!   sbc_tide            : 
-   !!   tide_init_potential :
-   !!----------------------------------------------------------------------
-
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   amp_pot, phi_pot
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   amp_load, phi_load
- 
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE sbc_tide( kt )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE sbc_tide  ***
-      !!----------------------------------------------------------------------      
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step
-      INTEGER               ::   jk     ! dummy loop index
-      INTEGER               ::   nsec_day_orig     ! Temporary variable
-      !!----------------------------------------------------------------------
-      
-      IF( nsec_day == NINT(0.5_wp * rdt) .OR. kt == nit000 ) THEN      ! start a new day
-         !
-         IF( kt == nit000 )THEN
-            ALLOCATE( amp_pot(jpi,jpj,nb_harmo),                      &
-               &      phi_pot(jpi,jpj,nb_harmo), pot_astro(jpi,jpj)   )
-            IF( ln_read_load )THEN
-               ALLOCATE( amp_load(jpi,jpj,nb_harmo), phi_load(jpi,jpj,nb_harmo) )
-               CALL tide_init_load
-            ENDIF
-         ENDIF
-         !
-         IF( ln_read_load )THEN
-            amp_pot(:,:,:) = amp_load(:,:,:)
-            phi_pot(:,:,:) = phi_load(:,:,:)
-         ELSE 
-            amp_pot(:,:,:) = 0._wp
-            phi_pot(:,:,:) = 0._wp
-         ENDIF
-         pot_astro(:,:) = 0._wp
-         !
-         ! If the run does not start from midnight then need to initialise tides
-         ! at the start of the current day (only occurs when kt==nit000)
-         ! Temporarily set nsec_day to beginning of day.
-         nsec_day_orig = nsec_day
-         IF ( nsec_day /= NINT(0.5_wp * rdt) ) THEN 
-            kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
-            nsec_day = NINT(0.5_wp * rdt)
-         ELSE
-            kt_tide = kt 
-         ENDIF
-         CALL tide_harmo( omega_tide, v0tide, utide, ftide, ntide, nb_harmo )
-         !
-         !
-         IF(lwp) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) 'sbc_tide : Update of the components and (re)Init. the potential at kt=', kt
-            WRITE(numout,*) '~~~~~~~~ '
-            DO jk = 1, nb_harmo
-               WRITE(numout,*) Wave(ntide(jk))%cname_tide, utide(jk), ftide(jk), v0tide(jk), omega_tide(jk)
-            END DO
-         ENDIF
-         !
-         IF( ln_tide_pot )   CALL tide_init_potential
-         !
-         ! Reset nsec_day
-         nsec_day = nsec_day_orig 
-      ENDIF
-      !
-   END SUBROUTINE sbc_tide
-
-
-   SUBROUTINE tide_init_potential
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE tide_init_potential  ***
-      !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zcons, ztmp1, ztmp2, zlat, zlon, ztmp, zamp, zcs   ! local scalar
-      !!----------------------------------------------------------------------
-
-      DO jk = 1, nb_harmo
-         zcons = 0.7_wp * Wave(ntide(jk))%equitide * ftide(jk)
-         DO ji = 1, jpi
-            DO jj = 1, jpj
-               ztmp1 =  ftide(jk) * amp_pot(ji,jj,jk) * COS( phi_pot(ji,jj,jk) + v0tide(jk) + utide(jk) )
-               ztmp2 = -ftide(jk) * amp_pot(ji,jj,jk) * SIN( phi_pot(ji,jj,jk) + v0tide(jk) + utide(jk) )
-               zlat = gphit(ji,jj)*rad !! latitude en radian
-               zlon = glamt(ji,jj)*rad !! longitude en radian
-               ztmp = v0tide(jk) + utide(jk) + Wave(ntide(jk))%nutide * zlon
-               ! le potentiel est composé des effets des astres:
-               IF    ( Wave(ntide(jk))%nutide == 1 )  THEN  ;  zcs = zcons * SIN( 2._wp*zlat )
-               ELSEIF( Wave(ntide(jk))%nutide == 2 )  THEN  ;  zcs = zcons * COS( zlat )**2
-               ELSE                                         ;  zcs = 0._wp
-               ENDIF
-               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(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
-         CALL iom_get  ( inum, jpdom_data,TRIM(Wave(ntide(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
-
-  !!======================================================================
-END MODULE sbctide
Index: MO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/SBC/tide.h90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/SBC/tide.h90	(revision 12116)
+++ 	(revision )
@@ -1,29 +1,0 @@
-   !!----------------------------------------------------------------------
-   !! History :  3.2  !  2007  (O. Le Galloudec)  Original code
-   !!----------------------------------------------------------------------
-
-   !             !! name_tide , equitide , nutide , nt , ns , nh , np , np1 , shift , nksi , nnu0 , nnu1 , nnu2 , R , formula !!
-   !             !!           !          !        !    !    !    !    !     !       !      !      !      !      !   !         !!
-   Wave( 1) = tide(  'M2'     , 0.242297 ,    2   ,  2 , -2 ,  2 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
-   Wave( 2) = tide(  'N2'     , 0.046313 ,    2   ,  2 , -3 ,  2 ,  1 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
-   Wave( 3) = tide( '2N2'     , 0.006184 ,    2   ,  2 , -4 ,  2 ,  2 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
-   Wave( 4) = tide(  'S2'     , 0.113572 ,    2   ,  2 ,  0 ,  0 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,     0   )
-   Wave( 5) = tide(  'K2'     , 0.030875 ,    2   ,  2 ,  0 ,  2 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   , -2   , 0 ,   235   )
-   !              !           !          !        !    !    !    !    !     !       !      !      !      !      !   !         !
-   Wave( 6) = tide(  'K1'     , 0.142408 ,    1   ,  1 ,  0 ,  1 ,  0 ,  0  ,  -90  ,  0   ,  0   , -1   ,  0   , 0 ,   227   )
-   Wave( 7) = tide(  'O1'     , 0.101266 ,    1   ,  1 , -2 ,  1 ,  0 ,  0  ,  +90  ,  2   , -1   ,  0   ,  0   , 0 ,    75   )
-   Wave( 8) = tide(  'Q1'     , 0.019387 ,    1   ,  1 , -3 ,  1 ,  1 ,  0  ,  +90  ,  2   , -1   ,  0   ,  0   , 0 ,    75   )
-   Wave( 9) = tide(  'P1'     , 0.047129 ,    1   ,  1 ,  0 , -1 ,  0 ,  0  ,  +90  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    )
-   !              !           !          !        !    !    !    !    !     !       !      !      !      !      !   !         !
-   Wave(10) = tide(  'M4'     , 0.000000 ,    4   ,  4 , -4 ,  4 ,  0 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,    1    )
-   !              !           !          !        !    !    !    !    !     !       !      !      !      !      !   !         !
-   Wave(11) = tide(  'Mf'     , 0.042017 ,    0   ,  0 ,  2 ,  0 ,  0 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    )
-   Wave(12) = tide(  'Mm'     , 0.022191 ,    0   ,  0 ,  1 ,  0 , -1 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,   73    )
-   Wave(13) = tide(  'Msqm'   , 0.000667 ,    0   ,  0 ,  4 , -2 ,  0 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    )
-   Wave(14) = tide(  'Mtm'    , 0.008049 ,    0   ,  0 ,  3 ,  0 , -1 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    )
-   !              !           !          !        !    !    !    !    !     !       !      !      !      !      !   !         !
-   Wave(15) = tide(  'S1'     , 0.000000 ,    1   ,  1 ,  0 ,  0 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    )   
-   Wave(16) = tide(  'MU2'    , 0.005841 ,    2   ,  2 , -4 ,  4 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,   78    )
-   Wave(17) = tide(  'NU2'    , 0.009094 ,    2   ,  2 , -3 ,  4 , -1 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,   78    ) 
-   Wave(18) = tide(  'L2'     , 0.006694 ,    2   ,  2 , -1 ,  2 , -1 ,  0  , +180  ,  2   , -2   ,  0   ,  0   , 0 ,  215    )
-   Wave(19) = tide(  'T2'     , 0.006614 ,    2   ,  2 ,  0 , -1 ,  0 ,  1  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    )
Index: MO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/SBC/tideini.F90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/SBC/tideini.F90	(revision 12116)
+++ 	(revision )
@@ -1,134 +1,0 @@
-MODULE tideini
-   !!======================================================================
-   !!                       ***  MODULE  tideini  ***
-   !! Initialization of tidal forcing
-   !!======================================================================
-   !! History :  1.0  !  2007  (O. Le Galloudec)  Original code
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and tracers variables
-   USE dom_oce        ! ocean space and time domain
-   USE phycst         ! physical constants
-   USE daymod         ! calendar
-   USE tide_mod       ! 
-   !
-   USE in_out_manager ! I/O units
-   USE iom            ! xIOs server
-   USE ioipsl         ! NetCDF IPSL library
-   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-
-   IMPLICIT NONE
-   PUBLIC
-
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   omega_tide   !:
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   v0tide       !:
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   utide        !:
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   ftide        !:
-
-   LOGICAL , PUBLIC ::   ln_tide         !:
-   LOGICAL , PUBLIC ::   ln_tide_pot     !:
-   LOGICAL , PUBLIC ::   ln_read_load    !:
-   LOGICAL , PUBLIC ::   ln_scal_load    !:
-   LOGICAL , PUBLIC ::   ln_tide_ramp    !:
-   INTEGER , PUBLIC ::   nb_harmo        !:
-   INTEGER , PUBLIC ::   kt_tide         !:
-   REAL(wp), PUBLIC ::   rdttideramp     !:
-   REAL(wp), PUBLIC ::   rn_scal_load    !:
-   CHARACTER(lc), PUBLIC ::   cn_tide_load   !: 
-
-   INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:) ::   ntide   !:
-
-   !!----------------------------------------------------------------------
-   !! 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) :: clname
-      INTEGER  ::   ios                 ! Local integer output status for namelist read
-      !
-      NAMELIST/nam_tide/ln_tide, ln_tide_pot, ln_scal_load, ln_read_load, cn_tide_load, &
-                  &     ln_tide_ramp, rn_scal_load, rdttideramp, clname
-      !!----------------------------------------------------------------------
-      !
-      ! 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' )
-      !
-      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' )
-      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,*) '         Apply astronomical potential            ln_tide_pot  = ', ln_tide_pot
-            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                 rdttideramp  = ', rdttideramp
-         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
-      !
-      CALL tide_init_Wave
-      !
-      nb_harmo=0
-      DO jk = 1, jpmax_harmo
-         DO ji = 1,jpmax_harmo
-            IF( TRIM(clname(jk)) == Wave(ji)%cname_tide )   nb_harmo = nb_harmo + 1
-         END DO
-      END DO
-      !       
-      ! 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( 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 < rdttideramp) )   &
-         &   CALL ctl_stop('rdttideramp must be lower than run duration')
-      IF( ln_tide_ramp.AND.(rdttideramp<0.) ) &
-         &   CALL ctl_stop('rdttideramp must be positive')
-      !
-      ALLOCATE( ntide(nb_harmo) )
-      DO jk = 1, nb_harmo
-         DO ji = 1, jpmax_harmo
-            IF( TRIM(clname(jk)) == Wave(ji)%cname_tide ) THEN
-               ntide(jk) = ji
-               EXIT
-            ENDIF
-         END DO
-      END DO
-      !
-      ALLOCATE( omega_tide(nb_harmo), v0tide    (nb_harmo),   &
-         &      utide     (nb_harmo), ftide     (nb_harmo)  )
-      kt_tide = nit000
-      !
-      IF (.NOT.ln_scal_load ) rn_scal_load = 0._wp
-      !
-   END SUBROUTINE tide_init
-     
-   !!======================================================================
-END MODULE tideini
Index: MO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/SBC/updtide.F90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/SBC/updtide.F90	(revision 12116)
+++ 	(revision )
@@ -1,81 +1,0 @@
-MODULE updtide
-   !!======================================================================
-   !!                       ***  MODULE  updtide  ***
-   !! Initialization of tidal forcing
-   !!======================================================================
-   !! History :  9.0  !  07  (O. Le Galloudec)  Original code
-   !!----------------------------------------------------------------------
-   !!   upd_tide       : update tidal potential
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain
-   USE in_out_manager  ! I/O units
-   USE phycst          ! physical constant
-   USE sbctide         ! tide potential variable
-   USE tideini, ONLY: ln_tide_ramp, rdttideramp
-
-   IMPLICIT NONE
-   PUBLIC
-
-   PUBLIC   upd_tide   ! called in dynspg_... modules
-  
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE upd_tide( kt, kit, kt_offset )
-      !!----------------------------------------------------------------------
-      !!                 ***  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
-      !!----------------------------------------------------------------------      
-      INTEGER, INTENT(in)           ::   kt      ! ocean time-step index
-      INTEGER, INTENT(in), OPTIONAL ::   kit     ! external mode sub-time-step index (lk_dynspg_ts=T)
-      INTEGER, INTENT(in), OPTIONAL ::   kt_offset ! time offset in number 
-                                                     ! of internal steps             (lk_dynspg_ts=F)
-                                                     ! of external steps             (lk_dynspg_ts=T)
-      !
-      INTEGER  ::   ioffset      ! local integer
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zt, zramp    ! local scalar
-      REAL(wp), DIMENSION(nb_harmo) ::   zwt 
-      !!----------------------------------------------------------------------      
-      !
-      !                               ! tide pulsation at model time step (or sub-time-step)
-      zt = ( kt - kt_tide ) * rdt
-      !
-      ioffset = 0
-      IF( PRESENT( kt_offset ) )   ioffset = kt_offset
-      !
-      IF( PRESENT( kit ) )   THEN
-         zt = zt + ( kit +  ioffset - 1 ) * rdt / REAL( nn_baro, wp )
-      ELSE
-         zt = zt + ioffset * rdt
-      ENDIF
-      !
-      zwt(:) = omega_tide(:) * zt
-
-      pot_astro(:,:) = 0._wp          ! update tidal potential (sum of all harmonics)
-      DO jk = 1, nb_harmo   
-         pot_astro(:,:) = pot_astro(:,:) + amp_pot(:,:,jk) * COS( zwt(jk) + phi_pot(:,:,jk) )      
-      END DO
-      !
-      IF( ln_tide_ramp ) THEN         ! linear increase if asked
-         zt = ( kt - nit000 ) * rdt
-         IF( PRESENT( kit ) )   zt = zt + ( kit + ioffset -1) * rdt / REAL( nn_baro, wp )
-         zramp = MIN(  MAX( zt / (rdttideramp*rday) , 0._wp ) , 1._wp  )
-         pot_astro(:,:) = zramp * pot_astro(:,:)
-      ENDIF
-      !
-   END SUBROUTINE upd_tide
-
-  !!======================================================================
-
-END MODULE updtide
Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/TDE/tide.h90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/TDE/tide.h90	(revision 12117)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/TDE/tide.h90	(revision 12117)
@@ -0,0 +1,166 @@
+   !!=====================================================================
+   !!                  ***  Include file tide.h90  ***
+   !!======================================================================
+   !! History :  3.2  !  2007  (O. Le Galloudec)  Original code
+   !!                 !  2019  (S. Mueller, N. Bruneau)  Alternative parameter set
+   !!----------------------------------------------------------------------
+   !!
+   !! ** Purpose : Inclusion of alternative variants of tidal-constituent
+   !!              parameter definitions during code preprocessing: the default
+   !!              variant includes the 34 constituents available in the FES2014
+   !!              version of the Finite Element Solution - Global Tide data
+   !!              product
+   !!              (https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/global-tide-fes.html);
+   !!              also available is the default parameter set available in
+   !!              previous NEMO versions
+   !!
+   !! ** References :
+   !!       S58) Schureman, P. (1958): Manual of Harmonic Analysis and
+   !!            Prediction of Tides (Revised (1940) Edition (Reprinted 1958
+   !!            with corrections). Reprinted June 2001). U.S. Department of
+   !!            Commerce, Coast and Geodetic Survey Special Publication
+   !!            No. 98. Washington DC, United States Government Printing
+   !!            Office. 317 pp. DOI: 10.25607/OBP-155.
+   !!      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
+   !!      CE73) Cartwright, D. E. and Edden, A. C. (1973): Corrected Tables of
+   !!            Tidal Harmonics. Geophys. J. R. astr. Soc. 33,
+   !!            pp. 253-264. DOI: 10.1111/j.1365-246X.1973.tb03420.x
+   !!   FES2014) FES (Finite element Solution) - Global
+   !!            tide. https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/global-tide-fes.html
+   !!
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2019)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+#ifndef TIDE_VAR_0
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   !                         | Name   | Equilibrium  | nt | ns | nh | np | np1 | Phase | nxi | nnu0 | nnu1 | nnu2 |  R | Nodal      | Equilibrium    | Parameters source     | Notes  |
+   !                         |        | tide         |    |    |    |    |     | shift |     |      |      |      |    | correction | tide           |                       |        |
+   !                         |        |              |    |    |    |    |     |       |     |      |      |      |    | formula    | source/comment |                       |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   !                         |                                                            Long-period tidal constituents                                                              |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components( 1) = tide( 'Mf'   ,  0.042054_wp ,  0 ,  2 ,  1 ,  0 ,   0 ,     0 ,  -2 ,    0 ,    0 ,    0 ,  0 ,  74 )      ! CE73           | S54 (Table 2, A6)     |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components( 2) = tide( 'Mm'   ,  0.022187_wp ,  0 ,  1 ,  0 , -1 ,   0 ,     0 ,   0 ,    0 ,    0 ,    0 ,  0 ,  73 )      ! CE73           | S54 (Table 2, A2)     |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components( 3) = tide( 'Ssa'  ,  0.019572_wp ,  0 ,  0 ,  2 ,  0 ,   0 ,     0 ,   0 ,    0 ,    0 ,    0 ,  0 ,   0 )      ! CE73           | S54 (Table 2, B6)     |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components( 4) = tide( 'Mtm'  ,  0.008052_wp ,  0 ,  3 ,  0 , -1 ,   0 ,     0 ,  -2 ,    0 ,    0 ,    0 ,  0 ,  74 )      ! CE73           | FES2014 (prediction   |        |
+   !                         |        |              |    |    |    |    |     |       |     |      |      |      |    |            |                | algorithm); S54       |        |
+   !                         |        |              |    |    |    |    |     |       |     |      |      |      |    |            |                | (Table 2, A7)         |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components( 5) = tide( 'Msf'  ,  0.003677_wp ,  0 ,  2 , -2 ,  0 ,   0 ,     0 ,   0 ,    0 ,    0 ,    0 ,  0 ,  73 )      ! CE73           | S54 (Table 2, A5)     |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components( 6) = tide( 'Msqm' ,  0.001287_wp ,  0 ,  4 , -2 ,  0 ,   0 ,     0 ,  -2 ,    0 ,    0 ,    0 ,  0 ,  74 )      ! CE73           | FES2014 (prediction   |        |
+   !                         |        |              |    |    |    |    |     |       |     |      |      |      |    |            |                | algorithm); S54       |        |
+   !                         |        |              |    |    |    |    |     |       |     |      |      |      |    |            |                | (Table 2, A12)        |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components( 7) = tide( 'Sa'   ,  0.000000_wp ,  0 ,  0 ,  1 ,  0 ,   0 ,     0 ,   0 ,    0 ,    0 ,    0 ,  0 ,   0 )      ! Meteorological | S54 (Table 2, B64)    |        |
+   !                         |        |              |    |    |    |    |     |       |     |      |      |      |    |            | tide only      |                       |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   !                         |                                                              Diurnal tidal constituents                                                                |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components( 8) = tide( 'K1'   , -0.142486_wp ,  1 ,  0 ,  1 ,  0 ,   0 ,    90 ,   0 ,    0 ,   -1 ,    0 ,  0 , 227 )      ! CE73           | S54 (Table 2)         | Note 1 |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components( 9) = tide( 'O1'   ,  0.101316_wp ,  1 , -2 ,  1 ,  0 ,   0 ,    90 ,   2 ,   -1 ,    0 ,    0 ,  0 ,  75 )      ! CE73           | S54 (Table 2, A14)    |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(10) = tide( 'P1'   ,  0.047152_wp ,  1 ,  0 , -1 ,  0 ,   0 ,    90 ,   0 ,    0 ,    0 ,    0 ,  0 ,   0 )      ! CE73           | S54 (Table 2, B14)    |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(11) = tide( 'Q1'   ,  0.019396_wp ,  1 , -3 ,  1 ,  1 ,   0 ,    90 ,   2 ,   -1 ,    0 ,    0 ,  0 ,  75 )      ! CE73           | S54 (Table 2, A15)    |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(12) = tide( 'J1'   , -0.007967_wp ,  1 ,  1 ,  1 , -1 ,   0 ,    90 ,   0 ,   -1 ,    0 ,    0 ,  0 ,  76 )      ! CE73           | S54 (Table 2, A24)    | Note 1 |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(13) = tide( 'S1'   ,  0.000000_wp ,  1 ,  0 ,  0 ,  0 ,   0 ,     0 ,   0 ,    0 ,    0 ,    0 ,  0 ,   0 )      ! Meteorological | S54 (Table 2, B71)    |        |
+   !                         |        |              |    |    |    |    |     |       |     |      |      |      |    |            | tide only      |                       |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   !                         |                                                            Semidiurnal tidal constituents                                                              |
+   !                         +--------+-------------+-----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(14) = tide( 'M2'   ,  0.244081_wp ,  2 , -2 ,  2 ,  0 ,   0 ,     0 ,   2 ,   -2 ,    0 ,    0 ,  0 ,  78 )      ! CE73           | S54 (Table 2, A39)    |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(15) = tide( 'S2'   ,  0.110242_wp ,  2 ,  0 ,  0 ,  0 ,   0 ,     0 ,   0 ,    0 ,    0 ,    0 ,  0 ,   0 )      ! CE73           | S54 (Table 2, B39)    |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(16) = tide( 'N2'   ,  0.046732_wp ,  2 , -3 ,  2 ,  1 ,   0 ,     0 ,   2 ,   -2 ,    0 ,    0 ,  0 ,  78 )      ! CE73           | S54 (Table 2, A40)    |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(17) = tide( 'K2'   ,  0.030905_wp ,  2 ,  0 ,  2 ,  0 ,   0 ,     0 ,   0 ,    0 ,    0 ,   -2 ,  0 , 235 )      ! CE73           | S54 (Table 2)         |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(18) = tide( 'nu2'  ,  0.008877_wp ,  2 , -3 ,  4 , -1 ,   0 ,     0 ,   2 ,   -2 ,    0 ,    0 ,  0 ,  78 )      ! CE73           | S54 (Table 2, A43)    |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(19) = tide( 'mu2'  ,  0.007463_wp ,  2 , -4 ,  4 ,  0 ,   0 ,     0 ,   2 ,   -2 ,    0 ,    0 ,  0 ,  78 )      ! CE73           | S54 (Table 2, A45)    |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(20) = tide( '2N2'  ,  0.006184_wp ,  2 , -4 ,  2 ,  2 ,   0 ,     0 ,   2 ,   -2 ,    0 ,    0 ,  0 ,  78 )      ! CE73           | S54 (Table 2, A42)    |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(21) = tide( 'L2'   , -0.006899_wp ,  2 , -1 ,  2 , -1 ,   0 ,     0 ,   2 ,   -2 ,    0 ,    0 , -1 , 215 )      ! CE73           | S54 (Table 2)         | Note 1 |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(22) = tide( 'T2'   ,  0.006655_wp ,  2 ,  0 , -1 ,  0 ,   1 ,     0 ,   0 ,    0 ,    0 ,    0 ,  0 ,   0 )      ! CE73           | S54 (Table 2, B40)    |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(23) = tide( 'eps2' ,  0.001804_wp ,  2 , -5 ,  4 ,  1 ,   0 ,     0 ,   2 ,   -2 ,    0 ,    0 ,  0 ,  78 )      ! CE73           | FES2014 (prediction   |        |
+   !                         |        |              |    |    |    |    |     |       |     |      |      |      |    |            |                | algorithm)            |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(24) = tide( 'lam2' , -0.001800_wp ,  2 , -1 ,  0 ,  1 ,   0 ,     0 ,   2 ,   -2 ,    0 ,    0 ,  0 ,  78 )      ! CE73           | S54 (Table 2, A44)    | Note 1 |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(25) = tide( 'R2'   , -0.000952_wp ,  2 ,  0 ,  1 ,  0 ,  -1 ,     0 ,   0 ,    0 ,    0 ,    0 ,  0 ,   0 )      ! CE73           | S54 (Table 2, B41)    | Note 1 |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   !                         |                                                            Terdiurnal tidal constituents                                                               |
+   !                         +--------+-------------+-----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(26) = tide( 'M3'   , -0.003192_wp ,  3 , -3 ,  3 ,  0 ,   0 ,   180 ,   3 ,   -3 ,    0 ,    0 ,  0 , 149 )      ! CT71           | S54 (Table 2, A82)    | Note 2 |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   !                         |                                                                    Compound tides                                                                      |
+   !                         +--------+-------------+-----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(27) = tide( 'MKS2' ,  0.000000_wp ,  2 , -2 ,  4 ,  0 ,   0 ,     0 ,   2 ,   -2 ,    0 ,   -2 ,  0 ,   4 )      ! Compound tide  | FES2014 (prediction   |        |
+   !                         |        |              |    |    |    |    |     |       |     |      |      |      |    |            |                | algorithm)            |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(28) = tide( 'MN4'  ,  0.000000_wp ,  4 , -5 ,  4 ,  1 ,   0 ,     0 ,   4 ,   -4 ,    0 ,    0 ,  0 ,   1 )      ! Compound tide  | S54 (Table 2a)        |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(29) = tide( 'MS4'  ,  0.000000_wp ,  4 , -2 ,  2 ,  0 ,   0 ,     0 ,   2 ,   -2 ,    0 ,    0 ,  0 ,  78 )      ! Compound tide  | FES2014 (prediction   | Note 3 |
+   !                         |        |              |    |    |    |    |     |       |     |      |      |      |    |            |                | algorithm); S54       |        |
+   !                         |        |              |    |    |    |    |     |       |     |      |      |      |    |            |                | (Table 2a)            |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   !                         |                                                                       Overtides                                                                        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(30) = tide( 'M4'   ,  0.000000_wp ,  4 , -4 ,  4 ,  0 ,   0 ,     0 ,   4 ,   -4 ,    0 ,    0 ,  0 ,   1 )      ! Overtide       | S54                   |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(31) = tide( 'N4'   ,  0.000000_wp ,  4 , -6 ,  4 ,  2 ,   0 ,     0 ,   4 ,   -4 ,    0 ,    0 ,  0 ,   1 )      ! Overtide       | FES2014 (prediction   |        |
+   !                         |        |              |    |    |    |    |     |       |     |      |      |      |    |            |                | algorithm)            |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(32) = tide( 'S4'   ,  0.000000_wp ,  4 ,  0 ,  0 ,  0 ,  0  ,     0 ,   0 ,    0 ,    0 ,    0 ,  0 ,   0 )      ! Overtide       | S54                   |        | 
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(33) = tide( 'M6'   ,  0.000000_wp ,  6 , -6 ,  6 ,  0 ,  0  ,     0 ,   6 ,   -6 ,    0 ,    0 ,  0 ,  18 )      ! Overtide       | S54                   |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   tide_components(34) = tide( 'M8'   ,  0.000000_wp ,  8 , -8 ,  8 ,  0 ,  0  ,     0 ,   8 ,   -8 ,    0 ,    0 ,  0 ,  20 )      ! Overtide       | S54                   |        |
+   !                         +--------+--------------+----+----+----+----+-----+-------+-----+------+------+------+----+------------+----------------+-----------------------+--------+
+   !   Note 1: the phase shift from Table 2 of S54 has been adjusted to accomodate the negative sign of the equilibrium-tide value derived from CE73.
+   !   Note 2: the phase shift from Table 2 of S54 has been adjusted to accomodate the negative sign of the equilibrium-tide value derived from CT71.
+   !   Note 3: the nodal correction factor formulas from FES2014 and S54 differ; here, the version from FES2014 has been selected.
+#else
+   !                        !! name_tide , equitide , nt , ns , nh , np , np1 , shift , nksi , nnu0 , nnu1 , nnu2 , R , formula !!
+   !                        !!           !          !    !    !    !    !     !       !      !      !      !      !   !         !!
+   tide_components( 1) = tide(  'M2'     , 0.242297 ,  2 , -2 ,  2 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
+   tide_components( 2) = tide(  'N2'     , 0.046313 ,  2 , -3 ,  2 ,  1 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
+   tide_components( 3) = tide( '2N2'     , 0.006184 ,  2 , -4 ,  2 ,  2 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
+   tide_components( 4) = tide(  'S2'     , 0.113572 ,  2 ,  0 ,  0 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,     0   )
+   tide_components( 5) = tide(  'K2'     , 0.030875 ,  2 ,  0 ,  2 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   , -2   , 0 ,   235   )
+   !                         !           !          !    !    !    !    !     !       !      !      !      !      !   !         !
+   tide_components( 6) = tide(  'K1'     , 0.142408 ,  1 ,  0 ,  1 ,  0 ,  0  ,  -90  ,  0   ,  0   , -1   ,  0   , 0 ,   227   )
+   tide_components( 7) = tide(  'O1'     , 0.101266 ,  1 , -2 ,  1 ,  0 ,  0  ,  +90  ,  2   , -1   ,  0   ,  0   , 0 ,    75   )
+   tide_components( 8) = tide(  'Q1'     , 0.019387 ,  1 , -3 ,  1 ,  1 ,  0  ,  +90  ,  2   , -1   ,  0   ,  0   , 0 ,    75   )
+   tide_components( 9) = tide(  'P1'     , 0.047129 ,  1 ,  0 , -1 ,  0 ,  0  ,  +90  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    )
+   !                         !           !          !    !    !    !    !     !       !      !      !      !      !   !         !
+   tide_components(10) = tide(  'M4'     , 0.000000 ,  4 , -4 ,  4 ,  0 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,    1    )
+   !                         !           !          !    !    !    !    !     !       !      !      !      !      !   !         !
+   tide_components(11) = tide(  'Mf'     , 0.042017 ,  0 ,  2 ,  0 ,  0 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    )
+   tide_components(12) = tide(  'Mm'     , 0.022191 ,  0 ,  1 ,  0 , -1 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,   73    )
+   tide_components(13) = tide(  'Msqm'   , 0.000667 ,  0 ,  4 , -2 ,  0 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    )
+   tide_components(14) = tide(  'Mtm'    , 0.008049 ,  0 ,  3 ,  0 , -1 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    )
+   !                         !           !          !    !    !    !    !     !       !      !      !      !      !   !         !
+   tide_components(15) = tide(  'S1'     , 0.000000 ,  1 ,  0 ,  0 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    )   
+   tide_components(16) = tide(  'MU2'    , 0.005841 ,  2 , -4 ,  4 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,   78    )
+   tide_components(17) = tide(  'NU2'    , 0.009094 ,  2 , -3 ,  4 , -1 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,   78    ) 
+   tide_components(18) = tide(  'L2'     , 0.006694 ,  2 , -1 ,  2 , -1 ,  0  , +180  ,  2   , -2   ,  0   ,  0   , 0 ,  215    )
+   tide_components(19) = tide(  'T2'     , 0.006614 ,  2 ,  0 , -1 ,  0 ,  1  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    )
+#endif
Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/TDE/tide_mod.F90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/TDE/tide_mod.F90	(revision 12117)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/TDE/tide_mod.F90	(revision 12117)
@@ -0,0 +1,763 @@
+MODULE tide_mod
+   !!======================================================================
+   !!                       ***  MODULE  tide_mod  ***
+   !! Compute nodal modulations corrections and pulsations
+   !!======================================================================
+   !! History :  1.0  !  2007  (O. Le Galloudec)  Original code
+   !!                 !  2019  (S. Mueller)
+   !!----------------------------------------------------------------------
+   !!
+   !! ** Reference :
+   !!       S58) Schureman, P. (1958): Manual of Harmonic Analysis and
+   !!            Prediction of Tides (Revised (1940) Edition (Reprinted 1958
+   !!            with corrections). Reprinted June 2001). U.S. Department of
+   !!            Commerce, Coast and Geodetic Survey Special Publication
+   !!            No. 98. Washington DC, United States Government Printing
+   !!            Office. 317 pp. DOI: 10.25607/OBP-155.
+   !!----------------------------------------------------------------------
+
+   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                       !
+
+   ! Longitudes on 1 Jan 1900, 00h and angular velocities (units of deg and
+   ! deg/h, respectively. The values of these module variables have been copied
+   ! from subroutine astronomic_angle of the version of this module used in
+   ! release version 4.0 of NEMO.
+   REAL(wp) ::   rlon00_N  =  259.1560564_wp               ! Longitude of ascending lunar node
+   REAL(wp) ::   romega_N  = -.0022064139_wp
+   REAL(wp) ::   rlon00_T  =  180.0_wp                     ! Mean solar angle (GMT)
+   REAL(wp) ::   romega_T  =  15.0_wp
+   REAL(wp) ::   rlon00_h  =  280.1895014_wp               ! Mean solar Longitude
+   REAL(wp) ::   romega_h  =  .0410686387_wp
+   REAL(wp) ::   rlon00_s  =  277.0256206_wp               ! Mean lunar Longitude
+   REAL(wp) ::   romega_s  =  .549016532_wp
+   REAL(wp) ::   rlon00_p1 =  281.2208569_wp               ! Longitude of solar perigee
+   REAL(wp) ::   romega_p1 =  .000001961_wp
+   REAL(wp) ::   rlon00_p  =  334.3837214_wp               ! Longitude of lunar perigee
+   REAL(wp) ::   romega_p  =  .004641834_wp
+   ! Values of cos(i)*cos(epsilon), rcice, and sin(incl)*sin(epsilon), rsise,
+   ! where i is the inclination of the orbit of the Moon w.r.t. the ecliptic and
+   ! epsilon the obliquity of the ecliptic on 1 January 1900, 00h. The values of
+   ! these module variables have been copied from subroutine astronomic_angle
+   ! (computation of the cosine of inclination of orbit of Moon to the celestial
+   ! equator) of the version of this module used in release version 4.0 of NEMO.
+   REAL(wp) ::   rcice    = 0.913694997_wp
+   REAL(wp) ::   rsise    = 0.035692561_wp
+   ! Coefficients used to compute sh_xi and sh_nu in subroutine astronomic_angle
+   ! according to two equations given in the explanation of Table 6 of S58
+   REAL(wp) ::   rxinu1, rxinu2
+
+   !!----------------------------------------------------------------------
+   !! 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' )
+      !
+      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' )
+      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')
+      !
+      ! Compute coefficients which are used in subroutine astronomic_angle to
+      ! compute sh_xi and sh_nu according to two equations given in the
+      ! explanation of Table 6 of S58
+      rxinu1 = COS( 0.5_wp * ( ABS( ACOS( rcice + rsise ) ) ) ) / COS( 0.5_wp * ( ACOS( rcice - rsise ) ) )
+      rxinu2 = SIN( 0.5_wp * ( ABS( ACOS( rcice + rsise ) ) ) ) / SIN( 0.5_wp * ( ACOS( rcice - rsise ) ) )
+      !
+      ! 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)
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE astronomic_angle  ***
+      !!                      
+      !! ** Purpose : Compute astronomic angles
+      !!----------------------------------------------------------------------
+      INTEGER  ::   psec_day !   Number of seconds from midnight
+      REAL(wp) ::   zp, zq, zt2, zs2, ztgI2, zP1, ztgn2, zat1, zat2
+      REAL(wp) ::   zqy , zsy, zday, zdj, zhfrac, zt
+      !!----------------------------------------------------------------------
+      !
+      ! Computation of the time from 1 Jan 1900, 00h in years
+      zqy = AINT( (nyear - 1901.0_wp) / 4.0_wp )
+      zsy = nyear - 1900.0_wp
+      !
+      zdj  = dayjul( nyear, nmonth, nday )
+      zday = zdj + zqy - 1.0_wp
+      !
+      zhfrac = psec_day / 3600.0_wp
+      !
+      zt = zsy * 365.0_wp * 24.0_wp + zday * 24.0_wp + zhfrac
+      !
+      ! Longitude of ascending lunar node
+      sh_N  = ( rlon00_N  + romega_N * zt     ) * rad
+      sh_N = MOD( sh_N,  2*rpi )
+      ! Mean solar angle (Greenwhich time)
+      sh_T  = ( rlon00_T  + romega_T * zhfrac ) * rad
+      ! Mean solar Longitude
+      sh_h  = ( rlon00_h  + romega_h * zt     ) * rad
+      sh_h = MOD( sh_h,  2*rpi )
+      ! Mean lunar Longitude
+      sh_s  = ( rlon00_s  + romega_s * zt     ) * rad
+      sh_s = MOD( sh_s,  2*rpi )
+      ! Longitude of solar perigee
+      sh_p1 = ( rlon00_p1 + romega_p1 * zt    ) * rad
+      sh_p1= MOD( sh_p1, 2*rpi )
+      ! Longitude of lunar perigee
+      sh_p  = ( rlon00_p  + romega_p  * zt    ) * rad
+      sh_p = MOD( sh_p,  2*rpi )
+      !
+      ! Inclination of the orbit of the moon w.r.t. the celestial equator, see
+      ! explanation of Table 6 of S58
+      sh_I = ACOS( rcice - rsise * COS( sh_N ) )
+      !
+      ! Computation of sh_xi and sh_nu, see explanation of Table 6 of S58
+      ztgn2 = TAN( sh_N / 2.0_wp )
+      zat1 = ATAN( rxinu1 * ztgn2 )
+      zat2 = ATAN( rxinu2 * ztgn2 )
+      sh_xi = sh_N - zat1 - zat2
+      IF( sh_N > rpi ) sh_xi = sh_xi - 2.0_wp * rpi
+      sh_nu = zat1 - zat2
+      !
+      ! Computation of sh_x1ra, sh_R, sh_nuprim, and sh_nusec used for tidal
+      ! constituents L2, K1, and K2
+      !
+      ! Computation of sh_x1ra and sh_R (Equations 204, 213, and 214 of S58)
+      ztgI2 = tan( sh_I / 2.0_wp )
+      zP1   = sh_p - sh_xi
+      zt2   = ztgI2 * ztgI2
+      sh_x1ra = SQRT( 1.0 - 12.0 * zt2 * COS( 2.0_wp * zP1 ) + 36.0_wp * zt2 * zt2 )
+      zp = SIN( 2.0_wp * zP1 )
+      zq = 1.0_wp / ( 6.0_wp * zt2 ) - COS( 2.0_wp * zP1 )
+      sh_R = ATAN( zp / zq )
+      !
+      ! Computation of sh_nuprim (Equation 224 of S58)
+      zp = SIN( 2.0_wp * sh_I ) * SIN( sh_nu )
+      zq = SIN( 2.0_wp * sh_I ) * COS( sh_nu ) + 0.3347_wp
+      sh_nuprim = ATAN( zp / zq )
+      !
+      ! Computation of sh_nusec  (Equation 232 of S58)
+      zs2 = SIN( sh_I ) * SIN( sh_I )
+      zp  = zs2 * SIN( 2.0_wp * sh_nu )
+      zq  = zs2 * COS( 2.0_wp * sh_nu ) + 0.0727_wp
+      sh_nusec = 0.5_wp * ATAN( zp / zq )
+      !
+   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
+      !!----------------------------------------------------------------------
+      !
+      zscale =  rad / 3600.0_wp
+      !
+      DO jh = 1, size(ptide_harmo)
+         ptide_harmo(jh)%omega = (  romega_T * ptide_comp( jh )%nT   &
+            &                     + romega_s * ptide_comp( jh )%ns   &
+            &                     + romega_h * ptide_comp( jh )%nh   &
+            &                     + romega_p * ptide_comp( jh )%np   &
+            &                     + romega_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 )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION nodal_factort  ***
+      !!
+      !! ** Purpose : Compute amplitude correction factors due to nodal motion
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kformula
+      !
+      REAL(wp)            ::   zf
+      REAL(wp)            ::   zs, zf1, zf2
+      CHARACTER(LEN=3)    ::   clformula
+      !!----------------------------------------------------------------------
+      !
+      SELECT CASE( kformula )
+      !
+      CASE( 0 )                  ! Formula 0, solar waves
+         zf = 1.0
+         !
+      CASE( 1 )                  ! Formula 1, compound waves (78 x 78)
+         zf=nodal_factort( 78 )
+         zf = zf * zf
+         !
+      CASE ( 4 )                 ! Formula 4,  compound waves (78 x 235) 
+         zf1 = nodal_factort( 78 )
+         zf  = nodal_factort(235)
+         zf  = zf1 * zf
+         !
+      CASE( 18 )                 ! Formula 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 )                 ! Formula 73 of S58
+         zs = SIN( sh_I )
+         zf = ( 2.0_wp / 3.0_wp - zs * zs ) / 0.5021_wp
+         !
+      CASE( 74 )                 ! Formula 74 of S58
+         zs = SIN(sh_I)
+         zf = zs * zs / 0.1578_wp
+         !
+      CASE( 75 )                 ! Formula 75 of S58
+         zs = COS( sh_I / 2.0_wp )
+         zf = SIN( sh_I ) * zs * zs / 0.3800_wp
+         !
+      CASE( 76 )                 ! Formula 76 of S58
+         zf = SIN( 2.0_wp * sh_I ) / 0.7214_wp
+         !
+      CASE( 78 )                 ! Formula 78 of S58
+         zs = COS( sh_I/2 )
+         zf = zs * zs * zs * zs / 0.9154_wp
+         !
+      CASE( 149 )                ! Formula 149 of S58
+         zs = COS( sh_I/2 )
+         zf = zs * zs * zs * zs * zs * zs / 0.8758_wp
+         !
+      CASE( 215 )                ! Formula 215 of S58 with typo correction (0.9154 instead of 0.9145)
+         zs = COS( sh_I/2 )
+         zf = zs * zs * zs * zs / 0.9154_wp * sh_x1ra
+         !
+      CASE( 227 )                ! Formula 227 of S58
+         zs = SIN( 2.0_wp * sh_I )
+         zf = SQRT( 0.8965_wp * zs * zs + 0.6001_wp * zs * COS( sh_nu ) + 0.1006_wp )
+         !
+      CASE ( 235 )               ! Formula 235 of S58
+         zs = SIN( sh_I )
+         zf = SQRT( 19.0444_wp * zs * zs * zs * zs + 2.7702_wp * zs * zs * cos( 2.0_wp * sh_nu ) + 0.0981_wp )
+         !
+      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 )
+      !!----------------------------------------------------------------------
+      !!                   ***  FUNCTION dayjul  ***
+      !!
+      !! Purpose :  compute the Julian day
+      !!----------------------------------------------------------------------
+      INTEGER,INTENT(in)    ::   kyr, kmonth, kday
+      !
+      INTEGER,DIMENSION(12) ::   idayt = (/ 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 /)
+      INTEGER,DIMENSION(12) ::   idays
+      INTEGER               ::   inc, ji, zyq
+      REAL(wp)              ::   dayjul
+      !!----------------------------------------------------------------------
+      !
+      idays(1) = 0
+      idays(2) = 31
+      inc = 0.0_wp
+      zyq = MOD( kyr - 1900 , 4 )
+      IF( zyq == 0 ) inc = 1
+      DO ji = 3, 12
+         idays(ji) = idayt(ji) + inc
+      END DO
+      dayjul = REAL( idays(kmonth) + kday, KIND=wp )
+      !
+   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
Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/nemogcm.F90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/nemogcm.F90	(revision 12116)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/nemogcm.F90	(revision 12117)
@@ -46,5 +46,5 @@
    USE closea         ! treatment of closed seas (for ln_closea)
    USE usrdef_nam     ! user defined configuration
-   USE tideini        ! tidal components initialization   (tide_ini routine)
+   USE tide_mod, ONLY : tide_init ! tidal components initialization   (tide_init routine)
    USE bdy_oce,  ONLY : ln_bdy
    USE bdyini         ! open boundary cond. setting       (bdy_init routine)
Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/step.F90	(revision 12116)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/step.F90	(revision 12117)
@@ -112,5 +112,5 @@
       ! Update external forcing (tides, open boundaries, and surface boundary condition (including sea-ice)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_tide    )   CALL sbc_tide( kstp )                   ! update tide potential
+      IF( ln_tide    )   CALL tide_update( kstp )                ! update tide potential
       IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                   ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 
       IF( ln_bdy     )   CALL bdy_dta ( kstp, kt_offset = +1 )   ! update dynamic & tracer data at open boundaries
Index: /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/step_oce.F90
===================================================================
--- /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/step_oce.F90	(revision 12116)
+++ /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/step_oce.F90	(revision 12117)
@@ -19,5 +19,5 @@
    USE sbccpl          ! surface boundary condition: coupled formulation (call send at end of step)
    USE sbcapr          ! surface boundary condition: atmospheric pressure
-   USE sbctide         ! Tide initialisation
+   USE tide_mod, ONLY : ln_tide, tide_update
    USE sbcwave         ! Wave intialisation
 
