bug in spitter call (NESP2pt9 and TRENDY_POP branches)

[har917]

At some point in the history of CABLE development the allocation of downwelling shortwave radiation into direct and diffuse components has gained bug for offline studies - this likely occurred at the time when the code considered the fractionation between visible and near-infrared radiation.

The direct-diffuse radiation (unless provided via other means) is evaluated through the use of the spitter function. Originally this assigned the fractions based on the ratio of the timestep's shortwave radiation to a clear sky value (at that location and time of year). When the partitioning between VIS:NIR occurred spitter was not adjusted and currently the call takes in the current VIS (or NIR) shortwave and compares to the clear sky value for total shortwave radiation (not the clear sky value for the VIS (or NIR) shortwave). The ratio evaluated is then smaller than expected - and the diffuse fraction diagnosed is greater than it would be.

If only for testing purposes - we need to retrofit the model to allow for the original intent of spitter to be used.

This implementation will need to be different in the current MAIN (CABLE3) and in the CABLE-POP-TRENDY/NESP2pt9/BLAZE branches because of the code structure.

see also #355

[har917]

@ccarouge There are a couple of up front question we need answer as part of this issue

• as a bug fix - should this go in on a switch?
• how do we manage the code structure difference between MAIN and the older branches?

[har917]

For CABLE-POP branch the suggestion is to work in the cable_radiation.f90, subroutine init_radiation.

In the short-term we can a logical switch - at around line 64 add
LOGICAL :: use_spitter_bugfix !to become a namelist option

The default value of the new switch would be (apply the bug fix - which is a change) via line 67 add
use_spitter_bugfix = .TRUE.

Finally at the calls to spitter (line 132) is then replaced

! Define beam fraction, fbeam:
IF (use_spitter_bugfix) THEN   
      !beam fraction defined on total shortwave, equal in VIS and NIR
      rad%fbeam(:,1) = spitter(met%doy, met%coszen, met%fsd(:,1)+met%fsd(:,2))
      rad%fbeam(:,2) = spitter(met%doy, met%coszen, met%fsd(:,1)+met%fsd(:,2))
ELSE 
    !CABLE2.x scheme
     rad%fbeam(:,1) = spitter(met%doy, met%coszen, met%fsd(:,1))
     rad%fbeam(:,2) = spitter(met%doy, met%coszen, met%fsd(:,2))
END IF

In the longer-term we would look to move the logic switch into a cable_users%namelist

edit: for efficiency it's probably better to do

! Define beam fraction, fbeam:
IF (use_spitter_bugfix) THEN   
      !beam fraction defined on total shortwave, equal in VIS and NIR
      rad%fbeam(:,1) = spitter(met%doy, met%coszen, met%fsd(:,1)+met%fsd(:,2))
      rad%fbeam(:,2) = rad%fbeam(:,1)
ELSE 
    !CABLE2.x scheme
     rad%fbeam(:,1) = spitter(met%doy, met%coszen, met%fsd(:,1))
     rad%fbeam(:,2) = spitter(met%doy, met%coszen, met%fsd(:,2))
END IF

[mcuntz]

@har917 Why should a bugfix have a switch? If it is a bug, then it should simply be corrected, isn't it? Why would I want to calculate the wrong code?

[har917]

@mcuntz Backwards compatibility/reproducibility for at least the short term. We're also going to assess the impact of this on GPP, NPP, equilibrium carbon pools etc. - it's highly likely that (since GPP has been one of the metrics we, recently meaning @juergenknauer, tune against) we should pair this bug fix with a recalibration of the model parameters.

[juergenknauer]

Well done on finding this bug @har917. Backward compatibility sounds reasonable to me in the short term (with the default flag being the corrected version).

Most likely we'd have to retune GPP. The extkn parameter (canopy N extinction coefficient) would be an obvious candidate.

[har917]

@rml599gh did the hard work not me. All part of the development of an offline-as-ACCESS capability under the CABLE4/ESM3 project.

[har917]

@AlisonBennett @juergenknauer @mcuntz - I think we pause this until the parallel work in #355 is completed. That will give us an initial idea as to how important the issue potentially is and whether a recalibration of the parameters is necessary. Working into MAIN will allow the use of benchcab and me.org.

[AlisonBennett]

I have performed a test on the reccap1000pts domain using the BLAZE_9184 branch with and without this fix. See attached some plots of the main carbon pools and stocks from the cable output. The black line (AGCD_1000pt_354_spitter_bugfix) is with the bug fix implemented, the orange line (AGCD_1000pts_BLAZE_off) is without the bug fix implemented.

[rml599gh]

Here are results from one of my tests using the main branch. Blue is with the fix, red without and both are compared to GPP from an ACCESS-ESM1.5 run. Plots are daily GPP for two grid-cells, one dominated by evergreen broadleaf (71.25W, 6.25S) and one by deciduous broadleaf (61.9W, 20S). Differences look to be larger for larger daily GPP - probably not surprising.

[har917]

see #357 Take home message: Applying this bug fix to MAIN without recalibration results in minor/modest impacts (mainly improvements) in the short-term biophysics.

The results seen in the biogeochemistry above likely reflect the accumulation of these small differences through the spin up.

[ccarouge]

To decide on a switch or not for this fix, I need more information. I have discussed this with various people and followed this discussion. For me, a switch would only make sense if we say the "fixed" version isn't yet ready for production runs. This means the default option should use the "original" (buggy) version that has been used so far for production runs. The switch enables the "fix" to be put in the code and continue work on to get a better set of parameters for it.

Once we have a full fix, i.e. code + parameters, we remove the switch and the fixed version becomes the production version.

Does that work for everyone @har917 @mcuntz @juergenknauer @rml599gh ? If you are going to use the fixed implementation for published runs now, I would not put it behind a switch.

[mcuntz]

My twopence worth:
From the plots by @har917 @rml599gh and @AlisonBennett, the fix does not seem to change the results very much (which is good in some respects).
From my experience, CABLE is much further off the observations than the changes by the fix. Global tuning of @juergenknauer is also not that precise but wants to get the ballpark right.
So I will just use the fixed version in the future.

[juergenknauer]

For my part it's not necessary to put it behind a switch. I'm happy to switch to the corrected version instantly.

Based on the changes shown in @AlisonBennett 's comment above it is also unnecessary to retune parameters. The changes are well within the observational uncertainty. But keep in mind that I always compared with the global evaluation datasets only, not specifically for Australia.