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

nutrient dynamics v2 of allocation and acquisition #880

Merged
merged 61 commits into from
Dec 2, 2022
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
61 commits
Select commit Hold shift + click to select a range
751490c
Setting target_m to dummy value to overcome IBM issues with base func…
rgknox Oct 18, 2021
9563af1
NCLMax = 3
rgknox Oct 28, 2021
7191021
Added nutrient storage limit and derivative to test effects on statur…
rgknox Oct 30, 2021
949e9ca
fixed logic limiting growth
rgknox Nov 3, 2021
5084a1d
Adjusted growth restrictions and uptake downregulation so that downre…
rgknox Nov 4, 2021
760e870
First complete pass-through on having a dynamic root response to nutr…
rgknox Nov 6, 2021
6f1b1ca
Work on dynamic root responses to nutrient availability
rgknox Nov 18, 2021
e1d88b1
version 2 of fates nutrient aquisition, bug fixes and enabling dynami…
rgknox Nov 24, 2021
561750f
Updates to fates nutrient coupling and feedback algorithm. Modified p…
rgknox Dec 3, 2021
f2a9b05
Adding l2fr to the restarts, updating soil root biomass diagnostic
rgknox Jan 1, 2022
9b50744
Updates to fates history variables for NP cycling
rgknox Jan 20, 2022
04c454f
Incremental work on l2fr searching
rgknox Feb 2, 2022
9f49427
Resolved conflicts between long-lived CNP update 2 branch and API 20
rgknox Feb 3, 2022
6bd5eb1
Updated l2fr diagnostics to be canopy and understory
rgknox Feb 3, 2022
b570617
Fixed l2fr and storage fraction history fields, updated some name con…
rgknox Feb 8, 2022
7a4a95a
Incremental work towards getting dynamic root response c/n balance to…
rgknox Feb 20, 2022
12f824d
Incremental changes towards fine-root biomass adaptation to N:C and P…
rgknox Feb 24, 2022
d017a31
merging in Charlies sym fix code
ckoven Oct 29, 2021
c4422f1
various updates to the cnp acquisition algorithm, including symbiotic…
rgknox Mar 15, 2022
1e34c39
Various updates to l2fr search algorithm
rgknox Apr 1, 2022
2326459
Added some parameters that were hard coded constants to the parameter…
rgknox Apr 15, 2022
5f44bfa
Added carbon storage fraction by sizexpft ustory/canopy
rgknox Apr 15, 2022
7b1a8b7
Started work on root reabsorption for cnp optimization
rgknox Apr 15, 2022
8fc9706
Merge branch 'cnp-dynamic-root-allom-dfdd-api21-symfix' of github.com…
rgknox Apr 15, 2022
7a1178f
CNP updates, small modifications to things like root trimming and lim…
rgknox May 6, 2022
8762093
cleaning up the cnp version 2 allocation code, adding in a diagnostic…
rgknox May 24, 2022
36491b9
Adding recruit l2fr and smoothing to the f_cn function
rgknox Jul 29, 2022
a75b6d6
Adding canopy layer x size x pft history diagnostics
rgknox Aug 3, 2022
434c0e0
Added in CNP limiter logic
rgknox Aug 8, 2022
893e227
Alternative PID error terms in CNP
rgknox Aug 9, 2022
681e4f5
Testing different PID controller methods and data prep methods.
rgknox Aug 14, 2022
43ee9b9
yet more reorganization and PID algorithm testing
rgknox Aug 19, 2022
8df96ae
various code modifications while exploring PID functions
rgknox Sep 21, 2022
75a7710
Updating parameter names to include PID. Cleaning up parteh CNP code
rgknox Sep 22, 2022
18bde56
Merge CNP branch up to API 24, resolve conflicts
rgknox Sep 22, 2022
13303d1
Updated parameter format to add PID terms, remove adaptive timescale …
rgknox Sep 22, 2022
2ae674c
Updating parameters for CNP, and various code cleaning for CNP
rgknox Sep 23, 2022
5d17995
Merge resolution fixes between CNP v2 and api24
rgknox Sep 23, 2022
26ee582
merge resolutions and parameter stuff
rgknox Sep 24, 2022
6e019d2
Updates to parameter file
rgknox Sep 24, 2022
06757ed
updating patch files
rgknox Sep 24, 2022
4465d7c
More merge resolutions between cnp v2 and api24
rgknox Sep 25, 2022
d490584
Cleanup in PRTAllometricCNP
rgknox Oct 11, 2022
68d71c2
More cleanup CNP growth code, removal of unused variables, labels, un…
rgknox Oct 11, 2022
cf89717
cnp v2, syntax cleanup
rgknox Oct 14, 2022
91bd865
cnp cleanup: change p uptake to p gain for consistency, remove tempor…
rgknox Oct 17, 2022
4ee19d2
Reverting nclmax
rgknox Oct 17, 2022
b194c1f
Clean up of fates cnp restart code and testing
rgknox Oct 19, 2022
9e69a0f
Added NPP boundary condition with ELM/CLM to enable npp hypothesis on…
rgknox Nov 8, 2022
f360032
first pass at resolving conflicts between nutrient cycling v2 and cro…
rgknox Nov 9, 2022
20e1c76
fixes to symbiotic fixation
rgknox Nov 9, 2022
fa8841e
update to bci xml parameter patch
rgknox Nov 15, 2022
c1f49d2
Merge branch 'master' into cnp-dynamic-root-allom-dfdd-api21-symfix
rgknox Nov 17, 2022
fb5a9de
Updated the fates default parameter file API update xmls
rgknox Nov 28, 2022
dd6f067
fix typo on carbon12_element
rgknox Nov 29, 2022
ab7af56
bug fix in fineroot allometry
rgknox Nov 29, 2022
9f6b8e9
Removed unused variable
rgknox Nov 29, 2022
b746edd
Reverting nclmax to 2
rgknox Nov 30, 2022
89cf208
fixed integration array indices for CNP
rgknox Nov 30, 2022
d3279d3
Removed the CLSZPF dimension temporarily to pass clm tests
rgknox Dec 1, 2022
3027e30
Merge resolution between v2 nutrients and the hydro pointer bug fix
rgknox Dec 1, 2022
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
yet more reorganization and PID algorithm testing
  • Loading branch information
rgknox committed Aug 19, 2022
commit 43ee9b93ead585992235d76531ec7434f87bab64
2 changes: 0 additions & 2 deletions biogeochem/EDPhysiologyMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -2647,8 +2647,6 @@ subroutine UpdateRecruitL2FR(csite)
rec_l2fr0(ft,cl) = rec_l2fr0(ft,cl) / rec_n(ft,cl)
csite%rec_l2fr(ft,cl) = &
(1._r8-smth_wgt)*csite%rec_l2fr(ft,cl) + smth_wgt*rec_l2fr0(ft,cl)

!print*,"REC_L2FR:",cl,csite%rec_l2fr(ft,cl)
end if
end do
end do
Expand Down
155 changes: 86 additions & 69 deletions parteh/PRTAllometricCNPMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -107,8 +107,8 @@ module PRTAllometricCNPMod

integer, parameter, private :: pid_spe_controller = storage_spe

real(r8), parameter, private :: pid_int_wgt = 1._r8/30._r8 ! n-day smoothing (K on the integral of PID)
real(r8), parameter, private :: pid_drv_wgt = 1._r8/30._r8 ! n-day smoothing (K on the derivative of PID)
real(r8), parameter, private :: pid_int_wgt = 1._r8/10._r8 ! n-day smoothing (K on the integral of PID)
real(r8), parameter, private :: pid_drv_wgt = 1._r8/10._r8 ! n-day smoothing (K on the derivative of PID)

! Global identifiers for the two stoichiometry values
integer,public, parameter :: stoich_growth_min = 1 ! Flag for stoichiometry associated with
Expand Down Expand Up @@ -481,7 +481,7 @@ subroutine DailyPRTAllometricCNP(this,co_num,nplant)

! This routine updates the l2fr (leaf 2 fine-root multiplier) variable
! It will also update the target
call this%CNPAdjustFRootTargets(target_c,target_dcdd,co_num,nplant)
!call this%CNPAdjustFRootTargets(target_c,target_dcdd,co_num,nplant)

! Remember the original C,N,P states to help with final
! evaluation of how much was allocated
Expand Down Expand Up @@ -582,6 +582,11 @@ subroutine DailyPRTAllometricCNP(this,co_num,nplant)
end do
call endrun(msg=errMsg(sourcefile, __LINE__))
end if


! This routine updates the l2fr (leaf 2 fine-root multiplier) variable
! It will also update the target
!call this%CNPAdjustFRootTargets(target_c,target_dcdd,co_num,nplant)

! ===================================================================================
! Step 3.
Expand All @@ -591,7 +596,7 @@ subroutine DailyPRTAllometricCNP(this,co_num,nplant)

call this%CNPAllocateRemainder(c_gain, n_gain, p_gain, &
c_gain0, n_gain0, p_gain0, &
c_efflux, n_efflux, p_efflux)
c_efflux, n_efflux, p_efflux,co_num,nplant,target_c,target_dcdd)


if(n_uptake_mode.ne.prescribed_n_uptake) then
Expand Down Expand Up @@ -730,6 +735,7 @@ subroutine CNPAdjustFRootTargets(this, target_c, target_dcdd,co_num,nplant)
real(r8) :: logi_k
real(r8) :: l2fr_mult
real(r8) :: l2fr_delta
real(r8) :: nc_ratio, pc_ratio
real(r8) :: dxcdt_ratio ! log change (derivative) of the maximum of the N/C and P/C storage ratio
real(r8) :: xc_ratio ! log Maximum of the N/C and P/C storage ratio
real(r8), pointer :: ema_xc ! The exponential moving average of the N-or-P versus C PID error function
Expand All @@ -739,6 +745,14 @@ subroutine CNPAdjustFRootTargets(this, target_c, target_dcdd,co_num,nplant)
real(r8), parameter :: max_l2fr_cgain_frac = 0.99_r8
real(r8), parameter :: xc_ratio_correction = 1.0_r8

integer, parameter :: pid_c_function = 0
integer, parameter :: pid_n_function = 1
integer, parameter :: pid_minnc_function = 2
integer, parameter :: pid_alogmaxnc_function = 3
integer, parameter :: pid_ncratio_function = 4

!integer, parameter :: pid_function = pid_c_function

! This is the relative scaling strength of the derivative term (K_d),
! compared to the combined proportion and integral term (K_p and K_i)
! Note the strength of the derivative is should be about half of the period
Expand Down Expand Up @@ -767,8 +781,8 @@ subroutine CNPAdjustFRootTargets(this, target_c, target_dcdd,co_num,nplant)

store_c_max = target_c(store_organ)

store_c_act = this%GetState(store_organ, carbon12_element) + &
this%bc_in(acnp_bc_in_id_netdc)%rval
store_c_act = max(0.001_r8*store_c_max,this%GetState(store_organ, carbon12_element) + &
this%bc_in(acnp_bc_in_id_netdc)%rval)

if(n_uptake_mode.ne.prescribed_n_uptake)then

Expand All @@ -780,15 +794,28 @@ subroutine CNPAdjustFRootTargets(this, target_c, target_dcdd,co_num,nplant)
store_nut_act = this%GetState(store_organ, nitrogen_element) + &
this%bc_inout(acnp_bc_inout_id_netdn)%rval

!n_ratio = xc_ratio_correction*min(50.0_r8,max(0.02_r8,(store_nut_act/store_nut_max)/(store_c_act/store_c_max)))
!n_ratio = xc_ratio_correction*min(50.0_r8,max(0.02_r8,(store_nut_act/store_nut_max)))
!n_ratio = xc_ratio_correction*min(50.0_r8,max(0.02_r8,1._r8/(store_c_act/store_c_max)))

if((store_nut_act/store_nut_max) > (store_c_act/store_c_max))then
n_ratio = xc_ratio_correction*min(50.0_r8,max(0.02_r8,(store_c_max/store_c_act)))
else
n_ratio = xc_ratio_correction*min(50.0_r8,max(0.02_r8,(store_nut_act/store_nut_max)))
end if
select case(nint(EDPftvarcon_inst%dev_arbitrary_pft(ipft)))
case(pid_c_function)
n_ratio = store_c_max/store_c_act
case(pid_n_function)
n_ratio = store_nut_act/store_nut_max
case(pid_minnc_function)
if((store_nut_act/store_nut_max) > (store_c_act/store_c_max))then
n_ratio = (store_c_max/store_c_act)
else
n_ratio = (store_nut_act/store_nut_max)
end if
case(pid_alogmaxnc_function)
if( abs(log(store_nut_act/store_nut_max)) < abs(log(store_c_act/store_c_max))) then
n_ratio = (store_c_max/store_c_act)
else
n_ratio = (store_nut_act/store_nut_max)
end if
case(pid_ncratio_function)
n_ratio = (store_nut_act/store_nut_max)/(store_c_act/store_c_max)
end select

nc_ratio = (store_nut_act/store_nut_max)/(store_c_act/store_c_max)

else
n_ratio = -1._r8
Expand All @@ -804,16 +831,28 @@ subroutine CNPAdjustFRootTargets(this, target_c, target_dcdd,co_num,nplant)
store_nut_act = this%GetState(store_organ, phosphorus_element) + &
this%bc_inout(acnp_bc_inout_id_netdp)%rval

!p_ratio = xc_ratio_correction*min(50.0_r8,max(0.02_r8,(store_nut_act/store_nut_max)/(store_c_act/store_c_max)))
p_ratio = xc_ratio_correction*min(50.0_r8,max(0.02_r8,(store_nut_act/store_nut_max)))
!p_ratio = xc_ratio_correction*min(50.0_r8,max(0.02_r8,1._r8/(store_c_act/store_c_max)))

if((store_nut_act/store_nut_max) > (store_c_act/store_c_max))then
p_ratio = xc_ratio_correction*min(50.0_r8,max(0.02_r8,(store_c_max/store_c_act)))
else
p_ratio = xc_ratio_correction*min(50.0_r8,max(0.02_r8,(store_nut_act/store_nut_max)))
end if

select case(nint(EDPftvarcon_inst%dev_arbitrary_pft(ipft)))
case(pid_c_function)
p_ratio = store_c_max/store_c_act
case(pid_n_function)
p_ratio = store_nut_act/store_nut_max
case(pid_minnc_function)
if((store_nut_act/store_nut_max) > (store_c_act/store_c_max))then
p_ratio = (store_c_max/store_c_act)
else
p_ratio = (store_nut_act/store_nut_max)
end if
case(pid_alogmaxnc_function)
if( abs(log(store_nut_act/store_nut_max)) < abs(log(store_c_act/store_c_max))) then
p_ratio = (store_c_max/store_c_act)
else
p_ratio = (store_nut_act/store_nut_max)
end if
case(pid_ncratio_function)
p_ratio = (store_nut_act/store_nut_max)/(store_c_act/store_c_max)
end select

pc_ratio = (store_nut_act/store_nut_max)/(store_c_act/store_c_max)

else
p_ratio = -1._r8
Expand Down Expand Up @@ -875,34 +914,7 @@ subroutine CNPAdjustFRootTargets(this, target_c, target_dcdd,co_num,nplant)
max(0._r8,target_c(struct_organ)-struct_c) - &
max(0._r8,target_c(store_organ)-store_c))

if(pid_method==pid_ratio) then

l2fr_delta_max = (c_fnrt_expand + target_c(fnrt_organ))/target_c(fnrt_organ)

! This value could be negative if there is no gain, or less gain
! than what can replace tissues, just ensure the multiplier is GT 1

l2fr_delta_max = max(1._r8,l2fr_delta_max)

! Determine the max change for the doubling timescale
! 2.0 = l2fr_delta_max^frnt_adapt_tscl

l2fr_delta_scale = 2._r8**(1._r8/prt_params%fnrt_adapt_tscale(ipft))-1.0_r8

! Calculate the un-regulated l2fr multiplier

logi_k = EDPftvarcon_inst%dev_arbitrary_pft(ipft)
l2fr_mult = l2fr_delta_scale*(2.0_r8/(1.0_r8 + exp(ema_xc)**logi_k)-1.0_r8)+1.0_r8


if(l2fr_mult>1.0_r8)then
l2fr_mult = min(l2fr_mult,l2fr_delta_max)
end if

l2fr = max(l2fr_min,l2fr * l2fr_mult)


elseif(pid_method==pid_logratio) then
if(pid_method==pid_logratio) then

! When using a log based (additive) PID search method
! we define 1/fnrt_adapt_tscale as delta we want
Expand All @@ -914,48 +926,49 @@ subroutine CNPAdjustFRootTargets(this, target_c, target_dcdd,co_num,nplant)

l2fr_delta_scale = (1._r8/prt_params%fnrt_adapt_tscale(ipft))/log(2.0_r8)

! log(2.0_r8)*l2fr_delta_scale = (1._r8/prt_params%fnrt_adapt_tscale(ipft))

! ema_xc is already in log form to allow for averaging

! Want the derivative to be strongest when storage is most disproportionate
l2fr_deriv_scale = 0.0_r8 !0.25*abs(ema_xc/ema_dxcdt)
l2fr_deriv_scale = 0.0_r8 !-20.0_r8 !0.25*abs(ema_xc/ema_dxcdt)
l2fr_int_scale = 0.0_r8


! To limit overshoot, when either positive and decending or negative and ascending
! we have already corrected enough to change the behavior, lets
! decrease the scaling

if( ((ema_xc > 0._r8) .and. (ema_dxcdt<0)) .or. ((ema_xc < 0._r8) .and. (ema_dxcdt>0))) then
l2fr_delta_scale = 0.05_r8 * l2fr_delta_scale
if( ((ema_xc > 0._r8) .and. (ema_dxcdt<0._r8)) .or. ((ema_xc < 0._r8) .and. (ema_dxcdt>0._r8))) then
l2fr_delta_scale = 0.1_r8 * l2fr_delta_scale
sup_flag = 1
else
sup_flag = 0

end if


!l2fr_delta = -l2fr_delta_scale*(min(0.2,max(-0.2,xc_ratio)) + ema_xc*l2fr_int_scale + ema_dxcdt*l2fr_deriv_scale)
l2fr_delta = -l2fr_delta_scale*(xc_ratio + ema_xc*l2fr_int_scale + ema_dxcdt*l2fr_deriv_scale)

! Cap growth and shrinkage to avoid large changes
! (currently capping at projected rate for a 2:1 ratio
l2fr_delta_minmax = l2fr_delta_scale*log(2.0)
l2fr_delta_minmax = l2fr_delta_scale*log(20.0)

! Don't allow more growth than we have carbon to pay for
l2fr_delta_max = min(l2fr_delta_minmax,l2fr*(c_fnrt_expand + target_c(fnrt_organ))/target_c(fnrt_organ))
!! l2fr_delta_max = min(l2fr_delta_minmax,l2fr*(c_fnrt_expand + target_c(fnrt_organ))/target_c(fnrt_organ))


l2fr_delta = max(-l2fr_delta_minmax,min(l2fr_delta,l2fr_delta_max))
!! l2fr_delta = max(-l2fr_delta_minmax,min(l2fr_delta,l2fr_delta_max))
!!l2fr_delta = max(-l2fr_delta_minmax,min(l2fr_delta,l2fr_delta_minmax))


! Apply the delta, also, avoid generating incredibly small l2fr's,
! super small l2frs will occur in plants that perpetually get almost
! now carbon gain, such as newly recruited plants in a dark understory

l2fr = max(l2fr_min, l2fr + l2fr_delta )

!if(co_num==1)
print*,'AAX1',co_num,hlm_current_year,hlm_day_of_year,dbh,nplant,sup_flag,xc_ratio,l2fr

if(co_num==1) print*,'AAX1',co_num,hlm_current_year,hlm_day_of_year,dbh,nplant,sup_flag,log(nc_ratio),l2fr

else

Expand Down Expand Up @@ -1906,7 +1919,7 @@ end subroutine CNPStatureGrowth

subroutine CNPAllocateRemainder(this, c_gain,n_gain,p_gain, &
c_gain0, n_gain0, p_gain0, &
c_efflux, n_efflux, p_efflux)
c_efflux, n_efflux, p_efflux,co_num,nplant,target_c,target_dcdd)

class(cnp_allom_prt_vartypes) :: this
real(r8), intent(inout) :: c_gain
Expand All @@ -1918,7 +1931,11 @@ subroutine CNPAllocateRemainder(this, c_gain,n_gain,p_gain, &
real(r8), intent(inout) :: c_efflux
real(r8), intent(inout) :: n_efflux
real(r8), intent(inout) :: p_efflux

integer,intent(in) :: co_num
real(r8),intent(in) :: nplant
real(r8) :: target_c(:)
real(r8) :: target_dcdd(:)

integer :: i
real(r8), dimension(num_organs) :: deficit_n
real(r8), dimension(num_organs) :: deficit_p
Expand Down Expand Up @@ -2040,7 +2057,7 @@ subroutine CNPAllocateRemainder(this, c_gain,n_gain,p_gain, &
deficit_p(i) = max(0._r8,this%GetDeficit(phosphorus_element,l2g_organ_list(i),target_p))

end do

! -----------------------------------------------------------------------------------
! Nutrient Fluxes proportionally to each pool (these should be fully actualized)
! (this also removes from the gain pools)
Expand All @@ -2054,8 +2071,8 @@ subroutine CNPAllocateRemainder(this, c_gain,n_gain,p_gain, &
call ProportionalNutrAllocation(this,deficit_p(1:num_organs), &
p_gain, phosphorus_element, l2g_organ_list(1:num_organs))



call this%CNPAdjustFRootTargets(target_c,target_dcdd,co_num,nplant)


! -----------------------------------------------------------------------------------
Expand Down