source: LMDZ4/branches/LMDZ4_AR5/libf/dyn3dpar/filtreg_p.F @ 1601

Last change on this file since 1601 was 1601, checked in by Ehouarn Millour, 13 years ago

Updates and upgrades of filter:

  • add minor corrections (r1591 of trunk) in filtreg_mod.F90: some arrays were oversized and the computation of the indexes from which the filter is applied could go wrong in extreme cases.
  • adapt filtreg_p.F (r1597 of trunk) so that BLAS routine DGEMM is only used if 'BLAS' preprocessing flag is set.
  • update FFT filter routines to match current 'trunk' versions (mostly cosmetic changes, exept for the implementation of use of FFTW in mod_fft_fftw.F90).
  • update "arch-PW6_VARGAS.fcm" to enable use of FFTW

NB: implemented FFTs assume working precision to be double precision ; trying to use them when working precision is single precision will clearly end in tragedy.
EM

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