source: LMDZ4/branches/IPSL-CM4_IPCC_branch/libf/filtrez/inifilr.F @ 5095

Last change on this file since 5095 was 524, checked in by lmdzadmin, 21 years ago

Initial revision

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