var/da/da_spectral/da_setlegpol_test.inc

   1 subroutine da_setlegpol_test (nj, max_wavenumber, alp_size, int_wgts, alp)
   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    real,    intent(in)  :: int_wgts(1:nj)  ! Integration weights.
  13    real,    intent(in)  :: alp(1:alp_size) ! Associated Legendre Polynomials.
  14
  15    real, parameter :: tolerance = 1.0e-6  ! warn if normalization error exceeds
  16
  17    integer :: m, l1, l2, j, j1    ! Loop counters.
  18    integer :: index_m, index_j    ! Markers.
  19    integer :: index1, index2      ! Markers.
  20    integer :: sign_switch1        ! Defined to make use of symmetry of ALPs.
  21    integer :: sign_switch2        ! Defined to make use of symmetry of ALPs.
  22    real    :: eq_coeff            ! 1 if equator point, 0 otherwise.
  23    real    :: alp_norm_test       ! Summation scalar.
  24    real    :: eq_term             ! Summation scalar.
  25    integer :: spec_unit
  26
  27    if (trace_use) call da_trace_entry("da_setlegpol_test")
  28
  29    call da_get_unit(spec_unit)
  30    open(unit=spec_unit,file="spec_pol",status="replace")
  31
  32    if ((nj+1) / 2 == nj/2 + 1) then
  33       eq_coeff = 1.0 ! Odd latitudes
  34    else
  35       eq_coeff = 0.0 ! Even latitudes
  36       eq_term  = 0.0
  37    end if
  38
  39    ! Test 0.5 * integral_-1^1 alp(j,l1,m) * alp(j,l2,m) = 1 if l1=l2,
  40    ! 0 otherwise:
  41
  42    do m = 0, max_wavenumber
  43       index_m = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 - m
  44       do l1 = m, max_wavenumber
  45          do l2 = m, max_wavenumber
  46
  47             sign_switch1 = (-1)**(l1 + m)
  48             sign_switch2 = (-1)**(l2 + m)
  49
  50             alp_norm_test = 0.0
  51             do j = 1, nj / 2
  52                index_j = (j - 1) * (max_wavenumber+1) * (max_wavenumber+2) /2
  53                index1 = index_j + index_m + l1
  54                index2 = index_j + index_m + l2
  55
  56                ! Sum first quadrant:
  57                alp_norm_test = alp_norm_test + int_wgts(j) * alp(index1) &
  58                   * alp(index2)
  59
  60                ! Add second quadrant (use symmetry ALP(-mu)=(-1)^{n+|m|}ALP(mu)):
  61                j1 = nj + 1 - j
  62                alp_norm_test = alp_norm_test + int_wgts(j1) * &
  63                   sign_switch1 * alp(index1) * sign_switch2 * alp(index2)
  64             end do
  65
  66             if (eq_coeff > 0.0) then
  67                ! Skip this step for even lats       R! Syed RH Rizvi! S
  68                ! Add equator term (wrong if even nj, but then eq_coeff = 0.0
  69                ! so OK):
  70                j = nj/2 + 1
  71                index_j = (j - 1) * (max_wavenumber+1) * (max_wavenumber+2) /2
  72                index1 = index_j + index_m + l1
  73                index2 = index_j + index_m + l2
  74
  75                eq_term = int_wgts(j) * alp(index1) * alp(index2)
  76             end if
  77             alp_norm_test = 0.5 * (alp_norm_test + eq_coeff * eq_term)
  78
  79             ! if (l1 /= l2 .and. abs(alp_norm_test) >= tolerance) then
  80             !    write(unit=stdout,fmt='(a,3i6,f15.10,a,f15.10)')
  81             !      ' warning: ALP normalization error (m, l1, l2) = ', !&
  82             !                                      m, l1, l2, alp_norm_test, &
  83             !                                      ', > tolerance = ', tolerance
  84             !            end if
  85             if (l1 == l2 .and. abs(alp_norm_test-1.0) >= tolerance) then
  86                write(unit=spec_unit,fmt='(a,3i6,f15.10,a,f15.10)') &
  87                  ' warning: ALP normalization error (m, l1, l2) = ', &
  88                  m, l1, l2, alp_norm_test - 1.0, &
  89                  ', > tolerance = ', tolerance
  90
  91             end if
  92          end do
  93       end do
  94    end do
  95
  96    close(spec_unit)
  97    call da_free_unit(spec_unit)
  98
  99    if (trace_use) call da_trace_exit("da_setlegpol_test")
 100
 101 end subroutine da_setlegpol_test
 102
 103