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 dimens_m, ONLY : iim, jjm
    USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
    use inifgn_m, only: inifgn
    use nr_util, only: pi

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