Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Local Solar Time Computation #323

Merged
merged 7 commits into from
Aug 30, 2019
Merged

Local Solar Time Computation #323

merged 7 commits into from
Aug 30, 2019

Conversation

dabail10
Copy link
Contributor

Fix the local solar time calculation as highlighted by Pedro Duarte in Issue #3. Basically this is a couple lines as follows:

! solar_time = mod(real(sec,kind=dbl_kind),secday)/c3600 &
! + c12sin(p5TLON(i,j))
solar_time = mod(real(sec,kind=dbl_kind),secday)/c3600 &
+ TLON(i,j)/(15._dbl_kind*deg2rad)
if (solar_time .ge. 24._dbl_kind) solar_time = solar_time - 24._dbl_kind

I cannot find any documentation about the original sinusoidal calculation here. This change only impacts the runs which use the Large-Yeager CORE forcing as in AOMIP. Currently this is the case in gx1 runs. Hence this only changes answers for gx1 cases in the test suite. Here are some plots that show the before and after computation of the incoming shortwave.

fswdn_old
fswdn_new

  • Developer(s): D. Bailey

  • Please suggest code Pull Request reviewers in the column at right.

  • Are the code changes bit for bit, different at roundoff level, or more substantial? This changes answers for all cases that use the LY CORE forcing.

  • Please include the link to test results or paste the summary block from the bottom of the testing output below.

FAIL cheyenne_intel_restart_gx1_40x4_droundrobin_medium compare cice6.1 21.27 3.70 6.03 different-data
FAIL cheyenne_intel_restart_gx1_4x2_bgcsklclim_medium compare cice6.1 148.87 26.58 71.15 different-data
FAIL cheyenne_intel_restart_gx1_8x1_bgczclim_medium compare cice6.1 348.01 16.23 176.79 different-data

  • Does this PR create or have dependencies on Icepack or any other models? N

  • Is the documentation being updated with this PR? (Y/N) N
    If not, does the documentation need to be updated separately at a later time? (Y/N) N?

Note: "Documentation" includes information on the wiki and .rst files in doc/source/,
which are used to create the online technical docs at https://readthedocs.org/projects/cice-consortium-cice/.

  • Other Relevant Details:

@dabail10 dabail10 changed the title Solar Local Solar Time Computation Jun 13, 2019
Copy link
Contributor

@eclare108213 eclare108213 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This appears to be written specifically assuming TLON > 0. Please check whether that condition is enforced all the time for the code, or if sometimes -180 < TLON < 180. (Sorry I don't remember!) It would probably be better to code this so that TLON can be defined either way.

@dabail10
Copy link
Contributor Author

Good point. I was actually assuming it was between -180 and 180.

Copy link
Contributor

@apcraig apcraig left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch on TLON @eclare108213. I'll approve the PR assuming that gets sorted out.

@dabail10
Copy link
Contributor Author

I am still working on this, although cheyenne is down this week. The longitude thing is kind of a can of worms. I understand now why someone might have chosen the sine function. For example, what if longitude goes from 180 to 540? This is something they are doing with the MOM domains. Thoughts?

@apcraig
Copy link
Contributor

apcraig commented Jun 24, 2019

I think it's generally dangerous to assume anything about the lon value. But one option might be to use the mod function to turn any lon into an equivalent lon in a different range. I think we could assume something like the lon will never be less than -3600 or something then do

mod (lon+3600,360) to make it 0-360
mod (lon+180+3600,360)-180 to make it -180 to 180

etc. If you don't want to deal with a hardwired "3600", you can always do

do while (lon < 0)
lon = lon + 360
enddo

then use the mod function to shift into your range of interest. we could even have a new array in CICE that would be TLON_bounded that would convert TLON into whatever lon range we wanted to define for TLON_bounded. Then in place of TLON, the code could use TLON_bounded. We could also convert TLON after reading to be in the bounds that CICE decides it wants. In that sense, users could define TLON with any values, but the CICE could would covert to 0_to_360 or -180_to_180 or whatever.

@eclare108213
Copy link
Contributor

It seems like the code already has some of this kind of thing in it. Maybe for history, or input forcing. Grep...

@dabail10
Copy link
Contributor Author

Good thinking. There is a line in ice_history_write.F90 like this:

workr2 = mod(tlon(:,:,1:nblocks)*rad_to_deg + c360, c360)

You learn something new everyday.

@dabail10
Copy link
Contributor Author

Ok. I thought about Tony's suggestions above. I think this will account for all reasonable longitude values.

! solar_time = mod(real(sec,kind=dbl_kind),secday)/c3600 &
! + c12sin(p5TLON(i,j))
solar_time = mod(real(sec,kind=dbl_kind),secday)/c3600 &
+ (mod(TLON(i,j)+c180+c3600,c360)-c180)/c15
if (solar_time .ge. 24._dbl_kind) solar_time = solar_time - 24._dbl_kind

I am running the test suite now and will commit later today.

@dabail10
Copy link
Contributor Author

Sorry to take so long on this. Much more complicated than I thought. Here is the latest version that I believe will work for all possible longitudes. First I use mod to convert longitudes to a range between -360 to 360. Then I shift the longitudes to -180 to 180, so I get +/-12 hours from GMT.

! Convert longitude to range of -180 to 180 for LST calculation.

    lontmp = mod(TLON(i,j)/deg2rad,c360)
    if (lontmp .gt. c180) lontmp = lontmp - c360
    if (lontmp .lt. -c180) lontmp = lontmp + c360

    solar_time = mod(real(sec,kind=dbl_kind),secday)/c3600 &
               + lontmp/c15

@dabail10
Copy link
Contributor Author

Just committed the new changes. @eclare108213 let me know if this is what you expected.

@apcraig
Copy link
Contributor

apcraig commented Jul 19, 2019

I believe this implementation should work, although I tend to avoid mods of negative numbers as it tends to confuse people. I am fine with this, but will let others review and comment as well.

@eclare108213
Copy link
Contributor

I am okay with this solution, although I'd like it better (for efficiency) if it just used mathematical functions instead of 'if' conditionals. Does it not work to do something like this (with suitable conversion factors, etc):

lontmp = mod(TLON+360,360)-180

? Maybe not -- I'm not fully aware of all the complications you're trying to work around.

It would be good to include this PR in the release. To do so, we'll need to do the runs to replace whatever sample output we have for the gx1 CORE-forced case. It would also be useful to run the QC tests -- I'd expect this to change the answers significantly but not change the sea ice 'climate'. If these runs can't be done quickly, maybe this waits for the next release.

@dabail10
Copy link
Contributor Author

My main concern was if the longitude was less than -360 or greater than 360. Would a longitude of -720 be completely ridiculous? Sure, but I wanted to make sure it was covered.

@dabail10
Copy link
Contributor Author

I'll leave the decision up to Tony when to merge this.

@eclare108213
Copy link
Contributor

If the longitude is outside the (-360,360) range, the mod function will bring it back in, no matter how large (or negative) the initial value is. My uncertainty regards which side of the globe you're on, for the solar time.

mod(A,P) = A - (INT(A/P) * P)

@apcraig
Copy link
Contributor

apcraig commented Jul 24, 2019

With regard to "when to merge". If we merge this to the trunk before the release, we should do a qc test as well as regenerate the gx1 model output "science" results. If we do not merge this, I think we can keep our 6.0.0 gx1 science output for this release and the release will be a bit cleaner and quicker. I can go either way. If we do merge, who is going to do the qc and the gx1 run. We probably want them both done before we start the release process?

@rgrumbine
Copy link
Contributor

The mod function approach is ok with me. I will say that we cannot expect values to lie inside [-720:720] either. iirc, the Hycom grid has longitudes over 1000. But mod takes care of that.

Elizabeth mentions a point different from algorithmic that this seems a good place to ask about on style. Namely, computational efficiency of mod vs. if (or other logical) test. I think at this point, it would be worth a timing test to see what is more computationally efficient. Processor architecture and compiler optimization skills have changed a lot since I was last using assembly language (early 1990s).

In my C/C++ codes, I usually have lines like
while (lon < -180.) lon += 360.;
while (lon > 180.) lon -= 360.;

F90 equivalent is less compact, but what's a little typing between friends?
DO WHILE (lon < -180.)
lon = lon + 360.
ENDDO
(and correspondingly for > 180.)

Or, for array-level 1-liner:
WHERE (lon < -180.) lon = lon + 360.
More generally:
DO WHILE (MINVAL(lon) < -180.)
WHERE (lon < -180.) lon = lon + 360.
ENDDO

In ancient days, this would have been bog-slow, not least because array operations weren't in the language. More generally, because logical tests were slow. These days, with array-orientation, vastly expanded instruction sets, and goodies like speculative out of order execution, ... as I said, I'm not so sure.

@eclare108213 eclare108213 self-assigned this Aug 2, 2019
@eclare108213
Copy link
Contributor

@dabail10 please try this

lontmp = mod(TLON+360,360)-180

and let us know whether it works. Either way, let's do the QC testing and update the sample output using your new run, and leave the efficiency question for a separate issue so we can get this PR into the next release. Thanks!

@dabail10
Copy link
Contributor Author

I looked at this approach and the problem is that if the longitudes are already between -180 and 180, this shifts things incorrectly. I believe the if conditions are the only robust way to do this.

@apcraig
Copy link
Contributor

apcraig commented Aug 22, 2019

On June 24, I pointed out a few other ways to deal with normalizing the lon values without ifs. What Dave has implemented should work, but there are lots of other ways to deal with it. I don't think what's there now is going to be a performance issue. I propose we defer that discussion at this point.

I think the main thing at this point is a qc test. We need that to merge this PR. @dabail10, I think you said you would take on that task. If you need some help, let me know.

@dabail10
Copy link
Contributor Author

Here are the results of the QC test:

Running QC test on the following directories:
/glade/scratch/dbailey/CICE_RUNS/qc_base
/glade/scratch/dbailey/CICE_RUNS/qc_test
Number of files: 1825
2 Stage Test Passed
Quadratic Skill Test Passed for Northern Hemisphere
Quadratic Skill Test Passed for Southern Hemisphere

Quality Control Test PASSED

@eclare108213
Copy link
Contributor

This result surprises me! Are you sure your test experiment is different from the base?

@dabail10
Copy link
Contributor Author

I did compare the ice thickness fields and there are definitely differences. However, according to the python script, the test passed. I though the python script automagically produced images as well? I think I am missing something.

@eclare108213
Copy link
Contributor

@dabail10
I'm ready to merge this PR, but I'd like some plots of the differences to be included in this PR conversation, for the record. There are a couple of options -- post your own plots, or use @mattdturner 's new python scripts. Please, would you do one or the other?
Thanks!

@dabail10
Copy link
Contributor Author

dabail10 commented Aug 30, 2019

Got the script working. Turns out I needed -Y for ssh on cheyenne. Very weird. Here are the plots from by base and test cases. There are some big differences, but I guess not climate changing overall according to the QC tests.

ice_thickness_qc_base_minus_qc_te
ice_thickness_qc_base
ice_thickness_qc_te

@eclare108213
Copy link
Contributor

Excellent. My rule of thumb is that any change less than roughly 10 cm is in the noise. The QC bears that out for this test. It'd good to know that this error was not actually climate changing.

@eclare108213 eclare108213 merged commit 804b703 into CICE-Consortium:master Aug 30, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants