Option to turn off internal stress in free drift conditions for B-grid EVP

[eclare108213]

• Short (1 sentence) summary of your PR:
  Adds a namelist option to set the ice internal stress term to zero when the ice area is less than 1.e-3, and applies this option for B-grid EVP in addition to the previous C- and CD-grid implementation.
• Developer(s):
  @eclare108213, building on C- and CD- grid modifications by @JFLemieux73 in Set rheology to zero for the C-grid when concentration is very small #1078 and Set rheology to zero for the C-grid when concentration is very small #1105
• 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.
  Quick suite results:

./cice.setup -m conda -e macos -s JRA55do --suite quick_suite --testid quick_rheofactU --bcmp quick_baseline_20260915

20 measured results of 20 total results
18 of 20 tests PASSED
0 of 20 tests PENDING
0 of 20 tests MISSING data
2 of 20 tests FAILED

The two failing tests are on the gbox128 test, for which the JRA55do forcing data's grid does not match. (I use -s JRA55do because I don't have the regular JRA55 data.)
Still need to run a QC test for free_drift_approx = .true., for completeness.

• How much do the PR code changes differ from the unmodified code?
  • bit for bit when free_drift_appox = .false.
  • different at roundoff level when free_drift_approx = .true., but the solutions are very similar.
  • more substantial
• Does this PR create or have dependencies on Icepack or any other models?
  • Yes
  • No
• Does this PR update the Icepack submodule? If so, the Icepack submodule must point to a hash on Icepack's main branch.
  • Yes
  • No
• Does this PR add any new test cases?
  • Yes
  • No - the new namelist flag is invoked as part of the C- and CD-grid tests
• 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/. A test build of the technical docs will be performed as part of the PR testing.)
  • Yes
  • No, does the documentation need to be updated at a later time?
    • Yes
    • No
• Please document the changes in detail, including why the changes are made. This will become part of the PR commit log.

EVP is stable in its default B-grid discretization, but slight changes to the discretization (such as using ice area at a different time level during the time step) can make it unstable. Because of this, the C- and CD-grid discretizations already set the internal stress term to zero for ice areas smaller than rheo_area_min = 1.e-3 by multiplying the internal stress by the factor rheofactE or rheofactN, equal to 0 (or 1, when the ice area is larger). The PR adds the analogous factor rheofactU to the B-grid discretization, and adds a namelist flag free_drift_approx to control it. The new namelist is also used for the C and CD grid cases.

To prevent a lot of "if present()" tests within the dynamics subcycling, rheofactU is not optional in the stepu routine, which meant that it also had to be added to the EAP dynamics module. However it is always 1 in that case -- it is not fully implemented there, pending future testing.

[apcraig]

I ran a QC test comparing gx1 baseline with gx1 baseline + flag_free_drift_approx=.true.. QC test passed,

INFO:__main__:Running QC test on the following directories:
INFO:__main__:  /glade/derecho/scratch/tcraig/CICE_RUNS/derecho_intel_smoke_gx1_256x1_medium_qc.fdbase
INFO:__main__:  /glade/derecho/scratch/tcraig/CICE_RUNS/derecho_intel_smoke_gx1_256x1_medium_qc.fdtest
INFO:__main__:Number of files: 1825
INFO:__main__:2 Stage Test Passed
INFO:__main__:Quadratic Skill Test Passed for Northern Hemisphere
INFO:__main__:Quadratic Skill Test Passed for Southern Hemisphere
INFO:__main__:Creating map of the data (ice_thickness_derecho_intel_smoke_gx1_256x1_medium_qc.fdbase.png)
INFO:__main__:Creating map of the data (ice_thickness_derecho_intel_smoke_gx1_256x1_medium_qc.fdte.png)
INFO:__main__:Creating map of the data (ice_thickness_derecho_intel_smoke_gx1_256x1_medium_qc.fdbase_minus_derecho_intel_smoke_gx1_256x1_medium_qc.fdte.png)
INFO:__main__:
INFO:__main__:Quality Control Test PASSED

I also ran a test suite on derecho gnu and intel including B, C, and CD tests, everything looks good. https://github.com/CICE-Consortium/Test-Results/wiki/cice_by_hash_forks#2c1f9dbd9d96ca95066ee6a75016e7c918d2e958

Finally, I wonder if it makes sense to have the free drift approx on by default for all cases and remove the namelist. My sense is that it's both more stable and makes more sense physically.

[eclare108213]

Thanks for testing, @apcraig.

Finally, I wonder if it makes sense to have the free drift approx on by default for all cases and remove the namelist. My sense is that it's both more stable and makes more sense physically.

I'm certainly open to turning the free drift approx on by default for EVP on all grids (it's not yet implemented/tested for EAP or VP). But I think we'd want to keep the namelist flag available for the operational centers to be able to turn it off easily, since it's not BFB. Please weigh in @TillRasmussen, @NickSzapiro-NOAA, @dabail10, @daveh150, @anton-seaice.

[TillRasmussen]

This looks fine. I agree with Tony that free drift should be default as I see it being more realistic. We can remove the other option. Alternatively keep it for a year or two until we have had a chance to test it more thorough in the coupled systems.

[NickSzapiro-NOAA]

If free_drift_approx is more realistic and robust, it is nice for that to be the default.

If it is only available for one dynamics scheme though, I think it needs to be configurable if folks should be able to switch among them

[eclare108213]

@dabail10, @daveh150, @anton-seaice: any heartburn about making this non-bit-for-bit change to EVP?

[eclare108213]

The free drift approximation is now set as the default in CICE. It is turned off for EAP and VP via namelist. The quick_suite test yields the same results as before.

I decided to keep the namelist option available for now, in case modeling centers need to test it both ways. Happy to remove that option in the near future (or now) if everyone is happy with the default behavior.

This PR needs a formal review to merge, @apcraig @TillRasmussen (or anyone). Thanks for your input - e

[apcraig]

The free drift approximation is now set as the default in CICE. It is turned off for EAP and VP via namelist. The quick_suite test yields the same results as before.

I decided to keep the namelist option available for now, in case modeling centers need to test it both ways. Happy to remove that option in the near future (or now) if everyone is happy with the default behavior.

Sounds like a good plan. Thanks.

[apcraig]

Hate to ask, is there an equivalent feature for eap and would it be easy to implement it?

[eclare108213]

Hate to ask, is there an equivalent feature for eap and would it be easy to implement it?

There is an equivalent feature, and at first glance it appears easy to implement. EAP calls both dyn_prep2, where rheofactU is computed (EAP ignores this) and stepu, where it is used (EAP rheofactU=1). It would be easy to get rheofactU from dyn_prep2 and use it in stepu, as in EVP. What I don't know is whether the structure tensor would also need to be altered to be consistent (subroutine stepa).