From fbde7c076c8ae14c0c33088040b5e5b38a1c84b7 Mon Sep 17 00:00:00 2001 From: John Alex Date: Wed, 11 Oct 2023 09:49:31 -0600 Subject: [PATCH 01/30] 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 b3eb3f1f2034c29b2c4eee8a82472f2f876669e9 Mon Sep 17 00:00:00 2001 From: jessica needham Date: Wed, 18 Oct 2023 16:57:30 -0700 Subject: [PATCH 02/30] remove unnecessary columns from inventory init --- main/FatesInventoryInitMod.F90 | 56 ++++++++++------------------------ 1 file changed, 16 insertions(+), 40 deletions(-) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index ec099860f1..edfa97cf59 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -16,7 +16,12 @@ module FatesInventoryInitMod ! See: https://github.com/EDmodel/ED2/blob/master/ED/src/io/ed_read_ed10_20_history.f90 ! At the time of writing this ED2 is unlicensed, and only concepts were borrowed with no direct ! code copied. - !----------------------------------------------------------------------------------------------- + ! + ! + ! Update: Jessica Needham October 2023 + ! As discussed in FATES issue #1062 we decided to remove columns not used in FATES from the + ! PSS and CSS files. + !----------------------------------------------------------------------------------------------- ! CIME GLOBALS @@ -744,14 +749,6 @@ subroutine set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name ! trk (integer) LU type index (0 non-forest, 1 secondary, 2 primary ! age (years) Time since this patch was disturbed (created) ! area (fraction) Fraction of the site occupied by this patch - ! water (NA) Water content of soil (NOT USED) - ! fsc (kg/m2) Fast Soil Carbon - ! stsc (kg/m2) Structural Soil Carbon - ! stsl (kg/m2) Structural Soil Lignin - ! ssc (kg/m2) Slow Soil Carbon - ! psc (NA) Passive Soil Carbon (NOT USED) - ! msn (kg/m2) Mineralized Soil Nitrogen - ! fsn (kg/m2) Fast Soil Nitrogen ! -------------------------------------------------------------------------------------------- use FatesSizeAgeTypeIndicesMod, only: get_age_class_index @@ -772,14 +769,6 @@ subroutine set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name character(len=patchname_strlen) :: p_name ! unique string identifier of patch real(r8) :: p_age ! Patch age [years] real(r8) :: p_area ! Patch area [fraction] - real(r8) :: p_water ! Patch water (unused) - real(r8) :: p_fsc ! Patch fast soil carbon - real(r8) :: p_stsc ! Patch structural soil carbon - real(r8) :: p_stsl ! Patch structural soil lignins - real(r8) :: p_ssc ! Patch slow soil carbon - real(r8) :: p_psc ! Patch P soil carbon - real(r8) :: p_msn ! Patch mean soil nitrogen - real(r8) :: p_fsn ! Patch fast soil nitrogen integer :: icwd ! index for counting CWD pools integer :: ipft ! index for counting PFTs real(r8) :: pftfrac ! the inverse of the total number of PFTs @@ -788,9 +777,7 @@ subroutine set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name '(F5.2,2X,A4,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2)' - read(pss_file_unit,fmt=*,iostat=ios) p_time, p_name, p_trk, p_age, p_area, & - p_water,p_fsc, p_stsc, p_stsl, p_ssc, & - p_psc, p_msn, p_fsn + read(pss_file_unit,fmt=*,iostat=ios) p_time, p_name, p_trk, p_age, p_area if (ios/=0) return @@ -798,9 +785,7 @@ subroutine set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name if( debug_inv) then write(*,fmt=wr_fmt) & - p_time, p_name, p_trk, p_age, p_area, & - p_water,p_fsc, p_stsc, p_stsl, p_ssc, & - p_psc, p_msn, p_fsn + p_time, p_name, p_trk, p_age, p_area end if ! Fill in the patch's memory structures @@ -859,12 +844,8 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & ! patch (string) patch id string associated with this cohort ! index (integer) cohort index ! dbh (cm) diameter at breast height - ! height (m) height of the tree - ! pft (integer) the plant functional type index (must be consistent with param file) + ! pft (integer) the plant functional type index (must be consistent with param file) ! n (/m2) The plant number density - ! bdead (kgC/plant)The dead biomass per indiv of this cohort (NOT USED) - ! balive (kgC/plant)The live biomass per indiv of this cohort (NOT USED) - ! avgRG (cm/yr?) Average Radial Growth (NOT USED) ! -------------------------------------------------------------------------------------------- use FatesAllometryMod , only : h_allom @@ -895,12 +876,8 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & character(len=patchname_strlen) :: p_name ! The patch associated with this cohort character(len=cohortname_strlen) :: c_name ! cohort index real(r8) :: c_dbh ! diameter at breast height (cm) - real(r8) :: c_height ! tree height (m) integer :: c_pft ! plant functional type index real(r8) :: c_nplant ! plant density (/m2) - real(r8) :: c_bdead ! dead biomass (kg) - real(r8) :: c_balive ! live biomass (kg) - real(r8) :: c_avgRG ! avg radial growth (NOT USED) real(r8) :: site_spread ! initial guess of site spread ! should be quickly re-calculated integer,parameter :: rstatus = 0 ! recruit status @@ -938,13 +915,12 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & integer, parameter :: recruitstatus = 0 - read(css_file_unit,fmt=*,iostat=ios) c_time, p_name, c_name, c_dbh, c_height, & - c_pft, c_nplant, c_bdead, c_balive, c_avgRG + read(css_file_unit,fmt=*,iostat=ios) c_time, p_name, c_name, c_dbh, & + c_pft, c_nplant if( debug_inv) then write(*,fmt=wr_fmt) & - c_time, p_name, c_name, c_dbh, c_height, & - c_pft, c_nplant, c_bdead, c_balive, c_avgRG + c_time, p_name, c_name, c_dbh, c_pft, c_nplant end if if (ios/=0) return @@ -961,8 +937,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & if(.not.matched_patch)then write(fates_log(), *) 'could not match a cohort with a patch' write(fates_log(),fmt=wr_fmt) & - c_time, p_name, c_name, c_dbh, c_height, & - c_pft, c_nplant, c_bdead, c_balive, c_avgRG + c_time, p_name, c_name, c_dbh, c_pft, c_nplant call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -1215,6 +1190,7 @@ subroutine write_inventory_type1(currentSite) ! a recommended file type for restarting a run. ! The files will have a lat/long tag added to their name, and will be ! generated in the run folder. + ! JFN Oct 2023 - updated to get rid of unused ED columns ! -------------------------------------------------------------------------------- use shr_file_mod, only : shr_file_getUnit @@ -1267,8 +1243,8 @@ subroutine write_inventory_type1(currentSite) open(unit=pss_file_out,file=trim(pss_name_out), status='UNKNOWN',action='WRITE',form='FORMATTED') open(unit=css_file_out,file=trim(css_name_out), status='UNKNOWN',action='WRITE',form='FORMATTED') - write(pss_file_out,*) 'time patch trk age area water fsc stsc stsl ssc psc msn fsn' - write(css_file_out,*) 'time patch cohort dbh height pft nplant bdead alive Avgrg' + write(pss_file_out,*) 'time patch trk age area' + write(css_file_out,*) 'time patch cohort dbh pft nplant' ipatch=0 currentpatch => currentSite%youngest_patch From 10a724c2de48083607cababc290deffdfe1b91ce Mon Sep 17 00:00:00 2001 From: jessica needham Date: Wed, 18 Oct 2023 21:09:23 -0700 Subject: [PATCH 03/30] remove cohort index from inventory init files --- main/FatesInventoryInitMod.F90 | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index edfa97cf59..d1af0be58f 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -842,7 +842,6 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & ! FILE FORMAT: ! time (year) year of measurement ! patch (string) patch id string associated with this cohort - ! index (integer) cohort index ! dbh (cm) diameter at breast height ! pft (integer) the plant functional type index (must be consistent with param file) ! n (/m2) The plant number density @@ -874,7 +873,6 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & class(prt_vartypes), pointer :: prt_obj real(r8) :: c_time ! Time patch was recorded character(len=patchname_strlen) :: p_name ! The patch associated with this cohort - character(len=cohortname_strlen) :: c_name ! cohort index real(r8) :: c_dbh ! diameter at breast height (cm) integer :: c_pft ! plant functional type index real(r8) :: c_nplant ! plant density (/m2) @@ -915,12 +913,12 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & integer, parameter :: recruitstatus = 0 - read(css_file_unit,fmt=*,iostat=ios) c_time, p_name, c_name, c_dbh, & + read(css_file_unit,fmt=*,iostat=ios) c_time, p_name, c_dbh, & c_pft, c_nplant if( debug_inv) then write(*,fmt=wr_fmt) & - c_time, p_name, c_name, c_dbh, c_pft, c_nplant + c_time, p_name, c_dbh, c_pft, c_nplant end if if (ios/=0) return @@ -937,7 +935,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & if(.not.matched_patch)then write(fates_log(), *) 'could not match a cohort with a patch' write(fates_log(),fmt=wr_fmt) & - c_time, p_name, c_name, c_dbh, c_pft, c_nplant + c_time, p_name, c_dbh, c_pft, c_nplant call endrun(msg=errMsg(sourcefile, __LINE__)) end if From 97707c0275d22f95ba8a15b2a87a345eb9f5dcc0 Mon Sep 17 00:00:00 2001 From: jessica needham Date: Wed, 18 Oct 2023 21:58:10 -0700 Subject: [PATCH 04/30] add script to convert ed to fates inventory init files --- tools/ed2_to_fates_inventory_init.py | 46 ++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 tools/ed2_to_fates_inventory_init.py diff --git a/tools/ed2_to_fates_inventory_init.py b/tools/ed2_to_fates_inventory_init.py new file mode 100644 index 0000000000..a8959134ee --- /dev/null +++ b/tools/ed2_to_fates_inventory_init.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python + +### This script takes a ED2 style inventory init file and converts it to a file compatible with FATES. +# It accepts the following flags: +# --type : patch or cohort +# --fin : input filename +# --fout : output file name + +import argparse +import pandas as pd +import sys + +def main(): + parser = argparse.ArgumentParser(description='Parse command line arguments to this script.') + # + parser.add_argument('--type', dest='fatestype', type=str, help="patch or cohort. Required.", required=True) + parser.add_argument('--fin', dest='fnamein', type=str, help="Input filename. Required.", required=True) + parser.add_argument('--fout', dest='fnameout', type=str, help="Output filename. Required.", required=True) + + args = parser.parse_args() + + # open the input data + dsin = pd.read_csv(args.fnamein, delim_whitespace=True) + + # if patch file delete unnecessary patch columns + if args.fatestype == 'patch' : + keep_col = ['time', 'patch', 'trk', 'age', 'area'] + newds = dsin[keep_col] + + + # if cohort file delete unnecessary cohort columns + elif args.fatestype == 'cohort' : + keep_col = ['time', 'patch', 'dbh', 'pft', 'nplant'] + newds = dsin[keep_col] + + else : + print("type must be one of patch or cohort") + + + newds.to_csv(args.fnameout , index=False, sep=' ') +# ======================================================================================================== +# This is the actual call to main + +if __name__ == "__main__": + main() + From c541a74d4ca57492ce146ea66ddc2a46c79be1d8 Mon Sep 17 00:00:00 2001 From: jessica needham Date: Thu, 19 Oct 2023 10:18:04 -0700 Subject: [PATCH 05/30] clean up python script -make simpler --- tools/ed2_to_fates_inventory_init.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/tools/ed2_to_fates_inventory_init.py b/tools/ed2_to_fates_inventory_init.py index a8959134ee..a4a00475f4 100644 --- a/tools/ed2_to_fates_inventory_init.py +++ b/tools/ed2_to_fates_inventory_init.py @@ -2,9 +2,7 @@ ### This script takes a ED2 style inventory init file and converts it to a file compatible with FATES. # It accepts the following flags: -# --type : patch or cohort # --fin : input filename -# --fout : output file name import argparse import pandas as pd @@ -13,23 +11,28 @@ def main(): parser = argparse.ArgumentParser(description='Parse command line arguments to this script.') # - parser.add_argument('--type', dest='fatestype', type=str, help="patch or cohort. Required.", required=True) parser.add_argument('--fin', dest='fnamein', type=str, help="Input filename. Required.", required=True) - parser.add_argument('--fout', dest='fnameout', type=str, help="Output filename. Required.", required=True) - + args = parser.parse_args() + # is it a pss or css file? + filetype = args.fnamein.split('.')[1] + + # make the new file name + base_filename = args.fnamein.split('.')[0] + output_filename = f"{base_filename}_{'fates'}.{filetype}" + # open the input data dsin = pd.read_csv(args.fnamein, delim_whitespace=True) # if patch file delete unnecessary patch columns - if args.fatestype == 'patch' : + if filetype == 'pss' : keep_col = ['time', 'patch', 'trk', 'age', 'area'] newds = dsin[keep_col] # if cohort file delete unnecessary cohort columns - elif args.fatestype == 'cohort' : + elif filetype == 'css' : keep_col = ['time', 'patch', 'dbh', 'pft', 'nplant'] newds = dsin[keep_col] @@ -37,7 +40,7 @@ def main(): print("type must be one of patch or cohort") - newds.to_csv(args.fnameout , index=False, sep=' ') + newds.to_csv(output_filename, index=False, sep=' ') # ======================================================================================================== # This is the actual call to main From 1108024d67d3a8bfc9775c66694508ced91c8b49 Mon Sep 17 00:00:00 2001 From: jessica needham Date: Mon, 23 Oct 2023 13:49:38 -0700 Subject: [PATCH 06/30] add height back in --- main/FatesInventoryInitMod.F90 | 10 ++++++---- tools/ed2_to_fates_inventory_init.py | 4 ++-- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index d1af0be58f..fdef68cf11 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -843,6 +843,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & ! time (year) year of measurement ! patch (string) patch id string associated with this cohort ! dbh (cm) diameter at breast height + ! height (m) height of vegetation in m. Currently not used. ! pft (integer) the plant functional type index (must be consistent with param file) ! n (/m2) The plant number density ! -------------------------------------------------------------------------------------------- @@ -874,6 +875,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & real(r8) :: c_time ! Time patch was recorded character(len=patchname_strlen) :: p_name ! The patch associated with this cohort real(r8) :: c_dbh ! diameter at breast height (cm) + real(r8) :: c_height ! tree height (m) integer :: c_pft ! plant functional type index real(r8) :: c_nplant ! plant density (/m2) real(r8) :: site_spread ! initial guess of site spread @@ -914,11 +916,11 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & read(css_file_unit,fmt=*,iostat=ios) c_time, p_name, c_dbh, & - c_pft, c_nplant + c_height, c_pft, c_nplant if( debug_inv) then write(*,fmt=wr_fmt) & - c_time, p_name, c_dbh, c_pft, c_nplant + c_time, p_name, c_dbh, c_height, c_pft, c_nplant end if if (ios/=0) return @@ -935,7 +937,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & if(.not.matched_patch)then write(fates_log(), *) 'could not match a cohort with a patch' write(fates_log(),fmt=wr_fmt) & - c_time, p_name, c_dbh, c_pft, c_nplant + c_time, p_name, c_dbh, c_height, c_pft, c_nplant call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -1242,7 +1244,7 @@ subroutine write_inventory_type1(currentSite) open(unit=css_file_out,file=trim(css_name_out), status='UNKNOWN',action='WRITE',form='FORMATTED') write(pss_file_out,*) 'time patch trk age area' - write(css_file_out,*) 'time patch cohort dbh pft nplant' + write(css_file_out,*) 'time patch cohort dbh height pft nplant' ipatch=0 currentpatch => currentSite%youngest_patch diff --git a/tools/ed2_to_fates_inventory_init.py b/tools/ed2_to_fates_inventory_init.py index a4a00475f4..b258341a08 100644 --- a/tools/ed2_to_fates_inventory_init.py +++ b/tools/ed2_to_fates_inventory_init.py @@ -33,11 +33,11 @@ def main(): # if cohort file delete unnecessary cohort columns elif filetype == 'css' : - keep_col = ['time', 'patch', 'dbh', 'pft', 'nplant'] + keep_col = ['time', 'patch', 'dbh', 'height', 'pft', 'nplant'] newds = dsin[keep_col] else : - print("type must be one of patch or cohort") + print("file type must be one of patch (pss) or cohort (css)") newds.to_csv(output_filename, index=False, sep=' ') From 3d0e54df146ec1dd005ffcad2987ebc8ffc6e50c Mon Sep 17 00:00:00 2001 From: John Alex Date: Wed, 25 Oct 2023 09:51:11 -0600 Subject: [PATCH 07/30] 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 701a47172f70ceca27d8015039218877283c42f0 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 31 Oct 2023 12:30:47 -0400 Subject: [PATCH 08/30] Removed some ed syntax conventions, enabled using height in the inventory data --- main/FatesInventoryInitMod.F90 | 51 +++++++++++++++++++++------------- 1 file changed, 31 insertions(+), 20 deletions(-) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index fdef68cf11..92ed690786 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -11,9 +11,7 @@ module FatesInventoryInitMod ! site, or a small collection of sparse/irregularly spaced group of sites ! ! Created: Ryan Knox June 2017 - ! This code borrows heavily in concept from what is done in ED2. We will also do our best to - ! maintain compatibility with the PSS/CSS file formats that were used in ED2. - ! See: https://github.com/EDmodel/ED2/blob/master/ED/src/io/ed_read_ed10_20_history.f90 + ! This code borrows heavily in concept from what is done in ED2. ! At the time of writing this ED2 is unlicensed, and only concepts were borrowed with no direct ! code copied. ! @@ -295,7 +293,7 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) if( inv_format_list(invsite) == 1 ) then - call set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name) + call set_inventory_patch_type1(newpatch,pss_file_unit,ipa,ios,patch_name) end if ! Add it to the site's patch list @@ -384,7 +382,7 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) invcohortloop: do if ( inv_format_list(invsite) == 1 ) then - call set_inventory_edcohort_type1(sites(s),bc_in(s),css_file_unit, & + call set_inventory_cohort_type1(sites(s),bc_in(s),css_file_unit, & npatches, patch_pointer_vec,patch_name_vec, ios) end if if ( ios/=0 ) exit @@ -621,8 +619,8 @@ subroutine assess_inventory_sites(sitelist_file_unit,nsites, inv_format_list, & ! ! type integer We will accomodate different file format with different ! field values as the need arises. format 1 will read in - ! datasets via "set_inventory_edpatch_type1()", - ! "set_inventory_edcohort_type1()" + ! datasets via "set_inventory_patch_type1()", + ! "set_inventory_cohort_type1()" ! ! latitude float The geographic latitude coordinate of the site ! longitude float The geogarphic longitude coordinate of the site @@ -734,7 +732,7 @@ end subroutine assess_inventory_sites ! ============================================================================================== - subroutine set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name) + subroutine set_inventory_patch_type1(newpatch,pss_file_unit,ipa,ios,patch_name) ! -------------------------------------------------------------------------------------------- ! This subroutine reads in a line of an inventory patch file (pss) @@ -823,12 +821,12 @@ subroutine set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name end do return - end subroutine set_inventory_edpatch_type1 + end subroutine set_inventory_patch_type1 ! ============================================================================================== - subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & + subroutine set_inventory_cohort_type1(csite,bc_in,css_file_unit,npatches, & patch_pointer_vec,patch_name_vec,ios) ! -------------------------------------------------------------------------------------------- @@ -842,8 +840,8 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & ! FILE FORMAT: ! time (year) year of measurement ! patch (string) patch id string associated with this cohort - ! dbh (cm) diameter at breast height - ! height (m) height of vegetation in m. Currently not used. + ! dbh (cm) diameter at breast height. Optional, set height to negative if used + ! height (m) height of vegetation in m. Optional, set dbh to negative if used ! pft (integer) the plant functional type index (must be consistent with param file) ! n (/m2) The plant number density ! -------------------------------------------------------------------------------------------- @@ -912,6 +910,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & real(r8), parameter :: abnormal_large_nplant = 1000.0_r8 ! Used to catch bad values real(r8), parameter :: abnormal_large_dbh = 500.0_r8 ! I've never heard of a tree > 3m + real(r8), parameter :: abnormal_large_height = 500.0_r8 ! I've never heard of a tree > 500m tall integer, parameter :: recruitstatus = 0 @@ -958,9 +957,10 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if (c_dbh <=0 ) then + if (c_dbh < 0._r8 .and. c_height < 0._r8) then write(fates_log(), *) 'inventory dbh: ', c_dbh - write(fates_log(), *) 'The inventory produced a cohort with <= 0 dbh' + write(fates_log(), *) 'and inventory height: ',c_height + write(fates_log(), *) 'are both zero. One must be positive.' call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -970,6 +970,12 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & call endrun(msg=errMsg(sourcefile, __LINE__)) end if + if (c_height > abnormal_large_height ) then + write(fates_log(), *) 'inventory height: ', c_height + write(fates_log(), *) 'The inventory produced a cohort with very large height [m]' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + if (c_nplant <=0 ) then write(fates_log(), *) 'inventory nplant: ', c_nplant write(fates_log(), *) 'The inventory produced a cohort with <= 0 density /m2' @@ -1005,10 +1011,17 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & endif temp_cohort%n = c_nplant * cpatch%area / real(ncohorts_to_create,r8) - temp_cohort%dbh = c_dbh + temp_cohort%crowndamage = 1 ! assume undamaged - call h_allom(c_dbh,temp_cohort%pft,temp_cohort%height) + if( c_dbh> 0._r8)then + temp_cohort%dbh = c_dbh + call h_allom(c_dbh,temp_cohort%pft,temp_cohort%height) + else + temp_cohort%height = c_height + call h2d_allom(c_height,temp_cohort%pft,temp_cohort%dbh) + end if + temp_cohort%canopy_trim = 1.0_r8 ! Determine the phenology status and the elongation factors. @@ -1177,7 +1190,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & end do return - end subroutine set_inventory_edcohort_type1 + end subroutine set_inventory_cohort_type1 ! ==================================================================================== @@ -1185,9 +1198,7 @@ subroutine write_inventory_type1(currentSite) ! -------------------------------------------------------------------------------- ! This subroutine writes the cohort/patch inventory type files in the "type 1" - ! format. Note that for compatibility with ED2, we chose an old type that has - ! both extra unused fields and is missing fields from FATES. THis is not - ! a recommended file type for restarting a run. + ! format. ! The files will have a lat/long tag added to their name, and will be ! generated in the run folder. ! JFN Oct 2023 - updated to get rid of unused ED columns From 9c556ff03413be1241ca40e08187d751c6669321 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 31 Oct 2023 12:37:23 -0400 Subject: [PATCH 09/30] added check in inventory init to make sure at least on of dbh or height is positive --- main/FatesInventoryInitMod.F90 | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 92ed690786..962e05212f 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -957,13 +957,20 @@ subroutine set_inventory_cohort_type1(csite,bc_in,css_file_unit,npatches, & call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if (c_dbh < 0._r8 .and. c_height < 0._r8) then + if (c_dbh < nearzero .and. c_height < nearzero) then write(fates_log(), *) 'inventory dbh: ', c_dbh write(fates_log(), *) 'and inventory height: ',c_height - write(fates_log(), *) 'are both zero. One must be positive.' + write(fates_log(), *) 'are both zero or negative. One must be positive.' call endrun(msg=errMsg(sourcefile, __LINE__)) end if + if (c_dbh > nearzero .and. c_height > nearzero) then + write(fates_log(), *) 'inventory dbh: ', c_dbh + write(fates_log(), *) 'and inventory height: ',c_height + write(fates_log(), *) 'are both positive. One must be zero or negative.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + if (c_dbh > abnormal_large_dbh ) then write(fates_log(), *) 'inventory dbh: ', c_nplant write(fates_log(), *) 'The inventory produced a cohort with very large diameter [cm]' From 6dfe90a595a0521bb9f08dabe2a834b6144b2a9b Mon Sep 17 00:00:00 2001 From: jessica needham Date: Wed, 1 Nov 2023 10:18:42 -0700 Subject: [PATCH 10/30] add nearzero to inventory initialization mod --- main/FatesInventoryInitMod.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 962e05212f..fef10084f3 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -29,6 +29,7 @@ module FatesInventoryInitMod use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : pi_const use FatesConstantsMod, only : itrue + use FatesConstantsMod, only : nearzero use FatesGlobals , only : endrun => fates_endrun use FatesGlobals , only : fates_log use EDParamsMod , only : regeneration_model From 3c7837755329238d98506ee59a6909e729362d56 Mon Sep 17 00:00:00 2001 From: John Alex Date: Wed, 1 Nov 2023 12:27:20 -0600 Subject: [PATCH 11/30] 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 f3c9b54e698bb951ea1e4f87c9b3175b2db92714 Mon Sep 17 00:00:00 2001 From: Gregory Lemieux Date: Wed, 1 Nov 2023 16:32:02 -0700 Subject: [PATCH 12/30] add timesince to correct "YEAR" calculation for landuse data tool --- tools/luh2/luh2.py | 6 ++++-- tools/luh2/luh2mod.py | 3 +++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/tools/luh2/luh2.py b/tools/luh2/luh2.py index d0cd91afec..c5f3983e3e 100644 --- a/tools/luh2/luh2.py +++ b/tools/luh2/luh2.py @@ -59,10 +59,12 @@ def main(): # Add additional required variables for the host land model # Add 'YEAR' as a variable. - # This is an old requirement of the HLM and should simply be a copy of the `time` dimension # If we are merging, we might not need to do this, so check to see if its there already + # This is a requirement of the HLM dyn_subgrid module and should be the actual year. + # Note that the time variable from the LUH2 data is 'years since ...' so we need to + # add the input data year if (not "YEAR" in list(regrid_luh2.variables)): - regrid_luh2["YEAR"] = regrid_luh2.time + regrid_luh2["YEAR"] = regrid_luh2.time + regrid_luh2.timesince regrid_luh2["LONGXY"] = ds_regrid_target["LONGXY"] # TO DO: double check if this is strictly necessary regrid_luh2["LATIXY"] = ds_regrid_target["LATIXY"] # TO DO: double check if this is strictly necessary diff --git a/tools/luh2/luh2mod.py b/tools/luh2/luh2mod.py index c8534d42a9..d0f21b9060 100644 --- a/tools/luh2/luh2mod.py +++ b/tools/luh2/luh2mod.py @@ -66,6 +66,9 @@ def PrepDataset(input_dataset,start=None,stop=None,merge_flag=False): # the start/stop is out of range input_dataset = input_dataset.sel(time=slice(years_since_start,years_since_stop)) + # Save the timesince as a variable for future use + input_dataset["timesince"] = time_since + # Correct the necessary variables for both datasets # We don't need to Prep the incoming dataset if it's being opened to merge if(not merge_flag): From e7e9367176c192b5f341f7be4860431eda4faaf4 Mon Sep 17 00:00:00 2001 From: jessica needham Date: Wed, 1 Nov 2023 20:58:41 -0700 Subject: [PATCH 13/30] fix LL logic to allow case where all patches and cohorts are identical --- main/FatesInventoryInitMod.F90 | 56 ++++++++++++++++++++++++++++++++-- 1 file changed, 54 insertions(+), 2 deletions(-) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index fef10084f3..544cb3fc7f 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -168,6 +168,11 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) character(len=patchname_strlen), allocatable :: patch_name_vec(:) ! vector of patch ID strings real(r8) :: basal_area_postf ! basal area before fusion (m2/ha) real(r8) :: basal_area_pref ! basal area after fusion (m2/ha) + real(r8) :: min_patch_age + real(r8) :: max_patch_age + real(r8) :: min_cohort_dbh + real(r8) :: max_cohort_dbh + ! I. Load the inventory list file, do some file handle checks ! ------------------------------------------------------------------------------------------ @@ -367,6 +372,7 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) enddo end if + ! OPEN THE CSS FILE ! --------------------------------------------------------------------------------------- css_file_unit = shr_file_getUnit() @@ -399,8 +405,53 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) deallocate(patch_pointer_vec,patch_name_vec) - ! now that we've read in the patch and cohort info, check to see if there is any real age info - if ( abs(sites(s)%youngest_patch%age - sites(s)%oldest_patch%age) <= nearzero .and. & + + ! if all patches are identical in age and biomass then don't change the order of the LL + min_patch_age = 0._r8 + max_patch_age = 0._r8 + min_cohort_dbh = 100000._r8 + max_cohort_dbh = 0._r8 + + ! get min and max patch age and cohort dbh + currentpatch => sites(s)%youngest_patch + do while(associated(currentpatch)) + + if ( currentpatch%age > max_patch_age ) then + max_patch_age = currentpatch%age + else if ( currentpatch%age < min_patch_age ) then + min_patch_age = currentpatch%age + end if + + currentcohort => currentpatch%tallest + do while(associated(currentcohort)) + + if ( currentcohort%dbh > max_cohort_dbh ) then + max_cohort_dbh = currentcohort%dbh + else if ( currentcohort%dbh < min_cohort_dbh ) then + min_cohort_dbh = currentcohort%dbh + end if + + currentcohort => currentcohort%shorter + end do + currentPatch => currentpatch%older + enddo + + if (debug_inv) then + write(fates_log(),*) 'min patch age', min_patch_age + write(fates_log(),*) 'max patch age', max_patch_age + write(fates_log(),*) 'min cohort dbh', min_cohort_dbh + write(fates_log(),*) 'max cohort dbh', max_cohort_dbh + end if + + if ( min_patch_age .eq. max_patch_age .and. min_cohort_dbh .eq. max_cohort_dbh ) then + + if(debug_inv)then + write(fates_log(), *) 'All patches and cohorts are identical' + end if + + + ! now that we've read in the patch and cohort info, check to see if there is any real age info + else if ( abs(sites(s)%youngest_patch%age - sites(s)%oldest_patch%age) <= nearzero .and. & associated(sites(s)%youngest_patch%older) ) then ! so there are at least two patches and the oldest and youngest are the same age. @@ -501,6 +552,7 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) ! ---------------------------------------------------------------------------------------- ipa=1 total_cohorts = 0 + currentpatch => sites(s)%youngest_patch do while(associated(currentpatch)) currentpatch%patchno = ipa From 8e2e9f74e83ef2bc9082a5a4e67fee9936b98f4d Mon Sep 17 00:00:00 2001 From: jessica needham Date: Thu, 2 Nov 2023 14:42:58 -0700 Subject: [PATCH 14/30] remove sorting of ll after init data read --- main/FatesInventoryInitMod.F90 | 254 +++++++++++++++++---------------- 1 file changed, 132 insertions(+), 122 deletions(-) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 544cb3fc7f..f713102998 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -168,10 +168,12 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) character(len=patchname_strlen), allocatable :: patch_name_vec(:) ! vector of patch ID strings real(r8) :: basal_area_postf ! basal area before fusion (m2/ha) real(r8) :: basal_area_pref ! basal area after fusion (m2/ha) - real(r8) :: min_patch_age - real(r8) :: max_patch_age - real(r8) :: min_cohort_dbh - real(r8) :: max_cohort_dbh + real(r8) :: n_pref + real(r8) :: n_postf + ! real(r8) :: min_patch_age +! real(r8) :: max_patch_age +! real(r8) :: min_cohort_dbh +! real(r8) :: max_cohort_dbh ! I. Load the inventory list file, do some file handle checks @@ -406,137 +408,139 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) deallocate(patch_pointer_vec,patch_name_vec) - ! if all patches are identical in age and biomass then don't change the order of the LL - min_patch_age = 0._r8 - max_patch_age = 0._r8 - min_cohort_dbh = 100000._r8 - max_cohort_dbh = 0._r8 - - ! get min and max patch age and cohort dbh - currentpatch => sites(s)%youngest_patch - do while(associated(currentpatch)) - - if ( currentpatch%age > max_patch_age ) then - max_patch_age = currentpatch%age - else if ( currentpatch%age < min_patch_age ) then - min_patch_age = currentpatch%age - end if - - currentcohort => currentpatch%tallest - do while(associated(currentcohort)) - - if ( currentcohort%dbh > max_cohort_dbh ) then - max_cohort_dbh = currentcohort%dbh - else if ( currentcohort%dbh < min_cohort_dbh ) then - min_cohort_dbh = currentcohort%dbh - end if - - currentcohort => currentcohort%shorter - end do - currentPatch => currentpatch%older - enddo - - if (debug_inv) then - write(fates_log(),*) 'min patch age', min_patch_age - write(fates_log(),*) 'max patch age', max_patch_age - write(fates_log(),*) 'min cohort dbh', min_cohort_dbh - write(fates_log(),*) 'max cohort dbh', max_cohort_dbh - end if + ! ! if all patches are identical in age and biomass then don't change the order of the LL + ! min_patch_age = 0._r8 + ! max_patch_age = 0._r8 + ! min_cohort_dbh = 100000._r8 + ! max_cohort_dbh = 0._r8 + + ! ! get min and max patch age and cohort dbh + ! currentpatch => sites(s)%youngest_patch + ! do while(associated(currentpatch)) + + ! if ( currentpatch%age > max_patch_age ) then + ! max_patch_age = currentpatch%age + ! else if ( currentpatch%age < min_patch_age ) then + ! min_patch_age = currentpatch%age + ! end if + + ! currentcohort => currentpatch%tallest + ! do while(associated(currentcohort)) + + ! if ( currentcohort%dbh > max_cohort_dbh ) then + ! max_cohort_dbh = currentcohort%dbh + ! else if ( currentcohort%dbh < min_cohort_dbh ) then + ! min_cohort_dbh = currentcohort%dbh + ! end if + + ! currentcohort => currentcohort%shorter + ! end do + ! currentPatch => currentpatch%older + ! enddo + + ! if (debug_inv) then + ! write(fates_log(),*) 'min patch age', min_patch_age + ! write(fates_log(),*) 'max patch age', max_patch_age + ! write(fates_log(),*) 'min cohort dbh', min_cohort_dbh + ! write(fates_log(),*) 'max cohort dbh', max_cohort_dbh + ! end if - if ( min_patch_age .eq. max_patch_age .and. min_cohort_dbh .eq. max_cohort_dbh ) then + ! if ( min_patch_age .eq. max_patch_age .and. min_cohort_dbh .eq. max_cohort_dbh ) then - if(debug_inv)then - write(fates_log(), *) 'All patches and cohorts are identical' - end if + ! if(debug_inv)then + ! write(fates_log(), *) 'All patches and cohorts are identical' + ! end if - ! now that we've read in the patch and cohort info, check to see if there is any real age info - else if ( abs(sites(s)%youngest_patch%age - sites(s)%oldest_patch%age) <= nearzero .and. & - associated(sites(s)%youngest_patch%older) ) then - - ! so there are at least two patches and the oldest and youngest are the same age. - ! this means that sorting by age wasn't very useful. try sorting by total biomass instead - - ! first calculate the biomass in each patch. simplest way is to use the patch fusion criteria - currentpatch => sites(s)%youngest_patch - do while(associated(currentpatch)) - call patch_pft_size_profile(currentPatch) - currentPatch => currentpatch%older - enddo - - ! now we need to sort them. - ! first generate a new head of the linked list. - head_of_unsorted_patch_list => sites(s)%youngest_patch%older - - ! reset the site-level patch linked list, keeping only the youngest patch. - sites(s)%youngest_patch%older => null() - sites(s)%youngest_patch%younger => null() - sites(s)%oldest_patch => sites(s)%youngest_patch - - ! loop through each patch in the unsorted LL, peel it off, - ! and insert it into the new, sorted LL - do while(associated(head_of_unsorted_patch_list)) - - ! first keep track of the next patch in the old (unsorted) linked list - next_in_unsorted_patch_list => head_of_unsorted_patch_list%older - - ! check the two end-cases - - ! Youngest Patch - if(sum(head_of_unsorted_patch_list%pft_agb_profile(:,:)) <= & - sum(sites(s)%youngest_patch%pft_agb_profile(:,:)))then - head_of_unsorted_patch_list%older => sites(s)%youngest_patch - head_of_unsorted_patch_list%younger => null() - sites(s)%youngest_patch%younger => head_of_unsorted_patch_list - sites(s)%youngest_patch => head_of_unsorted_patch_list - - ! Oldest Patch - else if(sum(head_of_unsorted_patch_list%pft_agb_profile(:,:)) > & - sum(sites(s)%oldest_patch%pft_agb_profile(:,:)))then - head_of_unsorted_patch_list%older => null() - head_of_unsorted_patch_list%younger => sites(s)%oldest_patch - sites(s)%oldest_patch%older => head_of_unsorted_patch_list - sites(s)%oldest_patch => head_of_unsorted_patch_list - - ! Somewhere in the middle - else - currentpatch => sites(s)%youngest_patch - do while(associated(currentpatch)) - olderpatch => currentpatch%older - if(associated(currentpatch%older)) then - if(sum(head_of_unsorted_patch_list%pft_agb_profile(:,:)) >= & - sum(currentpatch%pft_agb_profile(:,:)) .and. & - sum(head_of_unsorted_patch_list%pft_agb_profile(:,:)) < & - sum(olderpatch%pft_agb_profile(:,:))) then - ! Set the new patches pointers - head_of_unsorted_patch_list%older => currentpatch%older - head_of_unsorted_patch_list%younger => currentpatch - ! Fix the patch's older pointer - currentpatch%older => head_of_unsorted_patch_list - ! Fix the older patch's younger pointer - olderpatch%younger => head_of_unsorted_patch_list - ! Exit the loop once head sorted to avoid later re-sort - exit - end if - end if - currentPatch => olderpatch - enddo - end if - - ! now work through to the next element in the unsorted linked list - head_of_unsorted_patch_list => next_in_unsorted_patch_list - end do - endif + ! ! now that we've read in the patch and cohort info, check to see if there is any real age info + ! else if ( abs(sites(s)%youngest_patch%age - sites(s)%oldest_patch%age) <= nearzero .and. & + ! associated(sites(s)%youngest_patch%older) ) then + + ! ! so there are at least two patches and the oldest and youngest are the same age. + ! ! this means that sorting by age wasn't very useful. try sorting by total biomass instead + + ! ! first calculate the biomass in each patch. simplest way is to use the patch fusion criteria + ! currentpatch => sites(s)%youngest_patch + ! do while(associated(currentpatch)) + ! call patch_pft_size_profile(currentPatch) + ! currentPatch => currentpatch%older + ! enddo + + ! ! now we need to sort them. + ! ! first generate a new head of the linked list. + ! head_of_unsorted_patch_list => sites(s)%youngest_patch%older + + ! ! reset the site-level patch linked list, keeping only the youngest patch. + ! sites(s)%youngest_patch%older => null() + ! sites(s)%youngest_patch%younger => null() + ! sites(s)%oldest_patch => sites(s)%youngest_patch + + ! ! loop through each patch in the unsorted LL, peel it off, + ! ! and insert it into the new, sorted LL + ! do while(associated(head_of_unsorted_patch_list)) + + ! ! first keep track of the next patch in the old (unsorted) linked list + ! next_in_unsorted_patch_list => head_of_unsorted_patch_list%older + + ! ! check the two end-cases + + ! ! Youngest Patch + ! if(sum(head_of_unsorted_patch_list%pft_agb_profile(:,:)) <= & + ! sum(sites(s)%youngest_patch%pft_agb_profile(:,:)))then + ! head_of_unsorted_patch_list%older => sites(s)%youngest_patch + ! head_of_unsorted_patch_list%younger => null() + ! sites(s)%youngest_patch%younger => head_of_unsorted_patch_list + ! sites(s)%youngest_patch => head_of_unsorted_patch_list + + ! ! Oldest Patch + ! else if(sum(head_of_unsorted_patch_list%pft_agb_profile(:,:)) > & + ! sum(sites(s)%oldest_patch%pft_agb_profile(:,:)))then + ! head_of_unsorted_patch_list%older => null() + ! head_of_unsorted_patch_list%younger => sites(s)%oldest_patch + ! sites(s)%oldest_patch%older => head_of_unsorted_patch_list + ! sites(s)%oldest_patch => head_of_unsorted_patch_list + + ! ! Somewhere in the middle + ! else + ! currentpatch => sites(s)%youngest_patch + ! do while(associated(currentpatch)) + ! olderpatch => currentpatch%older + ! if(associated(currentpatch%older)) then + ! if(sum(head_of_unsorted_patch_list%pft_agb_profile(:,:)) >= & + ! sum(currentpatch%pft_agb_profile(:,:)) .and. & + ! sum(head_of_unsorted_patch_list%pft_agb_profile(:,:)) < & + ! sum(olderpatch%pft_agb_profile(:,:))) then + ! ! Set the new patches pointers + ! head_of_unsorted_patch_list%older => currentpatch%older + ! head_of_unsorted_patch_list%younger => currentpatch + ! ! Fix the patch's older pointer + ! currentpatch%older => head_of_unsorted_patch_list + ! ! Fix the older patch's younger pointer + ! olderpatch%younger => head_of_unsorted_patch_list + ! ! Exit the loop once head sorted to avoid later re-sort + ! exit + ! end if + ! end if + ! currentPatch => olderpatch + ! enddo + ! end if + + ! ! now work through to the next element in the unsorted linked list + ! head_of_unsorted_patch_list => next_in_unsorted_patch_list + ! end do + ! endif ! Report Basal Area (as a check on if things were read in) ! ------------------------------------------------------------------------------ basal_area_pref = 0.0_r8 + n_pref = 0.0_r8 currentpatch => sites(s)%youngest_patch do while(associated(currentpatch)) currentcohort => currentpatch%tallest do while(associated(currentcohort)) basal_area_pref = basal_area_pref + & - currentcohort%n*0.25*((currentcohort%dbh/100.0_r8)**2.0_r8)*pi_const + currentcohort%n*0.25*((currentcohort%dbh/100.0_r8)**2.0_r8)*pi_const + n_pref = n_pref + currentcohort%n currentcohort => currentcohort%shorter end do currentPatch => currentpatch%older @@ -546,6 +550,7 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) write(fates_log(),*) 'Basal Area from inventory, BEFORE fusion' write(fates_log(),*) 'Lat: ',sites(s)%lat,' Lon: ',sites(s)%lon write(fates_log(),*) basal_area_pref,' [m2/ha]' + write(fates_log(),*) 'number of plants: ', n_pref write(fates_log(),*) '-------------------------------------------------------' ! Update the patch index numbers and fuse the cohorts in the patches @@ -583,12 +588,14 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) ! ---------------------------------------------------------------------------------------- !call canopy_structure(sites(s),bc_in(s)) basal_area_postf = 0.0_r8 + n_postf = 0.0_r8 currentpatch => sites(s)%youngest_patch do while(associated(currentpatch)) currentcohort => currentpatch%tallest do while(associated(currentcohort)) basal_area_postf = basal_area_postf + & currentcohort%n*0.25*((currentcohort%dbh/100.0_r8)**2.0_r8)*pi_const + n_postf = n_postf + currentcohort%n currentcohort => currentcohort%shorter end do @@ -601,6 +608,7 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) write(fates_log(),*) 'Basal Area from inventory, AFTER fusion' write(fates_log(),*) 'Lat: ',sites(s)%lat,' Lon: ',sites(s)%lon write(fates_log(),*) basal_area_postf,' [m2/ha]' + write(fates_log(),*) 'jfn n post f :', n_postf write(fates_log(),*) '-------------------------------------------------------' ! If this is flagged as true, the post-fusion inventory will be written to file @@ -1077,9 +1085,11 @@ subroutine set_inventory_cohort_type1(csite,bc_in,css_file_unit,npatches, & if( c_dbh> 0._r8)then temp_cohort%dbh = c_dbh call h_allom(c_dbh,temp_cohort%pft,temp_cohort%height) + write(fates_log(),*) 'jfn - using dbh' else temp_cohort%height = c_height call h2d_allom(c_height,temp_cohort%pft,temp_cohort%dbh) + write(fates_log(),*) 'jfn - using height' end if temp_cohort%canopy_trim = 1.0_r8 From 35a6abd748c91f56f949d4852713e8fcfab36662 Mon Sep 17 00:00:00 2001 From: jessica needham Date: Thu, 2 Nov 2023 16:02:13 -0700 Subject: [PATCH 15/30] Clean writing of inventory files. Remove unused code. --- main/FatesInventoryInitMod.F90 | 142 ++------------------------------- 1 file changed, 5 insertions(+), 137 deletions(-) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index f713102998..c462c45112 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -170,11 +170,6 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) real(r8) :: basal_area_pref ! basal area after fusion (m2/ha) real(r8) :: n_pref real(r8) :: n_postf - ! real(r8) :: min_patch_age -! real(r8) :: max_patch_age -! real(r8) :: min_cohort_dbh -! real(r8) :: max_cohort_dbh - ! I. Load the inventory list file, do some file handle checks ! ------------------------------------------------------------------------------------------ @@ -407,129 +402,6 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) deallocate(patch_pointer_vec,patch_name_vec) - - ! ! if all patches are identical in age and biomass then don't change the order of the LL - ! min_patch_age = 0._r8 - ! max_patch_age = 0._r8 - ! min_cohort_dbh = 100000._r8 - ! max_cohort_dbh = 0._r8 - - ! ! get min and max patch age and cohort dbh - ! currentpatch => sites(s)%youngest_patch - ! do while(associated(currentpatch)) - - ! if ( currentpatch%age > max_patch_age ) then - ! max_patch_age = currentpatch%age - ! else if ( currentpatch%age < min_patch_age ) then - ! min_patch_age = currentpatch%age - ! end if - - ! currentcohort => currentpatch%tallest - ! do while(associated(currentcohort)) - - ! if ( currentcohort%dbh > max_cohort_dbh ) then - ! max_cohort_dbh = currentcohort%dbh - ! else if ( currentcohort%dbh < min_cohort_dbh ) then - ! min_cohort_dbh = currentcohort%dbh - ! end if - - ! currentcohort => currentcohort%shorter - ! end do - ! currentPatch => currentpatch%older - ! enddo - - ! if (debug_inv) then - ! write(fates_log(),*) 'min patch age', min_patch_age - ! write(fates_log(),*) 'max patch age', max_patch_age - ! write(fates_log(),*) 'min cohort dbh', min_cohort_dbh - ! write(fates_log(),*) 'max cohort dbh', max_cohort_dbh - ! end if - - ! if ( min_patch_age .eq. max_patch_age .and. min_cohort_dbh .eq. max_cohort_dbh ) then - - ! if(debug_inv)then - ! write(fates_log(), *) 'All patches and cohorts are identical' - ! end if - - - ! ! now that we've read in the patch and cohort info, check to see if there is any real age info - ! else if ( abs(sites(s)%youngest_patch%age - sites(s)%oldest_patch%age) <= nearzero .and. & - ! associated(sites(s)%youngest_patch%older) ) then - - ! ! so there are at least two patches and the oldest and youngest are the same age. - ! ! this means that sorting by age wasn't very useful. try sorting by total biomass instead - - ! ! first calculate the biomass in each patch. simplest way is to use the patch fusion criteria - ! currentpatch => sites(s)%youngest_patch - ! do while(associated(currentpatch)) - ! call patch_pft_size_profile(currentPatch) - ! currentPatch => currentpatch%older - ! enddo - - ! ! now we need to sort them. - ! ! first generate a new head of the linked list. - ! head_of_unsorted_patch_list => sites(s)%youngest_patch%older - - ! ! reset the site-level patch linked list, keeping only the youngest patch. - ! sites(s)%youngest_patch%older => null() - ! sites(s)%youngest_patch%younger => null() - ! sites(s)%oldest_patch => sites(s)%youngest_patch - - ! ! loop through each patch in the unsorted LL, peel it off, - ! ! and insert it into the new, sorted LL - ! do while(associated(head_of_unsorted_patch_list)) - - ! ! first keep track of the next patch in the old (unsorted) linked list - ! next_in_unsorted_patch_list => head_of_unsorted_patch_list%older - - ! ! check the two end-cases - - ! ! Youngest Patch - ! if(sum(head_of_unsorted_patch_list%pft_agb_profile(:,:)) <= & - ! sum(sites(s)%youngest_patch%pft_agb_profile(:,:)))then - ! head_of_unsorted_patch_list%older => sites(s)%youngest_patch - ! head_of_unsorted_patch_list%younger => null() - ! sites(s)%youngest_patch%younger => head_of_unsorted_patch_list - ! sites(s)%youngest_patch => head_of_unsorted_patch_list - - ! ! Oldest Patch - ! else if(sum(head_of_unsorted_patch_list%pft_agb_profile(:,:)) > & - ! sum(sites(s)%oldest_patch%pft_agb_profile(:,:)))then - ! head_of_unsorted_patch_list%older => null() - ! head_of_unsorted_patch_list%younger => sites(s)%oldest_patch - ! sites(s)%oldest_patch%older => head_of_unsorted_patch_list - ! sites(s)%oldest_patch => head_of_unsorted_patch_list - - ! ! Somewhere in the middle - ! else - ! currentpatch => sites(s)%youngest_patch - ! do while(associated(currentpatch)) - ! olderpatch => currentpatch%older - ! if(associated(currentpatch%older)) then - ! if(sum(head_of_unsorted_patch_list%pft_agb_profile(:,:)) >= & - ! sum(currentpatch%pft_agb_profile(:,:)) .and. & - ! sum(head_of_unsorted_patch_list%pft_agb_profile(:,:)) < & - ! sum(olderpatch%pft_agb_profile(:,:))) then - ! ! Set the new patches pointers - ! head_of_unsorted_patch_list%older => currentpatch%older - ! head_of_unsorted_patch_list%younger => currentpatch - ! ! Fix the patch's older pointer - ! currentpatch%older => head_of_unsorted_patch_list - ! ! Fix the older patch's younger pointer - ! olderpatch%younger => head_of_unsorted_patch_list - ! ! Exit the loop once head sorted to avoid later re-sort - ! exit - ! end if - ! end if - ! currentPatch => olderpatch - ! enddo - ! end if - - ! ! now work through to the next element in the unsorted linked list - ! head_of_unsorted_patch_list => next_in_unsorted_patch_list - ! end do - ! endif - ! Report Basal Area (as a check on if things were read in) ! ------------------------------------------------------------------------------ basal_area_pref = 0.0_r8 @@ -1085,11 +957,9 @@ subroutine set_inventory_cohort_type1(csite,bc_in,css_file_unit,npatches, & if( c_dbh> 0._r8)then temp_cohort%dbh = c_dbh call h_allom(c_dbh,temp_cohort%pft,temp_cohort%height) - write(fates_log(),*) 'jfn - using dbh' else temp_cohort%height = c_height call h2d_allom(c_height,temp_cohort%pft,temp_cohort%dbh) - write(fates_log(),*) 'jfn - using height' end if temp_cohort%canopy_trim = 1.0_r8 @@ -1313,9 +1183,9 @@ subroutine write_inventory_type1(currentSite) ilon_sign = 'W' end if - write(pss_name_out,'(A8,I2.2,A1,I5.5,A1,A1,I3.3,A1,I5.5,A1,A4)') & + write(pss_name_out,'(A8, I2.2, A1, I5.5, A1)') & 'pss_out_',ilat_int,'.',ilat_dec,ilat_sign,'_',ilon_int,'.',ilon_dec,ilon_sign,'.txt' - write(css_name_out,'(A8,I2.2,A1,I5.5,A1,A1,I3.3,A1,I5.5,A1,A4)') & + write(css_name_out,'(A8, I2.2, A1, A1, I3.3, A1)') & 'css_out_',ilat_int,'.',ilat_dec,ilat_sign,'_',ilon_int,'.',ilon_dec,ilon_sign,'.txt' pss_file_out = shr_file_getUnit() @@ -1334,16 +1204,14 @@ subroutine write_inventory_type1(currentSite) write(patch_str,'(A7,i4.4,A)') '' - write(pss_file_out,*) '0000 ',trim(patch_str),' 2 ',currentPatch%age,currentPatch%area/AREA, & - '0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000' + write(pss_file_out,*) '0000 ',trim(patch_str),' 2 ',currentPatch%age,currentPatch%area/AREA icohort=0 currentcohort => currentpatch%tallest do while(associated(currentcohort)) icohort=icohort+1 - write(cohort_str,'(A7,i4.4,A)') '' - write(css_file_out,*) '0000 ',trim(patch_str),' ',trim(cohort_str), & - currentCohort%dbh,0.0,currentCohort%pft,currentCohort%n/currentPatch%area,0.0,0.0,0.0 + write(css_file_out,*) '0000 ',trim(patch_str), & + currentCohort%dbh,currentCohort%height,currentCohort%pft,currentCohort%n/currentPatch%area currentcohort => currentcohort%shorter end do From 00c0ab81d1cbc06e0f94987dccda1971090f8084 Mon Sep 17 00:00:00 2001 From: jessica needham Date: Thu, 2 Nov 2023 16:14:25 -0700 Subject: [PATCH 16/30] remove another cohort name from write css --- main/FatesInventoryInitMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index c462c45112..b3d305d48d 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -1195,7 +1195,7 @@ subroutine write_inventory_type1(currentSite) open(unit=css_file_out,file=trim(css_name_out), status='UNKNOWN',action='WRITE',form='FORMATTED') write(pss_file_out,*) 'time patch trk age area' - write(css_file_out,*) 'time patch cohort dbh height pft nplant' + write(css_file_out,*) 'time patch dbh height pft nplant' ipatch=0 currentpatch => currentSite%youngest_patch From 44f9a4aa1e8eec5e75ba8f831c4fbea827e7a450 Mon Sep 17 00:00:00 2001 From: jessica needham Date: Thu, 2 Nov 2023 17:11:05 -0700 Subject: [PATCH 17/30] remove debugging checks --- main/FatesInventoryInitMod.F90 | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index b3d305d48d..304c40cfb0 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -168,9 +168,7 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) character(len=patchname_strlen), allocatable :: patch_name_vec(:) ! vector of patch ID strings real(r8) :: basal_area_postf ! basal area before fusion (m2/ha) real(r8) :: basal_area_pref ! basal area after fusion (m2/ha) - real(r8) :: n_pref - real(r8) :: n_postf - + ! I. Load the inventory list file, do some file handle checks ! ------------------------------------------------------------------------------------------ @@ -405,14 +403,12 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) ! Report Basal Area (as a check on if things were read in) ! ------------------------------------------------------------------------------ basal_area_pref = 0.0_r8 - n_pref = 0.0_r8 currentpatch => sites(s)%youngest_patch do while(associated(currentpatch)) currentcohort => currentpatch%tallest do while(associated(currentcohort)) basal_area_pref = basal_area_pref + & currentcohort%n*0.25*((currentcohort%dbh/100.0_r8)**2.0_r8)*pi_const - n_pref = n_pref + currentcohort%n currentcohort => currentcohort%shorter end do currentPatch => currentpatch%older @@ -422,7 +418,6 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) write(fates_log(),*) 'Basal Area from inventory, BEFORE fusion' write(fates_log(),*) 'Lat: ',sites(s)%lat,' Lon: ',sites(s)%lon write(fates_log(),*) basal_area_pref,' [m2/ha]' - write(fates_log(),*) 'number of plants: ', n_pref write(fates_log(),*) '-------------------------------------------------------' ! Update the patch index numbers and fuse the cohorts in the patches @@ -460,14 +455,12 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) ! ---------------------------------------------------------------------------------------- !call canopy_structure(sites(s),bc_in(s)) basal_area_postf = 0.0_r8 - n_postf = 0.0_r8 currentpatch => sites(s)%youngest_patch do while(associated(currentpatch)) currentcohort => currentpatch%tallest do while(associated(currentcohort)) basal_area_postf = basal_area_postf + & currentcohort%n*0.25*((currentcohort%dbh/100.0_r8)**2.0_r8)*pi_const - n_postf = n_postf + currentcohort%n currentcohort => currentcohort%shorter end do @@ -480,7 +473,6 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) write(fates_log(),*) 'Basal Area from inventory, AFTER fusion' write(fates_log(),*) 'Lat: ',sites(s)%lat,' Lon: ',sites(s)%lon write(fates_log(),*) basal_area_postf,' [m2/ha]' - write(fates_log(),*) 'jfn n post f :', n_postf write(fates_log(),*) '-------------------------------------------------------' ! If this is flagged as true, the post-fusion inventory will be written to file From 90f50ba399e8a7998874b74f0cae81b8dd13c2aa Mon Sep 17 00:00:00 2001 From: jessica needham Date: Fri, 3 Nov 2023 14:16:00 -0700 Subject: [PATCH 18/30] turn off hmort when soil is frozen --- biogeochem/EDMortalityFunctionsMod.F90 | 14 +++++++++----- main/EDParamsMod.F90 | 4 +++- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index bf47a5cce3..029ca907dc 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -48,7 +48,7 @@ module EDMortalityFunctionsMod contains - subroutine mortality_rates( cohort_in,bc_in,btran_ft, mean_temp, & + subroutine mortality_rates( cohort_in,bc_in, btran_ft, mean_temp, & cmort,hmort,bmort, frmort,smort,asmort,dgmort ) ! ============================================================================ @@ -56,9 +56,10 @@ subroutine mortality_rates( cohort_in,bc_in,btran_ft, mean_temp, & ! background and freezing and size and age dependent senescence ! ============================================================================ - use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm - use FatesConstantsMod, only : fates_check_param_set - use DamageMainMod, only : GetDamageMortality + use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm + use FatesConstantsMod, only : fates_check_param_set + use DamageMainMod, only : GetDamageMortality + use EDParamsmod, only : tsoil_thresh_hmort type (fates_cohort_type), intent(in) :: cohort_in type (bc_in_type), intent(in) :: bc_in @@ -156,9 +157,12 @@ subroutine mortality_rates( cohort_in,bc_in,btran_ft, mean_temp, & hmort = 0.0_r8 endif else - if(btran_ft(cohort_in%pft) <= hf_sm_threshold)then + if(btran_ft(cohort_in%pft) <= hf_sm_threshold .and. & + minval(bc_in%t_soisno_sl) - tfrz .ge. tsoil_thresh_hmort )then hmort = EDPftvarcon_inst%mort_scalar_hydrfailure(cohort_in%pft) else + write(fates_log(),*) 'jfn soil frozen - no hmort' + write(fates_log(),*) 'jfn min soil temp - ', minval(bc_in%t_soisno_sl) - tfrz hmort = 0.0_r8 endif endif diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 438d387213..489e0f209a 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -95,8 +95,10 @@ module EDParamsMod integer, public :: n_uptake_mode integer, public :: p_uptake_mode + real(r8), parameter, public :: tsoil_thresh_hmort = -2.0_r8 ! Soil temperature threshold below which hydraulic failure mortality is off (non-hydro only) + integer, parameter, public :: nclmax = 2 ! Maximum number of canopy layers - + ! parameters that govern the VAI (LAI+SAI) bins used in radiative transfer code integer, parameter, public :: nlevleaf = 30 ! number of leaf+stem layers in each canopy layer From ef1ee873eda060efc544bfeea20e73f4a74ccf03 Mon Sep 17 00:00:00 2001 From: jessica needham Date: Fri, 3 Nov 2023 15:26:08 -0700 Subject: [PATCH 19/30] remove checks --- biogeochem/EDMortalityFunctionsMod.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 029ca907dc..ea2b11b06d 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -161,8 +161,6 @@ subroutine mortality_rates( cohort_in,bc_in, btran_ft, mean_temp, & minval(bc_in%t_soisno_sl) - tfrz .ge. tsoil_thresh_hmort )then hmort = EDPftvarcon_inst%mort_scalar_hydrfailure(cohort_in%pft) else - write(fates_log(),*) 'jfn soil frozen - no hmort' - write(fates_log(),*) 'jfn min soil temp - ', minval(bc_in%t_soisno_sl) - tfrz hmort = 0.0_r8 endif endif From 86eaa052060bb940ad899975cf9a53b601742f5c Mon Sep 17 00:00:00 2001 From: jessica needham Date: Sat, 4 Nov 2023 11:46:44 -0700 Subject: [PATCH 20/30] Remove checks - rename soil temp threshold --- biogeochem/EDMortalityFunctionsMod.F90 | 6 +++--- biogeophys/EDBtranMod.F90 | 3 ++- main/EDParamsMod.F90 | 4 ++-- 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index ea2b11b06d..2291e56b39 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -59,7 +59,7 @@ subroutine mortality_rates( cohort_in,bc_in, btran_ft, mean_temp, & use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm use FatesConstantsMod, only : fates_check_param_set use DamageMainMod, only : GetDamageMortality - use EDParamsmod, only : tsoil_thresh_hmort + use EDParamsmod, only : soil_tfrz_thresh type (fates_cohort_type), intent(in) :: cohort_in type (bc_in_type), intent(in) :: bc_in @@ -157,8 +157,8 @@ subroutine mortality_rates( cohort_in,bc_in, btran_ft, mean_temp, & hmort = 0.0_r8 endif else - if(btran_ft(cohort_in%pft) <= hf_sm_threshold .and. & - minval(bc_in%t_soisno_sl) - tfrz .ge. tsoil_thresh_hmort )then + if( ( btran_ft(cohort_in%pft) <= hf_sm_threshold ) .and. & + ( ( minval(bc_in%t_soisno_sl) - tfrz ) >= soil_tfrz_thresh ) ) then hmort = EDPftvarcon_inst%mort_scalar_hydrfailure(cohort_in%pft) else hmort = 0.0_r8 diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index a785493d54..c0a84ac701 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -12,6 +12,7 @@ module EDBtranMod use EDTypesMod , only : ed_site_type use FatesPatchMod, only : fates_patch_type use EDParamsMod, only : maxpft + use EDParamsMod, only : soil_tfrz_thresh use FatesCohortMod, only : fates_cohort_type use shr_kind_mod , only : r8 => shr_kind_r8 use FatesInterfaceTypesMod , only : bc_in_type, & @@ -48,7 +49,7 @@ logical function check_layer_water(h2o_liq_vol, tempk) check_layer_water = .false. if ( h2o_liq_vol .gt. 0._r8 ) then - if ( tempk .gt. tfrz-2._r8) then + if ( tempk .gt. soil_tfrz_thresh) then check_layer_water = .true. end if end if diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 489e0f209a..b3d0bfe669 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -95,12 +95,12 @@ module EDParamsMod integer, public :: n_uptake_mode integer, public :: p_uptake_mode - real(r8), parameter, public :: tsoil_thresh_hmort = -2.0_r8 ! Soil temperature threshold below which hydraulic failure mortality is off (non-hydro only) + real(r8), parameter, public :: soil_tfrz_thresh = -2.0_r8 ! Soil temperature threshold below which hydraulic failure mortality is off (non-hydro only) in degrees C integer, parameter, public :: nclmax = 2 ! Maximum number of canopy layers ! parameters that govern the VAI (LAI+SAI) bins used in radiative transfer code - integer, parameter, public :: nlevleaf = 30 ! number of leaf+stem layers in each canopy layer + integer, parameter, public :: nlevleaf = 40 ! number of leaf+stem layers in each canopy layer real(r8), public :: dinc_vai(nlevleaf) = fates_unset_r8 ! VAI bin widths array real(r8), public :: dlower_vai(nlevleaf) = fates_unset_r8 ! lower edges of VAI bins From 63c2c47698522402119602537c1a85ea0a6519f4 Mon Sep 17 00:00:00 2001 From: jessica needham Date: Sat, 4 Nov 2023 16:54:12 -0700 Subject: [PATCH 21/30] revert nlevleaf to 30 --- main/EDParamsMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index b3d0bfe669..215d4c3ee7 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -100,7 +100,7 @@ module EDParamsMod integer, parameter, public :: nclmax = 2 ! Maximum number of canopy layers ! parameters that govern the VAI (LAI+SAI) bins used in radiative transfer code - integer, parameter, public :: nlevleaf = 40 ! number of leaf+stem layers in each canopy layer + integer, parameter, public :: nlevleaf = 30 ! number of leaf+stem layers in each canopy layer real(r8), public :: dinc_vai(nlevleaf) = fates_unset_r8 ! VAI bin widths array real(r8), public :: dlower_vai(nlevleaf) = fates_unset_r8 ! lower edges of VAI bins From 03aa382726fcc55e23c0020141b0206f33d1fa75 Mon Sep 17 00:00:00 2001 From: jessica needham Date: Sat, 4 Nov 2023 19:53:58 -0700 Subject: [PATCH 22/30] change >= to > to be consistent with btran calculation --- biogeochem/EDMortalityFunctionsMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 2291e56b39..d849121c0e 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -158,7 +158,7 @@ subroutine mortality_rates( cohort_in,bc_in, btran_ft, mean_temp, & endif else if( ( btran_ft(cohort_in%pft) <= hf_sm_threshold ) .and. & - ( ( minval(bc_in%t_soisno_sl) - tfrz ) >= soil_tfrz_thresh ) ) then + ( ( minval(bc_in%t_soisno_sl) - tfrz ) > soil_tfrz_thresh ) ) then hmort = EDPftvarcon_inst%mort_scalar_hydrfailure(cohort_in%pft) else hmort = 0.0_r8 From 696365076d3029991bc9267fc753c68425949245 Mon Sep 17 00:00:00 2001 From: jessica needham Date: Sun, 5 Nov 2023 13:21:53 -0800 Subject: [PATCH 23/30] update soil freezing in btran calculation --- biogeophys/EDBtranMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index c0a84ac701..3e2401a033 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -49,7 +49,7 @@ logical function check_layer_water(h2o_liq_vol, tempk) check_layer_water = .false. if ( h2o_liq_vol .gt. 0._r8 ) then - if ( tempk .gt. soil_tfrz_thresh) then + if ( tempk .gt. soil_tfrz_thresh + tfrz) then check_layer_water = .true. end if end if From 58732e5bda2d6881b28ee06a618903bca293bbb8 Mon Sep 17 00:00:00 2001 From: Gregory Lemieux Date: Wed, 8 Nov 2023 14:26:42 -0800 Subject: [PATCH 24/30] correct the way YEAR is updated Since regridding skips 'timesince' in the input dataset (and isn't necessary to regrid), simply pass the value from the original ds_luh2 dataset --- tools/luh2/luh2.py | 2 +- tools/luh2/luh2.sh | 2 +- tools/luh2/luh2mod.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tools/luh2/luh2.py b/tools/luh2/luh2.py index c5f3983e3e..7ea0af9fa3 100644 --- a/tools/luh2/luh2.py +++ b/tools/luh2/luh2.py @@ -64,7 +64,7 @@ def main(): # Note that the time variable from the LUH2 data is 'years since ...' so we need to # add the input data year if (not "YEAR" in list(regrid_luh2.variables)): - regrid_luh2["YEAR"] = regrid_luh2.time + regrid_luh2.timesince + regrid_luh2["YEAR"] = regrid_luh2.time + ds_luh2.timesince regrid_luh2["LONGXY"] = ds_regrid_target["LONGXY"] # TO DO: double check if this is strictly necessary regrid_luh2["LATIXY"] = ds_regrid_target["LATIXY"] # TO DO: double check if this is strictly necessary diff --git a/tools/luh2/luh2.sh b/tools/luh2/luh2.sh index 3aa246907d..1d088ba7e4 100755 --- a/tools/luh2/luh2.sh +++ b/tools/luh2/luh2.sh @@ -42,7 +42,7 @@ echo -e"storage status:\n" du -h ${OUTPUT_LOC} # Regrid the luh2 transitions data using the saved regridder weights file and merge into previous regrid output -python luh2.py -b ${START} -e ${END} -l ${TRANSITIONS} -s ${STATIC} -r ${REGRID_TARGET} -w ${REGRIDDER} \ +python luh2.py -b ${START} -e ${END}-1 -l ${TRANSITIONS} -s ${STATIC} -r ${REGRID_TARGET} -w ${REGRIDDER} \ -m ${OUTPUT_LOC}/states_regrid.nc -o ${OUTPUT_LOC}/states_trans_regrid.nc echo -e"storage status:\n" du -h ${OUTPUT_LOC} diff --git a/tools/luh2/luh2mod.py b/tools/luh2/luh2mod.py index d0f21b9060..1581f949fc 100644 --- a/tools/luh2/luh2mod.py +++ b/tools/luh2/luh2mod.py @@ -202,7 +202,7 @@ def RegridLoop(ds_to_regrid, regridder): for i in range(varlen-1): # Skip time variable - if (ds_varnames[i] != "time"): + if (not "time" in ds_varnames[i]): # Only regrid variables that match the lat/lon shape. if (ds_to_regrid[ds_varnames[i]][0].shape == (ds_to_regrid.lat.shape[0], ds_to_regrid.lon.shape[0])): From 0913be784986ee06151bea7fd0f6650891defc9b Mon Sep 17 00:00:00 2001 From: Gregory Lemieux Date: Wed, 8 Nov 2023 14:36:47 -0800 Subject: [PATCH 25/30] start setting up a refactor of the luh2 shell script --- tools/luh2/luh2.sh | 2 +- tools/luh2/luh2mod.py | 15 +++++++++------ 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/tools/luh2/luh2.sh b/tools/luh2/luh2.sh index 1d088ba7e4..3aa246907d 100755 --- a/tools/luh2/luh2.sh +++ b/tools/luh2/luh2.sh @@ -42,7 +42,7 @@ echo -e"storage status:\n" du -h ${OUTPUT_LOC} # Regrid the luh2 transitions data using the saved regridder weights file and merge into previous regrid output -python luh2.py -b ${START} -e ${END}-1 -l ${TRANSITIONS} -s ${STATIC} -r ${REGRID_TARGET} -w ${REGRIDDER} \ +python luh2.py -b ${START} -e ${END} -l ${TRANSITIONS} -s ${STATIC} -r ${REGRID_TARGET} -w ${REGRIDDER} \ -m ${OUTPUT_LOC}/states_regrid.nc -o ${OUTPUT_LOC}/states_trans_regrid.nc echo -e"storage status:\n" du -h ${OUTPUT_LOC} diff --git a/tools/luh2/luh2mod.py b/tools/luh2/luh2mod.py index 1581f949fc..9f3d26333a 100644 --- a/tools/luh2/luh2mod.py +++ b/tools/luh2/luh2mod.py @@ -30,7 +30,7 @@ def PrepDataset(input_dataset,start=None,stop=None,merge_flag=False): # 'years since' style format. if(not(dstype in ('static','regrid'))): - if (dstype == 'LUH2'): + if ('LUH2' in dstype): # Get the units to determine the file time # It is expected that the units of time is 'years since ...' time_since_array = input_dataset.time.units.split() @@ -80,7 +80,7 @@ def PrepDataset(input_dataset,start=None,stop=None,merge_flag=False): def PrepDataset_ESMF(input_dataset,dsflag,dstype): if (dsflag): - if(dstype == "LUH2"): + if("LUH2" in dstype): print("PrepDataset: LUH2") input_dataset = BoundsVariableFixLUH2(input_dataset) elif(dstype == "surface"): @@ -142,15 +142,18 @@ def CheckDataset(input_dataset): dsflag = False dsvars = list(input_dataset.variables) if('primf' in dsvars or - 'primf_to_secdn' in dsvars or any('irrig' in subname for subname in dsvars)): - dstype = 'LUH2' + if ('primf_to_secdn' in dsvars): + dstype = 'LUH2_transitions' + else: + dstype = 'LUH2' + dsflag = True - print("LUH2") + # print("LUH2") elif('natpft' in dsvars): dstype = 'surface' dsflag = True - print("Surface") + # print("Surface") elif('icwtr' in dsvars): dstype = 'static' dsflag = True From 17ed1b3410113d3dd2f946c0360238016dc638bb Mon Sep 17 00:00:00 2001 From: Gregory Lemieux Date: Wed, 8 Nov 2023 14:43:11 -0800 Subject: [PATCH 26/30] fix how luh2 data tool check luh2 file type --- tools/luh2/luh2mod.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/luh2/luh2mod.py b/tools/luh2/luh2mod.py index 9f3d26333a..801baa96fc 100644 --- a/tools/luh2/luh2mod.py +++ b/tools/luh2/luh2mod.py @@ -141,7 +141,7 @@ def CheckDataset(input_dataset): dsflag = False dsvars = list(input_dataset.variables) - if('primf' in dsvars or + if(any('primf' in subname for subname in dsvars) or any('irrig' in subname for subname in dsvars)): if ('primf_to_secdn' in dsvars): dstype = 'LUH2_transitions' From 77899da0b9bc6d5fd022a0b65d8a80b4580caab2 Mon Sep 17 00:00:00 2001 From: Gregory Lemieux Date: Wed, 8 Nov 2023 19:28:44 -0800 Subject: [PATCH 27/30] adding attribute copy to luh2 data tool to avoid losing unit data during regrid --- tools/luh2/luh2.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/tools/luh2/luh2.py b/tools/luh2/luh2.py index 7ea0af9fa3..ec111c88f6 100644 --- a/tools/luh2/luh2.py +++ b/tools/luh2/luh2.py @@ -72,6 +72,12 @@ def main(): if (not 'lsmlat' in list(regrid_luh2.dims)): regrid_luh2 = regrid_luh2.rename_dims({'lat':'lsmlat','lon':'lsmlon'}) + # Reapply the coordinate attributes. This is a workaround for an xarray bug (#8047) + # Currently only need time + regrid_luh2.time.attrs = ds_luh2.time.attrs + regrid_luh2.lat.attrs = ds_luh2.lat.attrs + regrid_luh2.lon.attrs = ds_luh2.lon.attrs + # Merge existing regrided luh2 file with merge input target # TO DO: check that the grid resolution # We could do this with an append during the write phase instead of the merge From c37d878e527eff9728b927423530dc8096f1b3c2 Mon Sep 17 00:00:00 2001 From: adrifoster Date: Tue, 14 Nov 2023 08:52:57 -0700 Subject: [PATCH 28/30] 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 29/30] 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 From bbb0b54851f59f0d2e4f9ecc76087eb352024eb5 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 17 Nov 2023 11:55:54 -0500 Subject: [PATCH 30/30] Fixed argument type declaration for the number of cwd pools, identified by Noel Keen --- biogeochem/FatesLitterMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/FatesLitterMod.F90 b/biogeochem/FatesLitterMod.F90 index e16cce69db..7d55dc9aab 100644 --- a/biogeochem/FatesLitterMod.F90 +++ b/biogeochem/FatesLitterMod.F90 @@ -447,7 +447,7 @@ subroutine adjust_SF_CWD_frac(dbh,ncwd,SF_val_CWD_frac,SF_val_CWD_frac_adj) !ARGUMENTS real(r8), intent(in) :: dbh !dbh of cohort [cm] - type(integer), intent(in) :: ncwd !number of cwd pools + integer, intent(in) :: ncwd !number of cwd pools real(r8), intent(in) :: SF_val_CWD_frac(:) !fates parameter specifying the !fraction of struct + sapw going !to each CWD class