/
GLQGridCoord.f95
121 lines (107 loc) · 3.8 KB
/
GLQGridCoord.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
subroutine GLQGridCoord(latglq, longlq, lmax, nlat, nlong, extend, exitstatus)
!------------------------------------------------------------------------------
!
! Given a maximum spherical harmonic degree lmax, this routine
! will determine the latitude and longitude coordinates associated with
! grids that are used in the Gauss-Legendre quadratue spherical harmonic
! expansion routines. The coordinates are output in DEGREES.
!
! Calling Parameters
! IN
! lmax Maximum spherical harmonic degree of the expansion.
!
! OUT
! latglq Array of latitude points used in Gauss-Legendre
! grids, in degrees.
! longlq Array of longitude points used in Gauss-Legendre
! grids, in degrees.
! nlat Number of latitude points.
! nlong Number of longitude points.
!
! OPTIONAL (IN)
! extend If 1, return a grid that contains an additional column
! for 360 E longitude.
!
! OPTIONAL (OUT)
! exitstatus If present, instead of executing a STOP when an error
! is encountered, the variable exitstatus will be
! returned describing the error.
! 0 = No errors;
! 1 = Improper dimensions of input array;
! 2 = Improper bounds for input variable;
! 3 = Error allocating memory;
! 4 = File IO error.
!
! Copyright (c) 2005-2019, SHTOOLS
! All rights reserved.
!
!------------------------------------------------------------------------------
use SHTOOLS, only: PreGLQ
use ftypes
implicit none
integer(int32), intent(in) :: lmax
integer(int32), intent(out) :: nlat, nlong
real(dp), intent(out) :: latglq(:), longlq(:)
integer(int32), intent(in), optional :: extend
integer(int32), intent(out), optional :: exitstatus
real(dp) :: pi, upper, lower, zero(lmax+1), w(lmax+1)
integer(int32) :: i, nlong_out
if (present(exitstatus)) exitstatus = 0
nlat = lmax + 1
nlong = 2 * lmax + 1
if (present(extend)) then
if (extend == 0) then
nlong_out = nlong
else if (extend == 1) then
nlong_out = nlong + 1
else
print*, "Error --- GLQGridCoord"
print*, "Optional parameter EXTEND must be 0 or 1."
print*, "Input value is ", extend
if (present(exitstatus)) then
exitstatus = 2
return
else
stop
end if
end if
else
nlong_out = nlong
end if
if (size(latglq) < nlat) then
print*, "Error --- GLQGridCoord"
print*, "LATGLQ must be dimensioned as (LMAX+1) where LMAX is ", lmax
print*, "Input array is dimensioned as ", size(latglq)
if (present(exitstatus)) then
exitstatus = 1
return
else
stop
end if
else if (size(longlq) < nlong_out) then
print*, "Error --- GLQGridCoord"
print*, "LONGLQ must be dimensioned as ", nlong_out
print*, "Input array is dimensioned as ", size(longlq)
if (present(exitstatus)) then
exitstatus = 1
return
else
stop
end if
end if
pi = acos(-1.0_dp)
upper = 1.0_dp
lower = -1.0_dp
if (present(exitstatus)) then
call PreGLQ(lower, upper, nlat, zero, w, exitstatus = exitstatus)
if (exitstatus /= 0) return
else
call PreGLQ(lower, upper, nlat, zero, w)
end if
do i = 1, nlong_out
longlq(i) = 360.0_dp * dble(i-1) / dble(nlong)
end do
do i = 1, nlat
latglq(i) = asin(zero(i)) * 180.0_dp / pi
end do
end subroutine GLQGridCoord