subroutine da_get_gausslats( nj, glats, gwgts, sinlat, coslat) !----------------------------------------------------------------------- ! Purpose: TBD !----------------------------------------------------------------------- implicit none ! Calculates nj Gaussian latitudes i.e. latitudes at which the Legendre ! polynomial Pn(sin(lat)) = 0.0, n=nj, m=0. ! The integral from -1 to +1 of f(x)*Pn(x) where f is a polynomial ! of degree <= 2n-1 can be calculated using ! 0.5 * sum(GaussWgts(:)*Pn(:)*f(:)) with the values at Gaussian latitudes. ! See eqns 77-79 of 'The Spectral Technique' M.Jarraud and A.J.Simmons ! (1983 ECMWF Seminar or 1990 ECMWF Lecture Notes). ! The orthogonality and normalisation of the Legendre polynomials ! checked in this way are very accurate on the Cray, but somewhat ! less accurate on the HPs(32-bit arithmetic). ! Starting with a regular latitude grid, use Newton-Raphson interpolation ! (with bisection steps to add robustness) ! to find the zeros of the Legendre polynomial Pn(x), 0 <= x < 1, ! the negative roots(-1 < x < 0) are set by symmetry. ! ASin(x) gives the Gaussian latitudes. ! This gives slightly better results than finding the roots of Pn(sin(lat)) ! (Algorithm from Numerical Recipies(Fortran version), 1989, p 258) integer, intent(in) :: nj ! Gridpoints in N-S direction. real, intent(out) :: glats(1:nj) ! Gaussian latitudes(S->N, radians). real, intent(out) :: gwgts(1:nj) ! Gaussian weights. real, intent(out), optional :: sinlat(1:nj) ! sin(Latitude). real, intent(out), optional :: coslat(1:nj) ! cos(Latitude). integer, parameter :: maxiter = 100 ! Maximum number of iterations. integer :: i, j, k ! Loop counters. real :: fj, fjold ! Pn(x) on search grid real :: xj, xjold ! search grid real :: x1, x2 ! bounds on root real :: x ! iterated values of x real :: z ! = sqrt(1-x*x) real :: fn ! Pn(x) real :: fn1 ! Pn-1(x) real :: dfr ! 1/Pn'(x) real :: dx, dxold ! step size, previous step if (trace_use) call da_trace_entry("da_get_gausslats") k =(nj + 2) / 2 xj = 0.0 z = 1.0 call da_asslegpol(nj, 0, xj, z, fj) if (mod(nj,2) == 1) then call da_asslegpol(nj-1,0,xj,z,fn1) glats(k) = 0.0 gwgts(k) = 2.0 *(1.0 - xj * xj) /(real(nj) * fn1)**2 k = k+1 end if ! Search interval 0 < x <= 1 for zeros of Legendre polynomials: do j = 2, nj * 2 xjold = xj fjold = fj ! Roots are approximately equally spaced in asin(x) xj = Sin(real(j)*Pi/real(nj*4)) z = sqrt(1.0-xj*xj) call da_asslegpol(nj, 0, xj, z, fj) if (fj >= 0.0 .AND. fjold < 0.0 .OR. fj < 0.0 .AND. fjold >= 0.0) then ! Perform simple interpolation to improve roots(find Gaussian latitudes) if (fjold < 0.0) then ! Orient the search so that fn(x1) < 0 x1 = xjold x2 = xj else x1 = xj x2 = xjold end if x = 0.5*(x1 + x2) ! Initialise the guess for the root dxold = ABS(x1 - x2) ! the step size before last dx = dxold ! and the last step z = sqrt(1.0-x*x) call da_asslegpol(nj, 0, x, z, fn) call da_asslegpol(nj-1,0,x,z,fn1) dfr =(1.0 - x * x) /(real(nj)*(fn1 - x * fn)) do i = 1, maxiter ! Bisect if Newton out of range or not decreasing fast enough if (((x-x1)-fn*dfr)*((x-x2)-fn*dfr) > 0.0 & .OR. ABS(2.0*fn) > ABS(dxold/dfr)) then dxold = dx dx = 0.5 *(x1 - x2) x = x2 + dx else ! Newton-Raphson step dxold = dx dx = fn * dfr x = x - dx end if if (ABS(dx) < 2.0*SPACinG(x)) exit z = sqrt(1.0-x*x) call da_asslegpol(nj,0,x,z,fn) call da_asslegpol(nj-1,0,x,z,fn1) dfr =(1.0 - x * x) /(real(nj)*(fn1 - x * fn)) if (fn < 0.0) then ! Maintain the bracket on the root x1 = x else x2 = x end if end do if (i >= MaxIter) then call da_error(__FILE__,__LINE__, & (/"No convergence finding Gaussian latitudes"/)) end if glats(k) = ASin(x) z = sqrt(1.0-x*x) call da_asslegpol(nj-1,0,x,z,fn1) gwgts(k) = 2.0*(1.0 - x * x) /(real(nj) * fn1)**2 glats(nj+1-k) = -glats(k) gwgts(nj+1-k) = gwgts(k) k=k+1 end if end do if (k /= nj+1) then call da_error(__FILE__,__LINE__,(/"Not all roots found"/)) end if ! Calculate sin, cosine: do j = 1, nj / 2 sinlat(j) = sin(glats(j)) coslat(j) = cos(glats(j)) ! use symmetry for northern hemisphere: sinlat(nj+1-j) = -sinlat(j) coslat(nj+1-j) = coslat(j) end do if ((nj+1) / 2 == nj/2 + 1) then ! Odd, then equator point: glats(nj/2+1) = 0.0 sinlat(nj/2+1) = 0.0 coslat(nj/2+1) = 1.0 end if if (trace_use) call da_trace_exit("da_get_gausslats") end subroutine da_get_gausslats