var/da/da_spectral/da_get_gausslats.inc

   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