Context Navigation

filtreg_p.F @ 2626

Last change on this file since 2626 was 2626, checked in by emillour, 3 years ago
Common dynamics: Some rather harmless OpenMP fixes in the filtering identified by recent versions of ifort. EM
File size: 13.7 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, only: use_filtre_fft, filtre_u_fft
7	USE mod_filtre_fft, only: filtre_v_fft, filtre_inv_fft
8	USE timer_filtre, only: init_timer, start_timer, stop_timer
9	USE filtreg_mod, only: matriceun, matricevn, matrinvn, matrinvs
10	USE filtreg_mod, only: matriceus, matricevs
11
12	IMPLICIT NONE
13
14	c=======================================================================
15	c
16	c Auteur: P. Le Van 07/10/97
17	c ------
18	c
19	c Objet: filtre matriciel longitudinal ,avec les matrices precalculees
20	c pour l'operateur Filtre .
21	c ------
22	c
23	c Arguments:
24	c ----------
25	c
26	c
27	c ibeg..iend lattitude a filtrer
28	c nlat nombre de latitudes du champ
29	c nbniv nombre de niveaux verticaux a filtrer
30	c champ(iip1,nblat,nbniv) en entree : champ a filtrer
31	c en sortie : champ filtre
32	c ifiltre +1 Transformee directe
33	c -1 Transformee inverse
34	c +2 Filtre directe
35	c -2 Filtre inverse
36	c
37	c iaire 1 si champ intensif
38	c 2 si champ extensif (pondere par les aires)
39	c
40	c iter 1 filtre simple
41	c
42	c=======================================================================
43	c
44	c
45	c Variable Intensive
46	c ifiltre = 1 filtre directe
47	c ifiltre =-1 filtre inverse
48	c
49	c Variable Extensive
50	c ifiltre = 2 filtre directe
51	c ifiltre =-2 filtre inverse
52	c
53	c
54	include "dimensions.h"
55	include "paramet.h"
56	include "coefils.h"
57	c
58	REAL,INTENT(INOUT) :: champ( iip1,nlat,nbniv)
59	INTEGER,INTENT(IN) :: ibeg,iend,nlat,nbniv,ifiltre,iter
60	LOGICAL,INTENT(IN) :: griscal
61	INTEGER,INTENT(IN) :: iaire
62
63	INTEGER i,j,l,k
64	INTEGER iim2,immjm
65	INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
66
67	INTEGER hemisph
68
69	REAL :: champ_fft(iip1,nlat,nbniv)
70	REAL :: champ_in(iip1,nlat,nbniv)
71
72	LOGICAL,SAVE :: first=.TRUE.
73	c$OMP THREADPRIVATE(first)
74
75	REAL, DIMENSION(iip1,nlat,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	nbniv_loc = ll_nb
206
207	IF( hemisph.EQ.1 ) THEN
208
209	IF( ifiltre.EQ.-2 ) THEN
210	DO j = jdfil,jffil
211	#ifdef BLAS
212	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
213	& matrinvn(1,1,j), iim,
214	& champ_loc(1,j,1), iip1*nlat, 0.0,
215	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
216	#else
217	champ_fft(:iim,j-jdfil+1,:nbniv_loc)
218	& =matmul(matrinvn(:,:,j),
219	& champ_loc(:iim,j,:nbniv_loc))
220	#endif
221	ENDDO
222
223	ELSE IF ( griscal ) THEN
224	DO j = jdfil,jffil
225	#ifdef BLAS
226	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
227	& matriceun(1,1,j), iim,
228	& champ_loc(1,j,1), iip1*nlat, 0.0,
229	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
230	#else
231	champ_fft(:iim,j-jdfil+1,:nbniv_loc)
232	& =matmul(matriceun(:,:,j),
233	& champ_loc(:iim,j,:nbniv_loc))
234	#endif
235	ENDDO
236
237	ELSE
238	DO j = jdfil,jffil
239	#ifdef BLAS
240	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
241	& matricevn(1,1,j), iim,
242	& champ_loc(1,j,1), iip1*nlat, 0.0,
243	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
244	#else
245	champ_fft(:iim,j-jdfil+1,:nbniv_loc)
246	& =matmul(matricevn(:,:,j),
247	& champ_loc(:iim,j,:nbniv_loc))
248	#endif
249	ENDDO
250
251	ENDIF ! of IF( ifiltre.EQ.-2 ) ... ELSE IF ( griscal ) ... ELSE
252
253	ELSE ! of IF( hemisph.EQ.1 )
254
255	IF( ifiltre.EQ.-2 ) THEN
256	DO j = jdfil,jffil
257	#ifdef BLAS
258	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
259	& matrinvs(1,1,j-jfiltsu+1), iim,
260	& champ_loc(1,j,1), iip1*nlat, 0.0,
261	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
262	#else
263	champ_fft(:iim,j-jdfil+1,:nbniv_loc)
264	& =matmul(matrinvs(:,:,j-jfiltsu+1),
265	& champ_loc(:iim,j,:nbniv_loc))
266	#endif
267	ENDDO
268
269	ELSE IF ( griscal ) THEN
270
271	DO j = jdfil,jffil
272	#ifdef BLAS
273	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
274	& matriceus(1,1,j-jfiltsu+1), iim,
275	& champ_loc(1,j,1), iip1*nlat, 0.0,
276	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
277	#else
278	champ_fft(:iim,j-jdfil+1,:nbniv_loc)
279	& =matmul(matriceus(:,:,j-jfiltsu+1),
280	& champ_loc(:iim,j,:nbniv_loc))
281	#endif
282	ENDDO
283
284	ELSE
285
286	DO j = jdfil,jffil
287	#ifdef BLAS
288	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
289	& matricevs(1,1,j-jfiltsv+1), iim,
290	& champ_loc(1,j,1), iip1*nlat, 0.0,
291	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
292	#else
293	champ_fft(:iim,j-jdfil+1,:nbniv_loc)
294	& =matmul(matricevs(:,:,j-jfiltsv+1),
295	& champ_loc(:iim,j,:nbniv_loc))
296	#endif
297	ENDDO
298
299	ENDIF ! of IF( ifiltre.EQ.-2 ) ... ELSE IF ( griscal ) ... ELSE
300
301	ENDIF ! of IF( hemisph.EQ.1 ) ... ELSE
302	! c
303	IF( ifiltre.EQ.2 ) THEN
304
305	c !-------------------------------------!
306	c ! Dés-agregation des niveau verticaux !
307	c ! uniquement necessaire pour une !
308	c ! execution OpenMP !
309	c !-------------------------------------!
310	ll_nb = 0
311	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
312	DO l = 1, nbniv
313	ll_nb = ll_nb + 1
314	DO j = jdfil,jffil
315	DO i = 1, iim
316	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
317	& + champ_fft(i,j-jdfil+1,ll_nb))
318	& * sdd12(i,sdd2_type)
319	ENDDO
320	ENDDO
321	ENDDO
322	c$OMP END DO NOWAIT
323
324	ELSE
325
326	ll_nb = 0
327	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
328	DO l = 1, nbniv_loc
329	ll_nb = ll_nb + 1
330	DO j = jdfil,jffil
331	DO i = 1, iim
332	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
333	& - champ_fft(i,j-jdfil+1,ll_nb))
334	& * sdd12(i,sdd2_type)
335	ENDDO
336	ENDDO
337	ENDDO
338	c$OMP END DO NOWAIT
339
340	ENDIF ! of IF( ifiltre.EQ.2 ) ... ELSE
341
342	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
343	DO l = 1, nbniv
344	DO j = jdfil,jffil
345	champ( iip1,j,l ) = champ( 1,j,l )
346	ENDDO
347	ENDDO
348	c$OMP END DO NOWAIT
349
350	ccccccccccccccccccccccccccccccccccccccccccccc
351	c Utilisation du filtre FFT
352	ccccccccccccccccccccccccccccccccccccccccccccc
353
354	ELSE ! of IF (.NOT. use_filtre_fft)
355
356	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
357	DO l=1,nbniv
358	DO j=jdfil,jffil
359	DO i = 1, iim
360	champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
361	champ_fft( i,j,l) = champ(i,j,l)
362	ENDDO
363	ENDDO
364	ENDDO
365	c$OMP END DO NOWAIT
366
367	IF (jdfil<=jffil) THEN
368	IF( ifiltre. EQ. -2 ) THEN
369	CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
370	ELSE IF ( griscal ) THEN
371	CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
372	ELSE
373	CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
374	ENDIF
375	ENDIF
376
377
378	IF( ifiltre.EQ. 2 ) THEN
379	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
380	DO l=1,nbniv
381	DO j=jdfil,jffil
382	DO i = 1, iim
383	champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
384	& *sdd12(i,sdd2_type)
385	ENDDO
386	ENDDO
387	ENDDO
388	c$OMP END DO NOWAIT
389	ELSE
390
391	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
392	DO l=1,nbniv
393	DO j=jdfil,jffil
394	DO i = 1, iim
395	champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
396	& *sdd12(i,sdd2_type)
397	ENDDO
398	ENDDO
399	ENDDO
400	c$OMP END DO NOWAIT
401	ENDIF
402	c
403	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
404	DO l=1,nbniv
405	DO j=jdfil,jffil
406	! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
407	champ( iip1,j,l ) = champ( 1,j,l )
408	ENDDO
409	ENDDO
410	c$OMP END DO NOWAIT
411	ENDIF ! of IF (.NOT. use_filtre_fft) ... ELSE
412	c Fin de la zone de filtrage
413
414
415	ENDDO ! of DO hemisph = 1, 2
416
417	! DO j=1,nlat
418	!
419	! PRINT *,"check FFT ----> Delta(",j,")=",
420	! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
421	! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
422	! ENDDO
423
424	! PRINT *,"check FFT ----> Delta(",j,")=",
425	! & sum(champ-champ_fft)/sum(champ)
426	!
427
428	c
429	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
430	& filtrer, sur la grille des scalaires'/)
431	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
432	& ltrer, sur la grille de V ou de Z'/)
433	c$OMP MASTER
434	CALL stop_timer
435	c$OMP END MASTER
436
437	END

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

Download in other formats:

Original Format