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

Last change on this file 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
Line 
1
2
3      SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv,
4     &     ifiltre, iaire, griscal ,iter)
5      USE Parallel, only : OMP_CHUNK
6      USE mod_filtre_fft
7      USE timer_filtre
8     
9      USE filtreg_mod
10     
11      IMPLICIT NONE
12     
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
61     
62      REAL  champ( iip1,nlat,nbniv)
63     
64      LOGICAL    griscal
65      INTEGER    hemisph, iaire
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)
72
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
85      IF (first) THEN
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.
93      ENDIF
94
95c$OMP MASTER     
96      CALL start_timer
97c$OMP END MASTER
98
99c-------------------------------------------------------c
100
101      IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
102     &     STOP'Pas de transformee simple dans cette version'
103     
104      IF( iter.EQ. 2 )  THEN
105         PRINT *,' Pas d iteration du filtre dans cette version !'
106     &        , ' Utiliser old_filtreg et repasser !'
107         STOP
108      ENDIF
109
110      IF( ifiltre.EQ. -2 .AND..NOT.griscal )     THEN
111         PRINT *,' Cette routine ne calcule le filtre inverse que '
112     &        , ' sur la grille des scalaires !'
113         STOP
114      ENDIF
115
116      IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 )  THEN
117         PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
118     &        , ' corriger et repasser !'
119         STOP
120      ENDIF
121c
122
123      iim2   = iim * iim
124      immjm  = iim * jjm
125c
126c
127      IF( griscal )   THEN
128         IF( nlat. NE. jjp1 )  THEN
129            PRINT  1111
130            STOP
131         ELSE
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
140c
141            jdfil1 = 2
142            jffil1 = jfiltnu
143            jdfil2 = jfiltsu
144            jffil2 = jjm
145         ENDIF
146      ELSE
147         IF( nlat.NE.jjm )  THEN
148            PRINT  2222
149            STOP
150         ELSE
151c
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
178
179
180cccccccccccccccccccccccccccccccccccccccccccc
181c Utilisation du filtre classique
182cccccccccccccccccccccccccccccccccccccccccccc
183
184         IF (.NOT. use_filtre_fft) THEN
185     
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
203
204            nbniv_loc = ll_nb
205
206            IF( hemisph.EQ.1 )      THEN
207               
208               IF( ifiltre.EQ.-2 )   THEN
209                  DO j = jdfil,jffil
210#ifdef BLAS
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)
215#else
216                     champ_fft(:,j-jdfil+1,:)
217     &                    =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
218#endif
219                  ENDDO
220                 
221               ELSE IF ( griscal )     THEN
222                  DO j = jdfil,jffil
223#ifdef BLAS
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)
228#else
229                     champ_fft(:,j-jdfil+1,:)
230     &                    =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
231#endif
232                  ENDDO
233                 
234               ELSE
235                  DO j = jdfil,jffil
236#ifdef BLAS
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)
241#else
242                     champ_fft(:,j-jdfil+1,:)
243     &                    =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
244#endif
245                  ENDDO
246                 
247               ENDIF
248               
249            ELSE
250               
251               IF( ifiltre.EQ.-2 )   THEN
252                  DO j = jdfil,jffil
253#ifdef BLAS
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)
258#else
259                     champ_fft(:,j-jdfil+1,:)
260     &                    =matmul(matrinvs(:,:,j-jfiltsu+1),
261     &                            champ_loc(:iim,j,:))
262#endif
263                  ENDDO
264                 
265               ELSE IF ( griscal )     THEN
266                 
267                  DO j = jdfil,jffil
268#ifdef BLAS
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)
273#else
274                     champ_fft(:,j-jdfil+1,:)
275     &                    =matmul(matriceus(:,:,j-jfiltsu+1),
276     &                            champ_loc(:iim,j,:))
277#endif
278                  ENDDO
279                 
280               ELSE
281                 
282                  DO j = jdfil,jffil
283#ifdef BLAS
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)
288#else
289                     champ_fft(:,j-jdfil+1,:)
290     &                    =matmul(matricevs(:,:,j-jfiltsv+1),
291     &                            champ_loc(:iim,j,:))
292#endif
293                  ENDDO
294                 
295               ENDIF
296               
297            ENDIF
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               
336            ENDIF
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
344c$OMP END DO NOWAIT
345           
346ccccccccccccccccccccccccccccccccccccccccccccc
347c Utilisation du filtre FFT
348ccccccccccccccccccccccccccccccccccccccccccccc
349       
350         ELSE
351       
352c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
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
360            ENDDO
361c$OMP END DO NOWAIT
362
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
372
373
374            IF( ifiltre.EQ. 2 )  THEN
375c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
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
384c$OMP END DO NOWAIT       
385            ELSE
386       
387c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
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
396c$OMP END DO NOWAIT         
397            ENDIF
398c
399c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
400            DO l=1,nbniv
401               DO j=jdfil,jffil
402!            champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
403                  champ( iip1,j,l ) = champ( 1,j,l )
404               ENDDO
405            ENDDO
406c$OMP END DO NOWAIT             
407         ENDIF
408c Fin de la zone de filtrage
409
410       
411      ENDDO
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     
424c
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'/)
429c$OMP MASTER     
430      CALL stop_timer
431c$OMP END MASTER
432      RETURN
433      END
Note: See TracBrowser for help on using the repository browser.