Context Navigation

filtreg_p.F @ 4084

Last change on this file since 4084 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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: 12.9 KB

Line
1
2
3	SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv,
4	& ifiltre, iaire, griscal ,iter)
5	USE parallel_lmdz, only : OMP_CHUNK
6	USE mod_filtre_fft
7	USE timer_filtre
8
9	USE filtreg_mod
10
11	IMPLICIT NONE
12
13	c=======================================================================
14	c
15	c Auteur: P. Le Van 07/10/97
16	c ------
17	c
18	c Objet: filtre matriciel longitudinal ,avec les matrices precalculees
19	c pour l'operateur Filtre .
20	c ------
21	c
22	c Arguments:
23	c ----------
24	c
25	c
26	c ibeg..iend lattitude a filtrer
27	c nlat nombre de latitudes du champ
28	c nbniv nombre de niveaux verticaux a filtrer
29	c champ(iip1,nblat,nbniv) en entree : champ a filtrer
30	c en sortie : champ filtre
31	c ifiltre +1 Transformee directe
32	c -1 Transformee inverse
33	c +2 Filtre directe
34	c -2 Filtre inverse
35	c
36	c iaire 1 si champ intensif
37	c 2 si champ extensif (pondere par les aires)
38	c
39	c iter 1 filtre simple
40	c
41	c=======================================================================
42	c
43	c
44	c Variable Intensive
45	c ifiltre = 1 filtre directe
46	c ifiltre =-1 filtre inverse
47	c
48	c Variable Extensive
49	c ifiltre = 2 filtre directe
50	c ifiltre =-2 filtre inverse
51	c
52	c
53	#include "dimensions.h"
54	#include "paramet.h"
55	#include "coefils.h"
56	c
57	INTEGER ibeg,iend,nlat,nbniv,ifiltre,iter
58	INTEGER i,j,l,k
59	INTEGER iim2,immjm
60	INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
61
62	REAL champ( iip1,nlat,nbniv)
63
64	LOGICAL griscal
65	INTEGER hemisph, iaire
66
67	REAL :: champ_fft(iip1,nlat,nbniv)
68	REAL :: champ_in(iip1,nlat,nbniv)
69
70	LOGICAL,SAVE :: first=.TRUE.
71	c$OMP THREADPRIVATE(first)
72
73	REAL, DIMENSION(iip1,nlat,nbniv) :: champ_loc
74	INTEGER :: ll_nb, nbniv_loc
75	REAL, SAVE :: sdd12(iim,4)
76	c$OMP THREADPRIVATE(sdd12)
77
78	INTEGER, PARAMETER :: type_sddu=1
79	INTEGER, PARAMETER :: type_sddv=2
80	INTEGER, PARAMETER :: type_unsddu=3
81	INTEGER, PARAMETER :: type_unsddv=4
82
83	INTEGER :: sdd1_type, sdd2_type
84
85	IF (first) THEN
86	sdd12(1:iim,type_sddu) = sddu(1:iim)
87	sdd12(1:iim,type_sddv) = sddv(1:iim)
88	sdd12(1:iim,type_unsddu) = unsddu(1:iim)
89	sdd12(1:iim,type_unsddv) = unsddv(1:iim)
90
91	CALL Init_timer
92	first=.FALSE.
93	ENDIF
94
95	c$OMP MASTER
96	CALL start_timer
97	c$OMP END MASTER
98
99	c-------------------------------------------------------c
100
101	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
102	& STOP'Pas de transformee simple dans cette version'
103
104	IF( iter.EQ. 2 ) THEN
105	PRINT *,' Pas d iteration du filtre dans cette version !'
106	& , ' Utiliser old_filtreg et repasser !'
107	STOP
108	ENDIF
109
110	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
111	PRINT *,' Cette routine ne calcule le filtre inverse que '
112	& , ' sur la grille des scalaires !'
113	STOP
114	ENDIF
115
116	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
117	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
118	& , ' corriger et repasser !'
119	STOP
120	ENDIF
121	c
122
123	iim2 = iim * iim
124	immjm = iim * jjm
125	c
126	c
127	IF( griscal ) THEN
128	IF( nlat. NE. jjp1 ) THEN
129	PRINT 1111
130	STOP
131	ELSE
132	c
133	IF( iaire.EQ.1 ) THEN
134	sdd1_type = type_sddv
135	sdd2_type = type_unsddv
136	ELSE
137	sdd1_type = type_unsddv
138	sdd2_type = type_sddv
139	ENDIF
140	c
141	jdfil1 = 2
142	jffil1 = jfiltnu
143	jdfil2 = jfiltsu
144	jffil2 = jjm
145	ENDIF
146	ELSE
147	IF( nlat.NE.jjm ) THEN
148	PRINT 2222
149	STOP
150	ELSE
151	c
152	IF( iaire.EQ.1 ) THEN
153	sdd1_type = type_sddu
154	sdd2_type = type_unsddu
155	ELSE
156	sdd1_type = type_unsddu
157	sdd2_type = type_sddu
158	ENDIF
159	c
160	jdfil1 = 1
161	jffil1 = jfiltnv
162	jdfil2 = jfiltsv
163	jffil2 = jjm
164	ENDIF
165	ENDIF
166	c
167	DO hemisph = 1, 2
168	c
169	IF ( hemisph.EQ.1 ) THEN
170	cym
171	jdfil = max(jdfil1,ibeg)
172	jffil = min(jffil1,iend)
173	ELSE
174	cym
175	jdfil = max(jdfil2,ibeg)
176	jffil = min(jffil2,iend)
177	ENDIF
178
179
180	cccccccccccccccccccccccccccccccccccccccccccc
181	c Utilisation du filtre classique
182	cccccccccccccccccccccccccccccccccccccccccccc
183
184	IF (.NOT. use_filtre_fft) THEN
185
186	c !---------------------------------!
187	c ! Agregation des niveau verticaux !
188	c ! uniquement necessaire pour une !
189	c ! execution OpenMP !
190	c !---------------------------------!
191	ll_nb = 0
192	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
193	DO l = 1, nbniv
194	ll_nb = ll_nb+1
195	DO j = jdfil,jffil
196	DO i = 1, iim
197	champ_loc(i,j,ll_nb) =
198	& champ(i,j,l) * sdd12(i,sdd1_type)
199	ENDDO
200	ENDDO
201	ENDDO
202	c$OMP END DO NOWAIT
203
204	nbniv_loc = ll_nb
205
206	IF( hemisph.EQ.1 ) THEN
207
208	IF( ifiltre.EQ.-2 ) THEN
209	DO j = jdfil,jffil
210	#ifdef BLAS
211	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
212	& matrinvn(1,1,j), iim,
213	& champ_loc(1,j,1), iip1*nlat, 0.0,
214	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
215	#else
216	champ_fft(:iim,j-jdfil+1,:)
217	& =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
218	#endif
219	ENDDO
220
221	ELSE IF ( griscal ) THEN
222	DO j = jdfil,jffil
223	#ifdef BLAS
224	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
225	& matriceun(1,1,j), iim,
226	& champ_loc(1,j,1), iip1*nlat, 0.0,
227	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
228	#else
229	champ_fft(:iim,j-jdfil+1,:)
230	& =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
231	#endif
232	ENDDO
233
234	ELSE
235	DO j = jdfil,jffil
236	#ifdef BLAS
237	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
238	& matricevn(1,1,j), iim,
239	& champ_loc(1,j,1), iip1*nlat, 0.0,
240	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
241	#else
242	champ_fft(:iim,j-jdfil+1,:)
243	& =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
244	#endif
245	ENDDO
246
247	ENDIF
248
249	ELSE
250
251	IF( ifiltre.EQ.-2 ) THEN
252	DO j = jdfil,jffil
253	#ifdef BLAS
254	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
255	& matrinvs(1,1,j-jfiltsu+1), iim,
256	& champ_loc(1,j,1), iip1*nlat, 0.0,
257	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
258	#else
259	champ_fft(:iim,j-jdfil+1,:)
260	& =matmul(matrinvs(:,:,j-jfiltsu+1),
261	& champ_loc(:iim,j,:))
262	#endif
263	ENDDO
264
265	ELSE IF ( griscal ) THEN
266
267	DO j = jdfil,jffil
268	#ifdef BLAS
269	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
270	& matriceus(1,1,j-jfiltsu+1), iim,
271	& champ_loc(1,j,1), iip1*nlat, 0.0,
272	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
273	#else
274	champ_fft(:iim,j-jdfil+1,:)
275	& =matmul(matriceus(:,:,j-jfiltsu+1),
276	& champ_loc(:iim,j,:))
277	#endif
278	ENDDO
279
280	ELSE
281
282	DO j = jdfil,jffil
283	#ifdef BLAS
284	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
285	& matricevs(1,1,j-jfiltsv+1), iim,
286	& champ_loc(1,j,1), iip1*nlat, 0.0,
287	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
288	#else
289	champ_fft(:iim,j-jdfil+1,:)
290	& =matmul(matricevs(:,:,j-jfiltsv+1),
291	& champ_loc(:iim,j,:))
292	#endif
293	ENDDO
294
295	ENDIF
296
297	ENDIF
298	! c
299	IF( ifiltre.EQ.2 ) THEN
300
301	c !-------------------------------------!
302	c ! Dés-agregation des niveau verticaux !
303	c ! uniquement necessaire pour une !
304	c ! execution OpenMP !
305	c !-------------------------------------!
306	ll_nb = 0
307	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
308	DO l = 1, nbniv
309	ll_nb = ll_nb + 1
310	DO j = jdfil,jffil
311	DO i = 1, iim
312	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
313	& + champ_fft(i,j-jdfil+1,ll_nb))
314	& * sdd12(i,sdd2_type)
315	ENDDO
316	ENDDO
317	ENDDO
318	c$OMP END DO NOWAIT
319
320	ELSE
321
322	ll_nb = 0
323	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
324	DO l = 1, nbniv_loc
325	ll_nb = ll_nb + 1
326	DO j = jdfil,jffil
327	DO i = 1, iim
328	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
329	& - champ_fft(i,j-jdfil+1,ll_nb))
330	& * sdd12(i,sdd2_type)
331	ENDDO
332	ENDDO
333	ENDDO
334	c$OMP END DO NOWAIT
335
336	ENDIF
337
338	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
339	DO l = 1, nbniv
340	DO j = jdfil,jffil
341	champ( iip1,j,l ) = champ( 1,j,l )
342	ENDDO
343	ENDDO
344	c$OMP END DO NOWAIT
345
346	ccccccccccccccccccccccccccccccccccccccccccccc
347	c Utilisation du filtre FFT
348	ccccccccccccccccccccccccccccccccccccccccccccc
349
350	ELSE
351
352	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
353	DO l=1,nbniv
354	DO j=jdfil,jffil
355	DO i = 1, iim
356	champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
357	champ_fft( i,j,l) = champ(i,j,l)
358	ENDDO
359	ENDDO
360	ENDDO
361	c$OMP END DO NOWAIT
362
363	IF (jdfil<=jffil) THEN
364	IF( ifiltre. EQ. -2 ) THEN
365	CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
366	ELSE IF ( griscal ) THEN
367	CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
368	ELSE
369	CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
370	ENDIF
371	ENDIF
372
373
374	IF( ifiltre.EQ. 2 ) THEN
375	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
376	DO l=1,nbniv
377	DO j=jdfil,jffil
378	DO i = 1, iim
379	champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
380	& *sdd12(i,sdd2_type)
381	ENDDO
382	ENDDO
383	ENDDO
384	c$OMP END DO NOWAIT
385	ELSE
386
387	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
388	DO l=1,nbniv
389	DO j=jdfil,jffil
390	DO i = 1, iim
391	champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
392	& *sdd12(i,sdd2_type)
393	ENDDO
394	ENDDO
395	ENDDO
396	c$OMP END DO NOWAIT
397	ENDIF
398	c
399	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
400	DO l=1,nbniv
401	DO j=jdfil,jffil
402	! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
403	champ( iip1,j,l ) = champ( 1,j,l )
404	ENDDO
405	ENDDO
406	c$OMP END DO NOWAIT
407	ENDIF
408	c Fin de la zone de filtrage
409
410
411	ENDDO
412
413	! DO j=1,nlat
414	!
415	! PRINT *,"check FFT ----> Delta(",j,")=",
416	! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
417	! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
418	! ENDDO
419
420	! PRINT *,"check FFT ----> Delta(",j,")=",
421	! & sum(champ-champ_fft)/sum(champ)
422	!
423
424	c
425	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
426	& filtrer, sur la grille des scalaires'/)
427	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
428	& ltrer, sur la grille de V ou de Z'/)
429	c$OMP MASTER
430	CALL stop_timer
431	c$OMP END MASTER
432	RETURN
433	END

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

Download in other formats:

Original Format