MODULE dynzdf_exp
   !!==============================================================================
   !!                     ***  MODULE  dynzdf_exp  ***
   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
   !!==============================================================================
   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal
   !!   NEMO     0.5  !  2002-08  (G. Madec)  F90: Free form and module
   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   dyn_zdf_exp  : update the momentum trend with the vertical diffu-
   !!                  sion using an explicit time-stepping scheme.
   !!----------------------------------------------------------------------
   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 sbc_oce         ! surface boundary condition: ocean
   USE lib_mpp         ! MPP library
   USE in_out_manager  ! I/O manager
   USE lib_mpp         ! MPP library

   IMPLICIT NONE
   PRIVATE

   PUBLIC   dyn_zdf_exp   ! called by step.F90

   !! * Control permutation of array indices
#  include "oce_ftrans.h90"
#  include "dom_oce_ftrans.h90"
#  include "zdf_oce_ftrans.h90"
#  include "sbc_oce_ftrans.h90"
   
   !! * Substitutions
#  include "domzgr_substitute.h90"
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
   !! $Id$
   !! Software governed by the CeCILL licence     (NEMOGCM/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 (averaged over kt-1/2 & kt+1/2)
      !!      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
      !!---------------------------------------------------------------------
      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
      USE oce     , ONLY:   zwx => ta       , zwy => sa         ! (ta,sa) used as 3D workspace
      USE wrk_nemo, ONLY:   zwz => wrk_3d_1 , zww => wrk_3d_2   ! 3D workspace
      !! DCSE_NEMO: need additional directives for renamed module variables
!FTRANS zwx :I :I :z
!FTRANS zwy :I :I :z
!FTRANS zwz :I :I :z
!FTRANS zww :I :I :z
      !
      INTEGER , INTENT(in) ::   kt     ! ocean time-step index
      REAL(wp), INTENT(in) ::   p2dt   ! time-step 
      !
      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
      REAL(wp) ::   zrau0r, zlavmr, zua, zva   ! local scalars
      !!----------------------------------------------------------------------

      IF( wrk_in_use(3, 1,2) ) THEN
         CALL ctl_stop('dyn_zdf_exp: requested workspace arrays unavailable')   ;   RETURN
      ENDIF

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

      zrau0r = 1. / rau0               ! Local constant initialization
      zlavmr = 1. / REAL( nn_zdfexp )


      DO jj = 2, jpjm1                 ! Surface boundary condition
         DO ji = 2, jpim1
            zwy(ji,jj,1) = ( utau_b(ji,jj) + utau(ji,jj) ) * zrau0r
            zww(ji,jj,1) = ( vtau_b(ji,jj) + vtau(ji,jj) ) * zrau0r
         END DO  
      END DO  
#if defined key_z_first
      DO jj = 2, jpjm1                 ! Initialization of x, z and contingently trends array
         DO ji = 2, jpim1
            DO jk = 1, jpk
#else
      DO jk = 1, jpk                   ! Initialization of x, z and contingently trends array
         DO jj = 2, jpjm1 
            DO ji = 2, jpim1
#endif
               zwx(ji,jj,jk) = ub(ji,jj,jk)
               zwz(ji,jj,jk) = vb(ji,jj,jk)
            END DO  
         END DO  
      END DO  
      !
      DO jl = 1, nn_zdfexp             ! Time splitting loop
         !
#if defined key_z_first
         DO jj = 2, jpjm1 
            DO ji = 2, jpim1
               DO jk = 2, jpk                ! First vertical derivative
#else
         DO jk = 2, jpk                      ! First vertical derivative
            DO jj = 2, jpjm1 
               DO ji = 2, jpim1
#endif
                  zwy(ji,jj,jk) = avmu(ji,jj,jk) * ( zwx(ji,jj,jk-1) - zwx(ji,jj,jk) ) / fse3uw(ji,jj,jk) 
                  zww(ji,jj,jk) = avmv(ji,jj,jk) * ( zwz(ji,jj,jk-1) - zwz(ji,jj,jk) ) / fse3vw(ji,jj,jk)
               END DO  
            END DO  
         END DO  
#if defined key_z_first
         DO jj = 2, jpjm1 
            DO ji = 2, jpim1
               DO jk = 1, jpkm1              ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp
#else
         DO jk = 1, jpkm1                    ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp
            DO jj = 2, jpjm1 
               DO ji = 2, jpim1
#endif
                  zua = zlavmr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) / fse3u(ji,jj,jk)
                  zva = zlavmr * ( zww(ji,jj,jk) - zww(ji,jj,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,jj,jk) = zwx(ji,jj,jk) + p2dt * zua * umask(ji,jj,jk)
                  zwz(ji,jj,jk) = zwz(ji,jj,jk) + p2dt * zva * vmask(ji,jj,jk)
               END DO  
            END DO  
         END DO  
         !
      END DO                           ! End of time splitting
      !
      IF( wrk_not_released(3, 1,2) )   CALL ctl_stop('dyn_zdf_exp: failed to release workspace arrays')
      !
   END SUBROUTINE dyn_zdf_exp

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