  SUBROUTINE bathy

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE bathy  ***
!&E
!&E ** Purpose : read the bathymetry and mean sea level file
!&E              in the historical DEL/AO format
!&E
!&E ** Description :
!&E
!&E ** Called by : main
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : hnivmoy,zlit
!&E
!&E ** Modified variables : hx_g, hy_g, h0_g, dhx, dhy
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2005-10-14
!&E       !  2009-04  (F. Dumas) generalized sigma
!&E       !  01-2011 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
!&E       !  09-2012 (S. Petton) key_siggen hc = MIN(h0_g(:,:),hc)
!&E       !  12-2016 (A. Fiandrino, B. Thouvenin) specificities key_lagoon
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE ionc4, ONLY : ionc4_openr,ionc4_read_xy,ionc4_close,&
                     ionc4_createfile,ionc4_createvar,ionc4_write_xy
   USE comvars2d


   USE toolmpi,   ONLY : MASTER
   USE obccombine,ONLY : l_obc_tide

   IMPLICIT NONE

!! * Local declarations
   LOGICAL                                       ::  l_wetdry
   CHARACTER(lchain)                             :: chain
   INTEGER                                       :: igg,idd,i,j,ierr_mpi
   REAL(kind=riosh),DIMENSION(imin:imax,jmin:jmax) :: zlit,hnmx,hnmy
   REAL(kind=riolg),dimension(imin:imax,jmin:jmax) :: tmp2drlg
   REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax) :: hx_t,hy_t

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


   CALL ionc4_openr(file_bathy)
   CALL ionc4_read_xy(file_bathy,name_in_longitude,tmp2drlg,imin,imax,jmin,jmax)
   lon2d_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,name_in_latitude,tmp2drlg,imin,imax,jmin,jmax)
   lat2d_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,name_in_longitude_u,tmp2drlg,imin,imax,jmin,jmax)
   lon2du_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,name_in_latitude_u,tmp2drlg,imin,imax,jmin,jmax)
   lat2du_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,name_in_longitude_v,tmp2drlg,imin,imax,jmin,jmax)
   lon2dv_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,name_in_latitude_v,tmp2drlg,imin,imax,jmin,jmax)
   lat2dv_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,name_in_longitude_f,tmp2drlg,imin,imax,jmin,jmax)
   lon2df_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,name_in_latitude_f,tmp2drlg,imin,imax,jmin,jmax)
   lat2df_g(:,:)=tmp2drlg(:,:)

!Read the grid characteristics from bathy_file
! Angle is defined from the x-axis (east) and positive anti-clockwise !
   CALL ionc4_read_xy(file_bathy,'angle',tmp2drlg,imin,imax,jmin,jmax)
   grid_angle2D_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,'dx',tmp2drlg,imin,imax,jmin,jmax)
   dx_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,'dy',tmp2drlg,imin,imax,jmin,jmax)
   dy_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,'dx_u',tmp2drlg,imin,imax,jmin,jmax)
   dxu_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,'dx_v',tmp2drlg,imin,imax,jmin,jmax)
   dxv_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,'dx_f',tmp2drlg,imin,imax,jmin,jmax)
   dxf_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,'dy_u',tmp2drlg,imin,imax,jmin,jmax)
   dyu_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,'dy_v',tmp2drlg,imin,imax,jmin,jmax)
   dyv_g(:,:)=tmp2drlg(:,:)
   CALL ionc4_read_xy(file_bathy,'dy_f',tmp2drlg,imin,imax,jmin,jmax)
   dyf_g(:,:)=tmp2drlg(:,:)

   CALL ionc4_close(file_bathy)



! Bathymetry of the current simulation
! ------------------------------------



   hnivmoy_g(:,:)=0.0_rsh

   CALL ionc4_openr(file_bathy)
   CALL ionc4_read_xy(file_bathy,name_bathy_hx,zlit,imin,imax,jmin,jmax, &
                      -rg_valmanq_io-REAL(h0fond,riosh))
   hx_g(:,:)=zlit(:,:)
   CALL ionc4_read_xy(file_bathy,name_bathy_hy,zlit,imin,imax,jmin,jmax, &
                      -rg_valmanq_io-REAL(h0fond,riosh))
   hy_g(:,:)=zlit(:,:)

hy_g(30:33,162)=-999
hx_g(30:33,162)=-999
hx_g(30,7:160)=-999 !kheili mohem
   IF (l_bathy_meanlev) THEN                                        ! read mean level
!CALL ionc4_read_xy(file_bathy,name_bathy_meanlev,zlit,imin,imax,jmin,jmax)
!hnivmoy_g(:,:)=zlit(:,:)
     CALL ionc4_read_xy(file_bathy,'hniv_moyU',hnmx,imin,imax,jmin,jmax)
     CALL ionc4_read_xy(file_bathy,'hniv_moyV',hnmy,imin,imax,jmin,jmax)
     hnmx=MAX(0.0,hnmx)
     hnmy=MAX(0.0,hnmy)
     hnivmoy_g(:,:)=0.5*(hnmx+hnmy)
   END IF




!#ifdef key_lagoon_spec
!
!! C'est ici qu'on peut creuser les rivieres dans la bathy
! ex :on creuse le Pallas sur Thau
! hy_g(137,64:65)=0.6_rsh
! h0_g(137,65:66)=0.6_rsh
!#endif



hnmx(10,175)=-999   !?????????
hx_g(463:464,69:71)=0.5_rsh   !?????????
hy_g(463:464,69:71)=0.5_rsh   !?????????
!hx_g(31,7)=15.
!hy_g(31,7)=15.
!hx_g(31,160)=15.
!hy_g(31,160)=15.


   WHERE (hx_g(:,:) < hminim) hx_g(:,:) = -valmanq -h0fond
   WHERE (hy_g(:,:) < hminim) hy_g(:,:) = -valmanq -h0fond
   hx_g(:,:)=hx_g(:,:)+hnmx(:,:)-mslshift
   hy_g(:,:)=hy_g(:,:)+hnmy(:,:)-mslshift
   WHERE (hx_g(:,:) < -valmanq+20.0_rsh) hx_g(:,:) = -valmanq -h0fond
   WHERE (hy_g(:,:) < -valmanq+20.0_rsh) hy_g(:,:) = -valmanq -h0fond




! modif bathy ou les mailles ne sont pas definies en lon & lat
   where (lon2du_g <= -valmanq) hx_g = -valmanq-h0fond
   where (lon2dv_g <= -valmanq) hy_g = -valmanq-h0fond
   h0_g(:,:)=-valmanq

   DO j=jmin+1,jmax
      DO i=imin+1,imax

        IF ( lon2du_g(i,j) > -valmanq .or. lat2dv_g(i,j) > -valmanq .or. &
             lon2du_g(i-1,j) > -valmanq .or. lat2dv_g(i,j-1) > -valmanq ) &

            h0_g(i,j)=MAX(hx_g(i,j),hx_g(i-1,j),hy_g(i,j),hy_g(i,j-1))+h0fond
     END DO
   END DO

!   IF (no==0) THEN
!      PRINT*, 'WARNING : wall at the eastern boundary for rank 0'
!      PRINT*, 'i.e. the Channel is closed'
!      PRINT*, 'this correction is available for the default rank 0 over the North East Atlantic'
!      PRINT*, 'See bathy.F90'
!      PRINT*, ' '
!      h0_g(imax,:)=-valmanq
!      hx_g(imax-1,:)=-valmanq-h0fond
!   END IF

   DO j=jmin+1,jmax
      h0_g(imin,j)=MAX(hx_g(imin,j),hy_g(imin,j),hy_g(imin,j-1))+h0fond
   END DO
   DO i=imin+1,imax
      h0_g(i,jmin)=MAX(hx_g(i,jmin),hx_g(i-1,jmin),hy_g(i,jmin))+h0fond
   END DO
   h0_g(imin,jmin)=MAX(hx_g(imin,jmin),hy_g(imin,jmin))+h0fond
   WHERE (h0_g(:,:) < -valmanq+20.0_rsh) h0_g(:,:) = -valmanq

!  redefine hx_g and hy_g
   DO i=imin,imax-1
      DO j=jmin,jmax-1
        IF (hx_g(i,j) >= hrdef) hx_g(i,j)=(h0_g(i,j)+h0_g(i+1,j))/2.0_rsh
        IF (hy_g(i,j) >= hrdef) hy_g(i,j)=(h0_g(i,j)+h0_g(i,j+1))/2.0_rsh
      END DO
   END DO



   hx_t(:,:)=100000.0_rsh
   hy_t(:,:)=100000.0_rsh
   WHERE (hx_g(:,:) /= -valmanq -h0fond) hx_t(:,:)=hx_g(:,:)
   WHERE (hy_g(:,:) /= -valmanq -h0fond) hy_t(:,:)=hy_g(:,:)
   dhxy_g(:,:)=1.0_rsh
   DO j=jmin+1,jmax
      DO i=imin+1,imax
         dhxy_g(i,j)=MAX(hx_g(i,j),hx_g(i-1,j),hy_g(i,j),hy_g(i,j-1))- &
                     MIN(hx_t(i,j),hx_t(i-1,j),hy_t(i,j),hy_t(i,j-1))+epsilon
      END DO
   END DO

   dhx_g(:,:)=0.0_rsh
   dhy_g(:,:)=0.0_rsh
   WHERE (hx_g(:,:) < hrdef) dhx_g(:,:)=dhj
   WHERE (hy_g(:,:) < hrdef) dhy_g(:,:)=dhj



! Make sure the bathymetric parameter are correctly initialized
! -------------------------------------------------------------

   l_wetdry=.False.
   DO j=jmin+1,jmax
     DO i=imin+1,imax
       IF (h0_g(i,j) <= 0.0_rsh .AND. h0_g(i,j) /= -valmanq) THEN
         l_wetdry=.true.
       ENDIF
     END DO
   END DO
   IF (l_wetdry .AND. h0fond <= 0.0_rsh) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR : There are wet-drying cells'
         WRITE(ierrorlog,*) 'you MUST specify h0fond > 0.'
         WRITE(ierrorlog,*) 'The simulation is stopped'
         CALL MPI_FINALIZE(IERR_MPI)
         STOP
   ENDIF
   IF (l_wetdry .AND. hm<h0fond) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR : There are wet-drying cells'
     WRITE(ierrorlog,*) 'you MUST specify hm >= h0fond, at the moment hm=',hm,' h0fond=',h0fond
     WRITE(ierrorlog,*) 'The simulation is stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   ENDIF
   IF (l_wetdry .AND. dhj<=0.0_rsh) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR : There are wet-drying cells'
     WRITE(ierrorlog,*) 'you MUST specify dhj > 0., at the moment dhj=',dhj
     WRITE(ierrorlog,*) 'The simulation is stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   ENDIF

! Save the bathymetry
! -------------------

   IF (l_bathy_save) THEN

     IF (MASTER) THEN
       CALL ionc4_createfile('bathy_verif.nc',lon2d_g,lat2d_g, &
                             imin,imax,1, &
                             jmin,jmax,1, &
                             1,kmax,1,kmax,0,l_out_nc4par=.false.)
       CALL ionc4_createvar('bathy_verif.nc','longitude_u',                            &
                            units='degrees_east',long_name='longitude at u location',   &
                            standard_name='longitude_at_u_location',location_in_cell='u', &
                            valid_min=-180.0_riolg,valid_max=180.0_riolg,              &
                            dims='xy',fill_value=dg_valmanq_io,l_auxcoordinate=.false.)
       tmp2drlg(:,:)=lon2du_g(:,:)
       CALL ionc4_write_xy('bathy_verif.nc','longitude_u',tmp2drlg,imin,imax,jmin,jmax)
       CALL ionc4_createvar('bathy_verif.nc','latitude_u',                             &
                            units='degrees_north',long_name='latitude at u location',   &
                            standard_name='latitude_at_u_location',location_in_cell='u',  &
                            valid_min=-90.0_riolg,valid_max=90.0_riolg,                &
                            dims='xy',fill_value=dg_valmanq_io,l_auxcoordinate=.false.)
       tmp2drlg(:,:)=lat2du_g(:,:)
       CALL ionc4_write_xy('bathy_verif.nc','latitude_u',tmp2drlg,imin,imax,jmin,jmax)
       CALL ionc4_createvar('bathy_verif.nc','longitude_v',                            &
                            units='degrees_east',long_name='longitude at v location',   &
                            standard_name='longitude_at_v_location',location_in_cell='v', &
                            valid_min=-180.0_riolg,valid_max=180.0_riolg,              &
                            dims='xy',fill_value=dg_valmanq_io,l_auxcoordinate=.false.)
       tmp2drlg(:,:)=lon2dv_g(:,:)
       CALL ionc4_write_xy('bathy_verif.nc','longitude_v',tmp2drlg,imin,imax,jmin,jmax)
       CALL ionc4_createvar('bathy_verif.nc','latitude_v',                             &
                            units='degrees_north',long_name='latitude at v location',   &
                            standard_name='latitude_at_v_location',location_in_cell='v',  &
                            valid_min=-90.0_riolg,valid_max=90.0_riolg,                &
                            dims='xy',fill_value=dg_valmanq_io,l_auxcoordinate=.false.)
       tmp2drlg(:,:)=lat2dv_g(:,:)
       CALL ionc4_write_xy('bathy_verif.nc','latitude_v',tmp2drlg,imin,imax,jmin,jmax)
       CALL ionc4_createvar('bathy_verif.nc','longitude_f',                            &
                            units='degrees_east',long_name='longitude at f location',   &
                            standard_name='longitude_at_f_location',location_in_cell='f', &
                            valid_min=-180.0_riolg,valid_max=180.0_riolg,              &
                            dims='xy',fill_value=dg_valmanq_io,l_auxcoordinate=.false.)
       tmp2drlg(:,:)=lon2du_g(:,:)
       CALL ionc4_write_xy('bathy_verif.nc','longitude_f',tmp2drlg,imin,imax,jmin,jmax)
       CALL ionc4_createvar('bathy_verif.nc','latitude_f',                             &
                            units='degrees_north',long_name='latitude at f location',   &
                            standard_name='latitude_at_f_location',location_in_cell='f',  &
                            valid_min=-90.0_riolg,valid_max=90.0_riolg,                &
                            dims='xy',fill_value=dg_valmanq_io,l_auxcoordinate=.false.)
       tmp2drlg(:,:)=lat2dv_g(:,:)
       CALL ionc4_write_xy('bathy_verif.nc','latitude_f',tmp2drlg,imin,imax,jmin,jmax)

       IF (l_obc_tide) THEN
         chain='model_sea_floor_depth_below_sea_level'
       ELSE
         chain='model_sea_floor_depth_below_geoid'
       END IF
       CALL ionc4_createvar('bathy_verif.nc','H0',units='m',                     &
                            long_name='bathymetry relative to the mean level',   &
                            standard_name = TRIM(chain),                         &
                            fill_value=-rg_valmanq_io,dims='xy')
       zlit(:,:)=h0_g(:,:)
       CALL ionc4_write_xy('bathy_verif.nc','H0',zlit,imin,imax,jmin,jmax)
       IF (l_obc_tide) THEN
         chain='model_sea_floor_depth_below_sea_level_at_u_location'
       ELSE
         chain='model_sea_floor_depth_below_geoid_at_u_location'
       END IF
       CALL ionc4_createvar('bathy_verif.nc','HX',units='m',location_in_cell='u',              &
                            long_name = 'bathymetry at u-location relative to the mean level', &
                            standard_name = TRIM(chain),                         &
                            fill_value=-rg_valmanq_io-REAL(h0fond,riosh),dims='xy')
       zlit(:,:)=hx_g(:,:)
       CALL ionc4_write_xy('bathy_verif.nc','HX',zlit,imin,imax,jmin,jmax)
       IF (l_obc_tide) THEN
         chain='model_sea_floor_depth_below_sea_level_at_v_location'
       ELSE
         chain='model_sea_floor_depth_below_geoid_at_v_location'
       END IF
       CALL ionc4_createvar('bathy_verif.nc','HY',units='m',location_in_cell='v',              &
                            long_name = 'bathymetry at v-location relative to the mean level', &
                            standard_name = TRIM(chain),                         &
                            fill_value=-rg_valmanq_io-REAL(h0fond,riosh),dims='xy')
       zlit(:,:)=hy_g(:,:)
       CALL ionc4_write_xy('bathy_verif.nc','HY',zlit,imin,imax,jmin,jmax)

       CALL ionc4_close('bathy_verif.nc')
     ENDIF
   END IF


! key_offline


   CALL ionc4_close(file_bathy)

   !PRINT*, 'exit BATHY routine'

  END SUBROUTINE bathy
