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.

Sign up for GitHub

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

Arctic Carbon Cycle Updates #947

Closed
wants to merge 11 commits into from
Closed

Arctic Carbon Cycle Updates #947

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
64aa4d9
Add files via upload
lmbirch89 Mar 21, 2020
d91b2ee
Add files via upload
lmbirch89 Mar 21, 2020
8a5252e
Add files via upload
lmbirch89 Mar 21, 2020
ea984f5
Merge pull request #2 from lmbirch89/abz_update
lmbirch89 Mar 21, 2020
6432dd2
Add files via upload
lmbirch89 Mar 24, 2020
35b7098
Add files via upload
lmbirch89 Mar 24, 2020
298235a
Add files via upload
lmbirch89 Mar 31, 2020
0f0305a
Add files via upload
lmbirch89 Mar 31, 2020
61cf54f
Merge remote-tracking branch 'escomp/master' into abz_edits324_bill
billsacks Apr 15, 2020
aa0f31f
Add files via upload
lmbirch89 Apr 13, 2020
dd0cbec
Add files via upload
lmbirch89 Apr 15, 2020
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
37 changes: 30 additions & 7 deletions src/biogeochem/CNPhenologyMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -252,7 +252,7 @@ subroutine CNPhenology (bounds, num_soilc, filter_soilc, num_soilp, &
cnveg_state_inst, cnveg_carbonstate_inst, cnveg_nitrogenstate_inst, cnveg_carbonflux_inst, cnveg_nitrogenflux_inst)

call CNSeasonDecidPhenology(num_soilp, filter_soilp, &
temperature_inst, cnveg_state_inst, dgvs_inst, &
temperature_inst, waterdiagnosticbulk_inst, cnveg_state_inst, dgvs_inst, &
cnveg_carbonstate_inst, cnveg_nitrogenstate_inst, cnveg_carbonflux_inst, cnveg_nitrogenflux_inst)

call CNStressDecidPhenology(num_soilp, filter_soilp, &
Expand Down Expand Up @@ -627,7 +627,7 @@ end subroutine CNEvergreenPhenology

!-----------------------------------------------------------------------
subroutine CNSeasonDecidPhenology (num_soilp, filter_soilp , &
temperature_inst, cnveg_state_inst, dgvs_inst , &
temperature_inst, waterdiagnosticbulk_inst, cnveg_state_inst, dgvs_inst , &
cnveg_carbonstate_inst, cnveg_nitrogenstate_inst, cnveg_carbonflux_inst, cnveg_nitrogenflux_inst)
!
! !DESCRIPTION:
Expand All @@ -644,6 +644,7 @@ subroutine CNSeasonDecidPhenology (num_soilp, filter_soilp , &
integer , intent(in) :: num_soilp ! number of soil patches in filter
integer , intent(in) :: filter_soilp(:) ! filter for soil patches
type(temperature_type) , intent(in) :: temperature_inst
type(waterdiagnosticbulk_type) , intent(in) :: waterdiagnosticbulk_inst
type(cnveg_state_type) , intent(inout) :: cnveg_state_inst
type(dgvs_type) , intent(inout) :: dgvs_inst
type(cnveg_carbonstate_type) , intent(inout) :: cnveg_carbonstate_inst
Expand All @@ -656,6 +657,8 @@ subroutine CNSeasonDecidPhenology (num_soilp, filter_soilp , &
integer :: fp !lake filter patch index
real(r8):: ws_flag !winter-summer solstice flag (0 or 1)
real(r8):: crit_onset_gdd !critical onset growing degree-day sum
real(r8):: crit_daylat !latitudinal light gradient in arctic-boreal
real(r8):: onset_thresh !flag onset threshold
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks like you're using this effectively as a logical value. Please change it to a logical, such as:

logical :: do_onset  ! flag for whether to do onset now

(I realize that some other flags in this routine are real-valued like this... it's possible that's needed for some other reason, or possible that that shouldn't have been done that way in the past. In either case, let's keep new code cleaner in this respect.)

real(r8):: soilt
!-----------------------------------------------------------------------

Expand All @@ -668,7 +671,10 @@ subroutine CNSeasonDecidPhenology (num_soilp, filter_soilp , &
season_decid => pftcon%season_decid , & ! Input: binary flag for seasonal-deciduous leaf habit (0 or 1)

t_soisno => temperature_inst%t_soisno_col , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) (-nlevsno+1:nlevgrnd)

soila10 => temperature_inst%soila10_patch , & ! Input: [real(r8) (:) ]
t_a5min => temperature_inst%t_a5min_patch , & ! input: [real(r8) (:) ]
snow_5day => waterdiagnosticbulk_inst%snow_5day_col , & ! input: [real(r8) (:) ]

pftmayexist => dgvs_inst%pftmayexist_patch , & ! Output: [logical (:) ] exclude seasonal decid patches from tropics

annavg_t2m => cnveg_state_inst%annavg_t2m_patch , & ! Input: [real(r8) (:) ] annual average 2m air temperature (K)
Expand Down Expand Up @@ -742,6 +748,8 @@ subroutine CNSeasonDecidPhenology (num_soilp, filter_soilp , &
)

! start patch loop


do fp = 1,num_soilp
p = filter_soilp(fp)
c = patch%column(p)
Expand Down Expand Up @@ -838,7 +846,7 @@ subroutine CNSeasonDecidPhenology (num_soilp, filter_soilp , &

! test for switching from dormant period to growth period
if (dormant_flag(p) == 1.0_r8) then

onset_thresh = 0.0_r8
! Test to turn on growing degree-day sum, if off.
! switch on the growing degree day sum on the winter solstice

Expand All @@ -865,13 +873,21 @@ subroutine CNSeasonDecidPhenology (num_soilp, filter_soilp , &
if (onset_gddflag(p) == 1.0_r8 .and. soilt > SHR_CONST_TKFRZ) then
onset_gdd(p) = onset_gdd(p) + (soilt-SHR_CONST_TKFRZ)*fracday
end if

!separate into Arctic boreal and lower latitudes
billsacks marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ekluzek & @billsacks What's the best way to wrap this into a namelist flag to maintain backwards compatibility with CLM5 but also make sure the code is readable?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  • should we have this hardcoded behavior for latitude or define this by PFT?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  • the second part of this conditional with onset_gddflag(p) == 1.0_r8 seems to be asking for something different, but as implemented it also seems to be triggered in lower latitudes too, leading to higher than intended LAI and GPP

if (onset_gdd(p) > crit_onset_gdd .and. abs(grc%latdeg(g))<45.0_r8) then
onset_thresh=1.0_r8
else if (onset_gddflag(p) == 1.0_r8 .and. soila10(p) > SHR_CONST_TKFRZ .and. &
t_a5min(p) > SHR_CONST_TKFRZ .and. ws_flag==1.0_r8 .and. &
dayl(g)>(crit_dayl/2.0_r8) .and. snow_5day(c)<0.1_r8) then
onset_thresh=1.0_r8
end if
! set onset_flag if critical growing degree-day sum is exceeded
if (onset_gdd(p) > crit_onset_gdd) then
if (onset_thresh == 1.0_r8) then
onset_flag(p) = 1.0_r8
dormant_flag(p) = 0.0_r8
onset_gddflag(p) = 0.0_r8
onset_gdd(p) = 0.0_r8
onset_thresh = 0.0_r8
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With my above request to initialize onset_thresh (renamed to do_onset) a bit above, you can delete this line.

onset_counter(p) = ndays_on * secspday

! move all the storage pools into transfer pools,
Expand Down Expand Up @@ -913,8 +929,15 @@ subroutine CNSeasonDecidPhenology (num_soilp, filter_soilp , &
if (days_active(p) > 355._r8) pftmayexist(p) = .false.
end if

! use 15 hr (54000 min) at ~65N from eitel 2019, to ~11hours in temperate regions
! 15hr-11hr/(65N-45N)=linear slope = 720 min/latitude
crit_daylat=54000-720*(65-abs(grc%latdeg(g)))
if (crit_daylat < crit_dayl) then
crit_daylat = crit_dayl !maintain previous offset from White 2001 as minimum
end if

! only begin to test for offset daylength once past the summer sol
if (ws_flag == 0._r8 .and. dayl(g) < crit_dayl) then
if (ws_flag == 0._r8 .and. dayl(g) < crit_daylat) then
offset_flag(p) = 1._r8
offset_counter(p) = ndays_off * secspday
prev_leafc_to_litter(p) = 0._r8
Expand Down
40 changes: 22 additions & 18 deletions src/biogeophys/LunaMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -306,7 +306,9 @@ subroutine Update_Photosynthesis_Capacity(bounds, fn, filterp, &
vcmx25_z => photosyns_inst%vcmx25_z_patch , & ! Output: [real(r8) (:,:) ] patch leaf Vc,max25 (umol/m2 leaf/s) for canopy layer
jmx25_z => photosyns_inst%jmx25_z_patch , & ! Output: [real(r8) (:,:) ] patch leaf Jmax25 (umol electron/m**2/s) for canopy layer
pnlc_z => photosyns_inst%pnlc_z_patch , & ! Output: [real(r8) (:,:) ] patch proportion of leaf nitrogen allocated for light capture for canopy layer
enzs_z => photosyns_inst%enzs_z_patch & ! Output: [real(r8) (:,:) ] enzyme decay status 1.0-fully active; 0-all decayed during stress
enzs_z => photosyns_inst%enzs_z_patch , & ! Output: [real(r8) (:,:) ] enzyme decay status 1.0-fully active; 0-all decayed during stress
vcmx_prevyr => photosyns_inst%vcmx_prevyr , & ! Output: [real(r8) (:,:) ] patch leaf Vc,max25 from previous year avg
jmx_prevyr => photosyns_inst%jmx_prevyr & ! Output: [real(r8) (:,:) ] patch leaf Jmax25 from previous year avg
)
!----------------------------------------------------------------------------------------------------------------------------------------------------------
!set timestep
Expand All @@ -332,7 +334,7 @@ subroutine Update_Photosynthesis_Capacity(bounds, fn, filterp, &
hourpd = dayl(g) / 3600._r8
tleafd10 = t_veg10_day(p) - tfrz
tleafn10 = t_veg10_night(p) - tfrz
tleaf10 = (dayl(g)*tleafd10 +(86400._r8-dayl(g)) * tleafd10)/86400._r8
tleaf10 = (dayl(g)*tleafd10 +(86400._r8-dayl(g)) * tleafn10)/86400._r8
ekluzek marked this conversation as resolved.
Show resolved Hide resolved
tair10 = t10(p)- tfrz
relh10 = min(1.0_r8, rh10_p(p))
rb10v = rb10_p(p)
Expand Down Expand Up @@ -402,18 +404,20 @@ subroutine Update_Photosynthesis_Capacity(bounds, fn, filterp, &
PNcbold = 0.0_r8
call NitrogenAllocation(FNCa,forc_pbot10(p), relh10, CO2a10, O2a10, PARi10, PARimx10, rb10v, hourpd, &
tair10, tleafd10, tleafn10, &
Jmaxb0, Jmaxb1, Wc2Wjb0, relhExp, PNlcold, PNetold, PNrespold, &
PNcbold, PNstoreopt, PNlcopt, PNetopt, PNrespopt, PNcbopt)
Jmaxb0, Jmaxb1, Wc2Wjb0, relhExp, PNlcold, PNetold, PNrespold, PNcbold, &
dayl_factor(p), PNstoreopt, PNlcopt, PNetopt, PNrespopt, PNcbopt)
vcmx25_opt= PNcbopt * FNCa * Fc25
jmx25_opt= PNetopt * FNCa * Fj25

chg = vcmx25_opt-vcmx25_z(p, z)
chg_constrn = min(abs(chg),vcmx25_z(p, z)*max_daily_pchg)
vcmx25_z(p, z) = vcmx25_z(p, z)+sign(1.0_r8,chg)*chg_constrn

vcmx_prevyr(p,z) = vcmx25_z(p,z)
ekluzek marked this conversation as resolved.
Show resolved Hide resolved

chg = jmx25_opt-jmx25_z(p, z)
chg_constrn = min(abs(chg),jmx25_z(p, z)*max_daily_pchg)
jmx25_z(p, z) = jmx25_z(p, z)+sign(1.0_r8,chg)*chg_constrn
jmx_prevyr(p,z) = jmx25_z(p,z)

PNlc_z(p, z)= PNlcopt

Expand Down Expand Up @@ -472,8 +476,8 @@ subroutine Update_Photosynthesis_Capacity(bounds, fn, filterp, &
endif !if not C3 plants
else
do z = 1 , nrad(p)
jmx25_z(p, z) = 85._r8
vcmx25_z(p, z) = 50._r8
jmx25_z(p, z) = jmx_prevyr(p,z)
vcmx25_z(p, z) = vcmx_prevyr(p,z)
ekluzek marked this conversation as resolved.
Show resolved Hide resolved
end do
endif !checking for LAI and LNC
endif !the first daycheck
Expand Down Expand Up @@ -791,7 +795,7 @@ end subroutine Clear24_Climate_LUNA
!Use the LUNA model to calculate the Nitrogen partioning
subroutine NitrogenAllocation(FNCa,forc_pbot10, relh10, CO2a10,O2a10, PARi10,PARimx10,rb10, hourpd, tair10, tleafd10, tleafn10, &
Jmaxb0, Jmaxb1, Wc2Wjb0, relhExp,&
PNlcold, PNetold, PNrespold, PNcbold, &
PNlcold, PNetold, PNrespold, PNcbold, dayl_factor,&
PNstoreopt, PNlcopt, PNetopt, PNrespopt, PNcbopt)
implicit none
real(r8), intent (in) :: FNCa !Area based functional nitrogen content (g N/m2 leaf)
Expand All @@ -814,12 +818,12 @@ subroutine NitrogenAllocation(FNCa,forc_pbot10, relh10, CO2a10,O2a10, PARi10,PAR
real(r8), intent (in) :: PNetold !old value of the proportion of nitrogen allocated to electron transport (unitless)
real(r8), intent (in) :: PNrespold !old value of the proportion of nitrogen allocated to respiration (unitless)
real(r8), intent (in) :: PNcbold !old value of the proportion of nitrogen allocated to carboxylation (unitless)
real(r8), intent (in) :: dayl_factor !daylight scale factor
real(r8), intent (out):: PNstoreopt !optimal proportion of nitrogen for storage
real(r8), intent (out):: PNlcopt !optimal proportion of nitrogen for light capture
real(r8), intent (out):: PNetopt !optimal proportion of nitrogen for electron transport
real(r8), intent (out):: PNrespopt !optimal proportion of nitrogen for respiration
real(r8), intent (out):: PNcbopt !optial proportion of nitrogen for carboxyaltion

!-------------------------------------------------------------------------------------------------------------------------------
!intermediate variables
real(r8) :: Carboncost1 !absolute amount of carbon cost associated with maintenance respiration due to deccrease in light capture nitrogen(g dry mass per day)
Expand Down Expand Up @@ -896,12 +900,12 @@ subroutine NitrogenAllocation(FNCa,forc_pbot10, relh10, CO2a10,O2a10, PARi10,PAR
jj = 1
tleafd10c = min(max(tleafd10, Trange1), Trange2) !constrain the physiological range
tleafn10c = min(max(tleafn10, Trange1), Trange2) !constrain the physiological range
ci = 0.7_r8 * CO2a10
JmaxCoef = Jmaxb1 * ((hourpd / 12.0_r8)**2.0_r8) * (1.0_r8 - exp(-relhExp * max(relh10 - minrelh, 0.0_r8) / &
ci = 0.7_r8 * CO2a10
JmaxCoef = Jmaxb1 * dayl_factor * (1.0_r8 - exp(-relhExp * max(relh10 - minrelh, 0.0_r8) / &
ekluzek marked this conversation as resolved.
Show resolved Hide resolved
(1.0_r8 - minrelh)))
do while (PNlcoldi .NE. PNlc .and. jj < 100)
Fc = VcmxTKattge(tair10, tleafd10c) * Fc25
Fj = JmxTKattge(tair10, tleafd10c) * Fj25
Fc = VcmxTLeuning(tair10, tleafd10c) * Fc25
ekluzek marked this conversation as resolved.
Show resolved Hide resolved
Fj = JmxTLeuning(tair10, tleafd10c) * Fj25
NUEr = Cv * NUEr25 * (RespTBernacchi(tleafd10c) * hourpd + RespTBernacchi(tleafn10c) * (24.0_r8 - hourpd)) !nitrogen use efficiency for respiration (g biomass/m2/day/g N)
!****************************************************
!Nitrogen Allocation Scheme: store the initial value
Expand Down Expand Up @@ -1054,7 +1058,7 @@ subroutine Nitrogen_investments (KcKjFlag, FNCa, Nlc, forc_pbot10, relh10, &
A = (1.0_r8 - theta_cj) * max(Wc, Wj) + theta_cj * min(Wc, Wj)
endif
PSN = Cv * A * hourpd
Vcmaxnight = VcmxTKattge(tair10, tleafn10) / VcmxTKattge(tair10, tleafd10) * Vcmax
Vcmaxnight = VcmxTLeuning(tair10, tleafn10) / VcmxTLeuning(tair10, tleafd10) * Vcmax
RESP = Cv * leaf_mr_vcm * (Vcmax * hourpd + Vcmaxnight * (24.0_r8 - hourpd))
Net = Jmax / Fj
Ncb = Vcmax / Fc
Expand Down Expand Up @@ -1209,8 +1213,8 @@ subroutine NUEref(NUEjref,NUEcref,Kj2Kcref)

tgrow = 25.0_r8
tleaf = 25.0_r8
Fc = VcmxTKattge(tgrow, tleaf) * Fc25
Fj = JmxTKattge(tgrow, tleaf) * Fj25
Fc = VcmxTLeuning(tgrow, tleaf) * Fc25
Fj = JmxTLeuning(tgrow, tleaf) * Fj25
CO2c = co2ref * forc_pbot_ref * 1.0e-6_r8 !pa
O2c = O2ref * forc_pbot_ref * 1.0e-6_r8 !pa
k_c = params_inst%kc25_coef * exp((79430.0_r8 / (rgas*1.e-3_r8 * (25.0_r8 + tfrz))) * (1.0_r8 - (tfrz + 25.0_r8) / (tfrz + tleaf)))
Expand Down Expand Up @@ -1250,8 +1254,8 @@ subroutine NUE(O2a, ci, tgrow, tleaf, NUEj,NUEc,Kj2Kc)
real(r8) :: awc !second deminator term for rubsico limited carboxylation rate based on Farquhar model
real(r8) :: c_p !CO2 compenstation point (Pa)

Fc = VcmxTKattge(tgrow, tleaf) * Fc25
Fj = JmxTKattge(tgrow, tleaf) * Fj25
Fc = VcmxTLeuning(tgrow, tleaf) * Fc25
Fj = JmxTLeuning(tgrow, tleaf) * Fj25
k_c = params_inst%kc25_coef * exp((79430.0_r8 / (rgas*1.e-3_r8 * (25.0_r8 + tfrz))) * (1.0_r8 - (tfrz + 25.0_r8) / (tfrz + tleaf)))
k_o = params_inst%ko25_coef * exp((36380.0_r8 / (rgas*1.e-3_r8 * (25.0_r8 + tfrz))) * (1.0_r8 - (tfrz + 25.0_r8) / (tfrz + tleaf)))
c_p = params_inst%cp25_yr2000 * exp((37830.0_r8 / (rgas*1.e-3_r8 * (25.0_r8 + tfrz))) * (1.0_r8 - (tfrz + 25.0_r8) / (tfrz + tleaf)))
Expand Down
Loading