Index: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 5802)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 5820)
@@ -46,5 +46,6 @@
    REAL(wp) ::   frc_wn_t, frc_wn_s    ! global forcing trends
    !
-   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf          , ssh_ini          !
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf 
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf_ini      , ssh_ini          !
    REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   !
    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
@@ -144,5 +145,5 @@
 
       ! volume variation (calculated with ssh)
-      zdiff_v1 = glob_sum( surf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) )
+      zdiff_v1 = glob_sum_full( surf(:,:) * sshn(:,:) - surf_ini(:,:) * ssh_ini(:,:) )
 
       ! heat & salt content variation (associated with ssh)
@@ -165,12 +166,12 @@
       DO jk = 1, jpkm1
          ! volume variation (calculated with scale factors)
-         zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * ( tmask(:,:,jk) &
-            &                           * fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
-         ! heat content variation
-         zdiff_hc = zdiff_hc + glob_sum(  surf(:,:) * ( tmask(:,:,jk) & 
-            &                           * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) )
+         zdiff_v2 = zdiff_v2 + glob_sum_full( surf    (:,:) * fse3t_n  (:,:,jk) * tmask(:,:,jk) &
+              &                             - surf_ini(:,:) *   e3t_ini(:,:,jk) )
+         ! heat content variation 
+         zdiff_hc = zdiff_hc + glob_sum_full( surf    (:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) &
+              &                             - surf_ini(:,:) * hc_loc_ini(:,:,jk) )
          ! salt content variation
-         zdiff_sc = zdiff_sc + glob_sum(  surf(:,:) * ( tmask(:,:,jk)   &
-            &                           * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) )
+         zdiff_sc = zdiff_sc + glob_sum_full( surf    (:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) & 
+              &                             - surf_ini(:,:) * sc_loc_ini(:,:,jk) )
       ENDDO
 
@@ -192,5 +193,5 @@
       zvol_tot = 0._wp                    ! total ocean volume (calculated with scale factors)
       DO jk = 1, jpkm1
-         zvol_tot  = zvol_tot + glob_sum( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
+         zvol_tot  = zvol_tot + glob_sum_full( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
       END DO
 
@@ -268,4 +269,5 @@
               CALL iom_get( numror, 'frc_wn_s', frc_wn_s )
            ENDIF
+           CALL iom_get( numror, jpdom_autoglo, 'surf_ini', surf_ini )
            CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini )
            CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini )
@@ -280,4 +282,5 @@
           IF(lwp) WRITE(numout,*) ' dia_hsb at initial state '
           IF(lwp) WRITE(numout,*) '~~~~~~~'
+          surf_ini(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)             ! initial ocean surface
           ssh_ini(:,:) = sshn(:,:)                                       ! initial ssh
           DO jk = 1, jpk
@@ -320,4 +323,5 @@
            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s )
         ENDIF
+        CALL iom_rstput( kt, nitrst, numrow, 'surf_ini', surf_ini )
         CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini )
         CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini )
@@ -387,6 +391,6 @@
       ! 1 - Allocate memory !
       ! ------------------- !
-      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), &
-         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror )
+      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), &
+         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror  )
       IF( ierror > 0 ) THEN
          CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
Index: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90
===================================================================
--- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90	(revision 5802)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90	(revision 5820)
@@ -62,5 +62,5 @@
       REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
       REAL(wp):: r1_tiscpl
-      REAL(wp):: zjip1_ratio, zjim1_ratio, zjjp1_ratio, zjjm1_ratio
+      REAL(wp):: zjip1_ratio  , zjim1_ratio  , zjjp1_ratio  , zjjm1_ratio
       !!
       REAL(wp), DIMENSION(:,:    ), POINTER :: zde3t, zdtem, zdsal
@@ -83,5 +83,5 @@
       pts_flx (:,:,:,:) = 0.0_wp
 
-      ! mask tsn and tsb (should be useless)
+      ! mask tsn and tsb 
       tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)*ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)*tmask(:,:,:);
       tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)*ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)*tmask(:,:,:);
@@ -98,64 +98,51 @@
                IF (tmask_h(ji,jj) == 1._wp) THEN
 
-               ! volume differences
-               zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
-
-               ! heat diff
-               zdtem(ji,jj) = tsn(ji,jj,jk,jp_tem) * fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)   &
+                  ! volume differences
+                  zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
+
+                  ! heat diff
+                  zdtem(ji,jj) = tsn(ji,jj,jk,jp_tem) * fse3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
                             - tsb(ji,jj,jk,jp_tem) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
-               ! salt diff
-               zdsal(ji,jj) = tsn(ji,jj,jk,jp_sal) * fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)   &
-                            - tsb(ji,jj,jk,jp_sal) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
+                  ! salt diff
+                  zdsal(ji,jj) = tsn(ji,jj,jk,jp_sal) * fse3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
+                               - tsb(ji,jj,jk,jp_sal) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
                
-               ! shh changes
-               IF ( ptmask_b(ji,jj,jk) == 1 .OR. tmask(ji,jj,jk) == 1 ) THEN 
-                  zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj)
-                  zssh0(ji,jj) = 0._wp
-               END IF
-
-               ! volume, heat and salt differences in each cell 
-               pvol_flx(ji,jj,jk)       =   pvol_flx(ji,jj,jk)        + zde3t(ji,jj) * r1_tiscpl
-               pts_flx (ji,jj,jk,jp_sal)=   pts_flx (ji,jj,jk,jp_sal) + zdsal(ji,jj) * r1_tiscpl 
-               pts_flx (ji,jj,jk,jp_tem)=   pts_flx (ji,jj,jk,jp_tem) + zdtem(ji,jj) * r1_tiscpl
-
-               IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
-                  ! case where we close a cell: check if the neighbour cells are wet 
-
-                  jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ;
-
-                  zsum =   e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) + e12t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) &
-                    &    + e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) + e12t(jip1,jj  ) * tmask(jip1,jj  ,jk) 
-
-                  IF ( zsum .NE. 0._wp ) THEN
-                     zjip1_ratio = e12t(jip1,jj  ) * tmask(jip1,jj  ,jk) / zsum
-                     zjim1_ratio = e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) / zsum
-                     zjjp1_ratio = e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) / zsum
-                     zjjm1_ratio = e12t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) / zsum
-
-                     pvol_flx(ji  ,jjp1,jk       ) = pvol_flx(ji  ,jjp1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjp1_ratio
-                     pvol_flx(ji  ,jjm1,jk       ) = pvol_flx(ji  ,jjm1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjm1_ratio
-                     pvol_flx(jip1,jj  ,jk       ) = pvol_flx(jip1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjip1_ratio
-                     pvol_flx(jim1,jj  ,jk       ) = pvol_flx(jim1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjim1_ratio
-                     pts_flx (ji  ,jjp1,jk,jp_sal) = pts_flx (ji  ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio
-                     pts_flx (ji  ,jjm1,jk,jp_sal) = pts_flx (ji  ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio
-                     pts_flx (jip1,jj  ,jk,jp_sal) = pts_flx (jip1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio
-                     pts_flx (jim1,jj  ,jk,jp_sal) = pts_flx (jim1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio
-                     pts_flx (ji  ,jjp1,jk,jp_tem) = pts_flx (ji  ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio
-                     pts_flx (ji  ,jjm1,jk,jp_tem) = pts_flx (ji  ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio
-                     pts_flx (jip1,jj  ,jk,jp_tem) = pts_flx (jip1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio
-                     pts_flx (jim1,jj  ,jk,jp_tem) = pts_flx (jim1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio
-
-                     ! set to 0 the cell we distributed over neigbourg cells
-                     pvol_flx(ji,jj,jk       ) = 0._wp
-                     pts_flx (ji,jj,jk,jp_sal) = 0._wp
-                     pts_flx (ji,jj,jk,jp_tem) = 0._wp
-
-                  ELSE IF (zsum .EQ. 0._wp ) THEN
-                     ! case where we close a cell and no adjacent cell open
-                     ! check if the cell beneath is wet
-                     IF ( tmask(ji,jj,jk+1) .EQ. 1._wp ) THEN
-                        pvol_flx(ji,jj,jk+1)       =  pvol_flx(ji,jj,jk+1)        + pvol_flx(ji,jj,jk)
-                        pts_flx (ji,jj,jk+1,jp_sal)=  pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal)
-                        pts_flx (ji,jj,jk+1,jp_tem)=  pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem)
+                  ! shh changes
+                  IF ( ptmask_b(ji,jj,jk) == 1._wp .OR. tmask(ji,jj,jk) == 1._wp ) THEN 
+                     zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj) ! zssh0 = 0 if vvl
+                     zssh0(ji,jj) = 0._wp
+                  END IF
+
+                  ! volume, heat and salt differences in each cell 
+                  pvol_flx(ji,jj,jk)       =   pvol_flx(ji,jj,jk)        + zde3t(ji,jj) * r1_tiscpl
+                  pts_flx (ji,jj,jk,jp_sal)=   pts_flx (ji,jj,jk,jp_sal) + zdsal(ji,jj) * r1_tiscpl 
+                  pts_flx (ji,jj,jk,jp_tem)=   pts_flx (ji,jj,jk,jp_tem) + zdtem(ji,jj) * r1_tiscpl
+
+                  IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
+                     ! case where we close a cell: check if the neighbour cells are wet 
+
+                     jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ;
+
+                     zsum =   e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) + e12t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) &
+                       &    + e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) + e12t(jip1,jj  ) * tmask(jip1,jj  ,jk)
+
+                     IF ( zsum .NE. 0._wp ) THEN
+                        zjip1_ratio   = e12t(jip1,jj  ) * tmask(jip1,jj  ,jk) / zsum
+                        zjim1_ratio   = e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) / zsum
+                        zjjp1_ratio   = e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) / zsum
+                        zjjm1_ratio   = e12t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) / zsum
+
+                        pvol_flx(ji  ,jjp1,jk       ) = pvol_flx(ji  ,jjp1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjp1_ratio
+                        pvol_flx(ji  ,jjm1,jk       ) = pvol_flx(ji  ,jjm1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjm1_ratio
+                        pvol_flx(jip1,jj  ,jk       ) = pvol_flx(jip1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjip1_ratio
+                        pvol_flx(jim1,jj  ,jk       ) = pvol_flx(jim1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjim1_ratio
+                        pts_flx (ji  ,jjp1,jk,jp_sal) = pts_flx (ji  ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio
+                        pts_flx (ji  ,jjm1,jk,jp_sal) = pts_flx (ji  ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio
+                        pts_flx (jip1,jj  ,jk,jp_sal) = pts_flx (jip1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio
+                        pts_flx (jim1,jj  ,jk,jp_sal) = pts_flx (jim1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio
+                        pts_flx (ji  ,jjp1,jk,jp_tem) = pts_flx (ji  ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio
+                        pts_flx (ji  ,jjm1,jk,jp_tem) = pts_flx (ji  ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio
+                        pts_flx (jip1,jj  ,jk,jp_tem) = pts_flx (jip1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio
+                        pts_flx (jim1,jj  ,jk,jp_tem) = pts_flx (jim1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio
 
                         ! set to 0 the cell we distributed over neigbourg cells
@@ -163,15 +150,28 @@
                         pts_flx (ji,jj,jk,jp_sal) = 0._wp
                         pts_flx (ji,jj,jk,jp_tem) = 0._wp
-                     ELSE
-                     ! case no adjacent cell on the horizontal and on the vertical
-                        PRINT *,        'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
-                        PRINT *,        '                     ',mig(ji),' ',mjg(jj),' ',jk
-                        PRINT *,        '                     ',ji,' ',jj,' ',jk,' ',narea
-                        PRINT *,        ' we are now looking for the closest wet cell on the horizontal '
-                     ! We deal with this points later.
+
+                     ELSE IF (zsum .EQ. 0._wp ) THEN
+                        ! case where we close a cell and no adjacent cell open
+                        ! check if the cell beneath is wet
+                        IF ( tmask(ji,jj,jk+1) .EQ. 1._wp ) THEN
+                           pvol_flx(ji,jj,jk+1)       =  pvol_flx(ji,jj,jk+1)        + pvol_flx(ji,jj,jk)
+                           pts_flx (ji,jj,jk+1,jp_sal)=  pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal)
+                           pts_flx (ji,jj,jk+1,jp_tem)=  pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem)
+
+                           ! set to 0 the cell we distributed over neigbourg cells
+                           pvol_flx(ji,jj,jk       ) = 0._wp
+                           pts_flx (ji,jj,jk,jp_sal) = 0._wp
+                           pts_flx (ji,jj,jk,jp_tem) = 0._wp
+                        ELSE
+                        ! case no adjacent cell on the horizontal and on the vertical
+                           WRITE(numout,*) 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
+                           WRITE(numout,*) '                     ',mig(ji),' ',mjg(jj),' ',jk
+                           WRITE(numout,*) '                     ',ji,' ',jj,' ',jk,' ',narea
+                           WRITE(numout,*) ' we are now looking for the closest wet cell on the horizontal '
+                        ! We deal with these points later.
+                        END IF
                      END IF
                   END IF
                END IF
-               END IF
             END DO
          END DO
@@ -181,11 +181,4 @@
       CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
       CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
-
-      zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt
-      IF (lwp) PRINT *, 'total volume correction 21 = ',zsum
-      zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
-      IF (lwp) PRINT *, 'total heat correction 21 = ',zsum
-      zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
-      IF (lwp) PRINT *, 'total salt correction 21 = ',zsum
 
       ! if no neighbour wet cell in case of 2close a cell", need to find the nearest wet point 
@@ -198,7 +191,6 @@
          DO ji = 2,jpi-1
             DO jj = 2,jpj-1
-               IF (          ptmask_b(ji,jj,jk)      == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 &
-               &   .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND.     tmask(ji,jj,jk+1)       == 0 &
-               &   .AND. tmask_h(ji,jj) == 1._wp ) THEN
+               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
+                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
                   vnpts(narea) = vnpts(narea) + 1 
                END IF
@@ -216,7 +208,7 @@
       ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts))
       ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20 ; zlat(:) = -1.0e20
-      zcorr_vol(:) = 0.0_wp
-      zcorr_sal(:) = 0.0_wp
-      zcorr_tem(:) = 0.0_wp
+      zcorr_vol(:) = -1.0e20
+      zcorr_sal(:) = -1.0e20
+      zcorr_tem(:) = -1.0e20
 
       ! fill new variable
@@ -225,9 +217,7 @@
          DO ji = 2,jpi-1
             DO jj = 2,jpj-1
-               IF (          ptmask_b(ji,jj,jk)      == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 &
-               &   .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND.     tmask(ji,jj,jk+1)       == 0 &
-               &   .AND. tmask_h(ji,jj) == 1 ) THEN
+               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
+                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
                   jpts = jpts + 1  ! positioning in the vnpts vector for the area narea
-                  PRINT *, 'corrected point ', narea, ji, jj, jk, jpts
                   ixpts(jpts) = ji           ; iypts(jpts) = jj ; izpts(jpts) = jk
                   zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj)
@@ -235,9 +225,9 @@
                   zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal)
                   zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem)
+
                   ! set flx to 0 (safer)
                   pvol_flx(ji,jj,jk       ) = 0.0_wp
                   pts_flx (ji,jj,jk,jp_sal) = 0.0_wp
                   pts_flx (ji,jj,jk,jp_tem) = 0.0_wp
-                  PRINT *, zcorr_vol(jpts)*rn_fiscpl*rn_rdt, zcorr_sal(jpts)*rn_fiscpl*rn_rdt, zcorr_tem(jpts)*rn_fiscpl*rn_rdt
                END IF
             END DO
@@ -246,22 +236,36 @@
 
       ! build array of total problematic point on each cpu (share to each cpu)
+      ! point coordinates
       CALL mpp_max(zlat ,npts)
       CALL mpp_max(zlon ,npts)
       CALL mpp_max(izpts,npts)
 
+      ! correction values 
+      CALL mpp_max(zcorr_vol,npts)
+      CALL mpp_max(zcorr_sal,npts)
+      CALL mpp_max(zcorr_tem,npts)
+
       ! put correction term in the closest cell          
-      PRINT *, 'corrected point1 ', narea, zlon, zlat, izpts
       DO jpts = 1,npts
          CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
-         PRINT *, 'corrected point2 ', narea, jpts, ixpts(jpts), iypts(jpts), izpts(jpts)
          DO jj = mj0(iypts(jpts)),mj1(iypts(jpts))
             DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts))
                jk = izpts(jpts)
-               pvol_flx(ji,jj,jk)        =  pvol_flx(ji,jj,jk       ) + zcorr_vol(jpts)
-               pts_flx (ji,jj,jk,jp_sal) =  pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts)
-               pts_flx (ji,jj,jk,jp_tem) =  pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts)
-            END DO
-         END DO
-      END DO
+
+               IF (tmask_h(ji,jj) == 1._wp) THEN
+                  ! correct the vol_flx in the closest cell
+                  pvol_flx(ji,jj,jk)        =  pvol_flx(ji,jj,jk       ) + zcorr_vol(jpts)
+                  pts_flx (ji,jj,jk,jp_sal) =  pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts)
+                  pts_flx (ji,jj,jk,jp_tem) =  pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts)
+
+                  ! set correction to 0
+                  zcorr_vol(jpts) = 0.0_wp
+                  zcorr_sal(jpts) = 0.0_wp
+                  zcorr_tem(jpts) = 0.0_wp
+               END IF
+            END DO
+         END DO
+      END DO
+
       ! deallocate variables 
       DEALLOCATE(vnpts)
@@ -273,18 +277,10 @@
       pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
 
-      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1.)
-      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
-      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
-
-      ! CHECK vol !!!!!!!!! warning tmask_i wrong if deals with before value, so glob_sum wrong for before value!!!!
-      zsum  = glob_sum_full( pvol_flx(:,:,:) )
-      IF (lwp) PRINT *, 'CHECK vol = ',zsum
-      ! CHECK salt
-      zsum  = glob_sum( pts_flx(:,:,:,jp_sal) ) 
-      IF (lwp) PRINT *, 'CHECK salt = ',zsum
-      ! CHECK heat
-      zsum  = glob_sum( pts_flx(:,:,:,jp_tem) )
-      IF (lwp) PRINT *, 'CHECK heat = ',zsum
-   !! 
+      ! compute sum over the halo and set it to 0.
+      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1._wp)
+      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1._wp)
+      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1._wp)
+
+      ! deallocate variables
       CALL wrk_dealloc(jpi,jpj,jpk,   ztmp3d ) 
       CALL wrk_dealloc(jpi,jpj,       zde3t  ) 
Index: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90
===================================================================
--- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 5802)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 5820)
@@ -97,5 +97,4 @@
          CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) )
          CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) )
-         CALL iom_rstput( 0, 0, inum0, 'e3t_diff', (e3t_0(:,:,:)*tmask(:,:,:) - ze3t_b(:,:,:)*ztmask_b(:,:,:))*tmask(:,:,:) )
          CALL iom_close ( inum0 )
       END IF
@@ -191,9 +190,8 @@
                summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
                IF (zdsmask(ji,jj)==1 .AND. summsk .NE. 0) THEN
-                  PRINT *, 'add ssh to 1 cell',ji,jj,narea
-                  sshn(ji,jj)=(zssh0(jip1,jj)*zsmask0(jip1,jj)     &
-                  &           +zssh0(jim1,jj)*zsmask0(jim1,jj)     &
-                  &           +zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
-                  &           +zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
+                  sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
+                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     &
+                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
+                  &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
                   zsmask1(ji,jj)=1
                END IF
@@ -209,14 +207,18 @@
 !=============================================================================
       IF ( lk_vvl ) THEN
-! compute fse3t_n
-         DO jk = 1,jpk
-            fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) )
-         END DO
-
-! compute fse3u/v ... (call interpolation vvl)
-      ! Reconstruction of all vertical scale factors at now and before time steps
+      ! Reconstruction of all vertical scale factors at now time steps
       ! =============================================================================
       ! Horizontal scale factor interpolations
       ! --------------------------------------
+         DO jk = 1,jpk
+            DO jj=1,jpj
+               DO ji=1,jpi
+                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
+                     fse3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) / ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) )
+                  ENDIF
+               END DO
+            END DO
+         END DO
+
          CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
          CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
@@ -295,5 +297,5 @@
       zbvn(:,:)   = SUM(ztrp,DIM=3)
 
-! Already know ????????
+      ! new water column
       hu1=0.0_wp ;
       hv1=0.0_wp ;
@@ -302,4 +304,6 @@
         hv1(:,:) = hv1(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
       END DO
+      
+      ! compute correction      
       zucorr = 0._wp
       zvcorr = 0._wp
@@ -314,4 +318,6 @@
          END DO
       END DO 
+      
+      ! update velocity
       DO jk  = 1,jpk
          un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
@@ -333,15 +339,15 @@
                       jip1=ji+1; jim1=ji-1;
                       jjp1=jj+1; jjm1=jj-1;
-                      summsk=(ztmask0(jip1,jj,jk)+ztmask0(jim1,jj,jk)+ztmask0(ji,jjp1,jk)+ztmask0(ji,jjm1,jk))
-                      !! horizontal extrapolation
+                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
+                      !! horizontal basic extrapolation
                       IF (zdmask(ji,jj)==1 .AND. summsk .NE. 0) THEN
-                         tsn(ji,jj,jk,1)=( zts0(jip1,jj,jk,1)*ztmask0(jip1,jj,jk)     &
-                         &                +zts0(jim1,jj,jk,1)*ztmask0(jim1,jj,jk)     &
-                         &                +zts0(ji,jjp1,jk,1)*ztmask0(ji,jjp1,jk)     &
-                         &                +zts0(ji,jjm1,jk,1)*ztmask0(ji,jjm1,jk))/summsk
-                         tsn(ji,jj,jk,2)=( zts0(jip1,jj,jk,2)*ztmask0(jip1,jj,jk)     &
-                         &                +zts0(jim1,jj,jk,2)*ztmask0(jim1,jj,jk)     &
-                         &                +zts0(ji,jjp1,jk,2)*ztmask0(ji,jjp1,jk)     &
-                         &                +zts0(ji,jjm1,jk,2)*ztmask0(ji,jjm1,jk))/summsk
+                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
+                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
+                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
+                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
+                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
+                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
+                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
+                         &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
                          ztmask1(ji,jj,jk)=1
                       END IF
@@ -358,21 +364,4 @@
                          END IF
                       END IF
-                      !! horizontal corner extrapolation if horizontal and vertical failed
-                      IF (zdmask(ji,jj)==1 .AND. summsk==0) THEN
-                         jip1=ji+1; jim1=ji-1;
-                         jjp1=jj+1; jjm1=jj-1;
-                         summsk=(ztmask0(jip1,jjp1,jk)+ztmask0(jim1,jjm1,jk)+ztmask0(jim1,jjp1,jk)+ztmask0(jip1,jjm1,jk))
-                         IF (zdmask(ji,jj)==1 .AND. summsk .NE. 0) THEN
-                            tsn(ji,jj,jk,1)=( zts0(jip1,jjp1,jk,1)*ztmask0(jip1,jjp1,jk)     &
-                            &                +zts0(jim1,jjm1,jk,1)*ztmask0(jim1,jjm1,jk)     &
-                            &                +zts0(jim1,jjp1,jk,1)*ztmask0(jim1,jjp1,jk)     &
-                            &                +zts0(jip1,jjm1,jk,1)*ztmask0(jip1,jjm1,jk))/summsk
-                            tsn(ji,jj,jk,2)=( zts0(jip1,jjp1,jk,2)*ztmask0(jip1,jjp1,jk)     &
-                            &                +zts0(jim1,jjm1,jk,2)*ztmask0(jim1,jjm1,jk)     &
-                            &                +zts0(jim1,jjp1,jk,2)*ztmask0(jim1,jjp1,jk)     &
-                            &                +zts0(jip1,jjm1,jk,2)*ztmask0(jip1,jjm1,jk))/summsk
-                            ztmask1(ji,jj,jk)=1
-                         END IF
-                      END IF
                   END DO
               END DO
@@ -392,5 +381,5 @@
             DO jj = 1,jpj
                DO ji = 1,jpi
-                  IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1.0_wp ) THEN
+                  IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1.0_wp .AND. (tmask(ji,jj,1)==0 .OR. ptmask_b(ji,jj,1)==0) ) THEN
                      !compute weight
                      zdzp1 = MAX(0._wp,fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
@@ -398,15 +387,15 @@
                      zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - fsdepw_n(ji,jj,jk  ))
                      IF (zdz .LT. 0.0_wp) THEN 
-                        PRINT *, 'ERROR dz n ', ji,jj,jk,zdz,fsdepw_n(ji,jj,jk+1),fsdepw_n(ji,jj,jk),fsdepw_n(ji,jj,jk-1)
-                        PRINT *, 'ERROR dz n             = ',fse3t_n (ji,jj,jk+1),fse3t_n (ji,jj,jk),fse3t_n (ji,jj,jk-1), sshn(ji,jj)
-                        PRINT *, 'ERROR dz b             = ',pdepw_b(ji,jj,jk+1),pdepw_b(ji,jj,jk),pdepw_b(ji,jj,jk-1)
-                        PRINT *, 'ERROR dz b             = ',fse3t_b (ji,jj,jk+1),fse3t_b (ji,jj,jk),fse3t_b (ji,jj,jk-1), sshb(ji,jj)
-                        PRINT *, 'ERROR dz 0             = ',  e3t_0 (ji,jj,jk+1),  e3t_0 (ji,jj,jk),  e3t_0 (ji,jj,jk-1)
-                        PRINT *, 'ERROR dz n             = ',  tmask (ji,jj,jk+1),  tmask (ji,jj,jk),  tmask (ji,jj,jk-1)
-                        PRINT *, 'ERROR dz n             = ', zwmaskn(ji,jj,jk+1), zwmaskn(ji,jj,jk), zwmaskn(ji,jj,jk-1)
-                        PRINT *, 'ERROR dz b             = ', ptmask_b(ji,jj,jk+1), ptmask_b(ji,jj,jk), ptmask_b(ji,jj,jk-1)
-                        PRINT *, 'ERROR dz b             = ', zwmaskb(ji,jj,jk+1), zwmaskb(ji,jj,jk), zwmaskb(ji,jj,jk-1)
-                        PRINT *, 'ERROR dz b             = ', gdepw_0(ji,jj,jk+1), gdepw_0(ji,jj,jk), gdepw_0(ji,jj,jk-1)
-                        STOP
+                        WRITE(numout,*) 'ERROR dz n ', ji,jj,jk,zdz,fsdepw_n(ji,jj,jk+1),fsdepw_n(ji,jj,jk),fsdepw_n(ji,jj,jk-1)
+                        WRITE(numout,*) 'ERROR dz n             = ',fse3t_n (ji,jj,jk+1),fse3t_n (ji,jj,jk),fse3t_n (ji,jj,jk-1), sshn(ji,jj)
+                        WRITE(numout,*) 'ERROR dz b             = ',pdepw_b(ji,jj,jk+1),pdepw_b(ji,jj,jk),pdepw_b(ji,jj,jk-1)
+                        WRITE(numout,*) 'ERROR dz b             = ',fse3t_b (ji,jj,jk+1),fse3t_b (ji,jj,jk),fse3t_b (ji,jj,jk-1), sshb(ji,jj)
+                        WRITE(numout,*) 'ERROR dz 0             = ',  e3t_0 (ji,jj,jk+1),  e3t_0 (ji,jj,jk),  e3t_0 (ji,jj,jk-1)
+                        WRITE(numout,*) 'ERROR dz n             = ',  tmask (ji,jj,jk+1),  tmask (ji,jj,jk),  tmask (ji,jj,jk-1)
+                        WRITE(numout,*) 'ERROR dz n             = ', zwmaskn(ji,jj,jk+1), zwmaskn(ji,jj,jk), zwmaskn(ji,jj,jk-1)
+                        WRITE(numout,*) 'ERROR dz b             = ', ptmask_b(ji,jj,jk+1), ptmask_b(ji,jj,jk), ptmask_b(ji,jj,jk-1)
+                        WRITE(numout,*) 'ERROR dz b             = ', zwmaskb(ji,jj,jk+1), zwmaskb(ji,jj,jk), zwmaskb(ji,jj,jk-1)
+                        WRITE(numout,*) 'ERROR dz b             = ', gdepw_0(ji,jj,jk+1), gdepw_0(ji,jj,jk), gdepw_0(ji,jj,jk-1)
+                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
                      END IF
                      tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
