Index: NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DOM/dom_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DOM/dom_oce.F90	(revision 12095)
+++ NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DOM/dom_oce.F90	(revision 12326)
@@ -150,4 +150,9 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
 
+!CEOD Scale of water column down to shllowest of neighbourinbg points over total
+!water depth
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: scaled_e3t_0_ik   , scaled_e3t_0_jk
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: scaled_e3t_0_ip1k , scaled_e3t_0_jp1k
+!CEOD
    INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1)
    INTEGER, PUBLIC ::   nlb10              !: shallowest W level Bellow ~10m (nla10 + 1) 
@@ -288,4 +293,8 @@
          !
       ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) )
+!CEOD 
+      ALLOCATE(  scaled_e3t_0_ik  (jpi,jpj) , scaled_e3t_0_jk  (jpi,jpj) ,      & 
+         &       scaled_e3t_0_ip1k(jpi,jpj) , scaled_e3t_0_jp1k(jpi,jpj) ,  STAT=ierr(13) )
+         !
       !
       dom_oce_alloc = MAXVAL(ierr)
Index: NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DOM/domain.F90	(revision 12095)
+++ NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DOM/domain.F90	(revision 12326)
@@ -79,4 +79,10 @@
       INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level
       REAL(wp), DIMENSION(jpi,jpj) ::   z1_hu_0, z1_hv_0
+!CEOD
+      REAL(wp), DIMENSION(jpi,jpj) ::  k_bot_i_min, k_bot_j_min
+      REAL(wp), DIMENSION(jpi,jpj) :: sum_e3t_0_min_ik
+      REAL(wp), DIMENSION(jpi,jpj) :: sum_e3t_0_min_jk
+      REAL(wp), DIMENSION(jpi,jpj) :: sum_e3t_0_min_ip1k
+      REAL(wp), DIMENSION(jpi,jpj) :: sum_e3t_0_min_jp1k
       !!----------------------------------------------------------------------
       !
@@ -155,4 +161,54 @@
          hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
       END DO
+!CEOD Here we need to work out whihc is the the shallowest point in terms of k
+!levels at neighbouring points in i and je directions.
+! then we need to wkr out the proprtion of the water column down to ht_0 that
+! that level is for both i points.
+! e.g point i could say go down to level 10 but point i+1 only level 5
+! then work out depth of bottom of level 5 at point i over depth at point i to
+! level 10
+! This is needed later to work out where a u point depth should be when using
+! enveloping coordinates, as it wil be 1/2 depth at t point at level 5 for i,i+1
+! points Because level 5 at i can change in time, thus u bed can change in time!
+! What about under ice shelves? Need rto return to that later
+
+! 1) get the shallowest k level at nighbrougin i and j points store in ! k_bot_i_min
+
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1                ! SPG with the application of W/D gravity filters
+                 k_bot_i_min(ji,jj) = MIN( ik_bot(ji,jj), ik_bot(ji+1,jj  ))
+                 k_bot_j_min(ji,jj) = MIN( ik_bot(ji,jj), ik_bot(ji,  jj+1))
+            ENDDO
+        ENDDO
+! 2) Work out the depth down to the shallowest k levle both at point i and i+1
+! (likewise for j) 
+
+        sum_e3t_0_min_ik(:,:) = 0
+        sum_e3t_0_min_jk(:,:) = 0
+        sum_e3t_0_min_ip1k(:,:) = 0
+        sum_e3t_0_min_jp1k(:,:) = 0
+
+        !Get sum of e3t_0s down to local min
+        DO jj = 1, jpjm1
+           DO ji = 1, jpim1               
+             DO jk = 1, k_bot_i_min(ji,jj)
+ 
+                sum_e3t_0_min_ik  (ji,jj) =  sum_e3t_0_min_ik  (ji,jj) + e3t_0(ji  ,jj,jk)*tmask(ji,jj,jk)
+                sum_e3t_0_min_ip1k(ji,jj) =  sum_e3t_0_min_ip1k(ji,jj) + e3t_0(ji+1,jj,jk)*tmask(ji+1,jj,jk)
+             ENDDO
+           
+             DO jk = 1, k_bot_j_min(ji,jj)
+                sum_e3t_0_min_jk  (ji,jj) =  sum_e3t_0_min_jk  (ji,jj) + e3t_0(ji,jj  ,jk)*tmask(ji,jj,jk)
+                sum_e3t_0_min_jp1k(ji,jj) =  sum_e3t_0_min_jp1k(ji,jj) + e3t_0(ji,jj+1,jk)*tmask(ji,jj+1,jk)
+             ENDDO
+! 3)  Then work out what fraction of the at rest water column that is, we later
+! multiply the now water depth by this scale to work out the bottom of the kth
+! level at time now in dynspg_ts 
+             scaled_e3t_0_ik  (ji,jj) = sum_e3t_0_min_ik  (ji,jj)/( ht_0(ji  ,jj  ) + 1._wp - ssmask(ji,jj))
+             scaled_e3t_0_jk  (ji,jj) = sum_e3t_0_min_jk  (ji,jj)/( ht_0(ji  ,jj  ) + 1._wp - ssmask(ji,jj))
+             scaled_e3t_0_ip1k(ji,jj) = sum_e3t_0_min_ip1k(ji,jj)/( ht_0(ji+1,jj  ) + 1._wp - ssmask(ji+1,jj))
+             scaled_e3t_0_jp1k(ji,jj) = sum_e3t_0_min_jp1k(ji,jj)/( ht_0(ji  ,jj+1) + 1._wp - ssmask(ji,jj+1))
+          ENDDO
+       ENDDO
       !
       !           !==  time varying part of coordinate system  ==!
Index: NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DYN/dynspg_ts.F90	(revision 12095)
+++ NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DYN/dynspg_ts.F90	(revision 12326)
@@ -152,14 +152,5 @@
       REAL(wp) ::   r1_2dt_b, z1_hu, z1_hv          ! local scalars
       REAL(wp) ::   za0, za1, za2, za3              !   -      -
-!CEOD
-      REAL(wp), DIMENSION(jpi,jpj) :: e3u_0_tot, e3v_0_tot
-      REAL(wp), DIMENSION(jpi,jpj) :: sum_e3t_0_min_ik
-      REAL(wp), DIMENSION(jpi,jpj) :: sum_e3t_0_min_jk
-      REAL(wp), DIMENSION(jpi,jpj) :: scaled_e3u_0_tot, scaled_e3v_0_tot
-      REAL(wp), DIMENSION(jpi,jpj) :: k_bot
-      REAL(wp), DIMENSION(jpi,jpj) :: k_bot_i_min
-      REAL(wp), DIMENSION(jpi,jpj) :: k_bot_j_min
       REAL(wp), DIMENSION(jpi,jpj) :: zdep_u, zdep_v
-!CEOD
       REAL(wp) ::   zmdi, zztmp, zldg               !   -      -
       REAL(wp) ::   zhu_bck, zhv_bck, zhdiv         !   -      -
@@ -191,53 +182,4 @@
 !     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
       !                                         ! inverse of baroclinic time step 
-!CEOD
-         e3u_0_tot(:,:)=0.0
-         e3v_0_tot(:,:)=0.0
-!CEOD
-         DO jk=1,jpk
-            e3u_0_tot(:,:) = e3u_0_tot(:,:) + e3u_0(:,:,jk)
-            e3v_0_tot(:,:) = e3v_0_tot(:,:) + e3v_0(:,:,jk)
-         ENDDO 
-         k_bot(:,:) = 1.
-         DO jk=1,jpkm1
-          DO jj = 1, jpj
-            DO ji = 1, jpi             ! SPG with the application of W/D gravity filters
-              IF ( tmask(ji,jj,jk) .ne. 0 ) THEN ! If its a wet point we keep going down
-                 k_bot(ji,jj) = jk
-              ENDIF
-            ENDDO
-           ENDDO
-         ENDDO 
-!work out k_bot min
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1                ! SPG with the application of W/D gravity filters
-                 k_bot_i_min(ji,jj) = MIN( k_bot(ji,jj), k_bot(ji+1,jj  ))
-                 k_bot_j_min(ji,jj) = MIN( k_bot(ji,jj), k_bot(ji,  jj+1))
-            ENDDO
-        ENDDO
-        sum_e3t_0_min_ik(:,:) = 0
-        sum_e3t_0_min_jk(:,:) = 0
-
-        !Get sum of e3t_0s sown to local min
-        DO jj = 1, jpjm1
-           DO ji = 1, jpim1               
-             DO jk = 1, k_bot_i_min(ji,jj)
-                sum_e3t_0_min_ik(ji,jj) =  sum_e3t_0_min_ik(ji,jj) + e3t_0(ji,jj,jk)
-             ENDDO
-           
-             DO jk = 1, k_bot_j_min(ji,jj)
-                sum_e3t_0_min_jk(ji,jj) =  sum_e3t_0_min_jk(ji,jj) + e3t_0(ji,jj,jk)
-             ENDDO
-             !CEOD Now scale that against ht_0 will be one in the case it matches of course(Potential for ht_0 = zero)
-             scaled_e3u_0_tot(ji,jj) = sum_e3t_0_min_ik(ji,jj)/( ht_0(ji,jj) + 1._wp - ssmask(ji,jj))
-             scaled_e3v_0_tot(ji,jj) = sum_e3t_0_min_jk(ji,jj)/( ht_0(ji,jj) + 1._wp - ssmask(ji,jj))
-!CEOD             WRITE(numout,*)'scaled values', ji,jj ,scaled_e3u_0_tot(ji,jj) , sum_e3t_0_min_ik(ji,jj), ht_0(ji,jj) + 1._wp - ssmask(ji,jj)
-          ENDDO
-       ENDDO
-!CEO to join things up
-         CALL lbc_lnk_multi( 'dynspg_ts', scaled_e3u_0_tot, 'U', -1._wp )
-         CALL lbc_lnk_multi( 'dynspg_ts', scaled_e3v_0_tot, 'V', -1._wp )
-              
-
       IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   r1_2dt_b = 1._wp / (         rdt )
       ELSE                                        ;   r1_2dt_b = 1._wp / ( 2._wp * rdt )
@@ -522,18 +464,12 @@
             DO jj = 1, jpj
                DO ji = 1, jpim1   ! not jpi-column
-!CEOD                  zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj)                        &
-!CEOD                       &                              * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
-!CEOD                       &                                 + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj)
-                  zhup2_e(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj) * zhtp2_e(ji  ,jj)*scaled_e3u_0_tot(ji  ,jj) &
-                                                 &         +  e1e2t(ji+1,jj) * zhtp2_e(ji+1,jj)*scaled_e3u_0_tot(ji+1,jj) ) *ssumask(ji,jj)
+                  zhup2_e(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj) * zhtp2_e(ji  ,jj)*scaled_e3t_0_ik  (ji,jj) &
+                                                 &         +  e1e2t(ji+1,jj) * zhtp2_e(ji+1,jj)*scaled_e3t_0_ip1k(ji,jj) ) *ssumask(ji,jj)
                END DO
             END DO
             DO jj = 1, jpjm1        ! not jpj-row
                DO ji = 1, jpi
-!CEOD                  zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj)                        &
-!CEOD                       &                              * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
-!CEOD                       &                                 + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj)
-                  zhvp2_e(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * zhtp2_e(ji,  jj)*scaled_e3v_0_tot(ji,jj  ) &
-                                                     &     +  e1e2t(ji,jj+1) * zhtp2_e(ji,jj+1)*scaled_e3v_0_tot(ji,jj+1) ) *ssvmask(ji,jj)
+                  zhvp2_e(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * zhtp2_e(ji,  jj)*scaled_e3t_0_jk  (ji,jj) &
+                                                     &     +  e1e2t(ji,jj+1) * zhtp2_e(ji,jj+1)*scaled_e3t_0_jp1k(ji,jj) ) *ssvmask(ji,jj)
                END DO
             END DO
@@ -714,16 +650,19 @@
                   !                    ! hu_e, hv_e hold depth at jn,  zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
                   !                    ! backward interpolated depth used in spg terms at jn+1/2
-                  !CEODzhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    &
-                  zhu_bck = (hu_0(ji,jj)/e3u_0_tot(ji,jj))*(e3u_0_tot(ji,jj)+ r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    &
-                       &                                          + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj) )
-                  !CEODzhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    &
-                  zhv_bck = (hv_0(ji,jj)/e3v_0_tot(ji,jj))*(e3v_0_tot(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    &
-                       &                                          + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj) )
-                  !                    ! inverse depth at jn+1
-                  !CEODz1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
-                  !CEODz1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
-                  z1_hu = ssumask(ji,jj) / ( (hu_0(ji,jj)/e3u_0_tot(ji,jj))*(e3u_0_tot(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj)) ) 
-                  z1_hv = ssvmask(ji,jj) / ( (hv_0(ji,jj)/e3v_0_tot(ji,jj))*(e3v_0_tot(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj)) )
-                  !
+                  zhu_bck = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj)*(ht_0(ji  ,jj) + zsshp2_e(ji  ,jj))*scaled_e3t_0_ik  (ji,jj) &
+                           &                         + e1e2t(ji+1,jj)*(ht_0(ji+1,jj) + zsshp2_e(ji+1,jj))*scaled_e3t_0_ip1k(ji,jj) &
+                           &        ) * ssumask(ji,jj) 
+
+                  zhv_bck = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  )*(ht_0(ji,jj  ) + zsshp2_e(ji,jj  ))*scaled_e3t_0_jk  (ji,jj) &
+                           &                         + e1e2t(ji,jj+1)*(ht_0(ji,jj+1) + zsshp2_e(ji,jj+1))*scaled_e3t_0_jp1k(ji,jj) &
+                           &        ) * ssvmask(ji,jj) 
+
+                  z1_hu = ssumask(ji,jj) / ( r1_2 * r1_e1e2u(ji,jj) *(e1e2t(ji  ,jj)*(ht_0(ji,jj)   + ssha_e(ji,jj)  )*scaled_e3t_0_ik(ji  ,jj) &
+                                            &      +e1e2t(ji+1,jj)*(ht_0(ji+1,jj) + ssha_e(ji+1,jj))*scaled_e3t_0_ip1k(ji,jj)) + 1._wp - ssumask(ji,jj) )
+
+                  z1_hv = ssvmask(ji,jj) / ( r1_2 * r1_e1e2v(ji,jj)* ( e1e2t(ji,jj  )*(ht_0(ji,jj)   + ssha_e(ji,jj)  )*scaled_e3t_0_jk(ji  ,jj) & 
+                                            &      + e1e2t(ji,jj+1)* (ht_0(ji,jj+1) + ssha_e(ji,jj+1))*scaled_e3t_0_jp1k(ji,jj)) + 1._wp - ssvmask(ji,jj) )
+
+
                   ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      & 
                        &            + rdtbt * (  zhu_bck        * zu_spg (ji,jj)  &   !
@@ -749,13 +688,10 @@
        
          IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only)
-!CEOD TBD
          DO jj = 1, jpjm1
             DO ji = 1, jpim1      ! NO Vector Opt.
-               zdep_u(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji,jj  ) * (ssha_e(ji,  jj)+ht_0(ji,  jj))*scaled_e3u_0_tot(ji  ,jj) &
-                                                 &     +  e1e2t(ji+1,jj) * (ssha_e(ji+1,jj)+ht_0(ji+1,jj))*scaled_e3u_0_tot(ji+1,jj) )*ssumask(ji,jj)
-               zdep_v(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * (ssha_e(ji,  jj)+ht_0(ji,  jj))*scaled_e3v_0_tot(ji,jj  ) &
-                                                 &     +  e1e2t(ji,jj+1) * (ssha_e(ji,jj+1)+ht_0(ji,jj+1))*scaled_e3v_0_tot(ji,jj+1) )*ssvmask(ji,jj)
-
-!CEOD             WRITE(numout,*)'zdep_u', ji,jj ,zdep_u(ji,jj) ,scaled_e3u_0_tot(ji,jj),ssha_e(ji,jj),jn
+               zdep_u(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji,jj  ) * (ssha_e(ji,  jj)+ht_0(ji,  jj))*scaled_e3t_0_ik  (ji,jj) &
+                                                 &     +  e1e2t(ji+1,jj) * (ssha_e(ji+1,jj)+ht_0(ji+1,jj))*scaled_e3t_0_ip1k(ji,jj) )*ssumask(ji,jj)
+               zdep_v(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * (ssha_e(ji,  jj)+ht_0(ji,  jj))*scaled_e3t_0_jk  (ji,jj  ) &
+                                                 &     +  e1e2t(ji,jj+1) * (ssha_e(ji,jj+1)+ht_0(ji,jj+1))*scaled_e3t_0_jp1k(ji,jj) )*ssvmask(ji,jj)
             ENDDO
         ENDDO
@@ -763,8 +699,6 @@
          CALL lbc_lnk_multi( 'dynspg_ts', zdep_v, 'V', -1._wp )
 
-            !CEODhu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
             hu_e (2:jpim1,2:jpjm1) = zdep_u(2:jpim1,2:jpjm1) 
             hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
-            !CEODhv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
             hv_e (2:jpim1,2:jpjm1) = zdep_v(2:jpim1,2:jpjm1)
             hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) )
@@ -865,8 +799,8 @@
                   &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )      &
                   &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
-               zdep_u(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji,jj  ) * (ssha(ji,  jj)+ht_0(ji,  jj))*scaled_e3u_0_tot(ji  ,jj) &
-                                                 &     +  e1e2t(ji+1,jj) * (ssha(ji+1,jj)+ht_0(ji+1,jj))*scaled_e3u_0_tot(ji+1,jj) )
-               zdep_v(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * (ssha(ji,  jj)+ht_0(ji,  jj))*scaled_e3v_0_tot(ji,jj  ) &
-                                                 &     +  e1e2t(ji,jj+1) * (ssha(ji,jj+1)+ht_0(ji,jj+1))*scaled_e3v_0_tot(ji,jj+1) )
+               zdep_u(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji,jj  ) * (ssha(ji,  jj)+ht_0(ji,  jj))*scaled_e3t_0_ik(ji  ,jj) &
+                                                 &     +  e1e2t(ji+1,jj) * (ssha(ji+1,jj)+ht_0(ji+1,jj))*scaled_e3t_0_ip1k(ji,jj) )
+               zdep_v(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * (ssha(ji,  jj)+ht_0(ji,  jj))*scaled_e3t_0_jk(ji,jj  ) &
+                                                 &     +  e1e2t(ji,jj+1) * (ssha(ji,jj+1)+ht_0(ji,jj+1))*scaled_e3t_0_jp1k(ji,jj) )
             END DO
          END DO
@@ -881,10 +815,7 @@
          END DO
          ! Save barotropic velocities not transport:
-!CEOD         ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
-!CEOD         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
          ua_b(:,:) =  ua_b(:,:) / ( zdep_u(:,:) + 1._wp - ssumask(:,:) )
          va_b(:,:) =  va_b(:,:) / ( zdep_v(:,:) + 1._wp - ssvmask(:,:) )
-      ENDIF
-
+      ENDIF 
 
       ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)  
