/[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

revision 40 by guez, Tue Feb 22 13:49:36 2011 UTC revision 54 by guez, Tue Dec 6 15:07:04 2011 UTC
# Line 1  Line 1 
1  SUBROUTINE inifilr  module inifilr_m
2    
3    ! From filtrez/inifilr.F,v 1.1.1.1 2004/05/19 12:53:09    use dimens_m, only: iim
   ! 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  
4    
5    IMPLICIT NONE    IMPLICIT NONE
6    
7    REAL dlonu(iim), dlatu(jjm)    INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv
8    REAL rlamda(iim), eignvl(iim)    INTEGER, PARAMETER:: nfilun=3, nfilus=2, nfilvn=2, nfilvs=2
9    
10    REAL lamdamax, pi, cof    real matriceun(iim,iim,nfilun), matriceus(iim,iim,nfilus)
11    INTEGER i, j, modemax, imx, k, kf, ii    real matricevn(iim,iim,nfilvn), matricevs(iim,iim,nfilvs)
12    REAL dymin, dxmin, colat0    real matrinvn(iim,iim,nfilun), matrinvs(iim,iim,nfilus)
13    REAL eignft(iim,iim), coff  
14    EXTERNAL inifgn    private iim, nfilun, nfilus, nfilvn, nfilvs
15    
16    !-----------------------------------------------------------              contains
17    
18      SUBROUTINE inifilr
19    pi = 2.*asin(1.)  
20        ! From filtrez/inifilr.F, version 1.1.1.1 2004/05/19 12:53:09
21    DO i = 1, iim      ! H. Upadhyaya, O. Sharma
22       dlonu(i) = xprimu(i)  
23    END DO      ! This routine computes the eigenfunctions of the laplacian on the
24        ! stretched grid, and the filtering coefficients.
25    CALL inifgn(eignvl)      ! We designate:
26        ! eignfn eigenfunctions of the discrete laplacian
27    PRINT *, ' EIGNVL '      ! eigenvl eigenvalues
28    PRINT 250, eignvl      ! jfiltn index of the last scalar line filtered in NH
29  250 FORMAT (1X,5E13.6)      ! jfilts index of the first line filtered in SH
30        ! modfrst index of the mode from where modes are filtered
31    ! compute eigenvalues and eigenfunctions                                      ! modemax maximum number of modes (im)
32        ! coefil filtering coefficients (lamda_max * cos(rlat) / lamda)
33        ! sdd SQRT(dx)
34    !.................................................................        
35        ! The modes are filtered from modfrst to modemax.
36    !  compute the filtering coefficients for scalar lines and                
37    !  meridional wind v-lines                                                    USE dimens_m, ONLY : iim, jjm
38        USE logic, ONLY : fxyhypb, ysinus
39    !  we filter all those latitude lines where coefil < 1                        USE comgeom, ONLY : rlatu, rlatv, xprimu
40    !  NO FILTERING AT POLES                                                      use nr_util, only: pi
41        USE serre, ONLY : alphax
42    !  colat0 is to be used  when alpha (stretching coefficient)                  USE coefils, ONLY : coefilu, coefilu2, coefilv, coefilv2, eignfnu, &
43    !  is set equal to zero for the regular grid case                                 eignfnv, modfrstu, modfrstv
44    
45    !    .......   Calcul  de  colat0   .........                                ! Local:
46    !     .....  colat0 = minimum de ( 0.5, min dy/ min dx )   ...                REAL dlonu(iim), dlatu(jjm)
47        REAL rlamda(2: iim), eignvl(iim)
48    
49    DO  j = 1, jjm      REAL lamdamax, cof
50       dlatu(j) = rlatu(j) - rlatu(j+1)      INTEGER i, j, modemax, imx, k, kf
51    END DO      REAL dymin, dxmin, colat0
52        REAL eignft(iim, iim), coff
53    dxmin = dlonu(1)      EXTERNAL inifgn
54    DO i = 2, iim  
55       dxmin = min(dxmin,dlonu(i))      !-----------------------------------------------------------
56    END DO  
57    dymin = dlatu(1)      print *, "Call sequence information: inifilr"
58    DO j = 2, jjm  
59       dymin = min(dymin,dlatu(j))      DO i = 1, iim
60    END DO         dlonu(i) = xprimu(i)
61        END DO
62    
63    colat0 = min(0.5,dymin/dxmin)      CALL inifgn(eignvl)
64    
65    IF ( .NOT. fxyhypb .AND. ysinus) THEN      PRINT *, 'EIGNVL '
66       colat0 = 0.6      PRINT "(1X, 5E13.6)", eignvl
67       !         ...... a revoir  pour  ysinus !   .......                      
68       alphax = 0.      ! compute eigenvalues and eigenfunctions
69    END IF      ! compute the filtering coefficients for scalar lines and
70        ! meridional wind v-lines
71    PRINT 50, colat0, alphax      ! we filter all those latitude lines where coefil < 1
72  50 FORMAT (/15X,' Inifilr colat0 alphax ',2E16.7)      ! NO FILTERING AT POLES
73        ! colat0 is to be used when alpha (stretching coefficient)
74    IF (alphax==1.) THEN      ! is set equal to zero for the regular grid case
75       PRINT *, ' Inifilr  alphax doit etre  <  a 1.  Corriger '  
76       STOP 1      ! Calcul de colat0
77    END IF  
78        DO j = 1, jjm
79    lamdamax = iim/(pi*colat0*(1.-alphax))         dlatu(j) = rlatu(j) - rlatu(j+1)
80        END DO
81    DO  i = 2, iim  
82       rlamda(i) = lamdamax/sqrt(abs(eignvl(i)))      dxmin = dlonu(1)
83    END DO      DO i = 2, iim
84           dxmin = min(dxmin, dlonu(i))
85        END DO
86    DO  j = 1, jjm      dymin = dlatu(1)
87       DO  i = 1, iim      DO j = 2, jjm
88          coefilu(i,j) = 0.0         dymin = min(dymin, dlatu(j))
89          coefilv(i,j) = 0.0      END DO
90          coefilu2(i,j) = 0.0  
91          coefilv2(i,j) = 0.0      colat0 = min(0.5, dymin/dxmin)
92       end DO  
93    END DO      IF (.NOT. fxyhypb .AND. ysinus) THEN
94           colat0 = 0.6
95           ! À revoir pour ysinus
96    !    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....                   alphax = 0.
97    !    .........................................................                END IF
98    
99    modemax = iim      PRINT *, 'colat0 = ', colat0
100        PRINT *, 'alphax = ', alphax
101    !ccc    imx = modemax - 4 * (modemax/iim)                                
102        IF (alphax == 1.) THEN
103    imx = iim         PRINT *, 'alphax doit etre < a 1. Corriger '
104           STOP 1
105    PRINT *, ' TRUNCATION AT ', imx      END IF
106    
107    DO  j = 2, jjm/2 + 1      lamdamax = iim / (pi * colat0 * (1. - alphax))
108       cof = cos(rlatu(j))/colat0      rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))
109       IF (cof<1.) THEN  
110          IF (rlamda(imx)*cos(rlatu(j))<1.) jfiltnu = j      DO j = 1, jjm
111       END IF         DO i = 1, iim
112              coefilu(i, j) = 0.
113       cof = cos(rlatu(jjp1-j+1))/colat0            coefilv(i, j) = 0.
114       IF (cof<1.) THEN            coefilu2(i, j) = 0.
115          IF (rlamda(imx)*cos(rlatu(jjp1-j+1))<1.) jfiltsu = jjp1 - j + 1            coefilv2(i, j) = 0.
116       END IF         end DO
117    END DO      END DO
118    
119    DO  j = 1, jjm/2      ! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
120       cof = cos(rlatv(j))/colat0  
121       IF (cof<1.) THEN      modemax = iim
122          IF (rlamda(imx)*cos(rlatv(j))<1.) jfiltnv = j      imx = iim
123       END IF  
124        PRINT *, 'TRUNCATION AT ', imx
125       cof = cos(rlatv(jjm-j+1))/colat0  
126       IF (cof<1.) THEN      DO j = 2, jjm / 2 + 1
127          IF (rlamda(imx)*cos(rlatv(jjm-j+1))<1.) jfiltsv = jjm - j + 1         IF (cos(rlatu(j)) / colat0 < 1. &
128       END IF              .and. rlamda(imx) * cos(rlatu(j)) < 1.) jfiltnu = j
129    END DO  
130           IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &
131                .and. rlamda(imx) * cos(rlatu(jjm - j + 2)) < 1.) &
132    IF (jfiltnu<=0) jfiltnu = 1              jfiltsu = jjm - j + 2
133    IF (jfiltnu>jjm/2+1) THEN      END DO
134       PRINT *, ' jfiltnu en dehors des valeurs acceptables ', jfiltnu  
135       STOP 1      DO j = 1, jjm/2
136    END IF         cof = cos(rlatv(j))/colat0
137           IF (cof < 1.) THEN
138    IF (jfiltsu<=0) jfiltsu = 1            IF (rlamda(imx)*cos(rlatv(j)) < 1.) jfiltnv = j
139    IF (jfiltsu>jjm+1) THEN         END IF
140       PRINT *, ' jfiltsu en dehors des valeurs acceptables ', jfiltsu  
141       STOP 1         cof = cos(rlatv(jjm-j+1))/colat0
142    END IF         IF (cof < 1.) THEN
143              IF (rlamda(imx)*cos(rlatv(jjm-j+1)) < 1.) jfiltsv = jjm - j + 1
144    IF (jfiltnv<=0) jfiltnv = 1         END IF
145    IF (jfiltnv>jjm/2) THEN      END DO
146       PRINT *, ' jfiltnv en dehors des valeurs acceptables ', jfiltnv  
147       STOP 1      IF (jfiltnu <= 0) jfiltnu = 1
148    END IF      IF (jfiltnu > jjm/2+1) THEN
149           PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu
150    IF (jfiltsv<=0) jfiltsv = 1         STOP 1
151    IF (jfiltsv>jjm) THEN      END IF
152       PRINT *, ' jfiltsv en dehors des valeurs acceptables ', jfiltsv  
153       STOP 1      IF (jfiltsu <= 0) jfiltsu = 1
154    END IF      IF (jfiltsu > jjm + 1) THEN
155           PRINT *, 'jfiltsu en dehors des valeurs acceptables ', jfiltsu
156    PRINT *, ' jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &         STOP 1
157         jfiltsu      END IF
158    
159        IF (jfiltnv <= 0) jfiltnv = 1
160    !   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....            IF (jfiltnv > jjm/2) THEN
161    !................................................................               PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv
162           STOP 1
163        END IF
164    DO  j = 1, jjm  
165       modfrstu(j) = iim      IF (jfiltsv <= 0) jfiltsv = 1
166       modfrstv(j) = iim      IF (jfiltsv > jjm) THEN
167    END DO         PRINT *, 'jfiltsv en dehors des valeurs acceptables ', jfiltsv
168           STOP 1
169    DO  j = 2, jfiltnu      END IF
170       DO  k = 2, modemax  
171          cof = rlamda(k)*cos(rlatu(j))      PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
172          IF (cof<1.) GO TO 82           jfiltsu
173       end DO  
174       cycle      ! Determination de coefilu, coefilv, n=modfrstu, modfrstv
175  82   modfrstu(j) = k  
176        DO j = 1, jjm
177       kf = modfrstu(j)         modfrstu(j) = iim
178       DO  k = kf, modemax         modfrstv(j) = iim
179          cof = rlamda(k)*cos(rlatu(j))      END DO
180          coefilu(k,j) = cof - 1.  
181          coefilu2(k,j) = cof*cof - 1.      DO j = 2, jfiltnu
182       end DO         DO k = 2, modemax
183    END DO            cof = rlamda(k) * cos(rlatu(j))
184              IF (cof < 1.) exit
185           end DO
186    DO j = 1, jfiltnv         if (k == modemax + 1) cycle
187       DO  k = 2, modemax         modfrstu(j) = k
188          cof = rlamda(k)*cos(rlatv(j))  
189          IF (cof<1.) GO TO 87         kf = modfrstu(j)
190       end DO         DO k = kf, modemax
191       cycle            cof = rlamda(k)*cos(rlatu(j))
192  87   modfrstv(j) = k            coefilu(k, j) = cof - 1.
193              coefilu2(k, j) = cof*cof - 1.
194       kf = modfrstv(j)         end DO
195       DO  k = kf, modemax      END DO
196          cof = rlamda(k)*cos(rlatv(j))  
197          coefilv(k,j) = cof - 1.      DO j = 1, jfiltnv
198          coefilv2(k,j) = cof*cof - 1.         DO k = 2, modemax
199       end DO            cof = rlamda(k)*cos(rlatv(j))
200    end DO            IF (cof < 1.) exit
201           end DO
202    DO  j = jfiltsu, jjm         if (k == modemax + 1) cycle
203       DO  k = 2, modemax         modfrstv(j) = k
204          cof = rlamda(k)*cos(rlatu(j))  
205          IF (cof<1.) GO TO 92         kf = modfrstv(j)
206       end DO         DO k = kf, modemax
207       cycle            cof = rlamda(k)*cos(rlatv(j))
208  92   modfrstu(j) = k            coefilv(k, j) = cof - 1.
209              coefilv2(k, j) = cof*cof - 1.
210       kf = modfrstu(j)         end DO
211       DO  k = kf, modemax      end DO
212          cof = rlamda(k)*cos(rlatu(j))  
213          coefilu(k,j) = cof - 1.      DO j = jfiltsu, jjm
214          coefilu2(k,j) = cof*cof - 1.         DO k = 2, modemax
215       end DO            cof = rlamda(k)*cos(rlatu(j))
216    end DO            IF (cof < 1.) exit
217           end DO
218    DO  j = jfiltsv, jjm         if (k == modemax + 1) cycle
219       DO  k = 2, modemax         modfrstu(j) = k
220          cof = rlamda(k)*cos(rlatv(j))  
221          IF (cof<1.) GO TO 97         kf = modfrstu(j)
222       end DO         DO k = kf, modemax
223       cycle            cof = rlamda(k)*cos(rlatu(j))
224  97   modfrstv(j) = k            coefilu(k, j) = cof - 1.
225              coefilu2(k, j) = cof*cof - 1.
226       kf = modfrstv(j)         end DO
227       DO  k = kf, modemax      end DO
228          cof = rlamda(k)*cos(rlatv(j))  
229          coefilv(k,j) = cof - 1.      DO j = jfiltsv, jjm
230          coefilv2(k,j) = cof*cof - 1.         DO k = 2, modemax
231       end DO            cof = rlamda(k)*cos(rlatv(j))
232    END DO            IF (cof < 1.) exit
233           end DO
234           if (k == modemax + 1) cycle
235    IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN         modfrstv(j) = k
236    
237       IF (jfiltnv==jfiltsv) jfiltsv = 1 + jfiltnv         kf = modfrstv(j)
238       IF (jfiltnu==jfiltsu) jfiltsu = 1 + jfiltnu         DO k = kf, modemax
239              cof = rlamda(k)*cos(rlatv(j))
240       PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &            coefilv(k, j) = cof - 1.
241            jfiltsu            coefilv2(k, j) = cof*cof - 1.
242    END IF         end DO
243        END DO
244    PRINT *, '   Modes premiers  v  '  
245    PRINT 334, modfrstv      IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN
246    PRINT *, '   Modes premiers  u  '         IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
247    PRINT 334, modfrstu         IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
248    
249           PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
250    IF (nfilun<jfiltnu) THEN              jfiltsu
251       PRINT *, ' le parametre nfilun utilise pour la matrice ', &      END IF
252            ' matriceun  est trop petit ! '  
253       PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnu      PRINT *, 'Modes premiers v '
254       PRINT *, 'Pour information, nfilun,nfilus,nfilvn,nfilvs ', &      PRINT 334, modfrstv
255            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &      PRINT *, 'Modes premiers u '
256            jfiltnv, jjm - jfiltsv + 1      PRINT 334, modfrstu
257       STOP 1  
258    END IF      IF (nfilun < jfiltnu) THEN
259    IF (nfilun>jfiltnu+2) THEN         PRINT *, 'le parametre nfilun utilise pour la matrice ', &
260       PRINT *, ' le parametre nfilun utilise pour la matrice ', &              'matriceun est trop petit ! '
261            ' matriceun est trop grand ! Gachis de memoire ! '         PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnu
262       PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnu         PRINT *, 'Pour information, nfilun, nfilus, nfilvn, nfilvs ', &
263       PRINT *, 'Pour information, nfilun,nfilus,nfilvn,nfilvs ', &              'doivent etre egaux successivement a ', jfiltnu, &
264            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &              jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1
265            jfiltnv, jjm - jfiltsv + 1         STOP 1
266    END IF      END IF
267    IF (nfilus<jjm-jfiltsu+1) THEN      IF (nfilun > jfiltnu+2) THEN
268       PRINT *, ' le parametre nfilus utilise pour la matrice ', &         PRINT *, 'le parametre nfilun utilise pour la matrice ', &
269            ' matriceus  est trop petit !  '              'matriceun est trop grand ! Gachis de memoire ! '
270       PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &         PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnu
271            jjm - jfiltsu + 1         PRINT *, 'Pour information, nfilun, nfilus, nfilvn, nfilvs ', &
272       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &              'doivent etre egaux successivement a ', &
273            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &              jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1
274            jfiltnv, jjm - jfiltsv + 1      END IF
275       STOP 1      IF (nfilus < jjm-jfiltsu+1) THEN
276    END IF         PRINT *, 'le parametre nfilus utilise pour la matrice ', &
277    IF (nfilus>jjm-jfiltsu+3) THEN              'matriceus est trop petit ! '
278       PRINT *, ' le parametre nfilus utilise pour la matrice ', &         PRINT *, 'Le changer dans parafilt.h et le mettre a ', &
279            ' matriceus  est trop grand ! '              jjm - jfiltsu + 1
280       PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &         PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &
281            jjm - jfiltsu + 1              'doivent etre egaux successivement a ', &
282       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &              jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1
283            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &         STOP 1
284            jfiltnv, jjm - jfiltsv + 1      END IF
285    END IF      IF (nfilus > jjm-jfiltsu+3) THEN
286    IF (nfilvn<jfiltnv) THEN         PRINT *, 'le parametre nfilus utilise pour la matrice ', &
287       PRINT *, ' le parametre nfilvn utilise pour la matrice ', &              'matriceus est trop grand ! '
288            ' matricevn  est trop petit ! '         PRINT *, 'Le changer dans parafilt.h et le mettre a ', &
289       PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnv              jjm - jfiltsu + 1
290       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &         PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &
291            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &              'doivent etre egaux successivement a ', &
292            jfiltnv, jjm - jfiltsv + 1              jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1
293       STOP 1      END IF
294    END IF      IF (nfilvn < jfiltnv) THEN
295    IF (nfilvn>jfiltnv+2) THEN         PRINT *, 'le parametre nfilvn utilise pour la matrice ', &
296       PRINT *, ' le parametre nfilvn utilise pour la matrice ', &              'matricevn est trop petit ! '
297            ' matricevn est trop grand !  Gachis de memoire ! '         PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnv
298       PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnv         PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &
299       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &              'doivent etre egaux successivement a ', &
300            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &              jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1
301            jfiltnv, jjm - jfiltsv + 1         STOP 1
302    END IF      END IF
303    IF (nfilvs<jjm-jfiltsv+1) THEN      IF (nfilvn > jfiltnv+2) THEN
304       PRINT *, ' le parametre nfilvs utilise pour la matrice ', &         PRINT *, 'le parametre nfilvn utilise pour la matrice ', &
305            ' matricevs  est trop petit !  Le changer dans parafilt.h '              'matricevn est trop grand ! Gachis de memoire ! '
306       PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &         PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnv
307            jjm - jfiltsv + 1         PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &
308       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &              'doivent etre egaux successivement a ', &
309            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &              jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1
310            jfiltnv, jjm - jfiltsv + 1      END IF
311       STOP 1      IF (nfilvs < jjm-jfiltsv+1) THEN
312    END IF         PRINT *, 'le parametre nfilvs utilise pour la matrice ', &
313    IF (nfilvs>jjm-jfiltsv+3) THEN              'matricevs est trop petit ! Le changer dans parafilt.h '
314       PRINT *, ' le parametre nfilvs utilise pour la matrice ', &         PRINT *, 'Le changer dans parafilt.h et le mettre a ', &
315            ' matricevs  est trop grand ! Gachis de memoire ! '              jjm - jfiltsv + 1
316       PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &         PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &
317            jjm - jfiltsv + 1              'doivent etre egaux successivement a ', &
318       PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &              jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1
319            'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &         STOP 1
320            jfiltnv, jjm - jfiltsv + 1      END IF
321    END IF      IF (nfilvs > jjm-jfiltsv+3) THEN
322           PRINT *, 'le parametre nfilvs utilise pour la matrice ', &
323    !   ... Calcul de la matrice filtre 'matriceu'  pour les champs situes                'matricevs est trop grand ! Gachis de memoire ! '
324    !                       sur la grille scalaire                 ........         PRINT *, 'Le changer dans parafilt.h et le mettre a ', &
325                jjm - jfiltsv + 1
326    DO j = 2, jfiltnu         PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &
327                'doivent etre egaux successivement a ', &
328       DO i = 1, iim              jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1
329          coff = coefilu(i,j)      END IF
330          IF (i<modfrstu(j)) coff = 0.  
331          DO k = 1, iim      ! Calcul de la matrice filtre 'matriceu' pour les champs situes
332             eignft(i,k) = eignfnv(k,i)*coff      ! sur la grille scalaire
333          END DO  
334       END DO      DO j = 2, jfiltnu
335       DO k = 1, iim         DO i = 1, iim
336          DO i = 1, iim            IF (i < modfrstu(j)) then
337             matriceun(i,k,j) = 0.0               coff = 0.
338             DO ii = 1, iim            else
339                matriceun(i,k,j) = matriceun(i,k,j) + eignfnv(i,ii)*eignft(ii,k)               coff = coefilu(i, j)
340             END DO            end IF
341          END DO            eignft(i, :) = eignfnv(:, i)*coff
342       END DO         END DO
343           matriceun(:, :, j) = matmul(eignfnv, eignft)
344    END DO      END DO
345    
346    DO j = jfiltsu, jjm      DO j = jfiltsu, jjm
347           DO i = 1, iim
348       DO i = 1, iim            IF (i < modfrstu(j)) then
349          coff = coefilu(i,j)               coff = 0.
350          IF (i<modfrstu(j)) coff = 0.            else
351          DO k = 1, iim               coff = coefilu(i, j)
352             eignft(i,k) = eignfnv(k,i)*coff            end IF
353          END DO            eignft(i, :) = eignfnv(:, i) * coff
354       END DO         END DO
355       DO k = 1, iim         matriceus(:, :, j - jfiltsu + 1) = matmul(eignfnv, eignft)
356          DO i = 1, iim      END DO
357             matriceus(i,k,j-jfiltsu+1) = 0.0  
358             DO ii = 1, iim      ! Calcul de la matrice filtre 'matricev' pour les champs situes
359                matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) + &      ! sur la grille de V ou de Z
360                     eignfnv(i,ii)*eignft(ii,k)  
361             END DO      DO j = 1, jfiltnv
362          END DO         DO i = 1, iim
363       END DO            IF (i < modfrstv(j)) then
364                 coff = 0.
365    END DO            else
366                 coff = coefilv(i, j)
367    !   ...................................................................            end IF
368              eignft(i, :) = eignfnu(:, i)*coff
369    !   ... Calcul de la matrice filtre 'matricev'  pour les champs situes           END DO
370    !                       sur la grille   de V ou de Z           ........         matricevn(:, :, j) = matmul(eignfnu, eignft)
371    !   ...................................................................      END DO
372    
373    DO j = 1, jfiltnv      DO j = jfiltsv, jjm
374           DO i = 1, iim
375       DO i = 1, iim            IF (i < modfrstv(j)) then
376          coff = coefilv(i,j)               coff = 0.
377          IF (i<modfrstv(j)) coff = 0.            else
378          DO k = 1, iim               coff = coefilv(i, j)
379             eignft(i,k) = eignfnu(k,i)*coff            end IF
380          END DO            eignft(i, :) = eignfnu(:, i)*coff
381       END DO         END DO
382       DO k = 1, iim         matricevs(:, :, j-jfiltsv+1) = matmul(eignfnu, eignft)
383          DO i = 1, iim      END DO
384             matricevn(i,k,j) = 0.0  
385             DO ii = 1, iim      ! Calcul de la matrice filtre 'matrinv' pour les champs situes
386                matricevn(i,k,j) = matricevn(i,k,j) + eignfnu(i,ii)*eignft(ii,k)      ! sur la grille scalaire , pour le filtre inverse
387             END DO  
388          END DO      DO j = 2, jfiltnu
389       END DO         DO i = 1, iim
390              IF (i < modfrstu(j)) then
391    END DO               coff = 0.
392              else
393    DO j = jfiltsv, jjm               coff = coefilu(i, j)/(1.+coefilu(i, j))
394              end IF
395       DO i = 1, iim            eignft(i, :) = eignfnv(:, i)*coff
396          coff = coefilv(i,j)         END DO
397          IF (i<modfrstv(j)) coff = 0.         matrinvn(:, :, j) = matmul(eignfnv, eignft)
398          DO k = 1, iim      END DO
399             eignft(i,k) = eignfnu(k,i)*coff  
400          END DO      DO j = jfiltsu, jjm
401       END DO         DO i = 1, iim
402       DO k = 1, iim            IF (i < modfrstu(j)) then
403          DO i = 1, iim               coff = 0.
404             matricevs(i,k,j-jfiltsv+1) = 0.0            else
405             DO ii = 1, iim               coff = coefilu(i, j)/(1.+coefilu(i, j))
406                matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) + &            end IF
407                     eignfnu(i,ii)*eignft(ii,k)            eignft(i, :) = eignfnv(:, i)*coff
408             END DO         END DO
409          END DO         matrinvs(:, :, j-jfiltsu+1) = matmul(eignfnv, eignft)
410       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  
411    
412    END DO  334 FORMAT (1X, 24I3)
413    
414  334 FORMAT (1X,24I3)    END SUBROUTINE inifilr
415    
416  END SUBROUTINE inifilr  end module inifilr_m
Legend:
Removed from v.40  
changed lines
  Added in v.54
  ViewVC Help
Powered by ViewVC 1.1.21