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

  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---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d, ONLY : hex,hey,hx,hy,goulet,govlet,dhx,dhy
#if defined key_siggen || defined key_gencoord
   USE comsiggen, ONLY : hzex,hzey
   USE toolgeom,  ONLY : f_hzu_x,f_hzu_y
#endif
   USE_MPI 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_DBG*, 'enter INIT_HEIGHT routine'




#ifdef key_MPI_2D
!$OMP DO schedule(runtime)
#else
!$OMP DO schedule(dynamic,3)
#endif

   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
#ifdef key_dyn_adv_linear
         hex(i,j)=MAX(0.0_rsh,hx(i,j))*goulet(i,j)
#if defined key_siggen || defined key_gencoord
         DO k=1,kmax
           hzex(k,i,j)=f_hzu_x(hx(i,j),dhx(i,j),k,i,j)*goulet(i,j)
         END DO
#endif
#else
         hex(i,j)=(htx+dhx(i,j))*goulet(i,j)
#if defined key_siggen || defined key_gencoord
         IF ( hx(i,j)+0.5_rsh*(xe(i,j)+xe(i+1,j)) < 0.0_rsh ) THEN
           DO k=1,kmax
             hzex(k,i,j)=f_hzu_x(hx(i,j),MAX(xe(i,j),xe(i+1,j))+dhx(i,j),k,i,j)*goulet(i,j)
           END DO
         ELSE
           DO k=1,kmax
             hzex(k,i,j)=f_hzu_x(hx(i,j),0.5_rsh*(xe(i,j)+xe(i+1,j))+dhx(i,j),k,i,j)*goulet(i,j)
           END DO
         END IF
#endif
#endif
       ELSE
         hex(i,j)=0.0_rsh
#if defined key_siggen || defined key_gencoord
         hzex(:,i,j)=0.0_rsh
#endif
       END IF


       IF (hty > 0.0_rsh) THEN
#ifdef key_dyn_adv_linear
         hey(i,j)=MAX(0.0_rsh,hy(i,j))*govlet(i,j)
#if defined key_siggen || defined key_gencoord
         DO k=1,kmax
           hzey(k,i,j)=f_hzu_y(hy(i,j),dhy(i,j),k,i,j)*govlet(i,j)
         END DO
#endif
#else
         hey(i,j)=(hty+dhy(i,j))*govlet(i,j)
#if defined key_siggen || defined key_gencoord
         IF ( hy(i,j)+0.5_rsh*(xe(i,j)+xe(i,j+1)) < 0.0_rsh ) THEN
           DO k=1,kmax
             hzey(k,i,j)=f_hzu_y(hy(i,j),MAX(xe(i,j),xe(i,j+1))+dhy(i,j),k,i,j)*govlet(i,j)
           END DO
         ELSE
           DO k=1,kmax
             hzey(k,i,j)=f_hzu_y(hy(i,j),0.5_rsh*(xe(i,j)+xe(i,j+1))+dhy(i,j),k,i,j)*govlet(i,j)
           END DO
         END IF
#endif
#endif
       ELSE
         hey(i,j)=0.0_rsh
#if defined key_siggen || defined key_gencoord
         hzey(:,i,j)=0.0_rsh
#endif
       END IF

     END DO
   END DO
!$OMP END DO

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

   PRINT_DBG*, 'exit INIT_HEIGHT routine'
 
  END SUBROUTINE init_height
