Index: NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/BDY/bdy_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/BDY/bdy_oce.F90	(revision 15736)
+++ NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/BDY/bdy_oce.F90	(revision 15737)
@@ -43,4 +43,7 @@
       INTEGER          , DIMENSION(2)   ::  nread
       LOGICAL                           ::  lneed_ssh
+      !--- davbyr
+      LOGICAL                           ::  lforced_ssh 
+	  !--- END davbyr
       LOGICAL                           ::  lneed_dyn2d
       LOGICAL                           ::  lneed_dyn3d
@@ -117,4 +120,9 @@
    REAL(wp), DIMENSION(jp_bdy) ::   rice_hpnd               !: pond thick. of incoming sea ice
    REAL(wp), DIMENSION(jp_bdy) ::   rice_hlid               !: pond lid thick. of incoming sea ice
+   
+   !  davbyr 
+   LOGICAL, DIMENSION(jp_bdy) ::   ln_ssh_bdy               !: =T USE SSH BDY - name list switch
+   REAL(wp), DIMENSION(jp_bdy) ::  rn_ssh_shift             !: =F SHIFT SSH AT A BORDER BY rn_ssh_shift m_
+   !  END davbyr
    !
    !!----------------------------------------------------------------------
Index: NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/BDY/bdydta.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/BDY/bdydta.F90	(revision 15736)
+++ NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/BDY/bdydta.F90	(revision 15737)
@@ -362,4 +362,19 @@
          ENDIF
       ENDIF
+      
+      ! davbyr - add a shift to the boundary + free elevation Enda, JT from NEMO RAN 3.6
+      DO jbdy = 1, nb_bdy
+         IF( dta_bdy(jbdy)%lneed_ssh ) THEN
+            igrd  = 1
+            DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd)   ! ssh is used only on the rim
+                ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                dta_bdy(jbdy)%ssh(ib) = dta_bdy(jbdy)%ssh(ib) + rn_ssh_shift(jbdy) * tmask(ii,ij,1)
+                IF( .NOT. dta_bdy(jbdy)%lforced_ssh ) dta_bdy(jbdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)
+             END DO
+         END IF
+      END DO
+      !--- END davbyr
+      
       !
       IF( ln_timing )   CALL timing_stop('bdy_dta')
Index: NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/BDY/bdyini.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/BDY/bdyini.F90	(revision 15736)
+++ NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/BDY/bdyini.F90	(revision 15737)
@@ -69,4 +69,8 @@
          &             ln_vol, nn_volctl, nn_rimwidth
          !
+!     davbyr Propagating ENDA's stuff from 3.6
+      NAMELIST/nambdy_ssh/ ln_ssh_bdy, rn_ssh_shift
+      INTEGER  ::   ib_bdy
+!     END davbyr
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!----------------------------------------------------------------------
@@ -97,4 +101,29 @@
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist' )
       IF(lwm) WRITE ( numond, nambdy )
+      
+      ! davbyr Propagating ENDA's stuff from 3.6
+      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries
+      READ  ( numnam_ref, nambdy_ssh, IOSTAT = ios, ERR = 905)
+905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_ssh in reference namelist' )
+
+      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries
+      READ  ( numnam_cfg, nambdy_ssh, IOSTAT = ios, ERR = 906)
+906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_ssh in configuration namelist' )
+      IF(lwm) WRITE ( numond, nambdy_ssh )
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'nambdy_ssh : use of ssh boundaries'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~'
+      IF(lwp) WRITE(numout,*) '      ln_ssh_bdy: '
+      DO ib_bdy = 1,nb_bdy
+        IF(lwp) WRITE(numout,*) '      ln_ssh_bdy  (',ib_bdy,'): ',ln_ssh_bdy(ib_bdy)
+      IF(lwp) WRITE(numout,*) '      rn_ssh_shift: '
+      ENDDO
+      DO ib_bdy = 1,nb_bdy
+        IF(lwp) WRITE(numout,*) '      rn_ssh_shift(',ib_bdy,'): ',rn_ssh_shift(ib_bdy)
+      ENDDO
+      IF(lwp) WRITE(numout,*) '~~~~~~~~'
+      IF(lwp) WRITE(numout,*)
+!     END davbyr
 
       IF( .NOT. Agrif_Root() ) ln_bdy = .FALSE.   ! forced for Agrif children
@@ -201,4 +230,14 @@
          dta_bdy(ib_bdy)%lneed_ssh   = cn_dyn2d(ib_bdy) == 'flather'
          dta_bdy(ib_bdy)%lneed_dyn2d = cn_dyn2d(ib_bdy) /= 'none'
+
+		 ! davbyr propagating JT override dta_bdy(ib_bdy)%ll_ssh with namelist value (ln_ssh_bdy)
+         dta_bdy(ib_bdy)%lforced_ssh = ln_ssh_bdy(ib_bdy)
+         IF(lwp) WRITE(numout,*) 'nambdy_ssh : use of ssh boundaries'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~'
+         IF(lwp) WRITE(numout,*) '      ib_bdy: ',ib_bdy
+         IF(lwp) WRITE(numout,*) '      dta_bdy(ib_bdy)%lneed_ssh  : ',dta_bdy(ib_bdy)%lneed_ssh
+         IF(lwp) WRITE(numout,*) '      dta_bdy(ib_bdy)%lforced_ssh: ',dta_bdy(ib_bdy)%lforced_ssh
+         IF(lwp) WRITE(numout,*) '~~~~~~~~'
+         ! END davbyr
 
          IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
@@ -598,5 +637,7 @@
             DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
                ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same weights
-               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 )      ! tanh formulation
+!               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 )      ! tanh formulation
+               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( ir - 1 ) * 0.5 &
+                                                & *(10./FLOAT(nn_rimwidth(ib_bdy))) ) ! JGraham:modified for rim=15
                !               idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.  ! quadratic
                !               idx_bdy(ib_bdy)%nbw(ib,igrd) =  REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy))       ! linear
Index: NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/BDY/bdytides.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/BDY/bdytides.F90	(revision 15736)
+++ NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/BDY/bdytides.F90	(revision 15737)
@@ -147,5 +147,4 @@
                !
                ! SSH fields
-               IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
                   clfile = TRIM(filtide)//'_grid_T.nc'
                   CALL iom_open( clfile , inum ) 
@@ -154,4 +153,5 @@
                      CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
                      CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 
+                     IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
                      DO ib = 1, SIZE(dta%ssh)
                         ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -160,10 +160,9 @@
                         td%ssh0(ib,itide,2) = zti(ii,ij)
                      END DO
+                     ENDIF
                   END DO
                   CALL iom_close( inum )
-               END IF
                !
                ! U fields
-               IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
                   clfile = TRIM(filtide)//'_grid_U.nc'
                   CALL iom_open( clfile , inum ) 
@@ -172,4 +171,5 @@
                      CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
                      CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
+                     IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
                      DO ib = 1, SIZE(dta%u2d)
                         ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -178,10 +178,9 @@
                         td%u0(ib,itide,2) = zti(ii,ij)
                      END DO
+                  END IF
                   END DO
-                  CALL iom_close( inum )
-               END IF
+               CALL iom_close( inum )
                !
                ! V fields
-               IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
                   clfile = TRIM(filtide)//'_grid_V.nc'
                   CALL iom_open( clfile , inum ) 
@@ -190,4 +189,5 @@
                      CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
                      CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
+                     IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
                      DO ib = 1, SIZE(dta%v2d)
                         ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -196,7 +196,7 @@
                         td%v0(ib,itide,2) = zti(ii,ij)
                      END DO
+                  END IF
                   END DO
-                  CALL iom_close( inum )
-               END IF
+               CALL iom_close( inum )
                !
                DEALLOCATE( ztr, zti ) 
Index: NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/DOM/dtatsd.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/DOM/dtatsd.F90	(revision 15736)
+++ NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/DOM/dtatsd.F90	(revision 15737)
@@ -22,5 +22,6 @@
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
-
+   USE iom
+   
    IMPLICIT NONE
    PRIVATE
@@ -31,11 +32,18 @@
    !                                  !!* namtsd  namelist : Temperature & Salinity Data *
    LOGICAL , PUBLIC ::   ln_tsd_init   !: T & S data flag
+   LOGICAL , PUBLIC ::   ln_tsd_interp !: vertical interpolation flag
    LOGICAL , PUBLIC ::   ln_tsd_dmp    !: internal damping toward input data flag
+   INTEGER , PUBLIC ::   nn_tsd_zone    !: 1=AMM7 July on-shelf, 2=AMM7 July off-shelf, 3=AMM7 Jan on-shelf, 4=AMM7 Jan off-shelf
+   REAL(wp), PUBLIC ::   rn_sal_sf, rn_mld_sf, rn_maxdep_sf, rn_c0_sf, rn_c1_sf
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read)
+   INTEGER                                 ::   jpk_init , inum_dta
+   INTEGER                                 ::   id ,linum   ! local integers
+   INTEGER                                 ::   zdim(4)
+
 
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
+   !! $Id: dtatsd.F90 10213 2018-10-23 14:40:09Z aumont $ 
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -53,15 +61,17 @@
       LOGICAL, INTENT(in), OPTIONAL ::   ld_tradmp   ! force the initialization when tradp is used
       !
-      INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3   ! local integers
+      INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3, ierr4, ierr5   ! local integers
       !!
       CHARACTER(len=100)            ::   cn_dir          ! Root directory for location of ssr files
-      TYPE(FLD_N), DIMENSION( jpts) ::   slf_i           ! array of namelist informations on the fields to read
-      TYPE(FLD_N)                   ::   sn_tem, sn_sal
+      TYPE(FLD_N), DIMENSION(jpts+2)::   slf_i           ! array of namelist informations on the fields to read
+      TYPE(FLD_N)                   ::   sn_tem, sn_sal, sn_dep, sn_msk
+      
       !!
-      NAMELIST/namtsd/   ln_tsd_init, ln_tsd_dmp, cn_dir, sn_tem, sn_sal
+      NAMELIST/namtsd/   ln_tsd_init, ln_tsd_interp, ln_tsd_dmp, cn_dir, sn_tem, sn_sal, sn_dep, sn_msk, nn_tsd_zone,   &
+                         rn_sal_sf, rn_mld_sf, rn_maxdep_sf, rn_c0_sf, rn_c1_sf 
       !!----------------------------------------------------------------------
       !
       !  Initialisation
-      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0
+      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0  ; ierr4 = 0  ;  ierr5 = 0 
       !
       REWIND( numnam_ref )              ! Namelist namtsd in reference namelist : 
@@ -80,6 +90,7 @@
          WRITE(numout,*) '~~~~~~~~~~~~ '
          WRITE(numout,*) '   Namelist namtsd'
-         WRITE(numout,*) '      Initialisation of ocean T & S with T &S input data   ln_tsd_init = ', ln_tsd_init
-         WRITE(numout,*) '      damping of ocean T & S toward T &S input data        ln_tsd_dmp  = ', ln_tsd_dmp
+         WRITE(numout,*) '      Initialisation of ocean T & S with T &S input data   ln_tsd_init   = ', ln_tsd_init
+         WRITE(numout,*) '      Interpolation of initial conditions in the vertical  ln_tsd_interp = ', ln_tsd_interp
+         WRITE(numout,*) '      damping of ocean T & S toward T &S input data        ln_tsd_dmp    = ', ln_tsd_dmp
          WRITE(numout,*)
          IF( .NOT.ln_tsd_init .AND. .NOT.ln_tsd_dmp ) THEN
@@ -94,19 +105,49 @@
          ln_tsd_init = .FALSE.
       ENDIF
+      IF( ln_tsd_interp .AND. ln_tsd_dmp ) THEN
+            CALL ctl_stop( 'dta_tsd_init: Tracer damping and vertical interpolation not yet configured' )   ;   RETURN
+      ENDIF
+      IF( ln_tsd_interp .AND. LEN(TRIM(sn_msk%wname)) > 0 ) THEN
+            CALL ctl_stop( 'dta_tsd_init: Using vertical interpolation and weights files not recommended' )   ;   RETURN
+      ENDIF
       !
       !                             ! allocate the arrays (if necessary)
       IF( ln_tsd_init .OR. ln_tsd_dmp ) THEN
          !
-         ALLOCATE( sf_tsd(jpts), STAT=ierr0 )
+         IF( ln_tsd_interp ) THEN
+           ALLOCATE( sf_tsd(jpts+2), STAT=ierr0 ) ! to carry the addtional depth information
+         ELSE
+           ALLOCATE( sf_tsd(jpts  ), STAT=ierr0 ) 
+         ENDIF 
          IF( ierr0 > 0 ) THEN
             CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' )   ;   RETURN
          ENDIF
          !
-                                ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 )
-         IF( sn_tem%ln_tint )   ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )
-                                ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 )
-         IF( sn_sal%ln_tint )   ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
-         !
-         IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN
+         slf_i(jp_tem) = sn_tem   ;   slf_i(jp_sal) = sn_sal
+         IF( ln_tsd_interp ) slf_i(jp_dep) = sn_dep   ;   slf_i(jp_msk) = sn_msk
+         CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd' )
+
+         IF( ln_tsd_interp ) THEN
+            CALL fld_clopn ( sf_tsd(jp_dep) ) 
+            IF(lwp) WRITE(numout,*) 'INFO: ', sf_tsd(jp_dep)%num, sn_dep%clvar
+            id = iom_varid( sf_tsd(jp_dep)%num, sn_dep%clvar, zdim )
+            jpk_init = zdim(3)
+            IF(lwp) WRITE(numout,*) 'Dimension of veritcal coordinate in ICs: ', jpk_init
+            !
+                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr0 )
+            IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr1 )
+                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr2 )
+            IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr3 )  
+                                 ALLOCATE( sf_tsd(jp_dep)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr4 )
+                                 ALLOCATE( sf_tsd(jp_msk)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr5 )
+         ELSE
+                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 )
+            IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )
+                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 )
+            IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )  
+         ENDIF ! ln_tsd_interp
+
+         !
+         IF( ierr0 + ierr1 + ierr2 + ierr3 + ierr4 + ierr5 > 0 ) THEN
             CALL ctl_stop( 'dta_tsd : unable to allocate T & S data arrays' )   ;   RETURN
          ENDIF
@@ -138,8 +179,7 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data
       !
-      INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies
+      INTEGER ::   ji, jj, jk, jl, jk_init   ! dummy loop indicies
       INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1   ! local integers
       REAL(wp)::   zl, zi                             ! local scalars
-      REAL(wp), DIMENSION(jpk) ::  ztp, zsp   ! 1D workspace
       !!----------------------------------------------------------------------
       !
@@ -176,47 +216,48 @@
 !!gm end
       !
-      ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)    ! NO mask
-      ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) 
-      !
-      IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==!
+      IF( ln_tsd_interp ) THEN
          !
          IF( kt == nit000 .AND. lwp )THEN
             WRITE(numout,*)
-            WRITE(numout,*) 'dta_tsd: interpolates T & S data onto the s- or mixed s-z-coordinate mesh'
-         ENDIF
-         !
-         DO jj = 1, jpj                         ! vertical interpolation of T & S
-            DO ji = 1, jpi
-               DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
+            WRITE(numout,*) 'dta_tsd: interpolates T & S data onto current mesh'
+         ENDIF
+         DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
+            DO jj= 1, jpj
+               DO ji= 1, jpi
                   zl = gdept_0(ji,jj,jk)
-                  IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data
-                     ztp(jk) =  ptsd(ji,jj,1    ,jp_tem)
-                     zsp(jk) =  ptsd(ji,jj,1    ,jp_sal)
-                  ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data
-                     ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem)
-                     zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal)
-                  ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
-                     DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
-                        IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
-                           zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
-                           ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi 
-                           zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi
+                  IF( zl < sf_tsd(jp_dep)%fnow(ji,jj,1) ) THEN                     ! above the first level of data
+                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,1) 
+                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,1)
+                  ELSEIF( zl > sf_tsd(jp_dep)%fnow(ji,jj,jpk_init) ) THEN          ! below the last level of data
+                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jpk_init)
+                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jpk_init)
+                  ELSE                                                             ! inbetween : vertical interpolation between jk_init & jk_init+1
+                     DO jk_init = 1, jpk_init-1                                    ! when  gdept(jk_init) < zl < gdept(jk_init+1)
+                        IF( sf_tsd(jp_msk)%fnow(ji,jj,jk_init+1) == 0 ) THEN       ! if there is no data fill down
+                           sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init)
+                           sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init)
+                        ENDIF
+                        IF( (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init)) * (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)) <= 0._wp ) THEN
+                           zi = ( zl - sf_tsd(jp_dep)%fnow(ji,jj,jk_init) ) / &
+                        &       (sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_dep)%fnow(ji,jj,jk_init))
+                           ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init) + &
+                        &                          (sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_tem)%fnow(ji,jj,jk_init)) * zi
+                           ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init) + &
+                        &                          (sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_sal)%fnow(ji,jj,jk_init)) * zi
                         ENDIF
                      END DO
                   ENDIF
-               END DO
-               DO jk = 1, jpkm1
-                  ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
-                  ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk)
-               END DO
-               ptsd(ji,jj,jpk,jp_tem) = 0._wp
-               ptsd(ji,jj,jpk,jp_sal) = 0._wp
-            END DO
+               ENDDO
+            ENDDO
          END DO
-         ! 
-      ELSE                                !==   z- or zps- coordinate   ==!
-         !                             
-         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) * tmask(:,:,:)    ! Mask
-         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) * tmask(:,:,:)
+         !
+         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) *tmask(:,:,:)
+         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) *tmask(:,:,:)
+      ELSE                                
+         !
+         CALL ctl_warn('dta_tsd: T & S data are assumed to be on the current mesh. No interpolation performed')
+         !                  
+         ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)  * tmask(:,:,:)  ! Mask
+         ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:)  * tmask(:,:,:)
          !
          IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
@@ -248,4 +289,6 @@
                                         DEALLOCATE( sf_tsd(jp_sal)%fnow )     ! S arrays in the structure
          IF( sf_tsd(jp_sal)%ln_tint )   DEALLOCATE( sf_tsd(jp_sal)%fdta )
+         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_dep)%fnow )     ! T arrays in the structure
+         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_msk)%fnow )     ! T arrays in the structure
                                         DEALLOCATE( sf_tsd              )     ! the structure itself
       ENDIF
Index: NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/SBC/sbctide.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/SBC/sbctide.F90	(revision 15736)
+++ NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/SBC/sbctide.F90	(revision 15737)
@@ -16,4 +16,6 @@
    USE ioipsl         ! NetCDF IPSL library
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+
+   USE bdytides ! davbyr - Access to love number
 
    IMPLICIT NONE
@@ -108,5 +110,7 @@
 
       DO jk = 1, nb_harmo
-         zcons = 0.7_wp * Wave(ntide(jk))%equitide * ftide(jk)
+         ! davbyr - Insert variable Love number where once was 0.7
+         zcons = dn_love_number * Wave(ntide(jk))%equitide * ftide(jk)
+         ! END davbyr
          DO ji = 1, jpi
             DO jj = 1, jpj
@@ -119,4 +123,7 @@
                IF    ( Wave(ntide(jk))%nutide == 1 )  THEN  ;  zcs = zcons * SIN( 2._wp*zlat )
                ELSEIF( Wave(ntide(jk))%nutide == 2 )  THEN  ;  zcs = zcons * COS( zlat )**2
+               ! davbyr - Include long period tidal forcing
+               ELSEIF( Wave(ntide(jk))%nutide == 0 )  THEN  ;  zcs = zcons * (0.5_wp-1.5_wp*SIN(zlat)**2._wp)
+               ! END - davbyr
                ELSE                                         ;  zcs = 0._wp
                ENDIF
Index: NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/SBC/tide.h90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/SBC/tide.h90	(revision 15736)
+++ NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/SBC/tide.h90	(revision 15737)
@@ -2,28 +2,113 @@
    !! History :  3.2  !  2007  (O. Le Galloudec)  Original code
    !!----------------------------------------------------------------------
-
+   !! TIDES ADDED  ! 2017 (Nico Bruneau) 
+   !! Following this document that seems to match implemented code 
+   !!    https://docs.lib.noaa.gov/rescue/cgs_specpubs/QB275U35no981924.pdf
+   !! see page 189 for some proposed values
+   !!
+   !! The convention which seems to have been chosen is the Shureman one and 
+   !! not the Cartwright and Tayer (1971)
+   !! This is probably due to the fact the Schureman has a solar calendar 
+   !! while Cartwright and Tayer is based on a lunar calendar
+   !!
+   !! Therefore the coefficient are not the Doodson number but the one 
+   !! defined by Schureman. For example :
+   !! M2 : Doodson   : 2  0 0 0 0 0
+   !!      Schureman : 2 -2 2 0 0 0
+   !!
+   !! Components 1-34 are for FES 2014
+   !! Components >= 35 are the one that were initially present in NEMO and not in FES14
+   !!                  keep in mind than equitide coefficient have been ajusted for the
+   !!                  34 FES 2014 constituents
+   !! 
+   !! The different coefficient are as follows
+   !!   - nt   = T  = Number of Julian centuries (36625 days) from Greenwich mean noon on December 31, 1899.
+   !!               = Hour angle of mean sun
+   !!   - ns   = s  = mean longitude of the moon
+   !!   - nh   = h  = mean longitude of the sun
+   !!   - np   = p  = mean longitude of the lunar perigee
+   !!   - np1  = p1 = mean longitude of the solar perigee
+   !!   - shift appears in table as a bias in degree
+   !!   - nksi Coefficient for the longitude in moon's orbit of lunar intersection
+   !!   - nu0 Coefficient for the right ascension of lunar intersection
+   !!   - nu1 Coefficient for the term in argument of lunisolar constituent K1
+   !!   - nu2 Coefficient for the term in argument of lunisolar constituent K2
+   !!   - R = ??? 
+   !!   - Formula = Nodal factor function; see the table of Schureman. Implemented in tide_mod.F90 
+   !! 
+   !! The equitide parameter seems to be the equilibrium tide amplitude corrected
+   !! with the C_n^m coefficient: see Cartwright and Tayer (1971) equation 12 
+   !! and Table 2
+   !! As an example in their Table 4c (p66), M2 (200000) has an amplitude of 
+   !! around 0.63186 m
+   !! Table 2, give us a correction of m = 2, n = 2 (semi-diurnal)
+   !! 0.63186*3*sqrt( 5 / 96 / pi ) = 0.24407 
+   !! very close to the one define originally here : 0.242297
+   !! Third order terms are neglected
+   !!
+   !! So to correct (to match what is implemented in sbctide.F90 - take care CT71 uses co-latitude):
+   !!    - long wave : Amplitude from CT71 * [ -1   * sqrt( 5 /  4 / pi ) ]
+   !!    - diurnal   : Amplitude from CT71 * [ -3/2 * sqrt( 5 / 24 / pi ) ]
+   !!    - semi-diur : Amplitude from CT71 * [  3   * sqrt( 5 / 96 / pi ) ]
+   !!
+   !! ATTENTION: convention seems to be to have a positive coefficient and a 180 shift to 
+   !!            represent negative value. to be confirmed though.
+   !!
+   !! All equtide were computed using the last epocs from Cartwright and Tayer (1971) multiply by
+   !! the corresponding coefficient of their table 2
+   !! 
+   !! nutide is used to compute tide potential - it uses a different formulation depending of nutide
+   !! see sbctide.F90 in function tide_init_potential
+   !!
+   !! Some random note
+   !! in cnes fes tool:
+   !!        Msf has nksi = 2 and nnu0 = -2 which is reverse from Schureman (I kept the Schureman one)
+   !!
+   !!----------------------------------------------------------------------
+   !
    !             !! 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    )
+   !
+   ! Long Period Tides
+   Wave( 1) = tide(  'SA'     , 0.003103 ,    0   ,  0 ,  0 ,  1 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    )
+   Wave( 2) = tide(  'SSA'    , 0.019523 ,    0   ,  0 ,  0 ,  2 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    )
+   Wave( 3) = tide(  'MM'     , 0.022191 ,    0   ,  0 ,  1 ,  0 , -1 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,   73    )
+   Wave( 4) = tide(  'MF'     , 0.042023 ,    0   ,  0 ,  2 ,  0 ,  0 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    )
+   Wave( 5) = tide(  'MTM'    , 0.008042 ,    0   ,  0 ,  3 ,  0 , -1 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    )
+   Wave( 6) = tide(  'MSF'    , 0.003671 ,    0   ,  0 ,  2 , -2 ,  0 ,  0  ,    0  , -2   ,  2   ,  0   ,  0   , 0 ,   78    )
+   Wave( 7) = tide(  'MSQM'   , 0.001293 ,    0   ,  0 ,  4 , -2 ,  0 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    )
+   !
+   ! Diurnal Tides
+   Wave( 8) = tide(  'K1'     , 0.142442 ,    1   ,  1 ,  0 ,  1 ,  0 ,  0  ,  -90  ,  0   ,  0   , -1   ,  0   , 0 ,   227   )
+   Wave( 9) = tide(  'O1'     , 0.101277 ,    1   ,  1 , -2 ,  1 ,  0 ,  0  ,  +90  ,  2   , -1   ,  0   ,  0   , 0 ,    75   )
+   Wave(10) = tide(  'Q1'     , 0.019383 ,    1   ,  1 , -3 ,  1 ,  1 ,  0  ,  +90  ,  2   , -1   ,  0   ,  0   , 0 ,    75   )
+   Wave(11) = tide(  'P1'     , 0.047145 ,    1   ,  1 ,  0 , -1 ,  0 ,  0  ,  +90  ,  0   ,  0   ,  0   ,  0   , 0 ,     0   )
+   Wave(12) = tide(  'S1'     ,-0.001116 ,    1   ,  1 ,  0 ,  0 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,     0   )
+   Wave(13) = tide(  'J1'     ,-0.007961 ,    1   ,  1 ,  1 ,  1 , -1 ,  0  ,  -90  ,  0   , -1   ,  0   ,  0   , 0 ,    76   )
+   !
+   ! Semi-Diurnal Tides
+   Wave(14) = tide(  'M2'     , 0.244083 ,    2   ,  2 , -2 ,  2 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
+   Wave(15) = tide(  'N2'     , 0.046720 ,    2   ,  2 , -3 ,  2 ,  1 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
+   Wave(16) = tide(  'S2'     , 0.113565 ,    2   ,  2 ,  0 ,  0 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,     0   )
+   Wave(17) = tide(  'K2'     , 0.030875 ,    2   ,  2 ,  0 ,  2 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   , -2   , 0 ,   235   )
+   Wave(18) = tide(  'L2'     , 0.006903 ,    2   ,  2 , -1 ,  2 , -1 ,  0  , +180  ,  2   , -2   ,  0   ,  0   , 0 ,   215   )
+   Wave(19) = tide(  'T2'     , 0.006644 ,    2   ,  2 ,  0 , -1 ,  0 ,  1  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,     0   )
+   Wave(20) = tide(  'R2'     , 0.000950 ,    2   ,  2 ,  0 ,  1 ,  0 , -1  , +180  ,  2   ,  0   ,  0   ,  0   , 0 ,     0   )
+   !
+   Wave(21) = tide(  'MU2'    , 0.007451 ,    2   ,  2 , -4 ,  4 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
+   Wave(22) = tide(  'NU2'    , 0.008873 ,    2   ,  2 , -3 ,  4 , -1 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
+   Wave(23) = tide( '2N2'     , 0.006176 ,    2   ,  2 , -4 ,  2 ,  2 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
+   Wave(24) = tide(  'MKS2'   , 0.000000 ,    2   ,  2 , -2 ,  4 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   , -2   , 0 ,     4   )
+   Wave(25) = tide(  'LA2'    , 0.001800 ,    2   ,  2 , -1 ,  0 ,  1 ,  0  , +180  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
+   Wave(26) = tide(  'EPS2'   , 0.001796 ,    2   ,  2 , -5 ,  4 ,  1 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
+   !
+   ! Harmonic and others
+   Wave(27) = tide(  'M3'     , 0.000000 ,    3   ,  3 , -3 ,  3 ,  0 ,  0  ,    0  ,  3   , -3   ,  0   ,  0   , 0 ,   149   )
+   Wave(28) = tide(  'M4'     , 0.000000 ,    4   ,  4 , -4 ,  4 ,  0 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,     1   )
+   Wave(29) = tide(  'M6'     , 0.000000 ,    6   ,  6 , -6 ,  6 ,  0 ,  0  ,    0  ,  6   , -6   ,  0   ,  0   , 0 ,    18   )
+   Wave(30) = tide(  'M8'     , 0.000000 ,    8   ,  8 , -8 ,  8 ,  0 ,  0  ,    0  ,  8   , -8   ,  0   ,  0   , 0 ,    20   )
+   Wave(31) = tide(  'N4'     , 0.000000 ,    4   ,  4 , -6 ,  4 ,  2 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,     1   )
+   Wave(32) = tide(  'S4'     , 0.000000 ,    4   ,  4 ,  0 ,  0 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,     0   )
+   Wave(33) = tide(  'MN4'    , 0.000000 ,    4   ,  4 , -5 ,  4 ,  1 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,     1   )
+   Wave(34) = tide(  'MS4'    , 0.000000 ,    4   ,  4 , -2 ,  2 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   )
+   ! 
Index: NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/SBC/tide_mod.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/SBC/tide_mod.F90	(revision 15736)
+++ NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/SBC/tide_mod.F90	(revision 15737)
@@ -16,5 +16,6 @@
    PUBLIC   tide_init_Wave   ! called by tideini and diaharm modules
 
-   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 19   !: maximum number of harmonic
+   ! davbyr: increase maximum number of harmonics from 19 to 34
+   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 34   !: maximum number of harmonic
 
    TYPE, PUBLIC ::    tide
@@ -331,4 +332,10 @@
          zf = zf * zf1 * zf1
          !
+         
+      !--- davbyr 11/2017
+      CASE( 20 )                 !==  formule 20,  compound waves ( 78 x 78 x 78 x 78 )
+         zf1 = nodal_factort(78)
+         zf  = zf1 * zf1 * zf1 * zf1
+      !--- END davbyr
       CASE( 73 )                 !==  formule 73
          zs = sin(sh_I)
Index: NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/SBC/tideini.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/SBC/tideini.F90	(revision 15736)
+++ NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/SBC/tideini.F90	(revision 15737)
@@ -34,4 +34,7 @@
    REAL(wp), PUBLIC ::   rdttideramp     !:
    REAL(wp), PUBLIC ::   rn_scal_load    !:
+   ! davbyr - read love number from namelist
+   REAL(wp), PUBLIC ::   dn_love_number  !:
+   ! END davbyr
    CHARACTER(lc), PUBLIC ::   cn_tide_load   !: 
 
@@ -54,5 +57,5 @@
       !
       NAMELIST/nam_tide/ln_tide, ln_tide_pot, ln_scal_load, ln_read_load, cn_tide_load, &
-                  &     ln_tide_ramp, rn_scal_load, rdttideramp, clname
+                  &     ln_tide_ramp, rn_scal_load, rdttideramp, dn_love_number, clname
       !!----------------------------------------------------------------------
       !
@@ -80,4 +83,6 @@
             WRITE(numout,*) '         Fraction of SSH used in scal. approx.   rn_scal_load = ', rn_scal_load
             WRITE(numout,*) '         Duration (days) of ramp                 rdttideramp  = ', rdttideramp
+            ! davbyr - Love number (one line)
+            WRITE(numout,*) '         Love Number                             dn_love_number = ', dn_love_number
          ENDIF
       ELSE
Index: NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/ZDF/zdfgls.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/ZDF/zdfgls.F90	(revision 15736)
+++ NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/ZDF/zdfgls.F90	(revision 15737)
@@ -110,5 +110,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
+   !! $Id: zdfgls.F90 13511 2020-09-24 08:55:10Z smasson $
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -188,6 +188,6 @@
          DO jj = 2, jpjm1                      ! bottom friction
             DO ji = fs_2, fs_jpim1   ! vector opt.         
-               zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
-               zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
+               zmsku = 0.5_wp * ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
+               zmskv = 0.5_wp * ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
                ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  &
                   &                                         + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  )
@@ -197,6 +197,6 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zmsku = ( 2._wp - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
-                  zmskv = ( 2._wp - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0)
+                  zmsku = 0.5_wp * ( 2._wp - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
+                  zmskv = 0.5_wp * ( 2._wp - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0)
                   ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  &
                      &                                         + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  )
@@ -413,4 +413,5 @@
                zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1)
                zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp
+               en   (ji,jj,ibot) = z_en
             END DO
          END DO
@@ -429,4 +430,5 @@
                   zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1)
                   zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
+                  en   (ji,jj,itop) = z_en
                END DO
             END DO
@@ -1236,3 +1238,2 @@
    !!======================================================================
 END MODULE zdfgls
-
Index: NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/par_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/par_oce.F90	(revision 15736)
+++ NEMO/branches/UKMO/NEMO_4.0.4_eiwpo12/src/OCE/par_oce.F90	(revision 15737)
@@ -64,5 +64,7 @@
    INTEGER, PUBLIC, PARAMETER ::   jp_tem = 1    !: indice for temperature
    INTEGER, PUBLIC, PARAMETER ::   jp_sal = 2    !: indice for salinity
-
+   !annkat
+   INTEGER, PUBLIC, PARAMETER ::   jp_dep = 3    !: indice for depth
+   INTEGER, PUBLIC, PARAMETER ::   jp_msk = 4    !: indice for mask
    !!----------------------------------------------------------------------
    !!   Domain decomposition
