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
merged 8 commits into from
Nov 13, 2025
Merged

Fixes for sitemptop, sitempbot, and sitempsnic. #1054

Show file tree
Hide file tree
Changes from 5 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
112 changes: 59 additions & 53 deletions cicecore/cicedyn/analysis/ice_history.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1535,17 +1535,17 @@ subroutine init_hist (dt)

call define_hist_field(n_sitemptop,"sitemptop","K",tstr2D, tcstr, &
"sea ice surface temperature", &
"none", c1, c0, &
"none", c1, Tffresh, &
ns1, f_sitemptop, avg_ice_present=.true., mask_ice_free_points=.true.)

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

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

call define_hist_field(n_siu,"siu","m/s",ustr2D, ucstr, &
Expand Down Expand Up @@ -2753,22 +2753,18 @@ 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) = aice(i,j,iblk)*trcr(i,j,nt_Tsfc,iblk)
enddo
enddo
call accum_hist_field(n_sitemptop, iblk, worka(:,:), a2D)
endif

! Tsnice is already multiplied by aicen in icepack.
if (f_sitempsnic(1:1) /= 'x') then
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)
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 +2774,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) = aice(i,j,iblk)*Tbot(i,j,iblk)
enddo
enddo
call accum_hist_field(n_sitempbot, iblk, worka(:,:), a2D)
Expand Down Expand Up @@ -3705,33 +3700,40 @@ subroutine accum_hist (dt)
do n = 1, num_avail_hist_fields_2D
if (avail_hist_fields(n)%vhistfreq == histfreq(ns)) then

do j = jlo, jhi
do i = ilo, ihi
if (.not. tmask(i,j,iblk)) then ! mask out land points
a2D(i,j,n,iblk) = spval_dbl
else ! convert units
a2D(i,j,n,iblk) = avail_hist_fields(n)%cona*a2D(i,j,n,iblk) &
* ravgct + avail_hist_fields(n)%conb
endif
enddo ! i
enddo ! j

! Only average for timesteps when ice present
! Only average when/where 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)
if (.not. tmask(i,j,iblk)) then
a2D(i,j,n,iblk) = spval_dbl
else ! convert units
a2D(i,j,n,iblk) = avail_hist_fields(n)%cona*a2D(i,j,n,iblk) &
* ravgip(i,j) + avail_hist_fields(n)%conb
endif
! Mask ice-free points
if (avail_hist_fields(n)%mask_ice_free_points) then
if (ravgip(i,j) == c0) a2D(i,j,n,iblk) = spval_dbl
enddo ! i
enddo ! j
else
do j = jlo, jhi
do i = ilo, ihi
if (.not. tmask(i,j,iblk)) then ! mask out land points
a2D(i,j,n,iblk) = spval_dbl
else ! convert units
a2D(i,j,n,iblk) = avail_hist_fields(n)%cona*a2D(i,j,n,iblk) &
* ravgct + avail_hist_fields(n)%conb
endif
enddo ! i
enddo ! j
endif

! Mask ice-free points
if (avail_hist_fields(n)%mask_ice_free_points) then
do j = jlo, jhi
do i = ilo, ihi
if (ravgip(i,j) == c0) a2D(i,j,n,iblk) = spval_dbl
enddo ! i
enddo ! j
endif

! CMIP albedo: also mask points below horizon
if (index(avail_hist_fields(n)%vname,'sialb') /= 0) then
do j = jlo, jhi
Expand Down Expand Up @@ -3838,30 +3840,33 @@ subroutine accum_hist (dt)
nn = n2D + n
if (avail_hist_fields(nn)%vhistfreq == histfreq(ns)) then

do k = 1, ncat_hist
do j = jlo, jhi
do i = ilo, ihi
if (.not. tmask(i,j,iblk)) then ! mask out land points
a3Dc(i,j,k,n,iblk) = spval_dbl
else ! convert units
a3Dc(i,j,k,n,iblk) = avail_hist_fields(nn)%cona*a3Dc(i,j,k,n,iblk) &
* ravgct + avail_hist_fields(nn)%conb
if (avail_hist_fields(nn)%avg_ice_present) then
do k = 1, ncat_hist
do j = jlo, jhi
do i = ilo, ihi
if (.not. tmask(i,j,iblk)) then ! mask out land points
a3Dc(i,j,k,n,iblk) = spval_dbl
else ! convert units
a3Dc(i,j,k,n,iblk) = avail_hist_fields(nn)%cona*a3Dc(i,j,k,n,iblk) &
* ravgipn(i,j,n) + avail_hist_fields(nn)%conb
endif
enddo ! i
enddo ! j
enddo ! k
else
do k = 1, ncat_hist
do j = jlo, jhi
do i = ilo, ihi
if (.not. tmask(i,j,iblk)) then ! mask out land points
a3Dc(i,j,k,n,iblk) = spval_dbl
else ! convert units
a3Dc(i,j,k,n,iblk) = avail_hist_fields(nn)%cona*a3Dc(i,j,k,n,iblk) &
* ravgct + avail_hist_fields(nn)%conb
endif
enddo ! i
enddo ! j
enddo ! k
endif
enddo ! i
enddo ! j
enddo ! k
if (avail_hist_fields(nn)%avg_ice_present) then
do k = 1, ncat_hist
do j = jlo, jhi
do i = ilo, ihi
if (tmask(i,j,iblk)) then
a3Dc(i,j,k,n,iblk) = &
a3Dc(i,j,k,n,iblk)*avgct(ns)*ravgipn(i,j,k)
endif
enddo ! i
enddo ! j
enddo ! k
endif

endif

Expand All @@ -3885,6 +3890,7 @@ subroutine accum_hist (dt)
enddo ! k
endif
enddo ! n

do n = 1, num_avail_hist_fields_3Db
nn = n3Dzcum + n
if (avail_hist_fields(nn)%vhistfreq == histfreq(ns)) then
Expand Down
12 changes: 9 additions & 3 deletions cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1287,9 +1287,15 @@ subroutine ice_hist_field_def(ncid, hfield, lprecision, dimids, ns)
.and.TRIM(hfield%vname(1:9))/='sistreave' &
.and.TRIM(hfield%vname(1:9))/='sistremax' &
.and.TRIM(hfield%vname(1:4))/='sigP') then
status = nf90_put_att(ncid,varid,'cell_methods','time: mean')
call ice_check_nc(status, subname// ' ERROR: defining cell methods for '//hfield%vname, &
file=__FILE__, line=__LINE__)
if (hfield%avg_ice_present) then
status = nf90_put_att(ncid,varid,'cell_methods','area: time: mean where sea ice (mask=siconc)')
call ice_check_nc(status, subname// ' ERROR: defining cell methods for '//hfield%vname, &
file=__FILE__, line=__LINE__)
else
status = nf90_put_att(ncid,varid,'cell_methods','time: mean')
call ice_check_nc(status, subname// ' ERROR: defining cell methods for '//hfield%vname, &
file=__FILE__, line=__LINE__)
endif
endif
endif

Expand Down
10 changes: 8 additions & 2 deletions cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1428,8 +1428,14 @@ subroutine ice_hist_field_def(File, hfield,lprecision, dimids, ns)
.and.TRIM(hfield%vname(1:9))/='sistreave' &
.and.TRIM(hfield%vname(1:9))/='sistremax' &
.and.TRIM(hfield%vname(1:4))/='sigP') then
call ice_pio_check(pio_put_att(File,varid,'cell_methods','time: mean'), &
subname//' ERROR: defining att cell_methods',file=__FILE__,line=__LINE__)
if (hfield%avg_ice_present) then
call ice_pio_check(pio_put_att(File,varid,'cell_methods', &
'area: time: mean where sea ice (mask=siconc)'), &
subname//' ERROR: defining att cell_methods',file=__FILE__,line=__LINE__)
else
call ice_pio_check(pio_put_att(File,varid,'cell_methods','time: mean'), &
subname//' ERROR: defining att cell_methods',file=__FILE__,line=__LINE__)
endif
endif
endif

Expand Down
2 changes: 1 addition & 1 deletion icepack
Open in desktop
Submodule icepack updated 7 files
+1 −2 columnphysics/icepack_algae.F90
+0 −2 columnphysics/icepack_itd.F90
+3 −4 columnphysics/icepack_therm_vertical.F90
+2 −4 columnphysics/icepack_zbgc_shared.F90
+1 −0 configuration/driver/icedrv_flux.F90
+3 −4 configuration/scripts/tests/poll_queue.csh
+5 −1 doc/source/user_guide/ug_running.rst