Context Navigation

filtreg_p.F90 @ 5115

Last change on this file since 5115 was 5113, checked in by abarral, 4 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: 11.7 KB

Line
1	! Amaury: on a ce fichier + lmdz_filtreg_p ! C'est totalement incompréhensible ;_;
2	! A minima il faut un nom clair pour chaque
3
4	SUBROUTINE filtreg_p( champ, ibeg, iend, nlat, nbniv, &
5	ifiltre, iaire, griscal ,iter)
6	USE parallel_lmdz, ONLY: OMP_CHUNK
7	USE mod_filtre_fft
8	USE lmdz_timer_filtre
9	USE lmdz_coefils, ONLY: jfiltnu, jfiltnv, jfiltsu, jfiltsv, sddu, sddv, unsddu, unsddv, modfrstv, modfrstu
10	USE lmdz_filtreg, ONLY: matrinvn, matrinvs, matriceun, matriceus, matricevn, matricevs
11
12	IMPLICIT NONE
13
14	!=======================================================================
15	!
16	! Auteur: P. Le Van 07/10/97
17	! ------
18	!
19	! Objet: filtre matriciel longitudinal ,avec les matrices precalculees
20	! pour l'operateur Filtre .
21	! ------
22	!
23	! Arguments:
24	! ----------
25	!
26	!
27	! ibeg..iend lattitude a filtrer
28	! nlat nombre de latitudes du champ
29	! nbniv nombre de niveaux verticaux a filtrer
30	! champ(iip1,nblat,nbniv) en entree : champ a filtrer
31	! en sortie : champ filtre
32	! ifiltre +1 Transformee directe
33	! -1 Transformee inverse
34	! +2 Filtre directe
35	! -2 Filtre inverse
36	!
37	! iaire 1 si champ intensif
38	! 2 si champ extensif (pondere par les aires)
39	!
40	! iter 1 filtre simple
41	!
42	!=======================================================================
43	!
44	!
45	! Variable Intensive
46	! ifiltre = 1 filtre directe
47	! ifiltre =-1 filtre inverse
48	!
49	! Variable Extensive
50	! ifiltre = 2 filtre directe
51	! ifiltre =-2 filtre inverse
52	!
53	!
54	INCLUDE "dimensions.h"
55	INCLUDE "paramet.h"
56
57	INTEGER :: ibeg,iend,nlat,nbniv,ifiltre,iter
58	INTEGER :: i,j,l,k
59	INTEGER :: iim2,immjm
60	INTEGER :: jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
61
62	REAL :: champ( iip1,nlat,nbniv)
63
64	LOGICAL :: griscal
65	INTEGER :: hemisph, iaire
66
67	REAL :: champ_fft(iip1,nlat,nbniv)
68	REAL :: champ_in(iip1,nlat,nbniv)
69
70	LOGICAL,SAVE :: first=.TRUE.
71	!$OMP THREADPRIVATE(first)
72
73	REAL, DIMENSION(iip1,nlat,nbniv) :: champ_loc
74	INTEGER :: ll_nb, nbniv_loc
75	REAL, SAVE :: sdd12(iim,4)
76	!$OMP THREADPRIVATE(sdd12)
77
78	INTEGER, PARAMETER :: type_sddu=1
79	INTEGER, PARAMETER :: type_sddv=2
80	INTEGER, PARAMETER :: type_unsddu=3
81	INTEGER, PARAMETER :: type_unsddv=4
82
83	INTEGER :: sdd1_type, sdd2_type
84
85	IF (first) THEN
86	sdd12(1:iim,type_sddu) = sddu(1:iim)
87	sdd12(1:iim,type_sddv) = sddv(1:iim)
88	sdd12(1:iim,type_unsddu) = unsddu(1:iim)
89	sdd12(1:iim,type_unsddv) = unsddv(1:iim)
90
91	CALL Init_timer
92	first=.FALSE.
93	ENDIF
94
95	!$OMP MASTER
96	CALL start_timer
97	!$OMP END MASTER
98
99	!-------------------------------------------------------c
100
101	IF(ifiltre==1.or.ifiltre==-1) &
102	CALL abort_gcm("fitreg_p","Pas de transformee simple&
103	&dans cette version",1)
104
105	IF( iter== 2 ) THEN
106	PRINT *,' Pas d iteration du filtre dans cette version !'&
107	& , ' Utiliser old_filtreg et repasser !'
108	CALL abort_gcm("fitreg_p","stopped",1)
109	ENDIF
110
111	IF( ifiltre== -2 .AND..NOT.griscal ) THEN
112	PRINT *,' Cette routine ne calcule le filtre inverse que ' &
113	, ' sur la grille des scalaires !'
114	CALL abort_gcm("fitreg_p","stopped",1)
115	ENDIF
116
117	IF( ifiltre/=2 .AND.ifiltre/= - 2 ) THEN
118	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2' &
119	, ' corriger et repasser !'
120	CALL abort_gcm("fitreg_p","stopped",1)
121	ENDIF
122	!
123
124	iim2 = iim * iim
125	immjm = iim * jjm
126	!
127	!
128	IF( griscal ) THEN
129	IF( nlat /= jjp1 ) THEN
130	CALL abort_gcm("fitreg_p","nlat. NE. jjp1",1)
131	ELSE
132	!
133	IF( iaire==1 ) THEN
134	sdd1_type = type_sddv
135	sdd2_type = type_unsddv
136	ELSE
137	sdd1_type = type_unsddv
138	sdd2_type = type_sddv
139	ENDIF
140	!
141	jdfil1 = 2
142	jffil1 = jfiltnu
143	jdfil2 = jfiltsu
144	jffil2 = jjm
145	ENDIF
146	ELSE
147	IF( nlat/=jjm ) THEN
148	CALL abort_gcm("fitreg_p","nlat. NE. jjm",1)
149	ELSE
150	!
151	IF( iaire==1 ) THEN
152	sdd1_type = type_sddu
153	sdd2_type = type_unsddu
154	ELSE
155	sdd1_type = type_unsddu
156	sdd2_type = type_sddu
157	ENDIF
158	!
159	jdfil1 = 1
160	jffil1 = jfiltnv
161	jdfil2 = jfiltsv
162	jffil2 = jjm
163	ENDIF
164	ENDIF
165	!
166	DO hemisph = 1, 2
167	!
168	IF ( hemisph==1 ) THEN
169	!ym
170	jdfil = max(jdfil1,ibeg)
171	jffil = min(jffil1,iend)
172	ELSE
173	!ym
174	jdfil = max(jdfil2,ibeg)
175	jffil = min(jffil2,iend)
176	ENDIF
177
178
179	!ccccccccccccccccccccccccccccccccccccccccccc
180	! Utilisation du filtre classique
181	!ccccccccccccccccccccccccccccccccccccccccccc
182
183	IF (.NOT. use_filtre_fft) THEN
184
185	!---------------------------------!
186	! Agregation des niveau verticaux !
187	! uniquement necessaire pour une !
188	! execution OpenMP !
189	!---------------------------------!
190	ll_nb = 0
191	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
192	DO l = 1, nbniv
193	ll_nb = ll_nb+1
194	DO j = jdfil,jffil
195	DO i = 1, iim
196	champ_loc(i,j,ll_nb) = &
197	champ(i,j,l) * sdd12(i,sdd1_type)
198	ENDDO
199	ENDDO
200	ENDDO
201	!$OMP END DO NOWAIT
202
203	nbniv_loc = ll_nb
204
205	IF( hemisph==1 ) THEN
206
207	IF( ifiltre==-2 ) THEN
208	DO j = jdfil,jffil
209	#ifdef BLAS
210	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
211	matrinvn(1,1,j), iim, &
212	champ_loc(1,j,1), iip1*nlat, 0.0, &
213	champ_fft(1,j-jdfil+1,1), iip1*nlat)
214	#else
215	champ_fft(:iim,j-jdfil+1,:) &
216	=matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
217	#endif
218	ENDDO
219
220	ELSE IF ( griscal ) THEN
221	DO j = jdfil,jffil
222	#ifdef BLAS
223	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
224	matriceun(1,1,j), iim, &
225	champ_loc(1,j,1), iip1*nlat, 0.0, &
226	champ_fft(1,j-jdfil+1,1), iip1*nlat)
227	#else
228	champ_fft(:iim,j-jdfil+1,:) &
229	=matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
230	#endif
231	ENDDO
232
233	ELSE
234	DO j = jdfil,jffil
235	#ifdef BLAS
236	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
237	matricevn(1,1,j), iim, &
238	champ_loc(1,j,1), iip1*nlat, 0.0, &
239	champ_fft(1,j-jdfil+1,1), iip1*nlat)
240	#else
241	champ_fft(:iim,j-jdfil+1,:) &
242	=matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
243	#endif
244	ENDDO
245
246	ENDIF
247
248	ELSE
249
250	IF( ifiltre==-2 ) THEN
251	DO j = jdfil,jffil
252	#ifdef BLAS
253	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
254	matrinvs(1,1,j-jfiltsu+1), iim, &
255	champ_loc(1,j,1), iip1*nlat, 0.0, &
256	champ_fft(1,j-jdfil+1,1), iip1*nlat)
257	#else
258	champ_fft(:iim,j-jdfil+1,:) &
259	=matmul(matrinvs(:,:,j-jfiltsu+1), &
260	champ_loc(:iim,j,:))
261	#endif
262	ENDDO
263
264	ELSE IF ( griscal ) THEN
265
266	DO j = jdfil,jffil
267	#ifdef BLAS
268	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
269	matriceus(1,1,j-jfiltsu+1), iim, &
270	champ_loc(1,j,1), iip1*nlat, 0.0, &
271	champ_fft(1,j-jdfil+1,1), iip1*nlat)
272	#else
273	champ_fft(:iim,j-jdfil+1,:) &
274	=matmul(matriceus(:,:,j-jfiltsu+1), &
275	champ_loc(:iim,j,:))
276	#endif
277	ENDDO
278
279	ELSE
280
281	DO j = jdfil,jffil
282	#ifdef BLAS
283	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
284	matricevs(1,1,j-jfiltsv+1), iim, &
285	champ_loc(1,j,1), iip1*nlat, 0.0, &
286	champ_fft(1,j-jdfil+1,1), iip1*nlat)
287	#else
288	champ_fft(:iim,j-jdfil+1,:) &
289	=matmul(matricevs(:,:,j-jfiltsv+1), &
290	champ_loc(:iim,j,:))
291	#endif
292	ENDDO
293
294	ENDIF
295
296	ENDIF
297	! c
298	IF( ifiltre==2 ) THEN
299
300	!-------------------------------------!
301	! Dés-agregation des niveau verticaux !
302	! uniquement necessaire pour une !
303	! execution OpenMP !
304	!-------------------------------------!
305	ll_nb = 0
306	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
307	DO l = 1, nbniv
308	ll_nb = ll_nb + 1
309	DO j = jdfil,jffil
310	DO i = 1, iim
311	champ( i,j,l ) = (champ_loc(i,j,ll_nb) &
312	+ champ_fft(i,j-jdfil+1,ll_nb)) &
313	* sdd12(i,sdd2_type)
314	ENDDO
315	ENDDO
316	ENDDO
317	!$OMP END DO NOWAIT
318
319	ELSE
320
321	ll_nb = 0
322	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
323	DO l = 1, nbniv_loc
324	ll_nb = ll_nb + 1
325	DO j = jdfil,jffil
326	DO i = 1, iim
327	champ( i,j,l ) = (champ_loc(i,j,ll_nb) &
328	- champ_fft(i,j-jdfil+1,ll_nb)) &
329	* sdd12(i,sdd2_type)
330	ENDDO
331	ENDDO
332	ENDDO
333	!$OMP END DO NOWAIT
334
335	ENDIF
336
337	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
338	DO l = 1, nbniv
339	DO j = jdfil,jffil
340	champ( iip1,j,l ) = champ( 1,j,l )
341	ENDDO
342	ENDDO
343	!$OMP END DO NOWAIT
344
345	!cccccccccccccccccccccccccccccccccccccccccccc
346	! Utilisation du filtre FFT
347	!cccccccccccccccccccccccccccccccccccccccccccc
348
349	ELSE
350
351	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
352	DO l=1,nbniv
353	DO j=jdfil,jffil
354	DO i = 1, iim
355	champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
356	champ_fft( i,j,l) = champ(i,j,l)
357	ENDDO
358	ENDDO
359	ENDDO
360	!$OMP END DO NOWAIT
361
362	IF (jdfil<=jffil) THEN
363	IF( ifiltre == -2 ) THEN
364	CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
365	ELSE IF ( griscal ) THEN
366	CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
367	ELSE
368	CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
369	ENDIF
370	ENDIF
371
372
373	IF( ifiltre== 2 ) THEN
374	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
375	DO l=1,nbniv
376	DO j=jdfil,jffil
377	DO i = 1, iim
378	champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l)) &
379	*sdd12(i,sdd2_type)
380	ENDDO
381	ENDDO
382	ENDDO
383	!$OMP END DO NOWAIT
384	ELSE
385
386	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
387	DO l=1,nbniv
388	DO j=jdfil,jffil
389	DO i = 1, iim
390	champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l)) &
391	*sdd12(i,sdd2_type)
392	ENDDO
393	ENDDO
394	ENDDO
395	!$OMP END DO NOWAIT
396	ENDIF
397	!
398	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
399	DO l=1,nbniv
400	DO j=jdfil,jffil
401	! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
402	champ( iip1,j,l ) = champ( 1,j,l )
403	ENDDO
404	ENDDO
405	!$OMP END DO NOWAIT
406	ENDIF
407	! Fin de la zone de filtrage
408
409
410	ENDDO
411
412	! DO j=1,nlat
413
414	! PRINT *,"check FFT ----> Delta(",j,")=",
415	! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
416	! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
417	! ENDDO
418
419	! PRINT *,"check FFT ----> Delta(",j,")=",
420	! & sum(champ-champ_fft)/sum(champ)
421
422	!
423	!$OMP MASTER
424	CALL stop_timer
425	!$OMP END MASTER
426
427	END SUBROUTINE filtreg_p

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

Download in other formats:

Original Format