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

Update 5-band dEdd to support test data #472

Merged
merged 4 commits into from
Oct 27, 2023

Conversation

apcraig
Copy link
Contributor

@apcraig apcraig commented Oct 20, 2023

PR checklist

Add an interpolation method in icepack_shortwave.F90 for rsnw

Refactor portions of icepack_shortwave to improve robustness

Add snicar and snicartest test cases

Update the 5-band snicar shortwave rsnw_snicar_tab interpolation. This changes answers but QC tests pass. The new implementation checks whether the rsnw_snicar_tab is an array of integer values with deltas of value 1 (and sets the rsnw_snicar_chk variable). If so, it uses a shortcut to find the rsnw_snicar_tab bounds for interpolation, otherwise it calls the shortwave_search routine. In testing, the shortcut did not change performance much, but that could change in the future.

Add an interpolation method in icepack_shortwave.F90 for rsnw

Refactor portions of icepack_shortwave to improve robustness

Add snicar and snicartest test cases
gs = asm_prm_ice_drc (ns,nr-1)*(c1-delr) &
+ asm_prm_ice_drc (ns,nr )* delr
else
if (trim(snw_ssp_table) == 'snicar') then
Copy link
Contributor

Choose a reason for hiding this comment

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

Adding some comments, just to make sure I understand. This block is specific to the snicar table, which has a minimum snow grain radius of 30 and values are +1 from there.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

correct.

if (trim(snw_ssp_table) == 'snicar') then
! NOTE: Assumes delta rsnw is 1 and rsnw are basically integers
nr = ceiling(rsnw(ksnow)) - nint(rsnw_snicar_min) + 1
if (ceiling(rsnw(ksnow)) - rsnw(ksnow) < 1.0e-3_dbl_kind) then
Copy link
Contributor

Choose a reason for hiding this comment

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

This block sets the ks, ws, gs values to those of the nearest integer only when rsnw is less than but within 1.e-3 of that integer value, not when rsnw is above it. Should this be if (abs( )) then?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

correct, don't know why the diff is just one side of rsnw(ksnow), don't know if it should be abs or not, don't know why 1.e-3 was chosen.

ks = ext_cff_mss_ice_drc(ns,nr)
ws = ss_alb_ice_drc (ns,nr)
gs = asm_prm_ice_drc (ns,nr)
else
Copy link
Contributor

Choose a reason for hiding this comment

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

this block linearly interpolates the values as long as rsnw is not within 1.e-3 of the next highest integer. We could remove the denominator that is always = 1.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

correct.

+ asm_prm_ice_drc (ns,nr )* delr
endif
else
! linear interpolation
Copy link
Contributor

Choose a reason for hiding this comment

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

Finally, this block does linear interpolation all the time, for data not read from the snicar table. It seems that we could just use this all the time and eliminate the if-else blocks above. What would be the downside of that? (other than extra testing because the answers change)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

My goal here was to leave the "snicar" snw table results bit-for-bit and add support for the "test" table. That meant leaving all the logic and math as implemented. Changing the logic or the math could/would change answers. I agree that we could (should!) clean this up, but I propose it be done as a separate PR which would allow us to more carefully test the snicar results. This PR, to some degree, was just to add support for the "test" table.

@eclare108213
Copy link
Contributor

@proteanplanet I invite you to take a look at this PR. @apcraig is setting up an Icepack test for the 5-band shortwave and needed to generalize the hardwired 5-band coding to do that. Personally, I think we should fix the hard-wired weirdnesses in here (see the code diff) to just do the linear interpolation all the time, but that will change answers. @apcraig has modified the code in a way that allows our testing without changing the 5-band shortwave answers. We can fix this later, but I'm pinging you here in case you'd like it done sooner. I'm not very happy with the current 5-band coding.

This changes answers but QC tests pass.  The new implementation checks
whether the rsnw_snicar_tab is an array of integer values with deltas
of value 1 (and sets the rsnw_snicar_chk variable).  If so, it uses a
shortcut to find the rsnw_snicar_tab bounds for interpolation, otherwise
it calls the shortwave_search routine.  In testing, the shortcut did not
change performance much, but that could change in the future.
@apcraig
Copy link
Contributor Author

apcraig commented Oct 26, 2023

I just updated the implementation to cleanup the standard snicar 5-band rsnw_snicar_tab interpolation. This changes answers, but QC tests pass and the differences cannot be seen on the QC diff plot (attached). I ran run a small test suite focused on the snicar changes and that all works as expected. I am running a full test suite now, but don't expect any problems.

ice_thickness_cheyenne_intel_smoke_gx1_72x1_medium_qc_snicar qst1

ice_thickness_cheyenne_intel_smoke_gx1_72x1_medium_qc_snicar qsb1_minus_cheyenne_intel_smoke_gx1_72x1_medium_qc_snicar qst1

@apcraig
Copy link
Contributor Author

apcraig commented Oct 26, 2023

Full test suite on latest changes passes.

endif
enddo
do n = 2,size(rsnw_snicar_tab)
if (abs(rsnw_snicar_tab(n) - rsnw_snicar_tab(n-1) - c1) > 1.0e-6) then
Copy link
Contributor

Choose a reason for hiding this comment

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

why do you need 2 loops, one with n=1,size() using real(nint()) and the other with n=2,size() without real(nint()) but with -1?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

These loops are checking too different conditions. The first is that the array is basically all integers (to roundoff-ish). The second is that the delta between two elements in the array is 1.0 (to roundoff-ish). Both of those need to be true for the "shortcut" to work.

Copy link
Contributor

Choose a reason for hiding this comment

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

ok thx

else
! check whether rsnw_snicar_chk are integers with delta of one, makes search faster
if (rsnw_snicar_chk == 'unknown') then
rsnw_snicar_chk = 'delta1'
Copy link
Contributor

Choose a reason for hiding this comment

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

why not initialize rsnw_snicar_chk = 'delta1' above? Or assume it's something other than delta1 and rearrange the logic?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

You only want to do this check once. So it's initialized to "unknown". Then on first entry, it decides if it's "delta1" or "mixed" and then these checks are never done again. This is like a logical (first_call), but we need more than two states. We could have a first_call logical and then we could have a string that defaults to "mixed" then resets on first_call to "delta1" if it meets the criteria. Would that be clearer?

Copy link
Contributor

Choose a reason for hiding this comment

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

that would be clearer. Even better, could this check be moved into an initialization routine somewhere in columnphysics?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

OK, let me see what I can do.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

OK, I refactored this part. I added a new variable, rsnw_datatype, to icepack_shortwave_data.F90, so that variable is set at initialization. Then the value of that variable can be used to manage the rsnw interpolation logic.

Copy link
Contributor

Choose a reason for hiding this comment

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

This is much clearer now, thank you.

Set a character string, rsnw_datatype, at initialization.
@apcraig apcraig merged commit 6ad0f44 into CICE-Consortium:main Oct 27, 2023
1 check passed
Copy link
Contributor

@dabail10 dabail10 left a comment

Choose a reason for hiding this comment

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

This all makes sense to me and seems fine given that it is bfb.

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.

3 participants