<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
<A NAME='DA_LEGTRA'><A href='../../html_code/spectral/da_legtra.inc.html#DA_LEGTRA' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_legtra (nj, max_wavenumber, alp_size, m, int_wgts, alp, r_leg, v) 2,2
!-----------------------------------------------------------------------
! Purpose: TBD
!-----------------------------------------------------------------------
implicit none
integer, intent(in) :: nj ! Number of latitudes.
integer, intent(in) :: max_wavenumber ! Maximum wavenumber.
integer, intent(in) :: alp_size ! Dimension of ALPs.
integer, intent(in) :: m ! Zonal wavenumber.
real, intent(in) :: int_wgts(1:nj) ! Integration weights.
real, intent(in) :: alp(1:alp_size) ! Associated Legendre Polynomials.
complex, intent(in) :: r_leg(1:nj) ! Field to transform.
complex, intent(out) :: v(m:max_wavenumber) ! Output spectral coefficient.
integer :: l, j, j1 ! Loop counters.
integer :: index_m, index_j, index ! Markers.
integer :: sign_switch ! make use of symmetry of ALPs.
real :: eq_coeff ! 1 if equator point, 0 otherwise.
complex :: sum_legtra ! Summation scalar.
complex :: eq_term ! Summation scalar.
if (trace_use) call da_trace_entry
("da_legtra")
index_m = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 - m
if ((nj+1) / 2 == nj/2 + 1) then
eq_coeff = 1.0 ! Odd latitudes
else
eq_coeff = 0.0 ! Even latitudes
eq_term = 0.0
end if
do l = m, max_wavenumber
sign_switch = (-1)**(l + m)
sum_legtra = da_zero_complex
do j = 1, nj / 2
index_j = (j - 1) * (max_wavenumber + 1) * (max_wavenumber + 2) / 2
index = index_j + index_m + l
! Sum first quadrant:
sum_legtra = sum_legtra + int_wgts(j) * r_leg(j) * alp(index)
! Add second quadrant (use symmetry ALP(-mu)=(-1)^{n+|m|}ALP(mu)):
j1 = nj + 1 - j
sum_legtra = sum_legtra + sign_switch * int_wgts(j1) * r_leg(j1) * &
alp(index)
end do
if (eq_coeff > 0.0) then
! Skip this step for Even lats ! Syed RH Rizvi
! Add equator term (wrong if even nj, but then eq_coeff = 0.0 so OK):
j = nj/2 + 1
index_j = (j - 1) * (max_wavenumber + 1) * (max_wavenumber+2) / 2
index = index_j + index_m + l
eq_term = int_wgts(j) * r_leg(j) * alp(index)
end if
v(l) = 0.5 * (sum_legtra + eq_coeff * eq_term)
end do
if (trace_use) call da_trace_exit("da_legtra")
end subroutine da_legtra