Context Navigation

lmdz_filtre_fft.F90 @ 5204

Last change on this file since 5204 was 5159, checked in by abarral, 3 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 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 fairhead $
3
4	MODULE 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
14	CONTAINS
15
16	SUBROUTINE init_filtre_fft(coeffu,modfrstu,jfiltnu,jfiltsu,coeffv,modfrstv,jfiltnv,jfiltsv)
17	USE lmdz_fft
18	USE 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
125	USE 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
186	USE 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
248	USE 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
303	END MODULE lmdz_filtre_fft
304

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format