source: LMDZ6/branches/Amaury_dev/libf/filtrez/lmdz_filtre_fft.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
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.0 KB
Line 
1
2! $Id: lmdz_filtre_fft.F90 5134 2024-07-26 15:56:37Z abarral $
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
18    IMPLICIT NONE
19    INCLUDE 'dimensions.h'
20    REAL,   INTENT(IN) :: coeffu(iim,jjm)
21    INTEGER,INTENT(IN) :: modfrstu(jjm)
22    INTEGER,INTENT(IN) :: jfiltnu
23    INTEGER,INTENT(IN) :: jfiltsu
24    REAL,   INTENT(IN) :: coeffv(iim,jjm)
25    INTEGER,INTENT(IN) :: modfrstv(jjm)
26    INTEGER,INTENT(IN) :: jfiltnv
27    INTEGER,INTENT(IN) :: jfiltsv
28   
29    INTEGER            :: index_vp(iim)
30    INTEGER            :: i,j
31    INTEGER            :: l,ll_nb
32
33    index_vp(1)=1
34    DO i=1,iim/2
35      index_vp(i+1)=i*2
36    ENDDO
37   
38    DO i=1,iim/2-1
39      index_vp(iim/2+i+1)=iim-2*i+1
40    ENDDO
41   
42    ALLOCATE(filtre_u(iim,jjm))
43    ALLOCATE(filtre_v(iim,jjm))
44    ALLOCATE(filtre_inv(iim,jjm))
45 
46   
47    DO j=2,jfiltnu
48      DO i=1,iim
49        IF (index_vp(i) < modfrstu(j)) THEN
50          filtre_u(i,j)=0
51        ELSE
52          filtre_u(i,j)=coeffu(index_vp(i),j)
53        ENDIF
54      ENDDO
55    ENDDO
56   
57    DO j=jfiltsu,jjm
58      DO i=1,iim
59        IF (index_vp(i) < modfrstu(j)) THEN
60          filtre_u(i,j)=0
61        ELSE
62          filtre_u(i,j)=coeffu(index_vp(i),j)
63        ENDIF
64      ENDDO
65    ENDDO
66 
67    DO j=1,jfiltnv
68      DO i=1,iim
69        IF (index_vp(i) < modfrstv(j)) THEN
70          filtre_v(i,j)=0
71        ELSE
72          filtre_v(i,j)=coeffv(index_vp(i),j)
73        ENDIF
74      ENDDO
75    ENDDO
76   
77    DO j=jfiltsv,jjm
78      DO i=1,iim
79        IF (index_vp(i) < modfrstv(j)) THEN
80          filtre_v(i,j)=0
81        ELSE
82          filtre_v(i,j)=coeffv(index_vp(i),j)
83        ENDIF
84      ENDDO
85    ENDDO
86         
87    DO j=2,jfiltnu
88      DO i=1,iim
89        IF (index_vp(i) < modfrstu(j)) THEN
90          filtre_inv(i,j)=0
91        ELSE
92          filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
93        ENDIF
94      ENDDO
95    ENDDO
96
97    DO j=jfiltsu,jjm
98      DO i=1,iim
99        IF (index_vp(i) < modfrstu(j)) THEN
100          filtre_inv(i,j)=0
101        ELSE
102          filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
103        ENDIF
104      ENDDO
105    ENDDO
106   
107#ifdef FFT_FFTW
108
109    WRITE (*,*)"COTH jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv"
110    WRITE (*,*)jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv
111    WRITE (*,*)MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1
112    CALL init_FFT(iim,(llm+1)*(MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1))
113#else   
114    CALL init_FFT(iim,(jjm+1)*(llm+1))
115#endif       
116   
117  END SUBROUTINE init_filtre_fft
118 
119  SUBROUTINE filtre_u_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
120    USE lmdz_fft
121#ifdef CPP_PARA
122    USE parallel_lmdz,ONLY: OMP_CHUNK
123#endif
124    IMPLICIT NONE
125    INCLUDE 'dimensions.h'
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  SUBROUTINE filtre_v_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
180    USE lmdz_fft
181#ifdef CPP_PARA
182    USE parallel_lmdz,ONLY: OMP_CHUNK
183#endif
184    IMPLICIT NONE
185    INCLUDE 'dimensions.h'
186    INTEGER,INTENT(IN) :: nlat
187    INTEGER,INTENT(IN) :: jj_begin
188    INTEGER,INTENT(IN) :: jj_end
189    INTEGER,INTENT(IN) :: nbniv
190    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
191
192    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
193    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
194    INTEGER            :: nb_vect
195    INTEGER :: i,j,l
196    INTEGER :: ll_nb
197   
198    ll_nb=0
199!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
200    DO l=1,nbniv
201      ll_nb=ll_nb+1
202      DO j=1,jj_end-jj_begin+1
203        DO i=1,iim+1
204          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
205        ENDDO
206      ENDDO
207    ENDDO
208!$OMP END DO NOWAIT
209
210   
211    nb_vect=(jj_end-jj_begin+1)*ll_nb
212
213    CALL FFT_forward(vect,TF_vect,nb_vect)
214 
215    DO l=1,ll_nb
216      DO j=1,jj_end-jj_begin+1
217        DO i=1,iim/2+1
218          TF_vect(i,j,l)=TF_vect(i,j,l)*filtre_v(i,jj_begin+j-1)
219        ENDDO
220      ENDDO
221    ENDDO
222 
223    CALL FFT_backward(TF_vect,vect,nb_vect)
224   
225   
226    ll_nb=0
227!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
228    DO l=1,nbniv
229      ll_nb=ll_nb+1
230      DO j=1,jj_end-jj_begin+1
231        DO i=1,iim+1
232          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
233        ENDDO
234      ENDDO
235    ENDDO
236!$OMP END DO NOWAIT
237 
238  END SUBROUTINE filtre_v_fft
239
240  SUBROUTINE filtre_inv_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
241    USE lmdz_fft
242#ifdef CPP_PARA
243    USE parallel_lmdz,ONLY: OMP_CHUNK
244#endif
245    IMPLICIT NONE
246    INCLUDE 'dimensions.h'
247    INTEGER,INTENT(IN) :: nlat
248    INTEGER,INTENT(IN) :: jj_begin
249    INTEGER,INTENT(IN) :: jj_end
250    INTEGER,INTENT(IN) :: nbniv
251    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
252
253     REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
254    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
255    INTEGER            :: nb_vect
256    INTEGER :: i,j,l
257    INTEGER :: ll_nb
258   
259    ll_nb=0
260!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
261    DO l=1,nbniv
262      ll_nb=ll_nb+1
263      DO j=1,jj_end-jj_begin+1
264        DO i=1,iim+1
265          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
266        ENDDO
267      ENDDO
268    ENDDO
269!$OMP END DO NOWAIT
270
271    nb_vect=(jj_end-jj_begin+1)*ll_nb
272
273    CALL FFT_forward(vect,TF_vect,nb_vect)
274 
275    DO l=1,ll_nb
276      DO j=1,jj_end-jj_begin+1
277        DO i=1,iim/2+1
278          TF_vect(i,j,l)=TF_vect(i,j,l)*filtre_inv(i,jj_begin+j-1)
279        ENDDO
280      ENDDO
281    ENDDO
282 
283    CALL FFT_backward(TF_vect,vect,nb_vect)
284
285    ll_nb=0
286!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
287    DO l=1,nbniv
288      ll_nb=ll_nb+1
289      DO j=1,jj_end-jj_begin+1
290        DO i=1,iim+1
291          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
292        ENDDO
293      ENDDO
294    ENDDO
295!$OMP END DO NOWAIT
296
297  END SUBROUTINE filtre_inv_fft 
298   
299END MODULE lmdz_filtre_fft
300 
Note: See TracBrowser for help on using the repository browser.