sfc_climo_gen 1.14.0
Loading...
Searching...
No Matches
source_grid.F90
Go to the documentation of this file.
1
5
12module 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
Read grid specs, date information and land/sea mask for the source data that will be interpolated to ...
Definition source_grid.F90:12
source_grid::source
character(len=75), public source
Original source of the data.
Definition source_grid.F90:22
source_grid::source_grid_cleanup
subroutine, public source_grid_cleanup
Free up memory associated with this module.
Definition source_grid.F90:404
source_grid::num_time_recs
integer, public num_time_recs
Number of time records.
Definition source_grid.F90:27
source_grid::num_records
integer, public num_records
Number of fields times time records.
Definition source_grid.F90:26
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:52
source_grid::field_names
character(len=50), dimension(:), allocatable, public field_names
Names of fields to be processed.
Definition source_grid.F90:21
source_grid::grid_src
type(esmf_grid), public grid_src
ESMF grid object for the source grid.
Definition source_grid.F90:34
source_grid::j_src
integer, public j_src
j dimension of the source grid.
Definition source_grid.F90:25
source_grid::num_fields
integer, public num_fields
Number of fields in the file.
Definition source_grid.F90:28
source_grid::i_src
integer, public i_src
i dimension of the source grid.
Definition source_grid.F90:24
source_grid::day_of_rec
integer, dimension(:), allocatable, public day_of_rec
Day of each time record with respect to Jan 1.
Definition source_grid.F90:31
utils
Utilities.
Definition utils.f90:8
utils::netcdf_err
subroutine, public netcdf_err(err, string)
Handle netCDF error codes.
Definition utils.f90:23
utils::error_handler
subroutine, public error_handler(string, rc)
Handle errors.
Definition utils.f90:49