Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscpldiv.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscpldiv.F90	(revision 9813)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscpldiv.F90	(revision 9813)
@@ -0,0 +1,188 @@
+MODULE iscpldiv
+   !!======================================================================
+   !!                       ***  MODULE  iscpldiv***
+   !! Divergence correction : compute and save the divergence correction due to ice sheet/ocean coupling
+   !!=====================================================================
+   !! History :  NEMO  ! 2018-06 A. Siahaan: original 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   iscpl_div_corr     : compute the divergence correction 
+   !!   iscpl_div          : update hdivn with the divergence correction
+   !!   iscpl_div_tra      : update tsa with contribution from the divergence correction
+   !!----------------------------------------------------------------------
+   USE dom_oce         ! ocean space and time domain
+   USE oce             ! global tra/dyn variable
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! MPP library
+   USE lib_fortran     ! MPP library
+   USE wrk_nemo        ! Memory allocation
+   USE lbclnk          !
+   USE iscplini
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC iscpl_div
+   PUBLIC iscpl_div_tra
+   PUBLIC iscpl_div_corr
+   !!
+   !! * Substitutions  
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE iscpl_div( phdivn )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE iscpl_div  ***
+      !!
+      !! ** Purpose :   update the horizontal divergence
+      !!
+      !! ** Method  :
+      !!                CAUTION : iscpl is positive (inflow) and expressed in
+      !m/s
+      !!
+      !! ** Action  :   phdivn   increase by the iscpl correction term
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
+      !!
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*) 'adding the divergence correction contribution at nit000'
+
+      phdivn(:,:,:) =  phdivn(:,:,:) + rhdivdiff(:,:,:)
+
+   END SUBROUTINE iscpl_div
+
+   SUBROUTINE iscpl_div_tra( ptsa )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE iscpl_div_tra  ***
+      !!
+      !! ** Purpose :   update the active tracer tendency
+      !!
+      !! ** Method  :
+      !!                CAUTION : iscpl is positive (inflow) and expressed in
+      !m/s
+      !!
+      !! ** Action  :   ptsa   increase by the iscpl correction term
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   ptsa   ! active tracer tendency
+      !!
+      !!----------------------------------------------------------------------
+      !
+      IF (lwp) WRITE(numout,*) 'add tendency tsa due to re-meshing at kt = nit000'
+
+      ptsa(:,:,:,1) = ptsa(:,:,:,1) + rhdivdiff_trac(:,:,:,1)*tmask(:,:,:)
+      ptsa(:,:,:,2) = ptsa(:,:,:,2) + rhdivdiff_trac(:,:,:,2)*tmask(:,:,:)
+
+   END SUBROUTINE iscpl_div_tra
+
+
+
+   SUBROUTINE iscpl_div_corr(ptmask_b, pe3t_b, pe3u_b, pe3v_b)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE iscpl_div_corr  ***
+      !!
+      !! ** Purpose :   compute the divergence correction
+      !!
+      !! ** Method  :
+      !!               compute the difference of hdivn*e3t_n between before and after remeshing    
+      !!               transfer the missing hdvn*e3t_n from top of isf to the new mikt
+      !!               transfer the missing hdvn*e3t_n from bottom of isf to the new mbkt 
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b !! mask before
+      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
+      !!
+      INTEGER :: ji, jj, jk, jtop         !! loop index
+      !!
+      REAL(wp):: zvol, zuflux_sum, zvflux_sum
+
+      REAL(wp), DIMENSION(:,:,:  ), POINTER :: zhdiv, zhdiv2
+      !!----------------------------------------------------------------------
+
+      !! allocate variables
+      CALL wrk_alloc(jpi,jpj,jpk, zhdiv, zhdiv2)
+
+      !-----initialise the divergence and tracer correction------
+      zhdiv(:,:,:) = 0._wp
+      rhdivdiff(:,:,:) = 0._wp
+      rhdivdiff_trac(:,:,:,:) = 0._wp      
+      zhdiv2(:,:,:) = 0._wp
+
+
+     !calculate the hdivn*e3tn with the last values before re-meshing (should only be different from the real hdivn for cavity only)
+     !hdivb and rot should remain the same as the previous timestep bcause it will be used in dyn_ldf, so we did not use div_cur to calculate rhdivdiff as we should only use div_cur once for div_cur(nit000) in order not to update hdivb and rot
+
+!----------------------------------------------hdivn before remeshing-------------------------------------
+   DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            DO jk = 1, jpk                                 ! Horizontal slab
+               zhdiv2(ji,jj,jk) =   &
+                  (  e2u(ji,jj)*pe3u_b(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  )*pe3u_b(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)       &
+                   + e1v(ji,jj)*pe3v_b(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1)*pe3v_b(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  )    &
+                  / ( e1t(ji,jj) * e2t(ji,jj) )*ptmask_b(ji,jj,jk)
+            END DO
+         END DO
+      ENDDO
+
+! velocity with the new mask
+      un(:,:,:)=un(:,:,:)*umask(:,:,:)
+      vn(:,:,:)=vn(:,:,:)*vmask(:,:,:)
+
+
+      !calling divcur to update hdivn (ingat: hdivb = hdivn(before re-meshing), yet that hdivn was not using the recent values, as it was not updated at the end of timestep  )
+      !CALL div_cur(nit000)
+!----------------------------------------------hdivn after remeshing-------------------------------------
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            DO jk = 1, jpk                                 ! Horizontal slab
+                 zhdiv(ji,jj,jk) =   &
+                  (  e2u(ji,jj)*fse3u_n(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  )*fse3u_n(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)       &
+                   + e1v(ji,jj)*fse3v_n(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1)*fse3v_n(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  )    &
+                  / ( e1t(ji,jj) * e2t(ji,jj) )*tmask(ji,jj,jk)
+            END DO
+         END DO
+      ENDDO
+
+      zhdiv2(:,:,:) = zhdiv2(:,:,:) - zhdiv(:,:,:)  ! non-cavity will be = 0
+      rhdivdiff(:,:,:) =  zhdiv2(:,:,:)  
+
+      !----transfers the divergence of closed cells to the new top cell------
+      !until this point, rhdivdiff has the same unit as hdivn*e3t
+      DO jj = 2, jpjm1
+        DO ji = 2, jpim1
+          jk = mikt(ji,jj) 
+         IF (jk > 1 .AND. ssmask(ji,jj) == 1 ) THEN 
+            rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + SUM(rhdivdiff(ji,jj,1:jk-1)*ptmask_b(ji,jj,1:jk-1))  
+         ENDIF
+
+         jk = mbkt(ji,jj)
+         IF (ssmask(ji,jj) == 1 ) THEN       
+            rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + SUM(rhdivdiff(ji,jj,jk+1:jpk)*ptmask_b(ji,jj,jk+1:jpk)) 
+        ENDIF
+
+        END DO
+      ENDDO
+ 
+
+      !----cancel all the divergence of closed cells with tmask, then convert it to hdivn by dividing with e3t_n (it is now safe to do this) 
+      ! also cancelling totally closed cavity column 
+      rhdivdiff(:,:,:) = rhdivdiff(:,:,:)*tmask(:,:,:)/fse3t_n(:,:,:)
+      !tracer correction trend associated with the divergence correction; opposite sign of divergence
+      rhdivdiff_trac(:,:,:,1) = - rhdivdiff(:,:,:) * tsb(:,:,:,1)*tmask(:,:,:)
+      rhdivdiff_trac(:,:,:,2) = - rhdivdiff(:,:,:) * tsb(:,:,:,2)*tmask(:,:,:)
+
+      CALL lbc_lnk(rhdivdiff,   'T',1._wp)
+      CALL lbc_lnk(rhdivdiff_trac(:,:,:,1),   'T',1._wp)
+      CALL lbc_lnk(rhdivdiff_trac(:,:,:,2),   'T',1._wp)
+
+     CALL wrk_dealloc(jpi,jpj,jpk, zhdiv, zhdiv2 )
+
+   END SUBROUTINE iscpl_div_corr
+
+END MODULE iscpldiv
Index: anches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90	(revision 9812)
+++ 	(revision )
@@ -1,327 +1,0 @@
-MODULE iscplhsb
-   !!======================================================================
-   !!                       ***  MODULE  iscplhsb***
-   !! Ocean forcing: ice sheet/ocean coupling (conservation)
-   !!=====================================================================
-   !! History :  NEMO  ! 2015-01 P. Mathiot: original 
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   iscpl_alloc    : variable allocation
-   !!   iscpl_hsb      : compute and store the input of heat/salt/volume 
-   !!                    into the system due to the coupling process
-   !!   iscpl_div      : correction of divergence to keep volume conservation
-   !!----------------------------------------------------------------------
-   USE dom_oce         ! ocean space and time domain
-   USE domwri          ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE sbc_oce         ! surface boundary condition variables
-   USE oce             ! global tra/dyn variable
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! MPP library
-   USE lib_fortran     ! MPP library
-   USE wrk_nemo        ! Memory allocation
-   USE lbclnk          !
-   USE domngb          !
-   USE iscplini
-
-   IMPLICIT NONE
-   PRIVATE
-   
-   PUBLIC   iscpl_div   
-   PUBLIC   iscpl_cons       
-   !! * Substitutions  
-#  include "domzgr_substitute.h90"  
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE iscpl_cons(ptmask_b, psmask_b, pe3t_b, pts_flx, pvol_flx, prdt_iscpl)
-      !!---------------------------------------------------------------------- 
-      !!                   ***  ROUTINE iscpl_cons  ***
-      !! 
-      !! ** Purpose :   compute input into the system during the coupling step
-      !!                compute the correction term
-      !!                compute where the correction have to be applied
-      !! 
-      !! ** Method  :   compute tsn*e3t-tsb*e3tb and e3t-e3t_b
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b    !! mask before
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b      !! scale factor before
-      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b    !! mask before
-      REAL(wp), DIMENSION(:,:,:,:), INTENT(out) :: pts_flx     !! corrective flux to have tracer conservation
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(out) :: pvol_flx    !! corrective flux to have volume conservation
-      REAL(wp),                     INTENT(in ) :: prdt_iscpl  !! coupling period 
-      !!
-      INTEGER :: ji, jj, jk                                    !! loop index
-      INTEGER :: jip1, jim1, jjp1, jjm1
-      !!
-      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
-      REAL(wp):: r1_rdtiscpl
-      REAL(wp):: zjip1_ratio  , zjim1_ratio  , zjjp1_ratio  , zjjm1_ratio
-      !!
-      REAL(wp):: zde3t, zdtem, zdsal
-      REAL(wp), DIMENSION(:,:), POINTER :: zdssh
-      !!
-      REAL(wp), DIMENSION(:), ALLOCATABLE :: zlon, zlat
-      REAL(wp), DIMENSION(:), ALLOCATABLE :: zcorr_vol, zcorr_tem, zcorr_sal
-      INTEGER , DIMENSION(:), ALLOCATABLE :: ixpts, iypts, izpts, inpts
-      INTEGER :: jpts, npts
-
-      CALL wrk_alloc(jpi,jpj, zdssh )
-
-      ! get imbalance (volume heat and salt)
-      ! initialisation difference
-      zde3t = 0.0_wp; zdsal = 0.0_wp ; zdtem = 0.0_wp
-
-      ! initialisation correction term
-      pvol_flx(:,:,:  ) = 0.0_wp
-      pts_flx (:,:,:,:) = 0.0_wp
-      
-      r1_rdtiscpl = 1._wp / prdt_iscpl 
-
-      ! 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(:,:,:);
-
-      !==============================================================================
-      ! diagnose the heat, salt and volume input and compute the correction variable
-      !==============================================================================
-
-      ! 
-      zdssh(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
-      IF ( lk_vvl ) zdssh = 0.0_wp ! already included in the levels by definition
-      
-!      DO jk = 1,jpk-1
-!         DO jj = 2,jpj-1
-!            DO ji = fs_2,fs_jpim1
-!               IF (tmask_h(ji,jj) == 1._wp) THEN
-!
-!                  ! volume differences
-!                  zde3t = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
-!
-!                  ! heat diff
-!                  zdtem = 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 = 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._wp .OR. tmask(ji,jj,jk) == 1._wp ) THEN 
-!                     zde3t = zde3t + zdssh(ji,jj) ! zdssh = 0 if vvl
-!                     zdssh(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 * r1_rdtiscpl
-!                  pts_flx (ji,jj,jk,jp_sal)=   pts_flx (ji,jj,jk,jp_sal) + zdsal * r1_rdtiscpl 
-!                  pts_flx (ji,jj,jk,jp_tem)=   pts_flx (ji,jj,jk,jp_tem) + zdtem * r1_rdtiscpl
-!
-!                  ! case where we close a cell: check if the neighbour cells are wet 
-!                  IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
-!
-!                     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 /= 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 == 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) == 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
-!                           IF ( lwp ) THEN   ! JMM : cAution this warning may occur on any mpp subdomain but numout is only
-!                                             ! open for narea== 1 (lwp=T)
-!                           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 '
-!                           ENDIF
-!                        ! We deal with these 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.)
-
-      ! 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(inpts(jpnij))
-      inpts(:)=0
-
-      ! fill narea location with the number of problematic point
-!      DO jk = 1,jpk-1
-!         DO jj = 2,jpj-1
-!            DO ji = fs_2,fs_jpim1
-!               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
-!                  inpts(narea) = inpts(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(inpts,jpnij) 
-
-      ! size of the new variable
-      npts  = SUM(inpts)    
-      
-      ! 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_wp ; zlat(:) = -1.0e20_wp
-      zcorr_vol(:) = -1.0e20_wp
-      zcorr_sal(:) = -1.0e20_wp
-      zcorr_tem(:) = -1.0e20_wp
-
-      ! fill new variable
-!      jpts = SUM(inpts(1:narea-1))
-!      DO jk = 1,jpk-1
-!         DO jj = 2,jpj-1
-!            DO ji = fs_2,fs_jpim1
-!               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 inpts vector for the area narea
-!                  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
-!               END IF
-!            END DO
-!         END DO
-!      END DO
-
-      ! 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          
-!      DO jpts = 1,npts
-!         CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
-!         DO jj = mj0(iypts(jpts)),mj1(iypts(jpts))
-!            DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts))
-!               jk = izpts(jpts)
-!
-!               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(inpts)
-      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(:,:,:)
-
-      ! 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, zdssh ) 
-
-   END SUBROUTINE iscpl_cons
-
-   SUBROUTINE iscpl_div( phdivn )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE iscpl_div  ***
-      !!
-      !! ** Purpose :   update the horizontal divergenc
-      !!
-      !! ** Method  :
-      !!                CAUTION : iscpl is positive (inflow) and expressed in m/s
-      !!
-      !! ** Action  :   phdivn   increase by the iscpl correction term
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
-      !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) WRITE(numout,*) 'adding the divergence correction contribution at nit000' 
-
-      !DO jk = 1, jpk
-      !   DO jj = 1, jpj
-      !      DO ji = 1, jpi
-      !         !phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / fse3t_n(ji,jj,jk)
-      !         phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + rhdivdiff(ji,jj,jk)
-      !      END DO
-      !   END DO
-      !END DO
-      !
-      phdivn(:,:,:) =  phdivn(:,:,:) + rhdivdiff(:,:,:)
-
-   END SUBROUTINE iscpl_div
-
-END MODULE iscplhsb
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90	(revision 9812)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90	(revision 9813)
@@ -28,6 +28,8 @@
    REAL(wp), PUBLIC                                        ::   rdt_iscpl
    !!                                                      !!* namsbc_iscpl namelist *
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) ::   hdiv_iscpl
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   htsc_iscpl
+   !!------array used for divergence correction
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  rhdivdiff
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  rhdivdiff_trac
+
    !! * Substitutions  
 #  include "domzgr_substitute.h90"  
@@ -43,5 +45,5 @@
       !!                ***  ROUTINE sbc_iscpl_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( htsc_iscpl(jpi,jpj,jpk,jpts) , hdiv_iscpl(jpi,jpj,jpk) , STAT=iscpl_alloc )
+      ALLOCATE( rhdivdiff(jpi,jpj,jpk), rhdivdiff_trac(jpi,jpj,jpk,2), STAT=iscpl_alloc )
          !
       IF( lk_mpp          )   CALL mpp_sum ( iscpl_alloc )
@@ -51,4 +53,5 @@
    SUBROUTINE iscpl_init()
       INTEGER ::   ios           ! Local integer output status for namelist read
+      INTEGER :: ierr
       NAMELIST/namsbc_iscpl/nn_fiscpl,ln_hsb,nn_drown
       !!----------------------------------------------------------------------
@@ -83,4 +86,6 @@
       END IF
 
+      IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'iscpl_init : unable to allocate arrays' )
+
    END SUBROUTINE iscpl_init
 
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 9812)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 9813)
@@ -24,5 +24,6 @@
    USE lbclnk          ! communication
    USE iscplini        ! ice sheet coupling: initialisation
-   USE iscplhsb        ! ice sheet coupling: conservation
+   !USE iscplhsb        ! ice sheet coupling: conservation
+   USE iscpldiv
 
    IMPLICIT NONE
@@ -80,22 +81,6 @@
       CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
 
-      !! compute correction if conservation needed
-      IF ( ln_hsb ) THEN
-         IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
-         CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl)
-      END IF
-         
       !! print mesh/mask
       IF( nmsh /= 0 .AND. ln_iscpl )   CALL dom_wri      ! Create a domain file
-
-      IF ( ln_hsb ) THEN
-         cfile='correction'
-         cfile = TRIM( cfile )
-         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
-         CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) )
-         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_close ( inum0 )
-      END IF
 
       CALL wrk_dealloc(jpi,jpj,jpk, ztmask_b,zumask_b,zvmask_b )  
@@ -150,5 +135,5 @@
       !!
       REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb, zhdivdiff, zfp_ui, zfm_ui, zfp_vi, zfm_vi
-      REAL(wp):: zdz, zdzm1, zdzp1, zvol
+      REAL(wp):: zdz, zdzm1, zdzp1, zvol, zuflux_sum, zvflux_sum
       REAL(wp),DIMENSION(2)::ztrdtsb
       !!
@@ -156,6 +141,6 @@
       REAL(wp), DIMENSION(:,:    ), POINTER :: zbub   , zbvb   , zbun  , zbvn
       REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  , zssh1, zhu1, zhv1
-      REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1
-      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp, zuflux, zvflux
+      REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1, isfmask
+      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp, zuflux, zvflux, zhdiv, zhdiv2
       REAL(wp), DIMENSION(:,:,:  ), POINTER :: zwmaskn, zwmaskb, ztmp3d
       REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0
@@ -165,6 +150,6 @@
       CALL wrk_alloc(jpi,jpj,jpk,2, zts0                                   )
       CALL wrk_alloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp, ztmp3d        ) 
-      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb, zuflux, zvflux                       ) 
-      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1                       ) 
+      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb, zuflux, zvflux, zhdiv, zhdiv2               ) 
+      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1, isfmask                       ) 
       CALL wrk_alloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
       CALL wrk_alloc(jpi,jpj,       zbub   , zbvb    , zbun , zbvn         ) 
@@ -277,95 +262,6 @@
       END IF
 
-!=============================================================================
-! compute velocity
-
-      ub(:,:,:)=un(:,:,:)
-      vb(:,:,:)=vn(:,:,:)
-
-      un(:,:,:)=ub(:,:,:)*umask(:,:,:)
-      vn(:,:,:)=vb(:,:,:)*vmask(:,:,:)
-    
-
-      !-----initialise the divergence and tracer correction------
-      rhdivdiff(:,:,:) = 0._wp
-      rhdivdiff_trac(:,:,:,:) = 0._wp      
-
-      !-----the change of flux in each velocity cell 
-      zuflux(:,:,:) = (pe3u_b(:,:,:)*pumask_b(:,:,:) - fse3u_n(:,:,:)*umask(:,:,:))*ub(:,:,:)
-      zvflux(:,:,:) = (pe3v_b(:,:,:)*pvmask_b(:,:,:) - fse3v_n(:,:,:)*vmask(:,:,:))*vb(:,:,:)
-     
-
-      ! compute divergence and tracer correction (whether the horizontal velocity cell is closed or open)    
-      DO jj = 1,jpj
-         DO ji = 1,jpim1   ! jpim1 rather than jpi for the correct lbc_lnk synchronisation
-            !--------u-cell is closed-----------  
-            DO jk  = 1,jpk
-                 !------the contribution to the west side of u-cell
-                 IF(tmask(ji,jj,jk) == 1) THEN   ! the west cell is still active 
-                    zvol = e1e2t(ji,jj)*fse3t(ji,jj,jk)
-                    rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + zuflux(ji,jj,jk)*e2u(ji,jj)/zvol 
-                    rhdivdiff_trac(ji,jj,jk,1:2) =  rhdivdiff_trac(ji,jj,jk,1:2) - zuflux(ji,jj,jk)*tsb(ji,jj,jk,1:2)*e2u(ji,jj)/zvol
-                 ELSEIF(mikt(ji,jj) > 1) THEN  ! the west cell becomes dry, put the waterflux to the top of the west column
-                    jtop=mikt(ji,jj)
-                    zvol=e1e2t(ji,jj)*fse3t(ji,jj,jtop)
-                    rhdivdiff(ji,jj,jtop) = rhdivdiff(ji,jj,jtop) + zuflux(ji,jj,jk)*e2u(ji,jj)/zvol
-                    rhdivdiff_trac(ji,jj,jtop,1:2) = rhdivdiff_trac(ji,jj,jtop,1:2)-zuflux(ji,jj,jk)*tsb(ji,jj,jtop,1:2)*e2u(ji,jj)/zvol
-                 !ELSE ! for the closing T-column : nothing to be done 
-                 ENDIF
-
-                 !------the contribution to the east side of u-cell
-                 IF(tmask(ji+1,jj,jk) == 1) THEN  ! the east cell is still active
-                    zvol = e1e2t(ji+1,jj)*fse3t(ji+1,jj,jk)
-                    rhdivdiff(ji+1,jj,jk) = rhdivdiff(ji+1,jj,jk) - zuflux(ji,jj,jk)*e2u(ji,jj)/zvol
-                    rhdivdiff_trac(ji+1,jj,jk,1:2) = rhdivdiff_trac(ji+1,jj,jk,1:2) + zuflux(ji,jj,jk)*tsb(ji+1,jj,jk,1:2)*e2u(ji,jj)/zvol
-                 ELSEIF(mikt(ji+1,jj) > 1) THEN    ! the east cell becomes dry, put the waterflux to the top of the east column
-                    jtop = mikt(ji+1,jj)
-                    zvol = e1e2t(ji+1,jj)*fse3t(ji+1,jj,jtop)
-                    rhdivdiff(ji+1,jj,jtop) = rhdivdiff(ji+1,jj,jtop) - zuflux(ji,jj,jk)*e2u(ji,jj)/zvol
-                    rhdivdiff_trac(ji+1,jj,jtop,1:2) = rhdivdiff_trac(ji+1,jj,jtop,1:2) + zuflux(ji,jj,jk)*tsb(ji+1,jj,jtop,1:2)*e2u(ji,jj)/zvol
-                 !ELSE  ! for the closing T-column : nothing to be done
-                 ENDIF
-            ENDDO
-         END DO
-      END DO
-
-
-      DO jj = 1,jpjm1  ! jpim1 rather than jpi for the correct lbc_lnk synchronisation
-         DO ji = 1,jpi   
-            !-----------v-cell is closed-------------
-            DO jk  = 1,jpk
-                 !------the contribution to the south side of v-cell
-                 IF(tmask(ji,jj,jk) == 1) THEN  ! the south cell is still active
-                   zvol = e1e2t(ji,jj)*fse3t(ji,jj,jk)
-                   rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + zvflux(ji,jj,jk)*e1v(ji,jj)/zvol
-                   rhdivdiff_trac(ji,jj,jk,1:2) =  rhdivdiff_trac(ji,jj,jk,1:2) - zvflux(ji,jj,jk)*tsb(ji,jj,jk,1:2)*e1v(ji,jj)/zvol
-                 ELSEIF(mikt(ji,jj) > 1) THEN  !the south cell becomes dry, put the waterflux to the top of the south column
-                   jtop = mikt(ji,jj)
-                   zvol = e1e2t(ji,jj)*fse3t(ji,jj,jtop)
-                   rhdivdiff(ji,jj,jtop) = rhdivdiff(ji,jj,jtop) + zvflux(ji,jj,jk)*e1v(ji,jj)/zvol
-                   rhdivdiff_trac(ji,jj,jtop,1:2) = rhdivdiff_trac(ji,jj,jtop,1:2) - zvflux(ji,jj,jk)*tsb(ji,jj,jtop,1:2)*e1v(ji,jj)/zvol
-                 ENDIF
-
-                 !------the contribution to the north side of v-cell
-                 IF(tmask(ji,jj+1,jk) == 1) THEN ! the north cell is still active
-                    zvol = e1e2t(ji,jj+1)*fse3t(ji,jj+1,jk)
-                    rhdivdiff(ji,jj+1,jk) = rhdivdiff(ji,jj+1,jk) - zvflux(ji,jj,jk)*e1v(ji,jj)/zvol
-                    rhdivdiff_trac(ji,jj+1,jk,1:2) = rhdivdiff_trac(ji,jj+1,jk,1:2) + zvflux(ji,jj,jk)*tsb(ji,jj+1,jk,1:2)*e1v(ji,jj)/zvol
-                 ELSEIF(mikt(ji,jj+1) > 1) THEN  ! the north cell becomes dry, put the waterflux to the top of the north column
-                    jtop = mikt(ji,jj+1)
-                    zvol = e1e2t(ji,jj+1)*fse3t(ji,jj+1,jtop)
-                    rhdivdiff(ji,jj+1,jtop) = rhdivdiff(ji,jj+1,jtop) - zvflux(ji,jj,jk)*e1v(ji,jj)/zvol
-                    rhdivdiff_trac(ji,jj+1,jtop,1:2) = rhdivdiff_trac(ji,jj+1,jtop,1:2) + zvflux(ji,jj,jk)*tsb(ji,jj+1,jtop,1:2)*e1v(ji,jj)/zvol
-                 ENDIF
-            ENDDO
-
-         END DO
-      END DO 
-      
-
-      CALL lbc_lnk(rhdivdiff,   'T',1._wp)
-      CALL lbc_lnk(rhdivdiff_trac(:,:,:,1),   'T',1._wp)
-      CALL lbc_lnk(rhdivdiff_trac(:,:,:,2),   'T',1._wp)
-
+      !----compute and save the divergence difference (before and after remeshing) 
+      CALL iscpl_div_corr(ptmask_b, pe3t_b, pe3u_b, pe3v_b)
 
       ! compute temp and salt
@@ -473,6 +369,6 @@
       CALL wrk_dealloc(jpi,jpj,jpk,2, zts0                                   )
       CALL wrk_dealloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp                ) 
-      CALL wrk_dealloc(jpi,jpj,jpk,   zwmaskn, zwmaskb , ztmp3d, zuflux, zvflux    ) 
-      CALL wrk_dealloc(jpi,jpj,       zsmask0, zsmask1                       ) 
+      CALL wrk_dealloc(jpi,jpj,jpk,   zwmaskn, zwmaskb , ztmp3d, zuflux, zvflux, zhdiv, zhdiv2    ) 
+      CALL wrk_dealloc(jpi,jpj,       zsmask0, zsmask1, isfmask                     ) 
       CALL wrk_dealloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
       CALL wrk_dealloc(jpi,jpj,       zbub   , zbvb    , zbun  , zbvn        ) 
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 9812)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 9813)
@@ -35,5 +35,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
-   USE iscplhsb
+   USE iscpldiv
 
    IMPLICIT NONE
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 9812)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 9813)
@@ -33,4 +33,5 @@
    USE timing          ! Timing
    USE eosbn2
+   USE iscpldiv
 #if defined key_asminc   
    USE asminc          ! Assimilation increment
@@ -223,7 +224,5 @@
       !---start the tendency with tracer source/sink from divergence correction
       IF ( ln_iscpl .AND. kt == nit000) THEN        
-         IF (lwp) WRITE(numout,*) 'add tendency tsa due to re-meshing at kt = nit000'
-         tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + rhdivdiff_trac(:,:,:,jp_tem)*tmask(:,:,:)
-         tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + rhdivdiff_trac(:,:,:,jp_sal)*tmask(:,:,:)
+         CALL iscpl_div_tra( tsa )
       ENDIF
 
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 9812)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 9813)
@@ -71,8 +71,4 @@
    !! Energy budget of the leads (open water embedded in sea ice)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fraqsr_1lev        !: fraction of solar net radiation absorbed in the first ocean level [-]
-
-  !!------array used for divergence correction
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  rhdivdiff
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  rhdivdiff_trac
  
    !! Arrays used in coupling when MEDUSA is present. These arrays need to be declared
@@ -133,5 +129,5 @@
          &     ge3rui(jpi,jpj)   , ge3rvi(jpi,jpj)   ,                     &
          &     grui(jpi,jpj)     , grvi(jpi,jpj)     ,                     &
-         &     riceload(jpi,jpj),  rhdivdiff(jpi,jpj,jpk), rhdivdiff_trac(jpi,jpj,jpk,2) ,  STAT=ierr(2) )
+         &     riceload(jpi,jpj),                             STAT=ierr(2) )
          !
       ALLOCATE( snwice_mass(jpi,jpj) , snwice_mass_b(jpi,jpj), snwice_fmass(jpi,jpj) , STAT=ierr(3) )
