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_itd.F90 b/columnphysics/icepack_itd.F90 index 013373ad..5343341c 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 diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index ed76b6b4..84ef0fd2 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 @@ -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 diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 30c247ac..22bdf455 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,54 @@ 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 - - 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 +422,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 +431,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 +553,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 @@ -578,6 +611,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)) @@ -666,10 +709,7 @@ 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 & + ferr = abs( (enew-einit+e_num)/dt & - (fcondtopn - fcondbot + fswint) ) ! factor of 0.9 allows for roundoff errors later @@ -699,76 +739,150 @@ 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) + 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) + + 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 e931629c..c7725abc 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 @@ -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,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 @@ -1511,8 +1573,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 @@ -1934,7 +1995,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 +2008,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 +2021,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 +2038,51 @@ subroutine conservation_check_vthermo(dt, & einp = (fsurfn - flatn + fswint - fhocnn & - fsnow*Lfresh - fadvocn) * dt ferr = abs(efinal-einit-einp) / dt - - if (ferr > 1.1_dbl_kind*ferrmax) then + + ! 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 + 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 + 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, '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 + 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) + 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 + 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 +2103,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 @@ -2687,7 +2794,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 @@ -2697,8 +2803,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))