Index: /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
===================================================================
--- /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 6448)
+++ /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 6448)
@@ -0,0 +1,225 @@
+MODULE iceini
+   !!======================================================================
+   !!                       ***  MODULE iceini   ***
+   !!   Sea-ice model : LIM Sea ice model Initialization
+   !!======================================================================
+   !! History :  3.0  ! 2008-03  (M. Vancoppenolle) LIM-3 original code
+   !!            3.3  ! 2010-12  (G. Madec) add call to lim_thd_init and lim_thd_sal_init
+   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                     LIM sea-ice model
+   !!----------------------------------------------------------------------
+   !!   ice_init      : sea-ice model initialization
+   !!----------------------------------------------------------------------
+   USE phycst         ! physical constants
+   USE dom_oce        ! ocean domain
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE sbc_ice        ! Surface boundary condition: ice   fields
+   USE ice            ! LIM variables
+   USE par_ice        ! LIM parameters
+   USE dom_ice        ! LIM domain
+   USE thd_ice        ! LIM thermodynamical variables
+   USE limitd_me      ! LIM ice thickness distribution
+   USE limmsh         ! LIM mesh
+   USE limistate      ! LIM initial state
+   USE limrst         ! LIM restart
+   USE limthd         ! LIM ice thermodynamics
+   USE limthd_sal     ! LIM ice thermodynamics: salinity
+   USE limvar         ! LIM variables
+   USE limsbc         ! LIM surface boundary condition
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_init   ! called by sbcice_lim.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_init  ***
+      !!
+      !! ** purpose :   Allocate all the dynamic arrays of the LIM-3 modules
+      !!----------------------------------------------------------------------
+      INTEGER :: ierr
+      !!----------------------------------------------------------------------
+
+      !                                ! Allocate the ice arrays
+      ierr =        ice_alloc        ()      ! ice variables
+      ierr = ierr + dom_ice_alloc    ()      ! domain
+      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing
+      ierr = ierr + thd_ice_alloc    ()      ! thermodynamics
+      ierr = ierr + lim_itd_me_alloc ()      ! ice thickness distribution - mechanics
+      !
+      IF( lk_mpp    )   CALL mpp_sum( ierr )
+      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
+      !
+      !                                ! adequation jpk versus ice/snow layers/categories
+      IF( jpl > jpk .OR. (nlay_i+1) > jpk .OR. nlay_s > jpk )   &
+         &      CALL ctl_stop( 'STOP',                     &
+         &     'ice_init: the 3rd dimension of workspace arrays is too small.',   &
+         &     'use more ocean levels or less ice/snow layers/categories.' )
+
+                                       ! Open the reference and configuration namelist files and namelist output file 
+      CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
+      CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
+      IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
+      !
+      CALL ice_run                     ! set some ice run parameters
+      !
+      CALL lim_thd_init                ! set ice thermodynics parameters
+      !
+      CALL lim_thd_sal_init            ! set ice salinity parameters
+      !
+      rdt_ice   = nn_fsbc * rdttra(1)  ! sea-ice timestep
+      r1_rdtice = 1._wp / rdt_ice      ! sea-ice timestep inverse
+      !
+      CALL lim_msh                     ! ice mesh initialization
+      !
+      CALL lim_itd_ini                 ! ice thickness distribution initialization
+      !
+      CALL lim_itd_me_init             ! ice thickness distribution initialization
+      !                                ! Initial sea-ice state
+      IF( .NOT. ln_rstart ) THEN              ! start from rest
+         numit = 0
+         numit = nit000 - 1
+         CALL lim_istate                        ! start from rest: sea-ice deduced from sst
+         CALL lim_var_agg(1)                    ! aggregate category variables in bulk variables
+         CALL lim_var_glo2eqv                   ! convert global variables in equivalent variables
+      ELSE                                   ! start from a restart file
+         CALL lim_rst_read                      ! read the restart file
+         numit = nit000 - 1
+         CALL lim_var_agg(1)                    ! aggregate ice variables
+         CALL lim_var_glo2eqv                   ! convert global var in equivalent variables
+      ENDIF
+      !
+      hi_max(jpl) = 99._wp  ! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl)
+      !
+      CALL lim_sbc_init                 ! ice surface boundary condition   
+      !
+      fr_i(:,:) = at_i(:,:)             ! initialisation of sea-ice fraction
+      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
+      !
+      nstart = numit  + nn_fsbc      
+      nitrun = nitend - nit000 + 1 
+      nlast  = numit  + nitrun 
+      !
+      IF( nstock == 0  )   nstock = nlast + 1
+      !
+   END SUBROUTINE ice_init
+
+
+   SUBROUTINE ice_run
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_run ***
+      !!                 
+      !! ** Purpose :   Definition some run parameter for ice model
+      !!
+      !! ** Method  :   Read the namicerun namelist and check the parameter 
+      !!              values called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namelist namicerun
+      !!-------------------------------------------------------------------
+      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, cai, cao, ln_nicep, ln_limdiahsb, ln_limdiaout
+      INTEGER  ::   ios                 ! Local integer output status for namelist read
+      !!-------------------------------------------------------------------
+      !                    
+      REWIND( numnam_ice_ref )              ! Namelist namicerun in reference namelist : Parameters for ice
+      READ  ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namicerun in configuration namelist : Parameters for ice
+      READ  ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namicerun )
+      !
+      !IF( lk_mpp .AND. ln_nicep ) THEN
+      !   ln_nicep = .FALSE.
+      !   CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' )
+      !ENDIF
+      !
+      IF(lwp) THEN                        ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice'
+         WRITE(numout,*) ' ~~~~~~'
+         WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn
+         WRITE(numout,*) '   maximum ice concentration                               = ', amax 
+         WRITE(numout,*) '   atmospheric drag over sea ice                           = ', cai
+         WRITE(numout,*) '   atmospheric drag over ocean                             = ', cao
+         WRITE(numout,*) '   Several ice points in the ice or not in ocean.output    = ', ln_nicep
+         WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb
+         WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout
+      ENDIF
+      !
+   END SUBROUTINE ice_run
+
+
+   SUBROUTINE lim_itd_ini
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE lim_itd_ini ***
+      !!
+      !! ** Purpose :   Initializes the ice thickness distribution
+      !! ** Method  :   ...
+      !!    Note    : hi_max(jpl) is here set up to a value close to 7 m for
+      !!              limistate (only) and is changed to 99 m in ice_init
+      !!------------------------------------------------------------------
+      INTEGER  ::   jl                   ! dummy loop index
+      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars
+      !!------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'lim_itd_ini : Initialization of ice thickness distribution '
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+
+      !------------------------------------------------------------------------------!
+      ! 1) Ice thickness distribution parameters initialization    
+      !------------------------------------------------------------------------------!
+      IF(lwp) THEN  
+         WRITE(numout,*) ' Number of ice categories jpl = ', jpl
+      ENDIF
+
+      !- Thickness categories boundaries 
+      !----------------------------------
+      hi_max(:) = 0._wp
+
+      zc1 =  3._wp / REAL( jpl, wp )
+      zc2 = 10._wp * zc1
+      zc3 =  3._wp
+
+      DO jl = 1, jpl
+         zx1 = REAL( jl-1, wp ) / REAL( jpl, wp )
+         hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) )
+      END DO
+
+      IF(lwp) WRITE(numout,*) ' Thickness category boundaries '
+      IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
+
+      !
+      DO jl = 1, jpl
+         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
+      END DO
+      !
+      !
+   END SUBROUTINE lim_itd_ini
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :        Empty module           NO LIM sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE ice_init        ! Empty routine
+   END SUBROUTINE ice_init
+#endif
+
+   !!======================================================================
+END MODULE iceini
Index: /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90
===================================================================
--- /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90	(revision 6448)
+++ /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90	(revision 6448)
@@ -0,0 +1,25 @@
+MODULE par_ice
+   !!======================================================================
+   !!                       ***  MODULE par_ice   ***
+   !! LIM-3 Sea Ice :   definition of parameters
+   !!======================================================================
+   !! History :  3.0  ! 2003-08  (M. Vancoppenolle)  LIM-3
+   !!----------------------------------------------------------------------
+   USE par_oce
+
+   IMPLICIT NONE
+   PUBLIC          
+
+   !                                             !!! ice thermodynamics
+   INTEGER, PUBLIC, PARAMETER ::   nlay_i   = 5   !: number of ice layers
+   INTEGER, PUBLIC, PARAMETER ::   nlay_s   = 1   !: number of snow layers
+
+   !                                             !!! ice mechanical redistribution
+   INTEGER, PUBLIC, PARAMETER ::   jpl      = 5   !: number of ice categories
+
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!======================================================================
+END MODULE par_ice
Index: /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OFF_SRC/domain.F90
===================================================================
--- /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OFF_SRC/domain.F90	(revision 6448)
+++ /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OFF_SRC/domain.F90	(revision 6448)
@@ -0,0 +1,428 @@
+MODULE domain
+   !!==============================================================================
+   !!                       ***  MODULE domain   ***
+   !! Ocean initialization : domain initialization
+   !!==============================================================================
+
+   !!----------------------------------------------------------------------
+   !!   dom_init       : initialize the space and time domain
+   !!   dom_nam        : read and contral domain namelists
+   !!   dom_ctl        : control print for the ocean domain
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! 
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distributed memory computing library
+
+   USE domstp          ! domain: set the time-step
+   USE domrea          ! domain: write the meshmask file
+   USE dommsk          ! domain : mask
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC dom_init       ! called by opa.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE dom_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_init  ***
+      !!                    
+      !! ** Purpose :   Domain initialization. Call the routines that are 
+      !!      required to create the arrays which define the space and time
+      !!      domain of the ocean model.
+      !!
+      !! ** Method  :
+      !!      - dom_stp: defined the model time step
+      !!      - dom_rea: read the meshmask file if nmsh=1
+      !!
+      !! History :
+      !!        !  90-10  (C. Levy - G. Madec)  Original code
+      !!        !  91-11  (G. Madec)
+      !!        !  92-01  (M. Imbard) insert time step initialization
+      !!        !  96-06  (G. Madec) generalized vertical coordinate 
+      !!        !  97-02  (G. Madec) creation of domwri.F
+      !!        !  01-05  (E.Durand - G. Madec) insert closed sea
+      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
+      !!----------------------------------------------------------------------
+      !! * Local declarations
+      INTEGER ::   jk                ! dummy loop argument
+      INTEGER ::   iconf = 0         ! temporary integers
+      !!----------------------------------------------------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_init : domain initialization'
+         WRITE(numout,*) '~~~~~~~~'
+      ENDIF
+
+      CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
+      CALL dom_zgr      ! Vertical mesh and bathymetry option
+      CALL dom_rea      ! Create a domain file
+
+     !
+      ! - ML - Used in dom_vvl_sf_nxt and lateral diffusion routines
+      !        but could be usefull in many other routines
+      e12t    (:,:) = e1t(:,:) * e2t(:,:)
+      e1e2t   (:,:) = e1t(:,:) * e2t(:,:)
+      e12u    (:,:) = e1u(:,:) * e2u(:,:)
+      e12v    (:,:) = e1v(:,:) * e2v(:,:)
+      e12f    (:,:) = e1f(:,:) * e2f(:,:)
+      r1_e12t (:,:) = 1._wp    / e12t(:,:)
+      r1_e12u (:,:) = 1._wp    / e12u(:,:)
+      r1_e12v (:,:) = 1._wp    / e12v(:,:)
+      r1_e12f (:,:) = 1._wp    / e12f(:,:)
+      re2u_e1u(:,:) = e2u(:,:) / e1u(:,:)
+      re1v_e2v(:,:) = e1v(:,:) / e2v(:,:)
+      !
+      hu(:,:) = 0._wp                          ! Ocean depth at U- and V-points
+      hv(:,:) = 0._wp
+      DO jk = 1, jpk
+         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+      END DO
+      !                                        ! Inverse of the local depth
+      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
+      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
+
+      CALL dom_stp      ! Time step
+      CALL dom_msk      ! Masks
+      CALL dom_ctl      ! Domain control
+
+   END SUBROUTINE dom_init
+
+   SUBROUTINE dom_nam
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dom_nam  ***
+      !!                    
+      !! ** Purpose :   read domaine namelists and print the variables.
+      !!
+      !! ** input   : - namrun namelist
+      !!              - namdom namelist
+      !!              - namcla namelist
+      !!----------------------------------------------------------------------
+      USE ioipsl
+      INTEGER  ::   ios                 ! Local integer output status for namelist read
+      NAMELIST/namrun/ nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
+         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
+         &             nn_write, ln_dimgnnn, ln_mskland  , ln_clobber   , nn_chunksz, nn_euler
+      NAMELIST/namdom/ nn_bathy , rn_bathy, rn_e3zps_min, rn_e3zps_rat, nn_msh    , rn_hmin,   &
+         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,            &
+         &             rn_rdtmax, rn_rdth     , nn_baro     , nn_closea , ln_crs, &
+         &             jphgr_msh, &
+         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
+         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
+         &             ppa2, ppkth2, ppacr2
+      NAMELIST/namcla/ nn_cla
+#if defined key_netcdf4
+      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
+#endif
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
+      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
+      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namrun )
+      !
+      IF(lwp) THEN                  ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
+         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) '   Namelist namrun'  
+         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
+         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
+         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
+         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
+         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
+         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
+         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
+         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
+         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
+         WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
+         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
+         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
+         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
+         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
+         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
+      ENDIF
+      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
+      cexper = cn_exp
+      nrstdt = nn_rstctl
+      nit000 = nn_it000
+      nitend = nn_itend
+      ndate0 = nn_date0
+      nleapy = nn_leapy
+      ninist = nn_istate
+      nstock = nn_stock
+      nwrite = nn_write
+
+
+      !                             ! control of output frequency
+      IF ( nstock == 0 .OR. nstock > nitend ) THEN
+         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
+         CALL ctl_warn( ctmp1 )
+         nstock = nitend
+      ENDIF
+      IF ( nwrite == 0 ) THEN
+         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
+         CALL ctl_warn( ctmp1 )
+         nwrite = nitend
+      ENDIF
+
+      ! parameters correspondting to nit000 - 1 (as we start the step loop with a call to day)
+      ndastp = ndate0 - 1        ! ndate0 read in the namelist in dom_nam, we assume that we start run at 00:00
+      adatrj = ( REAL( nit000-1, wp ) * rdttra(1) ) / rday
+
+#if defined key_agrif
+      IF( Agrif_Root() ) THEN
+#endif
+      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
+      CASE (  1 ) 
+         CALL ioconf_calendar('gregorian')
+         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
+      CASE (  0 )
+         CALL ioconf_calendar('noleap')
+         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
+      CASE ( 30 )
+         CALL ioconf_calendar('360d')
+         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
+      END SELECT
+#if defined key_agrif
+      ENDIF
+#endif
+
+      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
+      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
+903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
+      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
+904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namdom )
+
+      IF(lwp) THEN
+         WRITE(numout,*) 
+         WRITE(numout,*) '   Namelist namdom : space & time domain'
+         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
+         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
+         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
+         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
+         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
+         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
+         WRITE(numout,*) '           = 0   no file created                 '
+         WRITE(numout,*) '           = 1   mesh_mask                       '
+         WRITE(numout,*) '           = 2   mesh and mask                   '
+         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask      '
+         WRITE(numout,*) '      ocean time step                      rn_rdt    = ', rn_rdt
+         WRITE(numout,*) '      asselin time filter parameter        rn_atfp   = ', rn_atfp
+         WRITE(numout,*) '      time-splitting: nb of sub time-step  nn_baro   = ', nn_baro
+         WRITE(numout,*) '      acceleration of converge             nn_acc    = ', nn_acc
+         WRITE(numout,*) '        nn_acc=1: surface tracer rdt       rn_rdtmin = ', rn_rdtmin
+         WRITE(numout,*) '                  bottom  tracer rdt       rdtmax    = ', rn_rdtmax
+         WRITE(numout,*) '                  depth of transition      rn_rdth   = ', rn_rdth
+         WRITE(numout,*) '      suppression of closed seas (=0)      nn_closea = ', nn_closea
+         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
+         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
+         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
+         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
+         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
+         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
+         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
+         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
+         WRITE(numout,*) '                                        ppa0            = ', ppa0
+         WRITE(numout,*) '                                        ppa1            = ', ppa1
+         WRITE(numout,*) '                                        ppkth           = ', ppkth
+         WRITE(numout,*) '                                        ppacr           = ', ppacr
+         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
+         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
+         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
+         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
+         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
+         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
+      ENDIF
+
+      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
+      e3zps_min = rn_e3zps_min
+      e3zps_rat = rn_e3zps_rat
+      nmsh      = nn_msh
+      nacc      = nn_acc
+      atfp      = rn_atfp
+      rdt       = rn_rdt
+      rdtmin    = rn_rdtmin
+      rdtmax    = rn_rdtmin
+      rdth      = rn_rdth
+
+      REWIND( numnam_ref )              ! Namelist namcla in reference namelist : Cross land advection
+      READ  ( numnam_ref, namcla, IOSTAT = ios, ERR = 905)
+905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namcla in configuration namelist : Cross land advection
+      READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )
+906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )
+      IF(lwm) WRITE( numond, namcla )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '   Namelist namcla'
+         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
+      ENDIF
+
+#if defined key_netcdf4
+      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
+      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
+      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
+907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
+      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
+908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
+      IF(lwm) WRITE( numond, namnc4 )
+      IF(lwp) THEN                        ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
+         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
+         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
+         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
+         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
+      ENDIF
+
+      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
+      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
+      snc4set%ni   = nn_nchunks_i
+      snc4set%nj   = nn_nchunks_j
+      snc4set%nk   = nn_nchunks_k
+      snc4set%luse = ln_nc4zip
+#else
+      snc4set%luse = .FALSE.        ! No NetCDF 4 case
+#endif
+      !
+   END SUBROUTINE dom_nam
+
+   SUBROUTINE dom_zgr
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_zgr  ***
+      !!                   
+      !! ** Purpose :  set the depth of model levels and the resulting 
+      !!      vertical scale factors.
+      !!
+      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d)
+      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
+      !!              - vertical coordinate (gdep., e3.) depending on the 
+      !!                coordinate chosen :
+      !!                   ln_zco=T   z-coordinate  
+      !!                   ln_zps=T   z-coordinate with partial steps
+      !!                   ln_zco=T   s-coordinate 
+      !!
+      !! ** Action  :   define gdep., e3., mbathy and bathy
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio = 0   ! temporary integer
+      INTEGER ::   ios
+      !!
+      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam_ref )              ! Namelist namzgr in reference namelist : Vertical coordinate
+      READ  ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namzgr in configuration namelist : Vertical coordinate
+      READ  ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namzgr )
+
+      IF(lwp) THEN                     ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_zgr : vertical coordinate'
+         WRITE(numout,*) '~~~~~~~'
+         WRITE(numout,*) '          Namelist namzgr : set vertical coordinate'
+         WRITE(numout,*) '             z-coordinate - full steps      ln_zco    = ', ln_zco
+         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps    = ', ln_zps
+         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
+         WRITE(numout,*) '             ice shelf cavity               ln_isfcav = ', ln_isfcav
+      ENDIF
+
+      ioptio = 0                       ! Check Vertical coordinate options
+      IF( ln_zco ) ioptio = ioptio + 1
+      IF( ln_zps ) ioptio = ioptio + 1
+      IF( ln_sco ) ioptio = ioptio + 1
+      IF( ln_isfcav ) ioptio = 33
+      IF ( ioptio /= 1  )   CALL ctl_stop( ' none or several vertical coordinate options used' )
+      IF ( ioptio == 33 )   CALL ctl_stop( ' isf cavity with off line module not yet done    ' )
+
+   END SUBROUTINE dom_zgr
+
+   SUBROUTINE dom_ctl
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dom_ctl  ***
+      !!
+      !! ** Purpose :   Domain control.
+      !!
+      !! ** Method  :   compute and print extrema of masked scale factors
+      !!
+      !! History :
+      !!   8.5  !  02-08  (G. Madec)    Original code
+      !!----------------------------------------------------------------------
+      !! * Local declarations
+      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
+      INTEGER, DIMENSION(2) ::   iloc      ! 
+      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
+      !!----------------------------------------------------------------------
+
+      ! Extrema of the scale factors
+
+      IF(lwp)WRITE(numout,*)
+      IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
+      IF(lwp)WRITE(numout,*) '~~~~~~~'
+
+      IF (lk_mpp) THEN
+         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
+         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
+         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
+         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
+      ELSE
+         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+
+         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iimi1 = iloc(1) + nimpp - 1
+         ijmi1 = iloc(2) + njmpp - 1
+         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iimi2 = iloc(1) + nimpp - 1
+         ijmi2 = iloc(2) + njmpp - 1
+         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iima1 = iloc(1) + nimpp - 1
+         ijma1 = iloc(2) + njmpp - 1
+         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iima2 = iloc(1) + nimpp - 1
+         ijma2 = iloc(2) + njmpp - 1
+      ENDIF
+
+      IF(lwp) THEN
+         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
+         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
+         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
+         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
+      ENDIF
+
+   END SUBROUTINE dom_ctl
+
+   !!======================================================================
+END MODULE domain
Index: /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OFF_SRC/dommsk.F90
===================================================================
--- /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OFF_SRC/dommsk.F90	(revision 6448)
+++ /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OFF_SRC/dommsk.F90	(revision 6448)
@@ -0,0 +1,104 @@
+MODULE dommsk
+   !!======================================================================
+   !!                       ***  MODULE dommsk   ***
+   !! Ocean initialization : domain land/sea masks, off-line case 
+   !!======================================================================
+   !! History :  3.3  ! 2010-10  (C. Ethe)  adapted from OPA_SRC/DOM/dommsk
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dom_msk        : compute land/ocean mask
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE lib_mpp         ! MPP library
+   USE in_out_manager  ! I/O manager
+   USE wrk_nemo  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   dom_msk    ! routine called by inidom.F90
+
+   REAL(wp)        :: rn_shlat   = 2.   ! type of lateral boundary condition on velocity
+   LOGICAL, PUBLIC :: ln_vorlat  = .false.   !  consistency of vorticity boundary condition 
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dom_msk
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE dom_msk  ***
+      !!
+      !! ** Purpose :   Off-line case: defines the interior domain T-mask.
+      !!
+      !! ** Method  :   The interior ocean/land mask is computed from tmask
+      !!              setting to zero the duplicated row and lines due to
+      !!              MPP exchange halos, est-west cyclic and north fold
+      !!              boundary conditions.
+      !!
+      !! ** Action :   tmask_i  : interiorland/ocean mask at t-point
+      !!               tpol     : ???
+      !!----------------------------------------------------------------------
+      !
+      INTEGER  ::   ji, jk                   ! dummy loop indices
+      INTEGER  ::   iif, iil, ijf, ijl       ! local integers
+      INTEGER, POINTER, DIMENSION(:,:) ::  imsk 
+      !
+      !!---------------------------------------------------------------------
+      
+      CALL wrk_alloc( jpi, jpj, imsk )
+      !
+      ! Interior domain mask (used for global sum)
+      ! --------------------
+      ssmask(:,:)  = tmask(:,:,1)
+      tmask_i(:,:) = tmask(:,:,1)
+      iif = jpreci                        ! thickness of exchange halos in i-axis
+      iil = nlci - jpreci + 1
+      ijf = jprecj                        ! thickness of exchange halos in j-axis
+      ijl = nlcj - jprecj + 1
+      !
+      tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns
+      tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns)
+      tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows
+      tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows)
+      !
+      !                                   ! north fold mask
+      tpol(1:jpiglo) = 1._wp
+      !                                
+      IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot
+      IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot
+      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row
+         IF( mjg(ijl-1) == jpjglo-1 ) THEN
+            DO ji = iif+1, iil-1
+               tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji))
+            END DO
+         ENDIF
+      ENDIF 
+      !
+      IF( nprint == 1 .AND. lwp ) THEN    ! Control print
+         imsk(:,:) = INT( tmask_i(:,:) )
+         WRITE(numout,*) ' tmask_i : '
+         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
+         WRITE (numout,*)
+         WRITE (numout,*) ' dommsk: tmask for each level'
+         WRITE (numout,*) ' ----------------------------'
+         DO jk = 1, jpk
+            imsk(:,:) = INT( tmask(:,:,jk) )
+            WRITE(numout,*)
+            WRITE(numout,*) ' level = ',jk
+            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
+         END DO
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, imsk )
+      !
+   END SUBROUTINE dom_msk
+   !!======================================================================
+END MODULE dommsk
Index: /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OFF_SRC/istate.F90
===================================================================
--- /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OFF_SRC/istate.F90	(revision 6448)
+++ /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OFF_SRC/istate.F90	(revision 6448)
@@ -0,0 +1,51 @@
+MODULE istate
+   !!======================================================================
+   !!                     ***  MODULE  istate  ***
+   !! Ocean state   :  initial state setting, off-line case
+   !!=====================================================================
+   !! History :  3.3  ! 2010-10  (C. Ethe)  original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   istate_init   : initial state set to zero
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and active tracers 
+   USE dom_oce         ! ocean space and time domain 
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   istate_init   ! routine called by step.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!---------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE istate_init
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE istate_init  ***
+      !! 
+      !! ** Purpose :   Initialization to zero of the dynamics and tracers.
+      !!----------------------------------------------------------------------
+      !
+      !     now fields         !     after fields      !
+      un   (:,:,:)   = 0._wp   ;   ua(:,:,:) = 0._wp   !
+      vn   (:,:,:)   = 0._wp   ;   va(:,:,:) = 0._wp   !
+      wn   (:,:,:)   = 0._wp   !                       !
+      hdivn(:,:,:)   = 0._wp   !                       !
+      tsn  (:,:,:,:) = 0._wp   !                       !
+      !
+      rhd  (:,:,:) = 0.e0
+      rhop (:,:,:) = 0.e0
+      rn2  (:,:,:) = 0.e0 
+      !
+   END SUBROUTINE istate_init
+
+   !!=====================================================================
+END MODULE istate
Index: /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OFF_SRC/stpctl.F90
===================================================================
--- /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OFF_SRC/stpctl.F90	(revision 6448)
+++ /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OFF_SRC/stpctl.F90	(revision 6448)
@@ -0,0 +1,62 @@
+MODULE stpctl
+   !!======================================================================
+   !!                       ***  MODULE  stpctl  ***
+   !! Ocean run control :  Off-line case, only save the time step in numstp
+   !!======================================================================
+   !! History :  OPA  ! 1991-03  (G. Madec) Original code
+   !!            6.0  ! 1992-06  (M. Imbard)
+   !!            8.0  ! 1997-06  (A.M. Treguier)
+   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
+   !!            2.0  ! 2009-07  (G. Madec)  Add statistic for time-spliting
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   stp_ctl      : Control the run
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables 
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp         ! distributed memory computing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   stp_ctl    ! routine called by opa.F90
+   
+   !!----------------------------------------------------------------------
+   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE stp_ctl( kt, kindic )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE stp_ctl  ***
+      !!                     
+      !! ** Purpose :   Control the run
+      !!
+      !! ** Method  : - Save the time step in numstp
+      !!
+      !! ** Actions :   'time.step' file containing the last ocean time-step
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in   ) ::   kt      ! ocean time-step index
+      INTEGER, INTENT(inout) ::   kindic  ! indicator of solver convergence
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'stp_ctl : time-stepping control'
+         WRITE(numout,*) '~~~~~~~'
+         ! open time.step file
+         CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+      ENDIF
+      !
+      IF(lwp) WRITE ( numstp, '(1x, i8)' )   kt      !* save the current time step in numstp
+      IF(lwp) REWIND( numstp )                       !  --------------------------
+      !
+   END SUBROUTINE stp_ctl
+
+   !!======================================================================
+END MODULE stpctl
Index: /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 6447)
+++ /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 6448)
@@ -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
@@ -228,4 +229,12 @@
          CALL iom_put( "sbv", z2d )                ! bottom j-current
       ENDIF
+      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
+      IF( lk_zdfddm ) THEN
+         CALL iom_put( "avs" , fsavs(:,:,:)                          )    ! S vert. eddy diff. coef.
 #if defined key_dynspg_ts
       CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic j-current
Index: /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6447)
+++ /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6448)
@@ -104,5 +104,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'
@@ -147,4 +147,5 @@
       ! no need to output in offline mode
       IF( .NOT.lk_offline ) THEN   
+      IF( kt >= nit000 ) THEN               ! workaround for calls before SOMETHING reads the XIOS namelist
          IF ( iom_use("mldr10_1") ) THEN
             IF( ln_isfcav ) THEN
@@ -161,4 +162,5 @@
             END IF
          END IF
+      ENDIF
       ENDIF
       
Index: /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 6447)
+++ /branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 6448)
@@ -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), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
+   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
@@ -118,5 +122,7 @@
          &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          &
 #endif
-         &      htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) , STAT= zdf_tke_alloc      )
+         &      htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,                          &
+         &      efr  (jpi,jpj)     , e_niw(jpi,jpj,jpk) ,                          &
+         &      STAT= zdf_tke_alloc      )
          !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
@@ -428,7 +434,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
@@ -465,6 +486,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
@@ -730,7 +777,12 @@
          &                 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   
+      !!----------------------------------------------------------------------
+      !
+      ! NB. Default values of namelist parameters should be set in the reference namelist
+      !     but setting this one here because we aren't set up to use the reference namelist
+      !     properly yet. 
+      rn_c = 0.8_wp
+
       REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy
       READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
@@ -776,4 +828,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
@@ -786,5 +839,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' )
 
@@ -794,11 +847,9 @@
       ENDIF
       
-      IF( nn_etau == 2  ) THEN
-          ierr = zdf_mxl_alloc()
-          nmln(:,:) = nlb10           ! Initialization of nmln
-      ENDIF
+      IF( nn_etau == 2  .OR. ( nn_etau /= 0 .AND. nn_htau == 2 ) )   CALL zdf_mxl( nit000 - 1 )      ! Initialization of nmln and hmlp
 
       !                               !* 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)
@@ -806,5 +857,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
