Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

mod_filtreg_p.F90 @ 5248

Last change on this file since 5248 was 5246, checked in by abarral, 21 hours ago
Convert fixed-form to free-form sources .F -> .{f,F}90 (WIP: some .F remain, will be handled in subsequent commits)
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, 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	!=======================================================================
18	!
19	! Auteur: P. Le Van 07/10/97
20	! ------
21	!
22	! Objet: filtre matriciel longitudinal ,avec les matrices precalculees
23	! pour l'operateur Filtre .
24	! ------
25	!
26	! Arguments:
27	! ----------
28	!
29	!
30	! ibeg..iend lattitude a filtrer
31	! nlat nombre de latitudes du champ
32	! nbniv nombre de niveaux verticaux a filtrer
33	! champ(iip1,nblat,nbniv) en entree : champ a filtrer
34	! en sortie : champ filtre
35	! ifiltre +1 Transformee directe
36	! -1 Transformee inverse
37	! +2 Filtre directe
38	! -2 Filtre inverse
39	!
40	! iaire 1 si champ intensif
41	! 2 si champ extensif (pondere par les aires)
42	!
43	! iter 1 filtre simple
44	!
45	!=======================================================================
46	!
47	!
48	! Variable Intensive
49	! ifiltre = 1 filtre directe
50	! ifiltre =-1 filtre inverse
51	!
52	! Variable Extensive
53	! ifiltre = 2 filtre directe
54	! ifiltre =-2 filtre inverse
55	!
56	!
57	INCLUDE "dimensions.h"
58	INCLUDE "paramet.h"
59	INCLUDE "coefils.h"
60	!
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	!$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	!$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	CHARACTER (LEN=132) :: abort_message
88
89	IF (first) THEN
90	sdd12(1:iim,type_sddu) = sddu(1:iim)
91	sdd12(1:iim,type_sddv) = sddv(1:iim)
92	sdd12(1:iim,type_unsddu) = unsddu(1:iim)
93	sdd12(1:iim,type_unsddv) = unsddv(1:iim)
94
95	CALL Init_timer
96	first=.FALSE.
97	ENDIF
98
99	!$OMP MASTER
100	CALL start_timer
101	!$OMP END MASTER
102
103	!-------------------------------------------------------c
104
105	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1) &
106	CALL abort_gcm("mod_filtreg_p",'Pas de transformee&
107	&simple dans cette version',1)
108
109	IF( iter.EQ. 2 ) THEN
110	PRINT *,' Pas d iteration du filtre dans cette version !'&
111	& , ' Utiliser old_filtreg et repasser !'
112	CALL abort_gcm("mod_filtreg_p","stopped",1)
113	ENDIF
114
115	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
116	PRINT *,' Cette routine ne calcule le filtre inverse que ' &
117	, ' sur la grille des scalaires !'
118	CALL abort_gcm("mod_filtreg_p","stopped",1)
119	ENDIF
120
121	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
122	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2' &
123	, ' corriger et repasser !'
124	CALL abort_gcm("mod_filtreg_p","stopped",1)
125	ENDIF
126	!
127
128	iim2 = iim * iim
129	immjm = iim * jjm
130	!
131	!
132	IF( griscal ) THEN
133	IF( nlat.NE. jjp1 ) THEN
134	CALL abort_gcm("mod_filtreg_p"," nlat.NE. jjp1",1)
135	ELSE
136	!
137	IF( iaire.EQ.1 ) THEN
138	sdd1_type = type_sddv
139	sdd2_type = type_unsddv
140	ELSE
141	sdd1_type = type_unsddv
142	sdd2_type = type_sddv
143	ENDIF
144	!
145	jdfil1 = 2
146	jffil1 = jfiltnu
147	jdfil2 = jfiltsu
148	jffil2 = jjm
149	ENDIF
150	ELSE
151	IF( nlat.NE.jjm ) THEN
152	CALL abort_gcm("mod_filtreg_p"," nlat.NE. jjm",1)
153	ELSE
154	!
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	!
163	jdfil1 = 1
164	jffil1 = jfiltnv
165	jdfil2 = jfiltsv
166	jffil2 = jjm
167	ENDIF
168	ENDIF
169	!
170	DO hemisph = 1, 2
171	!
172	IF ( hemisph.EQ.1 ) THEN
173	!ym
174	jdfil = max(jdfil1,ibeg)
175	jffil = min(jffil1,iend)
176	ELSE
177	!ym
178	jdfil = max(jdfil2,ibeg)
179	jffil = min(jffil2,iend)
180	ENDIF
181
182
183	!ccccccccccccccccccccccccccccccccccccccccccc
184	! Utilisation du filtre classique
185	!ccccccccccccccccccccccccccccccccccccccccccc
186
187	IF (.NOT. use_filtre_fft) THEN
188
189	! !---------------------------------!
190	! ! Agregation des niveau verticaux !
191	! ! uniquement necessaire pour une !
192	! ! execution OpenMP !
193	! !---------------------------------!
194	ll_nb = 0
195	!$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	!$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	! !-------------------------------------!
308	! ! DÃ©s-agregation des niveau verticaux !
309	! ! uniquement necessaire pour une !
310	! ! execution OpenMP !
311	! !-------------------------------------!
312	ll_nb = 0
313	!$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	!$OMP END DO NOWAIT
325
326	ELSE
327
328	ll_nb = 0
329	!$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	!$OMP END DO NOWAIT
341
342	ENDIF
343
344	!$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	!$OMP END DO NOWAIT
352
353	!cccccccccccccccccccccccccccccccccccccccccccc
354	! Utilisation du filtre FFT
355	!cccccccccccccccccccccccccccccccccccccccccccc
356
357	ELSE
358
359	!$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	!$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	!$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	!$OMP END DO NOWAIT
392	ELSE
393
394	!$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	!$OMP END DO NOWAIT
404	ENDIF
405	!
406	!$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	!$OMP END DO NOWAIT
415	ENDIF
416	! 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	!
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	!$OMP MASTER
438	CALL stop_timer
439	!$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