Index: trunk/NEMO/OPA_SRC/DYN/dynhpg_jki.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynhpg_jki.F90	(revision 456)
+++ trunk/NEMO/OPA_SRC/DYN/dynhpg_jki.F90	(revision 456)
@@ -0,0 +1,417 @@
+MODULE dynhpg_jki
+   !!======================================================================
+   !!                       ***  MODULE  dynhpg_jki  ***
+   !! Ocean dynamics:  hydrostatic pressure gradient trend
+   !!======================================================================
+
+   !!----------------------------------------------------------------------
+   !!   hpg...o_jki : update the momentum trend with the horizontal
+   !!                  gradient of the hydrostatic pressure
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   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
+
+   IMPLICIT NONE
+   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
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE hpg_sco_jki( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE hpg_sco_jki  ***
+      !!        
+      !! ** Purpose :   Compute the now momentum trend due to the horizontal 
+      !!     gradient of the hydrostatic pressure. Add it to the general
+      !!     momentum trend.
+      !!
+      !! ** Method  :   The now hydrostatic pressure gradient at a given level
+      !!      jk is computed by taking the vertical integral of the in-situ 
+      !!      density gradient along the model level from the suface to that
+      !!      level. s-coordinate case (ln_sco=T): a corrective term is
+      !!      added to the horizontal pressure gradient :
+      !!         zhpi = grav .....   + 1/e1u mi(rhd) di[ grav dep3w ]
+      !!         zhpj = grav .....   + 1/e2v mj(rhd) dj[ grav dep3w ]
+      !!      add it to the general momentum trend (ua,va).
+      !!         ua = ua - 1/e1u * zhpi
+      !!         va = va - 1/e2v * zhpj
+      !!      j-k-i loop (j-slab) ('key_mpp_omp')
+      !!
+      !! ** 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
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'hpg_sco_jki : s-coordinate hydrostatic pressure gradient trend'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   OpenMP / NEC autotasking case (j-k-i loop)'
+      ENDIF
+
+      ! Save ua and va trends
+      IF( l_trddyn )   THEN
+         ztdua(:,:,:) = ua(:,:,:) 
+         ztdva(:,:,:) = va(:,:,:) 
+      ENDIF
+
+      ! 0. Local constant initialization
+      ! --------------------------------
+      zcoef0 = - grav * 0.5
+      zuap   = 0.e0
+      zvap   = 0.e0
+
+      !                                                ! ===============
+      DO jj = 2, jpjm1                                 !  Vertical slab
+         !                                             ! ===============
+         ! 1. Surface value
+         ! ----------------
+         DO ji = 2, jpim1
+            ! hydrostatic pressure gradient along s-surfaces
+            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj)   &
+                       * ( fse3w(ji+1,jj,1) * rhd(ji+1,jj,1) - fse3w(ji,jj,1) * rhd(ji,jj,1)  )
+            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj)   &
+                       * ( fse3w(ji,jj+1,1) * rhd(ji,jj+1,1) - fse3w(ji,jj,1) * rhd(ji,jj,1)  )
+            ! s-coordinate pressure gradient correction
+            zuap = -zcoef0 * ( rhd(ji+1,jj,1) + rhd(ji,jj,1) )   &
+                 * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
+            zvap = -zcoef0 * ( rhd(ji,jj+1,1) + rhd(ji,jj,1) )   &
+                 * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
+            ! add to the general momentum trend
+            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
+            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 ji = 2, jpim1
+               ! hydrostatic pressure gradient along s-surfaces
+               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
+                  &           * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   &
+                  &              -fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  )
+               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
+                  &           * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   &
+                  &              -fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  )
+               ! s-coordinate pressure gradient correction 
+               zuap = -zcoef0 * ( rhd(ji+1,jj  ,jk) + rhd(ji,jj,jk) )   &
+                    * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj)
+               zvap = -zcoef0 * ( rhd(ji  ,jj+1,jk) + rhd(ji,jj,jk) )   &
+                    * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
+               ! add to the general momentum trend
+               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
+               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
+            END DO
+         END DO
+         !                                             ! ===============
+      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
+
+
+   SUBROUTINE hpg_zps_jki( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE hpg_zps_jki  ***
+      !!  
+      !! ** Purpose :   Compute the now momentum trend due to the hor. gradient
+      !!      of the hydrostatic pressure. Add it to the general momentum trend.
+      !!
+      !! ** Method  :   The now hydrostatic pressure gradient at a given level
+      !!      jk is computed by taking the vertical integral of the in-situ
+      !!      density gradient along the model level from the suface to that 
+      !!      level:    zhpi = grav .....
+      !!                zhpj = grav .....
+      !!      add it to the general momentum trend (ua,va).
+      !!            ua = ua - 1/e1u * zhpi
+      !!            va = va - 1/e2v * zhpj
+      !!      j-k-i loop (j-slab) ('key_mpp_omp')
+      !!
+      !! ** 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
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'hpg_zps_jki : z-coord. partial steps hydrostatic pressure gradient trend'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   OpenMP / NEC autotasking case (j-k-i loop)'
+      ENDIF
+
+      ! Save ua and va trends
+      IF( l_trddyn )   THEN
+         ztdua(:,:,:) = ua(:,:,:) 
+         ztdva(:,:,:) = va(:,:,:) 
+      ENDIF
+
+      ! 0. Local constant initialization
+      ! --------------------------------
+      zcoef0 = - grav * 0.5
+      zuap   = 0.e0
+      zvap   = 0.e0
+      !                                                ! ===============
+      DO jj = 2, jpjm1                                 !  Vertical slab
+         !                                             ! ===============
+         ! 1. Surface value
+         ! ----------------
+         DO ji = 2, jpim1
+            zcoef1 = zcoef0 * fse3w(ji,jj,1)
+            ! hydrostatic pressure gradient
+            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
+            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
+            ! add to the general momentum trend
+            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
+            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 ji = 2, jpim1
+               zcoef1 = zcoef0 * fse3w(ji,jj,jk)
+               ! hydrostatic pressure gradient
+               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
+                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   &
+                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj)
+
+               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
+                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   &
+                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj)
+               ! add to the general momentum trend
+               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
+            END DO 
+         END DO
+
+         ! partial steps correction at the last level  (new gradient with  intgrd.F)
+         DO ji = 2, jpim1
+            iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1
+            ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1
+            zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj  ,iku) )
+            zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji  ,jj+1,ikv) )
+            ! on i-direction
+            IF ( iku > 2 ) THEN
+               ! subtract old value  
+               ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)
+               ! compute the new one   
+               zhpi (ji,jj,iku) = zhpi(ji,jj,iku-1)   &
+                  + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj)
+               ! add the new one to the general momentum trend
+               ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)
+            ENDIF
+            ! on j-direction
+            IF ( ikv > 2 ) THEN
+               ! subtract old value  
+               va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)
+               ! compute the new one   
+               zhpj (ji,jj,ikv) = zhpj(ji,jj,ikv-1)   &
+                  + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj)
+               ! add the new one to the general momentum trend
+               va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)
+            ENDIF
+         END DO
+         !                                             ! ===============
+      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
+
+
+   SUBROUTINE hpg_zco_jki( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE hpg_zco_jki  ***
+      !!   
+      !! ** Purpose :   Compute the now momentum trend due to the horizontal
+      !!      gradient of the hydrostatic pressure. Add it to the general
+      !!      momentum trend.
+      !!
+      !! ** Method  :   The now hydrostatic pressure gradient at a given level
+      !!      jk is computed by taking the vertical integral of the in-situ  
+      !!      density gradient along the model level from the suface to that
+      !!      level:    zhpi = grav .....
+      !!                zhpj = grav .....
+      !!      add it to the general momentum trend (ua,va).
+      !!            ua = ua - 1/e1u * zhpi
+      !!            va = va - 1/e2v * zhpj
+      !!      j-k-i loop (j-slab) ('key_mpp_omp')
+      !!
+      !! ** 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
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'hpg_zco_jki : z-coordinate hydrostatic pressure gradient trend'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   OpenMP / NEC auto-tasking (j-k-i loop)'
+      ENDIF
+
+      ! Save ua and va trends
+      IF( l_trddyn )   THEN
+         ztdua(:,:,:) = ua(:,:,:) 
+         ztdva(:,:,:) = va(:,:,:) 
+      ENDIF
+
+      ! 0. Local constant initialization
+      ! --------------------------------
+      zcoef0 = - grav * 0.5
+      zuap   = 0.e0
+      zvap   = 0.e0
+
+      !                                                ! ===============
+      DO jj = 2, jpjm1                                 !  Vertical slab
+         !                                             ! ===============
+         ! 1. Surface value
+         ! ----------------
+         
+
+         DO ji = 2, jpim1
+            zcoef1 = zcoef0 * fse3w(ji,jj,1)
+            ! hydrostatic pressure gradient
+            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
+            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
+            ! add to the general momentum trend
+            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
+            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 ji = 2, jpim1
+               zcoef1 = zcoef0 * fse3w(ji,jj,jk)
+               ! hydrostatic pressure gradient
+               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
+                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   &
+                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj)
+
+               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
+                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   &
+                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj)
+               ! add to the general momentum trend
+               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
+            END DO 
+         END DO
+         !                                             ! ===============
+      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
+
+   !!======================================================================
+END MODULE dynhpg_jki
Index: trunk/NEMO/OPA_SRC/DYN/dynldf.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 456)
+++ trunk/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 456)
@@ -0,0 +1,242 @@
+MODULE dynldf
+   !!======================================================================
+   !!                       ***  MODULE  dynldf  ***
+   !! Ocean physics:  lateral diffusivity trends 
+   !!=====================================================================
+
+   !!----------------------------------------------------------------------
+   !!   dyn_ldf     : update the dynamics trend with the lateral diffusion
+   !!   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)
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! *  Routine accessibility
+   PUBLIC dyn_ldf      ! called by step.F90 
+
+   !! * module variables
+   INTEGER ::                        & 
+      nldf = 0                         ! type of lateral diffusion used
+      !                                ! defined from ln_dynldf_... namlist logicals)
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!---------------------------------------------------------------------------------
+   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!---------------------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE dyn_ldf( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_ldf  ***
+      !! 
+      !! ** Purpose :   compute the lateral ocean dynamics physics.
+      !!
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
+
+      !! * local declarations
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
+         ztrdu, ztrdv                         ! 3D temporary workspace
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 )   CALL dyn_ldf_ctl      ! initialisation & control of options
+
+      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
+         ztrdu(:,:,:) = ua(:,:,:) 
+         ztrdv(:,:,:) = va(:,:,:) 
+      ENDIF
+
+      SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
+      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 )
+      END SELECT
+
+      !                                          ! save the horizontal diffusive trends for further diagnostics
+      IF( l_trddyn )   THEN
+         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         CALL trd_mod( ztrdu, ztrdv, jpdtdldf, '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
+
+   END SUBROUTINE dyn_ldf
+
+
+   SUBROUTINE dyn_ldf_ctl
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_ldf_ctl  ***
+      !! 
+      !! ** 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 
+      !!----------------------------------------------------------------------
+
+      !  Define the lateral dynamics physics parameters
+      ! =============================================
+    
+      ! Namelist nam_dynldf already read in ldfdyn module
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         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
+      ENDIF
+
+      ! control the consistency
+      ioptio = 0
+      IF( ln_dynldf_lap   )   ioptio = ioptio + 1
+      IF( ln_dynldf_bilap )   ioptio = ioptio + 1
+      IF( ioptio /= 1 )   THEN
+          IF(lwp) WRITE(numout,cform_err)
+          IF(lwp) WRITE(numout,*) '          use ONE of the 2 lap/bilap operator type on dynamics'
+          nstop = nstop + 1
+      ENDIF
+      ioptio = 0
+      IF( ln_dynldf_level )   ioptio = ioptio + 1
+      IF( ln_dynldf_hor   )   ioptio = ioptio + 1
+      IF( ln_dynldf_iso   )   ioptio = ioptio + 1
+      IF( ioptio /= 1 ) THEN
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) '          use only ONE direction (level/hor/iso)'
+         nstop = nstop + 1
+      ENDIF
+
+      ! defined the type of lateral diffusion from ln_dynldf_... logicals
+      ierr = 0
+      IF ( ln_dynldf_lap ) THEN      ! laplacian operator
+         IF ( ln_zco ) THEN                ! z-coordinate
+            IF ( ln_dynldf_level )   nldf = 0      ! iso-level  (no rotation)
+            IF ( ln_dynldf_hor   )   nldf = 0      ! horizontal (no rotation)
+            IF ( ln_dynldf_iso   )   nldf = 1      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_zps ) THEN             ! z-coordinate
+            IF ( ln_dynldf_level )   ierr = 1      ! iso-level not allowed
+            IF ( ln_dynldf_hor   )   nldf = 0      ! horizontal (no rotation)
+            IF ( ln_dynldf_iso   )   nldf = 1      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_sco ) THEN             ! z-coordinate
+            IF ( ln_dynldf_level )   nldf = 0      ! iso-level  (no rotation)
+            IF ( ln_dynldf_hor   )   nldf = 1      ! horizontal (   rotation)
+            IF ( ln_dynldf_iso   )   nldf = 1      ! isoneutral (   rotation)
+         ENDIF
+      ENDIF
+
+      IF( ln_dynldf_bilap ) THEN      ! bilaplacian operator
+         IF ( ln_zco ) THEN                ! z-coordinate
+            IF ( ln_dynldf_level )   nldf = 2      ! iso-level  (no rotation)
+            IF ( ln_dynldf_hor   )   nldf = 2      ! horizontal (no rotation)
+            IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_zps ) THEN             ! z-coordinate
+            IF ( ln_dynldf_level )   ierr = 1      ! iso-level not allowed 
+            IF ( ln_dynldf_hor   )   nldf = 2      ! horizontal (no rotation)
+            IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_sco ) THEN             ! z-coordinate
+            IF ( ln_dynldf_level )   nldf = 2      ! iso-level  (no rotation)
+            IF ( ln_dynldf_hor   )   nldf = 3      ! horizontal (   rotation)
+            IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation)
+         ENDIF
+      ENDIF
+
+      IF( ierr == 1 ) THEN
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) ' iso-level in z-coordinate - partial step, not allowed'
+         nstop = nstop + 1
+      ENDIF
+      IF( ierr == 2 ) THEN
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) ' isoneutral bilaplacian operator does not exist'
+         nstop = nstop + 1
+      ENDIF
+      IF( nldf == 1 .OR. nldf == 3 ) THEN      ! rotation
+         IF( .NOT.lk_ldfslp ) THEN
+            IF(lwp) WRITE(numout,cform_err)
+            IF(lwp) WRITE(numout,*) '          the rotation of the diffusive tensor require key_ldfslp'
+            nstop = nstop + 1
+         ENDIF
+      ENDIF
+
+      IF( lk_esopa ) THEN
+         IF(lwp ) WRITE(numout,*) '          esopa test: use all lateral physics options'
+         nldf = -1
+      ENDIF
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         IF( nldf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used'
+         IF( nldf ==  0 )   WRITE(numout,*) '              laplacian operator'
+         IF( nldf ==  1 )   WRITE(numout,*) '              Rotated laplacian operator'
+         IF( nldf ==  2 )   WRITE(numout,*) '              bilaplacian operator'
+         IF( nldf ==  3 )   WRITE(numout,*) '              Rotated bilaplacian'
+      ENDIF
+
+   END SUBROUTINE dyn_ldf_ctl
+
+   !!======================================================================
+END MODULE dynldf
Index: trunk/NEMO/OPA_SRC/DYN/dynzdf.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 456)
+++ trunk/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 456)
@@ -0,0 +1,178 @@
+MODULE dynzdf
+   !!==============================================================================
+   !!                 ***  MODULE  dynzdf  ***
+   !! Ocean dynamics :  vertical component of the momentum mixing trend
+   !!==============================================================================
+   !!----------------------------------------------------------------------
+   !!   dyn_zdf      : Update the momentum trend with the vertical diffusion
+   !!       zdf_ctl  : ???
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables 
+   USE zdf_oce         ! ocean vertical physics variables
+
+   USE dynzdf_exp      ! vertical diffusion: explicit (dyn_zdf_exp     routine)
+   USE dynzdf_imp      ! vertical diffusion: implicit (dyn_zdf_imp     routine)
+   USE dynzdf_imp_jki  ! vertical diffusion  implicit (dyn_zdf_imp_jki routine)
+
+   USE ldfdyn_oce      ! ocean dynamics: lateral physics
+   USE trdmod          ! ocean active dynamics and tracers trends 
+   USE trdmod_oce      ! ocean variables trends
+   USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC dyn_zdf   !  routine called by step.F90
+
+   !! * module variables
+   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
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "zdfddm_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!----------------------------------------------------------------------
+
+CONTAINS
+   
+   SUBROUTINE dyn_zdf( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_zdf  ***
+      !!
+      !! ** Purpose :   compute the vertical ocean dynamics physics.
+      !! ** Method  :
+      !! ** Action  :
+      !!
+      !! 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
+      !!---------------------------------------------------------------------
+
+      IF( kt == nit000 )   CALL zdf_ctl          ! initialisation & control of options
+
+      !                                          ! 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)
+      ENDIF
+
+
+      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
+         ztrdu(:,:,:) = ua(:,:,:)
+         ztrdv(:,:,:) = va(:,:,:)
+      ENDIF
+
+      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 )
+
+      END SELECT
+
+      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 )
+      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
+
+
+   SUBROUTINE zdf_ctl
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE zdf_ctl  ***
+      !!
+      !! ** Purpose :   initializations of the vertical ocean dynamics physics
+      !!
+      !! ** 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
+      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
+
+      ! OpenMP / NEC autotasking
+#if defined key_mpp_omp
+      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(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dyn:zdf_ctl : vertical dynamics physics scheme'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         IF( nzdf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used'
+         IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme'
+         IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme'
+         IF( nzdf ==  2 )   WRITE(numout,*) '              Implicit (euler backward) scheme with j-k-i loops'
+      ENDIF
+
+   END SUBROUTINE zdf_ctl
+
+   !!==============================================================================
+END MODULE dynzdf
Index: trunk/NEMO/OPA_SRC/DYN/dynzdf_imp_jki.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzdf_imp_jki.F90	(revision 456)
+++ trunk/NEMO/OPA_SRC/DYN/dynzdf_imp_jki.F90	(revision 456)
@@ -0,0 +1,274 @@
+MODULE dynzdf_imp_jki
+   !!==============================================================================
+   !!                    ***  MODULE  dynzdf_imp_jki  ***
+   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
+   !!==============================================================================
+
+   !!----------------------------------------------------------------------
+   !!   dyn_zdf_imp_jki : update the momentum trend with the vertical
+   !!                     diffusion using an implicit time-stepping and
+   !!                     j-k-i loops.
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE zdf_oce         ! ocean vertical physics
+   USE in_out_manager  ! I/O manager
+   USE taumod          ! surface ocean stress
+   USE prtctl          ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC dyn_zdf_imp_jki     ! called by step.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE dyn_zdf_imp_jki( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_zdf_imp_jki  ***
+      !!                   
+      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
+      !!      and the surface forcing, and add it to the general trend of 
+      !!      the momentum equations.
+      !!
+      !! ** Method  :   The vertical momentum mixing trend is given by :
+      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
+      !!      backward time stepping
+      !!      Surface boundary conditions: wind stress input
+      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
+      !!      Add this trend to the general trend ua :
+      !!         ua = ua + dz( avmu dz(u) )
+      !!
+      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive
+      !!               mixing trend.
+      !!
+      !! 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
+
+      !! * 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
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp_jki : vertical momentum diffusion implicit operator'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~   auto-task case (j-k-i loop)'
+      ENDIF
+
+
+      ! 0. Local constant initialization
+      ! --------------------------------
+      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
+
+      !                                                ! ===============
+      DO jj = 2, jpjm1                                 !  Vertical slab
+         !                                             ! ===============
+         ! 1. Vertical diffusion on u
+         ! ---------------------------
+
+         ! Matrix and second member construction
+         ! bottom boundary condition: only zws must be masked as avmu can take
+         ! non zero value at the ocean bottom depending on the bottom friction
+         ! used (see zdfmix.F)
+         DO jk = 1, jpkm1
+            DO ji = 2, jpim1
+               zcoef = - z2dt / 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)
+            END DO
+         END DO
+
+         ! Surface boudary conditions
+         DO ji = 2, jpim1
+            z2dtf = z2dt / ( fse3u(ji,jj,1)*rau0 )
+            zwi(ji,1) = 0.
+            zwd(ji,1) = 1. - zws(ji,1)
+            zwy(ji,1) = zwy(ji,1) + z2dtf * taux(ji,jj)
+         END DO
+
+         ! Matrix inversion starting from the first level
+         ikst = 1
+!!----------------------------------------------------------------------
+!!         ZDF.MATRIXSOLVER
+!!       ********************
+!!----------------------------------------------------------------------
+!! Matrix inversion
+!   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
+!
+!        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
+!        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
+!        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
+!        (        ...               )( ...  ) ( ...  )
+!        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
+!
+!   m is decomposed in the product of an upper and lower triangular
+!   matrix
+!   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
+!   The second member is in 2d array zwy
+!   The solution is in 2d array zwx
+!   The 2d arry zwt and zwz are work space arrays
+!
+!   N.B. the starting vertical index (ikst) is equal to 1 except for
+!   the resolution of tke matrix where surface tke value is prescribed
+!   so that ikstrt=2.
+!!----------------------------------------------------------------------
+
+         ikstp1 = ikst + 1
+         ikenm2 = jpk - 2
+         DO ji = 2, jpim1
+            zwt(ji,ikst) = zwd(ji,ikst)
+         END DO
+         DO jk = ikstp1, jpkm1
+            DO ji = 2, jpim1
+               zwt(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwt(ji,jk-1)
+            END DO
+         END DO
+         DO ji = 2, jpim1
+            zwz(ji,ikst) = zwy(ji,ikst)
+         END DO
+         DO jk = ikstp1, jpkm1
+            DO ji = 2, jpim1
+               zwz(ji,jk) = zwy(ji,jk) - zwi(ji,jk) / zwt(ji,jk-1) * zwz(ji,jk-1)
+            END DO
+         END DO
+         DO ji = 2, jpim1
+            zwx(ji,jpkm1) = zwz(ji,jpkm1) / zwt(ji,jpkm1)
+         END DO
+         DO jk = ikenm2, ikst, -1
+            DO ji = 2, jpim1
+            zwx(ji,jk) =( zwz(ji,jk) - zws(ji,jk) * zwx(ji,jk+1) ) / zwt(ji,jk)
+            END DO
+         END DO
+      
+         ! Normalization to obtain the general momentum trend ua
+         DO jk = 1, jpkm1
+            DO ji = 2, jpim1
+               ua(ji,jj,jk) = ( zwx(ji,jk) - ub(ji,jj,jk) ) / z2dt
+            END DO
+         END DO
+
+
+         ! 2. Vertical diffusion on v
+         ! ---------------------------
+         ! Matrix and second member construction
+         ! bottom boundary condition: only zws must be masked as avmv can take
+         ! non zero value at the ocean bottom depending on the bottom friction
+         ! used (see zdfmix.F)
+         DO jk = 1, jpkm1
+            DO ji = 2, jpim1
+               zcoef = -z2dt/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)
+            END DO
+         END DO
+
+         ! Surface boudary conditions
+         DO ji = 2, jpim1
+            z2dtf = z2dt / ( fse3v(ji,jj,1)*rau0 )
+            zwi(ji,1) = 0.e0
+            zwd(ji,1) = 1. - zws(ji,1)
+            zwy(ji,1) = zwy(ji,1) + z2dtf * tauy(ji,jj)
+         END DO
+
+         ! Matrix inversion starting from the first level
+         ikst = 1
+!!----------------------------------------------------------------------
+!!         ZDF.MATRIXSOLVER
+!!       ********************
+!!----------------------------------------------------------------------
+!! Matrix inversion
+!   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
+!
+!        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
+!        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
+!        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
+!        (        ...               )( ...  ) ( ...  )
+!        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
+!
+!   m is decomposed in the product of an upper and lower triangular
+!   matrix
+!   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
+!   The second member is in 2d array zwy
+!   The solution is in 2d array zwx
+!   The 2d arry zwt and zwz are work space arrays
+!
+!   N.B. the starting vertical index (ikst) is equal to 1 except for
+!   the resolution of tke matrix where surface tke value is prescribed
+!   so that ikstrt=2.
+!!----------------------------------------------------------------------
+
+         ikstp1 = ikst + 1
+         ikenm2 = jpk - 2
+         DO ji = 2, jpim1
+            zwt(ji,ikst) = zwd(ji,ikst)
+         END DO
+         DO jk = ikstp1, jpkm1
+            DO ji = 2, jpim1
+               zwt(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwt(ji,jk-1)
+            END DO
+         END DO
+         DO ji = 2, jpim1
+            zwz(ji,ikst) = zwy(ji,ikst)
+         END DO
+         DO jk = ikstp1, jpkm1
+            DO ji = 2, jpim1
+               zwz(ji,jk) = zwy(ji,jk) - zwi(ji,jk) / zwt(ji,jk-1) * zwz(ji,jk-1)
+            END DO
+         END DO
+         DO ji = 2, jpim1
+            zwx(ji,jpkm1) = zwz(ji,jpkm1) / zwt(ji,jpkm1)
+         END DO
+         DO jk = ikenm2, ikst, -1
+            DO ji = 2, jpim1
+               zwx(ji,jk) =( zwz(ji,jk) - zws(ji,jk) * zwx(ji,jk+1) ) / zwt(ji,jk)
+            END DO
+         END DO
+      
+         ! Normalization to obtain the general momentum trend va
+         DO jk = 1, jpkm1
+            DO ji = 2, jpim1
+               va(ji,jj,jk) = ( zwx(ji,jk) - vb(ji,jj,jk) ) / z2dt
+            END DO
+         END DO
+         !                                             ! ===============
+      END DO                                           !   End of slab
+      !                                                ! ===============
+
+   END SUBROUTINE dyn_zdf_imp_jki
+
+   !!==============================================================================
+END MODULE dynzdf_imp_jki
