!!======================================================================

  SUBROUTINE dyn3dy_epg(xec,xeuv)

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dy_epg  ***
!&E
!&E ** Purpose : estimate the trend term relative to the sea surface
!&E              pressure gradient for the meridional momentum forcing
!&E              (3D equation)
!&E
!&E ** Description : semi-implicit
!&E
!&E ** Called by : dyn3dyx.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : igdv,yz_exp,dy_v,xec,xeuv
!&E
!&E ** Modified variables : yz_exp
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  1987-02  (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       !  2009-04  (A. Mano) MPI2D
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters
   USE comvars2d,  ONLY : dy_v
   USE comcolig2d, ONLY : igdv
   USE comcolig3d, ONLY : yz_exp


   IMPLICIT NONE

!! * Arguments
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: xeuv
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: xec

!! * Local declarations
   INTEGER                                       :: i,j

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

!$OMP DO schedule(runtime),private(i,j)
   DO j=ljmin,ljmaxm1
     DO i=liminp2,limaxm1
       IF(igdv(i,j) == 1) THEN
         yz_exp(:,i,j)=yz_exp(:,i,j)-grav/dy_v(i,j)*(    &
                                     alphap*( xec(i,j+1)-xec(i,j) )   &
                                    +alphan*(xeuv(i,j+1)-xeuv(i,j)) )



       END IF
     END DO
   END DO
!$OMP end do

  END SUBROUTINE dyn3dy_epg
