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	module inifilr_hemisph_m
2
3	implicit none
4
5	contains
6
7	subroutine inifilr_hemisph(rlat, rlamda, unit, eignfn, jfilt, matrice, &
8	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_lat)
16	! 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	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	! matrice pour le filtre inverse
31
32	! Local:
33
34	integer n_lat, i, j
35	REAL eignft(iim, iim)
36
37	! Index of the mode from where modes are filtered:
38	integer modfrst ! in {2, ..., iim}
39
40	! Filtering coefficients (lamda_max * cos(rlat) / lamda):
41	real coefil(2:iim)
42
43	!-----------------------------------------------------------
44
45	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
50	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
56	forall (i = modfrst:iim) eignft(i, :) = eignfn(:, i) * coefil(i)
57	matrice(:, :, j - jfilt + 1) = matmul(eignfn, eignft)
58
59	if (present(matrinv)) then
60	forall (i = modfrst:iim) eignft(i, :) = eignfn(:, i) * coefil(i) &
61	/ (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