Open olyson opened 3 years ago
heroldn
commented
3 years ago On a related note, I was cross checking the code and CLM5 technote and saw that term b in equation 9.25 of the technote seems to be missing a minus sign out front, which exists in the code. Removing the minus sign gives negative conductance errors so I assume its inclusion is correct. But it might be worth checking the entire implementation of that quadratic with the source/author of the code.
djk2120
commented
3 years ago Would we expect the operational stomatal conductance to match the Medlyn equation in this way? I'll have to look through exactly how I coded this up, but I think of it more as Medlyn provides a maximal stomatal conductance, absent hydraulic stress, which is downregulated by PHS based on leaf water potential.
Whether this is the best approach is something I'm certainly interested in exploring, but I'm not sure this is a bug per se.
I can edit the the technical note, if this is not documented well.
EcoESM
commented
3 years ago Thanks for this note @djk2120. I had wondered how PHS was influencing this problem but hadn't gotten a chance to look into it yet. Given that the error check implemented by @heroldn is for maximum conductance and PHS downregulates this value, it is perhaps not surprising that the values don't match. Prior to PHS implementation, conductance was downregulated through the effect of btran on Vcmax (and also its effect on g0), so water stress was included more directly in the conductance equation.
I see that it is mentioned in the PHS chapter of the tech note, but I think it might be a good idea for us to update the stomatal conductance chapter to mention that the Medlyn conductance is downregulated by PHS.
danicalombardozzi
commented
3 years ago Thanks for this note @djk2120. I had wondered how PHS was influencing this problem but hadn't gotten a chance to look into it yet. Given that the error check implemented by @heroldn is for maximum conductance and PHS downregulates this value, it is perhaps not surprising that the values don't match. Prior to PHS implementation, conductance was downregulated through the effect of btran on Vcmax (and also its effect on g0), so water stress was included more directly in the conductance equation.
I see that it is mentioned in the PHS chapter of the tech note, but I think it might be a good idea for us to update the stomatal conductance chapter to mention that the Medlyn conductance is downregulated by PHS.
I was logged into a different GitHub account when I wrote this.
olyson
commented
3 years ago Thanks Daniel, that seems like a reasonable explanation. The stomatal conductance from the error check (the equation 2.9.1) always seems to be greater than the operational conductance in my single-point simulation so that seems to be consistent with your explanation. It might be good to confirm that they agree in a case with no water stress?
heroldn
commented
3 years ago Actually my values using eq. 9.1 seemed consistently lower than what CLM gave. And they still weren't matching even when PHS was off. But that was using the edits you made @djk2120 manually transferred into my own code, so it's possible I got something wrong when I combined. But still doesn't explain my lower values when using eq. 9.1.
heroldn
commented
3 years ago Digging a little further (with a healthy injection of print statements) I notice that the different gs values given by CLM (compared to eq. 9.1) are already present from the quadratic formula given by eq. 9.24 (and implemented in ci_func_PHS/ci_func). Should this make us suspect the quadratic formula? The technote says that the quadratic is arrived at by substituting eq. 9.23 for es into eq. 9.1, but I don't see where es factors into eq. 9.1, or how eq. 9.23 could be rearranged to substitute for something besides es. @djk2120, @danicalombardozzi , do you know who derived the numerical implementation?
A colleague here also noted that the technote doesn't explicitly define a for the quadratic (implicitly it's "1"), and I note a typo or two in this section as well, so I wonder if something got missed and I'd be keen to confirm these equations with the author. Maybe it's nothing but still.
olyson
commented
3 years ago I'm just copying in an offline conversation that happened over the last week or so, so that there is a record of it here.
Hi Rosie and Gordon
Thanks for this. And thanks Gordon your book looks like a great reference and it looks like it's resolved the issue too, but I'd appreciate all your eyes on this to confirm, and there are definitely some inconsistencies remaining (see below). Apologies for the long email.
Basically, the analytical version of Medlyn Gs, according to Gordon's book eq. 12.21, should use the leaf surface VPD. The VPD in the source code and that I've used in my check is the leaf-to-air VPD. The technote could use revising in this regard because it uses the same symbol D to refer to both types of VPD (but I see it does say in paragraph form just under eq. 9.1 that that equation uses leaf surface VPD, but people like me just look at symbols!). Once I implemented leaf surface VPD in my check the values from the analytical equation matched much more closely with what comes from the quadratic form, often to within a few decimal places, but occasionally it's still extremely different too (attached is grepped output from a 5 day F case that I've been eye-balling). Most errors seem < 50 and thus represent a small % error, since Gs is usually on the order of tens of thousands. So, IF there was going to be an internal check in the code that compares the analytical and numerical solutions (as there is for BB), the error tolerance would have to be much higher than that used for BB (which is 0.1) and maybe % based. It would still fail occasionally as I see extremely large differences once in a while (factor of 2 or 3 different and Gs values can go into the millions, but these look rare). Not sure what the tolerance in any check should be but of course your intuitions are better than mine.
Now, some small inconsistencies to note/clarify/revise in future tech notes: The technote has a slight mathematical error, in eq. 9.25 term c of the quadratic has (1-g^2)/D where it should be 1 - g^2/D as per the code and Gordon's book (eq. 12.23). Gordon, while your textbook seems to have resolved this issue, I note that your equation 12.21 says leaf surface VPD should be used in Medlyn but that the Medlyn 2011 paper itself seems to indicate that the leaf-to-air VPD should be used. Have I missed something or have they mistyped (I note there's a correction published for that paper)? As stated above, the technote should have the symbols for VPD updated so that the D in equation 9.1 is distinguished from the D used in the numerical implementation. This is the source of the inconsistency. The text probably needs checking in chapter 9.6. One mistake is the single sentence between eq. 9.2 and 9.22 which was complete in previous versions of the technote but chopped up a bit here and might be irrelevant. And eq. 9.23 of the technote probably needs replacing with eq. 12.16 in Gordon's book. Rosie, hmm ok I'll have to look into putting a couple of FATES runs on the standby queue! What is the performance hit (roughly) of activating FATES?
Thanks for everyone's input. If we can confirm this issue is resolved (am I correct that the occasional large discrepancy is OK?), then I can add the above details to the github report.
Cheers
Nick. On 15/07/2021 05:38 rosie fisher rosieafisher@googlemail.com wrote:
Indeed, the quadratic solver predates the Medlyn implementation, and Gordon's book is probably your best resource on that...
FATES isn't just for Paleo! We now have a 'satellite phenology' version just for looking into canopy physiology issues like this too ;)
Keep us in the loop with what you find.
-Rosie
Le mer. 14 juil. 2021 à 18:24, Gordon Bonan bonan@ucar.edu a écrit : Nick,
Look at my book Climate Change and Terrestrial Ecosystem Modeling (Cambridge University Press, 2019). Chapter 12 on stomatal conductance gives the derivation and numerical solution needed for the Ball-Berry and Medlyn models (see specifically eqs. 12.14-12.23). There is accompanying matlab code (explained in the book). The code can be found at https://github.com/gbonan/bonanmodeling
Gordon
On Wed, Jul 14, 2021 at 6:00 AM Nicholas Herold heroldn@mailbox.org wrote: Hi Rosie
Thanks very much for your reply. I can understand if this is stretching your memory so really appreciate your thoughts on it.
What you say about the technote consisting of old text from the BB method makes sense from my comparison of the different versions of that document. I actually have been looking at the code and there is actually a github issue logged for this here. I didn't include it initially to try and keep my email from being an information dump! If the es equation in the technote is redundant then I guess my question would be how does one arrive at the quadratic equation used in the code to get Gs (line 4159 of your link and described in eq. 9.24-9.26)? It's this quadratic that's giving different values compared to the Medlyn equation (eq. 9.1). You can see from the github comments Daniel suggested it could be because PHS downregulates Gs, but if I understand the code and methods correctly the Gs coming out of the quadratic should be prior to any downregulating and thus should match eq. 9.1 (I'm not a land modeller so have been learning many a new concept during all this!). I also tested this without PHS and still get inconsistent values, so I'm ruling out PHS at the moment, and that's led me to the quadratic getting used. Do you have any recollection of how this is derived by any chance?
Thanks for the heads up about FATES too, I wasn't aware. In my paleo CESM days this would have been right up my alley (I was using the immensely simple BIOME4!). I'll definitely keep this in mind.
Cheers
Nick. On 14/07/2021 16:39 rosie fisher rosieafisher@googlemail.com wrote:
-actually- cc-ing you guys!
---------- Forwarded message --------- From: rosie fisher rosieafisher@googlemail.com Date: Wed, 14 Jul 2021, 08:37 Subject: Re: CLM5 Medlyn Gs numerical implementation To: Nicholas Herold heroldn@mailbox.org
Dear Nick, (cc Gordon, Danica, Keith)
Thanks for this. On your first question (did I implement this?) I had to sit down for a long while to try and remember! I think I might have done, but my recollection appears to be a little foggy. A lot has happened since then, to put it mildly.
Either way, looking at the tech note, what I think might have happened somewhat prosaically, is that we updated the stomatal conducts ce équation to Medlyn, but then didn't realize that the later text on es refera to it. My suspicion is that this (ew 23) was relevant to the previous explanation of BB. Of course, Medlyn is using VPD and so that's not how it works any more.
My suggestion is thus probably to look directly at the code (e.g. working from line 3311 here) https://github.com/ESCOMP/CTSM/blob/master/src/biogeophys/PhotosynthesisMod.F90
As you can see there, the stomatal model is really a one-time switch, (with a couple of upstream calculations of vod, etc) so if we did something wrong it shouldn't be too hard to identify, hopefully.
Let me/us know if that throws up any questions. Also you might be interested to know (if you don't already) that we are also testing the implementation of Medlyn vs BB in the FATES model extension of CLM with Alister Rogers at BNL and his post doc Cherry. They are at the stage of writing up their work but I'm quite sure they have no plans (/funding) to look at the climate impacts. You might be interested in exchanging note with them at some point.
Best wishes, -Rosie
Dr Rosie A. Fisher. Scientist, Évolution & Diversité Biologique, University of Toulouse Paul Sabatier III, France. & Affiliate Scientist National Center for Atmospheric Research, Boulder, Colorado, USA
On Wed, 14 Jul 2021, 01:24 Nicholas Herold, heroldn@mailbox.org wrote: Hi Rosie
My name's Nicholas Herold and I'm using CLM5 to study the climatic impact of different stomatal conductance schemes, namely Ball-Berry and Medlyn. I've been conversing with Keith, Danica and Daniel from NCAR on an issue I'm having reconciling Gs values I get from CLM with ones I calculate with the Medlyn equation. And I'm trying to trace where the source of this difference comes from. I was told you may have implemented the numerical solution for Medlyn in CLM5. Is that correct?
If so I was wondering if you could help join some dots for me. The CLM5 technote says that eq. 9.23 is substituted into eq. 9.1 to get the quadratic equation described in equations 9.24 - 9.26. But I don't see how this substitution works (i.e. how es can substitute into eq. 9.1). I'm interested in this because the solution to the quadratic that CLM is coming up with is different to what I get from eq. 9.1, which, if I understand correctly, shouldn't be happening (but I'd be happy to be told otherwise). This happens whether the plant hydraulic stress module is active or not.
Anyway, any insight would be appreciated. Thanks for your time.
Cheers
Nick.
heroldn
commented
3 years ago Ok, a bit more progress. Below I'm implementing the analytical equation after calculating es and vpd. There are a couple of modifying factors applied to VPD and g0 in the CLM source which I've included below (a minimum constraint for VPD, and g0 being multiplied by bsun). Removing the latter brings gs values into even closer agreement. But still some large exceptions occur. Removing the minimum constraint for VPD seems to improve things a little further, i.e. agreement is to the point that the 99.4th percentile of errors is 8.38190317153931E-09. But with the VPD constraint removed I see the maximum discrepancy goes up to 16,797,683,955 and is due to the analytical solution going bonkers when the numerical solution is close to zero - so I imagine that the VPD constraint is there for good reason and should be kept (comment in the code implies Jinyun Tang might know something about this constraint)!
So, it seems removing the bsun multiplication is a genuine improvement (and I see this modification isn't applied in the Medlyn numerical solution so that makes sense) but maybe the VPD constraint should be kept. Can someone comment on this? Does the VPD minimum constraint need tweaking to improve treatment of edge cases? I can't see any other reason for the remaining large discrepancies. When only removing the bsun multiplication, the max discrepancy is 5,039,221 and this seems to be when the numerical solution gets big (6,841,128) compared to the analytical solution (1,801,906). And the 95th percentile of the errors is 91.35. So it seems at first glance that the VPD constraint improves the most extreme discrepancies (~16.7 billion down to ~5 million) but increases discrepancies elsewhere.
es = (ceair*gb_mol(p) + esat_tv(p)*gs_mol_sun(p,iv)) / (gs_mol_sun(p,iv)+gb_mol(p)) !max((esat_tv(p) - ceair), 50._r8) * 0.001_r8
vpd = max((esat_tv(p) - es), 50._r8) * 0.001_r8
gs_mol_err1 = max(bsun(p)*gsminsun, 1._r8) + 1.6_r8*(1._r8+(gs_slope_sun/sqrt(vpd))) * (an_sun(p,iv)/(cs_sun/forc_pbot(c)))I've tried a slightly different approach following some of Gordon's Matlab code. The following code for sunlit doesn't throw any errors in a 2 year run at the US-UMB tower site. No errors on the shaded side either. According to Gordon's code, the error check is only done (valid) when (esat_tv(p) - ceair) .gt. 50. I've removed the multiplication of g0 by bsun except when an < 0 (where gsminsun is multiplied by bsun in the photosynthesis code). hs and rh_leaf_sun already exist in the code because of the Ball-Berry error check, but this could be simplified using es = hs * esat_tv as you've done.
hs = (gb_mol(p)*ceair + gs_mol_sun(p,iv)*esat_tv(p)) / ((gb_mol(p)+gs_mol_sun(p,iv))*esat_tv(p))
rh_leaf_sun(p) = hs
if ((esat_tv(p) - ceair) .gt. 50) then
vpd = max(esat_tv(p) - rh_leaf_sun(p) * esat_tv(p), 0.1_r8) * 0.001_r8
if (an_sun(p,iv) .ge. 0._r8) then
gs_mol_err = gsminsun + 1.6_r8*(1._r8+(gs_slope_sun/sqrt(vpd))) * (max(an_sun(p,iv),0._r8)/(cs_sun/forc_pbot(c)))
else
gs_mol_err = max(bsun(p)*gsminsun, 1._r8) + 1.6_r8*(1._r8+(gs_slope_sun/sqrt(vpd))) * (max(an_sun(p,iv),0._r8)/(cs_sun/forc_pbot(c)))
end if
if (abs(gs_mol_sun(p,iv)-gs_mol_err) > 1.e-01_r8) then
write (iulog,*) 'Medlyn error check - SUNLIT stomatal conductance error:'
write (iulog,*)'gs_sun: ',gs_mol_sun(p,iv),'gs_err: ',gs_mol_err
end if
end ifThat makes sense to me, Keith.
Also, sorry if I missed something, but will this PR add a bsun/bsha factor for g0 with Medlyn?
That seems more like a science change than a numerical bug fix.
My impression was that many supported the change with the Medlyn implementation to not attentuate g0.
On Tue, Jul 20, 2021 at 9:55 AM Keith Oleson @.***> wrote:
I've tried a slightly different approach following some of Gordon's Matlab code. The following code for sunlit doesn't throw any errors in a 2 year run at the US-UMB tower site. No errors on the shaded side either. According to Gordon's code, the error check is only done (valid) when (esat_tv(p) - ceair) .gt. 50. I've removed the multiplication of g0 by bsun except when an < 0 (where gsminsun is multiplied by bsun in the photosynthesis code). hs and rh_leaf_sun already exist in the code because of the Ball-Berry error check, but this could be simplified using es = hs * esat_tv as you've done.
hs = (gb_mol(p)*ceair + gs_mol_sun(p,iv)*esat_tv(p)) / ((gb_mol(p)+gs_mol_sun(p,iv))*esat_tv(p)) rh_leaf_sun(p) = hs if ((esat_tv(p) - ceair) .gt. 50) then vpd = max(esat_tv(p) - rh_leaf_sun(p) * esat_tv(p), 0.1_r8) * 0.001_r8 if (an_sun(p,iv) .ge. 0._r8) then gs_mol_err = gsminsun + 1.6_r8*(1._r8+(gs_slope_sun/sqrt(vpd))) * (max(an_sun(p,iv),0._r8)/(cs_sun/forc_pbot(c))) else gs_mol_err = max(bsun(p)*gsminsun, 1._r8) + 1.6_r8*(1._r8+(gs_slope_sun/sqrt(vpd))) * (max(an_sun(p,iv),0._r8)/(cs_sun/forc_pbot(c))) end if if (abs(gs_mol_sun(p,iv)-gs_mol_err) > 1.e-01_r8) then write (iulog,*) 'Medlyn error check - SUNLIT stomatal conductance error:' write (iulog,*)'gs_sun: ',gs_mol_sun(p,iv),'gs_err: ',gs_mol_err end if end if— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/ESCOMP/CTSM/issues/1398#issuecomment-883546516, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADNHR5WY2FUMHV4UVN7P62LTYWTA5ANCNFSM46RKMPGQ .
olyson
commented
3 years ago The existing CLM5 release code, that has Medlyn with PHS, has attenuated g0, but only for an<0:
if (an_sun(p,iv) < 0._r8) gs_mol_sun(p,iv) = max( bsun(p)gsminsun, 1._r8 ) if (an_sha(p,iv) < 0._r8) gs_mol_sha(p,iv) = max( bsha(p)gsminsha, 1._r8 )
and so the error check has to accommodate that.
But my impression for the Medlyn-SMS PR was that we don't want to attenuate g0.
djk2120
commented
3 years ago That seems like a bug to me. I think the intent is to have the Medlyn model g0 unattenuated.
On Tue, Jul 20, 2021 at 10:35 AM Keith Oleson @.***> wrote:
The existing CLM5 release code, that has Medlyn with PHS, has attenuated g0, but only for an<0:
if (an_sun(p,iv) < 0._r8) gs_mol_sun(p,iv) = max( bsun(p)gsminsun, 1._r8 ) if (an_sha(p,iv) < 0._r8) gs_mol_sha(p,iv) = max( bsha(p)gsminsha, 1._r8 )
and so the error check has to accommodate that.
But my impression for the Medlyn-SMS PR was that we don't want to attenuate g0.
— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/ESCOMP/CTSM/issues/1398#issuecomment-883570882, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADNHR5XVQKIY4G65T4S7WVDTYWXVTANCNFSM46RKMPGQ .
olyson
commented
3 years ago I'm fine with changing that, it will make the error check simpler also.
djk2120
commented
3 years ago We'll just have to wrap a medlyn/BB if around those statements, because I think we do still want the bsun/bsha attenuation for ball-berry
On Tue, Jul 20, 2021 at 10:46 AM Keith Oleson @.***> wrote:
I'm fine with changing that, it will make the error check simpler also.
— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/ESCOMP/CTSM/issues/1398#issuecomment-883578100, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADNHR5XS3IZFLNJS2TJKCFDTYWY7DANCNFSM46RKMPGQ .
olyson
commented
3 years ago Good point. I'm also going to do a short global simulation to see if the error check is robust. Maybe Nick can review my changes to the error check and verify that it seems correct and works in his own simulation.
wwieder
commented
3 years ago I have to confess, I'm only following this discussion superficially, but I wondered if it would be helpful to discuss at an upcoming CLM meeting?
heroldn
commented
3 years ago Great. Ok a 10 day global case doesn't show any errors on my end using @olyson's formulation. And it's fine with PHS on or off. For the record I'm testing in my own modified CLM code but don't see why this wouldn't be robust. I realise in my implementation above I hadn't applied that minimum constraint of zero for an either.
To implement this check in the Photosynthesis routine (i.e. when PHS is off) I removed the if condition for an >= 0 which is what I understood from the above exchange. Combined BB and MED error check for PHS = off:
hs = (gb_mol(p)*ceair + gs_mol(p,iv)*esat_tv(p)) / ((gb_mol(p)+gs_mol(p,iv))*esat_tv(p))
rh_leaf(p) = hs
! Compare with Ball-Berry or Medlyn model.
! BB: gs_mol = m * an * hs/cs p + b
! MED: gs_mol = g0 + 1.6*(1+(g1/sqrt(D))) * (an/(cs/patm))
if ( stomatalcond_mtd == stomatalcond_mtd_bb1987 ) then
gs_mol_err = mbb(p)*max(an(p,iv), 0._r8)*hs/cs*forc_pbot(c) + bbb(p)
if (abs(gs_mol(p,iv)-gs_mol_err) > 1.e-01_r8) then
write (iulog,*) 'Ball-Berry error check - stomatal conductance error (CLM, mycalc, lat, lon, phase):'
write (iulog,*) gs_mol(p,iv), gs_mol_err, grc%latdeg(g), grc%londeg(g), phase
end if
else if (stomatalcond_mtd == stomatalcond_mtd_medlyn2011) then
if ((esat_tv(p) - ceair) .gt. 50) then
vpd = max(esat_tv(p) - rh_leaf(p) * esat_tv(p), 0.1_r8) * 0.001_r8
! if (an(p,iv) .ge. 0._r8) then
gs_mol_err = medlynintercept(patch%itype(p)) + 1.6_r8*(1._r8+(medlynslope(patch%itype(p))/sqrt(vpd))) * (max(an(p,iv),0._r8)/(cs/forc_pbot(c)))
! else
! gs_mol_err = max(bsun(p)*medlynintercept(patch%itype(p)), 1._r8) + 1.6_r8*(1._r8+(medlynslope(patch%itype(p))/sqrt(vpd))) * (max(an(p,iv),0._r8)/(cs/forc_pbot(c)))
! end if
write (iulog,*) '*************** MEDLYN CHECK (NUMERICAL, ANALYTICAL, DIFFERENCE): ',gs_mol(p,iv), gs_mol_err, gs_mol(p,iv)-gs_mol_err
if (abs(gs_mol(p,iv)-gs_mol_err) > 1.e-01_r8) then
write (iulog,*) 'Medlyn error check - stomatal conductance error:'
write (iulog,*)'gs: ',gs_mol(p,iv),'gs_err: ',gs_mol_err
end if
end if
end ifAnd to confirm what you were suggesting @djk2120. Should the following two lines inside PhotosynthesisHydraulicStress:
if (an_sun(p,iv) < 0._r8) gs_mol_sun(p,iv) = max( bsun(p)*gsminsun, 1._r8 )
if (an_sha(p,iv) < 0._r8) gs_mol_sha(p,iv) = max( bsha(p)*gsminsha, 1._r8 )be replaced with something like:
if (an_sun(p,iv) < 0._r8) then
if (stomatalcond_mtd == stomatalcond_mtd_bb1987) then
gs_mol_sun(p,iv) = max( bsun(p)*gsminsun, 1._r8 )
else if ( stomatalcond_mtd == stomatalcond_mtd_medlyn2011 )then
gs_mol_sun(p,iv) = max( gsminsun, 1._r8 )
end if
end if
if (an_sha(p,iv) < 0._r8) then
if (stomatalcond_mtd == stomatalcond_mtd_bb1987) then
gs_mol_sha(p,iv) = max( bsha(p)*gsminsha, 1._r8 )
else if ( stomatalcond_mtd == stomatalcond_mtd_medlyn2011 )then
gs_mol_sha(p,iv) = max( gsminsha, 1._r8 )
end if
end ifIs there any equivalent change required when PHS = off?
olyson
commented
3 years ago @heroldn would you mind reviewing the photosynthesis chapter in the technical note? I've updated the document per your suggestions:
Let me know if I missed anything.
heroldn
commented
3 years ago I've read the parts that I made suggestions on (certainly not qualified to review 2.9 in general).
heroldn
commented
3 years ago @djk2120 @olyson is my code sample above what you meant RE attenuating g0 under BB but not MED. And does this need to be reproduced when PHS is off?
olyson
commented
3 years ago Thanks for the review. I'm not quite sure what you mean by a formatting error with "a=1", unless you mean how these equations are aligned.
Yes, I think your code sample is right and this needs to be reproduced when PHS is off.
heroldn
commented
3 years ago Ah yes, that equation is just right-aligned, it looks inconsistent with other equations.
Ok thanks re g0 attenuation. And when PHS=off, I now see that the way variables are named ensures g0 attenuation is already only applied to BB (in my code anyway - I was taking edits to allow Medlyn from Daniel's files). Thanks for your time on all this.
Brief summary of bug
Implementing an error check on medlyn stomatal conductance, similar to the Ball-Berry error check, indicates an error
General bug information
CTSM version you are using: release-clm5.0.34
Does this bug cause significantly incorrect results in the model's science? Maybe
Configurations affected: All using PhotosynthesisHydraulicStress (use_hydrstress=.true.)
Details of bug
A user tried to implement an error check for Medlyn stomatal conductance (sunlit and shaded) in PhotosynthesisHydraulicStress, similar to the error check that is done for Ball-Berry in Photosynthesis. Specifically, the error check calculates the difference between stomatal conductance produced by the solution in PhotosynthesisHydraulicStress and the stomatal conductance produced by equation 2.9.1 in the technical note, and triggers an error statement if the difference is greater than 0.1. Here is the added code for sunlit stomatal conductance:
Here, gs_mol_err represents the stomatal conductance as calculated by equation 2.9.1 in the technical note. I've replicated this behavior in a single point simulation at US-UMB. Here is some sample output:
Medlyn error check - SUNLIT stomatal conductance error: 30503.1558514552 30136.0322013949
It looks to me as if the code for gs_mol_err is consistent with the tech note equation. rh_can is not actually relative humidity, it is the vapor pressure deficit at the leaf surface which is calculated correctly in the code. I don't know if the equation is actually wrong or if there is something wrong in how gs_mol_sun(sha) is calculated in the PhotosynthesisHydraulicStress. Also not sure how to figure this out.
Important details of your setup / configuration so we can reproduce the bug
This case demonstrates the problem:
/glade/work/oleson/release-clm5.0.34/cime/scripts/spinclm50conr34sp_US-UMB_chkmedlyn