Context Navigation

lmdz_filtreg_p.F90 @ 5139

Last change on this file since 5139 was 5128, checked in by abarral, 8 weeks ago
Correct bug in vlspltqs_loc.f90 from r2270 where we call SSUM with incorrect arguments. Merge the three different versions of abort_gcm into one Fix seq, para 3D compilation broken from r5107 onwards (lint) usual + Remove uneeded fixed-form continuations
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 lmdz_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	!
20	! Auteur: P. Le Van 07/10/97
21	! ------
22	!
23	! Objet: filtre matriciel longitudinal ,avec les matrices precalculees
24	! pour l'operateur Filtre .
25	! ------
26	!
27	! Arguments:
28	! ----------
29	!
30	!
31	! ibeg..iend lattitude a filtrer
32	! nlat nombre de latitudes du champ
33	! nbniv nombre de niveaux verticaux a filtrer
34	! champ(iip1,nblat,nbniv) en entree : champ a filtrer
35	! en sortie : champ filtre
36	! ifiltre +1 Transformee directe
37	! -1 Transformee inverse
38	! +2 Filtre directe
39	! -2 Filtre inverse
40	!
41	! iaire 1 si champ intensif
42	! 2 si champ extensif (pondere par les aires)
43	!
44	! iter 1 filtre simple
45	!
46	!=======================================================================
47	!
48	!
49	! Variable Intensive
50	! ifiltre = 1 filtre directe
51	! ifiltre =-1 filtre inverse
52	!
53	! Variable Extensive
54	! ifiltre = 2 filtre directe
55	! ifiltre =-2 filtre inverse
56	!
57	!
58	INCLUDE "dimensions.h"
59	INCLUDE "paramet.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==1.OR.ifiltre==-1) &
106	CALL abort_gcm("lmdz_filtreg_p", 'Pas de transformee&
107	&simple dans cette version', 1)
108
109	IF(iter== 2) THEN
110	PRINT *, ' Pas d iteration du filtre dans cette version !'&
111	&, ' Utiliser old_filtreg et repasser !'
112	CALL abort_gcm("lmdz_filtreg_p", "stopped", 1)
113	ENDIF
114
115	IF(ifiltre== -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("lmdz_filtreg_p", "stopped", 1)
119	ENDIF
120
121	IF(ifiltre/=2 .AND.ifiltre/= - 2) THEN
122	PRINT *, ' Probleme dans filtreg car ifiltre NE 2 et NE -2' &
123	, ' corriger et repasser !'
124	CALL abort_gcm("lmdz_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 /= jjp1) THEN
134	CALL abort_gcm("lmdz_filtreg_p", " nlat. NE. jjp1", 1)
135	ELSE
136	!
137	IF(iaire==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/=jjm) THEN
152	CALL abort_gcm("lmdz_filtreg_p", " nlat. NE. jjm", 1)
153	ELSE
154	!
155	IF(iaire==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==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==1) THEN
210
211	IF(ifiltre==-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==-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==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 == -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	IF(ifiltre== 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	ENDDO
418
419	! DO j=1,nlat
420
421	! PRINT *,"check FFT ----> Delta(",j,")=",
422	! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
423	! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
424	! ENDDO
425
426	! PRINT *,"check FFT ----> Delta(",j,")=",
427	! & sum(champ-champ_fft)/sum(champ)
428
429	!
430	!$OMP MASTER
431	CALL stop_timer
432	!$OMP END MASTER
433
434	END SUBROUTINE filtreg_p
435	END MODULE lmdz_filtreg_p
436

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

Download in other formats:

Original Format