From 57923db8786ed2c928bdcad08ab28063da4a92e0 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Fri, 16 Jun 2023 11:15:18 -0800 Subject: [PATCH] Fix for sigma_i going to zero in ITD landfast ice. (#189) - We can't divide by zero, so don't let sigma_i be zero - Version coming from the CICE Consortium - Update Dupont citation --- docs/zotero.bib | 5 +- src/SIS_dyn_cgrid.F90 | 97 ++++++++++++++++++++------------------- src/SIS_open_boundary.F90 | 14 +++--- 3 files changed, 61 insertions(+), 55 deletions(-) diff --git a/docs/zotero.bib b/docs/zotero.bib index 30a5107d..0becb5c7 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -1,9 +1,10 @@ @article{dupont2021, author = "F. Dupont and D. Dumont and J.-F. Lemieux and E. Dumas-Lefebvre and A. Caya", - year = 2021, + year = 2022, title = "A probabilistic seabed-ice keel interaction model", journal = "The Cryosphere", - note = "in review" + volume = "16", + pages = "1963--1977" } @article{lemieux2015, diff --git a/src/SIS_dyn_cgrid.F90 b/src/SIS_dyn_cgrid.F90 index dc11f222..cea51ac1 100644 --- a/src/SIS_dyn_cgrid.F90 +++ b/src/SIS_dyn_cgrid.F90 @@ -1853,8 +1853,9 @@ end subroutine basal_stress_coeff_C !! a normal distribution with sigma_b = 2.5d0. An improvement would !! be to provide the distribution based on high resolution data. !! -!! Dupont, F. Dumont, D., Lemieux, J.F., Dumas-Lefebvre, E., Caya, A. -!! in prep. +!! Dupont, F., D. Dumont, J.F. Lemieux, E. Dumas-Lefebvre, A. Caya (2022). +!! A probabilistic seabed-ice keel interaction model, The Cryosphere, 16, +!! 1963-1977. !! !! authors: D. Dumont, J.F. Lemieux, E. Dumas-Lefebvre, F. Dupont !! @@ -1902,8 +1903,8 @@ subroutine basal_stress_coeff_itd(G, IG, IST, sea_lev, CS) real :: rho_water ! water density [R ~> kg m-3] real :: pi ! [nondim] integer :: i, ii, j, isc, iec, jsc, jec, k, n, ncat - real :: ci_u ! Concentration at u-points [nondim] - real :: ci_v ! Concentration at u-points [nondim] + real :: ci_u ! Concentration at u-points [nondim] + real :: ci_v ! Concentration at u-points [nondim] isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ncat = IG%CatIce @@ -1956,50 +1957,54 @@ subroutine basal_stress_coeff_itd(G, IG, IST, sea_lev, CS) ! parameters for the log-normal mu_i = log(m_i/(CS%onemeter * sqrt(1.0 + v_i/m_i**2))) - sigma_i = max(sqrt(log(1.0 + v_i/m_i**2)), CS%puny) + sigma_i = sqrt(log(1.0 + v_i/m_i**2)) ! max thickness associated with percentile of log-normal PDF - ! x_kmax=x997 was obtained from an optimization procedure (Dupont et al.) - - x_kmax = CS%onemeter * exp(mu_i + sqrt(2.0*sigma_i)*CS%basal_stress_cutoff) - - ! Set x_kmax to hlev of the last category where there is ice - ! when there is no ice in the last category - cut = x_k(CS%ncat_i) - do n = ncat,-1,1 - if (acat(n) < CS%puny) then - cut = hin_max(n-1) - else - exit - endif - enddo - x_kmax = min(cut, x_kmax) - - g_k(:) = exp(-(log(x_k(:)/CS%onemeter) - mu_i) ** 2 / (2.0 * sigma_i ** 2)) / & - (x_k(:) * sigma_i * sqrt(2.0 * pi)) - - b_n(:) = exp(-(y_n(:) - mu_b) ** 2 / (2.0 * CS%sigma_b(i,j) ** 2)) / (CS%sigma_b(i,j) * sqrt(2.0*pi)) - - P_x(:) = g_k(:) * wid_i - P_y(:) = b_n(:) * wid_b - - do n =1, CS%ncat_i - if (x_k(n) > x_kmax) P_x(n)=0.0 - enddo - - ! calculate Tb factor at t-location - do n=1, CS%ncat_i - gt(:) = (y_n(:) <= rho_ice*x_k(n)/rho_water) - tmp(:) = merge(1,0,gt(:)) - ii = sum(tmp) - if (ii == 0) then - tb_tmp(n) = 0.0 - else - tb_tmp(n) = max(CS%basal_stress_mu_s * G%g_Earth * P_x(n) * & - sum(P_y(1:ii)*(rho_ice*x_k(n) - rho_water*y_n(1:ii))), 0.0) - endif - enddo - Tbt(i,j) = sum(tb_tmp) * exp(-CS%lemieux_alphab * (1.0 - atot)) + ! x_kmax=x997 was obtained from an optimization procedure (Dupont et al. 2022) + + if (sigma_i > 0) then + x_kmax = CS%onemeter * exp(mu_i + sqrt(2.0*sigma_i)*CS%basal_stress_cutoff) + + ! Set x_kmax to hlev of the last category where there is ice + ! when there is no ice in the last category + cut = x_k(CS%ncat_i) + do n = ncat,-1,1 + if (acat(n) < CS%puny) then + cut = hin_max(n-1) + else + exit + endif + enddo + x_kmax = min(cut, x_kmax) + + g_k(:) = exp(-(log(x_k(:)/CS%onemeter) - mu_i) ** 2 / (2.0 * sigma_i ** 2)) / & + (x_k(:) * sigma_i * sqrt(2.0 * pi)) + + b_n(:) = exp(-(y_n(:) - mu_b) ** 2 / (2.0 * CS%sigma_b(i,j) ** 2)) / (CS%sigma_b(i,j) * sqrt(2.0*pi)) + + P_x(:) = g_k(:) * wid_i + P_y(:) = b_n(:) * wid_b + + do n =1, CS%ncat_i + if (x_k(n) > x_kmax) P_x(n)=0.0 + enddo + + ! calculate Tb factor at t-location + do n=1, CS%ncat_i + gt(:) = (y_n(:) <= rho_ice*x_k(n)/rho_water) + tmp(:) = merge(1,0,gt(:)) + ii = sum(tmp) + if (ii == 0) then + tb_tmp(n) = 0.0 + else + tb_tmp(n) = max(CS%basal_stress_mu_s * G%g_Earth * P_x(n) * & + sum(P_y(1:ii)*(rho_ice*x_k(n) - rho_water*y_n(1:ii))), 0.0) + endif + enddo + Tbt(i,j) = sum(tb_tmp) * exp(-CS%lemieux_alphab * (1.0 - atot)) + else + Tbt(i,j) = 0.0 + endif endif enddo enddo diff --git a/src/SIS_open_boundary.F90 b/src/SIS_open_boundary.F90 index 13ccda0c..80445e01 100644 --- a/src/SIS_open_boundary.F90 +++ b/src/SIS_open_boundary.F90 @@ -270,7 +270,7 @@ subroutine open_boundary_config(G, US, param_file, OBC) (OBC%freeslip_vorticity .and. OBC%computed_vorticity) .or. & (OBC%freeslip_vorticity .and. OBC%specified_vorticity) .or. & (OBC%computed_vorticity .and. OBC%specified_vorticity)) & - call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config:\n"// & + call MOM_error(FATAL, "SIS_open_boundary.F90, open_boundary_config:\n"// & "Only one of OBC_ZERO_VORTICITY, OBC_FREESLIP_VORTICITY, OBC_COMPUTED_VORTICITY\n"// & "and OBC_IMPORTED_VORTICITY can be True at once.") call get_param(param_file, mdl, "OBC_ZERO_STRAIN", OBC%zero_strain, & @@ -364,7 +364,7 @@ subroutine open_boundary_config(G, US, param_file, OBC) elseif (segment_str(1:2) == 'J=') then call setup_v_point_obc(OBC, G, US, segment_str, l, param_file, reentrant_x) else - call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config: "//& + call MOM_error(FATAL, "SIS_open_boundary.F90, open_boundary_config: "//& "Unable to interpret "//segment_param_str//" = "//trim(segment_str)) endif enddo @@ -412,7 +412,7 @@ subroutine open_boundary_config(G, US, param_file, OBC) ! Safety check if ((OBC%open_u_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) .and. & .not.G%symmetric ) call MOM_error(FATAL, & - "MOM_open_boundary, open_boundary_config: "//& + "SIS_open_boundary, open_boundary_config: "//& "Symmetric memory must be used when using Flather OBCs.") if (.not.(OBC%specified_u_BCs_exist_globally .or. OBC%specified_v_BCs_exist_globally .or. & @@ -652,7 +652,7 @@ subroutine setup_u_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_y) OBC%segment(l_seg)%specified_grad = .true. OBC%segment(l_seg)%g_values_needed = .true. else - call MOM_error(FATAL, "MOM_open_boundary.F90, setup_u_point_obc: "//& + call MOM_error(FATAL, "SIS_open_boundary.F90, setup_u_point_obc: "//& "String '"//trim(action_str(a_loop))//"' not understood.") endif if (OBC%segment(l_seg)%nudged .or. OBC%segment(l_seg)%nudged_tan) then @@ -774,7 +774,7 @@ subroutine setup_v_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_x) OBC%segment(l_seg)%specified_grad = .true. OBC%segment(l_seg)%g_values_needed = .true. else - call MOM_error(FATAL, "MOM_open_boundary.F90, setup_v_point_obc: "//& + call MOM_error(FATAL, "SIS_open_boundary.F90, setup_v_point_obc: "//& "String '"//trim(action_str(a_loop))//"' not understood.") endif if (OBC%segment(l_seg)%nudged .or. OBC%segment(l_seg)%nudged_tan) then @@ -1039,14 +1039,14 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) do j=G%jsd,G%jed ; do i=G%isd,G%ied if (color(i,j) /= color2(i,j)) then fatal_error = .True. - write(mesg,'("MOM_open_boundary: problem with OBC segments specification at ",I5,",",I5," during\n", & + write(mesg,'("SIS_open_boundary: problem with OBC segments specification at ",I5,",",I5," during\n", & "the masking of the outside grid points.")') i, j call MOM_error(WARNING,"MOM register_tracer: "//mesg, all_print=.true.) endif if (color(i,j) == cout) G%bathyT(i,j) = min_depth enddo ; enddo if (fatal_error) call MOM_error(FATAL, & - "MOM_open_boundary: inconsistent OBC segments.") + "SIS_open_boundary: inconsistent OBC segments.") deallocate(color) deallocate(color2)