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 !---------------------------------------------------------------------- diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index f9485e32..5809c06e 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, e_num) real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -122,6 +122,8 @@ subroutine temperature_changes (dt, & zqsn , & ! snow layer enthalpy (J m-3) zTsn ! internal snow layer temperatures + real (kind=dbl_kind), intent(out):: & + e_num ! local variables integer (kind=int_kind), parameter :: & @@ -189,20 +191,24 @@ subroutine temperature_changes (dt, & Iswabs_tmp , & ! energy to melt through fraction frac of layer Sswabs_tmp , & ! same for snow dswabs , & ! difference in swabs and swabs_tmp - frac + frac , & + fcondtopn_reduction, & + fcondtopn_force, dqmat_sn logical (kind=log_kind) :: & - converged ! = true when local solution has converged + converged, Top_T_was_reset_last_time ! = true when local solution has converged logical (kind=log_kind) , dimension (nilyr) :: & - reduce_kh ! reduce conductivity when T exceeds Tmlt + reduce_kh ! reduce conductivity when T exceeds Tmlt character(len=*),parameter :: subname='(temperature_changes)' !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- - + fcondtopn_reduction = c0 + e_num = c0 + Top_T_was_reset_last_time = .false. converged = .false. l_snow = .false. l_cold = .true. @@ -428,6 +434,7 @@ subroutine temperature_changes (dt, & else + fcondtopn_force = fcondtopn - fcondtopn_reduction call get_matrix_elements_know_Tsfc ( & l_snow, Tbot, & Tin_init, Tsn_init, & @@ -436,7 +443,7 @@ subroutine temperature_changes (dt, & etai, etas, & sbdiag, diag, & spdiag, rhs, & - fcondtopn) + fcondtopn_force) if (icepack_warnings_aborted(subname)) return endif ! calc_Tsfc @@ -558,7 +565,31 @@ subroutine temperature_changes (dt, & else zTsn(k) = c0 endif - if (l_brine) zTsn(k) = min(zTsn(k), c0) + if ((l_brine) .and. zTsn(k)>c0) then + + ! Alex West: return this energy to the ocean + + dqmat_sn = (zTsn(k)*cp_ice - Lfresh)*rhos - zqsn(k) + + ! Alex West: If this is the second time in succession that Tsn(1) has been + ! reset, tell the solver to reduce the forcing at the top, and + ! pass the difference to the array enum where it will eventually + ! go into the ocean + ! This is done to avoid an 'infinite loop' whereby temp continually evolves + ! to the same point above zero, is reset, ad infinitum + if (l_snow .AND. k == 1) then + if (Top_T_was_reset_last_time) then + fcondtopn_reduction = fcondtopn_reduction + dqmat_sn*hslyr / dt + Top_T_was_reset_last_time = .false. + e_num = e_num + hslyr * dqmat_sn + else + Top_T_was_reset_last_time = .true. + endif + endif + + zTsn(k) = min(zTsn(k), c0) + + endif !----------------------------------------------------------------- ! If condition 1 or 2 failed, average new snow layer @@ -592,6 +623,16 @@ subroutine temperature_changes (dt, & dTmat(k) = zTin(k) - Tmlts(k) dqmat(k) = rhoi * dTmat(k) & * (cp_ice - Lfresh * Tmlts(k)/zTin(k)**2) + + if ((.not. l_snow) .and. (k == 1)) then + if (Top_T_was_reset_last_time) then + fcondtopn_reduction = fcondtopn_reduction + dqmat(k)*hilyr / dt + Top_T_was_reset_last_time = .false. + e_num = e_num + hilyr * dqmat(k) + else + Top_T_was_reset_last_time = .true. + endif + endif ! use this for the case that Tmlt changes by an amount dTmlt=Tmltnew-Tmlt(k) ! + rhoi * dTmlt & ! * (cp_ocn - cp_ice + Lfresh/zTin(k)) @@ -683,11 +724,15 @@ subroutine temperature_changes (dt, & fcondbot = kh(1+nslyr+nilyr) * & (zTin(nilyr) - Tbot) - ! Flux extra energy out of the ice - fcondbot = fcondbot + einex/dt - - ferr = abs( (enew-einit)/dt & + if (calc_Tsfc) then + ! Flux extra energy out of the ice + fcondbot = fcondbot + einex/dt + ferr = abs( (enew-einit)/dt & + - (fcondtopn - fcondbot + fswint) ) + else + ferr = abs( (enew-einit+e_num)/dt & - (fcondtopn - fcondbot + fswint) ) + end if ! factor of 0.9 allows for roundoff errors later if (ferr > 0.9_dbl_kind*ferrmax) then ! condition (5) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 00b3756d..bb6287fb 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -75,6 +75,49 @@ module icepack_therm_vertical contains !======================================================================= + + function cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) result(fcondtopn_solve) + + integer (kind=int_kind), intent(in) :: & + nilyr , & ! number of ice layers + nslyr ! number of snow layers + real (kind=dbl_kind), intent(in) :: fcondtopn + real (kind=dbl_kind), intent(in) :: hin + real (kind=dbl_kind), intent(in) :: zTin(nilyr) + real (kind=dbl_kind), intent(in) :: zTsn(nslyr) + real (kind=dbl_kind), intent(in) :: hslyr + + real (kind=dbl_kind) :: fcondtopn_solve + + real (kind=dbl_kind), parameter :: ratio_Wm2_m = 1000.0, cold_temp_flag = c0 - 60.0 + + ! AEW: New variables for cold-ice flux capping + real (kind=dbl_kind) :: top_layer_temp, & + reduce_ratio, & + reduce_amount + + + if (abs(fcondtopn) > ratio_Wm2_m * hin) then + fcondtopn_solve = sign(ratio_Wm2_m * hin,fcondtopn) + + else + fcondtopn_solve = fcondtopn + endif + + if (hslyr>hs_min) then + top_layer_temp = zTsn(1) + else + top_layer_temp = zTin(1) + endif + + if ((top_layer_temp < cold_temp_flag) .and. (fcondtopn_solve < c0)) then + reduce_ratio = (cold_temp_flag - top_layer_temp) / (100.0 + cold_temp_flag) + reduce_amount = reduce_ratio * fcondtopn_solve + fcondtopn_solve = fcondtopn_solve - reduce_amount + endif + + end function cap_conductive_flux + ! ! Driver for updating ice and snow internal temperatures and ! computing thermodynamic growth rates and atmospheric fluxes. @@ -246,6 +289,9 @@ subroutine thermo_vertical (dt, aicen, & real (kind=dbl_kind) :: & fadvocn, saltvol, dfsalt ! advective heat flux to ocean + real (kind=dbl_kind) :: & + fcondtopn_solve, fcondtopn_extra, e_num + character(len=*),parameter :: subname='(thermo_vertical)' !----------------------------------------------------------------- @@ -272,6 +318,8 @@ subroutine thermo_vertical (dt, aicen, & meltsliq= c0 massice(:) = c0 massliq(:) = c0 + e_num = c0 + fcondtopn_extra = c0 if (tr_pond) then dpnd_flush = c0 dpnd_expon = c0 @@ -284,6 +332,8 @@ subroutine thermo_vertical (dt, aicen, & fcondtopn = c0 endif + fcondtopn_solve = fcondtopn + !----------------------------------------------------------------- ! Compute variables needed for vertical thermo calculation !----------------------------------------------------------------- @@ -339,6 +389,12 @@ subroutine thermo_vertical (dt, aicen, & if (icepack_warnings_aborted(subname)) return else ! ktherm + fcondtopn_solve = cap_conductive_flux(nilyr, nslyr, fcondtopn, hin, zTsn, zTin, hslyr) + fcondtopn_extra = fcondtopn - fcondtopn_solve + + ! if (calc_Tsfc) then + ! fcondtopn = fcondtopn_solve + ! end if call temperature_changes(dt, & rhoa, flw, & @@ -353,8 +409,8 @@ subroutine thermo_vertical (dt, aicen, & Tsf, Tbot, & fsensn, flatn, & flwoutn, fsurfn, & - fcondtopn, fcondbotn, & - einit ) + fcondtopn_solve, fcondbotn, & + einit, e_num) if (icepack_warnings_aborted(subname)) return endif ! ktherm @@ -411,7 +467,8 @@ subroutine thermo_vertical (dt, aicen, & mlt_onset, frz_onset, & zSin, sss, & sst, & - dsnow, rsnw) + dsnow, rsnw, & + e_num, fcondtopn_extra ) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- @@ -425,7 +482,7 @@ subroutine thermo_vertical (dt, aicen, & fsnow, einit, & einter, efinal, & fcondtopn, fcondbotn, & - fadvocn, fbot ) + fadvocn, fbot, e_num, fcondtopn_extra ) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- @@ -1070,7 +1127,8 @@ subroutine thickness_changes (dt, yday, & mlt_onset, frz_onset,& zSin, sss, & sst, & - dsnow, rsnw) + dsnow, rsnw, & + e_num, fcondtopn_extra) real (kind=dbl_kind), intent(in) :: & dt , & ! time step @@ -1138,6 +1196,8 @@ subroutine thickness_changes (dt, yday, & sst , & ! sea surface temperature (C) sss ! ocean salinity (PSU) + real (kind=dbl_kind), intent(in) :: & + e_num, fcondtopn_extra ! local variables integer (kind=int_kind) :: & @@ -1271,7 +1331,7 @@ subroutine thickness_changes (dt, yday, & wk1 = (fsurfn - fcondtopn) * dt etop_mlt = max(wk1, c0) ! etop_mlt > 0 - wk1 = (fcondbotn - fbot) * dt + wk1 = (fcondbotn - fbot + fcondtopn_extra) * dt ebot_mlt = max(wk1, c0) ! ebot_mlt > 0 ebot_gro = min(wk1, c0) ! ebot_gro < 0 @@ -1561,8 +1621,7 @@ subroutine thickness_changes (dt, yday, & ! fhocn is the available ocean heat that is left after use by ice !----------------------------------------------------------------- - fhocnn = fbot & - + (esub + etop_mlt + ebot_mlt)/dt + fhocnn = fbot + (esub + etop_mlt + ebot_mlt + e_num)/dt !----------------------------------------------------------------- ! Add new snowfall at top surface @@ -1987,7 +2046,7 @@ subroutine conservation_check_vthermo(dt, & einit, einter, & efinal, & fcondtopn,fcondbotn, & - fadvocn, fbot ) + fadvocn, fbot, e_num, fcondtopn_extra) real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -2000,7 +2059,7 @@ subroutine conservation_check_vthermo(dt, & fsnow , & ! snowfall rate (kg m-2 s-1) fcondtopn , & fadvocn , & - fbot + fbot, e_num, fcondtopn_extra real (kind=dbl_kind), intent(in) :: & einit , & ! initial energy of melting (J m-2) @@ -2053,6 +2112,10 @@ subroutine conservation_check_vthermo(dt, & call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Input energy =', einp call icepack_warnings_add(warnstr) + write(warnstr,*) subname, 'Numerical energy =', e_num + call icepack_warnings_add(warnstr) + write(warnstr,*) subname, 'fcondtopn_extra energy =', fcondtopn_extra + call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'fbot,fcondbot:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, fbot,fcondbotn 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, &