sfc_climo_gen  1.13.0
 All Data Structures Files Functions Variables Pages
source_grid.F90
Go to the documentation of this file.
1 
5 
12 module source_grid
13 
14  use esmf
15  use utils
16 
17  implicit none
18 
19  private
20 
21  character(len=50), allocatable, public :: field_names(:)
22  character(len=75), public :: source
23 
24  integer, public :: i_src
25  integer, public :: j_src
26  integer, public :: num_records
27  integer, public :: num_time_recs
28  integer, public :: num_fields
31  integer, allocatable, public :: day_of_rec(:)
33 
34  type(esmf_grid), public :: grid_src
35 
36  public :: define_source_grid
37  public :: source_grid_cleanup
38 
39  contains
40 
51  subroutine define_source_grid(localpet, npets, input_file)
52 
53  use mpi
54  use netcdf
55 
56  implicit none
57 
58  character(len=*), intent(in) :: input_file
59 
60  integer, intent(in) :: localpet, npets
61 
62  character(len=50) :: field_names_save(100)
63  character(len=50) :: vname
64 
65  integer :: dimid, dims(1), ncid, status
66  integer :: count, num_vars, n
67  integer :: rc, varid, i, j, i_src_corner
68  integer(esmf_kind_i4), pointer :: mask_ptr(:,:)
69  integer :: clb(2), cub(2)
70  integer :: clb_corner(2), cub_corner(2)
71 
72  real(esmf_kind_r4), allocatable :: mask_global(:,:)
73  real(esmf_kind_r8), allocatable :: lat_global(:)
74  real(esmf_kind_r8), allocatable :: lon_global(:)
75  real(esmf_kind_r8), allocatable :: lat_corner_global(:)
76  real(esmf_kind_r8), allocatable :: lon_corner_global(:)
77  real(esmf_kind_r4), pointer :: mask_field_ptr(:,:)
78  real(esmf_kind_r8), pointer :: lat_ptr(:,:)
79  real(esmf_kind_r8), pointer :: lon_ptr(:,:)
80  real(esmf_kind_r8), pointer :: lat_corner_ptr(:,:)
81  real(esmf_kind_r8), pointer :: lon_corner_ptr(:,:)
82  real :: lon_extent
83  real(esmf_kind_r4) :: missing
84 
85  type(esmf_field) :: mask_field
86  type(esmf_polekind_flag) :: polekindflag(2)
87 
88  print*,"- OPEN FILE: ", trim(input_file)
89  status = nf90_open(input_file, nf90_nowrite, ncid)
90  call netcdf_err(status, "OPENING INPUT SOURCE FILE")
91 
92  status = nf90_get_att(ncid, nf90_global, 'source', source)
93  if (status /= nf90_noerr) source="unknown"
94 
95  if(localpet == 0) print*,'- SOURCE OF DATA IS: ', trim(source)
96 
97  status = nf90_inq_dimid(ncid, 'idim', dimid)
98  call netcdf_err(status, "READING I DIMENSION ID OF SOURCE DATA")
99  status = nf90_inquire_dimension(ncid, dimid, len=i_src)
100  call netcdf_err(status, "READING I DIMENSION OF SOURCE DATA")
101 
102  status = nf90_inq_dimid(ncid, 'jdim', dimid)
103  call netcdf_err(status, "READING J DIMENSION ID OF SOURCE DATA")
104  status = nf90_inquire_dimension(ncid, dimid, len=j_src)
105  call netcdf_err(status, "READING J DIMENSION OF SOURCE DATA")
106 
107  if(localpet == 0) print*,'- I/J DIMENSION OF SOURCE DATA: ', i_src, j_src
108 
109  allocate(lat_global(j_src))
110  status = nf90_inq_varid(ncid, 'lat', varid)
111  call netcdf_err(status, "READING SOURCE DATA LATITUDE ID")
112  status = nf90_get_var(ncid, varid, lat_global)
113  call netcdf_err(status, "READING SOURCE DATA LATITUDES")
114 
115  allocate(lon_global(i_src))
116  status = nf90_inq_varid(ncid, 'lon', varid)
117  call netcdf_err(status, "READING SOURCE DATA LONGITUDE ID")
118  status = nf90_get_var(ncid, varid, lon_global)
119  call netcdf_err(status, "READING SOURCE DATA LONGITUDE")
120 
121  allocate(lat_corner_global(j_src+1))
122  status = nf90_inq_varid(ncid, 'lat_corner', varid)
123  call netcdf_err(status, "READING SOURCE DATA CORNER LATITUDE ID")
124  status = nf90_get_var(ncid, varid, lat_corner_global)
125  call netcdf_err(status, "READING SOURCE DATA CORNER LATITUDE")
126 
127 !-----------------------------------------------------------------------
128 ! Dimension of lon_corner depends on whether input data is periodic
129 ! (global) or regional.
130 !-----------------------------------------------------------------------
131 
132  status = nf90_inq_varid(ncid, 'lon_corner', varid)
133  call netcdf_err(status, "READING SOURCE DATA CORNER LONGITUDE ID")
134  status = nf90_inquire_variable(ncid, varid, dimids=dims)
135  call netcdf_err(status, "READING SOURCE DATA CORNER LONGITUDE DIMIDS")
136  status = nf90_inquire_dimension(ncid, dims(1), len=i_src_corner)
137  call netcdf_err(status, "READING SOURCE DATA CORNER LONGITUDE DIMS")
138  allocate(lon_corner_global(i_src_corner))
139  status = nf90_get_var(ncid, varid, lon_corner_global)
140  call netcdf_err(status, "READING SOURCE DATA CORNER LONGITUDE")
141 
142  status = nf90_inq_dimid(ncid, 'time', dimid)
143  call netcdf_err(status, "READING SOURCE DATA TIME ID")
144  status = nf90_inquire_dimension(ncid, dimid, len=num_time_recs)
145  call netcdf_err(status, "READING SOURCE DATA NUM TIME PERIODS")
146 
147  allocate(day_of_rec(num_time_recs))
148  status = nf90_inq_varid(ncid, 'time', varid)
149  call netcdf_err(status, "READING SOURCE DATA TIME VARID")
150  status = nf90_get_var(ncid, varid, day_of_rec)
151  call netcdf_err(status, "READING SOURCE DATA DAY OF RECORD")
152 
153  print*,'- SOURCE DATA DAYS OF RECORD: ', day_of_rec
154 
155  status = nf90_inquire(ncid, nvariables=num_vars)
156  call netcdf_err(status, "READING NUMBER OF VARIABLES SOURCE DATA")
157 
158 !-----------------------------------------------------------------------
159 ! Assumes files contain records for time, lat, lon, lat_corner,
160 ! lon_corner. So number of fields processed will be the total
161 ! number of variables minus 5.
162 !-----------------------------------------------------------------------
163 
164 ! NOTE: the new BNU soil type data contains extra records for
165 ! sand and clay percentages. These extra records are not need yet,
166 ! so add logic to temporarily ignore them.
167 
168  count = 0
169  do n = 1, num_vars
170 
171  status = nf90_inquire_variable(ncid, n, name=vname)
172  call netcdf_err(status, "READING FIELD NAMES")
173 
174  if (trim(vname) == 'time') cycle
175  if (trim(vname) == 'lon') cycle
176  if (trim(vname) == 'lon_corner') cycle
177  if (trim(vname) == 'lat') cycle
178  if (trim(vname) == 'lat_corner') cycle
179  if (trim(vname) == 'clay_lev1') cycle
180  if (trim(vname) == 'clay_top') cycle
181  if (trim(vname) == 'sand_lev1') cycle
182  if (trim(vname) == 'sand_top') cycle
183 
184  count = count + 1
185  field_names_save(count) = vname
186 
187  enddo
188 
189  num_fields = count
190  num_records = num_vars * num_time_recs
191 
192  allocate(field_names(num_fields))
193 
194  field_names = field_names_save(1:num_fields)
195 
196  if(localpet==0) print*,'- FIELDS TO BE PROCESSED: ', field_names
197 
198  if (localpet == 0) then
199  allocate(mask_global(i_src,j_src))
200  status = nf90_inq_varid(ncid, field_names(1), varid)
201  call netcdf_err(status, "READING FIELD ID")
202  status = nf90_get_var(ncid, varid, mask_global)
203  call netcdf_err(status, "READING FIELD")
204  else
205  allocate(mask_global(0,0))
206  endif
207 
208 !--------------------------------------------------------------------------
209 ! Read in missing value. This is used to mask out data at non-land
210 ! points.
211 !--------------------------------------------------------------------------
212 
213  status = nf90_inq_varid(ncid, field_names(1), varid)
214  call netcdf_err(status, "READING FIELD 1 ID")
215  status=nf90_get_att(ncid, varid, 'missing_value', missing)
216  call netcdf_err(status, "READING MISSING VALUE")
217 
218  status = nf90_close(ncid)
219 
220 !--------------------------------------------------------------------------
221 ! Create ESMF grid object for the source data grid. Check if
222 ! data is periodic in the east/west direction.
223 !
224 ! Note: When using regional data, there is always the chance of
225 ! the target grid being located outside the input grid. ESMF
226 ! support recommends passing back the dstField (esmf_typekind_i4) from
227 ! all calls to FieldRegridStore and checking for the
228 ! ESMF_REGRIDSTATUS_OUTSIDE flag.
229 !--------------------------------------------------------------------------
230 
231  lon_extent = mod((lon_global(i_src)-lon_global(1))-1+3600,360.)+1.0
232 
233  if (lon_extent < 350.0) then
234 
235  print*,"- CALL GridCreateNoPeriDim FOR SOURCE GRID"
236  grid_src = esmf_gridcreatenoperidim(minindex=(/1,1/), &
237  maxindex=(/i_src,j_src/), &
238  coordsys=esmf_coordsys_sph_deg, &
239  regdecomp=(/1,npets/), &
240  name="source_grid", &
241  indexflag=esmf_index_global, rc=rc)
242  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
243  call error_handler("IN GridCreateNoPeriDim.", rc)
244 
245  else
246 
247  polekindflag(1:2) = esmf_polekind_monopole
248 
249  print*,"- CALL GridCreate1PeriDim FOR SOURCE GRID"
250  grid_src = esmf_gridcreate1peridim(minindex=(/1,1/), &
251  maxindex=(/i_src,j_src/), &
252  polekindflag=polekindflag, &
253  periodicdim=1, &
254  poledim=2, &
255  coordsys=esmf_coordsys_sph_deg, &
256  regdecomp=(/1,npets/), &
257  name="source_grid", &
258  indexflag=esmf_index_global, rc=rc)
259  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
260  call error_handler("IN GridCreate1PeriDim.", rc)
261 
262  endif
263 
264  print*,"- CALL FieldCreate FOR SOURCE GRID LANDMASK."
265  mask_field = esmf_fieldcreate(grid_src, &
266  typekind=esmf_typekind_r4, &
267  staggerloc=esmf_staggerloc_center, &
268  name="source grid land mask", &
269  rc=rc)
270  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
271  call error_handler("IN FieldCreate.", rc)
272 
273  print*,"- CALL FieldScatter FOR SOURCE GRID MASK."
274  call esmf_fieldscatter(mask_field, mask_global, rootpet=0, rc=rc)
275  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
276  call error_handler("IN FieldScatter.", rc)
277 
278  print*,"- CALL GridAddItem FOR SOURCE GRID MASK."
279  call esmf_gridadditem(grid_src, &
280  itemflag=esmf_griditem_mask, &
281  staggerloc=esmf_staggerloc_center, &
282  rc=rc)
283  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
284  call error_handler("IN GridAddItem.", rc)
285 
286  print*,"- CALL GridGetItem FOR SOURCE GRID MASK."
287  nullify(mask_ptr)
288  call esmf_gridgetitem(grid_src, &
289  itemflag=esmf_griditem_mask, &
290  farrayptr=mask_ptr, &
291  totallbound=clb, &
292  totalubound=cub, &
293  rc=rc)
294  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
295  call error_handler("IN GridGetItem", rc)
296 
297  print*,"- CALL FieldGet FOR SOURCE GRID LANDMASK."
298  nullify(mask_field_ptr)
299  call esmf_fieldget(mask_field, &
300  farrayptr=mask_field_ptr, &
301  rc=rc)
302  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
303  call error_handler("IN FieldGet.", rc)
304 
305  do j = clb(2), cub(2)
306  do i = clb(1), cub(1)
307  if ( abs(mask_field_ptr(i,j)-missing) < 0.001) then
308  mask_ptr(i,j) = 0
309  else
310  mask_ptr(i,j) = 1
311  endif
312  enddo
313  enddo
314 
315  deallocate(mask_global)
316 
317  print*,"- CALL FieldDestroy FOR SOURCE GRID LAND MASK."
318  call esmf_fielddestroy(mask_field,rc=rc)
319  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
320  call error_handler("IN FieldDestroy.", rc)
321 
322 ! Set lat/lons of grid points
323 
324  print*,"- CALL GridAddCoord FOR SOURCE GRID CENTER LOCATION."
325  call esmf_gridaddcoord(grid_src, &
326  staggerloc=esmf_staggerloc_center, rc=rc)
327  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
328  call error_handler("IN GridAddCoord.", rc)
329 
330  print*,"- CALL GridGetCoord FOR SOURCE GRID CENTER LONGITUDE."
331  nullify(lon_ptr)
332  call esmf_gridgetcoord(grid_src, &
333  staggerloc=esmf_staggerloc_center, &
334  coorddim=1, &
335  farrayptr=lon_ptr, rc=rc)
336  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
337  call error_handler("IN GridGetCoord.", rc)
338 
339  print*,"- CALL GridGetCoord FOR SOURCE GRID CENTER LATITUDE."
340  nullify(lat_ptr)
341  call esmf_gridgetcoord(grid_src, &
342  staggerloc=esmf_staggerloc_center, &
343  coorddim=2, &
344  farrayptr=lat_ptr, rc=rc)
345  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
346  call error_handler("IN GridGetCoord.", rc)
347 
348  do j = clb(2), cub(2)
349  lat_ptr(:,j) = lat_global(j)
350  enddo
351 
352  do i = clb(1), cub(1)
353  lon_ptr(i,:) = lon_global(i)
354  enddo
355 
356  print*,"- CALL GridAddCoord FOR SOURCE GRID CORNER LOCATION."
357  call esmf_gridaddcoord(grid_src, &
358  staggerloc=esmf_staggerloc_corner, &
359  rc=rc)
360  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
361  call error_handler("IN GridAddCoord.", rc)
362 
363  print*,"- CALL GridGetCoord FOR SOURCE GRID CORNER LONGITUDE."
364  nullify(lon_corner_ptr)
365  call esmf_gridgetcoord(grid_src, &
366  staggerloc=esmf_staggerloc_corner, &
367  coorddim=1, &
368  farrayptr=lon_corner_ptr, &
369  rc=rc)
370  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
371  call error_handler("IN GridGetCoord.", rc)
372 
373  print*,"- CALL GridGetCoord FOR SOURCE GRID CORNER LATITUDE."
374  nullify(lat_corner_ptr)
375  call esmf_gridgetcoord(grid_src, &
376  staggerloc=esmf_staggerloc_corner, &
377  coorddim=2, &
378  computationallbound=clb_corner, &
379  computationalubound=cub_corner, &
380  farrayptr=lat_corner_ptr, &
381  rc=rc)
382  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
383  call error_handler("IN GridGetCoord.", rc)
384 
385  do j = clb_corner(2), cub_corner(2)
386  lat_corner_ptr(:,j) = lat_corner_global(j)
387  enddo
388 
389  do i = clb_corner(1), cub_corner(1)
390  lon_corner_ptr(i,:) = lon_corner_global(i)
391  enddo
392 
393  deallocate(lat_global)
394  deallocate(lon_global)
395  deallocate(lat_corner_global)
396  deallocate(lon_corner_global)
397 
398  end subroutine define_source_grid
399 
403  subroutine source_grid_cleanup
404 
405  implicit none
406 
407  integer :: rc
408 
409  print*,"- CALL GridDestroy FOR SOURCE GRID."
410  call esmf_griddestroy(grid_src,rc=rc)
411 
412  deallocate(day_of_rec)
413  deallocate(field_names)
414 
415  end subroutine source_grid_cleanup
416 
417  end module source_grid
source_grid::source_grid_cleanup
subroutine, public source_grid_cleanup
Free up memory associated with this module.
Definition: source_grid.F90:403
source_grid::define_source_grid
subroutine, public define_source_grid(localpet, npets, input_file)
Defines esmf grid object for source grid.
Definition: source_grid.F90:51
utils::netcdf_err
subroutine, public netcdf_err(err, string)
Handle netCDF error codes.
Definition: utils.f90:22
utils::error_handler
subroutine, public error_handler(string, rc)
Handle errors.
Definition: utils.f90:48
source_grid
Read grid specs, date information and land/sea mask for the source data that will be interpolated to ...
Definition: source_grid.F90:12
utils
Utilities.
Definition: utils.f90:8