Add conductive flux limiting for vertical thermodynamics when calc_Tsfc=.false.
#557
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -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 | ||
| ratio_Wm2_m = 1000.0,&! max condutive flux/depth ratio (W m) | ||
| cold_temp_flag = c0 - 60.0 ! 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, 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 | ||
| 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 | ||
|
Comment on lines
+1212
to
+1213
Contributor
There was a problem hiding this comment. @kieranricardo has added these as configurable - but we don't plan to add them to the cice namelist (as we dont foresee changing them). Are we happy to change answers when calc_Tsfc = .false. - @apcraig ? If so, should we not make these two new parameters configurable through
Contributor
There was a problem hiding this comment. Are there values of these parameters that would preserve BFBness when calc_Tsfc = .false.? If so, I recommend setting the code defaults to those values and allow them to be changed via namelist.
Contributor
Author
There was a problem hiding this comment. We could pick extreme value to ensure the |
||
| 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, 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 | ||
| 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 | ||
|
Contributor
There was a problem hiding this comment. should these be within a
Contributor
Author
There was a problem hiding this comment. It looks like there's no other
Contributor
There was a problem hiding this comment. this is okay as long as they are not used elsewhere except when calc_Tsfc = .false. |
||
| write(iounit,*) " shortwave = ", trim(shortwave) | ||
| write(iounit,*) " albedo_type= ", trim(albedo_type) | ||
| write(iounit,*) " albicev = ", albicev | ||
|
|
||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -69,7 +69,7 @@ subroutine temperature_changes (dt, & | |
| fsensn, flatn, & | ||
| flwoutn, fsurfn, & | ||
| fcondtopn,fcondbot, & | ||
| 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 | ||
|
kieranricardo marked this conversation as resolved.
Outdated
|
||
| ! 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 , & | ||
| fcondtopn_reduction, & | ||
| fcondtopn_force, dqmat_sn | ||
|
|
||
| logical (kind=log_kind) :: & | ||
| 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 | ||
|
|
||
| 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_force) | ||
| if (icepack_warnings_aborted(subname)) return | ||
|
|
||
| endif ! calc_Tsfc | ||
|
|
@@ -558,7 +565,33 @@ subroutine temperature_changes (dt, & | |
| else | ||
| zTsn(k) = c0 | ||
| endif | ||
| if ((l_brine) .and. zTsn(k)>c0) then | ||
|
|
||
| if (.not. calc_Tsfc) then | ||
| ! Alex West: return this energy to the ocean | ||
|
kieranricardo marked this conversation as resolved.
Outdated
|
||
|
|
||
| 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 | ||
| end if | ||
|
|
||
| zTsn(k) = min(zTsn(k), c0) | ||
|
|
||
| endif | ||
|
|
||
| !----------------------------------------------------------------- | ||
| ! If condition 1 or 2 failed, average new snow layer | ||
|
|
@@ -592,6 +625,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. | ||
| 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 +726,15 @@ subroutine temperature_changes (dt, & | |
| fcondbot = kh(1+nslyr+nilyr) * & | ||
| (zTin(nilyr) - Tbot) | ||
|
|
||
| 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 & | ||
|
Contributor
There was a problem hiding this comment. How are einex and e_num related? Could you modify einex instead, just above, so that this block of code does not change?
Contributor
Author
There was a problem hiding this comment.
The main difference is what's done with this excess energy. I'm actually not sure what the best approach is for
Contributor
There was a problem hiding this comment. Are you suggesting not using e_num at all and seeing if it's stable using the already available einex ? Give it a try I guess, it looks like einex should be calculated fine when calc_Tsfc is false?
Contributor
Author
There was a problem hiding this comment. I meant use the current approach for calculating In the code this would add extra checks around einex, calculating it normally if
Contributor
Author
There was a problem hiding this comment. After some discussion it doesn't seem like this approach will work. With
Contributor
There was a problem hiding this comment. Thanks for looking into this. Then perhaps a compromise will work. In the PR comments at the top:
Since the excess in 1. is used for bottom melt, that's essentially supplementing
Contributor
Author
There was a problem hiding this comment. Good point! But that first approach could be a little involved so I think the second approach sounds good. The first approach would involve adding an extra inout argument to
Contributor
Author
There was a problem hiding this comment. I've updated the variable name to |
||
| - (fcondtopn - fcondbot + fswint) ) | ||
| end if | ||
|
|
||
| ! factor of 0.9 allows for roundoff errors later | ||
| if (ferr > 0.9_dbl_kind*ferrmax) then ! condition (5) | ||
|
|
||
Uh oh!
There was an error while loading. Please reload this page.