Snicar-ad port to Icepack from MPAS-SI

[eclare108213]

• Short (1 sentence) summary of your PR:
  This PR contains snicar-ad changes to dEdd for snow-covered surfaces, including 3 sets of bug fixes.
• Developer(s):
  Code ported from E3SM/MPAS-SI by @eclare108213 and @apcraig
• Suggest PR reviewers from list in the column to the right.
• Please copy the PR test results link or provide a summary of testing completed below.
  Full test suite with modified CICE run on cheyenne, https://github.com/CICE-Consortium/Test-Results/wiki/cice_by_hash_forks#8fbe45fab39f24fc6f136a57ce1475b7f20c64a6. All results are bit-for-bit with CICE Consortium version with same dEdd bug fixes except dynpicard which is expected.
• How much do the PR code changes differ from the unmodified code?
  • bit for bit in Icepack tests but not in CICE for reasons we understand (see below)
  • different at roundoff level
  • more substantial
• Does this PR create or have dependencies on CICE or any other models?
  • Yes, the interface with the driver has changed
  • No
• Does this PR add any new test cases?
  • Yes
  • No
• Is the documentation being updated? ("Documentation" includes information on the wiki or in the .rst files from doc/source/, which are used to create the online technical docs at https://readthedocs.org/projects/cice-consortium-cice/.)
  • Yes
  • No, does the documentation need to be updated at a later time?
    • Yes
    • No
• Please provide any additional information or relevant details below:

This (was) a draft PR. Still to be done:

• use Icepack rhoi parameter value rather than the hard-coded rhoi in the shortwave module
• understand non-BFB results in CICE
• include PR#13
• update documentation
• update snicar-ad reference in icepack_shortwave.F90 (in subroutine compute_dEdd_5bd)
• add new fswthru_ diagnostics to compute_dEdd_5bd
• clean up more of the column alignment issues for readability
• complete testing

The 5-band snicar-ad additions to dEdd shortwave are being ported from MPAS-SI to Icepack. The 'ad' designation for snicar is included in the user interface (icepack_in) for consistency with the E3SM land component, but this modifier is dropped inside Icepack since it is the only snicar option.

Bug fixes:
The cleanup in commit 64ab4a1 is massive but did reveal some previously undetected bugs. This commit also could cause answer changes due to some precision changes. These bugs were fixed in this PR:

• missing code at the end of some tr_aero code sections, lost when modal aerosols were implemented
• reciprocal error in modal aerosol code (need to multiply by rnilyr and rnslyr rather than dividing by them)
• rsnw interpolation weights were backwards for the nsky=1 direct beam calculation in the 5-band dEdd routine (this bug was in the E3SM/MPAS-SI code, not in the Consortium repo).

This commit will definitely cause answers to change for some cases relative to the original code, as in CICE-Consortium#400. There are going to be conflicts regardless, so I didn't try to pull the changes from that PR and instead implemented them here directly.

Significant changes:

• There is a new module, icepack_shortwave_data.F90, which consolidates data tables from icepack_shortwave.F90, ice_forcing_bgc.F90 (in CICE), and new data tables derived from the SNICAR snow model.
• The original 3-band data is still used for non-snow-covered surfaces, and therefore the original nspint (number of spectral intervals) has been replaced with nspint_3bd, and nspint_5bd has been added.
• A new initialization routine has been added to the Icepack interface: icepack_init_radiation.
• A new subroutine has been added for snicar-ad 5-band snow calculations, compute_dEdd_5bd. This might be combined with the original routine (now called compute_dEdd_3bd), but it will need to be refactored due to conflicts in the loop and if-then-else structures.
• Unnecessary subroutine arguments have been removed from the shortwave calling tree, using available module data instead. This change is not backwards compatible.
• [PR#13] Calculations of g, w0 and tau have been refactored and simplified to be more precise (denominators with +puny are not necessary). This change is not BFB in some CICE tests. See Update CICE to run with eclare108213/Icepack branch snicar apcraig/CICE#100.
• Some constants were explicitly made _dbl_kind (instead of integer) in the modal_aerosol section of the code, which might not be BFB when that case is (eventually) tested.
• rhoi was hard-coded and now uses module data instead. Results will not be BFB if the driver value is not 917 kg/m^3.
• Some parameters that moved from ice_forcing_bgc.F90 to icepack_shortwave_data.F90 were single precision, and are now double precision. This change is not BFB in some CICE tests. See Update CICE to run with eclare108213/Icepack branch snicar apcraig/CICE#100.
• dEdd-snicar tests have been added.

Other changes:

• Repeat ridging generates a diagnostic only when niter > 1.
• Cleaned up diagnostic formats.
• Removed the assumption that the final grid cell in Icepack is land (for testing with a single cell).
• Scripts have been adjusted to improve build/compile efficiency with the new, larger data look-up tables.

[eclare108213]

I am changing the snw_ssp_table defaults and also making icepack_warnings_getall public as requested by @akturner at https://acme-climate.atlassian.net/wiki/spaces/ICE/pages/3526492170/Simulation+results+and+comparisons.

@apcraig it looks like we have both ICE_SNICARHC and USE_SNICARHC available in the code and/or scripts. I'm guessing that ICE_ is desired in order to avoid potential conflicts with other components in a coupled system, and so I should change the USE_SNICARHC instances to ICE_SNICARHC. Is that right?

[apcraig]

Regarding ICE_SNICARHC and USE_SNICARHC. ICE_SNICARHC is the variable name in cice.settings. USE_SNICARHC is the CPP in the code. So we want to keep both as they are now. ICE_SNICARHC = true/false and that turns on -DUSE_SNICARHC in the Makefile. Hopefully that makes sense?

When you say you're changing the snw_ssp_table defaults, does that mean you are changing the default in Icepack to "snicar". That would cause some problems. That would mean ICE_SNICARHC should be default "true" and that means the big snicar table is compiled by default when it's only needed when we use "dEdd_snicar_ad". The reason we have USE_SNICARHC is to avoid building that big table by default. That increases compile time significantly (say from 1 minute to 2 minutes). It may not sound large, but when running a test suite or debugging, it does matter.

[eclare108213]

I was responding to your comment above, suggesting to document USE_SNICARHC and correct the snw_ssp_table default value:

Just noticed a couple things. We probably need to add documentation of "USE_SNICARHC" CPP in ug_case_settings.rst. Also, I think the default for snw_ssp_table is "test". These could be corrected in a PR later.

ICE_SNICARHC is documented in ug_case_settings (in Icepack, at least), and maybe you could correct/document the other parts so that we'll all know how this is supposed to work? :)

[eclare108213]

Actually, I changed snw_ssp_table in icepack_in from 'test' to 'unknown_snw_ssp_table', which should be fine.

[apcraig]

Sorry if I wasn't clear. In the user guide, we should add USE_SNICARHC to the CPP list (section 3.4.1), https://cice-consortium-icepack.readthedocs.io/en/latest/user_guide/ug_case_settings.html#table-of-c-preprocessor-cpp-macros. With regard to the default setting of snw_ssp_table, I think you have it defined as "snicar" in the user guide, it should be "test", see line 268 is ug_case_settings. Finally, it's possible that "unknown_snw_ssp_table" will cause an abort in the implementation, we should check that and fix if it does. Once that works, we should also check that we don't need to modify some of the other test settings. Let me know if I can help with anything.

[apcraig]

To make icepack_warnings_getall truly public, you also have to add

use icepack_warnings, only: icepack_warnings_getall

to icepack_intfc.F90. That's how outside codes will get access to it.

[eclare108213]

I haven't thoroughly tested the last 2 commits and might still need to work on the code some more, but I think this PR is getting close.

Commit 64ab4a1 is massive -- ignore the whitespace differences and it's still massive -- perhaps it would be better to back off on that for the Icepack merge, but this exercise did reveal some previously undetected bugs. This commit also could cause answer changes due to some precision changes, which I probably should have committed separately. The second commit (9ae618c) fixes the bugs that have been found so far in icepack_shortwave.F90:

• missing code at the end of some tr_aero code sections, lost when modal aerosols were implemented
• reciprocal error in modal aerosol code (need to multiply by rnilyr and rnslyr rather than dividing by them)
• rsnw interpolation weights were backwards for the nsky=1 direct beam calculation in the 5-band dEdd routine.

The (not-yet-fixed) bugs are flagged with 'BUG' in 64ab4a1. This commit will definitely cause answers to change for some cases, as in CICE-Consortium#400. (There are going to be conflicts regardless, so I didn't try to pull the changes from that PR and instead implemented them directly. Heads up for @apcraig.)

@njeffery, the rsnw weights bug might affect the answers in E3SM runs too.

[apcraig]

This looks great to me. I can test for bit-for-bit vs CICE-Consortium main. Are there bugs fixes in this version that need to be implemented in CICE-Consortium (i.e. the nsky=1 bug?). I can add those to CICE-Consortium#400.

[eclare108213]

The new bug only affects the 5-band dEdd code, so it’s not in consortium main. Some of the precision changes might be applicable but I don’t think there is a pressing need to fix those yet, unless you just want to. Thanks!

[njeffery]

@eclare108213 : could you point me to the weights file that had the bug? Thanks

[eclare108213]

@njeffery the bug is here, line 4889 of icepack_shortwave.F90 in hash 64ab4a1
and the following lines

[njeffery]

Thanks @eclare108213 . I'll set up some runs in E3SM to test the impact of this bug.

[apcraig]

Testing against our "dEddfix1" branches of CICE+Icepack in CICE-Consortium, the only difference seems to be in the alt04 again. I might try to get them bit-for-bit again. Any chance you can point out the new "precision changes". I'm also fine just leaving things as they are at this point, with alt04 producing different results on the Consortium and E3SM versions for now.

[eclare108213]

@apcraig I was afraid you might ask that. The easiest way to find them is to search the diffs for _dbl_kind. I also changed a 0.01 to p01, and there were a couple of places where 1. became c1 (although I think those were not actually precision changes). That should catch all of them, but you could make sure by searching the diff for all of the named constants at the top of the shortwave module (use icepack_constants only: p01, c1, etc). Sorry I didn't separate those into a unique commit - I was fixing things as I came across them.

[apcraig]

@eclare108213, no problem. I'll do some diffs and see if I can apply the changes to the Consortium version. Thanks for pointing out what to look for.

[apcraig]

My initial results were incorrect, too many versions of the code laying around. I believe the latest changes did not change the results for the CICE test suite. So I don't think we need to port any of the new precision changes to the Consortium version. But I may do it anyway and recheck.

[eclare108213]

@apcraig, please, can you link to your test results in the main PR template above?

[apcraig]

I updated the test results section in the PR template.

[eclare108213]

Thank you! I'll merge this now to keep things moving forward.