/[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 40 by guez, Tue Feb 22 13:49:36 2011 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
328    
329  END SUBROUTINE inifilr  end module inifilr_m
Legend:
Removed from v.40  
changed lines
  Added in v.139
  ViewVC Help
Powered by ViewVC 1.1.21