source: LMDZ4/branches/LMDZ4_AR5/libf/filtrez/mod_filtre_fft.F90 @ 5192

Last change on this file since 5192 was 1601, checked in by Ehouarn Millour, 13 years ago

Updates and upgrades of filter:

  • add minor corrections (r1591 of trunk) in filtreg_mod.F90: some arrays were oversized and the computation of the indexes from which the filter is applied could go wrong in extreme cases.
  • adapt filtreg_p.F (r1597 of trunk) so that BLAS routine DGEMM is only used if 'BLAS' preprocessing flag is set.
  • update FFT filter routines to match current 'trunk' versions (mostly cosmetic changes, exept for the implementation of use of FFTW in mod_fft_fftw.F90).
  • update "arch-PW6_VARGAS.fcm" to enable use of FFTW

NB: implemented FFTs assume working precision to be double precision ; trying to use them when working precision is single precision will clearly end in tragedy.
EM
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.9 KB
Line 
1!
2! $Id: mod_filtre_fft.F90 1601 2011-12-07 07:05:43Z abarral $
3!
4
5MODULE mod_filtre_fft
6
7  LOGICAL,SAVE :: use_filtre_fft
8  REAL,SAVE,ALLOCATABLE :: Filtre_u(:,:)
9  REAL,SAVE,ALLOCATABLE :: Filtre_v(:,:)
10  REAL,SAVE,ALLOCATABLE :: Filtre_inv(:,:)
11
12CONTAINS
13 
14  SUBROUTINE Init_filtre_fft(coeffu,modfrstu,jfiltnu,jfiltsu,coeffv,modfrstv,jfiltnv,jfiltsv)
15    USE mod_fft
16    IMPLICIT NONE
17    include 'dimensions.h'
18    REAL,   INTENT(IN) :: coeffu(iim,jjm)
19    INTEGER,INTENT(IN) :: modfrstu(jjm)
20    INTEGER,INTENT(IN) :: jfiltnu
21    INTEGER,INTENT(IN) :: jfiltsu
22    REAL,   INTENT(IN) :: coeffv(iim,jjm)
23    INTEGER,INTENT(IN) :: modfrstv(jjm)
24    INTEGER,INTENT(IN) :: jfiltnv
25    INTEGER,INTENT(IN) :: jfiltsv
26   
27    INTEGER            :: index_vp(iim)
28    INTEGER            :: i,j
29    INTEGER            :: l,ll_nb
30
31    index_vp(1)=1
32    DO i=1,iim/2
33      index_vp(i+1)=i*2
34    ENDDO
35   
36    DO i=1,iim/2-1
37      index_vp(iim/2+i+1)=iim-2*i+1
38    ENDDO
39   
40    ALLOCATE(Filtre_u(iim,jjm))
41    ALLOCATE(Filtre_v(iim,jjm))
42    ALLOCATE(Filtre_inv(iim,jjm))
43 
44   
45    DO j=2,jfiltnu
46      DO i=1,iim
47        IF (index_vp(i) < modfrstu(j)) THEN
48          Filtre_u(i,j)=0
49        ELSE
50          Filtre_u(i,j)=coeffu(index_vp(i),j)
51        ENDIF
52      ENDDO
53    ENDDO
54   
55    DO j=jfiltsu,jjm
56      DO i=1,iim
57        IF (index_vp(i) < modfrstu(j)) THEN
58          Filtre_u(i,j)=0
59        ELSE
60          Filtre_u(i,j)=coeffu(index_vp(i),j)
61        ENDIF
62      ENDDO
63    ENDDO
64 
65    DO j=1,jfiltnv
66      DO i=1,iim
67        IF (index_vp(i) < modfrstv(j)) THEN
68          Filtre_v(i,j)=0
69        ELSE
70          Filtre_v(i,j)=coeffv(index_vp(i),j)
71        ENDIF
72      ENDDO
73    ENDDO
74   
75    DO j=jfiltsv,jjm
76      DO i=1,iim
77        IF (index_vp(i) < modfrstv(j)) THEN
78          Filtre_v(i,j)=0
79        ELSE
80          Filtre_v(i,j)=coeffv(index_vp(i),j)
81        ENDIF
82      ENDDO
83    ENDDO
84         
85    DO j=2,jfiltnu
86      DO i=1,iim
87        IF (index_vp(i) < modfrstu(j)) THEN
88          Filtre_inv(i,j)=0
89        ELSE
90          Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
91        ENDIF
92      ENDDO
93    ENDDO
94
95    DO j=jfiltsu,jjm
96      DO i=1,iim
97        IF (index_vp(i) < modfrstu(j)) THEN
98          Filtre_inv(i,j)=0
99        ELSE
100          Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
101        ENDIF
102      ENDDO
103    ENDDO
104   
105#ifdef FFT_FFTW
106
107    WRITE (*,*)"COTH jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv"
108    WRITE (*,*)jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv
109    WRITE (*,*)MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1
110    CALL Init_FFT(iim,(llm+1)*(MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1))
111#else   
112    CALL Init_FFT(iim,(jjm+1)*(llm+1))
113#endif       
114   
115  END SUBROUTINE Init_filtre_fft
116 
117  SUBROUTINE Filtre_u_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
118    USE mod_fft
119#ifdef CPP_PARA
120    USE parallel,ONLY : OMP_CHUNK
121#endif
122    IMPLICIT NONE
123    include 'dimensions.h'
124    INTEGER,INTENT(IN) :: nlat
125    INTEGER,INTENT(IN) :: jj_begin
126    INTEGER,INTENT(IN) :: jj_end
127    INTEGER,INTENT(IN) :: nbniv
128    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
129
130    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
131    COMPLEX         :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
132    INTEGER            :: nb_vect
133    INTEGER :: i,j,l
134    INTEGER :: ll_nb
135   
136    ll_nb=0
137!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
138    DO l=1,nbniv
139      ll_nb=ll_nb+1
140      DO j=1,jj_end-jj_begin+1
141        DO i=1,iim+1
142          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
143        ENDDO
144      ENDDO
145    ENDDO
146!$OMP END DO NOWAIT
147
148    nb_vect=(jj_end-jj_begin+1)*ll_nb
149
150    CALL FFT_forward(vect,TF_vect,nb_vect)
151
152    DO l=1,ll_nb
153      DO j=1,jj_end-jj_begin+1
154        DO i=1,iim/2+1
155          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_u(i,jj_begin+j-1)
156        ENDDO
157      ENDDO
158    ENDDO
159 
160    CALL FFT_backward(TF_vect,vect,nb_vect)
161     
162     
163    ll_nb=0
164!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
165    DO l=1,nbniv
166      ll_nb=ll_nb+1
167      DO j=1,jj_end-jj_begin+1
168        DO i=1,iim+1
169          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
170        ENDDO
171      ENDDO
172    ENDDO
173!$OMP END DO NOWAIT
174
175  END SUBROUTINE Filtre_u_fft
176 
177
178  SUBROUTINE Filtre_v_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
179    USE mod_fft
180#ifdef CPP_PARA
181    USE parallel,ONLY : OMP_CHUNK
182#endif
183    IMPLICIT NONE
184    INCLUDE 'dimensions.h'
185    INTEGER,INTENT(IN) :: nlat
186    INTEGER,INTENT(IN) :: jj_begin
187    INTEGER,INTENT(IN) :: jj_end
188    INTEGER,INTENT(IN) :: nbniv
189    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
190
191    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
192    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
193    INTEGER            :: nb_vect
194    INTEGER :: i,j,l
195    INTEGER :: ll_nb
196   
197    ll_nb=0
198!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
199    DO l=1,nbniv
200      ll_nb=ll_nb+1
201      DO j=1,jj_end-jj_begin+1
202        DO i=1,iim+1
203          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
204        ENDDO
205      ENDDO
206    ENDDO
207!$OMP END DO NOWAIT
208
209   
210    nb_vect=(jj_end-jj_begin+1)*ll_nb
211
212    CALL FFT_forward(vect,TF_vect,nb_vect)
213 
214    DO l=1,ll_nb
215      DO j=1,jj_end-jj_begin+1
216        DO i=1,iim/2+1
217          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_v(i,jj_begin+j-1)
218        ENDDO
219      ENDDO
220    ENDDO
221 
222    CALL FFT_backward(TF_vect,vect,nb_vect)
223   
224   
225    ll_nb=0
226!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
227    DO l=1,nbniv
228      ll_nb=ll_nb+1
229      DO j=1,jj_end-jj_begin+1
230        DO i=1,iim+1
231          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
232        ENDDO
233      ENDDO
234    ENDDO
235!$OMP END DO NOWAIT
236 
237  END SUBROUTINE Filtre_v_fft
238
239
240  SUBROUTINE Filtre_inv_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
241    USE mod_fft
242#ifdef CPP_PARA
243    USE parallel,ONLY : OMP_CHUNK
244#endif
245    IMPLICIT NONE
246    INCLUDE 'dimensions.h'
247    INTEGER,INTENT(IN) :: nlat
248    INTEGER,INTENT(IN) :: jj_begin
249    INTEGER,INTENT(IN) :: jj_end
250    INTEGER,INTENT(IN) :: nbniv
251    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
252
253     REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
254    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
255    INTEGER            :: nb_vect
256    INTEGER :: i,j,l
257    INTEGER :: ll_nb
258   
259    ll_nb=0
260!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
261    DO l=1,nbniv
262      ll_nb=ll_nb+1
263      DO j=1,jj_end-jj_begin+1
264        DO i=1,iim+1
265          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
266        ENDDO
267      ENDDO
268    ENDDO
269!$OMP END DO NOWAIT
270
271    nb_vect=(jj_end-jj_begin+1)*ll_nb
272
273    CALL FFT_forward(vect,TF_vect,nb_vect)
274 
275    DO l=1,ll_nb
276      DO j=1,jj_end-jj_begin+1
277        DO i=1,iim/2+1
278          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_inv(i,jj_begin+j-1)
279        ENDDO
280      ENDDO
281    ENDDO
282 
283    CALL FFT_backward(TF_vect,vect,nb_vect)
284
285    ll_nb=0
286!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
287    DO l=1,nbniv
288      ll_nb=ll_nb+1
289      DO j=1,jj_end-jj_begin+1
290        DO i=1,iim+1
291          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
292        ENDDO
293      ENDDO
294    ENDDO
295!$OMP END DO NOWAIT
296
297  END SUBROUTINE Filtre_inv_fft 
298   
299END MODULE mod_filtre_fft
300 
Note: See TracBrowser for help on using the repository browser.