-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathglatwgt_module.f90
200 lines (158 loc) · 4.6 KB
/
glatwgt_module.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
module glatwgt_module
! calculates Gaussian points and weights.
use kind_module, only: i4b, dp
use math_module, only: pi=>math_pi
! implicit none
private
! Source: Based on Swartztrauber (2003)
! Author: T. Enomoto
! Usage:
! Calculate Gaussian points and weights
! subroutine f_scgaus(x, w, jMax)
! glat(jMax) cos(Gaussian colatitudes) = sin(Gaussian latitudes)
! gwt(jMax) Gaussian weights
! History:
! TE 27 Apr 2003 Fixed a bug in setting SH colatitudes
! TE 28 Apr 2003 Ported for AFES
! TE 24 Apr 2003 Implemented Fourier-Legendre formulation.
! xacc_min double minimum accuracy available
real(kind=dp), private, parameter :: xacc_min = 1.0e-15_dp
integer(kind=i4b), private :: jMid
real(kind=dp), private, dimension(:), allocatable :: an
private :: legendre_init, legendre_P, legendre_dp, legendre_clean, newton
public :: glatwgt_calc
contains
! Public procedures
subroutine glatwgt_calc(glat,gwgt)
implicit none
! calculates sin(Gaussian latitudes) between 1 and -1
! and Gaussian weights.
! NB. Gaussian colatitudes are used during calculation
real(kind=dp), dimension(:), intent(inout) :: glat, gwgt
integer(kind=i4b) :: l, j, jj=1, jMax
real(kind=dp) :: guess, pn, dpn, pn_max=0.0_dp, s
jMax = size(glat)
jMid = jMax/2
call legendre_init()
guess = 0.5_dp*pi*(1.0_dp-1.0_dp/(jMax + 1))
call newton(legendre_P, legendre_dP, guess, glat(jMid))
guess = 3.0_dp*glat(jMid) - pi
call newton(legendre_P, legendre_dP, guess, glat(jMid-1))
do l = jMid-2, 1, -1
guess = 2*glat(l+1) - glat(l+2)
call newton(legendre_P, legendre_dP, guess, glat(l))
end do
do j = 1, jMid
call legendre_dP(glat(j), dpn)
gwgt(j) = (2.0_dp*jMax + 1.0_dp)/(dpn)**2
end do
! ***** debug
! print *, "j, lat, pn, w"
! do j=1, jMid
! call legendre_P(glat(j), pn)
! print *, j, (pi/2-glat(j))*180/pi, pn, gwgt(j)
! pn = abs(pn)
! if (pn>pn_max) then
! pn_max = pn
! jj = j
! end if
! end do
! print *, "Largest error:", pn_max, " at", jj
s = sum(gwgt(1:jMid))
print *, "sum of weights:", s, " error=", abs(1.0_dp - s)
! *****
gwgt(jMax:jMid+1:-1) = gwgt(1:jMid)
glat(1:jMid) = cos(glat(1:jMid))
glat(jMax:jMid+1:-1) = -glat(1:jMid)
call legendre_clean()
end subroutine glatwgt_calc
! private procedures
subroutine legendre_init()
implicit none
real(kind=dp), parameter :: a0sq = 1.5_dp
integer(kind=i4b) :: k, l, n
n = jMid*2
allocate(an(0:jMid))
an(jMid) = a0sq
do k=2, n
an(jMid) = (1.0_dp-1.0_dp/(4.0_dp*k*k))*an(jMid)
end do
an(jMid) = sqrt(an(jMid))
do k=1, jMid
l = 2*k
an(jMid-k) = (l-1.0_dp)*(2.0_dp*n-l+2.0_dp)/(l*(2.0_dp*n-l+1.0_dp)) * an(jMid-k+1)
end do
an(0) = 0.5_dp*an(0)
end subroutine legendre_init
subroutine legendre_P(theta, pn)
implicit none
real(kind=dp), intent(in) :: theta
real(kind=dp), intent(out) :: pn
integer(kind=i4b) :: k, l
pn = 0.0
do l=0, jMid
k=l*2 ! k = l*2 + 1 if n odd
pn = pn + an(l) * cos(k*theta)
end do
end subroutine legendre_P
subroutine legendre_dP(theta, dpn)
implicit none
real(kind=dp), intent(in) :: theta
real(kind=dp), intent(out) :: dpn
integer(kind=i4b) :: k, l
dpn = 0.0
do l=1, jMid
k=l*2 ! k = l*2 + 1 if n odd
dpn = dpn - k * an(l) * sin(k*theta)
end do
end subroutine legendre_dP
subroutine legendre_clean()
implicit none
deallocate(an)
end subroutine legendre_clean
subroutine newton(f, df, x0, x, tolerance)
use kind_module, only: dp, i4b
implicit none
! finds the root u
interface
subroutine f(x, f_result)
use kind_module, only: dp
real(kind=dp), intent(in) :: x
real(kind=dp), intent(out) :: f_result
end subroutine f
subroutine df(x, f_result)
use kind_module, only: dp
real(kind=dp), intent(in) :: x
real(kind=dp), intent(out) :: f_result
end subroutine df
end interface
real(kind=dp), intent(in) :: x0
real(kind=dp), intent(out) :: x
real(kind=dp), optional, intent(in) :: tolerance
integer(kind=i4b), parameter :: newton_max = 500
real(kind=dp) :: xacc = xacc_min
real(kind=dp) :: y, dy
integer(kind=i4b) :: i
if (present(tolerance)) then
if (tolerance < xacc_min) then
print *, "### Error in newton: tolerance too small"
stop
else
xacc = tolerance
end if
else
xacc = xacc_min
end if
x = x0
do i = 1, newton_max
call f(x, y)
call df(x, dy)
y = y/dy
if (abs(y) < xacc) then
return
end if
x = x - y
end do
print *, "### Error in newton : Too many refinement."
end subroutine newton
end module glatwgt_module