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 k =(nj + 2) / 2
47 xj = 0.0
48 z = 1.0
49
50 call da_asslegpol(nj, 0, xj, z, fj)
51
52 if (mod(nj,2) == 1) then
53 call da_asslegpol(nj-1,0,xj,z,fn1)
54 glats(k) = 0.0
55 gwgts(k) = 2.0 *(1.0 - xj * xj) /(real(nj) * fn1)**2
56 k = k+1
57 end if
58
59 ! Search interval 0 < x <= 1 for zeros of Legendre polynomials:
60 do j = 2, nj * 2
61 xjold = xj
62 fjold = fj
63
64 ! Roots are approximately equally spaced in asin(x)
65 xj = Sin(real(j)*Pi/real(nj*4))
66 z = sqrt(1.0-xj*xj)
67 call da_asslegpol(nj, 0, xj, z, fj)
68
69 if (fj >= 0.0 .AND. fjold < 0.0 .OR. fj < 0.0 .AND. fjold >= 0.0) then
70
71 ! Perform simple interpolation to improve roots(find Gaussian latitudes)
72 if (fjold < 0.0) then ! Orient the search so that fn(x1) < 0
73 x1 = xjold
74 x2 = xj
75 else
76 x1 = xj
77 x2 = xjold
78 end if
79
80 x = 0.5*(x1 + x2) ! Initialise the guess for the root
81 dxold = ABS(x1 - x2) ! the step size before last
82 dx = dxold ! and the last step
83 z = sqrt(1.0-x*x)
84 call da_asslegpol(nj, 0, x, z, fn)
85 call da_asslegpol(nj-1,0,x,z,fn1)
86 dfr =(1.0 - x * x) /(real(nj)*(fn1 - x * fn))
87
88 do i = 1, maxiter
89
90 ! Bisect if Newton out of range or not decreasing fast enough
91 if (((x-x1)-fn*dfr)*((x-x2)-fn*dfr) > 0.0 &
92 .OR. ABS(2.0*fn) > ABS(dxold/dfr)) then
93 dxold = dx
94 dx = 0.5 *(x1 - x2)
95 x = x2 + dx
96 else ! Newton-Raphson step
97 dxold = dx
98 dx = fn * dfr
99 x = x - dx
100 end if
101
102 if (ABS(dx) < 2.0*SPACinG(x)) exit
103 z = sqrt(1.0-x*x)
104 call da_asslegpol(nj,0,x,z,fn)
105 call da_asslegpol(nj-1,0,x,z,fn1)
106 dfr =(1.0 - x * x) /(real(nj)*(fn1 - x * fn))
107
108 if (fn < 0.0) then ! Maintain the bracket on the root
109 x1 = x
110 else
111 x2 = x
112 end if
113 end do
114
115 if (i >= MaxIter) then
116 call da_error(__FILE__,__LINE__, &
117 (/"No convergence finding Gaussian latitudes"/))
118 end if
119
120 glats(k) = ASin(x)
121 z = sqrt(1.0-x*x)
122 call da_asslegpol(nj-1,0,x,z,fn1)
123 gwgts(k) = 2.0*(1.0 - x * x) /(real(nj) * fn1)**2
124 glats(nj+1-k) = -glats(k)
125 gwgts(nj+1-k) = gwgts(k)
126 k=k+1
127 end if
128 end do
129
130 if (k /= nj+1) then
131 call da_error(__FILE__,__LINE__,(/"Not all roots found"/))
132 end if
133
134 ! Calculate sin, cosine:
135
136 do j = 1, nj / 2
137 sinlat(j) = sin(glats(j))
138 coslat(j) = cos(glats(j))
139
140 ! use symmetry for northern hemisphere:
141 sinlat(nj+1-j) = -sinlat(j)
142 coslat(nj+1-j) = coslat(j)
143 end do
144
145 if ((nj+1) / 2 == nj/2 + 1) then ! Odd, then equator point:
146 glats(nj/2+1) = 0.0
147 sinlat(nj/2+1) = 0.0
148 coslat(nj/2+1) = 1.0
149 end if
150
151 end subroutine da_get_gausslats
152
153