source: LMDZ5/trunk/libf/dyn3dmem/mod_filtreg_p.F @ 1632

Last change on this file since 1632 was 1632, checked in by Laurent Fairhead, 12 years ago

Import initial du répertoire dyn3dmem

Attention! ceci n'est qu'une version préliminaire du code "basse mémoire":
le code contenu dans ce répertoire est basé sur la r1320 et a donc besoin
d'être mis à jour par rapport à la dynamique parallèle d'aujourd'hui.
Ce code est toutefois mis à disposition pour circonvenir à des problèmes
de mémoire que certaines configurations du modèle pourraient rencontrer.
Dans l'état, il compile et tourne sur vargas et au CCRT

Initial import of dyn3dmem

Warning! this is just a preliminary version of the memory light code:
it is based on r1320 of the code and thus needs to be updated before
it can replace the present dyn3dpar code. It is nevertheless put at your
disposal to circumvent some memory problems some LMDZ configurations may
encounter. In its present state, it will compile and run on vargas and CCRT

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