source: LMDZ6/branches/Amaury_dev/libf/filtrez/mod_filtre_fft.F90 @ 5112

Last change on this file since 5112 was 5107, checked in by abarral, 4 months ago

Turn coefils.h into lmdz_coefils.f90
Put filtreg.F90, inifgn.f90, jacobi.F90, eigen_sort.f90, acc.f90 inside lmdz_filtreg.F90
Turn mod_* into lmdz_* in filtrez
Delete obsolete parafilt.h*
(lint) remove spaces between routine name and args
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • 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 5107 2024-07-24 08:26:10Z abarral $
3
4MODULE mod_filtre_fft
5
6  LOGICAL,SAVE :: use_filtre_fft
7  REAL,SAVE,ALLOCATABLE :: Filtre_u(:,:)
8  REAL,SAVE,ALLOCATABLE :: Filtre_v(:,:)
9  REAL,SAVE,ALLOCATABLE :: Filtre_inv(:,:)
10
11CONTAINS
12 
13  SUBROUTINE Init_filtre_fft(coeffu,modfrstu,jfiltnu,jfiltsu,coeffv,modfrstv,jfiltnv,jfiltsv)
14    USE lmdz_fft
15    IMPLICIT NONE
16    include 'dimensions.h'
17    REAL,   INTENT(IN) :: coeffu(iim,jjm)
18    INTEGER,INTENT(IN) :: modfrstu(jjm)
19    INTEGER,INTENT(IN) :: jfiltnu
20    INTEGER,INTENT(IN) :: jfiltsu
21    REAL,   INTENT(IN) :: coeffv(iim,jjm)
22    INTEGER,INTENT(IN) :: modfrstv(jjm)
23    INTEGER,INTENT(IN) :: jfiltnv
24    INTEGER,INTENT(IN) :: jfiltsv
25   
26    INTEGER            :: index_vp(iim)
27    INTEGER            :: i,j
28    INTEGER            :: l,ll_nb
29
30    index_vp(1)=1
31    DO i=1,iim/2
32      index_vp(i+1)=i*2
33    ENDDO
34   
35    DO i=1,iim/2-1
36      index_vp(iim/2+i+1)=iim-2*i+1
37    ENDDO
38   
39    ALLOCATE(Filtre_u(iim,jjm))
40    ALLOCATE(Filtre_v(iim,jjm))
41    ALLOCATE(Filtre_inv(iim,jjm))
42 
43   
44    DO j=2,jfiltnu
45      DO i=1,iim
46        IF (index_vp(i) < modfrstu(j)) THEN
47          Filtre_u(i,j)=0
48        ELSE
49          Filtre_u(i,j)=coeffu(index_vp(i),j)
50        ENDIF
51      ENDDO
52    ENDDO
53   
54    DO j=jfiltsu,jjm
55      DO i=1,iim
56        IF (index_vp(i) < modfrstu(j)) THEN
57          Filtre_u(i,j)=0
58        ELSE
59          Filtre_u(i,j)=coeffu(index_vp(i),j)
60        ENDIF
61      ENDDO
62    ENDDO
63 
64    DO j=1,jfiltnv
65      DO i=1,iim
66        IF (index_vp(i) < modfrstv(j)) THEN
67          Filtre_v(i,j)=0
68        ELSE
69          Filtre_v(i,j)=coeffv(index_vp(i),j)
70        ENDIF
71      ENDDO
72    ENDDO
73   
74    DO j=jfiltsv,jjm
75      DO i=1,iim
76        IF (index_vp(i) < modfrstv(j)) THEN
77          Filtre_v(i,j)=0
78        ELSE
79          Filtre_v(i,j)=coeffv(index_vp(i),j)
80        ENDIF
81      ENDDO
82    ENDDO
83         
84    DO j=2,jfiltnu
85      DO i=1,iim
86        IF (index_vp(i) < modfrstu(j)) THEN
87          Filtre_inv(i,j)=0
88        ELSE
89          Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
90        ENDIF
91      ENDDO
92    ENDDO
93
94    DO j=jfiltsu,jjm
95      DO i=1,iim
96        IF (index_vp(i) < modfrstu(j)) THEN
97          Filtre_inv(i,j)=0
98        ELSE
99          Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
100        ENDIF
101      ENDDO
102    ENDDO
103   
104#ifdef FFT_FFTW
105
106    WRITE (*,*)"COTH jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv"
107    WRITE (*,*)jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv
108    WRITE (*,*)MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1
109    CALL Init_FFT(iim,(llm+1)*(MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1))
110#else   
111    CALL Init_FFT(iim,(jjm+1)*(llm+1))
112#endif       
113   
114  END SUBROUTINE Init_filtre_fft
115 
116  SUBROUTINE Filtre_u_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
117    USE lmdz_fft
118#ifdef CPP_PARA
119    USE parallel_lmdz,ONLY: OMP_CHUNK
120#endif
121    IMPLICIT NONE
122    include 'dimensions.h'
123    INTEGER,INTENT(IN) :: nlat
124    INTEGER,INTENT(IN) :: jj_begin
125    INTEGER,INTENT(IN) :: jj_end
126    INTEGER,INTENT(IN) :: nbniv
127    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
128
129    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
130    COMPLEX         :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
131    INTEGER            :: nb_vect
132    INTEGER :: i,j,l
133    INTEGER :: ll_nb
134   
135    ll_nb=0
136!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
137    DO l=1,nbniv
138      ll_nb=ll_nb+1
139      DO j=1,jj_end-jj_begin+1
140        DO i=1,iim+1
141          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
142        ENDDO
143      ENDDO
144    ENDDO
145!$OMP END DO NOWAIT
146
147    nb_vect=(jj_end-jj_begin+1)*ll_nb
148
149    CALL FFT_forward(vect,TF_vect,nb_vect)
150
151    DO l=1,ll_nb
152      DO j=1,jj_end-jj_begin+1
153        DO i=1,iim/2+1
154          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_u(i,jj_begin+j-1)
155        ENDDO
156      ENDDO
157    ENDDO
158 
159    CALL FFT_backward(TF_vect,vect,nb_vect)
160     
161     
162    ll_nb=0
163!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
164    DO l=1,nbniv
165      ll_nb=ll_nb+1
166      DO j=1,jj_end-jj_begin+1
167        DO i=1,iim+1
168          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
169        ENDDO
170      ENDDO
171    ENDDO
172!$OMP END DO NOWAIT
173
174  END SUBROUTINE Filtre_u_fft
175 
176
177  SUBROUTINE Filtre_v_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
178    USE lmdz_fft
179#ifdef CPP_PARA
180    USE parallel_lmdz,ONLY: OMP_CHUNK
181#endif
182    IMPLICIT NONE
183    INCLUDE 'dimensions.h'
184    INTEGER,INTENT(IN) :: nlat
185    INTEGER,INTENT(IN) :: jj_begin
186    INTEGER,INTENT(IN) :: jj_end
187    INTEGER,INTENT(IN) :: nbniv
188    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
189
190    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
191    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
192    INTEGER            :: nb_vect
193    INTEGER :: i,j,l
194    INTEGER :: ll_nb
195   
196    ll_nb=0
197!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
198    DO l=1,nbniv
199      ll_nb=ll_nb+1
200      DO j=1,jj_end-jj_begin+1
201        DO i=1,iim+1
202          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
203        ENDDO
204      ENDDO
205    ENDDO
206!$OMP END DO NOWAIT
207
208   
209    nb_vect=(jj_end-jj_begin+1)*ll_nb
210
211    CALL FFT_forward(vect,TF_vect,nb_vect)
212 
213    DO l=1,ll_nb
214      DO j=1,jj_end-jj_begin+1
215        DO i=1,iim/2+1
216          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_v(i,jj_begin+j-1)
217        ENDDO
218      ENDDO
219    ENDDO
220 
221    CALL FFT_backward(TF_vect,vect,nb_vect)
222   
223   
224    ll_nb=0
225!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
226    DO l=1,nbniv
227      ll_nb=ll_nb+1
228      DO j=1,jj_end-jj_begin+1
229        DO i=1,iim+1
230          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
231        ENDDO
232      ENDDO
233    ENDDO
234!$OMP END DO NOWAIT
235 
236  END SUBROUTINE Filtre_v_fft
237
238
239  SUBROUTINE Filtre_inv_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
240    USE lmdz_fft
241#ifdef CPP_PARA
242    USE parallel_lmdz,ONLY: OMP_CHUNK
243#endif
244    IMPLICIT NONE
245    INCLUDE 'dimensions.h'
246    INTEGER,INTENT(IN) :: nlat
247    INTEGER,INTENT(IN) :: jj_begin
248    INTEGER,INTENT(IN) :: jj_end
249    INTEGER,INTENT(IN) :: nbniv
250    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
251
252     REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
253    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
254    INTEGER            :: nb_vect
255    INTEGER :: i,j,l
256    INTEGER :: ll_nb
257   
258    ll_nb=0
259!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
260    DO l=1,nbniv
261      ll_nb=ll_nb+1
262      DO j=1,jj_end-jj_begin+1
263        DO i=1,iim+1
264          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
265        ENDDO
266      ENDDO
267    ENDDO
268!$OMP END DO NOWAIT
269
270    nb_vect=(jj_end-jj_begin+1)*ll_nb
271
272    CALL FFT_forward(vect,TF_vect,nb_vect)
273 
274    DO l=1,ll_nb
275      DO j=1,jj_end-jj_begin+1
276        DO i=1,iim/2+1
277          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_inv(i,jj_begin+j-1)
278        ENDDO
279      ENDDO
280    ENDDO
281 
282    CALL FFT_backward(TF_vect,vect,nb_vect)
283
284    ll_nb=0
285!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
286    DO l=1,nbniv
287      ll_nb=ll_nb+1
288      DO j=1,jj_end-jj_begin+1
289        DO i=1,iim+1
290          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
291        ENDDO
292      ENDDO
293    ENDDO
294!$OMP END DO NOWAIT
295
296  END SUBROUTINE Filtre_inv_fft 
297   
298END MODULE mod_filtre_fft
299 
Note: See TracBrowser for help on using the repository browser.