Sources/filtrez/inifilr_hemisph.f

module inifilr_hemisph_m

  implicit none

contains

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

    ! See notes.

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

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

    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_lat - jfilt + 1)
    ! matrice filtre

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

    ! Local:

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

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

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

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

    n_lat = size(rlat)
    jfilt = ifirstloc(cos(rlat) < 1. / rlamda(iim))
    allocate(matrice(iim, iim, n_lat - jfilt + 1))
    if (present(matrinv)) allocate(matrinv(iim, iim, n_lat - jfilt + 1))

    DO j = jfilt, n_lat
       modfrst = ifirstloc(rlamda < 1. / cos(rlat(j)), my_lbound = 2)
       write(unit, fmt = *) rlat(j) / pi * 180., modfrst
       coefil(modfrst:) = rlamda(modfrst:) * cos(rlat(j)) - 1.
       eignft(:modfrst - 1, :) = 0.

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

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