source: LMDZ6/trunk/libf/filtrez/filtreg.F @ 4544

Last change on this file since 4544 was 4440, checked in by fhourdin, 22 months ago

Deactivating filtre and check if iim=1 (2D y-z)
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.2 KB
RevLine 
[524]1!
2! $Header$
3!
4      SUBROUTINE filtreg ( champ, nlat, nbniv, ifiltre,iaire,
[1146]5     &     griscal ,iter)
6     
7      USE filtreg_mod
8     
[524]9      IMPLICIT NONE
10c=======================================================================
11c
12c   Auteur: P. Le Van        07/10/97
13c   ------
14c
15c   Objet: filtre matriciel longitudinal ,avec les matrices precalculees
16c                     pour l'operateur  Filtre    .
17c   ------
18c
19c   Arguments:
20c   ----------
21c
22c      nblat                 nombre de latitudes a filtrer
23c      nbniv                 nombre de niveaux verticaux a filtrer
24c      champ(iip1,nblat,nbniv)  en entree : champ a filtrer
25c                            en sortie : champ filtre
26c      ifiltre               +1  Transformee directe
27c                            -1  Transformee inverse
28c                            +2  Filtre directe
29c                            -2  Filtre inverse
30c
31c      iaire                 1   si champ intensif
32c                            2   si champ extensif (pondere par les aires)
33c
34c      iter                  1   filtre simple
35c
36c=======================================================================
37c
38c
39c                      Variable Intensive
40c                ifiltre = 1     filtre directe
41c                ifiltre =-1     filtre inverse
42c
43c                      Variable Extensive
44c                ifiltre = 2     filtre directe
45c                ifiltre =-2     filtre inverse
46c
47c
48#include "dimensions.h"
49#include "paramet.h"
50#include "coefils.h"
51
[1146]52      INTEGER    nlat,nbniv,ifiltre,iter
53      INTEGER    i,j,l,k
54      INTEGER    iim2,immjm
55      INTEGER    jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
56
57      REAL       champ( iip1,nlat,nbniv)
58
59      REAL       eignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim)
[524]60      LOGICAL    griscal
61      INTEGER    hemisph, iaire
62
[1146]63      LOGICAL,SAVE     :: first=.TRUE.
64
65      REAL, SAVE :: sdd12(iim,4)
66
67      INTEGER, PARAMETER :: type_sddu=1
68      INTEGER, PARAMETER :: type_sddv=2
69      INTEGER, PARAMETER :: type_unsddu=3
70      INTEGER, PARAMETER :: type_unsddv=4
71
72      INTEGER :: sdd1_type, sdd2_type
73
[4440]74      if ( iim == 1 ) return ! no filtre in 2D y-z
75
[1146]76      IF (first) THEN
77         sdd12(1:iim,type_sddu) = sddu(1:iim)
78         sdd12(1:iim,type_sddv) = sddv(1:iim)
79         sdd12(1:iim,type_unsddu) = unsddu(1:iim)
80         sdd12(1:iim,type_unsddv) = unsddv(1:iim)
81
82         first=.FALSE.
83      ENDIF
84
[524]85      IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
[1146]86     &     STOP'Pas de transformee simple dans cette version'
87     
[524]88      IF( iter.EQ. 2 )  THEN
[1146]89         PRINT *,' Pas d iteration du filtre dans cette version !'
90     &        , ' Utiliser old_filtreg et repasser !'
91         STOP
[524]92      ENDIF
[1146]93     
[524]94      IF( ifiltre.EQ. -2 .AND..NOT.griscal )     THEN
[1146]95         PRINT *,' Cette routine ne calcule le filtre inverse que '
96     &        , ' sur la grille des scalaires !'
97         STOP
[524]98      ENDIF
[1146]99     
[524]100      IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 )  THEN
[1146]101         PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
102     &        , ' corriger et repasser !'
103         STOP
[524]104      ENDIF
[1146]105     
[524]106      iim2   = iim * iim
107      immjm  = iim * jjm
[1146]108
[524]109      IF( griscal )   THEN
110         IF( nlat. NE. jjp1 )  THEN
[1146]111            PRINT  1111
112            STOP
[524]113         ELSE
[1146]114           
115            IF( iaire.EQ.1 )  THEN
[1279]116               sdd1_type = type_sddv
117               sdd2_type = type_unsddv
[1146]118            ELSE
[1279]119               sdd1_type = type_unsddv
120               sdd2_type = type_sddv
[1146]121            ENDIF
122
123c            IF( iaire.EQ.1 )  THEN
124c               CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 )
125c               CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
126c            ELSE
127c               CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
128c               CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
129c            END IF
130           
131            jdfil1 = 2
132            jffil1 = jfiltnu
133            jdfil2 = jfiltsu
134            jffil2 = jjm
135         END IF
[524]136      ELSE
[1146]137         IF( nlat.NE.jjm )  THEN
138            PRINT  2222
139            STOP
140         ELSE
141           
142            IF( iaire.EQ.1 )  THEN
143               sdd1_type = type_sddu
144               sdd2_type = type_unsddu
145            ELSE
146               sdd1_type = type_unsddu
147               sdd2_type = type_sddu
148            ENDIF
149
150c            IF( iaire.EQ.1 )  THEN
151c               CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 )
152c               CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
153c            ELSE
154c               CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
155c               CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
156c            END IF
157           
158            jdfil1 = 1
159            jffil1 = jfiltnv
160            jdfil2 = jfiltsv
161            jffil2 = jjm
162         END IF
[524]163      END IF
[1146]164     
165      DO hemisph = 1, 2
166         
167         IF ( hemisph.EQ.1 )  THEN
168            jdfil = jdfil1
169            jffil = jffil1
170         ELSE
171            jdfil = jdfil2
172            jffil = jffil2
173         END IF
174         
175         DO l = 1, nbniv
176            DO j = jdfil,jffil
177               DO i = 1, iim
178                  champ(i,j,l) = champ(i,j,l) * sdd12(i,sdd1_type) ! sdd1(i)
179               END DO
180            END DO
181         END DO
182         
183         IF( hemisph. EQ. 1 )      THEN
184           
185            IF( ifiltre. EQ. -2 )   THEN
186               
187               DO j = jdfil,jffil
[524]188#ifdef BLAS
[1705]189                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
[1146]190     &                 matrinvn(1,1,j),
191     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
192     &                 eignq(1,j-jdfil+1,1), iim*nlat)
[524]193#else
[1146]194                  eignq(:,j-jdfil+1,:)
195     $                 = matmul(matrinvn(:,:,j), champ(:iim,j,:))
[524]196#endif
[1146]197               END DO
198               
199            ELSE IF ( griscal )     THEN
200               
201               DO j = jdfil,jffil
[524]202#ifdef BLAS
[1705]203                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
[1146]204     &                 matriceun(1,1,j),
205     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
206     &                 eignq(1,j-jdfil+1,1), iim*nlat)
[524]207#else
[1146]208                  eignq(:,j-jdfil+1,:)
209     $                 = matmul(matriceun(:,:,j), champ(:iim,j,:))
[524]210#endif
[1146]211               END DO
212               
213            ELSE
214               
215               DO j = jdfil,jffil
[524]216#ifdef BLAS
[1705]217                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
[1146]218     &                 matricevn(1,1,j),
219     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
220     &                 eignq(1,j-jdfil+1,1), iim*nlat)
[524]221#else
[1146]222                  eignq(:,j-jdfil+1,:)
223     $                 = matmul(matricevn(:,:,j), champ(:iim,j,:))
[524]224#endif
[1146]225               END DO
226               
227            ENDIF
228           
229         ELSE
230           
231            IF( ifiltre. EQ. -2 )   THEN
232               
233               DO j = jdfil,jffil
[524]234#ifdef BLAS
[1705]235                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
[1146]236     &                 matrinvs(1,1,j-jfiltsu+1),
237     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
238     &                 eignq(1,j-jdfil+1,1), iim*nlat)
[524]239#else
[1146]240                  eignq(:,j-jdfil+1,:)
241     $                 = matmul(matrinvs(:,:,j-jfiltsu+1),
242     $                 champ(:iim,j,:))
[524]243#endif
[1146]244               END DO
245               
246               
247            ELSE IF ( griscal )     THEN
248               
249               DO j = jdfil,jffil
[524]250#ifdef BLAS
[1705]251                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
[1146]252     &                 matriceus(1,1,j-jfiltsu+1),
253     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
254     &                 eignq(1,j-jdfil+1,1), iim*nlat)
[524]255#else
[1146]256                  eignq(:,j-jdfil+1,:)
257     $                 = matmul(matriceus(:,:,j-jfiltsu+1),
258     $                 champ(:iim,j,:))
[524]259#endif
[1146]260               END DO
261                             
262            ELSE
263               
264               DO j = jdfil,jffil
[524]265#ifdef BLAS
[1705]266                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
[1146]267     &                 matricevs(1,1,j-jfiltsv+1),
268     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
269     &                 eignq(1,j-jdfil+1,1), iim*nlat)
[524]270#else
[1146]271                  eignq(:,j-jdfil+1,:)
272     $                 = matmul(matricevs(:,:,j-jfiltsv+1),
273     $                 champ(:iim,j,:))
[524]274#endif
[1146]275               END DO
276                             
277            ENDIF
278           
279         ENDIF
280         
281         IF( ifiltre.EQ. 2 )  THEN
282           
283            DO l = 1, nbniv
284               DO j = jdfil,jffil
285                  DO i = 1, iim
286                     champ( i,j,l ) =
287     &                    (champ(i,j,l) + eignq(i,j-jdfil+1,l))
288     &                    * sdd12(i,sdd2_type) ! sdd2(i)
289                  END DO
290               END DO
291            END DO
[524]292
[1146]293         ELSE
294
295            DO l = 1, nbniv
296               DO j = jdfil,jffil
297                  DO i = 1, iim
298                     champ( i,j,l ) =
299     &                    (champ(i,j,l) - eignq(i,j-jdfil+1,l))
300     &                    * sdd12(i,sdd2_type) ! sdd2(i)
301                  END DO
302               END DO
303            END DO
304
305         ENDIF
306
307         DO l = 1, nbniv
308            DO j = jdfil,jffil
309               champ( iip1,j,l ) = champ( 1,j,l )
310            END DO
311         END DO
312
313     
314      ENDDO
315
[524]3161111  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
[1146]317     &     filtrer, sur la grille des scalaires'/)
[524]3182222  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
[1146]319     &     ltrer, sur la grille de V ou de Z'/)
[524]320      RETURN
321      END
Note: See TracBrowser for help on using the repository browser.