Changeset 1611 for LMDZ5/trunk/libf/dyn3d

Timestamp:

Jan 25, 2012, 3:31:54 PM (12 years ago)

Author:

lguez

Message:

Avoid aliasing in "inidissip".
Correct misplaced continuation lines in "initdynav".

File:

• LMDZ5/trunk/libf/dyn3d/inidissip.F90 (modified) (3 diffs)

LMDZ5/trunk/libf/dyn3d/inidissip.F90

@@ -28 +28 @@
 ! Local variables:
   REAL fact,zvert(llm),zz
-  REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
+  REAL zh(ip1jmp1),zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1)
+  real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1,llm)
   REAL ullm,vllm,umin,vmin,zhmin,zhmax
-  REAL zllm,z1llm
+  REAL zllm
 
   INTEGER l,ij,idum,ii
@@ -78 +79 @@
   DO l = 1,50
      IF(lstardis) THEN
-        CALL divgrad2(1,zh,deltap,niterh,zh)
+        CALL divgrad2(1,zh,deltap,niterh,divgra)
      ELSE
-        CALL divgrad (1,zh,niterh,zh)
+        CALL divgrad (1,zh,niterh,divgra)
      ENDIF
 
-     CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
-
-     zllm  = ABS( zhmax )
-     z1llm = 1./zllm
-     DO ij = 1,ip1jmp1
-        zh(ij) = zh(ij)* z1llm
-     ENDDO
+     zllm  = ABS(maxval(divgra))
+     zh = divgra / zllm
   ENDDO
 
@@ -123 +119 @@
            !cccc             CALL covcont( 1,zu,zv,zu,zv )
            IF(lstardis) THEN
-              CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
+              CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy )
            ELSE
-              CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
+              CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy )
            ENDIF
         ELSE
            IF(lstardis) THEN
-              CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
+              CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy )
            ELSE
-              CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
+              CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy )
            ENDIF
         ENDIF
 
-        CALL minmax(iip1*jjp1,zu,umin,ullm )
-        CALL minmax(iip1*jjm, zv,vmin,vllm )
-
-        ullm = ABS  ( ullm )
-        vllm = ABS  ( vllm )
-
-        zllm  = MAX( ullm,vllm )
-        z1llm = 1./ zllm
-        DO ij = 1, ip1jmp1
-           zu(ij) = zu(ij)* z1llm
-        ENDDO
-        DO ij = 1, ip1jm
-           zv(ij) = zv(ij)* z1llm
-        ENDDO
+        zllm = max(abs(maxval(gx)), abs(maxval(gy)))
+        zu = gx / zllm
+        zv = gy / zllm
      end DO