da_setlegpol_test.inc
References to this file elsewhere.
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 call da_get_unit(spec_unit)
28 open(unit=spec_unit,file="spec_pol",status="replace")
29
30 if ((nj+1) / 2 == nj/2 + 1) then
31 eq_coeff = 1.0 ! Odd latitudes
32 else
33 eq_coeff = 0.0 ! Even latitudes
34 eq_term = 0.0
35 end if
36
37 ! Test 0.5 * integral_-1^1 alp(j,l1,m) * alp(j,l2,m) = 1 if l1=l2,
38 ! 0 otherwise:
39
40 do m = 0, max_wavenumber
41 index_m = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 - m
42 do l1 = m, max_wavenumber
43 do l2 = m, max_wavenumber
44
45 sign_switch1 = (-1)**(l1 + m)
46 sign_switch2 = (-1)**(l2 + m)
47
48 alp_norm_test = 0.0
49 do j = 1, nj / 2
50 index_j = (j - 1) * (max_wavenumber+1) * (max_wavenumber+2) /2
51 index1 = index_j + index_m + l1
52 index2 = index_j + index_m + l2
53
54 ! Sum first quadrant:
55 alp_norm_test = alp_norm_test + int_wgts(j) * alp(index1) &
56 * alp(index2)
57
58 ! Add second quadrant (use symmetry ALP(-mu)=(-1)^{n+|m|}ALP(mu)):
59 j1 = nj + 1 - j
60 alp_norm_test = alp_norm_test + int_wgts(j1) * &
61 sign_switch1 * alp(index1) * sign_switch2 * alp(index2)
62 end do
63
64 if (eq_coeff > 0.0) then
65 ! Skip this step for even lats R! Syed RH Rizvi! S
66 ! Add equator term (wrong if even nj, but then eq_coeff = 0.0
67 ! so OK):
68 j = nj/2 + 1
69 index_j = (j - 1) * (max_wavenumber+1) * (max_wavenumber+2) /2
70 index1 = index_j + index_m + l1
71 index2 = index_j + index_m + l2
72
73 eq_term = int_wgts(j) * alp(index1) * alp(index2)
74 end if
75 alp_norm_test = 0.5 * (alp_norm_test + eq_coeff * eq_term)
76
77 ! if (l1 /= l2 .and. abs(alp_norm_test) >= tolerance) then
78 ! write(unit=stdout,fmt='(a,3i6,f15.10,a,f15.10)')
79 ! ' warning: ALP normalization error (m, l1, l2) = ', !&
80 ! m, l1, l2, alp_norm_test, &
81 ! ', > tolerance = ', tolerance
82 ! end if
83 if (l1 == l2 .and. abs(alp_norm_test-1.0) >= tolerance) then
84 write(unit=spec_unit,fmt='(a,3i6,f15.10,a,f15.10)') &
85 ' warning: ALP normalization error (m, l1, l2) = ', &
86 m, l1, l2, alp_norm_test - 1.0, &
87 ', > tolerance = ', tolerance
88
89 end if
90 end do
91 end do
92 end do
93
94 close(spec_unit)
95 call da_free_unit(spec_unit)
96
97 end subroutine da_setlegpol_test
98
99