Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diaprod.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diaprod.F90	(revision 6491)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diaprod.F90	(revision 6491)
@@ -0,0 +1,231 @@
+MODULE diaprod
+! Requires key_iom_put
+# if defined key_iomput
+   !!======================================================================
+   !!                     ***  MODULE  diaprod  ***
+   !! Ocean diagnostics :  write ocean product diagnostics
+   !!=====================================================================
+   !! History :  3.4  ! 2012  (D. Storkey)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dia_prod      : calculate and write out product diagnostics
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers 
+   USE dom_oce         ! ocean space and time domain
+   USE domvvl          ! for thickness weighted diagnostics if key_vvl
+   USE eosbn2          ! equation of state  (eos call)
+   USE phycst          ! physical constants
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager  ! I/O manager
+   USE diadimg         ! dimg direct access file format output
+   USE iom
+   USE ioipsl
+   USE lib_mpp         ! MPP library
+   USE timing          ! preformance summary
+   USE wrk_nemo        ! working array
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   dia_prod                 ! routines called by step.F90
+
+   !! * Substitutions
+#  include "zdfddm_substitute.h90"
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.4 , NEMO Consortium (2012)
+   !! $Id $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_prod( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_prod  ***
+      !!                   
+      !! ** Purpose :   Write out product diagnostics (uT, vS etc.)
+      !!
+      !! ** Method  :  use iom_put
+      !!               Product diagnostics are not thickness-weighted in this routine.
+      !!               They should be thickness-weighted using XIOS if key_vvl is set. 
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      !!
+      INTEGER                      ::   ji, jj, jk              ! dummy loop indices
+      REAL(wp)                     ::   zztmp, zztmpx, zztmpy   ! 
+      !!
+      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zrhop    ! potential density
+      !!----------------------------------------------------------------------
+      ! 
+      IF( nn_timing == 1 )   CALL timing_start('dia_prod')
+      ! 
+      CALL wrk_alloc( jpi , jpj      , z2d )
+      CALL wrk_alloc( jpi , jpj, jpk , z3d )
+      CALL wrk_alloc( jpi , jpj, jpk , zrhop )
+      !
+
+      IF( iom_use("urhop") .OR. iom_use("vrhop") .OR. iom_use("wrhop") &
+#if ! defined key_diaar5
+     &  .OR. iom_use("rhop") &
+#endif
+     & ) THEN
+         CALL eos( tsn, z3d, zrhop )                 ! now in situ and potential density
+         zrhop(:,:,:) = zrhop(:,:,:)-1000.e0         ! reference potential density to 1000 to avoid precision issues in rhop2 calculation
+         zrhop(:,:,jpk) = 0._wp
+#if ! defined key_diaar5
+         CALL iom_put( 'rhop', zrhop )
+#else
+         ! If key_diaar5 set then there is already an iom_put call to output rhop.
+         ! Really should be a standard diagnostics option?
+#endif
+      ENDIF
+
+      IF( iom_use("ut") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "ut", z3d )                  ! product of temperature and zonal velocity at U points
+      ENDIF
+
+      IF( iom_use("vt") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "vt", z3d )                  ! product of temperature and meridional velocity at V points
+      ENDIF
+
+      IF( iom_use("wt") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z3d(ji,jj,1) = wn(ji,jj,1) * tsn(ji,jj,1,jp_tem)
+            END DO
+         END DO
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = wn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk-1,jp_tem) + tsn(ji,jj,jk,jp_tem) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "wt", z3d )                  ! product of temperature and vertical velocity at W points
+      ENDIF
+
+      IF( iom_use("us") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "us", z3d )                  ! product of salinity and zonal velocity at U points
+      ENDIF
+
+      IF( iom_use("vs") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "vs", z3d )                  ! product of salinity and meridional velocity at V points
+      ENDIF
+
+      IF( iom_use("ws") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z3d(ji,jj,1) = wn(ji,jj,1) * tsn(ji,jj,1,jp_sal)
+            END DO
+         END DO
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = wn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk-1,jp_sal) + tsn(ji,jj,jk,jp_sal) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "ws", z3d )                  ! product of salinity and vertical velocity at W points
+      ENDIF
+
+      IF( iom_use("urhop") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = un(ji,jj,jk) * 0.5 * ( zrhop(ji,jj,jk) + zrhop(ji+1,jj,jk) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "urhop", z3d )                  ! product of density and zonal velocity at U points
+      ENDIF
+
+      IF( iom_use("vrhop") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = vn(ji,jj,jk) * 0.5 * ( zrhop(ji,jj,jk) + zrhop(ji,jj+1,jk) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "vrhop", z3d )                  ! product of density and meridional velocity at V points
+      ENDIF
+
+      IF( iom_use("wrhop") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z3d(ji,jj,1) = wn(ji,jj,1) * zrhop(ji,jj,1)
+            END DO
+         END DO
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = wn(ji,jj,jk) * 0.5 * ( zrhop(ji,jj,jk-1) + zrhop(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "wrhop", z3d )                  ! product of density and vertical velocity at W points
+      ENDIF
+
+      !
+      CALL wrk_dealloc( jpi , jpj      , z2d )
+      CALL wrk_dealloc( jpi , jpj, jpk , z3d )
+      CALL wrk_dealloc( jpi , jpj, jpk , zrhop )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dia_prod')
+      !
+   END SUBROUTINE dia_prod
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :                                         NO diaprod
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER :: lk_diaprod = .FALSE.   ! coupled flag
+CONTAINS
+   SUBROUTINE dia_prod( kt )   ! Empty routine
+      INTEGER ::   kt
+      WRITE(*,*) 'dia_prod: You should not have seen this print! error?', kt
+   END SUBROUTINE dia_prod
+#endif
+   !!======================================================================
+END MODULE diaprod
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 6491)
@@ -39,4 +39,5 @@
    USE zdfmxl          ! mixed layer
    USE dianam          ! build name of file (routine)
+   USE zdftke          ! vertical physics: one-equation scheme 
    USE zdfddm          ! vertical  physics: double diffusion
    USE diahth          ! thermocline diagnostics
@@ -242,4 +243,8 @@
       CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef.
       CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef.
+      IF( lk_zdftke ) THEN   
+         CALL iom_put( "tke"      , en                               )    ! TKE budget: Turbulent Kinetic Energy   
+         CALL iom_put( "tke_niw"  , e_niw                            )    ! TKE budget: Near-inertial waves   
+      ENDIF 
       CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm)
 
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 6491)
@@ -36,5 +36,5 @@
    PUBLIC clo_bat      ! routine called in domzgr module
 
-   INTEGER, PUBLIC, PARAMETER          ::   jpncs   = 4      !: number of closed sea
+   INTEGER, PUBLIC, PARAMETER          ::   jpncs   = 10      !: number of closed sea
    INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncstt            !: Type of closed sea
    INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi1, ncsj1     !: south-west closed sea limits (i,j)
@@ -155,16 +155,56 @@
             ncsi2(4)   = 76  ;  ncsj2(4)   = 61
             ncsir(4,1) = 84  ;  ncsjr(4,1) = 59 
-            !                                        ! =======================
-         CASE ( 025 )                                ! ORCA_R025 configuration
-            !                                        ! =======================
-            ncsnr(1)   = 1    ; ncstt(1)   = 0               ! Caspian + Aral sea
-            ncsi1(1)   = 1330 ; ncsj1(1)   = 645
-            ncsi2(1)   = 1400 ; ncsj2(1)   = 795
+            !                                        ! ================================
+         CASE ( 025 )                                ! ORCA_R025 extended configuration
+            !                                        ! ================================
+            ncsnr(1)   = 1    ; ncstt(1)   = 0               ! Caspian sea
+            ncsi1(1)   = 1330 ; ncsj1(1)   = 831
+            ncsi2(1)   = 1375 ; ncsj2(1)   = 981
             ncsir(1,1) = 1    ; ncsjr(1,1) = 1
             !                                        
-            ncsnr(2)   = 1    ; ncstt(2)   = 0               ! Azov Sea 
-            ncsi1(2)   = 1284 ; ncsj1(2)   = 722
-            ncsi2(2)   = 1304 ; ncsj2(2)   = 747
+            ncsnr(2)   = 1    ; ncstt(2)   = 0               ! Aral sea
+            ncsi1(2)   = 1376 ; ncsj1(2)   = 900
+            ncsi2(2)   = 1400 ; ncsj2(2)   = 981
             ncsir(2,1) = 1    ; ncsjr(2,1) = 1
+            !                                        
+            ncsnr(3)   = 1    ; ncstt(3)   = 0               ! Azov Sea 
+            ncsi1(3)   = 1284 ; ncsj1(3)   = 908
+            ncsi2(3)   = 1304 ; ncsj2(3)   = 933
+            ncsir(3,1) = 1    ; ncsjr(3,1) = 1
+            !
+            ncsnr(4)   = 1    ; ncstt(4)   = 0               ! Lake Superior  
+            ncsi1(4)   = 781  ; ncsj1(4)   = 904 
+            ncsi2(4)   = 815  ; ncsj2(4)   = 926 
+            ncsir(4,1) = 1    ; ncsjr(4,1) = 1 
+            ! 
+            ncsnr(5)   = 1    ; ncstt(5)   = 0               ! Lake Michigan
+            ncsi1(5)   = 795  ; ncsj1(5)   = 871             
+            ncsi2(5)   = 813  ; ncsj2(5)   = 905 
+            ncsir(5,1) = 1    ; ncsjr(5,1) = 1 
+            ! 
+            ncsnr(6)   = 1    ; ncstt(6)   = 0               ! Lake Huron part 1
+            ncsi1(6)   = 814  ; ncsj1(6)   = 882             
+            ncsi2(6)   = 825  ; ncsj2(6)   = 905 
+            ncsir(6,1) = 1    ; ncsjr(6,1) = 1 
+            ! 
+            ncsnr(7)   = 1    ; ncstt(7)   = 0               ! Lake Huron part 2
+            ncsi1(7)   = 826  ; ncsj1(7)   = 889             
+            ncsi2(7)   = 833  ; ncsj2(7)   = 905 
+            ncsir(7,1) = 1    ; ncsjr(7,1) = 1 
+            ! 
+            ncsnr(8)   = 1    ; ncstt(8)   = 0               ! Lake Erie
+            ncsi1(8)   = 816  ; ncsj1(8)   = 871             
+            ncsi2(8)   = 837  ; ncsj2(8)   = 881 
+            ncsir(8,1) = 1    ; ncsjr(8,1) = 1 
+            ! 
+            ncsnr(9)   = 1    ; ncstt(9)   = 0               ! Lake Ontario
+            ncsi1(9)   = 831  ; ncsj1(9)   = 882             
+            ncsi2(9)   = 847  ; ncsj2(9)   = 889 
+            ncsir(9,1) = 1    ; ncsjr(9,1) = 1 
+            ! 
+            ncsnr(10)   = 1    ; ncstt(10)   = 0               ! Lake Victoria  
+            ncsi1(10)   = 1274 ; ncsj1(10)   = 672 
+            ncsi2(10)   = 1289 ; ncsj2(10)   = 687 
+            ncsir(10,1) = 1    ; ncsjr(10,1) = 1 
             !
          END SELECT
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6491)
@@ -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_rstart , ln_rstdate, 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
@@ -173,5 +173,6 @@
          WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
          WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
-         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
+         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/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 6491)
@@ -136,4 +136,5 @@
       INTEGER  ::   isrow                    ! index for ORCA1 starting row
       INTEGER , POINTER, DIMENSION(:,:) ::  imsk
+      REAL(wp) ::  zphi_drake_passage, zshlat_antarc
       REAL(wp), POINTER, DIMENSION(:,:) ::  zwf
       !!
@@ -445,4 +446,27 @@
       ENDIF
       !
+      IF( cp_cfg == "orca" .AND. jp_cfg == 025 .AND. rn_shlat == 0.0 ) THEN   
+         !                                              ! ORCA_R025 configuration
+         !                                              ! Increased lateral friction on parts of Antarctic coastline
+         !                                              ! for increased stability
+         !                                              ! NB. This only works to do this here if we have free slip 
+         !                                              ! generally, so fmask is zero at coast points.
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   orca_r025: increase friction in following regions : '
+         IF(lwp) WRITE(numout,*) '      whole Antarctic coastline: partial slip shlat=1 '
+
+         zphi_drake_passage = -58.0_wp
+         zshlat_antarc = 1.0_wp
+         zwf(:,:) = fmask(:,:,1)         
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               IF( gphif(ji,jj) .lt. zphi_drake_passage .and. fmask(ji,jj,1) == 0._wp ) THEN
+                  fmask(ji,jj,:) = zshlat_antarc * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
+                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
+               ENDIF
+            END DO
+         END DO
+      END IF
+      !
       CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
 
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 6491)
@@ -594,4 +594,15 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
       ENDIF
+
+      ! Write outputs
+      ! =============
+      CALL iom_put(     "e3t" , fse3t_n  (:,:,:) )
+      CALL iom_put(     "e3u" , fse3u_n  (:,:,:) )
+      CALL iom_put(     "e3v" , fse3v_n  (:,:,:) )
+      CALL iom_put(     "e3w" , fse3w_n  (:,:,:) )
+      CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) )
+      IF( iom_use("e3tdef") )   &
+         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
+
       !
       ! Time filter and swap of scale factors
@@ -666,13 +677,4 @@
       END DO
 
-      ! Write outputs
-      ! =============
-      CALL iom_put(     "e3t" , fse3t_n  (:,:,:) )
-      CALL iom_put(     "e3u" , fse3u_n  (:,:,:) )
-      CALL iom_put(     "e3v" , fse3v_n  (:,:,:) )
-      CALL iom_put(     "e3w" , fse3w_n  (:,:,:) )
-      CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) )
-      IF( iom_use("e3tdef") )   &
-         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
 
       ! write restart file
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90	(revision 6491)
@@ -371,4 +371,6 @@
       IF( .NOT. ln_bergdia )   RETURN            !!gm useless iom will control whether it is output or not
       !
+      CALL iom_put( "berg_total_melt"      , berg_grid%floating_melt(:,:)   )   ! Total melt flux to ocean      [kg/m2/s]
+      CALL iom_put( "berg_total_heat_flux" , berg_grid%calving_hflx(:,:)    )   ! Total iceberg-ocean heat flux [W/m2]
       CALL iom_put( "berg_melt"        , berg_melt   (:,:)   )   ! Melt rate of icebergs                     [kg/m2/s]
       CALL iom_put( "berg_buoy_melt"   , buoy_melt   (:,:)   )   ! Buoyancy component of iceberg melt rate   [kg/m2/s]
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 6491)
@@ -12,4 +12,6 @@
    !!            -    !                            Currently needs a fixed processor
    !!            -    !                            layout between restarts
+   !!            -    !  2015-11  Dave Storkey     Convert icb_rst_read to use IOM so can
+   !!                                              read single restart files
    !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
@@ -18,8 +20,11 @@
    !!----------------------------------------------------------------------
    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
+   USE iom
    USE icb_oce        ! define iceberg arrays
    USE icbutl         ! iceberg utility routines
@@ -57,10 +62,6 @@
       INTEGER                      ::   idim, ivar, iatt
       INTEGER                      ::   jn, iunlim_dim, ibergs_in_file
-      INTEGER                      ::   iclass
-      INTEGER, DIMENSION(1)        ::   istrt, ilngth, idata
-      INTEGER, DIMENSION(2)        ::   istrt2, ilngth2
-      INTEGER, DIMENSION(nkounts)  ::   idata2
-      REAL(wp), DIMENSION(1)       ::   zdata                                         ! need 1d array to read in with
-                                                                                            ! start and count arrays
+      INTEGER                      ::   ii,ij,iclass
+      REAL(wp), DIMENSION(nkounts) ::   zdata      
       LOGICAL                      ::   ll_found_restart
       CHARACTER(len=256)           ::   cl_path
@@ -71,138 +72,62 @@
       !!----------------------------------------------------------------------
 
-      ! Find a restart file. Assume iceberg restarts in same directory as ocean restarts. 
+      ! Find a restart file. Assume iceberg restarts in same directory as ocean restarts
+      ! and are called TRIM(cn_ocerst)//'_icebergs'
       cl_path = TRIM(cn_ocerst_indir)
       IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/'
-      cl_filename = ' '
-      IF ( lk_mpp ) THEN
-         cl_filename = ' '
-         WRITE( cl_filename, '("restart_icebergs_",I4.4,".nc")' ) narea-1
-         INQUIRE( file=TRIM(cl_path)//TRIM(cl_filename), exist=ll_found_restart )
-      ELSE
-         cl_filename = 'restart_icebergs.nc'
-         INQUIRE( file=TRIM(cl_path)//TRIM(cl_filename), exist=ll_found_restart )
-      ENDIF
-
-      IF ( .NOT. ll_found_restart) THEN                     ! only do the following if a file was found
-         CALL ctl_stop('icebergs: no restart file found')
-      ENDIF
-
-      IF (nn_verbose_level >= 0 .AND. lwp)  &
-         WRITE(numout,'(2a)') 'icebergs, read_restart_bergs: found restart file = ',TRIM(cl_path)//TRIM(cl_filename)
-
-      nret = NF90_OPEN(TRIM(cl_path)//TRIM(cl_filename), NF90_NOWRITE, ncid)
-      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_open failed')
-
-      nret = nf90_inquire(ncid, idim, ivar, iatt, iunlim_dim)
-      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_inquire failed')
-
-      IF( iunlim_dim .NE. -1) THEN
-
-         nret = nf90_inquire_dimension(ncid, iunlim_dim, cl_dname, ibergs_in_file)
-         IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_inq_dimlen failed')
-
-         nret = NF90_INQ_VARID(ncid, 'number', numberid)
-         nret = NF90_INQ_VARID(ncid, 'mass_scaling', nscaling_id)
-         nret = NF90_INQ_VARID(ncid, 'xi', nxid)
-         nret = NF90_INQ_VARID(ncid, 'yj', nyid)
-         nret = NF90_INQ_VARID(ncid, 'lon', nlonid)
-         nret = NF90_INQ_VARID(ncid, 'lat', nlatid)
-         nret = NF90_INQ_VARID(ncid, 'uvel', nuvelid)
-         nret = NF90_INQ_VARID(ncid, 'vvel', nvvelid)
-         nret = NF90_INQ_VARID(ncid, 'mass', nmassid)
-         nret = NF90_INQ_VARID(ncid, 'thickness', nthicknessid)
-         nret = NF90_INQ_VARID(ncid, 'width', nwidthid)
-         nret = NF90_INQ_VARID(ncid, 'length', nlengthid)
-         nret = NF90_INQ_VARID(ncid, 'year', nyearid)
-         nret = NF90_INQ_VARID(ncid, 'day', ndayid)
-         nret = NF90_INQ_VARID(ncid, 'mass_of_bits', nmass_of_bits_id)
-         nret = NF90_INQ_VARID(ncid, 'heat_density', nheat_density_id)
-
-         ilngth(1) = 1
-         istrt2(1) = 1
-         ilngth2(1) = nkounts
-         ilngth2(2) = 1
-         DO jn=1, ibergs_in_file
-
-            istrt(1) = jn
-            istrt2(2) = jn
-
-            nret = NF90_GET_VAR(ncid, numberid, idata2, istrt2, ilngth2 )
-            localberg%number(:) = idata2(:)
-
-            nret = NF90_GET_VAR(ncid, nscaling_id, zdata, istrt, ilngth )
-            localberg%mass_scaling = zdata(1)
-
-            nret = NF90_GET_VAR(ncid, nlonid, zdata, istrt, ilngth)
-            localpt%lon = zdata(1)
-            nret = NF90_GET_VAR(ncid, nlatid, zdata, istrt, ilngth)
-            localpt%lat = zdata(1)
-            IF (nn_verbose_level >= 2 .AND. lwp) THEN
-               WRITE(numout,'(a,i5,a,2f10.4,a,i5)') 'icebergs, read_restart_bergs: berg ',jn,' is at ', &
-                                              localpt%lon,localpt%lat,' on PE ',narea-1
+      cl_filename = TRIM(cn_ocerst_in)//'_icebergs'
+      CALL iom_open( TRIM(cl_path)//cl_filename, ncid )
+
+      IF( iom_file(ncid)%iduld .GE. 0) THEN
+
+         ibergs_in_file = iom_file(ncid)%lenuld
+         DO jn = 1,ibergs_in_file
+
+            ! iom_get treats the unlimited dimension as time. Here the unlimited dimension 
+            ! is the iceberg index, but we can still use the ktime keyword to get the iceberg we want. 
+
+            CALL iom_get( ncid, 'xi'     ,localpt%xi  , ktime=jn )
+            CALL iom_get( ncid, 'yj'     ,localpt%yj  , ktime=jn )
+
+            ii = INT( localpt%xi + 0.5 )
+            ij = INT( localpt%yj + 0.5 )
+            ! Only proceed if this iceberg is on the local processor (excluding halos).
+            IF ( ii .GE. nldi+nimpp-1 .AND. ii .LE. nlei+nimpp-1 .AND. &
+           &     ij .GE. nldj+njmpp-1 .AND. ij .LE. nlej+njmpp-1 ) THEN           
+
+               CALL iom_get( ncid, jpdom_unknown, 'number'       , zdata(:) , ktime=jn, kstart=(/1/), kcount=(/nkounts/) )
+               localberg%number(:) = INT(zdata(:))
+               CALL iom_get( ncid, 'mass_scaling' , localberg%mass_scaling, ktime=jn )
+               CALL iom_get( ncid, 'lon'          , localpt%lon           , ktime=jn )
+               CALL iom_get( ncid, 'lat'          , localpt%lat           , ktime=jn )
+               CALL iom_get( ncid, 'uvel'         , localpt%uvel          , ktime=jn )
+               CALL iom_get( ncid, 'vvel'         , localpt%vvel          , ktime=jn )
+               CALL iom_get( ncid, 'mass'         , localpt%mass          , ktime=jn )
+               CALL iom_get( ncid, 'thickness'    , localpt%thickness     , ktime=jn )
+               CALL iom_get( ncid, 'width'        , localpt%width         , ktime=jn )
+               CALL iom_get( ncid, 'length'       , localpt%length        , ktime=jn )
+               CALL iom_get( ncid, 'year'         , zdata(1)              , ktime=jn )
+               localpt%year = INT(zdata(1))
+               CALL iom_get( ncid, 'day'          , localpt%day           , ktime=jn )
+               CALL iom_get( ncid, 'mass_of_bits' , localpt%mass_of_bits  , ktime=jn )
+               CALL iom_get( ncid, 'heat_density' , localpt%heat_density  , ktime=jn )
+
+               !
+               CALL icb_utl_add( localberg, localpt )
+
             ENDIF
-            nret = NF90_GET_VAR(ncid, nxid, zdata, istrt, ilngth)
-            localpt%xi = zdata(1)
-            nret = NF90_GET_VAR(ncid, nyid, zdata, istrt, ilngth)
-            localpt%yj = zdata(1)
-            nret = NF90_GET_VAR(ncid, nuvelid, zdata, istrt, ilngth )
-            localpt%uvel = zdata(1)
-            nret = NF90_GET_VAR(ncid, nvvelid, zdata, istrt, ilngth )
-            localpt%vvel = zdata(1)
-            nret = NF90_GET_VAR(ncid, nmassid, zdata, istrt, ilngth )
-            localpt%mass = zdata(1)
-            nret = NF90_GET_VAR(ncid, nthicknessid, zdata, istrt, ilngth )
-            localpt%thickness = zdata(1)
-            nret = NF90_GET_VAR(ncid, nwidthid, zdata, istrt, ilngth )
-            localpt%width = zdata(1)
-            nret = NF90_GET_VAR(ncid, nlengthid, zdata, istrt, ilngth )
-            localpt%length = zdata(1)
-            nret = NF90_GET_VAR(ncid, nyearid, idata, istrt, ilngth )
-            localpt%year = idata(1)
-            nret = NF90_GET_VAR(ncid, ndayid, zdata, istrt, ilngth )
-            localpt%day = zdata(1)
-            nret = NF90_GET_VAR(ncid, nmass_of_bits_id, zdata, istrt, ilngth )
-            localpt%mass_of_bits = zdata(1)
-            nret = NF90_GET_VAR(ncid, nheat_density_id, zdata, istrt, ilngth )
-            localpt%heat_density = zdata(1)
-            !
-            CALL icb_utl_add( localberg, localpt )
+
          END DO
-         !
-      ENDIF
-
-      nret = NF90_INQ_DIMID( ncid, 'c', nc_dim )
-      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_inq_dimid c failed')
-
-      nret = NF90_INQUIRE_DIMENSION( ncid, nc_dim, cl_dname, iclass )
-      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_inquire_dimension failed')
-
-      nret = NF90_INQ_VARID(ncid, 'kount'       , nkountid)
-      nret = NF90_INQ_VARID(ncid, 'calving'     , ncalvid)
-      nret = NF90_INQ_VARID(ncid, 'calving_hflx', ncalvhid)
-      nret = NF90_INQ_VARID(ncid, 'stored_ice'  , nsiceid)
-      nret = NF90_INQ_VARID(ncid, 'stored_heat' , nsheatid)
-
-      nstrt3(1) = 1
-      nstrt3(2) = 1
-      nlngth3(1) = jpi
-      nlngth3(2) = jpj
-      nlngth3(3) = 1
-
-      DO jn = 1, iclass
-         nstrt3(3) = jn
-         nret      = NF90_GET_VAR( ncid, nsiceid , griddata, nstrt3, nlngth3 )
-         berg_grid%stored_ice(:,:,jn) = griddata(:,:,1)
-      END DO
-
-      nret = NF90_GET_VAR( ncid, ncalvid , src_calving          (:,:) )
-      nret = NF90_GET_VAR( ncid, ncalvhid, src_calving_hflx     (:,:) )
-      nret = NF90_GET_VAR( ncid, nsheatid, berg_grid%stored_heat(:,:) )
-      nret = NF90_GET_VAR( ncid, nkountid, idata2(:) )
-      num_bergs(:) = idata2(:)
-
-      ! Finish up
-      nret = NF90_CLOSE(ncid)
-      IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_close failed')
+
+      ENDIF 
+
+      ! Gridded variables
+      CALL iom_get( ncid, jpdom_autoglo,    'calving'     , src_calving  )
+      CALL iom_get( ncid, jpdom_autoglo,    'calving_hflx', src_calving_hflx  )
+      CALL iom_get( ncid, jpdom_autoglo,    'stored_heat' , berg_grid%stored_heat  )
+      CALL iom_get( ncid, jpdom_autoglo_xy, 'stored_ice'  , berg_grid%stored_ice, kstart=(/1,1,1/), kcount=(/1,1,nclasses/) )
+      
+      CALL iom_get( ncid, jpdom_unknown, 'kount' , zdata(:) )
+      num_bergs(:) = INT(zdata(:))
 
       ! Sanity check
@@ -211,5 +136,6 @@
          WRITE(numout,'(2(a,i5))') 'icebergs, read_restart_bergs: # bergs =',jn,' on PE',narea-1
       IF( lk_mpp ) THEN
-         CALL mpp_sum(ibergs_in_file)
+         ! Only mpp_sum ibergs_in_file if we are reading from multiple restart files. 
+         IF( INDEX(iom_file(ncid)%name,'icebergs.nc' ) .EQ. 0 ) CALL mpp_sum(ibergs_in_file)
          CALL mpp_sum(jn)
       ENDIF
@@ -217,4 +143,7 @@
          &                                    ' bergs in the restart file and', jn,' bergs have been read'
       !
+      ! Finish up
+      CALL iom_close( ncid )
+      !
       IF( lwp .and. nn_verbose_level >= 0)  WRITE(numout,'(a)') 'icebergs, read_restart_bergs: completed'
       !
@@ -231,6 +160,9 @@
       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
+      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 +172,16 @@
       cl_path = TRIM(cn_ocerst_outdir)
       IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/'
+      IF ( ln_rstdate ) THEN
+         CALL ju2ymds( fjulday + rdttra(1) / rday, 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/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90	(revision 6491)
@@ -18,4 +18,5 @@
    USE dom_oce        ! NEMO domain
    USE in_out_manager ! NEMO IO routines, numout in particular
+   USE iom
    USE lib_mpp        ! NEMO MPI routines, ctl_stop in particular
    USE phycst         ! NEMO physical constants
@@ -160,5 +161,6 @@
             zmelt    = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt   ! kg/s
             berg_grid%floating_melt(ii,ij) = berg_grid%floating_melt(ii,ij) + zmelt    * z1_e1e2    ! kg/m2/s
-            zheat = zmelt * pt%heat_density              ! kg/s x J/kg = J/s
+!            zheat = zmelt * pt%heat_density              ! kg/s x J/kg = J/s
+            zheat = zmelt * lfus                           !rma kg/s x J/kg (latent heat of fusion) = J/s
             berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + zheat    * z1_e1e2    ! W/m2
             CALL icb_dia_melt( ii, ij, zMnew, zheat, this%mass_scaling,       &
@@ -208,5 +210,5 @@
       IF(.NOT. ln_passive_mode ) THEN
          emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:)
-!!       qns (:,:) = qns (:,:) + berg_grid%calving_hflx (:,:)  !!gm heat flux not yet properly coded ==>> need it, SOLVE that!
+         qns (:,:) = qns (:,:) - berg_grid%calving_hflx (:,:)  
       ENDIF
       !
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 6491)
@@ -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/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 6491)
@@ -710,4 +710,5 @@
       CHARACTER(LEN=256)             ::   clname      ! file name
       CHARACTER(LEN=1)               ::   clrankpv, cldmspc      ! 
+      LOGICAL                        ::   ll_depth_spec ! T => if kstart, kcount present then *only* use values for 3rd spatial dimension.
       !---------------------------------------------------------------------
       !
@@ -722,5 +723,6 @@
       IF( PRESENT(kcount) .AND. (.NOT. PRESENT(kstart)) ) CALL ctl_stop(trim(clinfo), 'kcount present needs kstart present')
       IF( PRESENT(kstart) .AND. (.NOT. PRESENT(kcount)) ) CALL ctl_stop(trim(clinfo), 'kstart present needs kcount present')
-      IF( PRESENT(kstart) .AND. idom /= jpdom_unknown   ) CALL ctl_stop(trim(clinfo), 'kstart present needs kdom = jpdom_unknown')
+      IF( PRESENT(kstart) .AND. idom /= jpdom_unknown .AND.  idom /= jpdom_autoglo_xy  ) &
+     &           CALL ctl_stop(trim(clinfo), 'kstart present needs kdom = jpdom_unknown or kdom = jpdom_autoglo_xy')
 
       luse_jattr = .false.
@@ -755,4 +757,10 @@
          ! update idom definition...
          ! Identify the domain in case of jpdom_auto(glo/dta) definition
+         IF( idom == jpdom_autoglo_xy ) THEN
+            ll_depth_spec = .TRUE.
+            idom = jpdom_autoglo
+         ELSE
+            ll_depth_spec = .FALSE.
+         ENDIF
          IF( idom == jpdom_autoglo .OR. idom == jpdom_autodta ) THEN            
             IF( idom == jpdom_autoglo ) THEN   ;   idom = jpdom_global 
@@ -808,7 +816,7 @@
          istart(idmspc+1) = itime
 
-         IF(              PRESENT(kstart)       ) THEN ; istart(1:idmspc) = kstart(1:idmspc) ; icnt(1:idmspc) = kcount(1:idmspc)
+         IF( PRESENT(kstart) .AND. .NOT. ll_depth_spec ) THEN ; istart(1:idmspc) = kstart(1:idmspc) ; icnt(1:idmspc) = kcount(1:idmspc)
          ELSE
-            IF(           idom == jpdom_unknown ) THEN                                       ; icnt(1:idmspc) = idimsz(1:idmspc)
+            IF(           idom == jpdom_unknown ) THEN                                                ; icnt(1:idmspc) = idimsz(1:idmspc)
             ELSE 
                IF( .NOT. PRESENT(pv_r1d) ) THEN   !   not a 1D array
@@ -833,6 +841,7 @@
                   ENDIF
                   IF( PRESENT(pv_r3d) ) THEN
-                     IF( idom == jpdom_data ) THEN   ; icnt(3) = jpkdta
-                     ELSE                            ; icnt(3) = jpk
+                     IF( idom == jpdom_data ) THEN                                  ; icnt(3) = jpkdta
+                     ELSE IF( ll_depth_spec .AND. PRESENT(kstart) ) THEN            ; istart(3) = kstart(3); icnt(3) = kcount(3)
+                     ELSE                                                           ; icnt(3) = jpk
                      ENDIF
                   ENDIF
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90	(revision 6491)
@@ -26,5 +26,6 @@
    INTEGER, PARAMETER, PUBLIC ::   jpdom_unknown       = 7   !: No dimension checking
    INTEGER, PARAMETER, PUBLIC ::   jpdom_autoglo       = 8   !: 
-   INTEGER, PARAMETER, PUBLIC ::   jpdom_autodta       = 9   !: 
+   INTEGER, PARAMETER, PUBLIC ::   jpdom_autoglo_xy    = 9   !: Automatically set horizontal dimensions only
+   INTEGER, PARAMETER, PUBLIC ::   jpdom_autodta       = 10  !: 
 
    INTEGER, PARAMETER, PUBLIC ::   jpioipsl    = 100      !: Use ioipsl (fliocom only) library
@@ -57,4 +58,5 @@
       INTEGER                                   ::   nvars    !: number of identified varibles in the file
       INTEGER                                   ::   iduld    !: id of the unlimited dimension
+      INTEGER                                   ::   lenuld   !: length of the unlimited dimension (number of records in file)
       INTEGER                                   ::   irec     !: writing record position  
       CHARACTER(LEN=32)                         ::   uldname  !: name of the unlimited dimension
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90	(revision 6491)
@@ -154,6 +154,7 @@
          CALL iom_nf90_check(NF90_Inquire(if90id, unlimitedDimId = iom_file(kiomid)%iduld), clinfo)
          IF ( iom_file(kiomid)%iduld .GE. 0 ) THEN
-           CALL iom_nf90_check(NF90_Inquire_Dimension(if90id, iom_file(kiomid)%iduld,   &
-        &                                               name = iom_file(kiomid)%uldname), clinfo)
+           CALL iom_nf90_check(NF90_Inquire_Dimension(if90id, iom_file(kiomid)%iduld,     & 
+        &                                               name = iom_file(kiomid)%uldname,  &
+        &                                               len  = iom_file(kiomid)%lenuld ), clinfo )
          ENDIF
          IF(lwp) WRITE(numout,*) '                   ---> '//TRIM(cdname)//' OK'
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 6491)
@@ -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
@@ -55,8 +56,10 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! ocean time-step
+      INTEGER             ::   iyear, imonth, iday
+      REAL (wp)           ::   zsec
       !!
       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
       !!----------------------------------------------------------------------
       !
@@ -82,7 +85,12 @@
       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
+               CALL ju2ymds( fjulday + rdttra(1) / rday, 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/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6491)
@@ -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
@@ -33,7 +34,16 @@
    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, SAVE, DIMENSION(:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m] 
+   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
+
+   ! Namelist variables for  namzdf_mldzint
+   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 
 
    !! * Substitutions
@@ -52,5 +62,6 @@
       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),       &
+        &          ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc )
          !
          IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc )
@@ -90,5 +101,5 @@
       CALL wrk_alloc( jpi,jpj, imld )
 
-      IF( kt == nit000 ) THEN
+      IF( kt <= nit000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth'
@@ -134,8 +145,16 @@
       END DO
       IF( .NOT.lk_offline ) THEN            ! no need to output in offline mode
+      IF( kt >= nit000 ) THEN               ! workaround for calls before SOMETHING reads the XIOS namelist
          CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth
          CALL iom_put( "mldkz5"  , hmld )   ! turbocline depth
       ENDIF
+      ENDIF
       
+      ! Vertically-interpolated mixed-layer depth diagnostic
+      IF( iom_use( "mldzint" ) ) THEN
+         CALL zdf_mxl_zint( kt )
+         CALL iom_put( "mldzint" , hmld_zint )
+      ENDIF
+
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 )
       !
@@ -145,4 +164,189 @@
       !
    END SUBROUTINE zdf_mxl
+
+   SUBROUTINE zdf_mxl_zint( kt ) 
+      !!---------------------------------------------------------------------------------- 
+      !!                    ***  ROUTINE zdf_mxl_zint  *** 
+      !                                                                        
+      !   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                                      
+      ! 
+      !!----------------------------------------------------------------------------------- 
+
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
+      ! Local variables
+      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 
+      INTEGER :: ios
+
+      NAMELIST/namzdf_mldzint/ nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac
+
+      !!------------------------------------------------------------------------------------- 
+      !  
+      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 ) 
+ 
+      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 )
+
+         WRITE(numout,*) '===== Vertically-interpolated mixed layer ====='
+         WRITE(numout,*) 'nn_mld_type = ',nn_mld_type
+         WRITE(numout,*) 'rn_zref = ',rn_zref
+         WRITE(numout,*) 'rn_dT_crit = ',rn_dT_crit
+         WRITE(numout,*) 'rn_iso_frac = ',rn_iso_frac
+         WRITE(numout,*) '==============================================='
+      ENDIF
+ 
+      ! Set the mixed layer depth criterion at each grid point 
+      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 
+
+      ! 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(:,:) 
+      ll_found(:,:) = .false. 
+      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 
+
+      ! 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
 
    !!======================================================================
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 6491)
@@ -83,4 +83,5 @@
    INTEGER  ::   nn_htau   ! type of tke profile of penetration (=0/1)
    REAL(wp) ::   rn_efr    ! fraction of TKE surface value which penetrates in the ocean
+   REAL(wp) ::   rn_c      ! fraction of TKE added within the mixed layer by nn_etau
    LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
    REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells
@@ -88,8 +89,11 @@
    REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
    REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m]
+   REAL(wp) ::   rhtau                     ! coefficient to relate MLD to htau when nn_htau == 2
    REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3)
    REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
 
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
 #if defined key_c1d
@@ -114,4 +118,5 @@
       !!----------------------------------------------------------------------
       ALLOCATE(                                                                    &
+         &      efr  (jpi,jpj)     , e_niw(jpi,jpj,jpk) ,                         &      
 #if defined key_c1d
          &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          &
@@ -420,7 +425,22 @@
       END DO
 
+      !                                 ! Save TKE prior to nn_etau addition  
+      e_niw(:,:,:) = en(:,:,:)  
+      !  
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  TKE due to surface and internal wave breaking
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( nn_htau == 2 ) THEN           !* mixed-layer depth dependant length scale
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               htau(ji,jj) = rhtau * hmlp(ji,jj)
+            END DO
+         END DO
+      ENDIF
+#if defined key_iomput
+      !
+      CALL iom_put( "htau", htau(:,:) )  ! Check htau (even if constant in time)
+#endif
+      !
       IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
          DO jk = 2, jpkm1
@@ -457,6 +477,32 @@
             END DO
          END DO
+      ELSEIF( nn_etau == 4 ) THEN       !* column integral independant of htau (rn_efr must be scaled up)
+         IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau 
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  efr(ji,jj) = rn_efr / ( htau(ji,jj) * ( 1._wp - EXP( -bathy(ji,jj) / htau(ji,jj) ) ) )
+               END DO
+            END DO
+         ENDIF
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  en(ji,jj,jk) = en(ji,jj,jk) + efr(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                     &                                                   * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
       ENDIF
       CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      !
+      DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term  
+         DO jj = 2, jpjm1  
+            DO ji = fs_2, fs_jpim1   ! vector opt.  
+               e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk)  
+            END DO  
+         END DO  
+      END DO  
+      !  
+      CALL lbc_lnk( e_niw, 'W', 1. )  
       !
       CALL wrk_dealloc( jpi,jpj, imlc )    ! integer
@@ -722,7 +768,7 @@
          &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
          &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
-         &                 nn_etau , nn_htau  , rn_efr   
-      !!----------------------------------------------------------------------
-      !
+         &                 nn_etau , nn_htau  , rn_efr , rn_c   
+      !!----------------------------------------------------------------------
+
       REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy
       READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
@@ -757,4 +803,5 @@
          WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau
          WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr
+         WRITE(numout,*) '      fraction of TKE added within the mixed layer by nn_etau rn_c    = ', rn_c
          WRITE(numout,*)
          WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
@@ -767,5 +814,5 @@
       IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
       IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
-      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
+      IF( nn_htau < 0  .OR.  nn_htau > 5 )   CALL ctl_stop( 'bad flag: nn_htau is 0 to 5    ' )
       IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
 
@@ -780,6 +827,12 @@
       ENDIF
 
+      IF( nn_etau /= 0 .and. nn_htau == 2 ) THEN
+          ierr = zdf_mxl_alloc()
+          nmln(:,:) = nlb10           ! Initialization of nmln
+      ENDIF
+
       !                               !* depth of penetration of surface tke
       IF( nn_etau /= 0 ) THEN      
+         htau(:,:) = 0._wp
          SELECT CASE( nn_htau )             ! Choice of the depth of penetration
          CASE( 0 )                                 ! constant depth penetration (here 10 meters)
@@ -787,5 +840,37 @@
          CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
             htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )            
+         CASE( 2 )                                 ! fraction of depth-integrated TKE within mixed-layer
+            rhtau = -1._wp / LOG( 1._wp - rn_c )
+         CASE( 3 )                                 ! F(latitude) : 0.5m to 15m poleward of 20 degrees
+            htau(:,:) = MAX(  0.5_wp, MIN( 15._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )
+         CASE( 4 )                                 ! F(latitude) : 0.5m to 10m/30m poleward of 13/40 degrees north/south
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( gphit(ji,jj) <= 0._wp ) THEN
+                     htau(ji,jj) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
+                  ELSE
+                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
+                  ENDIF
+               END DO
+            END DO
+         CASE ( 5 )                                ! F(latitude) : 0.5m to 10m poleward of 13 degrees north/south,
+            DO jj = 2, jpjm1                       !               10m to 30m between 30/45 degrees south
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( gphit(ji,jj) <= -30._wp ) THEN
+                     htau(ji,jj) = MAX(  10._wp, MIN( 30._wp, 55._wp* ABS( SIN( rpi/120._wp * ( gphit(ji,jj) + 23._wp ) ) ) )   )
+                  ELSE
+                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
+                  ENDIF
+               END DO
+            END DO
          END SELECT
+         !
+         IF( nn_etau == 4 .AND. nn_htau /= 2 ) THEN            ! efr dependant on constant htau
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  efr(ji,jj) = rn_efr / ( htau(ji,jj) * ( 1._wp - EXP( -bathy(ji,jj) / htau(ji,jj) ) ) )
+               END DO
+            END DO
+         ENDIF
       ENDIF
       !                               !* set vertical eddy coef. to the background value
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6491)
@@ -234,4 +234,5 @@
       IF( lk_diaar5  )      CALL dia_ar5( kstp )         ! ar5 diag
       IF( lk_diaharm )      CALL dia_harm( kstp )        ! Tidal harmonic analysis
+                            CALL dia_prod( kstp )        ! ocean model: product diagnostics
                             CALL dia_wri( kstp )         ! ocean model: outputs
       !
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 6488)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 6491)
@@ -88,4 +88,5 @@
 
    USE diawri           ! Standard run outputs             (dia_wri routine)
+   USE diaprod          ! Product diagnostics              (dia_prod routine)
    USE diaptr           ! poleward transports              (dia_ptr routine)
    USE diadct           ! sections transports              (dia_dct routine)
