Sources/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 jumble, only: new_unit
    use nr_util, only: pi

    ! Local:
    REAL dlatu(jjm)
    REAL rlamda(2: iim)
    real eignvl(iim) ! eigenvalues
    REAL cof
    INTEGER i, j, k, kf, unit
    REAL 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)

    call new_unit(unit)
    open(unit, file = "eignvl.txt", status = "replace", action = "write") 
    write(unit, fmt = *) EIGNVL
    close(unit)

    ! 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
    forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
    colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
    PRINT *, 'colat0 = ', colat0

    rlamda = iim / (pi * colat0 / grossismx) / sqrt(abs(eignvl(2: iim)))
    coefilu = 0.
    coefilv = 0.
    coefilu2 = 0.
    coefilv2 = 0.

    ! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv

    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	143	use jumble, only: new_unit
36	guez	54	use nr_util, only: pi
37	guez	3
38	guez	54	! Local:
39	guez	132	REAL dlatu(jjm)
40	guez	140	REAL rlamda(2: iim)
41			real eignvl(iim) ! eigenvalues
42	guez	143	REAL cof
43			INTEGER i, j, k, kf, unit
44			REAL colat0
45	guez	54	REAL eignft(iim, iim), coff
46	guez	3
47	guez	140	! Filtering coefficients (lamda_max * cos(rlat) / lamda):
48			real coefilu(iim, jjm), coefilv(iim, jjm)
49			real coefilu2(iim, jjm), coefilv2(iim, jjm)
50
51			integer modfrstu(jjm), modfrstv(jjm)
52			! index of the mode from where modes are filtered
53
54	guez	54	!-----------------------------------------------------------
55	guez	3
56	guez	54	print *, "Call sequence information: inifilr"
57	guez	3
58	guez	54	CALL inifgn(eignvl)
59	guez	3
60	guez	143	call new_unit(unit)
61			open(unit, file = "eignvl.txt", status = "replace", action = "write")
62			write(unit, fmt = *) EIGNVL
63			close(unit)
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	143	forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
75			colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
76	guez	54	PRINT *, 'colat0 = ', colat0
77	guez	3
78	guez	143	rlamda = iim / (pi * colat0 / grossismx) / sqrt(abs(eignvl(2: iim)))
79			coefilu = 0.
80			coefilv = 0.
81			coefilu2 = 0.
82			coefilv2 = 0.
83	guez	3
84	guez	54	! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
85	guez	3
86	guez	54	DO j = 2, jjm / 2 + 1
87			IF (cos(rlatu(j)) / colat0 < 1. &
88	guez	140	.and. rlamda(iim) * cos(rlatu(j)) < 1.) jfiltnu = j
89	guez	3
90	guez	54	IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &
91	guez	140	.and. rlamda(iim) * cos(rlatu(jjm - j + 2)) < 1.) &
92	guez	54	jfiltsu = jjm - j + 2
93			END DO
94	guez	3
95	guez	140	DO j = 1, jjm / 2
96			IF (cos(rlatv(j)) / colat0 < 1. .and. rlamda(iim) * cos(rlatv(j)) < 1.) &
97			jfiltnv = j
98	guez	3
99	guez	140	IF (cos(rlatv(jjm - j + 1)) / colat0 < 1. .and. rlamda(iim) &
100			* cos(rlatv(jjm - j + 1)) < 1.) jfiltsv = jjm - j + 1
101	guez	54	END DO
102	guez	3
103	guez	54	IF (jfiltnu <= 0) jfiltnu = 1
104	guez	140	IF (jfiltnu > jjm / 2 + 1) THEN
105	guez	54	PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu
106			STOP 1
107			END IF
108	guez	3
109	guez	54	IF (jfiltsu <= 0) jfiltsu = 1
110			IF (jfiltsu > jjm + 1) THEN
111			PRINT *, 'jfiltsu en dehors des valeurs acceptables ', jfiltsu
112			STOP 1
113			END IF
114	guez	3
115	guez	54	IF (jfiltnv <= 0) jfiltnv = 1
116	guez	140	IF (jfiltnv > jjm / 2) THEN
117	guez	54	PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv
118			STOP 1
119			END IF
120	guez	3
121	guez	54	IF (jfiltsv <= 0) jfiltsv = 1
122			IF (jfiltsv > jjm) THEN
123			PRINT *, 'jfiltsv en dehors des valeurs acceptables ', jfiltsv
124			STOP 1
125			END IF
126	guez	3
127	guez	54	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
128			jfiltsu
129	guez	32
130	guez	54	! Determination de coefilu, coefilv, n=modfrstu, modfrstv
131	guez	32
132	guez	54	DO j = 1, jjm
133			modfrstu(j) = iim
134			modfrstv(j) = iim
135			END DO
136	guez	32
137	guez	54	DO j = 2, jfiltnu
138	guez	140	DO k = 2, iim
139			IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
140	guez	54	end DO
141	guez	140	if (k == iim + 1) cycle
142	guez	54	modfrstu(j) = k
143	guez	32
144	guez	54	kf = modfrstu(j)
145	guez	140	DO k = kf, iim
146			cof = rlamda(k) * cos(rlatu(j))
147	guez	54	coefilu(k, j) = cof - 1.
148	guez	140	coefilu2(k, j) = cof**2 - 1.
149	guez	54	end DO
150			END DO
151	guez	32
152	guez	54	DO j = 1, jfiltnv
153	guez	140	DO k = 2, iim
154			IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
155	guez	54	end DO
156	guez	140	if (k == iim + 1) cycle
157	guez	54	modfrstv(j) = k
158	guez	32
159	guez	54	kf = modfrstv(j)
160	guez	140	DO k = kf, iim
161			cof = rlamda(k) * cos(rlatv(j))
162	guez	54	coefilv(k, j) = cof - 1.
163	guez	140	coefilv2(k, j) = cof**2 - 1.
164	guez	54	end DO
165			end DO
166	guez	32
167	guez	54	DO j = jfiltsu, jjm
168	guez	140	DO k = 2, iim
169			IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
170	guez	54	end DO
171	guez	140	if (k == iim + 1) cycle
172	guez	54	modfrstu(j) = k
173	guez	32
174	guez	54	kf = modfrstu(j)
175	guez	140	DO k = kf, iim
176			cof = rlamda(k) * cos(rlatu(j))
177	guez	54	coefilu(k, j) = cof - 1.
178	guez	140	coefilu2(k, j) = cof**2 - 1.
179	guez	54	end DO
180			end DO
181	guez	32
182	guez	54	DO j = jfiltsv, jjm
183	guez	140	DO k = 2, iim
184			IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
185	guez	54	end DO
186	guez	140	if (k == iim + 1) cycle
187	guez	54	modfrstv(j) = k
188	guez	32
189	guez	54	kf = modfrstv(j)
190	guez	140	DO k = kf, iim
191			cof = rlamda(k) * cos(rlatv(j))
192	guez	54	coefilv(k, j) = cof - 1.
193	guez	140	coefilv2(k, j) = cof**2 - 1.
194	guez	54	end DO
195			END DO
196	guez	32
197	guez	140	IF (jfiltnv>=jjm / 2 .OR. jfiltnu>=jjm / 2) THEN
198	guez	54	IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
199			IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
200	guez	32
201	guez	54	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
202			jfiltsu
203			END IF
204	guez	32
205	guez	54	PRINT *, 'Modes premiers v '
206			PRINT 334, modfrstv
207			PRINT *, 'Modes premiers u '
208			PRINT 334, modfrstu
209	guez	32
210	guez	136	allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
211			allocate(matricevn(iim, iim, jfiltnv))
212			allocate(matricevs(iim, iim, jfiltsv:jjm))
213			allocate(matriceus(iim, iim, jfiltsu:jjm), matrinvs(iim, iim, jfiltsu:jjm))
214	guez	32
215	guez	54	! Calcul de la matrice filtre 'matriceu' pour les champs situes
216			! sur la grille scalaire
217	guez	32
218	guez	54	DO j = 2, jfiltnu
219			DO i = 1, iim
220			IF (i < modfrstu(j)) then
221			coff = 0.
222			else
223			coff = coefilu(i, j)
224			end IF
225	guez	140	eignft(i, :) = eignfnv(:, i) * coff
226	guez	54	END DO
227			matriceun(:, :, j) = matmul(eignfnv, eignft)
228			END DO
229	guez	32
230	guez	54	DO j = jfiltsu, jjm
231			DO i = 1, iim
232			IF (i < modfrstu(j)) then
233			coff = 0.
234			else
235			coff = coefilu(i, j)
236			end IF
237			eignft(i, :) = eignfnv(:, i) * coff
238			END DO
239	guez	136	matriceus(:, :, j) = matmul(eignfnv, eignft)
240	guez	54	END DO
241	guez	32
242	guez	54	! Calcul de la matrice filtre 'matricev' pour les champs situes
243			! sur la grille de V ou de Z
244	guez	32
245	guez	54	DO j = 1, jfiltnv
246			DO i = 1, iim
247			IF (i < modfrstv(j)) then
248			coff = 0.
249			else
250			coff = coefilv(i, j)
251			end IF
252	guez	140	eignft(i, :) = eignfnu(:, i) * coff
253	guez	54	END DO
254			matricevn(:, :, j) = matmul(eignfnu, eignft)
255			END DO
256	guez	32
257	guez	54	DO j = jfiltsv, jjm
258			DO i = 1, iim
259			IF (i < modfrstv(j)) then
260			coff = 0.
261			else
262			coff = coefilv(i, j)
263			end IF
264	guez	140	eignft(i, :) = eignfnu(:, i) * coff
265	guez	54	END DO
266	guez	136	matricevs(:, :, j) = matmul(eignfnu, eignft)
267	guez	54	END DO
268	guez	32
269	guez	54	! Calcul de la matrice filtre 'matrinv' pour les champs situes
270			! sur la grille scalaire , pour le filtre inverse
271	guez	32
272	guez	54	DO j = 2, jfiltnu
273			DO i = 1, iim
274			IF (i < modfrstu(j)) then
275			coff = 0.
276			else
277	guez	140	coff = coefilu(i, j) / (1. + coefilu(i, j))
278	guez	54	end IF
279	guez	140	eignft(i, :) = eignfnv(:, i) * coff
280	guez	54	END DO
281			matrinvn(:, :, j) = matmul(eignfnv, eignft)
282			END DO
283	guez	32
284	guez	54	DO j = jfiltsu, jjm
285			DO i = 1, iim
286			IF (i < modfrstu(j)) then
287			coff = 0.
288			else
289	guez	140	coff = coefilu(i, j) / (1. + coefilu(i, j))
290	guez	54	end IF
291	guez	140	eignft(i, :) = eignfnv(:, i) * coff
292	guez	54	END DO
293	guez	136	matrinvs(:, :, j) = matmul(eignfnv, eignft)
294	guez	54	END DO
295	guez	32
296	guez	54	334 FORMAT (1X, 24I3)
297	guez	32
298	guez	54	END SUBROUTINE inifilr
299	guez	32
300	guez	54	end module inifilr_m