My Project
 All Classes Files Functions Variables Macros
filtreg_p.F
Go to the documentation of this file.
1 
2 
3  SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv,
4  & ifiltre, iaire, griscal ,iter)
5  USE parallel, only : OMP_CHUNK
6  USE mod_filtre_fft
7  USE timer_filtre
8 
9  USE filtreg_mod
10 
11  IMPLICIT NONE
12 
13 c=======================================================================
14 c
15 c Auteur: P. Le Van 07/10/97
16 c ------
17 c
18 c Objet: filtre matriciel longitudinal ,avec les matrices precalculees
19 c pour l'operateur Filtre .
20 c ------
21 c
22 c Arguments:
23 c ----------
24 c
25 c
26 c ibeg..iend lattitude a filtrer
27 c nlat nombre de latitudes du champ
28 c nbniv nombre de niveaux verticaux a filtrer
29 c champ(iip1,nblat,nbniv) en entree : champ a filtrer
30 c en sortie : champ filtre
31 c ifiltre +1 Transformee directe
32 c -1 Transformee inverse
33 c +2 Filtre directe
34 c -2 Filtre inverse
35 c
36 c iaire 1 si champ intensif
37 c 2 si champ extensif (pondere par les aires)
38 c
39 c iter 1 filtre simple
40 c
41 c=======================================================================
42 c
43 c
44 c Variable Intensive
45 c ifiltre = 1 filtre directe
46 c ifiltre =-1 filtre inverse
47 c
48 c Variable Extensive
49 c ifiltre = 2 filtre directe
50 c ifiltre =-2 filtre inverse
51 c
52 c
53 #include "dimensions.h"
54 #include "paramet.h"
55 #include "coefils.h"
56 c
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 c$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 c$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 c$OMP MASTER
96  CALL start_timer
97 c$OMP END MASTER
98 
99 c-------------------------------------------------------c
100 
101  IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
102  & stop'Pas de transformee simple dans cette version'
103 
104  IF( iter.EQ. 2 ) THEN
105  print *,' Pas d iteration du filtre dans cette version !'
106  & , ' Utiliser old_filtreg et repasser !'
107  stop
108  ENDIF
109 
110  IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
111  print *,' Cette routine ne calcule le filtre inverse que '
112  & , ' sur la grille des scalaires !'
113  stop
114  ENDIF
115 
116  IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
117  print *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
118  & , ' corriger et repasser !'
119  stop
120  ENDIF
121 c
122 
123  iim2 = iim * iim
124  immjm = iim * jjm
125 c
126 c
127  IF( griscal ) THEN
128  IF( nlat. ne. jjp1 ) THEN
129  print 1111
130  stop
131  ELSE
132 c
133  IF( iaire.EQ.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 c
141  jdfil1 = 2
142  jffil1 = jfiltnu
143  jdfil2 = jfiltsu
144  jffil2 = jjm
145  ENDIF
146  ELSE
147  IF( nlat.NE.jjm ) THEN
148  print 2222
149  stop
150  ELSE
151 c
152  IF( iaire.EQ.1 ) THEN
153  sdd1_type = type_sddu
154  sdd2_type = type_unsddu
155  ELSE
156  sdd1_type = type_unsddu
157  sdd2_type = type_sddu
158  ENDIF
159 c
160  jdfil1 = 1
161  jffil1 = jfiltnv
162  jdfil2 = jfiltsv
163  jffil2 = jjm
164  ENDIF
165  ENDIF
166 c
167  DO hemisph = 1, 2
168 c
169  IF ( hemisph.EQ.1 ) THEN
170 cym
171  jdfil = max(jdfil1,ibeg)
172  jffil = min(jffil1,iend)
173  ELSE
174 cym
175  jdfil = max(jdfil2,ibeg)
176  jffil = min(jffil2,iend)
177  ENDIF
178 
179 
180 cccccccccccccccccccccccccccccccccccccccccccc
181 c Utilisation du filtre classique
182 cccccccccccccccccccccccccccccccccccccccccccc
183 
184  IF (.NOT. use_filtre_fft) THEN
185 
186 c !---------------------------------!
187 c ! Agregation des niveau verticaux !
188 c ! uniquement necessaire pour une !
189 c ! execution OpenMP !
190 c !---------------------------------!
191  ll_nb = 0
192 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
193  DO l = 1, nbniv
194  ll_nb = ll_nb+1
195  DO j = jdfil,jffil
196  DO i = 1, iim
197  champ_loc(i,j,ll_nb) =
198  & champ(i,j,l) * sdd12(i,sdd1_type)
199  ENDDO
200  ENDDO
201  ENDDO
202 c$OMP END DO NOWAIT
203 
204  nbniv_loc = ll_nb
205 
206  IF( hemisph.EQ.1 ) THEN
207 
208  IF( ifiltre.EQ.-2 ) THEN
209  DO j = jdfil,jffil
210 #ifdef BLAS
211  CALL dgemm("N", "N", iim, nbniv_loc, iim, 1.0,
212  & matrinvn(1,1,j), iim,
213  & champ_loc(1,j,1), iip1*nlat, 0.0,
214  & champ_fft(1,j-jdfil+1,1), iip1*nlat)
215 #else
216  champ_fft(:iim,j-jdfil+1,:)
217  & =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
218 #endif
219  ENDDO
220 
221  ELSE IF ( griscal ) THEN
222  DO j = jdfil,jffil
223 #ifdef BLAS
224  CALL dgemm("N", "N", iim, nbniv_loc, iim, 1.0,
225  & matriceun(1,1,j), iim,
226  & champ_loc(1,j,1), iip1*nlat, 0.0,
227  & champ_fft(1,j-jdfil+1,1), iip1*nlat)
228 #else
229  champ_fft(:iim,j-jdfil+1,:)
230  & =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
231 #endif
232  ENDDO
233 
234  ELSE
235  DO j = jdfil,jffil
236 #ifdef BLAS
237  CALL dgemm("N", "N", iim, nbniv_loc, iim, 1.0,
238  & matricevn(1,1,j), iim,
239  & champ_loc(1,j,1), iip1*nlat, 0.0,
240  & champ_fft(1,j-jdfil+1,1), iip1*nlat)
241 #else
242  champ_fft(:iim,j-jdfil+1,:)
243  & =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
244 #endif
245  ENDDO
246 
247  ENDIF
248 
249  ELSE
250 
251  IF( ifiltre.EQ.-2 ) THEN
252  DO j = jdfil,jffil
253 #ifdef BLAS
254  CALL dgemm("N", "N", iim, nbniv_loc, iim, 1.0,
255  & matrinvs(1,1,j-jfiltsu+1), iim,
256  & champ_loc(1,j,1), iip1*nlat, 0.0,
257  & champ_fft(1,j-jdfil+1,1), iip1*nlat)
258 #else
259  champ_fft(:iim,j-jdfil+1,:)
260  & =matmul(matrinvs(:,:,j-jfiltsu+1),
261  & champ_loc(:iim,j,:))
262 #endif
263  ENDDO
264 
265  ELSE IF ( griscal ) THEN
266 
267  DO j = jdfil,jffil
268 #ifdef BLAS
269  CALL dgemm("N", "N", iim, nbniv_loc, iim, 1.0,
270  & matriceus(1,1,j-jfiltsu+1), iim,
271  & champ_loc(1,j,1), iip1*nlat, 0.0,
272  & champ_fft(1,j-jdfil+1,1), iip1*nlat)
273 #else
274  champ_fft(:iim,j-jdfil+1,:)
275  & =matmul(matriceus(:,:,j-jfiltsu+1),
276  & champ_loc(:iim,j,:))
277 #endif
278  ENDDO
279 
280  ELSE
281 
282  DO j = jdfil,jffil
283 #ifdef BLAS
284  CALL dgemm("N", "N", iim, nbniv_loc, iim, 1.0,
285  & matricevs(1,1,j-jfiltsv+1), iim,
286  & champ_loc(1,j,1), iip1*nlat, 0.0,
287  & champ_fft(1,j-jdfil+1,1), iip1*nlat)
288 #else
289  champ_fft(:iim,j-jdfil+1,:)
290  & =matmul(matricevs(:,:,j-jfiltsv+1),
291  & champ_loc(:iim,j,:))
292 #endif
293  ENDDO
294 
295  ENDIF
296 
297  ENDIF
298 ! c
299  IF( ifiltre.EQ.2 ) THEN
300 
301 c !-------------------------------------!
302 c ! Dés-agregation des niveau verticaux !
303 c ! uniquement necessaire pour une !
304 c ! execution OpenMP !
305 c !-------------------------------------!
306  ll_nb = 0
307 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
308  DO l = 1, nbniv
309  ll_nb = ll_nb + 1
310  DO j = jdfil,jffil
311  DO i = 1, iim
312  champ( i,j,l ) = (champ_loc(i,j,ll_nb)
313  & + champ_fft(i,j-jdfil+1,ll_nb))
314  & * sdd12(i,sdd2_type)
315  ENDDO
316  ENDDO
317  ENDDO
318 c$OMP END DO NOWAIT
319 
320  ELSE
321 
322  ll_nb = 0
323 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
324  DO l = 1, nbniv_loc
325  ll_nb = ll_nb + 1
326  DO j = jdfil,jffil
327  DO i = 1, iim
328  champ( i,j,l ) = (champ_loc(i,j,ll_nb)
329  & - champ_fft(i,j-jdfil+1,ll_nb))
330  & * sdd12(i,sdd2_type)
331  ENDDO
332  ENDDO
333  ENDDO
334 c$OMP END DO NOWAIT
335 
336  ENDIF
337 
338 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
339  DO l = 1, nbniv
340  DO j = jdfil,jffil
341  champ( iip1,j,l ) = champ( 1,j,l )
342  ENDDO
343  ENDDO
344 c$OMP END DO NOWAIT
345 
346 ccccccccccccccccccccccccccccccccccccccccccccc
347 c Utilisation du filtre FFT
348 ccccccccccccccccccccccccccccccccccccccccccccc
349 
350  ELSE
351 
352 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
353  DO l=1,nbniv
354  DO j=jdfil,jffil
355  DO i = 1, iim
356  champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
357  champ_fft( i,j,l) = champ(i,j,l)
358  ENDDO
359  ENDDO
360  ENDDO
361 c$OMP END DO NOWAIT
362 
363  IF (jdfil<=jffil) THEN
364  IF( ifiltre. eq. -2 ) THEN
365  CALL filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
366  ELSE IF ( griscal ) THEN
367  CALL filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
368  ELSE
369  CALL filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
370  ENDIF
371  ENDIF
372 
373 
374  IF( ifiltre.EQ. 2 ) THEN
375 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
376  DO l=1,nbniv
377  DO j=jdfil,jffil
378  DO i = 1, iim
379  champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
380  & *sdd12(i,sdd2_type)
381  ENDDO
382  ENDDO
383  ENDDO
384 c$OMP END DO NOWAIT
385  ELSE
386 
387 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
388  DO l=1,nbniv
389  DO j=jdfil,jffil
390  DO i = 1, iim
391  champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
392  & *sdd12(i,sdd2_type)
393  ENDDO
394  ENDDO
395  ENDDO
396 c$OMP END DO NOWAIT
397  ENDIF
398 c
399 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
400  DO l=1,nbniv
401  DO j=jdfil,jffil
402 ! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
403  champ( iip1,j,l ) = champ( 1,j,l )
404  ENDDO
405  ENDDO
406 c$OMP END DO NOWAIT
407  ENDIF
408 c Fin de la zone de filtrage
409 
410 
411  ENDDO
412 
413 ! DO j=1,nlat
414 !
415 ! PRINT *,"check FFT ----> Delta(",j,")=",
416 ! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
417 ! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
418 ! ENDDO
419 
420 ! PRINT *,"check FFT ----> Delta(",j,")=",
421 ! & sum(champ-champ_fft)/sum(champ)
422 !
423 
424 c
425  1111 FORMAT(//20x,
426 'ERREUR dans le dimensionnement du tableau CHAMP a & filtrer, sur la grille des scalaires'/)
427  2222 FORMAT(//20x,
428 'ERREUR dans le dimensionnement du tableau CHAMP a fi & ltrer, sur la grille de V ou de Z'/)
429 c$OMP MASTER
430  CALL stop_timer
431 c$OMP END MASTER
432  RETURN
433  END