Index: NEMO/branches/2019/ENHANCE-03_closea/src/OCE/SBC/sbcclo.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-03_closea/src/OCE/SBC/sbcclo.F90	(revision 11629)
+++ NEMO/branches/2019/ENHANCE-03_closea/src/OCE/SBC/sbcclo.F90	(revision 11656)
@@ -22,14 +22,13 @@
    !!----------------------------------------------------------------------
    !
-   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
+   USE closea                                  ! closed sea 
+   USE in_out_manager                          ! I/O manager
+   !
+   USE dom_oce,     ONLY: e1e2t                ! ocean space and time domain
+   USE phycst ,     ONLY: rcp                  ! physical constants
+   USE sbc_oce,     ONLY: emp, qns, rnf, sst_m ! ocean surface boundary conditions
+   USE iom    ,     ONLY: iom_put              ! I/O routines
+   USE lib_fortran, ONLY: glob_sum             ! fortran library
+   USE lib_mpp    , ONLY: mpp_min, ctl_stop    ! MPP library
    !
    IMPLICIT NONE
@@ -270,22 +269,29 @@
       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
+      REAL(wp) :: zcsfw, zcsh        ! total fresh water and associated heat over one closed sea
+      REAL(wp) :: zcsfwf             ! mean fresh water flux over one closed sea
+      REAL(wp) :: zsurftrg, zsurfsrc ! total target surface area
       !!----------------------------------------------------------------------
       !
       DO jcs = 1, kncs  ! loop over closed seas
          !
-         !! 0. get mask of the closed sea
-         imsk_src(:,:) = 0
-         WHERE ( kmsk_src(:,:) == jcs ) imsk_src(:,:) = 1
-         !
-         !! 1. Work out net freshwater fluxes over the closed sea from EMP - RNF.
+         !! 0. get mask and surface of the closed sea
+         ! mask src
+         WHERE ( kmsk_src(:,:) == jcs ) 
+            imsk_src(:,:) = 1
+         ELSEWHERE
+            imsk_src(:,:) = 0
+         END WHERE
+         ! area src
+         zsurfsrc = psurfsrc(jcs)
+         !
+         !! 1. Work out net freshwater over the closed sea from EMP - RNF.
+         !!    Work out net heat associated with the correction (needed for conservation)
          !!    (PM: should we consider used delayed glob sum ?)
-         zcsfwf = glob_sum( 'closea', e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * imsk_src(:,:) )
+         zcsfw  = glob_sum( 'closea', e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * imsk_src(:,:) )
          !
          !! 2. Deal with runoff special case (net evaporation spread globally)
          !!    and compute trg mask
-         IF (cdcstype == 'rnf' .AND. zcsfwf > 0) THEN
+         IF (cdcstype == 'rnf' .AND. zcsfw  > 0._wp) THEN
             zsurftrg = psurf_opnsea(1)           ! change the target area surface
             imsk_trg = kcsgrp(jcs) * kmsk_opnsea ! trg mask is now the open sea mask
@@ -295,22 +301,17 @@
          END IF
          !
-         !! 3. Add residuals to target points
-         zcoef  = zcsfwf / zsurftrg
-         zcoef1 = rcp * zcoef
+         !! 3. Subtract residuals from source points
+         zcsfwf = zcsfw / zsurfsrc
+         pwcs(:,:) = pwcs(:,:) -       zcsfwf              * imsk_src(:,:)
+         pqcs(:,:) = pqcs(:,:) + rcp * zcsfwf * sst_m(:,:) * imsk_src(:,:)
+         !
+         !! 4. Add residuals to target points 
+         !!    Do not use pqcs(:,:) = pqcs(:,:) - rcp * zcsfw  * sst_m(:,:) / zsurftrg 
+         !!    as there is no reason for heat conservation with this formulation
+         zcsh   = glob_sum( 'closea', e1e2t(:,:) * rcp * zcsfwf * sst_m(:,:) * imsk_src(:,:) )
          WHERE( imsk_trg(:,:) == kcsgrp(jcs) )
-            pwcs(:,:) = pwcs(:,:) + zcoef
-            pqcs(:,:) = pqcs(:,:) - zcoef1 * sst_m(:,:)
+            pwcs(:,:) = pwcs(:,:) + zcsfw / zsurftrg
+            pqcs(:,:) = pqcs(:,:) - zcsh  / zsurftrg
          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
-         ! WARNING (PM): the correction is done as it was done in the previous version
-         !               this do no conserve heat as there is no reason that 
-         !               sum(zcoef1*sst_m) over the source (closed sea) (4) = sum(zcoef1*sst_m) over the target (river mouth) (3)
          !
       END DO ! jcs
