From 778d326c0a4404e3b6fc6403f2ce213636983ba3 Mon Sep 17 00:00:00 2001 From: Sam Hatfield Date: Tue, 25 Oct 2022 13:37:29 +0100 Subject: [PATCH 1/4] Fix floating-point exception in solar radiation subroutine There are cases where qsr can be extremely small (~10**-33) so when multiplied by r1_rau0_rcp (inverse of volumic mass of reference times ocean specific heat) it gives a subnormal number, which can trigger an overflow when qsr_hc divides by this value, or an underflow, which triggers a divide-by-zero. This commit augments the existing check for divide-by-zeros to account for this scenario. --- src/OCE/TRA/traqsr.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OCE/TRA/traqsr.F90 b/src/OCE/TRA/traqsr.F90 index 5f9c8e94..ae8b2a53 100644 --- a/src/OCE/TRA/traqsr.F90 +++ b/src/OCE/TRA/traqsr.F90 @@ -287,7 +287,7 @@ CONTAINS ! ! sea-ice: store the 1st ocean level attenuation coefficient DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) - IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) + IF( r1_rau0_rcp * qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) ELSE ; fraqsr_1lev(ji,jj) = 1._wp ENDIF END_2D -- GitLab From 520f45ee89f5c47485892d8379daf8a682377e60 Mon Sep 17 00:00:00 2001 From: Sam Hatfield Date: Tue, 25 Oct 2022 13:45:31 +0100 Subject: [PATCH 2/4] Prevent divide-by-zero errors in the TKE scheme --- src/OCE/ZDF/zdftke.F90 | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/OCE/ZDF/zdftke.F90 b/src/OCE/ZDF/zdftke.F90 index 8b6d9155..5680ec1e 100644 --- a/src/OCE/ZDF/zdftke.F90 +++ b/src/OCE/ZDF/zdftke.F90 @@ -219,6 +219,7 @@ CONTAINS REAL(wp), DIMENSION(A2D(nn_hls)) :: zice_fra, zhlc, zus3, zWlc2 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zpelc, zdiag, zd_up, zd_lw REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztmp ! for diags + REAL(wp) :: zdiv !!-------------------------------------------------------------------- ! zbbrau = rn_ebb / rho0 ! Local constant initialisation @@ -370,7 +371,16 @@ CONTAINS IF (rn2b(ji,jj,jk) <= 0.0_wp) then zri = 0.0_wp ELSE - zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) + ! This logic is to avoid divide-by-zero errors which can occur for single-precision + ! The actual value you choose for the denominator instead of zero doesn't really + ! matter, as long as it is very small and so triggers the same logic below with the + ! inverse Prandtl number + zdiv = p_sh2(ji,jj,jk) + rn_bshear + IF (zdiv == 0.0_wp) THEN + zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / rn_bshear + ELSE + zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / zdiv + ENDIF ENDIF ! ! inverse of Prandtl number apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) -- GitLab From f2332bff5e7f9d152c4962676f6f50989c252bef Mon Sep 17 00:00:00 2001 From: Sam Hatfield Date: Tue, 25 Oct 2022 13:55:12 +0100 Subject: [PATCH 3/4] Remove potential div/0 calculation in icethd.F90:ice_thd_1d2d Previously we just checked for zero ice thickness and snow depth. However we actually divide by the product of these two with ice concentration. When both multiplicands are very low, the product can be rounded to zero leading to a divide-by-zero floating-point exception. This bug occurred for a single-precision coupled run (eORCA025). It only happened after 15 years of a 5 member ensemble (a combined 75 years of simulation) so is quite rare. This bug could have potentially happened with double precision as well. I've also added an ELSEWHERE statement to handle the grid points where the denominator is zero. I just set the ice and snow enthalpy to zero. --- src/ICE/icethd.F90 | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/ICE/icethd.F90 b/src/ICE/icethd.F90 index 3971c441..d5348ae1 100644 --- a/src/ICE/icethd.F90 +++ b/src/ICE/icethd.F90 @@ -343,10 +343,18 @@ CONTAINS ! ! --- Change units of e_i, e_s from J/m2 to J/m3 --- ! DO jk = 1, nlay_i - WHERE( h_i_1d(1:npti)>0._wp ) e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i + WHERE( (h_i_1d(1:npti) * a_i_1d(1:npti)) > 0._wp ) + e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i + ELSEWHERE + e_i_1d(1:npti,jk) = 0._wp + ENDWHERE END DO DO jk = 1, nlay_s - WHERE( h_s_1d(1:npti)>0._wp ) e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s + WHERE( (h_s_1d(1:npti) * a_i_1d(1:npti)) > 0._wp ) + e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s + ELSEWHERE + e_s_1d(1:npti,jk) = 0._wp + ENDWHERE END DO ! ! !---------------------! -- GitLab From 213032fb6705056587257d6d02fa94ca7da0d850 Mon Sep 17 00:00:00 2001 From: Sam Hatfield Date: Tue, 25 Oct 2022 13:56:54 +0100 Subject: [PATCH 4/4] Secure J/m^2->J/m^3 conversion in SI3 against overflows This builds on 3485fd66f4e8512da32f7d0b13eeade350a3f69f. The ice/snow enthalpy is converted from J/m^2 to J/m^3 by dividing by the product of snow depth/ice thickness and sea-ice concentration. When one or both of the depth/thickness and concentration are close to zero this can cause arithmetic errors. The commit above secured this calculation against divide-by-zero errors by ensuring that it was never executed when the divisor was zero. However, it left open the possibility of overflows when the divisor was a subnormal number. In any case, such small numbers are physically meaningless so we don't have to think too hard about how to treat them. In this commit we check that the divisor is greater than 10^-20, and this prevents both overflows and divides-by-zero. --- src/ICE/icethd.F90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/ICE/icethd.F90 b/src/ICE/icethd.F90 index d5348ae1..fbd87f71 100644 --- a/src/ICE/icethd.F90 +++ b/src/ICE/icethd.F90 @@ -342,15 +342,18 @@ CONTAINS CALL tab_2d_1d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) ) ! ! --- Change units of e_i, e_s from J/m2 to J/m3 --- ! + ! Here we make sure that we don't divide by very small, but physically + ! meaningless, products of sea ice thicknesses/snow depths and sea ice + ! concentration DO jk = 1, nlay_i - WHERE( (h_i_1d(1:npti) * a_i_1d(1:npti)) > 0._wp ) + WHERE( (h_i_1d(1:npti) * a_i_1d(1:npti)) > epsi20 ) e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i ELSEWHERE e_i_1d(1:npti,jk) = 0._wp ENDWHERE END DO DO jk = 1, nlay_s - WHERE( (h_s_1d(1:npti) * a_i_1d(1:npti)) > 0._wp ) + WHERE( (h_s_1d(1:npti) * a_i_1d(1:npti)) > epsi20 ) e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s ELSEWHERE e_s_1d(1:npti,jk) = 0._wp -- GitLab