Merge branch '71-wrong-automatic-definition-of-the-laplacian-coefficient-2' into 'main'

Resolve "Wrong automatic definition of the laplacian coefficient"

Closes #71

See merge request nemo/nemo!298

src/OCE/LDF/ldfdyn.F90

@@ -103,7 +103,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
@@ -351,7 +351,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 )
@@ -373,20 +373,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)

src/OCE/LDF/ldftra.F90

@@ -384,7 +384,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
@@ -435,10 +435,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