/[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 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.
     ! 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.  
33    
34      USE dimens_m, ONLY : iim, jjm      USE dimens_m, ONLY : iim, jjm
35      use conf_gcm_m, ONLY : fxyhypb, ysinus      USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
36      USE comgeom, ONLY : rlatu, rlatv, xprimu      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 : alphax  
     USE coefils, ONLY : coefilu, coefilu2, coefilv, coefilv2, eignfnu, &  
          eignfnv, modfrstu, modfrstv  
39    
40      ! Local:      ! Local:
41      REAL dlonu(iim), 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, dxmin, colat0      REAL colat0 ! > 0
47      REAL eignft(iim, iim), coff      REAL eignft(iim, iim), coff
     EXTERNAL inifgn  
48    
49      !-----------------------------------------------------------      real eignfnu(iim, iim), eignfnv(iim, iim)
50        ! eigenvectors of the discrete laplacian
51    
52      print *, "Call sequence information: inifilr"      ! 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      DO i = 1, iim      !-----------------------------------------------------------
        dlonu(i) = xprimu(i)  
     END DO  
61    
62      CALL inifgn(eignvl)      print *, "Call sequence information: inifilr"
63    
64      PRINT *, 'EIGNVL '      CALL inifgn(eignvl, eignfnu, eignfnv)
     PRINT "(1X, 5E13.6)", eignvl  
65    
66      ! compute eigenvalues and eigenfunctions      ! compute eigenvalues and eigenvectors
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 74  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  
   
     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  
   
77      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  
78    
79         IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &      rlamda = iim / (pi * colat0 / grossismx) / sqrt(abs(eignvl(2: iim)))
             .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 360  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 387  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.82  
changed lines
  Added in v.156
  ViewVC Help
Powered by ViewVC 1.1.21