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	guez	54	module inifilr_m
2	guez	3
3	guez	32	IMPLICIT NONE
4	guez	3
5	guez	54	INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv
6	guez	3
7	guez	136	! North:
8			real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :)
9			! (iim, iim, 2:jfiltnu)
10	guez	3
11	guez	136	real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)
12	guez	3
13	guez	136	! South:
14			real, allocatable:: matriceus(:, :, :), matrinvs(:, :, :)
15			! (iim, iim, jfiltsu:jjm)
16
17			real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
18
19	guez	54	contains
20	guez	3
21	guez	54	SUBROUTINE inifilr
22	guez	3
23	guez	54	! From filtrez/inifilr.F, version 1.1.1.1 2004/05/19 12:53:09
24			! H. Upadhyaya, O. Sharma
25	guez	3
26	guez	54	! 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	guez	3
38	guez	54	! The modes are filtered from modfrst to modemax.
39	guez	3
40	guez	113	USE coefils, ONLY : coefilu, coefilu2, coefilv, coefilv2, eignfnu, &
41			eignfnv, modfrstu, modfrstv
42			USE comgeom, ONLY : rlatu, rlatv, xprimu
43	guez	54	USE dimens_m, ONLY : iim, jjm
44	guez	113	use inifgn_m, only: inifgn
45	guez	54	use nr_util, only: pi
46	guez	113	USE serre, ONLY : grossismx
47	guez	3
48	guez	54	! Local:
49	guez	132	REAL dlatu(jjm)
50	guez	54	REAL rlamda(2: iim), eignvl(iim)
51	guez	3
52	guez	54	REAL lamdamax, cof
53			INTEGER i, j, modemax, imx, k, kf
54	guez	132	REAL dymin, colat0
55	guez	54	REAL eignft(iim, iim), coff
56	guez	3
57	guez	54	!-----------------------------------------------------------
58	guez	3
59	guez	54	print *, "Call sequence information: inifilr"
60	guez	3
61	guez	54	CALL inifgn(eignvl)
62	guez	3
63	guez	54	PRINT *, 'EIGNVL '
64			PRINT "(1X, 5E13.6)", eignvl
65	guez	3
66	guez	54	! 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	guez	3
74	guez	54	! Calcul de colat0
75	guez	3
76	guez	54	DO j = 1, jjm
77			dlatu(j) = rlatu(j) - rlatu(j+1)
78			END DO
79	guez	3
80	guez	54	dymin = dlatu(1)
81			DO j = 2, jjm
82			dymin = min(dymin, dlatu(j))
83			END DO
84	guez	3
85	guez	132	colat0 = min(0.5, dymin / minval(xprimu(:iim)))
86	guez	3
87	guez	54	PRINT *, 'colat0 = ', colat0
88	guez	3
89	guez	113	lamdamax = iim / (pi * colat0 / grossismx)
90	guez	54	rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))
91	guez	3
92	guez	54	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	guez	3
101	guez	54	! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
102	guez	3
103	guez	54	modemax = iim
104			imx = iim
105	guez	3
106	guez	54	PRINT *, 'TRUNCATION AT ', imx
107	guez	3
108	guez	54	DO j = 2, jjm / 2 + 1
109			IF (cos(rlatu(j)) / colat0 < 1. &
110			.and. rlamda(imx) * cos(rlatu(j)) < 1.) jfiltnu = j
111	guez	3
112	guez	54	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	guez	3
117	guez	54	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	guez	3
123	guez	54	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	guez	3
129	guez	54	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	guez	3
135	guez	54	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	guez	3
141	guez	54	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	guez	3
147	guez	54	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	guez	3
153	guez	54	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
154			jfiltsu
155	guez	32
156	guez	54	! Determination de coefilu, coefilv, n=modfrstu, modfrstv
157	guez	32
158	guez	54	DO j = 1, jjm
159			modfrstu(j) = iim
160			modfrstv(j) = iim
161			END DO
162	guez	32
163	guez	54	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	guez	32
171	guez	54	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	guez	32
179	guez	54	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	guez	32
187	guez	54	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	guez	32
195	guez	54	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	guez	32
203	guez	54	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	guez	32
211	guez	54	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	guez	32
219	guez	54	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	guez	32
227	guez	54	IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN
228			IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
229			IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
230	guez	32
231	guez	54	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
232			jfiltsu
233			END IF
234	guez	32
235	guez	54	PRINT *, 'Modes premiers v '
236			PRINT 334, modfrstv
237			PRINT *, 'Modes premiers u '
238			PRINT 334, modfrstu
239	guez	32
240	guez	136	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	guez	32
245	guez	54	! Calcul de la matrice filtre 'matriceu' pour les champs situes
246			! sur la grille scalaire
247	guez	32
248	guez	54	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	guez	32
260	guez	54	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	guez	136	matriceus(:, :, j) = matmul(eignfnv, eignft)
270	guez	54	END DO
271	guez	32
272	guez	54	! Calcul de la matrice filtre 'matricev' pour les champs situes
273			! sur la grille de V ou de Z
274	guez	32
275	guez	54	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	guez	32
287	guez	54	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	guez	136	matricevs(:, :, j) = matmul(eignfnu, eignft)
297	guez	54	END DO
298	guez	32
299	guez	54	! Calcul de la matrice filtre 'matrinv' pour les champs situes
300			! sur la grille scalaire , pour le filtre inverse
301	guez	32
302	guez	54	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	guez	32
314	guez	54	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	guez	136	matrinvs(:, :, j) = matmul(eignfnv, eignft)
324	guez	54	END DO
325	guez	32
326	guez	54	334 FORMAT (1X, 24I3)
327	guez	32
328	guez	54	END SUBROUTINE inifilr
329	guez	32
330	guez	54	end module inifilr_m