#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  returnend subroutine phcs