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 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	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	54	USE dimens_m, ONLY : iim, jjm
38	guez	57	use conf_gcm_m, ONLY : fxyhypb, ysinus
39	guez	54	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	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			EXTERNAL inifgn
54	guez	3
55	guez	54	!-----------------------------------------------------------
56	guez	3
57	guez	54	print *, "Call sequence information: inifilr"
58	guez	3
59	guez	54	DO i = 1, iim
60			dlonu(i) = xprimu(i)
61			END DO
62	guez	3
63	guez	54	CALL inifgn(eignvl)
64	guez	3
65	guez	54	PRINT *, 'EIGNVL '
66			PRINT "(1X, 5E13.6)", eignvl
67	guez	3
68	guez	54	! 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	guez	3
76	guez	54	! Calcul de colat0
77	guez	3
78	guez	54	DO j = 1, jjm
79			dlatu(j) = rlatu(j) - rlatu(j+1)
80			END DO
81	guez	3
82	guez	54	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	guez	3
91	guez	54	colat0 = min(0.5, dymin/dxmin)
92	guez	3
93	guez	54	IF (.NOT. fxyhypb .AND. ysinus) THEN
94			colat0 = 0.6
95			! À revoir pour ysinus
96			alphax = 0.
97			END IF
98	guez	3
99	guez	54	PRINT *, 'colat0 = ', colat0
100			PRINT *, 'alphax = ', alphax
101	guez	3
102	guez	54	IF (alphax == 1.) THEN
103			PRINT *, 'alphax doit etre < a 1. Corriger '
104			STOP 1
105			END IF
106	guez	3
107	guez	54	lamdamax = iim / (pi * colat0 * (1. - alphax))
108			rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))
109	guez	3
110	guez	54	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	guez	3
119	guez	54	! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
120	guez	3
121	guez	54	modemax = iim
122			imx = iim
123	guez	3
124	guez	54	PRINT *, 'TRUNCATION AT ', imx
125	guez	3
126	guez	54	DO j = 2, jjm / 2 + 1
127			IF (cos(rlatu(j)) / colat0 < 1. &
128			.and. rlamda(imx) * cos(rlatu(j)) < 1.) jfiltnu = j
129	guez	3
130	guez	54	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	guez	3
135	guez	54	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	guez	3
141	guez	54	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	guez	3
147	guez	54	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	guez	3
153	guez	54	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	guez	3
159	guez	54	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	guez	3
165	guez	54	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	guez	3
171	guez	54	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
172			jfiltsu
173	guez	32
174	guez	54	! Determination de coefilu, coefilv, n=modfrstu, modfrstv
175	guez	32
176	guez	54	DO j = 1, jjm
177			modfrstu(j) = iim
178			modfrstv(j) = iim
179			END DO
180	guez	32
181	guez	54	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	guez	32
189	guez	54	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	guez	32
197	guez	54	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	guez	32
205	guez	54	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	guez	32
213	guez	54	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	guez	32
221	guez	54	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	guez	32
229	guez	54	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	guez	32
237	guez	54	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	guez	32
245	guez	54	IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN
246			IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
247			IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
248	guez	32
249	guez	54	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
250			jfiltsu
251			END IF
252	guez	32
253	guez	54	PRINT *, 'Modes premiers v '
254			PRINT 334, modfrstv
255			PRINT *, 'Modes premiers u '
256			PRINT 334, modfrstu
257	guez	32
258	guez	54	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	guez	32
331	guez	54	! Calcul de la matrice filtre 'matriceu' pour les champs situes
332			! sur la grille scalaire
333	guez	32
334	guez	54	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	guez	32
346	guez	54	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	guez	32
358	guez	54	! Calcul de la matrice filtre 'matricev' pour les champs situes
359			! sur la grille de V ou de Z
360	guez	32
361	guez	54	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	guez	32
373	guez	54	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	guez	32
385	guez	54	! Calcul de la matrice filtre 'matrinv' pour les champs situes
386			! sur la grille scalaire , pour le filtre inverse
387	guez	32
388	guez	54	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	guez	32
400	guez	54	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	guez	32
412	guez	54	334 FORMAT (1X, 24I3)
413	guez	32
414	guez	54	END SUBROUTINE inifilr
415	guez	32
416	guez	54	end module inifilr_m