/[lmdze]/trunk/Sources/filtrez/inifilr.f
ViewVC logotype

Annotation of /trunk/Sources/filtrez/inifilr.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21