From fbde7c076c8ae14c0c33088040b5e5b38a1c84b7 Mon Sep 17 00:00:00 2001 From: John Alex Date: Wed, 11 Oct 2023 09:49:31 -0600 Subject: [PATCH 1/5] Define new fates_param_reader_type type for abstracting param I/O, and copy FatesReadParameters() over from HLM. --- main/FatesInterfaceMod.F90 | 64 ++++++++++++++++++++++++++++--- main/FatesParametersInterface.F90 | 29 ++++++++++++++ 2 files changed, 88 insertions(+), 5 deletions(-) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 857c79336a..4f569fc140 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -64,9 +64,16 @@ module FatesInterfaceMod use EDParamsMod , only : ED_val_history_ageclass_bin_edges use EDParamsMod , only : ED_val_history_height_bin_edges use EDParamsMod , only : ED_val_history_coageclass_bin_edges - use CLMFatesParamInterfaceMod , only : FatesReadParameters - use EDParamsMod , only : p_uptake_mode - use EDParamsMod , only : n_uptake_mode + use FatesParametersInterface , only : fates_param_reader_type + use FatesParametersInterface , only : fates_parameters_type + use EDParamsMod , only : FatesRegisterParams, FatesReceiveParams + use SFParamsMod , only : SpitFireRegisterParams, SpitFireReceiveParams + use PRTInitParamsFATESMod , only : PRTRegisterParams, PRTReceiveParams + use FatesSynchronizedParamsMod, only : FatesSynchronizedParamsInst + ! TODO(jpalex): remove this direct reference to HLM code. + use CLMFatesParamInterfaceMod , only : HLM_FatesReadParameters => FatesReadParameters + use EDParamsMod , only : p_uptake_mode + use EDParamsMod , only : n_uptake_mode use EDTypesMod , only : ed_site_type use FatesConstantsMod , only : prescribed_p_uptake use FatesConstantsMod , only : prescribed_n_uptake @@ -173,6 +180,8 @@ module FatesInterfaceMod public :: set_bcs public :: UpdateFatesRMeansTStep public :: InitTimeAveragingGlobals + + private :: FatesReadParameters contains @@ -726,7 +735,7 @@ end subroutine set_bcs ! =================================================================================== - subroutine SetFatesGlobalElements1(use_fates,surf_numpft,surf_numcft) + subroutine SetFatesGlobalElements1(use_fates,surf_numpft,surf_numcft,param_reader) ! -------------------------------------------------------------------------------- ! @@ -741,13 +750,21 @@ subroutine SetFatesGlobalElements1(use_fates,surf_numpft,surf_numcft) logical, intent(in) :: use_fates ! Is fates turned on? integer, intent(in) :: surf_numpft ! Number of PFTs in surface dataset integer, intent(in) :: surf_numcft ! Number of CFTs in surface dataset + ! TODO(jpalex): make non-optional once all HLMs pass it in. + class(fates_param_reader_type), optional, intent(in) :: param_reader ! HLM-provided param file reader integer :: fates_numpft ! Number of PFTs tracked in FATES if (use_fates) then ! Self explanatory, read the fates parameter file - call FatesReadParameters() + if (present(param_reader)) then + ! new, Fates-side. + call FatesReadParameters(param_reader) + else + ! old, HLM-side. + call HLM_FatesReadParameters() + end if fates_numpft = size(prt_params%wood_density,dim=1) @@ -2142,5 +2159,42 @@ subroutine SeedlingParPatch(cpatch, & return end subroutine SeedlingParPatch + !----------------------------------------------------------------------- + ! TODO(jpalex): this belongs in FatesParametersInterface.F90, but would require + ! untangling the dependencies of the *RegisterParams methods below. + subroutine FatesReadParameters(param_reader) + implicit none + + class(fates_param_reader_type), intent(in) :: param_reader ! HLM-provided param file reader + + character(len=32) :: subname = 'FatesReadParameters' + class(fates_parameters_type), allocatable :: fates_params + logical :: is_host_file + + if ( hlm_masterproc == itrue ) then + write(fates_log(), *) 'FatesParametersInterface.F90::'//trim(subname)//' :: CLM reading ED/FATES '//' parameters ' + end if + + allocate(fates_params) + call fates_params%Init() ! fates_params class, in FatesParameterInterfaceMod + call FatesRegisterParams(fates_params) !EDParamsMod, only operates on fates_params class + call SpitFireRegisterParams(fates_params) !SpitFire Mod, only operates of fates_params class + call PRTRegisterParams(fates_params) ! PRT mod, only operates on fates_params class + call FatesSynchronizedParamsInst%RegisterParams(fates_params) !Synchronized params class in Synchronized params mod, only operates on fates_params class + + is_host_file = .false. + call param_reader%Read(is_host_file, fates_params) + + is_host_file = .true. + call param_reader%Read(is_host_file, fates_params) + + call FatesReceiveParams(fates_params) + call SpitFireReceiveParams(fates_params) + call PRTReceiveParams(fates_params) + call FatesSynchronizedParamsInst%ReceiveParams(fates_params) + + call fates_params%Destroy() + deallocate(fates_params) + end subroutine FatesReadParameters end module FatesInterfaceMod diff --git a/main/FatesParametersInterface.F90 b/main/FatesParametersInterface.F90 index b19817a091..3a220e1066 100644 --- a/main/FatesParametersInterface.F90 +++ b/main/FatesParametersInterface.F90 @@ -83,6 +83,35 @@ module FatesParametersInterface end type fates_parameters_type + ! Abstract class (to be implemented by host land models) to read in + ! parameter values. + type, abstract, public :: fates_param_reader_type + contains + ! Public functions + procedure(Read_interface), public, deferred :: Read + + end type fates_param_reader_type + + abstract interface + subroutine Read_interface(this, is_host_file, fates_params ) + ! + ! !DESCRIPTION: + ! Read 'fates_params' parameters from appropriate filename given 'is_host_file'. + ! + ! USES + import :: fates_param_reader_type + import :: fates_parameters_type + ! !ARGUMENTS: + class(fates_param_reader_type) :: this + logical, intent(in) :: is_host_file + class(fates_parameters_type), intent(inout) :: fates_params + !----------------------------------------------------------------------- + + end subroutine Read_interface + + !----------------------------------------------------------------------- + end interface + contains !----------------------------------------------------------------------- From 3d0e54df146ec1dd005ffcad2987ebc8ffc6e50c Mon Sep 17 00:00:00 2001 From: John Alex Date: Wed, 25 Oct 2023 09:51:11 -0600 Subject: [PATCH 2/5] Remove is_host_file from new fates_param_reader_type API. (because in practice, all callers use the FatesParametersInterface::RegisterParameter() default of false). Note this parameter is also called sync_with_host elsewhere. --- main/FatesInterfaceMod.F90 | 7 +------ main/FatesParametersInterface.F90 | 6 +++--- 2 files changed, 4 insertions(+), 9 deletions(-) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 4f569fc140..fc5a549cd5 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -2169,7 +2169,6 @@ subroutine FatesReadParameters(param_reader) character(len=32) :: subname = 'FatesReadParameters' class(fates_parameters_type), allocatable :: fates_params - logical :: is_host_file if ( hlm_masterproc == itrue ) then write(fates_log(), *) 'FatesParametersInterface.F90::'//trim(subname)//' :: CLM reading ED/FATES '//' parameters ' @@ -2182,11 +2181,7 @@ subroutine FatesReadParameters(param_reader) call PRTRegisterParams(fates_params) ! PRT mod, only operates on fates_params class call FatesSynchronizedParamsInst%RegisterParams(fates_params) !Synchronized params class in Synchronized params mod, only operates on fates_params class - is_host_file = .false. - call param_reader%Read(is_host_file, fates_params) - - is_host_file = .true. - call param_reader%Read(is_host_file, fates_params) + call param_reader%Read(fates_params) call FatesReceiveParams(fates_params) call SpitFireReceiveParams(fates_params) diff --git a/main/FatesParametersInterface.F90 b/main/FatesParametersInterface.F90 index 3a220e1066..c559ec4cb4 100644 --- a/main/FatesParametersInterface.F90 +++ b/main/FatesParametersInterface.F90 @@ -93,17 +93,17 @@ module FatesParametersInterface end type fates_param_reader_type abstract interface - subroutine Read_interface(this, is_host_file, fates_params ) + subroutine Read_interface(this, fates_params ) ! ! !DESCRIPTION: - ! Read 'fates_params' parameters from appropriate filename given 'is_host_file'. + ! Read 'fates_params' parameters from (HLM-provided) storage. Note this ignores + ! the legacy parameter_type.sync_with_host setting. ! ! USES import :: fates_param_reader_type import :: fates_parameters_type ! !ARGUMENTS: class(fates_param_reader_type) :: this - logical, intent(in) :: is_host_file class(fates_parameters_type), intent(inout) :: fates_params !----------------------------------------------------------------------- From 3c7837755329238d98506ee59a6909e729362d56 Mon Sep 17 00:00:00 2001 From: John Alex Date: Wed, 1 Nov 2023 12:27:20 -0600 Subject: [PATCH 3/5] Remove the old code path altogether; CESM will be updated at the same time, and E3SM will not reference this new FATES version until it's also updated. --- main/FatesInterfaceMod.F90 | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index fc5a549cd5..a158340ca5 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -70,8 +70,6 @@ module FatesInterfaceMod use SFParamsMod , only : SpitFireRegisterParams, SpitFireReceiveParams use PRTInitParamsFATESMod , only : PRTRegisterParams, PRTReceiveParams use FatesSynchronizedParamsMod, only : FatesSynchronizedParamsInst - ! TODO(jpalex): remove this direct reference to HLM code. - use CLMFatesParamInterfaceMod , only : HLM_FatesReadParameters => FatesReadParameters use EDParamsMod , only : p_uptake_mode use EDParamsMod , only : n_uptake_mode use EDTypesMod , only : ed_site_type @@ -750,21 +748,14 @@ subroutine SetFatesGlobalElements1(use_fates,surf_numpft,surf_numcft,param_reade logical, intent(in) :: use_fates ! Is fates turned on? integer, intent(in) :: surf_numpft ! Number of PFTs in surface dataset integer, intent(in) :: surf_numcft ! Number of CFTs in surface dataset - ! TODO(jpalex): make non-optional once all HLMs pass it in. - class(fates_param_reader_type), optional, intent(in) :: param_reader ! HLM-provided param file reader + class(fates_param_reader_type), intent(in) :: param_reader ! HLM-provided param file reader integer :: fates_numpft ! Number of PFTs tracked in FATES if (use_fates) then ! Self explanatory, read the fates parameter file - if (present(param_reader)) then - ! new, Fates-side. - call FatesReadParameters(param_reader) - else - ! old, HLM-side. - call HLM_FatesReadParameters() - end if + call FatesReadParameters(param_reader) fates_numpft = size(prt_params%wood_density,dim=1) From c37d878e527eff9728b927423530dc8096f1b3c2 Mon Sep 17 00:00:00 2001 From: adrifoster Date: Tue, 14 Nov 2023 08:52:57 -0700 Subject: [PATCH 4/5] fixing indents --- main/FatesInterfaceMod.F90 | 52 +++++++++++++++++++------------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 5f8ff245ac..2fa1ad8c39 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -2391,38 +2391,38 @@ end subroutine DetermineGridCellNeighbors ! ====================================================================================== - !----------------------------------------------------------------------- - ! TODO(jpalex): this belongs in FatesParametersInterface.F90, but would require - ! untangling the dependencies of the *RegisterParams methods below. - subroutine FatesReadParameters(param_reader) - implicit none - - class(fates_param_reader_type), intent(in) :: param_reader ! HLM-provided param file reader +!----------------------------------------------------------------------- +! TODO(jpalex): this belongs in FatesParametersInterface.F90, but would require +! untangling the dependencies of the *RegisterParams methods below. +subroutine FatesReadParameters(param_reader) + implicit none + + class(fates_param_reader_type), intent(in) :: param_reader ! HLM-provided param file reader - character(len=32) :: subname = 'FatesReadParameters' - class(fates_parameters_type), allocatable :: fates_params + character(len=32) :: subname = 'FatesReadParameters' + class(fates_parameters_type), allocatable :: fates_params - if ( hlm_masterproc == itrue ) then - write(fates_log(), *) 'FatesParametersInterface.F90::'//trim(subname)//' :: CLM reading ED/FATES '//' parameters ' - end if + if ( hlm_masterproc == itrue ) then + write(fates_log(), *) 'FatesParametersInterface.F90::'//trim(subname)//' :: CLM reading ED/FATES '//' parameters ' + end if - allocate(fates_params) - call fates_params%Init() ! fates_params class, in FatesParameterInterfaceMod - call FatesRegisterParams(fates_params) !EDParamsMod, only operates on fates_params class - call SpitFireRegisterParams(fates_params) !SpitFire Mod, only operates of fates_params class - call PRTRegisterParams(fates_params) ! PRT mod, only operates on fates_params class - call FatesSynchronizedParamsInst%RegisterParams(fates_params) !Synchronized params class in Synchronized params mod, only operates on fates_params class + allocate(fates_params) + call fates_params%Init() ! fates_params class, in FatesParameterInterfaceMod + call FatesRegisterParams(fates_params) !EDParamsMod, only operates on fates_params class + call SpitFireRegisterParams(fates_params) !SpitFire Mod, only operates of fates_params class + call PRTRegisterParams(fates_params) ! PRT mod, only operates on fates_params class + call FatesSynchronizedParamsInst%RegisterParams(fates_params) !Synchronized params class in Synchronized params mod, only operates on fates_params class - call param_reader%Read(fates_params) + call param_reader%Read(fates_params) - call FatesReceiveParams(fates_params) - call SpitFireReceiveParams(fates_params) - call PRTReceiveParams(fates_params) - call FatesSynchronizedParamsInst%ReceiveParams(fates_params) + call FatesReceiveParams(fates_params) + call SpitFireReceiveParams(fates_params) + call PRTReceiveParams(fates_params) + call FatesSynchronizedParamsInst%ReceiveParams(fates_params) - call fates_params%Destroy() - deallocate(fates_params) + call fates_params%Destroy() + deallocate(fates_params) end subroutine FatesReadParameters - + end module FatesInterfaceMod From 03446604fa78e80fdbfe22979212ac92a2e5bc86 Mon Sep 17 00:00:00 2001 From: Adrianna Foster Date: Thu, 16 Nov 2023 13:16:48 -0700 Subject: [PATCH 5/5] update diffs --- main/FatesInterfaceMod.F90 | 300 ++++++++++++++++++------------------- 1 file changed, 150 insertions(+), 150 deletions(-) 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