source: LMDZ5/branches/LMDZ5-DOFOCO/libf/dyn3dmem/mod_filtreg_p.F

Last change on this file 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.