/
CilmPlus.f95
320 lines (282 loc) · 10.7 KB
/
CilmPlus.f95
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
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
subroutine CilmPlus(cilm, gridin, lmax, nmax, mass, d, rho, gridtype, w, zero, plx, n, dref)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! This routine will compute the potential coefficients associated
! with the input gridded relief using the method of Wieczorek and
! Phillips (1998). The input grid must contain the degree-0 term
! and the computed coefficients will be referenced to the corresonding
! degree-0 radius. Note that the array plx is optional, and should not
! be precomputed when memory is an issue (i.e., lmax>360).
!
! Calling Parameters
! IN
! gridin: Input grid to be transformed to spherical
! harmonics.
! lmax: Maximum spherical harmonic degree to compute. For gridtype 3 and 4,
! this must be less than or equal to N/2 - 1.
! nmax: Order of expansion.
! mass: Mass of planet.
! rho: density of the relief.
! gridtype 1 = Gauss-Legendre quadrature grid corresponding to LMAX.
! 2 = N by N Driscoll and Healy grid corresponding to LMAX.
! 3 = N by 2N Driscoll and Healy grid corresponding to LMAX.
! (4 = 2D Cartesian using MakeGrid2D is not implemented).
! OUT
! cilm: Output spherical harmonic coefficients with dimensions (2, lmax+1, lmax+1).
! d: The radius that the coefficients are referenced
! to. This parameter corresponds to the degree zero term of the data.
! OPTIONAL
! w: Gauss-Legendre points used in the integrations of dimension lmax+1.
! zero: Array of dimension lmax+1 that contains the latitudinal
! gridpoints used in the Gauss-Legendre quadrature integration
! scheme. Only needed when PLX is not included.
! (Determined from a call to PreCompute).
! plx: Input array of Associated Legendre Polnomials computed
! at the Gauss points (determined from a call to
! PreCompute). If this is not included, then the optional
! array zero MUST be inlcuded.
! N: Number of latitude points in the Driscoll and Healy grid. Required for Gridtype 2 and 3
! dref: The reference radius used to be used when calculating the spherical
! harmonic coefficients. If not specified, this will be set equal to the
! mean radius of GRIDIN.
!
! All units assumed to be SI.
!
! Dependencies: NGLQSH, SHExpandGLQ, SHExpandDH
!
! Written by Mark Wieczorek 2003
! September 3, 2005. Modifed so that the array plx is now optional.
! April 2009. Added the option to use N by N and N by 2N Driscoll and Healy (1994) grids.
! GRIDTYPE must now be specified, and W is now defined as an optional parameter.
! Added the optional parameter DREF, which sets the reference radius of the expansion.
! April 2012. Modified to calculate GRID**k before sending in to the subroutine call in order to improve
! memory management
! November 2, 2012. Modified the way in which the grid was scaled before cacluating the spherical
! harmonic transform.
! Jan 2013. Modified FACT so that it returned a real*8 number, and not INTEGER. The integer
! would overflow for powers beyond about 13.
!
! Copyright (c) 2005-2012, Mark A. Wieczorek
! All rights reserved.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use SHTOOLS, only: NGLQSH, SHExpandGLQ, SHExpandDH
implicit none
real*8, intent(in) :: gridin(:,:), mass, rho
real*8, intent(in), optional :: w(:), zero(:), plx(:,:), dref
real*8, intent(out) :: cilm(:,:,:), d
integer, intent(in) :: lmax, nmax, gridtype
integer, intent(in), optional :: n
real*8 :: prod, pi, scalef
real*8, allocatable :: cilmn(:, :, :), grid(:, :)
integer :: j, l, k, nlat, nlong, astat(2), lmax_dh
pi = acos(-1.0d0)
if (size(cilm(:,1,1)) < 2 .or. size(cilm(1,:,1)) < lmax+1 .or. size(cilm(1,1,:)) < lmax+1) then
print*, "Error --- CilmPlus"
print*, "CILM must be dimensioned as (2, LMAX+1, LMAX+1) where LMAX is ", lmax
print*, "Input dimension is ", size(cilm(:,1,1)), size(cilm(1,:,1)), size(cilm(1,1,:))
stop
endif
if (gridtype == 4) then
print*, "Error --- CilmPlus"
print*, "GRIDTYPE 4 (Cartesian obtained from MakeGrid2D) is not allowed."
stop
elseif (gridtype == 1) then
if (present(n)) then
print*, "Error --- CilmPlus"
print*, "N can not be present when using GLQ grids."
stop
endif
if (present(w) .and. present(zero) .and. present(plx)) then
print*, "Error --- CilmPlus"
print*, "For GLQ grids, either W and ZERO or W and PLX must be present, but not all three."
stop
elseif (present(w) .and. present(zero)) then
if (size(w) < lmax + 1) then
print*, "Error --- CilmPlus"
print*, "W must be dimensioned as (LMAX+1) where LMAX is ", lmax
print*, "Input dimension is ", size(w)
stop
endif
if (size(zero) < lmax + 1) then
print*, "Error --- CilmPlus"
print*, "ZERO must be dimensioned as (LMAX+1) where LMAX is ", lmax
print*, "Input dimension is ", size(zero)
stop
endif
elseif (present(plx) .and. present(w)) then
if (size(w) < lmax + 1) then
print*, "Error --- CilmPlus"
print*, "W must be dimensioned as (LMAX+1) where LMAX is ", lmax
print*, "Input dimension is ", size(w)
stop
endif
if (size(plx(:,1)) < lmax+1 .or. size(plx(1,:)) < (lmax+1)*(lmax+2)/2) then
print*, "Error --- CilmPlus"
print*, "PLX must be dimensioned as (LMAX+1, (LMAX+1)*(LMAX+2)/2) where LMAX is ", lmax
print*, "Input dimension is ", size(plx(:,1)), size(plx(1,:))
stop
endif
else
print*, "Error --- CilmPlus"
print*, "For GLQ grids, either W and ZERO or W and PLX must be present"
stop
endif
elseif (gridtype == 2) then
if (present(w) .or. present(zero) .or. present(plx)) then
print*, "Error --- CilmPlus"
print*, "W, ZERO and PLX can not be present for Driscoll-Healy grids."
stop
elseif (.not.present(N)) then
print*, "Error --- CilmPlus"
print*, "N must be present when GRIDTYPE is 2 or 3."
stop
endif
elseif (gridtype == 3) then
if (present(w) .or. present(zero) .or. present(plx)) then
print*, "Error --- CilmPlus"
print*, "W, ZERO and PLX can not be present for Driscoll-Healy grids."
stop
elseif (.not.present(N)) then
print*, "Error --- CilmPlus"
print*, "N must be present when GRIDTYPE is 2 or 3."
stop
endif
else
print*, "Error --- CilmPlus"
print*, "GRIDTYPE must be 1 (GLQ), 2 (NxN) or 3 (Nx2N)"
print*, "Input value is ", gridtype
stop
endif
if (gridtype == 1) then
nlat = NGLQSH(lmax) ! lmax+1
nlong = 2*lmax+1
elseif (gridtype == 2) then
nlat = N
nlong = N
lmax_dh = N/2 - 1
if (lmax > lmax_dh) then
print*, "Error --- CilmPlus"
print*, "For Driscoll-Healy grids, LMAX must be less than or equal to N/2 -1, where N is ", N
print*, "Input value of LMAX is ", lmax
stop
endif
elseif (gridtype == 3) then
nlat = N
nlong = 2*N
lmax_dh = N/2 - 1
if (lmax > lmax_dh) then
print*, "Error --- CilmPlus"
print*, "For Driscoll-Healy grids, LMAX must be less than or equal to N/2 -1, where N is ", N
print*, "Input value of LMAX is ", lmax
stop
endif
endif
if (size(gridin(1,:)) < nlong .or. size(gridin(:,1)) < nlat) then
print*, "Error --- CilmPlus"
if (gridtype == 2) then
print*, "GRIDIN must be dimensioned as (LMAX+1, 2*LMAX+1) where LMAX is ", lmax
print*, "Input dimension is ", size(gridin(1,:)), size(gridin(:,1))
stop
elseif (gridtype == 3) then
print*, "GRIDIN must be dimensioned as (N, N) where N is ", n
print*, "Input dimension is ", size(gridin(1,:)), size(gridin(:,1))
stop
elseif (gridtype == 4) then
print*, "GRIDIN must be dimensioned as (N, 2N) where N is ", n
print*, "Input dimension is ", size(gridin(1,:)), size(gridin(:,1))
stop
endif
endif
allocate(cilmn(2, lmax+1, lmax+1), stat=astat(1))
allocate(grid(nlat, nlong), stat=astat(2))
if (astat(1) /= 0 .or. astat(2) /= 0) then
print*, "Error --- CilmPlus"
print*, "Problem allocating arrays CILMN and GRID", astat(1), astat(2)
stop
endif
cilm = 0.0d0
cilmn = 0.0d0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Do the expansion.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Do k = 1 terms first
grid(1:nlat, 1:nlong) = gridin(1:nlat,1:nlong)
select case(gridtype)
case(1)
if (present(plx)) then
call SHExpandGLQ(cilmn(1:2,1:lmax+1,1:lmax+1), lmax, grid(1:nlat,1:nlong), &
w, plx = plx, norm = 1, csphase = 1)
else
call SHExpandGLQ(cilmn(1:2,1:lmax+1,1:lmax+1), lmax, grid(1:nlat,1:nlong), &
w, zero=zero, norm = 1, csphase = 1)
endif
case(2); call SHExpandDH(grid(1:nlat,1:nlong), n, cilmn(1:2,1:lmax+1,1:lmax+1), &
lmax_dh, norm = 1, sampling = 1, csphase = 1, lmax_calc = lmax)
case(3); call SHExpandDH(grid(1:nlat,1:nlong), n, cilmn(1:2,1:lmax+1,1:lmax+1), &
lmax_dh, norm = 1, sampling = 2, csphase = 1, lmax_calc = lmax)
end select
if (present(dref)) then
d = dref
else
d = cilmn(1,1,1) ! use mean radius of relief for reference sphere
endif
cilmn(1,1,1) = cilmn(1,1,1) - d
do l=0, lmax
cilm(1:2,l+1,1:l+1) = 4.0d0*pi*rho*(d**2)*cilmn(1:2,l+1,1:l+1) /mass /dble(2*l+1)
enddo
grid(1:nlat, 1:nlong) = grid(1:nlat, 1:nlong) - d
scalef = maxval(abs(grid(1:nlat,1:nlong)))
grid(1:nlat, 1:nlong) = grid(1:nlat, 1:nlong) / scalef
do k=2, nmax
grid(1:nlat,1:nlong) = grid(1:nlat,1:nlong) * ((gridin(1:nlat,1:nlong)-d)/scalef)
select case(gridtype)
case(1)
if (present(plx)) then
call SHExpandGLQ(cilmn(1:2,1:lmax+1,1:lmax+1), lmax, grid(1:nlat,1:nlong), &
w, plx=plx, norm = 1, csphase = 1)
else
call SHExpandGLQ(cilmn(1:2,1:lmax+1,1:lmax+1), lmax, grid(1:nlat,1:nlong), &
w, zero=zero, norm = 1, csphase = 1)
endif
case(2); call SHExpandDH(grid(1:nlat,1:nlong), n, cilmn(1:2,1:lmax+1,1:lmax+1), &
lmax_dh, norm = 1, sampling = 1, csphase = 1, lmax_calc = lmax)
case(3); call SHExpandDH(grid(1:nlat,1:nlong), n, cilmn(1:2,1:lmax+1,1:lmax+1), &
lmax_dh, norm = 1, sampling = 2, csphase = 1, lmax_calc = lmax)
end select
do l = 0, lmax
prod = 4*pi*rho*(d**3)/mass * (scalef/d)**k
do j=1, k
prod = prod * (l+4-j)
enddo
prod = prod/( dble(2*l+1) * dble(l+3) * dble(fact(k)) )
cilm(1:2,l+1,1:l+1) = cilm(1:2,l+1,1:l+1) + cilmn(1:2,l+1,1:l+1)*prod
enddo
enddo
deallocate(cilmn)
deallocate(grid)
contains
function fact(i)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! This function computes the factorial of an integer
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
integer :: i, j
real*8 :: fact
if (i == 0) then
fact = 1.0d0
elseif (i .lt. 0) then
print*, "Argument to FACT must be positive"
stop
else
fact = 1.0d0
do j = 1, i
fact = fact * j
enddo
endif
end function fact
end subroutine CilmPlus