Context Navigation

lmdz_filtreg_p.F90 @ 5116

Last change on this file since 5116 was 5113, checked in by abarral, 2 months ago
Rename modules in misc from _mod > lmdz_ Put cbrt.f90, ch.f90, pch.f90 in new lmdz_libmath_pch.f90
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: 12.9 KB

Line
1	MODULE lmdz_filtreg_p
2	USE lmdz_filtreg, ONLY: matrinvn, matrinvs, matriceun, matriceus, matricevn, matricevs
3
4	IMPLICIT NONE; PRIVATE
5	PUBLIC filtreg_p
6
7	CONTAINS
8
9	SUBROUTINE filtreg_p(champ, jjb, jje, ibeg, iend, nlat, nbniv, &
10	ifiltre, iaire, griscal, iter)
11	USE parallel_lmdz, ONLY: OMP_CHUNK
12	USE mod_filtre_fft_loc, ONLY: use_filtre_fft, filtre_u_fft, &
13	filtre_v_fft, filtre_inv_fft
14	USE lmdz_timer_filtre, ONLY: init_timer, start_timer, stop_timer
15	USE lmdz_coefils, ONLY: jfiltnu, jfiltnv, jfiltsu, jfiltsv, sddu, sddv, unsddu, unsddv, modfrstv, modfrstu
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	!
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==1.or.ifiltre==-1) &
105	CALL abort_gcm("lmdz_filtreg_p", 'Pas de transformee&
106	&simple dans cette version', 1)
107
108	IF(iter== 2) THEN
109	PRINT *, ' Pas d iteration du filtre dans cette version !'&
110	&, ' Utiliser old_filtreg et repasser !'
111	CALL abort_gcm("lmdz_filtreg_p", "stopped", 1)
112	ENDIF
113
114	IF(ifiltre== -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("lmdz_filtreg_p", "stopped", 1)
118	ENDIF
119
120	IF(ifiltre/=2 .AND.ifiltre/= - 2) THEN
121	PRINT *, ' Probleme dans filtreg car ifiltre NE 2 et NE -2' &
122	, ' corriger et repasser !'
123	CALL abort_gcm("lmdz_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 /= jjp1) THEN
133	CALL abort_gcm("lmdz_filtreg_p", " nlat. NE. jjp1", 1)
134	ELSE
135	!
136	IF(iaire==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/=jjm) THEN
151	CALL abort_gcm("lmdz_filtreg_p", " nlat. NE. jjm", 1)
152	ELSE
153	!
154	IF(iaire==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==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==1) THEN
209
210	IF(ifiltre==-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==-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==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 == -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	IF(ifiltre== 2) THEN
380	!$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	!$OMP END DO NOWAIT
390	ELSE
391
392	!$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	!$OMP END DO NOWAIT
402	ENDIF
403	!
404	!$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	! ! add redundant longitude
409	champ(iip1, j, l) = champ(1, j, l)
410	ENDDO
411	ENDDO
412	!$OMP END DO NOWAIT
413	ENDIF
414	! Fin de la zone de filtrage
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	!$OMP MASTER
430	CALL stop_timer
431	!$OMP END MASTER
432
433	END SUBROUTINE filtreg_p
434	END MODULE lmdz_filtreg_p
435

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

Download in other formats:

Original Format