Index: /NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/cfgs/SHARED/field_def_nemo-oce.xml	(revision 14341)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/cfgs/SHARED/field_def_nemo-oce.xml	(revision 14342)
@@ -129,4 +129,10 @@
         <field id="beta"         long_name="haline contraction"                                                        unit="1e3"    grid_ref="grid_T_3D" />
         <field id="rhop"         long_name="potential density (sigma0)"        standard_name="sea_water_sigma_theta"   unit="kg/m3"  grid_ref="grid_T_3D" />
+        <field id="N_2d"         long_name="Depth-mean N"                                                              unit="m/s" />
+        <field id="modslp"       long_name="sqrt( slpi^2 + slpj^2 )"                                                   unit="1"      grid_ref="grid_T_3D" />
+        <field id="RossRad"      long_name="Rossby radius"                                                             unit="m" />
+        <field id="RossRadlim"   long_name="Rossby radius (limited)"                                                   unit="m" />
+        <field id="Tclinic_recip" long_name="recip of baroclinic timescale"                                            unit="s-1" />
+        <field id="RR_GS"        long_name="Rossby radius / min(dx,dy)"                                                unit="1" />
 
         <!-- Energy - horizontal divergence -->
Index: /NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/cfgs/SHARED/namelist_ref	(revision 14341)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/cfgs/SHARED/namelist_ref	(revision 14342)
@@ -838,8 +838,9 @@
       !                             !   = 30 F(i,j,k)  =ldf_c2d * ldf_c1d
       !                        !  time invariant coefficients:  aei0 = 1/2  Ue*Le 
-      rn_Ue        = 0.02           !  lateral diffusive velocity [m/s] (nn_aht_ijk_t= 0, 10, 20, 30)
-      rn_Le        = 200.e+3        !  lateral diffusive length   [m]   (nn_aht_ijk_t= 0, 10)
+      rn_Ue        = 0.02           !  lateral diffusive velocity [m/s] (nn_aei_ijk_t= 0, 10, 20, 30)
+      rn_Le        = 200.e+3        !  lateral diffusive length   [m]   (nn_aei_ijk_t= 0, 10)
       !
       ln_ldfeiv_dia =.false.   ! diagnose eiv stream function and velocities
+      nn_ldfeiv_shape = 0           !  shape of bounding coefficient    (nn_aei_ijk_t= 21 only)
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/OCE/DIA/diawri.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/OCE/DIA/diawri.F90	(revision 14341)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/OCE/DIA/diawri.F90	(revision 14342)
@@ -388,4 +388,12 @@
       ENDIF
 
+      IF( (.NOT.l_ldfeiv_time) .AND. ( iom_use('RossRad')  .OR. iom_use('RossRadlim') &
+            &                     .OR. iom_use('TclinicR') .OR. iom_use('RR_GG')      &
+            &                     .OR. iom_use('aeiu_2d')  .OR. iom_use('aeiv_2d') ) ) THEN
+         CALL ldf_eiv(kt, 75.0, z2d, z3d(:,:,1))
+         CALL iom_put('aeiu_2d', z2d)
+         CALL iom_put('aeiv_2d', z3d(:,:,1))
+      ENDIF
+
       IF( iom_use("tosmint_"//ttype) ) THEN
          z2d(:,:) = 0._wp
Index: /NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/OCE/LDF/ldftra.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/OCE/LDF/ldftra.F90	(revision 14341)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/OCE/LDF/ldftra.F90	(revision 14342)
@@ -72,4 +72,5 @@
    REAL(wp), PUBLIC ::      rn_Ue               !: lateral diffusive velocity  [m/s]
    REAL(wp), PUBLIC ::      rn_Le               !: lateral diffusive length    [m]
+   INTEGER,  PUBLIC ::   nn_ldfeiv_shape     !: shape of bounding coefficient (Treguier et al formulation only)
    
    !                                  ! Flag to control the type of lateral diffusive operator
@@ -409,5 +410,6 @@
       !!----------------------------------------------------------------------
       !
-      IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN       ! eddy induced velocity coefficients
+      IF( ln_ldfeiv .AND. ( nn_aei_ijk_t == 21 .OR. nn_aei_ijk_t == 22 ) ) THEN       
+         !                                ! eddy induced velocity coefficients
          !                                ! =F(growth rate of baroclinic instability)
          !                                ! max value aeiv_0 ; decreased to 0 within 20N-20S
@@ -502,5 +504,6 @@
       !!
       NAMELIST/namtra_eiv/ ln_ldfeiv   , ln_ldfeiv_dia,   &   ! eddy induced velocity (eiv)
-         &                 nn_aei_ijk_t, rn_Ue, rn_Le         ! eiv  coefficient
+         &                 nn_aei_ijk_t, rn_Ue, rn_Le,    &   ! eiv  coefficient
+         &                 nn_ldfeiv_shape
       !!----------------------------------------------------------------------
       !
@@ -525,5 +528,5 @@
          WRITE(numout,*) '      eiv streamfunction & velocity diag.        ln_ldfeiv_dia = ', ln_ldfeiv_dia
          WRITE(numout,*) '      coefficients :'
-         WRITE(numout,*) '         type of time-space variation            nn_aei_ijk_t  = ', nn_aht_ijk_t
+         WRITE(numout,*) '         type of time-space variation            nn_aei_ijk_t  = ', nn_aei_ijk_t
          WRITE(numout,*) '         lateral diffusive velocity (if cst)     rn_Ue         = ', rn_Ue, ' m/s'
          WRITE(numout,*) '         lateral diffusive length   (if cst)     rn_Le         = ', rn_Le, ' m'
@@ -595,4 +598,5 @@
             IF(lwp) WRITE(numout,*) '                                       = F( growth rate of baroclinic instability )'
             IF(lwp) WRITE(numout,*) '           maximum allowed value: aei0 = ', aei0, ' m2/s'
+            IF(lwp) WRITE(numout,*) '           shape of bounding coefficient : ',nn_ldfeiv_shape
             !
             l_ldfeiv_time = .TRUE.     ! will be calculated by call to ldf_tra routine in step.F90
@@ -638,16 +642,22 @@
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
-      REAL(wp)                        , INTENT(inout) ::   paei0          ! max value            [m2/s]
+      REAL(wp)                        , INTENT(in   ) ::   paei0          ! max value            [m2/s]
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   paeiu, paeiv   ! eiv coefficient      [m2/s]
       !
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
-      REAL(wp) ::   zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei    ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zn, zah, zhw, zRo, zaeiw   ! 2D workspace
+      REAL(wp) ::   zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei, z2_3    ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj) ::   zn, zah, zhw, zRo, zRo_lim, zTclinic_recip, zaeiw, zratio   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmodslp ! 3D workspace
       !!----------------------------------------------------------------------
       !
       zn (:,:) = 0._wp        ! Local initialization
+      zmodslp(:,:,:) = 0._wp
       zhw(:,:) = 5._wp
       zah(:,:) = 0._wp
       zRo(:,:) = 0._wp
+      zRo_lim(:,:) = 0._wp
+      zTclinic_recip(:,:) = 0._wp
+      zratio(:,:) = 0._wp 
+      zaeiw(:,:) = 0._wp    
       !                       ! Compute lateral diffusive coefficient at T-point
       IF( ln_traldf_triad ) THEN
@@ -682,6 +692,7 @@
                   ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
                   ze3w = e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
-                  zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
-                     &                            + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w
+                  zmodslp(ji,jj,jk) =  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
+                     &               + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)
+                  zah(ji,jj) = zah(ji,jj) + zn2 * zmodslp(ji,jj,jk) * ze3w
                   zhw(ji,jj) = zhw(ji,jj) + ze3w
                END DO
@@ -693,19 +704,88 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 )
-            ! Rossby radius at w-point taken betwenn 2 km and  40km
-            zRo(ji,jj) = MAX(  2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 )  )
+            ! Rossby radius at w-point taken between 2 km and  40km
+            zRo(ji,jj) = .4 * zn(ji,jj) / zfw
+            zRo_lim(ji,jj) = MAX(  2.e3 , MIN( zRo(ji,jj), 40.e3 )  )
             ! Compute aeiw by multiplying Ro^2 and T^-1
-            zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
+            zTclinic_recip(ji,jj) = SQRT( MAX(zah(ji,jj),0._wp) / zhw(ji,jj) ) * tmask(ji,jj,1)
+            zaeiw(ji,jj) = zRo_lim(ji,jj) * zRo_lim(ji,jj) * zTclinic_recip(ji,jj) 
          END DO
       END DO
-
+      IF( iom_use('N_2d') ) CALL iom_put('N_2d',zn(:,:)/zhw(:,:))
+      IF( iom_use('modslp') ) CALL iom_put('modslp',SQRT(zmodslp(:,:,:)) )
+      CALL iom_put('RossRad',zRo)
+      CALL iom_put('RossRadlim',zRo_lim)
+      CALL iom_put('Tclinic_recip',zTclinic_recip)
       !                                         !==  Bound on eiv coeff.  ==!
       z1_f20 = 1._wp / (  2._wp * omega * sin( rad * 20._wp )  )
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj)     ! tropical decrease
-            zaeiw(ji,jj) = MIN( zzaei , paei0 )                                  ! Max value = paei0
-         END DO
-      END DO
+      z2_3 = 2._wp/3._wp
+
+      SELECT CASE(nn_ldfeiv_shape)
+         CASE(0) !! Standard shape applied - decrease in tropics and cap. 
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj)     ! tropical decrease
+                  zaeiw(ji,jj) = MIN( zzaei, paei0 )
+               END DO
+            END DO
+
+         CASE(1) !! Abrupt cut-off on Rossby radius:
+! JD : modifications here to introduce scaling by local rossby radius of deformation vs local grid scale
+! arbitrary decision that GM is de-activated if local rossy radius larger than 2 times local grid scale
+! based on Hallberg (2013)
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF ( zRo(ji,jj) >= ( 2._wp * MIN( e1t(ji,jj), e2t(ji,jj) ) ) ) THEN
+! TODO : use a version of zRo that integrates over a few time steps ?
+                      zaeiw(ji,jj) = 0._wp
+                  ELSE
+                      zaeiw(ji,jj) = MIN( zaeiw(ji,jj), paei0 )
+                  ENDIF
+               END DO
+            END DO
+
+         CASE(2) !! Rossby radius ramp type 1:
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zratio(ji,jj) = zRo(ji,jj)/MIN(e1t(ji,jj),e2t(ji,jj))
+                  zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, z2_3*(2._wp - zratio(ji,jj)) ) ) * paei0 )
+               END DO
+            END DO
+            CALL iom_put('RR_GS',zratio)
+
+         CASE(3) !! Rossby radius ramp type 2:
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zratio(ji,jj) = MIN(e1t(ji,jj),e2t(ji,jj))/zRo(ji,jj)
+                  zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, z2_3*( zratio(ji,jj) - 0.5_wp ) ) ) * paei0 )
+               END DO
+            END DO
+
+         CASE(4) !! bathymetry ramp:
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, 0.001*(ht_0(ji,jj) - 2000._wp) ) ) * paei0 )
+               END DO
+            END DO
+
+         CASE(5) !! Rossby radius ramp type 1 applied to Treguier et al coefficient rather than cap:
+                 !! Note the ramp is RR/GS=[2.0,1.0] (not [2.0,0.5] as for cases 2,3) and we ramp up 
+                 !! to 5% of the Treguier et al coefficient, aiming for peak values of around 100m2/s
+                 !! at high latitudes rather than 2000m2/s which is what you get in eORCA025 with an 
+                 !! uncapped coefficient.
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zratio(ji,jj) = zRo(ji,jj)/MIN(e1t(ji,jj),e2t(ji,jj))
+                  zaeiw(ji,jj) = MAX( 0._wp, MIN( 1._wp, 2._wp - zratio(ji,jj) ) ) * 0.05 * zaeiw(ji,jj)
+                  zaeiw(ji,jj) = MIN( zaeiw(ji,jj), paei0 )
+               END DO
+            END DO
+            CALL iom_put('RR_GS',zratio)
+
+         CASE DEFAULT
+               CALL ctl_stop('ldf_eiv: Unrecognised option for nn_ldfeiv_shape.')         
+
+      END SELECT
+
       CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1. )       ! lateral boundary condition
       !               
@@ -722,5 +802,5 @@
          paeiv(:,:,jk) = paeiv(:,:,1) * vmask(:,:,jk)
       END DO
-      !  
+      !
    END SUBROUTINE ldf_eiv
 
