MODULE compute_geometry_mod
   USE geometry
   USE dimensions

   USE domain_mod, ONLY: t_domain, &
                         domain, ndomain, assigned_domain, &
                         domain_glo, ndomain_glo, domloc_glo_ind, swap_needed
   USE omp_para, ONLY: is_omp_level_master, is_master
   USE transfert_mod, ONLY: req_i0, req_i1, t_message, transfert_request, transfert_message, init_message

   USE spherical_geom_mod, ONLY: xyz2lonlat, circumcenter, compute_centroid, &
                                 surf_triangle, dist_cart, div_arc_bis, schmidt_transform
   USE vector, ONLY: norm, cross_product2

   IMPLICIT NONE
   PRIVATE
   SAVE

   PUBLIC :: compute_geometry

CONTAINS

   SUBROUTINE compute_geometry
      USE grid_param, ONLY: pmesh, grid_ico, grid_unst, grid_type, max_primal_deg, max_dual_deg, max_trisk_deg
      USE data_unstructured_mod
      USE compute_geometry_ico_mod, ONLY: set_geometry, update_domain
      USE set_bounds_mod, ONLY: set_bounds
      USE read_pmesh_mod, ONLY: angle_to_vector, pmesh_to_cellsets
      INTEGER :: ij

      CALL allocate_geometry(max_primal_deg, max_dual_deg, max_trisk_deg)

      SELECT CASE (grid_type)
      CASE (grid_ico)
         CALL set_geometry
         CALL set_bounds(domain_glo, .TRUE.)
         CALL set_bounds(domain, .FALSE.)
      CASE (grid_unst)
         swap_needed = .TRUE.
         CALL swap_geometry(1)
         swap_needed = .FALSE.
         ! after swap_geometry, XX is a pointer to data allocated in geom%XX by init_geometry
         ! copy data from pmesh%XX to XX
         Ai(:) = pmesh%Ai(:)
         Av(:) = pmesh%Av(:)
         fv(:) = pmesh%fv(:)
         le_de(:) = pmesh%le_de(:)
         le(:) = pmesh%le(:)
         de(:) = pmesh%de(:)
         lon_e(:) = pmesh%lon_e(:)
         lat_e(:) = pmesh%lat_e(:)
         lon_v(:) = pmesh%lon_v(:)
         lat_v(:) = pmesh%lat_v(:)
         R_iv2(:, :) = pmesh%Riv2(:, :)
         wee(:, :) = pmesh%wee(:, :)
         grad_ii = pmesh%grad(:, :, :)
         perot_ie = pmesh%perot(:, :, :)

         ! Workaround for apparent compiler bug in ifort 15
         DO ij = 1, SIZE(lon_i)
            lon_i(ij) = pmesh%lon_i(ij)
         END DO
         DO ij = 1, SIZE(lat_i)
            lat_i(ij) = pmesh%lat_i(ij)
         END DO
         IF (MAXVAL(ABS(lon_i(1:10) - pmesh%lon_i(1:10))) > 1e-6) THEN
            PRINT *, "** geom%XX    **", lon_i(1:3)
            PRINT *, "** pmesh%XX   **", pmesh%lon_i(1:3)
            PRINT *, "** geom-pmesh mismatch **"
            STOP
         END IF
         IF (MAXVAL(ABS(lat_i(1:10) - pmesh%lat_i(1:10))) > 1e-6) THEN
            PRINT *, "** geom%XX    **", lat_i(1:3)
            PRINT *, "** pmesh%XX   **", pmesh%lat_i(1:3)
            PRINT *, "** geom-pmesh mismatch **"
            STOP
         END IF
         ! primal_deg, etc. are POINTERS so we must either associate (=>) or allocate them
         ! indices in mesh file follow Fortran convention (first index is 1)
         primal_deg => pmesh%primal_deg
         primal_ne => pmesh%primal_ne
         primal_edge => pmesh%primal_edge
         primal_vertex => pmesh%primal_vertex
         dual_deg => pmesh%dual_deg
         dual_edge => pmesh%dual_edge
         dual_ne => pmesh%dual_ne
         dual_vertex => pmesh%dual_vertex
         left => pmesh%left
         right => pmesh%right
         up => pmesh%up
         down => pmesh%down
         trisk_deg => pmesh%trisk_deg
         trisk => pmesh%trisk
         CALL angle_to_vector(lon_e, lat_e, pmesh%angle_e, ep_e)
         CALL pmesh_to_cellsets(pmesh, &
                                domain(1)%primal_own, domain(1)%dual_own, domain(1)%edge_own)
         CALL pmesh_to_cellsets(pmesh, &
                                domain_glo(1)%primal_own, domain_glo(1)%dual_own, domain_glo(1)%edge_own)
      CASE DEFAULT
         STOP 'Invalid value of grid_type encountered in compute_geometry'
      END SELECT
   END SUBROUTINE compute_geometry

END MODULE compute_geometry_mod
