Sources/filtrez/inifilr.f

module inifilr_m

  IMPLICIT NONE

  INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv

  ! 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.
    ! 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 dlatu(jjm)
    REAL rlamda(2: iim), eignvl(iim)

    REAL lamdamax, cof
    INTEGER i, j, modemax, imx, k, kf
    REAL dymin, colat0
    REAL eignft(iim, iim), coff

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

    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

    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

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