source: LMDZ6/branches/contrails/libf/filtrez/mod_filtre_fft.F90 @ 5440

Last change on this file since 5440 was 5271, checked in by abarral, 2 months ago

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