Index: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 2986)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 3008)
@@ -24,9 +24,9 @@
    PRIVATE
 
-   PUBLIC   dom_vvl       ! called by domain.F90
-   PUBLIC   dom_vvl_alloc ! called by nemogcm.F90
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ee_t, ee_u, ee_v, ee_f   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mut , muu , muv , muf    !: ??? 
+   PUBLIC   dom_vvl         ! called by domain.F90
+   PUBLIC   dom_vvl_2       ! called by domain.F90
+   PUBLIC   dom_vvl_alloc   ! called by nemogcm.F90
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mut , muu , muv , muf    !: 1/H_0 at t-,u-,v-,f-points 
 
    REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:)     ::   r2dt   ! vertical profile time-step, = 2 rdttra 
@@ -49,5 +49,4 @@
       !
       ALLOCATE( mut (jpi,jpj,jpk) , muu (jpi,jpj,jpk) , muv (jpi,jpj,jpk) , muf (jpi,jpj,jpk) ,     &
-         &      ee_t(jpi,jpj)     , ee_u(jpi,jpj)     , ee_v(jpi,jpj)     , ee_f(jpi,jpj)     ,     &
          &      r2dt        (jpk)                                                             , STAT=dom_vvl_alloc )
          !
@@ -62,16 +61,18 @@
       !!                ***  ROUTINE dom_vvl  ***
       !!                   
-      !! ** Purpose :  compute coefficients muX at T-U-V-F points to spread
-      !!               ssh over the whole water column (scale factors)
+      !! ** Purpose :   compute mu coefficients at t-, u-, v- and f-points to 
+      !!              spread ssh over the whole water column (scale factors)
+      !!                set the before and now ssh at u- and v-points 
+      !!              (also f-point in now case)
       !!----------------------------------------------------------------------
       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3     ! 2D workspace
+      USE wrk_nemo, ONLY:   zee_t => wrk_2d_1, zee_u => wrk_2d_2, zee_v => wrk_2d_3, zee_f => wrk_2d_4   ! 2D workspace
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zcoefu , zcoefv   , zcoeff                   ! local scalars
-      REAL(wp) ::   zv_t_ij, zv_t_ip1j, zv_t_ijp1, zv_t_ip1jp1   !   -      -
-      !!----------------------------------------------------------------------
-
-      IF( wrk_in_use(2, 1,2,3) ) THEN
+      REAL(wp) ::   zcoefu, zcoefv , zcoeff                ! local scalars
+      REAL(wp) ::   zvt   , zvt_ip1, zvt_jp1, zvt_ip1jp1   !   -      -
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(2, 1,2,3,4) ) THEN
          CALL ctl_stop('dom_vvl: requested workspace arrays unavailable')   ;   RETURN
       ENDIF
@@ -97,35 +98,35 @@
 
       !                                 !==  mu computation  ==!
-      ee_t(:,:) = fse3t_0(:,:,1)                ! Lower bound : thickness of the first model level
-      ee_u(:,:) = fse3u_0(:,:,1)
-      ee_v(:,:) = fse3v_0(:,:,1)
-      ee_f(:,:) = fse3f_0(:,:,1)
+      zee_t(:,:) = fse3t_0(:,:,1)                ! Lower bound : thickness of the first model level
+      zee_u(:,:) = fse3u_0(:,:,1)
+      zee_v(:,:) = fse3v_0(:,:,1)
+      zee_f(:,:) = fse3f_0(:,:,1)
       DO jk = 2, jpkm1                          ! Sum of the masked vertical scale factors
-         ee_t(:,:) = ee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk)
-         ee_u(:,:) = ee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
-         ee_v(:,:) = ee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
+         zee_t(:,:) = zee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk)
+         zee_u(:,:) = zee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
+         zee_v(:,:) = zee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
          DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask
-            ee_f(:,jj) = ee_f(:,jj) + fse3f_0(:,jj,jk) *  umask(:,jj,jk) * umask(:,jj+1,jk)
+            zee_f(:,jj) = zee_f(:,jj) + fse3f_0(:,jj,jk) *  umask(:,jj,jk) * umask(:,jj+1,jk)
          END DO
       END DO  
       !                                         ! Compute and mask the inverse of the local depth at T, U, V and F points
-      ee_t(:,:) = 1. / ee_t(:,:) * tmask(:,:,1)
-      ee_u(:,:) = 1. / ee_u(:,:) * umask(:,:,1)
-      ee_v(:,:) = 1. / ee_v(:,:) * vmask(:,:,1)
+      zee_t(:,:) = 1._wp / zee_t(:,:) * tmask(:,:,1)
+      zee_u(:,:) = 1._wp / zee_u(:,:) * umask(:,:,1)
+      zee_v(:,:) = 1._wp / zee_v(:,:) * vmask(:,:,1)
       DO jj = 1, jpjm1                               ! f-point case fmask cannot be used 
-         ee_f(:,jj) = 1. / ee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1)
-      END DO
-      CALL lbc_lnk( ee_f, 'F', 1. )                  ! lateral boundary condition on ee_f
+         zee_f(:,jj) = 1._wp / zee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1)
+      END DO
+      CALL lbc_lnk( zee_f, 'F', 1. )                 ! lateral boundary condition on ee_f
       !
       DO jk = 1, jpk                            ! mu coefficients
-         mut(:,:,jk) = ee_t(:,:) * tmask(:,:,jk)     ! T-point at T levels
-         muu(:,:,jk) = ee_u(:,:) * umask(:,:,jk)     ! U-point at T levels
-         muv(:,:,jk) = ee_v(:,:) * vmask(:,:,jk)     ! V-point at T levels
+         mut(:,:,jk) = zee_t(:,:) * tmask(:,:,jk)     ! T-point at T levels
+         muu(:,:,jk) = zee_u(:,:) * umask(:,:,jk)     ! U-point at T levels
+         muv(:,:,jk) = zee_v(:,:) * vmask(:,:,jk)     ! V-point at T levels
       END DO
       DO jk = 1, jpk                                 ! F-point : fmask=shlat at coasts, use the product of umask
          DO jj = 1, jpjm1
-               muf(:,jj,jk) = ee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk)   ! at T levels
-         END DO
-         muf(:,jpj,jk) = 0.e0
+               muf(:,jj,jk) = zee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk)   ! at T levels
+         END DO
+         muf(:,jpj,jk) = 0._wp
       END DO
       CALL lbc_lnk( muf, 'F', 1. )                   ! lateral boundary condition
@@ -139,29 +140,23 @@
       END DO
       
-      ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
-      ! for ssh and scale factors
-      zs_t  (:,:) =         e1t(:,:) * e2t(:,:)
-      zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )
-      zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )
-
       DO jj = 1, jpjm1                          ! initialise before and now Sea Surface Height at u-, v-, f-points
          DO ji = 1, jpim1   ! NO vector opt.
-            zcoefu = umask(ji,jj,1) * zs_u_1(ji,jj)
-            zcoefv = vmask(ji,jj,1) * zs_v_1(ji,jj)
-            zcoeff = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) / ( e1f(ji,jj) * e2f(ji,jj) )
-            ! before fields
-            zv_t_ij       = zs_t(ji  ,jj  ) * sshb(ji  ,jj  )
-            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshb(ji+1,jj  )
-            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshb(ji  ,jj+1)
-            sshu_b(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j )
-            sshv_b(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 )
-            ! now fields
-            zv_t_ij       = zs_t(ji  ,jj  ) * sshn(ji  ,jj  )
-            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshn(ji+1,jj  )
-            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1)
-            zv_t_ip1jp1   = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1)
-            sshu_n(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j )
-            sshv_n(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 )
-            sshf_n(ji,jj) = zcoeff * ( zv_t_ij + zv_t_ip1j + zv_t_ijp1 + zv_t_ip1jp1 )
+            zcoefu = 0.50_wp / ( e1u(ji,jj) * e2u(ji,jj) ) * umask(ji,jj,1)
+            zcoefv = 0.50_wp / ( e1v(ji,jj) * e2v(ji,jj) ) * vmask(ji,jj,1)
+            zcoeff = 0.25_wp / ( e1f(ji,jj) * e2f(ji,jj) ) * umask(ji,jj,1) * umask(ji,jj+1,1)
+            !
+            zvt           = e1e2t(ji  ,jj  ) * sshb(ji  ,jj  )    ! before fields
+            zvt_ip1       = e1e2t(ji+1,jj  ) * sshb(ji+1,jj  )
+            zvt_jp1       = e1e2t(ji  ,jj+1) * sshb(ji  ,jj+1)
+            sshu_b(ji,jj) = zcoefu * ( zvt + zvt_ip1 )
+            sshv_b(ji,jj) = zcoefv * ( zvt + zvt_jp1 )
+            !
+            zvt           = e1e2t(ji  ,jj  ) * sshn(ji  ,jj  )    ! now fields
+            zvt_ip1       = e1e2t(ji+1,jj  ) * sshn(ji+1,jj  )
+            zvt_jp1       = e1e2t(ji  ,jj+1) * sshn(ji  ,jj+1)
+            zvt_ip1jp1    = e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1)
+            sshu_n(ji,jj) = zcoefu * ( zvt + zvt_ip1 )
+            sshv_n(ji,jj) = zcoefv * ( zvt + zvt_jp1 )
+            sshf_n(ji,jj) = zcoeff * ( zvt + zvt_ip1 + zvt_jp1 + zvt_ip1jp1 )
          END DO
       END DO
@@ -169,28 +164,298 @@
       CALL lbc_lnk( sshv_n, 'V', 1. )   ;   CALL lbc_lnk( sshv_b, 'V', 1. )
       CALL lbc_lnk( sshf_n, 'F', 1. )
-
-                                                ! initialise before scale factors at (u/v)-points
-      ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
-      DO jk = 1, jpkm1
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
-               zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
-               zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
-               fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
-               fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions
-      CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
-      ! Add initial scale factor to scale factor anomaly
-      fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
-      fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
-      !
-      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dom_vvl: failed to release workspace arrays')
+      !
+      IF( wrk_not_released(2, 1,2,3,4) )   CALL ctl_stop('dom_vvl: failed to release workspace arrays')
       !
    END SUBROUTINE dom_vvl
 
+
+   SUBROUTINE dom_vvl_2( kt, pe3u_b, pe3v_b )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_2  ***
+      !!                   
+      !! ** Purpose :   compute the vertical scale factors at u- and v-points
+      !!              in variable volume case.
+      !!
+      !! ** Method  :   In variable volume case (non linear sea surface) the 
+      !!              the vertical scale factor at velocity points is computed
+      !!              as the average of the cell surface weighted e3t.
+      !!                It uses the sea surface heigth so it have to be initialized
+      !!              after ssh is read/set
+      !!----------------------------------------------------------------------
+      INTEGER                   , INTENT(in   ) ::   kt               ! ocean time-step index
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe3u_b, pe3v_b   ! before vertical scale factor at u- & v-pts
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   iku, ikv     ! local integers    
+      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
+      REAL(wp) ::   zvt          ! local scalars
+      !!----------------------------------------------------------------------
+
+      IF( lwp .AND. kt == nit000 ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_vvl_2 : Variable volume, fse3t_b initialization'
+         WRITE(numout,*) '~~~~~~~~~ '
+         pe3u_b(:,:,jpk) = fse3u_0(:,:,jpk)
+         pe3v_b(:,:,jpk) = fse3u_0(:,:,jpk)
+      ENDIF
+      
+      DO jk = 1, jpkm1           ! set the before scale factors at u- & v-points
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zvt = fse3t_b(ji,jj,jk) * e1e2t(ji,jj)
+               pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1e2t(ji+1,jj) ) / ( e1u(ji,jj) * e2u(ji,jj) )
+               pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e1e2t(ji,jj+1) ) / ( e1v(ji,jj) * e2v(ji,jj) )
+            END DO
+         END DO
+      END DO
+
+      ! Correct scale factors at locations that have been individually modified in domhgr
+      ! Such modifications break the relationship between e1e2t and e1u*e2u etc. Recompute
+      ! scale factors ignoring the modified metric.
+      !                                                ! =====================
+      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
+         !                                             ! =====================
+         IF( nn_cla == 0 ) THEN
+            !
+            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
+            ij0 = 102   ;   ij1 = 102   
+            DO jk = 1, jpkm1                 ! set the before scale factors at u-points
+               DO jj = mj0(ij0), mj1(ij1)
+                  DO ji = mi0(ii0), mi1(ii1)
+                     zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
+                     pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
+                  END DO
+               END DO
+            END DO
+            !
+            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
+            ij0 =  88   ;   ij1 =  88   
+            DO jk = 1, jpkm1                 ! set the before scale factors at u-points
+               DO jj = mj0(ij0), mj1(ij1)
+                  DO ji = mi0(ii0), mi1(ii1)
+                     zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
+                     pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
+                  END DO
+               END DO
+            END DO
+            DO jk = 1, jpkm1                 ! set the before scale factors at v-points
+               DO jj = mj0(ij0), mj1(ij1)
+                  DO ji = mi0(ii0), mi1(ii1)
+                     zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
+                     pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
+                  END DO
+               END DO
+            END DO
+         ENDIF
+
+         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
+         ij0 = 116   ;   ij1 = 116   
+         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
+                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
+               END DO
+            END DO
+         END DO
+         !
+      ENDIF
+         !                                             ! =====================
+      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
+         !                                             ! =====================
+
+         ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified)
+         ij0 = 200   ;   ij1 = 200   
+         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
+                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
+               END DO
+            END DO
+         END DO
+
+         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
+         ij0 = 208   ;   ij1 = 208   
+         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
+                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
+               END DO
+            END DO
+         END DO
+
+         ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified)
+         ij0 = 124   ;   ij1 = 125   
+         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
+                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
+               END DO
+            END DO
+         END DO
+
+         ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
+         ij0 = 124   ;   ij1 = 125   
+         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
+                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
+               END DO
+            END DO
+         END DO
+
+         ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified)
+         ij0 = 124   ;   ij1 = 125   
+         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
+                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
+               END DO
+            END DO
+         END DO
+
+         ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified)
+         ij0 = 124   ;   ij1 = 125   
+         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
+                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
+               END DO
+            END DO
+         END DO
+
+         ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
+         ij0 = 141   ;   ij1 = 142   
+         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
+                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
+               END DO
+            END DO
+         END DO
+
+         ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
+         ij0 = 141   ;   ij1 = 142   
+         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
+                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
+               END DO
+            END DO
+         END DO
+
+         !
+      ENDIF
+      !                                                ! ======================
+      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
+         !                                             ! ======================
+         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
+         ij0 = 327   ;   ij1 = 327   
+         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
+                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
+               END DO
+            END DO
+         END DO
+         !
+         ii0 = 627   ;   ii1 = 628        ! Bosphore Strait (e2u was modified)
+         ij0 = 343   ;   ij1 = 343   
+         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
+                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
+               END DO
+            END DO
+         END DO
+         !
+         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
+         ij0 = 232   ;   ij1 = 232   
+         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
+                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
+               END DO
+            END DO
+         END DO
+         !
+         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
+         ij0 = 232   ;   ij1 = 232   
+         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
+                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
+               END DO
+            END DO
+         END DO
+         !
+         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
+         ij0 = 270   ;   ij1 = 270   
+         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
+                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
+               END DO
+            END DO
+         END DO
+         !
+         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
+         ij0 = 232   ;   ij1 = 233   
+         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
+                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
+               END DO
+            END DO
+         END DO
+         !
+         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
+         ij0 = 276   ;   ij1 = 276   
+         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
+                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
+               END DO
+            END DO
+         END DO
+         !
+      ENDIF
+      ! End of individual corrections to scale factors
+
+      IF( ln_zps ) THEN          ! minimum of the e3t at partial cell level
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               iku = mbku(ji,jj)
+               ikv = mbkv(ji,jj)
+               pe3u_b(ji,jj,iku) = MIN( fse3t_b(ji,jj,iku), fse3t_b(ji+1,jj  ,iku) ) 
+               pe3v_b(ji,jj,ikv) = MIN( fse3t_b(ji,jj,ikv), fse3t_b(ji  ,jj+1,ikv) ) 
+            END DO
+         END DO
+      ENDIF
+
+      pe3u_b(:,:,:) = pe3u_b(:,:,:) - fse3u_0(:,:,:)      ! anomaly to avoid zero along closed boundary/extra halos
+      pe3v_b(:,:,:) = pe3v_b(:,:,:) - fse3v_0(:,:,:)
+      CALL lbc_lnk( pe3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions
+      CALL lbc_lnk( pe3v_b(:,:,:), 'V', 1. )
+      pe3u_b(:,:,:) = pe3u_b(:,:,:) + fse3u_0(:,:,:)      ! recover the full scale factor
+      pe3v_b(:,:,:) = pe3v_b(:,:,:) + fse3v_0(:,:,:)
+      !
+   END SUBROUTINE dom_vvl_2
+   
 #else
    !!----------------------------------------------------------------------
@@ -200,4 +465,9 @@
    SUBROUTINE dom_vvl
    END SUBROUTINE dom_vvl
+   SUBROUTINE dom_vvl_2(kdum, pudum, pvdum )
+      USE par_kind
+      INTEGER                   , INTENT(in   ) ::   kdum       
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pudum, pvdum
+   END SUBROUTINE dom_vvl_2
 #endif
 
Index: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 2986)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 3008)
@@ -83,4 +83,6 @@
          neuler = 1                              ! Set time-step indicator at nit000 (leap-frog)
          CALL rst_read                           ! Read the restart file
+         !                                       ! define e3u_b, e3v_b from e3t_b read in restart file
+         CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) )
          CALL tra_swap                           ! swap 3D arrays (t,s)  in a 4D array (ts)
          CALL day_init                           ! model calendar (using both namelist and restart infos)
@@ -92,10 +94,13 @@
          CALL day_init                           ! model calendar (using both namelist and restart infos)
          !                                       ! Initialization of ocean to zero
-         !   before fields     !       now fields     
-         sshb (:,:)   = 0.e0   ;   sshn (:,:)   = 0.e0
-         ub   (:,:,:) = 0.e0   ;   un   (:,:,:) = 0.e0
-         vb   (:,:,:) = 0.e0   ;   vn   (:,:,:) = 0.e0  
-         rotb (:,:,:) = 0.e0   ;   rotn (:,:,:) = 0.e0
-         hdivb(:,:,:) = 0.e0   ;   hdivn(:,:,:) = 0.e0
+         !   before fields      !       now fields     
+         sshb (:,:)   = 0._wp   ;   sshn (:,:)   = 0._wp
+         ub   (:,:,:) = 0._wp   ;   un   (:,:,:) = 0._wp
+         vb   (:,:,:) = 0._wp   ;   vn   (:,:,:) = 0._wp  
+         rotb (:,:,:) = 0._wp   ;   rotn (:,:,:) = 0._wp
+         hdivb(:,:,:) = 0._wp   ;   hdivn(:,:,:) = 0._wp
+         !
+         !                                       ! define e3u_b, e3v_b from e3t_b initialized in domzgr
+         CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) )
          !
          IF( cp_cfg == 'eel' ) THEN
Index: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 2986)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 3008)
@@ -92,22 +92,16 @@
       !!               un,vn   now horizontal velocity of next time-step
       !!----------------------------------------------------------------------
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
       USE oce     , ONLY:   ze3u_f => ta       , ze3v_f => sa       ! (ta,sa) used as 3D workspace
-      USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3
       !
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   iku, ikv     ! local integers
 #if ! defined key_dynspg_flt
       REAL(wp) ::   z2dt         ! temporary scalar
 #endif
-      REAL(wp) ::   zue3a, zue3n, zue3b, zuf    ! local scalars
-      REAL(wp) ::   zve3a, zve3n, zve3b, zvf    !   -      -
-      REAL(wp) ::   zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1
+      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars
+      REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      -
       !!----------------------------------------------------------------------
-
-      IF( wrk_in_use(2, 1,2,3) ) THEN
-         CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN
-      ENDIF
 
       IF( kt == nit000 ) THEN
@@ -238,44 +232,20 @@
          ELSE                             ! Variable volume !
             !                             ! ================!
-            ! Before scale factor at t-points
-            ! -------------------------------
-            DO jk = 1, jpkm1
+            !
+            DO jk = 1, jpkm1                 ! Before scale factor at t-points
                fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
                   &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
-                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
-            ENDDO
-            ! Add volume filter correction only at the first level of t-point scale factors
-            zec = atfp * rdt / rau0
+                  &                         - 2._wp * fse3t_n(:,:,jk)            )
+            END DO
+            zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors
             fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
-            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
-            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
-            zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )
-            zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )
             !
-            IF( ln_dynadv_vec ) THEN
-               ! Before scale factor at (u/v)-points
-               ! -----------------------------------
-               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
-               DO jk = 1, jpkm1
-                  DO jj = 1, jpjm1
-                     DO ji = 1, jpim1
-                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
-                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
-                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
-                        fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
-                        fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
-                     END DO
-                  END DO
-               END DO
-               ! lateral boundary conditions
-               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
-               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
-               ! Add initial scale factor to scale factor anomaly
-               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
-               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
-               ! Leap-Frog - Asselin filter and swap: applied on velocity
-               ! -----------------------------------
-               DO jk = 1, jpkm1
-                  DO jj = 1, jpj
+            IF( ln_dynadv_vec ) THEN         ! vector invariant form (no thickness weighted calulation)
+               !
+               !                                      ! before scale factors at u- & v-pts (computed from fse3t_b)
+               CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) )
+               !
+               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap: applied on velocity
+                  DO jj = 1, jpj                      !                                                 --------
                      DO ji = 1, jpi
                         zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
@@ -290,30 +260,11 @@
                END DO
                !
-            ELSE
-               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
-               !-----------------------------------------------
-               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
-               DO jk = 1, jpkm1
-                  DO jj = 1, jpjm1
-                     DO ji = 1, jpim1
-                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
-                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
-                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
-                        ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
-                        ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
-                     END DO
-                  END DO
-               END DO
-               ! lateral boundary conditions
-               CALL lbc_lnk( ze3u_f, 'U', 1. )
-               CALL lbc_lnk( ze3v_f, 'V', 1. )
-               ! Add initial scale factor to scale factor anomaly
-               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
-               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
-               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
-               ! -----------------------------------             ===========================
-               DO jk = 1, jpkm1
-                  DO jj = 1, jpj
-                     DO ji = 1, jpim1
+            ELSE                             ! flux form (thickness weighted calulation)
+               !
+               CALL dom_vvl_2( kt, ze3u_f, ze3v_f )   ! before scale factors at u- & v-pts (computed from fse3t_b)
+               !
+               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap: 
+                  DO jj = 1, jpj                      !                   applied on thickness weighted velocity
+                     DO ji = 1, jpim1                 !                              ---------------------------
                         zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
                         zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
@@ -323,17 +274,17 @@
                         zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
                         !
-                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
-                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
+                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
+                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
                         !
-                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
+                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
                         vb(ji,jj,jk) = zvf
-                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
+                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
                         vn(ji,jj,jk) = va(ji,jj,jk)
                      END DO
                   END DO
                END DO
-               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
-               fse3v_b(:,:,:) = ze3v_f(:,:,:)
-               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
+               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
+               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
+               CALL lbc_lnk( ub, 'U', -1. )                    ! lateral boundary conditions
                CALL lbc_lnk( vb, 'V', -1. )
             ENDIF
@@ -346,6 +297,4 @@
          &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
       ! 
-      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
-      !
    END SUBROUTINE dyn_nxt
 
Index: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 2986)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 3008)
@@ -119,5 +119,5 @@
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
       INTEGER  ::   icycle           ! local scalar
-      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! local scalars
+      REAL(wp) ::   zraur, zcoef, z2dt_e, z1_2dt_b     ! local scalars
       REAL(wp) ::   z1_8, zx1, zy1                   !   -      -
       REAL(wp) ::   z1_4, zx2, zy2                   !   -      -
@@ -161,5 +161,6 @@
 
       !                                   !* Local constant initialization
-      z2dt_b = 2.0 * rdt                                    ! baroclinic time step
+      z1_2dt_b = 1._wp / ( 2.0_wp * rdt )                               ! baroclinic time step
+      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt_b = 1.0_wp / rdt    ! baroclinic time step (starting with euler timestep)
       z1_8 = 0.5 * 0.25                                     ! coefficient for vorticity estimates
       z1_4 = 0.5 * 0.5
@@ -195,6 +196,8 @@
                !
 #if defined key_vvl
-               ub_b(ji,jj) = ub_b(ji,jj) + (fse3u_0(ji,jj,jk)*(1.+sshu_b(ji,jj)*muu(ji,jj,jk)))* ub(ji,jj,jk) 
-               vb_b(ji,jj) = vb_b(ji,jj) + (fse3v_0(ji,jj,jk)*(1.+sshv_b(ji,jj)*muv(ji,jj,jk)))* vb(ji,jj,jk)   
+!              ub_b(ji,jj) = ub_b(ji,jj) + (fse3u_0(ji,jj,jk)*(1.+sshu_b(ji,jj)*muu(ji,jj,jk)))* ub(ji,jj,jk) 
+!              vb_b(ji,jj) = vb_b(ji,jj) + (fse3v_0(ji,jj,jk)*(1.+sshv_b(ji,jj)*muv(ji,jj,jk)))* vb(ji,jj,jk)   
+               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk)   *umask(ji,jj,jk) 
+               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk)   *vmask(ji,jj,jk)
 #else
                ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
@@ -278,5 +281,5 @@
       !                                             ! Remove barotropic contribution of bottom friction 
       !                                             ! from the barotropic transport trend
-      zcoef = -1. / z2dt_b
+      zcoef = -1._wp * z1_2dt_b
 # if defined key_vectopt_loop
       DO jj = 1, 1
@@ -565,6 +568,6 @@
       !                                   !* update the general momentum trend
       DO jk=1,jpkm1
-         ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b
-         va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b
+         ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b
+         va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b
       END DO
       un_b  (:,:) =  zu_sum(:,:) 
