#if defined key_siggen || defined key_gencoord
  SUBROUTINE dyn3d_wz(xe2,xe3,uz,vz,u,v)
#else
  SUBROUTINE dyn3d_wz(xe2,xe3,uz,vz)
#endif

   !&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---------------------------------------------------------------------

# include "toolcpp.h"

   !! * 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
#ifdef key_diag_physwz
   USE comvars2d, ONLY : hx,hy,dx_u,dy_v
   USE comvars3d, ONLY : wz_phys
   USE comcolig2d,ONLY : mn,ms,me,mo  
#endif
   USE comfluxdo, ONLY : phieau,fxeau,fyeau
#if defined key_tssub_wquickupwind || defined key_tssub_wcompact
   USE comfluxdo, ONLY : fzeau
#endif
#if defined key_siggen || defined key_gencoord
   USE comvars2d, ONLY : hex,hey
   USE comsiggen, ONLY : hzex,hzey
   USE toolgeom,  ONLY : f_hzw
#endif
   USE toolgeom,  ONLY : f_zw_x,f_zw_y

#ifdef key_MPI_2D
   USE toolmpi,   ONLY : ex_j_rsh, ex_i_rsh,ex_peri_x_rsh,ex_peri_y_rsh
#endif

#ifdef key_adapt_imp_vert_adv
   USE parameters,ONLY : rlg
   USE comvars3d, ONLY : f_alpha_T
   USE toolgeom,  ONLY : f_hzu
   USE toolmath,  ONLY : f_alpha 
#ifdef key_adapt_imp_vert_adv_suivi
   USE comvars3d, ONLY : adapt_imp_ts
#endif
#endif 

   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
#if defined key_siggen || defined key_gencoord
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT ( in ) :: u,v
#endif

   !! * Local declarations
   INTEGER                             :: k,i,j
   REAL(KIND=rsh)                      :: d,xet,div,divt
#if defined key_siggen || defined key_gencoord
   REAL(KIND=rsh)                      :: hz3_w
#endif

#ifdef key_adapt_imp_vert_adv
   REAL(kind=rsh)    :: alpha
#endif

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

   PRINT_DBG*, 'enter DYN3D_WZ routine'

!$OMP PARALLEL DEFAULT(NONE) PRIVATE(j,i,k,d,divt,xet,div)  &
#if defined key_siggen || defined key_gencoord
!$OMP& PRIVATE(hz3_w)                                       &
#endif
#ifdef key_adapt_imp_vert_adv
!$OMP& PRIVATE(alpha)                                       &
#endif
!$OMP& SHARED(ljminp2,ljmaxm1,ljmin,ljmax,                  &
!$OMP& limin,limax,l_obc_cycl_x,l_obc_cycl_y,               &
#ifdef key_MPI_2D
!$OMP& liminm1,ljminm1,limaxp1,ljmaxp1,                     &
#endif
!$OMP& ig,id,h0,xe3,xe2,uz,vz,dy_u,dx_v,fwet,               &
#if ! defined key_MPI_2D
!$OMP& jb,jh,                                               &
#endif
#ifdef key_diag_physwz
!$OMP& hx,hy,dx_u,dy_v,wz_phys,mn,ms,me,mo,                 &
#endif
!$OMP& fxeau,fyeau,phieau,dsigu,wz,sigw,surf,halfdt         &
#if defined key_siggen || defined key_gencoord
!$OMP& ,hm                                                  &
#ifdef key_debug
!$OMP& ,hzex,hzey,hex,hey,u,v                               &
#endif
#endif
#if defined key_tssub_wquickupwind || defined key_tssub_wcompact
!$OMP& ,fzeau                                               &
#endif
#ifdef key_agrif
!$OMP& ,imin,imax,jmin,jmax,ishift_obc_west,ishift_obc_east &
!$OMP& ,jshift_obc_south,jshift_obc_north                   &
#else
!$OMP& ,fwetp                                               &
#endif
#ifdef key_adapt_imp_vert_adv
!$OMP& ,f_alpha_T                                           &
#ifdef key_adapt_imp_vert_adv_suivi
!$OMP& ,adapt_imp_ts                                        &
#endif
#endif
!$OMP& )

#ifdef key_MPI_2D
!$OMP DO SCHEDULE(runtime)
#else
!$OMP DO SCHEDULE(dynamic,3)
#endif
   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
#ifdef key_agrif
       xet=-(SUM(dsigu(1:kmax)*fxeau(1:kmax,i  ,j))*dy_u(i  ,j) &
                                      -SUM(dsigu(1:kmax)*fxeau(1:kmax,i-1,j))*dy_u(i-1,j) &
                                      +SUM(dsigu(1:kmax)*fyeau(1:kmax,i  ,j))*dx_v(i  ,j) &
                                      -SUM(dsigu(1:kmax)*fyeau(1:kmax,i,j-1))*dx_v(i,j-1) &
                                      -SUM(phieau(1:kmax,i,j))*halfdt )/(surf(i,j)*halfdt)
                                      
#else
         xet=(xe3(i,j)-xe2(i,j)*fwetp(i,j)/fwet(i,j))/halfdt
#endif
!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)
#if defined key_1DV
#if defined key_siggen || defined key_gencoord
           hz3_w=f_hzw(h0(i,j),xe3(i,j),k,i,j)
!cval continuite OK car 2d mais peut etre que l on peut mieux faire d apres franck
           IF (h0(i,j)+xe3(i,j) <= hm) hz3_w=hm/REAL(kmax,rsh)
           div=-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)/hz3_w
#else
           div=-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

#endif
#else
#if defined key_siggen || defined key_gencoord
           hz3_w=f_hzw(h0(i,j),xe3(i,j),k,i,j)
!cval continuite OK car 2d mais peut etre que l on peut mieux faire d apres franck
           IF (h0(i,j)+xe3(i,j) <= hm) hz3_w=hm/REAL(kmax,rsh)
           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)/hz3_w

#else
           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

#endif
#endif

#ifdef key_diag_physwz
!           wz_phys(k,i,j)=wz(k,i,j)      ! to save vertical sigma coordinates velocity
           wz_phys(k,i,j)=-divt   &       ! to save vertical  "natural" velocity
                    + me(i,j)*mo(i,j)*((f_zw_x(hx(i,j),(xe3(i,j)+xe3(i+1,j))/2.0_rsh,k,i,j)- &
                                        f_zw_x(hx(i-1,j),(xe3(i,j)+xe3(i-1,j))/2.0_rsh,k,i-1,j))/dx_u(i,j))* &
                      (uz(k,i-1,j)+uz(k,i,j)+uz(k+1,i-1,j)+uz(k+1,i,j))/4.0_rsh &
                    + mn(i,j)*ms(i,j)*((f_zw_y(hy(i,j),(xe3(i,j)+xe3(i,j+1))/2.0_rsh,k,i,j)- &
                                        f_zw_y(hy(i,j-1),(xe3(i,j)+xe3(i,j-1))/2.0_rsh,k,i,j-1))/dy_v(i,j))* &
                      (vz(k,i,j-1)+vz(k,i,j)+vz(k+1,i,j-1)+vz(k+1,i,j))/4.0_rsh 

           if (k == kmax-1) then
               wz_phys(kmax,i,j) = xet &     ! wz_phys at the surface 
                    + me(i,j)*mo(i,j)*((xe3(i,j)-xe3(i-1,j))/dx_u(i,j))*(uz(kmax,i-1,j)+uz(kmax,i,j))/2.0_rsh &
                    + mn(i,j)*ms(i,j)*((xe3(i,j)-xe3(i,j-1))/dy_v(i,j))*(vz(kmax,i,j-1)+vz(kmax,i,j))/2.0_rsh
               wz_phys(0,i,j) =  - ( &       ! wz_phys at the bottom cell
                      me(i,j)*mo(i,j)*((hx(i,j)-hx(i-1,j))/dx_u(i,j))*(uz(1,i-1,j)+uz(1,i,j))/2.0_rsh &
                    + mn(i,j)*ms(i,j)*((hy(i,j)-hy(i,j-1))/dy_v(i,j))*(vz(1,i,j-1)+vz(1,i,j))/2.0_rsh )
           endif
#endif

! 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 key_tssub_wcompact)
#if defined key_tssub_wquickupwind || defined key_tssub_wcompact
#if defined key_siggen || defined key_gencoord
           fzeau(k,i,j)=hz3_w*wz(k,i,j)*halfdt
#else
           fzeau(k,i,j)=(h0(i,j)+xe3(i,j))*wz(k,i,j)*halfdt
#endif
#endif

         END DO   ! ends loop on k

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

#ifdef key_debug
#if defined key_siggen || defined key_gencoord
!!!!!!!-----------------------------------------!!!!!!!!!!!!!!!!
!!!!!!! ATTENTION A METTRE A JOUR AVEC LES BONNES VALEURS:
!!!!!!! hzex*uz = fxeau/halfdt   etc.....
!!!!!!! C est pas tres simple attention!!!!
!!!!!!!-----------------------------------------!!!!!!!!!!!!!!!!
       IF (i==2.AND.j==2650) THEN
       PRINT*,(SUM(dsigu(1:kmax)*hzex(1:kmax,i  ,j))*u(i  ,j)* dy_u(i  ,j) &
              -SUM(dsigu(1:kmax)*hzex(1:kmax,i-1,j))*u(i-1,j)* dy_u(i-1,j) &
              +SUM(dsigu(1:kmax)*hzey(1:kmax,i  ,j))*v(i  ,j)*dx_v(i,j  ) &
              -SUM(dsigu(1:kmax)*hzey(1:kmax,i,j-1))*v(i,j-1)*dx_v(i,j-1))/surf(i,j)+xet
       PRINT*,SUM(dsigu(1:kmax)*hzex(1:kmax,i  ,j)*u(i  ,j)*dy_u(i  ,j) &
                 -dsigu(1:kmax)*hzex(1:kmax,i-1,j)*u(i-1,j)*dy_u(i-1,j))
       PRINT*,hex(i  ,j)*u(i  ,j)*dy_u(i  ,j) &
             -hex(i-1,j)*u(i-1,j)*dy_u(i-1,j)
       PRINT*,(hey(i,j  )*v(i,j  )*dx_v(i,j  ) &
              -hey(i,j-1)*v(i,j-1)*dx_v(i,j-1))/surf(i,j)+xet
       PRINT*,SUM(dsigu(1:kmax)*hzex(1:kmax,i  ,j  )),hex(i  ,j  )
       PRINT*,SUM(dsigu(1:kmax)*hzex(1:kmax,i-1,j  )),hex(i-1,j  )
       PRINT*,SUM(dsigu(1:kmax)*hzey(1:kmax,i  ,j  )),hey(i  ,j  )
       PRINT*,SUM(dsigu(1:kmax)*hzey(1:kmax,i  ,j-1)),hey(i  ,j-1)
       END IF
#endif
#endif

#ifdef key_adapt_imp_vert_adv
     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)) &
#ifdef key_siggen
                / ( dsigu(k)*surf(i,j)*f_hzu(h0(i,j),xe2(i,j),k,i,j) )
#else
                / ( dsigu(k)*surf(i,j)*(h0(i,j)+xe2(i,j)) )
#endif
       f_alpha_T(k,i,j) = f_alpha( alpha )
#ifdef key_adapt_imp_vert_adv_suivi
       if (f_alpha_T(k,i,j) < 1.) adapt_imp_ts(k,i,j) = adapt_imp_ts(k,i,j) + 1
#endif
       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)) &
#ifdef key_siggen
                  / ( dsigu(k)*surf(i,j)*f_hzu(h0(i,j),xe2(i,j),k,i,j) )
#else
                  / ( dsigu(k)*surf(i,j)*(h0(i,j)+xe2(i,j)) )
#endif
         f_alpha_T(k,i,j) = f_alpha( alpha )
#ifdef key_adapt_imp_vert_adv_suivi
         if (f_alpha_T(k,i,j) < 1.) adapt_imp_ts(k,i,j) = adapt_imp_ts(k,i,j) + 1
#endif
       END DO
     ELSE
       f_alpha_T(:,i,j)=1.0_rsh
     ENDIF
#endif

     END DO   ! ends loop on i


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

!$OMP SINGLE
   IF (l_obc_cycl_x) THEN
#ifdef key_MPI_2D
     CALL ex_peri_x_rsh(1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wz(:,liminm1:limaxp1,ljminm1:ljmaxp1))
#else
     DO j=ljmin,ljmax
       wz(:,ig(j),j)=wz(:,id(j)-1,j)
       wz(:,id(j),j)=wz(:,ig(j)+1,j)
     END DO
#endif
   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
#ifdef key_MPI_2D
     CALL ex_peri_y_rsh(1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wz(:,liminm1:limaxp1,ljminm1:ljmaxp1))
#else
     DO i=limin,limax
       wz(:,i,jb(i))=wz(:,i,jh(i)-1)
       wz(:,i,jh(i))=wz(:,i,jb(i)+1)
     END DO
#endif
   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_MPI ex_i_rsh(-1,1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wz(:,liminm1:limaxp1,ljminm1:ljmaxp1))
   CALL_MPI ex_j_rsh(-1,1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wz(:,liminm1:limaxp1,ljminm1:ljmaxp1))


   PRINT_DBG*, 'exit DYN3D_WZ routine'

  END SUBROUTINE dyn3d_wz
