Index: /CONFIG/UNIFORM/v6/IPSLESM6/SOURCES/NEMO/domain.F90
===================================================================
--- /CONFIG/UNIFORM/v6/IPSLESM6/SOURCES/NEMO/domain.F90	(revision 3695)
+++ /CONFIG/UNIFORM/v6/IPSLESM6/SOURCES/NEMO/domain.F90	(revision 3695)
@@ -0,0 +1,480 @@
+MODULE domain
+   !!==============================================================================
+   !!                       ***  MODULE domain   ***
+   !! Ocean initialization : domain initialization
+   !!==============================================================================
+   !! History :  OPA  !  1990-10  (C. Levy - G. Madec)  Original code
+   !!                 !  1992-01  (M. Imbard) insert time step initialization
+   !!                 !  1996-06  (G. Madec) generalized vertical coordinate 
+   !!                 !  1997-02  (G. Madec) creation of domwri.F
+   !!                 !  2001-05  (E.Durand - G. Madec) insert closed sea
+   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module
+   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
+   !!            3.3  !  2010-11  (G. Madec)  initialisation in C1D configuration
+   !!            3.6  !  2013     ( J. Simeon, C. Calone, G. Madec, C. Ethe ) Online coarsening of outputs
+   !!----------------------------------------------------------------------
+   
+   !!----------------------------------------------------------------------
+   !!   dom_init       : initialize the space and time domain
+   !!   dom_nam        : read and contral domain namelists
+   !!   dom_ctl        : control print for the ocean domain
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean variables
+   USE dom_oce         ! domain: ocean
+   USE sbc_oce         ! surface boundary condition: ocean
+   USE trc_oce         ! shared ocean-passive tracers variables
+   USE phycst          ! physical constants
+   USE closea          ! closed seas
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distributed memory computing library
+
+   USE domhgr          ! domain: set the horizontal mesh
+   USE domzgr          ! domain: set the vertical mesh
+   USE domstp          ! domain: set the time-step
+   USE dommsk          ! domain: set the mask system
+   USE domwri          ! domain: write the meshmask file
+   USE domvvl          ! variable volume
+   USE c1d             ! 1D vertical configuration
+   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine)
+   USE timing          ! Timing
+   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   dom_init   ! called by opa.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!-------------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: domain.F90 7522 2017-01-02 10:06:49Z cetlod $
+   !! Software governed by the CeCILL licence        (NEMOGCM/NEMO_CeCILL.txt)
+   !!-------------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dom_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_init  ***
+      !!                    
+      !! ** Purpose :   Domain initialization. Call the routines that are 
+      !!              required to create the arrays which define the space 
+      !!              and time domain of the ocean model.
+      !!
+      !! ** Method  : - dom_msk: compute the masks from the bathymetry file
+      !!              - dom_hgr: compute or read the horizontal grid-point position
+      !!                         and scale factors, and the coriolis factor
+      !!              - dom_zgr: define the vertical coordinate and the bathymetry
+      !!              - dom_stp: defined the model time step
+      !!              - dom_wri: create the meshmask file if nmsh=1
+      !!              - 1D configuration, move Coriolis, u and v at T-point
+      !!----------------------------------------------------------------------
+      INTEGER ::   jk          ! dummy loop argument
+      INTEGER ::   iconf = 0   ! local integers
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('dom_init')
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_init : domain initialization'
+         WRITE(numout,*) '~~~~~~~~'
+      ENDIF
+      !
+                             CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
+                             CALL dom_clo      ! Closed seas and lake
+                             CALL dom_hgr      ! Horizontal mesh
+                             CALL dom_zgr      ! Vertical mesh and bathymetry
+                             CALL dom_msk      ! Masks
+      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
+      !
+      ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at T-points
+      hu_0(:,:) = 0.0_wp                       ! Reference ocean depth at U-points
+      hv_0(:,:) = 0.0_wp                       ! Reference ocean depth at V-points
+      DO jk = 1, jpk
+         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
+         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
+         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
+      END DO
+      !
+      IF( lk_c1d )           CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
+      !
+      IF( .NOT.lk_offline ) THEN
+        !
+        IF( lk_vvl )         CALL dom_vvl_init ! Vertical variable mesh
+        !
+        hu(:,:) = 0._wp                          ! Ocean depth at U-points
+        hv(:,:) = 0._wp                          ! Ocean depth at V-points
+        ht(:,:) = 0._wp                          ! Ocean depth at T-points
+        DO jk = 1, jpkm1
+           hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+           hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+           ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+        END DO
+        !                                        ! Inverse of the local depth
+        hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
+        hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
+        !
+      ENDIF
+
+                             CALL dom_stp      ! time step
+      IF( nmsh /= 0      )   CALL dom_wri      ! Create a domain file
+      IF( .NOT.ln_rstart )   CALL dom_ctl      ! Domain control
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
+      !
+   END SUBROUTINE dom_init
+
+
+   SUBROUTINE dom_nam
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dom_nam  ***
+      !!                    
+      !! ** Purpose :   read domaine namelists and print the variables.
+      !!
+      !! ** input   : - namrun namelist
+      !!              - namdom namelist
+      !!              - namcla namelist
+      !!              - namnc4 namelist   ! "key_netcdf4" only
+      !!----------------------------------------------------------------------
+      USE ioipsl
+      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
+         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , ln_rstart_ts, nn_rstctl,   &
+         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
+         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler, ln_mskutil
+      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   &
+         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  &
+         &             rn_rdtmax, rn_rdth     , nn_closea , ln_crs,    &
+         &             jphgr_msh, &
+         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
+         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
+         &             ppa2, ppkth2, ppacr2
+      NAMELIST/namcla/ nn_cla
+#if defined key_netcdf4
+      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
+#endif
+      INTEGER  ::   ios                 ! Local integer output status for namelist read
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
+      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
+      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namrun )
+      !
+      IF(lwp) THEN                  ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
+         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) '   Namelist namrun'
+         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
+         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
+         WRITE(numout,*) '      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,*) '      restart but zero for current for T ln_rstart_ts = ', ln_rstart_ts
+         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
+         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
+         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
+         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
+         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
+         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
+         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
+         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
+      ENDIF
+
+      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
+      cexper = cn_exp
+      nrstdt = nn_rstctl
+      nit000 = nn_it000
+      nitend = nn_itend
+      ndate0 = nn_date0
+      nleapy = nn_leapy
+      ninist = nn_istate
+      nstock = nn_stock
+      nstocklist = nn_stocklist
+      nwrite = nn_write
+      neuler = nn_euler
+      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
+         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
+         CALL ctl_warn( ctmp1 )
+         neuler = 0
+      ENDIF
+
+      !                             ! control of output frequency
+      IF ( nstock == 0 .OR. nstock > nitend ) THEN
+         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
+         CALL ctl_warn( ctmp1 )
+         nstock = nitend
+      ENDIF
+      IF ( nwrite == 0 ) THEN
+         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
+         CALL ctl_warn( ctmp1 )
+         nwrite = nitend
+      ENDIF
+
+#if defined key_agrif
+      IF( Agrif_Root() ) THEN
+#endif
+      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
+      CASE (  1 ) 
+         CALL ioconf_calendar('gregorian')
+         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
+      CASE (  0 )
+         CALL ioconf_calendar('noleap')
+         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
+      CASE ( 30 )
+         CALL ioconf_calendar('360d')
+         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
+      END SELECT
+#if defined key_agrif
+      ENDIF
+#endif
+
+      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
+      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
+903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
+  
+      !
+      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
+      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
+904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namdom )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '   Namelist namdom : space & time domain'
+         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
+         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
+         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
+         WRITE(numout,*) '      min number of ocean level (<0)       '
+         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
+         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
+         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
+         WRITE(numout,*) '           = 0   no file created           '
+         WRITE(numout,*) '           = 1   mesh_mask                 '
+         WRITE(numout,*) '           = 2   mesh and mask             '
+         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
+         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
+         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
+         WRITE(numout,*) '      acceleration of converge              nn_acc    = ', nn_acc
+         WRITE(numout,*) '        nn_acc=1: surface tracer rdt        rn_rdtmin = ', rn_rdtmin
+         WRITE(numout,*) '                  bottom  tracer rdt        rdtmax    = ', rn_rdtmax
+         WRITE(numout,*) '                  depth of transition       rn_rdth   = ', rn_rdth
+         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
+         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
+         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
+         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
+         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
+         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
+         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
+         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
+         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
+         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
+         WRITE(numout,*) '                                        ppa0            = ', ppa0
+         WRITE(numout,*) '                                        ppa1            = ', ppa1
+         WRITE(numout,*) '                                        ppkth           = ', ppkth
+         WRITE(numout,*) '                                        ppacr           = ', ppacr
+         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
+         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
+         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
+         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
+         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
+         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
+      ENDIF
+
+      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
+      e3zps_min = rn_e3zps_min
+      e3zps_rat = rn_e3zps_rat
+      nmsh      = nn_msh
+      nacc      = nn_acc
+      atfp      = rn_atfp
+      rdt       = rn_rdt
+      rdtmin    = rn_rdtmin
+      rdtmax    = rn_rdtmin
+      rdth      = rn_rdth
+
+      REWIND( numnam_ref )              ! Namelist namcla in reference namelist : Cross land advection
+      READ  ( numnam_ref, namcla, IOSTAT = ios, ERR = 905)
+905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namcla in configuration namelist : Cross land advection
+      READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )
+906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )
+      IF(lwm) WRITE( numond, namcla )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '   Namelist namcla'
+         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
+      ENDIF
+      IF ( nn_cla .EQ. 1 ) THEN
+         IF  ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2 
+            CONTINUE
+         ELSE
+            CALL ctl_stop( 'STOP', 'Cross land advation iplemented only for ORCA2 configuration: cp_cfg = "orca" and jp_cfg = 2 ' )
+         ENDIF
+      ENDIF
+
+#if defined key_netcdf4
+      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
+      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
+      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
+907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
+      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
+908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
+      IF(lwm) WRITE( numond, namnc4 )
+
+      IF(lwp) THEN                        ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
+         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
+         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
+         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
+         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
+      ENDIF
+
+      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
+      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
+      snc4set%ni   = nn_nchunks_i
+      snc4set%nj   = nn_nchunks_j
+      snc4set%nk   = nn_nchunks_k
+      snc4set%luse = ln_nc4zip
+#else
+      snc4set%luse = .FALSE.        ! No NetCDF 4 case
+#endif
+      !
+   END SUBROUTINE dom_nam
+
+
+   SUBROUTINE dom_ctl
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dom_ctl  ***
+      !!
+      !! ** Purpose :   Domain control.
+      !!
+      !! ** Method  :   compute and print extrema of masked scale factors
+      !!----------------------------------------------------------------------
+      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
+      INTEGER, DIMENSION(2) ::   iloc   ! 
+      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
+      !!----------------------------------------------------------------------
+      !
+      IF(lk_mpp) THEN
+         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
+         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
+         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
+         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
+      ELSE
+         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )    
+         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )    
+         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )    
+         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )    
+
+         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
+         iimi1 = iloc(1) + nimpp - 1
+         ijmi1 = iloc(2) + njmpp - 1
+         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
+         iimi2 = iloc(1) + nimpp - 1
+         ijmi2 = iloc(2) + njmpp - 1
+         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
+         iima1 = iloc(1) + nimpp - 1
+         ijma1 = iloc(2) + njmpp - 1
+         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
+         iima2 = iloc(1) + nimpp - 1
+         ijma2 = iloc(2) + njmpp - 1
+      ENDIF
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
+         WRITE(numout,*) '~~~~~~~'
+         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
+         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
+         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
+         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
+      ENDIF
+      !
+   END SUBROUTINE dom_ctl
+
+   SUBROUTINE dom_stiff
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_stiff  ***
+      !!                     
+      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
+      !!
+      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
+      !!                Save the maximum in the vertical direction
+      !!                (this number is only relevant in s-coordinates)
+      !!
+      !!                Haney, R. L., 1991: On the pressure gradient force
+      !!                over steep topography in sigma coordinate ocean models. 
+      !!                J. Phys. Oceanogr., 21, 610???619.
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk 
+      REAL(wp) ::   zrxmax
+      REAL(wp), DIMENSION(4) :: zr1
+      !!----------------------------------------------------------------------
+      rx1(:,:) = 0.e0
+      zrxmax   = 0.e0
+      zr1(:)   = 0.e0
+      
+      DO ji = 2, jpim1
+         DO jj = 2, jpjm1
+            DO jk = 1, jpkm1
+               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji-1,jj  ,jk  )  & 
+                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1)) &
+                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji-1,jj  ,jk  )  &
+                    &                         -gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1) + rsmall) )
+               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw_0(ji+1,jj  ,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
+                    &                         +gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
+                    &                        /(gdepw_0(ji+1,jj  ,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
+                    &                         -gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
+               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw_0(ji  ,jj+1,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
+                    &                         +gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
+                    &                        /(gdepw_0(ji  ,jj+1,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
+                    &                         -gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
+               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji  ,jj-1,jk  )  &
+                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1)) &
+                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji  ,jj-1,jk  )  &
+                    &                         -gdepw_0(ji,  jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1) + rsmall) )
+               zrxmax = MAXVAL(zr1(1:4))
+               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
+            END DO
+         END DO
+      END DO
+
+      CALL lbc_lnk( rx1, 'T', 1. )
+
+      zrxmax = MAXVAL(rx1)
+
+      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
+         WRITE(numout,*) '~~~~~~~~~'
+      ENDIF
+
+   END SUBROUTINE dom_stiff
+
+
+
+   !!======================================================================
+END MODULE domain
Index: /CONFIG/UNIFORM/v6/IPSLESM6/SOURCES/NEMO/in_out_manager.F90
===================================================================
--- /CONFIG/UNIFORM/v6/IPSLESM6/SOURCES/NEMO/in_out_manager.F90	(revision 3695)
+++ /CONFIG/UNIFORM/v6/IPSLESM6/SOURCES/NEMO/in_out_manager.F90	(revision 3695)
@@ -0,0 +1,161 @@
+MODULE in_out_manager   
+   !!======================================================================
+   !!                       ***  MODULE  in_out_manager  ***
+   !! I/O manager utilities : Defines run parameters together with logical units
+   !!=====================================================================
+   !! History :   1.0  !  2002-06  (G. Madec)   original code
+   !!             2.0  !  2006-07  (S. Masson)  iom, add ctl_stop, ctl_warn
+   !!             3.0  !  2008-06  (G. Madec)   add ctmp4 to ctmp10
+   !!             3.2  !  2009-08  (S. MAsson)  add new ctl_opn
+   !!             3.3  !  2010-10  (A. Coward)  add NetCDF4 usage
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   USE par_oce       ! ocean parameter
+   USE lib_print     ! formated print library
+   USE nc4interface  ! NetCDF4 interface
+
+   IMPLICIT NONE
+   PUBLIC
+
+ 
+   !
+   !!----------------------------------------------------------------------
+   !!                   namrun namelist parameters
+   !!----------------------------------------------------------------------
+   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_rstart_ts     !: start from (F) rest for current only or (T) 
+   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)
+   INTEGER       ::   nn_rstssh   = 0  !: hand made initilization of ssh or not (1/0)
+   INTEGER       ::   nn_it000         !: index of the first time step
+   INTEGER       ::   nn_itend         !: index of the last time step
+   INTEGER       ::   nn_date0         !: initial calendar date aammjj
+   INTEGER       ::   nn_leapy         !: Leap year calendar flag (0/1 or 30)
+   INTEGER       ::   nn_istate        !: initial state output flag (0/1)
+   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_mskutil       !: output data excluding halos
+   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)
+#if defined key_netcdf4
+   !!----------------------------------------------------------------------
+   !!                   namnc4 namelist parameters                         (key_netcdf4)
+   !!----------------------------------------------------------------------
+   ! The following four values determine the partitioning of the output fields
+   ! into netcdf4 chunks. They are unrelated to the nn_chunk_sz setting which is
+   ! for runtime optimisation. The individual netcdf4 chunks can be optionally 
+   ! gzipped (recommended) leading to significant reductions in I/O volumes 
+   !                         !!!**  variables only used with iom_nf90 routines and key_netcdf4 **
+   INTEGER ::   nn_nchunks_i   !: number of chunks required in the i-dimension 
+   INTEGER ::   nn_nchunks_j   !: number of chunks required in the j-dimension 
+   INTEGER ::   nn_nchunks_k   !: number of chunks required in the k-dimension 
+   INTEGER ::   nn_nchunks_t   !: number of chunks required in the t-dimension 
+   LOGICAL ::   ln_nc4zip      !: netcdf4 usage: (T) chunk and compress output using the HDF5 sublayers of netcdf4
+   !                           !                 (F) ignore chunking request and use the netcdf4 library 
+   !                           !                     to produce netcdf3-compatible files 
+#endif
+!$AGRIF_DO_NOT_TREAT
+   TYPE(snc4_ctl)     :: snc4set        !: netcdf4 chunking control structure (always needed for decision making)
+!$AGRIF_END_DO_NOT_TREAT
+
+
+   !! conversion of DOCTOR norm namelist name into model name
+   !! (this should disappear in a near futur)
+
+   CHARACTER(lc) ::   cexper                      !: experiment name used for output filename
+   INTEGER       ::   no                          !: job number
+   INTEGER       ::   nrstdt                      !: control of the time step (0, 1 or 2)
+   INTEGER       ::   nit000                      !: index of the first time step
+   INTEGER       ::   nitend                      !: index of the last time step
+   INTEGER       ::   ndate0                      !: initial calendar date aammjj
+   INTEGER       ::   nleapy                      !: Leap year calendar flag (0/1 or 30)
+   INTEGER       ::   ninist                      !: initial state output flag (0/1)
+   INTEGER       ::   nwrite                      !: model standard output frequency
+   INTEGER       ::   nstock                      !: restart file frequency
+   INTEGER, DIMENSION(10) :: nstocklist           !: restart dump times
+
+   !!----------------------------------------------------------------------
+   !! was in restart but moved here because of the OFF line... better solution should be found...
+   !!----------------------------------------------------------------------
+   INTEGER ::   nitrst                !: time step at which restart file should be written
+   LOGICAL ::   lrst_oce              !: logical to control the oce restart write 
+   INTEGER ::   numror = 0            !: logical unit for ocean restart (read). Init to 0 is needed for SAS (in daymod.F90)
+   INTEGER ::   numrow                !: logical unit for ocean restart (write)
+   INTEGER ::   nrst_lst              !: number of restart to output next
+
+   !!----------------------------------------------------------------------
+   !!                    output monitoring
+   !!----------------------------------------------------------------------
+   LOGICAL ::   ln_ctl       !: run control for debugging
+   INTEGER ::   nn_timing    !: run control for timing
+   INTEGER ::   nn_print     !: level of print (0 no print)
+   INTEGER ::   nn_ictls     !: Start i indice for the SUM control
+   INTEGER ::   nn_ictle     !: End   i indice for the SUM control
+   INTEGER ::   nn_jctls     !: Start j indice for the SUM control
+   INTEGER ::   nn_jctle     !: End   j indice for the SUM control
+   INTEGER ::   nn_isplt     !: number of processors following i
+   INTEGER ::   nn_jsplt     !: number of processors following j
+   INTEGER ::   nn_bench     !: benchmark parameter (0/1)
+   INTEGER ::   nn_bit_cmp   =    0    !: bit reproducibility  (0/1)
+
+   !                                          
+   INTEGER ::   nprint, nictls, nictle, njctls, njctle, isplt, jsplt, nbench    !: OLD namelist names
+
+   INTEGER ::   ijsplt     =    1      !: nb of local domain = nb of processors
+
+   !!----------------------------------------------------------------------
+   !!                        logical units
+   !!----------------------------------------------------------------------
+   INTEGER ::   numstp          =   -1      !: logical unit for time step
+   INTEGER ::   numtime         =   -1      !: logical unit for timing
+   INTEGER ::   numout          =    6      !: logical unit for output print; Set to stdout to ensure any early
+                                            !  output can be collected; do not change
+   INTEGER ::   numnam_ref      =   -1      !: logical unit for reference namelist
+   INTEGER ::   numnam_cfg      =   -1      !: logical unit for configuration specific namelist
+   INTEGER ::   numond          =   -1      !: logical unit for Output Namelist Dynamics
+   INTEGER ::   numnam_ice_ref  =   -1      !: logical unit for ice reference namelist
+   INTEGER ::   numnam_ice_cfg  =   -1      !: logical unit for ice reference namelist
+   INTEGER ::   numoni          =   -1      !: logical unit for Output Namelist Ice
+   INTEGER ::   numevo_ice      =   -1      !: logical unit for ice variables (temp. evolution)
+   INTEGER ::   numsol          =   -1      !: logical unit for solver statistics
+   INTEGER ::   numdct_in       =   -1      !: logical unit for transports computing
+   INTEGER ::   numdct_vol      =   -1      !: logical unit for voulume transports output
+   INTEGER ::   numdct_heat     =   -1      !: logical unit for heat    transports output
+   INTEGER ::   numdct_salt     =   -1      !: logical unit for salt    transports output
+   INTEGER ::   numfl           =   -1      !: logical unit for floats ascii output
+   INTEGER ::   numflo          =   -1      !: logical unit for floats ascii output
+
+   !!----------------------------------------------------------------------
+   !!                          Run control  
+   !!----------------------------------------------------------------------
+   INTEGER       ::   nstop = 0             !: error flag (=number of reason for a premature stop run)
+   INTEGER       ::   nwarn = 0             !: warning flag (=number of warning found during the run)
+   CHARACTER(lc) ::   ctmp1, ctmp2, ctmp3   !: temporary characters 1 to 3
+   CHARACTER(lc) ::   ctmp4, ctmp5, ctmp6   !: temporary characters 4 to 6
+   CHARACTER(lc) ::   ctmp7, ctmp8, ctmp9   !: temporary characters 7 to 9
+   CHARACTER(lc) ::   ctmp10                !: temporary character 10
+   CHARACTER(lc) ::   cform_err = "(/,' ===>>> : E R R O R',     /,'         ===========',/)"       !:
+   CHARACTER(lc) ::   cform_war = "(/,' ===>>> : W A R N I N G', /,'         ===============',/)"   !:
+   LOGICAL       ::   lwm      = .FALSE.    !: boolean : true on the 1st processor only (always)
+   LOGICAL       ::   lwp      = .FALSE.    !: boolean : true on the 1st processor only .OR. ln_ctl
+   LOGICAL       ::   lsp_area = .TRUE.     !: to make a control print over a specific area
+   CHARACTER(lc) ::   cxios_context         !: context name used in xios
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: in_out_manager.F90 5518 2015-06-30 13:11:42Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!=====================================================================
+END MODULE in_out_manager
Index: /CONFIG/UNIFORM/v6/IPSLESM6/SOURCES/NEMO/restart.F90
===================================================================
--- /CONFIG/UNIFORM/v6/IPSLESM6/SOURCES/NEMO/restart.F90	(revision 3695)
+++ /CONFIG/UNIFORM/v6/IPSLESM6/SOURCES/NEMO/restart.F90	(revision 3695)
@@ -0,0 +1,287 @@
+MODULE restart
+   !!======================================================================
+   !!                     ***  MODULE  restart  ***
+   !! Ocean restart :  write the ocean restart file
+   !!======================================================================
+   !! History :  OPA  !  1999-11  (M. Imbard)  Original code
+   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form
+   !!            2.0  !  2006-07  (S. Masson)  use IOM for restart
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
+   !!            - -  !  2010-10  (C. Ethe, G. Madec) TRC-TRA merge (T-S in 4D)
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   rst_opn    : open the ocean restart file
+   !!   rst_write  : write the ocean restart file
+   !!   rst_read   : read the ocean restart file
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers 
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O module
+   USE eosbn2          ! equation of state            (eos bn2 routine)
+   USE trdmxl_oce      ! ocean active mixed layer tracers trends variables
+   USE divcur          ! hor. divergence and curl      (div & cur routines)
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   rst_opn         ! routine called by step module
+   PUBLIC   rst_write       ! routine called by step module
+   PUBLIC   rst_read        ! routine called by istate module
+   PUBLIC   rst_read_open   ! routine called in rst_read and (possibly) in dom_vvl_init
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: restart.F90 5407 2015-06-11 19:13:22Z smasson $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE rst_opn( kt )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE rst_opn  ***
+      !!                     
+      !! ** Purpose : + initialization (should be read in the namelist) of nitrst 
+      !!              + open the restart when we are one time step before nitrst
+      !!                   - restart header is defined when kt = nitrst-1
+      !!                   - restart data  are written when kt = nitrst
+      !!              + define lrst_oce to .TRUE. when we need to define or write the restart
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt     ! ocean time-step
+      !!
+      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character
+      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
+      CHARACTER(lc)       ::   clpath   ! full path to ocean output restart file
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 ) THEN   ! default definitions
+         lrst_oce = .FALSE.   
+         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
+         IF( nitrst > nitend )   nitrst = nitend   ! make sure we write a restart at the end of the run
+      ENDIF
+      ! to get better performances with NetCDF format:
+      ! we open and define the ocean restart file one time step before writing the data (-> at nitrst - 1)
+      ! 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
+         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
+            ! 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
+      !
+   END SUBROUTINE rst_opn
+
+
+   SUBROUTINE rst_write( kt )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE rstwrite  ***
+      !!                     
+      !! ** Purpose :   Write restart fields in the format corresponding to jprstlib
+      !!
+      !! ** Method  :   Write in numrow when kt == nitrst in NetCDF
+      !!              file, save fields which are necessary for restart
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step
+      !!----------------------------------------------------------------------
+
+                     CALL iom_rstput( kt, nitrst, numrow, 'rdt'    , rdt       )   ! dynamics time step
+                     CALL iom_rstput( kt, nitrst, numrow, 'rdttra1', rdttra(1) )   ! surface tracer time step
+
+                     CALL iom_rstput( kt, nitrst, numrow, 'ub'     , ub        )     ! before fields
+                     CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vb        )
+                     CALL iom_rstput( kt, nitrst, numrow, 'tb'     , tsb(:,:,:,jp_tem) )
+                     CALL iom_rstput( kt, nitrst, numrow, 'sb'     , tsb(:,:,:,jp_sal) )
+                     CALL iom_rstput( kt, nitrst, numrow, 'rotb'   , rotb      )
+                     CALL iom_rstput( kt, nitrst, numrow, 'hdivb'  , hdivb     )
+                     CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb      )
+                     !
+                     CALL iom_rstput( kt, nitrst, numrow, 'un'     , un        )     ! now fields
+                     CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vn        )
+                     CALL iom_rstput( kt, nitrst, numrow, 'tn'     , tsn(:,:,:,jp_tem) )
+                     CALL iom_rstput( kt, nitrst, numrow, 'sn'     , tsn(:,:,:,jp_sal) )
+                     CALL iom_rstput( kt, nitrst, numrow, 'rotn'   , rotn      )
+                     CALL iom_rstput( kt, nitrst, numrow, 'hdivn'  , hdivn     )
+                     CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , sshn      )
+                     CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop      )
+#if defined key_zdfkpp
+                     CALL iom_rstput( kt, nitrst, numrow, 'rhd'    , rhd       )
+#endif
+      IF( kt == nitrst ) THEN
+         CALL iom_close( numrow )     ! close the restart file (only at last time step)
+!!gm         IF( .NOT. lk_trdmld )   lrst_oce = .FALSE.
+!!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
+      !
+   END SUBROUTINE rst_write
+
+
+   SUBROUTINE rst_read_open
+      !!---------------------------------------------------------------------- 
+      !!                   ***  ROUTINE rst_read_open  ***
+      !! 
+      !! ** Purpose :   Open read files for restart (format fixed by jprstlib )
+      !! 
+      !! ** Method  :   Use a non-zero, positive value of numror to assess whether or not
+      !!                the file has already been opened
+      !!----------------------------------------------------------------------
+      INTEGER        ::   jlibalt = jprstlib
+      LOGICAL        ::   llok
+      CHARACTER(lc)  ::   clpath   ! full path to ocean output restart file
+      !!----------------------------------------------------------------------
+      !
+      IF( numror <= 0 ) THEN
+         IF(lwp) THEN                                             ! Contol prints
+            WRITE(numout,*)
+            SELECT CASE ( jprstlib )
+            CASE ( jpnf90    )   ;   WRITE(numout,*) 'rst_read : read oce NetCDF restart file'
+            CASE ( jprstdimg )   ;   WRITE(numout,*) 'rst_read : read oce binary restart file'
+            END SELECT
+            IF ( snc4set%luse )      WRITE(numout,*) 'rst_read : configured with NetCDF4 support'
+            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_ocerst_in}.nc exists, then set jlibalt to jpnf90
+           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( TRIM(clpath)//cn_ocerst_in, numror, kiolib = jlibalt )
+      ENDIF
+   END SUBROUTINE rst_read_open
+
+   SUBROUTINE rst_read
+      !!---------------------------------------------------------------------- 
+      !!                   ***  ROUTINE rst_read  ***
+      !! 
+      !! ** Purpose :   Read files for restart (format fixed by jprstlib )
+      !! 
+      !! ** Method  :   Read in restart.nc file fields which are necessary for restart
+      !!----------------------------------------------------------------------
+      REAL(wp) ::   zrdt, zrdttra1
+      INTEGER  ::   jk
+      LOGICAL  ::   llok
+      !!----------------------------------------------------------------------
+
+      CALL rst_read_open           ! open restart for reading (if not already opened)
+
+      ! Check dynamics and tracer time-step consistency and force Euler restart if changed
+      IF( iom_varid( numror, 'rdt', ldstop = .FALSE. ) > 0 )   THEN
+         CALL iom_get( numror, 'rdt', zrdt )
+         IF( zrdt /= rdt )   neuler = 0
+      ENDIF
+      IF( iom_varid( numror, 'rdttra1', ldstop = .FALSE. ) > 0 )   THEN
+         CALL iom_get( numror, 'rdttra1', zrdttra1 )
+         IF( zrdttra1 /= rdttra(1) )   neuler = 0
+      ENDIF
+      ! 
+      IF( iom_varid( numror, 'ub', ldstop = .FALSE. ) > 0 ) THEN
+         CALL iom_get( numror, jpdom_autoglo, 'ub'     , ub      )   ! before fields
+         CALL iom_get( numror, jpdom_autoglo, 'vb'     , vb      )
+         CALL iom_get( numror, jpdom_autoglo, 'tb'     , tsb(:,:,:,jp_tem) )
+         CALL iom_get( numror, jpdom_autoglo, 'sb'     , tsb(:,:,:,jp_sal) )
+         CALL iom_get( numror, jpdom_autoglo, 'rotb'   , rotb    )
+         CALL iom_get( numror, jpdom_autoglo, 'hdivb'  , hdivb   )
+         CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb    )
+      ELSE
+         neuler = 0
+      ENDIF
+      !
+      CALL iom_get( numror, jpdom_autoglo, 'un'     , un      )   ! now    fields
+      CALL iom_get( numror, jpdom_autoglo, 'vn'     , vn      )
+      CALL iom_get( numror, jpdom_autoglo, 'tn'     , tsn(:,:,:,jp_tem) )
+      CALL iom_get( numror, jpdom_autoglo, 'sn'     , tsn(:,:,:,jp_sal) )
+      CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    )
+      IF( iom_varid( numror, 'rotn', ldstop = .FALSE. ) > 0 ) THEN
+         CALL iom_get( numror, jpdom_autoglo, 'rotn'   , rotn    )
+         CALL iom_get( numror, jpdom_autoglo, 'hdivn'  , hdivn   )
+      ELSE
+         CALL div_cur( 0 )                              ! Horizontal divergence & Relative vorticity
+      ENDIF
+      IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN
+         CALL iom_get( numror, jpdom_autoglo, 'rhop'   , rhop    )   ! now    potential density
+      ELSE
+         CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )   
+      ENDIF
+#if defined key_zdfkpp
+      IF( iom_varid( numror, 'rhd', ldstop = .FALSE. ) > 0 ) THEN
+         CALL iom_get( numror, jpdom_autoglo, 'rhd'    , rhd     )   ! now    in situ density anomaly
+      ELSE
+         CALL eos( tsn, rhd, fsdept_n(:,:,:) )   ! compute rhd
+      ENDIF
+#endif
+      IF( ln_rstart_ts ) THEN    ! ocean is forced to rest
+         neuler = 0
+         sshn (:,:)   = 0._wp
+         un   (:,:,:) = 0._wp
+         vn   (:,:,:) = 0._wp  
+         rotn (:,:,:) = 0._wp
+         hdivn(:,:,:) = 0._wp
+      ENDIF
+      !
+      IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0)
+         tsb  (:,:,:,:) = tsn  (:,:,:,:)                             ! all before fields set to now values
+         ub   (:,:,:)   = un   (:,:,:)
+         vb   (:,:,:)   = vn   (:,:,:)
+         rotb (:,:,:)   = rotn (:,:,:)
+         hdivb(:,:,:)   = hdivn(:,:,:)
+         sshb (:,:)     = sshn (:,:)
+
+         IF( lk_vvl ) THEN
+            DO jk = 1, jpk
+               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
+            END DO
+         ENDIF
+
+      ENDIF
+      !
+   END SUBROUTINE rst_read
+
+   !!=====================================================================
+END MODULE restart
Index: /CONFIG/UNIFORM/v6/IPSLESM6/SOURCES/NEMO/sbcrnf.F90
===================================================================
--- /CONFIG/UNIFORM/v6/IPSLESM6/SOURCES/NEMO/sbcrnf.F90	(revision 3695)
+++ /CONFIG/UNIFORM/v6/IPSLESM6/SOURCES/NEMO/sbcrnf.F90	(revision 3695)
@@ -0,0 +1,530 @@
+MODULE sbcrnf
+   !!======================================================================
+   !!                       ***  MODULE  sbcrnf  ***
+   !! Ocean forcing:  river runoff
+   !!=====================================================================
+   !! History :  OPA  ! 2000-11  (R. Hordoir, E. Durand)  NetCDF FORMAT
+   !!   NEMO     1.0  ! 2002-09  (G. Madec)  F90: Free form and module
+   !!            3.0  ! 2006-07  (G. Madec)  Surface module
+   !!            3.2  ! 2009-04  (B. Lemaire)  Introduce iom_put
+   !!            3.3  ! 2010-10  (R. Furner, G. Madec) runoff distributed over ocean levels
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_rnf      : monthly runoffs read in a NetCDF file
+   !!   sbc_rnf_init : runoffs initialisation
+   !!   rnf_mouth    : set river mouth mask
+   !!----------------------------------------------------------------------
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE sbc_oce         ! surface boundary condition variables
+   USE sbcisf          ! PM we could remove it I think
+   USE closea          ! closed seas
+   USE fldread         ! read input field at current time step
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O module
+   USE lib_mpp         ! MPP library
+   USE eosbn2
+   USE wrk_nemo        ! Memory allocation
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc_rnf       ! routine call in sbcmod module
+   PUBLIC   sbc_rnf_div   ! routine called in divcurl module
+   PUBLIC   sbc_rnf_alloc ! routine call in sbcmod module
+   PUBLIC   sbc_rnf_init  ! (PUBLIC for TAM)
+   !                                                     !!* namsbc_rnf namelist *
+   CHARACTER(len=100)         ::   cn_dir          !: Root directory for location of rnf files
+   LOGICAL                    ::   ln_rnf_depth      !: depth       river runoffs attribute specified in a file
+   LOGICAL                    ::   ln_rnf_depth_ini  !: depth       river runoffs  computed at the initialisation
+   REAL(wp)                   ::   rn_rnf_max        !: maximum value of the runoff climatologie ( ln_rnf_depth_ini = .true )
+   REAL(wp)                   ::   rn_dep_max        !: depth over which runoffs is spread ( ln_rnf_depth_ini = .true )
+   INTEGER                    ::   nn_rnf_depth_file !: create (=1) a runoff depth file or not (=0)
+   LOGICAL                    ::   ln_rnf_tem      !: temperature river runoffs attribute specified in a file
+   LOGICAL           , PUBLIC ::   ln_rnf_sal      !: salinity    river runoffs attribute specified in a file
+   TYPE(FLD_N)       , PUBLIC ::   sn_rnf          !: information about the runoff file to be read
+   TYPE(FLD_N)                ::   sn_cnf          !: information about the runoff mouth file to be read
+   TYPE(FLD_N)                ::   sn_s_rnf        !: information about the salinities of runoff file to be read
+   TYPE(FLD_N)                ::   sn_t_rnf        !: information about the temperatures of runoff file to be read
+   TYPE(FLD_N)                ::   sn_dep_rnf      !: information about the depth which river inflow affects
+   LOGICAL           , PUBLIC ::   ln_rnf_mouth    !: specific treatment in mouths vicinity
+   REAL(wp)                   ::   rn_hrnf         !: runoffs, depth over which enhanced vertical mixing is used
+   REAL(wp)          , PUBLIC ::   rn_avt_rnf      !: runoffs, value of the additional vertical mixing coef. [m2/s]
+   REAL(wp)          , PUBLIC ::   rn_rfact        !: multiplicative factor for runoff
+
+   LOGICAL           , PUBLIC ::   l_rnfcpl = .false.       ! runoffs recieved from oasis
+
+   INTEGER , PUBLIC  ::   nkrnf = 0         !: nb of levels over which Kz is increased at river mouths
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rnfmsk              !: river mouth mask (hori.)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     ::   rnfmsk_z            !: river mouth mask (vert.)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   h_rnf               !: depth of runoff in m
+   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nk_rnf              !: depth of runoff in model levels
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rnf_tsc_b, rnf_tsc  !: before and now T & S runoff contents   [K.m/s & PSU.m/s]   
+
+   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_rnf       ! structure: river runoff (file information, fields read)
+   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_s_rnf     ! structure: river runoff salinity (file information, fields read)  
+   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_t_rnf     ! structure: river runoff temperature (file information, fields read)  
+ 
+   !! * Substitutions  
+#  include "domzgr_substitute.h90"  
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: sbcrnf.F90 7963 2017-04-24 16:28:04Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION sbc_rnf_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE sbc_rnf_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( rnfmsk(jpi,jpj)         , rnfmsk_z(jpk)          ,     &
+         &      h_rnf (jpi,jpj)         , nk_rnf  (jpi,jpj)      ,     &
+         &      rnf_tsc_b(jpi,jpj,jpts) , rnf_tsc (jpi,jpj,jpts) , STAT=sbc_rnf_alloc )
+         !
+      IF( lk_mpp            )   CALL mpp_sum ( sbc_rnf_alloc )
+      IF( sbc_rnf_alloc > 0 )   CALL ctl_warn('sbc_rnf_alloc: allocation of arrays failed')
+   END FUNCTION sbc_rnf_alloc
+
+
+   SUBROUTINE sbc_rnf( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_rnf  ***
+      !!
+      !! ** Purpose :   Introduce a climatological run off forcing
+      !!
+      !! ** Method  :   Set each river mouth with a monthly climatology
+      !!                provided from different data.
+      !!                CAUTION : upward water flux, runoff forced to be < 0
+      !!
+      !! ** Action  :   runoff updated runoff field at time-step kt
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt          ! ocean time step
+      !
+      INTEGER  ::   ji, jj    ! dummy loop indices
+      INTEGER  ::   z_err = 0 ! dummy integer for error handling
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), POINTER       ::   ztfrz   ! freezing point used for temperature correction
+      !
+      CALL wrk_alloc( jpi,jpj, ztfrz)
+      !
+      !                                            !-------------------!
+      !                                            !   Update runoff   !
+      !                                            !-------------------!
+      !
+      IF( .NOT. l_rnfcpl )   CALL fld_read ( kt, nn_fsbc, sf_rnf   )    ! Read Runoffs data and provide it at kt
+      IF(   ln_rnf_tem   )   CALL fld_read ( kt, nn_fsbc, sf_t_rnf )    ! idem for runoffs temperature if required
+      IF(   ln_rnf_sal   )   CALL fld_read ( kt, nn_fsbc, sf_s_rnf )    ! idem for runoffs salinity    if required
+      !
+      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
+         !
+         IF( .NOT. l_rnfcpl )   rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )       ! updated runoff value at time step kt
+         !
+         !                                                     ! set temperature & salinity content of runoffs
+         IF( ln_rnf_tem ) THEN                                       ! use runoffs temperature data
+            rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
+            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp )             ! if missing data value use SST as runoffs temperature
+               rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0
+            END WHERE
+            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg
+               ztfrz(:,:) = -1.9 !tfreez( sss_m(:,:) ) !PM to be discuss (trouble if sensitivity study)
+               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * lfusisf * r1_rau0_rcp
+            END WHERE
+         ELSE                                                        ! use SST as runoffs temperature
+!            rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0
+            rnf_tsc(:,:,jp_tem) = MAX( sst_m(:,:), 0._wp ) * rnf(:,:) * r1_rau0
+         ENDIF
+         !                                                           ! use runoffs salinity data
+         IF( ln_rnf_sal )   rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
+         !                                                           ! else use S=0 for runoffs (done one for all in the init)
+         IF( iom_use('runoffs') )        CALL iom_put( 'runoffs'     , rnf(:,:)                         )   ! output runoff mass flux
+         IF( iom_use('hflx_rnf_cea') )   CALL iom_put( 'hflx_rnf_cea', rnf_tsc(:,:,jp_tem) * rau0 * rcp )   ! output runoff sensible heat (W/m2)
+      ENDIF
+      !
+      !                                                ! ---------------------------------------- !
+      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
+         !                                             ! ---------------------------------------- !
+         IF( ln_rstart .AND.    &                               !* Restart: read in restart file
+            & iom_varid( numror, 'rnf_b', ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields red in the restart file'
+            CALL iom_get( numror, jpdom_autoglo, 'rnf_b', rnf_b )     ! before runoff
+            CALL iom_get( numror, jpdom_autoglo, 'rnf_hc_b', rnf_tsc_b(:,:,jp_tem) )   ! before heat content of runoff
+            CALL iom_get( numror, jpdom_autoglo, 'rnf_sc_b', rnf_tsc_b(:,:,jp_sal) )   ! before salinity content of runoff
+         ELSE                                                   !* no restart: set from nit000 values
+            IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields set to nit000'
+            rnf_b    (:,:  ) = rnf    (:,:  )
+            rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
+         ENDIF
+      ENDIF
+      !                                                ! ---------------------------------------- !
+      IF( lrst_oce ) THEN                              !      Write in the ocean restart file     !
+         !                                             ! ---------------------------------------- !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'sbcrnf : runoff forcing fields written in ocean restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         CALL iom_rstput( kt, nitrst, numrow, 'rnf_b' , rnf )
+         CALL iom_rstput( kt, nitrst, numrow, 'rnf_hc_b', rnf_tsc(:,:,jp_tem) )
+         CALL iom_rstput( kt, nitrst, numrow, 'rnf_sc_b', rnf_tsc(:,:,jp_sal) )
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi,jpj, ztfrz)
+      !
+   END SUBROUTINE sbc_rnf
+
+
+   SUBROUTINE sbc_rnf_div( phdivn )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_rnf  ***
+      !!
+      !! ** Purpose :   update the horizontal divergence with the runoff inflow
+      !!
+      !! ** Method  :
+      !!                CAUTION : rnf is positive (inflow) decreasing the
+      !!                          divergence and expressed in m/s
+      !!
+      !! ** Action  :   phdivn   decreased by the runoff inflow
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zfact     ! local scalar
+      !!----------------------------------------------------------------------
+      !
+      zfact = 0.5_wp
+      !
+      IF( ln_rnf_depth .OR. ln_rnf_depth_ini ) THEN      !==   runoff distributed over several levels   ==!
+         IF( lk_vvl ) THEN             ! variable volume case
+            DO jj = 1, jpj                   ! update the depth over which runoffs are distributed
+               DO ji = 1, jpi
+                  h_rnf(ji,jj) = 0._wp
+                  DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres
+                     h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)   ! to the bottom of the relevant grid box
+                  END DO
+                  !                          ! apply the runoff input flow
+                  DO jk = 1, nk_rnf(ji,jj)
+                     phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rau0 / h_rnf(ji,jj)
+                  END DO
+               END DO
+            END DO
+         ELSE                          ! constant volume case : just apply the runoff input flow
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  DO jk = 1, nk_rnf(ji,jj)
+                     phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rau0 / h_rnf(ji,jj)
+                  END DO
+               END DO
+            END DO
+         ENDIF
+      ELSE                       !==   runoff put only at the surface   ==!
+         IF( lk_vvl ) THEN              ! variable volume case
+            h_rnf(:,:) = fse3t(:,:,1)   ! recalculate h_rnf to be depth of top box
+         ENDIF
+         phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rau0 / fse3t(:,:,1)
+      ENDIF
+      !
+   END SUBROUTINE sbc_rnf_div
+
+
+   SUBROUTINE sbc_rnf_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_rnf_init  ***
+      !!
+      !! ** Purpose :   Initialisation of the runoffs if (ln_rnf=T)
+      !!
+      !! ** Method  : - read the runoff namsbc_rnf namelist
+      !!
+      !! ** Action  : - read parameters
+      !!----------------------------------------------------------------------
+      CHARACTER(len=32) ::   rn_dep_file   ! runoff file name
+      INTEGER           ::   ji, jj, jk, jm    ! dummy loop indices
+      INTEGER           ::   ierror, inum  ! temporary integer
+      INTEGER           ::   ios           ! Local integer output status for namelist read
+      INTEGER           ::   nbrec         ! temporary integer
+      REAL(wp)          ::   zacoef  
+      REAL(wp), DIMENSION(12)                 :: zrec             ! times records
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zrnfcl    
+      REAL(wp), DIMENSION(:,:  ), ALLOCATABLE :: zrnf
+      !
+      NAMELIST/namsbc_rnf/ cn_dir            , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal,   &
+         &                 sn_rnf, sn_cnf    , sn_s_rnf    , sn_t_rnf  , sn_dep_rnf,   &
+         &                 ln_rnf_mouth      , rn_hrnf     , rn_avt_rnf, rn_rfact,     &
+         &                 ln_rnf_depth_ini  , rn_dep_max  , rn_rnf_max, nn_rnf_depth_file
+      !!----------------------------------------------------------------------
+      !
+      !                                         !==  allocate runoff arrays
+      IF( sbc_rnf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_rnf_alloc : unable to allocate arrays' )
+      !
+      IF( .NOT. ln_rnf ) THEN                      ! no specific treatment in vicinity of river mouths 
+         ln_rnf_mouth  = .FALSE.                   ! default definition needed for example by sbc_ssr or by tra_adv_muscl
+         nkrnf         = 0
+         rnf     (:,:) = 0.0_wp
+         rnf_b   (:,:) = 0.0_wp
+         rnfmsk  (:,:) = 0.0_wp
+         rnfmsk_z(:)   = 0.0_wp
+         RETURN
+      ENDIF
+      !
+      !                                   ! ============
+      !                                   !   Namelist
+      !                                   ! ============
+      !
+      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
+      READ  ( numnam_ref, namsbc_rnf, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_rnf in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
+      READ  ( numnam_cfg, namsbc_rnf, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_rnf in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namsbc_rnf )
+      !
+      !                                         ! Control print
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_rnf : runoff '
+         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) '   Namelist namsbc_rnf'
+         WRITE(numout,*) '      specific river mouths treatment            ln_rnf_mouth = ', ln_rnf_mouth
+         WRITE(numout,*) '      river mouth additional Kz                  rn_avt_rnf   = ', rn_avt_rnf
+         WRITE(numout,*) '      depth of river mouth additional mixing     rn_hrnf      = ', rn_hrnf
+         WRITE(numout,*) '      multiplicative factor for runoff           rn_rfact     = ', rn_rfact
+      ENDIF
+      !                                   ! ==================
+      !                                   !   Type of runoff
+      !                                   ! ==================
+      !
+      IF( .NOT. l_rnfcpl ) THEN                    
+         ALLOCATE( sf_rnf(1), STAT=ierror )         ! Create sf_rnf structure (runoff inflow)
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          runoffs inflow read in a file'
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'sbc_rnf: unable to allocate sf_rnf structure' )   ;   RETURN
+         ENDIF
+         ALLOCATE( sf_rnf(1)%fnow(jpi,jpj,1)   )
+         IF( sn_rnf%ln_tint ) ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,1,2) )
+         CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf' )
+      ENDIF
+      !
+      IF( ln_rnf_tem ) THEN                      ! Create (if required) sf_t_rnf structure
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          runoffs temperatures read in a file'
+         ALLOCATE( sf_t_rnf(1), STAT=ierror  )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_t_rnf structure' )   ;   RETURN
+         ENDIF
+         ALLOCATE( sf_t_rnf(1)%fnow(jpi,jpj,1)   )
+         IF( sn_t_rnf%ln_tint ) ALLOCATE( sf_t_rnf(1)%fdta(jpi,jpj,1,2) )
+         CALL fld_fill (sf_t_rnf, (/ sn_t_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf' )
+      ENDIF
+      !
+      IF( ln_rnf_sal  ) THEN                     ! Create (if required) sf_s_rnf and sf_t_rnf structures
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          runoffs salinities read in a file'
+         ALLOCATE( sf_s_rnf(1), STAT=ierror  )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_s_rnf structure' )   ;   RETURN
+         ENDIF
+         ALLOCATE( sf_s_rnf(1)%fnow(jpi,jpj,1)   )
+         IF( sn_s_rnf%ln_tint ) ALLOCATE( sf_s_rnf(1)%fdta(jpi,jpj,1,2) )
+         CALL fld_fill (sf_s_rnf, (/ sn_s_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf' )
+      ENDIF
+      !
+      IF( ln_rnf_depth ) THEN                    ! depth of runoffs set from a file
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          runoffs depth read in a file'
+         rn_dep_file = TRIM( cn_dir )//TRIM( sn_dep_rnf%clname )
+         IF( .NOT. sn_dep_rnf%ln_clim ) THEN   ;   WRITE(rn_dep_file, '(a,"_y",i4)' ) TRIM( rn_dep_file ), nyear    ! add year 
+            IF( sn_dep_rnf%cltype == 'monthly' )   WRITE(rn_dep_file, '(a,"m",i2)'  ) TRIM( rn_dep_file ), nmonth   ! add month 
+         ENDIF
+         CALL iom_open ( rn_dep_file, inum )                           ! open file
+         CALL iom_get  ( inum, jpdom_data, sn_dep_rnf%clvar, h_rnf )   ! read the river mouth array
+         CALL iom_close( inum )                                        ! close file
+         !
+         nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( h_rnf(ji,jj) > 0._wp ) THEN
+                  jk = 2
+                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
+                  END DO
+                  nk_rnf(ji,jj) = jk
+               ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1
+               ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj)
+               ELSE
+                  CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  )
+                  WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj)
+               ENDIF
+            END DO
+         END DO
+         DO jj = 1, jpj                                ! set the associated depth
+            DO ji = 1, jpi
+               h_rnf(ji,jj) = 0._wp
+               DO jk = 1, nk_rnf(ji,jj)
+                  h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
+      ELSE IF( ln_rnf_depth_ini ) THEN           ! runoffs applied at the surface
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '    depth of runoff computed once from max value of runoff'
+         IF(lwp) WRITE(numout,*) '    max value of the runoff climatologie (over global domain) rn_rnf_max = ', rn_rnf_max
+         IF(lwp) WRITE(numout,*) '    depth over which runoffs is spread                        rn_dep_max = ', rn_dep_max
+         IF(lwp) WRITE(numout,*) '     create (=1) a runoff depth file or not (=0)      nn_rnf_depth_file  = ', nn_rnf_depth_file
+
+         CALL iom_open( TRIM( sn_rnf%clname ), inum )    !  open runoff file
+         CALL iom_gettime( inum, zrec, kntime=nbrec)
+         ALLOCATE( zrnfcl(jpi,jpj,nbrec) )     ;      ALLOCATE( zrnf(jpi,jpj) )
+         DO jm = 1, nbrec
+            CALL iom_get( inum, jpdom_data, TRIM( sn_rnf%clvar ), zrnfcl(:,:,jm), jm )
+         END DO
+         CALL iom_close( inum )
+         zrnf(:,:) = MAXVAL( zrnfcl(:,:,:), DIM=3 )   !  maximum value in time
+         DEALLOCATE( zrnfcl )
+         !
+         h_rnf(:,:) = 1.
+         !
+         zacoef = rn_dep_max / rn_rnf_max            ! coef of linear relation between runoff and its depth (150m for max of runoff)
+         !
+         WHERE( zrnf(:,:) > 0._wp )  h_rnf(:,:) = zacoef * zrnf(:,:)   ! compute depth for all runoffs
+         !
+         DO jj = 1, jpj                     ! take in account min depth of ocean rn_hmin
+            DO ji = 1, jpi
+               IF( zrnf(ji,jj) > 0._wp ) THEN
+                  jk = mbkt(ji,jj)
+                  h_rnf(ji,jj) = MIN( h_rnf(ji,jj), gdept_0(ji,jj,jk ) )
+               ENDIF
+            END DO
+         END DO
+         !
+         nk_rnf(:,:) = 0                       ! number of levels on which runoffs are distributed
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( zrnf(ji,jj) > 0._wp ) THEN
+                  jk = 2
+                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
+                  END DO
+                  nk_rnf(ji,jj) = jk
+               ELSE
+                  nk_rnf(ji,jj) = 1
+               ENDIF
+            END DO
+         END DO
+         !
+         DEALLOCATE( zrnf )
+         !
+         DO jj = 1, jpj                                ! set the associated depth
+            DO ji = 1, jpi
+               h_rnf(ji,jj) = 0._wp
+               DO jk = 1, nk_rnf(ji,jj)
+                  h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
+         IF( nn_rnf_depth_file == 1 ) THEN      !  save  output nb levels for runoff
+            IF(lwp) WRITE(numout,*) '              create runoff depht file'
+            CALL iom_open  ( TRIM( sn_dep_rnf%clname ), inum, ldwrt = .TRUE., kiolib = jprstlib )
+            CALL iom_rstput( 0, 0, inum, 'rodepth', h_rnf )
+            CALL iom_close ( inum )
+         ENDIF
+      ELSE                                       ! runoffs applied at the surface
+         nk_rnf(:,:) = 1
+         h_rnf (:,:) = fse3t(:,:,1)
+      ENDIF
+      !
+      rnf(:,:) =  0._wp                         ! runoff initialisation
+      rnf_tsc(:,:,:) = 0._wp                    ! runoffs temperature & salinty contents initilisation
+      !
+      !                                   ! ========================
+      !                                   !   River mouth vicinity
+      !                                   ! ========================
+      !
+      IF( ln_rnf_mouth ) THEN                   ! Specific treatment in vicinity of river mouths :
+         !                                      !    - Increase Kz in surface layers ( rn_hrnf > 0 )
+         !                                      !    - set to zero SSS damping (ln_ssr=T)
+         !                                      !    - mixed upstream-centered (ln_traadv_cen2=T)
+         !
+         IF ( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   &
+            &                                              'be spread through depth by ln_rnf_depth'               )
+         !
+         nkrnf = 0                                  ! Number of level over which Kz increase
+         IF( rn_hrnf > 0._wp ) THEN
+            nkrnf = 2
+            DO WHILE( nkrnf /= jpkm1 .AND. gdepw_1d(nkrnf+1) < rn_hrnf )   ;   nkrnf = nkrnf + 1
+            END DO
+            IF( ln_sco )   CALL ctl_warn( 'sbc_rnf: number of levels over which Kz is increased is computed for zco...' )
+         ENDIF
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          Specific treatment used in vicinity of river mouths :'
+         IF(lwp) WRITE(numout,*) '             - Increase Kz in surface layers (if rn_hrnf > 0 )'
+         IF(lwp) WRITE(numout,*) '               by ', rn_avt_rnf,' m2/s  over ', nkrnf, ' w-levels'
+         IF(lwp) WRITE(numout,*) '             - set to zero SSS damping       (if ln_ssr=T)'
+         IF(lwp) WRITE(numout,*) '             - mixed upstream-centered       (if ln_traadv_cen2=T)'
+         !
+         CALL rnf_mouth                             ! set river mouth mask
+         !
+      ELSE                                      ! No treatment at river mouths
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          No specific treatment at river mouths'
+         rnfmsk  (:,:) = 0._wp
+         rnfmsk_z(:)   = 0._wp
+         nkrnf = 0
+      ENDIF
+      !
+   END SUBROUTINE sbc_rnf_init
+
+
+   SUBROUTINE rnf_mouth
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE rnf_mouth  ***
+      !!
+      !! ** Purpose :   define the river mouths mask
+      !!
+      !! ** Method  :   read the river mouth mask (=0/1) in the river runoff
+      !!                climatological file. Defined a given vertical structure.
+      !!                CAUTION, the vertical structure is hard coded on the
+      !!                first 5 levels.
+      !!                This fields can be used to:
+      !!                 - set an upstream advection scheme
+      !!                   (ln_rnf_mouth=T and ln_traadv_cen2=T)
+      !!                 - increase vertical on the top nn_krnf vertical levels
+      !!                   at river runoff input grid point (nn_krnf>=2, see step.F90)
+      !!                 - set to zero SSS restoring flux at river mouth grid points
+      !!
+      !! ** Action  :   rnfmsk   set to 1 at river runoff input, 0 elsewhere
+      !!                rnfmsk_z vertical structure
+      !!----------------------------------------------------------------------
+      INTEGER            ::   inum        ! temporary integers
+      CHARACTER(len=140) ::   cl_rnfile   ! runoff file name
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'rnf_mouth : river mouth mask'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~ '
+      !
+      cl_rnfile = TRIM( cn_dir )//TRIM( sn_cnf%clname )
+      IF( .NOT. sn_cnf%ln_clim ) THEN   ;   WRITE(cl_rnfile, '(a,"_y",i4)' ) TRIM( cl_rnfile ), nyear    ! add year
+         IF( sn_cnf%cltype == 'monthly' )   WRITE(cl_rnfile, '(a,"m",i2)'  ) TRIM( cl_rnfile ), nmonth   ! add month
+      ENDIF
+      !
+      ! horizontal mask (read in NetCDF file)
+      CALL iom_open ( cl_rnfile, inum )                           ! open file
+      CALL iom_get  ( inum, jpdom_data, sn_cnf%clvar, rnfmsk )    ! read the river mouth array
+      CALL iom_close( inum )                                      ! close file
+      !
+      IF( nn_closea == 1 )   CALL clo_rnf( rnfmsk )               ! closed sea inflow set as ruver mouth
+      !
+      rnfmsk_z(:)   = 0._wp                                       ! vertical structure
+      rnfmsk_z(1)   = 1.0
+      rnfmsk_z(2)   = 1.0                                         ! **********
+      rnfmsk_z(3)   = 0.5                                         ! HARD CODED on the 5 first levels
+      rnfmsk_z(4)   = 0.25                                        ! **********
+      rnfmsk_z(5)   = 0.125
+      !
+   END SUBROUTINE rnf_mouth
+
+   !!======================================================================
+END MODULE sbcrnf
