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

  SUBROUTINE obc_init_height(xe)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE obc_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 : htx,hzutx,hty,hzuty
   !&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
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d, ONLY : htx,hty,hx,hy
#if defined key_siggen || defined key_gencoord
   USE comsiggen, ONLY : hzutx,hzuty
   USE toolgeom,  ONLY : f_hzu_x,f_hzu_y
#endif

   IMPLICIT NONE

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

   !! * Local declarations
   INTEGER                                                     :: i,j,k

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

#if defined key_siggen || defined key_gencoord

#ifdef key_MPI_2D
!$OMP DO schedule(runtime)
#else
!$OMP DO schedule(static)
#endif
   DO j=ljminp1,ljmaxm1

   ! Western boundary
     DO i=liminp1,MIN0(obc_width,limaxm1)
       htx(i,j)=hx(i,j)+0.5_rsh*(xe(i,j)+xe(i+1,j))
       hty(i,j)=hy(i,j)+0.5_rsh*(xe(i,j)+xe(i,j+1))
       DO k=1,kmax
   ! dz/ds at location hx (i+1/2) and depth u or w
         hzutx(k,i,j)=f_hzu_x(hx(i,j),0.5_rsh*(xe(i,j)+xe(i+1,j)),k,i,j)
   ! dz/ds at location hy (j+1/2) and depth of v or w
         hzuty(k,i,j)=f_hzu_y(hy(i,j),0.5_rsh*(xe(i,j)+xe(i,j+1)),k,i,j)
       END DO
   ! dz/ds at location hx (i+1/2) and depth u or w
       IF ( htx(i,j)<=0.0_rsh .AND. ( hx(i,j)+xe(i,j) > 0.0_rsh .OR. hx(i,j)+xe(i+1,j) > 0.0_rsh) ) THEN
         DO k=1,kmax
           hzutx(k,i,j)=f_hzu_x(hx(i,j),MAX(xe(i,j),xe(i+1,j)),k,i,j)
         END DO
       END IF
   ! dz/ds at location hy (j+1/2) and depth of v or w
       IF ( hty(i,j)<=0.0_rsh .AND. ( hy(i,j)+xe(i,j) > 0.0_rsh .OR. hy(i,j)+xe(i,j+1) > 0.0_rsh )) THEN
         DO k=1,kmax
           hzuty(k,i,j)=f_hzu_y(hy(i,j),MAX(xe(i,j),xe(i,j+1)),k,i,j)
         END DO
       END IF
     END DO

   ! Western boundary following
     IF_MPI (limin==imin) THEN
       htx(imin,j)=hx(imin,j)+         xe(imin+1,j)
       hty(imin,j)=hy(imin,j)+0.5_rsh*(xe(imin+1,j)+xe(imin+1,j+1))
       DO k=1,kmax
   ! dz/ds at location hx (i+1/2) and depth u or w
         hzutx(k,imin,j)=f_hzu_x(hx(imin,j),         xe(imin+1,j)                ,k,imin,j)
   ! dz/ds at location hy (j+1/2) and depth of v or w
         hzuty(k,imin,j)=f_hzu_y(hy(imin,j),0.5_rsh*(xe(imin+1,j)+xe(imin+1,j+1)),k,imin,j)
       END DO
       IF ( hty(imin,j)<=0.0_rsh .AND. ( hy(imin,j)+xe(imin+1,j) > 0.0_rsh .OR. hy(imin,j)+xe(imin+1,j+1) > 0.0_rsh )) THEN
         DO k=1,kmax
           hzuty(k,imin,j)=f_hzu_y(hy(imin,j),MAX(xe(imin+1,j),xe(imin+1,j+1)),k,imin,j)
         END DO
       END IF
     ENDIF_MPI

   ! Eastern boundary
     DO i=MAX0(imax-obc_width+1,limin),limaxm1
       htx(i,j)=hx(i,j)+0.5_rsh*(xe(i,j)+xe(i+1,j))
       hty(i,j)=hy(i,j)+0.5_rsh*(xe(i,j)+xe(i,j+1))
       DO k=1,kmax
   ! dz/ds at location hx (i+1/2) and depth u or w
         hzutx(k,i,j)=f_hzu_x(hx(i,j),0.5_rsh*(xe(i,j)+xe(i+1,j)),k,i,j)
   ! dz/ds at location hy (j+1/2) and depth of v or w
         hzuty(k,i,j)=f_hzu_y(hy(i,j),0.5_rsh*(xe(i,j)+xe(i,j+1)),k,i,j)
       END DO
   ! dz/ds at location hx (i+1/2) and depth u or w
       IF ( htx(i,j)<=0.0_rsh .AND. ( hx(i,j)+xe(i,j) > 0.0_rsh .OR. hx(i,j)+xe(i+1,j) > 0.0_rsh) ) THEN
         DO k=1,kmax
           hzutx(k,i,j)=f_hzu_x(hx(i,j),MAX(xe(i,j),xe(i+1,j)),k,i,j)
         END DO
       END IF
   ! dz/ds at location hy (j+1/2) and depth of v or w
       IF ( hty(i,j)<=0.0_rsh .AND. ( hy(i,j)+xe(i,j) > 0.0_rsh .OR. hy(i,j)+xe(i,j+1) > 0.0_rsh )) THEN
         DO k=1,kmax
           hzuty(k,i,j)=f_hzu_y(hy(i,j),MAX(xe(i,j),xe(i,j+1)),k,i,j)
         END DO
       END IF
     END DO

   ! Eastern boundary following
     IF_MPI (limax==imax) THEN
       htx(imax,j)=hx(imax,j)+         xe(imax,j)
       hty(imax,j)=hy(imax,j)+0.5_rsh*(xe(imax,j)+xe(imax,j+1))
       DO k=1,kmax
   ! dz/ds at location hx (i+1/2) and depth u or w
         hzutx(k,imax,j)=f_hzu_x(hx(imax,j),         xe(imax,j)              ,k,imax,j)
   ! dz/ds at location hy (j+1/2) and depth of v or w
         hzuty(k,imax,j)=f_hzu_y(hy(imax,j),0.5_rsh*(xe(imax,j)+xe(imax,j+1)),k,imax,j)
       END DO
       IF ( hty(imax,j)<=0.0_rsh .AND. ( hy(imax,j)+xe(imax,j) > 0.0_rsh .OR. hy(imax,j)+xe(imax,j+1) > 0.0_rsh )) THEN
         DO k=1,kmax
           hzuty(k,imax,j)=f_hzu_y(hy(imax,j),MAX(xe(imax,j),xe(imax,j+1)),k,imax,j)
         END DO
       END IF
     ENDIF_MPI

   END DO
!$OMP END DO

#ifdef key_MPI_2D
!$OMP DO schedule(runtime)
#else
!$OMP DO schedule(static)
#endif
   DO i=liminp1,limaxm1

   ! Southern boundary
     DO j=ljminp1,MIN0(obc_width,ljmaxm1)
       htx(i,j)=hx(i,j)+0.5_rsh*(xe(i,j)+xe(i+1,j))
       hty(i,j)=hy(i,j)+0.5_rsh*(xe(i,j)+xe(i,j+1))
       DO k=1,kmax
   ! dz/ds at location hx (i+1/2) and depth u or w
         hzutx(k,i,j)=f_hzu_x(hx(i,j),0.5_rsh*(xe(i,j)+xe(i+1,j)),k,i,j)
   ! dz/ds at location hy (j+1/2) and depth of v or w
         hzuty(k,i,j)=f_hzu_y(hy(i,j),0.5_rsh*(xe(i,j)+xe(i,j+1)),k,i,j)
       END DO
   ! dz/ds at location hx (i+1/2) and depth u or w
       IF ( htx(i,j)<=0.0_rsh .AND. ( hx(i,j)+xe(i,j) > 0.0_rsh .OR. hx(i,j)+xe(i+1,j) > 0.0_rsh) ) THEN
         DO k=1,kmax
           hzutx(k,i,j)=f_hzu_x(hx(i,j),MAX(xe(i,j),xe(i+1,j)),k,i,j)
         END DO
       END IF
   ! dz/ds at location hy (j+1/2) and depth of v or w
       IF ( hty(i,j)<=0.0_rsh .AND. ( hy(i,j)+xe(i,j) > 0.0_rsh .OR. hy(i,j)+xe(i,j+1) > 0.0_rsh )) THEN
         DO k=1,kmax
           hzuty(k,i,j)=f_hzu_y(hy(i,j),MAX(xe(i,j),xe(i,j+1)),k,i,j)
         END DO
       END IF
     END DO

   ! Southern boundary following
   IF_MPI (ljmin==jmin) THEN
     htx(i,jmin)=hx(i,jmin)+0.5_rsh*(xe(i,jmin+1)+xe(i+1,jmin+1))
     hty(i,jmin)=hy(i,jmin)+         xe(i,jmin+1)
   ! dz/ds at location hx (i+1/2) and depth u or w
     DO k=1,kmax
       hzutx(k,i,jmin)=f_hzu_x(hx(i,jmin),0.5_rsh*(xe(i,jmin+1)+xe(i+1,jmin+1)),k,i,jmin)
       hzuty(k,i,jmin)=f_hzu_y(hy(i,jmin),         xe(i,jmin+1)                ,k,i,jmin)
     END DO
     IF ( htx(i,jmin)<=0.0_rsh .AND. ( hx(i,jmin)+xe(i,jmin+1) > 0.0_rsh .OR. hx(i,jmin)+xe(i+1,jmin+1) > 0.0_rsh) ) THEN
       DO k=1,kmax
         hzutx(k,i,jmin)=f_hzu_x(hx(i,jmin),MAX(xe(i,jmin+1),xe(i+1,jmin+1)),k,i,jmin)
       END DO
     END IF
   ENDIF_MPI

   ! Northern boundary
     DO j=MAX0(jmax-obc_width+1,ljmin),ljmaxm1
       htx(i,j)=hx(i,j)+0.5_rsh*(xe(i,j)+xe(i+1,j))
       hty(i,j)=hy(i,j)+0.5_rsh*(xe(i,j)+xe(i,j+1))
       DO k=1,kmax
   ! dz/ds at location hx (i+1/2) and depth u or w
         hzutx(k,i,j)=f_hzu_x(hx(i,j),0.5_rsh*(xe(i,j)+xe(i+1,j)),k,i,j)
   ! dz/ds at location hy (j+1/2) and depth of v or w
         hzuty(k,i,j)=f_hzu_y(hy(i,j),0.5_rsh*(xe(i,j)+xe(i,j+1)),k,i,j)
       END DO
   ! dz/ds at location hx (i+1/2) and depth u or w
       IF ( htx(i,j)<=0.0_rsh .AND. ( hx(i,j)+xe(i,j) > 0.0_rsh .OR. hx(i,j)+xe(i+1,j) > 0.0_rsh) ) THEN
         DO k=1,kmax
           hzutx(k,i,j)=f_hzu_x(hx(i,j),MAX(xe(i,j),xe(i+1,j)),k,i,j)
         END DO
       END IF
   ! dz/ds at location hy (j+1/2) and depth of v or w
       IF ( hty(i,j)<=0.0_rsh .AND. ( hy(i,j)+xe(i,j) > 0.0_rsh .OR. hy(i,j)+xe(i,j+1) > 0.0_rsh )) THEN
         DO k=1,kmax
           hzuty(k,i,j)=f_hzu_y(hy(i,j),MAX(xe(i,j),xe(i,j+1)),k,i,j)
         END DO
       END IF
     END DO

   ! Northern boundary following
   IF_MPI (ljmax==jmax) THEN
     htx(i,jmax)=hx(i,jmax)+0.5_rsh*(xe(i,jmax)+xe(i+1,jmax))
     hty(i,jmax)=hy(i,jmax)+         xe(i,jmax)
   ! dz/ds at location hx (i+1/2) and depth u or w
     DO k=1,kmax
       hzutx(k,i,jmax)=f_hzu_x(hx(i,jmax),0.5_rsh*(xe(i,jmax)+xe(i+1,jmax)),k,i,jmax)
       hzuty(k,i,jmax)=f_hzu_y(hy(i,jmax),         xe(i,jmax)              ,k,i,jmax)
     END DO
     IF ( htx(i,jmax)<=0.0_rsh .AND. ( hx(i,jmax)+xe(i,jmax) > 0.0_rsh .OR. hx(i,jmax)+xe(i+1,jmax) > 0.0_rsh) ) THEN
       DO k=1,kmax
         hzutx(k,i,jmax)=f_hzu_x(hx(i,jmax),MAX(xe(i,jmax),xe(i+1,jmax)),k,i,jmax)
       END DO
     END IF
   ENDIF_MPI

   END DO
!$OMP END DO

   ! corners

!$OMP SINGLE
   ! upper left
   IF_MPI (ljmax==jmax .AND. limin==imin) THEN
     htx(imin,jmax)=hx(imin,jmax)+         xe(imin+1,jmax)
     hty(imin,jmax)=hy(imin,jmax)+         xe(imin+1,jmax)
   ! dz/ds at location hx (i+1/2) and depth u or w
     DO k=1,kmax
       hzutx(k,imin,jmax)=f_hzu_x(hx(imin,jmax),         xe(imin+1,jmax),k,imin,jmax)
       hzuty(k,imin,jmax)=f_hzu_y(hy(imin,jmax),         xe(imin+1,jmax),k,imin,jmax)
     END DO
   ENDIF_MPI

   ! lower right
   IF_MPI (ljmin==jmin .AND. limax==imax) THEN
     htx(imax,jmin)=hx(imax,jmin)+         xe(imax,jmin+1)
     hty(imax,jmin)=hy(imax,jmin)+         xe(imax,jmin+1)
   ! dz/ds at location hx (i+1/2) and depth u or w
     DO k=1,kmax
       hzutx(k,imax,jmin)=f_hzu_x(hx(imax,jmin),         xe(imax,jmin+1),k,imax,jmin)
       hzuty(k,imax,jmin)=f_hzu_y(hy(imax,jmin),         xe(imax,jmin+1),k,imax,jmin)
     END DO
   ENDIF_MPI

   ! upper right
   IF_MPI (ljmax==jmax .AND. limax==imax) THEN
     htx(imax,jmax)=hx(imax,jmax)+         xe(imax,jmax)
     hty(imax,jmax)=hy(imax,jmax)+         xe(imax,jmax)
   ! dz/ds at location hx (i+1/2) and depth u or w
     DO k=1,kmax
       hzutx(k,imax,jmax)=f_hzu_x(hx(imax,jmax),         xe(imax,jmax),k,imax,jmax)
       hzuty(k,imax,jmax)=f_hzu_y(hy(imax,jmax),         xe(imax,jmax),k,imax,jmax)
     END DO
   ENDIF_MPI

   ! lower left
   IF_MPI (ljmin==jmin .AND. limin==imin) THEN
     htx(imin,jmin)=hx(imin,jmin)+         xe(imin+1,jmin+1)
     hty(imin,jmin)=hy(imin,jmin)+         xe(imin+1,jmin+1)
   ! dz/ds at location hx (i+1/2) and depth u or w
     DO k=1,kmax
       hzutx(k,imin,jmin)=f_hzu_x(hx(imin,jmin),         xe(imin+1,jmin+1),k,imin,jmin)
       hzuty(k,imin,jmin)=f_hzu_y(hy(imin,jmin),         xe(imin+1,jmin+1),k,imin,jmin)
     END DO
   ENDIF_MPI
!$OMP END SINGLE

#endif

  END SUBROUTINE obc_init_height
