Index: /trunk/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90
===================================================================
--- /trunk/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90	(revision 900)
+++ /trunk/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90	(revision 901)
@@ -36,5 +36,6 @@
 
       !! * Local variables
-      REAL(wp) ::   za00, zdx_max
+      INTEGER :: ji, jj
+      REAL(wp) ::   za00, zd_max, zetmax, zeumax, zefmax, zevmax
       !!----------------------------------------------------------------------
 
@@ -50,14 +51,19 @@
          ! (USER: modify ahm1 and ahm2 following your desiderata)
 
-         zdx_max = MAXVAL( e1t(:,:) )
-         IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain
+         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
+         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
 
          IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1'
-         IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...'
-         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0
-
-         za00 = ahm0 / zdx_max
-         ahm1(:,:) = za00 * e1t(:,:)
-         ahm2(:,:) = za00 * e1f(:,:)
+         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
+
+         za00 = ahm0 / zd_max
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
+               zefmax = MAX( e1f(ji,jj), e2f(ji,jj) )
+               ahm1(ji,jj) = za00 * zetmax
+               ahm2(ji,jj) = za00 * zefmax
+            END DO
+         END DO
 
          IF( ln_dynldf_iso ) THEN
@@ -90,14 +96,19 @@
          !       in the to horizontal direction
 
-         zdx_max = MAXVAL( e1u(:,:) )
-         IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain
+         zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) )
+         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
 
          IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 '
-         IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...'
-         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0
-
-         za00 = ahm0 / ( zdx_max * zdx_max * zdx_max )
-         ahm3(:,:) = za00 * e1u(:,:) * e1u(:,:) * e1u(:,:)
-         ahm4(:,:) = za00 * e1v(:,:) * e1v(:,:) * e1v(:,:)
+         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
+
+         za00 = ahm0 / ( zd_max * zd_max * zd_max )
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )
+               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )
+               ahm3(ji,jj) = za00 * zeumax * zeumax * zeumax
+               ahm4(ji,jj) = za00 * zevmax * zevmax * zevmax
+            END DO
+         END DO
 
          ! Control print
Index: /trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90
===================================================================
--- /trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90	(revision 900)
+++ /trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90	(revision 901)
@@ -38,6 +38,9 @@
          !               ! ( 0 < zr < 1 )
          zh = 500.,   &  ! depth of at which start the reduction ( > dept(1) )
-         zdx_max  ,   &  ! maximum grid spacing over the global domain
+         zd_max   ,   &  ! maximum grid spacing over the global domain
          za00, zc, zd    ! temporary scalars
+      REAL(wp) ::        &
+         zetmax, zefmax, &
+         zeumax, zevmax   
       REAL(wp), DIMENSION(jpk) ::   zcoef   ! temporary workspace
       !!----------------------------------------------------------------------
@@ -54,13 +57,11 @@
          ! (USER: modify ahm1 and ahm2 following your desiderata)
 
-         zdx_max = MAXVAL( e1t(:,:) )
-         IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain
+         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
+         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
 
          IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1'
-         IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...'
-         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0
-
-
-         za00 = ahm0 / zdx_max
+         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
+
+         za00 = ahm0 / zd_max
 
          IF( ln_dynldf_iso ) THEN
@@ -74,6 +75,12 @@
          CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm2 )   ! vertical profile
          DO jk = 1,jpk
-            ahm1(:,:,jk) = za00 * e1t(:,:) * ahm1(:,:,jk)
-            ahm2(:,:,jk) = za00 * e1f(:,:) * ahm2(:,:,jk)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
+                  zefmax = MAX( e1f(ji,jj), e2f(ji,jj) )
+                  ahm1(ji,jj,jk) = za00 * zetmax * ahm1(ji,jj,jk)
+                  ahm2(ji,jj,jk) = za00 * zefmax * ahm2(ji,jj,jk)
+               END DO
+            END DO
          END DO
 
@@ -109,14 +116,19 @@
       IF( ln_dynldf_bilap ) THEN
 
-         zdx_max = MAXVAL( e1u(:,:) )
-         IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain
+         zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) )
+         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
 
          IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 '
-         IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...'
-         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0
-
-         za00 = ahm0 / ( zdx_max * zdx_max * zdx_max )
-         ahm3(:,:,1) = za00 * e1u(:,:) * e1u(:,:) * e1u(:,:)
-         ahm4(:,:,1) = za00 * e1v(:,:) * e1v(:,:) * e1v(:,:)
+         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
+
+         za00 = ahm0 / ( zd_max * zd_max * zd_max )
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )
+               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )
+               ahm3(ji,jj,1) = za00 * zeumax * zeumax * zeumax
+               ahm4(ji,jj,1) = za00 * zevmax * zevmax * zevmax
+            END DO
+         END DO
 
          zh = MAX( zh, fsdept(1,1,1) )   ! at least the first reach ahm0
Index: /trunk/NEMO/OPA_SRC/LDF/ldftra_c2d.h90
===================================================================
--- /trunk/NEMO/OPA_SRC/LDF/ldftra_c2d.h90	(revision 900)
+++ /trunk/NEMO/OPA_SRC/LDF/ldftra_c2d.h90	(revision 901)
@@ -31,11 +31,9 @@
 
       !! * Local variables
-# if defined key_traldf_eiv && defined key_orca_r4
       INTEGER ::   ji, jj                  ! dummy loop indices
-# endif
 # if defined key_orca_r4
       INTEGER :: i1, i2, j1, j2
 # endif
-      REAL(wp) ::   za00, zdx_max
+      REAL(wp) ::   za00, zd_max, zeumax, zevmax, zetmax
       
       !!----------------------------------------------------------------------
@@ -53,6 +51,6 @@
       ENDIF
 
-      zdx_max = MAXVAL( e1t(:,:) )
-      IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain
+      zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
+      IF( lk_mpp ) CALL mpp_max( zd_max )   ! max over the global domain
 
       ! harmonic operator : (U-, V-, W-points)
@@ -60,18 +58,31 @@
       IF( ln_traldf_lap ) THEN
   
-         za00 = aht0 / zdx_max
+         za00 = aht0 / zd_max
   
-!RB         ahtu(:,:) = za00 * e1u(:,:)               ! set ahtu = ahtv at u- and v-points,
-!RB         ahtv(:,:) = za00 * e1f(:,:)               ! and ahtw at w-point (idem T-point)
-!RB         ahtw(:,:) = za00 * e1t(:,:)               ! 
-         ahtu(:,:) = aht0                ! set ahtu = ahtv at u- and v-points,
-         ahtv(:,:) = aht0                ! and ahtw at w-point (idem T-point)
-         ahtw(:,:) = aht0                !
-
-
+         DO jj = 1, jpj 
+            DO ji = 1, jpi 
+               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) ) 
+               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) ) 
+               zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
+               ahtu(ji,jj) = za00 * zeumax ! set ahtu = ahtv at u- and v-points, 
+               ahtv(ji,jj) = za00 * zevmax ! and ahtw at w-point (idem T-point) 
+               ahtw(ji,jj) = za00 * zetmax ! 
+            END DO
+         END DO
 
          CALL lbc_lnk( ahtu, 'U', 1. )   ! Lateral boundary conditions
          CALL lbc_lnk( ahtv, 'V', 1. )   ! (no change of sign)
          CALL lbc_lnk( ahtw, 'W', 1. )
+
+         ! Special case for ORCA R2 and R4 configurations (overwrite the value of ahtu ahtv ahtw)
+         ! ==============================================
+         IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) )   THEN
+            ahtu(:,:) = aht0              ! set ahtu = ahtv at u- and v-points,
+            ahtv(:,:) = aht0              ! and ahtw at w-point
+            ahtw(:,:) = aht0              ! (here : no space variation)
+            IF(lwp) WRITE(numout,*) '               ORCA R2 or R4 case'
+            IF(lwp) WRITE(numout,*) '               Constant values used for eddy diffusivity coefficients'
+            IF(lwp) WRITE(numout,*) '               Variation lat/lon only for eddy induced velocity coefficients'
+         ENDIF
 
          ! Control print
@@ -99,6 +110,14 @@
          !       in the to horizontal direction
 
-         za00 = aht0 / ( zdx_max * zdx_max * zdx_max )
-         ahtt(:,:) = za00 * e1t(:,:) * e1t(:,:) *e1t(:,:)      ! set ahtt at T-point 
+         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
+         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
+
+         za00 = aht0 / ( zd_max * zd_max * zd_max )
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
+               ahtt(ji,jj) = za00 * zetmax * zetmax * zetmax      ! set ahtt at T-point
+            END DO
+         END DO
 
          CALL lbc_lnk( ahtt, 'T', 1. )   ! Lateral boundary conditions on ( ahtt )
Index: /trunk/NEMO/OPA_SRC/LDF/ldftra_c3d.h90
===================================================================
--- /trunk/NEMO/OPA_SRC/LDF/ldftra_c3d.h90	(revision 900)
+++ /trunk/NEMO/OPA_SRC/LDF/ldftra_c3d.h90	(revision 901)
@@ -45,4 +45,5 @@
          IF(lwp) WRITE(numout,*) ' ldf_tra_c3d : 3D eddy diffusivity and eddy induced velocity coefficients'
          IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   --  '
+         IF(lwp) WRITE(numout,*) '               Coefficients set to constant'
          IF(lwp) WRITE(numout,*)
       ELSE
@@ -50,4 +51,5 @@
          IF(lwp) WRITE(numout,*) ' ldf_tra_c3d : 3D eddy diffusivity coefficient'
          IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   -- '
+         IF(lwp) WRITE(numout,*) '               Coefficients set to constant'
          IF(lwp) WRITE(numout,*)
       ENDIF
