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