diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 2fa1ad8c39..bb3021cb4a 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -2182,210 +2182,210 @@ end subroutine SeedlingParPatch subroutine DetermineGridCellNeighbors(neighbors,seeds,numg) - ! This subroutine utilizes information from the decomposition and domain types to determine - ! the set of grid cell neighbors within some maximum distance. It records the distance for each - ! neighbor for later use. This should be called after decompInit_lnd and surf_get_grid - ! as it relies on ldecomp and ldomain information. - - use decompMod , only : procinfo - use domainMod , only : ldomain - use spmdMod , only : MPI_REAL8, MPI_INTEGER, mpicom, npes, masterproc, iam - use perf_mod , only : t_startf, t_stopf - use FatesDispersalMod , only : neighborhood_type, neighbor_type, ProbabilityDensity, dispersal_type - use FatesUtilsMod , only : GetNeighborDistance - use FatesConstantsMod , only : fates_unset_int - use EDPftvarcon , only : EDPftvarcon_inst + ! This subroutine utilizes information from the decomposition and domain types to determine + ! the set of grid cell neighbors within some maximum distance. It records the distance for each + ! neighbor for later use. This should be called after decompInit_lnd and surf_get_grid + ! as it relies on ldecomp and ldomain information. - ! Arguments - type(neighborhood_type), intent(inout), pointer :: neighbors(:) ! land gridcell neighbor data structure - type(dispersal_type), intent(inout) :: seeds ! land gridcell neighbor data structure - integer , intent(in) :: numg ! number of land gridcells + use decompMod , only : procinfo + use domainMod , only : ldomain + use spmdMod , only : MPI_REAL8, MPI_INTEGER, mpicom, npes, masterproc, iam + use perf_mod , only : t_startf, t_stopf + use FatesDispersalMod , only : neighborhood_type, neighbor_type, ProbabilityDensity, dispersal_type + use FatesUtilsMod , only : GetNeighborDistance + use FatesConstantsMod , only : fates_unset_int + use EDPftvarcon , only : EDPftvarcon_inst - ! Local variables - type (neighbor_type), pointer :: current_neighbor - type (neighbor_type), pointer :: another_neighbor + ! Arguments + type(neighborhood_type), intent(inout), pointer :: neighbors(:) ! land gridcell neighbor data structure + type(dispersal_type), intent(inout) :: seeds ! land gridcell neighbor data structure + integer , intent(in) :: numg ! number of land gridcells - integer :: i, gi, gj, ni ! indices - integer :: ier, mpierr ! error status - integer :: ipft ! pft index + ! Local variables + type (neighbor_type), pointer :: current_neighbor + type (neighbor_type), pointer :: another_neighbor - integer, allocatable :: ncells_array(:), begg_array(:) ! number of cells and starting global grid cell index per process - real(r8), allocatable :: gclat(:), gclon(:) ! local array holding gridcell lat and lon + integer :: i, gi, gj, ni ! indices + integer :: ier, mpierr ! error status + integer :: ipft ! pft index - real(r8) :: g2g_dist ! grid cell distance (m) - real(r8) :: pdf ! probability density function output + integer, allocatable :: ncells_array(:), begg_array(:) ! number of cells and starting global grid cell index per process + real(r8), allocatable :: gclat(:), gclon(:) ! local array holding gridcell lat and lon - if(debug .and. hlm_is_restart .eq. itrue) write(fates_log(),*) 'gridcell initialization during restart' + real(r8) :: g2g_dist ! grid cell distance (m) + real(r8) :: pdf ! probability density function output - if(debug) write(fates_log(),*)'DGCN: npes, numg: ', npes, numg + if(debug .and. hlm_is_restart .eq. itrue) write(fates_log(),*) 'gridcell initialization during restart' - ! Allocate and initialize array neighbor type - allocate(neighbors(numg), stat=ier) - neighbors(:)%neighbor_count = 0 + if(debug) write(fates_log(),*)'DGCN: npes, numg: ', npes, numg - ! Allocate and initialize local lat and lon arrays - allocate(gclat(numg), stat=ier) - if(debug) write(fates_log(),*)'DGCN: gclat alloc: ', ier + ! Allocate and initialize array neighbor type + allocate(neighbors(numg), stat=ier) + neighbors(:)%neighbor_count = 0 - allocate(gclon(numg), stat=ier) - if(debug) write(fates_log(),*)'DGCN: gclon alloc: ', ier + ! Allocate and initialize local lat and lon arrays + allocate(gclat(numg), stat=ier) + if(debug) write(fates_log(),*)'DGCN: gclat alloc: ', ier - gclon(:) = nan - gclat(:) = nan + allocate(gclon(numg), stat=ier) + if(debug) write(fates_log(),*)'DGCN: gclon alloc: ', ier - ! Allocate and initialize MPI count and displacement values - allocate(ncells_array(0:npes-1), stat=ier) - if(debug) write(fates_log(),*)'DGCN: ncells alloc: ', ier + gclon(:) = nan + gclat(:) = nan - allocate(begg_array(0:npes-1), stat=ier) - if(debug) write(fates_log(),*)'DGCN: begg alloc: ', ier + ! Allocate and initialize MPI count and displacement values + allocate(ncells_array(0:npes-1), stat=ier) + if(debug) write(fates_log(),*)'DGCN: ncells alloc: ', ier - ncells_array(:) = fates_unset_int - begg_array(:) = fates_unset_int + allocate(begg_array(0:npes-1), stat=ier) + if(debug) write(fates_log(),*)'DGCN: begg alloc: ', ier - call t_startf('fates-seed-init-allgather') + ncells_array(:) = fates_unset_int + begg_array(:) = fates_unset_int - if(debug) write(fates_log(),*)'DGCN: procinfo%begg: ', procinfo%begg - if(debug) write(fates_log(),*)'DGCN: procinfo%ncells: ', procinfo%ncells + call t_startf('fates-seed-init-allgather') - ! Gather the sizes of the ldomain that each mpi rank is passing - call MPI_Allgather(procinfo%ncells,1,MPI_INTEGER,ncells_array,1,MPI_INTEGER,mpicom,mpierr) - if(debug) write(fates_log(),*)'DGCN: ncells mpierr: ', mpierr + if(debug) write(fates_log(),*)'DGCN: procinfo%begg: ', procinfo%begg + if(debug) write(fates_log(),*)'DGCN: procinfo%ncells: ', procinfo%ncells - ! Gather the starting gridcell index for each ldomain - call MPI_Allgather(procinfo%begg,1,MPI_INTEGER,begg_array,1,MPI_INTEGER,mpicom,mpierr) - if(debug) write(fates_log(),*)'DGCN: begg mpierr: ', mpierr + ! Gather the sizes of the ldomain that each mpi rank is passing + call MPI_Allgather(procinfo%ncells,1,MPI_INTEGER,ncells_array,1,MPI_INTEGER,mpicom,mpierr) + if(debug) write(fates_log(),*)'DGCN: ncells mpierr: ', mpierr - ! reduce the begg_array displacements by one as MPI collectives expect zero indexed arrays - begg_array = begg_array - 1 + ! Gather the starting gridcell index for each ldomain + call MPI_Allgather(procinfo%begg,1,MPI_INTEGER,begg_array,1,MPI_INTEGER,mpicom,mpierr) + if(debug) write(fates_log(),*)'DGCN: begg mpierr: ', mpierr - if(debug) write(fates_log(),*)'DGCN: ncells_array: ' , ncells_array - if(debug) write(fates_log(),*)'DGCN: begg_array: ' , begg_array + ! reduce the begg_array displacements by one as MPI collectives expect zero indexed arrays + begg_array = begg_array - 1 - ! Gather the domain information together into the neighbor type - ! Note that MPI_Allgatherv is only gathering a subset of ldomain - if(debug) write(fates_log(),*)'DGCN: gathering latc' - call MPI_Allgatherv(ldomain%latc,procinfo%ncells,MPI_REAL8,gclat,ncells_array,begg_array,MPI_REAL8,mpicom,mpierr) + if(debug) write(fates_log(),*)'DGCN: ncells_array: ' , ncells_array + if(debug) write(fates_log(),*)'DGCN: begg_array: ' , begg_array - if(debug) write(fates_log(),*)'DGCN: gathering lonc' - call MPI_Allgatherv(ldomain%lonc,procinfo%ncells,MPI_REAL8,gclon,ncells_array,begg_array,MPI_REAL8,mpicom,mpierr) + ! Gather the domain information together into the neighbor type + ! Note that MPI_Allgatherv is only gathering a subset of ldomain + if(debug) write(fates_log(),*)'DGCN: gathering latc' + call MPI_Allgatherv(ldomain%latc,procinfo%ncells,MPI_REAL8,gclat,ncells_array,begg_array,MPI_REAL8,mpicom,mpierr) - if (debug .and. iam .eq. 0) then - write(fates_log(),*)'DGCN: sum(gclat):, sum(gclon): ', sum(gclat), sum(gclon) - end if + if(debug) write(fates_log(),*)'DGCN: gathering lonc' + call MPI_Allgatherv(ldomain%lonc,procinfo%ncells,MPI_REAL8,gclon,ncells_array,begg_array,MPI_REAL8,mpicom,mpierr) - ! Save number of cells and begging index arrays to dispersal type - if(debug) write(fates_log(),*)'DGCN: save to seeds type' - if(debug) write(fates_log(),*)'DGCN: seeds ncells alloc: ', allocated(seeds%ncells_array) - if(debug) write(fates_log(),*)'DGCN: seeds begg alloc: ', allocated(seeds%begg_array) - seeds%ncells_array = ncells_array - seeds%begg_array = begg_array + if (debug .and. iam .eq. 0) then + write(fates_log(),*)'DGCN: sum(gclat):, sum(gclon): ', sum(gclat), sum(gclon) + end if - if (debug .and. iam .eq. 0) then - write(fates_log(),*)'DGCN: seeds%ncells_array: ', seeds%ncells_array - write(fates_log(),*)'DGCN: seeds%begg_array: ', seeds%begg_array - end if + ! Save number of cells and begging index arrays to dispersal type + if(debug) write(fates_log(),*)'DGCN: save to seeds type' + if(debug) write(fates_log(),*)'DGCN: seeds ncells alloc: ', allocated(seeds%ncells_array) + if(debug) write(fates_log(),*)'DGCN: seeds begg alloc: ', allocated(seeds%begg_array) + seeds%ncells_array = ncells_array + seeds%begg_array = begg_array - call t_stopf('fates-seed-init-allgather') + if (debug .and. iam .eq. 0) then + write(fates_log(),*)'DGCN: seeds%ncells_array: ', seeds%ncells_array + write(fates_log(),*)'DGCN: seeds%begg_array: ', seeds%begg_array + end if - call t_startf('fates-seed-init-decomp') + call t_stopf('fates-seed-init-allgather') - if(debug) write(fates_log(), *) 'DGCN: maxdist: ', EDPftvarcon_inst%seed_dispersal_max_dist + call t_startf('fates-seed-init-decomp') - ! Iterate through the grid cell indices and determine if any neighboring cells are in range - gc_loop: do gi = 1,numg-1 + if(debug) write(fates_log(), *) 'DGCN: maxdist: ', EDPftvarcon_inst%seed_dispersal_max_dist - ! Seach forward through all indices for neighbors to current grid cell index - neighbor_search: do gj = gi+1,numg + ! Iterate through the grid cell indices and determine if any neighboring cells are in range + gc_loop: do gi = 1,numg-1 - ! Determine distance to old grid cells to the current one - g2g_dist = GetNeighborDistance(gi,gj,gclat,gclon) + ! Seach forward through all indices for neighbors to current grid cell index + neighbor_search: do gj = gi+1,numg - if(debug) write(fates_log(), *) 'DGCN: gi,gj,g2g_dist: ', gi,gj,g2g_dist + ! Determine distance to old grid cells to the current one + g2g_dist = GetNeighborDistance(gi,gj,gclat,gclon) - ! - dist_check: if (any(EDPftvarcon_inst%seed_dispersal_max_dist .gt. g2g_dist)) then + if(debug) write(fates_log(), *) 'DGCN: gi,gj,g2g_dist: ', gi,gj,g2g_dist - ! Add neighbor index to current grid cell index list - allocate(current_neighbor) - current_neighbor%next_neighbor => null() + ! + dist_check: if (any(EDPftvarcon_inst%seed_dispersal_max_dist .gt. g2g_dist)) then - current_neighbor%gindex = gj + ! Add neighbor index to current grid cell index list + allocate(current_neighbor) + current_neighbor%next_neighbor => null() - current_neighbor%gc_dist = g2g_dist + current_neighbor%gindex = gj - allocate(current_neighbor%density_prob(numpft)) + current_neighbor%gc_dist = g2g_dist - do ipft = 1, numpft - call ProbabilityDensity(pdf, ipft, g2g_dist) - current_neighbor%density_prob(ipft) = pdf - end do + allocate(current_neighbor%density_prob(numpft)) - if (associated(neighbors(gi)%first_neighbor)) then - neighbors(gi)%last_neighbor%next_neighbor => current_neighbor - neighbors(gi)%last_neighbor => current_neighbor - else - neighbors(gi)%first_neighbor => current_neighbor - neighbors(gi)%last_neighbor => current_neighbor - end if + do ipft = 1, numpft + call ProbabilityDensity(pdf, ipft, g2g_dist) + current_neighbor%density_prob(ipft) = pdf + end do - neighbors(gi)%neighbor_count = neighbors(gi)%neighbor_count + 1 + if (associated(neighbors(gi)%first_neighbor)) then + neighbors(gi)%last_neighbor%next_neighbor => current_neighbor + neighbors(gi)%last_neighbor => current_neighbor + else + neighbors(gi)%first_neighbor => current_neighbor + neighbors(gi)%last_neighbor => current_neighbor + end if - ! Add current grid cell index to the neighbor's list as well - allocate(another_neighbor) - another_neighbor%next_neighbor => null() + neighbors(gi)%neighbor_count = neighbors(gi)%neighbor_count + 1 - another_neighbor%gindex = gi + ! Add current grid cell index to the neighbor's list as well + allocate(another_neighbor) + another_neighbor%next_neighbor => null() - another_neighbor%gc_dist = current_neighbor%gc_dist - allocate(another_neighbor%density_prob(numpft)) - do ipft = 1, numpft - another_neighbor%density_prob(ipft) = current_neighbor%density_prob(ipft) - end do + another_neighbor%gindex = gi - if (associated(neighbors(gj)%first_neighbor)) then - neighbors(gj)%last_neighbor%next_neighbor => another_neighbor - neighbors(gj)%last_neighbor => another_neighbor - else - neighbors(gj)%first_neighbor => another_neighbor - neighbors(gj)%last_neighbor => another_neighbor - end if + another_neighbor%gc_dist = current_neighbor%gc_dist + allocate(another_neighbor%density_prob(numpft)) + do ipft = 1, numpft + another_neighbor%density_prob(ipft) = current_neighbor%density_prob(ipft) + end do - neighbors(gj)%neighbor_count = neighbors(gj)%neighbor_count + 1 + if (associated(neighbors(gj)%first_neighbor)) then + neighbors(gj)%last_neighbor%next_neighbor => another_neighbor + neighbors(gj)%last_neighbor => another_neighbor + else + neighbors(gj)%first_neighbor => another_neighbor + neighbors(gj)%last_neighbor => another_neighbor + end if - end if dist_check - end do neighbor_search - end do gc_loop + neighbors(gj)%neighbor_count = neighbors(gj)%neighbor_count + 1 - ! Loop through the list and populate the grid cell index array for each gridcell - do gi = 1,numg + end if dist_check + end do neighbor_search + end do gc_loop - ! Start at the first neighbor of each neighborhood list - current_neighbor => neighbors(gi)%first_neighbor + ! Loop through the list and populate the grid cell index array for each gridcell + do gi = 1,numg - ! Allocate an array to hold the gridcell indices in each neighborhood - allocate(neighbors(gi)%neighbor_indices(neighbors(gi)%neighbor_count)) + ! Start at the first neighbor of each neighborhood list + current_neighbor => neighbors(gi)%first_neighbor - ! Walk through the neighborhood linked list and populate the array - ni = 1 - do while (associated(current_neighbor)) - neighbors(gi)%neighbor_indices(ni) = current_neighbor%gindex - ni = ni + 1 - current_neighbor => current_neighbor%next_neighbor - end do + ! Allocate an array to hold the gridcell indices in each neighborhood + allocate(neighbors(gi)%neighbor_indices(neighbors(gi)%neighbor_count)) - if (debug .and. iam .eq. 0) then - write(fates_log(), *) 'DGCN: g, lat, lon: ', gi, gclat(gi), gclon(gi) - write(fates_log(), *) 'DGCN: g, ncount: ', gi, neighbors(gi)%neighbor_count - do i = 1,neighbors(gi)%neighbor_count - write(fates_log(), *) 'DGCN: g, gilist: ', gi, neighbors(gi)%neighbor_indices(i) - end do - end if + ! Walk through the neighborhood linked list and populate the array + ni = 1 + do while (associated(current_neighbor)) + neighbors(gi)%neighbor_indices(ni) = current_neighbor%gindex + ni = ni + 1 + current_neighbor => current_neighbor%next_neighbor + end do - end do + if (debug .and. iam .eq. 0) then + write(fates_log(), *) 'DGCN: g, lat, lon: ', gi, gclat(gi), gclon(gi) + write(fates_log(), *) 'DGCN: g, ncount: ', gi, neighbors(gi)%neighbor_count + do i = 1,neighbors(gi)%neighbor_count + write(fates_log(), *) 'DGCN: g, gilist: ', gi, neighbors(gi)%neighbor_indices(i) + end do + end if + + end do - call t_stopf('fates-seed-init-decomp') + call t_stopf('fates-seed-init-decomp') end subroutine DetermineGridCellNeighbors