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	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			REAL dlonu(iim), dlatu(jjm)
47			REAL rlamda(2: iim), eignvl(iim)
48	guez	3
49	guez	54	REAL lamdamax, cof
50			INTEGER i, j, modemax, imx, k, kf
51			REAL dymin, dxmin, colat0
52			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	DO i = 1, iim
59			dlonu(i) = xprimu(i)
60			END DO
61	guez	3
62	guez	54	CALL inifgn(eignvl)
63	guez	3
64	guez	54	PRINT *, 'EIGNVL '
65			PRINT "(1X, 5E13.6)", eignvl
66	guez	3
67	guez	54	! 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	guez	3
75	guez	54	! Calcul de colat0
76	guez	3
77	guez	54	DO j = 1, jjm
78			dlatu(j) = rlatu(j) - rlatu(j+1)
79			END DO
80	guez	3
81	guez	54	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	guez	3
90	guez	54	colat0 = min(0.5, dymin/dxmin)
91	guez	3
92	guez	54	PRINT *, 'colat0 = ', colat0
93	guez	3
94	guez	113	lamdamax = iim / (pi * colat0 / grossismx)
95	guez	54	rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))
96	guez	3
97	guez	54	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	guez	3
106	guez	54	! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
107	guez	3
108	guez	54	modemax = iim
109			imx = iim
110	guez	3
111	guez	54	PRINT *, 'TRUNCATION AT ', imx
112	guez	3
113	guez	54	DO j = 2, jjm / 2 + 1
114			IF (cos(rlatu(j)) / colat0 < 1. &
115			.and. rlamda(imx) * cos(rlatu(j)) < 1.) jfiltnu = j
116	guez	3
117	guez	54	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	guez	3
122	guez	54	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	guez	3
128	guez	54	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	guez	3
134	guez	54	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	guez	3
140	guez	54	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	guez	3
146	guez	54	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	guez	3
152	guez	54	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	guez	3
158	guez	54	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
159			jfiltsu
160	guez	32
161	guez	54	! Determination de coefilu, coefilv, n=modfrstu, modfrstv
162	guez	32
163	guez	54	DO j = 1, jjm
164			modfrstu(j) = iim
165			modfrstv(j) = iim
166			END DO
167	guez	32
168	guez	54	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	guez	32
176	guez	54	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	guez	32
184	guez	54	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	guez	32
192	guez	54	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	guez	32
200	guez	54	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	guez	32
208	guez	54	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	guez	32
216	guez	54	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	guez	32
224	guez	54	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	guez	32
232	guez	54	IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN
233			IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
234			IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
235	guez	32
236	guez	54	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
237			jfiltsu
238			END IF
239	guez	32
240	guez	54	PRINT *, 'Modes premiers v '
241			PRINT 334, modfrstv
242			PRINT *, 'Modes premiers u '
243			PRINT 334, modfrstu
244	guez	32
245	guez	54	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	guez	32
318	guez	54	! Calcul de la matrice filtre 'matriceu' pour les champs situes
319			! sur la grille scalaire
320	guez	32
321	guez	54	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	guez	32
333	guez	54	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	guez	32
345	guez	54	! Calcul de la matrice filtre 'matricev' pour les champs situes
346			! sur la grille de V ou de Z
347	guez	32
348	guez	54	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	guez	32
360	guez	54	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	guez	32
372	guez	54	! Calcul de la matrice filtre 'matrinv' pour les champs situes
373			! sur la grille scalaire , pour le filtre inverse
374	guez	32
375	guez	54	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	guez	32
387	guez	54	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	guez	32
399	guez	54	334 FORMAT (1X, 24I3)
400	guez	32
401	guez	54	END SUBROUTINE inifilr
402	guez	32
403	guez	54	end module inifilr_m