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

Draft
kieranricardo wants to merge 21 commits into
base: main
Choose a base branch
from
Draft

Flux coupling edits #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
21 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
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

#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)

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

! 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

! 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
9 changes: 5 additions & 4 deletions configuration/driver/icedrv_init.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ subroutine input_data

real (kind=dbl_kind) :: ustar_min, albicev, albicei, albsnowv, albsnowi, &
ahmax, R_ice, R_pnd, R_snw, dT_mlt, rsnw_mlt, ksno, hi_min, Tliquidus_max, &
mu_rdg, hs0, dpscale, rfracmin, rfracmax, pndaspect, hs1, hp1, &
hs_min, mu_rdg, hs0, dpscale, rfracmin, rfracmax, pndaspect, hs1, hp1, &
apnd_sl, tscale_pnd_drain, itd_area_min, itd_mass_min, &
a_rapid_mode, Rac_rapid_mode, aspect_rapid_mode, dSdt_slow_mode, &
phi_c_slow_mode, phi_i_mushy, kalg, emissivity, floediam, hfrazilmin, &
Expand Down Expand Up @@ -161,7 +161,7 @@ subroutine input_data
a_rapid_mode, Rac_rapid_mode, aspect_rapid_mode, &
dSdt_slow_mode, phi_c_slow_mode, phi_i_mushy, &
floediam, hfrazilmin, Tliquidus_max, hi_min, &
tscale_pnd_drain
tscale_pnd_drain, hs_min

namelist /dynamics_nml/ &
kstrength, krdg_partic, krdg_redist, mu_rdg, &
Expand Down Expand Up @@ -232,7 +232,7 @@ subroutine input_data
krdg_redist_out=krdg_redist, mu_rdg_out=mu_rdg, &
atmbndy_out=atmbndy, calc_strair_out=calc_strair, &
formdrag_out=formdrag, highfreq_out=highfreq, &
emissivity_out=emissivity, &
emissivity_out=emissivity, hs_min_out=hs_min, &
kitd_out=kitd, kcatbound_out=kcatbound, hs0_out=hs0, &
dpscale_out=dpscale, frzpnd_out=frzpnd, &
rfracmin_out=rfracmin, rfracmax_out=rfracmax, &
Expand Down Expand Up @@ -854,6 +854,7 @@ subroutine input_data
write(nu_diag,1010) ' l_mpond_fresh = ', l_mpond_fresh
write(nu_diag,1005) ' ustar_min = ', ustar_min
write(nu_diag,1005) ' hi_min = ', hi_min
write(nu_diag,1005) ' hs_min = ', hs_min
write(nu_diag,1030) ' fbot_xfer_type = ', trim(fbot_xfer_type)
write(nu_diag,1010) ' oceanmixed_ice = ', oceanmixed_ice
write(nu_diag,1030) ' congel_freeze = ', trim(congel_freeze)
Expand Down Expand Up @@ -1035,7 +1036,7 @@ subroutine input_data
krdg_redist_in=krdg_redist, mu_rdg_in=mu_rdg, &
atmbndy_in=atmbndy, calc_strair_in=calc_strair, &
formdrag_in=formdrag, highfreq_in=highfreq, &
emissivity_in=emissivity, &
emissivity_in=emissivity, hs_min_in=hs_min, &
kitd_in=kitd, kcatbound_in=kcatbound, hs0_in=hs0, &
dpscale_in=dpscale, frzpnd_in=frzpnd, &
rfracmin_in=rfracmin, rfracmax_in=rfracmax, &
Expand Down
Loading