source: LMDZ6/trunk/libf/dyn3dmem/mod_filtreg_p.F90 @ 5296

Last change on this file since 5296 was 5285, checked in by abarral, 9 months ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules

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