Index: branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
===================================================================
--- branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 9021)
+++ branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 9078)
@@ -14,8 +14,8 @@
 
    !!----------------------------------------------------------------------
-   !!   dom_clo    : modification of the ocean domain for closed seas cases
-   !!   sbc_clo    : Special handling of closed seas
+   !!   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 domzrg)
+   !!   clo_bat    : set to zero a field over closed sea (see domzgr)
    !!----------------------------------------------------------------------
    USE oce             ! dynamics and tracers
@@ -38,9 +38,19 @@
    PUBLIC sbc_clo      ! called by sbcmod module
    PUBLIC clo_rnf      ! called by sbcrnf module
-   PUBLIC clo_bat      ! called in domain module
-
-   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)
+   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
@@ -61,5 +71,5 @@
 CONTAINS
 
-   SUBROUTINE dom_clo( k_bot )
+   SUBROUTINE dom_clo()
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE dom_clo  ***
@@ -70,5 +80,5 @@
       !!                just the thermodynamic processes are applied.
       !!
-      !! ** Action  :   Read closea mask fields from closea_masks.nc and infer
+      !! ** 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)
@@ -78,74 +88,94 @@
       !!                                    closed seas to a runoff area for net fluxes.
       !!
-      !!                Python code to generate the closea_masks.nc file from the old-style indices
+      !!                Python code to generate the closea_masks* fields from the old-style indices
       !!                definitions is available at TOOLS/DOMAINcfg/make_closea_masks.py
       !!----------------------------------------------------------------------
-      INTEGER, DIMENSION(:,:), INTENT(in)  ::   k_bot   ! ocean last level index (for masking input fields) 
       INTEGER ::   inum    ! input file identifier
       INTEGER ::   ierr    ! error code
       INTEGER ::   id      ! netcdf variable ID
+
       REAL(wp), POINTER, DIMENSION(:,:) :: zdata_in ! temporary real array for input
       !!----------------------------------------------------------------------
       !
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*)'dom_clo : closed seas '
+      IF(lwp) WRITE(numout,*)'dom_clo : read in masks to define closed seas '
       IF(lwp) WRITE(numout,*)'~~~~~~~'
       !
       CALL wrk_alloc( jpi,jpj, zdata_in)
       ! 
-      ! read the closed seas masks
-      ! --------------------------
-      !
-      ALLOCATE( closea_mask(jpi,jpj) , STAT=ierr )
-      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask array')
-      ALLOCATE( closea_mask_rnf(jpi,jpj) , STAT=ierr )
-      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask_rnf array')
-      ALLOCATE( closea_mask_empmr(jpi,jpj) , STAT=ierr )
-      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask_empmr array')
-
-      closea_mask(:,:) = 0
-      closea_mask_rnf(:,:) = 0
-      closea_mask_empmr(:,:) = 0
-
-      CALL iom_open( 'closea_masks.nc', inum )
-
-      zdata_in(:,:) = 0.0
-      CALL iom_get ( inum, jpdom_data, 'closea_mask', zdata_in )
-      closea_mask(:,:) = NINT(zdata_in(:,:))
-      ! necessary because fill values in input fields can confuse things
-      ! we can't multiply by tmask here because it isn't defined yet
-      WHERE( k_bot(:,:) == 0 ) closea_mask(:,:) = 0
-
-      id = iom_varid(inum, 'closea_mask_rnf', ldstop = .false.)
-      IF(lwp) WRITE(numout,*) 'closea_mask_rnf, id : ',id
+      ! 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 
-         CALL iom_get ( inum, jpdom_data, 'closea_mask_rnf', zdata_in )
-         closea_mask_rnf(:,:) = NINT(zdata_in(:,:))
-         WHERE( k_bot(:,:) == 0 ) closea_mask_rnf(:,:) = 0
-      ENDIF
+         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 
-         CALL iom_get ( inum, jpdom_data, 'closea_mask_empmr', zdata_in )
-         closea_mask_empmr(:,:) = NINT(zdata_in(:,:))
-         WHERE( k_bot(:,:) == 0 ) closea_mask_empmr(:,:) = 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 )
-
-      ! number of closed seas = global maximum value in closea_mask field
-      jncs = maxval(closea_mask(:,:))
-      IF( lk_mpp ) CALL mpp_max(jncs)
-      IF( lwp ) WRITE(numout,*) 'Number of closed seas : ',jncs
-      !
-      ! 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( lwp ) WRITE(numout,*) 'Number of closed seas rnf mappings : ',jncsr
-      !
-      ! 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( lwp ) WRITE(numout,*) 'Number of closed seas empmr mappings : ',jncse
       !
       CALL wrk_dealloc(jpi,jpj, zdata_in)
@@ -217,16 +247,20 @@
          !
          !                                        ! surface areas of rnf target areas
-         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
+         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
-         DO jce = 1, jncse
-            ztmp2d(:,:) = 0.e0_wp
-            WHERE( closea_mask_empmr(:,:) == jcr .and. closea_mask(:,:) == 0 ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:)
-            surfe(jcr) = glob_sum( ztmp2d(:,:) )
-         END DO
+         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)'
@@ -268,61 +302,65 @@
       zfwf_total = SUM(zfwf)
 
-      !
-      ! 2. Work out total FW fluxes over rnf source areas and add to rnf target areas. If jncsr is zero does not loop.
-      !    Where zfwf is negative add flux at specified runoff points and subtract from fluxes for global redistribution.
-      !    Where positive leave in global redistribution total.
-      !
       zfwfr(:) = 0.e0_wp           
-      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(:,:) )
-         IF(lwp) WRITE(numout,*) 'closea runoff output: jcr, zfwfr(jcr), ABS(zfwfr(jcr)/surf(jncs+1) : ',jcr,zfwfr(jcr),ABS(zfwfr(jcr)/surf(jncs+1))
-         !
-         ! 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(lwp) WRITE(numout,*) 'closea runoff output: passed roundoff test!'
-            IF( zfwfr(jcr) < 0.0 ) THEN 
-               zfwf_total = zfwf_total - zfwfr(jcr)
-               zcoef    = zfwfr(jcr) / surfr(jcr)
+      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
+                  IF(lwp) WRITE(numout,*) 'NEGATIVE FLUX zfwfr(jcr), zcoef : ',zfwfr(jcr), zcoef
+                  WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) == 0.0)
+                     emp(:,:) = emp(:,:) + zcoef
+                     qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
+                  ENDWHERE
+               ELSE
+                  IF(lwp) WRITE(numout,*) 'POSITIVE FLUX zfwfr(jcr) : ',zfwfr(jcr)
+               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. If jncse is zero does not loop. 
+         !
+         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
-               IF(lwp) WRITE(numout,*) 'closea runoff output: jcr, zcoef: ',jcr, zcoef
-               WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) == 0.0)
+               WHERE( closea_mask_empmr(:,:) == jce .and. closea_mask(:,:) == 0.0)
                   emp(:,:) = emp(:,:) + zcoef
                   qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
                ENDWHERE
+               !
             ENDIF
-            !
-         ENDIF
-      END DO
-    
-      !
-      ! 3. Work out total fluxes over empmr source areas and add to empmr target areas. If jncse is zero does not loop. 
-      !
-      zfwfe(:) = 0.e0_wp           
-      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
-    
+         END DO
+      ENDIF ! jncse > 0    
+
       !
       ! 4. Spread residual flux over global ocean. 
@@ -332,5 +370,4 @@
          zcoef    = zfwf_total / surf(jncs+1)
          zcoef1   = rcp * zcoef
-         IF(lwp) WRITE(numout,*) 'closea global addition: zfwf_total, zcoef: ', zfwf_total, zcoef
          WHERE( closea_mask(:,:) == 0 )
             emp(:,:) = emp(:,:) + zcoef
@@ -349,5 +386,4 @@
             zcoef    = zfwf(jc) / surf(jc)
             zcoef1   = rcp * zcoef
-            IF(lwp) WRITE(numout,*) 'closea subtract mean: jc, zfwf(jc), zcoef: ',jc, zfwf(jc), zcoef
             WHERE( closea_mask(:,:) == jc )
                emp(:,:) = emp(:,:) - zcoef
@@ -386,7 +422,15 @@
       !!----------------------------------------------------------------------
       !
-      WHERE( ( closea_mask_rnf(:,:) > 0 .or. closea_mask_empmr(:,:) > 0 ) .and. closea_mask(:,:) == 0 )
-         p_rnfmsk(:,:) = MAX( p_rnfmsk(:,:), 1.0_wp )
-      ENDWHERE
+      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
@@ -397,17 +441,49 @@
       !!                  ***  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
+      !! ** 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
-      !!----------------------------------------------------------------------
-      !
-      WHERE( closea_mask(:,:) > 0 )
-         k_top(:,:) = 0   
-         k_bot(:,:) = 0   
-      ENDWHERE
+      INTEGER                           :: inum, id
+      INTEGER,  POINTER, DIMENSION(:,:) :: closea_mask ! closea_mask field
+      REAL(wp), POINTER, DIMENSION(:,:) :: 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 wrk_alloc( jpi,jpj, zdata_in)
+         CALL wrk_alloc( jpi,jpj, closea_mask)
+         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
+         CALL wrk_dealloc( jpi,jpj, zdata_in)
+         CALL wrk_dealloc( jpi,jpj, closea_mask)
+      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
Index: branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 9021)
+++ branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 9078)
@@ -31,6 +31,4 @@
    !                                   !!* Namelist namdom : time & space domain *
    LOGICAL , PUBLIC ::   ln_linssh      !: =T  linear free surface ==>> model level are fixed in time
-   LOGICAL , PUBLIC ::   ln_closea      !: =T apply special treatment to closed seas depending on value of nn_closea.
-   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/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 9021)
+++ branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 9078)
@@ -130,9 +130,7 @@
       CALL dom_hgr                     ! Horizontal mesh
       CALL dom_zgr( ik_top, ik_bot )   ! Vertical mesh and bathymetry
-      IF( ln_closea ) THEN    
-          CALL dom_clo( ik_bot )  ! Closed seas and lake 
-          IF( nn_closea == 0 )   CALL clo_bat( ik_top, ik_bot )    !==  remove closed seas or lakes  ==!
-      ENDIF
       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
@@ -288,5 +286,5 @@
          &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     &
          &             ln_cfmeta, ln_iscpl
-      NAMELIST/namdom/ ln_linssh, ln_closea, 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
@@ -397,6 +395,4 @@
          WRITE(numout,*) '   Namelist namdom : space & time domain'
          WRITE(numout,*) '      linear free surface (=T)              ln_linssh  = ', ln_linssh
-         WRITE(numout,*) '      special treatment of closed seas      ln_closea  = ', ln_closea
-         IF(ln_closea) 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/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 9021)
+++ branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 9078)
@@ -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: ln_wd, ht_wd
@@ -118,4 +119,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/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 9021)
+++ branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 9078)
@@ -153,6 +153,4 @@
          WRITE(numout,*) '         runoff / runoff mouths                     ln_rnf        = ', ln_rnf
          WRITE(numout,*) '         iceshelf formulation                       ln_isf        = ', ln_isf
-         WRITE(numout,*) '         special treatment of closed seas (set in namdom)  ln_closea     = ', ln_closea
-         IF(ln_closea) WRITE(numout,*) '         closed sea treatment (=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
@@ -445,7 +443,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. ln_closea .AND. nn_closea == 1 )   CALL sbc_clo( kt )   
+      ! 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/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 9021)
+++ branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 9078)
@@ -513,5 +513,5 @@
       CALL iom_close( inum )                                      ! close file
       !
-      IF( ln_closea .and. nn_closea == 1 )   CALL clo_rnf( rnfmsk )   ! closed sea inflow set as river mouth
+      IF( l_clo_rnf )   CALL clo_rnf( rnfmsk )   ! closed sea inflow set as river mouth
       !
       rnfmsk_z(:)   = 0._wp                                       ! vertical structure
Index: branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 9021)
+++ branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 9078)
@@ -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,   &
          &             nn_timing, nn_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
       !!----------------------------------------------------------------------
       !
