da_legtra.inc

1 subroutine da_legtra (nj, max_wavenumber, alp_size, m, int_wgts, alp, r_leg, v)
2 
3    !-----------------------------------------------------------------------
4    ! Purpose: TBD
5    !-----------------------------------------------------------------------
6 
7    implicit none
8 
9    integer, intent(in)  :: nj                      ! Number of latitudes.
10    integer, intent(in)  :: max_wavenumber          ! Maximum wavenumber.
11    integer, intent(in)  :: alp_size                ! Dimension of ALPs.
12    integer, intent(in)  :: m                       ! Zonal wavenumber.
13    real,    intent(in)  :: int_wgts(1:nj)          ! Integration weights.
14    real,    intent(in)  :: alp(1:alp_size)         ! Associated Legendre Polynomials.
15    complex, intent(in)  :: r_leg(1:nj)             ! Field to transform.
16    complex, intent(out) :: v(m:max_wavenumber)     ! Output spectral coefficient.
17 
18    integer              :: l, j, j1                ! Loop counters.
19    integer              :: index_m, index_j, index ! Markers.
20    integer              :: sign_switch             ! make use of symmetry of ALPs.
21    real                 :: eq_coeff                ! 1 if equator point, 0 otherwise.
22    complex              :: sum_legtra              ! Summation scalar.
23    complex              :: eq_term                 ! Summation scalar.
24 
25    index_m = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 - m
26 
27    if ((nj+1) / 2 == nj/2 + 1) then
28       eq_coeff = 1.0 ! Odd latitudes
29    else
30       eq_coeff = 0.0 ! Even latitudes
31       eq_term  = 0.0
32    end if
33 
34    do l = m, max_wavenumber
35 
36       sign_switch = (-1)**(l + m)
37       sum_legtra = da_zero_complex
38 
39       do j = 1, nj / 2
40          index_j = (j - 1) * (max_wavenumber + 1) * (max_wavenumber + 2) / 2
41          index = index_j + index_m + l
42 
43          ! Sum first quadrant:
44          sum_legtra = sum_legtra + int_wgts(j) * r_leg(j) * alp(index)
45 
46          ! Add second quadrant (use symmetry ALP(-mu)=(-1)^{n+|m|}ALP(mu)):
47          j1 = nj + 1 - j
48          sum_legtra = sum_legtra + sign_switch * int_wgts(j1) * r_leg(j1) * &
49             alp(index)
50       end do
51      
52       if (eq_coeff > 0.0) then
53          ! Skip this step for Even lats    ! Syed RH Rizvi
54          ! Add equator term (wrong if even nj, but then eq_coeff = 0.0 so OK):
55          j = nj/2 + 1
56          index_j = (j - 1) * (max_wavenumber + 1) * (max_wavenumber+2) / 2
57          index = index_j + index_m + l
58          eq_term = int_wgts(j) * r_leg(j) * alp(index)
59       end if
60 
61       v(l) = 0.5 * (sum_legtra + eq_coeff * eq_term)
62    end do
63 
64 end subroutine da_legtra
65 
66