/[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/libf/filtrez/inifilr.f90 revision 32 by guez, Tue Apr 6 17:52:58 2010 UTC trunk/Sources/filtrez/inifilr.f revision 140 by guez, Fri Jun 5 18:58:06 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 nr_util, only: pi
36    
37    !  we filter all those latitude lines where coefil < 1                        ! Local:
38    !  NO FILTERING AT POLES                                                      REAL dlatu(jjm)
39        REAL rlamda(2: iim)
40    !  colat0 is to be used  when alpha (stretching coefficient)                  real eignvl(iim) ! eigenvalues
41    !  is set equal to zero for the regular grid case                            REAL lamdamax, cof
42        INTEGER i, j, k, kf
43    !    .......   Calcul  de  colat0   .........                                REAL dymin, colat0
44    !     .....  colat0 = minimum de ( 0.5, min dy/ min dx )   ...                REAL eignft(iim, iim), coff
45    
46        ! Filtering coefficients (lamda_max * cos(rlat) / lamda):
47    DO  j = 1, jjm      real coefilu(iim, jjm), coefilv(iim, jjm)
48       dlatu(j) = rlatu(j) - rlatu(j+1)      real coefilu2(iim, jjm), coefilv2(iim, jjm)
49    END DO  
50        integer modfrstu(jjm), modfrstv(jjm)
51    dxmin = dlonu(1)      ! index of the mode from where modes are filtered
52    DO i = 2, iim  
53       dxmin = min(dxmin,dlonu(i))      !-----------------------------------------------------------
54    END DO  
55    dymin = dlatu(1)      print *, "Call sequence information: inifilr"
56    DO j = 2, jjm  
57       dymin = min(dymin,dlatu(j))      CALL inifgn(eignvl)
58    END DO  
59        PRINT *, 'EIGNVL '
60        PRINT "(1X, 5E13.6)", eignvl
61    colat0 = min(0.5,dymin/dxmin)  
62        ! compute eigenvalues and eigenfunctions
63    IF ( .NOT. fxyhypb .AND. ysinus) THEN      ! compute the filtering coefficients for scalar lines and
64       colat0 = 0.6      ! meridional wind v-lines
65       !         ...... a revoir  pour  ysinus !   .......                          ! we filter all those latitude lines where coefil < 1
66       alphax = 0.      ! NO FILTERING AT POLES
67    END IF      ! colat0 is to be used when alpha (stretching coefficient)
68        ! is set equal to zero for the regular grid case
69    PRINT 50, colat0, alphax  
70  50 FORMAT (/15X,' Inifilr colat0 alphax ',2E16.7)      ! Calcul de colat0
71    
72    IF (alphax==1.) THEN      DO j = 1, jjm
73       PRINT *, ' Inifilr  alphax doit etre  <  a 1.  Corriger '         dlatu(j) = rlatu(j) - rlatu(j + 1)
74       STOP 1      END DO
75    END IF  
76        dymin = dlatu(1)
77    lamdamax = iim/(pi*colat0*(1.-alphax))      DO j = 2, jjm
78           dymin = min(dymin, dlatu(j))
79    DO  i = 2, iim      END DO
80       rlamda(i) = lamdamax/sqrt(abs(eignvl(i)))  
81    END DO      colat0 = min(0.5, dymin / minval(xprimu(:iim)))
82    
83        PRINT *, 'colat0 = ', colat0
84    DO  j = 1, jjm  
85       DO  i = 1, iim      lamdamax = iim / (pi * colat0 / grossismx)
86          coefilu(i,j) = 0.0      rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))
87          coefilv(i,j) = 0.0  
88          coefilu2(i,j) = 0.0      DO j = 1, jjm
89          coefilv2(i,j) = 0.0         DO i = 1, iim
90       end DO            coefilu(i, j) = 0.
91    END DO            coefilv(i, j) = 0.
92              coefilu2(i, j) = 0.
93              coefilv2(i, j) = 0.
94    !    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....                   end DO
95    !    .........................................................                END DO
96    
97    modemax = iim      ! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
98    
99    !ccc    imx = modemax - 4 * (modemax/iim)                                    PRINT *, 'TRUNCATION AT ', iim
100    
101    imx = iim      DO j = 2, jjm / 2 + 1
102           IF (cos(rlatu(j)) / colat0 < 1. &
103    PRINT *, ' TRUNCATION AT ', imx              .and. rlamda(iim) * cos(rlatu(j)) < 1.) jfiltnu = j
104    
105    DO  j = 2, jjm/2 + 1         IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &
106       cof = cos(rlatu(j))/colat0              .and. rlamda(iim) * cos(rlatu(jjm - j + 2)) < 1.) &
107       IF (cof<1.) THEN              jfiltsu = jjm - j + 2
108          IF (rlamda(imx)*cos(rlatu(j))<1.) jfiltnu = j      END DO
109       END IF  
110        DO j = 1, jjm / 2
111       cof = cos(rlatu(jjp1-j+1))/colat0         IF (cos(rlatv(j)) / colat0 < 1. .and. rlamda(iim) * cos(rlatv(j)) < 1.) &
112       IF (cof<1.) THEN              jfiltnv = j
113          IF (rlamda(imx)*cos(rlatu(jjp1-j+1))<1.) jfiltsu = jjp1 - j + 1  
114       END IF         IF (cos(rlatv(jjm - j + 1)) / colat0 < 1. .and. rlamda(iim) &
115    END DO              * cos(rlatv(jjm - j + 1)) < 1.) jfiltsv = jjm - j + 1
116        END DO
117    DO  j = 1, jjm/2  
118       cof = cos(rlatv(j))/colat0      IF (jfiltnu <= 0) jfiltnu = 1
119       IF (cof<1.) THEN      IF (jfiltnu > jjm / 2 + 1) THEN
120          IF (rlamda(imx)*cos(rlatv(j))<1.) jfiltnv = j         PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu
121       END IF         STOP 1
122        END IF
123       cof = cos(rlatv(jjm-j+1))/colat0  
124       IF (cof<1.) THEN      IF (jfiltsu <= 0) jfiltsu = 1
125          IF (rlamda(imx)*cos(rlatv(jjm-j+1))<1.) jfiltsv = jjm - j + 1      IF (jfiltsu > jjm + 1) THEN
126       END IF         PRINT *, 'jfiltsu en dehors des valeurs acceptables ', jfiltsu
127    END DO         STOP 1
128        END IF
129    
130    IF (jfiltnu<=0) jfiltnu = 1      IF (jfiltnv <= 0) jfiltnv = 1
131    IF (jfiltnu>jjm/2+1) THEN      IF (jfiltnv > jjm / 2) THEN
132       PRINT *, ' jfiltnu en dehors des valeurs acceptables ', jfiltnu         PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv
133       STOP 1         STOP 1
134    END IF      END IF
135    
136    IF (jfiltsu<=0) jfiltsu = 1      IF (jfiltsv <= 0) jfiltsv = 1
137    IF (jfiltsu>jjm+1) THEN      IF (jfiltsv > jjm) THEN
138       PRINT *, ' jfiltsu en dehors des valeurs acceptables ', jfiltsu         PRINT *, 'jfiltsv en dehors des valeurs acceptables ', jfiltsv
139       STOP 1         STOP 1
140    END IF      END IF
141    
142    IF (jfiltnv<=0) jfiltnv = 1      PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
143    IF (jfiltnv>jjm/2) THEN           jfiltsu
144       PRINT *, ' jfiltnv en dehors des valeurs acceptables ', jfiltnv  
145       STOP 1      ! Determination de coefilu, coefilv, n=modfrstu, modfrstv
146    END IF  
147        DO j = 1, jjm
148    IF (jfiltsv<=0) jfiltsv = 1         modfrstu(j) = iim
149    IF (jfiltsv>jjm) THEN         modfrstv(j) = iim
150       PRINT *, ' jfiltsv en dehors des valeurs acceptables ', jfiltsv      END DO
151       STOP 1  
152    END IF      DO j = 2, jfiltnu
153           DO k = 2, iim
154    PRINT *, ' jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &            IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
155         jfiltsu         end DO
156           if (k == iim + 1) cycle
157           modfrstu(j) = k
158    !   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....        
159    !................................................................               kf = modfrstu(j)
160           DO k = kf, iim
161              cof = rlamda(k) * cos(rlatu(j))
162    DO  j = 1, jjm            coefilu(k, j) = cof - 1.
163       modfrstu(j) = iim            coefilu2(k, j) = cof**2 - 1.
164       modfrstv(j) = iim         end DO
165    END DO      END DO
166    
167    DO  j = 2, jfiltnu      DO j = 1, jfiltnv
168       DO  k = 2, modemax         DO k = 2, iim
169          cof = rlamda(k)*cos(rlatu(j))            IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
170          IF (cof<1.) GO TO 82         end DO
171       end DO         if (k == iim + 1) cycle
172       cycle         modfrstv(j) = k
173  82   modfrstu(j) = k  
174           kf = modfrstv(j)
175       kf = modfrstu(j)         DO k = kf, iim
176       DO  k = kf, modemax            cof = rlamda(k) * cos(rlatv(j))
177          cof = rlamda(k)*cos(rlatu(j))            coefilv(k, j) = cof - 1.
178          coefilu(k,j) = cof - 1.            coefilv2(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 = jfiltsu, jjm
183           DO k = 2, iim
184    DO j = 1, jfiltnv            IF (rlamda(k) * cos(rlatu(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         modfrstu(j) = k
188       end DO  
189       cycle         kf = modfrstu(j)
190  87   modfrstv(j) = k         DO k = kf, iim
191              cof = rlamda(k) * cos(rlatu(j))
192       kf = modfrstv(j)            coefilu(k, j) = cof - 1.
193       DO  k = kf, modemax            coefilu2(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      DO j = jfiltsv, jjm
198    end DO         DO k = 2, iim
199              IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
200    DO  j = jfiltsu, jjm         end DO
201       DO  k = 2, modemax         if (k == iim + 1) cycle
202          cof = rlamda(k)*cos(rlatu(j))         modfrstv(j) = k
203          IF (cof<1.) GO TO 92  
204       end DO         kf = modfrstv(j)
205       cycle         DO k = kf, iim
206  92   modfrstu(j) = k            cof = rlamda(k) * cos(rlatv(j))
207              coefilv(k, j) = cof - 1.
208       kf = modfrstu(j)            coefilv2(k, j) = cof**2 - 1.
209       DO  k = kf, modemax         end DO
210          cof = rlamda(k)*cos(rlatu(j))      END DO
211          coefilu(k,j) = cof - 1.  
212          coefilu2(k,j) = cof*cof - 1.      IF (jfiltnv>=jjm / 2 .OR. jfiltnu>=jjm / 2) THEN
213       end DO         IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
214    end DO         IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
215    
216    DO  j = jfiltsv, jjm         PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
217       DO  k = 2, modemax              jfiltsu
218          cof = rlamda(k)*cos(rlatv(j))      END IF
219          IF (cof<1.) GO TO 97  
220       end DO      PRINT *, 'Modes premiers v '
221       cycle      PRINT 334, modfrstv
222  97   modfrstv(j) = k      PRINT *, 'Modes premiers u '
223        PRINT 334, modfrstu
224       kf = modfrstv(j)  
225       DO  k = kf, modemax      allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
226          cof = rlamda(k)*cos(rlatv(j))      allocate(matricevn(iim, iim, jfiltnv))
227          coefilv(k,j) = cof - 1.      allocate(matricevs(iim, iim, jfiltsv:jjm))
228          coefilv2(k,j) = cof*cof - 1.      allocate(matriceus(iim, iim, jfiltsu:jjm), matrinvs(iim, iim, jfiltsu:jjm))
229       end DO  
230    END DO      ! Calcul de la matrice filtre 'matriceu' pour les champs situes
231        ! sur la grille scalaire
232    
233    IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN      DO j = 2, jfiltnu
234           DO i = 1, iim
235       IF (jfiltnv==jfiltsv) jfiltsv = 1 + jfiltnv            IF (i < modfrstu(j)) then
236       IF (jfiltnu==jfiltsu) jfiltsu = 1 + jfiltnu               coff = 0.
237              else
238       PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &               coff = coefilu(i, j)
239            jfiltsu            end IF
240    END IF            eignft(i, :) = eignfnv(:, i) * coff
241           END DO
242    PRINT *, '   Modes premiers  v  '         matriceun(:, :, j) = matmul(eignfnv, eignft)
243    PRINT 334, modfrstv      END DO
244    PRINT *, '   Modes premiers  u  '  
245    PRINT 334, modfrstu      DO j = jfiltsu, jjm
246           DO i = 1, iim
247              IF (i < modfrstu(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 = coefilu(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, :) = eignfnv(:, i) * coff
253            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &         END DO
254            jfiltnv, jjm - jfiltsv + 1         matriceus(:, :, j) = matmul(eignfnv, eignft)
255       STOP 1      END DO
256    END IF  
257    IF (nfilun>jfiltnu+2) THEN      ! Calcul de la matrice filtre 'matricev' pour les champs situes
258       PRINT *, ' le parametre nfilun utilise pour la matrice ', &      ! sur la grille de V ou de Z
259            ' matriceun est trop grand ! Gachis de memoire ! '  
260       PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnu      DO j = 1, jfiltnv
261       PRINT *, 'Pour information, nfilun,nfilus,nfilvn,nfilvs ', &         DO i = 1, iim
262            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &            IF (i < modfrstv(j)) then
263            jfiltnv, jjm - jfiltsv + 1               coff = 0.
264    END IF            else
265    IF (nfilus<jjm-jfiltsu+1) THEN               coff = coefilv(i, j)
266       PRINT *, ' le parametre nfilus utilise pour la matrice ', &            end IF
267            ' matriceus  est trop petit !  '            eignft(i, :) = eignfnu(:, i) * coff
268       PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &         END DO
269            jjm - jfiltsu + 1         matricevn(:, :, j) = matmul(eignfnu, eignft)
270       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &      END DO
271            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &  
272            jfiltnv, jjm - jfiltsv + 1      DO j = jfiltsv, jjm
273       STOP 1         DO i = 1, iim
274    END IF            IF (i < modfrstv(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 = coefilv(i, j)
278       PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &            end IF
279            jjm - jfiltsu + 1            eignft(i, :) = eignfnu(:, i) * coff
280       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &         END DO
281            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &         matricevs(:, :, j) = matmul(eignfnu, eignft)
282            jfiltnv, jjm - jfiltsv + 1      END DO
283    END IF  
284    IF (nfilvn<jfiltnv) THEN      ! Calcul de la matrice filtre 'matrinv' pour les champs situes
285       PRINT *, ' le parametre nfilvn utilise pour la matrice ', &      ! sur la grille scalaire , pour le filtre inverse
286            ' matricevn  est trop petit ! '  
287       PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnv      DO j = 2, jfiltnu
288       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &         DO i = 1, iim
289            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &            IF (i < modfrstu(j)) then
290            jfiltnv, jjm - jfiltsv + 1               coff = 0.
291       STOP 1            else
292    END IF               coff = coefilu(i, j) / (1. + coefilu(i, j))
293    IF (nfilvn>jfiltnv+2) THEN            end IF
294       PRINT *, ' le parametre nfilvn utilise pour la matrice ', &            eignft(i, :) = eignfnv(:, i) * coff
295            ' matricevn est trop grand !  Gachis de memoire ! '         END DO
296       PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnv         matrinvn(:, :, j) = matmul(eignfnv, eignft)
297       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &      END DO
298            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &  
299            jfiltnv, jjm - jfiltsv + 1      DO j = jfiltsu, jjm
300    END IF         DO i = 1, iim
301    IF (nfilvs<jjm-jfiltsv+1) THEN            IF (i < modfrstu(j)) then
302       PRINT *, ' le parametre nfilvs utilise pour la matrice ', &               coff = 0.
303            ' matricevs  est trop petit !  Le changer dans parafilt.h '            else
304       PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &               coff = coefilu(i, j) / (1. + coefilu(i, j))
305            jjm - jfiltsv + 1            end IF
306       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &            eignft(i, :) = eignfnv(:, i) * coff
307            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &         END DO
308            jfiltnv, jjm - jfiltsv + 1         matrinvs(:, :, j) = matmul(eignfnv, eignft)
309       STOP 1      END DO
   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  
310    
311    END DO  334 FORMAT (1X, 24I3)
312    
313  334 FORMAT (1X,24I3)    END SUBROUTINE inifilr
 755 FORMAT (1X,6F10.3,I3)  
314    
315  END SUBROUTINE inifilr  end module inifilr_m
Legend:
Removed from v.32  
changed lines
  Added in v.140
  ViewVC Help
Powered by ViewVC 1.1.21