MODULE dynzdf_exp
   !!==============================================================================
   !!                     ***  MODULE  dynzdf_exp  ***
   !! 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
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   dyn_zdf_exp  : update the momentum trend with the vertical diffu-
   !!                  sion using an explicit time-stepping scheme.
   !!----------------------------------------------------------------------
   !! * 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

   IMPLICIT NONE
   PRIVATE

   !! * Routine accessibility
   PUBLIC dyn_zdf_exp    ! called by step.F90

   !! * Substitutions
#  include "domzgr_substitute.h90"
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
   !! $Header$ 
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE dyn_zdf_exp( kt, p2dt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dyn_zdf_exp  ***
      !!                   
      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
      !!
      !! ** Method  :   Explicit forward time stepping with a time splitting
      !!      technique. The vertical diffusion of momentum is given by:
      !!         diffu = dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ub) )
      !!      Surface boundary conditions: wind stress input
      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F90)
      !!      Add this trend to the general trend ua :
      !!         ua = ua + dz( avmu dz(u) )
      !!
      !! ** Action : - Update (ua,va) with the vertical diffusive trend
      !!---------------------------------------------------------------------
      !! * Arguments
      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, zua, zva                   ! temporary scalars
      REAL(wp), DIMENSION(jpi,jpk) ::   zwx, zwy, zwz, zww     ! temporary workspace arrays
      !!----------------------------------------------------------------------

      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'dyn_zdf_exp : vertical momentum diffusion explicit operator'
         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
      ENDIF

      ! Local constant initialization
      ! -----------------------------
      zrau0r = 1. / rau0                                   ! inverse of the reference density
      zlavmr = 1. / float( n_zdfexp )                      ! inverse of the number of sub time step

      !                                                ! ===============
      DO jj = 2, jpjm1                                 !  Vertical slab
         !                                             ! ===============

         ! Surface boundary condition
         DO ji = 2, jpim1
            zwy(ji,1) = taux(ji,jj) * zrau0r
            zww(ji,1) = tauy(ji,jj) * zrau0r
         END DO  

         ! Initialization of x, z and contingently trends array
         DO jk = 1, jpk
            DO ji = 2, jpim1
               zwx(ji,jk) = ub(ji,jj,jk)
               zwz(ji,jk) = vb(ji,jj,jk)
            END DO  
         END DO  

         ! Time splitting loop
         DO jl = 1, n_zdfexp

            ! First vertical derivative
            DO jk = 2, jpk
               DO ji = 2, jpim1
                  zwy(ji,jk) = avmu(ji,jj,jk) * ( zwx(ji,jk-1) - zwx(ji,jk) ) / fse3uw(ji,jj,jk) 
                  zww(ji,jk) = avmv(ji,jj,jk) * ( zwz(ji,jk-1) - zwz(ji,jk) ) / fse3vw(ji,jj,jk)
               END DO  
            END DO  

            ! Second vertical derivative and trend estimation at kt+l*rdt/n_zdfexp
            DO jk = 1, jpkm1
               DO ji = 2, jpim1
                  zua = zlavmr*( zwy(ji,jk) - zwy(ji,jk+1) ) / fse3u(ji,jj,jk)
                  zva = zlavmr*( zww(ji,jk) - zww(ji,jk+1) ) / fse3v(ji,jj,jk)
                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua
                  va(ji,jj,jk) = va(ji,jj,jk) + zva

                  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  

         END DO  

         !                                             ! ===============
      END DO                                           !   End of slab
      !                                                ! ===============

   END SUBROUTINE dyn_zdf_exp

   !!==============================================================================
END MODULE dynzdf_exp