diff --git a/src/ICE/icethd.F90 b/src/ICE/icethd.F90
index 3971c4415240950ed6423023821b958247370d3f..fbd87f717158afabb8540444462742d147a4925e 100644
--- a/src/ICE/icethd.F90
+++ b/src/ICE/icethd.F90
@@ -342,11 +342,22 @@ 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)>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)) > 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)>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)) > 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
+            ENDWHERE
          END DO
          !
          !                 !---------------------!
diff --git a/src/OCE/TRA/traqsr.F90 b/src/OCE/TRA/traqsr.F90
index 5f9c8e94e44f98c29f0cbe76fa1d50aa42332f65..ae8b2a53f8cdb91a4db2c2fe5d9c10528fcb2f5d 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
diff --git a/src/OCE/ZDF/zdftke.F90 b/src/OCE/ZDF/zdftke.F90
index 8b6d915587952429fdf3237b4bb8d2718b7fa85a..5680ec1e20176e01de75677426a40208fdf571f2 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 )  )