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_v(:)
      TYPE(t_field), POINTER :: lon_e(:)
      TYPE(t_field), POINTER :: lat_i(:)
      TYPE(t_field), POINTER :: lat_v(:)
      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 :: fv(:)
      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 :: ne(:)
      TYPE(t_field), POINTER :: bi(:)
      TYPE(t_field), POINTER :: S1(:)
      TYPE(t_field), POINTER :: S2(:)
      ! Weights Riv and Riv2 are for averaging to dual cells (triangles) from primal cells (hexagons).
      ! On the icosahedral mesh it is possible to arrange these weights as 6 weights per hexagon
      ! but on the unstructured mesh they must be arranged as 3 weights per triangle.
      ! => to move forward we define both Riv2 (icosahedral) and R_iv2 (unstructured)
      ! but we should change the ico-mesh implementation to match the unst-mesh one.
      TYPE(t_field), POINTER :: Riv(:)
      TYPE(t_field), POINTER :: Riv2(:)
      TYPE(t_field), POINTER :: R_iv2(:)
      ! Weights wee are for averaging from edges to edges (TRiSK perp operator).
      ! On the icosahedral mesh there are 10=2x5 such weights
      ! => array wee with shape (:,10).
      ! On the unstructured mesh there are (a-1)+(b-1) weights
      ! with a,b the number of edges of the two cells that share the edge
      ! => array wee with shape (:,max_trisk_deg)
      TYPE(t_field), POINTER :: wee(:)
      TYPE(t_field), POINTER :: perot_ie(:)
      TYPE(t_field), POINTER :: grad_ii(:)
   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_v(:)       ! longitude of the center of the dual cell (triangle)
!$OMP THREADPRIVATE(lon_v)
   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_v(:)       ! latitude of the center of the dual cell (triangle)
!$OMP THREADPRIVATE(lat_v)
   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)
   INTEGER, POINTER    :: ne(:, :)        ! convention for the way on the normal wind on an edge
!$OMP THREADPRIVATE(ne)
   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(:, :)       ! weights
!$OMP THREADPRIVATE(Riv)
   REAL(rstd), POINTER :: Riv2(:, :)      ! weights, shape (t,6)
!$OMP THREADPRIVATE(Riv2)
   REAL(rstd), POINTER :: R_iv2(:, :)     ! same weights, shape (z, 3)
!$OMP THREADPRIVATE(R_iv2)
   REAL(rstd), POINTER :: wee(:, :)      ! same weights, shape (e,10)
!$OMP THREADPRIVATE(wee)
   REAL(rstd), POINTER :: fv(:)          ! coriolis (evaluated on a vertex)
!$OMP THREADPRIVATE(fv)
   REAL(rstd), POINTER :: perot_ie(:, :, :)  ! coefficients of Perot vector reconstruction operator
!$OMP THREADPRIVATE(perot_ie)
   REAL(rstd), POINTER :: grad_ii(:, :, :)  ! coefficients of Green-Gauss gradient operator
!$OMP THREADPRIVATE(grad_ii)
   REAL(rstd), POINTER :: bi(:)          ! orographie
!$OMP THREADPRIVATE(bi)

   ! FIXME : these definitions are ico-specific, they belong elsewhere
   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(nvertex, ndual_vertex, ntrisk)
      USE field_mod
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: nvertex, ndual_vertex, ntrisk

      CALL allocate_field(geom%Ai, field_t, type_real, name='geom%Ai')
      CALL allocate_field(geom%xyz_i, field_t, type_real, 3)
      CALL allocate_field(geom%lon_i, field_t, type_real, name='geom%lon_i')
      CALL allocate_field(geom%lat_i, field_t, type_real, name='geom%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='geom%lon_e')
      CALL allocate_field(geom%lat_e, field_u, type_real, name='geom%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='geom%de')
      CALL allocate_field(geom%le, field_u, type_real, name='geom%le')
      CALL allocate_field(geom%le_de, field_u, type_real, name='geom%le_de')
      CALL allocate_field(geom%bi, field_t, type_real)

      CALL allocate_field(geom%Av, field_z, type_real, name='geom%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, ntrisk, name="geom%wee")
      CALL allocate_field(geom%fv, field_z, type_real, name='geom%fv')

      CALL allocate_field(geom%lon_v, field_z, type_real, name='lon_v')
      CALL allocate_field(geom%lat_v, field_z, type_real, name='lat_v')
      CALL allocate_field(geom%R_iv2, field_z, type_real, ndual_vertex, name='geom%R_iv2')
      CALL allocate_field(geom%perot_ie, field_t, type_real, nvertex, 3)
      CALL allocate_field(geom%grad_ii, field_t, type_real, nvertex, 3)

   END SUBROUTINE allocate_geometry

   SUBROUTINE swap_geometry(ind)
      USE field_mod
      USE domain_mod, ONLY: swap_needed
      INTEGER, INTENT(IN) :: ind

      IF (.NOT. swap_needed) RETURN

      Ai = geom%Ai(ind)
      lon_i = geom%lon_i(ind)
      lat_i = geom%lat_i(ind)
      ne = geom%ne(ind)
      xyz_i = geom%xyz_i(ind)
      elon_i = geom%elon_i(ind)
      elat_i = geom%elat_i(ind)
      centroid = geom%centroid(ind)
      grad_ii = geom%grad_ii(ind)
      perot_ie = geom%perot_ie(ind)
      xyz_e = geom%xyz_e(ind)
      elon_e = geom%elon_e(ind)
      elat_e = geom%elat_e(ind)
      ep_e = geom%ep_e(ind)
      et_e = geom%et_e(ind)

      de = geom%de(ind)
      le = geom%le(ind)
      le_de = geom%le_de(ind)
      lon_e = geom%lon_e(ind)
      lat_e = geom%lat_e(ind)
      wee = geom%wee(ind)

      Av = geom%Av(ind)
      fv = geom%fv(ind)
      lon_v = geom%lon_v(ind)
      lat_v = geom%lat_v(ind)
      xyz_v = geom%xyz_v(ind)

      S1 = geom%S1(ind)
      S2 = geom%S2(ind)
      Riv = geom%Riv(ind)
      Riv2 = geom%Riv2(ind)
      R_iv2 = geom%R_iv2(ind)
      bi = geom%bi(ind)
   END SUBROUTINE swap_geometry

END MODULE geometry
