trunk/filtrez/inifilr.f

module inifilr_m

  IMPLICIT NONE

  ! North:

  INTEGER jfiltnu, jfiltnv
  ! index of the last scalar line filtered in northern hemisphere

  real, allocatable:: matriceun(:, :, :) ! (iim, iim, 2:jfiltnu)
  ! matrice filtre pour les champs situes sur la grille scalaire

  real, allocatable:: matrinvn(:, :, :) ! (iim, iim, 2:jfiltnu)
  ! matrice filtre pour les champs situes sur la grille scalaire, pour
  ! le filtre inverse

  real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)
  ! matrice filtre pour les champs situes sur la grille de V ou de Z

  ! South:

  integer jfiltsu, jfiltsv
  ! index of the first line filtered in southern hemisphere

  real, allocatable:: matriceus(:, :, :) ! (iim, iim, jfiltsu:jjm)
  ! matrice filtre pour les champs situes sur la grille scalaire

  real, allocatable:: matrinvs(:, :, :) ! (iim, iim, jfiltsu:jjm)
  ! matrice filtre pour les champs situes sur la grille scalaire, pour
  ! le filtre inverse

  real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
  ! matrice filtre pour les champs situes sur la grille de V ou de Z

contains

  SUBROUTINE inifilr

    ! From filtrez/inifilr.F, version 1.1.1.1, 2004/05/19 12:53:09
    ! H. Upadhyaya, O. Sharma

    ! This procedure computes the filtering coefficients for scalar
    ! lines and meridional wind v lines. The modes are filtered from
    ! modfrst to iim. 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.

    USE dimens_m, ONLY : iim, jjm
    USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
    use inifgn_m, only: inifgn
    use jumble, only: new_unit
    use nr_util, only: pi

    ! Local:

    REAL dlatu(jjm)
    REAL rlamda(2:iim)
    real eignvl(iim) ! eigenvalues sorted in descending order (<= 0)
    INTEGER i, j, unit
    REAL colat0 ! > 0
    REAL eignft(iim, iim)

    real eignfnu(iim, iim), eignfnv(iim, iim)
    ! eigenvectors of the discrete second derivative with respect to longitude

    ! Filtering coefficients (lamda_max * cos(rlat) / lamda):
    real coefil(iim)

    ! Index of the mode from where modes are filtered:
    integer, allocatable:: modfrstnu(:) ! (2:jfiltnu)
    integer, allocatable:: modfrstsu(:) ! (jfiltsu:jjm)
    integer, allocatable:: modfrstnv(:) ! (jfiltnv)
    integer, allocatable:: modfrstsv(:) ! (jfiltsv:jjm)

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

    print *, "Call sequence information: inifilr"

    CALL inifgn(eignvl, eignfnu, eignfnv)

    ! Calcul de colat0
    forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
    colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
    PRINT *, 'colat0 = ', colat0

    rlamda = iim / (pi * colat0 / grossismx) / sqrt(- eignvl(2: iim))

    ! D\'etermination de jfilt[ns][uv] :

    jfiltnu = (jjm + 1) / 2
    do while (cos(rlatu(jfiltnu)) >= colat0 &
         .or. rlamda(iim) * cos(rlatu(jfiltnu)) >= 1.)
       jfiltnu = jfiltnu - 1
    end do

    jfiltsu = jjm / 2 + 2
    do while (cos(rlatu(jfiltsu)) >= colat0 &
         .or. rlamda(iim) * cos(rlatu(jfiltsu)) >= 1.)
       jfiltsu = jfiltsu + 1
    end do

    jfiltnv = jjm / 2
    do while ((cos(rlatv(jfiltnv)) >= colat0 &
         .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) .and. jfiltnv >= 2)
       jfiltnv = jfiltnv - 1
    end do

    if (cos(rlatv(jfiltnv)) >= colat0 &
         .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) then
       ! {jfiltnv == 1}
       PRINT *, 'Could not find jfiltnv.'
       STOP 1
    END IF

    jfiltsv = (jjm + 1)/ 2 + 1
    do while ((cos(rlatv(jfiltsv)) >= colat0 &
         .or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) .and. jfiltsv <= jjm - 1)
       jfiltsv = jfiltsv + 1
    end do

    IF (cos(rlatv(jfiltsv)) >= colat0 &
         .or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) THEN
       ! {jfiltsv == jjm}
       PRINT *, 'Could not find jfiltsv.'
       STOP 1
    END IF

    PRINT *, 'jfiltnu =', jfiltnu
    PRINT *, 'jfiltsu =', jfiltsu
    PRINT *, 'jfiltnv =', jfiltnv
    PRINT *, 'jfiltsv =', jfiltsv

    ! D\'etermination de modfrst[ns][uv] :

    allocate(modfrstnu(2:jfiltnu), modfrstsu(jfiltsu:jjm))
    allocate(modfrstnv(jfiltnv), modfrstsv(jfiltsv:jjm))

    DO j = 2, jfiltnu
       modfrstnu(j) = 2
       do while (rlamda(modfrstnu(j)) * cos(rlatu(j)) >= 1. &
            .and. modfrstnu(j) <= iim - 1)
          modfrstnu(j) = modfrstnu(j) + 1
       end do
    END DO

    DO j = 1, jfiltnv
       modfrstnv(j) = 2
       do while (rlamda(modfrstnv(j)) * cos(rlatv(j)) >= 1. &
            .and. modfrstnv(j) <= iim - 1)
          modfrstnv(j) = modfrstnv(j) + 1
       end do
    end DO

    DO j = jfiltsu, jjm
       modfrstsu(j) = 2
       do while (rlamda(modfrstsu(j)) * cos(rlatu(j)) >= 1. &
            .and. modfrstsu(j) <= iim - 1)
          modfrstsu(j) = modfrstsu(j) + 1
       end do
    end DO

    DO j = jfiltsv, jjm
       modfrstsv(j) = 2
       do while (rlamda(modfrstsv(j)) * cos(rlatv(j)) >= 1. &
            .and. modfrstsv(j) <= iim - 1)
          modfrstsv(j) = modfrstsv(j) + 1
       end do
    END DO

    call new_unit(unit)

    open(unit, file = "inifilr_out.txt", status = "replace", action = "write") 
    write(unit, fmt = *) '"EIGNVL"', eignvl
    close(unit)

    open(unit, file = "modfrstnu.csv", status = "replace", action = "write") 
    write(unit, fmt = *) '"rlatu (degrees north)" modfrstnu ' &
         // '"rlamda(modfrstnu) * cos(rlatu) < 1"'
    DO j = 2, jfiltnu
       write(unit, fmt = *) rlatu(j) / pi * 180., modfrstnu(j), &
            rlamda(modfrstnu(j)) * cos(rlatu(j)) < 1
    end DO
    close(unit)

    open(unit, file = "modfrstnv.csv", status = "replace", action = "write") 
    write(unit, fmt = *) '"rlatv (degrees north)" modfrstnv ' &
         // '"rlamda(modfrstnv) * cos(rlatv) < 1"'
    DO j = 1, jfiltnv
       write(unit, fmt = *) rlatv(j) / pi * 180., modfrstnv(j), &
            rlamda(modfrstnv(j)) * cos(rlatv(j)) < 1
    end DO
    close(unit)

    open(unit, file = "modfrstsu.csv", status = "replace", action = "write") 
     write(unit, fmt = *) '"rlatu (degrees north)" modfrstsu ' &
         // '"rlamda(modfrstsu) * cos(rlatu) < 1"'
   DO j = jfiltsu, jjm
       write(unit, fmt = *) rlatu(j) / pi * 180., modfrstsu(j), &
            rlamda(modfrstsu(j)) * cos(rlatu(j)) < 1
    end DO
    close(unit)

    open(unit, file = "modfrstsv.csv", status = "replace", action = "write") 
    write(unit, fmt = *) '"rlatv (degrees north)" modfrstsv ' &
         // '"rlamda(modfrstsv) * cos(rlatv) < 1"'
    DO j = jfiltsv, jjm
       write(unit, fmt = *) rlatv(j) / pi * 180., modfrstsv(j), &
            rlamda(modfrstsv(j)) * cos(rlatv(j)) < 1
    end DO
    close(unit)

    allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
    allocate(matricevn(iim, iim, jfiltnv))
    allocate(matricevs(iim, iim, jfiltsv:jjm))
    allocate(matriceus(iim, iim, jfiltsu:jjm), matrinvs(iim, iim, jfiltsu:jjm))

    ! Calcul de matriceu et matrinv

    DO j = 2, jfiltnu
       if (rlamda(modfrstnu(j)) * cos(rlatu(j)) < 1.) then
          DO i = modfrstnu(j), iim
             coefil(i) = rlamda(i) * cos(rlatu(j)) - 1.
          end DO

          eignft(:modfrstnu(j) - 1, :) = 0.

          forall (i = modfrstnu(j):iim) eignft(i, :) = eignfnv(:, i) * coefil(i)
          matriceun(:, :, j) = matmul(eignfnv, eignft)

          forall (i = modfrstnu(j):iim) eignft(i, :) = eignfnv(:, i) &
               * coefil(i) / (1. + coefil(i))
          matrinvn(:, :, j) = matmul(eignfnv, eignft)
       else
          matriceun(:, :, j) = 0.
          matrinvn(:, :, j) = 0.
       end if
    END DO

    DO j = jfiltsu, jjm
       if (rlamda(modfrstsu(j)) * cos(rlatu(j)) < 1.) then
          DO i = modfrstsu(j), iim
             coefil(i) = rlamda(i) * cos(rlatu(j)) - 1.
          end DO

          eignft(:modfrstsu(j) - 1, :) = 0.

          forall (i = modfrstsu(j):iim) eignft(i, :) = eignfnv(:, i) * coefil(i)
          matriceus(:, :, j) = matmul(eignfnv, eignft)

          forall (i = modfrstsu(j):iim) eignft(i, :) = eignfnv(:, i) &
               * coefil(i) / (1. + coefil(i))
          matrinvs(:, :, j) = matmul(eignfnv, eignft)
       else
          matriceus(:, :, j) = 0.
          matrinvs(:, :, j) = 0.
       end if
    END DO

    ! Calcul de matricev

    DO j = 1, jfiltnv
       if (rlamda(modfrstnv(j)) * cos(rlatv(j)) < 1.) then
          DO i = modfrstnv(j), iim
             coefil(i) = rlamda(i) * cos(rlatv(j)) - 1.
          end DO

          eignft(:modfrstnv(j) - 1, :) = 0.
          forall (i = modfrstnv(j):iim) eignft(i, :) = eignfnu(:, i) * coefil(i)
          matricevn(:, :, j) = matmul(eignfnu, eignft)
       else
          matricevn(:, :, j) = 0.
       end if
    END DO

    DO j = jfiltsv, jjm
       if (rlamda(modfrstsv(j)) * cos(rlatv(j)) < 1.) then
          DO i = modfrstsv(j), iim
             coefil(i) = rlamda(i) * cos(rlatv(j)) - 1.
          end DO

          eignft(:modfrstsv(j) - 1, :) = 0.
          forall (i = modfrstsv(j):iim) eignft(i, :) = eignfnu(:, i) * coefil(i)
          matricevs(:, :, j) = matmul(eignfnu, eignft)
       else
          matricevs(:, :, j) = 0.
       end if
    END DO

  END SUBROUTINE inifilr

end module inifilr_m
1	guez	54	module inifilr_m
2	guez	3
3	guez	32	IMPLICIT NONE
4	guez	3
5	guez	136	! North:
6	guez	156
7	guez	164	INTEGER jfiltnu, jfiltnv
8			! index of the last scalar line filtered in northern hemisphere
9
10	guez	165	real, allocatable:: matriceun(:, :, :) ! (iim, iim, 2:jfiltnu)
11			! matrice filtre pour les champs situes sur la grille scalaire
12	guez	3
13	guez	165	real, allocatable:: matrinvn(:, :, :) ! (iim, iim, 2:jfiltnu)
14			! matrice filtre pour les champs situes sur la grille scalaire, pour
15			! le filtre inverse
16
17	guez	136	real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)
18	guez	165	! matrice filtre pour les champs situes sur la grille de V ou de Z
19	guez	3
20	guez	136	! South:
21	guez	156
22	guez	164	integer jfiltsu, jfiltsv
23			! index of the first line filtered in southern hemisphere
24
25	guez	165	real, allocatable:: matriceus(:, :, :) ! (iim, iim, jfiltsu:jjm)
26			! matrice filtre pour les champs situes sur la grille scalaire
27	guez	136
28	guez	165	real, allocatable:: matrinvs(:, :, :) ! (iim, iim, jfiltsu:jjm)
29			! matrice filtre pour les champs situes sur la grille scalaire, pour
30			! le filtre inverse
31
32	guez	136	real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
33	guez	165	! matrice filtre pour les champs situes sur la grille de V ou de Z
34	guez	136
35	guez	54	contains
36	guez	3
37	guez	54	SUBROUTINE inifilr
38	guez	3
39	guez	164	! From filtrez/inifilr.F, version 1.1.1.1, 2004/05/19 12:53:09
40	guez	54	! H. Upadhyaya, O. Sharma
41	guez	3
42	guez	164	! This procedure computes the filtering coefficients for scalar
43			! lines and meridional wind v lines. The modes are filtered from
44			! modfrst to iim. We filter all those latitude lines where coefil
45			! < 1. No filtering at poles. colat0 is to be used when alpha
46			! (stretching coefficient) is set equal to zero for the regular
47			! grid case.
48	guez	3
49	guez	54	USE dimens_m, ONLY : iim, jjm
50	guez	139	USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
51	guez	154	use inifgn_m, only: inifgn
52	guez	143	use jumble, only: new_unit
53	guez	54	use nr_util, only: pi
54	guez	3
55	guez	54	! Local:
56	guez	164
57	guez	132	REAL dlatu(jjm)
58	guez	165	REAL rlamda(2:iim)
59	guez	164	real eignvl(iim) ! eigenvalues sorted in descending order (<= 0)
60			INTEGER i, j, unit
61	guez	151	REAL colat0 ! > 0
62	guez	164	REAL eignft(iim, iim)
63	guez	3
64	guez	154	real eignfnu(iim, iim), eignfnv(iim, iim)
65	guez	164	! eigenvectors of the discrete second derivative with respect to longitude
66	guez	154
67	guez	140	! Filtering coefficients (lamda_max * cos(rlat) / lamda):
68	guez	165	real coefil(iim)
69	guez	140
70	guez	151	! Index of the mode from where modes are filtered:
71	guez	164	integer, allocatable:: modfrstnu(:) ! (2:jfiltnu)
72			integer, allocatable:: modfrstsu(:) ! (jfiltsu:jjm)
73			integer, allocatable:: modfrstnv(:) ! (jfiltnv)
74			integer, allocatable:: modfrstsv(:) ! (jfiltsv:jjm)
75	guez	140
76	guez	54	!-----------------------------------------------------------
77	guez	3
78	guez	54	print *, "Call sequence information: inifilr"
79	guez	3
80	guez	154	CALL inifgn(eignvl, eignfnu, eignfnv)
81	guez	3
82	guez	54	! Calcul de colat0
83	guez	143	forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
84			colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
85	guez	54	PRINT *, 'colat0 = ', colat0
86	guez	3
87	guez	164	rlamda = iim / (pi * colat0 / grossismx) / sqrt(- eignvl(2: iim))
88	guez	3
89	guez	165	! D\'etermination de jfilt[ns][uv] :
90	guez	3
91	guez	151	jfiltnu = (jjm + 1) / 2
92			do while (cos(rlatu(jfiltnu)) >= colat0 &
93			.or. rlamda(iim) * cos(rlatu(jfiltnu)) >= 1.)
94			jfiltnu = jfiltnu - 1
95			end do
96	guez	3
97	guez	151	jfiltsu = jjm / 2 + 2
98			do while (cos(rlatu(jfiltsu)) >= colat0 &
99			.or. rlamda(iim) * cos(rlatu(jfiltsu)) >= 1.)
100			jfiltsu = jfiltsu + 1
101			end do
102	guez	3
103	guez	151	jfiltnv = jjm / 2
104			do while ((cos(rlatv(jfiltnv)) >= colat0 &
105			.or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) .and. jfiltnv >= 2)
106			jfiltnv = jfiltnv - 1
107			end do
108	guez	3
109	guez	151	if (cos(rlatv(jfiltnv)) >= colat0 &
110			.or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) then
111			! {jfiltnv == 1}
112			PRINT *, 'Could not find jfiltnv.'
113	guez	54	STOP 1
114			END IF
115	guez	3
116	guez	151	jfiltsv = (jjm + 1)/ 2 + 1
117			do while ((cos(rlatv(jfiltsv)) >= colat0 &
118			.or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) .and. jfiltsv <= jjm - 1)
119			jfiltsv = jfiltsv + 1
120			end do
121	guez	3
122	guez	151	IF (cos(rlatv(jfiltsv)) >= colat0 &
123			.or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) THEN
124			! {jfiltsv == jjm}
125			PRINT *, 'Could not find jfiltsv.'
126	guez	54	STOP 1
127			END IF
128	guez	3
129	guez	151	PRINT *, 'jfiltnu =', jfiltnu
130			PRINT *, 'jfiltsu =', jfiltsu
131			PRINT *, 'jfiltnv =', jfiltnv
132			PRINT *, 'jfiltsv =', jfiltsv
133	guez	3
134	guez	165	! D\'etermination de modfrst[ns][uv] :
135	guez	32
136	guez	151	allocate(modfrstnu(2:jfiltnu), modfrstsu(jfiltsu:jjm))
137			allocate(modfrstnv(jfiltnv), modfrstsv(jfiltsv:jjm))
138	guez	32
139	guez	54	DO j = 2, jfiltnu
140	guez	151	modfrstnu(j) = 2
141			do while (rlamda(modfrstnu(j)) * cos(rlatu(j)) >= 1. &
142			.and. modfrstnu(j) <= iim - 1)
143			modfrstnu(j) = modfrstnu(j) + 1
144			end do
145	guez	54	END DO
146	guez	32
147	guez	54	DO j = 1, jfiltnv
148	guez	151	modfrstnv(j) = 2
149			do while (rlamda(modfrstnv(j)) * cos(rlatv(j)) >= 1. &
150			.and. modfrstnv(j) <= iim - 1)
151			modfrstnv(j) = modfrstnv(j) + 1
152			end do
153	guez	54	end DO
154	guez	32
155	guez	54	DO j = jfiltsu, jjm
156	guez	151	modfrstsu(j) = 2
157			do while (rlamda(modfrstsu(j)) * cos(rlatu(j)) >= 1. &
158			.and. modfrstsu(j) <= iim - 1)
159			modfrstsu(j) = modfrstsu(j) + 1
160			end do
161	guez	54	end DO
162	guez	32
163	guez	54	DO j = jfiltsv, jjm
164	guez	151	modfrstsv(j) = 2
165			do while (rlamda(modfrstsv(j)) * cos(rlatv(j)) >= 1. &
166			.and. modfrstsv(j) <= iim - 1)
167			modfrstsv(j) = modfrstsv(j) + 1
168			end do
169	guez	54	END DO
170	guez	32
171	guez	151	call new_unit(unit)
172	guez	165
173	guez	151	open(unit, file = "inifilr_out.txt", status = "replace", action = "write")
174			write(unit, fmt = *) '"EIGNVL"', eignvl
175			close(unit)
176	guez	32
177	guez	165	open(unit, file = "modfrstnu.csv", status = "replace", action = "write")
178			write(unit, fmt = *) '"rlatu (degrees north)" modfrstnu ' &
179			// '"rlamda(modfrstnu) * cos(rlatu) < 1"'
180			DO j = 2, jfiltnu
181			write(unit, fmt = ) rlatu(j) / pi 180., modfrstnu(j), &
182			rlamda(modfrstnu(j)) * cos(rlatu(j)) < 1
183			end DO
184			close(unit)
185
186			open(unit, file = "modfrstnv.csv", status = "replace", action = "write")
187			write(unit, fmt = *) '"rlatv (degrees north)" modfrstnv ' &
188			// '"rlamda(modfrstnv) * cos(rlatv) < 1"'
189			DO j = 1, jfiltnv
190			write(unit, fmt = ) rlatv(j) / pi 180., modfrstnv(j), &
191			rlamda(modfrstnv(j)) * cos(rlatv(j)) < 1
192			end DO
193			close(unit)
194
195			open(unit, file = "modfrstsu.csv", status = "replace", action = "write")
196			write(unit, fmt = *) '"rlatu (degrees north)" modfrstsu ' &
197			// '"rlamda(modfrstsu) * cos(rlatu) < 1"'
198			DO j = jfiltsu, jjm
199			write(unit, fmt = ) rlatu(j) / pi 180., modfrstsu(j), &
200			rlamda(modfrstsu(j)) * cos(rlatu(j)) < 1
201			end DO
202			close(unit)
203
204			open(unit, file = "modfrstsv.csv", status = "replace", action = "write")
205			write(unit, fmt = *) '"rlatv (degrees north)" modfrstsv ' &
206			// '"rlamda(modfrstsv) * cos(rlatv) < 1"'
207			DO j = jfiltsv, jjm
208			write(unit, fmt = ) rlatv(j) / pi 180., modfrstsv(j), &
209			rlamda(modfrstsv(j)) * cos(rlatv(j)) < 1
210			end DO
211			close(unit)
212
213	guez	136	allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
214			allocate(matricevn(iim, iim, jfiltnv))
215			allocate(matricevs(iim, iim, jfiltsv:jjm))
216			allocate(matriceus(iim, iim, jfiltsu:jjm), matrinvs(iim, iim, jfiltsu:jjm))
217	guez	32
218	guez	165	! Calcul de matriceu et matrinv
219	guez	32
220	guez	54	DO j = 2, jfiltnu
221	guez	165	if (rlamda(modfrstnu(j)) * cos(rlatu(j)) < 1.) then
222			DO i = modfrstnu(j), iim
223			coefil(i) = rlamda(i) * cos(rlatu(j)) - 1.
224			end DO
225
226			eignft(:modfrstnu(j) - 1, :) = 0.
227
228			forall (i = modfrstnu(j):iim) eignft(i, :) = eignfnv(:, i) * coefil(i)
229			matriceun(:, :, j) = matmul(eignfnv, eignft)
230
231			forall (i = modfrstnu(j):iim) eignft(i, :) = eignfnv(:, i) &
232			* coefil(i) / (1. + coefil(i))
233			matrinvn(:, :, j) = matmul(eignfnv, eignft)
234			else
235			matriceun(:, :, j) = 0.
236			matrinvn(:, :, j) = 0.
237			end if
238	guez	54	END DO
239	guez	32
240	guez	54	DO j = jfiltsu, jjm
241	guez	165	if (rlamda(modfrstsu(j)) * cos(rlatu(j)) < 1.) then
242			DO i = modfrstsu(j), iim
243			coefil(i) = rlamda(i) * cos(rlatu(j)) - 1.
244			end DO
245
246			eignft(:modfrstsu(j) - 1, :) = 0.
247
248			forall (i = modfrstsu(j):iim) eignft(i, :) = eignfnv(:, i) * coefil(i)
249			matriceus(:, :, j) = matmul(eignfnv, eignft)
250
251			forall (i = modfrstsu(j):iim) eignft(i, :) = eignfnv(:, i) &
252			* coefil(i) / (1. + coefil(i))
253			matrinvs(:, :, j) = matmul(eignfnv, eignft)
254			else
255			matriceus(:, :, j) = 0.
256			matrinvs(:, :, j) = 0.
257			end if
258	guez	54	END DO
259	guez	32
260	guez	165	! Calcul de matricev
261	guez	32
262	guez	54	DO j = 1, jfiltnv
263	guez	165	if (rlamda(modfrstnv(j)) * cos(rlatv(j)) < 1.) then
264			DO i = modfrstnv(j), iim
265			coefil(i) = rlamda(i) * cos(rlatv(j)) - 1.
266			end DO
267
268			eignft(:modfrstnv(j) - 1, :) = 0.
269			forall (i = modfrstnv(j):iim) eignft(i, :) = eignfnu(:, i) * coefil(i)
270			matricevn(:, :, j) = matmul(eignfnu, eignft)
271			else
272			matricevn(:, :, j) = 0.
273			end if
274	guez	54	END DO
275	guez	32
276	guez	54	DO j = jfiltsv, jjm
277	guez	165	if (rlamda(modfrstsv(j)) * cos(rlatv(j)) < 1.) then
278			DO i = modfrstsv(j), iim
279			coefil(i) = rlamda(i) * cos(rlatv(j)) - 1.
280			end DO
281	guez	32
282	guez	165	eignft(:modfrstsv(j) - 1, :) = 0.
283			forall (i = modfrstsv(j):iim) eignft(i, :) = eignfnu(:, i) * coefil(i)
284			matricevs(:, :, j) = matmul(eignfnu, eignft)
285			else
286			matricevs(:, :, j) = 0.
287			end if
288	guez	54	END DO
289	guez	32
290	guez	54	END SUBROUTINE inifilr
291	guez	32
292	guez	54	end module inifilr_m