Context Navigation

mod_filtreg_p.F90 @ 5407

Last change on this file since 5407 was 5297, checked in by abarral, 7 weeks ago
Turn gradsdef.h coefils.h into a module
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	USE dimensions_mod, ONLY: iim, jjm, llm, ndm
16	USE paramet_mod_h
17	USE coefils_mod_h
18	IMPLICIT NONE
19
20	!=======================================================================
21	!
22	! Auteur: P. Le Van 07/10/97
23	! ------
24	!
25	! Objet: filtre matriciel longitudinal ,avec les matrices precalculees
26	! pour l'operateur Filtre .
27	! ------
28	!
29	! Arguments:
30	! ----------
31	!
32	!
33	! ibeg..iend lattitude a filtrer
34	! nlat nombre de latitudes du champ
35	! nbniv nombre de niveaux verticaux a filtrer
36	! champ(iip1,nblat,nbniv) en entree : champ a filtrer
37	! en sortie : champ filtre
38	! ifiltre +1 Transformee directe
39	! -1 Transformee inverse
40	! +2 Filtre directe
41	! -2 Filtre inverse
42	!
43	! iaire 1 si champ intensif
44	! 2 si champ extensif (pondere par les aires)
45	!
46	! iter 1 filtre simple
47	!
48	!=======================================================================
49	!
50	!
51	! Variable Intensive
52	! ifiltre = 1 filtre directe
53	! ifiltre =-1 filtre inverse
54	!
55	! Variable Extensive
56	! ifiltre = 2 filtre directe
57	! ifiltre =-2 filtre inverse
58	!
59	!
60	INTEGER,INTENT(IN) :: jjb,jje,ibeg,iend,nlat,nbniv,ifiltre,iter
61	INTEGER,INTENT(IN) :: iaire
62	LOGICAL,INTENT(IN) :: griscal
63	REAL,INTENT(INOUT) :: champ( iip1,jjb:jje,nbniv)
64
65	INTEGER :: i,j,l,k
66	INTEGER :: iim2,immjm
67	INTEGER :: jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
68	INTEGER :: hemisph
69	REAL :: champ_fft(iip1,jjb:jje,nbniv)
70	! REAL :: champ_in(iip1,jjb:jje,nbniv)
71
72	LOGICAL,SAVE :: first=.TRUE.
73	!$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	!$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	CHARACTER (LEN=132) :: abort_message
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	!$OMP MASTER
99	CALL start_timer
100	!$OMP END MASTER
101
102	!-------------------------------------------------------c
103
104	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1) &
105	CALL abort_gcm("mod_filtreg_p",'Pas de transformee&
106	&simple dans cette version',1)
107
108	IF( iter.EQ. 2 ) THEN
109	PRINT *,' Pas d iteration du filtre dans cette version !'&
110	& , ' Utiliser old_filtreg et repasser !'
111	CALL abort_gcm("mod_filtreg_p","stopped",1)
112	ENDIF
113
114	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
115	PRINT *,' Cette routine ne calcule le filtre inverse que ' &
116	, ' sur la grille des scalaires !'
117	CALL abort_gcm("mod_filtreg_p","stopped",1)
118	ENDIF
119
120	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
121	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2' &
122	, ' corriger et repasser !'
123	CALL abort_gcm("mod_filtreg_p","stopped",1)
124	ENDIF
125	!
126
127	iim2 = iim * iim
128	immjm = iim * jjm
129	!
130	!
131	IF( griscal ) THEN
132	IF( nlat.NE. jjp1 ) THEN
133	CALL abort_gcm("mod_filtreg_p"," nlat.NE. jjp1",1)
134	ELSE
135	!
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	!
144	jdfil1 = 2
145	jffil1 = jfiltnu
146	jdfil2 = jfiltsu
147	jffil2 = jjm
148	ENDIF
149	ELSE
150	IF( nlat.NE.jjm ) THEN
151	CALL abort_gcm("mod_filtreg_p"," nlat.NE. jjm",1)
152	ELSE
153	!
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	!
162	jdfil1 = 1
163	jffil1 = jfiltnv
164	jdfil2 = jfiltsv
165	jffil2 = jjm
166	ENDIF
167	ENDIF
168	!
169	DO hemisph = 1, 2
170	!
171	IF ( hemisph.EQ.1 ) THEN
172	!ym
173	jdfil = max(jdfil1,ibeg)
174	jffil = min(jffil1,iend)
175	ELSE
176	!ym
177	jdfil = max(jdfil2,ibeg)
178	jffil = min(jffil2,iend)
179	ENDIF
180
181
182	!ccccccccccccccccccccccccccccccccccccccccccc
183	! Utilisation du filtre classique
184	!ccccccccccccccccccccccccccccccccccccccccccc
185
186	IF (.NOT. use_filtre_fft) THEN
187
188	! !---------------------------------!
189	! ! Agregation des niveau verticaux !
190	! ! uniquement necessaire pour une !
191	! ! execution OpenMP !
192	! !---------------------------------!
193	ll_nb = 0
194	!$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	!$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(1:iim,j,1:nbniv_loc)= &
219	matmul(matrinvn(1:iim,1:iim,j), &
220	champ_loc(1: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(1:iim,j,1:nbniv_loc)= &
233	matmul(matriceun(1:iim,1:iim,j), &
234	champ_loc(1: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(1:iim,j,1:nbniv_loc)= &
247	matmul(matricevn(1:iim,1:iim,j), &
248	champ_loc(1: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(1:iim,j,1:nbniv_loc)= &
265	matmul(matrinvs(1:iim,1:iim,j-jfiltsu+1), &
266	champ_loc(1: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(1:iim,j,1:nbniv_loc)= &
280	matmul(matriceus(1:iim,1:iim,j-jfiltsu+1), &
281	champ_loc(1: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(1:iim,j,1:nbniv_loc)= &
295	matmul(matricevs(1:iim,1:iim,j-jfiltsv+1), &
296	champ_loc(1: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	! !-------------------------------------!
307	! ! DÃ©s-agregation des niveau verticaux !
308	! ! uniquement necessaire pour une !
309	! ! execution OpenMP !
310	! !-------------------------------------!
311	ll_nb = 0
312	!$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	!$OMP END DO NOWAIT
324
325	ELSE
326
327	ll_nb = 0
328	!$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	!$OMP END DO NOWAIT
340
341	ENDIF
342
343	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
344	DO l = 1, nbniv
345	DO j = jdfil,jffil
346	! ! add redundant longitude
347	champ( iip1,j,l ) = champ( 1,j,l )
348	ENDDO
349	ENDDO
350	!$OMP END DO NOWAIT
351
352	!cccccccccccccccccccccccccccccccccccccccccccc
353	! Utilisation du filtre FFT
354	!cccccccccccccccccccccccccccccccccccccccccccc
355
356	ELSE
357
358	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
359	DO l=1,nbniv
360	DO j=jdfil,jffil
361	DO i = 1, iim
362	champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
363	champ_fft( i,j,l) = champ(i,j,l)
364	ENDDO
365	ENDDO
366	ENDDO
367	!$OMP END DO NOWAIT
368
369	IF (jdfil<=jffil) THEN
370	IF( ifiltre.EQ. -2 ) THEN
371	CALL Filtre_inv_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
372	ELSE IF ( griscal ) THEN
373	CALL Filtre_u_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
374	ELSE
375	CALL Filtre_v_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
376	ENDIF
377	ENDIF
378
379
380	IF( ifiltre.EQ. 2 ) THEN
381	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
382	DO l=1,nbniv
383	DO j=jdfil,jffil
384	DO i = 1, iim
385	champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l)) &
386	*sdd12(i,sdd2_type)
387	ENDDO
388	ENDDO
389	ENDDO
390	!$OMP END DO NOWAIT
391	ELSE
392
393	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
394	DO l=1,nbniv
395	DO j=jdfil,jffil
396	DO i = 1, iim
397	champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l)) &
398	*sdd12(i,sdd2_type)
399	ENDDO
400	ENDDO
401	ENDDO
402	!$OMP END DO NOWAIT
403	ENDIF
404	!
405	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
406	DO l=1,nbniv
407	DO j=jdfil,jffil
408	! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
409	! ! add redundant longitude
410	champ( iip1,j,l ) = champ( 1,j,l )
411	ENDDO
412	ENDDO
413	!$OMP END DO NOWAIT
414	ENDIF
415	! Fin de la zone de filtrage
416
417
418	ENDDO
419
420	! DO j=1,nlat
421	!
422	! PRINT *,"check FFT ----> Delta(",j,")=",
423	! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
424	! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
425	! ENDDO
426
427	! PRINT *,"check FFT ----> Delta(",j,")=",
428	! & sum(champ-champ_fft)/sum(champ)
429	!
430
431	!
432	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a&
433	& filtrer, sur la grille des scalaires'/)
434	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi&
435	& ltrer, sur la grille de V ou de Z'/)
436	!$OMP MASTER
437	CALL stop_timer
438	!$OMP END MASTER
439	RETURN
440	END SUBROUTINE filtreg_p
441	END MODULE mod_filtreg_p
442

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

Download in other formats:

Original Format