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

  SUBROUTINE dyn3dx_epg(xec,xeuv)

   !&E---------------------------------------------------------------------
   !&E                 ***  dyn3dx_epg  ***
   !&E
   !&E ** Purpose : estimate the trend term relative to the sea surface
   !&E              pressure gradient for the zonal momentum forcing
   !&E              (3D equation)
   !&E
   !&E ** Description : semi-implicit
   !&E
   !&E ** Called by : dyn3dxy.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays : igdu,yz_exp,dx_u,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---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d, ONLY : ig,id,dx_u
   USE comcolig2d,ONLY : igdu
   USE comcolig3d,ONLY : yz_exp
#ifdef key_diagdyn
   USE diagdyn,   ONLY : mask_diagx,diagx_epg
#endif

   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

   PRINT_DBG*, 'enter DYN3DX_EPG routine'

!$OMP DO schedule(runtime) private(i,j)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)),MIN0(limax,id(j)-1)
       IF (igdu(i,j) == 1) THEN
         yz_exp(:,i,j)=yz_exp(:,i,j)-grav/dx_u(i,j)*(      &
                                           alphap*( xec(i+1,j)-xec(i,j) )   &
                                          +alphan*(xeuv(i+1,j)-xeuv(i,j)) )

#ifdef key_diagdyn
         IF (mask_diagx(i,j)) THEN
           diagx_epg(:,i,j)=diagx_epg(:,i,j)-grav/dx_u(i,j)*(      &
                                             alphap*( xec(i+1,j)-xec(i,j) )   &
                                            +alphan*(xeuv(i+1,j)-xeuv(i,j)) )
         END IF
#endif

       END IF
     END DO
   END DO
!$OMP END do

   PRINT_DBG*, 'exit DYN3DX_EPG routine'

  END SUBROUTINE dyn3dx_epg
