/[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

revision 134 by guez, Wed Apr 29 15:47:56 2015 UTC revision 156 by guez, Thu Jul 16 17:39:10 2015 UTC
# Line 1  Line 1 
1  module inifilr_m  module inifilr_m
2    
   use dimens_m, only: iim  
   
3    IMPLICIT NONE    IMPLICIT NONE
4    
5    INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv    INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv
6    INTEGER, PARAMETER:: nfilun=3, nfilus=2, nfilvn=2, nfilvs=2    ! jfiltn index of the last scalar line filtered in NH
7      ! jfilts index of the first line filtered in SH
8    
9      ! North:
10    
11      real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :)
12      ! (iim, iim, 2:jfiltnu)
13    
14      real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)
15    
16    real matriceun(iim,iim,nfilun), matriceus(iim,iim,nfilus)    ! South:
   real matricevn(iim,iim,nfilvn), matricevs(iim,iim,nfilvs)  
   real matrinvn(iim,iim,nfilun), matrinvs(iim,iim,nfilus)  
17    
18    private iim, nfilun, nfilus, nfilvn, nfilvs    real, allocatable:: matriceus(:, :, :), matrinvs(:, :, :)
19      ! (iim, iim, jfiltsu:jjm)
20    
21      real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
22    
23  contains  contains
24    
# Line 20  contains Line 27  contains
27      ! From filtrez/inifilr.F, version 1.1.1.1 2004/05/19 12:53:09      ! From filtrez/inifilr.F, version 1.1.1.1 2004/05/19 12:53:09
28      ! H. Upadhyaya, O. Sharma      ! H. Upadhyaya, O. Sharma
29    
30      ! This routine computes the eigenfunctions of the laplacian on the      ! This routine computes the eigenvectors of the laplacian on the
31      ! stretched grid, and the filtering coefficients.      ! stretched grid, and the filtering coefficients. The modes are
32      ! We designate:      ! filtered from modfrst to iim.
33      ! eignfn eigenfunctions of the discrete laplacian  
     ! eigenvl eigenvalues  
     ! jfiltn index of 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 coefils, ONLY : coefilu, coefilu2, coefilv, coefilv2, eignfnu, &  
          eignfnv, modfrstu, modfrstv  
     USE comgeom, ONLY : rlatu, rlatv, xprimu  
34      USE dimens_m, ONLY : iim, jjm      USE dimens_m, ONLY : iim, jjm
35        USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
36      use inifgn_m, only: inifgn      use inifgn_m, only: inifgn
37        use jumble, only: new_unit
38      use nr_util, only: pi      use nr_util, only: pi
     USE serre, ONLY : grossismx  
39    
40      ! Local:      ! Local:
41      REAL dlatu(jjm)      REAL dlatu(jjm)
42      REAL rlamda(2: iim), eignvl(iim)      REAL rlamda(2: iim)
43        real eignvl(iim) ! eigenvalues sorted in descending order
44      REAL lamdamax, cof      REAL cof
45      INTEGER i, j, modemax, imx, k, kf      INTEGER i, j, k, unit
46      REAL dymin, colat0      REAL colat0 ! > 0
47      REAL eignft(iim, iim), coff      REAL eignft(iim, iim), coff
48    
49        real eignfnu(iim, iim), eignfnv(iim, iim)
50        ! eigenvectors of the discrete laplacian
51    
52        ! Filtering coefficients (lamda_max * cos(rlat) / lamda):
53        real coefilu(iim, jjm), coefilv(iim, jjm)
54        real coefilu2(iim, jjm), coefilv2(iim, jjm)
55    
56        ! Index of the mode from where modes are filtered:
57        integer, allocatable:: modfrstnu(:), modfrstsu(:)
58        integer, allocatable:: modfrstnv(:), modfrstsv(:)
59    
60      !-----------------------------------------------------------      !-----------------------------------------------------------
61    
62      print *, "Call sequence information: inifilr"      print *, "Call sequence information: inifilr"
63    
64      CALL inifgn(eignvl)      CALL inifgn(eignvl, eignfnu, eignfnv)
65    
66      PRINT *, 'EIGNVL '      ! compute eigenvalues and eigenvectors
     PRINT "(1X, 5E13.6)", eignvl  
   
     ! compute eigenvalues and eigenfunctions  
67      ! compute the filtering coefficients for scalar lines and      ! compute the filtering coefficients for scalar lines and
68      ! meridional wind v-lines      ! meridional wind v-lines
69      ! we filter all those latitude lines where coefil < 1      ! we filter all those latitude lines where coefil < 1
# Line 69  contains Line 72  contains
72      ! is set equal to zero for the regular grid case      ! is set equal to zero for the regular grid case
73    
74      ! Calcul de colat0      ! Calcul de colat0
75        forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
76      DO j = 1, jjm      colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
        dlatu(j) = rlatu(j) - rlatu(j+1)  
     END DO  
   
     dymin = dlatu(1)  
     DO j = 2, jjm  
        dymin = min(dymin, dlatu(j))  
     END DO  
   
     colat0 = min(0.5, dymin / minval(xprimu(:iim)))  
   
77      PRINT *, 'colat0 = ', colat0      PRINT *, 'colat0 = ', colat0
78    
79      lamdamax = iim / (pi * colat0 / grossismx)      rlamda = iim / (pi * colat0 / grossismx) / sqrt(abs(eignvl(2: iim)))
     rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))  
   
     DO j = 1, jjm  
        DO i = 1, iim  
           coefilu(i, j) = 0.  
           coefilv(i, j) = 0.  
           coefilu2(i, j) = 0.  
           coefilv2(i, j) = 0.  
        end DO  
     END DO  
   
     ! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv  
   
     modemax = iim  
     imx = iim  
   
     PRINT *, 'TRUNCATION AT ', imx  
   
     DO j = 2, jjm / 2 + 1  
        IF (cos(rlatu(j)) / colat0 < 1. &  
             .and. rlamda(imx) * cos(rlatu(j)) < 1.) jfiltnu = j  
   
        IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &  
             .and. rlamda(imx) * cos(rlatu(jjm - j + 2)) < 1.) &  
             jfiltsu = jjm - j + 2  
     END DO  
   
     DO j = 1, jjm/2  
        cof = cos(rlatv(j))/colat0  
        IF (cof < 1.) THEN  
           IF (rlamda(imx)*cos(rlatv(j)) < 1.) jfiltnv = j  
        END IF  
   
        cof = cos(rlatv(jjm-j+1))/colat0  
        IF (cof < 1.) THEN  
           IF (rlamda(imx)*cos(rlatv(jjm-j+1)) < 1.) jfiltsv = jjm - j + 1  
        END IF  
     END DO  
   
     IF (jfiltnu <= 0) jfiltnu = 1  
     IF (jfiltnu > jjm/2+1) THEN  
        PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu  
        STOP 1  
     END IF  
80    
81      IF (jfiltsu <= 0) jfiltsu = 1      ! Determination de jfiltnu, jfiltsu, jfiltnv, jfiltsv
     IF (jfiltsu > jjm + 1) THEN  
        PRINT *, 'jfiltsu en dehors des valeurs acceptables ', jfiltsu  
        STOP 1  
     END IF  
82    
83      IF (jfiltnv <= 0) jfiltnv = 1      jfiltnu = (jjm + 1) / 2
84      IF (jfiltnv > jjm/2) THEN      do while (cos(rlatu(jfiltnu)) >= colat0 &
85         PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv           .or. rlamda(iim) * cos(rlatu(jfiltnu)) >= 1.)
86           jfiltnu = jfiltnu - 1
87        end do
88    
89        jfiltsu = jjm / 2 + 2
90        do while (cos(rlatu(jfiltsu)) >= colat0 &
91             .or. rlamda(iim) * cos(rlatu(jfiltsu)) >= 1.)
92           jfiltsu = jfiltsu + 1
93        end do
94    
95        jfiltnv = jjm / 2
96        do while ((cos(rlatv(jfiltnv)) >= colat0 &
97             .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) .and. jfiltnv >= 2)
98           jfiltnv = jfiltnv - 1
99        end do
100    
101        if (cos(rlatv(jfiltnv)) >= colat0 &
102             .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) then
103           ! {jfiltnv == 1}
104           PRINT *, 'Could not find jfiltnv.'
105         STOP 1         STOP 1
106      END IF      END IF
107    
108      IF (jfiltsv <= 0) jfiltsv = 1      jfiltsv = (jjm + 1)/ 2 + 1
109      IF (jfiltsv > jjm) THEN      do while ((cos(rlatv(jfiltsv)) >= colat0 &
110         PRINT *, 'jfiltsv en dehors des valeurs acceptables ', jfiltsv           .or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) .and. jfiltsv <= jjm - 1)
111           jfiltsv = jfiltsv + 1
112        end do
113    
114        IF (cos(rlatv(jfiltsv)) >= colat0 &
115             .or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) THEN
116           ! {jfiltsv == jjm}
117           PRINT *, 'Could not find jfiltsv.'
118         STOP 1         STOP 1
119      END IF      END IF
120    
121      PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &      PRINT *, 'jfiltnu =', jfiltnu
122           jfiltsu      PRINT *, 'jfiltsu =', jfiltsu
123        PRINT *, 'jfiltnv =', jfiltnv
124      ! Determination de coefilu, coefilv, n=modfrstu, modfrstv      PRINT *, 'jfiltsv =', jfiltsv
125    
126      DO j = 1, jjm      ! Determination de coefilu, coefilv, modfrst[ns][uv]:
127         modfrstu(j) = iim  
128         modfrstv(j) = iim      allocate(modfrstnu(2:jfiltnu), modfrstsu(jfiltsu:jjm))
129      END DO      allocate(modfrstnv(jfiltnv), modfrstsv(jfiltsv:jjm))
130        coefilu = 0.
131        coefilv = 0.
132        coefilu2 = 0.
133        coefilv2 = 0.
134    
135      DO j = 2, jfiltnu      DO j = 2, jfiltnu
136         DO k = 2, modemax         modfrstnu(j) = 2
137            cof = rlamda(k) * cos(rlatu(j))         do while (rlamda(modfrstnu(j)) * cos(rlatu(j)) >= 1. &
138            IF (cof < 1.) exit              .and. modfrstnu(j) <= iim - 1)
139         end DO            modfrstnu(j) = modfrstnu(j) + 1
140         if (k == modemax + 1) cycle         end do
141         modfrstu(j) = k  
142           if (rlamda(modfrstnu(j)) * cos(rlatu(j)) < 1.) then
143         kf = modfrstu(j)            DO k = modfrstnu(j), iim
144         DO k = kf, modemax               cof = rlamda(k) * cos(rlatu(j))
145            cof = rlamda(k)*cos(rlatu(j))               coefilu(k, j) = cof - 1.
146            coefilu(k, j) = cof - 1.               coefilu2(k, j) = cof**2 - 1.
147            coefilu2(k, j) = cof*cof - 1.            end DO
148         end DO         end if
149      END DO      END DO
150    
151      DO j = 1, jfiltnv      DO j = 1, jfiltnv
152         DO k = 2, modemax         modfrstnv(j) = 2
153            cof = rlamda(k)*cos(rlatv(j))         do while (rlamda(modfrstnv(j)) * cos(rlatv(j)) >= 1. &
154            IF (cof < 1.) exit              .and. modfrstnv(j) <= iim - 1)
155         end DO            modfrstnv(j) = modfrstnv(j) + 1
156         if (k == modemax + 1) cycle         end do
157         modfrstv(j) = k  
158           if (rlamda(modfrstnv(j)) * cos(rlatv(j)) < 1.) then
159         kf = modfrstv(j)            DO k = modfrstnv(j), iim
160         DO k = kf, modemax               cof = rlamda(k) * cos(rlatv(j))
161            cof = rlamda(k)*cos(rlatv(j))               coefilv(k, j) = cof - 1.
162            coefilv(k, j) = cof - 1.               coefilv2(k, j) = cof**2 - 1.
163            coefilv2(k, j) = cof*cof - 1.            end DO
164         end DO         end if
165      end DO      end DO
166    
167      DO j = jfiltsu, jjm      DO j = jfiltsu, jjm
168         DO k = 2, modemax         modfrstsu(j) = 2
169            cof = rlamda(k)*cos(rlatu(j))         do while (rlamda(modfrstsu(j)) * cos(rlatu(j)) >= 1. &
170            IF (cof < 1.) exit              .and. modfrstsu(j) <= iim - 1)
171         end DO            modfrstsu(j) = modfrstsu(j) + 1
172         if (k == modemax + 1) cycle         end do
173         modfrstu(j) = k  
174           if (rlamda(modfrstsu(j)) * cos(rlatu(j)) < 1.) then
175         kf = modfrstu(j)            DO k = modfrstsu(j), iim
176         DO k = kf, modemax               cof = rlamda(k) * cos(rlatu(j))
177            cof = rlamda(k)*cos(rlatu(j))               coefilu(k, j) = cof - 1.
178            coefilu(k, j) = cof - 1.               coefilu2(k, j) = cof**2 - 1.
179            coefilu2(k, j) = cof*cof - 1.            end DO
180         end DO         end if
181      end DO      end DO
182    
183      DO j = jfiltsv, jjm      DO j = jfiltsv, jjm
184         DO k = 2, modemax         modfrstsv(j) = 2
185            cof = rlamda(k)*cos(rlatv(j))         do while (rlamda(modfrstsv(j)) * cos(rlatv(j)) >= 1. &
186            IF (cof < 1.) exit              .and. modfrstsv(j) <= iim - 1)
187         end DO            modfrstsv(j) = modfrstsv(j) + 1
188         if (k == modemax + 1) cycle         end do
189         modfrstv(j) = k  
190           if (rlamda(modfrstsv(j)) * cos(rlatv(j)) < 1.) then
191         kf = modfrstv(j)            DO k = modfrstsv(j), iim
192         DO k = kf, modemax               cof = rlamda(k) * cos(rlatv(j))
193            cof = rlamda(k)*cos(rlatv(j))               coefilv(k, j) = cof - 1.
194            coefilv(k, j) = cof - 1.               coefilv2(k, j) = cof**2 - 1.
195            coefilv2(k, j) = cof*cof - 1.            end DO
196         end DO         end if
197      END DO      END DO
198    
199      IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN      call new_unit(unit)
200         IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv      open(unit, file = "inifilr_out.txt", status = "replace", action = "write")
201         IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu      write(unit, fmt = *) '"EIGNVL"', eignvl
202        write(unit, fmt = *) '"modfrstnu"', modfrstnu
203         PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &      write(unit, fmt = *) '"modfrstsu"', modfrstsu
204              jfiltsu      write(unit, fmt = *) '"modfrstnv"', modfrstnv
205      END IF      write(unit, fmt = *) '"modfrstsv"', modfrstsv
206        close(unit)
207      PRINT *, 'Modes premiers v '  
208      PRINT 334, modfrstv      allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
209      PRINT *, 'Modes premiers u '      allocate(matricevn(iim, iim, jfiltnv))
210      PRINT 334, modfrstu      allocate(matricevs(iim, iim, jfiltsv:jjm))
211        allocate(matriceus(iim, iim, jfiltsu:jjm), matrinvs(iim, iim, jfiltsu:jjm))
     IF (nfilun < jfiltnu) THEN  
        PRINT *, 'le parametre nfilun utilise pour la matrice ', &  
             'matriceun est trop petit ! '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnu  
        PRINT *, 'Pour information, nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', jfiltnu, &  
             jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
        STOP 1  
     END IF  
     IF (nfilun > jfiltnu+2) THEN  
        PRINT *, 'le parametre nfilun utilise pour la matrice ', &  
             'matriceun est trop grand ! Gachis de memoire ! '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnu  
        PRINT *, 'Pour information, nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
     END IF  
     IF (nfilus < jjm-jfiltsu+1) THEN  
        PRINT *, 'le parametre nfilus utilise pour la matrice ', &  
             'matriceus est trop petit ! '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', &  
             jjm - jfiltsu + 1  
        PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
        STOP 1  
     END IF  
     IF (nfilus > jjm-jfiltsu+3) THEN  
        PRINT *, 'le parametre nfilus utilise pour la matrice ', &  
             'matriceus est trop grand ! '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', &  
             jjm - jfiltsu + 1  
        PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
     END IF  
     IF (nfilvn < jfiltnv) THEN  
        PRINT *, 'le parametre nfilvn utilise pour la matrice ', &  
             'matricevn est trop petit ! '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnv  
        PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
        STOP 1  
     END IF  
     IF (nfilvn > jfiltnv+2) THEN  
        PRINT *, 'le parametre nfilvn utilise pour la matrice ', &  
             'matricevn est trop grand ! Gachis de memoire ! '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnv  
        PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
     END IF  
     IF (nfilvs < jjm-jfiltsv+1) THEN  
        PRINT *, 'le parametre nfilvs utilise pour la matrice ', &  
             'matricevs est trop petit ! Le changer dans parafilt.h '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', &  
             jjm - jfiltsv + 1  
        PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
        STOP 1  
     END IF  
     IF (nfilvs > jjm-jfiltsv+3) THEN  
        PRINT *, 'le parametre nfilvs utilise pour la matrice ', &  
             'matricevs est trop grand ! Gachis de memoire ! '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', &  
             jjm - jfiltsv + 1  
        PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
     END IF  
212    
213      ! Calcul de la matrice filtre 'matriceu' pour les champs situes      ! Calcul de la matrice filtre 'matriceu' pour les champs situes
214      ! sur la grille scalaire      ! sur la grille scalaire
215    
216      DO j = 2, jfiltnu      DO j = 2, jfiltnu
217         DO i = 1, iim         DO i = 1, iim
218            IF (i < modfrstu(j)) then            IF (i < modfrstnu(j)) then
219               coff = 0.               coff = 0.
220            else            else
221               coff = coefilu(i, j)               coff = coefilu(i, j)
222            end IF            end IF
223            eignft(i, :) = eignfnv(:, i)*coff            eignft(i, :) = eignfnv(:, i) * coff
224         END DO         END DO
225         matriceun(:, :, j) = matmul(eignfnv, eignft)         matriceun(:, :, j) = matmul(eignfnv, eignft)
226      END DO      END DO
227    
228      DO j = jfiltsu, jjm      DO j = jfiltsu, jjm
229         DO i = 1, iim         DO i = 1, iim
230            IF (i < modfrstu(j)) then            IF (i < modfrstsu(j)) then
231               coff = 0.               coff = 0.
232            else            else
233               coff = coefilu(i, j)               coff = coefilu(i, j)
234            end IF            end IF
235            eignft(i, :) = eignfnv(:, i) * coff            eignft(i, :) = eignfnv(:, i) * coff
236         END DO         END DO
237         matriceus(:, :, j - jfiltsu + 1) = matmul(eignfnv, eignft)         matriceus(:, :, j) = matmul(eignfnv, eignft)
238      END DO      END DO
239    
240      ! Calcul de la matrice filtre 'matricev' pour les champs situes      ! Calcul de la matrice filtre 'matricev' pour les champs situes
# Line 339  contains Line 242  contains
242    
243      DO j = 1, jfiltnv      DO j = 1, jfiltnv
244         DO i = 1, iim         DO i = 1, iim
245            IF (i < modfrstv(j)) then            IF (i < modfrstnv(j)) then
246               coff = 0.               coff = 0.
247            else            else
248               coff = coefilv(i, j)               coff = coefilv(i, j)
249            end IF            end IF
250            eignft(i, :) = eignfnu(:, i)*coff            eignft(i, :) = eignfnu(:, i) * coff
251         END DO         END DO
252         matricevn(:, :, j) = matmul(eignfnu, eignft)         matricevn(:, :, j) = matmul(eignfnu, eignft)
253      END DO      END DO
254    
255      DO j = jfiltsv, jjm      DO j = jfiltsv, jjm
256         DO i = 1, iim         DO i = 1, iim
257            IF (i < modfrstv(j)) then            IF (i < modfrstsv(j)) then
258               coff = 0.               coff = 0.
259            else            else
260               coff = coefilv(i, j)               coff = coefilv(i, j)
261            end IF            end IF
262            eignft(i, :) = eignfnu(:, i)*coff            eignft(i, :) = eignfnu(:, i) * coff
263         END DO         END DO
264         matricevs(:, :, j-jfiltsv+1) = matmul(eignfnu, eignft)         matricevs(:, :, j) = matmul(eignfnu, eignft)
265      END DO      END DO
266    
267      ! Calcul de la matrice filtre 'matrinv' pour les champs situes      ! Calcul de la matrice filtre 'matrinv' pour les champs situes
# Line 366  contains Line 269  contains
269    
270      DO j = 2, jfiltnu      DO j = 2, jfiltnu
271         DO i = 1, iim         DO i = 1, iim
272            IF (i < modfrstu(j)) then            IF (i < modfrstnu(j)) then
273               coff = 0.               coff = 0.
274            else            else
275               coff = coefilu(i, j)/(1.+coefilu(i, j))               coff = coefilu(i, j) / (1. + coefilu(i, j))
276            end IF            end IF
277            eignft(i, :) = eignfnv(:, i)*coff            eignft(i, :) = eignfnv(:, i) * coff
278         END DO         END DO
279         matrinvn(:, :, j) = matmul(eignfnv, eignft)         matrinvn(:, :, j) = matmul(eignfnv, eignft)
280      END DO      END DO
281    
282      DO j = jfiltsu, jjm      DO j = jfiltsu, jjm
283         DO i = 1, iim         DO i = 1, iim
284            IF (i < modfrstu(j)) then            IF (i < modfrstsu(j)) then
285               coff = 0.               coff = 0.
286            else            else
287               coff = coefilu(i, j)/(1.+coefilu(i, j))               coff = coefilu(i, j) / (1. + coefilu(i, j))
288            end IF            end IF
289            eignft(i, :) = eignfnv(:, i)*coff            eignft(i, :) = eignfnv(:, i) * coff
290         END DO         END DO
291         matrinvs(:, :, j-jfiltsu+1) = matmul(eignfnv, eignft)         matrinvs(:, :, j) = matmul(eignfnv, eignft)
292      END DO      END DO
293    
 334 FORMAT (1X, 24I3)  
   
294    END SUBROUTINE inifilr    END SUBROUTINE inifilr
295    
296  end module inifilr_m  end module inifilr_m
Legend:
Removed from v.134  
changed lines
  Added in v.156
  ViewVC Help
Powered by ViewVC 1.1.21