source: LMDZ6/branches/Amaury_dev/libf/filtrez/lmdz_filtre_fft.F90

Last change on this file was 5159, checked in by abarral, 7 weeks ago

Put dimensions.h and paramet.h into modules
  • 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: 7.1 KB
Line 
1
2! $Id: lmdz_filtre_fft.F90 5159 2024-08-02 19:58:25Z dcugnet $
3
4MODULE lmdz_filtre_fft
5  IMPLICIT NONE; PRIVATE
6  PUBLIC use_filtre_fft, filtre_u, filtre_v, filtre_inv, init_filtre_fft, filtre_u_fft, &
7          filtre_v_fft, filtre_inv_fft
8
9  LOGICAL :: use_filtre_fft
10  REAL,ALLOCATABLE :: filtre_u(:,:)
11  REAL,ALLOCATABLE :: filtre_v(:,:)
12  REAL,ALLOCATABLE :: filtre_inv(:,:)
13
14CONTAINS
15 
16  SUBROUTINE init_filtre_fft(coeffu,modfrstu,jfiltnu,jfiltsu,coeffv,modfrstv,jfiltnv,jfiltsv)
17    USE lmdz_fft
18USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
19    IMPLICIT NONE
20
21    REAL,   INTENT(IN) :: coeffu(iim,jjm)
22    INTEGER,INTENT(IN) :: modfrstu(jjm)
23    INTEGER,INTENT(IN) :: jfiltnu
24    INTEGER,INTENT(IN) :: jfiltsu
25    REAL,   INTENT(IN) :: coeffv(iim,jjm)
26    INTEGER,INTENT(IN) :: modfrstv(jjm)
27    INTEGER,INTENT(IN) :: jfiltnv
28    INTEGER,INTENT(IN) :: jfiltsv
29   
30    INTEGER            :: index_vp(iim)
31    INTEGER            :: i,j
32    INTEGER            :: l,ll_nb
33
34    index_vp(1)=1
35    DO i=1,iim/2
36      index_vp(i+1)=i*2
37    ENDDO
38   
39    DO i=1,iim/2-1
40      index_vp(iim/2+i+1)=iim-2*i+1
41    ENDDO
42   
43    ALLOCATE(filtre_u(iim,jjm))
44    ALLOCATE(filtre_v(iim,jjm))
45    ALLOCATE(filtre_inv(iim,jjm))
46 
47   
48    DO j=2,jfiltnu
49      DO i=1,iim
50        IF (index_vp(i) < modfrstu(j)) THEN
51          filtre_u(i,j)=0
52        ELSE
53          filtre_u(i,j)=coeffu(index_vp(i),j)
54        ENDIF
55      ENDDO
56    ENDDO
57   
58    DO j=jfiltsu,jjm
59      DO i=1,iim
60        IF (index_vp(i) < modfrstu(j)) THEN
61          filtre_u(i,j)=0
62        ELSE
63          filtre_u(i,j)=coeffu(index_vp(i),j)
64        ENDIF
65      ENDDO
66    ENDDO
67 
68    DO j=1,jfiltnv
69      DO i=1,iim
70        IF (index_vp(i) < modfrstv(j)) THEN
71          filtre_v(i,j)=0
72        ELSE
73          filtre_v(i,j)=coeffv(index_vp(i),j)
74        ENDIF
75      ENDDO
76    ENDDO
77   
78    DO j=jfiltsv,jjm
79      DO i=1,iim
80        IF (index_vp(i) < modfrstv(j)) THEN
81          filtre_v(i,j)=0
82        ELSE
83          filtre_v(i,j)=coeffv(index_vp(i),j)
84        ENDIF
85      ENDDO
86    ENDDO
87         
88    DO j=2,jfiltnu
89      DO i=1,iim
90        IF (index_vp(i) < modfrstu(j)) THEN
91          filtre_inv(i,j)=0
92        ELSE
93          filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
94        ENDIF
95      ENDDO
96    ENDDO
97
98    DO j=jfiltsu,jjm
99      DO i=1,iim
100        IF (index_vp(i) < modfrstu(j)) THEN
101          filtre_inv(i,j)=0
102        ELSE
103          filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
104        ENDIF
105      ENDDO
106    ENDDO
107   
108#ifdef FFT_FFTW
109
110    WRITE (*,*)"COTH jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv"
111    WRITE (*,*)jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv
112    WRITE (*,*)MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1
113    CALL init_FFT(iim,(llm+1)*(MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1))
114#else   
115    CALL init_FFT(iim,(jjm+1)*(llm+1))
116#endif       
117   
118  END SUBROUTINE init_filtre_fft
119 
120  SUBROUTINE filtre_u_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
121    USE lmdz_fft
122#ifdef CPP_PARA
123    USE parallel_lmdz,ONLY: OMP_CHUNK
124#endif
125USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
126    IMPLICIT NONE
127
128    INTEGER,INTENT(IN) :: nlat
129    INTEGER,INTENT(IN) :: jj_begin
130    INTEGER,INTENT(IN) :: jj_end
131    INTEGER,INTENT(IN) :: nbniv
132    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
133
134    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
135    COMPLEX         :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
136    INTEGER            :: nb_vect
137    INTEGER :: i,j,l
138    INTEGER :: ll_nb
139   
140    ll_nb=0
141!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
142    DO l=1,nbniv
143      ll_nb=ll_nb+1
144      DO j=1,jj_end-jj_begin+1
145        DO i=1,iim+1
146          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
147        ENDDO
148      ENDDO
149    ENDDO
150!$OMP END DO NOWAIT
151
152    nb_vect=(jj_end-jj_begin+1)*ll_nb
153
154    CALL FFT_forward(vect,TF_vect,nb_vect)
155
156    DO l=1,ll_nb
157      DO j=1,jj_end-jj_begin+1
158        DO i=1,iim/2+1
159          TF_vect(i,j,l)=TF_vect(i,j,l)*filtre_u(i,jj_begin+j-1)
160        ENDDO
161      ENDDO
162    ENDDO
163 
164    CALL FFT_backward(TF_vect,vect,nb_vect)
165     
166     
167    ll_nb=0
168!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
169    DO l=1,nbniv
170      ll_nb=ll_nb+1
171      DO j=1,jj_end-jj_begin+1
172        DO i=1,iim+1
173          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
174        ENDDO
175      ENDDO
176    ENDDO
177!$OMP END DO NOWAIT
178
179  END SUBROUTINE filtre_u_fft
180
181  SUBROUTINE filtre_v_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
182    USE lmdz_fft
183#ifdef CPP_PARA
184    USE parallel_lmdz,ONLY: OMP_CHUNK
185#endif
186USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
187    IMPLICIT NONE
188
189    INTEGER,INTENT(IN) :: nlat
190    INTEGER,INTENT(IN) :: jj_begin
191    INTEGER,INTENT(IN) :: jj_end
192    INTEGER,INTENT(IN) :: nbniv
193    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
194
195    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
196    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
197    INTEGER            :: nb_vect
198    INTEGER :: i,j,l
199    INTEGER :: ll_nb
200   
201    ll_nb=0
202!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
203    DO l=1,nbniv
204      ll_nb=ll_nb+1
205      DO j=1,jj_end-jj_begin+1
206        DO i=1,iim+1
207          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
208        ENDDO
209      ENDDO
210    ENDDO
211!$OMP END DO NOWAIT
212
213   
214    nb_vect=(jj_end-jj_begin+1)*ll_nb
215
216    CALL FFT_forward(vect,TF_vect,nb_vect)
217 
218    DO l=1,ll_nb
219      DO j=1,jj_end-jj_begin+1
220        DO i=1,iim/2+1
221          TF_vect(i,j,l)=TF_vect(i,j,l)*filtre_v(i,jj_begin+j-1)
222        ENDDO
223      ENDDO
224    ENDDO
225 
226    CALL FFT_backward(TF_vect,vect,nb_vect)
227   
228   
229    ll_nb=0
230!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
231    DO l=1,nbniv
232      ll_nb=ll_nb+1
233      DO j=1,jj_end-jj_begin+1
234        DO i=1,iim+1
235          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
236        ENDDO
237      ENDDO
238    ENDDO
239!$OMP END DO NOWAIT
240 
241  END SUBROUTINE filtre_v_fft
242
243  SUBROUTINE filtre_inv_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
244    USE lmdz_fft
245#ifdef CPP_PARA
246    USE parallel_lmdz,ONLY: OMP_CHUNK
247#endif
248USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
249    IMPLICIT 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 lmdz_filtre_fft
304 
Note: See TracBrowser for help on using the repository browser.