Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 3951)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 3953)
@@ -1,16 +1,16 @@
 MODULE diaharm 
-
-#if defined key_diaharm && defined key_tide
-   !!=================================================================================
+   !!======================================================================
    !!                       ***  MODULE  diaharm  ***
    !! Harmonic analysis of tidal constituents 
-   !!=================================================================================
-   !! * Modules used
+   !!======================================================================
+   !! History :  3.1  !  2007  (O. Le Galloudec, J. Chanut)  Original code
+   !!----------------------------------------------------------------------
+#if defined key_diaharm && defined key_tide
+   !!----------------------------------------------------------------------
+   !!   'key_diaharm'
+   !!   'key_tide'
+   !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain
-   USE in_out_manager  ! I/O units
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE ioipsl          ! NetCDF IPSL library
-   USE diadimg         ! To write dimg
    USE phycst
    USE dynspg_oce
@@ -18,5 +18,9 @@
    USE daymod
    USE tide_mod
-   USE iom 
+   USE in_out_manager  ! I/O units
+   USE iom             ! I/0 library
+   USE ioipsl          ! NetCDF IPSL library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE diadimg         ! To write dimg
    USE timing          ! preformance summary
    USE wrk_nemo        ! working arrays
@@ -30,36 +34,32 @@
    INTEGER, PARAMETER :: jpdimsparse  = jpincomax*300*24
 
-   INTEGER ::                            & !! namelist variables
-                         nit000_han = 1, & ! First time step used for harmonic analysis
-                         nitend_han = 1, & ! Last time step used for harmonic analysis
-                         nstep_han  = 1, & ! Time step frequency for harmonic analysis
-                         nb_ana            ! Number of harmonics to analyse
-
-   INTEGER , ALLOCATABLE, DIMENSION(:)       :: name
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp
-   REAL(wp), ALLOCATABLE, DIMENSION(:)       :: ana_freq, vt, ut, ft
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: out_eta, &
-                                                out_u  , &
-                                                out_v
-
-   INTEGER :: ninco, nsparse
-   INTEGER ,       DIMENSION(jpdimsparse)         :: njsparse, nisparse
-   INTEGER , SAVE, DIMENSION(jpincomax)           :: ipos1
-   REAL(wp),       DIMENSION(jpdimsparse)         :: valuesparse
-   REAL(wp),       DIMENSION(jpincomax)           :: ztmp4 , ztmp7
-   REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3 , zpilier
-   REAL(wp), SAVE, DIMENSION(jpincomax)           :: zpivot
-
-   CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   &
-       tname         ! Names of tidal constituents ('M2', 'K1',...)
-
-
-!! * Routine accessibility
-   PUBLIC  dia_harm    ! routine called by step.F90
-
-   !!---------------------------------------------------------------------------------
-   !!  
-   !!---------------------------------------------------------------------------------
-
+   !                            !!!namelist variables
+   INTEGER ::   nit000_han = 1   ! First time step used for harmonic analysis
+   INTEGER ::   nitend_han = 1   ! Last time step used for harmonic analysis
+   INTEGER ::   nstep_han  = 1   ! Time step frequency for harmonic analysis
+   INTEGER ::   nb_ana           ! Number of harmonics to analyse
+
+   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ana_temp
+   REAL(wp), ALLOCATABLE, DIMENSION(:)       ::   ana_freq, ut   , vt   , ft
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   out_eta , out_u, out_v
+
+   INTEGER ::   ninco, nsparse
+   INTEGER ,       DIMENSION(jpdimsparse)         ::   njsparse, nisparse
+   INTEGER , SAVE, DIMENSION(jpincomax)           ::   ipos1
+   REAL(wp),       DIMENSION(jpdimsparse)         ::   valuesparse
+   REAL(wp),       DIMENSION(jpincomax)           ::   ztmp4 , ztmp7
+   REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) ::   ztmp3 , zpilier
+   REAL(wp), SAVE, DIMENSION(jpincomax)           ::   zpivot
+
+   CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   tname   ! Names of tidal constituents ('M2', 'K1',...)
+
+   PUBLIC   dia_harm   ! routine called by step.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.5 , NEMO Consortium (2013)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -67,5 +67,4 @@
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE dia_harm_init  ***
-      !!----------------------------------------------------------------------
       !!         
       !! ** Purpose :   Initialization of tidal harmonic analysis
@@ -73,10 +72,7 @@
       !! ** Method  :   Initialize frequency array and  nodal factor for nit000_han
       !!
-      !! History :
-      !!   9.0  O. Le Galloudec and J. Chanut (Original)
-      !!--------------------------------------------------------------------
-      !! * Local declarations 
+      !!--------------------------------------------------------------------
       INTEGER :: jh, nhan, jk, ji
-
+      !
       NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, tname
       !!----------------------------------------------------------------------
@@ -92,7 +88,6 @@
       tname(:)=''
       !
-      ! Read Namelist nam_diaharm
-      REWIND ( numnam )
-      READ   ( numnam, nam_diaharm )
+      REWIND( numnam )                   ! Read Namelist nam_diaharm
+      READ  ( numnam, nam_diaharm )
       !
       IF(lwp) THEN
@@ -104,21 +99,13 @@
       ! Basic checks on harmonic analysis time window:
       ! ----------------------------------------------
-      IF (nit000 > nit000_han) THEN
-        IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nit000_han must be greater than nit000, stop'
-        IF(lwp) WRITE(numout,*) ' restart capability not implemented'
-        nstop = nstop + 1
-      ENDIF
-      IF (nitend < nitend_han) THEN
-        IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nitend_han must be lower than nitend, stop'
-        IF(lwp) WRITE(numout,*) ' restart capability not implemented'
-        nstop = nstop + 1
-      ENDIF
-
-      IF (MOD(nitend_han-nit000_han+1,nstep_han).NE.0) THEN
-        IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : analysis time span must be a multiple of nstep_han, stop'
-        nstop = nstop + 1
-      END IF
-
-      nb_ana=0
+      IF( nit000 > nit000_han )   CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000',   &
+         &                                       ' restart capability not implemented' )
+      IF( nitend < nitend_han )   CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend',   &
+         &                                       'restart capability not implemented' )
+
+      IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 )   &
+         &                        CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' )
+
+      nb_ana = 0
       DO jk=1,jpmax_harmo
          DO ji=1,jpmax_harmo
@@ -153,10 +140,7 @@
       ! Initialize frequency array:
       ! ---------------------------
-      ALLOCATE(ana_freq(nb_ana))
-      ALLOCATE(vt      (nb_ana))
-      ALLOCATE(ut      (nb_ana))
-      ALLOCATE(ft      (nb_ana))
-
-      CALL tide_harmo(ana_freq, vt, ut , ft, name ,nb_ana)
+      ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) )
+
+      CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana )
 
       IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
@@ -168,13 +152,13 @@
       ! Initialize temporary arrays:
       ! ----------------------------
-      ALLOCATE( ana_temp(jpi,jpj,nb_ana*2,3))
+      ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) )
       ana_temp(:,:,:,:) = 0.e0
 
    END SUBROUTINE dia_harm_init
-   
+
+
    SUBROUTINE dia_harm ( kt )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE dia_harm  ***
-      !!----------------------------------------------------------------------
       !!         
       !! ** Purpose :   Tidal harmonic analysis main routine
@@ -182,11 +166,7 @@
       !! ** Action  :   Sums ssh/u/v over time analysis [nit000_han,nitend_han]
       !!
-      !! History :
-      !!   9.0  O. Le Galloudec and J. Chanut (Original)
-      !!--------------------------------------------------------------------
-      !! * Argument:
+      !!--------------------------------------------------------------------
       INTEGER, INTENT( IN ) :: kt
-
-      !! * Local declarations
+      !
       INTEGER  :: ji, jj, jh, jc, nhc
       REAL(wp) :: ztime, ztemp
@@ -194,5 +174,5 @@
       IF( nn_timing == 1 )   CALL timing_start('dia_harm')
 
-      IF ( kt .EQ. nit000 ) CALL dia_harm_init
+      IF ( kt == nit000 ) CALL dia_harm_init
 
       IF ( ((kt.GE.nit000_han).AND.(kt.LE.nitend_han)).AND. &
@@ -211,13 +191,8 @@
               DO ji = 1,jpi
                 ! Elevation
-                ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1)                &
-                                        + ztemp*sshn(ji,jj)*tmask(ji,jj,1)        
+                ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)           *tmask(ji,jj,1)        
 #if defined key_dynspg_ts
-                ! ubar
-                ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2)                &
-                                        + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1)
-                ! vbar
-                ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3)                &
-                                        + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1)
+                ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1)
+                ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1)
 #endif
               END DO
@@ -229,5 +204,5 @@
       END IF
 
-      IF ( kt .EQ. nitend_han ) CALL dia_harm_end
+      IF ( kt == nitend_han )   CALL dia_harm_end
 
       IF( nn_timing == 1 )   CALL timing_stop('dia_harm')
@@ -235,8 +210,8 @@
    END SUBROUTINE dia_harm
 
+
    SUBROUTINE dia_harm_end
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE diaharm_end  ***
-      !!----------------------------------------------------------------------
       !!         
       !! ** Purpose :  Compute the Real and Imaginary part of tidal constituents
@@ -244,9 +219,5 @@
       !! ** Action  :  Decompose the signal on the harmonic constituents 
       !!
-      !! History :
-      !!   9.0  O. Le Galloudec and J. Chanut (Original)
-      !!--------------------------------------------------------------------
-
-      !! * Local declarations
+      !!--------------------------------------------------------------------
       INTEGER :: ji, jj, jh, jc, jn, nhan, jl
       INTEGER :: ksp, kun, keq
@@ -279,12 +250,11 @@
                nisparse(ksp) = keq
                njsparse(ksp) = kun
-               valuesparse(ksp)= &
-                   +(   MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &
-                   +(1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
+               valuesparse(ksp) = (   MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))   &
+                  &             + (1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)) )
             END DO
          END DO
       END DO
 
-      nsparse=ksp
+      nsparse = ksp
 
       ! Elevation:
@@ -292,11 +262,11 @@
          DO ji = 1, jpi
             ! Fill input array
-            kun=0
-            DO jh = 1,nb_ana
-               DO jc = 1,2
+            kun = 0
+            DO jh = 1, nb_ana
+               DO jc = 1, 2
                   kun = kun + 1
                   ztmp4(kun)=ana_temp(ji,jj,kun,1)
-               ENDDO
-            ENDDO
+               END DO
+            END DO
 
             CALL SUR_DETERMINE(jj)
@@ -310,15 +280,15 @@
       END DO
 
-      ALLOCATE(out_eta(jpi,jpj,2*nb_ana))
-      ALLOCATE(out_u  (jpi,jpj,2*nb_ana))
-      ALLOCATE(out_v  (jpi,jpj,2*nb_ana))
+      ALLOCATE( out_eta(jpi,jpj,2*nb_ana),   & 
+         &      out_u  (jpi,jpj,2*nb_ana),   &
+         &      out_v  (jpi,jpj,2*nb_ana)  )
 
       DO jj = 1, jpj
          DO ji = 1, jpi
             DO jh = 1, nb_ana 
-               X1=ana_amp(ji,jj,jh,1)
-               X2=-ana_amp(ji,jj,jh,2)
-               out_eta(ji,jj,jh)=X1 * tmask(ji,jj,1)
-               out_eta(ji,jj,nb_ana+jh)=X2 * tmask(ji,jj,1)
+               X1 = ana_amp(ji,jj,jh,1)
+               X2 =-ana_amp(ji,jj,jh,2)
+               out_eta(ji,jj,jh       ) = X1 * tmask(ji,jj,1)
+               out_eta(ji,jj,jh+nb_ana) = X2 * tmask(ji,jj,1)
             ENDDO
          ENDDO
@@ -398,17 +368,11 @@
    END SUBROUTINE dia_harm_end
 
+
    SUBROUTINE dia_wri_harm
       !!--------------------------------------------------------------------
       !!                 ***  ROUTINE dia_wri_harm  ***
-      !!--------------------------------------------------------------------
       !!         
       !! ** Purpose : Write tidal harmonic analysis results in a netcdf file
-      !!
-      !!
-      !! History :
-      !!   9.0  O. Le Galloudec and J. Chanut (Original)
-      !!--------------------------------------------------------------------
-
-      !! * Local declarations
+      !!--------------------------------------------------------------------
       CHARACTER(LEN=lc) :: cltext
       CHARACTER(LEN=lc) ::   &
@@ -468,10 +432,11 @@
 #else
       DO jh = 1, nb_ana
-      CALL iom_put( TRIM(tname(jh))//'x_v', out_u(:,:,jh) )
-      CALL iom_put( TRIM(tname(jh))//'y_v', out_u(:,:,nb_ana+jh) )
+         CALL iom_put( TRIM(tname(jh))//'x_v', out_u(:,:,jh       ) )
+         CALL iom_put( TRIM(tname(jh))//'y_v', out_u(:,:,jh+nb_ana) )
       END DO
 #endif
 
    END SUBROUTINE dia_wri_harm
+
 
    SUBROUTINE SUR_DETERMINE(init)
@@ -482,6 +447,6 @@
    !!       
    !!---------------------------------------------------------------------------------
-   INTEGER, INTENT(in) :: init 
-               
+   INTEGER, INTENT(in) ::   init 
+   !
    INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
    REAL(wp)                        :: zval1, zval2, zx1
@@ -492,27 +457,20 @@
    CALL wrk_alloc( jpincomax , ipos2 , ipivot        )
             
-   IF( init==1 )THEN
-
-      IF( nsparse .GT. jpdimsparse ) &
-         CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
-
-      IF( ninco .GT. jpincomax ) &
-         CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
-
-      ztmp3(:,:)=0.e0
-
+   IF( init == 1 ) THEN
+      IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
+      IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
+      !
+      ztmp3(:,:) = 0._wp
+      !
       DO jk1_sd = 1, nsparse
          DO jk2_sd = 1, nsparse
-
-            nisparse(jk2_sd)=nisparse(jk2_sd)
-            njsparse(jk2_sd)=njsparse(jk2_sd)
-
+            nisparse(jk2_sd) = nisparse(jk2_sd)
+            njsparse(jk2_sd) = njsparse(jk2_sd)
             IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
                ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  &
                                                         + valuesparse(jk1_sd)*valuesparse(jk2_sd)
             ENDIF
-
-         ENDDO
-      ENDDO
+         END DO
+      END DO
 
       DO jj_sd = 1 ,ninco
@@ -584,5 +542,4 @@
    ENDDO
 
-
    CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 )
    CALL wrk_dealloc( jpincomax , ipos2 , ipivot        )
@@ -590,5 +547,4 @@
   END SUBROUTINE SUR_DETERMINE
 
-
 #else
    !!----------------------------------------------------------------------
@@ -597,12 +553,10 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_diaharm = .FALSE.
 CONTAINS
-
    SUBROUTINE dia_harm ( kt )     ! Empty routine
       INTEGER, INTENT( IN ) :: kt  
       WRITE(*,*) 'dia_harm: you should not have seen this print'
    END SUBROUTINE dia_harm
-
-
-#endif
+#endif
+
    !!======================================================================
 END MODULE diaharm
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 3951)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 3953)
@@ -22,4 +22,6 @@
    USE dynspg_flt     ! surface pressure gradient     (dyn_spg_flt routine)
    USE dynadv         ! dynamics: vector invariant versus flux form
+   USE sbctide
+   USE updtide
    USE trdmod         ! ocean dynamics trends
    USE trdmod_oce     ! ocean variables trends
@@ -100,15 +102,57 @@
       ENDIF
 
-      IF( ln_apr_dyn ) THEN                   !==  Atmospheric pressure gradient  ==!
-         zg_2 = grav * 0.5
-         DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh
+      IF(      ln_apr_dyn                                                &   ! atmos. pressure
+         .OR.  ( .NOT.lk_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) )   &   ! tide potential (no time slitting)
+         .OR.  nn_ice_embd == 2  ) THEN                                      ! embedded sea-ice
+         !
+         DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               spgu(ji,jj) =  zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
-                  &                   + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) /e1u(ji,jj)
-               spgv(ji,jj) =  zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    &
-                  &                   + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj)
-            END DO
-         END DO
-         DO jk = 1, jpkm1                          ! Add the apg to the general trend
+               spgu(ji,jj) = 0._wp
+               spgv(ji,jj) = 0._wp
+            END DO
+         END DO         
+         !
+         IF( ln_apr_dyn ) THEN                !==  Atmospheric pressure gradient  ==!
+            zg_2 = grav * 0.5
+            DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
+                     &                      + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) /e1u(ji,jj)
+                  spgv(ji,jj) = spgv(ji,jj) + zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    &
+                     &                      + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj)
+               END DO
+            END DO
+         ENDIF
+         !
+         !                                    !==  tide potential forcing term  ==!
+         IF( .NOT.lk_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case
+            !
+            CALL upd_tide( kt )                      ! update tide potential
+            !
+            DO jj = 2, jpjm1                         ! add tide potential forcing
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  spgv(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
+                  spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
+               END DO 
+            END DO
+         ENDIF
+         !
+         IF( nn_ice_embd == 2 ) THEN          !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
+            CALL wrk_alloc( jpi, jpj, zpice )
+            !                                            
+            zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
+            zgrau0r     = - grav * r1_rau0
+            zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj)
+                  spgv(ji,jj) = spgu(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj)
+               END DO
+            END DO
+            !
+            CALL wrk_dealloc( jpi, jpj, zpice )         
+         ENDIF
+         !
+         DO jk = 1, jpkm1                     !== Add all terms to the general trend
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -117,31 +161,6 @@
                END DO
             END DO
-         END DO
-      ENDIF
-
-      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
-         CALL wrk_alloc( jpi, jpj, zpice )
-         !                                            
-         zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
-         zgrau0r     = - grav * r1_rau0
-         zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               spgu(ji,jj) = ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj)
-               spgv(ji,jj) = ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj)
-            END DO
-         END DO
-         DO jk = 1, jpkm1                             ! Add the surface pressure trend to the general trend
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
-                  va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
-               END DO
-            END DO
-         END DO
-         !
-         CALL wrk_dealloc( jpi, jpj, zpice )
-      ENDIF
-
+         END DO         
+      ENDIF
 
       SELECT CASE ( nspg )                       ! compute surf. pressure gradient trend and add it to the general trend
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90	(revision 3951)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90	(revision 3953)
@@ -91,5 +91,6 @@
                spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
             END DO 
-         END DO 
+         END DO
+         !
          DO jk = 1, jpkm1                    ! Add it to the general trend
             DO jj = 2, jpjm1
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 3951)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 3953)
@@ -396,7 +396,7 @@
          !                                                !* Update the forcing (BDY and tides)
          !                                                !  ------------------
-         IF( lk_obc )   CALL obc_dta_bt ( kt, jn   )
-         IF( lk_bdy )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 )
-         IF ( ln_tide_pot .AND. lk_tide) CALL upd_tide( kt, jn )
+         IF( lk_obc                    )   CALL obc_dta_bt( kt, jn   )
+         IF( lk_bdy                    )   CALL bdy_dta   ( kt, jit=jn, time_offset=1 )
+         IF( ln_tide_pot .AND. lk_tide )   CALL upd_tide  ( kt, kit=jn, kbaro=nn_baro )
 
          !                                                !* after ssh_e
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90	(revision 3951)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90	(revision 3953)
@@ -1,107 +1,115 @@
 MODULE sbctide
-  !!=================================================================================
-  !!                       ***  MODULE  sbctide  ***
-  !! Initialization of tidal forcing
-  !! History :  9.0  !  07  (O. Le Galloudec)  Original code
-  !!=================================================================================
-  !! * Modules used
-  USE oce             ! ocean dynamics and tracers variables
-  USE dom_oce         ! ocean space and time domain
-  USE in_out_manager  ! I/O units
-  USE ioipsl          ! NetCDF IPSL library
-  USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-  USE phycst
-  USE daymod
-  USE dynspg_oce
-  USE tideini
-  USE iom
+   !!======================================================================
+   !!                       ***  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
+   USE daymod
+   USE dynspg_oce
+   USE tideini
+   !
+   USE iom
+   USE in_out_manager  ! I/O units
+   USE ioipsl          ! NetCDF IPSL library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
 
-  IMPLICIT NONE
-  PUBLIC
+   IMPLICIT NONE
+   PUBLIC
 
-  REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: pot_astro
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   pot_astro   !
 
 #if defined key_tide
+   !!----------------------------------------------------------------------
+   !!   'key_tide' :                                        tidal potential
+   !!----------------------------------------------------------------------
+   !!   sbc_tide            : 
+   !!   tide_init_potential :
+   !!----------------------------------------------------------------------
 
-  LOGICAL, PUBLIC, PARAMETER ::   lk_tide  = .TRUE.
-  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: amp_pot,phi_pot
-  !!---------------------------------------------------------------------------------
-  !!   OPA 9.0 , LODYC-IPSL  (2003)
-  !!---------------------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_tide  = .TRUE.
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   amp_pot, phi_pot
 
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.5 , NEMO Consortium (2013)
+   !! $Id: $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
-  SUBROUTINE sbc_tide ( kt )
-    !!----------------------------------------------------------------------
-    !!                 ***  ROUTINE sbc_tide  ***
-    !!----------------------------------------------------------------------      
-    !! * Arguments
-    INTEGER, INTENT( in ) ::   kt     ! ocean time-step
-    !!----------------------------------------------------------------------
+   SUBROUTINE sbc_tide( kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE sbc_tide  ***
+      !!----------------------------------------------------------------------      
+      INTEGER, INTENT( in ) ::   kt     ! ocean time-step
+      !!----------------------------------------------------------------------
 
-    IF ( kt == nit000 .AND. .NOT. lk_dynspg_ts )  CALL ctl_stop( 'STOP', 'sbc_tide : tidal potential use only with time splitting' )
+      IF( kt /= nit000 )   CALL tide_init( kt )
 
-    IF ( nsec_day == NINT(0.5 * rdttra(1)) ) THEN
+      IF( nsec_day == NINT(0.5 * rdttra(1)) ) 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)   )
+            !
+            amp_pot(:,:,:) = 0._wp
+            phi_pot(:,:,:) = 0._wp
+            pot_astro(:,:) = 0._wp
+         ENDIF
+         !
+         CALL tide_harmo( omega_tide, v0tide, utide, ftide, ntide, nb_harmo )
+         !
+         kt_tide = kt
+         !
+         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
+         !
+      ENDIF
       !
-      kt_tide = kt
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'sbc_tide : (re)Initialization of the tidal potential at kt=',kt
-         WRITE(numout,*) '~~~~~~~ '
-      ENDIF
-
-      IF(lwp) THEN
-         IF ( kt == nit000 ) WRITE(numout,*) 'Apply astronomical potential : ln_tide_pot =', ln_tide_pot
-         CALL flush(numout)
-      ENDIF
-
-      IF ( kt == nit000 ) ALLOCATE(amp_pot(jpi,jpj,nb_harmo))
-      IF ( kt == nit000 ) ALLOCATE(phi_pot(jpi,jpj,nb_harmo))
-      IF ( kt == nit000 ) ALLOCATE(pot_astro(jpi,jpj))
-
-      amp_pot(:,:,:) = 0.e0
-      phi_pot(:,:,:) = 0.e0
-      pot_astro(:,:) = 0.e0
-
-      IF ( ln_tide_pot ) CALL tide_init_potential
-      !
-    ENDIF
-
-  END SUBROUTINE sbc_tide
-
-  SUBROUTINE tide_init_potential
-    !!----------------------------------------------------------------------
-    !!                 ***  ROUTINE tide_init_potential  ***
-    !!----------------------------------------------------------------------
-    !! * Local declarations
-    INTEGER  :: ji,jj,jk
-    REAL(wp) :: zcons,ztmp1,ztmp2,zlat,zlon
+   END SUBROUTINE sbc_tide
 
 
-    DO jk=1,nb_harmo
-       zcons=0.7*Wave(ntide(jk))%equitide*ftide(jk)
-       do ji=1,jpi
-          do jj=1,jpj
-             ztmp1 = amp_pot(ji,jj,jk)*COS(phi_pot(ji,jj,jk))
-             ztmp2 = -amp_pot(ji,jj,jk)*SIN(phi_pot(ji,jj,jk))
-             zlat = gphit(ji,jj)*rad !! latitude en radian
-             zlon = glamt(ji,jj)*rad !! longitude en radian
-             ! le potentiel est composé des effets des astres:
-             IF (Wave(ntide(jk))%nutide .EQ.1) THEN
-                ztmp1= ztmp1 + zcons*(SIN(2.*zlat))*COS(v0tide(jk)+utide(jk)+Wave(ntide(jk))%nutide*zlon)
-                ztmp2= ztmp2 - zcons*(SIN(2.*zlat))*SIN(v0tide(jk)+utide(jk)+Wave(ntide(jk))%nutide*zlon)
-             ENDIF
-             IF (Wave(ntide(jk))%nutide.EQ.2) THEN
-                ztmp1= ztmp1 + zcons*(COS(zlat)**2)*COS(v0tide(jk)+utide(jk)+Wave(ntide(jk))%nutide*zlon)
-                ztmp2= ztmp2 - zcons*(COS(zlat)**2)*SIN(v0tide(jk)+utide(jk)+Wave(ntide(jk))%nutide*zlon)
-             ENDIF
-             amp_pot(ji,jj,jk)=SQRT(ztmp1**2+ztmp2**2)
-             phi_pot(ji,jj,jk)=ATAN2(-ztmp2/MAX(1.E-10,SQRT(ztmp1**2+ztmp2**2)),ztmp1/MAX(1.E-10,SQRT(ztmp1**2+ztmp2**2)))
-          enddo
-       enddo
-    END DO
+   SUBROUTINE tide_init_potential
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE tide_init_potential  ***
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zcons, ztmp1, ztmp2, zlat, zlon, ztmp, zamp, zsc   ! local scalar
+      !!----------------------------------------------------------------------
 
-  END SUBROUTINE tide_init_potential
+      DO jk = 1, nb_harmo
+         zcons = 0.7 * Wave(ntide(jk))%equitide * ftide(jk)
+         DO ji = 1, jpi
+            DO jj = 1, jpj
+               ztmp1 =  amp_pot(ji,jj,jk) * COS( phi_pot(ji,jj,jk) )
+               ztmp2 = -amp_pot(ji,jj,jk) * SIN( phi_pot(ji,jj,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 )   zcs = zcons * SIN( 2.*zlat )
+               IF( Wave(ntide(jk))%nutide == 2 )   zcs = zcons * COS( zlat )**2
+               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 , zamp ) ,   &
+                  &                        ztmp1 / MAX( 1.e-10,  zamp )   )
+            END DO
+         END DO
+      END DO
+      !
+   END SUBROUTINE tide_init_potential
 
 #else
@@ -116,5 +124,5 @@
   END SUBROUTINE sbc_tide
 #endif
+
   !!======================================================================
-
 END MODULE sbctide
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tide.h90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tide.h90	(revision 3951)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tide.h90	(revision 3953)
@@ -1,28 +1,29 @@
-!! History :  9.0  !  07  (O. Le Galloudec)  Original code
+   !!----------------------------------------------------------------------
+   !! History :  3.2  !  2007  (O. Le Galloudec)  Original code
+   !!----------------------------------------------------------------------
 
-! Wave(1)= tide(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)
+   !             !! 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: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90	(revision 3951)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90	(revision 3953)
@@ -1,493 +1,411 @@
 MODULE tide_mod
-  !!=================================================================================
-  !!                       ***  MODULE  tide_mod  ***
-  !! Compute nodal modulations corrections and pulsations
-  !!=================================================================================
-  !!---------------------------------------------------------------------------------
-  !!   OPA 9.0 , LODYC-IPSL  (2003)
-  !!---------------------------------------------------------------------------------
-  USE dom_oce         ! ocean space and time domain
-  USE phycst
-  USE daymod
-
-  IMPLICIT NONE
-  PRIVATE
-
-  REAL(wp) :: sh_T, sh_s, sh_h, sh_p, sh_p1, &
-       sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R, &
-       sh_I, sh_x1ra, sh_N
-
-  INTEGER,PUBLIC, PARAMETER ::   &
-       jpmax_harmo = 19             ! maximum number of harmonic
-
-  TYPE, PUBLIC ::    tide
-     CHARACTER(LEN=4)  :: cname_tide
-     REAL(wp) :: equitide
-     INTEGER  :: nutide
-     INTEGER  ::  nt,ns,nh,np,np1,shift
-     INTEGER  ::  nksi,nnu0,nnu1,nnu2,R
-     INTEGER  :: nformula
-  END TYPE tide
-
-  TYPE(tide), PUBLIC, DIMENSION(jpmax_harmo) :: Wave
-
-  !! * Accessibility
-  PUBLIC tide_harmo
-  PUBLIC nodal_factort
-  PUBLIC tide_init_Wave
-
+   !!======================================================================
+   !!                       ***  MODULE  tide_mod  ***
+   !! Compute nodal modulations corrections and pulsations
+   !!======================================================================
+   !! History :  1.0  !  2007  (O. Le Galloudec)  Original code
+   !!----------------------------------------------------------------------
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constant
+   USE daymod         ! calendar
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   tide_harmo       ! called by tideini and diaharm modules
+   PUBLIC   tide_init_Wave   ! called by tideini and diaharm modules
+
+   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 19   !: maximum number of harmonic
+
+   TYPE, PUBLIC ::    tide
+      CHARACTER(LEN=4) ::   cname_tide
+      REAL(wp)         ::   equitide
+      INTEGER          ::   nutide
+      INTEGER          ::   nt, ns, nh, np, np1, shift
+      INTEGER          ::   nksi, nnu0, nnu1, nnu2, R
+      INTEGER          ::   nformula
+   END TYPE tide
+
+   TYPE(tide), PUBLIC, DIMENSION(jpmax_harmo) ::   Wave   !:
+
+   REAL(wp) ::   sh_T, sh_s, sh_h, sh_p, sh_p1             ! astronomic angles
+   REAL(wp) ::   sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R   !
+   REAL(wp) ::   sh_I, sh_x1ra, sh_N                       !
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
-  SUBROUTINE tide_init_Wave
-
-#  include "tide.h90"
-
-  END SUBROUTINE tide_init_Wave
-
-  SUBROUTINE tide_harmo( pomega, pvt, put , pcor, ktide ,kc)
-
-    INTEGER, DIMENSION(kc), INTENT( in ) ::   &
-         ktide      ! Indice of tidal constituents
-
-    INTEGER, INTENT( in ) :: &
-         kc         ! Total number of tidal constituents
-
-    REAL (wp), DIMENSION(kc), INTENT( out ) ::   &
-         pomega      ! pulsation in radians/s
-
-    REAL (wp), DIMENSION(kc), INTENT( out ) ::   &
-         pvt, &      !
-         put, &      !
-         pcor         !
-
-    CALL astronomic_angle
-    CALL tide_pulse(pomega, ktide ,kc)
-    CALL tide_vuf( pvt, put, pcor, ktide ,kc)
-
-  END SUBROUTINE tide_harmo
-
-  SUBROUTINE astronomic_angle
-
-    !!----------------------------------------------------------------------
-    !!
-    !!  tj is time elapsed since 1st January 1900, 0 hour, counted in julian
-    !!  century (e.g. time in days divide by 36525)
-    !!----------------------------------------------------------------------
-
-    REAL(wp) ::  cosI,p,q,t2,t4,sin2I,s2,tgI2,P1,sh_tgn2,at1,at2
-    REAL(wp) :: zqy,zsy,zday,zdj,zhfrac
-
-    zqy=AINT((nyear-1901.)/4.)
-    zsy=nyear-1900.
-
-    zdj=dayjul(nyear,nmonth,nday)
-    zday=zdj+zqy-1.
-
-    zhfrac=nsec_day/3600.
-
-    !----------------------------------------------------------------------
-    !  Sh_n Longitude of ascending lunar node
-    !----------------------------------------------------------------------
-
-    sh_N=(259.1560564-19.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*rad
-    !----------------------------------------------------------------------
-    ! T mean solar angle (Greenwhich time)
-    !----------------------------------------------------------------------
-    sh_T=(180.+zhfrac*(360./24.))*rad
-    !----------------------------------------------------------------------
-    ! h mean solar Longitude
-    !----------------------------------------------------------------------
-
-    sh_h=(280.1895014-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad
-    !----------------------------------------------------------------------
-    ! s mean lunar Longitude
-    !----------------------------------------------------------------------
-
-    sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad
-    !----------------------------------------------------------------------
-    ! p1 Longitude of solar perigee
-    !----------------------------------------------------------------------
-
-    sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad
-    !----------------------------------------------------------------------
-    ! p Longitude of lunar perigee
-    !----------------------------------------------------------------------
-
-    sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad
-
-    sh_N =mod(sh_N ,2*rpi)
-    sh_s =mod(sh_s ,2*rpi)
-    sh_h =mod(sh_h, 2*rpi)
-    sh_p =mod(sh_p, 2*rpi)
-    sh_p1=mod(sh_p1,2*rpi)
-
-    cosI=0.913694997 -0.035692561 *cos(sh_N)
-
-    sh_I=acos(cosI)
-
-    sin2I=sin(sh_I)
-    sh_tgn2=tan(sh_N/2.0)
-
-    at1=atan(1.01883*sh_tgn2)
-    at2=atan(0.64412*sh_tgn2)
-
-    sh_xi=-at1-at2+sh_N
-
-    if (sh_N > rpi) sh_xi=sh_xi-2.0*rpi
-
-    sh_nu=at1-at2
-
-    !----------------------------------------------------------------------
-    ! For constituents l2 k1 k2
-    !----------------------------------------------------------------------
-
-    tgI2=tan(sh_I/2.0)
-    P1=sh_p-sh_xi
-
-    t2=tgI2*tgI2
-    t4=t2*t2
-    sh_x1ra=sqrt(1.0-12.0*t2*cos(2.0*P1)+36.0*t4)
-
-    p=sin(2.0*P1)
-    q=1.0/(6.0*t2)-cos(2.0*P1)
-    sh_R=atan(p/q)
-
-    p=sin(2.0*sh_I)*sin(sh_nu)
-    q=sin(2.0*sh_I)*cos(sh_nu)+0.3347
-    sh_nuprim=atan(p/q)
-
-    s2=sin(sh_I)*sin(sh_I)
-    p=s2*sin(2.0*sh_nu)
-    q=s2*cos(2.0*sh_nu)+0.0727
-    sh_nusec=0.5*atan(p/q)
-
-  END SUBROUTINE astronomic_angle
-
-  SUBROUTINE tide_pulse( pomega, ktide ,kc)
-    !!----------------------------------------------------------------------
-    !!                     ***  ROUTINE tide_pulse  ***
-    !!                      
-    !! ** Purpose : Compute tidal frequencies
-    !!
-    !!----------------------------------------------------------------------
-    !! * Arguments
-    INTEGER, DIMENSION(kc), INTENT( in ) ::   &
-         ktide      ! Indice of tidal constituents
-
-    INTEGER, INTENT( in ) :: &
-         kc         ! Total number of tidal constituents
-
-    REAL (wp), DIMENSION(kc), INTENT( out ) ::   &
-         pomega      ! pulsation in radians/s
-
-    !! * Local declarations
-    INTEGER :: jh
-    REAL(wp) :: zscale  =  36525*24.0
-    REAL(wp) :: zomega_T=  13149000.0
-    REAL(wp) :: zomega_s=    481267.892
-    REAL(wp) :: zomega_h=     36000.76892
-    REAL(wp) :: zomega_p=      4069.0322056
-    REAL(wp) :: zomega_n=      1934.1423972
-    REAL(wp) :: zomega_p1=        1.719175
-    !!----------------------------------------------------------------------
-
-    DO jh=1,kc
-       pomega(jh) = zomega_T * Wave(ktide(jh))%nT &
-            + zomega_s * Wave(ktide(jh))%ns &
-            + zomega_h * Wave(ktide(jh))%nh &
-            + zomega_p * Wave(ktide(jh))%np &
-            + zomega_p1* Wave(ktide(jh))%np1
-       pomega(jh) = (pomega(jh)/zscale)*rad/3600.
-    END DO
-
-  END SUBROUTINE tide_pulse
-
-  SUBROUTINE tide_vuf( pvt, put, pcor, ktide ,kc)
-    !!----------------------------------------------------------------------
-    !!                     ***  ROUTINE tide_vuf  ***
-    !!                      
-    !! ** Purpose : Compute nodal modulation corrections
-    !!
-    !! ** Outputs :
-    !!          vt: Pase of tidal potential relative to Greenwich (radians)
-    !!          ut: Phase correction u due to nodal motion (radians)
-    !!          ft: Nodal correction factor
-    !!
-    !! ** Inputs :
-    !!          tname: array of constituents names (dimension<=nc) 
-    !!             nc: number of constituents
-    !!   
-    !!----------------------------------------------------------------------
-    !! * Arguments
-    INTEGER, DIMENSION(kc), INTENT( in ) ::   &
-         ktide      ! Indice of tidal constituents
-    INTEGER, INTENT( in ) :: &
-         kc         ! Total number of tidal constituents
-    REAL (wp), DIMENSION(kc), INTENT( out ) ::   &
-         pvt, &      !
-         put, &      !
-         pcor         !
-    !! * Local declarations
-    INTEGER :: jh
-    !!----------------------------------------------------------------------
-
-    DO jh =1,kc
-       !  Phase of the tidal potential relative to the Greenwhich 
-       !  meridian (e.g. the position of the fictuous celestial body). Units are
-       !  radian:
-       pvt(jh) = sh_T *Wave(ktide(jh))%nT        &
-            +sh_s *Wave(ktide(jh))%ns        &
-            +sh_h *Wave(ktide(jh))%nh        &
-            +sh_p *Wave(ktide(jh))%np        &
-            +sh_p1*Wave(ktide(jh))%np1       &
-            +Wave(ktide(jh))%shift*rad
+   SUBROUTINE tide_init_Wave
+#     include "tide.h90"
+   END SUBROUTINE tide_init_Wave
+
+
+   SUBROUTINE tide_harmo( pomega, pvt, put , pcor, ktide ,kc)
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide            ! Indice of tidal constituents
+      INTEGER                , INTENT(in ) ::   kc               ! Total number of tidal constituents
+      REAL(wp), DIMENSION(kc), INTENT(out) ::   pomega           ! pulsation in radians/s
+      REAL(wp), DIMENSION(kc), INTENT(out) ::   pvt, put, pcor   !
+      !!----------------------------------------------------------------------
+      !
+      CALL astronomic_angle
+      CALL tide_pulse( pomega, ktide ,kc )
+      CALL tide_vuf  ( pvt, put, pcor, ktide ,kc )
+      !
+   END SUBROUTINE tide_harmo
+
+
+   SUBROUTINE astronomic_angle
+      !!----------------------------------------------------------------------
+      !!  tj is time elapsed since 1st January 1900, 0 hour, counted in julian
+      !!  century (e.g. time in days divide by 36525)
+      !!----------------------------------------------------------------------
+      REAL(wp) ::   cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2
+      REAL(wp) ::   zqy , zsy, zday, zdj, zhfrac
+      !!----------------------------------------------------------------------
+      !
+      zqy = AINT( (nyear-1901.)/4. )
+      zsy = nyear - 1900.
+      !
+      zdj  = dayjul( nyear, nmonth, nday )
+      zday = zdj + zqy - 1.
+      !
+      zhfrac = nsec_day / 3600.
+      !
+      !----------------------------------------------------------------------
+      !  Sh_n Longitude of ascending lunar node
+      !----------------------------------------------------------------------
+      sh_N=(259.1560564-19.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*rad
+      !----------------------------------------------------------------------
+      ! T mean solar angle (Greenwhich time)
+      !----------------------------------------------------------------------
+      sh_T=(180.+zhfrac*(360./24.))*rad
+      !----------------------------------------------------------------------
+      ! h mean solar Longitude
+      !----------------------------------------------------------------------
+      sh_h=(280.1895014-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad
+      !----------------------------------------------------------------------
+      ! s mean lunar Longitude
+      !----------------------------------------------------------------------
+      sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad
+      !----------------------------------------------------------------------
+      ! p1 Longitude of solar perigee
+      !----------------------------------------------------------------------
+      sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad
+      !----------------------------------------------------------------------
+      ! p Longitude of lunar perigee
+      !----------------------------------------------------------------------
+      sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad
+
+      sh_N = MOD( sh_N ,2*rpi )
+      sh_s = MOD( sh_s ,2*rpi )
+      sh_h = MOD( sh_h, 2*rpi )
+      sh_p = MOD( sh_p, 2*rpi )
+      sh_p1= MOD( sh_p1,2*rpi )
+
+      cosI = 0.913694997 -0.035692561 *cos(sh_N)
+
+      sh_I = ACOS( cosI )
+
+      sin2I   = sin(sh_I)
+      sh_tgn2 = tan(sh_N/2.0)
+
+      at1=atan(1.01883*sh_tgn2)
+      at2=atan(0.64412*sh_tgn2)
+
+      sh_xi=-at1-at2+sh_N
+
+      IF( sh_N > rpi )   sh_xi=sh_xi-2.0*rpi
+
+      sh_nu = at1 - at2
+
+      !----------------------------------------------------------------------
+      ! For constituents l2 k1 k2
+      !----------------------------------------------------------------------
+
+      tgI2 = tan(sh_I/2.0)
+      P1   = sh_p-sh_xi
+
+      t2 = tgI2*tgI2
+      t4 = t2*t2
+      sh_x1ra = sqrt( 1.0-12.0*t2*cos(2.0*P1)+36.0*t4 )
+
+      p = sin(2.0*P1)
+      q = 1.0/(6.0*t2)-cos(2.0*P1)
+      sh_R = atan(p/q)
+
+      p = sin(2.0*sh_I)*sin(sh_nu)
+      q = sin(2.0*sh_I)*cos(sh_nu)+0.3347
+      sh_nuprim = atan(p/q)
+
+      s2 = sin(sh_I)*sin(sh_I)
+      p  = s2*sin(2.0*sh_nu)
+      q  = s2*cos(2.0*sh_nu)+0.0727
+      sh_nusec = 0.5*atan(p/q)
+      !
+   END SUBROUTINE astronomic_angle
+
+
+   SUBROUTINE tide_pulse( pomega, ktide ,kc )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE tide_pulse  ***
+      !!                      
+      !! ** Purpose : Compute tidal frequencies
+      !!----------------------------------------------------------------------
+      INTEGER                , INTENT(in ) ::   kc       ! Total number of tidal constituents
+      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide    ! Indice of tidal constituents
+      REAL(wp), DIMENSION(kc), INTENT(out) ::   pomega   ! pulsation in radians/s
+      !
+      INTEGER  ::   jh
+      REAL(wp) ::   zscale
+      REAL(wp) ::   zomega_T =  13149000.0_wp
+      REAL(wp) ::   zomega_s =    481267.892_wp
+      REAL(wp) ::   zomega_h =     36000.76892_wp
+      REAL(wp) ::   zomega_p =      4069.0322056_wp
+      REAL(wp) ::   zomega_n =      1934.1423972_wp
+      REAL(wp) ::   zomega_p1=         1.719175_wp
+      !!----------------------------------------------------------------------
+      !
+      zscale =  rad / ( 36525._wp * 86400._wp ) 
+      !
+      DO jh = 1, kc
+         pomega(jh) = (  zomega_T * Wave( ktide(jh) )%nT   &
+            &          + zomega_s * Wave( ktide(jh) )%ns   &
+            &          + zomega_h * Wave( ktide(jh) )%nh   &
+            &          + zomega_p * Wave( ktide(jh) )%np   &
+            &          + zomega_p1* Wave( ktide(jh) )%np1  ) * zscale
+      END DO
+      !
+   END SUBROUTINE tide_pulse
+
+
+   SUBROUTINE tide_vuf( pvt, put, pcor, ktide ,kc )
+      !!----------------------------------------------------------------------
+      !!                     ***  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
+      !!----------------------------------------------------------------------
+      INTEGER                , INTENT(in ) ::   kc               ! Total number of tidal constituents
+      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide            ! Indice of tidal constituents
+      REAL(wp), DIMENSION(kc), INTENT(out) ::   pvt, put, pcor   !
+      !
+      INTEGER ::   jh   ! dummy loop index
+      !!----------------------------------------------------------------------
+      !
+      DO jh = 1, kc
+         !  Phase of the tidal potential relative to the Greenwhich 
+         !  meridian (e.g. the position of the fictuous celestial body). Units are radian:
+         pvt(jh) = sh_T * Wave( ktide(jh) )%nT    &
+            &    + sh_s * Wave( ktide(jh) )%ns    &
+            &    + sh_h * Wave( ktide(jh) )%nh    &
+            &    + sh_p * Wave( ktide(jh) )%np    &
+            &    + sh_p1* Wave( ktide(jh) )%np1   &
+            &    +        Wave( ktide(jh) )%shift * rad
+         !
+         !  Phase correction u due to nodal motion. Units are radian:
+         put(jh) = sh_xi     * Wave( ktide(jh) )%nksi   &
+            &    + sh_nu     * Wave( ktide(jh) )%nnu0   &
+            &    + sh_nuprim * Wave( ktide(jh) )%nnu1   &
+            &    + sh_nusec  * Wave( ktide(jh) )%nnu2   &
+            &    + sh_R      * Wave( ktide(jh) )%R
+
+         !  Nodal correction factor:
+         pcor(jh) = nodal_factort( Wave( ktide(jh) )%nformula )
+      END DO
+      !
+   END SUBROUTINE tide_vuf
+
+
+   RECURSIVE FUNCTION nodal_factort( kformula ) RESULT( zf )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kformula
+      !
+      REAL(wp) :: zf
+      REAL(wp) :: zs, zf1, zf2
+      !!----------------------------------------------------------------------
+      !
+      SELECT CASE( kformula )
+      !
+      CASE( 0 )                  !==  formule 0, solar waves
+         zf = 1.0
+         !
+      CASE( 1 )                  !==  formule 1, compound waves (78 x 78)
+         zf=nodal_factort(78)
+         zf = zf * zf
+         !
+      CASE ( 2 )                 !==  formule 2, compound waves (78 x 0)  ===  (78) 
+       zf1= nodal_factort(78)
+       zf = nodal_factort( 0)
+       zf = zf1 * zf
        !
-       !  Phase correction u due to nodal motion. Units are radian:
-       put(jh) = sh_xi    *Wave(ktide(jh))%nksi  &
-            +sh_nu    *Wave(ktide(jh))%nnu0  &
-            +sh_nuprim*Wave(ktide(jh))%nnu1  &
-            +sh_nusec *Wave(ktide(jh))%nnu2  &
-            +sh_R     *Wave(ktide(jh))%R
-
-       !  Nodal correction factor:
-       pcor(jh) = nodal_factort(Wave(ktide(jh))%nformula)
-    END DO
-
-  END SUBROUTINE tide_vuf
-
-  recursive function nodal_factort(kformula) result (zf)
-    !!----------------------------------------------------------------------
-    INTEGER, INTENT(IN) :: kformula
-    REAL(wp) :: zf
-    REAL(wp) :: zs,zf1,zf2
-
-    SELECT CASE (kformula)
-
-       !!  formule 0, solar waves
-
-    case ( 0 )
-       zf=1.0
-
-       !! formule 1, compound waves (78 x 78)
-
-    case ( 1 )
-       zf=nodal_factort(78)
-       zf=zf*zf
-
-       !! formule 2, compound waves (78 x 0)  ===  (78) 
-
-    case ( 2 )
-       zf1=nodal_factort(78)
-       zf=nodal_factort(0)
-       zf=zf1*zf
-
-       !! formule 4,  compound waves (78 x 235) 
-
-    case ( 4 )
-       zf1=nodal_factort(78)
-       zf=nodal_factort(235)
-       zf=zf1*zf
-
-       !! formule 5,  compound waves (78 *78 x 235)
-
-    case ( 5 )
-       zf1=nodal_factort(78)
-       zf=nodal_factort(235)
-       zf=zf*zf1*zf1
-
-       !! formule 6,  compound waves (78 *78 x 0)
-
-    case ( 6 )
-       zf1=nodal_factort(78)
-       zf=nodal_factort(0)
-       zf=zf*zf1*zf1 
-
-       !! formule 7, compound waves (75 x 75)
-
-    case ( 7 )
-       zf=nodal_factort(75)
-       zf=zf*zf
-
-       !! formule 8,  compound waves (78 x 0 x 235)
-
-    case ( 8 )
-       zf=nodal_factort(78)
-       zf1=nodal_factort(0)
-       zf2=nodal_factort(235)
-       zf=zf*zf1*zf2
-
-       !! formule 9,  compound waves (78 x 0 x 227)
-
-    case ( 9 )
-       zf=nodal_factort(78)
-       zf1=nodal_factort(0)
-       zf2=nodal_factort(227)
-       zf=zf*zf1*zf2
-
-       !! formule 10,  compound waves (78 x 227)
-
-    case ( 10 )
-       zf=nodal_factort(78)
-       zf1=nodal_factort(227)
-       zf=zf*zf1
-
-       !! formule 11,  compound waves (75 x 0)
-
-    case ( 11 )
-       zf=nodal_factort(75)
-       zf=nodal_factort(0)
-       zf=zf*zf1
-
-       !! formule 12,  compound waves (78 x 78 x 78 x 0) 
-
-    case ( 12 )
-       zf1=nodal_factort(78)
-       zf=nodal_factort(0)
-       zf=zf*zf1*zf1*zf1
-
-       !! formule 13, compound waves (78 x 75)
-
-    case ( 13 )
-       zf1=nodal_factort(78)
-       zf=nodal_factort(75)
-       zf=zf*zf1
-
-       !! formule 14, compound waves (235 x 0)  ===  (235)
-
-    case ( 14 )
-       zf=nodal_factort(235)
-       zf1=nodal_factort(0)
-       zf=zf*zf1
-
-       !! formule 15, compound waves (235 x 75) 
-
-    case ( 15 )
-       zf=nodal_factort(235)
-       zf1=nodal_factort(75)
-       zf=zf*zf1
-
-       !! formule 16, compound waves (78 x 0 x 0)  ===  (78)
-
-    case ( 16 )
-       zf=nodal_factort(78)
-       zf1=nodal_factort(0)
-       zf=zf*zf1*zf1
-
-       !! formule 17,  compound waves (227 x 0) 
-
-    case ( 17 )
-       zf1=nodal_factort(227)
-       zf=nodal_factort(0)
-       zf=zf*zf1
-
-       !! formule 18,  compound waves (78 x 78 x 78 )
-
-    case ( 18 ) 
-       zf1=nodal_factort(78)
-       zf=zf1*zf1*zf1
-
-       !! formule 19, compound waves (78 x 0 x 0 x 0)  ===  (78)
-
-    case ( 19 )
-       zf=nodal_factort(78)
-       zf1=nodal_factort(0)
-       zf=zf*zf1*zf1
-
-       !! formule 73
-
-    case ( 73 )
-       zs=sin(sh_I)
-       zf=(2./3.-zs*zs)/0.5021
-
-       !! formule 74
-
-    case ( 74 )
-       zs=sin(sh_I)
-       zf=zs*zs/0.1578
-
-       !! formule 75
-
-    case ( 75 )
-       zs=cos (sh_I/2)
-       zf=sin (sh_I)*zs*zs/0.3800
-
-       !! formule 76
-
-    case ( 76 )
-       zf=sin (2*sh_I)/0.7214
-
-       !! formule 77
-
-    case ( 77 )
-       zs=sin (sh_I/2)
-       zf=sin (sh_I)*zs*zs/0.0164
-
-       !! formule 78
-
-    case ( 78 )
-       zs=cos (sh_I/2)
-       zf=zs*zs*zs*zs/0.9154
-
-       !! formule 79
-
-    case ( 79 )
-       zs=sin(sh_I)
-       zf=zs*zs/0.1565
-
-       !! formule 144
-
-    case ( 144 )
-       zs=sin (sh_I/2)
-       zf=(1-10*zs*zs+15*zs*zs*zs*zs)*cos(sh_I/2)/0.5873
-
-       !! formule 149
-
-    case ( 149 )
-       zs=cos (sh_I/2)
-       zf=zs*zs*zs*zs*zs*zs/0.8758
-
-       !! formule 215
-
-    case ( 215 )
-       zs=cos (sh_I/2)
-       zf=zs*zs*zs*zs/0.9154*sh_x1ra
-
-       !! formule 227 
-
-    case ( 227 )
-       zs=sin (2*sh_I)
-       zf=sqrt (0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006)
-
-       !! formule 235 
-
-    case ( 235 )
-       zs=sin (sh_I)
-       zf=sqrt (19.0444*zs*zs*zs*zs+2.7702*zs*zs*cos (2*sh_nu)+.0981)
-
-    END SELECT
-
-  end function nodal_factort
-
-  function dayjul(kyr,kmonth,kday)
-    !
-    !*** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE)
-    !
-    INTEGER,INTENT(IN) :: kyr,kmonth,kday
-    INTEGER,DIMENSION(12) ::  idayt,idays
-    INTEGER :: inc,ji
-    REAL(wp) :: dayjul,zyq
-
-    DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./
-    idays(1)=0.
-    idays(2)=31.
-    inc=0.
-    zyq=MOD((kyr-1900.),4.)
-    IF(zyq .eq. 0.) inc=1.
-    DO ji=3,12
-       idays(ji)=idayt(ji)+inc
-    END DO
-    dayjul=idays(kmonth)+kday
-
-  END FUNCTION dayjul
-
+      CASE ( 4 )                 !==  formule 4,  compound waves (78 x 235) 
+         zf1 = nodal_factort( 78)
+         zf  = nodal_factort(235)
+         zf  = zf1 * zf
+         !
+      CASE ( 5 )                 !==  formule 5,  compound waves (78 *78 x 235)
+         zf1 = nodal_factort( 78)
+         zf  = nodal_factort(235)
+         zf  = zf * zf1 * zf1
+         !
+      CASE ( 6 )                 !==  formule 6,  compound waves (78 *78 x 0)
+         zf1 = nodal_factort(78)
+         zf  = nodal_factort( 0)
+         zf  = zf * zf1 * zf1 
+         !
+      CASE( 7 )                  !==  formule 7, compound waves (75 x 75)
+         zf = nodal_factort(75)
+         zf = zf * zf
+         !
+      CASE( 8 )                  !==  formule 8,  compound waves (78 x 0 x 235)
+         zf  = nodal_factort( 78)
+         zf1 = nodal_factort(  0)
+         zf2 = nodal_factort(235)
+         zf  = zf * zf1 * zf2
+         !
+      CASE( 9 )                  !==  formule 9,  compound waves (78 x 0 x 227)
+         zf  = nodal_factort( 78)
+         zf1 = nodal_factort(  0)
+         zf2 = nodal_factort(227)
+         zf  = zf * zf1 * zf2
+         !
+      CASE( 10 )                 !==  formule 10,  compound waves (78 x 227)
+         zf  = nodal_factort( 78)
+         zf1 = nodal_factort(227)
+         zf  = zf * zf1
+         !
+      CASE( 11 )                 !==  formule 11,  compound waves (75 x 0)
+!!gm bug???? zf 2 fois !
+         zf = nodal_factort(75)
+         zf = nodal_factort( 0)
+         zf = zf * zf1
+         !
+      CASE( 12 )                 !==  formule 12,  compound waves (78 x 78 x 78 x 0) 
+         zf1 = nodal_factort(78)
+         zf  = nodal_factort( 0)
+         zf  = zf * zf1 * zf1 * zf1
+         !
+      CASE( 13 )                 !==  formule 13, compound waves (78 x 75)
+         zf1 = nodal_factort(78)
+         zf  = nodal_factort(75)
+         zf  = zf * zf1
+         !
+      CASE( 14 )                 !==  formule 14, compound waves (235 x 0)  ===  (235)
+         zf  = nodal_factort(235)
+         zf1 = nodal_factort(  0)
+         zf  = zf * zf1
+         !
+      CASE( 15 )                 !==  formule 15, compound waves (235 x 75) 
+         zf  = nodal_factort(235)
+         zf1 = nodal_factort( 75)
+         zf  = zf * zf1
+         !
+      CASE( 16 )                 !==  formule 16, compound waves (78 x 0 x 0)  ===  (78)
+         zf  = nodal_factort(78)
+         zf1 = nodal_factort( 0)
+         zf  = zf * zf1 * zf1
+         !
+      CASE( 17 )                 !==  formule 17,  compound waves (227 x 0) 
+         zf1 = nodal_factort(227)
+         zf  = nodal_factort(  0)
+         zf  = zf * zf1
+         !
+      CASE( 18 )                 !==  formule 18,  compound waves (78 x 78 x 78 )
+         zf1 = nodal_factort(78)
+         zf  = zf1 * zf1 * zf1
+         !
+      CASE( 19 )                 !==  formule 19, compound waves (78 x 0 x 0 x 0)  ===  (78)
+!!gm bug2 ==>>>   here identical to formule 16,  a third multiplication by zf1 is missing
+         zf  = nodal_factort(78)
+         zf1 = nodal_factort( 0)
+         zf = zf * zf1 * zf1
+         !
+      CASE( 73 )                 !==  formule 73
+         zs = sin(sh_I)
+         zf = (2./3.-zs*zs)/0.5021
+         !
+      CASE( 74 )                 !==  formule 74
+         zs = sin(sh_I)
+         zf = zs * zs / 0.1578
+         !
+      CASE( 75 )                 !==  formule 75
+         zs = cos(sh_I/2)
+         zf = sin(sh_I) * zs * zs / 0.3800
+         !
+      CASE( 76 )                 !==  formule 76
+         zf = sin(2*sh_I) / 0.7214
+         !
+      CASE( 77 )                 !==  formule 77
+         zs = sin(sh_I/2)
+         zf = sin(sh_I) * zs * zs / 0.0164
+         !
+      CASE( 78 )                 !==  formule 78
+         zs = cos(sh_I/2)
+         zf = zs * zs * zs * zs / 0.9154
+         !
+      CASE( 79 )                 !==  formule 79
+         zs = sin(sh_I)
+         zf = zs * zs / 0.1565
+         !
+      CASE( 144 )                !==  formule 144
+         zs = sin(sh_I/2)
+         zf = ( 1-10*zs*zs+15*zs*zs*zs*zs ) * cos(sh_I/2) / 0.5873
+         !
+      CASE( 149 )                !==  formule 149
+         zs = cos(sh_I/2)
+         zf = zs*zs*zs*zs*zs*zs / 0.8758
+         !
+      CASE( 215 )                !==  formule 215
+         zs = cos(sh_I/2)
+         zf = zs*zs*zs*zs / 0.9154 * sh_x1ra
+         !
+      CASE( 227 )                !==  formule 227 
+         zs = sin(2*sh_I)
+         zf = sqrt( 0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006 )
+         !
+      CASE ( 235 )               !==  formule 235 
+         zs = sin(sh_I)
+         zf = sqrt( 19.0444*zs*zs*zs*zs + 2.7702*zs*zs*cos(2*sh_nu) + .0981 )
+         !
+      END SELECT
+      !
+   END FUNCTION nodal_factort
+
+
+   FUNCTION dayjul( kyr, kmonth, kday )
+      !!----------------------------------------------------------------------
+      !!  *** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE)
+      !!----------------------------------------------------------------------
+      INTEGER,INTENT(in) ::   kyr, kmonth, kday
+      !
+      INTEGER,DIMENSION(12) ::  idayt, idays
+      INTEGER  ::   inc, ji
+      REAL(wp) ::   dayjul, zyq
+      !
+      DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./
+      !!----------------------------------------------------------------------
+      !
+      idays(1) = 0.
+      idays(2) = 31.
+      inc = 0.
+      zyq = MOD( kyr-1900. , 4. )
+      IF( zyq == 0.)   inc = 1.
+      DO ji = 3, 12
+         idays(ji)=idayt(ji)+inc
+      END DO
+      dayjul = idays(kmonth) + kday
+      !
+   END FUNCTION dayjul
+
+   !!======================================================================
 END MODULE tide_mod
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90	(revision 3951)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90	(revision 3953)
@@ -1,126 +1,126 @@
 MODULE tideini
-  !!=================================================================================
-  !!                       ***  MODULE  tideini  ***
-  !! Initialization of tidal forcing
-  !! History :  9.0  !  07  (O. Le Galloudec)  Original code
-  !!=================================================================================
-  !! * Modules used
-  USE oce             ! ocean dynamics and tracers variables
-  USE dom_oce         ! ocean space and time domain
-  USE in_out_manager  ! I/O units
-  USE ioipsl          ! NetCDF IPSL library
-  USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-  USE phycst
-  USE daymod
-  USE dynspg_oce
-  USE tide_mod
-  USE iom
+   !!======================================================================
+   !!                       ***  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
+   USE daymod
+   USE dynspg_oce
+   USE tide_mod
+   !
+   USE iom
+   USE in_out_manager  ! I/O units
+   USE ioipsl          ! NetCDF IPSL library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
 
-  IMPLICIT NONE
-  PUBLIC
+   IMPLICIT NONE
+   PUBLIC
 
-  REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::  &
-       omega_tide,  &
-       v0tide,      &
-       utide,       &
-       ftide
+   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_pot = .false., ln_tide_ramp = .false.
-  REAL(wp), PUBLIC :: rdttideramp 
-  INTEGER, PUBLIC :: nb_harmo
-  INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide
-  INTEGER, PUBLIC :: kt_tide
+   LOGICAL , PUBLIC ::   ln_tide_pot  = .FALSE.   !:
+   LOGICAL , PUBLIC ::   ln_tide_ramp = .FALSE.   !:
+   INTEGER , PUBLIC ::   nb_harmo                 !:
+   INTEGER , PUBLIC ::   kt_tide                  !:
+   REAL(wp), PUBLIC ::   rdttideramp              !:
+   
+   INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:) ::   ntide   !:
 
-  !!---------------------------------------------------------------------------------
-  !!   OPA 9.0 , LODYC-IPSL  (2003)
-  !!---------------------------------------------------------------------------------
-
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.5 , NEMO Consortium (2013)
+   !! $Id: $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
    
-  SUBROUTINE tide_init ( kt )
-    !!----------------------------------------------------------------------
-    !!                 ***  ROUTINE tide_init  ***
-    !!----------------------------------------------------------------------      
-    !! * Local declarations
-    INTEGER  :: ji, jk
-    INTEGER, INTENT( in ) ::   kt     ! ocean time-step
-    CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname
-    !
-    NAMELIST/nam_tide/ln_tide_pot, ln_tide_ramp, rdttideramp, clname
-    !!----------------------------------------------------------------------
+   SUBROUTINE tide_init ( kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE tide_init  ***
+      !!----------------------------------------------------------------------      
+      INTEGER  :: ji, jk
+      INTEGER, INTENT( in ) ::   kt     ! ocean time-step
+      CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname
+      !
+      NAMELIST/nam_tide/ln_tide_pot, ln_tide_ramp, rdttideramp, clname
+      !!----------------------------------------------------------------------
 
-    IF ( kt == nit000 ) THEN
-       !
-       IF(lwp) THEN
-          WRITE(numout,*)
-          WRITE(numout,*) 'tide_init : Initialization of the tidal components'
-          WRITE(numout,*) '~~~~~~~~~ '
-       ENDIF
-       !
-       CALL tide_init_Wave
-       !
-       clname(:)=''
-       !
-       ! Read Namelist nam_tide
-       REWIND ( numnam )
-       READ   ( numnam, nam_tide )
-       !
-       nb_harmo=0
-       DO jk=1,jpmax_harmo
-          DO ji=1,jpmax_harmo
-             IF(TRIM(clname(jk)) .eq. Wave(ji)%cname_tide) THEN
-                nb_harmo=nb_harmo+1
-             ENDIF
-          END DO
-       ENDDO
-       !
-       IF(lwp) THEN
-          WRITE(numout,*) '        Namelist nam_tide'
-          WRITE(numout,*) '        nb_harmo    = ', nb_harmo
-          WRITE(numout,*) '        ln_tide_ramp = ', ln_tide_ramp 
-          WRITE(numout,*) '        rdttideramp = ', rdttideramp
-          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')
-          CALL flush(numout)
-       ENDIF
-       !
-       ALLOCATE(ntide(nb_harmo))
-       DO jk=1,nb_harmo
-          DO ji=1,jpmax_harmo
-             IF (TRIM(clname(jk)) .eq. Wave(ji)%cname_tide) THEN
-                ntide(jk) = ji
-                EXIT
-             END IF
-          END DO
-       END DO
-       !
-       ALLOCATE(omega_tide(nb_harmo))
-       ALLOCATE(v0tide    (nb_harmo))
-       ALLOCATE(utide     (nb_harmo))
-       ALLOCATE(ftide     (nb_harmo))
-       kt_tide = kt
-       !
-    ENDIF
-
-    IF ( nsec_day == NINT(0.5 * rdttra(1)) ) THEN
-       !
-       IF(lwp) THEN
-          WRITE(numout,*)
-          WRITE(numout,*) 'tide_ini : Update of the tidal components at kt=',kt
-          WRITE(numout,*) '~~~~~~~~ '
-       ENDIF
-       CALL tide_harmo(omega_tide, v0tide, utide, ftide, ntide, nb_harmo)
-       DO jk =1,nb_harmo
-         IF(lwp) WRITE(numout,*) Wave(ntide(jk))%cname_tide,utide(jk),ftide(jk),v0tide(jk),omega_tide(jk)
-         call flush(numout)
-       END DO
-       !
-       kt_tide = kt
-       !
-    ENDIF
-
-  END SUBROUTINE tide_init
-   
+      IF( kt == nit000 ) THEN
+         !
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'tide_init : Initialization of the tidal components'
+            WRITE(numout,*) '~~~~~~~~~ '
+         ENDIF
+         !
+         CALL tide_init_Wave
+         !
+         clname(:)=''
+         !
+         REWIND( numnam )            ! Read Namelist nam_tide
+         READ  ( numnam, nam_tide )
+         !
+         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
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) '   Namelist nam_tide'
+            WRITE(numout,*) '      Apply astronomical potential : ln_tide_pot  =', ln_tide_pot
+            WRITE(numout,*) '                                     nb_harmo     = ', nb_harmo
+            WRITE(numout,*) '                                     ln_tide_ramp = ', ln_tide_ramp 
+            WRITE(numout,*) '                                     rdttideramp  = ', rdttideramp
+         ENDIF
+         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')
+         !
+         IF( .NOT. lk_dynspg_ts )   CALL ctl_warn( 'sbc_tide : use of time splitting is recommended' )
+         !
+         ALLOCATE( ntide(nb_harmo) )
+         DO jk = 1, nb_harmo
+            DO ji = 1, jpmax_harmo
+               IF( TRIM(clname(jk)) .eq. Wave(ji)%cname_tide ) THEN
+                  ntide(jk) = ji
+                  EXIT
+               END IF
+            END DO
+         END DO
+         !
+         ALLOCATE( omega_tide(nb_harmo), v0tide    (nb_harmo),   &
+            &      utide     (nb_harmo), ftide     (nb_harmo)  )
+         kt_tide = kt
+         !
+      ENDIF
+      !
+      IF( nsec_day == NINT(0.5 * rdttra(1)) ) THEN
+         !
+         CALL tide_harmo( omega_tide, v0tide, utide, ftide, ntide, nb_harmo )
+         !
+         kt_tide = kt
+         !
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'tide_ini : Update of the tidal components 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
+         !
+      ENDIF
+      !
+   END SUBROUTINE tide_init
+     
+   !!======================================================================
 END MODULE tideini
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90	(revision 3951)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90	(revision 3953)
@@ -1,65 +1,71 @@
 MODULE updtide
-  !!=================================================================================
-  !!                       ***  MODULE  updtide  ***
-  !! Initialization of tidal forcing
-  !! History :  9.0  !  07  (O. Le Galloudec)  Original code
-  !!=================================================================================
+   !!======================================================================
+   !!                       ***  MODULE  updtide  ***
+   !! Initialization of tidal forcing
+   !!======================================================================
+   !! History :  9.0  !  07  (O. Le Galloudec)  Original code
+   !!----------------------------------------------------------------------
 #if defined key_tide
-  !! * Modules used
-  USE oce             ! ocean dynamics and tracers variables
-  USE dom_oce         ! ocean space and time domain
-  USE in_out_manager  ! I/O units
-  USE phycst
-  USE sbctide
-  USE dynspg_oce
-  USE tideini, ONLY: ln_tide_ramp, rdttideramp
+   !!----------------------------------------------------------------------
+   !!   'key_tide' :                                        tidal potential
+   !!----------------------------------------------------------------------
+   !!   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
+   IMPLICIT NONE
+   PUBLIC
 
-  !! * Routine accessibility
-  PUBLIC upd_tide
-  !!---------------------------------------------------------------------------------
-  !!   OPA 9.0 , LODYC-IPSL  (2003)
-  !!---------------------------------------------------------------------------------
-
+   PUBLIC   upd_tide   ! called in dynspg_... modules
+  
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: sbcfwb.F90 3625 2012-11-21 13:19:18Z acc $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
-  SUBROUTINE upd_tide (kt,kit)
-    !!----------------------------------------------------------------------
-    !!                 ***  ROUTINE upd_tide  ***
-    !!----------------------------------------------------------------------      
-    !! * Local declarations
+   SUBROUTINE upd_tide( kt, kit, kbaro )
+      !!----------------------------------------------------------------------
+      !!                 ***  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 only)
+      INTEGER, INTENT(in), OPTIONAL ::   kbaro   ! number of sub-time-step           (lk_dynspg_ts=T only)
+      !
+      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
+      IF( PRESENT( kit ) .AND. PRESENT( kbaro ) )   zt = zt + kit * rdt / REAL( kbaro, wp )
+      zwt(:) = omega_tide(:) * zt
 
-    INTEGER, INTENT( in ) ::   kt,kit      ! ocean time-step index
-    INTEGER  :: ji,jj,jk
-    REAL (wp) :: zramp
-    REAL (wp), DIMENSION(nb_harmo) :: zwt 
-    !...............................................................................
-
-    pot_astro(:,:)=0.e0
-    zramp = 1.e0
-
-    IF (lk_dynspg_ts) THEN
-       zwt(:) = omega_tide(:)* ((kt-kt_tide)*rdt + kit*(rdt/REAL(nn_baro,wp)))
-       IF (ln_tide_ramp) THEN
-          zramp = MIN(MAX( ((kt-nit000)*rdt + kit*(rdt/REAL(nn_baro,wp)))/(rdttideramp*rday),0.),1.)
-       ENDIF
-    ELSE
-       zwt(:) = omega_tide(:)*(kt-kt_tide)*rdt
-       IF (ln_tide_ramp) THEN
-          zramp = MIN(MAX( ((kt-nit000)*rdt)/(rdttideramp*rday),0.),1.) 
-       ENDIF  
-    ENDIF
-
-    do jk=1,nb_harmo
-       do ji=1,jpi
-          do jj=1,jpj
-             pot_astro(ji,jj)=pot_astro(ji,jj) + zramp*(amp_pot(ji,jj,jk)*COS(zwt(jk)+phi_pot(ji,jj,jk)))      
-          enddo
-       enddo
-    enddo
-
-  END SUBROUTINE upd_tide
+      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 ) .AND. PRESENT( kbaro ) )   zt = zt + kit * rdt / REAL( kbaro, wp )
+         zramp = MIN(  MAX( zt / (rdttideramp*rday) , 0._wp ) , 1._wp  )
+         pot_astro(:,:) = zramp * pot_astro(:,:)
+      ENDIF
+      !
+   END SUBROUTINE upd_tide
 
 #else
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 3951)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 3953)
@@ -89,10 +89,9 @@
       ! Update data, open boundaries, surface boundary condition (including sea-ice)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
-      IF( lk_tide.AND.(kstp /= nit000 ))   CALL tide_init ( kstp )
+                         CALL sbc     ( kstp )        ! Sea Boundary Condition (including sea-ice)
       IF( lk_tide    )   CALL sbc_tide( kstp )
-      IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries
-      IF( lk_obc     )   CALL obc_rad( kstp )         ! compute phase velocities at open boundaries
-      IF( lk_bdy     )   CALL bdy_dta( kstp, time_offset=+1 ) ! update dynamic and tracer data at open boundaries
+      IF( lk_obc     )   CALL obc_dta ( kstp )        ! update dynamic and tracer data at open boundaries
+      IF( lk_obc     )   CALL obc_rad ( kstp )        ! compute phase velocities at open boundaries
+      IF( lk_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -117,5 +116,5 @@
       IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz
       IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz
-      IF( lk_zdfcst  )   THEN                            ! Constant Kz (reset avt, avm[uv] to the background value)
+      IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value)
          avt (:,:,:) = rn_avt0 * tmask(:,:,:)
          avmu(:,:,:) = rn_avm0 * umask(:,:,:)
