LMDZ
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_lmdz, only : omp_chunk
8  USE mod_filtre_fft_loc, ONLY: use_filtre_fft, filtre_u_fft,
9  & filtre_v_fft, filtre_inv_fft
10  USE timer_filtre, ONLY: init_timer, start_timer, stop_timer
11 
12  USE filtreg_mod, ONLY: matrinvn, matrinvs, matriceun, matriceus,
13  & matricevn, matricevs
14 
15  IMPLICIT NONE
16 
17 c=======================================================================
18 c
19 c Auteur: P. Le Van 07/10/97
20 c ------
21 c
22 c Objet: filtre matriciel longitudinal ,avec les matrices precalculees
23 c pour l'operateur Filtre .
24 c ------
25 c
26 c Arguments:
27 c ----------
28 c
29 c
30 c ibeg..iend lattitude a filtrer
31 c nlat nombre de latitudes du champ
32 c nbniv nombre de niveaux verticaux a filtrer
33 c champ(iip1,nblat,nbniv) en entree : champ a filtrer
34 c en sortie : champ filtre
35 c ifiltre +1 Transformee directe
36 c -1 Transformee inverse
37 c +2 Filtre directe
38 c -2 Filtre inverse
39 c
40 c iaire 1 si champ intensif
41 c 2 si champ extensif (pondere par les aires)
42 c
43 c iter 1 filtre simple
44 c
45 c=======================================================================
46 c
47 c
48 c Variable Intensive
49 c ifiltre = 1 filtre directe
50 c ifiltre =-1 filtre inverse
51 c
52 c Variable Extensive
53 c ifiltre = 2 filtre directe
54 c ifiltre =-2 filtre inverse
55 c
56 c
57 #include "dimensions.h"
58 #include "paramet.h"
59 #include "coefils.h"
60 c
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 c$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 c$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 
88  IF (first) THEN
89  sdd12(1:iim,type_sddu) = sddu(1:iim)
90  sdd12(1:iim,type_sddv) = sddv(1:iim)
91  sdd12(1:iim,type_unsddu) = unsddu(1:iim)
92  sdd12(1:iim,type_unsddv) = unsddv(1:iim)
93 
94  CALL init_timer
95  first=.false.
96  ENDIF
97 
98 c$OMP MASTER
99  CALL start_timer
100 c$OMP END MASTER
101 
102 c-------------------------------------------------------c
103 
104  IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
105  & stop'Pas de transformee simple dans cette version'
106 
107  IF( iter.EQ. 2 ) THEN
108  print *,' Pas d iteration du filtre dans cette version !'
109  & , ' Utiliser old_filtreg et repasser !'
110  stop
111  ENDIF
112 
113  IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
114  print *,' Cette routine ne calcule le filtre inverse que '
115  & , ' sur la grille des scalaires !'
116  stop
117  ENDIF
118 
119  IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
120  print *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
121  & , ' corriger et repasser !'
122  stop
123  ENDIF
124 c
125 
126  iim2 = iim * iim
127  immjm = iim * jjm
128 c
129 c
130  IF( griscal ) THEN
131  IF( nlat. ne. jjp1 ) THEN
132  print 1111
133  stop
134  ELSE
135 c
136  IF( iaire.EQ.1 ) THEN
137  sdd1_type = type_sddv
138  sdd2_type = type_unsddv
139  ELSE
140  sdd1_type = type_unsddv
141  sdd2_type = type_sddv
142  ENDIF
143 c
144  jdfil1 = 2
145  jffil1 = jfiltnu
146  jdfil2 = jfiltsu
147  jffil2 = jjm
148  ENDIF
149  ELSE
150  IF( nlat.NE.jjm ) THEN
151  print 2222
152  stop
153  ELSE
154 c
155  IF( iaire.EQ.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 c
163  jdfil1 = 1
164  jffil1 = jfiltnv
165  jdfil2 = jfiltsv
166  jffil2 = jjm
167  ENDIF
168  ENDIF
169 c
170  DO hemisph = 1, 2
171 c
172  IF ( hemisph.EQ.1 ) THEN
173 cym
174  jdfil = max(jdfil1,ibeg)
175  jffil = min(jffil1,iend)
176  ELSE
177 cym
178  jdfil = max(jdfil2,ibeg)
179  jffil = min(jffil2,iend)
180  ENDIF
181 
182 
183 cccccccccccccccccccccccccccccccccccccccccccc
184 c Utilisation du filtre classique
185 cccccccccccccccccccccccccccccccccccccccccccc
186 
187  IF (.NOT. use_filtre_fft) THEN
188 
189 c !---------------------------------!
190 c ! Agregation des niveau verticaux !
191 c ! uniquement necessaire pour une !
192 c ! execution OpenMP !
193 c !---------------------------------!
194  ll_nb = 0
195 c$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 c$OMP END DO NOWAIT
206 
207  nbniv_loc = ll_nb
208 
209  IF( hemisph.EQ.1 ) THEN
210 
211  IF( ifiltre.EQ.-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.EQ.-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.EQ.2 ) THEN
306 
307 c !-------------------------------------!
308 c ! Dés-agregation des niveau verticaux !
309 c ! uniquement necessaire pour une !
310 c ! execution OpenMP !
311 c !-------------------------------------!
312  ll_nb = 0
313 c$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 c$OMP END DO NOWAIT
325 
326  ELSE
327 
328  ll_nb = 0
329 c$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 c$OMP END DO NOWAIT
341 
342  ENDIF
343 
344 c$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 c$OMP END DO NOWAIT
352 
353 ccccccccccccccccccccccccccccccccccccccccccccc
354 c Utilisation du filtre FFT
355 ccccccccccccccccccccccccccccccccccccccccccccc
356 
357  ELSE
358 
359 c$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 c$OMP END DO NOWAIT
369 
370  IF (jdfil<=jffil) THEN
371  IF( ifiltre. eq. -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 
381  IF( ifiltre.EQ. 2 ) THEN
382 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
383  DO l=1,nbniv
384  DO j=jdfil,jffil
385  DO i = 1, iim
386  champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
387  & *sdd12(i,sdd2_type)
388  ENDDO
389  ENDDO
390  ENDDO
391 c$OMP END DO NOWAIT
392  ELSE
393 
394 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
395  DO l=1,nbniv
396  DO j=jdfil,jffil
397  DO i = 1, iim
398  champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
399  & *sdd12(i,sdd2_type)
400  ENDDO
401  ENDDO
402  ENDDO
403 c$OMP END DO NOWAIT
404  ENDIF
405 c
406 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
407  DO l=1,nbniv
408  DO j=jdfil,jffil
409 ! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
410  ! add redundant longitude
411  champ( iip1,j,l ) = champ( 1,j,l )
412  ENDDO
413  ENDDO
414 c$OMP END DO NOWAIT
415  ENDIF
416 c Fin de la zone de filtrage
417 
418 
419  ENDDO
420 
421 ! DO j=1,nlat
422 !
423 ! PRINT *,"check FFT ----> Delta(",j,")=",
424 ! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
425 ! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
426 ! ENDDO
427 
428 ! PRINT *,"check FFT ----> Delta(",j,")=",
429 ! & sum(champ-champ_fft)/sum(champ)
430 !
431 
432 c
433  1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
434  & filtrer, sur la grille des scalaires'/)
435  2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
436  & ltrer, sur la grille de V ou de Z'/)
437 c$OMP MASTER
438  CALL stop_timer
439 c$OMP END MASTER
440  RETURN
441  END SUBROUTINE filtreg_p
442  END MODULE mod_filtreg_p
443 
sddu
!$Id!COMMON coefils sddu(iim)
mod_filtre_fft_loc::filtre_u_fft
subroutine filtre_u_fft(vect_inout, jjb, jje, jj_begin, jj_end, nbniv)
Definition: mod_filtre_fft_loc.F90:107
jfiltsu
!$Id!COMMON coefils jfiltsu
Definition: coefils.h:4
unsddu
!$Id!COMMON coefils unsddu(iim)
filtreg_mod::matriceus
real, dimension(:,:,:), allocatable matriceus
Definition: filtreg_mod.F90:6
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
jfiltsv
!$Id!COMMON coefils jfiltsv
Definition: coefils.h:4
mod_filtre_fft_loc::use_filtre_fft
logical, save use_filtre_fft
Definition: mod_filtre_fft_loc.F90:3
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
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
mod_filtreg_p
Definition: mod_filtreg_p.F:1
parallel_lmdz::omp_chunk
integer, save omp_chunk
Definition: parallel_lmdz.F90:46
mod_filtre_fft_loc
Definition: mod_filtre_fft_loc.F90:1
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
timer_filtre
Definition: timer_filtre.F90:1