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