-
Notifications
You must be signed in to change notification settings - Fork 0
/
netcdf_output.f90
369 lines (290 loc) · 17.6 KB
/
netcdf_output.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
module netcdf_output
use netcdf
implicit none
contains
subroutine create_output_file(fname, x, y, xbounds, ybounds, flip_x, flip_y, fid, varid, nrec, earth_model, radius, flattening,&
grid_area, area_fraction, ifname)
use miscellaneous_functions, only: utc_time_string
use geography, only: get_earth_parameters
character(len=*), intent(in):: fname
double precision, dimension(:), intent(in):: x, y, xbounds, ybounds
logical, intent(in):: flip_x, flip_y
integer, intent(out):: fid, varid
integer, intent(in), optional:: nrec
character(len=*), intent(in), optional:: earth_model
double precision, intent(in), optional:: radius, flattening
double precision, dimension(:,:), intent(in), optional:: grid_area
logical, intent(in), optional:: area_fraction
character(len=*), optional:: ifname
double precision:: loc_radius, loc_flattening, eccentricity
logical:: fract, write_earth_param
character(len=1000):: command_line
integer:: dimid(4), loc_varid(7), nx, ny, ierr, i
if (present(area_fraction)) then
fract = area_fraction
else
fract = .false.
end if
! Data dimension:
! ---------------
nx = size(x)
ny = size(y)
! Create output file
! ------------------
ierr = nf90_create(fname, NF90_CLOBBER, fid)
call nf90check(ierr, 'Error while creating netCDF output file '//trim(fname))
! Global attributes
! ----------------------
ierr = nf90_put_att(fid, NF90_GLOBAL, 'title', 'Area fraction of polygons on a longitude-latitude grid')
call nf90check(ierr, 'Error while putting attribute "title" in '//trim(fname))
if (present(earth_model)) then
call get_earth_parameters(earth_model, loc_radius, loc_flattening, eccentricity)
ierr = nf90_put_att(fid, NF90_GLOBAL, 'earth_model', earth_model)
call nf90check(ierr, 'Error while putting attribute "earth_model" in '//trim(fname))
write_earth_param = .true.
elseif (present(radius) .and. present(flattening)) then
loc_radius = radius
loc_flattening = flattening
eccentricity = sqrt(1d0 - (1d0-loc_flattening)**2)
write_earth_param = .true.
else
write_earth_param = .false.
end if
if (write_earth_param) then
ierr = nf90_put_att(fid, NF90_GLOBAL, 'earth_equatorial_radius', loc_radius)
call nf90check(ierr, 'Error while putting attribute "earth_equatorial_radius" in '//trim(fname))
ierr = nf90_put_att(fid, NF90_GLOBAL, 'earth_equatorial_radius_units', 'm')
call nf90check(ierr, 'Error while putting attribute "earth_equatorial_radius_units" in '//trim(fname))
ierr = nf90_put_att(fid, NF90_GLOBAL, 'earth_polar_radius', loc_radius*(1-loc_flattening))
call nf90check(ierr, 'Error while putting attribute "earth_equatorial_radius" in '//trim(fname))
ierr = nf90_put_att(fid, NF90_GLOBAL, 'earth_polar_radius_units', 'm')
call nf90check(ierr, 'Error while putting attribute "earth_polar_radius_units" in '//trim(fname))
ierr = nf90_put_att(fid, NF90_GLOBAL, 'earth_flattening', loc_flattening)
call nf90check(ierr, 'Error while putting attribute "earth_flattening" in '//trim(fname))
ierr = nf90_put_att(fid, NF90_GLOBAL, 'earth_first_eccentricity', eccentricity)
call nf90check(ierr, 'Error while putting attribute "earth_eccentricity" in '//trim(fname))
end if
if (present(ifname)) then
ierr = nf90_put_att(fid, NF90_GLOBAL, 'source_shapefile', ifname)
call nf90check(ierr, 'Error while putting attribute "source_shapefile" in '//trim(fname))
end if
call get_command(command_line)
ierr = nf90_put_att(fid, NF90_GLOBAL, 'history', utc_time_string()//': '//trim(command_line))
call nf90check(ierr, 'Error while putting attribute "history" in '//trim(fname))
! Define dimensions and variables
! -------------------------------
! Dimensions
ierr = nf90_def_dim(fid, name='longitude', len=nx, dimid=dimid(1))
call nf90check(ierr, 'Error while defining dimension "longitude" in '//trim(fname))
ierr = nf90_def_dim(fid, name='latitude', len=ny, dimid=dimid(2))
call nf90check(ierr, 'Error while defining dimension "latitude" in '//trim(fname))
ierr = nf90_def_dim(fid, name='bounds', len=2, dimid=dimid(3))
call nf90check(ierr, 'Error while defining dimension "latitude" in '//trim(fname))
if (present(nrec)) then
ierr = nf90_def_dim(fid, name='records', len=nrec, dimid=dimid(4))
call nf90check(ierr, 'Error while defining dimension "records" in '//trim(fname))
end if
! Define lon variable
ierr = nf90_def_var(fid, 'longitude', NF90_REAL, dimid(1:1), loc_varid(1))
call nf90check(ierr, 'Error while defining variable "longitude" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(1), 'axis', 'X')
call nf90check(ierr, 'Error while putting attribute "axis" in variable "longitude" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(1), 'long_name', 'longitude')
call nf90check(ierr, 'Error while putting attribute "long_name" in variable "longitude" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(1), 'standard_name', 'longitude')
call nf90check(ierr, 'Error while putting attribute "standard_name" in variable "longitude" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(1), 'units', 'degrees_east')
call nf90check(ierr, 'Error while putting attribute "units" in variable "longitude" in '//trim(fname))
! Define lat variable
ierr = nf90_def_var(fid, 'latitude', NF90_REAL, dimid(2:2), loc_varid(2))
call nf90check(ierr, 'Error while defining variable "latitude" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(2), 'axis', 'Y')
call nf90check(ierr, 'Error while putting attribute "axis" in variable "latitude" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(2), 'long_name', 'latitude')
call nf90check(ierr, 'Error while putting attribute "long_name" in variable "latitude" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(2), 'standard_name', 'latitude')
call nf90check(ierr, 'Error while putting attribute "standard_name" in variable "latitude" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(2), 'units', 'degrees_north')
call nf90check(ierr, 'Error while putting attribute "units" in variable "latitude" in '//trim(fname))
! Define bounds variable
ierr = nf90_def_var(fid, 'bounds', NF90_INT, dimid(3:3), loc_varid(3))
call nf90check(ierr, 'Error while defining variable "bounds" in '//trim(fname))
! Define record variable
if (present(nrec)) then
ierr = nf90_def_var(fid, 'records', NF90_INT, dimid(4:4), loc_varid(4))
call nf90check(ierr, 'Error while defining variable "records" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(4), 'long_name', 'record number in source shapefile')
call nf90check(ierr, 'Error while putting attribute "long_name" in variable "records" in '//trim(fname))
end if
! Define lon bounds variable
ierr = nf90_def_var(fid, 'longitude_bounds', NF90_REAL, (/dimid(1),dimid(3)/), loc_varid(5))
call nf90check(ierr, 'Error while defining variable "longitude_bounds" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(5), 'long_name', 'longitude boundaries')
call nf90check(ierr, 'Error while putting attribute "long_name" in variable "longitude_bounds" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(5), 'standard_name', 'lon_bounds')
call nf90check(ierr, 'Error while putting attribute "standard_name" in variable "longitude_bounds" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(5), 'units', 'degrees_east')
call nf90check(ierr, 'Error while putting attribute "units" in variable "longitude_bounds" in '//trim(fname))
! Define lat bounds variable
ierr = nf90_def_var(fid, 'latitude_bounds', NF90_REAL, dimid(2:3), loc_varid(6))
call nf90check(ierr, 'Error while defining variable "latitude_bounds" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(6), 'long_name', 'latitude boundaries')
call nf90check(ierr, 'Error while putting attribute "long_name" in variable "latitude_bounds" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(6), 'standard_name', 'lat_bounds')
call nf90check(ierr, 'Error while putting attribute "standard_name" in variable "latitude_bounds" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(6), 'units', 'degrees_north')
call nf90check(ierr, 'Error while putting attribute "units" in variable "latitude_bounds" in ' &
//trim(fname))
! Define area variable
if (present(grid_area)) then
ierr = nf90_def_var(fid, 'grid_area', NF90_REAL, dimid(1:2), loc_varid(7))
call nf90check(ierr, 'Error while defining variable "grid_area" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(7), 'long_name', 'grid cell area')
call nf90check(ierr, 'Error while putting attribute "long_name" in variable "area" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(7), 'standard_name', 'area')
call nf90check(ierr, 'Error while putting attribute "standard_name" in variable "area" in '//trim(fname))
ierr = nf90_put_att(fid, loc_varid(7), 'units', 'm2')
call nf90check(ierr, 'Error while putting attribute "units" in variable "area" in '//trim(fname))
end if
! Define polygon fraction variable
if (fract) then
if (present(nrec)) then
ierr = nf90_def_var(fid, 'polygon_fraction', NF90_REAL, (/dimid(1:2),dimid(4)/), varid)
else
ierr = nf90_def_var(fid, 'polygon_fraction', NF90_REAL, dimid(1:2), varid)
end if
call nf90check(ierr, 'Error while defining variable "polygon_fraction" in '//trim(fname))
ierr = nf90_put_att(fid, varid, 'long_name', 'fraction of grid cell covered by polygon')
call nf90check(ierr, 'Error while putting attribute "long_name" in variable "polygon_fraction" in '//trim(fname))
ierr = nf90_put_att(fid, varid, 'units', '-')
call nf90check(ierr, 'Error while putting attribute "units" in variable "polygon_fraction" in '//trim(fname))
else
if (present(nrec)) then
ierr = nf90_def_var(fid, 'polygon_area', NF90_REAL, (/dimid(1:2),dimid(4)/), varid)
else
ierr = nf90_def_var(fid, 'polygon_area', NF90_REAL, dimid(1:2), varid)
end if
call nf90check(ierr, 'Error while defining variable "polygon_area" in '//trim(fname))
ierr = nf90_put_att(fid, varid, 'long_name', 'area grid cell covered by polygon')
call nf90check(ierr, 'Error while putting attribute "long_name" in variable "polygon_area" in '//trim(fname))
ierr = nf90_put_att(fid, varid, 'units', 'm2')
call nf90check(ierr, 'Error while putting attribute "units" in variable "polygon_area" in '//trim(fname))
end if
! End definition mode:
ierr = nf90_enddef(fid)
call nf90check(ierr, 'error while end of definition of file '//trim(fname))
! Put variable
! ------------
! Longitude
if (flip_x) then
ierr = nf90_put_var(fid, loc_varid(1), x(nx:1:-1))
else
ierr = nf90_put_var(fid, loc_varid(1), x)
end if
call nf90check(ierr, 'Error while putting variable "longitude" in '//trim(fname))
! Latitude
if (flip_y) then
ierr = nf90_put_var(fid, loc_varid(2), y(ny:1:-1))
else
ierr = nf90_put_var(fid, loc_varid(2), y)
end if
call nf90check(ierr, 'Error while putting variable "latitude" in '//trim(fname))
! Bounds
ierr = nf90_put_var(fid, loc_varid(3), (/1,2/))
call nf90check(ierr, 'Error while putting variable "bounds" in '//trim(fname))
! Records
if (present(nrec)) then
ierr = nf90_put_var(fid, loc_varid(4), (/(i,i=1,nrec)/))
call nf90check(ierr, 'Error while putting variable "record" in '//trim(fname))
end if
! Longitude bounds
if (flip_x) then
ierr = nf90_put_var(fid, loc_varid(5), xbounds(nx+1:2:-1), start=(/1,1/), count=(/nx,1/))
ierr = nf90_put_var(fid, loc_varid(5), xbounds(nx:1:-1), start=(/1,2/), count=(/nx,1/))
else
ierr = nf90_put_var(fid, loc_varid(5), xbounds(1:nx), start=(/1,1/), count=(/nx,1/))
ierr = nf90_put_var(fid, loc_varid(5), xbounds(2:nx+1), start=(/1,2/), count=(/nx,1/))
end if
call nf90check(ierr, 'Error while putting variable "longitude_bounds" in '//trim(fname))
! Latitude bounds
if (flip_y) then
ierr = nf90_put_var(fid, loc_varid(6), ybounds(ny+1:2:-1), start=(/1,1/), count=(/ny,1/))
ierr = nf90_put_var(fid, loc_varid(6), ybounds(ny:1:-1), start=(/1,2/), count=(/ny,1/))
else
ierr = nf90_put_var(fid, loc_varid(6), ybounds(1:ny), start=(/1,1/), count=(/ny,1/))
ierr = nf90_put_var(fid, loc_varid(6), ybounds(2:ny+1), start=(/1,2/), count=(/ny,1/))
end if
call nf90check(ierr, 'Error while putting variable "latitude_bounds" in '//trim(fname))
! Area
if (present(grid_area)) then
ierr = nf90_put_var(fid, loc_varid(7), grid_area)
call nf90check(ierr, 'Error while putting variable "area" in '//trim(fname))
end if
end subroutine
!======================================================================!
subroutine write_output(fid, varid, nx, ny, flip_x, flip_y, i0, outvar)
integer, intent(in):: fid, varid, nx, ny, i0
logical, intent(in):: flip_x, flip_y
double precision, dimension(nx,ny), intent(in):: outvar
integer:: ierr
if (flip_x .and. flip_y) then
ierr = nf90_put_var(fid, varid, outvar(i0-1:1:-1, ny:1:-1), start=(/1,1/), count=(/i0-1,ny/))
ierr = nf90_put_var(fid, varid, outvar(nx:i0:-1, ny:1:-1), start=(/i0,1/), count=(/nx+1-i0,ny/))
elseif (flip_x) then
ierr = nf90_put_var(fid, varid, outvar(i0-1:1:-1, :), start=(/1,1/), count=(/i0-1,ny/))
ierr = nf90_put_var(fid, varid, outvar(nx:i0:-1, :), start=(/i0,1/), count=(/nx+1-i0,ny/))
elseif (flip_y) then
ierr = nf90_put_var(fid, varid, outvar(i0:nx, ny:1:-1), start=(/1,1/), count=(/nx+1-i0,ny/))
ierr = nf90_put_var(fid, varid, outvar(1:i0-1, ny:1:-1), start=(/nx+2-i0,1/), count=(/i0-1,ny/))
else
ierr = nf90_put_var(fid, varid, outvar(i0:nx, :), start=(/1,1/), count=(/nx+1-i0,ny/))
ierr = nf90_put_var(fid, varid, outvar(1:i0-1, :), start=(/nx+2-i0,1/), count=(/i0-1,ny/))
end if
call nf90check(ierr, 'Error while putting variable')
end subroutine
!======================================================================!
subroutine write_output_slice(fid, varid, nx, ny, flip_x, flip_y, i0, outvar, recnum)
integer, intent(in):: fid, varid, recnum, nx, ny, i0
logical, intent(in):: flip_x, flip_y
double precision, dimension(nx,ny), intent(in):: outvar
integer:: ierr
if (flip_x .and. flip_y) then
ierr = nf90_put_var(fid, varid, outvar(i0-1:1:-1, ny:1:-1), start=(/1,1,recnum/), count=(/i0-1,ny,1/))
ierr = nf90_put_var(fid, varid, outvar(nx:i0:-1, ny:1:-1), start=(/i0,1,recnum/), count=(/nx+1-i0,ny,1/))
elseif (flip_x) then
ierr = nf90_put_var(fid, varid, outvar(i0-1:1:-1, :), start=(/1,1,recnum/), count=(/i0-1,ny,1/))
ierr = nf90_put_var(fid, varid, outvar(nx:i0:-1, :), start=(/i0,1,recnum/), count=(/nx+1-i0,ny,1/))
elseif (flip_y) then
ierr = nf90_put_var(fid, varid, outvar(i0:nx, ny:1:-1), start=(/1,1,recnum/), count=(/nx+1-i0,ny,1/))
ierr = nf90_put_var(fid, varid, outvar(1:i0-1, ny:1:-1), start=(/nx+2-i0,1,recnum/), count=(/i0-1,ny,1/))
else
ierr = nf90_put_var(fid, varid, outvar(i0:nx, :), start=(/1,1,recnum/), count=(/nx+1-i0,ny,1/))
ierr = nf90_put_var(fid, varid, outvar(1:i0-1, :), start=(/nx+2-i0,1,recnum/), count=(/i0-1,ny,1/))
end if
call nf90check(ierr, 'Error while putting variable slice')
end subroutine
!======================================================================!
subroutine close_file(fid)
integer, intent(in):: fid
integer:: ierr
ierr = nf90_close(fid)
call nf90check(ierr, 'Error while closing file')
end subroutine
!======================================================================!
subroutine nf90check(ierr, message, nokill)
integer, intent(in):: ierr
character(len=*), intent(in):: message
logical, intent(in), optional:: nokill
if (ierr/=NF90_NOERR) then
print *, message
print *, nf90_strerror(ierr)
if (present(nokill)) then
if (.not. nokill) then
stop 2
end if
else
stop 2
end if
end if
end subroutine
end module