From 20eb5992b1a639850857e1c7bb6cd3fd6605eb9a Mon Sep 17 00:00:00 2001 From: Jinyun Tang Date: Mon, 23 Jan 2017 13:23:34 -0800 Subject: [PATCH 1/3] Bug fix for the iterative solver of stomata conductance The original solver occasionaly leads to negative ci. This revision fixes this bug. Also, it was found that the old gs-An coupling did not guarantee a solution for all situations. For cases without solution, the revision forced the solution to converge with gs at this minimum value. A comparision between orginal and the fixed versions for 1850 SP simulations suggest that the difference is small, with a slight cooling in the tropics. --- .../clm/src/biogeophys/PhotosynthesisMod.F90 | 159 +++++++++--------- 1 file changed, 77 insertions(+), 82 deletions(-) diff --git a/components/clm/src/biogeophys/PhotosynthesisMod.F90 b/components/clm/src/biogeophys/PhotosynthesisMod.F90 index 1e161e6a51b1..f837f96632ac 100644 --- a/components/clm/src/biogeophys/PhotosynthesisMod.F90 +++ b/components/clm/src/biogeophys/PhotosynthesisMod.F90 @@ -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 ! @@ -733,7 +732,9 @@ 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))<1.e-14_r8)an(p,iv)=0._r8 ! Final estimates for cs and ci (needed for early exit of ci iteration when an < 0) @@ -1039,78 +1040,67 @@ subroutine hybrid(x0, p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z,& 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 + 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 + if(abs(f0) < 1.e-14_r8)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 - if(f1 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 + else + ! write(iulog,'(I8,13(X,E15.8),X,L8)')p,ci_val(1),fi_val(1),ci_val(2),fi_val(2),& + ! ci_val(3),fi_val(3),ci_val(4),fi_val(4),ci_val(5),fi_val(5),photosyns_vars%aj_patch(p,iv),& + ! photosyns_vars%ag_patch(p,iv),photosyns_vars%an_patch(p,iv), all(fi_val<0._r8) + ! write(iulog,*)'no solution for ci and gs_mol is forced to be the minimum for pft',p + ! write(iulog,*)'no solution found for ci' + ! call endrun(decomp_index=p, clmlevel=namep, msg=errmsg(__FILE__, __LINE__)) + + endif + end associate end subroutine hybrid !------------------------------------------------------------------------------ @@ -1142,7 +1132,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 @@ -1217,7 +1207,7 @@ 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 @@ -1407,25 +1397,30 @@ 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 +! 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 From 510644f51bd654b1b53993f2e140e1996579cb19 Mon Sep 17 00:00:00 2001 From: Jinyun Tang Date: Tue, 24 Jan 2017 09:42:33 -0800 Subject: [PATCH 2/3] Reset ag for negative an in photosynthesis calculation In cases when there is no solution for the gs-An copuled system, an is set to zero. This fix further sets ag to dark respiration to make it consistent with the use of zero an. --- components/clm/src/biogeophys/PhotosynthesisMod.F90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/components/clm/src/biogeophys/PhotosynthesisMod.F90 b/components/clm/src/biogeophys/PhotosynthesisMod.F90 index f837f96632ac..56eea7ec5898 100644 --- a/components/clm/src/biogeophys/PhotosynthesisMod.F90 +++ b/components/clm/src/biogeophys/PhotosynthesisMod.F90 @@ -734,7 +734,10 @@ subroutine Photosynthesis ( bounds, fn, filterp, & ! 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))<1.e-14_r8)an(p,iv)=0._r8 + if(abs(gs_mol(p,iv)-bbb(p))<1.e-14_r8)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) From 9f1764ba19604167bb339cf6a45e5ad14a79dd44 Mon Sep 17 00:00:00 2001 From: Jinyun Tang Date: Wed, 25 Jan 2017 13:33:59 -0800 Subject: [PATCH 3/3] limit ci re-calculation to minimum stomata conductance In the last fix, ci was still updated after obtaining stomata conductance from cally hybrid. This occasionally resulted in negative ci_z because of the limited precision that computer can realize. Now, ci_z is only calculated when stomata conductance is at its minimum. This ensures no negative ci_z in all cases. --- .../clm/src/biogeophys/PhotosynthesisMod.F90 | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/components/clm/src/biogeophys/PhotosynthesisMod.F90 b/components/clm/src/biogeophys/PhotosynthesisMod.F90 index 56eea7ec5898..83451601e03b 100644 --- a/components/clm/src/biogeophys/PhotosynthesisMod.F90 +++ b/components/clm/src/biogeophys/PhotosynthesisMod.F90 @@ -743,8 +743,15 @@ subroutine Photosynthesis ( bounds, fn, filterp, & 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))<1.e-14_r8)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 @@ -1094,6 +1101,9 @@ subroutine hybrid(x0, p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z,& lmr_z, par_z, rh_can, gs_mol, & atm2lnd_vars, photosyns_vars) x0=x + if(x0<0._r8)then + write(iulog,*)'negative after brent',x0 + endif else ! write(iulog,'(I8,13(X,E15.8),X,L8)')p,ci_val(1),fi_val(1),ci_val(2),fi_val(2),& ! ci_val(3),fi_val(3),ci_val(4),fi_val(4),ci_val(5),fi_val(5),photosyns_vars%aj_patch(p,iv),& @@ -1171,7 +1181,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 @@ -1216,7 +1226,7 @@ subroutine brent(x, x1,x2,f1, f2, tol, ip, iv, ic, gb_mol, je, cair, oair,& 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 return end subroutine brent @@ -1354,7 +1364,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