Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 2813)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 2814)
@@ -26,4 +26,5 @@
    USE zdfbfr          ! bottom friction
    USE dynvor          ! vorticity term
+   USE obc_par         ! for lk_obc
    USE obc_oce         ! Lateral open boundary condition
    USE obcdta          ! open boundary condition data     
@@ -528,11 +529,4 @@
       !                                                    ! ==================== !
 
-#if defined key_obc
-      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ?????
-      IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
-      IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
-      IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
-#endif
-
       ! -----------------------------------------------------------------------------
       ! Phase 3. update the general trend with the barotropic trend
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/FLO/floblk.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/FLO/floblk.F90	(revision 2813)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/FLO/floblk.F90	(revision 2814)
@@ -345,34 +345,35 @@
       ! more time.      
 # if defined key_obc
-      DO jfl = 1, jpnfl
-         IF( lp_obc_east ) THEN
-            IF( jped <=  zgjfl(jfl) .AND. zgjfl(jfl) <= jpef .AND. nieob-1 <=  zgifl(jfl) ) THEN
-               zgifl (jfl) = INT(zgifl(jfl)) + 0.5
-               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5
-               zagefl(jfl) = rdt
-            END IF
-         END IF
-         IF( lp_obc_west ) THEN
-            IF( jpwd <= zgjfl(jfl) .AND. zgjfl(jfl) <= jpwf .AND. niwob >=  zgifl(jfl) ) THEN
-               zgifl (jfl) = INT(zgifl(jfl)) + 0.5
-               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5
-               zagefl(jfl) = rdt
-            END IF
-         END IF
-         IF( lp_obc_north ) THEN
-            IF( jpnd <=  zgifl(jfl) .AND. zgifl(jfl) <= jpnf .AND. njnob-1 >=  zgjfl(jfl) ) THEN
-               zgifl (jfl) = INT(zgifl(jfl)) + 0.5
-               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5
-               zagefl(jfl) = rdt
-            END IF
-         END IF
-         IF( lp_obc_south ) THEN
-            IF( jpsd <=  zgifl(jfl) .AND. zgifl(jfl) <= jpsf .AND.  njsob >= zgjfl(jfl) ) THEN
-               zgifl (jfl) = INT(zgifl(jfl)) + 0.5
-               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5
-               zagefl(jfl) = rdt
-            END IF
-         END IF
-      END DO
+!!!!!!!! NEED TO SORT THIS OUT !!!!!!!!
+!!$      DO jfl = 1, jpnfl
+!!$         IF( lp_obc_east ) THEN
+!!$            IF( jped <=  zgjfl(jfl) .AND. zgjfl(jfl) <= jpef .AND. nieob-1 <=  zgifl(jfl) ) THEN
+!!$               zgifl (jfl) = INT(zgifl(jfl)) + 0.5
+!!$               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5
+!!$               zagefl(jfl) = rdt
+!!$            END IF
+!!$         END IF
+!!$         IF( lp_obc_west ) THEN
+!!$            IF( jpwd <= zgjfl(jfl) .AND. zgjfl(jfl) <= jpwf .AND. niwob >=  zgifl(jfl) ) THEN
+!!$               zgifl (jfl) = INT(zgifl(jfl)) + 0.5
+!!$               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5
+!!$               zagefl(jfl) = rdt
+!!$            END IF
+!!$         END IF
+!!$         IF( lp_obc_north ) THEN
+!!$            IF( jpnd <=  zgifl(jfl) .AND. zgifl(jfl) <= jpnf .AND. njnob-1 >=  zgjfl(jfl) ) THEN
+!!$               zgifl (jfl) = INT(zgifl(jfl)) + 0.5
+!!$               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5
+!!$               zagefl(jfl) = rdt
+!!$            END IF
+!!$         END IF
+!!$         IF( lp_obc_south ) THEN
+!!$            IF( jpsd <=  zgifl(jfl) .AND. zgifl(jfl) <= jpsf .AND.  njsob >= zgjfl(jfl) ) THEN
+!!$               zgifl (jfl) = INT(zgifl(jfl)) + 0.5
+!!$               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5
+!!$               zagefl(jfl) = rdt
+!!$            END IF
+!!$         END IF
+!!$      END DO
 #endif
 
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obc_oce.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obc_oce.F90	(revision 2813)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obc_oce.F90	(revision 2814)
@@ -56,5 +56,4 @@
    !                                                        !: =F read obc coordinates from namelist
    LOGICAL                    ::   ln_mask_file             !: =T read obcmask from file
-   LOGICAL, DIMENSION(jp_obc) ::   ln_tides                 !: =T apply tidal harmonic forcing along open boundaries
    LOGICAL                    ::   ln_vol                   !: =T volume correction             
    LOGICAL, DIMENSION(jp_obc) ::   ln_clim                  !: =T obc data files contain climatological data (time-cyclic)
@@ -64,4 +63,7 @@
    INTEGER, DIMENSION(jp_obc) ::   nn_dtactl      	    !: = 0 use the initial state as obc dta ; 
                                                             !: = 1 read it in a NetCDF file
+   INTEGER, DIMENSION(jp_obc) ::   nn_tides                 !: = 0 no tidal harmonic forcing
+                                                            !: = 1 apply ONLY tidal harmonic forcing for barotropic solution
+                                                            !: = 2 ADD tidal harmonic forcing to other barotropic boundary data
    INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P 
    !                                                        !  = 1 the volume will be constant during all the integration.
@@ -98,7 +100,7 @@
    !!----------------------------------------------------------------------
 
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)       ::   dta_global        !: workspace for reading in global data arrays
-   TYPE(OBC_INDEX), DIMENSION(jp_obc), TARGET    ::   idx_obc           !: obc indices (local process)
-   TYPE(OBC_DATA) , DIMENSION(jp_obc)            ::   dta_obc           !: obc external data (local process)
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays
+   TYPE(OBC_INDEX), DIMENSION(jp_obc), TARGET      ::   idx_obc           !: obc indices (local process)
+   TYPE(OBC_DATA) , DIMENSION(jp_obc)              ::   dta_obc           !: obc external data (local process)
 
    !!----------------------------------------------------------------------
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90	(revision 2813)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90	(revision 2814)
@@ -65,5 +65,5 @@
       !!----------------------------------------------------------------------
       USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY: wrk_2d_9, wrk_2d_10      ! 2D workspace
+      USE wrk_nemo, ONLY: wrk_2d_22, wrk_2d_23   ! 2D workspace
       !!
       INTEGER, INTENT( in )           ::   kt    ! ocean time-step index 
@@ -76,5 +76,5 @@
       !!---------------------------------------------------------------------------
 
-      IF(wrk_in_use(2, 9,10) ) THEN
+      IF(wrk_in_use(2, 22,23) ) THEN
          CALL ctl_stop('obc_dta: ERROR: requested workspace arrays are unavailable.')   ;   RETURN
       END IF
@@ -87,6 +87,6 @@
          ! Calculate depth-mean currents
          !-----------------------------
-         pu2d => wrk_2d_9
-         pu2d => wrk_2d_10
+         pu2d => wrk_2d_22
+         pu2d => wrk_2d_23
 
          pu2d(:,:) = 0.e0
@@ -206,11 +206,20 @@
                ! Update barotropic boundary conditions only
                ! jit is optional argument for fld_read
-               IF( nn_dyn2d(ib_obc) .gt. 0 ) THEN
+               IF( nn_dyn2d(ib_obc) .gt. 0 .and. nn_tides(ib_obc) .ne. 1 ) THEN
                   jend = jstart + 2
                   CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit )
                ENDIF
+               IF( nn_tides(ib_obc) .gt. 0 ) THEN
+                  CALL tide_update( kt=kt, jit=jit, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc) )
+               ENDIF
             ELSE
-               jend = jstart + nb_obc_fld(ib_obc) - 1
-               CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend ), map=nbmap_ptr(jstart:jend), timeshift=1 )
+               IF( nb_obc_fld(ib_obc) .gt. 0 ) THEN 
+                  jend = jstart + nb_obc_fld(ib_obc) - 1
+                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), timeshift=1 )
+               ENDIF
+               IF( nn_tides(ib_obc) .gt. 0 ) THEN
+                  !!! THINK ABOUT kt, jit VALUES !!!
+                  CALL tide_update( kt=kt, jit=0, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc) )
+               ENDIF
             ENDIF
             jstart = jend+1
@@ -257,5 +266,5 @@
       END DO  ! ib_obc
 
-      IF(wrk_not_released(2, 9,10) )    CALL ctl_stop('obc_dta: ERROR: failed to release workspace arrays.')
+      IF(wrk_not_released(2, 22,23) )    CALL ctl_stop('obc_dta: ERROR: failed to release workspace arrays.')
 
       END SUBROUTINE obc_dta
@@ -383,5 +392,5 @@
             IF( nn_dyn2d(ib_obc) .gt. 0 ) THEN 
 
-               IF( nn_dyn2d(ib_obc) .ne. jp_frs ) THEN
+               IF( nn_dyn2d(ib_obc) .ne. jp_frs .and. nn_tides(ib_obc) .ne. 1) THEN
                   jfld = jfld + 1
                   blf_i(jfld) = bn_ssh
@@ -392,5 +401,5 @@
                ENDIF
 
-               IF( .not. ln_full_vel_array(ib_obc) ) THEN
+               IF( .not. ln_full_vel_array(ib_obc) .and. nn_tides(ib_obc) .ne. 1 ) THEN
 
                   jfld = jfld + 1
@@ -604,21 +613,20 @@
 
             ! nn_dtactl = 1
-            ! Set boundary data arrays to point to relevant "fnow" arrays
-            !-----------------------------------------------------------
+            ! Set boundary data arrays to point to "fnow" arrays
+            !---------------------------------------------------
             IF (nn_dyn2d(ib_obc) .gt. 0) THEN
-               IF( nn_dyn2d(ib_obc) .ne. jp_frs ) THEN
+               IF( nn_dyn2d(ib_obc) .ne. jp_frs .and. nn_tides(ib_obc) .ne. 1 ) THEN
                   jfld = jfld + 1
                   dta_obc(ib_obc)%ssh => bf(jfld)%fnow(:,1,1)
                ENDIF
-               IF( ln_full_vel_array(ib_obc) ) THEN
+               IF( ln_full_vel_array(ib_obc) .or. nn_tides(ib_obc) .eq. 1 ) THEN
                   ! In this case we need space but we aren't reading it 
                   ! directly from the external file. 
                   IF( nn_dyn2d(ib_obc) .eq. jp_frs ) THEN
-                     ilen0(2) = nblen(2)
-                     ilen0(3) = nblen(3)
+                     ilen0(:) = nblen(:)
                   ELSE
-                     ilen0(2) = nblenrim(2)
-                     ilen0(3) = nblenrim(3)
+                     ilen0(:) = nblenrim(:)
                   ENDIF
+                  IF( nn_tides(ib_obc) .eq. 1 ) ALLOCATE( dta_obc(ib_obc)%ssh(ilen0(1)) )
                   ALLOCATE( dta_obc(ib_obc)%u2d(ilen0(2)) )
                   ALLOCATE( dta_obc(ib_obc)%v2d(ilen0(3)) )
@@ -631,5 +639,5 @@
             ENDIF
             IF (nn_dyn3d(ib_obc) .gt. 0 .or. &
-              &   ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 ) THEN
+              &  ( ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 ) ) THEN
                jfld = jfld + 1
                dta_obc(ib_obc)%u3d => bf(jfld)%fnow(:,1,:)
@@ -666,4 +674,6 @@
 CONTAINS
    SUBROUTINE obc_dta( kt, jit )              ! Empty routine
+      INTEGER, INTENT( in )           ::   kt    
+      INTEGER, INTENT( in ), OPTIONAL ::   jit   
       WRITE(*,*) 'obc_dta: You should not have seen this print! error?', kt
    END SUBROUTINE obc_dta
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn2d.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn2d.F90	(revision 2813)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn2d.F90	(revision 2814)
@@ -24,5 +24,4 @@
    USE phycst          ! physical constants
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE obctides        ! for tidal harmonic forcing at boundary
    USE in_out_manager  !
 
@@ -30,5 +29,5 @@
    PRIVATE
 
-   PUBLIC   obc_dyn2d     ! routine called in dynspg_ts (time splitting case ONLY)
+   PUBLIC   obc_dyn2d     ! routine called in dynspg_ts and dyn_nxt and dynspg_flt
 
    !!----------------------------------------------------------------------
@@ -77,5 +76,5 @@
       !!               topography. Tellus, 365-382.
       !!----------------------------------------------------------------------
-      INTEGER,         INTENT(in) :: kt
+      INTEGER,         INTENT(in) ::   kt
       TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
       TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
@@ -114,5 +113,4 @@
       !!             
       !!              - Apply Flather boundary conditions on normal barotropic velocities 
-      !!                (ln_dyn_fla=.true. or ln_tides=.true.)
       !!
       !! ** WARNINGS about FLATHER implementation:
@@ -143,6 +141,4 @@
       ! ---------------------------------! 
      
-!!!!!!!!!!!! SOME WORK TO DO ON THE TIDES HERE !!!!!!!!!!!!!
-
 !!! REPLACE spgu with nemo_wrk work space
 
@@ -154,5 +150,4 @@
          ij = idx%nbj(jb,igrd)
          spgu(ii, ij) = dta%ssh(jb)
-!!$         IF( ln_tides )   spgu(ii, ij) = spgu(ii, ij) + sshtide(jb)
       END DO
       !
@@ -167,5 +162,4 @@
          !
          zcorr = - idx%flagu(jb) * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
-!!$         zforc = dta%u2d(jb) + utide(jb)
          zforc = dta%u2d(jb)
          pu2d(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
@@ -181,5 +175,4 @@
          !
          zcorr = - idx%flagv(jb) * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
-!!$         zforc = dta%v2d(jb) + vtide(jb)
          zforc = dta%v2d(jb)
          pv2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1)
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90	(revision 2813)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90	(revision 2814)
@@ -22,5 +22,4 @@
    USE dom_oce         ! ocean space and time domain
    USE obc_oce         ! unstructured open boundary conditions
-   USE obctides        ! tides at open boundaries initialization (tide_init routine)
    USE in_out_manager  ! I/O units
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -73,5 +72,5 @@
          &             nn_ice_lim2,                                        &
 #endif
-         &             ln_tides, ln_vol, ln_clim, nn_dtactl, nn_volctl,    &
+         &             nn_tides, ln_vol, ln_clim, nn_dtactl, nn_volctl,    &
          &             nn_rimwidth, nn_dmp2d_in, nn_dmp2d_out,             &
          &             nn_dmp3d_in, nn_dmp3d_out
@@ -105,8 +104,8 @@
       nn_ice_lim2(:)    = 0
 #endif
-      ln_tides(:)       = .false.
       ln_vol            = .false.
       ln_clim(:)        = .false.
       nn_dtactl(:)      = -1  ! uninitialised flag
+      nn_tides(:)       =  0  ! default to no tidal forcing
       nn_volctl         = -1  ! uninitialised flag
       nn_rimwidth(:)    = -1  ! uninitialised flag
@@ -182,12 +181,17 @@
         IF(lwp) WRITE(numout,*)
 
-        IF( ln_tides(ib_obc) ) THEN
-          IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing at unstructured open boundaries'
+        SELECT CASE( nn_tides(ib_obc) ) 
+        CASE(0)
+          IF(lwp) WRITE(numout,*) 'No tidal harmonic forcing'
           IF(lwp) WRITE(numout,*)
-        ENDIF
-
-!!$        ! Presumably will need to read in a separate namelist for each boundary that includes tides???
-!!$        IF( ln_tides )   CALL tide_init( ib_obc )      ! Read tides namelist 
-
+        CASE(1)
+          IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ONLY for barotropic solution'
+          IF(lwp) WRITE(numout,*)
+        CASE(2)
+          IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ADDED to other barotropic boundary conditions'
+          IF(lwp) WRITE(numout,*)
+        CASE DEFAULT
+          CALL ctl_stop( 'obc_ini: ERROR: incorrect value for nn_tides ' )
+        END SELECT
 
       ENDDO
@@ -345,6 +349,6 @@
          ENDDO 
 
-         ! Compute rim weights
-         ! -------------------
+         ! Compute rim weights for FRS scheme
+         ! ----------------------------------
          DO igrd = 1, jpbgrd
             DO ib = 1, idx_obc(ib_obc)%nblen(igrd)
@@ -537,8 +541,4 @@
          IF( lk_mpp )   CALL mpp_sum( obcsurftot )      ! sum over the global domain
       END IF   
-
-      ! Read in tidal constituents and adjust for model start time
-      ! ----------------------------------------------------------
-!!$      IF( ln_tides )   CALL tide_data
       !
       ! Tidy up
@@ -550,5 +550,5 @@
 #else
    !!---------------------------------------------------------------------------------
-   !!   Dummy module                                   NO unstructured open boundaries
+   !!   Dummy module                                   NO open boundaries
    !!---------------------------------------------------------------------------------
 CONTAINS
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obctides.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obctides.F90	(revision 2813)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obctides.F90	(revision 2814)
@@ -10,709 +10,690 @@
    !!----------------------------------------------------------------------
 #if defined key_obc
-!!$   !!----------------------------------------------------------------------
-!!$   !!   'key_obc'     Unstructured Open Boundary Condition
-!!$   !!----------------------------------------------------------------------
-!!$   !!   PUBLIC
-!!$   !!      tide_init     : read of namelist
-!!$   !!      tide_data     : read in and initialisation of tidal constituents at boundary
-!!$   !!      tide_update   : calculation of tidal forcing at each timestep
-!!$   !!   PRIVATE
-!!$   !!      uvset         :\
-!!$   !!      vday          : | Routines to correct tidal harmonics forcing for
-!!$   !!      shpen         : | start time of integration
-!!$   !!      ufset         : |
-!!$   !!      vset          :/
-!!$   !!----------------------------------------------------------------------
-!!$   USE oce             ! ocean dynamics and tracers 
-!!$   USE dom_oce         ! ocean space and time domain
-!!$   USE iom
-!!$   USE in_out_manager  ! I/O units
-!!$   USE phycst          ! physical constants
-!!$   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-!!$   USE obc_par         ! Unstructured boundary parameters
-!!$   USE obc_oce         ! ocean open boundary conditions
-!!$   USE daymod          ! calendar
-!!$
-!!$   IMPLICIT NONE
-!!$   PRIVATE
-!!$
-!!$   PUBLIC   tide_init     ! routine called in obcini
-!!$   PUBLIC   tide_data     ! routine called in obcini
-!!$   PUBLIC   tide_update   ! routine called in obcdyn
-!!$
-!!$   LOGICAL, PUBLIC            ::   ln_tide_date          !: =T correct tide phases and amplitude for model start date
-!!$   INTEGER, PUBLIC, PARAMETER ::   jptides_max = 15      !: Max number of tidal contituents
-!!$   INTEGER, PUBLIC            ::   ntide                 !: Actual number of tidal constituents
-!!$
-!!$   CHARACTER(len=80), PUBLIC                         ::   filtide    !: Filename root for tidal input files
-!!$   CHARACTER(len= 4), PUBLIC, DIMENSION(jptides_max) ::   tide_cpt   !: Names of tidal components used.
-!!$
-!!$   INTEGER , PUBLIC, DIMENSION(jptides_max) ::   nindx        !: ???
-!!$   REAL(wp), PUBLIC, DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr)
-!!$   
-!!$   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   ssh1, ssh2   ! Tidal constituents : SSH
-!!$   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   u1  , u2     ! Tidal constituents : U
-!!$   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   v1  , v2     ! Tidal constituents : V
-!!$   
-!!$   !!----------------------------------------------------------------------
-!!$   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-!!$   !! $Id: obctides.F90 2528 2010-12-27 17:33:53Z rblod $ 
-!!$   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-!!$   !!----------------------------------------------------------------------
-!!$CONTAINS
-!!$
-!!$   SUBROUTINE tide_init
-!!$      !!----------------------------------------------------------------------
-!!$      !!                    ***  SUBROUTINE tide_init  ***
-!!$      !!                     
-!!$      !! ** Purpose : - Read in namelist for tides
-!!$      !!
-!!$      !!----------------------------------------------------------------------
-!!$      INTEGER ::   itide                  ! dummy loop index 
-!!$      !!
-!!$      NAMELIST/namobc_tide/ln_tide_date, filtide, tide_cpt, tide_speed
-!!$      !!----------------------------------------------------------------------
-!!$
-!!$      IF(lwp) WRITE(numout,*)
-!!$      IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries'
-!!$      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
-!!$
-!!$      ! Namelist namobc_tide : tidal harmonic forcing at open boundaries
-!!$      ln_tide_date = .false.
-!!$      filtide(:) = ''
-!!$      tide_speed(:) = 0.0
-!!$      tide_cpt(:) = ''
-!!$      REWIND( numnam )                                ! Read namelist parameters
-!!$      READ  ( numnam, namobc_tide )
-!!$      !                                               ! Count number of components specified
-!!$      ntide = jptides_max
-!!$      DO itide = 1, jptides_max
-!!$        IF( tide_cpt(itide) == '' ) THEN
-!!$           ntide = itide-1
-!!$           exit
-!!$        ENDIF
-!!$      END DO
-!!$
-!!$      !                                               ! find constituents in standard list
-!!$      DO itide = 1, ntide
-!!$         nindx(itide) = 0
-!!$         IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1
-!!$         IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2
-!!$         IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3
-!!$         IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4
-!!$         IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5
-!!$         IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6
-!!$         IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7
-!!$         IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8
-!!$         IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9
-!!$         IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10
-!!$         IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11
-!!$         IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12
-!!$         IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13
-!!$         IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14
-!!$         IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15
-!!$         IF( nindx(itide) == 0  .AND. lwp ) THEN
-!!$            WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list'
-!!$            CALL ctl_warn( ctmp1 )
-!!$         ENDIF
-!!$      END DO
-!!$      !                                               ! Parameter control and print
-!!$      IF( ntide < 1 ) THEN 
-!!$         CALL ctl_stop( '          Did not find any tidal components in namelist namobc_tide' )
-!!$      ELSE
-!!$         IF(lwp) WRITE(numout,*) '          Namelist namobc_tide : tidal harmonic forcing at open boundaries'
-!!$         IF(lwp) WRITE(numout,*) '             tidal components specified ', ntide
-!!$         IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:ntide)
-!!$         IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : '
-!!$         IF(lwp) WRITE(numout,*) '                ', tide_speed(1:ntide)
-!!$      ENDIF
-!!$
-!!$      ! Initialisation of tidal harmonics arrays
-!!$      sshtide(:) = 0.e0
-!!$      utide  (:) = 0.e0
-!!$      vtide  (:) = 0.e0
-!!$      !
-!!$   END SUBROUTINE tide_init
-!!$
-!!$
-!!$   SUBROUTINE tide_data
-!!$      !!----------------------------------------------------------------------
-!!$      !!                    ***  SUBROUTINE tide_data  ***
-!!$      !!                    
-!!$      !! ** Purpose : - Read in tidal harmonics data and adjust for the start 
-!!$      !!                time of the model run. 
-!!$      !!
-!!$      !!----------------------------------------------------------------------
-!!$      INTEGER ::   itide, igrd, ib        ! dummy loop indices
-!!$      CHARACTER(len=80) :: clfile         ! full file name for tidal input file 
-!!$      INTEGER ::   ipi, ipj, inum, idvar  ! temporary integers (netcdf read)
-!!$      INTEGER, DIMENSION(6) :: lendta=0   ! length of data in the file (note may be different from nblendta!)
-!!$      REAL(wp) ::  z_arg, z_atde, z_btde, z1t, z2t           
-!!$      REAL(wp), DIMENSION(jpbdta,1) ::   zdta   ! temporary array for data fields
-!!$      REAL(wp), DIMENSION(jptides_max) :: z_vplu, z_ftc
-!!$      !!------------------------------------------------------------------------------
-!!$
-!!$      ! Open files and read in tidal forcing data
-!!$      ! -----------------------------------------
-!!$
-!!$      ipj = 1
-!!$
-!!$      DO itide = 1, ntide
-!!$         !                                                              ! SSH fields
-!!$         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc'
-!!$         IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
-!!$         CALL iom_open( clfile, inum )
-!!$         igrd = 4
-!!$         IF( nblendta(igrd) <= 0 ) THEN 
-!!$            idvar = iom_varid( inum,'z1' )
-!!$            IF(lwp) WRITE(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar)
-!!$            nblendta(igrd) = iom_file(inum)%dimsz(1,idvar)
-!!$            WRITE(numout,*) 'Dim size for z1 is ', nblendta(igrd)
-!!$         ENDIF
-!!$         ipi = nblendta(igrd)
-!!$         CALL iom_get( inum, jpdom_unknown, 'z1', zdta(1:ipi,1:ipj) )
-!!$         DO ib = 1, nblenrim(igrd)
-!!$            ssh1(ib,itide) = zdta(nbmap(ib,igrd),1)
-!!$         END DO
-!!$         CALL iom_get( inum, jpdom_unknown, 'z2', zdta(1:ipi,1:ipj) )
-!!$         DO ib = 1, nblenrim(igrd)
-!!$            ssh2(ib,itide) = zdta(nbmap(ib,igrd),1)
-!!$         END DO
-!!$         CALL iom_close( inum )
-!!$         !
-!!$         !                                                              ! U fields
-!!$         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc'
-!!$         IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
-!!$         CALL iom_open( clfile, inum )
-!!$         igrd = 5
-!!$         IF( lendta(igrd) <= 0 ) THEN 
-!!$            idvar = iom_varid( inum,'u1' )
-!!$            lendta(igrd) = iom_file(inum)%dimsz(1,idvar)
-!!$            WRITE(numout,*) 'Dim size for u1 is ',lendta(igrd)
-!!$         ENDIF
-!!$         ipi = lendta(igrd)
-!!$         CALL iom_get( inum, jpdom_unknown, 'u1', zdta(1:ipi,1:ipj) )
-!!$         DO ib = 1, nblenrim(igrd)
-!!$            u1(ib,itide) = zdta(nbmap(ib,igrd),1)
-!!$         END DO
-!!$         CALL iom_get( inum, jpdom_unknown, 'u2', zdta(1:ipi,1:ipj) )
-!!$         DO ib = 1, nblenrim(igrd)
-!!$            u2(ib,itide) = zdta(nbmap(ib,igrd),1)
-!!$         END DO
-!!$         CALL iom_close( inum )
-!!$         !
-!!$         !                                                              ! V fields
-!!$         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc'
-!!$         if(lwp) write(numout,*) 'Reading data from file ', clfile
-!!$         CALL iom_open( clfile, inum )
-!!$         igrd = 6
-!!$         IF( lendta(igrd) <= 0 ) THEN 
-!!$            idvar = iom_varid( inum,'v1' )
-!!$            lendta(igrd) = iom_file(inum)%dimsz(1,idvar)
-!!$            WRITE(numout,*) 'Dim size for v1 is ', lendta(igrd)
-!!$         ENDIF
-!!$         ipi = lendta(igrd)
-!!$         CALL iom_get( inum, jpdom_unknown, 'v1', zdta(1:ipi,1:ipj) )
-!!$         DO ib = 1, nblenrim(igrd)
-!!$            v1(ib,itide) = zdta(nbmap(ib,igrd),1)
-!!$         END DO
-!!$         CALL iom_get( inum, jpdom_unknown, 'v2', zdta(1:ipi,1:ipj) )
-!!$         DO ib=1, nblenrim(igrd)
-!!$            v2(ib,itide) = zdta(nbmap(ib,igrd),1)
-!!$         END DO
-!!$         CALL iom_close( inum )
-!!$         !
-!!$      END DO ! end loop on tidal components
-!!$
-!!$      IF( ln_tide_date ) THEN      ! correct for date factors
-!!$
-!!$!! used nmonth, nyear and nday from daymod....
-!!$         ! Calculate date corrects for 15 standard consituents
-!!$         ! This is the initialisation step, so nday, nmonth, nyear are the 
-!!$         ! initial date/time of the integration.
-!!$           print *, nday,nmonth,nyear
-!!$           nyear  = int(ndate0 / 10000  )                           ! initial year
-!!$           nmonth = int((ndate0 - nyear * 10000 ) / 100 )          ! initial month
-!!$           nday   = int(ndate0 - nyear * 10000 - nmonth * 100)
-!!$
-!!$         CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu )
-!!$
-!!$         IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear
-!!$
-!!$         DO itide = 1, ntide       ! loop on tidal components
-!!$            !
-!!$            IF( nindx(itide) /= 0 ) THEN
-!!$!!gm use rpi  and rad global variable  
-!!$               z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0
-!!$               z_atde=z_ftc(nindx(itide))*cos(z_arg)
-!!$               z_btde=z_ftc(nindx(itide))*sin(z_arg)
-!!$               IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), tide_speed(itide),   &
-!!$                  &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide))
-!!$            ELSE
-!!$               z_atde = 1.0_wp
-!!$               z_btde = 0.0_wp
-!!$            ENDIF
-!!$            !                                         !  elevation         
-!!$            igrd = 4
-!!$            DO ib = 1, nblenrim(igrd)                
-!!$               z1t = z_atde * ssh1(ib,itide) + z_btde * ssh2(ib,itide)
-!!$               z2t = z_atde * ssh2(ib,itide) - z_btde * ssh1(ib,itide)
-!!$               ssh1(ib,itide) = z1t
-!!$               ssh2(ib,itide) = z2t
-!!$            END DO
-!!$            !                                         !  u       
-!!$            igrd = 5
-!!$            DO ib = 1, nblenrim(igrd)                
-!!$               z1t = z_atde * u1(ib,itide) + z_btde * u2(ib,itide)
-!!$               z2t = z_atde * u2(ib,itide) - z_btde * u1(ib,itide)
-!!$               u1(ib,itide) = z1t
-!!$               u2(ib,itide) = z2t
-!!$            END DO
-!!$            !                                         !  v       
-!!$            igrd = 6
-!!$            DO ib = 1, nblenrim(igrd)                
-!!$               z1t = z_atde * v1(ib,itide) + z_btde * v2(ib,itide)
-!!$               z2t = z_atde * v2(ib,itide) - z_btde * v1(ib,itide)
-!!$               v1(ib,itide) = z1t
-!!$               v2(ib,itide) = z2t
-!!$            END DO
-!!$            !
-!!$         END DO   ! end loop on tidal components
-!!$         !
-!!$      ENDIF ! date correction
-!!$      !
-!!$   END SUBROUTINE tide_data
-!!$
-!!$
-!!$   SUBROUTINE tide_update ( kt, jit )
-!!$      !!----------------------------------------------------------------------
-!!$      !!                 ***  SUBROUTINE tide_update  ***
-!!$      !!                
-!!$      !! ** Purpose : - Add tidal forcing to ssh_obc, ubt_obc and vbt_obc arrays. 
-!!$      !!                
-!!$      !!----------------------------------------------------------------------
-!!$      INTEGER, INTENT( in ) ::   kt      ! Main timestep counter
-!!$!!gm doctor jit ==> kit
-!!$      INTEGER, INTENT( in ) ::   jit     ! Barotropic timestep counter (for timesplitting option)
-!!$      !!
-!!$      INTEGER  ::   itide, igrd, ib      ! dummy loop indices
-!!$      REAL(wp) ::   z_arg, z_sarg            !            
-!!$      REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost
-!!$      !!----------------------------------------------------------------------
-!!$
-!!$      ! Note tide phase speeds are in deg/hour, so we need to convert the
-!!$      ! elapsed time in seconds to hours by dividing by 3600.0
-!!$      IF( jit == 0 ) THEN  
-!!$         z_arg = kt * rdt * rad / 3600.0
-!!$      ELSE                              ! we are in a barotropic subcycle (for timesplitting option)
-!!$!         z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,lwp) ) * rad / 3600.0
-!!$         z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,wp) ) * rad / 3600.0
-!!$      ENDIF
-!!$
-!!$      DO itide = 1, ntide
-!!$         z_sarg = z_arg * tide_speed(itide)
-!!$         z_cost(itide) = COS( z_sarg )
-!!$         z_sist(itide) = SIN( z_sarg )
-!!$      END DO
-!!$
-!!$      ! summing of tidal constituents into BDY arrays
-!!$      sshtide(:) = 0.0
-!!$      utide (:) = 0.0
-!!$      vtide (:) = 0.0
-!!$      !
-!!$      DO itide = 1, ntide
-!!$         igrd=4                              ! SSH on tracer grid.
-!!$         DO ib = 1, nblenrim(igrd)
-!!$            sshtide(ib) =sshtide(ib)+ ssh1(ib,itide)*z_cost(itide) + ssh2(ib,itide)*z_sist(itide)
-!!$            !    if(lwp) write(numout,*) 'z',ib,itide,sshtide(ib), ssh1(ib,itide),ssh2(ib,itide)
-!!$         END DO
-!!$         igrd=5                              ! U grid
-!!$         DO ib=1, nblenrim(igrd)
-!!$            utide(ib) = utide(ib)+ u1(ib,itide)*z_cost(itide) + u2(ib,itide)*z_sist(itide)
-!!$            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), u1(ib,itide),u2(ib,itide)
-!!$         END DO
-!!$         igrd=6                              ! V grid
-!!$         DO ib=1, nblenrim(igrd)
-!!$            vtide(ib) = vtide(ib)+ v1(ib,itide)*z_cost(itide) + v2(ib,itide)*z_sist(itide)
-!!$            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), v1(ib,itide),v2(ib,itide)
-!!$         END DO
-!!$      END DO
-!!$      !
-!!$   END SUBROUTINE tide_update
-!!$
-!!$!!gm  doctor naming of dummy argument variables!!!   and all local variables
-!!$
-!!$   SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu )
-!!$      !!----------------------------------------------------------------------
-!!$      !!                   ***  SUBROUTINE uvset  ***
-!!$      !!                     
-!!$      !! ** Purpose : - adjust tidal forcing for date factors
-!!$      !!                
-!!$      !!----------------------------------------------------------------------
-!!$      implicit none
-!!$      INTEGER, INTENT( in ) ::   ihs     ! Start time hours
-!!$      INTEGER, INTENT( in ) ::   iday    ! start time days
-!!$      INTEGER, INTENT( in ) ::   imnth   ! start time month
-!!$      INTEGER, INTENT( in ) ::   iyr     ! start time year
-!!$      !!
-!!$!!gm  nc is jptides_max ????
-!!$      INTEGER , PARAMETER     ::  nc =15    ! maximum number of constituents
-!!$      CHARACTER(len=8), DIMENSION(nc) :: cname
-!!$      INTEGER                 ::   year, vd, ivdy, ndc, i, k
-!!$      REAL(wp)                ::   ss, h, p, en, p1, rtd
-!!$      REAL(wp), DIMENSION(nc) ::   f                          ! nodal correction
-!!$      REAL(wp), DIMENSION(nc) ::   z_vplu                     ! phase correction
-!!$      REAL(wp), DIMENSION(nc) ::   u, v, zig
-!!$      !!
-!!$      DATA cname/  'q1'    ,    'o1'    ,     'p1'   ,    's1'    ,     'k1'    ,   &
-!!$         &         '2n2'   ,    'mu2'   ,     'n2'   ,    'nu2'   ,     'm2'    ,   &
-!!$         &         'l2'    ,    't2'    ,     's2'   ,    'k2'    ,     'm4'      /
-!!$      DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878,  .2625161701,   &
-!!$         &      .4868657873, .4881373225, .4963669182, .4976384533,  .5058680490,   &
-!!$         &      .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098   /
-!!$      !!----------------------------------------------------------------------
-!!$!
-!!$!     ihs             -  start time gmt on ...
-!!$!     iday/imnth/iyr  -  date   e.g.   12/10/87
-!!$!
-!!$      CALL vday(iday,imnth,iyr,ivdy)
-!!$      vd = ivdy
-!!$!
-!!$!rp   note change of year number for d. blackman shpen
-!!$!rp   if(iyr.ge.1000) year=iyr-1900
-!!$!rp   if(iyr.lt.1000) year=iyr
-!!$      year = iyr
-!!$!
-!!$!.....year  =  year of required data
-!!$!.....vd    =  day of required data..set up for 0000gmt day year
-!!$      ndc = nc
-!!$!.....ndc   =  number of constituents allowed
-!!$!!gm use rpi ?
-!!$      rtd = 360.0 / 6.2831852
-!!$      DO i = 1, ndc
-!!$         zig(i) = zig(i)*rtd
-!!$         ! sigo(i)= zig(i)
-!!$      END DO
-!!$
-!!$!!gm try to avoid RETURN  in F90
-!!$      IF( year == 0 )   RETURN
-!!$      
-!!$         CALL shpen( year, vd, ss, h , p , en, p1 )
-!!$         CALL ufset( p   , en, u , f )
-!!$         CALL vset ( ss  , h , p , en, p1, v )
-!!$         !
-!!$         DO k = 1, nc
-!!$            z_vplu(k) = v(k) + u(k)
-!!$            z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k)
-!!$            DO WHILE( z_vplu(k) < 0    )
-!!$               z_vplu(k) = z_vplu(k) + 360.0
-!!$            END DO
-!!$            DO WHILE( z_vplu(k) > 360. )
-!!$               z_vplu(k) = z_vplu(k) - 360.0
-!!$            END DO
-!!$         END DO
-!!$         !
-!!$      END SUBROUTINE uvset
-!!$
-!!$
-!!$      SUBROUTINE vday( iday, imnth, iy, ivdy )
-!!$      !!----------------------------------------------------------------------
-!!$      !!                   *** SUBROUTINE vday  ***
-!!$      !!                  
-!!$      !! ** Purpose : - adjust tidal forcing for date factors
-!!$      !!                
-!!$      !!----------------------------------------------------------------------
-!!$      INTEGER, INTENT(in   ) ::   iday, imnth, iy   ! ????
-!!$      INTEGER, INTENT(  out) ::   ivdy              ! ???
-!!$      !! 
-!!$      INTEGER ::   iyr
-!!$      !!------------------------------------------------------------------------------
-!!$
-!!$!!gm   nday_year in day mode is the variable compiuted here, no?
-!!$!!gm      nday_year ,   &  !: curent day counted from jan 1st of the current year
-!!$
-!!$      !calculate day number in year from day/month/year
-!!$       if(imnth.eq.1) ivdy=iday
-!!$       if(imnth.eq.2) ivdy=iday+31
-!!$       if(imnth.eq.3) ivdy=iday+59
-!!$       if(imnth.eq.4) ivdy=iday+90
-!!$       if(imnth.eq.5) ivdy=iday+120
-!!$       if(imnth.eq.6) ivdy=iday+151
-!!$       if(imnth.eq.7) ivdy=iday+181
-!!$       if(imnth.eq.8) ivdy=iday+212
-!!$       if(imnth.eq.9) ivdy=iday+243
-!!$       if(imnth.eq.10) ivdy=iday+273
-!!$       if(imnth.eq.11) ivdy=iday+304
-!!$       if(imnth.eq.12) ivdy=iday+334
-!!$       iyr=iy
-!!$       if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1
-!!$       if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1
-!!$       if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1
-!!$      !
-!!$   END SUBROUTINE vday
-!!$
-!!$!!doctor norme for dummy arguments
-!!$
-!!$   SUBROUTINE shpen( year, vd, s, h, p, en, p1 )
-!!$      !!----------------------------------------------------------------------
-!!$      !!                    ***  SUBROUTINE shpen  ***
-!!$      !!                   
-!!$      !! ** Purpose : - calculate astronomical arguments for tides
-!!$      !!                this version from d. blackman 30 nove 1990
-!!$      !!
-!!$      !!----------------------------------------------------------------------
-!!$!!gm add INTENT in, out or inout....    DOCTOR name....
-!!$!!gm please do not use variable name with a single letter:  impossible to search in a code
-!!$      INTEGER  ::   year,vd
-!!$      REAL(wp) ::   s,h,p,en,p1
-!!$      !!
-!!$      INTEGER  ::   yr,ilc,icent,it,iday,ild,ipos,nn,iyd
-!!$      REAL(wp) ::   cycle,t,td,delt(84),delta,deltat
-!!$      !!
-!!$      DATA delt /-5.04, -3.90, -2.87, -0.58,  0.71,  1.80,           &
-!!$         &        3.08,  4.63,  5.86,  7.21,  8.58, 10.50, 12.10,    &
-!!$         &       12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39,    &
-!!$         &       19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68,    &
-!!$         &       22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63,    &
-!!$         &       23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80,    &
-!!$         &       24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89,    &
-!!$         &       27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96,    &
-!!$         &       31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39,    &
-!!$         &       33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87,    &
-!!$         &       38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00,    &
-!!$         &       45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81,    &
-!!$         &       52.57 /
-!!$      !!----------------------------------------------------------------------
-!!$
-!!$      cycle = 360.0
-!!$      ilc = 0
-!!$      icent = year / 100
-!!$      yr = year - icent * 100
-!!$      t = icent - 20
-!!$!
-!!$!     for the following equations
-!!$!     time origin is fixed at 00 hr of jan 1st,2000.
-!!$!     see notes by cartwright
-!!$!
-!!$!!gm  old coding style, use CASE instead  and avoid GOTO (obsolescence in fortran 90)
-!!$!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
-!!$      it = icent - 20
-!!$      if (it) 1,2,2
-!!$    1 iday = it/4 -it
-!!$      go to 3
-!!$    2 iday = (it+3)/4 - it
-!!$!
-!!$!     t is in julian century
-!!$!     correction in gegorian calander where only century year divisible
-!!$!     by 4 is leap year.
-!!$!
-!!$    3 continue
-!!$!
-!!$      td = 0.0
-!!$!
-!!$!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
-!!$      if (yr) 4,5,4
-!!$!
-!!$    4 iyd = 365*yr
-!!$      ild = (yr-1)/4
-!!$      if((icent - (icent/4)*4) .eq. 0) ilc = 1
-!!$      td = iyd + ild + ilc
-!!$!
-!!$    5 td = td + iday + vd -1.0 - 0.5
-!!$      t = t + (td/36525.0)
-!!$!
-!!$      ipos=year-1899
-!!$      if (ipos .lt. 0) go to 7
-!!$      if (ipos .gt. 83) go to 6
-!!$!
-!!$      delta = (delt(ipos+1)+delt(ipos))/2.0
-!!$      go to 7
-!!$!
-!!$    6 delta= (65.0-50.5)/20.0*(year-1980)+50.5
-!!$!
-!!$    7 deltat = delta * 1.0e-6
-!!$!
-!!$!!gm   precision of the computation   :  example for s it should be replace by:
-!!$!!gm  s   = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat   ==>  more precise  modify the last digits results
-!!$      s   = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat
-!!$      h   = 280.4661 + 36000.7698 *t + 0.0003*t*t +  11.0*deltat
-!!$      p   =  83.3535 + 4069.0139  *t - 0.0103*t*t +       deltat
-!!$      en  = 234.9555 + 1934.1363  *t - 0.0021*t*t +       deltat
-!!$      p1  = 282.9384 + 1.7195     *t + 0.0005*t*t
-!!$      !
-!!$      nn = s / cycle
-!!$      s = s - nn * cycle
-!!$      IF( s < 0.e0 )   s = s + cycle
-!!$      !
-!!$      nn = h / cycle
-!!$      h = h - cycle * nn
-!!$      IF( h < 0.e0 )   h = h + cycle
-!!$      !
-!!$      nn = p / cycle
-!!$      p = p - cycle * nn
-!!$      IF( p < 0.e0)   p = p + cycle
-!!$      !
-!!$      nn = en / cycle
-!!$      en = en - cycle * nn
-!!$      IF( en < 0.e0 )   en = en + cycle
-!!$      en = cycle - en
-!!$      !
-!!$      nn = p1 / cycle
-!!$      p1 = p1 - nn * cycle
-!!$      !
-!!$   END SUBROUTINE shpen
-!!$
-!!$
-!!$   SUBROUTINE ufset( p, cn, b, a )
-!!$      !!----------------------------------------------------------------------
-!!$      !!                    ***  SUBROUTINE ufset  ***
-!!$      !!                    
-!!$      !! ** Purpose : - calculate nodal parameters for the tides
-!!$      !!                
-!!$      !!----------------------------------------------------------------------
-!!$!!gm  doctor naming of dummy argument variables!!!   and all local variables
-!!$!!gm  nc is jptides_max ????
-!!$      integer nc
-!!$      parameter (nc=15)
-!!$      REAL(wp) p,cn
-!!$      !!
-!!$      
-!!$!!gm  rad is already a public variable defined in phycst.F90 ....   ==> doctor norme  local real start with "z"
-!!$      REAL(wp) ::   w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad
-!!$      REAL(wp) ::   a(nc), b(nc)
-!!$      INTEGER  ::   ndc, k
-!!$      !!----------------------------------------------------------------------
-!!$
-!!$      ndc = nc
-!!$
-!!$!    a=f       ,  b =u
-!!$!    t is  zero as compared to tifa.
-!!$!! use rad defined in phycst   (i.e.  add a USE phycst at the begining of the module
-!!$      rad = 6.2831852d0/360.0
-!!$      pw = p  * rad
-!!$      nw = cn * rad
-!!$      w1 = cos(   nw )
-!!$      w2 = cos( 2*nw )
-!!$      w3 = cos( 3*nw )
-!!$      w4 = sin(   nw )
-!!$      w5 = sin( 2*nw )
-!!$      w6 = sin( 3*nw )
-!!$      w7 = 1. - 0.2505 * COS( 2*pw      ) - 0.1102 * COS(2*pw-nw )   &
-!!$         &    - 0.156  * COS( 2*pw-2*nw ) - 0.037  * COS(     nw )
-!!$      w8 =    - 0.2505 * SIN( 2*pw      ) - 0.1102 * SIN(2*pw-nw )   &
-!!$         &    - 0.0156 * SIN( 2*pw-2*nw ) - 0.037  * SIN(     nw )
-!!$      !
-!!$      a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3
-!!$      b(1) =          0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6
-!!$      !   q1
-!!$      a(2) = a(1)
-!!$      b(2) = b(1)
-!!$      !   o1
-!!$      a(3) = 1.0
-!!$      b(3) = 0.0
-!!$      !   p1
-!!$      a(4) = 1.0
-!!$      b(4) = 0.0
-!!$      !   s1
-!!$      a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3
-!!$      b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6
-!!$      !   k1
-!!$      a(6) =1.0004 -0.0373*w1+ 0.0002*w2
-!!$      b(6) = -0.0374*w4
-!!$      !  2n2
-!!$      a(7) = a(6)
-!!$      b(7) = b(6)
-!!$      !  mu2
-!!$      a(8) = a(6)
-!!$      b(8) = b(6)
-!!$      !   n2
-!!$      a(9) = a(6)
-!!$      b(9) = b(6)
-!!$      !  nu2
-!!$      a(10) = a(6)
-!!$      b(10) = b(6)
-!!$      !   m2
-!!$      a(11) = SQRT( w7 * w7 + w8 * w8 )
-!!$      b(11) = ATAN( w8 / w7 )
-!!$!!gmuse rpi  instead of 3.141992 ???   true pi is rpi=3.141592653589793_wp  .....   ????
-!!$      IF( w7 < 0.e0 )   b(11) = b(11) + 3.141992
-!!$      !   l2
-!!$      a(12) = 1.0
-!!$      b(12) = 0.0
-!!$      !   t2
-!!$      a(13)= a(12)
-!!$      b(13)= b(12)
-!!$      !   s2
-!!$      a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3
-!!$      b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6
-!!$      !   k2
-!!$      a(15) = a(6)*a(6)
-!!$      b(15) = 2*b(6)
-!!$      !   m4
-!!$!!gm  old coding,  remove GOTO and label of lines
-!!$!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
-!!$      DO 40 k = 1,ndc
-!!$         b(k) = b(k)/rad
-!!$32       if (b(k)) 34,35,35
-!!$34       b(k) = b(k) + 360.0
-!!$         go to 32
-!!$35       if (b(k)-360.0) 40,37,37
-!!$37       b(k) = b(k)-360.0
-!!$         go to 35
-!!$40    continue
-!!$      !
-!!$   END SUBROUTINE ufset
-!!$   
-!!$
-!!$   SUBROUTINE vset( s,h,p,en,p1,v)
-!!$      !!----------------------------------------------------------------------
-!!$      !!                    ***  SUBROUTINE vset  ***
-!!$      !!                   
-!!$      !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run
-!!$      !!                
-!!$      !!----------------------------------------------------------------------
-!!$!!gm  doctor naming of dummy argument variables!!!   and all local variables
-!!$!!gm  nc is jptides_max ????
-!!$!!gm  en argument is not used: suppress it ?
-!!$      integer nc
-!!$      parameter (nc=15)
-!!$      real(wp) s,h,p,en,p1
-!!$      real(wp)   v(nc)
-!!$      !!
-!!$      integer ndc, k
-!!$      !!----------------------------------------------------------------------
-!!$
-!!$      ndc = nc
-!!$      !   v s  are computed here.
-!!$      v(1) =-3*s +h +p +270      ! Q1
-!!$      v(2) =-2*s +h +270.0       ! O1
-!!$      v(3) =-h +270              ! P1
-!!$      v(4) =180                  ! S1
-!!$      v(5) =h +90.0              ! K1
-!!$      v(6) =-4*s +2*h +2*p       ! 2N2
-!!$      v(7) =-4*(s-h)             ! MU2
-!!$      v(8) =-3*s +2*h +p         ! N2
-!!$      v(9) =-3*s +4*h -p         ! MU2
-!!$      v(10) =-2*s +2*h           ! M2
-!!$      v(11) =-s +2*h -p +180     ! L2 
-!!$      v(12) =-h +p1              ! T2
-!!$      v(13) =0                   ! S2
-!!$      v(14) =h+h                 ! K2
-!!$      v(15) =2*v(10)             ! M4
-!!$!
-!!$!!gm  old coding,  remove GOTO and label of lines
-!!$!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
-!!$      do 72 k = 1, ndc
-!!$69       if( v(k) )   70,71,71
-!!$70       v(k) = v(k)+360.0
-!!$         go to 69
-!!$71       if( v(k) - 360.0 )   72,73,73
-!!$73       v(k) = v(k)-360.0
-!!$         go to 71
-!!$72    continue
-!!$      !
-!!$   END SUBROUTINE vset
-!!$
-!!$#else
+   !!----------------------------------------------------------------------
+   !!   'key_obc'     Open Boundary Condition
+   !!----------------------------------------------------------------------
+   !!   PUBLIC
+   !!      tide_init     : read of namelist and initialisation of tidal harmonics data
+   !!      tide_update   : calculation of tidal forcing at each timestep
+   !!   PRIVATE
+   !!      uvset         :\
+   !!      vday          : | Routines to correct tidal harmonics forcing for
+   !!      shpen         : | start time of integration
+   !!      ufset         : |
+   !!      vset          :/
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers 
+   USE dom_oce         ! ocean space and time domain
+   USE iom
+   USE in_out_manager  ! I/O units
+   USE phycst          ! physical constants
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE obc_par         ! Unstructured boundary parameters
+   USE obc_oce         ! ocean open boundary conditions
+   USE daymod          ! calendar
+   USE fldread, ONLY: fld_map
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   tide_init     ! routine called in nemo_init
+   PUBLIC   tide_update   ! routine called in obcdyn
+
+   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
+      INTEGER                                ::   ncpt       !: Actual number of tidal components
+      REAL(wp), POINTER, DIMENSION(:)        ::   speed      !: Phase speed of tidal constituent (deg/hr)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V
+   END TYPE TIDES_DATA
+
+   INTEGER, PUBLIC, PARAMETER                  ::   jptides_max = 15      !: Max number of tidal contituents
+
+   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_obc), TARGET ::   tides                 !: External tidal harmonics data
+   
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: obctides.F90 2528 2010-12-27 17:33:53Z rblod $ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE tide_init
+      !!----------------------------------------------------------------------
+      !!                    ***  SUBROUTINE tide_init  ***
+      !!                     
+      !! ** Purpose : - Read in namelist for tides and initialise external
+      !!                tidal harmonics data
+      !!
+      !!----------------------------------------------------------------------
+      !! namelist variables
+      !!-------------------
+      LOGICAL                                   ::   ln_tide_date !: =T correct tide phases and amplitude for model start date
+      CHARACTER(len=80)                         ::   filtide      !: Filename root for tidal input files
+      CHARACTER(len= 4), DIMENSION(jptides_max) ::   tide_cpt     !: Names of tidal components used.
+      REAL(wp),          DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr)
+      !!
+      INTEGER, DIMENSION(jptides_max)           ::   nindx              !: index of pre-set tidal components
+      INTEGER                                   ::   ib_obc, itide, ib  !: dummy loop indices
+      INTEGER                                   ::   inum, igrd
+      INTEGER, DIMENSION(3)                     ::   ilen0              !: length of boundary data (from OBC arrays)
+      CHARACTER(len=80)                         ::   clfile             !: full file name for tidal input file 
+      REAL(wp)                                  ::   z_arg, z_atde, z_btde, z1t, z2t           
+      REAL(wp),DIMENSION(jptides_max)           ::   z_vplu, z_ftc
+      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read           !: work space to read in tidal harmonics data
+      !!
+      TYPE(TIDES_DATA),  POINTER                ::   td                 !: local short cut   
+      !!
+      NAMELIST/namobc_tide/ln_tide_date, filtide, tide_cpt, tide_speed
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
+
+      REWIND(numnam)
+      DO ib_obc = 1, nb_obc
+         IF( nn_tides(ib_obc) .gt. 0 ) THEN
+
+            td => tides(ib_obc)
+
+            ! Namelist namobc_tide : tidal harmonic forcing at open boundaries
+            ln_tide_date = .false.
+            filtide(:) = ''
+            tide_speed(:) = 0.0
+            tide_cpt(:) = ''
+
+            ! Don't REWIND here - may need to read more than one of these namelists.
+            READ  ( numnam, namobc_tide )
+            !                                               ! Count number of components specified
+            td%ncpt = 0
+            DO itide = 1, jptides_max
+              IF( tide_cpt(itide) /= '' ) THEN
+                 td%ncpt = td%ncpt + 1
+              ENDIF
+            END DO
+
+            ! Fill in phase speeds from namelist
+            ALLOCATE( td%speed(td%ncpt) )
+            td%speed = tide_speed(1:td%ncpt)
+
+            ! Find constituents in standard list
+            DO itide = 1, td%ncpt
+               nindx(itide) = 0
+               IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1
+               IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2
+               IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3
+               IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4
+               IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5
+               IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6
+               IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7
+               IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8
+               IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9
+               IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10
+               IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11
+               IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12
+               IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13
+               IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14
+               IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15
+               IF( nindx(itide) == 0  .AND. lwp ) THEN
+                  WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list'
+                  CALL ctl_warn( ctmp1 )
+               ENDIF
+            END DO
+            !                                               ! Parameter control and print
+            IF( td%ncpt < 1 ) THEN 
+               CALL ctl_stop( '          Did not find any tidal components in namelist namobc_tide' )
+            ELSE
+               IF(lwp) WRITE(numout,*) '          Namelist namobc_tide : tidal harmonic forcing at open boundaries'
+               IF(lwp) WRITE(numout,*) '             tidal components specified ', td%ncpt
+               IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:td%ncpt)
+               IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : '
+               IF(lwp) WRITE(numout,*) '                ', tide_speed(1:td%ncpt)
+            ENDIF
+
+
+            ! Allocate space for tidal harmonics data - 
+            ! get size from OBC data arrays
+            ! ---------------------------------------
+
+            ilen0(1) = SIZE( dta_obc(ib_obc)%ssh ) 
+            ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) )
+
+            ilen0(2) = SIZE( dta_obc(ib_obc)%u2d ) 
+            ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) )
+
+            ilen0(3) = SIZE( dta_obc(ib_obc)%v2d ) 
+            ALLOCATE( td%v( ilen0(3), td%ncpt, 2 ) )
+
+            ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) )
+
+
+            ! Open files and read in tidal forcing data
+            ! -----------------------------------------
+
+            DO itide = 1, td%ncpt
+               !                                                              ! SSH fields
+               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc'
+               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
+               CALL iom_open( clfile, inum )
+               CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,1) )
+               td%ssh(:,itide,1) = dta_read(1:ilen0(1),1,1)
+               CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,1) )
+               td%ssh(:,itide,2) = dta_read(1:ilen0(1),1,1)
+               CALL iom_close( inum )
+               !                                                              ! U fields
+               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc'
+               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
+               CALL iom_open( clfile, inum )
+               CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,2) )
+               td%u(:,itide,1) = dta_read(1:ilen0(2),1,1)
+               CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,2) )
+               td%u(:,itide,2) = dta_read(1:ilen0(2),1,1)
+               CALL iom_close( inum )
+               !                                                              ! V fields
+               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc'
+               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
+               CALL iom_open( clfile, inum )
+               CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,3) )
+               td%v(:,itide,1) = dta_read(1:ilen0(3),1,1)
+               CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,3) )
+               td%v(:,itide,2) = dta_read(1:ilen0(3),1,1)
+               CALL iom_close( inum )
+               !
+            END DO ! end loop on tidal components
+
+            IF( ln_tide_date ) THEN      ! correct for date factors
+
+!! used nmonth, nyear and nday from daymod....
+               ! Calculate date corrects for 15 standard consituents
+               ! This is the initialisation step, so nday, nmonth, nyear are the 
+               ! initial date/time of the integration.
+                 print *, nday,nmonth,nyear
+                 nyear  = int(ndate0 / 10000  )                          ! initial year
+                 nmonth = int((ndate0 - nyear * 10000 ) / 100 )          ! initial month
+                 nday   = int(ndate0 - nyear * 10000 - nmonth * 100)
+
+               CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu )
+
+               IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear
+
+               DO itide = 1, td%ncpt       ! loop on tidal components
+                  !
+                  IF( nindx(itide) /= 0 ) THEN
+!!gm use rpi  and rad global variable  
+                     z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0
+                     z_atde=z_ftc(nindx(itide))*cos(z_arg)
+                     z_btde=z_ftc(nindx(itide))*sin(z_arg)
+                     IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), td%speed(itide),   &
+                        &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide))
+                  ELSE
+                     z_atde = 1.0_wp
+                     z_btde = 0.0_wp
+                  ENDIF
+                  !                                         !  elevation         
+                  igrd = 1
+                  DO ib = 1, ilen0(igrd)                
+                     z1t = z_atde * td%ssh(ib,itide,1) + z_btde * td%ssh(ib,itide,2)
+                     z2t = z_atde * td%ssh(ib,itide,2) - z_btde * td%ssh(ib,itide,1)
+                     td%ssh(ib,itide,1) = z1t
+                     td%ssh(ib,itide,2) = z2t
+                  END DO
+                  !                                         !  u       
+                  igrd = 2
+                  DO ib = 1, ilen0(igrd)                
+                     z1t = z_atde * td%u(ib,itide,1) + z_btde * td%u(ib,itide,2)
+                     z2t = z_atde * td%u(ib,itide,2) - z_btde * td%u(ib,itide,1)
+                     td%u(ib,itide,1) = z1t
+                     td%u(ib,itide,2) = z2t
+                  END DO
+                  !                                         !  v       
+                  igrd = 3
+                  DO ib = 1, ilen0(igrd)                
+                     z1t = z_atde * td%v(ib,itide,1) + z_btde * td%v(ib,itide,2)
+                     z2t = z_atde * td%v(ib,itide,2) - z_btde * td%v(ib,itide,1)
+                     td%v(ib,itide,1) = z1t
+                     td%v(ib,itide,2) = z2t
+                  END DO
+                  !
+               END DO   ! end loop on tidal components
+               !
+            ENDIF ! date correction
+            !
+         ENDIF ! nn_tides .gt. 0
+         !
+      END DO ! loop on ib_obc
+
+   END SUBROUTINE tide_init
+
+
+   SUBROUTINE tide_update ( kt, jit, idx, dta, td )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE tide_update  ***
+      !!                
+      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays. 
+      !!                
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in )          ::   kt      ! Main timestep counter
+!!gm doctor jit ==> kit
+      INTEGER, INTENT( in )          ::   jit     ! Barotropic timestep counter (for timesplitting option)
+      TYPE(OBC_INDEX), INTENT( in )  ::   idx     ! OBC indices
+      TYPE(OBC_DATA),  INTENT(inout) ::   dta     ! OBC external data
+      TYPE(TIDES_DATA),INTENT( in )  ::   td      ! tidal harmonics data
+      !!
+      INTEGER                          :: itide, igrd, ib      ! dummy loop indices
+      REAL(wp)                         :: z_arg, z_sarg      
+      REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost
+      !!----------------------------------------------------------------------
+
+      ! Note tide phase speeds are in deg/hour, so we need to convert the
+      ! elapsed time in seconds to hours by dividing by 3600.0
+      IF( jit == 0 ) THEN  
+         z_arg = kt * rdt * rad / 3600.0
+      ELSE                              ! we are in a barotropic subcycle (for timesplitting option)
+!         z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,lwp) ) * rad / 3600.0
+         z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,wp) ) * rad / 3600.0
+      ENDIF
+
+      DO itide = 1, td%ncpt
+         z_sarg = z_arg * td%speed(itide)
+         z_cost(itide) = COS( z_sarg )
+         z_sist(itide) = SIN( z_sarg )
+      END DO
+
+      !
+      DO itide = 1, td%ncpt
+         igrd=1                              ! SSH on tracer grid.
+         DO ib = 1, idx%nblenrim(igrd)
+            dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)
+            !    if(lwp) write(numout,*) 'z', ib, itide, dta%ssh(ib), td%ssh(ib,itide,1),td%ssh(ib,itide,2)
+         END DO
+         igrd=2                              ! U grid
+         DO ib=1, idx%nblenrim(igrd)
+            dta%u2d(ib) = dta%u2d(ib) + td%u(ib,itide,1)*z_cost(itide) + td%u(ib,itide,2)*z_sist(itide)
+            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), td%u(ib,itide,1),td%u(ib,itide,2)
+         END DO
+         igrd=3                              ! V grid
+         DO ib=1, idx%nblenrim(igrd)
+            dta%v2d(ib) = dta%v2d(ib) + td%v(ib,itide,1)*z_cost(itide) + td%v(ib,itide,2)*z_sist(itide)
+            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), td%v(ib,itide,1),td%v(ib,itide,2)
+         END DO
+      END DO
+      !
+   END SUBROUTINE tide_update
+
+!!gm  doctor naming of dummy argument variables!!!   and all local variables
+
+   SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu )
+      !!----------------------------------------------------------------------
+      !!                   ***  SUBROUTINE uvset  ***
+      !!                     
+      !! ** Purpose : - adjust tidal forcing for date factors
+      !!                
+      !!----------------------------------------------------------------------
+      implicit none
+      INTEGER, INTENT( in ) ::   ihs     ! Start time hours
+      INTEGER, INTENT( in ) ::   iday    ! start time days
+      INTEGER, INTENT( in ) ::   imnth   ! start time month
+      INTEGER, INTENT( in ) ::   iyr     ! start time year
+      !!
+!!gm  nc is jptides_max ????
+      INTEGER , PARAMETER     ::  nc =15    ! maximum number of constituents
+      CHARACTER(len=8), DIMENSION(nc) :: cname
+      INTEGER                 ::   year, vd, ivdy, ndc, i, k
+      REAL(wp)                ::   ss, h, p, en, p1, rtd
+      REAL(wp), DIMENSION(nc) ::   f                          ! nodal correction
+      REAL(wp), DIMENSION(nc) ::   z_vplu                     ! phase correction
+      REAL(wp), DIMENSION(nc) ::   u, v, zig
+      !!
+      DATA cname/  'q1'    ,    'o1'    ,     'p1'   ,    's1'    ,     'k1'    ,   &
+         &         '2n2'   ,    'mu2'   ,     'n2'   ,    'nu2'   ,     'm2'    ,   &
+         &         'l2'    ,    't2'    ,     's2'   ,    'k2'    ,     'm4'      /
+      DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878,  .2625161701,   &
+         &      .4868657873, .4881373225, .4963669182, .4976384533,  .5058680490,   &
+         &      .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098   /
+      !!----------------------------------------------------------------------
+!
+!     ihs             -  start time gmt on ...
+!     iday/imnth/iyr  -  date   e.g.   12/10/87
+!
+      CALL vday(iday,imnth,iyr,ivdy)
+      vd = ivdy
+!
+!rp   note change of year number for d. blackman shpen
+!rp   if(iyr.ge.1000) year=iyr-1900
+!rp   if(iyr.lt.1000) year=iyr
+      year = iyr
+!
+!.....year  =  year of required data
+!.....vd    =  day of required data..set up for 0000gmt day year
+      ndc = nc
+!.....ndc   =  number of constituents allowed
+!!gm use rpi ?
+      rtd = 360.0 / 6.2831852
+      DO i = 1, ndc
+         zig(i) = zig(i)*rtd
+         ! sigo(i)= zig(i)
+      END DO
+
+!!gm try to avoid RETURN  in F90
+      IF( year == 0 )   RETURN
+      
+         CALL shpen( year, vd, ss, h , p , en, p1 )
+         CALL ufset( p   , en, u , f )
+         CALL vset ( ss  , h , p , en, p1, v )
+         !
+         DO k = 1, nc
+            z_vplu(k) = v(k) + u(k)
+            z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k)
+            DO WHILE( z_vplu(k) < 0    )
+               z_vplu(k) = z_vplu(k) + 360.0
+            END DO
+            DO WHILE( z_vplu(k) > 360. )
+               z_vplu(k) = z_vplu(k) - 360.0
+            END DO
+         END DO
+         !
+      END SUBROUTINE uvset
+
+
+      SUBROUTINE vday( iday, imnth, iy, ivdy )
+      !!----------------------------------------------------------------------
+      !!                   *** SUBROUTINE vday  ***
+      !!                  
+      !! ** Purpose : - adjust tidal forcing for date factors
+      !!                
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in   ) ::   iday, imnth, iy   ! ????
+      INTEGER, INTENT(  out) ::   ivdy              ! ???
+      !! 
+      INTEGER ::   iyr
+      !!------------------------------------------------------------------------------
+
+!!gm   nday_year in day mode is the variable compiuted here, no?
+!!gm      nday_year ,   &  !: curent day counted from jan 1st of the current year
+
+      !calculate day number in year from day/month/year
+       if(imnth.eq.1) ivdy=iday
+       if(imnth.eq.2) ivdy=iday+31
+       if(imnth.eq.3) ivdy=iday+59
+       if(imnth.eq.4) ivdy=iday+90
+       if(imnth.eq.5) ivdy=iday+120
+       if(imnth.eq.6) ivdy=iday+151
+       if(imnth.eq.7) ivdy=iday+181
+       if(imnth.eq.8) ivdy=iday+212
+       if(imnth.eq.9) ivdy=iday+243
+       if(imnth.eq.10) ivdy=iday+273
+       if(imnth.eq.11) ivdy=iday+304
+       if(imnth.eq.12) ivdy=iday+334
+       iyr=iy
+       if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1
+       if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1
+       if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1
+      !
+   END SUBROUTINE vday
+
+!!doctor norme for dummy arguments
+
+   SUBROUTINE shpen( year, vd, s, h, p, en, p1 )
+      !!----------------------------------------------------------------------
+      !!                    ***  SUBROUTINE shpen  ***
+      !!                   
+      !! ** Purpose : - calculate astronomical arguments for tides
+      !!                this version from d. blackman 30 nove 1990
+      !!
+      !!----------------------------------------------------------------------
+!!gm add INTENT in, out or inout....    DOCTOR name....
+!!gm please do not use variable name with a single letter:  impossible to search in a code
+      INTEGER  ::   year,vd
+      REAL(wp) ::   s,h,p,en,p1
+      !!
+      INTEGER  ::   yr,ilc,icent,it,iday,ild,ipos,nn,iyd
+      REAL(wp) ::   cycle,t,td,delt(84),delta,deltat
+      !!
+      DATA delt /-5.04, -3.90, -2.87, -0.58,  0.71,  1.80,           &
+         &        3.08,  4.63,  5.86,  7.21,  8.58, 10.50, 12.10,    &
+         &       12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39,    &
+         &       19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68,    &
+         &       22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63,    &
+         &       23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80,    &
+         &       24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89,    &
+         &       27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96,    &
+         &       31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39,    &
+         &       33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87,    &
+         &       38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00,    &
+         &       45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81,    &
+         &       52.57 /
+      !!----------------------------------------------------------------------
+
+      cycle = 360.0
+      ilc = 0
+      icent = year / 100
+      yr = year - icent * 100
+      t = icent - 20
+!
+!     for the following equations
+!     time origin is fixed at 00 hr of jan 1st,2000.
+!     see notes by cartwright
+!
+!!gm  old coding style, use CASE instead  and avoid GOTO (obsolescence in fortran 90)
+!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
+      it = icent - 20
+      if (it) 1,2,2
+    1 iday = it/4 -it
+      go to 3
+    2 iday = (it+3)/4 - it
+!
+!     t is in julian century
+!     correction in gegorian calander where only century year divisible
+!     by 4 is leap year.
+!
+    3 continue
+!
+      td = 0.0
+!
+!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
+      if (yr) 4,5,4
+!
+    4 iyd = 365*yr
+      ild = (yr-1)/4
+      if((icent - (icent/4)*4) .eq. 0) ilc = 1
+      td = iyd + ild + ilc
+!
+    5 td = td + iday + vd -1.0 - 0.5
+      t = t + (td/36525.0)
+!
+      ipos=year-1899
+      if (ipos .lt. 0) go to 7
+      if (ipos .gt. 83) go to 6
+!
+      delta = (delt(ipos+1)+delt(ipos))/2.0
+      go to 7
+!
+    6 delta= (65.0-50.5)/20.0*(year-1980)+50.5
+!
+    7 deltat = delta * 1.0e-6
+!
+!!gm   precision of the computation   :  example for s it should be replace by:
+!!gm  s   = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat   ==>  more precise  modify the last digits results
+      s   = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat
+      h   = 280.4661 + 36000.7698 *t + 0.0003*t*t +  11.0*deltat
+      p   =  83.3535 + 4069.0139  *t - 0.0103*t*t +       deltat
+      en  = 234.9555 + 1934.1363  *t - 0.0021*t*t +       deltat
+      p1  = 282.9384 + 1.7195     *t + 0.0005*t*t
+      !
+      nn = s / cycle
+      s = s - nn * cycle
+      IF( s < 0.e0 )   s = s + cycle
+      !
+      nn = h / cycle
+      h = h - cycle * nn
+      IF( h < 0.e0 )   h = h + cycle
+      !
+      nn = p / cycle
+      p = p - cycle * nn
+      IF( p < 0.e0)   p = p + cycle
+      !
+      nn = en / cycle
+      en = en - cycle * nn
+      IF( en < 0.e0 )   en = en + cycle
+      en = cycle - en
+      !
+      nn = p1 / cycle
+      p1 = p1 - nn * cycle
+      !
+   END SUBROUTINE shpen
+
+
+   SUBROUTINE ufset( p, cn, b, a )
+      !!----------------------------------------------------------------------
+      !!                    ***  SUBROUTINE ufset  ***
+      !!                    
+      !! ** Purpose : - calculate nodal parameters for the tides
+      !!                
+      !!----------------------------------------------------------------------
+!!gm  doctor naming of dummy argument variables!!!   and all local variables
+!!gm  nc is jptides_max ????
+      integer nc
+      parameter (nc=15)
+      REAL(wp) p,cn
+      !!
+      
+!!gm  rad is already a public variable defined in phycst.F90 ....   ==> doctor norme  local real start with "z"
+      REAL(wp) ::   w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad
+      REAL(wp) ::   a(nc), b(nc)
+      INTEGER  ::   ndc, k
+      !!----------------------------------------------------------------------
+
+      ndc = nc
+
+!    a=f       ,  b =u
+!    t is  zero as compared to tifa.
+!! use rad defined in phycst   (i.e.  add a USE phycst at the begining of the module
+      rad = 6.2831852d0/360.0
+      pw = p  * rad
+      nw = cn * rad
+      w1 = cos(   nw )
+      w2 = cos( 2*nw )
+      w3 = cos( 3*nw )
+      w4 = sin(   nw )
+      w5 = sin( 2*nw )
+      w6 = sin( 3*nw )
+      w7 = 1. - 0.2505 * COS( 2*pw      ) - 0.1102 * COS(2*pw-nw )   &
+         &    - 0.156  * COS( 2*pw-2*nw ) - 0.037  * COS(     nw )
+      w8 =    - 0.2505 * SIN( 2*pw      ) - 0.1102 * SIN(2*pw-nw )   &
+         &    - 0.0156 * SIN( 2*pw-2*nw ) - 0.037  * SIN(     nw )
+      !
+      a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3
+      b(1) =          0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6
+      !   q1
+      a(2) = a(1)
+      b(2) = b(1)
+      !   o1
+      a(3) = 1.0
+      b(3) = 0.0
+      !   p1
+      a(4) = 1.0
+      b(4) = 0.0
+      !   s1
+      a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3
+      b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6
+      !   k1
+      a(6) =1.0004 -0.0373*w1+ 0.0002*w2
+      b(6) = -0.0374*w4
+      !  2n2
+      a(7) = a(6)
+      b(7) = b(6)
+      !  mu2
+      a(8) = a(6)
+      b(8) = b(6)
+      !   n2
+      a(9) = a(6)
+      b(9) = b(6)
+      !  nu2
+      a(10) = a(6)
+      b(10) = b(6)
+      !   m2
+      a(11) = SQRT( w7 * w7 + w8 * w8 )
+      b(11) = ATAN( w8 / w7 )
+!!gmuse rpi  instead of 3.141992 ???   true pi is rpi=3.141592653589793_wp  .....   ????
+      IF( w7 < 0.e0 )   b(11) = b(11) + 3.141992
+      !   l2
+      a(12) = 1.0
+      b(12) = 0.0
+      !   t2
+      a(13)= a(12)
+      b(13)= b(12)
+      !   s2
+      a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3
+      b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6
+      !   k2
+      a(15) = a(6)*a(6)
+      b(15) = 2*b(6)
+      !   m4
+!!gm  old coding,  remove GOTO and label of lines
+!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
+      DO 40 k = 1,ndc
+         b(k) = b(k)/rad
+32       if (b(k)) 34,35,35
+34       b(k) = b(k) + 360.0
+         go to 32
+35       if (b(k)-360.0) 40,37,37
+37       b(k) = b(k)-360.0
+         go to 35
+40    continue
+      !
+   END SUBROUTINE ufset
+   
+
+   SUBROUTINE vset( s,h,p,en,p1,v)
+      !!----------------------------------------------------------------------
+      !!                    ***  SUBROUTINE vset  ***
+      !!                   
+      !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run
+      !!                
+      !!----------------------------------------------------------------------
+!!gm  doctor naming of dummy argument variables!!!   and all local variables
+!!gm  nc is jptides_max ????
+!!gm  en argument is not used: suppress it ?
+      integer nc
+      parameter (nc=15)
+      real(wp) s,h,p,en,p1
+      real(wp)   v(nc)
+      !!
+      integer ndc, k
+      !!----------------------------------------------------------------------
+
+      ndc = nc
+      !   v s  are computed here.
+      v(1) =-3*s +h +p +270      ! Q1
+      v(2) =-2*s +h +270.0       ! O1
+      v(3) =-h +270              ! P1
+      v(4) =180                  ! S1
+      v(5) =h +90.0              ! K1
+      v(6) =-4*s +2*h +2*p       ! 2N2
+      v(7) =-4*(s-h)             ! MU2
+      v(8) =-3*s +2*h +p         ! N2
+      v(9) =-3*s +4*h -p         ! MU2
+      v(10) =-2*s +2*h           ! M2
+      v(11) =-s +2*h -p +180     ! L2 
+      v(12) =-h +p1              ! T2
+      v(13) =0                   ! S2
+      v(14) =h+h                 ! K2
+      v(15) =2*v(10)             ! M4
+!
+!!gm  old coding,  remove GOTO and label of lines
+!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
+      do 72 k = 1, ndc
+69       if( v(k) )   70,71,71
+70       v(k) = v(k)+360.0
+         go to 69
+71       if( v(k) - 360.0 )   72,73,73
+73       v(k) = v(k)-360.0
+         go to 71
+72    continue
+      !
+   END SUBROUTINE vset
+
+#else
    !!----------------------------------------------------------------------
    !!   Dummy module         NO Unstruct Open Boundary Conditions for tides
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 2813)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 2814)
@@ -24,4 +24,6 @@
    IMPLICIT NONE
    PRIVATE   
+ 
+   PUBLIC   fld_map    ! routine called by tides_init
 
    TYPE, PUBLIC ::   FLD_N      !: Namelist field informations
@@ -629,23 +631,28 @@
       !!                using a general mapping (for open boundaries)
       !!----------------------------------------------------------------------
+#if defined key_obc
       USE obc_oce, ONLY:  dta_global         ! workspace to read in global data arrays
+#endif 
 
       INTEGER                   , INTENT(in ) ::   num     ! stream number
       CHARACTER(LEN=*)          , INTENT(in ) ::   clvar   ! variable name
-      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta     ! output field on model grid
+      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta   ! output field on model grid (2 dimensional)
       INTEGER                   , INTENT(in ) ::   nrec    ! record number to read (ie time slice)
       INTEGER,  DIMENSION(:)    , INTENT(in ) ::   map     ! global-to-local mapping indices
       !!
-      INTEGER                  ::   ipi      ! length of boundary data on local process
-      INTEGER                  ::   ipj      ! length of dummy dimension ( = 1 )
-      INTEGER                  ::   ipk      ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
-      INTEGER                  ::   ilendta  ! length of data in file
-      INTEGER                  ::   idvar    ! variable ID
-      INTEGER                  ::   ib, ik   ! loop counters
-      INTEGER                  ::   ierr
-      !!
-      CHARACTER(len=80)                   :: zfile
+      INTEGER                                 ::   ipi      ! length of boundary data on local process
+      INTEGER                                 ::   ipj      ! length of dummy dimension ( = 1 )
+      INTEGER                                 ::   ipk      ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
+      INTEGER                                 ::   ilendta  ! length of data in file
+      INTEGER                                 ::   idvar    ! variable ID
+      INTEGER                                 ::   ib, ik   ! loop counters
+      INTEGER                                 ::   ierr
+      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read ! work space for global data
       !!---------------------------------------------------------------------
             
+#if defined key_obc
+      dta_read => dta_global
+#endif
+
       ipi = SIZE( dta, 1 )
       ipj = 1
@@ -655,14 +662,21 @@
       ilendta = iom_file(num)%dimsz(1,idvar)
       IF(lwp) WRITE(numout,*) 'Dim size for ',TRIM(clvar),' is ', ilendta
-
-      CALL iom_get ( num, jpdom_unknown, clvar, dta_global(1:ilendta,1:ipj,1:ipk), nrec )
+      IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk
+
+      SELECT CASE( ipk )
+      CASE(1)   
+         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec )
+      CASE DEFAULT
+         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )
+      END SELECT
       !
       DO ib = 1, ipi
          DO ik = 1, ipk
-            dta(ib,1,ik) =  dta_global(map(ib),1,ik)
+            dta(ib,1,ik) =  dta_read(map(ib),1,ik)
          END DO
       END DO
 
    END SUBROUTINE fld_map
+
 
    SUBROUTINE fld_rot( kt, sd )
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 2813)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 2814)
@@ -47,4 +47,5 @@
    USE obcini          ! open boundary cond. initialization (obc_init routine)
    USE obcdta          ! open boundary cond. initialization (obc_dta_init routine)
+   USE obctides        ! open boundary cond. initialization (tide_init routine)
    USE istate          ! initial state setting          (istate_init routine)
    USE ldfdyn          ! lateral viscosity setting      (ldfdyn_init routine)
@@ -296,4 +297,5 @@
       IF( lk_obc        )   CALL     obc_init       ! Open boundaries initialisation
       IF( lk_obc        )   CALL     obc_dta_init   ! Open boundaries initialisation of external data arrays
+      IF( lk_obc        )   CALL     tide_init      ! Open boundaries initialisation of tidal harmonic forcing
 
                             CALL  istate_init   ! ocean initial state (Dynamics and tracers)
