Add L2SM model into JULES

JULES
21

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

22

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.

23

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?

24

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.0k(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

25

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

26

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.

27

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:

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

28

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.

29

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

30

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.

31

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

!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
32

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).

33

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
cfplot_gpppftlyr_v7i1e1170×830 108 KB

cfplot_aparpftlyr_v7i1f
cfplot_aparpftlyr_v7i1f1170×830 93.3 KB

34

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?

35

Thanks, Tristan!
I will answer you privately about meeting up.
Patrick

36

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
cfplot_gpppftlyr_v7i1g1170×830 104 KB

GPP can_rad_mod=7 L2SM/Meador&Weaver
cfplot_gpppftlyr_v7i1h
cfplot_gpppftlyr_v7i1h1170×830 104 KB

APAR can_rad_mod=6 Sellers

cfplot_aparpftlyr_v7i1g
cfplot_aparpftlyr_v7i1g1170×830 84.1 KB

APAR can_rad_mod=7 L2SM/Meador&Weaver
cfplot_aparpftlyr_v7i1h
cfplot_aparpftlyr_v7i1h1170×830 92.9 KB

37

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.

38

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

39

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

1 Like
40

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.

1 Like