diff --git a/conda_package/recipe/build.sh b/conda_package/recipe/build.sh index cf69c237b..2bba17ccf 100644 --- a/conda_package/recipe/build.sh +++ b/conda_package/recipe/build.sh @@ -51,3 +51,14 @@ cmake \ .. cmake --build . cmake --install . + +# build and install WAVE MESH tools +cd ${SRC_DIR}/conda_package/ocean/wave_mesh_tools +mkdir build +cd build +cmake \ + -D CMAKE_INSTALL_PREFIX=${PREFIX} \ + -D CMAKE_BUILD_TYPE=Release \ + .. +cmake --build . +cmake --install . diff --git a/ocean/wave_mesh_tools/CMakeLists.txt b/ocean/wave_mesh_tools/CMakeLists.txt new file mode 100644 index 000000000..f20c3d0cf --- /dev/null +++ b/ocean/wave_mesh_tools/CMakeLists.txt @@ -0,0 +1,24 @@ +cmake_minimum_required (VERSION 3.0.2) +enable_language(Fortran) +project (ocean_wave_mesh_tools) + +list (APPEND CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmake) + +set (NETCDF_F90 "YES") +find_package (NetCDF REQUIRED) + +message (STATUS "NETCDF_INCLUDES=${NETCDF_INCLUDES}") +message (STATUS "NETCDF_LIBRARIES=${NETCDF_LIBRARIES}") + +include_directories(${NETCDF_INCLUDES}) + +add_executable (ocean_rotate_wave_mesh rotate.f90 rotate_mod.f90 read_write_gmsh.f90 grid_file_mod.f90 globals.f90) +target_link_libraries (ocean_rotate_wave_mesh ${NETCDF_LIBRARIES}) + +add_executable (ocean_cull_wave_mesh cull_wave_mesh.f90 in_cell_mod.f90 edge_connectivity_mod.f90 fix_elements.f90 globals.f90 grid_file_mod.f90 kdtree2.f90 read_write_gmsh.f90 write_vtk.f90) +target_link_libraries (ocean_cull_wave_mesh ${NETCDF_LIBRARIES}) + +add_executable (ocean_scrip_wave_mesh scrip.f90 wmscrpmd.f90 read_write_gmsh.f90) +target_link_libraries (ocean_scrip_wave_mesh ${NETCDF_LIBRARIES}) + +install (TARGETS ocean_rotate_wave_mesh ocean_cull_wave_mesh ocean_scrip_wave_mesh DESTINATION bin) diff --git a/ocean/wave_mesh_tools/cmake/FindNetCDF.cmake b/ocean/wave_mesh_tools/cmake/FindNetCDF.cmake new file mode 100644 index 000000000..4db4b0d38 --- /dev/null +++ b/ocean/wave_mesh_tools/cmake/FindNetCDF.cmake @@ -0,0 +1,71 @@ +# - Find NetCDF +# Find the native NetCDF includes and library +# +# NETCDF_INCLUDES - where to find netcdf.h, etc +# NETCDF_LIBRARIES - Link these libraries when using NetCDF +# NETCDF_FOUND - True if NetCDF found including required interfaces (see below) +# +# Your package can require certain interfaces to be FOUND by setting these +# +# NETCDF_CXX - require the C++ interface and link the C++ library +# NETCDF_F77 - require the F77 interface and link the fortran library +# NETCDF_F90 - require the F90 interface and link the fortran library +# +# The following are not for general use and are included in +# NETCDF_LIBRARIES if the corresponding option above is set. +# +# NETCDF_LIBRARIES_C - Just the C interface +# NETCDF_LIBRARIES_CXX - C++ interface, if available +# NETCDF_LIBRARIES_F77 - Fortran 77 interface, if available +# NETCDF_LIBRARIES_F90 - Fortran 90 interface, if available +# +# Normal usage would be: +# set (NETCDF_F90 "YES") +# find_package (NetCDF REQUIRED) +# target_link_libraries (uses_f90_interface ${NETCDF_LIBRARIES}) +# target_link_libraries (only_uses_c_interface ${NETCDF_LIBRARIES_C}) + +if (NETCDF_INCLUDES AND NETCDF_LIBRARIES) + # Already in cache, be silent + set (NETCDF_FIND_QUIETLY TRUE) +endif (NETCDF_INCLUDES AND NETCDF_LIBRARIES) + +find_path (NETCDF_INCLUDES netcdf.h + HINTS NETCDF_DIR ENV NETCDF_DIR) + +find_library (NETCDF_LIBRARIES_C NAMES netcdf) +mark_as_advanced(NETCDF_LIBRARIES_C) + +set (NetCDF_has_interfaces "YES") # will be set to NO if we're missing any interfaces +set (NetCDF_libs "${NETCDF_LIBRARIES_C}") + +get_filename_component (NetCDF_lib_dirs "${NETCDF_LIBRARIES_C}" PATH) + +macro (NetCDF_check_interface lang header libs) + if (NETCDF_${lang}) + find_path (NETCDF_INCLUDES_${lang} NAMES ${header} + HINTS "${NETCDF_INCLUDES}" NO_DEFAULT_PATH) + find_library (NETCDF_LIBRARIES_${lang} NAMES ${libs} + HINTS "${NetCDF_lib_dirs}" NO_DEFAULT_PATH) + mark_as_advanced (NETCDF_INCLUDES_${lang} NETCDF_LIBRARIES_${lang}) + if (NETCDF_INCLUDES_${lang} AND NETCDF_LIBRARIES_${lang}) + list (INSERT NetCDF_libs 0 ${NETCDF_LIBRARIES_${lang}}) # prepend so that -lnetcdf is last + else (NETCDF_INCLUDES_${lang} AND NETCDF_LIBRARIES_${lang}) + set (NetCDF_has_interfaces "NO") + message (STATUS "Failed to find NetCDF interface for ${lang}") + endif (NETCDF_INCLUDES_${lang} AND NETCDF_LIBRARIES_${lang}) + endif (NETCDF_${lang}) +endmacro (NetCDF_check_interface) + +NetCDF_check_interface (CXX netcdfcpp.h netcdf_c++) +NetCDF_check_interface (F77 netcdf.inc netcdff) +NetCDF_check_interface (F90 netcdf.mod netcdff) + +set (NETCDF_LIBRARIES "${NetCDF_libs}" CACHE STRING "All NetCDF libraries required for interface level") + +# handle the QUIETLY and REQUIRED arguments and set NETCDF_FOUND to TRUE if +# all listed variables are TRUE +include (FindPackageHandleStandardArgs) +find_package_handle_standard_args (NetCDF DEFAULT_MSG NETCDF_LIBRARIES NETCDF_INCLUDES NetCDF_has_interfaces) + +mark_as_advanced (NETCDF_LIBRARIES NETCDF_INCLUDES) diff --git a/ocean/wave_mesh_tools/cmake/LICENSE b/ocean/wave_mesh_tools/cmake/LICENSE new file mode 100644 index 000000000..49bac9892 --- /dev/null +++ b/ocean/wave_mesh_tools/cmake/LICENSE @@ -0,0 +1,23 @@ +Copyright $(git shortlog -s) +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +* Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +* Redistributions in binary form must reproduce the above copyright notice, this + list of conditions and the following disclaimer in the documentation and/or + other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/ocean/wave_mesh_tools/cmake/README b/ocean/wave_mesh_tools/cmake/README new file mode 100644 index 000000000..b4811763c --- /dev/null +++ b/ocean/wave_mesh_tools/cmake/README @@ -0,0 +1 @@ +These files are from https://github.com/jedbrown/cmake-modules diff --git a/ocean/wave_mesh_tools/cull_wave_mesh.f90 b/ocean/wave_mesh_tools/cull_wave_mesh.f90 new file mode 100644 index 000000000..a8712fe8d --- /dev/null +++ b/ocean/wave_mesh_tools/cull_wave_mesh.f90 @@ -0,0 +1,295 @@ + PROGRAM cull_waves_mesh + + USE netcdf + USE in_cell_mod, ONLY: in_cell_init, pt_in_cell, check, pi, lonlat2xyz, R + USE globals, ONLY: rp, nverts + USE write_vtk, ONLY: write_vtk_file + USE read_write_gmsh, ONLY: write_header, write_nodes, write_elements + USE edge_connectivity_mod, ONLY: elements_per_node, find_edge_pairs, find_interior_edges, & + find_element_edges, find_neighbor_elements, find_boundary_nodes + USE fix_elements, ONLY: fix_single_node_connections_across_islands, & + fix_element_node_orientation, & + create_new_ect, & + remove_unconnected_nodes, & + add_elements_to_ect + + IMPLICIT NONE + + INTEGER :: waves_ncid + INTEGER :: ocean_ncid + INTEGER :: nCells_dimid, nVertices_dimid + INTEGER :: cellsOnVertex_varid, lonCell_varid, latCell_varid + INTEGER :: bottomDepth_varid + INTEGER :: iceMask_varid + + CHARACTER(100) :: waves_mesh_file + CHARACTER(100) :: waves_mesh_culled_vtk + CHARACTER(100) :: waves_mesh_culled_gmsh + CHARACTER(100) :: ocean_mesh_file + + INTEGER :: ne_waves, nn_waves + INTEGER :: ne_new, nn_new + INTEGER, DIMENSION(:,:), ALLOCATABLE :: ect_waves + INTEGER, DIMENSION(:,:), ALLOCATABLE :: ect_new + REAL(rp), DIMENSION(:), ALLOCATABLE :: x_waves, y_waves + REAL(rp), DIMENSION(:,:), ALLOCATABLE :: xy_waves + REAL(rp), DIMENSION(:), ALLOCATABLE :: x_new, y_new + REAL(rp), DIMENSION(:), ALLOCATABLE :: depth_waves + REAL(rp), DIMENSION(:), ALLOCATABLE :: depth_new + REAL(rp) :: xy(2) + REAL(rp), DIMENSION(:,:), ALLOCATABLE :: xyz + REAL(rp), DIMENSION(:,:), ALLOCATABLE :: xy_new + REAL(rp) :: x1, x2, x3 + REAL(rp) :: y1, y2, y3 + REAL(rp) :: x2x1, x1x3, area + + INTEGER :: ncells_ocean + REAL(rp), DIMENSION(:), ALLOCATABLE :: depth_ocean + REAL(rp), DIMENSION(:), ALLOCATABLE :: icemask_ocean + + INTEGER :: el, nd, ged + INTEGER :: i,j + INTEGER :: in_cell + INTEGER :: tmp + INTEGER :: mnepn + INTEGER :: ned + INTEGER, DIMENSION(:), ALLOCATABLE :: nepn + INTEGER, DIMENSION(:,:), ALLOCATABLE :: epn + INTEGER, DIMENSION(:,:), ALLOCATABLE :: ged2el + INTEGER, DIMENSION(:,:), ALLOCATABLE :: ged2nn + INTEGER, DIMENSION(:,:), ALLOCATABLE :: ged2led + INTEGER, DIMENSION(:), ALLOCATABLE :: el_type + INTEGER :: nied + INTEGER, DIMENSION(:), ALLOCATABLE :: iedn + INTEGER, DIMENSION(:), ALLOCATABLE :: ed_type + INTEGER, DIMENSION(:), ALLOCATABLE :: recv_edge + INTEGER, DIMENSION(:,:), ALLOCATABLE :: el2ged + INTEGER, DIMENSION(:,:), ALLOCATABLE :: el2el + INTEGER :: nbed + INTEGER, DIMENSION(:), ALLOCATABLE :: bedn + INTEGER :: nbnd + INTEGER, DIMENSION(:), ALLOCATABLE :: bndn + INTEGER, DIMENSION(:), ALLOCATABLE :: new_node_numbers + INTEGER, DIMENSION(:), ALLOCATABLE :: nd_flag + INTEGER, DIMENSION(:), ALLOCATABLE :: elem_flag + INTEGER :: nfill_elements + INTEGER, DIMENSION(:,:), ALLOCATABLE :: fill_elements + INTEGER, DIMENSION(:), ALLOCATABLE :: ed_flag + INTEGER :: nd_count + INTEGER :: nd_start, nd_next + INTEGER, DIMENSION(:), ALLOCATABLE :: nd_found + + INTEGER, DIMENSION(:), ALLOCATABLE :: wave_node_in_ocean + INTEGER, DIMENSION(:), ALLOCATABLE :: wave_elements_keep + + NAMELIST / inputs / waves_mesh_file, ocean_mesh_file + NAMELIST / output / waves_mesh_culled_vtk, waves_mesh_culled_gmsh + + OPEN(UNIT=15, FILE='cull_waves_mesh.nml') + READ(UNIT=15, NML=inputs) + READ(UNIT=15, NML=output) + CLOSE(15) + + PRINT*, 'waves_mesh_file = ', TRIM(ADJUSTL(waves_mesh_file)) + PRINT*, 'ocean_mesh_file = ', TRIM(ADJUSTL(ocean_mesh_file)) + PRINT*, 'waves_mesh_file_vtk = ', TRIM(ADJUSTL(waves_mesh_culled_vtk)) + PRINT*, 'waves_mesh_file_gmsh = ', TRIM(ADJUSTL(waves_mesh_culled_gmsh)) + CALL check(NF90_OPEN(TRIM(ADJUSTL(waves_mesh_file)), NF90_NOWRITE, waves_ncid)) + CALL check(NF90_INQ_DIMID(waves_ncid, 'nCells', nCells_dimid)) + CALL check(NF90_INQ_DIMID(waves_ncid, 'nVertices', nVertices_dimid)) + + CALL check(NF90_INQUIRE_DIMENSION(waves_ncid, nCells_dimid, len=nn_waves)) + CALL check(NF90_INQUIRE_DIMENSION(waves_ncid, nVertices_dimid, len=ne_waves)) + + CALL check(NF90_INQ_VARID(waves_ncid, 'cellsOnVertex', cellsOnVertex_varid)) + CALL check(NF90_INQ_VARID(waves_ncid, 'lonCell', lonCell_varid)) + CALL check(NF90_INQ_VARID(waves_ncid, 'latCell', latCell_varid)) + + ALLOCATE(x_waves(nn_waves), y_waves(nn_waves)) + CALL check(NF90_GET_VAR(waves_ncid, lonCell_varid, x_waves)) + CALL check(NF90_GET_VAR(waves_ncid, latCell_varid, y_waves)) + ALLOCATE(ect_waves(3, ne_waves)) + CALL check(NF90_GET_VAR(waves_ncid, cellsOnVertex_varid, ect_waves)) + CALL check(NF90_CLOSE(waves_ncid)) + + CALL check(NF90_OPEN(ocean_mesh_file, NF90_NOWRITE, ocean_ncid)) + CALL check(NF90_INQ_DIMID(ocean_ncid, 'nCells', nCells_dimid)) + CALL check(NF90_INQUIRE_DIMENSION(ocean_ncid, nCells_dimid, len=ncells_ocean)) + CALL check(NF90_INQ_VARID(ocean_ncid, 'bottomDepth', bottomDepth_varid)) + ALLOCATE(depth_ocean(ncells_ocean)) + CALL check(NF90_GET_VAR(ocean_ncid, bottomDepth_varid, depth_ocean)) + CALL check(NF90_INQ_VARID(ocean_ncid, 'landIceMask', iceMask_varid)) + ALLOCATE(icemask_ocean(ncells_ocean)) + CALL check(NF90_GET_VAR(ocean_ncid, iceMask_varid, icemask_ocean)) + CALL check(NF90_CLOSE(ocean_ncid)) + + + + CALL in_cell_init(TRIM(ADJUSTL(ocean_mesh_file))) + + ALLOCATE(xy_waves(2,nn_waves)) + DO i = 1,nn_waves + xy_waves(1,i) = x_waves(i) + xy_waves(2,i) = y_waves(i) + ENDDO + + + ALLOCATE(wave_node_in_ocean(nn_waves)) + wave_node_in_ocean = 0 + ALLOCATE(wave_elements_keep(ne_waves)) + wave_elements_keep = 0 + ALLOCATE(depth_waves(nn_waves)) + depth_waves = -2d0 + + ! Determine which elements to keep based on nodes found in ocean mesh + ne_new = 0 +elems:DO el = 1,ne_waves + + IF (mod(el,1000) == 0) THEN + PRINT*, el + ENDIF + + nds: DO i = 1,3 + + nd = ect_waves(i,el) + + IF (nd == 0) THEN + wave_elements_keep(el) = 0 + EXIT nds + ENDIF + + IF (wave_node_in_ocean(nd) == 1) THEN + wave_elements_keep(el) = 1 + ne_new = ne_new + 1 + CYCLE nds + ENDIF + + xy(1) = x_waves(nd) + xy(2) = y_waves(nd) + + CALL pt_in_cell(xy,in_cell) + + IF (icemask_ocean(in_cell) == 0) THEN + wave_node_in_ocean(nd) = 1 + depth_waves(nd) = depth_ocean(in_cell) + IF (wave_elements_keep(el) == 0) THEN + wave_elements_keep(el) = 1 + ne_new = ne_new + 1 + ENDIF + ENDIF + + ENDDO nds + ENDDO elems + + + ! Create new element connectivity table + CALL create_new_ect(ne_waves,ect_waves,wave_elements_keep,ne_new,ect_new) + PRINT*, "" + PRINT*, "Elements in original waves mesh: ", ne_waves + PRINT*, "Elements in culled waves mesh:", ne_new + + ! Find number of elements per node + ALLOCATE(el_type(ne_new)) + el_type = 1 + CALL elements_per_node(ne_new,nn_waves,nverts,el_type,ect_new,nepn,mnepn,epn) + + ! Build new coordinate list + DEALLOCATE(ect_waves) + ALLOCATE(ect_waves(3,ne_new)) + ect_waves = ect_new + CALL remove_unconnected_nodes(nn_waves,nepn,xy_waves,ne_new,ect_waves,nn_new,xy_new,ect_new,depth_waves,depth_new) + PRINT*, "Nodes in original waves mesh: ", nn_waves + PRINT*, "Nodes in culled waves mesh:", nn_new + + + ! Convert to Cartesian coordinates for vtk output + ALLOCATE(xyz(3,nn_new)) + DO i = 1,nn_new + CALL lonlat2xyz(xy_new(1,i),xy_new(2,i),R,xyz(1,i),xyz(2,i),xyz(3,i)) + ENDDO + + ! Convert lon,lat to degrees and adjust lon to -180,180 range + DO i = 1,nn_new + xy_new(1,i) = xy_new(1,i)*180d0/pi + xy_new(2,i) = xy_new(2,i)*180d0/pi + IF (xy_new(1,i) > 180d0) THEN + xy_new(1,i) = xy_new(1,i) - 360d0 + ENDIF + ENDDO + + ! Compute edge connectivity and indentify single node connections + DEALLOCATE(nepn,epn) + CALL elements_per_node(ne_new,nn_new,nverts,el_type,ect_new,nepn,mnepn,epn) + CALL find_edge_pairs(ne_new,nverts,el_type,ect_new,nepn,epn,ned,ged2el,ged2nn,ged2led) + CALL find_interior_edges(ned,ged2el,nied,iedn,ed_type,recv_edge,nbed,bedn) + CALL find_element_edges(ne_new,ned,ged2el,ged2led,el2ged) + CALL find_neighbor_elements(ne_new,ned,ged2el,ged2led,el2el) + CALL find_boundary_nodes(nn_new,nbed,bedn,ged2nn,nbnd,bndn) + ALLOCATE(elem_flag(ne_new)) + CALL fix_single_node_connections_across_islands(nn_new,nbnd,bndn,nbed,bedn,nepn,epn,ged2nn,nd_flag,elem_flag) + + ! Find nodes in order to add elements to correct single node connections + PRINT*, "Fixing single node connections" + ALLOCATE(fill_elements(3,ne_new)) + ALLOCATE(nd_found(nbed)) + ALLOCATE(ed_flag(ned)) + nfill_elements = 0 + DO nd = 1,nn_new + IF (nd_flag(nd) > 2) THEN + nd_start = nd + nd_next = nd + nd_count = 1 + nd_found(nd_count) = nd_next + ed_flag = 0 + link: DO + search: DO i = 1,nbed + + ged = bedn(i) + IF (ed_flag(ged) == 1) THEN + CYCLE search + ENDIF + + IF (ged2nn(1,ged) == nd_next) THEN + nd_next = ged2nn(2,ged) + ed_flag(ged) = 1 + nd_count = nd_count + 1 + nd_found(nd_count) = nd_next + ELSE IF (ged2nn(2,ged) == nd_next) THEN + nd_next = ged2nn(1,ged) + ed_flag(ged) = 1 + nd_count = nd_count + 1 + nd_found(nd_count) = nd_next + ENDIF + + IF (ed_flag(ged) == 1) THEN + IF (nd_next == nd_start) THEN + EXIT link + ENDIF + EXIT search + ENDIF + + ENDDO search + ENDDO link + + nfill_elements = nfill_elements + 1 + fill_elements(1,nfill_elements) = nd_found(1) + fill_elements(2,nfill_elements) = nd_found(2) + fill_elements(3,nfill_elements) = nd_found(nd_count-1) + PRINT*, "Adding element: ", (fill_elements(j,nfill_elements), j = 1,3) + + ENDIF + ENDDO + + CALL add_elements_to_ect(ne_new,ect_new,nfill_elements,fill_elements) + CALL fix_element_node_orientation(ne_new,xy_new,ect_new) + + + + ! Write mesh + CALL write_vtk_file(TRIM(ADJUSTL(waves_mesh_culled_vtk)),nn_new,xyz,ne_new,ect_new,depth_new) + + CALL write_header(TRIM(ADJUSTL(waves_mesh_culled_gmsh))) + CALL write_nodes(nn_new,xy_new,depth_new) + CALL write_elements(ne_new,ect_new) + + END PROGRAM cull_waves_mesh diff --git a/ocean/wave_mesh_tools/edge_connectivity_mod.f90 b/ocean/wave_mesh_tools/edge_connectivity_mod.f90 new file mode 100644 index 000000000..228f9223e --- /dev/null +++ b/ocean/wave_mesh_tools/edge_connectivity_mod.f90 @@ -0,0 +1,804 @@ + MODULE edge_connectivity_mod + + USE globals, ONLY: rp + + IMPLICIT NONE + + INTEGER, SAVE :: nbed + INTEGER, SAVE, DIMENSION(:), ALLOCATABLE :: bedn ! allocated in find_interior_edges, deallocated in find_flow_edges + + CONTAINS + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE elements_per_node(ne,nn,nverts,el_type,vct,nepn,mnepn,epn) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ne + INTEGER, INTENT(IN) :: nn + INTEGER, DIMENSION(:), INTENT(IN) :: nverts + INTEGER, DIMENSION(:), INTENT(IN) :: el_type + INTEGER, DIMENSION(:,:), INTENT(IN) :: vct + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: nepn + INTEGER, INTENT(OUT) :: mnepn + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: epn + + INTEGER :: el + INTEGER :: nd + INTEGER :: alloc_status + INTEGER :: nvert + INTEGER :: n + + + + ! count the number of elements per node + + ALLOCATE(nepn(nn), STAT=alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + nepn(:) = 0 + DO el = 1,ne + nvert = nverts(el_type(el)) + DO nd = 1,nvert + n = vct(nd,el) + nepn(n) = nepn(n) + 1 + ENDDO + ENDDO + + + ! find the maximum number of elements per node + + mnepn = 0 + DO nd = 1,nn + IF(nepn(nd) > mnepn) THEN + mnepn = nepn(nd) + ENDIF + ENDDO + + + + ! find the elements associated with each node + + ALLOCATE(epn(mnepn,nn), STAT=alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + nepn(:) = 0 + DO el = 1,ne + nvert = nverts(el_type(el)) + DO nd = 1,nvert + n = vct(nd,el) + nepn(n) = nepn(n) + 1 + epn(nepn(n),n) = el + ENDDO + ENDDO + + RETURN + END SUBROUTINE elements_per_node + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE find_edge_pairs(ne,nverts,el_type,vct,nepn,epn,ned,ged2el,ged2nn,ged2led) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ne + INTEGER, DIMENSION(:), INTENT(IN) :: nverts + INTEGER, DIMENSION(:), INTENT(IN) :: el_type + INTEGER, DIMENSION(:,:), INTENT(IN) :: vct + INTEGER, DIMENSION(:), INTENT(IN) :: nepn + INTEGER, DIMENSION(:,:), INTENT(IN) :: epn + INTEGER, INTENT(OUT) :: ned + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: ged2el + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: ged2nn + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: ged2led + + INTEGER :: el + INTEGER :: nd + INTEGER :: alloc_status + INTEGER :: el1,el2 + INTEGER :: led1,led2 + INTEGER :: nvert1,nvert2 + INTEGER :: n1ed1,n2ed1 + INTEGER :: n1ed2,n2ed2 + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ged2nn_temp, ged2el_temp, ged2led_temp + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: edflag + + + ALLOCATE(ged2nn_temp(2,4*ne),ged2el_temp(2,4*ne),ged2led_temp(2,4*ne), STAT = alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + ALLOCATE(edflag(4,ne), STAT = alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + + ned = 0 + edflag = 0 + ged2nn_temp = 0 + ged2led_temp = 0 + ged2el_temp = 0 + + elem1: DO el1 = 1,ne ! loop through trial elements + + nvert1 = nverts(el_type(el1)) + + local_ed1: DO led1 = 1,nvert1 ! loop through trial edges + + IF(edflag(led1,el1) == 1) THEN ! skip if edge has already been flagged + CYCLE local_ed1 + ENDIF + + ned = ned + 1 ! increment edge number + + n1ed1 = vct(mod(led1+0,nvert1)+1,el1) ! find nodes on trial edge + n2ed1 = vct(mod(led1+1,nvert1)+1,el1) + + ged2nn_temp(1,ned) = n1ed1 ! set nodes on global edge # ned + ged2nn_temp(2,ned) = n2ed1 + + ged2led_temp(1,ned) = led1 ! set local edge number of first element sharing the edge + ged2el_temp(1,ned) = el1 ! set the first element that shares the edge + + edflag(led1,el1) = 1 ! flag the edge so it is not repeated + + elem2: DO el = 1,nepn(n1ed1) ! loop through test elements (only those that contain node n1ed1) + + el2 = epn(el,n1ed1) ! choose a test element that contains node n1ed1 + + IF(el2 == el1) THEN ! skip if the test element is the same as the trial element + CYCLE elem2 + ENDIF + + nvert2 = nverts(el_type(el2)) + + local_ed2: DO led2 = 1,nvert2 ! loop through local test edge numbers + + n1ed2 = vct(MOD(led2+0,nvert2)+1,el2) ! find nodes on test edge + n2ed2 = vct(MOD(led2+1,nvert2)+1,el2) + + IF(((n1ed1 == n1ed2) .AND. (n2ed1 == n2ed2)) .OR. & ! check if nodes on trial edge matches test edge + ((n1ed1 == n2ed2) .AND. (n2ed1 == n1ed2))) THEN + + ged2led_temp(2,ned) = led2 ! set local edge number of second element sharing the edge + ged2el_temp(2,ned) = el2 ! set the second element that shares the edge + + edflag(led2,el2) = 1 ! flag the edge so it is not repeated + + EXIT elem2 + ENDIF + + ENDDO local_ed2 + + ENDDO elem2 + + ENDDO local_ed1 + + ENDDO elem1 + + ALLOCATE(ged2nn(2,ned),ged2el(2,ned),ged2led(2,ned), STAT=alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + ged2nn(:,1:ned) = ged2nn_temp(:,1:ned) + ged2el(:,1:ned) = ged2el_temp(:,1:ned) + ged2led(:,1:ned) = ged2led_temp(:,1:ned) + + RETURN + END SUBROUTINE find_edge_pairs + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE find_interior_edges(ned,ged2el,nied,iedn,ed_type,recv_edge,nbou_edges,bou_edges) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ned + INTEGER, DIMENSION(:,:), INTENT(IN) :: ged2el + INTEGER, INTENT(OUT) :: nied + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: iedn + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: ed_type + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: recv_edge + INTEGER, INTENT(OUT), OPTIONAL :: nbou_edges + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT), OPTIONAL :: bou_edges + + INTEGER :: el + INTEGER :: nd + INTEGER :: alloc_status + INTEGER :: ged,ed + INTEGER :: el1,el2 + INTEGER, DIMENSION(:), ALLOCATABLE :: iedn_temp + + ALLOCATE(iedn_temp(ned), STAT = alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + ALLOCATE(recv_edge(ned),ed_type(ned), STAT = alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + IF (.NOT. ALLOCATED(bedn)) THEN + ALLOCATE(bedn(ned), STAT = alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + ENDIF + + recv_edge = 1 + ed_type = -999 + + nied = 0 + nbed = 0 + iedn_temp(:) = 0 + DO ged = 1,ned + el1 = ged2el(1,ged) + el2 = ged2el(2,ged) + IF ((el1 /= 0) .AND. (el2 /= 0)) THEN + nied = nied + 1 + iedn_temp(nied) = ged + recv_edge(ged) = 0 + ed_type(ged) = 0 + ELSE + nbed = nbed + 1 + bedn(nbed) = ged + ENDIF + ENDDO + + ALLOCATE(iedn(nied), STAT=alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + iedn(1:nied) = iedn_temp(1:nied) + + ! ed_type : interior edges = 0 + ! open edges = 1 + ! no normal flow edges = 10 + ! specified flow edges = 12 + ! recieve edges = -1 + + IF (PRESENT(nbou_edges).AND.PRESENT(bou_edges)) THEN + ALLOCATE(bou_edges(nbed)) + DO ed = 1,nbed + bou_edges(ed) = bedn(ed) + ENDDO + nbou_edges = nbed + ENDIF + + RETURN + END SUBROUTINE find_interior_edges + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE find_open_edges(nope,obseg,obnds,ged2nn,nobed,obedn,ed_type,recv_edge) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nope + INTEGER, DIMENSION(:), INTENT(IN) :: obseg + INTEGER, DIMENSION(:,:), INTENT(IN) :: obnds + INTEGER, DIMENSION(:,:), INTENT(IN) :: ged2nn + INTEGER, INTENT(OUT) :: nobed + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: obedn + INTEGER, DIMENSION(:), INTENT(INOUT) :: ed_type + INTEGER, DIMENSION(:), INTENT(INOUT) :: recv_edge + + + INTEGER :: el + INTEGER :: alloc_status + INTEGER :: seg,nd,ed,ged + INTEGER :: n1bed,n2bed + INTEGER :: n1ed2,n2ed2 + INTEGER, DIMENSION(:), ALLOCATABLE :: obedn_temp + + ALLOCATE(obedn_temp(nbed), STAT = alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + nobed = 0 + obedn_temp = 0 + + DO seg = 1,nope + DO nd = 1,obseg(seg)-1 + + n1bed = obnds(nd,seg) + n2bed = obnds(nd+1,seg) + + edges1: DO ed = 1,nbed + + ged = bedn(ed) + n1ed2 = ged2nn(1,ged) + n2ed2 = ged2nn(2,ged) + + IF(((n1ed2 == n1bed).AND.(n2ed2 == n2bed)).OR. & + ((n1ed2 == n2bed).AND.(n2ed2 == n1bed))) THEN + + nobed = nobed + 1 + obedn_temp(nobed) = ged + recv_edge(ged) = 0 + ed_type(ged) = 1 + + EXIT edges1 + ENDIF + + ENDDO edges1 + ENDDO + ENDDO + + ALLOCATE(obedn(nobed), STAT=alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + obedn(1:nobed) = obedn_temp(1:nobed) + + RETURN + END SUBROUTINE find_open_edges + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE find_flow_edges(nbou,fbseg,fbnds,ged2nn,nnfbed,nfbedn,nfbednn,nfbed,fbedn,recv_edge,ed_type) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nbou + INTEGER, DIMENSION(:,:), INTENT(IN) :: fbseg + INTEGER, DIMENSION(:,:), INTENT(IN) :: fbnds + INTEGER, DIMENSION(:,:), INTENT(IN) :: ged2nn + INTEGER, INTENT(OUT) :: nnfbed + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: nfbedn + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: nfbednn + INTEGER, INTENT(OUT) :: nfbed + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: fbedn + INTEGER, DIMENSION(:), INTENT(INOUT) :: recv_edge + INTEGER, DIMENSION(:), INTENT(INOUT) :: ed_type + + INTEGER :: el + INTEGER :: nd + INTEGER :: alloc_status + INTEGER :: seg,ed + INTEGER :: ged,segtype + INTEGER :: n1bed,n2bed + INTEGER :: n1ed2,n2ed2 + INTEGER :: found + INTEGER, DIMENSION(:), ALLOCATABLE :: nfbedn_temp + INTEGER, DIMENSION(:,:), ALLOCATABLE :: nfbednn_temp + INTEGER, DIMENSION(:), ALLOCATABLE :: fbedn_temp + + + ALLOCATE(nfbedn_temp(nbed),nfbednn_temp(nbed,3),fbedn_temp(nbed), STAT=alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + + nnfbed = 0 ! # no normal flow edge + nfbedn_temp = 0 + nfbednn_temp = 0 + + nfbed = 0 ! # flow specifed edges + fbedn_temp = 0 + + found = 0 + + DO seg = 1,nbou + + segtype = fbseg(2,seg) + + DO nd = 1,fbseg(1,seg)-1 + + n1bed = fbnds(nd,seg) + n2bed = fbnds(nd+1,seg) + found = 0 + + edges2: DO ed = 1,nbed + + ged = bedn(ed) + n1ed2 = ged2nn(1,ged) + n2ed2 = ged2nn(2,ged) + + IF(((n1ed2 == n1bed).AND.(n2ed2 == n2bed)).OR. & + ((n1ed2 == n2bed).AND.(n2ed2 == n1bed))) THEN + + ! no normal flow edges + IF( segtype == 0 .OR. segtype == 10 .OR. segtype == 20 .OR. & ! land boundaries + segtype == 1 .OR. segtype == 11 .OR. segtype == 21 ) THEN ! island boundaries + + nnfbed = nnfbed + 1 + nfbedn_temp(nnfbed) = ged + nfbednn_temp(nnfbed,1) = seg + nfbednn_temp(nnfbed,2) = n1bed + nfbednn_temp(nnfbed,3) = nd + recv_edge(ged) = 0 + ed_type(ged) = 10 + found = 1 + + EXIT edges2 + ENDIF + + ! specified normal flow edges + IF ( segtype == 2 .OR. segtype == 12 .OR. segtype == 22 ) THEN + + nfbed = nfbed + 1 + fbedn_temp(nfbed) = ged + recv_edge(ged) = 0 + ed_type(ged) = 12 + found = 1 + + EXIT edges2 + ENDIF + + ENDIF + ENDDO edges2 + + IF (found == 0) THEN +! PRINT*, seg, nd, fbseg(1,seg), segtype,myrank + PRINT "(4(A,I9))", "seg = ",seg, " nd = ",nd, " #nds in seg = ",fbseg(1,seg), " type = ",segtype + PRINT "(A)", " edge not found" + ENDIF + + ENDDO + ENDDO + + ALLOCATE(nfbedn(nnfbed),nfbednn(nnfbed,3),fbedn(nfbed), STAT=alloc_status) + DEALLOCATE(bedn, STAT=alloc_status) + + + nfbedn(1:nnfbed) = nfbedn_temp(1:nnfbed) + nfbednn(1:nnfbed,1:3) = nfbednn_temp(1:nnfbed,1:3) + fbedn(1:nfbed) = fbedn_temp(1:nfbed) + + END SUBROUTINE find_flow_edges + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE find_recieve_edges(ned,recv_edge,nred,redn,ed_type) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ned + INTEGER, DIMENSION(:), INTENT(IN) :: recv_edge + INTEGER, INTENT(OUT) :: nred + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: redn + INTEGER, DIMENSION(:), INTENT(INOUT) :: ed_type + + INTEGER :: el + INTEGER :: nd + INTEGER :: alloc_status + INTEGER :: ged + + nred = 0 + +#ifdef CMPI + DO ged = 1,ned + IF(recv_edge(ged) == 1) THEN + nred = nred + 1 + ENDIF + ENDDO + + ALLOCATE(redn(nred), STAT = alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + nred = 0 + DO ged = 1,ned + IF(recv_edge(ged) == 1) THEN + nred = nred + 1 + redn(nred) = ged + ed_type(ged) = -1 + ENDIF + ENDDO +#endif + + RETURN + END SUBROUTINE find_recieve_edges + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE find_element_edges(ne,ned,ged2el,ged2led,el2ged) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ne + INTEGER, INTENT(IN) :: ned + INTEGER, DIMENSION(:,:), INTENT(IN) :: ged2el + INTEGER, DIMENSION(:,:), INTENT(IN) :: ged2led + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: el2ged + + INTEGER :: ged + INTEGER :: el1,el2 + INTEGER :: led + INTEGER :: alloc_status + + ALLOCATE(el2ged(ne,4), STAT = alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + el2ged = 0 + + DO ged = 1,ned + + el1 = ged2el(1,ged) + el2 = ged2el(2,ged) + + IF (el1 /= 0) THEN + led = ged2led(1,ged) + el2ged(el1,led) = ged + ENDIF + + IF (el2 /= 0) THEN + led = ged2led(2,ged) + el2ged(el2,led) = ged + ENDIF + + ENDDO + + RETURN + END SUBROUTINE find_element_edges + + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE find_neighbor_elements(ne,ned,ged2el,ged2led,el2el) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ne + INTEGER, INTENT(IN) :: ned + INTEGER, DIMENSION(:,:), INTENT(IN) :: ged2el + INTEGER, DIMENSION(:,:), INTENT(IN) :: ged2led + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: el2el + + INTEGER :: ged + INTEGER :: el1,el2 + INTEGER :: led1,led2 + INTEGER :: alloc_status + + ALLOCATE(el2el(ne,4), STAT=alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + el2el = 0 + + DO ged = 1,ned + + el1 = ged2el(1,ged) + el2 = ged2el(2,ged) + + IF (el1 /= 0 .and. el2 /= 0 ) THEN + led1 = ged2led(1,ged) + led2 = ged2led(2,ged) + + el2el(el1,led1) = el2 + el2el(el2,led2) = el1 + ENDIF + + ENDDO + + + END SUBROUTINE + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE find_adjacent_nodes(nn,mnepn,ned,ged2nn,nadjnds,adjnds) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nn + INTEGER, INTENT(IN) :: mnepn + INTEGER, INTENT(IN) :: ned + INTEGER, DIMENSION(:,:), INTENT(IN) :: ged2nn + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: nadjnds + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: adjnds + + INTEGER :: ed + INTEGER :: n1,n2 + + ALLOCATE(nadjnds(nn)) + ALLOCATE(adjnds(2*mnepn,nn)) + nadjnds = 0 + + DO ed = 1,ned + n1 = ged2nn(1,ed) ! find the node numbers on each edge + n2 = ged2nn(2,ed) + + nadjnds(n1) = nadjnds(n1) + 1 ! count the nodes adjacent to node n1 + nadjnds(n2) = nadjnds(n2) + 1 ! count the nodes adjacent to node n2 + + adjnds(nadjnds(n1),n1) = n2 ! node n2 is adjacent to node n1 + adjnds(nadjnds(n2),n2) = n1 ! node n1 is adjacent to node n2 + ENDDO + + RETURN + END SUBROUTINE find_adjacent_nodes + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + SUBROUTINE find_boundary_nodes(nn,nbed,bedn,ged2nn,nbnd,bndn) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nn + INTEGER, INTENT(IN) :: nbed + INTEGER, DIMENSION(:), INTENT(IN) :: bedn + INTEGER, DIMENSION(:,:), INTENT(IN) :: ged2nn + INTEGER, INTENT(OUT) :: nbnd + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: bndn + + INTEGER :: ed,i + INTEGER :: ged,nd + + INTEGER, DIMENSION(:), ALLOCATABLE :: nd_flag + + ALLOCATE(bndn(2*nbed),nd_flag(nn)) + nd_flag = 0 + nbnd = 0 + DO ed = 1,nbed + ged = bedn(ed) + + DO i = 1,2 + nd = ged2nn(i,ged) + IF (nd_flag(nd) == 0) THEN + nbnd = nbnd + 1 + bndn(nbnd) = nd + nd_flag(nd) = 1 + ENDIF + ENDDO + + ENDDO + + RETURN + END SUBROUTINE find_boundary_nodes +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE min_edge_length(ned,ne,xy,ged2el,ged2nn,edlen,edlen_min) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ned + INTEGER, INTENT(IN) :: ne + REAL(rp), DIMENSION(:,:), INTENT(IN) :: xy + INTEGER, DIMENSION(:,:), INTENT(IN) :: ged2el + INTEGER, DIMENSION(:,:), INTENT(IN) :: ged2nn + REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: edlen + REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: edlen_min + + INTEGER :: ed + INTEGER :: el_in,el_ex + INTEGER :: n1,n2 + REAL(rp) :: x1,x2,y1,y2 + INTEGER :: alloc_status + + + ALLOCATE(edlen(ned),edlen_min(ne), STAT=alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + + + edlen_min = 1d10 + DO ed = 1,ned + n1 = ged2nn(1,ed) + n2 = ged2nn(2,ed) + + x1 = xy(1,n1) + x2 = xy(1,n2) + + y1 = xy(2,n1) + y2 = xy(2,n2) + + edlen(ed) = sqrt((x2-x1)**2 + (y2-y1)**2) + + el_in = ged2el(1,ed) + el_ex = ged2el(2,ed) + + IF (edlen(ed) < edlen_min(el_in)) THEN + edlen_min(el_in) = edlen(ed) + ENDIF + + IF (el_ex > 0 ) THEN + IF (edlen(ed) < edlen_min(el_ex)) THEN + edlen_min(el_ex) = edlen(ed) + ENDIF + ENDIF + + + ENDDO + + END SUBROUTINE min_edge_length + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE print_connect_info(mnepn,ned,nied,nobed,nfbed,nnfbed,nred) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: mnepn + INTEGER, INTENT(IN) :: ned + INTEGER, INTENT(IN) :: nied + INTEGER, INTENT(IN) :: nobed + INTEGER, INTENT(IN) :: nfbed + INTEGER, INTENT(IN) :: nnfbed + INTEGER, INTENT(IN) :: nred + + + + PRINT "(A)", "---------------------------------------------" + PRINT "(A)", " Edge Connectivity Information " + PRINT "(A)", "---------------------------------------------" + PRINT "(A)", " " + PRINT "(A,I7)", ' maximum elements per node:', mnepn + PRINT "(A)", ' ' + PRINT "(A,I7)", ' number of total edges:', ned + PRINT "(A)", ' ' + PRINT "(A,I7)", ' number of interior edges:', nied + PRINT "(A)", ' ' + PRINT "(A,I7)", ' number of open boundary edges:', nobed + PRINT "(A)", ' ' + PRINT "(A,I7)", ' number of specified normal boundary edges:', nfbed + PRINT "(A,I7)", ' number of no normal flow boundary edges:', nnfbed + PRINT "(A)", ' ' + PRINT "(A,I7)", ' number of recieve edges:', nred + PRINT "(A)", ' ' + PRINT "(A,I7)", ' number of missing edges:',ned-(nied+nobed+nfbed+nnfbed+nred) + PRINT "(A)", ' ' + + + RETURN + END SUBROUTINE print_connect_info + + END MODULE edge_connectivity_mod diff --git a/ocean/wave_mesh_tools/fix_elements.f90 b/ocean/wave_mesh_tools/fix_elements.f90 new file mode 100644 index 000000000..9cddd35b3 --- /dev/null +++ b/ocean/wave_mesh_tools/fix_elements.f90 @@ -0,0 +1,535 @@ + MODULE fix_elements + + USE globals, ONLY: rp + + IMPLICIT NONE + + CONTAINS + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE flag_problem_elements(ne,ect,nepn,el2el,keep_element) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ne + INTEGER, DIMENSION(:,:), INTENT(IN) :: ect + INTEGER, DIMENSION(:), INTENT(IN) :: nepn + INTEGER, DIMENSION(:,:), INTENT(IN) :: el2el + INTEGER, DIMENSION(:), INTENT(INOUT) :: keep_element + + INTEGER :: i,el,nd,led + INTEGER :: nneigh + INTEGER :: nconnds + + DO el = 1,ne + + ! Determine how many edge neighbors an element has + nneigh = 0 + DO led = 1,3 + IF (el2el(el,led) > 0) THEN + nneigh = nneigh + 1 + ENDIF + ENDDO + + ! Determine how many vertices are shared by more than one element + nconnds = 0 + DO i = 1,3 + nd = ect(i,el) + IF (nepn(nd) > 1) THEN + nconnds = nconnds + 1 + ENDIF + ENDDO + + ! Flag elements with no edge neighbors + IF (nneigh == 0) THEN + keep_element(el) = 0 + ENDIF + + ! Flag elements with two boundary edges and three shared verticies + IF (nneigh == 1 .and. nconnds == 3) THEN + keep_element(el) = 0 + ENDIF + + ENDDO + + RETURN + END SUBROUTINE flag_problem_elements + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE flag_single_element_passages(ne,nn,ect,nbnd,bndn,keep_element) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ne + INTEGER, INTENT(IN) :: nn + INTEGER, DIMENSION(:,:), INTENT(IN) :: ect + INTEGER :: nbnd + INTEGER, DIMENSION(:), INTENT(IN) :: bndn + INTEGER, DIMENSION(:), INTENT(INOUT) :: keep_element + + INTEGER :: el,nd + INTEGER :: bnd_count + INTEGER, DIMENSION(:), ALLOCATABLE :: bnd_flag + + ! Compute list of boundary node flags + ALLOCATE(bnd_flag(nn)) + bnd_flag = 0 + DO nd = 1,nbnd + bnd_flag(bndn(nd)) = 1 + ENDDO + + ! Flag elements with 3 boundary nodes + DO el = 1,ne + bnd_count = 0 + DO nd = 1,3 + bnd_count = bnd_count + bnd_flag(ect(nd,el)) + ENDDO + IF (bnd_count ==3) THEN + keep_element(el) = 0 + ENDIF + ENDDO + + RETURN + END SUBROUTINE flag_single_element_passages + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE flag_isolated_element_patches(ne,el2el,keep_element) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ne + INTEGER, DIMENSION(:,:), INTENT(IN) :: el2el + INTEGER, DIMENSION(:), INTENT(INOUT) :: keep_element + + INTEGER :: i,j,k,el + INTEGER :: el_next,el_check + INTEGER :: nel_found,found + INTEGER :: pos + INTEGER, DIMENSION(:), ALLOCATABLE :: el_found + + ALLOCATE(el_found(ne)) + DO el = 1,ne + + el_next = el + nel_found = 1 + pos = 1 + el_found(nel_found) = el_next + search:DO i = 1,200 + elems: DO j = 1,3 + + ! find elements connected to el + el_check = el2el(el_next,j) + + IF (el_check == 0 ) THEN + CYCLE elems + ENDIF + + IF (keep_element(el_check) == 0) THEN + CYCLE elems + ENDIF + + ! check if the element has been added to the queue + found = 0 + quen: DO k = 1,nel_found + IF (el_found(k) == el_check) THEN + found= 1 + EXIT quen + ENDIF + ENDDO quen + + ! add element to the queue if it hasn't been already + IF (found == 0 ) THEN + nel_found = nel_found + 1 + el_found(nel_found) = el_check + ENDIF + + ENDDO elems + + ! exit when all elements have been added + IF (pos == nel_found) THEN + EXIT search + ENDIF + + ! get ready to check the next element + pos = pos + 1 + el_next = el_found(pos) + + ENDDO search + + ! Flag elements + IF (nel_found < 200) THEN + DO i = 1,nel_found + keep_element(el_found(i)) = 0 + ENDDO + ENDIF + + ENDDO + + END SUBROUTINE flag_isolated_element_patches + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE fix_single_node_connections_across_islands(nn,nbnd,bndn,nbed,bedn,nepn,epn,ged2nn,nd_flag,keep_element) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nn + INTEGER, INTENT(IN) :: nbnd + INTEGER, DIMENSION(:), INTENT(IN) :: bndn + INTEGER, INTENT(IN) :: nbed + INTEGER, DIMENSION(:), INTENT(IN) :: bedn + INTEGER, DIMENSION(:), INTENT(IN) :: nepn + INTEGER, DIMENSION(:,:), INTENT(IN) :: epn + INTEGER, DIMENSION(:,:), INTENT(IN) :: ged2nn + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: nd_flag + INTEGER, DIMENSION(:), INTENT(INOUT) :: keep_element + + INTEGER :: i,j + INTEGER :: nd,ged + + ! Find nodes that are shared by more than two boundary edges + ALLOCATE(nd_flag(nn)) + nd_flag = 0 + DO i = 1,nbnd + nd = bndn(i) + DO j = 1,nbed + ged = bedn(j) + IF (ged2nn(1,ged) == nd .OR. ged2nn(2,ged) == nd) THEN + nd_flag(nd) = nd_flag(nd) + 1 + ENDIF + ENDDO + ENDDO + + ! Eliminate elements associated with these nodes + DO nd = 1,nn + IF (nd_flag(nd) > 2) THEN + DO i = 1,nepn(nd) + keep_element(epn(i,nd)) = 0 + ENDDO + ENDIF + ENDDO + + RETURN + END SUBROUTINE fix_single_node_connections_across_islands + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE flag_single_element_islands(ned,nbed,bedn,ged2nn,xy,nfill_elements,fill_elements) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ned + INTEGER, INTENT(IN) :: nbed + INTEGER, DIMENSION(:), INTENT(IN) :: bedn + INTEGER, DIMENSION(:,:), INTENT(IN) :: ged2nn + REAL(rp), DIMENSION(:,:), INTENT(IN) :: xy + INTEGER, INTENT(OUT) :: nfill_elements + INTEGER, DIMENSION(:,:), INTENT(OUT) :: fill_elements + + INTEGER :: i,j,ed,ged + INTEGER :: nd_start,nd_end,nd_next + INTEGER :: ed_count + INTEGER :: tmp + INTEGER, DIMENSION(:), ALLOCATABLE :: ed_flag + INTEGER, DIMENSION(:), ALLOCATABLE :: ed_skip + INTEGER, DIMENSION(5) :: ed_found,nd_found + REAL(rp) :: x1,x2,x3,y1,y2,y3 + REAL(rp) :: area + + ALLOCATE(ed_flag(ned)) + ALLOCATE(ed_skip(ned)) + ed_skip = 0 + nfill_elements = 0 + edge:DO ed = 1,nbed + ed_flag = 0 + ged = bedn(ed) + IF (ed_skip(ged) == 1) THEN + CYCLE edge + ENDIF + + ed_flag(ged) = 1 + nd_start = ged2nn(1,ged) + nd_end = ged2nn(2,ged) + nd_next = nd_start + nd_found = 0 + ed_found = 0 + nd_found(1) = nd_start + ed_found(1) = ged + ed_count = 1 + + ! Try to find 4 consecutive edges before wrapping around + link:DO j = 1,4 + search:DO i = 1,nbed + ged = bedn(i) + IF (ed_flag(ged) == 1) THEN + CYCLE search + ENDIF + + IF (ged2nn(1,ged) == nd_next) THEN + nd_next = ged2nn(2,ged) + ed_flag(ged) = 1 + ed_count = ed_count + 1 + nd_found(ed_count) = nd_next + ed_found(ed_count) = ged + ELSE IF (ged2nn(2,ged) == nd_next) THEN + nd_next = ged2nn(1,ged) + ed_flag(ged) = 1 + ed_count = ed_count + 1 + nd_found(ed_count) = nd_next + ed_found(ed_count) = ged + ENDIF + + IF (ed_flag(ged) == 1) THEN + EXIT search + ENDIF + IF (nd_next == nd_end) THEN + EXIT link + ENDIF + + ENDDO search + ENDDO link + + ! Fill element if only three edges/nodes were found + IF (ed_count == 3) THEN + DO i = 1,3 + ed_skip(ed_found(i)) = 1 + ENDDO + + ! Check to make sure node numbering is CCW + x1 = xy(1,nd_found(1)) + x2 = xy(1,nd_found(2)) + x3 = xy(1,nd_found(3)) + y1 = xy(2,nd_found(1)) + y2 = xy(2,nd_found(2)) + y3 = xy(2,nd_found(3)) + area = (x1-x3)*(y2-y3)+(x3-x2)*(y1-y3) + IF (area < 0d0) THEN + tmp = nd_found(2) + nd_found(2) = nd_found(3) + nd_found(3) = tmp + ENDIF + + ! Add element + nfill_elements = nfill_elements + 1 + DO i = 1,3 + fill_elements(i,nfill_elements) = nd_found(i) + ENDDO + ENDIF + + ENDDO edge + + RETURN + END SUBROUTINE flag_single_element_islands + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE create_new_ect(ne,ect,keep_element,ne_new,ect_new) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ne + INTEGER, DIMENSION(:,:), INTENT(IN) :: ect + INTEGER, DIMENSION(:), INTENT(IN) :: keep_element + INTEGER, INTENT(OUT) :: ne_new + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: ect_new + + INTEGER :: el,i + + ne_new = 0 + DO el = 1,ne + IF (keep_element(el) == 1) THEN + ne_new = ne_new + 1 + ENDIF + ENDDO + + ALLOCATE(ect_new(3,ne_new)) + ne_new = 0 + DO el = 1,ne + IF (keep_element(el) == 1) THEN + ne_new = ne_new + 1 + DO i = 1,3 + ect_new(i,ne_new) = ect(i,el) + ENDDO + ENDIF + ENDDO + PRINT*, ne,ne_new + + RETURN + END SUBROUTINE create_new_ect + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE add_elements_to_ect(ne,ect,nfill_elements,fill_elements) + + IMPLICIT NONE + + INTEGER, INTENT(INOUT) :: ne + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: ect + INTEGER, INTENT(IN) :: nfill_elements + INTEGER, DIMENSION(:,:), INTENT(IN) :: fill_elements + + INTEGER :: i,j + INTEGER, DIMENSION(:,:), ALLOCATABLE :: ect_tmp + + ALLOCATE(ect_tmp(3,2*ne)) + + DO i = 1,ne + DO j = 1,3 + ect_tmp(j,i) = ect(j,i) + ENDDO + ENDDO + + DEALLOCATE(ect) + ALLOCATE(ect(3,ne+nfill_elements)) + + DO i = 1,ne + DO j = 1,3 + ect(j,i) = ect_tmp(j,i) + ENDDO + ENDDO + + + DO i = 1,nfill_elements + DO j = 1,3 + ect(j,ne+i) = fill_elements(j,i) + ENDDO + ENDDO + PRINT*, ne,ne+nfill_elements + + ne = ne + nfill_elements + + RETURN + END SUBROUTINE add_elements_to_ect + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE remove_unconnected_nodes(nn,nepn,xy,ne,ect,nn_new,xy_new,ect_new,depth,depth_new) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nn + INTEGER, DIMENSION(:), INTENT(IN) :: nepn + REAL(rp), DIMENSION(:,:), INTENT(IN) :: xy + INTEGER, INTENT(IN) :: ne + INTEGER, DIMENSION(:,:), INTENT(IN) :: ect + INTEGER, INTENT(OUT) :: nn_new + REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: xy_new + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: ect_new + REAL(rp), DIMENSION(:), INTENT(IN), OPTIONAL :: depth + REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(OUT), OPTIONAL :: depth_new + + + INTEGER :: nd,el + INTEGER, DIMENSION(:), ALLOCATABLE :: keep_nd + INTEGER :: new_node_number + + ! Find nodes that aren't connencted to any elements + ALLOCATE(keep_nd(nn)) + nn_new = 0 + DO nd = 1,nn + IF (nepn(nd) == 0) THEN + keep_nd(nd) = 0 + ELSE + nn_new = nn_new + 1 + keep_nd(nd) = nn_new + ENDIF + ENDDO + + ! Remove nodes from node coordinate array + ALLOCATE(xy_new(2,nn_new)) + IF (PRESENT(depth)) THEN + ALLOCATE(depth_new(nn_new)) + ENDIF + DO nd = 1,nn + new_node_number = keep_nd(nd) + IF (new_node_number /= 0) THEN + xy_new(1,new_node_number) = xy(1,nd) + xy_new(2,new_node_number) = xy(2,nd) + IF (PRESENT(depth)) THEN + depth_new(new_node_number) = depth(nd) + ENDIF + ENDIF + ENDDO + + ! Renumber element connectivity table + ALLOCATE(ect_new(3,ne)) + DO el = 1,ne + DO nd = 1,3 + new_node_number = keep_nd(ect(nd,el)) + IF (new_node_number == 0) THEN + PRINT*, "Error: Node connected to element was removed" + STOP + ENDIF + ect_new(nd,el) = new_node_number + ENDDO + ENDDO + + PRINT*, nn,nn_new + + RETURN + END SUBROUTINE remove_unconnected_nodes + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE fix_element_node_orientation(ne,xy,ect) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ne + REAL(rp), DIMENSION(:,:), INTENT(IN) :: xy + INTEGER, DIMENSION(:,:), INTENT(INOUT) :: ect + + INTEGER :: el + INTEGER :: tmp + REAL(rp) :: x1,x2,x3 + REAL(rp) :: y1,y2,y3 + REAL(rp) :: x1x3,x2x1 + REAL(rp) :: area + + DO el = 1,ne + x1 = xy(1,ect(1,el)) + y1 = xy(2,ect(1,el)) + x2 = xy(1,ect(2,el)) + y2 = xy(2,ect(2,el)) + x3 = xy(1,ect(3,el)) + y3 = xy(2,ect(3,el)) + + x1x3 = x1-x3 + IF (abs(x1x3) > 180d0) THEN + x1x3 = x1x3 - SIGN(360d0,x1x3) + ENDIF + + x2x1 = x2-x1 + IF (abs(x2-x1) > 180d0) THEN + x2x1 = x2x1 - SIGN(360d0,x2x1) + ENDIF + + area = (y2-y1)*(x1x3)+(y3-y1)*(x2x1) + + IF (area < 0d0) THEN + tmp = ect(2,el) + ect(2,el) = ect(3,el) + ect(3,el) = tmp + ENDIF + + ENDDO + + RETURN + END SUBROUTINE fix_element_node_orientation + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + END MODULE fix_elements diff --git a/ocean/wave_mesh_tools/globals.f90 b/ocean/wave_mesh_tools/globals.f90 new file mode 100644 index 000000000..a15a8878d --- /dev/null +++ b/ocean/wave_mesh_tools/globals.f90 @@ -0,0 +1,16 @@ + MODULE globals + + IMPLICIT NONE + + INTEGER, PARAMETER :: rp = kind(1d0) + INTEGER, DIMENSION(4), PARAMETER :: nverts = (/3,4,3,4/) + + TYPE :: boundary_type + INTEGER :: nnds + INTEGER :: bou_type + INTEGER :: flag + REAL(rp), DIMENSION(:,:), ALLOCATABLE :: xy + INTEGER, DIMENSION(:), ALLOCATABLE :: nds + END TYPE + + END MODULE globals diff --git a/ocean/wave_mesh_tools/grid_file_mod.f90 b/ocean/wave_mesh_tools/grid_file_mod.f90 new file mode 100644 index 000000000..b6ca08296 --- /dev/null +++ b/ocean/wave_mesh_tools/grid_file_mod.f90 @@ -0,0 +1,763 @@ + MODULE grid_file_mod + + USE globals, ONLY: rp + + + IMPLICIT NONE + + TYPE :: grid_type + + CHARACTER(100) :: grid_name ! name of the grid + CHARACTER(100) :: grid_file ! name of fort.14 file + + INTEGER :: ne ! number of elements + INTEGER :: nn ! number of nodes + + INTEGER, ALLOCATABLE, DIMENSION(:) :: el_type ! element type flag + + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ect ! element connectivity table + REAL(rp), ALLOCATABLE, DIMENSION(:,:) :: xy ! x,y coordinates of nodes + REAL(rp), ALLOCATABLE, DIMENSION(:) :: depth ! depth at each node + + INTEGER :: nope ! number of open boundary segents + INTEGER, ALLOCATABLE, DIMENSION(:) :: obseg ! number of nodes in each open boundary segment + INTEGER :: neta ! total elevation specified boundary nodes + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: obnds ! open boundary nodes + + INTEGER :: nbou ! number of normal flow boundary segments + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: fbseg ! number of nodes and type of each normal flow boundary segment + INTEGER :: nvel ! total number of normal flow boundary nodes + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: fbnds ! normal flow boundary nodes + + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: inbconn ! back face node pair with front face + REAL(rp), ALLOCATABLE, DIMENSION(:,:) :: barinht ! internal barrier height + REAL(rp), ALLOCATABLE, DIMENSION(:,:) :: barincfsb ! sub-critical flow coefficient + REAL(rp), ALLOCATABLE, DIMENSION(:,:) :: barincfsp ! super-critical flow coefficient + + INTEGER :: mnepn + INTEGER, ALLOCATABLE, DIMENSION(:) :: nepn ! number of elements per node + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: epn ! elements per node + + END TYPE + + + CONTAINS + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE read_grid(mesh,option) + + IMPLICIT NONE + + TYPE(grid_type), INTENT(INOUT) :: mesh + INTEGER :: option + + + CALL read_header(0,mesh%grid_file,mesh%grid_name,mesh%ne,mesh%nn) + + CALL read_coords(mesh%nn,mesh%xy,mesh%depth) + + IF (option > 1) THEN + CALL read_connectivity(mesh%ne,mesh%ect,mesh%el_type) + ENDIF + + IF (option > 2) THEN + CALL read_open_boundaries(mesh%nope,mesh%neta,mesh%obseg,mesh%obnds) + CALL read_flow_boundaries(mesh%nbou,mesh%nvel,mesh%fbseg,mesh%fbnds, & + mesh%inbconn,mesh%barinht,mesh%barincfsb,mesh%barincfsp) + ENDIF + + CALL print_grid_info(mesh%grid_file,mesh%grid_name,mesh%ne,mesh%nn) + + RETURN + END SUBROUTINE + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE read_header(myrank,grid_file,grid_name,ne,nn) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: myrank + CHARACTER(*), INTENT(IN) :: grid_file + CHARACTER(100), INTENT(OUT) :: grid_name + INTEGER, INTENT(OUT) :: ne + INTEGER, INTENT(OUT) :: nn + + LOGICAL :: file_exists + + + ! open fort.14 grid file + INQUIRE(FILE=grid_file, EXIST=file_exists) + IF(file_exists .eqv. .FALSE.) THEN + PRINT*, "grid file does not exist" + STOP + ENDIF + + + IF (myrank == 0 ) PRINT "(A)", "Reading in grid file" + + OPEN(UNIT=14, FILE=grid_file) + + ! read in name of grid + READ(14,"(A)") grid_name + + ! read in number of elements and number of nodes + READ(14,*) ne, nn + + RETURN + END SUBROUTINE read_header + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE read_coords(nn,xy,depth,h0) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nn + REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: xy + REAL(rp), DIMENSION(:) , ALLOCATABLE, INTENT(OUT) :: depth + REAL(rp), INTENT(IN), OPTIONAL :: h0 + + + INTEGER :: i,j + INTEGER :: limit_depth + INTEGER :: alloc_status + alloc_status = 0 + + ALLOCATE(xy(2,nn),depth(nn), STAT=alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + ! read in node coordinates and depths + DO i = 1,nn + READ(14,*) j, xy(1,j), xy(2,j), depth(j) + ENDDO + + + limit_depth = 0 + IF (PRESENT(h0)) THEN + limit_depth = 1 + ENDIF + + IF (limit_depth == 1) THEN + DO i = 1,nn + IF (depth(i) < h0) THEN + depth(i) = h0 + ENDIF + ENDDO + ENDIF + +! PRINT "(A)", "Node coordinates and depth: " +! DO i = 1,nn +! PRINT "(I5,3(F11.3,3x))", i,xy(1,i), xy(2,i), depth(i) +! ENDDO +! PRINT*, " " + + + + RETURN + END SUBROUTINE read_coords + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE read_connectivity(ne,ect,el_type) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ne + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: ect + INTEGER, DIMENSION(:) , ALLOCATABLE, INTENT(OUT) :: el_type + + INTEGER :: i,j + INTEGER :: el + INTEGER :: nv + INTEGER :: alloc_status + + + ALLOCATE(ect(4,ne),el_type(ne), STAT=alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + ! read in element connectivity + DO i = 1,ne + READ(14,*) el,nv,(ect(j,el),j = 1,nv) + + IF (nv == 3) THEN + el_type(el) = 1 + ELSE IF (nv == 4) THEN + el_type(el) = 2 + ELSE + PRINT*, "Element type not supported" + STOP + ENDIF + + ENDDO + +! PRINT "(A)", "Element connectivity table: " +! DO i = 1,ne +! PRINT "(2(I5,3x),8x,4(I5,3x))", i,nelnds(i),(ect(j,i),j=1,nelnds(i)) +! ENDDO +! PRINT*, " " + + RETURN + END SUBROUTINE read_connectivity + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE read_open_boundaries(nope,neta,obseg,obnds) + + IMPLICIT NONE + + INTEGER, INTENT(OUT) :: nope + INTEGER, INTENT(OUT) :: neta + INTEGER, DIMENSION(:) , ALLOCATABLE, INTENT(OUT) :: obseg + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: obnds + + INTEGER :: i,j + INTEGER :: nbseg + INTEGER :: alloc_status + + + + READ(14,*) nope ! number of open boundaries + READ(14,*) neta ! number of total elevation specified boundary nodes + + ALLOCATE(obseg(nope),obnds(neta,nope) ,STAT=alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + DO i = 1,nope + READ(14,*) nbseg ! read in # of nodes in segment, boundary type + obseg(i) = nbseg + DO j = 1,nbseg + READ(14,*) obnds(j,i) ! read in open boundary node numbers + ENDDO + ENDDO + +! PRINT "(A)", "Open boundary segments:" +! DO i = 1,nope +! nbseg = obseg(i) +! PRINT "(A,I5,A,I5,A)", "Open boundary segment ",i," contains ",nbseg," nodes" +! DO j = 1,nbseg +! PRINT "(I5)",obnds(j,i) +! ENDDO +! ENDDO +! PRINT*, " " + + RETURN + END SUBROUTINE read_open_boundaries + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE read_flow_boundaries(nbou,nvel,fbseg,fbnds,inbconn,barinht,barincfsb,barincfsp) + + IMPLICIT NONE + + INTEGER, INTENT(OUT) :: nbou + INTEGER, INTENT(OUT) :: nvel + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: fbseg + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: fbnds + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT), OPTIONAL :: inbconn + REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(OUT), OPTIONAL :: barinht + REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(OUT), OPTIONAL :: barincfsb + REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(OUT), OPTIONAL :: barincfsp + + INTEGER :: i,j + INTEGER :: nbseg,btype + INTEGER :: alloc_status + + + READ(14,*) nbou ! number of normal flow boundaries + READ(14,*) nvel ! total number of normal flow nodes + + ALLOCATE(fbseg(2,nbou),fbnds(nvel,nbou), STAT=alloc_status) + IF (alloc_status /= 0) THEN + PRINT*, "Allocation error" + STOP + ENDIF + + DO i = 1,nbou + + READ(14,*) nbseg, btype ! read in # of nodes in segment, boundary type + fbseg(1,i) = nbseg + fbseg(2,i) = btype + + IF (btype == 0 .OR. btype == 1 .OR. btype == 2 .OR. & + btype == 10 .OR. btype == 11 .OR. btype == 12 .OR. & + btype == 20 .OR. btype == 21 .OR. btype == 22 .OR. & + btype == 30 ) THEN + + DO j = 1,nbseg + READ(14,*) fbnds(j,i) ! read in normal flow boundary node numbers + ENDDO + ENDIF + + IF (btype == 4 .OR. btype == 24) THEN + IF (.NOT.(ALLOCATED(inbconn))) THEN + ALLOCATE(inbconn(nvel,nbou),barinht(nvel,nbou),barincfsb(nvel,nbou),barincfsp(nvel,nbou)) + ENDIF + + DO j = 1,nbseg + READ(14,*) fbnds(j,i), inbconn(j,i), barinht(j,i), barincfsb(j,i), barincfsp(j,i) + ENDDO + ENDIF + + !IF (btype == 1 .OR. btype == 11 .OR. btype == 21) THEN + ! IF (fbnds(nbseg,i) /= fbnds(1,i)) THEN + ! fbnds(nbseg+1,i) = fbnds(1,i) ! close island boundaries + ! fbseg(1,i) = fbseg(1,i) + 1 + ! ENDIF + !ENDIF + ENDDO + + CLOSE(14) + +! PRINT "(A)", "Normal flow boundary segments: " +! DO i = 1,nbou +! nbseg = fbseg(1,i) +! btype = fbseg(2,i) +! PRINT "(A,I3,A,I3,A,I5,A)", "Normal flow boundary segment ",i," type ",btype, " contains ",nbseg," nodes" +! DO j = 1,nbseg +! PRINT "(I5)", fbnds(j,i) +! ENDDO +! ENDDO +! PRINT*, " " + + RETURN + END SUBROUTINE read_flow_boundaries + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE print_grid_info(grid_file,grid_name,ne,nn) + + IMPLICIT NONE + + CHARACTER(*), INTENT(IN) :: grid_file + CHARACTER(*), INTENT(IN) :: grid_name + INTEGER, INTENT(IN) :: ne + INTEGER, INTENT(IN) :: nn + + PRINT "(A)", "---------------------------------------------" + PRINT "(A)", " Grid Information " + PRINT "(A)", "---------------------------------------------" + PRINT "(A)", " " + PRINT "(A,A)", "Grid file: ", grid_file + PRINT "(A,A)", "Grid name: ", grid_name + PRINT "(A,I15)", "Number of elements: ", ne + PRINT "(A,I15)", "Number of nodes: ", nn + PRINT*, " " + + + RETURN + END SUBROUTINE print_grid_info + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE grid_size(ne,el_type,ect,xy,h) + + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ne + INTEGER, DIMENSION(:), INTENT(IN) :: el_type + INTEGER, DIMENSION(:,:), INTENT(IN) :: ect + REAL(rp), DIMENSION(:,:), INTENT(IN) :: xy + REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: h + + INTEGER :: el,i,et + INTEGER :: n1,n2,nv + REAL(rp) :: x(4),y(4) + REAL(rp) :: psml,pl + REAL(rp) :: alpha,gamma + REAL(rp) :: area,s,r,l(4) + + + ALLOCATE(h(ne)) + + DO el = 1,ne + + et = el_type(el) + + IF (mod(et,2) == 1) THEN + nv = 3 + ELSE IF (mod(et,2) == 0) THEN + nv = 4 + ENDIF + + ! calculate edge lengths + DO i = 1,nv + n1 = mod(i+0,nv) + 1 + n2 = mod(i+1,nv) + 1 + + x(n1) = xy(1,ect(n1,el)) + x(n2) = xy(1,ect(n2,el)) + + y(n1) = xy(2,ect(n1,el)) + y(n2) = xy(2,ect(n2,el)) + + l(i) = sqrt((x(n2)-x(n1))**2 + (y(n2)-y(n1))**2) + ENDDO + + ! calculate semiperimeter + s = 0 + DO i = 1,nv + s = s + l(i) + ENDDO + s = 0.5*s + + + psml = 1d0 ! product of semiperimeter - edge lengths + pl = 1d0 ! product of edge lengths + DO i = 1,nv + psml = psml*(s-l(i)) + pl = pl*l(i) + ENDDO + + IF (mod(et,2) == 1) THEN + + ! calculate radius of smallest inscribed circle for triangles + r = sqrt(psml/s) + h(el) = 2d0*r + + ELSE IF (mod(et,2) == 0) THEN + +! !calculate radius of circle with equal area for quadrilaterals +! alpha = angle(x(4),y(4),x(1),y(1),x(2),y(2)) +! gamma = angle(x(2),y(2),x(3),y(3),x(4),y(4)) +! area = sqrt(psml-.5d0*pl*(1d0+cos(alpha+gamma))) ! Bretschneider's formula +! +! r = sqrt(area/pi) +! h(el) = 2d0*r + + h(el) = minval(l) + + ENDIF + + + ENDDO + + RETURN + END SUBROUTINE grid_size + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE boundary_orientation(nbou,fbseg,fbnds,ne,ect,xy,nepn,epn,orient) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nbou + INTEGER, DIMENSION(:,:), INTENT(IN) :: fbseg + INTEGER, DIMENSION(:,:), INTENT(IN) :: fbnds + INTEGER, INTENT(IN) :: ne + INTEGER, DIMENSION(:,:), INTENT(IN) :: ect + REAL(rp), DIMENSION(:,:), INTENT(IN) :: xy + INTEGER, DIMENSION(:), INTENT(IN) :: nepn + INTEGER, DIMENSION(:,:), INTENT(IN) :: epn + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: orient + + INTEGER :: i,j,k,bou + INTEGER :: n1,n2,n3 + INTEGER :: v1,v2,v3 + INTEGER :: el + INTEGER :: cnt + REAL(rp) :: xc,yc + REAL(rp) :: x1,y1,x2,y2,x3,y3 + REAL(rp) :: cross_product + REAL(rp), PARAMETER :: t = -0.9d0 + + ALLOCATE(orient(nbou)) + cnt = 0 + + DO bou = 1,nbou + n1 = fbnds(1,bou) + n2 = fbnds(2,bou) + +elsrch: DO k = 1,nepn(n1) + + el = epn(k,n1) + + DO i = 1,3 + v1 = mod(i+0,3) + 1 + v2 = mod(i+1,3) + 1 + v3 = mod(i+2,3) + 1 + IF ((ect(v1,el) == n1 .and. ect(v2,el) == n2) .or. & + (ect(v1,el) == n2 .and. ect(v2,el) == n1)) THEN + + n3 = ect(v3,el) + + + x1 = xy(1,n1) + y1 = xy(2,n1) + x2 = xy(1,n2) + y2 = xy(2,n2) + x3 = xy(1,n3) + y3 = xy(2,n3) + + xc = 0.5d0*(1d0-t)*x2 + 0.5*(1d0+t)*x3 + yc = 0.5d0*(1d0-t)*y2 + 0.5*(1d0+t)*y3 + + cross_product = (x2-x1)*(yc-y1) - (y2-y1)*(xc-x1) + + IF (cross_product < 0d0) THEN + orient(bou) = -1 ! CCW + ELSE + orient(bou) = 1 ! CW + ENDIF + + EXIT elsrch + ENDIF + ENDDO + + ENDDO elsrch + + IF (fbseg(2,bou) /= 30) THEN + cnt = cnt + 1 + PRINT*, cnt,bou,fbseg(2,bou),orient(bou) + ENDIF + + ENDDO + + RETURN + END SUBROUTINE boundary_orientation + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE write_grid(mesh) + + IMPLICIT NONE + + TYPE(grid_type), INTENT(IN) :: mesh + + INTEGER :: nverts(4) + + nverts = (/ 3, 4, 3, 4 /) + + CALL write_header(mesh%grid_file,mesh%grid_name,mesh%ne,mesh%nn) + CALL write_coords(mesh%nn,mesh%xy,mesh%depth) + CALL write_connectivity(mesh%ne,mesh%ect,mesh%el_type,nverts) + CALL write_open_boundaries(mesh%nope,mesh%neta,mesh%obseg,mesh%obnds) + CALL write_flow_boundaries(mesh%nbou,mesh%nvel,mesh%fbseg,mesh%fbnds, & + mesh%inbconn,mesh%barinht,mesh%barincfsb,mesh%barincfsp) + + + RETURN + END SUBROUTINE write_grid + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE write_header(grid_file,grid_name,ne,nn) + + IMPLICIT NONE + + CHARACTER(*), INTENT(IN) :: grid_file + CHARACTER(*), INTENT(IN) :: grid_name + INTEGER, INTENT(IN) :: ne + INTEGER, INTENT(IN) :: nn + + ! open fort.14 grid file + OPEN(UNIT=14, FILE=TRIM(ADJUSTL(grid_file))) + + ! write out name of grid + WRITE(14,"(A)") TRIM(ADJUSTL(grid_name)) + + ! read in number of elements and number of nodes + WRITE(14,"(2(I8,1x))") ne, nn + + RETURN + END SUBROUTINE write_header + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE write_coords(nn,xy,depth) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nn + REAL(rp), DIMENSION(:,:), INTENT(IN) :: xy + REAL(rp), DIMENSION(:) , INTENT(IN) :: depth + + INTEGER :: i + + ! write out node coordinates and depths + DO i = 1,nn +! WRITE(14,"(I8,1X,3(D24.17,1X))") i, xy(1,i), xy(2,i), depth(i) + WRITE(14,"(I8,1X,3(E24.17,1X))") i, xy(1,i), xy(2,i), depth(i) + ENDDO + + RETURN + END SUBROUTINE write_coords + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE write_connectivity(ne,ect,el_type,nverts) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ne + INTEGER, DIMENSION(:,:), INTENT(IN) :: ect + INTEGER, DIMENSION(:) , INTENT(IN) :: el_type + INTEGER, DIMENSION(:) , INTENT(IN) :: nverts + + INTEGER :: i,j + INTEGER :: et,nv + + ! write out element connectivity + DO i = 1,ne + et = el_type(i) + nv = nverts(et) + WRITE(14,"(6(I8,1x))") i,nv,(ect(j,i),j = 1,nv) + ENDDO + + + RETURN + END SUBROUTINE write_connectivity + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE write_open_boundaries(nope,neta,obseg,obnds) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nope + INTEGER, INTENT(IN) :: neta + INTEGER, DIMENSION(:) , INTENT(IN) :: obseg + INTEGER, DIMENSION(:,:), INTENT(IN) :: obnds + + INTEGER :: i,j + INTEGER :: nbseg + + WRITE(14,"(I8,19x,A)") nope, "! number of open boundaries" + WRITE(14,"(I8,19x,A)") neta, "! number of total elevation specified boundary nodes" + + DO i = 1,nope + nbseg = obseg(i) + WRITE(14,"(I8,1x,I8,10x,A,1x,I8)") nbseg,0, "! number of nodes in open boundary", i + DO j = 1,nbseg + WRITE(14,"(I8)") obnds(j,i) ! write out open boundary node numbers + ENDDO + ENDDO + + RETURN + END SUBROUTINE write_open_boundaries + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + SUBROUTINE write_flow_boundaries(nbou,nvel,fbseg,fbnds,inbconn,barinht,barincfsb,barincfsp) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nbou + INTEGER, INTENT(IN) :: nvel + INTEGER, DIMENSION(:,:), INTENT(IN) :: fbseg + INTEGER, DIMENSION(:,:), INTENT(IN) :: fbnds + INTEGER, DIMENSION(:,:), INTENT(IN), OPTIONAL :: inbconn + REAL(rp), DIMENSION(:,:), INTENT(IN), OPTIONAL :: barinht + REAL(rp), DIMENSION(:,:), INTENT(IN), OPTIONAL :: barincfsb + REAL(rp), DIMENSION(:,:), INTENT(IN), OPTIONAL :: barincfsp + + INTEGER :: i,j + INTEGER :: nbseg,btype + + WRITE(14,"(I8,19x,A)") nbou, "! number of normal flow boundaries" + WRITE(14,"(I8,19x,A)") nvel, "! total number of normal flow nodes" + + DO i = 1,nbou + + nbseg = fbseg(1,i) + btype = fbseg(2,i) + WRITE(14,"(I8,1x,I8,10x,A,1x,I8)") nbseg, btype, "! number of nodes in normal flow boundary", i + + IF (btype == 0 .OR. btype == 1 .OR. btype == 2 .OR. & + btype == 10 .OR. btype == 11 .OR. btype == 12 .OR. & + btype == 20 .OR. btype == 21 .OR. btype == 22 ) THEN + + DO j = 1,nbseg + WRITE(14,"(I8)") fbnds(j,i) ! write out normal flow boundary node numbers + ENDDO + + ELSE IF (btype == 4 .OR. btype == 24) THEN + + DO j = 1,nbseg + WRITE(14,"(I8,1x,I8,1x,3(E24.17,1X))") fbnds(j,i), inbconn(j,i), barinht(j,i), barincfsb(j,i), barincfsp(j,i) + ENDDO + ENDIF + + ENDDO + + CLOSE(14) + + RETURN + END SUBROUTINE write_flow_boundaries + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE write_open_boundary_coords(nope,obseg,obnds,xy) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nope + INTEGER, DIMENSION(:), INTENT(IN) :: obseg + INTEGER, DIMENSION(:,:), INTENT(IN) :: obnds + REAL(rp), DIMENSION(:,:), INTENT(IN) :: xy + + INTEGER :: bou,nd + + OPEN(UNIT=25,FILE='obou.d') + DO bou = 1,nope + DO nd = 1,obseg(bou) + WRITE(25,"(F12.6,1x,F12.6)") xy(2,obnds(nd,bou)), xy(1,obnds(nd,bou)) + ENDDO + ENDDO + CLOSE(25) + + RETURN + END SUBROUTINE write_open_boundary_coords + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + END MODULE grid_file_mod diff --git a/ocean/wave_mesh_tools/in_cell_mod.f90 b/ocean/wave_mesh_tools/in_cell_mod.f90 new file mode 100644 index 000000000..93388d0c4 --- /dev/null +++ b/ocean/wave_mesh_tools/in_cell_mod.f90 @@ -0,0 +1,217 @@ + MODULE in_cell_mod + + USE globals, ONLY: rp + USE kdtree2_module + USE netcdf + + IMPLICIT NONE + + INTEGER :: ocean_ncid + INTEGER :: nCells_dimid, nVertices_dimid, maxEdges_dimid + INTEGER :: cellsOnVertex_varid, lonCell_varid, latCell_varid + INTEGER :: verticesOnCell_varid, nEdgesOnCell_varid + INTEGER :: lonVertex_varid, latVertex_varid + + INTEGER :: nCells, nVertices + INTEGER :: maxEdges + INTEGER, DIMENSION(:,:), ALLOCATABLE :: verticesOnCell + INTEGER, DIMENSION(:), ALLOCATABLE :: nEdgesOnCell + REAL(rp), DIMENSION(:), ALLOCATABLE :: lonCell, latCell + REAL(rp), DIMENSION(:), ALLOCATABLE :: lonVertex, latVertex + REAL(rp), DIMENSION(:,:), ALLOCATABLE :: xyzCell + + REAL(rp), DIMENSION(:,:), ALLOCATABLE :: lonlatCells + REAL(rp), DIMENSION(:), ALLOCATABLE :: area + + TYPE(kdtree2), POINTER :: tree_xy + TYPE(kdtree2_result), ALLOCATABLE, DIMENSION(:) :: closest + INTEGER :: srchdp + REAL(rp), PARAMETER :: tol = 1d-12 + REAL(rp), PARAMETER :: pi = 4d0*atan(1d0) + REAL(rp), PARAMETER :: R = 6371d0 + + CONTAINS + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE in_cell_init(ocean_mesh_file) + + IMPLICIT NONE + + CHARACTER(*), INTENT(IN) :: ocean_mesh_file + + INTEGER :: i + + PRINT*, "Initializing in_cell" + + ! Open ocean mesh + CALL check(NF90_OPEN(ocean_mesh_file, NF90_NOWRITE, ocean_ncid)) + + ! Get dimension IDs + CALL check(NF90_INQ_DIMID(ocean_ncid, 'nCells', nCells_dimid)) + CALL check(NF90_INQ_DIMID(ocean_ncid, 'nVertices', nVertices_dimid)) + CALL check(NF90_INQ_DIMID(ocean_ncid, 'maxEdges', maxEdges_dimid)) + + ! Get dimension values for ocean mesh + CALL check(NF90_INQUIRE_DIMENSION(ocean_ncid, nCells_dimid, len=nCells)) + CALL check(NF90_INQUIRE_DIMENSION(ocean_ncid, nVertices_dimid, len=nVertices)) + CALL check(NF90_INQUIRE_DIMENSION(ocean_ncid, maxEdges_dimid, len=maxEdges)) + + ! Get variables IDs + CALL check(NF90_INQ_VARID(ocean_ncid, 'verticesOnCell', verticesOnCell_varid)) + CALL check(NF90_INQ_VARID(ocean_ncid, 'nEdgesOnCell', nEdgesOnCell_varid)) + CALL check(NF90_INQ_VARID(ocean_ncid, 'lonCell', lonCell_varid)) + CALL check(NF90_INQ_VARID(ocean_ncid, 'latCell', latCell_varid)) + CALL check(NF90_INQ_VARID(ocean_ncid, 'lonVertex', lonVertex_varid)) + CALL check(NF90_INQ_VARID(ocean_ncid, 'latVertex', latVertex_varid)) + + ! Get variable values for oean mesh + ALLOCATE(verticesOnCell(maxEdges, nCells)) + CALL check(NF90_GET_VAR(ocean_ncid, verticesOnCell_varid, verticesOnCell)) + ALLOCATE(nEdgesOnCell(nCells)) + CALL check(NF90_GET_VAR(ocean_ncid, nEdgesOnCell_varid, nEdgesOnCell)) + ALLOCATE(lonCell(nCells),latCell(nCells)) + CALL check(NF90_GET_VAR(ocean_ncid, lonCell_varid, lonCell)) + CALL check(NF90_GET_VAR(ocean_ncid, latCell_varid, latCell)) + ALLOCATE(lonVertex(nVertices),latVertex(nVertices)) + CALL check(NF90_GET_VAR(ocean_ncid, lonVertex_varid, lonVertex)) + CALL check(NF90_GET_VAR(ocean_ncid, latVertex_varid, latVertex)) + CALL check(NF90_CLOSE(ocean_ncid)) + + ! Initialize kd-tree based on ocean cell centers + srchdp = 20 + ALLOCATE(closest(srchdp)) + ALLOCATE(xyzCell(3,nCells)) + DO i = 1,nCells + CALL lonlat2xyz(lonCell(i),latCell(i),R,xyzCell(1,i),xyzCell(2,i),xyzCell(3,i)) + ENDDO + tree_xy => kdtree2_create(xyzCell(1:3,1:nCells), rearrange=.true., sort=.true.) + + ! Compute ocean cell areas + ALLOCATE(area(nCells)) + DO i = 1,nCells + CALL compute_total_area(i,lonCell(i),latCell(i),area(i)) + ENDDO + + PRINT*, "Done" + + RETURN + END SUBROUTINE in_cell_init + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE pt_in_cell(lonlat,in_cell) + + IMPLICIT NONE + + REAL(rp), INTENT(IN) :: lonlat(2) + INTEGER, INTENT(OUT) :: in_cell + + INTEGER :: i + INTEGER :: cell + REAL(rp) :: area_sum + REAL(rp) :: xyz(3) + + CALL lonlat2xyz(lonlat(1),lonlat(2),R,xyz(1),xyz(2),xyz(3)) + + CALL kdtree2_n_nearest(tp=tree_xy,qv=xyz,nn=srchdp,results=closest) + + in_cell = 0 +srch: DO i = 1,srchdp + + cell = closest(i)%idx + CALL compute_total_area(cell,lonlat(1),lonlat(2),area_sum) + + IF (abs(area_sum-area(cell)) < tol) THEN + in_cell = cell + EXIT srch + ENDIF + + ENDDO srch + + RETURN + END SUBROUTINE pt_in_cell + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE compute_total_area(cell,x3,y3,total_area) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: cell + REAL(rp), INTENT(IN) :: x3 + REAL(rp), INTENT(IN) :: y3 + REAL(rp), INTENT(OUT) :: total_area + + INTEGER :: j + INTEGER :: n1,n2 + REAL(rp) :: x1,x2,x2x1,x3x1 + REAL(rp) :: y1,y2 + + total_area = 0d0 + DO j = 1,nEdgesOnCell(cell) + n1 = mod(j+0,nEdgesOnCell(cell))+1 + n2 = mod(j+1,nEdgesOnCell(cell))+1 + + x1 = lonVertex(verticesOnCell(n1,cell)) + y1 = latVertex(verticesOnCell(n1,cell)) + + x2 = lonVertex(verticesOnCell(n2,cell)) + y2 = latVertex(verticesOnCell(n2,cell)) + + x2x1 = x2-x1 + IF (abs(x2x1) > pi) THEN + x2x1 = x2x1 - SIGN(2d0*pi,x2x1) + ENDIF + x3x1 = x3-x1 + IF (abs(x3x1) > pi) THEN + x3x1 = x3x1 - SIGN(2d0*pi,x3x1) + ENDIF + + total_area = total_area + 0.5d0*abs((x2x1)*(y3-y1) - (x3x1)*(y2-y1)) + ENDDO + + RETURN + END SUBROUTINE compute_total_area + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE check(status) + + USE netcdf + + IMPLICIT NONE + INTEGER :: status + + IF (status /= NF90_NOERR) THEN + PRINT("(A, A)"), "fatal error from ", TRIM(NF90_STRERROR(status)) + ENDIF + + RETURN + END SUBROUTINE check + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE lonlat2xyz(lon,lat,R,x,y,z) + + IMPLICIT NONE + + REAL(rp), INTENT(IN) :: lon, lat, R + REAL(rp), INTENT(OUT) :: x, y, z + + x = R*cos(lat)*cos(lon) + y = R*cos(lat)*sin(lon) + z = R*sin(lat) + + RETURN + END SUBROUTINE lonlat2xyz + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + END MODULE in_cell_mod diff --git a/ocean/wave_mesh_tools/kdtree2.f90 b/ocean/wave_mesh_tools/kdtree2.f90 new file mode 100644 index 000000000..380897a1c --- /dev/null +++ b/ocean/wave_mesh_tools/kdtree2.f90 @@ -0,0 +1,1923 @@ +! +!(c) Matthew Kennel, Institute for Nonlinear Science (2004) +! +! Licensed under the Academic Free License version 1.1 found in file LICENSE +! with additional provisions found in that same file. +! +! Modified by Chris Massey, USACE-ARMY-MIL to use ADCIRC's SIZES, SZ +! parameter for the kdkind type for consistency within ADCIRC. +! Also converted to fixed format instead of free format (f90). + + + module kdtree2_precision_module + integer, parameter :: sp = kind(0.0) + integer, parameter :: dp = kind(0.0d0) + + private :: sp, dp + + ! + ! You must comment out exactly one + ! of the two lines. If you comment + ! out kdkind = sp then you get single precision + ! and if you comment out kdkind = dp + ! you get double precision. + ! + + !integer, parameter :: kdkind = sp + integer, parameter :: kdkind = dp + public :: kdkind + + end module kdtree2_precision_module + + module kdtree2_priority_queue_module + use kdtree2_precision_module + ! + ! maintain a priority queue (PQ) of data, pairs of 'priority/payload', + ! implemented with a binary heap. This is the type, and the 'dis' field + ! is the priority. + ! + type kdtree2_result + ! a pair of distances, indexes + real(kdkind) :: dis!=0.0 + integer :: idx!=-1 Initializers cause some bugs in compilers. + end type kdtree2_result + ! + ! A heap-based priority queue lets one efficiently implement the following + ! operations, each in log(N) time, as opposed to linear time. + ! + ! 1) add a datum (push a datum onto the queue, increasing its length) + ! 2) return the priority value of the maximum priority element + ! 3) pop-off (and delete) the element with the maximum priority, decreasing + ! the size of the queue. + ! 4) replace the datum with the maximum priority with a supplied datum + ! (of either higher or lower priority), maintaining the size of the + ! queue. + ! + ! + ! In the k-d tree case, the 'priority' is the square distance of a point in + ! the data set to a reference point. The goal is to keep the smallest M + ! distances to a reference point. The tree algorithm searches terminal + ! nodes to decide whether to add points under consideration. + ! + ! A priority queue is useful here because it lets one quickly return the + ! largest distance currently existing in the list. If a new candidate + ! distance is smaller than this, then the new candidate ought to replace + ! the old candidate. In priority queue terms, this means removing the + ! highest priority element, and inserting the new one. + ! + ! Algorithms based on Cormen, Leiserson, Rivest, _Introduction + ! to Algorithms_, 1990, with further optimization by the author. + ! + ! Originally informed by a C implementation by Sriranga Veeraraghavan. + ! + ! This module is not written in the most clear way, but is implemented such + ! for speed, as it its operations will be called many times during searches + ! of large numbers of neighbors. + ! + type pq + ! + ! The priority queue consists of elements + ! priority(1:heap_size), with associated payload(:). + ! + ! There are heap_size active elements. + ! Assumes the allocation is always sufficient. Will NOT increase it + ! to match. + integer :: heap_size = 0 + type(kdtree2_result), pointer :: elems(:) + end type pq + + public :: kdtree2_result + + public :: pq + public :: pq_create + public :: pq_delete, pq_insert + public :: pq_extract_max, pq_max, pq_replace_max, pq_maxpri + private + + contains + + + function pq_create(results_in) result(res) + ! + ! Create a priority queue from ALREADY allocated + ! array pointers for storage. NOTE! It will NOT + ! add any alements to the heap, i.e. any existing + ! data in the input arrays will NOT be used and may + ! be overwritten. + ! + ! usage: + ! real(kdkind), pointer :: x(:) + ! integer, pointer :: k(:) + ! allocate(x(1000),k(1000)) + ! pq => pq_create(x,k) + ! + type(kdtree2_result), target:: results_in(:) + type(pq) :: res + ! + ! + integer :: nalloc + + nalloc = size(results_in,1) + if (nalloc .lt. 1) then + write (*,*) 'PQ_CREATE: error, input arrays must be allocated.' + end if + res%elems => results_in + res%heap_size = 0 + return + end function pq_create + + ! + ! operations for getting parents and left + right children + ! of elements in a binary heap. + ! + + ! + ! These are written inline for speed. + ! + ! integer function parent(i) + ! integer, intent(in) :: i + ! parent = (i/2) + ! return + ! end function parent + + ! integer function left(i) + ! integer, intent(in) ::i + ! left = (2*i) + ! return + ! end function left + + ! integer function right(i) + ! integer, intent(in) :: i + ! right = (2*i)+1 + ! return + ! end function right + + ! logical function compare_priority(p1,p2) + ! real(kdkind), intent(in) :: p1, p2 + ! + ! compare_priority = (p1 .gt. p2) + ! return + ! end function compare_priority + + subroutine heapify(a,i_in) + ! + ! take a heap rooted at 'i' and force it to be in the + ! heap canonical form. This is performance critical + ! and has been tweaked a little to reflect this. + ! + type(pq),pointer :: a + integer, intent(in) :: i_in + ! + integer :: i, l, r, largest + + real(kdkind) :: pri_i, pri_l, pri_r, pri_largest + + type(kdtree2_result) :: temp + + i = i_in + + bigloop: do + l = 2*i ! left(i) + r = l+1 ! right(i) + ! + ! set 'largest' to the index of either i, l, r + ! depending on whose priority is largest. + ! + ! note that l or r can be larger than the heap size + ! in which case they do not count. + + + ! does left child have higher priority? + if (l .gt. a%heap_size) then + ! we know that i is the largest as both l and r are invalid. + exit + else + pri_i = a%elems(i)%dis + pri_l = a%elems(l)%dis + if (pri_l .gt. pri_i) then + largest = l + pri_largest = pri_l + else + largest = i + pri_largest = pri_i + endif + + ! + ! between i and l we have a winner + ! now choose between that and r. + ! + if (r .le. a%heap_size) then + pri_r = a%elems(r)%dis + if (pri_r .gt. pri_largest) then + largest = r + endif + endif + endif + + if (largest .ne. i) then + ! swap data in nodes largest and i, then heapify + + temp = a%elems(i) + a%elems(i) = a%elems(largest) + a%elems(largest) = temp + ! + ! Canonical heapify() algorithm has tail-ecursive call: + ! + ! call heapify(a,largest) + ! we will simulate with cycle + ! + i = largest + cycle bigloop ! continue the loop + else + return ! break from the loop + end if + enddo bigloop + return + end subroutine heapify + + subroutine pq_max(a,e) + ! + ! return the priority and its payload of the maximum priority element + ! on the queue, which should be the first one, if it is + ! in heapified form. + ! + type(pq),pointer :: a + type(kdtree2_result),intent(out) :: e + + if (a%heap_size .gt. 0) then + e = a%elems(1) + else + write (*,*) 'PQ_MAX: ERROR, heap_size < 1' + stop + endif + return + end subroutine pq_max + + real(kdkind) function pq_maxpri(a) + type(pq), pointer :: a + + if (a%heap_size .gt. 0) then + pq_maxpri = a%elems(1)%dis + else + write (*,*) 'PQ_MAX_PRI: ERROR, heapsize < 1' + stop + endif + return + end function pq_maxpri + + subroutine pq_extract_max(a,e) + ! + ! return the priority and payload of maximum priority + ! element, and remove it from the queue. + ! (equivalent to 'pop()' on a stack) + ! + type(pq),pointer :: a + type(kdtree2_result), intent(out) :: e + + if (a%heap_size .ge. 1) then + ! + ! return max as first element + ! + e = a%elems(1) + + ! + ! move last element to first + ! + a%elems(1) = a%elems(a%heap_size) + a%heap_size = a%heap_size-1 + call heapify(a,1) + return + else + write (*,*) 'PQ_EXTRACT_MAX: error,', & + ' attempted to pop non-posiive PQ' + stop + end if + + end subroutine pq_extract_max + + + real(kdkind) function pq_insert(a,dis,idx) + ! + ! Insert a new element and return the new maximum priority, + ! which may or may not be the same as the old maximum priority. + ! + type(pq),pointer :: a + real(kdkind), intent(in) :: dis + integer, intent(in) :: idx + ! type(kdtree2_result), intent(in) :: e + ! + integer :: i, isparent + real(kdkind) :: parentdis + ! + + ! if (a%heap_size .ge. a%max_elems) then + ! write (*,*) 'PQ_INSERT: error, attempt made to insert element on full PQ' + ! stop + ! else + a%heap_size = a%heap_size + 1 + i = a%heap_size + + do while (i .gt. 1) + isparent = int(i/2) + parentdis = a%elems(isparent)%dis + if (dis .gt. parentdis) then + ! move what was in i's parent into i. + a%elems(i)%dis = parentdis + a%elems(i)%idx = a%elems(isparent)%idx + i = isparent + else + exit + endif + end do + + ! insert the element at the determined position + a%elems(i)%dis = dis + a%elems(i)%idx = idx + + pq_insert = a%elems(1)%dis + return + ! end if + + end function pq_insert + + subroutine pq_adjust_heap(a,i) + type(pq),pointer :: a + integer, intent(in) :: i + ! + ! nominally arguments (a,i), but specialize for a=1 + ! + ! This routine assumes that the trees with roots 2 and 3 are already heaps, i.e. + ! the children of '1' are heaps. When the procedure is completed, the + ! tree rooted at 1 is a heap. + real(kdkind) :: prichild + integer :: parent, child, N + + type(kdtree2_result) :: e + + e = a%elems(i) + + parent = i + child = 2*i + N = a%heap_size + + do while (child .le. N) + if (child .lt. N) then + if (a%elems(child)%dis .lt. a%elems(child+1)%dis) then + child = child+1 + endif + endif + prichild = a%elems(child)%dis + if (e%dis .ge. prichild) then + exit + else + ! move child into parent. + a%elems(parent) = a%elems(child) + parent = child + child = 2*parent + end if + end do + a%elems(parent) = e + return + end subroutine pq_adjust_heap + + + real(kdkind) function pq_replace_max(a,dis,idx) + ! + ! Replace the extant maximum priority element + ! in the PQ with (dis,idx). Return + ! the new maximum priority, which may be larger + ! or smaller than the old one. + ! + type(pq),pointer :: a + real(kdkind), intent(in) :: dis + integer, intent(in) :: idx + ! type(kdtree2_result), intent(in) :: e + ! not tested as well! + + integer :: parent, child, N + real(kdkind) :: prichild, prichildp1 + + type(kdtree2_result) :: etmp + + if (.true.) then + N=a%heap_size + if (N .ge. 1) then + parent =1 + child=2 + + loop: do while (child .le. N) + prichild = a%elems(child)%dis + + ! + ! posibly child+1 has higher priority, and if + ! so, get it, and increment child. + ! + + if (child .lt. N) then + prichildp1 = a%elems(child+1)%dis + if (prichild .lt. prichildp1) then + child = child+1 + prichild = prichildp1 + endif + endif + + if (dis .ge. prichild) then + exit loop + ! we have a proper place for our new element, + ! bigger than either children's priority. + else + ! move child into parent. + a%elems(parent) = a%elems(child) + parent = child + child = 2*parent + end if + end do loop + a%elems(parent)%dis = dis + a%elems(parent)%idx = idx + pq_replace_max = a%elems(1)%dis + else + a%elems(1)%dis = dis + a%elems(1)%idx = idx + pq_replace_max = dis + endif + else + ! + ! slower version using elementary pop and push operations. + ! + call pq_extract_max(a,etmp) + etmp%dis = dis + etmp%idx = idx + pq_replace_max = pq_insert(a,dis,idx) + endif + return + end function pq_replace_max + + subroutine pq_delete(a,i) + ! + ! delete item with index 'i' + ! + type(pq),pointer :: a + integer :: i + + if ((i .lt. 1) .or. (i .gt. a%heap_size)) then + write (*,*) 'PQ_DELETE: error, attempt to remove', & + ' out of bounds element.' + stop + endif + + ! swap the item to be deleted with the last element + ! and shorten heap by one. + a%elems(i) = a%elems(a%heap_size) + a%heap_size = a%heap_size - 1 + + call heapify(a,i) + + end subroutine pq_delete + + end module kdtree2_priority_queue_module + + + module kdtree2_module + use kdtree2_precision_module + use kdtree2_priority_queue_module + ! K-D tree routines in Fortran 90 by Matt Kennel. + ! Original program was written in Sather by Steve Omohundro and + ! Matt Kennel. Only the Euclidean metric is supported. + ! + ! + ! This module is identical to 'kd_tree', except that the order + ! of subscripts is reversed in the data file. + ! In otherwords for an embedding of N D-dimensional vectors, the + ! data file is here, in natural Fortran order data(1:D, 1:N) + ! because Fortran lays out columns first, + ! + ! whereas conventionally (C-style) it is data(1:N,1:D) + ! as in the original kd_tree module. + ! + !-------------DATA TYPE, CREATION, DELETION--------------------- + public :: kdkind + public :: kdtree2, kdtree2_result, tree_node + public :: kdtree2_create, kdtree2_destroy + !--------------------------------------------------------------- + !-------------------SEARCH ROUTINES----------------------------- + public :: kdtree2_n_nearest,kdtree2_n_nearest_around_point + ! Return fixed number of nearest neighbors around arbitrary vector, + ! or extant point in dataset, with decorrelation window. + ! + public :: kdtree2_r_nearest, kdtree2_r_nearest_around_point + ! Return points within a fixed ball of arb vector/extant point + ! + public :: kdtree2_sort_results + ! Sort, in order of increasing distance, rseults from above. + ! + public :: kdtree2_r_count, kdtree2_r_count_around_point + ! Count points within a fixed ball of arb vector/extant point + ! + public :: kdtree2_n_nearest_brute_force + public :: kdtree2_r_nearest_brute_force + ! brute force of kdtree2_[n|r]_nearest + !---------------------------------------------------------------- + + + integer, parameter :: bucket_size = 12 + ! The maximum number of points to keep in a terminal node. + + type interval + real(kdkind) :: lower,upper + end type interval + + type :: tree_node + ! an internal tree node + private + integer :: cut_dim + ! the dimension to cut + real(kdkind) :: cut_val + ! where to cut the dimension + real(kdkind) :: cut_val_left, cut_val_right + ! improved cutoffs knowing the spread in child boxes. + integer :: l, u + type (tree_node), pointer :: left, right + type(interval), pointer :: box(:) => null() + ! child pointers + ! Points included in this node are indexes[k] with k \in [l,u] + + + end type tree_node + + type :: kdtree2 + ! Global information about the tree, one per tree + integer :: dimen=0, n=0 + ! dimensionality and total # of points + real(kdkind), pointer :: the_data(:,:) => null() + ! pointer to the actual data array + ! + ! IMPORTANT NOTE: IT IS DIMENSIONED the_data(1:d,1:N) + ! which may be opposite of what may be conventional. + ! This is, because in Fortran, the memory layout is such that + ! the first dimension is in sequential order. Hence, with + ! (1:d,1:N), all components of the vector will be in consecutive + ! memory locations. The search time is dominated by the + ! evaluation of distances in the terminal nodes. Putting all + ! vector components in consecutive memory location improves + ! memory cache locality, and hence search speed, and may enable + ! vectorization on some processors and compilers. + + integer, pointer :: ind(:) => null() + ! permuted index into the data, so that indexes[l..u] of some + ! bucket represent the indexes of the actual points in that + ! bucket. + logical :: sort = .false. + ! do we always sort output results? + logical :: rearrange = .false. + real(kdkind), pointer :: rearranged_data(:,:) => null() + ! if (rearrange .eqv. .true.) then rearranged_data has been + ! created so that rearranged_data(:,i) = the_data(:,ind(i)), + ! permitting search to use more cache-friendly rearranged_data, at + ! some initial computation and storage cost. + type (tree_node), pointer :: root => null() + ! root pointer of the tree + end type kdtree2 + + + type :: tree_search_record + ! + ! One of these is created for each search. + ! + private + ! + ! Many fields are copied from the tree structure, in order to + ! speed up the search. + ! + integer :: dimen + integer :: nn, nfound + real(kdkind) :: ballsize + integer :: centeridx=999, correltime=9999 + ! exclude points within 'correltime' of 'centeridx', iff centeridx >= 0 + integer :: nalloc ! how much allocated for results(:)? + logical :: rearrange ! are the data rearranged or original? + ! did the # of points found overflow the storage provided? + logical :: overflow + real(kdkind), pointer :: qv(:) ! query vector + type(kdtree2_result), pointer :: results(:) ! results + type(pq) :: pq + real(kdkind), pointer :: data(:,:) ! temp pointer to data + integer, pointer :: ind(:) ! temp pointer to indexes + end type tree_search_record + + private + ! everything else is private. + + type(tree_search_record), save, target :: sr ! A GLOBAL VARIABLE for search + + contains + + function kdtree2_create(input_data,idim2,sort,rearrange) result (mr) + ! + ! create the actual tree structure, given an input array of data. + ! + ! Note, input data is input_data(1:d,1:N), NOT the other way around. + ! THIS IS THE REVERSE OF THE PREVIOUS VERSION OF THIS MODULE. + ! The reason for it is cache friendliness, improving performance. + ! + ! Optional arguments: If 'idim2' is specified, then the tree + ! will only search the first 'idim2' components + ! of input_data, otherwise, idim2 is inferred + ! from SIZE(input_data,1). + ! + ! if sort .eqv. .true. then output results + ! will be sorted by increasing distance. + ! default=.false., as it is faster to not sort. + ! + ! if rearrange .eqv. .true. then an internal + ! copy of the data, rearranged by terminal node, + ! will be made for cache friendliness. + ! default=.true., as it speeds searches, but + ! building takes longer, and extra memory is used. + ! + ! .. Function Return Cut_value .. + type (kdtree2), pointer :: mr + integer, intent(in), optional :: idim2 + logical, intent(in), optional :: sort + logical, intent(in), optional :: rearrange + ! .. + ! .. Array Arguments .. + real(kdkind), target :: input_data(:,:) + ! + integer :: i + ! .. + allocate (mr) + mr%the_data => input_data + ! pointer assignment + + if (present(idim2)) then + mr%dimen = idim2 + else + mr%dimen = size(input_data,1) + end if + mr%n = size(input_data,2) + + if (mr%dimen > mr%n) then + ! unlikely to be correct + write (*,*) 'KD_TREE_TRANS: likely user error.' + write (*,*) 'KD_TREE_TRANS: You passed in matrix with D=', & + mr%dimen + write (*,*) 'KD_TREE_TRANS: and N=',mr%n + + write (*,*) 'KD_TREE_TRANS: note, that new format is', & + ' data(1:D,1:N)' + write (*,*) 'KD_TREE_TRANS: with usually N >> D. ', & + 'If N =approx= D, then a k-d tree' + write (*,*) 'KD_TREE_TRANS: is not an appropriate data', & + ' structure.' + stop + end if + + call build_tree(mr) + + if (present(sort)) then + mr%sort = sort + else + mr%sort = .false. + endif + + if (present(rearrange)) then + mr%rearrange = rearrange + else + mr%rearrange = .true. + endif + + if (mr%rearrange) then + allocate(mr%rearranged_data(mr%dimen,mr%n)) + do i=1,mr%n + mr%rearranged_data(:,i) = mr%the_data(:, & + mr%ind(i)) + enddo + else + nullify(mr%rearranged_data) + endif + + end function kdtree2_create + + subroutine build_tree(tp) + type (kdtree2), pointer :: tp + ! .. + integer :: j + type(tree_node), pointer :: dummy => null() + ! .. + allocate (tp%ind(tp%n)) + forall (j=1:tp%n) + tp%ind(j) = j + end forall + tp%root => build_tree_for_range(tp,1,tp%n, dummy) + end subroutine build_tree + + recursive function build_tree_for_range(tp,l,u,parent) result (res) + ! .. Function Return Cut_value .. + type (tree_node), pointer :: res + ! .. + ! .. Structure Arguments .. + type (kdtree2), pointer :: tp + type (tree_node),pointer :: parent + ! .. + ! .. Scalar Arguments .. + integer, intent (In) :: l, u + ! .. + ! .. Local Scalars .. + integer :: i, c, m, dimen + logical :: recompute + real(kdkind) :: average + + !!$ If (.False.) Then + !!$ If ((l .Lt. 1) .Or. (l .Gt. tp%n)) Then + !!$ Stop 'illegal L value in build_tree_for_range' + !!$ End If + !!$ If ((u .Lt. 1) .Or. (u .Gt. tp%n)) Then + !!$ Stop 'illegal u value in build_tree_for_range' + !!$ End If + !!$ If (u .Lt. l) Then + !!$ Stop 'U is less than L, thats illegal.' + !!$ End If + !!$ Endif + !!$ + ! first compute min and max + dimen = tp%dimen + allocate (res) + allocate(res%box(dimen)) + + ! First, compute an APPROXIMATE bounding box of all points associated with this node. + if ( u < l ) then + ! no points in this box + nullify(res) + return + end if + + if ((u-l)<=bucket_size) then + ! + ! always compute true bounding box for terminal nodes. + ! + do i=1,dimen + call spread_in_coordinate(tp,i,l,u,res%box(i)) + end do + res%cut_dim = 0 + res%cut_val = 0.0 + res%l = l + res%u = u + res%left =>null() + res%right => null() + else + ! + ! modify approximate bounding box. This will be an + ! overestimate of the true bounding box, as we are only recomputing + ! the bounding box for the dimension that the parent split on. + ! + ! Going to a true bounding box computation would significantly + ! increase the time necessary to build the tree, and usually + ! has only a very small difference. This box is not used + ! for searching but only for deciding which coordinate to split on. + ! + do i=1,dimen + recompute=.true. + if (associated(parent)) then + if (i .ne. parent%cut_dim) then + recompute=.false. + end if + endif + if (recompute) then + call spread_in_coordinate(tp,i,l,u,res%box(i)) + else + res%box(i) = parent%box(i) + endif + end do + + + c = maxloc(res%box(1:dimen)%upper-res%box(1:dimen)%lower,1) + ! + ! c is the identity of which coordinate has the greatest spread. + ! + + if (.false.) then + ! select exact median to have fully balanced tree. + m = (l+u)/2 + call select_on_coordinate(tp%the_data,tp%ind,c,m,l,u) + else + ! + ! select point halfway between min and max, as per A. Moore, + ! who says this helps in some degenerate cases, or + ! actual arithmetic average. + ! + if (.true.) then + ! actually compute average + average = sum(tp%the_data(c,tp%ind(l:u))) / real(u-l+1,kdkind) + else + average = (res%box(c)%upper + res%box(c)%lower)/2.0 + endif + + res%cut_val = average + m = select_on_coordinate_value( & + tp%the_data,tp%ind,c,average,l,u) + endif + + ! moves indexes around + res%cut_dim = c + res%l = l + res%u = u + ! res%cut_val = tp%the_data(c,tp%ind(m)) + + res%left => build_tree_for_range(tp,l,m,res) + res%right => build_tree_for_range(tp,m+1,u,res) + + if (associated(res%right) .eqv. .false.) then + res%box = res%left%box + res%cut_val_left = res%left%box(c)%upper + res%cut_val = res%cut_val_left + elseif (associated(res%left) .eqv. .false.) then + res%box = res%right%box + res%cut_val_right = res%right%box(c)%lower + res%cut_val = res%cut_val_right + else + res%cut_val_right = res%right%box(c)%lower + res%cut_val_left = res%left%box(c)%upper + res%cut_val = (res%cut_val_left + res%cut_val_right)/2 + + + ! now remake the true bounding box for self. + ! Since we are taking unions (in effect) of a tree structure, + ! this is much faster than doing an exhaustive + ! search over all points + res%box%upper = max(res%left%box%upper,res%right%box%upper) + res%box%lower = min(res%left%box%lower,res%right%box%lower) + endif + end if + end function build_tree_for_range + + integer function select_on_coordinate_value(v,ind,c,alpha,li,ui) & + result(res) + ! Move elts of ind around between l and u, so that all points + ! <= than alpha (in c cooordinate) are first, and then + ! all points > alpha are second. + + ! + ! Algorithm (matt kennel). + ! + ! Consider the list as having three parts: on the left, + ! the points known to be <= alpha. On the right, the points + ! known to be > alpha, and in the middle, the currently unknown + ! points. The algorithm is to scan the unknown points, starting + ! from the left, and swapping them so that they are added to + ! the left stack or the right stack, as appropriate. + ! + ! The algorithm finishes when the unknown stack is empty. + ! + ! .. Scalar Arguments .. + integer, intent (In) :: c, li, ui + real(kdkind), intent(in) :: alpha + ! .. + real(kdkind) :: v(1:,1:) + integer :: ind(1:) + integer :: tmp + ! .. + integer :: lb, rb + ! + ! The points known to be <= alpha are in + ! [l,lb-1] + ! + ! The points known to be > alpha are in + ! [rb+1,u]. + ! + ! Therefore we add new points into lb or + ! rb as appropriate. When lb=rb + ! we are done. We return the location of the last point <= alpha. + ! + ! + lb = li; rb = ui + + do while (lb < rb) + if ( v(c,ind(lb)) <= alpha ) then + ! it is good where it is. + lb = lb+1 + else + ! swap it with rb. + tmp = ind(lb); ind(lb) = ind(rb); ind(rb) = tmp + rb = rb-1 + endif + end do + + ! now lb .eq. ub + if (v(c,ind(lb)) <= alpha) then + res = lb + else + res = lb-1 + endif + + end function select_on_coordinate_value + + subroutine select_on_coordinate(v,ind,c,k,li,ui) + ! Move elts of ind around between l and u, so that the kth + ! element + ! is >= those below, <= those above, in the coordinate c. + ! .. Scalar Arguments .. + integer, intent (In) :: c, k, li, ui + ! .. + integer :: i, l, m, s, t, u + ! .. + real(kdkind) :: v(:,:) + integer :: ind(:) + ! .. + l = li + u = ui + do while (l=k) u = m - 1 + end do + end subroutine select_on_coordinate + + subroutine spread_in_coordinate(tp,c,l,u,interv) + ! the spread in coordinate 'c', between l and u. + ! + ! Return lower bound in 'smin', and upper in 'smax', + ! .. + ! .. Structure Arguments .. + type (kdtree2), pointer :: tp + type(interval), intent(out) :: interv + ! .. + ! .. Scalar Arguments .. + integer, intent (In) :: c, l, u + ! .. + ! .. Local Scalars .. + real(kdkind) :: last, lmax, lmin, t, smin,smax + integer :: i, ulocal + ! .. + ! .. Local Arrays .. + real(kdkind), pointer :: v(:,:) + integer, pointer :: ind(:) + ! .. + v => tp%the_data(1:,1:) + ind => tp%ind(1:) + smin = v(c,ind(l)) + smax = smin + + ulocal = u + + do i = l + 2, ulocal, 2 + lmin = v(c,ind(i-1)) + lmax = v(c,ind(i)) + if (lmin>lmax) then + t = lmin + lmin = lmax + lmax = t + end if + if (smin>lmin) smin = lmin + if (smaxlast) smin = last + if (smax qv + sr%nn = nn + sr%nfound = 0 + sr%centeridx = -1 + sr%correltime = 0 + sr%overflow = .false. + + sr%results => results + + sr%nalloc = nn ! will be checked + + sr%ind => tp%ind + sr%rearrange = tp%rearrange + if (tp%rearrange) then + sr%Data => tp%rearranged_data + else + sr%Data => tp%the_data + endif + sr%dimen = tp%dimen + + call validate_query_storage(nn) + sr%pq = pq_create(results) + + call search(tp%root) + + if (tp%sort) then + call kdtree2_sort_results(nn, results) + endif + ! deallocate(sr%pqp) + return + end subroutine kdtree2_n_nearest + + subroutine kdtree2_n_nearest_around_point(tp,idxin,correltime, & + nn,results) + ! Find the 'nn' vectors in the tree nearest to point 'idxin', + ! with correlation window 'correltime', returing results in + ! results(:), which must be pre-allocated upon entry. + type (kdtree2), pointer :: tp + integer, intent (In) :: idxin, correltime, nn + type(kdtree2_result), target :: results(:) + + allocate (sr%qv(tp%dimen)) + sr%qv = tp%the_data(:,idxin) ! copy the vector + sr%ballsize = huge(1.0) ! the largest real(kdkind) number + sr%centeridx = idxin + sr%correltime = correltime + + sr%nn = nn + sr%nfound = 0 + + sr%dimen = tp%dimen + sr%nalloc = nn + + sr%results => results + + sr%ind => tp%ind + sr%rearrange = tp%rearrange + + if (sr%rearrange) then + sr%Data => tp%rearranged_data + else + sr%Data => tp%the_data + endif + + call validate_query_storage(nn) + sr%pq = pq_create(results) + + call search(tp%root) + + if (tp%sort) then + call kdtree2_sort_results(nn, results) + endif + deallocate (sr%qv) + return + end subroutine kdtree2_n_nearest_around_point + + subroutine kdtree2_r_nearest(tp,qv,r2,nfound,nalloc,results) + ! find the nearest neighbors to point 'idxin', within SQUARED + ! Euclidean distance 'r2'. Upon ENTRY, nalloc must be the + ! size of memory allocated for results(1:nalloc). Upon + ! EXIT, nfound is the number actually found within the ball. + ! + ! Note that if nfound .gt. nalloc then more neighbors were found + ! than there were storage to store. The resulting list is NOT + ! the smallest ball inside norm r^2 + ! + ! Results are NOT sorted unless tree was created with sort option. + type (kdtree2), pointer :: tp + real(kdkind), target, intent (In) :: qv(:) + real(kdkind), intent(in) :: r2 + integer, intent(out) :: nfound + integer, intent (In) :: nalloc + type(kdtree2_result), target :: results(:) + + ! + sr%qv => qv + sr%ballsize = r2 + sr%nn = 0 ! flag for fixed ball search + sr%nfound = 0 + sr%centeridx = -1 + sr%correltime = 0 + + sr%results => results + + call validate_query_storage(nalloc) + sr%nalloc = nalloc + sr%overflow = .false. + sr%ind => tp%ind + sr%rearrange= tp%rearrange + + if (tp%rearrange) then + sr%Data => tp%rearranged_data + else + sr%Data => tp%the_data + endif + sr%dimen = tp%dimen + + ! + !sr%dsl = Huge(sr%dsl) ! set to huge positive values + !sr%il = -1 ! set to invalid indexes + ! + + call search(tp%root) + nfound = sr%nfound + if (tp%sort) then + call kdtree2_sort_results(nfound, results) + endif + + if (sr%overflow) then + write (*,*) 'KD_TREE_TRANS: warning! return from', & + ' kdtree2_r_nearest found more neighbors' + write (*,*) 'KD_TREE_TRANS: than storage was provided for.', & + ' Answer is NOT smallest ball' + write (*,*) 'KD_TREE_TRANS: with that number of neighbors!', & + ' I.e. it is wrong.' + endif + + return + end subroutine kdtree2_r_nearest + + subroutine kdtree2_r_nearest_around_point(tp,idxin,correltime,r2, & + nfound,nalloc,results) + ! + ! Like kdtree2_r_nearest, but around a point 'idxin' already existing + ! in the data set. + ! + ! Results are NOT sorted unless tree was created with sort option. + ! + type (kdtree2), pointer :: tp + integer, intent (In) :: idxin, correltime, nalloc + real(kdkind), intent(in) :: r2 + integer, intent(out) :: nfound + type(kdtree2_result), target :: results(:) + ! .. + ! .. Intrinsic Functions .. + intrinsic HUGE + ! .. + allocate (sr%qv(tp%dimen)) + sr%qv = tp%the_data(:,idxin) ! copy the vector + sr%ballsize = r2 + sr%nn = 0 ! flag for fixed r search + sr%nfound = 0 + sr%centeridx = idxin + sr%correltime = correltime + + sr%results => results + + sr%nalloc = nalloc + sr%overflow = .false. + + call validate_query_storage(nalloc) + + ! sr%dsl = HUGE(sr%dsl) ! set to huge positive values + ! sr%il = -1 ! set to invalid indexes + + sr%ind => tp%ind + sr%rearrange = tp%rearrange + + if (tp%rearrange) then + sr%Data => tp%rearranged_data + else + sr%Data => tp%the_data + endif + sr%rearrange = tp%rearrange + sr%dimen = tp%dimen + + ! + !sr%dsl = Huge(sr%dsl) ! set to huge positive values + !sr%il = -1 ! set to invalid indexes + ! + + call search(tp%root) + nfound = sr%nfound + if (tp%sort) then + call kdtree2_sort_results(nfound,results) + endif + + if (sr%overflow) then + write (*,*) 'KD_TREE_TRANS: warning! return from', & + ' kdtree2_r_nearest found more neighbors' + write (*,*) 'KD_TREE_TRANS: than storage was provided for.', & + ' Answer is NOT smallest ball' + write (*,*) 'KD_TREE_TRANS: with that number of neighbors!', & + ' I.e. it is wrong.' + endif + + deallocate (sr%qv) + return + end subroutine kdtree2_r_nearest_around_point + + function kdtree2_r_count(tp,qv,r2) result(nfound) + ! Count the number of neighbors within square distance 'r2'. + type (kdtree2), pointer :: tp + real(kdkind), target, intent (In) :: qv(:) + real(kdkind), intent(in) :: r2 + integer :: nfound + ! .. + ! .. Intrinsic Functions .. + intrinsic HUGE + ! .. + sr%qv => qv + sr%ballsize = r2 + + sr%nn = 0 ! flag for fixed r search + sr%nfound = 0 + sr%centeridx = -1 + sr%correltime = 0 + + nullify(sr%results) ! for some reason, FTN 95 chokes on '=> null()' + + sr%nalloc = 0 ! we do not allocate any storage but that's OK + ! for counting. + sr%ind => tp%ind + sr%rearrange = tp%rearrange + if (tp%rearrange) then + sr%Data => tp%rearranged_data + else + sr%Data => tp%the_data + endif + sr%dimen = tp%dimen + + ! + !sr%dsl = Huge(sr%dsl) ! set to huge positive values + !sr%il = -1 ! set to invalid indexes + ! + sr%overflow = .false. + + call search(tp%root) + + nfound = sr%nfound + + return + end function kdtree2_r_count + + function kdtree2_r_count_around_point(tp,idxin,correltime,r2) & + result(nfound) + ! Count the number of neighbors within square distance 'r2' around + ! point 'idxin' with decorrelation time 'correltime'. + ! + type (kdtree2), pointer :: tp + integer, intent (In) :: correltime, idxin + real(kdkind), intent(in) :: r2 + integer :: nfound + ! .. + ! .. + ! .. Intrinsic Functions .. + intrinsic HUGE + ! .. + allocate (sr%qv(tp%dimen)) + sr%qv = tp%the_data(:,idxin) + sr%ballsize = r2 + + sr%nn = 0 ! flag for fixed r search + sr%nfound = 0 + sr%centeridx = idxin + sr%correltime = correltime + nullify(sr%results) + + sr%nalloc = 0 ! we do not allocate any storage but that's OK + ! for counting. + + sr%ind => tp%ind + sr%rearrange = tp%rearrange + + if (sr%rearrange) then + sr%Data => tp%rearranged_data + else + sr%Data => tp%the_data + endif + sr%dimen = tp%dimen + + ! + !sr%dsl = Huge(sr%dsl) ! set to huge positive values + !sr%il = -1 ! set to invalid indexes + ! + sr%overflow = .false. + + call search(tp%root) + + nfound = sr%nfound + + return + end function kdtree2_r_count_around_point + + + subroutine validate_query_storage(n) + ! + ! make sure we have enough storage for n + ! + integer, intent(in) :: n + + if (size(sr%results,1) .lt. n) then + write (*,*) 'KD_TREE_TRANS: you did not provide enough', & + ' storage for results(1:n)' + stop + return + endif + + return + end subroutine validate_query_storage + + function square_distance(d, iv,qv) result (res) + ! distance between iv[1:n] and qv[1:n] + ! .. Function Return Value .. + ! re-implemented to improve vectorization. + real(kdkind) :: res + ! .. + ! .. + ! .. Scalar Arguments .. + integer :: d + ! .. + ! .. Array Arguments .. + real(kdkind) :: iv(:),qv(:) + ! .. + ! .. + res = sum( (iv(1:d)-qv(1:d))**2 ) + end function square_distance + + recursive subroutine search(node) + ! + ! This is the innermost core routine of the kd-tree search. Along + ! with "process_terminal_node", it is the performance bottleneck. + ! + ! This version uses a logically complete secondary search of + ! "box in bounds", whether the sear + ! + type (Tree_node), pointer :: node + ! .. + type(tree_node),pointer :: ncloser, nfarther + ! + integer :: cut_dim, i + ! .. + real(kdkind) :: qval, dis + real(kdkind) :: ballsize + real(kdkind), pointer :: qv(:) + type(interval), pointer :: box(:) + + if ((associated(node%left) .and. associated(node%right)) & + .eqv. .false.) then + ! we are on a terminal node + if (sr%nn .eq. 0) then + call process_terminal_node_fixedball(node) + else + call process_terminal_node(node) + endif + else + ! we are not on a terminal node + qv => sr%qv(1:) + cut_dim = node%cut_dim + qval = qv(cut_dim) + + if (qval < node%cut_val) then + ncloser => node%left + nfarther => node%right + dis = (node%cut_val_right - qval)**2 + ! extra = node%cut_val - qval + else + ncloser => node%right + nfarther => node%left + dis = (node%cut_val_left - qval)**2 + ! extra = qval- node%cut_val_left + endif + + if (associated(ncloser)) call search(ncloser) + + ! we may need to search the second node. + if (associated(nfarther)) then + ballsize = sr%ballsize + ! dis=extra**2 + if (dis <= ballsize) then + ! + ! we do this separately as going on the first cut dimen is often + ! a good idea. + ! note that if extra**2 < sr%ballsize, then the next + ! check will also be false. + ! + box => node%box(1:) + do i=1,sr%dimen + if (i .ne. cut_dim) then + dis = dis + dis2_from_bnd(qv(i),box(i)%lower,box(i)%upper) + if (dis > ballsize) then + return + endif + endif + end do + + ! + ! if we are still here then we need to search mroe. + ! + call search(nfarther) + endif + endif + end if + end subroutine search + + + real(kdkind) function dis2_from_bnd(x,amin,amax) result (res) + real(kdkind), intent(in) :: x, amin,amax + + if (x > amax) then + res = (x-amax)**2; + return + else + if (x < amin) then + res = (amin-x)**2; + return + else + res = 0.0 + return + endif + endif + return + end function dis2_from_bnd + + logical function box_in_search_range(node, sr) result(res) + ! + ! Return the distance from 'qv' to the CLOSEST corner of node's + ! bounding box + ! for all coordinates outside the box. Coordinates inside the box + ! contribute nothing to the distance. + ! + type (tree_node), pointer :: node + type (tree_search_record), pointer :: sr + + integer :: dimen, i + real(kdkind) :: dis, ballsize + real(kdkind) :: l, u + + dimen = sr%dimen + ballsize = sr%ballsize + dis = 0.0 + res = .true. + do i=1,dimen + l = node%box(i)%lower + u = node%box(i)%upper + dis = dis + (dis2_from_bnd(sr%qv(i),l,u)) + if (dis > ballsize) then + res = .false. + return + endif + end do + res = .true. + return + end function box_in_search_range + + + subroutine process_terminal_node(node) + ! + ! Look for actual near neighbors in 'node', and update + ! the search results on the sr data structure. + ! + type (tree_node), pointer :: node + ! + real(kdkind), pointer :: qv(:) + integer, pointer :: ind(:) + real(kdkind), pointer :: data(:,:) + ! + integer :: dimen, i, indexofi, k, centeridx, correltime + real(kdkind) :: ballsize, sd, newpri + logical :: rearrange + type(pq), pointer :: pqp + ! + ! copy values from sr to local variables + ! + ! + ! Notice, making local pointers with an EXPLICIT lower bound + ! seems to generate faster code. + ! why? I don't know. + qv => sr%qv(1:) + pqp => sr%pq + dimen = sr%dimen + ballsize = sr%ballsize + rearrange = sr%rearrange + ind => sr%ind(1:) + data => sr%Data(1:,1:) + centeridx = sr%centeridx + correltime = sr%correltime + + ! doing_correl = (centeridx >= 0) ! Do we have a decorrelation window? + ! include_point = .true. ! by default include all points + ! search through terminal bucket. + + mainloop: do i = node%l, node%u + if (rearrange) then + sd = 0.0 + do k = 1,dimen + sd = sd + (data(k,i) - qv(k))**2 + if (sd>ballsize) cycle mainloop + end do + indexofi = ind(i) ! only read it if we have not broken out + else + indexofi = ind(i) + sd = 0.0 + do k = 1,dimen + sd = sd + (data(k,indexofi) - qv(k))**2 + if (sd>ballsize) cycle mainloop + end do + endif + + if (centeridx > 0) then ! doing correlation interval? + if (abs(indexofi-centeridx) < correltime) cycle mainloop + endif + + + ! + ! two choices for any point. The list so far is either undersized, + ! or it is not. + ! + ! If it is undersized, then add the point and its distance + ! unconditionally. If the point added fills up the working + ! list then set the sr%ballsize, maximum distance bound (largest distance on + ! list) to be that distance, instead of the initialized +infinity. + ! + ! If the running list is full size, then compute the + ! distance but break out immediately if it is larger + ! than sr%ballsize, "best squared distance" (of the largest element), + ! as it cannot be a good neighbor. + ! + ! Once computed, compare to best_square distance. + ! if it is smaller, then delete the previous largest + ! element and add the new one. + + if (sr%nfound .lt. sr%nn) then + ! + ! add this point unconditionally to fill list. + ! + sr%nfound = sr%nfound +1 + newpri = pq_insert(pqp,sd,indexofi) + if (sr%nfound .eq. sr%nn) ballsize = newpri + ! we have just filled the working list. + ! put the best square distance to the maximum value + ! on the list, which is extractable from the PQ. + + else + ! + ! now, if we get here, + ! we know that the current node has a squared + ! distance smaller than the largest one on the list, and + ! belongs on the list. + ! Hence we replace that with the current one. + ! + ballsize = pq_replace_max(pqp,sd,indexofi) + endif + end do mainloop + ! + ! Reset sr variables which may have changed during loop + ! + sr%ballsize = ballsize + + end subroutine process_terminal_node + + subroutine process_terminal_node_fixedball(node) + ! + ! Look for actual near neighbors in 'node', and update + ! the search results on the sr data structure, i.e. + ! save all within a fixed ball. + ! + type (tree_node), pointer :: node + ! + real(kdkind), pointer :: qv(:) + integer, pointer :: ind(:) + real(kdkind), pointer :: data(:,:) + ! + integer :: nfound + integer :: dimen, i, indexofi, k + integer :: centeridx, correltime, nn + real(kdkind) :: ballsize, sd + logical :: rearrange + + ! + ! copy values from sr to local variables + ! + qv => sr%qv(1:) + dimen = sr%dimen + ballsize = sr%ballsize + rearrange = sr%rearrange + ind => sr%ind(1:) + data => sr%Data(1:,1:) + centeridx = sr%centeridx + correltime = sr%correltime + nn = sr%nn ! number to search for + nfound = sr%nfound + + ! search through terminal bucket. + mainloop: do i = node%l, node%u + + ! + ! two choices for any point. The list so far is either undersized, + ! or it is not. + ! + ! If it is undersized, then add the point and its distance + ! unconditionally. If the point added fills up the working + ! list then set the sr%ballsize, maximum distance bound (largest distance on + ! list) to be that distance, instead of the initialized +infinity. + ! + ! If the running list is full size, then compute the + ! distance but break out immediately if it is larger + ! than sr%ballsize, "best squared distance" (of the largest element), + ! as it cannot be a good neighbor. + ! + ! Once computed, compare to best_square distance. + ! if it is smaller, then delete the previous largest + ! element and add the new one. + + ! which index to the point do we use? + + if (rearrange) then + sd = 0.0 + do k = 1,dimen + sd = sd + (data(k,i) - qv(k))**2 + if (sd>ballsize) cycle mainloop + end do + indexofi = ind(i) ! only read it if we have not broken out + else + indexofi = ind(i) + sd = 0.0 + do k = 1,dimen + sd = sd + (data(k,indexofi) - qv(k))**2 + if (sd>ballsize) cycle mainloop + end do + endif + + if (centeridx > 0) then ! doing correlation interval? + if (abs(indexofi-centeridx) 1)then + ileft=ileft-1 + value=a(ileft); ivalue=ind(ileft) + else + value=a(iright); ivalue=ind(iright) + a(iright)=a(1); ind(iright)=ind(1) + iright=iright-1 + if (iright == 1) then + a(1)=value;ind(1)=ivalue + return + endif + endif + i=ileft + j=2*ileft + do while (j <= iright) + if(j < iright) then + if(a(j) < a(j+1)) j=j+1 + endif + if(value < a(j)) then + a(i)=a(j); ind(i)=ind(j) + i=j + j=j+j + else + j=iright+1 + endif + end do + a(i)=value; ind(i)=ivalue + end do + end subroutine heapsort + + subroutine heapsort_struct(a,n) + ! + ! Sort a(1:n) in ascending order + ! + ! + integer,intent(in) :: n + type(kdtree2_result),intent(inout) :: a(:) + + ! + ! + type(kdtree2_result) :: value ! temporary value + + integer :: i,j + integer :: ileft,iright + + ileft=n/2+1 + iright=n + + ! do i=1,n + ! ind(i)=i + ! Generate initial idum array + ! end do + + if(n.eq.1) return + + do + if(ileft > 1)then + ileft=ileft-1 + value=a(ileft) + else + value=a(iright) + a(iright)=a(1) + iright=iright-1 + if (iright == 1) then + a(1) = value + return + endif + endif + i=ileft + j=2*ileft + do while (j <= iright) + if(j < iright) then + if(a(j)%dis < a(j+1)%dis) j=j+1 + endif + if(value%dis < a(j)%dis) then + a(i)=a(j); + i=j + j=j+j + else + j=iright+1 + endif + end do + a(i)=value + end do + end subroutine heapsort_struct + + end module kdtree2_module + diff --git a/ocean/wave_mesh_tools/read_write_gmsh.f90 b/ocean/wave_mesh_tools/read_write_gmsh.f90 new file mode 100644 index 000000000..73159b2ff --- /dev/null +++ b/ocean/wave_mesh_tools/read_write_gmsh.f90 @@ -0,0 +1,194 @@ + MODULE read_write_gmsh + + USE grid_file_mod, ONLY: grid_type + USE globals, ONLY: rp + + + CONTAINS + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE read_gmsh_file(filename,mesh) + + IMPLICIT NONE + + CHARACTER(*), INTENT(IN) :: filename + TYPE(grid_type), INTENT(OUT) :: mesh + + CALL read_header(filename) + CALL read_nodes(mesh%nn,mesh%xy,mesh%depth) + CALL read_elements(mesh%ne,mesh%ect,mesh%el_type) + + RETURN + END SUBROUTINE read_gmsh_file + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE read_header(filename) + + IMPLICIT NONE + + CHARACTER(*), INTENT(IN) :: filename + + CHARACTER(20) :: comment + + OPEN(UNIT=14,FILE=filename) + READ(14,*) comment + READ(14,*) comment + READ(14,*) comment + + RETURN + END SUBROUTINE read_header + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE read_nodes(nn,xy,depth) + + IMPLICIT NONE + + INTEGER, INTENT(OUT) :: nn + REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: xy + REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: depth + + CHARACTER(20) :: comment + INTEGER :: i + INTEGER :: j + + READ(14,*) comment + READ(14,*) nn + + ALLOCATE(xy(2,nn)) + ALLOCATE(depth(nn)) + + DO i = 1,nn + READ(14,*), j, xy(1,j), xy(2,j), depth(j) + ENDDO + + READ(14,*) comment + + RETURN + END SUBROUTINE read_nodes + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE read_elements(ne,ect,el_type) + + IMPLICIT NONE + + INTEGER, INTENT(OUT) :: ne + INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: ect + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: el_type + + CHARACTER(20) :: comment + INTEGER :: i + INTEGER :: j + INTEGER :: k + INTEGER :: tag,ntags + INTEGER, DIMENSION(15) :: lookup + + lookup(15) = 1 + lookup(2) = 3 + + READ(14,*) comment + READ(14,*) ne + + ALLOCATE(ect(3,ne)) + ALLOCATE(el_type(ne)) + + DO i = 1,ne + READ(14,*) j,el_type(j),ntags,(tag,k=1,ntags),(ect(k,j),k=1,lookup(el_type(j))) + ENDDO + + READ(14,*) comment + + RETURN + END SUBROUTINE read_elements + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE write_gmsh_file(filename,mesh) + + IMPLICIT NONE + CHARACTER(*), INTENT(IN) :: filename + TYPE(grid_type), INTENT(IN) :: mesh + + CALL write_header(filename) + CALL write_nodes(mesh%nn,mesh%xy,mesh%depth) + CALL write_elements(mesh%ne,mesh%ect,mesh%el_type) + + RETURN + END SUBROUTINE write_gmsh_file + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE write_header(filename) + + IMPLICIT NONE + CHARACTER(*), INTENT(IN) :: filename + + OPEN(UNIT=14,FILE=filename) + WRITE(14,"(A)") "$MeshFormat" + WRITE(14,"(A)") "2 0 8" + WRITE(14,"(A)") "$EndMeshFormat" + + RETURN + END SUBROUTINE write_header + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE write_nodes(nn,xy,depth) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nn + REAL(rp), DIMENSION(:,:), INTENT(IN) :: xy + REAL(rp), DIMENSION(:), INTENT(IN) :: depth + + INTEGER :: i + + WRITE(14,"(A)") "$Nodes" + WRITE(14,"(I0)") nn + DO i = 1,nn + WRITE(14,"(I0,3(F16.8))") i, xy(1,i), xy(2,i), depth(i) + ENDDO + WRITE(14,"(A)") "$EndNodes" + + RETURN + END SUBROUTINE write_nodes + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE write_elements(ne,ect,el_type) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: ne + INTEGER, DIMENSION(:,:), INTENT(IN) :: ect + INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: el_type + + INTEGER :: i,j + + WRITE(14,"(A)") "$Elements" + WRITE(14,"(I0)") ne + DO i = 1,ne + WRITE(14,"(9(I0,2x))") i, 2, 3, 0, i, 0, ect(1,i), ect(2,i), ect(3,i) + ENDDO + WRITE(14,"(A)") "$EndElements" + + CLOSE(14) + + RETURN + END SUBROUTINE write_elements + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + END MODULE read_write_gmsh diff --git a/ocean/wave_mesh_tools/rotate.f90 b/ocean/wave_mesh_tools/rotate.f90 new file mode 100644 index 000000000..8926bb2c1 --- /dev/null +++ b/ocean/wave_mesh_tools/rotate.f90 @@ -0,0 +1,225 @@ + PROGRAM rotate + + USE rotate_mod + USE read_write_gmsh + + IMPLICIT NONE + + REAL(rp) :: LON_POLE,LAT_POLE + INTEGER :: NP + REAL(rp), DIMENSION(:), ALLOCATABLE :: LON,LAT + REAL(rp), DIMENSION(:), ALLOCATABLE :: LONR,LATR + REAL(rp), DIMENSION(:), ALLOCATABLE :: UVEC,VVEC + REAL(rp), DIMENSION(:), ALLOCATABLE :: UVECR,VVECR + REAL(rp), DIMENSION(:), ALLOCATABLE :: ANGLED + INTEGER :: ncid,lon_dimid,lat_dimid,t_dimid + INTEGER :: lon_varid,lat_varid + INTEGER :: u_varid,v_varid + INTEGER :: nlon,nlat,ntime + INTEGER :: lon_pole_varid,lat_pole_varid + REAL(rp), DIMENSION(:), ALLOCATABLE :: lon_nc,lat_nc + REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: u_nc,v_nc + INTEGER :: i,j,k,n + CHARACTER (:), ALLOCATABLE :: wind_file + CHARACTER (:), ALLOCATABLE :: mesh_file + CHARACTER (:), ALLOCATABLE :: mesh_file_out + TYPE(grid_type) :: grid + LOGICAL :: file_exists + + NAMELIST / inputs / LON_POLE, LAT_POLE, wind_file, mesh_file, mesh_file_out + + LON_POLE = -42.8906d0 + LAT_POLE = 72.3200d0 + wind_file = 'wnd10mx0.5.gdas.200506.ww3.nc' + mesh_file = 'mesh_shallow_4000_edit_mv_nd.msh' + mesh_file_out = 'mesh_shallow_4000_edit_mv_nd_RTD.msh' + + OPEN(UNIT=15, FILE='rotate.nml') + READ(UNIT=15, NML=inputs) + CLOSE(15) + + PRINT*, 'LON_POLE = ',LON_POLE + PRINT*, 'LAT_POLE = ',LAT_POLE + PRINT*, 'wind_file = ',wind_file + PRINT*, 'mesh_file = ',mesh_file + PRINT*, 'mesh_file_out = ',mesh_file_out + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Rotate mesh coordinates + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + CALL read_gmsh_file(mesh_file,grid) + ALLOCATE(LON(grid%nn),LAT(grid%nn)) + ALLOCATE(LONR(grid%nn),LATR(grid%nn)) + ALLOCATE(ANGLED(grid%nn)) + + LONR = 0d0 + LATR = 0d0 + ANGLED = 0d0 + DO i = 1,grid%nn + LON(i) = grid%xy(1,i) + LAT(i) = grid%xy(2,i) + ENDDO + + CALL W3LLTOEQ ( LAT , LON, LATR, LONR, & + ANGLED, LAT_POLE, LON_POLE, grid%nn ) + + DO i = 1,grid%nn + IF (LONR(i) > 180d0) THEN + grid%xy(1,i) = LONR(i) - 360d0 + ELSE + grid%xy(1,i) = LONR(i) + ENDIF + grid%xy(2,i) = LATR(i) + ENDDO + + !DO i = 1,grid%nn + ! PRINT("(I7,2(F12.6,F12.6,15x),F12.6)"), i, LON(i), LAT(i), LONR(i), LATR(i), ANGLED(i) + !ENDDO + + + CALL write_gmsh_file(mesh_file_out,grid) + + OPEN(UNIT=10,FILE='angled.d') + DO i = 1,grid%nn + WRITE(10,"(3(F12.6))") LON(i), LAT(i), ANGLED(i) + ENDDO + CLOSE(10) + + DEALLOCATE(LON,LAT) + DEALLOCATE(LONR,LATR) + DEALLOCATE(ANGLED) + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Read wind data + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + INQUIRE(FILE=wind_file,EXIST=file_exists) + IF (.not. file_exists) THEN + PRINT*, "Wind file does not exist" + STOP + ENDIF + + CALL check(NF90_OPEN(wind_file,NF90_WRITE,ncid)) + + IF (NF90_INQ_VARID(ncid,'LON_POLE',lon_pole_varid) .eq. NF90_NOERR) THEN + PRINT*, "Winds already rotated" + STOP + ENDIF + + !CALL check(NF90_INQ_DIMID(ncid,'lon',lon_dimid)) + !CALL check(NF90_INQ_DIMID(ncid,'lat',lat_dimid)) + CALL check(NF90_INQ_DIMID(ncid,'longitude',lon_dimid)) + CALL check(NF90_INQ_DIMID(ncid,'latitude',lat_dimid)) + CALL check(NF90_INQ_DIMID(ncid,'time',t_dimid)) + + CALL check(NF90_INQUIRE_DIMENSION(ncid,lon_dimid,len=nlon)) + CALL check(NF90_INQUIRE_DIMENSION(ncid,lat_dimid,len=nlat)) + CALL check(NF90_INQUIRE_DIMENSION(ncid,t_dimid,len=ntime)) + PRINT*, nlon,nlat,ntime + + CALL check(NF90_INQ_VARID(ncid,'longitude',lon_varid)) + CALL check(NF90_INQ_VARID(ncid,'latitude',lat_varid)) + !CALL check(NF90_INQ_VARID(ncid,'U_GRD_L103',u_varid)) + !CALL check(NF90_INQ_VARID(ncid,'V_GRD_L103',v_varid)) + CALL check(NF90_INQ_VARID(ncid,'u10',u_varid)) + CALL check(NF90_INQ_VARID(ncid,'v10',v_varid)) + + ALLOCATE(lon_nc(nlon),lat_nc(nlat)) + CALL check(NF90_GET_VAR(ncid,lon_varid,lon_nc)) + CALL check(NF90_GET_VAR(ncid,lat_varid,lat_nc)) + ALLOCATE(u_nc(nlon,nlat,ntime),v_nc(nlon,nlat,ntime)) + CALL check(NF90_GET_VAR(ncid,u_varid,u_nc)) + CALL check(NF90_GET_VAR(ncid,v_varid,v_nc)) + + NP = nlon*nlat + ALLOCATE(LON(NP), LAT(NP)) + ALLOCATE(LONR(NP), LATR(NP)) + ALLOCATE(ANGLED(NP)) + ALLOCATE(UVEC(NP),VVEC(NP)) + ALLOCATE(UVECR(NP),VVECR(NP)) + + + LONR = 0d0 + LATR = 0d0 + UVECR = 0d0 + VVECR = 0d0 + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Rotate wind coordinates + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + n = 1 + DO i = 1,nlat + DO j = 1,nlon + + IF (lon_nc(j) > 180d0) THEN + LON(n) = lon_nc(j) - 360d0 + ELSE + LON(n) = lon_nc(j) + ENDIF + + LAT(n) = lat_nc(i) + + n = n+1 + ENDDO + ENDDO + + + CALL W3LLTOEQ ( LAT , LON, LATR, LONR, & + ANGLED, LAT_POLE, LON_POLE, NP ) + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Rotate wind vectors + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + DO k = 1,ntime + PRINT*, 't = ',k + + n = 1 + DO i = 1,nlat + DO j = 1,nlon + + UVEC(n) = u_nc(j,i,k) + VVEC(n) = v_nc(j,i,k) + + n = n+1 + ENDDO + ENDDO + + UVECR = UVEC + VVECR = VVEC + CALL W3XYRTN ( NP, UVECR, VVECR, -ANGLED) + + n = 1 + DO i = 1,nlat + DO j = 1,nlon + + u_nc(j,i,k) = UVECR(n) + v_nc(j,i,k) = VVECR(n) + + if (k == 1) then + PRINT*, sqrt(UVEC(n)**2 + VVEC(n)**2), sqrt(UVECR(n)**2 + VVECR(n)**2) + endif + + n = n + 1 + ENDDO + ENDDO + + ENDDO + + + CALL check(NF90_PUT_VAR(ncid,u_varid,u_nc)) + CALL check(NF90_PUT_VAR(ncid,v_varid,v_nc)) + + CALL check(NF90_REDEF(ncid)) + CALL check(NF90_DEF_VAR(ncid,'LON_POLE',NF90_DOUBLE,lon_pole_varid)) + CALL check(NF90_DEF_VAR(ncid,'LAT_POLE',NF90_DOUBLE,lat_pole_varid)) + CALL check(NF90_ENDDEF(ncid)) + CALL check(NF90_PUT_VAR(ncid,lon_pole_varid,LON_POLE)) + CALL check(NF90_PUT_VAR(ncid,lat_pole_varid,LAT_POLE)) + + CALL check(NF90_CLOSE(ncid)) + + + END PROGRAM rotate diff --git a/ocean/wave_mesh_tools/rotate_mod.f90 b/ocean/wave_mesh_tools/rotate_mod.f90 new file mode 100644 index 000000000..d94319e05 --- /dev/null +++ b/ocean/wave_mesh_tools/rotate_mod.f90 @@ -0,0 +1,633 @@ + MODULE rotate_mod + + USE NETCDF + USE globals, ONLY: rp + + REAL (rp), PARAMETER :: PI = 4d0*atan(1d0) + REAL (rp), PARAMETER :: DERA = PI/180d0 + REAL (rp), PARAMETER :: RECIP_PI_OVER_180 = 180d0 / PI + REAL (rp), PARAMETER :: PI_OVER_180 = PI / 180d0 + REAL (rp), PARAMETER :: DEG2RAD = PI / 180d0 + REAL (rp), PARAMETER :: RAD2DEG = 180d0 / PI + REAL (rp) :: UNDEF = -999.9d0 + + CONTAINS + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! Given a radius vector (\varphi,\theta) (in rad) + ! of z', the new z after rotation, return + ! a transformation matrix: + ! RTOT: (x,y,z) -> (x',y',z') + SUBROUTINE GET_ROTMAT_ZNVEC( RTOT, VARPHIN, THETAN ) + IMPLICIT NONE + + REAL (rp):: RTOT(:,:) + REAL (rp):: VARPHIN, THETAN, PHIN, FPI + REAL (rp) :: XROT(3,3),YROT(3,3),ZROT(3,3) + + YROT(1,1) = sin( THETAN ) + YROT(2,1) = 0d0 + YROT(3,1) = -cos( THETAN ) + + YROT(1,2) = 0d0 + YROT(2,2) = 1d0 + YROT(3,2) = 0d0 + + YROT(1,3) = cos( THETAN ) + YROT(2,3) = 0d0 + YROT(3,3) = sin( THETAN ) + + !!!!!!!!!!!!!!!!!!!!!!!!!! + + PHIN = PI - VARPHIN + + ZROT(1,1) = cos( PHIN ) + ZROT(2,1) = sin( PHIN ) + ZROT(3,1) = 0d0 + + ZROT(1,2) = -sin( PHIN ) + ZROT(2,2) = cos( PHIN ) + ZROT(3,2) = 0d0 + + ZROT(1,3) = 0d0 + ZROT(2,3) = 0d0 + ZROT(3,3) = 1d0 + + RTOT = MATMUL(YROT,ZROT) + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !RTOT = 0.0_rp ; + + !RTOT(1,1) = sin( VARPHIN ) ; + !RTOT(2,1) = sin( THETAN )*cos( VARPHIN ) ; + !RTOT(3,1) = cos( THETAN )*cos( VARPHIN ) ; + + !RTOT(1,2) = -cos( VARPHIN ) ; + !RTOT(2,2) = sin( THETAN )*sin( VARPHIN ) ; + !RTOT(3,2) = cos( THETAN )*sin( VARPHIN ) ; + + !RTOT(1,3) = 0.0_rp ; + !RTOT(2,3) = -cos( THETAN ) ; + !RTOT(3,3) = sin( THETAN ) ; + + RETURN ; + END SUBROUTINE GET_ROTMAT_ZNVEC + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! + ! Transformation matrix mapping a vector from the Cartesian to spherical + ! coordinates + ! + ! [e_{r},e_{varphi},e_{theta}]^{T} = TM [ i, j, k]^{T} + FUNCTION CT2SP_VECTRANSMAT( VARPHI, THETA ) RESULT(TM) + IMPLICIT NONE + + REAL (rp) :: TM(3,3) + REAL (rp), intent(in):: VARPHI, THETA + + + TM(1,1) = cos( THETA )*cos( VARPHI ) ; + TM(2,1) = -sin( VARPHI ) ; + TM(3,1) = -sin( THETA )*cos( VARPHI ) ; + + TM(1,2) = cos( THETA )*sin( VARPHI ) ; + TM(2,2) = cos( VARPHI ) ; + TM(3,2) = -sin( THETA )*sin( VARPHI ) ; + + TM(1,3) = sin( THETA ) ; + TM(2,3) = 0.0_rp ; + TM(3,3) = cos( THETA ) ; + + + RETURN ; + END FUNCTION CT2SP_VECTRANSMAT + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE SPCOORSROTS1( RTOTS, LONR, LATR, LONO, LATO, NN ) + IMPLICIT NONE + + !c dummy c! + REAL (rp), INTENT(IN) :: RTOTS(3,3) + REAL (rp), dimension(:), INTENT(OUT) :: LONR, LATR + REAL (rp), dimension(:), INTENT(IN) :: LONO, LATO + INTEGER, INTENT(IN) :: NN + + !c local c! + INTEGER:: IP + REAL (rp):: XP(3), XPR(3), LLO, LTO + + + DO IP = 1, NN + LLO = LONO( IP ) ; + LTO = LATO( IP ) ; + + XP(1) = cos( LTO )*cos( LLO ) ; + XP(2) = cos( LTO )*sin( LLO ) ; + XP(3) = sin( LTO ) ; + + XPR = MATMUL( RTOTS, XP) ; + + LATR(IP) = atan2( XPR(3), SQRT(XPR(1)*XPR(1) & + + XPR(2)*XPR(2)) ) + LONR(IP) = atan2( XPR(2), XPR(1) ) ; + END DO + + RETURN ; + END SUBROUTINE SPCOORSROTS1 + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! (e_{r},e_{\varphi},e_{\theta})^{T} = T (i,j,k)^{T} + ! x' = R x + ! (e_{r'},e_{\varphi'},e_{\theta'})^{T} = T' (i',j',k')^{T} + ! + ! V'(r',\varphi',\theta') = T' R T^{T} V(r,\varphi,\theta) + ! + ! RVELF = T' R T^{T} + ! RVELI = (T')^{R} R^{T} T + ! + ! Evaluated the forward and inverse transformation matrices + ! between two spherical coordinate + ! RVELF: VEC(LONO,LATO) -> VEC(LONR,LATR) + ! RVELI: VEC(LONR,LATR) -> VEC(LONO,LATO) + ! + ! lonr,latr in rad + ! + SUBROUTINE SPVECROTSMAT( RVELF, RVELI, RTOTS, & + LONR, LATR, LONO, LATO, NN ) + IMPLICIT NONE + + ! Dummys ! + REAL (rp), DIMENSION(:,:,:), INTENT(OUT) :: RVELF, RVELI + + INTEGER:: NN + REAL (rp), INTENT(IN) :: RTOTS(3,3) + REAL (rp), DIMENSION(:), INTENT(IN) :: LONR, LATR, LONO, LATO + + ! local ! + INTEGER:: IP + REAL (rp), DIMENSION(3,3):: TP, TM, TRT + + + DO IP = 1, NN + TRT = 0.0_rp ; + + ! T' = T(varphi',theta') + TP = CT2SP_VECTRANSMAT( LONR(IP), LATR(IP) ) ; + + ! T = T(varphi,theta) + TM = CT2SP_VECTRANSMAT( LONO(IP), LATO(IP) ) ; + TM = TRANSPOSE( TM ) ; ! T^{T} + + ! T' R T^{T} + TRT = MATMUL( RTOTS, TM ) ; + TRT = MATMUL( TP, TRT ) ; + + RVELF(1:2,1:2,IP) = TRT(2:3,2:3) ; ! forward transform + RVELI(:,:,IP) = TRANSPOSE( RVELF(:,:,IP) ) ; ! inverse transform + END DO + + RETURN ; + END SUBROUTINE SPVECROTSMAT + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! + ! Map two vectors in spherical coordinates + ! + ! VR(:,II) = SPRV(:,:,II)*VO(:,II) + SUBROUTINE MAP2DSPVECS( VPHIR, VTHETAR, VPHIO, VTHETAO, SPRV, NN ) + IMPLICIT NONE + + ! dummy ! + INTEGER, INTENT(IN) :: NN + REAL (rp), INTENT(IN) :: SPRV(:,:,:) + REAL (rp), DIMENSION(:), INTENT(OUT) :: VPHIR, VTHETAR + REAL (rp), DIMENSION(:), INTENT(IN) :: VPHIO, VTHETAO + + ! local ! + INTEGER:: II + REAL (rp):: VO(2), VR(2) + + DO II = 1, NN + VO = (/ VPHIO(II), VTHETAO(II) /) ; + VR = MATMUL( SPRV(1:2,1:2,II), VO ) ; + + VPHIR(II) = VR(1) ; + VTHETAR(II) = VR(2) ; + END DO + + RETURN ; + END SUBROUTINE MAP2DSPVECS + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE CHECK_RTOTMAT( RTOTS, IERR ) + IMPLICIT NONE + + INTEGER:: IERR + REAL (rp), intent(in):: RTOTS(3,3) + + !c local c! + REAL (rp):: RTT(3,3), TR, SOFFD + + RTT = TRANSPOSE( RTOTS ) ; + + RTT = MATMUL(RTT, RTOTS ) ; + + ! trace of a matrix RTT + TR = RTT(1,1) + RTT(2,2) + RTT(3,3) ; + + RTT(1,1) = 0.0 ; + RTT(2,2) = 0.0 ; + RTT(3,3) = 0.0 ; + + SOFFD = SUM( MATMUL(ABS(RTT), (/ 1.0_rp, 1.0_rp, 1.0_rp /)) ) ; + + IERR = 0 ; + IF ( SOFFD > 1.0e-10_rp ) THEN + IERR = 1 ; + + PRINT*, "Error: a given rotation matrix is not orthgonal and thus invalid" ; + STOP + END IF + + RETURN ; + END SUBROUTINE CHECK_RTOTMAT + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +!Li +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!Li +!Li Merged UM source code for rotated grid, consisting the following +!Li original subroutines in UM 6.1 +!Li LLTOEQ1A WCOEFF1A and LBCROTWINDS1 +!Li The last subroutine is modified to process only one level winds +!Li cpp directives are removed and required header C_Pi.h inserted. +!Li Jian-Guo Li 26 May 2005 +!Li +!Li The WCOEFF1A subroutine is merged into LLTOEQ to reduce repetition +!Li of the same calculations. Subroutine interface changed to +!Li LLTOEQANGLE +!Li Jian-GUo Li 23 Aug 2005 +!Li +!Li Subroutine W3LLTOEQ -------------------------------------------- +!Li +!Li Purpose: Calculates latitude and longitude on equatorial +!Li latitude-longitude (eq) grid used in regional +!Li models from input arrays of latitude and +!Li longitude on standard grid. Both input and output +!Li latitudes and longitudes are in degrees. +!Li Also calculate rotation angle in degree to tranform +!Li standard wind velocity into equatorial wind. +!Li Valid for 0= 0.0) THEN + SIN_PHI_POLE = SIN(PI_OVER_180*PHI_POLE) + COS_PHI_POLE = COS(PI_OVER_180*PHI_POLE) + ELSE + SIN_PHI_POLE = -SIN(PI_OVER_180*PHI_POLE) + COS_PHI_POLE = -COS(PI_OVER_180*PHI_POLE) + ENDIF + +! 2. Transform from standard to equatorial latitude-longitude + + DO I= 1, POINTS + +! Scale longitude to range -180 to +180 degs + + A_LAMBDA=LAMBDA(I)-LAMBDA_ZERO + IF(A_LAMBDA.GT. 180.D0) A_LAMBDA=A_LAMBDA-360.D0 + IF(A_LAMBDA.LE.-180.D0) A_LAMBDA=A_LAMBDA+360.D0 + +! Convert latitude & longitude to radians + + A_LAMBDA=PI_OVER_180*A_LAMBDA + A_PHI=PI_OVER_180*PHI(I) + +! Compute eq latitude using equation (4.4) + + ARG=-COS_PHI_POLE*COS(A_PHI)*COS(A_LAMBDA) & + +SIN_PHI_POLE*SIN(A_PHI) + ARG=MIN(ARG, 1.D0) + ARG=MAX(ARG,-1.D0) + E_PHI=ASIN(ARG) + PHI_EQ(I)=RECIP_PI_OVER_180*E_PHI + +! Compute eq longitude using equation (4.6) + + TERM1 = SIN_PHI_POLE*COS(A_PHI)*COS(A_LAMBDA) & + +COS_PHI_POLE*SIN(A_PHI) + TERM2 = COS(E_PHI) + IF(TERM2 .LT. SMALL) THEN + E_LAMBDA=0.D0 + ELSE + ARG=TERM1/TERM2 + ARG=MIN(ARG, 1.D0) + ARG=MAX(ARG,-1.D0) + E_LAMBDA=RECIP_PI_OVER_180*ACOS(ARG) + E_LAMBDA=SIGN(E_LAMBDA,A_LAMBDA) + ENDIF + +! Scale longitude to range 0 to 360 degs + + IF(E_LAMBDA.GE.360.D0) E_LAMBDA=E_LAMBDA-360.D0 + IF(E_LAMBDA.LT. 0.D0) E_LAMBDA=E_LAMBDA+360.D0 + LAMBDA_EQ(I)=E_LAMBDA + +!Li Calculate turning angle for standard wind velocity + + E_LAMBDA=PI_OVER_180*LAMBDA_EQ(I) + +! Formulae used are from eqs (4.19) and (4.21) + + TERM2=SIN(E_LAMBDA) + ARG= SIN(A_LAMBDA)*TERM2*SIN_PHI_POLE & + +COS(A_LAMBDA)*COS(E_LAMBDA) + ARG=MIN(ARG, 1.D0) + ARG=MAX(ARG,-1.D0) + TERM1=RECIP_PI_OVER_180*ACOS(ARG) + ANGLED(I)=SIGN(TERM1,TERM2) +!Li + + ENDDO + +! Reset Lambda pole to the setting on entry to subroutine + LAMBDA_POLE=LAMBDA_POLE_KEEP + + RETURN + END SUBROUTINE W3LLTOEQ +!Li +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!Li +!Li Merged UM source code for rotated grid, consiting the following +!Li original subroutines in UM 6.1 +!Li EQTOLL1A WCOEFF1A and LBCROTWINDS1 +!Li The last subroutine is modified to process only one level winds +!Li cpp directives are removed and required header C_Pi.h inserted. +!Li Jian-Guo Li 26 May 2005 +!Li +!Li The WCOEFF1A subroutine is merged into EQTOLL to reduce repetition +!Li of the same calculations. Subroutine interface changed to +!Li EQTOLLANGLE +!Li First created: Jian-GUo Li 23 Aug 2005 +!Li Last modified: Jian-GUo Li 25 Feb 2008 +!Li +!Li Subroutine W3EQTOLL -------------------------------------------- +!Li +!Li Purpose: Calculates latitude and longitude on standard grid +!Li from input arrays of latitude and longitude on +!Li equatorial latitude-longitude (eq) grid used +!Li in regional models. Both input and output latitudes +!Li and longitudes are in degrees. +!Li Also calculate rotation angle in degree to tranform +!Li standard wind velocity into equatorial wind. +!Li Valid for 0= 0.0) THEN + SIN_PHI_POLE = SIN(PI_OVER_180*PHI_POLE) + COS_PHI_POLE = COS(PI_OVER_180*PHI_POLE) + ELSE + SIN_PHI_POLE = -SIN(PI_OVER_180*PHI_POLE) + COS_PHI_POLE = -COS(PI_OVER_180*PHI_POLE) + ENDIF + +! 2. Transform from equatorial to standard latitude-longitude + + DO I= 1, POINTS + +! Scale eq longitude to range -180 to +180 degs + + E_LAMBDA=LAMBDA_EQ(I) + IF(E_LAMBDA.GT. 180.0) E_LAMBDA=E_LAMBDA-360.D0 + IF(E_LAMBDA.LT.-180.0) E_LAMBDA=E_LAMBDA+360.D0 + +! Convert eq latitude & longitude to radians + + E_LAMBDA=PI_OVER_180*E_LAMBDA + E_PHI=PI_OVER_180*PHI_EQ(I) + +! Compute latitude using equation (4.7) + + ARG=COS_PHI_POLE*COS(E_PHI)*COS(E_LAMBDA) & + & +SIN_PHI_POLE*SIN(E_PHI) + ARG=MIN(ARG, 1.D0) + ARG=MAX(ARG,-1.D0) + A_PHI=ASIN(ARG) + PHI(I)=RECIP_PI_OVER_180*A_PHI + +! Compute longitude using equation (4.8) + + TERM1 = COS(E_PHI)*SIN_PHI_POLE*COS(E_LAMBDA) & + & -SIN(E_PHI)*COS_PHI_POLE + TERM2 = COS(A_PHI) + IF(TERM2.LT.SMALL) THEN + A_LAMBDA=0.D0 + ELSE + ARG=TERM1/TERM2 + ARG=MIN(ARG, 1.D0) + ARG=MAX(ARG,-1.D0) + A_LAMBDA=RECIP_PI_OVER_180*ACOS(ARG) + A_LAMBDA=SIGN(A_LAMBDA,E_LAMBDA) + A_LAMBDA=A_LAMBDA+LAMBDA_ZERO + END IF + +! Scale longitude to range 0 to 360 degs + + IF(A_LAMBDA.GE.360.0) A_LAMBDA=A_LAMBDA-360.D0 + IF(A_LAMBDA.LT. 0.0) A_LAMBDA=A_LAMBDA+360.D0 + LAMBDA(I)=A_LAMBDA + +!Li Calculate turning angle for standard wind velocity + + A_LAMBDA=PI_OVER_180*(LAMBDA(I)-LAMBDA_ZERO) + +! Formulae used are from eqs (4.19) and (4.21) + + TERM2=SIN(E_LAMBDA) + ARG=SIN(A_LAMBDA)*TERM2*SIN_PHI_POLE & + & +COS(A_LAMBDA)*COS(E_LAMBDA) + ARG=MIN(ARG, 1.D0) + ARG=MAX(ARG,-1.D0) + TERM1=RECIP_PI_OVER_180*ACOS(ARG) + ANGLED(I)=SIGN(TERM1,TERM2) +!Li + + ENDDO + + RETURN + END SUBROUTINE W3EQTOLL + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +!/ ------------------------------------------------------------------- / +!/ ------------------------------------------------------------------- / + SUBROUTINE W3XYRTN ( NSEA, XVEC, YVEC, AnglD ) +!/ +!/ +-----------------------------------+ +!/ | WAVEWATCH III NOAA/NMC | +!/ | A. Saulter | +!/ | FORTRAN 90 | +!/ | Last update : 01-Mar-2018 | +!/ +-----------------------------------+ +!/ +!/ 01-Mar-2018 : Added subroutine ( version 6.02 ) +! +! 1. Purpose : +! Subroutine to de-rotate x,y vectors from rotated to standard pole +! reference system +! +! 2. Method: +! Rotates x,y vectors anticlockwise by angle alpha in radians +! +!/ ------------------------------------------------------------------- / + IMPLICIT NONE +! +!/ ------------------------------------------------------------------- / +!/ Parameter list +!/ + INTEGER, INTENT(IN) :: NSEA ! Number of sea points + REAL(rp), INTENT(IN) :: AnglD(NSEA) ! Turning angle (degrees) + REAL(rp), INTENT(INOUT) :: XVEC(NSEA), YVEC(NSEA) +! +!/ ------------------------------------------------------------------- / +!/ Local parameters +!/ + INTEGER :: ISEA + REAL(rp) :: XVTMP, YVTMP +! +!/ ------------------------------------------------------------------- / +! Apply the rotation +! + DO ISEA=1, NSEA + IF (( XVEC(ISEA) .NE. UNDEF ) .AND. & + ( YVEC(ISEA) .NE. UNDEF )) THEN + XVTMP = XVEC(ISEA)*COS(AnglD(ISEA)*DERA) + & + YVEC(ISEA)*SIN(AnglD(ISEA)*DERA) + YVTMP = YVEC(ISEA)*COS(AnglD(ISEA)*DERA) - & + XVEC(ISEA)*SIN(AnglD(ISEA)*DERA) + XVEC(ISEA) = XVTMP + YVEC(ISEA) = YVTMP + END IF + END DO + + RETURN + END SUBROUTINE W3XYRTN + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE check(status) + + IMPLICIT NONE + INTEGER :: status + + IF(status /= NF90_NOERR) THEN + PRINT("(A,A)"), "fatal error from ", TRIM(NF90_STRERROR(status)) + ENDIF + + END SUBROUTINE check + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + END MODULE rotate_mod diff --git a/ocean/wave_mesh_tools/scrip.f90 b/ocean/wave_mesh_tools/scrip.f90 new file mode 100644 index 000000000..e67f79864 --- /dev/null +++ b/ocean/wave_mesh_tools/scrip.f90 @@ -0,0 +1,105 @@ +PROGRAM scrip + + USE NETCDF + USE read_write_gmsh + USE WMSCRPMD + + CHARACTER(100) :: waves_mesh_file + CHARACTER(100) :: waves_scrip_file + TYPE(grid_type) :: mesh + REAL*8, DIMENSION(:,:), ALLOCATABLE :: xyb + INTEGER, DIMENSION(:,:), ALLOCATABLE :: trigp + REAL*8, DIMENSION(:), ALLOCATABLE :: grid_center_lon + REAL*8, DIMENSION(:), ALLOCATABLE :: grid_center_lat + REAL*8, DIMENSION(:,:), ALLOCATABLE :: grid_corner_lon + REAL*8, DIMENSION(:,:), ALLOCATABLE :: grid_corner_lat + LOGICAL, DIMENSION(:), ALLOCATABLE :: grid_mask + INTEGER, DIMENSION(:), ALLOCATABLE :: grid_imask + INTEGER, DIMENSION(:), ALLOCATABLE :: grid_dims + INTEGER :: grid_size, grid_corners, grid_rank + INTEGER :: NCID, IERR + INTEGER :: grid_size_dimid, grid_rank_dimid, grid_corners_dimid + INTEGER :: grid_center_lat_varid, grid_center_lon_varid + INTEGER :: grid_corner_lat_varid, grid_corner_lon_varid + INTEGER :: grid_area_varid, grid_imask_varid + INTEGER :: grid_dims_varid + CHARACTER (:), ALLOCATABLE :: GRID_UNITS, GRID_NAME + + NAMELIST / inputs / waves_mesh_file + NAMELIST / outputs / waves_scrip_file + + OPEN(UNIT=15, FILE='scrip.nml') + READ(UNIT=15, NML=inputs) + READ(UNIT=15, NML=outputs) + CLOSE(15) + + CALL read_gmsh_file(TRIM(waves_mesh_file), mesh) + + ALLOCATE(xyb(mesh%nn,3)) + DO i = 1,mesh%nn + xyb(i,1) = mesh%xy(1,i) + xyb(i,2) = mesh%xy(2,i) + xyb(i,3) = mesh%depth(i) + ENDDO + + ALLOCATE(trigp(mesh%ne,3)) + DO i = 1,mesh%ne + DO j = 1,3 + trigp(i,j) = mesh%ect(j,i) + ENDDO + ENDDO + + + CALL GET_SCRIP_INFO(mesh%ne, mesh%nn, xyb, trigp, & + grid_center_lon, grid_center_lat, & + grid_corner_lon, grid_corner_lat, grid_mask, & + grid_dims, grid_size, grid_corners, grid_rank) + + GRID_UNITS='degrees' ! the other option is radians...we don't use this + GRID_NAME='src' ! this is not used, except for netcdf output + + + IERR = NF90_CREATE(TRIM(waves_scrip_file), NF90_NETCDF4, NCID) + IERR = NF90_DEF_DIM(NCID, 'grid_size', GRID_SIZE, grid_size_dimid) + IERR = NF90_DEF_DIM(NCID, 'grid_corners', GRID_CORNERS, grid_corners_dimid) + IERR = NF90_DEF_DIM(NCID, 'grid_rank', GRID_RANK, grid_rank_dimid) + + IERR = NF90_DEF_VAR(NCID, 'grid_center_lat', NF90_DOUBLE, & + (/grid_size_dimid/),grid_center_lat_varid) + IERR = NF90_DEF_VAR(NCID, 'grid_center_lon', NF90_DOUBLE, & + (/grid_size_dimid/),grid_center_lon_varid) + IERR = NF90_DEF_VAR(NCID, 'grid_corner_lat', NF90_DOUBLE, & + (/grid_corners_dimid,grid_size_dimid/), & + grid_corner_lat_varid) + IERR = NF90_DEF_VAR(NCID, 'grid_corner_lon', NF90_DOUBLE, & + (/grid_corners_dimid,grid_size_dimid/), & + grid_corner_lon_varid) + IERR = NF90_DEF_VAR(NCID, 'grid_imask', NF90_INT, & + (/grid_size_dimid/),grid_imask_varid) + IERR = NF90_DEF_VAR(NCID, 'grid_dims', NF90_INT, & + (/grid_rank_dimid/),grid_dims_varid) + IERR = NF90_ENDDEF(NCID) + + ALLOCATE(GRID_IMASK(GRID_DIMS(1))) + GRID_IMASK = 0 + DO I = 1,GRID_DIMS(1) + IF (GRID_MASK(I)) THEN + GRID_IMASK(I) = 1 + ENDIF + ENDDO + + IERR = NF90_PUT_ATT(NCID,grid_center_lat_varid,'units',GRID_UNITS) + IERR = NF90_PUT_ATT(NCID,grid_center_lon_varid,'units',GRID_UNITS) + IERR = NF90_PUT_ATT(NCID,grid_corner_lat_varid,'units',GRID_UNITS) + IERR = NF90_PUT_ATT(NCID,grid_corner_lon_varid,'units',GRID_UNITS) + IERR = NF90_PUT_ATT(NCID,grid_imask_varid,'units','unitless') + + IERR = NF90_PUT_VAR(NCID,grid_center_lat_varid,GRID_CENTER_LAT) + IERR = NF90_PUT_VAR(NCID,grid_center_lon_varid,GRID_CENTER_LON) + IERR = NF90_PUT_VAR(NCID,grid_corner_lat_varid,GRID_CORNER_LAT) + IERR = NF90_PUT_VAR(NCID,grid_corner_lon_varid,GRID_CORNER_LON) + IERR = NF90_PUT_VAR(NCID,grid_imask_varid,GRID_IMASK) + IERR = NF90_PUT_VAR(NCID,grid_dims_varid,GRID_DIMS) + IERR = NF90_CLOSE(NCID) + +END PROGRAM scrip diff --git a/ocean/wave_mesh_tools/wmscrpmd.f90 b/ocean/wave_mesh_tools/wmscrpmd.f90 new file mode 100644 index 000000000..19b08aeba --- /dev/null +++ b/ocean/wave_mesh_tools/wmscrpmd.f90 @@ -0,0 +1,1192 @@ +!/ ------------------------------------------------------------------- / + MODULE WMSCRPMD +!/ +!/ +-----------------------------------+ +!/ | WAVEWATCH III | +!/ | E. Rogers, M. Dutour, | +!/ | A. Roland, F. Ardhuin | +!/ | FORTRAN 90 | +!/ | Last update : 10-Dec-2014 | +!/ +-----------------------------------+ +!/ +!/ 06-Sep-2012 : Origination, transfer from WMGRIDMD ( version 4.08 ) +!/ 10-Dec-2014 : Add checks for allocate status ( version 5.04 ) +!/ +!/ Copyright 2012 National Weather Service (NWS), +!/ National Oceanic and Atmospheric Administration. All rights +!/ reserved. WAVEWATCH III is a trademark of the NWS. +!/ No unauthorized use without permission. +!/ +! 1. Purpose : +! +! Routines to determine and process grid dependencies in the +! multi-grid wave model. +! +! 2. Variables and types : +! +! 3. Subroutines and functions : +! +! Name Type Scope Description +! ---------------------------------------------------------------- +! scrip_wrapper Subr. Public as the name says ... +! get_scrip_info_structured Subr. Public as the name says ... +! get_scrip_info_unstructured Subr. Public as the name says ... +! get_scrip_info Subr. Public as the name says ... +! scrip_info_renormalization Subr. Public as the name says ... +! TRIANG_INDEXES Subr. Public as the name says ... +! get_unstructured_vertex_degree Subr. Public as the name says ... +! GET_BOUNDARY Subr. Public as the name says ... +! ---------------------------------------------------------------- +! +! 4. Subroutines and functions used : +! +! Name Type Module Description +! ---------------------------------------------------------------- +! get_unstructured_vertex_degree Subr. W3TRIAMD Manage data structures +! ---------------------------------------------------------------- +! +! 5. Remarks : +! +! - The subroutines TRIANG_INDEXES, get_unstructured_vertex_degree, and +! GET_BOUNDARY should probably be renamed and moved to the module w3triamd +! +! 6. Switches : +! +! +! !/S Enable subroutine tracing. +! !/Tn Enable test output. +! +! 7. Source code : +! +!/ ------------------------------------------------------------------- / +!/ +!/ Specify default accessibility +!/ + PUBLIC + REAL*8 , ALLOCATABLE :: GRID_AREA(:) +!/ +!/ Module private variable for checking error returns +!/ + INTEGER, PRIVATE :: ISTAT +!/ + CONTAINS +!/ ------------------------------------------------------------------- / + SUBROUTINE GET_SCRIP_INFO_UNSTRUCTURED (MNE, MNP, XYB, TRIGP, & + GRID_CENTER_LON, GRID_CENTER_LAT, & + GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK, & + GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK) +!/ +-----------------------------------+ +!/ | WAVEWATCH III | +!/ | M. Dutour, A. Roland | +!/ | FORTRAN 90 | +!/ | Last update : 10-Dec-2014 ! +!/ +-----------------------------------+ +!/ +! 1. Original author : +! +! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P +! +! 2. Last update : +! +! See revisions. +! +! 3. Revisions : +! +! 20-Feb-2012 : Origination +!/ 10-Dec-2014 : Add checks for allocate status ( version 5.04 ) +! +! 4. Copyright : +! +! 5. Purpose : +! +! Compute grid arrays for scrip for a specific unstructured grid +! For interior vertices, we select for every triangle the barycenter +! of the triangle. So to every node contained in N triangles we associate +! a domain with N corners. Every one of those corners is contained +! in 3 different domains. +! For nodes that are on the boundary, we have to proceed differently. +! For every such node, we have NEIGHBOR_PREV and NEIGHBOR_NEXT which +! are the neighbor on each side of the boundary. +! We put a corner on the middle of the edge. We also put a corner +! on the node itself. +! Note that instead of taking barycenters, we could have taken +! Voronoi vertices, but this carries danger since Voronoi vertices +! can be outside of the triangle. And it leaves open how to treat +! the boundary. +! +! 6. Method : +! +! 7. Parameters, Variables and types : +! +! 8. Called by : +! +! Subroutine get_scrip_info +! +! 9. Subroutines and functions used : +! +! 10. Error messages: +! +! 11. Remarks : +! +! 12. Structure : +! +! 13. Switches : +! +! 14. Source code : + IMPLICIT NONE + INTEGER, INTENT(IN) :: MNE + INTEGER, INTENT(IN) :: MNP + REAL*8, INTENT(IN), DIMENSION(:,:) :: XYB + INTEGER, INTENT(IN), DIMENSION(:,:) :: TRIGP + REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CENTER_LON(:) + REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CENTER_LAT(:) + LOGICAL, INTENT(OUT), ALLOCATABLE :: GRID_MASK(:) + REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CORNER_LON(:,:) + REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CORNER_LAT(:,:) + INTEGER, INTENT(OUT), ALLOCATABLE :: GRID_DIMS(:) + INTEGER, INTENT(OUT) :: GRID_SIZE, GRID_CORNERS, GRID_RANK + + INTEGER DIRAPPROACH, DUALAPPROACH, THEAPPROACH + INTEGER IE, IP, I, J + INTEGER NBPLUS, NBMINUS + INTEGER I1, I2, I3 + REAL*8 :: ELON1, ELON2, ELON3, ELON, ELONC + REAL*8 :: ELAT1, ELAT2, ELAT3, ELAT, ELATC + REAL *8 :: DELTALON12, DELTALON13, DELTALAT12, DELTALAT13 + REAL *8 :: THEDET + REAL*8 :: PT(3,2) + INTEGER, POINTER :: IOBP(:), TRIGINCD(:) + INTEGER, POINTER :: NEIGHBOR_PREV(:), NEIGHBOR_NEXT(:) + INTEGER, POINTER :: NBASSIGNEDCORNER(:), LISTNBCORNER(:) + INTEGER, POINTER :: STATUS(:), NEXTVERT(:), PREVVERT(:), FINALVERT(:) + INTEGER :: MAXCORNER, NBCORNER + INTEGER :: IDX, IPNEXT, IPPREV, NB, INEXT, IPREV + REAL*8, POINTER :: LON_CENT_TRIG(:), LAT_CENT_TRIG(:) + REAL*8 :: ELONIP, ELONNEXT, ELONPREV, ELONN, ELONP + REAL*8 :: ELATIP, ELATNEXT, ELATPREV, ELATN, ELATP + INTEGER :: ISFINISHED, ZPREV + INTEGER :: DODEBUG + GRID_RANK=2 + DIRAPPROACH=1 + DUALAPPROACH=2 + THEAPPROACH=DUALAPPROACH + IF (THEAPPROACH .EQ. DIRAPPROACH) THEN + ALLOCATE(GRID_CENTER_LON(MNE), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(GRID_CENTER_LAT(MNE), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(GRID_CORNER_LON(3,MNE), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(GRID_CORNER_LAT(3,MNE), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(GRID_MASK(MNE), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + DO IE=1,MNE + I1=TRIGP(IE,1) + I2=TRIGP(IE,2) + I3=TRIGP(IE,3) + ELON1=XYB(I1,1) + ELON2=XYB(I2,1) + ELON3=XYB(I3,1) + ELAT1=XYB(I1,2) + ELAT2=XYB(I2,2) + ELAT3=XYB(I3,2) + ELON=(ELON1 + ELON2 + ELON3)/3 + ELAT=(ELAT1 + ELAT2 + ELAT3)/3 + GRID_CENTER_LON(IE)=ELON + GRID_CENTER_LAT(IE)=ELAT + GRID_CORNER_LON(1,IE)=ELON1 + GRID_CORNER_LON(2,IE)=ELON2 + GRID_CORNER_LON(3,IE)=ELON3 + GRID_CORNER_LAT(1,IE)=ELAT1 + GRID_CORNER_LAT(2,IE)=ELAT2 + GRID_CORNER_LAT(3,IE)=ELAT3 + GRID_MASK(IE)=.TRUE. + END DO + GRID_CORNERS=3 + END IF + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + IF (THEAPPROACH .EQ. DUALAPPROACH) THEN +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ALLOCATE(TRIGINCD(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(IOBP(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(NEIGHBOR_NEXT(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(NEIGHBOR_PREV(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(NBASSIGNEDCORNER(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(LISTNBCORNER(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + + ALLOCATE(STATUS(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(NEXTVERT(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(PREVVERT(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(FINALVERT(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(LON_CENT_TRIG(MNE), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(LAT_CENT_TRIG(MNE), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + + CALL GET_UNSTRUCTURED_VERTEX_DEGREE (MNP, MNE, & + TRIGP, TRIGINCD) + CALL GET_BOUNDARY(MNP, MNE, TRIGP, IOBP, & + NEIGHBOR_PREV, NEIGHBOR_NEXT) + + ! Find max number of corners + MAXCORNER=0 + DO IP=1,MNP + !IF (IP == 53591) THEN + ! PRINT*, "NEIGHBOR_NEXT", NEIGHBOR_NEXT(IP) + ! PRINT*, XYB(IP,1), XYB(IP,2), XYB(IP,3) + ! PRINT*, TRIGINCD(IP) + !ENDIF + IF (NEIGHBOR_NEXT(IP) .EQ. 0) THEN + NBCORNER=TRIGINCD(IP) + ELSE + NBCORNER=TRIGINCD(IP) + 3 + PRINT*, IP, IOBP(IP) + END IF + LISTNBCORNER(IP)=NBCORNER + IF (NBCORNER .GT. MAXCORNER) THEN + MAXCORNER=NBCORNER + END IF + END DO + GRID_CORNERS=MAXCORNER + + ALLOCATE(GRID_CENTER_LON(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(GRID_CENTER_LAT(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(GRID_CORNER_LON(MAXCORNER,MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(GRID_CORNER_LAT(MAXCORNER,MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(GRID_MASK(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + + ! Add first three corners for boundaries + NBASSIGNEDCORNER(:)=0 + DO IP=1,MNP + GRID_MASK(IP)=.TRUE. + IF (NEIGHBOR_NEXT(IP) .GT. 0) THEN + IPNEXT=NEIGHBOR_NEXT(IP) + IPPREV=NEIGHBOR_PREV(IP) + ELONIP=DBLE(XYB(IP,1)) + ELATIP=DBLE(XYB(IP,2)) + ELONNEXT=DBLE(XYB(IPNEXT,1)) + ELATNEXT=DBLE(XYB(IPNEXT,2)) + ELONPREV=DBLE(XYB(IPPREV,1)) + ELATPREV=DBLE(XYB(IPPREV,2)) + + ! Periodicity fix for corner node + IF ( ABS(ELONIP - ELONNEXT) .GT. 180.0 ) THEN + ELONNEXT = ELONNEXT -SIGN(360.0d0,(ELONIP - ELONNEXT)) + ENDIF + IF ( ABS(ELONIP - ELONPREV) .GT. 180.0 ) THEN + ELONPREV = ELONPREV -SIGN(360.0d0,(ELONIP - ELONPREV)) + ENDIF + + ELONN=(ELONIP+ELONNEXT)/2.0 + ELATN=(ELATIP+ELATNEXT)/2.0 + ELONP=(ELONIP+ELONPREV)/2.0 + ELATP=(ELATIP+ELATPREV)/2.0 + + + GRID_CORNER_LON(1,IP)=ELONN + GRID_CORNER_LAT(1,IP)=ELATN + GRID_CORNER_LON(2,IP)=ELONIP + GRID_CORNER_LAT(2,IP)=ELATIP + GRID_CORNER_LON(3,IP)=ELONP + GRID_CORNER_LAT(3,IP)=ELATP + NBASSIGNEDCORNER(IP)=3 + END IF + END DO + + ! Compute centers + DO IP=1,MNP + GRID_CENTER_LON(IP)=DBLE(XYB(IP,1)) + GRID_CENTER_LAT(IP)=DBLE(XYB(IP,2)) + END DO + + ! Check triangle node orientation + ! Compute triangle centers + NBPLUS=0 + NBMINUS=0 + DO IE=1,MNE + I1=TRIGP(IE,1) + I2=TRIGP(IE,2) + I3=TRIGP(IE,3) + + PT(1,1)=DBLE(XYB(I1,1)) + PT(2,1)=DBLE(XYB(I2,1)) + PT(3,1)=DBLE(XYB(I3,1)) + PT(1,2)=DBLE(XYB(I1,2)) + PT(2,2)=DBLE(XYB(I2,2)) + PT(3,2)=DBLE(XYB(I3,2)) + + CALL FIX_PERIODCITY(PT) + + ELON1 = PT(1,1) + ELON2 = PT(2,1) + ELON3 = PT(3,1) + ELAT1 = PT(1,2) + ELAT2 = PT(2,2) + ELAT3 = PT(3,2) + + DELTALON12=ELON2 - ELON1 + DELTALON13=ELON3 - ELON1 + DELTALAT12=ELAT2 - ELAT1 + DELTALAT13=ELAT3 - ELAT1 + THEDET=DELTALON12*DELTALAT13 - DELTALON13*DELTALAT12 + IF (THEDET.GT.0) THEN + NBPLUS=NBPLUS+1 + END IF + IF (THEDET.LT.0) THEN + NBMINUS=NBMINUS+1 + END IF + ELON=(ELON1 + ELON2 + ELON3)/3.0 + ELAT=(ELAT1 + ELAT2 + ELAT3)/3.0 + + + LON_CENT_TRIG(IE)=ELON + LAT_CENT_TRIG(IE)=ELAT + + END DO + DODEBUG=0 + IF (DODEBUG.EQ.1) THEN + print *, 'nbplus=', nbplus, ' nbminus=', nbminus + END IF + + STATUS(:) = 0 + NEXTVERT(:) = 0 + PREVVERT(:) = 0 + DO IE=1,MNE + DO I=1,3 + CALL TRIANG_INDEXES(I, INEXT, IPREV) + IP=TRIGP(IE,I) + IPNEXT=TRIGP(IE,INEXT) + IPPREV=TRIGP(IE,IPREV) + IF (STATUS(IP).EQ.0) THEN + IF (NEIGHBOR_NEXT(IP).EQ.0) THEN + STATUS(IP)=1 + FINALVERT(IP)=IPPREV + PREVVERT(IP)=IPPREV + NEXTVERT(IP)=IPNEXT + ELSE + IF (NEIGHBOR_PREV(IP).EQ.IPPREV) THEN + STATUS(IP)=1 + PREVVERT(IP)=IPPREV + NEXTVERT(IP)=IPNEXT + FINALVERT(IP)=NEIGHBOR_NEXT(IP) + END IF + END IF + END IF + END DO + END DO + STATUS(:)=0 + DO + ISFINISHED=1 + DO IE=1,MNE + ELON=LON_CENT_TRIG(IE) + ELAT=LAT_CENT_TRIG(IE) + DO I=1,3 + CALL TRIANG_INDEXES(I, INEXT, IPREV) + IP=TRIGP(IE,I) + IPNEXT=TRIGP(IE,INEXT) + IPPREV=TRIGP(IE,IPREV) + IF (STATUS(IP).EQ.0) THEN + ISFINISHED=0 + ZPREV=PREVVERT(IP) + IF (ZPREV.EQ.IPPREV) THEN + IDX=NBASSIGNEDCORNER(IP) + IDX=IDX+1 + GRID_CORNER_LON(IDX,IP)=ELON + GRID_CORNER_LAT(IDX,IP)=ELAT + NBASSIGNEDCORNER(IP)=IDX + PREVVERT(IP)=IPNEXT + IF (IPNEXT.EQ.FINALVERT(IP)) THEN + STATUS(IP)=1 + END IF + END IF + END IF + END DO + END DO + IF (ISFINISHED.EQ.1) THEN + EXIT + END IF + END DO + DO IP=1,MNP + IF (NBASSIGNEDCORNER(IP).NE.LISTNBCORNER(IP)) THEN + WRITE(*,*) 'Incoherent number at IP=', IP + WRITE(*,*) ' NbAssignedCorner(IP)=', NbAssignedCorner(IP) + WRITE(*,*) ' ListNbCorner(IP)=', ListNbCorner(IP) + WRITE(*,*) ' N_N=', NEIGHBOR_NEXT(IP), 'N_P=', NEIGHBOR_PREV(IP) + WRITE(*,*) ' TrigIncd=', TrigIncd(IP) + !STOP 'wmscrpmd, case 2' + END IF + END DO + + ! if the number of corner is below threshold, we have to + ! add some more. + DO IP=1,MNP + NB=NBASSIGNEDCORNER(IP) + IF (NB .LT. MAXCORNER) THEN + ELON=GRID_CORNER_LON(NB,IP) + ELAT=GRID_CORNER_LAT(NB,IP) + DO IDX=NB+1,MAXCORNER + GRID_CORNER_LON(IDX,IP)=ELON + GRID_CORNER_LAT(IDX,IP)=ELAT + END DO + END IF + END DO + + ! Calculate area of cell + ALLOCATE(GRID_AREA(MNP)) + DO I1=1,MNP + GRID_AREA(I1) = 0.0 + NB=NBASSIGNEDCORNER(I1) + DO I2=1,NB + I3 = MOD(I2,NB)+1 + PT(1,1)=GRID_CENTER_LON(I1) + PT(2,1)=GRID_CORNER_LON(I2,I1) + PT(3,1)=GRID_CORNER_LON(I3,I1) + PT(1,2)=GRID_CENTER_LAT(I1) + PT(2,2)=GRID_CORNER_LAT(I2,I1) + PT(3,2)=GRID_CORNER_LON(I3,I1) + + CALL FIX_PERIODCITY(PT) + + ELON1 = PT(1,1) + ELON2 = PT(2,1) + ELON3 = PT(3,1) + ELAT1 = PT(1,2) + ELAT2 = PT(2,2) + ELAT3 = PT(3,2) + + DELTALON12=ELON2 - ELON1 + DELTALON13=ELON3 - ELON1 + DELTALAT12=ELAT2 - ELAT1 + DELTALAT13=ELAT3 - ELAT1 + THEDET=DELTALON12*DELTALAT13 - DELTALON13*DELTALAT12 + + GRID_AREA(I1) = GRID_AREA(I1) + ABS(THEDET) + ENDDO + ENDDO + + DEALLOCATE(NBASSIGNEDCORNER, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(LISTNBCORNER, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(TRIGINCD, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(IOBP, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(NEIGHBOR_PREV, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(NEIGHBOR_NEXT, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(STATUS, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(NEXTVERT, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(PREVVERT, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(FINALVERT, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(LON_CENT_TRIG, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(LAT_CENT_TRIG, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + + GRID_SIZE=MNP + GRID_RANK=1 + ALLOCATE(GRID_DIMS(GRID_RANK)) + GRID_DIMS(1) = GRID_SIZE + + DO I = 1,GRID_SIZE + IF (GRID_CENTER_LON(I) < 0.0) THEN + GRID_CENTER_LON(I) = GRID_CENTER_LON(I)+360.0 + ENDIF + DO J = 1,GRID_CORNERS + IF (GRID_CORNER_LON(J,I) < 0.0) THEN + GRID_CORNER_LON(J,I) = GRID_CORNER_LON(J,I)+360.0 + ENDIF + ENDDO + ENDDO + + END IF + END SUBROUTINE GET_SCRIP_INFO_UNSTRUCTURED + +!/ ------------------------------------------------------------------- / + SUBROUTINE GET_SCRIP_INFO(MNE, MNP, XYB, TRIGP, & + & GRID_CENTER_LON, GRID_CENTER_LAT, & + & GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK, & + & GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK) +! 1. Original author : +! +! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P +! +! 2. Last update : +! +! See revisions. +! +! 3. Revisions : +! +! 20-Feb-2012 : Origination, this is adapted from Erick Rogers +! code by splitting the code into sections. +! +! 4. Copyright : +! +! 5. Purpose : +! +! Compute grid for scrip for a specific structured grid. +! This is adapted from Erick Rogers code by making it cleaner. +! +! 6. Method : +! +! 7. Parameters, Variables and types : +! +! 8. Called by : +! +! Subroutine scrip_wrapper +! +! 9. Subroutines and functions used : +! +! 10. Error messages: +! +! 11. Remarks : +! +! 12. Structure : +! +! 13. Switches : +! +! 14. Source code : + IMPLICIT NONE + INTEGER, INTENT(IN) :: MNE + INTEGER, INTENT(IN) :: MNP + REAL*8, INTENT(IN), DIMENSION(:,:) :: XYB + INTEGER, INTENT(IN), DIMENSION(:,:) :: TRIGP + REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CENTER_LON(:) + REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CENTER_LAT(:) + LOGICAL, INTENT(OUT), ALLOCATABLE :: GRID_MASK(:) + REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CORNER_LON(:,:) + REAL*8, INTENT(OUT), ALLOCATABLE :: GRID_CORNER_LAT(:,:) + INTEGER, INTENT(OUT), ALLOCATABLE :: GRID_DIMS(:) + INTEGER, INTENT(OUT) :: GRID_SIZE, GRID_CORNERS, GRID_RANK + REAL*8 :: DLON1, DLAT1, DLON2, DLAT2, THEDET + INTEGER :: I, J + INTEGER :: IC, JC, IP, CHECKSIGNS, NBPLUS, NBMINUS, NBZERO + INTEGER :: PRINTDATA, PRINTMINMAX + REAL*8 :: MINLON, MINLAT, MAXLON, MAXLAT + REAL*8 :: MINLONCORNER, MAXLONCORNER, MINLATCORNER, MAXLATCORNER + REAL*8 :: PT(3,2) + CALL GET_SCRIP_INFO_UNSTRUCTURED (MNE, MNP, XYB, TRIGP, & + & GRID_CENTER_LON, GRID_CENTER_LAT, & + & GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK, & + & GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK) + CHECKSIGNS=1 + IF (CHECKSIGNS.EQ.1) THEN + NBPLUS=0 + NBMINUS=0 + NBZERO=0 + DO IP=1,GRID_SIZE + DO IC=1,GRID_CORNERS + IF (IC.EQ.GRID_CORNERS) THEN + JC=1 + ELSE + JC=IC+1 + END IF + + PT(1,1) = GRID_CENTER_LON(IP) + PT(1,2) = GRID_CENTER_LAT(IP) + PT(2,1) = GRID_CORNER_LON(IC,IP) + PT(2,2) = GRID_CORNER_LAT(IC,IP) + PT(3,1) = GRID_CORNER_LON(JC,IP) + PT(3,2) = GRID_CORNER_LAT(JC,IP) + + CALL FIX_PERIODCITY(PT) + + DLON1=PT(2,1)-PT(1,1) + DLON2=PT(3,1)-PT(1,1) + DLAT1=PT(2,2)-PT(1,2) + DLAT2=PT(3,2)-PT(1,2) + + THEDET=DLON1*DLAT2 - DLON2*DLAT1 + IF (THEDET.GT.1d-8) THEN + NBPLUS=NBPLUS+1 + ELSE IF (THEDET.LT.-1d-8) THEN + NBMINUS=NBMINUS+1 + ELSE + NBZERO=NBZERO+1 + END IF + END DO + END DO + + WRITE(*,*) 'SI nbPlus=', nbPlus, ' nbMinus=', nbMinus, ' nbZero=', nbZero + + END IF + END SUBROUTINE GET_SCRIP_INFO + +!/ ------------------------------------------------------------------- / + SUBROUTINE SCRIP_INFO_RENORMALIZATION( & + & GRID_CENTER_LON, GRID_CENTER_LAT, & + & GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK, & + & GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK, & + & CONV_DX, CONV_DY, OFFSET, GRIDSHIFT) +! 1. Original author : +! +! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P +! Adapted from Erick Rogers scrip_wrapper +! +! 2. Last update : +! +! See revisions. +! +! 3. Revisions : +! +! 20-Feb-2012 : Origination +! +! 4. Copyright : +! +! 5. Purpose : +! +! This is adapted from Erick Rogers scrip_wrapper +! Purpose is to rescale according to whether the grid is spherical +! or not and to adjust by some small shift to keep SCRIP happy +! in situations where nodes of different grids overlay +! +! 6. Method : +! +! We apply various transformations to the longitude latitude +! following here the transformations that were done only in +! finite difference meshes. +! +! 7. Parameters, Variables and types : +! +! 8. Called by : +! +! Subroutine WMGHGH +! +! 9. Subroutines and functions used : +! +! 10. Error messages: +! +! 11. Remarks : +! +! 12. Structure : +! +! 13. Switches : +! +! 14. Source code : + IMPLICIT NONE + REAL*8, INTENT(INOUT) :: GRID_CENTER_LON(:) + REAL*8, INTENT(INOUT) :: GRID_CENTER_LAT(:) + LOGICAL, INTENT(IN) :: GRID_MASK(:) + REAL*8, INTENT(INOUT) :: GRID_CORNER_LON(:,:) + REAL*8, INTENT(INOUT) :: GRID_CORNER_LAT(:,:) + INTEGER, INTENT(IN) :: GRID_DIMS(:) + INTEGER, INTENT(IN) :: GRID_SIZE, GRID_CORNERS, GRID_RANK + REAL*8 :: CONV_DX, CONV_DY, OFFSET, GRIDSHIFT + REAL*8 DEG2RAD + ! + INTEGER :: I, J, IP + REAL*8 :: MINLON, MINLAT, MAXLON, MAXLAT, HLON, HLAT + REAL*8 :: MINLONCORNER, MAXLONCORNER, MINLATCORNER, MAXLATCORNER + + DO I=1,GRID_SIZE + GRID_CENTER_LON(I)=(GRID_CENTER_LON(I)+OFFSET)/CONV_DX + & + & GRIDSHIFT + GRID_CENTER_LAT(I)=GRID_CENTER_LAT(I)/CONV_DY + & + & GRIDSHIFT + IF(GRID_CENTER_LON(I)>360.0) THEN + GRID_CENTER_LON(I)=GRID_CENTER_LON(I)-360.0 + END IF + IF(GRID_CENTER_LON(I)<000.0) THEN + GRID_CENTER_LON(I)=GRID_CENTER_LON(I)+360.0 + END IF + DO J=1,GRID_CORNERS + GRID_CORNER_LON(J, I)=(GRID_CORNER_LON(J, I)+OFFSET)/CONV_DX+ & + & GRIDSHIFT + GRID_CORNER_LAT(J, I)=GRID_CORNER_LAT(J, I)/CONV_DY + & + & GRIDSHIFT + IF(GRID_CORNER_LON(J,I)>360.0) THEN + GRID_CORNER_LON(J,I)=GRID_CORNER_LON(J,I)-360.0 + END IF + IF(GRID_CORNER_LON(J,I)<000.0) THEN + GRID_CORNER_LON(J,I)=GRID_CORNER_LON(J,I)+360.0 + END IF + END DO + END DO + + END SUBROUTINE SCRIP_INFO_RENORMALIZATION + +!/ ------------------------------------------------------------------- / + SUBROUTINE TRIANG_INDEXES(I, INEXT, IPREV) +! 1. Original author : +! +! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P +! + INTEGER, INTENT(IN) :: I + INTEGER, INTENT(OUT) :: INEXT, IPREV + IF (I.EQ.1) THEN + INEXT=3 + ELSE + INEXT=I-1 + END IF + IF (I.EQ.3) THEN + IPREV=1 + ELSE + IPREV=I+1 + END IF + END SUBROUTINE + +!/ ------------------------------------------------------------------- / + SUBROUTINE GET_UNSTRUCTURED_VERTEX_DEGREE (MNP, MNE, TRIGP, & + & TRIGINCD) +! Written: +! +! 20-Feb-2012 +! +! Author: +! +! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P +! +! Parameters: +! Input: +! MNP: number of nodes +! INE: list of nodes +! MNE: number of triangles +! Output: +! TrigIncd (number of triangles contained by vertices +! +! Description: +! this function returns the list of incidences +! + IMPLICIT NONE + INTEGER, INTENT(IN) :: MNP, MNE + INTEGER, INTENT(IN) :: TRIGP(:,:) + INTEGER, INTENT(OUT) :: TRIGINCD(:) + INTEGER :: IP, IE, I + TRIGINCD=0 + DO IE=1,MNE + DO I=1,3 + IP=TRIGP(IE,I) + TRIGINCD(IP)=TRIGINCD(IP) + 1 + END DO + END DO + END SUBROUTINE GET_UNSTRUCTURED_VERTEX_DEGREE + +!/ ------------------------------------------------------------------- / + SUBROUTINE GET_BOUNDARY(MNP, MNE, TRIGP, IOBP, NEIGHBOR_PREV, & + & NEIGHBOR_NEXT) +!/ +-----------------------------------+ +!/ | WAVEWATCH III | +!/ | M. Dutour, A. Roland | +!/ | FORTRAN 90 | +!/ | Last update : 10-Dec-2014 ! +!/ +-----------------------------------+ +!/ +! Written: +! +! 20-Feb-2012 +!/ 10-Dec-2014 : Add checks for allocate status ( version 5.04 ) +! +! Author: +! +! Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P +! +! Parameters: +! Input: +! MNP: number of nodes +! TRIGP: list of nodes +! MNE: number of triangles +! Output: +! NEIGHBOR +! +! Description: +! if a node belong to a boundary, the function +! returns the neighbor of this point on one side. +! if the point is interior then the value 0 is set. +! + IMPLICIT NONE + + INTEGER, INTENT(IN) :: MNP, MNE, TRIGP(MNE,3) + INTEGER, INTENT(INOUT) :: IOBP(MNP) + INTEGER, INTENT(INOUT) :: NEIGHBOR_PREV(MNP) + INTEGER, INTENT(INOUT) :: NEIGHBOR_NEXT(MNP) + + INTEGER, POINTER :: STATUS(:) + INTEGER, POINTER :: COLLECTED(:) + INTEGER, POINTER :: NEXTVERT(:) + INTEGER, POINTER :: PREVVERT(:) + + INTEGER :: IE, I, IP, IP2, IP3 + INTEGER :: ISFINISHED, INEXT, IPREV + INTEGER :: IPNEXT, IPPREV, ZNEXT, ZPREV + + ALLOCATE(STATUS(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(COLLECTED(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(PREVVERT(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + ALLOCATE(NEXTVERT(MNP), STAT=ISTAT) + CALL CHECK_ALLOC_STATUS ( ISTAT ) + IOBP = 0 + NEIGHBOR_NEXT = 0 + NEIGHBOR_PREV = 0 + +! Now computing the next items + STATUS = 0 + NEXTVERT = 0 + PREVVERT = 0 + DO IE=1,MNE + DO I=1,3 + CALL TRIANG_INDEXES(I, INEXT, IPREV) + IP=TRIGP(IE,I) + IPNEXT=TRIGP(IE,INEXT) + IPPREV=TRIGP(IE,IPREV) + IF (STATUS(IP).EQ.0) THEN + STATUS(IP)=1 + PREVVERT(IP)=IPPREV + NEXTVERT(IP)=IPNEXT + END IF + END DO + END DO + STATUS(:)=0 + DO + COLLECTED(:)=0 + DO IE=1,MNE + DO I=1,3 + CALL TRIANG_INDEXES(I, INEXT, IPREV) + IP=TRIGP(IE,I) + IPNEXT=TRIGP(IE,INEXT) + IPPREV=TRIGP(IE,IPREV) + IF (STATUS(IP).EQ.0) THEN + ZNEXT=NEXTVERT(IP) + IF (ZNEXT.EQ.IPPREV) THEN + COLLECTED(IP)=1 + NEXTVERT(IP)=IPNEXT + IF (NEXTVERT(IP).EQ.PREVVERT(IP)) THEN + STATUS(IP)=1 + END IF + END IF + END IF + END DO + END DO + + ISFINISHED=1 + DO IP=1,MNP + IF ((COLLECTED(IP).EQ.0).AND.(STATUS(IP).EQ.0)) THEN + STATUS(IP)=-1 + NEIGHBOR_NEXT(IP)=NEXTVERT(IP) + END IF + IF (STATUS(IP).EQ.0) THEN + ISFINISHED=0 + END IF + END DO + IF (ISFINISHED.EQ.1) THEN + EXIT + END IF + END DO + +! Now computing the prev items + STATUS = 0 + NEXTVERT = 0 + PREVVERT = 0 + DO IE=1,MNE + DO I=1,3 + CALL TRIANG_INDEXES(I, INEXT, IPREV) + IP=TRIGP(IE,I) + IPNEXT=TRIGP(IE,INEXT) + IPPREV=TRIGP(IE,IPREV) + IF (STATUS(IP).EQ.0) THEN + STATUS(IP)=1 + PREVVERT(IP)=IPPREV + NEXTVERT(IP)=IPNEXT + END IF + END DO + END DO + STATUS(:)=0 + DO + COLLECTED(:)=0 + DO IE=1,MNE + DO I=1,3 + CALL TRIANG_INDEXES(I, INEXT, IPREV) + IP=TRIGP(IE,I) + IPNEXT=TRIGP(IE,INEXT) + IPPREV=TRIGP(IE,IPREV) + IF (STATUS(IP).EQ.0) THEN + ZPREV=PREVVERT(IP) + IF (ZPREV.EQ.IPNEXT) THEN + COLLECTED(IP)=1 + PREVVERT(IP)=IPPREV + IF (PREVVERT(IP).EQ.NEXTVERT(IP)) THEN + STATUS(IP)=1 + END IF + END IF + END IF + END DO + END DO + + ISFINISHED=1 + DO IP=1,MNP + IF ((COLLECTED(IP).EQ.0).AND.(STATUS(IP).EQ.0)) THEN + STATUS(IP)=-1 + NEIGHBOR_PREV(IP)=PREVVERT(IP) ! new code + END IF + IF (STATUS(IP).EQ.0) THEN + ISFINISHED=0 + END IF + END DO + IF (ISFINISHED.EQ.1) THEN + EXIT + END IF + END DO +! Now making checks + DO IP=1,MNP + IP2=NEIGHBOR_NEXT(IP) + IF (IP2.GT.0) THEN + IP3=NEIGHBOR_PREV(IP2) + IF (ABS(IP3 - IP).GT.0) THEN + WRITE(*,*) 'IP=', IP, ' IP2=', IP2, ' IP3=', IP3 + WRITE(*,*) 'We have a dramatic inconsistency' + STOP 'wmscrpmd, case 3' + END IF + END IF + END DO +! Now assigning the boundary IOBP array + DO IP=1,MNP + IF (STATUS(IP).EQ.-1 .AND. IOBP(IP) .EQ. 0) THEN + IOBP(IP)=1 + END IF + END DO + + DEALLOCATE(STATUS, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(COLLECTED, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(NEXTVERT, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + DEALLOCATE(PREVVERT, STAT=ISTAT) + CALL CHECK_DEALLOC_STATUS ( ISTAT ) + + END SUBROUTINE +!/ ------------------------------------------------------------------- / + SUBROUTINE FIX_PERIODCITY(PT) + +!/ +!/ +-----------------------------------+ +!/ | WAVEWATCH III NOAA/NCEP | +!/ | Steven Brus | +!/ | Ali Abdolali | +!/ | FORTRAN 90 | +!/ | Last update : 21-May-2020 | +!/ +-----------------------------------+ +!/ +!/ 21-May-2020 : Origination. ( version 6.07 ) +!/ +!/ +! 1. Purpose : +! +! Adjust element longitude coordinates for elements straddling the +! dateline with distance of ~360 degrees +! +! 2. Method : +! +! Detect if element has nodes on both sides of dateline and adjust +! coordinates so that all nodes have the same sign +! +! 3. Parameters : +! +! Parameter list +! ---------------------------------------------------------------- + IMPLICIT NONE + REAL*8, INTENT(INOUT) :: PT(3,2) +! ---------------------------------------------------------------- +! +! Local variables. +! ---------------------------------------------------------------- + INTEGER :: I + INTEGER :: R1GT180, R2GT180, R3GT180 +! ---------------------------------------------------------------- +! +! 4. Subroutines used : +! + +! 5. Called by : +! +! Name Type Module Description +! ---------------------------------------------------------------- +! GET_SCRIP_INFO_UNSTRUCTURED Subr. WMSCRPMD Element center calculation +! GET_SCRIP_INFO Subr. WMSCRPMD Check signs +! ---------------------------------------------------------------- +! +! 6. Error messages : +! +! None. +! +! 7. Remarks : +! +! 8. Structure : +! +! 9. Switches : +! +! 10. Source code : +!/ ------------------------------------------------------------------- / + + R1GT180 = MERGE(1, 0, ABS(PT(3,1)-PT(2,1)).GT.180) + R2GT180 = MERGE(1, 0, ABS(PT(1,1)-PT(3,1)).GT.180) + R3GT180 = MERGE(1, 0, ABS(PT(2,1)-PT(1,1)).GT.180) +! if R1GT180+R2GT180+R3GT180 .eq. 0 the element does not cross the +! dateline +! if R1GT180+R2GT180+R3GT180 .eq. 1 the element contains the pole +! if R1GT180+R2GT180+R3GT180 .eq. 2 the element crosses the dateline + + IF ( R1GT180 + R2GT180 == 2 ) THEN + PT(3,1)=PT(3,1)-SIGN(360.0d0,(PT(3,1)-PT(2,1))) + ELSE IF ( R2GT180 + R3GT180 == 2 ) THEN + PT(1,1)=PT(1,1)-SIGN(360.0d0,(PT(1,1)-PT(2,1))) + ELSE IF ( R1GT180 + R3GT180 == 2 ) THEN + PT(2,1)=PT(2,1)-SIGN(360.0d0,(PT(2,1)-PT(3,1))) + ENDIF + + RETURN + END SUBROUTINE FIX_PERIODCITY + +!/ ------------------------------------------------------------------- / + SUBROUTINE EXTCDE ( IEXIT, UNIT, MSG, FILE, LINE, COMM ) +!/ +!/ +-----------------------------------+ +!/ | WAVEWATCH III NOAA/NCEP | +!/ | H. L. Tolman | +!/ | FORTRAN 90 | +!/ | Last update : 06-Jun-2018 | +!/ +-----------------------------------+ +!/ +!/ 06-Jan-1998 : Final FORTRAN 77 ( version 1.18 ) +!/ 23-Nov-1999 : Upgrade to FORTRAN 90 ( version 2.00 ) +!/ 10-Dec-2014 : Add checks for allocate status ( version 5.04 ) +!/ 11-Mar-2015 : Allow non-error exit (iexit=0) ( version 5.04 ) +!/ 20-Jan-2017 : Add optional MPI communicator arg ( version 6.02 ) +!/ 06-Jun-2018 : Add optional MPI ( version 6.04 ) +!/ +! 1. Purpose : +! +! Perform a program stop with an exit code. +! +! If exit code IEXIT=0, then it is not an error, but +! a stop has been requested by the calling routine: +! wait for other processes in communicator to catch up. +! +! If exit code IEXIT.ne.0, then abort program w/out +! waiting for other processes to catch up (important for example +! when not all processes are used by WW3). +! +! 2. Method : +! +! Machine dependent. +! +! 3. Parameters : +! +! Parameter list +! ---------------------------------------------------------------- +! IEXIT Int. I Exit code to be used. +! UNIT Int. I (optional) file unit to write error message +! MSG Str. I (optional) error message +! FILE Str. I (optional) name of source code file +! LINE Int. I (optional) line number in source code file +! COMM Int. I (optional) MPI communicator +! ---------------------------------------------------------------- +! +! 4. Subroutines used : +! +! 5. Called by : +! +! Any. +! +! 9. Switches : +! +! !/MPI MPI finalize interface if active +! +! 10. Source code : +! +!/ ------------------------------------------------------------------- / + IMPLICIT NONE +! +!/MPI INCLUDE "mpif.h" +!/ +!/ ------------------------------------------------------------------- / +!/ Parameter list +!/ + INTEGER, INTENT(IN) :: IEXIT + INTEGER, INTENT(IN), OPTIONAL :: UNIT + CHARACTER(*), INTENT(IN), OPTIONAL :: MSG + CHARACTER(*), INTENT(IN), OPTIONAL :: FILE + INTEGER, INTENT(IN), OPTIONAL :: LINE + INTEGER, INTENT(IN), OPTIONAL :: COMM +!/ +!/ ------------------------------------------------------------------- / +!/ +!/MPI INTEGER :: IERR_MPI +!/MPI LOGICAL :: RUN + INTEGER :: IUN + CHARACTER(256) :: LMSG = "" + CHARACTER(6) :: LSTR + CHARACTER(10) :: PREFIX = "WW3 ERROR:" +!/ +!/ Set file unit for error output +!/ + IUN = 0 + IF (PRESENT(UNIT)) IUN = UNIT +!/ +!/ Report error message +!/ + IF (PRESENT(MSG)) THEN + WRITE (IUN,"(A)") PREFIX//" "//TRIM(MSG) + END IF +!/ +!/ Report context +!/ + IF ( PRESENT(FILE) ) THEN + LMSG = TRIM(LMSG)//" FILE="//TRIM(FILE) + END IF + IF ( PRESENT(LINE) ) THEN + WRITE (LSTR,'(I0)') LINE + LMSG = TRIM(LMSG)//" LINE="//TRIM(LSTR) + END IF + IF ( LEN_TRIM(LMSG).GT.0 ) THEN + WRITE (IUN,"(A)") PREFIX//TRIM(LMSG) + END IF +!/F90 CALL EXIT ( IEXIT ) +!/DUM STOP +!/ +!/ End of EXTCDE ----------------------------------------------------- / +!/ + END SUBROUTINE EXTCDE + + SUBROUTINE CHECK_ALLOC_STATUS( STAT ) + INTEGER, INTENT(IN) :: STAT + IF ( STAT .NE. 0 ) & + CALL EXTCDE ( 99, MSG="ALLOCATE FAILED") + END SUBROUTINE CHECK_ALLOC_STATUS + + SUBROUTINE CHECK_DEALLOC_STATUS( STAT ) + INTEGER, INTENT(IN) :: STAT + IF ( STAT .NE. 0 ) & + CALL EXTCDE ( 99, MSG="DEALLOCATE FAILED") + END SUBROUTINE CHECK_DEALLOC_STATUS + +!/ +!/ End of module WMSCRPMD -------------------------------------------- / +!/ + END MODULE WMSCRPMD diff --git a/ocean/wave_mesh_tools/write_vtk.f90 b/ocean/wave_mesh_tools/write_vtk.f90 new file mode 100644 index 000000000..8df9c63d9 --- /dev/null +++ b/ocean/wave_mesh_tools/write_vtk.f90 @@ -0,0 +1,205 @@ + MODULE write_vtk + + USE globals, ONLY: rp + + IMPLICIT NONE + + CONTAINS + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE write_vtk_file(filename,nn,xy,ne,ect,depth,nope,obseg,obnds,nbou,fbseg,fbnds,nsnap,solution) + + IMPLICIT NONE + + CHARACTER(*), INTENT(IN) :: filename + INTEGER, INTENT(IN) :: nn + REAL(rp), DIMENSION(:,:), INTENT(IN) :: xy + INTEGER, INTENT(IN) :: ne + INTEGER, DIMENSION(:,:), INTENT(IN) :: ect + REAL(rp), DIMENSION(:), INTENT(IN), OPTIONAL :: depth + INTEGER, INTENT(IN),OPTIONAL :: nope + INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: obseg + INTEGER, DIMENSION(:,:), INTENT(IN), OPTIONAL :: obnds + INTEGER, INTENT(IN), OPTIONAL :: nbou + INTEGER, INTENT(IN), DIMENSION(:,:), OPTIONAL :: fbseg + INTEGER, INTENT(IN), DIMENSION(:,:), OPTIONAL :: fbnds + INTEGER, INTENT(IN), OPTIONAL :: nsnap + REAL(rp), INTENT(IN), DIMENSION(:,:), OPTIONAL :: solution + + INTEGER :: open_bou_flag + INTEGER :: flow_bou_flag + INTEGER :: depth_flag + INTEGER :: solution_flag + + CHARACTER(100) :: tmp1,tmp2,tmp3 + CHARACTER(3) :: tmp4 + INTEGER :: i,j + INTEGER :: tbnds + INTEGER :: snap + INTEGER :: R3 + + open_bou_flag = 0 + IF (PRESENT(nope) .and. PRESENT(obseg) .and. PRESENT(obnds)) THEN + open_bou_flag = 1 + ENDIF + + flow_bou_flag = 0 + IF (PRESENT(nbou) .and. PRESENT(fbseg) .and. PRESENT(fbnds)) THEN + flow_bou_flag = 1 + ENDIF + + depth_flag = 0 + IF (PRESENT(depth)) THEN + depth_flag = 1 + ENDIF + + solution_flag = 0 + IF (PRESENT(solution)) THEN + solution_flag = 1 + ENDIF + + OPEN(UNIT=11,FILE=TRIM(filename)) + + ! Write header + WRITE(11,"(A)") "# vtk DataFile Version 3.0" + WRITE(11,"(A)") filename + WRITE(11,"(A)") "ASCII" + WRITE(11,"(A)") "DATASET UNSTRUCTURED_GRID" + WRITE(tmp1,"(I10)") nn + WRITE(11,"(A,1x,A,1x,A)") "POINTS", TRIM(ADJUSTL(tmp1)), "double" + + IF (SIZE(xy,1) == 3) THEN + R3 = 1 + ELSE IF (SIZE(xy,1) == 2) THEN + R3 = 0 + ELSE + PRINT*, "Incorrect xy size" + STOP + ENDIF + + ! Write point coordinates + WRITE(tmp3,"(I4)") 0 + DO i = 1,nn + WRITE(tmp1,"(F20.14)") xy(1,i) + WRITE(tmp2,"(F20.14)") xy(2,i) + IF (R3 == 1) THEN + WRITE(tmp3,"(F20.14)") xy(3,i) + ENDIF + WRITE(11,"(A,1x,A,1x,A)") TRIM(ADJUSTL(tmp1)),TRIM(ADJUSTL(tmp2)),TRIM(ADJUSTL(tmp3)) + ENDDO + + + ! Coutn total boundary nodes + tbnds = 0 + IF (open_bou_flag == 1) THEN + DO i = 1,nope + tbnds = tbnds + obseg(i)-1 + ENDDO + ENDIF + IF (flow_bou_flag == 1) THEN + DO i = 1,nbou + tbnds = tbnds + fbseg(1,i)-1 + ENDDO + ENDIF + + ! Write boundary and element connectivity + WRITE(tmp1,"(I10)") ne+tbnds + WRITE(tmp2,"(I10)") 4*(ne+tbnds) + WRITE(11,"(A,1x,A,1x,A)") "CELLS",TRIM(ADJUSTL(tmp1)), TRIM(ADJUSTL(tmp2)) + IF (open_bou_flag == 1) THEN + DO i = 1,nope + DO j = 1,obseg(i)-1 + WRITE(tmp1,"(I10)") obnds(j,i)-1 + WRITE(tmp2,"(I10)") obnds(j+1,i)-1 + WRITE(11,"(I1,1x,A,1x,A)") 2,TRIM(ADJUSTL(tmp1)),TRIM(ADJUSTL(tmp2)) + ENDDO + ENDDO + ENDIF + IF (flow_bou_flag == 1) THEN + DO i = 1,nbou + DO j = 1,fbseg(1,i)-1 + WRITE(tmp1,"(I10)") fbnds(j,i)-1 + WRITE(tmp2,"(I10)") fbnds(j+1,i)-1 + WRITE(11,"(I1,1x,A,1x,A)") 2,TRIM(ADJUSTL(tmp1)),TRIM(ADJUSTL(tmp2)) + ENDDO + ENDDO + ENDIF + DO i = 1,ne + WRITE(tmp1,"(I10)") ect(1,i)-1 + WRITE(tmp2,"(I10)") ect(2,i)-1 + WRITE(tmp3,"(I10)") ect(3,i)-1 + WRITE(11,"(I1,1x,A,1x,A,1x,A)") 3,TRIM(ADJUSTL(tmp1)), TRIM(ADJUSTL(tmp2)), TRIM(ADJUSTL(tmp3)) + ENDDO + + ! Write cell types + WRITE(tmp1,"(I10)") ne+tbnds + WRITE(11,"(A,1x,A)") "CELL_TYPES", TRIM(ADJUSTL(tmp1)) + IF (open_bou_flag == 1) THEN + DO i = 1,nope + DO j = 1,obseg(i)-1 + WRITE(11,"(I1)") 3 + ENDDO + ENDDO + ENDIF + IF (flow_bou_flag == 1) THEN + DO i = 1,nbou + DO j = 1,fbseg(1,i)-1 + WRITE(11,"(I1)") 3 + ENDDO + ENDDO + ENDIF + DO i = 1,ne + WRITE(11,"(I1)") 5 + ENDDO + + ! Write depth + IF (depth_flag == 1) THEN + WRITE(tmp1,"(I10)") nn + WRITE(11,"(A,1x,A)") "POINT_DATA", TRIM(ADJUSTL(tmp1)) + WRITE(11,"(A)") "SCALARS bathymetry float 1" + WRITE(11,"(A)") "LOOKUP_TABLE default" + DO i = 1,nn + WRITE(tmp1,"(F24.7)") depth(i) + WRITE(11,"(A)") TRIM(ADJUSTL(tmp1)) + ENDDO + ENDIF + + IF (solution_flag == 1) THEN + WRITE(tmp1,"(I10)") nn + WRITE(11,"(A,1x,A)") "POINT_DATA", TRIM(ADJUSTL(tmp1)) + DO snap = 1,nsnap + WRITE(tmp1,"(I10)") snap + tmp4 = '000' + WRITE(tmp4,"(I3.3)") snap + WRITE(11,"(A)") "SCALARS t="//TRIM(ADJUSTL(tmp4))//" float 1" + WRITE(11,"(A)") "LOOKUP_TABLE default" + DO i = 1,nn + WRITE(tmp1,"(F24.7)") solution(i,snap) + WRITE(11,"(A)") TRIM(ADJUSTL(tmp1)) + ENDDO + ENDDO + !WRITE(tmp1,"(I10)") nsnap + !WRITE(11,"(A,1x,A)") "FIELD solution", TRIM(ADJUSTL(tmp1)) + !DO snap = 1,nsnap + ! WRITE(tmp1,"(I10)") snap + ! WRITE(tmp2,"(I10)") nn + ! WRITE(11,"(A)") "t="//TRIM(ADJUSTL(tmp1))//" 1 "//TRIM(ADJUSTL(tmp2))//" float" + ! DO i = 1,nn + ! WRITE(tmp1,"(F24.7)") solution(i,snap) + ! WRITE(11,"(A)") TRIM(ADJUSTL(tmp1)) + ! ENDDO + !ENDDO + ENDIF + + CLOSE(11) + + RETURN + END SUBROUTINE write_vtk_file + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + END MODULE write_vtk +