Context Navigation

mod_filtreg_p.F @ 5308

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

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

Download in other formats:

Original Format