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

  SUBROUTINE init_height(xe)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE init_height  ***
!&E
!&E ** Purpose : Evaluation of water column height and layer depths
!&E              at u,v locations
!&E
!&E ** Description :
!&E
!&E ** Called by : main.F90,initfromfile.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : hx,hy,xe
!&E
!&E ** Modified variables : hex,hey,hzex,hzey
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2008-07  (V. Garnier)  Original code from dyn3d_height_mask
!&E       !  2011-03  (F. Dumas) remove case 2 and 3 on tidal flats
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2014-11  (V. Garnier) Estimate of hzex/y instead of hzutx/y to estimate barotropic velocities
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d, ONLY : hex,hey,hx,hy,goulet,govlet,dhx,dhy

   USE toolmpi, ONLY : ex_i_rsh,ex_j_rsh,ex_i_int,ex_j_int

   IMPLICIT NONE

!! * Arguments
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: xe

!! * Local declarations
   INTEGER                                                     :: i,j,k,i1,i2,j1,j2
   REAL(KIND=rsh)                                              :: htx,hty

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

   !PRINT*, 'enter INIT_HEIGHT routine'





!$OMP DO schedule(runtime)


   DO j=ljmin,ljmaxm1
     DO i=limin,limaxm1


       htx=hx(i,j)+0.5_rsh*(xe(i,j)+xe(i+1,j))
       hty=hy(i,j)+0.5_rsh*(xe(i,j)+xe(i,j+1))
       IF ( htx<=0.0_rsh .AND. ( hx(i,j)+xe(i,j) > 0.0_rsh .OR.  hx(i,j)+xe(i+1,j) > 0.0_rsh) ) &
         htx=hx(i,j)+MAX(xe(i,j),xe(i+1,j))
       IF ( hty<=0.0_rsh .AND. ( hy(i,j)+xe(i,j) > 0.0_rsh .OR.  hy(i,j)+xe(i,j+1) > 0.0_rsh )) &
         hty=hy(i,j)+MAX(xe(i,j),xe(i,j+1))


       IF (htx > 0.0_rsh) THEN

         hex(i,j)=(htx+dhx(i,j))*goulet(i,j)


       ELSE
         hex(i,j)=0.0_rsh

       END IF


       IF (hty > 0.0_rsh) THEN

         hey(i,j)=(hty+dhy(i,j))*govlet(i,j)


       ELSE
         hey(i,j)=0.0_rsh

       END IF

     END DO
   END DO
!$OMP END DO

!!$OMP barrier
!!$OMP master
!   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljmin,ljmax,hex(liminm2:limaxp2,ljmin:ljmax))
!   CALL ex_j_rsh(-2,2,1,limin,limax,ljminm2,ljmaxp2,hey(limin:limax,ljminm2:ljmaxp2))
!#if defined key_siggen || defined key_gencoord
!   CALL ex_i_rsh(-1,0,kmax,liminm1,limax,ljmin,ljmax,hzex(:,liminm1:limax,ljmin:ljmax))
!   CALL ex_j_rsh(-1,0,kmax,limin,limax,ljminm1,ljmax,hzey(:,limin:limax,ljminm1:ljmax))
!#endif
!!$OMP end master
!!$OMP barrier
!!$OMP flush(hex,hey,whx,why,igdu,igdv,mn,ms,mo,me)
!#if defined key_siggen || defined key_gencoord
!!$OMP flush(hzex,hzey)
!#endif

   !PRINT*, 'exit INIT_HEIGHT routine'
 
  END SUBROUTINE init_height
