# Add L2SM model into JULES

**URL:** <https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928>\
**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:** 1

<div class="post-metadata">

**Author:** ![tquaife](https://avatars.discourse-cdn.com/v4/letter/t/22d042/32.png) [@tquaife](https://cms-helpdesk.ncas.ac.uk/u/tquaife)\
**Post date:** [1 March 2023 13:20 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/1 "2023-03-01T13:20:51Z")

</div>

Request to add the L2SM model into JULES as an optional replacement for the current Sellers scheme. Ideally the new branch should be submitted to the trunk.

The master version of L2SM is here:

> **[GitHub - tquaife/layeredCanopyTwoStream: Radiative transfer code for...](https://github.com/tquaife/layeredCanopyTwoStream)**
>
> Radiative transfer code for vertically inhomogeneous canopies using the Meador and Weavor two stream approximation with adding - GitHub - tquaife/layeredCanopyTwoStream: Radiative transfer code for...

but this is a very generalised form of the model, written in Python, and we (probably Tristan) would need to boil it down into a single Fortran function.

Note - a key task will be deciding how we input vertically varying canopy properties to JULES. Probably, in the first instance, this should be done on a per-PFT basis in the appropriate NML file.

---

<div class="post-metadata">

**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\
**Post date:** [1 March 2023 13:49 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/2 "2023-03-01T13:49:13Z")

</div>

Thanks, Tristan:  
I will get working on this.  
Patrick

---

<div class="post-metadata">

**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\
**Post date:** [16 March 2023 21:31 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/3 "2023-03-16T21:31:35Z")

</div>

Hi Tristan:  
Just a quick check-in.  
I was able to fork your layeredCanopyTwoStream repository, then run both layeredCanopy\_wrapper\_example.py and leafGeometry.py with python2 on JASMIN. The first program produced numbers and the second program produced plots, both without error.  
The checked-out code is now in:  
`/gws/nopw/j04/odanceo/pmcguire/layeredCanopyTwoStream_pcm/layeredCanopyTwoStream`

I guess you want me to translate the function `get_layered_canopy_fluxes( )` into FORTRAN. Let me know if this is not the case. Otherwise, I will try to do that next week.  
Patrick

---

<div class="post-metadata">

**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\
**Post date:** [29 March 2023 13:55 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/4 "2023-03-29T13:55:30Z")

</div>

Hi Tristan:  
I have studied a bit of the JULES radiative-transfer (RT) FORTRAN source code, which uses the Sellers scheme. I am working on JASMIN, currently with the layered-output-canopy branch of JULES 5.2 trunk, which I call `vn5.2_sif1`. The RT code I looked at is here: `~pmcguire/jules/vn5.2_sif1/src/science/radiation/albpft_jls_mod.F90`

I had the idea to try to compare the (Python) code changes between your L2SM model (using the Meador&Weavor two-stream scheme) and the Sellers scheme. This led to me finding another repository of yours:

> **[GitHub - tquaife/pySellersTwoStream: Python implementation of the Sellers...](https://github.com/tquaife/pySellersTwoStream/)**
>
> Python implementation of the Sellers two-stream canopy radiative transfer model - GitHub - tquaife/pySellersTwoStream: Python implementation of the Sellers two-stream canopy radiative transfer model

The `albpft( )` FORTRAN function in `albpft_jls_mod.F90` isn’t a standalone function for vegetation. It also has code in it for urban surfaces, for snow-covered canopies, and for matching observed albedos with `l_albedo_obs=TRUE`. But I guess, we mainly just want an L2SM version of albpft( )?

Patrick

---

<div class="post-metadata">

**Author:** ![system](https://europe1.discourse-cdn.com/flex013/uploads/cms_support/original/1X/1fd2411499ffcbc299fe756cd5cdf26e44956558.png) [@system](https://cms-helpdesk.ncas.ac.uk/u/system)\
**Post date:** [28 April 2023 13:55 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/5 "2023-04-28T13:55:56Z")

</div>

This topic was automatically closed 30 days after the last reply. New replies are no longer allowed.

---

<div class="post-metadata">

**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\
**Post date:** [2 May 2023 12:43 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/6 "2023-05-02T12:43:28Z")

</div>



---

<div class="post-metadata">

**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\
**Post date:** [2 May 2023 12:55 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/7 "2023-05-02T12:55:12Z")

</div>

Hi Tristan & Meg:  
I have been able to add the L2SM model to JULES. It’s a 1st draft version of the code, and I am sure there are mistakes of various sorts. It is a version of the Python code at the GitHub repository tquaife/layeredCanopyTwoStream, which I manually translated from Python to FORTRAN.

The new code is checked into a MOSRS branch at [https://code.metoffice.gov.uk/svn/jules/main/branches/dev/patrickmcguire/r25125\_vn5.2\_layeredCanopyTwoStream](https://code.metoffice.gov.uk/svn/jules/main/branches/dev/patrickmcguire/r25125_vn5.2_layeredCanopyTwoStream)

I modified a GL7 (global land) N96 suite so that it uses this branch of JULES 5.2. The new suite is u-cw591, and it is currently running. It has compiled the code, and it has passed the RECON step, and it is more than 1 year through a 20 year spinup cycle at the moment. Probably the output values will show mistakes of some sorts, but it is running! I added a can\_rad\_mod=7 option which is used to switch from can\_rad\_mod=6 (Sellers) to can\_rad\_mod=7 (L2SM), using the can\_rad\_mod=6 options otherwise.  
Patrick

---

<div class="post-metadata">

**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\
**Post date:** [10 May 2023 09:04 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/8 "2023-05-10T09:04:19Z")

</div>

Hi Tristan & Meg:  
The L2SM GL7 (global land) N96 suite did run on JASMIN for a few years (from 1860-1863), before crashing on 1864 sometime. It turned out that it was running low on memory. I increased the memory available in SLURM on the JASMIN LOTUS batch-processing nodes up to near the normal maximum of 64GB per core, and it has now run for 10 years (1860-1869).

It still has many years to go to reach the 21st century. The LS2M version currently takes about 90 minutes per year, whereas the original Sellers version took about 4 minutes per year.

Also, the LS2M version gives annual-average GPP values in many places that are negative, whereas the original Sellers version gives GPP values that are non-negative everywhere. There certainly is something that needs to be fixed.  
Patrick

---

<div class="post-metadata">

**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\
**Post date:** [25 May 2023 19:18 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/9 "2023-05-25T19:18:18Z")

</div>

Hi Tristan & Meg:  
I started trying to figure out why the GPP carbon-flux values were negative in places. They should always be non-negative.

It looks like some of my code changes for putting the L2SM into JULES were never being reached, when I thought they should be reached. This turns out to be because `l_getprofile` is hard wired as `.false.` in the code, so that the `can_rad_mod == 7` code changes aren’t being reached in the file `~pmcguire/jules/r25125_vn5.2_layeredCanopyTwoStream/src/science/radiation/albpft_jls_mod.F90` .

I am trying out hardwiring it with `l_getprofile=.true.` in `~pmcguire/jules/r25125_vn5.2_layeredCanopyTwoStream/src/science/radiation/tile_albedo_jls_mod.F90`.

There is a comment that says:

```auto
LOGICAL, INTENT(IN) :: l_getprofile
                                   ! TRUE means profiles through the
! canopy are calculated
! ! FALSE means no profiles calculated

```

This l\_getprofile LOGICAL variable is there both in the JULES 5.2 trunk and in the canopy-layering branch  
[https://code.metoffice.gov.uk/svn/jules/main/branches/dev/patrickmcguire/vn5.2\_sif1](https://code.metoffice.gov.uk/svn/jules/main/branches/dev/patrickmcguire/vn5.2_sif1)

The new suite that uses the latest revised version of the L2SM JULES branch is u-cx323.  
It is running right now. We should have some results about if it runs better with `l_getprofile=.true.` maybe by tomorrow.  
Patrick

---

<div class="post-metadata">

**Author:** ![tquaife](https://avatars.discourse-cdn.com/v4/letter/t/22d042/32.png) [@tquaife](https://cms-helpdesk.ncas.ac.uk/u/tquaife)\
**Post date:** [26 May 2023 09:34 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/10 "2023-05-26T09:34:32Z")

</div>

Thanks Patrick.

---

<div class="post-metadata">

**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\
**Post date:** [1 June 2023 19:57 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/11 "2023-06-01T19:57:10Z")

</div>

Hi Tristan:  
It turns out that setting `l_getprofile=.true.` was the wrong thing to do. This gets set to true when the `albpft_jls` subroutine is called in an alternative way.

It turns out that I had a mistake in the `can_rad_mod == 7` logic, so this is why the `can_rad_mod == 7` section wasn’t getting reached in the `albpft_jls` subroutine. I fixed this mistake, and I am rerunning now, so hopefully it works better this time.  
Patrick

---

<div class="post-metadata">

**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\
**Post date:** [5 June 2023 14:28 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/12 "2023-06-05T14:28:00Z")

</div>

Hi Tristan:  
I have been through several more rounds of debugging now. With the `can_rad_mod == 7` logic fixed, and a couple of other bugs fixed, the next bug is that the code fails when `cosz_gb = mu = 0.0`, due to division by zero. I have put in an if statement to skip these cases.

I had some print statements embedded, and prior to adding this `if (mu > 0.01)` statement, the L2SM code in JULES was producing a fair number of NAN’s in the recon & spinup phases of the suite. But when I added this statement, there were no NAN’s in the recon phase; the L2SM code in JULES was producing various float values that didn’t seem unreasonable. However, I ran out of disk space because I was printing things out to the log files, so the recon phase failed. I have now commented out the print statements in the L2SM code in JULES and now, the recon phase succeeds.

Maybe now, the spinup phase will succeed, and the main-run phase of the suite will start producing output that doesn’t have the GPP carbon flux as being \< 0.0.  
Patrick

---

<div class="post-metadata">

**Author:** ![tquaife](https://avatars.discourse-cdn.com/v4/letter/t/22d042/32.png) [@tquaife](https://cms-helpdesk.ncas.ac.uk/u/tquaife)\
**Post date:** [5 June 2023 14:43 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/13 "2023-06-05T14:43:34Z")

</div>

Thanks Patrick.

I assume the Seller’s code must also need to guard against mu=0.0 as it also has to divide by mu. I’ll check the code later, but from memory, it sets a variables K=G/mu.

---

<div class="post-metadata">

**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\
**Post date:** [6 June 2023 11:31 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/14 "2023-06-06T11:31:43Z")

</div>

Hi Tristan  
Yes. Both the Sellers code and the L2SM code need to guard against division by mu=0. I think mu=0 when the sun doesn’t rise, due to the eternal night of winter in the arctic.  
Patrick

---

<div class="post-metadata">

**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\
**Post date:** [6 June 2023 12:18 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/15 "2023-06-06T12:18:42Z")

</div>

Hi Tristan:  
The code now seems to run the L2SM sub-code (unlike before) and get to the main\_run of the JULES suite u-cx323, after spinning up.

For some reason, it previously crashed during spin up for some mysterious reason. I think it was running out of memory; the normal maximum memory for a JASMIN node is 64GB, and I had it set to use the maximum. But for the global JULES suite, it does domain decomposition, normally with 10 processors. And with the par-multi queue setup that we use for this suite, it assigned 3 processors to the same node of 1/10th the single-processor global size, and so I think because there were 3 such JULES processes running on the same node, it ran out of memory. My initial solution is to double the number of processors from 10 to 20, which maybe reduces the memory usage per processor by up to a factor of 2. I tried this, and this time it is running without crashing mysteriously.

But the GPP is still often \< 0, so I still need to work on that. I will try to do some print statements for after spinup, within the L2SM code.  
Patrick

---

<div class="post-metadata">

**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\
**Post date:** [6 June 2023 12:50 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/16 "2023-06-06T12:50:43Z")

</div>

Hi Tristan:  
I suspect one reason that GPP\<0 with the L2SM code in JULES, is that I don’t currently return all the output variables from the albpft routine:

Currently, there are values returned in `can_rad_mod == 5` normal non-L2SM usage:

```auto
            rnet_dir(i) = rdird - rdiru
            rnet_dif(i) = rdifd - rdifu

            fapar_dir(l,n,i) = (rnet_dir(i-1) - rnet_dir(i)) / dlai
            fapar_dif(l,n,i) = (rnet_dif(i-1) - rnet_dif(i)) / dlai

            fapar_dir2dir(l,n,i) = (1.0 - om) * dabeer_dla(i)
            fapar_dir2dif(l,n,i) = om * dabeer_dla(i) + fapar_dir(l,n,i)
            fapar_dif2dif(l,n,i) = fapar_dif(l,n,i)

            fsun(l,n,i) = exp_minus_kla * (EXP(k * dlai) - 1.0) / (k * dlai)

            fapar_dir(l,n,i) = fapar_dir(l,n,i) * can_struct_a(n)
            fapar_dif(l,n,i) = fapar_dif(l,n,i) * can_struct_a(n)
            fapar_dir2dir(l,n,i) = fapar_dir2dir(l,n,i) * can_struct_a(n)
            fapar_dir2dif(l,n,i) = fapar_dir2dif(l,n,i) * can_struct_a(n)
            fapar_dif2dif(l,n,i) = fapar_dif2dif(l,n,i) * can_struct_a(n)

```

But in the L2SM version with can\_rad\_mod == 7, I currently only return a couple of the variables (fapar\_dir & fapar\_dif). I need to go through the math and return the other variables besides those 2).

```auto
            !Original L2SM Python code from Tristan Quaife:
            !Idn=self.layers[layer].Idn_dir
            !Iup=self.layers[layer].Iup_dir
            !self.layers[layer].Iab_dir=np.abs(Iup_below-Iup+Idn_above-Idn)
            !self.layers[layer].Iab_dir=np.abs(Iup_below-Iup+Idn_above-Idn)
            !Idn=self.layers[layer].Idn_dif
            !Iup=self.layers[layer].Iup_dif
            !self.layers[layer].Iab_dif=np.abs(Iup_below-Iup+Idn_above-Idn)

            !initial L2SM FORTRAN code using JULES variable names, by Patrick McGuire
            rdn = rdirdn(i)
            rup = rdirup(i)
            fapar_dir(l,n,i) = ABS(rdirup_below-rup+rdirdn_above-rdn)
            rdn = rdifdn(i)
            rup = rdifup(i)
            fapar_dif(l,n,i) = ABS(rdifup_below-rup+rdifdn_above-rdn)

```

Patrick

---

<div class="post-metadata">

**Author:** ![tquaife](https://avatars.discourse-cdn.com/v4/letter/t/22d042/32.png) [@tquaife](https://cms-helpdesk.ncas.ac.uk/u/tquaife)\
**Post date:** [7 June 2023 08:55 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/17 "2023-06-07T08:55:50Z")

</div>

Hi Patrick,

It might be quicker for me to take you through what is required here. You should be able to pull those calculations from other parts of the code _but_ when we have vertical variability in the code, they will need to be slightly different.

The basis of the equations JULES uses is given in this paper:

> **[A Two-Big-Leaf Model for Canopy Temperature, Photosynthesis, and Stomatal...](https://journals.ametsoc.org/view/journals/clim/17/12/1520-0442_2004_017_2281_atmfct_2.0.co_2.xml?tab_body=fulltext-display)**
>
> Abstract The energy exchange, evapotranspiration, and carbon exchange by plant canopies depend on leaf stomatal control. The treatment of this control has been required by land components of climate and carbon models. Physiological models can be used...

Thanks,  
Tristan.

---

<div class="post-metadata">

**Author:** ![pmcguire](https://avatars.discourse-cdn.com/v4/letter/p/b782af/32.png) [@pmcguire](https://cms-helpdesk.ncas.ac.uk/u/pmcguire)\
**Post date:** [9 June 2023 14:44 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/18 "2023-06-09T14:44:31Z")

</div>

Hi Tristan:  
Thanks.  
For the can\_rad\_mod == 7 version (L2SM), I did lift the equations for fsun, fapar\_dir2dir, fapar\_dir2dif, etc. verbatim from the can\_rad\_mod == 5 version (as quoted in my previous post here). I assumed that can\_struct\_a(NPFT) == 1.0, as it is in the GL7 suite, which can be revisited if necessary. The new suite (suite number u-cx605) compiles, and is starting to run now.  
Patrick

---

<div class="post-metadata">

**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 June 2023 21:23 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/19 "2023-06-13T21:23:23Z")

</div>

Hi Tristan:  
The suite `u-cx605` did start running, but it soon failed, in the RECON step. I had some trouble tracing this down. Using `jules_print()` statements wasn’t enough, so I also switched from `JULES_BUILD=normal` to `JULES_BUILD=debug`. This allowed the logs to show that in the `albpft( )` code, I was addressing one of the L2SM arrays with a 0 index when the lowest index was actually 1.

I also found a division-by-zero-snowdepth error in an unrelated file, `relayersnow_mod.F90`, which I have temporarily disabled.

And I found that my setting of a lower-limit on `mu = cos_solar_zenith >= 0.01` was allowing overflows, in the `albpft( )` code, with `tau/mu >= 38` and `EXP(tau/mu)` overflowing in some cases. So I changed the lower limit of `mu` to 0.05. This allowed the code to get through the RECON step.

Now, it gets into the SPINUP step, but there are multiple different errors on the different processors after running for up to 6 months in the year 1860. I am currently using 40 processors (for domain decomposition), increasing first from 10 to 20 to avoid running out of memory and then from 20 to 40 to try to get rid of the aforementioned overflow. But making it to the SPINUP step is progress, as I can now look at the end-dump from the RECON step, to make sure the numbers appear legitimate.

I do note that with `JULES_BUILD=debug`, when it makes it through the `RECON` step, it fails with overflow in the `compactsnow( )` routine, and possibly other places. So I am now running with the maximum number of snow layers changed from 3 to 0 in the rose cylc suite `nsmax` variable, in order to avoid running the `compactsnow( )` routine’s overflow.

With these settings, in the SPINUP step, there is another overflow in the `albpft( )` code at the same spot. So I changed the minimum `mu=cos_solar_zenith` further, from 0.05 to 0.10. Now it can run through at least 2 spinup years for each of the 40 processors running on the subdomains. Previously, it was limited to no more than 6 to 10 months of the 1st spinup year (1860), prior to overflowing. It should start producing output after 10 spinup years. Maybe the output now has the `GPP carbon flux >=0.0`, like it should.

Patrick

---

<div class="post-metadata">

**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 14:28 UTC](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928/20 "2023-06-14T14:28:11Z")

</div>

Thanks Patrick.

Long term, I am not sure that having such a big threshold on mu is appropriate. It is probably better to work out what is going on. The L2SM model itself should be fine with any value \>0.0 and hence the numbers it returns to the rest of JULES should be OK.

Tristan.

[Next page](https://cms-helpdesk.ncas.ac.uk/t/add-l2sm-model-into-jules/928.md?page=2)