source: LMDZ6/trunk/libf/filtrez/filtreg.F90 @ 5297

Last change on this file since 5297 was 5297, checked in by abarral, 35 hours ago

Turn gradsdef.h coefils.h into a module

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