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

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

Mise a jour de dyn3dpar par rapport a dyn3d, inclusion OpenMP et filtre FFT YM
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.0 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      IMPLICIT NONE
9
10c=======================================================================
11c
12c   Auteur: P. Le Van        07/10/97
13c   ------
14c
15c   Objet: filtre matriciel longitudinal ,avec les matrices precalculees
16c                     pour l'operateur  Filtre    .
17c   ------
18c
19c   Arguments:
20c   ----------
21c
22c     
23c      ibeg..iend            lattitude a filtrer
24c      nlat                  nombre de latitudes du champ
25c      nbniv                 nombre de niveaux verticaux a filtrer
26c      champ(iip1,nblat,nbniv)  en entree : champ a filtrer
27c                            en sortie : champ filtre
28c      ifiltre               +1  Transformee directe
29c                            -1  Transformee inverse
30c                            +2  Filtre directe
31c                            -2  Filtre inverse
32c
33c      iaire                 1   si champ intensif
34c                            2   si champ extensif (pondere par les aires)
35c
36c      iter                  1   filtre simple
37c
38c=======================================================================
39c
40c
41c                      Variable Intensive
42c                ifiltre = 1     filtre directe
43c                ifiltre =-1     filtre inverse
44c
45c                      Variable Extensive
46c                ifiltre = 2     filtre directe
47c                ifiltre =-2     filtre inverse
48c
49c
50#include "dimensions.h"
51#include "paramet.h"
52#include "parafilt.h"
53#include "coefils.h"
54c
55      INTEGER ibeg,iend,nlat,nbniv,ifiltre,iter
56      INTEGER i,j,l,k
57      INTEGER iim2,immjm
58      INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
59
60      REAL  champ( iip1,nlat,nbniv)
61      REAL matriceun,matriceus,matricevn,matricevs,matrinvn,matrinvs
62      COMMON/matrfil/matriceun(iim,iim,nfilun),matriceus(iim,iim,nfilus)
63     ,             , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs)
64     ,             ,  matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus)
65cym      REAL  eignq(iim), sdd1(iim),sdd2(iim)
66
67      REAL  eignq(iim)
68      REAL :: sdd1(iim),sdd2(iim)
69     
70      LOGICAL    griscal
71      INTEGER    hemisph, iaire
72     
73      REAL :: champ_fft(iip1,nlat,nbniv)
74      REAL :: champ_in(iip1,nlat,nbniv)
75     
76      REAL,SAVE,TARGET :: sddu_loc(iim)
77      REAL,SAVE,TARGET :: sddv_loc(iim)
78      REAL,SAVE,TARGET :: unsddu_loc(iim)
79      REAL,SAVE,TARGET :: unsddv_loc(iim)
80c$OMP THREADPRIVATE(sddu_loc,sddv_loc,unsddu_loc,unsddv_loc)
81      LOGICAL,SAVE     :: first=.TRUE.
82c$OMP THREADPRIVATE(first)
83
84      IF (first) THEN
85        sddu_loc(1:iim)=sddu(1:iim)
86        sddv_loc(1:iim)=sddv(1:iim)
87        unsddu_loc(1:iim)=unsddu(1:iim)
88        unsddv_loc(1:iim)=unsddv(1:iim)
89        CALL Init_timer
90        first=.FALSE.
91c       PRINT *,"----> sddu_loc=",sddu_loc
92c       PRINT *,"----> sddv_loc=",sddv_loc
93c       PRINT *,"----> unsddu_loc=",unsddu_loc
94c       PRINT *,"----> unsddv_loc=",unsddv_loc
95      ENDIF
96
97c$OMP MASTER     
98      CALL start_timer
99c$OMP END MASTER
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
134cym                CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 )
135cym                CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
136cym               sdd1=>sddv_loc
137cym               sdd2=>unsddv_loc
138               sdd1(1:iim)=sddv_loc(1:iim)
139               sdd2(1:iim)=unsddv_loc(1:iim)
140             ELSE
141cym                CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
142cym                CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
143               sdd1(1:iim)=unsddv_loc(1:iim)
144               sdd2(1:iim)=sddv_loc(1:iim)
145             END IF
146c
147             jdfil1 = 2
148             jffil1 = jfiltnu
149             jdfil2 = jfiltsu
150             jffil2 = jjm
151          END IF
152      ELSE
153          IF( nlat.NE.jjm )  THEN
154             PRINT  2222
155             STOP
156          ELSE
157c
158             IF( iaire.EQ.1 )  THEN
159cym                CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 )
160cym                CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
161cym                sdd1=>sddu_loc
162cym                sdd2=>unsddu_loc
163                sdd1(1:iim)=sddu_loc(1:iim)
164                sdd2(1:iim)=unsddu_loc(1:iim)
165
166             ELSE
167cym                CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
168cym                CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
169cym               sdd1=>unsddu_loc
170cym               sdd2=>sddu_loc
171               sdd1(1:iim)=unsddu_loc(1:iim)
172               sdd2(1:iim)=sddu_loc(1:iim)
173             END IF
174c
175             jdfil1 = 1
176             jffil1 = jfiltnv
177             jdfil2 = jfiltsv
178             jffil2 = jjm
179          END IF
180      END IF
181
182c      PRINT *,"APPEL a filtreg --> sdd1=",sdd1
183c      PRINT *,"APPEL a filtreg --> sdd2=",sdd2
184c      PRINT *,"----> sddu_loc=",sddu_loc
185c       PRINT *,"----> sddv_loc=",sddv_loc
186c       PRINT *,"----> unsddu_loc=",unsddu_loc
187c       PRINT *,"----> unsddv_loc=",unsddv_loc
188 
189c
190c
191      DO 100  hemisph = 1, 2
192c
193      IF ( hemisph.EQ.1 )  THEN
194c ym
195          jdfil = max(jdfil1,ibeg)
196          jffil = min(jffil1,iend)
197      ELSE
198c ym
199          jdfil = max(jdfil2,ibeg)
200          jffil = min(jffil2,iend)
201      END IF
202
203
204cccccccccccccccccccccccccccccccccccccccccccc
205c Utilisation du filtre classique
206cccccccccccccccccccccccccccccccccccccccccccc
207
208      IF (.NOT. use_filtre_fft) THEN
209     
210c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
211      DO 50  l = 1, nbniv
212        DO 30  j = jdfil,jffil
213 
214 
215          DO  5  i = 1, iim
216            champ(i,j,l) = champ(i,j,l) * sdd1(i)
217   5      CONTINUE
218c
219
220          IF( hemisph. EQ. 1 )      THEN
221
222            IF( ifiltre. EQ. -2 )   THEN
223
224
225              CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
226     .                     champ(1,j,l), 1, 0.0, eignq, 1)
227
228
229            ELSE IF ( griscal )     THEN
230
231              CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
232     .                    champ(1,j,l), 1, 0.0, eignq, 1)
233
234            ELSE
235
236              CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
237     .                   champ(1,j,l), 1, 0.0, eignq, 1)
238            ENDIF
239
240          ELSE
241
242            IF( ifiltre. EQ. -2 )   THEN
243     
244              CALL SGEMV("N",iim,iim,1.0, matrinvs(1,1,j-jfiltsu+1),iim,
245     .                   champ(1,j,l), 1, 0.0, eignq, 1)
246     
247            ELSE IF ( griscal )     THEN
248     
249              CALL SGEMV("N",iim,iim,1.0,matriceus(1,1,j-jfiltsu+1),iim,
250     .                   champ(1,j,l), 1, 0.0, eignq, 1)
251            ELSE
252         
253              CALL SGEMV("N",iim,iim,1.0,matricevs(1,1,j-jfiltsv+1),iim,
254     .                    champ(1,j,l), 1, 0.0, eignq, 1)
255            ENDIF
256
257          ENDIF
258
259
260c
261          IF( ifiltre.EQ. 2 )  THEN
262         
263            DO 15 i = 1, iim
264              champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
265  15        CONTINUE
266         
267          ELSE
268       
269            DO 16 i=1,iim
270               champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
27116          CONTINUE
272         
273          ENDIF
274c
275          champ( iip1,j,l ) = champ( 1,j,l )
276c
277  30    CONTINUE
278c
279  50  CONTINUE
280c$OMP END DO NOWAIT
281
282ccccccccccccccccccccccccccccccccccccccccccccc
283c Utilisation du filtre FFT
284ccccccccccccccccccccccccccccccccccccccccccccc
285       
286       ELSE
287       
288c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
289          DO l=1,nbniv
290            DO j=jdfil,jffil
291              DO  i = 1, iim
292                champ( i,j,l)= champ(i,j,l)*sdd1(i)
293                champ_fft( i,j,l) = champ(i,j,l)
294              ENDDO
295            ENDDO
296          ENDDO
297c$OMP END DO NOWAIT
298
299      IF (jdfil<=jffil) THEN
300        IF( ifiltre. EQ. -2 )   THEN
301          CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
302        ELSE IF ( griscal )     THEN
303          CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
304        ELSE
305          CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
306        ENDIF
307      ENDIF
308
309
310        IF( ifiltre.EQ. 2 )  THEN
311c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
312          DO l=1,nbniv
313            DO j=jdfil,jffil
314              DO  i = 1, iim
315                champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
316     .                             *sdd2(i)
317              ENDDO
318            ENDDO
319          ENDDO
320c$OMP END DO NOWAIT       
321        ELSE
322       
323c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
324          DO l=1,nbniv
325            DO j=jdfil,jffil
326              DO  i = 1, iim
327                champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
328     .                            *sdd2(i)
329              ENDDO
330            ENDDO
331          ENDDO
332c$OMP END DO NOWAIT         
333        ENDIF
334c
335c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
336        DO l=1,nbniv
337          DO j=jdfil,jffil
338!            champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
339            champ( iip1,j,l ) = champ( 1,j,l )
340          ENDDO
341        ENDDO
342c$OMP END DO NOWAIT             
343      ENDIF
344c Fin de la zone de filtrage
345
346       
347 100  CONTINUE
348
349!      DO j=1,nlat
350!     
351!          PRINT *,"check FFT ----> Delta(",j,")=",
352!     &            sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
353!     &            sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
354!      ENDDO
355     
356!          PRINT *,"check FFT ----> Delta(",j,")=",
357!     &            sum(champ-champ_fft)/sum(champ)
358!     
359     
360c
3611111  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
362     *filtrer, sur la grille des scalaires'/)
3632222  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
364     *ltrer, sur la grille de V ou de Z'/)
365c$OMP MASTER     
366      CALL stop_timer
367c$OMP END MASTER
368      RETURN
369      END
Note: See TracBrowser for help on using the repository browser.