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 dlonu(iim), dlatu(jjm)
    REAL rlamda(2: iim), eignvl(iim)

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

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

    print *, "Call sequence information: inifilr"

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

    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

    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)

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