Additional cleanup of Stieglitz snow code

[gmao-rreichle]

Miscellaneous, 0-diff cleanup of Stieglitz snow code:

• Print WARNING when encountering positive snow temperatures (PosSnowHeat) [see item 4 below].
• Pass tile lat/lon into StieglitzSnow_snowrt().
• Clean up logic of if-block for (tzero, ice1) conditions.
• Added comment about value of CPW (spec heat of ice) [see item 2 below].

This PR started out as an attempt to address occasionally large garbage fluxes during snow compaction, which would be a non-0-diff change [see item 1 below]. This particular change should be included in a future, non-0-diff PR.

Additional non-0-diff changes were attempted but found to be problematic or undesirable.

Below is the summary of all changes that were implemented or considered and ultimately.

---

Non-0-diff fixes of StieglitzSnow():

1. Maintain SWE after snow compaction (avoid large garbage fluxes)

For LAND tiles only: When snow compaction reaches the max allowed snow density, keep the snow mass (or snow water equivalent; SWE) unchanged and adjust the snow depth to limit the density to the max allowed (as proposed by @rdkoster).
Previously, after compaction reached the max allowed density, snow mass was removed to cut the density back to the max allowed value. This excess snow mass was added into a garbage flux term ("excs"). This mechanism resulted in large excs residual terms because of repeated densification beyond the max density and subsequent removal of snow mass excs. (This mechanism is retained for LANDICE tiles.)

A global, land-only simulation on the 36-km EASEv2 grid with hourly output for one winter season (~Oct-May) confirmed that the standard snow-related output (snow mass, snow depth, snow surface temperature, sublimation) all remain reasonable with the currently implemented bug fix.

The downside of the bug fix is that the snow mass and snow depth can grow indefinitely for land tiles that are "faux landice" (where snow does not melt off completely during the warm season). It is not clear if among the landice tiles there are "faux land" tiles (where snow/ice melt off completely during the warm season) and what, if anything, could be done to minimize the total number of "faux" tiles.

For reference, after snow densification and the aforementioned check on the max density, there is another check on the max snow depth here:

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/StieglitzSnow.F90

Line 890 in 3d764c6

if(sndzsum > maxsndepth) then

For land tiles, however, this snow depth check is essentially disabled (CATCH_MAXSNDEPTH=1.e20 [m]).
For landice tiles, the threshold is MAXSNDZ=15. [m].

Summary: Non-0-diff change implemented was originally implemented as discussed above. Rolled back in this commit to avoid non-0-diff changes. This particular change should be considered in a future non-0-diff change.

2. Custom value of CPW (spec heat of ice)

StieglitzSnow_CPW=2065.22 J/kg/K is hard-coded here:

https://github.com/GEOS-ESM/GEOSgcm_GridComp/blob/65465f6dbe9a524876225a10fe9b21b8294be362/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/StieglitzSnow.F90#L62C31-L62C43

Brief research suggests that this value is good near 0 deg C whereas MAPL_CAPICE=2000. J/kg/K is more appropriate near -10 deg C.

Summary: Added comment to source code, but did not change numerical values.

3. Make calc_tpsnow calculations consistent for scalar and vector inputs.

The legacy code uses hardwired 1.0000001 and 0.99999999 factors in the if condition when diagnosing snow temperature in the diagnostic routine for scalars [originally get_tf0d()]. These factors are not used in the diagnostic routine for vectors [originally get_tf_nd()].

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/StieglitzSnow.F90

Line 1363 in 65465f6

if (use_threshold_fac) then

Legacy threshold factors were disabled in this commit

Summary: Non-0-diff change was originally implemented as discussed above. Rolled back in this commit because the change resulted in bad snow simulations.

4. Fix bug in assigned snow temperature of partially frozen pack at 0 deg

The calc_tpsnow subroutine could produce a positive snow temperature if the input snow heat content is positive.

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/StieglitzSnow.F90

Line 1409 in 65465f6

t = hbw*tfac ! t>0. ??????

Summary: Check for positive snow heat content in StieglitzSnow_calc_tpsnow() added in this set of commits. Note that one call to StieglitzSnow_calc_tpsnow() in StieglitzSnow_snowrt() uses a tentative snow heat content and may produce positive snow temperatures. The check is turned off for this call.

5. Snow heat content adjustments after relayer()

After snow redistribution in relayer(), the thermal conditions of each new layer i are checked against those of the corresponding old layer i. Partially frozen conditions at 0 deg are not allowed to turn into fully frozen conditions below 0 deg (and vice versa), and the snow heat content and temperature are adjusted accordingly. It is unclear, however, if this makes sense when the redistribution is such that the old and new layers have minimal overlap. This is a particular concern when relayer is used in Landice (N_snow=15) or after adding snow mass increments from a land analysis.

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/StieglitzSnow.F90

Line 1279 in 65465f6

do i=1,N_snow

Summary: After initially abandoning the (ice10,tzero) check and snow heat content adjustment in this commit, the check was found to be important for maintaining a correct energy balance and reinstated in this commit (plus the subsequent 2 fixes).

6. Snow density check in relayer()?

After (or near the end of) relayer(), do we need a check of the snow density against the min/max allowed?

Summary: Snow density will be checked again at next time step. Leave source code unchanged.

7. Minimum value for snow cover area fraction

In a couple of places, the snow cover area fraction is set to a minimum of small=1.e6. Is there an inconsistency with the minimum of 0. implied elsewhere?

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/StieglitzSnow.F90

Line 428 in 65465f6

sndz = wesn/(max(areasc0,small)*rhofs)

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/StieglitzSnow.F90

Line 894 in 65465f6

areasc0 = max(small, areasc0 )

Summary: Slight inconsistency is acceptable. Leave source code unchange.

cc: @rdkoster @lcandre2

[gmao-rreichle]

@lcandre2, @zhaobin74: I'm wondering why the "tiletype" isn't passed into StieglitzSnow_snowrt() by using the simple constant MAPL_LANDICE, as in call SNOWRT( 1, NUM_SNOW_LAYERS, MAPL_LANDICE, ...

Instead, GEOS_LandIceGridComp.F90 uses TileTypes:

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSlandice_GridComp/GEOS_LandIceGridComp.F90

Line 3089 in 9e8ad9a

call SNOWRT(1,NUM_SNOW_LAYERS,TILETYPES(k), &

which is obtained from MAPL here:

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSlandice_GridComp/GEOS_LandIceGridComp.F90

Line 2183 in 9e8ad9a

TILETYPES = TILETYPES, &

Is it possible that TileTypes is a vector for all tiles (incl. land, lake, landice, and ocean)? I assume TileTypes is somehow zoomed into just the landice tiles, but then it would be equivalent to using the scalar constant MAPL_LANDICE as suggested above, which is much simpler.

If TileTypes is not zoomed into just the landice tiles, could StieglitzSnow_snowrt() be getting the wrong tile type? Again, I can't imagine this is the case, but I thought I should flag it just in case. Sorry I'm so ignorant about MAPL's inner workings.

[zhaobin74]

@lcandre2, @zhaobin74: I'm wondering why the "tiletype" isn't passed into StieglitzSnow_snowrt() by using the simple constant MAPL_LANDICE, as in call SNOWRT( 1, NUM_SNOW_LAYERS, MAPL_LANDICE, ...

Instead, GEOS_LandIceGridComp.F90 uses TileTypes:

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSlandice_GridComp/GEOS_LandIceGridComp.F90

Line 3089 in 9e8ad9a

call SNOWRT(1,NUM_SNOW_LAYERS,TILETYPES(k), &

which is obtained from MAPL here:

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSlandice_GridComp/GEOS_LandIceGridComp.F90

Line 2183 in 9e8ad9a

TILETYPES = TILETYPES, &

Is it possible that TileTypes is a vector for all tiles (incl. land, lake, landice, and ocean)? I assume TileTypes is somehow zoomed into just the landice tiles, but then it would be equivalent to using the scalar constant MAPL_LANDICE as suggested above, which is much simpler.

If TileTypes is not zoomed into just the landice tiles, could StieglitzSnow_snowrt() be getting the wrong tile type? Again, I can't imagine this is the case, but I thought I should flag it just in case. Sorry I'm so ignorant about MAPL's inner workings.

@gmao-rreichle, you are right about TileTypes being zoomed into just the landice tiles. Yes, it is equivalent to using MAPL_LANDICE. I don't recall the details as to why the call to snow_rt used TileTypes. I think using MAPL_LANDICE in the call could simply it.

[gmao-rreichle]

@sdrabenh : This PR is finally ready for you. Thanks for your patience!

@biljanaorescanin ran the GCM tests. As expected for a non-0-diff change, the tests runs finished ok but failed the comparisons against the Baseline results.

@lcandre2 approved the PR for the landice team, but after her approval I made a very minor and trivially 0-diff commit (removed a few comment lines), which may have formally invalidated her approval.

When you merge the PR, you may want to include the 5 bullet points of the PR "Summary" above in the commit message.

@weiyuan-jiang @rdkoster

[gmao-rreichle]

@biljanaorescanin: This PR is now ready for you and Scott. The CI tests all pass, and I also successfully ran the full suite of nightly LDAS tests. Please let me know if you have any questions, thanks, Rolf

[biljanaorescanin]

@sdrabenh PR was tested for 1day AMIP; REPLAY and Inc REPLAY and all checkpoint/restart files were zero diff between exp and baseline.
Regress job returned success status.

[sdrabenh]

@gmao-rreichle please approve this PR if it is ready

[gmao-rreichle]

@gmao-rreichle please approve this PR if it is ready

@sdrabenh : I initiated the PR and cannot formally approve it. @biljanaorescanin approved it for the land group. It is ready for merging.

[gmao-rreichle]

I guess we'll still need approval from @GEOS-ESM/landice-team. @lcandre2, @rcullath, can one of you please take a look and approve, or let me know if you have any questions or concerns? Thanks!

[sdrabenh]

@gmao-rreichle please approve this PR if it is ready

@sdrabenh : I initiated the PR and cannot formally approve it. @biljanaorescanin approved it for the land group. It is ready for merging.

@gmao-rreichle you are right of course. My oversight. I'll wait for @lcandre2 or @rcullath approval