orog_mask_tools  1.13.0
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
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
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::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_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
Module containing utilities that read and write data.
Definition: io_utils.F90:9
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
lake_frac
program lake_frac
This program computes lake fraction and depth numbers for FV3 cubed sphere grid cells, from a high resolution lat/lon data set.
Definition: lakefrac.F90:21
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
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
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
orog_utils
Module containing utilites used by the orog program.
Definition: orog_utils.F90:8
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::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
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