var/da/da_spectral/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    if (trace_use) call da_trace_entry("da_legtra")
  26
  27    index_m = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 - m
  28
  29    if ((nj+1) / 2 == nj/2 + 1) then
  30       eq_coeff = 1.0 ! Odd latitudes
  31    else
  32       eq_coeff = 0.0 ! Even latitudes
  33       eq_term  = 0.0
  34    end if
  35
  36    do l = m, max_wavenumber
  37
  38       sign_switch = (-1)**(l + m)
  39       sum_legtra = da_zero_complex
  40
  41       do j = 1, nj / 2
  42          index_j = (j - 1) * (max_wavenumber + 1) * (max_wavenumber + 2) / 2
  43          index = index_j + index_m + l
  44
  45          ! Sum first quadrant:
  46          sum_legtra = sum_legtra + int_wgts(j) * r_leg(j) * alp(index)
  47
  48          ! Add second quadrant (use symmetry ALP(-mu)=(-1)^{n+|m|}ALP(mu)):
  49          j1 = nj + 1 - j
  50          sum_legtra = sum_legtra + sign_switch * int_wgts(j1) * r_leg(j1) * &
  51             alp(index)
  52       end do
  53
  54       if (eq_coeff > 0.0) then
  55          ! Skip this step for Even lats    ! Syed RH Rizvi
  56          ! Add equator term (wrong if even nj, but then eq_coeff = 0.0 so OK):
  57          j = nj/2 + 1
  58          index_j = (j - 1) * (max_wavenumber + 1) * (max_wavenumber+2) / 2
  59          index = index_j + index_m + l
  60          eq_term = int_wgts(j) * r_leg(j) * alp(index)
  61       end if
  62
  63       v(l) = 0.5 * (sum_legtra + eq_coeff * eq_term)
  64    end do
  65
  66    if (trace_use) call da_trace_exit("da_legtra")
  67
  68 end subroutine da_legtra
  69
  70