Changeset 162 for trunk/NEMO/OPA_SRC/LDF

Timestamp:

2004-11-03T15:47:11+01:00 (20 years ago)

Author:

opalod

Message:

CT + CL : BUGFIX108 : Add latitudinal variation of viscosity for the ORCA 4� configuration to damp instabilities

File:

• trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90 (modified) (6 diffs)

trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90

@@ -190 +190 @@
       !! * local variables
       INTEGER ::   ji, jj, jk, jn      ! dummy loop indices
+      INTEGER ::   ii0, ii1, ij0, ij1  ! temporary integers
       INTEGER ::   inum = 11           ! temporary logical unit
       INTEGER ::   iost, iim, ijm
@@ -196 +197 @@
       INTEGER, DIMENSION(jpi   , jpj   ) ::   icof
 
-      REAL(wp) ::   zahmeq, zcoff, zcoft, zmsk
+      REAL(wp) ::   &
+         zahmeq, zcoff, zcoft, zmsk,   & ! ???
+         zemax, zemin, zeref, zahmm
+      REAL(wp), DIMENSION(jpi,jpj) ::   zahm0
       REAL(wp), DIMENSION(jpk) ::   zcoef
 
@@ -288 +292 @@
       ! symmetric in the south hemisphere)
       
-      IF( jp_cfg == 4 )   zahmeq = 5.0 * aht0
-      IF( jp_cfg == 2 )   zahmeq =       aht0
+      IF( jp_cfg == 4 )   THEN
+         zahmeq = 5.0 * aht0
+         zahmm  = min( 160000.0, ahm0)
+         zemax = MAXVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
+         zemin = MINVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
+         zeref = MAXVAL ( e1t(:,:) * e2t(:,:),   &
+             &   tmask(:,:,1) .GE. 0.5 .AND. ABS(gphit(:,:)) .GT. 50. )
+ 
+         DO jj = 1, jpj
+           DO ji = 1, jpi
+              zmsk = e1t(ji,jj) * e2t(ji,jj)
+              IF( abs(gphit(ji,jj)) .LE. 15 ) THEN 
+                 zahm0(ji,jj) = ahm0
+              ELSE 
+                 IF( zmsk .GE. zeref ) THEN 
+                    zahm0(ji,jj) = ahm0
+                 ELSE 
+                    zahm0(ji,jj) = zahmm + (ahm0-zahmm)*(1.0 -   &
+                        &          cos((rpi*0.5*(zmsk-zemin)/(zeref-zemin))))
+                 ENDIF 
+              ENDIF 
+           END DO 
+         END DO 
+      ENDIF
+
+      IF( jp_cfg == 2 )   THEN
+         zahmeq     = aht0
+         zahmm      = ahm0
+         zahm0(:,:) = ahm0
+      ENDIF
 
       DO jj = 1, jpj
          DO ji = 1, jpi
             IF( ABS(gphif(ji,jj)) >= 20.) THEN
-               ahm2(ji,jj,1) =  ahm0
+               ahm2(ji,jj,1) =  zahm0(ji,jj)
             ELSEIF( ABS(gphif(ji,jj)) <= 2.5) THEN
                ahm2(ji,jj,1) =  zahmeq
             ELSE
-               ahm2(ji,jj,1) = zahmeq    &
-                             + (ahm0-zahmeq)/2.*(1.-COS( rad*(ABS(gphif(ji,jj))-2.5)*180./17.5 ) )
+               ahm2(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2.   &
+                  &            *(1.-COS( rad*(ABS(gphif(ji,jj))-2.5)*180./17.5 ) )
             ENDIF
             IF( ABS(gphit(ji,jj)) >= 20.) THEN
-               ahm1(ji,jj,1) =  ahm0
+               ahm1(ji,jj,1) =  zahm0(ji,jj)
             ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
                ahm1(ji,jj,1) =  zahmeq
             ELSE
-               ahm1(ji,jj,1) = zahmeq    &
-                             + (ahm0-zahmeq)/2.*(1.-COS( rad*(ABS(gphit(ji,jj))-2.5)*180./17.5 ) )
+               ahm1(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2.   &
+                  &            *(1.-COS( rad*(ABS(gphit(ji,jj))-2.5)*180./17.5 ) )
             ENDIF
          END DO
@@ -317 +349 @@
          DO ji = 1, jpim1
             zcoft = float( icof(ji,jj) ) / 100.
-            ahm1(ji,jj,1) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj,1)
+            ahm1(ji,jj,1) = zcoft * zahm0(ji,jj) + (1.-zcoft) * ahm1(ji,jj,1)
          END DO
       END DO
@@ -331 +363 @@
                      / (zmsk * 100.)
             ENDIF
-            ahm2(ji,jj,1) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj,1)
+            ahm2(ji,jj,1) = zcoff * zahm0(ji,jj) + (1.-zcoff) * ahm2(ji,jj,1)
          END DO
       END DO
@@ -350 +382 @@
       
       DO jk = 2, jpk
-         ahm1(:,:,jk) = MIN( ahm0, zcoef(jk) * ahm1(:,:,1) )
-         ahm2(:,:,jk) = MIN( ahm0, zcoef(jk) * ahm2(:,:,1) )
-      END DO
+         ahm1(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm1(:,:,1) )
+         ahm2(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm2(:,:,1) )
+      END DO
+
+      IF( jp_cfg == 4 )   THEN
+         ! Limit AHM in Gibraltar strait
+         ij0 = 50   ;   ij1 = 53
+         ii0 = 69   ;   ii1 = 71
+         DO jk = 1, jpk
+            ahm1(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) = min( zahmm, ahm1 (mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) ) 
+            ahm2(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) = min( zahmm, ahm2 (mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) )
+         END DO
+      ENDIF
       
       ! Lateral boundary conditions on ( ahm1, ahm2 )