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

Fix negative leaf intracellular CO2 computed in PhotosynthesisMod #1272

Merged
merged 6 commits into from
Apr 25, 2017
+81 −89
Merged

Fix negative leaf intracellular CO2 computed in PhotosynthesisMod #1272

Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
34e30d0
Bug fix for the iterative solver of stomata conductance
jinyuntang Jan 23, 2017
7a0ec08
Reset ag for negative an in photosynthesis calculation
jinyuntang Jan 24, 2017
286e51b
limit ci re-calculation to minimum stomata conductance
jinyuntang Jan 25, 2017
c1e586b
clean up unnecssary comments
jinyuntang Apr 19, 2017
d3b8fd7
replace 1.e-14_r8 with tiny_val
jinyuntang Apr 20, 2017
28d854c
remove some commented lines
jinyuntang Apr 20, 2017
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
170 changes: 81 additions & 89 deletions components/clm/src/biogeophys/PhotosynthesisMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,6 @@ module PhotosynthesisMod
use clm_varctl , only : cnallocate_carbon_only
use clm_varctl , only : cnallocate_carbonnitrogen_only
use clm_varctl , only : cnallocate_carbonphosphorus_only
use clm_varctl , only : iulog
use pftvarcon , only : noveg

!
Expand Down Expand Up @@ -209,6 +208,7 @@ subroutine Photosynthesis ( bounds, fn, filterp, &
real(r8) :: lpc(bounds%begp:bounds%endp) ! leaf P concentration (gP leaf/m^2)
real(r8) :: sum_nscaler
real(r8) :: total_lai
real(r8),parameter :: tiny_val=1.e-14_r8
!------------------------------------------------------------------------------

! Temperature and soil water response functions
Expand Down Expand Up @@ -728,14 +728,26 @@ subroutine Photosynthesis ( bounds, fn, filterp, &

! End of ci iteration. Check for an < 0, in which case gs_mol = bbb

if (an(p,iv) < 0._r8) gs_mol(p,iv) = bbb(p)
! if (an(p,iv) < 0._r8) gs_mol(p,iv) = bbb(p)
! Brutely force an to zero if gs_mol is at its minimal value
if(abs(gs_mol(p,iv)-bbb(p))<tiny_val)then
an(p,iv)=0._r8
ag(p,iv)=lmr_z(p,iv)
endif

! Final estimates for cs and ci (needed for early exit of ci iteration when an < 0)

cs = cair(p) - 1.4_r8/gb_mol(p) * an(p,iv) * forc_pbot(c)
cs = max(cs,1.e-06_r8)
ci_z(p,iv) = cair(p) - an(p,iv) * forc_pbot(c) * (1.4_r8*gs_mol(p,iv)+1.6_r8*gb_mol(p)) / (gb_mol(p)*gs_mol(p,iv))

if(abs(gs_mol(p,iv)-bbb(p))<tiny_val)then
ci_z(p,iv) = cair(p) - an(p,iv) * forc_pbot(c) * (1.4_r8*gs_mol(p,iv)+1.6_r8*gb_mol(p)) / (gb_mol(p)*gs_mol(p,iv))
else
ci_z(p,iv) = ciold
endif
if(ci_z(p,iv)<0._r8)then
write(iulog,*)'negative ci',ci_z(p,iv),an(p,iv),gs_mol(p,iv)-bbb(p)
call endrun(decomp_index=p, clmlevel=namep, msg=errmsg(__FILE__, __LINE__))
endif
! Convert gs_mol (umol H2O/m**2/s) to gs (m/s) and then to rs (s/m)

gs = gs_mol(p,iv) / cf
Expand Down Expand Up @@ -1030,82 +1042,62 @@ subroutine hybrid(x0, p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z,&
real(r8) :: fa, fb
real(r8) :: x1, f0, f1
real(r8) :: x, dx
real(r8), parameter :: eps = 1.e-2_r8 !relative accuracy
real(r8), parameter :: eps1= 1.e-4_r8
integer, parameter :: itmax = 40 !maximum number of iterations
real(r8) :: tol,minx,minf

real(r8) :: ci_val(5)
real(r8) :: fi_val(5)
integer :: ii, mi
real(r8) :: tiny_val=1.e-14_r8
associate(&
c3flag => photosyns_vars%c3flag_patch , & ! Output: [logical (:) ] true if C3 and false if C4
cp => photosyns_vars%cp_patch & !Intput: [real(r8) (:) ] CO2 compensation point (Pa)
)
iter=0
call ci_func(x0, f0, p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z, rh_can, gs_mol, &
atm2lnd_vars, photosyns_vars)

if(f0 == 0._r8)return

minx=x0
minf=f0
x1 = x0 * 0.99_r8

call ci_func(x1,f1, p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z, rh_can, gs_mol, &
atm2lnd_vars, photosyns_vars)

if(f1==0._r8)then
x0 = x1
return
endif
if(f1<minf)then
minx=x1
minf=f1
if(abs(f0) < tiny_val)return
ci_val(3)=x0
fi_val(3)=f0
!compute the minimum ci value
if(c3flag(p))then
ci_val(1)=cp(p)+1.e-6_r8
else
ci_val(1)=1.e-6_r8
endif
call ci_func(ci_val(1), fi_val(1), p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z, rh_can, gs_mol, &
atm2lnd_vars, photosyns_vars)

!first use the secant approach, then use the brent approach as a backup
iter = 0
do
iter = iter + 1
dx = - f1 * (x1-x0)/(f1-f0)
x = x1 + dx
tol = abs(x) * eps
if(abs(dx)<tol)then
x0 = x
exit
endif
x0 = x1
f0 = f1
x1 = x

call ci_func(x1,f1, p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z, rh_can, gs_mol, &
atm2lnd_vars, photosyns_vars)
ci_val(2)=(ci_val(1)+ci_val(3))*0.5
call ci_func(ci_val(2), fi_val(2), p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z, rh_can, gs_mol, &
atm2lnd_vars, photosyns_vars)

if(f1<minf)then
minx=x1
minf=f1
endif
if(abs(f1)<=eps1)then
x0 = x1
exit
endif
ci_val(4)=(cair+ci_val(3))*0.5
call ci_func(ci_val(4), fi_val(4), p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z, rh_can, gs_mol, &
atm2lnd_vars, photosyns_vars)

!if a root zone is found, use the brent method for a robust backup strategy
if(f1 * f0 < 0._r8)then
!compute the maximum ci value
ci_val(5)=cair*0.999_r8
call ci_func(ci_val(5), fi_val(5), p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z, rh_can, gs_mol, &
atm2lnd_vars, photosyns_vars)

call brent(x, x0,x1,f0,f1, tol, p, iv, c, gb_mol, je, cair, oair, &
mi = -1
do ii = 1, 4
if(fi_val(ii)*fi_val(ii+1)<0._r8)then
mi = ii
endif
enddo
if(mi > 0)then
x0 = ci_val(mi)
f0 = fi_val(mi)
x1 = ci_val(mi+1)
f1 = ci_val(mi+1)
call brent(x, x0,x1,f0,f1, tol, p, iv, c, gb_mol, je, cair, oair, &
lmr_z, par_z, rh_can, gs_mol, &
atm2lnd_vars, photosyns_vars)

x0=x
exit
endif
if(iter>itmax)then
!in case of failing to converge within itmax iterations
!stop at the minimum function
!this happens because of some other issues besides the stomatal conductance calculation
!and it happens usually in very dry places and more likely with c4 plants.

call ci_func(minx,f1, p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z, rh_can, gs_mol, &
atm2lnd_vars, photosyns_vars)

exit
endif
enddo

x0=x
endif
end associate
end subroutine hybrid

!------------------------------------------------------------------------------
Expand Down Expand Up @@ -1137,7 +1129,7 @@ subroutine brent(x, x1,x2,f1, f2, tol, ip, iv, ic, gb_mol, je, cair, oair,&
type(photosyns_type), intent(inout) :: photosyns_vars
!
!!LOCAL VARIABLES:
integer, parameter :: ITMAX=20 !maximum number of iterations
integer, parameter :: ITMAX=30 !maximum number of iterations
real(r8), parameter :: EPS=1.e-2_r8 !relative error tolerance
integer :: iter
real(r8) :: a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm
Expand Down Expand Up @@ -1173,7 +1165,7 @@ subroutine brent(x, x1,x2,f1, f2, tol, ip, iv, ic, gb_mol, je, cair, oair,&
endif
tol1=2._r8*EPS*abs(b)+0.5_r8*tol !Convergence check.
xm=0.5_r8*(c-b)
if(abs(xm) <= tol1 .or. fb == 0.)then
if(abs(xm) <= tol1 .or. fb == 0._r8)then
x=b
return
endif
Expand Down Expand Up @@ -1212,13 +1204,13 @@ subroutine brent(x, x1,x2,f1, f2, tol, ip, iv, ic, gb_mol, je, cair, oair,&
call ci_func(b, fb, ip, iv, ic, gb_mol, je, cair, oair, lmr_z, par_z, rh_can, gs_mol, &
atm2lnd_vars, photosyns_vars)

if(fb==0._r8)exit
if(abs(fb)<1.e-5_r8)exit

enddo

if(iter==ITMAX)write(iulog,*) 'brent exceeding maximum iterations', b, fb
x=b

if(x<0._r8)write(iulog,*)'negative in brent',x,x1,x2
Copy link
Contributor

Choose a reason for hiding this comment

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

It seems that this write statement and the write statement at L1105 are duplicate.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It is now removed

return
end subroutine brent

Expand Down Expand Up @@ -1356,7 +1348,6 @@ subroutine ci_func(ci, fval, p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z,&
bbb => photosyns_vars%bbb_patch , & ! Output: [real(r8) (:) ] Ball-Berry minimum leaf conductance (umol H2O/m**2/s)
mbb => photosyns_vars%mbb_patch & ! Output: [real(r8) (:) ] Ball-Berry slope of conductance-photosynthesis relationship
)

! Miscellaneous parameters, from Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593
fnps = 0.15_r8
theta_psii = 0.7_r8
Expand Down Expand Up @@ -1402,25 +1393,26 @@ subroutine ci_func(ci, fval, p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z,&
! Net photosynthesis. Exit iteration if an < 0

an(p,iv) = ag(p,iv) - lmr_z
if (an(p,iv) < 0._r8) then
fval = 0._r8
return
endif
! Quadratic gs_mol calculation with an known. Valid for an >= 0.
! With an <= 0, then gs_mol = bbb

cs = cair - 1.4_r8/gb_mol * an(p,iv) * forc_pbot(c)
cs = max(cs,1.e-06_r8)
aquad = cs
bquad = cs*(gb_mol - bbb(p)) - mbb(p)*an(p,iv)*forc_pbot(c)
cquad = -gb_mol*(cs*bbb(p) + mbb(p)*an(p,iv)*forc_pbot(c)*rh_can)
call quadratic (aquad, bquad, cquad, r1, r2)
gs_mol = max(r1,r2)

if(an(p,iv)<=0.0)then
gs_mol=bbb(p)
if(aj(p,iv)<=1.e-20_r8)then
fval=0._r8
else
fval = ci-cair
endif
else
cs = cair - 1.4_r8/gb_mol * an(p,iv) * forc_pbot(c)
cs = max(cs,1.e-06_r8)
aquad = cs
bquad = cs*(gb_mol - bbb(p)) - mbb(p)*an(p,iv)*forc_pbot(c)
cquad = -gb_mol*(cs*bbb(p) + mbb(p)*an(p,iv)*forc_pbot(c)*rh_can)
call quadratic (aquad, bquad, cquad, r1, r2)
gs_mol = max(r1,r2)
! Derive new estimate for ci

fval =ci - cair + an(p,iv) * forc_pbot(c) * (1.4_r8*gs_mol+1.6_r8*gb_mol) / (gb_mol*gs_mol)

fval =ci - cair + an(p,iv) * forc_pbot(c) * (1.4_r8*gs_mol+1.6_r8*gb_mol) / (gb_mol*gs_mol)
endif
end associate

end subroutine ci_func
Expand Down