/[lmdze]/trunk/filtrez/inifilr.f90
ViewVC logotype

Diff of /trunk/filtrez/inifilr.f90

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