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

Last change on this file since 5280 was 5272, checked in by abarral, 7 weeks ago

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