Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 5382)
@@ -24,5 +24,7 @@
    !                                                   !!* namicerun read in iceini  *
    CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in     !: suffix of ice restart name (input)
+   CHARACTER(len=256)    , PUBLIC ::   cn_icerst_indir  !: ice restart in directory
    CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out    !: suffix of ice restart name (output)
+   CHARACTER(len=256)    , PUBLIC ::   cn_icerst_outdir !: ice restart out directory
    LOGICAL               , PUBLIC ::   ln_limdyn        !: flag for ice dynamics (T) or not (F)
    LOGICAL               , PUBLIC ::   ln_limdmp        !: Ice damping
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90	(revision 5382)
@@ -40,5 +40,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/LIM2 4.0 , UCL - NEMO Consortium (2011)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -123,5 +123,6 @@
       !! ** input   :   Namelist namicerun
       !!-------------------------------------------------------------------
-      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, ln_limdmp, acrit, hsndif, hicdif
+      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir, &
+                          ln_limdyn, ln_limdmp, acrit, hsndif, hicdif
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!-------------------------------------------------------------------
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_2/limrst_2.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_2/limrst_2.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_2/limrst_2.F90	(revision 5382)
@@ -50,4 +50,5 @@
       CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character
       CHARACTER(LEN=50)   ::   clname   ! ice output restart file name
+      CHARACTER(len=150)  ::   clpath   ! full path to ice output restart file
       !!----------------------------------------------------------------------
       !
@@ -58,24 +59,30 @@
       ! except if we write ice restart files every ice time step or if an ice restart file was writen at nitend - 2*nn_fsbc + 1
       IF( kt == nitrst - 2*nn_fsbc + 1 .OR. nstock == nn_fsbc .OR. ( kt == nitend - nn_fsbc + 1 .AND. .NOT. lrst_ice ) ) THEN
-         ! beware of the format used to write kt (default is i8.8, that should be large enough...)
-         IF( nitrst > 99999999 ) THEN   ;   WRITE(clkt, *       ) nitrst
-         ELSE                           ;   WRITE(clkt, '(i8.8)') nitrst
+         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 > 99999999 ) THEN   ;   WRITE(clkt, *       ) nitrst
+            ELSE                           ;   WRITE(clkt, '(i8.8)') nitrst
+            ENDIF
+            ! create the file
+            clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_icerst_out)
+            clpath = TRIM(cn_icerst_outdir) 
+            IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath)//'/' 
+            IF(lwp) THEN
+               WRITE(numout,*)
+               SELECT CASE ( jprstlib )
+               CASE ( jprstdimg )
+                  WRITE(numout,*) '             open ice restart binary file: ',TRIM(clpath)//clname
+               CASE DEFAULT
+                  WRITE(numout,*) '             open ice restart NetCDF file: ',TRIM(clpath)//clname
+               END SELECT
+               IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN   
+                  WRITE(numout,*)         '             kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp
+               ELSE   ;   WRITE(numout,*) '             kt = '                         , kt,' date= ', ndastp
+               ENDIF
+            ENDIF
+
+            CALL iom_open( TRIM(clpath)//TRIM(clname), numriw, ldwrt = .TRUE., kiolib = jprstlib )
+            lrst_ice = .TRUE.
          ENDIF
-         ! create the file
-         clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_icerst_out)
-         IF(lwp) THEN
-            WRITE(numout,*)
-            SELECT CASE ( jprstlib )
-            CASE ( jprstdimg )   ;   WRITE(numout,*) '             open ice restart binary file: '//clname
-            CASE DEFAULT         ;   WRITE(numout,*) '             open ice restart NetCDF file: '//clname
-            END SELECT
-            IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN   
-               WRITE(numout,*)         '             kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp
-            ELSE   ;   WRITE(numout,*) '             kt = '                         , kt,' date= ', ndastp
-            ENDIF
-         ENDIF
-
-         CALL iom_open( clname, numriw, ldwrt = .TRUE., kiolib = jprstlib )
-         lrst_ice = .TRUE.
       ENDIF
       !
@@ -188,9 +195,9 @@
         ! eventually read netcdf file (monobloc)  for restarting on different number of processors
         ! if {cn_icerst_in}.nc exists, then set jlibalt to jpnf90
-        INQUIRE( FILE = TRIM(cn_icerst_in)//'.nc', EXIST = llok )
+        INQUIRE( FILE = TRIM(cn_icerst_indir)//'/'//TRIM(cn_icerst_in)//'.nc', EXIST = llok )
         IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF
       ENDIF
 
-      CALL iom_open ( cn_icerst_in, numrir, kiolib = jlibalt )
+      CALL iom_open ( TRIM(cn_icerst_indir)//'/'//TRIM(cn_icerst_in), numrir, kiolib = jlibalt )
 
       CALL iom_get( numrir, 'kt_ice' , ziter )
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5382)
@@ -373,5 +373,7 @@
    INTEGER          , PUBLIC ::   nlay_s          !: number of snow layers 
    CHARACTER(len=32), PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input)
+   CHARACTER(len=256), PUBLIC ::   cn_icerst_indir !: ice restart input directory
    CHARACTER(len=32), PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output)
+   CHARACTER(len=256), PUBLIC ::   cn_icerst_outdir!: ice restart output directory
    LOGICAL          , PUBLIC ::   ln_limdyn       !: flag for ice dynamics (T) or not (F)
    LOGICAL          , PUBLIC ::   ln_icectl       !: flag for sea-ice points output (T) or not (F)
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 5382)
@@ -55,4 +55,5 @@
       CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
       CHARACTER(LEN=50)   ::   clname   ! ice output restart file name
+      CHARACTER(len=256)  ::   clpath   ! full path to ice output restart file 
       !!----------------------------------------------------------------------
       !
@@ -64,24 +65,30 @@
       IF( kt == nitrst - 2*nn_fsbc + 1 .OR. nstock == nn_fsbc    &
          &                             .OR. ( kt == nitend - nn_fsbc + 1 .AND. .NOT. lrst_ice ) ) THEN
-         ! beware of the format used to write kt (default is i8.8, that should be large enough...)
-         IF( nitrst > 99999999 ) THEN   ;   WRITE(clkt, *       ) nitrst
-         ELSE                           ;   WRITE(clkt, '(i8.8)') nitrst
+         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 > 99999999 ) THEN   ;   WRITE(clkt, *       ) nitrst
+            ELSE                           ;   WRITE(clkt, '(i8.8)') nitrst
+            ENDIF
+            ! create the file
+            clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_icerst_out)
+            clpath = TRIM(cn_icerst_outdir) 
+            IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath)//'/'
+            IF(lwp) THEN
+               WRITE(numout,*)
+               SELECT CASE ( jprstlib )
+               CASE ( jprstdimg )
+                  WRITE(numout,*) '             open ice restart binary file: ',TRIM(clpath)//clname
+               CASE DEFAULT
+                  WRITE(numout,*) '             open ice restart NetCDF file: ',TRIM(clpath)//clname
+               END SELECT
+               IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN   
+                  WRITE(numout,*)         '             kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp
+               ELSE   ;   WRITE(numout,*) '             kt = '                         , kt,' date= ', ndastp
+               ENDIF
+            ENDIF
+            !
+            CALL iom_open( TRIM(clpath)//TRIM(clname), numriw, ldwrt = .TRUE., kiolib = jprstlib )
+            lrst_ice = .TRUE.
          ENDIF
-         ! create the file
-         clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_icerst_out)
-         IF(lwp) THEN
-            WRITE(numout,*)
-            SELECT CASE ( jprstlib )
-            CASE ( jprstdimg )   ;   WRITE(numout,*) '             open ice restart binary file: '//clname
-            CASE DEFAULT         ;   WRITE(numout,*) '             open ice restart NetCDF file: '//clname
-            END SELECT
-            IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN   
-               WRITE(numout,*)         '             kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp
-            ELSE   ;   WRITE(numout,*) '             kt = '                         , kt,' date= ', ndastp
-            ENDIF
-         ENDIF
-         !
-         CALL iom_open( clname, numriw, ldwrt = .TRUE., kiolib = jprstlib )
-         lrst_ice = .TRUE.
       ENDIF
       !
@@ -143,5 +150,5 @@
          CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
       END DO
-      
+
       DO jl = 1, jpl 
          WRITE(zchar,'(I1)') jl
@@ -327,9 +334,9 @@
         ! eventually read netcdf file (monobloc)  for restarting on different number of processors
         ! if {cn_icerst_in}.nc exists, then set jlibalt to jpnf90
-        INQUIRE( FILE = TRIM(cn_icerst_in)//'.nc', EXIST = llok )
+        INQUIRE( FILE = TRIM(cn_icerst_indir)//'/'//TRIM(cn_icerst_in)//'.nc', EXIST = llok )
         IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF
       ENDIF
 
-      CALL iom_open ( cn_icerst_in, numrir, kiolib = jprstlib )
+      CALL iom_open ( TRIM(cn_icerst_indir)//'/'//cn_icerst_in, numrir, kiolib = jprstlib )
 
       CALL iom_get( numrir, 'nn_fsbc', zfice )
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OFF_SRC/domain.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OFF_SRC/domain.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OFF_SRC/domain.F90	(revision 5382)
@@ -116,7 +116,8 @@
       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,   &
+      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_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
-         &             nn_write, ln_dimgnnn, ln_mskland  , ln_clobber   , nn_chunksz, nn_euler
+         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , 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 ,            &
@@ -159,4 +160,5 @@
          WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
          WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
+         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
          WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
          WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
@@ -171,4 +173,5 @@
       ninist = nn_istate
       nstock = nn_stock
+      nstocklist = nn_stocklist
       nwrite = nn_write
 
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90	(revision 5382)
@@ -1883,5 +1883,5 @@
       CALL crs_lbc_lnk( p_surf_crs_msk, cd_type, 1.0, pval=1.0 )
 
-      CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk )
+      CALL wrk_dealloc( jpi, jpj, jpk, zsurfmsk, zsurf )
 
    END SUBROUTINE crs_dom_sfc
@@ -2275,6 +2275,4 @@
       ENDDO
      
-      CALL wrk_alloc( jpi_crs, jpj_crs, zmbk )
-
       zmbk(:,:) = 0.0
       zmbk(:,:) = REAL( mbathy_crs(:,:), wp ) ;   CALL crs_lbc_lnk(zmbk,'T',1.0)   ;   mbathy_crs(:,:) = INT( zmbk(:,:) )
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90	(revision 5382)
@@ -245,5 +245,4 @@
       CALL iom_put( "mldr10_3", zrho10_3     )   ! MLD delta rho(10m) = 0.03
       CALL iom_put( "pycndep" , zpycn        )   ! MLD delta rho equi. delta T(10m) = 0.2
-      CALL iom_put( "BLT"     , ztm2 - zpycn )   ! Barrier Layer Thickness
       CALL iom_put( "tinv"    , ztinv        )   ! max. temp. inv. (t10 ref) 
       CALL iom_put( "depti"   , zdepinv      )   ! depth of max. temp. inv. (t10 ref) 
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 5382)
@@ -135,7 +135,8 @@
       !!----------------------------------------------------------------------
       USE ioipsl
-      NAMELIST/namrun/ nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
+      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_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
-         &             nn_write, ln_dimgnnn, ln_mskland  , ln_clobber   , nn_chunksz, nn_euler
+         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , 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 ,                  &
@@ -169,5 +170,7 @@
          WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
          WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in
+         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir
          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,*) '      start with forward time step    nn_euler   = ', nn_euler
@@ -178,8 +181,13 @@
          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
+         IF( ln_rst_list ) THEN
+            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
+         ELSE
+            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
+         ENDIF
          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,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
          WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
          WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
@@ -195,7 +203,8 @@
       ninist = nn_istate
       nstock = nn_stock
+      nstocklist = nn_stocklist
       nwrite = nn_write
       neuler = nn_euler
-      IF ( neuler == 1 .AND. .NOT.ln_rstart ) THEN
+      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
          WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
          CALL ctl_warn( ctmp1 )
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5382)
@@ -472,25 +472,4 @@
          risfdep(:,:)=0.e0
          misfdep(:,:)=1
-         !
-         ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code
-         IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN 
-            risfdep(:,:)=200.e0 
-            misfdep(:,:)=1 
-            ij0 = 1 ; ij1 = 40 
-            DO jj = mj0(ij0), mj1(ij1) 
-               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
-            END DO 
-            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
-         ! 
-         ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN
-         ! 
-            risfdep(:,:)=0.e0
-            misfdep(:,:)=1
-            ij0 = 1 ; ij1 = 40
-            DO jj = mj0(ij0), mj1(ij1)
-               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp
-            END DO
-            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp
-         END IF
          !
          DEALLOCATE( idta, zdta )
@@ -969,5 +948,5 @@
       !!
       INTEGER  ::   ji, jj, jk       ! dummy loop indices
-      INTEGER  ::   ik, it           ! temporary integers
+      INTEGER  ::   ik, it, ikb, ikt ! temporary integers
       LOGICAL  ::   ll_print         ! Allow  control print for debugging
       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
@@ -1152,17 +1131,18 @@
       IF ( ln_isfcav ) THEN
       ! (ISF) define e3uw (adapted for 2 cells in the water column)
-      ! Need to test if the modification of only mikt and mbkt levels is enough
-         DO jk = 2,jpk                          
-            DO jj = 1, jpjm1 
-               DO ji = 1, fs_jpim1   ! vector opt. 
-                  e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) &
-                    &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) )
-                  e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) &
-                    &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) )
-               END DO 
-            END DO 
+         DO jj = 2, jpjm1 
+            DO ji = 2, fs_jpim1   ! vector opt. 
+               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
+               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
+               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) &
+                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) )
+               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
+               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
+               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) &
+                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) )
+            END DO
          END DO
       END IF
-      
+
       CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
       CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
@@ -1538,5 +1518,5 @@
  
  ! remove single point "bay" on isf coast line in the ice shelf draft'
-         DO jk = 1, jpk
+         DO jk = 2, jpk
             WHERE (misfdep==0) misfdep=jpk
             zmask=0
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 5382)
@@ -110,9 +110,4 @@
          ELSEIF( cp_cfg == 'gyre' ) THEN         
             CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields
-        ELSEIF( cp_cfg == 'isomip' .OR. cp_cfg == 'isomip2') THEN
-            IF(lwp) WRITE(numout,*) 'Initialization of T+S for ISOMIP domain' 
-            tsn(:,:,:,jp_tem)=-1.9*tmask(:,:,:)          ! ISOMIP configuration : start from constant T+S fields 
-            tsn(:,:,:,jp_sal)=34.4*tmask(:,:,:)
-            tsb(:,:,:,:)=tsn(:,:,:,:)  
          ELSE                                    ! Initial T-S, U-V fields read in files
             IF ( ln_tsd_init ) THEN              ! read 3D T and S data at nit000
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 5382)
@@ -5,5 +5,6 @@
    !!==============================================================================
    !! History :  1.0  !  2006-11  (G. Madec)  Original code
-   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!            3.3  !  2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
+   !!            3.6  !  2015-05  (N. Ducousso, G. Madec)  add Hollingsworth scheme as an option 
    !!----------------------------------------------------------------------
 
@@ -17,4 +18,5 @@
    USE dynkeg          ! kinetic energy gradient          (dyn_keg      routine)
    USE dynzad          ! vertical advection               (dyn_zad      routine)
+   !
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
@@ -25,7 +27,9 @@
 
    PUBLIC dyn_adv       ! routine called by step module
-   PUBLIC dyn_adv_init  ! routine called by opa module
+   PUBLIC dyn_adv_init  ! routine called by opa  module
  
+   !                                    !* namdyn_adv namelist *
    LOGICAL, PUBLIC ::   ln_dynadv_vec   !: vector form flag
+   INTEGER, PUBLIC ::   nn_dynkeg       !: scheme of kinetic energy gradient: =0 C2 ; =1 Hollingsworth
    LOGICAL, PUBLIC ::   ln_dynadv_cen2  !: flux form - 2nd order centered scheme flag
    LOGICAL, PUBLIC ::   ln_dynadv_ubs   !: flux form - 3rd order UBS scheme flag
@@ -38,5 +42,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.6 , NEMO Consortium (2015)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -63,16 +67,16 @@
       SELECT CASE ( nadv )                  ! compute advection trend and add it to general trend
       CASE ( 0 )     
-                      CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy
-                      CALL dyn_zad     ( kt )    ! vector form : vertical advection
+                      CALL dyn_keg     ( kt, nn_dynkeg )    ! vector form : horizontal gradient of kinetic energy
+                      CALL dyn_zad     ( kt )               ! vector form : vertical advection
       CASE ( 1 )     
-                      CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy
-                      CALL dyn_zad_zts ( kt )    ! vector form : vertical advection with sub-timestepping
+                      CALL dyn_keg     ( kt, nn_dynkeg )    ! vector form : horizontal gradient of kinetic energy
+                      CALL dyn_zad_zts ( kt )               ! vector form : vertical advection with sub-timestepping
       CASE ( 2 ) 
-                      CALL dyn_adv_cen2( kt )    ! 2nd order centered scheme
+                      CALL dyn_adv_cen2( kt )               ! 2nd order centered scheme
       CASE ( 3 )   
-                      CALL dyn_adv_ubs ( kt )    ! 3rd order UBS      scheme
+                      CALL dyn_adv_ubs ( kt )               ! 3rd order UBS      scheme
       !
-      CASE (-1 )                                 ! esopa: test all possibility with control print
-                      CALL dyn_keg     ( kt )
+      CASE (-1 )                                            ! esopa: test all possibility with control print
+                      CALL dyn_keg     ( kt, nn_dynkeg )
                       CALL dyn_zad     ( kt )
                       CALL dyn_adv_cen2( kt )
@@ -92,10 +96,9 @@
       !!              momentum advection formulation & scheme and set nadv
       !!----------------------------------------------------------------------
-      INTEGER ::   ioptio
-      INTEGER  ::   ios                 ! Local integer output status for namelist read
-      !!
-      NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts
+      INTEGER ::   ioptio, ios   ! Local integer
+      !
+      NAMELIST/namdyn_adv/ ln_dynadv_vec, nn_dynkeg, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts
       !!----------------------------------------------------------------------
-
+      !
       REWIND( numnam_ref )              ! Namelist namdyn_adv in reference namelist : Momentum advection scheme
       READ  ( numnam_ref, namdyn_adv, IOSTAT = ios, ERR = 901)
@@ -112,8 +115,9 @@
          WRITE(numout,*) '~~~~~~~~~~~'
          WRITE(numout,*) '       Namelist namdyn_adv : chose a advection formulation & scheme for momentum'
-         WRITE(numout,*) '          Vector/flux form (T/F)             ln_dynadv_vec  = ', ln_dynadv_vec
-         WRITE(numout,*) '          2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2
-         WRITE(numout,*) '          3rd order UBS advection scheme     ln_dynadv_ubs  = ', ln_dynadv_ubs
-         WRITE(numout,*) '      Sub timestepping of vertical advection ln_dynzad_zts  = ', ln_dynzad_zts
+         WRITE(numout,*) '          Vector/flux form (T/F)                           ln_dynadv_vec  = ', ln_dynadv_vec
+         WRITE(numout,*) '          = 0 standard scheme  ; =1 Hollingsworth scheme   nn_dynkeg      = ', nn_dynkeg
+         WRITE(numout,*) '          2nd order centred advection scheme               ln_dynadv_cen2 = ', ln_dynadv_cen2
+         WRITE(numout,*) '          3rd order UBS advection scheme                   ln_dynadv_ubs  = ', ln_dynadv_ubs
+         WRITE(numout,*) '          Sub timestepping of vertical advection           ln_dynzad_zts  = ', ln_dynzad_zts
       ENDIF
 
@@ -126,7 +130,7 @@
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' )
       IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec )   &
-          CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' )
-      IF( ln_dynzad_zts .AND. ln_isfcav )   &
-          CALL ctl_stop( 'Sub timestepping of vertical advection does not work with ln_isfcav = .TRUE.' )
+         CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' )
+      IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_HW )   &  
+         CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' )
 
       !                               ! Set nadv
@@ -139,6 +143,10 @@
       IF(lwp) THEN                    ! Print the choice
          WRITE(numout,*)
-         IF( nadv ==  0 )   WRITE(numout,*) '         vector form : keg + zad + vor is used'
+         IF( nadv ==  0 )   WRITE(numout,*) '         vector form : keg + zad + vor is used' 
          IF( nadv ==  1 )   WRITE(numout,*) '         vector form : keg + zad_zts + vor is used'
+         IF( nadv ==  0 .OR. nadv ==  1 ) THEN
+            IF( nn_dynkeg == nkeg_C2  )   WRITE(numout,*) 'with Centered standard keg scheme'
+            IF( nn_dynkeg == nkeg_HW  )   WRITE(numout,*) 'with Hollingsworth keg scheme'
+         ENDIF
          IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used'
          IF( nadv ==  3 )   WRITE(numout,*) '         flux form   : UBS       scheme is used'
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 5382)
@@ -4,7 +4,8 @@
    !! Ocean dynamics:  kinetic energy gradient trend
    !!======================================================================
-   !! History :  1.0  !  87-09  (P. Andrich, m.-a. Foujols)  Original code
-   !!            7.0  !  97-05  (G. Madec)  Split dynber into dynkeg and dynhpg
-   !!            9.0  !  02-07  (G. Madec)  F90: Free form and module
+   !! History :  1.0  !  1987-09  (P. Andrich, M.-A. Foujols)  Original code
+   !!            7.0  !  1997-05  (G. Madec)  Split dynber into dynkeg and dynhpg
+   !!  NEMO      1.0  !  2002-07  (G. Madec)  F90: Free form and module
+   !!            3.6  !  2015-05  (N. Ducousso, G. Madec)  add Hollingsworth scheme as an option 
    !!----------------------------------------------------------------------
    
@@ -18,4 +19,5 @@
    !
    USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp         ! MPP library
    USE prtctl          ! Print control
@@ -28,8 +30,13 @@
    PUBLIC   dyn_keg    ! routine called by step module
    
+   INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2  = 0   !: 2nd order centered scheme (standard scheme)
+   INTEGER, PARAMETER, PUBLIC  ::   nkeg_HW  = 1   !: Hollingsworth et al., QJRMS, 1983
+   !
+   REAL(wp) ::   r1_48 = 1._wp / 48._wp   !: =1/(4*2*6)
+   
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.6 , NEMO Consortium (2015)
    !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -37,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_keg( kt )
+   SUBROUTINE dyn_keg( kt, kscheme )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_keg  ***
@@ -45,6 +52,12 @@
       !!      general momentum trend.
       !!
-      !! ** Method  :   Compute the now horizontal kinetic energy 
+      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that 
+      !!      conserve kinetic energy. Compute the now horizontal kinetic energy 
       !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
+      !!              * kscheme = nkeg_HW : Hollingsworth correction following
+      !!      Arakawa (2001). The now horizontal kinetic energy is given by:
+      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((un(j+1)+un(j-1))/2)^2  )
+      !!                    + mj-1(  2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2  ) ]
+      !!      
       !!      Take its horizontal gradient and add it to the general momentum
       !!      trend (ua,va).
@@ -54,6 +67,10 @@
       !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
       !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
+      !!
+      !! ** References : Arakawa, A., International Geophysics 2001.
+      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
+      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme 
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -63,40 +80,60 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_keg')
+      IF( nn_timing == 1 )   CALL timing_start('dyn_keg')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zhke )
+      CALL wrk_alloc( jpi,jpj,jpk,   zhke )
       !
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend'
+         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
          IF(lwp) WRITE(numout,*) '~~~~~~~'
       ENDIF
 
       IF( l_trddyn ) THEN           ! Save ua and va trends
-         CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
+         CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
          ztrdu(:,:,:) = ua(:,:,:) 
          ztrdv(:,:,:) = va(:,:,:) 
       ENDIF
       
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         DO jj = 2, jpj         ! Horizontal kinetic energy at T-point
-            DO ji = fs_2, jpi   ! vector opt.
-               zu = 0.25 * (  un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
-                  &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)  )
-               zv = 0.25 * (  vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
-                  &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)  )
-               zhke(ji,jj,jk) = zv + zu
-!!gm simplier coding  ==>>   ~ faster
-!    don't forget to suppress local zu zv scalars
-!               zhke(ji,jj,jk) = 0.25 * (   un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
-!                  &                      + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)   &
-!                  &                      + vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
-!                  &                      + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   )
-!!gm end <<==
-            END DO  
-         END DO  
-         DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends
+      zhke(:,:,jpk) = 0._wp
+      
+      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
+      !
+      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
+         DO jk = 1, jpkm1
+            DO jj = 2, jpj
+               DO ji = fs_2, jpi   ! vector opt.
+                  zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
+                     &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
+                  zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
+                     &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
+                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
+               END DO  
+            END DO
+         END DO
+         !
+      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1       
+               DO ji = fs_2, jpim1   ! vector opt.
+                  zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
+                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
+                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
+                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
+                     !
+                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
+                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
+                     &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
+                     &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
+                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
+               END DO  
+            END DO
+         END DO
+         CALL lbc_lnk( zhke, 'T', 1. )
+         !
+      END SELECT
+      !
+      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
+         DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
                ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
@@ -104,34 +141,11 @@
             END DO 
          END DO
-!!gm idea to be tested  ==>>   is it faster on scalar computers ?
-!         DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends
-!            DO ji = fs_2, fs_jpim1   ! vector opt.
-!               ua(ji,jj,jk) = ua(ji,jj,jk) - 0.25 * ( + un(ji+1,jj  ,jk) * un(ji+1,jj  ,jk)   &
-!                  &                                   + vn(ji+1,jj-1,jk) * vn(ji+1,jj-1,jk)   &
-!                  &                                   + vn(ji+1,jj  ,jk) * vn(ji+1,jj  ,jk)   &
-!                  !
-!                  &                                   - un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
-!                  &                                   - vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
-!                  &                                   - vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   ) / e1u(ji,jj)
-!                  !
-!               va(ji,jj,jk) = va(ji,jj,jk) - 0.25 * (   un(ji-1,jj+1,jk) * un(ji-1,jj+1,jk)   &
-!                  &                                   + un(ji  ,jj+1,jk) * un(ji  ,jj+1,jk)   &
-!                  &                                   + vn(ji  ,jj+1,jk) * vn(ji  ,jj+1,jk)   &
-!                  !
-!                  &                                   - un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
-!                  &                                   - un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)   &
-!                  &                                   - vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   ) / e2v(ji,jj)
-!            END DO 
-!         END DO
-!!gm en idea            <<==
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-
-      IF( l_trddyn ) THEN      ! save the Kinetic Energy trends for diagnostic
+      END DO
+      !
+      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
-         CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
+         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv )
       ENDIF
       !
@@ -139,7 +153,7 @@
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zhke )
+      CALL wrk_dealloc( jpi,jpj,jpk,   zhke )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_keg')
+      IF( nn_timing == 1 )   CALL timing_stop('dyn_keg')
       !
    END SUBROUTINE dyn_keg
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 5382)
@@ -64,4 +64,5 @@
                                                                                             ! start and count arrays
       LOGICAL                      ::   ll_found_restart
+      CHARACTER(len=256)           ::   cl_path
       CHARACTER(len=256)           ::   cl_filename
       CHARACTER(len=NF90_MAX_NAME) ::   cl_dname
@@ -70,13 +71,15 @@
       !!----------------------------------------------------------------------
 
-      ! Find a restart file
+      ! Find a restart file. Assume iceberg restarts in same directory as ocean restarts. 
+      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_filename), exist=ll_found_restart )
+         INQUIRE( file=TRIM(cl_path)//TRIM(cl_filename), exist=ll_found_restart )
       ELSE
          cl_filename = 'restart_icebergs.nc'
-         INQUIRE( file=TRIM(cl_filename), exist=ll_found_restart )
+         INQUIRE( file=TRIM(cl_path)//TRIM(cl_filename), exist=ll_found_restart )
       ENDIF
 
@@ -86,7 +89,7 @@
 
       IF (nn_verbose_level >= 0 .AND. lwp)  &
-         WRITE(numout,'(2a)') 'icebergs, read_restart_bergs: found restart file = ',TRIM(cl_filename)
-
-      nret = NF90_OPEN(TRIM(cl_filename), NF90_NOWRITE, ncid)
+         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')
 
@@ -228,4 +231,5 @@
       INTEGER ::   jn   ! dummy loop index
       INTEGER ::   ix_dim, iy_dim, ik_dim, in_dim
+      CHARACTER(len=256)     :: cl_path
       CHARACTER(len=256)     :: cl_filename
       TYPE(iceberg), POINTER :: this
@@ -233,12 +237,15 @@
       !!----------------------------------------------------------------------
 
+      ! Assume we write iceberg restarts to same directory as ocean restarts.
+      cl_path = TRIM(cn_ocerst_outdir)
+      IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/'
       IF( lk_mpp ) THEN
-         WRITE(cl_filename,'("icebergs_",I8.8,"_restart_",I4.4,".nc")') kt, narea-1
+         WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart_",I4.4,".nc")') TRIM(cexper), kt, narea-1
       ELSE
-         WRITE(cl_filename,'("icebergs_",I8.8,"_restart.nc")') kt
-      ENDIF
-      IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM(cl_filename)
-
-      nret = NF90_CREATE(TRIM(cl_filename), NF90_CLOBBER, ncid)
+         WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart.nc")') TRIM(cexper), kt
+      ENDIF
+      IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM(cl_path)//TRIM(cl_filename)
+
+      nret = NF90_CREATE(TRIM(cl_path)//TRIM(cl_filename), NF90_CLOBBER, ncid)
       IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_create failed')
 
@@ -372,5 +379,5 @@
          ENDIF
       ENDDO
-      IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_filename),' var: stored_ice  written'
+      IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_path)//TRIM(cl_filename),' var: stored_ice  written'
 
       nret = NF90_PUT_VAR( ncid, nkountid, num_bergs(:) )
@@ -379,5 +386,5 @@
       nret = NF90_PUT_VAR( ncid, nsheatid, berg_grid%stored_heat(:,:) )
       IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var stored_heat failed')
-      IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_filename),' var: stored_heat written'
+      IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_path)//TRIM(cl_filename),' var: stored_heat written'
 
       nret = NF90_PUT_VAR( ncid, ncalvid , src_calving(:,:) )
@@ -385,5 +392,5 @@
       nret = NF90_PUT_VAR( ncid, ncalvhid, src_calving_hflx(:,:) )
       IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, write_restart: nf_put_var calving_hflx failed')
-      IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_filename),' var: calving written'
+      IF( lwp ) WRITE(numout,*) 'file: ',TRIM(cl_path)//TRIM(cl_filename),' var: calving written'
 
       IF ( ASSOCIATED(first_berg) ) THEN
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 5382)
@@ -26,6 +26,9 @@
    CHARACTER(lc) ::   cn_exp           !: experiment name used for output filename
    CHARACTER(lc) ::   cn_ocerst_in     !: suffix of ocean restart name (input)
+   CHARACTER(lc) ::   cn_ocerst_indir  !: restart input directory
    CHARACTER(lc) ::   cn_ocerst_out    !: suffix of ocean restart name (output)
+   CHARACTER(lc) ::   cn_ocerst_outdir !: restart output directory
    LOGICAL       ::   ln_rstart        !: start from (F) rest or (T) a restart file
+   LOGICAL       ::   ln_rst_list      !: output restarts at list of times (T) or by frequency (F)
    INTEGER       ::   nn_no            !: job number
    INTEGER       ::   nn_rstctl        !: control of the time step (0, 1 or 2)
@@ -38,7 +41,9 @@
    INTEGER       ::   nn_write         !: model standard output frequency
    INTEGER       ::   nn_stock         !: restart file frequency
+   INTEGER, DIMENSION(10) :: nn_stocklist  !: restart dump times
    LOGICAL       ::   ln_dimgnnn       !: type of dimgout. (F): 1 file for all proc
                                                        !:                  (T): 1 file per proc
    LOGICAL       ::   ln_mskland       !: mask land points in NetCDF outputs (costly: + ~15%)
+   LOGICAL       ::   ln_cfmeta        !: output additional data to netCDF files required for compliance with the CF metadata standard
    LOGICAL       ::   ln_clobber       !: clobber (overwrite) an existing file
    INTEGER       ::   nn_chunksz       !: chunksize (bytes) for NetCDF file (works only with iom_nf90 routines)
@@ -78,4 +83,5 @@
    INTEGER       ::   nwrite                      !: model standard output frequency
    INTEGER       ::   nstock                      !: restart file frequency
+   INTEGER, DIMENSION(10) :: nstocklist           !: restart dump times
 
    !!----------------------------------------------------------------------
@@ -85,4 +91,5 @@
    LOGICAL ::   lrst_oce              !: logical to control the oce restart write 
    INTEGER ::   numror, numrow        !: logical unit for cean restart (read and write)
+   INTEGER ::   nrst_lst              !: number of restart to output next
 
    !!----------------------------------------------------------------------
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 5382)
@@ -61,5 +61,5 @@
 #if defined key_iomput
    PRIVATE iom_set_domain_attr, iom_set_axis_attr, iom_set_field_attr, iom_set_file_attr, iom_get_file_attr, iom_set_grid_attr
-   PRIVATE set_grid, set_scalar, set_xmlatt, set_mooring, iom_update_file_name, iom_sdate
+   PRIVATE set_grid, set_grid_bounds, set_scalar, set_xmlatt, set_mooring, iom_update_file_name, iom_sdate
 # endif
 
@@ -139,4 +139,10 @@
       ENDIF
 
+      IF ( ln_cfmeta ) THEN   ! add cell bounds
+         CALL set_grid_bounds( "T", cdname ) 
+         CALL set_grid_bounds( "U", cdname )
+         CALL set_grid_bounds( "V", cdname )
+         CALL set_grid_bounds( "W", cdname )
+      ENDIF
 
       ! vertical grid definition
@@ -152,4 +158,6 @@
 #endif
       CALL iom_set_axis_attr( "icbcla", class_num )
+      CALL iom_set_axis_attr( "iax_20C", (/ REAL(20,wp) /) )
+      CALL iom_set_axis_attr( "iax_28C", (/ REAL(28,wp) /) )
       
       ! automatic definitions of some of the xml attributs
@@ -1107,11 +1115,13 @@
 
    SUBROUTINE iom_set_domain_attr( cdid, ni_glo, nj_glo, ibegin, jbegin, ni, nj, zoom_ibegin, zoom_jbegin, zoom_ni, zoom_nj,   &
-      &                                    data_dim, data_ibegin, data_ni, data_jbegin, data_nj, lonvalue, latvalue, mask )
-      CHARACTER(LEN=*)                 , INTENT(in) ::   cdid
-      INTEGER                , OPTIONAL, INTENT(in) ::   ni_glo, nj_glo, ibegin, jbegin, ni, nj
-      INTEGER                , OPTIONAL, INTENT(in) ::   data_dim, data_ibegin, data_ni, data_jbegin, data_nj
-      INTEGER                , OPTIONAL, INTENT(in) ::   zoom_ibegin, zoom_jbegin, zoom_ni, zoom_nj
-      REAL(wp), DIMENSION(:) , OPTIONAL, INTENT(in) ::   lonvalue, latvalue
-      LOGICAL, DIMENSION(:,:), OPTIONAL, INTENT(in) ::   mask
+      &                                    data_dim, data_ibegin, data_ni, data_jbegin, data_nj, lonvalue, latvalue, mask,     &
+      &                                    nvertex, bounds_lon, bounds_lat )
+      CHARACTER(LEN=*)                   , INTENT(in) ::   cdid
+      INTEGER                  , OPTIONAL, INTENT(in) ::   ni_glo, nj_glo, ibegin, jbegin, ni, nj
+      INTEGER                  , OPTIONAL, INTENT(in) ::   data_dim, data_ibegin, data_ni, data_jbegin, data_nj
+      INTEGER                  , OPTIONAL, INTENT(in) ::   zoom_ibegin, zoom_jbegin, zoom_ni, zoom_nj, nvertex
+      REAL(wp), DIMENSION(:)   , OPTIONAL, INTENT(in) ::   lonvalue, latvalue
+      REAL(wp), DIMENSION(:,:) , OPTIONAL, INTENT(in) ::   bounds_lon, bounds_lat
+      LOGICAL,  DIMENSION(:,:) , OPTIONAL, INTENT(in) ::   mask
 
       IF ( xios_is_valid_domain     (cdid) ) THEN
@@ -1119,5 +1129,5 @@
             &    data_dim=data_dim, data_ibegin=data_ibegin, data_ni=data_ni, data_jbegin=data_jbegin, data_nj=data_nj ,   &
             &    zoom_ibegin=zoom_ibegin, zoom_jbegin=zoom_jbegin, zoom_ni=zoom_ni, zoom_nj=zoom_nj,                       &
-            &    lonvalue=lonvalue, latvalue=latvalue,mask=mask )
+            &    lonvalue=lonvalue, latvalue=latvalue, mask=mask, nvertex=nvertex, bounds_lon=bounds_lon, bounds_lat=bounds_lat )
       ENDIF
 
@@ -1126,5 +1136,5 @@
             &    data_dim=data_dim, data_ibegin=data_ibegin, data_ni=data_ni, data_jbegin=data_jbegin, data_nj=data_nj ,   &
             &    zoom_ibegin=zoom_ibegin, zoom_jbegin=zoom_jbegin, zoom_ni=zoom_ni, zoom_nj=zoom_nj,                       &
-            &    lonvalue=lonvalue, latvalue=latvalue,mask=mask )
+            &    lonvalue=lonvalue, latvalue=latvalue, mask=mask, nvertex=nvertex, bounds_lon=bounds_lon, bounds_lat=bounds_lat )
       ENDIF
       CALL xios_solve_inheritance()
@@ -1251,4 +1261,125 @@
       
    END SUBROUTINE set_grid
+
+
+   SUBROUTINE set_grid_bounds( cdgrd, cdname )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE set_grid_bounds  ***
+      !!
+      !! ** Purpose :   define horizontal grid corners
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(LEN=1) , INTENT(in) :: cdgrd
+      CHARACTER(LEN=*) , INTENT(in) :: cdname
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)   :: z_bnds      ! Lat/lon coordinates of the vertices of cell (i,j)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)     :: z_cnr       ! Lat/lon coordinates of a contiguous vertex of cell (i,j)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)     :: z_pnt       ! Lat/lon coordinates of the point of cell (i,j)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)       :: z_fld       ! Working array to determine where to rotate cells
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)       :: z_rot       ! Lat/lon working array for rotation of cells
+      !
+      INTEGER :: icnr, jcnr                                      ! Offset such that the vertex coordinate (i+icnr,j+jcnr)
+      !                                                          ! represents the bottom-left corner of cell (i,j)
+      INTEGER :: ji, jj, jn, ni, nj
+
+      ALLOCATE( z_bnds(4,jpi,jpj,2), z_cnr(jpi,jpj,2), z_pnt(jpi,jpj,2), z_fld(jpi,jpj), z_rot(4,2)  )
+
+      ! Grid information
+      SELECT CASE ( TRIM(cdgrd) )
+         CASE ('T', 'W')
+            icnr = -1 ; jcnr = -1
+            IF( TRIM(cdname) == "nemo_crs" ) THEN
+               z_cnr(:,:,1) = gphif_crs ; z_cnr(:,:,2) = glamf_crs
+               z_pnt(:,:,1) = gphit_crs ; z_pnt(:,:,2) = glamt_crs
+            ELSE
+               z_cnr(:,:,1) = gphif ; z_cnr(:,:,2) = glamf
+               z_pnt(:,:,1) = gphit ; z_pnt(:,:,2) = glamt
+            ENDIF
+         CASE ('U')
+            icnr =  0 ; jcnr = -1
+            IF( TRIM(cdname) == "nemo_crs" ) THEN
+               z_cnr(:,:,1) = gphiv_crs ; z_cnr(:,:,2) = glamv_crs
+               z_pnt(:,:,1) = gphiu_crs ; z_pnt(:,:,2) = glamu_crs
+            ELSE
+               z_cnr(:,:,1) = gphiv ; z_cnr(:,:,2) = glamv
+               z_pnt(:,:,1) = gphiu ; z_pnt(:,:,2) = glamu
+            ENDIF
+         CASE ('V')
+            icnr = -1 ; jcnr =  0
+            IF( TRIM(cdname) == "nemo_crs" ) THEN
+               z_cnr(:,:,1) = gphiu_crs ; z_cnr(:,:,2) = glamu_crs
+               z_pnt(:,:,1) = gphiv_crs ; z_pnt(:,:,2) = glamv_crs
+            ELSE
+               z_cnr(:,:,1) = gphiu ; z_cnr(:,:,2) = glamu
+               z_pnt(:,:,1) = gphiv ; z_pnt(:,:,2) = glamv
+            ENDIF
+      END SELECT
+
+      ni = nlei-nldi+1 ; nj = nlej-nldj+1  ! Dimensions of subdomain interior
+
+      z_fld(:,:) = 1._wp
+      CALL lbc_lnk( z_fld, cdgrd, -1. )    ! Working array for location of northfold
+
+      ! Cell vertices that can be defined
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            z_bnds(1,ji,jj,:) = z_cnr(ji+icnr,  jj+jcnr  ,:) ! Bottom-left
+            z_bnds(2,ji,jj,:) = z_cnr(ji+icnr+1,jj+jcnr  ,:) ! Bottom-right
+            z_bnds(3,ji,jj,:) = z_cnr(ji+icnr+1,jj+jcnr+1,:) ! Top-right
+            z_bnds(4,ji,jj,:) = z_cnr(ji+icnr,  jj+jcnr+1,:) ! Top-left
+         END DO
+      END DO
+
+      ! Cell vertices on boundries
+      DO jn = 1, 4
+         CALL lbc_lnk( z_bnds(jn,:,:,1), cdgrd, 1., pval=999._wp )
+         CALL lbc_lnk( z_bnds(jn,:,:,2), cdgrd, 1., pval=999._wp )
+      END DO
+
+      ! Zero-size cells at closed boundaries
+      IF( (nbondi == -1 .OR. nbondi == 2) .AND. .NOT. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6) ) THEN
+         DO jn = 1, 4  ;  z_bnds(jn,1,:,:) = z_pnt(1,:,:)        ;  END DO      ! (West or jpni = 1), closed E-W
+      ENDIF
+      IF( (nbondi == 1 .OR. nbondi == 2) .AND. .NOT. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6) ) THEN
+         DO jn = 1, 4  ;  z_bnds(jn,nlci,:,:) = z_pnt(nlci,:,:)  ;  END DO      ! (East or jpni = 1), closed E-W
+      ENDIF
+      IF( nbondj == -1 .OR. (nbondj == 2 .AND. jperio /= 2) ) THEN
+         DO jn = 1, 4  ;  z_bnds(jn,:,1,:) = z_pnt(:,1,:)        ;  END DO      ! South or (jpnj = 1, not symmetric)
+      ENDIF
+      IF( (nbondj == 1 .OR. nbondj == 2) .AND. jperio  < 3 ) THEN
+         DO jn = 1, 4  ;  z_bnds(jn,:,nlcj,:) = z_pnt(:,nlcj,:)  ;  END DO      ! (North or jpnj = 1), no north fold
+      ENDIF
+
+! =====================================================================================================
+!             Do we need to set zero-size cells at Mediterranean / Persian Gulf region?
+! =====================================================================================================
+
+      ! Rotate cells at the north fold
+      IF( (nbondj == 1 .OR. nbondj == 2) .AND. jperio >= 3 ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( z_fld(ji,jj) == -1. ) THEN
+                  z_rot(1,:) = z_bnds(3,ji,jj,:) ; z_rot(2,:) = z_bnds(4,ji,jj,:)
+                  z_rot(3,:) = z_bnds(1,ji,jj,:) ; z_rot(4,:) = z_bnds(2,ji,jj,:)
+                  z_bnds(:,ji,jj,:) = z_rot(:,:)
+               ENDIF
+            END DO
+         END DO
+
+      ! Invert cells at the symmetric equator
+      ELSE IF( nbondj == 2 .AND. jperio == 2 ) THEN
+         DO ji = 1, jpi
+            z_rot(1:2,:) = z_bnds(3:4,ji,1,:)
+            z_rot(3:4,:) = z_bnds(1:2,ji,1,:)
+            z_bnds(:,ji,1,:) = z_rot(:,:)
+         END DO
+      ENDIF
+
+      CALL iom_set_domain_attr("grid_"//cdgrd, bounds_lat = RESHAPE(z_bnds(:,nldi:nlei,nldj:nlej,1),(/ 4,ni*nj /)),           &
+                                               bounds_lon = RESHAPE(z_bnds(:,nldi:nlei,nldj:nlej,2),(/ 4,ni*nj /)), nvertex=4 )
+
+      DEALLOCATE( z_bnds, z_cnr, z_pnt, z_fld, z_rot ) 
+
+   END SUBROUTINE set_grid_bounds
 
 
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90	(revision 5382)
@@ -61,6 +61,6 @@
       INTEGER, DIMENSION(2,5), INTENT(in   ), OPTIONAL ::   kdompar     ! domain parameters: 
 
-      CHARACTER(LEN=100) ::   clinfo           ! info character
-      CHARACTER(LEN=100) ::   cltmp            ! temporary character
+      CHARACTER(LEN=256) ::   clinfo           ! info character
+      CHARACTER(LEN=256) ::   cltmp            ! temporary character
       INTEGER            ::   iln              ! lengths of character
       INTEGER            ::   istop            ! temporary storage of nstop
@@ -393,5 +393,5 @@
       INTEGER, DIMENSION(4) :: idimsz               ! dimensions size  
       INTEGER, DIMENSION(4) :: idimid               ! dimensions id
-      CHARACTER(LEN=100)    :: clinfo               ! info character
+      CHARACTER(LEN=256)    :: clinfo               ! info character
       CHARACTER(LEN= 12), DIMENSION(4) :: cltmp     ! temporary character
       INTEGER               :: if90id               ! nf90 file identifier
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 5382)
@@ -56,12 +56,20 @@
       !!
       CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character
-      CHARACTER(LEN=50)   ::   clname   ! ice output restart file name
+      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
+      CHARACTER(lc)       ::   clpath   ! full path to ocean output restart file
       !!----------------------------------------------------------------------
       !
       IF( kt == nit000 ) THEN   ! default definitions
          lrst_oce = .FALSE.   
-         nitrst = nitend
-      ENDIF
-      IF( MOD( kt - 1, nstock ) == 0 ) THEN   
+         IF( ln_rst_list ) THEN
+            nrst_lst = 1
+            nitrst = nstocklist( nrst_lst )
+         ELSE
+            nitrst = nitend
+         ENDIF
+      ENDIF
+
+      ! frequency-based restart dumping (nn_stock)
+      IF( .NOT. ln_rst_list .AND. MOD( kt - 1, nstock ) == 0 ) THEN   
          ! we use kt - 1 and not kt - nit000 to keep the same periodicity from the beginning of the experiment
          nitrst = kt + nstock - 1                  ! define the next value of nitrst for restart writing
@@ -72,24 +80,30 @@
       ! except if we write ocean restart files every time step or if an ocean restart file was writen at nitend - 1
       IF( kt == nitrst - 1 .OR. nstock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) 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
-         ENDIF
-         ! create the file
-         clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_ocerst_out)
-         IF(lwp) THEN
-            WRITE(numout,*)
-            SELECT CASE ( jprstlib )
-            CASE ( jprstdimg )   ;   WRITE(numout,*) '             open ocean restart binary file: '//clname
-            CASE DEFAULT         ;   WRITE(numout,*) '             open ocean restart NetCDF file: '//clname
-            END SELECT
-            IF ( snc4set%luse )      WRITE(numout,*) '             opened for NetCDF4 chunking and compression'
-            IF( kt == nitrst - 1 ) THEN   ;   WRITE(numout,*) '             kt = nitrst - 1 = ', kt
-            ELSE                          ;   WRITE(numout,*) '             kt = '             , kt
+         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
             ENDIF
-         ENDIF
-         !
-         CALL iom_open( clname, numrow, ldwrt = .TRUE., kiolib = jprstlib )
-         lrst_oce = .TRUE.
+            ! create the file
+            clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_ocerst_out)
+            clpath = TRIM(cn_ocerst_outdir)
+            IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
+            IF(lwp) THEN
+               WRITE(numout,*)
+               SELECT CASE ( jprstlib )
+               CASE ( jprstdimg )   ;   WRITE(numout,*)                            &
+                   '             open ocean restart binary file: ',TRIM(clpath)//clname
+               CASE DEFAULT         ;   WRITE(numout,*)                            &
+                   '             open ocean restart NetCDF file: ',TRIM(clpath)//clname
+               END SELECT
+               IF ( snc4set%luse )      WRITE(numout,*) '             opened for NetCDF4 chunking and compression'
+               IF( kt == nitrst - 1 ) THEN   ;   WRITE(numout,*) '             kt = nitrst - 1 = ', kt
+               ELSE                          ;   WRITE(numout,*) '             kt = '             , kt
+               ENDIF
+            ENDIF
+            !
+            CALL iom_open( TRIM(clpath)//TRIM(clname), numrow, ldwrt = .TRUE., kiolib = jprstlib )
+            lrst_oce = .TRUE.
+         ENDIF
       ENDIF
       !
@@ -136,4 +150,9 @@
 !!gm  not sure what to do here   ===>>>  ask to Sebastian
          lrst_oce = .FALSE.
+            IF( ln_rst_list ) THEN
+               nrst_lst = MIN(nrst_lst + 1, SIZE(nstocklist,1))
+               nitrst = nstocklist( nrst_lst )
+            ENDIF
+            lrst_oce = .FALSE.
       ENDIF
       !
@@ -150,6 +169,7 @@
       !!                the file has already been opened
       !!----------------------------------------------------------------------
-      INTEGER  ::   jlibalt = jprstlib
-      LOGICAL  ::   llok
+      INTEGER        ::   jlibalt = jprstlib
+      LOGICAL        ::   llok
+      CHARACTER(lc)  ::   clpath   ! full path to ocean output restart file
       !!----------------------------------------------------------------------
       !
@@ -165,11 +185,13 @@
          ENDIF
 
+         clpath = TRIM(cn_ocerst_indir)
+         IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
          IF ( jprstlib == jprstdimg ) THEN
            ! eventually read netcdf file (monobloc)  for restarting on different number of processors
            ! if {cn_ocerst_in}.nc exists, then set jlibalt to jpnf90
-           INQUIRE( FILE = TRIM(cn_ocerst_in)//'.nc', EXIST = llok )
+           INQUIRE( FILE = TRIM(cn_ocerst_indir)//'/'//TRIM(cn_ocerst_in)//'.nc', EXIST = llok )
            IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF
          ENDIF
-         CALL iom_open( cn_ocerst_in, numror, kiolib = jlibalt )
+         CALL iom_open( TRIM(clpath)//cn_ocerst_in, numror, kiolib = jlibalt )
       ENDIF
    END SUBROUTINE rst_read_open
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 5382)
@@ -1440,5 +1440,5 @@
       IF( srcv(jpr_cal)%laction ) THEN 
          zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
-         CALL iom_put( 'calving', frcv(jpr_cal)%z3(:,:,1) )
+         CALL iom_put( 'calving_cea', frcv(jpr_cal)%z3(:,:,1) )
       ENDIF
 
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5382)
@@ -382,5 +382,5 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_out,   &
+      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir,  &
          &                ln_limdyn, rn_amax, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt  
       !!-------------------------------------------------------------------
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 5382)
@@ -561,6 +561,5 @@
       CALL iom_put('isfgammat', zgammat2d)
       CALL iom_put('isfgammas', zgammas2d)
-         !
-      !CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zqisf, zfwfisf  )
+      !
       CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
       !
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 5382)
@@ -47,4 +47,6 @@
    USE lbclnk         ! ocean lateral boundary conditions
    USE timing          ! Timing
+   USE stopar          ! Stochastic T/S fluctuations
+   USE stopts          ! Stochastic T/S fluctuations
 
    IMPLICIT NONE
@@ -313,7 +315,9 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m]
       !
-      INTEGER  ::   ji, jj, jk                ! dummy loop indices
-      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars
-      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
+      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices
+      INTEGER  ::   jdof
+      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars
+      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      -
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors
       !!----------------------------------------------------------------------
       !
@@ -324,43 +328,108 @@
       CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  !
-                  zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
-                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
-                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-                  ztm = tmask(ji,jj,jk)                                         ! tmask
-                  !
-                  zn3 = EOS013*zt   &
-                     &   + EOS103*zs+EOS003
-                     !
-                  zn2 = (EOS022*zt   &
-                     &   + EOS112*zs+EOS012)*zt   &
-                     &   + (EOS202*zs+EOS102)*zs+EOS002
-                     !
-                  zn1 = (((EOS041*zt   &
-                     &   + EOS131*zs+EOS031)*zt   &
-                     &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
-                     &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
-                     &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
-                     !
-                  zn0 = (((((EOS060*zt   &
-                     &   + EOS150*zs+EOS050)*zt   &
-                     &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
-                     &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
-                     &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
-                     &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
-                     &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
-                     !
-                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
-                  !
-                  prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface
-                  !
-                  prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked)
+         ! Stochastic equation of state
+         IF ( ln_sto_eos ) THEN
+            ALLOCATE(zn0_sto(1:2*nn_sto_eos))
+            ALLOCATE(zn_sto(1:2*nn_sto_eos))
+            ALLOCATE(zsign(1:2*nn_sto_eos))
+            DO jsmp = 1, 2*nn_sto_eos, 2
+              zsign(jsmp)   = 1._wp
+              zsign(jsmp+1) = -1._wp
+            END DO
+            !
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     !
+                     ! compute density (2*nn_sto_eos) times:
+                     ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts)
+                     ! (2) for t-dt, s-ds (with the opposite fluctuation)
+                     DO jsmp = 1, nn_sto_eos*2
+                        jdof   = (jsmp + 1) / 2
+                        zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth
+                        zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature
+                        zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp)
+                        zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity
+                        ztm    = tmask(ji,jj,jk)                                         ! tmask
+                        !
+                        zn3 = EOS013*zt   &
+                           &   + EOS103*zs+EOS003
+                           !
+                        zn2 = (EOS022*zt   &
+                           &   + EOS112*zs+EOS012)*zt   &
+                           &   + (EOS202*zs+EOS102)*zs+EOS002
+                           !
+                        zn1 = (((EOS041*zt   &
+                           &   + EOS131*zs+EOS031)*zt   &
+                           &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
+                           &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
+                           &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+                           !
+                        zn0_sto(jsmp) = (((((EOS060*zt   &
+                           &   + EOS150*zs+EOS050)*zt   &
+                           &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+                           &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+                           &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+                           &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+                           &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+                           !
+                        zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp)
+                     END DO
+                     !
+                     ! compute stochastic density as the mean of the (2*nn_sto_eos) densities
+                     prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp
+                     DO jsmp = 1, nn_sto_eos*2
+                        prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface
+                        !
+                        prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rau0 - 1._wp  )   ! density anomaly (masked)
+                     END DO
+                     prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos
+                     prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos
+                  END DO
                END DO
             END DO
-         END DO
-         !
+            DEALLOCATE(zn0_sto,zn_sto,zsign)
+         ! Non-stochastic equation of state
+         ELSE
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     !
+                     zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
+                     zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
+                     zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+                     ztm = tmask(ji,jj,jk)                                         ! tmask
+                     !
+                     zn3 = EOS013*zt   &
+                        &   + EOS103*zs+EOS003
+                        !
+                     zn2 = (EOS022*zt   &
+                        &   + EOS112*zs+EOS012)*zt   &
+                        &   + (EOS202*zs+EOS102)*zs+EOS002
+                        !
+                     zn1 = (((EOS041*zt   &
+                        &   + EOS131*zs+EOS031)*zt   &
+                        &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
+                        &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
+                        &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+                        !
+                     zn0 = (((((EOS060*zt   &
+                        &   + EOS150*zs+EOS050)*zt   &
+                        &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+                        &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+                        &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+                        &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+                        &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+                        !
+                     zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+                     !
+                     prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface
+                     !
+                     prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked)
+                  END DO
+               END DO
+            END DO
+         ENDIF
+         
       CASE( 1 )                !==  simplified EOS  ==!
          !
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl_rst.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl_rst.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl_rst.F90	(revision 5382)
@@ -27,5 +27,5 @@
    !!---------------------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!---------------------------------------------------------------------------------
@@ -43,5 +43,6 @@
       INTEGER ::   jk                 ! loop indice
       CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character
-      CHARACTER(LEN=50)   ::   clname   ! ice output restart file name
+      CHARACTER(LEN=50)   ::   clname   ! output restart file name
+      CHARACTER(LEN=256)  ::   clpath   ! full path to restart file
       !!--------------------------------------------------------------------------------
 
@@ -56,4 +57,6 @@
          ! create the file
          clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_trdrst_out)
+         clpath = TRIM(cn_ocerst_outdir)
+         IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
          IF(lwp) THEN
             WRITE(numout,*)
@@ -67,5 +70,5 @@
          ENDIF
 
-         CALL iom_open( clname, nummxlw, ldwrt = .TRUE., kiolib = jprstlib )
+         CALL iom_open( TRIM(clpath)//TRIM(clname), nummxlw, ldwrt = .TRUE., kiolib = jprstlib )
       ENDIF
 
@@ -133,4 +136,5 @@
       INTEGER ::   jlibalt = jprstlib
       LOGICAL ::   llok
+      CHARACTER(LEN=256)  ::   clpath   ! full path to restart file
       !!-----------------------------------------------------------------------------
 
@@ -140,8 +144,12 @@
          WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
       ENDIF
+
+      clpath = TRIM(cn_ocerst_indir)
+      IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
+
       IF ( jprstlib == jprstdimg ) THEN
          ! eventually read netcdf file (monobloc)  for restarting on different number of processors
          ! if {cn_trdrst_in}.nc exists, then set jlibalt to jpnf90
-         INQUIRE( FILE = TRIM(cn_trdrst_in)//'.nc', EXIST = llok )
+         INQUIRE( FILE = TRIM(clpath)//TRIM(cn_trdrst_in)//'.nc', EXIST = llok )
          IF ( llok ) THEN   ;   jlibalt = jpnf90   
          ELSE               ;   jlibalt = jprstlib   
@@ -149,5 +157,5 @@
       ENDIF
 
-      CALL iom_open( cn_trdrst_in, inum, kiolib = jlibalt ) 
+      CALL iom_open( TRIM(clpath)//TRIM(cn_trdrst_in), inum, kiolib = jlibalt ) 
 
       IF( ln_trdmxl_instant ) THEN 
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 5382)
@@ -171,10 +171,12 @@
             END DO
          END DO
+         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition
+
          IF ( ln_isfcav ) THEN
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
                   ! (ISF) ========================================================================
-                  ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points 
-                  ikbv = mikv(ji,jj)         ! (deepest ocean u- and v-points)
+                  ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
+                  ikbv = mikv(ji,jj)         ! (1st wet ocean u- and v-points)
                   !
                   zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     &
@@ -183,6 +185,6 @@
                      &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  )
               !
-                  zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 )
-                  zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 )
+                  zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_tfeb2 )
+                  zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_tfeb2 )
               !
                   tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu * (1._wp - umask(ji,jj,1))
@@ -202,7 +204,7 @@
                END DO
             END DO
-         END IF
-         !
-         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition
+            CALL lbc_lnk( tfrua, 'U', 1. )   ;   CALL lbc_lnk( tfrva, 'V', 1. )      ! Lateral boundary condition
+         END IF
+         !
          !
          IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        &
@@ -277,9 +279,11 @@
             IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien
          ENDIF
-         IF(lwp) WRITE(numout,*) '      top    friction coef.   rn_bfri1  = ', rn_bfri1
-         IF( ln_tfr2d ) THEN
-            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d
-            IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien
-         ENDIF
+         IF ( ln_isfcav ) THEN
+            IF(lwp) WRITE(numout,*) '      top    friction coef.   rn_bfri1  = ', rn_tfri1
+            IF( ln_tfr2d ) THEN
+               IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d
+               IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien
+            ENDIF
+         END IF
          !
          IF(ln_bfr2d) THEN
@@ -295,4 +299,19 @@
          bfrua(:,:) = - bfrcoef2d(:,:)
          bfrva(:,:) = - bfrcoef2d(:,:)
+         !
+         IF ( ln_isfcav ) THEN
+            IF(ln_tfr2d) THEN
+               ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
+               CALL iom_open('tfr_coef.nc',inum)
+               CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array
+               CALL iom_close(inum)
+               tfrcoef2d(:,:) = rn_tfri1 * ( 1 + rn_tfrien * tfrcoef2d(:,:) )
+            ELSE
+               tfrcoef2d(:,:) = rn_tfri1  ! initialize tfrcoef2d to the namelist variable
+            ENDIF
+            !
+            tfrua(:,:) = - tfrcoef2d(:,:)
+            tfrva(:,:) = - tfrcoef2d(:,:)
+         END IF
          !
       CASE( 2 )
@@ -311,18 +330,20 @@
             IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien
          ENDIF
-         IF(lwp) WRITE(numout,*) '      quadratic top    friction'
-         IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri2  = ', rn_tfri2
-         IF(lwp) WRITE(numout,*) '      Max. coef. (log case)   rn_tfri2_max  = ', rn_tfri2_max
-         IF(lwp) WRITE(numout,*) '      background tke   rn_tfeb2  = ', rn_tfeb2
-         IF(lwp) WRITE(numout,*) '      log formulation   ln_tfr2d = ', ln_loglayer
-         IF(lwp) WRITE(numout,*) '      bottom roughness  rn_tfrz0 [m] = ', rn_tfrz0
-         IF( rn_tfrz0<=0.e0 ) THEN
-            WRITE(ctmp1,*) '      bottom roughness must be strictly positive'
-            CALL ctl_stop( ctmp1 )
-         ENDIF
-         IF( ln_tfr2d ) THEN
-            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d
-            IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien
-         ENDIF
+         IF ( ln_isfcav ) THEN
+            IF(lwp) WRITE(numout,*) '      quadratic top    friction'
+            IF(lwp) WRITE(numout,*) '      friction coef.    rn_tfri2     = ', rn_tfri2
+            IF(lwp) WRITE(numout,*) '      Max. coef. (log case)   rn_tfri2_max  = ', rn_tfri2_max
+            IF(lwp) WRITE(numout,*) '      background tke    rn_tfeb2     = ', rn_tfeb2
+            IF(lwp) WRITE(numout,*) '      log formulation   ln_tfr2d     = ', ln_loglayer
+            IF(lwp) WRITE(numout,*) '      top roughness     rn_tfrz0 [m] = ', rn_tfrz0
+            IF( rn_tfrz0<=0.e0 ) THEN
+               WRITE(ctmp1,*) '      top roughness must be strictly positive'
+               CALL ctl_stop( ctmp1 )
+            ENDIF
+            IF( ln_tfr2d ) THEN
+               IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d
+               IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien
+            ENDIF
+         END IF
          !
          IF(ln_bfr2d) THEN
@@ -336,4 +357,17 @@
             bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable
          ENDIF
+         
+         IF ( ln_isfcav ) THEN
+            IF(ln_tfr2d) THEN
+               ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
+               CALL iom_open('tfr_coef.nc',inum)
+               CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array
+               CALL iom_close(inum)
+               !
+               tfrcoef2d(:,:) = rn_tfri2 * ( 1 + rn_tfrien * tfrcoef2d(:,:) )
+            ELSE
+               tfrcoef2d(:,:) = rn_tfri2  ! initialize tfrcoef2d to the namelist variable
+            ENDIF
+         END IF
          !
          IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all
@@ -346,4 +380,14 @@
                END DO
             END DO
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ikbt = mikt(ji,jj)
+                     ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_tfrz0 ))**2._wp
+                     tfrcoef2d(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp)
+                     tfrcoef2d(ji,jj) = MIN(tfrcoef2d(ji,jj), rn_tfri2_max)
+                  END DO
+               END DO
+            END IF
          ENDIF
          !
@@ -398,4 +442,28 @@
              zminbfr = MIN(  zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) )  )
              zmaxbfr = MAX(  zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) )  )
+! (ISF)
+             IF ( ln_isfcav ) THEN
+                ikbu = miku(ji,jj)       ! 1st wet ocean level at u- and v-points
+                ikbv = mikv(ji,jj)
+                zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt
+                zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt
+                IF( ABS( tfrcoef2d(ji,jj) ) > zfru ) THEN
+                   IF( ln_ctl ) THEN
+                      WRITE(numout,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbu
+                      WRITE(numout,*) 'TFR ', ABS( tfrcoef2d(ji,jj) ), zfru
+                   ENDIF
+                   ictu = ictu + 1
+                ENDIF
+                IF( ABS( tfrcoef2d(ji,jj) ) > zfrv ) THEN
+                   IF( ln_ctl ) THEN
+                      WRITE(numout,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbv
+                      WRITE(numout,*) 'TFR ', tfrcoef2d(ji,jj), zfrv
+                   ENDIF
+                   ictv = ictv + 1
+                ENDIF
+                zmintfr = MIN(  zmintfr, MIN( zfru, ABS( tfrcoef2d(ji,jj) ) )  )
+                zmaxtfr = MAX(  zmaxtfr, MIN( zfrv, ABS( tfrcoef2d(ji,jj) ) )  )
+             END IF
+! END ISF
          END DO
       END DO
@@ -405,12 +473,14 @@
          CALL mpp_min( zminbfr )
          CALL mpp_max( zmaxbfr )
+         IF ( ln_isfcav) CALL mpp_min( zmintfr )
+         IF ( ln_isfcav) CALL mpp_max( zmaxtfr )
       ENDIF
       IF( .NOT.ln_bfrimp) THEN
       IF( lwp .AND. ictu + ictv > 0 ) THEN
-         WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points '
-         WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points '
+         WRITE(numout,*) ' Bottom/Top friction stability check failed at ', ictu, ' U-points '
+         WRITE(numout,*) ' Bottom/Top friction stability check failed at ', ictv, ' V-points '
          WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr
-         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr
-         WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary'
+         IF ( ln_isfcav ) WRITE(numout,*) ' Top friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr
+         WRITE(numout,*) ' Bottom/Top friction coefficient will be reduced where necessary'
       ENDIF
       ENDIF
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5382)
@@ -83,4 +83,6 @@
    USE lbcnfd, ONLY: isendto, nsndto, nfsloop, nfeloop ! Setup of north fold exchanges 
    USE sbc_oce, ONLY: lk_oasis
+   USE stopar
+   USE stopts
 
    IMPLICIT NONE
@@ -437,4 +439,6 @@
       IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_init       ! Cross Land Advection
                             CALL icb_init( rdt, nit000)   ! initialise icebergs instance
+                            CALL sto_par_init   ! Stochastic parametrization
+      IF( ln_sto_eos     )  CALL sto_pts_init   ! RRandom T/S fluctuations
      
 #if defined key_top
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5382)
@@ -107,4 +107,9 @@
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update stochastic parameters and random T/S fluctuations
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+                        CALL sto_par( kstp )          ! Stochastic parameters
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Ocean physics update                (ua, va, tsa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -146,4 +151,5 @@
       !
       IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
+         IF(ln_sto_eos ) CALL sto_pts( tsn )          ! Random T/S fluctuations
                          CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density
          IF( ln_zps .AND. .NOT. ln_isfcav)                               &
@@ -181,4 +187,5 @@
           ! Note that the computation of vertical velocity above, hence "after" sea level
           ! is necessary to compute momentum advection for the rhs of barotropic loop:
+            IF(ln_sto_eos ) CALL sto_pts( tsn )                             ! Random T/S fluctuations
                             CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation
             IF( ln_zps .AND. .NOT. ln_isfcav)                               &
@@ -262,4 +269,5 @@
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
+            IF( ln_sto_eos ) CALL sto_pts( tsn )                 ! Random T/S fluctuations
                              CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation
             IF( ln_zps .AND. .NOT. ln_isfcav)                                &
@@ -272,4 +280,5 @@
       ELSE                                                  ! centered hpg  (eos then time stepping)
          IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case
+            IF( ln_sto_eos ) CALL sto_pts( tsn )    ! Random T/S fluctuations
                              CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation
          IF( ln_zps .AND. .NOT. ln_isfcav)                                   &
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5382)
@@ -53,4 +53,7 @@
 
    USE dynnxt           ! time-stepping                    (dyn_nxt routine)
+
+   USE stopar           ! Stochastic parametrization       (sto_par routine)
+   USE stopts 
 
    USE bdy_par          ! for lk_bdy
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/stopar.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/stopar.F90	(revision 5382)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/stopar.F90	(revision 5382)
@@ -0,0 +1,909 @@
+MODULE stopar
+   !!======================================================================
+   !!                       ***  MODULE  stopar  ***
+   !! Stochastic parameters : definition and time stepping
+   !!=====================================================================
+   !! History :  3.3  ! 2011-10 (J.-M. Brankart)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sto_par       : update the stochastic parameters
+   !!   sto_par_init  : define the stochastic parameterization
+   !!   sto_rst_read  : read restart file for stochastic parameters
+   !!   sto_rst_write : write restart file for stochastic parameters
+   !!   sto_par_white : fill input array with white Gaussian noise
+   !!   sto_par_flt   : apply horizontal Laplacian filter to input array
+   !!----------------------------------------------------------------------
+   USE storng          ! random number generator (external module)
+   USE par_oce         ! ocean parameters
+   USE dom_oce         ! ocean space and time domain variables
+   USE lbclnk          ! lateral boundary conditions (or mpp link)
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O module
+   USE lib_mpp
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sto_par_init    ! called by nemogcm.F90
+   PUBLIC   sto_par         ! called by step.F90
+   PUBLIC   sto_rst_write   ! called by step.F90
+
+   LOGICAL           :: ln_rststo = .FALSE.  ! restart stochastic parameters from restart file
+   LOGICAL           :: ln_rstseed = .FALSE. ! read seed of RNG from restart file
+   CHARACTER(len=32) :: cn_storst_in = "restart_sto"     ! suffix of sto restart name (input)
+   CHARACTER(len=32) :: cn_storst_out = "restart_sto"    ! suffix of sto restart name (output)
+   INTEGER           :: numstor, numstow     ! logical unit for restart (read and write)
+
+   INTEGER           :: jpsto2d = 0          ! number of 2D stochastic parameters
+   INTEGER           :: jpsto3d = 0          ! number of 3D stochastic parameters
+
+   REAL(wp), PUBLIC, DIMENSION(:,:,:),   ALLOCATABLE :: sto2d      ! 2D stochastic parameters
+   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: sto3d      ! 3D stochastic parameters
+   REAL(wp),         DIMENSION(:,:),     ALLOCATABLE :: sto_tmp    ! temporary workspace
+   REAL(wp),         DIMENSION(:,:),     ALLOCATABLE :: sto2d_abc  ! a, b, c parameters (for 2D arrays)
+   REAL(wp),         DIMENSION(:,:),     ALLOCATABLE :: sto3d_abc  ! a, b, c parameters (for 3D arrays)
+   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto2d_ave  ! mean value (for 2D arrays)
+   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto3d_ave  ! mean value (for 3D arrays)
+   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto2d_std  ! standard deviation (for 2D arrays)
+   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto3d_std  ! standard deviation (for 3D arrays)
+   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto2d_lim  ! limitation factor (for 2D arrays)
+   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto3d_lim  ! limitation factor (for 3D arrays)
+   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto2d_tcor ! time correlation (for 2D arrays)
+   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto3d_tcor ! time correlation (for 3D arrays)
+   INTEGER,          DIMENSION(:),       ALLOCATABLE :: sto2d_ord  ! order of autoregressive process
+   INTEGER,          DIMENSION(:),       ALLOCATABLE :: sto3d_ord  ! order of autoregressive process
+
+   CHARACTER(len=1), DIMENSION(:),       ALLOCATABLE :: sto2d_typ  ! nature of grid point (T, U, V, W, F, I)
+   CHARACTER(len=1), DIMENSION(:),       ALLOCATABLE :: sto3d_typ  ! nature of grid point (T, U, V, W, F, I)
+   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto2d_sgn  ! control of the sign accross the north fold
+   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto3d_sgn  ! control of the sign accross the north fold
+   INTEGER,          DIMENSION(:),       ALLOCATABLE :: sto2d_flt  ! number of passes of Laplacian filter
+   INTEGER,          DIMENSION(:),       ALLOCATABLE :: sto3d_flt  ! number of passes of Laplacian filter
+   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto2d_fac  ! factor to restore std after filtering
+   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto3d_fac  ! factor to restore std after filtering
+
+   LOGICAL, PUBLIC :: ln_sto_ldf = .FALSE.    ! stochastic lateral diffusion
+   INTEGER, PUBLIC :: jsto_ldf                ! index of lateral diffusion stochastic parameter
+   REAL(wp)        :: rn_ldf_std              ! lateral diffusion standard deviation (in percent)
+   REAL(wp)        :: rn_ldf_tcor             ! lateral diffusion correlation timescale (in timesteps)
+
+   LOGICAL, PUBLIC :: ln_sto_hpg = .FALSE.    ! stochastic horizontal pressure gradient
+   INTEGER, PUBLIC :: jsto_hpgi               ! index of stochastic hpg parameter (i direction)
+   INTEGER, PUBLIC :: jsto_hpgj               ! index of stochastic hpg parameter (j direction)
+   REAL(wp)        :: rn_hpg_std              ! density gradient standard deviation (in percent)
+   REAL(wp)        :: rn_hpg_tcor             ! density gradient correlation timescale (in timesteps)
+
+   LOGICAL, PUBLIC :: ln_sto_pstar = .FALSE.  ! stochastic ice strength
+   INTEGER, PUBLIC :: jsto_pstar              ! index of stochastic ice strength
+   REAL(wp), PUBLIC:: rn_pstar_std            ! ice strength standard deviation (in percent)
+   REAL(wp)        :: rn_pstar_tcor           ! ice strength correlation timescale (in timesteps)
+   INTEGER         :: nn_pstar_flt = 0        ! number of passes of Laplacian filter
+   INTEGER         :: nn_pstar_ord = 1        ! order of autoregressive processes
+
+   LOGICAL, PUBLIC :: ln_sto_trd = .FALSE.    ! stochastic model trend
+   INTEGER, PUBLIC :: jsto_trd                ! index of stochastic trend parameter
+   REAL(wp)        :: rn_trd_std              ! trend standard deviation (in percent)
+   REAL(wp)        :: rn_trd_tcor             ! trend correlation timescale (in timesteps)
+
+   LOGICAL, PUBLIC :: ln_sto_eos = .FALSE.    ! stochastic equation of state
+   INTEGER, PUBLIC :: nn_sto_eos = 1          ! number of degrees of freedom in stochastic equation of state
+   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE :: jsto_eosi ! index of stochastic eos parameter (i direction)
+   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE :: jsto_eosj ! index of stochastic eos parameter (j direction)
+   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE :: jsto_eosk ! index of stochastic eos parameter (k direction)
+   REAL(wp)        :: rn_eos_stdxy            ! random walk horz. standard deviation (in grid points)
+   REAL(wp)        :: rn_eos_stdz             ! random walk vert. standard deviation (in grid points)
+   REAL(wp)        :: rn_eos_tcor             ! random walk correlation timescale (in timesteps)
+   REAL(wp)        :: rn_eos_lim = 3.0_wp     ! limitation factor
+   INTEGER         :: nn_eos_flt = 0          ! number of passes of Laplacian filter
+   INTEGER         :: nn_eos_ord = 1          ! order of autoregressive processes
+
+   LOGICAL, PUBLIC :: ln_sto_trc = .FALSE.    ! stochastic tracer dynamics
+   INTEGER, PUBLIC :: nn_sto_trc = 1          ! number of degrees of freedom in stochastic tracer dynamics
+   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE :: jsto_trci ! index of stochastic trc parameter (i direction)
+   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE :: jsto_trcj ! index of stochastic trc parameter (j direction)
+   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE :: jsto_trck ! index of stochastic trc parameter (k direction)
+   REAL(wp)        :: rn_trc_stdxy            ! random walk horz. standard deviation (in grid points)
+   REAL(wp)        :: rn_trc_stdz             ! random walk vert. standard deviation (in grid points)
+   REAL(wp)        :: rn_trc_tcor             ! random walk correlation timescale (in timesteps)
+   REAL(wp)        :: rn_trc_lim = 3.0_wp     ! limitation factor
+   INTEGER         :: nn_trc_flt = 0          ! number of passes of Laplacian filter
+   INTEGER         :: nn_trc_ord = 1          ! order of autoregressive processes
+
+   ! Public array with density correction
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: drho_ran
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: dynhpg.F90 2528 2010-12-27 17:33:53Z rblod $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE sto_par( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sto_par  ***
+      !!
+      !! ** Purpose :   update the stochastic parameters
+      !!
+      !! ** Method  :   model basic stochastic parameters
+      !!                as a first order autoregressive process AR(1),
+      !!                governed by the equation:
+      !!                   X(t) = a * X(t-1) + b * w + c
+      !!                where the parameters a, b and c are related
+      !!                to expected value, standard deviation
+      !!                and time correlation (all stationary in time) by:
+      !!                   E   [X(t)]        = c / ( 1 - a )
+      !!                   STD [X(t)]        = b / SQRT( 1 - a * a )
+      !!                   COR [X(t),X(t-k)] = a ** k
+      !!                and w is a Gaussian white noise.
+      !!
+      !!                Higher order autoregressive proces can be optionally generated
+      !!                by replacing the white noise by a lower order process.
+      !!
+      !!                1) The statistics of the stochastic parameters (X) are assumed
+      !!                constant in space (homogeneous) and time (stationary).
+      !!                This could be generalized by replacing the constant
+      !!                a, b, c parameters by functions of space and time.
+      !!
+      !!                2) The computation is performed independently for every model
+      !!                grid point, which corresponds to assume that the stochastic
+      !!                parameters are uncorrelated in space.
+      !!                This could be generalized by including a spatial filter: Y = Filt[ X ]
+      !!                (possibly non-homgeneous and non-stationary) in the computation,
+      !!                or by solving an elliptic equation: L[ Y ] = X.
+      !!
+      !!                3) The stochastic model for the parameters could also
+      !!                be generalized to depend on the current state of the ocean (not done here).
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      !!
+      INTEGER  :: ji, jj, jk, jsto, jflt
+      REAL(wp) :: stomax
+
+      !
+      ! Update 2D stochastic arrays
+      !
+      DO jsto = 1, jpsto2d
+        ! Store array from previous time step
+        sto_tmp(:,:) = sto2d(:,:,jsto)
+
+        IF ( sto2d_ord(jsto) == 1 ) THEN
+          ! Draw new random numbers from N(0,1) --> w
+          CALL sto_par_white( sto2d(:,:,jsto) )
+          ! Apply horizontal Laplacian filter to w
+          DO jflt = 1, sto2d_flt(jsto)
+            CALL lbc_lnk( sto2d(:,:,jsto), sto2d_typ(jsto), sto2d_sgn(jsto) )
+            CALL sto_par_flt( sto2d(:,:,jsto) )
+          END DO
+          ! Factor to restore standard deviation after filtering
+          sto2d(:,:,jsto) = sto2d(:,:,jsto) * sto2d_fac(jsto)
+        ELSE
+          ! Use previous process (one order lower) instead of white noise
+          sto2d(:,:,jsto) = sto2d(:,:,jsto-1)
+        ENDIF
+
+        ! Multiply white noise (or lower order process) by b --> b * w
+        sto2d(:,:,jsto) = sto2d(:,:,jsto) * sto2d_abc(jsto,2)
+        ! Update autoregressive processes --> a * X(t-1) + b * w
+        sto2d(:,:,jsto) = sto2d(:,:,jsto) + sto_tmp(:,:) * sto2d_abc(jsto,1)
+        ! Add parameter c --> a * X(t-1) + b * w + c
+        sto2d(:,:,jsto) = sto2d(:,:,jsto) + sto2d_abc(jsto,3)
+        ! Limit random parameter anomalies to std times the limitation factor
+        stomax = sto2d_std(jsto) * sto2d_lim(jsto)
+        sto2d(:,:,jsto) = sto2d(:,:,jsto) - sto2d_ave(jsto)
+        sto2d(:,:,jsto) = SIGN(MIN(stomax,ABS(sto2d(:,:,jsto))),sto2d(:,:,jsto))
+        sto2d(:,:,jsto) = sto2d(:,:,jsto) + sto2d_ave(jsto)
+
+        ! Lateral boundary conditions on sto2d
+        CALL lbc_lnk( sto2d(:,:,jsto), sto2d_typ(jsto), sto2d_sgn(jsto) )
+      END DO
+      !
+      ! Update 3D stochastic arrays
+      !
+      DO jsto = 1, jpsto3d
+         DO jk = 1, jpk
+           ! Store array from previous time step
+           sto_tmp(:,:) = sto3d(:,:,jk,jsto)
+
+           IF ( sto3d_ord(jsto) == 1 ) THEN
+             ! Draw new random numbers from N(0,1) --> w
+             CALL sto_par_white( sto3d(:,:,jk,jsto) )
+             ! Apply horizontal Laplacian filter to w
+             DO jflt = 1, sto3d_flt(jsto)
+               CALL lbc_lnk( sto3d(:,:,jk,jsto), sto3d_typ(jsto), sto3d_sgn(jsto) )
+               CALL sto_par_flt( sto3d(:,:,jk,jsto) )
+             END DO
+             ! Factor to restore standard deviation after filtering
+             sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) * sto3d_fac(jsto)
+           ELSE
+             ! Use previous process (one order lower) instead of white noise
+             sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto-1)
+           ENDIF
+
+           ! Multiply white noise by b --> b * w
+           sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) * sto3d_abc(jsto,2)
+           ! Update autoregressive processes --> a * X(t-1) + b * w
+           sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) + sto_tmp(:,:) * sto3d_abc(jsto,1)
+           ! Add parameter c --> a * X(t-1) + b * w + c
+           sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) + sto3d_abc(jsto,3)
+           ! Limit random parameters anomalies to std times the limitation factor
+           stomax = sto3d_std(jsto) * sto3d_lim(jsto)
+           sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) - sto3d_ave(jsto)
+           sto3d(:,:,jk,jsto) = SIGN(MIN(stomax,ABS(sto3d(:,:,jk,jsto))),sto3d(:,:,jk,jsto))
+           sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) + sto3d_ave(jsto)
+         END DO
+         ! Lateral boundary conditions on sto3d
+         CALL lbc_lnk( sto3d(:,:,:,jsto), sto3d_typ(jsto), sto3d_sgn(jsto) )
+      END DO
+
+   END SUBROUTINE sto_par
+
+
+   SUBROUTINE sto_par_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sto_par_init  ***
+      !!
+      !! ** Purpose :   define the stochastic parameterization
+      !!----------------------------------------------------------------------
+      NAMELIST/namsto/ ln_sto_ldf, rn_ldf_std, rn_ldf_tcor, &
+        &              ln_sto_hpg, rn_hpg_std, rn_hpg_tcor, &
+        &              ln_sto_pstar, rn_pstar_std, rn_pstar_tcor, nn_pstar_flt, nn_pstar_ord, &
+        &              ln_sto_trd, rn_trd_std, rn_trd_tcor, &
+        &              ln_sto_eos, nn_sto_eos, rn_eos_stdxy, rn_eos_stdz, &
+        &              rn_eos_tcor, nn_eos_ord, nn_eos_flt, rn_eos_lim, &
+        &              ln_sto_trc, nn_sto_trc, rn_trc_stdxy, rn_trc_stdz, &
+        &              rn_trc_tcor, nn_trc_ord, nn_trc_flt, rn_trc_lim, &
+        &              ln_rststo, ln_rstseed, cn_storst_in, cn_storst_out
+      !!----------------------------------------------------------------------
+      INTEGER :: jsto, jmem, jarea, jdof, jord, jordm1, jk, jflt
+      INTEGER(KIND=8) :: zseed1, zseed2, zseed3, zseed4
+      REAL(wp) :: rinflate
+      INTEGER  ::   ios                 ! Local integer output status for namelist read
+
+      ! Read namsto namelist : stochastic parameterization
+      REWIND( numnam_ref )              ! Namelist namdyn_adv in reference namelist : Momentum advection scheme
+      READ  ( numnam_ref, namsto, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsto in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namdyn_adv in configuration namelist : Momentum advection scheme
+      READ  ( numnam_cfg, namsto, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsto in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namsto )
+
+      !IF(ln_ens_rst_in) cn_storst_in = cn_mem//cn_storst_in
+
+      ! Parameter print
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sto_par_init : stochastic parameterization'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namsto : stochastic parameterization'
+         WRITE(numout,*) '      restart stochastic parameters           ln_rststo     = ', ln_rststo
+         WRITE(numout,*) '      read seed of RNG from restart file      ln_rstseed    = ', ln_rstseed
+         WRITE(numout,*) '      suffix of sto restart name (input)      cn_storst_in  = ', cn_storst_in
+         WRITE(numout,*) '      suffix of sto restart name (output)     cn_storst_out = ', cn_storst_out
+
+         ! WRITE(numout,*) '      stochastic lateral diffusion            ln_sto_ldf    = ', ln_sto_ldf
+         ! WRITE(numout,*) '      lateral diffusion std (in percent)      rn_ldf_std    = ', rn_ldf_std
+         ! WRITE(numout,*) '      lateral diffusion tcor (in timesteps)   rn_ldf_tcor   = ', rn_ldf_tcor
+
+         ! WRITE(numout,*) '      stochastic horizontal pressure gradient ln_sto_hpg    = ', ln_sto_hpg
+         ! WRITE(numout,*) '      density gradient std (in percent)       rn_hpg_std    = ', rn_hpg_std
+         ! WRITE(numout,*) '      density gradient tcor (in timesteps)    rn_hpg_tcor   = ', rn_hpg_tcor
+
+         ! WRITE(numout,*) '      stochastic ice strength                 ln_sto_pstar  = ', ln_sto_pstar
+         ! WRITE(numout,*) '      ice strength std (in percent)           rn_pstar_std  = ', rn_pstar_std
+         ! WRITE(numout,*) '      ice strength tcor (in timesteps)        rn_pstar_tcor = ', rn_pstar_tcor
+         ! WRITE(numout,*) '      order of autoregressive  processes      nn_pstar_ord  = ', nn_pstar_ord
+         ! WRITE(numout,*) '      passes of Laplacian filter              nn_pstar_flt  = ', nn_pstar_flt
+
+         !WRITE(numout,*) '      stochastic trend                        ln_sto_trd    = ', ln_sto_trd
+         !WRITE(numout,*) '      trend std (in percent)                  rn_trd_std    = ', rn_trd_std
+         !WRITE(numout,*) '      trend tcor (in timesteps)               rn_trd_tcor   = ', rn_trd_tcor
+
+         WRITE(numout,*) '      stochastic equation of state            ln_sto_eos    = ', ln_sto_eos
+         WRITE(numout,*) '      number of degrees of freedom            nn_sto_eos    = ', nn_sto_eos
+         WRITE(numout,*) '      random walk horz. std (in grid points)  rn_eos_stdxy  = ', rn_eos_stdxy
+         WRITE(numout,*) '      random walk vert. std (in grid points)  rn_eos_stdz   = ', rn_eos_stdz
+         WRITE(numout,*) '      random walk tcor (in timesteps)         rn_eos_tcor   = ', rn_eos_tcor
+         WRITE(numout,*) '      order of autoregressive  processes      nn_eos_ord    = ', nn_eos_ord
+         WRITE(numout,*) '      passes of Laplacian filter              nn_eos_flt    = ', nn_eos_flt
+         WRITE(numout,*) '      limitation factor                       rn_eos_lim    = ', rn_eos_lim
+
+         ! WRITE(numout,*) '      stochastic tracers dynamics             ln_sto_trc    = ', ln_sto_trc
+         ! WRITE(numout,*) '      number of degrees of freedom            nn_sto_trc    = ', nn_sto_trc
+         ! WRITE(numout,*) '      random walk horz. std (in grid points)  rn_trc_stdxy  = ', rn_trc_stdxy
+         ! WRITE(numout,*) '      random walk vert. std (in grid points)  rn_trc_stdz   = ', rn_trc_stdz
+         ! WRITE(numout,*) '      random walk tcor (in timesteps)         rn_trc_tcor   = ', rn_trc_tcor
+         ! WRITE(numout,*) '      order of autoregressive  processes      nn_trc_ord    = ', nn_trc_ord
+         ! WRITE(numout,*) '      passes of Laplacian filter              nn_trc_flt    = ', nn_trc_flt
+         ! WRITE(numout,*) '      limitation factor                       rn_trc_lim    = ', rn_trc_lim
+
+      ENDIF
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '   stochastic parameterization :'
+
+      ! Set number of 2D stochastic arrays
+      jpsto2d = 0
+      IF( ln_sto_ldf ) THEN
+         IF(lwp) WRITE(numout,*) '       - stochastic lateral diffusion'
+         jpsto2d   = jpsto2d + 1
+         jsto_ldf  = jpsto2d
+      ENDIF
+      IF( ln_sto_pstar ) THEN
+         IF(lwp) WRITE(numout,*) '       - stochastic ice strength'
+         jpsto2d    = jpsto2d + 1 * nn_pstar_ord
+         jsto_pstar = jpsto2d
+      ENDIF
+      IF( ln_sto_eos ) THEN
+         IF(lwp) WRITE(numout,*) '       - stochastic equation of state'
+         ALLOCATE(jsto_eosi(nn_sto_eos))
+         ALLOCATE(jsto_eosj(nn_sto_eos))
+         ALLOCATE(jsto_eosk(nn_sto_eos))
+         DO jdof = 1, nn_sto_eos
+            jpsto2d   = jpsto2d + 3 * nn_eos_ord
+            jsto_eosi(jdof) = jpsto2d - 2 * nn_eos_ord
+            jsto_eosj(jdof) = jpsto2d - 1 * nn_eos_ord
+            jsto_eosk(jdof) = jpsto2d
+         END DO
+      ELSE
+         nn_sto_eos = 0
+      ENDIF
+      IF( ln_sto_trc ) THEN
+         IF(lwp) WRITE(numout,*) '       - stochastic tracers dynamics'
+         ALLOCATE(jsto_trci(nn_sto_trc))
+         ALLOCATE(jsto_trcj(nn_sto_trc))
+         ALLOCATE(jsto_trck(nn_sto_trc))
+         DO jdof = 1, nn_sto_trc
+            jpsto2d   = jpsto2d + 3 * nn_trc_ord
+            jsto_trci(jdof) = jpsto2d - 2 * nn_trc_ord
+            jsto_trcj(jdof) = jpsto2d - 1 * nn_trc_ord
+            jsto_trck(jdof) = jpsto2d
+         END DO
+      ELSE
+         nn_sto_trc = 0
+      ENDIF
+
+      ! Set number of 3D stochastic arrays
+      jpsto3d = 0
+      IF( ln_sto_hpg ) THEN
+         IF(lwp) WRITE(numout,*) '       - stochastic horizontal pressure gradient'
+         jpsto3d   = jpsto3d + 2
+         jsto_hpgi = jpsto3d - 1
+         jsto_hpgj = jpsto3d
+      ENDIF
+      IF( ln_sto_trd ) THEN
+         IF(lwp) WRITE(numout,*) '       - stochastic trend'
+         jpsto3d   = jpsto3d + 1
+         jsto_trd  = jpsto3d
+      ENDIF
+
+      ! Allocate 2D stochastic arrays
+      IF ( jpsto2d > 0 ) THEN
+         ALLOCATE ( sto2d(jpi,jpj,jpsto2d) )
+         ALLOCATE ( sto2d_abc(jpsto2d,3) )
+         ALLOCATE ( sto2d_ave(jpsto2d) )
+         ALLOCATE ( sto2d_std(jpsto2d) )
+         ALLOCATE ( sto2d_lim(jpsto2d) )
+         ALLOCATE ( sto2d_tcor(jpsto2d) )
+         ALLOCATE ( sto2d_ord(jpsto2d) )
+         ALLOCATE ( sto2d_typ(jpsto2d) )
+         ALLOCATE ( sto2d_sgn(jpsto2d) )
+         ALLOCATE ( sto2d_flt(jpsto2d) )
+         ALLOCATE ( sto2d_fac(jpsto2d) )
+      ENDIF
+
+      ! Allocate 3D stochastic arrays
+      IF ( jpsto3d > 0 ) THEN
+         ALLOCATE ( sto3d(jpi,jpj,jpk,jpsto3d) )
+         ALLOCATE ( sto3d_abc(jpsto3d,3) )
+         ALLOCATE ( sto3d_ave(jpsto3d) )
+         ALLOCATE ( sto3d_std(jpsto3d) )
+         ALLOCATE ( sto3d_lim(jpsto3d) )
+         ALLOCATE ( sto3d_tcor(jpsto3d) )
+         ALLOCATE ( sto3d_ord(jpsto3d) )
+         ALLOCATE ( sto3d_typ(jpsto3d) )
+         ALLOCATE ( sto3d_sgn(jpsto3d) )
+         ALLOCATE ( sto3d_flt(jpsto3d) )
+         ALLOCATE ( sto3d_fac(jpsto3d) )
+      ENDIF
+
+      ! Allocate temporary workspace
+      IF ( jpsto2d > 0 .OR. jpsto3d > 0 ) THEN
+         ALLOCATE ( sto_tmp(jpi,jpj) ) ; sto_tmp(:,:) = 0._wp
+      ENDIF
+
+      ! 1) For every stochastic parameter:
+      ! ----------------------------------
+      ! - set nature of grid point and control of the sign
+      !       across the north fold (sto2d_typ, sto2d_sgn)
+      ! - set number of passes of Laplacian filter (sto2d_flt)
+      ! - set order of every autoregressive process (sto2d_ord)
+      DO jsto = 1, jpsto2d
+         sto2d_typ(jsto) = 'T'
+         sto2d_sgn(jsto) = 1._wp
+         sto2d_flt(jsto) = 0
+         sto2d_ord(jsto) = 1
+         DO jord = 0, nn_pstar_ord-1
+            IF ( jsto+jord == jsto_pstar ) THEN ! Stochastic ice strength (ave=1)
+               sto2d_ord(jsto) = nn_pstar_ord - jord
+               sto2d_flt(jsto) = nn_pstar_flt
+            ENDIF
+         ENDDO
+         DO jdof = 1, nn_sto_eos
+         DO jord = 0, nn_eos_ord-1
+            IF ( jsto+jord == jsto_eosi(jdof) ) THEN ! Stochastic equation of state i (ave=0)
+               sto2d_ord(jsto) = nn_eos_ord - jord
+               sto2d_sgn(jsto) = -1._wp
+               sto2d_flt(jsto) = nn_eos_flt
+            ENDIF
+            IF ( jsto+jord == jsto_eosj(jdof) ) THEN ! Stochastic equation of state j (ave=0)
+               sto2d_ord(jsto) = nn_eos_ord - jord
+               sto2d_sgn(jsto) = -1._wp
+               sto2d_flt(jsto) = nn_eos_flt
+            ENDIF
+            IF ( jsto+jord == jsto_eosk(jdof) ) THEN ! Stochastic equation of state k (ave=0)
+               sto2d_ord(jsto) = nn_eos_ord - jord
+               sto2d_flt(jsto) = nn_eos_flt
+            ENDIF
+         END DO
+         END DO
+         DO jdof = 1, nn_sto_trc
+         DO jord = 0, nn_trc_ord-1
+            IF ( jsto+jord == jsto_trci(jdof) ) THEN ! Stochastic tracers dynamics i (ave=0)
+               sto2d_ord(jsto) = nn_trc_ord - jord
+               sto2d_sgn(jsto) = -1._wp
+               sto2d_flt(jsto) = nn_trc_flt
+            ENDIF
+            IF ( jsto+jord == jsto_trcj(jdof) ) THEN ! Stochastic tracers dynamics j (ave=0)
+               sto2d_ord(jsto) = nn_trc_ord - jord
+               sto2d_sgn(jsto) = -1._wp
+               sto2d_flt(jsto) = nn_trc_flt
+            ENDIF
+            IF ( jsto+jord == jsto_trck(jdof) ) THEN ! Stochastic tracers dynamics k (ave=0)
+               sto2d_ord(jsto) = nn_trc_ord - jord
+               sto2d_flt(jsto) = nn_trc_flt
+            ENDIF
+         END DO
+         END DO
+
+         sto2d_fac(jsto) = sto_par_flt_fac ( sto2d_flt(jsto) )
+      END DO
+      !
+      DO jsto = 1, jpsto3d
+         sto3d_typ(jsto) = 'T'
+         sto3d_sgn(jsto) = 1._wp
+         sto3d_flt(jsto) = 0
+         sto3d_ord(jsto) = 1
+         IF ( jsto == jsto_hpgi ) THEN ! Stochastic density gradient i (ave=1)
+            sto3d_typ(jsto) = 'U'
+         ENDIF
+         IF ( jsto == jsto_hpgj ) THEN ! Stochastic density gradient j (ave=1)
+            sto3d_typ(jsto) = 'V'
+         ENDIF
+         sto3d_fac(jsto) = sto_par_flt_fac ( sto3d_flt(jsto) )
+      END DO
+
+      ! 2) For every stochastic parameter:
+      ! ----------------------------------
+      ! set average, standard deviation and time correlation
+      DO jsto = 1, jpsto2d
+         sto2d_ave(jsto)  = 0._wp
+         sto2d_std(jsto)  = 1._wp
+         sto2d_tcor(jsto) = 1._wp
+         sto2d_lim(jsto)  = 3._wp
+         IF ( jsto == jsto_ldf  ) THEN ! Stochastic lateral diffusion (ave=1)
+            sto2d_ave(jsto)  = 1._wp
+            sto2d_std(jsto)  = rn_ldf_std
+            sto2d_tcor(jsto) = rn_ldf_tcor
+         ENDIF
+         DO jord = 0, nn_pstar_ord-1
+            IF ( jsto+jord == jsto_pstar ) THEN ! Stochastic ice strength (ave=1)
+               sto2d_std(jsto) = 1._wp
+               sto2d_tcor(jsto) = rn_pstar_tcor
+            ENDIF
+         ENDDO
+         DO jdof = 1, nn_sto_eos
+         DO jord = 0, nn_eos_ord-1
+            IF ( jsto+jord == jsto_eosi(jdof) ) THEN ! Stochastic equation of state i (ave=0)
+               sto2d_std(jsto)  = rn_eos_stdxy
+               sto2d_tcor(jsto) = rn_eos_tcor
+               sto2d_lim(jsto)  = rn_eos_lim
+            ENDIF
+            IF ( jsto+jord == jsto_eosj(jdof) ) THEN ! Stochastic equation of state j (ave=0)
+               sto2d_std(jsto)  = rn_eos_stdxy
+               sto2d_tcor(jsto) = rn_eos_tcor
+               sto2d_lim(jsto)  = rn_eos_lim
+            ENDIF
+            IF ( jsto+jord == jsto_eosk(jdof) ) THEN ! Stochastic equation of state k (ave=0)
+               sto2d_std(jsto)  = rn_eos_stdz
+               sto2d_tcor(jsto) = rn_eos_tcor
+               sto2d_lim(jsto)  = rn_eos_lim
+            ENDIF
+         END DO
+         END DO
+         DO jdof = 1, nn_sto_trc
+         DO jord = 0, nn_trc_ord-1
+            IF ( jsto+jord == jsto_trci(jdof) ) THEN ! Stochastic tracer dynamics i (ave=0)
+               sto2d_std(jsto)  = rn_trc_stdxy
+               sto2d_tcor(jsto) = rn_trc_tcor
+               sto2d_lim(jsto)  = rn_trc_lim
+            ENDIF
+            IF ( jsto+jord == jsto_trcj(jdof) ) THEN ! Stochastic tracer dynamics j (ave=0)
+               sto2d_std(jsto)  = rn_trc_stdxy
+               sto2d_tcor(jsto) = rn_trc_tcor
+               sto2d_lim(jsto)  = rn_trc_lim
+            ENDIF
+            IF ( jsto+jord == jsto_trck(jdof) ) THEN ! Stochastic tracer dynamics k (ave=0)
+               sto2d_std(jsto)  = rn_trc_stdz
+               sto2d_tcor(jsto) = rn_trc_tcor
+               sto2d_lim(jsto)  = rn_trc_lim
+            ENDIF
+         END DO
+         END DO
+
+      END DO
+      !
+      DO jsto = 1, jpsto3d
+         sto3d_ave(jsto)  = 0._wp
+         sto3d_std(jsto)  = 1._wp
+         sto3d_tcor(jsto) = 1._wp
+         sto3d_lim(jsto)  = 3._wp
+         IF ( jsto == jsto_hpgi ) THEN ! Stochastic density gradient i (ave=1)
+            sto3d_ave(jsto)  = 1._wp
+            sto3d_std(jsto)  = rn_hpg_std
+            sto3d_tcor(jsto) = rn_hpg_tcor
+         ENDIF
+         IF ( jsto == jsto_hpgj ) THEN ! Stochastic density gradient j (ave=1)
+            sto3d_ave(jsto)  = 1._wp
+            sto3d_std(jsto)  = rn_hpg_std
+            sto3d_tcor(jsto) = rn_hpg_tcor
+         ENDIF
+         IF ( jsto == jsto_trd ) THEN ! Stochastic trend (ave=1)
+            sto3d_ave(jsto)  = 1._wp
+            sto3d_std(jsto)  = rn_trd_std
+            sto3d_tcor(jsto) = rn_trd_tcor
+         ENDIF
+      END DO
+
+      ! 3) For every stochastic parameter:
+      ! ----------------------------------
+      ! - compute parameters (a, b, c) of the AR1 autoregressive process
+      !   from expected value (ave), standard deviation (std)
+      !   and time correlation (tcor):
+      !     a = EXP ( - 1 / tcor )           --> sto2d_abc(:,1)
+      !     b = std * SQRT( 1 - a * a )      --> sto2d_abc(:,2)
+      !     c = ave * ( 1 - a )              --> sto2d_abc(:,3)
+      ! - for higher order processes (ARn, n>1), use approximate formula
+      !   for the b parameter (valid for tcor>>1 time step)
+      DO jsto = 1, jpsto2d
+         IF ( sto2d_tcor(jsto) == 0._wp ) THEN
+            sto2d_abc(jsto,1) = 0._wp
+         ELSE
+            sto2d_abc(jsto,1) = EXP ( - 1._wp / sto2d_tcor(jsto) )
+         ENDIF
+         IF ( sto2d_ord(jsto) == 1 ) THEN      ! Exact formula for 1st order process
+            rinflate = sto2d_std(jsto)
+         ELSE
+            ! Approximate formula, valid for tcor >> 1
+            jordm1 = sto2d_ord(jsto) - 1
+            rinflate = SQRT ( REAL( jordm1 , wp ) / REAL( 2*(2*jordm1-1) , wp ) )
+         ENDIF
+         sto2d_abc(jsto,2) = rinflate * SQRT ( 1._wp - sto2d_abc(jsto,1) &
+                                                     * sto2d_abc(jsto,1) )
+         sto2d_abc(jsto,3) = sto2d_ave(jsto) * ( 1._wp - sto2d_abc(jsto,1) )
+      END DO
+      !
+      DO jsto = 1, jpsto3d
+         IF ( sto3d_tcor(jsto) == 0._wp ) THEN
+            sto3d_abc(jsto,1) = 0._wp
+         ELSE
+            sto3d_abc(jsto,1) = EXP ( - 1._wp / sto3d_tcor(jsto) )
+         ENDIF
+         IF ( sto3d_ord(jsto) == 1 ) THEN      ! Exact formula for 1st order process
+            rinflate = sto3d_std(jsto)
+         ELSE
+            ! Approximate formula, valid for tcor >> 1
+            jordm1 = sto3d_ord(jsto) - 1
+            rinflate = SQRT ( REAL( jordm1 , wp ) / REAL( 2*(2*jordm1-1) , wp ) )
+         ENDIF
+         sto3d_abc(jsto,2) = rinflate * SQRT ( 1._wp - sto3d_abc(jsto,1) &
+                                                     * sto3d_abc(jsto,1) )
+         sto3d_abc(jsto,3) = sto3d_ave(jsto) * ( 1._wp - sto3d_abc(jsto,1) )
+      END DO
+
+      ! 4) Initialize seeds for random number generator
+      ! -----------------------------------------------
+      ! using different seeds for different processors (jarea)
+      ! and different ensemble members (jmem)
+      CALL kiss_reset( )
+      DO jarea = 1, narea
+         !DO jmem = 0, nmember
+         zseed1 = kiss() ; zseed2 = kiss() ; zseed3 = kiss() ; zseed4 = kiss()
+         !END DO
+      END DO
+      CALL kiss_seed( zseed1, zseed2, zseed3, zseed4 )
+
+      ! 5) Initialize stochastic parameters to: ave + std * w
+      ! -----------------------------------------------------
+      DO jsto = 1, jpsto2d
+         ! Draw random numbers from N(0,1) --> w
+         CALL sto_par_white( sto2d(:,:,jsto) )
+         ! Apply horizontal Laplacian filter to w
+         DO jflt = 1, sto2d_flt(jsto)
+            CALL lbc_lnk( sto2d(:,:,jsto), sto2d_typ(jsto), sto2d_sgn(jsto) )
+            CALL sto_par_flt( sto2d(:,:,jsto) )
+         END DO
+         ! Factor to restore standard deviation after filtering
+         sto2d(:,:,jsto) = sto2d(:,:,jsto) * sto2d_fac(jsto)
+         ! Limit random parameter to the limitation factor
+         sto2d(:,:,jsto) = SIGN(MIN(sto2d_lim(jsto),ABS(sto2d(:,:,jsto))),sto2d(:,:,jsto))
+         ! Multiply by standard devation and add average value
+         sto2d(:,:,jsto) = sto2d(:,:,jsto) * sto2d_std(jsto) + sto2d_ave(jsto)
+      END DO
+      !
+      DO jsto = 1, jpsto3d
+         DO jk = 1, jpk
+            ! Draw random numbers from N(0,1) --> w
+            CALL sto_par_white( sto3d(:,:,jk,jsto) )
+            ! Apply horizontal Laplacian filter to w
+            DO jflt = 1, sto3d_flt(jsto)
+               CALL lbc_lnk( sto3d(:,:,jk,jsto), sto3d_typ(jsto), sto3d_sgn(jsto) )
+               CALL sto_par_flt( sto3d(:,:,jk,jsto) )
+            END DO
+            ! Factor to restore standard deviation after filtering
+            sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) * sto3d_fac(jsto)
+            ! Limit random parameter to the limitation factor
+            sto3d(:,:,jk,jsto) = SIGN(MIN(sto3d_lim(jsto),ABS(sto3d(:,:,jk,jsto))),sto3d(:,:,jk,jsto))
+            ! Multiply by standard devation and add average value
+            sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) * sto3d_std(jsto) + sto3d_ave(jsto)
+         END DO
+      END DO
+
+      ! 6) Restart stochastic parameters from file
+      ! ------------------------------------------
+      IF( ln_rststo ) CALL sto_rst_read
+
+      ! Allocate drho_ran
+      ALLOCATE(drho_ran(jpi,jpj,jpk))
+
+   END SUBROUTINE sto_par_init
+
+
+   SUBROUTINE sto_rst_read
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sto_rst_read  ***
+      !!
+
+      !! ** Purpose :   read stochastic parameters from restart file
+      !!----------------------------------------------------------------------
+
+      INTEGER  :: jsto, jseed
+      INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type
+      REAL(KIND=8)        ::   zrseed(4)           ! RNG seeds in real type (with same bits to save in restart)
+      CHARACTER(LEN=9)    ::   clsto2d='sto2d_000' ! stochastic parameter variable name
+      CHARACTER(LEN=9)    ::   clsto3d='sto3d_000' ! stochastic parameter variable name
+      CHARACTER(LEN=10)   ::   clseed='seed0_0000' ! seed variable name
+
+      IF ( jpsto2d > 0 .OR. jpsto3d > 0 ) THEN
+
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'sto_rst_read : read stochastic parameters from restart file'
+            WRITE(numout,*) '~~~~~~~~~~~~'
+         ENDIF
+
+         ! Open the restart file
+         CALL iom_open( cn_storst_in, numstor, kiolib = jprstlib )
+
+         ! Get stochastic parameters from restart file:
+         ! 2D stochastic parameters
+         DO jsto = 1 , jpsto2d
+            WRITE(clsto2d(7:9),'(i3.3)') jsto
+            CALL iom_get( numstor, jpdom_autoglo, clsto2d , sto2d(:,:,jsto) )
+         END DO
+         ! 3D stochastic parameters
+         DO jsto = 1 , jpsto3d
+            WRITE(clsto3d(7:9),'(i3.3)') jsto
+            CALL iom_get( numstor, jpdom_autoglo, clsto3d , sto3d(:,:,:,jsto) )
+         END DO
+
+         IF (ln_rstseed) THEN
+            ! Get saved state of the random number generator
+            DO jseed = 1 , 4
+               WRITE(clseed(5:5) ,'(i1.1)') jseed
+               WRITE(clseed(7:10),'(i4.4)') narea
+               CALL iom_get( numstor, clseed , zrseed(jseed) )
+            END DO
+            ziseed = TRANSFER( zrseed , ziseed)
+            CALL kiss_seed( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
+         ENDIF
+
+         ! Close the restart file
+         CALL iom_close( numstor )
+
+      ENDIF
+
+   END SUBROUTINE sto_rst_read
+
+
+   SUBROUTINE sto_rst_write( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sto_rst_write  ***
+      !!
+      !! ** Purpose :   write stochastic parameters in restart file
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt     ! ocean time-step
+      !!
+      INTEGER  :: jsto, jseed
+      INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type
+      REAL(KIND=8)        ::   zrseed(4)           ! RNG seeds in real type (with same bits to save in restart)
+      CHARACTER(LEN=20)   ::   clkt                ! ocean time-step defined as a character
+      CHARACTER(LEN=50)   ::   clname              ! restart file name
+      CHARACTER(LEN=9)    ::   clsto2d='sto2d_000' ! stochastic parameter variable name
+      CHARACTER(LEN=9)    ::   clsto3d='sto3d_000' ! stochastic parameter variable name
+      CHARACTER(LEN=10)   ::   clseed='seed0_0000' ! seed variable name
+
+      IF ( jpsto2d > 0 .OR. jpsto3d > 0 ) THEN
+
+         IF( kt == nitrst .OR. kt == nitend ) THEN
+            IF(lwp) THEN
+               WRITE(numout,*)
+               WRITE(numout,*) 'sto_rst_write : write stochastic parameters in restart file'
+               WRITE(numout,*) '~~~~~~~~~~~~~'
+            ENDIF
+         ENDIF
+
+         ! Put stochastic parameters in restart files
+         ! (as opened at previous timestep, see below)
+         IF( kt > nit000) THEN
+         IF( kt == nitrst .OR. kt == nitend ) THEN
+            ! get and save current state of the random number generator
+            CALL kiss_state( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
+            zrseed = TRANSFER( ziseed , zrseed)
+            DO jseed = 1 , 4
+               WRITE(clseed(5:5) ,'(i1.1)') jseed
+               WRITE(clseed(7:10),'(i4.4)') narea
+               CALL iom_rstput( kt, nitrst, numstow, clseed , zrseed(jseed) )
+            END DO
+            ! 2D stochastic parameters
+            DO jsto = 1 , jpsto2d
+               WRITE(clsto2d(7:9),'(i3.3)') jsto
+               CALL iom_rstput( kt, nitrst, numstow, clsto2d , sto2d(:,:,jsto) )
+            END DO
+            ! 3D stochastic parameters
+            DO jsto = 1 , jpsto3d
+               WRITE(clsto3d(7:9),'(i3.3)') jsto
+               CALL iom_rstput( kt, nitrst, numstow, clsto3d , sto3d(:,:,:,jsto) )
+            END DO
+            ! Save drho_ran in restart file
+            CALL iom_rstput( kt, nitrst, numstow, 'drho' , drho_ran(:,:,:) )
+            ! close the restart file
+            CALL iom_close( numstow )
+         ENDIF
+         ENDIF
+
+         ! Open the restart file one timestep before writing restart
+         IF( kt < nitend) THEN
+         IF( kt == nitrst - 1 .OR. nstock == 1 .OR. kt == nitend-1 ) THEN
+            ! create the filename
+            IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst
+            ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst
+            ENDIF
+            clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_storst_out)
+            ! print information
+            IF(lwp) THEN
+               WRITE(numout,*) '             open stochastic parameters restart file: '//clname
+               IF( kt == nitrst - 1 ) THEN
+                  WRITE(numout,*) '             kt = nitrst - 1 = ', kt
+               ELSE
+                  WRITE(numout,*) '             kt = '             , kt
+               ENDIF
+            ENDIF
+            ! open the restart file
+            CALL iom_open( clname, numstow, ldwrt = .TRUE., kiolib = jprstlib )
+         ENDIF
+         ENDIF
+
+      ENDIF
+
+   END SUBROUTINE sto_rst_write
+
+
+   SUBROUTINE sto_par_white( psto )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sto_par_white  ***
+      !!
+      !! ** Purpose :   fill input array with white Gaussian noise
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out)           ::   psto
+      !!
+      INTEGER  :: ji, jj
+      REAL(KIND=8) :: gran   ! Gaussian random number (forced KIND=8 as in kiss_gaussian)
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            CALL kiss_gaussian( gran )
+            psto(ji,jj) = gran
+         END DO
+      END DO
+
+   END SUBROUTINE sto_par_white
+
+
+   SUBROUTINE sto_par_flt( psto )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sto_par_flt  ***
+      !!
+      !! ** Purpose :   apply horizontal Laplacian filter to input array
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out)           ::   psto
+      !!
+      INTEGER  :: ji, jj
+
+      DO jj = 2, jpj-1
+         DO ji = 2, jpi-1
+            psto(ji,jj) = 0.5_wp * psto(ji,jj) + 0.125_wp * &
+                              &  ( psto(ji-1,jj) + psto(ji+1,jj) +  &
+                              &    psto(ji,jj-1) + psto(ji,jj+1) )
+         END DO
+      END DO
+
+   END SUBROUTINE sto_par_flt
+
+
+   REAL(wp) FUNCTION sto_par_flt_fac( kpasses )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION sto_par_flt_fac  ***
+      !!
+      !! ** Purpose :   compute factor to restore standard deviation
+      !!                as a function of the number of passes
+      !!                of the Laplacian filter
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kpasses
+      !!
+      INTEGER :: jpasses, ji, jj, jflti, jfltj
+      INTEGER, DIMENSION(-1:1,-1:1) :: pflt0
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pfltb
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pflta
+      REAL(wp) :: ratio
+
+      pflt0(-1,-1) = 0 ; pflt0(-1,0) = 1 ; pflt0(-1,1) = 0
+      pflt0( 0,-1) = 1 ; pflt0( 0,0) = 4 ; pflt0( 0,1) = 1
+      pflt0( 1,-1) = 0 ; pflt0( 1,0) = 1 ; pflt0( 1,1) = 0
+
+      ALLOCATE(pfltb(-kpasses-1:kpasses+1,-kpasses-1:kpasses+1))
+      ALLOCATE(pflta(-kpasses-1:kpasses+1,-kpasses-1:kpasses+1))
+
+      pfltb(:,:) = 0
+      pfltb(0,0) = 1
+      DO jpasses = 1, kpasses
+        pflta(:,:) = 0
+        DO jflti= -1, 1
+        DO jfltj= -1, 1
+          DO ji= -kpasses, kpasses
+          DO jj= -kpasses, kpasses
+            pflta(ji,jj) = pflta(ji,jj) + pfltb(ji+jflti,jj+jfltj) * pflt0(jflti,jfltj)
+          ENDDO
+          ENDDO
+        ENDDO
+        ENDDO
+        pfltb(:,:) = pflta(:,:)
+      ENDDO
+
+      ratio = SUM(pfltb(:,:))
+      ratio = ratio * ratio / SUM(pfltb(:,:)*pfltb(:,:))
+      ratio = SQRT(ratio)
+
+      DEALLOCATE(pfltb,pflta)
+
+      sto_par_flt_fac = ratio
+
+   END FUNCTION sto_par_flt_fac
+
+
+END MODULE stopar
+
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/stopts.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/stopts.F90	(revision 5382)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/stopts.F90	(revision 5382)
@@ -0,0 +1,146 @@
+MODULE stopts
+   !!==============================================================================
+   !!                       ***  MODULE  stopts  ***
+   !! Stochastic parameterization: compute stochastic tracer fluctuations
+   !!==============================================================================
+   !! History :  3.3  ! 2011-12 (J.-M. Brankart)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sto_pts        : compute current stochastic tracer fluctuations
+   !!   sto_pts_init   : initialisation for stochastic tracer fluctuations
+   !!----------------------------------------------------------------------
+   USE dom_oce         ! ocean space and time domain
+   USE lbclnk          ! lateral boundary conditions (or mpp link)
+   USE phycst          ! physical constants
+   USE stopar          ! stochastic parameterization
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sto_pts         ! called by step.F90
+   PUBLIC   sto_pts_init    ! called by nemogcm.F90
+
+   ! Public array with random tracer fluctuations
+   REAL(wp), PUBLIC, DIMENSION(:,:,:,:,:), ALLOCATABLE :: pts_ran
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: eosbn2.F90 2528 2010-12-27 17:33:53Z rblod $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE sto_pts( pts )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sto_pts  ***
+      !!
+      !! ** Purpose :   Compute current stochastic tracer fluctuations
+      !!
+      !! ** Method  :   Compute tracer differences from a random walk
+      !!                around every model grid point
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout) ::   pts   ! 1 : potential temperature  [Celsius]
+      !                                                               ! 2 : salinity               [psu]
+      INTEGER  ::   ji, jj, jk, jts, jdof ! dummy loop indices
+      INTEGER  ::   jim1, jjm1, jkm1  ! incremented indices
+      INTEGER  ::   jip1, jjp1, jkp1  !     -          -
+      REAL(wp) ::   zdtsim, zdtsjm, zdtskm         ! temporary scalars
+      REAL(wp) ::   zdtsip, zdtsjp, zdtskp, zdts   !     -        -
+      !!----------------------------------------------------------------------
+
+      DO jts = 1, jpts
+        CALL lbc_lnk( pts(:,:,:,jts), 'T' , 1._wp )
+      ENDDO
+
+      DO jdof = 1, nn_sto_eos
+        DO jts = 1, jpts
+           DO jk = 1, jpkm1
+              jkm1 = MAX(jk-1,1) ; jkp1 = MIN(jk+1,jpkm1)
+              DO jj = 1, jpj
+                 jjm1 = MAX(jj-1,1) ; jjp1 = MIN(jj+1,jpj)
+                 DO ji = 1, jpi
+                    jim1 = MAX(ji-1,1) ; jip1 = MIN(ji+1,jpi)
+                    !
+                    ! compute tracer gradient
+                    zdtsip = ( pts(jip1,jj,jk,jts) - pts(ji,jj,jk,jts) ) * tmask(jip1,jj,jk)
+                    zdtsim = ( pts(ji,jj,jk,jts) - pts(jim1,jj,jk,jts) ) * tmask(jim1,jj,jk)
+                    zdtsjp = ( pts(ji,jjp1,jk,jts) - pts(ji,jj,jk,jts) ) * tmask(ji,jjp1,jk)
+                    zdtsjm = ( pts(ji,jj,jk,jts) - pts(ji,jjm1,jk,jts) ) * tmask(ji,jjm1,jk)
+                    zdtskp = ( pts(ji,jj,jkp1,jts) - pts(ji,jj,jk,jts) ) * tmask(ji,jj,jkp1)
+                    zdtskm = ( pts(ji,jj,jk,jts) - pts(ji,jj,jkm1,jts) ) * tmask(ji,jj,jkm1)
+                    !
+                    ! compute random tracer fluctuation (zdts)
+                    zdts   = ( zdtsip + zdtsim ) * sto2d(ji,jj,jsto_eosi(jdof)) + &
+                           & ( zdtsjp + zdtsjm ) * sto2d(ji,jj,jsto_eosj(jdof)) + &
+                           & ( zdtskp + zdtskm ) * sto2d(ji,jj,jsto_eosk(jdof))
+!                   zdts   = zdtsip * MAX(sto2d(ji,jj,jsto_eosi),0._wp) + &
+!                          & zdtsim * MIN(sto2d(ji,jj,jsto_eosi),0._wp) + &
+!                          & zdtsjp * MAX(sto2d(ji,jj,jsto_eosj),0._wp) + &
+!                          & zdtsjm * MIN(sto2d(ji,jj,jsto_eosj),0._wp) + &
+!                          & zdtskp * MAX(sto2d(ji,jj,jsto_eosk),0._wp) + &
+!                          & zdtskm * MIN(sto2d(ji,jj,jsto_eosk),0._wp)
+                    zdts   = zdts  * tmask(ji,jj,jk) *SIN( gphit(ji,jj) * rad )
+                    pts_ran(ji,jj,jk,jts,jdof) = zdts * 0.5_wp
+                    !
+                  END DO
+               END DO
+            END DO
+         END DO
+      END DO
+
+      ! Eliminate any possible negative salinity
+      DO jdof = 1, nn_sto_eos
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  pts_ran(ji,jj,jk,jp_sal,jdof) = MIN( ABS(pts_ran(ji,jj,jk,jp_sal,jdof)) ,  &
+                                                &      MAX(pts(ji,jj,jk,jp_sal),0._wp) )     &
+                                                &  * SIGN(1._wp,pts_ran(ji,jj,jk,jp_sal,jdof))
+               END DO
+            END DO
+         END DO
+      END DO
+
+      ! Eliminate any temperature lower than -2 degC
+!     DO jdof = 1, nn_sto_eos
+!        DO jk = 1, jpkm1
+!           DO jj = 1, jpj
+!              DO ji = 1, jpi
+!                 pts_ran(ji,jj,jk,jp_tem,jdof) = MIN( ABS(pts_ran(ji,jj,jk,jp_tem,jdof)) ,    &
+!                                               &      MAX(pts(ji,jj,jk,jp_tem)+2._wp,0._wp) ) &
+!                                               &  * SIGN(1._wp,pts_ran(ji,jj,jk,jp_tem,jdof))
+!              END DO
+!           END DO
+!        END DO
+!     END DO
+
+
+      ! Lateral boundary conditions on pts_ran
+      DO jdof = 1, nn_sto_eos
+         DO jts = 1, jpts
+            CALL lbc_lnk( pts_ran(:,:,:,jts,jdof), 'T' , 1._wp )
+         END DO
+      END DO
+
+   END SUBROUTINE sto_pts
+
+
+   SUBROUTINE sto_pts_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sto_pts_init  ***
+      !!
+      !! ** Purpose :   Initialisation for stochastic tracer fluctuations
+      !!
+      !! ** Method  :   Allocate required array
+      !!
+      !!----------------------------------------------------------------------
+
+      ALLOCATE(pts_ran(jpi,jpj,jpk,jpts,nn_sto_eos))
+
+   END SUBROUTINE sto_pts_init
+
+END MODULE stopts
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/storng.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/storng.F90	(revision 5382)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/storng.F90	(revision 5382)
@@ -0,0 +1,400 @@
+MODULE storng
+!$AGRIF_DO_NOT_TREAT
+   !!======================================================================
+   !!                       ***  MODULE  storng  ***
+   !! Random number generator, used in NEMO stochastic parameterization
+   !!
+   !!=====================================================================
+   !! History :  3.3  ! 2011-10 (J.-M. Brankart)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !! The module is based on (and includes) the
+   !! 64-bit KISS (Keep It Simple Stupid) random number generator
+   !! distributed by George Marsaglia :
+   !! http://groups.google.com/group/comp.lang.fortran/
+   !!        browse_thread/thread/a85bf5f2a97f5a55
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   kiss          : 64-bit KISS random number generator (period ~ 2^250)
+   !!   kiss_seed     : Define seeds for KISS random number generator
+   !!   kiss_state    : Get current state of KISS random number generator
+   !!   kiss_save     : Save current state of KISS (for future restart)
+   !!   kiss_load     : Load the saved state of KISS
+   !!   kiss_reset    : Reset the default seeds
+   !!   kiss_check    : Check the KISS pseudo-random sequence
+   !!   kiss_uniform  : Real random numbers with uniform distribution in [0,1]
+   !!   kiss_gaussian : Real random numbers with Gaussian distribution N(0,1)
+   !!   kiss_gamma    : Real random numbers with Gamma distribution Gamma(k,1)
+   !!   kiss_sample   : Select a random sample from a set of integers
+   !!
+   !!   ---CURRENTLY NOT USED IN NEMO :
+   !!   kiss_save, kiss_load, kiss_check, kiss_gamma, kiss_sample
+   !!----------------------------------------------------------------------
+   IMPLICIT NONE
+   PRIVATE
+
+   ! Public functions/subroutines
+   PUBLIC :: kiss, kiss_seed, kiss_state, kiss_save, kiss_load, kiss_reset, kiss_check
+   PUBLIC :: kiss_uniform, kiss_gaussian, kiss_gamma, kiss_sample
+
+   ! Default/initial seeds
+   INTEGER(KIND=8) :: x=1234567890987654321_8
+   INTEGER(KIND=8) :: y=362436362436362436_8
+   INTEGER(KIND=8) :: z=1066149217761810_8
+   INTEGER(KIND=8) :: w=123456123456123456_8
+
+   ! Parameters to generate real random variates
+   REAL(KIND=8), PARAMETER :: huge64=9223372036854775808.0_8  ! +1
+   REAL(KIND=8), PARAMETER :: zero=0.0_8, half=0.5_8, one=1.0_8, two=2.0_8
+
+   ! Variables to store 2 Gaussian random numbers with current index (ig)
+   INTEGER(KIND=8), SAVE :: ig=1
+   REAL(KIND=8), SAVE :: gran1, gran2
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: dynhpg.F90 2528 2010-12-27 17:33:53Z rblod $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   FUNCTION kiss()
+      !! --------------------------------------------------------------------
+      !!                  ***  FUNCTION kiss  ***
+      !!
+      !! ** Purpose :   64-bit KISS random number generator
+      !!
+      !! ** Method  :   combine several random number generators:
+      !!                (1) Xorshift (XSH), period 2^64-1,
+      !!                (2) Multiply-with-carry (MWC), period (2^121+2^63-1)
+      !!                (3) Congruential generator (CNG), period 2^64.
+      !!
+      !!                overall period:
+      !!                (2^250+2^192+2^64-2^186-2^129)/6
+      !!                            ~= 2^(247.42) or 10^(74.48)
+      !!
+      !!                set your own seeds with 'kiss_seed'
+      ! --------------------------------------------------------------------
+      IMPLICIT NONE
+      INTEGER(KIND=8) :: kiss, t
+
+      t = ISHFT(x,58) + w
+      IF (s(x).eq.s(t)) THEN
+         w = ISHFT(x,-6) + s(x)
+      ELSE
+         w = ISHFT(x,-6) + 1 - s(x+t)
+      ENDIF
+      x = t + x
+      y = m( m( m(y,13_8), -17_8 ), 43_8 )
+      z = 6906969069_8 * z + 1234567_8
+
+      kiss = x + y + z
+
+      CONTAINS
+
+         FUNCTION s(k)
+            INTEGER(KIND=8) :: s, k
+            s = ISHFT(k,-63)
+         END FUNCTION s
+
+         FUNCTION m(k, n)
+            INTEGER(KIND=8) :: m, k, n
+            m =  IEOR(k, ISHFT(k, n) )
+         END FUNCTION m
+
+   END FUNCTION kiss
+
+
+   SUBROUTINE kiss_seed(ix, iy, iz, iw)
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_seed  ***
+      !!
+      !! ** Purpose :   Define seeds for KISS random number generator
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      INTEGER(KIND=8) :: ix, iy, iz, iw
+
+      x = ix
+      y = iy
+      z = iz
+      w = iw
+
+   END SUBROUTINE kiss_seed
+
+
+   SUBROUTINE kiss_state(ix, iy, iz, iw)
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_state  ***
+      !!
+      !! ** Purpose :   Get current state of KISS random number generator
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      INTEGER(KIND=8) :: ix, iy, iz, iw
+
+      ix = x
+      iy = y
+      iz = z
+      iw = w
+
+   END SUBROUTINE kiss_state
+
+
+   SUBROUTINE kiss_reset()
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_reset  ***
+      !!
+      !! ** Purpose :   Reset the default seeds for KISS random number generator
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+
+      x=1234567890987654321_8
+      y=362436362436362436_8
+      z=1066149217761810_8
+      w=123456123456123456_8
+
+    END SUBROUTINE kiss_reset
+
+
+    SUBROUTINE kiss_check(check_type)
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_check  ***
+      !!
+      !! ** Purpose :   Check the KISS pseudo-random sequence
+      !!
+      !! ** Method  :   Check that it reproduces the correct sequence
+      !!                from the default seed
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      INTEGER(KIND=8) :: iter, niter, correct, iran
+      CHARACTER(LEN=*) :: check_type
+      LOGICAL :: print_success
+
+      ! Save current state of KISS
+      CALL kiss_save()
+      ! Reset the default seed
+      CALL kiss_reset()
+
+      ! Select check type
+      SELECT CASE(check_type)
+      CASE('short')
+         niter = 5_8
+         correct = 542381058189297533_8
+         print_success = .FALSE.
+      CASE('long')
+         niter = 100000000_8
+         correct = 1666297717051644203_8 ! Check provided by G. Marsaglia
+         print_success = .TRUE.
+      CASE('default')
+      CASE DEFAULT
+         STOP 'Bad check type in kiss_check'
+      END SELECT
+
+      ! Run kiss for the required number of iterations (niter)
+      DO iter=1,niter
+         iran = kiss()
+      ENDDO
+
+      ! Check that last iterate is correct
+      IF (iran.NE.correct) THEN
+         STOP 'Check failed: KISS internal error !!'
+      ELSE
+         IF (print_success) PRINT *, 'Check successful: 100 million calls to KISS OK'
+      ENDIF
+
+      ! Reload the previous state of KISS
+      CALL kiss_load()
+
+   END SUBROUTINE kiss_check
+
+
+   SUBROUTINE kiss_save
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_save  ***
+      !!
+      !! ** Purpose :   Save current state of KISS random number generator
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+
+      OPEN(UNIT=30,FILE='.kiss_restart')
+      WRITE(30,*) x
+      WRITE(30,*) y
+      WRITE(30,*) z
+      WRITE(30,*) w
+      CLOSE(30)
+
+    END SUBROUTINE kiss_save
+
+
+    SUBROUTINE kiss_load
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_load  ***
+      !!
+      !! ** Purpose :   Load the saved state of KISS random number generator
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      LOGICAL :: filexists
+
+      INQUIRE(FILE='.kiss_restart',EXIST=filexists)
+      IF (filexists) THEN
+         OPEN(UNIT=30,FILE='.kiss_restart')
+         READ(30,*) x
+         READ(30,*) y
+         READ(30,*) z
+         READ(30,*) w
+         CLOSE(30)
+      ENDIF
+
+   END SUBROUTINE kiss_load
+
+
+   SUBROUTINE kiss_uniform(uran)
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_uniform  ***
+      !!
+      !! ** Purpose :   Real random numbers with uniform distribution in [0,1]
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(KIND=8) :: uran
+
+      uran = half * ( one + REAL(kiss(),8) / huge64 )
+
+   END SUBROUTINE kiss_uniform
+
+
+   SUBROUTINE kiss_gaussian(gran)
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_gaussian  ***
+      !!
+      !! ** Purpose :   Real random numbers with Gaussian distribution N(0,1)
+      !!
+      !! ** Method  :   Generate 2 new Gaussian draws (gran1 and gran2)
+      !!                from 2 uniform draws on [-1,1] (u1 and u2),
+      !!                using the Marsaglia polar method
+      !!                (see Devroye, Non-Uniform Random Variate Generation, p. 235-236)
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(KIND=8) :: gran, u1, u2, rsq, fac
+
+      IF (ig.EQ.1) THEN
+         rsq = two
+         DO WHILE ( (rsq.GE.one).OR. (rsq.EQ.zero) )
+            u1 = REAL(kiss(),8) / huge64
+            u2 = REAL(kiss(),8) / huge64
+            rsq = u1*u1 + u2*u2
+         ENDDO
+         fac = SQRT(-two*LOG(rsq)/rsq)
+         gran1 = u1 * fac
+         gran2 = u2 * fac
+      ENDIF
+
+      ! Output one of the 2 draws
+      IF (ig.EQ.1) THEN
+         gran = gran1 ; ig = 2
+      ELSE
+         gran = gran2 ; ig = 1
+      ENDIF
+
+   END SUBROUTINE kiss_gaussian
+
+
+   SUBROUTINE kiss_gamma(gamr,k)
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_gamma  ***
+      !!
+      !! ** Purpose :   Real random numbers with Gamma distribution Gamma(k,1)
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(KIND=8), PARAMETER :: p1 = 4.5_8
+      REAL(KIND=8), PARAMETER :: p2 = 2.50407739677627_8  ! 1+LOG(9/2)
+      REAL(KIND=8), PARAMETER :: p3 = 1.38629436111989_8  ! LOG(4)
+      REAL(KIND=8) :: gamr, k, u1, u2, b, c, d, xx, yy, zz, rr, ee
+      LOGICAL :: accepted
+
+      IF (k.GT.one) THEN
+         ! Cheng's rejection algorithm
+         ! (see Devroye, Non-Uniform Random Variate Generation, p. 413)
+         b = k - p3 ; d = SQRT(two*k-one) ; c = k + d
+
+         accepted=.FALSE.
+         DO WHILE (.NOT.accepted)
+            CALL kiss_uniform(u1)
+            yy = LOG(u1/(one-u1)) / d  ! Mistake in Devroye: "* k" instead of "/ d"
+            xx = k * EXP(yy)
+            rr = b + c * yy - xx
+            CALL kiss_uniform(u2)
+            zz = u1 * u1 * u2
+
+            accepted = rr .GE. (zz*p1-p2)
+            IF (.NOT.accepted) accepted =  rr .GE. LOG(zz)
+         ENDDO
+
+         gamr = xx
+
+      ELSEIF (k.LT.one) THEN
+        ! Rejection from the Weibull density
+        ! (see Devroye, Non-Uniform Random Variate Generation, p. 415)
+        c = one/k ; d = (one-k) * EXP( (k/(one-k)) * LOG(k) )
+
+        accepted=.FALSE.
+        DO WHILE (.NOT.accepted)
+           CALL kiss_uniform(u1)
+           zz = -LOG(u1)
+           xx = EXP( c * LOG(zz) )
+           CALL kiss_uniform(u2)
+           ee = -LOG(u2)
+
+           accepted = (zz+ee) .GE. (d+xx)  ! Mistake in Devroye: "LE" instead of "GE"
+        ENDDO
+
+        gamr = xx
+
+      ELSE
+         ! Exponential distribution
+         CALL kiss_uniform(u1)
+         gamr = -LOG(u1)
+
+      ENDIF
+
+   END SUBROUTINE kiss_gamma
+
+
+   SUBROUTINE kiss_sample(a,n,k)
+      !! --------------------------------------------------------------------
+      !!                  ***  ROUTINE kiss_sample  ***
+      !!
+      !! ** Purpose :   Select a random sample of size k from a set of n integers
+      !!
+      !! ** Method  :   The sample is output in the first k elements of a
+      !!                Set k equal to n to obtain a random permutation
+      !!                  of the whole set of integers
+      !!
+      !! --------------------------------------------------------------------
+      IMPLICIT NONE
+      INTEGER(KIND=8), DIMENSION(:) :: a
+      INTEGER(KIND=8) :: n, k, i, j, atmp
+      REAL(KIND=8) :: uran
+
+      ! Select the sample using the swapping method
+      ! (see Devroye, Non-Uniform Random Variate Generation, p. 612)
+      DO i=1,k
+         ! Randomly select the swapping element between i and n (inclusive)
+         CALL kiss_uniform(uran)
+         j = i - 1 + CEILING( REAL(n-i+1,8) * uran )
+         ! Swap elements i and j
+         atmp = a(i) ; a(i) = a(j) ; a(j) = atmp
+      ENDDO
+
+   END SUBROUTINE kiss_sample
+!$AGRIF_END_DO_NOT_TREAT
+END MODULE storng
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc_rst.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc_rst.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc_rst.F90	(revision 5382)
@@ -39,5 +39,6 @@
       !
       CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character
-      CHARACTER(LEN=50)   ::   clname   ! ice output restart file name
+      CHARACTER(LEN=50)   ::   clname   ! output restart file name
+      CHARACTER(LEN=256)  ::   clpath   ! full path to restart file
       CHARACTER (len=35) :: charout
       INTEGER :: jl,  jk, jn               ! loop indice
@@ -51,6 +52,8 @@
          ENDIF
          clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_trdrst_trc_out)
-         IF(lwp) WRITE(numout,*) '             open ocean restart_mld_trc NetCDF  '//clname
-         CALL iom_open( clname, nummldw_trc, ldwrt = .TRUE., kiolib = jprstlib )
+         clpath = TRIM(cn_trcrst_outdir)
+         IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
+         IF(lwp) WRITE(numout,*) '             open ocean restart_mld_trc NetCDF  'TRIM(clpath)//TRIM(clname)
+         CALL iom_open( TRIM(clpath)//TRIM(clname), nummldw_trc, ldwrt = .TRUE., kiolib = jprstlib )
       ENDIF
 
@@ -133,4 +136,5 @@
       INTEGER ::  jlibalt = jprstlib
       LOGICAL ::  llok
+      CHARACTER(LEN=256)  ::   clpath   ! full path to restart file
       !!-----------------------------------------------------------------------------
       
@@ -141,12 +145,15 @@
       ENDIF
       
+      clpath = TRIM(cn_trcrst_indir)
+      IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
+
       IF ( jprstlib == jprstdimg ) THEN
         ! eventually read netcdf file (monobloc)  for restarting on different number of processors
         ! if {cn_trdrst_trc_in}.nc exists, then set jlibalt to jpnf90
-        INQUIRE( FILE = TRIM(cn_trdrst_trc_in)//'.nc', EXIST = llok )
+        INQUIRE( FILE = TRIM(clpath)//TRIM(cn_trdrst_trc_in)//'.nc', EXIST = llok )
         IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF
       ENDIF
 
-      CALL iom_open( cn_trdrst_trc_in, inum, kiolib = jlibalt ) 
+      CALL iom_open( TRIM(clpath)//TRIM(cn_trdrst_trc_in), inum, kiolib = jlibalt ) 
       
       IF( ln_trdmxl_trc_instant ) THEN 
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/TOP_SRC/trc.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 5382)
@@ -54,5 +54,7 @@
    INTEGER             , PUBLIC                                    ::  nn_rsttr       !: control of the time step ( 0 or 1 ) for pass. tr.
    CHARACTER(len = 80) , PUBLIC                                    ::  cn_trcrst_in   !: suffix of pass. tracer restart name (input)
+   CHARACTER(len = 256), PUBLIC                                    ::  cn_trcrst_indir  !: restart input directory
    CHARACTER(len = 80) , PUBLIC                                    ::  cn_trcrst_out  !: suffix of pass. tracer restart name (output)
+   CHARACTER(len = 256), PUBLIC                                    ::  cn_trcrst_outdir  !: restart output directory
    REAL(wp)            , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::  rdttrc         !: vertical profile of passive tracer time step
    LOGICAL             , PUBLIC                                    ::  ln_top_euler  !: boolean term for euler integration 
@@ -172,5 +174,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3.1 , NEMO Consortium (2010)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/TOP_SRC/trcnam.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/TOP_SRC/trcnam.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/TOP_SRC/trcnam.F90	(revision 5382)
@@ -39,5 +39,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -175,5 +175,6 @@
       !!---------------------------------------------------------------------
       NAMELIST/namtrc_run/ nn_dttrc, nn_writetrc, ln_rsttr, nn_rsttr, ln_top_euler, &
-        &                  cn_trcrst_in, cn_trcrst_out
+        &                  cn_trcrst_indir, cn_trcrst_outdir, cn_trcrst_in, cn_trcrst_out
+
 
       INTEGER  ::   ios                 ! Local integer output status for namelist read
@@ -339,5 +340,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!======================================================================
Index: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 5381)
+++ branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 5382)
@@ -51,4 +51,5 @@
       CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
       CHARACTER(LEN=50)   ::   clname   ! trc output restart file name
+      CHARACTER(LEN=256)  ::   clpath   ! full path to ocean output restart file
       !!----------------------------------------------------------------------
       !
@@ -56,8 +57,13 @@
          IF( kt == nittrc000 ) THEN
             lrst_trc = .FALSE.
-            nitrst = nitend
-         ENDIF
-
-         IF( MOD( kt - 1, nstock ) == 0 ) THEN
+            IF( ln_rst_list ) THEN
+               nrst_lst = 1
+               nitrst = nstocklist( nrst_lst )
+            ELSE
+               nitrst = nitend
+            ENDIF
+         ENDIF
+
+         IF( .NOT. ln_rst_list .AND. MOD( kt - 1, nstock ) == 0 ) THEN
             ! we use kt - 1 and not kt - nittrc000 to keep the same periodicity from the beginning of the experiment
             nitrst = kt + nstock - 1                  ! define the next value of nitrst for restart writing
@@ -79,6 +85,9 @@
          IF(lwp) WRITE(numout,*)
          clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_trcrst_out)
-         IF(lwp) WRITE(numout,*) '             open trc restart.output NetCDF file: '//clname
-         CALL iom_open( clname, numrtw, ldwrt = .TRUE., kiolib = jprstlib )
+         clpath = TRIM(cn_trcrst_outdir)
+         IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
+         IF(lwp) WRITE(numout,*) &
+             '             open trc restart.output NetCDF file: ',TRIM(clpath)//clname
+         CALL iom_open( TRIM(clpath)//TRIM(clname), numrtw, ldwrt = .TRUE., kiolib = jprstlib )
          lrst_trc = .TRUE.
       ENDIF
@@ -140,4 +149,8 @@
           lrst_trc = .FALSE.
 #endif
+          IF( lk_offline .AND. ln_rst_list ) THEN
+             nrst_lst = nrst_lst + 1
+             nitrst = nstocklist( nrst_lst )
+          ENDIF
       ENDIF
       !
@@ -190,9 +203,9 @@
            ! eventually read netcdf file (monobloc)  for restarting on different number of processors
            ! if {cn_trcrst_in}.nc exists, then set jlibalt to jpnf90 
-           INQUIRE( FILE = TRIM(cn_trcrst_in)//'.nc', EXIST = llok )
+           INQUIRE( FILE = TRIM(cn_trcrst_indir)//'/'//TRIM(cn_trcrst_in)//'.nc', EXIST = llok )
            IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF
          ENDIF
 
-         CALL iom_open( cn_trcrst_in, numrtr, kiolib = jlibalt )
+         CALL iom_open( TRIM(cn_trcrst_indir)//'/'//cn_trcrst_in, numrtr, kiolib = jlibalt )
 
          CALL iom_get ( numrtr, 'kt', zkt )   ! last time-step of previous run
@@ -306,5 +319,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!======================================================================
