Updated Numerical Solution for Wave Attenuation in sea ice

[erinethomas]

Pull Request Summary

This PR provides code modifications to apply a new MISC namelist parameter called "ICNUMERICS". Setting ICNUMERICS to be "TRUE" allows the user to place the 'sea ice' source terms with the other source terms, thus reducing the errors introduced with the time splitting method for wave attenuation in sea ice resulting in more accurate wave action in sea ice.

Description

The default behavior for wave attenuation in sea ice uses 'time splitting' for the wave action equation. This PR and the use of the new namelist parameter places the 'sea ice' source terms with the other source terms, thus reducing the errors introduced with the time splitting method for wave attenuation in sea ice resulting in more accurate wave action in sea ice.

Currently, this code is intended to be used in tandem with the newly proposed Meylan et al 2021 attenuation scheme : IC4M10 (see PR: #1293 ). This combination (IC4M10 plus ICNUMERICS=TRUE) has been tested in coupled CESM and E3SM runs, yet we note, at this time, this new numerics calculation for wave attenuation has not yet been tested with any other IC switch settings or other IC4 methods (1-9), although, we expect similar improvements for all IC methods.

This PR is the original work of @cmbitz and Vince Cooper (University of Washington). Collaborators include @erinethomas, @dabail10, @NickSzapiro-NOAA

Issue(s) addressed

This PR fixes the following Issue: #738

Commit Message

Introduction of a new MISC namelist parameter (ICNUMERICS) that allows for more accurate numerics scheme for wave attenuation in sea ice. Coauthors: Cecilia Bitz, Vince Cooper, Erin Thomas, David Bailey, Nick Szapiro.

Check list

• Branch is up to date with the authoritative repository (NOAA-EMC) develop branch.
• Checked the checklist for a developer submitting to develop.
• If a version number update is required, checked the updating version number checklist.
• If a new feature was added, a regression test for testing the new feature is added.

Testing

• How were these changes tested?
  This code has been tested by @cmbitz (details to be posted below), by @erinethomas in E3SM configurations with wave-sea ice coupling, and by @dabail10 in CESM configurations with wave-sea ice coupling.
• Are the changes covered by regression tests? (If not, why? Do new tests need to be added?)
• Have the matrix regression tests been run (if yes, please note HPC and compiler)?
• Please indicate the expected changes in the regression test output, (Note the list of known non-identical tests.)
• Please provide the summary output of matrix.comp (matrix.Diff.txt, matrixCompFull.txt and matrixCompSummary.txt):

[NickSzapiro-NOAA]

Maybe it is possible to progress with some discussion/review during this pause, @MatthewMasarik-NOAA ? I don't know if any IC4 developers in particular have feedback

[MatthewMasarik-NOAA]

Hi @NickSzapiro-NOAA, others should feel free to discuss. If this is not needed for either gfsv17 or gefsv13 though, than Avichal has reiterated that the waves team needs to pause any reviewing since we don't have the manpower currently.

[MatthewMasarik-NOAA]

@erinethomas, we will also process this in ~November. Thank you for the very informative PR header.

[cmbitz]

I've made a github repository about this issue at https://github.com/cmbitz/WW3Numerics where I've written up some more detailed notes about this issue and show some results in slides from a number of integrations to test it. I also provide a Python Notebook with some idealized examples that can be solved exactly where I compare the various solution options. In short, the modifications I'm suggesting can be implemented so they only matter in the sea ice. They can work for any IC/IS option combination, though I confess I've only tested with IC4Method10. I'm 100% sure it would improve wave amplitudes in the sea ice for any IC option. Happy to talk to anyone who wants to listen. I understand that this issue is on the back-burner, which is okay with me. I'll wait.

[erinethomas]

thanks for the extra details @cmbitz!

[sbanihash]

In a case of IC0 for example, where there is no IC1-IC5 defined, then VDIC has not been defined (look at line 744 in this same script) and the build will give this error: This name does not have a type, and must have an explicit type. [VDIC]
ATT=EXP(ICE*VDIC(IS)*DTG)
suggested changes is:
#ifdef WW3_IC4
#ifndef WW3_IC4_ACCURATE_NUMERICS
ATT = EXP(ICE * VDIC(IS) * DTG)
#endif
#endif

[sbanihash]

@erinethomas @NickSzapiro-NOAA do you agree with this change? If so could you please go ahead and make the changes?

[NickSzapiro-NOAA]

Keeping both makes sense to me to protect VDIC

[ErickRogers]

Why is this implemented as a switch? Could a logical variable in the IC4 namelist do the same thing?

[dabail10]

Why is this implemented as a switch? Could a logical variable in the IC4 namelist do the same thing?

Hi Erick. We were following the example of the W3_IC? options here. Happy to change it to a namelist option if this is preferred.

[ErickRogers]

Hi David.
Yes, I'm not sure if this is covered in the guidelines for developers, but the switches tend to be used for things that are pretty major, like selecting a propagation scheme, selecting a physics parameterization, selecting parallel computing method, or compiling in extra software (e.g., SCRIP). When a switch is changed, it requires a re-compile, whereas changing a namelist variable in the ww3_grid.inp file does not require a re-compile. Any new switch needs to be accommodated in the build system (/bin/build_utils.sh), which is an argument for "fewer switches=better". If this is limited to IC4, you can use a new logical variable in the IC4 namelist. If your option adds to the computation time, e.g., by forcing a small time step, then the original setting should be the default. (This will also make the discussion re: NCEP's "matrix" testing straightforward.) If this is intended to be used with any of the other IC# methods, you can put the new logical variable in the MISC namelist instead of the IC4 namelist.

[dabail10]

Ok. I will work with @erinethomas on this.

[ErickRogers]

Adding a namelist variable is very easy. I just did this the other day, to put in a new option for hacking the ice thickness. You can look at an existing variable to see how it is done, e.g. "grep IC4 *.F90 | grep -i ki". Any new variable will have 2 names. One is used in ww3_grid.inp and the other is used internally, e.g. IC4KI and IC4_KI. There are 3-4 files that need minor editing to add the new variable, as you will see with the output from above grep command. Keep in mind that after making this change, any old mod-def files that you have will no longer work.

[erinethomas]

This seems like a good way forward (to go with namelist options rather than switches) - I will update the PR as soon as I can!

[sbanihash]

@erinethomas @dabail10 Do we have any updates to this PR?

[erinethomas]

yes -I have worked out how to move this from a switch setting to a namelist parameter - I am still doing some final tests - code coming soon.

[erinethomas]

@cmbitz @dabail10 @NickSzapiro-NOAA
UPDATE: I have modified the code to use a logical Namelist setting called "ICNUMERICS"/"IC_NUMERICS". This namelist setting should be placed under "MISC" namelist category. I placed it under MISC (as opposed to IC4) since, in the future this numerical scheme should be applicable (and tested) for any Sea Ice dissipation scheme (not just IC4). This code should be applicable for users to use in any IC/IS combination. Keeping in mind this flexibility, would you all take a look at the modifications, and check if the logic I implemented is correct ... especially with respect to the changes I made in w3srcemd.F90 and any "IC/IS" switches used.... (for example Line ~1330, Line ~1416, Line 1460, etc.)

[sbanihash]

Thank you @erinethomas for your updates! Please let me know when you are ready for the second round of review.

[ErickRogers]

The recent updates by Erin look good to me. I don't have any other concerns or suggestions. I suppose it needs to get through the matrix testing now.

[sbanihash]

Thank you @ErickRogers for your review. @erinethomas I will be starting my second round of review soon. Please let me know if there are still changes or checks you need to do on you side before I do so. Thanks

[dabail10]

@cmbitz was actually advocating that we should remove this timestep limitation under the ice. However, I don't know if that is only when OFFSET=0.5 in w3srcemod.F90

[dabail10]

These changes work in the CESM code.

[sbanihash]

@erinethomas are you done with the changes that you were planning on making? Is the PR ready for testing again?

[erinethomas]

@sbanihash - I believe this PR is ready for review again.

[sbanihash]

yes -I have worked out how to move this from a switch setting to a namelist parameter - I am still doing some final tests - code coming soon.

@erinethomas could you edit your PR description to reflect this change?

[sbanihash]

@erinethomas just double checking to see if you have that one last change ready to be committed for me to submit my final review? Thank you in advance!

[erinethomas]

@sbanihash - thanks for the reminder - change to line 704 (removed the pointer) has now been made. Review ready.

[sbanihash]

I have reviewed the changes in this PR and tested the regression tests successfully. This PR will write an additional term in all ww3_grid.out files for ICNUMERICS and result in ww3_grid.out files being different for all regression tests. For that reason the regtests were first ran without any changes in the original PR and the tests were successful. In an attempt to compare matrixes, the write statements in ww3gridmd.F90 were taken out before the second round of regression tests and the matrix diff result is as below:

---

********************* non-identical cases ****************************

---

mww3_test_03/./work_PR2_UNO_MPI_e (1 files differ)
mww3_test_03/./work_PR2_UQ_MPI_d2 (16 files differ)
mww3_test_03/./work_PR3_UQ_MPI_e_c (1 files differ)
mww3_test_03/./work_PR2_UNO_MPI_d2 (16 files differ)
mww3_test_03/./work_PR3_UQ_MPI_d2_c (15 files differ)
mww3_test_03/./work_PR2_UQ_MPI_e (1 files differ)
mww3_test_03/./work_PR3_UNO_MPI_e_c (1 files differ)
mww3_test_03/./work_PR1_MPI_d2 (11 files differ)
mww3_test_03/./work_PR3_UNO_MPI_d2_c (15 files differ)
mww3_test_03/./work_PR3_UNO_MPI_d2 (16 files differ)
mww3_test_03/./work_PR3_UQ_MPI_d2 (15 files differ)
mww3_test_09/./work_MPI_ASCII (0 files differ)
ww3_tp2.10/./work_MPI_OMPH (7 files differ)
ww3_tp2.16/./work_MPI_OMPH (4 files differ)
ww3_tp2.6/./work_ST4_ASCII (0 files differ)
ww3_ufs1.3/./work_a (2 files differ)

---

All diff's are as expected. I approve this PR.
matrixCompSummary.txt
matrixCompFull.txt

matrixDiff.txt

[sbanihash]

@erinethomas I had already made that change and ran the tests, so when I double checked your latest commit and made sure it is the same as the change I had made in my runs, I went ahead and approved this PR. Thank you again and also thanks to @dabail10, @cmbitz and @NickSzapiro-NOAA for your work on this PR and to @ErickRogers for his review.
@erinethomas I will be creating a pull request next week with the ICNUMERICS truned on in IC4M10 so that this new feature is tested in our regtests as well. Stay tuned!

[sbanihash]

@erinethomas @cmbitz I have been working on adding a regtest for this feature, by changing ICNUMERICS=T in the input for the IC4M10 case in tic1.1 , but when comparing the results to a case where ICNUMERICS=F, I do not see any changes in output. I know you had mentioned that your tests have been successful, could you confirm that what is already in the IC4_M10 case is sufficient to test this numerical solution option?
Thanks

[dabail10]

These changes should only impact ice covered areas. Do you have coupling with the sea ice model?

[sbanihash]

Thanks @dabail10 for your prompt response. So you are confirming that we would need to test it in a coupled setting to see the results? If so I would have to look at setting up a test within the ufs-weather-model setup.

[dabail10]

You should definitely test it in the coupled setup. Or do you have a method in WW3 to read in sea ice data?

[sbanihash]

@dabail10 yes, we read an ice input for example this is how the ice is being read in the regtest that I am working on:

WW3/regtests/ww3_tic1.1/input_IC4_M10/ww3_prep_icecon.inp

Line 34 in 0a9e10f

2345 '../input_IC2_nondisp/icecon.156x3.txt'

[NickSzapiro-NOAA]

@sbanihash For coupled test, I have stale PR (ufs-community/ufs-weather-model#2072) once this is synced to dev/ufs-weather-model.

Idea is to switch on wave-ice coupling in cpld_control_pdlib_p8 and cpld_debug_pdlib_p8 RTs

[sbanihash]

Thanks @NickSzapiro-NOAA was just looking at your PR. Let me see if I can do a test within your branch with the addition of the changes in this PR. I was under the impression that as long as I am reading an ice input, turing on ICNUMERICS=T should result in answer changes even in a WW3 regtest, am I wrong?

[dabail10]

I agree. If you are reading in a sea ice field, then the numerics should be doing something in the presence of ice.

[NickSzapiro-NOAA]

Do you have plots for your ww3 regtest? Are there waves in the ice?

[sbanihash]

@NickSzapiro-NOAA I can definitely generate some and share.

[sbanihash]

@NickSzapiro-NOAA I have opened an issue of what I think to be the problem with this PR not working as expected when turning IC_NUMERICS on + a PR with the fixes which worked for the test case I have added. Still need to wait for all the regtests to finish to have the full matrix comparisons so that I could officially have the PR reviewed and tested by others.