source: LMDZ6/trunk/libf/dyn3dmem/filtreg_p.F @ 4660

Last change on this file since 4660 was 4593, checked in by yann meurdesoif, 18 months ago

Replace #include (c preprocessor) by INCLUDE (fortran keyword)

in phylmd (except rrtm and ecrad) filtrez, dy3dmem and dyn3dcommon

Other directories will follow
YM

  • 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: 13.1 KB
RevLine 
[1632]1
2
3      SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv,
4     &     ifiltre, iaire, griscal ,iter)
[1823]5      USE parallel_lmdz, only : OMP_CHUNK
[1632]6      USE mod_filtre_fft
7      USE timer_filtre
8     
9      USE filtreg_mod
10     
11      IMPLICIT NONE
12     
13c=======================================================================
14c
15c   Auteur: P. Le Van        07/10/97
16c   ------
17c
18c   Objet: filtre matriciel longitudinal ,avec les matrices precalculees
19c                     pour l'operateur  Filtre    .
20c   ------
21c
22c   Arguments:
23c   ----------
24c
25c     
26c      ibeg..iend            lattitude a filtrer
27c      nlat                  nombre de latitudes du champ
28c      nbniv                 nombre de niveaux verticaux a filtrer
29c      champ(iip1,nblat,nbniv)  en entree : champ a filtrer
30c                            en sortie : champ filtre
31c      ifiltre               +1  Transformee directe
32c                            -1  Transformee inverse
33c                            +2  Filtre directe
34c                            -2  Filtre inverse
35c
36c      iaire                 1   si champ intensif
37c                            2   si champ extensif (pondere par les aires)
38c
39c      iter                  1   filtre simple
40c
41c=======================================================================
42c
43c
44c                      Variable Intensive
45c                ifiltre = 1     filtre directe
46c                ifiltre =-1     filtre inverse
47c
48c                      Variable Extensive
49c                ifiltre = 2     filtre directe
50c                ifiltre =-2     filtre inverse
51c
52c
[4593]53      INCLUDE "dimensions.h"
54      INCLUDE "paramet.h"
55      INCLUDE "coefils.h"
[1632]56c
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.
71c$OMP THREADPRIVATE(first)
72
73      REAL, DIMENSION(iip1,nlat,nbniv) :: champ_loc
74      INTEGER :: ll_nb, nbniv_loc
75      REAL, SAVE :: sdd12(iim,4)
76c$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
95c$OMP MASTER     
96      CALL start_timer
97c$OMP END MASTER
98
99c-------------------------------------------------------c
100
[4469]101      IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
102     & CALL abort_gcm("fitreg_p","Pas de transformee simple
103     &dans cette version",1)
[1632]104     
105      IF( iter.EQ. 2 )  THEN
106         PRINT *,' Pas d iteration du filtre dans cette version !'
107     &        , ' Utiliser old_filtreg et repasser !'
[4469]108         CALL abort_gcm("fitreg_p","stopped",1)
[1632]109      ENDIF
110
111      IF( ifiltre.EQ. -2 .AND..NOT.griscal )     THEN
112         PRINT *,' Cette routine ne calcule le filtre inverse que '
113     &        , ' sur la grille des scalaires !'
[4469]114         CALL abort_gcm("fitreg_p","stopped",1)
[1632]115      ENDIF
116
117      IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 )  THEN
118         PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
119     &        , ' corriger et repasser !'
[4469]120         CALL abort_gcm("fitreg_p","stopped",1)
121       ENDIF
[1632]122c
123
124      iim2   = iim * iim
125      immjm  = iim * jjm
126c
127c
128      IF( griscal )   THEN
129         IF( nlat. NE. jjp1 )  THEN
[4469]130            CALL abort_gcm("fitreg_p","nlat. NE. jjp1",1)
[1632]131         ELSE
132c     
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
140c
141            jdfil1 = 2
142            jffil1 = jfiltnu
143            jdfil2 = jfiltsu
144            jffil2 = jjm
145         ENDIF
146      ELSE
147         IF( nlat.NE.jjm )  THEN
[4469]148            CALL abort_gcm("fitreg_p","nlat. NE. jjm",1)
[1632]149         ELSE
150c
151            IF( iaire.EQ.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
158c     
159            jdfil1 = 1
160            jffil1 = jfiltnv
161            jdfil2 = jfiltsv
162            jffil2 = jjm
163         ENDIF
164      ENDIF
165c     
166      DO hemisph = 1, 2
167c     
168         IF ( hemisph.EQ.1 )  THEN
169cym
170            jdfil = max(jdfil1,ibeg)
171            jffil = min(jffil1,iend)
172         ELSE
173cym
174            jdfil = max(jdfil2,ibeg)
175            jffil = min(jffil2,iend)
176         ENDIF
177
178
179cccccccccccccccccccccccccccccccccccccccccccc
180c Utilisation du filtre classique
181cccccccccccccccccccccccccccccccccccccccccccc
182
183         IF (.NOT. use_filtre_fft) THEN
184     
185c     !---------------------------------!
186c     ! Agregation des niveau verticaux !
187c     ! uniquement necessaire pour une  !
188c     ! execution OpenMP                !
189c     !---------------------------------!
190            ll_nb = 0
191c$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
201c$OMP END DO NOWAIT
202
203            nbniv_loc = ll_nb
204
205            IF( hemisph.EQ.1 )      THEN
206               
207               IF( ifiltre.EQ.-2 )   THEN
208                  DO j = jdfil,jffil
[1673]209#ifdef BLAS
[1632]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)
[1673]214#else
[1680]215                     champ_fft(:iim,j-jdfil+1,:)
[1673]216     &                    =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
217#endif
[1632]218                  ENDDO
219                 
220               ELSE IF ( griscal )     THEN
221                  DO j = jdfil,jffil
[1673]222#ifdef BLAS
[1632]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)
[1673]227#else
[1680]228                     champ_fft(:iim,j-jdfil+1,:)
[1673]229     &                    =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
230#endif
[1632]231                  ENDDO
232                 
233               ELSE
234                  DO j = jdfil,jffil
[1673]235#ifdef BLAS
[1632]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)
[1673]240#else
[1680]241                     champ_fft(:iim,j-jdfil+1,:)
[1673]242     &                    =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
243#endif
[1632]244                  ENDDO
245                 
246               ENDIF
247               
248            ELSE
249               
250               IF( ifiltre.EQ.-2 )   THEN
251                  DO j = jdfil,jffil
[1673]252#ifdef BLAS
[1632]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)
[1673]257#else
[1680]258                     champ_fft(:iim,j-jdfil+1,:)
[1673]259     &                    =matmul(matrinvs(:,:,j-jfiltsu+1),
260     &                            champ_loc(:iim,j,:))
261#endif
[1632]262                  ENDDO
263                 
264               ELSE IF ( griscal )     THEN
265                 
266                  DO j = jdfil,jffil
[1673]267#ifdef BLAS
[1632]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)
[1673]272#else
[1680]273                     champ_fft(:iim,j-jdfil+1,:)
[1673]274     &                    =matmul(matriceus(:,:,j-jfiltsu+1),
275     &                            champ_loc(:iim,j,:))
276#endif
[1632]277                  ENDDO
278                 
279               ELSE
280                 
281                  DO j = jdfil,jffil
[1673]282#ifdef BLAS
[1632]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)
[1673]287#else
[1680]288                     champ_fft(:iim,j-jdfil+1,:)
[1673]289     &                    =matmul(matricevs(:,:,j-jfiltsv+1),
290     &                            champ_loc(:iim,j,:))
291#endif
[1632]292                  ENDDO
293                 
294               ENDIF
295               
296            ENDIF
297!     c     
298            IF( ifiltre.EQ.2 )  THEN
299               
300c     !-------------------------------------!
301c     ! Dés-agregation des niveau verticaux !
302c     ! uniquement necessaire pour une      !
303c     ! execution OpenMP                    !
304c     !-------------------------------------!
305               ll_nb = 0
306c$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
317c$OMP END DO NOWAIT
318               
319            ELSE
320               
321               ll_nb = 0
322c$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
333c$OMP END DO NOWAIT
334               
335            ENDIF
336           
337c$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
343c$OMP END DO NOWAIT
344           
345ccccccccccccccccccccccccccccccccccccccccccccc
346c Utilisation du filtre FFT
347ccccccccccccccccccccccccccccccccccccccccccccc
348       
349         ELSE
350       
351c$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
360c$OMP END DO NOWAIT
361
362            IF (jdfil<=jffil) THEN
363               IF( ifiltre. EQ. -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.EQ. 2 )  THEN
374c$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
383c$OMP END DO NOWAIT       
384            ELSE
385       
386c$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
395c$OMP END DO NOWAIT         
396            ENDIF
397c
398c$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
405c$OMP END DO NOWAIT             
406         ENDIF
407c 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     
423c
424 1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
425     &     filtrer, sur la grille des scalaires'/)
426 2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
427     &     ltrer, sur la grille de V ou de Z'/)
428c$OMP MASTER     
429      CALL stop_timer
430c$OMP END MASTER
431      RETURN
432      END
Note: See TracBrowser for help on using the repository browser.