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