  SUBROUTINE dyn2d_uvstar(u,v,xe,ustar,vstar)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE flx_cd  ***
   !&E
   !&E ** Purpose : evaluation of a first guess for u and v
   !&E
   !&E ** Description : explicit estimation, except for bottom stress
   !&E                  (semi-implicit according to cfroe,cfroi)
   !&E
   !&E ** Called by : dyn2dxy.F90, dyn2dyx.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays : u,v,xe,htx,hty,dx,dy,dx_u,dy_u,dx_v,dy_v,
   !&E                     dx_f,dy_f,fu,fv,igdu,igdv,tfiusr,tfivsr,
   !&E                     vis2d,tosx,tosy,pa,potgene,buoy2d
   !&E                     me,mo,ms,mn,whx,why
   !&E
   !&E ** Modified variables : ustar, vstar
   !&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        !  2008-07  (L. Debreu, T. Duhaut, M. Honnorat) time scheme
   !&E        !  2008-07  (T. Odaka) MPI
   !&E        !  2009-01  (V. Garnier) OMP
   !&E        !  2009-01  (A. Mano T. Odaka) MPI2D
   !&E        !  2009-03  (V. Garnier) modification of vdusdy and udvsdx terms (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        !  2011-03  (F. Dumas)  bottom stress parametrisation + remove case 2 and 3
   !&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                                       plus big bug line 140
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d,  ONLY : botstressmax,z0b,tfiusr,tfivsr,halfdt,    &
                          htx,hx,hy,hty,hey,hex,fu,fv,vis2d,        & 
                          dy,dx,dy_f,dy_u,dx_u,dx_v,dy_v,dx_f,      &
                          jb,jh,ig,id,hminfrot,hmaxfrot,icon,frot_lin
   USE comvarp2d,  ONLY : u_vtmp,v_utmp
   USE comcolig2d, ONLY : igdu,igdv,mn,ms,me,mo,whx,why
   USE  comfluxdo, ONLY : nbfl,buoy2d
   USE comaree,    ONLY : l_genpot,potgene
   USE cometdyn,   ONLY : tosx,pa,tosy
   USE comvars3d,  ONLY : rhoref
#ifdef key_bottom_stress_wave
   USE waveforcing,ONLY : CfH,Orb_u,Orb_v
#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 ) :: xe
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in ) :: u,v
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( out ) :: ustar,vstar

   !! * Local declarations
   INTEGER              :: i,j,k,i1,i2,j1,j2
   REAL(KIND=rsh)       :: vit,factft,ft,udusdx,vdusdy,dusdy, &
                           udvsdx,vdvsdy,dvsdx
#if defined key_dyn_adv_quick || defined key_dyn_adv_quickest
   REAL(KIND=rsh)       ::uc,vc,um,up,vm,vp 
#endif

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

   PRINT_DBG*, 'enter DYN2D_UVSTAR routine'

   ! -----------------
   ! Estimate of ustar
   ! -----------------

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp1,ljmax
     DO i=MAX0(ig(j)-1,limin),MIN0(id(j)-1,limax)
       v_utmp(i,j)=0.25_rsh*(v(i  ,j)+v(i+1,j)+v(i  ,j-1)+v(i+1,j-1))
     END DO
   END DO
!$OMP END DO

#ifdef key_waveforcing
OMPMPI barrier
OMPMPI master
   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))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(v_utmp)
#endif

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(ig(j),limin),MIN0(id(j)-1,limax)

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

     ! bottom stress factor
#ifdef key_no_bottom_stress
         ft=0.0_rsh
#elif key_linear_bottom_stress
         !ft=frot_lin/htx(i,j)
         ft=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*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

     ! non-linear terms (udu/dy and vdu/dy) : explicit
#if ! defined key_dyn_adv_none
#if ! defined key_dyn_adv_linear
#ifdef key_dyn_adv_quickest
         uc=u(i,j)
         up=0.5_rsh*(uc+ABS(uc))
         um=0.5_rsh*(uc-ABS(uc))
         vp=0.5_rsh*(v_utmp(i,j)+ABS(v_utmp(i,j)))
         vm=0.5_rsh*(v_utmp(i,j)-ABS(v_utmp(i,j)))
         IF (i >= imin+2.AND.i <= imax-2.AND.j <= jmax-2.AND.j >= jmin+3) THEN
           udusdx=( (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)
           vdusdy=( (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)
         ELSE
   ! upwind scheme
           udusdx  = up*(u(i,j)-u(i-1,j))/dx(i,j)+um*(u(i+1,j)-u(i,j))/dx(i+1,j)
           vdusdy  = 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)
         END IF
#elif key_dyn_adv_quick
         uc=u(i,j)
         up=0.5_rsh*(uc+ABS(uc))
         um=0.5_rsh*(uc-ABS(uc))
         vp=0.5_rsh*(v_utmp(i,j)+ABS(v_utmp(i,j)))
         vm=0.5_rsh*(v_utmp(i,j)-ABS(v_utmp(i,j)))
         IF (i >= imin+2.AND.i <= imax-2.AND.j <= jmax-2.AND.j >= jmin+3) THEN
           udusdx  =(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))
           vdusdy  =(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))
         ELSE
   ! upwind scheme
           udusdx  = up*(u(i,j)-u(i-1,j))/dx(i,j)+um*(u(i+1,j)-u(i,j))/dx(i+1,j)
           vdusdy  = 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)
         END IF
#else
   ! default : centered scheme
         udusdx=u(i,j)*(u(i+1,j)-u(i-1,j))/(dx(i,j)+dx(i+1,j))
         vdusdy=0.5_rsh*v_utmp(i,j)*                                         &
                (  ms(i,j)*(u(i,j  )-u(i,j-1))/dy_f(i,j-1)                   &
                 + mn(i,j)*(u(i,j+1)-u(i,j  ))/dy_f(i,j  ) )
#endif
#endif
#endif

         ustar(i,j)=u(i,j)-halfdt*(                                          &
#if ! defined key_dyn_adv_none
#if ! defined key_dyn_adv_linear
                   ! non-linear advection : explicit
                   udusdx + vdusdy - u(i,j)*v_utmp(i,j)*tfiusr(i,j)          &
#endif
#endif
                   ! bottom stress : explicit part
                  +ft*cfroe*u(i,j)                                           & 
                   ! viscosity (along x) : explicit
                  -((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)  &
                   ! sea surface height gradient (along x) : explicit
                  +grav*(xe(i+1,j)-xe(i,j))/dx_u(i,j)                        &
                   ! coriolis term : explicit
                  -fu(i,j)*v_utmp(i,j)                                       &
                   ! viscosity (along y) : explicit
                  +( vis2d(i,j-1)*ms(i,j)*(u(i,j  )-u(i,j-1))/dy_f(i,j-1)    &
                    -vis2d(i,j  )*mn(i,j)*(u(i,j+1)-u(i,j  ))/dy_f(i,j  ) )  &
                    /dy_u(i,j) )

      ! wind stress and atmospherical pressure gradient
         ustar(i,j)=ustar(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) &
           ustar(i,j)=ustar(i,j)+halfdt/dx_u(i,j)*(potgene(i+1,j)-potgene(i,j))

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

       ! bottom stress : implicit part
         ustar(i,j) = ustar(i,j)/(1.0_rlg + ft*cfroi*halfdt)

#ifdef key_waveforcing
         IF (l_wave_brgm_phillips) THEN
           ustar(i,j) = ustar(i,j) - halfdt * temp_radu(i,j) / MAX(htx(i,j),hminfrot)/rhoref
           ustar(i,j) = ustar(i,j) - halfdt * u  (i,j) * (Qx(i+1,j)-Qx(i,j))/dx_u(i,j) / MAX(htx(i,j),hminfrot)
           ustar(i,j) = ustar(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) THEN
             ustar(i,j) = ustar(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))
           END IF
           IF ((hty(i,j  )+hty(i+1,j  )) >= 0.0_rsh) THEN
             ustar(i,j) = ustar(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))  
           END IF
         ELSE IF (l_wave_brgm_smith) THEN
           ustar(i,j) = ustar(i,j) - halfdt * (JSMITH(i+1,j)-JSMITH(i,j)) /dx_u(i,j)
           ustar(i,j) = ustar(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
           ustar(i,j) = ustar(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) &
             ustar(i,j) = ustar(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) &
             ustar(i,j) = ustar(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
           ustar(i,j) = ustar(i,j) - halfdt * temp_radu(i,j) / MAX(htx(i,j),hminfrot)/rhoref
           ustar(i,j) = ustar(i,j) + ft*Qx_u(i,j)/MAX(htx(i,j),hminfrot)
         END IF
#endif

       ENDIF   ! ends test on j

     END DO
   END DO
!$OMP END DO

   ! -----------------
   ! Estimate of vstar
   ! -----------------

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljmin,ljmax
     DO i=liminp1,limax
       IF (j >= jb(i)-1 .AND. j <= jh(i)-1) &
         u_vtmp(i,j)=0.25_rsh*(u(i,j  )+u(i,j+1)+u(i-1,j  )+u(i-1,j+1))
     END DO
   END DO
!$OMP END DO

#ifdef key_waveforcing
OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,u_vtmp(liminm1:limaxp1,ljminm1:ljmaxp1))
   CALL_MPI ex_j_rsh(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,u_vtmp(liminm1:limaxp1,ljminm1:ljmaxp1))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(u_vtmp)
#endif

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp1,ljmax
     DO i=liminp2,limaxm1
       IF (j >= jb(i) .AND. j <= jh(i)-1) THEN
         IF (igdv(i,j) == 1) THEN

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

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

#endif

     ! non-linear terms (udu/dy and vdu/dy) : explicit
#if ! defined key_dyn_adv_none
#if ! defined key_dyn_adv_linear
#ifdef key_dyn_adv_quickest
   ! QUICKEST scheme
           vc=v(i,j)
           vp=0.5_rsh*(vc+ABS(vc))
           vm=0.5_rsh*(vc-ABS(vc))
           up=0.5_rsh*(u_vtmp(i,j)+ABS(u_vtmp(i,j)))
           um=0.5_rsh*(u_vtmp(i,j)-ABS(u_vtmp(i,j)))
           IF (i >= imin+3.AND.i <= imax-2.AND.j <= jmax-2.AND.j >= jmin+2) THEN
             vdvsdy=( (vp+vm)*( 0.5_rsh*(v(i,j+1)-v(i,j-1))                              &
                     -0.5_rsh*(vp+vm)*(v(i,j+1)+v(i,j-1)-2.0_rsh*v(i,j)) )               &
                     +vp/6.0_rsh*(1-vp*vp)*(v(i,j-2)*why(i,j-1)-3.0_rsh*v(i,j-1)+3.0_rsh*v(i,j  )-v(i,j+1))          &
                     +vm/6.0_rsh*(1-vm*vm)*(v(i,j-1)-3.0_rsh*v(i,j  )+3.0_rsh*v(i,j+1)-v(i,j+2)*why(i,j+1)) )/dy(i,j)
             udvsdx=( (up+um)*( 0.5_rsh*(v(i+1,j)-v(i-1,j))                              &
                     -0.5_rsh*(up+um)*(v(i+1,j)+v(i-1,j)-2.0_rsh*v(i,j)) )               &
                     +up/6.0_rsh*(1-up*up)*(v(i-2,j)*mo(i-1,j)*mo(i,j)-3.0_rsh*v(i-1,j)+3.0_rsh*v(i  ,j)-v(i+1,j))   &
                     +um/6.0_rsh*(1-um*um)*(v(i-1,j)-3.0_rsh*v(i  ,j)+3.0_rsh*v(i+1,j)-v(i+2,j)*me(i+1,j)*me(i,j)) ) &
                    /dx_v(i,j)
           ELSE
   ! upwind scheme
             vdvsdy  = vp*(v(i,j)-v(i,j-1))/dy(i,j)+vm*(v(i,j+1)-v(i,j))/dy(i,j+1)
             udvsdx  = up*(v(i,j)-v(i-1,j))/dx_f(i-1,j)+um*(v(i+1,j)-v(i,j))/dx_f(i,j)
           END IF
#elif key_dyn_adv_quick
   ! QUICK scheme
           vc=v(i,j)
           vp=0.5_rsh*(vc+ABS(vc))
           vm=0.5_rsh*(vc-ABS(vc))
           up=0.5_rsh*(u_vtmp(i,j)+ABS(u_vtmp(i,j)))
           um=0.5_rsh*(u_vtmp(i,j)-ABS(u_vtmp(i,j)))
           IF (i >= imin+3.AND.i <= imax-2.AND.j <= jmax-2.AND.j >= jmin+2) THEN
             vdvsdy=(vp*v(i,j-2)*why(i,j-1)+(-7.0_rsh*vp-3.0_rsh*vm)*v(i,j-1)+         &
                     (3.0_rsh*vp-3.0_rsh*vm)*v(i,j)+(3.0_rsh*vp+7.0_rsh*vm)*v(i,j+1)-vm*v(i,j+2)*why(i,j+1)) &
                    /(8.0_rsh*dy_v(i,j))
             udvsdx=(up*v(i-2,j)*mo(i-1,j)*mo(i,j)+(-7.0_rsh*up-3.0_rsh*um)*v(i-1,j)+  &
                     (3.0_rsh*up-3.0_rsh*um)*v(i,j)+(3.0_rsh*up+7.0_rsh*um)*v(i+1,j)-um*v(i+2,j)*me(i+1,j)*me(i,j)) &
                    /(8.0_rsh*dx_v(i,j))
           ELSE
   ! upwind scheme
             vdvsdy  = vp*(v(i,j)-v(i,j-1))/dy(i,j)+vm*(v(i,j+1)-v(i,j))/dy(i,j+1)
             udvsdx  = up*(v(i,j)-v(i-1,j))/dx_f(i-1,j)+um*(v(i+1,j)-v(i,j))/dx_f(i,j)
           END IF
#else
   ! default : centered scheme
           vdvsdy=v(i,j)*(v(i,j+1)-v(i,j-1))/(dy(i,j+1)+dy(i,j))
           udvsdx=0.5_rsh*u_vtmp(i,j)*                                       &
                  (  mo(i,j)*(v(i  ,j)-v(i-1,j))/dx_f(i-1,j)                 &
                   + me(i,j)*(v(i+1,j)-v(i  ,j))/dx_f(i  ,j) )
#endif
#endif
#endif

           vstar(i,j)=v(i,j)-halfdt*(                                        &
#if ! defined key_dyn_adv_none
#if ! defined key_dyn_adv_linear
                    ! non-linear advection : explicit
                    vdvsdy + udvsdx + tfivsr(i,j)*u_vtmp(i,j)**2             &
#endif
#endif
                    ! bottom stress : explicit
                   +ft*cfroe*v(i,j)                                          &
                    ! viscosity (along y) : explicit
                   -( vis2d(i,j+1)*(v(i,j+1)-v(i,j  ))/dy(i,j+1)             &
                     -vis2d(i,j  )*(v(i,j  )-v(i,j-1))/dy(i,j  ) )/dy_v(i,j) &
                    ! sea surface height gradient (along y) : explicit
                   +grav*(xe(i,j+1)-xe(i,j))/dy_v(i,j)                       &
                    ! coriolis term : explicit
                   +fv(i,j)*u_vtmp(i,j)                                      &
                    ! viscosity (along x) : explicit
                   +( vis2d(i-1,j)*mo(i,j)*(v(i  ,j)-v(i-1,j))/dx_f(i-1,j)   &
                     -vis2d(i  ,j)*me(i,j)*(v(i+1,j)-v(i  ,j))/dx_f(i  ,j) ) &
                     /dx_v(i,j) )

        ! wind stress and atmospherical pressure gradient
           vstar(i,j)=vstar(i,j)                                     &
                  ! wind stress in y-direction : explicit
                     +halfdt*tosy(i,j)/MAX(hty(i,j),hminfrot)/rhoref &
                  ! atmospherical pressure gradient : explicit
                     -halfdt*(pa(i,j+1)-pa(i,j))/(dy_v(i,j)*rhoref)

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

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

           vstar(i,j) = vstar(i,j)/(1.0_rlg + ft*cfroi*halfdt)

#ifdef key_waveforcing
           IF (l_wave_brgm_phillips) THEN
             vstar(i,j) = vstar(i,j) - halfdt *( temp_radv(i,j) / rhoref + v  (i,j) * (Qy(i,j+1)-Qy(i,j))/dy_v(i,j) &
                                               +  2.0_rsh*Qy_v(i,j) * (v(i,j+1)-v(i,j-1))/(dy(i,j)+dy(i,j+1)) )     &
                                              / MAX(hty(i,j),hminfrot)
             IF ( htx(i-1,j)+htx(i-1,j+1) >= 0.0_rsh )  &
               vstar(i,j) = vstar(i,j) - halfdt * ( Qy_v(i,j) * (u_vtmp(i,j )-u_vtmp(i-1,j)) &
                                                  + Qx_v(i,j) * (v (i,j  )-v   (i-1,j))      &
                                                  + u_vtmp(i,j) * (Qy_v(i,j)-Qy_v(i-1,j))  ) &
                                                / (2.0_rsh*dx_u(i-1,j)*MAX(hty(i,j),hminfrot))
             IF ( htx(i  ,j)+htx(i  ,j+1) >= 0.0_rsh ) &
               vstar(i,j) = vstar(i,j) - halfdt * (  Qy_v(i,j) * (u_vtmp (i+1,j)-u_vtmp (i,j)) &
                                                  +  Qx_v(i,j) * (v   (i+1,j)-v   (i,j))       &
                                                  +  u_vtmp(i,j) * (Qy_v(i+1,j)-Qy_v(i,j)) )   &
                                                / (2.0_rsh*dx_u(i,j)*MAX(hty(i,j),hminfrot))

           ELSE IF (l_wave_brgm_smith) THEN
             vstar(i,j) = vstar(i,j) - halfdt * (JSMITH(i,j+1)-JSMITH(i,j)) /dy_v(i,j)
             vstar(i,j) = vstar(i,j) &
                          + halfdt * (k_y(i,j)+k_y(i,j+1))/2.0_rsh * (Dsmith(i,j)+Dsmith(i,j+1))/2.0_rsh &
                                   / MAX(hty(i,j),hminfrot)/rhoref
             vstar(i,j) = vstar(i,j) + halfdt * Qx_v(i,j) * (u_vtmp(i,j+1)-u_vtmp(i,j-1))/(dy(i,j)+dy(i,j+1)) &
                                              / MAX(hty(i,j),hminfrot)
             IF ((htx(i-1,j)+htx(i-1,j+1)) >= 0.0_rsh) &
               vstar(i,j) = vstar(i,j) - halfdt * Qx_v(i,j) * (v  (i,j)-v  (i-1,j))/(2.0_rsh*dx_u(i-1,j)) &
                                                / MAX(hty(i,j),hminfrot)
             IF ((htx(i,j)+htx(i,j+1)) >= 0.0_rsh) &
               vstar(i,j) = vstar(i,j) - halfdt * Qx_v(i,j) * (v  (i+1,j)-v  (i,j))/(2.0_rsh*dx_u(i,j )) &
                                                / MAX(hty(i,j),hminfrot)
           ELSE
             vstar(i,j) = vstar(i,j) - halfdt * temp_radv(i,j) / MAX(hty(i,j),hminfrot)/rhoref
             vstar(i,j) = vstar(i,j) + ft*Qy_u(i,j)/MAX(hty(i,j),hminfrot)
           END IF
#endif 

         END IF

       ENDIF   ! ends test on j
     END DO
   END DO
!$OMP END DO

   PRINT_DBG*, 'exit DYN2D_UVSTAR routine'

  END SUBROUTINE dyn2d_uvstar
