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 dimensions, 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. - 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	guez	54	module inifilr_m
2	guez	3
3	guez	32	IMPLICIT NONE
4	guez	3
5	guez	164	INTEGER jfiltnu, jfiltnv
6	guez	169	! index of the last line filtered in northern hemisphere at rlat[uv]
7			! latitudes
8	guez	164
9			integer jfiltsu, jfiltsv
10	guez	169	! index of the first line filtered in southern hemisphere at
11			! rlat[uv] latitudes
12	guez	164
13	guez	169	! Filtre pour les champs situes sur la grille scalaire (longitudes
14			! rlonv, latitudes rlatu) :
15			real, pointer:: matriceun(:, :, :) ! (iim, iim, jfiltnu - 1)
16	guez	166	real, pointer:: matriceus(:, :, :) ! (iim, iim, jjm - jfiltsu + 1)
17	guez	136
18	guez	169	! 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	guez	166	real, pointer:: matrinvs(:, :, :) ! (iim, iim, jjm - jfiltsu + 1)
22	guez	165
23	guez	169	! 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	guez	166	real, pointer:: matricevs(:, :, :) ! (iim, iim, jjm - jfiltsv + 1)
27	guez	136
28	guez	54	contains
29	guez	3
30	guez	54	SUBROUTINE inifilr
31	guez	3
32	guez	164	! From filtrez/inifilr.F, version 1.1.1.1, 2004/05/19 12:53:09
33	guez	54	! H. Upadhyaya, O. Sharma
34	guez	3
35	guez	164	! 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	guez	3
42	guez	265	USE dimensions, ONLY : iim, jjm
43	guez	139	USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
44	guez	154	use inifgn_m, only: inifgn
45	guez	166	use inifilr_hemisph_m, only: inifilr_hemisph
46	guez	143	use jumble, only: new_unit
47	guez	169	use nr_util, only: pi, ifirstloc, assert
48	guez	3
49	guez	54	! Local:
50	guez	164
51	guez	132	REAL dlatu(jjm)
52	guez	166	REAL rlamda(2:iim) ! > 0, in descending order
53	guez	168	real eignvl(iim) ! eigenvalues (<= 0) sorted in descending order
54	guez	166	INTEGER j, unit
55	guez	151	REAL colat0 ! > 0
56	guez	169	integer j1 ! index of negative latitude closest to the equator
57	guez	3
58	guez	154	real eignfnu(iim, iim), eignfnv(iim, iim)
59	guez	169	! eigenvectors of the discrete second derivative with respect to
60			! longitude, at rlon[uv] longitudes
61	guez	154
62	guez	54	!-----------------------------------------------------------
63	guez	3
64	guez	54	print *, "Call sequence information: inifilr"
65	guez	3
66	guez	154	CALL inifgn(eignvl, eignfnu, eignfnv)
67	guez	3
68	guez	54	! Calcul de colat0
69	guez	143	forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
70			colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
71	guez	54	PRINT *, 'colat0 = ', colat0
72	guez	3
73	guez	169	rlamda = iim / pi / colat0 * grossismx / sqrt(- eignvl(2: iim))
74	guez	168	print *, "1 / rlamda(iim) = ", 1. / rlamda(iim)
75	guez	169	! This is demonstrated in the notes but just to be sure:
76	guez	175	call assert(rlamda(iim) * colat0 >= 1. - 2. * epsilon(0.), &
77	guez	169	"inifilr rlamda(iim) * colat0")
78
79	guez	166	call new_unit(unit)
80			open(unit, file = "modfrst.csv", status = "replace", action = "write")
81			write(unit, fmt = *) '"rlat (degrees)" modfrst' ! title line
82	guez	3
83	guez	167	j1 = ifirstloc(rlatu <= 0.)
84	guez	3
85	guez	169	call inifilr_hemisph(rlatu(j1 - 1:2:- 1), rlamda, unit, eignfnv, jfiltnu, &
86			matriceun, matrinvn)
87	guez	167	jfiltnu = j1 - jfiltnu
88	guez	166	matriceun = matriceun(:, :, jfiltnu - 1:1:- 1)
89			matrinvn = matrinvn(:, :, jfiltnu - 1:1:- 1)
90	guez	3
91	guez	169	call inifilr_hemisph(- rlatu(j1:jjm), rlamda, unit, eignfnv, jfiltsu, &
92			matriceus, matrinvs)
93	guez	167	jfiltsu = j1 - 1 + jfiltsu
94	guez	3
95	guez	167	j1 = ifirstloc(rlatv <= 0.)
96	guez	3
97	guez	169	call inifilr_hemisph(rlatv(j1 - 1:1:- 1), rlamda, unit, eignfnu, jfiltnv, &
98			matricevn)
99	guez	167	jfiltnv = j1 - jfiltnv
100	guez	166	matricevn = matricevn(:, :, jfiltnv:1:- 1)
101	guez	3
102	guez	169	call inifilr_hemisph(- rlatv(j1:jjm), rlamda, unit, eignfnu, jfiltsv, &
103			matricevs)
104	guez	167	jfiltsv = j1 - 1 + jfiltsv
105	guez	3
106	guez	166	close(unit)
107	guez	151	PRINT *, 'jfiltnu =', jfiltnu
108			PRINT *, 'jfiltsu =', jfiltsu
109			PRINT *, 'jfiltnv =', jfiltnv
110			PRINT *, 'jfiltsv =', jfiltsv
111	guez	3
112	guez	54	END SUBROUTINE inifilr
113	guez	32
114	guez	54	end module inifilr_m