Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 9161)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 9161)
@@ -0,0 +1,471 @@
+MODULE closea
+   !!======================================================================
+   !!                   ***  MODULE  closea  ***
+   !!
+   !! User define : specific treatments associated with closed seas
+   !!======================================================================
+   !! History :   8.2  !  2000-05  (O. Marti)  Original code
+   !!   NEMO      1.0  !  2002-06  (E. Durand, G. Madec)  F90
+   !!             3.0  !  2006-07  (G. Madec)  add clo_rnf, clo_ups, clo_bat
+   !!             3.4  !  2014-12  (P.G. Fogli) sbc_clo bug fix & mpp reproducibility
+   !!             4.0  !  2016-06  (G. Madec)  move to usrdef_closea, remove clo_ups
+   !!             4.0  !  2017-12  (D. Storkey) new formulation based on masks read from file
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dom_clo    : read in masks which define closed seas and runoff areas
+   !!   sbc_clo    : Special handling of freshwater fluxes over closed seas
+   !!   clo_rnf    : set close sea outflows as river mouths (see sbcrnf)
+   !!   clo_bat    : set to zero a field over closed sea (see domzgr)
+   !!----------------------------------------------------------------------
+   USE oce             ! dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE sbc_oce         ! ocean surface boundary conditions
+   USE iom             ! I/O routines
+   !
+   USE in_out_manager  ! I/O manager
+   USE lib_fortran,    ONLY: glob_sum
+   USE lbclnk          ! lateral boundary condition - MPP exchanges
+   USE lib_mpp         ! MPP library
+   USE timing          ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC dom_clo      ! called by domain module
+   PUBLIC sbc_clo      ! called by sbcmod module
+   PUBLIC clo_rnf      ! called by sbcrnf module
+   PUBLIC clo_bat      ! called in domzgr module
+
+   LOGICAL, PUBLIC :: ln_closea  !:  T => keep closed seas (defined by closea_mask field) in the domain and apply
+                                 !:       special treatment of freshwater fluxes.
+                                 !:  F => suppress closed seas (defined by closea_mask field) from the bathymetry
+                                 !:       at runtime.
+                                 !:  If there is no closea_mask field in the domain_cfg file or we do not use
+                                 !:  a domain_cfg file then this logical does nothing.
+                                 !:
+   LOGICAL, PUBLIC :: l_sbc_clo  !: T => Closed seas defined, apply special treatment of freshwater fluxes.
+                                 !: F => No closed seas defined (closea_mask field not found).
+   LOGICAL, PUBLIC :: l_clo_rnf  !: T => Some closed seas output freshwater (RNF or EMPMR) to specified runoff points.
+   INTEGER, PUBLIC :: jncs       !: number of closed seas (inferred from closea_mask field)
+   INTEGER, PUBLIC :: jncsr      !: number of closed seas rnf mappings (inferred from closea_mask_rnf field)
+   INTEGER, PUBLIC :: jncse      !: number of closed seas empmr mappings (inferred from closea_mask_empmr field)
+   
+   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::  closea_mask       !: mask of integers defining closed seas
+   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::  closea_mask_rnf   !: mask of integers defining closed seas rnf mappings
+   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::  closea_mask_empmr !: mask of integers defining closed seas empmr mappings
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)  ::   surf         !: closed sea surface areas 
+                                                                  !: (and residual global surface area) 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)  ::   surfr        !: closed sea target rnf surface areas 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)  ::   surfe        !: closed sea target empmr surface areas 
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
+   !! $Id: usrdef_closea.F90 9124 2017-12-19 08:26:25Z gm $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dom_clo()
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_clo  ***
+      !!        
+      !! ** Purpose :   Closed sea domain initialization
+      !!
+      !! ** Method  :   if a closed sea is located only in a model grid point
+      !!                just the thermodynamic processes are applied.
+      !!
+      !! ** Action  :   Read closea_mask* fields (if they exist) from domain_cfg file and infer
+      !!                number of closed seas from closea_mask field.
+      !!                closea_mask       : integer values defining closed seas (or groups of closed seas)
+      !!                closea_mask_rnf   : integer values defining mappings from closed seas or groups of
+      !!                                    closed seas to a runoff area for downwards flux only.
+      !!                closea_mask_empmr : integer values defining mappings from closed seas or groups of
+      !!                                    closed seas to a runoff area for net fluxes.
+      !!
+      !!                Python code to generate the closea_masks* fields from the old-style indices
+      !!                definitions is available at TOOLS/DOMAINcfg/make_closea_masks.py
+      !!----------------------------------------------------------------------
+      INTEGER ::   inum    ! input file identifier
+      INTEGER ::   ierr    ! error code
+      INTEGER ::   id      ! netcdf variable ID
+
+      REAL(wp), DIMENSION(jpi,jpj) :: zdata_in ! temporary real array for input
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*)'dom_clo : read in masks to define closed seas '
+      IF(lwp) WRITE(numout,*)'~~~~~~~'
+      !
+      ! read the closed seas masks (if they exist) from domain_cfg file
+      ! ---------------------------------------------------------------
+      !
+      CALL iom_open( cn_domcfg, inum )
+
+      id = iom_varid(inum, 'closea_mask', ldstop = .false.)
+      IF( id > 0 ) THEN 
+         l_sbc_clo = .true.
+         ALLOCATE( closea_mask(jpi,jpj) , STAT=ierr )
+         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask array')
+         zdata_in(:,:) = 0.0
+         CALL iom_get ( inum, jpdom_data, 'closea_mask', zdata_in )
+         closea_mask(:,:) = NINT(zdata_in(:,:)) * tmask(:,:,1)
+         ! number of closed seas = global maximum value in closea_mask field
+         jncs = maxval(closea_mask(:,:))
+         IF( lk_mpp ) CALL mpp_max(jncs)
+         IF( jncs > 0 ) THEN
+            IF( lwp ) WRITE(numout,*) 'Number of closed seas : ',jncs
+         ELSE
+            CALL ctl_stop( 'Problem with closea_mask field in domain_cfg file. Has no values > 0 so no closed seas defined.')
+         ENDIF
+      ELSE 
+         IF( lwp ) WRITE(numout,*) 'closea_mask field not found in domain_cfg file. No closed seas defined.'
+         l_sbc_clo = .false.
+         jncs = 0 
+      ENDIF
+
+      l_clo_rnf = .false.
+
+      IF( l_sbc_clo ) THEN ! No point reading in closea_mask_rnf or closea_mask_empmr fields if no closed seas defined.
+
+         id = iom_varid(inum, 'closea_mask_rnf', ldstop = .false.)
+         IF( id > 0 ) THEN 
+            l_clo_rnf = .true.            
+            ALLOCATE( closea_mask_rnf(jpi,jpj) , STAT=ierr )
+            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask_rnf array')
+            CALL iom_get ( inum, jpdom_data, 'closea_mask_rnf', zdata_in )
+            closea_mask_rnf(:,:) = NINT(zdata_in(:,:)) * tmask(:,:,1)
+            ! number of closed seas rnf mappings = global maximum in closea_mask_rnf field
+            jncsr = maxval(closea_mask_rnf(:,:))
+            IF( lk_mpp ) CALL mpp_max(jncsr)
+            IF( jncsr > 0 ) THEN
+               IF( lwp ) WRITE(numout,*) 'Number of closed seas rnf mappings : ',jncsr
+            ELSE
+               CALL ctl_stop( 'Problem with closea_mask_rnf field in domain_cfg file. Has no values > 0 so no closed seas rnf mappings defined.')
+            ENDIF
+         ELSE 
+            IF( lwp ) WRITE(numout,*) 'closea_mask_rnf field not found in domain_cfg file. No closed seas rnf mappings defined.'
+            jncsr = 0
+         ENDIF
+ 
+         id = iom_varid(inum, 'closea_mask_empmr', ldstop = .false.)
+         IF( id > 0 ) THEN 
+            l_clo_rnf = .true.            
+            ALLOCATE( closea_mask_empmr(jpi,jpj) , STAT=ierr )
+            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask_empmr array')
+            CALL iom_get ( inum, jpdom_data, 'closea_mask_empmr', zdata_in )
+            closea_mask_empmr(:,:) = NINT(zdata_in(:,:)) * tmask(:,:,1)
+            ! number of closed seas empmr mappings = global maximum value in closea_mask_empmr field
+            jncse = maxval(closea_mask_empmr(:,:))
+            IF( lk_mpp ) CALL mpp_max(jncse)
+            IF( jncse > 0 ) THEN 
+               IF( lwp ) WRITE(numout,*) 'Number of closed seas empmr mappings : ',jncse
+            ELSE
+               CALL ctl_stop( 'Problem with closea_mask_empmr field in domain_cfg file. Has no values > 0 so no closed seas empmr mappings defined.')
+            ENDIF
+         ELSE 
+            IF( lwp ) WRITE(numout,*) 'closea_mask_empmr field not found in domain_cfg file. No closed seas empmr mappings defined.'
+            jncse = 0
+         ENDIF
+
+      ENDIF ! l_sbc_clo
+      !
+      CALL iom_close( inum )
+      !
+   END SUBROUTINE dom_clo
+
+
+   SUBROUTINE sbc_clo( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_clo  ***
+      !!                    
+      !! ** Purpose :   Special handling of closed seas
+      !!
+      !! ** Method  :   Water flux is forced to zero over closed sea
+      !!      Excess is shared between remaining ocean, or
+      !!      put as run-off in open ocean.
+      !!
+      !! ** Action  :   emp updated surface freshwater fluxes and associated heat content at kt
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in   ) ::   kt       ! ocean model time step
+      !
+      INTEGER             ::   ierr
+      INTEGER             ::   jc, jcr, jce   ! dummy loop indices
+      REAL(wp), PARAMETER ::   rsmall = 1.e-20_wp    ! Closed sea correction epsilon
+      REAL(wp)            ::   zfwf_total, zcoef, zcoef1         ! 
+      REAL(wp), DIMENSION(jncs)    ::   zfwf      !:
+      REAL(wp), DIMENSION(jncsr+1) ::   zfwfr     !: freshwater fluxes over closed seas
+      REAL(wp), DIMENSION(jncse+1) ::   zfwfe     !: 
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp2d   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )  CALL timing_start('sbc_clo')
+      !
+      !                                                   !------------------! 
+      IF( kt == nit000 ) THEN                             !  Initialisation  !
+         !                                                !------------------!
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*)'sbc_clo : closed seas '
+         IF(lwp) WRITE(numout,*)'~~~~~~~'
+
+         ALLOCATE( surf(jncs+1) , STAT=ierr )
+         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surf array')
+         surf(:) = 0.e0_wp
+         !
+         ! jncsr can be zero so add 1 to avoid allocating zero-length array
+         ALLOCATE( surfr(jncsr+1) , STAT=ierr )
+         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surfr array')
+         surfr(:) = 0.e0_wp
+         !
+         ! jncse can be zero so add 1 to avoid allocating zero-length array
+         ALLOCATE( surfe(jncse+1) , STAT=ierr )
+         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surfe array')
+         surfe(:) = 0.e0_wp
+         !
+         surf(jncs+1) = glob_sum( e1e2t(:,:) )   ! surface of the global ocean
+         !
+         !                                        ! surface areas of closed seas 
+         DO jc = 1, jncs
+            ztmp2d(:,:) = 0.e0_wp
+            WHERE( closea_mask(:,:) == jc ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:)
+            surf(jc) = glob_sum( ztmp2d(:,:) )
+         END DO
+         !
+         ! jncs+1 : surface area of global ocean, closed seas excluded
+         surf(jncs+1) = surf(jncs+1) - SUM(surf(1:jncs))
+         !
+         !                                        ! surface areas of rnf target areas
+         IF( jncsr > 0 ) THEN
+            DO jcr = 1, jncsr
+               ztmp2d(:,:) = 0.e0_wp
+               WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) == 0 ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:)
+               surfr(jcr) = glob_sum( ztmp2d(:,:) )
+            END DO
+         ENDIF
+         !
+         !                                        ! surface areas of empmr target areas
+         IF( jncse > 0 ) THEN
+            DO jce = 1, jncse
+               ztmp2d(:,:) = 0.e0_wp
+               WHERE( closea_mask_empmr(:,:) == jce .and. closea_mask(:,:) == 0 ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:)
+               surfe(jce) = glob_sum( ztmp2d(:,:) )
+            END DO
+         ENDIF
+         !
+         IF(lwp) WRITE(numout,*)'     Closed sea surface areas (km2)'
+         DO jc = 1, jncs
+            IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jc, surf(jc) * 1.0e-6
+         END DO
+         IF(lwp) WRITE(numout,FMT='(A,ES12.2)') 'Global surface area excluding closed seas (km2): ', surf(jncs+1) * 1.0e-6
+         !
+         IF(jncsr > 0) THEN
+            IF(lwp) WRITE(numout,*)'     Closed sea target rnf surface areas (km2)'
+            DO jcr = 1, jncsr
+               IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jcr, surfr(jcr) * 1.0e-6
+            END DO
+         ENDIF
+         !
+         IF(jncse > 0) THEN
+            IF(lwp) WRITE(numout,*)'     Closed sea target empmr surface areas (km2)'
+            DO jce = 1, jncse
+               IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jce, surfe(jce) * 1.0e-6
+            END DO
+         ENDIF
+      ENDIF
+      !
+      !                                                      !--------------------!
+      !                                                      !  update emp        !
+      !                                                      !--------------------!
+
+      zfwf_total = 0._wp
+
+      !
+      ! 1. Work out total freshwater fluxes over closed seas from EMP - RNF.
+      !
+      zfwf(:) = 0.e0_wp           
+      DO jc = 1, jncs
+         ztmp2d(:,:) = 0.e0_wp
+         WHERE( closea_mask(:,:) == jc ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:)
+         zfwf(jc) = glob_sum( ztmp2d(:,:) )
+      END DO
+      zfwf_total = SUM(zfwf)
+
+      zfwfr(:) = 0.e0_wp           
+      IF( jncsr > 0 ) THEN
+         !
+         ! 2. Work out total FW fluxes over rnf source areas and add to rnf target areas. 
+         !    Where zfwf is negative add flux at specified runoff points and subtract from fluxes for global redistribution.
+         !    Where positive leave in global redistribution total.
+         !
+         DO jcr = 1, jncsr
+            !
+            ztmp2d(:,:) = 0.e0_wp
+            WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) > 0 ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:)
+            zfwfr(jcr) = glob_sum( ztmp2d(:,:) )
+            !
+            ! The following if avoids the redistribution of the round off
+            IF ( ABS(zfwfr(jcr) / surf(jncs+1) ) > rsmall) THEN
+               !
+               ! Add residuals to target runoff points if negative and subtract from total to be added globally
+               IF( zfwfr(jcr) < 0.0 ) THEN 
+                  zfwf_total = zfwf_total - zfwfr(jcr)
+                  zcoef    = zfwfr(jcr) / surfr(jcr)
+                  zcoef1   = rcp * zcoef
+                  WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) == 0.0)
+                     emp(:,:) = emp(:,:) + zcoef
+                     qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
+                  ENDWHERE
+               ENDIF
+               !
+            ENDIF
+         END DO
+      ENDIF  ! jncsr > 0    
+      !
+      zfwfe(:) = 0.e0_wp           
+      IF( jncse > 0 ) THEN
+         !
+         ! 3. Work out total fluxes over empmr source areas and add to empmr target areas. 
+         !
+         DO jce = 1, jncse
+            !
+            ztmp2d(:,:) = 0.e0_wp
+            WHERE( closea_mask_empmr(:,:) == jce .and. closea_mask(:,:) > 0 ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:)
+            zfwfe(jce) = glob_sum( ztmp2d(:,:) )
+            !
+            ! The following if avoids the redistribution of the round off
+            IF ( ABS( zfwfe(jce) / surf(jncs+1) ) > rsmall ) THEN
+               !
+               ! Add residuals to runoff points and subtract from total to be added globally
+               zfwf_total = zfwf_total - zfwfe(jce)
+               zcoef    = zfwfe(jce) / surfe(jce)
+               zcoef1   = rcp * zcoef
+               WHERE( closea_mask_empmr(:,:) == jce .and. closea_mask(:,:) == 0.0)
+                  emp(:,:) = emp(:,:) + zcoef
+                  qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
+               ENDWHERE
+               !
+            ENDIF
+         END DO
+      ENDIF ! jncse > 0    
+
+      !
+      ! 4. Spread residual flux over global ocean. 
+      !
+      ! The following if avoids the redistribution of the round off
+      IF ( ABS(zfwf_total / surf(jncs+1) ) > rsmall) THEN
+         zcoef    = zfwf_total / surf(jncs+1)
+         zcoef1   = rcp * zcoef
+         WHERE( closea_mask(:,:) == 0 )
+            emp(:,:) = emp(:,:) + zcoef
+            qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
+         ENDWHERE
+      ENDIF
+
+      !
+      ! 5. Subtract area means from emp (and qns) over closed seas to give zero mean FW flux over each sea.
+      !
+      DO jc = 1, jncs
+         ! The following if avoids the redistribution of the round off
+         IF ( ABS(zfwf(jc) / surf(jncs+1) ) > rsmall) THEN
+            !
+            ! Subtract residuals from fluxes over closed sea
+            zcoef    = zfwf(jc) / surf(jc)
+            zcoef1   = rcp * zcoef
+            WHERE( closea_mask(:,:) == jc )
+               emp(:,:) = emp(:,:) - zcoef
+               qns(:,:) = qns(:,:) + zcoef1 * sst_m(:,:)
+            ENDWHERE
+            !
+         ENDIF
+      END DO
+      !
+      emp (:,:) = emp (:,:) * tmask(:,:,1)
+      !
+      CALL lbc_lnk( emp , 'T', 1._wp )
+      !
+   END SUBROUTINE sbc_clo
+
+   SUBROUTINE clo_rnf( p_rnfmsk )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_rnf  ***
+      !!                    
+      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
+      !!                to be the same as river mouth grid-points
+      !!
+      !! ** Method  :   set to 1 the runoff mask (mskrnf, see sbcrnf module)
+      !!                at the closed sea outflow grid-point.
+      !!
+      !! ** Action  :   update (p_)mskrnf (set 1 at closed sea outflow)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array)
+      !!----------------------------------------------------------------------
+      !
+      IF( jncsr > 0 ) THEN
+         WHERE( closea_mask_rnf(:,:) > 0 .and. closea_mask(:,:) == 0 )
+            p_rnfmsk(:,:) = MAX( p_rnfmsk(:,:), 1.0_wp )
+         ENDWHERE
+      ENDIF
+      !
+      IF( jncse > 0 ) THEN
+         WHERE( closea_mask_empmr(:,:) > 0 .and. closea_mask(:,:) == 0 )
+            p_rnfmsk(:,:) = MAX( p_rnfmsk(:,:), 1.0_wp )
+         ENDWHERE
+      ENDIF
+      !
+   END SUBROUTINE clo_rnf
+   
+      
+   SUBROUTINE clo_bat( k_top, k_bot )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE clo_bat  ***
+      !!                    
+      !! ** Purpose :   Suppress closed sea from the domain
+      !!
+      !! ** Method  :   Read in closea_mask field (if it exists) from domain_cfg file.
+      !!                Where closea_mask > 0 set first and last ocean level to 0
+      !!                (As currently coded you can't define a closea_mask field in 
+      !!                usr_def_zgr).
+      !!
+      !! ** Action  :   set k_top=0 and k_bot=0 over closed seas
+      !!----------------------------------------------------------------------
+      INTEGER, DIMENSION(:,:), INTENT(inout) ::   k_top, k_bot   ! ocean first and last level indices
+      INTEGER                           :: inum, id
+      INTEGER,  DIMENSION(jpi,jpj) :: closea_mask ! closea_mask field
+      REAL(wp), DIMENSION(jpi,jpj) :: zdata_in ! temporary real array for input
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN                     ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'clo_bat : suppression of closed seas'
+         WRITE(numout,*) '~~~~~~~'
+      ENDIF
+      !
+      CALL iom_open( cn_domcfg, inum )
+      !
+      id = iom_varid(inum, 'closea_mask', ldstop = .false.)      
+      IF( id > 0 ) THEN
+         IF( lwp ) WRITE(numout,*) 'Suppressing closed seas in bathymetry based on closea_mask field,'
+         CALL iom_get ( inum, jpdom_data, 'closea_mask', zdata_in )
+         closea_mask(:,:) = NINT(zdata_in(:,:))
+         WHERE( closea_mask(:,:) > 0 )
+            k_top(:,:) = 0   
+            k_bot(:,:) = 0   
+         ENDWHERE
+      ELSE
+         IF( lwp ) WRITE(numout,*) 'No closea_mask field found in domain_cfg file. No suppression of closed seas.'
+      ENDIF
+      !
+      ! Initialise l_sbc_clo and l_clo_rnf for this case (ln_closea=.false.)
+      l_sbc_clo = .false.
+      l_clo_rnf = .false.
+      !
+      CALL iom_close(inum)
+      !
+   END SUBROUTINE clo_bat
+
+   !!======================================================================
+END MODULE closea
+
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 9161)
@@ -31,5 +31,4 @@
    !                                   !!* Namelist namdom : time & space domain *
    LOGICAL , PUBLIC ::   ln_linssh      !: =T  linear free surface ==>> model level are fixed in time
-   INTEGER , PUBLIC ::   nn_closea      !: =0 suppress closed sea/lake from the ORCA domain or not (=1)
    INTEGER , PUBLIC ::   nn_msh         !: >0  create a mesh-mask file (mesh_mask.nc)
    REAL(wp), PUBLIC ::   rn_isfhmin     !: threshold to discriminate grounded ice to floating ice
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 9161)
@@ -30,5 +30,5 @@
    USE trc_oce        ! shared ocean & passive tracers variab
    USE phycst         ! physical constants
-   USE usrdef_closea  ! closed seas
+   USE closea         ! closed seas
    USE domhgr         ! domain: set the horizontal mesh
    USE domzgr         ! domain: set the vertical mesh
@@ -118,9 +118,9 @@
       CALL dom_glo                     ! global domain versus local domain
       CALL dom_nam                     ! read namelist ( namrun, namdom )
-      CALL dom_clo( cn_cfg, nn_cfg )   ! Closed seas and lake
       CALL dom_hgr                     ! Horizontal mesh
       CALL dom_zgr( ik_top, ik_bot )   ! Vertical mesh and bathymetry
-      IF( nn_closea == 0 )   CALL clo_bat( ik_top, ik_bot )    !==  remove closed seas or lakes  ==!
       CALL dom_msk( ik_top, ik_bot )   ! Masks
+      IF( ln_closea ) CALL dom_clo     ! ln_closea=T : closed seas included in the simulation
+                                       ! Read in masks to define closed seas and lakes 
       !
       DO jj = 1, jpj                   ! depth of the iceshelves
@@ -274,5 +274,5 @@
          &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     &
          &             ln_cfmeta, ln_iscpl
-      NAMELIST/namdom/ ln_linssh, nn_closea, nn_msh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs
+      NAMELIST/namdom/ ln_linssh, nn_msh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs
 #if defined key_netcdf4
       NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
@@ -383,5 +383,4 @@
          WRITE(numout,*) '   Namelist namdom : space & time domain'
          WRITE(numout,*) '      linear free surface (=T)              ln_linssh  = ', ln_linssh
-         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea  = ', nn_closea
          WRITE(numout,*) '      create mesh/mask file(s)              nn_msh     = ', nn_msh
          WRITE(numout,*) '           = 0   no file created           '
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 9161)
@@ -29,4 +29,5 @@
    USE dom_oce        ! ocean domain
    USE usrdef_zgr     ! user defined vertical coordinate system
+   USE closea         ! closed seas
    USE depth_e3       ! depth <=> e3
    USE wet_dry,   ONLY: ll_wd, ssh_ref  ! Wetting and drying
@@ -114,4 +115,8 @@
          gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
       END DO
+      !
+      ! Any closed seas (defined by closea_mask > 0 in domain_cfg file) to be filled 
+      ! in at runtime if ln_closea=.false.
+      IF( .not. ln_closea ) CALL clo_bat( k_top, k_bot )
       !
       IF(lwp) THEN                     ! Control print
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 9161)
@@ -50,5 +50,5 @@
    USE bdy_oce   , ONLY: ln_bdy
    USE usrdef_sbc     ! user defined: surface boundary condition
-   USE usrdef_closea  ! user defined: closed sea
+   USE closea         ! closed sea
    !
    USE prtctl         ! Print control                    (prt_ctl routine)
@@ -151,5 +151,4 @@
          WRITE(numout,*) '         runoff / runoff mouths                     ln_rnf        = ', ln_rnf
          WRITE(numout,*) '         iceshelf formulation                       ln_isf        = ', ln_isf
-         WRITE(numout,*) '         closed sea (=0/1) (set in namdom)          nn_closea     = ', nn_closea
          WRITE(numout,*) '         nb of iterations if land-sea-mask applied  nn_lsm        = ', nn_lsm
          WRITE(numout,*) '         surface wave                               ln_wave       = ', ln_wave
@@ -435,7 +434,7 @@
       IF( nn_fwb    /= 0 )   CALL sbc_fwb( kt, nn_fwb, nn_fsbc )  ! control the freshwater budget
 
-      ! treatment of closed sea in the model domain   (update freshwater fluxes)
-      ! Should not be ran if ln_diurnal_only
-      IF( .NOT.ln_diurnal_only .AND. nn_closea == 1 )   CALL sbc_clo( kt, cn_cfg, nn_cfg )   
+      ! Special treatment of freshwater fluxes over closed seas in the model domain
+      ! Should not be run if ln_diurnal_only
+      IF( l_sbc_clo .AND. (.NOT. ln_diurnal_only) )   CALL sbc_clo( kt )   
 
 !RBbug do not understand why see ticket 667
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 9161)
@@ -21,5 +21,5 @@
    USE sbcisf         ! PM we could remove it I think
    USE eosbn2         ! Equation Of State
-   USE usrdef_closea  ! closed seas
+   USE closea         ! closed seas
    !
    USE in_out_manager ! I/O manager
@@ -510,5 +510,5 @@
       CALL iom_close( inum )                                      ! close file
       !
-      IF( nn_closea == 1 )   CALL clo_rnf( rnfmsk )               ! closed sea inflow set as ruver mouth
+      IF( l_clo_rnf )   CALL clo_rnf( rnfmsk )   ! closed sea inflow set as river mouth
       !
       rnfmsk_z(:)   = 0._wp                                       ! vertical structure
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_closea.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_closea.F90	(revision 9125)
+++ 	(revision )
@@ -1,398 +1,0 @@
-MODULE usrdef_closea
-   !!======================================================================
-   !!                   ***  MODULE  usrdef_closea  ***
-   !!
-   !!                      ===  ORCA configuration  ===
-   !!                         (2, 1 and 1/4 degrees)
-   !!
-   !! User define : specific treatments associated with closed seas
-   !!======================================================================
-   !! History :   8.2  !  2000-05  (O. Marti)  Original code
-   !!   NEMO      1.0  !  2002-06  (E. Durand, G. Madec)  F90
-   !!             3.0  !  2006-07  (G. Madec)  add clo_rnf, clo_ups, clo_bat
-   !!             3.4  !  2014-12  (P.G. Fogli) sbc_clo bug fix & mpp reproducibility
-   !!             4.0  !  2016-06  (G. Madec)  move to usrdef_closea, remove clo_ups
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   dom_clo    : modification of the ocean domain for closed seas cases
-   !!   sbc_clo    : Special handling of closed seas
-   !!   clo_rnf    : set close sea outflows as river mouths (see sbcrnf)
-   !!   clo_bat    : set to zero a field over closed sea (see domzrg)
-   !!----------------------------------------------------------------------
-   USE oce             ! dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE sbc_oce         ! ocean surface boundary conditions
-   !
-   USE in_out_manager  ! I/O manager
-   USE lib_fortran,    ONLY: glob_sum, DDPDD
-   USE lbclnk          ! lateral boundary condition - MPP exchanges
-   USE lib_mpp         ! MPP library
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC dom_clo      ! called by domain module
-   PUBLIC sbc_clo      ! called by step module
-   PUBLIC clo_rnf      ! called by sbcrnf module
-   PUBLIC clo_bat      ! called in domzgr module
-
-   INTEGER, PUBLIC, PARAMETER          ::   jpncs   = 4      !: number of closed sea
-   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncstt            !: Type of closed sea
-   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi1, ncsj1     !: south-west closed sea limits (i,j)
-   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi2, ncsj2     !: north-east closed sea limits (i,j)
-   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsnr            !: number of point where run-off pours
-   INTEGER, PUBLIC, DIMENSION(jpncs,4) ::   ncsir, ncsjr     !: Location of runoff
-
-   REAL(wp), DIMENSION (jpncs+1)       ::   surf             ! closed sea surface
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE dom_clo( cd_cfg, kcfg )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE dom_clo  ***
-      !!        
-      !! ** Purpose :   Closed sea domain initialization
-      !!
-      !! ** Method  :   if a closed sea is located only in a model grid point
-      !!                just the thermodynamic processes are applied.
-      !!
-      !! ** Action  :   ncsi1(), ncsj1() : south-west closed sea limits (i,j)
-      !!                ncsi2(), ncsj2() : north-east Closed sea limits (i,j)
-      !!                ncsir(), ncsjr() : Location of runoff
-      !!                ncsnr            : number of point where run-off pours
-      !!                ncstt            : Type of closed sea
-      !!                                   =0 spread over the world ocean
-      !!                                   =2 put at location runoff
-      !!----------------------------------------------------------------------
-      CHARACTER(len=*), INTENT(in   ) ::   cd_cfg   ! configuration name
-      INTEGER         , INTENT(in   ) ::   kcfg     ! configuration identifier 
-      !
-      INTEGER ::   jc      ! dummy loop indices
-      INTEGER ::   isrow   ! local index
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*)'dom_clo : closed seas '
-      IF(lwp) WRITE(numout,*)'~~~~~~~'
-      !
-      ! initial values
-      ncsnr(:) = 1  ;  ncsi1(:) = 1  ;  ncsi2(:) = 1  ;  ncsir(:,:) = 1
-      ncstt(:) = 0  ;  ncsj1(:) = 1  ;  ncsj2(:) = 1  ;  ncsjr(:,:) = 1
-      !
-      ! set the closed seas (in data domain indices)
-      ! -------------------
-      !
-      IF( cd_cfg == "orca" ) THEN      !==  ORCA configuration  ==!
-         !
-         SELECT CASE ( kcfg )
-         !                                           ! =======================
-         CASE ( 1 )                                  !  ORCA_R1 configuration
-            !                                        ! =======================
-            IF(lwp) WRITE(numout,*)'   ORCA_R1 closed seas :  only the Caspian Sea'
-            ! This dirty section will be suppressed by simplification process:
-            ! all this will come back in input files
-            ! Currently these hard-wired indices relate to configuration with
-            ! extend grid (jpjglo=332)
-            isrow = 332 - jpjglo
-            !
-            ncsnr(1)   = 1    ; ncstt(1)   = 0           ! Caspian Sea  (spread over the globe)
-            ncsi1(1)   = 332  ; ncsj1(1)   = 243 - isrow
-            ncsi2(1)   = 344  ; ncsj2(1)   = 275 - isrow
-            ncsir(1,1) = 1    ; ncsjr(1,1) = 1
-            !                                        
-            !                                        ! =======================
-         CASE ( 2 )                                  !  ORCA_R2 configuration
-            !                                        ! =======================
-            IF(lwp) WRITE(numout,*)'   ORCA_R2 closed seas and lakes : '
-            !                                            ! Caspian Sea
-            IF(lwp) WRITE(numout,*)'      Caspian Sea  '
-            ncsnr(1)   =   1  ;  ncstt(1)   =   0           ! spread over the globe
-            ncsi1(1)   =  11  ;  ncsj1(1)   = 103
-            ncsi2(1)   =  17  ;  ncsj2(1)   = 112
-            ncsir(1,1) =   1  ;  ncsjr(1,1) =   1 
-            !                                            ! Great North American Lakes
-            IF(lwp) WRITE(numout,*)'      Great North American Lakes  '
-            ncsnr(2)   =   1  ;  ncstt(2)   =   2           ! put at St Laurent mouth
-            ncsi1(2)   =  97  ;  ncsj1(2)   = 107
-            ncsi2(2)   = 103  ;  ncsj2(2)   = 111
-            ncsir(2,1) = 110  ;  ncsjr(2,1) = 111           
-            !                                            ! Black Sea (crossed by the cyclic boundary condition)
-            IF(lwp) WRITE(numout,*)'      Black Sea  '
-            ncsnr(3:4) =   4  ;  ncstt(3:4) =   2           ! put in Med Sea (north of Aegean Sea)
-            ncsir(3:4,1) = 171;  ncsjr(3:4,1) = 106         !
-            ncsir(3:4,2) = 170;  ncsjr(3:4,2) = 106 
-            ncsir(3:4,3) = 171;  ncsjr(3:4,3) = 105 
-            ncsir(3:4,4) = 170;  ncsjr(3:4,4) = 105 
-            ncsi1(3)   = 174  ;  ncsj1(3)   = 107           ! 1 : west part of the Black Sea      
-            ncsi2(3)   = 181  ;  ncsj2(3)   = 112           !            (ie west of the cyclic b.c.)
-            ncsi1(4)   =   2  ;  ncsj1(4)   = 107           ! 2 : east part of the Black Sea 
-            ncsi2(4)   =   6  ;  ncsj2(4)   = 112           !           (ie east of the cyclic b.c.)
-            !
-            !                                        ! =========================
-         CASE ( 025 )                                !  ORCA_R025 configuration
-            !                                        ! =========================
-            IF(lwp) WRITE(numout,*)'   ORCA_R025 closed seas : '
-            !                                            ! Caspian Sea
-            IF(lwp) WRITE(numout,*)'      Caspian Sea  '
-            ncsnr(1)   = 1    ; ncstt(1)   = 0               ! Caspian + Aral sea
-            ncsi1(1)   = 1330 ; ncsj1(1)   = 645
-            ncsi2(1)   = 1400 ; ncsj2(1)   = 795
-            ncsir(1,1) = 1    ; ncsjr(1,1) = 1
-            !                                        
-            IF(lwp) WRITE(numout,*)'      Azov Sea  '
-            ncsnr(2)   = 1    ; ncstt(2)   = 0               ! Azov Sea 
-            ncsi1(2)   = 1284 ; ncsj1(2)   = 722
-            ncsi2(2)   = 1304 ; ncsj2(2)   = 747
-            ncsir(2,1) = 1    ; ncsjr(2,1) = 1
-            !
-         END SELECT
-         !
-      ELSE                             !==  No closed sea in the configuration  ==!
-         !
-         IF(lwp) WRITE(numout,*)'   No closed seas or lakes in the configuration '
-         !
-      ENDIF
-
-      ! convert the position in local domain indices
-      ! --------------------------------------------
-      DO jc = 1, jpncs
-         ncsi1(jc)   = mi0( ncsi1(jc) )
-         ncsj1(jc)   = mj0( ncsj1(jc) )
-         !
-         ncsi2(jc)   = mi1( ncsi2(jc) )   
-         ncsj2(jc)   = mj1( ncsj2(jc) )  
-      END DO
-      !
-   END SUBROUTINE dom_clo
-
-
-   SUBROUTINE sbc_clo( kt, cd_cfg, kcfg )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_clo  ***
-      !!                    
-      !! ** Purpose :   Special handling of closed seas
-      !!
-      !! ** Method  :   Water flux is forced to zero over closed sea
-      !!      Excess is shared between remaining ocean, or
-      !!      put as run-off in open ocean.
-      !!
-      !! ** Action  :   emp updated surface freshwater fluxes and associated heat content at kt
-      !!----------------------------------------------------------------------
-      INTEGER         , INTENT(in   ) ::   kt       ! ocean model time step
-      CHARACTER(len=*), INTENT(in   ) ::   cd_cfg   ! configuration name
-      INTEGER         , INTENT(in   ) ::   kcfg     ! configuration identifier 
-      !
-      INTEGER             ::   ji, jj, jc, jn   ! dummy loop indices
-      REAL(wp), PARAMETER ::   rsmall = 1.e-20_wp    ! Closed sea correction epsilon
-      REAL(wp)            ::   zze2, ztmp, zcorr     ! 
-      REAL(wp)            ::   zcoef, zcoef1         ! 
-      COMPLEX(wp)         ::   ctmp 
-      REAL(wp), DIMENSION(jpncs) ::   zfwf   ! 1D workspace
-      !!----------------------------------------------------------------------
-      !
-      !                                                   !------------------!
-      IF( kt == nit000 ) THEN                             !  Initialisation  !
-         !                                                !------------------!
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*)'sbc_clo : closed seas '
-         IF(lwp) WRITE(numout,*)'~~~~~~~'
-
-         surf(:) = 0._wp
-         !
-         surf(jpncs+1) = glob_sum( e1e2t(:,:) )   ! surface of the global ocean
-         !
-         !                                        ! surface of closed seas 
-         DO jc = 1, jpncs
-            ctmp = CMPLX( 0.e0, 0.e0, wp )
-            DO jj = ncsj1(jc), ncsj2(jc)
-               DO ji = ncsi1(jc), ncsi2(jc)
-                  ztmp = e1e2t(ji,jj) * tmask_i(ji,jj)
-                  CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
-               END DO 
-            END DO 
-            IF( lk_mpp )   CALL mpp_sum( ctmp )
-            surf(jc) = REAL(ctmp,wp)
-         END DO
-
-         IF(lwp) WRITE(numout,*)'     Closed sea surfaces'
-         DO jc = 1, jpncs
-            IF(lwp)WRITE(numout,FMT='(1I3,4I4,5X,F16.2)') jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc)
-         END DO
-
-         ! jpncs+1 : surface of sea, closed seas excluded
-         DO jc = 1, jpncs
-            surf(jpncs+1) = surf(jpncs+1) - surf(jc)
-         END DO           
-         !
-      ENDIF
-      !                                                   !--------------------!
-      !                                                   !  update emp        !
-      zfwf = 0.e0_wp                                      !--------------------!
-      DO jc = 1, jpncs
-         ctmp = CMPLX( 0.e0, 0.e0, wp )
-         DO jj = ncsj1(jc), ncsj2(jc)
-            DO ji = ncsi1(jc), ncsi2(jc)
-               ztmp = e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj)
-               CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
-            END DO  
-         END DO 
-         IF( lk_mpp )   CALL mpp_sum( ctmp )
-         zfwf(jc) = REAL(ctmp,wp)
-      END DO
-
-      IF( cd_cfg == "orca" .AND. kcfg == 2 ) THEN      ! Black Sea case for ORCA_R2 configuration
-         zze2    = ( zfwf(3) + zfwf(4) ) * 0.5_wp
-         zfwf(3) = zze2
-         zfwf(4) = zze2
-      ENDIF
-
-      zcorr = 0._wp
-
-      DO jc = 1, jpncs
-         !
-         ! The following if avoids the redistribution of the round off
-         IF ( ABS(zfwf(jc) / surf(jpncs+1) ) > rsmall) THEN
-            !
-            IF( ncstt(jc) == 0 ) THEN           ! water/evap excess is shared by all open ocean
-               zcoef    = zfwf(jc) / surf(jpncs+1)
-               zcoef1   = rcp * zcoef
-               emp(:,:) = emp(:,:) + zcoef
-               qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
-               ! accumulate closed seas correction
-               zcorr    = zcorr    + zcoef
-               !
-            ELSEIF( ncstt(jc) == 1 ) THEN       ! Excess water in open sea, at outflow location, excess evap shared
-               IF ( zfwf(jc) <= 0.e0_wp ) THEN 
-                   DO jn = 1, ncsnr(jc)
-                     ji = mi0(ncsir(jc,jn))
-                     jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
-                     IF (      ji > 1 .AND. ji < jpi   &
-                         .AND. jj > 1 .AND. jj < jpj ) THEN 
-                         zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) )
-                         zcoef1     = rcp * zcoef
-                         emp(ji,jj) = emp(ji,jj) + zcoef
-                         qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)
-                     ENDIF 
-                   END DO 
-               ELSE 
-                   zcoef    = zfwf(jc) / surf(jpncs+1)
-                   zcoef1   = rcp * zcoef
-                   emp(:,:) = emp(:,:) + zcoef
-                   qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
-                   ! accumulate closed seas correction
-                   zcorr    = zcorr    + zcoef
-               ENDIF
-            ELSEIF( ncstt(jc) == 2 ) THEN       ! Excess e-p-r (either sign) goes to open ocean, at outflow location
-               DO jn = 1, ncsnr(jc)
-                  ji = mi0(ncsir(jc,jn))
-                  jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
-                  IF(      ji > 1 .AND. ji < jpi    &
-                     .AND. jj > 1 .AND. jj < jpj ) THEN 
-                     zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) *  e1e2t(ji,jj) )
-                     zcoef1     = rcp * zcoef
-                     emp(ji,jj) = emp(ji,jj) + zcoef
-                     qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)
-                  ENDIF 
-               END DO 
-            ENDIF 
-            !
-            DO jj = ncsj1(jc), ncsj2(jc)
-               DO ji = ncsi1(jc), ncsi2(jc)
-                  zcoef      = zfwf(jc) / surf(jc)
-                  zcoef1     = rcp * zcoef
-                  emp(ji,jj) = emp(ji,jj) - zcoef
-                  qns(ji,jj) = qns(ji,jj) + zcoef1 * sst_m(ji,jj)
-               END DO  
-            END DO 
-            !
-         END IF
-      END DO 
-
-      IF ( ABS(zcorr) > rsmall ) THEN      ! remove the global correction from the closed seas
-         DO jc = 1, jpncs                  ! only if it is large enough
-            DO jj = ncsj1(jc), ncsj2(jc)
-               DO ji = ncsi1(jc), ncsi2(jc)
-                  emp(ji,jj) = emp(ji,jj) - zcorr
-                  qns(ji,jj) = qns(ji,jj) + rcp * zcorr * sst_m(ji,jj)
-               END DO  
-             END DO 
-          END DO
-      ENDIF
-      !
-      emp (:,:) = emp (:,:) * tmask(:,:,1)
-      !
-      CALL lbc_lnk( emp , 'T', 1._wp )
-      !
-   END SUBROUTINE sbc_clo
-
-
-   SUBROUTINE clo_rnf( p_rnfmsk )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_rnf  ***
-      !!                    
-      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
-      !!                to be the same as river mouth grid-points
-      !!
-      !! ** Method  :   set to 1 the runoff mask (mskrnf, see sbcrnf module)
-      !!                at the closed sea outflow grid-point.
-      !!
-      !! ** Action  :   update (p_)mskrnf (set 1 at closed sea outflow)
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array)
-      !
-      INTEGER  ::   jc, jn, ji, jj      ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      DO jc = 1, jpncs
-         IF( ncstt(jc) >= 1 ) THEN            ! runoff mask set to 1 at closed sea outflows
-             DO jn = 1, 4
-                DO jj =    mj0( ncsjr(jc,jn) ), mj1( ncsjr(jc,jn) )
-                   DO ji = mi0( ncsir(jc,jn) ), mi1( ncsir(jc,jn) )
-                      p_rnfmsk(ji,jj) = MAX( p_rnfmsk(ji,jj), 1.0_wp )
-                   END DO
-                END DO
-            END DO 
-         ENDIF 
-      END DO 
-      !
-   END SUBROUTINE clo_rnf
-   
-      
-   SUBROUTINE clo_bat( k_top, k_bot )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE clo_bat  ***
-      !!                    
-      !! ** Purpose :   suppress closed sea from the domain
-      !!
-      !! ** Method  :   set first and last ocean level to 0 over the closed seas.
-      !!
-      !! ** Action  :   set pbat=0 and kbat=0 over closed seas
-      !!----------------------------------------------------------------------
-      INTEGER, DIMENSION(:,:), INTENT(inout) ::   k_top, k_bot   ! ocean first and last level indices
-      !
-      INTEGER  ::   jc, ji, jj      ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      DO jc = 1, jpncs
-         DO jj = ncsj1(jc), ncsj2(jc)
-            DO ji = ncsi1(jc), ncsi2(jc)
-               k_top(ji,jj) = 0   
-               k_bot(ji,jj) = 0   
-            END DO 
-         END DO 
-       END DO 
-       !
-   END SUBROUTINE clo_bat
-
-   !!======================================================================
-END MODULE usrdef_closea
-
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 9161)
@@ -48,4 +48,5 @@
    USE phycst         ! physical constant                  (par_cst routine)
    USE domain         ! domain initialization   (dom_init & dom_cfg routines)
+   USE closea         ! treatment of closed seas (for ln_closea)
    USE usrdef_nam     ! user defined configuration
    USE tideini        ! tidal components initialization   (tide_ini routine)
@@ -242,5 +243,5 @@
          &             nn_isplt , nn_jsplt, nn_jctls, nn_jctle,   &
          &             ln_timing, ln_diacfl
-      NAMELIST/namcfg/ ln_read_cfg, cn_domcfg, ln_write_cfg, cn_domcfg_out, ln_use_jattr
+      NAMELIST/namcfg/ ln_read_cfg, cn_domcfg, ln_write_cfg, cn_domcfg_out, ln_use_jattr, ln_closea
       !!----------------------------------------------------------------------
       !
