MODULE wzvmod
   !!==============================================================================
   !!                       ***  MODULE  wzvmod  ***
   !! Ocean diagnostic variable : vertical velocity
   !!==============================================================================
   !! History :   5.0  !  90-10  (C. Levy, G. Madec)  Original code
   !!             7.0  !  96-01  (G. Madec)  Statement function for e3
   !!             8.5  !  02-07  (G. Madec)  Free form, F90
   !!----------------------------------------------------------------------
   !!   wzv        : Compute the vertical velocity
   !!----------------------------------------------------------------------
   !! * Modules used
   USE oce             ! ocean dynamics and tracers variables
   USE dom_oce         ! ocean space and time domain variables 
   USE sbc_oce         ! surface boundary condition: ocean
   USE domvvl          ! Variable volume
   USE in_out_manager  ! I/O manager
   USE prtctl          ! Print control
   USE phycst
   USE lbclnk          ! ocean lateral boundary condition (or mpp link)

   IMPLICIT NONE
   PRIVATE

   !! * Routine accessibility
   PUBLIC wzv          ! routine called by step.F90 and inidtr.F90

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

CONTAINS

   SUBROUTINE wzv( kt )
      !!----------------------------------------------------------------------
      !!                    ***  ROUTINE wzv  ***
      !!
      !! ** Purpose :   Compute the now vertical velocity after the array swap
      !!
      !! ** Method  :   Using the incompressibility hypothesis, the vertical
      !!      velocity is computed by integrating the horizontal divergence 
      !!      from the bottom to the surface.
      !!        The boundary conditions are w=0 at the bottom (no flux) and,
      !!      in regid-lid case, w=0 at the sea surface.
      !!
      !! ** action  :   wn array : the now vertical velocity
      !!----------------------------------------------------------------------
      !! * Arguments
      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index

      !! * Local declarations
      INTEGER  ::           jk           ! dummy loop indices
      !! Variable volume
      INTEGER  ::   ji, jj               ! dummy loop indices
      REAL(wp) ::   z2dt, zraur          ! temporary scalar
      REAL(wp), DIMENSION (jpi,jpj) ::   zssha, zun, zvn, zhdiv
      !!----------------------------------------------------------------------

      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'wzv     : vertical velocity from continuity eq.'
         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 

         ! bottom boundary condition: w=0 (set once for all)
         wn(:,:,jpk) = 0.e0
      ENDIF

      IF( lk_vvl ) THEN                ! Variable volume
         !
         z2dt = 2. * rdt                                       ! time step: leap-frog
         IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest
         zraur  = 1. / rauw

         ! Vertically integrated quantities
         ! --------------------------------
         zun(:,:) = 0.e0
         zvn(:,:) = 0.e0
         !
         DO jk = 1, jpkm1             ! Vertically integrated transports (now)
            zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk)
            zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
         END DO

         ! Horizontal divergence of barotropic transports
         !--------------------------------------------------
         DO jj = 2, jpjm1
            DO ji = 2, jpim1   ! vector opt.
               zhdiv(ji,jj) = (  e2u(ji  ,jj  ) * zun(ji  ,jj  )     &
                  &            - e2u(ji-1,jj  ) * zun(ji-1,jj  )     &
                  &            + e1v(ji  ,jj  ) * zvn(ji  ,jj  )     &
                  &            - e1v(ji  ,jj-1) * zvn(ji  ,jj-1) )   &
                  &           / ( e1t(ji,jj) * e2t(ji,jj) )
            END DO
         END DO

#if defined key_obc && ( key_dynspg_exp || key_dynspg_ts )
         ! open boundaries (div must be zero behind the open boundary)
         !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east
         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west
         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
#endif

         CALL lbc_lnk( zhdiv, 'T', 1. )

         ! Sea surface elevation time stepping
         ! -----------------------------------
         zssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1)

         ! Vertical velocity computed from bottom
         ! --------------------------------------
         DO jk = jpkm1, 1, -1
            wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) &
              &        - ( zssha(:,:) - sshb(:,:) ) * fsve3t(:,:,jk) * mut(:,:,jk) / z2dt
         END DO

      ELSE                             ! Fixed volume 

         ! Vertical velocity computed from bottom
         ! --------------------------------------
         DO jk = jpkm1, 1, -1
            wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk)
         END DO
      
      ENDIF 

      IF(ln_ctl)   CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 -   : ', mask1=wn)

   END SUBROUTINE wzv

   !!======================================================================
END MODULE wzvmod