# Add L2SM model into JULES **URL:** \ **Category:** JULES\ **Created:** [1 March 2023 13:20 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928 "2023-03-01T13:20:51Z")\ **Posts on this page:** 20\ **Page:** 2
**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\ **Post date:** [14 June 2023 15:18 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/21 "2023-06-14T15:18:26Z")
Hi Tristan: Well, for `tau/mu = 38` (at least in one case, with `0.01 <= mu <= 0.05`), `EXP(tau/mu)` is a huge number (`3.2 * 10^(18)`), so I am not surprised that there are issues with overflow. There are probably other cases with larger values of tau/mu. The natural log of `10^99` is `ln(10^99) = 228`, which might be the upper limit on what values of `tau/mu` that can get through without overflow. I am not sure if such high values of `EXP(tau/mu)` (`=10^18`) are meaningful. Here’s a table of `mu=cos(solar zenith angle)` versus `solar zenith angle`. mu` `SZA 0.01 89.43deg 0.05 87.13deg 0.10 84.26deg I can try to see if there are overflows of `EXP(tau/mu)` when `mu = 0.07` or something like that, but the angle for `mu = 0.10` is 84.26deg, which is very low on the horizon already. So, maybe the limit of `mu=0.10` is a sufficiently-low threshold. Patrick ---
**Author:** ![tquaife](https://avatars.discourse-cdn.com/v4/letter/t/22d042/32.png) [@tquaife](https://cms-helpdesk.ncas.ac.uk/u/tquaife)\ **Post date:** [14 June 2023 16:01 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/22 "2023-06-14T16:01:33Z")
Hi Patrick, My version of the code (in Python) can handle mu=0.01 and tau=5.0. Is it really the case the Fortran is _less_ capable than Python numerically? Also worth noting that value of tau in a single layer is unlikely to ever exceed 1.0 using the default ten layer canopy, so that allows us to drop to even lower values of mu. I appreciate what the consequences are for sun angle, but JULES still needs a way of returning values of absorbed radiation and the canopy albedo at low sun angles. If it worked previously with the existing Sellers canopy RT in JULES it should also work for L2SM. They are, ultimately, both derived from the same 2-stream solution. Tristan. ---
**Author:** ![tquaife](https://avatars.discourse-cdn.com/v4/letter/t/22d042/32.png) [@tquaife](https://cms-helpdesk.ncas.ac.uk/u/tquaife)\ **Post date:** [14 June 2023 17:09 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/23 "2023-06-14T17:09:28Z")
Hi Patrick, I have a solution to propose. Assuming the JULES Fortran code really can’t handle these values, then we still need to be able to provide values for albedo and absorption/fAPAR. To do this, I have derived expressions for the reflectance (R) when mu=0. We also know that transmission (T) must be 0, and therefore fAPAR=1-R (in the visible/PAR part of the spectrum) for the top layer, and 0 for all subsequent layers (because T=0). How does that sound? ---
**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\ **Post date:** [14 June 2023 17:20 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/24 "2023-06-14T17:20:15Z")
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: ```auto 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 ---
**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\ **Post date:** [14 June 2023 17:23 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/25 "2023-06-14T17:23:15Z")
Hi Tristan: I was writing the last post to respond to your 2nd to last post when you sent your last post. Maybe you can read my last post first before we decide what to do. Your solution for mu=0 sounds very useful, and it might be needed! But maybe there are other solutions. Patrick ---
**Author:** ![tquaife](https://avatars.discourse-cdn.com/v4/letter/t/22d042/32.png) [@tquaife](https://cms-helpdesk.ncas.ac.uk/u/tquaife)\ **Post date:** [14 June 2023 17:45 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/26 "2023-06-14T17:45:28Z")
Hi Patrick, coszm in JULES is the average inverse diffuse optical depth per unit leaf area. The fact it is always equal to one is indeed because the the limited selection of leaf geometry archetypes, but it does not represent the leaf normals (it just happens the coszm is unity for both uniform and horizontal leaf distributions). However - you are correct that it’s not mu as in the cosine of the sun angle, and so that is why the Seller’s implementation dodges the same overflow problem (at least in the JULES implementation; coszm can vary in other implementation, e.g. in CLM). The exp(tau/mu) is correct. The minus sign really should be absent. It appears in the direct transmission term, which is equation 15 in M&W. I don’t totally follow the need for your approximation to T with mu-\>0. It is just T=0, because exp(-tau/mu)-\>0 as mu-\>0. Like I mentioned, I also have an exact solution for reflection and absorption when mu=0, so we can use those for low sun angles. Thanks for all the effort you are putting into this, it is much appreciated. Tristan. ---
**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\ **Post date:** [14 June 2023 21:46 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/27 "2023-06-14T21:46:41Z")
Thanks, Tristan, for helping to sort this out. I see equation 15 now. Yes, the approximations for when mu-\>0 would be useful to have. As far as my reasoning goes, here is the spelt-out maths: ```auto H = 2.0*k*(g4+a1*mu)*EXP(tau/mu) \def H' * EXP(tau/mu) Tu = EXP(-tau/mu)*(1.0-w/D*(F-G-H)) For abs(H) >> abs(G) and abs(H) >> abs(F) and H >> 1.0: Tu \approx EXP(-tau/mu)*(w/D * H) = EXP(-tau/mu)*(w/D * H' * EXP(tau/mu)) = w/D * H' So Tu is not tending to 0 in these cases (unless w or H' tend to zero). ``` Or maybe I made a maths mistake or improper approximation somewhere? Patrick ---
**Author:** ![tquaife](https://avatars.discourse-cdn.com/v4/letter/t/22d042/32.png) [@tquaife](https://cms-helpdesk.ncas.ac.uk/u/tquaife)\ **Post date:** [15 June 2023 07:04 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/28 "2023-06-15T07:04:03Z")
But also exp(-tau/mu) will go to zero for small mu, so hence T will also be zero. I am in an ESA workshop for the next two days. I’ll check my working on the reflectance term over the weekend (I think it’s fairly simple) and up date my Python version with it. ---
**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\ **Post date:** [15 June 2023 07:48 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/29 "2023-06-15T07:48:32Z")
Hi Tristan: Yes. EXP(-tau/mu) goes to zero for small mu. But it is multiplied by the H term that includes an EXP(+tau/mu) term. The product of these two terms is a constant that doesn’t go to zero with small mu. Or maybe I am misunderstanding something? Patrick ---
**Author:** ![tquaife](https://avatars.discourse-cdn.com/v4/letter/t/22d042/32.png) [@tquaife](https://cms-helpdesk.ncas.ac.uk/u/tquaife)\ **Post date:** [15 June 2023 11:26 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/30 "2023-06-15T11:26:35Z")
Hi Patrick, My apologies - you are correct. I have made some changes in my Python version for T, based on your suggestion, and corresponding ones to R, and run some tests. Works to a _very_ high order of approximation. Haven’t pushed to github yet, but will let you know when I do. Tristan. ---
**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\ **Post date:** [21 June 2023 19:36 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/31 "2023-06-21T19:36:17Z")
Hi Tristan I’m changing the code to the following. This will avoid the `EXP(tau/mu)` for small `mu`, so `mu` as low as 0.001 should work. We could make this part of the code work even for `mu = 0`, since (`EXP(-1/x)->0` as `x->0+ `, so we could put the limiting case in. But there are other divisions by `mu` elsewhere in the code that would need handling somehow. Patrick ```auto !IF ( mu .GE. 0.1 ) THEN !setting this limit to 0.01 causes overflows with EXP(tau/mu) IF ( mu .GE. 0.001 ) THEN DO i = ilayers,1,-1 . . . 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) H2 = 2.0*k*(g4+a1*mu) ! Tu = EXP(-tau/mu)*(1.0-w/D*(F-G-H)) Tu = EXP(-tau/mu)*(1.0-w/D*(F-G))+w*H2/D . . . END DO END IF ``` ---
**Author:** ![tquaife](https://avatars.discourse-cdn.com/v4/letter/t/22d042/32.png) [@tquaife](https://cms-helpdesk.ncas.ac.uk/u/tquaife)\ **Post date:** [22 June 2023 08:34 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/32 "2023-06-22T08:34:03Z")
Perfect. Thanks Patrick. I forgot to commit my changes to github, but it is essentially the same thing I think (I haven’t double checked). ---
**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\ **Post date:** [7 September 2023 21:13 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/33 "2023-09-07T21:13:19Z")
Hi Tristan: I did JULES/L2SM simulations of the region around the western UK, covering 5x16 grid cells. The grid was an N96 global grid. I produced layered-canopy output, and I also added JULES switches to produce monthly averages of the layered-canopy output. The simulation had the JULES GL7 settings, and went from 1860-1947. The simulation was run as Rose/Cylc suite `u-cx688e_UK`, which I checked into MOSRS as `u-cz631`. I plot below the monthly averages of the GPP and the APAR in the different canopy layers of the Deciduous Broadleaf PFT for the year 1947. This is for a single N96 grid cell, which is centered at Lat = 51.875 & Lon = -2.8125 degrees, near the border of England and southern Wales. The plots make sense to me, with decreasing values with canopy depth. And there are some shadowing-like effects in the summer time, where possibly the increased leaf-areas cause reduced values of GPP and APAR for the deeper layers. I am also running this globally, on an N96 grid using 16 processors, as a Rose/Cylc suite `u-cx688d`, checked into MOSRS as `u-cx688`. The data for the UK suite has been copied from `scratch-pw2` more permanently to the `landsurf_rdg` group workspace at the path `/gws/nopw/j04/landsurf_rdg/pmcguire/TristanQ/u-cx688e_UK`. Patrick ![cfplot_gpppftlyr_v7i1e](https://europe1.discourse-cdn.com/flex013/uploads/cms_support/original/1X/dbcf08d008d6ee17b7b9a606af32537dc766a86b.png) ![cfplot_aparpftlyr_v7i1f](https://europe1.discourse-cdn.com/flex013/uploads/cms_support/original/1X/9d16f142f6d98890ebfc0107253bbf6851294b6c.png) ---
**Author:** ![tquaife](https://avatars.discourse-cdn.com/v4/letter/t/22d042/32.png) [@tquaife](https://cms-helpdesk.ncas.ac.uk/u/tquaife)\ **Post date:** [8 September 2023 13:15 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/34 "2023-09-08T13:15:43Z")
Amazing stuff. Thanks Patrick. I am away for the next two weeks at meetings, but we should have a chat when I get back. How is your availability for the 27th and 29th? ---
**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\ **Post date:** [8 September 2023 14:05 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/35 "2023-09-08T14:05:37Z")
Thanks, Tristan! I will answer you privately about meeting up. Patrick ---
**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\ **Post date:** [27 September 2023 12:43 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/36 "2023-09-27T12:43:11Z")
Hi Tristan: It was nice talking with you today. The global JULES runs with can\_rad\_mod=6 and canopy-layered output that I started yesterday have made it from the year 1860 to about 1890 (and are continuing on to 1947). This would be with the Sellers radiative transfer for JULES. The previous global runs with the new can\_rad\_mod=7 L2SM option for JULES also finished a while back. This L2SM option is with my JULES FORTRAN code from your L2SM Python code for the Meador and Weaver equations. I compare the annual layered curves below for my JULES with either can\_rad\_mod=7 L2SM or can\_rad\_mod=6 Sellers below (suites u-cx688 and u-da046, respectively). This is for one grid cell in western England and eastern Wales, for the year 1869. I had plotted the year 1947 in a previous post here, for the can\_rad\_mod=7 L2SM only run for a subgrid around the western UK. I plot the GPP and the APAR in layers. There is indeed more work to be done to get the model output for the can\_rad\_mod=7 L2SM/Meador&Weaver code and the can\_rad\_mod=6 Sellers code to match up. The APAR is higher in the top layers for the new code, which results in the GPP being lower for all but the top layer with the new code. Patrick **GPP can\_rad\_mod=6 Sellers** ![cfplot_gpppftlyr_v7i1g](https://europe1.discourse-cdn.com/flex013/uploads/cms_support/original/1X/ef418c7f07ffc48d9fcbfc44604149ebb4ed461e.png) **GPP can\_rad\_mod=7 L2SM/Meador&Weaver** ![cfplot_gpppftlyr_v7i1h](https://europe1.discourse-cdn.com/flex013/uploads/cms_support/original/1X/48e6a7a8113327638144b6cb4904545c60883a17.png) **APAR can\_rad\_mod=6 Sellers** ![cfplot_aparpftlyr_v7i1g](https://europe1.discourse-cdn.com/flex013/uploads/cms_support/original/1X/a5725a1e8b6beb3b852f288f946ebfe3c1abea09.png) **APAR can\_rad\_mod=7 L2SM/Meador&Weaver** ![cfplot_aparpftlyr_v7i1h](https://europe1.discourse-cdn.com/flex013/uploads/cms_support/original/1X/96f71f1ca4cd65555de3e36847cff0ea2763ed86.png) ---
**Author:** ![tquaife](https://avatars.discourse-cdn.com/v4/letter/t/22d042/32.png) [@tquaife](https://cms-helpdesk.ncas.ac.uk/u/tquaife)\ **Post date:** [27 September 2023 13:04 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/37 "2023-09-27T13:04:30Z")
Thanks Patrick. This is very informative. I’ll have a think what might be causing it, but we probably need to sit down and look at the code. ---
**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\ **Post date:** [13 November 2023 12:15 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/38 "2023-11-13T12:15:23Z")
Hi Tristan: I have been working on adding the output variable for electron transport, `j`, to the JULES/L2SM code. I am calling this variable `ej` since the CABLE folks have an `ej_max` hook for the `j_max` variable in the code. I noticed when I do this, that the other output variables `gpp` and `gpp_lyr` are greater than 0 at night time in the UK in summer, when the sun isn’t shining for `can_rad_mod=4`. I would have expected that `gpp` be equal to 0. (The option `can_rad_mod=4` doesn’t use sunlit/shaded leaves, BTW.) This doesn’t make sense to me, so I am going through the antecedent variables in the code now. But if that’s what is expected for some reason, please let me know. Patrick ---
**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\ **Post date:** [13 November 2023 13:37 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/39 "2023-11-13T13:37:33Z")
Hi again, Tristan: So, I now have added these output variables to the suite u-da046k\_UK over the UK region: `ej_lyr` for the `j` electron transport, and `tstar` and `tstar_gb` for the leaf temperature (not layered). The non-`gb` variables `tstar` & `ej_lyr` are indexed for PFT. The temperature variables already existed in JULES5.2, they just needed to be added to the suite. The `ej_lyr` variable is = `4.0*gpp_lyr/cmpf_lyr`. This means that `ej_lyr` didn’t really need to be added to the FORTRAN code (but it was), since it could have been computed from `gpp_lyr` and `cmpf_lyr` in post-processing. I have put a sample of 1 year’s of monthly & 6-hourly NETCDF output in `~pmcguire/for_tristan/` on JASMIN. The new variables are in the 2 files with the substring `sif_vars`. You’ll note that the `gpp_lyr` and `gpp` variables are often greater than zero at night in the summer. I still need to sort that out. Patrick ---
**Author:** ![tquaife](https://avatars.discourse-cdn.com/v4/letter/t/22d042/32.png) [@tquaife](https://cms-helpdesk.ncas.ac.uk/u/tquaife)\ **Post date:** [14 November 2023 11:50 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/40 "2023-11-14T11:50:37Z")
Brilliant. Thanks Patrick. As discussed, would also be great to get wlite per layer as an output, but there is no immediate rush on that. [Previous page](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928.md?page=1) [Next page](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928.md?page=3)