source: LMDZ4/trunk/libf/filtrez/inifilr.F @ 1088

Last change on this file since 1088 was 986, checked in by Laurent Fairhead, 16 years ago

Inclusion du filtre FFT YM
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.9 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE inifilr
5c
6c    ... H. Upadhyaya, O.Sharma   ...
7c
8      USE mod_filtre_fft
9      IMPLICIT NONE
10c
11c     version 3 .....
12
13c     Correction  le 28/10/97    P. Le Van .
14c  -------------------------------------------------------------------
15#include "dimensions.h"
16#include "paramet.h"
17#include "parafilt.h"
18c  -------------------------------------------------------------------
19#include "comgeom.h"
20#include "coefils.h"
21#include "logic.h"
22#include "serre.h"
23
24      REAL  dlonu(iim),dlatu(jjm)
25      REAL  rlamda( iim ),  eignvl( iim )
26c
27
28      REAL    lamdamax,pi,cof
29      INTEGER i,j,modemax,imx,k,kf,ii
30      REAL dymin,dxmin,colat0
31      REAL eignft(iim,iim), coff
32      REAL matriceun,matriceus,matricevn,matricevs,matrinvn,matrinvs
33      COMMON/matrfil/matriceun(iim,iim,nfilun),matriceus(iim,iim,nfilus)
34     ,             , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs)
35     ,             ,  matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus)
36#ifdef CRAY
37      INTEGER   ISMIN
38      EXTERNAL  ISMIN
39      INTEGER iymin
40      INTEGER ixmineq
41#endif
42      EXTERNAL  inifgn
43c
44c ------------------------------------------------------------
45c   This routine computes the eigenfunctions of the laplacien
46c   on the stretched grid, and the filtering coefficients
47c     
48c  We designate:
49c   eignfn   eigenfunctions of the discrete laplacien
50c   eigenvl  eigenvalues
51c   jfiltn   indexof the last scalar line filtered in NH
52c   jfilts   index of the first line filtered in SH
53c   modfrst  index of the mode from where modes are filtered
54c   modemax  maximum number of modes ( im )
55c   coefil   filtering coefficients ( lamda_max*cos(rlat)/lamda )
56c   sdd      SQRT( dx )
57c     
58c     the modes are filtered from modfrst to modemax
59c     
60c-----------------------------------------------------------
61c
62
63       pi       = 2. * ASIN( 1. )
64
65       DO i = 1,iim
66        dlonu(i) = xprimu( i )
67       ENDDO
68c
69       CALL inifgn(eignvl)
70c
71        print *,' EIGNVL '
72        PRINT 250,eignvl
73250     FORMAT( 1x,5e13.6)
74c
75c compute eigenvalues and eigenfunctions
76c
77c
78c.................................................................
79c
80c  compute the filtering coefficients for scalar lines and
81c  meridional wind v-lines
82c
83c  we filter all those latitude lines where coefil < 1
84c  NO FILTERING AT POLES
85c
86c  colat0 is to be used  when alpha (stretching coefficient)
87c  is set equal to zero for the regular grid case
88c
89c    .......   Calcul  de  colat0   .........
90c     .....  colat0 = minimum de ( 0.5, min dy/ min dx )   ...
91c
92c
93      DO 45 j = 1,jjm
94         dlatu( j ) = rlatu( j ) - rlatu( j+1 )
95 45   CONTINUE
96c
97#ifdef CRAY
98      iymin   = ISMIN( jjm, dlatu, 1 )
99      ixmineq = ISMIN( iim, dlonu, 1 )
100      dymin   = dlatu( iymin )
101      dxmin   = dlonu( ixmineq )
102#else
103      dxmin   =  dlonu(1)
104       DO  i  = 2, iim
105        dxmin = MIN( dxmin,dlonu(i) )
106       ENDDO
107      dymin  = dlatu(1)
108       DO j  = 2, jjm
109        dymin = MIN( dymin,dlatu(j) )
110       ENDDO
111#endif
112c
113c
114      colat0  =  MIN( 0.5, dymin/dxmin )
115c
116      IF( .NOT.fxyhypb.AND.ysinus )  THEN
117           colat0 = 0.6
118c         ...... a revoir  pour  ysinus !   .......
119           alphax = 0.
120      ENDIF
121c
122      PRINT 50, colat0,alphax
123  50  FORMAT(/15x,' Inifilr colat0 alphax ',2e16.7)
124c
125      IF(alphax.EQ.1. )  THEN
126        PRINT *,' Inifilr  alphax doit etre  <  a 1.  Corriger '
127         STOP
128      ENDIF
129c
130      lamdamax = iim / ( pi * colat0 * ( 1. - alphax ) )
131
132cc                        ... Correction  le 28/10/97  ( P.Le Van ) ..
133c
134      DO 71 i = 2,iim
135       rlamda( i ) = lamdamax/ SQRT( ABS( eignvl(i) ) )
136 71   CONTINUE
137c
138
139      DO 72 j = 1,jjm
140            DO 73 i = 1,iim
141            coefilu( i,j )  = 0.0
142            coefilv( i,j )  = 0.0
143            coefilu2( i,j ) = 0.0
144            coefilv2( i,j ) = 0.0
145 73     CONTINUE
146 72   CONTINUE
147
148c
149c    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....
150c    .........................................................
151c
152       modemax = iim
153
154cccc    imx = modemax - 4 * (modemax/iim)
155
156       imx  = iim
157c
158       PRINT *,' TRUNCATION AT ',imx
159c
160      DO 75 j = 2, jjm/2+1
161       cof = COS( rlatu(j) )/ colat0
162            IF ( cof .LT. 1. ) THEN
163          IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) jfiltnu= j
164        ENDIF
165
166       cof = COS( rlatu(jjp1-j+1) )/ colat0
167            IF ( cof .LT. 1. ) THEN
168          IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. )
169     $      jfiltsu= jjp1-j+1
170        ENDIF
171 75   CONTINUE
172c
173      DO 76 j = 1, jjm/2
174       cof = COS( rlatv(j) )/ colat0
175            IF ( cof .LT. 1. ) THEN
176          IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) jfiltnv= j
177        ENDIF
178
179       cof = COS( rlatv(jjm-j+1) )/ colat0
180            IF ( cof .LT. 1. ) THEN
181          IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. )
182     $       jfiltsv= jjm-j+1
183        ENDIF
184 76   CONTINUE
185c                                 
186
187      if ( jfiltnu.LE.0 ) jfiltnu=1
188      IF( jfiltnu.GT. jjm/2 +1 )  THEN
189        PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
190        STOP
191      ENDIF
192
193      IF( jfiltsu.LE.0) jfiltsu=1
194      IF( jfiltsu.GT.  jjm  +1 )  THEN
195        PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
196        STOP
197      ENDIF
198
199      IF( jfiltnv.LE.0) jfiltnv=1
200      IF( jfiltnv.GT. jjm/2    )  THEN
201        PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
202        STOP
203      ENDIF
204
205      IF( jfiltsv.LE.0) jfiltsv=1
206      IF( jfiltsv.GT.     jjm  )  THEN
207        PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
208        STOP
209      ENDIF
210
211       PRINT *,' jfiltnv jfiltsv jfiltnu jfiltsu ' ,
212     *           jfiltnv,jfiltsv,jfiltnu,jfiltsu
213
214c                                 
215c   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
216c................................................................
217c
218c
219      DO 77 j = 1,jjm
220          modfrstu( j ) = iim
221          modfrstv( j ) = iim
222 77   CONTINUE
223c
224      DO 84 j = 2,jfiltnu
225       DO 81 k = 2,modemax
226             cof = rlamda(k) * COS( rlatu(j) )
227         IF ( cof .LT. 1. ) GOTO 82
228 81    CONTINUE
229      GOTO 84
230 82   modfrstu( j ) = k
231c
232          kf = modfrstu( j )
233           DO 83 k = kf , modemax
234        cof = rlamda(k) * COS( rlatu(j) )
235            coefilu(k,j) = cof - 1.
236            coefilu2(k,j) = cof*cof - 1.
237 83    CONTINUE
238 84   CONTINUE
239c                                 
240c
241      DO 89 j = 1,jfiltnv
242c
243       DO 86 k = 2,modemax
244            cof = rlamda(k) * COS( rlatv(j) )
245         IF ( cof .LT. 1. ) GOTO 87
246 86    CONTINUE
247      GOTO 89
248 87   modfrstv( j ) = k
249c
250           kf = modfrstv( j )
251           DO 88 k = kf , modemax
252        cof = rlamda(k) * COS( rlatv(j) )
253            coefilv(k,j) = cof - 1.
254            coefilv2(k,j) = cof*cof - 1.
255 88    CONTINUE
256c
257 89    CONTINUE
258c
259      DO 94 j = jfiltsu,jjm
260       DO 91 k = 2,modemax
261            cof = rlamda(k) * COS( rlatu(j) )
262         IF ( cof .LT. 1. ) GOTO 92
263 91    CONTINUE
264      GOTO 94
265 92   modfrstu( j ) = k
266c
267        kf = modfrstu( j )
268         DO 93 k = kf , modemax
269          cof = rlamda(k) * COS( rlatu(j) )
270          coefilu(k,j) = cof - 1.
271          coefilu2(k,j) = cof*cof - 1.
272 93      CONTINUE
273 94    CONTINUE
274c                                 
275      DO 99 j = jfiltsv,jjm
276       DO 96 k = 2,modemax
277             cof = rlamda(k) * COS( rlatv(j) )
278         IF ( cof .LT. 1. ) GOTO 97
279 96    CONTINUE
280      GOTO 99
281 97   modfrstv( j ) = k
282c
283       kf = modfrstv( j )
284           DO 98 k = kf , modemax
285        cof = rlamda(k) * COS( rlatv(j) )
286            coefilv(k,j) = cof - 1.
287            coefilv2(k,j) = cof*cof - 1.
288 98    CONTINUE
289 99   CONTINUE
290c
291
292       IF(jfiltnv.GE.jjm/2 .OR. jfiltnu.GE.jjm/2)THEN
293
294         IF(jfiltnv.EQ.jfiltsv)jfiltsv=1+jfiltnv
295         IF(jfiltnu.EQ.jfiltsu)jfiltsu=1+jfiltnu
296
297          PRINT *,'jfiltnv jfiltsv jfiltnu jfiltsu' ,
298     *        jfiltnv,jfiltsv,jfiltnu,jfiltsu
299       ENDIF
300
301       PRINT *,'   Modes premiers  v  '
302       PRINT 334,modfrstv
303       PRINT *,'   Modes premiers  u  '
304       PRINT 334,modfrstu
305
306     
307      IF( nfilun.LT. jfiltnu )  THEN
308       PRINT *,' le parametre nfilun utilise pour la matrice ',
309     *   ' matriceun  est trop petit ! '
310       PRINT *,'Le changer dans parafilt.h et le mettre a  ',jfiltnu
311        PRINT *,' Pour information, nfilun,nfilus,nfilvn,nfilvs '
312     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
313     *  ,jfiltnv,jjm-jfiltsv+1
314               STOP
315      ENDIF
316      IF( nfilun.GT. jfiltnu+ 2 )  THEN
317           PRINT *,' le parametre nfilun utilise pour la matrice ',
318     *' matriceun est trop grand ! Gachis de memoire ! '
319       PRINT *,'Le changer dans parafilt.h et le mettre a  ',jfiltnu
320        PRINT *,' Pour information, nfilun,nfilus,nfilvn,nfilvs '
321     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
322     *  ,jfiltnv,jjm-jfiltsv+1
323c              STOP
324      ENDIF
325      IF( nfilus.LT. jjm - jfiltsu +1 )  THEN
326            PRINT *,' le parametre nfilus utilise pour la matrice ',
327     *   ' matriceus  est trop petit !  '
328       PRINT *,' Le changer dans parafilt.h et le mettre a  ',
329     * jjm - jfiltsu + 1
330        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs '
331     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
332     *  ,jfiltnv,jjm-jfiltsv+1
333               STOP
334      ENDIF
335      IF( nfilus.GT. jjm - jfiltsu + 3 )  THEN
336           PRINT *,' le parametre nfilus utilise pour la matrice ',
337     * ' matriceus  est trop grand ! '
338       PRINT *,' Le changer dans parafilt.h et le mettre a  ' ,
339     * jjm - jfiltsu + 1
340        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs '
341     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
342     *  ,jfiltnv,jjm-jfiltsv+1
343c              STOP
344      ENDIF
345      IF( nfilvn.LT. jfiltnv )  THEN
346            PRINT *,' le parametre nfilvn utilise pour la matrice ',
347     *   ' matricevn  est trop petit ! ' 
348       PRINT *,'Le changer dans parafilt.h et le mettre a  ',jfiltnv
349        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs '
350     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
351     *  ,jfiltnv,jjm-jfiltsv+1
352               STOP
353      ENDIF
354      IF( nfilvn.GT. jfiltnv+ 2 )  THEN
355           PRINT *,' le parametre nfilvn utilise pour la matrice ',
356     *' matricevn est trop grand !  Gachis de memoire ! '
357       PRINT *,'Le changer dans parafilt.h et le mettre a  ',jfiltnv
358        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs '
359     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
360     *  ,jfiltnv,jjm-jfiltsv+1
361c              STOP
362      ENDIF
363      IF( nfilvs.LT. jjm - jfiltsv +1 )  THEN
364            PRINT *,' le parametre nfilvs utilise pour la matrice ',
365     *   ' matricevs  est trop petit !  Le changer dans parafilt.h '
366       PRINT *,' Le changer dans parafilt.h et le mettre a  '
367     * , jjm - jfiltsv + 1
368        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs '
369     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
370     *  ,jfiltnv,jjm-jfiltsv+1
371               STOP
372      ENDIF
373      IF( nfilvs.GT. jjm - jfiltsv + 3 )  THEN
374           PRINT *,' le parametre nfilvs utilise pour la matrice ',
375     * ' matricevs  est trop grand ! Gachis de memoire ! '
376       PRINT *,' Le changer dans parafilt.h et le mettre a  '
377     *   ,  jjm - jfiltsv + 1
378        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs '
379     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
380     *  ,jfiltnv,jjm-jfiltsv+1
381c              STOP
382      ENDIF
383
384
385c   ...................................................................
386c
387c   ... Calcul de la matrice filtre 'matriceu'  pour les champs situes
388c                       sur la grille scalaire                 ........
389c   ...................................................................
390c
391        DO j = 2, jfiltnu
392
393         DO i=1,iim
394          coff = coefilu(i,j)
395          IF( i.LT.modfrstu(j) ) coff = 0.
396           DO k=1,iim
397            eignft(i,k) = eignfnv(k,i) * coff
398           ENDDO
399         ENDDO
400#ifdef CRAY
401         CALL MXM( eignfnv,iim,eignft,iim,matriceun(1,1,j),iim )
402#else
403#ifdef BLAS
404         CALL SGEMM ('N', 'N', iim, iim, iim, 1.0,
405     $   eignfnv, iim, eignft, iim, 0.0, matriceun(1,1,j), iim)
406#else
407         DO k = 1, iim
408         DO i = 1, iim
409            matriceun(i,k,j) = 0.0
410            DO ii = 1, iim
411               matriceun(i,k,j) = matriceun(i,k,j)
412     .                          + eignfnv(i,ii)*eignft(ii,k)
413            ENDDO
414         ENDDO
415         ENDDO
416#endif
417#endif
418
419        ENDDO
420
421        DO j = jfiltsu, jjm
422
423         DO i=1,iim
424          coff = coefilu(i,j)
425          IF( i.LT.modfrstu(j) ) coff = 0.
426            DO k=1,iim
427             eignft(i,k) = eignfnv(k,i) * coff
428            ENDDO
429         ENDDO
430#ifdef CRAY
431         CALL MXM(eignfnv,iim,eignft,iim,matriceus(1,1,j-jfiltsu+1),iim)
432#else
433#ifdef BLAS
434      CALL SGEMM ('N', 'N', iim, iim, iim, 1.0,
435     $           eignfnv, iim, eignft, iim, 0.0,
436     $           matriceus(1,1,j-jfiltsu+1), iim)
437#else
438         DO k = 1, iim
439         DO i = 1, iim
440            matriceus(i,k,j-jfiltsu+1) = 0.0
441            DO ii = 1, iim
442               matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1)
443     .                                    + eignfnv(i,ii)*eignft(ii,k)
444            ENDDO
445         ENDDO
446         ENDDO
447#endif
448#endif
449
450        ENDDO
451
452c   ...................................................................
453c
454c   ... Calcul de la matrice filtre 'matricev'  pour les champs situes
455c                       sur la grille   de V ou de Z           ........
456c   ...................................................................
457c
458        DO j = 1, jfiltnv
459
460         DO i = 1, iim
461          coff = coefilv(i,j)
462          IF( i.LT.modfrstv(j) ) coff = 0.
463           DO k = 1, iim
464            eignft(i,k) = eignfnu(k,i) * coff
465           ENDDO
466         ENDDO
467#ifdef CRAY
468          CALL MXM( eignfnu,iim,eignft,iim,matricevn(1,1,j),iim )
469#else
470#ifdef BLAS
471      CALL SGEMM ('N', 'N', iim, iim, iim, 1.0,
472     $     eignfnu, iim, eignft, iim, 0.0, matricevn(1,1,j), iim)
473#else
474         DO k = 1, iim
475         DO i = 1, iim
476            matricevn(i,k,j) = 0.0
477            DO ii = 1, iim
478               matricevn(i,k,j) = matricevn(i,k,j)
479     .                          + eignfnu(i,ii)*eignft(ii,k)
480            ENDDO
481         ENDDO
482         ENDDO
483#endif
484#endif
485
486        ENDDO
487
488        DO j = jfiltsv, jjm
489
490         DO i = 1, iim
491          coff = coefilv(i,j)
492          IF( i.LT.modfrstv(j) ) coff = 0.
493            DO k = 1, iim
494             eignft(i,k) = eignfnu(k,i) * coff
495            ENDDO
496         ENDDO
497#ifdef CRAY
498         CALL MXM(eignfnu,iim,eignft,iim,matricevs(1,1,j-jfiltsv+1),iim)
499#else
500#ifdef BLAS
501      CALL SGEMM ('N', 'N', iim, iim, iim, 1.0,
502     $           eignfnu, iim, eignft, iim, 0.0,
503     $           matricevs(1,1,j-jfiltsv+1), iim)
504#else
505         DO k = 1, iim
506         DO i = 1, iim
507            matricevs(i,k,j-jfiltsv+1) = 0.0
508            DO ii = 1, iim
509               matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1)
510     .                          + eignfnu(i,ii)*eignft(ii,k)
511            ENDDO
512         ENDDO
513         ENDDO
514#endif
515#endif
516
517        ENDDO
518
519c   ...................................................................
520c
521c   ... Calcul de la matrice filtre 'matrinv'  pour les champs situes
522c              sur la grille scalaire , pour le filtre inverse ........
523c   ...................................................................
524c
525        DO j = 2, jfiltnu
526
527         DO i = 1,iim
528          coff = coefilu(i,j)/ ( 1. + coefilu(i,j) )
529          IF( i.LT.modfrstu(j) ) coff = 0.
530           DO k=1,iim
531            eignft(i,k) = eignfnv(k,i) * coff
532           ENDDO
533         ENDDO
534#ifdef CRAY
535          CALL MXM( eignfnv,iim,eignft,iim,matrinvn(1,1,j),iim )
536#else
537#ifdef BLAS
538      CALL SGEMM ('N', 'N', iim, iim, iim, 1.0,
539     $     eignfnv, iim, eignft, iim, 0.0, matrinvn(1,1,j), iim)
540#else
541         DO k = 1, iim
542         DO i = 1, iim
543            matrinvn(i,k,j) = 0.0
544            DO ii = 1, iim
545               matrinvn(i,k,j) = matrinvn(i,k,j)
546     .                          + eignfnv(i,ii)*eignft(ii,k)
547            ENDDO
548         ENDDO
549         ENDDO
550#endif
551#endif
552
553        ENDDO
554
555        DO j = jfiltsu, jjm
556
557         DO i = 1,iim
558          coff = coefilu(i,j) / ( 1. + coefilu(i,j) )
559          IF( i.LT.modfrstu(j) ) coff = 0.
560            DO k=1,iim
561             eignft(i,k) = eignfnv(k,i) * coff
562            ENDDO
563         ENDDO
564#ifdef CRAY
565         CALL MXM(eignfnv,iim,eignft,iim,matrinvs(1,1,j-jfiltsu+1),iim)
566#else
567#ifdef BLAS
568      CALL SGEMM ('N', 'N', iim, iim, iim, 1.0,
569     $ eignfnv, iim, eignft, iim, 0.0, matrinvs(1,1,j-jfiltsu+1), iim)
570#else
571         DO k = 1, iim
572         DO i = 1, iim
573            matrinvs(i,k,j-jfiltsu+1) = 0.0
574            DO ii = 1, iim
575               matrinvs(i,k,j-jfiltsu+1) = matrinvs(i,k,j-jfiltsu+1)
576     .                          + eignfnv(i,ii)*eignft(ii,k)
577            ENDDO
578         ENDDO
579         ENDDO
580#endif
581#endif
582
583        ENDDO
584
585c   ...................................................................
586
587c
588       IF (use_filtre_fft) THEN
589       CALL Init_filtre_fft(coefilu,modfrstu,jfiltnu,jfiltsu,
590     &                      coefilv,modfrstv,jfiltnv,jfiltsv)
591       ENDIF
592       
593334    FORMAT(1x,24i3)
594755    FORMAT(1x,6f10.3,i3)
595
596       RETURN
597       END
Note: See TracBrowser for help on using the repository browser.