Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8486)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8498)
@@ -48,14 +48,9 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
-
+   !
    REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier
    REAL(wp), PARAMETER ::   kraft   = 0.5_wp    ! rafting multipliyer
-
-!!gm Cp is  1) not DOCTOR, 
-!!          2) misleading name: heat capacity instead of a constant,
-!!          3) recomputed at each time-step, whereas it is stored in the module memory...
-!!             ===>>> compute it one for all inside the IF( kt == nit000 ) (i.e. without the ".AND. lwp")
-   REAL(wp)            ::   Cp                  ! ??? !!gm  Not doctor ! 
-   
+   !
+   REAL(wp)            ::   zdrho               ! 
    !
    !
@@ -111,6 +106,5 @@
       INTEGER  ::   niter              ! local integer 
       INTEGER  ::   iterate_ridging    ! if =1, repeat the ridging
-      REAL(wp) ::   za, zfac, zcs_2    ! local scalar
-      CHARACTER (len = 15) ::   fieldid
+      REAL(wp) ::   za                 ! local scalar
       REAL(wp), DIMENSION(jpi,jpj) ::   closing_net     ! net rate at which area is removed    (1/s)
       !                                                 ! (ridging ice area - area of new ridges) / dt
@@ -125,22 +119,12 @@
       IF( nn_timing == 1 )  CALL timing_start('icerdgrft')
 
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*)'icerdgrft : ice ridging and rafting'
-         WRITE(numout,*)'~~~~~~~~~~'
-      ENDIF
-!!gm should be:      
-!      IF( kt == nit000 ) THEN
-!         IF(lwp) WRITE(numout,*)
-!         IF(lwp) WRITE(numout,*)'icerdgrft : ???'
-!         IF(lwp) WRITE(numout,*)'~~~~~~~~~~'
-!         !
-!         Cp    = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE
-!         !
-!!gm why not adding here zcs_2 computation
-!         !
-!      ENDIF
-!!gm end
-      
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*)'icerdgrft : ice ridging and rafting'
+         IF(lwp) WRITE(numout,*)'~~~~~~~~~~'
+         !
+         zdrho = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE
+         !
+      ENDIF      
       !                    ! conservation test
       IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
@@ -149,6 +133,4 @@
       ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
       !-----------------------------------------------------------------------------!
-      Cp    = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE
-      zcs_2 = rn_cs * 0.5_wp
       !
       CALL ice_rdgrft_ridgeprep                             ! prepare ridging
@@ -175,5 +157,5 @@
             !  (thick, newly ridged ice).
 
-            closing_net(ji,jj) = zcs_2 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
+            closing_net(ji,jj) = rn_cs * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
 
             ! 2.2 divu_adv
@@ -220,8 +202,6 @@
                za   = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice
                IF    ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN                  ! would lead to negative ato_i
-                  zfac          = - ato_i(ji,jj) / za
                   opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 
                ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN  ! would lead to ato_i > asum
-                  zfac          = ( asum(ji,jj) - ato_i(ji,jj) ) / za
                   opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 
                ENDIF
@@ -238,6 +218,5 @@
                   za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
                   IF( za  >  a_i(ji,jj,jl) ) THEN
-                     zfac = a_i(ji,jj,jl) / za
-                     closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
+                     closing_gross(ji,jj) = closing_gross(ji,jj) * a_i(ji,jj,jl) / za
                   ENDIF
                END DO
@@ -254,5 +233,5 @@
          asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
 
-         ! 3.5 Do we keep on iterating ???
+         ! 3.5 Do we keep on iterating?
          !-----------------------------------------------------------------------------!
          ! Check whether asum = 1.  If not (because the closing and opening
@@ -316,34 +295,17 @@
       !!---------------------------------------------------------------------!
       INTEGER  ::   ji, jj, jl    ! dummy loop indices
-      REAL(wp) ::   Gstari, astari, hrmean, zdummy   ! local scalar     !!gm DOCTOR norme should start with z !!!!!
-      REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n
+      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zdummy   ! local scalar
+      REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::   zGsum     ! zGsum(n) = sum of areas in categories 0 to n
       !------------------------------------------------------------------------------!
 
-      Gstari        = 1._wp / rn_gstar    
-      astari        = 1._wp / rn_astar    
-      aksum(:,:)    = 0._wp
-      athorn(:,:,:) = 0._wp
-      aridge(:,:,:) = 0._wp
-      araft (:,:,:) = 0._wp
+      z1_gstar      = 1._wp / rn_gstar    
+      z1_astar      = 1._wp / rn_astar    
 
       CALL ice_var_zapsmall   ! Zero out categories with very small areas
 
-!      DO jl = 1, jpl          ! Ice thickness needed for rafting
-!         DO jj = 1, jpj
-!            DO ji = 1, jpi
-!!gm               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
-!!gm               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
-!               IF( a_i(ji,jj,jl) >= epsi20 ) THEN   ;   ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl)
-!               ELSE                                 ;   ht_i(ji,jj,jl) = 0._wp
-!               ENDIF
-!           END DO
-!         END DO
-!      END DO
-!!gm or even better :
-!     !                       ! Ice thickness needed for rafting
+      !                       ! Ice thickness needed for rafting
       WHERE( a_i(:,:,:) >= epsi20 )   ;   ht_i(:,:,:) = v_i (:,:,:) / a_i(:,:,:)
       ELSEWHERE                       ;   ht_i(:,:,:) = 0._wp
       END WHERE
-!!gm end
 
       !------------------------------------------------------------------------------!
@@ -356,17 +318,11 @@
       !
       ! Compute cumulative thickness distribution function
-      ! Compute the cumulative thickness distribution function Gsum,
-      ! where Gsum(n) is the fractional area in categories 0 to n.
+      ! Compute the cumulative thickness distribution function zGsum,
+      ! where zGsum(n) is the fractional area in categories 0 to n.
       ! initial value (in h = 0) equals open water area
-      Gsum(:,:,-1) = 0._wp
-      Gsum(:,:,0 ) = ato_i(:,:)
-      ! for each value of h, you have to add ice concentration then
+      zGsum(:,:,-1) = 0._wp
+      zGsum(:,:,0 ) = ato_i(:,:) / asum(:,:)
       DO jl = 1, jpl
-         Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl)
-      END DO
-      !
-      ! Normalize the cumulative distribution to 1
-      DO jl = 0, jpl
-         Gsum(:,:,jl) = Gsum(:,:,jl) / asum(:,:)
+         zGsum(:,:,jl) = ( ato_i(:,:) + SUM( a_i(:,:,1:jl), dim=3 ) ) / asum(:,:) ! sum(1:jl) is ok (and not jpl)
       END DO
 
@@ -388,10 +344,10 @@
             DO jj = 1, jpj 
                DO ji = 1, jpi
-                  IF    ( Gsum(ji,jj,jl)   < rn_gstar ) THEN
-                     athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * &
-                        &                        ( 2._wp - ( Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari )
-                  ELSEIF( Gsum(ji,jj,jl-1) < rn_gstar ) THEN
-                     athorn(ji,jj,jl) = Gstari * ( rn_gstar       - Gsum(ji,jj,jl-1) ) *  &
-                        &                        ( 2._wp - ( Gsum(ji,jj,jl-1) + rn_gstar       ) * Gstari )
+                  IF    ( zGsum(ji,jj,jl)   < rn_gstar ) THEN
+                     athorn(ji,jj,jl) = z1_gstar * ( zGsum(ji,jj,jl) - zGsum(ji,jj,jl-1) ) * &
+                        &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + zGsum(ji,jj,jl) ) * z1_gstar )
+                  ELSEIF( zGsum(ji,jj,jl-1) < rn_gstar ) THEN
+                     athorn(ji,jj,jl) = z1_gstar * ( rn_gstar        - zGsum(ji,jj,jl-1) ) *  &
+                        &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + rn_gstar        ) * z1_gstar )
                   ELSE
                      athorn(ji,jj,jl) = 0._wp
@@ -403,10 +359,10 @@
       ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007)
          !                        
-         zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array
+         zdummy = 1._wp / ( 1._wp - EXP(-z1_astar) )        ! precompute exponential terms using zGsum as a work array
          DO jl = -1, jpl
-            Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
+            zGsum(:,:,jl) = EXP( -zGsum(:,:,jl) * z1_astar ) * zdummy
          END DO
          DO jl = 0, jpl
-             athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
+            athorn(:,:,jl) = zGsum(:,:,jl-1) - zGsum(:,:,jl)
          END DO
          !
@@ -418,18 +374,18 @@
             DO jj = 1, jpj 
                DO ji = 1, jpi
-                  zdummy           = TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) )
-                  aridge(ji,jj,jl) = ( 1._wp + zdummy ) * 0.5_wp * athorn(ji,jj,jl)
+                  aridge(ji,jj,jl) = ( 1._wp + TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) ) * 0.5_wp * athorn(ji,jj,jl)
                   araft (ji,jj,jl) = athorn(ji,jj,jl) - aridge(ji,jj,jl)
                END DO
             END DO
          END DO
-      ELSEIF( ln_ridging .AND. .NOT.ln_rafting ) THEN   !- ridging alone
-         DO jl = 1, jpl
-            aridge(:,:,jl) = athorn(:,:,jl)
-         END DO
-      ELSEIF( ln_rafting .AND. .NOT.ln_ridging ) THEN   !- rafting alone   
-         DO jl = 1, jpl
-            araft(:,:,jl) = athorn(:,:,jl)
-         END DO
+      ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone
+         aridge(:,:,:) = athorn(:,:,1:jpl)
+         araft (:,:,:) = 0._wp
+      ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone   
+        aridge(:,:,:) = 0._wp
+         araft (:,:,:) = athorn(:,:,1:jpl)
+      ELSE                                               !- no ridging & no rafting
+         aridge(:,:,:) = 0._wp
+         araft (:,:,:) = 0._wp         
       ENDIF
 
@@ -441,13 +397,13 @@
       ! 
       ! This parameterization is a modified version of Hibler (1980).
-      ! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
+      ! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
       !  and for very thick ridging ice must be >= krdgmin*hi
       !
       ! The minimum ridging thickness, hrmin, is equal to 2*hi 
       !  (i.e., rafting) and for very thick ridging ice is
-      !  constrained by hrmin <= (hrmean + hi)/2.
+      !  constrained by hrmin <= (zhmean + hi)/2.
       ! 
       ! The maximum ridging thickness, hrmax, is determined by
-      !  hrmean and hrmin.
+      !  zhmean and hrmin.
       !
       ! These modifications have the effect of reducing the ice strength
@@ -466,9 +422,9 @@
             DO ji = 1, jpi
                IF ( athorn(ji,jj,jl) > 0._wp ) THEN
-                  hrmean          = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin )
-                  hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( hrmean + ht_i(ji,jj,jl) ) )
-                  hrmax(ji,jj,jl) = 2._wp * hrmean - hrmin(ji,jj,jl)
+                  zhmean          = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin )
+                  hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( zhmean + ht_i(ji,jj,jl) ) )
+                  hrmax(ji,jj,jl) = 2._wp * zhmean - hrmin(ji,jj,jl)
                   hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft
-                  krdg (ji,jj,jl) = ht_i(ji,jj,jl) / MAX( hrmean, epsi20 )
+                  krdg (ji,jj,jl) = ht_i(ji,jj,jl) / MAX( zhmean, epsi20 )
                   !
                   ! Normalization factor : aksum, ensures mass conservation
@@ -498,16 +454,13 @@
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges
-      !
-      CHARACTER (len=80) ::   fieldid   ! field identifier     !!gm DOCTOR name wrong
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area retreats, excluding area of new ridges
       !
       INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
       INTEGER ::   ij                ! horizontal index, combines i and j loops
       INTEGER ::   icells            ! number of cells with a_i > puny
-      REAL(wp) ::   hL, hR, farea    ! left and right limits of integration
-      REAL(wp) ::   zwfx_snw         ! snow mass flux increment
+      REAL(wp) ::   hL, hR, farea    ! left and right limits of integration and new area going to jl2
 
       INTEGER , DIMENSION(jpij) ::   indxi, indxj   ! compressed indices
-      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol   ! new ridge volume going to n2
+      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol   ! new ridge volume going to jl2
 
       REAL(wp), DIMENSION(jpij) ::   afrac            ! fraction of category area ridged 
@@ -656,8 +609,8 @@
             !  During the next time step, thermo_rates will determine whether
             !  the ocean cools or new ice grows.
-            zwfx_snw           = ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
-               &                 + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
-
-            wfx_snw_dyn(ji,jj)  =   wfx_snw_dyn(ji,jj) + zwfx_snw
+            wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
+               &                                      + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
+
+!!!            wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + zwfx_snw
 
             hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
@@ -722,5 +675,5 @@
 
                ! update jl1
-               e_i  (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
+               e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
 
             END DO
@@ -806,6 +759,6 @@
       !
       INTEGER             ::   ji,jj, jl   ! dummy loop indices
-      INTEGER             ::   ksmooth     ! smoothing the resistance to deformation    !!gm not DOCTOR : start with i !!!
-      INTEGER             ::   numts_rm    ! number of time steps for the P smoothing    !!gm not DOCTOR : start with i !!!
+      INTEGER             ::   ismooth     ! smoothing the resistance to deformation
+      INTEGER             ::   itframe     ! number of time steps for the P smoothing
       REAL(wp)            ::   zp, z1_3    ! local scalars
       REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
@@ -813,68 +766,39 @@
       !!----------------------------------------------------------------------
 
-      !------------------------------------------------------------------------------!
-      ! 1) Initialize
-      !------------------------------------------------------------------------------!
-      strength(:,:) = 0._wp
-
-      !------------------------------------------------------------------------------!
-      ! 2) Compute thickness distribution of ridging and ridged ice
-      !------------------------------------------------------------------------------!
-      CALL ice_rdgrft_ridgeprep
-
-      !------------------------------------------------------------------------------!
-      ! 3) Rothrock(1975)'s method
-      !------------------------------------------------------------------------------!
-      IF( kstrngth == 1 ) THEN
+      !                              !--------------------------------------------------!
+      CALL ice_rdgrft_ridgeprep      ! Thickness distribution of ridging and ridged ice !
+      !                              !--------------------------------------------------!
+
+      !                              !--------------------------------------------------!
+      IF( kstrngth == 1 ) THEN       ! Ice strength => Rothrock (1975) method           !
+      !                              !--------------------------------------------------!
          z1_3 = 1._wp / 3._wp
          DO jl = 1, jpl
-            DO jj= 1, jpj
-               DO ji = 1, jpi
-                  !
-                  IF( athorn(ji,jj,jl) > 0._wp ) THEN
-                     !----------------------------
-                     ! PE loss from deforming ice
-                     !----------------------------
-                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
-
-                     !--------------------------
-                     ! PE gain from rafting ice
-                     !--------------------------
-                     strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
-
-                     !----------------------------
-                     ! PE gain from ridging ice
-                     !----------------------------
-                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) * krdg(ji,jj,jl) * z1_3 *  &
-                        &                              ( hrmax(ji,jj,jl) * hrmax(ji,jj,jl) +         &
-                        &                                hrmin(ji,jj,jl) * hrmin(ji,jj,jl) +         &
-                        &                                hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )  
-                        !!(a**3-b**3)/(a-b) = a*a+ab+b*b                      
-                  ENDIF
-                  !
-               END DO
-            END DO
-         END DO
-   
-         strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1)
-                         ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
-         ksmooth = 1
-
-      !------------------------------------------------------------------------------!
-      ! 4) Hibler (1979)' method
-      !------------------------------------------------------------------------------!
-      ELSE                      ! kstrngth ne 1:  Hibler (1979) form
-         !
+            WHERE( athorn(:,:,jl) > 0._wp )
+               strength(:,:) =  -         athorn(:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE loss from deforming ice
+                  &             + 2._wp * araft (:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE gain from rafting ice
+                  &             +         aridge(:,:,jl) * krdg(:,:,jl) * z1_3 *   &        ! PE gain from ridging ice
+                  &                      ( hrmax(:,:,jl) * hrmax(:,:,jl) +         &
+                  &                        hrmin(:,:,jl) * hrmin(:,:,jl) +         &
+                  &                        hrmax(:,:,jl) * hrmin(:,:,jl) )
+            ELSEWHERE
+               strength(:,:) = 0._wp
+            END WHERE
+         END DO
+         strength(:,:) = rn_pe_rdg * zdrho * strength(:,:) / aksum(:,:) * tmask(:,:,1)
+                         ! where zdrho = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
+         ismooth = 1
+      !                              !--------------------------------------------------!
+      ELSE                           ! Ice strength => Hibler (1979) method             !
+      !                              !--------------------------------------------------!
          strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1)
          !
-         ksmooth = 1
-         !
-      ENDIF                     ! kstrngth
-      !
-      !------------------------------------------------------------------------------!
-      ! 5) Impact of brine volume
-      !------------------------------------------------------------------------------!
+         ismooth = 1
+         !
+      ENDIF
+      !                              !--------------------------------------------------!
+      IF( ln_icestr_bvf ) THEN       ! Impact of brine volume                           !
+      !                              !--------------------------------------------------!
       ! CAN BE REMOVED
-      IF( ln_icestr_bvf ) THEN
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -883,12 +807,9 @@
          END DO
       ENDIF
-      !
-      !------------------------------------------------------------------------------!
-      ! 6) Smoothing ice strength
-      !------------------------------------------------------------------------------!
-      SELECT CASE( ksmooth )
-      !                       !-------------------
-      CASE( 1 )               ! Spatial smoothing
-         !                    !-------------------
+      !                              !--------------------------------------------------!
+      SELECT CASE( ismooth )         ! Smoothing ice strength                           !
+      !                              !--------------------------------------------------!
+      !
+      CASE( 1 )               !--- Spatial smoothing
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
@@ -903,5 +824,5 @@
             END DO
          END DO
-
+         
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
@@ -911,7 +832,5 @@
          CALL lbc_lnk( strength, 'T', 1. )
          !
-         !                    !--------------------
-      CASE( 2 )               ! Temporal smoothing
-         !                    !--------------------
+      CASE( 2 )               !--- Temporal smoothing
          IF ( kt_ice == nit000 ) THEN
             zstrp1(:,:) = 0._wp
@@ -922,8 +841,8 @@
             DO ji = 2, jpim1
                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN 
-                  numts_rm = 1 ! number of time steps for the running mean
-                  IF ( zstrp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
-                  IF ( zstrp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
-                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / numts_rm
+                  itframe = 1 ! number of time steps for the running mean
+                  IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
+                  IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
+                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
                   zstrp2  (ji,jj) = zstrp1  (ji,jj)
                   zstrp1  (ji,jj) = strength(ji,jj)
@@ -932,5 +851,5 @@
             END DO
          END DO
-         CALL lbc_lnk( strength, 'T', 1. )      ! Boundary conditions
+         CALL lbc_lnk( strength, 'T', 1. )
          !
       END SELECT
@@ -970,20 +889,21 @@
       IF (lwp) THEN                          ! control print
          WRITE(numout,*)
-         WRITE(numout,*)'ice_rdgrft_init : ice parameters for mechanical ice redistribution '
-         WRITE(numout,*)'~~~~~~~~~~~~~~~'
-         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
-         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
-         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
-         WRITE(numout,*)'   Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
-         WRITE(numout,*)'   Ridging of ice sheets or not                            ln_ridging  = ', ln_ridging
-         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
-         WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
-         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
-         WRITE(numout,*)'   Fraction of pond volume conserved during ridging        rn_fpondrdg = ', rn_fpondrdg 
-         WRITE(numout,*)'   Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
-         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
-         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft  
-         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
-         WRITE(numout,*)'   Fraction of pond volume conserved during rafting        rn_fpondrft = ', rn_fpondrft 
+         WRITE(numout,*) 'ice_rdgrft_init : ice parameters for mechanical ice redistribution '
+         WRITE(numout,*) '~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namiceitdme'
+         WRITE(numout,*) '      Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
+         WRITE(numout,*) '      Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
+         WRITE(numout,*) '      Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
+         WRITE(numout,*) '      Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
+         WRITE(numout,*) '      Ridging of ice sheets or not                            ln_ridging  = ', ln_ridging
+         WRITE(numout,*) '      Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
+         WRITE(numout,*) '      Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
+         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
+         WRITE(numout,*) '      Fraction of pond volume conserved during ridging        rn_fpondrdg = ', rn_fpondrdg 
+         WRITE(numout,*) '      Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
+         WRITE(numout,*) '      Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
+         WRITE(numout,*) '      Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft  
+         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
+         WRITE(numout,*) '      Fraction of pond volume conserved during rafting        rn_fpondrft = ', rn_fpondrft 
       ENDIF
       !
