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