/[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 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
314    
315  END SUBROUTINE inifilr  end module inifilr_m
Legend:
Removed from v.40  
changed lines
  Added in v.140
  ViewVC Help
Powered by ViewVC 1.1.21