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

Last change on this file since 5297 was 5297, checked in by abarral, 29 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
RevLine 
[524]1!
2! $Header$
3!
[5246]4SUBROUTINE filtreg ( champ, nlat, nbniv, ifiltre,iaire, &
5        griscal ,iter)
[524]6
[5246]7  USE filtreg_mod
[1146]8
[5271]9  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5297]10  USE paramet_mod_h
11  USE coefils_mod_h
[5271]12IMPLICIT NONE
[5246]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  !
[5271]50
[5246]51  INTEGER :: nlat,nbniv,ifiltre,iter
52  INTEGER :: i,j,l,k
53  INTEGER :: iim2,immjm
54  INTEGER :: jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
[524]55
[5246]56  REAL :: champ( iip1,nlat,nbniv)
[1146]57
[5246]58  REAL :: eignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim)
59  LOGICAL :: griscal
60  INTEGER :: hemisph, iaire
[1146]61
[5246]62  LOGICAL,SAVE     :: first=.TRUE.
[1146]63
[5246]64  REAL, SAVE :: sdd12(iim,4)
[1146]65
[5246]66  INTEGER, PARAMETER :: type_sddu=1
67  INTEGER, PARAMETER :: type_sddv=2
68  INTEGER, PARAMETER :: type_unsddu=3
69  INTEGER, PARAMETER :: type_unsddv=4
[4440]70
[5246]71  INTEGER :: sdd1_type, sdd2_type
[1146]72
[5246]73  if ( iim == 1 ) return ! no filtre in 2D y-z
[1146]74
[5246]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)
[1146]80
[5246]81     first=.FALSE.
82  ENDIF
[1146]83
[5246]84  IF(ifiltre.EQ.1.or.ifiltre.EQ.-1) &
85        STOP 'Pas de transformee simple dans cette version'
[1146]86
[5246]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
[524]187#ifdef BLAS
[5246]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)
[524]192#else
[5246]193              eignq(:,j-jdfil+1,:) &
194                    = matmul(matrinvn(:,:,j), champ(:iim,j,:))
[524]195#endif
[5246]196           END DO
197
198        ELSE IF ( griscal )     THEN
199
200           DO j = jdfil,jffil
[524]201#ifdef BLAS
[5246]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)
[524]206#else
[5246]207              eignq(:,j-jdfil+1,:) &
208                    = matmul(matriceun(:,:,j), champ(:iim,j,:))
[524]209#endif
[5246]210           END DO
211
212        ELSE
213
214           DO j = jdfil,jffil
[524]215#ifdef BLAS
[5246]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)
[524]220#else
[5246]221              eignq(:,j-jdfil+1,:) &
222                    = matmul(matricevn(:,:,j), champ(:iim,j,:))
[524]223#endif
[5246]224           END DO
225
226        ENDIF
227
228     ELSE
229
230        IF( ifiltre.EQ. -2 )   THEN
231
232           DO j = jdfil,jffil
[524]233#ifdef BLAS
[5246]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)
[524]238#else
[5246]239              eignq(:,j-jdfil+1,:) &
240                    = matmul(matrinvs(:,:,j-jfiltsu+1), &
241                    champ(:iim,j,:))
[524]242#endif
[5246]243           END DO
244
245
246        ELSE IF ( griscal )     THEN
247
248           DO j = jdfil,jffil
[524]249#ifdef BLAS
[5246]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)
[524]254#else
[5246]255              eignq(:,j-jdfil+1,:) &
256                    = matmul(matriceus(:,:,j-jfiltsu+1), &
257                    champ(:iim,j,:))
[524]258#endif
[5246]259           END DO
260
261        ELSE
262
263           DO j = jdfil,jffil
[524]264#ifdef BLAS
[5246]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)
[524]269#else
[5246]270              eignq(:,j-jdfil+1,:) &
271                    = matmul(matricevs(:,:,j-jfiltsv+1), &
272                    champ(:iim,j,:))
[524]273#endif
[5246]274           END DO
[524]275
[5246]276        ENDIF
[1146]277
[5246]278     ENDIF
[1146]279
[5246]280     IF( ifiltre.EQ. 2 )  THEN
[1146]281
[5246]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
[1146]291
[5246]292     ELSE
[1146]293
[5246]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.