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

  SUBROUTINE grid_metric

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE grid_metric  ***
   !&E
   !&E ** Purpose : estimate of grid parametrization
   !&E
   !&E ** Description : estimate of surface cell, begining and ending cells,
   !&E                  sigma repartition,...
   !&E
   !&E ** Called by : main.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !         (J.-C. Salomon, P. Lazure)  Original code
   !&E       !  2007-09 (V. Garnier) homogeneneity of dx dy...
   !&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
   !&E       !  01-2011 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
   !&E       !  04-2015 (R. Le Gendre) Estimate of the circle distance from 2 longitude/latitude locations (dist_circle)
   !&E       !  06-2015 (V. Garnier) MPILARGERLOOP for grid variables to avoid problem with sleeping cpus (and 3D exchanges in future?) 
   !&E       !  12-2016 (A. Fiandrino, B. Thouvenin. ) introducton of specificities key_lagoon (govlet modified by user) 
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE ionc4,     ONLY : ionc4_open,ionc4_createvar,ionc4_write_xy,ionc4_close
   USE comvars2d, ONLY : dx,dy,dfiac,goulet,govlet,cosg,sing,cos2g,sin2g
   USE comvars2d, ONLY : dx_u,dx_v,dx_f,dy_v,dy_f,dy_u,l_corio_cst,surf,ig,id,h0,radearth
   USE comvars2d, ONLY : sinfiu,sin2fi,coscafi,fv,tfivsr,l_modele2d,iwarnlog
   USE comvars2d, ONLY : cosfiu,tfiusr,fu,fisac,dgac,gwac,lonmin,latmin,lonmax,latmax
   USE comvars2d, ONLY : ierrorlog,l_bathy_save,lon2du,lon2dv,lat2du,lat2dv,lon2d,lat2d
   USE toolgeom,  ONLY : dist_circle
#ifdef key_lagoon
   USE comvars2d,  ONLY: nbgrau,igrau,jgrau,l_open_grau,type_grau,dxac
#endif
#ifdef key_grid_curv
  USE comvars2d, ONLY : dx_g,dy_g,dxu_g,dyu_g,dxv_g,dyv_g,dxf_g,dyf_g
#endif
#if defined key_ts_diffh_geo || defined key_ts_diffh_rho
   USE comvars2d, ONLY : jbound_omp
#endif
   USE comvars2d, ONLY : lon2d_g,lat2d_g,lon2du_g,lat2du_g,lon2dv_g,lat2dv_g,lon2df_g,lat2df_g
   USE comvars3d, ONLY : l_equisig,sig,dsigu,dsigw,sigw
#if defined key_siggen || defined key_gencoord
   USE comsiggen, ONLY : hc_sig,hcx_sig,hcy_sig,Csu_sig,Csw_sig,dcusds,dcwsds,b_sig,theta_sig,dcsigdxy
#if ! defined key_siggen_g2
   USE comsiggen, ONLY : hzdcusds,hzdcwsds,hzdcusds_x,hzdcwsds_x,hzdcusds_y,hzdcwsds_y
#endif
   USE comvars2d, ONLY : hx,hy
#endif
#if defined key_profil_back
   USE comvars2d, ONLY : comdir
   USE comvars3d, ONLY : rho2,sali_back,temp_back,l_stateeq_lin
   USE toolgeom,  ONLY : f_zu
#endif
#if defined key_critere_inconsistence
   USE comvars2d, ONLY : h0_g,hx_g,hy_g,h0fond
   USE toolgeom,  ONLY : f_zw
#endif

   USE_MPI toolmpi,   ONLY : MASTER,ex_i_rsh,ex_j_rsh

   IMPLICIT NONE

   !! * Arguments

   !! * Local declarations
   INTEGER        :: i,j,k,ixl,iyl,nb
#if defined key_ts_diffh_geo || defined key_ts_diffh_rho
   INTEGER        :: nb_wetpoint,OMP_GET_NUM_THREADS,nb_point,nb_thread
#endif
   REAL(KIND=rlg) :: rhox,rhoy
   REAL(KIND=rsh) :: fi,g,radeg
#ifdef key_grid_verif
   INTEGER        :: IERR_MPI
   REAL(kind=rsh),dimension(limin:limax,ljmin:ljmax) :: tmp2drsh
#endif

#if defined key_profil_back
   REAL(KIND=rsh),dimension(300) :: sal_profil,temp_profil,depth
   REAL(KIND=rsh) :: sal_profill,temp_profill,depthl,zr
   INTEGER :: kk,kback,eof_back
   REAL(KIND=rsh)           :: TP,S,P
   REAL(KIND=rlg)           :: row,rost,c,pcc
#endif
#if defined key_critere_inconsistence
   REAL(KIND=rsh) :: zikp12,zikm12,zip1kp12,zip1km12 
   REAL(KIND=rsh) :: zjkp12,zjkm12,zjp1kp12,zjp1km12 
   REAL(KIND=rsh) :: crit 
#endif
#ifdef key_lagoon
   REAL(KIND=rsh) :: etr
#endif

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

   PRINT_DBG*, 'enter GRIDMETRIC routine'

#ifdef key_agrif
   IF (.NOT.agrif_root()) THEN
     rhox=Agrif_Rhox()
     rhoy=Agrif_Rhoy()

     ixl=Agrif_ix()
     iyl=Agrif_iy()

     dfiac = (Agrif_Parent(dfiac))/rhoy
     fisac = Agrif_Parent(fisac)+(iyl-0.5_rlg-0.5_rlg/rhoy) &
             *Agrif_Parent(dfiac)
     dgac = ( Agrif_Parent(dgac) )/rhox
     gwac = Agrif_Parent(gwac)+(ixl-0.5_rlg-0.5_rlg/rhox) &
             *Agrif_Parent(dgac)
   END IF
#endif

!$OMP PARALLEL DEFAULT(none) PRIVATE(i,j,k,radeg,g,           &
#ifdef key_lagoon
!$OMP& etr,                                                   &
#endif
!$OMP& fi)                                                    &
!$OMP& SHARED(dx,dy,govlet,goulet,cosg,sing,cos2g,sin2g       &
!$OMP& ,lonmin,latmin,lonmax,latmax                           &
#ifdef key_agrif
!$OMP& ,imin,imax,jmin,jmax                                   &
#endif
#ifdef key_lagoon
!$OMP& ,nbgrau,igrau,jgrau,type_grau,l_open_grau,dxac         &
#endif
!$OMP& ,l_corio_cst,fu,tfiusr,cosfiu,sinfiu,sin2fi,coscafi    &
!$OMP& ,dx_u,dx_v,dx_f,dy_v,dy_f,dy_u                         &
!$OMP& ,lon2du,lat2du,lat2dv,lon2d,lat2d                      &
!$OMP& ,lon2d_g,lat2d_g,lon2du_g,lat2du_g,lon2dv_g            &
!$OMP& ,lat2dv_g,lon2df_g,lat2df_g                            &
!$OMP& ,liminm2,limaxp2,ljminm2,ljmaxp2                       &
#if defined key_MPI_2D
!$OMP& ,MASTER                                                &
#endif
!$OMP& ,fv,tfivsr,surf,sig,dsigu,l_modele2d,iwarnlog          &
!$OMP& ,l_equisig,dsigw,sigw,fisac                            &
#if defined key_siggen || defined key_gencoord || defined key_profil_back
!$OMP& ,h0                                                    &
#endif
#if defined key_siggen || defined key_gencoord
!$OMP& ,hc_sig,Csu_sig,Csw_sig,dcusds,dcwsds,b_sig,theta_sig  &
!$OMP& ,hcx_sig,hcy_sig,hx,hy,hzdcusds,hzdcwsds,dcsigdxy      &
!$OMP& ,hzdcusds_x,hzdcwsds_x,hzdcusds_y,hzdcwsds_y           &
#endif
#if defined key_ts_diffh_geo || defined key_ts_diffh_rho
!$OMP& ,nb_thread,ig,id                                       &
#endif
#ifdef key_grid_verif
!$OMP ,lon2dv                                                 &
#endif
#ifdef key_grid_curv
!$OMP ,h0,dx_g,dy_g,dxu_g,dyu_g,dxv_g,dyv_g,dxf_g,dyf_g          &
#endif
!$OMP& ,limin,limax,ljmin,ljmax,limaxp1,ljmaxp1,              &
!$OMP& ljminp1,liminm1,ljminm1,ierrorlog)

   radeg = pi/180.0_rsh

   ! Initially, we know the location of points hx and hy at each cell i,j :
   ! The longitudes and latitudes are read from file_bathy
   !
   !       ---- V(i,j) --- phi(i,j) 
   !       ^                  |                        <-- dx_v(i,j)--->
   !       |                  |                       |
   !       |                  |                       |
   !      dy(    XE(i,j)     U(i,j) and hx(i,j)       |<----dx(i,j)---->
   !       |                  |                       |
   !       |                  |                       |
   !       <---dx_v(i,j-1)--->|                       |-----------------
   !
   !
   !       |<-- dx_v(i,j)--->|            |--------|--------|
   !       |                 |            |        |        |
   !       |                 |            |        |        |
   !       |<----dx(i,j)---->|            |     dy(i,j)   dy_u(i,j)
   !       |                 |            |        |        |
   !       |                 |            |        |        |
   !       |-----------------|            |--------|--------| 
   !
   ! so we know dx and dy at the center of the cell and we estimate :
   ! dx_u and dy_u at location u or hx
   ! dx_v and dy_v at location v or hy

!-----------------------------------------------------------
! coriolis parameter, estimates of dx,dy at the different locations, etc...
!-----------------------------------------------------------

!MPILARGERLOOP
!$OMP DO SCHEDULE(STATIC)
   DO j=ljminm2,ljmaxp2
     DO i = liminm2,limaxp2

       goulet(i,j)=1.0_rsh
       govlet(i,j)=1.0_rsh

       ! T location
       g = REAL(lon2d(i,j),rsh)
       fi = REAL(lat2d(i,j),rsh)

       cosg(i,j) = COS(g*radeg)
       sing(i,j) = SIN(g*radeg)
       cos2g(i,j) = COS(2.0_rsh*g*radeg)
       sin2g(i,j) = SIN(2.0_rsh*g*radeg)

       IF (l_corio_cst) fi=fisac
       cosfiu(i,j) = COS(fi*radeg)
       sinfiu(i,j) = SIN(fi*radeg)
       sin2fi(i,j) = SIN(2.0_rsh*fi*radeg)
       coscafi(i,j) = cosfiu(i,j)*cosfiu(i,j)

#if !defined key_grid_curv
       IF (j==jmin) THEN
         dy(i,j)=dist_circle(lat2dv_g(i,j+1),lat2dv_g(i,j),lon2dv_g(i,j+1),lon2dv_g(i,j))
       ELSE
         dy(i,j)=dist_circle(lat2dv_g(i,j),lat2dv_g(i,j-1),lon2dv_g(i,j),lon2dv_g(i,j-1))
       END IF
       IF (i==imin) THEN
         dx(i,j)=dist_circle(lat2du_g(i+1,j),lat2du_g(i,j),lon2du_g(i+1,j),lon2du_g(i,j))
       ELSE
         dx(i,j)=dist_circle(lat2du_g(i,j),lat2du_g(i-1,j),lon2du_g(i,j),lon2du_g(i-1,j))
       END IF
       IF_CURV (dx(i,j)==0.0_rsh) dx(i,j)=dxmin
       IF_CURV (dy(i,j)==0.0_rsh) dy(i,j)=dxmin
#endif

       
       ! U location
       fi = REAL(lat2du(i,j),rsh)
       IF (l_corio_cst) fi=fisac
       ! 0.0001458 = 2 *omega with omega=2*pi/86164 (86164 = sidereal day)
       fu(i,j)=0.0001458_rsh*SIN(fi*radeg)
       tfiusr(i,j)=TAN(fi*radeg)/radearth
       IF (l_corio_cst) tfiusr(i,j)=0.0_rsh

#if !defined key_grid_curv
       IF (j==jmin) THEN
         dy_u(i,j)=dist_circle(lat2df_g(i,j+1),lat2df_g(i,j),lon2df_g(i,j+1),lon2df_g(i,j))
       ELSE
         dy_u(i,j)=dist_circle(lat2df_g(i,j),lat2df_g(i,j-1),lon2df_g(i,j),lon2df_g(i,j-1))
       END IF
       IF (i==imax) THEN
         dx_u(i,j)=dist_circle(lat2d_g(i,j),lat2d_g(i-1,j),lon2d_g(i,j),lon2d_g(i-1,j))
       ELSE
         dx_u(i,j)=dist_circle(lat2d_g(i+1,j),lat2d_g(i,j),lon2d_g(i+1,j),lon2d_g(i,j))
       END IF
       IF_CURV (dx_u(i,j)==0.0_rsh) dx_u(i,j)=dxmin
       IF_CURV (dy_u(i,j)==0.0_rsh) dy_u(i,j)=dxmin
#endif


       ! V location
       fi = REAL(lat2dv(i,j),rsh)
       IF (l_corio_cst) fi=fisac
       fv(i,j)=0.0001458_rsh*SIN(fi*radeg)
       tfivsr(i,j)=TAN(fi*radeg)/radearth
       IF (l_corio_cst) tfivsr(i,j)=0.0_rsh
#if !defined key_grid_curv
       IF (j==jmax) THEN
         dy_v(i,j)=dist_circle(lat2d_g(i,j),lat2d_g(i,j-1),lon2d_g(i,j),lon2d_g(i,j-1))
       ELSE
         dy_v(i,j)=dist_circle(lat2d_g(i,j+1),lat2d_g(i,j),lon2d_g(i,j+1),lon2d_g(i,j))
       END IF
       IF (i==imin) THEN
         dx_v(i,j)=dist_circle(lat2df_g(i+1,j),lat2df_g(i,j),lon2df_g(i+1,j),lon2df_g(i,j))
       ELSE
         dx_v(i,j)=dist_circle(lat2df_g(i,j),lat2df_g(i-1,j),lon2df_g(i,j),lon2df_g(i-1,j))
       END IF

       ! F location
       IF (j==jmax) THEN
         dy_f(i,j)=dist_circle(lat2du_g(i,j),lat2du_g(i,j-1),lon2du_g(i,j),lon2du_g(i,j-1))
       ELSE
         dy_f(i,j)=dist_circle(lat2du_g(i,j+1),lat2du_g(i,j),lon2du_g(i,j+1),lon2du_g(i,j))
       END IF
       IF (i==imax) THEN
         dx_f(i,j)=dist_circle(lat2dv_g(i,j),lat2dv_g(i-1,j),lon2dv_g(i,j),lon2dv_g(i-1,j))
       ELSE
         dx_f(i,j)=dist_circle(lat2dv_g(i+1,j),lat2dv_g(i,j),lon2dv_g(i+1,j),lon2dv_g(i,j))
       END IF
       IF_CURV (dx_f(i,j)==0.0_rsh) dx_f(i,j)=dxmin
       IF_CURV (dy_f(i,j)==0.0_rsh) dy_f(i,j)=dxmin
#endif

#ifdef key_grid_curv
       dx(i,j)=max(dxmin,dx_g(i,j))
       dy(i,j)=max(dxmin,dy_g(i,j))
       dx_u(i,j)=max(dxmin,dxu_g(i,j))
       dy_u(i,j)=max(dxmin,dyu_g(i,j))
       dx_v(i,j)=max(dxmin,dxv_g(i,j))
       dy_v(i,j)=max(dxmin,dyv_g(i,j))
       dx_f(i,j)=max(dxmin,dxf_g(i,j))
       dy_f(i,j)=max(dxmin,dyf_g(i,j))
#endif

     END DO
   END DO
!$OMP END DO

#ifdef key_lagoon
! le 14/10/09 : Retrecissement du canal et betonnage bord
! pour limiter les débits dans le grau : etr de base =40/dxac
!$OMP SINGLE
   etr = 30.0_rsh/ dxac   
! canal droit pour thau  
   i=igrau(1)
   DO j=1,10
     IF_MPI ( (j >= ljminm2) .AND. (j <= ljmaxp2) .AND. (i >= liminm2) .AND. (i <= limaxp2) ) &
     govlet(i,j)=etr
   END DO
!$OMP END SINGLE
#endif

#ifdef key_nocoriolis
!$OMP SINGLE
   fu(:,:)=0.0_rsh
   fv(:,:)=0.0_rsh
!$OMP END SINGLE
#endif

   !-------------------------------------------------------
   ! boundaries for meteorological forcings :
   !--------------------------------------------------------
   latmax=maxval(lat2dv,lat2dv>-valmanq)
   latmin=minval(lat2d,lat2d>-valmanq)
   lonmin=minval(lon2d,lon2d>-valmanq)
   lonmax=maxval(lon2du,lon2du>-valmanq)
   !print*,"lonmin,lonmax,latmin,latmax",lonmin,lonmax,latmin,latmax

   IF_MPI (MASTER) THEN
!cval si idiot que cela... Il me semble qu il fallait s assurer que les valeurs
!ne soient pas trop petites pour les envoyer a l interpolation ou pour reperer
!indices ou... 
!   IF ( ABS(latmax) < epsilon .OR. &
!        ABS(latmin) < epsilon .OR. &
!        ABS(lonmax) < epsilon .OR. &
!        ABS(lonmin) < epsilon ) THEN
!$OMP SINGLE
   IF ( (latmax-latmin) < epsilon .OR. &
        (lonmax-lonmin) < epsilon ) THEN
     WRITE(iwarnlog,*) " "
     WRITE(iwarnlog,*) "ERROR: routine BATHY.F90 ?"
     WRITE(iwarnlog,*) " "
     WRITE(iwarnlog,*) 'lonmin=',lonmin
     WRITE(iwarnlog,*) 'lonmax=',lonmax
     WRITE(iwarnlog,*) 'latmax=',latmax
     WRITE(iwarnlog,*) 'latmin=',latmin
     WRITE(iwarnlog,*) "Make sure the max/min latitude and longitude are correct"
     WRITE(iwarnlog,*) " "
   END IF
!$OMP END SINGLE
   ENDIF_MPI

!$OMP DO SCHEDULE(STATIC) 
   DO j = ljminp1,ljmax
     DO i = limin,limax
       IF (l_corio_cst) THEN
         surf(i,j)=dy(i,j)*govlet(i,j)*dx(i,j)*goulet(i,j)
       ELSE
#ifdef key_grid_curv
         IF (h0(i,j)== -valmanq) THEN
           surf(i,j)=-valmanq
         ELSE
           !surf(i,j)=goulet(i,j)*dy(i,j)*(govlet(i,j)*dx_v(i,j)+govlet(i,j-1)*dx_v(i,j-1))*0.5_rsh 
           surf(i,j)=goulet(i,j)*dy(i,j)*(govlet(i,j)*dx_v(i,j)+govlet(i,j-1)*dx_v(i,j-1))*0.5_rsh  &
               *SQRT(1.0_rsh-(0.5_rsh*(govlet(i,j)*dx_v(i,j)-govlet(i,j-1)*dx_v(i,j-1))/dy(i,j))**2)
         END IF
#else
           surf(i,j)=goulet(i,j)*dy(i,j)*(govlet(i,j)*dx_v(i,j)+govlet(i,j-1)*dx_v(i,j-1))*0.5_rsh  &
               *SQRT(1.0_rsh-(0.5_rsh*(govlet(i,j)*dx_v(i,j)-govlet(i,j-1)*dx_v(i,j-1))/dy(i,j))**2)

#endif
       END IF
     END DO
   END DO
!$OMP END DO
!$OMP SINGLE
   IF_MPI (ljmin==jmin) THEN
     surf(limin:limax,jmin)=surf(limin:limax,jmin+1)
   ENDIF_MPI
!$OMP END SINGLE
OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,surf(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,surf(liminm2:limaxp2,ljminm2:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(surf)

   ! sigma repartition
   !------------------

!$OMP single
#if defined key_siggen || defined key_gencoord
  ! pour l instant on calcule les coefs de contraction avec sig regulier mais a terme on pourra choisir regulier ou non
   l_equisig=.TRUE.
#endif
   IF (l_equisig) THEN
     DO k=1,kmax
       sig(k)=(REAL(k-kmax,rsh)-0.5_rsh)/REAL(kmax,rsh)
     END DO
   END IF

   sig(kmax+1)=-sig(kmax)
   sig(0)=ABS(sig(1))-2.0_rsh

   DO k=1,kmax
     IF (.NOT.(sig(k) > sig(k-1) .AND. sig(k) < sig(k+1))) THEN
       WRITE(ierrorlog,*) ' '
       WRITE(ierrorlog,*) 'ERROR: paraspec.txt input file'
       WRITE(ierrorlog,*) 'ERROR: test in grid_metric.F90 routine'
       WRITE(ierrorlog,*) 'Sigma level must be defined with increasing number.'
       WRITE(ierrorlog,*) 'There are not in paraspec.txt namelist.'
       STOP
     END IF
   END DO

   DO k=1,kmax
     !cval  dsigu(k)=0.5*(sig(k+1)-sig(k-1))
     dsigw(k)=sig(k+1)-sig(k)
     sigw(k)=0.5_rsh*(sig(k)+sig(k+1))
   END DO
   sigw(kmax+1)=-sigw(kmax)
   sigw(0)=-1.0_rsh
   IF (l_modele2d) THEN
     dsigu(1)=1.0_rsh
     dsigu(2:kmax)=0.0_rsh
   END IF

   IF (.NOT.l_modele2d) THEN
     DO k=1,kmax
       dsigu(k)=sigw(k)-sigw(k-1)
     END DO
   END IF

#if defined key_siggen || defined key_gencoord
#ifdef key_siggen_g2
   DO k=0,kmax+1
     IF (theta_sig > 0.0_rsh) THEN
       Csu_sig(k)=(1.0_rsh-COSH(theta_sig*sig(k)))/(COSH(theta_sig)-1.0_rsh)
     ELSE
       Csu_sig(k)=-sig(k)*sig(k)
     END IF
     IF (b_sig > 0.0_rsh) &
       Csu_sig(k)=(EXP(b_sig*Csu_sig(k))-1.0_rsh)/(1.0_rsh-EXP(-b_sig))
   !END DO
   !DO k=0,kmax
     IF (theta_sig > 0.0_rsh) THEN
       Csw_sig(k)=(1.0_rsh-COSH(theta_sig*sigw(k)))/(COSH(theta_sig)-1.0_rsh)
     ELSE
       Csw_sig(k)=-sigw(k)*sigw(k)
     END IF
     IF (b_sig > 0.0_rsh) &
       Csw_sig(k)=(EXP(b_sig*Csw_sig(k))-1.0_rsh)/(1.0_rsh-EXP(-b_sig))
   END DO
#else
   DO k=0,kmax+1
!   Csu_sig(k)=(1.0_rsh-b_sig)*SINH(theta_sig*(sig(k)-1.0_rsh))/SINH(theta_sig) &
!         + 0.5_rsh*b_sig*(TANH(theta_sig*(sig(k)-0.5_rsh))-TANH(0.5_rsh*theta_sig))/TANH(theta_sig/2.0_rsh)
     Csu_sig(k)=(1.0_rsh-b_sig)*SINH(theta_sig*sig(k))/SINH(theta_sig) &
         + 0.5_rsh*b_sig*(TANH(theta_sig*(sig(k)+0.5_rsh))-TANH(0.5_rsh*theta_sig))/TANH(theta_sig/2.0_rsh)
   !END DO
   !DO k=0,kmax
!   Csw_sig(k)=(1.0_rsh-b_sig)*SINH(theta_sig*(sigw(k)-1.0_rsh))/SINH(theta_sig) &
!         + 0.5_rsh*b_sig*(TANH(theta_sig*(sigw(k)-0.5_rsh))-TANH(0.5_rsh*theta_sig))/TANH(theta_sig/2.0_rsh)
     Csw_sig(k)=(1.0_rsh-b_sig)*SINH(theta_sig*sigw(k))/SINH(theta_sig) &
         + 0.5_rsh*b_sig*(TANH(theta_sig*(sigw(k)+0.5_rsh))-TANH(0.5_rsh*theta_sig))/TANH(theta_sig/2.0_rsh)

     dcsigdxy(k) = (1.0_rsh-b_sig)*theta_sig*COSH(theta_sig*sig(k))/SINH(theta_sig) &
         + 0.5_rsh*theta_sig*b_sig*(1.0_rsh-(TANH(theta_sig*(sigw(k)+0.5_rsh)))**2)/TANH(theta_sig/2.0_rsh)
   END DO
#endif

   DO k=1,kmax
     dcusds(k)=(Csw_sig(k)-Csw_sig(k-1))/(sigw(k)-sigw(k-1))
   END DO

   DO k=1,kmax-1
     dcwsds(k)=(Csu_sig(k+1)-Csu_sig(k))/(sig(k+1)-sig(k))
   END DO
   dcwsds(kmax)=(0.0_rsh-csu_sig(kmax))/(0.0_rsh-sig(kmax))
   IF_MPI (MASTER) THEN
   DO k=1,kmax
     PRINT '(i4,1x,5(f6.3,1x))',k,sig(k),sigw(k),csw_sig(k)+1,csu_sig(k)+1,dcusds(k)
   END DO
   PRINT*,'sum(dcusds(1:kmax)*dsigu(1:kmax)) must be equal to 1'
   PRINT*,SUM(dcusds(1:kmax)*dsigu(1:kmax))
   ENDIF_MPI

#if ! defined key_siggen_g2
   DO j=ljminm1,ljmaxp1
     DO i=liminm1,limaxp1
       DO k=1,kmax
         hzdcusds(k,i,j)=hc_sig(i,j)+(h0(i,j)-hc_sig(i,j))*dcusds(k)
         hzdcusds_x(k,i,j)=hcx_sig(i,j)+(hx(i,j)-hcx_sig(i,j))*dcusds(k)
         hzdcusds_y(k,i,j)=hcy_sig(i,j)+(hy(i,j)-hcy_sig(i,j))*dcusds(k)
       END DO
     END DO
   END DO
   DO j=ljmin,ljmaxp1
     DO i=limin,limaxp1
       DO k=1,kmax
         hzdcwsds(k,i,j)=hc_sig(i,j)+(h0(i,j)-hc_sig(i,j))*dcwsds(k)
         hzdcwsds_x(k,i,j)=hcx_sig(i,j)+(hx(i,j)-hcx_sig(i,j))*dcwsds(k)
         hzdcwsds_y(k,i,j)=hcy_sig(i,j)+(hy(i,j)-hcy_sig(i,j))*dcwsds(k)
       END DO
     END DO
   END DO
#endif

#endif

!$OMP END single

#if defined key_ts_diffh_geo || defined key_ts_diffh_rho
!$ nb_thread=OMP_GET_NUM_THREADS()
#endif

!$OMP END PARALLEL

#if defined key_ts_diffh_geo || defined key_ts_diffh_rho
  ! estimate limits for switching thread with OMP : jbound_omp
  ! with balancing the number of wet cells of each thread

! initialization of locales variables + jbound_omp (public)
!$ ALLOCATE( jbound_omp(0:nb_thread) )
!$ jbound_omp(0)=0
!$ j=jmin
!$ nb=0
!$ nb_point=0
! number of wet cells by thread
!$ nb_wetpoint=INT( SUM(id(:)-ig(:)+1)/nb_thread ) + 1

! set jbound_omp(num thread) to limit of switching
!$ DO WHILE(j<=jmax)
!$   IF ( j==jmax .OR. nb_point>=nb_wetpoint ) THEN
!$     nb=nb+1
!$     jbound_omp(nb)=j
!$     nb_point=0
!$   ELSE IF ( nb_point<nb_wetpoint ) THEN
!$     nb_point=nb_point+id(j)-ig(j)+1
!$   END IF
!$   j=j+1
!$ END DO
!$ PRINT*,jbound_omp,'jbound_omp'
#endif

#ifdef key_profil_back
  IF_MPI(MASTER) THEN
    print*,'Internal pressure gradient will be estimated with density anomaly (real '
    print*,'density substract from a typical profil "profil_back.txt"'
  ENDIF_MPI
#ifdef key_oasis
  open(623,file="./"//TRIM(comdir)//"/inputs/profil_back.txt",form="formatted")
#else
  open(623,file="../"//TRIM(comdir)//"/inputs/profil_back.txt",form="formatted")
#endif
   eof_back=0
   kback=0
   do while(eof_back==0)
    read(623,*,iostat=eof_back)depthl,sal_profill,temp_profill
    if(eof_back==0) then
    kback=kback+1
    depth(kback)=depthl
    sal_profil(kback)=sal_profill
    temp_profil(kback)=temp_profill
    endif
   end do
  close(623)

   DO j=ljminm1,ljmaxp2
     DO i=liminm1,limaxp2
       IF (h0(i,j) > -valmanq) THEN
         DO k=1,kmax
           zr=-f_zu(h0(i,j),0.0_rsh,k,i,j) 
           temp_back(k,i,j)=temp_profil(1)
           sali_back(k,i,j)=sal_profil(1)
           kk=1
           IF (zr > depth(kk)) THEN
             do  while(depth(kk) < zr .and. kk < kback)
               kk=kk+1
             end do
             temp_back(k,i,j) = ((zr-depth(kk-1))*temp_profil(kk)+(depth(kk)-zr)*temp_profil(kk-1))/(depth(kk)-depth(kk-1))
             sali_back(k,i,j) = ((zr-depth(kk-1))*sal_profil(kk) +(depth(kk)-zr)* sal_profil(kk-1))/(depth(kk)-depth(kk-1))
           END IF
         END DO
       ELSE
         temp_back(:,i,j) = temp_profil(1)
         sali_back(:,i,j)  = sal_profil(1)
       END IF
       DO k=1,kmax
         P = 10.0_rsh
         tp=temp_back(k,i,j)
         S =sali_back(k,i,j)
   ! density of pure water
         row = 999.842594_rlg + (  6.793952d-02 + ( - 9.095290d-03 +                &
               (1.001685d-04 - 1.120083d-06*TP + 6.536332d-09*TP*TP)*TP )*TP  )*TP
   ! density of sea water at the surface
         rost = row + S * (   0.824493_rlg + (  -4.0899d-03 + ( 7.6438d-05 +        &
                          (-8.2467d-07+5.3875d-09*TP)*TP )*TP  )*TP   ) +           &
                S**1.5_rlg * ( -5.72466d-03 + (1.0227d-04-1.6546d-06*TP)*TP ) +     &
                S*S*4.8314d-04
   ! density of sea water at pressure p (decibars)
         c = 1449.2_rlg + 1.34_rlg*(S-35.0_rlg) + 4.55_rlg*TP - 0.045_rlg*TP*TP +   &
             0.00821_rlg*p + 15.0d-09*p*p
         pcc = p/(c*c)
         rho2(k,i,j) = rost + 1.0d4*(1.0_rlg-0.2_rlg*pcc)*pcc
       END DO
     END DO
   END DO
   IF (l_stateeq_lin) THEN
      DEALLOCATE(rho2)
   ELSE
      DEALLOCATE(sali_back,temp_back)
   END IF 
#endif /* profil_back */

!! Verification du critere d inconsistence hydrostatique !!
#if defined key_critere_inconsistence
#if defined key_siggen || defined key_gencoord
 IF_MPI (MASTER) THEN
  DO j=jmin+1,jmax-1
    DO i=imin+1,imax-1
     if(hx_g(i,j) /= -valmanq -h0fond .and. h0_g(i,j) > 50.0) then
       do k=1,kmax
         zikp12   = f_zw(h0_g(i  ,j),0.0_rsh,k  ,i  ,j) 
         zikm12   = f_zw(h0_g(i  ,j),0.0_rsh,k-1,i  ,j) 
         zip1kp12 = f_zw(h0_g(i+1,j),0.0_rsh,k  ,i+1,j)
         zip1km12 = f_zw(h0_g(i+1,j),0.0_rsh,k-1,i+1,j)
         if(zip1kp12-zip1km12+zikp12-zikm12/=0.0 ) crit=abs((zip1kp12+zip1km12-(zikp12+zikm12))/(zip1kp12-zip1km12+zikp12-zikm12))
         if(crit > 1) write(566,'(3i4,1x,f6.3,1x,3(e10.4,1x))') k,i,j,crit,h0_g(i,j),h0_g(i+1,j),hx_g(i,j)
       end do
     end if
     if(hy_g(i,j) /= -valmanq -h0fond .and. h0_g(i,j) > 50.0) then
       do k=1,kmax
         zjkp12   = f_zw(h0_g(i,j  ),0.0_rsh,k  ,i,j  ) 
         zjkm12   = f_zw(h0_g(i,j  ),0.0_rsh,k-1,i,j  ) 
         zjp1kp12 = f_zw(h0_g(i,j+1),0.0_rsh,k  ,i,j+1)
         zjp1km12 = f_zw(h0_g(i,j+1),0.0_rsh,k-1,i,j+1)
         if(zjp1kp12-zjp1km12+zjkp12-zjkm12/=0.0 ) crit=abs((zjp1kp12+zjp1km12-(zjkp12+zjkm12))/(zjp1kp12-zjp1km12+zjkp12-zjkm12))
         if(crit > 1) write(567,'(3i4,1x,f6.3,1x,3(e10.4,1x))') k,i,j,crit,h0_g(i,j),h0_g(i,j+1),hy_g(i,j)
       end do
     end if
    END DO
  END DO
 ENDIF_MPI
#endif
#endif


#ifdef key_grid_verif
   ! saving of dx, dy variables for verification. Optional (in sequential only)
   IF (l_bathy_save) THEN

#ifdef key_MPI_2D
     WRITE(ierrorlog,*) 'ERROR: routine GRDIMETRICS.F90'
     WRITE(ierrorlog,*) 'Check the grid (dx,dy,dx_u...) in sequential only'
     CALL_MPI MPI_FINALIZE(IERR_MPI) ; STOP
#endif

     write(1231,*) 'j,MINVAL(lon2d(i,ljmin:ljmax)),MAXVAL(lon2d(i,ljmin:ljmax))'
     DO i = imin,imax
       write(1231,*) j,MINVAL(lon2d(i,ljmin:ljmax)),MAXVAL(lon2d(i,ljmin:ljmax))
     END DO
     write(1231,*) 'j,MINVAL(lon2du(i,ljmin:ljmax)),MAXVAL(lon2du(i,ljmin:ljmax))'
     DO i = imin,imax
       write(1231,*) j,MINVAL(lon2du(i,ljmin:ljmax)),MAXVAL(lon2du(i,ljmin:ljmax))
     END DO
     write(1231,*) 'j,MINVAL(lon2dv(i,ljmin:ljmax)),MAXVAL(lon2dv(i,ljmin:ljmax))'
     DO i = imin,imax
       write(1231,*) j,MINVAL(lon2dv(i,ljmin:ljmax)),MAXVAL(lon2dv(i,ljmin:ljmax))
     END DO

     write(1232,*) 'j,MINVAL(lat2d(limin:limax,j)),MAXVAL(lat2d(limin:limax,j))'
     DO j = jmin,jmax
       write(1232,*) j,MINVAL(lat2d(limin:limax,j)),MAXVAL(lat2d(limin:limax,j))
     END DO
     write(1232,*) 'j,MINVAL(lat2du(limin:limax,j)),MAXVAL(lat2du(limin:limax,j))'
     DO j = jmin,jmax
       write(1232,*) j,MINVAL(lat2du(limin:limax,j)),MAXVAL(lat2du(limin:limax,j))
     END DO
     write(1232,*) 'j,MINVAL(lat2dv(limin:limax,j)),MAXVAL(lat2dv(limin:limax,j))'
     DO j = jmin,jmax
       write(1232,*) j,MINVAL(lat2dv(limin:limax,j)),MAXVAL(lat2dv(limin:limax,j))
     END DO

     write(1233,*) 'j,MINVAL(dx(limin:limax,j)),MAXVAL(dx(limin:limax,j))'
     DO j = jmin,jmax
       write(1233,*) j,MINVAL(dx(limin:limax,j)),MAXVAL(dx(limin:limax,j))
     END DO
     write(1233,*) 'j,MINVAL(dx_u(limin:limax,j)),MAXVAL(dx_u(limin:limax,j))'
     DO j = jmin,jmax
       write(1233,*) j,MINVAL(dx_u(limin:limax,j)),MAXVAL(dx_u(limin:limax,j))
     END DO
     write(1233,*) 'j,MINVAL(dx_v(limin:limax,j)),MAXVAL(dx_v(limin:limax,j))'
     DO j = jmin,jmax
       write(1233,*) j,MINVAL(dx_v(limin:limax,j)),MAXVAL(dx_v(limin:limax,j))
     END DO
     write(1233,*) 'j,MINVAL(dx_f(limin:limax,j)),MAXVAL(dx_f(limin:limax,j))'
     DO j = jmin,jmax
       write(1233,*) j,MINVAL(dx_f(limin:limax,j)),MAXVAL(dx_f(limin:limax,j))
     END DO

     write(1234,*) 'dy etc = cst or not ?'
     write(1234,*) 'dy',MINVAL(dy(:,:)),MAXVAL(dy(:,:))
     write(1234,*) 'dy_u',MINVAL(dy_u(:,:)),MAXVAL(dy_u(:,:))
     write(1234,*) 'dy_v',MINVAL(dy_v(:,:)),MAXVAL(dy_v(:,:))
     write(1234,*) 'dy_f',MINVAL(dy_f(:,:)),MAXVAL(dy_f(:,:))
     write(1234,*) 'j,MINVAL(dy(limin:limax,j)),MAXVAL(dy(limin:limax,j))'
     DO j = jmin,jmax
       write(1234,*) j,MINVAL(dy(limin:limax,j)),MAXVAL(dy(limin:limax,j))
     END DO
     write(1234,*) 'j,MINVAL(dy_u(limin:limax,j)),MAXVAL(dy_u(limin:limax,j))'
     DO j = jmin,jmax
       write(1234,*) j,MINVAL(dy_u(limin:limax,j)),MAXVAL(dy_u(limin:limax,j))
     END DO
     write(1234,*) 'j,MINVAL(dy_v(limin:limax,j)),MAXVAL(dy_v(limin:limax,j))'
     DO j = jmin,jmax
       write(1234,*) j,MINVAL(dy_v(limin:limax,j)),MAXVAL(dy_v(limin:limax,j))
     END DO
     write(1234,*) 'j,MINVAL(dy_f(limin:limax,j)),MAXVAL(dy_f(limin:limax,j))'
     DO j = jmin,jmax
       write(1234,*) j,MINVAL(dy_f(limin:limax,j)),MAXVAL(dy_f(limin:limax,j))
     END DO


     IF_MPI (MASTER) THEN
     CALL ionc4_open('bathy_verif.nc')
     CALL ionc4_createvar('bathy_verif.nc','dx','degres','dx', &
                          dims='xy',fill_value=dg_valmanq_io)
     tmp2drsh(:,:)=dx(limin:limax,ljmin:ljmax)
     CALL ionc4_write_xy('bathy_verif.nc','dx',tmp2drsh,imin,imax,jmin,jmax)
     CALL ionc4_createvar('bathy_verif.nc','dx_u','degres','dx_u', &
                          dims='xy',fill_value=dg_valmanq_io)
     tmp2drsh(:,:)=dx_u(limin:limax,ljmin:ljmax)
     CALL ionc4_write_xy('bathy_verif.nc','dx_u',tmp2drsh,imin,imax,jmin,jmax)
     CALL ionc4_createvar('bathy_verif.nc','dx_v','degres','dx_v', &
                          dims='xy',fill_value=dg_valmanq_io)
     tmp2drsh(:,:)=dx_v(limin:limax,ljmin:ljmax)
     CALL ionc4_write_xy('bathy_verif.nc','dx_v',tmp2drsh,imin,imax,jmin,jmax)
     CALL ionc4_createvar('bathy_verif.nc','dx_f','degres','dx_f', &
                          dims='xy',fill_value=dg_valmanq_io)
     tmp2drsh(:,:)=dx_u(limin:limax,ljmin:ljmax)
     CALL ionc4_write_xy('bathy_verif.nc','dx_f',tmp2drsh,imin,imax,jmin,jmax)

     CALL ionc4_createvar('bathy_verif.nc','dy','degres','dy', &
                          dims='xy',fill_value=dg_valmanq_io)
     tmp2drsh(:,:)=dy(limin:limax,ljmin:ljmax)
     CALL ionc4_write_xy('bathy_verif.nc','dy',tmp2drsh,imin,imax,jmin,jmax)
     CALL ionc4_createvar('bathy_verif.nc','dy_u','degres','dy_u', &
                          dims='xy',fill_value=dg_valmanq_io)
     tmp2drsh(:,:)=dy_u(limin:limax,ljmin:ljmax)
     CALL ionc4_write_xy('bathy_verif.nc','dy_u',tmp2drsh,imin,imax,jmin,jmax)
     CALL ionc4_createvar('bathy_verif.nc','dy_v','degres','dy_v', &
                          dims='xy',fill_value=dg_valmanq_io)
     tmp2drsh(:,:)=dy_v(limin:limax,ljmin:ljmax)
     CALL ionc4_write_xy('bathy_verif.nc','dy_v',tmp2drsh,imin,imax,jmin,jmax)
     CALL ionc4_createvar('bathy_verif.nc','dy_f','degres','dy_f', &
                          dims='xy',fill_value=dg_valmanq_io)
     tmp2drsh(:,:)=dy_v(limin:limax,ljmin:ljmax)
     CALL ionc4_write_xy('bathy_verif.nc','dy_f',tmp2drsh,imin,imax,jmin,jmax)
     CALL ionc4_close('bathy_verif.nc')
     ENDIF_MPI
   ENDIF
#endif

   PRINT_DBG*, 'exit GRID_METRIC routine'
   
  END SUBROUTINE grid_metric
