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 sorted in descending order
    REAL cof
    INTEGER i, j, k, unit
    REAL colat0 ! > 0
    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)

    ! 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)

    ! 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)))

    ! 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	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 sorted in descending order
42	REAL cof
43	INTEGER i, j, k, unit
44	REAL colat0 ! > 0
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	! Index of the mode from where modes are filtered:
52	integer, allocatable:: modfrstnu(:), modfrstsu(:)
53	integer, allocatable:: modfrstnv(:), modfrstsv(:)
54
55	!-----------------------------------------------------------
56
57	print *, "Call sequence information: inifilr"
58
59	CALL inifgn(eignvl)
60
61	! compute eigenvalues and eigenfunctions
62	! compute the filtering coefficients for scalar lines and
63	! meridional wind v-lines
64	! we filter all those latitude lines where coefil < 1
65	! NO FILTERING AT POLES
66	! colat0 is to be used when alpha (stretching coefficient)
67	! is set equal to zero for the regular grid case
68
69	! Calcul de colat0
70	forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
71	colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
72	PRINT *, 'colat0 = ', colat0
73
74	rlamda = iim / (pi * colat0 / grossismx) / sqrt(abs(eignvl(2: iim)))
75
76	! Determination de jfiltnu, jfiltsu, jfiltnv, jfiltsv
77
78	jfiltnu = (jjm + 1) / 2
79	do while (cos(rlatu(jfiltnu)) >= colat0 &
80	.or. rlamda(iim) * cos(rlatu(jfiltnu)) >= 1.)
81	jfiltnu = jfiltnu - 1
82	end do
83
84	jfiltsu = jjm / 2 + 2
85	do while (cos(rlatu(jfiltsu)) >= colat0 &
86	.or. rlamda(iim) * cos(rlatu(jfiltsu)) >= 1.)
87	jfiltsu = jfiltsu + 1
88	end do
89
90	jfiltnv = jjm / 2
91	do while ((cos(rlatv(jfiltnv)) >= colat0 &
92	.or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) .and. jfiltnv >= 2)
93	jfiltnv = jfiltnv - 1
94	end do
95
96	if (cos(rlatv(jfiltnv)) >= colat0 &
97	.or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) then
98	! {jfiltnv == 1}
99	PRINT *, 'Could not find jfiltnv.'
100	STOP 1
101	END IF
102
103	jfiltsv = (jjm + 1)/ 2 + 1
104	do while ((cos(rlatv(jfiltsv)) >= colat0 &
105	.or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) .and. jfiltsv <= jjm - 1)
106	jfiltsv = jfiltsv + 1
107	end do
108
109	IF (cos(rlatv(jfiltsv)) >= colat0 &
110	.or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) THEN
111	! {jfiltsv == jjm}
112	PRINT *, 'Could not find jfiltsv.'
113	STOP 1
114	END IF
115
116	PRINT *, 'jfiltnu =', jfiltnu
117	PRINT *, 'jfiltsu =', jfiltsu
118	PRINT *, 'jfiltnv =', jfiltnv
119	PRINT *, 'jfiltsv =', jfiltsv
120
121	! Determination de coefilu, coefilv, modfrst[ns][uv]:
122
123	allocate(modfrstnu(2:jfiltnu), modfrstsu(jfiltsu:jjm))
124	allocate(modfrstnv(jfiltnv), modfrstsv(jfiltsv:jjm))
125	coefilu = 0.
126	coefilv = 0.
127	coefilu2 = 0.
128	coefilv2 = 0.
129
130	DO j = 2, jfiltnu
131	modfrstnu(j) = 2
132	do while (rlamda(modfrstnu(j)) * cos(rlatu(j)) >= 1. &
133	.and. modfrstnu(j) <= iim - 1)
134	modfrstnu(j) = modfrstnu(j) + 1
135	end do
136
137	if (rlamda(modfrstnu(j)) * cos(rlatu(j)) < 1.) then
138	DO k = modfrstnu(j), iim
139	cof = rlamda(k) * cos(rlatu(j))
140	coefilu(k, j) = cof - 1.
141	coefilu2(k, j) = cof**2 - 1.
142	end DO
143	end if
144	END DO
145
146	DO j = 1, jfiltnv
147	modfrstnv(j) = 2
148	do while (rlamda(modfrstnv(j)) * cos(rlatv(j)) >= 1. &
149	.and. modfrstnv(j) <= iim - 1)
150	modfrstnv(j) = modfrstnv(j) + 1
151	end do
152
153	if (rlamda(modfrstnv(j)) * cos(rlatv(j)) < 1.) then
154	DO k = modfrstnv(j), iim
155	cof = rlamda(k) * cos(rlatv(j))
156	coefilv(k, j) = cof - 1.
157	coefilv2(k, j) = cof**2 - 1.
158	end DO
159	end if
160	end DO
161
162	DO j = jfiltsu, jjm
163	modfrstsu(j) = 2
164	do while (rlamda(modfrstsu(j)) * cos(rlatu(j)) >= 1. &
165	.and. modfrstsu(j) <= iim - 1)
166	modfrstsu(j) = modfrstsu(j) + 1
167	end do
168
169	if (rlamda(modfrstsu(j)) * cos(rlatu(j)) < 1.) then
170	DO k = modfrstsu(j), iim
171	cof = rlamda(k) * cos(rlatu(j))
172	coefilu(k, j) = cof - 1.
173	coefilu2(k, j) = cof**2 - 1.
174	end DO
175	end if
176	end DO
177
178	DO j = jfiltsv, jjm
179	modfrstsv(j) = 2
180	do while (rlamda(modfrstsv(j)) * cos(rlatv(j)) >= 1. &
181	.and. modfrstsv(j) <= iim - 1)
182	modfrstsv(j) = modfrstsv(j) + 1
183	end do
184
185	if (rlamda(modfrstsv(j)) * cos(rlatv(j)) < 1.) then
186	DO k = modfrstsv(j), iim
187	cof = rlamda(k) * cos(rlatv(j))
188	coefilv(k, j) = cof - 1.
189	coefilv2(k, j) = cof**2 - 1.
190	end DO
191	end if
192	END DO
193
194	call new_unit(unit)
195	open(unit, file = "inifilr_out.txt", status = "replace", action = "write")
196	write(unit, fmt = *) '"EIGNVL"', eignvl
197	write(unit, fmt = *) '"modfrstnu"', modfrstnu
198	write(unit, fmt = *) '"modfrstsu"', modfrstsu
199	write(unit, fmt = *) '"modfrstnv"', modfrstnv
200	write(unit, fmt = *) '"modfrstsv"', modfrstsv
201	close(unit)
202
203	allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
204	allocate(matricevn(iim, iim, jfiltnv))
205	allocate(matricevs(iim, iim, jfiltsv:jjm))
206	allocate(matriceus(iim, iim, jfiltsu:jjm), matrinvs(iim, iim, jfiltsu:jjm))
207
208	! Calcul de la matrice filtre 'matriceu' pour les champs situes
209	! sur la grille scalaire
210
211	DO j = 2, jfiltnu
212	DO i = 1, iim
213	IF (i < modfrstnu(j)) then
214	coff = 0.
215	else
216	coff = coefilu(i, j)
217	end IF
218	eignft(i, :) = eignfnv(:, i) * coff
219	END DO
220	matriceun(:, :, j) = matmul(eignfnv, eignft)
221	END DO
222
223	DO j = jfiltsu, jjm
224	DO i = 1, iim
225	IF (i < modfrstsu(j)) then
226	coff = 0.
227	else
228	coff = coefilu(i, j)
229	end IF
230	eignft(i, :) = eignfnv(:, i) * coff
231	END DO
232	matriceus(:, :, j) = matmul(eignfnv, eignft)
233	END DO
234
235	! Calcul de la matrice filtre 'matricev' pour les champs situes
236	! sur la grille de V ou de Z
237
238	DO j = 1, jfiltnv
239	DO i = 1, iim
240	IF (i < modfrstnv(j)) then
241	coff = 0.
242	else
243	coff = coefilv(i, j)
244	end IF
245	eignft(i, :) = eignfnu(:, i) * coff
246	END DO
247	matricevn(:, :, j) = matmul(eignfnu, eignft)
248	END DO
249
250	DO j = jfiltsv, jjm
251	DO i = 1, iim
252	IF (i < modfrstsv(j)) then
253	coff = 0.
254	else
255	coff = coefilv(i, j)
256	end IF
257	eignft(i, :) = eignfnu(:, i) * coff
258	END DO
259	matricevs(:, :, j) = matmul(eignfnu, eignft)
260	END DO
261
262	! Calcul de la matrice filtre 'matrinv' pour les champs situes
263	! sur la grille scalaire , pour le filtre inverse
264
265	DO j = 2, jfiltnu
266	DO i = 1, iim
267	IF (i < modfrstnu(j)) then
268	coff = 0.
269	else
270	coff = coefilu(i, j) / (1. + coefilu(i, j))
271	end IF
272	eignft(i, :) = eignfnv(:, i) * coff
273	END DO
274	matrinvn(:, :, j) = matmul(eignfnv, eignft)
275	END DO
276
277	DO j = jfiltsu, jjm
278	DO i = 1, iim
279	IF (i < modfrstsu(j)) then
280	coff = 0.
281	else
282	coff = coefilu(i, j) / (1. + coefilu(i, j))
283	end IF
284	eignft(i, :) = eignfnv(:, i) * coff
285	END DO
286	matrinvs(:, :, j) = matmul(eignfnv, eignft)
287	END DO
288
289	END SUBROUTINE inifilr
290
291	end module inifilr_m