diff --git a/columnphysics/icepack_flux.F90 b/columnphysics/icepack_flux.F90 index d5578076..082f6eff 100644 --- a/columnphysics/icepack_flux.F90 +++ b/columnphysics/icepack_flux.F90 @@ -351,7 +351,7 @@ subroutine set_sfcflux (aicen, & raicen = c1 -#ifdef CICE_IN_NEMO +#if defined(CICE_IN_NEMO) || defined(CICE_ACCESS3) !---------------------------------------------------------------------- ! Convert fluxes from GBM values to per ice area values when ! running in NEMO environment. (When in standalone mode, fluxes diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index fe6ce051..d75909d1 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -136,7 +136,9 @@ module icepack_parameters aspect_rapid_mode = 1.0_dbl_kind,&! aspect ratio (larger is wider) dSdt_slow_mode = -1.5e-7_dbl_kind,&! slow mode drainage strength (m s-1 K-1) phi_c_slow_mode = 0.05_dbl_kind,&! critical liquid fraction porosity cutoff - phi_i_mushy = 0.85_dbl_kind ! liquid fraction of congelation ice + phi_i_mushy = 0.85_dbl_kind,&! liquid fraction of congelation ice + ratio_Wm2_m = 1000.0_dbl_kind,&! max condutive flux/depth ratio (W m) + cold_temp_flag = -60.0_dbl_kind ! min temp used to limit the conductive flux (C) integer (kind=int_kind), public :: & ktherm = 1 ! type of thermodynamics @@ -585,8 +587,8 @@ subroutine icepack_init_parameters( & cpl_frazil_in, semi_implicit_Tsfc_in, vapor_flux_correction_in, & Rac_rapid_mode_in, aspect_rapid_mode_in, & dSdt_slow_mode_in, phi_c_slow_mode_in, & - phi_i_mushy_in, shortwave_in, albedo_type_in, albsnowi_in, & - albicev_in, albicei_in, albsnowv_in, & + phi_i_mushy_in, ratio_Wm2_m_in, cold_temp_flag_in, shortwave_in, albedo_type_in, & + albsnowi_in, albicev_in, albicei_in, albsnowv_in, & ahmax_in, R_ice_in, R_pnd_in, R_snw_in, dT_mlt_in, rsnw_mlt_in, & kalg_in, R_gC2molC_in, kstrength_in, krdg_partic_in, krdg_redist_in, mu_rdg_in, & atmbndy_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, & @@ -732,7 +734,9 @@ subroutine icepack_init_parameters( & aspect_rapid_mode_in , & ! aspect ratio for rapid drainage mode (larger=wider) dSdt_slow_mode_in , & ! slow mode drainage strength (m s-1 K-1) phi_c_slow_mode_in , & ! liquid fraction porosity cutoff for slow mode - phi_i_mushy_in ! liquid fraction of congelation ice + phi_i_mushy_in , & ! liquid fraction of congelation ice + ratio_Wm2_m_in , & ! max condutive flux/depth ratio (W m) + cold_temp_flag_in ! min temp used to limit the conductive flux (C) character(len=*), intent(in), optional :: & congel_freeze_in ! congelation computation @@ -1205,6 +1209,8 @@ subroutine icepack_init_parameters( & if (present(dSdt_slow_mode_in) ) dSdt_slow_mode = dSdt_slow_mode_in if (present(phi_c_slow_mode_in) ) phi_c_slow_mode = phi_c_slow_mode_in if (present(phi_i_mushy_in) ) phi_i_mushy = phi_i_mushy_in + if (present(ratio_Wm2_m_in) ) ratio_Wm2_m = ratio_Wm2_m_in + if (present(cold_temp_flag_in) ) cold_temp_flag = cold_temp_flag_in if (present(shortwave_in) ) shortwave = shortwave_in if (present(albedo_type_in) ) albedo_type = albedo_type_in if (present(albicev_in) ) albicev = albicev_in @@ -1598,8 +1604,9 @@ subroutine icepack_query_parameters( & Lfresh_out, cprho_out, Cp_out, ustar_min_out, hi_min_out, a_rapid_mode_out, & ktherm_out, conduct_out, fbot_xfer_type_out, calc_Tsfc_out, & Rac_rapid_mode_out, aspect_rapid_mode_out, dSdt_slow_mode_out, & - phi_c_slow_mode_out, phi_i_mushy_out, shortwave_out, semi_implicit_Tsfc_out, & - albedo_type_out, albicev_out, albicei_out, albsnowv_out, vapor_flux_correction_out, & + phi_c_slow_mode_out, phi_i_mushy_out, ratio_Wm2_m_out, cold_temp_flag_out, & + shortwave_out, semi_implicit_Tsfc_out, albedo_type_out, albicev_out, & + albicei_out, albsnowv_out, vapor_flux_correction_out, & albsnowi_out, ahmax_out, R_ice_out, R_pnd_out, R_snw_out, dT_mlt_out, & rsnw_mlt_out, dEdd_algae_out, & kalg_out, R_gC2molC_out, kstrength_out, krdg_partic_out, krdg_redist_out, mu_rdg_out, & @@ -1754,7 +1761,9 @@ subroutine icepack_query_parameters( & aspect_rapid_mode_out , & ! aspect ratio for rapid drainage mode (larger=wider) dSdt_slow_mode_out , & ! slow mode drainage strength (m s-1 K-1) phi_c_slow_mode_out , & ! liquid fraction porosity cutoff for slow mode - phi_i_mushy_out ! liquid fraction of congelation ice + phi_i_mushy_out , & ! liquid fraction of congelation ice + ratio_Wm2_m_out , & ! max condutive flux/depth ratio (W m) + cold_temp_flag_out ! min temp used to limit the conductive flux (C) character(len=*), intent(out), optional :: & congel_freeze_out ! congelation computation @@ -2261,6 +2270,8 @@ subroutine icepack_query_parameters( & if (present(dSdt_slow_mode_out) ) dSdt_slow_mode_out = dSdt_slow_mode if (present(phi_c_slow_mode_out) ) phi_c_slow_mode_out = phi_c_slow_mode if (present(phi_i_mushy_out) ) phi_i_mushy_out = phi_i_mushy + if (present(ratio_Wm2_m_out) ) ratio_Wm2_m_out = ratio_Wm2_m + if (present(cold_temp_flag_out) ) cold_temp_flag_out = cold_temp_flag if (present(shortwave_out) ) shortwave_out = shortwave if (present(albedo_type_out) ) albedo_type_out = albedo_type if (present(albicev_out) ) albicev_out = albicev @@ -2575,6 +2586,8 @@ subroutine icepack_write_parameters(iounit) write(iounit,*) " dSdt_slow_mode = ", dSdt_slow_mode write(iounit,*) " phi_c_slow_mode = ", phi_c_slow_mode write(iounit,*) " phi_i_mushy= ", phi_i_mushy + write(iounit,*) " ratio_Wm2_m= ", ratio_Wm2_m + write(iounit,*) " cold_temp_flag = ", cold_temp_flag write(iounit,*) " shortwave = ", trim(shortwave) write(iounit,*) " albedo_type= ", trim(albedo_type) write(iounit,*) " albicev = ", albicev diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index f9485e32..7dc0e309 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -69,7 +69,7 @@ subroutine temperature_changes (dt, & fsensn, flatn, & flwoutn, fsurfn, & fcondtopn,fcondbot, & - einit ) + einit, einex_sfc_flux) real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -122,6 +122,9 @@ subroutine temperature_changes (dt, & zqsn , & ! snow layer enthalpy (J m-3) zTsn ! internal snow layer temperatures + real (kind=dbl_kind), intent(out):: & + einex_sfc_flux ! excess energy from conductive flux (J m-2) + ! local variables integer (kind=int_kind), parameter :: & @@ -149,14 +152,14 @@ subroutine temperature_changes (dt, & enew ! new energy of melting after temp change (J m-2) real (kind=dbl_kind) :: & - dTsf_prev , & ! dTsf from previous iteration - dTi1_prev , & ! dTi1 from previous iteration - dfsens_dT , & ! deriv of fsens wrt Tsf (W m-2 deg-1) - dflat_dT , & ! deriv of flat wrt Tsf (W m-2 deg-1) - dflwout_dT , & ! deriv of flwout wrt Tsf (W m-2 deg-1) - dt_rhoi_hlyr, & ! dt/(rhoi*hilyr) - einex , & ! excess energy from dqmat to ocean - ferr ! energy conservation error (W m-2) + dTsf_prev , & ! dTsf from previous iteration + dTi1_prev , & ! dTi1 from previous iteration + dfsens_dT , & ! deriv of fsens wrt Tsf (W m-2 deg-1) + dflat_dT , & ! deriv of flat wrt Tsf (W m-2 deg-1) + dflwout_dT , & ! deriv of flwout wrt Tsf (W m-2 deg-1) + dt_rhoi_hlyr , & ! dt/(rhoi*hilyr) + einex_sfc_calc, & ! excess energy from dqmat to ocean (J m-2) + ferr ! energy conservation error (W m-2) real (kind=dbl_kind), dimension (nilyr) :: & Tin_init , & ! zTin at beginning of time step @@ -184,15 +187,19 @@ subroutine temperature_changes (dt, & kh ! effective conductivity at interfaces (W m-2 deg-1) real (kind=dbl_kind) :: & - ci , & ! specific heat of sea ice (J kg-1 deg-1) - avg_Tsf , & ! = 1. if Tsf averaged w/Tsf_start, else = 0. - Iswabs_tmp , & ! energy to melt through fraction frac of layer - Sswabs_tmp , & ! same for snow - dswabs , & ! difference in swabs and swabs_tmp - frac + ci , & ! specific heat of sea ice (J kg-1 deg-1) + avg_Tsf , & ! = 1. if Tsf averaged w/Tsf_start, else = 0. + Iswabs_tmp , & ! energy to melt through fraction frac of layer + Sswabs_tmp , & ! same for snow + dswabs , & ! difference in swabs and swabs_tmp + frac , & + fcondtopn_reduction, & ! reduction in downward cond flux at top surface (W m-2) + fcondtopn_force , & ! reduced downward cond flux at top surface (W m-2) + dqmat_sn ! associated enthalpy difference at top snow layer (J m-3) logical (kind=log_kind) :: & - converged ! = true when local solution has converged + converged , & ! = true when local solution has converged + Top_T_was_reset_last_time ! = true when surface temperature was reset in prev iteration logical (kind=log_kind) , dimension (nilyr) :: & reduce_kh ! reduce conductivity when T exceeds Tmlt @@ -202,7 +209,9 @@ subroutine temperature_changes (dt, & !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- - + fcondtopn_reduction = c0 + einex_sfc_flux = c0 + Top_T_was_reset_last_time = .false. converged = .false. l_snow = .false. l_cold = .true. @@ -212,7 +221,7 @@ subroutine temperature_changes (dt, & dfsens_dT = c0 dflat_dT = c0 dflwout_dT = c0 - einex = c0 + einex_sfc_calc = c0 if (semi_implicit_Tsfc) then ! initialize dfsurf_dT = dfsurfdTs_cpl dflat_dT = dflatdTs_cpl @@ -337,7 +346,7 @@ subroutine temperature_changes (dt, & if (.not.semi_implicit_Tsfc) dfsurf_dT = c0 avg_Tsi = c0 enew = c0 - einex = c0 + einex_sfc_calc = c0 !----------------------------------------------------------------- ! Update specific heat of ice layers. @@ -428,6 +437,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 +446,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 +568,33 @@ 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 + + if (.not. calc_Tsfc) then + ! return this energy to the ocean + + dqmat_sn = (zTsn(k)*cp_ice - Lfresh)*rhos - zqsn(k) + + ! 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. + einex_sfc_flux = einex_sfc_flux + hslyr * dqmat_sn + else + Top_T_was_reset_last_time = .true. + endif + endif + end if + + 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. calc_Tsfc) .and. (.not. l_snow) .and. (k == 1)) then + if (Top_T_was_reset_last_time) then + fcondtopn_reduction = fcondtopn_reduction + dqmat(k)*hilyr / dt + Top_T_was_reset_last_time = .false. + einex_sfc_flux = einex_sfc_flux + 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)) @@ -637,7 +683,7 @@ subroutine temperature_changes (dt, & zqin(k) = -rhoi * (-cp_ice*zTin(k) + Lfresh) endif enew = enew + hilyr * zqin(k) - einex = einex + hilyr * dqmat(k) + einex_sfc_calc = einex_sfc_calc + hilyr * dqmat(k) Tin_start(k) = zTin(k) ! for next iteration @@ -683,11 +729,15 @@ subroutine temperature_changes (dt, & fcondbot = kh(1+nslyr+nilyr) * & (zTin(nilyr) - Tbot) - ! Flux extra energy out of the ice - fcondbot = fcondbot + einex/dt - - ferr = abs( (enew-einit)/dt & + if (calc_Tsfc) then + ! Flux extra energy out of the ice + fcondbot = fcondbot + einex_sfc_calc/dt + ferr = abs( (enew-einit)/dt & + - (fcondtopn - fcondbot + fswint) ) + else + ferr = abs( (enew-einit+einex_sfc_flux)/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) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 1a2be844..15f6ac71 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -23,10 +23,11 @@ module icepack_therm_vertical use icepack_fsd, only: floe_rad_c, floe_binwidth - use icepack_parameters, only: c0, c1, c2, p001, p5, puny + use icepack_parameters, only: c0, c1, c2, p001, p5, puny, c100 use icepack_parameters, only: pi, depressT, Lvap, hs_min, cp_ice, min_salin use icepack_parameters, only: cp_ocn, rhow, rhoi, rhos, Lfresh, rhofresh, ice_ref_salinity use icepack_parameters, only: ktherm, calc_Tsfc, rsnw_fall, rsnw_tmax + use icepack_parameters, only: ratio_Wm2_m, cold_temp_flag use icepack_parameters, only: ustar_min, fbot_xfer_type, formdrag, calc_strair use icepack_parameters, only: rfracmin, rfracmax, dpscale, frzpnd, snwgrain, snwlvlfac use icepack_parameters, only: phi_i_mushy, floeshape, floediam, use_smliq_pnd, snwredist @@ -75,6 +76,48 @@ module icepack_therm_vertical contains !======================================================================= + +! Function for limiting the conductive flux supplied from an atmosphere model when calc_tsfc=.false. + + function cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) result(fcondtopn_solve) + + integer (kind=int_kind), intent(in) :: & + nilyr , & ! number of ice layers + nslyr ! number of snow layers + real (kind=dbl_kind), intent(in) :: fcondtopn ! downward cond flux at top surface (W m-2) + real (kind=dbl_kind), intent(in) :: hin ! ice thickness (m) + real (kind=dbl_kind), intent(in) :: zTin(nilyr) ! internal ice layer temperatures (C) + real (kind=dbl_kind), intent(in) :: zTsn(nslyr) ! internal snow layer temperatures (C) + real (kind=dbl_kind), intent(in) :: hslyr ! snow layer thickness (m) + + real (kind=dbl_kind) :: fcondtopn_solve ! limited downward cond flux at top surface (W m-2) + real (kind=dbl_kind) :: & + top_layer_temp, & ! top layer temperature (C) + reduce_ratio , & ! reduction ratio of downward cond flux at top surface + reduce_amount ! reduction in downward cond flux at top surface (W m-2) + + if (abs(fcondtopn) > ratio_Wm2_m * hin) then + fcondtopn_solve = sign(ratio_Wm2_m * hin,fcondtopn) + else + fcondtopn_solve = fcondtopn + endif + + 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) / (c100 + 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 +289,11 @@ subroutine thermo_vertical (dt, aicen, & real (kind=dbl_kind) :: & fadvocn, saltvol, dfsalt ! advective heat flux to ocean + real (kind=dbl_kind) :: & + fcondtopn_solve, & ! limited downward cond flux at top surface (W m-2) + fcondtopn_extra, & ! excess downward cond flux at top surface (W m-2) + einex_sfc_flux ! excess energy removed during temperature_changes (J m-2) + character(len=*),parameter :: subname='(thermo_vertical)' !----------------------------------------------------------------- @@ -272,6 +320,8 @@ subroutine thermo_vertical (dt, aicen, & meltsliq= c0 massice(:) = c0 massliq(:) = c0 + einex_sfc_flux = 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 !----------------------------------------------------------------- @@ -340,6 +392,13 @@ subroutine thermo_vertical (dt, aicen, & else ! ktherm + if (calc_Tsfc) then + fcondtopn_solve = fcondtopn + else + fcondtopn_solve = cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) + end if + fcondtopn_extra = fcondtopn - fcondtopn_solve + call temperature_changes(dt, & rhoa, flw, & potT, Qa, & @@ -353,10 +412,12 @@ subroutine thermo_vertical (dt, aicen, & Tsf, Tbot, & fsensn, flatn, & flwoutn, fsurfn, & - fcondtopn, fcondbotn, & - einit ) + fcondtopn_solve, fcondbotn, & + einit, einex_sfc_flux) if (icepack_warnings_aborted(subname)) return + fcondtopn = fcondtopn_solve + fcondtopn_extra + endif ! ktherm ! mass of ice and liquid water in snow @@ -411,7 +472,8 @@ subroutine thermo_vertical (dt, aicen, & mlt_onset, frz_onset, & zSin, sss, & sst, & - dsnow, rsnw) + dsnow, rsnw, & + einex_sfc_flux, fcondtopn_extra ) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- @@ -425,7 +487,7 @@ subroutine thermo_vertical (dt, aicen, & fsnow, einit, & einter, efinal, & fcondtopn, fcondbotn, & - fadvocn, fbot ) + fadvocn, fbot, einex_sfc_flux, fcondtopn_extra ) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- @@ -1070,7 +1132,8 @@ subroutine thickness_changes (dt, yday, & mlt_onset, frz_onset,& zSin, sss, & sst, & - dsnow, rsnw) + dsnow, rsnw, & + einex_sfc_flux, fcondtopn_extra) real (kind=dbl_kind), intent(in) :: & dt , & ! time step @@ -1138,6 +1201,9 @@ subroutine thickness_changes (dt, yday, & sst , & ! sea surface temperature (C) sss ! ocean salinity (PSU) + real (kind=dbl_kind), intent(in) :: & + einex_sfc_flux, & ! excess energy removed during temperature_changes (J m-2) + fcondtopn_extra ! excess downward cond flux at top surface (W m-2) ! local variables integer (kind=int_kind) :: & @@ -1271,7 +1337,7 @@ subroutine thickness_changes (dt, yday, & wk1 = (fsurfn - fcondtopn) * dt etop_mlt = max(wk1, c0) ! etop_mlt > 0 - wk1 = (fcondbotn - fbot) * dt + wk1 = (fcondbotn - fbot + fcondtopn_extra) * dt ebot_mlt = max(wk1, c0) ! ebot_mlt > 0 ebot_gro = min(wk1, c0) ! ebot_gro < 0 @@ -1561,8 +1627,7 @@ subroutine thickness_changes (dt, yday, & ! fhocn is the available ocean heat that is left after use by ice !----------------------------------------------------------------- - fhocnn = fbot & - + (esub + etop_mlt + ebot_mlt)/dt + fhocnn = fbot + (esub + etop_mlt + ebot_mlt + einex_sfc_flux)/dt !----------------------------------------------------------------- ! Add new snowfall at top surface @@ -1988,20 +2053,22 @@ subroutine conservation_check_vthermo(dt, & einit, einter, & efinal, & fcondtopn,fcondbotn, & - fadvocn, fbot ) + fadvocn, fbot, einex_sfc_flux, fcondtopn_extra) real (kind=dbl_kind), intent(in) :: & dt ! time step real (kind=dbl_kind), intent(in) :: & - fsurfn , & ! net flux to top surface, excluding fcondtopn - flatn , & ! surface downward latent heat (W m-2) - fhocnn , & ! fbot, corrected for any surplus energy - fswint , & ! SW absorbed in ice interior, below surface (W m-2) - fsnow , & ! snowfall rate (kg m-2 s-1) - fcondtopn , & - fadvocn , & - fbot + fsurfn , & ! net flux to top surface, excluding fcondtopn + flatn , & ! surface downward latent heat (W m-2) + fhocnn , & ! fbot, corrected for any surplus energy + fswint , & ! SW absorbed in ice interior, below surface (W m-2) + fsnow , & ! snowfall rate (kg m-2 s-1) + fcondtopn , & ! downward cond flux at top surface (W m-2) + fadvocn , & ! advective heat flux to ocean (W m-2) + fbot , & ! ice-ocean heat flux at bottom surface (W/m^2) + einex_sfc_flux, & ! excess energy removed during temperature_changes (J m-2) + fcondtopn_extra ! excess downward cond flux at top surface (W m-2) real (kind=dbl_kind), intent(in) :: & einit , & ! initial energy of melting (J m-2) @@ -2054,6 +2121,12 @@ subroutine conservation_check_vthermo(dt, & call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Input energy =', einp call icepack_warnings_add(warnstr) + if (.not. calc_Tsfc) then + write(warnstr,*) subname, 'Numerical energy =', einex_sfc_flux + call icepack_warnings_add(warnstr) + write(warnstr,*) subname, 'fcondtopn_extra energy =', fcondtopn_extra + call icepack_warnings_add(warnstr) + end if write(warnstr,*) subname, 'fbot,fcondbot:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, fbot,fcondbotn diff --git a/configuration/scripts/machines/Macros.gadi_oneapi b/configuration/scripts/machines/Macros.gadi_oneapi new file mode 100644 index 00000000..eacc6cb0 --- /dev/null +++ b/configuration/scripts/machines/Macros.gadi_oneapi @@ -0,0 +1,44 @@ +#============================================================================== +# Makefile macros for NCI Gadi, intel compiler +#============================================================================== + +CPP := fpp +CPPDEFS := -DFORTRANUNDERSCORE -DREPRODUCIBLE ${ICE_CPPDEFS} +CFLAGS := -c -O2 -fp-model precise -Wno-unused-variable -Wno-unused-parameter + +FIXEDFLAGS := -132 +FREEFLAGS := -FR + +NCI_INTEL_FLAGS := -r8 -i4 -traceback -w -fpe0 -ftz -convert big_endian -assume byterecl -check noarg_temp_created +NCI_REPRO_FLAGS := -fp-model precise -fp-model source -align all + +ifeq ($(ICE_BLDDEBUG), true) + NCI_DEBUG_FLAGS := -g3 -O0 -debug all -check all -no-vec -assume nobuffered_io + FFLAGS := $(NCI_INTEL_FLAGS) $(NCI_REPRO_FLAGS) $(NCI_DEBUG_FLAGS) + CPPDEFS := $(CPPDEFS) -DDEBUG=$(DEBUG) +else + NCI_OPTIM_FLAGS := -g3 -O2 + FFLAGS := $(NCI_INTEL_FLAGS) $(NCI_REPRO_FLAGS) $(NCI_OPTIM_FLAGS) +endif + +SCC := icx +SFC := ifx +MPICC := mpicc +MPIFC := mpifort + +ifeq ($(ICE_COMMDIR), mpi) + FC := $(MPIFC) + CC := $(MPICC) +else + FC := $(SFC) + CC := $(SCC) +endif +LD:= $(FC) + +SLIBS := $(SLIBS) +INCLDIR := $(INCLDIR) + +# use system netcdf build +SLIBS += -L$(NETCDF)/lib -lnetcdf -lnetcdff +INCLDIR += -I$(NETCDF)/include + diff --git a/configuration/scripts/machines/env.gadi_oneapi b/configuration/scripts/machines/env.gadi_oneapi new file mode 100644 index 00000000..b3c387ce --- /dev/null +++ b/configuration/scripts/machines/env.gadi_oneapi @@ -0,0 +1,33 @@ +#!/bin/csh -f + +set inp = "undefined" +if ($#argv == 1) then + set inp = $1 +endif + +if ("$inp" != "-nomodules") then + + source /etc/profile.d/modules.csh + + module load intel-compiler-llvm + module load netcdf + module load openmpi + +endif + +setenv ICE_MACHINE_MACHNAME gadi +setenv ICE_MACHINE_MACHINFO "Intel Xeon Scalable" +setenv ICE_MACHINE_ENVNAME intel +setenv ICE_MACHINE_ENVINFO INTEL_COMPILER_VERSION=$INTEL_COMPILER_LLVM_BASE +setenv ICE_MACHINE_MAKE gmake +setenv ICE_MACHINE_WKDIR /scratch/$PROJECT/$USER/ICEPACK_RUNS +setenv ICE_MACHINE_INPUTDATA /g/data/ik11/inputs/CICE_data/icepack-dirs/input +setenv ICE_MACHINE_BASELINE /scratch/$PROJECT/$USER/ICEPACK_BASELINE +setenv ICE_MACHINE_SUBMIT "qsub" +setenv ICE_MACHINE_PROJ $PROJECT +setenv ICE_MACHINE_ACCT $USER +setenv ICE_MACHINE_QUEUE "normal" +setenv ICE_MACHINE_TPNODE 48 +setenv ICE_MACHINE_BLDTHRDS 4 +setenv ICE_MACHINE_QSTAT "qstat" +setenv ICE_CPPDEFS -DUSE_NETCDF \ No newline at end of file diff --git a/doc/source/science_guide/sg_boundary_forcing.rst b/doc/source/science_guide/sg_boundary_forcing.rst index e73a13a9..c4cb7e30 100755 --- a/doc/source/science_guide/sg_boundary_forcing.rst +++ b/doc/source/science_guide/sg_boundary_forcing.rst @@ -81,9 +81,9 @@ subroutine, *thermo_vertical*. At the end of the time step, the surface temperature and effective conductivity (i.e., thermal conductivity divided by thickness) of the top ice/snow layer in each category are returned to the atmosphere model via the coupler. Since the ice surface -temperature is treated explicitly, the effective conductivity may need -to be limited to ensure stability. As a result, accuracy may be -significantly reduced, especially for thin ice or snow layers. A more +temperature is treated explicitly, the conductive flux is limited to ensure stability, +with any excess energy either used for bottom melting or transferred to the ocean as a heat flux. +As a result, accuracy may be significantly reduced, especially for thin ice or snow layers. A more stable and accurate procedure would be to compute the temperature profiles for both the atmosphere and ice, together with the surface fluxes, in a single implicit calculation. This was judged impractical, @@ -97,6 +97,8 @@ implicitly. The resultant surface temperature change is passed back to the atmosphere model via coupler to complete the full update of its temperature profiles. This middle-ground approach, enabled by ``vapor_flux_correction=true``, does not sacrifice accuracy because it does not need effective conductivity limiting as in the explicit case. +The GEOS approach is available for either mushy or BL99 thermodynamics, while the +flux-limiting approaches are only available for BL99. diff --git a/doc/source/science_guide/sg_thermo.rst b/doc/source/science_guide/sg_thermo.rst index bf3d3c25..36355312 100755 --- a/doc/source/science_guide/sg_thermo.rst +++ b/doc/source/science_guide/sg_thermo.rst @@ -1375,9 +1375,15 @@ is that the solution scheme is no longer unconditionally stable. Instead, the effective conductivity in the top layer must satisfy a diffusive CFL condition: -.. math:: +.. math:: K^* \le {\rho ch \over \Delta t}. +To improve stability when ``calc_Tsfc=false``, the conductive flux is limited before +the thermodynamic solve. Any excess flux is then used for bottom melt. The +conductive flux may also be limited during the thermodynamic solve if the +surface temperature exceeds the melting point, in which case the excess flux +is passed to the ocean. + For thin layers and typical coupling intervals (:math:`\sim 1` hr), :math:`K^*` may need to be limited before being passed to the atmosphere via the coupler. Otherwise, the fluxes that are returned to the host sea ice model may