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

Diff of /trunk/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 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
329    
330  END SUBROUTINE inifilr  end module inifilr_m
Legend:
Removed from v.40  
changed lines
  Added in v.136
  ViewVC Help
Powered by ViewVC 1.1.21