Context Navigation

filtreg_p.F @ 4553

Last change on this file since 4553 was 4482, checked in by lguez, 20 months ago
Sync latest trunk changes to branch LMDZ_ECRad
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.0 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	& CALL abort_gcm("fitreg_p","Pas de transformee simple
103	&dans cette version",1)
104
105	IF( iter.EQ. 2 ) THEN
106	PRINT *,' Pas d iteration du filtre dans cette version !'
107	& , ' Utiliser old_filtreg et repasser !'
108	CALL abort_gcm("fitreg_p","stopped",1)
109	ENDIF
110
111	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
112	PRINT *,' Cette routine ne calcule le filtre inverse que '
113	& , ' sur la grille des scalaires !'
114	CALL abort_gcm("fitreg_p","stopped",1)
115	ENDIF
116
117	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
118	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
119	& , ' corriger et repasser !'
120	CALL abort_gcm("fitreg_p","stopped",1)
121	ENDIF
122	c
123
124	iim2 = iim * iim
125	immjm = iim * jjm
126	c
127	c
128	IF( griscal ) THEN
129	IF( nlat. NE. jjp1 ) THEN
130	CALL abort_gcm("fitreg_p","nlat. NE. jjp1",1)
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	CALL abort_gcm("fitreg_p","nlat. NE. jjm",1)
149	ELSE
150	c
151	IF( iaire.EQ.1 ) THEN
152	sdd1_type = type_sddu
153	sdd2_type = type_unsddu
154	ELSE
155	sdd1_type = type_unsddu
156	sdd2_type = type_sddu
157	ENDIF
158	c
159	jdfil1 = 1
160	jffil1 = jfiltnv
161	jdfil2 = jfiltsv
162	jffil2 = jjm
163	ENDIF
164	ENDIF
165	c
166	DO hemisph = 1, 2
167	c
168	IF ( hemisph.EQ.1 ) THEN
169	cym
170	jdfil = max(jdfil1,ibeg)
171	jffil = min(jffil1,iend)
172	ELSE
173	cym
174	jdfil = max(jdfil2,ibeg)
175	jffil = min(jffil2,iend)
176	ENDIF
177
178
179	cccccccccccccccccccccccccccccccccccccccccccc
180	c Utilisation du filtre classique
181	cccccccccccccccccccccccccccccccccccccccccccc
182
183	IF (.NOT. use_filtre_fft) THEN
184
185	c !---------------------------------!
186	c ! Agregation des niveau verticaux !
187	c ! uniquement necessaire pour une !
188	c ! execution OpenMP !
189	c !---------------------------------!
190	ll_nb = 0
191	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
192	DO l = 1, nbniv
193	ll_nb = ll_nb+1
194	DO j = jdfil,jffil
195	DO i = 1, iim
196	champ_loc(i,j,ll_nb) =
197	& champ(i,j,l) * sdd12(i,sdd1_type)
198	ENDDO
199	ENDDO
200	ENDDO
201	c$OMP END DO NOWAIT
202
203	nbniv_loc = ll_nb
204
205	IF( hemisph.EQ.1 ) THEN
206
207	IF( ifiltre.EQ.-2 ) THEN
208	DO j = jdfil,jffil
209	#ifdef BLAS
210	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
211	& matrinvn(1,1,j), iim,
212	& champ_loc(1,j,1), iip1*nlat, 0.0,
213	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
214	#else
215	champ_fft(:iim,j-jdfil+1,:)
216	& =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
217	#endif
218	ENDDO
219
220	ELSE IF ( griscal ) THEN
221	DO j = jdfil,jffil
222	#ifdef BLAS
223	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
224	& matriceun(1,1,j), iim,
225	& champ_loc(1,j,1), iip1*nlat, 0.0,
226	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
227	#else
228	champ_fft(:iim,j-jdfil+1,:)
229	& =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
230	#endif
231	ENDDO
232
233	ELSE
234	DO j = jdfil,jffil
235	#ifdef BLAS
236	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
237	& matricevn(1,1,j), iim,
238	& champ_loc(1,j,1), iip1*nlat, 0.0,
239	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
240	#else
241	champ_fft(:iim,j-jdfil+1,:)
242	& =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
243	#endif
244	ENDDO
245
246	ENDIF
247
248	ELSE
249
250	IF( ifiltre.EQ.-2 ) THEN
251	DO j = jdfil,jffil
252	#ifdef BLAS
253	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
254	& matrinvs(1,1,j-jfiltsu+1), iim,
255	& champ_loc(1,j,1), iip1*nlat, 0.0,
256	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
257	#else
258	champ_fft(:iim,j-jdfil+1,:)
259	& =matmul(matrinvs(:,:,j-jfiltsu+1),
260	& champ_loc(:iim,j,:))
261	#endif
262	ENDDO
263
264	ELSE IF ( griscal ) THEN
265
266	DO j = jdfil,jffil
267	#ifdef BLAS
268	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
269	& matriceus(1,1,j-jfiltsu+1), iim,
270	& champ_loc(1,j,1), iip1*nlat, 0.0,
271	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
272	#else
273	champ_fft(:iim,j-jdfil+1,:)
274	& =matmul(matriceus(:,:,j-jfiltsu+1),
275	& champ_loc(:iim,j,:))
276	#endif
277	ENDDO
278
279	ELSE
280
281	DO j = jdfil,jffil
282	#ifdef BLAS
283	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
284	& matricevs(1,1,j-jfiltsv+1), iim,
285	& champ_loc(1,j,1), iip1*nlat, 0.0,
286	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
287	#else
288	champ_fft(:iim,j-jdfil+1,:)
289	& =matmul(matricevs(:,:,j-jfiltsv+1),
290	& champ_loc(:iim,j,:))
291	#endif
292	ENDDO
293
294	ENDIF
295
296	ENDIF
297	! c
298	IF( ifiltre.EQ.2 ) THEN
299
300	c !-------------------------------------!
301	c ! Dés-agregation des niveau verticaux !
302	c ! uniquement necessaire pour une !
303	c ! execution OpenMP !
304	c !-------------------------------------!
305	ll_nb = 0
306	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
307	DO l = 1, nbniv
308	ll_nb = ll_nb + 1
309	DO j = jdfil,jffil
310	DO i = 1, iim
311	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
312	& + champ_fft(i,j-jdfil+1,ll_nb))
313	& * sdd12(i,sdd2_type)
314	ENDDO
315	ENDDO
316	ENDDO
317	c$OMP END DO NOWAIT
318
319	ELSE
320
321	ll_nb = 0
322	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
323	DO l = 1, nbniv_loc
324	ll_nb = ll_nb + 1
325	DO j = jdfil,jffil
326	DO i = 1, iim
327	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
328	& - champ_fft(i,j-jdfil+1,ll_nb))
329	& * sdd12(i,sdd2_type)
330	ENDDO
331	ENDDO
332	ENDDO
333	c$OMP END DO NOWAIT
334
335	ENDIF
336
337	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
338	DO l = 1, nbniv
339	DO j = jdfil,jffil
340	champ( iip1,j,l ) = champ( 1,j,l )
341	ENDDO
342	ENDDO
343	c$OMP END DO NOWAIT
344
345	ccccccccccccccccccccccccccccccccccccccccccccc
346	c Utilisation du filtre FFT
347	ccccccccccccccccccccccccccccccccccccccccccccc
348
349	ELSE
350
351	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
352	DO l=1,nbniv
353	DO j=jdfil,jffil
354	DO i = 1, iim
355	champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
356	champ_fft( i,j,l) = champ(i,j,l)
357	ENDDO
358	ENDDO
359	ENDDO
360	c$OMP END DO NOWAIT
361
362	IF (jdfil<=jffil) THEN
363	IF( ifiltre. EQ. -2 ) THEN
364	CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
365	ELSE IF ( griscal ) THEN
366	CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
367	ELSE
368	CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
369	ENDIF
370	ENDIF
371
372
373	IF( ifiltre.EQ. 2 ) THEN
374	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
375	DO l=1,nbniv
376	DO j=jdfil,jffil
377	DO i = 1, iim
378	champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
379	& *sdd12(i,sdd2_type)
380	ENDDO
381	ENDDO
382	ENDDO
383	c$OMP END DO NOWAIT
384	ELSE
385
386	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
387	DO l=1,nbniv
388	DO j=jdfil,jffil
389	DO i = 1, iim
390	champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
391	& *sdd12(i,sdd2_type)
392	ENDDO
393	ENDDO
394	ENDDO
395	c$OMP END DO NOWAIT
396	ENDIF
397	c
398	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
399	DO l=1,nbniv
400	DO j=jdfil,jffil
401	! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
402	champ( iip1,j,l ) = champ( 1,j,l )
403	ENDDO
404	ENDDO
405	c$OMP END DO NOWAIT
406	ENDIF
407	c Fin de la zone de filtrage
408
409
410	ENDDO
411
412	! DO j=1,nlat
413	!
414	! PRINT *,"check FFT ----> Delta(",j,")=",
415	! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
416	! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
417	! ENDDO
418
419	! PRINT *,"check FFT ----> Delta(",j,")=",
420	! & sum(champ-champ_fft)/sum(champ)
421	!
422
423	c
424	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
425	& filtrer, sur la grille des scalaires'/)
426	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
427	& ltrer, sur la grille de V ou de Z'/)
428	c$OMP MASTER
429	CALL stop_timer
430	c$OMP END MASTER
431	RETURN
432	END

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

Download in other formats:

Original Format