Sources/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(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
       DO i = modfrst(j), iim
          coefil(i) = rlamda(i) * cos(rlat(j)) - 1.
       end DO

       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	guez	166	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(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			DO i = modfrst(j), iim
58			coefil(i) = rlamda(i) * cos(rlat(j)) - 1.
59			end DO
60
61			eignft(:modfrst(j) - 1, :) = 0.
62
63			forall (i = modfrst(j):iim) eignft(i, :) = eignfn(:, i) * coefil(i)
64			matrice(:, :, j - jfilt + 1) = matmul(eignfn, eignft)
65
66			if (present(matrinv)) then
67			forall (i = modfrst(j):iim) eignft(i, :) = eignfn(:, i) * coefil(i) &
68			/ (1. + coefil(i))
69			matrinv(:, :, j - jfilt + 1) = matmul(eignfn, eignft)
70			end if
71			END DO
72
73			end subroutine inifilr_hemisph
74
75			end module inifilr_hemisph_m