da_get_gausslats.inc
References to this file elsewhere.
1 subroutine da_get_gausslats( nj, glats, gwgts, sinlat, coslat)
2
3 !-----------------------------------------------------------------------
4 ! Purpose: TBD
5 !-----------------------------------------------------------------------
6
7 implicit none
8
9 ! Calculates nj Gaussian latitudes i.e. latitudes at which the Legendre
10 ! polynomial Pn(sin(lat)) = 0.0, n=nj, m=0.
11 ! The integral from -1 to +1 of f(x)*Pn(x) where f is a polynomial
12 ! of degree <= 2n-1 can be calculated using
13 ! 0.5 * sum(GaussWgts(:)*Pn(:)*f(:)) with the values at Gaussian latitudes.
14 ! See eqns 77-79 of 'The Spectral Technique' M.Jarraud and A.J.Simmons
15 ! (1983 ECMWF Seminar or 1990 ECMWF Lecture Notes).
16 ! The orthogonality and normalisation of the Legendre polynomials
17 ! checked in this way are very accurate on the Cray, but somewhat
18 ! less accurate on the HPs(32-bit arithmetic).
19 ! Starting with a regular latitude grid, use Newton-Raphson interpolation
20 ! (with bisection steps to add robustness)
21 ! to find the zeros of the Legendre polynomial Pn(x), 0 <= x < 1,
22 ! the negative roots(-1 < x < 0) are set by symmetry.
23 ! ASin(x) gives the Gaussian latitudes.
24 ! This gives slightly better results than finding the roots of Pn(sin(lat))
25 ! (Algorithm from Numerical Recipies(Fortran version), 1989, p 258)
26
27 integer, intent(in) :: nj ! Gridpoints in N-S direction.
28 real, intent(out) :: glats(1:nj) ! Gaussian latitudes(S->N, radians).
29 real, intent(out) :: gwgts(1:nj) ! Gaussian weights.
30 real, intent(out), optional :: sinlat(1:nj) ! sin(Latitude).
31 real, intent(out), optional :: coslat(1:nj) ! cos(Latitude).
32
33 integer, parameter :: maxiter = 100 ! Maximum number of iterations.
34 integer :: i, j, k ! Loop counters.
35
36 real :: fj, fjold ! Pn(x) on search grid
37 real :: xj, xjold ! search grid
38 real :: x1, x2 ! bounds on root
39 real :: x ! iterated values of x
40 real :: z ! = sqrt(1-x*x)
41 real :: fn ! Pn(x)
42 real :: fn1 ! Pn-1(x)
43 real :: dfr ! 1/Pn'(x)
44 real :: dx, dxold ! step size, previous step
45
46 if (trace_use) call da_trace_entry("da_get_gausslats")
47
48 k =(nj + 2) / 2
49 xj = 0.0
50 z = 1.0
51
52 call da_asslegpol(nj, 0, xj, z, fj)
53
54 if (mod(nj,2) == 1) then
55 call da_asslegpol(nj-1,0,xj,z,fn1)
56 glats(k) = 0.0
57 gwgts(k) = 2.0 *(1.0 - xj * xj) /(real(nj) * fn1)**2
58 k = k+1
59 end if
60
61 ! Search interval 0 < x <= 1 for zeros of Legendre polynomials:
62 do j = 2, nj * 2
63 xjold = xj
64 fjold = fj
65
66 ! Roots are approximately equally spaced in asin(x)
67 xj = Sin(real(j)*Pi/real(nj*4))
68 z = sqrt(1.0-xj*xj)
69 call da_asslegpol(nj, 0, xj, z, fj)
70
71 if (fj >= 0.0 .AND. fjold < 0.0 .OR. fj < 0.0 .AND. fjold >= 0.0) then
72
73 ! Perform simple interpolation to improve roots(find Gaussian latitudes)
74 if (fjold < 0.0) then ! Orient the search so that fn(x1) < 0
75 x1 = xjold
76 x2 = xj
77 else
78 x1 = xj
79 x2 = xjold
80 end if
81
82 x = 0.5*(x1 + x2) ! Initialise the guess for the root
83 dxold = ABS(x1 - x2) ! the step size before last
84 dx = dxold ! and the last step
85 z = sqrt(1.0-x*x)
86 call da_asslegpol(nj, 0, x, z, fn)
87 call da_asslegpol(nj-1,0,x,z,fn1)
88 dfr =(1.0 - x * x) /(real(nj)*(fn1 - x * fn))
89
90 do i = 1, maxiter
91
92 ! Bisect if Newton out of range or not decreasing fast enough
93 if (((x-x1)-fn*dfr)*((x-x2)-fn*dfr) > 0.0 &
94 .OR. ABS(2.0*fn) > ABS(dxold/dfr)) then
95 dxold = dx
96 dx = 0.5 *(x1 - x2)
97 x = x2 + dx
98 else ! Newton-Raphson step
99 dxold = dx
100 dx = fn * dfr
101 x = x - dx
102 end if
103
104 if (ABS(dx) < 2.0*SPACinG(x)) exit
105 z = sqrt(1.0-x*x)
106 call da_asslegpol(nj,0,x,z,fn)
107 call da_asslegpol(nj-1,0,x,z,fn1)
108 dfr =(1.0 - x * x) /(real(nj)*(fn1 - x * fn))
109
110 if (fn < 0.0) then ! Maintain the bracket on the root
111 x1 = x
112 else
113 x2 = x
114 end if
115 end do
116
117 if (i >= MaxIter) then
118 call da_error(__FILE__,__LINE__, &
119 (/"No convergence finding Gaussian latitudes"/))
120 end if
121
122 glats(k) = ASin(x)
123 z = sqrt(1.0-x*x)
124 call da_asslegpol(nj-1,0,x,z,fn1)
125 gwgts(k) = 2.0*(1.0 - x * x) /(real(nj) * fn1)**2
126 glats(nj+1-k) = -glats(k)
127 gwgts(nj+1-k) = gwgts(k)
128 k=k+1
129 end if
130 end do
131
132 if (k /= nj+1) then
133 call da_error(__FILE__,__LINE__,(/"Not all roots found"/))
134 end if
135
136 ! Calculate sin, cosine:
137
138 do j = 1, nj / 2
139 sinlat(j) = sin(glats(j))
140 coslat(j) = cos(glats(j))
141
142 ! use symmetry for northern hemisphere:
143 sinlat(nj+1-j) = -sinlat(j)
144 coslat(nj+1-j) = coslat(j)
145 end do
146
147 if ((nj+1) / 2 == nj/2 + 1) then ! Odd, then equator point:
148 glats(nj/2+1) = 0.0
149 sinlat(nj/2+1) = 0.0
150 coslat(nj/2+1) = 1.0
151 end if
152
153 if (trace_use) call da_trace_exit("da_get_gausslats")
154
155 end subroutine da_get_gausslats
156
157