Add equation of time correction to local time for sw radiation calculation

[weiwangncar]

This PR adds equation of time correction to the local time for computing shortwave radiation. The EOT adds a max of 16 min in Nov, and a min of -14 min in Feb. This correction has been added to WRF since 2014.

The PR improves accuracy of shortwave radiation, but relatively small overall effects are expected.

[mgduda]

@weiwangncar It seems like there are quite a few magic numbers in the new code:

       da=6.2831853071795862*(julian-1.)/365.
       eot=(0.000075+0.001868*cos(da)-0.032077*sin(da) &
            -0.014615*cos(2*da)-0.04089*sin(2*da))*(229.18)
       xt24 = mod(xtime+radt*0.5,1440.)+eot

Could some of these be defined as parameters with descriptive names? For example, it seems like 6.2831853071795862 could be twopi or 2.0 * pi. I have no idea what 0.032077, 0.001868, etc. represent.

[weiwangncar]

@mgduda I will see if we can find the original contributor on the WRF side.

[dudhia]

Yes, that is 2*pi. The other numbers are a fit to an annual cycle. They describe the phase.

[mgduda]

Yes, that is 2*pi. The other numbers are a fit to an annual cycle. They describe the phase.

If we can't create descriptively named parameters for the numbers, would it at least be possible to lay out the equation with an explanation of the parameters in a comment block?

[dudhia]

It's a simple harmonic function fit based on the day of the year. The numbers are empirical fit coefficients.

[weiwangncar]

Will revisit at a later time.

[weiwangncar]

Here is a reference for the equation used, the second equation on the first page. Would this help with this PR?