Hi Tristan:

It looks like your Python version of the code can handle `tau/mu = 500`

, for your stated values of `tau`

& `mu`

. That is a `tau/mu`

ratio that my calculator app can’t do `EXP(tau/mu)`

for. I can switch to DOUBLE precision, which can handle `EXP(tau/mu) = 10^308`

or `tau/mu = 709`

, if you’d like. But the Sellers code only uses SINGLE precision (“REAL”).

I do see two cos variables defined for the Sellers code in JULES. One is a local variable called `coszm`

, which is set to 1.0, and the comment says that it is the ‘mean value of COSZ’. This must be the angle of the normals of the leaves, and it must assume the leaves are all horizontal. The other is cosz_gb, which is the one I am using for `mu = cosz_gb`

in the `EXP(tau/mu)`

calculation for the L2SM code.

I don’t see anything like `EXP(tau/mu)`

with `mu = cosz_gb`

in the Sellers code. The closest thing I see for exponentials with positive arguments is `exp_hla = EXP(h * la)`

, where `la`

is the leaf area index, and where ` h = SQRT(b * b - c * c) / coszm`

, with `coszm = 1.0`

. This value for `h`

is not inversely proportial to cosz_gb.

The code in question is:

```
F = (1.0+k*mu)*(a1+k*g4)*EXP(k*tauD)
G = (1.0-k*mu)*(a1-k*g4)*EXP(-k*tauD)
H = 2.0*k*(g4+a1*mu)*EXP(tau/mu)
Tu = EXP(-tau/mu)*(1.0-w/D*(F-G-H))
```

If `EXP(tau/mu)`

is a huge number like 10^18 and if it can be represented in the precision of the compiler, then H is also likely going to be a huge number. Let’s call H’ = 2.0*k*(g4+a1*mu). If F & G are relatively small, then Tu is approximately w * H’ / D, since the exponentials cancel out.

It’s possible that I have things confused. But then again, maybe not.

One question I haven’t completely checked is if this `EXP(tau/mu)`

shouldn’t actually be `EXP(-tau/mu)`

. The `EXP(tau/mu)`

is what is in your original repository for L2SM. I can’t immediately/quickly figure out which of the equations in the Meador & Weaver 1980 paper correspond to these equations, though.

In other news, the suite `u-cx602`

(with `mu_min = 0.10`

) failed in the MAIN_RUN step, having a division by zero when it was trying to write out the NETCDF files for the ANNUAL data.

So, I decided that it might be crashing because I had `JULES_BUILD`

set to `'debug'`

(and hence properly crashing when there are divisions by zero, or the like). So I set `JULES_BUILD`

back to `'normal'`

in a new suite `u-cx563`

, and this suite runs properly. It even produces monthly GPP carbon fluxes that are >= 0.0, instead of often < 0.0 as the suite did before. So, a bit of progress has been made.

Patrick