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	module inifilr_m
2
3	IMPLICIT NONE
4
5	! North:
6
7	INTEGER jfiltnu, jfiltnv
8	! index of the last scalar line filtered in northern hemisphere
9
10	real, allocatable:: matriceun(:, :, :) ! (iim, iim, 2:jfiltnu)
11	! matrice filtre pour les champs situes sur la grille scalaire
12
13	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	real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)
18	! matrice filtre pour les champs situes sur la grille de V ou de Z
19
20	! South:
21
22	integer jfiltsu, jfiltsv
23	! index of the first line filtered in southern hemisphere
24
25	real, allocatable:: matriceus(:, :, :) ! (iim, iim, jfiltsu:jjm)
26	! matrice filtre pour les champs situes sur la grille scalaire
27
28	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	real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
33	! matrice filtre pour les champs situes sur la grille de V ou de Z
34
35	contains
36
37	SUBROUTINE inifilr
38
39	! From filtrez/inifilr.F, version 1.1.1.1, 2004/05/19 12:53:09
40	! H. Upadhyaya, O. Sharma
41
42	! 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
49	USE dimens_m, ONLY : iim, jjm
50	USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
51	use inifgn_m, only: inifgn
52	use jumble, only: new_unit
53	use nr_util, only: pi
54
55	! Local:
56
57	REAL dlatu(jjm)
58	REAL rlamda(2:iim)
59	real eignvl(iim) ! eigenvalues sorted in descending order (<= 0)
60	INTEGER i, j, unit
61	REAL colat0 ! > 0
62	REAL eignft(iim, iim)
63
64	real eignfnu(iim, iim), eignfnv(iim, iim)
65	! eigenvectors of the discrete second derivative with respect to longitude
66
67	! Filtering coefficients (lamda_max * cos(rlat) / lamda):
68	real coefil(iim)
69
70	! Index of the mode from where modes are filtered:
71	integer, allocatable:: modfrstnu(:) ! (2:jfiltnu)
72	integer, allocatable:: modfrstsu(:) ! (jfiltsu:jjm)
73	integer, allocatable:: modfrstnv(:) ! (jfiltnv)
74	integer, allocatable:: modfrstsv(:) ! (jfiltsv:jjm)
75
76	!-----------------------------------------------------------
77
78	print *, "Call sequence information: inifilr"
79
80	CALL inifgn(eignvl, eignfnu, eignfnv)
81
82	! Calcul de colat0
83	forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
84	colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
85	PRINT *, 'colat0 = ', colat0
86
87	rlamda = iim / (pi * colat0 / grossismx) / sqrt(- eignvl(2: iim))
88
89	! D\'etermination de jfilt[ns][uv] :
90
91	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
97	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
103	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
109	if (cos(rlatv(jfiltnv)) >= colat0 &
110	.or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) then
111	! {jfiltnv == 1}
112	PRINT *, 'Could not find jfiltnv.'
113	STOP 1
114	END IF
115
116	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
122	IF (cos(rlatv(jfiltsv)) >= colat0 &
123	.or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) THEN
124	! {jfiltsv == jjm}
125	PRINT *, 'Could not find jfiltsv.'
126	STOP 1
127	END IF
128
129	PRINT *, 'jfiltnu =', jfiltnu
130	PRINT *, 'jfiltsu =', jfiltsu
131	PRINT *, 'jfiltnv =', jfiltnv
132	PRINT *, 'jfiltsv =', jfiltsv
133
134	! D\'etermination de modfrst[ns][uv] :
135
136	allocate(modfrstnu(2:jfiltnu), modfrstsu(jfiltsu:jjm))
137	allocate(modfrstnv(jfiltnv), modfrstsv(jfiltsv:jjm))
138
139	DO j = 2, jfiltnu
140	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	END DO
146
147	DO j = 1, jfiltnv
148	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	end DO
154
155	DO j = jfiltsu, jjm
156	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	end DO
162
163	DO j = jfiltsv, jjm
164	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	END DO
170
171	call new_unit(unit)
172
173	open(unit, file = "inifilr_out.txt", status = "replace", action = "write")
174	write(unit, fmt = *) '"EIGNVL"', eignvl
175	close(unit)
176
177	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	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
218	! Calcul de matriceu et matrinv
219
220	DO j = 2, jfiltnu
221	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	END DO
239
240	DO j = jfiltsu, jjm
241	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	END DO
259
260	! Calcul de matricev
261
262	DO j = 1, jfiltnv
263	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	END DO
275
276	DO j = jfiltsv, jjm
277	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
282	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	END DO
289
290	END SUBROUTINE inifilr
291
292	end module inifilr_m