ROMS
Loading...
Searching...
No Matches
nf_fread3d.F
Go to the documentation of this file.
1#include "cppdefs.h"
3!
4!git $Id$
5!================================================== Hernan G. Arango ===
6! Copyright (c) 2002-2025 The ROMS Group !
7! Licensed under a MIT/X style license !
8! See License_ROMS.md !
9!=======================================================================
10! !
11! This function reads in a generic floating point 3D array from an !
12! input NetCDF file. !
13! !
14! On Input: !
15! !
16! ng Nested grid number (integer) !
17! model Calling model identifier (integer) !
18! ncname NetCDF file name (string) !
19! ncid NetCDF file ID (integer) !
20#if defined PIO_LIB && defined DISTRIBUTE
21!or pioFile PIO file descriptor structure, TYPE(file_desc_t) !
22! pioFile%fh file handler !
23! pioFile%iosystem IO system descriptor (struct) !
24#endif
25! ncvname NetCDF variable name (string) !
26! ncvarid NetCDF variable ID (integer) !
27#if defined PIO_LIB && defined DISTRIBUTE
28!or pioVar PIO variable descriptor structure, TYPE(My_VarDesc) !
29! pioVar%vd variable descriptor TYPE(Var_Desc_t)!
30! pioVar%dkind variable data kind !
31! pioVar%gtype variable C-gridtype !
32#endif
33! tindex NetCDF time record index to read (integer) !
34! gtype C-grid type (integer) !
35#if defined PIO_LIB && defined DISTRIBUTE
36!or pioDesc IO data decomposition descriptor, TYPE(IO_desc_t) !
37#endif
38! Vsize Variable dimensions in NetCDF file (integer 1D array) !
39! LBi I-dimension Lower bound (integer) !
40! UBi I-dimension Upper bound (integer) !
41! LBj J-dimension Lower bound (integer) !
42! UBj J-dimension Upper bound (integer) !
43! LBk K-dimension Lower bound (integer) !
44! UBk K-dimension Upper bound (integer) !
45! Ascl Factor to scale field after reading (real). !
46! Amask Land/Sea mask, if any (real 3D array) !
47! !
48! On Output: !
49! !
50! Amin Field minimum value (real) !
51! Amax Field maximum value (real) !
52! Adat Field to read in (real 3D array) !
53! checksum Field checksum value (32-bit integer; OPTIONAL) !
54! status Result Error flag (integer) !
55! !
56!=======================================================================
57!
58 USE mod_param
59 USE mod_parallel
60 USE mod_iounits
61 USE mod_ncparam
62 USE mod_scalars
63!
64 USE get_hash_mod, ONLY : get_hash
65 USE strings_mod, ONLY : founderror
66!
67 implicit none
68!
69 INTERFACE nf_fread3d
70 MODULE PROCEDURE nf90_fread3d
71#if defined PIO_LIB && defined DISTRIBUTE
72 MODULE PROCEDURE pio_fread3d
73#endif
74 END INTERFACE nf_fread3d
75!
76 CONTAINS
77
78#if defined PARALLEL_IO && defined DISTRIBUTE
79!
80!***********************************************************************
81 FUNCTION nf90_fread3d (ng, model, ncname, ncid, &
82 & ncvname, ncvarid, &
83 & tindex, gtype, Vsize, &
84 & LBi, UBi, LBj, UBj, LBk, UBk, &
85 & Ascl, Amin, Amax, &
86# ifdef MASKING
87 & Amask, &
88# endif
89 & Adat, &
90 & checksum) RESULT (status)
91!***********************************************************************
92!
93 USE mod_netcdf
94!
96# if defined MASKING && defined READ_WATER
97 USE distribute_mod, ONLY : mp_collect
98# endif
100!
101! Imported variable declarations.
102!
103 integer, intent(in) :: ng, model, ncid, ncvarid, tindex, gtype
104 integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk
105 integer, intent(in) :: vsize(4)
106!
107 integer(i8b), intent(out), optional :: checksum
108!
109 real(dp), intent(in) :: ascl
110 real(r8), intent(out) :: amin
111 real(r8), intent(out) :: amax
112!
113 character (len=*), intent(in) :: ncname
114 character (len=*), intent(in) :: ncvname
115!
116# ifdef ASSUMED_SHAPE
117# ifdef MASKING
118 real(r8), intent(in) :: amask(lbi:,lbj:)
119# endif
120 real(r8), intent(out) :: adat(lbi:,lbj:,lbk:)
121# else
122# ifdef MASKING
123 real(r8), intent(in) :: amask(lbi:ubi,lbj:ubj)
124# endif
125 real(r8), intent(out) :: adat(lbi:ubi,lbj:ubj,lbk:ubk)
126# endif
127!
128! Local variable declarations.
129!
130 logical :: lchecksum
131 logical, dimension(3) :: foundit
132!
133 integer :: i, ic, ij, j, jc, k, kc, np, mynpts, npts
134 integer :: imin, imax, isize, jmin, jmax, jsize, ijsize
135 integer :: istr, iend
136 integer :: ioff, joff, koff
137 integer :: ilen, jlen, klen, ijlen
138 integer :: cgrid, mytype, ghost, status, wtype
139
140 integer, dimension(4) :: start, total
141!
142 real(r8) :: afactor, aoffset, aspval
143
144 real(r8), parameter :: inival= 0.0_r8
145
146 real(r8), dimension(2) :: rbuffer
147 real(r8), dimension(3) :: attvalue
148
149# if defined MASKING && defined READ_WATER
150 real(r8), allocatable :: a2d(:)
151# endif
152 real(r8), allocatable :: wrk(:)
153!
154 character (len= 3), dimension(2) :: op_handle
155 character (len=12), dimension(3) :: attname
156
157 character (len=*), parameter :: myfile = &
158 & __FILE__//", nf90_fread3d"
159!
160!-----------------------------------------------------------------------
161! Set starting and ending indices to process.
162!-----------------------------------------------------------------------
163!
164 status=nf90_noerr
165!
166! Set first and last grid point according to staggered C-grid
167! classification. Set the offsets for variables with starting
168! zero-index. Recall the NetCDF does not support a zero-index.
169!
170! Notice that (Imin,Jmin) and (Imax,Jmax) are the corner of the
171! computational tile. If ghost=0, ghost points are not processed.
172! They will be processed elsewhere by the appropriate call to any
173! of the routines in "mp_exchange.F". If ghost=1, the ghost points
174! are read.
175!
176# ifdef NO_READ_GHOST
177 ghost=0 ! non-overlapping, no ghost points
178# else
179 IF (model.eq.iadm) THEN
180 ghost=0 ! non-overlapping, no ghost points
181 ELSE
182 ghost=1 ! overlapping, read ghost points
183 END IF
184# endif
185
186 mytype=gtype
187
188 SELECT CASE (abs(mytype))
189 CASE (p2dvar, p3dvar)
190 cgrid=1
191 isize=iobounds(ng)%xi_psi
192 jsize=iobounds(ng)%eta_psi
193 CASE (r2dvar, r3dvar, w3dvar)
194 cgrid=2
195 isize=iobounds(ng)%xi_rho
196 jsize=iobounds(ng)%eta_rho
197 CASE (u2dvar, u3dvar)
198 cgrid=3
199 isize=iobounds(ng)%xi_u
200 jsize=iobounds(ng)%eta_u
201 CASE (v2dvar, v3dvar)
202 cgrid=4
203 isize=iobounds(ng)%xi_v
204 jsize=iobounds(ng)%eta_v
205 CASE DEFAULT
206 cgrid=2
207 isize=iobounds(ng)%xi_rho
208 jsize=iobounds(ng)%eta_rho
209 END SELECT
210
211 imin=bounds(ng)%Imin(cgrid,ghost,myrank)
212 imax=bounds(ng)%Imax(cgrid,ghost,myrank)
213 jmin=bounds(ng)%Jmin(cgrid,ghost,myrank)
214 jmax=bounds(ng)%Jmax(cgrid,ghost,myrank)
215
216 ilen=imax-imin+1
217 jlen=jmax-jmin+1
218 klen=ubk-lbk+1
219!
220! Check if the following attributes: "scale_factor", "add_offset", and
221! "_FillValue" are present in the input NetCDF variable:
222!
223! If the "scale_value" attribute is present, the data is multiplied by
224! this factor after reading.
225! If the "add_offset" attribute is present, this value is added to the
226! data after reading.
227! If both "scale_factor" and "add_offset" attributes are present, the
228! data are first scaled before the offset is added.
229! If the "_FillValue" attribute is present, the data having this value
230! is treated as missing and it is replaced with zero. This feature it
231! is usually related with the land/sea masking.
232!
233 attname(1)='scale_factor'
234 attname(2)='add_offset '
235 attname(3)='_FillValue '
236
237 CALL netcdf_get_fatt (ng, model, ncname, ncvarid, attname, &
238 & attvalue, foundit, &
239 & ncid = ncid)
240 IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
241 status=ioerror
242 RETURN
243 END IF
244
245 IF (.not.foundit(1)) THEN
246 afactor=1.0_r8
247 ELSE
248 afactor=attvalue(1)
249 END IF
250
251 IF (.not.foundit(2)) THEN
252 aoffset=0.0_r8
253 ELSE
254 aoffset=attvalue(2)
255 END IF
256
257 IF (.not.foundit(3)) THEN
258 aspval=spval_check
259 ELSE
260 aspval=attvalue(3)
261 END IF
262!
263! Initialize checsum value.
264!
265 IF (PRESENT(checksum)) THEN
266 lchecksum=.true.
267 checksum=0_i8b
268 ELSE
269 lchecksum=.false.
270 END IF
271!
272!-----------------------------------------------------------------------
273! Parallel I/O: Read in tile data from requested field and scale it.
274! Processing both water and land points.
275!-----------------------------------------------------------------------
276!
277 IF (gtype.gt.0) THEN
278!
279! Set offsets due the NetCDF dimensions. Recall that some output
280! variables not always start at one.
281!
282 SELECT CASE (abs(mytype))
283 CASE (p2dvar, p3dvar)
284 ioff=0
285 joff=0
286 CASE (r2dvar, r3dvar, w3dvar)
287 ioff=1
288 joff=1
289 CASE (u2dvar, u3dvar)
290 ioff=0
291 joff=1
292 CASE (v2dvar, v3dvar)
293 ioff=1
294 joff=0
295 CASE DEFAULT
296 ioff=1
297 joff=1
298 END SELECT
299
300 IF (lbk.eq.0) THEN
301 koff=1
302 ELSE
303 koff=0
304 END IF
305
306 npts=ilen*jlen*klen
307!
308! Allocate scratch work array.
309!
310 IF (.not.allocated(wrk)) THEN
311 allocate ( wrk(npts) )
312 wrk=0.0_r8
313 END IF
314!
315! Read in data: all parallel nodes read their own tile data.
316!
317 start(1)=imin+ioff
318 total(1)=ilen
319 start(2)=jmin+joff
320 total(2)=jlen
321 start(3)=lbk+koff
322 total(3)=klen
323 start(4)=tindex
324 total(4)=1
325
326 status=nf90_get_var(ncid, ncvarid, wrk, start, total)
327!
328! Scale read data and process fill values, if any. Compute minimum
329! and maximum values.
330!
331 IF (status.eq.nf90_noerr) THEN
332 amin=spval
333 amax=-spval
334 DO i=1,npts
335 IF (abs(wrk(i)).ge.abs(aspval)) THEN
336 wrk(i)=0.0_r8 ! masked with _FillValue
337 ELSE
338 wrk(i)=ascl*(afactor*wrk(i)+aoffset)
339 amin=min(amin,wrk(i))
340 amax=max(amax,wrk(i))
341 END IF
342 END DO
343 IF ((abs(amin).ge.abs(aspval)).and. &
344 & (abs(amax).ge.abs(aspval))) THEN
345 amin=0.0_r8 ! the entire data is all
346 amax=0.0_r8 ! field value, _FillValue
347 END IF
348!
349! Set minimum and maximum values: global reduction.
350!
351 rbuffer(1)=amin
352 op_handle(1)='MIN'
353 rbuffer(2)=amax
354 op_handle(2)='MAX'
355 CALL mp_reduce (ng, model, 2, rbuffer, op_handle)
356 amin=rbuffer(1)
357 amax=rbuffer(2)
358!
359! Unpack read data.
360!
361 ic=0
362 DO k=lbk,ubk
363 DO j=jmin,jmax
364 DO i=imin,imax
365 ic=ic+1
366 adat(i,j,k)=wrk(ic)
367 END DO
368 END DO
369 END DO
370 ELSE
371 exit_flag=2
372 ioerror=status
373 END IF
374 END IF
375
376# if defined MASKING && defined READ_WATER
377!
378!-----------------------------------------------------------------------
379! Parallel I/O: Read in tile data from requested field and scale it.
380! Processing water points only.
381!-----------------------------------------------------------------------
382!
383 IF (gtype.lt.0) THEN
384!
385! Set number of points to process, grid type switch, and offsets due
386! array packing into 1D array in column-major order.
387!
388 SELECT CASE (abs(mytype))
389 CASE (p3dvar)
390 ijlen=iobounds(ng)%xy_psi
391 wtype=p2dvar
392 ioff=0
393 joff=1
394 CASE (r3dvar, w3dvar)
395 ijlen=iobounds(ng)%xy_rho
396 wtype=r2dvar
397 ioff=1
398 joff=0
399 CASE (u3dvar)
400 ijlen=iobounds(ng)%xy_u
401 wtype=u2dvar
402 ioff=0
403 joff=0
404 CASE (v3dvar)
405 ijlen=iobounds(ng)%xy_v
406 wtype=v2dvar
407 ioff=1
408 joff=1
409 CASE DEFAULT
410 ijlen=iobounds(ng)%xy_rho
411 wtype=r2dvar
412 ioff=1
413 joff=0
414 END SELECT
415
416 IF (lbk.eq.0) THEN
417 koff=0
418 ELSE
419 koff=1
420 END IF
421
422 npts=ijlen*klen
423 ijsize=isize*jsize
424!
425! Allocate scratch work arrays.
426!
427 IF (.not.allocated(a2d)) THEN
428 allocate ( a2d(ijsize) )
429 END IF
430 IF (.not.allocated(wrk)) THEN
431 allocate ( wrk(npts) )
432 wrk=inival
433 END IF
434!
435! Read in data: all parallel nodes read a segment of the 1D data.
436! Recall that water points are pack in the NetCDF file in a single
437! dimension.
438!
439 CALL tile_bounds_1d (ng, myrank, npts, istr, iend)
440
441 start(1)=istr
442 total(1)=iend-istr+1
443 start(2)=1
444 total(2)=tindex
445
446 status=nf90_get_var(ncid, ncvarid, wrk(istr:), start, total)
447!
448! Global reduction of work array. We need this because the packing
449! of the water point only affects the model tile partition.
450!
451 IF (status.eq.nf90_noerr) THEN
452 CALL mp_collect (ng, model, npts, inival, wrk)
453!
454! Scale read data and process fill values, if any. Compute minimum
455! and maximum values.
456!
457 amin=spval
458 amax=-spval
459 DO i=1,npts
460 IF (abs(wrk(i)).ge.abs(aspval)) THEN
461 wrk(i)=0.0_r8 ! set _FillValue to zero
462 ELSE
463 wrk(i)=ascl*(afactor*wrk(i)+aoffset)
464 amin=min(amin,wrk(i))
465 amax=max(amax,wrk(i))
466 END IF
467 END DO
468 IF ((abs(amin).ge.abs(aspval)).and. &
469 & (abs(amax).ge.abs(aspval))) THEN
470 amin=0.0_r8 ! the entire data is all
471 amax=0.0_r8 ! field value, _FillValue
472 END IF
473!
474! Unpack read data. This is tricky in parallel I/O. The cheapeast
475! thing to do is reconstruct a packed 2D global array and then select
476! the appropriate values for the tile.
477!
478 DO k=lbk,ubk
479 kc=(k-koff)*ijlen
480 a2d=inival
481 DO np=1,ijlen
482 ij=scalars(ng)%IJwater(np,wtype)
483 a2d(ij)=wrk(np+kc)
484 END DO
485 DO j=jmin,jmax
486 jc=(j-joff)*isize
487 DO i=imin,imax
488 ij=i+ioff+jc
489 adat(i,j,k)=a2d(ij)
490 END DO
491 END DO
492 END DO
493 ELSE
494 exit_flag=2
495 ioerror=status
496 END IF
497 END IF
498# endif
499!
500!-----------------------------------------------------------------------
501! Deallocate scratch work vector.
502!-----------------------------------------------------------------------
503!
504# if defined MASKING && defined READ_WATER
505 IF (allocated(a2d)) THEN
506 deallocate (a2d)
507 END IF
508# endif
509
510 IF (allocated(wrk)) THEN
511 deallocate (wrk)
512 END IF
513!
514 RETURN
515 END FUNCTION nf90_fread3d
516
517#else
518
519!
520!***********************************************************************
521 FUNCTION nf90_fread3d (ng, model, ncname, ncid, &
522 & ncvname, ncvarid, &
523 & tindex, gtype, Vsize, &
524 & LBi, UBi, LBj, UBj, LBk, UBk, &
525 & Ascl, Amin, Amax, &
526# ifdef MASKING
527 & Amask, &
528# endif
529 & Adat, &
530 & checksum) RESULT (status)
531!***********************************************************************
532!
533 USE mod_netcdf
534!
535# ifdef DISTRIBUTE
536 USE distribute_mod, ONLY : mp_bcasti
537# ifdef INLINE_2DIO
538 USE distribute_mod, ONLY : mp_scatter2d
539# else
540 USE distribute_mod, ONLY : mp_scatter3d
541# endif
542# endif
543!
544! Imported variable declarations.
545!
546 integer, intent(in) :: ng, model, ncid, ncvarid, tindex, gtype
547 integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk
548 integer, intent(in) :: vsize(4)
549!
550 integer(i8b), intent(out), optional :: checksum
551!
552 real(dp), intent(in) :: ascl
553 real(r8), intent(out) :: amin
554 real(r8), intent(out) :: amax
555!
556 character (len=*), intent(in) :: ncname
557 character (len=*), intent(in) :: ncvname
558!
559# ifdef ASSUMED_SHAPE
560# ifdef MASKING
561 real(r8), intent(in) :: amask(lbi:,lbj:)
562# endif
563 real(r8), intent(out) :: adat(lbi:,lbj:,lbk:)
564# else
565# ifdef MASKING
566 real(r8), intent(in) :: amask(lbi:ubi,lbj:ubj)
567# endif
568 real(r8), intent(out) :: adat(lbi:ubi,lbj:ubj,lbk:ubk)
569# endif
570!
571! Local variable declarations.
572!
573 logical :: lchecksum
574 logical, dimension(3) :: foundit
575!
576 integer :: i, j, k, ic, npts, nwpts, status, wtype
577 integer :: is, ie, js, je
578 integer :: imin, imax, jmin, jmax, koff
579 integer :: ilen, jlen, klen, ijlen
580 integer :: cgrid, mytype, ghost
581# ifdef DISTRIBUTE
582 integer :: nghost
583# endif
584 integer, dimension(4) :: start, total
585!
586 real(r8) :: afactor, aoffset, aspval
587
588 real(r8), dimension(3) :: attvalue
589!
590 real(r8), allocatable :: cwrk(:) ! used for checksum
591
592# if defined INLINE_2DIO && defined DISTRIBUTE
593 real(r8), dimension(2+(Lm(ng)+2)*(Mm(ng)+2)) :: wrk
594# else
595 real(r8), dimension(2+(Lm(ng)+2)*(Mm(ng)+2)*(UBk-LBk+1)) :: wrk
596# endif
597!
598 character (len=12), dimension(3) :: attname
599
600 character (len=*), parameter :: myfile = &
601 & __FILE__//", nf90_fread3d"
602!
603!-----------------------------------------------------------------------
604! Set starting and ending indices to process.
605!-----------------------------------------------------------------------
606!
607! Set global (interior plus boundary) starting and ending grid cell
608! indices in the I- and J-directions according to staggered C-grid
609! classification.
610!
611 mytype=gtype
612
613 SELECT CASE (abs(mytype))
614 CASE (p2dvar, p3dvar)
615 cgrid=1 ! PSI-points
616 is=iobounds(ng)%ILB_psi
617 ie=iobounds(ng)%IUB_psi
618 js=iobounds(ng)%JLB_psi
619 je=iobounds(ng)%JUB_psi
620 CASE (r2dvar, r3dvar, w3dvar)
621 cgrid=2 ! RHO-points
622 is=iobounds(ng)%ILB_rho
623 ie=iobounds(ng)%IUB_rho
624 js=iobounds(ng)%JLB_rho
625 je=iobounds(ng)%JUB_rho
626 CASE (u2dvar, u3dvar)
627 cgrid=3 ! U-points
628 is=iobounds(ng)%ILB_u
629 ie=iobounds(ng)%IUB_u
630 js=iobounds(ng)%JLB_u
631 je=iobounds(ng)%JUB_u
632 CASE (v2dvar, v3dvar)
633 cgrid=4 ! V-points
634 is=iobounds(ng)%ILB_v
635 ie=iobounds(ng)%IUB_v
636 js=iobounds(ng)%JLB_v
637 je=iobounds(ng)%JUB_v
638 CASE DEFAULT
639 cgrid=2 ! RHO-points
640 is=iobounds(ng)%ILB_rho
641 ie=iobounds(ng)%IUB_rho
642 js=iobounds(ng)%JLB_rho
643 je=iobounds(ng)%JUB_rho
644 END SELECT
645
646 ilen=ie-is+1
647 jlen=je-js+1
648 klen=ubk-lbk+1
649 ijlen=ilen*jlen
650
651 IF (lbk.eq.0) THEN
652 koff=0
653 ELSE
654 koff=1
655 END IF
656!
657! Set the tile computational I- and J-bounds (no ghost points).
658!
659 ghost=0
660 imin=bounds(ng)%Imin(cgrid,ghost,myrank)
661 imax=bounds(ng)%Imax(cgrid,ghost,myrank)
662 jmin=bounds(ng)%Jmin(cgrid,ghost,myrank)
663 jmax=bounds(ng)%Jmax(cgrid,ghost,myrank)
664!
665! Check if the following attributes: "scale_factor", "add_offset", and
666! "_FillValue" are present in the input NetCDF variable:
667!
668! If the "scale_value" attribute is present, the data is multiplied by
669! this factor after reading.
670! If the "add_offset" attribute is present, this value is added to the
671! data after reading.
672! If both "scale_factor" and "add_offset" attributes are present, the
673! data are first scaled before the offset is added.
674! If the "_FillValue" attribute is present, the data having this value
675! is treated a missing and it is replaced with zero. This feature it is
676! usually related with the land/sea masking.
677!
678 attname(1)='scale_factor'
679 attname(2)='add_offset '
680 attname(3)='_FillValue '
681
682 CALL netcdf_get_fatt (ng, model, ncname, ncvarid, attname, &
683 & attvalue, foundit, &
684 & ncid = ncid)
685 IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
686 status=ioerror
687 RETURN
688 END IF
689
690 IF (.not.foundit(1)) THEN
691 afactor=1.0_r8
692 ELSE
693 afactor=attvalue(1)
694 END IF
695
696 IF (.not.foundit(2)) THEN
697 aoffset=0.0_r8
698 ELSE
699 aoffset=attvalue(2)
700 END IF
701
702 IF (.not.foundit(3)) THEN
703 aspval=spval_check
704 ELSE
705 aspval=attvalue(3)
706 END IF
707
708# ifdef DISTRIBUTE
709!
710! Set the number of tile ghost points, Nghost, to scatter in
711! distributed-memory applications. If Nghost=0, the ghost points
712! are not processed. They will be processed elsewhere by the
713! appropriate call to any of the routines in "mp_exchange.F".
714!
715# ifdef NO_READ_GHOST
716 nghost=0 ! no ghost points exchange
717# else
718 IF (model.eq.iadm) THEN
719 nghost=0 ! no ghost points exchange
720 ELSE
721 nghost=nghostpoints ! do ghost points exchange
722 END IF
723# endif
724# endif
725# if defined READ_WATER && defined MASKING
726!
727! If processing water points only, set number of points and type
728! switch.
729!
730 SELECT CASE (abs(mytype))
731 CASE (p3dvar)
732 npts=iobounds(ng)%xy_psi
733 wtype=p2dvar
734 CASE (r3dvar, w3dvar)
735 npts=iobounds(ng)%xy_rho
736 wtype=r2dvar
737 CASE (u3dvar)
738 npts=iobounds(ng)%xy_u
739 wtype=u2dvar
740 CASE (v3dvar)
741 npts=iobounds(ng)%xy_v
742 wtype=v2dvar
743 CASE DEFAULT
744 npts=iobounds(ng)%xy_rho
745 wtype=r2dvar
746 END SELECT
747 nwpts=(lm(ng)+2)*(mm(ng)+2)
748# if !(defined INLINE_2DIO && defined DISTRIBUTE)
749 npts=npts*klen
750# endif
751# endif
752!
753! Set NetCDF dimension counters for processing requested field.
754!
755 IF (mytype.gt.0) THEN
756 start(1)=1
757 total(1)=ilen
758 start(2)=1
759 total(2)=jlen
760 start(3)=1
761 total(3)=klen
762 start(4)=tindex
763 total(4)=1
764 npts=ijlen
765# if !(defined INLINE_2DIO && defined DISTRIBUTE)
766 npts=npts*klen
767# endif
768# if defined READ_WATER && defined MASKING
769 ELSE
770 start(1)=1
771 total(1)=npts
772 start(2)=1
773 total(2)=tindex
774# endif
775 END IF
776!
777! Initialize local array to avoid denormalized numbers. This
778! facilitates processing and debugging.
779!
780 wrk=0.0_r8
781!
782! Initialize checsum value.
783!
784 IF (PRESENT(checksum)) THEN
785 lchecksum=.true.
786 checksum=0_i8b
787 ELSE
788 lchecksum=.false.
789 END IF
790!
791!-----------------------------------------------------------------------
792! Serial I/O: Read in requested field and scale it.
793!-----------------------------------------------------------------------
794!
795 amin=spval
796 amax=-spval
797
798# if defined INLINE_2DIO && defined DISTRIBUTE
799!
800! If appropriate, process 3D data level by level to reduce memory
801! requirements.
802!
803 DO k=lbk,ubk
804 start(3)=k-koff+1
805 total(3)=1
806# endif
807 status=nf90_noerr
808 IF (inpthread) THEN
809 status=nf90_get_var(ncid, ncvarid, wrk, start, total)
810 IF (status.eq.nf90_noerr) THEN
811 DO i=1,npts
812 IF (abs(wrk(i)).ge.abs(aspval)) THEN
813 wrk(i)=0.0_r8 ! masked with _FillValue
814 ELSE
815 wrk(i)=ascl*(afactor*wrk(i)+aoffset)
816 amin=min(amin,wrk(i))
817 amax=max(amax,wrk(i))
818 END IF
819 END DO
820 IF ((abs(amin).ge.abs(aspval)).and. &
821 & (abs(amax).ge.abs(aspval))) THEN
822 amin=0.0_r8 ! the entire data is all
823 amax=0.0_r8 ! field value, _FillValue
824 END IF
825 END IF
826 END IF
827# ifdef DISTRIBUTE
828 CALL mp_bcasti (ng, model, status)
829# endif
830 IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
831 exit_flag=2
832 ioerror=status
833 RETURN
834 END IF
835!
836!-----------------------------------------------------------------------
837! Serial I/O: Unpack read field.
838!-----------------------------------------------------------------------
839
840# ifdef DISTRIBUTE
841!
842! Scatter read data over the distributed memory tiles.
843!
844# ifdef INLINE_2DIO
845 CALL mp_scatter2d (ng, model, lbi, ubi, lbj, ubj, &
846 & nghost, mytype, amin, amax, &
847# if defined READ_WATER && defined MASKING
848 & nwpts, scalars(ng)%IJwater(:,wtype), &
849# endif
850 & npts, wrk, adat(:,:,k))
851 END DO
852# else
853 CALL mp_scatter3d (ng, model, lbi, ubi, lbj, ubj, lbk, ubk, &
854 & nghost, mytype, amin, amax, &
855# if defined READ_WATER && defined MASKING
856 & nwpts, scalars(ng)%IJwater(:,wtype), &
857# endif
858 & npts, wrk, adat)
859# endif
860# else
861!
862! Unpack data into the global array: serial, serial with partitions,
863! and shared-memory applications.
864!
865 IF (mytype.gt.0) THEN
866 ic=0
867 DO k=lbk,ubk
868 DO j=js,je
869 DO i=is,ie
870 ic=ic+1
871 adat(i,j,k)=wrk(ic)
872 END DO
873 END DO
874 END DO
875# if defined MASKING || defined READ_WATER
876 ELSE
877 ic=0
878 DO k=lbk,ubk
879 DO j=js,je
880 DO i=is,ie
881 IF (amask(i,j).gt.0.0_r8) THEN
882 ic=ic+1
883 adat(i,j,k)=wrk(ic)
884 ELSE
885 adat(i,j,k)=0.0_r8
886 END IF
887 END DO
888 END DO
889 END DO
890# endif
891 END IF
892# endif
893!
894!-----------------------------------------------------------------------
895! If requested, compute data checksum value.
896!-----------------------------------------------------------------------
897!
898 IF (lchecksum) THEN
899# ifdef DISTRIBUTE
900 npts=(imax-imin+1)*(jmax-jmin+1)*(ubk-lbk+1)
901 IF (.not.allocated(cwrk)) allocate ( cwrk(npts) )
902 cwrk=pack(adat(imin:imax, jmin:jmax, lbk:ubk), .true.)
903 CALL get_hash (cwrk, npts, checksum, .true.)
904# else
905 npts=(ie-is+1)*(je-js+1)*(ubk-lbk+1)
906 IF (.not.allocated(cwrk)) allocate ( cwrk(npts) )
907 cwrk=pack(adat(is:ie, js:je, lbk:ubk), .true.)
908 CALL get_hash (cwrk, npts, checksum)
909# endif
910 IF (allocated(cwrk)) deallocate (cwrk)
911 END IF
912!
913 RETURN
914 END FUNCTION nf90_fread3d
915#endif
916
917#if defined PIO_LIB && defined DISTRIBUTE
918!
919!***********************************************************************
920 FUNCTION pio_fread3d (ng, model, ncname, pioFile, &
921 & ncvname, pioVar, &
922 & tindex, pioDesc, Vsize, &
923 & LBi, UBi, LBj, UBj, LBk, UBk, &
924 & Ascl, Amin, Amax, &
925# ifdef MASKING
926 & Amask, &
927# endif
928 & Adat, &
929 & checksum) RESULT (status)
930!***********************************************************************
931!
933!
934 USE distribute_mod, ONLY : mp_reduce
935!
936! Imported variable declarations.
937!
938 integer, intent(in) :: ng, model, tindex
939 integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk
940 integer, intent(in) :: vsize(4)
941!
942 integer(i8b), intent(out), optional :: checksum
943!
944 real(dp), intent(in) :: ascl
945 real(r8), intent(out) :: amin
946 real(r8), intent(out) :: amax
947!
948 character (len=*), intent(in) :: ncname
949 character (len=*), intent(in) :: ncvname
950!
951# ifdef ASSUMED_SHAPE
952# ifdef MASKING
953 real(r8), intent(in) :: amask(lbi:,lbj:)
954# endif
955 real(r8), intent(out) :: adat(lbi:,lbj:,lbk:)
956# else
957# ifdef MASKING
958 real(r8), intent(in) :: amask(lbi:ubi,lbj:ubj)
959# endif
960 real(r8), intent(out) :: adat(lbi:ubi,lbj:ubj,lbk:ubk)
961# endif
962!
963 TYPE (file_desc_t), intent(inout) :: piofile
964 TYPE (io_desc_t), intent(inout) :: piodesc
965 TYPE (my_vardesc), intent(inout) :: piovar
966!
967! Local variable declarations.
968!
969 logical :: lchecksum
970 logical, dimension(3) :: foundit
971!
972 integer :: i, j, k, npts, status
973 integer :: is, ie, js, je
974 integer :: imin, imax, jmin, jmax
975 integer :: cgrid, ghost, dkind, gtype
976
977 integer, dimension(4) :: start, total
978!
979 real(r8) :: afactor, aoffset, aspval, avalue
980 real(r8) :: my_amin, my_amax
981
982 real(r8), dimension(3) :: attvalue
983 real(r8), dimension(2) :: rbuffer
984!
985 real(r4), pointer :: awrk4(:,:,:) ! single precision
986 real(r8), pointer :: awrk8(:,:,:) ! double precision
987 real(r8), allocatable :: cwrk(:) ! used for checksum
988!
989 character (len=12), dimension(3) :: attname
990 character (len= 3), dimension(2) :: op_handle
991
992 character (len=*), parameter :: myfile = &
993 & __FILE__//", pio_fread3d"
994!
995!-----------------------------------------------------------------------
996! Set starting and ending indices to process.
997!-----------------------------------------------------------------------
998!
999 status=pio_noerr
1000 amin=spval
1001 amax=-spval
1002 my_amin=spval
1003 my_amax=-spval
1004!
1005! Set global (interior plus boundary) starting and ending grid cell
1006! indices in the I- and J-directions according to staggered C-grid
1007! classification.
1008!
1009 dkind=piovar%dkind
1010 gtype=piovar%gtype
1011!
1012 SELECT CASE (abs(gtype))
1013 CASE (p2dvar, p3dvar)
1014 cgrid=1 ! PSI-points
1015 is=iobounds(ng)%ILB_psi
1016 ie=iobounds(ng)%IUB_psi
1017 js=iobounds(ng)%JLB_psi
1018 je=iobounds(ng)%JUB_psi
1019 CASE (b3dvar, l3dvar, r2dvar, r3dvar, w3dvar)
1020 cgrid=2 ! RHO-points
1021 is=iobounds(ng)%ILB_rho
1022 ie=iobounds(ng)%IUB_rho
1023 js=iobounds(ng)%JLB_rho
1024 je=iobounds(ng)%JUB_rho
1025 CASE (u2dvar, u3dvar)
1026 cgrid=3 ! U-points
1027 is=iobounds(ng)%ILB_u
1028 ie=iobounds(ng)%IUB_u
1029 js=iobounds(ng)%JLB_u
1030 je=iobounds(ng)%JUB_u
1031 CASE (v2dvar, v3dvar)
1032 cgrid=4 ! V-points
1033 is=iobounds(ng)%ILB_v
1034 ie=iobounds(ng)%IUB_v
1035 js=iobounds(ng)%JLB_v
1036 je=iobounds(ng)%JUB_v
1037 CASE DEFAULT
1038 cgrid=2 ! RHO-points
1039 is=iobounds(ng)%ILB_rho
1040 ie=iobounds(ng)%IUB_rho
1041 js=iobounds(ng)%JLB_rho
1042 je=iobounds(ng)%JUB_rho
1043 END SELECT
1044!
1045! Set the tile computational I- and J-bounds (no ghost points).
1046!
1047 ghost=0
1048 imin=bounds(ng)%Imin(cgrid,ghost,myrank)
1049 imax=bounds(ng)%Imax(cgrid,ghost,myrank)
1050 jmin=bounds(ng)%Jmin(cgrid,ghost,myrank)
1051 jmax=bounds(ng)%Jmax(cgrid,ghost,myrank)
1052!
1053! Check if the following attributes: "scale_factor", "add_offset", and
1054! "_FillValue" are present in the input NetCDF variable:
1055!
1056! If the "scale_value" attribute is present, the data is multiplied by
1057! this factor after reading.
1058! If the "add_offset" attribute is present, this value is added to the
1059! data after reading.
1060! If both "scale_factor" and "add_offset" attributes are present, the
1061! data are first scaled before the offset is added.
1062! If the "_FillValue" attribute is present, the data having this value
1063! is treated a missing and it is replaced with zero. This feature it is
1064! usually related with the land/sea masking.
1065!
1066 attname(1)='scale_factor'
1067 attname(2)='add_offset '
1068 attname(3)='_FillValue '
1069
1070 CALL pio_netcdf_get_fatt (ng, model, ncname, piovar%vd, attname, &
1071 & attvalue, foundit, &
1072 & piofile = piofile)
1073 IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
1074 status=ioerror
1075 RETURN
1076 END IF
1077
1078 IF (.not.foundit(1)) THEN
1079 afactor=1.0_r8
1080 ELSE
1081 afactor=attvalue(1)
1082 END IF
1083
1084 IF (.not.foundit(2)) THEN
1085 aoffset=0.0_r8
1086 ELSE
1087 aoffset=attvalue(2)
1088 END IF
1089
1090 IF (.not.foundit(3)) THEN
1091 aspval=spval_check
1092 ELSE
1093 aspval=attvalue(3)
1094 END IF
1095!
1096! Initialize checsum value.
1097!
1098 IF (PRESENT(checksum)) THEN
1099 lchecksum=.true.
1100 checksum=0_i8b
1101 ELSE
1102 lchecksum=.false.
1103 END IF
1104!
1105!-----------------------------------------------------------------------
1106! Read in requested field and scale it.
1107!-----------------------------------------------------------------------
1108!
1109! Allocate and initialize local array used for reading. The local array
1110! needs to be of the same precision as "A" and its IO decomposition
1111! descriptor "pioDesc".
1112!
1113 IF (dkind.eq.pio_double) THEN ! double precision
1114 IF (.not.associated(awrk8)) THEN
1115 allocate ( awrk8(imin:imax, jmin:jmax, lbk:ubk) )
1116 END IF
1117 awrk8=0.0_r8
1118 ELSE ! single precision
1119 IF (.not.associated(awrk4)) THEN
1120 allocate ( awrk4(imin:imax, jmin:jmax, lbk:ubk) )
1121 END IF
1122 awrk4=0.0_r4
1123 END IF
1124!
1125! Set unlimited time dimension record to read, if any.
1126!
1127 IF (tindex.gt.0) THEN
1128 CALL pio_setframe (piofile, &
1129 & piovar%vd, &
1130 & int(tindex, kind=pio_offset_kind))
1131 END IF
1132!
1133! Read in double precision data from NetCDF file.
1134!
1135 IF (dkind.eq.pio_double) THEN
1136 CALL pio_read_darray (piofile, &
1137 & piovar%vd, &
1138 & piodesc, &
1139 & awrk8(imin:,jmin:,lbk:), &
1140 & status)
1141!
1142 DO k=lbk,ubk
1143 DO j=jmin,jmax
1144 DO i=imin,imax
1145 IF (abs(awrk8(i,j,k)).ge.abs(aspval)) THEN
1146 adat(i,j,k)=0.0_r8 ! masked with _FillValue
1147 ELSE
1148 avalue=ascl*(afactor*awrk8(i,j,k)+aoffset)
1149 adat(i,j,k)=avalue
1150 my_amin=min(my_amin,avalue)
1151 my_amax=max(my_amax,avalue)
1152 END IF
1153 END DO
1154 END DO
1155 END DO
1156 IF (associated(awrk8)) deallocate (awrk8)
1157!
1158! Read in single precision data from NetCDF file.
1159!
1160 ELSE
1161 CALL pio_read_darray (piofile, &
1162 & piovar%vd, &
1163 & piodesc, &
1164 & awrk4(imin:,jmin:,lbk:), &
1165 & status)
1166!
1167 DO k=lbk,ubk
1168 DO j=jmin,jmax
1169 DO i=imin,imax
1170 IF (abs(awrk4(i,j,k)).ge.abs(aspval)) THEN
1171 adat(i,j,k)=0.0_r8 ! masked with _FillValue
1172 ELSE
1173 avalue=real(ascl*(afactor*awrk4(i,j,k)+aoffset),r8)
1174 adat(i,j,k)=avalue
1175 my_amin=real(min(my_amin,avalue),r8)
1176 my_amax=real(max(my_amax,avalue),r8)
1177 END IF
1178 END DO
1179 END DO
1180 END DO
1181 IF (associated(awrk4)) deallocate (awrk4)
1182 END IF
1183!
1184! If requested, compute input data checksum value.
1185!
1186 IF (lchecksum) THEN
1187 npts=(imax-imin+1)*(jmax-jmin+1)*(ubk-lbk+1)
1188 IF (.not.allocated(cwrk)) allocate ( cwrk(npts) )
1189 cwrk=pack(adat(imin:imax, jmin:jmax, lbk:ubk), .true.)
1190 CALL get_hash (cwrk, npts, checksum, .true.)
1191 IF (allocated(cwrk)) deallocate (cwrk)
1192 END IF
1193!
1194! Compute global minimum and maximum values.
1195!
1196 rbuffer(1)=my_amin
1197 rbuffer(2)=my_amax
1198 op_handle(1)='MIN'
1199 op_handle(2)='MAX'
1200 CALL mp_reduce (ng, model, 2, rbuffer, op_handle)
1201 amin=rbuffer(1)
1202 amax=rbuffer(2)
1203!
1204 IF ((abs(amin).ge.abs(spval)).and. &
1205 & (abs(amax).ge.abs(spval))) THEN
1206 amin=0.0_r8 ! the entire data is all
1207 amax=0.0_r8 ! field value, _FillValue
1208 END IF ! and was zeroth out
1209!
1210 RETURN
1211 END FUNCTION pio_fread3d
1212#endif
1213 END MODULE nf_fread3d_mod
subroutine mp_scatter3d(ng, model, lbi, ubi, lbj, ubj, lbk, ubk, nghost, gtype, amin, amax, nwpts, ij_water, npts, a, awrk)
subroutine mp_scatter2d(ng, model, lbi, ubi, lbj, ubj, nghost, gtype, amin, amax, nwpts, ij_water, npts, a, awrk)
subroutine, public tile_bounds_1d(ng, tile, imax, istr, iend)
Definition get_bounds.F:921
subroutine, public get_hash(a, asize, hash, lreduce)
Definition get_hash.F:72
integer ioerror
logical inpthread
type(t_bounds), dimension(:), allocatable bounds
Definition mod_param.F:232
integer, parameter b3dvar
Definition mod_param.F:725
integer, parameter r3dvar
Definition mod_param.F:721
type(t_iobounds), dimension(:), allocatable iobounds
Definition mod_param.F:282
integer nghostpoints
Definition mod_param.F:710
integer, parameter iadm
Definition mod_param.F:665
integer, parameter u3dvar
Definition mod_param.F:722
integer, dimension(:), allocatable lm
Definition mod_param.F:455
integer, parameter u2dvar
Definition mod_param.F:718
integer, parameter w3dvar
Definition mod_param.F:724
integer, parameter p2dvar
Definition mod_param.F:716
integer, dimension(:), allocatable mm
Definition mod_param.F:456
integer, parameter r2dvar
Definition mod_param.F:717
integer, parameter l3dvar
Definition mod_param.F:726
integer, parameter v2dvar
Definition mod_param.F:719
integer, parameter p3dvar
Definition mod_param.F:720
integer, parameter v3dvar
Definition mod_param.F:723
real(dp), parameter spval
real(dp), parameter spval_check
integer exit_flag
type(t_scalars), dimension(:), allocatable scalars
Definition mod_scalars.F:65
integer noerror
integer function pio_fread3d(ng, model, ncname, piofile, ncvname, piovar, tindex, piodesc, vsize, lbi, ubi, lbj, ubj, lbk, ubk, ascl, amin, amax, amask, adat, checksum)
Definition nf_fread3d.F:930
integer function nf90_fread3d(ng, model, ncname, ncid, ncvname, ncvarid, tindex, gtype, vsize, lbi, ubi, lbj, ubj, lbk, ubk, ascl, amin, amax, amask, adat, checksum)
Definition nf_fread3d.F:91
logical function, public founderror(flag, noerr, line, routine)
Definition strings.F:52