source: LMDZ5/trunk/libf/dyn3dmem/mod_filtreg_p.F @ 1994

Last change on this file since 1994 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
File size: 13.2 KB
Line 
1      MODULE mod_filtreg_p
2     
3      CONTAINS
4     
5      SUBROUTINE filtreg_p ( champ,jjb,jje, ibeg, iend, nlat, nbniv,
6     &     ifiltre, iaire, griscal ,iter)
7      USE parallel_lmdz, only : OMP_CHUNK
8      USE mod_filtre_fft_loc
9      USE timer_filtre
10     
11      USE filtreg_mod
12     
13      IMPLICIT NONE
14     
15c=======================================================================
16c
17c   Auteur: P. Le Van        07/10/97
18c   ------
19c
20c   Objet: filtre matriciel longitudinal ,avec les matrices precalculees
21c                     pour l'operateur  Filtre    .
22c   ------
23c
24c   Arguments:
25c   ----------
26c
27c     
28c      ibeg..iend            lattitude a filtrer
29c      nlat                  nombre de latitudes du champ
30c      nbniv                 nombre de niveaux verticaux a filtrer
31c      champ(iip1,nblat,nbniv)  en entree : champ a filtrer
32c                            en sortie : champ filtre
33c      ifiltre               +1  Transformee directe
34c                            -1  Transformee inverse
35c                            +2  Filtre directe
36c                            -2  Filtre inverse
37c
38c      iaire                 1   si champ intensif
39c                            2   si champ extensif (pondere par les aires)
40c
41c      iter                  1   filtre simple
42c
43c=======================================================================
44c
45c
46c                      Variable Intensive
47c                ifiltre = 1     filtre directe
48c                ifiltre =-1     filtre inverse
49c
50c                      Variable Extensive
51c                ifiltre = 2     filtre directe
52c                ifiltre =-2     filtre inverse
53c
54c
55#include "dimensions.h"
56#include "paramet.h"
57#include "coefils.h"
58c
59      INTEGER jjb,jje,ibeg,iend,nlat,nbniv,ifiltre,iter
60      INTEGER i,j,l,k
61      INTEGER iim2,immjm
62      INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
63     
64      REAL  champ( iip1,jjb:jje,nbniv)
65     
66      LOGICAL    griscal
67      INTEGER    hemisph, iaire
68     
69      REAL :: champ_fft(iip1,jjb:jje,nbniv)
70      REAL :: champ_in(iip1,jjb:jje,nbniv)
71     
72      LOGICAL,SAVE     :: first=.TRUE.
73c$OMP THREADPRIVATE(first)
74
75      REAL, DIMENSION(iip1,jjb:jje,nbniv) :: champ_loc
76      INTEGER :: ll_nb, nbniv_loc
77      REAL, SAVE :: sdd12(iim,4)
78c$OMP THREADPRIVATE(sdd12)
79
80      INTEGER, PARAMETER :: type_sddu=1
81      INTEGER, PARAMETER :: type_sddv=2
82      INTEGER, PARAMETER :: type_unsddu=3
83      INTEGER, PARAMETER :: type_unsddv=4
84
85      INTEGER :: sdd1_type, sdd2_type
86
87      IF (first) THEN
88         sdd12(1:iim,type_sddu) = sddu(1:iim)
89         sdd12(1:iim,type_sddv) = sddv(1:iim)
90         sdd12(1:iim,type_unsddu) = unsddu(1:iim)
91         sdd12(1:iim,type_unsddv) = unsddv(1:iim)
92
93         CALL Init_timer
94         first=.FALSE.
95      ENDIF
96
97c$OMP MASTER     
98      CALL start_timer
99c$OMP END MASTER
100
101c-------------------------------------------------------c
102
103      IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
104     &     STOP'Pas de transformee simple dans cette version'
105     
106      IF( iter.EQ. 2 )  THEN
107         PRINT *,' Pas d iteration du filtre dans cette version !'
108     &        , ' Utiliser old_filtreg et repasser !'
109         STOP
110      ENDIF
111
112      IF( ifiltre.EQ. -2 .AND..NOT.griscal )     THEN
113         PRINT *,' Cette routine ne calcule le filtre inverse que '
114     &        , ' sur la grille des scalaires !'
115         STOP
116      ENDIF
117
118      IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 )  THEN
119         PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
120     &        , ' corriger et repasser !'
121         STOP
122      ENDIF
123c
124
125      iim2   = iim * iim
126      immjm  = iim * jjm
127c
128c
129      IF( griscal )   THEN
130         IF( nlat. NE. jjp1 )  THEN
131            PRINT  1111
132            STOP
133         ELSE
134c     
135            IF( iaire.EQ.1 )  THEN
136               sdd1_type = type_sddv
137               sdd2_type = type_unsddv
138            ELSE
139               sdd1_type = type_unsddv
140               sdd2_type = type_sddv
141            ENDIF
142c
143            jdfil1 = 2
144            jffil1 = jfiltnu
145            jdfil2 = jfiltsu
146            jffil2 = jjm
147         ENDIF
148      ELSE
149         IF( nlat.NE.jjm )  THEN
150            PRINT  2222
151            STOP
152         ELSE
153c
154            IF( iaire.EQ.1 )  THEN
155               sdd1_type = type_sddu
156               sdd2_type = type_unsddu
157            ELSE
158               sdd1_type = type_unsddu
159               sdd2_type = type_sddu
160            ENDIF
161c     
162            jdfil1 = 1
163            jffil1 = jfiltnv
164            jdfil2 = jfiltsv
165            jffil2 = jjm
166         ENDIF
167      ENDIF
168c     
169      DO hemisph = 1, 2
170c     
171         IF ( hemisph.EQ.1 )  THEN
172cym
173            jdfil = max(jdfil1,ibeg)
174            jffil = min(jffil1,iend)
175         ELSE
176cym
177            jdfil = max(jdfil2,ibeg)
178            jffil = min(jffil2,iend)
179         ENDIF
180
181
182cccccccccccccccccccccccccccccccccccccccccccc
183c Utilisation du filtre classique
184cccccccccccccccccccccccccccccccccccccccccccc
185
186         IF (.NOT. use_filtre_fft) THEN
187     
188c     !---------------------------------!
189c     ! Agregation des niveau verticaux !
190c     ! uniquement necessaire pour une  !
191c     ! execution OpenMP                !
192c     !---------------------------------!
193            ll_nb = 0
194c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
195            DO l = 1, nbniv
196               ll_nb = ll_nb+1
197               DO j = jdfil,jffil
198                  DO i = 1, iim
199                     champ_loc(i,j,ll_nb) =
200     &                    champ(i,j,l) * sdd12(i,sdd1_type)
201                  ENDDO
202               ENDDO
203            ENDDO
204c$OMP END DO NOWAIT
205
206            nbniv_loc = ll_nb
207
208            IF( hemisph.EQ.1 )      THEN
209               
210               IF( ifiltre.EQ.-2 )   THEN
211                  DO j = jdfil,jffil
212#ifdef BLAS
213                     CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
214     &                    matrinvn(1,1,j), iim,
215     &                    champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
216     &                    champ_fft(1,j,1), iip1*(jje-jjb+1))
217#else
218                     champ_fft(:,j,1:nbniv_loc)=
219     &                    matmul(matrinvn(:,:,j),
220     &                    champ_loc(:iim,j,1:nbniv_loc))
221#endif
222                  ENDDO
223                 
224               ELSE IF ( griscal )     THEN
225                  DO j = jdfil,jffil
226#ifdef BLAS
227                     CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
228     &                    matriceun(1,1,j), iim,
229     &                    champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
230     &                    champ_fft(1,j,1), iip1*(jje-jjb+1))
231#else
232                     champ_fft(:,j,1:nbniv_loc)=
233     &                    matmul(matriceun(:,:,j),
234     &                           champ_loc(:iim,j,1:nbniv_loc))
235#endif
236                  ENDDO
237                 
238               ELSE
239                  DO j = jdfil,jffil
240#ifdef BLAS
241                     CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
242     &                    matricevn(1,1,j), iim,
243     &                    champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
244     &                    champ_fft(1,j,1), iip1*(jje-jjb+1))
245#else
246                     champ_fft(:,j,1:nbniv_loc)=
247     &                    matmul(matricevn(:,:,j),           
248     &                           champ_loc(:iim,j,1:nbniv_loc))
249#endif
250                  ENDDO
251                 
252               ENDIF
253               
254            ELSE
255               
256               IF( ifiltre.EQ.-2 )   THEN
257                  DO j = jdfil,jffil
258#ifdef BLAS
259                     CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
260     &                    matrinvs(1,1,j-jfiltsu+1), iim,
261     &                    champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
262     &                    champ_fft(1,j,1), iip1*(jje-jjb+1))
263#else
264                     champ_fft(:,j,1:nbniv_loc)=
265     &                    matmul(matrinvs(:,:,j-jfiltsu+1),
266     &                           champ_loc(:iim,j,1:nbniv_loc))
267#endif
268                  ENDDO
269                 
270               ELSE IF ( griscal )     THEN
271                 
272                  DO j = jdfil,jffil
273#ifdef BLAS
274                     CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
275     &                    matriceus(1,1,j-jfiltsu+1), iim,
276     &                    champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
277     &                    champ_fft(1,j,1), iip1*(jje-jjb+1))
278#else
279                     champ_fft(:,j,1:nbniv_loc)=
280     &                    matmul(matriceus(:,:,j-jfiltsu+1),
281     &                           champ_loc(:iim,j,1:nbniv_loc))
282#endif
283                  ENDDO
284                 
285               ELSE
286                 
287                  DO j = jdfil,jffil
288#ifdef BLAS
289                     CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
290     &                    matricevs(1,1,j-jfiltsv+1), iim,
291     &                    champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
292     &                    champ_fft(1,j,1), iip1*(jje-jjb+1))
293#else
294                     champ_fft(:,j,1:nbniv_loc)=
295     &                    matmul(matricevs(:,:,j-jfiltsv+1),
296     &                           champ_loc(:iim,j,1:nbniv_loc))
297#endif
298                  ENDDO
299                 
300               ENDIF
301               
302            ENDIF
303!     c     
304            IF( ifiltre.EQ.2 )  THEN
305               
306c     !-------------------------------------!
307c     ! Dés-agregation des niveau verticaux !
308c     ! uniquement necessaire pour une      !
309c     ! execution OpenMP                    !
310c     !-------------------------------------!
311               ll_nb = 0
312c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
313               DO l = 1, nbniv
314                  ll_nb = ll_nb + 1
315                  DO j = jdfil,jffil
316                     DO i = 1, iim
317                        champ( i,j,l ) = (champ_loc(i,j,ll_nb)
318     &                       + champ_fft(i,j,ll_nb))
319     &                       * sdd12(i,sdd2_type)
320                     ENDDO
321                  ENDDO
322               ENDDO
323c$OMP END DO NOWAIT
324               
325            ELSE
326               
327               ll_nb = 0
328c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
329               DO l = 1, nbniv
330                  ll_nb = ll_nb + 1
331                  DO j = jdfil,jffil
332                     DO i = 1, iim
333                        champ( i,j,l ) = (champ_loc(i,j,ll_nb)
334     &                       - champ_fft(i,j,ll_nb))
335     &                       * sdd12(i,sdd2_type)
336                     ENDDO
337                  ENDDO
338               ENDDO
339c$OMP END DO NOWAIT
340               
341            ENDIF
342           
343c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
344            DO l = 1, nbniv
345               DO j = jdfil,jffil
346                  champ( iip1,j,l ) = champ( 1,j,l )
347               ENDDO
348            ENDDO
349c$OMP END DO NOWAIT
350           
351ccccccccccccccccccccccccccccccccccccccccccccc
352c Utilisation du filtre FFT
353ccccccccccccccccccccccccccccccccccccccccccccc
354       
355         ELSE
356       
357c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
358            DO l=1,nbniv
359               DO j=jdfil,jffil
360                  DO  i = 1, iim
361                     champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
362                     champ_fft( i,j,l) = champ(i,j,l)
363                  ENDDO
364               ENDDO
365            ENDDO
366c$OMP END DO NOWAIT
367
368            IF (jdfil<=jffil) THEN
369               IF( ifiltre. EQ. -2 )   THEN
370                CALL Filtre_inv_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
371               ELSE IF ( griscal )     THEN
372                  CALL Filtre_u_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
373               ELSE
374                  CALL Filtre_v_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
375               ENDIF
376            ENDIF
377
378
379            IF( ifiltre.EQ. 2 )  THEN
380c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
381               DO l=1,nbniv
382                  DO j=jdfil,jffil
383                     DO  i = 1, iim
384                        champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
385     &                       *sdd12(i,sdd2_type)
386                     ENDDO
387                  ENDDO
388               ENDDO
389c$OMP END DO NOWAIT       
390            ELSE
391       
392c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
393               DO l=1,nbniv
394                  DO j=jdfil,jffil
395                     DO  i = 1, iim
396                        champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
397     &                       *sdd12(i,sdd2_type)
398                     ENDDO
399                  ENDDO
400               ENDDO
401c$OMP END DO NOWAIT         
402            ENDIF
403c
404c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
405            DO l=1,nbniv
406               DO j=jdfil,jffil
407!            champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
408                  champ( iip1,j,l ) = champ( 1,j,l )
409               ENDDO
410            ENDDO
411c$OMP END DO NOWAIT             
412         ENDIF
413c Fin de la zone de filtrage
414
415       
416      ENDDO
417
418!      DO j=1,nlat
419!     
420!          PRINT *,"check FFT ----> Delta(",j,")=",
421!     &            sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
422!     &            sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
423!      ENDDO
424     
425!          PRINT *,"check FFT ----> Delta(",j,")=",
426!     &            sum(champ-champ_fft)/sum(champ)
427!     
428     
429c
430 1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
431     &     filtrer, sur la grille des scalaires'/)
432 2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
433     &     ltrer, sur la grille de V ou de Z'/)
434c$OMP MASTER     
435      CALL stop_timer
436c$OMP END MASTER
437      RETURN
438      END SUBROUTINE filtreg_p
439      END MODULE mod_filtreg_p
440
Note: See TracBrowser for help on using the repository browser.