From 21a8603b58e597f76dd2a02aba7ccf35aee907ab Mon Sep 17 00:00:00 2001
From: clem <clement.rousset@locean.ipsl.fr>
Date: Fri, 7 Oct 2022 11:38:58 +0200
Subject: [PATCH] change rswitch into zswitch to gain 25% performance in
 rheology

---
 src/ICE/icedyn_adv_pra.F90 | 26 +++++++++++----------
 src/ICE/icedyn_rhg_eap.F90 | 46 ++++++++++++++++++++------------------
 src/ICE/icedyn_rhg_evp.F90 | 46 ++++++++++++++++++++------------------
 src/ICE/icesbc.F90         |  9 ++++----
 src/ICE/icethd.F90         | 10 +++++----
 src/ICE/icethd_dh.F90      | 10 ++++-----
 src/ICE/icethd_do.F90      | 30 +++++++++++++------------
 7 files changed, 94 insertions(+), 83 deletions(-)

diff --git a/src/ICE/icedyn_adv_pra.F90 b/src/ICE/icedyn_adv_pra.F90
index 763f71a66..b87b34a5b 100644
--- a/src/ICE/icedyn_adv_pra.F90
+++ b/src/ICE/icedyn_adv_pra.F90
@@ -377,6 +377,7 @@ CONTAINS
       REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
       REAL(wp) ::   zpsm, zps0
       REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy
+      REAL(wp) ::   zswitch
       REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
       REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      -
@@ -408,14 +409,14 @@ CONTAINS
             zs1max  = 1.5 * zslpmax
             zs1new  = MIN( zs1max, MAX( -zs1max, zpsx ) )
             zs2new  = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) )
-            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
+            zswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
             
             zps0  = zslpmax
-            zpsx  = zs1new  * rswitch
-            zpsxx = zs2new  * rswitch
-            zpsy  = zpsy    * rswitch
-            zpsyy = zpsyy   * rswitch
-            zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
+            zpsx  = zs1new  * zswitch
+            zpsxx = zs2new  * zswitch
+            zpsy  = zpsy    * zswitch
+            zpsyy = zpsyy   * zswitch
+            zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * zswitch
             
             !  Calculate fluxes and moments between boxes i<-->i+1
             !                                !  Flux from i to i+1 WHEN u GT 0
@@ -566,6 +567,7 @@ CONTAINS
       REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
       REAL(wp) ::   zpsm, zps0
       REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy
+      REAL(wp) ::   zswitch
       REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
       REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
@@ -597,14 +599,14 @@ CONTAINS
             zs1max  = 1.5 * zslpmax
             zs1new  = MIN( zs1max, MAX( -zs1max, zpsy ) )
             zs2new  = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new ) - zslpmax, zpsyy ) )
-            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
+            zswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
             !
             zps0  = zslpmax
-            zpsx  = zpsx    * rswitch
-            zpsxx = zpsxx   * rswitch
-            zpsy  = zs1new  * rswitch
-            zpsyy = zs2new  * rswitch
-            zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
+            zpsx  = zpsx    * zswitch
+            zpsxx = zpsxx   * zswitch
+            zpsy  = zs1new  * zswitch
+            zpsyy = zs2new  * zswitch
+            zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * zswitch
 
             !  Calculate fluxes and moments between boxes j<-->j+1
             !                                !  Flux from j to j+1 WHEN v GT 0
diff --git a/src/ICE/icedyn_rhg_eap.F90 b/src/ICE/icedyn_rhg_eap.F90
index e87e879a3..cf2fd163f 100644
--- a/src/ICE/icedyn_rhg_eap.F90
+++ b/src/ICE/icedyn_rhg_eap.F90
@@ -153,6 +153,8 @@ CONTAINS
       REAL(wp) ::   zalphar, zalphas                                    ! for mechanical redistribution
       REAL(wp) ::   zmresult11, zmresult12, z1dtevpkth, zp5kth, z1_dtevp_A  ! for structure tensor evolution
       !
+      REAL(wp) ::   zswitch
+      !
       REAL(wp), DIMENSION(jpi,jpj) ::   zstress12tmp                    ! anisotropic stress tensor component for regridding
       REAL(wp), DIMENSION(A2D(0))  ::   zyield11, zyield22, zyield12    ! yield surface tensor for history
       REAL(wp), DIMENSION(jpi,jpj) ::   zdelta, zp_delt                 ! delta and P/delta at T points
@@ -551,23 +553,23 @@ CONTAINS
                !
                !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
                !                                         1 = sliding friction : TauB < RHS
-               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+               zswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
                !
                IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
                   zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
-                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
+                  v_ice(ji,jj) = ( (          zswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
                      &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &
+                     &            + ( 1._wp - zswitch ) * (  v_ice_b(ji,jj)                                                   &
                      &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
                      &                                    ) / ( zbetav + 1._wp )                                              &
                      &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &           )   * zmsk00y(ji,jj)
                ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
-                  v_ice(ji,jj) = ( (         rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
+                  v_ice(ji,jj) = ( (         zswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
                      &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
+                     &            + ( 1._wp - zswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
                      &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                  & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &            )   * zmsk00y(ji,jj)
                ENDIF
@@ -597,23 +599,23 @@ CONTAINS
                !
                !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
                !                                         1 = sliding friction : TauB < RHS
-               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+               zswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
                !
                IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
                   zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
-                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
+                  u_ice(ji,jj) = ( (          zswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
                      &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
+                     &            + ( 1._wp - zswitch ) * (  u_ice_b(ji,jj)                                                   &
                      &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
                      &                                    ) / ( zbetau + 1._wp )                                              &
                      &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &           )   * zmsk00x(ji,jj)
                ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
-                  u_ice(ji,jj) = ( (         rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
+                  u_ice(ji,jj) = ( (         zswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
                      &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
+                     &            + ( 1._wp - zswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
                      &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                  & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &            )   * zmsk00x(ji,jj)
                ENDIF
@@ -647,23 +649,23 @@ CONTAINS
                !
                !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
                !                                         1 = sliding friction : TauB < RHS
-               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+               zswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
                !
                IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
                   zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
-                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
+                  u_ice(ji,jj) = ( (          zswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
                      &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
+                     &            + ( 1._wp - zswitch ) * (  u_ice_b(ji,jj)                                                   &
                      &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
                      &                                    ) / ( zbetau + 1._wp )                                              &
                      &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &           )   * zmsk00x(ji,jj)
                ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
-                  u_ice(ji,jj) = ( (         rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
+                  u_ice(ji,jj) = ( (         zswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
                      &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
+                     &            + ( 1._wp - zswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
                      &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                  & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &            )   * zmsk00x(ji,jj)
                ENDIF
@@ -693,23 +695,23 @@ CONTAINS
                !
                !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
                !                                         1 = sliding friction : TauB < RHS
-               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+               zswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
                !
                IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
                   zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
-                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
+                  v_ice(ji,jj) = ( (          zswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
                      &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &
+                     &            + ( 1._wp - zswitch ) * (  v_ice_b(ji,jj)                                                   &
                      &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
                      &                                    ) / ( zbetav + 1._wp )                                              &
                      &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &           )   * zmsk00y(ji,jj)
                ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
-                  v_ice(ji,jj) = ( (         rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
+                  v_ice(ji,jj) = ( (         zswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
                      &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
+                     &            + ( 1._wp - zswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
                      &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                  & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &            )   * zmsk00y(ji,jj)
                ENDIF
@@ -780,8 +782,8 @@ CONTAINS
          zdelta(ji,jj)   = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta
 
          ! delta at T points
-         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0
-         pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * rswitch  
+         zswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0
+         pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * zswitch  
                            ! it seems that deformation used for advection and mech redistribution is delta*
                            ! MV in principle adding creep limit is a regularization for viscosity not for delta
                            ! delta_star should not (in my view) be used in a replacement for delta
diff --git a/src/ICE/icedyn_rhg_evp.F90 b/src/ICE/icedyn_rhg_evp.F90
index 136f13d12..a2b075d47 100644
--- a/src/ICE/icedyn_rhg_evp.F90
+++ b/src/ICE/icedyn_rhg_evp.F90
@@ -135,6 +135,8 @@ CONTAINS
       !
       REAL(wp) ::   zfac_x, zfac_y
       !
+      REAL(wp) ::   zswitch
+      !
       REAL(wp), DIMENSION(jpi,jpj)       ::   zdelta, zp_delt                 ! delta, P/delta at T points
       REAL(wp), DIMENSION(jpi,jpj)       ::   zbeta                           ! beta coef from Kimmritz 2017
       !
@@ -519,23 +521,23 @@ CONTAINS
                !
                !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
                !                                         1 = sliding friction : TauB < RHS
-               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+               zswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
                !
                IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
                   zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
-                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
+                  v_ice(ji,jj) = ( (          zswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
                      &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &
+                     &            + ( 1._wp - zswitch ) * (  v_ice_b(ji,jj)                                                   &
                      &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
                      &                                    ) / ( zbetav + 1._wp )                                              &
                      &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &           )   * zmsk00y(ji,jj)
                ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
-                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
+                  v_ice(ji,jj) = ( (          zswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
                      &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
+                     &            + ( 1._wp - zswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
                      &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &            )  * zmsk00y(ji,jj)
                ENDIF
@@ -565,23 +567,23 @@ CONTAINS
                !
                !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
                !                                         1 = sliding friction : TauB < RHS
-               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+               zswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
                !
                IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
                   zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
-                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
+                  u_ice(ji,jj) = ( (          zswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
                      &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
+                     &            + ( 1._wp - zswitch ) * (  u_ice_b(ji,jj)                                                   &
                      &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
                      &                                    ) / ( zbetau + 1._wp )                                              &
                      &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &           )   * zmsk00x(ji,jj)
                ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
-                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
+                  u_ice(ji,jj) = ( (          zswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
                      &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
+                     &            + ( 1._wp - zswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
                      &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &           )   * zmsk00x(ji,jj)
                ENDIF
@@ -615,23 +617,23 @@ CONTAINS
                !
                !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
                !                                         1 = sliding friction : TauB < RHS
-               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+               zswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
                !
                IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
                   zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
-                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
+                  u_ice(ji,jj) = ( (          zswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
                      &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
+                     &            + ( 1._wp - zswitch ) * (  u_ice_b(ji,jj)                                                   &
                      &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
                      &                                    ) / ( zbetau + 1._wp )                                              &
                      &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &           )   * zmsk00x(ji,jj)
                ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
-                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
+                  u_ice(ji,jj) = ( (          zswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
                      &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
+                     &            + ( 1._wp - zswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
                      &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &           )   * zmsk00x(ji,jj)
                ENDIF
@@ -661,23 +663,23 @@ CONTAINS
                !
                !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
                !                                         1 = sliding friction : TauB < RHS
-               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+               zswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
                !
                IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
                   zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
-                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
+                  v_ice(ji,jj) = ( (          zswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
                      &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &
+                     &            + ( 1._wp - zswitch ) * (  v_ice_b(ji,jj)                                                   &
                      &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
                      &                                    ) / ( zbetav + 1._wp )                                              &
                      &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &           )   * zmsk00y(ji,jj)
                ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
-                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
+                  v_ice(ji,jj) = ( (          zswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
                      &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
-                     &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
+                     &            + ( 1._wp - zswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
                      &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                      &           )   * zmsk00y(ji,jj)
                ENDIF
@@ -773,8 +775,8 @@ CONTAINS
          zdelta(ji,jj)   = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) * zmsk(ji,jj) ! delta
 
          ! delta* at T points (pdelta_i)
-         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0
-         pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * rswitch  
+         zswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0
+         pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * zswitch  
                            ! it seems that deformation used for advection and mech redistribution is delta*
                            ! MV in principle adding creep limit is a regularization for viscosity not for delta
                            ! delta_star should not (in my view) be used in a replacement for delta
diff --git a/src/ICE/icesbc.F90 b/src/ICE/icesbc.F90
index 1efec2099..7e357b47d 100644
--- a/src/ICE/icesbc.F90
+++ b/src/ICE/icesbc.F90
@@ -295,6 +295,7 @@ CONTAINS
       !!-----------------------------------------------------------------------
       INTEGER  ::   ji, jj             ! dummy loop indices
       REAL(wp) ::   zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos, zu_io, zv_io, zu_iom1, zv_iom1
+      REAL(wp) ::   zswitch
       REAL(wp), PARAMETER ::   zfric_umin = 0._wp       ! lower bound for the friction velocity (cice value=5.e-04)
       REAL(wp), PARAMETER ::   zch        = 0.0057_wp   ! heat transfer coefficient
       REAL(wp), DIMENSION(A2D(0)) ::  zfric, zvel       ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
@@ -324,7 +325,7 @@ CONTAINS
       ! Partial computation of forcing for the thermodynamic sea ice model
       !--------------------------------------------------------------------!
       DO_2D( 0, 0, 0, 0 )   ! needed for qlead
-         rswitch  = smask0(ji,jj) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
+         zswitch  = smask0(ji,jj) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
          !
          ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- !
          zqld =  smask0(ji,jj) * rDt_ice *  &
@@ -340,13 +341,13 @@ CONTAINS
          ! --- Sensible ocean-to-ice heat flux (W/m2) --- !
          !     (mostly>0 but <0 if supercooling)
          zfric_u            = MAX( SQRT( zfric(ji,jj) ), zfric_umin )
-         qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) )
+         qsb_ice_bot(ji,jj) = zswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) )
 
          ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach
          !                              the freezing point, so that we do not have SST < T_freeze
          !                              This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg
          !                              The following formulation is ok for both normal conditions and supercooling
-         qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) )
+         qsb_ice_bot(ji,jj) = zswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) )
 
          ! If conditions are always supercooled (such as at the mouth of ice-shelves), then ice grows continuously
          ! ==> stop ice formation by artificially setting up the turbulent fluxes to 0 when volume > 20m (arbitrary)
@@ -365,7 +366,7 @@ CONTAINS
             !                        but we have to make sure that this heat will not make the sst drop below the freezing point
             !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos
             !                        The following formulation is ok for both normal conditions and supercooling
-            fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) &  ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
+            fhld (ji,jj) = zswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) &  ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
                &                                 - qsb_ice_bot(ji,jj) )
             qlead(ji,jj) = 0._wp
          ELSE
diff --git a/src/ICE/icethd.F90 b/src/ICE/icethd.F90
index 448ae5aa9..a0b8f0c70 100644
--- a/src/ICE/icethd.F90
+++ b/src/ICE/icethd.F90
@@ -202,6 +202,7 @@ CONTAINS
       !!-------------------------------------------------------------------
       INTEGER  ::   ji, jk   ! dummy loop indices
       REAL(wp) ::   ztmelts, zbbb, zccc  ! local scalar
+      REAL(wp) ::   zswitch
       !!-------------------------------------------------------------------
       ! Recover ice temperature
       DO jk = 1, nlay_i
@@ -213,8 +214,8 @@ CONTAINS
             t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_rcpi
 
             ! mask temperature
-            rswitch       = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) )
-            t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0
+            zswitch       = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) )
+            t_i_1d(ji,jk) = zswitch * t_i_1d(ji,jk) + ( 1._wp - zswitch ) * rt0
          END DO
       END DO
       !
@@ -232,6 +233,7 @@ CONTAINS
       REAL(wp) ::   zhi_bef            ! ice thickness before thermo
       REAL(wp) ::   zdh_mel, zda_mel   ! net melting
       REAL(wp) ::   zvi, zvs           ! ice/snow volumes
+      REAL(wp) ::   zswitch
       !!-----------------------------------------------------------------------
       !
       DO ji = 1, npti
@@ -241,8 +243,8 @@ CONTAINS
             zvs          = a_i_1d(ji) * h_s_1d(ji)
             ! lateral melting = concentration change
             zhi_bef     = h_i_1d(ji) - zdh_mel
-            rswitch     = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) )
-            zda_mel     = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) )
+            zswitch     = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) )
+            zda_mel     = zswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) )
             a_i_1d(ji)  = MAX( epsi20, a_i_1d(ji) + zda_mel )
             ! adjust thickness
             h_i_1d(ji) = zvi / a_i_1d(ji)
diff --git a/src/ICE/icethd_dh.F90 b/src/ICE/icethd_dh.F90
index a0d40b7c8..0d7497f62 100644
--- a/src/ICE/icethd_dh.F90
+++ b/src/ICE/icethd_dh.F90
@@ -96,7 +96,7 @@ CONTAINS
       REAL(wp), DIMENSION(0:nlay_i+1) ::   zh_i_old  ! old thickness
       REAL(wp), DIMENSION(0:nlay_i+1) ::   ze_i_old  ! old enthalpy
 
-      REAL(wp) ::   zswitch_sal
+      REAL(wp) ::   zswitch, zswitch_sal
 
       INTEGER  ::   num_iter_max      ! Heat conservation
       !!------------------------------------------------------------------
@@ -202,8 +202,8 @@ CONTAINS
          DO jk = 0, nlay_s
             IF( zh_s(jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN
                !
-               rswitch = MAX( 0._wp , SIGN( 1._wp , ze_s(jk) - epsi20 ) )
-               zdum    = - rswitch * zq_top(ji) / MAX( ze_s(jk), epsi20 )   ! thickness change
+               zswitch = MAX( 0._wp , SIGN( 1._wp , ze_s(jk) - epsi20 ) )
+               zdum    = - zswitch * zq_top(ji) / MAX( ze_s(jk), epsi20 )   ! thickness change
                zdum    = MAX( zdum , - zh_s(jk) )                           ! bound melting
 
                hfx_snw_1d    (ji) = hfx_snw_1d    (ji) - ze_s(jk) * zdum * a_i_1d(ji) * r1_Dt_ice   ! heat used to melt snow(W.m-2, >0)
@@ -324,8 +324,8 @@ CONTAINS
             ! record which layers have disappeared (for bottom melting)
             !    => icount=0 : no layer has vanished
             !    => icount=5 : 5 layers have vanished
-            rswitch    = MAX( 0._wp , SIGN( 1._wp , - zh_i(jk) ) )
-            icount(jk) = NINT( rswitch )
+            zswitch    = MAX( 0._wp , SIGN( 1._wp , - zh_i(jk) ) )
+            icount(jk) = NINT( zswitch )
 
          END DO
 
diff --git a/src/ICE/icethd_do.F90 b/src/ICE/icethd_do.F90
index 063c3e903..6545a9708 100644
--- a/src/ICE/icethd_do.F90
+++ b/src/ICE/icethd_do.F90
@@ -88,6 +88,7 @@ CONTAINS
       REAL(wp) ::   zdv_res     ! residual volume in case of excessive heat budget
       REAL(wp) ::   zda_res     ! residual area in case of excessive heat budget
       REAL(wp) ::   zv_frazb    ! accretion of frazil ice at the ice bottom
+      REAL(wp) ::   zswitch
       !
       REAL(wp), DIMENSION(jpl) ::   zv_b    ! old volume of ice in category jl
       REAL(wp), DIMENSION(jpl) ::   za_b    ! old area of ice in category jl
@@ -205,9 +206,9 @@ CONTAINS
             sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice * rhoi * zs_newice(ji) * r1_Dt_ice
          
             ! A fraction fraz_frac of frazil ice is accreted at the ice bottom
-            rswitch   = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
-            zv_frazb  =           zfraz_frac_1d(ji) * rswitch   * zv_newice
-            zv_newice = ( 1._wp - zfraz_frac_1d(ji) * rswitch ) * zv_newice
+            zswitch   = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
+            zv_frazb  =           zfraz_frac_1d(ji) * zswitch   * zv_newice
+            zv_newice = ( 1._wp - zfraz_frac_1d(ji) * zswitch ) * zv_newice
          
             ! --- Area of new ice --- !
             za_newice = zv_newice / zh_newice(ji)
@@ -245,10 +246,10 @@ CONTAINS
 
             DO jk = 1, nlay_i
                jl = jcat
-               rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) )
+               zswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) )
                e_i_2d(ji,jk,jl) = zswinew   *   ze_newice +                                            &
                   &       ( 1.0 - zswinew ) * ( ze_newice * zv_newice + e_i_2d(ji,jk,jl) * zv_b(jl) )  &
-                  &        * rswitch / MAX( v_i_2d(ji,jl), epsi20 )
+                  &        * zswitch / MAX( v_i_2d(ji,jl), epsi20 )
             END DO
          
             ! --- bottom ice growth + ice enthalpy remapping --- !
@@ -263,9 +264,9 @@ CONTAINS
                END DO
 
                ! new volumes including lateral/bottom accretion + residual
-               rswitch       = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) )
-               zv_newfra     = rswitch * ( zdv_res + zv_frazb ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 )
-               a_i_2d(ji,jl) = rswitch * a_i_2d(ji,jl)               
+               zswitch       = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) )
+               zv_newfra     = zswitch * ( zdv_res + zv_frazb ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 )
+               a_i_2d(ji,jl) = zswitch * a_i_2d(ji,jl)               
                v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newfra
                ! for remapping
                zh_i_old(nlay_i+1) = zv_newfra
@@ -322,6 +323,7 @@ CONTAINS
       INTEGER  ::   ji, jj             ! dummy loop indices
       INTEGER  ::   iter
       REAL(wp) ::   zvfrx, zvgx, ztaux, zf, ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp
+      REAL(wp) ::   zswitch
       REAL(wp), PARAMETER ::   zcai    = 1.4e-3_wp                       ! ice-air drag (clem: should be dependent on coupling/forcing used)
       REAL(wp), PARAMETER ::   zhicrit = 0.04_wp                         ! frazil ice thickness
       REAL(wp), PARAMETER ::   zsqcd   = 1.0_wp / SQRT( 1.3_wp * zcai )  ! 1/SQRT(airdensity*drag)
@@ -357,20 +359,20 @@ CONTAINS
                ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
 
                ! -- Frazil ice velocity -- !
-               rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
-               zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
-               zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
+               zswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
+               zvfrx   = zswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
+               zvfry   = zswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
 
                ! -- Pack ice velocity -- !
                zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
                zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
 
                ! -- Relative frazil/pack ice velocity -- !
-               rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
-               zvrel2  = MAX( (zvfrx - zvgx)*(zvfrx - zvgx) + (zvfry - zvgy)*(zvfry - zvgy), 0.15_wp*0.15_wp ) * rswitch
+               zswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
+               zvrel2  = MAX( (zvfrx - zvgx)*(zvfrx - zvgx) + (zvfry - zvgy)*(zvfry - zvgy), 0.15_wp*0.15_wp ) * zswitch
 
                ! -- fraction of frazil ice -- !
-               fraz_frac(ji,jj) = rswitch * ( TANH( rn_Cfraz * ( SQRT(zvrel2) - rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz
+               fraz_frac(ji,jj) = zswitch * ( TANH( rn_Cfraz * ( SQRT(zvrel2) - rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz
                
                ! -- new ice thickness (iterative loop) -- !
                ht_i_new(ji,jj) = zhicrit +   ( zhicrit + 0.1_wp )    &
-- 
GitLab