Sources/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	guez	54	module inifilr_m
2	guez	3
3	guez	54	use dimens_m, only: iim
4	guez	3
5	guez	32	IMPLICIT NONE
6	guez	3
7	guez	54	INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv
8			INTEGER, PARAMETER:: nfilun=3, nfilus=2, nfilvn=2, nfilvs=2
9	guez	3
10	guez	54	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	guez	3
14	guez	54	private iim, nfilun, nfilus, nfilvn, nfilvs
15	guez	3
16	guez	54	contains
17	guez	3
18	guez	54	SUBROUTINE inifilr
19	guez	3
20	guez	54	! From filtrez/inifilr.F, version 1.1.1.1 2004/05/19 12:53:09
21			! H. Upadhyaya, O. Sharma
22	guez	3
23	guez	54	! 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	guez	3
35	guez	54	! The modes are filtered from modfrst to modemax.
36	guez	3
37	guez	113	USE coefils, ONLY : coefilu, coefilu2, coefilv, coefilv2, eignfnu, &
38			eignfnv, modfrstu, modfrstv
39			USE comgeom, ONLY : rlatu, rlatv, xprimu
40	guez	54	USE dimens_m, ONLY : iim, jjm
41	guez	113	use inifgn_m, only: inifgn
42	guez	54	use nr_util, only: pi
43	guez	113	USE serre, ONLY : grossismx
44	guez	3
45	guez	54	! Local:
46	guez	132	REAL dlatu(jjm)
47	guez	54	REAL rlamda(2: iim), eignvl(iim)
48	guez	3
49	guez	54	REAL lamdamax, cof
50			INTEGER i, j, modemax, imx, k, kf
51	guez	132	REAL dymin, colat0
52	guez	54	REAL eignft(iim, iim), coff
53	guez	3
54	guez	54	!-----------------------------------------------------------
55	guez	3
56	guez	54	print *, "Call sequence information: inifilr"
57	guez	3
58	guez	54	CALL inifgn(eignvl)
59	guez	3
60	guez	54	PRINT *, 'EIGNVL '
61			PRINT "(1X, 5E13.6)", eignvl
62	guez	3
63	guez	54	! 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	guez	3
71	guez	54	! Calcul de colat0
72	guez	3
73	guez	54	DO j = 1, jjm
74			dlatu(j) = rlatu(j) - rlatu(j+1)
75			END DO
76	guez	3
77	guez	54	dymin = dlatu(1)
78			DO j = 2, jjm
79			dymin = min(dymin, dlatu(j))
80			END DO
81	guez	3
82	guez	132	colat0 = min(0.5, dymin / minval(xprimu(:iim)))
83	guez	3
84	guez	54	PRINT *, 'colat0 = ', colat0
85	guez	3
86	guez	113	lamdamax = iim / (pi * colat0 / grossismx)
87	guez	54	rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))
88	guez	3
89	guez	54	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	guez	3
98	guez	54	! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
99	guez	3
100	guez	54	modemax = iim
101			imx = iim
102	guez	3
103	guez	54	PRINT *, 'TRUNCATION AT ', imx
104	guez	3
105	guez	54	DO j = 2, jjm / 2 + 1
106			IF (cos(rlatu(j)) / colat0 < 1. &
107			.and. rlamda(imx) * cos(rlatu(j)) < 1.) jfiltnu = j
108	guez	3
109	guez	54	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	guez	3
114	guez	54	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	guez	3
120	guez	54	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	guez	3
126	guez	54	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	guez	3
132	guez	54	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	guez	3
138	guez	54	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	guez	3
144	guez	54	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	guez	3
150	guez	54	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
151			jfiltsu
152	guez	32
153	guez	54	! Determination de coefilu, coefilv, n=modfrstu, modfrstv
154	guez	32
155	guez	54	DO j = 1, jjm
156			modfrstu(j) = iim
157			modfrstv(j) = iim
158			END DO
159	guez	32
160	guez	54	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	guez	32
168	guez	54	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	guez	32
176	guez	54	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	guez	32
184	guez	54	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	guez	32
192	guez	54	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	guez	32
200	guez	54	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	guez	32
208	guez	54	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	guez	32
216	guez	54	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	guez	32
224	guez	54	IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN
225			IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
226			IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
227	guez	32
228	guez	54	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
229			jfiltsu
230			END IF
231	guez	32
232	guez	54	PRINT *, 'Modes premiers v '
233			PRINT 334, modfrstv
234			PRINT *, 'Modes premiers u '
235			PRINT 334, modfrstu
236	guez	32
237	guez	54	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	guez	32
310	guez	54	! Calcul de la matrice filtre 'matriceu' pour les champs situes
311			! sur la grille scalaire
312	guez	32
313	guez	54	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	guez	32
325	guez	54	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	guez	32
337	guez	54	! Calcul de la matrice filtre 'matricev' pour les champs situes
338			! sur la grille de V ou de Z
339	guez	32
340	guez	54	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	guez	32
352	guez	54	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	guez	32
364	guez	54	! Calcul de la matrice filtre 'matrinv' pour les champs situes
365			! sur la grille scalaire , pour le filtre inverse
366	guez	32
367	guez	54	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	guez	32
379	guez	54	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	guez	32
391	guez	54	334 FORMAT (1X, 24I3)
392	guez	32
393	guez	54	END SUBROUTINE inifilr
394	guez	32
395	guez	54	end module inifilr_m