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

Last change on this file since 1783 was 1705, checked in by Ehouarn Millour, 12 years ago

Added arch files for ADA (IDRIS IBMx3750) and made the following code modifications:
phylmd/printflag.F : removed "print" of unset variable (radpas0)
dyn3dmem/integrd_loc.F : removed unecessary "include mpif.h"
dyn3dmem/leapfrog_loc.F : removed unecessary "include mpif.h" and allocate saved variables at first call
dyn3dmem/mod_filtreg_p.F : added matmul() alternatives to call to BLAS routine SGEMM (which was incorectly set as DGEMM; which would fail if running with -r4)
filtrez/filtreg.F: changed calls to DGEMM into calls to SGEMM, so that code works with either -r4 or -r8 (the later being used in conjunction with "BLAS SGEMV=DGEMV SGEMM=DGEMM" preprocessing statements)
EM
File size: 13.0 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, 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,:)=
219     &                    matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
220#endif
221                  ENDDO
222                 
223               ELSE IF ( griscal )     THEN
224                  DO j = jdfil,jffil
225#ifdef BLAS
226                     CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
227     &                    matriceun(1,1,j), iim,
228     &                    champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
229     &                    champ_fft(1,j,1), iip1*(jje-jjb+1))
230#else
231                     champ_fft(:,j,:)=
232     &                    matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
233#endif
234                  ENDDO
235                 
236               ELSE
237                  DO j = jdfil,jffil
238#ifdef BLAS
239                     CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
240     &                    matricevn(1,1,j), iim,
241     &                    champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
242     &                    champ_fft(1,j,1), iip1*(jje-jjb+1))
243#else
244                     champ_fft(:,j,:)=
245     &                    matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
246#endif
247                  ENDDO
248                 
249               ENDIF
250               
251            ELSE
252               
253               IF( ifiltre.EQ.-2 )   THEN
254                  DO j = jdfil,jffil
255#ifdef BLAS
256                     CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
257     &                    matrinvs(1,1,j-jfiltsu+1), iim,
258     &                    champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
259     &                    champ_fft(1,j,1), iip1*(jje-jjb+1))
260#else
261                     champ_fft(:,j,:)=
262     &                    matmul(matrinvs(:,:,j-jfiltsu+1),
263     &                           champ_loc(:iim,j,:))
264#endif
265                  ENDDO
266                 
267               ELSE IF ( griscal )     THEN
268                 
269                  DO j = jdfil,jffil
270#ifdef BLAS
271                     CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
272     &                    matriceus(1,1,j-jfiltsu+1), iim,
273     &                    champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
274     &                    champ_fft(1,j,1), iip1*(jje-jjb+1))
275#else
276                     champ_fft(:,j,:)=
277     &                    matmul(matriceus(:,:,j-jfiltsu+1),
278     &                           champ_loc(:iim,j,:))
279#endif
280                  ENDDO
281                 
282               ELSE
283                 
284                  DO j = jdfil,jffil
285#ifdef BLAS
286                     CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
287     &                    matricevs(1,1,j-jfiltsv+1), iim,
288     &                    champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
289     &                    champ_fft(1,j,1), iip1*(jje-jjb+1))
290#else
291                     champ_fft(:,j,:)=
292     &                    matmul(matricevs(:,:,j-jfiltsv+1),
293     &                           champ_loc(:iim,j,:))
294#endif
295                  ENDDO
296                 
297               ENDIF
298               
299            ENDIF
300!     c     
301            IF( ifiltre.EQ.2 )  THEN
302               
303c     !-------------------------------------!
304c     ! Dés-agregation des niveau verticaux !
305c     ! uniquement necessaire pour une      !
306c     ! execution OpenMP                    !
307c     !-------------------------------------!
308               ll_nb = 0
309c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
310               DO l = 1, nbniv_loc
311                  ll_nb = ll_nb + 1
312                  DO j = jdfil,jffil
313                     DO i = 1, iim
314                        champ( i,j,l ) = (champ_loc(i,j,ll_nb)
315     &                       + champ_fft(i,j,ll_nb))
316     &                       * sdd12(i,sdd2_type)
317                     ENDDO
318                  ENDDO
319               ENDDO
320c$OMP END DO NOWAIT
321               
322            ELSE
323               
324               ll_nb = 0
325c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
326               DO l = 1, nbniv_loc
327                  ll_nb = ll_nb + 1
328                  DO j = jdfil,jffil
329                     DO i = 1, iim
330                        champ( i,j,l ) = (champ_loc(i,j,ll_nb)
331     &                       - champ_fft(i,j,ll_nb))
332     &                       * sdd12(i,sdd2_type)
333                     ENDDO
334                  ENDDO
335               ENDDO
336c$OMP END DO NOWAIT
337               
338            ENDIF
339           
340c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
341            DO l = 1, nbniv
342               DO j = jdfil,jffil
343                  champ( iip1,j,l ) = champ( 1,j,l )
344               ENDDO
345            ENDDO
346c$OMP END DO NOWAIT
347           
348ccccccccccccccccccccccccccccccccccccccccccccc
349c Utilisation du filtre FFT
350ccccccccccccccccccccccccccccccccccccccccccccc
351       
352         ELSE
353       
354c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
355            DO l=1,nbniv
356               DO j=jdfil,jffil
357                  DO  i = 1, iim
358                     champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
359                     champ_fft( i,j,l) = champ(i,j,l)
360                  ENDDO
361               ENDDO
362            ENDDO
363c$OMP END DO NOWAIT
364
365            IF (jdfil<=jffil) THEN
366               IF( ifiltre. EQ. -2 )   THEN
367                CALL Filtre_inv_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
368               ELSE IF ( griscal )     THEN
369                  CALL Filtre_u_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
370               ELSE
371                  CALL Filtre_v_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
372               ENDIF
373            ENDIF
374
375
376            IF( ifiltre.EQ. 2 )  THEN
377c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
378               DO l=1,nbniv
379                  DO j=jdfil,jffil
380                     DO  i = 1, iim
381                        champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
382     &                       *sdd12(i,sdd2_type)
383                     ENDDO
384                  ENDDO
385               ENDDO
386c$OMP END DO NOWAIT       
387            ELSE
388       
389c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
390               DO l=1,nbniv
391                  DO j=jdfil,jffil
392                     DO  i = 1, iim
393                        champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
394     &                       *sdd12(i,sdd2_type)
395                     ENDDO
396                  ENDDO
397               ENDDO
398c$OMP END DO NOWAIT         
399            ENDIF
400c
401c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
402            DO l=1,nbniv
403               DO j=jdfil,jffil
404!            champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
405                  champ( iip1,j,l ) = champ( 1,j,l )
406               ENDDO
407            ENDDO
408c$OMP END DO NOWAIT             
409         ENDIF
410c Fin de la zone de filtrage
411
412       
413      ENDDO
414
415!      DO j=1,nlat
416!     
417!          PRINT *,"check FFT ----> Delta(",j,")=",
418!     &            sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
419!     &            sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
420!      ENDDO
421     
422!          PRINT *,"check FFT ----> Delta(",j,")=",
423!     &            sum(champ-champ_fft)/sum(champ)
424!     
425     
426c
427 1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
428     &     filtrer, sur la grille des scalaires'/)
429 2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
430     &     ltrer, sur la grille de V ou de Z'/)
431c$OMP MASTER     
432      CALL stop_timer
433c$OMP END MASTER
434      RETURN
435      END SUBROUTINE filtreg_p
436      END MODULE mod_filtreg_p
437
Note: See TracBrowser for help on using the repository browser.