diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 196b40ba15..696a88514c 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -3815,156 +3815,104 @@ subroutine InsertPatch(currentSite, newPatch) ! !LOCAL VARIABLES: type (fates_patch_type), pointer :: currentPatch - integer :: insert_method ! Temporary dev - logical :: found_landuselabel_match - integer, parameter :: unordered_lul_groups= 1 - integer, parameter :: primaryland_oldest_group = 2 - integer, parameter :: numerical_order_lul_groups = 3 - integer, parameter :: age_order_only = 4 - - ! Insert new patch case options: - ! Option 1: Group the landuse types together, but the group order doesn't matter - ! Option 2: Option 1, but primarylands are forced to be the oldest group - ! Option 3: Option 1, but groups are in numerical order according to land use label index integer - ! (i.e. primarylands=1, secondarylands=2, ..., croplands=5) - ! Option 4: Don't group the patches by land use label. Simply add new patches to the youngest end. - - ! Hardcode the default insertion method. The options developed during FATES V1 land use are - ! currently being held for potential future usage. - insert_method = primaryland_oldest_group - - ! Start from the youngest patch and work to oldest, regarless of insertion_method - currentPatch => currentSite%youngest_patch + logical :: patch_inserted - ! For the three grouped cases, if the land use label of the youngest patch on the site - ! is a match to the new patch land use label, simply insert it as the new youngest. - ! This is applicable to the non-grouped option 4 method as well. - if (currentPatch%land_use_label .eq. newPatch%land_use_label ) then - newPatch%older => currentPatch - newPatch%younger => null() - currentPatch%younger => newPatch - currentSite%youngest_patch => newPatch - else - - ! If the current site youngest patch land use label doesn't match the new patch - ! land use label then work through the list until you find the matching type. - ! Since we've just checked the youngest patch, move to the next patch and - ! initialize the match flag to false. - found_landuselabel_match = .false. - currentPatch => currentPatch%older - select case(insert_method) + ! The goal here is to have patches ordered in a specific way. That way is to have them + ! looped as the following, where LU refers to the land use label, PFT refers to the + ! nocomp PFT label, and Y and O refer to continuous patch ages. + ! + ! LU1 ---- LU2 ---- LU3 -- etc + ! / \ / \ / \ + ! PFT1 --- PFT2 | PFT1 --- PFT2 | PFT1 --- PFT2 -- etc + ! / \ / \ / \ / \ / \ / \ + ! O - Y O - Y O - Y O - Y O - Y O - Y -- etc - ! Option 1 - order of land use label groups does not matter - case (unordered_lul_groups) + ! I.e. to treat land use as the outermost loop element, then nocomp PFT as next loop element, + ! and then age as the innermost loop element. Visualizing the above as a linked list patches: - do while(associated(currentPatch) .and. .not. found_landuselabel_match) - if (currentPatch%land_use_label .eq. newPatch%land_use_label) then - found_landuselabel_match = .true. - else - currentPatch => currentPatch%older - end if - end do + ! LU1/PFT1/O <-> LU1/PFT1/Y <-> LU1/PFT2/O <- ... -> LU3/PFT2/O <-> LU3/PFT2/Y - ! In the case where we've found a land use label matching the new patch label, - ! insert the newPatch will as the youngest patch for that land use type. - if (associated(currentPatch)) then - newPatch%older => currentPatch - newPatch%younger => currentPatch%younger - currentPatch%younger%older => newPatch - currentPatch%younger => newPatch - else - ! In the case in which we get to the end of the list and haven't found - ! a landuse label match simply add the new patch to the youngest end. - newPatch%older => currentSite%youngest_patch - newPatch%younger => null() - currentSite%youngest_patch%younger => newPatch - currentSite%youngest_patch => newPatch - endif - - ! Option 2 - primaryland group must be on the oldest end - case (primaryland_oldest_group) + ! Mapping this setup onto the existing "older/younger" scheme means that lower number + ! land use and pft labels are considered "older". Note that this matches the current + ! initialization scheme in which patches are created and linked in increasing pft + ! numerical order starting from 1. This also aligns with the current set_patchno scheme + ! in which patches are given an indexable number for the API iteration loops. + + ! The way to accomplsh this most simply is to define a pseudo-age that includes all of the + ! above info and sort the patches based on the pseudo-age. i.e. take some number larger + ! than any patch will ever reach in actual age. Then take the LU, multiply it by the big + ! number squared, add it to the pft number multiplied by the big number, and add to the age. + ! And lastly to sort using that instead of the actual age. - do while(associated(currentPatch) .and. .not. found_landuselabel_match) - if (currentPatch%land_use_label .eq. newPatch%land_use_label) then - found_landuselabel_match = .true. - else - currentPatch => currentPatch%older - end if - end do + ! If land use is turned off or nocomp is turned off, then this should devolve to the prior + ! behavior of just age sorting. - ! In the case where we've found a land use label matching the new patch label, - ! insert the newPatch will as the youngest patch for that land use type. - if (associated(currentPatch)) then - newPatch%older => currentPatch - newPatch%younger => currentPatch%younger - currentPatch%younger%older => newPatch - currentPatch%younger => newPatch - else - ! In the case in which we get to the end of the list and haven't found - ! a landuse label match. - - ! If the new patch is primaryland add it to the oldest end of the list - if (newPatch%land_use_label .eq. primaryland) then - newPatch%older => null() - newPatch%younger => currentSite%oldest_patch - currentSite%oldest_patch%older => newPatch - currentSite%oldest_patch => newPatch - else - ! If the new patch land use type is not primaryland and we are at the - ! oldest end of the list, add it to the youngest end - newPatch%older => currentSite%youngest_patch - newPatch%younger => null() - currentSite%youngest_patch%younger => newPatch - currentSite%youngest_patch => newPatch - endif - endif + patch_inserted = .false. + + if (GetPseudoPatchAge(newPatch) .le. GetPseudoPatchAge(currentSite%youngest_patch)) then - ! Option 3 - groups are numerically ordered with primaryland group starting at oldest end. - case (numerical_order_lul_groups) + ! insert new patch at the head of the linked list + newPatch%older => currentSite%youngest_patch + newPatch%younger => null() + currentSite%youngest_patch%younger => newPatch + currentSite%youngest_patch => newPatch - ! If the youngest patch landuse label number is greater than the new - ! patch land use label number, the new patch must be inserted somewhere - ! in between oldest and youngest - do while(associated(currentPatch) .and. .not. found_landuselabel_match) - if (currentPatch%land_use_label .eq. newPatch%land_use_label .or. & - currentPatch%land_use_label .lt. newPatch%land_use_label) then - found_landuselabel_match = .true. - else - currentPatch => currentPatch%older - endif - end do + patch_inserted = .true. + else if (GetPseudoPatchAge(newPatch) .ge. GetPseudoPatchAge(currentSite%oldest_patch)) then - ! In the case where we've found a landuse label matching the new patch label - ! insert the newPatch will as the youngest patch for that land use type. - if (associated(currentPatch)) then + ! insert new patch at the end of the linked list + newPatch%younger => currentSite%oldest_patch + newPatch%older => null() + currentSite%oldest_patch%older => newPatch + currentSite%oldest_patch => newPatch - newPatch%older => currentPatch - newPatch%younger => currentPatch%younger + patch_inserted = .true. + else + ! new patch has a pseudo-age somewhere within the linked list. find the first patch which + ! has a pseudo age older than it, and put it ahead of that patch + currentPatch => currentSite%youngest_patch + do while (associated(currentPatch) .and. ( .not. patch_inserted) ) + if (GetPseudoPatchAge(newPatch) .lt. GetPseudoPatchAge(currentPatch)) then + newPatch%older => currentPatch + newPatch%younger => currentPatch%younger currentPatch%younger%older => newPatch - currentPatch%younger => newPatch - - else - - ! In the case were we get to the end, the new patch - ! must be numerically the smallest, so put it at the oldest position - newPatch%older => null() - newPatch%younger => currentSite%oldest_patch - currentSite%oldest_patch%older => newPatch - currentSite%oldest_patch => newPatch + currentPatch%younger => newPatch + patch_inserted = .true. endif + currentPatch => currentPatch%older + end do + end if - ! Option 4 - always add the new patch as the youngest regardless of land use label - case (age_order_only) - ! Set the current patch to the youngest patch - newPatch%older => currentSite%youngest_patch - newPatch%younger => null() - currentSite%youngest_patch%younger => newPatch - currentSite%youngest_patch => newPatch - end select + if ( .not. patch_inserted) then + ! something has gone wrong. abort. + write(fates_log(),*) 'something has gone wrong in the patch insertion, because no place to put the new patch was found' + call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end subroutine InsertPatch + end subroutine InsertPatch + + ! ===================================================================================== + + function GetPseudoPatchAge(CurrentPatch) result(pseudo_age) + + ! Purpose: we want to sort the patches in a way that takes into account both their + ! continuous and categorical variables. Calculate a pseudo age that does this, by taking + ! the integer labels, multiplying these by large numbers, and adding to the continuous age. + ! Note to ensure that lower integer land use label and pft label numbers are considered + ! "younger" (i.e higher index patchno) in the linked list, they are summed and multiplied by + ! negative one. The patch age is still added normally to this negative pseudoage calculation + ! as a higher age will result in a less negative number correlating with an "older" patch. + + type (fates_patch_type), intent(in), pointer :: CurrentPatch + real(r8) :: pseudo_age + real(r8), parameter :: max_actual_age = 1.e4 ! hard to imagine a patch older than 10,000 years + real(r8), parameter :: max_actual_age_squared = 1.e8 + + pseudo_age = -1.0_r8 * (real(CurrentPatch%land_use_label,r8) * max_actual_age_squared + & + real(CurrentPatch%nocomp_pft_label,r8) * max_actual_age) + CurrentPatch%age + + end function GetPseudoPatchAge ! ===================================================================================== diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index d529a0e123..f56e3e0dd2 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -706,6 +706,12 @@ subroutine init_patches( nsites, sites, bc_in) do el=1,num_elements call SiteMassStock(sites(s),el,sites(s)%mass_balance(el)%old_stock, & biomass_stock,litter_stock,seed_stock) + ! Initialize the integrated flux balance diagnostics + ! No need to initialize the instantaneous states, those are re-calculated + sites(s)%iflux_balance(el)%iflux_liveveg = & + (biomass_stock + seed_stock)*area_inv + sites(s)%iflux_balance(el)%iflux_litter = litter_stock * area_inv + end do call set_patchno(sites(s)) enddo @@ -972,7 +978,7 @@ subroutine init_patches( nsites, sites, bc_in) do el=1,num_elements call SiteMassStock(sites(s),el,sites(s)%mass_balance(el)%old_stock, & biomass_stock,litter_stock,seed_stock) - + ! Initialize the integrated flux balance diagnostics ! No need to initialize the instantaneous states, those are re-calculated sites(s)%iflux_balance(el)%iflux_liveveg = & @@ -983,7 +989,7 @@ subroutine init_patches( nsites, sites, bc_in) call set_patchno(sites(s)) - enddo sites_loop !s + enddo sites_loop end if ! zero all the patch fire variables for the first timestep diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index e69d1e9bcb..1e7ee14e13 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -834,7 +834,10 @@ subroutine ed_update_site( currentSite, bc_in, bc_out, is_restarting ) enddo ! Check to see if the time integrated fluxes match the state - call CheckIntegratedMassPools(currentSite) + ! Dont call this if we are restarting, it will double count the flux + if(.not.is_restarting)then + call CheckIntegratedMassPools(currentSite) + end if ! The HLMs need to know about nutrient demand, and/or ! root mass and affinities diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 739ea30503..364ad816cf 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -53,6 +53,7 @@ module FatesHistoryInterfaceMod use FatesInterfaceTypesMod , only : numpft use FatesInterfaceTypesMod , only : hlm_freq_day use FatesInterfaceTypesMod , only : hlm_parteh_mode + use FatesInterfaceTypesMod , only : hlm_use_sp use EDParamsMod , only : ED_val_comp_excln use EDParamsMod , only : ED_val_phen_coldtemp use EDParamsMod , only : nlevleaf @@ -3101,9 +3102,7 @@ subroutine update_history_dyn2(this,nc,nsites,sites,bc_in) real(r8), parameter :: reallytalltrees = 1000. ! some large number (m) -! - associate( hio_err_fates_elem => this%hvars(ih_err_fates_elem)%r82d, & - hio_biomass_si_pft => this%hvars(ih_biomass_si_pft)%r82d, & + associate( hio_biomass_si_pft => this%hvars(ih_biomass_si_pft)%r82d, & hio_biomass_sec_si_pft => this%hvars(ih_biomass_sec_si_pft)%r82d, & hio_leafbiomass_si_pft => this%hvars(ih_leafbiomass_si_pft)%r82d, & hio_storebiomass_si_pft => this%hvars(ih_storebiomass_si_pft)%r82d, & @@ -3315,11 +3314,8 @@ subroutine update_history_dyn2(this,nc,nsites,sites,bc_in) hio_nplant_understory_si_scag => this%hvars(ih_nplant_understory_si_scag)%r82d, & hio_disturbance_rate_si_lulu => this%hvars(ih_disturbance_rate_si_lulu)%r82d, & hio_cstarvmortality_continuous_carbonflux_si_pft => this%hvars(ih_cstarvmortality_continuous_carbonflux_si_pft)%r82d, & - hio_interr_liveveg_elem => this%hvars(ih_interr_liveveg_elem)%r82d, & - hio_interr_litter_elem => this%hvars(ih_interr_litter_elem)%r82d, & hio_transition_matrix_si_lulu => this%hvars(ih_transition_matrix_si_lulu)%r82d) - model_day_int = nint(hlm_model_day) ! --------------------------------------------------------------------------------- @@ -3363,13 +3359,17 @@ subroutine update_history_dyn2(this,nc,nsites,sites,bc_in) do el = 1, num_elements - ! Total model error [kg/day -> kg/s] (all elements) - hio_err_fates_elem(io_si,el) = sites(s)%mass_balance(el)%err_fates / sec_per_day + if((hlm_use_ed_st3 .eq. ifalse) .and. (hlm_use_sp .eq. ifalse)) then - hio_interr_liveveg_elem(io_si,el) = sites(s)%flux_diags%elem(el)%err_liveveg + ! Total model error [kg/day -> kg/s] (all elements) + this%hvars(ih_err_fates_elem)%r82d(io_si,el) = sites(s)%mass_balance(el)%err_fates / sec_per_day + + this%hvars(ih_interr_liveveg_elem)%r82d(io_si,el) = sites(s)%flux_diags%elem(el)%err_liveveg - hio_interr_litter_elem(io_si,el) = sites(s)%flux_diags%elem(el)%err_litter - + this%hvars(ih_interr_litter_elem)%r82d(io_si,el) = sites(s)%flux_diags%elem(el)%err_litter + + end if + ! Total element lost to atmosphere from burning (kg/site/day -> kg/m2/s) hio_burn_flux_elem(io_si,el) = & sites(s)%mass_balance(el)%burn_flux_to_atm * ha_per_m2 * & @@ -8462,25 +8462,26 @@ subroutine define_history_vars(this, initialize_variables) hlms='CLM:ALM', upfreq=group_dyna_complx, ivar=ivar, & initialize=initialize_variables, index = ih_burn_flux_elem) - - call this%set_history_var(vname='FATES_ERROR_EL', units='kg s-1', & - long='total mass-balance error in kg per second by element', & - use_default='active', avgflag='A', vtype=site_elem_r8, & - hlms='CLM:ALM', upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & - index = ih_err_fates_elem) - - call this%set_history_var(vname='FATES_INTERR_LIVEVEG_EL',units='kg m-2', & - long='Bias error between integrated flux and (minus) state in live vegetation ', & - use_default='active', avgflag='A', vtype=site_elem_r8, hlms='CLM:ALM', & - upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & - index = ih_interr_liveveg_elem) - - call this%set_history_var(vname='FATES_INTERR_LITTER_EL',units='kg m-2', & - long='Bias error between integrated flux and (minus) state in litter ', & - use_default='active', avgflag='A', vtype=site_elem_r8, hlms='CLM:ALM', & - upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & - index = ih_interr_litter_elem) - + if((hlm_use_ed_st3 .eq. ifalse) .and. (hlm_use_sp .eq. ifalse))then + call this%set_history_var(vname='FATES_ERROR_EL', units='kg s-1', & + long='total mass-balance error in kg per second by element', & + use_default='active', avgflag='A', vtype=site_elem_r8, & + hlms='CLM:ALM', upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & + index = ih_err_fates_elem) + + call this%set_history_var(vname='FATES_INTERR_LIVEVEG_EL',units='kg m-2', & + long='Bias error between integrated flux and (minus) state in live vegetation ', & + use_default='active', avgflag='A', vtype=site_elem_r8, hlms='CLM:ALM', & + upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & + index = ih_interr_liveveg_elem) + + call this%set_history_var(vname='FATES_INTERR_LITTER_EL',units='kg m-2', & + long='Bias error between integrated flux and (minus) state in litter ', & + use_default='active', avgflag='A', vtype=site_elem_r8, hlms='CLM:ALM', & + upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & + index = ih_interr_litter_elem) + end if + call this%set_history_var(vname='FATES_LITTER_AG_FINE_EL', units='kg m-2', & long='mass of aboveground litter in fines (leaves, nonviable seed) by element', & use_default='active', avgflag='A', vtype=site_elem_r8, &