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.

    ! Libraries:
    use jumble, only: new_unit
    use nr_util, only: pi, ifirstloc, assert

    USE dimensions, ONLY: iim, jjm
    USE dynetat0_m, ONLY: rlatu, rlatv, xprimu
    USE dynetat0_chosen_m, ONLY: grossismx
    use inifgn_m, only: inifgn
    use inifilr_hemisph_m, only: inifilr_hemisph

    ! 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. - 2. * 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	! Libraries:
43	use jumble, only: new_unit
44	use nr_util, only: pi, ifirstloc, assert
45
46	USE dimensions, ONLY: iim, jjm
47	USE dynetat0_m, ONLY: rlatu, rlatv, xprimu
48	USE dynetat0_chosen_m, ONLY: grossismx
49	use inifgn_m, only: inifgn
50	use inifilr_hemisph_m, only: inifilr_hemisph
51
52	! Local:
53
54	REAL dlatu(jjm)
55	REAL rlamda(2:iim) ! > 0, in descending order
56	real eignvl(iim) ! eigenvalues (<= 0) sorted in descending order
57	INTEGER j, unit
58	REAL colat0 ! > 0
59	integer j1 ! index of negative latitude closest to the equator
60
61	real eignfnu(iim, iim), eignfnv(iim, iim)
62	! eigenvectors of the discrete second derivative with respect to
63	! longitude, at rlon[uv] longitudes
64
65	!-----------------------------------------------------------
66
67	print *, "Call sequence information: inifilr"
68
69	CALL inifgn(eignvl, eignfnu, eignfnv)
70
71	! Calcul de colat0
72	forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
73	colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
74	PRINT *, 'colat0 = ', colat0
75
76	rlamda = iim / pi / colat0 * grossismx / sqrt(- eignvl(2: iim))
77	print *, "1 / rlamda(iim) = ", 1. / rlamda(iim)
78	! This is demonstrated in the notes but just to be sure:
79	call assert(rlamda(iim) * colat0 >= 1. - 2. * epsilon(0.), &
80	"inifilr rlamda(iim) * colat0")
81
82	call new_unit(unit)
83	open(unit, file = "modfrst.csv", status = "replace", action = "write")
84	write(unit, fmt = *) '"rlat (degrees)" modfrst' ! title line
85
86	j1 = ifirstloc(rlatu <= 0.)
87
88	call inifilr_hemisph(rlatu(j1 - 1:2:- 1), rlamda, unit, eignfnv, jfiltnu, &
89	matriceun, matrinvn)
90	jfiltnu = j1 - jfiltnu
91	matriceun = matriceun(:, :, jfiltnu - 1:1:- 1)
92	matrinvn = matrinvn(:, :, jfiltnu - 1:1:- 1)
93
94	call inifilr_hemisph(- rlatu(j1:jjm), rlamda, unit, eignfnv, jfiltsu, &
95	matriceus, matrinvs)
96	jfiltsu = j1 - 1 + jfiltsu
97
98	j1 = ifirstloc(rlatv <= 0.)
99
100	call inifilr_hemisph(rlatv(j1 - 1:1:- 1), rlamda, unit, eignfnu, jfiltnv, &
101	matricevn)
102	jfiltnv = j1 - jfiltnv
103	matricevn = matricevn(:, :, jfiltnv:1:- 1)
104
105	call inifilr_hemisph(- rlatv(j1:jjm), rlamda, unit, eignfnu, jfiltsv, &
106	matricevs)
107	jfiltsv = j1 - 1 + jfiltsv
108
109	close(unit)
110	PRINT *, 'jfiltnu =', jfiltnu
111	PRINT *, 'jfiltsu =', jfiltsu
112	PRINT *, 'jfiltnv =', jfiltnv
113	PRINT *, 'jfiltsv =', jfiltsv
114
115	END SUBROUTINE inifilr
116
117	end module inifilr_m