Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

compatability changes for clm bgc #999

Merged
merged 10 commits into from
Aug 14, 2023
128 changes: 63 additions & 65 deletions biogeochem/FatesSoilBGCFluxMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -255,7 +255,6 @@ subroutine PrepCH4BCs(csite,bc_in,bc_out)
integer :: fp ! patch index of the site
real(r8) :: agnpp ! Above ground daily npp
real(r8) :: bgnpp ! Below ground daily npp
real(r8) :: site_npp ! Site level NPP gC/m2/year
real(r8) :: plant_area ! crown area (m2) of all plants in patch
real(r8) :: woody_area ! corwn area (m2) of woody plants in patch
real(r8) :: fnrt_c ! Fine root carbon [kg/plant]
Expand All @@ -268,8 +267,9 @@ subroutine PrepCH4BCs(csite,bc_in,bc_out)

real(r8), parameter :: ema_npp_tscale = 10._r8 ! 10 day


! Exit if we need not communicate with the hlm's ch4 module
if(.not.(hlm_use_ch4==itrue) .and. .not.(hlm_parteh_mode==prt_cnp_flex_allom_hyp) ) return
! if(.not.(hlm_use_ch4==itrue) .and. .not.(hlm_parteh_mode==prt_cnp_flex_allom_hyp) ) return

! Initialize to zero
bc_out%annavg_agnpp_pa(:) = 0._r8
Expand All @@ -279,12 +279,15 @@ subroutine PrepCH4BCs(csite,bc_in,bc_out)
bc_out%frootc_pa(:) = 0._r8
bc_out%root_resp(:) = 0._r8
bc_out%woody_frac_aere_pa(:) = 0._r8
site_npp = 0._r8
bc_out%ema_npp = 0._r8

! Process CH4 variables first
!if(.not.(hlm_use_ch4==itrue) .and. .not.(hlm_parteh_mode==prt_cnp_flex_allom_hyp) )

fp = 0
cpatch => csite%oldest_patch
do while (associated(cpatch))

! Patch ordering when passing boundary conditions
! always goes from oldest to youngest, following
! the convention of EDPatchDynamics::set_patchno()
Expand All @@ -298,55 +301,58 @@ subroutine PrepCH4BCs(csite,bc_in,bc_out)

ccohort => cpatch%tallest
do while (associated(ccohort))

! For consistency, only apply calculations to non-new
! cohorts. New cohorts will not have respiration rates
! at this point in the call sequence.

if(.not.ccohort%isnew) then

pft = ccohort%pft

call set_root_fraction(csite%rootfrac_scr, pft, csite%zi_soil, &
bc_in%max_rooting_depth_index_col )

fnrt_c = ccohort%prt%GetState(fnrt_organ, carbon12_element)

! Fine root fraction over depth

bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) = &
bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) + &
csite%rootfrac_scr(1:bc_in%nlevsoil)

! Fine root carbon, convert [kg/plant] -> [g/m2]
bc_out%frootc_pa(fp) = &
bc_out%frootc_pa(fp) + &
fnrt_c*ccohort%n/cpatch%area * g_per_kg


! [kgC/day]
sapw_net_alloc = ccohort%prt%GetNetAlloc(sapw_organ, carbon12_element) * days_per_sec
store_net_alloc = ccohort%prt%GetNetAlloc(store_organ, carbon12_element) * days_per_sec
leaf_net_alloc = ccohort%prt%GetNetAlloc(leaf_organ, carbon12_element) * days_per_sec
fnrt_net_alloc = ccohort%prt%GetNetAlloc(fnrt_organ, carbon12_element) * days_per_sec
struct_net_alloc = ccohort%prt%GetNetAlloc(struct_organ, carbon12_element) * days_per_sec
repro_net_alloc = ccohort%prt%GetNetAlloc(repro_organ, carbon12_element) * days_per_sec

! [kgC/plant/day] -> [gC/m2/s]
agnpp = agnpp + ccohort%n/cpatch%area * (leaf_net_alloc + repro_net_alloc + &
prt_params%allom_agb_frac(pft)*(sapw_net_alloc+store_net_alloc+struct_net_alloc)) * g_per_kg

! [kgC/plant/day] -> [gC/m2/s]
bgnpp = bgnpp + ccohort%n/cpatch%area * (fnrt_net_alloc + &
(1._r8-prt_params%allom_agb_frac(pft))*(sapw_net_alloc+store_net_alloc+struct_net_alloc)) * g_per_kg

! (gC/m2/s) root respiration (fine root MR + total root GR)
! RGK: We do not save root respiration and average over the day. Until we do
! this is a best (bad) guess at fine root MR + total root GR
! (kgC/indiv/yr) -> gC/m2/s
bc_out%root_resp(1:bc_in%nlevsoil) = bc_out%root_resp(1:bc_in%nlevsoil) + &
ccohort%resp_acc_hold*years_per_day*g_per_kg*days_per_sec* &
ccohort%n*area_inv*(1._r8-prt_params%allom_agb_frac(pft)) * csite%rootfrac_scr(1:bc_in%nlevsoil)


if(hlm_use_ch4==itrue)then

! Fine root fraction over depth
bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) = &
bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) + &
csite%rootfrac_scr(1:bc_in%nlevsoil)

! Fine root carbon, convert [kg/plant] -> [g/m2]
bc_out%frootc_pa(fp) = &
bc_out%frootc_pa(fp) + &
fnrt_c*ccohort%n/cpatch%area * g_per_kg

! (gC/m2/s) root respiration (fine root MR + total root GR)
! RGK: We do not save root respiration and average over the day. Until we do
! this is a best (bad) guess at fine root MR + total root GR
! (kgC/indiv/yr) -> gC/m2/s
bc_out%root_resp(1:bc_in%nlevsoil) = bc_out%root_resp(1:bc_in%nlevsoil) + &
ccohort%resp_acc_hold*years_per_day*g_per_kg*days_per_sec* &
ccohort%n*area_inv*(1._r8-prt_params%allom_agb_frac(pft)) * csite%rootfrac_scr(1:bc_in%nlevsoil)

end if

if( prt_params%woody(pft)==itrue ) then
woody_area = woody_area + ccohort%c_area
end if
Expand All @@ -357,44 +363,36 @@ subroutine PrepCH4BCs(csite,bc_in,bc_out)

ccohort => ccohort%shorter
end do

if( sum(bc_out%rootfr_pa(fp,1:bc_in%nlevsoil)) > nearzero) then
bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) = &
bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) / &
sum(bc_out%rootfr_pa(fp,1:bc_in%nlevsoil))
end if

! RGK: These averages should switch to the new patch averaging methods
! when available. Right now we are not doing any time averaging
! because it would be mixing the memory of patches, which
! would be arguably worse than just using the instantaneous value

! gC/m2/s
bc_out%annavg_agnpp_pa(fp) = agnpp
bc_out%annavg_bgnpp_pa(fp) = bgnpp
! gc/m2/yr
bc_out%annsum_npp_pa(fp) = (bgnpp+agnpp)*days_per_year*sec_per_day

site_npp = site_npp + bc_out%annsum_npp_pa(fp)*cpatch%area*area_inv

if(plant_area>nearzero) then
bc_out%woody_frac_aere_pa(fp) = woody_area/plant_area
if(hlm_use_ch4==itrue)then
if( sum(bc_out%rootfr_pa(fp,1:bc_in%nlevsoil)) > nearzero) then
bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) = &
bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) / &
sum(bc_out%rootfr_pa(fp,1:bc_in%nlevsoil))
end if

! RGK: These averages should switch to the new patch averaging methods
! when available. Right now we are not doing any time averaging
! because it would be mixing the memory of patches, which
! would be arguably worse than just using the instantaneous value

! gC/m2/s
bc_out%annavg_agnpp_pa(fp) = agnpp
bc_out%annavg_bgnpp_pa(fp) = bgnpp
! gc/m2/yr
bc_out%annsum_npp_pa(fp) = (bgnpp+agnpp)*days_per_year*sec_per_day

if(plant_area>nearzero) then
bc_out%woody_frac_aere_pa(fp) = woody_area/plant_area
end if

end if

cpatch => cpatch%younger
end do

! Smoothed [gc/m2/yr]
if(csite%ema_npp<-9000._r8)then
! For cold starts, we initialize the ema_npp value at -9999, so that
! it becomes memoryless and uses the first calculated value
csite%ema_npp = site_npp
else
csite%ema_npp = (1._r8-1._r8/ema_npp_tscale)*csite%ema_npp + (1._r8/ema_npp_tscale)*site_npp
end if

bc_out%ema_npp = csite%ema_npp

return
end subroutine PrepCH4BCs

Expand Down Expand Up @@ -700,7 +698,7 @@ subroutine FluxIntoLitterPools(csite, bc_in, bc_out)
do id = 1,nlev_eff_decomp
surface_prof(id) = surface_prof(id)/surface_prof_tot
end do

! Loop over the different elements.
do el = 1, num_elements

Expand All @@ -710,9 +708,9 @@ subroutine FluxIntoLitterPools(csite, bc_in, bc_out)

select case (element_list(el))
case (carbon12_element)
bc_out%litt_flux_cel_c_si(:) = 0._r8
bc_out%litt_flux_lig_c_si(:) = 0._r8
bc_out%litt_flux_lab_c_si(:) = 0._r8
bc_out%litt_flux_cel_c_si(:) = 0.0_r8
bc_out%litt_flux_lig_c_si(:) = 0.0_r8
bc_out%litt_flux_lab_c_si(:) = 0.0_r8
flux_cel_si => bc_out%litt_flux_cel_c_si(:)
flux_lab_si => bc_out%litt_flux_lab_c_si(:)
flux_lig_si => bc_out%litt_flux_lig_c_si(:)
Expand Down
6 changes: 3 additions & 3 deletions main/EDMainMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -736,9 +736,9 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
enddo


! Before we start messing with the patch areas, and before we start removing
! trees, this is a good time to pass fragmentation litter fluxes and
! plant-to-soil fluxes (such as efflux and fixation fluxes)
! RGK: This call is unecessary for CLM coupling. I believe we
! can remove it completely if/when this call is added in ELM to
! subroutine UpdateLitterFluxes(this,bounds_clump) in elmfates_interfaceMod.F90

call FluxIntoLitterPools(currentsite, bc_in, bc_out)

Expand Down
53 changes: 38 additions & 15 deletions main/FatesInterfaceMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ module FatesInterfaceMod
use FatesConstantsMod , only : sec_per_day
use FatesConstantsMod , only : days_per_year
use FatesConstantsMod , only : TRS_regeneration
use FatesConstantsMod , only : g_per_kg
use FatesGlobals , only : fates_global_verbose
use FatesGlobals , only : fates_log
use FatesGlobals , only : endrun => fates_endrun
Expand Down Expand Up @@ -328,6 +329,12 @@ subroutine zero_bcs(fates,s)
fates%bc_out(s)%litt_flux_cel_c_si(:) = 0._r8
fates%bc_out(s)%litt_flux_lig_c_si(:) = 0._r8
fates%bc_out(s)%litt_flux_lab_c_si(:) = 0._r8
! Yes, zero out N flux arrays for c-only runs.
! This is because we want these on (and zero)
! with CLM.
fates%bc_out(s)%litt_flux_cel_n_si(:) = 0._r8
fates%bc_out(s)%litt_flux_lig_n_si(:) = 0._r8
fates%bc_out(s)%litt_flux_lab_n_si(:) = 0._r8
case(prt_cnp_flex_allom_hyp)

fates%bc_in(s)%plant_nh4_uptake_flux(:,:) = 0._r8
Expand Down Expand Up @@ -458,7 +465,6 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats,
allocate(bc_in%plant_p_uptake_flux(1,1))
end if


allocate(bc_in%zi_sisl(0:nlevsoil_in))
allocate(bc_in%dz_sisl(nlevsoil_in))
allocate(bc_in%z_sisl(nlevsoil_in))
Expand Down Expand Up @@ -615,7 +621,7 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in)

! Include the bare-ground patch for these patch-level boundary conditions
! (it will always be zero for all of these)
if(hlm_use_ch4.eq.itrue) then
!if(hlm_use_ch4.eq.itrue) then
allocate(bc_out%annavg_agnpp_pa(0:maxpatch_total));bc_out%annavg_agnpp_pa(:)=nan
allocate(bc_out%annavg_bgnpp_pa(0:maxpatch_total));bc_out%annavg_bgnpp_pa(:)=nan
allocate(bc_out%annsum_npp_pa(0:maxpatch_total));bc_out%annsum_npp_pa(:)=nan
Expand All @@ -627,19 +633,23 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in)

! Give the bare-ground root fractions a nominal fraction of unity over depth
bc_out%rootfr_pa(0,1:nlevsoil_in)=1._r8/real(nlevsoil_in,r8)
end if
!end if

bc_out%ema_npp = -9999.9_r8
bc_out%ema_npp = nan

! Fates -> BGC fragmentation mass fluxes
select case(hlm_parteh_mode)
case(prt_carbon_allom_hyp)
allocate(bc_out%litt_flux_cel_c_si(nlevdecomp_in))
allocate(bc_out%litt_flux_lig_c_si(nlevdecomp_in))
allocate(bc_out%litt_flux_lab_c_si(nlevdecomp_in))
! Yes, allocate N flux arrays for c-only runs.
! This is because we want these on (and zero)
! with CLM.
allocate(bc_out%litt_flux_cel_n_si(nlevdecomp_in))
allocate(bc_out%litt_flux_lig_n_si(nlevdecomp_in))
allocate(bc_out%litt_flux_lab_n_si(nlevdecomp_in))
case(prt_cnp_flex_allom_hyp)


allocate(bc_out%litt_flux_cel_c_si(nlevdecomp_in))
allocate(bc_out%litt_flux_lig_c_si(nlevdecomp_in))
allocate(bc_out%litt_flux_lab_c_si(nlevdecomp_in))
Expand All @@ -649,10 +659,8 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in)
allocate(bc_out%litt_flux_cel_p_si(nlevdecomp_in))
allocate(bc_out%litt_flux_lig_p_si(nlevdecomp_in))
allocate(bc_out%litt_flux_lab_p_si(nlevdecomp_in))

allocate(bc_out%source_nh4(nlevdecomp_in))
allocate(bc_out%source_p(nlevdecomp_in))

case default
write(fates_log(), *) 'An unknown parteh hypothesis was passed'
write(fates_log(), *) 'to the site level output boundary conditions'
Expand Down Expand Up @@ -1956,6 +1964,7 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in)
type(fates_patch_type), pointer :: cpatch
type(fates_cohort_type), pointer :: ccohort
integer :: s, ifp, io_si, pft
real(r8) :: site_npp ! Site level NPP gC/m2/year
real(r8) :: new_seedling_layer_par ! seedling layer par in the current timestep
real(r8) :: new_seedling_layer_smp ! seedling layer smp in the current timestep
real(r8) :: new_seedling_mdd ! seedling layer moisture deficit days in the current timestep
Expand All @@ -1965,10 +1974,12 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in)
real(r8) :: seedling_par_low ! lower intensity par for seedlings (par under the undergrowth) [W/m2]
real(r8) :: par_low_frac ! fraction of ground where PAR is low
integer,parameter :: ipar = 1 ! solar radiation in the shortwave band (i.e. par)
real(r8), parameter :: ema_npp_tscale = 480._r8 ! 10 day (10*48 steps)

do s = 1,size(sites,dim=1)

ifp=0
site_npp = 0._r8
cpatch => sites(s)%oldest_patch
do while(associated(cpatch))
if (cpatch%patchno .ne. 0) then
Expand All @@ -1977,8 +1988,6 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in)
call cpatch%tveg_lpa%UpdateRMean(bc_in(s)%t_veg_pa(ifp))
call cpatch%tveg_longterm%UpdateRMean(bc_in(s)%t_veg_pa(ifp))



! Update the seedling layer par running means
if ( regeneration_model == TRS_regeneration ) then

Expand Down Expand Up @@ -2024,15 +2033,29 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in)

end if

!ccohort => cpatch%tallest
!do while (associated(ccohort))
! call ccohort%tveg_lpa%UpdateRMean(bc_in(s)%t_veg_pa(ifp))
! ccohort => ccohort%shorter
!end do
ccohort => cpatch%tallest
do while (associated(ccohort))
! call ccohort%tveg_lpa%UpdateRMean(bc_in(s)%t_veg_pa(ifp))

! [kgC/plant/yr] -> [gC/m2/s]
site_npp = site_npp + ccohort%npp_acc_hold * ccohort%n*area_inv * &
g_per_kg * hlm_days_per_year / sec_per_day

ccohort => ccohort%shorter
end do

end if

cpatch => cpatch%younger
enddo

! Smoothed [gc/m2/yr]
if(sites(s)%ema_npp<-9000._r8)then
sites(s)%ema_npp = site_npp
else
sites(s)%ema_npp = (1._r8-1._r8/ema_npp_tscale)*sites(s)%ema_npp + (1._r8/ema_npp_tscale)*site_npp
end if

end do

return
Expand Down