Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 8059)
@@ -71,4 +71,11 @@
       REAL, POINTER, DIMENSION(:,:)   ::  ht_s  !: now snow thickness
 #endif
+#if defined key_top
+      CHARACTER(LEN=20)                   :: cn_obc  !: type of boundary condition to apply
+      REAL(wp)                            :: rn_fac  !: multiplicative scaling factor
+      REAL(wp), POINTER, DIMENSION(:,:)   :: trc     !: now field of the tracer
+      LOGICAL                             :: dmp     !: obc damping term
+#endif
+
    END TYPE OBC_DATA
 
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 8059)
@@ -37,4 +37,8 @@
 #endif
    USE sbcapr
+#if defined key_top
+   USE par_trc
+   USE trc, ONLY: trn
+#endif
 
    IMPLICIT NONE
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 8059)
@@ -61,4 +61,8 @@
          CASE('zero')
             CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+         CASE('zerograd') 
+            CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )       
+         CASE('neumann') 
+            CALL bdy_dyn3d_nmn( idx_bdy(ib_bdy), ib_bdy )
          CASE('orlanski')
             CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. )
@@ -117,4 +121,54 @@
 
    END SUBROUTINE bdy_dyn3d_spe
+
+   SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt , ib_bdy ) 
+      !!---------------------------------------------------------------------- 
+      !!                  ***  SUBROUTINE bdy_dyn3d_zgrad  *** 
+      !! 
+      !! ** Purpose : - Enforce a zero gradient of normal velocity 
+      !! 
+      !!---------------------------------------------------------------------- 
+      INTEGER                     ::   kt 
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
+      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index 
+      !! 
+      INTEGER  ::   jb, jk         ! dummy loop indices 
+      INTEGER  ::   ii, ij, igrd   ! local integers 
+      REAL(wp) ::   zwgt           ! boundary weight 
+      INTEGER  ::   fu, fv 
+      !!---------------------------------------------------------------------- 
+      ! 
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zgrad') 
+      ! 
+      igrd = 2                      ! Copying tangential velocity into bdy points 
+      DO jb = 1, idx%nblenrim(igrd) 
+         DO jk = 1, jpkm1 
+            ii   = idx%nbi(jb,igrd) 
+            ij   = idx%nbj(jb,igrd) 
+            fu   = ABS( ABS (NINT( idx%flagu(jb,igrd) ) ) - 1 ) 
+            ua(ii,ij,jk) = ua(ii,ij,jk) * REAL( 1 - fu) + ( ua(ii,ij+fu,jk) * umask(ii,ij+fu,jk) & 
+                        &+ ua(ii,ij-fu,jk) * umask(ii,ij-fu,jk) ) * umask(ii,ij,jk) * REAL( fu ) 
+         END DO 
+      END DO 
+      ! 
+      igrd = 3                      ! Copying tangential velocity into bdy points 
+      DO jb = 1, idx%nblenrim(igrd) 
+         DO jk = 1, jpkm1 
+            ii   = idx%nbi(jb,igrd) 
+            ij   = idx%nbj(jb,igrd) 
+            fv   = ABS( ABS (NINT( idx%flagv(jb,igrd) ) ) - 1 ) 
+            va(ii,ij,jk) = va(ii,ij,jk) * REAL( 1 - fv ) + ( va(ii+fv,ij,jk) * vmask(ii+fv,ij,jk) & 
+                        &+ va(ii-fv,ij,jk) * vmask(ii-fv,ij,jk) ) * vmask(ii,ij,jk) * REAL( fv ) 
+         END DO 
+      END DO 
+      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ! Boundary points should be updated   
+      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )    
+      ! 
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
+ 
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zgrad') 
+ 
+   END SUBROUTINE bdy_dyn3d_zgrad 
 
    SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy )
@@ -303,4 +357,39 @@
    END SUBROUTINE bdy_dyn3d_dmp
 
+   SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy ) 
+      !!---------------------------------------------------------------------- 
+      !!                 ***  SUBROUTINE bdy_dyn3d_nmn  *** 
+      !!              
+      !!              - Apply Neumann condition to baroclinic velocities.  
+      !!              - Wrapper routine for bdy_nmn 
+      !!  
+      !! 
+      !!---------------------------------------------------------------------- 
+      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices 
+      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index 
+ 
+      INTEGER  ::   jb, igrd                               ! dummy loop indices 
+      !!---------------------------------------------------------------------- 
+ 
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_nmn') 
+      ! 
+      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities.  
+      ! 
+      igrd = 2      ! Neumann bc on u-velocity;  
+      !             
+      CALL bdy_nmn( idx, igrd, ua ) 
+ 
+      igrd = 3      ! Neumann bc on v-velocity 
+      !   
+      CALL bdy_nmn( idx, igrd, va ) 
+      ! 
+      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated 
+      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy ) 
+      ! 
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_nmn') 
+      ! 
+   END SUBROUTINE bdy_dyn3d_nmn 
+ 
+
 #else
    !!----------------------------------------------------------------------
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 8059)
@@ -213,4 +213,12 @@
              dta_bdy(ib_bdy)%ll_u3d = .true.
              dta_bdy(ib_bdy)%ll_v3d = .true.
+          CASE('neumann') 
+             IF(lwp) WRITE(numout,*) '      Neumann conditions' 
+             dta_bdy(ib_bdy)%ll_u3d = .false. 
+             dta_bdy(ib_bdy)%ll_v3d = .false. 
+          CASE('zerograd') 
+             IF(lwp) WRITE(numout,*) '      Zero gradient for baroclinic velocities' 
+             dta_bdy(ib_bdy)%ll_u3d = .false. 
+             dta_bdy(ib_bdy)%ll_v3d = .false.
           CASE('zero')
              IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
@@ -1087,5 +1095,7 @@
             DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
                nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
-               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
+               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) * 0.5 &
+                                            &  *(10./ FLOAT(nn_rimwidth(ib_bdy))) ) ! JGraham:modified for rim=15
+!               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
 !               idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.  ! quadratic
 !               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy))       ! linear
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90	(revision 8059)
@@ -26,4 +26,5 @@
    PUBLIC   bdy_orlanski_2d     ! routine called where?
    PUBLIC   bdy_orlanski_3d     ! routine called where?
+   PUBLIC   bdy_nmn     ! routine called where? 
 
    !!----------------------------------------------------------------------
@@ -354,4 +355,77 @@
    END SUBROUTINE bdy_orlanski_3d
 
+   SUBROUTINE bdy_nmn( idx, igrd, phia ) 
+      !!---------------------------------------------------------------------- 
+      !!                 ***  SUBROUTINE bdy_nmn  *** 
+      !!                     
+      !! ** Purpose : Duplicate the value at open boundaries, zero gradient. 
+      !!  
+      !!---------------------------------------------------------------------- 
+      INTEGER,                    INTENT(in)     ::   igrd     ! grid index 
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated) 
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
+      !!  
+      REAL(wp) ::   zcoef, zcoef1, zcoef2 
+      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field 
+      REAL(wp), POINTER, DIMENSION(:,:)        :: bdypmask      ! land/sea mask for field 
+      INTEGER  ::   ib, ik   ! dummy loop indices 
+      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses 
+      !!---------------------------------------------------------------------- 
+      ! 
+      IF( nn_timing == 1 ) CALL timing_start('bdy_nmn') 
+      ! 
+      SELECT CASE(igrd) 
+         CASE(1) 
+            pmask => tmask(:,:,:) 
+            bdypmask => bdytmask(:,:) 
+         CASE(2) 
+            pmask => umask(:,:,:) 
+            bdypmask => bdyumask(:,:) 
+         CASE(3) 
+            pmask => vmask(:,:,:) 
+            bdypmask => bdyvmask(:,:) 
+         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' ) 
+      END SELECT 
+      DO ib = 1, idx%nblenrim(igrd) 
+         ii = idx%nbi(ib,igrd) 
+         ij = idx%nbj(ib,igrd) 
+         DO ik = 1, jpkm1 
+            ! search the sense of the gradient 
+            zcoef1 = bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik) +  bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik) 
+            zcoef2 = bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik) +  bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik) 
+            IF ( nint(zcoef1+zcoef2) == 0) THEN 
+               ! corner **** we probably only want to set the tangentail component for the dynamics here 
+               zcoef = pmask(ii-1,ij,ik) + pmask(ii+1,ij,ik) +  pmask(ii,ij-1,ik) +  pmask(ii,ij+1,ik) 
+               IF (zcoef > .5_wp) THEN ! Only set none isolated points. 
+                 phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik) + & 
+                   &              phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik) + & 
+                   &              phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik) + & 
+                   &              phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik) 
+                 phia(ii,ij,ik) = ( phia(ii,ij,ik) / zcoef ) * pmask(ii,ij,ik) 
+               ELSE 
+                 phia(ii,ij,ik) = phia(ii,ij  ,ik) * pmask(ii,ij  ,ik) 
+               ENDIF 
+            ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN 
+               ! oblique corner **** we probably only want to set the normal component for the dynamics here 
+               zcoef = pmask(ii-1,ij,ik)*bdypmask(ii-1,ij  ) + pmask(ii+1,ij,ik)*bdypmask(ii+1,ij  ) + & 
+                   &   pmask(ii,ij-1,ik)*bdypmask(ii,ij -1 ) +  pmask(ii,ij+1,ik)*bdypmask(ii,ij+1  ) 
+               phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik)*bdypmask(ii-1,ij  ) + & 
+                   &            phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik)*bdypmask(ii+1,ij  )  + & 
+                   &            phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik)*bdypmask(ii,ij -1 ) + & 
+                   &            phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik)*bdypmask(ii,ij+1  ) 
+   
+               phia(ii,ij,ik) = ( phia(ii,ij,ik) / MAX(1._wp, zcoef) ) * pmask(ii,ij,ik) 
+            ELSE 
+               ip = nint(bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik) - bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik)) 
+               jp = nint(bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik) - bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik)) 
+               phia(ii,ij,ik) = phia(ii+ip,ij+jp,ik) * pmask(ii+ip,ij+jp,ik) 
+            ENDIF 
+         END DO 
+      END DO 
+      ! 
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_nmn') 
+      ! 
+   END SUBROUTINE bdy_nmn 
+
 
 #else
@@ -366,4 +440,7 @@
       WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt
    END SUBROUTINE bdy_orlanski_3d
+   SUBROUTINE bdy_nmn( idx, igrd, phia )      ! Empty routine 
+      WRITE(*,*) 'bdy_nmn: You should not have seen this print! error?', kt 
+   END SUBROUTINE bdy_nmn 
 #endif
 
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 8059)
@@ -50,4 +50,9 @@
       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V    (after nodal cor.)
    END TYPE TIDES_DATA
+   INTEGER, PUBLIC, PARAMETER                  ::   jptides_max = 15      !: Max number of tidal contituents
+      LOGICAL, PUBLIC                           ::   ln_harm_ana_store    !: =T Stores data for  harmonic Analysis
+      LOGICAL, PUBLIC                           ::   ln_harm_ana_compute     !: =T  Compute harmonic Analysis
+      LOGICAL, PUBLIC                           ::   ln_harmana_read         !: =T  Decide to do the analysis 
+                                                                             !from scratch or continue previous run
 
 !$AGRIF_DO_NOT_TREAT
@@ -90,5 +95,5 @@
       TYPE(MAP_POINTER), DIMENSION(jpbgrd)      ::   ibmap_ptr           !: array of pointers to nbmap
       !!
-      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
+      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj, ln_harm_ana_store, ln_harm_ana_compute, ln_harmana_read
       !!----------------------------------------------------------------------
 
@@ -102,4 +107,5 @@
 
       REWIND(numnam_cfg)
+      REWIND(numnam_ref)   ! slwa
 
       DO ib_bdy = 1, nb_bdy
@@ -125,4 +131,7 @@
             IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj
             IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo
+            IF(lwp) WRITE(numout,*) '             Use PCOMS harmonic ananalysis or not: ', ln_harm_ana_store
+            IF(lwp) WRITE(numout,*) '             Compute Final  harmonic ananalysis or not: ', ln_harm_ana_compute
+            IF(lwp) WRITE(numout,*) '             Read in previous days harmonic data or start afresh: ', ln_harmana_read
             IF(lwp) THEN 
                     WRITE(numout,*) '             Tidal components: ' 
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 8059)
@@ -91,6 +91,7 @@
       !! 
       REAL(wp) ::   zwgt           ! boundary weight
-      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
-      INTEGER  ::   ii, ij         ! 2D addresses
+      REAL(wp) ::   zcoef, zcoef1,zcoef2 
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
+      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses 
       !!----------------------------------------------------------------------
       !
@@ -160,6 +161,7 @@
       !! 
       REAL(wp) ::   zwgt           ! boundary weight
-      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
-      INTEGER  ::   ii, ij,zcoef, zcoef1,zcoef2, ip, jp   ! 2D addresses
+      REAL(wp) ::   zcoef, zcoef1,zcoef2 
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
+      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses 
       !!----------------------------------------------------------------------
       !
@@ -174,5 +176,5 @@
             zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  )
             zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1)
-            IF ( zcoef1+zcoef2 == 0) THEN
+            IF ( NINT(zcoef1+zcoef2) == 0) THEN 
                ! corner
                zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik)
@@ -181,13 +183,13 @@
                  &                    tsa(ii  ,ij-1,ik,jp_tem) * tmask(ii  ,ij-1,ik) + &
                  &                    tsa(ii  ,ij+1,ik,jp_tem) * tmask(ii  ,ij+1,ik)
-               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1, zcoef) ) * tmask(ii,ij,ik)
+               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik) 
                tsa(ii,ij,ik,jp_sal) = tsa(ii-1,ij  ,ik,jp_sal) * tmask(ii-1,ij  ,ik) + &
                  &                    tsa(ii+1,ij  ,ik,jp_sal) * tmask(ii+1,ij  ,ik) + &
                  &                    tsa(ii  ,ij-1,ik,jp_sal) * tmask(ii  ,ij-1,ik) + &
                  &                    tsa(ii  ,ij+1,ik,jp_sal) * tmask(ii  ,ij+1,ik)
-               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1, zcoef) ) * tmask(ii,ij,ik)
+               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik)
             ELSE
-               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
-               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
+               ip = NINT(bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )) 
+               jp = NINT(bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)) 
                tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik)
                tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik)
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 8059)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 8059)
@@ -0,0 +1,344 @@
+MODULE dia25h 
+   !!======================================================================
+   !!                       ***  MODULE  diaharm  ***
+   !! Harmonic analysis of tidal constituents 
+   !!======================================================================
+   !! History :  3.6  !  2014  (E O'Dea)  Original code
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain
+   USE diainsitutem, ONLY: rinsitu_t, theta2t
+   USE in_out_manager  ! I/O units
+   USE iom             ! I/0 library
+   USE wrk_nemo        ! working arrays
+#if defined key_zdftke 
+   USE zdftke, ONLY: en
+#endif
+   USE zdf_oce, ONLY: avt, avm
+#if defined key_zdfgls
+   USE zdfgls, ONLY: mxln
+   USE zdf_oce, ONLY: en
+#endif
+
+   IMPLICIT NONE
+   PRIVATE
+
+   LOGICAL , PUBLIC ::   ln_dia25h     !:  25h mean output
+   PUBLIC   dia_25h_init               ! routine called by nemogcm.F90
+   PUBLIC   dia_25h                    ! routine called by diawri.F90
+
+  !! * variables for calculating 25-hourly means
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_25h  , sn_25h, rinsitu_t_25h  
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_25h 
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_25h  , vn_25h  , wn_25h
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_25h , avm_25h
+#if defined key_zdfgls || key_zdftke
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_25h
+#endif
+#if defined key_zdfgls 
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   rmxln_25h
+#endif
+   INTEGER, SAVE :: cnt_25h     ! Counter for 25 hour means
+
+
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_25h_init 
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_25h_init  ***
+      !!     
+      !! ** Purpose: Initialization of 25h mean namelist 
+      !!        
+      !! ** Method : Read namelist
+      !!   History
+      !!   3.6  !  08-14  (E. O'Dea) Routine to initialize dia_25h
+      !!---------------------------------------------------------------------------
+      !!
+      INTEGER ::   ios                 ! Local integer output status for namelist read
+      INTEGER ::   ierror              ! Local integer for memory allocation
+      !
+      NAMELIST/nam_dia25h/ ln_dia25h
+      !!----------------------------------------------------------------------
+      !
+      REWIND ( numnam_ref )              ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics
+      READ   ( numnam_ref, nam_dia25h, IOSTAT=ios, ERR= 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist nam_dia25h in configuration namelist  25hour diagnostics
+      READ  ( numnam_cfg, nam_dia25h, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nam_dia25h )
+
+      IF(lwp) THEN                   ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_25h_init : Output 25 hour mean diagnostics'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) 'Namelist nam_dia25h : set 25h outputs '
+         WRITE(numout,*) 'Switch for 25h diagnostics (T) or not (F)  ln_dia25h  = ', ln_dia25h
+      ENDIF
+      IF( .NOT. ln_dia25h )   RETURN
+      ! ------------------- !
+      ! 1 - Allocate memory !
+      ! ------------------- !
+      ALLOCATE( tn_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate tn_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( sn_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate sn_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( rinsitu_t_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate rinsitu_t_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( un_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate un_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( vn_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate vn_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( wn_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate wn_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( avt_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate avt_25h' )   ;   RETURN
+      ENDIF
+      ALLOCATE( avm_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate avm_25h' )   ;   RETURN
+      ENDIF
+# if defined key_zdfgls || defined key_zdftke
+      ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate en_25h' )   ;   RETURN
+      ENDIF
+#endif
+# if defined key_zdfgls 
+      ALLOCATE( rmxln_25h(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate rmxln_25h' )   ;   RETURN
+      ENDIF
+#endif
+      ALLOCATE( sshn_25h(jpi,jpj), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_25h: unable to allocate sshn_25h' )   ;   RETURN
+      ENDIF
+      ! ------------------------- !
+      ! 2 - Assign Initial Values !
+      ! ------------------------- !
+      cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible) 
+      tn_25h(:,:,:) = tsb(:,:,:,jp_tem)
+      sn_25h(:,:,:) = tsb(:,:,:,jp_sal)
+      CALL theta2t
+      rinsitu_t_25h(:,:,:) = rinsitu_t(:,:,:)
+      sshn_25h(:,:) = sshb(:,:)
+      un_25h(:,:,:) = ub(:,:,:)
+      vn_25h(:,:,:) = vb(:,:,:)
+      wn_25h(:,:,:) = wn(:,:,:)
+      avt_25h(:,:,:) = avt(:,:,:)
+      avm_25h(:,:,:) = avm(:,:,:)
+# if defined key_zdfgls || defined key_zdftke
+         en_25h(:,:,:) = en(:,:,:)
+#endif
+# if defined key_zdfgls
+         rmxln_25h(:,:,:) = mxln(:,:,:)
+#endif
+#if defined key_lim3 || defined key_lim2
+         CALL ctl_stop('STOP', 'dia_25h not setup yet to do tidemean ice')
+#endif 
+
+      ! -------------------------- !
+      ! 3 - Return to dia_wri      !
+      ! -------------------------- !
+
+
+   END SUBROUTINE dia_25h_init
+
+
+   SUBROUTINE dia_25h( kt )  
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE dia_25h  ***
+      !!         
+      !!
+      !!--------------------------------------------------------------------
+      !!                   
+      !! ** Purpose :   Write diagnostics with M2/S2 tide removed
+      !!
+      !! ** Method  :   
+      !!      25hr mean outputs for shelf seas
+      !!
+      !! History :
+      !!   ?.0  !  07-04  (A. Hines) New routine, developed from dia_wri_foam
+      !!   3.4  !  02-13  (J. Siddorn) Routine taken from old dia_wri_foam
+      !!   3.6  !  08-14  (E. O'Dea) adapted for VN3.6
+      !!----------------------------------------------------------------------
+      !! * Modules used
+
+      IMPLICIT NONE
+
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+
+
+      !! * Local declarations
+      INTEGER ::   ji, jj, jk
+
+      LOGICAL ::   ll_print = .FALSE.    ! =T print and flush numout
+      REAL(wp)                         ::   zsto, zout, zmax, zjulian, zdt, zmdi  ! temporary reals
+      INTEGER                          ::   i_steps                               ! no of timesteps per hour
+      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                    ! temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                                  ! temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                                 ! temporary workspace
+      INTEGER                          ::   iyear0, nimonth0,iday0                ! start year,imonth,day
+
+      !!----------------------------------------------------------------------
+
+      ! 0. Initialisation
+      ! -----------------
+      ! Define frequency of summing to create 25 h mean
+      zdt = rdt
+      IF( nacc == 1 ) zdt = rdtmin
+
+      IF( MOD( 3600,INT(zdt) ) == 0 ) THEN
+         i_steps = 3600/INT(zdt)
+      ELSE
+         CALL ctl_stop('STOP', 'dia_wri_tide: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible')
+      ENDIF
+
+#if defined key_lim3 || defined key_lim2
+      CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ice')
+#endif
+
+      ! local variable for debugging
+      ll_print = ll_print .AND. lwp
+
+      ! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours
+      ! every day
+      IF( MOD( kt, i_steps ) == 0  .and. kt .ne. nn_it000 ) THEN
+
+         IF (lwp) THEN
+              WRITE(numout,*) 'dia_wri_tide : Summing instantaneous hourly diagnostics at timestep ',kt
+              WRITE(numout,*) '~~~~~~~~~~~~ '
+         ENDIF
+
+         tn_25h(:,:,:)        = tn_25h(:,:,:) + tsn(:,:,:,jp_tem)
+         sn_25h(:,:,:)        = sn_25h(:,:,:) + tsn(:,:,:,jp_sal)
+         CALL theta2t
+         rinsitu_t_25h(:,:,:)  = rinsitu_t_25h(:,:,:) + rinsitu_t(:,:,:)
+         sshn_25h(:,:)        = sshn_25h(:,:) + sshn (:,:)
+         un_25h(:,:,:)        = un_25h(:,:,:) + un(:,:,:)
+         vn_25h(:,:,:)        = vn_25h(:,:,:) + vn(:,:,:)
+         wn_25h(:,:,:)        = wn_25h(:,:,:) + wn(:,:,:)
+         avt_25h(:,:,:)       = avt_25h(:,:,:) + avt(:,:,:)
+         avm_25h(:,:,:)       = avm_25h(:,:,:) + avm(:,:,:)
+# if defined key_zdfgls || defined key_zdftke
+         en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:)
+#endif
+# if defined key_zdfgls
+         rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + mxln(:,:,:)
+#endif
+         cnt_25h = cnt_25h + 1
+
+         IF (lwp) THEN
+            WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far',cnt_25h
+         ENDIF
+
+      ENDIF ! MOD( kt, i_steps ) == 0
+
+         ! Write data for 25 hour mean output streams
+      IF( cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN
+
+            IF(lwp) THEN
+               WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt
+               WRITE(numout,*) '~~~~~~~~~~~~ '
+            ENDIF
+
+            tn_25h(:,:,:)        = tn_25h(:,:,:) / 25.0_wp
+            sn_25h(:,:,:)        = sn_25h(:,:,:) / 25.0_wp
+            rinsitu_t_25h(:,:,:)  = rinsitu_t_25h(:,:,:) / 25.0_wp
+            sshn_25h(:,:)        = sshn_25h(:,:) / 25.0_wp
+            un_25h(:,:,:)        = un_25h(:,:,:) / 25.0_wp
+            vn_25h(:,:,:)        = vn_25h(:,:,:) / 25.0_wp
+            wn_25h(:,:,:)        = wn_25h(:,:,:) / 25.0_wp
+            avt_25h(:,:,:)       = avt_25h(:,:,:) / 25.0_wp
+            avm_25h(:,:,:)       = avm_25h(:,:,:) / 25.0_wp
+# if defined key_zdfgls || defined key_zdftke
+            en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp
+#endif
+# if defined key_zdfgls
+            rmxln_25h(:,:,:)       = rmxln_25h(:,:,:) / 25.0_wp
+#endif
+
+            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output'
+            zmdi=1.e+20 !missing data indicator for masking
+            ! write tracers (instantaneous)
+            zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("temper25h", zw3d)   ! potential temperature
+            CALL theta2t                                                                    ! calculate insitu temp
+            zw3d(:,:,:) = rinsitu_t_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("tempis25h", zw3d)   ! in-situ temperature
+            zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put( "salin25h", zw3d  )   ! salinity
+            zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+            CALL iom_put( "ssh25h", zw2d )   ! sea surface 
+
+
+            ! Write velocities (instantaneous)
+            zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:))
+            CALL iom_put("vozocrtx25h", zw3d)    ! i-current
+            zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))
+            CALL iom_put("vomecrty25h", zw3d  )   ! j-current
+
+            zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("vomecrtz25h", zw3d )   ! k-current
+            zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("avt25h", zw3d )   ! diffusivity
+            zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("avm25h", zw3d)   ! viscosity
+#if defined key_zdftke || defined key_zdfgls 
+!           zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+!           CALL iom_put("tke25h", zw3d)   ! tke
+#endif
+#if defined key_zdfgls 
+!           zw3d(:,:,:) = rmxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+!           CALL iom_put( "mxln25h",zw3d)
+#endif
+
+            ! After the write reset the values to cnt=1 and sum values equal current value 
+            tn_25h(:,:,:) = tsn(:,:,:,jp_tem)
+            sn_25h(:,:,:) = tsn(:,:,:,jp_sal)
+            CALL theta2t
+            rinsitu_t_25h(:,:,:) = rinsitu_t(:,:,:)
+            sshn_25h(:,:) = sshn (:,:)
+            un_25h(:,:,:) = un(:,:,:)
+            vn_25h(:,:,:) = vn(:,:,:)
+            wn_25h(:,:,:) = wn(:,:,:)
+            avt_25h(:,:,:) = avt(:,:,:)
+            avm_25h(:,:,:) = avm(:,:,:)
+# if defined key_zdfgls || defined key_zdftke
+            en_25h(:,:,:) = en(:,:,:)
+#endif
+# if defined key_zdfgls
+            rmxln_25h(:,:,:) = mxln(:,:,:)
+#endif
+            cnt_25h = 1
+            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average',cnt_25h
+
+      ENDIF !  cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps * 24) == 0 .AND. kt .NE. nn_it000
+
+
+   END SUBROUTINE dia_25h 
+
+   !!======================================================================
+END MODULE dia25h
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diaharmana.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diaharmana.F90	(revision 8059)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diaharmana.F90	(revision 8059)
@@ -0,0 +1,667 @@
+MODULE harmonic_analysis
+   !!======================================================================
+   !!                       ***  MODULE  example  ***
+   !! Ocean physics:  On line harmonic analyser
+   !!                 
+   !!=====================================================================
+
+#if defined key_harm_ana
+   !!----------------------------------------------------------------------
+   !!   'key_harm_ana'  :                Calculate harmonic analysis
+   !!----------------------------------------------------------------------
+   !!   harm_ana        :
+   !!   harm_ana_init   :
+   !!----------------------------------------------------------------------
+
+   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_par         ! Unstructured boundary parameters
+   USE bdy_oce         ! ocean open boundary conditions
+   USE bdytides        ! tidal bdy forcing
+   USE daymod          ! calendar
+   USE tideini
+   USE restart
+   USE ioipsl, ONLY : ju2ymds    ! for calendar
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! *  Routine accessibility
+   PUBLIC harm_ana    ! routine called in step.F90 module
+
+   !! * Module variables
+   INTEGER, PARAMETER ::  nharm_max  =   jptides_max ! max number of harmonics to be analysed 
+   INTEGER, PARAMETER ::  nhm_max    =   2*jptides_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
+   REAL(wp) :: ccau,ccav,ssau,ssav
+   REAL(wp), PUBLIC, ALLOCATABLE,DIMENSION(:) :: anau, anav, anaf   ! nodel/phase corrections used by diaharmana
+!
+   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:) ::   &
+      bssh,bubar,bvbar        ! work array for ssh anaylsis
+   REAL(wp), ALLOCATABLE,SAVE, DIMENSION(:,:,:) ::   &
+      cosampz, cosampu,cosampv,      &
+      sinampz, sinampu,sinampv
+   REAL(WP), ALLOCATABLE,SAVE,DIMENSION(:,:) :: cc,a
+!
+   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:) ::   &
+      gout,hout        ! arrays for output
+   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:) ::   &
+      huout,hvout,guout,gvout        ! arrays for output
+   REAL(wp), PUBLIC ::   fjulday_startharm       !: Julian Day since start of harmonic analysis
+
+   !! * Substitutions
+
+   !!----------------------------------------------------------------------
+   !!  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 harm_ana( 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, jj, jk,  &        ! dummy loop arguments
+                   ih,i1,i2,iv,jgrid
+
+
+      !!--------------------------------------------------------------------
+
+
+
+      IF( kt == nit000  )   CALL harm_ana_init    ! Initialization (first time-step only)
+
+! this bit done every time step
+        !nharm=ncpt_anal
+        nharm=nb_harmo
+        nhm=2*nb_harmo+1
+       
+        c(1)=1.0
+
+        do ih=1,nb_harmo
+             c(2*ih)=    cos((fjulday-fjulday_startharm)*86400._wp*om_tide(ih)  )
+             c(2*ih+1)=  sin((fjulday-fjulday_startharm)*86400._wp*om_tide(ih)  )
+        enddo
+
+
+         do ji=1,jpi
+          do jj=1,jpj
+            do ih=1,nhm
+            bssh (ih,ji,jj)=bssh (ih,ji,jj)+c(ih)*sshn(ji,jj)
+            bubar(ih,ji,jj)=bubar(ih,ji,jj)+c(ih)*un_b(ji,jj)
+            bvbar(ih,ji,jj)=bvbar(ih,ji,jj)+c(ih)*vn_b(ji,jj)
+            enddo
+          enddo
+        enddo
+!  
+!        IF(lwp) WRITE(666,'(4E15.5)') adatrj*86400._wp, sshn(10,10),bssh(2,10,10),bssh(3,10,10)
+!
+         do i1=1,nhm
+           do i2=1,nhm
+             cc(i1,i2)=cc(i1,i2)+c(i1)*c(i2)
+            enddo
+         enddo
+!
+! 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,*) '~~~~~~~~~'
+      CALL harm_rst_write(kt)     ! Dump out data for a restarted run
+
+
+      if(ln_harm_ana_compute) then
+               WRITE(numout,*) "Computing harmonics at last step"
+!
+        cosampu=0.0_wp
+        sinampu=0.0_wp
+        cosampv=0.0_wp
+        sinampv=0.0_wp
+        cosampz=0.0_wp
+        sinampz=0.0_wp
+!
+     do jgrid=1,3 ! elevation, Ubar,Vbar
+        do ji=1,jpi
+            do jj=1,jpj
+               bt=1.0
+               do ih=1,nhm
+                   if(jgrid .eq. 1) then
+                     bzz(ih)=bssh(ih,ji,jj)
+                   endif
+                   if(jgrid .eq. 2) then
+                     bzz(ih)=bubar(ih,ji,jj)
+                   endif
+                   if(jgrid .eq. 3) then
+                     bzz(ih)=bvbar(ih,ji,jj)
+                   endif
+                   bt=bt*bzz(ih)
+               enddo
+
+               do i1=1,nhm
+                   do i2=1,nhm
+                       a(i1,i2)=cc(i1,i2)
+                   enddo
+               enddo
+!	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.0d0
+              if(bt.ne.0.) then
+                  if(lwp .and. ji == 10 .and. ji == 10) then
+              endif
+              call gelim(a,bzz,x,nhm)
+
+              do ih=1,nb_harmo
+                if(jgrid .eq. 1) then
+                   cosampz(ih,ji,jj)=x(ih*2)
+                   sinampz(ih,ji,jj)=x(ih*2+1)
+                endif
+                if(jgrid .eq. 2) then
+                   cosampu(ih,ji,jj)=x(ih*2)
+                   sinampu(ih,ji,jj)=x(ih*2+1)
+                endif
+                if(jgrid .eq. 3) then
+                   cosampv(ih,ji,jj)=x(ih*2)
+                   sinampv(ih,ji,jj)=x(ih*2+1)
+                endif
+              enddo
+                if(jgrid .eq. 1) then
+                  cosampz(0,ji,jj)=x(1)
+                  sinampz(0,ji,jj)=0.0_wp
+                endif
+                if(jgrid .eq. 2) then
+                  cosampu(0,ji,jj)=x(1)
+                  sinampu(0,ji,jj)=0.0_wp
+                endif
+                if(jgrid .eq. 3) then
+                  cosampv(0,ji,jj)=x(1)
+                  sinampv(0,ji,jj)=0.0_wp
+                endif
+              endif
+           enddo
+        enddo
+      enddo ! jgrid
+!
+!
+        CALL harm_ana_out     ! output analysis (last time step)
+     ELSE ! ln_harmana_compute 
+         WRITE(numout,*) "Skipping Computing harmonics at last step"
+     ENDIF
+       
+    ENDIF
+   END SUBROUTINE harm_ana
+
+   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, jj, jk, jit,ih   ! dummy loop indices
+
+      !!----------------------------------------------------------------------
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'harm_init : initialization of harmonic analysis of tides'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
+     ! DO ALLOCATIONS
+
+      ALLOCATE( bubar(nb_harmo*2+1,jpi,jpj) )
+      ALLOCATE( bvbar(nb_harmo*2+1,jpi,jpj) )
+      ALLOCATE( bssh (nb_harmo*2+1,jpi,jpj) )
+
+      ALLOCATE( cosampz(0:nb_harmo*2+1,jpi,jpj))
+      ALLOCATE( cosampu(0:nb_harmo*2+1,jpi,jpj))
+      ALLOCATE( cosampv(0:nb_harmo*2+1,jpi,jpj))
+
+      ALLOCATE( sinampz(0:nb_harmo*2+1,jpi,jpj))
+      ALLOCATE( sinampu(0:nb_harmo*2+1,jpi,jpj))
+      ALLOCATE( sinampv(0:nb_harmo*2+1,jpi,jpj))
+
+      ALLOCATE( cc(nb_harmo*2+1,nb_harmo*2+1) )
+
+      ALLOCATE( a(nb_harmo*2+1,nb_harmo*2+1) )
+
+      ALLOCATE( anau(nb_harmo) )
+      ALLOCATE( anav(nb_harmo) )
+      ALLOCATE( anaf(nb_harmo) )
+
+      ALLOCATE( bzz(nb_harmo*2+1) )
+      ALLOCATE( x(nb_harmo*2+1) )
+      ALLOCATE( c(nb_harmo*2+1) )
+
+      ALLOCATE( gout(jpi,jpj))
+      ALLOCATE( hout(jpi,jpj))
+
+      ALLOCATE( guout(jpi,jpj))
+      ALLOCATE( huout(jpi,jpj))
+
+      ALLOCATE( gvout(jpi,jpj))
+      ALLOCATE( hvout(jpi,jpj))
+
+        ! END ALLOCATE 
+
+                 do ih=1,nb_harmo
+                  om_tide(ih)= omega_tide(ih) 
+                 enddo
+       if(ln_harmana_read) then
+               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
+               WRITE(numout,*) "Starting harmonic analysis from Fresh "
+         bssh(:,:,:)  = 0.0_wp
+         bubar(:,:,:) = 0.0_wp
+         bvbar(:,:,:) = 0.0_wp
+
+         cc=0.0_wp
+
+         anau(:) = 0.0_wp
+         anav(:) = 0.0_wp
+         anaf(:) = 0.0_wp
+         bzz(:) = 0.0_wp
+         x(:) = 0.0_wp
+         c(:) = 0.0_wp
+
+         DO ih = 1, nb_harmo
+            anau(ih) = utide(ih)
+            anav(ih) = v0tide(ih)
+            anaf(ih) = ftide(ih)
+
+         END DO
+         fjulday_startharm=fjulday !Set this at very start and store
+      endif
+
+   END SUBROUTINE harm_ana_init
+!
+   SUBROUTINE gelim (a,b,x,n)
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE harm_ana  ***
+      !!
+      !! ** Purpose :   Guassian elimination
+      !!
+      !! ** Method  :   
+      !!
+      !! ** 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_harmo*2+1),a(nb_harmo*2+1,nb_harmo*2+1)
+   REAL(WP) :: x(nb_harmo*2+1)
+        INTEGER  :: row,col,prow,pivrow,rrow,ntemp
+   REAL(WP) ::  atemp
+   REAL(WP) ::  pivot
+   REAL(WP) ::  m
+   REAL(WP) ::  tempsol(nb_harmo*2+1)
+
+
+        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,jit, ih, jjk   ! dummy loop indices
+      INTEGER :: nh_T
+      INTEGER :: nid_harm
+      CHARACTER (len=40) ::           &
+         clhstnamt, clop1, clop2          ! temporary names 
+      CHARACTER (len=40) ::           &
+         clhstnamu, clhstnamv   ! temporary names 
+      REAL(wp) ::   &
+         zsto1, zsto2, zout, zmax,            &  ! temporary scalars
+         zjulian, zdt, zmdi  
+         
+        do jgrid=1,3
+          do ih=1,nb_harmo
+            hout = 0.0
+            gout = 0.0
+            do jj=1,nlcj
+                do ji=1,nlci
+                    if(jgrid .eq. 1) then ! SSH 
+                       cca=cosampz(ih,ji,jj)
+                       ssa=sinampz(ih,ji,jj)
+                    endif
+                    if(jgrid .eq. 2) then ! UBAR
+                       cca=cosampu(ih,ji,jj)
+                       ssa=sinampu(ih,ji,jj)
+                    endif
+                    if(jgrid .eq. 3) then ! VBAR
+                       cca=cosampv(ih,ji,jj)
+                       ssa=sinampv(ih,ji,jj)
+                    endif
+
+                    hout(ji,jj)=sqrt(cca**2+ssa**2)
+        
+                    if (cca.ne.0.0) then
+                         gout(ji,jj)=(180.0/rpi)*atan(ssa/cca)
+                    else
+                        if (ssa.ne.0.0) then
+                            if (ssa.gt. 0.0) then
+                                 gout(ji,jj)=90.0
+                            else
+                                 gout(ji,jj)=270.0
+                            endif
+                        else
+                            gout(ji,jj)=0.0
+                        endif
+                    endif
+
+                    if (gout(ji,jj).lt.0) then
+                        gout(ji,jj)=gout(ji,jj)+180.0
+                    endif
+                    if (ssa.lt.0) then
+                        gout(ji,jj)=gout(ji,jj)+180.0
+                    endif
+                    
+                    if (hout(ji,jj).ne.0) then
+                        hout(ji,jj)=hout(ji,jj)/anaf(ih)
+                    endif
+                    if (gout(ji,jj).ne.0) then
+                                                !Convert Rad to degree and take
+                                                !modulus
+                         gout(ji,jj)=gout(ji,jj)+MOD( (anau(ih)+anav(ih))/rad , 360.0)
+                        if (gout(ji,jj).gt.360.0) then
+                            gout(ji,jj)=gout(ji,jj)-360.0
+                        else if (gout(ji,jj).lt.0) then
+                            gout(ji,jj)=gout(ji,jj)+360.0
+                        endif
+                    endif
+                enddo
+            enddo
+!NETCDF OUTPUT
+            if(jgrid==1) then
+             WRITE(numout,*) TRIM(Wave(ntide(ih))%cname_tide)//'amp'
+             WRITE(numout,*) TRIM(Wave(ntide(ih))%cname_tide)//'phase'
+             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp', hout(:,:) )
+             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase', gout(:,:) )
+            endif
+            if(jgrid==2) then
+             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp_u', hout(:,:) )
+             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase_u', gout(:,:) )
+            endif
+            if(jgrid==3) then
+             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp_v', hout(:,:) )
+             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase_v', gout(:,:) )
+            endif
+          enddo
+        enddo
+ 
+!
+   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             ::   iyear, imonth, iday,ih
+      REAL (wp)           ::   zsec
+      !!
+      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
+
+      CALL ju2ymds( fjulday , iyear, imonth, iday, zsec)
+
+      WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
+      ! create the file
+      clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana"
+      clpath = TRIM(cn_ocerst_outdir)
+      IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
+      WRITE(numout,*) 'Open tidal harmonics restart file for writing: ',TRIM(clpath)//clname
+
+
+
+
+!!      clhstnam=TRIM(cexper)//'.restart_harm_ana'
+      ! Open a file to write the data into
+         write (clfinal,'(a,''_'',i4.4)') trim(clpath)//trim(clname),nproc
+         open(66,file=TRIM(clfinal))
+         write(66,'(1e20.9)') cc
+         !These Three which contain the most data can be moved to a regular
+         !restart file
+
+            CALL iom_rstput( kt, nitrst, numrow, 'Mean_bssh'     , bssh(1,:,:)       )
+            CALL iom_rstput( kt, nitrst, numrow, 'Mean_bubar'    , bubar(1,:,:)       )
+            CALL iom_rstput( kt, nitrst, numrow, 'Mean_bvbar'    , bvbar(1,:,:)       )
+         do ih=1,nb_harmo
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_cos'    , bssh (ih*2  , : , : )     )
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_sin'    , bssh (ih*2+1, : , : )     )
+
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_cos'   , bubar(ih*2  , : , : )    )
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_sin'   , bubar(ih*2+1, : , : )    )
+
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_cos'   , bvbar(ih*2  , : , : )    )
+            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_sin'   , bvbar(ih*2+1, : , : )    )
+         enddo
+
+         write(66,'(1e20.9)') anau
+         write(66,'(1e20.9)') anav
+         write(66,'(1e20.9)') anaf
+         write(66,'(1e25.20)') fjulday_startharm
+         close(66)
+ 
+   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=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             ::   ih
+
+      ! create the file
+      clname = "restart_harm_ana"
+      clpath = TRIM(cn_ocerst_indir)
+      IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
+      WRITE(numout,*) 'Open tidal harmonics restart file for reading: ',TRIM(clpath)//clname
+
+
+
+
+!!      clhstnam=TRIM(cexper)//'.restart_harm_ana'
+      ! Open a file to read the data ifrom
+         write (clfinal,'(a,''_'',i4.4)') trim(clpath)//trim(clname),nproc
+         open(66,file=TRIM(clfinal),status='old')
+         read(66,'(1e20.9)') cc
+!2D regular fields can be read from normal restart this saves space and handy to
+!view in netcdf format also.
+         CALL iom_get( numror,jpdom_autoglo, 'Mean_bssh'     , bssh(1,:,:)       )
+         CALL iom_get( numror,jpdom_autoglo, 'Mean_bubar'    , bubar(1,:,:)       )
+         CALL iom_get( numror,jpdom_autoglo, 'Mean_bvbar'    , bvbar(1,:,:)       )
+          
+         do ih=1,nb_harmo
+
+            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_cos'    , bssh (ih*2  , : , : )     )
+            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_sin'    , bssh (ih*2+1, : , : )     )
+
+            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_cos'   , bubar(ih*2  , : , : )    )
+            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_sin'   , bubar(ih*2+1, : , : )    )
+
+            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_cos'   , bvbar(ih*2  , : , : )    )
+            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_sin'   , bvbar(ih*2+1, : , : )    )
+         enddo
+         read(66,'(1e20.9)') anau
+         read(66,'(1e20.9)') anav
+         read(66,'(1e20.9)') anaf
+         read(66,'(1e25.20)') fjulday_startharm
+      WRITE(numout,*) 'Checking anaf is correct'
+      WRITE(numout,*) anaf
+      WRITE(numout,*) '---end of anaf checking----'
+
+         close(66)
+
+   END SUBROUTINE harm_rst_read
+
+   !!======================================================================
+#else
+!!---------------------------------------------------------------------------------
+!!   Dummy module                                   NO harmonic Analysis
+!!---------------------------------------------------------------------------------
+        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 )
+           END SUBROUTINE harm_ana_init
+           SUBROUTINE gelim (a,b,x,n)
+           END SUBROUTINE gelim (a,b,x,n)
+           
+#endif
+
+END MODULE harmonic_analysis
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90	(revision 8059)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90	(revision 8059)
@@ -0,0 +1,124 @@
+MODULE diainsitutem
+   !!======================================================================
+   !!                       ***  MODULE  diaharm  ***
+   !! Harmonic analysis of tidal constituents
+   !!======================================================================
+   !! Initial version: D. Acreman June 2006 Original Code
+   !! History :  3.6  !  2014  (E O'Dea) for 3.6 25hour mean output
+   !!----------------------------------------------------------------------
+
+  USE dom_oce  ! ocean space and time domain
+  USE oce,            ONLY: tsn
+  USE par_oce,        ONLY: jpi, jpj, jpk, jpkm1
+  USE lbclnk,         ONLY: lbc_lnk
+  USE lib_mpp         ! MPP library
+
+  IMPLICIT NONE
+  PRIVATE
+
+  PUBLIC theta2t
+
+  !! * Accessibility
+  PUBLIC insitu_tem_alloc          ! routines called by step.F90
+
+  REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: rinsitu_t
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+
+CONTAINS
+
+   REAL FUNCTION insitu_tem_alloc()
+      !!----------------------------------------------------------------------
+      INTEGER, DIMENSION(2) :: ierr
+      !!----------------------------------------------------------------------
+      !
+      ierr = 0
+      !
+      ALLOCATE( rinsitu_t(jpi,jpj,jpk), STAT=ierr(1) )
+         !
+      insitu_tem_alloc = MAXVAL(ierr)
+      IF( lk_mpp )   CALL mpp_sum( insitu_tem_alloc )
+      !
+  END FUNCTION insitu_tem_alloc
+
+!-----------------------------------------------------------------------------------
+!
+! Calculate the insitu temperature by integrating the adiabatic lapse rate from zero 
+! to pressure at depth of tracer level. Based on UM subroutine POTTEM and function ATG
+!
+! Initial version: D. Acreman June 2006
+!
+!-----------------------------------------------------------------------------------
+
+  SUBROUTINE theta2t()
+
+    INTEGER, PARAMETER :: inum_steps=10                ! number of steps in integration
+    INTEGER            :: jstep                        ! iteration counter
+    INTEGER            :: ji, jj, jk                  ! loop indices
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zP            ! pressure (decibars)
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zT            ! temperature at pressure p
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zTB           ! temperature at p-dp
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zTA           ! temperature at p+dp
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zDP           ! pressure step
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zSS           ! salinity(PSU) - 35.0
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zLAPSE        ! adiabatic lapse rate
+
+!CDIR IEXPAND (ATG)
+
+! Set the pressure interval for the integration. The integration is carried out from 
+! zero (pressure at the surface) to the pressure at the depth of the tracer point. The
+! pressure in decibars is represented by the depth in metres. Pressures are "Oceanographic"
+! pressures equal to absolute pressure minus one atmosphere
+     zDP(:,:,:) = 0.0
+     DO jk = 1, jpkm1
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+              ! These loops expanded for case where fsdept may be 1D
+              zDP(ji,jj,jk) = fsdept(ji,jj,jk) / real(inum_steps)
+           END DO
+        END DO
+     END DO
+
+! Salinity at each point
+     zSS(:,:,:)      = tsn(:,:,:,jp_sal) - 35.0 
+
+! Set initial values of temperature and pressure. zT is the temperature at pressure zP, 
+! zTB is the temperature at pressure zP-zdP and zTA is the temperature at pressure zP+zdP
+     zT(:,:,:)  = tsn(:,:,:,jp_tem)
+     zP(:,:,:)  = 0.0         ! Pressure at surface
+     CALL ATG(zP, zT, zSS, zLAPSE)
+     zTB(:,:,:) = zT(:,:,:)  - zLAPSE(:,:,:) * zDP(:,:,:)
+
+     interation: DO jstep=1, inum_steps
+        ! Calculate lapse rate (dT/dP) and hence TA
+        CALL ATG(zP, zT, zSS, zLAPSE)
+        zTA(:,:,:) = zTB(:,:,:) + 2.0 * zLAPSE(:,:,:) * zDP(:,:,:)
+        ! Have calculated TB, T and TA for this pressure, now advance solution by dP
+        zP(:,:,:)  = zP(:,:,:)  + zDP
+        zTB(:,:,:) = zT(:,:,:)
+        zT(:,:,:)  = zTA(:,:,:)
+
+     END DO interation
+
+     rinsitu_t(:,:,:) = zT(:,:,:) * tmask(:,:,:)
+     CALL lbc_lnk( rinsitu_t,  'T', 1.0)
+
+   END SUBROUTINE theta2t
+
+   SUBROUTINE ATG(P,T,DS,zLAPSE)
+
+     REAL, INTENT(IN   ) :: P(jpi,jpj,jpk)        ! PRESSURE (DECIBARS)
+     REAL, INTENT(IN   ) :: T(jpi,jpj,jpk)        ! TEMPERATURE (DEG C)
+     REAL, INTENT(IN   ) :: DS(jpi,jpj,jpk)       ! SALINITY (PSU) -35.0
+     REAL, INTENT(  OUT) :: zLAPSE(jpi,jpj,jpk)   ! LAPSE RATE
+
+       zLAPSE = ((( -2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P &
+           + ((2.7759E-12*T-1.1351E-10)*DS+((-5.4481E-14*T    &
+           + 8.733E-12)*T-6.7795E-10)*T+1.8741E-8))*P         &
+           + (-4.2393E-8*T+1.8932E-6)*DS                      &
+           + ((6.6228E-10*T-6.836E-8)*T+8.5258E-6)*T+3.5803E-5 
+
+   END SUBROUTINE ATG
+
+END MODULE diainsitutem
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90	(revision 8059)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90	(revision 8059)
@@ -0,0 +1,169 @@
+MODULE diatmb 
+   !!======================================================================
+   !!                       ***  MODULE  diaharm  ***
+   !! Harmonic analysis of tidal constituents 
+   !!======================================================================
+   !! History :  3.6  !  2014  (E O'Dea)  Original code
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain
+   USE in_out_manager  ! I/O units
+   USE iom             ! I/0 library
+   USE wrk_nemo        ! working arrays
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   LOGICAL , PUBLIC ::   ln_diatmb     !: Top Middle and Bottom output 
+   PUBLIC   dia_tmb_init            ! routine called by nemogcm.F90
+   PUBLIC   dia_tmb                 ! routine called by diawri.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_tmb_init 
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_tmb_init  ***
+      !!     
+      !! ** Purpose: Initialization of tmb namelist 
+      !!        
+      !! ** Method : Read namelist
+      !!   History
+      !!   3.6  !  08-14  (E. O'Dea) Routine to initialize dia_tmb
+      !!---------------------------------------------------------------------------
+      !!
+      INTEGER ::   ios                 ! Local integer output status for namelist read
+      !
+      NAMELIST/nam_diatmb/ ln_diatmb
+      !!----------------------------------------------------------------------
+      !
+      REWIND ( numnam_ref )              ! Read Namelist nam_diatmb in reference namelist : TMB diagnostics
+      READ   ( numnam_ref, nam_diatmb, IOSTAT=ios, ERR= 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diatmb in reference namelist', lwp )
+ 
+      REWIND( numnam_cfg )              ! Namelist nam_diatmb in configuration namelist  TMB diagnostics
+      READ  ( numnam_cfg, nam_diatmb, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diatmb in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nam_diatmb )
+
+      IF(lwp) THEN                   ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_tmb_init : Output Top, Middle, Bottom Diagnostics'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) 'Namelist nam_diatmb : set tmb outputs '
+         WRITE(numout,*) 'Switch for TMB diagnostics (T) or not (F)  ln_diatmb  = ', ln_diatmb
+      ENDIF
+
+   END SUBROUTINE dia_tmb_init
+
+   SUBROUTINE dia_calctmb( pinfield,pouttmb )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_tmb  ***
+      !!                   
+      !! ** Purpose :    Find the Top, Mid and Bottom fields of water Column
+      !!
+      !! ** Method  :   
+      !!      use mbathy to find surface, mid and bottom of model levels
+      !!
+      !! History :
+      !!   3.6  !  08-14  (E. O'Dea) Routine based on dia_wri_foam
+      !!----------------------------------------------------------------------
+      !! * Modules used
+
+      ! Routine to map 3d field to top, middle, bottom
+      IMPLICIT NONE
+
+
+      ! Routine arguments
+      REAL(wp), DIMENSION(jpi, jpj, jpk), INTENT(IN   ) :: pinfield    ! Input 3d field and mask
+      REAL(wp), DIMENSION(jpi, jpj, 3  ), INTENT(  OUT) :: pouttmb     ! Output top, middle, bottom
+
+
+
+      ! Local variables
+      INTEGER :: ji,jj,jk  ! Dummy loop indices
+
+      ! Local Real
+      REAL(wp)                         ::   zmdi  !  set masked values
+
+      zmdi=1.e+20 !missing data indicator for masking
+
+      ! Calculate top
+      pouttmb(:,:,1) = pinfield(:,:,1)*tmask(:,:,1)  + zmdi*(1.0-tmask(:,:,1))
+
+      ! Calculate middle
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            jk              = max(1,mbathy(ji,jj)/2)
+            pouttmb(ji,jj,2) = pinfield(ji,jj,jk)*tmask(ji,jj,jk)  + zmdi*(1.0-tmask(ji,jj,jk))
+         END DO
+      END DO
+
+      ! Calculate bottom
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            jk              = max(1,mbathy(ji,jj) )
+            pouttmb(ji,jj,3) = pinfield(ji,jj,jk)*tmask(ji,jj,jk)  + zmdi*(1.0-tmask(ji,jj,jk))
+         END DO
+      END DO
+
+   END SUBROUTINE dia_calctmb
+
+
+
+   SUBROUTINE dia_tmb
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE dia_tmb  ***
+      !! ** Purpose :   Write diagnostics for Top, Mid and Bottom of water Column
+      !!
+      !! ** Method  :   
+      !!      use mbathy to find surface, mid and bottom of model levels
+      !!      calls calctmb to retrieve TMB values before sending to iom_put
+      !!
+      !! History :
+      !!   3.6  !  08-14  (E. O'Dea) 
+      !!         
+      !!--------------------------------------------------------------------
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmb    ! temporary workspace 
+      REAL(wp)                         ::   zmdi      ! set masked values
+
+      zmdi=1.e+20 !missing data indicator for maskin
+
+      IF (ln_diatmb) THEN
+         CALL wrk_alloc( jpi , jpj, 3 , zwtmb )
+         CALL dia_calctmb(  tsn(:,:,:,jp_tem),zwtmb )
+         !ssh already output but here we output it masked
+         CALL iom_put( "sshnmasked" , sshn(:,:)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) )   ! tmb Temperature
+         CALL iom_put( "top_temp" , zwtmb(:,:,1) )    ! tmb Temperature
+         CALL iom_put( "mid_temp" , zwtmb(:,:,2) )    ! tmb Temperature
+         CALL iom_put( "bot_temp" , zwtmb(:,:,3) )    ! tmb Temperature
+!         CALL iom_put( "sotrefml" , hmld_tref(:,:) )    ! "T criterion Mixed Layer Depth
+
+         CALL dia_calctmb(  tsn(:,:,:,jp_sal),zwtmb )
+         CALL iom_put( "top_sal" , zwtmb(:,:,1) )    ! tmb Salinity 
+         CALL iom_put( "mid_sal" , zwtmb(:,:,2) )    ! tmb Salinity
+         CALL iom_put( "bot_sal" , zwtmb(:,:,3) )    ! tmb Salinity
+
+         CALL dia_calctmb(  un(:,:,:),zwtmb )
+         CALL iom_put( "top_u" , zwtmb(:,:,1) )    ! tmb  U Velocity
+         CALL iom_put( "mid_u" , zwtmb(:,:,2) )    ! tmb  U Velocity
+         CALL iom_put( "bot_u" , zwtmb(:,:,3) )    ! tmb  U Velocity
+!Called in  dynspg_ts.F90        CALL iom_put( "baro_u" , un_b )    ! Barotropic  U Velocity
+
+         CALL dia_calctmb(  vn(:,:,:),zwtmb )
+         CALL iom_put( "top_v" , zwtmb(:,:,1) )    ! tmb  V Velocity
+         CALL iom_put( "mid_v" , zwtmb(:,:,2) )    ! tmb  V Velocity
+         CALL iom_put( "bot_v" , zwtmb(:,:,3) )    ! tmb  V Velocity
+!Called in  dynspg_ts.F90       CALL iom_put( "baro_v" , vn_b )    ! Barotropic  V Velocity
+      ELSE
+         CALL ctl_warn('dia_tmb: tmb diagnostic is set to false you should not have seen this')
+      ENDIF
+
+   END SUBROUTINE dia_tmb
+   !!======================================================================
+END MODULE diatmb
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 8059)
@@ -44,7 +44,11 @@
    USE in_out_manager  ! I/O manager
    USE diadimg         ! dimg direct access file format output
+   USE diatmb          ! Top,middle,bottom output
+   USE dia25h          ! 25h Mean output
    USE iom
    USE ioipsl
    USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities     
+   USE eosbn2         ! equation of state                (eos_bn2 routine)
+
 
 #if defined key_lim2
@@ -132,4 +136,5 @@
       REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zrhd , zrhop  ! 3D workspace
       !!----------------------------------------------------------------------
       ! 
@@ -138,4 +143,5 @@
       CALL wrk_alloc( jpi , jpj      , z2d )
       CALL wrk_alloc( jpi , jpj, jpk , z3d )
+      CALL wrk_alloc( jpi , jpj, jpk , zrhd , zrhop )
       !
       ! Output the initial state and forcings
@@ -376,7 +382,24 @@
          CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction
       ENDIF
+
+      IF( iom_use("rhop") ) THEN
+         CALL eos( tsn, zrhd, zrhop, fsdept_n(:,:,:) )       ! now in situ and potential density
+         zrhop(:,:,jpk) = 0._wp
+         CALL iom_put( 'rhop', zrhop )
+      ENDIF
+
       !
       CALL wrk_dealloc( jpi , jpj      , z2d )
       CALL wrk_dealloc( jpi , jpj, jpk , z3d )
+      CALL wrk_dealloc( jpi , jpj, jpk , zrhd , zrhop )
+      !
+      ! If we want tmb values 
+
+      IF (ln_diatmb) THEN
+         CALL dia_tmb
+      ENDIF
+      IF (ln_dia25h) THEN
+         CALL dia_25h( kt )
+      ENDIF
       !
       IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 8059)
@@ -136,5 +136,5 @@
       USE ioipsl
       NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
-         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
+         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstdate, ln_rstart , nn_rstctl,   &
          &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
          &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler
@@ -174,4 +174,5 @@
          WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
          WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
+         WRITE(numout,*) '      datestamping of restarts        ln_rstdate  = ', ln_rstdate
          WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
          WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 8059)
@@ -31,4 +31,5 @@
    USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
+   USE iom    ! slwa
 
    IMPLICIT NONE
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 8059)
@@ -209,5 +209,5 @@
       ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) 
       !
-      IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==!
+      IF( ln_sco .and. 1==0) THEN                   !==   s- or mixed s-zps-coordinate   ==!
          !
          CALL wrk_alloc( jpk, ztp, zsp )
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 8059)
@@ -24,4 +24,7 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+#if defined key_bdy 
+   USE bdy_oce         ! ocean open boundary conditions 
+#endif 
 
    IMPLICIT NONE
@@ -78,4 +81,9 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
       REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
+#if defined key_bdy 
+      INTEGER  ::   jb                 ! dummy loop indices 
+      INTEGER  ::   ii, ij, igrd, ib_bdy   ! local integers 
+      INTEGER  ::   fu, fv 
+#endif 
       !!----------------------------------------------------------------------
       !
@@ -97,5 +105,33 @@
       
       zhke(:,:,jpk) = 0._wp
-      
+
+#if defined key_bdy 
+      ! Maria Luneva & Fred Wobus: July-2016 
+      ! compensate for lack of turbulent kinetic energy on liquid bdy points 
+      DO ib_bdy = 1, nb_bdy 
+         IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
+            igrd = 2           ! Copying normal velocity into points outside bdy 
+            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
+               DO jk = 1, jpkm1 
+                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
+                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
+                  fu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) 
+                  un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk) 
+               END DO 
+            END DO 
+            ! 
+            igrd = 3           ! Copying normal velocity into points outside bdy 
+            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
+               DO jk = 1, jpkm1 
+                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
+                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
+                  fv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) 
+                  vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk) 
+               END DO 
+            END DO 
+         ENDIF 
+      ENDDO  
+#endif       
+
       SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
       !
@@ -134,4 +170,11 @@
       END SELECT
       !
+
+#if defined key_bdy 
+      ! restore velocity masks at points outside boundary 
+      un(:,:,:) = un(:,:,:) * umask(:,:,:) 
+      vn(:,:,:) = vn(:,:,:) * vmask(:,:,:) 
+#endif  
+
       DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
          DO jj = 2, jpjm1
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 8059)
@@ -41,4 +41,5 @@
    USE timing          ! Timing    
    USE sbcapr          ! surface boundary condition: atmospheric pressure
+   USE diatmb          ! Top,middle,bottom output
    USE dynadv, ONLY: ln_dynadv_vec
 #if defined key_agrif
@@ -144,4 +145,5 @@
       INTEGER  ::   ji, jj, jk, jn   		! dummy loop indices
       INTEGER  ::   ikbu, ikbv, noffset    	! local integers
+      REAL(wp) ::   zmdi
       REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf  	! local scalars
       REAL(wp) ::   zx1, zy1, zx2, zy2         !   -      -
@@ -169,4 +171,5 @@
       CALL wrk_alloc( jpi, jpj, zhf )
       !
+      zmdi=1.e+20                               !  missing data indicator for masking
       !                                         !* Local constant initialization
       z1_12 = 1._wp / 12._wp 
@@ -926,4 +929,8 @@
       CALL wrk_dealloc( jpi, jpj, zhf )
       !
+      IF ( ln_diatmb ) THEN
+         CALL iom_put( "baro_u" , un_b*umask(:,:,1)+zmdi*(1-umask(:,:,1 ) ) )  ! Barotropic  U Velocity
+         CALL iom_put( "baro_v" , vn_b*vmask(:,:,1)+zmdi*(1-vmask(:,:,1 ) ) )  ! Barotropic  V Velocity
+      ENDIF
       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
       !
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 8059)
@@ -18,6 +18,8 @@
    !!----------------------------------------------------------------------
    USE par_oce        ! NEMO parameters
+   USE phycst         ! for rday
    USE dom_oce        ! NEMO domain
    USE in_out_manager ! NEMO IO routines
+   USE ioipsl, ONLY : ju2ymds    ! for calendar
    USE lib_mpp        ! NEMO MPI library, lk_mpp in particular
    USE netcdf         ! netcdf routines for IO
@@ -231,6 +233,10 @@
       INTEGER ::   jn   ! dummy loop index
       INTEGER ::   ix_dim, iy_dim, ik_dim, in_dim
-      CHARACTER(len=256)     :: cl_path
-      CHARACTER(len=256)     :: cl_filename
+      INTEGER             ::   iyear, imonth, iday
+      REAL (wp)           ::   zsec
+      REAL (wp)           ::   zfjulday
+      CHARACTER(len=256)  :: cl_path
+      CHARACTER(len=256)  :: cl_filename
+      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character
       TYPE(iceberg), POINTER :: this
       TYPE(point)  , POINTER :: pt
@@ -240,8 +246,18 @@
       cl_path = TRIM(cn_ocerst_outdir)
       IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/'
+      IF ( ln_rstdate ) THEN
+         zfjulday = fjulday + rdttra(1) / rday
+         IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error
+         CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )           
+         WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
+      ELSE
+         IF( kt > 999999999 ) THEN   ;   WRITE(clkt, *       ) kt
+         ELSE                        ;   WRITE(clkt, '(i8.8)') kt
+         ENDIF
+      ENDIF
       IF( lk_mpp ) THEN
-         WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart_",I4.4,".nc")') TRIM(cexper), kt, narea-1
+         WRITE(cl_filename,'(A,"_icebergs_",A,"_restart_",I4.4,".nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)), narea-1
       ELSE
-         WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart.nc")') TRIM(cexper), kt
+         WRITE(cl_filename,'(A,"_icebergs_",A,"_restart.nc")') TRIM(cexper), TRIM(ADJUSTL(clkt))
       ENDIF
       IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM(cl_path)//TRIM(cl_filename)
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 8059)
@@ -30,4 +30,5 @@
    CHARACTER(lc) ::   cn_ocerst_outdir !: restart output directory
    LOGICAL       ::   ln_rstart        !: start from (F) rest or (T) a restart file
+   LOGICAL       ::   ln_rstdate       !: datestamping of restarts 
    LOGICAL       ::   ln_rst_list      !: output restarts at list of times (T) or by frequency (F)
    INTEGER       ::   nn_no            !: job number
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 8059)
@@ -21,4 +21,5 @@
    USE in_out_manager  ! I/O manager
    USE iom             ! I/O module
+   USE ioipsl, ONLY : ju2ymds    ! for calendar
    USE eosbn2          ! equation of state            (eos bn2 routine)
    USE trdmxl_oce      ! ocean active mixed layer tracers trends variables
@@ -54,8 +55,11 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! ocean time-step
+      INTEGER             ::   iyear, imonth, iday
+      REAL (wp)           ::   zsec
+      REAL (wp)           ::   zfjulday
       !!
       CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character
       CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
-      CHARACTER(lc)       ::   clpath   ! full path to ocean output restart file
+      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file
       !!----------------------------------------------------------------------
       !
@@ -81,7 +85,14 @@
       IF( kt == nitrst - 1 .OR. nstock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN
          IF( nitrst <= nitend .AND. nitrst > 0 ) THEN 
-            ! beware of the format used to write kt (default is i8.8, that should be large enough...)
-            IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst
-            ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst
+            IF ( ln_rstdate ) THEN
+               zfjulday = fjulday + rdttra(1) / rday
+               IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error
+               CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )           
+               WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
+            ELSE
+               ! beware of the format used to write kt (default is i8.8, that should be large enough...)
+               IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst
+               ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst
+               ENDIF
             ENDIF
             ! create the file
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 8059)
@@ -121,4 +121,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sprecip           !: solid precipitation                          [Kg/m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fr_i              !: ice fraction = 1 - lead fraction      (between 0 to 1)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   pressnow          !: UKMO SHELF pressure
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   apgu              !: UKMO SHELF pressure forcing
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   apgv              !: UKMO SHELF pressure forcing
 #if defined key_cpl_carbon_cycle
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm]
@@ -171,4 +174,5 @@
 #endif
          &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      &
+         &      pressnow(jpi,jpj), apgu(jpi,jpj)    , apgv(jpi,jpj) ,     &
          &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) )
          !
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 8059)
@@ -28,5 +28,5 @@
    PUBLIC sbc_flx       ! routine called by step.F90
 
-   INTEGER , PARAMETER ::   jpfld   = 5   ! maximum number of files to read 
+   INTEGER , PARAMETER ::   jpfld   = 6   ! maximum number of files to read 
    INTEGER , PARAMETER ::   jp_utau = 1   ! index of wind stress (i-component) file
    INTEGER , PARAMETER ::   jp_vtau = 2   ! index of wind stress (j-component) file
@@ -34,5 +34,8 @@
    INTEGER , PARAMETER ::   jp_qsr  = 4   ! index of solar heat file
    INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file
+   INTEGER , PARAMETER ::   jp_press = 6  ! index of pressure for UKMO shelf fluxes
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
+   LOGICAL , PUBLIC    ::   ln_shelf_flx = .FALSE. ! UKMO SHELF specific flux flag
+   INTEGER             ::   jpfld_local   ! maximum number of files to read (locally modified depending on ln_shelf_flx) 
 
    !! * Substitutions
@@ -82,9 +85,16 @@
       REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
       REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
+      REAL     ::   cs                    ! UKMO SHELF: Friction co-efficient at surface
+      REAL     ::   totwindspd            ! UKMO SHELF: Magnitude of wind speed vector
+
+      REAL(wp) ::   rhoa  = 1.22         ! Air density kg/m3
+      REAL(wp) ::   cdrag = 1.5e-3       ! drag coefficient 
       !!
       CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files
       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                    ! array of namelist information structures
-      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp  ! informations about the fields to be read
-      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp
+      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp, sn_press  !  informations about the fields to be read
+      LOGICAL     ::   ln_foam_flx  = .FALSE.                     ! UKMO FOAM specific flux flag
+      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp,   &
+      &                    ln_foam_flx, sn_press, ln_shelf_flx
       !!---------------------------------------------------------------------
       !
@@ -109,9 +119,18 @@
          slf_i(jp_emp ) = sn_emp
          !
-         ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure
+            ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure
+            IF( ln_shelf_flx ) slf_i(jp_press) = sn_press
+   
+            ! define local jpfld depending on shelf_flx logical
+            IF( ln_shelf_flx ) THEN
+               jpfld_local = jpfld
+            ELSE
+               jpfld_local = jpfld-1
+            ENDIF
+            !
          IF( ierror > 0 ) THEN   
             CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN  
          ENDIF
-         DO ji= 1, jpfld
+         DO ji= 1, jpfld_local
             ALLOCATE( sf(ji)%fnow(jpi,jpj,1) )
             IF( slf_i(ji)%ln_tint ) ALLOCATE( sf(ji)%fdta(jpi,jpj,1,2) )
@@ -132,10 +151,39 @@
          ENDIF
 !CDIR COLLAPSE
+            !!UKMO SHELF effect of atmospheric pressure on SSH
+            ! If using ln_apr_dyn, this is done there so don't repeat here.
+            IF( ln_shelf_flx .AND. .NOT. ln_apr_dyn) THEN
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1
+                     apgu(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji+1,jj,1)-sf(jp_press)%fnow(ji,jj,1))/e1u(ji,jj)
+                     apgv(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji,jj+1,1)-sf(jp_press)%fnow(ji,jj,1))/e2v(ji,jj)
+                  END DO
+               END DO
+            ENDIF ! ln_shelf_flx
+      
          DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
             DO ji = 1, jpi
-               utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
-               vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
-               qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1)
-               emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1)
+                   IF( ln_shelf_flx ) THEN
+                      !! UKMO SHELF - need atmospheric pressure to calculate Haney forcing
+                      pressnow(ji,jj) = sf(jp_press)%fnow(ji,jj,1)
+                      !! UKMO SHELF flux files contain wind speed not wind stress
+                      totwindspd = sqrt((sf(jp_utau)%fnow(ji,jj,1))**2.0 + (sf(jp_vtau)%fnow(ji,jj,1))**2.0)
+                      cs = 0.63 + (0.066 * totwindspd)
+                      utau(ji,jj) = cs * (rhoa/rau0) * sf(jp_utau)%fnow(ji,jj,1) * totwindspd
+                      vtau(ji,jj) = cs * (rhoa/rau0) * sf(jp_vtau)%fnow(ji,jj,1) * totwindspd
+                   ELSE
+                      utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
+                      vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
+                   ENDIF
+                   qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1)
+                   IF( ln_foam_flx .OR. ln_shelf_flx ) THEN
+                      !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot)
+                      qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1)
+                      !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P
+                      emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1)
+                   ELSE
+                      qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1)
+                      emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1)
+                   ENDIF
             END DO
          END DO
@@ -143,4 +191,11 @@
          qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
          !
+   
+            !! UKMO FOAM wind fluxes need lbc_lnk calls owing to a bug in interp.exe
+            IF( ln_foam_flx ) THEN
+               CALL lbc_lnk( utau(:,:), 'U', -1. )
+               CALL lbc_lnk( vtau(:,:), 'V', -1. )
+            ENDIF
+    
          !                                                        ! module of wind stress and wind speed at T-point
          zcoef = 1. / ( zrhoa * zcdrag )
@@ -162,5 +217,5 @@
             WRITE(numout,*) 
             WRITE(numout,*) '        read daily momentum, heat and freshwater fluxes OK'
-            DO jf = 1, jpfld
+            DO jf = 1, jpfld_local
                IF( jf == jp_utau .OR. jf == jp_vtau )   zfact =     1.
                IF( jf == jp_qtot .OR. jf == jp_qsr  )   zfact =     0.1
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 8059)
@@ -42,4 +42,5 @@
    LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term 
    REAL(wp)        ::   rn_sssr_bnd     ! ABS(Max./Min.) value of erp term [mm/day]
+   LOGICAL         ::   ln_UKMO_haney   ! UKMO specific flag to calculate Haney forcing  
 
    REAL(wp) , ALLOCATABLE, DIMENSION(:) ::   buffer   ! Temporary buffer for exchange
@@ -79,4 +80,12 @@
       INTEGER  ::   ierror   ! return error code
       !!
+      REAL(wp) ::   sst1,sst2                      ! sea surface temperature
+      REAL(wp) ::   e_sst1, e_sst2                 ! saturation vapour pressure
+      REAL(wp) ::   qs1,qs2                        ! specific humidity
+      REAL(wp) ::   pr_tmp                         ! temporary variable for pressure
+ 
+      REAL(wp), DIMENSION(jpi,jpj) ::  hny_frc1    ! Haney forcing for sensible heat, correction for latent heat   
+      REAL(wp), DIMENSION(jpi,jpj) ::  hny_frc2    ! Haney forcing for sensible heat, correction for latent heat   
+      !!
       CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
       TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read
@@ -95,11 +104,32 @@
             !
             IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) )
-                     qns(ji,jj) = qns(ji,jj) + zqrp
-                     qrp(ji,jj) = zqrp
-                  END DO
-               END DO
+                  IF( ln_UKMO_haney ) THEN
+                     DO jj = 1, jpj
+                        DO ji = 1, jpi
+                           sst1   =  sst_m(ji,jj)
+                           sst2   =  sf_sst(1)%fnow(ji,jj,1)   
+                           e_sst1 = 10**((0.7859+0.03477*sst1)/(1.+0.00412*sst1))
+                           e_sst2 = 10**((0.7859+0.03477*sst2)/(1.+0.00412*sst2))         
+                           pr_tmp = 0.01*pressnow(ji,jj)  !pr_tmp = 1012.0
+                           qs1    = (0.62197*e_sst1)/(pr_tmp-0.378*e_sst1)
+                           qs2    = (0.62197*e_sst2)/(pr_tmp-0.378*e_sst2)
+                           hny_frc1(ji,jj) = sst1-sst2                   
+                           hny_frc2(ji,jj) = qs1-qs2                     
+                          !Might need to mask off land points.
+                           hny_frc1(ji,jj)=-hny_frc1(ji,jj)*wndm(ji,jj)*1.42
+                           hny_frc2(ji,jj)=-hny_frc2(ji,jj)*wndm(ji,jj)*4688.0
+                           qns(ji,jj)=qns(ji,jj)+hny_frc1(ji,jj)+hny_frc2(ji,jj)   
+                           qrp(ji,jj) = 0.e0
+                        END DO
+                     END DO
+                  ELSE
+                     DO jj = 1, jpj
+                        DO ji = 1, jpi
+                           zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) )
+                           qns(ji,jj) = qns(ji,jj) + zqrp
+                           qrp(ji,jj) = zqrp
+                        END DO
+                     END DO
+                  ENDIF
                CALL iom_put( "qrp", qrp )                             ! heat flux damping
             ENDIF
@@ -163,5 +193,5 @@
       CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
       TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read
-      NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, sn_sss, ln_sssr_bnd, rn_sssr_bnd
+      NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, sn_sss, ln_sssr_bnd, rn_sssr_bnd, ln_UKMO_haney
       INTEGER     ::  ios
       !!----------------------------------------------------------------------
@@ -189,4 +219,5 @@
          WRITE(numout,*) '      flag to bound erp term                 ln_sssr_bnd = ', ln_sssr_bnd
          WRITE(numout,*) '      ABS(Max./Min.) erp threshold           rn_sssr_bnd = ', rn_sssr_bnd, ' mm/day'
+         WRITE(numout,*) '      Haney forcing                          ln_UKMO_haney = ', ln_UKMO_haney
       ENDIF
       !
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/tide.h90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/tide.h90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/tide.h90	(revision 8059)
@@ -28,2 +28,16 @@
    Wave(18) = tide(  'L2'     , 0.006694 ,    2   ,  2 , -1 ,  2 , -1 ,  0  , +180  ,  2   , -2   ,  0   ,  0   , 0 ,  215    )
    Wave(19) = tide(  'T2'     , 0.006614 ,    2   ,  2 ,  0 , -1 ,  0 ,  1  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    )
+   !
+   !             !! name_tide , equitide , nutide , nt , ns , nh , np , np1 , shift , nksi , nnu0 , nnu1 , nnu2 , R , formula !!
+   Wave(20) = tide(  'MNS2'   , 0.000000 ,    2   ,  2 , -5 ,  4 ,  1 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,    6    )
+   Wave(21) = tide(  'Lam2'   , 0.001760 ,    2   ,  2 , -1 ,  0 ,  1 ,  0  , +180  ,  2   , -2   ,  0   ,  0   , 0 ,   78    )
+   Wave(22) = tide(  'MSN2'   , 0.000000 ,    2   ,  2 ,  1 ,  0 ,  1 ,  0  ,    0  ,  2   , -2   ,  0   ,  2   , 0 ,    6    )
+   Wave(23) = tide(  '2SM2'   , 0.000000 ,    2   ,  2 ,  2 , -2 ,  0 ,  0  ,    0  , -2   ,  2   ,  0   ,  0   , 0 ,   16    )
+   Wave(24) = tide(  'MO3'    , 0.000000 ,    3   ,  3 , -4 ,  1 ,  0 ,  0  ,  +90  ,  2   , -2   ,  0   ,  0   , 0 ,   13    )
+   Wave(25) = tide(  'MK3'    , 0.000000 ,    3   ,  3 , -2 ,  3 ,  0 ,  0  ,  -90  ,  2   , -2   , -1   ,  0   , 0 ,   10    )
+   Wave(26) = tide(  'MN4'    , 0.000000 ,    4   ,  4 , -5 ,  4 ,  1 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,    1    )
+   Wave(27) = tide(  'MS4'    , 0.000000 ,    4   ,  4 , -2 ,  2 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    2    )
+   Wave(28) = tide(  'M6'     , 0.000000 ,    6   ,  6 , -6 ,  6 ,  0 ,  0  ,    0  ,  6   , -6   ,  0   ,  0   , 0 ,    4    )
+   Wave(29) = tide(  '2MS6'   , 0.000000 ,    6   ,  6 , -4 ,  4 ,  0 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,    6    )
+   Wave(30) = tide(  '2MK6'   , 0.000000 ,    6   ,  6 , -4 ,  6 ,  0 ,  0  ,    0  ,  4   , -4   ,  0   , -2   , 0 ,    5    )
+   Wave(31) = tide(  '3M2S2'  , 0.000000 ,    2   , 2  , -6 ,  6 ,  0 ,  0  ,    0  ,  6   , -6   ,  0   ,  0   , 0 ,   12    )
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90	(revision 8059)
@@ -16,8 +16,8 @@
    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
-      CHARACTER(LEN=4) ::   cname_tide
+      CHARACTER(LEN=5) ::   cname_tide
       REAL(wp)         ::   equitide
       INTEGER          ::   nutide
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 8059)
@@ -1241,5 +1241,6 @@
       IF(lwm) WRITE( numond, nameos )
       !
-      rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3]
+      rau0        = 1020._wp                 !: volumic mass of reference     [kg/m3]
+!     rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3]
       rcp         = 3991.86795711963_wp      !: heat capacity     [J/K]
       !
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 8059)
@@ -100,4 +100,8 @@
          IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
       ENDIF
+! slwa unless you use l_trdtra too, the above switches off trend calculations for l_trdtrc
+         l_trd = .FALSE.
+         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
+!slwa
       !
       IF( l_trd )  THEN
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/tradwl.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/tradwl.F90	(revision 8059)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/tradwl.F90	(revision 8059)
@@ -0,0 +1,238 @@
+MODULE tradwl
+   !!======================================================================
+   !!                       ***  MODULE  tradwl  ***
+   !! Ocean physics: solar radiation penetration in the top ocean levels
+   !!======================================================================
+   !! History :  POLCOMS  !  1996-10  (J. Holt)  Original code
+   !!   NEMO     3.2  !  2010-03  (E. O'Dea)  Import to Nemo for use in Shelf Model
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   tra_dwl      : trend due to the solar radiation penetration
+   !!   tra_dwl_init : solar radiation penetration initialization
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and active tracers
+   USE dom_oce         ! ocean space and time domain
+   USE sbc_oce         ! surface boundary condition: ocean
+   USE trc_oce         ! share SMS/Ocean variables
+   USE trd_oce         ! trends: ocean variables
+   USE trdtra          ! ocean active tracers trends 
+   USE in_out_manager  ! I/O manager
+   USE phycst          ! physical constants
+   USE prtctl          ! Print control
+   USE iom             ! I/O manager
+   USE fldread         ! read input fields
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   tra_dwl        ! routine called by step.F90 (ln_tradwl=T)
+
+   !                                           !!* Namelist namtra_qsr: penetrative solar radiation
+   LOGICAL , PUBLIC ::   ln_tradwl  = .TRUE.    ! light absorption (dwl) flag
+   LOGICAL , PUBLIC ::   ln_vary_lambda  = .TRUE.    ! vary Lambda or not flag
+   
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE tra_dwl( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_qsr  ***
+      !!
+      !! ** Purpose :   Compute the temperature trend due to the solar radiation
+      !!      penetration and add it to the general temperature trend.
+      !!
+      !! ** Method  : The profile of the solar radiation within the ocean is defined
+      !!
+      !!          Jason Holt Oct 96
+      !!
+      !!          Calculates change in temperature due to penetrating
+      !!          radiation, with cooling at the surface layer
+      !!
+      !!          rad=rad0*exp(lambda*z)
+      !!
+      !!       Heat input into box is between z=K and z=K+1 is RAD(K)-RAD(K+1)
+      !!
+      !!
+      !! ** Action  : - update ta with the penetrative solar radiation trend
+      !!              - save the trend in ttrd ('key_trdtra')
+      !!
+      !!----------------------------------------------------------------------
+      USE oce, ONLY :   ztrdt => ua   ! use ua as 3D workspace   
+      USE oce, ONLY :   ztrds => va   ! use va as 3D workspace   
+      !!
+      INTEGER, INTENT(in) ::   kt     ! ocean time-step
+      !!
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      INTEGER  ::   irgb                 ! temporary integers
+      REAL(wp) ::   zchl, zcoef, zsi0r   ! temporary scalars
+      REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr            ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0, ze1 , ze2, ze3, zea    ! 3D workspace
+      !!----------------------------------------------------------------------
+      !! HERE GO VARIABLES USED IN POLCOMS CLEAN UP LATER
+
+      integer i,j,k
+!      real*8 dtmp(n-1)
+      real*8 dtmp(jpkm1)
+      real*8 z1,z2,Rad0,Rad1,Rad2,rD,SurfOut,cp
+      logical first
+      save first
+      data first/.true./
+      !!--------------------------End of POLCOMS variables Note instead of using saves
+      !!--------------------------Could shift this into initial code
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_dwl : penetration of the surface solar radiation'
+         IF(lwp) WRITE(numout,*) '~~~~~~~'
+         CALL tra_dwl_init
+         IF( .NOT.ln_tradwl )   RETURN
+      ENDIF
+
+      IF( l_trdtra ) THEN      ! Save ta and sa trends
+         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
+         ztrds(:,:,:) = 0.e0
+      ENDIF
+!--------------------------------------------------------------------
+!  Set transmissivity
+!--------------------------------------------------------------------
+!
+!  Normal value
+!
+!---------------------------------------------------------------------------
+!
+!  Convert Heat fluxes to units used in POL subroutine dwl
+!
+!---------------------------------------------------------------------------
+      cp=3986.0d0
+
+    DO jj = 2, jpj
+         DO ji = fs_2, fs_jpim1
+           qsr(ji,jj)  = qsr(ji,jj)  * (r1_rau0_rcp)
+         ENDDO       !ji
+    ENDDO            !jj
+!--------------------------------------------------------------------------------
+      if ( first ) then
+        rlambda2(:,:) = 0.0
+        first=.false.
+           if ( ln_vary_lambda ) then
+!        do j=1,jesub    ! Original Polcoms style Loop
+!          do i=1,iesub  ! Original Polcoms style Loop
+
+        do jj=2,jpjm1
+          do ji = fs_2, fs_jpim1   ! vector opt.
+
+!            if(ipexb(i,j).ne. 0) then  (Mask, use Tmask instead)
+
+              rlambda2(ji,jj)=-0.033*log(hbatt(ji,jj))+0.2583    ! JIAs formula
+              rlambda2(ji,jj)=max(0.05,rlambda2(ji,jj))     ! limit in deep water
+          enddo ! ji
+        enddo ! jj
+      rlambda = 0.0
+       else
+      rLambda=0.154
+       endif ! If vary lambda
+      endif ! If first
+
+!      do j=1,jesub      ! Original Polcoms Style Loop
+!        do i=1,iesub    ! Original Polcoms Style Loop
+              do jk=2,jpk
+        DO jj=2,jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+
+!--------------------------------------------------------------------
+! Calculate change in temperature
+!--------------------------------------------------------------------
+!
+!        rad0 = hfl_in(i,j)   ! change hfl_in to qsr I assume
+
+               rad0 = qsr(ji,jj)
+                     rD = rLambda2(ji,jj)  +rLambda      !  Transmissivity to be used here
+!       if rlambda 0 then rlambda2 not zer and vica versa 
+
+                z2=gdepw_0(ji,jj,jk-1)    ! grid box is from z=z1 to z=z2
+                z1=gdepw_0(ji,jj,jk)
+
+                Rad2=Rad0*(exp(-z2*rD)) ! radiation entering box
+                Rad1=Rad0*(exp(-z1*rD)) ! radiation leaving box
+
+
+                dtmp(jk)=1.0/(e3t_0(ji,jj,jk))*(Rad2-Rad1) !change in temperature
+                tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + dtmp(jk)
+         enddo  ! ji
+         enddo  ! jj
+             enddo !jk
+
+
+      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
+         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
+         !CEODCALL trd_mod( ztrdt, ztrds, jptra_trd_qsr, 'TRA', kt )
+      ENDIF
+      !                       ! print mean trends (used for debugging)
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
+      !
+   END SUBROUTINE tra_dwl
+
+
+   SUBROUTINE tra_dwl_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_dwl_init  ***
+      !!
+      !! ** Purpose :   Initialization for the penetrative solar radiation for Downwell routine
+      !!
+      !! ** Method  :   The profile of solar radiation within the ocean is set
+      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
+      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
+      !!      default values correspond to clear water (type I in Jerlov' 
+      !!      (1968) classification.
+      !!         called by tra_qsr at the first timestep (nit000)
+      !!
+      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
+      !!
+      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk            ! dummy loop indices
+      INTEGER  ::   ios                   ! Local integer output status for namelist read
+      INTEGER  ::   irgb, ierror          ! temporary integer
+      INTEGER  ::   ioptio, nqsr          ! temporary integer
+      REAL(wp) ::   zc0  , zc1            ! temporary scalars
+      REAL(wp) ::   zc2  , zc3  , zchl    !    -         -
+      REAL(wp) ::   zsi0r, zsi1r, zcoef   !    -         -
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr              ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0 , ze1 , ze2 , ze3 , zea   ! 3D workspace
+      !!
+      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
+      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
+      NAMELIST/namtra_dwl/  ln_tradwl, ln_vary_lambda
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam_ref )            ! Read Namelist namtra_dwl in reference namelist :
+      READ  ( numnam_ref, namtra_dwl, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )
+      
+      REWIND( numnam_cfg )            ! Read Namelist namtra_dwl in configuration namelist :
+      READ  ( numnam_cfg, namtra_dwl, IOSTAT = ios, ERR = 902)
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp )
+      !
+      IF(lwp) THEN                ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'tra_dwl_init : '
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namtra_dwl : set the parameter of penetration'
+         WRITE(numout,*) '      Light penetration (T) or not (F)         ln_tradwl  = ', ln_tradwl
+         WRITE(numout,*) '      Vary Lambda  (T) or not (F))             ln_vary_lambda  = ', ln_vary_lambda
+      ENDIF
+
+   END SUBROUTINE tra_dwl_init
+
+   !!======================================================================
+END MODULE tradwl
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 8059)
@@ -58,4 +58,5 @@
 
    REAL(wp) ::   rbcp   ! Brown & Campana parameters for semi-implicit hpg
+   INTEGER  ::   warn_1, warn_2   ! indicators for warning statement
 
    !! * Substitutions
@@ -93,6 +94,6 @@
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
-      INTEGER  ::   jk, jn    ! dummy loop indices
-      REAL(wp) ::   zfact     ! local scalars
+      INTEGER  ::   jk, jn, ji, jj     ! dummy loop indices
+      REAL(wp) ::   zfact, zfreeze     ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
       !!----------------------------------------------------------------------
@@ -120,5 +121,34 @@
       IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
 #endif
- 
+
+#if ( ! defined key_lim3 && ! defined key_lim2 && ! key_cice )
+      IF ( kt == nit000 ) warn_1=0
+      warn_2=0
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( tsa(ji,jj,jk,jp_tem) .lt. 0.0 ) THEN
+                  ! calculate freezing point
+                  zfreeze = ( -0.0575_wp + 1.710523E-3 * Sqrt(Abs(tsn(ji,jj,jk,jp_sal)))   & 
+                            - 2.154996E-4 * tsn(ji,jj,jk,jp_sal) ) * tsn(ji,jj,jk,jp_sal) - 7.53E-4 * ( 10.0_wp + fsdept(ji,jj,jk) )
+                  IF ( tsa(ji,jj,jk,jp_tem) .lt. zfreeze ) THEN
+                     tsa(ji,jj,jk,jp_tem)=zfreeze
+                     warn_2=1
+                  ENDIF
+               ENDIF
+            END DO
+         END DO
+      END DO
+      CALL mpp_max(warn_1)
+      CALL mpp_max(warn_2)
+      IF ( (warn_1 == 0) .and. (warn_2 /= 0) ) THEN
+         IF(lwp) THEN
+            CALL ctl_warn( ' Temperatures dropping below freezing point, ', &
+                      &    ' being forced to freezing point, no longer conservative' ) 
+         ENDIF
+         warn_1=1
+      ENDIF
+#endif
+
       ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler)
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 8059)
@@ -46,4 +46,5 @@
    LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem)
    INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0)
+   INTEGER , PUBLIC ::   nn_kd490dta  !: use kd490dta data (=1) or not (=0)
    REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands)
    REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
@@ -54,4 +55,5 @@
    REAL(wp) ::   xsi1r                           !: inverse of rn_si1
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_kd490 ! structure of input kd490 (file informations, fields read)
    INTEGER, PUBLIC ::   nksr              ! levels below which the light cannot penetrate ( depth larger than 391 m)
    REAL(wp), DIMENSION(3,61) ::   rkrgb   !: tabulated attenuation coefficients for RGB absorption
@@ -306,4 +308,48 @@
             !
          ENDIF
+! slwa
+         IF( nn_kd490dta == 1 ) THEN                      !  use KD490 data read in   !
+            !                                             ! ------------------------- !
+               nksr = jpk - 1
+               !
+               CALL fld_read( kt, 1, sf_kd490 )     ! Read kd490 data and provide it at the current time step
+               !
+               zcoef  = ( 1. - rn_abs )
+               ze0(:,:,1) = rn_abs  * qsr(:,:)
+               ze1(:,:,1) = zcoef * qsr(:,:)
+               zea(:,:,1) =         qsr(:,:)
+               !
+               DO jk = 2, nksr+1
+!CDIR NOVERRCHK
+                  DO jj = 1, jpj
+!CDIR NOVERRCHK   
+                     DO ji = 1, jpi
+                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     )
+                        zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * sf_kd490(1)%fnow(ji,jj,1) )
+                        ze0(ji,jj,jk) = zc0
+                        ze1(ji,jj,jk) = zc1
+                        zea(ji,jj,jk) = ( zc0 + zc1 ) * tmask(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO
+               ! clem: store attenuation coefficient of the first ocean level
+               IF ( ln_qsr_ice ) THEN
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi
+                        zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r     )
+                        zzc1 = zcoef  * EXP( - fse3t(ji,jj,1) * sf_kd490(1)%fnow(ji,jj,1) )
+                        fraqsr_1lev(ji,jj) = 1.0 - ( zzc0 + zzc1 ) * tmask(ji,jj,2) 
+                     END DO
+                  END DO
+               ENDIF
+               !
+               DO jk = 1, nksr                                        ! compute and add qsr trend to ta
+                  qsr_hc(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) )
+               END DO
+               zea(:,:,nksr+1:jpk) = 0.e0     ! 
+               CALL iom_put( 'qsr3d', zea )   ! Shortwave Radiation 3D distribution
+               !
+        ENDIF   ! use KD490 data
+!slwa
          !
          !                                        Add to the general trend
@@ -374,7 +420,8 @@
       CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
       TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
-      !!
-      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  &
-         &                  nn_chldta, rn_abs, rn_si0, rn_si1
+      TYPE(FLD_N)        ::   sn_kd490 ! informations about the kd490 field to be read
+      !!
+      NAMELIST/namtra_qsr/  sn_chl, sn_kd490, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  &
+         &                  nn_chldta, rn_abs, rn_si0, rn_si1, nn_kd490dta
       !!----------------------------------------------------------------------
 
@@ -409,4 +456,5 @@
          WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0
          WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1 = ', rn_si1
+         WRITE(numout,*) '      read in KD490 data                       nn_kd490dta  = ', nn_kd490dta
       ENDIF
 
@@ -422,4 +470,5 @@
          IF( ln_qsr_2bd  )   ioptio = ioptio + 1
          IF( ln_qsr_bio  )   ioptio = ioptio + 1
+         IF( nn_kd490dta == 1 )   ioptio = ioptio + 1
          !
          IF( ioptio /= 1 ) &
@@ -431,4 +480,5 @@
          IF( ln_qsr_2bd                      )   nqsr =  3
          IF( ln_qsr_bio                      )   nqsr =  4
+         IF( nn_kd490dta == 1                )   nqsr =  5
          !
          IF(lwp) THEN                   ! Print the choice
@@ -438,4 +488,5 @@
             IF( nqsr ==  3 )   WRITE(numout,*) '         2 bands light penetration'
             IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration'
+            IF( nqsr ==  5 )   WRITE(numout,*) '         KD490 light penetration'
          ENDIF
          !
@@ -447,6 +498,18 @@
          xsi0r = 1.e0 / rn_si0
          xsi1r = 1.e0 / rn_si1
+         IF( nn_kd490dta == 1 ) THEN           !* KD490 data : set sf_kd490 structure
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '        KD490 read in a file'
+            ALLOCATE( sf_kd490(1), STAT=ierror )
+            IF( ierror > 0 ) THEN
+               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_kd490 structure' )   ;   RETURN
+            ENDIF
+            ALLOCATE( sf_kd490(1)%fnow(jpi,jpj,1)   )
+            IF( sn_kd490%ln_tint )ALLOCATE( sf_kd490(1)%fdta(jpi,jpj,1,2) )
+            !                                        ! fill sf_kd490 with sn_kd490 and control print
+            CALL fld_fill( sf_kd490, (/ sn_kd490 /), cn_dir, 'tra_qsr_init',   &
+               &                                         'Solar penetration function of read KD490', 'namtra_qsr' )
          !                                ! ---------------------------------- !
-         IF( ln_qsr_rgb ) THEN            !  Red-Green-Blue light penetration  !
+         ELSEIF( ln_qsr_rgb ) THEN            !  Red-Green-Blue light penetration  !
             !                             ! ---------------------------------- !
             !
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 8059)
@@ -25,4 +25,5 @@
    USE trd_oce         ! trends: ocean variables
    USE trdtra          ! trends manager: tracers 
+   USE tradwl          ! solar radiation penetration (downwell method)
    !
    USE in_out_manager  ! I/O manager
@@ -138,5 +139,5 @@
 
 !!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
-      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
+      IF( .NOT.ln_traqsr .and. .NOT.ln_tradwl ) THEN     ! no solar radiation penetration
          qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
          qsr(:,:) = 0.e0                     ! qsr set to zero
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 8059)
@@ -203,7 +203,12 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
+#if defined key_tracer_budget
+!              ptrd(ji,jj,jk) = - (     pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)  )  * tmask(ji,jj,jk)
+               ptrd(ji,jj,jk) = -      pf (ji,jj,jk) * tmask(ji,jj,jk)
+#else
                ptrd(ji,jj,jk) = - (     pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)                        &
                  &                  - ( pun(ji,jj,jk) - pun(ji-ii,jj-ij,jk-ik) ) * ptn(ji,jj,jk)  )   &
                  &              / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )  * tmask(ji,jj,jk)
+#endif
             END DO
          END DO
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 8059)
@@ -18,4 +18,5 @@
    USE phycst          ! physical constants
    USE iom             ! I/O library
+   USE eosbn2          ! for zdf_mxl_zint
    USE lib_mpp         ! MPP library
    USE wrk_nemo        ! work arrays
@@ -27,5 +28,4 @@
 
    PUBLIC   zdf_mxl       ! called by step.F90
-   PUBLIC   zdf_mxl_alloc ! Used in zdf_tke_init
 
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP)
@@ -33,7 +33,18 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m]
+   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m] 
+   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   htc_mld    ! Heat content of hmld_zint
+   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: ll_found   ! Is T_b to be found by interpolation ? 
+   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: ll_belowml ! Flag points below mixed layer when ll_found=F
 
    REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
    REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth
+
+   TYPE, PUBLIC :: MXL_ZINT   !: Structure for MLD defs
+      INTEGER   :: mld_type   ! mixed layer type     
+      REAL(wp)  :: zref       ! depth of initial T_ref
+      REAL(wp)  :: dT_crit    ! Critical temp diff
+      REAL(wp)  :: iso_frac   ! Fraction of rn_dT_crit used
+   END TYPE MXL_ZINT
 
    !! * Substitutions
@@ -52,5 +63,7 @@
       zdf_mxl_alloc = 0      ! set to zero if no array to be allocated
       IF( .NOT. ALLOCATED( nmln ) ) THEN
-         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc )
+         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj),       &
+        &          htc_mld(jpi,jpj),                                                                    &
+        &          ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc )
          !
          IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc )
@@ -128,7 +141,7 @@
             iikn = nmln(ji,jj)
             imkt = mikt(ji,jj)
-            hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! Turbocline depth 
-            hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! Mixed layer depth
-            hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
+            hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! Turbocline depth 
+            hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj)    ! Mixed layer depth
+            hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
          END DO
       END DO
@@ -138,4 +151,7 @@
       ENDIF
       
+      ! Vertically-interpolated mixed-layer depth diagnostic
+      CALL zdf_mxl_zint( kt )
+
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 )
       !
@@ -145,4 +161,344 @@
       !
    END SUBROUTINE zdf_mxl
+
+   SUBROUTINE zdf_mxl_zint_mld( sf ) 
+      !!---------------------------------------------------------------------------------- 
+      !!                    ***  ROUTINE zdf_mxl_zint_mld  *** 
+      !                                                                        
+      !   Calculate vertically-interpolated mixed layer depth diagnostic. 
+      !            
+      !   This routine can calculate the mixed layer depth diagnostic suggested by
+      !   Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate
+      !   vertically-interpolated mixed-layer depth diagnostics with other parameter
+      !   settings set in the namzdf_mldzint namelist.  
+      ! 
+      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the  
+      !   density has increased by an amount equivalent to a temperature difference of  
+      !   0.8C at the surface. 
+      ! 
+      !   For other values of mld_type the mixed layer is calculated as the depth at  
+      !   which the temperature differs by 0.8C from the surface temperature.  
+      !                                                                        
+      !   David Acreman, Daley Calvert                                      
+      ! 
+      !!----------------------------------------------------------------------------------- 
+
+      TYPE(MXL_ZINT), INTENT(in)  :: sf
+
+      ! Diagnostic criteria
+      INTEGER   :: nn_mld_type   ! mixed layer type     
+      REAL(wp)  :: rn_zref       ! depth of initial T_ref
+      REAL(wp)  :: rn_dT_crit    ! Critical temp diff
+      REAL(wp)  :: rn_iso_frac   ! Fraction of rn_dT_crit used
+
+      ! Local variables
+      REAL(wp), PARAMETER :: zepsilon = 1.e-30          ! local small value
+      INTEGER, POINTER, DIMENSION(:,:) :: ikmt          ! number of active tracer levels 
+      INTEGER, POINTER, DIMENSION(:,:) :: ik_ref        ! index of reference level 
+      INTEGER, POINTER, DIMENSION(:,:) :: ik_iso        ! index of last uniform temp level 
+      REAL, POINTER, DIMENSION(:,:,:)  :: zT            ! Temperature or density 
+      REAL, POINTER, DIMENSION(:,:)    :: ppzdep        ! depth for use in calculating d(rho) 
+      REAL, POINTER, DIMENSION(:,:)    :: zT_ref        ! reference temperature 
+      REAL    :: zT_b                                   ! base temperature 
+      REAL, POINTER, DIMENSION(:,:,:)  :: zdTdz         ! gradient of zT 
+      REAL, POINTER, DIMENSION(:,:,:)  :: zmoddT        ! Absolute temperature difference 
+      REAL    :: zdz                                    ! depth difference 
+      REAL    :: zdT                                    ! temperature difference 
+      REAL, POINTER, DIMENSION(:,:)    :: zdelta_T      ! difference critereon 
+      REAL, POINTER, DIMENSION(:,:)    :: zRHO1, zRHO2  ! Densities 
+      INTEGER :: ji, jj, jk                             ! loop counter 
+
+      !!------------------------------------------------------------------------------------- 
+      !  
+      CALL wrk_alloc( jpi, jpj, ikmt, ik_ref, ik_iso) 
+      CALL wrk_alloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 ) 
+      CALL wrk_alloc( jpi, jpj, jpk, zT, zdTdz, zmoddT ) 
+
+      ! Unpack structure
+      nn_mld_type = sf%mld_type
+      rn_zref     = sf%zref
+      rn_dT_crit  = sf%dT_crit
+      rn_iso_frac = sf%iso_frac
+
+      ! Set the mixed layer depth criterion at each grid point 
+      IF( nn_mld_type == 0 ) THEN
+         zdelta_T(:,:) = rn_dT_crit
+         zT(:,:,:) = rhop(:,:,:)
+      ELSE IF( nn_mld_type == 1 ) THEN
+         ppzdep(:,:)=0.0 
+         call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) ) 
+! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST 
+! [assumes number of tracers less than number of vertical levels] 
+         zT(:,:,1:jpts)=tsn(:,:,1,1:jpts) 
+         zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit 
+         CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) ) 
+         zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 
+         ! RHO from eos (2d version) doesn't calculate north or east halo: 
+         CALL lbc_lnk( zdelta_T, 'T', 1. ) 
+         zT(:,:,:) = rhop(:,:,:) 
+      ELSE 
+         zdelta_T(:,:) = rn_dT_crit                      
+         zT(:,:,:) = tsn(:,:,:,jp_tem)                           
+      END IF 
+
+      ! Calculate the gradient of zT and absolute difference for use later 
+      DO jk = 1 ,jpk-2 
+         zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1) 
+         zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 
+      END DO 
+
+      ! Find density/temperature at the reference level (Kara et al use 10m).          
+      ! ik_ref is the index of the box centre immediately above or at the reference level 
+      ! Find rn_zref in the array of model level depths and find the ref    
+      ! density/temperature by linear interpolation.                                   
+      DO jk = jpkm1, 2, -1 
+         WHERE ( fsdept(:,:,jk) > rn_zref ) 
+           ik_ref(:,:) = jk - 1 
+           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - fsdept(:,:,jk-1) ) 
+         END WHERE 
+      END DO 
+
+      ! If the first grid box centre is below the reference level then use the 
+      ! top model level to get zT_ref 
+      WHERE ( fsdept(:,:,1) > rn_zref )  
+         zT_ref = zT(:,:,1) 
+         ik_ref = 1 
+      END WHERE 
+
+      ! The number of active tracer levels is 1 less than the number of active w levels 
+      ikmt(:,:) = mbathy(:,:) - 1 
+
+      ! Initialize / reset
+      ll_found(:,:) = .false.
+
+      IF ( rn_iso_frac - zepsilon > 0. ) THEN
+         ! Search for a uniform density/temperature region where adjacent levels          
+         ! differ by less than rn_iso_frac * deltaT.                                      
+         ! ik_iso is the index of the last level in the uniform layer  
+         ! ll_found indicates whether the mixed layer depth can be found by interpolation 
+         ik_iso(:,:)   = ik_ref(:,:) 
+         DO jj = 1, nlcj 
+            DO ji = 1, nlci 
+!CDIR NOVECTOR 
+               DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1 
+                  IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN 
+                     ik_iso(ji,jj)   = jk 
+                     ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 
+                     EXIT 
+                  END IF 
+               END DO 
+            END DO 
+         END DO 
+
+         ! Use linear interpolation to find depth of mixed layer base where possible 
+         hmld_zint(:,:) = rn_zref 
+         DO jj = 1, jpj 
+            DO ji = 1, jpi 
+               IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN 
+                  zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 
+                  hmld_zint(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz 
+               END IF 
+            END DO 
+         END DO 
+      END IF
+
+      ! If ll_found = .false. then calculate MLD using difference of zdelta_T    
+      ! from the reference density/temperature 
+ 
+! Prevent this section from working on land points 
+      WHERE ( tmask(:,:,1) /= 1.0 ) 
+         ll_found = .true. 
+      END WHERE 
+ 
+      DO jk=1, jpk 
+         ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)  
+      END DO 
+ 
+! Set default value where interpolation cannot be used (ll_found=false)  
+      DO jj = 1, jpj 
+         DO ji = 1, jpi 
+            IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = fsdept(ji,jj,ikmt(ji,jj)) 
+         END DO 
+      END DO 
+
+      DO jj = 1, jpj 
+         DO ji = 1, jpi 
+!CDIR NOVECTOR 
+            DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj) 
+               IF ( ll_found(ji,jj) ) EXIT 
+               IF ( ll_belowml(ji,jj,jk) ) THEN                
+                  zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) ) 
+                  zdT  = zT_b - zT(ji,jj,jk-1)                                      
+                  zdz  = zdT / zdTdz(ji,jj,jk-1)                                       
+                  hmld_zint(ji,jj) = fsdept(ji,jj,jk-1) + zdz 
+                  EXIT                                                   
+               END IF 
+            END DO 
+         END DO 
+      END DO 
+
+      hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1) 
+      !  
+      CALL wrk_dealloc( jpi, jpj, ikmt, ik_ref, ik_iso) 
+      CALL wrk_dealloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 ) 
+      CALL wrk_dealloc( jpi,jpj, jpk, zT, zdTdz, zmoddT ) 
+      ! 
+   END SUBROUTINE zdf_mxl_zint_mld
+
+   SUBROUTINE zdf_mxl_zint_htc( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_mxl_zint_htc  ***
+      !! 
+      !! ** Purpose :   
+      !!
+      !! ** Method  :   
+      !!----------------------------------------------------------------------
+
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+
+      INTEGER :: ji, jj, jk
+      INTEGER :: ikmax
+      REAL(wp) :: zc, zcoef
+      !
+      INTEGER,  ALLOCATABLE, DIMENSION(:,:) ::   ilevel
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zthick_0, zthick
+
+      !!----------------------------------------------------------------------
+
+      IF( .NOT. ALLOCATED(ilevel) ) THEN
+         ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), &
+         &         zthick(jpi,jpj), STAT=ji )
+         IF( lk_mpp  )   CALL mpp_sum(ji)
+         IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' )
+      ENDIF
+
+      ! Find last whole model T level above the MLD
+      ilevel(:,:)   = 0
+      zthick_0(:,:) = 0._wp
+
+      DO jk = 1, jpkm1  
+         DO jj = 1, jpj
+            DO ji = 1, jpi                    
+               zthick_0(ji,jj) = zthick_0(ji,jj) + fse3t(ji,jj,jk)
+               IF( zthick_0(ji,jj) < hmld_zint(ji,jj) )   ilevel(ji,jj) = jk
+            END DO
+         END DO
+         WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2)
+         WRITE(numout,*) 'fsdepw(jk+1 =',jk+1,') =',fsdepw(2,2,jk+1)
+      END DO
+
+      ! Surface boundary condition
+      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc_mld(:,:) = 0._wp                                   
+      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)   
+      ENDIF
+
+      ! Deepest whole T level above the MLD
+      ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 )
+
+      ! Integration down to last whole model T level
+      DO jk = 1, ikmax
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zc = fse3t(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1  )  )    ! 0 below ilevel
+               zthick(ji,jj) = zthick(ji,jj) + zc
+               htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      ! Subsequent partial T level
+      zthick(:,:) = hmld_zint(:,:) - zthick(:,:)   !   remaining thickness to reach MLD
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem)  & 
+      &                      * MIN( fse3t(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1)
+         END DO
+      END DO
+
+      WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2)
+
+      ! Convert to heat content
+      zcoef = rau0 * rcp
+      htc_mld(:,:) = zcoef * htc_mld(:,:)
+
+   END SUBROUTINE zdf_mxl_zint_htc
+
+   SUBROUTINE zdf_mxl_zint( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_mxl_zint  ***
+      !! 
+      !! ** Purpose :   
+      !!
+      !! ** Method  :   
+      !!----------------------------------------------------------------------
+
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+
+      INTEGER :: ios
+      INTEGER :: jn
+
+      INTEGER :: nn_mld_diag = 0    ! number of diagnostics
+
+      CHARACTER(len=1) :: cmld
+
+      TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
+      TYPE(MXL_ZINT), SAVE, DIMENSION(5) ::   mld_diags
+
+      NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
+
+      !!----------------------------------------------------------------------
+      
+      IF( kt == nit000 ) THEN
+         REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist 
+         READ  ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901)
+901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist', lwp )
+
+         REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist 
+         READ  ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 )
+902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist', lwp )
+         IF(lwm) WRITE ( numond, namzdf_mldzint )
+
+         IF( nn_mld_diag > 5 )   CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' )
+
+         mld_diags(1) = sn_mld1
+         mld_diags(2) = sn_mld2
+         mld_diags(3) = sn_mld3
+         mld_diags(4) = sn_mld4
+         mld_diags(5) = sn_mld5
+
+         IF( nn_mld_diag > 0 ) THEN
+            WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================'
+            WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)'
+            DO jn = 1, nn_mld_diag
+               WRITE(numout,*) 'MLD criterion',jn,':'
+               WRITE(numout,*) '    nn_mld_type =', mld_diags(jn)%mld_type
+               WRITE(numout,*) '    rn_zref ='    , mld_diags(jn)%zref
+               WRITE(numout,*) '    rn_dT_crit =' , mld_diags(jn)%dT_crit
+               WRITE(numout,*) '    rn_iso_frac =', mld_diags(jn)%iso_frac
+            END DO
+            WRITE(numout,*) '===================================================================='
+         ENDIF
+      ENDIF
+
+      IF( nn_mld_diag > 0 ) THEN
+         DO jn = 1, nn_mld_diag
+            WRITE(cmld,'(I1)') jn
+            IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN
+               CALL zdf_mxl_zint_mld( mld_diags(jn) )
+
+               IF( iom_use( "mldzint_"//cmld ) ) THEN
+                  CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) )
+               ENDIF
+
+               IF( iom_use( "mldhtc_"//cmld ) )  THEN
+                  CALL zdf_mxl_zint_htc( kt )
+                  CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:)   )
+               ENDIF
+            ENDIF
+         END DO
+      ENDIF
+
+   END SUBROUTINE zdf_mxl_zint
 
    !!======================================================================
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 8059)
@@ -85,4 +85,6 @@
    USE stopar
    USE stopts
+   USE diatmb          ! Top,middle,bottom output
+   USE dia25h          ! 25h mean output
 
    IMPLICIT NONE
@@ -475,4 +477,6 @@
       IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments
       IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler
+                            CALL dia_tmb_init  ! TMB outputs
+                            CALL dia_25h_init  ! 25h mean  outputs
       !
    END SUBROUTINE nemo_init
@@ -630,4 +634,5 @@
       USE ldftra_oce, ONLY: ldftra_oce_alloc
       USE trc_oce   , ONLY: trc_oce_alloc
+      USE diainsitutem, ONLY: insitu_tem_alloc
 #if defined key_diadct 
       USE diadct    , ONLY: diadct_alloc 
@@ -646,4 +651,5 @@
       ierr = ierr + ldftra_oce_alloc()          ! ocean lateral  physics : tracers
       ierr = ierr + zdf_oce_alloc   ()          ! ocean vertical physics
+      ierr = ierr + insitu_tem_alloc()
       !
       ierr = ierr + trc_oce_alloc   ()          ! shared TRC / TRA arrays
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8059)
@@ -255,4 +255,5 @@
                              CALL tra_sbc    ( kstp )       ! surface boundary condition
       IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
+      IF( ln_tradwl      )   CALL tra_dwl    ( kstp )       ! Polcoms Style Short Wave Radiation 
       IF( ln_trabbc      )   CALL tra_bbc    ( kstp )       ! bottom heat flux
       IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme
@@ -337,5 +338,4 @@
       IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
       !
-      IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file
 
 #if defined key_agrif
@@ -361,4 +361,5 @@
                                CALL dia_wri_state( 'output.abort', kstp )
       ENDIF
+      IF( ln_harm_ana_store   )   CALL harm_ana( kstp )        ! Harmonic analysis of tides 
       IF( kstp == nit000   )   THEN
                  CALL iom_close( numror )     ! close input  ocean restart file
@@ -366,4 +367,5 @@
          IF( lwm.AND.numoni /= -1 ) CALL FLUSH    ( numoni )     ! flush output namelist ice
       ENDIF
+      IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 8059)
@@ -25,4 +25,5 @@
    USE sbcrnf           ! surface boundary condition: runoff variables
    USE sbccpl           ! surface boundary condition: coupled formulation (call send at end of step)
+   USE sbcflx           ! surface boundary condition: Fluxes
    USE sbc_oce          ! surface boundary condition: ocean
    USE sbctide          ! Tide initialisation
@@ -30,4 +31,5 @@
 
    USE traqsr           ! solar radiation penetration      (tra_qsr routine)
+   USE tradwl           ! POLCOMS style solar radiation    (tra_dwl routine) 
    USE trasbc           ! surface boundary condition       (tra_sbc routine)
    USE trabbc           ! bottom boundary condition        (tra_bbc routine)
@@ -106,4 +108,6 @@
    USE prtctl           ! Print control                    (prt_ctl routine)
 
+   USE harmonic_analysis ! harmonic analysis of tides (harm_ana routine) 
+   USE bdytides          ! harmonic analysis of tides (harm_ana routine) 
    USE diaobs           ! Observation operator
 
Index: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90	(revision 8058)
+++ branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90	(revision 8059)
@@ -27,4 +27,6 @@
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   etot3         !: light absortion coefficient
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   facvol        !: volume for degraded regions
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   rlambda2      !: Lambda2 for downwell version of Short wave Radiation
+   REAL(wp), PUBLIC                                      ::   rlambda       !: Lambda  for downwell version of Short wave Radiation
 
 #if defined key_top 
@@ -78,12 +80,14 @@
       !!                  ***  trc_oce_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER ::   ierr(2)        ! Local variables
+      INTEGER ::   ierr(3)        ! Local variables
       !!----------------------------------------------------------------------
       ierr(:) = 0
                      ALLOCATE( etot3 (jpi,jpj,jpk), STAT=ierr(1) )
       IF( lk_degrad) ALLOCATE( facvol(jpi,jpj,jpk), STAT=ierr(2) )
+                    ALLOCATE( rlambda2(jpi,jpj),   STAT=ierr(3) )
       trc_oce_alloc  = MAXVAL( ierr )
       !
       IF( trc_oce_alloc /= 0 )   CALL ctl_warn('trc_oce_alloc: failed to allocate etot3 array')
+      IF( trc_oce_alloc /= 0 )   CALL ctl_warn('trc_oce_alloc: failed to allocate etot3, facvol or rlambda2 array')
    END FUNCTION trc_oce_alloc
 
