
  SUBROUTINE dyn3d_wz(xe2,xe3,uz,vz)


!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE flx_cd  ***
!&E
!&E ** Purpose : estimate of the vertical velocity from the vertical
!&E              integration (from the bottom to the surface)
!&E              of the local continuity equation
!&E
!&E ** Description : prognostic estimate
!&E
!&E ** Called by : step3dxy.F90,step3dyx.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : h0,xe3,xe2,fxeau,fyeau,phieau,dy_u,dx_v,surf
!&E
!&E ** Modified variables : wz,fzeau
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  19887-02 (J.-C. Salomon, P. Lazure)  Original code
!&E       !  2007-09  (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
!&E       !  2008-08  (T. Odaka) MPI
!&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
!&E       !  2009-04  (V. Garnier) MPI2D and removing twisted obc part
!&E       !  2010-04  (T. Odaka) MPI2D cyclic
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2012-06  (P. Garreau) Estimate of physical vertical velocity (key_diag_physwz)
!&E       !  2015-03  (P. Garreau) Correction of physical vertical velocity (key_diag_physwz)
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters,ONLY : imin,imax,jmin,jmax,kmax,ljmin,ljmax,rsh,            &
                         ljminm1,ljminp2,ljmaxm1,limin,limax,limaxp1,ljmaxp1, &
                         liminm2,limaxp2,ljminm2,ljmaxp2,liminm1,ljminm1,     &
                         ishift_obc_west,ishift_obc_east,jshift_obc_south,jshift_obc_north
   USE comvars2d, ONLY : ig,id,h0,dy_u,dx_v,l_obc_cycl_y,l_obc_cycl_x,jb,jh,surf,halfdt,hm
   USE comvars2d, ONLY : fwet,fwetp
   USE comvars3d, ONLY : dsigu,wz,sigw

   USE comfluxdo, ONLY : phieau,fxeau,fyeau

   USE comfluxdo, ONLY : fzeau


   USE toolgeom,  ONLY : f_zw_x,f_zw_y


   USE toolmpi,   ONLY : ex_j_rsh, ex_i_rsh,ex_peri_x_rsh,ex_peri_y_rsh



   USE parameters,ONLY : rlg
   USE comvars3d, ONLY : f_alpha_T
   USE toolgeom,  ONLY : f_hzu
   USE toolmath,  ONLY : f_alpha 



   IMPLICIT NONE

!! * Arguments
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT ( in ) :: xe2
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT ( in ) :: xe3
   REAL(KIND=rsh), DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2), INTENT ( in )   :: uz,vz


!! * Local declarations
   INTEGER                             :: k,i,j
   REAL(KIND=rsh)                      :: d,xet,div,divt



   REAL(kind=rsh)    :: alpha


!!----------------------------------------------------------------------
!! * Executable part

   !PRINT*, 'enter DYN3D_WZ routine'

!$OMP PARALLEL DEFAULT(NONE) PRIVATE(j,i,k,d,divt,xet,div)  &


!$OMP& PRIVATE(alpha)                                       &

!$OMP& SHARED(ljminp2,ljmaxm1,ljmin,ljmax,                  &
!$OMP& limin,limax,l_obc_cycl_x,l_obc_cycl_y,               &

!$OMP& liminm1,ljminm1,limaxp1,ljmaxp1,                     &

!$OMP& ig,id,h0,xe3,xe2,uz,vz,dy_u,dx_v,fwet,               &


!$OMP& fxeau,fyeau,phieau,dsigu,wz,sigw,surf,halfdt         &


!$OMP& ,fzeau                                               &


!$OMP& ,fwetp                                               &


!$OMP& ,f_alpha_T                                           &


!$OMP& )


!$OMP DO SCHEDULE(runtime)

   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)

       d=h0(i,j)+xe3(i,j)

       IF (d > 0.0_rsh) THEN
         divt=0.0_rsh

         xet=(xe3(i,j)-xe2(i,j)*fwetp(i,j)/fwet(i,j))/halfdt

!cval         DO  k=1,kmax
! Inutile de calculer wz(kmax). On neglige la vitesse verticale induite par la pluie
! et les fleuves. Ailleurs dans le code, les discretisation suppose que wz(kmax)=0
! et wz(kmax) n est jamais appele
         DO  k=1,kmax-1

! estimates vertical velocity (sigma)


           div=(fxeau(k,i,j)*dy_u(i,j)-fxeau(k,i-1,j)*dy_u(i-1,j))/surf(i,j)/halfdt &
              +(fyeau(k,i,j)*dx_v(i,j)-fyeau(k,i,j-1)*dx_v(i,j-1))/surf(i,j)/halfdt &
              -phieau(k,i,j)/(surf(i,j)*fwet(i,j))/dsigu(k)
           divt=divt+div*dsigu(k)
           wz(k,i,j)=-(xet*(sigw(k)+1.0_rsh)+divt)/d






! Rq cval : xet*(sigw(k)+1.0_rsh) correspond a intregation verticale de -1 a sigw du terme xet

! time integration of the vertical velocity wz ( hz*wz if key_tssub_wquickupwind or 1)


           fzeau(k,i,j)=(h0(i,j)+xe3(i,j))*wz(k,i,j)*halfdt



         END DO   ! ends loop on k

       ELSE
         wz(:,i,j)=0.0_rsh
       END IF   ! ends test on water depth




     IF (h0(i,j)+xe2(i,j) > 0.0_rsh) THEN
       k=1
       alpha = ((max(fxeau(k,i,j),0.0_rlg)*dy_u(i,j)-min(fxeau(k,i-1,j),0.0_rlg)*dy_u(i-1,j) &
                +max(fyeau(k,i,j),0.0_rlg)*dx_v(i,j)-min(fyeau(k,i,j-1),0.0_rlg)*dx_v(i,j-1) ) * dsigu(k) &
                +max(fzeau(k,i,j),0.0_rlg)*surf(i,j)) &

                / ( dsigu(k)*surf(i,j)*(h0(i,j)+xe2(i,j)) )

       f_alpha_T(k,i,j) = f_alpha( alpha )

       DO k=2,kmax
         alpha = ((max(fxeau(k,i,j),0.0_rlg)*dy_u(i,j)-min(fxeau(k,i-1,j),0.0_rlg)*dy_u(i-1,j) &
                  +max(fyeau(k,i,j),0.0_rlg)*dx_v(i,j)-min(fyeau(k,i,j-1),0.0_rlg)*dx_v(i,j-1) ) * dsigu(k) &
                  +(max(fzeau(k,i,j),0.0_rlg)-min(fzeau(k-1,i,j),0.0_rlg))*surf(i,j)) &

                  / ( dsigu(k)*surf(i,j)*(h0(i,j)+xe2(i,j)) )

         f_alpha_T(k,i,j) = f_alpha( alpha )

       END DO
     ELSE
       f_alpha_T(:,i,j)=1.0_rsh
     ENDIF


     END DO   ! ends loop on i


   END DO   ! ends loop on j
!$OMP END DO

!$OMP SINGLE
   IF (l_obc_cycl_x) THEN

     CALL ex_peri_x_rsh(1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wz(:,liminm1:limaxp1,ljminm1:ljmaxp1))

   ELSE
     IF (limax==imax) wz(:,imax-ishift_obc_east,:)=wz(:,imax-1-ishift_obc_east,:)
     IF (limin==imin) wz(:,imin+1+ishift_obc_west,:)=wz(:,imin+2+ishift_obc_west,:)
   END IF
   IF (l_obc_cycl_y) THEN

     CALL ex_peri_y_rsh(1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wz(:,liminm1:limaxp1,ljminm1:ljmaxp1))

   ELSE
     IF (ljmin==jmin) wz(:,limin:limax,jmin+1+jshift_obc_south)=wz(:,limin:limax,jmin+2+jshift_obc_south)
     IF (ljmax==jmax) wz(:,limin:limax,jmax-jshift_obc_north)=wz(:,limin:limax,jmax-1-jshift_obc_north)
   END IF
!$OMP END SINGLE

!$OMP END PARALLEL

   CALL ex_i_rsh(-1,1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wz(:,liminm1:limaxp1,ljminm1:ljmaxp1))
   CALL ex_j_rsh(-1,1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wz(:,liminm1:limaxp1,ljminm1:ljmaxp1))


   !PRINT*, 'exit DYN3D_WZ routine'

  END SUBROUTINE dyn3d_wz
