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

Diff of /trunk/filtrez/inifilr.f90

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

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