Line data Source code
1 : module ncdio_atm
2 :
3 : !-----------------------------------------------------------------------
4 : !BOP
5 : !
6 : ! !MODULE: ncdio_atm
7 : !
8 : ! !DESCRIPTION:
9 : ! Generic interfaces to write fields to PIO files
10 : !
11 : ! !USES:
12 :
13 : use pio, only: pio_offset_kind, file_desc_t, var_desc_t, pio_double, &
14 : pio_inq_dimid, pio_max_var_dims, io_desc_t, pio_setframe
15 : use shr_kind_mod, only: r8 => shr_kind_r8
16 : use shr_sys_mod, only: shr_sys_flush ! Standardized system subroutines
17 : use shr_scam_mod, only: shr_scam_getCloseLatLon ! Standardized system subroutines
18 : use spmd_utils, only: masterproc
19 : use cam_abortutils, only: endrun
20 : use scamMod, only: scmlat,scmlon,single_column
21 : use cam_logfile, only: iulog
22 : use string_utils, only: to_lower
23 : !
24 : ! !PUBLIC TYPES:
25 : implicit none
26 :
27 : PRIVATE
28 :
29 : save
30 :
31 : logical :: debug = .false.
32 :
33 : !
34 : !EOP
35 : !
36 : interface infld
37 : module procedure infld_real_1d_2d
38 : module procedure infld_real_2d_2d
39 : module procedure infld_real_2d_3d
40 : module procedure infld_real_3d_3d
41 : end interface
42 :
43 :
44 : public :: infld
45 :
46 : integer STATUS
47 : real(r8) surfdat
48 : !-----------------------------------------------------------------------
49 :
50 : contains
51 :
52 : !-----------------------------------------------------------------------
53 : !BOP
54 : !
55 : ! !IROUTINE: infld_real_1d_2d
56 : !
57 : ! !INTERFACE:
58 204288 : subroutine infld_real_1d_2d(varname, ncid, dimname1, &
59 146688 : dim1b, dim1e, dim2b, dim2e, field, readvar, gridname, timelevel, &
60 : fillvalue)
61 : !
62 : ! !DESCRIPTION:
63 : ! Netcdf I/O of initial real field from netCDF file
64 : ! Read a 1-D field (or slice) into a 2-D variable
65 : !
66 : ! !USES
67 : !
68 :
69 : use pio, only: pio_get_var, pio_read_darray, pio_setdebuglevel
70 : use pio, only: PIO_MAX_NAME, pio_inquire, pio_inq_dimname
71 : use cam_grid_support, only: cam_grid_check, cam_grid_get_decomp, cam_grid_id, &
72 : cam_grid_dimensions
73 : use cam_pio_utils, only: cam_pio_check_var, cam_pio_inq_var_fill
74 :
75 : !
76 : ! !ARGUMENTS:
77 : implicit none
78 : character(len=*), intent(in) :: varname ! variable name
79 : type(file_desc_t), intent(inout) :: ncid ! input unit
80 : character(len=*), intent(in) :: dimname1 ! name of 1st array dimensions of field on file (array order)
81 : integer, intent(in) :: dim1b ! start of first dimension of array to be returned
82 : integer, intent(in) :: dim1e ! end of first dimension of array to be returned
83 : integer, intent(in) :: dim2b ! start of second dimension of array to be returned
84 : integer, intent(in) :: dim2e ! end of second dimension of array to be returned
85 : real(r8), target, intent(out) :: field(dim1b:dim1e,dim2b:dim2e) ! array to be returned (decomposed or global)
86 : logical, intent(out) :: readvar ! true => variable is on initial dataset
87 : character(len=*), optional, intent(in) :: gridname ! Name of variable's grid
88 : integer, optional, intent(in) :: timelevel
89 : real(r8), optional, intent(out) :: fillvalue
90 : !
91 : !EOP
92 : !
93 : ! !LOCAL VARIABLES:
94 : type(io_desc_t), pointer :: iodesc
95 : integer :: grid_id ! grid ID for data mapping
96 : integer :: i, j ! indices
97 : integer :: ierr ! error status
98 : type(var_desc_t) :: varid ! variable id
99 : integer :: no_fill
100 : integer :: arraydimsize(2) ! field dimension lengths
101 :
102 : integer :: ndims ! number of dimensions
103 : integer :: dimids(PIO_MAX_VAR_DIMS) ! file variable dims
104 : integer :: dimlens(PIO_MAX_VAR_DIMS) ! file variable shape
105 : integer :: grid_dimlens(2)
106 :
107 : ! Offsets for reading global variables
108 : integer :: strt(1) = 1 ! start ncol index for netcdf 1-d
109 : integer :: cnt (1) = 1 ! ncol count for netcdf 1-d
110 : character(len=PIO_MAX_NAME) :: tmpname
111 : character(len=128) :: errormsg
112 :
113 : logical :: readvar_tmp ! if true, variable is on tape
114 : character(len=*), parameter :: subname='INFLD_REAL_1D_2D' ! subroutine name
115 :
116 : ! For SCAM
117 : real(r8) :: closelat, closelon
118 : integer :: lonidx, latidx
119 :
120 146688 : nullify(iodesc)
121 :
122 : !
123 : !-----------------------------------------------------------------------
124 : !
125 : ! call pio_setdebuglevel(3)
126 :
127 : !
128 : ! Error conditions
129 : !
130 146688 : if (present(gridname)) then
131 146688 : grid_id = cam_grid_id(trim(gridname))
132 : else
133 0 : grid_id = cam_grid_id('physgrid')
134 : end if
135 146688 : if (.not. cam_grid_check(grid_id)) then
136 0 : if(masterproc) then
137 0 : if (present(gridname)) then
138 0 : write(errormsg, *)': invalid gridname, "',trim(gridname),'", specified for field ',trim(varname)
139 : else
140 0 : write(errormsg, *)': Internal error, no "physgrid" gridname'
141 : end if
142 : end if
143 0 : call endrun(trim(subname)//errormsg)
144 : end if
145 :
146 : ! Get the number of columns in the global grid.
147 146688 : call cam_grid_dimensions(grid_id, grid_dimlens)
148 :
149 146688 : if (debug .and. masterproc) then
150 0 : if (present(gridname)) then
151 0 : write(iulog, '(5a)') trim(subname),': field = ',trim(varname),', grid = ',trim(gridname)
152 : else
153 0 : write(iulog, '(4a)') trim(subname),': field = ',trim(varname),', grid = physgrid'
154 : end if
155 0 : call shr_sys_flush(iulog)
156 : end if
157 : !
158 : ! Read netCDF file
159 : !
160 : !
161 : ! Check if field is on file; get netCDF variable id
162 : !
163 146688 : call cam_pio_check_var(ncid, varname, varid, ndims, dimids, dimlens, readvar_tmp)
164 : !
165 : ! If field is on file:
166 : !
167 146688 : if (readvar_tmp) then
168 144384 : if (debug .and. masterproc) then
169 0 : write(iulog, '(2a,5(i0,a))') trim(subname),': field(', &
170 0 : dim1b,':',dim1e,',',dim2b,':',dim2e, '), file(',dimlens(1),')'
171 0 : call shr_sys_flush(iulog)
172 : end if
173 : !
174 : ! Get array dimension id's and sizes
175 : !
176 144384 : arraydimsize(1) = (dim1e - dim1b + 1)
177 144384 : arraydimsize(2) = (dim2e - dim2b + 1)
178 433152 : do j = 1, 2
179 433152 : if (arraydimsize(j) /= size(field, j)) then
180 0 : write(errormsg, *) ': Mismatch between array bounds and field size for ', &
181 0 : trim(varname), ', dimension', j
182 0 : call endrun(trim(subname)//errormsg)
183 : end if
184 : end do
185 :
186 144384 : if (ndims > 2) then
187 0 : call endrun(trim(subname)//': too many dimensions for '//trim(varname))
188 144384 : else if (ndims < 1) then
189 0 : call endrun(trim(subname)//': too few dimensions for '//trim(varname))
190 : else
191 : ! Check that the number of columns in the file matches the number of
192 : ! columns in the grid object.
193 144384 : if (dimlens(1) /= grid_dimlens(1)) then
194 0 : readvar = .false.
195 0 : return
196 : end if
197 :
198 : ! Check to make sure that the second dimension is time
199 144384 : if (ndims == 2) then
200 135936 : ierr = pio_inq_dimname(ncid, dimids(2), tmpname)
201 135936 : if (trim(tmpname) /= 'time') then
202 0 : call endrun(trim(subname)//': dimension mismatch for '//trim(varname))
203 : end if
204 : end if
205 : end if
206 :
207 144384 : if(ndims == 2) then
208 135936 : if(present(timelevel)) then
209 135168 : call pio_setframe(ncid, varid, int(timelevel,kind=pio_offset_kind))
210 : else
211 768 : call pio_setframe(ncid, varid, int(1,kind=pio_offset_kind))
212 : end if
213 135936 : ndims = ndims - 1
214 : end if
215 :
216 : ! NB: strt and cnt were initialized to 1
217 144384 : if (single_column) then
218 : !!XXgoldyXX: Clearly, this will not work for an unstructured dycore
219 0 : call endrun(trim(subname)//': SCAM not supported in this configuration')
220 : else
221 : ! All distributed array processing
222 0 : call cam_grid_get_decomp(grid_id, arraydimsize, dimlens(1:ndims), &
223 144384 : pio_double, iodesc)
224 144384 : call pio_read_darray(ncid, varid, iodesc, field, ierr)
225 144384 : if (present(fillvalue)) then
226 5376 : ierr = cam_pio_inq_var_fill(ncid, varid, fillvalue)
227 : end if
228 : end if
229 :
230 :
231 144384 : if (masterproc) write(iulog,*) subname//': read field '//trim(varname)
232 :
233 : end if ! end of readvar_tmp
234 :
235 146688 : readvar = readvar_tmp
236 :
237 146688 : return
238 :
239 146688 : end subroutine infld_real_1d_2d
240 :
241 : !-----------------------------------------------------------------------
242 : !BOP
243 : !
244 : ! !IROUTINE: infld_real_2d_2d
245 : !
246 : ! !INTERFACE:
247 141312 : subroutine infld_real_2d_2d(varname, ncid, dimname1, dimname2, &
248 141312 : dim1b, dim1e, dim2b, dim2e, field, readvar, gridname, timelevel, &
249 : fillvalue)
250 : !
251 : ! !DESCRIPTION:
252 : ! Netcdf I/O of initial real field from netCDF file
253 : ! Read a 2-D field (or slice) into a 2-D variable
254 : !
255 : ! !USES
256 : !
257 :
258 146688 : use pio, only: pio_get_var, pio_read_darray, pio_setdebuglevel
259 : use pio, only: PIO_MAX_NAME, pio_inquire, pio_inq_dimname
260 : use cam_grid_support, only: cam_grid_check, cam_grid_get_decomp, cam_grid_id
261 : use cam_pio_utils, only: cam_permute_array, calc_permutation
262 : use cam_pio_utils, only: cam_pio_check_var, cam_pio_inq_var_fill
263 :
264 : !
265 : ! !ARGUMENTS:
266 : implicit none
267 : character(len=*), intent(in) :: varname ! variable name
268 : type(file_desc_t), intent(inout) :: ncid ! input unit
269 : character(len=*), intent(in) :: dimname1 ! name of 1st array dimensions of field on file (array order)
270 : character(len=*), intent(in) :: dimname2 ! name of 2nd array dimensions of field on file (array order)
271 : integer, intent(in) :: dim1b ! start of first dimension of array to be returned
272 : integer, intent(in) :: dim1e ! end of first dimension of array to be returned
273 : integer, intent(in) :: dim2b ! start of second dimension of array to be returned
274 : integer, intent(in) :: dim2e ! end of second dimension of array to be returned
275 : real(r8), target, intent(out) :: field(dim1b:dim1e,dim2b:dim2e) ! array to be returned (decomposed or global)
276 : logical, intent(out) :: readvar ! true => variable is on initial dataset
277 : character(len=*), optional, intent(in) :: gridname ! Name of variable's grid
278 : integer, optional, intent(in) :: timelevel
279 : real(r8), optional, intent(out) :: fillvalue
280 : !
281 : !EOP
282 : !
283 : ! !LOCAL VARIABLES:
284 : type(io_desc_t), pointer :: iodesc
285 : integer :: grid_id ! grid ID for data mapping
286 : integer :: i, j ! indices
287 : integer :: ierr ! error status
288 : type(var_desc_t) :: varid ! variable id
289 :
290 : integer :: arraydimsize(2) ! field dimension lengths
291 : integer :: arraydimids(2) ! Dimension IDs
292 : integer :: permutation(2)
293 : logical :: ispermuted
294 :
295 : integer :: ndims ! number of dimensions on file
296 : integer :: dimids(PIO_MAX_VAR_DIMS) ! file variable dims
297 : integer :: dimlens(PIO_MAX_VAR_DIMS) ! file variable shape
298 :
299 : ! Offsets for reading global variables
300 : integer :: strt(2) ! start lon, lat indices for netcdf 2-d
301 : integer :: cnt (2) ! lon, lat counts for netcdf 2-d
302 : character(len=PIO_MAX_NAME) :: tmpname
303 : character(len=128) :: errormsg
304 :
305 141312 : real(r8), pointer :: tmp2d(:,:) ! input data for permutation
306 :
307 : logical :: readvar_tmp ! if true, variable is on tape
308 : character(len=*), parameter :: subname='INFLD_REAL_2D_2D' ! subroutine name
309 : character(len=PIO_MAX_NAME) :: field_dnames(2)
310 :
311 : ! For SCAM
312 : real(r8) :: closelat, closelon
313 : integer :: lonidx, latidx
314 :
315 141312 : nullify(iodesc)
316 :
317 : !
318 : !-----------------------------------------------------------------------
319 : !
320 : ! call pio_setdebuglevel(3)
321 :
322 : ! Should we be using a different interface?
323 141312 : if ((trim(dimname1) == trim(dimname2)) .or. (len_trim(dimname2) == 0)) then
324 : call infld(varname, ncid, dimname1, dim1b, dim1e, dim2b, dim2e, &
325 141312 : field, readvar, gridname, timelevel)
326 : else
327 :
328 : !
329 : ! Error conditions
330 : !
331 0 : if (present(gridname)) then
332 0 : grid_id = cam_grid_id(trim(gridname))
333 : else
334 0 : grid_id = cam_grid_id('physgrid')
335 : end if
336 0 : if (.not. cam_grid_check(grid_id)) then
337 0 : if(masterproc) then
338 0 : if (present(gridname)) then
339 0 : write(errormsg, *)': invalid gridname, "',trim(gridname),'", specified for field ',trim(varname)
340 : else
341 0 : write(errormsg, *)': Internal error, no "physgrid" gridname'
342 : end if
343 : end if
344 0 : call endrun(trim(subname)//errormsg)
345 : end if
346 :
347 0 : if (debug .and. masterproc) then
348 0 : if (present(gridname)) then
349 0 : write(iulog, '(5a)') trim(subname),': field = ',trim(varname),', grid = ',trim(gridname)
350 : else
351 0 : write(iulog, '(4a)') trim(subname),': field = ',trim(varname),', grid = physgrid'
352 : end if
353 0 : call shr_sys_flush(iulog)
354 : end if
355 :
356 : !
357 : ! Read netCDF file
358 : !
359 : !
360 : ! Check if field is on file; get netCDF variable id
361 : !
362 0 : call cam_pio_check_var(ncid, varname, varid, ndims, dimids, dimlens, readvar_tmp)
363 : !
364 : ! If field is on file:
365 : !
366 0 : if (readvar_tmp) then
367 0 : if (debug .and. masterproc) then
368 0 : write(iulog, '(2a,6(i0,a))') trim(subname),': field(', &
369 0 : dim1b,':',dim1e,',',dim2b,':',dim2e, &
370 0 : '), file(',dimlens(1),',',dimlens(2),')'
371 0 : call shr_sys_flush(iulog)
372 : end if
373 : !
374 : ! Get array dimension id's and sizes
375 : !
376 0 : ierr = PIO_inq_dimid(ncid, dimname1, arraydimids(1))
377 0 : ierr = PIO_inq_dimid(ncid, dimname2, arraydimids(2))
378 0 : arraydimsize(1) = (dim1e - dim1b + 1)
379 0 : arraydimsize(2) = (dim2e - dim2b + 1)
380 0 : do j = 1, 2
381 0 : if (arraydimsize(j) /= size(field, j)) then
382 0 : write(errormsg, *) ': Mismatch between array bounds and field size for ', &
383 0 : trim(varname), ', dimension', j
384 0 : call endrun(trim(subname)//errormsg)
385 : end if
386 : end do
387 :
388 0 : if (ndims > 3) then
389 0 : call endrun(trim(subname)//': too many dimensions for '//trim(varname))
390 0 : else if (ndims < 2) then
391 0 : call endrun(trim(subname)//': too few dimensions for '//trim(varname))
392 : else
393 : ! Check to make sure that the third dimension is time
394 0 : if (ndims == 3) then
395 0 : ierr = pio_inq_dimname(ncid, dimids(3), tmpname)
396 0 : if (trim(tmpname) /= 'time') then
397 0 : call endrun(trim(subname)//': dimension mismatch for '//trim(varname))
398 : end if
399 : end if
400 : end if
401 :
402 0 : if(ndims == 3) then
403 0 : if(present(timelevel)) then
404 0 : call pio_setframe(ncid, varid, int(timelevel,kind=pio_offset_kind))
405 : else
406 0 : call pio_setframe(ncid, varid, int(1,kind=pio_offset_kind))
407 : end if
408 : end if
409 :
410 0 : field_dnames(1) = dimname1
411 0 : field_dnames(2) = dimname2
412 0 : if (single_column) then
413 : ! This could be generalized but for now only handles a single point
414 0 : strt(1) = dim1b
415 0 : strt(2) = dim2b
416 0 : cnt = arraydimsize
417 0 : call shr_scam_getCloseLatLon(ncid,scmlat,scmlon,closelat,closelon,latidx,lonidx)
418 0 : if (trim(field_dnames(1)) == 'lon') then
419 0 : strt(1) = lonidx ! First dim always lon for Eulerian dycore
420 : else
421 0 : call endrun(trim(subname)//': lon should be first dimension for '//trim(varname))
422 : end if
423 0 : if (trim(field_dnames(2)) == 'lat') then
424 0 : strt(2) = latidx
425 : else
426 0 : call endrun(trim(subname)//': lat dimension not found for '//trim(varname))
427 : end if
428 :
429 : ! Check for permuted dimensions ('out of order' array)
430 0 : call calc_permutation(dimids, arraydimids, permutation, ispermuted)
431 0 : if (ispermuted) then
432 0 : call cam_permute_array(strt, permutation)
433 0 : call cam_permute_array(cnt, permutation)
434 0 : allocate(tmp2d(1:cnt(1), 1:cnt(2)))
435 0 : ierr = pio_get_var(ncid, varid, strt, cnt, tmp2d)
436 0 : do j = dim2b, dim2e
437 0 : do i = dim1b, dim1e
438 : ! We don't need strt anymore, reuse it
439 0 : strt(1) = i - dim1b + 1
440 0 : strt(2) = j - dim2b + 1
441 0 : call cam_permute_array(strt, permutation)
442 0 : field(i,j) = tmp2d(strt(1), strt(2))
443 : end do
444 : end do
445 0 : deallocate(tmp2d)
446 : else
447 0 : ierr = pio_get_var(ncid, varid, strt, cnt, field)
448 : end if
449 : else
450 : ! All distributed array processing
451 : call cam_grid_get_decomp(grid_id, arraydimsize, dimlens(1:2), &
452 0 : pio_double, iodesc, field_dnames=field_dnames)
453 0 : call pio_read_darray(ncid, varid, iodesc, field, ierr)
454 0 : if (present(fillvalue)) then
455 0 : ierr = cam_pio_inq_var_fill(ncid, varid, fillvalue)
456 : end if
457 : end if
458 :
459 0 : if (masterproc) write(iulog,*) subname//': read field '//trim(varname)
460 :
461 : end if ! end of readvar_tmp
462 :
463 0 : readvar = readvar_tmp
464 :
465 : end if ! end of call infld_real_1d_2d instead
466 :
467 282624 : end subroutine infld_real_2d_2d
468 :
469 :
470 : !-----------------------------------------------------------------------
471 : !BOP
472 : !
473 : ! !IROUTINE: infld_real_2d_3d
474 : !
475 : ! !INTERFACE:
476 125184 : subroutine infld_real_2d_3d(varname, ncid, dimname1, dimname2, &
477 : dim1b, dim1e, dim2b, dim2e, dim3b, dim3e, &
478 125184 : field, readvar, gridname, timelevel, fillvalue)
479 : !
480 : ! !DESCRIPTION:
481 : ! Netcdf I/O of initial real field from netCDF file
482 : ! Read a 2-D field (or slice) into a 3-D variable
483 : !
484 : ! !USES
485 : !
486 :
487 141312 : use pio, only: pio_get_var, pio_read_darray, pio_setdebuglevel
488 : use pio, only: PIO_MAX_NAME, pio_inquire, pio_inq_dimname
489 : use cam_grid_support, only: cam_grid_check, cam_grid_get_decomp, cam_grid_id, &
490 : cam_grid_dimensions
491 : use cam_pio_utils, only: cam_permute_array, calc_permutation
492 : use cam_pio_utils, only: cam_pio_check_var, cam_pio_inq_var_fill
493 :
494 : !
495 : ! !ARGUMENTS:
496 : implicit none
497 : character(len=*), intent(in) :: varname ! variable name
498 : type(file_desc_t), intent(inout) :: ncid ! input unit
499 : character(len=*), intent(in) :: dimname1 ! name of 1st array dimensions of field on file (array order)
500 : character(len=*), intent(in) :: dimname2 ! name of 2nd array dimensions of field on file (array order)
501 : integer, intent(in) :: dim1b ! start of first dimension of array to be returned
502 : integer, intent(in) :: dim1e ! end of first dimension of array to be returned
503 : integer, intent(in) :: dim2b ! start of second dimension of array to be returned
504 : integer, intent(in) :: dim2e ! end of second dimension of array to be returned
505 : integer, intent(in) :: dim3b ! start of third dimension of array to be returned
506 : integer, intent(in) :: dim3e ! end of third dimension of array to be returned
507 : real(r8), target, intent(out) :: field(dim1b:dim1e,dim2b:dim2e,dim3b:dim3e) ! array to be returned (decomposed or global)
508 : logical, intent(out) :: readvar ! true => variable is on initial dataset
509 : character(len=*), optional, intent(in) :: gridname ! Name of variable's grid
510 : integer, optional, intent(in) :: timelevel
511 : real(r8), optional, intent(out) :: fillvalue
512 : !
513 : !EOP
514 : !
515 : ! !LOCAL VARIABLES:
516 : type(io_desc_t), pointer :: iodesc
517 : integer :: grid_id ! grid ID for data mapping
518 : integer :: i, j, k ! indices
519 : integer :: ierr ! error status
520 : type(var_desc_t) :: varid ! variable id
521 :
522 : integer :: arraydimsize(3) ! field dimension lengths
523 : integer :: arraydimids(2) ! Dimension IDs
524 : integer :: permutation(2)
525 : logical :: ispermuted
526 :
527 : integer :: ndims ! number of dimensions
528 : integer :: dimids(PIO_MAX_VAR_DIMS) ! file variable dims
529 : integer :: dimlens(PIO_MAX_VAR_DIMS) ! file variable shape
530 : integer :: grid_dimlens(2)
531 :
532 : ! Offsets for reading global variables
533 : integer :: strt(3) = 1 ! start ncol, lev indices for netcdf 2-d
534 : integer :: cnt (3) = 1 ! ncol, lev counts for netcdf 2-d
535 : character(len=PIO_MAX_NAME) :: tmpname
536 :
537 : real(r8), pointer :: tmp3d(:,:,:) ! input data for permutation
538 :
539 : logical :: readvar_tmp ! if true, variable is on tape
540 : character(len=*), parameter :: subname='INFLD_REAL_2D_3D' ! subroutine name
541 : character(len=128) :: errormsg
542 : character(len=PIO_MAX_NAME) :: field_dnames(2)
543 : character(len=PIO_MAX_NAME) :: file_dnames(3)
544 :
545 : ! For SCAM
546 : real(r8) :: closelat, closelon
547 : integer :: lonidx, latidx
548 :
549 125184 : nullify(iodesc)
550 :
551 : !
552 : !-----------------------------------------------------------------------
553 : !
554 : ! call pio_setdebuglevel(3)
555 :
556 : !
557 : ! Error conditions
558 : !
559 125184 : if (present(gridname)) then
560 125184 : grid_id = cam_grid_id(trim(gridname))
561 : else
562 0 : grid_id = cam_grid_id('physgrid')
563 : end if
564 125184 : if (.not. cam_grid_check(grid_id)) then
565 0 : if(masterproc) then
566 0 : if (present(gridname)) then
567 0 : write(errormsg, *)': invalid gridname, "',trim(gridname),'", specified for field ',trim(varname)
568 : else
569 0 : write(errormsg, *)': Internal error, no "physgrid" gridname'
570 : end if
571 : end if
572 0 : call endrun(trim(subname)//errormsg)
573 : end if
574 :
575 : ! Get the number of columns in the global grid.
576 125184 : call cam_grid_dimensions(grid_id, grid_dimlens)
577 :
578 125184 : if (debug .and. masterproc) then
579 0 : if (present(gridname)) then
580 0 : write(iulog, '(5a)') trim(subname),': field = ',trim(varname),', grid = ',trim(gridname)
581 : else
582 0 : write(iulog, '(4a)') trim(subname),': field = ',trim(varname),', grid = physgrid'
583 : end if
584 0 : call shr_sys_flush(iulog)
585 : end if
586 :
587 : !
588 : ! Read netCDF file
589 : !
590 : !
591 : ! Check if field is on file; get netCDF variable id
592 : !
593 : call cam_pio_check_var(ncid, varname, varid, ndims, dimids, dimlens, &
594 125184 : readvar_tmp, dimnames=file_dnames)
595 :
596 : ! If field is on file:
597 : !
598 125184 : if (readvar_tmp) then
599 121344 : if (debug .and. masterproc) then
600 0 : write(iulog, '(2a,8(i0,a))') trim(subname),': field(', &
601 0 : dim1b,':',dim1e,',',dim2b,':',dim2e,',',dim3b,':',dim3e, &
602 0 : '), file(',dimlens(1),',',dimlens(2),')'
603 0 : call shr_sys_flush(iulog)
604 : end if
605 : !
606 : ! Get array dimension id's and sizes
607 : !
608 121344 : arraydimsize(1) = (dim1e - dim1b + 1)
609 121344 : arraydimsize(2) = (dim2e - dim2b + 1)
610 121344 : arraydimsize(3) = (dim3e - dim3b + 1)
611 485376 : do j = 1, 3
612 485376 : if (arraydimsize(j) /= size(field, j)) then
613 0 : write(errormsg, *) ': Mismatch between array bounds and field size for ', &
614 0 : trim(varname), ', dimension', j
615 0 : call endrun(trim(subname)//errormsg)
616 : end if
617 : end do
618 :
619 121344 : if (ndims > 3) then
620 0 : call endrun(trim(subname)//': too many dimensions for '//trim(varname))
621 121344 : else if (ndims < 2) then
622 0 : call endrun(trim(subname)//': too few dimensions for '//trim(varname))
623 : else
624 : ! Check that the number of columns in the file matches the number of
625 : ! columns in the grid object.
626 121344 : if (dimlens(1) /= grid_dimlens(1) .and. dimlens(2) /= grid_dimlens(1)) then
627 0 : readvar = .false.
628 0 : return
629 : end if
630 :
631 : ! Check to make sure that the 3rd dimension is time
632 121344 : if (ndims == 3) then
633 112128 : ierr = pio_inq_dimname(ncid, dimids(3), tmpname)
634 224256 : if (to_lower(trim(tmpname)) /= 'time') then
635 112128 : call endrun(trim(subname)//': dimension mismatch for '//trim(varname))
636 : end if
637 : end if
638 : end if
639 :
640 121344 : if(ndims == 3) then
641 112128 : if(present(timelevel)) then
642 78336 : call pio_setframe(ncid, varid, int(timelevel,kind=pio_offset_kind))
643 : else
644 33792 : call pio_setframe(ncid, varid, int(1,kind=pio_offset_kind))
645 : end if
646 112128 : ndims = ndims - 1
647 : end if
648 :
649 121344 : field_dnames(1) = dimname1
650 121344 : field_dnames(2) = dimname2
651 : ! NB: strt and cnt were initialized to 1
652 121344 : if (single_column) then
653 : !!XXgoldyXX: Clearly, this will not work for an unstructured dycore
654 : ! Check for permuted dimensions ('out of order' array)
655 : ! call calc_permutation(dimids(1:2), arraydimids, permutation, ispermuted)
656 0 : call endrun(trim(subname)//': SCAM not supported in this configuration')
657 : else
658 : ! All distributed array processing
659 : call cam_grid_get_decomp(grid_id, arraydimsize, dimlens(1:2), &
660 : pio_double, iodesc, field_dnames=field_dnames, &
661 121344 : file_dnames=file_dnames(1:2))
662 121344 : call pio_read_darray(ncid, varid, iodesc, field, ierr)
663 121344 : if (present(fillvalue)) then
664 33792 : ierr = cam_pio_inq_var_fill(ncid, varid, fillvalue)
665 : end if
666 : end if
667 :
668 121344 : if (masterproc) write(iulog,*) subname//': read field '//trim(varname)
669 :
670 : end if ! end of readvar_tmp
671 :
672 125184 : readvar = readvar_tmp
673 :
674 125184 : return
675 :
676 250368 : end subroutine infld_real_2d_3d
677 :
678 : !-----------------------------------------------------------------------
679 : !BOP
680 : !
681 : ! !IROUTINE: infld_real_3d_3d
682 : !
683 : ! !INTERFACE:
684 57600 : subroutine infld_real_3d_3d(varname, ncid, dimname1, dimname2, dimname3, &
685 : dim1b, dim1e, dim2b, dim2e, dim3b, dim3e, &
686 57600 : field, readvar, gridname, timelevel, fillvalue)
687 : !
688 : ! !DESCRIPTION:
689 : ! Netcdf I/O of initial real field from netCDF file
690 : ! Read a 3-D field (or slice) into a 3-D variable
691 : !
692 : ! !USES
693 : !
694 :
695 125184 : use pio, only: pio_get_var, pio_read_darray, pio_setdebuglevel
696 : use pio, only: PIO_MAX_NAME, pio_inquire, pio_inq_dimname
697 : use cam_grid_support, only: cam_grid_check, cam_grid_get_decomp, cam_grid_id
698 : use cam_pio_utils, only: cam_permute_array, calc_permutation
699 : use cam_pio_utils, only: cam_pio_check_var, cam_pio_inq_var_fill
700 :
701 : !
702 : ! !ARGUMENTS:
703 : implicit none
704 : character(len=*), intent(in) :: varname ! variable name
705 : type(file_desc_t), intent(inout) :: ncid ! input unit
706 : character(len=*), intent(in) :: dimname1 ! name of 1st array dimensions of field on file (array order)
707 : character(len=*), intent(in) :: dimname2 ! name of 2nd array dimensions of field on file (array order)
708 : character(len=*), intent(in) :: dimname3 ! name of 3rd array dimensions of field on file (array order)
709 : integer, intent(in) :: dim1b ! start of first dimension of array to be returned
710 : integer, intent(in) :: dim1e ! end of first dimension of array to be returned
711 : integer, intent(in) :: dim2b ! start of second dimension of array to be returned
712 : integer, intent(in) :: dim2e ! end of second dimension of array to be returned
713 : integer, intent(in) :: dim3b ! start of third dimension of array to be returned
714 : integer, intent(in) :: dim3e ! end of third dimension of array to be returned
715 : real(r8), target, intent(out) :: field(dim1b:dim1e,dim2b:dim2e,dim3b:dim3e) ! array to be returned (decomposed or global)
716 : logical, intent(out) :: readvar ! true => variable is on initial dataset
717 : character(len=*), optional, intent(in) :: gridname ! Name of variable's grid
718 : integer, optional, intent(in) :: timelevel
719 : real(r8), optional, intent(out) :: fillvalue
720 : !
721 : !EOP
722 : !
723 : ! !LOCAL VARIABLES:
724 : type(io_desc_t), pointer :: iodesc
725 : integer :: grid_id ! grid ID for data mapping
726 : integer :: i, j, k ! indices
727 : integer :: ierr ! error status
728 : type(var_desc_t) :: varid ! variable id
729 :
730 : integer :: arraydimsize(3) ! field dimension lengths
731 : integer :: arraydimids(3) ! Dimension IDs
732 : integer :: permutation(3)
733 : logical :: ispermuted
734 :
735 : integer :: ndims ! number of dimensions
736 : integer :: pdims ! number of dimensions w/o timeslice
737 : integer :: dimids(PIO_MAX_VAR_DIMS) ! file variable dims
738 : integer :: dimlens(PIO_MAX_VAR_DIMS) ! file variable shape
739 :
740 : ! Offsets for reading global variables
741 : integer :: strt(3) ! start lon, lev, lat indices for netcdf 3-d
742 : integer :: cnt (3) ! lon, lat counts for netcdf 3-d
743 : character(len=PIO_MAX_NAME) :: tmpname
744 :
745 57600 : real(r8), pointer :: tmp3d(:,:,:) ! input data for permutation
746 :
747 : logical :: readvar_tmp ! if true, variable is on tape
748 : character(len=*), parameter :: subname='INFLD_REAL_3D_3D' ! subroutine name
749 : character(len=128) :: errormsg
750 : character(len=PIO_MAX_NAME) :: field_dnames(3)
751 : character(len=PIO_MAX_NAME) :: file_dnames(4)
752 :
753 : ! For SCAM
754 : real(r8) :: closelat, closelon
755 : integer :: lonidx, latidx
756 :
757 57600 : nullify(iodesc)
758 :
759 : !
760 : !-----------------------------------------------------------------------
761 : !
762 : ! call pio_setdebuglevel(3)
763 :
764 : ! Should we be using a different interface?
765 57600 : if ((trim(dimname1) == trim(dimname2)) .or. (len_trim(dimname2) == 0)) then
766 : call infld(varname, ncid, dimname1, dimname3, &
767 : dim1b, dim1e, dim2b, dim2e, dim3b, dim3e, &
768 0 : field, readvar, gridname, timelevel)
769 57600 : else if ((trim(dimname1) == trim(dimname3)) .or. (len_trim(dimname3) == 0)) then
770 : call infld(varname, ncid, dimname1, dimname2, &
771 : dim1b, dim1e, dim2b, dim2e, dim3b, dim3e, &
772 57600 : field, readvar, gridname, timelevel)
773 : else
774 :
775 : !
776 : ! Error conditions
777 : !
778 0 : if (present(gridname)) then
779 0 : grid_id = cam_grid_id(trim(gridname))
780 : else
781 0 : grid_id = cam_grid_id('physgrid')
782 : end if
783 0 : if (.not. cam_grid_check(grid_id)) then
784 0 : if(masterproc) then
785 0 : if (present(gridname)) then
786 0 : write(errormsg, *)': invalid gridname, "',trim(gridname),'", specified for field ',trim(varname)
787 : else
788 0 : write(errormsg, *)': Internal error, no "physgrid" gridname'
789 : end if
790 : end if
791 0 : call endrun(trim(subname)//errormsg)
792 : end if
793 :
794 0 : if (debug .and. masterproc) then
795 0 : if (present(gridname)) then
796 0 : write(iulog, '(5a)') trim(subname),': field = ',trim(varname),', grid = ',trim(gridname)
797 : else
798 0 : write(iulog, '(4a)') trim(subname),': field = ',trim(varname),', grid = physgrid'
799 : end if
800 0 : call shr_sys_flush(iulog)
801 : end if
802 :
803 : !
804 : ! Read netCDF file
805 : !
806 : !
807 : ! Check if field is on file; get netCDF variable id
808 : !
809 : call cam_pio_check_var(ncid, varname, varid, ndims, dimids, dimlens, &
810 0 : readvar_tmp, dimnames=file_dnames)
811 : !
812 : ! If field is on file:
813 : !
814 0 : if (readvar_tmp) then
815 0 : if (debug .and. masterproc) then
816 0 : write(iulog, '(2a,9(i0,a))') trim(subname),': field(', &
817 0 : dim1b,':',dim1e,',',dim2b,':',dim2e,',',dim3b,':',dim3e, &
818 0 : '), file(',dimlens(1),',',dimlens(2),',',dimlens(3),')'
819 0 : call shr_sys_flush(iulog)
820 : end if
821 : !
822 : ! Get array dimension id's and sizes
823 : !
824 0 : ierr = PIO_inq_dimid(ncid, dimname1, arraydimids(1))
825 0 : ierr = PIO_inq_dimid(ncid, dimname2, arraydimids(2))
826 0 : ierr = PIO_inq_dimid(ncid, dimname3, arraydimids(3))
827 0 : arraydimsize(1) = (dim1e - dim1b + 1)
828 0 : arraydimsize(2) = (dim2e - dim2b + 1)
829 0 : arraydimsize(3) = (dim3e - dim3b + 1)
830 :
831 0 : do j = 1, 3
832 0 : if (arraydimsize(j) /= size(field, j)) then
833 0 : write(errormsg, *) ': Mismatch between array bounds and field size for ', &
834 0 : trim(varname), ', dimension', j
835 0 : call endrun(trim(subname)//errormsg)
836 : end if
837 : end do
838 :
839 0 : pdims = ndims
840 0 : if (ndims > 4) then
841 0 : call endrun(trim(subname)//': too many dimensions for '//trim(varname))
842 0 : else if (ndims < 3) then
843 0 : call endrun(trim(subname)//': too few dimensions for '//trim(varname))
844 : else
845 : ! Check to make sure that the fourth dimension is time
846 0 : if (ndims == 4) then
847 0 : ierr = pio_inq_dimname(ncid, dimids(4), tmpname)
848 0 : if (trim(tmpname) /= 'time') then
849 0 : call endrun(trim(subname)//': dimension mismatch for '//trim(varname))
850 : end if
851 0 : pdims = 3
852 : end if
853 : end if
854 :
855 0 : if(ndims == 4) then
856 0 : if(present(timelevel)) then
857 0 : call pio_setframe(ncid, varid, int(timelevel,kind=pio_offset_kind))
858 : else
859 0 : call pio_setframe(ncid, varid, int(1,kind=pio_offset_kind))
860 : end if
861 : end if
862 :
863 0 : field_dnames(1) = dimname1
864 0 : field_dnames(2) = dimname2
865 0 : field_dnames(3) = dimname3
866 :
867 0 : if (single_column) then
868 : ! This could be generalized but for now only handles a single point
869 0 : strt(1) = dim1b
870 0 : strt(2) = dim2b
871 0 : strt(3) = dim3b
872 0 : cnt = arraydimsize
873 0 : call shr_scam_getCloseLatLon(ncid,scmlat,scmlon,closelat,closelon,latidx,lonidx)
874 0 : if (trim(field_dnames(1)) == 'lon') then
875 0 : strt(1) = lonidx ! First dim always lon for Eulerian dycore
876 : else
877 0 : call endrun(trim(subname)//': lon should be first dimension for '//trim(varname))
878 : end if
879 0 : if (trim(field_dnames(2)) == 'lat') then
880 0 : strt(2) = latidx
881 0 : else if (trim(field_dnames(3)) == 'lat') then
882 0 : strt(3) = latidx
883 : else
884 0 : call endrun(trim(subname)//': lat dimension not found for '//trim(varname))
885 : end if
886 :
887 : ! Check for permuted dimensions ('out of order' array)
888 0 : call calc_permutation(dimids, arraydimids, permutation, ispermuted)
889 0 : if (ispermuted) then
890 0 : call cam_permute_array(strt, permutation)
891 0 : call cam_permute_array(cnt, permutation)
892 0 : allocate(tmp3d(1:cnt(1), 1:cnt(2), 1:cnt(3)))
893 0 : ierr = pio_get_var(ncid, varid, strt, cnt, tmp3d)
894 0 : do k = dim3b, dim3e
895 0 : do j = dim2b, dim2e
896 0 : do i = dim1b, dim1e
897 : ! We don't need strt anymore, reuse it
898 0 : strt(1) = i - dim1b + 1
899 0 : strt(2) = j - dim2b + 1
900 0 : strt(3) = k - dim3b + 1
901 0 : call cam_permute_array(strt, permutation)
902 0 : field(i,j,k) = tmp3d(strt(1), strt(2), strt(3))
903 : end do
904 : end do
905 : end do
906 0 : deallocate(tmp3d)
907 : else
908 0 : ierr = pio_get_var(ncid, varid, strt, cnt, field)
909 : end if
910 : else
911 : ! All distributed array processing
912 0 : call cam_grid_get_decomp(grid_id, arraydimsize, dimlens(1:pdims), &
913 : pio_double, iodesc, field_dnames=field_dnames, &
914 0 : file_dnames=file_dnames(1:3))
915 0 : call pio_read_darray(ncid, varid, iodesc, field, ierr)
916 0 : if (present(fillvalue)) then
917 0 : ierr = cam_pio_inq_var_fill(ncid, varid, fillvalue)
918 : end if
919 : end if ! end of single column
920 :
921 0 : if (masterproc) write(iulog,*) subname//': read field '//trim(varname)
922 :
923 : end if ! end of readvar_tmp
924 :
925 0 : readvar = readvar_tmp
926 :
927 : end if ! end of call infld_real_2d_3d instead
928 :
929 115200 : end subroutine infld_real_3d_3d
930 :
931 :
932 : end module ncdio_atm
|
