/[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/libf/filtrez/inifilr.f90 revision 40 by guez, Tue Feb 22 13:49:36 2011 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
299    
300  END SUBROUTINE inifilr  end module inifilr_m
Legend:
Removed from v.40  
changed lines
  Added in v.143
  ViewVC Help
Powered by ViewVC 1.1.21