A small bug in ana_srflux.h

[takashi]

Hello all,
I'm now using ALBEDO function.
I found a small bug in line 175 of the file ana_srflux.h.
The original version is:

Code: Select all

          zenith=cff1+cff2*COS(Hangle-lonr(i,j)*deg2rad/15.0_r8)

I think the lonr should not be divided by 15.

Code: Select all

          zenith=cff1+cff2*COS(Hangle-lonr(i,j)*deg2rad)

Timings of sunrise and sunset had been incorrect, but the local diurnal daylight pattern has been well reconstructed after this modification.

[kate]

Great, thanks for the fix!

[wilkin]

Wait! I believe the change you made is only justified if in your model the hour that defines Hangle is local time.

As coded it should be correct for a model assuming Universal Time. The hour angle in UT (or GMT) must depend on longitude:
http://en.wikipedia.org/wiki/Hour_angle

Users assuming UT or GMT will introduce a bug if they make this change.

[takashi]

Thank you for your reply.
I found this bug when I use UT (or GMT).
I’ve not removed the longitude (just removed the 15) by this modification.
I think this modification is justified when we use UT (or GMT).
But as you mentioned about the compatibility between UT and local time, the compatibility problem is still remaining when we use local time.

One solution for this problem is, for example,
(1) to define the new CPP flag in the project header file (XXXX.h) as following, if we use local time:

Code: Select all

# define LOCAL_TIME +9.0

(+9.0 is Japan time case, this value is the difference between local time and UT),

(2) to change the line 149 of ana_srflux.h:

Code: Select all

      Hangle=(12.0_r8-hour)*pi/12.0_r8

to

Code: Select all

# ifdef LOCAL_TIME
      Hangle=(12.0_r8-hour+(LOCAL_TIME))*pi/12.0_r8
# else
      Hangle=(12.0_r8-hour)*pi/12.0_r8
# endif

I have tested this modified code under Cygwin+gfortran, and it has been run with no problem.

[mitya]

the dvision by 15 is included into deg2rad already:

a locally corrected Hangle is:
(12.0_r8-hour-(lonr(i,j)/15.0_r8))*pi/12.0_r8 =
Hangle – lonr(i,j)*pi/180 = Hangle-lon*deg2rad

deg2rad = pi/180

Wait! I believe the change you made is only justified if in your model the hour that defines Hangle is local time.

As coded it should be correct for a model assuming Universal Time. The hour angle in UT (or GMT) must depend on longitude:
http://en.wikipedia.org/wiki/Hour_angle

Users assuming UT or GMT will introduce a bug if they make this change.

[lcbernardo]

Dear ROMS users,

I'd like to ask if this issue regarding the possible bug in ana_srflux.h has already been resolved. I ask because I'm using ROMS code updated up to circa 2018, and I find that if I convert the model outputs for solar radiation flux and surface net heat flux to local time, the peaks are occurring close to midnight. My model domain is in the subtropical northern hemisphere so it doesn't make physical sense.

Thanks,
Lawrence

[pwallhead]

For me the division by 15.0 on line 176 is a definite bug:
1) It breaks periodicity over 360 degrees change in lonr
2) It is inconsistent with lines 205 and 213
3) As mitya points out, the division is in effect already accounted for in deg2rad

Changes to account for local time definitions should occur where "hour" enters the calculation of Hangle, as takashi suggests. These will add a constant phase offset to the argument of cos(). The division by 15.0 on line 176 seems to be wrong in any case.

One other thing. I understand that the Rutgers code assumes input of NET downward shortwave irradiance as srflx, since there are no internal corrections for ocean albedo (unlike the Hedstrom code, which assumes input of downwelling irradiance). Therefore, shouldn't the Rutgers version of ana_srflux.h include a correction for ocean albedo (at least the simple 6% reduction)? From Niemela et al. (2001) it seems that the Zillman (1972) parameterization used in ana_srflux.h is for the downwelling irradiance, not the net downward.

Best regards
Phil

[wilkin]

Can you clarify which lines you are referring to in which version of the code:

For me the division by 15.0 on line 176 is a definite bug:
1) It breaks periodicity over 360 degrees change in lonr
2) It is inconsistent with lines 205 and 213

They don't align with the version of ana_srflux.h I'm looking at from myroms.org.

By net downward radiation we mean the NET (downward minus upward reflected) entering the ocean i.e. incident upon the ocean ocean surface in the downward direction. Perhaps net downward is confusing, but net radiation leaves the sign convention ambiguous.

So it is assumed this is given with the effect of ocean albedo included. I compute this from WRF dswrfsfc minus uswrfsfc. In most instances of WRF met model output I see ocean uswrfsfc = 0.06 dswrfsfc as you note.

[arango]

Both ana_srflux.h and ana_specir.h were coded longtime ago when we were testing the EcoSim bio-optical model in the late 90's. It was like twenty years ago. So I don't remember the details. At that time, both routines were intended for idealized toy problems for debugging purposes. The EcoSim model is quite complex and may have up to 60 spectral bands, and we needed surface solar downwelling spectral irradiance. One needs to be very precise in the terminology when dealing with bio-optical models.

In practical applications, we get the shortwave, and longwave radiation fluxes from atmospheric models and the data need to resolve the daily cycle. The turbulent boundary layer between atmosphere and ocean is complicated with various spatial and temporal scales (see CBLAST cartoon below), and the upward and downward fluxes depend if you are looking in terms of the atmosphere or the ocean. It is due to the monotonic ascending (atmosphere) or descending (ocean) discretized vertical coordinate.

I want to be very clear about the nomenclature in the ocean (ROMS) about surface fluxes:

Heat Flux: Downward Flux = heating (positive); Upward Flux = cooling (negative)
E - P Flux: Downward Flux = salting (positive); Upward Flux = freshening (negative)

The heat flux in the atmosphere is reversed and the moisture flux is actually P-E. I imagine lots of ROMS users processing the atmospheric forcing incorrectly. It will take a long time in the simulation to feel the wrong direction of the flux. Recall that these fluxes are applied as surface conditions to the vertical diffusion equations, and its kinematic values are small (around E-05). I highly recommend using the downward fluxes for shortwave and longwave radiation and make the corrections to the shortwave due to the mean absorption in the cool-skin layer, and removing the outgoing infrared radiation from the longwave due to the sea surface temperature (-emiss*StBolt*SST^4; where the SST is in Kelvin). Now, the surface fluxes get more complicated when we have sea ice. It is tricky, if not impossible, for ROMS to check these fluxes and figure out what the user provided. It will require a rigorous metadata model, which most everybody ignores. Recall that the user controls the metadata that goes into ROMS via the customizable input file varinfo.dat.

Now, ROMS clock was reworked recently (see trac ticket 724). The time in ROMS is always with respect to the Greenwich Mean Time (GMT), also known as Coordinated Universal Time (UTC). If you are using any other time zone, you will likely get into trouble when dealing with atmospheric forcing, tidal forcing, coupled Earth Modeling Systems (ESMs), analysis against observations, data assimilation, and any other application with forcing sensitive to the time clock.

Therefore, I haven't reread the literature about the shortwave to check if we have a bug or not. But one needs to be aware of the time clock and default GMT time zone in ROMS. It will be insane using any other time zone in a ROMS application.

[pwallhead]

Hi

Strange that our line numbers don't coincide - I refer to the Rutgers code, svn-updated 05/03/2019, and the ana_srflux.h is from ROMS/Functionals/. Anyway here is line 176:
zenith=cff1+cff2*COS(Hangle-lonr(i,j)*deg2rad/15.0_r8)
and lines 205, 213:
& (cff1+cff2*COS(Hangle-lonr(i,j)*deg2rad)))
& (cff1+cff2*COS(Hangle-lonr(i,j)*deg2rad)))

I fully agree that it is better to use atmospheric model outputs to set srflx (and usually I do). However, you should be aware that many groups do in fact use the analytical formula approach (ANA_SRFLUX+ALBEDO), whether as a provisional, rough solution, and I think in some cases as a final solution. If you believe that this approach is not usable under any circumstances, you should probably remove it from your code. I personally think it can provide a reasonable approximation, if you remove the erroneous division by 15.0 and add a 6% correction for ocean albedo. This latter should come immediately after the IF END IF within the # if defined ALBEDO section (the 6% correction should not be applied if you have DIURNAL_SRFLUX because input of net downward radiation is assumed in Rutgers ROMS --- the comments could be improved to clarify this).

Phil

[arango]

Finally, I have time to browse the literature, Googling, and write a Matlab script to check the issue mentioned above. Indeed, the 15.0 factor used to compute the zenith is incorrect We need to have instead:

Code: Select all

        zenith=cff1+cff2*COS(Hangle-lonr(i,j)*deg2rad)

In my Matlab script, I compared against a WRF-ROMS coupled solution for the US West Coast at a particular point in the grid:

Here, the X-axis is the time since initialization and it is in GMT time zone. We need to subtract 8 hours for local Pacific Standard Time (PST) to analyze the results. As you can see, the WRF-ROMS coupled solution values agree well with the red curve that computes zenith with updated above equation. I am now satisfied that I proved to myself that it is a bug.

Many thanks to all that brought this bug to my attention.

[mitya]

sorry for not insisting back then, I had quite good visual proof, as our domain was huge over all Arctic. I benchmarked "ecodynamo" (by Pedro Duarte from the Norwegian polar instutute) code which used own swrad calculation with ROMS' Powell model using ana_srflx from ROMS. The match became ok when I corrected ana_srflx the way I described above:
https://source.uit.no/mitya/results_sha ... ynamo_roms

[lcbernardo]

Thank you for the bug fix. With this modification, the peaks in swrad in my application now occur at the expected times.