Skip to content

Add conductive flux limiting for vertical thermodynamics when calc_Tsfc=.false. #557

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

Open
kieranricardo wants to merge 37 commits into
base: main
Choose a base branch
from
Open

Add conductive flux limiting for vertical thermodynamics when calc_Tsfc=.false. #557

Show file tree
Hide file tree
Changes from 34 commits
Commits
Show all changes
37 commits
Select commit Hold shift + click to select a range
bc088ac
logging, tweaked newton solver tolerances, re-add einex
kieranricardo Apr 29, 2024
ea4cde3
fix argument ordering in set_sfcflux
kieranricardo May 6, 2024
d4b4a31
remove cm2 style conductive flux limiting
kieranricardo May 7, 2024
a15fc03
use cm2 style conductive flux limiting
kieranricardo May 20, 2024
65259a0
change flux energy error threshold, and reduce fsurf before thermo solve
kieranricardo Jul 22, 2024
e3b9746
don't reduce fsurf before thermo solve - has no effect
kieranricardo Jul 23, 2024
2cf33b1
convert GBM ice fluxes to local
kieranricardo Aug 20, 2024
af43525
remove scaling by inverse ice fraction
kieranricardo Aug 26, 2024
d1ff164
increase ice and snow thickness minimum values to match CM2, revert c…
kieranricardo Sep 17, 2024
e6e429c
scale ice fluxes by changing ice area
kieranricardo Sep 17, 2024
07352ba
Remove stray deprecated zsal flag
blimlim Mar 31, 2025
e8ef69f
Change default hi_min to p1
blimlim Apr 9, 2025
9bd1b3f
Remove unnecessary comments
blimlim Apr 23, 2025
140e8cb
Review suggestions: remove unnecessary include and comment
blimlim Apr 24, 2025
bf0995b
white space changes, remove potentiall unescessary calc_Tsfc statement
kieranricardo Feb 3, 2026
5d2160d
fix edit
kieranricardo Feb 3, 2026
6782499
calculate energy error differently for calc_tsfc options
kieranricardo Feb 3, 2026
ed20c15
white space change
kieranricardo Feb 3, 2026
7036d9c
revert snow and ice thickness limits to default values
kieranricardo Feb 3, 2026
10faac6
white space change
kieranricardo Feb 3, 2026
789f982
make hs_min configurable
kieranricardo Feb 3, 2026
d2c2a09
revert configurable hs_min
kieranricardo Feb 17, 2026
1fb991d
remove tab
kieranricardo Feb 17, 2026
8516838
Update columnphysics/icepack_therm_vertical.F90
kieranricardo Feb 24, 2026
594f9ff
add extra calc_Tsfc, and add cap_conductive_flux parameters to icepac…
kieranricardo Feb 24, 2026
6ad9e4c
add ACCESS flag
kieranricardo Feb 24, 2026
e92ca42
add ACCESS flag
kieranricardo Feb 24, 2026
bfd3c4c
extra .not. calc_Tsfc check
kieranricardo Feb 25, 2026
7e86b82
extra .not. calc_Tsfc check
kieranricardo Feb 25, 2026
9b8a3ff
fix fcondtop issue
kieranricardo Feb 25, 2026
9e0e644
Update columnphysics/icepack_parameters.F90
kieranricardo Mar 24, 2026
79cef70
add .not. calc_Tsfc check for thermo conservation error logging
kieranricardo Mar 24, 2026
ebed003
formatting and cleanup
kieranricardo Apr 14, 2026
ed835b8
update einex and e_num variable names
kieranricardo May 11, 2026
cda0f14
comments with units for new variables + formatting
kieranricardo Jun 16, 2026
edbb4ce
comments with units for new variables + formatting
kieranricardo Jun 16, 2026
a1d5fac
formatting
kieranricardo Jun 16, 2026
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
2 changes: 1 addition & 1 deletion columnphysics/icepack_flux.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -351,7 +351,7 @@ subroutine set_sfcflux (aicen, &

raicen = c1

#ifdef CICE_IN_NEMO
#if defined(CICE_IN_NEMO) || defined(CICE_ACCESS3)
!----------------------------------------------------------------------
! Convert fluxes from GBM values to per ice area values when
! running in NEMO environment. (When in standalone mode, fluxes
Expand Down
27 changes: 20 additions & 7 deletions columnphysics/icepack_parameters.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,9 @@ module icepack_parameters
aspect_rapid_mode = 1.0_dbl_kind,&! aspect ratio (larger is wider)
dSdt_slow_mode = -1.5e-7_dbl_kind,&! slow mode drainage strength (m s-1 K-1)
phi_c_slow_mode = 0.05_dbl_kind,&! critical liquid fraction porosity cutoff
phi_i_mushy = 0.85_dbl_kind ! liquid fraction of congelation ice
phi_i_mushy = 0.85_dbl_kind,&! liquid fraction of congelation ice
ratio_Wm2_m = 1000.0_dbl_kind ,&! max condutive flux/depth ratio (W m)
cold_temp_flag = -60.0_dbl_kind ! min temp used to limit the conductive flux (C)

integer (kind=int_kind), public :: &
ktherm = 1 ! type of thermodynamics
Expand Down Expand Up @@ -585,8 +587,8 @@ subroutine icepack_init_parameters( &
cpl_frazil_in, semi_implicit_Tsfc_in, vapor_flux_correction_in, &
Rac_rapid_mode_in, aspect_rapid_mode_in, &
dSdt_slow_mode_in, phi_c_slow_mode_in, &
phi_i_mushy_in, shortwave_in, albedo_type_in, albsnowi_in, &
albicev_in, albicei_in, albsnowv_in, &
phi_i_mushy_in, ratio_Wm2_m_in, cold_temp_flag_in, shortwave_in, albedo_type_in, &
albsnowi_in, albicev_in, albicei_in, albsnowv_in, &
ahmax_in, R_ice_in, R_pnd_in, R_snw_in, dT_mlt_in, rsnw_mlt_in, &
kalg_in, R_gC2molC_in, kstrength_in, krdg_partic_in, krdg_redist_in, mu_rdg_in, &
atmbndy_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, &
Expand Down Expand Up @@ -732,7 +734,9 @@ subroutine icepack_init_parameters( &
aspect_rapid_mode_in , & ! aspect ratio for rapid drainage mode (larger=wider)
dSdt_slow_mode_in , & ! slow mode drainage strength (m s-1 K-1)
phi_c_slow_mode_in , & ! liquid fraction porosity cutoff for slow mode
phi_i_mushy_in ! liquid fraction of congelation ice
phi_i_mushy_in , & ! liquid fraction of congelation ice
ratio_Wm2_m_in , & ! max condutive flux/depth ratio (W m)
cold_temp_flag_in ! min temp used to limit the conductive flux (C)

character(len=*), intent(in), optional :: &
congel_freeze_in ! congelation computation
Expand Down Expand Up @@ -1205,6 +1209,8 @@ subroutine icepack_init_parameters( &
if (present(dSdt_slow_mode_in) ) dSdt_slow_mode = dSdt_slow_mode_in
if (present(phi_c_slow_mode_in) ) phi_c_slow_mode = phi_c_slow_mode_in
if (present(phi_i_mushy_in) ) phi_i_mushy = phi_i_mushy_in
if (present(ratio_Wm2_m_in) ) ratio_Wm2_m = ratio_Wm2_m_in
if (present(cold_temp_flag_in) ) cold_temp_flag = cold_temp_flag_in
Comment on lines +1212 to +1213

@anton-seaice anton-seaice Mar 25, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

@kieranricardo has added these as configurable - but we don't plan to add them to the cice namelist (as we dont foresee changing them).

Are we happy to change answers when calc_Tsfc = .false. - @apcraig ? If so, should we not make these two new parameters configurable through icepack_init_parameters ?

@eclare108213 eclare108213 Apr 6, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Are there values of these parameters that would preserve BFBness when calc_Tsfc = .false.? If so, I recommend setting the code defaults to those values and allow them to be changed via namelist.

@kieranricardo kieranricardo May 4, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

We could pick extreme value to ensure the cap_conductive_fluxroutine doesn't do anything, but because the flux is also limited in bl99 with e_num this will still change answers when calc_Tsfc = .false..

if (present(shortwave_in) ) shortwave = shortwave_in
if (present(albedo_type_in) ) albedo_type = albedo_type_in
if (present(albicev_in) ) albicev = albicev_in
Expand Down Expand Up @@ -1598,8 +1604,9 @@ subroutine icepack_query_parameters( &
Lfresh_out, cprho_out, Cp_out, ustar_min_out, hi_min_out, a_rapid_mode_out, &
ktherm_out, conduct_out, fbot_xfer_type_out, calc_Tsfc_out, &
Rac_rapid_mode_out, aspect_rapid_mode_out, dSdt_slow_mode_out, &
phi_c_slow_mode_out, phi_i_mushy_out, shortwave_out, semi_implicit_Tsfc_out, &
albedo_type_out, albicev_out, albicei_out, albsnowv_out, vapor_flux_correction_out, &
phi_c_slow_mode_out, phi_i_mushy_out, ratio_Wm2_m_out, cold_temp_flag_out, &
shortwave_out, semi_implicit_Tsfc_out, albedo_type_out, albicev_out, &
albicei_out, albsnowv_out, vapor_flux_correction_out, &
albsnowi_out, ahmax_out, R_ice_out, R_pnd_out, R_snw_out, dT_mlt_out, &
rsnw_mlt_out, dEdd_algae_out, &
kalg_out, R_gC2molC_out, kstrength_out, krdg_partic_out, krdg_redist_out, mu_rdg_out, &
Expand Down Expand Up @@ -1754,7 +1761,9 @@ subroutine icepack_query_parameters( &
aspect_rapid_mode_out , & ! aspect ratio for rapid drainage mode (larger=wider)
dSdt_slow_mode_out , & ! slow mode drainage strength (m s-1 K-1)
phi_c_slow_mode_out , & ! liquid fraction porosity cutoff for slow mode
phi_i_mushy_out ! liquid fraction of congelation ice
phi_i_mushy_out , & ! liquid fraction of congelation ice
ratio_Wm2_m_out , & ! max condutive flux/depth ratio (W m)
cold_temp_flag_out ! min temp used to limit the conductive flux (C)

character(len=*), intent(out), optional :: &
congel_freeze_out ! congelation computation
Expand Down Expand Up @@ -2261,6 +2270,8 @@ subroutine icepack_query_parameters( &
if (present(dSdt_slow_mode_out) ) dSdt_slow_mode_out = dSdt_slow_mode
if (present(phi_c_slow_mode_out) ) phi_c_slow_mode_out = phi_c_slow_mode
if (present(phi_i_mushy_out) ) phi_i_mushy_out = phi_i_mushy
if (present(ratio_Wm2_m_out) ) ratio_Wm2_m_out = ratio_Wm2_m
if (present(cold_temp_flag_out) ) cold_temp_flag_out = cold_temp_flag
if (present(shortwave_out) ) shortwave_out = shortwave
if (present(albedo_type_out) ) albedo_type_out = albedo_type
if (present(albicev_out) ) albicev_out = albicev
Expand Down Expand Up @@ -2575,6 +2586,8 @@ subroutine icepack_write_parameters(iounit)
write(iounit,*) " dSdt_slow_mode = ", dSdt_slow_mode
write(iounit,*) " phi_c_slow_mode = ", phi_c_slow_mode
write(iounit,*) " phi_i_mushy= ", phi_i_mushy
write(iounit,*) " ratio_Wm2_m= ", ratio_Wm2_m
write(iounit,*) " cold_temp_flag = ", cold_temp_flag

@anton-seaice anton-seaice Mar 24, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

should these be within a if ( .not. calc_Tsfc) ?

@kieranricardo kieranricardo Mar 24, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

It looks like there's no other if checks in icepack_write_parameters, so I think all parameters are output regardless of whether they're used. It might be inconsistent to add an if statement here.

@eclare108213 eclare108213 Apr 6, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

this is okay as long as they are not used elsewhere except when calc_Tsfc = .false.

write(iounit,*) " shortwave = ", trim(shortwave)
write(iounit,*) " albedo_type= ", trim(albedo_type)
write(iounit,*) " albicev = ", albicev
Expand Down
86 changes: 63 additions & 23 deletions columnphysics/icepack_therm_bl99.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ subroutine temperature_changes (dt, &
fsensn, flatn, &
flwoutn, fsurfn, &
fcondtopn,fcondbot, &
einit )
einit, einex_sfc_flux)

real (kind=dbl_kind), intent(in) :: &
dt ! time step
Expand Down Expand Up @@ -122,6 +122,9 @@ subroutine temperature_changes (dt, &
zqsn , & ! snow layer enthalpy (J m-3)
zTsn ! internal snow layer temperatures

real (kind=dbl_kind), intent(out):: &
einex_sfc_flux ! excess energy from conductive flux (J m-2)

! local variables

integer (kind=int_kind), parameter :: &
Expand Down Expand Up @@ -155,7 +158,7 @@ subroutine temperature_changes (dt, &
dflat_dT , & ! deriv of flat wrt Tsf (W m-2 deg-1)
dflwout_dT , & ! deriv of flwout wrt Tsf (W m-2 deg-1)
dt_rhoi_hlyr, & ! dt/(rhoi*hilyr)
einex , & ! excess energy from dqmat to ocean
einex_sfc_calc, & ! excess energy from dqmat to ocean
ferr ! energy conservation error (W m-2)

real (kind=dbl_kind), dimension (nilyr) :: &
Expand Down Expand Up @@ -189,20 +192,24 @@ subroutine temperature_changes (dt, &
Iswabs_tmp , & ! energy to melt through fraction frac of layer
Sswabs_tmp , & ! same for snow
dswabs , & ! difference in swabs and swabs_tmp
frac
frac , &
fcondtopn_reduction, &
fcondtopn_force, dqmat_sn

logical (kind=log_kind) :: &
converged ! = true when local solution has converged
converged, Top_T_was_reset_last_time ! = true when local solution has converged

logical (kind=log_kind) , dimension (nilyr) :: &
reduce_kh ! reduce conductivity when T exceeds Tmlt
reduce_kh ! reduce conductivity when T exceeds Tmlt

character(len=*),parameter :: subname='(temperature_changes)'

!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------

fcondtopn_reduction = c0
einex_sfc_flux = c0
Top_T_was_reset_last_time = .false.
converged = .false.
l_snow = .false.
l_cold = .true.
Expand All @@ -212,7 +219,7 @@ subroutine temperature_changes (dt, &
dfsens_dT = c0
dflat_dT = c0
dflwout_dT = c0
einex = c0
einex_sfc_calc = c0
if (semi_implicit_Tsfc) then ! initialize
dfsurf_dT = dfsurfdTs_cpl
dflat_dT = dflatdTs_cpl
Expand Down Expand Up @@ -337,7 +344,7 @@ subroutine temperature_changes (dt, &
if (.not.semi_implicit_Tsfc) dfsurf_dT = c0
avg_Tsi = c0
enew = c0
einex = c0
einex_sfc_calc = c0

!-----------------------------------------------------------------
! Update specific heat of ice layers.
Expand Down Expand Up @@ -428,6 +435,7 @@ subroutine temperature_changes (dt, &

else

fcondtopn_force = fcondtopn - fcondtopn_reduction
call get_matrix_elements_know_Tsfc ( &
l_snow, Tbot, &
Tin_init, Tsn_init, &
Expand All @@ -436,7 +444,7 @@ subroutine temperature_changes (dt, &
etai, etas, &
sbdiag, diag, &
spdiag, rhs, &
fcondtopn)
fcondtopn_force)
if (icepack_warnings_aborted(subname)) return

endif ! calc_Tsfc
Expand Down Expand Up @@ -558,7 +566,33 @@ subroutine temperature_changes (dt, &
else
zTsn(k) = c0
endif
if (l_brine) zTsn(k) = min(zTsn(k), c0)
if ((l_brine) .and. zTsn(k)>c0) then

if (.not. calc_Tsfc) then
! return this energy to the ocean

dqmat_sn = (zTsn(k)*cp_ice - Lfresh)*rhos - zqsn(k)

! If this is the second time in succession that Tsn(1) has been
! reset, tell the solver to reduce the forcing at the top, and
! pass the difference to the array enum where it will eventually
! go into the ocean
! This is done to avoid an 'infinite loop' whereby temp continually evolves
! to the same point above zero, is reset, ad infinitum
if (l_snow .AND. k == 1) then
if (Top_T_was_reset_last_time) then
fcondtopn_reduction = fcondtopn_reduction + dqmat_sn*hslyr / dt
Top_T_was_reset_last_time = .false.
einex_sfc_flux = einex_sfc_flux + hslyr * dqmat_sn
else
Top_T_was_reset_last_time = .true.
endif
endif
end if

zTsn(k) = min(zTsn(k), c0)

endif

!-----------------------------------------------------------------
! If condition 1 or 2 failed, average new snow layer
Expand Down Expand Up @@ -592,6 +626,16 @@ subroutine temperature_changes (dt, &
dTmat(k) = zTin(k) - Tmlts(k)
dqmat(k) = rhoi * dTmat(k) &
* (cp_ice - Lfresh * Tmlts(k)/zTin(k)**2)

if ((.not. calc_Tsfc) .and. (.not. l_snow) .and. (k == 1)) then
if (Top_T_was_reset_last_time) then
fcondtopn_reduction = fcondtopn_reduction + dqmat(k)*hilyr / dt
Top_T_was_reset_last_time = .false.
einex_sfc_flux = einex_sfc_flux + hilyr * dqmat(k)
else
Top_T_was_reset_last_time = .true.
endif
endif
! use this for the case that Tmlt changes by an amount dTmlt=Tmltnew-Tmlt(k)
! + rhoi * dTmlt &
! * (cp_ocn - cp_ice + Lfresh/zTin(k))
Expand Down Expand Up @@ -637,7 +681,7 @@ subroutine temperature_changes (dt, &
zqin(k) = -rhoi * (-cp_ice*zTin(k) + Lfresh)
endif
enew = enew + hilyr * zqin(k)
einex = einex + hilyr * dqmat(k)
einex_sfc_calc = einex_sfc_calc + hilyr * dqmat(k)

Tin_start(k) = zTin(k) ! for next iteration

Expand Down Expand Up @@ -683,11 +727,15 @@ subroutine temperature_changes (dt, &
fcondbot = kh(1+nslyr+nilyr) * &
(zTin(nilyr) - Tbot)

! Flux extra energy out of the ice
fcondbot = fcondbot + einex/dt

ferr = abs( (enew-einit)/dt &
if (calc_Tsfc) then
! Flux extra energy out of the ice
fcondbot = fcondbot + einex_sfc_calc/dt
ferr = abs( (enew-einit)/dt &
- (fcondtopn - fcondbot + fswint) )
else
ferr = abs( (enew-einit+einex_sfc_flux)/dt &
- (fcondtopn - fcondbot + fswint) )
end if

! factor of 0.9 allows for roundoff errors later
if (ferr > 0.9_dbl_kind*ferrmax) then ! condition (5)
Expand Down Expand Up @@ -763,14 +811,6 @@ subroutine temperature_changes (dt, &
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, (zTsn(k),k=1,nslyr)
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'Matrix ice temperature diff:'
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, (dTmat(k),k=1,nilyr)
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'dqmat*hilyr/dt:'
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, (hilyr*dqmat(k)/dt,k=1,nilyr)
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'Final ice temperatures:'
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, (zTin(k),k=1,nilyr)
Expand Down
Loading
Loading