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