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 12094)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DYN/dynspg_ts.F90	(revision 12095)
@@ -152,4 +152,14 @@
       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         !   -      -
@@ -181,4 +191,53 @@
 !     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 )
@@ -463,14 +522,18 @@
             DO jj = 1, jpj
                DO ji = 1, jpim1   ! not jpi-column
-                  zhup2_e(ji,jj) = hu_0(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)
+!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)
                END DO
             END DO
             DO jj = 1, jpjm1        ! not jpj-row
                DO ji = 1, jpi
-                  zhvp2_e(ji,jj) = hv_0(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)
+!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)
                END DO
             END DO
@@ -651,11 +714,15 @@
                   !                    ! 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
-                  zhu_bck = hu_0(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)
-                  zhv_bck = hv_0(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)
+                  !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
-                  z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
-                  z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
+                  !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)) )
                   !
                   ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      & 
@@ -682,7 +749,23 @@
        
          IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only)
-            hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
+!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
+            ENDDO
+        ENDDO
+         CALL lbc_lnk_multi( 'dynspg_ts', zdep_u, 'U', -1._wp )
+         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) )
-            hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(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) )
             CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  &
@@ -782,6 +865,13 @@
                   &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )      &
                   &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
-            END DO
-         END DO
+               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) )
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'dynspg_ts', zdep_u, 'U', -1._wp )
+         CALL lbc_lnk_multi( 'dynspg_ts', zdep_v, 'V', -1._wp )
+
          CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
          !
@@ -791,6 +881,8 @@
          END DO
          ! Save barotropic velocities not transport:
-         ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
-         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
+!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
 
