module remapper use mpas_mesh use scan_input, only : input_field_type, FIELD_TYPE_REAL, FIELD_TYPE_DOUBLE, FIELD_TYPE_INTEGER use target_mesh integer, parameter :: max_queue_length = 2700 ! should be at least (earth circumference / minimum grid distance) integer, parameter :: max_dictionary_size = 82000 ! should be at least (nCells/32) integer :: queue_head = 0 integer :: queue_tail = 0 integer :: queue_size = 0 integer, dimension(0:max_queue_length-1) :: queue_array integer :: int_size = 32 integer, dimension(max_dictionary_size) :: dictionary_array type remap_info_type integer :: method = -1 type (mpas_mesh_type), pointer :: src_mesh type (target_mesh_type), pointer :: dst_mesh ! For nearest-neighbor integer, dimension(:,:), pointer :: nearestCell => null() integer, dimension(:,:), pointer :: nearestVertex => null() integer, dimension(:,:), pointer :: nearestEdge => null() ! For Wachspress interpolation real, dimension(:,:,:), pointer :: cellWeights => null() real, dimension(:,:,:), pointer :: vertexWeights => null() real, dimension(:,:,:), pointer :: edgeWeights => null() integer, dimension(:,:,:), pointer :: sourceCells => null() integer, dimension(:,:,:), pointer :: sourceVertices => null() integer, dimension(:,:,:), pointer :: sourceEdges => null() ! For masked interpolation real, dimension(:,:,:), pointer :: cellMaskedWeights => null() integer, dimension(:,:,:), pointer :: sourceMaskedCells => null() end type remap_info_type type target_field_type character (len=64) :: name integer :: ndims = -1 integer :: xtype = -1 logical :: isTimeDependent = .false. integer, dimension(:), pointer :: dimlens => null() character (len=64), dimension(:), pointer :: dimnames => null() ! Members to store field data real :: array0r real, dimension(:), pointer :: array1r => null() real, dimension(:,:), pointer :: array2r => null() real, dimension(:,:,:), pointer :: array3r => null() real, dimension(:,:,:,:), pointer :: array4r => null() double precision :: array0d double precision, dimension(:), pointer :: array1d => null() double precision, dimension(:,:), pointer :: array2d => null() double precision, dimension(:,:,:), pointer :: array3d => null() double precision, dimension(:,:,:,:), pointer :: array4d => null() integer :: array0i integer, dimension(:), pointer :: array1i => null() integer, dimension(:,:), pointer :: array2i => null() integer, dimension(:,:,:), pointer :: array3i => null() integer, dimension(:,:,:,:), pointer :: array4i => null() end type target_field_type private :: nearest_cell, & nearest_vertex, & sphere_distance, & mpas_arc_length, & mpas_triangle_signed_area_sphere, & mpas_wachspress_coordinates, & convert_lx, & index2d contains integer function remap_info_setup(src_mesh, dst_mesh, remap_info) result(stat) implicit none type (mpas_mesh_type), intent(in), target :: src_mesh type (target_mesh_type), intent(in), target :: dst_mesh type (remap_info_type), intent(out) :: remap_info integer :: idx integer :: j integer :: nn integer :: ix, iy integer :: irank real :: sumWeights real, dimension(:,:), allocatable :: vertCoords real, dimension(3) :: pointInterp stat = 0 remap_info % method = 1 remap_info % src_mesh => src_mesh remap_info % dst_mesh => dst_mesh ! ! For nearest neighbor ! allocate(remap_info % nearestCell(dst_mesh % nlon, dst_mesh % nlat)) allocate(remap_info % nearestVertex(dst_mesh % nlon, dst_mesh % nlat)) allocate(remap_info % nearestEdge(dst_mesh % nlon, dst_mesh % nlat)) irank = dst_mesh % irank idx = 1 do iy=1,dst_mesh % nlat do ix=1,dst_mesh % nlon idx = nearest_cell(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & src_mesh % nCells, src_mesh % maxEdges, & src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, src_mesh % latCell, src_mesh % lonCell) remap_info % nearestCell(ix, iy) = idx end do end do idx = 1 do iy=1,dst_mesh % nlat do ix=1,dst_mesh % nlon idx = nearest_vertex(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & src_mesh % nCells, src_mesh % nVertices, src_mesh % maxEdges, 3, & src_mesh % nEdgesOnCell, src_mesh % verticesOnCell, & src_mesh % cellsOnVertex, src_mesh % latCell, src_mesh % lonCell, & src_mesh % latVertex, src_mesh % lonVertex ) remap_info % nearestVertex(ix, iy) = idx end do end do idx = 1 do iy=1,dst_mesh % nlat do ix=1,dst_mesh % nlon idx = nearest_vertex(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & src_mesh % nCells, src_mesh % nEdges, src_mesh % maxEdges, 2, & src_mesh % nEdgesOnCell, src_mesh % edgesOnCell, & src_mesh % cellsOnEdge, src_mesh % latCell, src_mesh % lonCell, & src_mesh % latEdge, src_mesh % lonEdge ) remap_info % nearestEdge(ix, iy) = idx end do end do ! ! For Wachspress interpolation ! allocate(vertCoords(3,3)) allocate(remap_info % cellWeights(3, dst_mesh % nlon, dst_mesh % nlat)) allocate(remap_info % sourceCells(3, dst_mesh % nlon, dst_mesh % nlat)) idx = 1 do iy=1,dst_mesh % nlat do ix=1,dst_mesh % nlon idx = nearest_vertex(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & src_mesh % nCells, src_mesh % nVertices, src_mesh % maxEdges, 3, & src_mesh % nEdgesOnCell, src_mesh % verticesOnCell, & src_mesh % cellsOnVertex, src_mesh % latCell, src_mesh % lonCell, & src_mesh % latVertex, src_mesh % lonVertex ) remap_info % sourceCells(:,ix,iy) = src_mesh % cellsOnVertex(:,idx) pointInterp(:) = convert_lx(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), 6371229.0) do j=1,3 vertCoords(:,j) = convert_lx(src_mesh % latCell(src_mesh % cellsOnVertex(j,idx)), & src_mesh % lonCell(src_mesh % cellsOnVertex(j,idx)), & 6371229.0) end do remap_info % cellWeights(:,ix,iy) = mpas_wachspress_coordinates(3, vertCoords, pointInterp) end do end do deallocate(vertCoords) allocate(vertCoords(3,src_mesh % maxEdges)) allocate(remap_info % vertexWeights(src_mesh % maxEdges, dst_mesh % nlon, dst_mesh % nlat)) allocate(remap_info % sourceVertices(src_mesh % maxEdges, dst_mesh % nlon, dst_mesh % nlat)) idx = 1 do iy=1,dst_mesh % nlat do ix=1,dst_mesh % nlon idx = nearest_cell(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & src_mesh % nCells, src_mesh % maxEdges, & src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, src_mesh % latCell, src_mesh % lonCell) nn = src_mesh % nEdgesOnCell(idx) remap_info % sourceVertices(:,ix,iy) = 1 remap_info % sourceVertices(1:nn,ix,iy) = src_mesh % verticesOnCell(1:nn,idx) pointInterp(:) = convert_lx(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), 6371229.0) do j=1,nn vertCoords(:,j) = convert_lx(src_mesh % latVertex(src_mesh % verticesOnCell(j,idx)), & src_mesh % lonVertex(src_mesh % verticesOnCell(j,idx)), & 6371229.0) end do remap_info % vertexWeights(:,ix,iy) = 0.0 remap_info % vertexWeights(1:nn,ix,iy) = mpas_wachspress_coordinates(3, vertCoords(:,1:nn), pointInterp) end do end do deallocate(vertCoords) allocate(vertCoords(3,src_mesh % maxEdges)) allocate(remap_info % edgeWeights(src_mesh % maxEdges, dst_mesh % nlon, dst_mesh % nlat)) allocate(remap_info % sourceEdges(src_mesh % maxEdges, dst_mesh % nlon, dst_mesh % nlat)) idx = 1 do iy=1,dst_mesh % nlat do ix=1,dst_mesh % nlon idx = nearest_cell(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & src_mesh % nCells, src_mesh % maxEdges, & src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, src_mesh % latCell, src_mesh % lonCell) nn = src_mesh % nEdgesOnCell(idx) remap_info % sourceEdges(:,ix,iy) = 1 remap_info % sourceEdges(1:nn,ix,iy) = src_mesh % edgesOnCell(1:nn,idx) pointInterp(:) = convert_lx(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), 6371229.0) do j=1,nn vertCoords(:,j) = convert_lx(src_mesh % latEdge(src_mesh % edgesOnCell(j,idx)), & src_mesh % lonEdge(src_mesh % edgesOnCell(j,idx)), & 6371229.0) end do remap_info % edgeWeights(:,ix,iy) = 0.0 remap_info % edgeWeights(1:nn,ix,iy) = mpas_wachspress_coordinates(3, vertCoords(:,1:nn), pointInterp) end do end do deallocate(vertCoords) ! ! For masked interpolation ! allocate(vertCoords(3,3)) allocate(remap_info % cellMaskedWeights(3, dst_mesh % nlon, dst_mesh % nlat)) allocate(remap_info % sourceMaskedCells(3, dst_mesh % nlon, dst_mesh % nlat)) idx = 1 do iy=1,dst_mesh % nlat do ix=1,dst_mesh % nlon idx = nearest_vertex(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), idx, & src_mesh % nCells, src_mesh % nVertices, src_mesh % maxEdges, 3, & src_mesh % nEdgesOnCell, src_mesh % verticesOnCell, & src_mesh % cellsOnVertex, src_mesh % latCell, src_mesh % lonCell, & src_mesh % latVertex, src_mesh % lonVertex ) remap_info % sourceMaskedCells(:,ix,iy) = src_mesh % cellsOnVertex(:,idx) pointInterp(:) = convert_lx(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), 6371229.0) do j=1,3 vertCoords(:,j) = convert_lx(src_mesh % latCell(src_mesh % cellsOnVertex(j,idx)), & src_mesh % lonCell(src_mesh % cellsOnVertex(j,idx)), & 6371229.0) end do remap_info % cellMaskedWeights(:,ix,iy) = mpas_wachspress_coordinates(3, vertCoords, pointInterp) sumWeights = 0.0 do j=1,3 if (src_mesh % landmask(remap_info % sourceMaskedCells(j,ix,iy)) == 1) then sumWeights = sumWeights + remap_info % cellMaskedWeights(j,ix,iy) else remap_info % cellMaskedWeights(j,ix,iy) = 0.0 end if end do if (sumWeights > 0.0) then remap_info % cellMaskedWeights(:,ix,iy) = remap_info % cellMaskedWeights(:,ix,iy) / sumWeights else idx = nearest_cell(dst_mesh % lats(index2d(irank,ix),iy), dst_mesh % lons(ix,index2d(irank,iy)), & remap_info % sourceMaskedCells(1,ix,iy), & src_mesh % nCells, src_mesh % maxEdges, & src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, src_mesh % latCell, src_mesh % lonCell) call search_for_cells(src_mesh % nCells, src_mesh % maxEdges, src_mesh % nEdgesOnCell, src_mesh % cellsOnCell, & src_mesh % landmask, dst_mesh % lats(index2d(irank,ix),iy), & dst_mesh % lons(ix,index2d(irank,iy)), & src_mesh % latCell, src_mesh % lonCell, idx, remap_info % sourceMaskedCells(:,ix,iy), & remap_info % cellMaskedWeights(:,ix,iy)) end if end do end do deallocate(vertCoords) end function remap_info_setup integer function remap_info_free(remap_info) result(stat) implicit none type (remap_info_type), intent(inout) :: remap_info stat = 0 remap_info % method = -1 nullify(remap_info % src_mesh) nullify(remap_info % dst_mesh) if (associated(remap_info % nearestCell)) then deallocate(remap_info % nearestCell) end if if (associated(remap_info % nearestVertex)) then deallocate(remap_info % nearestVertex) end if if (associated(remap_info % nearestEdge)) then deallocate(remap_info % nearestEdge) end if if (associated(remap_info % cellWeights)) then deallocate(remap_info % cellWeights) end if if (associated(remap_info % vertexWeights)) then deallocate(remap_info % vertexWeights) end if if (associated(remap_info % edgeWeights)) then deallocate(remap_info % edgeWeights) end if if (associated(remap_info % sourceCells)) then deallocate(remap_info % sourceCells) end if if (associated(remap_info % sourceVertices)) then deallocate(remap_info % sourceVertices) end if if (associated(remap_info % sourceEdges)) then deallocate(remap_info % sourceEdges) end if if (associated(remap_info % cellMaskedWeights)) then deallocate(remap_info % cellMaskedWeights) end if if (associated(remap_info % sourceMaskedCells)) then deallocate(remap_info % sourceMaskedCells) end if end function remap_info_free logical function can_remap_field(field) implicit none type (input_field_type), intent(in) :: field integer :: decompDim can_remap_field = .true. if (field % xtype /= FIELD_TYPE_INTEGER .and. & field % xtype /= FIELD_TYPE_REAL .and. & field % xtype /= FIELD_TYPE_DOUBLE) then can_remap_field = .false. return end if if (field % ndims == 0 .or. & (field % ndims == 1 .and. field % isTimeDependent)) then can_remap_field = .false. return end if decompDim = field % ndims if (field % isTimeDependent) then decompDim = decompDim - 1 end if if (trim(field % dimnames(decompDim)) /= 'nCells' .and. & trim(field % dimnames(decompDim)) /= 'nVertices' .and. & trim(field % dimnames(decompDim)) /= 'nEdges') then can_remap_field = .false. return end if end function can_remap_field integer function remap_field_dryrun(remap_info, src_field, dst_field) result(stat) implicit none type (remap_info_type), intent(in) :: remap_info type (input_field_type), intent(in) :: src_field type (target_field_type), intent(out) :: dst_field integer :: idim stat = 0 dst_field % name = src_field % name dst_field % xtype = src_field % xtype if (src_field % isTimeDependent) then ! Single horizontal dimension becomes two horizontal dimensions, nlat and nlon, ! but the time dimension is not counted in the target field dst_field % ndims = src_field % ndims else ! Single horizontal dimension becomes two horizontal dimensions, nlat and nlon dst_field % ndims = src_field % ndims + 1 end if allocate(dst_field % dimnames(dst_field % ndims)) allocate(dst_field % dimlens(dst_field % ndims)) dst_field % isTimeDependent = src_field % isTimeDependent do idim=1,dst_field % ndims-2 dst_field % dimlens(idim) = src_field % dimlens(idim) dst_field % dimnames(idim) = src_field % dimnames(idim) end do dst_field % dimlens(dst_field % ndims-1) = remap_info % dst_mesh % nlon dst_field % dimnames(dst_field % ndims-1) = 'nLons' dst_field % dimlens(dst_field % ndims) = remap_info % dst_mesh % nlat dst_field % dimnames(dst_field % ndims) = 'nLats' end function remap_field_dryrun integer function remap_field(remap_info, src_field, dst_field, masked) result(stat) implicit none type (remap_info_type), intent(in) :: remap_info type (input_field_type), intent(in) :: src_field type (target_field_type), intent(out) :: dst_field logical, intent(in), optional :: masked integer :: decompDim integer :: idim integer :: j integer :: iy, ix logical :: local_masked integer, dimension(:,:), pointer :: nearestIndex integer, dimension(:,:,:), pointer :: sourceNodes real, dimension(:,:,:), pointer :: nodeWeights stat = 0 decompDim = src_field % ndims if (src_field % isTimeDependent) then decompDim = decompDim - 1 end if local_masked = .false. if (present(masked)) then local_masked = masked end if if (trim(src_field % dimnames(decompDim)) == 'nCells') then nearestIndex => remap_info % nearestCell if (local_masked) then sourceNodes => remap_info % sourceMaskedCells nodeWeights => remap_info % cellMaskedWeights else sourceNodes => remap_info % sourceCells nodeWeights => remap_info % cellWeights end if else if (trim(src_field % dimnames(decompDim)) == 'nVertices') then nearestIndex => remap_info % nearestVertex sourceNodes => remap_info % sourceVertices nodeWeights => remap_info % vertexWeights else if (trim(src_field % dimnames(decompDim)) == 'nEdges') then nearestIndex => remap_info % nearestEdge sourceNodes => remap_info % sourceEdges nodeWeights => remap_info % edgeWeights else write(0,*) 'Remap exception: unhandled decomposed dim' stat = 1 return end if dst_field % name = src_field % name dst_field % xtype = src_field % xtype if (src_field % isTimeDependent) then ! Single horizontal dimension becomes two horizontal dimensions, nlat and nlon, ! but the time dimension is not counted in the target field dst_field % ndims = src_field % ndims else ! Single horizontal dimension becomes two horizontal dimensions, nlat and nlon dst_field % ndims = src_field % ndims + 1 end if allocate(dst_field % dimnames(dst_field % ndims)) allocate(dst_field % dimlens(dst_field % ndims)) dst_field % isTimeDependent = src_field % isTimeDependent do idim=1,dst_field % ndims-2 dst_field % dimlens(idim) = src_field % dimlens(idim) dst_field % dimnames(idim) = src_field % dimnames(idim) end do dst_field % dimlens(dst_field % ndims-1) = remap_info % dst_mesh % nlon dst_field % dimnames(dst_field % ndims-1) = 'nLons' dst_field % dimlens(dst_field % ndims) = remap_info % dst_mesh % nlat dst_field % dimnames(dst_field % ndims) = 'nLats' if (src_field % xtype == FIELD_TYPE_REAL) then if (dst_field % ndims == 2) then allocate(dst_field % array2r(dst_field % dimlens(1), & dst_field % dimlens(2))) #ifdef NEAREST_NEIGHBOR do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array2r(ix,iy) = src_field % array1r(nearestIndex(ix,iy)) end do end do #else do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array2r(ix,iy) = 0.0 do j=1,size(sourceNodes, dim=1) dst_field % array2r(ix,iy) = dst_field % array2r(ix,iy) + & nodeWeights(j,ix,iy) & * src_field % array1r(sourceNodes(j,ix,iy)) end do end do end do #endif else if (dst_field % ndims == 3) then allocate(dst_field % array3r(dst_field % dimlens(1), & dst_field % dimlens(2), & dst_field % dimlens(3))) #ifdef NEAREST_NEIGHBOR do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array3r(:,ix,iy) = src_field % array2r(:,nearestIndex(ix,iy)) end do end do #else do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array3r(:,ix,iy) = 0.0 do j=1,3 dst_field % array3r(:,ix,iy) = dst_field % array3r(:,ix,iy) + & nodeWeights(j,ix,iy) & * src_field % array2r(:,sourceNodes(j,ix,iy)) end do end do end do #endif else if (dst_field % ndims == 4) then allocate(dst_field % array4r(dst_field % dimlens(1), & dst_field % dimlens(2), & dst_field % dimlens(3), & dst_field % dimlens(4))) #ifdef NEAREST_NEIGHBOR do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array4r(:,:,ix,iy) = src_field % array3r(:,:,nearestIndex(ix,iy)) end do end do #endif do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array4r(:,:,ix,iy) = 0.0 do j=1,3 dst_field % array4r(:,:,ix,iy) = dst_field % array4r(:,:,ix,iy) + & remap_info % cellWeights(j,ix,iy) & * src_field % array3r(:,:,remap_info % sourceCells(j,ix,iy)) end do end do end do else write(0,*) 'Remap exception: unhandled dimension for real ', dst_field % ndims end if else if (src_field % xtype == FIELD_TYPE_DOUBLE) then if (dst_field % ndims == 2) then allocate(dst_field % array2d(dst_field % dimlens(1), & dst_field % dimlens(2))) #ifdef NEAREST_NEIGHBOR do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array2d(ix,iy) = src_field % array1d(nearestIndex(ix,iy)) end do end do #endif do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array2d(ix,iy) = 0.0 do j=1,size(sourceNodes, dim=1) dst_field % array2d(ix,iy) = dst_field % array2d(ix,iy) + & nodeWeights(j,ix,iy) & * src_field % array1d(sourceNodes(j,ix,iy)) end do end do end do else if (dst_field % ndims == 3) then allocate(dst_field % array3d(dst_field % dimlens(1), & dst_field % dimlens(2), & dst_field % dimlens(3))) #ifdef NEAREST_NEIGHBOR do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array3d(:,ix,iy) = src_field % array2d(:,nearestIndex(ix,iy)) end do end do #endif do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array3d(:,ix,iy) = 0.0 do j=1,3 dst_field % array3d(:,ix,iy) = dst_field % array3d(:,ix,iy) + & remap_info % cellWeights(j,ix,iy) & * src_field % array2d(:,remap_info % sourceCells(j,ix,iy)) end do end do end do else if (dst_field % ndims == 4) then allocate(dst_field % array4d(dst_field % dimlens(1), & dst_field % dimlens(2), & dst_field % dimlens(3), & dst_field % dimlens(4))) #ifdef NEAREST_NEIGHBOR do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array4d(:,:,ix,iy) = src_field % array3d(:,:,nearestIndex(ix,iy)) end do end do #endif do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array4d(:,:,ix,iy) = 0.0 do j=1,3 dst_field % array4d(:,:,ix,iy) = dst_field % array4d(:,:,ix,iy) + & remap_info % cellWeights(j,ix,iy) & * src_field % array3d(:,:,remap_info % sourceCells(j,ix,iy)) end do end do end do else write(0,*) 'Remap exception: unhandled dimension for dbl ', dst_field % ndims end if else if (src_field % xtype == FIELD_TYPE_INTEGER) then if (dst_field % ndims == 2) then allocate(dst_field % array2i(dst_field % dimlens(1), & dst_field % dimlens(2))) do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array2i(ix,iy) = src_field % array1i(nearestIndex(ix,iy)) end do end do else if (dst_field % ndims == 3) then allocate(dst_field % array3i(dst_field % dimlens(1), & dst_field % dimlens(2), & dst_field % dimlens(3))) do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array3i(:,ix,iy) = src_field % array2i(:,nearestIndex(ix,iy)) end do end do else if (dst_field % ndims == 4) then allocate(dst_field % array4i(dst_field % dimlens(1), & dst_field % dimlens(2), & dst_field % dimlens(3), & dst_field % dimlens(4))) do iy=1,size(nearestIndex, dim=2) do ix=1,size(nearestIndex, dim=1) dst_field % array4i(:,:,ix,iy) = src_field % array3i(:,:,nearestIndex(ix,iy)) end do end do else write(0,*) 'Remap exception: unhandled dimension for int ', dst_field % ndims end if else write(0,*) 'Remap exception: unhandled type' end if end function remap_field integer function remap_get_target_latitudes(remap_info, lat_field) result(stat) implicit none type (remap_info_type), intent(in) :: remap_info type (target_field_type), intent(out) :: lat_field real, parameter :: rad2deg = 90.0 / asin(1.0) stat = 0 lat_field % name = 'lat' lat_field % xtype = FIELD_TYPE_REAL lat_field % ndims = 1 lat_field % isTimeDependent = .false. allocate(lat_field % dimnames(lat_field % ndims)) allocate(lat_field % dimlens(lat_field % ndims)) lat_field % dimnames(1) = 'nLats' lat_field % dimlens(1) = remap_info % dst_mesh % nlat allocate(lat_field % array1r(lat_field % dimlens(1))) lat_field % array1r(:) = remap_info % dst_mesh % lats(1,:) * rad2deg end function remap_get_target_latitudes integer function remap_get_target_longitudes(remap_info, lon_field) result(stat) implicit none type (remap_info_type), intent(in) :: remap_info type (target_field_type), intent(out) :: lon_field real, parameter :: rad2deg = 90.0 / asin(1.0) stat = 0 lon_field % name = 'lon' lon_field % xtype = FIELD_TYPE_REAL lon_field % ndims = 1 lon_field % isTimeDependent = .false. allocate(lon_field % dimnames(lon_field % ndims)) allocate(lon_field % dimlens(lon_field % ndims)) lon_field % dimnames(1) = 'nLons' lon_field % dimlens(1) = remap_info % dst_mesh % nlon allocate(lon_field % array1r(lon_field % dimlens(1))) lon_field % array1r(:) = remap_info % dst_mesh % lons(:,1) * rad2deg end function remap_get_target_longitudes integer function free_target_field(field) result(stat) implicit none type (target_field_type), intent(inout) :: field stat = 0 if (associated(field % dimlens)) then deallocate(field % dimlens) end if if (associated(field % dimnames)) then deallocate(field % dimnames) end if if (associated(field % array1r)) then deallocate(field % array1r) end if if (associated(field % array2r)) then deallocate(field % array2r) end if if (associated(field % array3r)) then deallocate(field % array3r) end if if (associated(field % array4r)) then deallocate(field % array4r) end if if (associated(field % array1d)) then deallocate(field % array1d) end if if (associated(field % array2d)) then deallocate(field % array2d) end if if (associated(field % array3d)) then deallocate(field % array3d) end if if (associated(field % array4d)) then deallocate(field % array4d) end if if (associated(field % array1i)) then deallocate(field % array1i) end if if (associated(field % array2i)) then deallocate(field % array2i) end if if (associated(field % array3i)) then deallocate(field % array3i) end if if (associated(field % array4i)) then deallocate(field % array4i) end if end function free_target_field integer function nearest_cell(target_lat, target_lon, start_cell, nCells, maxEdges, & nEdgesOnCell, cellsOnCell, latCell, lonCell) implicit none real, intent(in) :: target_lat, target_lon integer, intent(in) :: start_cell integer, intent(in) :: nCells, maxEdges integer, dimension(nCells), intent(in) :: nEdgesOnCell integer, dimension(maxEdges,nCells), intent(in) :: cellsOnCell real, dimension(nCells), intent(in) :: latCell, lonCell integer :: i integer :: iCell integer :: current_cell real :: current_distance, d real :: nearest_distance nearest_cell = start_cell current_cell = -1 do while (nearest_cell /= current_cell) current_cell = nearest_cell current_distance = sphere_distance(latCell(current_cell), lonCell(current_cell), target_lat, & target_lon, 1.0) nearest_cell = current_cell nearest_distance = current_distance do i = 1, nEdgesOnCell(current_cell) iCell = cellsOnCell(i,current_cell) if (iCell <= nCells) then d = sphere_distance(latCell(iCell), lonCell(iCell), target_lat, target_lon, 1.0) if (d < nearest_distance) then nearest_cell = iCell nearest_distance = d end if end if end do end do end function nearest_cell integer function nearest_vertex( target_lat, target_lon, & start_vertex, & nCells, nVertices, maxEdges, vertexDegree, & nEdgesOnCell, verticesOnCell, & cellsOnVertex, latCell, lonCell, & latVertex, lonVertex ) implicit none real, intent(in) :: target_lat, target_lon integer, intent(in) :: start_vertex integer, intent(in) :: nCells, nVertices, maxEdges, vertexDegree integer, dimension(nCells), intent(in) :: nEdgesOnCell integer, dimension(maxEdges,nCells), intent(in) :: verticesOnCell integer, dimension(vertexDegree,nVertices), intent(in) :: cellsOnVertex real, dimension(nCells), intent(in) :: latCell, lonCell real, dimension(nVertices), intent(in) :: latVertex, lonVertex integer :: i, cell1, cell2, cell3, iCell integer :: iVtx integer :: current_vertex real :: cell1_dist, cell2_dist, cell3_dist real :: current_distance, d real :: nearest_distance nearest_vertex = start_vertex current_vertex = -1 do while (nearest_vertex /= current_vertex) current_vertex = nearest_vertex current_distance = sphere_distance(latVertex(current_vertex), lonVertex(current_vertex), & target_lat, target_lon, 1.0) nearest_vertex = current_vertex nearest_distance = current_distance cell1 = cellsOnVertex(1,current_vertex) cell1_dist = sphere_distance(latCell(cell1), lonCell(cell1), target_lat, target_lon, 1.0) cell2 = cellsOnVertex(2,current_vertex) cell2_dist = sphere_distance(latCell(cell2), lonCell(cell2), target_lat, target_lon, 1.0) if (vertexDegree == 3) then cell3 = cellsOnVertex(3,current_vertex) cell3_dist = sphere_distance(latCell(cell3), lonCell(cell3), target_lat, target_lon, 1.0) end if if (vertexDegree == 3) then if (cell1_dist < cell2_dist) then if (cell1_dist < cell3_dist) then iCell = cell1 else iCell = cell3 end if else if (cell2_dist < cell3_dist) then iCell = cell2 else iCell = cell3 end if end if else if (cell1_dist < cell2_dist) then iCell = cell1 else iCell = cell2 end if end if do i = 1, nEdgesOnCell(iCell) iVtx = verticesOnCell(i,iCell) d = sphere_distance(latVertex(iVtx), lonVertex(iVtx), target_lat, target_lon, 1.0) if (d < nearest_distance) then nearest_vertex = iVtx nearest_distance = d end if end do end do end function nearest_vertex real function sphere_distance(lat1, lon1, lat2, lon2, radius) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a ! sphere with given radius. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! implicit none real, intent(in) :: lat1, lon1, lat2, lon2, radius real :: arg1 arg1 = sqrt( sin(0.5*(lat2-lat1))**2 + & cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 ) sphere_distance = 2.*radius*asin(arg1) end function sphere_distance !*********************************************************************** ! ! function mpas_wachspress_coordinates ! !> \brief Compute the barycentric Wachspress coordinates for a polygon !> \author Phillip Wolfram !> \date 01/26/2015 !> \details !> Computes the barycentric Wachspress coordinates for a polygon with nVertices !> points in R3, vertCoords for a particular pointInterp with normalized radius. !> Follows Gillette, A., Rand, A., Bajaj, C., 2011. !> Error estimates for generalized barycentric interpolation. !> Advances in computational mathematics 37 (3), 417–439. !> Optimized version of mpas_wachspress_coordinates uses optional cached B_i areas !------------------------------------------------------------------------ function mpas_wachspress_coordinates(nVertices, vertCoords, pointInterp, areaBin) implicit none ! input points integer, intent(in) :: nVertices real, dimension(3, nVertices), intent(in) :: vertCoords real, dimension(3), intent(in) :: pointInterp real, dimension(nVertices), optional, intent(in) :: areaBin ! output real, dimension(nVertices) :: mpas_wachspress_coordinates ! computational intermediates real, dimension(nVertices) :: wach ! The wachpress area-product real :: wach_total ! The wachpress total weight integer :: i, j ! Loop indices integer :: im1, i0, ip1 ! im1 = (i-1), i0 = i, ip1 = (i+1) ! triangle areas to compute wachspress coordinate real, dimension(nVertices) :: areaA real, dimension(nVertices) :: areaB real :: radiusLocal radiusLocal = sqrt(sum(vertCoords(:,1)**2)) if (.not. present(areaBin)) then ! compute areas do i = 1, nVertices ! compute first area B_i ! get vertex indices im1 = mod(nVertices + i - 2, nVertices) + 1 i0 = mod(nVertices + i - 1, nVertices) + 1 ip1 = mod(nVertices + i , nVertices) + 1 ! precompute B_i areas ! always the same because B_i independent of xp,yp,zp ! (COULD CACHE AND USE RESULT FROM ARRAY FOR FURTHER OPTIMIZATION) areaB(i) = mpas_triangle_signed_area_sphere(vertCoords(:, im1), vertCoords(:, i0), vertCoords(:, ip1), radiusLocal) end do else ! assign areas do i = 1, nVertices areaB(i) = areaBin(i) end do end if ! compute areas do i = 1, nVertices ! compute first area B_i ! get vertex indices im1 = mod(nVertices + i - 2, nVertices) + 1 i0 = mod(nVertices + i - 1, nVertices) + 1 ip1 = mod(nVertices + i , nVertices) + 1 ! compute A_ij areas ! must be computed each time areaA(i0) = mpas_triangle_signed_area_sphere(pointInterp, vertCoords(:, i0), vertCoords(:, ip1), radiusLocal) ! precomputed B_i areas, cached end do ! for each vertex compute wachpress coordinate do i = 1, nVertices wach(i) = areaB(i) do j = (i + 1), (i + nVertices - 2) i0 = mod(nVertices + j - 1, nVertices) + 1 ! accumulate products for A_ij subareas wach(i) = wach(i) * areaA(i0) end do end do ! get summed weights for normalization wach_total = 0 do i = 1, nVertices wach_total = wach_total + wach(i) end do ! compute lambda mpas_wachspress_coordinates= 0.0 do i = 1, nVertices mpas_wachspress_coordinates(i) = wach(i)/wach_total end do end function mpas_wachspress_coordinates !*********************************************************************** ! ! routine mpas_triangle_signed_area_sphere ! !> \brief Calculates area of a triangle on a sphere !> \author Matthew Hoffman !> \date 13 January 2015 !> \details !> This routine calculates the area of a triangle on the surface of a sphere. !> Uses the spherical analog of Heron's formula. !> Copied from mesh generator. A CCW winding angle is positive. !----------------------------------------------------------------------- real function mpas_triangle_signed_area_sphere(a, b, c, radius) implicit none !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- real, dimension(3), intent(in) :: a, b, c !< Input: 3d (x,y,z) points forming the triangle in which to calculate the bary weights real, intent(in) :: radius !< sphere radius !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- real :: ab, bc, ca, semiperim, tanqe real, dimension(3) :: ablen, aclen, Dlen ab = mpas_arc_length(a(1), a(2), a(3), b(1), b(2), b(3))/radius bc = mpas_arc_length(b(1), b(2), b(3), c(1), c(2), c(3))/radius ca = mpas_arc_length(c(1), c(2), c(3), a(1), a(2), a(3))/radius semiperim = 0.5 * (ab + bc + ca) tanqe = sqrt(max(0.0,tan(0.5 * semiperim) * tan(0.5 * (semiperim - ab)) & * tan(0.5 * (semiperim - bc)) * tan(0.5 * (semiperim - ca)))) mpas_triangle_signed_area_sphere = 4.0 * radius * radius * atan(tanqe) ! computing correct signs (in similar fashion to mpas_sphere_angle) ablen(1) = b(1) - a(1) ablen(2) = b(2) - a(2) ablen(3) = b(3) - a(3) aclen(1) = c(1) - a(1) aclen(2) = c(2) - a(2) aclen(3) = c(3) - a(3) dlen(1) = (ablen(2) * aclen(3)) - (ablen(3) * aclen(2)) dlen(2) = -((ablen(1) * aclen(3)) - (ablen(3) * aclen(1))) dlen(3) = (ablen(1) * aclen(2)) - (ablen(2) * aclen(1)) if ((Dlen(1)*a(1) + Dlen(2)*a(2) + Dlen(3)*a(3)) < 0.0) then mpas_triangle_signed_area_sphere = -mpas_triangle_signed_area_sphere end if end function mpas_triangle_signed_area_sphere !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! FUNCTION MPAS_ARC_LENGTH ! ! Returns the length of the great circle arc from A=(ax, ay, az) to ! B=(bx, by, bz). It is assumed that both A and B lie on the surface of the ! same sphere centered at the origin. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real function mpas_arc_length(ax, ay, az, bx, by, bz) implicit none real, intent(in) :: ax, ay, az, bx, by, bz real :: r, c real :: cx, cy, cz cx = bx - ax cy = by - ay cz = bz - az r = sqrt(ax*ax + ay*ay + az*az) c = sqrt(cx*cx + cy*cy + cz*cz) mpas_arc_length = r * 2.0 * asin(c/(2.0*r)) end function mpas_arc_length function convert_lx(lat, lon, radius) result(vec) implicit none real, intent(in) :: lat, lon, radius real, dimension(3) :: vec vec(1) = radius * cos(lon) * cos(lat) vec(2) = radius * sin(lon) * cos(lat) vec(3) = radius * sin(lat) end function convert_lx subroutine search_for_cells(nCells, maxEdges, nEdgesOnCell, cellsOnCell, cellMask, & targetLat, targetLon, latCell, lonCell, & startIdx, sourceMaskedCells, cellMaskedWeights) implicit none integer, intent(in) :: nCells integer, intent(in) :: maxEdges integer, dimension(nCells), intent(in) :: nEdgesOnCell integer, dimension(maxEdges,nCells), intent(in) :: cellsOnCell integer, dimension(nCells), intent(in) :: cellMask real, intent(in) :: targetLat real, intent(in) :: targetLon real, dimension(nCells), intent(in) :: latCell real, dimension(nCells), intent(in) :: lonCell integer, intent(in) :: startIdx integer, dimension(3), intent(inout) :: sourceMaskedCells real, dimension(3), intent(inout) :: cellMaskedWeights integer :: i integer :: scan_cell integer :: neighbor integer :: unscanned logical :: no_more_queueing real :: d, d_min ! ! Reset data structures ! call queue_reset() call dictionary_reset() no_more_queueing = .false. unscanned = 0 d_min = 1000.0 sourceMaskedCells(:) = 1 cellMaskedWeights(:) = 0.0 ! ! Insert the origin cell into the queue for processing ! call queue_insert(startIdx) call dictionary_insert(startIdx) do while (queue_size > 0) scan_cell = queue_remove() ! ! Each cell index removed from the queue represents a unique grid cell ! that falls within the specified radius of the origin point ! Here, we can do any processing we like for these cells ! if (cellMask(scan_cell) == 1) then d = sphere_distance(targetLat, targetLon, latCell(scan_cell), lonCell(scan_cell), 1.0) if (d < d_min) then d_min = d sourceMaskedCells(1) = scan_cell cellMaskedWeights(1) = 1.0 if (.not. no_more_queueing) then unscanned = queue_size end if no_more_queueing = .true. end if end if ! ! Add any neighbors of scan_cell within specified radius to the queue for processing ! if (.not. no_more_queueing .or. unscanned > 0) then do i=1,nEdgesOnCell(scan_cell) neighbor = cellsOnCell(i,scan_cell) if (.not. dictionary_search(neighbor)) then call queue_insert(neighbor) call dictionary_insert(neighbor) end if end do end if if (no_more_queueing) then unscanned = unscanned - 1 end if end do end subroutine search_for_cells !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Insert a new integer into the queue !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine queue_insert(i) implicit none integer, intent(in) :: i if (queue_size == max_queue_length) then write(0,*) 'Error: queue overrun' return end if queue_size = queue_size + 1 queue_array(queue_head) = i queue_head = mod(queue_head + 1, max_queue_length) end subroutine queue_insert !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Remove the oldest integer from the queue !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! integer function queue_remove() implicit none if (queue_size <= 0) then write(0,*) 'Error: queue underrun' queue_remove = -1 return end if queue_size = queue_size - 1 queue_remove = queue_array(queue_tail) queue_tail = mod(queue_tail + 1, max_queue_length) end function queue_remove !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Reset the queue to an empty state !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine queue_reset() implicit none queue_head = 0 queue_tail = 0 queue_size = 0 end subroutine queue_reset !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Insert an integer into the dictionary !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine dictionary_insert(i) implicit none integer, intent(in) :: i integer :: n_integer integer :: n_bit n_integer = ((i-1) / int_size) + 1 n_bit = mod((i-1), int_size) if (n_integer > max_dictionary_size) then write(0,*) 'Error: dictionary insert out of bounds' return end if dictionary_array(n_integer) = ibset(dictionary_array(n_integer), n_bit) end subroutine dictionary_insert !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Search for an integer in the dictionary !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! logical function dictionary_search(i) implicit none integer, intent(in) :: i integer :: n_integer integer :: n_bit n_integer = ((i-1) / int_size) + 1 n_bit = mod((i-1), int_size) if (n_integer > max_dictionary_size) then write(0,*) 'Error: dictionary search out of bounds' dictionary_search = .false. return end if dictionary_search = btest(dictionary_array(n_integer), n_bit) end function dictionary_search !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Reset the dictionary to an empty state !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine dictionary_reset() implicit none dictionary_array(:) = 0 end subroutine dictionary_reset function index2d(irank, idx) result(i) implicit none integer, intent(in) :: irank, idx integer :: i i = irank * (idx - 1) + 1 end function index2d end module remapper