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