Index: NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rdgrft.F90	(revision 15570)
+++ NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rdgrft.F90	(revision 15742)
@@ -43,7 +43,11 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   opning          ! rate of opening due to divergence/shear
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   closing_gross   ! rate at which area removed, not counting area of new ridges
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   zaksum          ! normalisation factor
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   zasum           ! sum of ice area plus open water   
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   apartf          ! participation function; fraction of ridging/closing associated w/ category n
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   zhi             ! ice thickness in category n
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrmin           ! minimum ridge thickness
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrmax           ! maximum ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrexp           ! e-folding ridge thickness
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hraft           ! thickness of rafted ice
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hi_hrdg         ! thickness of ridging ice / mean ridge thickness
@@ -53,10 +57,23 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_s_2d
    !
+   ! For ridging diagnostics
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   airdg1          ! Ridging ice area loss
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   airft1          ! Rafting ice area loss
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   airdg2          ! New ridged ice area gain
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   airft2          ! New rafted ice area gain
+   !
    REAL(wp), PARAMETER ::   hrdg_hi_min = 1.1_wp    ! min ridge thickness multiplier: min(hrdg/hi)
-   REAL(wp), PARAMETER ::   hi_hrft     = 0.5_wp    ! rafting multipliyer: (hi/hraft)
+   REAL(wp), PARAMETER ::   hi_hrft     = 0.5_wp    ! rafting multiplier: (hi/hraft)
    !
    ! ** namelist (namdyn_rdgrft) **
    LOGICAL  ::   ln_str_H79       ! ice strength parameterization (Hibler79)
-   REAL(wp) ::   rn_pstar         ! determines ice strength, Hibler JPO79
+   REAL(wp) ::   rn_pstar         !    determines ice strength, Hibler JPO79
+   LOGICAL  ::   ln_str_R75       ! ice strength parameterization (Rothrock75)
+   REAL(wp) ::   rn_Cf            !    ratio of ridging work to PE change in ridging, R75
+   INTEGER  ::   ismooth          ! smoothing the resistance to deformation (strength)
+   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            
    LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975))
@@ -86,8 +103,10 @@
       !!                ***  ROUTINE ice_dyn_rdgrft_alloc ***
       !!-------------------------------------------------------------------
-      ALLOCATE( closing_net(jpij)  , opning(jpij)      , closing_gross(jpij) ,               &
-         &      apartf(jpij,0:jpl) , hrmin  (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), &
-         &      hrmax (jpij,jpl)   , hi_hrdg(jpij,jpl) , araft(jpij,jpl) ,                   &
-         &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc )
+      ALLOCATE( closing_net(jpij)  , opning(jpij)      , closing_gross(jpij) ,  &
+         &      apartf(jpij,0:jpl) , hrmin(jpij,jpl)   , hrmax(jpij,jpl)     , zhi(jpij,jpl),  &
+         &      zaksum(jpij)       , hraft(jpij,jpl)   , hrexp(jpij,jpl)     ,  &
+         &      hi_hrdg(jpij,jpl)  , araft(jpij,jpl)   , aridge(jpij,jpl)    , zasum(jpij), &
+         &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), &
+         &      airdg1(jpij), airft1(jpij), airdg2(jpij), airft2(jpij), STAT=ice_dyn_rdgrft_alloc )
 
       CALL mpp_sum ( 'icedyn_rdgrft', ice_dyn_rdgrft_alloc )
@@ -139,4 +158,15 @@
       REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i
       REAL(wp), DIMENSION(jpij) ::   zconv         ! 1D rdg_conv (if EAP rheology)
+
+      REAL(wp), DIMENSION(jpi,jpj) ::   opning_2d  ! Lead opening diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   airdg1_2d  ! Ridging ice area loss diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   airft1_2d  ! Rafting ice area loss diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   airdg2_2d  ! New ridged ice area gain diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   airft2_2d  ! New rafted ice area gain diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   dairdg1dt  ! Ridging ice area loss rate diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   dairft1dt  ! Rafting ice area loss rate diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   dairdg2dt  ! New ridged ice area gain rate diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   dairft2dt  ! New rafted ice area gain rate diagnostic
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00     ! Mask for ice presence
       !
       INTEGER, PARAMETER ::   jp_itermax = 20    
@@ -151,5 +181,16 @@
          IF(lwp) WRITE(numout,*)'ice_dyn_rdgrft: ice ridging and rafting'
          IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~~'
-      ENDIF      
+      ENDIF     
+
+      ! Initialise ridging diagnostics
+      opning_2d(:,:) = 0.0_wp
+      dairdg1dt(:,:) = 0.0_wp
+      dairft1dt(:,:) = 0.0_wp
+      dairdg2dt(:,:) = 0.0_wp
+      dairft2dt(:,:) = 0.0_wp
+      airdg1_2d(:,:) = 0.0_wp
+      airft1_2d(:,:) = 0.0_wp
+      airdg2_2d(:,:) = 0.0_wp
+      airft2_2d(:,:) = 0.0_wp
 
       !--------------------------------
@@ -162,4 +203,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
+            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice, 0 if no ice. epsi06 used here rather than epsi10 as below for consistency with ice masking elsewhere
             IF ( at_i(ji,jj) > epsi10 ) THEN
                npti           = npti + 1
@@ -188,9 +230,10 @@
             IF( ln_rhg_EVP )  closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
             IF( ln_rhg_EAP )  closing_net(ji) = zconv(ji)
-            !
+            !     
             IF( zdivu(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) )   ! make sure the closing rate is large enough
             !                                                                                ! to give asum = 1.0 after ridging
             ! Opening rate (non-negative) that will give asum = 1.0 after ridging.
             opning(ji) = closing_net(ji) + zdivu(ji)
+            !
          END DO
          !
@@ -266,6 +309,60 @@
          CALL ice_dyn_1d2d( 2 )            ! --- Move to 2D arrays --- !
 
+! --- Ridging diagnostics --- !
+
+         IF( iom_use('opning') ) THEN
+           CALL tab_1d_2d( npti, nptidx(1:npti), opning(1:npti), opning_2d(:,:) )
+         END IF
+
+         IF( iom_use('dairdg1dt') .OR. iom_use('dairft1dt') .OR. iom_use('dairdg2dt') .OR. iom_use('dairft2dt') ) THEN
+           !
+           CALL tab_1d_2d( npti, nptidx(1:npti), airdg1(1:npti), airdg1_2d(:,:) )
+           CALL tab_1d_2d( npti, nptidx(1:npti), airft1(1:npti), airft1_2d(:,:) )
+           CALL tab_1d_2d( npti, nptidx(1:npti), airdg2(1:npti), airdg2_2d(:,:) )
+           CALL tab_1d_2d( npti, nptidx(1:npti), airft2(1:npti), airft2_2d(:,:) )
+           !
+           DO jj = 1, jpj
+              DO ji = 1, jpi
+              !
+              dairdg1dt(ji,jj) = airdg1_2d(ji,jj) * r1_rdtice
+              dairft1dt(ji,jj) = airft1_2d(ji,jj) * r1_rdtice                      
+              dairdg2dt(ji,jj) = airdg2_2d(ji,jj) * r1_rdtice
+              dairft2dt(ji,jj) = airft2_2d(ji,jj) * r1_rdtice 
+              !
+              END DO
+           END DO    
+          !         
+         ENDIF
+        ! 
+      ENDIF  ! npti>0
+
+
+      IF( iom_use('opning') ) THEN
+        CALL lbc_lnk( 'icedyn_rdgrft', opning_2d(:,:), 'T', 1. )
+        CALL iom_put( 'opning', opning_2d(:,:) * zmsk00(:,:) )     ! Lead area opening rate
       ENDIF
-   
+
+      IF( iom_use('dairdg1dt') ) THEN
+        CALL lbc_lnk( 'icedyn_rdgrft', dairdg1dt(:,:), 'T', 1. )
+        CALL iom_put( 'dairdg1dt', dairdg1dt(:,:) * zmsk00(:,:) )  ! Ridging ice area loss rate
+      ENDIF
+
+      IF( iom_use('dairft1dt') ) THEN
+        CALL lbc_lnk( 'icedyn_rdgrft', dairft1dt(:,:), 'T', 1. )
+        CALL iom_put( 'dairft1dt', dairft1dt(:,:) * zmsk00(:,:) )  ! Rafting ice area loss rate
+      ENDIF
+
+      IF( iom_use('dairdg2dt') ) THEN
+        CALL lbc_lnk( 'icedyn_rdgrft', dairdg2dt(:,:), 'T', 1. )
+        CALL iom_put( 'dairdg2dt', dairdg2dt(:,:) * zmsk00(:,:) )  ! New ridged ice area gain rate
+      ENDIF
+
+      IF( iom_use('dairft2dt') ) THEN
+        CALL lbc_lnk( 'icedyn_rdgrft', dairft2dt(:,:), 'T', 1. )
+        CALL iom_put( 'dairft2dt', dairft2dt(:,:) * zmsk00(:,:) )  ! New rafted ice area gain rate
+      ENDIF
+
+! ------- !
+
       CALL ice_var_agg( 1 ) 
 
@@ -280,5 +377,6 @@
 
 
-   SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i, pclosing_net )
+
+   SUBROUTINE rdgrft_prep( pa_i_in, pv_i, pato_i, pclosing_net )
       !!-------------------------------------------------------------------
       !!                ***  ROUTINE rdgrft_prep ***
@@ -290,10 +388,10 @@
       !!-------------------------------------------------------------------
       REAL(wp), DIMENSION(:)  , INTENT(in) ::   pato_i, pclosing_net 
-      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pa_i, pv_i 
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pa_i_in, pv_i 
+      REAL(wp), DIMENSION(jpij,jpl)        ::   pa_i                     ! area adjusted for minimum zhi
       !!
       INTEGER  ::   ji, jl                     ! dummy loop indices
       REAL(wp) ::   z1_gstar, z1_astar, zhmean, zfac   ! local scalar
-      REAL(wp), DIMENSION(jpij)        ::   zasum, z1_asum, zaksum   ! sum of a_i+ato_i and reverse 
-      REAL(wp), DIMENSION(jpij,jpl)    ::   zhi                      ! ice thickness
+      REAL(wp), DIMENSION(jpij)        ::   z1_asum                  ! inverse of sum of a_i+ato_i
       REAL(wp), DIMENSION(jpij,-1:jpl) ::   zGsum                    ! zGsum(n) = sum of areas in categories 0 to n
       !--------------------------------------------------------------------
@@ -301,4 +399,6 @@
       z1_gstar = 1._wp / rn_gstar
       z1_astar = 1._wp / rn_astar
+
+      pa_i(:,:) = pa_i_in(:,:)
 
       !                       ! Ice thickness needed for rafting
@@ -306,4 +406,16 @@
       ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp
       END WHERE
+      ! Make sure ice thickness is not below the minimum (from iceitd.F90; no minimum was causing issues for R75 
+      ! strength)
+      WHERE( pa_i(1:npti,:) > epsi10 .AND. zhi(1:npti,:) < rn_himin ) 
+        pa_i(1:npti,:) = pa_i(1:npti,:) * zhi(1:npti,:) / rn_himin
+        zhi(1:npti,:) = rn_himin
+      END WHERE
+
+      ! Set tiny values of pa_i to zero
+      ! If this isn't done, can end up with zhi=0 but apartf>0
+      WHERE( pa_i(1:npti,:) <= epsi10 )   ;   pa_i(1:npti,:) = 0._wp
+      END WHERE
+
 
       ! 1) Participation function (apartf): a(h) = b(h).g(h)
@@ -320,10 +432,10 @@
       ! Compute total area of ice plus open water.
       ! This is in general not equal to one because of divergence during transport
-      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
-      !
+      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,1:jpl), dim=2 )
+
       WHERE( zasum(1:npti) > epsi10 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
       ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp
       END WHERE
-      !
+
       ! Compute cumulative thickness distribution function
       ! Compute the cumulative thickness distribution function zGsum,
@@ -333,5 +445,5 @@
       zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti)
       DO jl = 1, jpl
-         zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti)  ! sum(1:jl) is ok (and not jpl)
+         zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti)  ! sum(1:jl) is correct (and not jpl) 
       END DO
       !
@@ -398,8 +510,12 @@
       ENDIF
 
+
+
       ! 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).
@@ -416,7 +532,35 @@
       ! (relative to the Hibler formulation) when very thick ice is ridging.
       !
-      ! zaksum = net area removed/ total area removed
-      ! where total area removed = area of ice that ridges
-      !         net area removed = total area removed - area of new ridges
+      !-----------------------------------------------------------------
+      ! 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
+      ! where total area participating = area of ice that ridges
+      !         net area removed = total area participating - area of new ridges
       !-----------------------------------------------------------------
       zfac = 1._wp / hi_hrft
@@ -428,7 +572,13 @@
                zhmean         = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min )
                hrmin  (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) )
-               hrmax  (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
                hraft  (ji,jl) = zhi(ji,jl) * zfac
-               hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 )
+               !
+               IF ( ln_redist_ufm ) THEN
+                 hrmax  (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
+                 hi_hrdg(ji,jl) = zhi(ji,jl) / zhmean
+               ELSE IF ( ln_redist_exp ) THEN
+                 hrexp  (ji,jl) = rn_murdg * SQRT( zhi(ji,jl) )               
+                 hi_hrdg(ji,jl) = zhi(ji,jl) / ( hrmin(ji,jl) + hrexp(ji,jl) )     
+               END IF    
                !
                ! Normalization factor : zaksum, ensures mass conservation
@@ -439,5 +589,7 @@
                hrmax  (ji,jl) = 0._wp 
                hraft  (ji,jl) = 0._wp 
+               hrexp  (ji,jl) = 0._wp
                hi_hrdg(ji,jl) = 1._wp
+               zaksum (ji)    = 0._wp
             ENDIF
          END DO
@@ -493,10 +645,11 @@
       INTEGER  ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
       REAL(wp) ::   hL, hR, farea              ! left and right limits of integration and new area going to jl2
+      REAL(wp) ::   expL, expR                 ! exponentials involving hL, hR
       REAL(wp) ::   vsw                        ! vol of water trapped into ridges
       REAL(wp) ::   afrdg, afrft               ! fraction of category area ridged/rafted 
-      REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1
-      REAL(wp)                  ::   airft1, oirft1, aprft1
-      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg  ! area etc of new ridges
-      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft  ! area etc of rafted ice
+      REAL(wp)                  ::   oirdg1, aprdg1, virdg1, sirdg1
+      REAL(wp)                  ::   oirft1, aprft1
+      REAL(wp), DIMENSION(jpij) ::   oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg  ! area etc of new ridges
+      REAL(wp), DIMENSION(jpij) ::   oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft  ! area etc of rafted ice
       !
       REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges
@@ -538,13 +691,13 @@
                
                ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2)
-               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
-               airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice
-
-               airdg2(ji) = airdg1 * hi_hrdg(ji,jl1)
-               airft2(ji) = airft1 * hi_hrft
+               airdg1(ji) = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
+               airft1(ji) = araft (ji,jl1) * closing_gross(ji) * rdt_ice
+
+               airdg2(ji) = airdg1(ji) * hi_hrdg(ji,jl1)
+               airft2(ji) = airft1(ji) * hi_hrft
 
                ! ridging /rafting fractions
-               afrdg = airdg1 * z1_ai(ji)
-               afrft = airft1 * z1_ai(ji)
+               afrdg = airdg1(ji) * z1_ai(ji)
+               afrft = airft1(ji) * z1_ai(ji)
 
                ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
@@ -602,5 +755,5 @@
                ! Remove area, volume of new ridge to each category jl1
                !------------------------------------------------------
-               a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1    - airft1
+               a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1(ji)- airft1(ji)
                v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1    - virft(ji)
                v_s_2d (ji,jl1) = v_s_2d (ji,jl1) - vsrdg(ji) - vsrft(ji)
@@ -661,11 +814,15 @@
          itest_rdg(1:npti) = 0
          itest_rft(1:npti) = 0
+         !
          DO jl2  = 1, jpl 
             !
             DO ji = 1, npti
 
-               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
-
-                  ! Compute the fraction of ridged ice area and volume going to thickness category jl2
+              IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
+                 !
+                 ! Compute the fraction of ridged ice area and volume going to thickness category jl2
+                 !
+                 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
                      hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
@@ -679,4 +836,32 @@
                      fvol(ji) = 0._wp                  
                   ENDIF
+                  !
+                 ELSE IF ( ln_redist_exp ) THEN ! Lipscomb et al. (2007) exponential formulation
+                  !
+                  IF (jl2 < jpl) THEN
+                  !
+                    IF (hrmin(ji,jl1) <= hi_max(jl2)) THEN
+                       hL       = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
+                       hR       = hi_max(jl2)
+                       expL     = EXP( -( hL - hrmin(ji,jl1) ) / hrexp(ji,jl1) )
+                       expR     = EXP( -( hR - hrmin(ji,jl1) ) / hrexp(ji,jl1) )
+                       farea    = expL - expR
+                       fvol(ji) = ( ( hL + hrexp(ji,jl1) ) * expL  &
+                                  - ( hR + hrexp(ji,jl1) ) * expR ) / ( hrmin(ji,jl1) + hrexp(ji,jl1) )
+                     ELSE
+                       farea    = 0._wp 
+                       fvol(ji) = 0._wp  
+                     END IF             
+                     !                 
+                  ELSE             ! jl2 = jpl
+                    !
+                    hL       = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
+                    expL     = EXP(-( hL - hrmin(ji,jl1) ) / hrexp(ji,jl1) )
+                    farea    = expL
+                    fvol(ji) = ( hL + hrexp(ji,jl1) ) * expL / ( hrmin(ji,jl1) + hrexp(ji,jl1) )  
+                    !
+                  END IF            ! jl2 < jpl
+                  ! 
+                 END IF             ! ridge redistribution
 
                   ! Compute the fraction of rafted ice area and volume going to thickness category jl2
@@ -757,28 +942,130 @@
       !!              dissipated per unit area removed from the ice pack under compression,
       !!              and assumed proportional to the change in potential energy caused
-      !!              by ridging. Note that only Hibler's formulation is stable and that
-      !!              ice strength has to be smoothed
+      !!              by ridging. Note that ice strength using Hibler's formulation must be
+      !!              smoothed.
       !!----------------------------------------------------------------------
       INTEGER             ::   ji, jj, jl  ! dummy loop indices
-      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)            ::   Cp          ! proportional constant for PE
+      REAL(wp)            ::   h2rdg       ! mean value of h^2 for new ridge
+      REAL(wp)            ::   dh2rdg      ! change in mean value of h^2 per unit area consumed by ridging 
       REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
       REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
+      REAL(wp), DIMENSION(jpij)     ::   strength_1d      ! 1-dimensional strength
+      REAL(wp), DIMENSION(jpij,jpl) ::   strength_pe      ! PE component of R75 strength
+      REAL(wp), DIMENSION(jpij)     ::   strength_pe_sum  ! PE component of R75 strength 
       !!----------------------------------------------------------------------
       !                              !--------------------------------------------------!
       IF( ln_str_H79 ) THEN          ! Ice strength => Hibler (1979) method             !
       !                              !--------------------------------------------------!
-         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
-         ismooth = 1
+      ! Recommend using ismooth = 1 (spatial smoothing)
+        strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
+      !                              !--------------------------------------------------!
+      ELSE IF( ln_str_R75 ) THEN     ! Ice strength => Rothrock (1975) method           !
+      !                              !--------------------------------------------------!
+      ! Recommend using ismooth = 0 (no smoothing needed)
+      ! Initialise the strength field at each timestep
+        strength(:,:) = 0._wp
+      !
+        !--------------------------------
+        ! 0) Identify grid cells with ice
+        !--------------------------------
+        at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
+        vt_i(:,:) = SUM( v_i(:,:,:), dim=3 )
+        ato_i(:,:) = 1 - at_i(:,:)
+        !
+        npti = 0   ;   nptidx(:) = 0
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+              IF ( at_i(ji,jj) > epsi10 ) THEN
+                 npti           = npti + 1
+                 nptidx( npti ) = (jj - 1) * jpi + ji
+              ENDIF
+           END DO
+        END DO
+
+        CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i   )
+        CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i   )
+        CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i )
+
+        !-----------------------------------
+        ! Get thickness distribution of ice
+        !-----------------------------------
+
+        Cp = 0.5_wp * grav * (rhow-rhoi) * rhoi / rhow  ! EF: This should go in an initialisation routine but unclear where, won't let me set as a parameter above
+
+        CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
+        ! EF: Note that closing rate is not needed here, could split this function up, and refactor as initial rdgrft_prep is 
+        ! called here, before going into ridging iterations
+
+        ! Initialise
+        strength_pe_sum(:) = 0._wp
+
+        IF ( npti > 0 ) THEN
+
+          DO jl = 1, jpl              ! categories
+            DO ji = 1, npti           ! grid cells
+              IF ( apartf(ji,jl) > 0._wp ) THEN
+                !
+                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 redistribution of ridged ice
+	 	  h2rdg = hrmin(ji,jl) * hrmin(ji,jl) & 
+	 	          + 2._wp * hrmin(ji,jl) * hrexp(ji,jl) & 
+	 	          + 2._wp * hrexp(ji,jl) * hrexp(ji,jl) 
+                END IF
+                
+                dh2rdg = -zhi(ji,jl)*zhi(ji,jl) + h2rdg * hi_hrdg(ji,jl)
+   
+                strength_pe_sum(ji) = strength_pe_sum(ji) + apartf(ji,jl) * dh2rdg
+
+              ELSE
+
+                strength_pe_sum(ji) = strength_pe_sum(ji) + 0._wp
+
+              END IF
+
+            END DO                    ! grid cells
+          END DO                      ! categories
+
+
+          DO ji = 1, npti
+            IF ( zaksum(ji) > epsi10 ) THEN
+              strength_1d(ji) = rn_Cf * Cp * strength_pe_sum(ji) / zaksum(ji)
+            ELSE
+              strength_1d(ji) = 0._wp
+            END IF
+          END DO
+
+
+          ! Enforce a maximum for strength
+          ! 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
+
+
+          ! Convert strength_1d to 2d
+          CALL tab_1d_2d( npti, nptidx(1:npti), strength_1d(1:npti), strength(:,:) )
+
+         END IF  ! number of ice points
+
+
          !                           !--------------------------------------------------!
-      ELSE                           ! Zero strength                                    !
+      ELSE                           ! Zero strength      (this just crashes)           !
          !                           !--------------------------------------------------!
          strength(:,:) = 0._wp
-         ismooth = 0
-      ENDIF
-      !                              !--------------------------------------------------!
+         !
+      ENDIF  ! Ice strength formulation
+      !      
+                                     !--------------------------------------------------!
       SELECT CASE( ismooth )         ! Smoothing ice strength                           !
       !                              !--------------------------------------------------!
+      CASE( 0 )               !--- No smoothing
+
+         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
+
       CASE( 1 )               !--- Spatial smoothing
          DO jj = 2, jpjm1
@@ -918,6 +1205,9 @@
       INTEGER :: ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
-         &                    rn_csrdg  ,                    &
+      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,        &
          &                    ln_partf_exp, rn_astar,        & 
@@ -939,7 +1229,14 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namdyn_rdgrft:'
-         WRITE(numout,*) '      ice strength parameterization Hibler (1979)              ln_str_H79   = ', ln_str_H79 
+         WRITE(numout,*) '      ice strength parameterization Hibler (1979)              ln_str_H79   = ', ln_str_H79
          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_R75   = ', ln_str_R75
+         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
+         WRITE(numout,*) '            e-folding scale of ridged ice for ln_redist_exp    rn_murdg     = ', rn_murdg
          WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_csrdg     = ', rn_csrdg 
          WRITE(numout,*) '      linear ridging participation function                    ln_partf_lin = ', ln_partf_lin
@@ -959,7 +1256,20 @@
       ENDIF
       !
+      IF ( ( ln_str_H79 .AND. ln_str_R75 ) .OR. ( .NOT.ln_str_H79 .AND. .NOT.ln_str_R75 ) ) THEN
+         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one ice strength formulation (ln_str_H79 or ln_str_R75)' )
+      ENDIF
+      !
       IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
          CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
       ENDIF
+      !
+      IF ( ( ln_redist_ufm .AND. ln_redist_exp ) .OR. ( .NOT.ln_redist_ufm .AND. .NOT.ln_redist_exp ) ) THEN
+         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one ice redistribution function (ln_redist_ufm or ln_redist_exp)' )
+      ENDIF
+      !
+      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
