Index: /NEMO/branches/UKMO/NEMO_4.0.1_GM_rossby_radius_cutoff/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_GM_rossby_radius_cutoff/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12860)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_GM_rossby_radius_cutoff/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12861)
@@ -70,4 +70,5 @@
         <field id="bn2"          long_name="squared Brunt-Vaisala frequency"                                           unit="s-1"    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="RoRad"        long_name="Rossby radius"                                                             unit="m" />
 
         <!-- Energy - horizontal divergence -->
Index: /NEMO/branches/UKMO/NEMO_4.0.1_GM_rossby_radius_cutoff/src/OCE/LDF/ldftra.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_GM_rossby_radius_cutoff/src/OCE/LDF/ldftra.F90	(revision 12860)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_GM_rossby_radius_cutoff/src/OCE/LDF/ldftra.F90	(revision 12861)
@@ -84,4 +84,5 @@
    LOGICAL , PUBLIC ::   l_ldftra_time = .FALSE.   !: flag for time variation of the lateral eddy diffusivity coef.
    LOGICAL , PUBLIC ::   l_ldfeiv_time = .FALSE.   !: flag for time variation of the eiv coef.
+   LOGICAL , PUBLIC ::   l_Ro_cutoff   = .FALSE.   !: flag for Rossby-radius cutoff in calcn of eiv coeff
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahtu, ahtv   !: eddy diffusivity coef. at U- and V-points   [m2/s]
@@ -408,8 +409,9 @@
       !!----------------------------------------------------------------------
       !
-      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
-         CALL ldf_eiv( kt, aei0, aeiu, aeiv )
+         CALL ldf_eiv( kt, aei0, aeiu, aeiv, l_Ro_cutoff )
       ENDIF
       !
@@ -424,5 +426,5 @@
             ahtv(:,:,1) = aeiv(:,:,1)
          ELSE                                            ! compute aht. 
-            CALL ldf_eiv( kt, aht0, ahtu, ahtv )
+            CALL ldf_eiv( kt, aht0, ahtu, ahtv, .false. )
          ENDIF
          !
@@ -524,5 +526,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'
@@ -559,4 +561,5 @@
          zah_max = zUfac * (ra*rad)**inn        ! maximum reachable coefficient (value at the Equator)
 
+         l_Ro_cutoff = .FALSE.
          SELECT CASE( nn_aei_ijk_t )         !* Specification of space-time variations
          !
@@ -597,4 +600,13 @@
             l_ldfeiv_time = .TRUE.     ! will be calculated by call to ldf_tra routine in step.F90
             !
+         CASE(  22  )                        !--  time varying 2D field with cutoff --!
+            IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity coef. = F( latitude, longitude, time )'
+            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,*) '           Applied only where Rossby radius resolved by 2.dx'
+            !
+            l_Ro_cutoff   = .TRUE.
+            l_ldfeiv_time = .TRUE.     ! will be calculated by call to ldf_tra routine in step.F90
+            !
          CASE( -30  )                        !-- fixed 3D shape read in file  --!
             IF(lwp) WRITE(numout,*) '   ==>>>   eddy induced velocity coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'
@@ -625,5 +637,5 @@
 
 
-   SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv )
+   SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv, l_Ro_cutoff )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_eiv  ***
@@ -639,8 +651,9 @@
       REAL(wp)                        , INTENT(inout) ::   paei0          ! max value            [m2/s]
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   paeiu, paeiv   ! eiv coefficient      [m2/s]
+      LOGICAL                         , INTENT(in   ) ::   l_Ro_cutoff    ! Rossby radius cutoff yes/no
       !
       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), DIMENSION(jpi,jpj) ::   zn, zah, zhw, zRo, zRo_lim, zaeiw   ! 2D workspace
       !!----------------------------------------------------------------------
       !
@@ -649,4 +662,5 @@
       zah(:,:) = 0._wp
       zRo(:,:) = 0._wp
+      zRo_lim(:,:) = 0._wp
       !                       ! Compute lateral diffusive coefficient at T-point
       IF( ln_traldf_triad ) THEN
@@ -692,11 +706,12 @@
          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)
-         END DO
-      END DO
-
+            zaeiw(ji,jj) = zRo_lim(ji,jj) * zRo_lim(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
+         END DO
+      END DO
+      CALL iom_put('RossRad',zRo)
       !                                         !==  Bound on eiv coeff.  ==!
       z1_f20 = 1._wp / (  2._wp * omega * sin( rad * 20._wp )  )
@@ -704,5 +719,13 @@
          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
+! 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)
+            IF ( l_Ro_cutoff .AND. 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( zzaei, paei0 )
+            ENDIF
          END DO
       END DO
