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

Last change on this file since 3367 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
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
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               
[985]264            ENDIF
[1146]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               
[985]303            ENDIF
[1146]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
[763]311c$OMP END DO NOWAIT
[1146]312           
[985]313ccccccccccccccccccccccccccccccccccccccccccccc
314c Utilisation du filtre FFT
315ccccccccccccccccccccccccccccccccccccccccccccc
316       
[1146]317         ELSE
[985]318       
319c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1146]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
[985]327            ENDDO
328c$OMP END DO NOWAIT
329
[1146]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
[985]339
340
[1146]341            IF( ifiltre.EQ. 2 )  THEN
[985]342c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
[1146]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
[985]351c$OMP END DO NOWAIT       
[1146]352            ELSE
[985]353       
354c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
[1146]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
[985]363c$OMP END DO NOWAIT         
[1146]364            ENDIF
[985]365c
366c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1146]367            DO l=1,nbniv
368               DO j=jdfil,jffil
[985]369!            champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
[1146]370                  champ( iip1,j,l ) = champ( 1,j,l )
371               ENDDO
372            ENDDO
[985]373c$OMP END DO NOWAIT             
[1146]374         ENDIF
[985]375c Fin de la zone de filtrage
376
377       
[1146]378      ENDDO
[985]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     
[763]391c
[1146]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'/)
[985]396c$OMP MASTER     
397      CALL stop_timer
398c$OMP END MASTER
[763]399      RETURN
400      END
Note: See TracBrowser for help on using the repository browser.