LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - mesh_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 534 | 0.0 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 0 | 29 | 0.0 %

              Line data    Source code
           1             : module mesh_mod
           2             :
           3             :   use shr_kind_mod,   only: r8=>shr_kind_r8
           4             :   use physconst,      only: PI
           5             :   use control_mod,    only: MAX_FILE_LEN
           6             :   use cam_abortutils, only: endrun
           7             :
           8             :   use netcdf,         only: nf90_strerror, nf90_open, nf90_close
           9             :   use netcdf,         only: NF90_NOWRITE, nf90_NoErr
          10             :   use netcdf,         only: nf90_inq_dimid, nf90_inquire_dimension
          11             :   use netcdf,         only: nf90_inq_varid, nf90_get_var
          12             :
          13             :   implicit none
          14             :   logical, public           :: MeshUseMeshFile = .false.
          15             :
          16             :   public  :: MeshOpen           ! Must be called first
          17             :
          18             :   integer, parameter :: MXSTLN = 32
          19             :
          20             :   ! ===============================
          21             :   ! Public methods for mesh_mod
          22             :   ! ===============================
          23             :
          24             :   public  :: MeshCubeEdgeCount  ! called anytime afer MeshOpen
          25             :   public  :: MeshCubeElemCount  ! called anytime afer MeshOpen
          26             :   public  :: MeshCubeTopology   ! called afer MeshOpen
          27             :   public  :: MeshSetCoordinates ! called after MeshCubeTopology
          28             :   public  :: MeshPrint          ! show the contents of the Mesh after it has been loaded into the module
          29             :   public  :: MeshClose
          30             :   ! ===============================
          31             :   ! Private members
          32             :   ! ===============================
          33             :
          34             :   integer,private,parameter :: nfaces          = 6 ! number of faces on the cube
          35             :   integer,private,parameter :: nInnerElemEdge  = 8 ! number of edges for an interior element
          36             :
          37             :   character (len=MAX_FILE_LEN), private              :: p_mesh_file_name
          38             :   integer                     , private              :: p_ncid
          39             :   integer                     , private              :: p_number_elements
          40             :   integer                     , private              :: p_number_elements_per_face
          41             :   integer                     , private              :: p_number_blocks
          42             :   integer                     , private              :: p_number_nodes
          43             :   integer                     , private              :: p_number_dimensions
          44             :   integer                     , private              :: p_number_neighbor_edges
          45             :    real(kind=r8)       , private, allocatable :: p_node_coordinates(:,:)
          46             :   integer                     , private, allocatable :: p_connectivity(:,:)
          47             :
          48             :   ! ===============================
          49             :   ! Private methods
          50             :   ! ===============================
          51             :
          52             :   private :: create_index_table
          53             :   private :: find_side_neighbors
          54             :   private :: find_corner_neighbors
          55             :   private :: get_node_coordinates
          56             :   private :: get_2D_sub_coordinate_indexes
          57             :   private :: mesh_connectivity
          58             :   private :: cube_face_element_centroids
          59             :   private :: smallest_diameter_element
          60             :   private :: cube_to_cube_coordinates
          61             :   private :: sphere_to_cube_coordinates
          62             :   private :: initialize_space_filling_curve
          63             :
          64             :   private :: handle_error
          65             :   private :: open_mesh_file
          66             :   private :: close_mesh_file
          67             :   private :: get_number_of_elements
          68             :   private :: get_number_of_dimensions
          69             :   private :: get_number_of_elements_per_face
          70             :   private :: get_number_of_nodes
          71             :   private :: get_number_of_element_blocks
          72             :   private :: get_node_multiplicity
          73             :   private :: get_face_connectivity
          74             :
          75             : CONTAINS
          76             :
          77             : !======================================================================
          78             : !  subroutine handle_error
          79             : !======================================================================
          80           0 :   subroutine handle_error (status, file, line)
          81             :
          82             :     integer,            intent(in) :: status
          83             :     character (len=*),  intent(in) :: file
          84             :     integer,            intent(in) :: line
          85           0 :     print *, file,':', line, ': ', trim(nf90_strerror(status))
          86           0 :     call endrun("Terminating program due to netcdf error while obtaining mesh information, please see message in standard output.")
          87           0 :   end subroutine handle_error
          88             :
          89             : !======================================================================
          90             : !  open_mesh_file()
          91             : !
          92             : !> Open the netcdf file containing the mesh.
          93             : !! Assign the holder to the file to p_ncid so everyone else knows
          94             : !! how to use it without passing the argument around.
          95             : !======================================================================
          96           0 :   subroutine open_mesh_file()
          97             :     implicit none
          98             :     integer                        :: status
          99             :
         100           0 :     status = nf90_open(p_mesh_file_name, NF90_NOWRITE, p_ncid)
         101           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         102             :
         103           0 :     MeshUseMeshFile = .true.
         104             :
         105           0 :   end subroutine open_mesh_file
         106             :
         107             : !======================================================================
         108             : ! subroutine close_mesh_file()
         109             : !======================================================================
         110             :
         111           0 :   subroutine close_mesh_file()
         112             :     implicit none
         113             :     integer              :: status
         114             :
         115           0 :     status = nf90_close(p_ncid)
         116           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         117             :
         118           0 :   end subroutine close_mesh_file
         119             :
         120             : !======================================================================
         121             : ! function get_number_of_dimensions()
         122             : !======================================================================
         123             :
         124           0 :   function get_number_of_dimensions() result(number_dimensions)
         125             :     implicit none
         126             :     integer              :: number_dimensions
         127             :
         128             :      ! local variables
         129             :      integer              :: status, number_of_dim_id
         130             :
         131             :      ! Get the id of 'num_elem', if such dimension is not there panic and quit :P
         132           0 :     status = nf90_inq_dimid(p_ncid, "num_dim", number_of_dim_id)
         133           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         134             :
         135             :     ! How many values for 'num_elem' are there?
         136           0 :     status = nf90_inquire_dimension(p_ncid, number_of_dim_id, len = number_dimensions)
         137           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         138             :
         139           0 :   end function get_number_of_dimensions
         140             :
         141             : !======================================================================
         142             : ! function get_number_of_elements()
         143             : !======================================================================
         144             :
         145           0 :   function get_number_of_elements() result(number_elements)
         146             :     implicit none
         147             :     integer              :: number_elements
         148             :     ! local variables
         149             :     integer              :: status, number_of_elements_id
         150             :
         151             :     ! Get the id of 'num_elem', if such dimension is not there panic and quit :P
         152           0 :     status = nf90_inq_dimid(p_ncid, "num_elem", number_of_elements_id)
         153           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         154             :
         155             :     ! How many values for 'num_elem' are there?
         156           0 :     status = nf90_inquire_dimension(p_ncid, number_of_elements_id, len = number_elements)
         157           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         158             :
         159           0 :   end function get_number_of_elements
         160             :
         161             : !======================================================================
         162             : !  function get_number_of_nodes()
         163             : !======================================================================
         164           0 :   function get_number_of_nodes() result(number_nodes)
         165             :     implicit none
         166             :     integer              :: number_nodes
         167             :     ! local variables
         168             :     integer              :: status, number_of_nodes_id
         169             :
         170             :     ! Get the id of 'num_nodes', if such dimension is not there panic and quit :P
         171           0 :     status = nf90_inq_dimid(p_ncid, "num_nodes", number_of_nodes_id)
         172           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         173             :
         174             :     ! How many values for 'num_nodes' are there?
         175           0 :     status = nf90_inquire_dimension(p_ncid, number_of_nodes_id, len = number_nodes)
         176           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         177             :
         178           0 :   end function get_number_of_nodes
         179             :
         180             :
         181             : !======================================================================
         182             : !  function get_number_of_element_blocks()
         183             : !======================================================================
         184           0 :   function get_number_of_element_blocks() result(number_element_blocks)
         185             :
         186             :     integer              :: number_element_blocks
         187             :     ! local variables
         188             :     integer              :: status, number_of_element_blocks_id
         189             :
         190             :     ! Get the id of 'num_el_blk', if such dimension is not there panic and quit :P
         191           0 :     status = nf90_inq_dimid(p_ncid, "num_el_blk", number_of_element_blocks_id)
         192           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         193             :
         194             :     ! How many values for 'num_el_blk' are there?
         195           0 :     status = nf90_inquire_dimension(p_ncid, number_of_element_blocks_id, len = number_element_blocks)
         196           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         197             :
         198           0 :      if (number_element_blocks /= 1) then
         199           0 :         if (number_element_blocks /= 6  ) then
         200           0 :             call endrun('Reading cube-sphere from input file is not supported')
         201             :         else
         202           0 :            call endrun('Number of elements blocks not exactly 1 (sphere) or 6 (cube)')
         203             :         endif
         204             :      endif
         205             :
         206           0 :   end function get_number_of_element_blocks
         207             :
         208             : !======================================================================
         209             : !  function get_number_of_elements_per_face()
         210             : !======================================================================
         211           0 :     function get_number_of_elements_per_face() result(number_elements_per_face)
         212             :
         213             :     integer             :: number_elements_per_face
         214             :
         215             :     integer               :: face_num ! For each of the face, we get the information
         216             :     character(len=MXSTLN) :: element_type ! Each face is composed of elements of certain type
         217             :     integer               :: number_elements_in_face ! How many elements in this face
         218             :     integer               :: num_nodes_per_elem ! How many nodes in each element
         219             :     integer               :: number_of_attributes ! How many attributes in the face
         220             :
         221             :     integer               :: status, dimension_id
         222             :
         223           0 :     if (p_number_blocks == 0)  then
         224           0 :        call endrun('get_number_of_elements_per_face called before MeshOpen')
         225           0 :     else if (p_number_blocks == 1) then ! we are in the presence of a sphere
         226             :        ! First we get sure the number of nodes per element is four
         227           0 :        status = nf90_inq_dimid(p_ncid, "num_nod_per_el1", dimension_id)
         228           0 :        if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         229           0 :        status = nf90_inquire_dimension(p_ncid, dimension_id, len =  num_nodes_per_elem)
         230           0 :        if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         231           0 :        if (num_nodes_per_elem /= 4)  call endrun('Number of nodes per element is not four')
         232             :        ! now we check how many elements there are in the face
         233           0 :        status = nf90_inq_dimid(p_ncid, "num_el_in_blk1", dimension_id)
         234           0 :        if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         235           0 :        status = nf90_inquire_dimension(p_ncid, dimension_id, len = number_elements_in_face)
         236           0 :        if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         237           0 :        number_elements_per_face =  number_elements_in_face
         238           0 :     else if (p_number_blocks == 6) then ! we are in the presence of a cube-sphere
         239           0 :        call endrun('Reading a mesh for a cube-sphere is not supported')
         240             :     else
         241           0 :        call endrun('Number of elements blocks not exactly 1 (sphere) or 6 (cube)')
         242             :     end if
         243             :
         244           0 :   end function get_number_of_elements_per_face
         245             :
         246             : !======================================================================
         247             : ! subroutine get_face_connectivity
         248             : !======================================================================
         249           0 :   subroutine get_face_connectivity()
         250             :
         251             :     integer              :: var_id, status
         252             :
         253           0 :     status = nf90_inq_varid(p_ncid, "connect1", var_id)
         254           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         255           0 :     status = nf90_get_var(p_ncid, var_id, p_connectivity)
         256           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         257           0 :   end subroutine get_face_connectivity
         258             :
         259             : !======================================================================
         260             : ! subroutine get_node_multiplicity
         261             : !======================================================================
         262           0 :   subroutine get_node_multiplicity(node_multiplicity)
         263             :     use dimensions_mod, only : max_elements_attached_to_node
         264             :
         265             :     integer, intent(out) :: node_multiplicity(:)
         266             :     integer              :: node_num(4)
         267             :
         268             :     integer              :: k, number_nodes
         269             :
         270           0 :     node_multiplicity(:) = 0
         271           0 :     number_nodes = SIZE(node_multiplicity)
         272             :     ! check this external buffer was allocated correctly
         273           0 :     if (number_nodes /= p_number_nodes) call endrun('Number of nodes does not matches size of node multiplicity array')
         274             :     ! for each node, we have for four other nodes
         275             :
         276           0 :     if (minval(p_connectivity) < 1 .or. number_nodes < maxval(p_connectivity)) then
         277           0 :        call endrun('get_node_multiplicity: Node number less than 1 or greater than max.')
         278             :     end if
         279             :
         280           0 :     do k=1,p_number_elements_per_face
         281           0 :        node_num = p_connectivity(:,k)
         282           0 :        node_multiplicity(node_num) = node_multiplicity(node_num) + 1
         283             :     enddo
         284             :
         285           0 :     if (minval(node_multiplicity) < 3 .or. max_elements_attached_to_node < maxval(node_multiplicity)) then
         286           0 :       print *, 'minval(node_multiplicity)', minval(node_multiplicity)
         287           0 :       print *, 'maxval(node_multiplicity)', maxval(node_multiplicity),&
         288           0 :            ' and max_elements_attached_to_node ',max_elements_attached_to_node
         289           0 :       call endrun('get_node_multiplicity: Number of elements attached to node less than 3 or greater than maximum.')
         290             :     endif
         291             :
         292           0 :   end subroutine get_node_multiplicity
         293             :
         294             : !======================================================================
         295             : !  subroutine get_node_coordinates ()
         296             : !======================================================================
         297           0 :   subroutine get_node_coordinates ()
         298             :
         299             :     integer              :: var_id, status
         300             :
         301           0 :     status = nf90_inq_varid(p_ncid, "coord", var_id)
         302           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         303           0 :     status = nf90_get_var(p_ncid, var_id, p_node_coordinates)
         304           0 :     if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
         305           0 :   end subroutine get_node_coordinates
         306             :
         307             :   ! ================================================================================
         308             :   !
         309             :   ! -----------------Internal private routines that do not use netCDF IO -----------
         310             :   !
         311             :   ! ================================================================================
         312             :
         313             : !======================================================================
         314             : ! subroutine get_2D_sub_coordinate_indexes
         315             : !======================================================================
         316           0 :   subroutine get_2D_sub_coordinate_indexes(x, y, sgnx, sgny, face_no)
         317             :      implicit none
         318             :     integer, intent(in)              :: face_no
         319             :     integer, intent(out)             :: x,y
         320             :     integer, intent(out)             :: sgnx, sgny
         321           0 :     if (face_no == 1 .or. face_no == 3) then
         322           0 :        x = 2
         323           0 :        y = 3
         324           0 :     else if (face_no == 2 .or. face_no == 4) then
         325           0 :        x = 1
         326           0 :        y = 3
         327             :     else
         328           0 :        x = 2
         329           0 :        y = 1
         330             :     endif
         331           0 :     if (face_no == 1 .or. face_no == 4 .or. face_no == 5) then
         332           0 :        sgnx =  1
         333           0 :        sgny =  1
         334           0 :     else if (face_no == 2 .or. face_no == 3) then
         335           0 :        sgnx = -1
         336           0 :        sgny =  1
         337             :     else
         338           0 :        sgnx =  1
         339           0 :        sgny = -1
         340             :     endif
         341           0 :   end subroutine get_2D_sub_coordinate_indexes
         342             :
         343             :
         344             :
         345             : !======================================================================
         346             : ! subroutine mesh_connectivity(connect)
         347             : !
         348             : ! puts the transpose of p_connectivity into connect
         349             : !======================================================================
         350             :
         351           0 :   subroutine  mesh_connectivity (connect)
         352             :
         353             :     integer,  intent(out) :: connect(p_number_elements,4)
         354             :
         355             :     integer :: k, j
         356             :
         357           0 :     if (0 == p_number_blocks)  call endrun('mesh_connectivity called before MeshOpen')
         358           0 :     j=0
         359           0 :     do k=1, p_number_elements_per_face
         360           0 :        j=j+1
         361           0 :        connect(j,:) = p_connectivity(:,k)
         362             :     enddo
         363             :
         364           0 :     if (j /= p_number_elements) call endrun('mesh_connectivity: Number of elements in side sets not equal to total elements')
         365             :
         366           0 :     if (minval(connect) < 1 .or. maxval(connect) > p_number_nodes) then
         367           0 :        call endrun('mesh_connectivity: Node number out of bounds')
         368             :     end if
         369             :
         370           0 :   end subroutine mesh_connectivity
         371             : !======================================================================
         372             : ! subroutine create_index_table()
         373             : !
         374             : ! this is needed to detremine side and corner neighbors
         375             : !======================================================================
         376             :
         377           0 :   subroutine create_index_table(index_table, element_nodes)
         378             :     use dimensions_mod, only : max_elements_attached_to_node
         379             :
         380             :     integer, allocatable, intent(inout)  :: index_table(:,:)
         381             :     integer             ,  intent(in)    :: element_nodes(p_number_elements, 4)
         382             :     integer                              :: cnt, cnt_index, node
         383             :     integer                              :: k, ll
         384             :
         385             :     !Create an index table so that we can find neighbors on O(n)
         386             :     ! so for each node, we want to know which elements it is part of
         387           0 :     allocate(index_table(p_number_nodes, max_elements_attached_to_node + 1))
         388             :
         389             :     !the last column in the index table is a count of the number of elements
         390           0 :     index_table = 0
         391             :
         392           0 :     cnt_index =  max_elements_attached_to_node + 1
         393             :
         394           0 :     do k=1,p_number_elements
         395           0 :         do ll=1,4
         396           0 :            node = element_nodes(k, ll) !the node
         397           0 :            cnt = index_table(node, cnt_index)  !how many elements for that node already in table
         398           0 :            cnt = cnt + 1 !increment since we are adding an element
         399           0 :            if (cnt >  max_elements_attached_to_node) then
         400           0 :               call endrun('Found a node in too many elements.')
         401             :            endif
         402           0 :            index_table(node, cnt_index) = cnt
         403           0 :            index_table(node, cnt) = k !put the element in the indextable
         404             :         enddo
         405             :     enddo
         406             :
         407           0 :   end subroutine create_index_table
         408             :
         409             : !======================================================================
         410             : ! subroutine find_side_neighbors()
         411             : !
         412             : ! find the element neighbors to the n,s,e,w and put them in GridVertex_t
         413             : !  (only 1 neighbor to the n,s,e,w)
         414             : !======================================================================
         415           0 :   subroutine find_side_neighbors (GridVertex, normal_to_homme_ordering, element_nodes, edge_wgt, index_table)
         416             :     use coordinate_systems_mod, only : cartesian3D_t
         417             :     use gridgraph_mod, only   : GridVertex_t
         418             :     use dimensions_mod, only : max_elements_attached_to_node
         419             :
         420             :     integer             ,  intent(in)    :: normal_to_homme_ordering(8)
         421             :     integer             ,  intent(in)    :: element_nodes(p_number_elements, 4)
         422             :     integer             ,  intent(in)    :: edge_wgt
         423             :     integer             ,  intent(in)     :: index_table(:,:)
         424             :     type (GridVertex_t) ,  intent(inout) :: GridVertex(:)
         425             :
         426             :     integer                              :: i_node(2), my_node(2)
         427             :     integer                              :: neighbor, direction
         428             :     integer                              :: j,k,ll,i, m
         429             :     integer                              :: i_elem, jump, end_i
         430             :     integer                              :: loc, cnt_index, a_count(2)
         431             :     logical                              :: found
         432           0 :     if (0 == p_number_blocks)  call endrun('find_side_neighbors called before MeshOpen')
         433             :
         434             :
         435             :     !the last column in the index table is a count of the number of elements
         436           0 :     cnt_index =  max_elements_attached_to_node + 1
         437             :
         438             :     !use index table to find neighbors
         439           0 :     do k=1,p_number_elements  ! for each element k
         440             :        !set the side weights
         441           0 :        GridVertex(k)%nbrs_wgt(1:4) = edge_wgt
         442           0 :        do ll=1,4 !loop through the four sides
         443             :
         444           0 :           jump = normal_to_homme_ordering(ll)
         445           0 :           loc = GridVertex(k)%nbrs_ptr(jump)
         446             :
         447           0 :           if (GridVertex(k)%nbrs(loc) == 0) then  !if side is not set yet, then
         448             :              !look for side element
         449           0 :              found = .false.
         450           0 :              neighbor = 0
         451             :
         452           0 :              my_node(1) = element_nodes(k, ll)
         453           0 :              a_count(1) = index_table(my_node(1), cnt_index)
         454           0 :              my_node(2) = element_nodes(k, mod(ll,4)+1)
         455           0 :              a_count(2) = index_table(my_node(2), cnt_index)
         456             :
         457             :              !loop through the elements that are in the index table for each node
         458             :              !and find the element number and direction of the side neighbor
         459           0 :              do m = 1,2
         460           0 :                 if (found) exit
         461           0 :                 end_i = a_count(m)
         462           0 :                 do i = 1, end_i
         463           0 :                    if (found) exit
         464           0 :                    i_elem = index_table(my_node(m),i)
         465           0 :                    if (i_elem /= k) then !k is the element we are setting sides for
         466           0 :                       do j=1,4 !loop through each of i_elem's four sides
         467           0 :                          i_node(1) = element_nodes(i_elem, j)
         468           0 :                          i_node(2) = element_nodes(i_elem, mod(j,4)+1)
         469           0 :                          if ( (i_node(1) == my_node(2) .and. i_node(2) == my_node(1)) .or. &
         470           0 :                               (i_node(1) == my_node(1) .and. i_node(2) == my_node(2)) ) then
         471             :                             neighbor = i_elem
         472             :                             direction = j
         473             :                             found = .true.
         474             :                             !found a match
         475             :                             exit
         476             :                          end if
         477             :                       end do ! j loop
         478             :                    end if
         479             :                 enddo ! i loop
         480             :              enddo !m loop
         481             :
         482           0 :              if (neighbor == 0) call endrun('find_side_neighbor: Neighbor not found! Every side should have a neighbor.')
         483             :
         484           0 :              GridVertex(k)%nbrs(loc) = neighbor
         485           0 :              jump = normal_to_homme_ordering(direction)
         486           0 :              loc = GridVertex(neighbor)%nbrs_ptr(jump)
         487           0 :              GridVertex(neighbor)%nbrs(loc)= k
         488             :           endif
         489             :        enddo !  ll loop => 4 sides
         490             :     enddo ! k loop: each element
         491             :
         492           0 :     do k=1,p_number_elements
         493           0 :        do ll=1,4
         494           0 :           if ( 0 == GridVertex(k)%nbrs(ll)) then
         495           0 :              call endrun('Found one side of one element witout a neighbor.  Bummer!')
         496             :           end if
         497             :        end do
         498             :     end do
         499             :
         500           0 :   end subroutine find_side_neighbors
         501             :
         502             : !======================================================================
         503             : ! function smallest_diameter_element
         504             : !======================================================================
         505             :
         506           0 :   function smallest_diameter_element(element_nodes) result(min_diameter)
         507             :
         508             :     integer             ,intent(in)  :: element_nodes(:,:)
         509             :
         510             :     integer                          :: i, j
         511             :     integer                          :: node_numbers(4)
         512             :     real(kind=r8)             :: coordinates (4,3)
         513             :     real(kind=r8)             :: x(3), y(3), r(3), d, min_diameter
         514             :
         515           0 :     if (SIZE(element_nodes,dim=1) /= p_number_elements) then
         516             :        call endrun('smallest_diameter_element:Element count check failed in &
         517           0 :             &exodus_mesh. Connectivity array length not equal to number of elements.')
         518             :     end if
         519           0 :     if ( p_number_elements_per_face /= p_number_elements) then
         520             :        call endrun('smallest_diameter_element: Element count check failed in &
         521           0 :             &exodus_mesh. Element array length not equal to sum of face.')
         522             :     end if
         523             :
         524           0 :     min_diameter = 9999999.0_r8
         525           0 :     do i=1, p_number_elements
         526           0 :        node_numbers = element_nodes(i,:)
         527           0 :        coordinates = p_node_coordinates(node_numbers,:)
         528             :        ! smallest side length
         529           0 :        do j=1,4
         530           0 :           x = coordinates(j         ,:)
         531           0 :           y = coordinates(1+MOD(j,4),:)
         532           0 :           r = x-y
         533           0 :           d  = dot_product(r,r)
         534           0 :           if (d < min_diameter ) then
         535           0 :              min_diameter = d
         536             :           end if
         537             :        end do
         538             :        ! smallest diameter length
         539           0 :        do j=1,2
         540           0 :           x = coordinates(j         ,:)
         541           0 :           y = coordinates(2+MOD(j,4),:)
         542           0 :           r = x-y
         543           0 :           d  = dot_product(r,r)
         544           0 :           if (d < min_diameter ) then
         545           0 :              min_diameter = d
         546             :           end if
         547             :        end do
         548             :     enddo
         549           0 :     min_diameter = SQRT(min_diameter)
         550           0 :   end function smallest_diameter_element
         551             :
         552             : !======================================================================
         553             : !  subroutine cube_to_cube_coordinates
         554             : !======================================================================
         555             :
         556           0 :   subroutine cube_to_cube_coordinates (cube_coor, node_coor, face_number)
         557             :
         558             :     real(kind=r8), intent(in)  :: node_coor(4,3)
         559             :     integer,       intent(in)  :: face_number
         560             :     real(kind=r8), intent(out) :: cube_coor(4,2)
         561             :
         562             :     integer                              :: x_index, y_index, sgnx, sgny
         563           0 :     call get_2D_sub_coordinate_indexes(x_index, y_index, sgnx, sgny, face_number)
         564           0 :     cube_coor(:,1) = sgnx*node_coor(:,x_index)
         565           0 :     cube_coor(:,2) = sgny*node_coor(:,y_index)
         566           0 :   end subroutine cube_to_cube_coordinates
         567             :
         568             :
         569             : !======================================================================
         570             : !  subroutine sphere_to_cube_coordinates
         571             : !======================================================================
         572             :
         573           0 :   subroutine sphere_to_cube_coordinates (cube_coor, node_coor, face_number)
         574             :     use coordinate_systems_mod, only   : cartesian2d_t, change_coordinates, sphere2cubedsphere
         575             :     implicit none
         576             :     real(kind=r8),    intent(in)   :: node_coor(4,3)
         577             :     integer,                 intent(in)   :: face_number
         578             :     real(kind=r8),    intent(out)  :: cube_coor(4,2)
         579             :     integer                               :: i
         580             :     type(cartesian2d_t)                   :: cart(4)
         581             :
         582           0 :     do i=1,4
         583           0 :        cart(i) = sphere2cubedsphere(change_coordinates(node_coor(i,:)), face_number)
         584             :     end do
         585           0 :     cube_coor(:,1) = cart(:)%x
         586           0 :     cube_coor(:,2) = cart(:)%y
         587           0 :   end subroutine sphere_to_cube_coordinates
         588             :
         589             :
         590             : !======================================================================
         591             : !  subroutine cube_face_element_centroids
         592             : !======================================================================
         593             :
         594           0 :   subroutine cube_face_element_centroids(centroids, face_numbers, element_nodes)
         595             :
         596             :     integer            , intent(in)  :: element_nodes(:,:)
         597             :     integer,             intent(in)  :: face_numbers    (p_number_elements)
         598             :     real(kind=r8),intent(out) :: centroids       (p_number_elements,2)
         599             :     real(kind=r8)             :: coordinates(4,3)
         600             :     real(kind=r8)             :: cube_coor  (4,2)
         601             :     integer                          :: i, node_numbers(4)
         602             :
         603           0 :     if (0 == p_number_blocks)  call endrun('cube_face_element_centroids called before MeshOpen')
         604           0 :     if (SIZE(element_nodes,dim=1) /= p_number_elements) then
         605             :        call endrun('cube_face_element_centroids:Element count check failed in &
         606           0 :             &exodus_mesh. Connectivity array length not equal to number of elements.')
         607             :     end if
         608           0 :     if ( p_number_elements_per_face /= p_number_elements ) then
         609             :        call endrun('cube_face_element_centroids: Element count check failed in &
         610           0 :             &exodus_mesh. Element array length not equal to sum of face.')
         611             :     end if
         612             :
         613           0 :     do i=1, p_number_elements
         614           0 :        node_numbers = element_nodes(i,:)
         615           0 :        coordinates = p_node_coordinates(node_numbers,:)
         616           0 :        if (6 == p_number_blocks) then
         617           0 :           call cube_to_cube_coordinates   (cube_coor, coordinates, face_numbers(i))
         618             :        else
         619           0 :           call sphere_to_cube_coordinates (cube_coor, coordinates, face_numbers(i))
         620             :        end if
         621           0 :        centroids(i,:) = SUM(cube_coor,dim=1)/4.0_r8
         622             :     enddo
         623           0 :   end subroutine cube_face_element_centroids
         624             :
         625             : !======================================================================
         626             : ! subroutine initialize_space_filling_curve
         627             : !======================================================================
         628           0 :   subroutine initialize_space_filling_curve(GridVertex, element_nodes)
         629             :     use gridgraph_mod, only   : GridVertex_t
         630             :     use spacecurve_mod, only  : GenspaceCurve
         631             :
         632             :     type (GridVertex_t), intent(inout) :: GridVertex(:)
         633             :     integer            , intent(in)    :: element_nodes(:,:)
         634             :
         635           0 :     integer,allocatable                :: Mesh2(:,:),Mesh2_map(:,:),sfcij(:,:)
         636             :
         637           0 :     real(kind=r8)               :: centroids(p_number_elements,2)
         638           0 :     integer                            :: face_numbers(p_number_elements)
         639             :     real(kind=r8)               :: x, y, h
         640             :     integer                            :: i, j, i2, j2, ne, ne2
         641             :     integer                            :: sfc_index, face
         642             :
         643           0 :     if (SIZE(GridVertex) /= p_number_elements) then
         644             :        call endrun('initialize_space_filling_curve:Element count check failed &
         645           0 :             &in exodus_mesh. Vertex array length not equal to number of elements.')
         646             :     end if
         647           0 :     if (SIZE(element_nodes,dim=1) /= p_number_elements) then
         648             :        call endrun('initialize_space_filling_curve:Element count check failed &
         649           0 :             &in exodus_mesh. Connectivity array length not equal to number of elements.')
         650             :     end if
         651             :
         652           0 :     face_numbers(:) = GridVertex(:)%face_number
         653           0 :     h = smallest_diameter_element    (                         element_nodes)
         654             :
         655           0 :     call cube_face_element_centroids (centroids, face_numbers, element_nodes)
         656             :
         657           0 :     if (h<.00001_r8) then
         658           0 :       call endrun('initialize_space_filling_curve: Unreasonably small element found. less than .00001')
         659             :     end if
         660             :
         661           0 :     ne = CEILING(0.5_r8*PI/(h/2));
         662             :
         663             :     ! find the smallest ne2 which is a power of 2 and ne2>ne
         664           0 :     ne2=2**ceiling( log(real(ne))/log(2._r8) )
         665           0 :     if (ne2<ne) call endrun('initialize_space_filling_curve: Fatel SFC error')
         666             :
         667           0 :     allocate(Mesh2(ne2,ne2))
         668           0 :     allocate(Mesh2_map(ne2,ne2))
         669           0 :     allocate(sfcij(0:ne2*ne2,2))
         670             :
         671             :     ! create a reverse index array for Mesh2
         672             :     ! j = Mesh2(i,j)
         673             :     ! (i,j) = (sfcij(j,1),sfci(j,2))
         674           0 :     call GenspaceCurve(Mesh2)  ! SFC partition for ne2
         675           0 :     do j2=1,ne2
         676           0 :        do i2=1,ne2
         677           0 :           j=Mesh2(i2,j2)
         678           0 :           sfcij(j,1)=i2
         679           0 :           sfcij(j,2)=j2
         680             :        enddo
         681             :     enddo
         682             :
         683             :
         684           0 :     GridVertex(:)%SpaceCurve=-1
         685             :     sfc_index   = 0
         686           0 :     do face = 1,nfaces
         687             :        ! associate every element on the ne x ne mesh (Mesh)
         688             :        ! with its closest element on the ne2 x ne2 mesh (Mesh2)
         689             :        ! Store this as a map from Mesh2 -> Mesh in Mesh2_map.
         690             :        ! elements in Mesh2 which are not mapped get assigned a value of 0
         691           0 :        Mesh2_map=0
         692           0 :        do i=1,p_number_elements
         693           0 :           if (face_numbers(i) == face ) then
         694           0 :              x = centroids(i,1)
         695           0 :              y = centroids(i,2)
         696             :              ! map this element to an (i2,j2) element
         697             :              ! [ -PI/4, PI/4 ]  -> [ 0, ne2 ]
         698           0 :              i2=nint( (0.5_r8 + 2.0_r8*x/PI)*ne2 + 0.5_r8 )
         699           0 :              j2=nint( (0.5_r8 + 2.0_r8*y/PI)*ne2 + 0.5_r8 )
         700           0 :              if (face == 4 .or. face == 6 )               i2 = ne2-i2+1
         701           0 :              if (face == 1 .or. face == 2 .or. face == 6) j2 = ne2-j2+1
         702           0 :              if (i2<1  ) i2=1
         703           0 :              if (i2>ne2) i2=ne2
         704           0 :              if (j2<1  ) j2=1
         705           0 :              if (j2>ne2) j2=ne2
         706           0 :              Mesh2_map(i2,j2)=i
         707             :           end if
         708             :        end do
         709             :
         710             :        ! generate a SFC for Mesh with the same ordering as the
         711             :        ! elements in Mesh2 which map to Mesh.
         712           0 :        do j=0,ne2*ne2-1
         713           0 :           i2=sfcij(j,1)
         714           0 :           j2=sfcij(j,2)
         715           0 :           i=Mesh2_map(i2,j2)
         716           0 :           if (i/=0) then
         717             :              ! (i2,j2) element maps to element
         718           0 :              GridVertex(i)%SpaceCurve=sfc_index
         719           0 :              sfc_index=sfc_index+1
         720             :           endif
         721             :        enddo
         722             :     enddo
         723           0 :     deallocate(Mesh2)
         724           0 :     deallocate(Mesh2_map)
         725           0 :     deallocate(sfcij)
         726             :
         727           0 :     if (minval(GridVertex(:)%SpaceCurve) == -1) then
         728           0 :        do i=1,p_number_elements
         729           0 :           if (-1==GridVertex(i)%SpaceCurve) then
         730           0 :              write (*,*) " Error in projecting element ",i," to space filling curve."
         731           0 :              write (*,*) " Face:",face_numbers(i)
         732           0 :              write (*,*) " Centroid:",centroids(i,:)
         733             :           end if
         734             :        end do
         735           0 :        call endrun('initialize_space_filling_curve: Vertex not on SpaceCurve')
         736             :     end if
         737             :
         738           0 :   end subroutine initialize_space_filling_curve
         739             :
         740             : !======================================================================
         741             : ! subroutine find_corner_neighbors
         742             : !======================================================================
         743             :
         744           0 :   subroutine find_corner_neighbors  (GridVertex, normal_to_homme_ordering, element_nodes, corner_wgt, index_table)
         745             :     use gridgraph_mod,          only : GridVertex_t
         746             :     use dimensions_mod,         only : max_elements_attached_to_node, max_corner_elem
         747             :     use control_mod, only: north, south, east, west, neast,seast, nwest,swest
         748             :
         749             :     type (GridVertex_t), intent(inout) :: GridVertex(:)
         750             :     integer            , intent(in)    :: normal_to_homme_ordering(8)
         751             :     integer            , intent(in)    :: element_nodes(p_number_elements, 4)
         752             :     integer            , intent(in)    :: corner_wgt
         753             :     integer            , intent(in)    :: index_table(:,:)
         754             :
         755           0 :     integer                          :: node_elements (2*max_elements_attached_to_node)
         756           0 :     integer                          :: elem_neighbor (4*max_elements_attached_to_node)
         757             :     integer                          :: nbr_cnt(4)
         758             :     integer                          :: elem_nbr_start, start
         759             :     integer                          :: i, j, k, ll, jj, kk
         760             :     integer                          :: node, loc, cnt, cnt_index
         761           0 :     integer                          :: corner_array(max_corner_elem), orig_pos(max_corner_elem)
         762           0 :     integer                          :: face_array(max_corner_elem), a_corner_elems(max_corner_elem)
         763             :     integer                          :: corner_sides(2)
         764             :     integer                          :: side_elem, corner_elem, tmp_s
         765             :
         766             :     !the last column in the index table is a count of the number of elements
         767           0 :     cnt_index =  max_elements_attached_to_node + 1
         768             :
         769           0 :     do i=1, p_number_elements  !loop through all elements
         770           0 :        node_elements(:) = 0
         771           0 :        elem_neighbor(:) = 0
         772           0 :        elem_nbr_start = 0
         773           0 :        nbr_cnt(:) = 0
         774             :
         775           0 :        do j=1,4  !check each of the 4 nodes at the element corners
         776           0 :           node = element_nodes(i,j)
         777           0 :           cnt = index_table(node, cnt_index)
         778           0 :           if (cnt < 3 .or. max_elements_attached_to_node < cnt) then
         779           0 :              call endrun('find_corner_neighbors: Number of elements attached to node less than 3 or greater than maximum.')
         780             :           endif
         781             :
         782           0 :           node_elements(1:cnt) = index_table(node, 1:cnt)
         783             :
         784             :           !now node_elements contains the element neighbors to that node - so grab the
         785             :           ! corner neighbors - these are the ones that are not already side neighbors (or myself)
         786           0 :           k = 0
         787           0 :           do ll=1,cnt
         788           0 :              if ( i /= node_elements(ll) .and. & !not me
         789           0 :                   GridVertex(i)%nbrs(1) /= node_elements(ll) .and. & !not side 1
         790           0 :                   GridVertex(i)%nbrs(2) /= node_elements(ll) .and. & ! etc ...
         791           0 :                   GridVertex(i)%nbrs(3) /= node_elements(ll) .and. &
         792           0 :                   GridVertex(i)%nbrs(4) /= node_elements(ll)) then
         793           0 :                 k = k + 1
         794           0 :                 elem_neighbor(elem_nbr_start + k) = node_elements(ll)
         795             :              end if
         796             :           end do ! end of ll loop for multiplicity
         797             :
         798             :           !keep track of where we are starting in elem_neighbor for each corner j
         799           0 :           elem_nbr_start = elem_nbr_start + k
         800           0 :           nbr_cnt(j) = k !how many neighbors in this corner
         801             :        end do ! end of j loop through 4 nodes
         802             :
         803             :
         804             :        ! now that we have done the 4 corners we can populate nbrs and nbrs_ptr
         805             :        ! with the corners in the proper order (clockwise) in neighbors
         806             :        ! also we can add the corner weight
         807             :
         808           0 :        do j=5,8 !loop through 4 corners
         809             :           elem_nbr_start = 1
         810             :           !easiest to do the corner in ascending order - find loc
         811           0 :           do jj = 5,8
         812           0 :              ll = normal_to_homme_ordering(jj)
         813           0 :              if (j == ll) then
         814             :                 loc = jj
         815             :                 exit
         816             :              end if
         817           0 :              elem_nbr_start = elem_nbr_start + nbr_cnt(jj-4)
         818             :           end do
         819             :
         820           0 :           start =  GridVertex(i)%nbrs_ptr(j)
         821           0 :           cnt = nbr_cnt(loc - 4)
         822           0 :           GridVertex(i)%nbrs_ptr(j+1) = start + cnt
         823             :
         824           0 :           if (cnt > 0) then
         825           0 :              GridVertex(i)%nbrs(start : start + cnt-1) = &
         826           0 :                   elem_neighbor(elem_nbr_start : elem_nbr_start + cnt -1)
         827           0 :              GridVertex(i)%nbrs_face(start : start + cnt - 1) = &
         828           0 :                   GridVertex(elem_neighbor(elem_nbr_start : elem_nbr_start + cnt -1))%face_number
         829           0 :              GridVertex(i)%nbrs_wgt(start : start + cnt-1)  = corner_wgt
         830             :
         831             :           end if
         832             :
         833             :           ! within each corner neighbor, lets list the corners in clockwise order
         834           0 :           if (cnt > 1) then !cnt is the number of neighbors in this corner j
         835             :                             !there can be at most max_corner element of these
         836             :
         837           0 :              a_corner_elems = 0
         838           0 :              a_corner_elems(1:cnt) = elem_neighbor(elem_nbr_start : elem_nbr_start + cnt -1)
         839             :              !corner-sides(2) is clockwise of corner_side(1)
         840           0 :              corner_array= 0
         841           0 :              orig_pos = 0
         842           0 :              select case (j)
         843             :              case(neast)
         844           0 :                 corner_sides(1) = north
         845           0 :                 corner_sides(2) = east
         846             :              case(seast)
         847           0 :                 corner_sides(1) = east
         848           0 :                 corner_sides(2) = south
         849             :              case(swest)
         850           0 :                 corner_sides(1) = south
         851           0 :                 corner_sides(2) = west
         852             :              case(nwest)
         853           0 :                 corner_sides(1) = west
         854           0 :                 corner_sides(2) = north
         855             :              end select
         856             :
         857             :              !so the first element to list touches  corner_sides(1) element
         858           0 :              side_elem = GridVertex(i)%nbrs(corner_sides(1))
         859             :
         860             :              !loop though the corner elements and see if any have a side neighbor
         861             :              !that = side_elem
         862           0 :              do k = 1,cnt !number of corner elements
         863           0 :                 corner_elem = a_corner_elems(k)
         864           0 :                 do kk = 1,4 !number of sides to check
         865           0 :                    loc = GridVertex(corner_elem)%nbrs_ptr(kk)
         866           0 :                    tmp_s = GridVertex(corner_elem)%nbrs(loc)
         867           0 :                    if (tmp_s == side_elem) then
         868           0 :                       corner_array(1) = corner_elem
         869           0 :                       orig_pos(1) = k
         870           0 :                       exit
         871             :                    endif
         872             :                 enddo
         873           0 :                 if  (corner_array(1)> 0) exit
         874             :              enddo
         875           0 :              if (corner_array(1)==0) then
         876           0 :                 print *, i, cnt
         877           0 :                 call endrun('find_corner_neighbors (1) : mistake finding corner neighbor order')
         878             :              endif
         879             :
         880             :              !if cnt == 2, we are done (we know the order of neighbors)
         881           0 :              if (cnt ==2) then
         882           0 :                 if (corner_array(1) ==  a_corner_elems(1)) then
         883           0 :                    corner_array(2) =  a_corner_elems(2)
         884           0 :                    orig_pos(2) = 2
         885             :                 else
         886           0 :                    corner_array(2) =  a_corner_elems(1)
         887           0 :                    orig_pos(2) = 1
         888             :                 end if
         889             :              else !cnt = 3 or 4
         890             :                 !find which corner element borders corner_sides(2)
         891           0 :                 side_elem = GridVertex(i)%nbrs(corner_sides(2))
         892           0 :                 do k = 1,cnt
         893           0 :                    corner_elem = a_corner_elems(k)
         894           0 :                    do kk = 1,4
         895           0 :                       loc = GridVertex(corner_elem)%nbrs_ptr(kk)
         896           0 :                       tmp_s = GridVertex(corner_elem)%nbrs(loc)
         897           0 :                       if (tmp_s == side_elem) then
         898           0 :                          corner_array(4) = corner_elem
         899           0 :                          orig_pos(4) = k
         900           0 :                          exit
         901             :                       endif
         902             :                    enddo
         903           0 :                    if  (corner_array(4)> 0) exit
         904             :                 enddo
         905           0 :                 if (corner_array(4)==0 .or. corner_array(4) == corner_array(1)) then
         906           0 :                    print *, i, cnt
         907           0 :                    call endrun('find_corner_neighbors (2) : mistake finding corner neighbor order')
         908             :                 endif
         909             :
         910             :                 !now if cnt = 3 then we are done
         911           0 :                 if (cnt ==3) then
         912           0 :                    corner_array(3) = corner_array(4)
         913           0 :                    orig_pos(3) = orig_pos(4)
         914             :
         915           0 :                    do k = 1,cnt !find the "middle" element
         916           0 :                       if (k /= orig_pos(1) .and. k /= orig_pos(3)) then
         917           0 :                          orig_pos(2) = k
         918           0 :                          corner_array(2) = a_corner_elems(k)
         919           0 :                          exit
         920             :                       endif
         921             :                    enddo
         922             :                 else  !cnt = 4
         923             :                    !which of the two unassigned elements borders the element in
         924             :                    !corner_array(1) => put in corner_array(2)
         925           0 :                    side_elem = corner_array(1)
         926             :
         927           0 :                    do k = 1,cnt
         928           0 :                       corner_elem = a_corner_elems(k)
         929           0 :                       if (corner_elem == corner_array(4) .or. corner_elem == corner_array(1)) then
         930             :                          cycle
         931             :                       else
         932           0 :                          do kk = 1,4 !check each side
         933           0 :                             loc = GridVertex(corner_elem)%nbrs_ptr(kk)
         934           0 :                             tmp_s = GridVertex(corner_elem)%nbrs(loc)
         935           0 :                             if (tmp_s == side_elem) then
         936           0 :                                corner_array(2) = corner_elem
         937           0 :                                orig_pos(2) = k
         938           0 :                                exit
         939             :                             endif
         940             :                          enddo
         941             :                       endif
         942           0 :                       if  (corner_array(2)> 0) exit
         943             :                    enddo
         944             :                    !now put the remaining one in pos 3
         945           0 :                    do k = 1,cnt
         946           0 :                       corner_elem = a_corner_elems(k)
         947             :                       if (corner_elem /= corner_array(4) .and. corner_elem /= &
         948           0 :                            corner_array(2) .and. corner_elem /= corner_array(1)) then
         949           0 :                          corner_array(3) = corner_elem
         950           0 :                          orig_pos(3) = k
         951           0 :                          exit
         952             :                       endif
         953             :                    enddo
         954             :                 endif ! end of cnt=4
         955             :              endif! end of not cnt=2
         956             :
         957             :              !now re-set the elements in this corner
         958           0 :              GridVertex(i)%nbrs(start : start + cnt-1) = corner_array(1:cnt)
         959             :              !nbrs_wgt are the same - nothing to do
         960             :              !fix neighbors face
         961           0 :              do k = 1,cnt
         962           0 :                 face_array(k) = GridVertex(i)%nbrs_face(start + orig_pos(k) - 1)
         963             :              end do
         964           0 :              GridVertex(i)%nbrs_face(start : start + cnt - 1) = face_array(1:cnt)
         965             :           endif !end of cnt > 1 loop for corners
         966             :
         967             :        end do !j loop through each corner
         968             :
         969             :     end do ! end of i loop through elements
         970           0 :   end subroutine find_corner_neighbors
         971             :
         972             :   ! ================================================================================
         973             :   !
         974             :   ! -------------------------------Public Methods-----------------------------------
         975             :   !
         976             :   ! ================================================================================
         977             :
         978             : !======================================================================
         979             : ! subroutine MeshOpen
         980             : !======================================================================
         981             :
         982           0 :   subroutine MeshOpen(mesh_file_name, par)
         983             :     use parallel_mod, only: parallel_t
         984             :     use cam_logfile,  only: iulog
         985             :
         986             :     character (len=*), intent(in) :: mesh_file_name
         987             :     type (parallel_t), intent(in) :: par
         988             :
         989           0 :     integer, allocatable :: node_multiplicity(:)
         990             :     integer              :: k
         991             :
         992           0 :     p_mesh_file_name    = mesh_file_name
         993           0 :     call open_mesh_file ()
         994             :
         995           0 :     p_number_elements   = get_number_of_elements       ()
         996           0 :     p_number_nodes      = get_number_of_nodes          ()
         997           0 :     p_number_blocks     = get_number_of_element_blocks ()
         998           0 :     p_number_dimensions = get_number_of_dimensions     ()
         999             :
        1000           0 :     if (p_number_dimensions /= 3) then
        1001           0 :        call endrun('The number of dimensions must be 3, otherwise the mesh algorithms will not work')
        1002             :     endif
        1003             :
        1004             :     ! Only spheres are allowed in input files.
        1005           0 :     if (par%masterproc) then
        1006           0 :        if (p_number_blocks == 1) then
        1007           0 :           write(iulog,*) "Since the mesh file has only one block, it is assumed to be a sphere."
        1008             :        endif
        1009             :     end if
        1010             :
        1011           0 :     if (p_number_blocks /= 1) then
        1012           0 :        call endrun('Number of elements blocks not exactly 1 (sphere)')
        1013             :     end if
        1014             :
        1015           0 :     p_number_elements_per_face = get_number_of_elements_per_face()
        1016             :     ! Because all elements are in one face, this value must match  p_number_elements
        1017           0 :     if ( p_number_elements /= p_number_elements_per_face) then
        1018           0 :        call endrun('The value of the total number of elements does not match all the elements found in face 1')
        1019             :     end if
        1020             :
        1021           0 :     allocate( p_connectivity(4,p_number_elements_per_face) )
        1022           0 :     p_connectivity(:,:)=0
        1023             :     ! extract the connectivity from the netcdf file
        1024           0 :     call get_face_connectivity()
        1025             :
        1026           0 :     allocate(node_multiplicity(p_number_nodes))
        1027           0 :     call get_node_multiplicity(node_multiplicity)
        1028             :
        1029             :     ! tricky:  For each node with multiplicity n, there are n(n-1) neighbor links
        1030             :     ! created.  But this counts each edge twice, so:  n(n-1) -n
        1031             :     ! Should be the same as SUM(SIZE(GridVertex(i)%nbrs(j)%n),i=1:p_number_elements,j=1:8)
        1032             :     ! p_number_neighbor_edges = dot_product(mult,mult) - 2*sum(mult)
        1033           0 :     p_number_neighbor_edges = 0
        1034           0 :     do k=1,p_number_nodes
        1035           0 :       p_number_neighbor_edges = p_number_neighbor_edges + node_multiplicity(k)*(node_multiplicity(k)-2)
        1036             :     end do
        1037             :
        1038           0 :     deallocate(node_multiplicity)
        1039             :
        1040             :     ! allocate the space for the coordinates, this is used in many functions
        1041           0 :     allocate(p_node_coordinates(p_number_nodes, p_number_dimensions))
        1042           0 :     call get_node_coordinates()
        1043             :
        1044           0 :     if (p_number_elements_per_face /= p_number_elements) then
        1045           0 :        call endrun('MeshOpen: Total number of elements not equal to the number of elements on face 1!')
        1046             :     end if
        1047             :
        1048           0 :   end subroutine MeshOpen
        1049             :
        1050             : !======================================================================
        1051             : ! subroutine MeshClose
        1052             : !
        1053             : ! This routine acts as a destructor cleaning the memory allocated in MeshOpen
        1054             : ! which acts as a constructor allocated dynamical memory for the nodes coordinates.
        1055             : !======================================================================
        1056             :
        1057           0 :   subroutine MeshClose
        1058             :
        1059             :     ! release memory
        1060           0 :     deallocate(p_node_coordinates)
        1061           0 :     deallocate(p_connectivity)
        1062             :     ! let the file go
        1063           0 :     call close_mesh_file ()
        1064             :
        1065           0 :   end subroutine MeshClose
        1066             :
        1067             :
        1068             : !======================================================================
        1069             : ! subroutine MeshPrint
        1070             : !======================================================================
        1071             :
        1072             :
        1073           0 :   subroutine MeshPrint(par)
        1074             :     use parallel_mod, only: parallel_t
        1075             :     use cam_logfile,  only: iulog
        1076             :
        1077             :     type (parallel_t), intent(in) :: par
        1078           0 :     if (par%masterproc) then
        1079           0 :       write(iulog,*) 'This are the values for file ', trim(p_mesh_file_name)
        1080           0 :       write(iulog,*) 'The value for the number of dimensions (num_dim) is ', p_number_dimensions
        1081           0 :       write(iulog,*) 'The number of elements in the mesh file is ', p_number_elements
        1082           0 :       write(iulog,*) 'The number of nodes in the mesh file is ', p_number_nodes
        1083           0 :       write(iulog,*) 'The number of blocks in the mesh file is ',  p_number_blocks
        1084           0 :       write(iulog,*) 'The number of elements in the face 1 (sphere) is ',  p_number_elements_per_face
        1085             :        if ( p_number_elements == p_number_elements) then
        1086           0 :          write(iulog,*) 'The value of the total number of elements does match all the elements found in face 1 (the only face)'
        1087             :        else
        1088             :          write(iulog,*) 'The value of the total number of elements does not match all the elements found in face 1'
        1089             :          write(iulog,*) 'This message should not be appearing, there is something wrong in the code'
        1090             :        endif
        1091           0 :       write(iulog,*) 'The number of neighbor edges ', p_number_neighbor_edges
        1092             :     end if
        1093             :
        1094           0 :   end subroutine MeshPrint
        1095             :
        1096             : !======================================================================
        1097             : ! subroutine MeshCubeTopology
        1098             : !======================================================================
        1099           0 :    subroutine MeshCubeTopology(GridEdge, GridVertex)
        1100             :     use dimensions_mod,         only : np
        1101             :     use coordinate_systems_mod, only : cartesian3D_t, cube_face_number_from_cart
        1102             :     use gridgraph_mod,          only : GridVertex_t
        1103             :     use gridgraph_mod,          only : GridEdge_t
        1104             :     use cube_mod,               only : CubeSetupEdgeIndex
        1105             :     use gridgraph_mod,          only : initgridedge, num_neighbors
        1106             :     use control_mod,            only : north, south, east, west, neast, seast, swest, nwest
        1107             :
        1108             :     type (GridEdge_t),   intent(inout), target :: GridEdge(:)
        1109             :     type (GridVertex_t), intent(inout), target :: GridVertex(:)
        1110             :
        1111             :     real(kind=r8)             :: coordinates(4,3)
        1112             :     real(kind=r8)             :: centroid(3)
        1113             :     type (cartesian3D_t)             :: face_center
        1114             :
        1115             :     integer                          :: i, j, k, ll, loc
        1116           0 :     integer                          :: element_nodes(p_number_elements, 4)
        1117             :     integer                          :: EdgeWgtP,CornerWgt
        1118             :     integer                          :: normal_to_homme_ordering(8)
        1119             :     integer                          :: node_numbers(4)
        1120           0 :     integer, allocatable             :: index_table(:,:)
        1121             :
        1122           0 :     normal_to_homme_ordering(1) = south
        1123           0 :     normal_to_homme_ordering(2) =  east
        1124           0 :     normal_to_homme_ordering(3) = north
        1125           0 :     normal_to_homme_ordering(4) =  west
        1126           0 :     normal_to_homme_ordering(5) = swest
        1127           0 :     normal_to_homme_ordering(6) = seast
        1128           0 :     normal_to_homme_ordering(7) = neast
        1129           0 :     normal_to_homme_ordering(8) = nwest
        1130             :
        1131           0 :     if (SIZE(GridVertex) /= p_number_elements) then
        1132             :        call endrun('MeshCubeTopology: Element count check failed in exodus_mesh. &
        1133           0 :             &Vertex array length not equal to number of elements.')
        1134             :     end if
        1135           0 :     if (p_number_elements_per_face /= p_number_elements) then
        1136             :        call endrun('MeshCubeTopology: Element count check failed in exodus_mesh. &
        1137           0 :             &Element array length not equal to sum of face.')
        1138             :     end if
        1139             :
        1140           0 :     EdgeWgtP = np
        1141           0 :     CornerWgt = 1
        1142             :
        1143             :
        1144           0 :     call mesh_connectivity (element_nodes)
        1145             :
        1146           0 :     do i=1, p_number_elements
        1147           0 :        GridVertex(i)%number           = i
        1148           0 :        GridVertex(i)%face_number      = 0
        1149           0 :        GridVertex(i)%processor_number = 0
        1150           0 :        GridVertex(i)%SpaceCurve       = 0
        1151             :
        1152           0 :        GridVertex(i)%nbrs(:) = 0
        1153           0 :        GridVertex(i)%nbrs_face(:) = 0
        1154           0 :        GridVertex(i)%nbrs_wgt(:) = 0
        1155           0 :        GridVertex(i)%nbrs_wgt_ghost(:) = 1
        1156             :
        1157             :        !each elements has one side neighbor (first 4)
        1158           0 :        GridVertex(i)%nbrs_ptr(1) = 1
        1159           0 :        GridVertex(i)%nbrs_ptr(2) = 2
        1160           0 :        GridVertex(i)%nbrs_ptr(3) = 3
        1161           0 :        GridVertex(i)%nbrs_ptr(4) = 4
        1162             :        !don't know about corners yet
        1163           0 :        GridVertex(i)%nbrs_ptr(5:num_neighbors+1) = 5
        1164             :
        1165             :     end do
        1166             :
        1167             :     !create index table to find neighbors
        1168           0 :     call create_index_table(index_table, element_nodes)
        1169             :
        1170             :     ! side neighbors
        1171           0 :     call find_side_neighbors(GridVertex, normal_to_homme_ordering, element_nodes, EdgeWgtP, index_table)
        1172             :
        1173             :     ! set vertex faces
        1174           0 :     do i=1, p_number_elements
        1175           0 :        node_numbers = element_nodes(i,:)
        1176           0 :        coordinates = p_node_coordinates(node_numbers,:)
        1177           0 :        centroid = SUM(coordinates, dim=1)/4.0_r8
        1178           0 :        face_center%x = centroid(1)
        1179           0 :        face_center%y = centroid(2)
        1180           0 :        face_center%z = centroid(3)
        1181           0 :        GridVertex(i)%face_number = cube_face_number_from_cart(face_center)
        1182             :     end do
        1183             :
        1184             :     ! set side neighbor faces
        1185           0 :     do i=1, p_number_elements
        1186           0 :        do j=1,4 !look at each side
        1187           0 :           k = normal_to_homme_ordering(j)
        1188           0 :           loc =  GridVertex(i)%nbrs_ptr(k)
        1189           0 :           ll = GridVertex(i)%nbrs(loc)
        1190           0 :           GridVertex(i)%nbrs_face(loc) = GridVertex(ll)%face_number
        1191             :        end do
        1192             :     end do
        1193             :
        1194             :     ! find corner neighbor and faces (weights added also)
        1195           0 :     call find_corner_neighbors  (GridVertex, normal_to_homme_ordering, element_nodes, CornerWgt, index_table)
        1196             :
        1197             :     !done with the index table
        1198           0 :     deallocate(index_table)
        1199             :
        1200             :
        1201           0 :     call initgridedge(GridEdge,GridVertex)
        1202           0 :     do i=1,SIZE(GridEdge)
        1203           0 :        call CubeSetupEdgeIndex(GridEdge(i))
        1204             :     enddo
        1205             :
        1206           0 :     call initialize_space_filling_curve(GridVertex, element_nodes)
        1207           0 :   end subroutine MeshCubeTopology
        1208             :
        1209             : !======================================================================
        1210             : ! subroutine MeshSetCoordinates(elem)
        1211             : !======================================================================
        1212             :
        1213           0 :   subroutine MeshSetCoordinates(elem)
        1214             :     use element_mod, only: element_t
        1215             :
        1216             :     type (element_t), intent(inout) :: elem(:)
        1217             :
        1218           0 :     integer                         :: connectivity(p_number_elements,4)
        1219           0 :     integer                         :: node_multiplicity(p_number_nodes)
        1220             :     integer                         :: face_no, i, k, l
        1221             :     integer                         :: number
        1222             :     integer                         :: node_num(4)
        1223             :     real(kind=r8)                   :: coordinates(4,3)
        1224             :     real(kind=r8)                   :: cube_coor  (4,2)
        1225             :
        1226           0 :     connectivity     =0
        1227           0 :     node_multiplicity=0
        1228             :
        1229           0 :     call mesh_connectivity (connectivity)
        1230             :
        1231           0 :     do k=1,p_number_elements
        1232           0 :        node_num = connectivity(k,:)
        1233           0 :        node_multiplicity(node_num(:)) = node_multiplicity(node_num(:)) + 1
        1234             :     end do
        1235             :
        1236           0 :     do k=1,SIZE(elem)
        1237           0 :        number  = elem(k)%vertex%number
        1238           0 :        face_no = elem(k)%vertex%face_number
        1239           0 :        node_num = connectivity(number,:)
        1240           0 :        coordinates = p_node_coordinates(node_num,:)
        1241             :
        1242           0 :       if (6 == p_number_blocks) then
        1243           0 :          call cube_to_cube_coordinates   (cube_coor, coordinates, face_no)
        1244             :       else
        1245           0 :          call sphere_to_cube_coordinates (cube_coor, coordinates, face_no)
        1246             :       end if
        1247             : !      elem(k)%node_numbers         = node_num
        1248             : !      elem(k)%node_multiplicity(:) = node_multiplicity(node_num(:))
        1249           0 :       elem(k)%corners(:)%x         = cube_coor(:,1)
        1250           0 :       elem(k)%corners(:)%y         = cube_coor(:,2)
        1251             :     end do
        1252           0 :   end subroutine MeshSetCoordinates
        1253             :
        1254             : !======================================================================
        1255             : !function MeshCubeEdgeCount()
        1256             : !======================================================================
        1257           0 :   function MeshCubeEdgeCount()  result(nedge)
        1258             :
        1259             :     integer                     :: nedge
        1260           0 :     if (0 == p_number_blocks)  call endrun('MeshCubeEdgeCount called before MeshOpenMesh')
        1261           0 :     if (MeshUseMeshFile) then
        1262             :       ! should be the same as SUM(SIZE(GridVertex(i)%nbrs(j)%n),i=1:p_number_elements,j=1:nInnerElemEdge)
        1263             :       ! the total number of neighbors.
        1264           0 :       nedge = p_number_neighbor_edges
        1265             :     else
        1266           0 :       call endrun('Error in MeshCubeEdgeCount: Should not call for non-exodus mesh file.')
        1267             :     endif
        1268             :
        1269           0 :   end function MeshCubeEdgeCount
        1270             :
        1271           0 :   function MeshCubeElemCount()  result(nelem)
        1272             :
        1273             :     integer                     :: nelem
        1274           0 :     if (0 == p_number_blocks)  call endrun('MeshCubeElemCount called before MeshOpenMesh')
        1275           0 :     if (MeshUseMeshFile) then
        1276           0 :       nelem = p_number_elements
        1277             :     else
        1278           0 :       call endrun('Error in MeshCubeElemCount: Should not call for non-exodus mesh file.')
        1279             :     end if
        1280           0 :   end function MeshCubeElemCount
        1281             :
        1282           0 :   subroutine test_private_methods
        1283             :     implicit none
        1284           0 :     integer                          :: element_nodes(p_number_elements, 4)
        1285           0 :     call mesh_connectivity (element_nodes)
        1286           0 :   end subroutine test_private_methods
        1287             :
        1288             :
        1289             : end module mesh_mod

Generated by: LCOV version 1.14