Context Navigation

lmdz_filtreg_p.F90 @ 5209

Last change on this file since 5209 was 5159, checked in by abarral, 3 months ago
Put dimensions.h and paramet.h into modules
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.8 KB

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

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

Download in other formats:

Original Format