MODULE geometry
  USE field_mod
  IMPLICIT NONE

  TYPE t_geometry
    TYPE(t_field),POINTER :: centroid(:)
    TYPE(t_field),POINTER :: xyz_i(:)
    TYPE(t_field),POINTER :: xyz_e(:)
    TYPE(t_field),POINTER :: xyz_v(:)
    TYPE(t_field),POINTER :: lon_i(:)
    TYPE(t_field),POINTER :: lon_e(:)
    TYPE(t_field),POINTER :: lat_i(:)
    TYPE(t_field),POINTER :: lat_e(:)
    TYPE(t_field),POINTER :: ep_e(:)
    TYPE(t_field),POINTER :: et_e(:)
    TYPE(t_field),POINTER :: elon_i(:)
    TYPE(t_field),POINTER :: elat_i(:)
    TYPE(t_field),POINTER :: elon_e(:)
    TYPE(t_field),POINTER :: elat_e(:)
    TYPE(t_field),POINTER :: Ai(:)
    TYPE(t_field),POINTER :: Av(:)
    TYPE(t_field),POINTER :: de(:)
    TYPE(t_field),POINTER :: le(:)
    TYPE(t_field),POINTER :: le_de(:) ! le/de, 0. if de=0.
    TYPE(t_field),POINTER :: Riv(:)
    TYPE(t_field),POINTER :: S1(:)
    TYPE(t_field),POINTER :: S2(:)
    TYPE(t_field),POINTER :: Riv2(:)
    TYPE(t_field),POINTER :: ne(:)
    TYPE(t_field),POINTER :: Wee(:)
    TYPE(t_field),POINTER :: bi(:)
    TYPE(t_field),POINTER :: fv(:)
  END TYPE t_geometry    

  TYPE(t_geometry),SAVE,TARGET :: geom

  
  REAL(rstd),POINTER :: Ai(:)          ! area of a cell
!$OMP THREADPRIVATE(Ai)
  REAL(rstd),POINTER :: centroid(:,:)  ! coordinate of the centroid of the cell
!$OMP THREADPRIVATE(centroid)
  REAL(rstd),POINTER :: xyz_i(:,:)     ! coordinate of the center of the cell (voronoi)
!$OMP THREADPRIVATE(xyz_i)
  REAL(rstd),POINTER :: xyz_e(:,:)     ! coordinate of a wind point on the cell on a edge
!$OMP THREADPRIVATE(xyz_e)
  REAL(rstd),POINTER :: xyz_v(:,:)     ! coordinate of a vertex (center of the dual mesh)
!$OMP THREADPRIVATE(xyz_v)
  REAL(rstd),POINTER :: lon_i(:)       ! longitude of the center of the cell (voronoi)
!$OMP THREADPRIVATE(lon_i)
  REAL(rstd),POINTER :: lon_e(:)       ! longitude of a wind point on the cell on a edge
!$OMP THREADPRIVATE(lon_e)
  REAL(rstd),POINTER :: lat_i(:)       ! latitude of the center of the cell (voronoi)
!$OMP THREADPRIVATE(lat_i)
  REAL(rstd),POINTER :: lat_e(:)       ! latitude of a wind point on the cell on a edge
!$OMP THREADPRIVATE(lat_e)
  REAL(rstd),POINTER :: ep_e(:,:)      ! perpendicular unit vector of a edge (outsider)
!$OMP THREADPRIVATE(ep_e)
  REAL(rstd),POINTER :: et_e(:,:)      ! tangeantial unit vector of a edge
!$OMP THREADPRIVATE(et_e)
  REAL(rstd),POINTER :: elon_i(:,:)    ! unit longitude vector on the center 
!$OMP THREADPRIVATE(elon_i)
  REAL(rstd),POINTER :: elat_i(:,:)    ! unit latitude vector on the center 
!$OMP THREADPRIVATE(elat_i)
  REAL(rstd),POINTER :: elon_e(:,:)    ! unit longitude vector on a wind point 
!$OMP THREADPRIVATE(elon_e)
  REAL(rstd),POINTER :: elat_e(:,:)    ! unit latitude vector on a wind point
!$OMP THREADPRIVATE(elat_e)
  REAL(rstd),POINTER :: Av(:)          ! area of dual mesk cell
!$OMP THREADPRIVATE(Av)
  REAL(rstd),POINTER :: de(:)          ! distance from a neighbour == lenght of an edge of the dual mesh
!$OMP THREADPRIVATE(de)
  REAL(rstd),POINTER :: le(:)          ! lenght of a edge
!$OMP THREADPRIVATE(le)
  REAL(rstd),POINTER :: le_de(:)       ! le/de
!$OMP THREADPRIVATE(le_de)
  REAL(rstd),POINTER :: S1(:,:)        ! area of sub-triangle
!$OMP THREADPRIVATE(S1)
  REAL(rstd),POINTER :: S2(:,:)        ! area of sub-tirangle
!$OMP THREADPRIVATE(S2)
  REAL(rstd),POINTER :: Riv(:,:)       ! weight
!$OMP THREADPRIVATE(Riv)
  REAL(rstd),POINTER :: Riv2(:,:)      ! weight
!$OMP THREADPRIVATE(Riv2)
  INTEGER,POINTER    :: ne(:,:)        ! convention for the way on the normal wind on an edge 
!$OMP THREADPRIVATE(ne)
  REAL(rstd),POINTER :: Wee(:,:,:)     ! weight
!$OMP THREADPRIVATE(Wee)
  REAL(rstd),POINTER :: bi(:)          ! orographie
!$OMP THREADPRIVATE(bi)
  REAL(rstd),POINTER :: fv(:)          ! coriolis (evaluted on a vertex)
!$OMP THREADPRIVATE(fv)
     

  INTEGER, PARAMETER :: ne_right=1
  INTEGER, PARAMETER :: ne_rup=-1
  INTEGER, PARAMETER :: ne_lup=1
  INTEGER, PARAMETER :: ne_left=-1
  INTEGER, PARAMETER :: ne_ldown=1
  INTEGER, PARAMETER :: ne_rdown=-1

CONTAINS

  SUBROUTINE allocate_geometry
  USE field_mod
  IMPLICIT NONE
  INTEGER, PARAMETER :: nvertex=6    ! FIXME unstructured

    CALL allocate_field(geom%Ai,field_t,type_real,name='Ai')

    CALL allocate_field(geom%xyz_i,field_t,type_real,3)
    CALL allocate_field(geom%lon_i,field_t,type_real, name='lon_i')
    CALL allocate_field(geom%lat_i,field_t,type_real, name='lat_i')
    CALL allocate_field(geom%elon_i,field_t,type_real,3)
    CALL allocate_field(geom%elat_i,field_t,type_real,3)
    CALL allocate_field(geom%centroid,field_t,type_real,3)

    CALL allocate_field(geom%xyz_e,field_u,type_real,3)
    CALL allocate_field(geom%lon_e,field_u,type_real, name='lon_e')
    CALL allocate_field(geom%lat_e,field_u,type_real, name='lat_e')
    CALL allocate_field(geom%elon_e,field_u,type_real,3)
    CALL allocate_field(geom%elat_e,field_u,type_real,3)
    CALL allocate_field(geom%ep_e,field_u,type_real,3)
    CALL allocate_field(geom%et_e,field_u,type_real,3)

    CALL allocate_field(geom%xyz_v,field_z,type_real,3)
    CALL allocate_field(geom%de,field_u,type_real, name='de')
    CALL allocate_field(geom%le,field_u,type_real, name='le')
    CALL allocate_field(geom%le_de,field_u,type_real, name='le_de')
    CALL allocate_field(geom%bi,field_t,type_real)
    CALL allocate_field(geom%Av,field_z,type_real, name='Av')
    CALL allocate_field(geom%S1,field_t,type_real,nvertex)
    CALL allocate_field(geom%S2,field_t,type_real,nvertex)
    CALL allocate_field(geom%Riv,field_t,type_real,nvertex)
    CALL allocate_field(geom%Riv2,field_t,type_real,nvertex)
    CALL allocate_field(geom%ne,field_t,type_integer,nvertex)
    CALL allocate_field(geom%Wee,field_u,type_real,5,2) ! FIXME unstructured
    CALL allocate_field(geom%bi,field_t,type_real)
    CALL allocate_field(geom%fv,field_z,type_real, name='fv')

  END SUBROUTINE allocate_geometry
  
  
  SUBROUTINE swap_geometry(ind)
  USE field_mod
  USE domain_mod, ONLY : swap_needed
  IMPLICIT NONE
    INTEGER,INTENT(IN) :: ind

    IF(.NOT. swap_needed) RETURN

!!$OMP MASTER
    Ai=geom%Ai(ind)
    xyz_i=geom%xyz_i(ind)
    centroid=geom%centroid(ind)
    xyz_e=geom%xyz_e(ind)
    ep_e=geom%ep_e(ind)
    et_e=geom%et_e(ind)
    lon_i=geom%lon_i(ind)
    lat_i=geom%lat_i(ind)
    lon_e=geom%lon_e(ind)
    lat_e=geom%lat_e(ind)
    elon_i=geom%elon_i(ind)
    elat_i=geom%elat_i(ind)
    elon_e=geom%elon_e(ind)
    elat_e=geom%elat_e(ind)
    xyz_v=geom%xyz_v(ind)
    de=geom%de(ind)
    le=geom%le(ind)
    le_de=geom%le_de(ind)
    Av=geom%Av(ind)
    S1=geom%S1(ind)
    S2=geom%S2(ind)
    Riv=geom%Riv(ind)
    Riv2=geom%Riv2(ind)
    ne=geom%ne(ind)
    Wee=geom%Wee(ind)
    bi=geom%bi(ind)
    fv=geom%fv(ind)
!!$OMP END MASTER
!!$OMP BARRIER    
  END SUBROUTINE swap_geometry

END MODULE geometry