source: LMDZ4/trunk/libf/filtrez/filtreg.F @ 5021

Last change on this file since 5021 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210
  • 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
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
[524]83      IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
[1146]84     &     STOP'Pas de transformee simple dans cette version'
85     
[524]86      IF( iter.EQ. 2 )  THEN
[1146]87         PRINT *,' Pas d iteration du filtre dans cette version !'
88     &        , ' Utiliser old_filtreg et repasser !'
89         STOP
[524]90      ENDIF
[1146]91     
[524]92      IF( ifiltre.EQ. -2 .AND..NOT.griscal )     THEN
[1146]93         PRINT *,' Cette routine ne calcule le filtre inverse que '
94     &        , ' sur la grille des scalaires !'
95         STOP
[524]96      ENDIF
[1146]97     
[524]98      IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 )  THEN
[1146]99         PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
100     &        , ' corriger et repasser !'
101         STOP
[524]102      ENDIF
[1146]103     
[524]104      iim2   = iim * iim
105      immjm  = iim * jjm
[1146]106
[524]107      IF( griscal )   THEN
108         IF( nlat. NE. jjp1 )  THEN
[1146]109            PRINT  1111
110            STOP
[524]111         ELSE
[1146]112           
113            IF( iaire.EQ.1 )  THEN
[1279]114               sdd1_type = type_sddv
115               sdd2_type = type_unsddv
[1146]116            ELSE
[1279]117               sdd1_type = type_unsddv
118               sdd2_type = type_sddv
[1146]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
[524]134      ELSE
[1146]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
[524]161      END IF
[1146]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
[524]186#ifdef BLAS
[1146]187                  CALL DGEMM("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)
[524]191#else
[1146]192                  eignq(:,j-jdfil+1,:)
193     $                 = matmul(matrinvn(:,:,j), champ(:iim,j,:))
[524]194#endif
[1146]195               END DO
196               
197            ELSE IF ( griscal )     THEN
198               
199               DO j = jdfil,jffil
[524]200#ifdef BLAS
[1146]201                  CALL DGEMM("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)
[524]205#else
[1146]206                  eignq(:,j-jdfil+1,:)
207     $                 = matmul(matriceun(:,:,j), champ(:iim,j,:))
[524]208#endif
[1146]209               END DO
210               
211            ELSE
212               
213               DO j = jdfil,jffil
[524]214#ifdef BLAS
[1146]215                  CALL DGEMM("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)
[524]219#else
[1146]220                  eignq(:,j-jdfil+1,:)
221     $                 = matmul(matricevn(:,:,j), champ(:iim,j,:))
[524]222#endif
[1146]223               END DO
224               
225            ENDIF
226           
227         ELSE
228           
229            IF( ifiltre. EQ. -2 )   THEN
230               
231               DO j = jdfil,jffil
[524]232#ifdef BLAS
[1146]233                  CALL DGEMM("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)
[524]237#else
[1146]238                  eignq(:,j-jdfil+1,:)
239     $                 = matmul(matrinvs(:,:,j-jfiltsu+1),
240     $                 champ(:iim,j,:))
[524]241#endif
[1146]242               END DO
243               
244               
245            ELSE IF ( griscal )     THEN
246               
247               DO j = jdfil,jffil
[524]248#ifdef BLAS
[1146]249                  CALL DGEMM("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)
[524]253#else
[1146]254                  eignq(:,j-jdfil+1,:)
255     $                 = matmul(matriceus(:,:,j-jfiltsu+1),
256     $                 champ(:iim,j,:))
[524]257#endif
[1146]258               END DO
259                             
260            ELSE
261               
262               DO j = jdfil,jffil
[524]263#ifdef BLAS
[1146]264                  CALL DGEMM("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)
[524]268#else
[1146]269                  eignq(:,j-jdfil+1,:)
270     $                 = matmul(matricevs(:,:,j-jfiltsv+1),
271     $                 champ(:iim,j,:))
[524]272#endif
[1146]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
[524]290
[1146]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
[524]3141111  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
[1146]315     &     filtrer, sur la grille des scalaires'/)
[524]3162222  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
[1146]317     &     ltrer, sur la grille de V ou de Z'/)
[524]318      RETURN
319      END
Note: See TracBrowser for help on using the repository browser.