Fixes for sitemptop, sitempbot, and sitempsnic. #1054

dabail10 · 2025-09-24T20:19:18Z

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

This PR will also require an icepack commit later on. For now, this is the changes for CICE only. Start of addressing issue #1038 and #1033

NickSzapiro-NOAA · 2025-09-25T19:51:37Z

This isn't bit-for-bit in UFS for sitempsnic and sitempbot. Not that it should be...I just see that box checked.

More substantially, I don't know that we should separate the CICE and Icepack changes for these variables. But maybe I just don't see your plan

dabail10 · 2025-09-25T21:18:29Z

This isn't bit-for-bit in UFS for sitempsnic and sitempbot. Not that it should be...I just see that box checked.

More substantially, I don't know that we should separate the CICE and Icepack changes for these variables. But maybe I just don't see your plan

I guess I think of bfb meaning only the internal variables and not the history variables. In terms of CICE and Icepack, this is the workflow where we have to do each separately and then they come in together.

anton-seaice

This looks ok - but CICE-Consortium/Icepack#542 has changed Tsnice (in Celsius) to already be aice weighted, so in combination with this change , it would get weighted by aice twice (sortof)

dabail10 · 2025-09-26T13:40:36Z

This looks ok - but CICE-Consortium/Icepack#542 has changed Tsnice (in Celsius) to already be aice weighted, so in combination with this change , it would get weighted by aice twice (sortof)

I would be happy to have a discussion about this over Google Meet. Icepack is producing the gridcell Tsnice. So, it is:

Tsnice = sum(aicen(n)*Tsnice).

Technically this can be n = 0,ncat where n=0 is the open water category and Tsnice is zero here. So, this is like dividing by sum(aicen(n)), n = 0,ncat = 1. Then in ice_history.F90 we do sum(aice(t)*Tsnice(t)) / sum(aice(t)) over each timestep in the averaging interval. The sum(aice(t)) division is done later on. So this code:

                    if (tmask(i,j,iblk)) then
                          a2D(i,j,n,iblk) = &
                          a2D(i,j,n,iblk)*avgct(ns)*ravgip(i,j)
                    endif

eclare108213 · 2025-09-26T13:58:52Z

The category and time averaging are confusing enough for everyone (including me) that I think it's worth writing it out mathematically and making that available in the docs. I'll work on that.

My question for @dabail10 is whether the CMIP data request asks for Tsnice over the grid cell or over just the ice. I think it's the latter, and the confusing part is why you're only computing a grid-cell value in Icepack.

I think of the time averaging this way:

Average volume V_avg = sum(V) / M, where the sum is over all time steps in the averaging period and M is the number of time steps with ice present.
Average content, e.g. salt_avg = sum(V*S) / M where S is salinity i.e. salt tracer
Average tracer, e.g. salinity_avg = sum(V*S) / sum(V) = salt_avg / V_avg

The quantities being averaged in time here have already been integrated or averaged over categories, and it's always better to use 'content' quantities for conservation and/or accuracy. I think that's what @dabail10 is doing but I want to write it out with units to check. The temperatures are surface quantities so the integrals would be over the surfaces rather than the volumes.

Let me know if you see errors in my thinking here!

dabail10 · 2025-09-26T16:22:22Z

The category and time averaging are confusing enough for everyone (including me) that I think it's worth writing it out mathematically and making that available in the docs. I'll work on that.

My question for @dabail10 is whether the CMIP data request asks for Tsnice over the grid cell or over just the ice. I think it's the latter, and the confusing part is why you're only computing a grid-cell value in Icepack.

I think of the time averaging this way:

Average volume V_avg = sum(V) / M, where the sum is over all time steps in the averaging period and M is the number of time steps with ice present. Average content, e.g. salt_avg = sum(V*S) / M where S is salinity i.e. salt tracer Average tracer, e.g. salinity_avg = sum(V*S) / sum(V) = salt_avg / V_avg

The quantities being averaged in time here have already been integrated or averaged over categories, and it's always better to use 'content' quantities for conservation and/or accuracy. I think that's what @dabail10 is doing but I want to write it out with units to check. The temperatures are surface quantities so the integrals would be over the surfaces rather than the volumes.

Let me know if you see errors in my thinking here!

The SIMIP request only asks for the temperature within / over / under the ice. However, you point is valid in terms of the definition of intensive versus extensive. I think the point is we want to avoid averaging in zeroes that are not physical. So, the aggregate Tsnice is:

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

From the Notz et al. 2016 paper:

"For all variables that are proportional to area fraction, i.e. extensive variables such as volume, mass, or area fraction, a zero should always be averaged in for all time steps where no sea ice is present. This is because the extensive variables naturally approach zero as area fraction approaches zero."

So in this sense, this definition of Tsnice is indeed extensive. It becomes intensive if we do:

Tsnice_cell = sum(aicen(n)*Tsnicen(n)) / aice(n), n=1,5

this is a weighted mean of Tsnicen, again only over the ice. However, it does not go to zero as aice goes to zero. So, I agree that I should just do a straight time average of Tsnice as it is defined. This is similar for trcr(:,:,nt_Tsfc) and Tbot. For Tbot though since it is all the same for each category, then Tbot_cell = Tbot.

dabail10 · 2025-10-01T16:15:21Z

Here are some images with the updated history averaging and fix to Tsnice in Icepack.

dabail10 · 2025-10-01T16:17:37Z

@eclare108213 @anton-seaice @NickSzapiro-NOAA Hopefully I have addressed your concerns on this as well.

eclare108213

One question about average_ice_present...

eclare108213 · 2025-10-06T17:46:47Z

cicecore/cicedyn/analysis/ice_history.F90

             "sea ice surface temperature", &
             "none", c1, c0,           &
-             ns1, f_sitemptop, avg_ice_present=.true., mask_ice_free_points=.true.)
+             ns1, f_sitemptop, avg_ice_present=.false., mask_ice_free_points=.true.)


Can you explain why avg_ice_present should be false? Since these temperature variables are supposed to be only over the ice area, it seems to me that avg_ice_present should be true. If it is false, then the accumulated field (e.g. a2D) is not multiplied by ravgip, so the field written out is the average over all space and time for that accumulation interval, right?

I keep going around and around on this one. The variable sitemptop is trcr(:,:,nt_Tsfc), so the category aggregate surface temperature. As aicen goes to zero, then trcr(:,:,nt_Tsfc) goes to zero and hence is "extensive". However, this is only true in Celsius and 0C is actually a valid temperature. I would like others interpretation here.

The "extensive" and "intensive" nomenclature is nonintuitive in my opinion. Since nearly every quantity in the ice model is aggregated over categories, ice area is always going to be a factor and those will all go to zero as area goes to zero. The exceptions will be tracers (or similar) for which you want the "actual" value over the ice. I would think that skin temperature would be one of those, and so the aggregated value ought to be divided by ice area. Can you work through an example case that is easy to understand, by hand, following the code? E.g. let's say the ice area changes from 0 to 0.5 to 1 over 3 time steps, and the temperature changes from nothing to something... What should the time averaged value be over the average ice, and does the code get it right? If not, is the problem with avg_ice_present or somewhere else?

However, this is only true in Celsius and 0C is actually a valid temperature.

I think this is kind of coincidental. (We report temperature in K after all)

It clear in Notz 2016 / the CMIP7 data request that temperatures should be masked / averaged only where there is sea ice

anton-seaice · 2025-10-06T22:56:43Z

cicecore/cicedyn/analysis/ice_history.F90

             "sea ice surface temperature", &
             "none", c1, c0,           &
-             ns1, f_sitemptop, avg_ice_present=.true., mask_ice_free_points=.true.)
+             ns1, f_sitemptop, avg_ice_present=.false., mask_ice_free_points=.true.)


Does removing the avg_ice_present and removing if (aice(i,j,iblk) > puny) later on mean that this field is not masked over land anymore ?

Nevermind, tmask is applied seperately

anton-seaice · 2025-10-07T04:39:13Z

Sorry @dabail10 - I am still confused, well probably need to have a meeting !

In my very simply example, with one ice thickness category:

If ncat = 1 and aicen = aice = 0.5
And there is no snow
And surface temperature is 2°C

Then icepack produces a grid cell mean temperature (Tsnice) of 1°C (as sum(aicen*Tsf))

If we assume that there is no change during a day, then the implementation in this PR calculates the sitemptop daily average as 274.15K (for either setting of avg_ice_present)

But the correct answer is 275.15K ? (As its a diagnostic of surface temperature over sea ice)

anton-seaice · 2025-10-10T03:15:32Z

I realised in our discussion this morning we forgot to talk about the right step the C->K conversion should happen (:

eclare108213 · 2025-10-10T12:28:26Z

the right step the C->K conversion should happen

I think that should be last, to just do the calculation once rather than over and over during the time accumulation.

dabail10 · 2025-10-10T13:10:28Z

the right step the C->K conversion should happen

I think that should be last, to just do the calculation once rather than over and over during the time accumulation.

I always forget that we can add a conversion factor in the variable definition. This is added at the end.

anton-seaice · 2025-10-14T01:22:20Z

cicecore/cicedyn/analysis/ice_history.F90

           do i = ilo, ihi
-              if (aice(i,j,iblk) > puny) &
-              worka(i,j) = aice(i,j,iblk)*(trcr(i,j,nt_Tsfc,iblk)+Tffresh)
+              worka(i,j) = trcr(i,j,nt_Tsfc,iblk)+Tffresh


Is it correct to use trcr(:,:,nt_Tsfc,:) here?

I believe trcr(:,:,nt_Tsfc,:) is a grid cell average, where the open water area is Tfrz.

But what we want is a ice area aggregate only of Tsfc?

e.g.

sum(aicen(n)*trcrn(nt_Tsfc,n)) for n=1,ncat

without including the open water temperature

This is a good question. My assumption was that all trcr tracers were just sum(aicen*trcrn) from n=1,5. Is this not the case?

I found this quite hard to figure out, but I think Tf is averaged in.

There is this note for the Tsfc diagnostic:

CICE/cicecore/cicedyn/analysis/ice_history.F90

Lines 757 to 760 in bd30d63

call define_hist_field(n_Tsfc,"Tsfc","C",tstr2D, tcstr, &

"snow/ice surface temperature", &

"averaged with Tf if no ice is present", c1, c0, &

ns1, f_Tsfc)

And this special handling for Tsfc:

if (aicen > puny) then trcrn(it) = atrcrn(it) / aicen else trcrn(it) = c0 if (it == nt_Tsfc) then trcrn(it) = Tf ! surface temperature endif

from https://github.com/CICE-Consortium/Icepack/blob/7ee2266089f80effb5fee8335c4a1056d7da11e7/columnphysics/icepack_tracers.F90#L1288

Right, so this means trcrn(:,:,nt_Tsfc) is only set to Tf when aicen(n) < puny. So, this would not be averaged in. Here is the code in icepack_aggregate. I suppose this might still average values of Tf when aicen(n) < puny. @eclare108213 is this correct?

do n = 1, ncat aice = aice + aicen(n) vice = vice + vicen(n) vsno = vsno + vsnon(n) do it = 1, ntrcr atrcrn = trcrn(it,n)*(trcr_base(it,1) * aicen(n) & + trcr_base(it,2) * vicen(n) & + trcr_base(it,3) * vsnon(n)) if (n_trcr_strata(it) > 0) then ! additional tracer layers do itl = 1, n_trcr_strata(it) ntr = nt_strata(it,itl) atrcrn = atrcrn * trcrn(ntr,n) enddo endif atrcr(it) = atrcr(it) + atrcrn enddo ! ntrcr enddo ! ncat

Oh ok - not as problematic as I was thinking. Probably ok then !

anton-seaice · 2025-10-21T04:41:21Z

the right step the C->K conversion should happen

I think that should be last, to just do the calculation once rather than over and over during the time accumulation.

I always forget that we can add a conversion factor in the variable definition. This is added at the end.

I just realised we might need to change the order of operations to do this, the masking for "avg_ice_present" being true occurs after "conb" (the constant to add to a diagnostic) is added

see

CICE/cicecore/cicedyn/analysis/ice_history.F90

Line 3720 in bd30d63

if (avail_hist_fields(n)%avg_ice_present) then

dabail10 · 2025-10-21T19:43:31Z

I remember this code originally. I think I need to change the logic. The following two lines are accumulating normally:

                    a2D(i,j,n,iblk) = avail_hist_fields(n)%cona*a2D(i,j,n,iblk) &
                                   * ravgct + avail_hist_fields(n)%conb

where the multiplication of ravgct is a division by the number of steps in the averaging period. However, the next set of lines in the avg_ice_present if block are:

                          a2D(i,j,n,iblk) = &
                          a2D(i,j,n,iblk)*avgct(ns)*ravgip(i,j)

Are undoing the multiplication by ravgct by mutiplying by avgct(ns) and then dividing by the accumulated ice area (ravgip). I think this is a bug. I should probably only do the first part when avg_ice_present is .false. and then change the second line to:

                      a2D(i,j,n,iblk) = avail_hist_fields(n)%cona*a2D(i,j,n,iblk) &
                                   * ravgip(i,j) + avail_hist_fields(n)%conb

this was how the pond variables were coded. Again, I think I did this. Not sure why I did it this way.

…hods for variables when avg_ice_present is true.

dabail10 · 2025-10-23T15:27:32Z

Ok. I am happy with the averaging when ice present for sitemptop, sitempsnic, and sitempbot. I also added the updated cell_methods for variables when avg_ice_present = .true. I will work on the aice, aice_mid, aice_init stuff in the next PR.

dabail10 · 2025-10-23T15:28:05Z

@eclare108213 @anton-seaice @NickSzapiro-NOAA @apcraig This is ready for review again.

anton-seaice

Thanks Dave !

cicecore/cicedyn/analysis/ice_history.F90

cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90

cicecore/cicedyn/analysis/ice_diagnostics.F90

icepack

cicecore/cicedyn/analysis/ice_history.F90

apcraig

Why is Icepack updated in this PR? The PR says Icepack is not being updated. Does this PR need the latest Icepack and/or CICE-Consortium/Icepack#542?

Are there any changes needed for io_binary (noticed only io_netcdf and io_pio2) were updated?

I have not checked the science of the averaging, just the technical details.

apcraig · 2025-10-28T18:55:04Z

cicecore/cicedyn/analysis/ice_history.F90

           do i = ilo, ihi
-              if (aice_init(i,j,iblk) > puny) &
-                 worka(i,j) = aice(i,j,iblk)*(Tbot(i,j,iblk)/aice_init(i,j,iblk)+Tffresh)
+               worka(i,j) = aice(i,j,iblk)*Tbot(i,j,iblk)


could shift this line 1 character left

This is just metadata in the netcdf files. We don't have this for binary files?

dabail10 · 2025-10-28T19:06:09Z

We need the Icepack fix for the computation of Tsnice. I assumed this woul come in later when the Icepack PR has been merged. Am I missing something?

apcraig · 2025-10-28T22:10:28Z

We need the Icepack fix for the computation of Tsnice. I assumed this woul come in later when the Icepack PR has been merged. Am I missing something?

That's fine, so we do need to merge the Icepack PR first. Then you will want to update this PR, right?

Fixes for sitemptop, sitempbot, and sitempsnic.

ae5b8fd

dabail10 requested review from NickSzapiro-NOAA, anton-seaice, apcraig and eclare108213 September 24, 2025 20:19

NickSzapiro-NOAA approved these changes Sep 25, 2025

View reviewed changes

anton-seaice reviewed Sep 25, 2025

View reviewed changes

Change the averaging for the sitemp variables

0dfcba1

eclare108213 mentioned this pull request Sep 30, 2025

documentation for averaging #1056

Open

18 tasks

dabail10 added the Type: Bug label Oct 1, 2025

dabail10 self-assigned this Oct 1, 2025

Fix comment

13de839

eclare108213 reviewed Oct 6, 2025

View reviewed changes

anton-seaice reviewed Oct 6, 2025

View reviewed changes

anton-seaice reviewed Oct 14, 2025

View reviewed changes

Address issues with averaging when ice present. Also, change cell_met…

1bd92f4

…hods for variables when avg_ice_present is true.

anton-seaice suggested changes Oct 23, 2025

View reviewed changes

Cleanup the code a bit more

069d89c

apcraig requested changes Oct 28, 2025

View reviewed changes

	call define_hist_field(n_Tsfc,"Tsfc","C",tstr2D, tcstr, &
	"snow/ice surface temperature", &
	"averaged with Tf if no ice is present", c1, c0, &
	ns1, f_Tsfc)

Uh oh!

Fixes for sitemptop, sitempbot, and sitempsnic. #1054

Are you sure you want to change the base?

Fixes for sitemptop, sitempbot, and sitempsnic. #1054

Conversation

dabail10 commented Sep 24, 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 25, 2025

Uh oh!

dabail10 commented Sep 25, 2025

Uh oh!

anton-seaice left a comment • 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 26, 2025

Uh oh!

eclare108213 commented Sep 26, 2025

Uh oh!

dabail10 commented Sep 26, 2025

Uh oh!

dabail10 commented Oct 1, 2025

Uh oh!

dabail10 commented Oct 1, 2025

Uh oh!

eclare108213 left a comment

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 Oct 6, 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!

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

Uh oh!

Uh oh!

anton-seaice commented Oct 10, 2025

Uh oh!

eclare108213 commented Oct 10, 2025

Uh oh!

dabail10 commented Oct 10, 2025

Uh oh!

anton-seaice Oct 14, 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!

anton-seaice Oct 14, 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!

anton-seaice commented Oct 21, 2025

Uh oh!

dabail10 commented Oct 21, 2025

Uh oh!

dabail10 commented Oct 23, 2025

Uh oh!

dabail10 commented Oct 23, 2025

Uh oh!

anton-seaice left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

dabail10 commented Sep 24, 2025 •

edited

Loading

anton-seaice left a comment •

edited

Loading

anton-seaice Oct 6, 2025 •

edited

Loading

anton-seaice commented Oct 7, 2025 •

edited

Loading

anton-seaice Oct 14, 2025 •

edited

Loading

anton-seaice Oct 14, 2025 •

edited

Loading