trunk/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	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	guez	156
11	guez	136	real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :)
12			! (iim, iim, 2:jfiltnu)
13	guez	3
14	guez	136	real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)
15	guez	3
16	guez	136	! South:
17	guez	156
18	guez	136	real, allocatable:: matriceus(:, :, :), matrinvs(:, :, :)
19			! (iim, iim, jfiltsu:jjm)
20
21			real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
22
23	guez	54	contains
24	guez	3
25	guez	54	SUBROUTINE inifilr
26	guez	3
27	guez	54	! From filtrez/inifilr.F, version 1.1.1.1 2004/05/19 12:53:09
28			! H. Upadhyaya, O. Sharma
29	guez	3
30	guez	156	! This routine computes the eigenvectors of the laplacian on the
31	guez	140	! stretched grid, and the filtering coefficients. The modes are
32			! filtered from modfrst to iim.
33	guez	3
34	guez	54	USE dimens_m, ONLY : iim, jjm
35	guez	139	USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
36	guez	154	use inifgn_m, only: inifgn
37	guez	143	use jumble, only: new_unit
38	guez	54	use nr_util, only: pi
39	guez	3
40	guez	54	! Local:
41	guez	132	REAL dlatu(jjm)
42	guez	140	REAL rlamda(2: iim)
43	guez	151	real eignvl(iim) ! eigenvalues sorted in descending order
44	guez	143	REAL cof
45	guez	151	INTEGER i, j, k, unit
46			REAL colat0 ! > 0
47	guez	54	REAL eignft(iim, iim), coff
48	guez	3
49	guez	154	real eignfnu(iim, iim), eignfnv(iim, iim)
50	guez	156	! eigenvectors of the discrete laplacian
51	guez	154
52	guez	140	! 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	guez	151	! Index of the mode from where modes are filtered:
57			integer, allocatable:: modfrstnu(:), modfrstsu(:)
58			integer, allocatable:: modfrstnv(:), modfrstsv(:)
59	guez	140
60	guez	54	!-----------------------------------------------------------
61	guez	3
62	guez	54	print *, "Call sequence information: inifilr"
63	guez	3
64	guez	154	CALL inifgn(eignvl, eignfnu, eignfnv)
65	guez	3
66	guez	156	! compute eigenvalues and eigenvectors
67	guez	54	! 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	guez	3
74	guez	54	! Calcul de colat0
75	guez	143	forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
76			colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
77	guez	54	PRINT *, 'colat0 = ', colat0
78	guez	3
79	guez	143	rlamda = iim / (pi * colat0 / grossismx) / sqrt(abs(eignvl(2: iim)))
80	guez	3
81	guez	151	! Determination de jfiltnu, jfiltsu, jfiltnv, jfiltsv
82	guez	3
83	guez	151	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	guez	3
89	guez	151	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	guez	3
95	guez	151	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	guez	3
101	guez	151	if (cos(rlatv(jfiltnv)) >= colat0 &
102			.or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) then
103			! {jfiltnv == 1}
104			PRINT *, 'Could not find jfiltnv.'
105	guez	54	STOP 1
106			END IF
107	guez	3
108	guez	151	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	guez	3
114	guez	151	IF (cos(rlatv(jfiltsv)) >= colat0 &
115			.or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) THEN
116			! {jfiltsv == jjm}
117			PRINT *, 'Could not find jfiltsv.'
118	guez	54	STOP 1
119			END IF
120	guez	3
121	guez	151	PRINT *, 'jfiltnu =', jfiltnu
122			PRINT *, 'jfiltsu =', jfiltsu
123			PRINT *, 'jfiltnv =', jfiltnv
124			PRINT *, 'jfiltsv =', jfiltsv
125	guez	3
126	guez	151	! Determination de coefilu, coefilv, modfrst[ns][uv]:
127	guez	32
128	guez	151	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	guez	32
135	guez	54	DO j = 2, jfiltnu
136	guez	151	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	guez	32
142	guez	151	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	guez	54	END DO
150	guez	32
151	guez	54	DO j = 1, jfiltnv
152	guez	151	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	guez	32
158	guez	151	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	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			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	guez	54	end DO
182	guez	32
183	guez	54	DO j = jfiltsv, jjm
184	guez	151	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	guez	32
190	guez	151	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	guez	54	END DO
198	guez	32
199	guez	151	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	guez	32
208	guez	136	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	guez	32
213	guez	54	! Calcul de la matrice filtre 'matriceu' pour les champs situes
214			! sur la grille scalaire
215	guez	32
216	guez	54	DO j = 2, jfiltnu
217			DO i = 1, iim
218	guez	151	IF (i < modfrstnu(j)) then
219	guez	54	coff = 0.
220			else
221			coff = coefilu(i, j)
222			end IF
223	guez	140	eignft(i, :) = eignfnv(:, i) * coff
224	guez	54	END DO
225			matriceun(:, :, j) = matmul(eignfnv, eignft)
226			END DO
227	guez	32
228	guez	54	DO j = jfiltsu, jjm
229			DO i = 1, iim
230	guez	151	IF (i < modfrstsu(j)) then
231	guez	54	coff = 0.
232			else
233			coff = coefilu(i, j)
234			end IF
235			eignft(i, :) = eignfnv(:, i) * coff
236			END DO
237	guez	136	matriceus(:, :, j) = matmul(eignfnv, eignft)
238	guez	54	END DO
239	guez	32
240	guez	54	! Calcul de la matrice filtre 'matricev' pour les champs situes
241			! sur la grille de V ou de Z
242	guez	32
243	guez	54	DO j = 1, jfiltnv
244			DO i = 1, iim
245	guez	151	IF (i < modfrstnv(j)) then
246	guez	54	coff = 0.
247			else
248			coff = coefilv(i, j)
249			end IF
250	guez	140	eignft(i, :) = eignfnu(:, i) * coff
251	guez	54	END DO
252			matricevn(:, :, j) = matmul(eignfnu, eignft)
253			END DO
254	guez	32
255	guez	54	DO j = jfiltsv, jjm
256			DO i = 1, iim
257	guez	151	IF (i < modfrstsv(j)) then
258	guez	54	coff = 0.
259			else
260			coff = coefilv(i, j)
261			end IF
262	guez	140	eignft(i, :) = eignfnu(:, i) * coff
263	guez	54	END DO
264	guez	136	matricevs(:, :, j) = matmul(eignfnu, eignft)
265	guez	54	END DO
266	guez	32
267	guez	54	! Calcul de la matrice filtre 'matrinv' pour les champs situes
268			! sur la grille scalaire , pour le filtre inverse
269	guez	32
270	guez	54	DO j = 2, jfiltnu
271			DO i = 1, iim
272	guez	151	IF (i < modfrstnu(j)) then
273	guez	54	coff = 0.
274			else
275	guez	140	coff = coefilu(i, j) / (1. + coefilu(i, j))
276	guez	54	end IF
277	guez	140	eignft(i, :) = eignfnv(:, i) * coff
278	guez	54	END DO
279			matrinvn(:, :, j) = matmul(eignfnv, eignft)
280			END DO
281	guez	32
282	guez	54	DO j = jfiltsu, jjm
283			DO i = 1, iim
284	guez	151	IF (i < modfrstsu(j)) then
285	guez	54	coff = 0.
286			else
287	guez	140	coff = coefilu(i, j) / (1. + coefilu(i, j))
288	guez	54	end IF
289	guez	140	eignft(i, :) = eignfnv(:, i) * coff
290	guez	54	END DO
291	guez	136	matrinvs(:, :, j) = matmul(eignfnv, eignft)
292	guez	54	END DO
293	guez	32
294	guez	54	END SUBROUTINE inifilr
295	guez	32
296	guez	54	end module inifilr_m