source: LMDZ6/branches/IPSLCM6.0.14/libf/filtrez/filtreg.F @ 3553

Last change on this file since 3553 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

  • 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
Line 
1!
2! $Header$
3!
4      SUBROUTINE filtreg ( champ, nlat, nbniv, ifiltre,iaire,
5     &     griscal ,iter)
6     
7      USE filtreg_mod
8     
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
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)
60      LOGICAL    griscal
61      INTEGER    hemisph, iaire
62
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
74      IF (first) THEN
75         sdd12(1:iim,type_sddu) = sddu(1:iim)
76         sdd12(1:iim,type_sddv) = sddv(1:iim)
77         sdd12(1:iim,type_unsddu) = unsddu(1:iim)
78         sdd12(1:iim,type_unsddv) = unsddv(1:iim)
79
80         first=.FALSE.
81      ENDIF
82
83      IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
84     &     STOP'Pas de transformee simple dans cette version'
85     
86      IF( iter.EQ. 2 )  THEN
87         PRINT *,' Pas d iteration du filtre dans cette version !'
88     &        , ' Utiliser old_filtreg et repasser !'
89         STOP
90      ENDIF
91     
92      IF( ifiltre.EQ. -2 .AND..NOT.griscal )     THEN
93         PRINT *,' Cette routine ne calcule le filtre inverse que '
94     &        , ' sur la grille des scalaires !'
95         STOP
96      ENDIF
97     
98      IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 )  THEN
99         PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
100     &        , ' corriger et repasser !'
101         STOP
102      ENDIF
103     
104      iim2   = iim * iim
105      immjm  = iim * jjm
106
107      IF( griscal )   THEN
108         IF( nlat. NE. jjp1 )  THEN
109            PRINT  1111
110            STOP
111         ELSE
112           
113            IF( iaire.EQ.1 )  THEN
114               sdd1_type = type_sddv
115               sdd2_type = type_unsddv
116            ELSE
117               sdd1_type = type_unsddv
118               sdd2_type = type_sddv
119            ENDIF
120
121c            IF( iaire.EQ.1 )  THEN
122c               CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 )
123c               CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
124c            ELSE
125c               CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
126c               CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
127c            END IF
128           
129            jdfil1 = 2
130            jffil1 = jfiltnu
131            jdfil2 = jfiltsu
132            jffil2 = jjm
133         END IF
134      ELSE
135         IF( nlat.NE.jjm )  THEN
136            PRINT  2222
137            STOP
138         ELSE
139           
140            IF( iaire.EQ.1 )  THEN
141               sdd1_type = type_sddu
142               sdd2_type = type_unsddu
143            ELSE
144               sdd1_type = type_unsddu
145               sdd2_type = type_sddu
146            ENDIF
147
148c            IF( iaire.EQ.1 )  THEN
149c               CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 )
150c               CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
151c            ELSE
152c               CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
153c               CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
154c            END IF
155           
156            jdfil1 = 1
157            jffil1 = jfiltnv
158            jdfil2 = jfiltsv
159            jffil2 = jjm
160         END IF
161      END IF
162     
163      DO hemisph = 1, 2
164         
165         IF ( hemisph.EQ.1 )  THEN
166            jdfil = jdfil1
167            jffil = jffil1
168         ELSE
169            jdfil = jdfil2
170            jffil = jffil2
171         END IF
172         
173         DO l = 1, nbniv
174            DO j = jdfil,jffil
175               DO i = 1, iim
176                  champ(i,j,l) = champ(i,j,l) * sdd12(i,sdd1_type) ! sdd1(i)
177               END DO
178            END DO
179         END DO
180         
181         IF( hemisph. EQ. 1 )      THEN
182           
183            IF( ifiltre. EQ. -2 )   THEN
184               
185               DO j = jdfil,jffil
186#ifdef BLAS
187                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
188     &                 matrinvn(1,1,j),
189     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
190     &                 eignq(1,j-jdfil+1,1), iim*nlat)
191#else
192                  eignq(:,j-jdfil+1,:)
193     $                 = matmul(matrinvn(:,:,j), champ(:iim,j,:))
194#endif
195               END DO
196               
197            ELSE IF ( griscal )     THEN
198               
199               DO j = jdfil,jffil
200#ifdef BLAS
201                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
202     &                 matriceun(1,1,j),
203     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
204     &                 eignq(1,j-jdfil+1,1), iim*nlat)
205#else
206                  eignq(:,j-jdfil+1,:)
207     $                 = matmul(matriceun(:,:,j), champ(:iim,j,:))
208#endif
209               END DO
210               
211            ELSE
212               
213               DO j = jdfil,jffil
214#ifdef BLAS
215                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
216     &                 matricevn(1,1,j),
217     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
218     &                 eignq(1,j-jdfil+1,1), iim*nlat)
219#else
220                  eignq(:,j-jdfil+1,:)
221     $                 = matmul(matricevn(:,:,j), champ(:iim,j,:))
222#endif
223               END DO
224               
225            ENDIF
226           
227         ELSE
228           
229            IF( ifiltre. EQ. -2 )   THEN
230               
231               DO j = jdfil,jffil
232#ifdef BLAS
233                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
234     &                 matrinvs(1,1,j-jfiltsu+1),
235     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
236     &                 eignq(1,j-jdfil+1,1), iim*nlat)
237#else
238                  eignq(:,j-jdfil+1,:)
239     $                 = matmul(matrinvs(:,:,j-jfiltsu+1),
240     $                 champ(:iim,j,:))
241#endif
242               END DO
243               
244               
245            ELSE IF ( griscal )     THEN
246               
247               DO j = jdfil,jffil
248#ifdef BLAS
249                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
250     &                 matriceus(1,1,j-jfiltsu+1),
251     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
252     &                 eignq(1,j-jdfil+1,1), iim*nlat)
253#else
254                  eignq(:,j-jdfil+1,:)
255     $                 = matmul(matriceus(:,:,j-jfiltsu+1),
256     $                 champ(:iim,j,:))
257#endif
258               END DO
259                             
260            ELSE
261               
262               DO j = jdfil,jffil
263#ifdef BLAS
264                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
265     &                 matricevs(1,1,j-jfiltsv+1),
266     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
267     &                 eignq(1,j-jdfil+1,1), iim*nlat)
268#else
269                  eignq(:,j-jdfil+1,:)
270     $                 = matmul(matricevs(:,:,j-jfiltsv+1),
271     $                 champ(:iim,j,:))
272#endif
273               END DO
274                             
275            ENDIF
276           
277         ENDIF
278         
279         IF( ifiltre.EQ. 2 )  THEN
280           
281            DO l = 1, nbniv
282               DO j = jdfil,jffil
283                  DO i = 1, iim
284                     champ( i,j,l ) =
285     &                    (champ(i,j,l) + eignq(i,j-jdfil+1,l))
286     &                    * sdd12(i,sdd2_type) ! sdd2(i)
287                  END DO
288               END DO
289            END DO
290
291         ELSE
292
293            DO l = 1, nbniv
294               DO j = jdfil,jffil
295                  DO i = 1, iim
296                     champ( i,j,l ) =
297     &                    (champ(i,j,l) - eignq(i,j-jdfil+1,l))
298     &                    * sdd12(i,sdd2_type) ! sdd2(i)
299                  END DO
300               END DO
301            END DO
302
303         ENDIF
304
305         DO l = 1, nbniv
306            DO j = jdfil,jffil
307               champ( iip1,j,l ) = champ( 1,j,l )
308            END DO
309         END DO
310
311     
312      ENDDO
313
3141111  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
315     &     filtrer, sur la grille des scalaires'/)
3162222  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
317     &     ltrer, sur la grille de V ou de Z'/)
318      RETURN
319      END
Note: See TracBrowser for help on using the repository browser.