source: LMDZ4/trunk/libf/dyn3dpar/filtreg_p.F @ 4005

Last change on this file since 4005 was 1146, checked in by Laurent Fairhead, 16 years ago

Réintegration dans le tronc des modifications issues de la branche LMDZ-dev
comprises entre la révision 1074 et 1145
Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes
résultats entre la trunk et la dev
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.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                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
211     &                    matrinvn(1,1,j), iim,
212     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
213     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
214                  ENDDO
215                 
216               ELSE IF ( griscal )     THEN
217                  DO j = jdfil,jffil
218                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
219     &                    matriceun(1,1,j), iim,
220     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
221     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
222                  ENDDO
223                 
224               ELSE
225                  DO j = jdfil,jffil
226                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
227     &                    matricevn(1,1,j), iim,
228     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
229     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
230                  ENDDO
231                 
232               ENDIF
233               
234            ELSE
235               
236               IF( ifiltre.EQ.-2 )   THEN
237                  DO j = jdfil,jffil
238                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
239     &                    matrinvs(1,1,j-jfiltsu+1), iim,
240     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
241     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
242                  ENDDO
243                 
244               ELSE IF ( griscal )     THEN
245                 
246                  DO j = jdfil,jffil
247                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
248     &                    matriceus(1,1,j-jfiltsu+1), iim,
249     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
250     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
251                  ENDDO
252                 
253               ELSE
254                 
255                  DO j = jdfil,jffil
256                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
257     &                    matricevs(1,1,j-jfiltsv+1), iim,
258     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
259     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
260                  ENDDO
261                 
262               ENDIF
263               
264            ENDIF
265!     c     
266            IF( ifiltre.EQ.2 )  THEN
267               
268c     !-------------------------------------!
269c     ! Dés-agregation des niveau verticaux !
270c     ! uniquement necessaire pour une      !
271c     ! execution OpenMP                    !
272c     !-------------------------------------!
273               ll_nb = 0
274c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
275               DO l = 1, nbniv
276                  ll_nb = ll_nb + 1
277                  DO j = jdfil,jffil
278                     DO i = 1, iim
279                        champ( i,j,l ) = (champ_loc(i,j,ll_nb)
280     &                       + champ_fft(i,j-jdfil+1,ll_nb))
281     &                       * sdd12(i,sdd2_type)
282                     ENDDO
283                  ENDDO
284               ENDDO
285c$OMP END DO NOWAIT
286               
287            ELSE
288               
289               ll_nb = 0
290c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
291               DO l = 1, nbniv_loc
292                  ll_nb = ll_nb + 1
293                  DO j = jdfil,jffil
294                     DO i = 1, iim
295                        champ( i,j,l ) = (champ_loc(i,j,ll_nb)
296     &                       - champ_fft(i,j-jdfil+1,ll_nb))
297     &                       * sdd12(i,sdd2_type)
298                     ENDDO
299                  ENDDO
300               ENDDO
301c$OMP END DO NOWAIT
302               
303            ENDIF
304           
305c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
306            DO l = 1, nbniv
307               DO j = jdfil,jffil
308                  champ( iip1,j,l ) = champ( 1,j,l )
309               ENDDO
310            ENDDO
311c$OMP END DO NOWAIT
312           
313ccccccccccccccccccccccccccccccccccccccccccccc
314c Utilisation du filtre FFT
315ccccccccccccccccccccccccccccccccccccccccccccc
316       
317         ELSE
318       
319c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
320            DO l=1,nbniv
321               DO j=jdfil,jffil
322                  DO  i = 1, iim
323                     champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
324                     champ_fft( i,j,l) = champ(i,j,l)
325                  ENDDO
326               ENDDO
327            ENDDO
328c$OMP END DO NOWAIT
329
330            IF (jdfil<=jffil) THEN
331               IF( ifiltre. EQ. -2 )   THEN
332                  CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
333               ELSE IF ( griscal )     THEN
334                  CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
335               ELSE
336                  CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
337               ENDIF
338            ENDIF
339
340
341            IF( ifiltre.EQ. 2 )  THEN
342c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
343               DO l=1,nbniv
344                  DO j=jdfil,jffil
345                     DO  i = 1, iim
346                        champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
347     &                       *sdd12(i,sdd2_type)
348                     ENDDO
349                  ENDDO
350               ENDDO
351c$OMP END DO NOWAIT       
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)-champ_fft(i,j,l))
359     &                       *sdd12(i,sdd2_type)
360                     ENDDO
361                  ENDDO
362               ENDDO
363c$OMP END DO NOWAIT         
364            ENDIF
365c
366c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
367            DO l=1,nbniv
368               DO j=jdfil,jffil
369!            champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
370                  champ( iip1,j,l ) = champ( 1,j,l )
371               ENDDO
372            ENDDO
373c$OMP END DO NOWAIT             
374         ENDIF
375c Fin de la zone de filtrage
376
377       
378      ENDDO
379
380!      DO j=1,nlat
381!     
382!          PRINT *,"check FFT ----> Delta(",j,")=",
383!     &            sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
384!     &            sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
385!      ENDDO
386     
387!          PRINT *,"check FFT ----> Delta(",j,")=",
388!     &            sum(champ-champ_fft)/sum(champ)
389!     
390     
391c
392 1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
393     &     filtrer, sur la grille des scalaires'/)
394 2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
395     &     ltrer, sur la grille de V ou de Z'/)
396c$OMP MASTER     
397      CALL stop_timer
398c$OMP END MASTER
399      RETURN
400      END
Note: See TracBrowser for help on using the repository browser.