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 5790)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 5802)
@@ -200,5 +200,11 @@
 !      ENDIF
 !!gm end
-      
+      IF ( lk_vvl ) THEN
+         IF (lwp) PRINT *, 'cons heat : ', kt, zdiff_hc / zvol_tot, zdiff_sc / zvol_tot 
+         IF (lwp) PRINT *, 'cons volu : ', kt, zdiff_v2 * 1.e-9      
+      ELSE
+         IF (lwp) PRINT *, 'cons heat : ', kt, zdiff_hc1 * 1.e-20 * rau0 * rcp, zdiff_sc1 * 1.e-9
+         IF (lwp) PRINT *, 'cons vol  : ', kt, zdiff_v1 * 1.e-9      
+      END IF
       IF( lk_vvl ) THEN
         CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature variation (C) 
Index: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 5790)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 5802)
@@ -253,6 +253,7 @@
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt               !: vertical index of the bottom last T- ocean level
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbku, mbkv         !: vertical index of the bottom last U- and W- ocean level
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy                     !: ocean depth (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i                   !: interior domain T-point mask
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i            !: interior domain T-point mask
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_h            !: internal domain T-point mask (Figure 8.5 NEMO book)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask                     !: land/ocean mask of barotropic stream function
 
@@ -391,5 +392,5 @@
 
       ALLOCATE( mbathy(jpi,jpj) , bathy  (jpi,jpj) ,                                       &
-         &     tmask_i(jpi,jpj) ,                                                          & 
+         &     tmask_i(jpi,jpj) , tmask_h(jpi, jpj),                                       & 
          &     ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , &
          &     bmask(jpi,jpj)   ,                                                          &
Index: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 5790)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 5802)
@@ -222,4 +222,6 @@
       ! --------------------
       tmask_i(:,:) = ssmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf
+
+      tmask_h(:,:) = 1._wp                 ! 0 on the halo and 1 elsewhere
       iif = jpreci                         ! ???
       iil = nlci - jpreci + 1
@@ -227,8 +229,8 @@
       ijl = nlcj - jprecj + 1
 
-      tmask_i( 1 :iif,   :   ) = 0._wp      ! first columns
-      tmask_i(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
-      tmask_i(   :   , 1 :ijf) = 0._wp      ! first rows
-      tmask_i(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
+      tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns
+      tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
+      tmask_h(   :   , 1 :ijf) = 0._wp      ! first rows
+      tmask_h(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
 
       ! north fold mask
@@ -241,8 +243,11 @@
          IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row
             DO ji = iif+1, iil-1
-               tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
+               tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji))
             END DO
          ENDIF
       ENDIF
+     
+      tmask_i(:,:) = tmask_i(:,:) * tmask_h(:,:)
+
       IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
          tpol(     1    :jpiglo) = 0._wp
Index: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5790)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5802)
@@ -1286,11 +1286,14 @@
 
          ! split last cell if possible (only where water column is 2 cell or less)
-         !DO jk = jpkm1, 1, -1
-         !   zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
-         !   WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
-         !      mbathy(:,:) = jk
-         !      bathy(:,:)  = zmax
-         !   END WHERE
-         !END DO
+         ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss).
+         IF ( .NOT. ln_iscpl) THEN
+            DO jk = jpkm1, 1, -1
+               zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
+               WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
+                  mbathy(:,:) = jk
+                  bathy(:,:)  = zmax
+               END WHERE
+            END DO
+         END IF
  
          ! split top cell if possible (only where water column is 2 cell or less)
@@ -1310,15 +1313,24 @@
                ! test bathy
                IF (risfdep(ji,jj) .GT. 1) THEN
-               zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
-               zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
+                  IF ( .NOT. ln_iscpl ) THEN
+                     zbathydiff  =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) &
+                         &            + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) &
+                         &            - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
  
-                  IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
-!                     IF (zbathydiff .LE. zrisfdepdiff) THEN
-!                        bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
-!                        mbathy(ji,jj)= mbathy(ji,jj) + 1
-!                     ELSE
+                     IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
+                        IF (zbathydiff .LE. zrisfdepdiff) THEN
+                           bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
+                           mbathy(ji,jj)= mbathy(ji,jj) + 1
+                        ELSE
+                           risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
+                           misfdep(ji,jj) = misfdep(ji,jj) - 1
+                        END IF
+                     ENDIF
+                  ELSE
+                     IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
                         risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
                         misfdep(ji,jj) = misfdep(ji,jj) - 1
-!                     END IF
+                     END IF
                   END IF
                END IF
@@ -1332,13 +1344,20 @@
                ! test bathy
                IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
-!                  zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
-!                  zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
-!                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-!                     mbathy(ji,jj) = mbathy(ji,jj) + 1
-!                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
-!                  ELSE
+                  IF ( .NOT. ln_iscpl ) THEN 
+                     zbathydiff  =ABS(bathy(ji,jj)   - ( gdepw_1d(mbathy (ji,jj)+1) &
+                         &                             + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj)  ) & 
+                         &                             - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
+                     IF (zbathydiff .LE. zrisfdepdiff) THEN
+                        mbathy(ji,jj) = mbathy(ji,jj) + 1
+                        bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
+                     ELSE
+                        misfdep(ji,jj)= misfdep(ji,jj) - 1
+                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
+                     END IF
+                  ELSE
                      misfdep(ji,jj)= misfdep(ji,jj) - 1
-                     risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
-!                  END IF
+                     risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
+                  END IF
                ENDIF
             END DO
@@ -1349,13 +1368,20 @@
             DO ji = 1, jpim1
                IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
-!                  zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat )))
-!                  zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
-!                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-!                     mbathy(ji,jj) = mbathy(ji,jj) + 1
-!                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat )
-!                  ELSE
+                  IF ( .NOT. ln_iscpl ) THEN 
+                     zbathydiff  =ABS(bathy(ji,jj    ) - ( gdepw_1d(mbathy (ji,jj)+1) &
+                          &                              + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) &
+                          &                              - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
+                     IF (zbathydiff .LE. zrisfdepdiff) THEN
+                        mbathy(ji,jj) = mbathy(ji,jj) + 1
+                        bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat )
+                     ELSE
+                        misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1
+                        risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
+                     END IF
+                  ELSE
                      misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1
                      risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
-!                  END IF
+                  END IF
                ENDIF
             END DO
@@ -1376,13 +1402,20 @@
             DO ji = 1, jpim1
                IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN
-!                  zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
-!                  zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat )))
-!                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-!                     mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
-!                     bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
-!                  ELSE
+                  IF ( .NOT. ln_iscpl ) THEN 
+                     zbathydiff  =ABS(  bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) &
+                           &                             + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji,jj  ) - ( gdepw_1d(misfdep(ji,jj  )  ) &
+                           &                             - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat )))
+                     IF (zbathydiff .LE. zrisfdepdiff) THEN
+                        mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
+                        bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
+                     ELSE
+                        misfdep(ji,jj)   = misfdep(ji,jj) - 1
+                        risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat )
+                     END IF
+                  ELSE
                      misfdep(ji,jj)   = misfdep(ji,jj) - 1
                      risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat )
-!                  END IF
+                  END IF
                ENDIF
             END DO
@@ -1405,13 +1438,20 @@
             DO ji = 1, jpim1
                IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
-!                  zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat )))
-!                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
-!                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-!                     mbathy(ji,jj) = mbathy(ji,jj) + 1
-!                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
-!                  ELSE
+                  IF ( .NOT. ln_iscpl ) THEN 
+                  zbathydiff  =ABS(  bathy(ji  ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) &
+                       &                              + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat )))
+                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) &
+                       &                              - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
+                  IF (zbathydiff .LE. zrisfdepdiff) THEN
+                     mbathy(ji,jj) = mbathy(ji,jj) + 1
+                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
+                  ELSE
                      misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
                      risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
-!                  END IF
+                  END IF
+                  ELSE
+                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
+                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
+                  ENDIF
                ENDIF
             ENDDO
@@ -1433,13 +1473,20 @@
             DO ji = 1, jpim1
                IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
-!                  zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
-!                  zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  ) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat )))
-!                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-!                     mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1
-!                     bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
-!                  ELSE
+                  IF ( .NOT. ln_iscpl ) THEN 
+                     zbathydiff  =ABS(  bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) &
+                          &                              + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji  ,jj) - ( gdepw_1d(misfdep(ji  ,jj)  ) &
+                          &                              - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat )))
+                     IF (zbathydiff .LE. zrisfdepdiff) THEN
+                        mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1
+                        bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
+                     ELSE
+                        misfdep(ji,jj)   = misfdep(ji  ,jj) - 1
+                        risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat )
+                     END IF
+                  ELSE
                      misfdep(ji,jj)   = misfdep(ji  ,jj) - 1
                      risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat )
-!                  END IF
+                  ENDIF
                ENDIF
             ENDDO
@@ -1530,8 +1577,8 @@
                ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  )
                ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1)
+               IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk
+               IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk
+               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk
+               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk
                ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
                IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN
@@ -1562,5 +1609,5 @@
                ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  )
                ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1)
-               IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1)
+               IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0
                IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj  ) ) ibtestip1 = 0
                IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj-1) ) ibtestjm1 = 0
@@ -1756,5 +1803,5 @@
                ENDIF 
             !       ... on ik / ik-1 
-               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
+               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) 
                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 
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 5790)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90	(revision 5802)
@@ -64,5 +64,5 @@
       REAL(wp):: zjip1_ratio, zjim1_ratio, zjjp1_ratio, zjjm1_ratio
       !!
-      REAL(wp), DIMENSION(:,:    ), POINTER :: zde3t
+      REAL(wp), DIMENSION(:,:    ), POINTER :: zde3t, zdtem, zdsal
       REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  
       REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmp3d
@@ -74,246 +74,220 @@
 
       CALL wrk_alloc(jpi,jpj,jpk,   ztmp3d )
-      CALL wrk_alloc(jpi,jpj,       zde3t  )
+      CALL wrk_alloc(jpi,jpj,       zde3t , zdtem, zdsal )
       CALL wrk_alloc(jpi,jpj,       zssh0  )
 
     ! get unbalance (volume heat and salt)
     ! initialisation
-       zde3t   (:,:)     = 0.0_wp
-       pvol_flx(:,:,:  ) = 0.0_wp
-       pts_flx (:,:,:,:) = 0.0_wp
-
-       zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total volume correction 0 = ',zsum
-       zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total heat correction 0 = ',zsum
-       zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total salt correction 0 = ',zsum
-
-       ! mask tsn and tsb (should be useless)
-       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(:,:,:);
-
-       ! diagnose non conservation of heat, salt and volume 
-       r1_tiscpl = 1._wp / (prdt_iscpl * rn_rdt) 
-       zssh0(:,:)        = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
-       IF ( lk_vvl ) zssh0 = 0.0_wp
-       DO jk = 1,jpk-1
-          DO ji = 2,jpi-1
-             DO jj = 2,jpj-1
-                ! volume differences
-                zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - 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
-
-                ! ocean cell now
-                ! case where we open, enlarge or thin a cell :
-                pvol_flx(ji,jj,jk)       =                          zde3t(ji,jj) * r1_tiscpl
-                pts_flx (ji,jj,jk,jp_sal)=   tsn(ji,jj,jk,jp_sal) * zde3t(ji,jj) * r1_tiscpl 
-                pts_flx (ji,jj,jk,jp_tem)=   tsn(ji,jj,jk,jp_tem) * zde3t(ji,jj) * r1_tiscpl
-             END DO
-          END DO
-       END DO
-       ! glob_sum_full because with glob summ some data can be masked. WARNING the halo have to be set at 0
-       PRINT *, 'test ', narea, SUM(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt, SUM(pvol_flx(2:jpi-1,2:jpj-1,:)) * rn_fiscpl * rn_rdt
-       zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total volume correction 1 = ',zsum
-       zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total heat correction 1 = ',zsum
-       zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total salt correction 1 = ',zsum
-
-       zssh0(:,:)        = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
-       IF ( lk_vvl ) zssh0 = 0.0_wp
-       DO jk = 1,jpk-1
-          DO ji = 2,jpi-1
-             DO jj = 2,jpj-1
-                ! volume differences
-                zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - 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
-
-                ! ocean cell before and mask cell now
-                IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
-                   ! case where we close a cell and adjacent cell open
-                   pvol_flx(ji,jj,jk)       =                         zde3t(ji,jj) * r1_tiscpl
-                   pts_flx (ji,jj,jk,jp_sal)=  tsb(ji,jj,jk,jp_sal) * zde3t(ji,jj) * r1_tiscpl 
-                   pts_flx (ji,jj,jk,jp_tem)=  tsb(ji,jj,jk,jp_tem) * zde3t(ji,jj) * r1_tiscpl
-
-                   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)
-
-                         ! 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
-                         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.
-                      END IF
-                   END IF
-                END IF
-             END DO
-          END DO
-       END DO
-
-       zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total volume correction 2 = ',zsum
-       zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total heat correction 2 = ',zsum
-       zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total salt correction 2 = ',zsum
-
-       ! allocation and initialisation of the list of problematic point
-       ALLOCATE(vnpts(jpnij))
-       vnpts(:)=0
-
-       ! fill narea location with the number of problematic point
-       DO jk = 1,jpk-1
-          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 ) THEN
-                   vnpts(narea) = vnpts(narea) + 1 
-                END IF
-             END DO
-          END DO
-       END DO
-
-       ! build array of total problematic point on each cpu (share to each cpu)
-       CALL mpp_max(vnpts,jpnij) 
-
-       ! size of the new variable
-       npts  = SUM(vnpts)    
-       
-       ! allocation of the coordinates, correction, index vector for the problematic points 
-       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
-
-       ! fill new variable
-       jpts = SUM(vnpts(1:narea-1))
-       DO jk = 1,jpk-1
-          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 ) 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)
-                   zcorr_vol(jpts) = pvol_flx(ji,jj,jk)
-                   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
-          END DO
-       END DO
-
-       ! build array of total problematic point on each cpu (share to each cpu)
-       CALL mpp_max(zlat ,npts)
-       CALL mpp_max(zlon ,npts)
-       CALL mpp_max(izpts,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
-       ! deallocate variables 
-       DEALLOCATE(vnpts)
-       DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
-     
-       ! add contribution store on the hallo (lbclnk remove one of the contribution)
-       pvol_flx(:,:,:       ) = pvol_flx(:,:,:       ) * tmask(:,:,:)
-       pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
-       pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * 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!!!!
-       zsumn = glob_sum     ( fse3t_n(:,:,:) * tmask  (:,:,:)) - glob_sum(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt
-       ztmp3d(:,:,:) = 0.0
-       ztmp3d(2:jpi-1,2:jpj-1,:) = pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
-       zsumb = glob_sum_full(ztmp3d) 
-       zsum  = glob_sum     ( pvol_flx(:,:,:) * rn_fiscpl * rn_rdt)
-       IF (lwp) PRINT *, 'CHECK vol = ',zsumn, zsumb, zsumn - zsumb, zsum
-       ! CHECK salt
-       zsumn = glob_sum( tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) *  tmask  (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
-       zsumb = glob_sum( tsb(:,:,:,jp_sal) *  pe3t_b(:,:,:) * ptmask_b(:,:,:))
-       zsum  = glob_sum( pts_flx(:,:,:,jp_sal)*rn_fiscpl * rn_rdt)
-       IF (lwp) PRINT *, 'CHECK salt = ',zsumn, zsumb, zsumn - zsumb, zsum
-       ! CHECK heat
-       zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) *  tmask  (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
-       zsumb = glob_sum( tsb(:,:,:,jp_tem) *  pe3t_b(:,:,:) * ptmask_b(:,:,:))
-       zsum  = glob_sum( pts_flx(:,:,:,jp_tem)*rn_fiscpl * rn_rdt)
-       IF (lwp) PRINT *, 'CHECK heat = ',zsumn, zsumb, zsumn - zsumb, zsum
-    !! 
-    CALL wrk_dealloc(jpi,jpj,jpk,   ztmp3d ) 
-    CALL wrk_dealloc(jpi,jpj,       zde3t  ) 
-    CALL wrk_dealloc(jpi,jpj,       zssh0  ) 
+      zde3t   (:,:)     = 0.0_wp
+      pvol_flx(:,:,:  ) = 0.0_wp
+      pts_flx (:,:,:,:) = 0.0_wp
+
+      ! mask tsn and tsb (should be useless)
+      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(:,:,:);
+
+      ! diagnose non conservation of heat, salt and volume 
+      r1_tiscpl = 1._wp / (prdt_iscpl * rn_rdt) 
+
+      zssh0(:,:)        = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
+      IF ( lk_vvl ) zssh0 = 0.0_wp ! already include in the levels by definition
+
+      DO jk = 1,jpk-1
+         DO ji = 2,jpi-1
+            DO jj = 2,jpj-1
+               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)   &
+                            - 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)
+               
+               ! 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)
+
+                        ! 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
+                        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.
+                     END IF
+                  END IF
+               END IF
+               END IF
+            END DO
+         END DO
+      END DO
+
+      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.)
+
+      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 
+      ! allocation and initialisation of the list of problematic point
+      ALLOCATE(vnpts(jpnij))
+      vnpts(:)=0
+
+      ! fill narea location with the number of problematic point
+      DO jk = 1,jpk-1
+         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
+                  vnpts(narea) = vnpts(narea) + 1 
+               END IF
+            END DO
+         END DO
+      END DO
+
+      ! build array of total problematic point on each cpu (share to each cpu)
+      CALL mpp_max(vnpts,jpnij) 
+
+      ! size of the new variable
+      npts  = SUM(vnpts)    
+      
+      ! allocation of the coordinates, correction, index vector for the problematic points 
+      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
+
+      ! fill new variable
+      jpts = SUM(vnpts(1:narea-1))
+      DO jk = 1,jpk-1
+         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
+                  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)
+                  zcorr_vol(jpts) = pvol_flx(ji,jj,jk)
+                  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
+         END DO
+      END DO
+
+      ! build array of total problematic point on each cpu (share to each cpu)
+      CALL mpp_max(zlat ,npts)
+      CALL mpp_max(zlon ,npts)
+      CALL mpp_max(izpts,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
+      ! deallocate variables 
+      DEALLOCATE(vnpts)
+      DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
+    
+      ! add contribution store on the hallo (lbclnk remove one of the contribution)
+      pvol_flx(:,:,:       ) = pvol_flx(:,:,:       ) * tmask(:,:,:)
+      pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
+      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
+   !! 
+      CALL wrk_dealloc(jpi,jpj,jpk,   ztmp3d ) 
+      CALL wrk_dealloc(jpi,jpj,       zde3t  ) 
+      CALL wrk_dealloc(jpi,jpj,       zssh0  ) 
    END SUBROUTINE iscpl_cons
 
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 5790)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 5802)
@@ -59,6 +59,6 @@
       CHARACTER(20) :: cfile
 
-      CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! need this new variables for evoving isf cavity
-      CALL wrk_alloc(jpi,jpj,jpk, ze3t_b  , ze3u_b  , ze3v_b  ) ! need of this variables for interpolation)
+      CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! mask before
+      CALL wrk_alloc(jpi,jpj,jpk, ze3t_b  , ze3u_b  , ze3v_b  ) ! e3   before
       CALL wrk_alloc(jpi,jpj,jpk, zdepw_b )
       CALL wrk_alloc(jpi,jpj,     zsmask_b                    )
@@ -105,14 +105,26 @@
       CALL wrk_dealloc(jpi,jpj,jpk, zdepw_b                    )
       CALL wrk_dealloc(jpi,jpj,     zsmask_b                   )
- 
-      neuler = 0              ! Euler restart (neuler=0)
-
-      tsb    (:,:,:,:) = tsn    (:,:,:,:)  ! all before fields set to now values
-      ub     (:,:,:)   = un     (:,:,:)
-      vb     (:,:,:)   = vn     (:,:,:)
-      rotb   (:,:,:)   = rotn   (:,:,:)
-      hdivb  (:,:,:)   = hdivn  (:,:,:)
-      fse3t_b(:,:,:)   = fse3t_n(:,:,:)
-      sshb   (:,:)     = sshn   (:,:)
+
+      !! next step is an euler time step
+      neuler = 0
+      !! set _b and _n variables equal
+      tsb (:,:,:,:) = tsn (:,:,:,:)
+      ub  (:,:,:  ) = un  (:,:,:  )
+      vb  (:,:,:  ) = vn  (:,:,:  )
+      sshb(:,:    ) = sshn(:,:)
+      !! set _b and _n vertical scale factor equal
+      fse3t_b (:,:,:) = fse3t_n (:,:,:)
+      fse3u_b (:,:,:) = fse3u_n (:,:,:)
+      fse3v_b (:,:,:) = fse3v_n (:,:,:)
+      IF ( lk_vvl ) THEN
+         fse3uw_b(:,:,:) = fse3uw_n(:,:,:)
+         fse3vw_b(:,:,:) = fse3vw_n(:,:,:)
+         fsdept_b(:,:,:) = fsdept_n(:,:,:)
+         fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
+         hu_b (:,:) = hu(:,:)
+         hv_b (:,:) = hv(:,:)
+         hur_b(:,:) = hur(:,:)
+         hvr_b(:,:) = hvr(:,:)
+      END IF
       !
    END SUBROUTINE iscpl_stp
@@ -166,11 +178,9 @@
       END DO
            
-      ! compute ssh     
-      zsmask0(:,:) = psmask_b(:,:)
-      zsmask1(:,:) = psmask_b(:,:)
-      
       ! compute new ssh if we open a full water column (average of the closest neigbourgs)  
       sshb (:,:)=sshn(:,:)
       zssh0(:,:)=sshn(:,:)
+      zsmask0(:,:) = psmask_b(:,:)
+      zsmask1(:,:) = psmask_b(:,:)
       DO iz = 1,10    ! need to be tuned (configuration dependent)
          zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
@@ -198,20 +208,9 @@
 
 !=============================================================================
-! WARNING we cannot used glob_sum for before time variable (because glob_sum use tmask_i). Need to mask the halo and glob_sum_full
-      ztmp3d(:,:,:) = 0.0
-      ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
-      &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
-      zsum = glob_sum_full(ztmp3d) 
-      IF (lwp) PRINT *, 'total volume correction 00 = ',zsum
       IF ( lk_vvl ) THEN
 ! compute fse3t_n
          DO jk = 1,jpk
-            fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1 + sshn(:,:) / ( ht_0(:,:) + 1. - ssmask(:,:) ) * tmask(:,:,jk) )
-         END DO
-         ztmp3d(:,:,:) = 0.0
-         ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
-         &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
-         zsum = glob_sum_full(ztmp3d) 
-         IF (lwp) PRINT *, 'total volume correction 01 = ',zsum
+            fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) )
+         END DO
 
 ! compute fse3u/v ... (call interpolation vvl)
@@ -269,9 +268,4 @@
 
       END IF
-      ztmp3d(:,:,:) = 0.0
-      ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
-      &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
-      zsum = glob_sum_full(ztmp3d) 
-      IF (lwp) PRINT *, 'total volume correction 02 = ',zsum
 
 !=============================================================================
@@ -289,4 +283,5 @@
          END DO
       END DO
+
 ! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
 ! compute barotropic velocity now and after 
@@ -299,4 +294,5 @@
       ztrp(:,:,:) = vn(:,:,:)*fse3v_n(:,:,:); 
       zbvn(:,:)   = SUM(ztrp,DIM=3)
+
 ! Already know ????????
       hu1=0.0_wp ;
@@ -391,12 +387,6 @@
       tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
 
-      ! CHECK heat
-      zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) * zwmaskb(:,:,:) * zwmaskn(:,:,:))
-      zsumb = glob_sum( tsb(:,:,:,jp_tem) *  pe3t_b(:,:,:) * zwmaskb(:,:,:) * zwmaskn(:,:,:))
-      IF (lwp) PRINT *, 'CHECK tsn = ',zsumn, zsumb
-
-      ! compute new T/S (interpolation) because of vvl
-      IF ( lk_vvl .AND. .FALSE. ) THEN
-      !IF ( lk_vvl ) THEN
+      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
+      IF ( lk_vvl ) THEN
          DO jk = 2,jpk-1
             DO jj = 1,jpj
@@ -432,11 +422,5 @@
       END IF
 
-      ! CHECK heat
-      ztmp3d(:,:,:) = 0.0
-      ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
-      &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
-      zsum = glob_sum_full(ztmp3d) 
-      IF (lwp) PRINT *, 'total volume correction 03 = ',zsum
-
+      ! Special value for closed pool and set the mask to 0 to run
       WHERE (tmask(:,:,:) == 1.0 .AND. tsn(:,:,:,2) == 0._wp) 
          tsn(:,:,:,2)=  -99._wp
@@ -453,27 +437,6 @@
       ! compute new tn and sn if we close cell 
       ! nothing to do
-
-      !! next step is an euler time step
-      neuler = 0 
-      !! set _b and _n variables equal
-      tsb (:,:,:,:)=tsn (:,:,:,:)
-      ub  (:,:,:  )=un  (:,:,:  )
-      vb  (:,:,:  )=vn  (:,:,:  )
-      sshb(:,:    )=sshn(:,:)
-      !! set _b and _n vertical scale factor equal
-      fse3t_b (:,:,:)=fse3t_n (:,:,:)
-      fse3u_b (:,:,:)=fse3u_n (:,:,:)
-      fse3v_b (:,:,:)=fse3v_n (:,:,:)
-      IF ( lk_vvl ) THEN
-         fse3uw_b(:,:,:) = fse3uw_n(:,:,:)
-         fse3vw_b(:,:,:) = fse3vw_n(:,:,:)
-         fsdept_b(:,:,:) = fsdept_n(:,:,:)
-         fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
-         hu_b(:,:) = hu(:,:)
-         hv_b(:,:) = hv(:,:)
-         hur_b(:,:) = hur(:,:)
-         hvr_b(:,:) = hvr(:,:)
-      END IF
-      !! 
+      ! 
+      ! deallocation tmp arrays
       CALL wrk_dealloc(jpi,jpj,jpk,2, zts0                                   )
       CALL wrk_dealloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp                ) 
Index: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90
===================================================================
--- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90	(revision 5790)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90	(revision 5802)
@@ -172,5 +172,5 @@
       !!-----------------------------------------------------------------------
       !
-      glob_sum_full_2d = SUM( ptab(:,:) )
+      glob_sum_full_2d = SUM( ptab(:,:) * tmask_h(:,:) )
       IF( lk_mpp )   CALL mpp_sum( glob_sum_full_2d )
       !
@@ -194,5 +194,5 @@
       glob_sum_full_3d = 0.e0
       DO jk = 1, ijpk
-         glob_sum_full_3d = glob_sum_full_3d + SUM( ptab(:,:,jk) )
+         glob_sum_full_3d = glob_sum_full_3d + SUM( ptab(:,:,jk) * tmask_h(:,:) )
       END DO
       IF( lk_mpp )   CALL mpp_sum( glob_sum_full_3d )
@@ -376,5 +376,5 @@
       DO jj = 1, jpj
          DO ji =1, jpi
-         ztmp =  ptab(ji,jj)
+         ztmp =  ptab(ji,jj) * tmask_h(ji,jj) 
          CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
          END DO
@@ -407,5 +407,5 @@
          DO jj = 1, jpj
             DO ji =1, jpi
-            ztmp =  ptab(ji,jj,jk)
+            ztmp =  ptab(ji,jj,jk) * tmask_h(ji,jj)
             CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
             END DO
