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