From bc088acf63b5bbfa04745bb7bba80258749134ce Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 29 Apr 2024 14:33:19 +1000 Subject: [PATCH 01/44] logging, tweaked newton solver tolerances, re-add einex --- columnphysics/icepack_therm_bl99.F90 | 132 +++++++++++++++-------- columnphysics/icepack_therm_shared.F90 | 2 +- columnphysics/icepack_therm_vertical.F90 | 89 +++++++++++++-- 3 files changed, 168 insertions(+), 55 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index f9485e32f..8b798e59e 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -12,6 +12,7 @@ module icepack_therm_bl99 use icepack_kinds + use ESMF use icepack_parameters, only: c0, c1, c2, p1, p5, puny use icepack_parameters, only: rhoi, rhos, hs_min, cp_ice, cp_ocn, depressT, Lfresh, ksno, kice use icepack_parameters, only: conduct, calc_Tsfc, semi_implicit_Tsfc @@ -69,7 +70,7 @@ subroutine temperature_changes (dt, & fsensn, flatn, & flwoutn, fsurfn, & fcondtopn,fcondbot, & - einit ) + einit, e_num) real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -122,6 +123,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 :: & @@ -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 + e_num = c0 + Top_T_was_reset_last_time = .false. converged = .false. l_snow = .false. l_cold = .true. @@ -263,53 +270,56 @@ subroutine temperature_changes (dt, & ! has already computed fsurf. (Unless we adjust fsurf here) !----------------------------------------------------------------- !mclaren: Should there be an if calc_Tsfc statement here then?? + if (calc_Tsfc) then + if (sw_redist) then - if (sw_redist) then + if (solve_zsal) sw_dtemp = p1 ! lower tolerance with dynamic salinity - do k = 1, nilyr + do k = 1, nilyr - Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc - if (Tin_init(k) <= Tmlts(k) - sw_dtemp) then - if (l_brine) then - ci = cp_ice - Lfresh * Tmlts(k) / (Tin_init(k)**2) - Iswabs_tmp = min(Iswabs(k), & - sw_frac*(Tmlts(k)-Tin_init(k))*ci/dt_rhoi_hlyr) - else - ci = cp_ice - Iswabs_tmp = min(Iswabs(k), & - sw_frac*(-Tin_init(k))*ci/dt_rhoi_hlyr) + Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc + if (Tin_init(k) <= Tmlts(k) - sw_dtemp) then + if (l_brine) then + ci = cp_ice - Lfresh * Tmlts(k) / (Tin_init(k)**2) + Iswabs_tmp = min(Iswabs(k), & + sw_frac*(Tmlts(k)-Tin_init(k))*ci/dt_rhoi_hlyr) + else + ci = cp_ice + Iswabs_tmp = min(Iswabs(k), & + sw_frac*(-Tin_init(k))*ci/dt_rhoi_hlyr) + endif endif - endif - if (Iswabs_tmp < puny) Iswabs_tmp = c0 + if (Iswabs_tmp < puny) Iswabs_tmp = c0 - dswabs = min(Iswabs(k) - Iswabs_tmp, fswint) + dswabs = min(Iswabs(k) - Iswabs_tmp, fswint) - fswsfc = fswsfc + dswabs - fswint = fswint - dswabs - Iswabs(k) = Iswabs_tmp + fswsfc = fswsfc + dswabs + fswint = fswint - dswabs + Iswabs(k) = Iswabs_tmp - enddo + enddo - do k = 1, nslyr - if (l_snow) then + do k = 1, nslyr + if (l_snow) then - Sswabs_tmp = c0 - if (Tsn_init(k) <= -sw_dtemp) then - Sswabs_tmp = min(Sswabs(k), & - -sw_frac*Tsn_init(k)/etas(k)) - endif - if (Sswabs_tmp < puny) Sswabs_tmp = c0 + Sswabs_tmp = c0 + if (Tsn_init(k) <= -sw_dtemp) then + Sswabs_tmp = min(Sswabs(k), & + -sw_frac*Tsn_init(k)/etas(k)) + endif + if (Sswabs_tmp < puny) Sswabs_tmp = c0 - dswabs = min(Sswabs(k) - Sswabs_tmp, fswint) + dswabs = min(Sswabs(k) - Sswabs_tmp, fswint) - fswsfc = fswsfc + dswabs - fswint = fswint - dswabs - Sswabs(k) = Sswabs_tmp + fswsfc = fswsfc + dswabs + fswint = fswint - dswabs + Sswabs(k) = Sswabs_tmp - endif - enddo + endif + enddo - endif + endif + endif ! calc_Tsfc if (semi_implicit_Tsfc) then fsurfn = fsurfn + fswsfc ! this is the total heat flux @@ -428,6 +438,7 @@ subroutine temperature_changes (dt, & else + fcondtopn_force = fcondtopn - fcondtopn_reduction call get_matrix_elements_know_Tsfc ( & l_snow, Tbot, & Tin_init, Tsn_init, & @@ -436,7 +447,7 @@ subroutine temperature_changes (dt, & etai, etas, & sbdiag, diag, & spdiag, rhs, & - fcondtopn) + fcondtopn_force) if (icepack_warnings_aborted(subname)) return endif ! calc_Tsfc @@ -558,7 +569,32 @@ subroutine temperature_changes (dt, & else zTsn(k) = c0 endif - if (l_brine) zTsn(k) = min(zTsn(k), c0) + ! 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 @@ -592,6 +628,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)) @@ -686,11 +732,11 @@ subroutine temperature_changes (dt, & ! Flux extra energy out of the ice fcondbot = fcondbot + einex/dt - ferr = abs( (enew-einit)/dt & + ferr = abs( (enew-einit+e_num)/dt & - (fcondtopn - fcondbot + fswint) ) ! factor of 0.9 allows for roundoff errors later - if (ferr > 0.9_dbl_kind*ferrmax) then ! condition (5) + if (ferr > 0.9_dbl_kind*ferrmax*0.05_dbl_kind) then ! condition (5) converged = .false. diff --git a/columnphysics/icepack_therm_shared.F90 b/columnphysics/icepack_therm_shared.F90 index 73d560b60..8c1116fe7 100644 --- a/columnphysics/icepack_therm_shared.F90 +++ b/columnphysics/icepack_therm_shared.F90 @@ -39,7 +39,7 @@ module icepack_therm_shared adjust_enthalpy real (kind=dbl_kind), parameter, public :: & - ferrmax = 1.0e-3_dbl_kind ! max allowed energy flux error (W m-2) + ferrmax = 2.0e-2_dbl_kind ! max allowed energy flux error (W m-2) ! recommend ferrmax < 0.01 W m-2 real (kind=dbl_kind), parameter, public :: & diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 00b3756d6..758ed0d30 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -65,6 +65,8 @@ module icepack_therm_vertical use icepack_meltpond_sealvl, only: compute_ponds_sealvl use icepack_snow, only: drain_snow + use ESMF + implicit none private @@ -75,6 +77,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. @@ -246,6 +291,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)' !----------------------------------------------------------------- @@ -272,6 +320,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 @@ -284,6 +334,8 @@ subroutine thermo_vertical (dt, aicen, & fcondtopn = c0 endif + fcondtopn_solve = fcondtopn + !----------------------------------------------------------------- ! Compute variables needed for vertical thermo calculation !----------------------------------------------------------------- @@ -339,6 +391,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, & @@ -353,8 +411,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 @@ -400,7 +458,7 @@ subroutine thermo_vertical (dt, aicen, & smliq, massliq, & fbot, Tbot, & flatn, fsurfn, & - fcondtopn, fcondbotn, & + fcondtopn_solve, fcondbotn, & fsnow, hsn_new, & fhocnn, evapn, & evapsn, evapin, & @@ -411,7 +469,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 !----------------------------------------------------------------- @@ -425,7 +484,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 !----------------------------------------------------------------- @@ -1070,7 +1129,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 @@ -1138,6 +1198,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) :: & @@ -1271,7 +1333,8 @@ 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 + ! wk1 = (fcondbotn - fbot) * dt ebot_mlt = max(wk1, c0) ! ebot_mlt > 0 ebot_gro = min(wk1, c0) ! ebot_gro < 0 @@ -1561,8 +1624,8 @@ 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 + ! fhocnn = fbot + (esub + etop_mlt + ebot_mlt)/dt !----------------------------------------------------------------- ! Add new snowfall at top surface @@ -1987,7 +2050,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 @@ -2000,7 +2063,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) @@ -2053,6 +2116,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 From ea4cde3c6324abedaba625e56afa7882fb069650 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 6 May 2024 12:21:12 +1000 Subject: [PATCH 02/44] fix argument ordering in set_sfcflux --- columnphysics/icepack_therm_bl99.F90 | 2 +- columnphysics/icepack_therm_vertical.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 8b798e59e..e93da4034 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -730,7 +730,7 @@ subroutine temperature_changes (dt, & (zTin(nilyr) - Tbot) ! Flux extra energy out of the ice - fcondbot = fcondbot + einex/dt + ! fcondbot = fcondbot + einex/dt ferr = abs( (enew-einit+e_num)/dt & - (fcondtopn - fcondbot + fswint) ) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 758ed0d30..2bd9b9bb1 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -1331,6 +1331,7 @@ subroutine thickness_changes (dt, yday, & econ = min(wk1, c0) ! energy for condensation, < 0 wk1 = (fsurfn - fcondtopn) * dt + ! wk1 = fsurfn * dt etop_mlt = max(wk1, c0) ! etop_mlt > 0 wk1 = (fcondbotn - fbot + fcondtopn_extra) * dt @@ -2900,7 +2901,6 @@ subroutine icepack_step_therm1(dt, & !----------------------------------------------------------------- ! Vertical thermodynamics: Heat conduction, growth and melting. !----------------------------------------------------------------- - if (.not.(calc_Tsfc)) then ! If not calculating surface temperature and fluxes, set From d4b4a31fb525d6da5d7b1ca7c580a6052443691f Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 7 May 2024 16:16:53 +1000 Subject: [PATCH 03/44] remove cm2 style conductive flux limiting --- columnphysics/icepack_therm_bl99.F90 | 6 +++--- columnphysics/icepack_therm_vertical.F90 | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index e93da4034..f2eab66f5 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -586,12 +586,12 @@ subroutine temperature_changes (dt, & 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 + e_num = e_num + hslyr * dqmat_sn else Top_T_was_reset_last_time = .true. endif endif - + zTsn(k) = min(zTsn(k), c0) endif @@ -730,7 +730,7 @@ subroutine temperature_changes (dt, & (zTin(nilyr) - Tbot) ! Flux extra energy out of the ice - ! fcondbot = fcondbot + einex/dt + fcondbot = fcondbot + einex/dt ferr = abs( (enew-einit+e_num)/dt & - (fcondtopn - fcondbot + fswint) ) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 2bd9b9bb1..80c04f71f 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -391,8 +391,8 @@ 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 + ! fcondtopn_solve = cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) + ! fcondtopn_extra = fcondtopn - fcondtopn_solve ! if (calc_Tsfc) then ! fcondtopn = fcondtopn_solve @@ -411,7 +411,7 @@ subroutine thermo_vertical (dt, aicen, & Tsf, Tbot, & fsensn, flatn, & flwoutn, fsurfn, & - fcondtopn_solve, fcondbotn, & + fcondtopn, fcondbotn, & einit, e_num) if (icepack_warnings_aborted(subname)) return @@ -458,7 +458,7 @@ subroutine thermo_vertical (dt, aicen, & smliq, massliq, & fbot, Tbot, & flatn, fsurfn, & - fcondtopn_solve, fcondbotn, & + fcondtopn, fcondbotn, & fsnow, hsn_new, & fhocnn, evapn, & evapsn, evapin, & From a15fc03b538c8ba50ea16c4595b398a131dd0358 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 20 May 2024 11:46:28 +1000 Subject: [PATCH 04/44] use cm2 style conductive flux limiting --- columnphysics/icepack_therm_bl99.F90 | 2 +- columnphysics/icepack_therm_vertical.F90 | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index f2eab66f5..e40c9dd16 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -730,7 +730,7 @@ subroutine temperature_changes (dt, & (zTin(nilyr) - Tbot) ! Flux extra energy out of the ice - fcondbot = fcondbot + einex/dt + ! fcondbot = fcondbot + einex/dt ferr = abs( (enew-einit+e_num)/dt & - (fcondtopn - fcondbot + fswint) ) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 80c04f71f..d0f064575 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -391,8 +391,8 @@ 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 + fcondtopn_solve = cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) + fcondtopn_extra = fcondtopn - fcondtopn_solve ! if (calc_Tsfc) then ! fcondtopn = fcondtopn_solve @@ -411,7 +411,7 @@ subroutine thermo_vertical (dt, aicen, & Tsf, Tbot, & fsensn, flatn, & flwoutn, fsurfn, & - fcondtopn, fcondbotn, & + fcondtopn_solve, fcondbotn, & einit, e_num) if (icepack_warnings_aborted(subname)) return From 65259a0e0acdd6afe707c98936dc462d5db95dc6 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 22 Jul 2024 13:56:26 +1000 Subject: [PATCH 05/44] change flux energy error threshold, and reduce fsurf before thermo solve --- columnphysics/icepack_therm_bl99.F90 | 2 +- columnphysics/icepack_therm_shared.F90 | 2 +- columnphysics/icepack_therm_vertical.F90 | 5 +++-- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index e40c9dd16..008f66a9c 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -736,7 +736,7 @@ subroutine temperature_changes (dt, & - (fcondtopn - fcondbot + fswint) ) ! factor of 0.9 allows for roundoff errors later - if (ferr > 0.9_dbl_kind*ferrmax*0.05_dbl_kind) then ! condition (5) + if (ferr > 0.9_dbl_kind*ferrmax) then ! condition (5) converged = .false. diff --git a/columnphysics/icepack_therm_shared.F90 b/columnphysics/icepack_therm_shared.F90 index 8c1116fe7..73d560b60 100644 --- a/columnphysics/icepack_therm_shared.F90 +++ b/columnphysics/icepack_therm_shared.F90 @@ -39,7 +39,7 @@ module icepack_therm_shared adjust_enthalpy real (kind=dbl_kind), parameter, public :: & - ferrmax = 2.0e-2_dbl_kind ! max allowed energy flux error (W m-2) + ferrmax = 1.0e-3_dbl_kind ! max allowed energy flux error (W m-2) ! recommend ferrmax < 0.01 W m-2 real (kind=dbl_kind), parameter, public :: & diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index d0f064575..e9ce400f5 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -292,7 +292,7 @@ subroutine thermo_vertical (dt, aicen, & fadvocn, saltvol, dfsalt ! advective heat flux to ocean real (kind=dbl_kind) :: & - fcondtopn_solve, fcondtopn_extra, e_num + fcondtopn_solve, fcondtopn_extra, e_num, fsurfn_solve character(len=*),parameter :: subname='(thermo_vertical)' @@ -393,6 +393,7 @@ subroutine thermo_vertical (dt, aicen, & else ! ktherm fcondtopn_solve = cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) fcondtopn_extra = fcondtopn - fcondtopn_solve + fsurfn_solve = fsurfn - fcondtopn_extra ! if (calc_Tsfc) then ! fcondtopn = fcondtopn_solve @@ -410,7 +411,7 @@ subroutine thermo_vertical (dt, aicen, & zSin, & Tsf, Tbot, & fsensn, flatn, & - flwoutn, fsurfn, & + flwoutn, fsurfn_solve, & fcondtopn_solve, fcondbotn, & einit, e_num) if (icepack_warnings_aborted(subname)) return From e3b9746d93ff33d2d3c65f8f84120a9abe058eb1 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 23 Jul 2024 10:10:23 +1000 Subject: [PATCH 06/44] don't reduce fsurf before thermo solve - has no effect --- columnphysics/icepack_therm_vertical.F90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index e9ce400f5..d0f064575 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -292,7 +292,7 @@ subroutine thermo_vertical (dt, aicen, & fadvocn, saltvol, dfsalt ! advective heat flux to ocean real (kind=dbl_kind) :: & - fcondtopn_solve, fcondtopn_extra, e_num, fsurfn_solve + fcondtopn_solve, fcondtopn_extra, e_num character(len=*),parameter :: subname='(thermo_vertical)' @@ -393,7 +393,6 @@ subroutine thermo_vertical (dt, aicen, & else ! ktherm fcondtopn_solve = cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) fcondtopn_extra = fcondtopn - fcondtopn_solve - fsurfn_solve = fsurfn - fcondtopn_extra ! if (calc_Tsfc) then ! fcondtopn = fcondtopn_solve @@ -411,7 +410,7 @@ subroutine thermo_vertical (dt, aicen, & zSin, & Tsf, Tbot, & fsensn, flatn, & - flwoutn, fsurfn_solve, & + flwoutn, fsurfn, & fcondtopn_solve, fcondbotn, & einit, e_num) if (icepack_warnings_aborted(subname)) return From 2cf33b149697245d48b1ad0ab9b1963465ddd1b1 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 20 Aug 2024 14:38:39 +1000 Subject: [PATCH 07/44] convert GBM ice fluxes to local --- columnphysics/icepack_flux.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/columnphysics/icepack_flux.F90 b/columnphysics/icepack_flux.F90 index d55780763..c441a7cc3 100644 --- a/columnphysics/icepack_flux.F90 +++ b/columnphysics/icepack_flux.F90 @@ -349,7 +349,7 @@ subroutine set_sfcflux (aicen, & character(len=*),parameter :: subname='(set_sfcflux)' - raicen = c1 + raicen = c1 / aicen #ifdef CICE_IN_NEMO !---------------------------------------------------------------------- From af43525cd7c92d7efc5471a66df3c5d3ba1eebbd Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 26 Aug 2024 14:24:18 +1000 Subject: [PATCH 08/44] remove scaling by inverse ice fraction --- columnphysics/icepack_flux.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/columnphysics/icepack_flux.F90 b/columnphysics/icepack_flux.F90 index c441a7cc3..bca86fc10 100644 --- a/columnphysics/icepack_flux.F90 +++ b/columnphysics/icepack_flux.F90 @@ -349,7 +349,7 @@ subroutine set_sfcflux (aicen, & character(len=*),parameter :: subname='(set_sfcflux)' - raicen = c1 / aicen + raicen = c1 !/ aicen #ifdef CICE_IN_NEMO !---------------------------------------------------------------------- From d1ff164290f7ba7d5a8880e7cc08190a33fd9bb4 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 17 Sep 2024 11:54:44 +1000 Subject: [PATCH 09/44] increase ice and snow thickness minimum values to match CM2, revert conductive flux limits to CM2 --- columnphysics/icepack_itd.F90 | 2 +- columnphysics/icepack_parameters.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index 8802a943b..df17671b9 100644 --- a/columnphysics/icepack_itd.F90 +++ b/columnphysics/icepack_itd.F90 @@ -26,7 +26,7 @@ module icepack_itd use icepack_kinds - use icepack_parameters, only: c0, c1, c2, c3, c15, c25, c100, p1, p01, p001, p5, puny + use icepack_parameters, only: c0, c1, c2, c3, c15, c25, c100, p1, p01, p001, p5, puny, p2 use icepack_parameters, only: Lfresh, rhos, ice_ref_salinity, hs_min, cp_ice, rhoi use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin, rsnw_fall, rhosnew use icepack_tracers, only: ncat, nilyr, nslyr, nblyr, ntrcr, nbtrcr, n_aero diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index fe6ce051e..69e61116c 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -120,7 +120,7 @@ module icepack_parameters kice = 2.03_dbl_kind ,&! thermal conductivity of fresh ice(W/m/deg) ! kice is only used with ktherm=1 (BL99) and conduct='MU71' ksno = 0.30_dbl_kind ,&! thermal conductivity of snow (W/m/deg) - hs_min = 1.e-4_dbl_kind ,&! min snow thickness for computing zTsn (m) + hs_min = 0.10_dbl_kind ,&! min snow thickness for computing zTsn (m) snowpatch = 0.02_dbl_kind ,&! parameter for fractional snow area (m) saltmax = 3.2_dbl_kind ,&! max salinity at ice base for BL99 (ppt) ! phi_init, dSin0_frazil are for mushy thermo From e6e429c4a29e3abcb3f9f9ef752a17a02fa3ddc7 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 17 Sep 2024 12:11:48 +1000 Subject: [PATCH 10/44] scale ice fluxes by changing ice area --- columnphysics/icepack_flux.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/columnphysics/icepack_flux.F90 b/columnphysics/icepack_flux.F90 index bca86fc10..c441a7cc3 100644 --- a/columnphysics/icepack_flux.F90 +++ b/columnphysics/icepack_flux.F90 @@ -349,7 +349,7 @@ subroutine set_sfcflux (aicen, & character(len=*),parameter :: subname='(set_sfcflux)' - raicen = c1 !/ aicen + raicen = c1 / aicen #ifdef CICE_IN_NEMO !---------------------------------------------------------------------- From 07352ba4b07385e9da79def88126697f83c20ced Mon Sep 17 00:00:00 2001 From: Spencer Wong Date: Mon, 31 Mar 2025 22:15:46 +1100 Subject: [PATCH 11/44] Remove stray deprecated zsal flag --- columnphysics/icepack_therm_bl99.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 008f66a9c..9ac9ed314 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -273,8 +273,6 @@ subroutine temperature_changes (dt, & if (calc_Tsfc) then if (sw_redist) then - if (solve_zsal) sw_dtemp = p1 ! lower tolerance with dynamic salinity - do k = 1, nilyr Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc From e8ef69f123ed829d45cb12e863498da1dba9eb7d Mon Sep 17 00:00:00 2001 From: Spencer Wong Date: Wed, 9 Apr 2025 15:11:37 +1000 Subject: [PATCH 12/44] Change default hi_min to p1 --- columnphysics/icepack_parameters.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index 69e61116c..9a9ea8180 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -129,7 +129,7 @@ module icepack_parameters Tliquidus_max = c0 ,&! maximum liquidus temperature of mush (C) dSin0_frazil = c3 ,&! bulk salinity reduction of newly formed frazil ustar_min = 0.005_dbl_kind ,&! minimum friction velocity for ocean heat flux (m/s) - hi_min = p01 ,&! minimum ice thickness allowed (m) for thermo + hi_min = p1 ,&! minimum ice thickness allowed (m) for thermo ! mushy thermo a_rapid_mode = 0.5e-3_dbl_kind,&! channel radius for rapid drainage mode (m) Rac_rapid_mode = 10.0_dbl_kind,&! critical Rayleigh number From 9bd1b3f24cc94f06423baba3852fc0d2e4558387 Mon Sep 17 00:00:00 2001 From: Spencer Wong <88933912+blimlim@users.noreply.github.com> Date: Thu, 24 Apr 2025 09:43:04 +1000 Subject: [PATCH 13/44] Remove unnecessary comments Co-authored-by: Kieran Ricardo --- columnphysics/icepack_therm_bl99.F90 | 1 - columnphysics/icepack_therm_vertical.F90 | 3 --- 2 files changed, 4 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 9ac9ed314..55d35cae8 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -567,7 +567,6 @@ 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 diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index d0f064575..a022dc51b 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -1331,11 +1331,9 @@ subroutine thickness_changes (dt, yday, & econ = min(wk1, c0) ! energy for condensation, < 0 wk1 = (fsurfn - fcondtopn) * dt - ! wk1 = fsurfn * dt etop_mlt = max(wk1, c0) ! etop_mlt > 0 wk1 = (fcondbotn - fbot + fcondtopn_extra) * dt - ! wk1 = (fcondbotn - fbot) * dt ebot_mlt = max(wk1, c0) ! ebot_mlt > 0 ebot_gro = min(wk1, c0) ! ebot_gro < 0 @@ -1626,7 +1624,6 @@ subroutine thickness_changes (dt, yday, & !----------------------------------------------------------------- fhocnn = fbot + (esub + etop_mlt + ebot_mlt + e_num)/dt - ! fhocnn = fbot + (esub + etop_mlt + ebot_mlt)/dt !----------------------------------------------------------------- ! Add new snowfall at top surface From 140e8cb68aa1f0280fa394529fb7104aecc5a827 Mon Sep 17 00:00:00 2001 From: Spencer Wong Date: Thu, 24 Apr 2025 10:49:21 +1000 Subject: [PATCH 14/44] Review suggestions: remove unnecessary include and comment --- columnphysics/icepack_itd.F90 | 2 +- columnphysics/icepack_therm_bl99.F90 | 3 --- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index df17671b9..8802a943b 100644 --- a/columnphysics/icepack_itd.F90 +++ b/columnphysics/icepack_itd.F90 @@ -26,7 +26,7 @@ module icepack_itd use icepack_kinds - use icepack_parameters, only: c0, c1, c2, c3, c15, c25, c100, p1, p01, p001, p5, puny, p2 + use icepack_parameters, only: c0, c1, c2, c3, c15, c25, c100, p1, p01, p001, p5, puny use icepack_parameters, only: Lfresh, rhos, ice_ref_salinity, hs_min, cp_ice, rhoi use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin, rsnw_fall, rhosnew use icepack_tracers, only: ncat, nilyr, nslyr, nblyr, ntrcr, nbtrcr, n_aero diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 55d35cae8..825bbcf3a 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -726,9 +726,6 @@ 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+e_num)/dt & - (fcondtopn - fcondbot + fswint) ) From bf0995b4cc275203589227f6f4ebef0be3de4f86 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 3 Feb 2026 15:10:13 +1100 Subject: [PATCH 15/44] white space changes, remove potentiall unescessary calc_Tsfc statement --- columnphysics/icepack_therm_bl99.F90 | 63 ++++++++---------------- columnphysics/icepack_therm_vertical.F90 | 3 +- 2 files changed, 21 insertions(+), 45 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 825bbcf3a..68c4e2940 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -12,7 +12,6 @@ module icepack_therm_bl99 use icepack_kinds - use ESMF use icepack_parameters, only: c0, c1, c2, p1, p5, puny use icepack_parameters, only: rhoi, rhos, hs_min, cp_ice, cp_ocn, depressT, Lfresh, ksno, kice use icepack_parameters, only: conduct, calc_Tsfc, semi_implicit_Tsfc @@ -270,54 +269,32 @@ subroutine temperature_changes (dt, & ! has already computed fsurf. (Unless we adjust fsurf here) !----------------------------------------------------------------- !mclaren: Should there be an if calc_Tsfc statement here then?? - if (calc_Tsfc) then - if (sw_redist) then - - do k = 1, nilyr - - Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc - if (Tin_init(k) <= Tmlts(k) - sw_dtemp) then - if (l_brine) then - ci = cp_ice - Lfresh * Tmlts(k) / (Tin_init(k)**2) - Iswabs_tmp = min(Iswabs(k), & - sw_frac*(Tmlts(k)-Tin_init(k))*ci/dt_rhoi_hlyr) - else - ci = cp_ice - Iswabs_tmp = min(Iswabs(k), & - sw_frac*(-Tin_init(k))*ci/dt_rhoi_hlyr) - endif - endif - if (Iswabs_tmp < puny) Iswabs_tmp = c0 - - dswabs = min(Iswabs(k) - Iswabs_tmp, fswint) - - fswsfc = fswsfc + dswabs - fswint = fswint - dswabs - Iswabs(k) = Iswabs_tmp - enddo + if (sw_redist) then - do k = 1, nslyr - if (l_snow) then - - Sswabs_tmp = c0 - if (Tsn_init(k) <= -sw_dtemp) then - Sswabs_tmp = min(Sswabs(k), & - -sw_frac*Tsn_init(k)/etas(k)) - endif - if (Sswabs_tmp < puny) Sswabs_tmp = c0 + do k = 1, nilyr - dswabs = min(Sswabs(k) - Sswabs_tmp, fswint) + Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc + if (Tin_init(k) <= Tmlts(k) - sw_dtemp) then + if (l_brine) then + ci = cp_ice - Lfresh * Tmlts(k) / (Tin_init(k)**2) + Iswabs_tmp = min(Iswabs(k), & + sw_frac*(Tmlts(k)-Tin_init(k))*ci/dt_rhoi_hlyr) + else + ci = cp_ice + Iswabs_tmp = min(Iswabs(k), & + sw_frac*(-Tin_init(k))*ci/dt_rhoi_hlyr) + endif + endif + if (Iswabs_tmp < puny) Iswabs_tmp = c0 - fswsfc = fswsfc + dswabs - fswint = fswint - dswabs - Sswabs(k) = Sswabs_tmp + dswabs = min(Iswabs(k) - Iswabs_tmp, fswint) - endif - enddo + fswsfc = fswsfc + dswabs + fswint = fswint - dswabs + Iswabs(k) = Iswabs_tmp - endif - endif ! calc_Tsfc + enddo if (semi_implicit_Tsfc) then fsurfn = fsurfn + fswsfc ! this is the total heat flux diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index a022dc51b..bb6287fba 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -65,8 +65,6 @@ module icepack_therm_vertical use icepack_meltpond_sealvl, only: compute_ponds_sealvl use icepack_snow, only: drain_snow - use ESMF - implicit none private @@ -2898,6 +2896,7 @@ subroutine icepack_step_therm1(dt, & !----------------------------------------------------------------- ! Vertical thermodynamics: Heat conduction, growth and melting. !----------------------------------------------------------------- + if (.not.(calc_Tsfc)) then ! If not calculating surface temperature and fluxes, set From 5d2160dac883f7e0ed5af6edb3c9bc2018ff344d Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 3 Feb 2026 15:11:18 +1100 Subject: [PATCH 16/44] fix edit --- columnphysics/icepack_therm_bl99.F90 | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 68c4e2940..3ec445911 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -296,6 +296,27 @@ subroutine temperature_changes (dt, & enddo + do k = 1, nslyr + if (l_snow) then + + Sswabs_tmp = c0 + if (Tsn_init(k) <= -sw_dtemp) then + Sswabs_tmp = min(Sswabs(k), & + -sw_frac*Tsn_init(k)/etas(k)) + endif + if (Sswabs_tmp < puny) Sswabs_tmp = c0 + + dswabs = min(Sswabs(k) - Sswabs_tmp, fswint) + + fswsfc = fswsfc + dswabs + fswint = fswint - dswabs + Sswabs(k) = Sswabs_tmp + + endif + enddo + + endif + if (semi_implicit_Tsfc) then fsurfn = fsurfn + fswsfc ! this is the total heat flux endif From 6782499eaaa73e5f430d7ad9ec6339d801c49e3e Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 3 Feb 2026 15:13:49 +1100 Subject: [PATCH 17/44] calculate energy error differently for calc_tsfc options --- columnphysics/icepack_therm_bl99.F90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 3ec445911..d2212282f 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -720,12 +720,18 @@ subroutine temperature_changes (dt, & ! Condition 5: check for energy conservation error ! Change in internal ice energy should equal net energy input. !----------------------------------------------------------------- - fcondbot = kh(1+nslyr+nilyr) * & (zTin(nilyr) - Tbot) - ferr = abs( (enew-einit+e_num)/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) From ed20c15e6aae55c440bc3b2e15591a48e50812ad Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 3 Feb 2026 15:14:17 +1100 Subject: [PATCH 18/44] white space change --- columnphysics/icepack_therm_bl99.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index d2212282f..5809c06e0 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -720,6 +720,7 @@ subroutine temperature_changes (dt, & ! Condition 5: check for energy conservation error ! Change in internal ice energy should equal net energy input. !----------------------------------------------------------------- + fcondbot = kh(1+nslyr+nilyr) * & (zTin(nilyr) - Tbot) From 7036d9cf792c793b6589e943af3db2884b3f0797 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 3 Feb 2026 15:19:57 +1100 Subject: [PATCH 19/44] revert snow and ice thickness limits to default values --- columnphysics/icepack_parameters.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index 9a9ea8180..301d05bba 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -120,7 +120,7 @@ module icepack_parameters kice = 2.03_dbl_kind ,&! thermal conductivity of fresh ice(W/m/deg) ! kice is only used with ktherm=1 (BL99) and conduct='MU71' ksno = 0.30_dbl_kind ,&! thermal conductivity of snow (W/m/deg) - hs_min = 0.10_dbl_kind ,&! min snow thickness for computing zTsn (m) + hs_min = 1.e-4_dbl_kind ,&! min snow thickness for computing zTsn (m) snowpatch = 0.02_dbl_kind ,&! parameter for fractional snow area (m) saltmax = 3.2_dbl_kind ,&! max salinity at ice base for BL99 (ppt) ! phi_init, dSin0_frazil are for mushy thermo @@ -129,7 +129,7 @@ module icepack_parameters Tliquidus_max = c0 ,&! maximum liquidus temperature of mush (C) dSin0_frazil = c3 ,&! bulk salinity reduction of newly formed frazil ustar_min = 0.005_dbl_kind ,&! minimum friction velocity for ocean heat flux (m/s) - hi_min = p1 ,&! minimum ice thickness allowed (m) for thermo + hi_min = p01 ,&! minimum ice thickness allowed (m) for thermo ! mushy thermo a_rapid_mode = 0.5e-3_dbl_kind,&! channel radius for rapid drainage mode (m) Rac_rapid_mode = 10.0_dbl_kind,&! critical Rayleigh number From 10faac69c2bc3babc0e392d551d4fbb8afa4060c Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 3 Feb 2026 15:20:20 +1100 Subject: [PATCH 20/44] white space change --- columnphysics/icepack_parameters.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index 301d05bba..fe6ce051e 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -120,7 +120,7 @@ module icepack_parameters kice = 2.03_dbl_kind ,&! thermal conductivity of fresh ice(W/m/deg) ! kice is only used with ktherm=1 (BL99) and conduct='MU71' ksno = 0.30_dbl_kind ,&! thermal conductivity of snow (W/m/deg) - hs_min = 1.e-4_dbl_kind ,&! min snow thickness for computing zTsn (m) + hs_min = 1.e-4_dbl_kind ,&! min snow thickness for computing zTsn (m) snowpatch = 0.02_dbl_kind ,&! parameter for fractional snow area (m) saltmax = 3.2_dbl_kind ,&! max salinity at ice base for BL99 (ppt) ! phi_init, dSin0_frazil are for mushy thermo @@ -129,7 +129,7 @@ module icepack_parameters Tliquidus_max = c0 ,&! maximum liquidus temperature of mush (C) dSin0_frazil = c3 ,&! bulk salinity reduction of newly formed frazil ustar_min = 0.005_dbl_kind ,&! minimum friction velocity for ocean heat flux (m/s) - hi_min = p01 ,&! minimum ice thickness allowed (m) for thermo + hi_min = p01 ,&! minimum ice thickness allowed (m) for thermo ! mushy thermo a_rapid_mode = 0.5e-3_dbl_kind,&! channel radius for rapid drainage mode (m) Rac_rapid_mode = 10.0_dbl_kind,&! critical Rayleigh number From 789f98280fc80b18527be3e04f61761d764c7215 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 3 Feb 2026 15:45:20 +1100 Subject: [PATCH 21/44] make hs_min configurable --- configuration/driver/icedrv_init.F90 | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90 index 133683dd6..e5ab9800b 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/configuration/driver/icedrv_init.F90 @@ -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, & @@ -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, & @@ -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, & @@ -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) @@ -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, & From d2c2a0926973e9adb1aa0c194d8eb8a16aa0096b Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 17 Feb 2026 13:38:45 +1100 Subject: [PATCH 22/44] revert configurable hs_min --- configuration/driver/icedrv_init.F90 | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90 index e5ab9800b..133683dd6 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/configuration/driver/icedrv_init.F90 @@ -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, & - hs_min, mu_rdg, hs0, dpscale, rfracmin, rfracmax, pndaspect, hs1, hp1, & + 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, & @@ -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, hs_min + tscale_pnd_drain namelist /dynamics_nml/ & kstrength, krdg_partic, krdg_redist, mu_rdg, & @@ -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, hs_min_out=hs_min, & + emissivity_out=emissivity, & kitd_out=kitd, kcatbound_out=kcatbound, hs0_out=hs0, & dpscale_out=dpscale, frzpnd_out=frzpnd, & rfracmin_out=rfracmin, rfracmax_out=rfracmax, & @@ -854,7 +854,6 @@ 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) @@ -1036,7 +1035,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, hs_min_in=hs_min, & + emissivity_in=emissivity, & kitd_in=kitd, kcatbound_in=kcatbound, hs0_in=hs0, & dpscale_in=dpscale, frzpnd_in=frzpnd, & rfracmin_in=rfracmin, rfracmax_in=rfracmax, & From 1fb991d66eee4c67e62ee3662fed91312971fa0a Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 17 Feb 2026 16:13:52 +1100 Subject: [PATCH 23/44] remove tab --- columnphysics/icepack_therm_bl99.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 5809c06e0..f72e6af1f 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -581,7 +581,7 @@ subroutine temperature_changes (dt, & 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 + e_num = e_num + hslyr * dqmat_sn else Top_T_was_reset_last_time = .true. endif From 8516838c3bd332d34a9202a7c0e7105d1d7f3901 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 24 Feb 2026 16:17:00 +1100 Subject: [PATCH 24/44] Update columnphysics/icepack_therm_vertical.F90 Co-authored-by: Anton Steketee <79179784+anton-seaice@users.noreply.github.com> --- columnphysics/icepack_therm_vertical.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index bb6287fba..c7cfe65c1 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -111,7 +111,7 @@ function cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) re 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_ratio = (cold_temp_flag - top_layer_temp) / (c100 + cold_temp_flag) reduce_amount = reduce_ratio * fcondtopn_solve fcondtopn_solve = fcondtopn_solve - reduce_amount endif From 594f9ff54383b5a8889d3ca571ee2ff49803e40d Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 24 Feb 2026 16:57:48 +1100 Subject: [PATCH 25/44] add extra calc_Tsfc, and add cap_conductive_flux parameters to icepack parameters --- columnphysics/icepack_parameters.F90 | 27 +++++++++++----- columnphysics/icepack_therm_bl99.F90 | 40 +++++++++++++----------- columnphysics/icepack_therm_vertical.F90 | 22 ++++++++----- 3 files changed, 55 insertions(+), 34 deletions(-) diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index fe6ce051e..83d709bc7 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -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,&! max condutive flux/depth ratio (W m) + cold_temp_flag = c0 - 60.0 ! min temp used to limit the conductive flux (C) integer (kind=int_kind), public :: & ktherm = 1 ! type of thermodynamics @@ -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, & @@ -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 @@ -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 if (present(shortwave_in) ) shortwave = shortwave_in if (present(albedo_type_in) ) albedo_type = albedo_type_in if (present(albicev_in) ) albicev = albicev_in @@ -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, & @@ -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 @@ -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 @@ -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 write(iounit,*) " shortwave = ", trim(shortwave) write(iounit,*) " albedo_type= ", trim(albedo_type) write(iounit,*) " albicev = ", albicev diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index f72e6af1f..f54ee1296 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -567,26 +567,28 @@ subroutine temperature_changes (dt, & endif 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. + if (.not. calc_Tsfc) 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 - endif - + end if + zTsn(k) = min(zTsn(k), c0) endif diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index c7cfe65c1..3db9ac34c 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -23,10 +23,11 @@ module icepack_therm_vertical use icepack_fsd, only: floe_rad_c, floe_binwidth - use icepack_parameters, only: c0, c1, c2, p001, p5, puny + use icepack_parameters, only: c0, c1, c2, p001, p5, puny, c100 use icepack_parameters, only: pi, depressT, Lvap, hs_min, cp_ice, min_salin use icepack_parameters, only: cp_ocn, rhow, rhoi, rhos, Lfresh, rhofresh, ice_ref_salinity use icepack_parameters, only: ktherm, calc_Tsfc, rsnw_fall, rsnw_tmax + use icepack_parameters, only: ratio_Wm2_m, cold_temp_flag use icepack_parameters, only: ustar_min, fbot_xfer_type, formdrag, calc_strair use icepack_parameters, only: rfracmin, rfracmax, dpscale, frzpnd, snwgrain, snwlvlfac use icepack_parameters, only: phi_i_mushy, floeshape, floediam, use_smliq_pnd, snwredist @@ -76,6 +77,10 @@ module icepack_therm_vertical !======================================================================= +! +! Function for limiting the conductive flux supplied from an atmosphere model when calc_tsfc=.false. +! + function cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) result(fcondtopn_solve) integer (kind=int_kind), intent(in) :: & @@ -89,8 +94,6 @@ function cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) re 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, & @@ -389,12 +392,15 @@ 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 + if (calc_Tsfc) then + fcondtopn_solve = fcondtopn + else + fcondtopn_solve = cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) + end if + + + fcondtopn_extra = fcondtopn - fcondtopn_solve call temperature_changes(dt, & rhoa, flw, & From 6ad9e4c449b2dc2736076c5943447e5125b737a2 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 24 Feb 2026 17:17:03 +1100 Subject: [PATCH 26/44] add ACCESS flag --- columnphysics/icepack_flux.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/columnphysics/icepack_flux.F90 b/columnphysics/icepack_flux.F90 index c441a7cc3..9ec452e0b 100644 --- a/columnphysics/icepack_flux.F90 +++ b/columnphysics/icepack_flux.F90 @@ -349,9 +349,9 @@ subroutine set_sfcflux (aicen, & character(len=*),parameter :: subname='(set_sfcflux)' - raicen = c1 / aicen + raicen = c1 -#ifdef CICE_IN_NEMO +#if defined(CICE_IN_NEMO) || defined(ACCESS) !---------------------------------------------------------------------- ! Convert fluxes from GBM values to per ice area values when ! running in NEMO environment. (When in standalone mode, fluxes From e92ca4217ccfe89c955b551411cf911cd81f24cf Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 24 Feb 2026 17:27:05 +1100 Subject: [PATCH 27/44] add ACCESS flag --- columnphysics/icepack_flux.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/columnphysics/icepack_flux.F90 b/columnphysics/icepack_flux.F90 index 9ec452e0b..082f6effa 100644 --- a/columnphysics/icepack_flux.F90 +++ b/columnphysics/icepack_flux.F90 @@ -351,7 +351,7 @@ subroutine set_sfcflux (aicen, & raicen = c1 -#if defined(CICE_IN_NEMO) || defined(ACCESS) +#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 From bfd3c4c85b0172785ac15ebd484fab7211a0ce26 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Wed, 25 Feb 2026 15:48:22 +1100 Subject: [PATCH 28/44] extra .not. calc_Tsfc check --- columnphysics/icepack_therm_bl99.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index f54ee1296..3879b5cc9 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -621,7 +621,7 @@ subroutine temperature_changes (dt, & zTin(k) = Tmat(k+1+nslyr) - if (l_brine .and. zTin(k) > Tmlts(k) - puny) then + if (.not. calc_Tsfc .and. l_brine .and. zTin(k) > Tmlts(k) - puny) then dTmat(k) = zTin(k) - Tmlts(k) dqmat(k) = rhoi * dTmat(k) & * (cp_ice - Lfresh * Tmlts(k)/zTin(k)**2) From 7e86b82c6d8315be8d8fbeb2651afc05cc5e6e58 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Wed, 25 Feb 2026 16:25:31 +1100 Subject: [PATCH 29/44] extra .not. calc_Tsfc check --- columnphysics/icepack_therm_bl99.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 3879b5cc9..5c24e6756 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -621,12 +621,12 @@ subroutine temperature_changes (dt, & zTin(k) = Tmat(k+1+nslyr) - if (.not. calc_Tsfc .and. l_brine .and. zTin(k) > Tmlts(k) - puny) then + if (l_brine .and. zTin(k) > Tmlts(k) - puny) then 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 ((.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. From 9b8a3ff55a8869f1138a2a422bfb06d873501b3e Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Wed, 25 Feb 2026 17:07:26 +1100 Subject: [PATCH 30/44] fix fcondtop issue --- columnphysics/icepack_therm_vertical.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 3db9ac34c..6c4b8a897 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -419,6 +419,8 @@ subroutine thermo_vertical (dt, aicen, & einit, e_num) if (icepack_warnings_aborted(subname)) return + fcondtopn = fcondtopn_solve + fcondtopn_extra + endif ! ktherm ! mass of ice and liquid water in snow From 9e0e644ba3e4b519d4e1262a241845e5a5adb900 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Wed, 25 Mar 2026 10:28:50 +1100 Subject: [PATCH 31/44] Update columnphysics/icepack_parameters.F90 Co-authored-by: Anton Steketee <79179784+anton-seaice@users.noreply.github.com> --- columnphysics/icepack_parameters.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index 83d709bc7..fca969297 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -137,8 +137,8 @@ module icepack_parameters 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 - ratio_Wm2_m = 1000.0,&! max condutive flux/depth ratio (W m) - cold_temp_flag = c0 - 60.0 ! min temp used to limit the conductive flux (C) + 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 From 79cef70e21bec3289cc60f9ee5de1c4ce3fd34b0 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Wed, 25 Mar 2026 10:29:25 +1100 Subject: [PATCH 32/44] add .not. calc_Tsfc check for thermo conservation error logging --- columnphysics/icepack_therm_vertical.F90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 6c4b8a897..7a3096d54 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -2120,10 +2120,12 @@ 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) + if (.not. calc_Tsfc) then + 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) + end if write(warnstr,*) subname, 'fbot,fcondbot:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, fbot,fcondbotn From ebed003d82c1900241aacee7a1d4318f18a7f92a Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 14 Apr 2026 16:12:48 +1000 Subject: [PATCH 33/44] formatting and cleanup --- columnphysics/icepack_therm_bl99.F90 | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 5c24e6756..cd9539507 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -123,7 +123,8 @@ subroutine temperature_changes (dt, & zTsn ! internal snow layer temperatures real (kind=dbl_kind), intent(out):: & - e_num + e_num ! excess energy from conductive flux (J m-2) + ! local variables integer (kind=int_kind), parameter :: & @@ -568,11 +569,11 @@ subroutine temperature_changes (dt, & if ((l_brine) .and. zTsn(k)>c0) then if (.not. calc_Tsfc) then - ! Alex West: return this energy to the ocean + ! 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 + ! 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 @@ -810,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) From ed835b8f8f1568490613109a29077bdd78aedd8a Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 11 May 2026 15:58:34 +1000 Subject: [PATCH 34/44] update einex and e_num variable names --- columnphysics/icepack_therm_bl99.F90 | 22 +++++++++++----------- columnphysics/icepack_therm_vertical.F90 | 22 +++++++++++----------- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index cd9539507..17b1f09a1 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -69,7 +69,7 @@ subroutine temperature_changes (dt, & fsensn, flatn, & flwoutn, fsurfn, & fcondtopn,fcondbot, & - einit, e_num) + einit, einex_sfc_flux) real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -123,7 +123,7 @@ subroutine temperature_changes (dt, & zTsn ! internal snow layer temperatures real (kind=dbl_kind), intent(out):: & - e_num ! excess energy from conductive flux (J m-2) + einex_sfc_flux ! excess energy from conductive flux (J m-2) ! local variables @@ -158,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) :: & @@ -208,7 +208,7 @@ subroutine temperature_changes (dt, & ! Initialize !----------------------------------------------------------------- fcondtopn_reduction = c0 - e_num = c0 + einex_sfc_flux = c0 Top_T_was_reset_last_time = .false. converged = .false. l_snow = .false. @@ -219,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 @@ -344,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. @@ -583,7 +583,7 @@ subroutine temperature_changes (dt, & 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 + einex_sfc_flux = einex_sfc_flux + hslyr * dqmat_sn else Top_T_was_reset_last_time = .true. endif @@ -631,7 +631,7 @@ subroutine temperature_changes (dt, & 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) + einex_sfc_flux = einex_sfc_flux + hilyr * dqmat(k) else Top_T_was_reset_last_time = .true. endif @@ -681,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 @@ -729,11 +729,11 @@ subroutine temperature_changes (dt, & if (calc_Tsfc) then ! Flux extra energy out of the ice - fcondbot = fcondbot + einex/dt + fcondbot = fcondbot + einex_sfc_calc/dt ferr = abs( (enew-einit)/dt & - (fcondtopn - fcondbot + fswint) ) else - ferr = abs( (enew-einit+e_num)/dt & + ferr = abs( (enew-einit+einex_sfc_flux)/dt & - (fcondtopn - fcondbot + fswint) ) end if diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 7a3096d54..a61b71ac5 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -293,7 +293,7 @@ subroutine thermo_vertical (dt, aicen, & fadvocn, saltvol, dfsalt ! advective heat flux to ocean real (kind=dbl_kind) :: & - fcondtopn_solve, fcondtopn_extra, e_num + fcondtopn_solve, fcondtopn_extra, einex_sfc_flux character(len=*),parameter :: subname='(thermo_vertical)' @@ -321,7 +321,7 @@ subroutine thermo_vertical (dt, aicen, & meltsliq= c0 massice(:) = c0 massliq(:) = c0 - e_num = c0 + einex_sfc_flux = c0 fcondtopn_extra = c0 if (tr_pond) then dpnd_flush = c0 @@ -416,7 +416,7 @@ subroutine thermo_vertical (dt, aicen, & fsensn, flatn, & flwoutn, fsurfn, & fcondtopn_solve, fcondbotn, & - einit, e_num) + einit, einex_sfc_flux) if (icepack_warnings_aborted(subname)) return fcondtopn = fcondtopn_solve + fcondtopn_extra @@ -476,7 +476,7 @@ subroutine thermo_vertical (dt, aicen, & zSin, sss, & sst, & dsnow, rsnw, & - e_num, fcondtopn_extra ) + einex_sfc_flux, fcondtopn_extra ) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- @@ -490,7 +490,7 @@ subroutine thermo_vertical (dt, aicen, & fsnow, einit, & einter, efinal, & fcondtopn, fcondbotn, & - fadvocn, fbot, e_num, fcondtopn_extra ) + fadvocn, fbot, einex_sfc_flux, fcondtopn_extra ) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- @@ -1136,7 +1136,7 @@ subroutine thickness_changes (dt, yday, & zSin, sss, & sst, & dsnow, rsnw, & - e_num, fcondtopn_extra) + einex_sfc_flux, fcondtopn_extra) real (kind=dbl_kind), intent(in) :: & dt , & ! time step @@ -1205,7 +1205,7 @@ subroutine thickness_changes (dt, yday, & sss ! ocean salinity (PSU) real (kind=dbl_kind), intent(in) :: & - e_num, fcondtopn_extra + einex_sfc_flux, fcondtopn_extra ! local variables integer (kind=int_kind) :: & @@ -1629,7 +1629,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 + e_num)/dt + fhocnn = fbot + (esub + etop_mlt + ebot_mlt + einex_sfc_flux)/dt !----------------------------------------------------------------- ! Add new snowfall at top surface @@ -2054,7 +2054,7 @@ subroutine conservation_check_vthermo(dt, & einit, einter, & efinal, & fcondtopn,fcondbotn, & - fadvocn, fbot, e_num, fcondtopn_extra) + fadvocn, fbot, einex_sfc_flux, fcondtopn_extra) real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -2067,7 +2067,7 @@ subroutine conservation_check_vthermo(dt, & fsnow , & ! snowfall rate (kg m-2 s-1) fcondtopn , & fadvocn , & - fbot, e_num, fcondtopn_extra + fbot, einex_sfc_flux, fcondtopn_extra real (kind=dbl_kind), intent(in) :: & einit , & ! initial energy of melting (J m-2) @@ -2121,7 +2121,7 @@ subroutine conservation_check_vthermo(dt, & write(warnstr,*) subname, 'Input energy =', einp call icepack_warnings_add(warnstr) if (.not. calc_Tsfc) then - write(warnstr,*) subname, 'Numerical energy =', e_num + write(warnstr,*) subname, 'Numerical energy =', einex_sfc_flux call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'fcondtopn_extra energy =', fcondtopn_extra call icepack_warnings_add(warnstr) From cda0f140b1256b10eaaff6993c3f032279d4e06c Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 16 Jun 2026 16:56:24 +1000 Subject: [PATCH 35/44] comments with units for new variables + formatting --- columnphysics/icepack_therm_bl99.F90 | 20 +++++---- columnphysics/icepack_therm_vertical.F90 | 56 ++++++++++++------------ 2 files changed, 39 insertions(+), 37 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 17b1f09a1..844fb6ba3 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -187,17 +187,19 @@ subroutine temperature_changes (dt, & kh ! effective conductivity at interfaces (W m-2 deg-1) real (kind=dbl_kind) :: & - ci , & ! specific heat of sea ice (J kg-1 deg-1) - avg_Tsf , & ! = 1. if Tsf averaged w/Tsf_start, else = 0. - Iswabs_tmp , & ! energy to melt through fraction frac of layer - Sswabs_tmp , & ! same for snow - dswabs , & ! difference in swabs and swabs_tmp - frac , & - fcondtopn_reduction, & - fcondtopn_force, dqmat_sn + ci , & ! specific heat of sea ice (J kg-1 deg-1) + avg_Tsf , & ! = 1. if Tsf averaged w/Tsf_start, else = 0. + Iswabs_tmp , & ! energy to melt through fraction frac of layer + Sswabs_tmp , & ! same for snow + dswabs , & ! difference in swabs and swabs_tmp + frac , & + fcondtopn_reduction, & ! reduction in downward cond flux at top surface (W m-2) + fcondtopn_force , & ! reduced downward cond flux at top surface (W m-2) + dqmat_sn ! associated enthalpy difference at top snow layer (J m-3) logical (kind=log_kind) :: & - converged, Top_T_was_reset_last_time ! = true when local solution has converged + converged , & ! = true when local solution has converged + Top_T_was_reset_last_time ! = true when surface temperature was reset in prev iteration logical (kind=log_kind) , dimension (nilyr) :: & reduce_kh ! reduce conductivity when T exceeds Tmlt diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index a61b71ac5..52919c5da 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -77,32 +77,27 @@ module icepack_therm_vertical !======================================================================= -! ! Function for limiting the conductive flux supplied from an atmosphere model when calc_tsfc=.false. -! 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 - - ! AEW: New variables for cold-ice flux capping - real (kind=dbl_kind) :: top_layer_temp, & - reduce_ratio, & - reduce_amount - + real (kind=dbl_kind), intent(in) :: fcondtopn ! downward cond flux at top surface (W m-2) + real (kind=dbl_kind), intent(in) :: hin ! ice thickness (m) + real (kind=dbl_kind), intent(in) :: zTin(nilyr) ! internal ice layer temperatures (C) + real (kind=dbl_kind), intent(in) :: zTsn(nslyr) ! internal snow layer temperatures (C) + real (kind=dbl_kind), intent(in) :: hslyr ! snow layer thickness (m) + + real (kind=dbl_kind) :: fcondtopn_solve ! limited downward cond flux at top surface (W m-2) + real (kind=dbl_kind) :: & + top_layer_temp, & ! top layer temperature (C) + reduce_ratio , & ! reduction ratio of downward cond flux at top surface + reduce_amount ! reduction in downward cond flux at top surface (W m-2) if (abs(fcondtopn) > ratio_Wm2_m * hin) then fcondtopn_solve = sign(ratio_Wm2_m * hin,fcondtopn) - else fcondtopn_solve = fcondtopn endif @@ -121,6 +116,8 @@ function cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) re end function cap_conductive_flux +!======================================================================= + ! ! Driver for updating ice and snow internal temperatures and ! computing thermodynamic growth rates and atmospheric fluxes. @@ -293,7 +290,9 @@ subroutine thermo_vertical (dt, aicen, & fadvocn, saltvol, dfsalt ! advective heat flux to ocean real (kind=dbl_kind) :: & - fcondtopn_solve, fcondtopn_extra, einex_sfc_flux + fcondtopn_solve, & ! limited downward cond flux at top surface (W m-2) + fcondtopn_extra, & ! excess downward cond flux at top surface (W m-2) + einex_sfc_flux ! excess energy removed during temperature_changes (J m-2) character(len=*),parameter :: subname='(thermo_vertical)' @@ -398,8 +397,6 @@ subroutine thermo_vertical (dt, aicen, & else fcondtopn_solve = cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) end if - - fcondtopn_extra = fcondtopn - fcondtopn_solve call temperature_changes(dt, & @@ -1205,7 +1202,8 @@ subroutine thickness_changes (dt, yday, & sss ! ocean salinity (PSU) real (kind=dbl_kind), intent(in) :: & - einex_sfc_flux, fcondtopn_extra + einex_sfc_flux, & ! excess energy removed during temperature_changes (J m-2) + fcondtopn_extra ! excess downward cond flux at top surface (W m-2) ! local variables integer (kind=int_kind) :: & @@ -2060,14 +2058,16 @@ subroutine conservation_check_vthermo(dt, & dt ! time step real (kind=dbl_kind), intent(in) :: & - fsurfn , & ! net flux to top surface, excluding fcondtopn - flatn , & ! surface downward latent heat (W m-2) - fhocnn , & ! fbot, corrected for any surplus energy - fswint , & ! SW absorbed in ice interior, below surface (W m-2) - fsnow , & ! snowfall rate (kg m-2 s-1) - fcondtopn , & - fadvocn , & - fbot, einex_sfc_flux, fcondtopn_extra + fsurfn , & ! net flux to top surface, excluding fcondtopn + flatn , & ! surface downward latent heat (W m-2) + fhocnn , & ! fbot, corrected for any surplus energy + fswint , & ! SW absorbed in ice interior, below surface (W m-2) + fsnow , & ! snowfall rate (kg m-2 s-1) + fcondtopn , & ! downward cond flux at top surface (W m-2) + fadvocn , & ! advective heat flux to ocean (W m-2) + fbot , & ! ice-ocean heat flux at bottom surface (W/m^2) + einex_sfc_flux, & ! excess energy removed during temperature_changes (J m-2) + fcondtopn_extra ! excess downward cond flux at top surface (W m-2) real (kind=dbl_kind), intent(in) :: & einit , & ! initial energy of melting (J m-2) From edbb4ce7af9a76eef816f973c234bbcf0f55944c Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 16 Jun 2026 16:58:08 +1000 Subject: [PATCH 36/44] comments with units for new variables + formatting --- columnphysics/icepack_therm_bl99.F90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 844fb6ba3..a5f9f1cee 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -152,14 +152,14 @@ subroutine temperature_changes (dt, & enew ! new energy of melting after temp change (J m-2) real (kind=dbl_kind) :: & - dTsf_prev , & ! dTsf from previous iteration - dTi1_prev , & ! dTi1 from previous iteration - dfsens_dT , & ! deriv of fsens wrt Tsf (W m-2 deg-1) - 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_sfc_calc, & ! excess energy from dqmat to ocean - ferr ! energy conservation error (W m-2) + dTsf_prev , & ! dTsf from previous iteration + dTi1_prev , & ! dTi1 from previous iteration + dfsens_dT , & ! deriv of fsens wrt Tsf (W m-2 deg-1) + 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_sfc_calc, & ! excess energy from dqmat to ocean (J m-2) + ferr ! energy conservation error (W m-2) real (kind=dbl_kind), dimension (nilyr) :: & Tin_init , & ! zTin at beginning of time step From a1d5fac4ccc0e777871cc76bbb53050fc11346d4 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 16 Jun 2026 16:59:29 +1000 Subject: [PATCH 37/44] formatting --- columnphysics/icepack_therm_bl99.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index a5f9f1cee..36b8fb540 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -202,7 +202,7 @@ subroutine temperature_changes (dt, & Top_T_was_reset_last_time ! = true when surface temperature was reset in prev iteration 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)' From 630c13d8415c7d640ae65890fc65ddf11e8b4e9a Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Wed, 1 Jul 2026 18:09:00 +1000 Subject: [PATCH 38/44] add conductive flux limiting documentation --- doc/source/science_guide/sg_boundary_forcing.rst | 5 ++--- doc/source/science_guide/sg_thermo.rst | 8 +++++++- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/doc/source/science_guide/sg_boundary_forcing.rst b/doc/source/science_guide/sg_boundary_forcing.rst index e73a13a93..17b56f25e 100755 --- a/doc/source/science_guide/sg_boundary_forcing.rst +++ b/doc/source/science_guide/sg_boundary_forcing.rst @@ -81,9 +81,8 @@ subroutine, *thermo_vertical*. At the end of the time step, the surface temperature and effective conductivity (i.e., thermal conductivity divided by thickness) of the top ice/snow layer in each category are returned to the atmosphere model via the coupler. Since the ice surface -temperature is treated explicitly, the effective conductivity may need -to be limited to ensure stability. As a result, accuracy may be -significantly reduced, especially for thin ice or snow layers. A more +temperature is treated explicitly, the conductive flux is limited to ensure stability. +As a result, accuracy may be significantly reduced, especially for thin ice or snow layers. A more stable and accurate procedure would be to compute the temperature profiles for both the atmosphere and ice, together with the surface fluxes, in a single implicit calculation. This was judged impractical, diff --git a/doc/source/science_guide/sg_thermo.rst b/doc/source/science_guide/sg_thermo.rst index bf3d3c251..ea30a5038 100755 --- a/doc/source/science_guide/sg_thermo.rst +++ b/doc/source/science_guide/sg_thermo.rst @@ -1375,9 +1375,15 @@ is that the solution scheme is no longer unconditionally stable. Instead, the effective conductivity in the top layer must satisfy a diffusive CFL condition: -.. math:: +.. math:: K^* \le {\rho ch \over \Delta t}. +To improve stability, the conductive flux is limited before +the thermodynamic solve. Any excess flux is then used for bottom melt. The +conductive flux may also be limited during the thermodynamic solve if the +surface temperature exceeds the melting point, in which case the excess flux +is passed to the ocean. + For thin layers and typical coupling intervals (:math:`\sim 1` hr), :math:`K^*` may need to be limited before being passed to the atmosphere via the coupler. Otherwise, the fluxes that are returned to the host sea ice model may From 0cc7ff49b03106fcd0b7e1ffc05e3ac725fc18ca Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 13 Jul 2026 17:07:24 +1000 Subject: [PATCH 39/44] add warnings back in --- columnphysics/icepack_therm_bl99.F90 | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 36b8fb540..f3e7bfebd 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -812,6 +812,14 @@ subroutine temperature_changes (dt, & write(warnstr,*) subname, 'Final snow temperatures:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (zTsn(k),k=1,nslyr) + 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) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Final ice temperatures:' call icepack_warnings_add(warnstr) From b88bc33839cead1f585400991d79feb3e807f13a Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 13 Jul 2026 17:07:59 +1000 Subject: [PATCH 40/44] add warnings back in --- columnphysics/icepack_therm_bl99.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index f3e7bfebd..0c8d555b6 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -812,6 +812,7 @@ subroutine temperature_changes (dt, & write(warnstr,*) subname, 'Final snow temperatures:' 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) From 0005fff428307aa137a7279bed3d8b91d01a5707 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 13 Jul 2026 17:08:09 +1000 Subject: [PATCH 41/44] add warnings back in --- columnphysics/icepack_therm_bl99.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 0c8d555b6..7dc0e3093 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -821,7 +821,6 @@ subroutine temperature_changes (dt, & call icepack_warnings_add(warnstr) write(warnstr,*) subname, (hilyr*dqmat(k)/dt,k=1,nilyr) call icepack_warnings_add(warnstr) - call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Final ice temperatures:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (zTin(k),k=1,nilyr) From 49b04c6143f89e0e014efa7b051cebf638ba343f Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 13 Jul 2026 17:20:43 +1000 Subject: [PATCH 42/44] documentation changes --- doc/source/science_guide/sg_boundary_forcing.rst | 5 ++++- doc/source/science_guide/sg_thermo.rst | 2 +- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/doc/source/science_guide/sg_boundary_forcing.rst b/doc/source/science_guide/sg_boundary_forcing.rst index 17b56f25e..c4cb7e308 100755 --- a/doc/source/science_guide/sg_boundary_forcing.rst +++ b/doc/source/science_guide/sg_boundary_forcing.rst @@ -81,7 +81,8 @@ subroutine, *thermo_vertical*. At the end of the time step, the surface temperature and effective conductivity (i.e., thermal conductivity divided by thickness) of the top ice/snow layer in each category are returned to the atmosphere model via the coupler. Since the ice surface -temperature is treated explicitly, the conductive flux is limited to ensure stability. +temperature is treated explicitly, the conductive flux is limited to ensure stability, +with any excess energy either used for bottom melting or transferred to the ocean as a heat flux. As a result, accuracy may be significantly reduced, especially for thin ice or snow layers. A more stable and accurate procedure would be to compute the temperature profiles for both the atmosphere and ice, together with the surface @@ -96,6 +97,8 @@ implicitly. The resultant surface temperature change is passed back to the atmosphere model via coupler to complete the full update of its temperature profiles. This middle-ground approach, enabled by ``vapor_flux_correction=true``, does not sacrifice accuracy because it does not need effective conductivity limiting as in the explicit case. +The GEOS approach is available for either mushy or BL99 thermodynamics, while the +flux-limiting approaches are only available for BL99. diff --git a/doc/source/science_guide/sg_thermo.rst b/doc/source/science_guide/sg_thermo.rst index ea30a5038..363553128 100755 --- a/doc/source/science_guide/sg_thermo.rst +++ b/doc/source/science_guide/sg_thermo.rst @@ -1378,7 +1378,7 @@ diffusive CFL condition: .. math:: K^* \le {\rho ch \over \Delta t}. -To improve stability, the conductive flux is limited before +To improve stability when ``calc_Tsfc=false``, the conductive flux is limited before the thermodynamic solve. Any excess flux is then used for bottom melt. The conductive flux may also be limited during the thermodynamic solve if the surface temperature exceeds the melting point, in which case the excess flux From 19ff472799b8beff721c2c3f91374d7b4345bbb9 Mon Sep 17 00:00:00 2001 From: anton-seaice Date: Tue, 21 Jul 2026 13:43:04 +1000 Subject: [PATCH 43/44] line up formatting --- columnphysics/icepack_parameters.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index fca969297..d75909d17 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -137,8 +137,8 @@ module icepack_parameters 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 - 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) + 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 From 6cf77f7ab107b65ec4d8d470a6f7607fe7dc98e3 Mon Sep 17 00:00:00 2001 From: anton-seaice Date: Tue, 21 Jul 2026 13:43:22 +1000 Subject: [PATCH 44/44] Test setup --- .../scripts/machines/Macros.gadi_oneapi | 44 +++++++++++++++++++ .../scripts/machines/env.gadi_oneapi | 33 ++++++++++++++ 2 files changed, 77 insertions(+) create mode 100644 configuration/scripts/machines/Macros.gadi_oneapi create mode 100644 configuration/scripts/machines/env.gadi_oneapi diff --git a/configuration/scripts/machines/Macros.gadi_oneapi b/configuration/scripts/machines/Macros.gadi_oneapi new file mode 100644 index 000000000..eacc6cb0a --- /dev/null +++ b/configuration/scripts/machines/Macros.gadi_oneapi @@ -0,0 +1,44 @@ +#============================================================================== +# Makefile macros for NCI Gadi, intel compiler +#============================================================================== + +CPP := fpp +CPPDEFS := -DFORTRANUNDERSCORE -DREPRODUCIBLE ${ICE_CPPDEFS} +CFLAGS := -c -O2 -fp-model precise -Wno-unused-variable -Wno-unused-parameter + +FIXEDFLAGS := -132 +FREEFLAGS := -FR + +NCI_INTEL_FLAGS := -r8 -i4 -traceback -w -fpe0 -ftz -convert big_endian -assume byterecl -check noarg_temp_created +NCI_REPRO_FLAGS := -fp-model precise -fp-model source -align all + +ifeq ($(ICE_BLDDEBUG), true) + NCI_DEBUG_FLAGS := -g3 -O0 -debug all -check all -no-vec -assume nobuffered_io + FFLAGS := $(NCI_INTEL_FLAGS) $(NCI_REPRO_FLAGS) $(NCI_DEBUG_FLAGS) + CPPDEFS := $(CPPDEFS) -DDEBUG=$(DEBUG) +else + NCI_OPTIM_FLAGS := -g3 -O2 + FFLAGS := $(NCI_INTEL_FLAGS) $(NCI_REPRO_FLAGS) $(NCI_OPTIM_FLAGS) +endif + +SCC := icx +SFC := ifx +MPICC := mpicc +MPIFC := mpifort + +ifeq ($(ICE_COMMDIR), mpi) + FC := $(MPIFC) + CC := $(MPICC) +else + FC := $(SFC) + CC := $(SCC) +endif +LD:= $(FC) + +SLIBS := $(SLIBS) +INCLDIR := $(INCLDIR) + +# use system netcdf build +SLIBS += -L$(NETCDF)/lib -lnetcdf -lnetcdff +INCLDIR += -I$(NETCDF)/include + diff --git a/configuration/scripts/machines/env.gadi_oneapi b/configuration/scripts/machines/env.gadi_oneapi new file mode 100644 index 000000000..b3c387cea --- /dev/null +++ b/configuration/scripts/machines/env.gadi_oneapi @@ -0,0 +1,33 @@ +#!/bin/csh -f + +set inp = "undefined" +if ($#argv == 1) then + set inp = $1 +endif + +if ("$inp" != "-nomodules") then + + source /etc/profile.d/modules.csh + + module load intel-compiler-llvm + module load netcdf + module load openmpi + +endif + +setenv ICE_MACHINE_MACHNAME gadi +setenv ICE_MACHINE_MACHINFO "Intel Xeon Scalable" +setenv ICE_MACHINE_ENVNAME intel +setenv ICE_MACHINE_ENVINFO INTEL_COMPILER_VERSION=$INTEL_COMPILER_LLVM_BASE +setenv ICE_MACHINE_MAKE gmake +setenv ICE_MACHINE_WKDIR /scratch/$PROJECT/$USER/ICEPACK_RUNS +setenv ICE_MACHINE_INPUTDATA /g/data/ik11/inputs/CICE_data/icepack-dirs/input +setenv ICE_MACHINE_BASELINE /scratch/$PROJECT/$USER/ICEPACK_BASELINE +setenv ICE_MACHINE_SUBMIT "qsub" +setenv ICE_MACHINE_PROJ $PROJECT +setenv ICE_MACHINE_ACCT $USER +setenv ICE_MACHINE_QUEUE "normal" +setenv ICE_MACHINE_TPNODE 48 +setenv ICE_MACHINE_BLDTHRDS 4 +setenv ICE_MACHINE_QSTAT "qstat" +setenv ICE_CPPDEFS -DUSE_NETCDF \ No newline at end of file