Skip to content
Snippets Groups Projects
Commit 16cf1ca6 authored by Jérôme Chanut's avatar Jérôme Chanut Committed by Guillaume Samson
Browse files

Resolve "Wrong automatic definition of the laplacian coefficient"

parent 22ee99e9
No related branches found
No related tags found
No related merge requests found
......@@ -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)
......
......@@ -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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment