Index: /NEMO/branches/2019/ENHANCE-03_domcfg/namelist_ref
===================================================================
--- /NEMO/branches/2019/ENHANCE-03_domcfg/namelist_ref	(revision 11200)
+++ /NEMO/branches/2019/ENHANCE-03_domcfg/namelist_ref	(revision 11201)
@@ -42,6 +42,5 @@
    nn_bathy    =    1      !  compute (=0) or read (=1) the bathymetry file
    rn_bathy    =    0.     !  value of the bathymetry. if (=0) bottom flat at jpkm1
-   nn_closea   =    0      !  remove (=0) or keep (=1) closed seas and lakes (ORCA)
-   nn_msh      =    1      !  create (=1) a mesh file or not (=0)
+   nn_msh      =    0      !  create (=1) a mesh file or not (=0)
    rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0)
    rn_isfhmin  =    1.00   !  treshold (m) to discriminate grounding ice to floating ice
@@ -93,6 +92,4 @@
    jpiglo      =     10    !  1st dimension of global domain --> i =jpidta
    jpjglo      =     12    !  2nd    -                  -    --> j =jpjdta
-   jpizoom     =      1    !  left bottom (i,j) indices of the zoom
-   jpjzoom     =      1    !  in data domain indices
    jperio      =      0    !  lateral cond. type (between 0 and 6)
                                  !  = 0 closed                 ;   = 1 cyclic East-West
@@ -103,4 +100,5 @@
    ln_use_jattr = .false.  !  use (T) the file attribute: open_ocean_jstart, if present
                            !  in netcdf input files, as the start j-row for reading
+   ln_domclo = .true.      ! computation of closed sea masks (see namclo)
 /
 !-----------------------------------------------------------------------
@@ -137,4 +135,24 @@
                         !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
    rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1)
+/
+!-----------------------------------------------------------------------
+&namclo ! (closed sea : need ln_domclo = .true. in namcfg)
+!-----------------------------------------------------------------------
+   rn_lon_opnsea = -2.0     ! longitude seed of open ocean
+   rn_lat_opnsea = -2.0     ! latitude  seed of open ocean
+   nn_closea = 8           ! number of closed seas ( = 0; only the open_sea mask will be computed)
+   !                name   ! lon_src ! lat_src ! lon_trg ! lat_trg ! river mouth area   ! net evap/precip correction scheme ! radius tgt   ! id trg
+   !                       ! (degree)! (degree)! (degree)! (degree)! local/coast/global ! (glo/rnf/emp)                     !     (m)      !
+   ! North American lakes
+   sn_lake(1) = 'superior' ,  -86.57 ,  47.30  , -66.49  , 50.45   , 'local'            , 'rnf'                             ,   550000.0 , 2    
+   sn_lake(2) = 'michigan' ,  -87.06 ,  42.74  , -66.49  , 50.45   , 'local'            , 'rnf'                             ,   550000.0 , 2    
+   sn_lake(3) = 'huron'    ,  -82.51 ,  44.74  , -66.49  , 50.45   , 'local'            , 'rnf'                             ,   550000.0 , 2    
+   sn_lake(4) = 'erie'     ,  -81.13 ,  42.25  , -66.49  , 50.45   , 'local'            , 'rnf'                             ,   550000.0 , 2    
+   sn_lake(5) = 'ontario'  ,  -77.72 ,  43.62  , -66.49  , 50.45   , 'local'            , 'rnf'                             ,   550000.0 , 2    
+   ! African Lake
+   sn_lake(6) = 'victoria' ,   32.93 ,  -1.08  ,  30.44  , 31.37   , 'coast'            , 'emp'                             ,   100000.0 , 3    
+   ! Asian Lakes
+   sn_lake(7) = 'caspian'  ,   50.0  ,  44.0   ,   0.0   ,  0.0    , 'global'           , 'glo'                             ,        0.0 , 1     
+   sn_lake(8) = 'aral'     ,   60.0  ,  45.0   ,   0.0   ,  0.0    , 'global'           , 'glo'                             ,        0.0 , 1    
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2019/ENHANCE-03_domcfg/src/dom_oce.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-03_domcfg/src/dom_oce.F90	(revision 11200)
+++ /NEMO/branches/2019/ENHANCE-03_domcfg/src/dom_oce.F90	(revision 11201)
@@ -36,5 +36,4 @@
    REAL(wp), PUBLIC ::   rn_e3zps_rat    !: minimum thickness ration for partial steps
    INTEGER , PUBLIC ::   nn_msh          !: = 1 create a mesh-mask file
-   INTEGER , PUBLIC ::   nn_closea       !: =0 suppress closed sea/lake from the ORCA domain or not (=1)
 
    INTEGER, PUBLIC :: nn_interp
@@ -50,5 +49,5 @@
    LOGICAL, PUBLIC ::   lzoom_n    =  .FALSE.   !: North zoom type flag
 
-   LOGICAL, PUBLIC ::   ln_closea  =  .FALSE.
+   LOGICAL, PUBLIC ::   ln_domclo  =  .FALSE.
 
    INTEGER       ::   jphgr_msh          !: type of horizontal mesh
@@ -210,12 +209,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdepw_0 !: w- depth              [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gde3w_0 !: w- depth (sum of e3w) [m]
-   ! 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0  !: t-depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  !: u-depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  !: v-depth              [m]
-
+   !
    INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1)
    INTEGER, PUBLIC ::   nlb10              !: shallowest W level Bellow ~10m (nla10 + 1) 
-
+   !
    !! 1D reference  vertical coordinate
    !! =-----------------====------
@@ -243,5 +238,9 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask             !: surface mask at T-,U-, V- and F-pts
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask                        !: land/ocean mask at W- pts               
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: msk_opnsea, msk_closea                 !: open ocean mask, closed sea mask (all of them)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: msk_glo  , msk_rnf  , msk_emp                !: closed sea masks
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: msk_gloid, msk_rnfid, msk_empid              !: closed sea masks
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4)
@@ -302,5 +301,5 @@
    INTEGER FUNCTION dom_oce_alloc()
       !!----------------------------------------------------------------------
-      INTEGER, DIMENSION(12) :: ierr
+      INTEGER, DIMENSION(11) :: ierr
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -328,26 +327,25 @@
          &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , STAT=ierr(5) )                       
          !
-      ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) , STAT=ierr(6)  )
-         !
-         !
-      ALLOCATE( gdept_1d(jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) ,gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(7) )
+      ALLOCATE( gdept_1d(jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) ,gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(6) )
          !
       ALLOCATE( bathy(jpi,jpj),mbathy(jpi,jpj), tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        & 
          &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) ,     &
-         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(9) )
+         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(7) )
          !
       ALLOCATE( misfdep(jpi,jpj) , mikt(jpi,jpj) , miku(jpi,jpj) ,     &
-         &      risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(10) )
+         &      risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(8) )
          !
       ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     & 
-         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , STAT=ierr(11) )
-         !
-      ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) )
-
+         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , wmask(jpi,jpj,jpk) , STAT=ierr(9) )
+         !
       ALLOCATE( hbatv (jpi,jpj) , hbatf (jpi,jpj) ,     &
          &      hbatt (jpi,jpj) , hbatu (jpi,jpj) ,     &
          &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     &
          &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     &
-         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , STAT=ierr(8) )
+         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , STAT=ierr(10) )
+
+      ALLOCATE( msk_opnsea(jpi,jpj), msk_closea(jpi,jpj),                   &
+         &      msk_glo  (jpi,jpj), msk_rnf  (jpi,jpj), msk_emp  (jpi,jpj), &
+         &      msk_gloid(jpi,jpj), msk_rnfid(jpi,jpj), msk_empid(jpi,jpj), STAT=ierr(11) )
       !
       dom_oce_alloc = MAXVAL(ierr)
Index: /NEMO/branches/2019/ENHANCE-03_domcfg/src/domain.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-03_domcfg/src/domain.F90	(revision 11200)
+++ /NEMO/branches/2019/ENHANCE-03_domcfg/src/domain.F90	(revision 11201)
@@ -23,9 +23,8 @@
    USE dom_oce         ! domain: ocean
    USE phycst          ! physical constants
-   !  USE closea          ! closed seas
    USE domhgr          ! domain: set the horizontal mesh
    USE domzgr          ! domain: set the vertical mesh
    USE dommsk          ! domain: set the mask system
-   USE domwri          ! domain: write the meshmask file
+   USE domclo          ! domain: set closed sea mask
    !
    USE in_out_manager  ! I/O manager
@@ -75,26 +74,15 @@
       !                       !==  Reference coordinate system  ==!
       !
-      CALL dom_nam               ! read namelist ( namrun, namdom )
-      !
-      !   CALL dom_clo               ! Closed seas and lake
+      CALL dom_nam                  ! read namelist ( namrun, namdom )
+      !
+      CALL dom_hgr                  ! Horizontal mesh
+      !
+      CALL dom_zgr                  ! Vertical mesh and bathymetry
+      !
+      IF ( ln_domclo .OR. nmsh > 0 ) CALL dom_msk                  ! compute mask (needed by dom_clo
+      !
+      IF ( ln_domclo ) CALL dom_clo ! Closed seas and lake
       ! 
-      CALL dom_hgr               ! Horizontal mesh
-      !
-      CALL dom_zgr               ! Vertical mesh and bathymetry
-      !
-      CALL dom_msk               ! Masks
-      !
-      ht_0(:,:) = e3t_0(:,:,1) * tmask(:,:,1)   ! Reference ocean thickness
-      hu_0(:,:) = e3u_0(:,:,1) * umask(:,:,1)
-      hv_0(:,:) = e3v_0(:,:,1) * vmask(:,:,1)
-      DO jk = 2, 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
-      !
-      CALL cfg_write             ! create the configuration file
-      !
-      CALL dom_wri
+      CALL cfg_write                ! create the configuration file
       !
    END SUBROUTINE dom_init
@@ -118,5 +106,5 @@
       NAMELIST/namdom/ nn_bathy, cn_topo, cn_bath, cn_lon, cn_lat, nn_interp,                        &
          &             rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin,           &
-         &             rn_atfp , rn_rdt   , nn_closea   , ln_crs      , jphgr_msh ,                  &
+         &             rn_atfp , rn_rdt   , ln_crs      , jphgr_msh ,                  &
          &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         &
          &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  &
@@ -142,30 +130,9 @@
          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,*) '      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,*) '      initial time of day in hhmm     nn_time0   = ', nn_time0
-         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,*) '      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
-         WRITE(numout,*) '      IS coupling at the restart step ln_iscpl   = ', ln_iscpl
       ENDIF
 
@@ -245,8 +212,4 @@
          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,*) '      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
@@ -453,10 +416,20 @@
       END DO
       CALL iom_rstput( 0, 0, inum, 'bathy_metry'   , z2d , ktype = jp_r4 )
-
-      !
-      IF( ln_sco ) THEN             ! s-coordinate: store grid stiffness ratio  (Not required anyway)
-         CALL dom_stiff( z2d )
-         CALL iom_rstput( 0, 0, inum, 'stiffness', z2d )        !    ! Max. grid stiffness ratio
-      ENDIF
+      !
+      !                              !== closed sea ==!
+      IF (ln_domclo) THEN
+         ! mask for the open sea
+         CALL iom_rstput( 0, 0, inum, 'mask_opensea', msk_opnsea, ktype = jp_i4 )
+         ! mask for all the under closed sea
+         CALL iom_rstput( 0, 0, inum, 'mask_csundef', msk_closea, ktype = jp_i4 )
+         ! mask for global, local net precip, local net precip and evaporation correction
+         CALL iom_rstput( 0, 0, inum, 'mask_csglo', msk_glo, ktype = jp_i4 )
+         CALL iom_rstput( 0, 0, inum, 'mask_csemp', msk_emp, ktype = jp_i4 )
+         CALL iom_rstput( 0, 0, inum, 'mask_csrnf', msk_rnf, ktype = jp_i4 )
+         ! mask for the various river mouth (in case multiple lake in the same outlet)
+         CALL iom_rstput( 0, 0, inum, 'mask_csgrpglo', msk_gloid, ktype = jp_i4 )
+         CALL iom_rstput( 0, 0, inum, 'mask_csgrpemp', msk_empid, ktype = jp_i4 )
+         CALL iom_rstput( 0, 0, inum, 'mask_csgrprnf', msk_rnfid, ktype = jp_i4 )
+      END IF
       !
       !                                ! ============================
Index: /NEMO/branches/2019/ENHANCE-03_domcfg/src/dombat.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-03_domcfg/src/dombat.F90	(revision 11200)
+++ /NEMO/branches/2019/ENHANCE-03_domcfg/src/dombat.F90	(revision 11201)
@@ -8,5 +8,4 @@
    USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp           ! distributed memory computing library
-   USE wrk_nemo          ! Memory allocation
    USE agrif_modutil
    USE bilinear_interp
Index: /NEMO/branches/2019/ENHANCE-03_domcfg/src/domclo.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-03_domcfg/src/domclo.F90	(revision 11201)
+++ /NEMO/branches/2019/ENHANCE-03_domcfg/src/domclo.F90	(revision 11201)
@@ -0,0 +1,312 @@
+MODULE domclo
+   !!======================================================================
+   !!                       ***  MODULE  domclo  ***
+   !! domclo : definition of closed sea mask needed by NEMO
+   !!======================================================================
+   !! History : NEMO 4.0 ! 07-2019 (P. Mathiot) original version
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   dom_clo    : definition of closed sea mask needed by NEMO
+   !!----------------------------------------------------------------------
+   USE dom_oce         ! ocean space and time domain
+   USE domngb          ! closest point algorithm
+   USE phycst          ! rpi, rad, ra
+   USE domutil         ! flood filling algorithm (fill_pool)
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! lateral boundary condition - MPP exchanges
+   USE lib_mpp
+   USE lib_fortran
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC dom_clo      ! routine called by domain module
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: closea.F90 7213 2016-11-09 09:07:42Z timgraham $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dom_clo
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_clo  ***
+      !!        
+      !! ** Purpose :   Closed sea mask definition
+      !!
+      !! ** Method  :   A flood filling algorithm is used to detect lake and open ocean
+      !!                Namelist provide seed for each lake and condition on how to
+      !!                manage net evap/runoff.
+      !!
+      !! ** Action  :  - compute open ocean mask 
+      !!               - loop over all the lake in namelist 
+      !!               - define lake msk (msk_csglo/emp/rnf)
+      !!               - define river mouth mask for each lake or group of lake with the same river mouth
+      !!                 (msk_csgrpglo/emp/rnf) 
+      !!               
+      !!----------------------------------------------------------------------
+
+      TYPE closea 
+         CHARACTER(256) :: cname                       ! name
+         REAL(wp)       :: rlonsrc                     ! seed src location (lon)
+         REAL(wp)       :: rlatsrc                     ! seed src location (lat)
+         REAL(wp)       :: rlontrg                     ! seed trg location (lon)
+         REAL(wp)       :: rlattrg                     ! seed trg location (lat)
+         CHARACTER(256) :: cloctrg                     ! where water is spread
+         CHARACTER(256) :: cschtrg                     ! how   water is spread
+         REAL(wp)       :: radtrg                      ! radius of closed sea river mouth
+         INTEGER        :: idtrg                       ! target id in case multiple lakes for the same river mouth
+      END TYPE
+         
+      INTEGER, PARAMETER :: jp_closea = 98 ! number maximal of closea
+      INTEGER :: ji, jj, jcs, jsch  ! loop indexes
+      INTEGER :: jglo, jemp, jrnf   ! closed sea indexes for global, emp, rnf spreading
+      INTEGER :: jiseed, jjseed     ! seed indexes
+      INTEGER :: jglo, jloc, jcoast ! counter
+      INTEGER :: ios
+      INTEGER :: nn_closea          ! number of closed seas
+
+      REAL(wp) :: zdistseed         ! distance to seed
+      REAL(wp) :: zarea             ! river mouth area
+      REAL(wp) :: rn_lon_opnsea, rn_lat_opnsea ! open sea seed
+      REAL(wp), DIMENSION(1)       :: zchk, zradtrg
+      REAL(wp), DIMENSION(jpi,jpj) :: zmsksrc, zmsktrg, zmsk_coastline, zdist
+
+      CHARACTER(256) :: csch, cloc        ! scheme name for water spreading (glo, rnf, emp)
+      TYPE(closea)  , DIMENSION(jp_closea)   :: sn_lake   ! lake properties
+
+      LOGICAL :: lskip     ! flag in case lake seed on land or already filled (...)
+
+      !!----------------------------------------------------------------------
+      !! 0 : Read namelist for closed sea definition
+      !!----------------------------------------------------------------------
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*)'dom_clo : closed seas '
+      IF(lwp) WRITE(numout,*)'~~~~~~~'
+
+      NAMELIST/namclo/ rn_lon_opnsea, rn_lat_opnsea, nn_closea, sn_lake
+      !!---------------------------------------------------------------------
+      PRINT *, rn_lon_opnsea, rn_lat_opnsea, nn_closea
+      !!
+      REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
+      READ  ( numnam_ref, namclo, IOSTAT = ios, ERR = 901 )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namclo in reference namelist', lwp )
+      REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
+      READ  ( numnam_cfg, namclo, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namclo in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namclo )
+
+      !!----------------------------------------------------------------------
+      !! 1 : Define open ocean and closed sea masks
+      !!----------------------------------------------------------------------
+
+      !! 1.1 get closest point to the lat/lon seed
+      msk_opnsea(:,:) = ssmask(:,:)
+      CALL dom_ngb(rn_lon_opnsea, rn_lat_opnsea, jiseed, jjseed, zdistseed, 'T')
+
+      !! 1.2 fill select cell to -1
+      CALL fill_pool( jiseed, jjseed, msk_opnsea, -1._wp )
+
+      !! sanity check on the msk value (closea mask should be 99 otherwise, lake already processed)
+      !! check if some lake are not connected
+      zchk = 0._wp
+      IF (mi0(jiseed) == mi1(jiseed) .AND. mj0(jjseed) == mj1(jjseed)) zchk = ssmask(mi0(jiseed),mj0(jjseed))
+      CALL mpp_max('domclo',zchk)
+      IF (zchk(1) == 0._wp) CALL ctl_stop( 'STOP', 'open sea seed is on land, please update namelist (rn_lon_opnsea,rn_lat_opnsea)' ) 
+
+      !! 1.3 set to 0 everything >0 and revert mask
+      IF (lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)'open ocean: '
+         WRITE(numout,'(a,2f7.2)')'    lon/lat seed to detect main ocean is: ', rn_lon_opnsea, rn_lat_opnsea
+         WRITE(numout,'(a,2i7)'  )'    i/j     seed to detect main ocean is: ', jiseed, jjseed
+         WRITE(numout,'(a,f8.0)' )'    distance (m) between namelist seed location and model seed location: ', zdistseed
+         WRITE(numout,*)
+      END IF
+      
+      WHERE (msk_opnsea(:,:) > 0) msk_opnsea(:,:) = 0 ! mask all closed seas
+      WHERE (msk_opnsea(:,:) < 0) msk_opnsea(:,:) = 1 ! restore mask value
+
+      !! 1.4 Define closed sea mask (all of them, ie defined in the namelist or not)
+      !! needed to remove the undefined closed seas at the end
+      msk_closea = (ssmask - msk_opnsea) * 99
+
+      !!----------------------------------------------------------------------
+      !! 2 : compute closed sea mask for global, rnf and emp cases
+      !!----------------------------------------------------------------------
+
+      !! 2.1 : initialisation of masks and lake indexes
+      jglo = 1 ; jrnf = 1 ; jemp = 1
+      !! mask used to group lake by net evap/precip distribution technics
+      msk_glo = msk_closea
+      msk_rnf = msk_closea
+      msk_emp = msk_closea
+
+      !! mask used to group multiple lake with the same river mouth (great lake for example)
+      msk_gloid = 0.0_wp
+      msk_rnfid = 0.0_wp
+      msk_empid = 0.0_wp
+
+      IF (lwp) WRITE(numout,*)'closed seas: '
+ 
+      DO jcs = 1,nn_closea
+
+         !! 2.2 : print out and sanity check
+         !! flag changed in sanity checks
+         lskip = .FALSE.
+
+         !! define seed to extract the closed sea jcs
+         CALL dom_ngb(sn_lake(jcs)%rlonsrc, sn_lake(jcs)%rlatsrc, jiseed, jjseed, zdistseed, 'T')
+
+         !! how the excess of the closed seas is spread out:
+         IF (lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,'(a,a10,a,2f7.2,a)')'    processing lake ',TRIM(sn_lake(jcs)%cname)             &
+                 &        ,' ( lat/lon : ',sn_lake(jcs)%rlatsrc, sn_lake(jcs)%rlonsrc,' )'
+         END IF
+         cloc = sn_lake(jcs)%cloctrg
+         csch = sn_lake(jcs)%cschtrg
+  
+         !! set up indexes and mask  ! SELECT CASE ...
+         SELECT CASE (csch)
+         CASE('glo')
+            jsch = jglo
+            zmsksrc = msk_glo
+            zmsktrg = msk_gloid
+            IF (lwp) WRITE(numout,*)'        net evap/precip will be spread globally'
+         CASE('rnf') 
+            jsch = jrnf
+            zmsksrc = msk_rnf
+            zmsktrg = msk_rnfid
+            IF (lwp) WRITE(numout,*)'        net precip will be spread locally and net evap globally'
+         CASE('emp')
+            jsch = jemp 
+            zmsksrc = msk_emp
+            zmsktrg = msk_empid
+            IF (lwp) WRITE(numout,*)'        net precip will be spread locally'
+         CASE DEFAULT
+            CALL ctl_stop( 'STOP', 'domclo: unknown target type for lake (should be glo, emp, rnf)' )
+         END SELECT
+
+         !! sanity check on the msk value (closea mask should be 99 otherwise, lake already processed)
+         !! check if some lake are not connected
+         zchk = 0._wp
+         IF (mi0(jiseed) == mi1(jiseed) .AND. mj0(jjseed) == mj1(jjseed)) zchk = msk_closea(mi0(jiseed),mj0(jjseed))
+         CALL mpp_max('domclo',zchk)
+         
+         IF (zchk(1) /= 99._wp) THEN
+            IF (lwp) THEN
+               WRITE(numout,*)
+               WRITE(numout,*)'W A R N I N G: lake ',TRIM(sn_lake(jcs)%cname),' already processed or seed on land, we skip it.'
+               WRITE(numout,*)'It is possible this lake is connected to another one already processed.'
+               WRITE(numout,*)'It is possible the seed is on land because the lake is not present in the bathymetry file'
+               WRITE(numout,*)'If concerned about it, check input file and namelist definition of this lake.'
+            END IF
+            lskip = .TRUE.
+         END IF
+
+         !! 2.3 : compute mask for global, rnf, emp case
+         IF (.NOT. lskip) THEN
+            !! fill close sea mask with counter value
+            CALL fill_pool( jiseed, jjseed, zmsksrc, REAL(jsch  ,8))
+            WHERE (zmsksrc == REAL(jsch,8))
+               msk_closea = 0._wp
+               zmsktrg    = sn_lake(jcs)%idtrg
+            END WHERE
+
+            !! compute location of river mouth and distance to river mouth
+            IF (cloc /= 'global') THEN
+               !! set a minimum value for radius of the river influence
+               zradtrg = 0._wp
+               IF (mi0(jiseed) == mi1(jiseed) .AND. mj0(jjseed) == mj1(jjseed)) THEN
+                  zradtrg(1) = MAX( sn_lake(jcs)%radtrg, e1t( mi0(jiseed), mj0(jjseed) ), e2t( mi0(jiseed), mj0(jjseed) ) )
+               END IF
+               CALL mpp_max('domclo',zradtrg)
+              
+               !! compute seed location for print 
+               CALL dom_ngb(sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg, jiseed, jjseed, zdistseed, 'T')
+               !! compute distance to river mouth
+               zdist(:,:) = dist(sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg, glamt, gphit)
+            END IF
+
+            !! define estuary
+            !! deal with local/coast cases
+            SELECT CASE (cloc)
+            CASE ('global')
+               WHERE (msk_opnsea(:,:) == 1) zmsktrg = sn_lake(jcs)%idtrg
+
+            CASE ('local')
+               !! compute mask
+               WHERE (zdist(:,:) < zradtrg(1) .AND. msk_opnsea(:,:) == 1 ) zmsktrg = sn_lake(jcs)%idtrg
+               !! print
+               IF (lwp) WRITE(numout,'(a,f7.0,a,2f7.2,a,2i7,a)')'         river mouth area is defined by         points within ',zradtrg(1)         &
+                             &        ,' m of lat/lon ', sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg                                  &
+                             &        ,' (closest point is i/j ',jiseed, jjseed,')'
+  
+            CASE ('coast')
+               !! define coastline mask
+               zmsk_coastline = 0._wp
+               DO jj=2,jpj-1
+                  DO ji=2,jpi-1
+                     IF ( ssmask(ji,jj) == 1._wp .AND. SUM(ssmask(ji-1:ji+1,jj-1:jj+1)) < 9 ) zmsk_coastline(ji,jj) = 1._wp
+                  END DO
+               END DO
+               CALL lbc_lnk('domclo', zmsk_coastline,'T',1._wp)
+   
+               !! compute mask
+               WHERE ( zdist(:,:) < zradtrg(1) .AND. zmsk_coastline(:,:) == 1 .AND. msk_opnsea(:,:) == 1 ) zmsktrg = sn_lake(jcs)%idtrg
+               !! print
+               IF (lwp) WRITE(numout,'(a,f7.0,a,2f7.2,a,2i7,a)')'         river mouth area is defined by coastal points within ',zradtrg(1) &
+                             &        ,' m of lat/lon ', sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg                                  &
+                             &        ,' (closest point is i/j ',jiseed, jjseed,')'
+
+            CASE DEFAULT
+               CALL ctl_stop( 'STOP', 'domclo: unknown target type for lake (should be global, coast or local)' )
+
+            END SELECT
+   
+            !! sanity check
+            zarea = glob_sum('domclo',zmsktrg * msk_opnsea)
+            IF (zarea == 0._wp) CALL ctl_stop( 'STOP', 'river mouth area is 0, tune lon/lat trg or radtrg for this lake') 
+
+            !! set up indexes and mask 
+            SELECT CASE (csch)
+            CASE ('glo')
+               jglo = jglo + 1
+               msk_glo   = zmsksrc
+               msk_gloid = zmsktrg
+            CASE ('rnf')
+               jrnf = jrnf + 1
+               msk_rnf   = zmsksrc
+               msk_rnfid = zmsktrg
+            CASE ('emp')
+               jemp = jemp + 1
+               msk_emp   = zmsksrc
+               msk_empid = zmsktrg
+            END SELECT
+
+         END IF ! lskip
+
+      END DO ! nn_closea
+
+      !!----------------------------------------------------------------------
+      !! 3 : clean the masks of possible remaining undefined closed seas
+      !!----------------------------------------------------------------------
+
+      !! mask all the cells not defined as closed sea
+      WHERE ( msk_glo == 99 ) msk_glo = 0
+      WHERE ( msk_rnf == 99 ) msk_rnf = 0
+      WHERE ( msk_emp == 99 ) msk_emp = 0
+
+      !!  non defined closed sea
+      WHERE (msk_closea > 0) msk_closea = 1
+
+   END SUBROUTINE dom_clo
+
+   !!======================================================================
+
+END MODULE domclo
+
Index: /NEMO/branches/2019/ENHANCE-03_domcfg/src/domhgr.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-03_domcfg/src/domhgr.F90	(revision 11200)
+++ /NEMO/branches/2019/ENHANCE-03_domcfg/src/domhgr.F90	(revision 11201)
@@ -28,5 +28,4 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
-   USE timing         ! Timing
 
    IMPLICIT NONE
Index: /NEMO/branches/2019/ENHANCE-03_domcfg/src/dommsk.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-03_domcfg/src/dommsk.F90	(revision 11200)
+++ /NEMO/branches/2019/ENHANCE-03_domcfg/src/dommsk.F90	(revision 11201)
@@ -24,4 +24,5 @@
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
+   USE domwri         ! domain: write the meshmask file
    USE bdy_oce        ! open boundary
    !
@@ -60,5 +61,5 @@
       !!      and ko_bot, the indices of the fist and last ocean t-levels which 
       !!      are either defined in usrdef_zgr or read in zgr_read.
-      !!                The velocity masks (umask, vmask, wmask, wumask, wvmask) 
+      !!                The velocity masks (umask, vmask) 
       !!      are deduced from a product of the two neighboring tmask.
       !!                The vorticity mask (fmask) is deduced from tmask taking
@@ -75,5 +76,5 @@
       !!                due to cyclic or North Fold boundaries as well as MPP halos.
       !!
-      !! ** Action :   tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask 
+      !! ** Action :   tmask, umask, vmask, wmask : land/ocean mask 
       !!                         at t-, u-, v- w, wu-, and wv-points (=0. or 1.)
       !!               fmask   : land/ocean mask at f-point (=0., or =1., or 
@@ -188,10 +189,6 @@
       !-----------------------------------------
       wmask (:,:,1) = tmask(:,:,1)     ! surface
-      wumask(:,:,1) = umask(:,:,1)
-      wvmask(:,:,1) = vmask(:,:,1)
       DO jk = 2, jpk                   ! interior values
          wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
-         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
-         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
       END DO
 
@@ -202,5 +199,4 @@
       ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 )
       ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 )
-
 
       ! Interior domain mask  (used for global sum)
@@ -289,4 +285,7 @@
       ENDIF
       !
+      ! write out mesh mask
+      IF ( nn_msh > 0 ) CALL dom_wri
+      !
    END SUBROUTINE dom_msk
    
Index: /NEMO/branches/2019/ENHANCE-03_domcfg/src/domngb.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-03_domcfg/src/domngb.F90	(revision 11200)
+++ /NEMO/branches/2019/ENHANCE-03_domcfg/src/domngb.F90	(revision 11201)
@@ -11,4 +11,5 @@
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
+   USE phycst
    !
    USE in_out_manager ! I/O manager
@@ -18,5 +19,6 @@
    PRIVATE
 
-   PUBLIC   dom_ngb   ! routine called in iom.F90 module
+   PUBLIC   dom_ngb   ! routine called in iom.F90 and domclo.F90 module
+   PUBLIC   dist
 
    !!----------------------------------------------------------------------
@@ -27,5 +29,5 @@
 CONTAINS
 
-   SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, kkk )
+   SUBROUTINE dom_ngb( plon, plat, kii, kjj, rdist, cdgrid, kkk )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dom_ngb  ***
@@ -37,4 +39,5 @@
       !!----------------------------------------------------------------------
       REAL(wp)        , INTENT(in   ) ::   plon, plat   ! longitude,latitude of the point
+      REAL(wp)        , INTENT(  out) ::   rdist        ! distance between the located point and the source
       INTEGER         , INTENT(  out) ::   kii, kjj     ! i-,j-index of the closes grid point
       INTEGER         , INTENT(in   ), OPTIONAL :: kkk  ! k-index of the mask level used
@@ -43,5 +46,4 @@
       INTEGER :: ik         ! working level
       INTEGER , DIMENSION(2) ::   iloc
-      REAL(wp)               ::   zlon, zmini
       REAL(wp), DIMENSION(jpi,jpj) ::   zglam, zgphi, zmask, zdist
       !!--------------------------------------------------------------------
@@ -57,18 +59,12 @@
       END SELECT
 
-      zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360
-      zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360
-      IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270
-      IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180
-      zglam(:,:) = zglam(:,:) - zlon
-
-      zgphi(:,:) = zgphi(:,:) - plat
-      zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:)
+      zdist = dist(plon, plat, zglam, zgphi)
       
       IF( lk_mpp ) THEN  
-         CALL mpp_minloc( 'domngb', zdist(:,:), zmask, zmini, iloc)
+         CALL mpp_minloc( 'domngb', zdist(:,:), zmask, rdist, iloc)
          kii = iloc(1) ; kjj = iloc(2)
       ELSE
          iloc(:) = MINLOC( zdist(:,:), mask = zmask(:,:) == 1.e0 )
+         rdist = MINVAL( zdist(:,:) )
          kii = iloc(1) + nimpp - 1
          kjj = iloc(2) + njmpp - 1
@@ -77,4 +73,29 @@
    END SUBROUTINE dom_ngb
 
+   FUNCTION dist(plonsrc, platsrc, plontrg, plattrg)
+      REAL(wp), INTENT(in) :: plonsrc, platsrc                     ! lat/lon of the source point
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: plontrg, plattrg ! lat/lon of the target (2d in this case)
+
+      REAL(wp) :: zxs, zys, zzs
+      REAL(wp), DIMENSION(jpi,jpj) :: zxt, zyt, zzt
+
+      REAL(wp), DIMENSION(jpi,jpj) :: dist ! distance between src and trg
+
+      zxs = COS( rad * platsrc ) * COS( rad * plonsrc )
+      zys = COS( rad * platsrc ) * SIN( rad * plonsrc )
+      zzs = SIN( rad * platsrc )
+
+      zxt = COS( rad * plattrg ) * COS( rad * plontrg )
+      zyt = COS( rad * plattrg ) * SIN( rad * plontrg )
+      zzt = SIN( rad * plattrg )
+
+      dist(:,:) = ( zxs - zxt(:,:) )**2   &
+         &      + ( zys - zyt(:,:) )**2   &
+         &      + ( zzs - zzt(:,:) )**2
+
+      dist = ra * SQRT( dist )
+
+   END FUNCTION dist
+
    !!======================================================================
 END MODULE domngb
Index: /NEMO/branches/2019/ENHANCE-03_domcfg/src/domutl.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-03_domcfg/src/domutl.F90	(revision 11201)
+++ /NEMO/branches/2019/ENHANCE-03_domcfg/src/domutl.F90	(revision 11201)
@@ -0,0 +1,271 @@
+MODULE domutil
+   !!
+   !!======================================================================
+   !!                       ***  MODULE  closea  ***
+   !! modutl : specific module to do
+   !!              - 2d flood filling (closea)
+   !!              - 3d flood filling (zgr_isf)
+   !!
+   !!======================================================================
+   !! History :   4.0  !  03-18  (P. Mathiot) Original code
+   !!----------------------------------------------------------------------
+   !!
+   !!----------------------------------------------------------------------
+   !!   fillpool    : replace all point connected to the seed by a fillvalue 
+   !!----------------------------------------------------------------------
+   !!
+   USE dom_oce         ! ocean space and time domain
+   USE domngb
+   USE in_out_manager    ! I/O manager
+   USE lib_mpp
+   USE lbclnk
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC fill_pool      ! routine called by domain module
+
+   INTERFACE fill_pool
+      MODULE PROCEDURE FillPool2D, FillPool3D
+   END INTERFACE
+
+   TYPE idx
+      INTEGER :: i,j,k
+      TYPE(idx), POINTER :: next
+   END TYPE idx
+
+CONTAINS
+   SUBROUTINE FillPool3D(kiseed, kjseed, kkseed, rdta, rfill)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE FillPool3D  ***
+      !!
+      !! ** Purpose :  Replace all area surrounding by mask value by kifill value
+      !!
+      !! ** Method  :  flood fill algorithm
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER,                          INTENT(in)    :: kiseed, kjseed, kkseed
+      REAL(wp),                         INTENT(in)    :: rfill    ! filling value
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: rdta     ! input data
+      REAL(wp), DIMENSION(jpi,jpj,jpk)                :: rseedmap, rseedmap_b  !
+ 
+      INTEGER :: ip=0                 ! size of the pile
+      INTEGER :: ii  , ij  , ik   , kii, kjj, jj, kk    ! working integer
+      INTEGER :: iip1, ijp1, ikp1
+      INTEGER :: iim1, ijm1, ikm1
+      INTEGER :: nseed
+      TYPE (idx), POINTER :: seed
+      !!---------------------------------------------------------------------- 
+      !
+      ! initialisation seed
+      NULLIFY(seed)
+      !
+      ! create the first seed and update nseed for all processors
+      nseed = 0
+      rseedmap = 0.
+      ii=kiseed; ij=kjseed ; ik=kkseed
+      IF ((mi0(ii) == mi1(ii)) .AND. (mj0(ij) == mj1(ij))) THEN   ! T if seed on the local domain
+         ii = mi0(ii) ; ij = mj0(ij);                             ! conversion to local index
+         IF (rdta(ii,ij,ik) > 0 ) THEN                            ! create seed if not on land
+            CALL add_head_idx(seed,ii,ij,ik)
+            rdta    (ii,ij,ik) = rfill
+            rseedmap(ii,ij,ik) = 1.
+         END IF
+      END IF
+      nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! nseed =0 means on land => WARNING later on
+      !
+      ! loop until the pile size is 0 or if the pool is larger than the critical size
+      IF (nseed > 0) THEN
+         ! seed on ocean continue
+         DO WHILE ( nseed /= 0 )
+            DO WHILE ( ASSOCIATED(seed) )
+               ip=ip+1
+               ii=seed%i; ij=seed%j ; ik=seed%k ; rseedmap(ii,ij,ik)=1.
+               ! 
+               ! update bathy and update pile size
+               CALL del_head_idx(seed) 
+               !
+               ! check neighbour cells
+               iip1=ii+1 ; IF (iip1 == jpi+1 ) iip1=ii
+               iim1=ii-1 ; IF (iim1 == 0     ) iim1=ii
+               ijp1=ij+1 ; IF (ijp1 == jpj+1 ) ijp1=ij
+               ijm1=ij-1 ; IF (ijm1 == 0     ) ijm1=ij
+               ikp1=ik+1 ; IF (ikp1 == jpk+1 ) ikp1=ik
+               ikm1=ik-1 ; IF (ikm1 == 0     ) ikm1=ik
+               !
+               ! ji direction
+               IF (rdta(ii,ijp1,ik) > 0 .AND. rdta(ii,ijp1,ik) /= rfill ) THEN ; CALL add_head_idx(seed,ii,ijp1,ik) ; rdta(ii,ijp1,ik) = rfill ; rseedmap(ii,ijp1,ik)=1. ; END IF
+               IF (rdta(ii,ijm1,ik) > 0 .AND. rdta(ii,ijm1,ik) /= rfill ) THEN ; CALL add_head_idx(seed,ii,ijm1,ik) ; rdta(ii,ijm1,ik) = rfill ; rseedmap(ii,ijm1,ik)=1. ; END IF
+               !
+               ! jj direction
+               IF (rdta(iip1,ij,ik) > 0 .AND. rdta(iip1,ij,ik) /= rfill ) THEN ; CALL add_head_idx(seed,iip1,ij,ik) ; rdta(iip1,ij,ik) = rfill ; rseedmap(iip1,ij,ik)=1. ; END IF
+               IF (rdta(iim1,ij,ik) > 0 .AND. rdta(iim1,ij,ik) /= rfill ) THEN ; CALL add_head_idx(seed,iim1,ij,ik) ; rdta(iim1,ij,ik) = rfill ; rseedmap(iim1,ij,ik)=1. ; END IF
+               !
+               ! jk direction
+               IF (rdta(ii,ij,ikp1) > 0 .AND. rdta(ii,ij,ikp1) /= rfill ) THEN ; CALL add_head_idx(seed,ii,ij,ikp1) ; rdta(ii,ij,ikp1) = rfill ; rseedmap(ii,ij,ik)=1. ; END IF
+               IF (rdta(ii,ij,ikm1) > 0 .AND. rdta(ii,ij,ikm1) /= rfill ) THEN ; CALL add_head_idx(seed,ii,ij,ikm1) ; rdta(ii,ij,ikm1) = rfill ; rseedmap(ii,ij,ik)=1. ; END IF
+            END DO
+            !
+            ! exchange seed
+            nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! this is the sum of all the point check this iteration
+            !
+            rseedmap_b(:,:,:)=rseedmap(:,:,:)
+            CALL lbc_lnk('domutil', rseedmap, 'T', 1.)
+            !
+            ! build new list of seed
+            DO ii=1,jpi
+               DO jj=1,jpj
+                  DO kk=1,jpk
+                     IF (rseedmap(ii,jj,kk) > 0.0 .AND. rseedmap(ii,jj,kk) /= rseedmap_b(ii,jj,kk) .AND. rdta(ii,jj,kk) > 0 .AND. rdta(ii,jj,kk) /= rfill) THEN
+                        CALL add_head_idx(seed,ii,jj,kk)
+                     END IF
+                  END DO
+               END DO
+            END DO
+            !
+            ! reset map of seed
+            rseedmap(:,:,:)=0.0
+            !
+         END DO
+      ELSE
+         WRITE(numout,*) 'W A R N I N G : SEED (',ii,ij,') is on land, nothing to do'
+      END IF
+ 
+   END SUBROUTINE FillPool3D
+
+   SUBROUTINE FillPool2D(kiseed, kjseed, rdta, rfill)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE FillPool3D  ***
+      !!
+      !! ** Purpose :  Replace all area surrounding by mask value by kifill value
+      !!
+      !! ** Method  :  flood fill algorithm
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER,                    INTENT(in)    :: kiseed, kjseed
+      REAL(wp),                   INTENT(in)    :: rfill    ! filling value
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: rdta     ! input data
+      REAL(wp), DIMENSION(jpi,jpj) :: rseedmap, rseedmap_b 
+ 
+      INTEGER :: ip=0                     ! size of the pile
+      INTEGER :: ii  , ij  , jj, kii, kjj     ! working integer
+      INTEGER :: iip1, ijp1               ! working integer
+      INTEGER :: iim1, ijm1
+      INTEGER :: nseed
+      TYPE (idx), POINTER :: seed
+      !!----------------------------------------------------------------------
+      !
+      ! initialisation seed
+      NULLIFY(seed)
+      !
+      ! create the first seed and update nseed for all processors
+      nseed = 0
+      rseedmap = 0.
+      ii = kiseed ; ij = kjseed ;
+      IF ((mi0(ii) == mi1(ii)) .AND. (mj0(ij) == mj1(ij))) THEN   ! T if seed on the local domain
+         ii = mi0(ii) ; ij = mj0(ij)                              ! conversion to local index
+         IF (rdta(ii,ij) > 0 ) THEN                               ! create seed if not on land
+            CALL add_head_idx(seed,ii,ij,1)
+            rdta    (ii,ij) = rfill
+            rseedmap(ii,ij) = 1.
+         END IF
+      END IF
+      nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! nseed =0 means on land => WARNING later on
+      !
+      ! loop until the pile size is 0 or if the pool is larger than the critical size
+      IF (nseed > 0) THEN
+         ! seed on ocean continue
+         DO WHILE ( nseed .NE. 0 )
+            DO WHILE ( ASSOCIATED(seed) )
+               ip=ip+1
+               ii=seed%i; ij=seed%j ; rseedmap(ii,ij)=1.
+               ! update pile size
+               CALL del_head_idx(seed) 
+               ! 
+               ! check neighbour cells
+               iip1=ii+1 ; IF (iip1 == jpi+1 ) iip1=ii
+               iim1=ii-1 ; IF (iim1 == 0     ) iim1=ii
+               ijp1=ij+1 ; IF (ijp1 == jpj+1 ) ijp1=ij
+               ijm1=ij-1 ; IF (ijm1 == 0     ) ijm1=ij
+               !
+               ! ji direction
+               IF (rdta(ii  ,ijp1) > 0 .AND. rdta(ii  ,ijp1) /= rfill ) THEN ; CALL add_head_idx(seed,ii  ,ijp1,1) ; rdta(ii  ,ijp1) = rfill ; rseedmap(ii  ,ijp1)=1. ; END IF
+               IF (rdta(ii  ,ijm1) > 0 .AND. rdta(ii  ,ijm1) /= rfill ) THEN ; CALL add_head_idx(seed,ii  ,ijm1,1) ; rdta(ii  ,ijm1) = rfill ; rseedmap(ii  ,ijm1)=1. ; END IF
+               !
+               ! jj direction
+               IF (rdta(iip1,ij  ) > 0 .AND. rdta(iip1,ij  ) /= rfill ) THEN ; CALL add_head_idx(seed,iip1,ij  ,1) ; rdta(iip1,ij  ) = rfill ; rseedmap(iip1,ij  )=1. ; END IF
+               IF (rdta(iim1,ij  ) > 0 .AND. rdta(iim1,ij  ) /= rfill ) THEN ; CALL add_head_idx(seed,iim1,ij  ,1) ; rdta(iim1,ij  ) = rfill ; rseedmap(iim1,ij  )=1. ; END IF
+            END DO
+            !
+            ! exchange seed
+            nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! this is the sum of all the point check this iteration
+            !
+            rseedmap_b(:,:)=rseedmap(:,:)
+            CALL lbc_lnk('domutil', rseedmap, 'T', 1.)
+            !
+            ! build new list of seed
+            DO ii=1,jpi
+               DO jj=1,jpj
+                  IF (rseedmap(ii,jj) > 0.0 .AND. rseedmap(ii,jj) /= rseedmap_b(ii,jj) .AND. rdta(ii,jj) > 0 .AND. rdta(ii,jj) /= rfill) THEN
+                     CALL add_head_idx(seed, ii, jj, 1)
+                  END IF
+               END DO
+            END DO 
+            !
+            ! reset map of seed
+            rseedmap(:,:)=0.0
+            !
+         END DO
+      ELSE
+         WRITE(numout,*) 'W A R N I N G : SEED (',ii,ij,') is on land, nothing to do', narea
+      END IF
+ 
+   END SUBROUTINE FillPool2D
+   !
+   !
+   ! subroutine to deals with link list
+   !
+   SUBROUTINE create_idx(pt_idx, ki, kj, kk)
+      TYPE (idx), POINTER :: pt_idx
+      INTEGER, INTENT(in) :: ki, kj, kk
+      !
+      ! initialised all field to NULL()
+      NULLIFY(pt_idx)
+      !
+      ! allocate new element
+      ALLOCATE(pt_idx)
+      pt_idx%i=ki ; pt_idx%j=kj ; pt_idx%k=kk ;
+      pt_idx%next => NULL()  
+   END SUBROUTINE create_idx
+
+   SUBROUTINE add_head_idx(pt_idx, ki, kj, kk)
+      TYPE (idx), POINTER :: pt_idx
+      TYPE (idx), POINTER :: zpt_new
+      INTEGER, INTENT(in) :: ki, kj, kk
+      !
+      ! allocate new element
+      ALLOCATE(zpt_new)
+      zpt_new%i=ki ; zpt_new%j=kj ; zpt_new%k=kk ;
+      zpt_new%next => NULL()
+      !
+      ! linked new element to main list
+      zpt_new%next => pt_idx
+      pt_idx => zpt_new
+   END SUBROUTINE add_head_idx
+
+   SUBROUTINE del_head_idx(pt_idx)
+      TYPE (idx), POINTER :: pt_idx
+      TYPE (idx), POINTER :: zpt_tmp
+      !
+      ! delete first element
+      IF (ASSOCIATED(pt_idx%next)) THEN
+         zpt_tmp => pt_idx%next
+         DEALLOCATE(pt_idx)
+         NULLIFY(pt_idx)
+         pt_idx => zpt_tmp
+      ELSE
+         NULLIFY(pt_idx)
+      ENDIF
+   END SUBROUTINE del_head_idx
+   !
+END MODULE domutil
Index: /NEMO/branches/2019/ENHANCE-03_domcfg/src/domwri.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-03_domcfg/src/domwri.F90	(revision 11200)
+++ /NEMO/branches/2019/ENHANCE-03_domcfg/src/domwri.F90	(revision 11201)
@@ -177,6 +177,4 @@
       ENDIF
       !
-   !   IF( ll_wd ) CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 )
-
       !                                     ! ============================
       CALL iom_close( inum )                !        close the files 
Index: /NEMO/branches/2019/ENHANCE-03_domcfg/src/domzgr.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-03_domcfg/src/domzgr.F90	(revision 11200)
+++ /NEMO/branches/2019/ENHANCE-03_domcfg/src/domzgr.F90	(revision 11201)
@@ -36,5 +36,4 @@
    !!---------------------------------------------------------------------
    USE dom_oce           ! ocean domain
-!   USE closea            ! closed seas
    !
    USE in_out_manager    ! I/O manager
@@ -164,5 +163,4 @@
       ! final adjustment of mbathy & check 
       ! -----------------------------------
-      IF( lzoom       )   CALL zgr_bat_zoom     ! correct mbathy in case of zoom subdomain
                           CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points
                           CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points
@@ -622,6 +620,4 @@
       ENDIF
       !
-    !  IF( nn_closea == 0 )   CALL clo_bat( bathy, mbathy )    !==  NO closed seas or lakes  ==!
-      !                       
       IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
          IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
@@ -636,57 +632,4 @@
       !
    END SUBROUTINE zgr_bat
-
-
-   SUBROUTINE zgr_bat_zoom
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE zgr_bat_zoom  ***
-      !!
-      !! ** Purpose : - Close zoom domain boundary if necessary
-      !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom
-      !!
-      !! ** Method  : 
-      !!
-      !! ** Action  : - update mbathy: level bathymetry (in level index)
-      !!----------------------------------------------------------------------
-      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain'
-      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
-      !
-      ! Zoom domain
-      ! ===========
-      !
-      ! Forced closed boundary if required
-      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) )      = 0
-      IF( lzoom_w )   mbathy(      mi0(jpizoom):mi1(jpizoom) , :  ) = 0
-      IF( lzoom_e )   mbathy(      mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0
-      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) )      = 0
-      !
-      ! Configuration specific domain modifications
-      ! (here, ORCA arctic configuration: suppress Med Sea)
-      IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
-         SELECT CASE ( jp_cfg )
-         !                                        ! =======================
-         CASE ( 2 )                               !  ORCA_R2 configuration
-            !                                     ! =======================
-            IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea'
-            ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices
-            ij0 =  98   ;   ij1 = 110
-            !                                     ! =======================
-         CASE ( 05 )                              !  ORCA_R05 configuration
-            !                                     ! =======================
-            IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea'
-            ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe
-            ij0 = 314   ;   ij1 = 370 
-         END SELECT
-         !
-         mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe
-         !
-      ENDIF
-      !
-   END SUBROUTINE zgr_bat_zoom
-
 
    SUBROUTINE zgr_bat_ctl
Index: /NEMO/branches/2019/ENHANCE-03_domcfg/src/nemogcm.F90
===================================================================
--- /NEMO/branches/2019/ENHANCE-03_domcfg/src/nemogcm.F90	(revision 11200)
+++ /NEMO/branches/2019/ENHANCE-03_domcfg/src/nemogcm.F90	(revision 11201)
@@ -149,5 +149,5 @@
       NAMELIST/namcfg/ ln_e3_dep,                                &
          &             cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, &
-         &             jpizoom, jpjzoom, jperio, ln_use_jattr
+         &             jperio, ln_use_jattr, ln_domclo
       !!----------------------------------------------------------------------
       !
@@ -315,5 +315,4 @@
       jsplt     = nn_jsplt
 
-     !  IF( .NOT.ln_read_cfg )   ln_closea = .false.   ! dealing possible only with a domcfg file
       !
       !                             ! Parameter control
Index: /NEMO/branches/2019/ENHANCE-03_domcfg/src/par_oce.f90
===================================================================
--- /NEMO/branches/2019/ENHANCE-03_domcfg/src/par_oce.f90	(revision 11200)
+++ /NEMO/branches/2019/ENHANCE-03_domcfg/src/par_oce.f90	(revision 11201)
@@ -13,9 +13,5 @@
    PUBLIC
 
-   ! zoom starting position
-   INTEGER       ::   jpizoom          !: left bottom (i,j) indices of the zoom
-   INTEGER       ::   jpjzoom          !: in data domain indices
-
-  CHARACTER(lc) ::   cp_cfg           !: name of the configuration
+   CHARACTER(lc) ::   cp_cfg           !: name of the configuration
    CHARACTER(lc) ::   cp_cfz           !: name of the zoom of configuration
    INTEGER       ::   jp_cfg           !: resolution of the configuration
