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

    ! Local:
    REAL dlatu(jjm)
    REAL rlamda(2: iim)
    real eignvl(iim) ! eigenvalues sorted in descending order
    REAL cof
    INTEGER i, j, k, unit
    REAL colat0 ! > 0
    REAL eignft(iim, iim), coff

    real eignfnu(iim, iim), eignfnv(iim, iim)
    ! eigenvectors of the discrete laplacian

    ! Filtering coefficients (lamda_max * cos(rlat) / lamda):
    real coefilu(iim, jjm), coefilv(iim, jjm)
    real coefilu2(iim, jjm), coefilv2(iim, jjm)

    ! Index of the mode from where modes are filtered:
    integer, allocatable:: modfrstnu(:), modfrstsu(:)
    integer, allocatable:: modfrstnv(:), modfrstsv(:)

    !-----------------------------------------------------------

    print *, "Call sequence information: inifilr"

    CALL inifgn(eignvl, eignfnu, eignfnv)

    ! compute eigenvalues and eigenvectors
    ! 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)))

    ! Determination de jfiltnu, jfiltsu, jfiltnv, jfiltsv

    jfiltnu = (jjm + 1) / 2
    do while (cos(rlatu(jfiltnu)) >= colat0 &
         .or. rlamda(iim) * cos(rlatu(jfiltnu)) >= 1.)
       jfiltnu = jfiltnu - 1
    end do

    jfiltsu = jjm / 2 + 2
    do while (cos(rlatu(jfiltsu)) >= colat0 &
         .or. rlamda(iim) * cos(rlatu(jfiltsu)) >= 1.)
       jfiltsu = jfiltsu + 1
    end do

    jfiltnv = jjm / 2
    do while ((cos(rlatv(jfiltnv)) >= colat0 &
         .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) .and. jfiltnv >= 2)
       jfiltnv = jfiltnv - 1
    end do

    if (cos(rlatv(jfiltnv)) >= colat0 &
         .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) then
       ! {jfiltnv == 1}
       PRINT *, 'Could not find jfiltnv.'
       STOP 1
    END IF

    jfiltsv = (jjm + 1)/ 2 + 1
    do while ((cos(rlatv(jfiltsv)) >= colat0 &
         .or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) .and. jfiltsv <= jjm - 1)
       jfiltsv = jfiltsv + 1
    end do

    IF (cos(rlatv(jfiltsv)) >= colat0 &
         .or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) THEN
       ! {jfiltsv == jjm}
       PRINT *, 'Could not find jfiltsv.'
       STOP 1
    END IF

    PRINT *, 'jfiltnu =', jfiltnu
    PRINT *, 'jfiltsu =', jfiltsu
    PRINT *, 'jfiltnv =', jfiltnv
    PRINT *, 'jfiltsv =', jfiltsv

    ! Determination de coefilu, coefilv, modfrst[ns][uv]:

    allocate(modfrstnu(2:jfiltnu), modfrstsu(jfiltsu:jjm))
    allocate(modfrstnv(jfiltnv), modfrstsv(jfiltsv:jjm))
    coefilu = 0.
    coefilv = 0.
    coefilu2 = 0.
    coefilv2 = 0.

    DO j = 2, jfiltnu
       modfrstnu(j) = 2
       do while (rlamda(modfrstnu(j)) * cos(rlatu(j)) >= 1. &
            .and. modfrstnu(j) <= iim - 1)
          modfrstnu(j) = modfrstnu(j) + 1
       end do

       if (rlamda(modfrstnu(j)) * cos(rlatu(j)) < 1.) then
          DO k = modfrstnu(j), iim
             cof = rlamda(k) * cos(rlatu(j))
             coefilu(k, j) = cof - 1.
             coefilu2(k, j) = cof**2 - 1.
          end DO
       end if
    END DO

    DO j = 1, jfiltnv
       modfrstnv(j) = 2
       do while (rlamda(modfrstnv(j)) * cos(rlatv(j)) >= 1. &
            .and. modfrstnv(j) <= iim - 1)
          modfrstnv(j) = modfrstnv(j) + 1
       end do

       if (rlamda(modfrstnv(j)) * cos(rlatv(j)) < 1.) then
          DO k = modfrstnv(j), iim
             cof = rlamda(k) * cos(rlatv(j))
             coefilv(k, j) = cof - 1.
             coefilv2(k, j) = cof**2 - 1.
          end DO
       end if
    end DO

    DO j = jfiltsu, jjm
       modfrstsu(j) = 2
       do while (rlamda(modfrstsu(j)) * cos(rlatu(j)) >= 1. &
            .and. modfrstsu(j) <= iim - 1)
          modfrstsu(j) = modfrstsu(j) + 1
       end do

       if (rlamda(modfrstsu(j)) * cos(rlatu(j)) < 1.) then
          DO k = modfrstsu(j), iim
             cof = rlamda(k) * cos(rlatu(j))
             coefilu(k, j) = cof - 1.
             coefilu2(k, j) = cof**2 - 1.
          end DO
       end if
    end DO

    DO j = jfiltsv, jjm
       modfrstsv(j) = 2
       do while (rlamda(modfrstsv(j)) * cos(rlatv(j)) >= 1. &
            .and. modfrstsv(j) <= iim - 1)
          modfrstsv(j) = modfrstsv(j) + 1
       end do

       if (rlamda(modfrstsv(j)) * cos(rlatv(j)) < 1.) then
          DO k = modfrstsv(j), iim
             cof = rlamda(k) * cos(rlatv(j))
             coefilv(k, j) = cof - 1.
             coefilv2(k, j) = cof**2 - 1.
          end DO
       end if
    END DO

    call new_unit(unit)
    open(unit, file = "inifilr_out.txt", status = "replace", action = "write") 
    write(unit, fmt = *) '"EIGNVL"', eignvl
    write(unit, fmt = *) '"modfrstnu"', modfrstnu
    write(unit, fmt = *) '"modfrstsu"', modfrstsu
    write(unit, fmt = *) '"modfrstnv"', modfrstnv
    write(unit, fmt = *) '"modfrstsv"', modfrstsv
    close(unit)

    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 < modfrstnu(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 < modfrstsu(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 < modfrstnv(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 < modfrstsv(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 < modfrstnu(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 < modfrstsu(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

  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
11	real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :)
12	! (iim, iim, 2:jfiltnu)
13
14	real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)
15
16	! South:
17
18	real, allocatable:: matriceus(:, :, :), matrinvs(:, :, :)
19	! (iim, iim, jfiltsu:jjm)
20
21	real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
22
23	contains
24
25	SUBROUTINE inifilr
26
27	! From filtrez/inifilr.F, version 1.1.1.1 2004/05/19 12:53:09
28	! H. Upadhyaya, O. Sharma
29
30	! This routine computes the eigenvectors of the laplacian on the
31	! stretched grid, and the filtering coefficients. The modes are
32	! filtered from modfrst to iim.
33
34	USE dimens_m, ONLY : iim, jjm
35	USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
36	use inifgn_m, only: inifgn
37	use jumble, only: new_unit
38	use nr_util, only: pi
39
40	! Local:
41	REAL dlatu(jjm)
42	REAL rlamda(2: iim)
43	real eignvl(iim) ! eigenvalues sorted in descending order
44	REAL cof
45	INTEGER i, j, k, unit
46	REAL colat0 ! > 0
47	REAL eignft(iim, iim), coff
48
49	real eignfnu(iim, iim), eignfnv(iim, iim)
50	! eigenvectors of the discrete laplacian
51
52	! Filtering coefficients (lamda_max * cos(rlat) / lamda):
53	real coefilu(iim, jjm), coefilv(iim, jjm)
54	real coefilu2(iim, jjm), coefilv2(iim, jjm)
55
56	! Index of the mode from where modes are filtered:
57	integer, allocatable:: modfrstnu(:), modfrstsu(:)
58	integer, allocatable:: modfrstnv(:), modfrstsv(:)
59
60	!-----------------------------------------------------------
61
62	print *, "Call sequence information: inifilr"
63
64	CALL inifgn(eignvl, eignfnu, eignfnv)
65
66	! compute eigenvalues and eigenvectors
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
74	! Calcul de colat0
75	forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
76	colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
77	PRINT *, 'colat0 = ', colat0
78
79	rlamda = iim / (pi * colat0 / grossismx) / sqrt(abs(eignvl(2: iim)))
80
81	! Determination de jfiltnu, jfiltsu, jfiltnv, jfiltsv
82
83	jfiltnu = (jjm + 1) / 2
84	do while (cos(rlatu(jfiltnu)) >= colat0 &
85	.or. rlamda(iim) * cos(rlatu(jfiltnu)) >= 1.)
86	jfiltnu = jfiltnu - 1
87	end do
88
89	jfiltsu = jjm / 2 + 2
90	do while (cos(rlatu(jfiltsu)) >= colat0 &
91	.or. rlamda(iim) * cos(rlatu(jfiltsu)) >= 1.)
92	jfiltsu = jfiltsu + 1
93	end do
94
95	jfiltnv = jjm / 2
96	do while ((cos(rlatv(jfiltnv)) >= colat0 &
97	.or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) .and. jfiltnv >= 2)
98	jfiltnv = jfiltnv - 1
99	end do
100
101	if (cos(rlatv(jfiltnv)) >= colat0 &
102	.or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) then
103	! {jfiltnv == 1}
104	PRINT *, 'Could not find jfiltnv.'
105	STOP 1
106	END IF
107
108	jfiltsv = (jjm + 1)/ 2 + 1
109	do while ((cos(rlatv(jfiltsv)) >= colat0 &
110	.or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) .and. jfiltsv <= jjm - 1)
111	jfiltsv = jfiltsv + 1
112	end do
113
114	IF (cos(rlatv(jfiltsv)) >= colat0 &
115	.or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) THEN
116	! {jfiltsv == jjm}
117	PRINT *, 'Could not find jfiltsv.'
118	STOP 1
119	END IF
120
121	PRINT *, 'jfiltnu =', jfiltnu
122	PRINT *, 'jfiltsu =', jfiltsu
123	PRINT *, 'jfiltnv =', jfiltnv
124	PRINT *, 'jfiltsv =', jfiltsv
125
126	! Determination de coefilu, coefilv, modfrst[ns][uv]:
127
128	allocate(modfrstnu(2:jfiltnu), modfrstsu(jfiltsu:jjm))
129	allocate(modfrstnv(jfiltnv), modfrstsv(jfiltsv:jjm))
130	coefilu = 0.
131	coefilv = 0.
132	coefilu2 = 0.
133	coefilv2 = 0.
134
135	DO j = 2, jfiltnu
136	modfrstnu(j) = 2
137	do while (rlamda(modfrstnu(j)) * cos(rlatu(j)) >= 1. &
138	.and. modfrstnu(j) <= iim - 1)
139	modfrstnu(j) = modfrstnu(j) + 1
140	end do
141
142	if (rlamda(modfrstnu(j)) * cos(rlatu(j)) < 1.) then
143	DO k = modfrstnu(j), iim
144	cof = rlamda(k) * cos(rlatu(j))
145	coefilu(k, j) = cof - 1.
146	coefilu2(k, j) = cof**2 - 1.
147	end DO
148	end if
149	END DO
150
151	DO j = 1, jfiltnv
152	modfrstnv(j) = 2
153	do while (rlamda(modfrstnv(j)) * cos(rlatv(j)) >= 1. &
154	.and. modfrstnv(j) <= iim - 1)
155	modfrstnv(j) = modfrstnv(j) + 1
156	end do
157
158	if (rlamda(modfrstnv(j)) * cos(rlatv(j)) < 1.) then
159	DO k = modfrstnv(j), iim
160	cof = rlamda(k) * cos(rlatv(j))
161	coefilv(k, j) = cof - 1.
162	coefilv2(k, j) = cof**2 - 1.
163	end DO
164	end if
165	end DO
166
167	DO j = jfiltsu, jjm
168	modfrstsu(j) = 2
169	do while (rlamda(modfrstsu(j)) * cos(rlatu(j)) >= 1. &
170	.and. modfrstsu(j) <= iim - 1)
171	modfrstsu(j) = modfrstsu(j) + 1
172	end do
173
174	if (rlamda(modfrstsu(j)) * cos(rlatu(j)) < 1.) then
175	DO k = modfrstsu(j), 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 if
181	end DO
182
183	DO j = jfiltsv, jjm
184	modfrstsv(j) = 2
185	do while (rlamda(modfrstsv(j)) * cos(rlatv(j)) >= 1. &
186	.and. modfrstsv(j) <= iim - 1)
187	modfrstsv(j) = modfrstsv(j) + 1
188	end do
189
190	if (rlamda(modfrstsv(j)) * cos(rlatv(j)) < 1.) then
191	DO k = modfrstsv(j), iim
192	cof = rlamda(k) * cos(rlatv(j))
193	coefilv(k, j) = cof - 1.
194	coefilv2(k, j) = cof**2 - 1.
195	end DO
196	end if
197	END DO
198
199	call new_unit(unit)
200	open(unit, file = "inifilr_out.txt", status = "replace", action = "write")
201	write(unit, fmt = *) '"EIGNVL"', eignvl
202	write(unit, fmt = *) '"modfrstnu"', modfrstnu
203	write(unit, fmt = *) '"modfrstsu"', modfrstsu
204	write(unit, fmt = *) '"modfrstnv"', modfrstnv
205	write(unit, fmt = *) '"modfrstsv"', modfrstsv
206	close(unit)
207
208	allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
209	allocate(matricevn(iim, iim, jfiltnv))
210	allocate(matricevs(iim, iim, jfiltsv:jjm))
211	allocate(matriceus(iim, iim, jfiltsu:jjm), matrinvs(iim, iim, jfiltsu:jjm))
212
213	! Calcul de la matrice filtre 'matriceu' pour les champs situes
214	! sur la grille scalaire
215
216	DO j = 2, jfiltnu
217	DO i = 1, iim
218	IF (i < modfrstnu(j)) then
219	coff = 0.
220	else
221	coff = coefilu(i, j)
222	end IF
223	eignft(i, :) = eignfnv(:, i) * coff
224	END DO
225	matriceun(:, :, j) = matmul(eignfnv, eignft)
226	END DO
227
228	DO j = jfiltsu, jjm
229	DO i = 1, iim
230	IF (i < modfrstsu(j)) then
231	coff = 0.
232	else
233	coff = coefilu(i, j)
234	end IF
235	eignft(i, :) = eignfnv(:, i) * coff
236	END DO
237	matriceus(:, :, j) = matmul(eignfnv, eignft)
238	END DO
239
240	! Calcul de la matrice filtre 'matricev' pour les champs situes
241	! sur la grille de V ou de Z
242
243	DO j = 1, jfiltnv
244	DO i = 1, iim
245	IF (i < modfrstnv(j)) then
246	coff = 0.
247	else
248	coff = coefilv(i, j)
249	end IF
250	eignft(i, :) = eignfnu(:, i) * coff
251	END DO
252	matricevn(:, :, j) = matmul(eignfnu, eignft)
253	END DO
254
255	DO j = jfiltsv, jjm
256	DO i = 1, iim
257	IF (i < modfrstsv(j)) then
258	coff = 0.
259	else
260	coff = coefilv(i, j)
261	end IF
262	eignft(i, :) = eignfnu(:, i) * coff
263	END DO
264	matricevs(:, :, j) = matmul(eignfnu, eignft)
265	END DO
266
267	! Calcul de la matrice filtre 'matrinv' pour les champs situes
268	! sur la grille scalaire , pour le filtre inverse
269
270	DO j = 2, jfiltnu
271	DO i = 1, iim
272	IF (i < modfrstnu(j)) then
273	coff = 0.
274	else
275	coff = coefilu(i, j) / (1. + coefilu(i, j))
276	end IF
277	eignft(i, :) = eignfnv(:, i) * coff
278	END DO
279	matrinvn(:, :, j) = matmul(eignfnv, eignft)
280	END DO
281
282	DO j = jfiltsu, jjm
283	DO i = 1, iim
284	IF (i < modfrstsu(j)) then
285	coff = 0.
286	else
287	coff = coefilu(i, j) / (1. + coefilu(i, j))
288	end IF
289	eignft(i, :) = eignfnv(:, i) * coff
290	END DO
291	matrinvs(:, :, j) = matmul(eignfnv, eignft)
292	END DO
293
294	END SUBROUTINE inifilr
295
296	end module inifilr_m