Index: /NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE/icedyn_rdgrft.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE/icedyn_rdgrft.F90	(revision 15734)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE/icedyn_rdgrft.F90	(revision 15735)
@@ -74,4 +74,5 @@
    LOGICAL  ::   ln_redist_ufm    ! uniform redistribution of ridged ice (Hibler, 1980) 
    LOGICAL  ::   ln_redist_exp    ! exponential redistribution of ridged ice
+   REAL(wp) ::   max_strength     !    maximum strength cap   
    REAL(wp) ::   rn_murdg         !    gives e-folding scale of ridged ice (m^.5) for ln_redist_exp
    REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging            
@@ -513,6 +514,8 @@
       ! 2) Transfer function
       !-----------------------------------------------------------------
+      ! If assuming ridged ice is uniformly distributed between hrmin and 
+      ! hrmax (ln_redist_ufm):
+      !
       ! Compute max and min ridged ice thickness for each ridging category.
-      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
       ! 
       ! This parameterization is a modified version of Hibler (1980).
@@ -528,4 +531,32 @@
       ! These modifications have the effect of reducing the ice strength
       ! (relative to the Hibler formulation) when very thick ice is ridging.
+      !
+      !-----------------------------------------------------------------
+      ! If assuming ridged ice ITD is a negative exponential 
+      ! (ln_redist_exp) and following CICE implementation: 
+      ! 
+      !  g(h) ~ exp[-(h-hrmin)/hrexp], h >= hrmin 
+      ! 
+      ! where hrmin is the minimum thickness of ridging ice and 
+      ! hrexp is the e-folding thickness.
+      ! 
+      ! Here, assume as above that hrmin = min(2*hi, hi+maxraft).
+      ! That is, the minimum ridge thickness results from rafting,
+      !  unless the ice is thicker than maxraft.
+      !
+      ! Also, assume that hrexp = mu_rdg*sqrt(hi).
+      ! The parameter mu_rdg is tuned to give e-folding scales mostly
+      !  in the range 2-4 m as observed by upward-looking sonar.
+      !
+      ! Values of mu_rdg in the right column give ice strengths
+      !  roughly equal to values of Hstar in the left column
+      !  (within ~10 kN/m for typical ITDs):
+      !
+      !   Hstar     mu_rdg
+      !
+      !     25        3.0
+      !     50        4.0
+      !     75        5.0
+      !    100        6.0
       !
       ! zaksum = net area removed/ total area participating
@@ -792,5 +823,5 @@
                  ! Compute the fraction of ridged ice area and volume going to thickness category jl2
                  !
-                 IF ( ln_redist_ufm ) THEN ! Hibler 1980 uniform formulation
+                 IF ( ln_redist_ufm ) THEN ! Hibler (1980) uniform formulation
                  !
                   IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
@@ -806,5 +837,5 @@
                   ENDIF
                   !
-                 ELSE IF ( ln_redist_exp ) THEN ! exponential formulation
+                 ELSE IF ( ln_redist_exp ) THEN ! Lipscomb et al. (2007) exponential formulation
                   !
                   IF (jl2 < jpl) THEN
@@ -929,14 +960,13 @@
       IF( ln_str_H79 ) THEN          ! Ice strength => Hibler (1979) method             !
       !                              !--------------------------------------------------!
+      ! Recommend using ismooth = 1 (spatial smoothing)
         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
-      ! Recommend ismooth = 1
       !                              !--------------------------------------------------!
       ELSE IF( ln_str_R75 ) THEN     ! Ice strength => Rothrock (1975) method           !
       !                              !--------------------------------------------------!
-      !
-        ! Initialise the strength field at each timestep
+      ! Recommend using ismooth = 0 (no smoothing needed)
+      ! Initialise the strength field at each timestep
         strength(:,:) = 0._wp
-      ! Recommend ismooth = 0
-
+      !
         !--------------------------------
         ! 0) Identify grid cells with ice
@@ -979,9 +1009,9 @@
               IF ( apartf(ji,jl) > 0._wp ) THEN
                 !
-                IF ( ln_redist_ufm ) THEN   ! Uniform distribution of ridged ice               
+                IF ( ln_redist_ufm ) THEN   ! Uniform redistribution of ridged ice               
                   h2rdg = (1._wp/3._wp) * (hrmax(ji,jl)**3._wp - hrmin(ji,jl)**3._wp)  &
                            / (hrmax(ji,jl) - hrmin(ji,jl))
                 !
-                ELSE IF ( ln_redist_exp ) THEN   ! Exponential distribution of ridged ice
+                ELSE IF ( ln_redist_exp ) THEN   ! Exponential redistribution of ridged ice
 	 	  h2rdg = hrmin(ji,jl) * hrmin(ji,jl) & 
 	 	          + 2._wp * hrmin(ji,jl) * hrexp(ji,jl) & 
@@ -1004,5 +1034,5 @@
 
           DO ji = 1, npti
-            IF ( zaksum(ji) > epsi10 ) THEN  ! Unclear why is this check should be needed? Occasional x10E-24...
+            IF ( zaksum(ji) > epsi10 ) THEN
               strength_1d(ji) = rn_Cf * Cp * strength_pe_sum(ji) / zaksum(ji)
             ELSE
@@ -1013,5 +1043,6 @@
 
           ! Enforce a maximum for strength
-          WHERE( strength_1d(:) > 2.0e5) ; strength_1d(:) = 2.0e5
+          ! Richter-Menge and Elder (1998) estimate maximum in Beaufort Sea in wintertime of the order 150 kN/m
+          WHERE( strength_1d(:) > max_strength ) ; strength_1d(:) = max_strength
           END WHERE
 
@@ -1035,5 +1066,5 @@
       CASE( 0 )               !--- No smoothing
 
-      CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
+         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
 
       CASE( 1 )               !--- Spatial smoothing
@@ -1174,7 +1205,8 @@
       INTEGER :: ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_rdgrft/ ismooth, ln_str_H79, rn_pstar,          &
-         &                    ln_str_R75, rn_Cf,                      &
-         &                    ln_redist_ufm, ln_redist_exp, rn_murdg, &
+      NAMELIST/namdyn_rdgrft/ ismooth, ln_str_H79, rn_pstar, &
+         &                    ln_str_R75, rn_Cf,             &
+         &                    ln_redist_ufm, ln_redist_exp,  &
+         &                    max_strength, rn_murdg,        &
          &                    rn_crhg, rn_csrdg,             &
          &                    ln_partf_lin, rn_gstar,        &
@@ -1203,4 +1235,5 @@
          WRITE(numout,*) '            ratio of ridging work to PE change in ridging      rn_Cf        = ', rn_Cf
          WRITE(numout,*) '      smoothing the resistance to deformation (strength)       ismooth      = ', ismooth
+         WRITE(numout,*) '      maximum strength cap                                     max_strength = ', max_strength
          WRITE(numout,*) '      uniform redistribution of ridged ice (Hibler, 1980)     ln_redist_ufm = ', ln_redist_ufm
          WRITE(numout,*) '      exponential redistribution of ridged ice                ln_redist_exp = ', ln_redist_exp
@@ -1235,6 +1268,9 @@
       ENDIF
       !
-      !! Add check for ismooth too, =0, 1 or 2
-
+      IF ( ( ismooth .NE. 0 ) .AND. ( ismooth .NE. 1 ) .AND. (ismooth .NE. 2 ) ) THEN
+         CALL ctl_stop( 'ice_dyn_rdgrft_init: ismooth must be one of 0 (no smoothing), 1 (spatial), or 2 (temporal)' )
+      ENDIF
+      !
+      !
       IF( .NOT. ln_icethd ) THEN
          rn_porordg = 0._wp
Index: /NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE/icedyn_rhg_evp.F90	(revision 15734)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE/icedyn_rhg_evp.F90	(revision 15735)
@@ -792,5 +792,5 @@
       END DO
       CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., &
-         &                                  zs1     , 'T', 1., zs2    , 'T', 1., zs12    , 'F', 1. )
+         &                                  zs1     , 'T', 1., zs2    , 'T', 1., zs12    , 'F', 1., zshear, 'T', 1. )
       
       ! --- Store the stress tensor for the next time step --- !
@@ -885,6 +885,6 @@
          END DO               
          !
-         CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00(:,:) ) 
-         CALL iom_put( 'sig2_pnorm' , zsig2_p(:,:) * zmsk00(:,:) ) 
+         IF( iom_use('sig1_pnorm') )  CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00(:,:) )  
+         IF( iom_use('sig2_pnorm') )  CALL iom_put( 'sig2_pnorm' , zsig2_p(:,:) * zmsk00(:,:) ) 
       
          DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
@@ -1009,5 +1009,5 @@
 
       ! time
-      it = ( kt - 1 ) * kitermax + kiter
+      it = ( kt - nit000 ) * kitermax + kiter
       
       ! convergence
@@ -1029,5 +1029,5 @@
          istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) )
          ! close file
-         IF( kt == nitend - nn_fsbc + 1 )   istatus = NF90_CLOSE(ncvgid)
+         IF( kt == nitend - nn_fsbc + 1 .AND. kiter == kitermax )   istatus = NF90_CLOSE(ncvgid)
       ENDIF
       
