da_setlegpol.inc
References to this file elsewhere.
1 subroutine da_setlegpol(nj, max_wavenumber, alp_size, sinlat, coslat, alp)
2
3 !-----------------------------------------------------------------------
4 ! Purpose: TBD
5 !-----------------------------------------------------------------------
6
7 implicit none
8
9 ! Method:
10 ! Uses ECMWF recursion relation as opposed to Num Rec. one which can only go
11 ! to m about 12 with single precision). However, still use NumRec one for
12 ! m = 0 and 1 as ECMWF one requires m-2 for asslegpol of order m.
13 ! Reference: Jarraud and Simmons (1990) and Belousov (1962).
14
15 integer, intent(in) :: nj ! Latitude dimension.
16 integer, intent(in) :: max_wavenumber ! Maximum zonal wavenumber.
17 integer, intent(in) :: alp_size ! Ass. Leg. Pol. array size.
18 real, intent(in) :: sinlat(1:nj) ! sine(latitude).
19 real, intent(in) :: coslat(1:nj) ! cosine(latitude).
20 real, intent(out) :: alp(1:alp_size) ! Associated Legendre Polynomial.
21
22 integer :: j, l, m, mm ! Loop counters.
23 integer :: l1, l2, m2 ! Markers.
24 integer :: index_j, index_m, index ! Indexing counters for ALP array.
25 integer :: index_m2, index_l2m2 ! Indexing counters for ALP array.
26 integer :: index_l1m2, index_l1m ! Indexing counters for ALP array.
27 real :: s
28 real :: c(2:max_wavenumber,2:max_wavenumber) ! Recursion coefficient.
29 real :: d(2:max_wavenumber,2:max_wavenumber) ! Recursion coefficient.
30 real :: e(2:max_wavenumber,2:max_wavenumber) ! Recursion coefficient.
31
32 alp(:) = 0.0
33
34 ! Compute Associated Legendre polynomials for latitude range:
35
36 do j = 1, (nj + 1) / 2
37 index_j = (j - 1) * (max_wavenumber + 1) * (max_wavenumber + 2) / 2
38
39 ! use Num. Rec. recursion relation for m = 0, 1:
40
41 do m = 0, 1
42 index_m = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 - m
43 do l = m, max_wavenumber
44 index = index_m + index_j + l
45 call da_asslegpol(l, m, sinlat(j), coslat(j), alp(index))
46
47 ! Multiply by normalization constant
48 ! (to ensure 1/2 integral^1_-1 Pml Pml1 = 1 if l = l1):
49
50 s = 1.0
51 do mm = l-m+1, l+m
52 s = s * real(mm)
53 end do
54 alp(index) = sqrt(real(2*l+1) / s) * alp(index)
55 end do
56 end do
57 end do
58
59 ! Jarraud recursion relation coefficients:
60
61 do m = 2, max_wavenumber
62 do l = m, max_wavenumber
63 c(l,m) = sqrt ((real(2*l+1)/real(2*l-3)) * (real(l+m-1)/real(l+m)) * &
64 (real(l+m-3)/real(l+m-2)))
65 d(l,m) = sqrt ((real(2*l+1)/real(2*l-1)) * (real(l+m-1)/real(l+m)) * &
66 (real(l-m+1)/real(l+m-2)))
67 e(l,m) = sqrt ((real(2*l+1)/real(2*l-1)) * (real(l-m) /real(l+m)))
68 end do
69 end do
70
71 ! use Jarraud recursion relation for m>=2:
72
73 do j = 1, (nj + 1) / 2
74 index_j = (j - 1) * (max_wavenumber + 1) * (max_wavenumber + 2) / 2
75
76 do m = 2, max_wavenumber
77 index_m = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 - m
78 m2 = m - 2
79 index_m2 = m2 * (max_wavenumber + 1 - m2) + m2 * (m2+1) / 2 + 1 - m2
80
81 do l = m, max_wavenumber
82 l1 = l - 1
83 l2 = l - 2
84 index = index_j + index_m + l
85 index_l2m2 = index_j + index_m2 + l2
86 index_l1m2 = index_j + index_m2 + l1
87 index_l1m = index_j + index_m + l1
88
89 alp(index) = c(l,m) * alp(index_l2m2) - d(l,m) * sinlat(j) * &
90 alp(index_l1m2) + e(l,m) * sinlat(j) * alp(index_l1m)
91 end do
92 end do
93 end do
94
95 end subroutine da_setlegpol
96
97