diff --git a/src/OCE/LDF/ldfdyn.F90 b/src/OCE/LDF/ldfdyn.F90
index 92623ba08b2e2271515a09c0a65ee12dd7ce2111..66c797b44c050d524be43abaf5213624e695790d 100644
--- a/src/OCE/LDF/ldfdyn.F90
+++ b/src/OCE/LDF/ldfdyn.F90
@@ -102,7 +102,7 @@ CONTAINS
!! = 20 = F(i,j) = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
!! =-30 => = F(i,j,k) = shape read in 'eddy_viscosity_3D.nc' file
!! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10)
- !! = 31 = F(i,j,k,t) = F(local velocity) ( |u|e /12 laplacian operator
+ !! = 31 = F(i,j,k,t) = F(local velocity) ( |u|e / 2 laplacian operator
!! or |u|e^3/12 bilaplacian operator )
!! = 32 = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky)
!! ( L^2|D| laplacian operator
@@ -357,7 +357,7 @@ CONTAINS
!! ** Method : time varying eddy viscosity coefficients:
!!
!! nn_ahm_ijk_t = 31 ahmt, ahmf = F(i,j,k,t) = F(local velocity)
- !! ( |u|e /12 or |u|e^3/12 for laplacian or bilaplacian operator )
+ !! ( |u|e / 2 or |u|e^3/2 for laplacian, u|e^3/12 for bilaplacian operator )
!!
!! nn_ahm_ijk_t = 32 ahmt, ahmf = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky)
!! ( L^2|D| or L^4|D|/8 for laplacian or bilaplacian operator )
@@ -379,20 +379,20 @@ CONTAINS
!
CASE( 31 ) !== time varying 3D field ==! = F( local velocity )
!
- IF( ln_dynldf_lap ) THEN ! laplacian operator : |u| e /12 = |u/144| e
- DO jk = 1, jpkm1
+ IF( ln_dynldf_lap ) THEN ! laplacian operator : |u| e / 2 = |u/ 4| e
+ DO jk = 1, jpkm1 ! r1_8 = 1 / (2*2 * 2)
DO_2D( 0, 0, 0, 0 )
zu2pv2_ij = uu(ji ,jj ,jk,Kbb) * uu(ji ,jj ,jk,Kbb) + vv(ji ,jj ,jk,Kbb) * vv(ji ,jj ,jk,Kbb)
zu2pv2_ij_m1 = uu(ji-1,jj ,jk,Kbb) * uu(ji-1,jj ,jk,Kbb) + vv(ji ,jj-1,jk,Kbb) * vv(ji ,jj-1,jk,Kbb)
zu2pv2_ij_p1 = uu(ji ,jj+1,jk,Kbb) * uu(ji ,jj+1,jk,Kbb) + vv(ji+1,jj ,jk,Kbb) * vv(ji+1,jj ,jk,Kbb)
zemax = MAX( e1t(ji,jj) , e2t(ji,jj) )
- ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax * tmask(ji,jj,jk) ! 288= 12*12 * 2
+ ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_8 ) * zemax * tmask(ji,jj,jk)
zemax = MAX( e1f(ji,jj) , e2f(ji,jj) )
- ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax * fmask(ji,jj,jk) ! 288= 12*12 * 2
+ ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_8 ) * zemax * fmask(ji,jj,jk)
END_2D
END DO
ELSEIF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( |u| e^3 /12 ) = sqrt( |u/144| e ) * e
- DO jk = 1, jpkm1
+ DO jk = 1, jpkm1 ! r1_288 = 1 / (12*12 * 2)
DO_2D( 0, 0, 0, 0 )
zu2pv2_ij = uu(ji ,jj ,jk,Kbb) * uu(ji ,jj ,jk,Kbb) + vv(ji ,jj ,jk,Kbb) * vv(ji ,jj ,jk,Kbb)
zu2pv2_ij_m1 = uu(ji-1,jj ,jk,Kbb) * uu(ji-1,jj ,jk,Kbb) + vv(ji ,jj-1,jk,Kbb) * vv(ji ,jj-1,jk,Kbb)
diff --git a/src/OCE/LDF/ldftra.F90 b/src/OCE/LDF/ldftra.F90
index cfa3f05aa0fd4cc4c0a7e42958e4e3879120c0de..510166d0955b9fd418b04aae9ba2965af7799c4b 100644
--- a/src/OCE/LDF/ldftra.F90
+++ b/src/OCE/LDF/ldftra.F90
@@ -389,7 +389,7 @@ CONTAINS
!! with a reduction to 0 in vicinity of the Equator
!! nn_aht_ijk_t = 21 ahtu, ahtv = F(i,j, t) = F(growth rate of baroclinic instability)
!!
- !! = 31 ahtu, ahtv = F(i,j,k,t) = F(local velocity) ( |u|e /12 laplacian operator
+ !! = 31 ahtu, ahtv = F(i,j,k,t) = F(local velocity) ( |u|e / 2 laplacian operator
!! or |u|e^3/12 bilaplacian operator )
!!
!! * time varying EIV coefficients: call to ldf_eiv routine
@@ -440,10 +440,10 @@ CONTAINS
END DO
!
CASE( 31 ) !== time varying 3D field ==! = F( local velocity )
- IF( ln_traldf_lap ) THEN ! laplacian operator |u| e /12
+ IF( ln_traldf_lap ) THEN ! laplacian operator |u| e / 2
DO jk = 1, jpkm1
- ahtu(:,:,jk) = ABS( uu(:,:,jk,Kbb) ) * e1u(:,:) * r1_12 ! n.b. uu,vv are masked
- ahtv(:,:,jk) = ABS( vv(:,:,jk,Kbb) ) * e2v(:,:) * r1_12
+ ahtu(:,:,jk) = ABS( uu(:,:,jk,Kbb) ) * e1u(:,:) * r1_2 ! n.b. uu,vv are masked
+ ahtv(:,:,jk) = ABS( vv(:,:,jk,Kbb) ) * e2v(:,:) * r1_2
END DO
ELSEIF( ln_traldf_blp ) THEN ! bilaplacian operator sqrt( |u| e^3 /12 ) = sqrt( |u| e /12 ) * e
DO jk = 1, jpkm1