Fix Tsnice calculation #542

dabail10 · 2025-09-19T17:54:54Z

For detailed information about submitting Pull Requests (PRs) to the CICE-Consortium,
please refer to: https://github.com/CICE-Consortium/About-Us/wiki/Resource-Index#information-for-developers

PR checklist

NickSzapiro-NOAA · 2025-09-22T06:12:12Z

The range of sitempsnic is certainly better in 24-hour history from ufs-weather-model cpld_control_gfsv17 RT, especially no (huge) negatives. Left is current, right is this PR (with different colorbars):

NickSzapiro-NOAA · 2025-09-22T06:22:54Z

columnphysics/icepack_therm_vertical.F90

      if ((hslyr+hilyr) > puny) then
         if (hslyr > puny) then
-            Tsnice = (hslyr*zTsn(nslyr) + hilyr*zTin(1)) / (hslyr+hilyr)
+            Tsnice = Tsnice + aicen*((hslyr*zTsn(nslyr) + hilyr*zTin(1)) / (hslyr+hilyr))


Can I check that I understand? Does this aicen factor get "canceled" in ice_history in:

! Only average for timesteps when ice present if (avail_hist_fields(n)%avg_ice_present) then do j = jlo, jhi do i = ilo, ihi if (tmask(i,j,iblk)) then a2D(i,j,n,iblk) = & a2D(i,j,n,iblk)*avgct(ns)*ravgip(i,j) endif

No, this is aggregating the category Tsnice into the grid cell Tsnice. The time averaging in ice_history.F90 will come in a CICE PR!

Thanks! Then isn't a Tsnice = Tsnice/sum(aicen) factor missing? For example, for single category aicen almost puny

Sorry this is so confusing. I should really call this:

Tsnice = sum(aicen(n)*Tsnicen(n)) for n=1,5

this is done for every point in space within icepack. The history averaging in CICE is then:

sum(aice(t)*Tsnice(t)) / sum(aice(t)), for t over the averaging period. The CICE changes are coming in a future PR!

I'm also confused :)

Tsnice = sum(aicen(n)*Tsnicen(n)) for n=1,5

sum(aicen(n)) is not 1, so Tsnice ends up weighted by aice here right ?

... The history averaging in CICE is then:

sum(aice(t)*Tsnice(t)) / sum(aice(t)), for t over the averaging period. The CICE changes are coming in a future PR!

Doesn't this mean its weighted by aice twice ?

Oh - maybe this line would be correct now:

worka(i,j) = aice(i,j,iblk)*(Tsnice(i,j,iblk)/aice_init(i,j,iblk)+Tffresh)

https://github.com/CICE-Consortium/CICE/blob/b01032506d71d86ca51cef8c3c28c1ab87532833/cicecore/cicedyn/analysis/ice_history.F90#L2768

the divide by aice_init corrects it for the aice weighting applied during thermo vertical

No, this is aggregating the category Tsnice into the grid cell Tsnice.

This is true only if Tsnice=0 for the open water area. Assuming that, then the grid cell average Tavg = sum(aicen(n)*Tsnice(n))/sum(aicen(n)) =sum(aicen(n)*Tsnice(n))/1 where the sums are for n=0,ncat mathematically but only need to be taken for n=1,ncat in the code. Confusing, yes, but possibly okay here. (Explicit code comments to this effect would be helpful.) A better assumption for Tsnice for open water might be the ocean temperature. On the other hand, if the data request is for Tsnice only over the ice-covered area, then it had better be
sum(aicen(n)*Tsnice(n))/sum(aicen(n)) for n=1,ncat, where aicen might actually be aicen_init at this point in the code.

I'm generally suspicious of any 'mixed' ice area scaling using both aice and aice_init. Those will need to be considered on a case-by-case basis, in my opinion.

dabail10 · 2025-09-22T15:59:57Z

The range of sitempsnic is certainly better in 24-hour history from ufs-weather-model cpld_control_gfsv17 RT, especially no (huge) negatives. Left is current, right is this PR (with different colorbars):

Great. There is another fix coming in CICE to account for a case when aice_init is zero, but aice is not.

apcraig · 2025-09-23T19:53:29Z

Is this ready for review and merge? Are there any test results verifying bit-for-bit (except for the diagnostic)?

dabail10 · 2025-09-23T19:54:57Z

Is this ready for review and merge? Are there any test results verifying bit-for-bit (except for the diagnostic)?

I posted the test results in the comments above.

apcraig · 2025-09-23T20:04:10Z

Is this ready for review and merge? Are there any test results verifying bit-for-bit (except for the diagnostic)?

I posted the test results in the comments above.

where?

dabail10 · 2025-09-23T21:51:41Z

Is this ready for review and merge? Are there any test results verifying bit-for-bit (except for the diagnostic)?

I posted the test results in the comments above.

where?

I guess I forgot to save it earlier. Now in the first comment.

https://github.com/CICE-Consortium/Test-Results/wiki/icepack_by_hash_forks#e9d626f0e5b743e143a2e87248a1aa22ee4f3751

anton-seaice · 2025-09-24T00:04:18Z

columnphysics/icepack_therm_vertical.F90

         einter = einter + hilyr * zqin(k)
      enddo ! k

      Tsnice = c0


I think this Tsnice = c0 is in the wrong location.

thermo_vertical is called for once for each thickness category, which means Tsnice=0 gets reset every time and the result would be just the temperature for the ncat = 5

Good catch. Fixed.

Tsnice is now initialized in icedrv_flux.F90 and ice_flux.F90.

columnphysics/icepack_therm_vertical.F90

apcraig · 2025-09-24T17:31:34Z

As part of these updates, it makes sense to me to do some manual checks of the updated diagnostic output. Maybe run the model for a day, write out a bunch of data either in history or via stdout, and verify that the new diagnostics are generating results consistent with the raw data. If we do that, these errors should get caught, and we'll have more confidence in the results as well.

dabail10 · 2025-09-24T17:36:43Z

As part of these updates, it makes sense to me to do some manual checks of the updated diagnostic output. Maybe run the model for a day, write out a bunch of data either in history or via stdout, and verify that the new diagnostics are generating results consistent with the raw data. If we do that, these errors should get caught, and we'll have more confidence in the results as well.

Ok. This is part one of a set of fixes for Icepack and CICE. I am working on the CICE PR now. In the meantime, here is a plot from the CICE run.

apcraig

Needs further review

dabail10 · 2025-10-01T16:05:30Z

@anton-seaice @eclare108213 @NickSzapiro-NOAA Have I addressed all of your concerns here? I know Nick is not able to look at this right now.

NickSzapiro-NOAA · 2025-10-01T16:12:30Z

Thanks @dabail10 ! fyi, EMC contract and federal staff are still working. Great to hear! I have a friend, Todd Arbetter who is a contractor working on UFS now. More in the hurricane group.

eclare108213

I've looked at this code (the PRs in both Icepack and CICE) fairly closely and I think it is okay. I'll approve both.

Generally speaking, the history code is VERY confusing because the averaging calculations are spread between Icepack and CICE and also among CICE history subroutines. I'm not 100% convinced that it's outputting what we think it is for every field (not just those in this PR).

Some additional comment statements to explain what's happening when would be helpful. E.g. I think this comment

          ! Only average for timesteps when ice present

ought to be

          ! Only average when/where ice is present

because the multiplication by ravgip divides by both time and area, whereas avgct is only for time (and ravgip is not defined in the declaration comments!).

I'll suggest some unit tests in a new issue, which should impart a bit more confidence in general.

apcraig · 2025-10-28T18:51:02Z

Should we merge this? Or do we want to wait for the CICE PR to be ready too?

dabail10 · 2025-10-28T19:02:00Z

I have one more question for the others. To get the snow-ice interface temperature I am doing a thickness weighted average:

Tsnice = Tsnice + aicen*((hslyr*zTsn(nslyr) + hilyr*zTin(1)) / (hslyr+hilyr))

So, this assumes a straight line from the bottom snow level to the top ice level. However, I am values above 0C in my history variable. I tried using both aice_init and aice in the averaging, but this did not help. I don't think my weighted average could be generating values above 0C, but I wonder if I should use ksno and kice here? Or is it possible that the mushy physics is generating temperatures above 0C?

Dave

dabail10 · 2025-10-28T19:34:55Z

So, for example, say zTsn=-10C and zTin = -1C and hslyr=0.1 and hilyr=0.25. Then this formula is:

(0.1*(-3) + 0.25*(-1C)) / 0.35 = (-0.3-0.25) / 0.35 = -1.57C

Whereas a conductivity based thing would be something like this?

ksno * (zTsn - Tsnice) / hslyr = kice * (Tsnice - zTin) / hilyr

In case you were curious, with the numbers above this comes out around -3.5C.

dabail10 · 2025-10-28T19:55:33Z

Also, I did a search again in case there was already an interface temperature. There is iTin in the icepack_brine.F90 module, but I couldn't find one otherwise.

eclare108213 · 2025-10-29T20:09:29Z

@dabail10 I don't think your thickness-weighted temperature formula gives you the temperature at the interface. That is probably the average temperature over the two layers, which isn't the same thing as the interface temperature. Setting up a grid with temperature on the y axis and distance on the x axis (i.e. the snow-ice column lying sideways), I get the interface temperature as -7.42857C in your example above with -10C, 0.1 m snow and -1C, 0.25 m ice layers.

dabail10 · 2025-10-29T21:03:48Z

I agree. I had naively assumed it was a weighted average of the two layers. Instead maybe the straight line fit of:

T = (Tice - Tsno) / 0.5 *(hslyr + hilyr) * x + Tsno

So, when x = 0, T = Tsno. when x = 0.5*(hslyr+hilyr), T=Tice. With this in mind, x = 0.5*hslyr is the snow-ice interface and T = (Tice-Tsno)*hslyr / (hilyr+hslyr) + Tsno.

Using Tsno = -10 and Tice = -1, hslyr=0.1, hilyr=0.25, this gives -7.429 as you found.

eclare108213 · 2025-10-29T21:29:32Z

T = (Tice-Tsno)*hslyr / (hilyr+hslyr) + Tsno

reduces to

T = (Tice * hslyr + Tsno * hilyr) / (hslyr + hilyr)

Same number of flops but the interpolant is a little easier to recognize.

dabail10 · 2025-10-30T16:59:43Z

T = (Tice-Tsno)*hslyr / (hilyr+hslyr) + Tsno

reduces to

T = (Tice * hslyr + Tsno * hilyr) / (hslyr + hilyr)

Same number of flops but the interpolant is a little easier to recognize.

Which is funny, because I was wondering if I was doing the weighted average correctly. This expression weights the snow temperature with hilyr, whereas my version was weighting it with hslyr.

dabail10 · 2025-10-30T17:52:21Z

Ok. Fixed the computation of Tsnice here. This should be ready to go.

apcraig · 2025-10-31T17:27:22Z

Should I merge, do we need to wait for the CICE PR to also be ready? If we merge this, can we pull this into CICE even if the corresponding CICE PR isn't ready?

dabail10 · 2025-10-31T17:29:15Z

I guess strictly speaking this does not change the icepack_step_therm1 interface and should still work with CICE main.

Fix Tsnice calculation

e9d626f

dabail10 requested review from NickSzapiro-NOAA, apcraig and eclare108213 September 19, 2025 17:54

dabail10 added the Type: Bug label Sep 20, 2025

NickSzapiro-NOAA reviewed Sep 22, 2025

View reviewed changes

apcraig approved these changes Sep 23, 2025

View reviewed changes

anton-seaice reviewed Sep 24, 2025

View reviewed changes

eclare108213 reviewed Sep 24, 2025

View reviewed changes

columnphysics/icepack_therm_vertical.F90 Show resolved Hide resolved

eclare108213 self-assigned this Sep 24, 2025

anton-seaice added a commit to ACCESS-NRI/cice5 that referenced this pull request Sep 24, 2025

Tsnice based on CICE-Consortium/Icepack#542

a8c2d26

dabail10 added 2 commits September 24, 2025 11:01

Fix some bug with Tsnice

1b47cf1

Add zero initialization for Tsnice in driver

453743d

anton-seaice mentioned this pull request Sep 25, 2025

Fixes for sitemptop, sitempbot, and sitempsnic. CICE-Consortium/CICE#1054

Open

18 tasks

apcraig reviewed Sep 30, 2025

View reviewed changes

apcraig self-requested a review September 30, 2025 21:44

NickSzapiro-NOAA approved these changes Oct 1, 2025

View reviewed changes

eclare108213 approved these changes Oct 6, 2025

View reviewed changes

Fix the computation of Tsnice to be a straight line

ba7a2d3

Fix Tsnice calculation #542

Are you sure you want to change the base?

Fix Tsnice calculation #542

Uh oh!

Conversation

dabail10 commented Sep 19, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

PR checklist

Uh oh!

NickSzapiro-NOAA commented Sep 22, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

anton-seaice Sep 24, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

eclare108213 Sep 24, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

dabail10 commented Sep 22, 2025

Uh oh!

apcraig commented Sep 23, 2025

Uh oh!

dabail10 commented Sep 23, 2025

Uh oh!

apcraig commented Sep 23, 2025

Uh oh!

dabail10 commented Sep 23, 2025

Uh oh!

anton-seaice Sep 24, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

apcraig commented Sep 24, 2025

Uh oh!

dabail10 commented Sep 24, 2025

Uh oh!

apcraig left a comment

Choose a reason for hiding this comment

Uh oh!

dabail10 commented Oct 1, 2025

Uh oh!

NickSzapiro-NOAA commented Oct 1, 2025 • edited by dabail10 Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

eclare108213 left a comment

Choose a reason for hiding this comment

Uh oh!

apcraig commented Oct 28, 2025

Uh oh!

dabail10 commented Oct 28, 2025

Uh oh!

dabail10 commented Oct 28, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

dabail10 commented Oct 28, 2025

Uh oh!

eclare108213 commented Oct 29, 2025

Uh oh!

dabail10 commented Oct 29, 2025

Uh oh!

eclare108213 commented Oct 29, 2025

dabail10 commented Sep 19, 2025 •

edited

Loading

NickSzapiro-NOAA commented Sep 22, 2025 •

edited

Loading

anton-seaice Sep 24, 2025 •

edited

Loading

eclare108213 Sep 24, 2025 •

edited

Loading

anton-seaice Sep 24, 2025 •

edited

Loading

NickSzapiro-NOAA commented Oct 1, 2025 •

edited by dabail10

Loading

dabail10 commented Oct 28, 2025 •

edited

Loading