Skip to content

Flux coupling edits #8

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

Closed
kieranricardo wants to merge 30 commits into from
Closed

Flux coupling edits #8

Show file tree
Hide file tree
Changes from 23 commits
Commits
Show all changes
30 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
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 @@ -349,7 +349,7 @@ subroutine set_sfcflux (aicen, &

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

raicen = c1
raicen = c1 / aicen
Copy link
Copy Markdown
Collaborator

@anton-seaice anton-seaice Feb 18, 2026

Choose a reason for hiding this comment

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

So this is the same as line 360 now?

The comment says (When in standalone mode, fluxes ! are input as per ice area.) does that mean there is a standalone icepack using this subroutine, or the comment is out of date ?

Copy link
Copy Markdown
Author

@kieranricardo kieranricardo Feb 25, 2026 •
edited
Loading

Choose a reason for hiding this comment

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

I've changed this to use an #idfef(ACCESS3_CICE) (which is now defined in CICE)


#ifdef CICE_IN_NEMO
!----------------------------------------------------------------------
Expand Down
67 changes: 56 additions & 11 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, e_num)

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

real (kind=dbl_kind), intent(out):: &
e_num
! local variables

integer (kind=int_kind), parameter :: &
Expand Down Expand Up @@ -189,20 +191,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
e_num = c0
Top_T_was_reset_last_time = .false.
converged = .false.
l_snow = .false.
l_cold = .true.
Expand Down Expand Up @@ -428,6 +434,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 +443,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 +565,31 @@ 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

! Alex West: return this energy to the ocean

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

! Alex West: 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.
e_num = e_num + hslyr * dqmat_sn
else
Top_T_was_reset_last_time = .true.
endif
endif

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

endif

!-----------------------------------------------------------------
! If condition 1 or 2 failed, average new snow layer
Expand Down Expand Up @@ -592,6 +623,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. 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.
e_num = e_num + 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 @@ -683,11 +724,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/dt
ferr = abs( (enew-einit)/dt &
- (fcondtopn - fcondbot + fswint) )
else
ferr = abs( (enew-einit+e_num)/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
83 changes: 73 additions & 10 deletions columnphysics/icepack_therm_vertical.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,49 @@ module icepack_therm_vertical
contains

!=======================================================================

function cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) result(fcondtopn_solve)
Copy link
Copy Markdown
Collaborator

@anton-seaice anton-seaice Feb 18, 2026

Choose a reason for hiding this comment

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

should this be a subroutine? is there a particular reason to have a function ?

Copy link
Copy Markdown
Author

@kieranricardo kieranricardo Feb 24, 2026

Choose a reason for hiding this comment

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

this just returns a single float, so a function felt cleaner to me here.


Copy link
Copy Markdown
Collaborator

@anton-seaice anton-seaice Feb 18, 2026

Choose a reason for hiding this comment

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

This will need a description of some form

integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
real (kind=dbl_kind), intent(in) :: fcondtopn
real (kind=dbl_kind), intent(in) :: hin
real (kind=dbl_kind), intent(in) :: zTin(nilyr)
real (kind=dbl_kind), intent(in) :: zTsn(nslyr)
real (kind=dbl_kind), intent(in) :: hslyr

real (kind=dbl_kind) :: fcondtopn_solve

real (kind=dbl_kind), parameter :: ratio_Wm2_m = 1000.0, cold_temp_flag = c0 - 60.0
Copy link
Copy Markdown
Collaborator

@anton-seaice anton-seaice Feb 18, 2026

Choose a reason for hiding this comment

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

I guess we should put these in icepack_parameter.F90 , and add them to icepack_init_parameters


! AEW: New variables for cold-ice flux capping
real (kind=dbl_kind) :: top_layer_temp, &
reduce_ratio, &
reduce_amount


if (abs(fcondtopn) > ratio_Wm2_m * hin) then
fcondtopn_solve = sign(ratio_Wm2_m * hin,fcondtopn)

else
fcondtopn_solve = fcondtopn
endif

if (hslyr>hs_min) then
top_layer_temp = zTsn(1)
else
top_layer_temp = zTin(1)
endif

if ((top_layer_temp < cold_temp_flag) .and. (fcondtopn_solve < c0)) then
reduce_ratio = (cold_temp_flag - top_layer_temp) / (100.0 + cold_temp_flag)
reduce_amount = reduce_ratio * fcondtopn_solve
fcondtopn_solve = fcondtopn_solve - reduce_amount
endif

end function cap_conductive_flux

!
! Driver for updating ice and snow internal temperatures and
! computing thermodynamic growth rates and atmospheric fluxes.
Expand Down Expand Up @@ -246,6 +289,9 @@ subroutine thermo_vertical (dt, aicen, &
real (kind=dbl_kind) :: &
fadvocn, saltvol, dfsalt ! advective heat flux to ocean

real (kind=dbl_kind) :: &
fcondtopn_solve, fcondtopn_extra, e_num

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

!-----------------------------------------------------------------
Expand All @@ -272,6 +318,8 @@ subroutine thermo_vertical (dt, aicen, &
meltsliq= c0
massice(:) = c0
massliq(:) = c0
e_num = c0
fcondtopn_extra = c0
if (tr_pond) then
dpnd_flush = c0
dpnd_expon = c0
Expand All @@ -284,6 +332,8 @@ subroutine thermo_vertical (dt, aicen, &
fcondtopn = c0
endif

fcondtopn_solve = fcondtopn

!-----------------------------------------------------------------
! Compute variables needed for vertical thermo calculation
!-----------------------------------------------------------------
Expand Down Expand Up @@ -339,6 +389,12 @@ subroutine thermo_vertical (dt, aicen, &
if (icepack_warnings_aborted(subname)) return

else ! ktherm
fcondtopn_solve = cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr)
fcondtopn_extra = fcondtopn - fcondtopn_solve
Copy link
Copy Markdown
Collaborator

@anton-seaice anton-seaice Feb 18, 2026

Choose a reason for hiding this comment

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

Does this need an if ( .not. calc_Tsfc ) ?

If CICE is internally calculating the conductive flux then there is hopefully no issue that would mean the flux needs capping ?

Copy link
Copy Markdown
Author

@kieranricardo kieranricardo Feb 24, 2026

Choose a reason for hiding this comment

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

good catch. this shouldn't be doing anything if calc_Tsfc=true bc fcondtop will be calculated in the following routine, but we should be needlessly capping a flux that potentially isn't even defined.


! if (calc_Tsfc) then
! fcondtopn = fcondtopn_solve
! end if

call temperature_changes(dt, &
rhoa, flw, &
Expand All @@ -353,8 +409,8 @@ subroutine thermo_vertical (dt, aicen, &
Tsf, Tbot, &
fsensn, flatn, &
flwoutn, fsurfn, &
fcondtopn, fcondbotn, &
einit )
fcondtopn_solve, fcondbotn, &
einit, e_num)
if (icepack_warnings_aborted(subname)) return

endif ! ktherm
Expand Down Expand Up @@ -411,7 +467,8 @@ subroutine thermo_vertical (dt, aicen, &
mlt_onset, frz_onset, &
zSin, sss, &
sst, &
dsnow, rsnw)
dsnow, rsnw, &
e_num, fcondtopn_extra )
if (icepack_warnings_aborted(subname)) return

!-----------------------------------------------------------------
Expand All @@ -425,7 +482,7 @@ subroutine thermo_vertical (dt, aicen, &
fsnow, einit, &
einter, efinal, &
fcondtopn, fcondbotn, &
fadvocn, fbot )
fadvocn, fbot, e_num, fcondtopn_extra )
if (icepack_warnings_aborted(subname)) return

!-----------------------------------------------------------------
Expand Down Expand Up @@ -1070,7 +1127,8 @@ subroutine thickness_changes (dt, yday, &
mlt_onset, frz_onset,&
zSin, sss, &
sst, &
dsnow, rsnw)
dsnow, rsnw, &
e_num, fcondtopn_extra)

real (kind=dbl_kind), intent(in) :: &
dt , & ! time step
Expand Down Expand Up @@ -1138,6 +1196,8 @@ subroutine thickness_changes (dt, yday, &
sst , & ! sea surface temperature (C)
sss ! ocean salinity (PSU)

real (kind=dbl_kind), intent(in) :: &
e_num, fcondtopn_extra
! local variables

integer (kind=int_kind) :: &
Expand Down Expand Up @@ -1271,7 +1331,7 @@ subroutine thickness_changes (dt, yday, &
wk1 = (fsurfn - fcondtopn) * dt
etop_mlt = max(wk1, c0) ! etop_mlt > 0

wk1 = (fcondbotn - fbot) * dt
wk1 = (fcondbotn - fbot + fcondtopn_extra) * dt
ebot_mlt = max(wk1, c0) ! ebot_mlt > 0
ebot_gro = min(wk1, c0) ! ebot_gro < 0

Expand Down Expand Up @@ -1561,8 +1621,7 @@ subroutine thickness_changes (dt, yday, &
! fhocn is the available ocean heat that is left after use by ice
!-----------------------------------------------------------------

fhocnn = fbot &
+ (esub + etop_mlt + ebot_mlt)/dt
fhocnn = fbot + (esub + etop_mlt + ebot_mlt + e_num)/dt

!-----------------------------------------------------------------
! Add new snowfall at top surface
Expand Down Expand Up @@ -1987,7 +2046,7 @@ subroutine conservation_check_vthermo(dt, &
einit, einter, &
efinal, &
fcondtopn,fcondbotn, &
fadvocn, fbot )
fadvocn, fbot, e_num, fcondtopn_extra)

real (kind=dbl_kind), intent(in) :: &
dt ! time step
Expand All @@ -2000,7 +2059,7 @@ subroutine conservation_check_vthermo(dt, &
fsnow , & ! snowfall rate (kg m-2 s-1)
fcondtopn , &
fadvocn , &
fbot
fbot, e_num, fcondtopn_extra

real (kind=dbl_kind), intent(in) :: &
einit , & ! initial energy of melting (J m-2)
Expand Down Expand Up @@ -2053,6 +2112,10 @@ subroutine conservation_check_vthermo(dt, &
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'Input energy =', einp
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'Numerical energy =', e_num
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'fcondtopn_extra energy =', fcondtopn_extra
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'fbot,fcondbot:'
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, fbot,fcondbotn
Expand Down
Loading