Index: branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 2111)
+++ branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 2148)
@@ -4,7 +4,8 @@
    !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
    !!==============================================================================
-   !! History :      !  90-10  (B. Blanke)  Original code
-   !!                !  97-05  (G. Madec)  vertical component of isopycnal
-   !!           8.5  !  02-08  (G. Madec)  F90: Free form and module
+   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
+   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal
+   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!----------------------------------------------------------------------
 
@@ -13,9 +14,4 @@
    !!                  sion using a implicit time-stepping.
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
-   !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -28,6 +24,5 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dyn_zdf_imp    ! called by step.F90
+   PUBLIC   dyn_zdf_imp   ! called by step.F90
 
    !! * Substitutions
@@ -35,11 +30,10 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
-
 
    SUBROUTINE dyn_zdf_imp( kt, p2dt )
@@ -54,25 +48,21 @@
       !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
       !!      backward time stepping
-      !!      Surface boundary conditions: wind stress input
+      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
       !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
       !!      Add this trend to the general trend ua :
       !!         ua = ua + dz( avmu dz(u) )
       !!
-      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive
-      !!               mixing trend.
+      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
       !!---------------------------------------------------------------------
-      !! * Modules used
-      USE oce, ONLY :  zwd   => ta,   &                ! use ta as workspace
-                       zws   => sa                     ! use sa as workspace
-
-      !! * Arguments
-      INTEGER , INTENT( in ) ::   kt                   ! ocean time-step index
-      REAL(wp), INTENT( in ) ::  p2dt                  ! vertical profile of tracer time-step
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk                          ! dummy loop indices
-      REAL(wp) ::   zrau0r, z2dtf, zcoef, zzws, zrhs   ! temporary scalars
-      REAL(wp) ::   zzwi                               ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! temporary workspace arrays
+      USE oce, ONLY :  zwd   => ta      ! use ta as workspace
+      USE oce, ONLY :  zws   => sa      ! use sa as workspace
+      !!
+      INTEGER , INTENT( in ) ::   kt    ! ocean time-step index
+      REAL(wp), INTENT( in ) ::  p2dt   ! vertical profile of tracer time-step
+      !!
+      INTEGER  ::   ji, jj, jk             ! dummy loop indices
+      REAL(wp) ::   zrau0r, zcoef         ! temporary scalars
+      REAL(wp) ::   zzwi, zzws, zrhs       ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! 3D workspace
       !!----------------------------------------------------------------------
 
@@ -95,19 +85,18 @@
       ! friction velocity in dyn_bfr using values from the previous timestep. There
       ! is no need to include these in the implicit calculation.
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1        ! Matrix
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zcoef = - p2dt / fse3u(ji,jj,jk)
-               zzwi          = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
-               zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk)
-               zzws          = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
-               zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)
+               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
+               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
+               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
+               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
                zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
             END DO
          END DO
       END DO
-
-      ! Surface boudary conditions
-      DO jj = 2, jpjm1 
+      DO jj = 2, jpjm1        ! Surface boudary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0.
@@ -130,7 +119,6 @@
       !   The solution (the after velocity) is in ua
       !-----------------------------------------------------------------------
-
-      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
-      DO jk = 2, jpkm1
+      !
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -139,14 +127,9 @@
          END DO
       END DO
-
-      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
-      DO jj = 2, jpjm1   
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
-!!!         ua(ji,jj,1) = ub(ji,jj,1)  &
-!!!                      + p2dt * ( ua(ji,jj,1) + utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) )
-            z2dtf = p2dt / ( fse3u(ji,jj,1)*rau0 )
-            ua(ji,jj,1) = ub(ji,jj,1)  &
-                         + p2dt *  ua(ji,jj,1) + z2dtf * utau(ji,jj)
+      !
+      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ua(ji,jj,1) = ub(ji,jj,1)   &
+               &        + p2dt * (  ua(ji,jj,1) + 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 )  )
          END DO
       END DO
@@ -159,7 +142,6 @@
          END DO
       END DO
-
-      ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk
-      DO jj = 2, jpjm1   
+      !
+      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
             ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
@@ -173,7 +155,6 @@
          END DO
       END DO
-
-      ! Normalization to obtain the general momentum trend ua
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1        !==  Normalization to obtain the general momentum trend ua  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -192,11 +173,12 @@
       ! friction velocity in dyn_bfr using values from the previous timestep. There
       ! is no need to include these in the implicit calculation.
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1        ! Matrix
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zcoef = -p2dt / fse3v(ji,jj,jk)
-               zzwi          = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
+               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
                zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
-               zzws       = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
+               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
                zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
                zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
@@ -204,7 +186,5 @@
          END DO
       END DO
-
-      ! Surface boudary conditions
-      DO jj = 2, jpjm1   
+      DO jj = 2, jpjm1        ! Surface boudary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0.e0
@@ -223,12 +203,10 @@
       !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
       !
-      !   m is decomposed in the product of an upper and lower triangular
-      !   matrix
+      !   m is decomposed in the product of an upper and lower triangular matrix
       !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
       !   The solution (after velocity) is in 2d array va
       !-----------------------------------------------------------------------
-
-      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
-      DO jk = 2, jpkm1
+      !
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -237,14 +215,9 @@
          END DO
       END DO
-
-      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
-!!!         va(ji,jj,1) = vb(ji,jj,1)  &
-!!!                      + p2dt * ( va(ji,jj,1) + vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) )
-            z2dtf = p2dt / ( fse3v(ji,jj,1)*rau0 )
-            va(ji,jj,1) = vb(ji,jj,1)  &
-                         + p2dt * va(ji,jj,1) + z2dtf * vtau(ji,jj)
+      !
+      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            va(ji,jj,1) = vb(ji,jj,1)   &
+               &        + p2dt * (  va(ji,jj,1) + 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 )  )
          END DO
       END DO
@@ -257,7 +230,6 @@
          END DO
       END DO
-
-      ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
-      DO jj = 2, jpjm1   
+      !
+      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
             va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
@@ -271,7 +243,6 @@
          END DO
       END DO
-
-      ! Normalization to obtain the general momentum trend va
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1     !==  Normalization to obtain the general momentum trend va  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -280,5 +251,5 @@
          END DO
       END DO
-
+      !
    END SUBROUTINE dyn_zdf_imp
 
