libf/filtrez/inifilr.f90

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 dimens_m, ONLY : iim, jjm
    use conf_gcm_m, ONLY : fxyhypb, ysinus
    USE comgeom, ONLY : rlatu, rlatv, xprimu
    use nr_util, only: pi
    USE serre, ONLY : alphax
    USE coefils, ONLY : coefilu, coefilu2, coefilv, coefilv2, eignfnu, &
         eignfnv, modfrstu, modfrstv

    ! 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
    EXTERNAL inifgn

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

    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)

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

    PRINT *, 'colat0 = ', colat0
    PRINT *, 'alphax = ', alphax

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

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