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