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