Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/BDY/bdyini.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/BDY/bdyini.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/BDY/bdyini.F90	(revision 15422)
@@ -167,6 +167,6 @@
       ! Check and write out namelist parameters
       ! -----------------------------------------
-      IF( jperio /= 0 )   CALL ctl_stop( 'bdy_segs: Cyclic or symmetric,',   &
-         &                               ' and general open boundary condition are not compatible' )
+!      IF( jperio /= 0 )   CALL ctl_stop( 'bdy_segs: Cyclic or symmetric,',   &
+!         &                               ' and general open boundary condition are not compatible' )
 
       IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DIA/diaharm.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DIA/diaharm.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DIA/diaharm.F90	(revision 15422)
@@ -5,4 +5,8 @@
    !!======================================================================
    !! History :  3.1  !  2007  (O. Le Galloudec, J. Chanut)  Original code
+   !!  
+   !!   NB: 2017-12 : add 3D harmonic analysis of velocities  
+   !!                 integration of Maria Luneva's development  
+   !!   'key_3Ddiaharm
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
@@ -13,4 +17,8 @@
    USE sbctide         ! Tidal forcing or not
    !
+# if defined key_3Ddiaharm  
+   USE zdf_oce  
+#endif  
+   !
    USE in_out_manager  ! I/O units
    USE iom             ! I/0 library
@@ -33,8 +41,13 @@
    INTEGER         ::   nb_ana        ! Number of harmonics to analyse
 
-   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ana_temp
-   REAL(wp), ALLOCATABLE, DIMENSION(:)       ::   ana_freq, ut, vt, ft
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   out_eta, out_u, out_v
+   INTEGER , ALLOCATABLE, DIMENSION(:)           ::   name 
+   REAL(wp), ALLOCATABLE, DIMENSION(:)           ::   ana_freq, ut, vt, ft  
+# if defined key_3Ddiaharm  
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:)   ::   ana_temp  
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)     ::   out_eta, out_u, out_v, out_w, out_dzi  
+# else  
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)     ::   ana_temp  
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)       ::   out_eta, out_u, out_v  
+# endif
 
    INTEGER  ::   ninco, nsparse
@@ -76,4 +89,7 @@
          WRITE(numout,*)
          WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization'
+# if defined key_3Ddiaharm  
+         WRITE(numout,*) '  - 3D harmonic analysis of currents activated (key_3Ddiaharm)'  
+#endif
          WRITE(numout,*) '~~~~~~~~~~~~~ '
       ENDIF
@@ -155,6 +171,11 @@
          ! Initialize temporary arrays:
          ! ----------------------------
-         ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) )
-         ana_temp(:,:,:,:) = 0._wp
+# if defined key_3Ddiaharm  
+         ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 5, jpk ) )  
+         ana_temp(:,:,:,:,:) = 0._wp  
+# else  
+         ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 3      ) )  
+         ana_temp(:,:,:,:  ) = 0._wp  
+#endif
 
       ENDIF
@@ -175,4 +196,7 @@
       !
       INTEGER  ::   ji, jj, jh, jc, nhc
+# if defined key_3Ddiaharm  
+      INTEGER  :: jk  
+# endif
       REAL(wp) ::   ztime, ztemp
       !!--------------------------------------------------------------------
@@ -190,11 +214,43 @@
                   &    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
                   !
+               ! ssh, ub, vb are stored at the last level of 5d array
                DO jj = 2, jpjm1
                   DO ji = 2, jpim1
-                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp * sshn(ji,jj) * ssmask (ji,jj) ! elevation      
-                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp * un_b(ji,jj) * ssumask(ji,jj) ! u-vel
-                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp * vn_b(ji,jj) * ssvmask(ji,jj) ! v-vel
+                     ! Elevation and currents  
+# if defined key_3Ddiaharm  
+                     ana_temp(ji,jj,nhc,1,jpk) = ana_temp(ji,jj,nhc,1,jpk) + ztemp*sshn(ji,jj)*ssmask (ji,jj)          
+                     ana_temp(ji,jj,nhc,2,jpk) = ana_temp(ji,jj,nhc,2,jpk) + ztemp*un_b(ji,jj)*ssumask(ji,jj)  
+                     ana_temp(ji,jj,nhc,3,jpk) = ana_temp(ji,jj,nhc,3,jpk) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)  
+  
+                     ana_temp(ji,jj,nhc,5,jpk) = ana_temp(ji,jj,nhc,5,jpk)                               &  
+                   &                              + ztemp*bfrva(ji,jj)*vn(ji,jj,mbkv(ji,jj))*ssvmask(ji,jj)  
+                     ana_temp(ji,jj,nhc,4,jpk) = ana_temp(ji,jj,nhc,4,jpk)                               &   
+                   &                              + ztemp*bfrua(ji,jj)*un(ji,jj,mbku(ji,jj))*ssumask(ji,jj)  
+# else  
+                      ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj)          
+                      ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj)  
+                      ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)  
+# endif
                   END DO
                END DO
+               ! 
+# if defined key_3Ddiaharm  
+! 3d velocity and density:  
+             DO jk=1,jpk-1  
+               DO jj = 1,jpj  
+                  DO ji = 1,jpi  
+                     ! density and velocity  
+                     ana_temp(ji,jj,nhc,1,jk) = ana_temp(ji,jj,nhc,1,jk) + ztemp*rhd(ji,jj,jk)  
+                     ana_temp(ji,jj,nhc,2,jk) = ana_temp(ji,jj,nhc,2,jk) + ztemp*(un(ji,jj,jk)-un_b(ji,jj)) &  
+                &                                          *umask(ji,jj,jk)  
+                     ana_temp(ji,jj,nhc,3,jk) = ana_temp(ji,jj,nhc,3,jk) + ztemp*(vn(ji,jj,jk)-vn_b(ji,jj)) &  
+                &                                          *vmask(ji,jj,jk)   
+                     ana_temp(ji,jj,nhc,4,jk) = ana_temp(ji,jj,nhc,4,jk) + ztemp*wn(ji,jj,jk)  
+   
+                     ana_temp(ji,jj,nhc,5,jk) = ana_temp(ji,jj,nhc,5,jk) - 0.5*grav*ztemp*(rhd(ji,jj,jk)+rhd(ji,jj,jk+1))/max(rn2(ji,jj,jk),1.e-8_wp)  
+                  END DO  
+               END DO  
+             ENDDO  
+# endif
             END DO
          END DO
@@ -218,4 +274,7 @@
       !!--------------------------------------------------------------------
       INTEGER  ::   ji, jj, jh, jc, jn, nhan
+# if defined key_3Ddiaharm  
+      INTEGER :: jk  
+# endif 
       INTEGER  ::   ksp, kun, keq
       REAL(wp) ::   ztime, ztime_ini, ztime_end, z1_han
@@ -226,6 +285,4 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
       
-      ALLOCATE( out_eta(jpi,jpj,2*nb_ana), out_u(jpi,jpj,2*nb_ana), out_v(jpi,jpj,2*nb_ana) )
-
       ztime_ini = nit000_han*rdt                 ! Initial time in seconds at the beginning of analysis
       ztime_end = nitend_han*rdt                 ! Final time in seconds at the end of analysis
@@ -233,4 +290,20 @@
       z1_han = 1._wp / REAL(nhan-1) 
       
+# if defined key_3Ddiaharm  
+      ALLOCATE( out_eta(jpi,jpj,jpk,2*nb_ana),   &  
+         &      out_u  (jpi,jpj,jpk,2*nb_ana),   &  
+         &      out_v  (jpi,jpj,jpk,2*nb_ana),   &  
+         &      out_w  (jpi,jpj,jpk,2*nb_ana),   &  
+         &      out_dzi(jpi,jpj,jpk,2*nb_ana) )  
+# else  
+      ALLOCATE( out_eta(jpi,jpj,2*nb_ana),   &  
+         &      out_u  (jpi,jpj,2*nb_ana),   &  
+         &      out_v  (jpi,jpj,2*nb_ana)  )  
+# endif  
+  
+      IF(lwp) WRITE(numout,*) 'ANA F OLD', ft   
+      IF(lwp) WRITE(numout,*) 'ANA U OLD', ut  
+      IF(lwp) WRITE(numout,*) 'ANA V OLD', vt
+
       ninco = 2*nb_ana
 
@@ -260,16 +333,28 @@
       CALL SUR_DETERMINE_INIT
 
-      ! Elevation:
+      ! Density and Elevation:  
+# if defined key_3Ddiaharm  
+    DO jk=1,jpk  
+# endif
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
 
             ! Fill input array
+# if defined key_3Ddiaharm  
+            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,1,jk)  
+# else
             ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,1)
+# endif
             CALL SUR_DETERMINE
             
             ! Fill output array
             DO jh = 1, nb_ana
-               out_eta(ji,jj,jh       ) =  ztmp7((jh-1)*2+1) * ssmask(ji,jj)
-               out_eta(ji,jj,jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj)
+# if defined key_3Ddiaharm  
+               out_eta(ji,jj,jk,jh       ) =  ztmp7((jh-1)*2+1) * ssmask(ji,jj)
+               out_eta(ji,jj,jk,jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj)
+# else  
+               out_eta(ji,jj,   jh       ) =  ztmp7((jh-1)*2+1) * ssmask(ji,jj)
+               out_eta(ji,jj,   jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj)
+# endif
             END DO
          END DO
@@ -281,11 +366,21 @@
 
             ! Fill input array
+# if defined key_3Ddiaharm  
+            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,2,jk)
+# else 
             ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,2)
+# endif
             CALL SUR_DETERMINE
 
             ! Fill output array
             DO jh = 1, nb_ana
-               out_u(ji,jj,       jh) =  ztmp7((jh-1)*2+1) * ssumask(ji,jj)
-               out_u(ji,jj,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj)
+# if defined key_3Ddiaharm  
+               out_u(ji,jj,jk,       jh) =  ztmp7((jh-1)*2+1) * ssumask(ji,jj)
+               out_u(ji,jj,jk,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj)
+# else  
+               out_u(ji,jj,          jh) =  ztmp7((jh-1)*2+1) * ssumask(ji,jj)
+               out_u(ji,jj,   nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj)
+# endif
+
             END DO
 
@@ -298,15 +393,96 @@
 
             ! Fill input array
+# if defined key_3Ddiaharm  
+            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,3,jk)
+# else
             ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,3)
+# endif
             CALL SUR_DETERMINE
 
             ! Fill output array
             DO jh = 1, nb_ana
-               out_v(ji,jj,       jh) =  ztmp7((jh-1)*2+1) * ssvmask(ji,jj)
-               out_v(ji,jj,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj)
+# if defined key_3Ddiaharm  
+               out_v(ji,jj,jk,       jh) =  ztmp7((jh-1)*2+1) * ssvmask(ji,jj)
+               out_v(ji,jj,jk,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj)
+# else  
+               out_v(ji,jj,          jh) =  ztmp7((jh-1)*2+1) * ssvmask(ji,jj)
+               out_v(ji,jj,   nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj)
+# endif
             END DO
 
          END DO
       END DO
+
+# if defined key_3Ddiaharm  
+      ! w- velocity  
+      DO jj = 1, jpj  
+         DO ji = 1, jpi  
+            ! Fill input array  
+            kun=0  
+            DO jh = 1,nb_ana  
+               DO jc = 1,2  
+                  kun = kun + 1  
+                  ztmp4(kun)=ana_temp(ji,jj,kun,4,jk)  
+               END DO  
+            END DO  
+  
+            CALL SUR_DETERMINE(jj+1)  
+  
+            ! Fill output array  
+            DO jh = 1, nb_ana  
+               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)  
+               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)  
+            END DO  
+  
+         END DO  
+      END DO  
+  
+      DO jj = 1, jpj  
+         DO ji = 1, jpi  
+            DO jh = 1, nb_ana  
+               X1=ana_amp(ji,jj,jh,1)  
+               X2=-ana_amp(ji,jj,jh,2)  
+               out_w(ji,jj,jk,       jh)=X1 * tmask_i(ji,jj)  
+               out_w(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj)  
+            END DO  
+         END DO  
+      END DO  
+  
+       ! dzi- isopycnal displacements  
+      DO jj = 1, jpj  
+         DO ji = 1, jpi  
+            ! Fill input array  
+            kun=0  
+            DO jh = 1,nb_ana  
+               DO jc = 1,2  
+                  kun = kun + 1  
+                  ztmp4(kun)=ana_temp(ji,jj,kun,5,jk)  
+               END DO  
+            END DO  
+  
+            CALL SUR_DETERMINE(jj+1)  
+  
+            ! Fill output array  
+            DO jh = 1, nb_ana  
+               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)  
+               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)  
+            END DO  
+  
+         END DO  
+      END DO  
+  
+      DO jj = 1, jpj  
+         DO ji = 1, jpi  
+            DO jh = 1, nb_ana  
+               X1=ana_amp(ji,jj,jh,1)  
+               X2=-ana_amp(ji,jj,jh,2)  
+               out_dzi(ji,jj,jk,       jh)=X1 * tmask_i(ji,jj)  
+               out_dzi(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj)  
+            END DO  
+         END DO  
+      END DO  
+  
+   ENDDO ! jk  
+# endif
       !
       ! clem: we could avoid this call if all the loops were from 1:jpi and 1:jpj
@@ -328,31 +504,163 @@
       !!--------------------------------------------------------------------
       INTEGER  ::   jh
-      !!----------------------------------------------------------------------
+ 
+# if defined key_3Ddiaharm  
+      CHARACTER(LEN=lc) :: cdfile_name_W         ! name of the file created (W-points)  
+      INTEGER  :: jk  
+      REAL(WP), ALLOCATABLE, DIMENSION (:,:,:) :: z3real, z3im   
+      REAL(WP), ALLOCATABLE, DIMENSION (:,:)   :: z2real, z2im        
+# endif  
+!!----------------------------------------------------------------------  
+  
+#if defined key_dimgout  
+      cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc'  
+      cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc'  
+      cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc'  
+#   if defined key_3Ddiaharm  
+      cdfile_name_W = TRIM(cexper)//'_Tidal_harmonics_gridW.dimgproc'  
+#   endif  
+#endif
 
       IF(lwp) WRITE(numout,*) '  '
       IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
+#if defined key_dimgout  
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~  Output files: ', TRIM(cdfile_name_T)  
+      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_U)  
+      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_V)  
+#   if defined key_3Ddiaharm  
+      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_W)  
+#   endif  
+#endif
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
 
-      ! A) Elevation
+# if defined key_3Ddiaharm  
+      ALLOCATE(z3real(jpi,jpj,jpk),z3im(jpi,jpj,jpk),z2real(jpi,jpj),z2im(jpi,jpj))  
+# endif  
+  
+      ! A) density and elevation
       !/////////////
+#if defined key_dimgout  
+      cltext='density amplitude and phase; elevation is level=jpk '  
+      CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2')  
+#else  
+#   if defined key_3Ddiaharm  
+      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp  
+#   endif
       DO jh = 1, nb_ana
+#   if defined key_3Ddiaharm  
+        DO jk=1,jpkm1  
+          z3real(:,:,jk)=out_eta(:,:,jk,jh)  
+          z3im  (:,:,jk)=out_eta(:,:,jk,jh+nb_ana)  
+        ENDDO  
+      z2real(:,:)=out_eta(:,:,jpk,jh); z2im(:,:)=out_eta(:,:,jpk,jh+nb_ana)  
+      CALL iom_put( TRIM(tname(jh))//'x_ro', z3real(:,:,:) )  
+      CALL iom_put( TRIM(tname(jh))//'y_ro', z3im  (:,:,:) )  
+      CALL iom_put( TRIM(tname(jh))//'x'   , z2real(:,:  ) )  
+      CALL iom_put( TRIM(tname(jh))//'y'   , z2im  (:,:  ) )  
+#   else   
+      WRITE(numout,*) "OUTPUT ORI: ", TRIM(tname(jh))//'x', ' & ', TRIM(tname(jh))//'y', MAXVAL(out_eta(:,:,jh))
       CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
       CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,jh+nb_ana) )
-      END DO
-
-      ! B) ubar
+#   endif  
+      END DO  
+#endif  
+  
+      ! B) u
       !/////////
+#if defined key_dimgout  
+      cltext='3d u amplitude and phase; ubar is the last level'  
+      CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2')  
+#else  
+#   if defined key_3Ddiaharm  
+      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp  
+#   endif
       DO jh = 1, nb_ana
-      CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) )
-      CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,jh+nb_ana) )
-      END DO
-
-      ! C) vbar
+#   if defined key_3Ddiaharm  
+        DO jk=1,jpkm1  
+          z3real(:,:,jk)=out_u(:,:,jk,jh)  
+          z3im  (:,:,jk)=out_u(:,:,jk,jh+nb_ana)  
+        ENDDO  
+        z2real(:,:)=out_u(:,:,jpk,jh); z2im(:,:)=out_u(:,:,jpk,jh+nb_ana)  
+        CALL iom_put( TRIM(tname(jh))//'x_u3d', z3real(:,:,:) )  
+        CALL iom_put( TRIM(tname(jh))//'y_u3d', z3im (:,:,:)  )  
+        CALL iom_put( TRIM(tname(jh))//'x_u2d', z2real(:,:) )  
+        CALL iom_put( TRIM(tname(jh))//'y_u2d', z2im (:,:)  )  
+        z2real(:,:)=out_w(:,:,jpk,jh); z2im(:,:)=out_w(:,:,jpk,jh+nb_ana)  
+        CALL iom_put( TRIM(tname(jh))//'x_tabx', z2real(:,:) )  
+        CALL iom_put( TRIM(tname(jh))//'y_tabx', z2im (:,:)  )  
+#   else  
+        CALL iom_put( TRIM(tname(jh))//'x_u2d', out_u(:,:,jh) )  
+        CALL iom_put( TRIM(tname(jh))//'y_u2d', out_u(:,:,nb_ana+jh) )  
+#   endif  
+      END DO  
+#endif  
+  
+      ! C) v
       !/////////
+#if defined key_dimgout  
+      cltext='3d v amplitude and phase; vbar is the last level'  
+      CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2')  
+#else  
+#   if defined key_3Ddiaharm  
+      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp  
+#   endif
       DO jh = 1, nb_ana
-         CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh       ) )
-         CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) )
-      END DO
-      !
+#   if defined key_3Ddiaharm  
+        DO jk=1,jpkm1  
+          z3real(:,:,jk)=out_v(:,:,jk,jh)  
+          z3im  (:,:,jk)=out_v(:,:,jk,jh+nb_ana)  
+        ENDDO  
+        z2real(:,:)=out_v(:,:,jpk,jh); z2im(:,:)=out_v(:,:,jpk,jh+nb_ana)  
+        CALL iom_put( TRIM(tname(jh))//'x_v3d', z3real(:,:,:) )  
+        CALL iom_put( TRIM(tname(jh))//'y_v3d', z3im (:,:,:)  )  
+        CALL iom_put( TRIM(tname(jh))//'x_v2d'  , z2real(:,:) )  
+        CALL iom_put( TRIM(tname(jh))//'y_v2d'  , z2im (:,:)  )  
+        z2real(:,:)=out_dzi(:,:,jpk,jh); z2im(:,:)=out_dzi(:,:,jpk,jh+nb_ana)  
+        CALL iom_put( TRIM(tname(jh))//'x_taby', z2real(:,:) )  
+        CALL iom_put( TRIM(tname(jh))//'y_taby', z2im (:,:)  )  
+#   else  
+         CALL iom_put( TRIM(tname(jh))//'x_v2d', out_v(:,:,jh ) )  
+         CALL iom_put( TRIM(tname(jh))//'y_v2d', out_v(:,:,jh+nb_ana) )  
+#   endif  
+       END DO  
+  
+#endif  
+      ! D) w  
+# if defined key_3Ddiaharm  
+#   if defined key_dimgout  
+      cltext='3d w amplitude and phase; vort_baro is the last level'  
+      CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2')  
+#   else  
+      DO jh = 1, nb_ana  
+        DO jk=1,jpkm1  
+         z3real(:,:,jk)=out_w(:,:,jk,jh)  
+         z3im(:,:,jk)=out_w(:,:,jk,jh+nb_ana)  
+        ENDDO  
+        CALL iom_put( TRIM(tname(jh))//'x_w3d', z3real(:,:,:) )  
+        CALL iom_put( TRIM(tname(jh))//'y_w3d', z3im(:,:,:) )  
+      END DO  
+#   endif  
+  
+!       E) dzi + tau_bot  
+#   if defined key_dimgout  
+      cltext='dzi=g*ro/N2 amplitude and phase'  
+      CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2')  
+#   else  
+      DO jh = 1, nb_ana  
+        DO jk=1,jpkm1  
+         z3real(:,:,jk)=out_dzi(:,:,jk,jh)  
+         z3im(:,:,jk)=out_dzi(:,:,jk,jh+nb_ana)  
+        ENDDO  
+        CALL iom_put( TRIM(tname(jh))//'x_dzi', z3real(:,:,:) )  
+        CALL iom_put( TRIM(tname(jh))//'y_dzi', z3im(:,:,:) )  
+      END DO  
+#   endif  
+# endif   
+  
+      !  
+# if defined key_3Ddiaharm  
+   DEALLOCATE(z3real, z3im, z2real,z2im)  
+# endif
+
    END SUBROUTINE dia_wri_harm
 
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DIA/diaharm_fast.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DIA/diaharm_fast.F90	(revision 15422)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DIA/diaharm_fast.F90	(revision 15422)
@@ -0,0 +1,855 @@
+MODULE diaharm_fast 
+   !!======================================================================
+   !!                       ***  MODULE  example  ***
+   !! Ocean physics:  On line harmonic analyser
+   !!                 
+   !!=====================================================================
+
+#if defined key_diaharm_fast 
+
+   !!----------------------------------------------------------------------
+   !!   'key_harm_ana'  :                Calculate harmonic analysis
+   !!----------------------------------------------------------------------
+   !!   harm_ana        :
+   !!   harm_ana_init   :
+   !!   NB: 2017-12 : add 3D harmonic analysis of velocities
+   !!                 integration of Maria Luneva's development
+   !!   'key_3Ddiaharm'
+   !!----------------------------------------------------------------------
+
+   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 bdy_oce         ! ocean open boundary conditions
+   USE bdytides        ! tidal bdy forcing
+   USE zdfdrg  , ONLY : rCdU_bot   ! bottom friction
+   USE daymod          ! calendar
+   USE tideini
+   USE restart
+   USE ioipsl, ONLY : ju2ymds    ! for calendar
+   !
+   !
+   USE timing          ! preformance summary
+   USE zdf_oce
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! *  Routine accessibility
+   PUBLIC dia_harm_fast                                      ! routine called in step.F90 module
+   LOGICAL, PUBLIC, PARAMETER :: lk_diaharm_fast  = .TRUE.   ! to be run or not
+   LOGICAL, PUBLIC :: lk_diaharm_2D   ! = .TRUE.   ! to run 2d
+   LOGICAL, PUBLIC :: lk_diaharm_3D   ! = .TRUE.   ! to run 3d
+
+   !! * Module variables
+   INTEGER, PARAMETER ::  nharm_max  = jpmax_harmo  ! max number of harmonics to be analysed 
+   INTEGER, PARAMETER ::  nhm_max    = 2*nharm_max+1 
+   INTEGER, PARAMETER ::  nvab       = 2 ! number of 3D variables
+   INTEGER            ::  nharm
+   INTEGER            ::  nhm 
+   INTEGER ::                 & !!! ** toto namelist (namtoto) **
+      nflag  =  1                ! default value of nflag 
+   REAL(wp), DIMENSION(nharm_max) ::                & 
+      om_tide                     ! tidal frequencies ( rads/sec)
+   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:)   ::                & 
+      bzz,c,x    ! work arrays
+   REAL(wp) :: cca,ssa,zm,bt,dd_cumul
+!
+   REAL(wp), PUBLIC ::   fjulday_startharm       !: Julian Day since start of harmonic analysis
+   REAL(wp), PUBLIC, ALLOCATABLE,DIMENSION(:) :: anau, anav, anaf   ! nodel/phase corrections used by diaharmana
+   REAL(WP), ALLOCATABLE,SAVE,DIMENSION(:,:)   :: cc,a
+!
+   INTEGER ::  nvar_2d, nvar_3d    !: number of 2d and 3d variables to analyse
+   INTEGER, ALLOCATABLE,DIMENSION(:) :: m_posi_2d, m_posi_3d
+
+!  Name of variables used in the restart
+   CHARACTER( LEN = 10 ), DIMENSION(5), PARAMETER :: m_varName2d = (/'ssh','u2d','v2d','ubfr','vbfr'/)
+   CHARACTER( LEN = 10 ), DIMENSION(4), PARAMETER :: m_varName3d = (/'rho','u3d','v3d','w3d'/)
+!
+   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,:  ) :: g_cosamp2D, g_sinamp2D, g_cumul_var2D
+   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,:,:) :: g_cosamp3D, g_sinamp3D, g_cumul_var3D  
+!
+   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:)       :: g_out2D,h_out2D  ! arrays for output
+   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:)     :: g_out3D,h_out3D  ! arrays for 3D output
+!
+!  NAMELIST
+   LOGICAL, PUBLIC :: ln_diaharm_store           !: =T  Stores data for harmonic Analysis
+   LOGICAL, PUBLIC :: ln_diaharm_compute         !: =T  Compute harmonic Analysis
+   LOGICAL, PUBLIC :: ln_diaharm_read_restart   !: =T  Read restart from a previous run 
+   LOGICAL, PUBLIC :: ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d
+   INTEGER ::   nb_ana        ! Number of harmonics to analyse
+   CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   tname   ! Names of tidal constituents ('M2', 'K1',...)
+   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   ntide_all ! INDEX within the full set of constituents (tide.h90)
+   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   ntide_sub ! INDEX within the subset of constituents pass in input
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! or LIM 2.0 , UCL-LOCEAN-IPSL (2005)
+   !! or  TOP 1.0 , LOCEAN-IPSL (2005)
+   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/module_example,v 1.3 2005/03/27 18:34:47 opalod Exp $ 
+   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE dia_harm_fast( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE harm_ana  ***
+      !!
+      !! ** Purpose :   Harmonic analyser
+      !!
+      !! ** Method  :   
+      !!
+      !! ** Action  : - first action (share memory array/varible modified
+      !!                in this routine
+      !!              - second action .....
+      !!              - .....
+      !!
+      !! References :
+      !!   Give references if exist otherwise suppress these lines
+      !!
+      !! History :
+      !!   9.0  !  03-08  (Autor Names)  Original code
+      !!        !  02-08  (Author names)  brief description of modifications
+      !!----------------------------------------------------------------------
+      !! * Modules used
+      
+      !! * arguments
+      INTEGER, INTENT( in  ) ::   &  
+         kt                          ! describe it!!!
+
+      !! * local declarations
+      INTEGER  :: ji, jk, jj          ! dummy loop arguments
+      INTEGER  :: jh, i1, i2, jgrid
+      INTEGER  :: j2d, j3d
+      REAL(WP) :: sec2start
+      !!--------------------------------------------------------------------
+
+      IF( ln_timing )   CALL timing_start( 'dia_harm_fast' )
+      IF( kt == nit000   )   CALL harm_ana_init    ! Initialization (first time-step only)
+
+     IF ( ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
+
+      ! this bit done every time step
+      nhm=2*nb_ana+1
+      c(1) = 1.0
+
+      sec2start = nint( (fjulday-fjulday_startharm)*86400._wp ) 
+
+      DO jh=1,nb_ana
+         c(2*jh  ) = anaf(jh)*cos( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
+         c(2*jh+1) = anaf(jh)*sin( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
+      ENDDO 
+
+      ! CUMULATE
+       DO jj = 2, jpjm1
+          DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jh=1,nhm   ! loop harmonic
+
+               DO j2d=1,nvar_2d
+                  IF ( m_posi_2d(j2d) .eq. 1 ) dd_cumul = c(jh) * sshn(ji,jj) * ssmask (ji,jj)             ! analysis elevation
+                  IF ( m_posi_2d(j2d) .eq. 2 ) dd_cumul = c(jh) * un_b(ji,jj) * ssumask(ji,jj)             ! analysis depth average velocities 
+                  IF ( m_posi_2d(j2d) .eq. 3 ) dd_cumul = c(jh) * vn_b(ji,jj) * ssvmask(ji,jj)
+                  IF ( m_posi_2d(j2d) .eq. 4 ) dd_cumul = c(jh) * 0.5*(rCdU_bot(ji+1,jj) + rCdU_bot(ji,jj)) * &
+                                                          un(ji,jj,mbku(ji,jj)) * ssumask(ji,jj) ! analysis bottom friction
+                  IF ( m_posi_2d(j2d) .eq. 5 ) dd_cumul = c(jh) * 0.5*(rCdU_bot(ji,jj+1) + rCdU_bot(ji,jj)) * &
+                                                          vn(ji,jj,mbkv(ji,jj)) * ssvmask(ji,jj)
+                  g_cumul_var2D(jh,ji,jj,j2d) = g_cumul_var2D(jh,ji,jj,j2d) + dd_cumul
+               ENDDO
+
+               DO j3d=1,nvar_3d
+                  DO jk=1,jpkm1
+                     IF ( m_posi_3d(j3d) .eq. 1 ) dd_cumul = c(jh) *  rhd(ji,jj,jk)               * tmask(ji,jj,jk)   
+                     IF ( m_posi_3d(j3d) .eq. 2 ) dd_cumul = c(jh) * ( un(ji,jj,jk)-un_b(ji,jj) ) * umask(ji,jj,jk) 
+                     IF ( m_posi_3d(j3d) .eq. 3 ) dd_cumul = c(jh) * ( vn(ji,jj,jk)-vn_b(ji,jj) ) * vmask(ji,jj,jk)
+                     IF ( m_posi_3d(j3d) .eq. 4 ) dd_cumul = c(jh) *   wn(ji,jj,jk)               * wmask(ji,jj,jk)
+                     g_cumul_var3D(jh,ji,jj,jk,j3d) = g_cumul_var3D(jh,ji,jj,jk,j3d) + dd_cumul
+                  ENDDO
+               ENDDO
+
+            ENDDO     ! end loop harmonic
+         ENDDO        ! end loop lat
+      ENDDO           ! end loop lon
+
+      ! Compute nodal factor cumulative cross-product
+      DO i1=1,nhm
+         DO i2=1,nhm
+            cc(i1,i2)=cc(i1,i2)+c(i1)*c(i2)
+         ENDDO
+      ENDDO
+
+      ! Output RESTART
+      IF( kt == nitrst ) THEN
+         CALL harm_rst_write(kt) ! Dump out data for a restarted run 
+      ENDIF
+
+      ! At End of run
+      IF ( kt ==  nitend ) THEN
+
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'harm_ana : harmonic analysis of tides at end of run'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~'
+
+         IF( ln_diaharm_compute ) THEN
+
+             ! INITIALISE TABLE TO 0
+             IF ( nvar_2d .gt. 0 ) THEN
+                g_cosamp2D = 0.0_wp
+                g_sinamp2D = 0.0_wp
+             ENDIF
+             IF ( nvar_3d .gt. 0 ) THEN
+                g_cosamp3D = 0.0_wp
+                g_sinamp3D = 0.0_wp
+             ENDIF
+
+             ! FIRST OUTPUT 2D VARIABLES
+             DO jgrid=1,nvar_2d    ! loop number of 2d variables (ssh, U2d, V2d, UVfric) to analyse harmonically
+                DO ji=1,jpi        ! loop lon
+                   DO jj=1,jpj     ! loop lat
+                      bt = 1.0_wp; bzz(:) = 0.0_wp
+                      DO jh=1,nhm  ! loop harmonic
+                         bzz(jh) = g_cumul_var2D(jh,ji,jj,jgrid)
+                         bt = bt*bzz(jh)
+                      ENDDO
+                      ! Copy back original cumulated nodal factor
+                      a(:,:) = cc(:,:)
+!                     now do gaussian elimination of the system
+!                     a * x = b
+!                     the matrix x is (a0,a1,b1,a2,b2 ...)
+!                     the matrix a and rhs b solved here for x
+                      x=0.0_wp
+                      IF(bt.ne.0.) THEN
+                        CALL gelim( a, bzz, x, nhm )
+!                       Backup output in variables
+                        DO jh=1,nb_ana
+                           g_cosamp2D(jh,ji,jj,jgrid) = x(jh*2  )
+                           g_sinamp2D(jh,ji,jj,jgrid) = x(jh*2+1)
+                        ENDDO
+                        g_cosamp2D( 0,ji,jj,jgrid) = x(1)
+                        g_sinamp2D( 0,ji,jj,jgrid) = 0.0_wp
+                      ENDIF     ! bt.ne.0.
+                   ENDDO        ! jj
+                ENDDO           ! ji
+             ENDDO              ! jgrid
+
+             ! SECOND OUTPUT 3D VARIABLES
+             DO jgrid=1,nvar_3d     ! loop number of 3d variables rho, U, V, W
+                DO jk=1,jpkm1       ! loop over vertical level
+                   DO ji=1,jpi      ! loop over lon
+                      DO jj=1,jpj   ! loop over lat
+                         bt = 1.0_wp; bzz(:) = 0.0_wp
+                         DO jh=1,nhm
+                            bzz(jh) = g_cumul_var3D(jh,ji,jj,jk,jgrid)
+                            bt = bt*bzz(jh)
+                         ENDDO
+                         ! Copy back original cumulated nodal factor
+                         a(:,:) = cc(:,:)                      
+!                        now do gaussian elimination of the system
+!                        a * x = b
+!                        the matrix x is (a0,a1,b1,a2,b2 ...)
+!                        the matrix a and rhs b solved here for x
+                         x=0.0_wp
+                         IF(bt.ne.0.) THEN
+                           CALL gelim( a, bzz, x, nhm )
+!                          Backup output in variables
+                           DO jh=1,nb_ana
+                              g_cosamp3D(jh,ji,jj,jk,jgrid) = x(jh*2  )
+                              g_sinamp3D(jh,ji,jj,jk,jgrid) = x(jh*2+1)
+                           ENDDO
+                           g_cosamp3D   ( 0,ji,jj,jk,jgrid) = x(1)
+                           g_sinamp3D   ( 0,ji,jj,jk,jgrid) = 0.0_wp
+                        ENDIF     ! bt.ne.0.
+                      ENDDO       ! jj
+                   ENDDO          ! ji
+                ENDDO             ! jk
+             ENDDO                ! jgrid
+
+             CALL harm_ana_out     ! output analysis (last time step)
+
+         ELSE    ! ln_harmana_compute = False 
+             IF(lwp) WRITE(numout,*) " Skipping Computing harmonics at last step"
+
+         ENDIF   ! ln_harmana_compute 
+      ENDIF      ! kt ==  nitend
+
+     ENDIF
+
+      IF( ln_timing )   CALL timing_stop( 'dia_harm_fast' )
+
+   END SUBROUTINE dia_harm_fast 
+
+   SUBROUTINE harm_ana_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE harm_ana_init  ***
+      !!                   
+      !! ** Purpose :   initialization of ....
+      !!
+      !! ** Method  :   blah blah blah ...
+      !!
+      !! ** input   :   Namlist namexa
+      !!
+      !! ** Action  :   ...  
+      !!
+      !! history :
+      !!   9.0  !  03-08  (Autor Names)  Original code
+      !!----------------------------------------------------------------------
+      !! * local declarations
+      INTEGER ::   ji, jk, jh  ! dummy loop indices
+      INTEGER ::   ios                  ! Local integer output status for namelist read
+      INTEGER ::   k2d, k3d             ! dummy number of analysis
+      NAMELIST/nam_diaharm_fast/ ln_diaharm_store, ln_diaharm_compute, ln_diaharm_read_restart, ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d, tname
+      !!----------------------------------------------------------------------
+
+      lk_diaharm_2D    = .TRUE.   ! to run 2d
+      lk_diaharm_3D    = .TRUE.   ! to run 3d
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'harm_init : initialization of harmonic analysis of tides'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
+
+      ! GET NAMELIST DETAILS
+      REWIND( numnam_ref )              ! Namelist nam_diaharm_fast in reference namelist : Tidal harmonic analysis
+      READ  ( numnam_ref, nam_diaharm_fast, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in reference namelist' )
+
+      REWIND( numnam_cfg )              ! Namelist nam_diaharm_fast in configuration namelist : Tidal harmonic analysis
+      READ  ( numnam_cfg, nam_diaharm_fast, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in configuration namelist' )
+      IF(lwm) WRITE ( numond, nam_diaharm_fast )
+
+      ! GET NUMBER OF HARMONIC TO ANALYSE - from diaharm.F90
+      nb_ana = 0
+      DO jk=1,jpmax_harmo
+         DO ji=1,nb_harmo
+            IF(TRIM(tname(jk)) == Wave( ntide(ji) )%cname_tide ) THEN
+               nb_ana=nb_ana+1
+            ENDIF
+         END DO
+      END DO
+      !
+      IF(lwp) THEN
+         WRITE(numout,*) '        Namelist nam_diaharm_fast'
+         WRITE(numout,*) '        nb_ana    = ', nb_ana
+         CALL flush(numout)
+      ENDIF
+      !
+      IF (nb_ana > nharm_max) THEN
+        IF(lwp) WRITE(numout,*) ' E R R O R harm_ana : nb_ana must be lower than nharm_max, stop'
+        IF(lwp) WRITE(numout,*) ' nharm_max = ', nharm_max
+        nstop = nstop + 1
+      ENDIF
+
+      ALLOCATE(ntide_all(nb_ana))
+      ALLOCATE(ntide_sub(nb_ana))
+
+      DO jk=1,nb_ana
+       DO ji=1,nb_harmo
+          IF (TRIM(tname(jk)) .eq. Wave( ntide(ji) )%cname_tide ) THEN
+             ntide_sub(jk) = ji
+             ntide_all(jk) = ntide(ji)
+             EXIT
+          END IF
+       END DO
+      END DO
+
+      ! SEARCH HOW MANY VARIABLES 2D AND 3D TO PROCESS
+      nvar_2d = 0; nvar_3d = 0
+      IF ( ln_ana_ssh   ) nvar_2d = nvar_2d + 1       ! analysis elevation
+      IF ( ln_ana_uvbar ) nvar_2d = nvar_2d + 2       ! analysis depth-averaged velocity
+      IF ( ln_ana_bfric ) nvar_2d = nvar_2d + 2       ! analysis bottom friction 
+            
+      IF ( ln_ana_rho   ) nvar_3d = nvar_3d + 1       ! analysis density
+      IF ( ln_ana_uv3d  ) nvar_3d = nvar_3d + 2       ! analysis 3D horizontal velocities
+      IF ( ln_ana_w3d   ) nvar_3d = nvar_3d + 1       ! analysis 3D vertical velocity
+
+      ! CHECK IF SOMETHING TO RUN
+      IF ( nvar_2d .eq. 0 ) lk_diaharm_2D = .FALSE.   ! no 2d to run
+      IF ( nvar_3d .eq. 0 ) lk_diaharm_3D = .FALSE.   ! no 3d to run
+
+      IF ( ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
+
+         ! DO ALLOCATIONS
+         IF ( lk_diaharm_2D ) THEN
+            ALLOCATE( g_cumul_var2D(nb_ana*2+1,jpi,jpj,    nvar_2d) )
+            ALLOCATE( g_cosamp2D( 0:nb_ana*2+1,jpi,jpj,    nvar_2d) )
+            ALLOCATE( g_sinamp2D( 0:nb_ana*2+1,jpi,jpj,    nvar_2d) )
+            ALLOCATE( g_out2D (jpi,jpj) )
+            ALLOCATE( h_out2D (jpi,jpj) )
+            ALLOCATE( m_posi_2d( nvar_2d ) ); m_posi_2d(:)=0
+         ENDIF
+ 
+         IF ( lk_diaharm_3D ) THEN
+            ALLOCATE( g_cumul_var3D(nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
+            ALLOCATE( g_cosamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
+            ALLOCATE( g_sinamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
+            ALLOCATE( g_out3D (jpi,jpj,jpk) )
+            ALLOCATE( h_out3D (jpi,jpj,jpk) )
+            ALLOCATE( m_posi_3d( nvar_3d ) ); m_posi_3d(:)=0
+         ENDIF
+
+         ALLOCATE( cc(nb_ana*2+1,nb_ana*2+1) )
+         ALLOCATE( a (nb_ana*2+1,nb_ana*2+1) )
+         ALLOCATE( bzz(nb_ana*2+1) )
+         ALLOCATE( x  (nb_ana*2+1) )
+         ALLOCATE( c  (nb_ana*2+1) )
+         ALLOCATE( anau(nb_ana) )
+         ALLOCATE( anav(nb_ana) )
+         ALLOCATE( anaf(nb_ana) )
+         ! END ALLOCATE 
+
+         ! STORE INDEX OF WHAT TO PRODUCE DEPENDING ON ACTIVATED LOGICAL
+         ! MAKES THINGS EASIER AND FASTER LATER
+         ! !!! UGLY !!!
+         jh = 1; k2d = 0; 
+         IF ( ln_ana_ssh   ) THEN
+            k2d = k2d + 1; m_posi_2d(k2d) = jh
+            IF(lwp) WRITE(numout,*) "   - ssh harmonic analysis activated (ln_ana_ssh)"
+         ENDIF
+         jh = jh + 1
+         IF ( ln_ana_uvbar ) THEN
+            k2d = k2d + 1; m_posi_2d(k2d) = jh
+            jh  = jh  + 1 
+            k2d = k2d + 1; m_posi_2d(k2d) = jh
+            IF(lwp) WRITE(numout,*) "   - barotropic currents harmonic analysis activated (ln_ana_uvbar)"
+         ELSE
+            jh  = jh  + 1
+         ENDIF
+         jh = jh + 1
+         IF ( ln_ana_bfric ) THEN
+            k2d = k2d + 1; m_posi_2d(k2d) = jh
+            jh  = jh  + 1; 
+            k2d = k2d + 1; m_posi_2d(k2d) = jh
+            IF(lwp) WRITE(numout,*) "   - bottom friction harmonic analysis activated (ln_ana_vbfr)"
+         ELSE
+            jh  = jh  + 1
+         ENDIF
+
+         ! and for 3D
+         jh = 1; k3d = 0; 
+         IF ( ln_ana_rho  ) THEN
+            k3d = k3d + 1; m_posi_3d(k3d) = jh
+            IF(lwp) WRITE(numout,*) "   - 3D density harmonic analysis activated (ln_ana_rho)"
+         ENDIF
+         jh = jh + 1
+         IF ( ln_ana_uv3d )  THEN
+            k3d = k3d + 1; m_posi_3d(k3d) = jh
+            jh  = jh  + 1 
+            k3d = k3d + 1; m_posi_3d(k3d) = jh
+            IF(lwp) WRITE(numout,*) "   - 3D horizontal currents harmonic analysis activated (ln_ana_uv3d)"
+         ELSE
+            jh  = jh  + 1
+         ENDIF
+         jh = jh + 1
+         IF ( ln_ana_w3d ) THEN
+            k3d = k3d + 1; m_posi_3d(k3d) = jh
+            IF(lwp) WRITE(numout,*) "   - 3D vertical currents harmonic analysis activated (ln_ana_w3d)"
+         ENDIF
+
+         ! SELECT AND STORE FREQUENCIES
+         IF(lwp)    WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
+         DO jh=1,nb_ana
+            om_tide(jh) = omega_tide( ntide_sub(jh) ) 
+            IF(lwp) WRITE(numout,*) '        - ',tname(jh),' ',om_tide(jh)
+         ENDDO
+
+         ! READ RESTART IF 
+         IF ( ln_diaharm_read_restart ) THEN
+            IF (lwp) WRITE(numout,*) "Reading previous harmonic data from previous run"
+            ! Need to read in  bssh bz, cc anau anav and anaf 
+            call harm_rst_read  ! This reads in from the previous day
+                                ! Currrently the data in in assci format
+         ELSE 
+
+            IF (lwp) WRITE(numout,*) "Starting harmonic analysis from Fresh "
+ 
+            IF ( lk_diaharm_2D ) g_cumul_var2D(:,:,:,:  ) = 0.0_wp
+            IF ( lk_diaharm_3D ) g_cumul_var3D(:,:,:,:,:) = 0.0_wp
+            cc           = 0.0_wp
+            a    (:,:)   = 0.0_wp ! NB
+            bzz  (:)     = 0.0_wp
+            x    (:)     = 0.0_wp
+            c    (:)     = 0.0_wp
+            anau (:)     = 0.0_wp
+            anav (:)     = 0.0_wp
+            anaf (:)     = 0.0_wp
+
+            DO jh = 1, nb_ana
+               anau(jh) = utide ( ntide_sub(jh) )
+               anav(jh) = v0tide( ntide_sub(jh) )
+               anaf(jh) = ftide ( ntide_sub(jh) )
+            END DO
+
+            fjulday_startharm=fjulday !Set this at very start and store
+
+            IF (lwp) THEN
+               WRITE(numout,*) '--------------------------'
+               WRITE(numout,*) '   - Output anaf for check'
+               WRITE(numout,*) 'ANA F', anaf
+               WRITE(numout,*) 'ANA U', anau
+               WRITE(numout,*) 'ANA V', anav
+               WRITE(numout,*) fjulday_startharm
+               WRITE(numout,*) '--------------------------'
+            ENDIF
+
+         ENDIF
+
+      ELSE
+
+         IF (lwp) WRITE(numout,*) "No variable setup for harmonic analysis"
+
+      ENDIF
+
+   END SUBROUTINE harm_ana_init
+!
+   SUBROUTINE gelim (a,b,x,n)
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE harm_ana  ***
+      !!
+      !! ** Purpose :   Guassian elimination
+      !!
+      !!
+      !! ** Action  : - first action (share memory array/varible modified
+      !!                in this routine
+      !!              - second action .....
+      !!              - .....
+      !!
+      !! References :
+      !!   Give references if exist otherwise suppress these lines
+      !!
+      !! History :
+        implicit none
+!
+        integer  :: n
+        REAL(WP) :: b(nb_ana*2+1), a(nb_ana*2+1,nb_ana*2+1)
+        REAL(WP) :: x(nb_ana*2+1)
+        INTEGER  :: row,col,prow,pivrow,rrow
+        REAL(WP) :: atemp
+        REAL(WP) :: pivot
+        REAL(WP) :: m
+
+        do row=1,n-1
+           pivrow=row
+           pivot=a(row,n-row+1)
+           do prow=row+1,n
+              if (abs(a(prow,n-row+1)).gt.abs(pivot)  ) then
+                 pivot=a(prow,n-row+1)
+                 pivrow=prow
+              endif
+           enddo
+!	swap row and prow
+           if ( pivrow .ne. row ) then
+              atemp=b(pivrow)
+              b(pivrow)=b(row)
+              b(row)=atemp
+              do col=1,n
+                 atemp=a(pivrow,col)
+                 a(pivrow,col)=a(row,col)
+                 a(row,col)=atemp
+              enddo
+           endif
+
+           do rrow=row+1,n
+              if (a(row,row).ne.0) then
+   
+                 m=-a(rrow,n-row+1)/a(row,n-row+1)
+                 do col=1,n
+                    a(rrow,col)=m*a(row,col)+a(rrow,col)
+                 enddo
+                 b(rrow)=m*b(row)+b(rrow)
+              endif
+           enddo
+        enddo
+!	back substitution now
+
+        x(1)=b(n)/a(n,1)
+        do row=n-1,1,-1
+           x(n-row+1)=b(row)
+           do col=1,(n-row)
+              x(n-row+1)=(x(n-row+1)-a(row,col)*x(col)) 
+           enddo
+
+           x(n-row+1)=(x(n-row+1)/a(row,(n-row)+1))
+        enddo
+
+        return
+   END SUBROUTINE gelim
+
+   SUBROUTINE harm_ana_out
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE harm_ana_init  ***
+      !!                   
+      !! ** Purpose :   initialization of ....
+      !!
+      !! ** Method  :   blah blah blah ...
+      !!
+      !! ** input   :   Namlist namexa
+      !!
+      !! ** Action  :   ...  
+      !!
+      !! history :
+      !!   9.0  !  03-08  (Autor Names)  Original code
+      !!----------------------------------------------------------------------
+        USE dianam          ! build name of file (routine)
+ 
+      !! * local declarations
+      INTEGER :: ji, jj, jk, jgrid, jh    ! dummy loop indices
+!      INTEGER :: nh_T
+!      INTEGER :: nid_harm
+!      CHARACTER (len=40) :: clhstnamt, clop1, clop2 ! temporary names 
+!      CHARACTER (len=40) :: clhstnamu, clhstnamv    ! temporary names 
+      CHARACTER (len=40) :: suffix
+!      REAL(wp) :: zsto1, zsto2, zout, zmax, zjulian, zdt, zmdi  ! temporary scalars
+
+      do jgrid=1,nvar_2d
+          do jh=1,nb_ana
+             h_out2D = 0.0
+             g_out2D = 0.0
+             do jj=1,nlcj
+                do ji=1,nlci
+                   cca=g_cosamp2D(jh,ji,jj,jgrid)
+                   ssa=g_sinamp2D(jh,ji,jj,jgrid)
+                   h_out2D(ji,jj)=sqrt(cca**2+ssa**2)
+                   IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN 
+                      g_out2D(ji,jj)= 0.0_wp
+                   ELSE
+                      g_out2D(ji,jj)=(180.0/rpi)*atan2(ssa,cca)       
+                   ENDIF 
+                   IF (h_out2D(ji,jj).ne.0) THEN
+                       h_out2D(ji,jj)=h_out2D(ji,jj)/anaf(jh)
+                   ENDIF
+                   IF (g_out2D(ji,jj).ne.0) THEN  !Correct and take modulus
+                       g_out2D(ji,jj) = g_out2D(ji,jj) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
+                       if (g_out2D(ji,jj).gt.360.0) then
+                           g_out2D(ji,jj)=g_out2D(ji,jj)-360.0
+                       else if (g_out2D(ji,jj).lt.0.0) then
+                           g_out2D(ji,jj)=g_out2D(ji,jj)+360.0
+                       endif
+                   ENDIF
+                enddo
+             enddo
+             !
+             ! NETCDF OUTPUT
+             suffix = TRIM( m_varName2d( m_posi_2d(jgrid) ) )
+             CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix), h_out2D(:,:) )
+             CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix), g_out2D(:,:) )
+
+          enddo
+      enddo
+!
+! DO THE SAME FOR 3D VARIABLES
+!
+      do jgrid=1,nvar_3d
+          do jh=1,nb_ana
+             h_out3D = 0.0
+             g_out3D = 0.0
+             DO jk=1,jpkm1
+                do jj=1,nlcj
+                   do ji=1,nlci
+                      cca=g_cosamp3D(jh,ji,jj,jk,jgrid)
+                      ssa=g_sinamp3D(jh,ji,jj,jk,jgrid)
+                      h_out3D(ji,jj,jk)=sqrt(cca**2+ssa**2)
+                      IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
+                         g_out3D(ji,jj,jk) = 0.0_wp
+                      ELSE
+                         g_out3D(ji,jj,jk) = (180.0/rpi)*atan2(ssa,cca)
+                      ENDIF
+                      IF (h_out3D(ji,jj,jk).ne.0) THEN
+                          h_out3D(ji,jj,jk) = h_out3D(ji,jj,jk)/anaf(jh)
+                      ENDIF
+                      IF (g_out3D(ji,jj,jk).ne.0) THEN  !Correct and take modulus
+                          g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
+                          if      (g_out3D(ji,jj,jk).gt.360.0) then
+                                   g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)-360.0
+                          else if (g_out3D(ji,jj,jk).lt.0.0) then
+                                   g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)+360.0
+                          endif
+                      ENDIF
+                   enddo    ! ji
+                enddo       ! jj
+             ENDDO          ! jk
+             !
+             ! NETCDF OUTPUT
+             suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) )
+             IF(lwp) WRITE(numout,*) "harm_ana_out", suffix
+             CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix), h_out3D(:,:,:) )
+             CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix), g_out3D(:,:,:) )
+          enddo             ! jh 
+      enddo                 ! jgrid
+!
+   END SUBROUTINE harm_ana_out
+!
+   SUBROUTINE harm_rst_write(kt)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE harm_ana_init  ***
+      !!                   
+      !! ** Purpose :  To write out cummulated Tidal Harmomnic data to file for
+      !!               restarting
+      !!
+      !! ** Method  :   restart files will be dated by default
+      !!
+      !! ** input   :   
+      !!
+      !! ** Action  :   ...  
+      !!
+      !! history :
+      !!   0.0  !  01-16  (Enda O'Dea)  Original code
+      !! ASSUMES  dated file for rose  , can change later to be more generic
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt     ! ocean time-step
+      !!
+      INTEGER             ::   jh, j2d, j3d
+      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
+      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
+      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file
+      CHARACTER(LEN=250)  ::   clfinal   ! full name
+
+      !restart file
+      DO j2d=1,nvar_2d
+         CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
+         DO jh=1,nb_ana
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2  , :, :, j2d ) )
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) )
+         ENDDO
+      ENDDO
+
+      DO j3d=1,nvar_3d
+         CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
+         DO jh=1,nb_ana
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2  , :, :, :, j3d ) )
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) )
+         ENDDO
+      ENDDO
+
+      IF(lwp) THEN
+        IF( kt > 999999999 ) THEN ; WRITE(clkt, *       ) kt
+        ELSE                      ; WRITE(clkt, '(i8.8)') kt
+        ENDIF
+        clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
+        clpath = TRIM(cn_ocerst_outdir)
+        IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
+        IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for writing: ',TRIM(clpath)//clname
+
+        WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
+        OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
+        WRITE(66) cc
+        WRITE(66) anau
+        WRITE(66) anav
+        WRITE(66) anaf
+        WRITE(66) fjulday_startharm
+        CLOSE(66)
+        WRITE(numout,*) '----------------------------'
+        WRITE(numout,*) '   harm_rst_write: DONE '
+        WRITE(numout,*) cc
+        WRITE(numout,*) anaf
+        WRITE(numout,*) fjulday_startharm
+        WRITE(numout,*) '----------------------------'
+      ENDIF
+ 
+   END SUBROUTINE harm_rst_write
+
+   SUBROUTINE harm_rst_read
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE harm_ana_init  ***
+      !!                   
+      !! ** Purpose :  To read in  cummulated Tidal Harmomnic data to file for
+      !!               restarting
+      !!
+      !! ** Method  :   
+      !!
+      !! ** input   :   
+      !!
+      !! ** Action  :   ...  
+      !!
+      !! history :
+      !!   0.0  !  01-16  (Enda O'Dea)  Original code
+      !! ASSUMES  dated file for rose  , can change later to be more generic
+      !!----------------------------------------------------------------------
+      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
+      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
+      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file
+      CHARACTER(LEN=250)  ::   clfinal   ! full name
+      INTEGER             ::   jh, j2d, j3d
+
+      IF( nit000 > 999999999 ) THEN ; WRITE(clkt, *       ) nit000-1
+      ELSE                      ; WRITE(clkt, '(i8.8)') nit000-1
+      ENDIF
+      clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
+      clpath = TRIM(cn_ocerst_outdir)
+      IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
+
+      IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for reading: ',TRIM(clpath)//clname
+
+      DO j2d=1,nvar_2d
+         CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
+         IF(lwp) WRITE(numout,*) "2D", j2d, m_posi_2d(j2d), m_varName2d( m_posi_2d(j2d) )
+         DO jh=1,nb_ana
+            CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2  , :, :, j2d ) )
+            CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) )
+         ENDDO
+      ENDDO
+
+      DO j3d=1,nvar_3d
+         CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
+         IF(lwp) WRITE(numout,*) "3D", j3d,  m_posi_3d(j3d), m_varName3d( m_posi_3d(j3d) )
+
+         DO jh=1,nb_ana
+            CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2  , :, :, :, j3d ) )
+            CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) )
+         ENDDO
+      ENDDO
+
+      WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
+      OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
+      READ(66) cc
+      READ(66) anau
+      READ(66) anav
+      READ(66) anaf
+      READ(66) fjulday_startharm
+      CLOSE(66)
+
+      IF(lwp) THEN
+        WRITE(numout,*) '----------------------------'
+        WRITE(numout,*) '   Checking anaf is correct'
+        WRITE(numout,*) cc
+        WRITE(numout,*) anaf
+        WRITE(numout,*) fjulday_startharm
+        WRITE(numout,*) '----------------------------'
+      ENDIF
+ 
+   END SUBROUTINE harm_rst_read
+
+   !!======================================================================
+#else
+!!---------------------------------------------------------------------------------
+!!   Dummy module                                   NO harmonic Analysis
+!!---------------------------------------------------------------------------------
+        LOGICAL, PUBLIC, PARAMETER :: lk_diaharm_fast  = .FALSE.   ! to be run or not
+
+        CONTAINS
+           SUBROUTINE harm_rst_write(kt)     ! Dummy routine
+           END SUBROUTINE harm_rst_write
+           SUBROUTINE harm_rst_read    ! Dummy routine
+           END SUBROUTINE harm_rst_read
+           SUBROUTINE harm_ana_out      ! Dummy routine
+           END SUBROUTINE harm_ana_out
+           SUBROUTINE harm_ana_init
+           END SUBROUTINE harm_ana_init
+           SUBROUTINE harm_ana( kt )
+!--- NB : end call not properly written
+           END SUBROUTINE harm_ana
+!           END SUBROUTINE harm_ana_init
+!--- END NB
+           SUBROUTINE gelim (a,b,x,n)
+!--- NB : end call not properly written
+           END SUBROUTINE gelim
+!           END SUBROUTINE gelim (a,b,x,n)
+!--- END NB           
+#endif
+
+END MODULE diaharm_fast 
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DOM/dtatsd.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DOM/dtatsd.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DOM/dtatsd.F90	(revision 15422)
@@ -21,4 +21,5 @@
    !
    USE in_out_manager  ! I/O manager
+   USE iom             ! IOM library
    USE lib_mpp         ! MPP library
 
@@ -31,7 +32,11 @@
    !                                  !!* 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
 
    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)
 
    !!----------------------------------------------------------------------
@@ -53,9 +58,9 @@
       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
@@ -63,5 +68,5 @@
       !
       !  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 +85,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,23 +100,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
+         IF( ln_tsd_interp ) THEN  
+            CALL iom_open ( trim(cn_dir) // trim(sn_dep%clname), inum_dta )   
+            id = iom_varid( inum_dta, sn_dep%clvar, zdim )  
+            jpk_init = zdim(3)  
+            IF(lwp) WRITE(numout,*) 'Dimension of vertical coordinate in ICs: ', jpk_init  
+            CALL iom_close( inum_dta )   ! Close the input file  
+            !  
+                                 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
          !                         ! fill sf_tsd with sn_tem & sn_sal and control print
          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', no_print )
          !
@@ -138,8 +170,7 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data
       !
-      INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies
-      INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1   ! local integers
+      INTEGER ::   ji, jj, jk, jl, jk_init   ! dummy loop indices  
+      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 +207,47 @@
 !!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( 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
+      IF( kt == nit000 .AND. lwp )THEN  
+         WRITE(numout,*)  
+         WRITE(numout,*) 'dta_tsd: interpolates T & S data onto current mesh'  
+      ENDIF  
+      !  
+      IF( ln_tsd_interp ) THEN                 ! probably should use pointers in the following to make more readable  
+      !  
+         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
          END DO
          ! 
+         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) *tmask(:,:,:)  
+         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) *tmask(:,:,:)
       ELSE                                !==   z- or zps- coordinate   ==!
          !                             
-         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) * tmask(:,:,:)    ! Mask
-         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) * tmask(:,:,:)
+         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 +279,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/r14075_India_uncoupled/src/OCE/DYN/dynnxt.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DYN/dynnxt.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DYN/dynnxt.F90	(revision 15422)
@@ -33,4 +33,5 @@
    USE dynadv         ! dynamics: vector invariant versus flux form
    USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
+   USE dynspg
    USE domvvl         ! variable volume
    USE bdy_oce , ONLY : ln_bdy
@@ -58,4 +59,6 @@
    PUBLIC    dyn_nxt   ! routine called by step.F90
 
+   !! Substitution 
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -182,4 +185,8 @@
             vn(:,:,jk) = va(:,:,jk)
          END DO
+         ! limit velocities  
+         IF (ln_ulimit) THEN  
+            call dyn_limit_velocity (kt)  
+         ENDIF
          IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
 !!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a  
@@ -214,4 +221,8 @@
                END DO
             END DO
+            ! limit velocities  
+            IF (ln_ulimit) THEN  
+               call dyn_limit_velocity (kt)  
+            ENDIF
             !                             ! ================!
          ELSE                             ! Variable volume !
@@ -263,4 +274,8 @@
                   END DO
                END DO
+               ! limit velocities  
+               IF (ln_ulimit) THEN  
+                  call dyn_limit_velocity (kt)  
+               ENDIF
                !
             ELSE                          ! Asselin filter applied on thickness weighted velocity
@@ -290,4 +305,8 @@
                   END DO
                END DO
+               ! limit velocities  
+               IF (ln_ulimit) THEN  
+                  call dyn_limit_velocity (kt)  
+               ENDIF
                e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
                e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
@@ -401,4 +420,59 @@
    END SUBROUTINE dyn_nxt
 
+   SUBROUTINE dyn_limit_velocity (kt)  
+      ! limits maximum values of un and vn by volume courant number  
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index  
+      !  
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices  
+      REAL(wp) ::   zzu,zplim,zmlim,isp,ism,zcn,ze3e1,zzcn,zcnn,idivp,idivm  
+   
+      ! limit fluxes  
+      zcn =cn_ulimit !0.9 ! maximum velocity inverse courant number  
+      zcnn = cnn_ulimit !0.54 ! how much to reduce cn by in divergen flow  
+   
+      DO jk = 1, jpkm1  
+         DO jj = 2, jpjm1  
+            DO ji = fs_2, fs_jpim1   ! vect. opt. 
+               ! U direction  
+               zzu = un(ji,jj,jk)  
+               ze3e1 = e3u_n(ji  ,jj,jk) * e2u(ji  ,jj)  
+               ! ips is 1 if flow is positive othersize zero  
+               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )  
+               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )  
+               ! idev is 1 if divergent flow otherwise zero  
+               idivp = -isp * 0.5 * (sign(1.0_wp, un(ji-1,jj,jk)) - 1.0_wp )  
+               idivm =  ism * 0.5 * (sign(1.0_wp, un(ji+1,jj,jk)) + 1.0_wp )  
+               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp  
+               zzcn = zzcn * zcn  
+               zplim =  zzcn * (e3t_n(ji  ,jj,jk) * e1t(ji  ,jj) * e2t(ji  ,jj)) / &  
+                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)  
+               zmlim = -zzcn * (e3t_n(ji+1,jj,jk) * e1t(ji+1,jj) * e2t(ji+1,jj)) / &  
+                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)  
+               ! limit currents  
+               un(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism  
+   
+               ! V  direction  
+               zzu = vn(ji,jj,jk)  
+               ze3e1 = e3v_n(ji  ,jj,jk) * e1v(ji  ,jj)  
+               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )  
+               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )  
+               ! idev is 1 if divergent flow otherwise zero  
+               idivp = -isp * 0.5 * (sign(1.0_wp, vn(ji,jj-1,jk)) - 1.0_wp )  
+               idivm =  ism * 0.5 * (sign(1.0_wp, vn(ji,jj+1,jk)) + 1.0_wp )  
+               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp  
+               zzcn = zzcn * zcn  
+               zplim =  zzcn * (e3t_n(ji,jj  ,jk) * e1t(ji,jj  ) * e2t(ji,jj  )) / &  
+                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)  
+               zmlim = -zzcn * (e3t_n(ji,jj+1,jk) * e1t(ji,jj+1) * e2t(ji,jj+1)) / &  
+                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)  
+               ! limit currents  
+               vn(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism  
+            ENDDO  
+         ENDDO  
+      ENDDO  
+       CALL lbc_lnk_multi( 'dynnxt', un(:,:,:), 'U',  -1., vn(:,:,:), 'V',  -1. ) 
+   
+    END SUBROUTINE dyn_limit_velocity
+
    !!=========================================================================
 END MODULE dynnxt
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DYN/dynspg.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DYN/dynspg.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DYN/dynspg.F90	(revision 15422)
@@ -39,4 +39,7 @@
    INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_... 
 
+   LOGICAL, PUBLIC ::  ln_ulimit  
+   REAL(wp), PUBLIC :: cn_ulimit,cnn_ulimit
+
    !                       ! Parameter to control the surface pressure gradient scheme
    INTEGER, PARAMETER ::   np_TS  = 1   ! split-explicit time stepping (Time-Splitting)
@@ -191,5 +194,6 @@
       NAMELIST/namdyn_spg/ ln_dynspg_exp       , ln_dynspg_ts,   &
       &                    ln_bt_fw, ln_bt_av  , ln_bt_auto  ,   &
-      &                    nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha
+      &                    nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha, & 
+      &                    ln_ulimit, cn_ulimit, cnn_ulimit
       !!----------------------------------------------------------------------
       !
@@ -213,4 +217,7 @@
          WRITE(numout,*) '      Explicit free surface                  ln_dynspg_exp = ', ln_dynspg_exp
          WRITE(numout,*) '      Free surface with time splitting       ln_dynspg_ts  = ', ln_dynspg_ts
+         WRITE(numout,*) '     Limit velocities ln_ulimit     = ', ln_ulimit  
+         WRITE(numout,*) '     Limit velocities      max inverse Courant number     = ', cn_ulimit  
+         WRITE(numout,*) '     Limit velocities   multiplier for divergant flow     = ', cnn_ulimit
       ENDIF
       !                          ! Control of surface pressure gradient scheme options
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/fldread.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/fldread.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/fldread.F90	(revision 15422)
@@ -901,6 +901,6 @@
             IF(     zdepth(jk) <= pdta_read_z(jb,1,1)      ) THEN   ! above the first level of external data
                pdta(jb,1,jk) = pdta_read(jb,1,1)
-            ELSEIF( zdepth(jk) >  pdta_read_z(jb,1,ipkmax) ) THEN   ! below the last level of external data /= FillValue
-               pdta(jb,1,jk) = pdta_read(jb,1,ipkmax)
+            ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkmax) ) THEN    ! below the last level of external data  
+               pdta(jb,1,jk) =  pdta_read(jb,1,ipkmax)
             ELSE                                                    ! inbetween: vertical interpolation between jkb & jkb+1
                DO jkb = 1, ipkmax-1
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/sbctide.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/sbctide.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/sbctide.F90	(revision 15422)
@@ -16,4 +16,5 @@
    USE ioipsl         ! NetCDF IPSL library
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE bdytides
 
    IMPLICIT NONE
@@ -108,5 +109,6 @@
 
       DO jk = 1, nb_harmo
-         zcons = 0.7_wp * Wave(ntide(jk))%equitide * ftide(jk)
+         ! love number now provides in tide namelist  
+         zcons = dn_love_number * Wave(ntide(jk))%equitide * ftide(jk)
          DO ji = 1, jpi
             DO jj = 1, jpj
@@ -119,4 +121,6 @@
                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
+               ! Add tide potential for long period tides  
+               ELSEIF( Wave(ntide(jk))%nutide == 0 )  THEN  ;  zcs = zcons * (0.5_wp-1.5_wp*SIN(zlat)**2._wp)
                ELSE                                         ;  zcs = 0._wp
                ENDIF
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/tide_FES14.h90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/tide_FES14.h90	(revision 15422)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/tide_FES14.h90	(revision 15422)
@@ -0,0 +1,114 @@
+   !!----------------------------------------------------------------------
+   !! 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 !!
+   !             !!           !          !        !    !    !    !    !     !       !      !      !      !      !   !         !!
+   !
+   ! 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/r14075_India_uncoupled/src/OCE/SBC/tide_mod.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/tide_mod.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/tide_mod.F90	(revision 15422)
@@ -16,5 +16,5 @@
    PUBLIC   tide_init_Wave   ! called by tideini and diaharm modules
 
-   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 19   !: maximum number of harmonic
+   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 31   !: maximum number of harmonic
 
    TYPE, PUBLIC ::    tide
@@ -41,5 +41,9 @@
 
    SUBROUTINE tide_init_Wave
+# if defined key_FES14_tides  
+#     include "tide_FES14.h90"  
+# else
 #     include "tide.h90"
+# endif
    END SUBROUTINE tide_init_Wave
 
@@ -331,4 +335,8 @@
          zf = zf * zf1 * zf1
          !
+      CASE( 20 )                 !==  formule 20,  compound waves ( 78 x 78 x 78 x 78 )  
+         zf1 = nodal_factort(78)  
+         zf  = zf1 * zf1 * zf1 * zf1  
+         ! 
       CASE( 73 )                 !==  formule 73
          zs = sin(sh_I)
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/tideini.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/tideini.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/tideini.F90	(revision 15422)
@@ -33,4 +33,5 @@
    INTEGER , PUBLIC ::   kt_tide         !:
    REAL(wp), PUBLIC ::   rdttideramp     !:
+   REAL(wp), PUBLIC ::   dn_love_number  !:
    REAL(wp), PUBLIC ::   rn_scal_load    !:
    CHARACTER(lc), PUBLIC ::   cn_tide_load   !: 
@@ -54,5 +55,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
       !!----------------------------------------------------------------------
       !
@@ -78,4 +79,5 @@
             WRITE(numout,*) '         Read load potential from file           ln_read_load = ', ln_read_load
             WRITE(numout,*) '         Apply ramp on tides at startup          ln_tide_ramp = ', ln_tide_ramp
+            WRITE(numout,*) '                                              dn_love_number  = ', dn_love_number
             WRITE(numout,*) '         Fraction of SSH used in scal. approx.   rn_scal_load = ', rn_scal_load
             WRITE(numout,*) '         Duration (days) of ramp                 rdttideramp  = ', rdttideramp
@@ -99,4 +101,6 @@
       END DO
       !       
+      IF (ln_tide .and.lwp) WRITE(numout,*) 'nb_harmo     = ', nb_harmo
+
       ! Ensure that tidal components have been set in namelist_cfg
       IF( nb_harmo == 0 )   CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' )
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/USR/usrdef_istate.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/USR/usrdef_istate.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/USR/usrdef_istate.F90	(revision 15422)
@@ -62,18 +62,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               pts(ji,jj,jk,jp_tem) =  (  (  16. - 12. * TANH( (pdept(ji,jj,jk) - 400) / 700 ) )   &
-                    &           * (-TANH( (500. - pdept(ji,jj,jk)) / 150. ) + 1.) / 2.             &
-                    &           + ( 15. * ( 1. - TANH( (pdept(ji,jj,jk)-50.) / 1500.) )            &
-                    &           - 1.4 * TANH((pdept(ji,jj,jk)-100.) / 100.)                        &
-                    &           + 7.  * (1500. - pdept(ji,jj,jk) ) / 1500.)                        &
-                    &           * (-TANH( (pdept(ji,jj,jk) - 500.) / 150.) + 1.) / 2.  ) * ptmask(ji,jj,jk)
-
-               pts(ji,jj,jk,jp_sal) =  (  (  36.25 - 1.13 * TANH( (pdept(ji,jj,jk) - 305) / 460 ) )  &
-                    &         * (-TANH((500. - pdept(ji,jj,jk)) / 150.) + 1.) / 2                  &
-                    &         + ( 35.55 + 1.25 * (5000. - pdept(ji,jj,jk)) / 5000.                 &
-                    &         - 1.62 * TANH( (pdept(ji,jj,jk) - 60.  ) / 650. )                    &
-                    &         + 0.2  * TANH( (pdept(ji,jj,jk) - 35.  ) / 100. )                    &
-                    &         + 0.2  * TANH( (pdept(ji,jj,jk) - 1000.) / 5000.) )                  &
-                    &         * (-TANH( (pdept(ji,jj,jk) - 500.) / 150.) + 1.) / 2  ) * ptmask(ji,jj,jk)
+               pts(ji,jj,jk,jp_tem) = 20._wp * ptmask(ji,jj,jk)  
+               pts(ji,jj,jk,jp_sal) = 36.25_wp * ptmask(ji,jj,jk)
             END DO
          END DO
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/USR/usrdef_sbc.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/USR/usrdef_sbc.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/USR/usrdef_sbc.F90	(revision 15422)
@@ -1,32 +1,44 @@
 MODULE usrdef_sbc
    !!======================================================================
-   !!                     ***  MODULE  usrdef_sbc  ***
-   !!
-   !!                     ===  GYRE configuration  ===
+   !!                       ***  MODULE usrdef_sbc  *** 
+   !!  
+   !!                  ===  AMM7_SURGE configuration  ===
    !!
    !! User defined :   surface forcing of a user configuration
    !!======================================================================
    !! History :  4.0   ! 2016-03  (S. Flavoni, G. Madec)  user defined interface
+   !!            4.0   ! 2017-12  (C. O'Neill)  add necessary options for surge work - either no fluxes  
+   !!                                           (for tide-only run) or wind and pressure only
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   usrdef_sbc    : user defined surface bounday conditions in GYRE case
+   !!   usrdef_sbc    : user defined surface bounday conditions in LOCK_EXCHANGE case
    !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and tracers
-   USE dom_oce        ! ocean space and time domain
-   USE sbc_oce        ! Surface boundary condition: ocean fields
-   USE phycst         ! physical constants
+   USE oce             ! ocean dynamics and tracers 
+   USE dom_oce         ! ocean space and time domain 
+   USE sbc_oce         ! Surface boundary condition: ocean fields 
+   USE sbc_ice         ! Surface boundary condition: ocean fields 
+   USE fldread         ! read input fields 
+   USE phycst          ! physical constants 
+   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
    !
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! distribued memory computing library
-   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE lib_fortran    !
+   USE in_out_manager  ! I/O manager 
+   USE iom 
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
+   USE lib_mpp         ! distribued memory computing library 
+   !USE wrk_nemo       ! work arrays 
+   USE timing         ! Timing 
+   USE prtctl         ! Print control
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   usrdef_sbc_oce       ! routine called in sbcmod module
-   PUBLIC   usrdef_sbc_ice_tau   ! routine called by icestp.F90 for ice dynamics
-   PUBLIC   usrdef_sbc_ice_flx   ! routine called by icestp.F90 for ice thermo
+   PUBLIC   usrdef_sbc_oce    ! routine called in sbcmod module 
+   PUBLIC   usrdef_sbc_ice_tau  ! routine called by sbcice_lim.F90 for ice dynamics 
+   PUBLIC   usrdef_sbc_ice_flx  ! routine called by sbcice_lim.F90 for ice thermo 
+   !                                  !!* Namelist namsbc_usr 
+   REAL(wp) ::   rn_vfac     ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem) 
+   REAL(wp) ::   rn_charn_const  
+   LOGICAL  ::   ln_use_sbc  ! Surface fluxes on or not
 
    !! * Substitutions
@@ -43,182 +55,46 @@
       !!                    ***  ROUTINE usrdef_sbc  ***
       !!              
-      !! ** Purpose :   provide at each time-step the GYRE surface boundary
+      !! ** Purpose :   provide at each time-step the surface boundary
       !!              condition, i.e. the momentum, heat and freshwater fluxes.
       !!
-      !! ** Method  :   analytical seasonal cycle for GYRE configuration.
+      !! ** Method  :   all 0 fields, for AMM7_SURGE case
       !!                CAUTION : never mask the surface stress field !
       !!
-      !! ** Action  : - set the ocean surface boundary condition, i.e.   
+      !! ** Action  : - if tide-only case - set to ZERO all the ocean surface boundary condition, i.e.
       !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
-      !!
-      !! Reference : Hazeleger, W., and S. Drijfhout, JPO, 30, 677-695, 2000.
+      !!              - if tide+surge case - read in wind and air pressure      !!
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
-      !!
-      INTEGER  ::   ji, jj                 ! dummy loop indices
-      INTEGER  ::   zyear0                 ! initial year 
-      INTEGER  ::   zmonth0                ! initial month
-      INTEGER  ::   zday0                  ! initial day
-      INTEGER  ::   zday_year0             ! initial day since january 1st
-      REAL(wp) ::   ztau     , ztau_sais   ! wind intensity and of the seasonal cycle
-      REAL(wp) ::   ztime                  ! time in hour
-      REAL(wp) ::   ztimemax , ztimemin    ! 21th June, and 21th decem. if date0 = 1st january
-      REAL(wp) ::   ztimemax1, ztimemin1   ! 21th June, and 21th decem. if date0 = 1st january
-      REAL(wp) ::   ztimemax2, ztimemin2   ! 21th June, and 21th decem. if date0 = 1st january
-      REAL(wp) ::   ztaun                  ! intensity
-      REAL(wp) ::   zemp_s, zemp_n, zemp_sais, ztstar
-      REAL(wp) ::   zcos_sais1, zcos_sais2, ztrp, zconv, t_star
-      REAL(wp) ::   zsumemp, zsurf
-      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
-      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
-      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
-      REAL(wp) ::   zyydd                 ! number of days in one year
+
+      INTEGER  ::   ios      ! Local integer output status for namelist read 
+      ! 
+      CHARACTER(len=100) ::  cn_dir   !   Root directory for location of flux files 
+      TYPE(FLD_N) ::   sn_wndi, sn_wndj                        ! informations about the fields to be read 
+
+      NAMELIST/namsbc_usr/ ln_use_sbc, cn_dir , rn_vfac,  & 
+         &                   sn_wndi, sn_wndj, rn_charn_const
       !!---------------------------------------------------------------------
-      zyydd = REAL(nyear_len(1),wp)
+      ! 
+      IF( kt == nit000 ) THEN 
+          
+          
+         REWIND( numnam_cfg )              ! Namelist namsbc_usr in configuration namelist 
+         READ  ( numnam_cfg, namsbc_usr, IOSTAT = ios, ERR = 902 ) 
+902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_surge in configuration namelist' )
 
-      ! ---------------------------- !
-      !  heat and freshwater fluxes  !
-      ! ---------------------------- !
-      !same temperature, E-P as in HAZELEGER 2000
+         IF(lwm) WRITE( numond, namsbc_usr )
+         IF(lwp) WRITE(numout,*)' usr_sbc : AMM7_SURGE tide only case: NO surface forcing' 
+         IF(lwp) WRITE(numout,*)' ~~~~~~~~~~~   utau = vtau = taum = wndm = qns = qsr = emp = sfx = 0'
 
-      zyear0     =   ndate0 / 10000                             ! initial year
-      zmonth0    = ( ndate0 - zyear0 * 10000 ) / 100            ! initial month
-      zday0      =   ndate0 - zyear0 * 10000 - zmonth0 * 100    ! initial day betwen 1 and 30
-      zday_year0 = ( zmonth0 - 1 ) * 30.+zday0                  ! initial day betwen 1 and 360
-
-      ! current day (in hours) since january the 1st of the current year
-      ztime = REAL( kt ) * rdt / (rmmss * rhhmm)   &       !  total incrementation (in hours)
-         &      - (nyear  - 1) * rjjhh * zyydd             !  minus years since beginning of experiment (in hours)
-
-      ztimemax1 = ((5.*30.)+21.)* 24.                      ! 21th june     at 24h in hours
-      ztimemin1 = ztimemax1 + rjjhh * zyydd / 2            ! 21th december        in hours
-      ztimemax2 = ((6.*30.)+21.)* 24.                      ! 21th july     at 24h in hours
-      ztimemin2 = ztimemax2 - rjjhh * zyydd / 2            ! 21th january         in hours
-      !                                                    ! NB: rjjhh * zyydd / 4 = one seasonal cycle in hours
-
-      ! amplitudes
-      zemp_S    = 0.7       ! intensity of COS in the South
-      zemp_N    = 0.8       ! intensity of COS in the North
-      zemp_sais = 0.1
-      zTstar    = 28.3      ! intemsity from 28.3 a -5 deg
-
-      ! 1/2 period between 21th June and 21th December and between 21th July and 21th January
-      zcos_sais1 = COS( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi ) 
-      zcos_sais2 = COS( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi )
-
-      ztrp= - 40.e0        ! retroaction term on heat fluxes (W/m2/K)
-      zconv = 3.16e-5      ! convertion factor: 1 m/yr => 3.16e-5 mm/s
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3
-            ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
-            ! 64.5 in summer, 42.5 in winter
-            t_star = zTstar * ( 1. + 1. / 50. * zcos_sais2 )                &
-               &                    * COS( rpi * (gphit(ji,jj) - 5.)               &
-               &                    / ( 53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) )
-            ! 23.5 deg : tropics
-            qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )
-            qns (ji,jj) = ztrp * ( tsb(ji,jj,1,jp_tem) - t_star ) - qsr(ji,jj)
-            IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj) ) THEN    ! zero at 37.8 deg, max at 24.6 deg
-               emp  (ji,jj) =   zemp_S * zconv   &
-                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  &
-                  &         * ( 1 - zemp_sais / zemp_S * zcos_sais1)
-            ELSE
-               emp (ji,jj) =  - zemp_N * zconv   &
-                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  &
-                  &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 )
-            ENDIF
-         END DO
-      END DO
-
-      zsumemp = GLOB_SUM( 'usrdef_sbc', emp  (:,:)   ) 
-      zsurf   = GLOB_SUM( 'usrdef_sbc', tmask(:,:,1) ) 
-      zsumemp = zsumemp / zsurf         ! Default GYRE configuration
-
-      ! freshwater (mass flux) and update of qns with heat content of emp
-      emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)        ! freshwater flux (=0 in domain average)
-      sfx (:,:) = 0.0_wp                                   ! no salt flux
-      qns (:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp   ! evap and precip are at SST
-
-
-      ! ---------------------------- !
-      !       momentum fluxes        !
-      ! ---------------------------- !
-      ! same wind as in Wico
-      !test date0 : ndate0 = 010203
-      zyear0  =   ndate0 / 10000
-      zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100
-      zday0   =   ndate0 - zyear0 * 10000 - zmonth0 * 100
-      !Calculates nday_year, day since january 1st
-      zday_year0 = (zmonth0-1)*30.+zday0
-
-      !accumulates days of previous months of this year
-      ! day (in hours) since january the 1st
-      ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm)  &  ! incrementation in hour
-         &     - (nyear - 1) * rjjhh * zyydd          !  - nber of hours the precedent years
-      ztimemax = ((5.*30.)+21.)* 24.               ! 21th june     in hours
-      ztimemin = ztimemax + rjjhh * zyydd / 2      ! 21th december in hours
-      !                                            ! NB: rjjhh * zyydd / 4 = 1 seasonal cycle in hours
-
-      ! mean intensity at 0.105 ; srqt(2) because projected with 45deg angle
-      ztau = 0.105 / SQRT( 2. )
-      ! seasonal oscillation intensity
-      ztau_sais = 0.015
-      ztaun = ztau - ztau_sais * COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi )
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-           ! domain from 15deg to 50deg and 1/2 period along 14deg
-           ! so 5/4 of half period with seasonal cycle
-           utau(ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) )
-           vtau(ji,jj) =   ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) )
-         END DO
-      END DO
-
-      ! module of wind stress and wind speed at T-point
-      zcoef = 1. / ( zrhoa * zcdrag ) 
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vect. opt.
-            ztx = utau(ji-1,jj  ) + utau(ji,jj) 
-            zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
-            zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
-            taum(ji,jj) = zmod
-            wndm(ji,jj) = SQRT( zmod * zcoef )
-         END DO
-      END DO
-      CALL lbc_lnk_multi( 'usrdef_sbc', taum(:,:), 'T', 1. , wndm(:,:), 'T', 1. )
-
-      ! ---------------------------------- !
-      !  control print at first time-step  !
-      ! ---------------------------------- !
-      IF( kt == nit000 .AND. lwp ) THEN 
-         WRITE(numout,*)
-         WRITE(numout,*)'usrdef_sbc_oce : analytical surface fluxes for GYRE configuration'               
-         WRITE(numout,*)'~~~~~~~~~~~ ' 
-         WRITE(numout,*)'           nyear      = ', nyear
-         WRITE(numout,*)'           nmonth     = ', nmonth
-         WRITE(numout,*)'           nday       = ', nday
-         WRITE(numout,*)'           nday_year  = ', nday_year
-         WRITE(numout,*)'           ztime      = ', ztime
-         WRITE(numout,*)'           ztimemax   = ', ztimemax
-         WRITE(numout,*)'           ztimemin   = ', ztimemin
-         WRITE(numout,*)'           ztimemax1  = ', ztimemax1
-         WRITE(numout,*)'           ztimemin1  = ', ztimemin1
-         WRITE(numout,*)'           ztimemax2  = ', ztimemax2
-         WRITE(numout,*)'           ztimemin2  = ', ztimemin2
-         WRITE(numout,*)'           zyear0     = ', zyear0
-         WRITE(numout,*)'           zmonth0    = ', zmonth0
-         WRITE(numout,*)'           zday0      = ', zday0
-         WRITE(numout,*)'           zday_year0 = ', zday_year0
-         WRITE(numout,*)'           zyydd      = ', zyydd
-         WRITE(numout,*)'           zemp_S     = ', zemp_S
-         WRITE(numout,*)'           zemp_N     = ', zemp_N
-         WRITE(numout,*)'           zemp_sais  = ', zemp_sais
-         WRITE(numout,*)'           zTstar     = ', zTstar
-         WRITE(numout,*)'           zsumemp    = ', zsumemp
-         WRITE(numout,*)'           zsurf      = ', zsurf
-         WRITE(numout,*)'           ztrp       = ', ztrp
-         WRITE(numout,*)'           zconv      = ', zconv
-         WRITE(numout,*)'           ndastp     = ', ndastp
-         WRITE(numout,*)'           adatrj     = ', adatrj
+         utau(:,:) = 0._wp 
+         vtau(:,:) = 0._wp 
+         taum(:,:) = 0._wp 
+         wndm(:,:) = 0._wp 
+         ! 
+         emp (:,:) = 0._wp 
+         sfx (:,:) = 0._wp 
+         qns (:,:) = 0._wp 
+         qsr (:,:) = 0._wp 
+         !
       ENDIF
       !
@@ -231,8 +107,6 @@
 
 
-   SUBROUTINE usrdef_sbc_ice_flx( kt, phs, phi )
+   SUBROUTINE usrdef_sbc_ice_flx( kt )
       INTEGER, INTENT(in) ::   kt   ! ocean time step
-      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
-      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
    END SUBROUTINE usrdef_sbc_ice_flx
 
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/par_oce.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/par_oce.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/par_oce.F90	(revision 15422)
@@ -64,4 +64,6 @@
    INTEGER, PUBLIC, PARAMETER ::   jp_tem = 1    !: indice for temperature
    INTEGER, PUBLIC, PARAMETER ::   jp_sal = 2    !: indice for salinity
+   INTEGER, PUBLIC, PARAMETER ::   jp_dep = 3    !: indice for depth  
+   INTEGER, PUBLIC, PARAMETER ::   jp_msk = 4    !: indice for mask
 
    !!----------------------------------------------------------------------
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/step.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/step.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/step.F90	(revision 15422)
@@ -206,4 +206,6 @@
                          CALL dia_ar5 ( kstp )        ! ar5 diag
       IF( ln_diaptr  )   CALL dia_ptr                 ! Poleward adv/ldf TRansports diagnostics
+      IF( lk_diaharm_fast )                           &  
+            &            CALL dia_harm_fast( kstp )   ! Tidal harmonic analysis - restart and faster version 
       IF( ln_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
                          CALL dia_wri ( kstp )        ! ocean model: outputs
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/step_oce.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/step_oce.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/step_oce.F90	(revision 15422)
@@ -79,4 +79,5 @@
    USE diahth          ! thermocline depth                (dia_hth routine)
    USE diahsb          ! heat, salt and volume budgets    (dia_hsb routine)
+   USE diaharm_fast    ! harmonic analysis of tides (harm_ana routine)
    USE diaharm
    USE diacfl
Index: /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/stpctl.F90
===================================================================
--- /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/stpctl.F90	(revision 15421)
+++ /NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/stpctl.F90	(revision 15422)
@@ -126,5 +126,5 @@
          zmax(1) = MAXVAL(  ABS( sshn(:,:) )  )                               ! ssh max
       ENDIF
-      zmax(2) = MAXVAL(  ABS( un(:,:,:) )  )                                  ! velocity max (zonal only)
+      zmax(2) = MAXVAL(  un(:,:,:)*un(:,:,:) + vn(:,:,:)*vn(:,:,:)  )         ! velocity max
       llmsk(:,:,:) = tmask(:,:,:) == 1._wp
       IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors...      
@@ -167,5 +167,5 @@
       !                                   !==  error handling  ==!
       IF(   zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m )
-         &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s)
+         &  zmax(2) >  100._wp .OR.   &                    ! too large velocity ( > 10 m/s)
          &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity
          &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 )
@@ -177,5 +177,5 @@
             IF( lwm .AND. kt /= nitend )   istatus = NF90_CLOSE(idrun)
             CALL mpp_maxloc( 'stpctl', ABS(sshn)        , ssmask(:,:)  , zzz, ih(1:2)  )   ;   ih(3) = 0
-            CALL mpp_maxloc( 'stpctl', ABS(un)          , umask (:,:,:), zzz, iu  )
+            CALL mpp_maxloc( 'stpctl', un*un + vn*vn    , umask (:,:,:), zzz, iu  )
             CALL mpp_minloc( 'stpctl', tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, is1 )
             CALL mpp_maxloc( 'stpctl', tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, is2 )
@@ -192,5 +192,5 @@
          ELSE
             ih(1:2)= MAXLOC( ABS( sshn(:,:)   )                              ) + (/ nimpp - 1, njmpp - 1    /)   ;   ih(3) = 0
-            iu(:)  = MAXLOC( ABS( un  (:,:,:) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /)
+            iu(:)  = MAXLOC( un(:,:,:)*un(:,:,:) + vn(:,:,:)*vn(:,:,:)       ) + (/ nimpp - 1, njmpp - 1, 0 /)
             is1(:) = MINLOC( tsn(:,:,:,jp_sal), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
             is2(:) = MAXLOC( tsn(:,:,:,jp_sal), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
@@ -200,5 +200,5 @@
          WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
          CALL wrt_line(ctmp2, kt, ' |ssh| max ',   zmax(1), ih , iareasum(1), iareamin(1), iareamax(1) ) 
-         CALL wrt_line(ctmp3, kt, ' |U|   max ',   zmax(2), iu , iareasum(2), iareamin(2), iareamax(2) ) 
+         CALL wrt_line(ctmp3, kt, ' Vel2  max ',   zmax(2), iu , iareasum(2), iareamin(2), iareamax(2) ) 
          CALL wrt_line(ctmp4, kt, ' Sal   min ', - zmax(3), is1, iareasum(3), iareamin(3), iareamax(3) ) 
          CALL wrt_line(ctmp5, kt, ' Sal   max ',   zmax(4), is2, iareasum(4), iareamin(4), iareamax(4) ) 
@@ -226,5 +226,5 @@
       ENDIF
       !
-9500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16)
+9500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' Vel2_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16)
       !
    END SUBROUTINE stp_ctl
