-
Notifications
You must be signed in to change notification settings - Fork 103
/
Curve2Mask.f95
310 lines (268 loc) · 10.2 KB
/
Curve2Mask.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
subroutine Curve2Mask(dhgrid, n, sampling, profile, nprofile, NP, extend, &
exitstatus)
!------------------------------------------------------------------------------
!
! Given a list of coordinates of a single closed curve, this routine
! will output a Driscoll and Healy sampled grid where the interior and
! exterior of the curve are filled with zeros and ones. The value at the
! north pole (either 0 or 1) is specified by the input parameter NP.
!
! Longitudes can span the range from -360 to 720 degrees. If the longitudes
! of two adjacent points differ by more than 180 degrees, it will be assumed
! that the curve passes from 360 to 0 degrees, or from -180 to 180 degrees.
!
! Calling Parameters
!
! IN
! profile Latitude (:,1) and longitude (:,2) coordinates of a
! single close curve having dimension (nprofile, 2).
! nprofile Number of coordinates in the vector profile.
! np Value of the output function at the North Pole, which
! can be either 0 or 1.
! n Number of latitude bands in the output Driscoll and
! Healy sampled grid.
! sampling 1 sets the number of longitude bands equal to 1,
! whereas 2 sets the number to twice that of n.
!
! OUT
! dhgrid A Driscoll and Healy sampled grid specifiying whether
! the point is in the curve (1), or outside of it (0).
!
! OPTIONAL (IN)
! extend If 1, return a grid that contains an additional column
! and row corresponding to 360 E longitude and 90 S
! latitude, respectively.
!
! 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 ftypes
implicit none
integer, intent(out) :: dhgrid(:,:)
real(dp), intent(in) :: profile(:,:)
integer(int32), intent(in) :: n, sampling, nprofile, np
integer(int32), intent(in), optional :: extend
integer(int32), intent(out), optional :: exitstatus
integer(int32), parameter :: maxcross = 2000
integer(int32) :: i, j, k, k_loc, nlat, nlong, numcross, next, ind1, &
ind2, nlat_out, nlong_out, extend_grid
real(dp) :: lat_int, long_int, lon, cross(maxcross), &
cross_sort(maxcross), lat1, lat2, lon1, lon2
if (present(exitstatus)) exitstatus = 0
nlat = n
lat_int = 180.0_dp / dble(nlat)
dhgrid = 0
if (mod(n,2) /= 0) then
print*, "Error --- Curve2Mask"
print*, "N must be even"
print*, "N = ", n
if (present(exitstatus)) then
exitstatus = 2
return
else
stop
end if
end if
if (sampling == 1) then
nlong = nlat
long_int = 2.0_dp * lat_int
else if (sampling == 2) then
nlong = 2 * nlat
long_int = lat_int
else
print*, "Error --- Curve2Mask"
print*, "SAMPLING of DHGRID must be 1 (equally sampled) or 2 (equally spaced)."
print*, "SAMPLING = ", sampling
if (present(exitstatus)) then
exitstatus = 2
return
else
stop
end if
end if
if (present(extend)) then
if (extend == 0) then
extend_grid = 0
nlat_out = nlat
nlong_out = nlong
else if (extend == 1) then
extend_grid = 1
nlat_out = nlat + 1
nlong_out = nlong + 1
else
print*, "Error --- Curve2Mask"
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
extend_grid = 0
nlat_out = nlat
nlong_out = nlong
end if
if (NP /= 1 .and. NP /= 0) then
print*, "Error --- Curve2Mask"
print*, "NP must be 0 if the North pole is outside of curve,"
print*, "or 1 if the North pole is inside of the curve."
print*, "NP = ", np
if (present(exitstatus)) then
exitstatus = 2
return
else
stop
end if
end if
if (size(dhgrid(1,:)) < nlong_out .or. size(dhgrid(:,1)) < nlat_out) then
print*, "Error --- Curve2Mask"
print*, "DHGRID must be dimensioned as (NLAT, NLONG)."
print*, "NLAT = ", nlat_out
print*, "NLONG = ", nlong_out
print*, "Size of GRID = ", size(dhgrid(:,1)), size(dhgrid(1,:))
if (present(exitstatus)) then
exitstatus = 1
return
else
stop
end if
end if
if (size(profile(:,1)) < nprofile .or. size(profile(1,:)) < 2) then
print*, "Error --- Curve2Mask"
print*, "PROFILE must be dimensioned as (NPROFILE, 2)."
print*, "Dimension of NPROFILE = ", size(profile(:,1)), &
size(profile(1,:))
if (present(exitstatus)) then
exitstatus = 1
return
else
stop
end if
end if
!--------------------------------------------------------------------------
!
! Start at 90N and 0E. Determine where the curve crosses in this
! longitude band, sort the values, and then set the pixels between zero
! crossings to either 0 or 1.
!
!--------------------------------------------------------------------------
do j = 1, nlong
lon = dble(j-1) * long_int
numcross = 0
do i = 1, nprofile
lat1 = profile(i, 1)
lon1 = profile(i, 2)
if (i /= nprofile) then
lat2 = profile(i+1, 1)
lon2 = profile(i+1, 2)
else
lat2 = profile(1, 1)
lon2 = profile(1, 2)
end if
! The output grid will be for longitudes between 0 and 360 degrees.
! Check if both longitudes are less than 0 or greater than 360,
! and then convert to this range.
if (lon1 < 0._dp .and. lon2 < 0._dp) then
lon1 = lon1 + 360._dp
lon2 = lon2 + 360._dp
end if
if (lon1 > 360._dp .and. lon2 > 360._dp) then
lon1 = lon1 - 360._dp
lon2 = lon2 - 360._dp
end if
! Depending on how the longitude coordinates are defined, adjacent
! points may differ by almost 360 degrees. Check for such jumps by
! looking for differences between adjacent points of more than
! 180 degrees.
if (abs(lon1 - lon2) > 180._dp) then
if (lon1 < 0._dp .or. lon2 < 0._dp) then
! There is a jump from -180 to 180
if (lon1 < 0._dp) lon1 = lon1 + 360._dp
if (lon2 < 0._dp) lon2 = lon2 + 360._dp
elseif (lon1 >= 0._dp .and. lon2 >= 0._dp) then
! There is a jump from 360 to 0. Set the largest value
! to a negative longitude.
if (lon1 > lon2) then
lon1 = lon1 - 360._dp
else
lon2 = lon2 - 360._dp
end if
end if
end if
! Find the latitude crossing by interpolating between the two
! points.
if (lon1 <= lon .and. lon2 > lon) then
numcross = numcross + 1
if (numcross > maxcross) then
print*, "Error --- Curve2Mask"
print*, "Internal variable MAXCROSS needs to be increased."
print*, "MAXCROSS = ", maxcross
if (present(exitstatus)) then
exitstatus = 5
return
else
stop
end if
end if
cross(numcross) = lat1 + (lat2 - lat1) / (lon2 - lon1) &
* (lon - lon1)
else if (lon1 > lon .and. lon2 <= lon) then
numcross = numcross + 1
if (numcross > maxcross) then
print*, "Error --- Curve2Mask"
print*, "Internal variable MAXCROSS needs to be increased."
print*, "MAXCROSS = ", maxcross
if (present(exitstatus)) then
exitstatus = 5
return
else
stop
end if
end if
cross(numcross) = lat2 + (lat1 - lat2) / (lon1 - lon2) &
* (lon - lon2)
end if
end do
if (numcross == 0) then
dhgrid(1:nlat_out, j) = np
cycle
else ! sort crossings by decreasing latitude
do k = 1, numcross
k_loc = maxloc(cross(1:numcross), 1)
cross_sort(k) = cross(k_loc)
cross(k_loc) = -999.0_dp
end do
end if
ind1 = 1
next = np
do k = 1, numcross
ind2 = nint( (90.0_dp - cross_sort(k)) / lat_int) + 1
dhgrid(ind1:ind2, j) = next
if (next == 0) then
next = 1
else
next = 0
end if
ind1 = ind2 + 1
end do
if (ind1 <= nlat_out) dhgrid(ind1:nlat_out, j) = next
end do
if (extend_grid == 1) then
dhgrid(nlat_out, 1) = dhgrid(nlat_out-1, 1)
dhgrid(1:nlat_out, nlong_out) = dhgrid(1:nlat_out, 1)
end if
end subroutine Curve2Mask