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

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





























































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


  USE comvars2d, ONLY : dx_g,dy_g,dxu_g,dyu_g,dxv_g,dyv_g,dxf_g,dyf_g


   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




   USE toolmpi,   ONLY : MASTER,ex_i_rsh,ex_j_rsh

   IMPLICIT NONE

!! * Arguments

!! * Local declarations
   INTEGER        :: i,j,k,ixl,iyl,nb

   REAL(KIND=rlg) :: rhox,rhoy
   REAL(KIND=rsh) :: fi,g,radeg






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

   !PRINT*, 'enter GRIDMETRIC routine'



!$OMP PARALLEL DEFAULT(none) PRIVATE(i,j,k,radeg,g,           &

!$OMP& fi)                                                    &
!$OMP& SHARED(dx,dy,govlet,goulet,cosg,sing,cos2g,sin2g       &
!$OMP& ,lonmin,latmin,lonmax,latmax                           &


!$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                       &

!$OMP& ,MASTER                                                &

!$OMP& ,fv,tfivsr,surf,sig,dsigu,l_modele2d,iwarnlog          &
!$OMP& ,l_equisig,dsigw,sigw,fisac                            &





!$OMP ,h0,dx_g,dy_g,dxu_g,dyu_g,dxv_g,dyv_g,dxf_g,dyf_g          &

!$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)



       
! 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




! 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



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


     END DO
   END DO
!$OMP END DO





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

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

         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

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

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

!$OMP single

   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



!$OMP END single



!$OMP END PARALLEL





!! Verification du critere d inconsistence hydrostatique !!





   !PRINT*, 'exit GRID_METRIC routine'
   
  END SUBROUTINE grid_metric
