trunk/filtrez/inifilr.f

module inifilr_m

  IMPLICIT NONE

  INTEGER jfiltnu, jfiltnv
  ! index of the last line filtered in northern hemisphere at rlat[uv]
  ! latitudes

  integer jfiltsu, jfiltsv
  ! index of the first line filtered in southern hemisphere at
  ! rlat[uv] latitudes

  ! Filtre pour les champs situes sur la grille scalaire (longitudes
  ! rlonv, latitudes rlatu) :
  real, pointer:: matriceun(:, :, :) ! (iim, iim, jfiltnu - 1)
  real, pointer:: matriceus(:, :, :) ! (iim, iim, jjm - jfiltsu + 1)

  ! Filtre pour les champs situes sur la grille scalaire (longitudes
  ! rlonv, latitudes rlatu), pour le filtre inverse :
  real, pointer:: matrinvn(:, :, :) ! (iim, iim, jfiltnu - 1)
  real, pointer:: matrinvs(:, :, :) ! (iim, iim, jjm - jfiltsu + 1)

  ! Filtre pour les champs situes sur la grille de la vorticit\'e
  ! (longitudes rlonu, latitudes rlatv)
  real, pointer:: matricevn(:, :, :) ! (iim, iim, jfiltnv) matrice
  real, pointer:: matricevs(:, :, :) ! (iim, iim, jjm - jfiltsv + 1)

contains

  SUBROUTINE inifilr

    ! From filtrez/inifilr.F, version 1.1.1.1, 2004/05/19 12:53:09
    ! H. Upadhyaya, O. Sharma

    ! This procedure computes the filtering coefficients for scalar
    ! lines and meridional wind v lines. The modes are filtered from
    ! modfrst to iim. 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.

    USE dimens_m, ONLY : iim, jjm
    USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
    use inifgn_m, only: inifgn
    use inifilr_hemisph_m, only: inifilr_hemisph
    use jumble, only: new_unit
    use nr_util, only: pi, ifirstloc, assert

    ! Local:

    REAL dlatu(jjm)
    REAL rlamda(2:iim) ! > 0, in descending order
    real eignvl(iim) ! eigenvalues (<= 0) sorted in descending order
    INTEGER j, unit
    REAL colat0 ! > 0
    integer j1 ! index of negative latitude closest to the equator

    real eignfnu(iim, iim), eignfnv(iim, iim)
    ! eigenvectors of the discrete second derivative with respect to
    ! longitude, at rlon[uv] longitudes

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

    print *, "Call sequence information: inifilr"

    CALL inifgn(eignvl, eignfnu, eignfnv)

    ! 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(- eignvl(2: iim))
    print *, "1 / rlamda(iim) = ", 1. / rlamda(iim)

    ! This is demonstrated in the notes but just to be sure:
    call assert(rlamda(iim) * colat0 >= 1. - epsilon(0.), &
         "inifilr rlamda(iim) * colat0")

    call new_unit(unit)
    open(unit, file = "modfrst.csv", status = "replace", action = "write") 
    write(unit, fmt = *) '"rlat (degrees)" modfrst' ! title line

    j1 = ifirstloc(rlatu <= 0.)

    call inifilr_hemisph(rlatu(j1 - 1:2:- 1), rlamda, unit, eignfnv, jfiltnu, &
         matriceun, matrinvn)
    jfiltnu = j1 - jfiltnu
    matriceun = matriceun(:, :, jfiltnu - 1:1:- 1)
    matrinvn = matrinvn(:, :, jfiltnu - 1:1:- 1)

    call inifilr_hemisph(- rlatu(j1:jjm), rlamda, unit, eignfnv, jfiltsu, &
         matriceus, matrinvs)
    jfiltsu = j1 - 1 + jfiltsu

    j1 = ifirstloc(rlatv <= 0.)

    call inifilr_hemisph(rlatv(j1 - 1:1:- 1), rlamda, unit, eignfnu, jfiltnv, &
         matricevn)
    jfiltnv = j1 - jfiltnv
    matricevn = matricevn(:, :, jfiltnv:1:- 1)

    call inifilr_hemisph(- rlatv(j1:jjm), rlamda, unit, eignfnu, jfiltsv, &
         matricevs)
    jfiltsv = j1 - 1 + jfiltsv

    close(unit)
    PRINT *, 'jfiltnu =', jfiltnu
    PRINT *, 'jfiltsu =', jfiltsu
    PRINT *, 'jfiltnv =', jfiltnv
    PRINT *, 'jfiltsv =', jfiltsv

  END SUBROUTINE inifilr

end module inifilr_m
1	module inifilr_m
2
3	IMPLICIT NONE
4
5	INTEGER jfiltnu, jfiltnv
6	! index of the last line filtered in northern hemisphere at rlat[uv]
7	! latitudes
8
9	integer jfiltsu, jfiltsv
10	! index of the first line filtered in southern hemisphere at
11	! rlat[uv] latitudes
12
13	! Filtre pour les champs situes sur la grille scalaire (longitudes
14	! rlonv, latitudes rlatu) :
15	real, pointer:: matriceun(:, :, :) ! (iim, iim, jfiltnu - 1)
16	real, pointer:: matriceus(:, :, :) ! (iim, iim, jjm - jfiltsu + 1)
17
18	! Filtre pour les champs situes sur la grille scalaire (longitudes
19	! rlonv, latitudes rlatu), pour le filtre inverse :
20	real, pointer:: matrinvn(:, :, :) ! (iim, iim, jfiltnu - 1)
21	real, pointer:: matrinvs(:, :, :) ! (iim, iim, jjm - jfiltsu + 1)
22
23	! Filtre pour les champs situes sur la grille de la vorticit\'e
24	! (longitudes rlonu, latitudes rlatv)
25	real, pointer:: matricevn(:, :, :) ! (iim, iim, jfiltnv) matrice
26	real, pointer:: matricevs(:, :, :) ! (iim, iim, jjm - jfiltsv + 1)
27
28	contains
29
30	SUBROUTINE inifilr
31
32	! From filtrez/inifilr.F, version 1.1.1.1, 2004/05/19 12:53:09
33	! H. Upadhyaya, O. Sharma
34
35	! This procedure computes the filtering coefficients for scalar
36	! lines and meridional wind v lines. The modes are filtered from
37	! modfrst to iim. We filter all those latitude lines where coefil
38	! < 1. No filtering at poles. colat0 is to be used when alpha
39	! (stretching coefficient) is set equal to zero for the regular
40	! grid case.
41
42	USE dimens_m, ONLY : iim, jjm
43	USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
44	use inifgn_m, only: inifgn
45	use inifilr_hemisph_m, only: inifilr_hemisph
46	use jumble, only: new_unit
47	use nr_util, only: pi, ifirstloc, assert
48
49	! Local:
50
51	REAL dlatu(jjm)
52	REAL rlamda(2:iim) ! > 0, in descending order
53	real eignvl(iim) ! eigenvalues (<= 0) sorted in descending order
54	INTEGER j, unit
55	REAL colat0 ! > 0
56	integer j1 ! index of negative latitude closest to the equator
57
58	real eignfnu(iim, iim), eignfnv(iim, iim)
59	! eigenvectors of the discrete second derivative with respect to
60	! longitude, at rlon[uv] longitudes
61
62	!-----------------------------------------------------------
63
64	print *, "Call sequence information: inifilr"
65
66	CALL inifgn(eignvl, eignfnu, eignfnv)
67
68	! Calcul de colat0
69	forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
70	colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
71	PRINT *, 'colat0 = ', colat0
72
73	rlamda = iim / pi / colat0 * grossismx / sqrt(- eignvl(2: iim))
74	print *, "1 / rlamda(iim) = ", 1. / rlamda(iim)
75
76	! This is demonstrated in the notes but just to be sure:
77	call assert(rlamda(iim) * colat0 >= 1. - epsilon(0.), &
78	"inifilr rlamda(iim) * colat0")
79
80	call new_unit(unit)
81	open(unit, file = "modfrst.csv", status = "replace", action = "write")
82	write(unit, fmt = *) '"rlat (degrees)" modfrst' ! title line
83
84	j1 = ifirstloc(rlatu <= 0.)
85
86	call inifilr_hemisph(rlatu(j1 - 1:2:- 1), rlamda, unit, eignfnv, jfiltnu, &
87	matriceun, matrinvn)
88	jfiltnu = j1 - jfiltnu
89	matriceun = matriceun(:, :, jfiltnu - 1:1:- 1)
90	matrinvn = matrinvn(:, :, jfiltnu - 1:1:- 1)
91
92	call inifilr_hemisph(- rlatu(j1:jjm), rlamda, unit, eignfnv, jfiltsu, &
93	matriceus, matrinvs)
94	jfiltsu = j1 - 1 + jfiltsu
95
96	j1 = ifirstloc(rlatv <= 0.)
97
98	call inifilr_hemisph(rlatv(j1 - 1:1:- 1), rlamda, unit, eignfnu, jfiltnv, &
99	matricevn)
100	jfiltnv = j1 - jfiltnv
101	matricevn = matricevn(:, :, jfiltnv:1:- 1)
102
103	call inifilr_hemisph(- rlatv(j1:jjm), rlamda, unit, eignfnu, jfiltsv, &
104	matricevs)
105	jfiltsv = j1 - 1 + jfiltsv
106
107	close(unit)
108	PRINT *, 'jfiltnu =', jfiltnu
109	PRINT *, 'jfiltsu =', jfiltsu
110	PRINT *, 'jfiltnv =', jfiltnv
111	PRINT *, 'jfiltsv =', jfiltsv
112
113	END SUBROUTINE inifilr
114
115	end module inifilr_m