Index: NEMO/branches/2019/ENHANCE-03_closea/src/OCE/DOM/closea.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-03_closea/src/OCE/DOM/closea.F90	(revision 11207)
+++ NEMO/branches/2019/ENHANCE-03_closea/src/OCE/DOM/closea.F90	(revision 11295)
@@ -33,34 +33,27 @@
 
    IMPLICIT NONE
-   PRIVATE
+   PRIVATE read_csmask
+   PRIVATE alloc_csmask
 
    PUBLIC dom_clo      ! called by domain module
-   PUBLIC sbc_clo_init ! called by sbcmod 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_maskcs  !: mask all closed sea
-   LOGICAL, PUBLIC :: ln_mask_csundef  !: mask all closed sea
-   LOGICAL, PUBLIC :: ln_clo_rnf !: closed sea treated as runoff (update rnf mask)
+   LOGICAL, PUBLIC :: ln_maskcs        !: logical to mask all closed sea
+   LOGICAL, PUBLIC :: ln_mask_csundef  !: logical to mask all undefined closed sea
+   LOGICAL, PUBLIC :: ln_clo_rnf       !: closed sea treated as runoff (update rnf mask)
 
    LOGICAL, PUBLIC :: l_sbc_clo  !: T => net evap/precip over closed seas spread outover the globe/river mouth
    LOGICAL, PUBLIC :: l_clo_rnf  !: T => Some closed seas output freshwater (RNF) to specified runoff points.
 
-   INTEGER, PUBLIC :: jncsg      !: number of closed seas global mappings (inferred from closea_mask_glo 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_emp field)
-
-   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)   :: jcsgrpg, jcsgrpr, jcsgrpe
-  
-   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_opnsea, mask_csundef 
+   INTEGER, PUBLIC :: ncsg      !: number of closed seas global mappings (inferred from closea_mask_glo field)
+   INTEGER, PUBLIC :: ncsr      !: number of closed seas rnf    mappings (inferred from closea_mask_rnf field)
+   INTEGER, PUBLIC :: ncse      !: number of closed seas empmr  mappings (inferred from closea_mask_emp field)
+
+   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_opnsea, mask_csundef !: mask defining the open sea and the undefined closed sea
  
-   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) ::  mask_csglo, mask_csgrpglo !: mask of integers defining closed seas
-   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) ::  mask_csrnf, mask_csgrprnf !: mask of integers defining closed seas rnf mappings
-   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) ::  mask_csemp, mask_csgrpemp !: mask of integers defining closed seas empmr mappings
-
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: rsurfsrcg, rsurftrgg !: closed sea target glo surface areas 
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: rsurfsrcr, rsurftrgr !: closed sea target rnf surface areas 
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: rsurfsrce, rsurftrge !: closed sea target emp surface areas 
+   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_csglo, mask_csgrpglo !: mask of integers defining closed seas
+   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_csrnf, mask_csgrprnf !: mask of integers defining closed seas rnf mappings
+   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_csemp, mask_csgrpemp !: mask of integers defining closed seas empmr mappings
 
    !! * Substitutions
@@ -154,7 +147,7 @@
          !
          ! compute number of cs for each cases
-         jncsg = MAXVAL( mask_csglo(:,:) ) ; CALL mpp_max( 'closea', jncsg )
-         jncsr = MAXVAL( mask_csrnf(:,:) ) ; CALL mpp_max( 'closea', jncsr )
-         jncse = MAXVAL( mask_csemp(:,:) ) ; CALL mpp_max( 'closea', jncse )
+         ncsg = MAXVAL( mask_csglo(:,:) ) ; CALL mpp_max( 'closea', ncsg )
+         ncsr = MAXVAL( mask_csrnf(:,:) ) ; CALL mpp_max( 'closea', ncsr )
+         ncse = MAXVAL( mask_csemp(:,:) ) ; CALL mpp_max( 'closea', ncse )
          !
          ! allocate closed sea group masks
@@ -168,227 +161,6 @@
          CALL read_csmask( cn_domcfg, 'mask_csgrpemp', mask_csgrpemp )
          !
-         ! Allocate cs variables (surf)
-         CALL alloc_cssurf( jncsg, rsurfsrcg, rsurftrgg ) 
-         CALL alloc_cssurf( jncsr, rsurfsrcr, rsurftrgr )
-         CALL alloc_cssurf( jncse, rsurfsrce, rsurftrge )
-         !
-         ! Allocate cs group variables (jcsgrp)
-         CALL alloc_csgrp( jncsg, jcsgrpg )
-         CALL alloc_csgrp( jncsr, jcsgrpr )
-         CALL alloc_csgrp( jncse, jcsgrpe )
-         !
       END IF
    END SUBROUTINE dom_clo
-
-   SUBROUTINE sbc_clo_init
-
-      ! compute source surface area
-      CALL get_cssrcsurf( jncsg, mask_csglo, rsurfsrcg )
-      CALL get_cssrcsurf( jncsr, mask_csrnf, rsurfsrcr )
-      CALL get_cssrcsurf( jncse, mask_csemp, rsurfsrce )
-      !
-      ! compute target surface area and group number (jcsgrp) for all cs and cases 
-      ! glo could be simpler but for lisibility, all treated the same way
-      ! It is only done once, so not a big deal
-      CALL get_cstrgsurf( jncsg, mask_csglo, mask_csgrpglo, rsurftrgg, jcsgrpg )
-      CALL get_cstrgsurf( jncsr, mask_csrnf, mask_csgrprnf, rsurftrgr, jcsgrpr )
-      CALL get_cstrgsurf( jncse, mask_csemp, mask_csgrpemp, rsurftrge, jcsgrpe )
-      ! 
-      ! print out in ocean.ouput
-      CALL prt_csctl( jncsg, rsurfsrcg, rsurftrgg, jcsgrpg, 'glo' )
-      CALL prt_csctl( jncsr, rsurfsrcr, rsurftrgr, jcsgrpr, 'rnf' )
-      CALL prt_csctl( jncse, rsurfsrce, rsurftrge, jcsgrpe, 'emp' )
-
-   END SUBROUTINE sbc_clo_init
-
-   SUBROUTINE get_cssrcsurf(kncs, pmaskcs, psurfsrc)
-
-      ! subroutine parameters
-      INTEGER,                 INTENT(in   ) :: kncs
-      INTEGER, DIMENSION(:,:), INTENT(in   ) :: pmaskcs
-
-      REAL(wp), DIMENSION(:)  , INTENT(inout) :: psurfsrc
-
-      ! local variables
-      INTEGER, DIMENSION(jpi,jpj) :: zmsksrc
-      INTEGER :: jcs
-
-      DO jcs = 1,kncs
-         !
-         ! build river mouth mask for this lake
-         WHERE ( pmaskcs == jcs )
-            zmsksrc = 1
-         ELSE WHERE
-            zmsksrc = 0
-         END WHERE
-         !
-         ! compute target area
-         psurfsrc(jcs) = glob_sum('closea', e1e2t(:,:) * zmsksrc(:,:) )
-         !
-      END DO
-
-   END SUBROUTINE
-
-   SUBROUTINE get_cstrgsurf(kncs, pmaskcs, pmaskcsgrp, psurftrg, kcsgrp )
-
-      ! subroutine parameters
-      INTEGER,               INTENT(in   ) :: kncs
-      INTEGER, DIMENSION(:), INTENT(inout) :: kcsgrp
-      INTEGER, DIMENSION(:,:), INTENT(in   ) :: pmaskcs, pmaskcsgrp
-
-      REAL(wp), DIMENSION(:)  , INTENT(inout) :: psurftrg
-
-      ! local variables
-      INTEGER, DIMENSION(jpi,jpj) :: zmskgrp, zmsksrc, zmsktrg
-      INTEGER :: jcs, jtmp
-
-      DO jcs = 1,kncs
-         !
-         ! find group number
-         zmskgrp = pmaskcsgrp
-         zmsksrc = pmaskcs
-         !
-         ! set value where cs is
-         zmsktrg = HUGE(1)
-         WHERE ( zmsksrc == jcs ) zmsktrg = jcs
-         !
-         ! zmsk = HUGE outside the cs number jcs
-         ! ktmp = jcs - group number
-         ! jgrp = group corresponding to the cs jcs
-         zmsktrg = zmsktrg - zmskgrp
-         jtmp = MINVAL(zmsktrg) ; CALL mpp_min('closea',jtmp)
-         kcsgrp(jcs) = jcs - jtmp
-         !
-         ! build river mouth mask for this lake
-         WHERE ( zmskgrp * mask_opnsea == kcsgrp(jcs) )
-            zmsktrg = 1
-         ELSE WHERE
-            zmsktrg = 0
-         END WHERE
-         !
-         ! compute target area
-         psurftrg(jcs) = glob_sum('closea', e1e2t(:,:) * zmsktrg(:,:) )
-         !
-      END DO
-
-   END SUBROUTINE
-
-   SUBROUTINE prt_csctl(kncs, psurfsrc, psurftrg, kcsgrp, pcstype)
-      ! subroutine parameters
-      INTEGER, INTENT(in   ) :: kncs
-      INTEGER, DIMENSION(:), INTENT(in   ) :: kcsgrp
-      !
-      REAL(wp), DIMENSION(:), INTENT(in   ) :: psurfsrc, psurftrg
-      !
-      CHARACTER(256), INTENT(in   ) :: pcstype
-      !
-      ! local variable
-      INTEGER :: jcs
-      
-      IF ( lwp .AND. kncs > 0 ) THEN
-         WRITE(numout,*)''
-         !
-         WRITE(numout,*)'Closed sea target ',TRIM(pcstype),' : '
-         !
-         DO jcs = 1,kncs
-            WRITE(numout,FMT='(3a,i3,a,i3)') ' ',TRIM(pcstype),' closed sea id is ',jcs,' and trg id is : ', kcsgrp(jcs)
-            WRITE(numout,FMT='(a,f12.2)'   ) ' src surface areas (km2) : ', psurfsrc(jcs) * 1.0e-6
-            WRITE(numout,FMT='(a,f12.2)'   ) ' trg surface areas (km2) : ', psurftrg(jcs) * 1.0e-6
-         END DO
-         !
-         WRITE(numout,*)''
-      END IF
-   END SUBROUTINE
-
-   SUBROUTINE sbc_clo( kt ) ! to be move in SBC in a file sbcclo ???
-      !!---------------------------------------------------------------------
-      !!                  ***  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
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_timing )  CALL timing_start('sbc_clo')
-      !
-      ! update emp and qns
-      CALL sbc_csupdate( jncsg, jcsgrpg, mask_csglo, mask_csgrpglo, rsurfsrcg, rsurftrgg, 'glo', mask_opnsea, rsurftrgg )
-      CALL sbc_csupdate( jncsr, jcsgrpr, mask_csrnf, mask_csgrprnf, rsurfsrcr, rsurftrgr, 'rnf', mask_opnsea, rsurftrgg )
-      CALL sbc_csupdate( jncse, jcsgrpe, mask_csemp, mask_csgrpemp, rsurfsrce, rsurftrge, 'emp', mask_opnsea, rsurftrgg )
-      !
-      ! is this really useful ??????
-      emp(:,:) = emp(:,:) * tmask(:,:,1)
-      qns(:,:) = qns(:,:) * tmask(:,:,1)
-      !
-      ! is this really useful ??????
-      CALL lbc_lnk( 'closea', emp , 'T', 1._wp )
-      CALL lbc_lnk( 'closea', qns , 'T', 1._wp )
-      !
-   END SUBROUTINE sbc_clo
-   
-   SUBROUTINE sbc_csupdate(kncs, kcsgrp, pmsk_src, pmsk_trg, psurfsrc, psurftrg, pcstype, pmsk_opnsea, psurf_opnsea)
-
-      ! subroutine parameters
-      INTEGER, INTENT(in   ) :: kncs
-      INTEGER, DIMENSION(:  ), INTENT(in   ) :: kcsgrp
-      INTEGER, DIMENSION(:,:), INTENT(in   ) :: pmsk_src, pmsk_trg, pmsk_opnsea
-      
-      REAL(wp), DIMENSION(:), INTENT(inout) :: psurfsrc, psurftrg, psurf_opnsea
-
-      CHARACTER(256), INTENT(in   ) :: pcstype
-
-      ! local variables
-      INTEGER :: jcs
-      INTEGER, DIMENSION(jpi,jpj) :: zmsk_src, zmsk_trg
-      
-      REAL(wp) :: zcoef, zcoef1, ztmp
-      REAL(wp) :: zcsfwf
-      REAL(wp) :: zsurftrg
-     
-      DO jcs = 1, kncs
-         !!
-         !! 0. get mask of each closed sea
-         zmsk_src(:,:) = 0
-         WHERE ( pmsk_src(:,:) == jcs ) zmsk_src = 1
-         !!
-         !! 1. Work out net freshwater fluxes over each closed seas from EMP - RNF.
-         zcsfwf = glob_sum( 'closea', e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * zmsk_src )
-         !!
-         !! 2. Deal with runoff special case (net evaporation spread globally)
-         IF (pcstype == 'rnf' .AND. zcsfwf > 0) THEN
-            zsurftrg = psurf_opnsea(1)
-            zmsk_trg = pmsk_opnsea
-         ELSE
-            zsurftrg = psurftrg(jcs)
-            zmsk_trg = pmsk_trg
-         END IF
-         zmsk_trg = zmsk_trg * pmsk_opnsea
-         !!
-         !! 3. Add residuals to target points
-         zcoef    = zcsfwf / zsurftrg
-         zcoef1   = rcp * zcoef
-         WHERE( zmsk_trg(:,:) == kcsgrp(jcs) )
-            emp(:,:) = emp(:,:) + zcoef
-            qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
-         ENDWHERE
-         !!
-         !! 4. Subtract residuals from source points
-         zcoef    = zcsfwf / psurfsrc(jcs)
-         zcoef1   = rcp * zcoef
-         WHERE( pmsk_src(:,:) == jcs )
-            emp(:,:) = emp(:,:) - zcoef
-            qns(:,:) = qns(:,:) + zcoef1 * sst_m(:,:)
-         ENDWHERE
-         !!
-      END DO ! jcs
-
-   END SUBROUTINE
-
 
    SUBROUTINE clo_rnf( p_rnfmsk )
@@ -420,5 +192,5 @@
    END SUBROUTINE clo_rnf
       
-   SUBROUTINE clo_bat( k_top, k_bot, p_mask, p_prt )
+   SUBROUTINE clo_bat( k_top, k_bot, k_mask, cd_prt )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE clo_bat  ***
@@ -436,6 +208,6 @@
       !! subroutine parameter
       INTEGER, DIMENSION(:,:), INTENT(inout) ::   k_top, k_bot   ! ocean first and last level indices
-      INTEGER, DIMENSION(:,:), INTENT(in   ) ::   p_mask         ! mask used to mask ktop and k_bot
-      CHARACTER(256),          INTENT(in   ) ::   p_prt          ! text for control print
+      INTEGER, DIMENSION(:,:), INTENT(in   ) ::   k_mask         ! mask used to mask ktop and k_bot
+      CHARACTER(256),          INTENT(in   ) ::   cd_prt         ! text for control print
       !!
       !! local variables
@@ -444,79 +216,43 @@
       IF ( lwp ) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'clo_bat : Suppression closed seas based on ',TRIM(p_prt),' field.'
+         WRITE(numout,*) 'clo_bat : Suppression closed seas based on ',TRIM(cd_prt),' field.'
          WRITE(numout,*) '~~~~~~~'
          WRITE(numout,*)
       ENDIF
       !!
-      k_top(:,:) = k_top(:,:) * p_mask(:,:)
-      k_bot(:,:) = k_bot(:,:) * p_mask(:,:)
+      k_top(:,:) = k_top(:,:) * k_mask(:,:)
+      k_bot(:,:) = k_bot(:,:) * k_mask(:,:)
       !!
    END SUBROUTINE clo_bat
 
-   SUBROUTINE read_csmask(p_file, p_var, p_mskout)
+   SUBROUTINE read_csmask(cd_file, cd_var, k_mskout)
       !
       ! subroutine parameter
-      CHARACTER(256), INTENT(in   ) :: p_file, p_var
-      INTEGER, DIMENSION(:,:), INTENT(inout) :: p_mskout
+      CHARACTER(256), INTENT(in   ) :: cd_file, cd_var     ! netcdf file and variable name
+      INTEGER, DIMENSION(:,:), INTENT(inout) :: k_mskout ! output mask variable
       !
       ! local variables
-      INTEGER :: ics
-      REAL(wp), DIMENSION(jpi,jpj) :: zdta
-      !
-      CALL iom_open ( p_file, ics )
-      CALL iom_get  ( ics, jpdom_data, TRIM(p_var), zdta )
+      INTEGER :: ics                       ! netcdf id
+      REAL(wp), DIMENSION(jpi,jpj) :: zdta ! netcdf data
+      !
+      CALL iom_open ( cd_file, ics )
+      CALL iom_get  ( ics, jpdom_data, TRIM(cd_var), zdta )
       CALL iom_close( ics )
-      p_mskout(:,:) = NINT(zdta(:,:))
+      k_mskout(:,:) = NINT(zdta(:,:))
       !
    END SUBROUTINE read_csmask
 
-   SUBROUTINE alloc_csmask( pmask )
+   SUBROUTINE alloc_csmask( kmask )
       !
       ! subroutine parameter
-      INTEGER, ALLOCATABLE, DIMENSION(:,:), INTENT(inout) :: pmask
+      INTEGER, ALLOCATABLE, DIMENSION(:,:), INTENT(inout) :: kmask
       !
       ! local variables
       INTEGER :: ierr
       !
-      ALLOCATE( pmask(jpi,jpj) , STAT=ierr )
+      ALLOCATE( kmask(jpi,jpj) , STAT=ierr )
       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'alloc_csmask: failed to allocate surf array')
       !
    END SUBROUTINE
 
-
-   SUBROUTINE alloc_cssurf( klen, pvarsrc, pvartrg )
-      !
-      ! subroutine parameter
-      INTEGER, INTENT(in) :: klen
-      REAL(wp), ALLOCATABLE, DIMENSION(:), INTENT(inout) :: pvarsrc, pvartrg 
-      !
-      ! local variables
-      INTEGER :: ierr
-      !
-      ! klen (number of lake) can be zero so use MAX(klen,1) to avoid 0 length array
-      ALLOCATE( pvarsrc(MAX(klen,1)) , pvartrg(MAX(klen,1)) , STAT=ierr )
-      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surf array')
-      ! initialise to 0
-      pvarsrc(:) = 0.e0_wp
-      pvartrg(:) = 0.e0_wp
-   END SUBROUTINE
-
-   SUBROUTINE alloc_csgrp( klen, kvar )
-      !
-      ! subroutine parameter
-      INTEGER, INTENT(in) :: klen
-      INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(inout) :: kvar 
-      !
-      ! local variables
-      INTEGER :: ierr
-      !
-      ! klen (number of lake) can be zero so use MAX(klen,1) to avoid 0 length array
-      ALLOCATE( kvar(MAX(klen,1)) , STAT=ierr )
-      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate group array')
-      ! initialise to 0
-      kvar(:) = 0
-   END SUBROUTINE
-
-   !!======================================================================
 END MODULE closea
-
Index: NEMO/branches/2019/ENHANCE-03_closea/src/OCE/SBC/sbcclo.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-03_closea/src/OCE/SBC/sbcclo.F90	(revision 11295)
+++ NEMO/branches/2019/ENHANCE-03_closea/src/OCE/SBC/sbcclo.F90	(revision 11295)
@@ -0,0 +1,364 @@
+MODULE sbcclo
+   !!======================================================================
+   !!                       ***  MODULE  sbcclo  ***
+   !! Ocean forcing:  closea sea correction
+   !!=====================================================================
+   !! History :  4.1  ! 2019-09  (P. Mathiot) original
+   !!   NEMO 
+   !!----------------------------------------------------------------------
+   !
+   !!----------------------------------------------------------------------
+   !!   sbc_clo       : update emp and qns over target area and source area
+   !!   sbc_clo_init  : initialise all variable needed for closed sea correction
+   !!
+   !!   alloc_cssurf  : allocate closed sea surface array
+   !!   alloc_csgrp   : allocate closed sea group array
+   !!   get_cssrcsurf : compute source surface area
+   !!   get_cstrgsurf : compute target surface area
+   !!   prt_csctl     : closed sea control print
+   !!   sbc_csupdate  : compute net fw from closed sea
+   !!----------------------------------------------------------------------
+   !
+   USE oce             ! dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE closea          ! closed sea 
+   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 lib_mpp         ! MPP library
+   !
+   IMPLICIT NONE
+   !
+   PRIVATE alloc_cssurf
+   PRIVATE alloc_csgrp
+   PRIVATE get_cssrcsurf
+   PRIVATE get_cstrgsurf
+   PRIVATE prt_csctl
+   PRIVATE sbc_csupdate
+   !
+   PUBLIC sbc_clo
+   PUBLIC sbc_clo_init
+   !
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: rsurfsrcg, rsurftrgg !: closed sea target glo surface areas 
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: rsurfsrcr, rsurftrgr !: closed sea target rnf surface areas 
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: rsurfsrce, rsurftrge !: closed sea target emp surface areas 
+   !
+   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: mcsgrpg, mcsgrpr, mcsgrpe !: closed sea group for glo, rnf and emp
+   !
+   CONTAINS
+   !
+   !!----------------------------------------------------------------------
+   !!  Public subroutines
+   !!----------------------------------------------------------------------
+   !
+   SUBROUTINE sbc_clo_init
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_clo_init  ***
+      !!                    
+      !! ** Purpose :  Initialisation of the net fw closed sea correction
+      !!
+      !! ** Method  :  - compute source surface area for each closed sea
+      !!               - defined the group of each closed sea 
+      !!                    (needed to manage multiple closed sea and one target area like great lakes / St Laurent outlet)
+      !!               - compute target surface area and corresponding group for each closed sea
+      !!----------------------------------------------------------------------
+      !
+      ! 0. Allocate cs variables (surf)
+      CALL alloc_cssurf( ncsg, rsurfsrcg, rsurftrgg ) 
+      CALL alloc_cssurf( ncsr, rsurfsrcr, rsurftrgr )
+      CALL alloc_cssurf( ncse, rsurfsrce, rsurftrge )
+      !
+      ! 1. compute source surface area
+      CALL get_cssrcsurf( ncsg, mask_csglo, rsurfsrcg )
+      CALL get_cssrcsurf( ncsr, mask_csrnf, rsurfsrcr )
+      CALL get_cssrcsurf( ncse, mask_csemp, rsurfsrce )
+      !
+      ! 2. Allocate cs group variables (mcsgrp)
+      CALL alloc_csgrp( ncsg, mcsgrpg )
+      CALL alloc_csgrp( ncsr, mcsgrpr )
+      CALL alloc_csgrp( ncse, mcsgrpe )
+      !
+      ! 3. compute target surface area and group number (mcsgrp) for all cs and cases 
+      ! glo could be simpler but for lisibility, all treated the same way
+      ! It is only done once, so not a big deal
+      CALL get_cstrgsurf( ncsg, mask_csglo, mask_csgrpglo, rsurftrgg, mcsgrpg )
+      CALL get_cstrgsurf( ncsr, mask_csrnf, mask_csgrprnf, rsurftrgr, mcsgrpr )
+      CALL get_cstrgsurf( ncse, mask_csemp, mask_csgrpemp, rsurftrge, mcsgrpe )
+      ! 
+      ! 4. print out in ocean.ouput
+      CALL prt_csctl( ncsg, rsurfsrcg, rsurftrgg, mcsgrpg, 'glo' )
+      CALL prt_csctl( ncsr, rsurfsrcr, rsurftrgr, mcsgrpr, 'rnf' )
+      CALL prt_csctl( ncse, rsurfsrce, rsurftrge, mcsgrpe, 'emp' )
+
+   END SUBROUTINE sbc_clo_init
+
+   SUBROUTINE sbc_clo( kt ) ! to be move in SBC in a file sbcclo ???
+      !!---------------------------------------------------------------------
+      !!                  ***  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  : - compute surface freshwater fluxes and associated heat content flux at kt
+      !!              - output closed sea contribution to fw and heat budget
+      !!              - update emp and qns
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in   ) ::   kt       ! ocean model time step
+      !
+      REAL(wp), DIMENSION(jpi,jpj) :: zwcs, zqcs    ! water flux and heat flux correction due to closed seas
+      !!----------------------------------------------------------------------
+      !
+      ! 0. initialisation
+      zwcs(:,:) = 0._wp ; zqcs(:,:) = 0._wp
+      !
+      ! 1. update emp and qns
+      CALL sbc_csupdate( ncsg, mcsgrpg, mask_csglo, mask_csgrpglo, rsurfsrcg, rsurftrgg, 'glo', mask_opnsea, rsurftrgg, zwcs, zqcs )
+      CALL sbc_csupdate( ncsr, mcsgrpr, mask_csrnf, mask_csgrprnf, rsurfsrcr, rsurftrgr, 'rnf', mask_opnsea, rsurftrgg, zwcs, zqcs )
+      CALL sbc_csupdate( ncse, mcsgrpe, mask_csemp, mask_csgrpemp, rsurfsrce, rsurftrge, 'emp', mask_opnsea, rsurftrgg, zwcs, zqcs )
+      !
+      ! 2. ouput closed sea contributions
+      CALL iom_put('wclosea',zwcs)
+      CALL iom_put('qclosea',zqcs)
+      !
+      ! 3. update emp and qns
+      emp(:,:) = emp(:,:) + zwcs(:,:)
+      qns(:,:) = qns(:,:) + zqcs(:,:)
+      !
+   END SUBROUTINE sbc_clo
+   !
+   !!----------------------------------------------------------------------
+   !!  Private subroutines
+   !!----------------------------------------------------------------------
+   !
+   SUBROUTINE get_cssrcsurf(kncs, kmaskcs, psurfsrc)
+      !!-----------------------------------------------------------------------
+      !!                  ***  routine get_cssrcsurf  ***
+      !!
+      !! ** Purpose : compute closed sea (source) surface area
+      !!----------------------------------------------------------------------
+      ! subroutine parameters
+      INTEGER,                 INTENT(in   ) :: kncs           ! closed sea number
+      INTEGER, DIMENSION(:,:), INTENT(in   ) :: kmaskcs        ! closed sea mask
+
+      REAL(wp), DIMENSION(:)  , INTENT(inout) :: psurfsrc      ! source surface area
+
+      ! local variables
+      INTEGER :: jcs                                           ! loop index
+      INTEGER, DIMENSION(jpi,jpj) :: imsksrc                   ! source mask
+      !!----------------------------------------------------------------------
+      !
+      DO jcs = 1,kncs  ! loop over closed seas
+         !
+         ! 0. build river mouth mask for this lake
+         WHERE ( kmaskcs == jcs )
+            imsksrc = 1
+         ELSE WHERE
+            imsksrc = 0
+         END WHERE
+         !
+         ! 1. compute target area
+         psurfsrc(jcs) = glob_sum('closea', e1e2t(:,:) * imsksrc(:,:) )
+         !
+      END DO  ! jcs
+
+   END SUBROUTINE
+
+   SUBROUTINE get_cstrgsurf(kncs, kmaskcs, kmaskcsgrp, psurftrg, kcsgrp )
+      !!-----------------------------------------------------------------------
+      !!                  ***  routine get_cstrgsurf  ***
+      !!
+      !! ** Purpose : compute closed sea (target) surface area
+      !!----------------------------------------------------------------------
+      ! subroutine parameters
+      INTEGER,                 INTENT(in   ) :: kncs                 ! closed sea number
+      INTEGER, DIMENSION(:)  , INTENT(inout) :: kcsgrp               ! closed sea group number
+      INTEGER, DIMENSION(:,:), INTENT(in   ) :: kmaskcs, kmaskcsgrp  ! closed sea and group mask
+
+      REAL(wp), DIMENSION(:)  , INTENT(inout) :: psurftrg            ! target surface area
+
+      ! local variables
+      INTEGER :: jcs, jtmp                                           ! tmp
+      INTEGER, DIMENSION(jpi,jpj) :: imskgrp, imsksrc, imsktrg       ! tmp group, source and target mask
+      !!----------------------------------------------------------------------
+      !
+      DO jcs = 1,kncs  ! loop over closed seas
+         !
+         !! 0. find group number for cs number jcs
+         imskgrp = kmaskcsgrp
+         imsksrc = kmaskcs
+         !
+         ! set cs value where cs is
+         imsktrg = HUGE(1)
+         WHERE ( imsksrc == jcs ) imsktrg = jcs
+         !
+         ! zmsk = HUGE outside the cs number jcs
+         ! ktmp = jcs - group number
+         ! jgrp = group corresponding to the cs jcs
+         imsktrg = imsktrg - imskgrp
+         jtmp = MINVAL(imsktrg) ; CALL mpp_min('closea',jtmp)
+         kcsgrp(jcs) = jcs - jtmp
+         !
+         !! 1. build river mouth mask for this lake
+         WHERE ( imskgrp * mask_opnsea == kcsgrp(jcs) )
+            imsktrg = 1
+         ELSE WHERE
+            imsktrg = 0
+         END WHERE
+         !
+         !! 2. compute target area
+         psurftrg(jcs) = glob_sum('closea', e1e2t(:,:) * imsktrg(:,:) )
+         !
+      END DO ! jcs
+
+   END SUBROUTINE
+
+   SUBROUTINE prt_csctl(kncs, psurfsrc, psurftrg, kcsgrp, cdcstype)
+      !!-----------------------------------------------------------------------
+      !!                  ***  routine prt_csctl  ***
+      !!
+      !! ** Purpose : output information about each closed sea (src id, trg id, src area and trg area)
+      !!----------------------------------------------------------------------
+      ! subroutine parameters
+      INTEGER,               INTENT(in   ) :: kncs                 ! closed sea number                
+      INTEGER, DIMENSION(:), INTENT(in   ) :: kcsgrp               ! closed sea group number
+      !
+      REAL(wp), DIMENSION(:), INTENT(in   ) :: psurfsrc, psurftrg  ! source and target surface area
+      !
+      CHARACTER(256), INTENT(in   ) :: cdcstype  ! closed sea scheme used for redistribution
+      !
+      ! local variable
+      INTEGER :: jcs
+      !!----------------------------------------------------------------------
+      !
+      IF ( lwp .AND. kncs > 0 ) THEN
+         WRITE(numout,*)''
+         !
+         WRITE(numout,*)'Closed sea target ',TRIM(cdcstype),' : '
+         !
+         DO jcs = 1,kncs
+            WRITE(numout,FMT='(3a,i3,a,i3)') ' ',TRIM(cdcstype),' closed sea id is ',jcs,' and trg id is : ', kcsgrp(jcs)
+            WRITE(numout,FMT='(a,f12.2)'   ) ' src surface areas (km2) : ', psurfsrc(jcs) * 1.0e-6
+            WRITE(numout,FMT='(a,f12.2)'   ) ' trg surface areas (km2) : ', psurftrg(jcs) * 1.0e-6
+         END DO
+         !
+         WRITE(numout,*)''
+      END IF
+
+   END SUBROUTINE
+
+   SUBROUTINE sbc_csupdate(kncs, kcsgrp, kmsk_src, kmsk_trg, psurfsrc, psurftrg, cdcstype, kmsk_opnsea, psurf_opnsea, pwcs, pqcs)
+      !!-----------------------------------------------------------------------
+      !!                  ***  routine sbc_csupdate  ***
+      !!
+      !! ** Purpose : - compute the net freshwater fluxes over each closed seas
+      !!              - apply correction to closed sea source/target net fwf accordingly
+      !!----------------------------------------------------------------------
+      ! subroutine parameters
+      INTEGER,                 INTENT(in) :: kncs                              ! closed sea number
+      INTEGER, DIMENSION(:  ), INTENT(in) :: kcsgrp                            ! closed sea group number
+      INTEGER, DIMENSION(:,:), INTENT(in) :: kmsk_src, kmsk_trg, kmsk_opnsea   ! source, target, open ocean mask
+      
+      REAL(wp), DIMENSION(:)  , INTENT(in   ) :: psurfsrc, psurftrg, psurf_opnsea ! source, target and open ocean surface area
+      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pwcs, pqcs                       ! water and heat flux correction due to closed seas
+
+      CHARACTER(256), INTENT(in   ) :: cdcstype  ! closed sea scheme used for redistribution
+
+      ! local variables
+      INTEGER :: jcs  ! loop index over closed sea 
+      INTEGER, DIMENSION(jpi,jpj) :: imsk_src, imsk_trg  ! tmp array source and target closed sea masks
+      
+      REAL(wp) :: zcoef, zcoef1, ztmp ! tmp
+      REAL(wp) :: zcsfwf              ! tmp net fwf over one closed sea
+      REAL(wp) :: zsurftrg            ! tmp target surface area
+      !!----------------------------------------------------------------------
+      !
+      DO jcs = 1, kncs  ! loop over closed seas
+         !
+         !! 0. get mask of each closed sea
+         imsk_src(:,:) = 0
+         WHERE ( kmsk_src(:,:) == jcs ) imsk_src(:,:) = 1
+         !
+         !! 1. Work out net freshwater fluxes over each closed seas from EMP - RNF.
+         zcsfwf = glob_sum( 'closea', e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * imsk_src(:,:) )
+         !
+         !! 2. Deal with runoff special case (net evaporation spread globally)
+         IF (cdcstype == 'rnf' .AND. zcsfwf > 0) THEN
+            zsurftrg = psurf_opnsea(1)
+            imsk_trg = kmsk_opnsea * kcsgrp(jcs) ! set imsk_trg value to the corresponding group id
+         ELSE
+            zsurftrg = psurftrg(jcs)
+            imsk_trg = kmsk_trg
+         END IF
+         imsk_trg = imsk_trg * kmsk_opnsea
+         !
+         !! 3. Add residuals to target points
+         zcoef    = zcsfwf / zsurftrg
+         zcoef1   = rcp * zcoef
+         WHERE( imsk_trg(:,:) == kcsgrp(jcs) )
+            pwcs(:,:) = pwcs(:,:) + zcoef
+            pqcs(:,:) = pqcs(:,:) - zcoef1 * sst_m(:,:)
+         ENDWHERE
+         !
+         !! 4. Subtract residuals from source points
+         zcoef    = zcsfwf / psurfsrc(jcs)
+         zcoef1   = rcp * zcoef
+         WHERE( kmsk_src(:,:) == jcs )
+            pwcs(:,:) = pwcs(:,:) - zcoef
+            pqcs(:,:) = pqcs(:,:) + zcoef1 * sst_m(:,:)
+         ENDWHERE
+         !
+      END DO ! jcs
+
+   END SUBROUTINE
+
+   SUBROUTINE alloc_cssurf( klen, pvarsrc, pvartrg )
+      !!-----------------------------------------------------------------------
+      !!                  ***  routine alloc_cssurf  ***
+      !!
+      !! ** Purpose : allocate closed sea surface array (source)
+      !!----------------------------------------------------------------------
+      ! subroutine parameters
+      INTEGER, INTENT(in) :: klen
+      REAL(wp), ALLOCATABLE, DIMENSION(:), INTENT(inout) :: pvarsrc, pvartrg 
+      !
+      ! local variables
+      INTEGER :: ierr
+      !!----------------------------------------------------------------------
+      !
+      ! klen (number of lake) can be zero so use MAX(klen,1) to avoid 0 length array
+      ALLOCATE( pvarsrc(MAX(klen,1)) , pvartrg(MAX(klen,1)) , STAT=ierr )
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surf array')
+      !
+      ! initialise to 0
+      pvarsrc(:) = 0.e0_wp
+      pvartrg(:) = 0.e0_wp
+   END SUBROUTINE
+
+   SUBROUTINE alloc_csgrp( klen, kvar )
+      !!-----------------------------------------------------------------------
+      !!                  ***  routine alloc_csgrp  ***
+      !!
+      !! ** Purpose : allocate closed sea group surface array
+      !!----------------------------------------------------------------------
+      ! subroutine parameters
+      INTEGER, INTENT(in) :: klen
+      INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(inout) :: kvar 
+      !
+      ! local variables
+      INTEGER :: ierr
+      !!----------------------------------------------------------------------
+      !
+      ! klen (number of lake) can be zero so use MAX(klen,1) to avoid 0 length array
+      ALLOCATE( kvar(MAX(klen,1)) , STAT=ierr )
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate group array')
+      ! initialise to 0
+      kvar(:) = 0
+   END SUBROUTINE
+
+END MODULE
Index: NEMO/branches/2019/ENHANCE-03_closea/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-03_closea/src/OCE/SBC/sbcmod.F90	(revision 11207)
+++ NEMO/branches/2019/ENHANCE-03_closea/src/OCE/SBC/sbcmod.F90	(revision 11295)
@@ -41,4 +41,5 @@
    USE sbccpl         ! surface boundary condition: coupled formulation
    USE cpl_oasis3     ! OASIS routines for coupling
+   USE sbcclo         ! surface boundary condition: closed sea correction
    USE sbcssr         ! surface boundary condition: sea surface restoring
    USE sbcrnf         ! surface boundary condition: runoffs
