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	module inifilr_m
2
3	IMPLICIT NONE
4
5	INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv
6	! jfiltn index of the last scalar line filtered in NH
7	! jfilts index of the first line filtered in SH
8
9	! North:
10	real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :)
11	! (iim, iim, 2:jfiltnu)
12
13	real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)
14
15	! South:
16	real, allocatable:: matriceus(:, :, :), matrinvs(:, :, :)
17	! (iim, iim, jfiltsu:jjm)
18
19	real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
20
21	contains
22
23	SUBROUTINE inifilr
24
25	! From filtrez/inifilr.F, version 1.1.1.1 2004/05/19 12:53:09
26	! H. Upadhyaya, O. Sharma
27
28	! This routine computes the eigenfunctions of the laplacian on the
29	! stretched grid, and the filtering coefficients. The modes are
30	! filtered from modfrst to iim.
31
32	USE dimens_m, ONLY : iim, jjm
33	USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
34	use inifgn_m, only: inifgn, eignfnu, eignfnv
35	use jumble, only: new_unit
36	use nr_util, only: pi
37
38	! Local:
39	REAL dlatu(jjm)
40	REAL rlamda(2: iim)
41	real eignvl(iim) ! eigenvalues
42	REAL cof
43	INTEGER i, j, k, kf, unit
44	REAL colat0
45	REAL eignft(iim, iim), coff
46
47	! 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	!-----------------------------------------------------------
55
56	print *, "Call sequence information: inifilr"
57
58	CALL inifgn(eignvl)
59
60	call new_unit(unit)
61	open(unit, file = "eignvl.txt", status = "replace", action = "write")
62	write(unit, fmt = *) EIGNVL
63	close(unit)
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	forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
75	colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
76	PRINT *, 'colat0 = ', colat0
77
78	rlamda = iim / (pi * colat0 / grossismx) / sqrt(abs(eignvl(2: iim)))
79	coefilu = 0.
80	coefilv = 0.
81	coefilu2 = 0.
82	coefilv2 = 0.
83
84	! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
85
86	DO j = 2, jjm / 2 + 1
87	IF (cos(rlatu(j)) / colat0 < 1. &
88	.and. rlamda(iim) * cos(rlatu(j)) < 1.) jfiltnu = j
89
90	IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &
91	.and. rlamda(iim) * cos(rlatu(jjm - j + 2)) < 1.) &
92	jfiltsu = jjm - j + 2
93	END DO
94
95	DO j = 1, jjm / 2
96	IF (cos(rlatv(j)) / colat0 < 1. .and. rlamda(iim) * cos(rlatv(j)) < 1.) &
97	jfiltnv = j
98
99	IF (cos(rlatv(jjm - j + 1)) / colat0 < 1. .and. rlamda(iim) &
100	* cos(rlatv(jjm - j + 1)) < 1.) jfiltsv = jjm - j + 1
101	END DO
102
103	IF (jfiltnu <= 0) jfiltnu = 1
104	IF (jfiltnu > jjm / 2 + 1) THEN
105	PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu
106	STOP 1
107	END IF
108
109	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
115	IF (jfiltnv <= 0) jfiltnv = 1
116	IF (jfiltnv > jjm / 2) THEN
117	PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv
118	STOP 1
119	END IF
120
121	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
127	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
128	jfiltsu
129
130	! Determination de coefilu, coefilv, n=modfrstu, modfrstv
131
132	DO j = 1, jjm
133	modfrstu(j) = iim
134	modfrstv(j) = iim
135	END DO
136
137	DO j = 2, jfiltnu
138	DO k = 2, iim
139	IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
140	end DO
141	if (k == iim + 1) cycle
142	modfrstu(j) = k
143
144	kf = modfrstu(j)
145	DO k = kf, iim
146	cof = rlamda(k) * cos(rlatu(j))
147	coefilu(k, j) = cof - 1.
148	coefilu2(k, j) = cof**2 - 1.
149	end DO
150	END DO
151
152	DO j = 1, jfiltnv
153	DO k = 2, iim
154	IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
155	end DO
156	if (k == iim + 1) cycle
157	modfrstv(j) = k
158
159	kf = modfrstv(j)
160	DO k = kf, iim
161	cof = rlamda(k) * cos(rlatv(j))
162	coefilv(k, j) = cof - 1.
163	coefilv2(k, j) = cof**2 - 1.
164	end DO
165	end DO
166
167	DO j = jfiltsu, jjm
168	DO k = 2, iim
169	IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
170	end DO
171	if (k == iim + 1) cycle
172	modfrstu(j) = k
173
174	kf = modfrstu(j)
175	DO k = kf, iim
176	cof = rlamda(k) * cos(rlatu(j))
177	coefilu(k, j) = cof - 1.
178	coefilu2(k, j) = cof**2 - 1.
179	end DO
180	end DO
181
182	DO j = jfiltsv, jjm
183	DO k = 2, iim
184	IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
185	end DO
186	if (k == iim + 1) cycle
187	modfrstv(j) = k
188
189	kf = modfrstv(j)
190	DO k = kf, iim
191	cof = rlamda(k) * cos(rlatv(j))
192	coefilv(k, j) = cof - 1.
193	coefilv2(k, j) = cof**2 - 1.
194	end DO
195	END DO
196
197	IF (jfiltnv>=jjm / 2 .OR. jfiltnu>=jjm / 2) THEN
198	IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
199	IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
200
201	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
202	jfiltsu
203	END IF
204
205	PRINT *, 'Modes premiers v '
206	PRINT 334, modfrstv
207	PRINT *, 'Modes premiers u '
208	PRINT 334, modfrstu
209
210	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
215	! Calcul de la matrice filtre 'matriceu' pour les champs situes
216	! sur la grille scalaire
217
218	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	eignft(i, :) = eignfnv(:, i) * coff
226	END DO
227	matriceun(:, :, j) = matmul(eignfnv, eignft)
228	END DO
229
230	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	matriceus(:, :, j) = matmul(eignfnv, eignft)
240	END DO
241
242	! Calcul de la matrice filtre 'matricev' pour les champs situes
243	! sur la grille de V ou de Z
244
245	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	eignft(i, :) = eignfnu(:, i) * coff
253	END DO
254	matricevn(:, :, j) = matmul(eignfnu, eignft)
255	END DO
256
257	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	eignft(i, :) = eignfnu(:, i) * coff
265	END DO
266	matricevs(:, :, j) = matmul(eignfnu, eignft)
267	END DO
268
269	! Calcul de la matrice filtre 'matrinv' pour les champs situes
270	! sur la grille scalaire , pour le filtre inverse
271
272	DO j = 2, jfiltnu
273	DO i = 1, iim
274	IF (i < modfrstu(j)) then
275	coff = 0.
276	else
277	coff = coefilu(i, j) / (1. + coefilu(i, j))
278	end IF
279	eignft(i, :) = eignfnv(:, i) * coff
280	END DO
281	matrinvn(:, :, j) = matmul(eignfnv, eignft)
282	END DO
283
284	DO j = jfiltsu, jjm
285	DO i = 1, iim
286	IF (i < modfrstu(j)) then
287	coff = 0.
288	else
289	coff = coefilu(i, j) / (1. + coefilu(i, j))
290	end IF
291	eignft(i, :) = eignfnv(:, i) * coff
292	END DO
293	matrinvs(:, :, j) = matmul(eignfnv, eignft)
294	END DO
295
296	334 FORMAT (1X, 24I3)
297
298	END SUBROUTINE inifilr
299
300	end module inifilr_m