LMDZ
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_lmdz, 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,'ERREUR dans le dimensionnement du tableau CHAMP a
426  & filtrer, sur la grille des scalaires'/)
427  2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
428  & 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
sddu
!$Id!COMMON coefils sddu(iim)
jfiltsu
!$Id!COMMON coefils jfiltsu
Definition: coefils.h:4
unsddu
!$Id!COMMON coefils unsddu(iim)
mod_filtre_fft::filtre_u_fft
subroutine filtre_u_fft(vect_inout, nlat, jj_begin, jj_end, nbniv)
Definition: mod_filtre_fft.F90:118
filtreg_mod::matriceus
real, dimension(:,:,:), allocatable matriceus
Definition: filtreg_mod.F90:6
mod_filtre_fft::use_filtre_fft
logical, save use_filtre_fft
Definition: mod_filtre_fft.F90:7
filtreg_mod::matrinvs
real, dimension(:,:,:), allocatable matrinvs
Definition: filtreg_mod.F90:7
filtreg_p
subroutine filtreg_p(champ, ibeg, iend, nlat, nbniv, ifiltre, iaire, griscal, iter)
Definition: filtreg_p.F:5
timer_filtre::start_timer
subroutine, public start_timer
Definition: timer_filtre.F90:15
mod_filtre_fft::filtre_v_fft
subroutine filtre_v_fft(vect_inout, nlat, jj_begin, jj_end, nbniv)
Definition: mod_filtre_fft.F90:179
jfiltsv
!$Id!COMMON coefils jfiltsv
Definition: coefils.h:4
jfiltnu
!$Id!COMMON coefils jfiltnu
Definition: coefils.h:4
timer_filtre::stop_timer
subroutine, public stop_timer
Definition: timer_filtre.F90:22
false
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
mod_filtre_fft
Definition: mod_filtre_fft.F90:5
parallel_lmdz
Definition: parallel_lmdz.F90:4
jjp1
!$Header jjp1
Definition: paramet.h:14
filtreg_mod
Definition: filtreg_mod.F90:4
x
!$Header!c c INCLUDE fxyprim h c c c Fonctions in line c c REAL fyprim REAL rj c c il faut la calculer avant d appeler ces fonctions c c c Fonctions a changer selon x(x) et y(y) choisis.c-----------------------------------------------------------------c c.....ici
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
sddv
!$Id!COMMON coefils sddv(iim)&&
unsddv
!$Id!COMMON coefils unsddv(iim)
filtreg_mod::matrinvn
real, dimension(:,:,:), allocatable matrinvn
Definition: filtreg_mod.F90:7
parallel_lmdz::omp_chunk
integer, save omp_chunk
Definition: parallel_lmdz.F90:46
jfiltnv
!$Id!COMMON coefils jfiltnv
Definition: coefils.h:4
timer_filtre::init_timer
subroutine, public init_timer
Definition: timer_filtre.F90:10
iim
c c zjulian c cym CALL iim cym klev iim
Definition: ini_bilKP_ave.h:24
filtreg_mod::matriceun
real, dimension(:,:,:), allocatable matriceun
Definition: filtreg_mod.F90:6
filtreg_mod::matricevn
real, dimension(:,:,:), allocatable matricevn
Definition: filtreg_mod.F90:6
filtreg_mod::matricevs
real, dimension(:,:,:), allocatable matricevs
Definition: filtreg_mod.F90:7
timer_filtre
Definition: timer_filtre.F90:1
mod_filtre_fft::filtre_inv_fft
subroutine filtre_inv_fft(vect_inout, nlat, jj_begin, jj_end, nbniv)
Definition: mod_filtre_fft.F90:241