Changeset 901 for trunk/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90

Timestamp:

2008-04-22T21:16:18+02:00 (16 years ago)

Author:

rblod

Message:

Change ldf computation to take into account the grid anisotropy, see ticket #115

File:

• trunk/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 (modified) (3 diffs)

trunk/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90

@@ -36 +36 @@
 
       !! * Local variables
-      REAL(wp) ::   za00, zdx_max
+      INTEGER :: ji, jj
+      REAL(wp) ::   za00, zd_max, zetmax, zeumax, zefmax, zevmax
       !!----------------------------------------------------------------------
 
@@ -50 +51 @@
          ! (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 +96 @@
          !       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