Context Navigation

mod_filtreg_p.F @ 3365

Last change on this file since 3365 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: 13.2 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
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,1:nbniv_loc)=
219	& matmul(matrinvn(:,:,j),
220	& champ_loc(:iim,j,1:nbniv_loc))
221	#endif
222	ENDDO
223
224	ELSE IF ( griscal ) THEN
225	DO j = jdfil,jffil
226	#ifdef BLAS
227	CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
228	& matriceun(1,1,j), iim,
229	& champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
230	& champ_fft(1,j,1), iip1*(jje-jjb+1))
231	#else
232	champ_fft(:,j,1:nbniv_loc)=
233	& matmul(matriceun(:,:,j),
234	& champ_loc(:iim,j,1:nbniv_loc))
235	#endif
236	ENDDO
237
238	ELSE
239	DO j = jdfil,jffil
240	#ifdef BLAS
241	CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
242	& matricevn(1,1,j), iim,
243	& champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
244	& champ_fft(1,j,1), iip1*(jje-jjb+1))
245	#else
246	champ_fft(:,j,1:nbniv_loc)=
247	& matmul(matricevn(:,:,j),
248	& champ_loc(:iim,j,1:nbniv_loc))
249	#endif
250	ENDDO
251
252	ENDIF
253
254	ELSE
255
256	IF( ifiltre.EQ.-2 ) THEN
257	DO j = jdfil,jffil
258	#ifdef BLAS
259	CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
260	& matrinvs(1,1,j-jfiltsu+1), iim,
261	& champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
262	& champ_fft(1,j,1), iip1*(jje-jjb+1))
263	#else
264	champ_fft(:,j,1:nbniv_loc)=
265	& matmul(matrinvs(:,:,j-jfiltsu+1),
266	& champ_loc(:iim,j,1:nbniv_loc))
267	#endif
268	ENDDO
269
270	ELSE IF ( griscal ) THEN
271
272	DO j = jdfil,jffil
273	#ifdef BLAS
274	CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
275	& matriceus(1,1,j-jfiltsu+1), iim,
276	& champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
277	& champ_fft(1,j,1), iip1*(jje-jjb+1))
278	#else
279	champ_fft(:,j,1:nbniv_loc)=
280	& matmul(matriceus(:,:,j-jfiltsu+1),
281	& champ_loc(:iim,j,1:nbniv_loc))
282	#endif
283	ENDDO
284
285	ELSE
286
287	DO j = jdfil,jffil
288	#ifdef BLAS
289	CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
290	& matricevs(1,1,j-jfiltsv+1), iim,
291	& champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0,
292	& champ_fft(1,j,1), iip1*(jje-jjb+1))
293	#else
294	champ_fft(:,j,1:nbniv_loc)=
295	& matmul(matricevs(:,:,j-jfiltsv+1),
296	& champ_loc(:iim,j,1:nbniv_loc))
297	#endif
298	ENDDO
299
300	ENDIF
301
302	ENDIF
303	! c
304	IF( ifiltre.EQ.2 ) THEN
305
306	c !-------------------------------------!
307	c ! DÃ©s-agregation des niveau verticaux !
308	c ! uniquement necessaire pour une !
309	c ! execution OpenMP !
310	c !-------------------------------------!
311	ll_nb = 0
312	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
313	DO l = 1, nbniv
314	ll_nb = ll_nb + 1
315	DO j = jdfil,jffil
316	DO i = 1, iim
317	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
318	& + champ_fft(i,j,ll_nb))
319	& * sdd12(i,sdd2_type)
320	ENDDO
321	ENDDO
322	ENDDO
323	c$OMP END DO NOWAIT
324
325	ELSE
326
327	ll_nb = 0
328	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
329	DO l = 1, nbniv
330	ll_nb = ll_nb + 1
331	DO j = jdfil,jffil
332	DO i = 1, iim
333	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
334	& - champ_fft(i,j,ll_nb))
335	& * sdd12(i,sdd2_type)
336	ENDDO
337	ENDDO
338	ENDDO
339	c$OMP END DO NOWAIT
340
341	ENDIF
342
343	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
344	DO l = 1, nbniv
345	DO j = jdfil,jffil
346	champ( iip1,j,l ) = champ( 1,j,l )
347	ENDDO
348	ENDDO
349	c$OMP END DO NOWAIT
350
351	ccccccccccccccccccccccccccccccccccccccccccccc
352	c Utilisation du filtre FFT
353	ccccccccccccccccccccccccccccccccccccccccccccc
354
355	ELSE
356
357	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
358	DO l=1,nbniv
359	DO j=jdfil,jffil
360	DO i = 1, iim
361	champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
362	champ_fft( i,j,l) = champ(i,j,l)
363	ENDDO
364	ENDDO
365	ENDDO
366	c$OMP END DO NOWAIT
367
368	IF (jdfil<=jffil) THEN
369	IF( ifiltre. EQ. -2 ) THEN
370	CALL Filtre_inv_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
371	ELSE IF ( griscal ) THEN
372	CALL Filtre_u_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
373	ELSE
374	CALL Filtre_v_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
375	ENDIF
376	ENDIF
377
378
379	IF( ifiltre.EQ. 2 ) THEN
380	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
381	DO l=1,nbniv
382	DO j=jdfil,jffil
383	DO i = 1, iim
384	champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
385	& *sdd12(i,sdd2_type)
386	ENDDO
387	ENDDO
388	ENDDO
389	c$OMP END DO NOWAIT
390	ELSE
391
392	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
393	DO l=1,nbniv
394	DO j=jdfil,jffil
395	DO i = 1, iim
396	champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
397	& *sdd12(i,sdd2_type)
398	ENDDO
399	ENDDO
400	ENDDO
401	c$OMP END DO NOWAIT
402	ENDIF
403	c
404	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
405	DO l=1,nbniv
406	DO j=jdfil,jffil
407	! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
408	champ( iip1,j,l ) = champ( 1,j,l )
409	ENDDO
410	ENDDO
411	c$OMP END DO NOWAIT
412	ENDIF
413	c Fin de la zone de filtrage
414
415
416	ENDDO
417
418	! DO j=1,nlat
419	!
420	! PRINT *,"check FFT ----> Delta(",j,")=",
421	! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
422	! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
423	! ENDDO
424
425	! PRINT *,"check FFT ----> Delta(",j,")=",
426	! & sum(champ-champ_fft)/sum(champ)
427	!
428
429	c
430	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
431	& filtrer, sur la grille des scalaires'/)
432	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
433	& ltrer, sur la grille de V ou de Z'/)
434	c$OMP MASTER
435	CALL stop_timer
436	c$OMP END MASTER
437	RETURN
438	END SUBROUTINE filtreg_p
439	END MODULE mod_filtreg_p
440

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

Download in other formats:

Original Format