From 0b5be471f844671a11bb09c901385ce5256a9db0 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 29 Apr 2024 14:33:19 +1000 Subject: [PATCH 01/19] logging, tweaked newton solver tolerances, re-add einex --- columnphysics/icepack_therm_bl99.F90 | 168 +++++++++++++++++------ columnphysics/icepack_therm_shared.F90 | 2 +- columnphysics/icepack_therm_vertical.F90 | 135 +++++++++++++++--- 3 files changed, 244 insertions(+), 61 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 30c247ac..135eccd2 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 @@ -67,7 +68,7 @@ subroutine temperature_changes (dt, & fsensn, flatn, & flwoutn, fsurfn, & fcondtopn,fcondbot, & - einit ) + einit, e_num) real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -119,7 +120,9 @@ subroutine temperature_changes (dt, & real (kind=dbl_kind), dimension (nslyr), intent(inout) :: & 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 :: & @@ -187,20 +190,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. @@ -255,53 +262,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 !----------------------------------------------------------------- ! Solve for new temperatures. @@ -414,6 +424,7 @@ subroutine temperature_changes (dt, & else + fcondtopn_force = fcondtopn - fcondtopn_reduction call get_matrix_elements_know_Tsfc ( & l_snow, Tbot, & Tin_init, Tsn_init, & @@ -422,7 +433,7 @@ subroutine temperature_changes (dt, & etai, etas, & sbdiag, diag, & spdiag, rhs, & - fcondtopn) + fcondtopn_force) if (icepack_warnings_aborted(subname)) return endif ! calc_Tsfc @@ -544,7 +555,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 @@ -578,6 +614,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)) @@ -669,11 +715,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. @@ -699,76 +745,110 @@ subroutine temperature_changes (dt, & if (.not.converged) then write(warnstr,*) subname, 'Thermo iteration does not converge,' call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Ice thickness:', hilyr*nilyr call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Snow thickness:', hslyr*nslyr call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'dTsf, Tsf_errmax:',dTsf_prev, & Tsf_errmax call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Tsf:', Tsf call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'fsurf:', fsurfn call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'fcondtop, fcondbot, fswint', & fcondtopn, fcondbot, fswint call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'fswsfc', fswsfc call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Iswabs',(Iswabs(k),k=1,nilyr) call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Flux conservation error =', ferr call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Initial snow temperatures:' call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, (Tsn_init(k),k=1,nslyr) call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Initial ice temperatures:' call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, (Tin_init(k),k=1,nilyr) call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Matrix ice temperature diff:' call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, (dTmat(k),k=1,nilyr) call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'dqmat*hilyr/dt:' call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, (hilyr*dqmat(k)/dt,k=1,nilyr) call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Final snow temperatures:' call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, (zTsn(k),k=1,nslyr) call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Matrix ice temperature diff:' call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, (dTmat(k),k=1,nilyr) call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'dqmat*hilyr/dt:' call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, (hilyr*dqmat(k)/dt,k=1,nilyr) call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Final ice temperatures:' call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, (zTin(k),k=1,nilyr) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Ice melting temperatures:' call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, (Tmlts(k),k=1,nilyr) call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Ice bottom temperature:', Tbot call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'dT initial:' call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, (Tmlts(k)-Tin_init(k),k=1,nilyr) call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'dT final:' call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, (Tmlts(k)-zTin(k),k=1,nilyr) call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'zSin' call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, (zSin(k),k=1,nilyr) call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//" temperature_changes: Thermo iteration does not converge" ) return diff --git a/columnphysics/icepack_therm_shared.F90 b/columnphysics/icepack_therm_shared.F90 index 858e6c56..7643de80 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 e931629c..494fba60 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -62,6 +62,8 @@ module icepack_therm_vertical use icepack_meltpond_topo, only: compute_ponds_topo use icepack_snow, only: drain_snow + use ESMF + implicit none private @@ -72,6 +74,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. @@ -240,6 +285,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)' !----------------------------------------------------------------- @@ -266,6 +314,8 @@ subroutine thermo_vertical (dt, aicen, & meltsliq= c0 massice(:) = c0 massliq(:) = c0 + e_num = c0 + fcondtopn_extra = c0 if (calc_Tsfc) then fsensn = c0 @@ -274,6 +324,8 @@ subroutine thermo_vertical (dt, aicen, & fcondtopn = c0 endif + fcondtopn_solve = fcondtopn + !----------------------------------------------------------------- ! Compute variables needed for vertical thermo calculation !----------------------------------------------------------------- @@ -328,6 +380,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, & @@ -342,8 +400,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 @@ -390,7 +448,7 @@ subroutine thermo_vertical (dt, aicen, & smliq, massliq, & fbot, Tbot, & flatn, fsurfn, & - fcondtopn, fcondbotn, & + fcondtopn_solve, fcondbotn, & fsnow, hsn_new, & fhocnn, evapn, & evapsn, evapin, & @@ -401,7 +459,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 !----------------------------------------------------------------- @@ -415,7 +474,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 !----------------------------------------------------------------- @@ -1059,8 +1118,9 @@ subroutine thickness_changes (dt, yday, & congel, snoice, & mlt_onset, frz_onset,& zSin, sss, & - sst, & - dsnow, rsnw) + sst, & + dsnow, rsnw, & + e_num, fcondtopn_extra) real (kind=dbl_kind), intent(in) :: & dt , & ! time step @@ -1127,7 +1187,9 @@ subroutine thickness_changes (dt, yday, & real (kind=dbl_kind), intent(in) :: & 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) :: & @@ -1259,7 +1321,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 @@ -1511,8 +1574,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 @@ -1934,7 +1997,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 @@ -1947,7 +2010,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) @@ -1960,7 +2023,7 @@ subroutine conservation_check_vthermo(dt, & real (kind=dbl_kind) :: & einp , & ! energy input during timestep (J m-2) ferr ! energy conservation error (W m-2) - + character(len=*),parameter :: subname='(conservation_check_vthermo)' !---------------------------------------------------------------- @@ -1977,9 +2040,45 @@ subroutine conservation_check_vthermo(dt, & einp = (fsurfn - flatn + fswint - fhocnn & - fsnow*Lfresh - fadvocn) * dt ferr = abs(efinal-einit-einp) / dt + + ! if (ferr > 10.0_dbl_kind*ferrmax) then + ! call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + ! call icepack_warnings_add(subname//" conservation_check_vthermo: Thermo energy conservation error" ) + ! end if + + if (ferr > 10.0_dbl_kind*ferrmax) then + + write(warnstr,*) subname, 'Thermo energy conservation error' + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Flux error (W/m^2) =', ferr + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Energy error (J) =', ferr*dt + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Initial energy =', einit + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Final energy =', efinal + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'efinal - einit =', efinal-einit + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'fsurfn,flatn,fswint,fhocn, fsnow*Lfresh:' + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, fsurfn,flatn,fswint,fhocnn, fsnow*Lfresh + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Input energy =', einp + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Numerical energy =', e_num + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'fcondtopn_extra energy =', fcondtopn_extra + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'fbot,fcondbot:' + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, fbot,fcondbotn + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + end if - if (ferr > 1.1_dbl_kind*ferrmax) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + ! if (ferr > 1.1_dbl_kind*ferrmax) then + if (.false.) then + ! call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//" conservation_check_vthermo: Thermo energy conservation error" ) write(warnstr,*) subname, 'Thermo energy conservation error' @@ -2000,6 +2099,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 03d523508b0e8a37462b683a796c06a0a0bb27c9 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 30 Apr 2024 16:42:51 +1000 Subject: [PATCH 02/19] extra logging --- columnphysics/icepack_therm_vertical.F90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 494fba60..6ea04a6b 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -2047,7 +2047,7 @@ subroutine conservation_check_vthermo(dt, & ! end if if (ferr > 10.0_dbl_kind*ferrmax) then - + call ESMF_LogWrite("\n\n----------------------------", ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Thermo energy conservation error' call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Flux error (W/m^2) =', ferr @@ -2066,6 +2066,10 @@ subroutine conservation_check_vthermo(dt, & call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Input energy =', einp call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'fcondtopn =', fcondtopn + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Intermediate energy =', einter + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'Numerical energy =', e_num call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, 'fcondtopn_extra energy =', fcondtopn_extra @@ -2074,6 +2078,7 @@ subroutine conservation_check_vthermo(dt, & call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) write(warnstr,*) subname, fbot,fcondbotn call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + call ESMF_LogWrite("----------------------------\n\n", ESMF_LOGMSG_INFO) end if ! if (ferr > 1.1_dbl_kind*ferrmax) then From 80dc4024af9677fed7fa82a8c155d8087bef29ba Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 6 May 2024 12:21:12 +1000 Subject: [PATCH 03/19] fix argument ordering in set_sfcflux --- columnphysics/icepack_therm_bl99.F90 | 2 +- columnphysics/icepack_therm_vertical.F90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 135eccd2..e77cd80e 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -713,7 +713,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 6ea04a6b..4556446d 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -1319,6 +1319,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 @@ -2795,7 +2796,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 @@ -2805,8 +2805,8 @@ subroutine icepack_step_therm1(dt, & ! hadgem routine sets fluxes to default values in ice-only mode call set_sfcflux(aicen (n), & flatn_f (n), fsensn_f (n), & - fcondtopn_f(n), & fsurfn_f (n), & + fcondtopn_f(n), & flatn (n), fsensn (n), & fsurfn (n), & fcondtopn (n)) From 3d50087639d9101c92bf7892f57c93d005d1e669 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 7 May 2024 16:16:53 +1000 Subject: [PATCH 04/19] remove cm2 style conductive flux limiting --- columnphysics/icepack_therm_bl99.F90 | 20 ++++++++++---------- columnphysics/icepack_therm_vertical.F90 | 8 ++++---- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index e77cd80e..521cfadb 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -568,15 +568,15 @@ subroutine temperature_changes (dt, & ! 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 + ! 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) @@ -713,7 +713,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 4556446d..dd00d680 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -380,8 +380,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 @@ -400,7 +400,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 @@ -448,7 +448,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 1b539cf9ee659b2230bccce3cdbef7319ce5f10a Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 20 May 2024 11:46:28 +1000 Subject: [PATCH 05/19] use cm2 style conductive flux limiting --- columnphysics/icepack_therm_bl99.F90 | 20 ++++++++++---------- columnphysics/icepack_therm_vertical.F90 | 6 +++--- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 521cfadb..e77cd80e 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -568,15 +568,15 @@ subroutine temperature_changes (dt, & ! 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 + 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) @@ -713,7 +713,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 dd00d680..7ff2bb2f 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -380,8 +380,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 @@ -400,7 +400,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 4d4c7b1287bddfc9b0df6224e60a7f6545e2d699 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 22 Jul 2024 13:56:26 +1000 Subject: [PATCH 06/19] 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 e77cd80e..725e458f 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -719,7 +719,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 7643de80..858e6c56 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 7ff2bb2f..ee10ef75 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -286,7 +286,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)' @@ -382,6 +382,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 @@ -399,7 +400,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 862c39c50567b4d23043c7b7e13e64cc05d47e4a Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 23 Jul 2024 10:10:23 +1000 Subject: [PATCH 07/19] 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 ee10ef75..7ff2bb2f 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -286,7 +286,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)' @@ -382,7 +382,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 @@ -400,7 +399,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 8f336c5673f0254585294719cd171d869e8e10b5 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 20 Aug 2024 14:38:39 +1000 Subject: [PATCH 08/19] convert GBM ice fluxes to local --- columnphysics/icepack_flux.F90 | 2 +- columnphysics/icepack_therm_bl99.F90 | 36 ++++++++++++++++++++++++ columnphysics/icepack_therm_vertical.F90 | 1 + 3 files changed, 38 insertions(+), 1 deletion(-) diff --git a/columnphysics/icepack_flux.F90 b/columnphysics/icepack_flux.F90 index 27a9f0f8..7a25bae0 100644 --- a/columnphysics/icepack_flux.F90 +++ b/columnphysics/icepack_flux.F90 @@ -299,7 +299,7 @@ subroutine set_sfcflux (aicen, & character(len=*),parameter :: subname='(set_sfcflux)' - raicen = c1 + raicen = c1 / aicen #ifdef CICE_IN_NEMO !---------------------------------------------------------------------- diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 725e458f..3fc36dda 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -746,109 +746,145 @@ subroutine temperature_changes (dt, & write(warnstr,*) subname, 'Thermo iteration does not converge,' call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Ice thickness:', hilyr*nilyr call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Snow thickness:', hslyr*nslyr call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'dTsf, Tsf_errmax:',dTsf_prev, & Tsf_errmax call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Tsf:', Tsf call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'fsurf:', fsurfn call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'fcondtop, fcondbot, fswint', & fcondtopn, fcondbot, fswint call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'fswsfc', fswsfc call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Iswabs',(Iswabs(k),k=1,nilyr) call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Flux conservation error =', ferr call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Initial snow temperatures:' call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, (Tsn_init(k),k=1,nslyr) call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Initial ice temperatures:' call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, (Tin_init(k),k=1,nilyr) call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Matrix ice temperature diff:' call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, (dTmat(k),k=1,nilyr) call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'dqmat*hilyr/dt:' call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, (hilyr*dqmat(k)/dt,k=1,nilyr) call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Final snow temperatures:' call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, (zTsn(k),k=1,nslyr) call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Matrix ice temperature diff:' call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, (dTmat(k),k=1,nilyr) call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'dqmat*hilyr/dt:' call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, (hilyr*dqmat(k)/dt,k=1,nilyr) call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Final ice temperatures:' call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, (zTin(k),k=1,nilyr) call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Ice melting temperatures:' call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, (Tmlts(k),k=1,nilyr) call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'Ice bottom temperature:', Tbot call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'dT initial:' call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, (Tmlts(k)-Tin_init(k),k=1,nilyr) call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'dT final:' call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, (Tmlts(k)-zTin(k),k=1,nilyr) call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'zSin' call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, (zSin(k),k=1,nilyr) call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//" temperature_changes: Thermo iteration does not converge" ) return diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 7ff2bb2f..76dff861 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -2080,6 +2080,7 @@ subroutine conservation_check_vthermo(dt, & write(warnstr,*) subname, fbot,fcondbotn call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) call ESMF_LogWrite("----------------------------\n\n", ESMF_LOGMSG_INFO) + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) end if ! if (ferr > 1.1_dbl_kind*ferrmax) then From ca72f1b9e414576626f6e296ab781ed50198e9eb Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Mon, 26 Aug 2024 14:24:18 +1000 Subject: [PATCH 09/19] 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 7a25bae0..e3ed16df 100644 --- a/columnphysics/icepack_flux.F90 +++ b/columnphysics/icepack_flux.F90 @@ -299,7 +299,7 @@ subroutine set_sfcflux (aicen, & character(len=*),parameter :: subname='(set_sfcflux)' - raicen = c1 / aicen + raicen = c1 !/ aicen #ifdef CICE_IN_NEMO !---------------------------------------------------------------------- From e6419abb7fc15e73eb6de6a153b708b3ca4c0e2d Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Wed, 28 Aug 2024 09:34:45 +1000 Subject: [PATCH 10/19] reduce fcondtop limiting threshold --- 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 76dff861..65057a60 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -88,7 +88,7 @@ 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 + real (kind=dbl_kind), parameter :: ratio_Wm2_m = 800.0, cold_temp_flag = c0 - 60.0 ! AEW: New variables for cold-ice flux capping real (kind=dbl_kind) :: top_layer_temp, & From a2d88f3dd9d478006ba7be1b401e0ea407cb33a2 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 17 Sep 2024 11:54:44 +1000 Subject: [PATCH 11/19] 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 +- columnphysics/icepack_therm_bl99.F90 | 4 ++++ columnphysics/icepack_therm_vertical.F90 | 2 +- 4 files changed, 7 insertions(+), 3 deletions(-) diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index 013373ad..669dcced 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 ed76b6b4..1ebf3322 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 diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 3fc36dda..95a1fdad 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -885,6 +885,10 @@ subroutine temperature_changes (dt, & call icepack_warnings_add(warnstr) call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + write(warnstr,*) subname, 'enum:', e_num + call icepack_warnings_add(warnstr) + call ESMF_LogWrite(warnstr, ESMF_LOGMSG_INFO) + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//" temperature_changes: Thermo iteration does not converge" ) return diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 65057a60..76dff861 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -88,7 +88,7 @@ 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 = 800.0, cold_temp_flag = c0 - 60.0 + 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, & From 0a4650f63aafeb4aa16739bc78d1f5fe58403781 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 17 Sep 2024 12:11:48 +1000 Subject: [PATCH 12/19] 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 e3ed16df..7a25bae0 100644 --- a/columnphysics/icepack_flux.F90 +++ b/columnphysics/icepack_flux.F90 @@ -299,7 +299,7 @@ subroutine set_sfcflux (aicen, & character(len=*),parameter :: subname='(set_sfcflux)' - raicen = c1 !/ aicen + raicen = c1 / aicen #ifdef CICE_IN_NEMO !---------------------------------------------------------------------- From 9d7ef628d3e4289c6fb38ea1cd1e41f4400e3b69 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 17 Sep 2024 12:31:55 +1000 Subject: [PATCH 13/19] more agressively zap small ice areas --- columnphysics/icepack_itd.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index 669dcced..940e03b3 100644 --- a/columnphysics/icepack_itd.F90 +++ b/columnphysics/icepack_itd.F90 @@ -1080,7 +1080,7 @@ subroutine zap_small_areas (dt, & call icepack_warnings_add(subname//' Zap ice: negative ice area') return elseif (abs(aicen(n)) /= c0 .and. & - abs(aicen(n)) <= puny) then + abs(aicen(n)) <= 1.0e-2_dbl_kind) then !----------------------------------------------------------------- ! Account for tracers important for conservation From d5296fe42b2b91ec83d1ac3a676929e3f947a979 Mon Sep 17 00:00:00 2001 From: Kieran Ricardo Date: Tue, 15 Oct 2024 11:42:26 +1100 Subject: [PATCH 14/19] revert zapping back to original values - possible water conservation isse --- columnphysics/icepack_itd.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index 940e03b3..669dcced 100644 --- a/columnphysics/icepack_itd.F90 +++ b/columnphysics/icepack_itd.F90 @@ -1080,7 +1080,7 @@ subroutine zap_small_areas (dt, & call icepack_warnings_add(subname//' Zap ice: negative ice area') return elseif (abs(aicen(n)) /= c0 .and. & - abs(aicen(n)) <= 1.0e-2_dbl_kind) then + abs(aicen(n)) <= puny) then !----------------------------------------------------------------- ! Account for tracers important for conservation From 5cb072df687cbaf0631520e0dc810aea2b88318e Mon Sep 17 00:00:00 2001 From: Spencer Wong Date: Thu, 14 Nov 2024 22:33:05 +1100 Subject: [PATCH 15/19] Swap hi_min to 0.2 and add aicenmin=1.e-5 --- columnphysics/icepack_itd.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index 669dcced..242151c7 100644 --- a/columnphysics/icepack_itd.F90 +++ b/columnphysics/icepack_itd.F90 @@ -1053,6 +1053,7 @@ subroutine zap_small_areas (dt, & first_ice ! For bgc tracers. Set to true if zapping ice ! local variables + real (kind=dbl_kind) :: aicenmin = 1.0e-5 ! Cutoff concentration for zapping integer (kind=int_kind) :: & n, k, it, & !counting indices @@ -1080,7 +1081,7 @@ subroutine zap_small_areas (dt, & call icepack_warnings_add(subname//' Zap ice: negative ice area') return elseif (abs(aicen(n)) /= c0 .and. & - abs(aicen(n)) <= puny) then + abs(aicen(n)) <= aicenmin) then !----------------------------------------------------------------- ! Account for tracers important for conservation From 0cfa7fdc83069746d5651612340737d7c7d5ee8f Mon Sep 17 00:00:00 2001 From: Spencer Wong Date: Mon, 31 Mar 2025 22:15:46 +1100 Subject: [PATCH 16/19] 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 95a1fdad..4ff835eb 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -265,8 +265,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 31305ce7ebd30688bfb5d6929bc09a31f89e6c67 Mon Sep 17 00:00:00 2001 From: Spencer Wong Date: Wed, 9 Apr 2025 15:11:37 +1000 Subject: [PATCH 17/19] 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 1ebf3322..84ef0fd2 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -131,7 +131,7 @@ module icepack_parameters dSin0_frazil = c3 ,&! bulk salinity reduction of newly formed frazil dts_b = 50._dbl_kind ,&! zsalinity timestep 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 05e415da0c17697054f51279d01aa62cb7a32521 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 18/19] 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 4ff835eb..d64e7e91 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -553,7 +553,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 76dff861..c7725abc 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -1319,11 +1319,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 @@ -1576,7 +1574,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 33bbca5a61e9449f77f451e8809f27e2304734f4 Mon Sep 17 00:00:00 2001 From: Spencer Wong Date: Thu, 24 Apr 2025 10:49:21 +1000 Subject: [PATCH 19/19] 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 242151c7..5343341c 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 d64e7e91..22bdf455 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -709,9 +709,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) )