Skip to content

Fixes for sitemptop, sitempbot, and sitempsnic. #1054

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.

Sign up for GitHub

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

Merged
apcraig merged 8 commits into from
Nov 13, 2025
Merged

Fixes for sitemptop, sitempbot, and sitempsnic. #1054

Changes from 2 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
ae5b8fd
Fixes for sitemptop, sitempbot, and sitempsnic.
dabail10 Sep 24, 2025
0dfcba1
Change the averaging for the sitemp variables
dabail10 Sep 26, 2025
13de839
Fix comment
dabail10 Oct 6, 2025
1bd92f4
Address issues with averaging when ice present. Also, change cell_met…
dabail10 Oct 23, 2025
069d89c
Cleanup the code a bit more
dabail10 Oct 24, 2025
fea53d7
Fix indent
dabail10 Oct 30, 2025
7552bd1
Update icepack with Tsnice fix
dabail10 Nov 3, 2025
58de585
Fix syntax error in ice_history.F90
dabail10 Nov 5, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 6 additions & 12 deletions cicecore/cicedyn/analysis/ice_history.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1536,17 +1536,17 @@ subroutine init_hist (dt)
call define_hist_field(n_sitemptop,"sitemptop","K",tstr2D, tcstr, &
"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.)
Copy link
Contributor

@eclare108213 eclare108213 Oct 6, 2025

Choose a reason for hiding this comment

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

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?

Copy link
Contributor Author

@dabail10 dabail10 Oct 6, 2025

Choose a reason for hiding this comment

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

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.

Copy link
Contributor

@eclare108213 eclare108213 Oct 6, 2025

Choose a reason for hiding this comment

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

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?

Copy link
Contributor

@anton-seaice anton-seaice Oct 6, 2025

Choose a reason for hiding this comment

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

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

Copy link
Contributor

@anton-seaice anton-seaice Oct 6, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

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 ?

Copy link
Contributor

@anton-seaice anton-seaice Oct 7, 2025

Choose a reason for hiding this comment

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

Nevermind, tmask is applied seperately


call define_hist_field(n_sitempsnic,"sitempsnic","K",tstr2D, tcstr, &
"snow ice interface temperature", &
"surface temperature when no snow present", c1, c0, &
ns1, f_sitempsnic, avg_ice_present=.true., mask_ice_free_points=.true.)
ns1, f_sitempsnic, avg_ice_present=.false., mask_ice_free_points=.true.)

call define_hist_field(n_sitempbot,"sitempbot","K",tstr2D, tcstr, &
"sea ice bottom temperature", &
"none", c1, c0, &
ns1, f_sitempbot, avg_ice_present=.true., mask_ice_free_points=.true.)
ns1, f_sitempbot, avg_ice_present=.false., mask_ice_free_points=.true.)

call define_hist_field(n_siu,"siu","m/s",ustr2D, ucstr, &
"ice x velocity component", &
Expand Down Expand Up @@ -2753,8 +2753,7 @@ subroutine accum_hist (dt)
worka(:,:) = c0
do j = jlo, jhi
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
Copy link
Contributor

@anton-seaice anton-seaice Oct 14, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

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

Copy link
Contributor Author

@dabail10 dabail10 Oct 14, 2025

Choose a reason for hiding this comment

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

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?

Copy link
Contributor

@anton-seaice anton-seaice Oct 14, 2025

Choose a reason for hiding this comment

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

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

Copy link
Contributor Author

@dabail10 dabail10 Oct 14, 2025

Choose a reason for hiding this comment

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

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

Copy link
Contributor

@anton-seaice anton-seaice Oct 14, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

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

enddo
enddo
call accum_hist_field(n_sitemptop, iblk, worka(:,:), a2D)
Expand All @@ -2764,11 +2763,7 @@ subroutine accum_hist (dt)
worka(:,:) = c0
do j = jlo, jhi
do i = ilo, ihi
if (vsno(i,j,iblk) > puny .and. aice_init(i,j,iblk) > puny) then
worka(i,j) = aice(i,j,iblk)*(Tsnice(i,j,iblk)/aice_init(i,j,iblk)+Tffresh)
else
worka(i,j) = aice(i,j,iblk)*(trcr(i,j,nt_Tsfc,iblk)+Tffresh)
endif
worka(i,j) = Tsnice(i,j,iblk)+Tffresh
enddo
enddo
call accum_hist_field(n_sitempsnic, iblk, worka(:,:), a2D)
Comment on lines 2763 to 2770
Copy link
Contributor

@anton-seaice anton-seaice Oct 23, 2025

Choose a reason for hiding this comment

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

Suggested change
call accum_hist_field(n_sitempsnic, iblk, worka(:,:), a2D)
if (f_sitempsnic(1:1) /= 'x') &
call accum_hist_field(n_sitempsnic, iblk, Tsnice(:,:,iblk), a2D)

is it useful to add a comment here,

e.g. sitempsnic is intensive, and Tsnice is already weighted by aice in icepack

Copy link
Contributor

@eclare108213 eclare108213 Nov 3, 2025

Choose a reason for hiding this comment

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

Why is Tsnice treated differently from the other surface temperatures? This is confusing. Could the multiplication by aice be done here instead?

Copy link
Contributor Author

@dabail10 dabail10 Nov 3, 2025

Choose a reason for hiding this comment

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

I guess I am not sure what you are saying. Tsnice is computed in step_therm1 and aggregated using aicen_init. Whereas Tsfc is advected and changed after the dynamics. Also, trcr(:,:,nt_Tsfc) is divided by sum(aicen(:)).

Copy link
Contributor

@eclare108213 eclare108213 Nov 3, 2025

Choose a reason for hiding this comment

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

In this section of code, the Tsfc tracer and Tbot are multiplied by aice, but Tsnice is not. I think it would be less confusing if all three were sent out of Icepack in the same form. I know they are computed differently in Icepack, but shouldn't the end result that is sent to the driver be something like "a surface temperature (top, ice-snow interface, or bottom) that is averaged only over the ice that's present" for all three? At the moment it looks like Tsfc and Tbot are averaged over the ice and Tsnice is averaged over the grid cell when they are sent out of Icepack.

Copy link
Contributor Author

@dabail10 dabail10 Nov 3, 2025

Choose a reason for hiding this comment

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

In this section of code, the Tsfc tracer and Tbot are multiplied by aice, but Tsnice is not. I think it would be less confusing if all three were sent out of Icepack in the same form. I know they are computed differently in Icepack, but shouldn't the end result that is sent to the driver be something like "a surface temperature (top, ice-snow interface, or bottom) that is averaged only over the ice that's present" for all three? At the moment it looks like Tsfc and Tbot are averaged over the ice and Tsnice is averaged over the grid cell when they are sent out of Icepack.

I feel like I am going around in circles on this one. The idea was that we did not want Tsnice as:

Tsnice = sum(aicen(n)*Tsnicen(n)) / sum(aicen(n))

and then multiply by aice again in the history averaging. So, Icepack is just sending:

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

thus we do not need to mutilply again by aice when accumulating in ice_history.F90. However Tsfc and Tbot do require multiplication by aice.

Copy link
Contributor

@eclare108213 eclare108213 Nov 3, 2025

Choose a reason for hiding this comment

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

Then Tsfc and Tbot are the ones that ought to be fixed. Let's plan to clean this up when we fix the general aice and aice_init issues rather than in this PR. Thanks for the clarification.
There are a lot of "will do so-and-so in a later PR" in the comments here -- need to keep track of them in an issue (new or updates to an existing one)...

Expand All @@ -2778,8 +2773,7 @@ subroutine accum_hist (dt)
worka(:,:) = c0
do j = jlo, jhi
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) = Tbot(i,j,iblk)+Tffresh
enddo
enddo
call accum_hist_field(n_sitempbot, iblk, worka(:,:), a2D)
Expand Down