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	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	guez	54	USE dimens_m, ONLY : iim, jjm
43	guez	139	USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
44	guez	113	use inifgn_m, only: inifgn
45	guez	54	use nr_util, only: pi
46	guez	3
47	guez	54	! Local:
48	guez	132	REAL dlatu(jjm)
49	guez	54	REAL rlamda(2: iim), eignvl(iim)
50	guez	3
51	guez	54	REAL lamdamax, cof
52			INTEGER i, j, modemax, imx, k, kf
53	guez	132	REAL dymin, colat0
54	guez	54	REAL eignft(iim, iim), coff
55	guez	3
56	guez	54	!-----------------------------------------------------------
57	guez	3
58	guez	54	print *, "Call sequence information: inifilr"
59	guez	3
60	guez	54	CALL inifgn(eignvl)
61	guez	3
62	guez	54	PRINT *, 'EIGNVL '
63			PRINT "(1X, 5E13.6)", eignvl
64	guez	3
65	guez	54	! 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	guez	3
73	guez	54	! Calcul de colat0
74	guez	3
75	guez	54	DO j = 1, jjm
76			dlatu(j) = rlatu(j) - rlatu(j+1)
77			END DO
78	guez	3
79	guez	54	dymin = dlatu(1)
80			DO j = 2, jjm
81			dymin = min(dymin, dlatu(j))
82			END DO
83	guez	3
84	guez	132	colat0 = min(0.5, dymin / minval(xprimu(:iim)))
85	guez	3
86	guez	54	PRINT *, 'colat0 = ', colat0
87	guez	3
88	guez	113	lamdamax = iim / (pi * colat0 / grossismx)
89	guez	54	rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))
90	guez	3
91	guez	54	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	guez	3
100	guez	54	! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
101	guez	3
102	guez	54	modemax = iim
103			imx = iim
104	guez	3
105	guez	54	PRINT *, 'TRUNCATION AT ', imx
106	guez	3
107	guez	54	DO j = 2, jjm / 2 + 1
108			IF (cos(rlatu(j)) / colat0 < 1. &
109			.and. rlamda(imx) * cos(rlatu(j)) < 1.) jfiltnu = j
110	guez	3
111	guez	54	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	guez	3
116	guez	54	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	guez	3
122	guez	54	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	guez	3
128	guez	54	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	guez	3
134	guez	54	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	guez	3
140	guez	54	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	guez	3
146	guez	54	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	guez	3
152	guez	54	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
153			jfiltsu
154	guez	32
155	guez	54	! Determination de coefilu, coefilv, n=modfrstu, modfrstv
156	guez	32
157	guez	54	DO j = 1, jjm
158			modfrstu(j) = iim
159			modfrstv(j) = iim
160			END DO
161	guez	32
162	guez	54	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	guez	32
170	guez	54	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	guez	32
178	guez	54	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	guez	32
186	guez	54	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	guez	32
194	guez	54	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	guez	32
202	guez	54	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	guez	32
210	guez	54	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	guez	32
218	guez	54	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	guez	32
226	guez	54	IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN
227			IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
228			IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
229	guez	32
230	guez	54	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
231			jfiltsu
232			END IF
233	guez	32
234	guez	54	PRINT *, 'Modes premiers v '
235			PRINT 334, modfrstv
236			PRINT *, 'Modes premiers u '
237			PRINT 334, modfrstu
238	guez	32
239	guez	136	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	guez	32
244	guez	54	! Calcul de la matrice filtre 'matriceu' pour les champs situes
245			! sur la grille scalaire
246	guez	32
247	guez	54	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	guez	32
259	guez	54	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	guez	136	matriceus(:, :, j) = matmul(eignfnv, eignft)
269	guez	54	END DO
270	guez	32
271	guez	54	! Calcul de la matrice filtre 'matricev' pour les champs situes
272			! sur la grille de V ou de Z
273	guez	32
274	guez	54	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	guez	32
286	guez	54	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	guez	136	matricevs(:, :, j) = matmul(eignfnu, eignft)
296	guez	54	END DO
297	guez	32
298	guez	54	! Calcul de la matrice filtre 'matrinv' pour les champs situes
299			! sur la grille scalaire , pour le filtre inverse
300	guez	32
301	guez	54	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	guez	32
313	guez	54	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	guez	136	matrinvs(:, :, j) = matmul(eignfnv, eignft)
323	guez	54	END DO
324	guez	32
325	guez	54	334 FORMAT (1X, 24I3)
326	guez	32
327	guez	54	END SUBROUTINE inifilr
328	guez	32
329	guez	54	end module inifilr_m