Sources/filtrez/inifilr.f

module inifilr_m

  IMPLICIT NONE

  INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv
  ! jfiltn index of the last scalar line filtered in NH
  ! jfilts index of the first line filtered in SH

  ! North:
  real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :) 
  ! (iim, iim, 2:jfiltnu)

  real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)

  ! South:
  real, allocatable:: matriceus(:, :, :), matrinvs(:, :, :)
  ! (iim, iim, jfiltsu:jjm)

  real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)

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. The modes are
    ! filtered from modfrst to iim.

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

    ! Local:
    REAL dlatu(jjm)
    REAL rlamda(2: iim)
    real eignvl(iim) ! eigenvalues
    REAL lamdamax, cof
    INTEGER i, j, k, kf
    REAL dymin, colat0
    REAL eignft(iim, iim), coff

    ! Filtering coefficients (lamda_max * cos(rlat) / lamda):
    real coefilu(iim, jjm), coefilv(iim, jjm)
    real coefilu2(iim, jjm), coefilv2(iim, jjm)

    integer modfrstu(jjm), modfrstv(jjm)
    ! index of the mode from where modes are filtered

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

    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

    PRINT *, 'TRUNCATION AT ', iim

    DO j = 2, jjm / 2 + 1
       IF (cos(rlatu(j)) / colat0 < 1. &
            .and. rlamda(iim) * cos(rlatu(j)) < 1.) jfiltnu = j

       IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &
            .and. rlamda(iim) * cos(rlatu(jjm - j + 2)) < 1.) &
            jfiltsu = jjm - j + 2
    END DO

    DO j = 1, jjm / 2
       IF (cos(rlatv(j)) / colat0 < 1. .and. rlamda(iim) * cos(rlatv(j)) < 1.) &
            jfiltnv = j

       IF (cos(rlatv(jjm - j + 1)) / colat0 < 1. .and. rlamda(iim) &
            * cos(rlatv(jjm - j + 1)) < 1.) jfiltsv = jjm - j + 1
    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, iim
          IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
       end DO
       if (k == iim + 1) cycle
       modfrstu(j) = k

       kf = modfrstu(j)
       DO k = kf, iim
          cof = rlamda(k) * cos(rlatu(j))
          coefilu(k, j) = cof - 1.
          coefilu2(k, j) = cof**2 - 1.
       end DO
    END DO

    DO j = 1, jfiltnv
       DO k = 2, iim
          IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
       end DO
       if (k == iim + 1) cycle
       modfrstv(j) = k

       kf = modfrstv(j)
       DO k = kf, iim
          cof = rlamda(k) * cos(rlatv(j))
          coefilv(k, j) = cof - 1.
          coefilv2(k, j) = cof**2 - 1.
       end DO
    end DO

    DO j = jfiltsu, jjm
       DO k = 2, iim
          IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
       end DO
       if (k == iim + 1) cycle
       modfrstu(j) = k

       kf = modfrstu(j)
       DO k = kf, iim
          cof = rlamda(k) * cos(rlatu(j))
          coefilu(k, j) = cof - 1.
          coefilu2(k, j) = cof**2 - 1.
       end DO
    end DO

    DO j = jfiltsv, jjm
       DO k = 2, iim
          IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
       end DO
       if (k == iim + 1) cycle
       modfrstv(j) = k

       kf = modfrstv(j)
       DO k = kf, iim
          cof = rlamda(k) * cos(rlatv(j))
          coefilv(k, j) = cof - 1.
          coefilv2(k, j) = cof**2 - 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

    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 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) = 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) = 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) = matmul(eignfnv, eignft)
    END DO

334 FORMAT (1X, 24I3)

  END SUBROUTINE inifilr

end module inifilr_m
1	module inifilr_m
2
3	IMPLICIT NONE
4
5	INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv
6	! jfiltn index of the last scalar line filtered in NH
7	! jfilts index of the first line filtered in SH
8
9	! North:
10	real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :)
11	! (iim, iim, 2:jfiltnu)
12
13	real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)
14
15	! South:
16	real, allocatable:: matriceus(:, :, :), matrinvs(:, :, :)
17	! (iim, iim, jfiltsu:jjm)
18
19	real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
20
21	contains
22
23	SUBROUTINE inifilr
24
25	! From filtrez/inifilr.F, version 1.1.1.1 2004/05/19 12:53:09
26	! H. Upadhyaya, O. Sharma
27
28	! This routine computes the eigenfunctions of the laplacian on the
29	! stretched grid, and the filtering coefficients. The modes are
30	! filtered from modfrst to iim.
31
32	USE dimens_m, ONLY : iim, jjm
33	USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
34	use inifgn_m, only: inifgn, eignfnu, eignfnv
35	use nr_util, only: pi
36
37	! Local:
38	REAL dlatu(jjm)
39	REAL rlamda(2: iim)
40	real eignvl(iim) ! eigenvalues
41	REAL lamdamax, cof
42	INTEGER i, j, k, kf
43	REAL dymin, colat0
44	REAL eignft(iim, iim), coff
45
46	! Filtering coefficients (lamda_max * cos(rlat) / lamda):
47	real coefilu(iim, jjm), coefilv(iim, jjm)
48	real coefilu2(iim, jjm), coefilv2(iim, jjm)
49
50	integer modfrstu(jjm), modfrstv(jjm)
51	! index of the mode from where modes are filtered
52
53	!-----------------------------------------------------------
54
55	print *, "Call sequence information: inifilr"
56
57	CALL inifgn(eignvl)
58
59	PRINT *, 'EIGNVL '
60	PRINT "(1X, 5E13.6)", eignvl
61
62	! compute eigenvalues and eigenfunctions
63	! compute the filtering coefficients for scalar lines and
64	! meridional wind v-lines
65	! we filter all those latitude lines where coefil < 1
66	! NO FILTERING AT POLES
67	! colat0 is to be used when alpha (stretching coefficient)
68	! is set equal to zero for the regular grid case
69
70	! Calcul de colat0
71
72	DO j = 1, jjm
73	dlatu(j) = rlatu(j) - rlatu(j + 1)
74	END DO
75
76	dymin = dlatu(1)
77	DO j = 2, jjm
78	dymin = min(dymin, dlatu(j))
79	END DO
80
81	colat0 = min(0.5, dymin / minval(xprimu(:iim)))
82
83	PRINT *, 'colat0 = ', colat0
84
85	lamdamax = iim / (pi * colat0 / grossismx)
86	rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))
87
88	DO j = 1, jjm
89	DO i = 1, iim
90	coefilu(i, j) = 0.
91	coefilv(i, j) = 0.
92	coefilu2(i, j) = 0.
93	coefilv2(i, j) = 0.
94	end DO
95	END DO
96
97	! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
98
99	PRINT *, 'TRUNCATION AT ', iim
100
101	DO j = 2, jjm / 2 + 1
102	IF (cos(rlatu(j)) / colat0 < 1. &
103	.and. rlamda(iim) * cos(rlatu(j)) < 1.) jfiltnu = j
104
105	IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &
106	.and. rlamda(iim) * cos(rlatu(jjm - j + 2)) < 1.) &
107	jfiltsu = jjm - j + 2
108	END DO
109
110	DO j = 1, jjm / 2
111	IF (cos(rlatv(j)) / colat0 < 1. .and. rlamda(iim) * cos(rlatv(j)) < 1.) &
112	jfiltnv = j
113
114	IF (cos(rlatv(jjm - j + 1)) / colat0 < 1. .and. rlamda(iim) &
115	* cos(rlatv(jjm - j + 1)) < 1.) jfiltsv = jjm - j + 1
116	END DO
117
118	IF (jfiltnu <= 0) jfiltnu = 1
119	IF (jfiltnu > jjm / 2 + 1) THEN
120	PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu
121	STOP 1
122	END IF
123
124	IF (jfiltsu <= 0) jfiltsu = 1
125	IF (jfiltsu > jjm + 1) THEN
126	PRINT *, 'jfiltsu en dehors des valeurs acceptables ', jfiltsu
127	STOP 1
128	END IF
129
130	IF (jfiltnv <= 0) jfiltnv = 1
131	IF (jfiltnv > jjm / 2) THEN
132	PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv
133	STOP 1
134	END IF
135
136	IF (jfiltsv <= 0) jfiltsv = 1
137	IF (jfiltsv > jjm) THEN
138	PRINT *, 'jfiltsv en dehors des valeurs acceptables ', jfiltsv
139	STOP 1
140	END IF
141
142	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
143	jfiltsu
144
145	! Determination de coefilu, coefilv, n=modfrstu, modfrstv
146
147	DO j = 1, jjm
148	modfrstu(j) = iim
149	modfrstv(j) = iim
150	END DO
151
152	DO j = 2, jfiltnu
153	DO k = 2, iim
154	IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
155	end DO
156	if (k == iim + 1) cycle
157	modfrstu(j) = k
158
159	kf = modfrstu(j)
160	DO k = kf, iim
161	cof = rlamda(k) * cos(rlatu(j))
162	coefilu(k, j) = cof - 1.
163	coefilu2(k, j) = cof**2 - 1.
164	end DO
165	END DO
166
167	DO j = 1, jfiltnv
168	DO k = 2, iim
169	IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
170	end DO
171	if (k == iim + 1) cycle
172	modfrstv(j) = k
173
174	kf = modfrstv(j)
175	DO k = kf, iim
176	cof = rlamda(k) * cos(rlatv(j))
177	coefilv(k, j) = cof - 1.
178	coefilv2(k, j) = cof**2 - 1.
179	end DO
180	end DO
181
182	DO j = jfiltsu, jjm
183	DO k = 2, iim
184	IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
185	end DO
186	if (k == iim + 1) cycle
187	modfrstu(j) = k
188
189	kf = modfrstu(j)
190	DO k = kf, iim
191	cof = rlamda(k) * cos(rlatu(j))
192	coefilu(k, j) = cof - 1.
193	coefilu2(k, j) = cof**2 - 1.
194	end DO
195	end DO
196
197	DO j = jfiltsv, jjm
198	DO k = 2, iim
199	IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
200	end DO
201	if (k == iim + 1) cycle
202	modfrstv(j) = k
203
204	kf = modfrstv(j)
205	DO k = kf, iim
206	cof = rlamda(k) * cos(rlatv(j))
207	coefilv(k, j) = cof - 1.
208	coefilv2(k, j) = cof**2 - 1.
209	end DO
210	END DO
211
212	IF (jfiltnv>=jjm / 2 .OR. jfiltnu>=jjm / 2) THEN
213	IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
214	IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
215
216	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
217	jfiltsu
218	END IF
219
220	PRINT *, 'Modes premiers v '
221	PRINT 334, modfrstv
222	PRINT *, 'Modes premiers u '
223	PRINT 334, modfrstu
224
225	allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
226	allocate(matricevn(iim, iim, jfiltnv))
227	allocate(matricevs(iim, iim, jfiltsv:jjm))
228	allocate(matriceus(iim, iim, jfiltsu:jjm), matrinvs(iim, iim, jfiltsu:jjm))
229
230	! Calcul de la matrice filtre 'matriceu' pour les champs situes
231	! sur la grille scalaire
232
233	DO j = 2, jfiltnu
234	DO i = 1, iim
235	IF (i < modfrstu(j)) then
236	coff = 0.
237	else
238	coff = coefilu(i, j)
239	end IF
240	eignft(i, :) = eignfnv(:, i) * coff
241	END DO
242	matriceun(:, :, j) = matmul(eignfnv, eignft)
243	END DO
244
245	DO j = jfiltsu, jjm
246	DO i = 1, iim
247	IF (i < modfrstu(j)) then
248	coff = 0.
249	else
250	coff = coefilu(i, j)
251	end IF
252	eignft(i, :) = eignfnv(:, i) * coff
253	END DO
254	matriceus(:, :, j) = matmul(eignfnv, eignft)
255	END DO
256
257	! Calcul de la matrice filtre 'matricev' pour les champs situes
258	! sur la grille de V ou de Z
259
260	DO j = 1, jfiltnv
261	DO i = 1, iim
262	IF (i < modfrstv(j)) then
263	coff = 0.
264	else
265	coff = coefilv(i, j)
266	end IF
267	eignft(i, :) = eignfnu(:, i) * coff
268	END DO
269	matricevn(:, :, j) = matmul(eignfnu, eignft)
270	END DO
271
272	DO j = jfiltsv, jjm
273	DO i = 1, iim
274	IF (i < modfrstv(j)) then
275	coff = 0.
276	else
277	coff = coefilv(i, j)
278	end IF
279	eignft(i, :) = eignfnu(:, i) * coff
280	END DO
281	matricevs(:, :, j) = matmul(eignfnu, eignft)
282	END DO
283
284	! Calcul de la matrice filtre 'matrinv' pour les champs situes
285	! sur la grille scalaire , pour le filtre inverse
286
287	DO j = 2, jfiltnu
288	DO i = 1, iim
289	IF (i < modfrstu(j)) then
290	coff = 0.
291	else
292	coff = coefilu(i, j) / (1. + coefilu(i, j))
293	end IF
294	eignft(i, :) = eignfnv(:, i) * coff
295	END DO
296	matrinvn(:, :, j) = matmul(eignfnv, eignft)
297	END DO
298
299	DO j = jfiltsu, jjm
300	DO i = 1, iim
301	IF (i < modfrstu(j)) then
302	coff = 0.
303	else
304	coff = coefilu(i, j) / (1. + coefilu(i, j))
305	end IF
306	eignft(i, :) = eignfnv(:, i) * coff
307	END DO
308	matrinvs(:, :, j) = matmul(eignfnv, eignft)
309	END DO
310
311	334 FORMAT (1X, 24I3)
312
313	END SUBROUTINE inifilr
314
315	end module inifilr_m