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