-
Notifications
You must be signed in to change notification settings - Fork 1
/
create_ss_grid.f90
415 lines (280 loc) · 12.7 KB
/
create_ss_grid.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
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
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
program create_ss_grid
! reads in a GMT formatted file containing the shear stress domains, output from QGIS,
! and creates a regular grid from a given input file.
! heavily modified from the version of create_ss_grid that I was using before
use global_parameters
! reads in the file with polygons "gmt_file.txt" and creates a grid with a given resolution set by the user
implicit none
integer :: grid_spacing, minimum_y, minimum_x, maximum_x, maximum_y, x_int, y_int
double precision :: x, y, x1, y1, x2, y2, slope, intercept, x_offset, y_offset
integer :: local_minimum_x, local_maximum_x, local_minimum_y, local_maximum_y
double precision :: local_x, local_y
character(len=255) :: input_parameter, dummy, bin_file, polygon_file, domain_max_file, domain_adjust_file, info
character(len=255) :: domain_min_file
character(len=255), parameter :: output_domains_file = "domains.txt"
character(len=255), parameter :: grid_parameters_file="ss_parameters.txt" ! parameter file used in ICESHEET
character(len=255), parameter :: output_grid_file = "shear_stress_grid.txt"
character(len=1) :: bracket
integer :: istat, number_polygons, ss_polygon_counter, point_counter, maximum_points, number_y, number_x, counter
integer :: start_y_index, start_x_index, end_y_index, end_x_index, domain_id, shear_stress, max_domain_id, current_time
integer :: local_y_counter, local_x_counter, y_counter, x_counter
integer :: time_of_maximum_ss, time_of_minimum_ss, minimum_ss
integer, parameter :: polygon_file_unit=10, max_polygons = 10000, output_domains_file_unit=20, grid_parameters_unit=30
integer, parameter :: shear_stress_max_unit=40, adjust_unit=50, output_grid_unit = 60, shear_stress_min_unit=70
integer, dimension(max_polygons) :: polygon_point_size, polygon_domain_id ! if you have more than 10000 polygons, you are ambitious!
double precision, allocatable, dimension(:,:) :: y_array, x_array
integer, dimension(max_polygons) :: shear_stress_value_array, shear_stress_min
integer, allocatable, dimension(:,:) :: grid
logical :: inside, use_max_file, use_adjust_file, use_min_file
integer, parameter :: nominal_shear_stress = 5000 ! default shear stress value
character(len=50), parameter :: output_format = "(I9,1X,I9,1X,I6)"
! requires up to three files from the command line, the GMT file is manditory
call getarg(1,polygon_file)
! read(input_parameter,*)
call getarg(2,domain_max_file)
! read(input_parameter,*)
call getarg(3,domain_adjust_file)
domain_min_file = ""
if(polygon_file == "") THEN
write(6,*) "you need to include the shear stress GMT formatted polygon file as a command line argument"
stop
end if
if(domain_max_file == "") THEN
write(6,*) "using default shear stress values from the GMT file"
use_max_file = .false.
else
write(6,*) "using user derived shear stress values"
use_max_file = .true.
endif
if(domain_adjust_file == "") THEN
write(6,*) "no time variable adjustments to shear stress requested"
use_adjust_file = .false.
else
write(6,*) "using user derived time variable shear stress adjustments"
use_adjust_file = .true.
! read in the time value (in thosands of years before present so it is an integer)
call getarg(4,dummy)
read(dummy,*) current_time
call getarg(5,domain_min_file)
endif
if(domain_min_file == "") THEN
write(6,*) "using default minimum shear stress value:", nominal_shear_stress
use_min_file = .false.
else
write(6,*) "using user derived minimum shear stress values"
use_min_file = .true.
endif
! the program expects a file called ss_parameters.txt that contains the grid parameters, x and y min/max and grid spacing
open(unit=grid_parameters_unit, file=grid_parameters_file, access="sequential", form="formatted", status="old")
read(grid_parameters_unit,*) bin_file ! not actually used at this point, will be created afterwards with GMT
read(grid_parameters_unit,*) minimum_x
read(grid_parameters_unit,*) maximum_x
read(grid_parameters_unit,*) minimum_y
read(grid_parameters_unit,*) maximum_y
read(grid_parameters_unit,*) grid_spacing
! if there is an offset to the grid (i.e. QGIS had the projection center at the center coordinate using the Lambert projection, while
! it is more convenient to use the bottom left corner as is default in GMT to keep it a regular grid
! if there is no offset, do not include it in ss_parameters.txt
read(grid_parameters_unit,*, iostat=istat) x_offset, y_offset
if(istat /= 0) THEN
x_offset = 0.0
y_offset = 0.0
endif
close(unit=grid_parameters_unit)
! first, find the number of polygons, and the maximum amount of points
number_polygons = 0
polygon_point_size = 0
shear_stress_value_array = nominal_shear_stress
max_domain_id = 0
open(unit=polygon_file_unit, file=polygon_file, access="sequential", form="formatted", status="old")
initial_read: do
read(polygon_file_unit,'(A1,A)', iostat=istat) bracket, dummy
if(istat /=0) THEN !hopefully end of file
exit initial_read
endif
if(bracket == ">") then
number_polygons = number_polygons + 1
if(number_polygons > max_polygons) THEN
write(6,*) "you have to increase the max_polygons parameter in create_ss_grid.f90, and recompile"
stop
endif
! the QGIS output file will immediately afterwards have a line containing the domain number and shear stress
read(polygon_file_unit,*) dummy, info
! replace the | character with spaces
do counter=1,len_trim(info)
if (info(counter:counter) == "|") info(counter:counter) = " "
end do
read(info,*) dummy, domain_id, shear_stress
if(domain_id > max_domain_id) THEN
max_domain_id = domain_id
end if
shear_stress_value_array(domain_id) = shear_stress
polygon_domain_id(number_polygons) = domain_id
else if (bracket == "#") then
! skip
else
polygon_point_size(number_polygons) = polygon_point_size(number_polygons) + 1
endif
end do initial_read
rewind(unit=polygon_file_unit)
maximum_points = maxval(polygon_point_size)
allocate(y_array(number_polygons,maximum_points), x_array(number_polygons,maximum_points),&
stat=istat)
if(istat /=0) THEN
write(6,*) "problem allocating arrays"
stop
endif
! read in the polygons
do ss_polygon_counter = 1, number_polygons, 1
!write(6,*) ss_polygon_counter
break_out: do
read(polygon_file_unit,*) bracket, dummy
if(bracket == ">" .or. bracket == "#") then
cycle break_out
else
backspace polygon_file_unit
exit break_out
end if
end do break_out
do point_counter = 1, polygon_point_size(ss_polygon_counter), 1
read(polygon_file_unit,*) x_array(ss_polygon_counter,point_counter), &
y_array(ss_polygon_counter,point_counter)
x_array(ss_polygon_counter,point_counter) = x_array(ss_polygon_counter,point_counter) + x_offset
y_array(ss_polygon_counter,point_counter) = y_array(ss_polygon_counter,point_counter) + y_offset
end do
end do
close(polygon_file_unit)
if(use_max_file) THEN
open(unit=shear_stress_max_unit, file=domain_max_file, access="sequential", form="formatted", status="old")
read_max: do
read(shear_stress_max_unit,*, iostat=istat) domain_id, shear_stress
if(istat /=0) THEN
exit read_max
endif
shear_stress_value_array(domain_id) = shear_stress
end do read_max
close(shear_stress_max_unit)
end if
shear_stress_min = nominal_shear_stress
if(use_min_file) THEN
open(unit=shear_stress_min_unit, file=domain_min_file, access="sequential", form="formatted", status="old")
read_min: do
read(shear_stress_min_unit,*, iostat=istat) domain_id, shear_stress
if(istat /=0) THEN
exit read_min
endif
shear_stress_min(domain_id) = shear_stress
end do read_min
close(shear_stress_min_unit)
end if
if(use_adjust_file) THEN
! adjust the maximum shear stress values based on a format
! time_of_maximum_ss time_of_minimum_ss minimum_ss
! this will linearly go between those times and adjust the shear stress based on these values, the minimum time should probably represent an ice free time
! or zero time
! if minimum_ss is set to zero, the domains_min.txt value is used
! if minimum_ss is set to a value less than zero, a minimal value is used for the whole time period. Use this to flag off times that are ice free if you want.
open(unit=adjust_unit, file=domain_adjust_file, access="sequential", form="formatted", status="old")
read_adjust: do
read(adjust_unit,*, iostat=istat) domain_id, time_of_maximum_ss, time_of_minimum_ss, minimum_ss ! must be integers
if(istat /=0) THEN
exit read_adjust
endif
! are adjustments to the current time necessary?
if(current_time >= min(time_of_maximum_ss,time_of_minimum_ss) .and. &
current_time <= max(time_of_maximum_ss,time_of_minimum_ss)) THEN
if(minimum_ss >= 0) THEN
! originally I thought of making the decay a
x1 = dble(time_of_minimum_ss)
if (minimum_ss == 0) THEN
y1 = dble(shear_stress_min(domain_id))
else
y1 = dble(minimum_ss)
endif
x2 = dble(time_of_maximum_ss)
y2 = dble(shear_stress_value_array(domain_id))
slope = (y2 - y1) / (x2 - x1)
intercept = y1 - slope * x1
shear_stress_value_array(domain_id) = nint(slope * current_time + intercept)
else ! assume that minimum value is requested
shear_stress_value_array(domain_id) = shear_stress_min(domain_id)
endif
end if
end do read_adjust
close(adjust_unit)
end if
! write out the domains file, not really used for anything else, just for checking
open(unit=output_domains_file_unit, file=output_domains_file, access="sequential", form="formatted", status="replace")
do counter = 1, max_domain_id, 1
write(output_domains_file_unit,*) counter, shear_stress_value_array(counter)
end do
! allocate the output grid
number_y = nint(dble(maximum_y - minimum_y) / dble(grid_spacing)) + 1
number_x = nint(dble(maximum_x - minimum_x) / dble(grid_spacing)) + 1
allocate(grid(number_x,number_y), stat=istat)
if(istat /=0) THEN
write(6,*) "problem allocating arrays"
stop
endif
! note that it is entirely possible that polgons will be skipped if your lat/long spacing is not dense enough
! also, if your input shapefile has any gaps between polygons, they could also be skipped.
! If they polygons overlap, the values could be overwritten
grid = nominal_shear_stress
cycle_polygon: do ss_polygon_counter = 1, number_polygons, 1
local_minimum_x = int(minval(x_array(ss_polygon_counter,1:polygon_point_size(ss_polygon_counter)))/dble(grid_spacing))&
*grid_spacing
local_maximum_x = ceiling(maxval(x_array(ss_polygon_counter,1:polygon_point_size(ss_polygon_counter)))/dble(grid_spacing))&
*grid_spacing
local_minimum_y = int(minval(y_array(ss_polygon_counter,1:polygon_point_size(ss_polygon_counter)))/dble(grid_spacing))&
*grid_spacing
local_maximum_y = ceiling(maxval(y_array(ss_polygon_counter,1:polygon_point_size(ss_polygon_counter)))/dble(grid_spacing))&
*grid_spacing
start_y_index = (local_minimum_y - minimum_y)/grid_spacing + 1
start_x_index = (local_minimum_x - minimum_x)/grid_spacing + 1
end_y_index = (local_maximum_y - minimum_y)/grid_spacing + 1
end_x_index = (local_maximum_x - minimum_x)/grid_spacing + 1
if(end_x_index < 1 .or. start_x_index > number_x .or. end_y_index < 1 .or. start_y_index > number_y) THEN ! polygon is outside of the grid
cycle cycle_polygon
end if
if (start_x_index < 1) THEN
start_x_index = 1
endif
if (end_x_index > number_x) THEN
end_x_index = number_x
endif
if (start_y_index < 1) THEN
start_y_index = 1
endif
if (end_y_index > number_y) THEN
end_y_index = number_y
endif
do local_x_counter = start_x_index, end_x_index
do local_y_counter = start_y_index, end_y_index
local_x = dble(minimum_x + (local_x_counter-1) * grid_spacing)
local_y = dble(minimum_y + (local_y_counter-1) * grid_spacing)
inside = point_in_polygon(x_array(ss_polygon_counter,1:polygon_point_size(ss_polygon_counter)), &
y_array(ss_polygon_counter,1:polygon_point_size(ss_polygon_counter)), local_x, &
local_y, polygon_point_size(ss_polygon_counter))
if(inside) THEN
grid(local_x_counter,local_y_counter) = shear_stress_value_array(polygon_domain_id(ss_polygon_counter))
endif
end do
end do
end do cycle_polygon
! output
open(unit=output_grid_unit, file=output_grid_file, access="sequential", form="formatted", status="replace")
do x_counter = 1, number_x, 1
do y_counter = 1, number_y, 1
x_int = minimum_x + (x_counter-1) * grid_spacing
y_int = minimum_y + (y_counter-1) * grid_spacing
write(output_grid_unit,output_format) x_int, y_int, grid(x_counter,y_counter)
end do
end do
close(output_grid_unit)
deallocate(y_array, x_array, grid, stat=istat)
if(istat /=0) THEN
write(6,*) "problem deallocating arrays"
stop
endif
write(6,*) "finished create_ss_grid"
end program create_ss_grid