  SUBROUTINE dyn2dxctymtm2(u,xeu,xeuv,v,uv,ustar,vstar)

   !&E---------------------------------------------------------------------
   !&E                 ***  dyn2dxctymtm2  ***
   !&E
   !&E ** Purpose : calcul des coefficients de la matrice tridiagonale
   !&E              du systeme lineaire et du second membre relatifs
   !&E              a l equation de conservation de la quantite de mouvement
   !&E              pour le calcul en ligne.
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : dyn2dyx.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays : xeu,xeuv,u,v,uv,ustar,vstar,igdu,hx,htx,hty,
   !&E                     dx,dx_u,dx_v,dy_f,dy_u,fu,hey,hex,surf,phieau,
   !&E                     tfiusr,vis2d,tosx,pa,potgene,buoy2d,
   !&E                     ms,mn,whx
   !&E
   !&E ** Modified variables : y2dxe,y2dv,axe,bxe,cxe,cv,bv,av
   !&E
   !&E ** Reference :
   !&E                T. Duhaut, M. Honnorat, L. Debreu, 2008 : Rapport
   !&E                contrat PREVIMER - Ref : 06/2 210 290 (pages 100-125)
   !&E
   !&E ** History :
   !&E        !  2007-09  (V. Garnier) Commentaires et uniformisation dx,dy
   !&E        !  2008-07  (L. Debreu, T. Duhaut, M. Honnorat) time scheme
   !&E        !  2009-01  (A. Mano, T. Odaka) MPI2D
   !&E        !  2009-03  (V. Garnier) modification of vdusdy and udvsdx term (centered scheme)
   !&E        !  2009-03  (T. Duhaut, V. Garnier) QUICK and QUICKEST scheme (with masks)
   !&E        !  2009-07  (N. Bruneau, F. Dumas) coupling with swan
   !&E        !  2010-07  (L. Debreu, F. Dumas) rewritting of ADI
   !&E        !  2011-03  (F. Dumas) bottom stress parametrisation + remove case 2 and 3 on tidal flats
   !&E        !  2011-10  (L. Pineau, F. Dumas, V. garnier) correction bottom stress parametrisation
   !&E        !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E        !  2011-11  (F. Vandermeirsch) correct spatial extension on which we can estimate quick and quickest scheme
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d,  ONLY : id,ig,surf,dx_u,dx_v,dx,dy_u,dy_v,dy_f,icon,halfdt, &
                          hx,hex,hey,htx,hty,fu,tfiusr,vis2d,dt,h0,frot_lin,  &
                          hminfrot,botstressmax,hmaxfrot,z0b,fwet,fwetp
   USE comvarp2d,  ONLY : v_utmp
   USE comvars3d,  ONLY : rhoref
   USE comcolig2d, ONLY : y2dv,cv,bv,av,igdu,axe,bxe,cxe,y2dxe,mn,ms,whx
   USE comaree,    ONLY : l_genpot,potgene
   USE cometdyn,   ONLY : pa,tosx
   USE comfluxdo,  ONLY : phieau,buoy2d,nbfl
#ifdef key_bottom_stress_wave
   USE waveforcing,ONLY : CfH,Orb_u
#endif
#ifdef key_waveforcing
   USE waveforcing,ONLY : Qx,Qy,Jsmith,Dsmith,Qy_u,k_x,k_y,hlim_calc, &
                          l_wave_brgm_phillips,l_wave_brgm_smith,     &
                          Qx_u,Qx_v,Qy_u,Qy_v,temp_radu,temp_radv
   USE_MPI toolmpi,ONLY : ex_i_rsh,ex_j_rsh
#endif

   IMPLICIT NONE

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in ) :: u,v,uv,ustar,vstar
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in ) :: xeu
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in ) ::xeuv

   !! * Local declarations
   INTEGER                                      :: i,j,k
   REAL(KIND=rsh)                               :: gdt,vit,factft,ft,dusdy,dusdy2,          &
                                                   udusdx,ududx,ududx2,vdusdy,vdudy,vdudy2, &
                                                   v_utmp2,phieautot
#if defined key_dyn_adv_quick || defined key_dyn_adv_quickest
   REAL(KIND=rsh)                               :: uc,um,up,vm,vp
#endif

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

   gdt = grav*dt

#ifdef key_MPI_2D
!$OMP DO schedule(runtime)
#else
!$OMP DO schedule(dynamic,5)
#endif
   DO j=ljminp1,ljmax
    DO i = MAX0(limin,ig(j)-1),MIN0(limax,id(j)-1)
      v_utmp(i,j)=0.25_rsh*(uv(i  ,j  )+uv(i+1,j  )+uv(i  ,j-1)+uv(i+1,j-1))
    END DO
   END DO
!$OMP END DO

#ifdef key_waveforcing
   CALL_MPI ex_i_rsh(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,v_utmp(liminm1:limaxp1,ljminm1:ljmaxp1))
   CALL_MPI ex_j_rsh(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,v_utmp(liminm1:limaxp1,ljminm1:ljmaxp1))
#endif

#ifdef key_MPI_2D
!$OMP DO schedule(runtime)
#else
!$OMP DO schedule(dynamic,5)
#endif

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

   ! Evaluation of axe, bxe, cxe and y2dxe
   !--------------------------------------

   ! 1- terme de diagonale de la matrice tridiagonale
   !    du modele barotrope
   !--------------------------------------------------

   ! 2- terme de sous-diagonale de la matrice tridiagonale
   !    issu de l implicitation de la divergence de flux zonale
   !    (flux a l ouest du point xe(i,j) considere). Remplissage
   !    du tableau hex qui ressert dans les flux de traceurs
   !------------------------------------------------------------

       axe(i,j)=-alphan*halfdt*dy_u(i-1,j)/surf(i,j)*hex(i-1,j)
       bxe(i,j)=fwet(i,j)

   ! 3- terme de sous-diagonale de la matrice tridiagonale
   !    issu de l implicitation de la divergence de flux zonale
   !    (flux a l est du point xe(i,j) considere). Remplissage
   !    du tableau hex qui ressert dans les flux de traceurs
   !-----------------------------------------------------------

       cxe(i,j)=alphan*halfdt*dy_u(i,j)/surf(i,j)*hex(i,j)

   ! 4- terme du vecteur second membre du systeme lineaire
   !    issu de la divergence de flux meridienne
   !    (flux au nord et sud  du point xe(i,j) considere)
   !    Prise en compte de tous les apports d eau dans
   !    la colonne d eau. Remplissage
   !    du tableau hey qui ressert dans les flux de traceurs
   !----------------------------------------------------------

       phieautot=0.0_rsh
       DO k=1,kmax
         phieautot=phieautot+phieau(k,i,j)
       END DO

       y2dxe(i,j)=xeu(i,j)*fwetp(i,j)-halfdt/surf(i,j)*(                 &
                            alphap*( hey(i,j  )*v(i,j  )*dx_v(i,j  )    &
                                   - hey(i,j-1)*v(i,j-1)*dx_v(i,j-1) )  &
                           +alphan*( hey(i,j  )*uv(i,j  )*dx_v(i,j  )   &
                                   - hey(i,j-1)*uv(i,j-1)*dx_v(i,j-1) ) &
                           +alphap*( hex(i  ,j)*u(i  ,j)*dy_u(i  ,j)    &
                                   - hex(i-1,j)*u(i-1,j)*dy_u(i-1,j) )  &
                                   - phieautot )

#ifdef key_waveforcing
       IF (l_wave_brgm_smith .OR. l_wave_brgm_phillips) THEN
         IF (h0(i  ,j  )+xeu(i  ,j  ) >= hlim_calc .AND. &
             h0(i+1,j  )+xeu(i+1,j  ) >= hlim_calc .AND. &
             h0(i-1,j  )+xeu(i-1,j  ) >= hlim_calc)      &
           y2dxe(i,j) = y2dxe(i,j) - halfdt * (Qx(i+1,j)-Qx(i-1,j))/(dx_u(i,j)+dx_u(i-1,j))
         IF (h0(i  ,j  )+xeu(i  ,j  ) >= hlim_calc .AND. &
             h0(i  ,j+1)+xeu(i  ,j+1) >= hlim_calc .AND. &
             h0(i-1,j-1)+xeu(i-1,j-1) >= hlim_calc)      &
           y2dxe(i,j) = y2dxe(i,j) - halfdt * (Qy(i,j+1)-Qy(i,j-1))/(dy_v(i,j)+dy_v(i-1,j))
       END IF
! Attention forme non conservative a priori
#endif

   ! Evaluation of av, bv, cv and y2dv
   !----------------------------------

   ! cas d un point u sans eau de part ni d autre
   ! drying u-point : no water on the right and left

       IF (igdu(i,j) == 4) THEN
         av(i,j)=0.0_rlg
         bv(i,j)=1.0_rlg
         cv(i,j)=0.0_rlg
         y2dv(i,j)=0.0_rlg


   ! cas d un point u avec de l eau de part et d autre
   ! u-point surrounded by water

       ELSE IF (igdu(i,j) == 1) THEN

         v_utmp2=0.25_rsh*(vstar(i,j)+vstar(i+1,j)+vstar(i,j-1)+vstar(i+1,j-1))

       ! sea surface pressure gradient : implicit part
         av(i,j)=-alphan*0.5_rlg*gdt/dx_u(i,j)                ! i-1
         cv(i,j)= alphan*0.5_rlg*gdt/dx_u(i,j)                ! i

#ifdef key_dyn_adv_linear

       ! time evolution (u at n+1/2)
         bv(i,j)=1.0_rlg

       ! rhs : explicit terms
         y2dv(i,j)=                                       &
                 ! time evolution (u at n-1/2)
                 u(i,j)                                 &
                 ! coriolis
                 +fu(i,j)*v_utmp(i,j)*halfdt                 &
                 ! sea surface pressure gradient
                 -0.5_rlg*gdt/dx_u(i,j)*alphap*(xeu(i+1,j)-xeu(i,j))      &
                 ! viscosity (along x)
                 +halfdt*( (vis2d(i+1,j)*(u(i+1,j)-u(i  ,j)))/dx(i+1,j)   &
                          -(vis2d(i  ,j)*(u(i  ,j)-u(i-1,j)))/dx(i  ,j) ) &
                           /dx_u(i,j)

#else

       ! bottom stress
#ifdef key_no_bottom_stress
         ft=0.0_rsh
#elif key_linear_bottom_stress
         !ft=frot_lin/htx(i,j)
         ft=halfdt*frot_lin/MAX(htx(i,j),hminfrot)
#else

#ifdef key_bottom_stress_wave
         ft = halfdt * 2.0_rsh/pi* Orb_u(i,j) * CfH(i,j) / MAX(htx(i,j),hminfrot)
#else
         vit=SQRT(u(i,j)**2+v_utmp(i,j)**2)
         factft=1.0_rsh+vit-botstressmax+ABS(vit-botstressmax)
         ft=factft*0.5_rsh*dt*vit/htx(i,j)*MIN((0.4_rsh/LOG(MAX(MIN(htx(i,j),hmaxfrot),hminfrot)/(exp1*z0b(i,j))))**2,0.1_rsh)
#endif

#endif

#if ! defined key_dyn_adv_none
     ! non-linear terms (udu/dy and vdu/dy) : explicit
#ifdef key_dyn_adv_quickest
         uc=ustar(i,j)
         up=0.5_rsh*(uc+ABS(uc))
         um=0.5_rsh*(uc-ABS(uc))
         vp=0.5_rsh*(v_utmp2+ABS(v_utmp2))
         vm=0.5_rsh*(v_utmp2-ABS(v_utmp2))
         IF (i >= imin+2.AND.i <= imax-2.AND.j <= jmax-2.AND.j >= jmin+3) THEN
           ududx =( (up+um)*( 0.5_rsh*(u(i+1,j)-u(i-1,j))                                                 &
                   -0.5_rsh*(up+um)*(u(i+1,j)+u(i-1,j)-2.0_rsh*u(i,j)) )                                  &
                   +up/6.0_rsh*(1-up*up)*(u(i-2,j)*whx(i-1,j)-3.0_rsh*u(i-1,j)+3.0_rsh*u(i  ,j)-u(i+1,j)) &
                   +um/6.0_rsh*(1-um*um)*(u(i-1,j)-3.0_rsh*u(i  ,j)+3.0_rsh*u(i+1,j)-u(i+2,j)*whx(i+1,j)) )/dx_u(i,j)
           ududx2=( (up+um)*(0.5_rsh*(ustar(i+1,j)-ustar(i-1,j))                       &
                   -0.5_rsh*(up+um)*(ustar(i+1,j)+ustar(i-1,j)-2.0_rsh*ustar(i,j)) )   &
                   +up/6.0_rsh*(1-up*up)*(ustar(i-2,j)*whx(i-1,j)-3.0_rsh*ustar(i-1,j) &
                                          +3.0_rsh*ustar(i  ,j)-ustar(i+1,j))          &
                   +um/6.0_rsh*(1-um*um)*(ustar(i-1,j)-3.0_rsh*ustar(i  ,j)            &
                                          +3.0_rsh*ustar(i+1,j)-ustar(i+2,j)*whx(i+1,j)) )/dx_u(i,j)
           udusdx = (alphan*ududx2 + alphap*ududx)
           vdudy =( (vp+vm)*(0.5_rsh*(u(i,j+1)-u(i,j-1))                                                         & 
                   -0.5_rsh*(vp+vm)*(u(i,j+1)+u(i,j-1)-2.0_rsh*u(i,j)) )                                         &
                   +vp/6.0_rsh*(1-vp*vp)*(u(i,j-2)*ms(i,j-1)*ms(i,j)-3.0_rsh*u(i,j-1)+3.0_rsh*u(i,j  )-u(i,j+1)) &
                   +vm/6.0_rsh*(1-vm*vm)*(u(i,j-1)-3.0_rsh*u(i,j  )+3.0_rsh*u(i,j+1)                             &
                                         -u(i,j+2)*mn(i,j+1)*mn(i,j)) )/dy_u(i,j) 
           vdudy2=( (vp+vm)*( 0.5_rsh*(ustar(i,j+1)-ustar(i,j-1))                                              & 
                   -0.5_rsh*(vp+vm)*(ustar(i,j+1)+ustar(i,j-1)-2.0_rsh*ustar(i,j)) )                           &
                   +vp/6.0_rsh*(1-vp*vp)*(ustar(i,j-2)*ms(i,j-1)*ms(i,j)-3.0_rsh*ustar(i,j-1)                  &
                                          +3.0_rsh*ustar(i,j  )-ustar(i,j+1))                                  &
                   +vm/6.0_rsh*(1-vm*vm)*(ustar(i,j-1)-3.0_rsh*ustar(i,j  )+3.0_rsh*ustar(i,j+1)               &
                                         -ustar(i,j+2)*mn(i,j+1)*mn(i,j)) )/dy_u(i,j) 
           vdusdy = (alphan*vdudy2 + alphap*vdudy)
         ELSE
   ! upwind scheme
           ududx  = up*(u(i,j)-u(i-1,j))/dx(i,j)+um*(u(i+1,j)-u(i,j))/dx(i+1,j)
           ududx2 = up*(ustar(i,j)-ustar(i-1,j))/dx(i,j)+um*(ustar(i+1,j)-ustar(i,j))/dx(i+1,j)
           udusdx = (alphan*ududx2 + alphap*ududx)
           vdudy  = vp*(u(i,j)-u(i,j-1))/dy_f(i,j-1)+vm*(u(i,j+1)-u(i,j))/dy_f(i,j)
           vdudy2 = vp*(ustar(i,j)-ustar(i,j-1))/dy_f(i,j-1)+vm*(ustar(i,j+1)-ustar(i,j))/dy_f(i,j)
           vdusdy = (alphan*vdudy2 + alphap*vdudy)
         END IF
#elif key_dyn_adv_quick
         uc=ustar(i,j)
         up=0.5_rsh*(uc+ABS(uc))
         um=0.5_rsh*(uc-ABS(uc))
         vp=0.5_rsh*(v_utmp2+ABS(v_utmp2))
         vm=0.5_rsh*(v_utmp2-ABS(v_utmp2))
         IF (i >= imin+2.AND.i <= imax-2.AND.j <= jmax-2.AND.j >= jmin+3) THEN
           ududx = (up*u(i-2,j)*whx(i-1,j)+(-7.0_rsh*up-3.0_rsh*um)*u(i-1,j)+                               &
                    (3.0_rsh*up-3.0_rsh*um)*u(i,j)+(3.0_rsh*up+7.0_rsh*um)*u(i+1,j)-um*u(i+2,j)*whx(i+1,j)) &
                   /(8.0_rsh*dx_u(i,j))
           ududx2 = (up*ustar(i-2,j)*whx(i-1,j)+(-7.0_rsh*up-3.0_rsh*um)*ustar(i-1,j)+                                  &
                    (3.0_rsh*up-3.0_rsh*um)*ustar(i,j)+(3.0_rsh*up+7.0_rsh*um)*ustar(i+1,j)-um*ustar(i+2,j)*whx(i+1,j)) &
                 /(8.0_rsh*dx_u(i,j))
           udusdx = (alphan*ududx2 + alphap*ududx)
           vdudy  = (vp*u(i,j-2)*ms(i,j-1)*ms(i,j)+(-7.0_rsh*vp-3.0_rsh*vm)*u(i,j-1)+                              &
                    (3.0_rsh*vp-3.0_rsh*vm)*u(i,j)+(3.0_rsh*vp+7.0_rsh*vm)*u(i,j+1)-vm*u(i,j+2)*mn(i,j+1)*mn(i,j)) &
                    /(8.0_rsh*dy_u(i,j)) 
           vdudy2 = (vp*ustar(i,j-2)*ms(i,j-1)*ms(i,j)+(-7.0_rsh*vp-3.0_rsh*vm)*ustar(i,j-1)+ &
                    (3.0_rsh*vp-3.0_rsh*vm)*ustar(i,j)+(3.0_rsh*vp+7.0_rsh*vm)*ustar(i,j+1)   &
                     -vm*ustar(i,j+2)*mn(i,j+1)*mn(i,j))/(8.0_rsh*dy_u(i,j)) 
           vdusdy = (alphan*vdudy2 + alphap*vdudy)
         ELSE
   ! upwind scheme
           ududx  = up*(u(i,j)-u(i-1,j))/dx(i,j)+um*(u(i+1,j)-u(i,j))/dx(i+1,j)
           ududx2 = up*(ustar(i,j)-ustar(i-1,j))/dx(i,j)+um*(ustar(i+1,j)-ustar(i,j))/dx(i+1,j)
           udusdx = (alphan*ududx2 + alphap*ududx)
           vdudy  = vp*(u(i,j)-u(i,j-1))/dy_f(i,j-1)+vm*(u(i,j+1)-u(i,j))/dy_f(i,j)
           vdudy2 = vp*(ustar(i,j)-ustar(i,j-1))/dy_f(i,j-1)+vm*(ustar(i,j+1)-ustar(i,j))/dy_f(i,j)
           vdusdy = (alphan*vdudy2 + alphap*vdudy)
         END IF
#else
   ! default : centered scheme
         udusdx=ustar(i,j)*                                              &
                ( alphan*(ustar(i+1,j)-ustar(i-1,j))/(dx(i+1,j)+dx(i,j)) &
                 +alphap*(u(i+1,j)-u(i-1,j))/(dx(i+1,j)+dx(i,j)) )
         vdusdy=0.5_rsh*v_utmp2*(                                        &
                ms(i,j)*(  alphan*(ustar(i,j  )-ustar(i,j-1))              &
                         + alphap*(u(i,j  )-u(i,j-1)) )/dy_f(i,j-1)        &
               +mn(i,j)*(  alphan*(ustar(i,j+1)-ustar(i,j  ))              &
                         + alphap*(u(i,j+1)-u(i,j  )) )/dy_f(i,j) )
#endif
#endif

         bv(i,j)=1.0_rlg/(                              &
               ! time evolution (u at n+1/2)
               1.0_rlg                                  &
               ! bottom stress : implicit part
               +ft*cfroi)

       ! rhs : explicit terms
         y2dv(i,j)=                                       &
                 ! time evolution (u at n-1/2)
                 u(i,j)                                 &
                 ! bottom stress (u at n-1/2)
                 -ft*cfroe*u(i,j)                       &
                 ! coriolis
                 +fu(i,j)*v_utmp(i,j)*halfdt                 &
                 ! sea surface pressure gradient
                 -0.5_rlg*gdt/dx_u(i,j)*alphap*(xeu(i+1,j)-xeu(i,j))      &
#if ! defined key_dyn_adv_none
                 ! NL advection : u*v*tan(phi)/R (spherical coordinates)
                 +halfdt*u(i,j)*v_utmp(i,j)*tfiusr(i,j)      &
                 ! non-linear (udu/dx) : explicit
                 -halfdt*udusdx                         &
                 ! non-linear (vdu/dy) : explicit
                 -halfdt*vdusdy                         &
#endif
                 ! viscosity (along x)
                 +halfdt*( (vis2d(i+1,j)*(u(i+1,j)-u(i  ,j)))/dx(i+1,j)   &
                          -(vis2d(i  ,j)*(u(i  ,j)-u(i-1,j)))/dx(i  ,j) ) &
                           /dx_u(i,j)                                     &
                 ! viscosity (along y)
                 +halfdt*( vis2d(i,j  )*mn(i,j)*(u(i,j+1)-u(i,j  ))/dy_f(i,j)     &
                          -vis2d(i,j-1)*ms(i,j)*(u(i,j  )-u(i,j-1))/dy_f(i,j-1) ) &
                           /dy_u(i,j)
#endif

       ! wind stress and atmospherical pressure gradient
         y2dv(i,j)=y2dv(i,j)                              &
                 ! wind stress in x-direction : explicit
                 +halfdt*tosx(i,j)/MAX(htx(i,j),hminfrot)/rhoref &
                 ! atmospherical pressure gradient : explicit
                 -halfdt*(pa(i+1,j)-pa(i,j))/(dx_u(i,j)*rhoref)

       ! generator potential : explicit
         IF (l_genpot) &
           y2dv(i,j)=y2dv(i,j)+halfdt/dx_u(i,j)*(potgene(i+1,j)-potgene(i,j))

       ! internal pressure gradient : explicit
         IF (nbfl /= 0 .AND. icon==2) &
           y2dv(i,j)=y2dv(i,j)+halfdt/dx_u(i,j)*htx(i,j)*(buoy2d(i+1,j)-buoy2d(i,j))

#ifdef key_waveforcing
         IF (l_wave_brgm_phillips) THEN
           y2dv(i,j) = y2dv(i,j) - halfdt * temp_radu(i,j) / MAX(htx(i,j),hminfrot)/rhoref
           y2dv(i,j) = y2dv(i,j) - halfdt * u  (i,j) * (Qx(i+1,j)-Qx(i,j))/dx_u(i,j) / MAX(htx(i,j),hminfrot)
           y2dv(i,j) = y2dv(i,j) - halfdt * 2.0_rsh*Qx_u(i,j) * (u(i+1,j)-u(i-1,j))/(dx(i,j)+dx(i+1,j)) &
                                          / MAX(htx(i,j),hminfrot)
           IF ((hty(i,j-1)+hty(i+1,j-1)) >= 0.0_rsh) &
             y2dv(i,j) = y2dv(i,j) - halfdt *( Qy_u(i,j) * (u  (i,j  )-u  (i,j-1))                           &
                                             + Qx_u(i,j) * (v_utmp(i,j)    -v_utmp(i,j-1)  )                 &
                                             + v_utmp (i,j) * (Qx_u(i,j)-Qx_u(i,j-1)) )/(2.0_rsh*dy_v(i,j-1) &
                                                            * MAX(htx(i,j),hminfrot))
           IF ((hty(i,j  )+hty(i+1,j  )) >= 0.0_rsh) &
             y2dv(i,j) = y2dv(i,j) - halfdt * ( Qy_u(i,j) * (u  (i,j+1)-u  (i,j  ))                         &
                                              + Qx_u(i,j) * (v_utmp(i,j+1)-v_utmp(i,j))                     &
                                              + v_utmp (i,j) * (Qx_u(i,j+1)-Qx_u(i,j)) )/(2.0_rsh*dy_v(i,j) &
                                              * MAX(htx(i,j),hminfrot))  
         ELSE IF (l_wave_brgm_smith) THEN
           y2dv(i,j) = y2dv(i,j) - halfdt * (JSMITH(i+1,j)-JSMITH(i,j)) /dx_u(i,j)
           y2dv(i,j) = y2dv(i,j) +                                                                  &
                       halfdt * (k_x(i,j)+k_x(i+1,j))/2.0_rsh * (Dsmith(i,j)+Dsmith(i+1,j))/2.0_rsh &
                              / MAX(htx(i,j),hminfrot)/rhoref
           y2dv(i,j) = y2dv(i,j) + halfdt * Qy_u(i,j) * (v_utmp(i+1,j)-v_utmp(i-1,j))/(dx(i,j)+dx(i+1,j)) &
                                          / MAX(htx(i,j),hminfrot)
           IF ((hty(i,j-1)+hty(i+1,j-1)) >= 0.0_rsh) &
             y2dv(i,j) = y2dv(i,j) - halfdt * Qy_u(i,j) * (u(i,j  )-u(i,j-1)  )/(2.0_rsh*dy_v(i,j-1)) &
                                            / MAX(htx(i,j),hminfrot)
           IF ((hty(i,j  )+hty(i+1,j  )) >= 0.0_rsh) &
             y2dv(i,j) = y2dv(i,j) - halfdt * Qy_u(i,j) * (u(i,j+1)-u(i,j  )  )/(2.0_rsh*dy_v(i,j)) &
                                            / MAX(htx(i,j),hminfrot)
         ELSE
           y2dv(i,j) = y2dv(i,j) - halfdt * temp_radu(i,j) / MAX(htx(i,j),hminfrot)/rhoref
           y2dv(i,j) = y2dv(i,j) + ft*Qx_u(i,j)/MAX(htx(i,j),hminfrot)
         END IF
#endif

       END IF   ! ends test on mask for wetting and drying

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

  END SUBROUTINE dyn2dxctymtm2
