phcs.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
#include <params.h>

subroutine phcs(pmn     ,hmn     ,ix      ,x1)

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Compute associated Legendre functions of the first kind of order m and
! degree n, and the associated derivatives for arg x1.

! Method: 
! Compute associated Legendre functions of the first kind of order m and
! degree n, and the associated derivatives for arg x1.  The associated
! Legendre functions are evaluated using relationships contained in
! "Tables of Normalized Associated Legendre Polynomials",
! S. L. Belousov (1962).  Both the functions and their derivatives are
! ordered in a linear stored rectangular array (with a large enough
! domain to contain the particular wavenumber truncation defined in the
! pspect common block) by column. m = 0->ptrm, and  n = m->ptrn + m 
!                m
! The functions P (x) are normalized such that 
!                n
!                          /   m     2
!                          | [P  (x)] dx = 1/2
!                          /   n
!                             __
! and must be multiplied by  |2  to match Belousov tables.
!                           \|
!                  m
! The derivatives H (x) are defined as 
!                  n        m           2    m
!                          H (x) = -(1-x ) dP (x)/dx
!                           n                n
!
! and are evaluated using the recurrence relationship
!                          _________________________
!      m          m       |  2   2                     m
!     H (x) = nx P (x) -  |(n - m )(2n + 1)/(2n - 1)  P   (x)
!      n          n      \|                            n-1
!
! Modified 1/23/97 by Jim Rosinski to use real*16 arithmetic in order to 
! achieve (nearly) identical values on all machines.
! 
! Author: CCM1
! 
!-----------------------------------------------------------------------
!
! $Id: phcs.F90,v 1.1 2001/11/06 18:42:49 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

  use precision
  use pmgrid
  use pspect
  implicit none

!------------------------------Arguments--------------------------------
  integer , intent(in)  :: ix       ! Dimension of Legendre funct arrays
  real(r8), intent(in)  :: x1       ! sin of latitude, [sin(phi), or mu]
  real(r8), intent(out) ::  pmn(ix) ! Legendre function array
  real(r8), intent(out) ::  hmn(ix) ! Derivative array
!-----------------------------------------------------------------------

!---------------------------Local variables-----------------------------
#ifdef PGF90
  integer, parameter :: r16 = selected_real_kind(12)
#else
  integer, parameter :: r16 = selected_real_kind(17)
#endif
  integer jmax       ! Loop limit (N+1=> 2D wavenumber limit +1)
  integer nmax       ! Large enough n to envelope truncation
  integer n          ! 2-D wavenumber index (up/down column)
  integer ml         ! intermediate scratch variable
  integer k          ! counter on terms in trig series expansion
  integer n2         ! 2*n
  integer m          ! zonal wavenumber index
  integer nto        ! intermediate scratch variable
  integer mto        ! intermediate scratch variable
  integer j          ! 2-D wavenumber index in recurrence evaluation
  integer nmaxm      ! loop limit in recurrence evaluation

  real(r16) xtemp(3,pmmax+ptrn+1)  ! Workspace for evaluating recurrence
!                                    ! relation where xtemp(m-2,n) and
!                                    ! xtemp(m-1,n) contain Pmn's required 
!                                    ! to evaluate xtemp(m,n) (i.e.,always
!                                    ! contains three adjacent columns of
!                                    ! the Pmn data structure)
!
  real(r16) xx1        ! x1 in extended precision
  real(r16) xte        ! cosine latitude [cos(phi)]
  real(r16) teta       ! pi/2 - latitute (colatitude)
  real(r16) an         ! coefficient on trig. series expansion
  real(r16) sinpar     ! accumulator in trig. series expansion 
  real(r16) cospar     ! accumulator in trig. series expansion 
  real(r16) p          ! 2-D wavenumber (series expansion)
  real(r16) q          ! intermediate variable in series expansion 
  real(r16) r          ! zonal wavenumber (recurrence evaluation)
  real(r16) p2         ! intermediate variable in series expansion 
  real(r16) rr         ! twice the zonal wavenumber (recurrence)
  real(r16) sqp        ! intermediate variable in series expansion 
  real(r16) cosfak     ! coef. on cos term in series expansion
  real(r16) sinfak     ! coef. on sin term in series expansion
  real(r16) ateta      ! intermediate variable in series expansion 
  real(r16) costet     ! cos term in trigonometric series expansion
  real(r16) sintet     ! sin term in trigonometric series expansion
!
  real(r16) t          ! intermediate variable (recurrence evaluation)
  real(r16) wm2        ! intermediate variable (recurrence evaluation)
  real(r16) wmq2       ! intermediate variable (recurrence evaluation)
  real(r16) w          ! intermediate variable (recurrence evaluation)
  real(r16) wq         ! intermediate variable (recurrence evaluation)
  real(r16) q2         ! intermediate variable (recurrence evaluation)
  real(r16) wt         ! intermediate variable (recurrence evaluation)
  real(r16) q2d        ! intermediate variable (recurrence evaluation)
  real(r16) cmn        ! cmn  recurrence coefficient (see Belousov)
  real(r16) xdmn       ! dmn  recurrence coefficient (see Belousov)
  real(r16) emn        ! emn  recurrence coefficient (see Belousov)
  real(r16) n2m1       ! n2 - 1 in extended precision
  real(r16) n2m3       ! n2 - 3 in extended precision
  real(r16) n2p1nnm1   ! (n2+1)*(n*n-1) in extended precision
  real(r16) twopmq     ! p + p - q in extended precision
!-----------------------------------------------------------------------
!
! Begin procedure by evaluating the first two columns of the Legendre
! function matrix (i.e., all n for m=0,1) via a trigonometric series 
! expansion (see eqs. 19 and 21 in Belousov, 1962).  Note that indexing
! is offset by one (e.g., m index for wavenumber m=0 is 1 and so on)
! Setup first ...
!
  xx1  = x1
  jmax = ptrn + 1
  nmax = pmmax + jmax
  xte = sqrt(1.-xx1*xx1)
  teta = acos(xx1)
  an = 1.
  xtemp(1,1) = 0.5    ! P00
!
! begin loop over n (2D wavenumber, or degree of associated Legendre 
! function) beginning with n=1 (i.e., P00 was assigned above)
! note n odd/even distinction yielding 2 results per n cycle
!
  do n=2,nmax
     sinpar = 0.
     cospar = 0.
     ml = n
     p = n - 1
     p2 = p*p
     sqp = 1./sqrt(p2+p)
     an = an*sqrt(1. - 1./(4.*p2))
     cosfak = 1.
     sinfak = p*sqp
     do k=1,ml,2
        q = k - 1
        twopmq = p + p - q
        ateta = (p-q)*teta
        costet = cos(ateta)
        sintet = sin(ateta)
        if (n==k) costet = costet*0.5
        if (k/=1) then
           cosfak = (q-1.)/q*(twopmq+2.)/(twopmq+1.)*cosfak
           sinfak = cosfak*(p-q)*sqp
        end if
        cospar = cospar + costet*cosfak
        sinpar = sinpar + sintet*sinfak
     end do
     xtemp(1,n)   = an*cospar      ! P0n vector
     xtemp(2,n-1) = an*sinpar      ! P1n vector
  end do
!
! Assign Legendre functions and evaluate derivatives for all n and m=0,1
!
  pmn(1) = 0.5
  pmn(1+jmax) = xtemp(2,1)
  hmn(1) = 0.
  hmn(1+jmax) = xx1*xtemp(2,1)
  do n=2,jmax
     pmn(n) = xtemp(1,n)
     pmn(n+jmax) = xtemp(2,n)
     n2 = n + n
     n2m1 = n2 - 1
     n2m3 = n2 - 3
     n2p1nnm1 = (n2+1)*(n*n-1)
     hmn(n) = (n-1)*(xx1*xtemp(1,n)-sqrt(n2m1/n2m3)*xtemp(1,n-1))
     hmn(n+jmax) = n*xx1*xtemp(2,n)-sqrt(n2p1nnm1/n2m1)*xtemp(2,n-1)
  end do
!
! Evaluate recurrence relationship for remaining Legendre functions
! (i.e., m=2 ... PTRM) and associated derivatives (see eq 17, Belousov)
!
  do m=3,pmmax
     r = m - 1
     rr = r + r
     xtemp(3,1) = sqrt(1.+1./rr)*xte*xtemp(2,1)
     nto = (m-1)*jmax
     pmn(nto+1) = xtemp(3,1)
     hmn(nto+1) = r*xx1*xtemp(3,1)
     nmaxm = nmax - m
!
! Loop over 2-D wavenumber (i.e., degree of Legendre function)
! Pmn's and Hmn's for current zonal wavenumber, r
!
     do j=2,nmaxm
        mto = nto + j
        t = j - 1
        q = rr + t - 1
        wm2 = q + t
        w = wm2 + 2
        wq = w*q
        q2 = q*q - 1
        wmq2 = wm2*q2
        wt = w*t
        q2d = q2 + q2
        cmn = sqrt((wq*(q-2.))/(wmq2-q2d))
        xdmn = sqrt((wq*(t+1.))/wmq2)
        emn = sqrt(wt/((q+1.)*wm2))
        xtemp(3,j) = cmn*xtemp(1,j) - xx1*(xdmn*xtemp(1,j+1)-emn*xtemp(3,j-1))
        pmn(mto) = xtemp(3,j)
        hmn(mto) = (r+t)*xx1*xtemp(3,j) - sqrt(wt*(q+1.)/wm2)*xtemp(3,j-1)
     end do
!
! shift Pmn's to left in workspace (setup for next recurrence pass)
!
!++pjr
! not initialized above
     xtemp(2,nmax) = 0.
     do j=nmaxm,nmax
        xtemp(3,j) = 0.
     end do
!--pjr
     do n=1,nmax
        xtemp(1,n) = xtemp(2,n)
        xtemp(2,n) = xtemp(3,n)
     end do
  end do

  return
end subroutine phcs