Context Navigation

mod_filtreg_p.F @ 1707

Last change on this file since 1707 was 1707, checked in by Laurent Fairhead, 11 years ago

Version testing basée sur la r1706

Testing release based on r1706

File size: 13.0 KB

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

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

Download in other formats:

Original Format