Changeset 1665 for LMDZ5/branches/testing/libf/filtrez/filtreg_mod.F90

Timestamp:

Oct 9, 2012, 3:35:26 PM (12 years ago)

Author:

Laurent Fairhead

Message:

Version testing basée sur la r1628

http://lmdz.lmd.jussieu.fr/utilisateurs/distribution-du-modele/versions-intermediaires

---

Testing release based on r1628

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/filtrez/filtreg_mod.F90 (modified) (20 diffs)

LMDZ5/branches/testing/libf/filtrez/filtreg_mod.F90

@@ +1 @@
+!
+! $Id $
+!
 MODULE filtreg_mod
 
@@ -42 +45 @@
     INTEGER ixmineq
 #endif
-    EXTERNAL  inifgn
     !
     ! ------------------------------------------------------------
@@ -71 +73 @@
     CALL inifgn(eignvl)
     !
-    PRINT *,' EIGNVL '
+    PRINT *,'inifilr: EIGNVL '
     PRINT 250,eignvl
-250 FORMAT( 1x,5e13.6)
+250 FORMAT( 1x,5e14.6)
     !
     ! compute eigenvalues and eigenfunctions
@@ -113 +115 @@
 #endif
     !
+    ! For a regular grid, we want the filter to start at latitudes
+    ! corresponding to lengths dx of the same size as dy (in terms
+    ! of angles: dx=2*dy) => at colat0=0.5 (i.e. colatitude=30 degrees
+    !  <=> latitude=60 degrees).
+    ! Same idea for the zoomed grid: start filtering polewards as soon
+    ! as length dx becomes of the same size as dy 
     !
     colat0  =  MIN( 0.5, dymin/dxmin )
@@ -158 +166 @@
     imx  = iim
     !
-    PRINT *,' TRUNCATION AT ',imx
-    !
+    PRINT *,'inifilr: TRUNCATION AT ',imx
+    !
+! Ehouarn: set up some defaults
+    jfiltnu=2 ! avoid north pole
+    jfiltsu=jjm ! avoid south pole (which is at jjm+1)
+    jfiltnv=1 ! NB: no poles on the V grid
+    jfiltsv=jjm
+
     DO j = 2, jjm/2+1
        cof = COS( rlatu(j) )/ colat0
        IF ( cof .LT. 1. ) THEN
-          IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) jfiltnu= j
+          IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) THEN
+            jfiltnu= j
+          ENDIF
        ENDIF
 
        cof = COS( rlatu(jjp1-j+1) )/ colat0
        IF ( cof .LT. 1. ) THEN
-          IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) &
+          IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) THEN
                jfiltsu= jjp1-j+1
+          ENDIF
        ENDIF
     ENDDO
@@ -176 +193 @@
        cof = COS( rlatv(j) )/ colat0
        IF ( cof .LT. 1. ) THEN
-          IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) jfiltnv= j
+          IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) THEN
+            jfiltnv= j
+          ENDIF
        ENDIF
 
        cof = COS( rlatv(jjm-j+1) )/ colat0
        IF ( cof .LT. 1. ) THEN
-          IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) &
+          IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) THEN
                jfiltsv= jjm-j+1
+          ENDIF
        ENDIF
     ENDDO
     !                                 
 
-    IF ( jfiltnu.LE.0 ) jfiltnu=1
     IF( jfiltnu.GT. jjm/2 +1 )  THEN
        PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
@@ -193 +212 @@
     ENDIF
 
-    IF( jfiltsu.LE.0) jfiltsu=1
     IF( jfiltsu.GT.  jjm  +1 )  THEN
        PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
@@ -199 +217 @@
     ENDIF
 
-    IF( jfiltnv.LE.0) jfiltnv=1
     IF( jfiltnv.GT. jjm/2    )  THEN
        PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
@@ -205 +222 @@
     ENDIF
 
-    IF( jfiltsv.LE.0) jfiltsv=1
     IF( jfiltsv.GT.     jjm  )  THEN
        PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
@@ -211 +227 @@
     ENDIF
 
-    PRINT *,' jfiltnv jfiltsv jfiltnu jfiltsu ' , &
+    PRINT *,'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ' , &
          jfiltnv,jfiltsv,jfiltnu,jfiltsu
 
     IF(first_call_inifilr) THEN
        ALLOCATE(matriceun(iim,iim,jfiltnu))
-       ALLOCATE(matriceus(iim,iim,jfiltsu))
+       ALLOCATE(matriceus(iim,iim,jjm-jfiltsu+1))
        ALLOCATE(matricevn(iim,iim,jfiltnv))
-       ALLOCATE(matricevs(iim,iim,jfiltsv))
+       ALLOCATE(matricevs(iim,iim,jjm-jfiltsv+1))
        ALLOCATE( matrinvn(iim,iim,jfiltnu))
-       ALLOCATE( matrinvs(iim,iim,jfiltsu))
+       ALLOCATE( matrinvs(iim,iim,jjm-jfiltsu+1))
        first_call_inifilr = .FALSE.
     ENDIF
@@ -230 +246 @@
     !
     DO j = 1,jjm
+    !default initialization: all modes are retained (i.e. no filtering)
        modfrstu( j ) = iim
        modfrstv( j ) = iim
@@ -306 +323 @@
 
     IF(jfiltnv.GE.jjm/2 .OR. jfiltnu.GE.jjm/2)THEN
-
+! Ehouarn: and what are these for??? Trying to handle a limit case
+!          where filters extend to and meet at the equator?
        IF(jfiltnv.EQ.jfiltsv)jfiltsv=1+jfiltnv
        IF(jfiltnu.EQ.jfiltsu)jfiltsu=1+jfiltnu
@@ -334 +352 @@
              eignft(i,k) = eignfnv(k,i) * coff
           ENDDO
-       ENDDO
+       ENDDO ! of DO i=1,iim
 #ifdef CRAY
        CALL MXM( eignfnv,iim,eignft,iim,matriceun(1,1,j),iim )
@@ -350 +368 @@
              ENDDO
           ENDDO
-       ENDDO
-#endif
-#endif
-
-    ENDDO
+       ENDDO ! of DO k = 1, iim
+#endif
+#endif
+
+    ENDDO ! of DO j = 2, jfiltnu
 
     DO j = jfiltsu, jjm
@@ -364 +382 @@
              eignft(i,k) = eignfnv(k,i) * coff
           ENDDO
-       ENDDO
+       ENDDO ! of DO i=1,iim
 #ifdef CRAY
        CALL MXM(eignfnv,iim,eignft,iim,matriceus(1,1,j-jfiltsu+1),iim)
@@ -381 +399 @@
              ENDDO
           ENDDO
-       ENDDO
-#endif
-#endif
-
-    ENDDO
+       ENDDO ! of DO k = 1, iim
+#endif
+#endif
+
+    ENDDO ! of DO j = jfiltsu, jjm
 
     !   ...................................................................
@@ -421 +439 @@
 #endif
 
-    ENDDO
+    ENDDO ! of DO j = 1, jfiltnv
 
     DO j = jfiltsv, jjm
@@ -452 +470 @@
 #endif
 
-    ENDDO
+    ENDDO ! of DO j = jfiltsv, jjm
 
     !   ...................................................................
@@ -488 +506 @@
 #endif
 
-    ENDDO
+    ENDDO ! of DO j = 2, jfiltnu
 
     DO j = jfiltsu, jjm
@@ -518 +536 @@
 #endif
 
-    ENDDO
+    ENDDO ! of DO j = jfiltsu, jjm
 
     IF (use_filtre_fft) THEN