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 8506)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8512)
@@ -32,5 +32,5 @@
 
    PUBLIC   ice_rdgrft               ! called by ice_stp
-   PUBLIC   ice_rdgrft_icestrength   ! called by icerhg_evp
+   PUBLIC   ice_rdgrft_strength      ! called by icerhg_evp
    PUBLIC   ice_rdgrft_init          ! called by ice_stp
    PUBLIC   ice_rdgrft_alloc         ! called by ice_init 
@@ -51,19 +51,20 @@
    REAL(wp)            ::   zdrho               ! 
    !
-   ! ** namelist (namiceitdme) **
-   REAL(wp) ::   rn_cs            !: fraction of shearing energy contributing to ridging            
-   INTEGER  ::   nn_partfun       !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007)
-   REAL(wp) ::   rn_gstar         !: fractional area of young ice contributing to ridging
-   REAL(wp) ::   rn_astar         !: equivalent of G* for an exponential participation function
-   LOGICAL  ::   ln_ridging       !: ridging of ice or not                        
-   REAL(wp) ::   rn_hstar         !: thickness that determines the maximal thickness of ridged ice
-   REAL(wp) ::   rn_por_rdg       !: initial porosity of ridges (0.3 regular value)
-   REAL(wp) ::   rn_fsnowrdg      !: fractional snow loss to the ocean during ridging
-   REAL(wp) ::   rn_fpondrdg      !: fractional melt pond loss to the ocean during ridging
-   LOGICAL  ::   ln_rafting       !: rafting of ice or not                        
-   REAL(wp) ::   rn_hraft         !: threshold thickness (m) for rafting / ridging 
-   REAL(wp) ::   rn_craft         !: coefficient for smoothness of the hyperbolic tangent in rafting
-   REAL(wp) ::   rn_fsnowrft      !: fractional snow loss to the ocean during ridging
-   REAL(wp) ::   rn_fpondrft      !: fractional snow loss to the ocean during rafting
+   ! ** namelist (namice_rdgrft) **
+   REAL(wp) ::   rn_cs            ! fraction of shearing energy contributing to ridging            
+   LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975))
+   REAL(wp) ::   rn_gstar         !    fractional area of young ice contributing to ridging
+   LOGICAL  ::   ln_partf_exp     ! participation function exponential (Lipscomb et al. (2007))
+   REAL(wp) ::   rn_astar         !    equivalent of G* for an exponential participation function
+   LOGICAL  ::   ln_ridging       ! ridging of ice or not                        
+   REAL(wp) ::   rn_hstar         !    thickness that determines the maximal thickness of ridged ice
+   REAL(wp) ::   rn_porordg       !    initial porosity of ridges (0.3 regular value)
+   REAL(wp) ::   rn_fsnwrdg       !    fractional snow loss to the ocean during ridging
+   REAL(wp) ::   rn_fpndrdg       !    fractional pond loss to the ocean during ridging
+   LOGICAL  ::   ln_rafting       ! rafting of ice or not                        
+   REAL(wp) ::   rn_hraft         !    threshold thickness (m) for rafting / ridging 
+   REAL(wp) ::   rn_craft         !    coefficient for smoothness of the hyperbolic tangent in rafting
+   REAL(wp) ::   rn_fsnwrft       !    fractional snow loss to the ocean during rafting
+   REAL(wp) ::   rn_fpndrft       !    fractional pond loss to the ocean during rafting
    !
    !!----------------------------------------------------------------------
@@ -133,5 +134,5 @@
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*)'icerdgrft : ice ridging and rafting'
+         IF(lwp) WRITE(numout,*)'ice_rdgrft: ice ridging and rafting'
          IF(lwp) WRITE(numout,*)'~~~~~~~~~~'
          !
@@ -146,5 +147,5 @@
       !-----------------------------------------------------------------------------!
       !
-      CALL ice_rdgrft_ridgeprep                             ! prepare ridging
+      CALL ice_rdgrft_prep                             ! prepare ridging
       !
       DO jj = 1, jpj                                        ! Initialize arrays.
@@ -272,5 +273,5 @@
          !
          IF( iterate_ridging == 1 ) THEN
-            CALL ice_rdgrft_ridgeprep
+            CALL ice_rdgrft_prep
             IF( niter  >  nitermax ) THEN
                WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
@@ -297,7 +298,7 @@
 
 
-   SUBROUTINE ice_rdgrft_ridgeprep
+   SUBROUTINE ice_rdgrft_prep
       !!---------------------------------------------------------------------!
-      !!                ***  ROUTINE ice_rdgrft_ridgeprep ***
+      !!                ***  ROUTINE ice_rdgrft_prep ***
       !!
       !! ** Purpose :   preparation for ridging and strength calculations
@@ -352,5 +353,5 @@
       !-----------------------------------------------------------------
       !
-      IF( nn_partfun == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975)
+      IF( ln_partf_lin ) THEN          !--- Linear formulation (Thorndike et al., 1975)
          DO jl = 0, jpl    
             DO jj = 1, jpj 
@@ -369,5 +370,5 @@
          END DO
          !
-      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007)
+      ELSEIF( ln_partf_exp ) THEN      !--- Exponential, more stable formulation (Lipscomb et al, 2007)
          !                        
          zdummy = 1._wp / ( 1._wp - EXP(-z1_astar) )        ! precompute exponential terms using zGsum as a work array
@@ -453,10 +454,10 @@
       END DO
       !
-   END SUBROUTINE ice_rdgrft_ridgeprep
+   END SUBROUTINE ice_rdgrft_prep
 
 
    SUBROUTINE ice_rdgrft_ridgeshift( opning, closing_gross )
       !!----------------------------------------------------------------------
-      !!                ***  ROUTINE ice_rdgrft_icestrength ***
+      !!                ***  ROUTINE ice_rdgrft_strength ***
       !!
       !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
@@ -565,6 +566,6 @@
             !--------------------------------------------------------------------------
             vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij)
-            vrdg2(ij) = vrdg1(ij) * ( 1. + rn_por_rdg )
-            vsw  (ij) = vrdg1(ij) * rn_por_rdg
+            vrdg2(ij) = vrdg1(ij) * ( 1. + rn_porordg )
+            vsw  (ij) = vrdg1(ij) * rn_porordg
 
             vsrdg (ij) = v_s  (ji,jj,  jl1) * afrac(ij)
@@ -621,9 +622,9 @@
             !  During the next time step, thermo_rates will determine whether
             !  the ocean cools or new ice grows.
-            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
-
-            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
-               &                                - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2)
+            wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnwrdg )   & 
+               &                                      + rhosn * vsrft(ij) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice  ! fresh water source for ocean
+
+            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnwrdg )         & 
+               &                                - esrft(ij) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2)
 
             ! MV MP 2016
@@ -633,6 +634,6 @@
             !  Place part of the melt pond volume into the ocean. 
             IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN
-               wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( rhofw * vprdg(ij) * ( 1._wp - rn_fpondrdg )   & 
-               &                                 + rhofw * vprft(ij) * ( 1._wp - rn_fpondrft ) ) * r1_rdtice  ! fresh water source for ocean
+               wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( rhofw * vprdg(ij) * ( 1._wp - rn_fpndrdg )   & 
+               &                                 + rhofw * vprft(ij) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice  ! fresh water source for ocean
             ENDIF
             ! END MV MP 2016
@@ -719,14 +720,14 @@
                v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) )
                smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) )
-               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + ( vsrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
-                  &                                        vsrft (ij) * rn_fsnowrft * zswitch(ij) )
-               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
-                  &                                        esrft (ij) * rn_fsnowrft * zswitch(ij) )
+               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + ( vsrdg (ij) * rn_fsnwrdg * fvol(ij)  +  &
+                  &                                        vsrft (ij) * rn_fsnwrft * zswitch(ij) )
+               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnwrdg * fvol(ij)  +  &
+                  &                                        esrft (ij) * rn_fsnwrft * zswitch(ij) )
                ! MV MP 2016
                IF ( nn_pnd_scheme > 0 ) THEN
-                  v_ip (ji,jj,jl2) = v_ip(ji,jj,jl2) + (   vprdg (ij) * rn_fpondrdg * fvol   (ij)   &
-                     &                                   + vprft (ij) * rn_fpondrft * zswitch(ij)   )
-                  a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2) + (   aprdg2(ij) * rn_fpondrdg * farea         & 
-                     &                                   + aprft2(ij) * rn_fpondrft * zswitch(ji)   )
+                  v_ip (ji,jj,jl2) = v_ip(ji,jj,jl2) + (   vprdg (ij) * rn_fpndrdg * fvol   (ij)   &
+                     &                                   + vprft (ij) * rn_fpndrft * zswitch(ij)   )
+                  a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2) + (   aprdg2(ij) * rn_fpndrdg * farea         & 
+                     &                                   + aprft2(ij) * rn_fpndrft * zswitch(ji)   )
                ENDIF
                ! END MV MP 2016
@@ -748,7 +749,7 @@
 
 
-   SUBROUTINE ice_rdgrft_icestrength( kstrngth )
+   SUBROUTINE ice_rdgrft_strength
       !!----------------------------------------------------------------------
-      !!                ***  ROUTINE ice_rdgrft_icestrength ***
+      !!                ***  ROUTINE ice_rdgrft_strength ***
       !!
       !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
@@ -759,9 +760,5 @@
       !!              by ridging. Note that only Hibler's formulation is stable and that
       !!              ice strength has to be smoothed
-      !!
-      !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979)
-      !
       INTEGER             ::   ji,jj, jl   ! dummy loop indices
       INTEGER             ::   ismooth     ! smoothing the resistance to deformation
@@ -773,9 +770,9 @@
 
       !                              !--------------------------------------------------!
-      CALL ice_rdgrft_ridgeprep      ! Thickness distribution of ridging and ridged ice !
+      CALL ice_rdgrft_prep           ! Thickness distribution of ridging and ridged ice !
       !                              !--------------------------------------------------!
 
       !                              !--------------------------------------------------!
-      IF( kstrngth == 1 ) THEN       ! Ice strength => Rothrock (1975) method           !
+      IF( ln_str_Rot ) THEN          ! Ice strength => Rothrock (1975) method           !
       !                              !--------------------------------------------------!
          z1_3 = 1._wp / 3._wp
@@ -792,9 +789,9 @@
             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
+         strength(:,:) = rn_perdg * zdrho * strength(:,:) / aksum(:,:) * tmask(:,:,1)
+                         ! where zdrho = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_perdg accounts for frictional dissipation
          ismooth = 1
       !                              !--------------------------------------------------!
-      ELSE                           ! Ice strength => Hibler (1979) method             !
+      ELSEIF( ln_str_Hib ) THEN      ! Ice strength => Hibler (1979) method             !
       !                              !--------------------------------------------------!
          strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1)
@@ -802,19 +799,8 @@
          ismooth = 1
          !
-      ENDIF
-      !                              !--------------------------------------------------!
-      IF( ln_icestr_bvf ) THEN       ! Impact of brine volume                           !
-      !                              !--------------------------------------------------!
-      ! CAN BE REMOVED
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0)))
-            END DO
-         END DO
       ENDIF
       !                              !--------------------------------------------------!
       SELECT CASE( ismooth )         ! Smoothing ice strength                           !
       !                              !--------------------------------------------------!
-      !
       CASE( 1 )               !--- Spatial smoothing
          DO jj = 2, jpjm1
@@ -861,5 +847,5 @@
       END SELECT
       !
-   END SUBROUTINE ice_rdgrft_icestrength
+   END SUBROUTINE ice_rdgrft_strength
 
 
@@ -871,48 +857,66 @@
       !!                to the mechanical ice redistribution
       !!
-      !! ** Method  :   Read the namiceitdme namelist 
+      !! ** Method  :   Read the namice_rdgrft namelist 
       !!                and check the parameters values 
       !!                called at the first timestep (nit000)
       !!
-      !! ** input   :   Namelist namiceitdme
+      !! ** input   :   Namelist namice_rdgrft
       !!-------------------------------------------------------------------
       INTEGER :: ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namiceitdme/ rn_cs     , nn_partfun, rn_gstar  , rn_astar   ,                & 
-         &                  ln_ridging, rn_hstar  , rn_por_rdg, rn_fsnowrdg, rn_fpondrdg,   & 
-         &                  ln_rafting, rn_hraft  , rn_craft  , rn_fsnowrft, rn_fpondrft
+      NAMELIST/namice_rdgrft/ ln_str_Hib, rn_pstar, rn_crhg, &
+         &                    ln_str_Rot, rn_perdg,          &
+         &                    rn_cs     ,                    &
+         &                    ln_partf_lin, rn_gstar,        &
+         &                    ln_partf_exp, rn_astar,        & 
+         &                    ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg,  & 
+         &                    ln_rafting, rn_hraft, rn_craft  , rn_fsnwrft, rn_fpndrft
       !!-------------------------------------------------------------------
       !
       REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
-      READ  ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
-      !
-      REWIND( numnam_ice_cfg )              ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
-      READ  ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namiceitdme )
+      READ  ( numnam_ice_ref, namice_rdgrft, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_rdgrft in reference namelist', lwp )
+      !
+      REWIND( numnam_ice_cfg )              ! Namelist namice_rdgrft in configuration namelist : Ice mechanical ice redistribution
+      READ  ( numnam_ice_cfg, namice_rdgrft, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_rdgrft in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namice_rdgrft )
       !
       IF (lwp) THEN                          ! control print
          WRITE(numout,*)
-         WRITE(numout,*) 'ice_rdgrft_init : ice parameters for mechanical ice redistribution '
+         WRITE(numout,*) 'ice_rdgrft_init : ice parameters for ridging/rafting '
          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 
+         WRITE(numout,*) '   Namelist namice_rdgrft'
+         WRITE(numout,*) '      ice strength parameterization Hibler (1979)              ln_str_Hib   = ', ln_str_Hib 
+         WRITE(numout,*) '            1st bulk-rheology parameter                        rn_pstar     = ', rn_pstar
+         WRITE(numout,*) '            2nd bulk-rhelogy parameter                         rn_crhg      = ', rn_crhg
+         WRITE(numout,*) '      ice strength parameterization Rothrock (1975)            ln_str_Rot   = ', ln_str_Rot 
+         WRITE(numout,*) '            Ratio of ridging work to PotEner change in ridging rn_perdg     = ', rn_perdg 
+         WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_cs        = ', rn_cs 
+         WRITE(numout,*) '      linear ridging participation function                    ln_partf_lin = ', ln_partf_lin
+         WRITE(numout,*) '            Fraction of ice coverage contributing to ridging   rn_gstar     = ', rn_gstar
+         WRITE(numout,*) '      Exponential ridging participation function               ln_partf_exp = ', ln_partf_exp
+         WRITE(numout,*) '            Equivalent to G* for an exponential function       rn_astar     = ', rn_astar
+         WRITE(numout,*) '      Ridging of ice sheets or not                             ln_ridging   = ', ln_ridging
+         WRITE(numout,*) '            max ridged ice thickness                           rn_hstar     = ', rn_hstar
+         WRITE(numout,*) '            Initial porosity of ridges                         rn_porordg   = ', rn_porordg
+         WRITE(numout,*) '            Fraction of snow volume conserved during ridging   rn_fsnwrdg   = ', rn_fsnwrdg 
+         WRITE(numout,*) '            Fraction of pond volume conserved during ridging   rn_fpndrdg   = ', rn_fpndrdg 
+         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 rafting   rn_fsnwrft   = ', rn_fsnwrft 
+         WRITE(numout,*) '            Fraction of pond volume conserved during rafting   rn_fpndrft   = ', rn_fpndrft 
       ENDIF
       !
-   END SUBROUTINE ice_rdgrft_init
+      IF ( ( ln_str_Hib .AND. ln_str_Rot ) .OR. ( .NOT.ln_str_Hib .AND. .NOT.ln_str_Rot ) ) THEN
+         CALL ctl_stop( 'ice_rdgrft_init: choose one and only one formulation for ice strength (ln_str_Hib or ln_str_Rot)' )
+      ENDIF
+      !
+      IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
+         CALL ctl_stop( 'ice_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
+      ENDIF
+      !
+  END SUBROUTINE ice_rdgrft_init
 
 #else
