/[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 113 by guez, Thu Sep 18 19:56:46 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.
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 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
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 73  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)  
   
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  
80    
81      modemax = iim      ! Determination de jfiltnu, jfiltsu, jfiltnv, jfiltsv
     imx = iim  
82    
83      PRINT *, 'TRUNCATION AT ', imx      jfiltnu = (jjm + 1) / 2
84        do while (cos(rlatu(jfiltnu)) >= colat0 &
85      DO j = 2, jjm / 2 + 1           .or. rlamda(iim) * cos(rlatu(jfiltnu)) >= 1.)
86         IF (cos(rlatu(j)) / colat0 < 1. &         jfiltnu = jfiltnu - 1
87              .and. rlamda(imx) * cos(rlatu(j)) < 1.) jfiltnu = j      end do
88    
89         IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &      jfiltsu = jjm / 2 + 2
90              .and. rlamda(imx) * cos(rlatu(jjm - j + 2)) < 1.) &      do while (cos(rlatu(jfiltsu)) >= colat0 &
91              jfiltsu = jjm - j + 2           .or. rlamda(iim) * cos(rlatu(jfiltsu)) >= 1.)
92      END DO         jfiltsu = jfiltsu + 1
93        end do
94      DO j = 1, jjm/2  
95         cof = cos(rlatv(j))/colat0      jfiltnv = jjm / 2
96         IF (cof < 1.) THEN      do while ((cos(rlatv(jfiltnv)) >= colat0 &
97            IF (rlamda(imx)*cos(rlatv(j)) < 1.) jfiltnv = j           .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) .and. jfiltnv >= 2)
98         END IF         jfiltnv = jfiltnv - 1
99        end do
100         cof = cos(rlatv(jjm-j+1))/colat0  
101         IF (cof < 1.) THEN      if (cos(rlatv(jfiltnv)) >= colat0 &
102            IF (rlamda(imx)*cos(rlatv(jjm-j+1)) < 1.) jfiltsv = jjm - j + 1           .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) then
103         END IF         ! {jfiltnv == 1}
104      END DO         PRINT *, 'Could not find jfiltnv.'
   
     IF (jfiltnu <= 0) jfiltnu = 1  
     IF (jfiltnu > jjm/2+1) THEN  
        PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu  
105         STOP 1         STOP 1
106      END IF      END IF
107    
108      IF (jfiltsu <= 0) jfiltsu = 1      jfiltsv = (jjm + 1)/ 2 + 1
109      IF (jfiltsu > jjm + 1) THEN      do while ((cos(rlatv(jfiltsv)) >= colat0 &
110         PRINT *, 'jfiltsu en dehors des valeurs acceptables ', jfiltsu           .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      IF (jfiltnv <= 0) jfiltnv = 1      PRINT *, 'jfiltnu =', jfiltnu
122      IF (jfiltnv > jjm/2) THEN      PRINT *, 'jfiltsu =', jfiltsu
123         PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv      PRINT *, 'jfiltnv =', jfiltnv
124         STOP 1      PRINT *, 'jfiltsv =', jfiltsv
125      END IF  
126        ! Determination de coefilu, coefilv, modfrst[ns][uv]:
127      IF (jfiltsv <= 0) jfiltsv = 1  
128      IF (jfiltsv > jjm) THEN      allocate(modfrstnu(2:jfiltnu), modfrstsu(jfiltsu:jjm))
129         PRINT *, 'jfiltsv en dehors des valeurs acceptables ', jfiltsv      allocate(modfrstnv(jfiltnv), modfrstsv(jfiltsv:jjm))
130         STOP 1      coefilu = 0.
131      END IF      coefilv = 0.
132        coefilu2 = 0.
133      PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &      coefilv2 = 0.
          jfiltsu  
   
     ! Determination de coefilu, coefilv, n=modfrstu, modfrstv  
   
     DO j = 1, jjm  
        modfrstu(j) = iim  
        modfrstv(j) = iim  
     END DO  
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 347  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 374  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.113  
changed lines
  Added in v.156
  ViewVC Help
Powered by ViewVC 1.1.21