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

  SUBROUTINE init_xe(xeu)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE init_xe  ***
   !&E
   !&E ** Purpose : Initialize water levels.
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : main
   !&E
   !&E ** External calls : obc2d_rank, tide_astroarg, tide_harmcompo
   !&E
   !&E ** Used ij-arrays : h0
   !&E
   !&E ** Modified variables :  xeu, xeadv
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2004-10-04
   !&E       !  1987-02  (J.-C. Salomon, P. Lazure)  Original code
   !&E       !  2009-04  (V. Garnier) separation of OBC origin (OGCM/tide)
   !&E                                use of obc_tide_ssh*
   !&E       !  2009-10  (F. Dumas, E. Romieu) add obc with and without wind to tidal signal
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2012-05  (V. Garnier, L. Pineau) SSH to mslshift at imin and jmin if western or southern open boundary
   !&E       !  2016-12  (A. Fiandrino,B.Thouvenin),V. Garnier) integration of key_lagoon
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d, ONLY : h0,h0_g,ig,id,jb,jh,hminim,mslshift,l_initfromfile,l_tide_surge
   USE comvars2d, ONLY : fwet,fwetp,fwetmin,dhxy
#ifdef key_lagoon
   USE comvars2d, ONLY : igrau,jgrau,type_grau
#ifdef key_MPI_2D
   USE comvars2d, ONLY : type_mpi_rsh
#endif
#endif
   USE comaree,   ONLY : l_tide_harmcompo
   USE comfluxdo, ONLY : xeadv
   USE obcogcm,   ONLY : obc_ogcm
   USE obccombine,ONLY : l_obc_ogcm_rssh,obc_combine,obc_tide_sshs,obc_tide_sshn, &
                         obc_tide_sshe,obc_tide_sshw,l_obc_tide
   USE obcapplyssh
   USE_MPI toolmpi, ONLY : exchange_rsh,ex_i_rsh,ex_j_rsh,MASTER

   IMPLICIT NONE
   
#ifdef key_MPI_2D
   INCLUDE 'mpif.h'
#endif

   !! * Arguments
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(out) :: xeu

   !! * Local declarations
   REAL(kind=rsh),DIMENSION(imin:imax,jmin:jmax) :: som
   REAL(kind=rsh),DIMENSION(imin:imax,jmin:jmax) :: xe
   REAL(kind=rsh)                                :: d2
   INTEGER                                       :: i,j,ii,jj,i1,j1
#ifdef key_lagoon
   REAL(kind=rsh)      :: xe_init,value
   INTEGER             :: IERR_MPI
#endif

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

   PRINT_DBG*, 'enter INIT_XE routine'

   som(:,:)=0.0_rsh
   xe(:,:)=0.0_rsh
   xeu(:,:)=valmanq

   ! Intermediate calculations
   ! -------------------------

   IF (l_tide_surge) THEN
   ! add obc with and without wind to tidal signal
     CALL tide_harmcompo
     CALL obc2d_rank 
   ELSE
     IF (l_obc_tide) THEN
       IF (l_tide_harmcompo) THEN
#if ! defined key_lagoon
         CALL tide_astroarg
#endif
         CALL tide_harmcompo
       ELSE
#if ! defined key_offline
         CALL obc2d_rank
#endif
       END IF
     END IF
   END IF

   IF (l_obc_tide ) THEN
    IF_MPI(MASTER) print*,"INIT XE WITH OBC"
#ifdef key_lagoon
   ! for key_lagoon : l_obc_tide=.true. always
   ! CLIM donnee par le niveau en mer au grau numero 1
   !  pas de conditions limites en x interpolables
   !  une haut eau / zero ngf.
   ! le niveau moyen est ramene a zero et non plus a xemoy_m
   ! dans limah
   xe_init=0.0_rsh
   i1=igrau(1)
   j1=jgrau(1)
   IF_MPI (liminm2<=i1 .AND. i1<=limaxp2 .AND. ljminm2<=j1 .AND. j1<=ljmaxp2) THEN
    IF (type_grau(1)==1) xe_init=obc_tide_sshn(i1)
    IF(type_grau(1)==2) xe_init=obc_tide_sshe(j1)
    IF(type_grau(1)==3) xe_init=obc_tide_sshs(i1)
    IF (type_grau(1)==4) xe_init=obc_tide_sshw(j1)
   ENDIF_MPI 

#ifdef key_MPI_2D
   CALL MPI_ALLREDUCE(xe_init,value,1,type_mpi_rsh,MPI_SUM,MPI_COMM_WORLD,IERR_MPI)
   xe_init=value
#endif
   xeu(:,:)=xe_init
   WHERE(h0(liminm2:limaxp2,ljminm2:ljmaxp2)== -valmanq) xeu(:,:)=valmanq

!cval$OMP do schedule(runtime)
     DO j=ljmin,ljmax
       DO i=limin,limax
         IF (dhxy(i,j) == 0.0_rsh) THEN
           fwet(i,j)=1.0_rsh
           fwetp(i,j)=1.0_rsh
         ELSE
           fwet(i,j)=MIN(MAX(fwetmin,(xeu(i,j)+h0(i,j))/dhxy(i,j)),1.0_rsh)
           fwetp(i,j)=MIN(MAX(fwetmin,(xeu(i,j)+h0(i,j))/dhxy(i,j)),1.0_rsh)
         ENDIF
       END DO
     END DO
!cval$OMP end do
!cval$OMP end parallel

!  end key_lagoon
#else   
    
!cval$OMP parallel default(none) private(i,j,d2,jj,i1,ii,j1) &
!cval$OMP& shared(ljmin,ljmax,h0,h0_g,ig,id,jb,jh,hminim,som,xeu,xeadv,mslshift)

   ! initialize the open boundaries
     IF_MPI (limin==imin) &
     !xe(imin+1,ljminp1:ljmax)=obc_tide_sshw(ljminp1:ljmax)
     xe(imin+31,ljminp1:ljmax)=obc_tide_sshw(ljminp1:ljmax) ! iman 4mai2015 
     IF_MPI (limax==imax) &
     xe(imax  ,ljminp1:ljmax)=obc_tide_sshe(ljminp1:ljmax)
     IF_MPI (ljmin==jmin) &
     xe(liminp1:limax,jmin+1)=obc_tide_sshs(liminp1:limax)
     IF_MPI (ljmax==jmax) &
     xe(liminp1:limax,jmax  )=obc_tide_sshn(liminp1:limax)

OMPMPI barrier
OMPMPI master
     CALL_MPI exchange_rsh(xe,1)
OMPMPI end master
OMPMPI barrier
OMPMPI flush(xe)

     DO jj=jmin+2,jmax-1,6

   ! from the western open boundary

       i1=ig(jj)

       IF (h0_g(i1,jj) > hminim) THEN
!cval$OMP do schedule(runtime)
         DO j=jmin+2,jmax-1
           DO i=ig(j)+1,id(j)-1
             IF (h0_g(i,j) == -valmanq) THEN
               som(i,j)=0.0_rsh
               xe(i,j)=valmanq
             ELSE
               IF (i1==i.AND.jj==j) THEN
                 xe(i,j)=xe(i,j)+xe(i1,jj)
                 som(i,j)=som(i,j)+1.0_rsh
               ELSE
                 d2=(i-i1)*(i-i1)+(j-jj)*(j-jj)
                 xe(i,j)=xe(i,j)+xe(i1,jj)/d2
                 som(i,j)=som(i,j)+1.0_rsh/d2
               END IF
             END IF
           END DO
         END DO
!cval$OMP end do
       END IF

   ! from the eastern open boundary

       i1=id(jj)
       IF (h0_g(i1,jj) > hminim) THEN
!cval$OMP do schedule(runtime)
         DO j=jmin+2,jmax-1
           DO i=ig(j)+1,id(j)-1
             IF (h0_g(i,j) == -valmanq) THEN
               som(i,j)=0.0_rsh
               xe(i,j)=valmanq
             ELSE
               IF (i1==i.AND.jj==j) THEN
                 xe(i,j)=xe(i,j)+xe(i1,jj)
                 som(i,j)=som(i,j)+1.0_rsh
               ELSE
                 d2=(i-i1)*(i-i1)+(j-jj)*(j-jj)
                 xe(i,j)=xe(i,j)+xe(i1,jj)/d2
                 som(i,j)=som(i,j)+1.0_rsh/d2
               END IF
             END IF
           END DO
         END DO
!cval$OMP end do
       END IF

     END DO

     DO ii=imin+2,imax-1,6

   ! from the southern open boundary

       j1=jb(ii)
       IF (h0_g(ii,j1) > hminim) THEN
!cval$OMP do schedule(runtime)
         DO j=jmin,jmax
           DO i=imin+2,imax-1
             IF ((j.GE.jb(i)+1).AND.(j.LE.jh(i)-1)) THEN
!           do j=jb(i)+1,jh(i)-1
               IF (h0_g(i,j) == -valmanq) THEN
                 som(i,j)=0.0_rsh
                 xe(i,j)=valmanq
               ELSE
                 IF (ii==i.AND.j==j1) THEN
                   xe(i,j)=xe(i,j)+xe(ii,j1)
                   som(i,j)=som(i,j)+1.0_rsh
                 ELSE
                   d2=(i-ii)*(i-ii)+(j-j1)*(j-j1)
                   xe(i,j)=xe(i,j)+xe(ii,j1)/d2
                   som(i,j)=som(i,j)+1.0_rsh/d2
                 END IF
               END IF
             END IF
           END DO
         END DO
!cval$OMP end do
       END IF

   ! from the northern open boundary

       j1=jh(ii)
       IF (h0_g(ii,j1) > hminim) THEN
!cval$OMP do schedule(runtime)
         DO j=jmin,jmax
           DO i=imin+2,imax-1
!           do j=jb(i)+1,jh(i)-1
             IF ((j.GE.jb(i)+1).AND.(j.LE.jh(i)-1)) THEN
               IF (h0_g(i,j) == -valmanq) THEN
                 som(i,j)=0.0_rsh
                 xe(i,j)=valmanq
               ELSE
                 IF (ii==i.AND.j==j1) THEN
                   xe(i,j)=xe(i,j)+xe(ii,j1)
                   som(i,j)=som(i,j)+1.0_rsh
                 ELSE
                   d2=(i-ii)*(i-ii)+(j-j1)*(j-j1)
                   xe(i,j)=xe(i,j)+xe(ii,j1)/d2
                   som(i,j)=som(i,j)+1.0_rsh/d2
                 END IF
               END IF
             END IF
           END DO
         END DO
!cval$OMP end do
       END IF

     END DO


   ! initialization of xeu and xeadv
   ! -------------------------------

   ! initialize the open boundaries

     IF_MPI (limin==imin) &
     !xeu(imin+1,ljmin:ljmaxp1)=obc_tide_sshw(ljmin:ljmaxp1)
     xeu(imin+31,ljminp1:ljmaxp1)=obc_tide_sshw(ljminp1:ljmaxp1) ! iman 4mai2015 
     IF_MPI (limax==imax) &
     xeu(imax  ,ljmin:ljmaxp1)=obc_tide_sshe(ljmin:ljmaxp1)
     IF_MPI (ljmin==jmin) &
     xeu(limin:limaxp1,jmin+1)=obc_tide_sshs(limin:limaxp1)
     IF_MPI (ljmax==jmax) &
     xeu(limin:limaxp1,jmax  )=obc_tide_sshn(limin:limaxp1)

! do i=imin+2,imax-1
!  do j=jb(i)+1,jh(i)-1
!   if(som(i,j) /= 0.0_rsh) then
!    xeu(i,j)=max(xe(i,j)/som(i,j),-h0_g(i,j))
!   else
!    xeu(i,j)=valmanq
!   endif
!~  enddo
! enddo
!

!cval$OMP do schedule(runtime)
     DO j=ljmin,ljmax
       DO i=limin,limax
         IF ((i.GE.imin+2.AND.i.LE.imax-1).AND.(j.GE.jb(i)+1).AND.(j.LE.jh(i)-1)) THEN
           IF (som(i,j) /= 0.0_rsh) THEN
             xeu(i,j)=MAX(xe(i,j)/som(i,j),-h0(i,j))
           ELSE
             xeu(i,j)=valmanq
           END IF
         END IF
         !IF (xeu(i,j) == 0.0_rsh) xeu(i,j)=valmanq
         IF ((xeu(i,j) == valmanq) .AND. (h0(i,j)/= -valmanq)) xeu(i,j)=mslshift
         IF (j.EQ.jmin .AND. h0(i,j)== -valmanq) xeu(i,jmin)=valmanq
         IF (i.EQ.imin .AND. h0(i,j)== -valmanq) xeu(imin,j)=valmanq
         IF (dhxy(i,j) == 0.0_rsh) THEN
            fwet(i,j)=1.0_rsh
           fwetp(i,j)=1.0_rsh
         ELSE
            fwet(i,j)=MIN(MAX(fwetmin,(xeu(i,j)+h0(i,j))/dhxy(i,j)),1.0_rsh)
           fwetp(i,j)=MIN(MAX(fwetmin,(xeu(i,j)+h0(i,j))/dhxy(i,j)),1.0_rsh)
         ENDIF
       END DO
     END DO
!cval$OMP end do
!cval$OMP end parallel

    DO j=ljmin,ljmax
       DO i=limin,limax
       if (i==31 .and. j==100) print*,'init xe',xeu(i,j),xeu(i+1,j)
        
     END DO
     END DO

     CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeu(liminm2:limaxp2,ljminm2:ljmaxp2))
     CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,xeu(liminm2:limaxp2,ljminm2:ljmaxp2))
#endif

   ELSE

     WHERE (h0(liminm2:limaxp2,ljminm2:ljmaxp2) /= -valmanq) xeu(liminm2:limaxp2,ljminm2:ljmaxp2)=0.0_rsh

   END IF

   xeadv(limin:limax,ljmin:ljmax)=xeu(limin:limax,ljmin:ljmax)

   PRINT_DBG*, 'exit INITXE routine'

  END SUBROUTINE init_xe
