trunk/filtrez/inifilr.f

module inifilr_m

  IMPLICIT NONE

  ! North:

  INTEGER jfiltnu, jfiltnv
  ! index of the last scalar line filtered in northern hemisphere

  real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :) 
  ! (iim, iim, 2:jfiltnu)

  real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)

  ! South:

  integer jfiltsu, jfiltsv
  ! index of the first line filtered in southern hemisphere

  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 procedure computes the filtering coefficients for scalar
    ! lines and meridional wind v lines. The modes are filtered from
    ! modfrst to iim. 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.

    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 (<= 0)
    INTEGER i, j, unit
    REAL colat0 ! > 0
    REAL eignft(iim, iim)

    real eignfnu(iim, iim), eignfnv(iim, iim)
    ! eigenvectors of the discrete second derivative with respect to longitude

    ! Filtering coefficients (lamda_max * cos(rlat) / lamda):
    real, allocatable:: coefilnu(:, :) ! (iim, 2:jfiltnu)
    real, allocatable:: coefilsu(:, :) ! (iim, jfiltsu:jjm)
    real, allocatable:: coefilnv(:, :) ! (iim, jfiltnv)
    real, allocatable:: coefilsv(:, :) ! (iim, jfiltsv:jjm)

    ! Index of the mode from where modes are filtered:
    integer, allocatable:: modfrstnu(:) ! (2:jfiltnu)
    integer, allocatable:: modfrstsu(:) ! (jfiltsu:jjm)
    integer, allocatable:: modfrstnv(:) ! (jfiltnv)
    integer, allocatable:: modfrstsv(:) ! (jfiltsv:jjm)

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

    print *, "Call sequence information: inifilr"

    CALL inifgn(eignvl, eignfnu, eignfnv)

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

    ! D\'etermination de coefil[ns][uv], modfrst[ns][uv]:

    allocate(modfrstnu(2:jfiltnu), modfrstsu(jfiltsu:jjm))
    allocate(modfrstnv(jfiltnv), modfrstsv(jfiltsv:jjm))
    allocate(coefilnu(iim, 2:jfiltnu), coefilsu(iim, jfiltsu:jjm))
    allocate(coefilnv(iim, jfiltnv), coefilsv(iim, jfiltsv:jjm))

    coefilnu = 0.
    coefilnv = 0.
    coefilsu = 0.
    coefilsv = 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 i = modfrstnu(j), iim
             coefilnu(i, j) = rlamda(i) * cos(rlatu(j)) - 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 i = modfrstnv(j), iim
             coefilnv(i, j) = rlamda(i) * cos(rlatv(j)) - 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 i = modfrstsu(j), iim
             coefilsu(i, j) = rlamda(i) * cos(rlatu(j)) - 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 i = modfrstsv(j), iim
             coefilsv(i, j) = rlamda(i) * cos(rlatv(j)) - 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
       eignft(:modfrstnu(j) - 1, :) = 0.
       forall (i = modfrstnu(j):iim) eignft(i, :) = eignfnv(:, i) &
            * coefilnu(i, j)
       matriceun(:, :, j) = matmul(eignfnv, eignft)
    END DO

    DO j = jfiltsu, jjm
       eignft(:modfrstsu(j) - 1, :) = 0.
       forall (i = modfrstsu(j):iim) eignft(i, :) = eignfnv(:, i) &
            * coefilsu(i, j)
       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
       eignft(:modfrstnv(j) - 1, :) = 0.
       forall (i = modfrstnv(j): iim) eignft(i, :) = eignfnu(:, i) &
            * coefilnv(i, j)
       matricevn(:, :, j) = matmul(eignfnu, eignft)
    END DO

    DO j = jfiltsv, jjm
       eignft(:modfrstsv(j) - 1, :) = 0.
       forall (i = modfrstsv(j):iim) eignft(i, :) = eignfnu(:, i) &
            * coefilsv(i, j)
       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
       eignft(:modfrstnu(j) - 1, :) = 0.
       forall (i = modfrstnu(j):iim) eignft(i, :) = eignfnv(:, i) &
            * coefilnu(i, j) / (1. + coefilnu(i, j))
       matrinvn(:, :, j) = matmul(eignfnv, eignft)
    END DO

    DO j = jfiltsu, jjm
       eignft(:modfrstsu(j) - 1, :) = 0.
       forall (i = modfrstsu(j):iim) eignft(i, :) = eignfnv(:, i) &
            * coefilsu(i, j) / (1. + coefilsu(i, j))
       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	136	! North:
6	guez	156
7	guez	164	INTEGER jfiltnu, jfiltnv
8			! index of the last scalar line filtered in northern hemisphere
9
10	guez	136	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	guez	156
17	guez	164	integer jfiltsu, jfiltsv
18			! index of the first line filtered in southern hemisphere
19
20	guez	136	real, allocatable:: matriceus(:, :, :), matrinvs(:, :, :)
21			! (iim, iim, jfiltsu:jjm)
22
23			real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
24
25	guez	54	contains
26	guez	3
27	guez	54	SUBROUTINE inifilr
28	guez	3
29	guez	164	! From filtrez/inifilr.F, version 1.1.1.1, 2004/05/19 12:53:09
30	guez	54	! H. Upadhyaya, O. Sharma
31	guez	3
32	guez	164	! This procedure computes the filtering coefficients for scalar
33			! lines and meridional wind v lines. The modes are filtered from
34			! modfrst to iim. We filter all those latitude lines where coefil
35			! < 1. No filtering at poles. colat0 is to be used when alpha
36			! (stretching coefficient) is set equal to zero for the regular
37			! grid case.
38	guez	3
39	guez	54	USE dimens_m, ONLY : iim, jjm
40	guez	139	USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
41	guez	154	use inifgn_m, only: inifgn
42	guez	143	use jumble, only: new_unit
43	guez	54	use nr_util, only: pi
44	guez	3
45	guez	54	! Local:
46	guez	164
47	guez	132	REAL dlatu(jjm)
48	guez	140	REAL rlamda(2: iim)
49	guez	164	real eignvl(iim) ! eigenvalues sorted in descending order (<= 0)
50			INTEGER i, j, unit
51	guez	151	REAL colat0 ! > 0
52	guez	164	REAL eignft(iim, iim)
53	guez	3
54	guez	154	real eignfnu(iim, iim), eignfnv(iim, iim)
55	guez	164	! eigenvectors of the discrete second derivative with respect to longitude
56	guez	154
57	guez	140	! Filtering coefficients (lamda_max * cos(rlat) / lamda):
58	guez	164	real, allocatable:: coefilnu(:, :) ! (iim, 2:jfiltnu)
59			real, allocatable:: coefilsu(:, :) ! (iim, jfiltsu:jjm)
60			real, allocatable:: coefilnv(:, :) ! (iim, jfiltnv)
61			real, allocatable:: coefilsv(:, :) ! (iim, jfiltsv:jjm)
62	guez	140
63	guez	151	! Index of the mode from where modes are filtered:
64	guez	164	integer, allocatable:: modfrstnu(:) ! (2:jfiltnu)
65			integer, allocatable:: modfrstsu(:) ! (jfiltsu:jjm)
66			integer, allocatable:: modfrstnv(:) ! (jfiltnv)
67			integer, allocatable:: modfrstsv(:) ! (jfiltsv:jjm)
68	guez	140
69	guez	54	!-----------------------------------------------------------
70	guez	3
71	guez	54	print *, "Call sequence information: inifilr"
72	guez	3
73	guez	154	CALL inifgn(eignvl, eignfnu, eignfnv)
74	guez	3
75	guez	54	! Calcul de colat0
76	guez	143	forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
77			colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
78	guez	54	PRINT *, 'colat0 = ', colat0
79	guez	3
80	guez	164	rlamda = iim / (pi * colat0 / grossismx) / sqrt(- eignvl(2: iim))
81	guez	3
82	guez	151	! Determination de jfiltnu, jfiltsu, jfiltnv, jfiltsv
83	guez	3
84	guez	151	jfiltnu = (jjm + 1) / 2
85			do while (cos(rlatu(jfiltnu)) >= colat0 &
86			.or. rlamda(iim) * cos(rlatu(jfiltnu)) >= 1.)
87			jfiltnu = jfiltnu - 1
88			end do
89	guez	3
90	guez	151	jfiltsu = jjm / 2 + 2
91			do while (cos(rlatu(jfiltsu)) >= colat0 &
92			.or. rlamda(iim) * cos(rlatu(jfiltsu)) >= 1.)
93			jfiltsu = jfiltsu + 1
94			end do
95	guez	3
96	guez	151	jfiltnv = jjm / 2
97			do while ((cos(rlatv(jfiltnv)) >= colat0 &
98			.or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) .and. jfiltnv >= 2)
99			jfiltnv = jfiltnv - 1
100			end do
101	guez	3
102	guez	151	if (cos(rlatv(jfiltnv)) >= colat0 &
103			.or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) then
104			! {jfiltnv == 1}
105			PRINT *, 'Could not find jfiltnv.'
106	guez	54	STOP 1
107			END IF
108	guez	3
109	guez	151	jfiltsv = (jjm + 1)/ 2 + 1
110			do while ((cos(rlatv(jfiltsv)) >= colat0 &
111			.or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) .and. jfiltsv <= jjm - 1)
112			jfiltsv = jfiltsv + 1
113			end do
114	guez	3
115	guez	151	IF (cos(rlatv(jfiltsv)) >= colat0 &
116			.or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) THEN
117			! {jfiltsv == jjm}
118			PRINT *, 'Could not find jfiltsv.'
119	guez	54	STOP 1
120			END IF
121	guez	3
122	guez	151	PRINT *, 'jfiltnu =', jfiltnu
123			PRINT *, 'jfiltsu =', jfiltsu
124			PRINT *, 'jfiltnv =', jfiltnv
125			PRINT *, 'jfiltsv =', jfiltsv
126	guez	3
127	guez	164	! D\'etermination de coefil[ns][uv], modfrst[ns][uv]:
128	guez	32
129	guez	151	allocate(modfrstnu(2:jfiltnu), modfrstsu(jfiltsu:jjm))
130			allocate(modfrstnv(jfiltnv), modfrstsv(jfiltsv:jjm))
131	guez	164	allocate(coefilnu(iim, 2:jfiltnu), coefilsu(iim, jfiltsu:jjm))
132			allocate(coefilnv(iim, jfiltnv), coefilsv(iim, jfiltsv:jjm))
133	guez	32
134	guez	164	coefilnu = 0.
135			coefilnv = 0.
136			coefilsu = 0.
137			coefilsv = 0.
138
139	guez	54	DO j = 2, jfiltnu
140	guez	151	modfrstnu(j) = 2
141			do while (rlamda(modfrstnu(j)) * cos(rlatu(j)) >= 1. &
142			.and. modfrstnu(j) <= iim - 1)
143			modfrstnu(j) = modfrstnu(j) + 1
144			end do
145	guez	32
146	guez	151	if (rlamda(modfrstnu(j)) * cos(rlatu(j)) < 1.) then
147	guez	164	DO i = modfrstnu(j), iim
148			coefilnu(i, j) = rlamda(i) * cos(rlatu(j)) - 1.
149	guez	151	end DO
150			end if
151	guez	54	END DO
152	guez	32
153	guez	54	DO j = 1, jfiltnv
154	guez	151	modfrstnv(j) = 2
155			do while (rlamda(modfrstnv(j)) * cos(rlatv(j)) >= 1. &
156			.and. modfrstnv(j) <= iim - 1)
157			modfrstnv(j) = modfrstnv(j) + 1
158			end do
159	guez	32
160	guez	151	if (rlamda(modfrstnv(j)) * cos(rlatv(j)) < 1.) then
161	guez	164	DO i = modfrstnv(j), iim
162			coefilnv(i, j) = rlamda(i) * cos(rlatv(j)) - 1.
163	guez	151	end DO
164			end if
165	guez	54	end DO
166	guez	32
167	guez	54	DO j = jfiltsu, jjm
168	guez	151	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	guez	32
174	guez	151	if (rlamda(modfrstsu(j)) * cos(rlatu(j)) < 1.) then
175	guez	164	DO i = modfrstsu(j), iim
176			coefilsu(i, j) = rlamda(i) * cos(rlatu(j)) - 1.
177	guez	151	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	guez	164	DO i = modfrstsv(j), iim
190			coefilsv(i, j) = rlamda(i) * cos(rlatv(j)) - 1.
191	guez	151	end DO
192			end if
193	guez	54	END DO
194	guez	32
195	guez	151	call new_unit(unit)
196			open(unit, file = "inifilr_out.txt", status = "replace", action = "write")
197			write(unit, fmt = *) '"EIGNVL"', eignvl
198			write(unit, fmt = *) '"modfrstnu"', modfrstnu
199			write(unit, fmt = *) '"modfrstsu"', modfrstsu
200			write(unit, fmt = *) '"modfrstnv"', modfrstnv
201			write(unit, fmt = *) '"modfrstsv"', modfrstsv
202			close(unit)
203	guez	32
204	guez	136	allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
205			allocate(matricevn(iim, iim, jfiltnv))
206			allocate(matricevs(iim, iim, jfiltsv:jjm))
207			allocate(matriceus(iim, iim, jfiltsu:jjm), matrinvs(iim, iim, jfiltsu:jjm))
208	guez	32
209	guez	54	! Calcul de la matrice filtre 'matriceu' pour les champs situes
210			! sur la grille scalaire
211	guez	32
212	guez	54	DO j = 2, jfiltnu
213	guez	164	eignft(:modfrstnu(j) - 1, :) = 0.
214			forall (i = modfrstnu(j):iim) eignft(i, :) = eignfnv(:, i) &
215			* coefilnu(i, j)
216	guez	54	matriceun(:, :, j) = matmul(eignfnv, eignft)
217			END DO
218	guez	32
219	guez	54	DO j = jfiltsu, jjm
220	guez	164	eignft(:modfrstsu(j) - 1, :) = 0.
221			forall (i = modfrstsu(j):iim) eignft(i, :) = eignfnv(:, i) &
222			* coefilsu(i, j)
223	guez	136	matriceus(:, :, j) = matmul(eignfnv, eignft)
224	guez	54	END DO
225	guez	32
226	guez	54	! Calcul de la matrice filtre 'matricev' pour les champs situes
227			! sur la grille de V ou de Z
228	guez	32
229	guez	54	DO j = 1, jfiltnv
230	guez	164	eignft(:modfrstnv(j) - 1, :) = 0.
231			forall (i = modfrstnv(j): iim) eignft(i, :) = eignfnu(:, i) &
232			* coefilnv(i, j)
233	guez	54	matricevn(:, :, j) = matmul(eignfnu, eignft)
234			END DO
235	guez	32
236	guez	54	DO j = jfiltsv, jjm
237	guez	164	eignft(:modfrstsv(j) - 1, :) = 0.
238			forall (i = modfrstsv(j):iim) eignft(i, :) = eignfnu(:, i) &
239			* coefilsv(i, j)
240	guez	136	matricevs(:, :, j) = matmul(eignfnu, eignft)
241	guez	54	END DO
242	guez	32
243	guez	54	! Calcul de la matrice filtre 'matrinv' pour les champs situes
244			! sur la grille scalaire , pour le filtre inverse
245	guez	32
246	guez	54	DO j = 2, jfiltnu
247	guez	164	eignft(:modfrstnu(j) - 1, :) = 0.
248			forall (i = modfrstnu(j):iim) eignft(i, :) = eignfnv(:, i) &
249			* coefilnu(i, j) / (1. + coefilnu(i, j))
250	guez	54	matrinvn(:, :, j) = matmul(eignfnv, eignft)
251			END DO
252	guez	32
253	guez	54	DO j = jfiltsu, jjm
254	guez	164	eignft(:modfrstsu(j) - 1, :) = 0.
255			forall (i = modfrstsu(j):iim) eignft(i, :) = eignfnv(:, i) &
256			* coefilsu(i, j) / (1. + coefilsu(i, j))
257	guez	136	matrinvs(:, :, j) = matmul(eignfnv, eignft)
258	guez	54	END DO
259	guez	32
260	guez	54	END SUBROUTINE inifilr
261	guez	32
262	guez	54	end module inifilr_m