trunk/filtrez/inifilr_hemisph.f

module inifilr_hemisph_m

  implicit none

contains

  subroutine inifilr_hemisph(rlat, colat0, rlamda, unit, eignfn, jfilt, &
       matrice, matrinv)

    ! See notes.

    USE dimens_m, ONLY : iim
    use nr_util, only: pi, ifirstloc

    real, intent(in):: rlat(:) ! (n)
    ! latitudes, in rad, in [0, pi / 2[, in strictly ascending order

    REAL, intent(in):: colat0 ! > 0
    REAL, intent(in):: rlamda(2:) ! (2:iim) > 0, in descending order
    integer, intent(in):: unit

    real, intent(in):: eignfn(:, :) ! (iim, iim)
    ! eigenvectors of the discrete second derivative with respect to longitude

    integer, intent(out):: jfilt
    real, pointer:: matrice(:, :, :) ! (iim, iim, n - jfilt + 1) matrice filtre

    real, pointer, optional:: matrinv(:, :, :) ! (iim, iim, n - jfilt + 1) 
    ! matrice pour le filtre inverse

    ! Local:

    integer n, i, j
    REAL eignft(iim, iim)

    ! Index of the mode from where modes are filtered:
    integer, allocatable:: modfrst(:) ! (jfilt:n) in {2, ..., iim}

    ! Filtering coefficients (lamda_max * cos(rlat) / lamda):
    real coefil(2:iim)

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

    n = size(rlat)
    jfilt = ifirstloc(cos(rlat) < min(colat0, 1. / rlamda(iim)))
    allocate(modfrst(jfilt:n))

    DO j = jfilt, n
       modfrst(j) = ifirstloc(rlamda < 1. / cos(rlat(j)), my_lbound = 2)
       write(unit, fmt = *) rlat(j) / pi * 180., modfrst(j)
    END DO

    allocate(matrice(iim, iim, n - jfilt + 1))
    if (present(matrinv)) allocate(matrinv(iim, iim, n - jfilt + 1))

    DO j = jfilt, n
       coefil(modfrst(j):) = rlamda(modfrst(j):) * cos(rlat(j)) - 1.
       eignft(:modfrst(j) - 1, :) = 0.

       forall (i = modfrst(j):iim) eignft(i, :) = eignfn(:, i) * coefil(i)
       matrice(:, :, j - jfilt + 1) = matmul(eignfn, eignft)

       if (present(matrinv)) then
          forall (i = modfrst(j):iim) eignft(i, :) = eignfn(:, i) * coefil(i) &
               / (1. + coefil(i))
          matrinv(:, :, j - jfilt + 1) = matmul(eignfn, eignft)
       end if
    END DO

  end subroutine inifilr_hemisph

end module inifilr_hemisph_m
1	module inifilr_hemisph_m
2
3	implicit none
4
5	contains
6
7	subroutine inifilr_hemisph(rlat, colat0, rlamda, unit, eignfn, jfilt, &
8	matrice, matrinv)
9
10	! See notes.
11
12	USE dimens_m, ONLY : iim
13	use nr_util, only: pi, ifirstloc
14
15	real, intent(in):: rlat(:) ! (n)
16	! latitudes, in rad, in [0, pi / 2[, in strictly ascending order
17
18	REAL, intent(in):: colat0 ! > 0
19	REAL, intent(in):: rlamda(2:) ! (2:iim) > 0, in descending order
20	integer, intent(in):: unit
21
22	real, intent(in):: eignfn(:, :) ! (iim, iim)
23	! eigenvectors of the discrete second derivative with respect to longitude
24
25	integer, intent(out):: jfilt
26	real, pointer:: matrice(:, :, :) ! (iim, iim, n - jfilt + 1) matrice filtre
27
28	real, pointer, optional:: matrinv(:, :, :) ! (iim, iim, n - jfilt + 1)
29	! matrice pour le filtre inverse
30
31	! Local:
32
33	integer n, i, j
34	REAL eignft(iim, iim)
35
36	! Index of the mode from where modes are filtered:
37	integer, allocatable:: modfrst(:) ! (jfilt:n) in {2, ..., iim}
38
39	! Filtering coefficients (lamda_max * cos(rlat) / lamda):
40	real coefil(2:iim)
41
42	!-----------------------------------------------------------
43
44	n = size(rlat)
45	jfilt = ifirstloc(cos(rlat) < min(colat0, 1. / rlamda(iim)))
46	allocate(modfrst(jfilt:n))
47
48	DO j = jfilt, n
49	modfrst(j) = ifirstloc(rlamda < 1. / cos(rlat(j)), my_lbound = 2)
50	write(unit, fmt = ) rlat(j) / pi 180., modfrst(j)
51	END DO
52
53	allocate(matrice(iim, iim, n - jfilt + 1))
54	if (present(matrinv)) allocate(matrinv(iim, iim, n - jfilt + 1))
55
56	DO j = jfilt, n
57	coefil(modfrst(j):) = rlamda(modfrst(j):) * cos(rlat(j)) - 1.
58	eignft(:modfrst(j) - 1, :) = 0.
59
60	forall (i = modfrst(j):iim) eignft(i, :) = eignfn(:, i) * coefil(i)
61	matrice(:, :, j - jfilt + 1) = matmul(eignfn, eignft)
62
63	if (present(matrinv)) then
64	forall (i = modfrst(j):iim) eignft(i, :) = eignfn(:, i) * coefil(i) &
65	/ (1. + coefil(i))
66	matrinv(:, :, j - jfilt + 1) = matmul(eignfn, eignft)
67	end if
68	END DO
69
70	end subroutine inifilr_hemisph
71
72	end module inifilr_hemisph_m