source: LMDZ.3.3/trunk/libf/filtrez/inifilr.F

Last change on this file was 2, checked in by lmdz, 25 years ago

Initial revision

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