From 5d892a1dbbc1ead2eadedabaa901202a13fe4aea Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 10 Dec 2024 11:37:37 -0500 Subject: [PATCH] Updated patch allocation initialization to include discussion points with ckoven, and to have a b4b preservation bypass. --- main/EDInitMod.F90 | 69 +++++++++++++++++++++----------------- main/FatesInterfaceMod.F90 | 60 +++++++++++++++++++++------------ 2 files changed, 77 insertions(+), 52 deletions(-) diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 9fc11491c1..81a53b8f52 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -481,42 +481,42 @@ subroutine set_site_properties( nsites, sites,bc_in ) use_fates_luh_if: if (hlm_use_luh .eq. itrue) then ! MAPPING OF FATES PFTs on to HLM_PFTs with land use ! add up the area associated with each FATES PFT - ! where pft_areafrac_lu is the area of land in each HLM PFT and land use type (from surface dataset) - ! hlm_pft_map is the area of that land in each FATES PFT (from param file) - - ! First check for NaNs in bc_in(s)%pft_areafrac_lu. If so, make everything bare ground. - if ( .not. (any( isnan( bc_in(s)%pft_areafrac_lu (:,:) )) .or. isnan( bc_in(s)%baregroundfrac))) then + ! where pft_areafrac_lu is the area of land in each HLM + ! PFT and land use type (from surface dataset) + ! hlm_pft_map is the area of that land in each FATES + ! PFT (from param file) + + ! First check for NaNs in bc_in(s)%pft_areafrac_lu. + ! If so, make everything bare ground. + if ( .not. (any( isnan( bc_in(s)%pft_areafrac_lu (:,:) )) .or. & + isnan( bc_in(s)%baregroundfrac))) then do i_landusetype = 1, n_landuse_cats if (.not. is_crop(i_landusetype)) then do hlm_pft = 1,size( EDPftvarcon_inst%hlm_pft_map,2) do fates_pft = 1,numpft ! loop round all fates pfts for all hlm pfts - sites(s)%area_pft(fates_pft,i_landusetype) = sites(s)%area_pft(fates_pft,i_landusetype) + & - EDPftvarcon_inst%hlm_pft_map(fates_pft,hlm_pft) * bc_in(s)%pft_areafrac_lu(hlm_pft,i_landusetype) + sites(s)%area_pft(fates_pft,i_landusetype) = & + sites(s)%area_pft(fates_pft,i_landusetype) + & + EDPftvarcon_inst%hlm_pft_map(fates_pft,hlm_pft) * & + bc_in(s)%pft_areafrac_lu(hlm_pft,i_landusetype) end do end do !hlm_pft else - ! for crops, we need to use different logic because the bc_in(s)%pft_areafrac_lu() information only exists for natural PFTs + ! for crops, we need to use different logic because + ! the bc_in(s)%pft_areafrac_lu() information only exists for natural PFTs sites(s)%area_pft(crop_lu_pft_vector(i_landusetype),i_landusetype) = 1._r8 endif end do sites(s)%area_bareground = bc_in(s)%baregroundfrac else - !if ( all( isnan( bc_in(s)%pft_areafrac_lu (:,:))) .and. isnan(bc_in(s)%baregroundfrac)) then - ! if given all NaNs, then make everything bare ground - sites(s)%area_bareground = 1._r8 - sites(s)%area_pft(:,:) = 0._r8 - write(fates_log(),*) 'Nan values for pftareafrac. dumping site info.' - call dump_site(sites(s)) - !else - ! ! if only some things are NaN but not all, then something terrible has probably happened. crash. - ! write(fates_log(),*) 'some but, not all, of the data in the PFT by LU matrix at this site is NaN.' - ! write(fates_log(),*) 'recommend checking the dataset to see what has happened.' - ! call endrun(msg=errMsg(sourcefile, __LINE__)) - !endif + + sites(s)%area_bareground = 1._r8 + sites(s)%area_pft(:,:) = 0._r8 + endif - else + else ! use_fates_luh_if + ! MAPPING OF FATES PFTs on to HLM_PFTs ! add up the area associated with each FATES PFT ! where pft_areafrac is the area of land in each HLM PFT and (from surface dataset) @@ -539,7 +539,8 @@ subroutine set_site_properties( nsites, sites,bc_in ) ! remove tiny patches to prevent numerical errors in terminate patches if (sites(s)%area_pft(ft, i_landusetype) .lt. min_nocomp_pftfrac_perlanduse & .and. sites(s)%area_pft(ft, i_landusetype) .gt. nearzero) then - if(debug) write(fates_log(),*) 'removing small numbers in site%area_pft',s,ft,i_landusetype,sites(s)%area_pft(ft, i_landusetype) + if(debug) write(fates_log(),*) 'removing small numbers in site%area_pft', & + s,ft,i_landusetype,sites(s)%area_pft(ft, i_landusetype) sites(s)%area_pft(ft, i_landusetype)=0.0_r8 endif @@ -551,17 +552,24 @@ subroutine set_site_properties( nsites, sites,bc_in ) end do end do - ! if in nocomp mode, and the number of nocomp PFTs of a given land use type is greater than the maximum number of patches - ! allowed to be allocated for that land use type, then only keep the number of PFTs correspondign to the number of patches - ! allowed on that land use type, starting with the PFTs with greatest area coverage and working down + ! if in nocomp mode, and the number of nocomp PFTs of a given + ! land use type is greater than the maximum number of patches + ! allowed to be allocated for that land use type, then only + ! keep the number of PFTs correspondign to the number of patches + ! allowed on that land use type, starting with the PFTs with + ! greatest area coverage and working down if (hlm_use_nocomp .eq. itrue) then do i_landusetype = 1, n_landuse_cats - ! count how many PFTs have areas greater than zero and compare to the number of patches allowed - if (COUNT(sites(s)%area_pft(:, i_landusetype) .gt. 0._r8) > max_nocomp_pfts_by_landuse(i_landusetype)) then + ! count how many PFTs have areas greater than zero and + ! compare to the number of patches allowed + if (COUNT(sites(s)%area_pft(:, i_landusetype) .gt. 0._r8) > & + max_nocomp_pfts_by_landuse(i_landusetype)) then ! write current vector to log file - if(debug) write(fates_log(),*) 'too many PFTs for LU type ', i_landusetype, sites(s)%area_pft(:, i_landusetype) + if(debug) write(fates_log(),*) 'too many PFTs for LU type ', & + i_landusetype, sites(s)%area_pft(:, i_landusetype) - ! start from largest area, put that PFT's area into a temp vector, and then work down to successively smaller-area PFTs, + ! start from largest area, put that PFT's area into a temp vector, + ! and then work down to successively smaller-area PFTs, ! at the end replace the original vector with the temp vector temp_vec(:) = 0._r8 do i_pftcount = 1, max_nocomp_pfts_by_landuse(i_landusetype) @@ -572,7 +580,8 @@ subroutine set_site_properties( nsites, sites,bc_in ) sites(s)%area_pft(:, i_landusetype) = temp_vec(:) ! write adjusted vector to log file - if(debug) write(fates_log(),*) 'new PFT vector for LU type', i_landusetype, sites(s)%area_pft(:, i_landusetype) + if(debug) write(fates_log(),*) 'new PFT vector for LU type', & + i_landusetype, sites(s)%area_pft(:, i_landusetype) endif end do end if diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 3650e4a179..793a7b1f5b 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -184,6 +184,8 @@ module FatesInterfaceMod public :: DetermineGridCellNeighbors logical :: debug = .false. ! for debugging this module + + logical, parameter :: preserve_b4b = .true. contains @@ -810,32 +812,46 @@ subroutine SetFatesGlobalElements1(use_fates,surf_numpft,surf_numcft,param_reade else - ! If we are using fixed biogeography or no-comp then we - ! can also apply those constraints to maxpatch_primaryland and secondary - ! and that value will match fates_maxPatchesPerSite - - if(hlm_use_luh==ifalse) then - maxpatches_by_landuse(secondaryland:n_landuse_cats) = 0 - end if - - if(hlm_use_nocomp==itrue) then - maxpatches_by_landuse(primaryland) = max(maxpatches_by_landuse(primaryland),fates_numpft) - !! REVIEWERS: NEED THIS NEXT LINE? - !maxpatches_by_landuse(secondaryland) = max(maxpatches_by_landuse(secondaryland),fates_numpft) - - ! We hold the bareground area in one of the patches - maxpatch_total = sum(maxpatches_by_landuse(:))+1 + if_preserve_b4b: if(preserve_b4b)then + if(hlm_use_nocomp==itrue) then + maxpatches_by_landuse(primaryland) = max(maxpatches_by_landuse(primaryland),fates_numpft) + maxpatch_total = sum(maxpatches_by_landuse(:))+1 + else + maxpatch_total = sum(maxpatches_by_landuse(:)) + end if else - maxpatch_total = sum(maxpatches_by_landuse(:)) - end if + + ! If we are using fixed biogeography or no-comp then we + ! can also apply those constraints to maxpatch_primaryland and secondary + ! and that value will match fates_maxPatchesPerSite + + if(hlm_use_luh==ifalse) then + maxpatches_by_landuse(secondaryland:n_landuse_cats) = 0 + end if + + if(hlm_use_nocomp==itrue) then + + ! We do not modify maxpatches_by_landuse here, and/or ensure + ! that it is large enough to accomodate all of our pfts, because + ! we will compare it to the max_nocomp_pfts_by_landuse() array to + ! make sure all indices are larger. + + ! We hold the bareground area in one of the patches + maxpatch_total = sum(maxpatches_by_landuse(:))+1 + else + maxpatch_total = sum(maxpatches_by_landuse(:)) + end if + + ! maxpatch_total does not include HLM-side bare ground (so add 1) + ! Note: Yes, there are two bare-ground patch indices needed for no-comp, where + ! index 1 on the fates side will donate all its area to index 0 on the host + ! side, and index 1 on the host side will have zero area. + + end if if_preserve_b4b - ! maxpatch_total does not include HLM-side bare ground (so add 1) - ! Note: Yes, there are two bare-ground patch indices needed for no-comp, where - ! index 1 on the fates side will donate all its area to index 0 on the host - ! side, and index 1 on the host side will have zero area. - fates_maxPatchesPerSite = maxpatch_total+1 + end if end if