source: LMDZ6/branches/LMDZ_ECRad/libf/filtrez/filtreg.F @ 4890

Last change on this file since 4890 was 4727, checked in by idelkadi, 14 months ago

Merged trunk changes -r4488:4726 LMDZ_ECRad branch

  • 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.3 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 ( iim == 1 ) return ! no filtre in 2D y-z
75
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
85      IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
86     &     STOP'Pas de transformee simple dans cette version'
87     
88      IF( iter.EQ. 2 )  THEN
89         PRINT *,' Pas d iteration du filtre dans cette version !'
90     &        , ' Utiliser old_filtreg et repasser !'
91         STOP
92      ENDIF
93     
94      IF( ifiltre.EQ. -2 .AND..NOT.griscal )     THEN
95         PRINT *,' Cette routine ne calcule le filtre inverse que '
96     &        , ' sur la grille des scalaires !'
97         STOP
98      ENDIF
99     
100      IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 )  THEN
101         PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
102     &        , ' corriger et repasser !'
103         STOP
104      ENDIF
105     
106      iim2   = iim * iim
107      immjm  = iim * jjm
108
109      IF( griscal )   THEN
110         IF( nlat. NE. jjp1 )  THEN
111            PRINT  1111
112            STOP
113         ELSE
114           
115            IF( iaire.EQ.1 )  THEN
116               sdd1_type = type_sddv
117               sdd2_type = type_unsddv
118            ELSE
119               sdd1_type = type_unsddv
120               sdd2_type = type_sddv
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
136      ELSE
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
163      END IF
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
188#ifdef BLAS
189                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
190     &                 matrinvn(1,1,j),
191     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
192     &                 eignq(1,j-jdfil+1,1), iim*nlat)
193#else
194                  eignq(:,j-jdfil+1,:)
195     $                 = matmul(matrinvn(:,:,j), champ(:iim,j,:))
196#endif
197               END DO
198               
199            ELSE IF ( griscal )     THEN
200               
201               DO j = jdfil,jffil
202#ifdef BLAS
203                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
204     &                 matriceun(1,1,j),
205     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
206     &                 eignq(1,j-jdfil+1,1), iim*nlat)
207#else
208                  eignq(:,j-jdfil+1,:)
209     $                 = matmul(matriceun(:,:,j), champ(:iim,j,:))
210#endif
211               END DO
212               
213            ELSE
214               
215               DO j = jdfil,jffil
216#ifdef BLAS
217                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
218     &                 matricevn(1,1,j),
219     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
220     &                 eignq(1,j-jdfil+1,1), iim*nlat)
221#else
222                  eignq(:,j-jdfil+1,:)
223     $                 = matmul(matricevn(:,:,j), champ(:iim,j,:))
224#endif
225               END DO
226               
227            ENDIF
228           
229         ELSE
230           
231            IF( ifiltre. EQ. -2 )   THEN
232               
233               DO j = jdfil,jffil
234#ifdef BLAS
235                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
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)
239#else
240                  eignq(:,j-jdfil+1,:)
241     $                 = matmul(matrinvs(:,:,j-jfiltsu+1),
242     $                 champ(:iim,j,:))
243#endif
244               END DO
245               
246               
247            ELSE IF ( griscal )     THEN
248               
249               DO j = jdfil,jffil
250#ifdef BLAS
251                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
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)
255#else
256                  eignq(:,j-jdfil+1,:)
257     $                 = matmul(matriceus(:,:,j-jfiltsu+1),
258     $                 champ(:iim,j,:))
259#endif
260               END DO
261                             
262            ELSE
263               
264               DO j = jdfil,jffil
265#ifdef BLAS
266                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
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)
270#else
271                  eignq(:,j-jdfil+1,:)
272     $                 = matmul(matricevs(:,:,j-jfiltsv+1),
273     $                 champ(:iim,j,:))
274#endif
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
292
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
3161111  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
317     &     filtrer, sur la grille des scalaires'/)
3182222  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
319     &     ltrer, sur la grille de V ou de Z'/)
320      RETURN
321      END
Note: See TracBrowser for help on using the repository browser.