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