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