source: LMDZ6/branches/Amaury_dev/libf/filtrez/lmdz_filtre_fft_loc.F90 @ 5224

Last change on this file since 5224 was 5159, checked in by abarral, 5 months 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
File size: 8.7 KB
Line 
1MODULE lmdz_filtre_fft_loc
2  IMPLICIT NONE; PRIVATE
3  PUBLIC use_filtre_fft, filtre_u, filtre_v, filtre_inv, init_filtre_fft, filtre_u_fft, &
4          filtre_v_fft, filtre_inv_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
15USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
16    IMPLICIT NONE
17
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
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    DO j = 2, jfiltnu
44      DO i = 1, iim
45        IF (index_vp(i) < modfrstu(j)) THEN
46          Filtre_u(i, j) = 0
47        ELSE
48          Filtre_u(i, j) = coeffu(index_vp(i), j)
49        ENDIF
50      ENDDO
51    ENDDO
52
53    DO j = jfiltsu, jjm
54      DO i = 1, iim
55        IF (index_vp(i) < modfrstu(j)) THEN
56          Filtre_u(i, j) = 0
57        ELSE
58          Filtre_u(i, j) = coeffu(index_vp(i), j)
59        ENDIF
60      ENDDO
61    ENDDO
62
63    DO j = 1, jfiltnv
64      DO i = 1, iim
65        IF (index_vp(i) < modfrstv(j)) THEN
66          Filtre_v(i, j) = 0
67        ELSE
68          Filtre_v(i, j) = coeffv(index_vp(i), j)
69        ENDIF
70      ENDDO
71    ENDDO
72
73    DO j = jfiltsv, jjm
74      DO i = 1, iim
75        IF (index_vp(i) < modfrstv(j)) THEN
76          Filtre_v(i, j) = 0
77        ELSE
78          Filtre_v(i, j) = coeffv(index_vp(i), j)
79        ENDIF
80      ENDDO
81    ENDDO
82
83    DO j = 2, jfiltnu
84      DO i = 1, iim
85        IF (index_vp(i) < modfrstu(j)) THEN
86          Filtre_inv(i, j) = 0
87        ELSE
88          Filtre_inv(i, j) = coeffu(index_vp(i), j) / (1. + coeffu(index_vp(i), j))
89        ENDIF
90      ENDDO
91    ENDDO
92
93    DO j = jfiltsu, jjm
94      DO i = 1, iim
95        IF (index_vp(i) < modfrstu(j)) THEN
96          Filtre_inv(i, j) = 0
97        ELSE
98          Filtre_inv(i, j) = coeffu(index_vp(i), j) / (1. + coeffu(index_vp(i), j))
99        ENDIF
100      ENDDO
101    ENDDO
102
103
104    !    CALL Init_FFT(iim,(jjm+1)*(llm+1))
105
106  END SUBROUTINE Init_filtre_fft
107
108  SUBROUTINE Filtre_u_fft(vect_inout, jjb, jje, jj_begin, jj_end, nbniv)
109    USE lmdz_fft
110#ifdef CPP_PARA
111    USE parallel_lmdz,ONLY: OMP_CHUNK
112#endif
113USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
114    IMPLICIT NONE
115
116    INTEGER, INTENT(IN) :: jjb
117    INTEGER, INTENT(IN) :: jje
118    INTEGER, INTENT(IN) :: jj_begin
119    INTEGER, INTENT(IN) :: jj_end
120    INTEGER, INTENT(IN) :: nbniv
121    REAL, INTENT(INOUT) :: vect_inout(iim + 1, jjb:jje, nbniv)
122
123    REAL :: vect(iim + inc, jj_end - jj_begin + 1, nbniv)
124    COMPLEX :: TF_vect(iim / 2 + 1, jj_end - jj_begin + 1, nbniv)
125    INTEGER :: nb_vect
126    INTEGER :: i, j, l
127    INTEGER :: ll_nb
128    !    REAL               :: vect_tmp(iim+inc,jj_end-jj_begin+1,nbniv)
129
130    ll_nb = 0
131    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
132    DO l = 1, nbniv
133      ll_nb = ll_nb + 1
134      DO j = 1, jj_end - jj_begin + 1
135        DO i = 1, iim + 1
136          vect(i, j, ll_nb) = vect_inout(i, j + jj_begin - 1, l)
137        ENDDO
138      ENDDO
139    ENDDO
140    !$OMP END DO NOWAIT
141
142    nb_vect = (jj_end - jj_begin + 1) * ll_nb
143
144    !    vect_tmp=vect
145
146    CALL FFT_forward(vect, TF_vect, nb_vect)
147
148    !    CALL FFT_forward(vect,TF_vect_test,nb_vect)
149    !      PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx"
150    !      DO j=1,jj_end-jj_begin+1
151    !      DO i=1,iim/2+1
152    !         PRINT *,"====",i,j,"----->",TF_vect_test(i,j,1)
153    !       ENDDO
154    !      ENDDO
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    !    CALL FFT_backward(TF_vect_test,vect_test,nb_vect)
166
167    !      PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx"
168    !      DO j=1,jj_end-jj_begin+1
169    !         DO i=1,iim
170    !           PRINT *,"====",i,j,"----->",vect_test(i,j,1)
171    !         ENDDO
172    !      ENDDO
173
174    ll_nb = 0
175    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
176    DO l = 1, nbniv
177      ll_nb = ll_nb + 1
178      DO j = 1, jj_end - jj_begin + 1
179        DO i = 1, iim + 1
180          vect_inout(i, j + jj_begin - 1, l) = vect(i, j, ll_nb)
181        ENDDO
182      ENDDO
183    ENDDO
184    !$OMP END DO NOWAIT
185
186  END SUBROUTINE Filtre_u_fft
187
188
189  SUBROUTINE Filtre_v_fft(vect_inout, jjb, jje, jj_begin, jj_end, nbniv)
190    USE lmdz_fft
191#ifdef CPP_PARA
192    USE parallel_lmdz,ONLY: OMP_CHUNK
193#endif
194USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
195    IMPLICIT NONE
196
197    INTEGER, INTENT(IN) :: jjb
198    INTEGER, INTENT(IN) :: jje
199    INTEGER, INTENT(IN) :: jj_begin
200    INTEGER, INTENT(IN) :: jj_end
201    INTEGER, INTENT(IN) :: nbniv
202    REAL, INTENT(INOUT) :: vect_inout(iim + 1, jjb:jje, nbniv)
203
204    REAL :: vect(iim + inc, jj_end - jj_begin + 1, nbniv)
205    COMPLEX :: TF_vect(iim / 2 + 1, jj_end - jj_begin + 1, nbniv)
206    INTEGER :: nb_vect
207    INTEGER :: i, j, l
208    INTEGER :: ll_nb
209
210    ll_nb = 0
211    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
212    DO l = 1, nbniv
213      ll_nb = ll_nb + 1
214      DO j = 1, jj_end - jj_begin + 1
215        DO i = 1, iim + 1
216          vect(i, j, ll_nb) = vect_inout(i, j + jj_begin - 1, l)
217        ENDDO
218      ENDDO
219    ENDDO
220    !$OMP END DO NOWAIT
221
222    nb_vect = (jj_end - jj_begin + 1) * ll_nb
223
224    CALL FFT_forward(vect, TF_vect, nb_vect)
225
226    DO l = 1, ll_nb
227      DO j = 1, jj_end - jj_begin + 1
228        DO i = 1, iim / 2 + 1
229          TF_vect(i, j, l) = TF_vect(i, j, l) * Filtre_v(i, jj_begin + j - 1)
230        ENDDO
231      ENDDO
232    ENDDO
233
234    CALL FFT_backward(TF_vect, vect, nb_vect)
235
236    ll_nb = 0
237    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
238    DO l = 1, nbniv
239      ll_nb = ll_nb + 1
240      DO j = 1, jj_end - jj_begin + 1
241        DO i = 1, iim + 1
242          vect_inout(i, j + jj_begin - 1, l) = vect(i, j, ll_nb)
243        ENDDO
244      ENDDO
245    ENDDO
246    !$OMP END DO NOWAIT
247
248  END SUBROUTINE Filtre_v_fft
249
250
251  SUBROUTINE Filtre_inv_fft(vect_inout, jjb, jje, jj_begin, jj_end, nbniv)
252    USE lmdz_fft
253#ifdef CPP_PARA
254    USE parallel_lmdz,ONLY: OMP_CHUNK
255#endif
256USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
257    IMPLICIT NONE
258
259    INTEGER, INTENT(IN) :: jjb
260    INTEGER, INTENT(IN) :: jje
261    INTEGER, INTENT(IN) :: jj_begin
262    INTEGER, INTENT(IN) :: jj_end
263    INTEGER, INTENT(IN) :: nbniv
264    REAL, INTENT(INOUT) :: vect_inout(iim + 1, jjb:jje, nbniv)
265
266    REAL :: vect(iim + inc, jj_end - jj_begin + 1, nbniv)
267    COMPLEX :: TF_vect(iim / 2 + 1, jj_end - jj_begin + 1, nbniv)
268    INTEGER :: nb_vect
269    INTEGER :: i, j, l
270    INTEGER :: ll_nb
271
272    ll_nb = 0
273    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
274    DO l = 1, nbniv
275      ll_nb = ll_nb + 1
276      DO j = 1, jj_end - jj_begin + 1
277        DO i = 1, iim + 1
278          vect(i, j, ll_nb) = vect_inout(i, j + jj_begin - 1, l)
279        ENDDO
280      ENDDO
281    ENDDO
282    !$OMP END DO NOWAIT
283
284    nb_vect = (jj_end - jj_begin + 1) * ll_nb
285
286    CALL FFT_forward(vect, TF_vect, nb_vect)
287
288    DO l = 1, ll_nb
289      DO j = 1, jj_end - jj_begin + 1
290        DO i = 1, iim / 2 + 1
291          TF_vect(i, j, l) = TF_vect(i, j, l) * Filtre_inv(i, jj_begin + j - 1)
292        ENDDO
293      ENDDO
294    ENDDO
295
296    CALL FFT_backward(TF_vect, vect, nb_vect)
297
298    ll_nb = 0
299    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
300    DO l = 1, nbniv
301      ll_nb = ll_nb + 1
302      DO j = 1, jj_end - jj_begin + 1
303        DO i = 1, iim + 1
304          vect_inout(i, j + jj_begin - 1, l) = vect(i, j, ll_nb)
305        ENDDO
306      ENDDO
307    ENDDO
308    !$OMP END DO NOWAIT
309
310  END SUBROUTINE Filtre_inv_fft
311
312
313  !  SUBROUTINE get_ll_index(nbniv,ll_index,ll_nb)
314  !  IMPLICIT NONE
315  !    INTEGER,INTENT(IN)  :: nbniv
316  !    INTEGER,INTENT(OUT) :: ll_index(nbniv)
317  !    INTEGER,INTENT(OUT) :: ll_nb
318
319  !    INTEGER :: l,ll_begin, ll_end
320  !   INTEGER :: omp_rank,omp_size
321  !   INTEGER :: OMP_GET_NUM_THREADS
322  !   INTEGER :: omp_chunk
323  !   EXTERNAL OMP_GET_NUM_THREADS
324  !   INTEGER :: OMP_GET_THREAD_NUM
325  !   EXTERNAL OMP_GET_THREAD_NUM
326
327
328  !   omp_size=OMP_GET_NUM_THREADS()
329  !   omp_rank=OMP_GET_THREAD_NUM()
330  !   omp_chunk=nbniv/omp_size+min(1,MOD(nbniv,omp_size))
331
332  !   ll_begin=omp_rank*OMP_CHUNK+1
333  !   ll_nb=0
334  !   DO WHILE (ll_begin<=nbniv)
335  !     ll_end=min(ll_begin+OMP_CHUNK-1,nbniv)
336  !     DO l=ll_begin,ll_end
337  !       ll_nb=ll_nb+1
338  !       ll_index(ll_nb)=l
339  !     ENDDO
340  !     ll_begin=ll_begin+omp_size*OMP_CHUNK
341  !   ENDDO
342
343  !  END SUBROUTINE get_ll_index
344
345END MODULE lmdz_filtre_fft_loc
346 
Note: See TracBrowser for help on using the repository browser.