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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 25 - (hide annotations)
Fri Mar 5 16:43:45 2010 UTC (14 years, 2 months ago) by guez
Original Path: trunk/libf/filtrez/inifilr.f
File size: 15424 byte(s)
Simplified "etat0_lim.sh" and "gcm.sh" because the full versions
depended on personal arrangements for directories and machines.

Translated included files into modules. Encapsulated procedures into modules.

Moved variables from module "comgeom" to local variables of
"inigeom". Deleted some unused variables in "comgeom".

Moved variable "day_ini" from module "temps" to module "dynetat0_m".

Removed useless test on variable "time" and useless "close" statement
in procedure "leapfrog".

Removed useless call to "inigeom" in procedure "limit".

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 guez 5 use parafilt
14 guez 25 use coefils
15 guez 3 IMPLICIT NONE
16     c
17     c version 3 .....
18    
19     c Correction le 28/10/97 P. Le Van .
20     c -------------------------------------------------------------------
21     c -------------------------------------------------------------------
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