Index: NEMO/branches/UKMO/NEMO_4.0.4_ice_strength/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_ice_strength/src/ICE/icedyn_rdgrft.F90	(revision 15537)
+++ NEMO/branches/UKMO/NEMO_4.0.4_ice_strength/src/ICE/icedyn_rdgrft.F90	(revision 15542)
@@ -49,4 +49,5 @@
    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
@@ -64,5 +65,6 @@
    LOGICAL  ::   ln_str_R75       ! ice strength parameterization (Rothrock75)
    REAL(wp) ::   rn_Cf            !    ratio of ridging work to PE change in ridging, R75
-   REAL(wp) ::   rn_murdg         !    gives e-folding scale of ridged ice (m^.5), R75
+   INTEGER  ::   rdg_redistr      ! redistribution of ridged ice (0 = uniform (Hibler 1980), 1 = exponential)
+   REAL(wp) ::   rn_murdg         !    gives e-folding scale of ridged ice (m^.5) for rdg_redistr = 1
    REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging            
    LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975))
@@ -94,5 +96,5 @@
       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)   ,  &
+         &      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), STAT=ice_dyn_rdgrft_alloc )
@@ -450,7 +452,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) / zhmean
+               !
+               IF (rdg_redistr == 0) THEN
+                 hrmax  (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
+                 hi_hrdg(ji,jl) = zhi(ji,jl) / zhmean
+               ELSE IF (rdg_redistr == 1) 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
@@ -461,4 +469,5 @@
                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
@@ -516,4 +525,5 @@
       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 
@@ -684,11 +694,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 (rdg_redistr == 0) 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) )
@@ -702,4 +716,32 @@
                      fvol(ji) = 0._wp                  
                   ENDIF
+                  !
+                 ELSE IF (rdg_redistr == 1) THEN ! 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             ! rdg_redistr
 
                   ! Compute the fraction of rafted ice area and volume going to thickness category jl2
@@ -848,9 +890,15 @@
             DO ji = 1, npti           ! grid cells
               IF ( apartf(ji,jl) > 0._wp ) THEN
-
-                ! Uniform distribution of ridged ice
-                h2rdg = (1._wp/3._wp) * (hrmax(ji,jl)**3._wp - hrmin(ji,jl)**3._wp)  &
-                         / (hrmax(ji,jl) - hrmin(ji,jl))
-
+                !
+                IF (rdg_redistr == 0) THEN   ! Uniform distribution 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 (rdg_redistr == 1) THEN   ! Exponential distribution 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)
    
@@ -1039,5 +1087,6 @@
       !!
       NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar,          &
-         &                    ln_str_R75, rn_Cf, rn_murdg,   &
+         &                    ln_str_R75, rn_Cf,             &
+         &                    rdg_redistr, rn_murdg,         &
          &                    rn_crhg, rn_csrdg,             &
          &                    ln_partf_lin, rn_gstar,        &
@@ -1065,5 +1114,7 @@
          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,*) '            gives e-folding scale of ridged ice (m^.5)         rn_murdg     = ', rn_murdg
+         WRITE(numout,*) '      redistribution of ridged ice                             rdg_redistr  = ', rdg_redistr
+         WRITE(numout,*) '            (0 = uniform (Hibler 1980), 1 = exponential)                      '
+         WRITE(numout,*) '            e-folding scale of ridged ice for rdg_redistr=1    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
@@ -1091,4 +1142,8 @@
       ENDIF
       !
+!      IF ( ( rdg_redistr .NE. 0 ) .OR. ( rdg_redistr .NE. 1 ) ) THEN  ! why doesn't this work?
+!         CALL ctl_stop( 'ice_dyn_rdgrft_init: ice redistribution function rdg_redistr must be 0 (uniform) or 1 (exponential)' )
+!      ENDIF
+      !
       IF( .NOT. ln_icethd ) THEN
          rn_porordg = 0._wp
