orog_mask_tools 1.14.0
Loading...
Searching...
No Matches
io_utils.F90
Go to the documentation of this file.
1
4
8
9 module io_utils
10
11 implicit none
12
13 private
14
15 public :: qc_orog_by_ramp
16 public :: read_global_mask
17 public :: read_global_orog
18 public :: read_mask
19 public :: read_mdl_dims
20 public :: read_mdl_grid_file
21 public :: write_mask_netcdf
22 public :: write_netcdf
23
24 contains
25
41 subroutine write_netcdf(im, jm, slm, land_frac, oro, hprime, ntiles, tile, geolon, geolat, lon, lat)
42 use netcdf
43 implicit none
44 integer, intent(in):: im, jm, ntiles, tile
45 real, intent(in) :: lon(im), lat(jm)
46 real, intent(in), dimension(im,jm) :: slm, oro, geolon, geolat, land_frac
47 real, intent(in), dimension(im,jm,14):: hprime
48 character(len=128) :: outfile
49 integer :: error, ncid
50 integer :: header_buffer_val = 16384
51 integer :: fsize=65536, inital = 0
52 integer :: dim1, dim2
53 integer :: dim_lon, dim_lat
54 integer :: id_geolon,id_geolat
55 integer :: id_slmsk,id_orog_raw,id_orog_filt,id_land_frac
56 integer :: id_stddev,id_convex
57 integer :: id_oa1,id_oa2,id_oa3,id_oa4
58 integer :: id_ol1,id_ol2,id_ol3,id_ol4
59 integer :: id_theta,id_gamma,id_sigma,id_elvmax
60
61 if(ntiles > 1) then
62 write(outfile, '(a,i4.4,a)') 'out.oro.tile', tile, '.nc'
63 else
64 outfile = "out.oro.nc"
65 endif
66
67 dim1=size(lon,1)
68 dim2=size(lat,1)
69
70 !--- open the file
71 error = nf90_create(outfile, ior(nf90_netcdf4,nf90_classic_model), ncid, &
72 initialsize=inital, chunksize=fsize)
73 call netcdf_err(error, 'Creating file '//trim(outfile) )
74 !--- define dimension
75 error = nf90_def_dim(ncid, 'lon', dim1, dim_lon)
76 call netcdf_err(error, 'define dimension lon for file='//trim(outfile) )
77 error = nf90_def_dim(ncid, 'lat', dim2, dim_lat)
78 call netcdf_err(error, 'define dimension lat for file='//trim(outfile) )
79
80 !--- define field
81!---geolon
82 error = nf90_def_var(ncid, 'geolon', nf90_float, (/dim_lon,dim_lat/), id_geolon)
83 call netcdf_err(error, 'define var geolon for file='//trim(outfile) )
84 error = nf90_put_att(ncid, id_geolon, "long_name", "Longitude")
85 call netcdf_err(error, 'define geolon name for file='//trim(outfile) )
86 error = nf90_put_att(ncid, id_geolon, "units", "degrees_east")
87 call netcdf_err(error, 'define geolon units for file='//trim(outfile) )
88!---geolat
89 error = nf90_def_var(ncid, 'geolat', nf90_float, (/dim_lon,dim_lat/), id_geolat)
90 call netcdf_err(error, 'define var geolat for file='//trim(outfile) )
91 error = nf90_put_att(ncid, id_geolat, "long_name", "Latitude")
92 call netcdf_err(error, 'define geolat name for file='//trim(outfile) )
93 error = nf90_put_att(ncid, id_geolat, "units", "degrees_north")
94 call netcdf_err(error, 'define geolat units for file='//trim(outfile) )
95!---slmsk
96 error = nf90_def_var(ncid, 'slmsk', nf90_float,(/dim_lon,dim_lat/), id_slmsk)
97 call netcdf_err(error, 'define var slmsk for file='//trim(outfile) )
98 error = nf90_put_att(ncid, id_slmsk, "coordinates", "geolon geolat")
99 call netcdf_err(error, 'define slmsk coordinates for file='//trim(outfile) )
100!--- land_frac
101 error = nf90_def_var(ncid, 'land_frac', nf90_float, (/dim_lon,dim_lat/), id_land_frac)
102 call netcdf_err(error, 'define var land_frac for file='//trim(outfile) )
103 error = nf90_put_att(ncid, id_land_frac, "coordinates", "geolon geolat")
104 call netcdf_err(error, 'define land_frac coordinates for file='//trim(outfile) )
105!---orography - raw
106 error = nf90_def_var(ncid, 'orog_raw', nf90_float, (/dim_lon,dim_lat/), id_orog_raw)
107 call netcdf_err(error, 'define var orog_raw for file='//trim(outfile) )
108 error = nf90_put_att(ncid, id_orog_raw, "coordinates", "geolon geolat")
109 call netcdf_err(error, 'define orog_raw coordinates for file='//trim(outfile) )
110!---orography - filtered
111 error = nf90_def_var(ncid, 'orog_filt', nf90_float, (/dim_lon,dim_lat/), id_orog_filt)
112 call netcdf_err(error, 'define var orog_filt for file='//trim(outfile) )
113 error = nf90_put_att(ncid, id_orog_filt, "coordinates", "geolon geolat")
114 call netcdf_err(error, 'define orog_filt coordinates for file='//trim(outfile) )
115!---stddev
116 error = nf90_def_var(ncid, 'stddev', nf90_float, (/dim_lon,dim_lat/), id_stddev)
117 call netcdf_err(error, 'define var stddev for file='//trim(outfile) )
118 error = nf90_put_att(ncid, id_stddev, "coordinates", "geolon geolat")
119 call netcdf_err(error, 'define stddev coordinates for file='//trim(outfile) )
120!---convexity
121 error = nf90_def_var(ncid, 'convexity', nf90_float, (/dim_lon,dim_lat/), id_convex)
122 call netcdf_err(error, 'define var convexity for file='//trim(outfile) )
123 error = nf90_put_att(ncid, id_convex, "coordinates", "geolon geolat")
124 call netcdf_err(error, 'define convexity coordinates for file='//trim(outfile) )
125!---oa1 -> oa4
126 error = nf90_def_var(ncid, 'oa1', nf90_float, (/dim_lon,dim_lat/), id_oa1)
127 call netcdf_err(error, 'define var oa1 for file='//trim(outfile) )
128 error = nf90_put_att(ncid, id_oa1, "coordinates", "geolon geolat")
129 call netcdf_err(error, 'define oa1 coordinates for file='//trim(outfile) )
130 error = nf90_def_var(ncid, 'oa2', nf90_float, (/dim_lon,dim_lat/), id_oa2)
131 call netcdf_err(error, 'define var oa2 for file='//trim(outfile) )
132 error = nf90_put_att(ncid, id_oa2, "coordinates", "geolon geolat")
133 call netcdf_err(error, 'define oa2 coordinates for file='//trim(outfile) )
134 error = nf90_def_var(ncid, 'oa3', nf90_float, (/dim_lon,dim_lat/), id_oa3)
135 call netcdf_err(error, 'define var oa3 for file='//trim(outfile) )
136 error = nf90_put_att(ncid, id_oa3, "coordinates", "geolon geolat")
137 call netcdf_err(error, 'define oa3 coordinates for file='//trim(outfile) )
138 error = nf90_def_var(ncid, 'oa4', nf90_float, (/dim_lon,dim_lat/), id_oa4)
139 call netcdf_err(error, 'define var oa4 for file='//trim(outfile) )
140 error = nf90_put_att(ncid, id_oa4, "coordinates", "geolon geolat")
141 call netcdf_err(error, 'define oa4 coordinates for file='//trim(outfile) )
142!---ol1 -> ol4
143 error = nf90_def_var(ncid, 'ol1', nf90_float, (/dim_lon,dim_lat/), id_ol1)
144 call netcdf_err(error, 'define var ol1 for file='//trim(outfile) )
145 error = nf90_put_att(ncid, id_ol1, "coordinates", "geolon geolat")
146 call netcdf_err(error, 'define ol1 coordinates for file='//trim(outfile) )
147 error = nf90_def_var(ncid, 'ol2', nf90_float, (/dim_lon,dim_lat/), id_ol2)
148 call netcdf_err(error, 'define var ol2 for file='//trim(outfile) )
149 error = nf90_put_att(ncid, id_ol2, "coordinates", "geolon geolat")
150 call netcdf_err(error, 'define ol2 coordinates for file='//trim(outfile) )
151 error = nf90_def_var(ncid, 'ol3', nf90_float, (/dim_lon,dim_lat/), id_ol3)
152 call netcdf_err(error, 'define var ol3 for file='//trim(outfile) )
153 error = nf90_put_att(ncid, id_ol3, "coordinates", "geolon geolat")
154 call netcdf_err(error, 'define ol3 coordinates for file='//trim(outfile) )
155 error = nf90_def_var(ncid, 'ol4', nf90_float, (/dim_lon,dim_lat/), id_ol4)
156 call netcdf_err(error, 'define var ol4 for file='//trim(outfile) )
157 error = nf90_put_att(ncid, id_ol4, "coordinates", "geolon geolat")
158 call netcdf_err(error, 'define ol4 coordinates for file='//trim(outfile) )
159!---theta gamma sigma elvmax
160 error = nf90_def_var(ncid, 'theta', nf90_float, (/dim_lon,dim_lat/), id_theta)
161 call netcdf_err(error, 'define var theta for file='//trim(outfile) )
162 error = nf90_put_att(ncid, id_theta, "coordinates", "geolon geolat")
163 call netcdf_err(error, 'define theta coordinates for file='//trim(outfile) )
164 error = nf90_def_var(ncid, 'gamma', nf90_float, (/dim_lon,dim_lat/), id_gamma)
165 call netcdf_err(error, 'define var gamma for file='//trim(outfile) )
166 error = nf90_put_att(ncid, id_gamma, "coordinates", "geolon geolat")
167 call netcdf_err(error, 'define gamma coordinates for file='//trim(outfile) )
168 error = nf90_def_var(ncid, 'sigma', nf90_float, (/dim_lon,dim_lat/), id_sigma)
169 call netcdf_err(error, 'define var sigma for file='//trim(outfile) )
170 error = nf90_put_att(ncid, id_sigma, "coordinates", "geolon geolat")
171 call netcdf_err(error, 'define sigma coordinates for file='//trim(outfile) )
172 error = nf90_def_var(ncid, 'elvmax', nf90_float, (/dim_lon,dim_lat/), id_elvmax)
173 call netcdf_err(error, 'define var elvmax for file='//trim(outfile) )
174 error = nf90_put_att(ncid, id_elvmax, "coordinates", "geolon geolat")
175 call netcdf_err(error, 'define elvmax coordinates for file='//trim(outfile) )
176
177 error = nf90_enddef(ncid, header_buffer_val,4,0,4)
178 call netcdf_err(error, 'end meta define for file='//trim(outfile) )
179
180 !--- write out data
181 error = nf90_put_var( ncid, id_geolon, geolon(:dim1,:dim2))
182 call netcdf_err(error, 'write var geolon for file='//trim(outfile) )
183 error = nf90_put_var( ncid, id_geolat, geolat(:dim1,:dim2))
184 call netcdf_err(error, 'write var geolat for file='//trim(outfile) )
185
186 error = nf90_put_var( ncid, id_slmsk, slm(:dim1,:dim2))
187 call netcdf_err(error, 'write var slmsk for file='//trim(outfile) )
188 error = nf90_put_var( ncid, id_land_frac, land_frac(:dim1,:dim2))
189 call netcdf_err(error, 'write var land_frac for file='//trim(outfile) )
190
191 error = nf90_put_var( ncid, id_orog_raw, oro(:dim1,:dim2))
192 call netcdf_err(error, 'write var orog_raw for file='//trim(outfile) )
193! We no longer filter the orog, so the raw and filtered records are the same.
194 error = nf90_put_var( ncid, id_orog_filt, oro(:dim1,:dim2))
195 call netcdf_err(error, 'write var orog_filt for file='//trim(outfile) )
196
197 error = nf90_put_var( ncid, id_stddev, hprime(:dim1,:dim2,1))
198 call netcdf_err(error, 'write var stddev for file='//trim(outfile) )
199 error = nf90_put_var( ncid, id_convex, hprime(:dim1,:dim2,2))
200 call netcdf_err(error, 'write var convex for file='//trim(outfile) )
201
202 error = nf90_put_var( ncid, id_oa1, hprime(:dim1,:dim2,3))
203 call netcdf_err(error, 'write var oa1 for file='//trim(outfile) )
204 error = nf90_put_var( ncid, id_oa2, hprime(:dim1,:dim2,4))
205 call netcdf_err(error, 'write var oa2 for file='//trim(outfile) )
206 error = nf90_put_var( ncid, id_oa3, hprime(:dim1,:dim2,5))
207 call netcdf_err(error, 'write var oa3 for file='//trim(outfile) )
208 error = nf90_put_var( ncid, id_oa4, hprime(:dim1,:dim2,6))
209 call netcdf_err(error, 'write var oa4 for file='//trim(outfile) )
210
211 error = nf90_put_var( ncid, id_ol1, hprime(:dim1,:dim2,7))
212 call netcdf_err(error, 'write var ol1 for file='//trim(outfile) )
213 error = nf90_put_var( ncid, id_ol2, hprime(:dim1,:dim2,8))
214 call netcdf_err(error, 'write var ol2 for file='//trim(outfile) )
215 error = nf90_put_var( ncid, id_ol3, hprime(:dim1,:dim2,9))
216 call netcdf_err(error, 'write var ol3 for file='//trim(outfile) )
217 error = nf90_put_var( ncid, id_ol4, hprime(:dim1,:dim2,10))
218 call netcdf_err(error, 'write var ol4 for file='//trim(outfile) )
219
220 error = nf90_put_var( ncid, id_theta, hprime(:dim1,:dim2,11))
221 call netcdf_err(error, 'write var theta for file='//trim(outfile) )
222 error = nf90_put_var( ncid, id_gamma, hprime(:dim1,:dim2,12))
223 call netcdf_err(error, 'write var gamma for file='//trim(outfile) )
224 error = nf90_put_var( ncid, id_sigma, hprime(:dim1,:dim2,13))
225 call netcdf_err(error, 'write var sigma for file='//trim(outfile) )
226 error = nf90_put_var( ncid, id_elvmax, hprime(:dim1,:dim2,14))
227 call netcdf_err(error, 'write var elvmax for file='//trim(outfile) )
228
229 error = nf90_close(ncid)
230 call netcdf_err(error, 'close file='//trim(outfile) )
231
232 end subroutine write_netcdf
233
239 subroutine netcdf_err( err, string )
240 use netcdf
241 implicit none
242 integer, intent(in) :: err
243 character(len=*), intent(in) :: string
244 character(len=256) :: errmsg
245
246 if( err.EQ.nf90_noerr )return
247 errmsg = nf90_strerror(err)
248 print*, 'FATAL ERROR: ', trim(string), ': ', trim(errmsg)
249 call abort
250
251 return
252 end subroutine netcdf_err
253
265
266 subroutine write_mask_netcdf(im, jm, slm, land_frac, ntiles, tile, geolon, geolat)
267 use netcdf
268 implicit none
269 integer, intent(in):: im, jm, ntiles, tile
270 real, intent(in), dimension(im,jm) :: slm, geolon, geolat, land_frac
271 character(len=128) :: outfile
272 integer :: error, ncid
273 integer :: header_buffer_val = 16384
274 integer :: fsize=65536, inital = 0
275 integer :: dim1, dim2
276 integer :: dim_lon, dim_lat
277 integer :: id_geolon,id_geolat
278 integer :: id_slmsk,id_land_frac
279
280 if(ntiles > 1) then
281 write(outfile, '(a,i4.4,a)') 'out.oro.tile', tile, '.nc'
282 else
283 outfile = "out.oro.nc"
284 endif
285
286 dim1=im
287 dim2=jm
288
289 !--- open the file
290 error = nf90_create(outfile, ior(nf90_netcdf4,nf90_classic_model), ncid, &
291 initialsize=inital, chunksize=fsize)
292 call netcdf_err(error, 'Creating file '//trim(outfile) )
293 !--- define dimension
294 error = nf90_def_dim(ncid, 'lon', dim1, dim_lon)
295 call netcdf_err(error, 'define dimension lon for file='//trim(outfile) )
296 error = nf90_def_dim(ncid, 'lat', dim2, dim_lat)
297 call netcdf_err(error, 'define dimension lat for file='//trim(outfile) )
298
299 !--- define field
300!---geolon
301 error = nf90_def_var(ncid, 'geolon', nf90_float, (/dim_lon,dim_lat/), id_geolon)
302 call netcdf_err(error, 'define var geolon for file='//trim(outfile) )
303 error = nf90_put_att(ncid, id_geolon, "long_name", "Longitude")
304 call netcdf_err(error, 'define geolon name for file='//trim(outfile) )
305 error = nf90_put_att(ncid, id_geolon, "units", "degrees_east")
306 call netcdf_err(error, 'define geolon units for file='//trim(outfile) )
307!---geolat
308 error = nf90_def_var(ncid, 'geolat', nf90_float, (/dim_lon,dim_lat/), id_geolat)
309 call netcdf_err(error, 'define var geolat for file='//trim(outfile) )
310 error = nf90_put_att(ncid, id_geolat, "long_name", "Latitude")
311 call netcdf_err(error, 'define geolat name for file='//trim(outfile) )
312 error = nf90_put_att(ncid, id_geolat, "units", "degrees_north")
313 call netcdf_err(error, 'define geolat units for file='//trim(outfile) )
314!---slmsk
315 error = nf90_def_var(ncid, 'slmsk', nf90_float, (/dim_lon,dim_lat/), id_slmsk)
316 call netcdf_err(error, 'define var slmsk for file='//trim(outfile) )
317 error = nf90_put_att(ncid, id_slmsk, "coordinates", "geolon geolat")
318 call netcdf_err(error, 'define slmsk coordinates for file='//trim(outfile) )
319!--- land_frac
320 error = nf90_def_var(ncid, 'land_frac', nf90_float, (/dim_lon,dim_lat/), id_land_frac)
321 call netcdf_err(error, 'define var land_frac for file='//trim(outfile) )
322 error = nf90_put_att(ncid, id_land_frac, "coordinates", "geolon geolat")
323 call netcdf_err(error, 'define land_frac coordinates for file='//trim(outfile) )
324
325 error = nf90_enddef(ncid, header_buffer_val,4,0,4)
326 call netcdf_err(error, 'end meta define for file='//trim(outfile) )
327
328 !--- write out data
329 error = nf90_put_var( ncid, id_geolon, geolon(:dim1,:dim2))
330 call netcdf_err(error, 'write var geolon for file='//trim(outfile) )
331
332 error = nf90_put_var( ncid, id_geolat, geolat(:dim1,:dim2))
333 call netcdf_err(error, 'write var geolat for file='//trim(outfile) )
334
335 error = nf90_put_var( ncid, id_slmsk, slm(:dim1,:dim2))
336 call netcdf_err(error, 'write var slmsk for file='//trim(outfile) )
337
338 error = nf90_put_var( ncid, id_land_frac, land_frac(:dim1,:dim2))
339 call netcdf_err(error, 'write var land_frac for file='//trim(outfile) )
340
341 error = nf90_close(ncid)
342 call netcdf_err(error, 'close file='//trim(outfile) )
343
344 end subroutine write_mask_netcdf
345
355
356 subroutine read_mask(merge_file,slm,land_frac,lake_frac,im,jm)
357
358 use netcdf
359
360 implicit none
361
362 character(len=*), intent(in) :: merge_file
363
364 integer, intent(in) :: im, jm
365
366 real, intent(out) :: land_frac(im,jm)
367 real, intent(out) :: lake_frac(im,jm)
368 real, intent(out) :: slm(im,jm)
369
370 integer ncid, error, id_var
371
372 print*,'- READ IN EXTERNAL LANDMASK FILE: ',trim(merge_file)
373 error=nf90_open(merge_file,nf90_nowrite,ncid)
374 call netcdf_err(error, 'Open file '//trim(merge_file) )
375
376 error=nf90_inq_varid(ncid, 'land_frac', id_var)
377 call netcdf_err(error, 'inquire varid of land_frac')
378 error=nf90_get_var(ncid, id_var, land_frac)
379 call netcdf_err(error, 'inquire data of land_frac')
380
381 error=nf90_inq_varid(ncid, 'slmsk', id_var)
382 call netcdf_err(error, 'inquire varid of slmsk')
383 error=nf90_get_var(ncid, id_var, slm)
384 call netcdf_err(error, 'inquire data of slmsk')
385
386 error=nf90_inq_varid(ncid, 'lake_frac', id_var)
387 call netcdf_err(error, 'inquire varid of lake_frac')
388 error=nf90_get_var(ncid, id_var, lake_frac)
389 call netcdf_err(error, 'inquire data of lake_frac')
390
391 error = nf90_close(ncid)
392
393 end subroutine read_mask
394
401 subroutine read_mdl_dims(mdl_grid_file, im, jm)
402
403 use netcdf
404
405 implicit none
406
407 character(len=*), intent(in) :: mdl_grid_file
408
409 integer, intent(out) :: im, jm
410
411 integer ncid, error, id_dim, nx, ny
412
413 print*, "- READ MDL GRID DIMENSIONS FROM= ", trim(mdl_grid_file)
414
415 error=nf90_open(mdl_grid_file, nf90_nowrite, ncid)
416 call netcdf_err(error, 'Opening file '//trim(mdl_grid_file) )
417
418 error=nf90_inq_dimid(ncid, 'nx', id_dim)
419 call netcdf_err(error, 'inquire dimension nx from file '// trim(mdl_grid_file) )
420 error=nf90_inquire_dimension(ncid, id_dim, len=nx)
421 call netcdf_err(error, 'inquire nx from file '//trim(mdl_grid_file) )
422
423 error=nf90_inq_dimid(ncid, 'ny', id_dim)
424 call netcdf_err(error, 'inquire dimension ny from file '// trim(mdl_grid_file) )
425 error=nf90_inquire_dimension(ncid, id_dim, len=ny)
426 call netcdf_err(error, 'inquire ny from file '//trim(mdl_grid_file) )
427
428 error=nf90_close(ncid)
429
430 im = nx/2
431 jm = ny/2
432
433 print*,"- MDL GRID DIMENSIONS ", im, jm
434
435 end subroutine read_mdl_dims
436
451 subroutine read_mdl_grid_file(mdl_grid_file, im, jm, &
452 geolon, geolon_c, geolat, geolat_c, dx, dy, &
453 is_north_pole, is_south_pole)
454
455 use netcdf
456
457 use orog_utils, only : find_poles, find_nearest_pole_points
458
459 implicit none
460
461 character(len=*), intent(in) :: mdl_grid_file
462
463 integer, intent(in) :: im, jm
464
465 logical, intent(out) :: is_north_pole(im,jm)
466 logical, intent(out) :: is_south_pole(im,jm)
467
468 real, intent(out) :: geolat(im,jm)
469 real, intent(out) :: geolat_c(im+1,jm+1)
470 real, intent(out) :: geolon(im,jm)
471 real, intent(out) :: geolon_c(im+1,jm+1)
472 real, intent(out) :: dx(im,jm), dy(im,jm)
473
474 integer :: i, j
475 integer :: ncid, error, id_var, nx, ny
476 integer :: i_south_pole,j_south_pole
477 integer :: i_north_pole,j_north_pole
478
479 real, allocatable :: tmpvar(:,:)
480
481 nx = 2*im
482 ny = 2*jm
483
484 allocate(tmpvar(nx+1,ny+1))
485
486 print*, "- OPEN AND READ= ", trim(mdl_grid_file)
487
488 error=nf90_open(mdl_grid_file, nf90_nowrite, ncid)
489 call netcdf_err(error, 'Opening file '//trim(mdl_grid_file) )
490
491 error=nf90_inq_varid(ncid, 'x', id_var)
492 call netcdf_err(error, 'inquire varid of x from file ' // trim(mdl_grid_file))
493 error=nf90_get_var(ncid, id_var, tmpvar)
494 call netcdf_err(error, 'inquire data of x from file ' // trim(mdl_grid_file))
495
496! Adjust lontitude to be between 0 and 360.
497 do j = 1,ny+1
498 do i = 1,nx+1
499 if(tmpvar(i,j) .GT. 360) tmpvar(i,j) = tmpvar(i,j) - 360
500 if(tmpvar(i,j) .LT. 0) tmpvar(i,j) = tmpvar(i,j) + 360
501 enddo
502 enddo
503
504 geolon(1:im,1:jm) = tmpvar(2:nx:2,2:ny:2)
505 geolon_c(1:im+1,1:jm+1) = tmpvar(1:nx+1:2,1:ny+1:2)
506
507 error=nf90_inq_varid(ncid, 'y', id_var)
508 call netcdf_err(error, 'inquire varid of y from file ' // trim(mdl_grid_file))
509 error=nf90_get_var(ncid, id_var, tmpvar)
510 call netcdf_err(error, 'inquire data of y from file ' // trim(mdl_grid_file))
511
512 geolat(1:im,1:jm) = tmpvar(2:nx:2,2:ny:2)
513 geolat_c(1:im+1,1:jm+1) = tmpvar(1:nx+1:2,1:ny+1:2)
514
515 call find_poles(tmpvar, nx, ny, i_north_pole, j_north_pole, &
516 i_south_pole, j_south_pole)
517
518 deallocate(tmpvar)
519
520 call find_nearest_pole_points(i_north_pole, j_north_pole, &
521 i_south_pole, j_south_pole, im, jm, is_north_pole, &
522 is_south_pole)
523
524 allocate(tmpvar(nx,ny))
525
526 error=nf90_inq_varid(ncid, 'area', id_var)
527 call netcdf_err(error, 'inquire varid of area from file ' // trim(mdl_grid_file))
528 error=nf90_get_var(ncid, id_var, tmpvar)
529 call netcdf_err(error, 'inquire data of area from file ' // trim(mdl_grid_file))
530
531 error = nf90_close(ncid)
532
533 do j = 1, jm
534 do i = 1, im
535 dx(i,j) = sqrt(tmpvar(2*i-1,2*j-1)+tmpvar(2*i,2*j-1) &
536 + tmpvar(2*i-1,2*j )+tmpvar(2*i,2*j ))
537 dy(i,j) = dx(i,j)
538 enddo
539 enddo
540
541 deallocate(tmpvar)
542
543 end subroutine read_mdl_grid_file
544
551 subroutine read_global_orog(imn,jmn,glob)
552
553 use orog_utils, only : transpose_orog
554 use netcdf
555
556 implicit none
557
558 integer, intent(in) :: imn, jmn
559 integer*2, intent(out) :: glob(imn,jmn)
560
561 integer :: ncid, error, id_dim, id_var, idim, jdim
562
563 print*,"- OPEN AND READ ./topography.gmted2010.30s.nc"
564
565 error=nf90_open("./topography.gmted2010.30s.nc", &
566 nf90_nowrite, ncid)
567 call netcdf_err(error, 'Opening file topography.gmted2010.30s.nc' )
568
569 error=nf90_inq_dimid(ncid, 'idim', id_dim)
570 call netcdf_err(error, 'Inquire dimid of idim' )
571
572 error=nf90_inquire_dimension(ncid,id_dim,len=idim)
573 call netcdf_err(error, 'Reading idim' )
574
575 if (imn /= idim) then
576 print*,"FATAL ERROR: i-dimensions do not match."
577 endif
578
579 error=nf90_inq_dimid(ncid, 'jdim', id_dim)
580 call netcdf_err(error, 'Inquire dimid of jdim' )
581
582 error=nf90_inquire_dimension(ncid,id_dim,len=jdim)
583 call netcdf_err(error, 'Reading jdim' )
584
585 if (jmn /= jdim) then
586 print*,"FATAL ERROR: j-dimensions do not match."
587 endif
588
589 error=nf90_inq_varid(ncid, 'topo', id_var)
590 call netcdf_err(error, 'Inquire varid of topo')
591
592 error=nf90_get_var(ncid, id_var, glob)
593 call netcdf_err(error, 'Reading topo')
594
595 error = nf90_close(ncid)
596
597 print*,"- MAX/MIN OF OROGRAPHY DATA ",maxval(glob),minval(glob)
598
599 call transpose_orog(imn,jmn,glob)
600
601 return
602 end subroutine read_global_orog
603
610 subroutine read_global_mask(imn, jmn, mask)
611
612 use orog_utils, only : transpose_mask
613 use netcdf
614
615 implicit none
616
617 integer, intent(in) :: imn, jmn
618
619 integer(1), intent(out) :: mask(imn,jmn)
620
621 integer :: ncid, id_var, id_dim, error, idim, jdim
622
623 print*,"- OPEN AND READ ./landcover.umd.30s.nc"
624
625 error=nf90_open("./landcover.umd.30s.nc",nf90_nowrite,ncid)
626 call netcdf_err(error, 'Opening file landcover.umd.30s.nc' )
627
628 error=nf90_inq_dimid(ncid, 'idim', id_dim)
629 call netcdf_err(error, 'Inquire dimid of idim' )
630
631 error=nf90_inquire_dimension(ncid,id_dim,len=idim)
632 call netcdf_err(error, 'Reading idim' )
633
634 if (imn /= idim) then
635 print*,"FATAL ERROR: i-dimensions do not match."
636 endif
637
638 error=nf90_inq_dimid(ncid, 'jdim', id_dim)
639 call netcdf_err(error, 'Inquire dimid of jdim' )
640
641 error=nf90_inquire_dimension(ncid,id_dim,len=jdim)
642 call netcdf_err(error, 'Reading jdim' )
643
644 if (jmn /= jdim) then
645 print*,"FATAL ERROR: j-dimensions do not match."
646 endif
647
648 error=nf90_inq_varid(ncid, 'land_mask', id_var)
649 call netcdf_err(error, 'Inquire varid of land_mask')
650
651 error=nf90_get_var(ncid, id_var, mask)
652 call netcdf_err(error, 'Inquire data of land_mask')
653
654 error = nf90_close(ncid)
655
656 call transpose_mask(imn,jmn,mask)
657
658 end subroutine read_global_mask
659
668 subroutine qc_orog_by_ramp(imn, jmn, zavg, zslm)
669
670 use netcdf
671
672 implicit none
673
674 integer, intent(in) :: imn, jmn
675 integer, intent(inout) :: zavg(imn,jmn)
676 integer, intent(inout) :: zslm(imn,jmn)
677
678 integer :: i, j, error, ncid, id_var, id_dim, jramp
679
680 real(4), allocatable :: gice(:,:)
681
682! Read 30-sec Antarctica RAMP data. Points scan from South
683! to North, and from Greenwich to Greenwich.
684
685 print*,"- OPEN/READ RAMP DATA ./topography.antarctica.ramp.30s.nc"
686
687 error=nf90_open("./topography.antarctica.ramp.30s.nc", &
688 nf90_nowrite, ncid)
689 call netcdf_err(error, 'Opening RAMP topo file' )
690
691 error=nf90_inq_dimid(ncid, 'jdim', id_dim)
692 call netcdf_err(error, 'Inquire dimid of jdim' )
693
694 error=nf90_inquire_dimension(ncid, id_dim, len=jramp)
695 call netcdf_err(error, 'Reading jdim' )
696
697 allocate (gice(imn+1,jramp))
698
699 error=nf90_inq_varid(ncid, 'topo', id_var)
700 call netcdf_err(error, 'Inquire varid of RAMP topo')
701
702 error=nf90_get_var(ncid, id_var, gice)
703 call netcdf_err(error, 'Inquire data of RAMP topo')
704
705 error = nf90_close(ncid)
706
707 print*,"- QC GLOBAL OROGRAPHY DATA WITH RAMP."
708
709! If RAMP values are valid, replace the global value with the RAMP
710! value. Invalid values are less than or equal to 0 (0, -1, or -99).
711
712 do j = 1, jramp
713 do i = 1, imn
714 if ( gice(i,j) .gt. 0.) then
715 zavg(i,j) = int( gice(i,j) + 0.5 )
716 zslm(i,j) = 1
717 endif
718 enddo
719 enddo
720
721 deallocate (gice)
722
723 end subroutine qc_orog_by_ramp
724
725 end module io_utils
io_utils
Module containing utilities that read and write data.
Definition io_utils.F90:9
io_utils::read_global_mask
subroutine, public read_global_mask(imn, jmn, mask)
Read input global 30-arc second land mask data.
Definition io_utils.F90:611
io_utils::netcdf_err
subroutine netcdf_err(err, string)
Check NetCDF error code and output the error message.
Definition io_utils.F90:240
io_utils::read_global_orog
subroutine, public read_global_orog(imn, jmn, glob)
Read input global 30-arc second orography data.
Definition io_utils.F90:552
io_utils::read_mdl_grid_file
subroutine, public read_mdl_grid_file(mdl_grid_file, im, jm, geolon, geolon_c, geolat, geolat_c, dx, dy, is_north_pole, is_south_pole)
Read the grid dimensions from the model 'grid' file.
Definition io_utils.F90:454
io_utils::qc_orog_by_ramp
subroutine, public qc_orog_by_ramp(imn, jmn, zavg, zslm)
Quality control the global orography and landmask data over Antarctica using RAMP data.
Definition io_utils.F90:669
io_utils::read_mdl_dims
subroutine, public read_mdl_dims(mdl_grid_file, im, jm)
Read the grid dimensions from the model 'grid' file.
Definition io_utils.F90:402
io_utils::read_mask
subroutine, public read_mask(merge_file, slm, land_frac, lake_frac, im, jm)
Read the land mask file.
Definition io_utils.F90:357
io_utils::write_mask_netcdf
subroutine, public write_mask_netcdf(im, jm, slm, land_frac, ntiles, tile, geolon, geolat)
Write the land mask file.
Definition io_utils.F90:267
io_utils::write_netcdf
subroutine, public write_netcdf(im, jm, slm, land_frac, oro, hprime, ntiles, tile, geolon, geolat, lon, lat)
Write out orography file in netcdf format.
Definition io_utils.F90:42
orog_utils
Module containing utilites used by the orog program.
Definition orog_utils.F90:8
orog_utils::transpose_orog
subroutine, public transpose_orog(imn, jmn, glob)
Transpose the global orography data by flipping the poles and moving the starting longitude to Greenw...
Definition orog_utils.F90:220
orog_utils::find_poles
subroutine, public find_poles(geolat, nx, ny, i_north_pole, j_north_pole, i_south_pole, j_south_pole)
Find the point on the model grid tile closest to the north and south pole.
Definition orog_utils.F90:411
orog_utils::find_nearest_pole_points
subroutine, public find_nearest_pole_points(i_north_pole, j_north_pole, i_south_pole, j_south_pole, im, jm, is_north_pole, is_south_pole)
Find the point on the model grid tile closest to the north and south pole.
Definition orog_utils.F90:487
orog_utils::transpose_mask
subroutine, public transpose_mask(imn, jmn, mask)
Transpose the global landmask by flipping the poles and moving the starting longitude to Greenwich.
Definition orog_utils.F90:177