! Original file is in jsrc
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
    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
    SELECT CASE(grid_type)
    CASE(grid_ico)
       CALL set_geometry
       ! copy position of generators and vertices back into domain_glo(:)%xyz/vertex
       ! so that write_field gets the right values
       CALL update_domain
       CALL set_bounds(domain_glo, .TRUE.)
       CALL set_bounds(domain, .FALSE.)
    CASE(grid_unst)
       swap_needed = .TRUE.
       CALL swap_geometry(1)
       swap_needed = .FALSE.
       ! now XX points to data allocated to geom%XX by init_geometry
       ! copy data from pmesh%XX to XX
       {%- for XX in 'Ai', 'Av', 'fv', 'le_de', 'le', 'de',
       'lon_e', 'lat_e', 'lon_v', 'lat_v' %}
       {{XX}}(:) = pmesh%{{XX}}(:)
       {%- endfor %}
       R_iv2(:,:) = pmesh%Riv2(:,:)
       wee(:,:)  = pmesh%wee(:,:)
       ! Workaround for apparent compiler bug in ifort 15
       {%- for XX in 'lon_i', 'lat_i' %}
       DO ij=1, SIZE({{XX}})
          {{XX}}(ij) = pmesh%{{XX}}(ij)
       END DO
       {%- endfor %}
       {%- for XX in 'lon_i', 'lat_i' %}
       IF(MAXVAL(ABS({{XX}}(1:10)-pmesh%{{XX}}(1:10)))>1e-6) THEN
          PRINT *, "** geom%XX    **", {{XX}}(1:3)
          PRINT *, "** pmesh%XX   **", pmesh%{{XX}}(1:3)
          PRINT *, "** geom-pmesh mismatch **"
          STOP
       END IF
       {%- endfor %}
       ! primal_deg, etc. are POINTERS so we must either associate (=>) or allocate them
       ! indices in mesh file follow Fortran convention (first index is 1)
       {%- for member in
           'primal_deg', 'primal_ne', 'primal_edge', 'primal_vertex',
           'dual_deg', 'dual_edge', 'dual_ne', 'dual_vertex',
           'left', 'right', 'up', 'down', 'trisk_deg', 'trisk' %}
       {{member}} => pmesh%{{member}}
       {%- endfor %}
       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
