Cosine solar zenith angle forcing is off by a time-step.

[ekluzek]

From some analysis by Sean Swenson, we see that the cosine solar zenith angle forcing is off by a time-step. See the PPT file attached.

I think the problem is in shr_tInterpMod.F90 in shr_tInterp_getAvgCosz, where a time update happens and then the solar-zenith angle is calculated. I think that should be reversed the solar zenith angle is calculated and then the time-update happens.

This affects CLM standalone simulations for all forcing types (I compsets). But, also CORE_IAF and CORE2_IAF forcing (JRA and IAF C and G compsets as well as the DIAF compset).

We want this to go on the maint-5.6 branch as well as master. But, do NOT expect this for the cesm2.1.2 release that is imminent. We do need to do some simulations with it to make sure it doesn't do anything unexpected. We think the climate signal should be minimal.

coszen_check.pdf

[ekluzek]

@jgfouca and @rljacob you should look into this issue for E3SM. Let us know if you want the e3sm version to preserve current answers.

[ekluzek]

@alperaltuntas you should consider this problem for the ocean compsets that it will impact.

[rljacob]

This has no impact on B or F cases?

[ekluzek]

This has no impact on B or F cases?

Yes correct.

[ekluzek]

OK Sean Swenson (@swensosc) has done more work on this. It looks like lower is doing the wrong thing, while upper and linear are correct. nearest is harder to judge.

Looking at the code, it looks like both lower and upper are wrong as they compare itime1 and time2 (the lower and upper bounds) rather than for the current time. We think upper would be correct for all but time 0, and lower would be wrong when itimein == itime2. This code as it is goes back at least 15 years. But, in general lower and upper aren't used, and there are no unit tests for them.

Sean is going to try out the proposed code changes and we'll see if we get the results we expect.

[rljacob]

@jonbob see above regarding G-cases.

[jonbob]

@rljacob - yes, I had been following this PR. I'm guessing we should wait and see what they decide after their test...

[olyson]

Sean asked me to review this bug and the proposed changes. I do agree that it is a bug and I agree with his fix. I've run an IHIST case with the bug fix and compared to a control and ran diagnostics on it here:

http://webext.cgd.ucar.edu/I20TR/clm50_ctsm10d084_2deg_GSWP3V1_shrinterp_hist/lnd/clm50_ctsm10d084_2deg_GSWP3V1_shrinterp_hist.1995_2014-clm50_ctsm10d084_2deg_GSWP3V1_hist.1995_2014/setsIndex.html

I don't see any significant changes in key variables. More than roundoff but not climate-changing. In terms of atmospheric forcing, it's just incoming solar radiation that changes. All other forcing fields appear to be bfb.

I'm not sure how changes to cime are vetted or implemented, but Sean and I are recommending that this be fixed. The code changes can be found here:

From /glade/work/oleson/ctsm_shrinterp/cime: git status HEAD detached at branch_tags/cime5.8.15_a01 modified: src/share/streams/shr_dmodel_mod.F90 modified: src/share/streams/shr_tInterp_mod.F90

[ekluzek]

Here's a note from @swensosc on what he found...

my code analysis showed a couple of possible bugs in the time interpolation code. First, the ‘upper’ and ‘lower’ conditionals appeared to never be triggered. Second, the calculation of the average coszen value for the forcing period advanced the time value used in the coszen calculation before the calculation occurred, rather than afterwards. Finally, the conditional used to determine when the forcing data should be advanced was delayed by one time step. (This may be why the coszen average subroutine was advanced. By advancing prior to calculating the coszen, it matched the way the forcing was advanced and therefore gave the correct value of forcing for the time period).

The upper/lower conditions are found in:

shr_tInterp_getFactors (shr_tInterp_mod.F90)

I think the following conditional is not needed:

if (itime1 < itime2) then

The average coszen value is calculated in:

shr_tInterp_getAvgCosz (shr_tInterp_mod.F90)

I think these statements:

call shr_cal_advDateInt(ldt,'seconds',ymd0,tod0,ymd,tod,calendar)

call shr_cal_timeSet(reday,ymd,tod,calendar)

should be moved after:

call shr_tInterp_getCosz(cosz,lonr,latr,ymd,tod,eccen,mvelpp,lambm0,obliqr,calendar)

n = n + ldt

tavCosz = tavCosz + cosz*real(ldt,SHR_KIND_R8) ! add to partial sum

The determination of when to advance the forcing data is in:

shr_dmodel_readLBUB (shr_dmodel_mod.F90)

I think the ‘>’ should be ‘>=’, i.e.:

old:

if (rDateM < rDateLB .or. rDateM > rDateUB) then

new:

if (rDateM < rDateLB .or. rDateM >= rDateUB) then

[github-actions[bot]]

This issue is stale because it has been open 30 days with no activity. Remove stale label or comment or this will be closed in 5 days.

[github-actions[bot]]

This issue was closed because it has been stalled for 5 days with no activity.