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

Diff of /trunk/filtrez/inifilr.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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