source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/lmdz_filtreg_p.F90 @ 5107

Last change on this file since 5107 was 5107, checked in by abarral, 3 months ago

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