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

Updates to advanced snow physics implementation #449

Merged
merged 5 commits into from
Aug 18, 2023

Conversation

eclare108213
Copy link
Contributor

@eclare108213 eclare108213 commented Jul 31, 2023

Draft PR while @njeffery tests the code in E3SM. There is still some code cleanup to do, too.

  • Short (1 sentence) summary of your PR:
    Updates to advanced snow physics implementation from E3SM.

  • Developer(s):
    @eclare108213 @njeffery

  • 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.

The answers change for the advanced snow physics tests due to changes in logic and a bug fix

Icepack base_suite:

20 measured results of 20 total results
16 of 20 tests PASSED
0 of 20 tests PENDING
0 of 20 tests MISSING data
4 of 20 tests FAILED
icepack/testsuite.test_snow_suite2> results.csh | grep FAIL
FAIL chicoma_intel_smoke_col_1x1_run1year_snw30percent compare base_snow_suite2 different-data
FAIL chicoma_intel_smoke_col_1x1_run1year_snwitdrdg compare base_snow_suite2 different-data
FAIL chicoma_intel_smoke_col_1x1_run1year_snwgrain compare base_snow_suite2 different-data
FAIL chicoma_intel_smoke_col_1x1_run1year_snwgrain_snwitdrdg compare base_snow_suite2 different-data

CICE base_suite:
The only changes to CICE other than the icepack submodule are the chicoma machine files.
The two jra55do runs fail because I haven't yet downloaded the data.

cice.snow/testsuite.test_snow2> results.csh | grep FAIL
FAIL chicoma_intel_smoke_gx3_8x2_jra55do run -1 -1 -1
FAIL chicoma_intel_smoke_gx3_8x2_jra55do test
FAIL chicoma_intel_smoke_gx3_4x1_debug_icdefault_snwgrain_snwitdrdg complog base_snow different-data
FAIL chicoma_intel_smoke_gx3_4x1_debug_icdefault_snwgrain_snwitdrdg compare base_snow 41.57 13.19 15.29 different-data
FAIL chicoma_intel_smoke_gx3_4x1_debug_icdefault_snw30percent compare base_snow 41.13 13.36 14.66 different-data
FAIL chicoma_intel_restart_gx3_8x2_icdefault_snwgrain_snwitdrdg complog base_snow different-data
FAIL chicoma_intel_restart_gx3_8x2_icdefault_snwgrain_snwitdrdg compare base_snow 12.29 2.38 4.55 different-data
FAIL - test failed
7 of 345 tests FAILED

CICE QC test:

INFO:main:Running QC test on the following directories:
INFO:main: /lustre/scratch5/.mdt1/eclare/CICE_RUNS/chicoma_intel_smoke_gx1_44x1_medium_qc.qc_base
INFO:main: /lustre/scratch5/.mdt1/eclare/CICE_RUNS/chicoma_intel_smoke_gx1_44x1_medium_qc.qc_test
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

  • How much do the PR code changes differ from the unmodified code?

    • bit for bit
    • different at roundoff level
    • more substantial
  • Does this PR create or have dependencies on CICE or any other models?

    • Yes
    • 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:

  • Add chicoma machine files.

  • Use snwredist and snwgrain instead of tr_snow in column physics (tr_snow is still used in Icepack's driver but is not required for other drivers).

  • Correct tr_snow logic to use snwredist and snwgrain, as needed. This change is not BFB due to tracer values when the various options are NOT in use.

  • Fix snwgrain bug: the tracers smicen and smliqn were incorrectly calculated, missing the factor of aice needed to translate between grid-cell-mean thickness and "actual" ice thickness for the thermodynamics.

@apcraig
Copy link
Contributor

apcraig commented Jul 31, 2023

This looks fine to me. Is there a reason it's a Draft at this point? Do you want me to do any additional testing (with the caveat that cheyenne is down most of this week, so it might be a few days)?

@eclare108213
Copy link
Contributor Author

I am satisfied with the testing that I've done so far, although I still plan to look more closely at the output differences. I compared base_suite and the QC -- are there other test suites that need to be run in this case?

I'd like to keep this as a draft for now, in case Nicole turns up something that needs to be changed.

Copy link
Contributor

Choose a reason for hiding this comment

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

@eclare108213 : For the snwredist options, it seems we're only supporting ITDrdg and none. Are bulk, 30percent and ITDsd defunct?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

We're supporting ITDrdg, bulk and none. The old 30percent is available - it's set using 'bulk' with a (new) namelist parameter that can be any percentage. ITDsd is defunct.

Copy link
Contributor

Choose a reason for hiding this comment

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

Good to know. I'll have to look into how the bulk option would work in mpas-seaice. Does 'bulk' use the effective snow density?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I think bulk works the same way as 30percent did before, just generalized the percentage. So it only rearranges the snow without changing the density. And yes the density tracer is still available in that case - it should remain constant, although I haven't looked at that yet after this latest round of changes.

Copy link
Contributor

Choose a reason for hiding this comment

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

I'm concerned that the way it's implemented, we are requiring the density tracer for the bulk option. Otherwise it will fail. And it doesn't do anything.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It doesn't fail in Icepack (or CICE). Can MPAS-SI not implement it that way? With it, analysis scripts looking for it (for comparison with ITDrdg) won't fail, and its values can be verified to be constant.
Have you tried running this icepack version in your E3SM tests? I think it's coded so that if the driver doesn't have the tracer, icepack will still handle it gracefully. But I'm not sure about that. If it doesn't work, then I can change it.

Copy link
Contributor

Choose a reason for hiding this comment

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

I haven't tested (recently) snwredist = bulk and snow effective density off, but I will. But tracers are passed to icepack in the tracer array which in mpas is only made up of active tracers. The index for nonactive tracers is set to 0. So in something like

do k = nt_rhos, nt_rhos+nslyr-1
trcrn(k,n) = max(trcrn(k,n)-rhosmin, c0)

I suspect F90 won't like that nt_rhos = 0, and the rest of the code will overwrite something valuable.

I could be wrong. I've almost got your changes implemented. I'll try running with different options.

Copy link
Contributor

Choose a reason for hiding this comment

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

The changes 141, 267 and 268 would be a problem with other ktherm options and the snow grain tracer active. Since the internal assignments of massice and massliq are wrapped in a ktherm == 2 statement, they would be sent as zeros to thickness_changes.

Copy link
Contributor

Choose a reason for hiding this comment

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

My first look didn't catch this. Probably the simplest fix is to just pull the assignment statements out of the ktherm == 2 if-else and move to after both the temperature_changes routines.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Do we need to support advanced snow physics on BL99 thermo? We could abort if that combination is attempted. Or we could initialize massice and massliq according to the BL99 assumptions (use rhos for massice and 0 for massliq?).

Copy link
Contributor

Choose a reason for hiding this comment

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

I think just moving it out of the if statement is easiest. It has it's own if snwgrain ... wrapper. I'm certainly not in favor of supporting a massive number of config combinations, but the snow stuff doesn't really care about mushy or BL99. That's not to say I've tested it with BL99...

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Do you mean line 333ff?

Copy link
Contributor

Choose a reason for hiding this comment

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

Yes

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Okay, I agree that should be changed.

…em in columnphysics/. Compute massice and massliq for either mushy or BL99 thermo.
Copy link
Contributor

Choose a reason for hiding this comment

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

@eclare108213 : line 1294 in icepack_tracers.F90 should also be changed to rhosnew.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

done

Copy link
Contributor

Choose a reason for hiding this comment

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

This else part of this statement in icepack_snow.F90 is a problem. rhos_cmpn can't be used at all if snow effective density is false otherwise it will overwrite ice area.
Lines 318 to 320
else
rhos_cmpn(:,:) = rhos
endif

Copy link
Contributor

Choose a reason for hiding this comment

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

And the first part of the if statement
if (trim(snwredist) /= 'none') then
should be changed to
if(snwredist(1:3) = 'ITD') 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.

Agreed. Will it work to just remove the else part? This is part of the reason for defining the all of the tracers when tr_snow=T.

Copy link
Contributor

Choose a reason for hiding this comment

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

yes, deleting L318 and 319 works

@apcraig
Copy link
Contributor

apcraig commented Aug 4, 2023

Full test suite on cheyenne looks fine. The snw tests changed answers, but otherwise, all results pass and are bit-for-bit. This testing was done up to #2e12190. If you want me to test the latest version as well, let me know. Test results are here, https://github.com/CICE-Consortium/Test-Results/wiki/cice_by_hash_forks#993d058f5042a29d7d8b645313d9352632affa88.

Copy link
Contributor

@njeffery njeffery Aug 7, 2023

Choose a reason for hiding this comment

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

@eclare108213 : This expression (L1245 and 1246) forces the ice mass tracer to maintain it's current ice density by adding/subtracting mass from the snow liquid tracer.

1240  if (hsn > puny) then    ! add snow with enthalpy zqsn(1)                                                            
1241        dhs = econ / (zqsn(1) - rhos*Lvap) ! econ < 0, dhs > 0  
1242         mass  = massice(1) + massliq(1)                                                                                  
1243         massi = c0                                                                                                       
1244         if (dzs(1) > puny) massi = c1 + dhs/dzs(1)
1245         massice(1) = massice(1) * massi                                                                                  
1246        massliq(1) = max(c0, mass + rhos*dhs - massice(1)) ! conserve new total mass 

It's only going to be round off level changes at this point because the ice mass is always near rho, but it's not what we want. It should just be

massice(1) = massice(1) + dhs * rhos

[edited formatting only]

Copy link
Contributor Author

Choose a reason for hiding this comment

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

To clarify, @njeffery, you are suggesting that we should only change massice(1) and not change massliq(1)? Having thought through this, I think you are right. It's pretty confusing, though, so I'm recording my thought process here for later reference:
Condensation could be done either way, with all of the condensation becoming solid ice or dividing it between ice and liquid. I think the massi code above divides it evenly based on what the tracers were to begin with. But the latent heat values would probably be different in that case. Is that the issue here? This calculation is using Lvaporization, which is between liquid and gas, but if we're condensing from gas directly to ice, shouldn't it be Lsublimation? Answer: No, because the enthalpy definition that we use in the model includes the heat needed to convert between solid and liquid phases (and to bring the liquid to 0 C). Here, we are only changing the ice and liquid mass tracers, which shouldn't affect the solution at this point in the code -- the liquid mass may be used later for melt ponds, which would change the solution in those configurations.

Copy link
Contributor

Choose a reason for hiding this comment

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

@eclare108213 : yes. Just massice. New snow formation has massice = rhos and this is most consistent with both our energy and mass exchanges. The points you're making above about vapor to liquid and vapor to ice are also worth thinking more about, especially if we ever start changing snow density in the thermodynamics. But we okay for conservation. You'll see differences through the melt ponds and the snow grain radius.

@dabail10
Copy link
Contributor

Are the bulk changes not in this PR?

@njeffery
Copy link
Contributor

@dabail10 : No.

@apcraig
Copy link
Contributor

apcraig commented Aug 18, 2023

Full test suite of Icepack and CICE run on cheyenne. Plus a QC test. Everything looks good.

Icepack is fully bit-for-bit except for snw test cases, https://github.com/CICE-Consortium/Test-Results/wiki/icepack_by_hash_forks#8fe89cdbd3ed62f14baade3170bacdfb6d11cc60

CICE is bit-for-bit except for snw test cases, https://github.com/CICE-Consortium/Test-Results/wiki/cice_by_hash_forks#993d058f5042a29d7d8b645313d9352632affa88. In addition, 3 results are still pending, but those should be fine. Also, there is a problem with some part of the test scripts on cheyenne, but manual comparison of a few of those cases show the results are bit-for-bit. Will need to look into the test script problem.

The QC test passes,

INFO:main:Running QC test on the following directories:
INFO:main: /glade/scratch/tcraig/CICE_RUNS/cheyenne_intel_smoke_gx1_72x1_icdefault_medium_qc_snwgrain_snwitdrdg.qcb
INFO:main: /glade/scratch/tcraig/CICE_RUNS/cheyenne_intel_smoke_gx1_72x1_icdefault_medium_qc_snwgrain_snwitdrdg.qc
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_cheyenne_intel_smoke_gx1_72x1_icdefault_medium_qc_snwgrain_snwitdrdg.qcb.png)
INFO:main:Creating map of the data (ice_thickness_cheyenne_intel_smoke_gx1_72x1_icdefault_medium_qc_snwgrain_snwitdrdg.qc.png)
INFO:main:Creating map of the data (ice_thickness_cheyenne_intel_smoke_gx1_72x1_icdefault_medium_qc_snwgrain_snwitdrdg.qcb_minus_cheyenne_intel_smoke_gx1_72x1_icdefault_medium_qc_snwgrain_snwitdrdg.qc.png)
INFO:main:
INFO:main:Quality Control Test PASSED

ice_thickness_cheyenne_intel_smoke_gx1_72x1_icdefault_medium_qc_snwgrain_snwitdrdg qcb_minus_cheyenne_intel_smoke_gx1_72x1_icdefault_medium_qc_snwgrain_snwitdrdg qc
ice_thickness_cheyenne_intel_smoke_gx1_72x1_icdefault_medium_qc_snwgrain_snwitdrdg qc
ice_thickness_cheyenne_intel_smoke_gx1_72x1_icdefault_medium_qc_snwgrain_snwitdrdg qcb

@apcraig
Copy link
Contributor

apcraig commented Aug 18, 2023

@eclare108213, looks like this could be merged. Then the CICE PR needs to be updated with this version of Icepack, and it can be merged. Are you able to do both quickly? I can do the CICE merge after a quick final review and after the Github actions runs. Maybe we can get both in before weekend testing. If not, that's OK too.

@eclare108213 eclare108213 marked this pull request as ready for review August 18, 2023 22:23
@eclare108213 eclare108213 merged commit 86cae16 into CICE-Consortium:main Aug 18, 2023
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.

4 participants