trunk/filtrez/inifilr.f

module inifilr_m

  use dimens_m, only: iim

  IMPLICIT NONE

  INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv
  INTEGER, PARAMETER:: nfilun=3, nfilus=2, nfilvn=2, nfilvs=2

  real matriceun(iim,iim,nfilun), matriceus(iim,iim,nfilus)
  real matricevn(iim,iim,nfilvn), matricevs(iim,iim,nfilvs)
  real matrinvn(iim,iim,nfilun), matrinvs(iim,iim,nfilus)

  private iim, nfilun, nfilus, nfilvn, nfilvs

contains

  SUBROUTINE inifilr

    ! From filtrez/inifilr.F, version 1.1.1.1 2004/05/19 12:53:09
    ! H. Upadhyaya, O. Sharma

    ! This routine computes the eigenfunctions of the laplacian on the
    ! stretched grid, and the filtering coefficients.
    ! We designate:
    ! eignfn eigenfunctions of the discrete laplacian
    ! eigenvl eigenvalues
    ! jfiltn index of 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 coefils, ONLY : coefilu, coefilu2, coefilv, coefilv2, eignfnu, &
         eignfnv, modfrstu, modfrstv
    USE comgeom, ONLY : rlatu, rlatv, xprimu
    USE dimens_m, ONLY : iim, jjm
    use inifgn_m, only: inifgn
    use nr_util, only: pi
    USE serre, ONLY : grossismx

    ! Local:
    REAL dlatu(jjm)
    REAL rlamda(2: iim), eignvl(iim)

    REAL lamdamax, cof
    INTEGER i, j, modemax, imx, k, kf
    REAL dymin, colat0
    REAL eignft(iim, iim), coff

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

    print *, "Call sequence information: inifilr"

    CALL inifgn(eignvl)

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

    ! 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

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

    dymin = dlatu(1)
    DO j = 2, jjm
       dymin = min(dymin, dlatu(j))
    END DO

    colat0 = min(0.5, dymin / minval(xprimu(:iim)))

    PRINT *, 'colat0 = ', colat0

    lamdamax = iim / (pi * colat0 / grossismx)
    rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))

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

    ! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv

    modemax = iim
    imx = iim

    PRINT *, 'TRUNCATION AT ', imx

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

       IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &
            .and. rlamda(imx) * cos(rlatu(jjm - j + 2)) < 1.) &
            jfiltsu = jjm - j + 2
    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.) exit
       end DO
       if (k == modemax + 1) cycle
       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.) exit
       end DO
       if (k == modemax + 1) cycle
       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.) exit
       end DO
       if (k == modemax + 1) cycle
       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.) exit
       end DO
       if (k == modemax + 1) cycle
       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
          IF (i < modfrstu(j)) then
             coff = 0.
          else
             coff = coefilu(i, j)
          end IF
          eignft(i, :) = eignfnv(:, i)*coff
       END DO
       matriceun(:, :, j) = matmul(eignfnv, eignft)
    END DO

    DO j = jfiltsu, jjm
       DO i = 1, iim
          IF (i < modfrstu(j)) then
             coff = 0.
          else
             coff = coefilu(i, j)
          end IF
          eignft(i, :) = eignfnv(:, i) * coff
       END DO
       matriceus(:, :, j - jfiltsu + 1) = matmul(eignfnv, eignft)
    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
          IF (i < modfrstv(j)) then
             coff = 0.
          else
             coff = coefilv(i, j)
          end IF
          eignft(i, :) = eignfnu(:, i)*coff
       END DO
       matricevn(:, :, j) = matmul(eignfnu, eignft)
    END DO

    DO j = jfiltsv, jjm
       DO i = 1, iim
          IF (i < modfrstv(j)) then
             coff = 0.
          else
             coff = coefilv(i, j)
          end IF
          eignft(i, :) = eignfnu(:, i)*coff
       END DO
       matricevs(:, :, j-jfiltsv+1) = matmul(eignfnu, eignft)
    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
          IF (i < modfrstu(j)) then
             coff = 0.
          else
             coff = coefilu(i, j)/(1.+coefilu(i, j))
          end IF
          eignft(i, :) = eignfnv(:, i)*coff
       END DO
       matrinvn(:, :, j) = matmul(eignfnv, eignft)
    END DO

    DO j = jfiltsu, jjm
       DO i = 1, iim
          IF (i < modfrstu(j)) then
             coff = 0.
          else
             coff = coefilu(i, j)/(1.+coefilu(i, j))
          end IF
          eignft(i, :) = eignfnv(:, i)*coff
       END DO
       matrinvs(:, :, j-jfiltsu+1) = matmul(eignfnv, eignft)
    END DO

334 FORMAT (1X, 24I3)

  END SUBROUTINE inifilr

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