diff --git a/README.md b/README.md index 3357943468..05d6bd710e 100644 --- a/README.md +++ b/README.md @@ -33,5 +33,8 @@ https://github.com/E3SM-Project/E3SM https://github.com/ESCOMP/cesm -The FATES, E3SM and CTSM teams maintain compatability of the NGEET/FATES master branch with the E3SM master and CTSM master branches respectively. There may be some modest lag time in which the latest commit on the FATES master branch is available to these host land models (HLM) by default. This is typically correlated with FATES development updates forcing necessary changes to the FATES API. See the table of [FATES API/HLM compatibility](https://fates-users-guide.readthedocs.io/en/latest/user/release-tags-compat-table.html) for information on which fates tag corresponds to which HLM tag or commit. +The FATES, E3SM and CTSM teams maintain compatability of the NGEET/FATES master branch with the E3SM master and CTSM master branches respectively. There may be some modest lag time in which the latest commit on the FATES master branch is available to these host land models (HLM) by default. This is typically correlated with FATES development updates forcing necessary changes to the FATES API. See the table of [FATES API/HLM compatibility](https://fates-users-guide.readthedocs.io/en/latest/user/release-tags-compat-table.html) for information on which fates tag corresponds to which HLM tag or commit. + +## TODO take out this section, INFO from 8/7-2026 +This branch represents a rejuvenation of a branch that had been merged but needed to be reverted because regression tests and runs from old finidats were failing. This needs to be reinvestigated and solved later / after the summer diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 8294953f59..43d0b4ad04 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -85,18 +85,6 @@ module EDCanopyStructureMod real(r8), parameter :: co_area_target_precision = 1.0E-9_r8 - ! Relative precision target used alongside the absolute targets above. - ! Area conservation checks compare differences of two large, similarly - ! sized areas (e.g. cohort crown areas, layer areas, patch areas). For - ! such differences a fixed absolute threshold becomes meaningless once - ! the operands are large, because a single rounding of an operand of - ! magnitude x already produces a residual of order ulp(x) ~ 2.2e-16*x. - ! To avoid spuriously tripping error checks on rounding noise, the - ! comparisons below use max(absolute_target, rel_area_precision*|operand|). - ! rel_area_precision is set to relative machine precision for r8 so that - ! a difference below machine precision is never treated as a real error. - real(r8), parameter :: rel_area_precision = 1.0E-15_r8 - integer, parameter :: demotion_phase = 1 integer, parameter :: promotion_phase = 2 @@ -577,7 +565,7 @@ subroutine PromoteOrDemote(site,patch,target_layer,phase,target_area) sumpd_area = 0._r8 ic = 1 do while( ic<=n_layer .and. (promdem_area-sumpd_area) > & - max(co_area_target_precision, rel_area_precision*promdem_area)) + max(co_area_target_precision, rsnbl_math_prec*promdem_area)) cohort => layer_co(ic)%p @@ -625,7 +613,7 @@ subroutine PromoteOrDemote(site,patch,target_layer,phase,target_area) whole_or_part: if( ((layer_co(ic)%pd_area - cohort%c_area) > & - max(co_area_target_precision, rel_area_precision*cohort%c_area) ) .or. & + max(co_area_target_precision, rsnbl_math_prec*cohort%c_area) ) .or. & (layer_co(ic)%pd_area < 0._r8) ) then write(fates_log(),*) 'negative,or more area than the cohort has is being promoted/demoted' write(fates_log(),*) 'change: ',layer_co(ic)%pd_area @@ -635,7 +623,7 @@ subroutine PromoteOrDemote(site,patch,target_layer,phase,target_area) elseif ( abs(layer_co(ic)%pd_area - cohort%c_area) < & - max(co_area_target_precision, rel_area_precision*cohort%c_area) ) then + max(co_area_target_precision, rsnbl_math_prec*cohort%c_area) ) then ! Whole cohort promotion/demotion cohort%canopy_layer = cohort%canopy_layer + ilyr_change @@ -894,7 +882,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) end if if (currentPatch%total_canopy_area - (1._r8-imperfect_fraction)*currentPatch%area > & - max(area_error_1, rel_area_precision*currentPatch%area)) then + max(area_error_1, rsnbl_math_prec*currentPatch%area)) then write(fates_log(),*) 'too much canopy in summary', s, & currentPatch%nocomp_pft_label, currentPatch%total_canopy_area - (1._r8-imperfect_fraction)*currentPatch%area call endrun(msg=errMsg(sourcefile, __LINE__)) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 745eae67d6..f8436b4d4f 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -34,6 +34,7 @@ Module EDCohortDynamicsMod use EDParamsMod , only : max_cohort_per_patch use EDTypesMod , only : min_npm2, min_nppatch use EDTypesMod , only : min_n_safemath + use EDTypesMod , only : min_patch_area use EDParamsMod , only : nlevleaf use FatesConstantsMod , only : ican_upper use EDTypesMod , only : elem_diag_type @@ -86,6 +87,7 @@ Module EDCohortDynamicsMod use DamageMainMod, only : undamaged_class use FatesConstantsMod, only : i_term_mort_type_cstarv use FatesConstantsMod, only : i_term_mort_type_canlev + use FatesConstantsMod, only : i_term_mort_type_nocomp use FatesConstantsMod, only : i_term_mort_type_numdens use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) @@ -287,7 +289,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ ! terminates all cohorts when they get too small ! ! !USES: - + use FatesInterfaceTypesMod , only : hlm_use_nocomp ! ! !ARGUMENTS type (ed_site_type) , intent(inout) :: currentSite @@ -336,7 +338,27 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ struct_c = currentCohort%prt%GetState(struct_organ, carbon12_element) repro_c = currentCohort%prt%GetState(repro_organ, carbon12_element) - ! Check if number density is so low is breaks math (level 1) + ! ---------------------------------------------------------------------- + ! Numerically dangerous states. These are checked at all levels and for + ! recruits as well, so that degenerate cohorts are removed before cohort + ! fusion (which runs between the level-1 and level-2 termination calls) + ! can smear their pathology into otherwise healthy cohorts. + ! ---------------------------------------------------------------------- + + ! Number density is so low it breaks math + if (currentCohort%n < min_n_safemath .and. level == 1) then + terminate = itrue + termination_type = i_term_mort_type_numdens + if ( debug ) then + write(fates_log(),*) 'terminating cohorts 0',currentCohort%n/currentPatch%area, & + currentCohort%dbh,currentCohort%pft,call_index + endif + endif + + + + ! The rest of these are biological starvation checks, only allowed if we + ! are not dealing with a recruit (level 2) if (currentcohort%n < min_n_safemath .and. level == 1) then terminate = itrue termination_type = i_term_mort_type_numdens @@ -410,6 +432,14 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ write(fates_log(),*) 'terminating cohorts 2', currentCohort%canopy_layer,currentCohort%pft,call_index endif endif + if (hlm_use_nocomp == itrue .and. level == 4 .and. & + currentPatch%nocomp_pft_label .ne. currentCohort%pft) then + terminate = itrue + termination_type = i_term_mort_type_nocomp + if ( debug ) then + write(fates_log(),*) 'terminating cohorts 7', currentPatch%nocomp_pft_label,currentCohort%pft,call_index + endif + endif if (terminate == itrue) then call terminate_cohort(currentSite, currentPatch, currentCohort, bc_in, termination_type) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 0c7642b94a..a916caca44 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -43,6 +43,8 @@ module EDPatchDynamicsMod use EDTypesMod , only : AREA_INV use EDTypesMod , only : dump_site use FatesConstantsMod , only : rsnbl_math_prec + use FatesConstantsMod , only : rel_patch_area_floor + use FatesConstantsMod , only : min_nocomp_pftfrac_perlanduse use FatesConstantsMod , only : fates_tiny use FatesConstantsMod , only : nocomp_bareground use FatesInterfaceTypesMod , only : hlm_use_planthydro @@ -558,6 +560,9 @@ subroutine spawn_patches( currentSite, bc_in ) integer :: which_pft_allowed logical :: buffer_patch_used logical :: clear_all + real(r8) :: total_realloc_area ! total area being reallocated for this land-use label [m2] + real(r8) :: area_tol ! relative area tolerance for this land-use label [m2] + integer :: dominant_pft ! pft with the largest area fraction on this land-use label !--------------------------------------------------------------------- if (hlm_use_nocomp .eq. itrue) then @@ -1445,13 +1450,28 @@ subroutine spawn_patches( currentSite, bc_in ) currentPatch => currentPatch%younger end do + ! Total area being reallocated across this land-use label, and the relative-area + ! tolerance derived from it. Using a tolerance that scales with the area being + ! reallocated (rather than an absolute rsnbl_math_prec) prevents rounding noise on + ! differences of large areas from carving off sub-threshold sliver patches. + total_realloc_area = sum(nocomp_pft_area_vector(:)) + area_tol = max(rsnbl_math_prec, total_realloc_area * rel_patch_area_floor) + ! figure out how may PFTs on each land use type. if only 1, then the next calculation is much simpler: we just need to know which PFT is allowed. + ! Admission is gated on min_nocomp_pftfrac_perlanduse (the same floor applied at + ! initialization in EDInitMod) rather than nearzero, so the two modules agree on + ! which PFTs have meaningful area on a given land-use type. n_pfts_by_landuse = 0 + dominant_pft = 1 do i_pft = 1,numpft - if ( currentSite%area_pft(i_pft,i_land_use_label) .gt. nearzero) then + if ( currentSite%area_pft(i_pft,i_land_use_label) .ge. min_nocomp_pftfrac_perlanduse) then n_pfts_by_landuse = n_pfts_by_landuse + 1 which_pft_allowed = i_pft end if + if ( currentSite%area_pft(i_pft,i_land_use_label) .gt. & + currentSite%area_pft(dominant_pft,i_land_use_label) ) then + dominant_pft = i_pft + end if end do if ( n_pfts_by_landuse .ne. 1) then which_pft_allowed = fates_unset_int @@ -1494,7 +1514,7 @@ subroutine spawn_patches( currentSite, bc_in ) newp_area = currentPatch%area - area_to_keep fraction_to_keep = area_to_keep / currentPatch%area - if (fraction_to_keep .le. nearzero .or. area_to_keep .lt. rsnbl_math_prec) then + if (fraction_to_keep .le. nearzero .or. area_to_keep .lt. area_tol) then ! we don't want any patch area with this PFT identity at all anymore. Fuse it into the buffer patch. currentPatch%nocomp_pft_label = 0 @@ -1509,7 +1529,7 @@ subroutine spawn_patches( currentSite, bc_in ) buffer_patch_used = .true. - elseif ( area_to_keep .ge. rsnbl_math_prec .and. newp_area .ge. rsnbl_math_prec) then + elseif ( area_to_keep .ge. area_tol .and. newp_area .ge. area_tol) then ! we have more patch are of this PFT than we want, but we do want to keep some of it. ! we want to split the patch into two here. leave one patch as-is, and put the rest into the buffer patch. @@ -1541,8 +1561,15 @@ subroutine spawn_patches( currentSite, bc_in ) end do buffer_patch_used_if: if ( buffer_patch_used ) then + ! The buffer patch has now absorbed all of the area carved off the changed + ! patches and is about to be redistributed into patches with explicit, valid + ! nocomp PFT labels for this land use type. Clear the land-use-change flag it + ! inherited (via fuse_2_patches) from its donors so the remapped output patches + ! are not spuriously reprocessed on a later timestep. + buffer_patch%changed_landuse_this_ts = .false. + ! at this point, lets check that the total patch area remaining to be relabelled equals what we think that it is. - if (abs(sum(nocomp_pft_area_vector(:) - nocomp_pft_area_vector_filled(:)) - buffer_patch%area) .gt. rsnbl_math_prec) then + if (abs(sum(nocomp_pft_area_vector(:) - nocomp_pft_area_vector_filled(:)) - buffer_patch%area) .gt. area_tol) then write(fates_log(),*) 'midway through patch reallocation and things are already not adding up.', i_land_use_label write(fates_log(),*) currentSite%area_pft(:,i_land_use_label) write(fates_log(),*) '-----' @@ -1570,10 +1597,10 @@ subroutine spawn_patches( currentSite, bc_in ) max_val = maxval(newp_area_buffer_frac) ! If the max value is the only value in the array then loop through the array to find the max value pft index and insert buffer - if (abs(sum(newp_area_buffer_frac(:)) - max_val) .le. nearzero) then + if (abs(sum(newp_area_buffer_frac(:)) - max_val) .le. rel_patch_area_floor) then i_pft = 1 do while(.not. buffer_patch_in_linked_list) - if (abs(newp_area_buffer_frac(i_pft) - max_val) .le. nearzero) then + if (abs(newp_area_buffer_frac(i_pft) - max_val) .le. rel_patch_area_floor) then ! give the buffer patch the intended nocomp PFT label buffer_patch%nocomp_pft_label = i_pft @@ -1596,8 +1623,8 @@ subroutine spawn_patches( currentSite, bc_in ) nocomp_pft_loop_2: do i_pft = 1, numpft ! Check the area fraction to makes sure that this pft should have area. Also make sure that the buffer patch hasn't been - ! added to the linked list already - if ( currentSite%area_pft(i_pft,i_land_use_label) .gt. nearzero .and. .not. buffer_patch_in_linked_list) then + ! added to the linked list already. Admission uses min_nocomp_pftfrac_perlanduse (consistent with EDInitMod) rather than nearzero. + if ( currentSite%area_pft(i_pft,i_land_use_label) .ge. min_nocomp_pftfrac_perlanduse .and. .not. buffer_patch_in_linked_list) then ! Slightly complicated way of making sure that the same pfts are subtracted from each other which may help to avoid precision ! errors due to differencing between very large and very small areas @@ -1610,10 +1637,10 @@ subroutine spawn_patches( currentSite, bc_in ) area_to_keep = buffer_patch%area - newp_area fraction_to_keep = area_to_keep / buffer_patch%area - ! only bother doing this if the new new patch area needed is greater than some tiny amount - if ( newp_area .gt. rsnbl_math_prec * 0.01_r8) then + ! only bother carving off a new patch if the area needed exceeds the relative-area tolerance + if ( newp_area .gt. area_tol) then - if (area_to_keep .gt. rsnbl_math_prec) then + if (area_to_keep .gt. area_tol) then ! split buffer patch in two, keeping the smaller buffer patch to put into new patches allocate(temp_patch) @@ -1659,12 +1686,16 @@ subroutine spawn_patches( currentSite, bc_in ) call endrun(msg=errMsg(sourcefile, __LINE__)) endif else - write(fates_log(),*) 'Buffer patch still has area and it wasnt put into the linked list' - write(fates_log(),*) 'buffer_patch%area', buffer_patch%area - write(fates_log(),*) sum(nocomp_pft_area_vector_filled(:)), sum(nocomp_pft_area_vector(:)) - write(fates_log(),*) sum(nocomp_pft_area_vector_filled(:) - nocomp_pft_area_vector(:)) - - call endrun(msg=errMsg(sourcefile, __LINE__)) + ! The buffer still holds a (typically sub-tolerance) remnant of area that was + ! not carved into its own patch above, because the relative-area gates declined + ! to spawn a sliver. Rather than crash, place the entire remnant into the + ! dominant (largest-area) PFT for this land-use label. This guarantees area + ! conservation and prevents an unplaced buffer or an orphan sliver patch. + buffer_patch%nocomp_pft_label = dominant_pft + nocomp_pft_area_vector_filled(dominant_pft) = & + nocomp_pft_area_vector_filled(dominant_pft) + buffer_patch%area + call InsertPatch(currentSite, buffer_patch) + buffer_patch_in_linked_list = .true. end if end if else @@ -1678,7 +1709,7 @@ subroutine spawn_patches( currentSite, bc_in ) end if buffer_patch_used_if ! check that the area we have added is the same as the area we have taken away. if not, crash. - if ( abs(sum(nocomp_pft_area_vector_filled(:) - nocomp_pft_area_vector(:))) .gt. rsnbl_math_prec) then + if ( abs(sum(nocomp_pft_area_vector_filled(:) - nocomp_pft_area_vector(:))) .gt. area_tol) then write(fates_log(),*) 'patch reallocation logic doesnt add up. difference is: ', sum(nocomp_pft_area_vector_filled(:) - nocomp_pft_area_vector(:)) write(fates_log(),*) nocomp_pft_area_vector_filled write(fates_log(),*) nocomp_pft_area_vector @@ -1759,6 +1790,11 @@ subroutine split_patch(currentSite, currentPatch, new_patch, fraction_to_keep, a num_swb, numpft, currentSite%nlevsoil, hlm_current_tod, & hlm_regeneration_model) + ! Conserve the land-use-change flag: a piece split off a patch that still needs its + ! nocomp PFT identity remapped to the new land use type must inherit the flag, so it + ! is not left behind as an orphan patch. + new_patch%changed_landuse_this_ts = currentPatch%changed_landuse_this_ts + ! Initialize the litter pools to zero, these ! pools will be populated shortly do el=1,num_elements @@ -2933,6 +2969,9 @@ subroutine fuse_patches( csite, bc_in ) ! !LOCAL VARIABLES: type(ed_site_type) , pointer :: currentSite type(fates_patch_type), pointer :: currentPatch,tpp,tmpptr + type(fates_patch_type), pointer :: largest_patch ! disposal target for an orphan patch + type(fates_patch_type), pointer :: smallest_patch ! smallest (blocking/orphan) patch to dispose of + character(len=1024) :: warn_msg ! for defining a warning message integer :: ft,z !counters for pft and height class real(r8) :: norm !normalized difference between biomass profiles real(r8) :: profiletol !tolerance of patch fusion routine. Starts off high and is reduced if there are too many patches. @@ -3190,21 +3229,109 @@ subroutine fuse_patches( csite, bc_in ) ! against relative biomass size, it shoudnt ever get above 2 (which would mean fusing ! a zero with a nonzero biomass in a given category) if (profiletol .gt. 100._r8) then - write(fates_log(),*) 'profile tolerance is too big, this shouldnt happen.' - write(fates_log(),*) 'probably this means there are too many distinct categorical ' - write(fates_log(),*) 'patch types for the maximum number of patches' - call dump_site(currentSite) - write(fates_log(),*) 'currentSite%area_bareground', currentSite%area_bareground - do i = 1, n_landuse_cats - write(fates_log(),*) 'i, currentSite%area_pft(:,i)',i, currentSite%area_pft(:,i) + + ! We have relaxed the patch-fusion tolerance all the way and still cannot get + ! nopatches(i_lulabel) down to maxpatches_by_landuse(i_lulabel). This happens when + ! there is an irreducible "orphan" patch in this land-use label whose + ! nocomp_pft_label does not match any other patch in the same land-use label, so + ! the categorical fusion above can never merge it. Rather than crash, gracefully + ! dispose of the smallest blocking patch by relabelling its nocomp_pft_label to + ! that of the largest patch in the same land-use label and fusing it into that + ! patch. This branch only executes in the previously-fatal case, so healthy runs + ! are unaffected. + + ! find the largest patch in this land-use label (the disposal target) + largest_patch => null() + + ! Prefer a non-bareground recipient when nocomp is enabled; fall back to + ! bareground only if no other target exists. + tmpptr => currentSite%youngest_patch + do while(associated(tmpptr)) + if (tmpptr%land_use_label .eq. i_lulabel .and. & + (hlm_use_nocomp .eq. ifalse .or. tmpptr%nocomp_pft_label .ne. nocomp_bareground)) then + if (.not.associated(largest_patch)) then + largest_patch => tmpptr + else if (tmpptr%area .gt. largest_patch%area) then + largest_patch => tmpptr + endif + endif + tmpptr => tmpptr%older end do + + if (.not.associated(largest_patch)) then + tmpptr => currentSite%youngest_patch + do while(associated(tmpptr)) + if (tmpptr%land_use_label .eq. i_lulabel) then + if (.not.associated(largest_patch)) then + largest_patch => tmpptr + else if (tmpptr%area .gt. largest_patch%area) then + largest_patch => tmpptr + endif + endif + tmpptr => tmpptr%older + end do + endif + + ! find the smallest patch in this land-use label that is not the disposal + ! target itself - this is the blocking/orphan patch to remove + smallest_patch => null() tmpptr => currentSite%youngest_patch do while(associated(tmpptr)) - write(fates_log(),*) tmpptr%area, tmpptr%nocomp_pft_label, tmpptr%land_use_label + if (tmpptr%land_use_label .eq. i_lulabel .and. & + .not.associated(tmpptr,largest_patch)) then + if (.not.associated(smallest_patch)) then + smallest_patch => tmpptr + else if (tmpptr%area .lt. smallest_patch%area) then + smallest_patch => tmpptr + endif + endif tmpptr => tmpptr%older end do - - call endrun(msg=errMsg(sourcefile, __LINE__)) + + if (associated(largest_patch) .and. associated(smallest_patch)) then + + warn_msg = 'fuse_patches could not reduce patch count for a land-use label; '// & + 'relabelling and fusing the smallest (orphan) patch into the '// & + 'largest patch of the same land-use label. '// & + 'lu label: '//trim(I2S(i_lulabel))// & + ' orphan nocomp pft: '//trim(I2S(smallest_patch%nocomp_pft_label))// & + ' orphan area: '//trim(N2S(smallest_patch%area)) + call FatesWarn(warn_msg,index=5) + + ! relabel BEFORE fusing so that the nocomp_pft_label check in + ! fuse_2_patches passes (donor and recipient labels must match) + smallest_patch%nocomp_pft_label = largest_patch%nocomp_pft_label + smallest_patch%age_since_anthro_disturbance = largest_patch%age_since_anthro_disturbance + call fuse_2_patches(csite, smallest_patch, largest_patch) + call fuse_cohorts(csite, largest_patch, bc_in) + call largest_patch%SortCohorts() + call largest_patch%ValidateCohorts() + + ! reset the tolerance and keep iterating; we removed one patch so the + ! count will eventually fall to/below the maximum + profiletol = ED_val_patch_fusion_tol + + else + + ! we could not identify a disposal target/donor pair; fall back to the + ! original diagnostic dump and crash + write(fates_log(),*) 'profile tolerance is too big, this shouldnt happen.' + write(fates_log(),*) 'probably this means there are too many distinct categorical ' + write(fates_log(),*) 'patch types for the maximum number of patches' + call dump_site(currentSite) + write(fates_log(),*) 'currentSite%area_bareground', currentSite%area_bareground + do i = 1, n_landuse_cats + write(fates_log(),*) 'i, currentSite%area_pft(:,i)',i, currentSite%area_pft(:,i) + end do + tmpptr => currentSite%youngest_patch + do while(associated(tmpptr)) + write(fates_log(),*) tmpptr%area, tmpptr%nocomp_pft_label, tmpptr%land_use_label + tmpptr => tmpptr%older + end do + + call endrun(msg=errMsg(sourcefile, __LINE__)) + + endif endif else iterate = 0 @@ -3322,6 +3449,12 @@ subroutine fuse_2_patches(csite, dp, rp) rp%c_stomata = (dp%c_stomata*dp%area + rp%c_stomata*rp%area) * inv_sum_area rp%c_lblayer = (dp%c_lblayer*dp%area + rp%c_lblayer*rp%area) * inv_sum_area + ! Conserve the land-use-change flag: if either patch still needs its nocomp PFT + ! identity remapped to the new land use type (changed_landuse_this_ts), the fused + ! patch must remain flagged so it is not silently left as an orphan (a patch whose + ! nocomp_pft_label is not valid for its land_use_label). + rp%changed_landuse_this_ts = rp%changed_landuse_this_ts .or. dp%changed_landuse_this_ts + rp%area = rp%area + dp%area !THIS MUST COME AT THE END! !insert donor cohorts into recipient patch @@ -3407,8 +3540,9 @@ subroutine terminate_patches(currentSite, bc_in) ! ! !DESCRIPTION: ! Terminate Patches if they are too small - ! - ! + ! USES: + use EDCohortDynamicsMod , only : terminate_cohorts + use FatesCohortMod , only : fates_cohort_type ! !ARGUMENTS: type(ed_site_type), target, intent(inout) :: currentSite type(bc_in_type), intent(in) :: bc_in @@ -3460,13 +3594,54 @@ subroutine terminate_patches(currentSite, bc_in) end do if ( .not. gotfused ) then - !! somehow didn't find a patch to fuse with. - warn_msg = 'small nocomp patch wasnt able to find '// & - 'another patch to fuse with. '// & - 'nocomp pft: '//trim(I2S(currentPatch%nocomp_pft_label))// & - 'lu label: '//trim(I2S(currentPatch%land_use_label))// & - 'area: '//trim(N2S(currentPatch%area)) - call FatesWarn(warn_msg,index=5) + !! somehow didn't find a patch of the same (nocomp_pft_label, land_use_label) + !! kind to fuse with. Left in place, this tiny patch becomes an irreducible + !! "orphan" that later blocks fuse_patches (its categorical labels match no + !! other patch in its land-use label). Dispose of it here by relabelling its + !! nocomp_pft_label to that of the largest patch in the same land-use label + !! and fusing that patch into it (keeping currentPatch as the recipient so the + !! enclosing loop's pointer remains valid). + + ! find the largest patch in the same land-use label (disposal partner) + largest_patch => null() + patchpointer => currentSite%youngest_patch + do while(associated(patchpointer)) + if ( .not.associated(patchpointer,currentPatch) .and. & + patchpointer%land_use_label .eq. currentPatch%land_use_label .and. & + patchpointer%nocomp_pft_label .ne. nocomp_bareground ) then + if (.not.associated(largest_patch)) then + largest_patch => patchpointer + else if (patchpointer%area .gt. largest_patch%area) then + largest_patch => patchpointer + endif + endif + patchpointer => patchpointer%older + end do + + if ( associated(largest_patch) ) then + warn_msg = 'small nocomp patch wasnt able to find a matching patch to '// & + 'fuse with; relabelling and fusing into the largest patch of '// & + 'the same land-use label. '// & + 'nocomp pft: '//trim(I2S(currentPatch%nocomp_pft_label))// & + 'lu label: '//trim(I2S(currentPatch%land_use_label))// & + 'area: '//trim(N2S(currentPatch%area)) + call FatesWarn(warn_msg,index=5) + + ! relabel currentPatch to the allowed pft before fusing so the + ! nocomp_pft_label check in fuse_2_patches passes + currentPatch%nocomp_pft_label = largest_patch%nocomp_pft_label + currentPatch%age_since_anthro_disturbance = largest_patch%age_since_anthro_disturbance + call fuse_2_patches(currentSite, largest_patch, currentPatch) + gotfused = .true. + else + !! truly nothing else in this land-use label to fuse into; warn as before + warn_msg = 'small nocomp patch wasnt able to find '// & + 'another patch to fuse with. '// & + 'nocomp pft: '//trim(I2S(currentPatch%nocomp_pft_label))// & + 'lu label: '//trim(I2S(currentPatch%land_use_label))// & + 'area: '//trim(N2S(currentPatch%area)) + call FatesWarn(warn_msg,index=5) + endif endif else nocomp_if @@ -3632,7 +3807,7 @@ subroutine terminate_patches(currentSite, bc_in) patchpointer%nocomp_pft_label = largest_patch%nocomp_pft_label patchpointer%age_since_anthro_disturbance = largest_patch%age_since_anthro_disturbance call fuse_2_patches(currentSite, patchpointer, largest_patch) - + ! start over in the loop to make sure we are removing every patch with the targeted land use type patchpointer => currentSite%youngest_patch @@ -3649,7 +3824,7 @@ subroutine terminate_patches(currentSite, bc_in) else write(fates_log(),*) 'this isnt because the land use was less than allowed' - !call endrun(msg=errMsg(sourcefile, __LINE__)) + call endrun(msg=errMsg(sourcefile, __LINE__)) ! Note to user. If you DO decide to remove the end-run above this line ! Make sure that you keep the pointer below this line, or you will get diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index bdf572039a..f5a26f86ba 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -262,6 +262,7 @@ module FatesPatchMod procedure :: Create procedure :: CountCohorts procedure :: ValidateCohorts + procedure :: CheckCohortsPfts procedure :: InsertCohort procedure :: SortCohorts procedure :: UpdateTreeGrassArea @@ -1068,13 +1069,14 @@ subroutine ValidateCohorts(this) ! DESCRIPTION: ! Validates a patch's cohort linked list ! - + use FatesInterfaceTypesMod , only : hlm_use_nocomp ! ARGUMENTS: class(fates_patch_type), intent(in), target :: this ! patch ! LOCALS: type(fates_cohort_type), pointer :: currentCohort ! cohort object integer :: forward_count, backward_count ! forwards and backwards counts of cohorts + logical :: debug=.false. ! check initial conditions if (.not. associated(this%shortest) .and. .not. associated(this%tallest)) then @@ -1139,9 +1141,52 @@ subroutine ValidateCohorts(this) additional_msg=errMsg(sourcefile, __LINE__)) return end if - + end subroutine ValidateCohorts + !=========================================================================== + + subroutine CheckCohortsPfts(this, mode, rc) + ! + ! DESCRIPTION: + ! IF we decide that Nocomp patches MUST have cohorts with %pft == nocom_pft_label + ! Here we either check or cohort%pfts + ! + use FatesInterfaceTypesMod , only : hlm_use_nocomp + + ! ARGUMENTS + class(fates_patch_type), intent(in), target :: this ! patch + integer, intent(in) :: mode ! mode: 0 skip check, 1 check, 2 assign + integer, intent(out), optional :: rc ! return code + + ! LOCAL VARIABLES + class(fates_cohort_type), pointer :: currentCohort + integer :: l_rc=0 + if ( (hlm_use_nocomp .eq. ifalse)) return + + currentCohort => this%shortest + do while (associated(currentCohort)) + if (this%nocomp_pft_label .ne. currentCohort%pft) then + l_rc = 1 + !if (mode == 1) then + write(fates_log(),*)'ERROR use_nocomp is true but a cohorts pft does not match nocomp_label (cpft,ncpft,c_area):', & + currentCohort%pft,this%nocomp_pft_label,this%land_use_label + + if (mode .ne. 0) then + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + if (mode == 2) then + currentCohort%pft = this%nocomp_pft_label + endif + end if + currentCohort => currentCohort%taller + end do + if (present(rc)) then + rc= l_rc + endif + return + end subroutine CheckCohortsPfts + !=========================================================================== subroutine CountCohorts(this) diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index 177557aedf..359189fbe0 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -9,6 +9,7 @@ module EDBtranMod use FatesConstantsMod , only : tfrz => t_water_freeze_k_1atm use FatesConstantsMod , only : itrue,ifalse,nearzero use FatesConstantsMod , only : nocomp_bareground + use FatesConstantsMod , only : rsnbl_math_prec use EDTypesMod , only : ed_site_type use FatesPatchMod, only : fates_patch_type use EDParamsMod, only : maxpft @@ -30,6 +31,7 @@ module EDBtranMod logical, parameter :: debug = .false. + character(len=*), parameter :: sourcefile = __FILE__ public :: btran_ed public :: get_active_suction_layers @@ -120,7 +122,9 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) real(r8) :: smp_node ! matrix potential real(r8) :: rresis ! suction limitation to transpiration independent ! of root density - real(r8) :: pftgs(maxpft) ! pft weighted stomatal conductance m/s + real(r8) :: pftgs(1:numpft) ! pft weighted stomatal conductance m/s + logical :: valid_pft(1:numpft) ! pft mask + real(r8) :: num_valid_pfts ! real(r8) :: temprootr real(r8) :: sum_pftgs ! sum of weighted conductances (for normalization) real(r8), allocatable :: root_resis(:,:) ! Root resistance in each pft x layer @@ -142,41 +146,89 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) ifp = cpatch%patchno - if_bare: if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then ! only for veg patches - - ! THIS SHOULD REALLY BE A COHORT LOOP ONCE WE HAVE rootfr_ft FOR COHORTS (RGK) + if_bare: if(cpatch%nocomp_pft_label.ne.nocomp_bareground .and. (cpatch%num_cohorts > 0))then ! only for veg patches + valid_pft(1:numpft) = .false. + pftgs(1:numpft)=0.0_r8 + sum_pftgs = 0.0_r8 + ccohort => cpatch%tallest + do while(associated(ccohort)) + if (ccohort%g_sb_laweight > rsnbl_math_prec) then + sum_pftgs = sum_pftgs + ccohort%g_sb_laweight + pftgs(ccohort%pft) = pftgs(ccohort%pft) + ccohort%g_sb_laweight + valid_pft(ccohort%pft) = .true. + endif + ccohort => ccohort%shorter + enddo + num_valid_pfts = 0.0_r8 do ft = 1,numpft - - call set_root_fraction(sites(s)%rootfrac_scr, ft, sites(s)%zi_soil, & - bc_in(s)%max_rooting_depth_index_col ) - - cpatch%btran_ft(ft) = 0.0_r8 - do j = 1,bc_in(s)%nlevsoil - - ! Calculations are only relevant where liquid water exists - ! see clm_fates%wrap_btran for calculation with CLM/ALM - - if ( check_layer_water(bc_in(s)%h2o_liqvol_sl(j),bc_in(s)%tempk_sl(j)) ) then - - smp_node = max(smpsc(ft), bc_in(s)%smp_sl(j)) - - rresis = min( (bc_in(s)%eff_porosity_sl(j)/bc_in(s)%watsat_sl(j))* & - (smp_node - smpsc(ft)) / (smpso(ft) - smpsc(ft)), 1._r8) - - root_resis(ft,j) = sites(s)%rootfrac_scr(j)*rresis - - ! root water uptake is not linearly proportional to root density, - ! to allow proper deep root funciton. Replace with equations from SPA/Newman. FIX(RF,032414) - - cpatch%btran_ft(ft) = cpatch%btran_ft(ft) + root_resis(ft,j) - + if(valid_pft(ft)) num_valid_pfts = num_valid_pfts +1._r8 + enddo + + do ft = 1,numpft + if (valid_pft(ft)) then + if (sum_pftgs < num_valid_pfts * rsnbl_math_prec) then + pftgs(ft) = 1._r8/num_valid_pfts ! condition above already means 1 valid pft else - root_resis(ft,j) = 0._r8 - end if - - end do !j + pftgs(ft) = pftgs(ft)/sum_pftgs + endif + else + pftgs(ft) = 0._r8 + endif + end do + + ! THIS SHOULD REALLY BE A COHORT LOOP ONCE WE HAVE rootfr_ft FOR COHORTS (RGK) + cpatch%btran_ft(:) = 0.0_r8 + root_resis(:,:) = 0.0_r8 + do ft = 1,numpft + if (valid_pft(ft)) then + call set_root_fraction(sites(s)%rootfrac_scr, ft, sites(s)%zi_soil, & + bc_in(s)%max_rooting_depth_index_col ) + + do j = 1,bc_in(s)%nlevsoil + + ! Calculations are only relevant where liquid water exists + ! see clm_fates%wrap_btran for calculation with CLM/ALM + if ( check_layer_water(bc_in(s)%h2o_liqvol_sl(j),bc_in(s)%tempk_sl(j)) ) then + + smp_node = max(smpsc(ft), bc_in(s)%smp_sl(j)) + + rresis = min( (bc_in(s)%eff_porosity_sl(j)/bc_in(s)%watsat_sl(j))* & + (smp_node - smpsc(ft)) / (smpso(ft) - smpsc(ft)), 1._r8) + + root_resis(ft,j) = root_resis(ft,j) + rresis * sites(s)%rootfrac_scr(j) + + ! root water uptake is not linearly proportional to root density, + ! to allow proper deep root funciton. Replace with equations from SPA/Newman. FIX(RF,032414) + + cpatch%btran_ft(ft) = cpatch%btran_ft(ft) + root_resis(ft,j) + + else + root_resis(ft,j) = 0._r8 + end if + end do !j + else + root_resis(ft,1:bc_in(s)%nlevsoil) = 0.0_r8 + cpatch%btran_ft(ft) = 0.0_r8 + endif ! valid_pfts + end do ! PFT + + ! remove this check when merging to noresm + if (debug) then + do ft=1,numpft + if (sum(root_resis(ft,1:bc_in(s)%nlevsoil)) .ne. cpatch%btran_ft(ft)) then + write(fates_log(),*) 'btran_ft not equal rresist ',cpatch%nocomp_pft_label,ft + write(fates_log(),*)' btran sum', cpatch%btran_ft(ft),sum(root_resis(ft,1:bc_in(s)%nlevsoil)) + write(fates_log(),*)' rootr ',root_resis(ft,1:bc_in(s)%nlevsoil) + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + enddo + endif + ! PFT-averaged point level root fraction for extraction purposese. + ! The cohort's conductance g_sb_laweighted, contains a weighting factor + ! based on the cohort's leaf area. units: [m/s] * [m2] + do ft = 1,numpft ! Normalize root resistances to get layer contribution to ET do j = 1,bc_in(s)%nlevsoil if (cpatch%btran_ft(ft) > nearzero) then @@ -185,63 +237,34 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) root_resis(ft,j) = 0._r8 end if end do - - end do !PFT - - ! PFT-averaged point level root fraction for extraction purposese. - ! The cohort's conductance g_sb_laweighted, contains a weighting factor - ! based on the cohort's leaf area. units: [m/s] * [m2] - - pftgs(1:maxpft) = 0._r8 - ccohort => cpatch%tallest - do while(associated(ccohort)) - pftgs(ccohort%pft) = pftgs(ccohort%pft) + ccohort%g_sb_laweight - ccohort => ccohort%shorter enddo - ! Process the boundary output, this is necessary for calculating the soil-moisture ! sink term across the different layers in driver/host. Photosynthesis will ! pass the host a total transpiration for the patch. This needs rootr to be ! distributed over the soil layers. - sum_pftgs = sum(pftgs(1:numpft)) - + bc_out(s)%rootr_pasl(ifp,:) = 0.0_r8 + bc_out(s)%btran_pa(ifp) = 0.0_r8 do j = 1, bc_in(s)%nlevsoil - bc_out(s)%rootr_pasl(ifp,j) = 0._r8 do ft = 1,numpft - if( sum_pftgs > 0._r8)then !prevent problem with the first timestep - might fail - !bit-retart test as a result? FIX(RF,032414) - bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j) + & - root_resis(ft,j) * pftgs(ft)/sum_pftgs - else - bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j) + & - root_resis(ft,j) * 1._r8/real(numpft,r8) - end if + bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j) + & + root_resis(ft,j) * pftgs(ft) enddo enddo - - ! Calculate the BTRAN that is passed back to the HLM - ! used only for diagnostics. If plant hydraulics is turned off - ! we are using the patchxpft level btran calculation - - if(hlm_use_planthydro.eq.ifalse) then - !weight patch level output BTRAN for the - bc_out(s)%btran_pa(ifp) = 0.0_r8 + if (hlm_use_planthydro.eq.ifalse) then do ft = 1,numpft - if( sum_pftgs > 0._r8)then !prevent problem with the first timestep - might fail - !bit-retart test as a result? FIX(RF,032414) - bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * pftgs(ft)/sum_pftgs - else - bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * 1./numpft - end if - enddo - end if + bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) *pftgs(ft) + end do + endif temprootr = sum(bc_out(s)%rootr_pasl(ifp,1:bc_in(s)%nlevsoil)) - if(abs(1.0_r8-temprootr) > 1.0e-10_r8 .and. temprootr > 1.0e-10_r8)then - - if(debug) write(fates_log(),*) 'error with rootr in canopy fluxes',temprootr,sum_pftgs - + if(abs(1.0_r8-temprootr) > rsnbl_math_prec .and. abs(temprootr) > rsnbl_math_prec)then + if (debug) then + write(fates_log(),*) 'error with rootr in canopy fluxes',temprootr,sum_pftgs,bc_out(s)%btran_pa(ifp) + ! remove this endrun later. + !call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + temprootr = abs(temprootr) do j = 1,bc_in(s)%nlevsoil bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j)/temprootr enddo diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index b599bd0961..92270a7eca 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -198,6 +198,12 @@ module FatesConstantsMod ! in nocomp simulations, what is the minimum PFT fraction for any given land use type? real(fates_r8), parameter, public :: min_nocomp_pftfrac_perlanduse = 0.01_fates_r8 + ! Relative (to the notional site area, or to a local area scale) floor below which a + ! patch/disturbance area is considered negligible. Used in place of absolute thresholds + ! (nearzero/rsnbl_math_prec) for patch-area comparisons so that the test scales with the + ! magnitude of the areas being compared. 1e-9 -> 1e-5 m2 at area=1e4. + real(fates_r8), parameter, public :: rel_patch_area_floor = 1.0e-9_fates_r8 + ! This is the precision of 8byte reals (~1e-308) real(fates_r8), parameter, public :: tinyr8 = tiny(1.0_fates_r8) @@ -378,9 +384,10 @@ module FatesConstantsMod ! some integers related to termination mortality - integer, parameter, public :: n_term_mort_types = 3 + integer, parameter, public :: n_term_mort_types = 4 integer, parameter, public :: i_term_mort_type_cstarv = 1 integer, parameter, public :: i_term_mort_type_canlev = 2 integer, parameter, public :: i_term_mort_type_numdens = 3 + integer, parameter, public :: i_term_mort_type_nocomp = 4 end module FatesConstantsMod diff --git a/parameter_files/fates_params_noresm.json b/parameter_files/fates_params_noresm.json index f375640c15..1c9ce04001 100644 --- a/parameter_files/fates_params_noresm.json +++ b/parameter_files/fates_params_noresm.json @@ -891,7 +891,7 @@ "dims": ["fates_pft"], "long_name": "background mortality rate", "units": "1/yr", - "data": [0.02, 0.02, 0.02, 0.004, 0.004, 0.004, 0.004, 0.004, 0.004, 0.004, 0.004, 0.014, 0.014, 0.05, 0.014] + "data": [0.02, 0.02, 0.02, 0.014, 0.004, 0.014, 0.004, 0.004, 0.004, 0.004, 0.004, 0.014, 0.014, 0.05, 0.014] }, "fates_mort_freezetol":{ "dtype": "float", @@ -926,7 +926,7 @@ "dims": ["fates_pft"], "long_name": "Mortality dbh senescence inflection point. If _ this mortality term is off. Setting this value turns on size dependent mortality", "units": "dbh cm", - "data": [140, 155, 175, 175, 140, 135, 175, null, null, null, null, null, null, null, null] + "data": [122, 135, 175, 130, 140, 120, 175, null, null, null, null, null, null, null, null] }, "fates_mort_prescribed_canopy":{ "dtype": "float", @@ -954,7 +954,7 @@ "dims": ["fates_pft"], "long_name": "Mortality dbh senescence rate of change. Sensible range is around 0.03-0.06. Larger values give steeper mortality curves.", "units": "mortality rate dbh^-1", - "data": [0.055, 0.055, 0.04, 0.04, 0.055, 0.055, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04] + "data": [0.059, 0.061, 0.04, 0.057, 0.055, 0.059, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04] }, "fates_mort_scalar_coldstress":{ "dtype": "float",