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

Last change on this file since 5158 was 5134, checked in by abarral, 8 weeks ago

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