Index: trunk/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 474)
+++ trunk/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 503)
@@ -4,4 +4,14 @@
    !! Ocean dynamics:  hydrostatic pressure gradient trend
    !!======================================================================
+   !! History :  1.0  !  87-09  (P. Andrich, M.-A. Foujols)  hpg_zco: Original code
+   !!            5.0  !  91-11  (G. Madec)
+   !!            7.0  !  96-01  (G. Madec)  hpg_sco: Original code for s-coordinates
+   !!            8.0  !  97-05  (G. Madec)  split dynber into dynkeg and dynhpg
+   !!            8.5  !  02-07  (G. Madec)  F90: Free form and module
+   !!            8.5  !  02-08  (A. Bozec)  hpg_zps: Original code
+   !!            9.0  !  05-10  (A. Beckmann, B.W. An)  various s-coordinate options
+   !!                           Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot 
+   !!            9.0  !  05-11  (G. Madec) style & small optimisation
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -18,5 +28,4 @@
    !!       hpg_rot  : s-coordinate (ROTated axes scheme)
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -32,6 +41,5 @@
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC dyn_hpg                ! routine called by step.F90
+   PUBLIC   dyn_hpg    ! routine called by step module
 
 #if defined key_mpp_omp
@@ -49,20 +57,18 @@
 #endif
 
-   !! * Share module variables
-   LOGICAL  ::               & !!! ** nam_dynhpg **   hpg flags
-      ln_hpg_zco = .TRUE. ,  &  ! z-coordinate - full steps
-      ln_hpg_zps = .FALSE.,  &  ! z-coordinate - partial steps (interpolation)
-      ln_hpg_sco = .FALSE.,  &  ! s-coordinate (standard jacobian formulation)
-      ln_hpg_hel = .FALSE.,  &  ! s-coordinate (helsinki modification)
-      ln_hpg_wdj = .FALSE.,  &  ! s-coordinate (weighted density jacobian)
-      ln_hpg_djc = .FALSE.,  &  ! s-coordinate (Density Jacobian with Cubic polynomial)
-      ln_hpg_rot = .FALSE.      ! s-coordinate (ROTated axes scheme)
-
-   REAL(wp) ::               & !!! ** nam_dynhpg ** 
-      gamm = 0.e0               ! weighting coefficient
-
-   INTEGER  ::               &  ! 
-      nhpg  =  0                 ! = 0 to 6, type of pressure gradient scheme used
-      !                         ! (deduced from ln_hpg_... flags)
+   !!* Namelist nam_dynhpg : Choice of horizontal pressure gradient computation
+   LOGICAL  ::   ln_hpg_zco = .TRUE.    ! z-coordinate - full steps
+   LOGICAL  ::   ln_hpg_zps = .FALSE.   ! z-coordinate - partial steps (interpolation)
+   LOGICAL  ::   ln_hpg_sco = .FALSE.   ! s-coordinate (standard jacobian formulation)
+   LOGICAL  ::   ln_hpg_hel = .FALSE.   ! s-coordinate (helsinki modification)
+   LOGICAL  ::   ln_hpg_wdj = .FALSE.   ! s-coordinate (weighted density jacobian)
+   LOGICAL  ::   ln_hpg_djc = .FALSE.   ! s-coordinate (Density Jacobian with Cubic polynomial)
+   LOGICAL  ::   ln_hpg_rot = .FALSE.   ! s-coordinate (ROTated axes scheme)
+   REAL(wp) ::   gamm       = 0.e0      ! weighting coefficient
+   NAMELIST/nam_dynhpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_hel,   &
+      &                 ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, gamm
+
+   INTEGER  ::   nhpg  =  0             ! = 0 to 6, type of pressure gradient scheme used
+      !                                 ! (deduced from ln_hpg_... flags)
 
    !! * Substitutions
@@ -72,5 +78,5 @@
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! 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)
    !!----------------------------------------------------------------------
 
@@ -82,25 +88,17 @@
       !!
       !! ** Method  :   Call the hydrostatic pressure gradient routine 
-      !!      using the scheme defined in the namelist (nhpg parameter)
+      !!              using the scheme defined in the namelist
       !!   
       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
       !!             - Save the trend (l_trddyn=T)
-      !!             - Control print  (ln_ctl)
-      !!
-      !! History :
-      !!   9.0  !  05-10  (A. Beckmann, G. Madec) various s-coordinate options
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-
-      !! * local declarations
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         ztrdu, ztrdv                         ! 3D temporary workspace
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !!
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrdu, ztrdv   ! 3D temporary workspace
       !!----------------------------------------------------------------------
    
       IF( kt == nit000 )   CALL hpg_ctl      ! initialisation & control of options
 
-      ! Temporary saving of ua and va trends (l_trddyn)
-      IF( l_trddyn )   THEN
+      IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn)
          ztrdu(:,:,:) = ua(:,:,:)  
          ztrdv(:,:,:) = va(:,:,:) 
@@ -108,38 +106,25 @@
 
       SELECT CASE ( nhpg )      ! Hydrastatic pressure gradient computation
-      CASE ( 0 )                  ! z-coordinate 
-         CALL hpg_zco( kt )
-      CASE ( 1 )                  ! z-coordinate plus partial steps (interpolation)
-         CALL hpg_zps( kt )
-      CASE ( 2 )                  ! s-coordinate (standard jacobian formulation)
-         CALL hpg_sco( kt )
-      CASE ( 3 )                  ! s-coordinate (helsinki modification)
-         CALL hpg_hel( kt )
-      CASE ( 4 )                  ! s-coordinate (weighted density jacobian)
-         CALL hpg_wdj( kt )
-      CASE ( 5 )                  ! s-coordinate (Density Jacobian with Cubic polynomial)
-         CALL hpg_djc( kt )
-      CASE ( 6 )                  ! s-coordinate (ROTated axes scheme)
-         CALL hpg_rot( kt )
-      CASE ( 10 )                  ! z-coordinate
-         CALL hpg_zco_jki( kt )
-      CASE ( 11 )                  ! z-coordinate plus partial steps (interpolation)
-         CALL hpg_zps_jki( kt )
-      CASE ( 12 )                  ! s-coordinate (standard jacobian formulation)
-         CALL hpg_sco_jki( kt )
+      CASE (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate
+      CASE (  1 )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation)
+      CASE (  2 )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation)
+      CASE (  3 )   ;   CALL hpg_hel    ( kt )      ! s-coordinate (helsinki modification)
+      CASE (  4 )   ;   CALL hpg_wdj    ( kt )      ! s-coordinate (weighted density jacobian)
+      CASE (  5 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial)
+      CASE (  6 )   ;   CALL hpg_rot    ( kt )      ! s-coordinate (ROTated axes scheme)
+      CASE ( 10 )   ;   CALL hpg_zco_jki( kt )      ! z-coordinate (k-j-i)
+      CASE ( 11 )   ;   CALL hpg_zps_jki( kt )      ! z-coordinate plus partial steps (interpolation) (k-j-i)
+      CASE ( 12 )   ;   CALL hpg_sco_jki( kt )      ! s-coordinate (standard jacobian formulation) (k-j-i)
       END SELECT
 
-      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
-      IF( l_trddyn )   THEN
+      IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-         CALL trd_mod( ztrdu, ztrdv, jpdtdhpg, 'DYN', kt )
+         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )
       ENDIF          
-      
-      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
-         CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, &
-            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-      ENDIF
-      
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask,   &
+         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      !
    END SUBROUTINE dyn_hpg
 
@@ -154,21 +139,12 @@
       !! ** Action  :   Read the namelist namdynhpg and check the consistency
       !!      with the type of vertical coordinate used (zco, zps, sco)
-      !!
-      !! History :
-      !!   9.0  !  05-10  (A. Beckmann)  Original code
       !!----------------------------------------------------------------------
       INTEGER ::   ioptio = 0      ! temporary integer
-
-      NAMELIST/nam_dynhpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,               &
-         &                 ln_hpg_hel, ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot,   &
-         &                 gamm
-      !!----------------------------------------------------------------------
-
-      ! Read Namelist nam_dynhpg : pressure gradient calculation options
-      REWIND ( numnam )
+      !!----------------------------------------------------------------------
+
+      REWIND ( numnam )               ! Read Namelist nam_dynhpg : pressure gradient calculation options
       READ   ( numnam, nam_dynhpg )
 
-      ! Control print
-      IF(lwp) THEN
+      IF(lwp) THEN                    ! Control print
          WRITE(numout,*)
          WRITE(numout,*) 'dyn:hpg_ctl : hydrostatic pressure gradient control'
@@ -185,5 +161,5 @@
       ENDIF
 
-      ! set nhpg from ln_hpg_... flags
+      !                               ! Set nhpg from ln_hpg_... flags
       IF( ln_hpg_zco )   nhpg = 0
       IF( ln_hpg_zps )   nhpg = 1
@@ -194,5 +170,5 @@
       IF( ln_hpg_rot )   nhpg = 6
 
-      ! Consitency check
+      !                               ! Consitency check
       ioptio = 0 
       IF( ln_hpg_zco )   ioptio = ioptio + 1
@@ -203,6 +179,5 @@
       IF( ln_hpg_djc )   ioptio = ioptio + 1
       IF( ln_hpg_rot )   ioptio = ioptio + 1
-      IF ( ioptio > 1 )   &
-           &   CALL ctl_stop( ' several hydrostatic pressure gradient options used' )
+      IF ( ioptio /= 1 )   CALL ctl_stop( ' NO or several hydrostatic pressure gradient options used' )
 
       IF( lk_dynhpg_jki ) THEN
@@ -211,5 +186,5 @@
          IF(lwp) WRITE(numout,*) '          Autotasking or OPENMP: use j-k-i loops (i.e. _jki routines)'
       ENDIF
-
+      !
    END SUBROUTINE hpg_ctl
 
@@ -230,23 +205,12 @@
       !! 
       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
-      !!
-      !! History :
-      !!   1.0  !  87-09  (P. Andrich, M.-A. Foujols)  Original code
-      !!   5.0  !  91-11  (G. Madec)
-      !!   7.0  !  96-01  (G. Madec)
-      !!   8.0  !  97-05  (G. Madec)  split dynber into dynkeg and dynhpg
-      !!   8.5  !  02-07  (G. Madec)  F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * modules used
-      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace
-         &              zhpj => sa      ! use sa as 3D workspace
-      
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-      
-      !! * local declarations
-      INTEGER ::   ji, jj, jk           ! dummy loop indices
-      REAL(wp) ::   &   
-         zcoef0, zcoef1                 ! temporary scalars
+      !!----------------------------------------------------------------------
+      USE oce, ONLY :   zhpi => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   zhpj => sa   ! use sa as 3D workspace
+      !!
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk       ! dummy loop indices
+      REAL(wp) ::   zcoef0, zcoef1   ! temporary scalars
       !!----------------------------------------------------------------------
       
@@ -257,11 +221,8 @@
       ENDIF
       
-      
       ! Local constant initialization 
-      ! -----------------------------
       zcoef0 = - grav * 0.5
 
       ! Surface value
-      ! -------------
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -275,7 +236,6 @@
          END DO
       END DO
-
+      !
       ! interior value (2=<jk=<jpkm1)
-      ! --------------
       DO jk = 2, jpkm1
          DO jj = 2, jpjm1
@@ -296,5 +256,5 @@
          END DO
       END DO
-
+      !
    END SUBROUTINE hpg_zco
 
@@ -307,22 +267,13 @@
       !! 
       !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend
-      !!    
-      !! History :
-      !!   8.5  !  02-08  (A. Bozec)  Original code
-      !!   9.0  !  04-08  (G. Madec)  F90
       !!---------------------------------------------------------------------- 
-      !! * modules used                                 
-      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace
-         &              zhpj => sa      ! use sa as 3D workspace
-               
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-
-      !! * local declarations
-      INTEGER ::   ji, jj, jk           ! dummy loop indices
-      INTEGER ::   iku, ikv             ! temporary integers
-      REAL(wp) ::   &
-         zcoef0, zcoef1,             &  ! temporary scalars
-         zcoef2, zcoef3                 !    "         "
+      USE oce, ONLY :   zhpi => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   zhpj => sa   ! use sa as 3D workspace
+      !!
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk                       ! dummy loop indices
+      INTEGER  ::   iku, ikv                         ! temporary integers
+      REAL(wp) ::   zcoef0, zcoef1, zcoef2, zcoef3   ! temporary scalars
       !!----------------------------------------------------------------------
 
@@ -330,15 +281,11 @@
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate with partial steps'
-         IF(lwp) WRITE(numout,*) '              vector optimization'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate with partial steps - vector optimization'
       ENDIF
 
-
-      ! 0. Local constant initialization
-      ! --------------------------------
+      ! Local constant initialization
       zcoef0 = - grav * 0.5
 
-      ! 1. Surface value
-      ! ----------------
+      !  Surface value
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -353,6 +300,5 @@
       END DO
 
-      ! 2. interior value (2=<jk=<jpkm1)
-      ! -----------------
+      ! interior value (2=<jk=<jpkm1)
       DO jk = 2, jpkm1
          DO jj = 2, jpjm1
@@ -410,5 +356,5 @@
 # endif
       END DO
-
+      !
    END SUBROUTINE hpg_zps
 
@@ -431,23 +377,12 @@
       !!
       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
-      !!
-      !! History :
-      !!   7.0  !  96-01  (G. Madec)  s-coordinates
-      !!        !  97-05  (G. Madec)  split dynber into dynkeg and dynhpg
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module, vector opt.
-      !!   9.0  !  04-08  (C. Talandier) New trends organization
-      !!   9.0  !  05-10  (A. Beckmann) various s-coordinate options
-      !!----------------------------------------------------------------------
-      !! * modules used
-      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace
-         &              zhpj => sa      ! use sa as 3D workspace
-
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk           ! dummy loop indices
-      REAL(wp) ::   &
-         zcoef0, zuap, zvap             ! temporary scalars
+      !!----------------------------------------------------------------------
+      USE oce, ONLY :   zhpi => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   zhpj => sa   ! use sa as 3D workspace
+      !!
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      REAL(wp) ::   zcoef0, zuap, zvap   ! temporary scalars
       !!----------------------------------------------------------------------
 
@@ -458,12 +393,8 @@
       ENDIF
 
-
-      ! 0. Local constant initialization
-      ! --------------------------------
+      ! Local constant initialization
       zcoef0 = - grav * 0.5
 
-
-      ! 1. Surface value
-      ! ----------------                          
+      ! Surface value
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.   
@@ -484,7 +415,5 @@
       END DO   
             
-               
-      ! 2. interior value (2=<jk=<jpkm1)
-      ! -----------------     
+      ! interior value (2=<jk=<jpkm1)
       DO jk = 2, jpkm1                                  
          DO jj = 2, jpjm1     
@@ -508,5 +437,5 @@
          END DO
       END DO
-
+      !
    END SUBROUTINE hpg_sco
 
@@ -530,19 +459,12 @@
       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
       !!             - Save the trend (l_trddyn=T)
-      !!
-      !! History :
-      !!   9.0  !  05-10  (A. Beckmann)  Original code
-      !!----------------------------------------------------------------------
-      !! * modules used
-      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace
-         &              zhpj => sa      ! use sa as 3D workspace
-
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-      
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk           ! dummy loop indices
-      REAL(wp) ::   &
-         zcoef0, zuap, zvap             ! temporary scalars
+      !!----------------------------------------------------------------------
+      USE oce, ONLY :   zhpi => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   zhpj => sa   ! use sa as 3D workspace
+      !!
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      REAL(wp) ::   zcoef0, zuap, zvap   ! temporary scalars
       !!----------------------------------------------------------------------
 
@@ -553,12 +475,8 @@
       ENDIF
 
-
-      ! 0. Local constant initialization
-      ! --------------------------------
+      ! Local constant initialization
       zcoef0 = - grav * 0.5
-
  
-      ! 1. Surface value
-      ! ----------------
+      ! Surface value
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -578,7 +496,6 @@
          END DO
       END DO
-
-      ! 2. interior value (2=<jk=<jpkm1)
-      ! -----------------
+      !
+      ! interior value (2=<jk=<jpkm1)
       DO jk = 2, jpkm1
          DO jj = 2, jpjm1
@@ -606,5 +523,5 @@
          END DO
       END DO
-
+      !
    END SUBROUTINE hpg_hel
 
@@ -619,21 +536,13 @@
       !!
       !! Reference : Song, Mon. Wea. Rev., 126, 3213-3230, 1998.
-      !!
-      !! History :
-      !!   9.0  !  05-05  (B.W. An)  Original code
-      !!        !  05-10  (G. Madec) style & small optimisation
-      !!----------------------------------------------------------------------
-      !! * modules used
-      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace
-         &              zhpj => sa      ! use sa as 3D workspace
-
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk           ! dummy loop indices
-      REAL(wp) ::   &
-         zcoef0, zuap, zvap,         &  ! temporary scalars
-         zalph , zbeta                  !    "         "
+      !!----------------------------------------------------------------------
+      USE oce, ONLY :   zhpi => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   zhpj => sa   ! use sa as 3D workspace
+      !!
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      REAL(wp) ::   zcoef0, zuap, zvap   ! temporary scalars
+      REAL(wp) ::   zalph , zbeta        !    "         "
       !!----------------------------------------------------------------------
 
@@ -644,7 +553,5 @@
       ENDIF
 
-
       ! Local constant initialization
-      ! -----------------------------
       zcoef0 = - grav * 0.5
       zalph  = 0.5 - gamm        ! weighting coefficients (alpha=0.5-gamm)
@@ -652,5 +559,4 @@
 
       ! Surface value (no ponderation)
-      ! -------------
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -672,5 +578,4 @@
 
       ! Interior value (2=<jk=<jpkm1) (weighted with zalph & zbeta)
-      ! --------------
       DO jk = 2, jpkm1
          DO jj = 2, jpjm1
@@ -700,5 +605,5 @@
          END DO
       END DO
-
+      !
    END SUBROUTINE hpg_wdj
 
@@ -710,26 +615,18 @@
       !! ** Method  :   Density Jacobian with Cubic polynomial scheme
       !! 
-      !! Reference: Shchepetkin, A.F. & J.C. McWilliams, J. Geophys. Res.,
-      !!            108(C3), 3090, 2003
-      !! History :
-      !!   9.0  !  05-05  (B.W. An)  Original code
-      !!----------------------------------------------------------------------
-      !! * modules used
-      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace
-         &              zhpj => sa      ! use sa as 3D workspace
-
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-      
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk           ! dummy loop indices
-      REAL(wp) ::   &
-         zcoef0, z1_10, cffu, cffx,  &  ! temporary scalars
-                 z1_12, cffv, cffy,  &  !    "         "
-         zep   , cffw                   !    "         "
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &  ! 3D workspace
-         drhox, dzx, drhou, dzu, rho_i,     &
-         drhoy, dzy, drhov, dzv, rho_j,     &
-         drhoz, dzz, drhow, dzw, rho_k
+      !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003
+      !!----------------------------------------------------------------------
+      USE oce, ONLY :   zhpi => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   zhpj => sa   ! use sa as 3D workspace
+      !!
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk          ! dummy loop indices
+      REAL(wp) ::   zcoef0, zep, cffw   ! temporary scalars
+      REAL(wp) ::   z1_10, cffu, cffx   !    "         "
+      REAL(wp) ::   z1_12, cffv, cffy   !    "         "
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   drhox, dzx, drhou, dzu, rho_i   ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   drhoy, dzy, drhov, dzv, rho_j   !  "      "
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   drhoz, dzz, drhow, dzw, rho_k   !  "      "
       !!----------------------------------------------------------------------
 
@@ -741,6 +638,5 @@
 
 
-      ! 0. Local constant initialization
-      ! --------------------------------
+      ! Local constant initialization
       zcoef0 = - grav * 0.5
       z1_10  = 1.0 / 10.0
@@ -771,6 +667,6 @@
       zep = 1.e-15
 
-      !!bug  gm  drhoz not defined at level 1 and used (jk-1 with jk=2)
-      !!bug  gm  idem for drhox, drhoy et ji=jpi and jj=jpj
+!!bug  gm  drhoz not defined at level 1 and used (jk-1 with jk=2)
+!!bug  gm  idem for drhox, drhoy et ji=jpi and jj=jpj
 
       DO jk = 2, jpkm1
@@ -930,5 +826,5 @@
          END DO
       END DO
-
+      !
    END SUBROUTINE hpg_djc
 
@@ -941,25 +837,16 @@
       !!
       !! Reference: Thiem & Berntsen, Ocean Modelling, In press, 2005.
-      !! History :
-      !!   9.0  !  05-07 (B.W. An)
-      !!   9.0  !  05-10 (A. Beckmann) adapted to non-equidistant and masked grids
-      !!----------------------------------------------------------------------
-      !! * modules used
-      USE oce, ONLY :   zhpi => ta,         &  ! use ta as 3D workspace
-         &              zhpj => sa             ! use sa as 3D workspace
-
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt            ! ocean time-step index
-      
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk                  ! dummy loop indices
-      REAL(wp) ::   &
-         zforg, zcoef0, zuap, zmskd1, zmskd1m,                             &
-         zfrot        , zvap, zmskd2, zmskd2m
-      REAL(wp), DIMENSION(jpi,jpj) ::       &  ! 2D temporary workspace
-         zdistr, zsina, zcosa             
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &  ! 3D temporary workspace
-         zhpiorg, zhpirot, zhpitra, zhpine, & 
-         zhpjorg, zhpjrot, zhpjtra, zhpjne
+      !!----------------------------------------------------------------------
+      USE oce, ONLY :   zhpi => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   zhpj => sa   ! use sa as 3D workspace
+      !!
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk          ! dummy loop indices
+      REAL(wp) ::   zforg, zcoef0, zuap, zmskd1, zmskd1m   ! temporary scalar
+      REAL(wp) ::   zfrot        , zvap, zmskd2, zmskd2m   !    "         "
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zdistr, zsina, zcosa                ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpiorg, zhpirot, zhpitra, zhpine   ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpjorg, zhpjrot, zhpjtra, zhpjne   !  "      "
       !!----------------------------------------------------------------------
 
@@ -1118,5 +1005,5 @@
          END DO
       END DO
-
+      !
    END SUBROUTINE hpg_rot
 
Index: trunk/NEMO/OPA_SRC/DYN/dynhpg_jki.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynhpg_jki.F90	(revision 474)
+++ trunk/NEMO/OPA_SRC/DYN/dynhpg_jki.F90	(revision 503)
@@ -4,17 +4,19 @@
    !! Ocean dynamics:  hydrostatic pressure gradient trend
    !!======================================================================
-
-   !!----------------------------------------------------------------------
-   !!   hpg...o_jki : update the momentum trend with the horizontal
-   !!                  gradient of the hydrostatic pressure
-   !!----------------------------------------------------------------------
-   !! * Modules used
+   !! History :  9.0  !  06-09  (G. Madec) From dynhpg module
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   hpg_sco_jki : update the momentum trend with the horizontal
+   !!                 gradient of the hydrostatic pressure (s-coordinate)
+   !!   hpg_zps_jki : update the momentum trend with the horizontal
+   !!                 gradient of the hydrostatic pressure (partial step)
+   !!   hpg_zco_jki : update the momentum trend with the horizontal
+   !!                 gradient of the hydrostatic pressure (z-coordinate)
+   !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
    USE in_out_manager  ! I/O manager
-   USE trdmod          ! ocean dynamics trends
-   USE trdmod_oce      ! ocean variables trends
-   USE prtctl          ! Print control
    USE lbclnk          ! lateral boundary condition
 
@@ -22,8 +24,7 @@
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC hpg_sco_jki ! routine called by step.F90
-   PUBLIC hpg_zps_jki ! routine called by step.F90
-   PUBLIC hpg_zco_jki ! routine called by step.F90
+   PUBLIC   hpg_sco_jki   ! routine called by step.F90
+   PUBLIC   hpg_zps_jki   ! routine called by step.F90
+   PUBLIC   hpg_zco_jki   ! routine called by step.F90
 
    !! * Substitutions
@@ -32,4 +33,6 @@
    !!----------------------------------------------------------------------
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! $Header$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -57,25 +60,13 @@
       !!
       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
-      !!             - Save the trend in (utrd,vtrd) ('key_trddyn')
-      !!
-      !! History :
-      !!   7.0  !  96-01  (G. Madec)  s-coordinates
-      !!        !  97-05  (G. Madec)  split dynber into dynkeg and dynhpg
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
-      !!   9.0  !  04-08  (C. Talandier) New trends organization
-      !!----------------------------------------------------------------------
-      !! * modules used
-      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace
-         &              zhpj => sa      ! use sa as 3D workspace
-
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-      
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk           ! dummy loop indices
-      REAL(wp) ::   &
-         zcoef0, zuap, zvap     ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         ztdua, ztdva                   ! temporary scalars
+      !!             - Save the trend in (ztrdu,ztrdv) ('key_trddyn')
+      !!----------------------------------------------------------------------
+      USE oce, ONLY :   zhpi => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   zhpj => sa   ! use sa as 3D workspace
+      !!
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      REAL(wp) ::   zcoef0, zuap, zvap   ! temporary scalars
       !!----------------------------------------------------------------------
 
@@ -86,12 +77,5 @@
       ENDIF
 
-      ! Save ua and va trends
-      IF( l_trddyn )   THEN
-         ztdua(:,:,:) = ua(:,:,:) 
-         ztdva(:,:,:) = va(:,:,:) 
-      ENDIF
-
-      ! 0. Local constant initialization
-      ! --------------------------------
+      ! Local constant initialization
       zcoef0 = - grav * 0.5
       zuap   = 0.e0
@@ -101,7 +85,5 @@
       DO jj = 2, jpjm1                                 !  Vertical slab
          !                                             ! ===============
-         ! 1. Surface value
-         ! ----------------
-         DO ji = 2, jpim1
+         DO ji = 2, jpim1         !  Surface value
             ! hydrostatic pressure gradient along s-surfaces
             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj)   &
@@ -118,8 +100,6 @@
             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
          END DO  
-
-         ! 2. interior value (2=<jk=<jpkm1)
-         ! -----------------
-         DO jk = 2, jpkm1
+         !
+         DO jk = 2, jpkm1         ! interior value (2=<jk=<jpkm1)
             DO ji = 2, jpim1
                ! hydrostatic pressure gradient along s-surfaces
@@ -143,18 +123,4 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-
-      ! save the hydrostatic pressure gradient trends for diagnostic
-      ! momentum trends
-      IF( l_trddyn )   THEN
-         zhpi(:,:,:) = ua(:,:,:) - ztdua(:,:,:)
-         zhpj(:,:,:) = va(:,:,:) - ztdva(:,:,:)
-         CALL trd_mod(zhpi, zhpj, jpdtdhpg, 'DYN', kt)
-      ENDIF
-
-      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
-         CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, &
-            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
-      ENDIF
-
    END SUBROUTINE hpg_sco_jki
 
@@ -178,24 +144,15 @@
       !!
       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
-      !!             - Save the trend in (utrd,vtrd) ('key_trddyn')
-      !!
-      !! History :
-      !!   8.5  !  02-08  (A. Bozec)  Original code
-      !!----------------------------------------------------------------------
-      !! * modules used
-      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace
-         &              zhpj => sa      ! use sa as 3D workspace
-
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-
-      !! * local declarations
-      INTEGER ::   ji, jj, jk           ! dummy loop indices
-      INTEGER ::   iku, ikv             ! temporary integers
-      REAL(wp) ::   &
-         zcoef0, zcoef1, zuap,       &  ! temporary scalars
-         zcoef2, zcoef3, zvap           !    "         "
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         ztdua, ztdva                   ! temporary scalars
+      !!             - Save the trend in (ztrdu,ztrdv) ('key_trddyn')
+      !!----------------------------------------------------------------------
+      USE oce, ONLY :   zhpi => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   zhpj => sa   ! use sa as 3D workspace
+      !!
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk             ! dummy loop indices
+      INTEGER  ::   iku, ikv               ! temporary integers
+      REAL(wp) ::   zcoef0, zcoef1, zuap   ! temporary scalars
+      REAL(wp) ::   zcoef2, zcoef3, zvap   !    "         "
       !!----------------------------------------------------------------------
 
@@ -206,12 +163,5 @@
       ENDIF
 
-      ! Save ua and va trends
-      IF( l_trddyn )   THEN
-         ztdua(:,:,:) = ua(:,:,:) 
-         ztdva(:,:,:) = va(:,:,:) 
-      ENDIF
-
-      ! 0. Local constant initialization
-      ! --------------------------------
+      ! Local constant initialization
       zcoef0 = - grav * 0.5
       zuap   = 0.e0
@@ -220,7 +170,5 @@
       DO jj = 2, jpjm1                                 !  Vertical slab
          !                                             ! ===============
-         ! 1. Surface value
-         ! ----------------
-         DO ji = 2, jpim1
+         DO ji = 2, jpim1         ! Surface value
             zcoef1 = zcoef0 * fse3w(ji,jj,1)
             ! hydrostatic pressure gradient
@@ -231,8 +179,6 @@
             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
          END DO
-
-         ! 2. interior value (2=<jk=<jpkm1)
-         ! -----------------
-         DO jk = 2, jpkm1
+         !
+         DO jk = 2, jpkm1         ! interior value (2=<jk=<jpkm1)
             DO ji = 2, jpim1
                zcoef1 = zcoef0 * fse3w(ji,jj,jk)
@@ -250,5 +196,5 @@
             END DO 
          END DO
-
+         !
          ! partial steps correction at the last level  (new gradient with  intgrd.F)
          DO ji = 2, jpim1
@@ -281,18 +227,4 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-
-      ! save the hydrostatic pressure gradient trends for diagnostic
-      ! momentum trends
-      IF( l_trddyn )   THEN
-         zhpi(:,:,:) = ua(:,:,:) - ztdua(:,:,:)
-         zhpj(:,:,:) = va(:,:,:) - ztdva(:,:,:)
-         CALL trd_mod(zhpi, zhpj, jpdtdhpg, 'DYN', kt)
-      ENDIF
-
-      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
-         CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, &
-            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
-      ENDIF   
-
    END SUBROUTINE hpg_zps_jki
 
@@ -317,26 +249,13 @@
       !!
       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
-      !!             - Save the trend in (utrd,vtrd) ('key_trddyn')
-      !!
-      !! History :
-      !!   1.0  !  87-09  (P. Andrich, m.-a. Foujols)  Original code
-      !!        !  91-11  (G. Madec)
-      !!        !  96-01  (G. Madec)  s-coordinates
-      !!        !  97-05  (G. Madec)  split dynber into dynkeg and dynhpg
-      !!   8.5  !  02-07  (G. Madec)  F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * modules used
-      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace
-         &              zhpj => sa      ! use sa as 3D workspace
-
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-
-      !! * local declarations
-      INTEGER ::   ji, jj, jk           ! dummy loop indices
-      REAL(wp) ::   &
-         zcoef0, zcoef1, zuap, zvap     ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         ztdua, ztdva                   ! temporary scalars
+      !!             - Save the trend in (ztrdu,ztrdv) ('key_trddyn')
+      !!----------------------------------------------------------------------
+      USE oce, ONLY :   zhpi => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   zhpj => sa   ! use sa as 3D workspace
+      !!
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
+      REAL(wp) ::   zcoef0, zcoef1, zuap, zvap   ! temporary scalars
       !!----------------------------------------------------------------------
 
@@ -347,12 +266,5 @@
       ENDIF
 
-      ! Save ua and va trends
-      IF( l_trddyn )   THEN
-         ztdua(:,:,:) = ua(:,:,:) 
-         ztdva(:,:,:) = va(:,:,:) 
-      ENDIF
-
-      ! 0. Local constant initialization
-      ! --------------------------------
+      ! Local constant initialization
       zcoef0 = - grav * 0.5
       zuap   = 0.e0
@@ -362,9 +274,5 @@
       DO jj = 2, jpjm1                                 !  Vertical slab
          !                                             ! ===============
-         ! 1. Surface value
-         ! ----------------
-         
-
-         DO ji = 2, jpim1
+         DO ji = 2, jpim1         ! Surface value
             zcoef1 = zcoef0 * fse3w(ji,jj,1)
             ! hydrostatic pressure gradient
@@ -375,8 +283,6 @@
             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
          END DO
-
-         ! 2. interior value (2=<jk=<jpkm1)
-         ! -----------------
-         DO jk = 2, jpkm1
+         !
+         DO jk = 2, jpkm1         ! interior value (2=<jk=<jpkm1)
             DO ji = 2, jpim1
                zcoef1 = zcoef0 * fse3w(ji,jj,jk)
@@ -397,19 +303,4 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-
-      ! save the hydrostatic pressure gradient trends for diagnostic
-      ! momentum trends
-      IF( l_trddyn )   THEN
-         zhpi(:,:,:) = ua(:,:,:) - ztdua(:,:,:)
-         zhpj(:,:,:) = va(:,:,:) - ztdva(:,:,:)
-
-         CALL trd_mod(zhpi, zhpj, jpdtdhpg, 'DYN', kt)
-      ENDIF
-
-      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
-         CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, &
-            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
-      ENDIF
-
    END SUBROUTINE hpg_zco_jki
 
Index: trunk/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 474)
+++ trunk/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 503)
@@ -4,9 +4,12 @@
    !! Ocean dynamics:  kinetic energy gradient trend
    !!======================================================================
+   !! History :  1.0  !  87-09  (P. Andrich, m.-a. Foujols)  Original code
+   !!            7.0  !  97-05  (G. Madec)  Split dynber into dynkeg and dynhpg
+   !!            9.0  !  02-07  (G. Madec)  F90: Free form and module
+   !!----------------------------------------------------------------------
    
    !!----------------------------------------------------------------------
    !!   dyn_keg      : update the momentum trend with the horizontal tke
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -19,14 +22,13 @@
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC dyn_keg                ! routine called by step.F90
+   PUBLIC   dyn_keg    ! routine called by step module
    
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
-   !!---------------------------------------------------------------------------------
+   !!----------------------------------------------------------------------
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! 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
@@ -40,5 +42,5 @@
       !!      general momentum trend.
       !!
-      !! ** Method  :   Compute the now horizontal kinetic energy:
+      !! ** Method  :   Compute the now horizontal kinetic energy 
       !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
       !!      Take its horizontal gradient and add it to the general momentum
@@ -48,24 +50,14 @@
       !!
       !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
-      !!             - Save the trends in (utrd,vtrd) ('key_trddyn')
+      !!             - save this trends (l_trddyn=T) for post-processing
+      !!----------------------------------------------------------------------
+      USE oce, ONLY :   ztrdu => ta   ! use ta as 3D workspace   
+      USE oce, ONLY :   ztrdv => sa   ! use sa as 3D workspace   
       !!
-      !! History :
-      !!   1.0  !  87-09  (P. Andrich, m.-a. Foujols)  Original code
-      !!   7.0  !  97-05  (G. Madec)  Split dynber into dynkeg and dynhpg
-      !!   9.0  !  02-07  (G. Madec)  F90: Free form and module
-      !!    "   !  04-08  (C. Talandier) New trends organization
-      !!----------------------------------------------------------------------
-      !! * Modules used     
-      USE oce, ONLY :    ztdua => ta,   & ! use ta as 3D workspace   
-                         ztdva => sa      ! use sa as 3D workspace   
-
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk          ! dummy loop indices
-      REAL(wp) ::   zua, zva, zu, zv    ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         zhke                           ! temporary workspace
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zu, zv       ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhke   ! temporary 3D workspace
       !!----------------------------------------------------------------------
 
@@ -76,8 +68,7 @@
       ENDIF
 
-      ! Save ua and va trends
-      IF( l_trddyn )   THEN
-         ztdua(:,:,:) = ua(:,:,:) 
-         ztdva(:,:,:) = va(:,:,:) 
+      IF( l_trddyn ) THEN           ! Save ua and va trends
+         ztrdu(:,:,:) = ua(:,:,:) 
+         ztrdv(:,:,:) = va(:,:,:) 
       ENDIF
       
@@ -85,24 +76,17 @@
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
-         ! Horizontal kinetic energy at T-point
-         DO jj = 2, jpj
+         DO jj = 2, jpj         ! Horizontal kinetic energy at T-point
             DO ji = fs_2, jpi   ! vector opt.
+               zu = 0.25 * (  un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
+                  &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)  )
                zv = 0.25 * (  vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
-                            + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)  )
-               zu = 0.25 * (  un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
-                            + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)  )
+                  &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)  )
                zhke(ji,jj,jk) = zv + zu
             END DO  
          END DO  
-         
-         ! Horizontal gradient of Horizontal kinetic energy
-         DO jj = 2, jpjm1
+         DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ! gradient of kinetic energy
-               zua = -( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
-               zva = -( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
-               ! add to the general momentum trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
-               va(ji,jj,jk) = va(ji,jj,jk) + zva
+               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
+               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
             END DO 
          END DO
@@ -111,18 +95,13 @@
       !                                                ! ===============
 
-      ! save the Kinetic Energy trends for diagnostic
-      ! momentum trends
-      IF( l_trddyn )   THEN
-         ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:)
-         ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:)
-
-         CALL trd_mod(ztdua, ztdva, jpdtdkeg, 'DYN', kt)
+      IF( l_trddyn ) THEN      ! save the Kinetic Energy trends for diagnostic
+         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_keg, 'DYN', kt )
       ENDIF
-
-      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
-         CALL prt_ctl(tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask, &
-            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
-      ENDIF
-
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
+         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      !
    END SUBROUTINE dyn_keg
 
Index: trunk/NEMO/OPA_SRC/DYN/dynldf.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 474)
+++ trunk/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 503)
@@ -4,4 +4,6 @@
    !! Ocean physics:  lateral diffusivity trends 
    !!=====================================================================
+   !! History :  9.0  !  05-11  (G. Madec)  Original code (new step architecture)
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -9,31 +11,26 @@
    !!   dyn_ldf_ctl : initialization, namelist read, and parameters control
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE ldfdyn_oce      ! ocean dynamics lateral physics
-   USE ldfslp          ! ???
-   USE dynldf_bilapg   ! lateral mixing            (dyn_ldf_bilapg routine)
-   USE dynldf_bilap    ! lateral mixing             (dyn_ldf_bilap routine)
-   USE dynldf_iso      ! lateral mixing               (dyn_ldf_iso routine)
-   USE dynldf_lap      ! lateral mixing               (dyn_ldf_lap routine)
-   USE trdmod          ! ocean dynamics and tracer trends
-   USE trdmod_oce      ! ocean variables trends
-   USE prtctl          ! Print control
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distribued memory computing library
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE ldfdyn_oce     ! ocean dynamics lateral physics
+   USE ldfslp         ! lateral mixing: slopes of mixing orientation
+   USE dynldf_bilapg  ! lateral mixing            (dyn_ldf_bilapg routine)
+   USE dynldf_bilap   ! lateral mixing            (dyn_ldf_bilap  routine)
+   USE dynldf_iso     ! lateral mixing            (dyn_ldf_iso    routine)
+   USE dynldf_lap     ! lateral mixing            (dyn_ldf_lap    routine)
+   USE trdmod         ! ocean dynamics and tracer trends
+   USE trdmod_oce     ! ocean variables trends
+   USE prtctl         ! Print control
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! distribued memory computing library
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
 
    IMPLICIT NONE
    PRIVATE
 
-   !! *  Routine accessibility
-   PUBLIC dyn_ldf      ! called by step.F90 
+   PUBLIC   dyn_ldf   ! called by step module 
 
-   !! * module variables
-   INTEGER ::                        & 
-      nldf = 0                         ! type of lateral diffusion used
-      !                                ! defined from ln_dynldf_... namlist logicals)
+   INTEGER ::   nldf = 0   ! type of lateral diffusion used defined from ln_dynldf_... namlist logicals)
 
    !! * Substitutions
@@ -42,5 +39,7 @@
    !!---------------------------------------------------------------------------------
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
-   !!---------------------------------------------------------------------------------
+   !! $Header$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 
 CONTAINS
@@ -51,12 +50,8 @@
       !! 
       !! ** Purpose :   compute the lateral ocean dynamics physics.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !!
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-
-      !! * local declarations
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         ztrdu, ztrdv                         ! 3D temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrdu, ztrdv   ! 3D workspace
       !!----------------------------------------------------------------------
 
@@ -69,48 +64,34 @@
 
       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
+      !
+      CASE ( 0 )    ;   CALL dyn_ldf_lap    ( kt )      ! iso-level laplacian
+      CASE ( 1 )    ;   CALL dyn_ldf_iso    ( kt )      ! rotated laplacian (except dk[ dk[.] ] part)
+      CASE ( 2 )    ;   CALL dyn_ldf_bilap  ( kt )      ! iso-level bilaplacian
+      CASE ( 3 )    ;   CALL dyn_ldf_bilapg ( kt )      ! s-coord. horizontal bilaplacian
+      !
       CASE ( -1 )                                       ! esopa: test all possibility with control print
-         CALL dyn_ldf_lap    ( kt )
-         IF(ln_ctl) THEN         ! print sum trends (used for debugging)
-            CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf0 - Ua: ', mask1=umask, &
-               &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         ENDIF
-         CALL dyn_ldf_iso    ( kt )
-         IF(ln_ctl) THEN         ! print sum trends (used for debugging)
-            CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf1 - Ua: ', mask1=umask, &
-               &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         ENDIF
-         CALL dyn_ldf_bilap  ( kt )
-         IF(ln_ctl) THEN         ! print sum trends (used for debugging)
-            CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf2 - Ua: ', mask1=umask, &
-               &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         ENDIF
-         CALL dyn_ldf_bilapg ( kt )
-         IF(ln_ctl) THEN         ! print sum trends (used for debugging)
-            CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf3 - Ua: ', mask1=umask, &
-               &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         ENDIF
-
-      CASE ( 0 )                                       ! iso-level laplacian
-         CALL dyn_ldf_lap    ( kt )
-      CASE ( 1 )                                       ! rotated laplacian (except dk[ dk[.] ] part)
-         CALL dyn_ldf_iso    ( kt )
-      CASE ( 2 )                                       ! iso-level bilaplacian
-         CALL dyn_ldf_bilap  ( kt )
-      CASE ( 3 )                                       ! s-coord. horizontal bilaplacian
-         CALL dyn_ldf_bilapg ( kt )
+         ;              CALL dyn_ldf_lap    ( kt )
+         ;              CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf0 - Ua: ', mask1=umask,   &
+            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         ;              CALL dyn_ldf_iso    ( kt )
+         ;              CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf1 - Ua: ', mask1=umask,   &
+            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         ;              CALL dyn_ldf_bilap  ( kt )
+         ;              CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf2 - Ua: ', mask1=umask,   &
+            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         ;              CALL dyn_ldf_bilapg ( kt )
+         ;              CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf3 - Ua: ', mask1=umask,   &
+            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       END SELECT
 
-      !                                          ! save the horizontal diffusive trends for further diagnostics
-      IF( l_trddyn )   THEN
+      IF( l_trddyn ) THEN                        ! save the horizontal diffusive trends for further diagnostics
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-         CALL trd_mod( ztrdu, ztrdv, jpdtdldf, 'DYN', kt )
+         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_ldf, 'DYN', kt )
       ENDIF
-
-      IF(ln_ctl) THEN                            ! print mean trends (used for debugging)
-         CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf  - Ua: ', mask1=umask, &
-            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-      ENDIF
-
+      !                                          ! print sum trends (used for debugging)
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf  - Ua: ', mask1=umask,   &
+         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      !
    END SUBROUTINE dyn_ldf
 
@@ -121,32 +102,23 @@
       !! 
       !! ** Purpose :   initializations of the horizontal ocean dynamics physics
-      !!
-      !! ** Method :
-      !!
-      !! History :
-      !!   9.0  !  05-11  (G. Madec)  Original code
       !!----------------------------------------------------------------------
-      !! * Local declarations
       INTEGER ::   ioptio, ierr         ! temporary integers 
       !!----------------------------------------------------------------------
+    
+      !                                   ! Namelist nam_dynldf: already read in ldfdyn module
 
-      !  Define the lateral dynamics physics parameters
-      ! =============================================
-    
-      ! Namelist nam_dynldf already read in ldfdyn module
-
-      IF(lwp) THEN
+      IF(lwp) THEN                        ! Namelist print
          WRITE(numout,*)
-         WRITE(numout,*) 'dyn:ldf_ctl : Choice of the lateral diffusive operator on dynamics'
+         WRITE(numout,*) 'dyn_ldf_ctl : Choice of the lateral diffusive operator on dynamics'
          WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*) '          Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)'
-         WRITE(numout,*) '             laplacian operator          ln_dynldf_lap   = ', ln_dynldf_lap
-         WRITE(numout,*) '             bilaplacian operator        ln_dynldf_bilap = ', ln_dynldf_bilap
-         WRITE(numout,*) '             iso-level                   ln_dynldf_level = ', ln_dynldf_level
-         WRITE(numout,*) '             horizontal (geopotential)   ln_dynldf_hor   = ', ln_dynldf_hor
-         WRITE(numout,*) '             iso-neutral                 ln_dynldf_iso   = ', ln_dynldf_iso
+         WRITE(numout,*) '       Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)'
+         WRITE(numout,*) '          laplacian operator          ln_dynldf_lap   = ', ln_dynldf_lap
+         WRITE(numout,*) '          bilaplacian operator        ln_dynldf_bilap = ', ln_dynldf_bilap
+         WRITE(numout,*) '          iso-level                   ln_dynldf_level = ', ln_dynldf_level
+         WRITE(numout,*) '          horizontal (geopotential)   ln_dynldf_hor   = ', ln_dynldf_hor
+         WRITE(numout,*) '          iso-neutral                 ln_dynldf_iso   = ', ln_dynldf_iso
       ENDIF
 
-      ! control the consistency
+      !                                   ! control the consistency
       ioptio = 0
       IF( ln_dynldf_lap   )   ioptio = ioptio + 1
@@ -159,5 +131,5 @@
       IF( ioptio /= 1 ) CALL ctl_stop( '          use only ONE direction (level/hor/iso)' )
 
-      ! defined the type of lateral diffusion from ln_dynldf_... logicals
+      !                                   ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals
       ierr = 0
       IF ( ln_dynldf_lap ) THEN      ! laplacian operator
@@ -196,17 +168,11 @@
          ENDIF
       ENDIF
+      
+      IF( lk_esopa )                 nldf = -1     ! esopa test
 
-      IF( ierr == 1 )   &
-           &   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' )
-      IF( ierr == 2 )   &
-           &   CALL ctl_stop( ' isoneutral bilaplacian operator does not exist' )
+      IF( ierr == 1 )   CALL ctl_stop( 'iso-level in z-coordinate - partial step, not allowed' )
+      IF( ierr == 2 )   CALL ctl_stop( 'isoneutral bilaplacian operator does not exist' )
       IF( nldf == 1 .OR. nldf == 3 ) THEN      ! rotation
-         IF( .NOT.lk_ldfslp )   &
-           &   CALL ctl_stop(  '          the rotation of the diffusive tensor require key_ldfslp' )
-      ENDIF
-
-      IF( lk_esopa ) THEN
-         IF(lwp ) WRITE(numout,*) '          esopa test: use all lateral physics options'
-         nldf = -1
+         IF( .NOT.lk_ldfslp )   CALL ctl_stop( 'the rotation of the diffusive tensor require key_ldfslp' )
       ENDIF
 
@@ -219,5 +185,5 @@
          IF( nldf ==  3 )   WRITE(numout,*) '              Rotated bilaplacian'
       ENDIF
-
+      !
    END SUBROUTINE dyn_ldf_ctl
 
Index: trunk/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 474)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 503)
@@ -4,4 +4,7 @@
    !! Ocean dynamics:  surface pressure gradient control
    !!======================================================================
+   !! History :  9.0  !  05-12  (C. Talandier, G. Madec)  Original code
+   !!            9.0  !  05-12  (V. Garnier)  dyn_spg_ctl: Original code
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -9,5 +12,4 @@
    !!   dyn_spg_ctl : initialization, namelist read, and parameters control
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce            ! ocean dynamics and tracers variables
    USE dom_oce        ! ocean space and time domain variables
@@ -29,11 +31,8 @@
    PRIVATE
 
-   !! * Accessibility
    PUBLIC dyn_spg         ! routine called by step module
 
    !! * module variables
-   INTEGER ::                        &
-      nspg = 0                         ! type of surface pressure gradient scheme
-      !                                ! defined from lk_dynspg_... 
+   INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_... 
 
    !! * Substitutions
@@ -43,5 +42,5 @@
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! 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)
    !!----------------------------------------------------------------------
 
@@ -53,16 +52,10 @@
       !!
       !! ** Purpose :   compute the lateral ocean dynamics physics.
-      !!
-      !! History :
-      !!   9.0  !  05-12  (C. Talandier, G. Madec)  Original code
-      !!----------------------------------------------------------------------
-      !! * Arguments
+      !!----------------------------------------------------------------------
       INTEGER, INTENT( in  ) ::   kt     ! ocean time-step index
       INTEGER, INTENT( out ) ::   kindic ! solver flag
-
-      !! * local declarations
+      !!
       REAL(wp) ::   z2dt                      ! temporary scalar
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         ztrdu, ztrdv                         ! 3D temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrdu, ztrdv   ! 3D workspace
       !!----------------------------------------------------------------------
 
@@ -75,42 +68,36 @@
 
       SELECT CASE ( nspg )                       ! compute surf. pressure gradient trend and add it to the general trend
+      !                                                      ! k-j-i loops
+      CASE (  0 )   ;   CALL dyn_spg_exp    ( kt )              ! explicit
+      CASE (  1 )   ;   CALL dyn_spg_ts     ( kt )              ! time-splitting
+      CASE (  2 )   ;   CALL dyn_spg_flt    ( kt, kindic )      ! filtered
+      CASE (  3 )   ;   CALL dyn_spg_rl     ( kt, kindic )      ! rigid lid
+      !                                                      ! j-k-i loops
+      CASE ( 10 )   ;   CALL dyn_spg_exp_jki( kt )              ! explicit with j-k-i loop
+      CASE ( 11 )   ;   CALL dyn_spg_ts_jki ( kt )              ! time-splitting with j-k-i loop
+      CASE ( 12 )   ;   CALL dyn_spg_flt_jki( kt, kindic )      ! filtered with j-k-i loop
+      !
       CASE ( -1 )                                       ! esopa: test all possibility with control print
-         CALL dyn_spg_exp    ( kt )
-         IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' spg0 - Ua: ', mask1=umask, &
-               &                    tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL dyn_spg_ts     ( kt )
-         IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' spg1 - Ua: ', mask1=umask, &
-               &                    tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL dyn_spg_flt  ( kt, kindic )
-         IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' spg2 - Ua: ', mask1=umask, &
-               &                    tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL dyn_spg_exp_jki( kt )
-         IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' spg10- Ua: ', mask1=umask, &
-               &                    tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL dyn_spg_ts_jki ( kt )
-         IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' spg12- Ua: ', mask1=umask, &
-               &                    tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL dyn_spg_flt_jki( kt, kindic )
-         IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' spg13- Ua: ', mask1=umask, &
-               &                    tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-      CASE (  0 )                                       ! explicit
-         CALL dyn_spg_exp    ( kt )
-      CASE (  1 )                                       ! time-splitting
-         CALL dyn_spg_ts     ( kt )
-      CASE (  2 )                                       ! filtered
-         CALL dyn_spg_flt    ( kt, kindic )
-      CASE (  3 )                                       ! rigid lid
-         CALL dyn_spg_rl     ( kt, kindic )
-
-      CASE ( 10 )                                       ! explicit with j-k-i loop
-         CALL dyn_spg_exp_jki( kt )
-      CASE ( 11 )                                       ! time-splitting with j-k-i loop
-         CALL dyn_spg_ts_jki ( kt )
-      CASE ( 12 )                                       ! filtered with j-k-i loop
-         CALL dyn_spg_flt_jki( kt, kindic )
+         ;              CALL dyn_spg_exp    ( kt )
+         ;              CALL prt_ctl( tab3d_1=ua, clinfo1=' spg0 - Ua: ', mask1=umask, &
+            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         ;              CALL dyn_spg_ts     ( kt )
+         ;              CALL prt_ctl( tab3d_1=ua, clinfo1=' spg1 - Ua: ', mask1=umask, &
+            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         ;              CALL dyn_spg_flt  ( kt, kindic )
+         ;              CALL prt_ctl( tab3d_1=ua, clinfo1=' spg2 - Ua: ', mask1=umask, &
+            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         ;              CALL dyn_spg_exp_jki( kt )
+         ;              CALL prt_ctl( tab3d_1=ua, clinfo1=' spg10- Ua: ', mask1=umask, &
+            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         ;              CALL dyn_spg_ts_jki ( kt )
+         ;              CALL prt_ctl( tab3d_1=ua, clinfo1=' spg12- Ua: ', mask1=umask, &
+            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         ;              CALL dyn_spg_flt_jki( kt, kindic )
+         ;              CALL prt_ctl( tab3d_1=ua, clinfo1=' spg13- Ua: ', mask1=umask, &
+            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       END SELECT
-
-      !                                          ! save the horizontal diffusive trends for further diagnostics
-      IF( l_trddyn )   THEN
+      !                    
+      IF( l_trddyn )   THEN                      ! save the horizontal diffusive trends for further diagnostics
          SELECT CASE ( nspg )
          CASE ( 0, 1, 3, 10, 11 )
@@ -123,11 +110,10 @@
             ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:)
          END SELECT
-         CALL trd_mod( ztrdu, ztrdv, jpdtdspg, 'DYN', kt )
-      ENDIF
-
+         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_spg, 'DYN', kt )
+      ENDIF
       !                                          ! print mean trends (used for debugging)
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' spg  - Ua: ', mask1=umask, &
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-
+      !
    END SUBROUTINE dyn_spg
 
@@ -139,7 +125,4 @@
       !! ** Purpose :   Control the consistency between cpp options for 
       !!      surface pressure gradient schemes
-      !!
-      !! History :
-      !!   9.0  !  05-10  (V. Garnier)  Original code : spg re-organization
       !!----------------------------------------------------------------------
       !! * Local declarations
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90	(revision 474)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90	(revision 503)
@@ -9,13 +9,9 @@
    !!   'key_mpp_omp'                              j-k-i loop (vector opt.)
    !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
    !!   dyn_spg_flt_jki : Update the momentum trend with the surface pressure
    !!                     gradient for the free surf. constant volume case
    !!                     with auto-tasking (j-slab) (no vectior opt.)
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
-   !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers 
    USE dom_oce         ! ocean space and time domain
@@ -51,5 +47,5 @@
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! 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)
    !!----------------------------------------------------------------------
 
@@ -100,21 +96,9 @@
       !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
       !!
-      !! References :
-      !!      Roullet and Madec 1999, JGR.
-      !!
-      !! History :
-      !!        !  98-05 (G. Roullet)  Original code
-      !!        !  98-10 (G. Madec, M. Imbard)  release 8.2
-      !!   8.5  !  02-08 (G. Madec)  F90: Free form and module
-      !!        !  02-11 (C. Talandier, A-M Treguier) Open boundaries
-      !!   9.0  !  04-08 (C. Talandier) New trends organization
-      !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
+      !! References : Roullet and Madec 1999, JGR.
       !!---------------------------------------------------------------------
-      !! * Arguments
       INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
-      INTEGER, INTENT( out ) ::   kindic     ! solver convergence flag
-                                             ! if the solver doesn t converge
-                                             ! the flag is < 0
-      !! * Local declarations
+      INTEGER, INTENT( out ) ::   kindic     ! solver convergence flag, <0 if the solver doesn t converge
+      !!
       INTEGER  ::   ji, jj, jk               ! dummy loop indices
       REAL(wp) ::   &              ! temporary scalars
@@ -154,5 +138,4 @@
             spgv(ji,jj) = - grav * ( sshn(ji  ,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
          END DO 
-
          !  Add the surface pressure trend to the general trend
          DO jk = 1, jpkm1
@@ -162,5 +145,4 @@
             END DO
          END DO
-
          ! Evaluate the masked next velocity (effect of the additional force not included)
          DO jk = 1, jpkm1
@@ -170,5 +152,4 @@
             END DO
          END DO
-
          !                                             ! ===============
       END DO                                           !   End of slab
@@ -202,5 +183,4 @@
             spgv(ji,jj) = 0.e0
          END DO
-
          ! vertical sum
          DO jk = 1, jpk
@@ -210,5 +190,4 @@
             END DO
          END DO
-
          ! transport: multiplied by the horizontal scale factor
          DO ji = 2, jpim1
@@ -216,5 +195,4 @@
             spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
          END DO
-
          !                                             ! ===============
       END DO                                           !   End of slab
@@ -302,10 +280,9 @@
             WRITE(ctmp1,*) ' ~~~~~~~~~~~~~~~~                not = ', nsolv
             CALL ctl_stop( ' dyn_spg_flt_jki : e r r o r, nsolv = 1, 2, 3 or 4', ctmp1 )
-      ENDIF
-
-      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-
-!CDIR PARALLEL DO
-!$OMP PARALLEL DO
+         ENDIF
+      ENDIF
+
+      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
+
       !                                                ! ===============
       DO jj = 2, jpjm1                                 !  Vertical slab
@@ -353,6 +330,5 @@
          ! 8. Sea surface elevation time stepping
          ! --------------------------------------
-         ! Euler (forward) time stepping, no time filter
-         IF( neuler == 0 .AND. kt == nit000 ) THEN
+         IF( neuler == 0 .AND. kt == nit000 ) THEN         ! Euler (forward) time stepping, no time filter
             DO ji = 1, jpi
                ! after free surface elevation
@@ -362,6 +338,5 @@
                sshn(ji,jj) = zssha
             END DO
-         ELSE
-            ! Leap-frog time stepping and time filter
+         ELSE                                              ! Leap-frog time stepping and time filter
             DO ji = 1, jpi
                ! after free surface elevation
@@ -384,6 +359,5 @@
          CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg  - ssh: ', mask1=tmask)
       ENDIF
-
-
+      !
    END SUBROUTINE dyn_spg_flt_jki
 
Index: trunk/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 474)
+++ trunk/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 503)
@@ -5,4 +5,12 @@
    !!                 planetary vorticity trends
    !!======================================================================
+   !! History :  1.0  !  89-12  (P. Andrich)  vor_ens: Original code
+   !!            5.0  !  91-11  (G. Madec) vor_ene, vor_mix: Original code
+   !!            6.0  !  96-01  (G. Madec)  s-coord, suppress work arrays
+   !!            8.5  !  02-08  (G. Madec)  F90: Free form and module
+   !!            8.5  !  04-02  (G. Madec)  vor_een: Original code
+   !!            9.0  !  03-08  (G. Madec)  vor_ctl: Original code
+   !!            9.0  !  05-11  (G. Madec)  dyn_vor: Original code (new step architecture)
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -14,28 +22,25 @@
    !!       vor_ctl : control of the different vorticity option
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE in_out_manager  ! I/O manager
-   USE trdmod          ! ocean dynamics trends 
-   USE trdmod_oce      ! ocean variables trends
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE prtctl          ! Print control
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE trdmod         ! ocean dynamics trends 
+   USE trdmod_oce     ! ocean variables trends
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl         ! Print control
+   USE in_out_manager ! I/O manager
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dyn_vor      ! routine called by step.F90
-
-   !! * Shared module variables
+   PUBLIC   dyn_vor   ! routine called by step.F90
+
+   !!* Namelist nam_dynvor: vorticity term
    LOGICAL, PUBLIC ::   ln_dynvor_ene = .FALSE.   !: energy conserving scheme
    LOGICAL, PUBLIC ::   ln_dynvor_ens = .TRUE.    !: enstrophy conserving scheme
    LOGICAL, PUBLIC ::   ln_dynvor_mix = .FALSE.   !: mixed scheme
-   LOGICAL, PUBLIC ::   ln_dynvor_een = .FALSE.    !: energy and enstrophy conserving scheme
-
-   !! * module variables
-   INTEGER ::                        &
-      nvor = 0                         ! type of vorticity trend used
+   LOGICAL, PUBLIC ::   ln_dynvor_een = .FALSE.   !: energy and enstrophy conserving scheme
+   NAMELIST/nam_dynvor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
+
+   INTEGER ::   nvor = 0   ! type of vorticity trend used
 
    !! * Substitutions
@@ -43,5 +48,7 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
+   !! $Header$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -54,15 +61,11 @@
       !!
       !! ** Action : - Update (ua,va) with the now vorticity term trend
-      !!             - save the trends in (utrd,vtrd) in 2 parts (relative
+      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
       !!               and planetary vorticity trends) ('key_trddyn')
-      !!
-      !! History :
-      !!   9.0  !  05-11  (G. Madec)  Original code
-      !!----------------------------------------------------------------------
-      USE oce, ONLY :    ztrdu => ta,     & ! use ta as 3D workspace
-                         ztrdv => sa        ! use sa as 3D workspace
-
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
+      !!----------------------------------------------------------------------
+      USE oce, ONLY :   ztrdu => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   ztrdv => sa   ! use sa as 3D workspace
+      !!
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
       !!----------------------------------------------------------------------
 
@@ -74,15 +77,15 @@
       CASE ( -1 )                                      ! esopa: test all possibility with control print
          CALL vor_ene( kt, 'TOT', ua, va )
-         IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
-            &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
+            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
          CALL vor_ens( kt, 'TOT', ua, va )
-         IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
-            &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
+            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
          CALL vor_mix( kt )
-         IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
-            &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
+            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
          CALL vor_een( kt, 'TOT', ua, va )
-         IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
-            &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
+            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
 
       CASE ( 0 )                                       ! energy conserving scheme
@@ -93,5 +96,5 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt )
+            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
@@ -99,6 +102,6 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt )
-            CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt )
+            CALL trd_mod( ztrdu, ztrdu, jpdyn_trd_pvo, 'DYN', kt )
+            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
          ELSE
             CALL vor_ene( kt, 'TOT', ua, va )                ! total vorticity
@@ -112,5 +115,5 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt )
+            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
@@ -118,6 +121,6 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt )
-            CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt )
+            CALL trd_mod( ztrdu, ztrdu, jpdyn_trd_pvo, 'DYN', kt )
+            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
          ELSE
             CALL vor_ens( kt, 'TOT', ua, va )                ! total vorticity
@@ -131,5 +134,5 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt )
+            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
@@ -137,6 +140,6 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt )
-            CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt )
+            CALL trd_mod( ztrdu, ztrdu, jpdyn_trd_pvo, 'DYN', kt )
+            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
          ELSE
             CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene)
@@ -150,5 +153,5 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt )
+            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
@@ -156,6 +159,6 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt )
-            CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt )
+            CALL trd_mod( ztrdu, ztrdu, jpdyn_trd_pvo, 'DYN', kt )
+            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
          ELSE
             CALL vor_een( kt, 'TOT', ua, va )                ! total vorticity
@@ -192,30 +195,19 @@
       !!
       !! ** Action : - Update (ua,va) with the now vorticity term trend
-      !!             - save the trends in (utrd,vtrd) in 2 parts (relative
+      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
       !!               and planetary vorticity trends) ('key_trddyn')
       !!
-      !! References :
-      !!      Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
-      !! History :
-      !!   5.0  !  91-11  (G. Madec)  Original code
-      !!   6.0  !  96-01  (G. Madec)  s-coord, suppress work arrays
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
-      !!   9.0  !  04-08  (C. Talandier)  New trends organization
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
-      CHARACTER(len=3) , INTENT( in ) ::   &
-         cd_vor        ! define the vorticity considered
-         !             ! ='COR' (planetary) ; ='VOR' (relative) ; ='TOT' (total)
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   &
-         pua, pva                      ! ????
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk               ! dummy loop indices
-      REAL(wp) ::   &
-         zfact2,                         &  ! temporary scalars
-         zx1, zx2, zy1, zy2                 !    "         "
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zwx, zwy, zwz                      ! temporary workspace
+      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
+      !!----------------------------------------------------------------------
+      INTEGER          , INTENT(in   )                         ::   kt       ! ocean time-step index
+      CHARACTER(len=3) , INTENT(in   )                         ::   cd_vor   ! ='COR' (planetary) ; ='VOR' (relative)
+         !                                                                   ! ='TOT' (total)
+      REAL(wp)         , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua      ! total u-trend
+      REAL(wp)         , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva      ! total v-trend
+      !!
+      INTEGER  ::   ji, jj, jk         ! dummy loop indices
+      REAL(wp) ::   zx1, zy1, zfact2   ! temporary scalars
+      REAL(wp) ::   zx2, zy2           !    "         "
+      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! temporary 2D workspace
       !!----------------------------------------------------------------------
 
@@ -234,14 +226,10 @@
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
-         
          ! Potential vorticity and horizontal fluxes
          ! -----------------------------------------
          SELECT CASE( cd_vor )      ! vorticity considered
-         CASE ( 'COR' )                   ! planetary vorticcity (Coriolis)
-            zwz(:,:) = ff(:,:)
-         CASE ( 'VOR' )                   ! relative  vorticity
-            zwz(:,:) = rotn(:,:,jk)
-         CASE ( 'TOT' )                   ! total     vorticity (planetary + relative)
-            zwz(:,:) = rotn(:,:,jk) + ff(:,:)
+         CASE ( 'COR' )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
+         CASE ( 'VOR' )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
+         CASE ( 'TOT' )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total     vorticity 
          END SELECT
 
@@ -270,5 +258,4 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-
    END SUBROUTINE vor_ene
 
@@ -300,25 +287,15 @@
       !!
       !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
-      !!             - Save the trends in (utrd,vtrd) in 2 parts (relative
+      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
       !!               and planetary vorticity trends) ('key_trddyn')
       !!
-      !! References :
-      !!      Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
-      !! History :
-      !!   5.0  !  91-11  (G. Madec) Original code, enstrophy-energy-combined schemes
-      !!   6.0  !  96-01  (G. Madec) s-coord, suppress work arrays
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
-      !!   9.0  !  04-08  (C. Talandier)  New trends organization
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt         ! ocean timestep index
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk               ! dummy loop indices
-      REAL(wp) ::   &
-         zfact1, zfact2, zua, zva,       &  ! temporary scalars
-         zcua, zcva, zx1, zx2, zy1, zy2
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zwx, zwy, zwz, zww                 ! temporary workspace
+      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean timestep index
+      !!
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! temporary scalars
+      REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !    "         "
+      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! temporary 3D workspace
       !!----------------------------------------------------------------------
 
@@ -367,5 +344,5 @@
                zcua = zfact2 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
                zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
-
+               ! mixed vorticity trend added to the momentum trends
                ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
                va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
@@ -375,5 +352,4 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-
    END SUBROUTINE vor_mix
 
@@ -400,29 +376,18 @@
       !!
       !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
-      !!             - Save the trends in (utrd,vtrd) in 2 parts (relative 
+      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative 
       !!               and planetary vorticity trends) ('key_trddyn')
       !!
-      !! References :
-      !!      Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
-      !! History :
-      !!   5.0  !  91-11  (G. Madec)  Original code
-      !!   6.0  !  96-01  (G. Madec)  s-coord, suppress work arrays
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
-      !!   9.0  !  04-08  (C. Talandier)  New trends organization
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt    ! ocean timestep
-      CHARACTER(len=3) , INTENT( in ) ::   &
-         cd_vor        ! define the vorticity considered
-         !             ! ='COR' (planetary) ; ='VOR' (relative) ; ='TOT' (total)
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   &
-         pua, pva                      ! ????
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk          ! dummy loop indices
-      REAL(wp) ::   &
-         zfact1, zuav, zvau            ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zwx, zwy, zwz                 ! temporary workspace
+      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
+      !!----------------------------------------------------------------------
+      INTEGER          , INTENT(in   )                         ::   kt       ! ocean time-step index
+      CHARACTER(len=3) , INTENT(in   )                         ::   cd_vor   ! ='COR' (planetary) ; ='VOR' (relative)
+         !                                                                   ! ='TOT' (total)
+      REAL(wp)         , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua      ! total u-trend
+      REAL(wp)         , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva      ! total v-trend
+      !!
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! temporary 3D workspace
       !!----------------------------------------------------------------------
       
@@ -441,14 +406,10 @@
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
-
          ! Potential vorticity and horizontal fluxes
          ! -----------------------------------------
          SELECT CASE( cd_vor )      ! vorticity considered
-         CASE ( 'COR' )                   ! planetary vorticcity (Coriolis)
-            zwz(:,:) = ff(:,:)
-         CASE ( 'VOR' )                   ! relative  vorticity
-            zwz(:,:) = rotn(:,:,jk)
-         CASE ( 'TOT' )                   ! total     vorticity (planetary + relative)
-            zwz(:,:) = rotn(:,:,jk) + ff(:,:)
+         CASE ( 'COR' )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
+         CASE ( 'VOR' )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
+         CASE ( 'TOT' )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total     vorticity 
          END SELECT
 
@@ -470,5 +431,4 @@
          ENDIF
 
-
          ! Compute and add the vorticity term trend
          ! ----------------------------------------
@@ -476,8 +436,7 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
-                                            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
+                  &                         + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
                zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
-                                            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
-
+                  &                         + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
                pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
                pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
@@ -487,5 +446,4 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-
    END SUBROUTINE vor_ens
 
@@ -509,33 +467,20 @@
       !!
       !! ** Action : - Update (ua,va) with the now vorticity term trend
-      !!             - save the trends in (utrd,vtrd) in 2 parts (relative
+      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
       !!               and planetary vorticity trends) ('key_trddyn')
       !!
-      !! References :
-      !!      Arakawa and Lamb 1980, A potential enstrophy and energy conserving
-      !!                             scheme for the Shallow water equations, 
-      !!                             Monthly Weather Review, vol. 109, p 18-36
-      !!
-      !! History :
-      !!   8.5  !  04-02  (G. Madec)  Original code
-      !!   9.0  !  04-08  (C. Talandier)  New trends organization
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
-      CHARACTER(len=3) , INTENT( in ) ::   &
-         cd_vor        ! define the vorticity considered
-         !             ! ='COR' (planetary) ; ='VOR' (relative) ; ='TOT' (total)
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   &
-         pua, pva                      ! ????
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk              ! dummy loop indices
-      REAL(wp) ::   &
-         zfac12, zua, zva                  ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zwx, zwy, zwz,                 &  ! temporary workspace
-         ztnw, ztne, ztsw, ztse            !    "           "
-      REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE ::   &
-         ze3f
+      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
+      !!----------------------------------------------------------------------
+      INTEGER          , INTENT(in   )                         ::   kt       ! ocean time-step index
+      CHARACTER(len=3) , INTENT(in   )                         ::   cd_vor   ! ='COR' (planetary) ; ='VOR' (relative)
+         !                                                                   ! ='TOT' (total)
+      REAL(wp)         , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua      ! total u-trend
+      REAL(wp)         , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva      ! total v-trend
+      !!
+      INTEGER ::   ji, jj, jk          ! dummy loop indices
+      REAL(wp) ::   zfac12, zua, zva   ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz            ! temporary 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztnw, ztne, ztsw, ztse   ! temporary 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE ::   ze3f
       !!----------------------------------------------------------------------
 
@@ -570,10 +515,7 @@
          ! -----------------------------------------
          SELECT CASE( cd_vor )      ! vorticity considered
-         CASE ( 'COR' )                   ! planetary vorticcity (Coriolis)
-            zwz(:,:) = ff(:,:) * ze3f(:,:,jk)
-         CASE ( 'VOR' )                   ! relative  vorticity
-            zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
-         CASE ( 'TOT' )                   ! total     vorticity (planetary + relative)
-            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
+         CASE ( 'COR' )   ;   zwz(:,:) =                  ff(:,:)   * ze3f(:,:,jk)   ! planetary vorticity (Coriolis)
+         CASE ( 'VOR' )   ;   zwz(:,:) =   rotn(:,:,jk)             * ze3f(:,:,jk)   ! relative  vorticity
+         CASE ( 'TOT' )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)   ! total     vorticity 
          END SELECT
 
@@ -599,5 +541,4 @@
             END DO
          END DO
-         
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -613,5 +554,4 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-
    END SUBROUTINE vor_een
 
@@ -623,67 +563,49 @@
       !! ** Purpose :   Control the consistency between cpp options for
       !!      tracer advection schemes
-      !!
-      !! History :
-      !!   9.0  !  03-08  (G. Madec)  Original code
-      !!----------------------------------------------------------------------
-      !! * Local declarations
+      !!----------------------------------------------------------------------
       INTEGER ::   ioptio          ! temporary integer
-
-      NAMELIST/nam_dynvor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
-      !!----------------------------------------------------------------------
-
-      ! Read Namelist nam_dynvor : Vorticity scheme options
-      ! ------------------------
-      REWIND ( numnam )
+      !!----------------------------------------------------------------------
+
+      REWIND ( numnam )               ! Read Namelist nam_dynvor : Vorticity scheme options
       READ   ( numnam, nam_dynvor )
 
-      ! Control of vorticity scheme options
-      ! -----------------------------------
-      ! Control print
-      IF(lwp) THEN
+      IF(lwp) THEN                    ! Namelist print
          WRITE(numout,*)
          WRITE(numout,*) 'dyn:vor_ctl : vorticity term : read namelist and control the consistency'
          WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*) '              Namelist nam_dynvor : oice of the vorticity term scheme'
-         WRITE(numout,*) '                 energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
-         WRITE(numout,*) '                 enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
-         WRITE(numout,*) '                 mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
-         WRITE(numout,*) '                 enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
+         WRITE(numout,*) '        Namelist nam_dynvor : oice of the vorticity term scheme'
+         WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
+         WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
+         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
+         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
       ENDIF
 
-      ioptio = 0
-      IF( ln_dynvor_ene ) THEN
-         nvor = 0
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '              vorticity term : energy conserving scheme'
-         ioptio = ioptio + 1
+      ioptio = 0                     ! Control of vorticity scheme options
+      IF( ln_dynvor_ene )   ioptio = ioptio + 1
+      IF( ln_dynvor_ens )   ioptio = ioptio + 1
+      IF( ln_dynvor_mix )   ioptio = ioptio + 1
+      IF( ln_dynvor_een )   ioptio = ioptio + 1
+      IF( lk_esopa      )   ioptio =          1
+
+      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
+
+      !                              ! Set nvor
+      IF( ln_dynvor_ene )   nvor =  0
+      IF( ln_dynvor_ens )   nvor =  1
+      IF( ln_dynvor_mix )   nvor =  2
+      IF( ln_dynvor_een )   nvor =  3
+      IF( lk_esopa      )   nvor = -1
+      
+      IF(lwp) THEN                   ! Print the choice
+         WRITE(numout,*)
+         IF( nvor ==  0 )   WRITE(numout,*) '         vorticity term used : energy conserving scheme'
+         IF( nvor ==  1 )   WRITE(numout,*) '         vorticity term used : enstrophy conserving scheme'
+         IF( nvor ==  2 )   WRITE(numout,*) '         vorticity term used : mixed enstrophy/energy conserving scheme'
+         IF( nvor ==  3 )   WRITE(numout,*) '         vorticity term used : energy and enstrophy conserving scheme'
+         IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options'
       ENDIF
-      IF( ln_dynvor_ens ) THEN
-         nvor = 1
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '              vorticity term : enstrophy conserving scheme'
-         ioptio = ioptio + 1
-      ENDIF
-      IF( ln_dynvor_mix ) THEN
-         nvor = 2
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '              vorticity term : mixed enstrophy/energy conserving scheme'
-         ioptio = ioptio + 1
-      ENDIF
-      IF( ln_dynvor_een ) THEN
-         nvor = 3
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '              vorticity term : energy and enstrophy conserving scheme'
-         ioptio = ioptio + 1
-      ENDIF
-      IF( ioptio /= 1 .AND. .NOT. lk_esopa ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
-      IF( lk_esopa ) THEN
-         nvor = -1
-         IF(lwp ) WRITE(numout,*) '          esopa test: use all lateral physics options'
-      ENDIF
-      IF(lwp) WRITE(numout,*) '             choice of vor_...                nvor   = ', nvor
-
+      !
    END SUBROUTINE vor_ctl
 
-!!==============================================================================
+   !!==============================================================================
 END MODULE dynvor
Index: trunk/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 474)
+++ trunk/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 503)
@@ -4,22 +4,26 @@
    !! Ocean dynamics : vertical advection trend
    !!======================================================================
+   !! History :  6.0  !  91-01  (G. Madec) Original code
+   !!            7.0  !  91-11  (G. Madec)
+   !!            7.5  !  96-01  (G. Madec) statement function for e3
+   !!            8.5  !  02-07  (G. Madec) j-k-i case: Original code
+   !!            8.5  !  02-07  (G. Madec) Free form, F90
+   !!----------------------------------------------------------------------
    
    !!----------------------------------------------------------------------
-   !!   dyn_zad      : vertical advection momentum trend
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE in_out_manager  ! I/O manager
-   USE trdmod          ! ocean dynamics trends 
-   USE trdmod_oce      ! ocean variables trends
-   USE flxrnf          ! ocean runoffs
-   USE prtctl          ! Print control
+   !!   dyn_zad       : vertical advection momentum trend
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE in_out_manager ! I/O manager
+   USE trdmod         ! ocean dynamics trends 
+   USE trdmod_oce     ! ocean variables trends
+   USE flxrnf         ! ocean runoffs
+   USE prtctl         ! Print control
 
    IMPLICIT NONE
    PRIVATE
    
-   !! * Accessibility
-   PUBLIC dyn_zad                ! routine called by step.F90
+   PUBLIC   dyn_zad   ! routine called by step.F90
 
    !! * Substitutions
@@ -29,5 +33,5 @@
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! 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)
    !!----------------------------------------------------------------------
 
@@ -54,27 +58,15 @@
       !!
       !! ** Action  : - Update (ua,va) with the vert. momentum advection trends
-      !!              - Save the trends in (utrd,vtrd) ('key_trddyn')
-      !!
-      !! History :
-      !!   6.0  !  91-01  (G. Madec) Original code
-      !!   7.0  !  91-11  (G. Madec)
-      !!   7.5  !  96-01  (G. Madec) statement function for e3
-      !!   8.5  !  02-07  (G. Madec) Free form, F90
-      !!   9.0  !  04-08  (C. Talandier) New trends organization
-      !!----------------------------------------------------------------------
-      !! * modules used
-      USE oce, ONLY:   zwuw => ta,   & ! use ta as 3D workspace
-                       zwvw => sa      ! use sa as 3D workspace
-
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt    ! ocean time-step inedx
-      
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk         ! dummy loop indices
-      REAL(wp) ::   zvn, zua, zva      ! temporary scalars
-      REAL(wp), DIMENSION(jpi) ::   &
-         zww                           ! temporary workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         ztdua, ztdva                  ! temporary workspace
+      !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn')
+      !!----------------------------------------------------------------------
+      USE oce, ONLY:   zwuw => ta   ! use ta as 3D workspace
+      USE oce, ONLY:   zwvw => sa   ! use sa as 3D workspace
+      !!
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
+      !!
+      INTEGER  ::   ji, jj, jk      ! dummy loop indices
+      REAL(wp) ::   zvn, zua, zva   ! temporary scalars
+      REAL(wp), DIMENSION(jpi)         ::   zww            ! 1D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrdu, ztrdv   ! 3D workspace
       !!----------------------------------------------------------------------
       
@@ -85,8 +77,7 @@
       ENDIF
 
-      ! Save ua and va trends
-      IF( l_trddyn )   THEN
-         ztdua(:,:,:) = ua(:,:,:) 
-         ztdva(:,:,:) = va(:,:,:) 
+      IF( l_trddyn )   THEN         ! Save ua and va trends
+         ztrdu(:,:,:) = ua(:,:,:) 
+         ztrdv(:,:,:) = va(:,:,:) 
       ENDIF
 
@@ -94,14 +85,9 @@
       DO jj = 2, jpjm1                                 !  Vertical slab
          !                                             ! ===============
-
-         ! Vertical momentum advection at level w and u- and v- vertical
-         ! ----------------------------------------------------------------
-         DO jk = 2, jpkm1
-            ! vertical fluxes 
-            DO ji = 2, jpi
+         DO jk = 2, jpkm1         ! Vertical momentum advection at uw and vw-pts
+            DO ji = 2, jpi              ! vertical fluxes 
                zww(ji) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
             END DO
-            ! vertical momentum advection at w-point
-            DO ji = 2, jpim1
+            DO ji = 2, jpim1            ! vertical momentum advection at w-point
                zvn = 0.25 * e1t(ji,jj+1) * e2t(ji,jj+1) * wn(ji,jj+1,jk)
                zwuw(ji,jj,jk) = ( zww(ji+1) + zww(ji) ) * ( un(ji,jj,jk-1)-un(ji,jj,jk) )
@@ -109,7 +95,5 @@
             END DO  
          END DO   
-
-         ! Surface and bottom values set to zero
-         DO ji = 2, jpim1
+         DO ji = 2, jpim1               ! Surface and bottom values set to zero
             zwuw(ji,jj, 1 ) = 0.e0
             zwvw(ji,jj, 1 ) = 0.e0
@@ -117,13 +101,11 @@
             zwvw(ji,jj,jpk) = 0.e0
          END DO  
-
-         ! Vertical momentum advection at u- and v-points
-         ! ----------------------------------------------
-         DO jk = 1, jpkm1
+         !
+         DO jk = 1, jpkm1         ! Vertical momentum advection at u- and v-points
             DO ji = 2, jpim1
-               ! vertical momentum advective trends
+               !                        ! vertical momentum advective trends
                zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
                zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
-               ! add the trends to the general momentum trends
+               !                        ! add the trends to the general momentum trends
                ua(ji,jj,jk) = ua(ji,jj,jk) + zua
                va(ji,jj,jk) = va(ji,jj,jk) + zva
@@ -133,19 +115,14 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-
-      ! save the vertical advection trends for diagnostic
-      ! momentum trends
-      IF( l_trddyn )   THEN
-         ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:)
-         ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:)
-
-         CALL trd_mod(ztdua, ztdva, jpdtdzad, 'DYN', kt)
-      ENDIF
-
-      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
-         CALL prt_ctl(tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask, &
-            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
-      ENDIF
-
+      !
+      IF( l_trddyn ) THEN         ! save the vertical advection trends for diagnostic
+         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt )
+      ENDIF
+      !                           ! Control print
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
+         &                       tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn' )
+      !
    END SUBROUTINE dyn_zad
 
@@ -169,22 +146,15 @@
       !!
       !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
-      !!              - Save the trends in (utrd,vtrd) ('key_trddyn')
-      !!
-      !! History :
-      !!   8.5  !  02-07  (G. Madec)  Original code
-      !!----------------------------------------------------------------------
-      !! * modules used
-      USE oce, ONLY:   zwuw => ta,   & ! use ta as 3D workspace
-                       zwvw => sa      ! use sa as 3D workspace
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt    ! ocean time-step inedx
-      
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk         ! dummy loop indices
-      REAL(wp) ::   zua, zva           ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zww                           ! temporary  workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         ztdua, ztdva                  ! temporary workspace
+      !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn')
+     !!----------------------------------------------------------------------
+      USE oce, ONLY:   zwuw => ta   ! use ta as 3D workspace
+      USE oce, ONLY:   zwvw => sa   ! use sa as 3D workspace
+      !!
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
+      !!
+      INTEGER  ::   ji, jj, jk      ! dummy loop indices
+      REAL(wp) ::   zua, zva        ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zww            ! 2D  workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrdu, ztrdv   ! 3D workspace
       !!----------------------------------------------------------------------
       
@@ -195,22 +165,17 @@
       ENDIF
 
-      ! Save ua and va trends
-      IF( l_trddyn )   THEN
-         ztdua(:,:,:) = ua(:,:,:) 
-         ztdva(:,:,:) = va(:,:,:) 
+      IF( l_trddyn )   THEN         ! Save ua and va trends
+         ztrdu(:,:,:) = ua(:,:,:) 
+         ztrdv(:,:,:) = va(:,:,:) 
       ENDIF
       
-      ! Vertical momentum advection at level w and u- and v- vertical
-      ! -------------------------------------------------------------
-      DO jk = 2, jpkm1
-         ! vertical fluxes 
-         DO jj = 2, jpj
-            DO ji = fs_2, jpi   ! vector opt.
+      DO jk = 2, jpkm1              ! Vertical momentum advection at level w and u- and v- vertical
+         DO jj = 2, jpj                   ! vertical fluxes 
+            DO ji = fs_2, jpi             ! vector opt.
                zww(ji,jj) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
             END DO
          END DO
-         ! vertical momentum advection at w-point
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
+            DO ji = fs_2, fs_jpim1        ! vector opt.
                zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1)-un(ji,jj,jk) )
                zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) )
@@ -218,8 +183,6 @@
          END DO   
       END DO
-
-      ! Surface and bottom values set to zero
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
+      DO jj = 2, jpjm1              ! Surface and bottom values set to zero
+         DO ji = fs_2, fs_jpim1           ! vector opt.
             zwuw(ji,jj, 1 ) = 0.e0
             zwvw(ji,jj, 1 ) = 0.e0
@@ -229,14 +192,11 @@
       END DO
 
-
-      ! Vertical momentum advection at u- and v-points
-      ! ----------------------------------------------
-      DO jk = 1, jpkm1
+      DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points
          DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ! vertical momentum advective trends
+            DO ji = fs_2, fs_jpim1       ! vector opt.
+               !                         ! vertical momentum advective trends
                zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
                zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
-               ! add the trends to the general momentum trends
+               !                         ! add the trends to the general momentum trends
                ua(ji,jj,jk) = ua(ji,jj,jk) + zua
                va(ji,jj,jk) = va(ji,jj,jk) + zva
@@ -245,21 +205,16 @@
       END DO
 
-      ! save the vertical advection trends for diagnostic
-      ! momentum trends
-      IF( l_trddyn )   THEN
-         ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:)
-         ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:)
-
-         CALL trd_mod(ztdua, ztdva, jpdtdzad, 'DYN', kt)
-      ENDIF
-
-      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
-         CALL prt_ctl(tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask, &
-            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
-      ENDIF
-
+      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
+         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         CALL trd_mod(ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt)
+      ENDIF
+      !                             ! Control print
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
+         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      !
    END SUBROUTINE dyn_zad
 #endif
 
-!!======================================================================
+   !!======================================================================
 END MODULE dynzad
Index: trunk/NEMO/OPA_SRC/DYN/dynzdf.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 474)
+++ trunk/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 503)
@@ -4,9 +4,11 @@
    !! Ocean dynamics :  vertical component of the momentum mixing trend
    !!==============================================================================
+   !! History :  9.0  !  05-11  (G. Madec)  Original code
+   !!----------------------------------------------------------------------
+
    !!----------------------------------------------------------------------
    !!   dyn_zdf      : Update the momentum trend with the vertical diffusion
-   !!       zdf_ctl  : ???
+   !!       zdf_ctl  : initializations of the vertical diffusion scheme
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain variables 
@@ -26,15 +28,11 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dyn_zdf   !  routine called by step.F90
+   PUBLIC   dyn_zdf    !  routine called by step.F90
 
-   !! * module variables
-   INTEGER ::                        &
-      nzdf = 0                         ! type vertical diffusion algorithm used
+   INTEGER  ::   nzdf = 0              ! type vertical diffusion algorithm used 
       !                                ! defined from ln_zdf...  namlist logicals)
 
-   !! * Module variables
-   REAL(wp) ::   &
-      r2dt                          ! time-step, = 2 rdttra except at nit000 (=rdttra) if neuler=0
+   REAL(wp) ::   r2dt                  ! time-step, = 2 rdttra
+      !                                ! except at nit000 (=rdttra) if neuler=0
 
    !! * Substitutions
@@ -44,4 +42,6 @@
    !!----------------------------------------------------------------------
    !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! $Header$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -53,16 +53,8 @@
       !!
       !! ** Purpose :   compute the vertical ocean dynamics physics.
-      !! ** Method  :
-      !! ** Action  :
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
       !!
-      !! History :
-      !!   9.0  !  05-11  (G. Madec)  Original code
-      !!---------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-
-      !! * local declarations
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         ztrdu, ztrdv                    ! 3D temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrdu, ztrdv   ! 3D workspace
       !!---------------------------------------------------------------------
 
@@ -70,10 +62,7 @@
 
       !                                          ! set time step
-      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
-         r2dt =  rdt                                    ! = rdtra (restarting with Euler time stepping)
-      ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1
-         r2dt = 2. * rdt                                ! = 2 rdttra (leapfrog)
+      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (restarting with Euler time stepping)
+      ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdttra (leapfrog)
       ENDIF
-
 
       IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
@@ -83,36 +72,30 @@
 
       SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( -1 )                                       ! esopa: test all possibility with control print
-         CALL dyn_zdf_exp    ( kt )
-         CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf0 - Ua: ', mask1=umask,               &
-            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL dyn_zdf_imp    ( kt )
-         CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf1 - Ua: ', mask1=umask,               &
-            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL dyn_zdf_imp_jki( kt )
-         CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf2 - Ua: ', mask1=umask,               &
-            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-
-      CASE ( 0 )                                       ! explicit scheme
-         CALL dyn_zdf_exp    ( kt )
-
-      CASE ( 1 )                                       ! implicit scheme (k-j-i loop)
-         CALL dyn_zdf_imp    ( kt )
-
-      CASE ( 2 )                                       ! implicit scheme (j-k-i loop)
-         CALL dyn_zdf_imp_jki( kt )
-
+      !
+      CASE ( 0 )   ;   CALL dyn_zdf_exp    ( kt, r2dt )      ! explicit scheme
+      CASE ( 1 )   ;   CALL dyn_zdf_imp    ( kt, r2dt )      ! implicit scheme (k-j-i loop)
+      CASE ( 2 )   ;   CALL dyn_zdf_imp_jki( kt, r2dt )      ! implicit scheme (j-k-i loop)
+      !
+      CASE ( -1 )                                      ! esopa: test all possibility with control print
+         ;             CALL dyn_zdf_exp    ( kt, r2dt )
+         ;             CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf0 - Ua: ', mask1=umask,               &
+            &                        tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         ;             CALL dyn_zdf_imp    ( kt, r2dt )
+         ;             CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf1 - Ua: ', mask1=umask,               &
+            &                        tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         ;             CALL dyn_zdf_imp_jki( kt, r2dt )
+         ;             CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf2 - Ua: ', mask1=umask,               &
+            &                        tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       END SELECT
 
-      IF( l_trddyn )   THEN                         ! save the vertical diffusive trends for further diagnostics
+      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-         CALL trd_mod( ztrdu, ztrdv, jpttdzdf, 'DYN', kt )
+         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_zdf, 'DYN', kt )
       ENDIF
-
       !                                          ! print mean trends (used for debugging)
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
             &                    tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-
+      !
    END SUBROUTINE dyn_zdf
 
@@ -122,46 +105,32 @@
       !!                 ***  ROUTINE zdf_ctl  ***
       !!
-      !! ** Purpose :   initializations of the vertical ocean dynamics physics
+      !! ** Purpose :   initializations of the vertical diffusion scheme
       !!
       !! ** Method  :   implicit (euler backward) scheme (default)
       !!                explicit (time-splitting) scheme if ln_zdfexp=T
       !!                OpenMP / NEC autotasking: use j-k-i loops
-      !!
-      !! History :
-      !!   9.0  !  05-11  (G. Madec)  Original code
       !!----------------------------------------------------------------------
-      !! * Module used
       USE zdftke
       USE zdfkpp
       !!----------------------------------------------------------------------
 
-      !  Define the vertical dynamics physics scheme
-      ! ==========================================
-
       ! Choice from ln_zdfexp read in namelist in zdfini
-      IF( ln_zdfexp ) THEN               ! use explicit scheme
-         nzdf = 0
-      ELSE                               ! use implicit scheme
-         nzdf = 1
+      IF( ln_zdfexp ) THEN   ;   nzdf = 0           ! use explicit scheme
+      ELSE                   ;   nzdf = 1           ! use implicit scheme
       ENDIF
 
       ! Force implicit schemes
-      IF( lk_zdftke .OR. lk_zdfkpp   )   nzdf = 1      ! TKE or KPP physics
-      IF( ln_dynldf_iso              )   nzdf = 1      ! iso-neutral lateral physics
-      IF( ln_dynldf_hor .AND. ln_sco )   nzdf = 1      ! horizontal lateral physics in s-coordinate
+      IF( lk_zdftke .OR. lk_zdfkpp   )   nzdf = 1   ! TKE or KPP physics
+      IF( ln_dynldf_iso              )   nzdf = 1   ! iso-neutral lateral physics
+      IF( ln_dynldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate
 
       ! OpenMP / NEC autotasking
 #if defined key_mpp_omp
-      IF( nzdf == 1 )   nzdf = 2                       ! j-k-i loop
+      IF( nzdf == 1 )   nzdf = 2                    ! j-k-i loop
 #endif
 
-!!bug
-!     IF( ln_dynldf_iso              )   nzdf = 3      ! iso-neutral lateral physics
-!     IF( ln_dynldf_hor .AND. ln_sco )   nzdf = 3      ! horizontal lateral physics in s-coordinate
-!!bug
-      ! Test: esopa
-      IF( lk_esopa )    nzdf = -1                      ! All schemes used
+      IF( lk_esopa )    nzdf = -1                   ! Esopa key: All schemes used
 
-      IF(lwp) THEN
+      IF(lwp) THEN                                  ! Print the choice
          WRITE(numout,*)
          WRITE(numout,*) 'dyn:zdf_ctl : vertical dynamics physics scheme'
@@ -172,5 +141,5 @@
          IF( nzdf ==  2 )   WRITE(numout,*) '              Implicit (euler backward) scheme with j-k-i loops'
       ENDIF
-
+      !
    END SUBROUTINE zdf_ctl
 
Index: trunk/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzdf_exp.F90	(revision 474)
+++ trunk/NEMO/OPA_SRC/DYN/dynzdf_exp.F90	(revision 503)
@@ -4,4 +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
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -16,5 +20,4 @@
    USE in_out_manager  ! I/O manager
    USE taumod          ! surface ocean stress
-   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -30,10 +33,10 @@
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! 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_exp( kt )
+   SUBROUTINE dyn_zdf_exp( kt, p2dt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_zdf_exp  ***
@@ -50,25 +53,13 @@
       !!
       !! ** Action : - Update (ua,va) with the vertical diffusive trend
-      !!             - Save the trends in (ztdua,ztdva) ('key_trddyn')
-      !!
-      !! 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
-      !!   9.0  !  04-08  (C. Talandier)  New trends organization
       !!---------------------------------------------------------------------
-      !! * Modules used     
-      USE oce, ONLY :    ztdua => ta,    & ! use ta as 3D workspace   
-                         ztdva => sa       ! use sa as 3D workspace   
       !! * Arguments
-      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
+      INTEGER , INTENT( in ) ::   kt                           ! ocean time-step index
+      REAL(wp), INTENT( in ) ::   p2dt                         ! time-step 
 
       !! * Local declarations
-      INTEGER ::   &
-         ji, jj, jk, jl                    ! dummy loop indices
-      REAL(wp) ::   &
-         zrau0r, zlavmr, z2dt, zua, zva    ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpk) ::    &
-         zwx, zwy, zwz, zww                ! temporary workspace arrays
+      INTEGER ::   ji, jj, jk, jl                              ! dummy loop indices
+      REAL(wp) ::   zrau0r, zlavmr, zua, zva                   ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpk) ::   zwx, zwy, zwz, zww     ! temporary workspace arrays
       !!----------------------------------------------------------------------
 
@@ -83,7 +74,4 @@
       zrau0r = 1. / rau0                                   ! inverse of the reference density
       zlavmr = 1. / float( n_zdfexp )                      ! inverse of the number of sub time step
-      z2dt = 2. * rdt                                      ! Leap-frog environnement
-
-      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt    ! Euler time stepping when starting from rest
 
       !                                                ! ===============
@@ -124,6 +112,6 @@
                   va(ji,jj,jk) = va(ji,jj,jk) + zva
 
-                  zwx(ji,jk) = zwx(ji,jk) + z2dt*zua*umask(ji,jj,jk)
-                  zwz(ji,jk) = zwz(ji,jk) + z2dt*zva*vmask(ji,jj,jk)
+                  zwx(ji,jk) = zwx(ji,jk) + p2dt*zua*umask(ji,jj,jk)
+                  zwz(ji,jk) = zwz(ji,jk) + p2dt*zva*vmask(ji,jj,jk)
                END DO  
             END DO  
Index: trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 474)
+++ trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 503)
@@ -4,4 +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
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -11,5 +15,5 @@
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! 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)
    !!----------------------------------------------------------------------
    !! * Modules used
@@ -20,5 +24,4 @@
    USE in_out_manager  ! I/O manager
    USE taumod          ! surface ocean stress
-   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -40,5 +43,5 @@
 
 
-   SUBROUTINE dyn_zdf_imp( kt )
+   SUBROUTINE dyn_zdf_imp( kt, p2dt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_zdf_imp  ***
@@ -58,25 +61,17 @@
       !! ** Action : - Update (ua,va) arrays with the after vertical diffusive
       !!               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
-      !!   9.0  !  04-08  (C. Talandier)  New trends organization
       !!---------------------------------------------------------------------
       !! * Modules used
-      USE oce, ONLY :  zwd   => ta,   &  ! use ta as workspace
-                       zws   => sa       ! use sa as workspace
+      USE oce, ONLY :  zwd   => ta,   &                ! use ta as workspace
+                       zws   => sa                     ! use sa as workspace
 
       !! * Arguments
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      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, z2dt,                &  ! temporary scalars
-         z2dtf, zcoef, zzws, zrhs        !    "         "
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         zwi                             ! temporary workspace arrays
+      INTEGER ::   ji, jj, jk                          ! dummy loop indices
+      REAL(wp) ::   zrau0r, z2dtf, zcoef, zzws, zrhs   ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! temporary workspace arrays
       !!----------------------------------------------------------------------
 
@@ -90,8 +85,4 @@
       ! --------------------------------
       zrau0r = 1. / rau0      ! inverse of the reference density
-      z2dt   = 2. * rdt       ! Leap-frog environnement
-
-      ! Euler time stepping when starting from rest
-      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
 
       ! 1. Vertical diffusion on u
@@ -104,5 +95,5 @@
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcoef = - z2dt / fse3u(ji,jj,jk)
+               zcoef = - p2dt / fse3u(ji,jj,jk)
                zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
                zzws          = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
@@ -150,14 +141,14 @@
 !!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
 !!!         ua(ji,jj,1) = ub(ji,jj,1)  &
-!!!                      + z2dt * ( ua(ji,jj,1) + taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) )
-            z2dtf = z2dt / ( fse3u(ji,jj,1)*rau0 )
+!!!                      + p2dt * ( ua(ji,jj,1) + taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) )
+            z2dtf = p2dt / ( fse3u(ji,jj,1)*rau0 )
             ua(ji,jj,1) = ub(ji,jj,1)  &
-                         + z2dt *  ua(ji,jj,1) + z2dtf * taux(ji,jj)
-         END DO
-      END DO
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zrhs = ub(ji,jj,jk) + z2dt * ua(ji,jj,jk)   ! zrhs=right hand side
+                         + p2dt *  ua(ji,jj,1) + z2dtf * taux(ji,jj)
+         END DO
+      END DO
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1   
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
                ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
             END DO
@@ -183,5 +174,5 @@
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) / z2dt
+               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) / p2dt
             END DO
          END DO
@@ -198,5 +189,5 @@
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcoef = -z2dt / fse3v(ji,jj,jk)
+               zcoef = -p2dt / fse3v(ji,jj,jk)
                zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
                zzws       = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
@@ -245,8 +236,8 @@
 !!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
 !!!         va(ji,jj,1) = vb(ji,jj,1)  &
-!!!                      + z2dt * ( va(ji,jj,1) + tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) )
-            z2dtf = z2dt / ( fse3v(ji,jj,1)*rau0 )
+!!!                      + p2dt * ( va(ji,jj,1) + tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) )
+            z2dtf = p2dt / ( fse3v(ji,jj,1)*rau0 )
             va(ji,jj,1) = vb(ji,jj,1)  &
-                         + z2dt * va(ji,jj,1) + z2dtf * tauy(ji,jj)
+                         + p2dt * va(ji,jj,1) + z2dtf * tauy(ji,jj)
          END DO
       END DO
@@ -254,5 +245,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zrhs = vb(ji,jj,jk) + z2dt * va(ji,jj,jk)   ! zrhs=right hand side
+               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
                va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
             END DO
@@ -274,31 +265,9 @@
       END DO
 
-! flux de surface doit etre calcule dans trdmod et boootom stress
-! deduit par integration verticale dans trmod pour jpdtdzdf
-!RB      IF( l_trddyn )  THEN 
-!         ! diagnose surface and bottom momentum fluxes
-!         DO jj = 2, jpjm1   
-!            DO ji = fs_2, fs_jpim1   ! vector opt.
-!               ! save the surface forcing momentum fluxes
-!               ztsy(ji,jj) = tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 )
-!               ztsx(ji,jj) = taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 )
-!               ! save bottom friction momentum fluxes
-!               ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
-!               ikbvm1 = MAX( ikbv-1, 1 )
-!               ztby(ji,jj) = - avmv(ji,jj,ikbv) * va(ji,jj,ikbvm1)   &
-!                  / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) )
-!               ! subtract surface forcing and bottom friction trend from vertical
-!               ! diffusive momentum trend
-!               ztdva(ji,jj,1     ) = ztdva(ji,jj,1     ) - ztsy(ji,jj)
-!               ztdva(ji,jj,ikbvm1) = ztdva(ji,jj,ikbvm1) - ztby(ji,jj)
-!            END DO
-!         END DO
-!      ENDIF
-
       ! Normalization to obtain the general momentum trend va
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) / z2dt
+               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) / p2dt
             END DO
          END DO
Index: trunk/NEMO/OPA_SRC/DYN/dynzdf_imp_jki.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzdf_imp_jki.F90	(revision 474)
+++ trunk/NEMO/OPA_SRC/DYN/dynzdf_imp_jki.F90	(revision 503)
@@ -4,4 +4,6 @@
    !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
    !!==============================================================================
+   !! History : 8.5  !  02-08  (G. Madec)  auto-tasking option
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -17,5 +19,4 @@
    USE in_out_manager  ! I/O manager
    USE taumod          ! surface ocean stress
-   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -30,9 +31,11 @@
    !!----------------------------------------------------------------------
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! $Header$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
 
-   SUBROUTINE dyn_zdf_imp_jki( kt )
+   SUBROUTINE dyn_zdf_imp_jki( kt, p2dt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_zdf_imp_jki  ***
@@ -52,24 +55,15 @@
       !! ** Action : - Update (ua,va) arrays with the after vertical diffusive
       !!               mixing trend.
-      !!
-      !! History :
-      !!   8.5  !  02-08  (G. Madec)  auto-tasking option
-      !!   9.0  !  04-08  (C. Talandier)  New trends organization
       !!---------------------------------------------------------------------
-      !! * Modules used     
-      USE oce, ONLY :    ztdua => ta,  & ! use ta as 3D workspace   
-                         ztdva => sa     ! use sa as 3D workspace   
-
       !! * Arguments
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      INTEGER , INTENT( in ) ::   kt                          ! ocean time-step index
+      REAL(wp), INTENT( in ) ::   p2dt                        ! ocean time-step index
 
       !! * Local declarations
-      INTEGER ::   ji, jj, jk            ! dummy loop indices
-      INTEGER ::   ikst, ikenm2, ikstp1  ! temporary integers
-      REAL(wp) ::  zrau0r, z2dt,       & !temporary scalars
-         &         z2dtf, zcoef, zzws
-      REAL(wp), DIMENSION(jpi,jpk) ::  &
-         zwx, zwy, zwz,                & ! workspace
-         zwd, zws, zwi, zwt
+      INTEGER ::   ji, jj, jk                                 ! dummy loop indices
+      INTEGER ::   ikst, ikenm2, ikstp1                       ! temporary integers
+      REAL(wp) ::  zrau0r, z2dtf, zcoef, zzws                 ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpk) ::  zwx, zwy, zwz,       & ! workspace
+         &                             zwd, zws, zwi, zwt
       !!----------------------------------------------------------------------
 
@@ -84,8 +78,4 @@
       ! --------------------------------
       zrau0r = 1. / rau0      ! inverse of the reference density
-      z2dt   = 2. * rdt       ! Leap-frog environnement
-
-      ! Euler time stepping when starting from rest
-      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
 
       !                                                ! ===============
@@ -101,10 +91,10 @@
          DO jk = 1, jpkm1
             DO ji = 2, jpim1
-               zcoef = - z2dt / fse3u(ji,jj,jk)
+               zcoef = - p2dt / fse3u(ji,jj,jk)
                zwi(ji,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
                zzws       = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
                zws(ji,jk) = zzws * umask(ji,jj,jk+1)
                zwd(ji,jk) = 1. - zwi(ji,jk) - zzws
-               zwy(ji,jk) = ub(ji,jj,jk) + z2dt * ua(ji,jj,jk)
+               zwy(ji,jk) = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
             END DO
          END DO
@@ -112,5 +102,5 @@
          ! Surface boudary conditions
          DO ji = 2, jpim1
-            z2dtf = z2dt / ( fse3u(ji,jj,1)*rau0 )
+            z2dtf = p2dt / ( fse3u(ji,jj,1)*rau0 )
             zwi(ji,1) = 0.
             zwd(ji,1) = 1. - zws(ji,1)
@@ -175,5 +165,5 @@
          DO jk = 1, jpkm1
             DO ji = 2, jpim1
-               ua(ji,jj,jk) = ( zwx(ji,jk) - ub(ji,jj,jk) ) / z2dt
+               ua(ji,jj,jk) = ( zwx(ji,jk) - ub(ji,jj,jk) ) / p2dt
             END DO
          END DO
@@ -188,10 +178,10 @@
          DO jk = 1, jpkm1
             DO ji = 2, jpim1
-               zcoef = -z2dt/fse3v(ji,jj,jk)
+               zcoef = -p2dt/fse3v(ji,jj,jk)
                zwi(ji,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
                zzws       = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
                zws(ji,jk) =  zzws * vmask(ji,jj,jk+1)
                zwd(ji,jk) = 1. - zwi(ji,jk) - zzws
-               zwy(ji,jk) = vb(ji,jj,jk) + z2dt * va(ji,jj,jk)
+               zwy(ji,jk) = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
             END DO
          END DO
@@ -199,5 +189,5 @@
          ! Surface boudary conditions
          DO ji = 2, jpim1
-            z2dtf = z2dt / ( fse3v(ji,jj,1)*rau0 )
+            z2dtf = p2dt / ( fse3v(ji,jj,1)*rau0 )
             zwi(ji,1) = 0.e0
             zwd(ji,1) = 1. - zws(ji,1)
@@ -262,5 +252,5 @@
          DO jk = 1, jpkm1
             DO ji = 2, jpim1
-               va(ji,jj,jk) = ( zwx(ji,jk) - vb(ji,jj,jk) ) / z2dt
+               va(ji,jj,jk) = ( zwx(ji,jk) - vb(ji,jj,jk) ) / p2dt
             END DO
          END DO
