Context Navigation

mod_filtreg_p.F90 @ 5280

Last change on this file since 5280 was 5272, checked in by abarral, 7 weeks ago
Turn paramet.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.3 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, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
17	ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm
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
61
62	INCLUDE "coefils.h"
63	!
64	INTEGER,INTENT(IN) :: jjb,jje,ibeg,iend,nlat,nbniv,ifiltre,iter
65	INTEGER,INTENT(IN) :: iaire
66	LOGICAL,INTENT(IN) :: griscal
67	REAL,INTENT(INOUT) :: champ( iip1,jjb:jje,nbniv)
68
69	INTEGER :: i,j,l,k
70	INTEGER :: iim2,immjm
71	INTEGER :: jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
72	INTEGER :: hemisph
73	REAL :: champ_fft(iip1,jjb:jje,nbniv)
74	! REAL :: champ_in(iip1,jjb:jje,nbniv)
75
76	LOGICAL,SAVE :: first=.TRUE.
77	!$OMP THREADPRIVATE(first)
78
79	REAL, DIMENSION(iip1,jjb:jje,nbniv) :: champ_loc
80	INTEGER :: ll_nb, nbniv_loc
81	REAL, SAVE :: sdd12(iim,4)
82	!$OMP THREADPRIVATE(sdd12)
83
84	INTEGER, PARAMETER :: type_sddu=1
85	INTEGER, PARAMETER :: type_sddv=2
86	INTEGER, PARAMETER :: type_unsddu=3
87	INTEGER, PARAMETER :: type_unsddv=4
88
89	INTEGER :: sdd1_type, sdd2_type
90	CHARACTER (LEN=132) :: abort_message
91
92	IF (first) THEN
93	sdd12(1:iim,type_sddu) = sddu(1:iim)
94	sdd12(1:iim,type_sddv) = sddv(1:iim)
95	sdd12(1:iim,type_unsddu) = unsddu(1:iim)
96	sdd12(1:iim,type_unsddv) = unsddv(1:iim)
97
98	CALL Init_timer
99	first=.FALSE.
100	ENDIF
101
102	!$OMP MASTER
103	CALL start_timer
104	!$OMP END MASTER
105
106	!-------------------------------------------------------c
107
108	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1) &
109	CALL abort_gcm("mod_filtreg_p",'Pas de transformee&
110	&simple dans cette version',1)
111
112	IF( iter.EQ. 2 ) THEN
113	PRINT *,' Pas d iteration du filtre dans cette version !'&
114	& , ' Utiliser old_filtreg et repasser !'
115	CALL abort_gcm("mod_filtreg_p","stopped",1)
116	ENDIF
117
118	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
119	PRINT *,' Cette routine ne calcule le filtre inverse que ' &
120	, ' sur la grille des scalaires !'
121	CALL abort_gcm("mod_filtreg_p","stopped",1)
122	ENDIF
123
124	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
125	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2' &
126	, ' corriger et repasser !'
127	CALL abort_gcm("mod_filtreg_p","stopped",1)
128	ENDIF
129	!
130
131	iim2 = iim * iim
132	immjm = iim * jjm
133	!
134	!
135	IF( griscal ) THEN
136	IF( nlat.NE. jjp1 ) THEN
137	CALL abort_gcm("mod_filtreg_p"," nlat.NE. jjp1",1)
138	ELSE
139	!
140	IF( iaire.EQ.1 ) THEN
141	sdd1_type = type_sddv
142	sdd2_type = type_unsddv
143	ELSE
144	sdd1_type = type_unsddv
145	sdd2_type = type_sddv
146	ENDIF
147	!
148	jdfil1 = 2
149	jffil1 = jfiltnu
150	jdfil2 = jfiltsu
151	jffil2 = jjm
152	ENDIF
153	ELSE
154	IF( nlat.NE.jjm ) THEN
155	CALL abort_gcm("mod_filtreg_p"," nlat.NE. jjm",1)
156	ELSE
157	!
158	IF( iaire.EQ.1 ) THEN
159	sdd1_type = type_sddu
160	sdd2_type = type_unsddu
161	ELSE
162	sdd1_type = type_unsddu
163	sdd2_type = type_sddu
164	ENDIF
165	!
166	jdfil1 = 1
167	jffil1 = jfiltnv
168	jdfil2 = jfiltsv
169	jffil2 = jjm
170	ENDIF
171	ENDIF
172	!
173	DO hemisph = 1, 2
174	!
175	IF ( hemisph.EQ.1 ) THEN
176	!ym
177	jdfil = max(jdfil1,ibeg)
178	jffil = min(jffil1,iend)
179	ELSE
180	!ym
181	jdfil = max(jdfil2,ibeg)
182	jffil = min(jffil2,iend)
183	ENDIF
184
185
186	!ccccccccccccccccccccccccccccccccccccccccccc
187	! Utilisation du filtre classique
188	!ccccccccccccccccccccccccccccccccccccccccccc
189
190	IF (.NOT. use_filtre_fft) THEN
191
192	! !---------------------------------!
193	! ! Agregation des niveau verticaux !
194	! ! uniquement necessaire pour une !
195	! ! execution OpenMP !
196	! !---------------------------------!
197	ll_nb = 0
198	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
199	DO l = 1, nbniv
200	ll_nb = ll_nb+1
201	DO j = jdfil,jffil
202	DO i = 1, iim
203	champ_loc(i,j,ll_nb) = &
204	champ(i,j,l) * sdd12(i,sdd1_type)
205	ENDDO
206	ENDDO
207	ENDDO
208	!$OMP END DO NOWAIT
209
210	nbniv_loc = ll_nb
211
212	IF( hemisph.EQ.1 ) THEN
213
214	IF( ifiltre.EQ.-2 ) THEN
215	DO j = jdfil,jffil
216	#ifdef BLAS
217	CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
218	matrinvn(1,1,j), iim, &
219	champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, &
220	champ_fft(1,j,1), iip1*(jje-jjb+1))
221	#else
222	champ_fft(1:iim,j,1:nbniv_loc)= &
223	matmul(matrinvn(1:iim,1:iim,j), &
224	champ_loc(1:iim,j,1:nbniv_loc))
225	#endif
226	ENDDO
227
228	ELSE IF ( griscal ) THEN
229	DO j = jdfil,jffil
230	#ifdef BLAS
231	CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
232	matriceun(1,1,j), iim, &
233	champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, &
234	champ_fft(1,j,1), iip1*(jje-jjb+1))
235	#else
236	champ_fft(1:iim,j,1:nbniv_loc)= &
237	matmul(matriceun(1:iim,1:iim,j), &
238	champ_loc(1:iim,j,1:nbniv_loc))
239	#endif
240	ENDDO
241
242	ELSE
243	DO j = jdfil,jffil
244	#ifdef BLAS
245	CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
246	matricevn(1,1,j), iim, &
247	champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, &
248	champ_fft(1,j,1), iip1*(jje-jjb+1))
249	#else
250	champ_fft(1:iim,j,1:nbniv_loc)= &
251	matmul(matricevn(1:iim,1:iim,j), &
252	champ_loc(1:iim,j,1:nbniv_loc))
253	#endif
254	ENDDO
255
256	ENDIF
257
258	ELSE
259
260	IF( ifiltre.EQ.-2 ) THEN
261	DO j = jdfil,jffil
262	#ifdef BLAS
263	CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
264	matrinvs(1,1,j-jfiltsu+1), iim, &
265	champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, &
266	champ_fft(1,j,1), iip1*(jje-jjb+1))
267	#else
268	champ_fft(1:iim,j,1:nbniv_loc)= &
269	matmul(matrinvs(1:iim,1:iim,j-jfiltsu+1), &
270	champ_loc(1:iim,j,1:nbniv_loc))
271	#endif
272	ENDDO
273
274	ELSE IF ( griscal ) THEN
275
276	DO j = jdfil,jffil
277	#ifdef BLAS
278	CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
279	matriceus(1,1,j-jfiltsu+1), iim, &
280	champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, &
281	champ_fft(1,j,1), iip1*(jje-jjb+1))
282	#else
283	champ_fft(1:iim,j,1:nbniv_loc)= &
284	matmul(matriceus(1:iim,1:iim,j-jfiltsu+1), &
285	champ_loc(1:iim,j,1:nbniv_loc))
286	#endif
287	ENDDO
288
289	ELSE
290
291	DO j = jdfil,jffil
292	#ifdef BLAS
293	CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
294	matricevs(1,1,j-jfiltsv+1), iim, &
295	champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, &
296	champ_fft(1,j,1), iip1*(jje-jjb+1))
297	#else
298	champ_fft(1:iim,j,1:nbniv_loc)= &
299	matmul(matricevs(1:iim,1:iim,j-jfiltsv+1), &
300	champ_loc(1:iim,j,1:nbniv_loc))
301	#endif
302	ENDDO
303
304	ENDIF
305
306	ENDIF
307	! c
308	IF( ifiltre.EQ.2 ) THEN
309
310	! !-------------------------------------!
311	! ! DÃ©s-agregation des niveau verticaux !
312	! ! uniquement necessaire pour une !
313	! ! execution OpenMP !
314	! !-------------------------------------!
315	ll_nb = 0
316	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
317	DO l = 1, nbniv
318	ll_nb = ll_nb + 1
319	DO j = jdfil,jffil
320	DO i = 1, iim
321	champ( i,j,l ) = (champ_loc(i,j,ll_nb) &
322	+ champ_fft(i,j,ll_nb)) &
323	* sdd12(i,sdd2_type)
324	ENDDO
325	ENDDO
326	ENDDO
327	!$OMP END DO NOWAIT
328
329	ELSE
330
331	ll_nb = 0
332	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
333	DO l = 1, nbniv
334	ll_nb = ll_nb + 1
335	DO j = jdfil,jffil
336	DO i = 1, iim
337	champ( i,j,l ) = (champ_loc(i,j,ll_nb) &
338	- champ_fft(i,j,ll_nb)) &
339	* sdd12(i,sdd2_type)
340	ENDDO
341	ENDDO
342	ENDDO
343	!$OMP END DO NOWAIT
344
345	ENDIF
346
347	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
348	DO l = 1, nbniv
349	DO j = jdfil,jffil
350	! ! add redundant longitude
351	champ( iip1,j,l ) = champ( 1,j,l )
352	ENDDO
353	ENDDO
354	!$OMP END DO NOWAIT
355
356	!cccccccccccccccccccccccccccccccccccccccccccc
357	! Utilisation du filtre FFT
358	!cccccccccccccccccccccccccccccccccccccccccccc
359
360	ELSE
361
362	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
363	DO l=1,nbniv
364	DO j=jdfil,jffil
365	DO i = 1, iim
366	champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
367	champ_fft( i,j,l) = champ(i,j,l)
368	ENDDO
369	ENDDO
370	ENDDO
371	!$OMP END DO NOWAIT
372
373	IF (jdfil<=jffil) THEN
374	IF( ifiltre.EQ. -2 ) THEN
375	CALL Filtre_inv_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
376	ELSE IF ( griscal ) THEN
377	CALL Filtre_u_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
378	ELSE
379	CALL Filtre_v_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv)
380	ENDIF
381	ENDIF
382
383
384	IF( ifiltre.EQ. 2 ) THEN
385	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
386	DO l=1,nbniv
387	DO j=jdfil,jffil
388	DO i = 1, iim
389	champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l)) &
390	*sdd12(i,sdd2_type)
391	ENDDO
392	ENDDO
393	ENDDO
394	!$OMP END DO NOWAIT
395	ELSE
396
397	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
398	DO l=1,nbniv
399	DO j=jdfil,jffil
400	DO i = 1, iim
401	champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l)) &
402	*sdd12(i,sdd2_type)
403	ENDDO
404	ENDDO
405	ENDDO
406	!$OMP END DO NOWAIT
407	ENDIF
408	!
409	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
410	DO l=1,nbniv
411	DO j=jdfil,jffil
412	! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
413	! ! add redundant longitude
414	champ( iip1,j,l ) = champ( 1,j,l )
415	ENDDO
416	ENDDO
417	!$OMP END DO NOWAIT
418	ENDIF
419	! Fin de la zone de filtrage
420
421
422	ENDDO
423
424	! DO j=1,nlat
425	!
426	! PRINT *,"check FFT ----> Delta(",j,")=",
427	! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
428	! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
429	! ENDDO
430
431	! PRINT *,"check FFT ----> Delta(",j,")=",
432	! & sum(champ-champ_fft)/sum(champ)
433	!
434
435	!
436	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a&
437	& filtrer, sur la grille des scalaires'/)
438	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi&
439	& ltrer, sur la grille de V ou de Z'/)
440	!$OMP MASTER
441	CALL stop_timer
442	!$OMP END MASTER
443	RETURN
444	END SUBROUTINE filtreg_p
445	END MODULE mod_filtreg_p
446

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

Download in other formats:

Original Format