Sources/filtrez/inifilr.f

SUBROUTINE inifilr

  ! 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

  IMPLICIT NONE

  REAL dlonu(iim), dlatu(jjm)
  REAL rlamda(iim), eignvl(iim)

  REAL lamdamax, pi, cof
  INTEGER i, j, modemax, imx, k, kf, ii
  REAL dymin, dxmin, colat0
  REAL eignft(iim,iim), coff
  EXTERNAL inifgn

  !-----------------------------------------------------------            


  pi = 2.*asin(1.)

  DO i = 1, iim
     dlonu(i) = xprimu(i)
  END DO

  CALL inifgn(eignvl)

  PRINT *, ' EIGNVL '
  PRINT 250, eignvl
250 FORMAT (1X,5E13.6)

  ! compute eigenvalues and eigenfunctions                                


  !.................................................................      

  !  compute the filtering coefficients for scalar lines and              
  !  meridional wind v-lines                                              

  !  we filter all those latitude lines where coefil < 1                  
  !  NO FILTERING AT POLES                                                

  !  colat0 is to be used  when alpha (stretching coefficient)            
  !  is set equal to zero for the regular grid case                       

  !    .......   Calcul  de  colat0   .........                           
  !     .....  colat0 = minimum de ( 0.5, min dy/ min dx )   ...          


  DO  j = 1, jjm
     dlatu(j) = rlatu(j) - rlatu(j+1)
  END DO

  dxmin = dlonu(1)
  DO i = 2, iim
     dxmin = min(dxmin,dlonu(i))
  END DO
  dymin = dlatu(1)
  DO j = 2, jjm
     dymin = min(dymin,dlatu(j))
  END DO


  colat0 = min(0.5,dymin/dxmin)

  IF ( .NOT. fxyhypb .AND. ysinus) THEN
     colat0 = 0.6
     !         ...... a revoir  pour  ysinus !   .......                     
     alphax = 0.
  END IF

  PRINT 50, colat0, alphax
50 FORMAT (/15X,' Inifilr colat0 alphax ',2E16.7)

  IF (alphax==1.) THEN
     PRINT *, ' Inifilr  alphax doit etre  <  a 1.  Corriger '
     STOP 1
  END IF

  lamdamax = iim/(pi*colat0*(1.-alphax))

  DO  i = 2, iim
     rlamda(i) = lamdamax/sqrt(abs(eignvl(i)))
  END DO


  DO  j = 1, jjm
     DO  i = 1, iim
        coefilu(i,j) = 0.0
        coefilv(i,j) = 0.0
        coefilu2(i,j) = 0.0
        coefilv2(i,j) = 0.0
     end DO
  END DO


  !    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....          
  !    .........................................................          

  modemax = iim

  !ccc    imx = modemax - 4 * (modemax/iim)                               

  imx = iim

  PRINT *, ' TRUNCATION AT ', imx

  DO  j = 2, jjm/2 + 1
     cof = cos(rlatu(j))/colat0
     IF (cof<1.) THEN
        IF (rlamda(imx)*cos(rlatu(j))<1.) jfiltnu = j
     END IF

     cof = cos(rlatu(jjp1-j+1))/colat0
     IF (cof<1.) THEN
        IF (rlamda(imx)*cos(rlatu(jjp1-j+1))<1.) jfiltsu = jjp1 - j + 1
     END IF
  END DO

  DO  j = 1, jjm/2
     cof = cos(rlatv(j))/colat0
     IF (cof<1.) THEN
        IF (rlamda(imx)*cos(rlatv(j))<1.) jfiltnv = j
     END IF

     cof = cos(rlatv(jjm-j+1))/colat0
     IF (cof<1.) THEN
        IF (rlamda(imx)*cos(rlatv(jjm-j+1))<1.) jfiltsv = jjm - j + 1
     END IF
  END DO


  IF (jfiltnu<=0) jfiltnu = 1
  IF (jfiltnu>jjm/2+1) THEN
     PRINT *, ' jfiltnu en dehors des valeurs acceptables ', jfiltnu
     STOP 1
  END IF

  IF (jfiltsu<=0) jfiltsu = 1
  IF (jfiltsu>jjm+1) THEN
     PRINT *, ' jfiltsu en dehors des valeurs acceptables ', jfiltsu
     STOP 1
  END IF

  IF (jfiltnv<=0) jfiltnv = 1
  IF (jfiltnv>jjm/2) THEN
     PRINT *, ' jfiltnv en dehors des valeurs acceptables ', jfiltnv
     STOP 1
  END IF

  IF (jfiltsv<=0) jfiltsv = 1
  IF (jfiltsv>jjm) THEN
     PRINT *, ' jfiltsv en dehors des valeurs acceptables ', jfiltsv
     STOP 1
  END IF

  PRINT *, ' jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
       jfiltsu


  !   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....       
  !................................................................       


  DO  j = 1, jjm
     modfrstu(j) = iim
     modfrstv(j) = iim
  END DO

  DO  j = 2, jfiltnu
     DO  k = 2, modemax
        cof = rlamda(k)*cos(rlatu(j))
        IF (cof<1.) GO TO 82
     end DO
     cycle
82   modfrstu(j) = k

     kf = modfrstu(j)
     DO  k = kf, modemax
        cof = rlamda(k)*cos(rlatu(j))
        coefilu(k,j) = cof - 1.
        coefilu2(k,j) = cof*cof - 1.
     end DO
  END DO


  DO j = 1, jfiltnv
     DO  k = 2, modemax
        cof = rlamda(k)*cos(rlatv(j))
        IF (cof<1.) GO TO 87
     end DO
     cycle
87   modfrstv(j) = k

     kf = modfrstv(j)
     DO  k = kf, modemax
        cof = rlamda(k)*cos(rlatv(j))
        coefilv(k,j) = cof - 1.
        coefilv2(k,j) = cof*cof - 1.
     end DO
  end DO

  DO  j = jfiltsu, jjm
     DO  k = 2, modemax
        cof = rlamda(k)*cos(rlatu(j))
        IF (cof<1.) GO TO 92
     end DO
     cycle
92   modfrstu(j) = k

     kf = modfrstu(j)
     DO  k = kf, modemax
        cof = rlamda(k)*cos(rlatu(j))
        coefilu(k,j) = cof - 1.
        coefilu2(k,j) = cof*cof - 1.
     end DO
  end DO

  DO  j = jfiltsv, jjm
     DO  k = 2, modemax
        cof = rlamda(k)*cos(rlatv(j))
        IF (cof<1.) GO TO 97
     end DO
     cycle
97   modfrstv(j) = k

     kf = modfrstv(j)
     DO  k = kf, modemax
        cof = rlamda(k)*cos(rlatv(j))
        coefilv(k,j) = cof - 1.
        coefilv2(k,j) = cof*cof - 1.
     end DO
  END DO


  IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN

     IF (jfiltnv==jfiltsv) jfiltsv = 1 + jfiltnv
     IF (jfiltnu==jfiltsu) jfiltsu = 1 + jfiltnu

     PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
          jfiltsu
  END IF

  PRINT *, '   Modes premiers  v  '
  PRINT 334, modfrstv
  PRINT *, '   Modes premiers  u  '
  PRINT 334, modfrstu


  IF (nfilun<jfiltnu) THEN
     PRINT *, ' le parametre nfilun utilise pour la matrice ', &
          ' matriceun  est trop petit ! '
     PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnu
     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 (nfilun>jfiltnu+2) THEN
     PRINT *, ' le parametre nfilun utilise pour la matrice ', &
          ' matriceun est trop grand ! Gachis de memoire ! '
     PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnu
     PRINT *, 'Pour information, nfilun,nfilus,nfilvn,nfilvs ', &
          'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
          jfiltnv, jjm - jfiltsv + 1
  END IF
  IF (nfilus<jjm-jfiltsu+1) THEN
     PRINT *, ' le parametre nfilus utilise pour la matrice ', &
          ' matriceus  est trop petit !  '
     PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &
          jjm - jfiltsu + 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 (nfilus>jjm-jfiltsu+3) THEN
     PRINT *, ' le parametre nfilus utilise pour la matrice ', &
          ' matriceus  est trop grand ! '
     PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &
          jjm - jfiltsu + 1
     PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &
          'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
          jfiltnv, jjm - jfiltsv + 1
  END IF
  IF (nfilvn<jfiltnv) THEN
     PRINT *, ' le parametre nfilvn utilise pour la matrice ', &
          ' matricevn  est trop petit ! '
     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
     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

334 FORMAT (1X,24I3)
755 FORMAT (1X,6F10.3,I3)

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