From bc088acf63b5bbfa04745bb7bba80258749134ce Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 29 Apr 2024 14:33:19 +1000 Subject: [PATCH 01/21] 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 f9485e32..8b798e59 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 73d560b6..8c1116fe 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 00b3756d..758ed0d3 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/21] 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 8b798e59..e93da403 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 758ed0d3..2bd9b9bb 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/21] 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 e93da403..f2eab66f 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 2bd9b9bb..80c04f71 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/21] 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 f2eab66f..e40c9dd1 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 80c04f71..d0f06457 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/21] 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 e40c9dd1..008f66a9 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 8c1116fe..73d560b6 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 d0f06457..e9ce400f 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/21] 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 e9ce400f..d0f06457 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/21] 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 d5578076..c441a7cc 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/21] 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 c441a7cc..bca86fc1 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/21] 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 8802a943..df17671b 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 fe6ce051..69e61116 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/21] 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 bca86fc1..c441a7cc 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/21] 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 008f66a9..9ac9ed31 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/21] 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 69e61116..9a9ea818 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/21] 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 9ac9ed31..55d35cae 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 d0f06457..a022dc51 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/21] 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 df17671b..8802a943 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 55d35cae..825bbcf3 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/21] 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 825bbcf3..68c4e294 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 a022dc51..bb6287fb 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/21] 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 68c4e294..3ec44591 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/21] 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 3ec44591..d2212282 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/21] 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 d2212282..5809c06e 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/21] 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 9a9ea818..301d05bb 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/21] 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 301d05bb..fe6ce051 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/21] 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 133683dd..e5ab9800 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, &