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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/filtrez/inifilr.f90 revision 32 by guez, Tue Apr 6 17:52:58 2010 UTC trunk/Sources/filtrez/inifilr.f revision 139 by guez, Tue May 26 17:46:03 2015 UTC
# Line 1  Line 1 
1  SUBROUTINE inifilr  module inifilr_m
   
   ! From filtrez/inifilr.F,v 1.1.1.1 2004/05/19 12:53:09  
   ! H. Upadhyaya, O.Sharma  
   
   !   This routine computes the eigenfunctions of the laplacien            
   !   on the stretched grid, and the filtering coefficients                
   
   !  We designate:                                                          
   !   eignfn   eigenfunctions of the discrete laplacien                    
   !   eigenvl  eigenvalues                                                  
   !   jfiltn   indexof the last scalar line filtered in NH                  
   !   jfilts   index of the first line filtered in SH                      
   !   modfrst  index of the mode from where modes are filtered              
   !   modemax  maximum number of modes ( im )                              
   !   coefil   filtering coefficients ( lamda_max*cos(rlat)/lamda )        
   !   sdd      SQRT( dx )                                                  
   
   !     the modes are filtered from modfrst to modemax                      
   
   USE dimens_m  
   USE paramet_m  
   USE logic  
   USE comgeom  
   USE serre  
   USE parafilt  
   USE coefils  
2    
3    IMPLICIT NONE    IMPLICIT NONE
4    
5    REAL dlonu(iim), dlatu(jjm)    INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv
   REAL rlamda(iim), eignvl(iim)  
6    
7    REAL lamdamax, pi, cof    ! North:
8    INTEGER i, j, modemax, imx, k, kf, ii    real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :)
9    REAL dymin, dxmin, colat0    ! (iim, iim, 2:jfiltnu)
10    REAL eignft(iim,iim), coff  
11    EXTERNAL inifgn    real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)
12    
13    !-----------------------------------------------------------                ! South:
14      real, allocatable:: matriceus(:, :, :), matrinvs(:, :, :)
15      ! (iim, iim, jfiltsu:jjm)
16    pi = 2.*asin(1.)  
17      real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
18    DO i = 1, iim  
19       dlonu(i) = xprimu(i)  contains
20    END DO  
21      SUBROUTINE inifilr
22    CALL inifgn(eignvl)  
23        ! From filtrez/inifilr.F, version 1.1.1.1 2004/05/19 12:53:09
24    PRINT *, ' EIGNVL '      ! H. Upadhyaya, O. Sharma
25    PRINT 250, eignvl  
26  250 FORMAT (1X,5E13.6)      ! This routine computes the eigenfunctions of the laplacian on the
27        ! stretched grid, and the filtering coefficients.
28    ! compute eigenvalues and eigenfunctions                                      ! We designate:
29        ! eignfn eigenfunctions of the discrete laplacian
30        ! eigenvl eigenvalues
31    !.................................................................            ! jfiltn index of the last scalar line filtered in NH
32        ! jfilts index of the first line filtered in SH
33    !  compute the filtering coefficients for scalar lines and                    ! modfrst index of the mode from where modes are filtered
34    !  meridional wind v-lines                                                    ! modemax maximum number of modes (im)
35        ! coefil filtering coefficients (lamda_max * cos(rlat) / lamda)
36    !  we filter all those latitude lines where coefil < 1                        ! sdd SQRT(dx)
37    !  NO FILTERING AT POLES                                                  
38        ! The modes are filtered from modfrst to modemax.
39    !  colat0 is to be used  when alpha (stretching coefficient)              
40    !  is set equal to zero for the regular grid case                            USE coefils, ONLY : coefilu, coefilu2, coefilv, coefilv2, eignfnu, &
41             eignfnv, modfrstu, modfrstv
42    !    .......   Calcul  de  colat0   .........                                USE dimens_m, ONLY : iim, jjm
43    !     .....  colat0 = minimum de ( 0.5, min dy/ min dx )   ...                USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
44        use inifgn_m, only: inifgn
45        use nr_util, only: pi
46    DO  j = 1, jjm  
47       dlatu(j) = rlatu(j) - rlatu(j+1)      ! Local:
48    END DO      REAL dlatu(jjm)
49        REAL rlamda(2: iim), eignvl(iim)
50    dxmin = dlonu(1)  
51    DO i = 2, iim      REAL lamdamax, cof
52       dxmin = min(dxmin,dlonu(i))      INTEGER i, j, modemax, imx, k, kf
53    END DO      REAL dymin, colat0
54    dymin = dlatu(1)      REAL eignft(iim, iim), coff
55    DO j = 2, jjm  
56       dymin = min(dymin,dlatu(j))      !-----------------------------------------------------------
57    END DO  
58        print *, "Call sequence information: inifilr"
59    
60    colat0 = min(0.5,dymin/dxmin)      CALL inifgn(eignvl)
61    
62    IF ( .NOT. fxyhypb .AND. ysinus) THEN      PRINT *, 'EIGNVL '
63       colat0 = 0.6      PRINT "(1X, 5E13.6)", eignvl
64       !         ...... a revoir  pour  ysinus !   .......                      
65       alphax = 0.      ! compute eigenvalues and eigenfunctions
66    END IF      ! compute the filtering coefficients for scalar lines and
67        ! meridional wind v-lines
68    PRINT 50, colat0, alphax      ! we filter all those latitude lines where coefil < 1
69  50 FORMAT (/15X,' Inifilr colat0 alphax ',2E16.7)      ! NO FILTERING AT POLES
70        ! colat0 is to be used when alpha (stretching coefficient)
71    IF (alphax==1.) THEN      ! is set equal to zero for the regular grid case
72       PRINT *, ' Inifilr  alphax doit etre  <  a 1.  Corriger '  
73       STOP 1      ! Calcul de colat0
74    END IF  
75        DO j = 1, jjm
76    lamdamax = iim/(pi*colat0*(1.-alphax))         dlatu(j) = rlatu(j) - rlatu(j+1)
77        END DO
78    DO  i = 2, iim  
79       rlamda(i) = lamdamax/sqrt(abs(eignvl(i)))      dymin = dlatu(1)
80    END DO      DO j = 2, jjm
81           dymin = min(dymin, dlatu(j))
82        END DO
83    DO  j = 1, jjm  
84       DO  i = 1, iim      colat0 = min(0.5, dymin / minval(xprimu(:iim)))
85          coefilu(i,j) = 0.0  
86          coefilv(i,j) = 0.0      PRINT *, 'colat0 = ', colat0
87          coefilu2(i,j) = 0.0  
88          coefilv2(i,j) = 0.0      lamdamax = iim / (pi * colat0 / grossismx)
89       end DO      rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))
90    END DO  
91        DO j = 1, jjm
92           DO i = 1, iim
93    !    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....                      coefilu(i, j) = 0.
94    !    .........................................................                      coefilv(i, j) = 0.
95              coefilu2(i, j) = 0.
96    modemax = iim            coefilv2(i, j) = 0.
97           end DO
98    !ccc    imx = modemax - 4 * (modemax/iim)                                    END DO
99    
100    imx = iim      ! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
101    
102    PRINT *, ' TRUNCATION AT ', imx      modemax = iim
103        imx = iim
104    DO  j = 2, jjm/2 + 1  
105       cof = cos(rlatu(j))/colat0      PRINT *, 'TRUNCATION AT ', imx
106       IF (cof<1.) THEN  
107          IF (rlamda(imx)*cos(rlatu(j))<1.) jfiltnu = j      DO j = 2, jjm / 2 + 1
108       END IF         IF (cos(rlatu(j)) / colat0 < 1. &
109                .and. rlamda(imx) * cos(rlatu(j)) < 1.) jfiltnu = j
110       cof = cos(rlatu(jjp1-j+1))/colat0  
111       IF (cof<1.) THEN         IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &
112          IF (rlamda(imx)*cos(rlatu(jjp1-j+1))<1.) jfiltsu = jjp1 - j + 1              .and. rlamda(imx) * cos(rlatu(jjm - j + 2)) < 1.) &
113       END IF              jfiltsu = jjm - j + 2
114    END DO      END DO
115    
116    DO  j = 1, jjm/2      DO j = 1, jjm/2
117       cof = cos(rlatv(j))/colat0         cof = cos(rlatv(j))/colat0
118       IF (cof<1.) THEN         IF (cof < 1.) THEN
119          IF (rlamda(imx)*cos(rlatv(j))<1.) jfiltnv = j            IF (rlamda(imx)*cos(rlatv(j)) < 1.) jfiltnv = j
120       END IF         END IF
121    
122       cof = cos(rlatv(jjm-j+1))/colat0         cof = cos(rlatv(jjm-j+1))/colat0
123       IF (cof<1.) THEN         IF (cof < 1.) THEN
124          IF (rlamda(imx)*cos(rlatv(jjm-j+1))<1.) jfiltsv = jjm - j + 1            IF (rlamda(imx)*cos(rlatv(jjm-j+1)) < 1.) jfiltsv = jjm - j + 1
125       END IF         END IF
126    END DO      END DO
127    
128        IF (jfiltnu <= 0) jfiltnu = 1
129    IF (jfiltnu<=0) jfiltnu = 1      IF (jfiltnu > jjm/2+1) THEN
130    IF (jfiltnu>jjm/2+1) THEN         PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu
131       PRINT *, ' jfiltnu en dehors des valeurs acceptables ', jfiltnu         STOP 1
132       STOP 1      END IF
133    END IF  
134        IF (jfiltsu <= 0) jfiltsu = 1
135    IF (jfiltsu<=0) jfiltsu = 1      IF (jfiltsu > jjm + 1) THEN
136    IF (jfiltsu>jjm+1) THEN         PRINT *, 'jfiltsu en dehors des valeurs acceptables ', jfiltsu
137       PRINT *, ' jfiltsu en dehors des valeurs acceptables ', jfiltsu         STOP 1
138       STOP 1      END IF
139    END IF  
140        IF (jfiltnv <= 0) jfiltnv = 1
141    IF (jfiltnv<=0) jfiltnv = 1      IF (jfiltnv > jjm/2) THEN
142    IF (jfiltnv>jjm/2) THEN         PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv
143       PRINT *, ' jfiltnv en dehors des valeurs acceptables ', jfiltnv         STOP 1
144       STOP 1      END IF
145    END IF  
146        IF (jfiltsv <= 0) jfiltsv = 1
147    IF (jfiltsv<=0) jfiltsv = 1      IF (jfiltsv > jjm) THEN
148    IF (jfiltsv>jjm) THEN         PRINT *, 'jfiltsv en dehors des valeurs acceptables ', jfiltsv
149       PRINT *, ' jfiltsv en dehors des valeurs acceptables ', jfiltsv         STOP 1
150       STOP 1      END IF
151    END IF  
152        PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
153    PRINT *, ' jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &           jfiltsu
154         jfiltsu  
155        ! Determination de coefilu, coefilv, n=modfrstu, modfrstv
156    
157    !   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....            DO j = 1, jjm
158    !................................................................               modfrstu(j) = iim
159           modfrstv(j) = iim
160        END DO
161    DO  j = 1, jjm  
162       modfrstu(j) = iim      DO j = 2, jfiltnu
163       modfrstv(j) = iim         DO k = 2, modemax
164    END DO            cof = rlamda(k) * cos(rlatu(j))
165              IF (cof < 1.) exit
166    DO  j = 2, jfiltnu         end DO
167       DO  k = 2, modemax         if (k == modemax + 1) cycle
168          cof = rlamda(k)*cos(rlatu(j))         modfrstu(j) = k
169          IF (cof<1.) GO TO 82  
170       end DO         kf = modfrstu(j)
171       cycle         DO k = kf, modemax
172  82   modfrstu(j) = k            cof = rlamda(k)*cos(rlatu(j))
173              coefilu(k, j) = cof - 1.
174       kf = modfrstu(j)            coefilu2(k, j) = cof*cof - 1.
175       DO  k = kf, modemax         end DO
176          cof = rlamda(k)*cos(rlatu(j))      END DO
177          coefilu(k,j) = cof - 1.  
178          coefilu2(k,j) = cof*cof - 1.      DO j = 1, jfiltnv
179       end DO         DO k = 2, modemax
180    END DO            cof = rlamda(k)*cos(rlatv(j))
181              IF (cof < 1.) exit
182           end DO
183    DO j = 1, jfiltnv         if (k == modemax + 1) cycle
184       DO  k = 2, modemax         modfrstv(j) = k
185          cof = rlamda(k)*cos(rlatv(j))  
186          IF (cof<1.) GO TO 87         kf = modfrstv(j)
187       end DO         DO k = kf, modemax
188       cycle            cof = rlamda(k)*cos(rlatv(j))
189  87   modfrstv(j) = k            coefilv(k, j) = cof - 1.
190              coefilv2(k, j) = cof*cof - 1.
191       kf = modfrstv(j)         end DO
192       DO  k = kf, modemax      end DO
193          cof = rlamda(k)*cos(rlatv(j))  
194          coefilv(k,j) = cof - 1.      DO j = jfiltsu, jjm
195          coefilv2(k,j) = cof*cof - 1.         DO k = 2, modemax
196       end DO            cof = rlamda(k)*cos(rlatu(j))
197    end DO            IF (cof < 1.) exit
198           end DO
199    DO  j = jfiltsu, jjm         if (k == modemax + 1) cycle
200       DO  k = 2, modemax         modfrstu(j) = k
201          cof = rlamda(k)*cos(rlatu(j))  
202          IF (cof<1.) GO TO 92         kf = modfrstu(j)
203       end DO         DO k = kf, modemax
204       cycle            cof = rlamda(k)*cos(rlatu(j))
205  92   modfrstu(j) = k            coefilu(k, j) = cof - 1.
206              coefilu2(k, j) = cof*cof - 1.
207       kf = modfrstu(j)         end DO
208       DO  k = kf, modemax      end DO
209          cof = rlamda(k)*cos(rlatu(j))  
210          coefilu(k,j) = cof - 1.      DO j = jfiltsv, jjm
211          coefilu2(k,j) = cof*cof - 1.         DO k = 2, modemax
212       end DO            cof = rlamda(k)*cos(rlatv(j))
213    end DO            IF (cof < 1.) exit
214           end DO
215    DO  j = jfiltsv, jjm         if (k == modemax + 1) cycle
216       DO  k = 2, modemax         modfrstv(j) = k
217          cof = rlamda(k)*cos(rlatv(j))  
218          IF (cof<1.) GO TO 97         kf = modfrstv(j)
219       end DO         DO k = kf, modemax
220       cycle            cof = rlamda(k)*cos(rlatv(j))
221  97   modfrstv(j) = k            coefilv(k, j) = cof - 1.
222              coefilv2(k, j) = cof*cof - 1.
223       kf = modfrstv(j)         end DO
224       DO  k = kf, modemax      END DO
225          cof = rlamda(k)*cos(rlatv(j))  
226          coefilv(k,j) = cof - 1.      IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN
227          coefilv2(k,j) = cof*cof - 1.         IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
228       end DO         IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
229    END DO  
230           PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
231                jfiltsu
232    IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN      END IF
233    
234       IF (jfiltnv==jfiltsv) jfiltsv = 1 + jfiltnv      PRINT *, 'Modes premiers v '
235       IF (jfiltnu==jfiltsu) jfiltsu = 1 + jfiltnu      PRINT 334, modfrstv
236        PRINT *, 'Modes premiers u '
237       PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &      PRINT 334, modfrstu
238            jfiltsu  
239    END IF      allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
240        allocate(matricevn(iim, iim, jfiltnv))
241    PRINT *, '   Modes premiers  v  '      allocate(matricevs(iim, iim, jfiltsv:jjm))
242    PRINT 334, modfrstv      allocate(matriceus(iim, iim, jfiltsu:jjm), matrinvs(iim, iim, jfiltsu:jjm))
243    PRINT *, '   Modes premiers  u  '  
244    PRINT 334, modfrstu      ! Calcul de la matrice filtre 'matriceu' pour les champs situes
245        ! sur la grille scalaire
246    
247    IF (nfilun<jfiltnu) THEN      DO j = 2, jfiltnu
248       PRINT *, ' le parametre nfilun utilise pour la matrice ', &         DO i = 1, iim
249            ' matriceun  est trop petit ! '            IF (i < modfrstu(j)) then
250       PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnu               coff = 0.
251       PRINT *, 'Pour information, nfilun,nfilus,nfilvn,nfilvs ', &            else
252            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &               coff = coefilu(i, j)
253            jfiltnv, jjm - jfiltsv + 1            end IF
254       STOP 1            eignft(i, :) = eignfnv(:, i)*coff
255    END IF         END DO
256    IF (nfilun>jfiltnu+2) THEN         matriceun(:, :, j) = matmul(eignfnv, eignft)
257       PRINT *, ' le parametre nfilun utilise pour la matrice ', &      END DO
258            ' matriceun est trop grand ! Gachis de memoire ! '  
259       PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnu      DO j = jfiltsu, jjm
260       PRINT *, 'Pour information, nfilun,nfilus,nfilvn,nfilvs ', &         DO i = 1, iim
261            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &            IF (i < modfrstu(j)) then
262            jfiltnv, jjm - jfiltsv + 1               coff = 0.
263    END IF            else
264    IF (nfilus<jjm-jfiltsu+1) THEN               coff = coefilu(i, j)
265       PRINT *, ' le parametre nfilus utilise pour la matrice ', &            end IF
266            ' matriceus  est trop petit !  '            eignft(i, :) = eignfnv(:, i) * coff
267       PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &         END DO
268            jjm - jfiltsu + 1         matriceus(:, :, j) = matmul(eignfnv, eignft)
269       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &      END DO
270            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &  
271            jfiltnv, jjm - jfiltsv + 1      ! Calcul de la matrice filtre 'matricev' pour les champs situes
272       STOP 1      ! sur la grille de V ou de Z
273    END IF  
274    IF (nfilus>jjm-jfiltsu+3) THEN      DO j = 1, jfiltnv
275       PRINT *, ' le parametre nfilus utilise pour la matrice ', &         DO i = 1, iim
276            ' matriceus  est trop grand ! '            IF (i < modfrstv(j)) then
277       PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &               coff = 0.
278            jjm - jfiltsu + 1            else
279       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &               coff = coefilv(i, j)
280            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &            end IF
281            jfiltnv, jjm - jfiltsv + 1            eignft(i, :) = eignfnu(:, i)*coff
282    END IF         END DO
283    IF (nfilvn<jfiltnv) THEN         matricevn(:, :, j) = matmul(eignfnu, eignft)
284       PRINT *, ' le parametre nfilvn utilise pour la matrice ', &      END DO
285            ' matricevn  est trop petit ! '  
286       PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnv      DO j = jfiltsv, jjm
287       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &         DO i = 1, iim
288            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &            IF (i < modfrstv(j)) then
289            jfiltnv, jjm - jfiltsv + 1               coff = 0.
290       STOP 1            else
291    END IF               coff = coefilv(i, j)
292    IF (nfilvn>jfiltnv+2) THEN            end IF
293       PRINT *, ' le parametre nfilvn utilise pour la matrice ', &            eignft(i, :) = eignfnu(:, i)*coff
294            ' matricevn est trop grand !  Gachis de memoire ! '         END DO
295       PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnv         matricevs(:, :, j) = matmul(eignfnu, eignft)
296       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &      END DO
297            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &  
298            jfiltnv, jjm - jfiltsv + 1      ! Calcul de la matrice filtre 'matrinv' pour les champs situes
299    END IF      ! sur la grille scalaire , pour le filtre inverse
300    IF (nfilvs<jjm-jfiltsv+1) THEN  
301       PRINT *, ' le parametre nfilvs utilise pour la matrice ', &      DO j = 2, jfiltnu
302            ' matricevs  est trop petit !  Le changer dans parafilt.h '         DO i = 1, iim
303       PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &            IF (i < modfrstu(j)) then
304            jjm - jfiltsv + 1               coff = 0.
305       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &            else
306            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &               coff = coefilu(i, j)/(1.+coefilu(i, j))
307            jfiltnv, jjm - jfiltsv + 1            end IF
308       STOP 1            eignft(i, :) = eignfnv(:, i)*coff
309    END IF         END DO
310    IF (nfilvs>jjm-jfiltsv+3) THEN         matrinvn(:, :, j) = matmul(eignfnv, eignft)
311       PRINT *, ' le parametre nfilvs utilise pour la matrice ', &      END DO
312            ' matricevs  est trop grand ! Gachis de memoire ! '  
313       PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &      DO j = jfiltsu, jjm
314            jjm - jfiltsv + 1         DO i = 1, iim
315       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &            IF (i < modfrstu(j)) then
316            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &               coff = 0.
317            jfiltnv, jjm - jfiltsv + 1            else
318    END IF               coff = coefilu(i, j)/(1.+coefilu(i, j))
319              end IF
320    !   ... Calcul de la matrice filtre 'matriceu'  pour les champs situes              eignft(i, :) = eignfnv(:, i)*coff
321    !                       sur la grille scalaire                 ........         END DO
322           matrinvs(:, :, j) = matmul(eignfnv, eignft)
323    DO j = 2, jfiltnu      END DO
   
      DO i = 1, iim  
         coff = coefilu(i,j)  
         IF (i<modfrstu(j)) coff = 0.  
         DO k = 1, iim  
            eignft(i,k) = eignfnv(k,i)*coff  
         END DO  
      END DO  
      DO k = 1, iim  
         DO i = 1, iim  
            matriceun(i,k,j) = 0.0  
            DO ii = 1, iim  
               matriceun(i,k,j) = matriceun(i,k,j) + eignfnv(i,ii)*eignft(ii,k)  
            END DO  
         END DO  
      END DO  
   
   END DO  
   
   DO j = jfiltsu, jjm  
   
      DO i = 1, iim  
         coff = coefilu(i,j)  
         IF (i<modfrstu(j)) coff = 0.  
         DO k = 1, iim  
            eignft(i,k) = eignfnv(k,i)*coff  
         END DO  
      END DO  
      DO k = 1, iim  
         DO i = 1, iim  
            matriceus(i,k,j-jfiltsu+1) = 0.0  
            DO ii = 1, iim  
               matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) + &  
                    eignfnv(i,ii)*eignft(ii,k)  
            END DO  
         END DO  
      END DO  
   
   END DO  
   
   !   ...................................................................  
   
   !   ... Calcul de la matrice filtre 'matricev'  pour les champs situes    
   !                       sur la grille   de V ou de Z           ........  
   !   ...................................................................  
   
   DO j = 1, jfiltnv  
   
      DO i = 1, iim  
         coff = coefilv(i,j)  
         IF (i<modfrstv(j)) coff = 0.  
         DO k = 1, iim  
            eignft(i,k) = eignfnu(k,i)*coff  
         END DO  
      END DO  
      DO k = 1, iim  
         DO i = 1, iim  
            matricevn(i,k,j) = 0.0  
            DO ii = 1, iim  
               matricevn(i,k,j) = matricevn(i,k,j) + eignfnu(i,ii)*eignft(ii,k)  
            END DO  
         END DO  
      END DO  
   
   END DO  
   
   DO j = jfiltsv, jjm  
   
      DO i = 1, iim  
         coff = coefilv(i,j)  
         IF (i<modfrstv(j)) coff = 0.  
         DO k = 1, iim  
            eignft(i,k) = eignfnu(k,i)*coff  
         END DO  
      END DO  
      DO k = 1, iim  
         DO i = 1, iim  
            matricevs(i,k,j-jfiltsv+1) = 0.0  
            DO ii = 1, iim  
               matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) + &  
                    eignfnu(i,ii)*eignft(ii,k)  
            END DO  
         END DO  
      END DO  
   
   END DO  
   
   !   ...................................................................  
   
   !   ... Calcul de la matrice filtre 'matrinv'  pour les champs situes    
   !              sur la grille scalaire , pour le filtre inverse ........  
   !   ...................................................................  
   
   DO j = 2, jfiltnu  
   
      DO i = 1, iim  
         coff = coefilu(i,j)/(1.+coefilu(i,j))  
         IF (i<modfrstu(j)) coff = 0.  
         DO k = 1, iim  
            eignft(i,k) = eignfnv(k,i)*coff  
         END DO  
      END DO  
      DO k = 1, iim  
         DO i = 1, iim  
            matrinvn(i,k,j) = 0.0  
            DO ii = 1, iim  
               matrinvn(i,k,j) = matrinvn(i,k,j) + eignfnv(i,ii)*eignft(ii,k)  
            END DO  
         END DO  
      END DO  
   
   END DO  
   
   DO j = jfiltsu, jjm  
   
      DO i = 1, iim  
         coff = coefilu(i,j)/(1.+coefilu(i,j))  
         IF (i<modfrstu(j)) coff = 0.  
         DO k = 1, iim  
            eignft(i,k) = eignfnv(k,i)*coff  
         END DO  
      END DO  
      DO k = 1, iim  
         DO i = 1, iim  
            matrinvs(i,k,j-jfiltsu+1) = 0.0  
            DO ii = 1, iim  
               matrinvs(i,k,j-jfiltsu+1) = matrinvs(i,k,j-jfiltsu+1) + &  
                    eignfnv(i,ii)*eignft(ii,k)  
            END DO  
         END DO  
      END DO  
324    
325    END DO  334 FORMAT (1X, 24I3)
326    
327  334 FORMAT (1X,24I3)    END SUBROUTINE inifilr
 755 FORMAT (1X,6F10.3,I3)  
328    
329  END SUBROUTINE inifilr  end module inifilr_m
Legend:
Removed from v.32  
changed lines
  Added in v.139
  ViewVC Help
Powered by ViewVC 1.1.21