nf_fread4d.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE nf_fread4d_mod
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This function reads in a generic floating point 4D array from an    !
!  input NetCDF file.                                                  !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Nested grid number (integer)                          !
!     model      Calling model identifier (integer)                    !
!     ncname     NetCDF file name (string)                             !
!     ncid       NetCDF file ID (integer)                              !
#if defined PIO_LIB && defined DISTRIBUTE
!or   pioFile    PIO file descriptor structure, TYPE(file_desc_t)      !
!                  pioFile%fh         file handler                     !
!                  pioFile%iosystem   IO system descriptor (struct)    !
#endif
!     ncvname    NetCDF variable name (string)                         !
!     ncvarid    NetCDF variable ID (integer)                          !
#if defined PIO_LIB && defined DISTRIBUTE
!or   pioVar     PIO variable descriptor structure, TYPE(My_VarDesc)   !
!                  pioVar%vd     variable descriptor TYPE(Var_Desc_t)  !
!                  pioVar%dkind  variable data kind                    !
!                  pioVar%gtype  variable C-gridtype                   !
#endif
!     tindex     NetCDF time record index to read (integer)            !
!     gtype      C-grid type (integer)                                 !
#if defined PIO_LIB && defined DISTRIBUTE
!or   pioDesc    IO data decomposition descriptor, TYPE(IO_desc_t)     !
#endif
!     Vsize      Variable dimensions in NetCDF file (integer 1D array) !
!     LBi        I-dimension Lower bound (integer)                     !
!     UBi        I-dimension Upper bound (integer)                     !
!     LBj        J-dimension Lower bound (integer)                     !
!     UBj        J-dimension Upper bound (integer)                     !
!     LBk        K-dimension Lower bound (integer)                     !
!     UBk        K-dimension Upper bound (integer)                     !
!     LBt        Time-dimension Lower bound (integer)                  !
!     UBt        Time-dimension Upper bound (integer)                  !
!     Ascl       Factor to scale field after reading (real).           !
!     Amask      Land/Sea mask, if any (real 4D array)                 !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     Amin       Field minimum value (real)                            !
!     Amax       Field maximum value (real)                            !
!     Adat       Field to read in (real 4D array)                      !
!     checksum   Field checksum value (32-bit integer; OPTIONAL)       !
!     status     Result error flag (integer)                           !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_iounits
      USE mod_ncparam
      USE mod_scalars
!
      USE get_hash_mod,   ONLY : get_hash
      USE strings_mod,    ONLY : founderror
!
      implicit none
!
      INTERFACE nf_fread4d
        MODULE PROCEDURE nf90_fread4d
#if defined PIO_LIB && defined DISTRIBUTE
        MODULE PROCEDURE pio_fread4d
#endif
      END INTERFACE nf_fread4d
!
      CONTAINS
 
#if defined PARALLEL_IO && defined DISTRIBUTE
!
!***********************************************************************
      FUNCTION nf90_fread4d (ng, model, ncname, ncid,                   &
     &                       ncvname, ncvarid,                          &
     &                       tindex, gtype, Vsize,                      &
     &                       LBi, UBi, LBj, UBj, LBk, UBk, LBt, UBt,    &
     &                       Ascl, Amin, Amax,                          &
# ifdef MASKING
     &                       Amask,                                     &
# endif
     &                       Adat,                                      &
     &                       checksum) RESULT (status)
!***********************************************************************
!
      USE mod_netcdf
!
      USE distribute_mod, ONLY : mp_bcasti, mp_reduce
# if defined MASKING && defined READ_WATER
      USE distribute_mod, ONLY : mp_collect
# endif
      USE get_bounds_mod, ONLY : tile_bounds_1d
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, model, ncid, ncvarid, tindex, gtype
      integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt
      integer, intent(in) :: vsize(4)
!
      integer(i8b), intent(out), optional :: checksum
!
      real(dp), intent(in)  :: ascl
      real(r8), intent(out) :: amin
      real(r8), intent(out) :: amax
!
      character (len=*), intent(in) :: ncname
      character (len=*), intent(in) :: ncvname
!
# ifdef ASSUMED_SHAPE
#  ifdef MASKING
      real(r8), intent(in) :: amask(lbi:,lbj:)
#  endif
      real(r8), intent(out) :: adat(lbi:,lbj:,lbk:,lbt:)
# else
#  ifdef MASKING
      real(r8), intent(in) :: amask(lbi:ubi,lbj:ubj)
#  endif
      real(r8), intent(out) :: adat(lbi:ubi,lbj:ubj,lbk:ubk,lbt:ubt)
# endif
!
!  Local variable declarations.
!
      logical :: lchecksum
      logical, dimension(3) :: foundit
!
      integer :: i, ic, ij, j, jc, k, kc, l, lc, np, npts
      integer :: imin, imax, isize, jmin, jmax, jsize, ijsize
      integer :: istr, iend
      integer :: ioff, joff, koff, loff
      integer :: ilen, jlen, klen, llen, ijlen
      integer :: cgrid, mytype, ghost, status, wtype
 
      integer, dimension(5) :: start, total
!
      real(r8) :: afactor, aoffset, aspval
 
      real(r8), parameter :: inival= 0.0_r8
 
      real(r8), dimension(2) :: rbuffer
      real(r8), dimension(3) :: attvalue
 
# if defined MASKING && defined READ_WATER
      real(r8), allocatable :: a2d(:)
# endif
      real(r8), allocatable :: wrk(:)
!
      character (len= 3), dimension(2) :: op_handle
      character (len=12), dimension(3) :: attname
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", nf90_fread4d"
!
!-----------------------------------------------------------------------
!  Set starting and ending indices to process.
!-----------------------------------------------------------------------
!
      status=nf90_noerr
!
!  Set first and last grid point according to staggered C-grid
!  classification. Set the offsets for variables with starting
!  zero-index.  Recall the NetCDF does not support a zero-index.
!
!  Notice that (Imin,Jmin) and (Imax,Jmax) are the corner of the
!  computational tile. If ghost=0, ghost points are not processed.
!  They will be processed elsewhere by the appropriate call to any
!  of the routines in "mp_exchange.F".  If ghost=1, the ghost points
!  are read.
!
# ifdef NO_READ_GHOST
      ghost=0                        ! non-overlapping, no ghost points
# else
      IF (model.eq.iadm) THEN
        ghost=0                      ! non-overlapping, no ghost points
      ELSE
        ghost=1                      ! overlapping, read ghost points
      END IF
# endif
 
      mytype=gtype
 
      SELECT CASE (abs(mytype))
        CASE (p2dvar, p3dvar)
          cgrid=1
          isize=iobounds(ng)%xi_psi
          jsize=iobounds(ng)%eta_psi
        CASE (r2dvar, r3dvar, w3dvar)
          cgrid=2
          isize=iobounds(ng)%xi_rho
          jsize=iobounds(ng)%eta_rho
        CASE (u2dvar, u3dvar)
          cgrid=3
          isize=iobounds(ng)%xi_u
          jsize=iobounds(ng)%eta_u
        CASE (v2dvar, v3dvar)
          cgrid=4
          isize=iobounds(ng)%xi_v
          jsize=iobounds(ng)%eta_v
        CASE DEFAULT
          cgrid=2
          isize=iobounds(ng)%xi_rho
          jsize=iobounds(ng)%eta_rho
      END SELECT
 
      imin=bounds(ng)%Imin(cgrid,ghost,myrank)
      imax=bounds(ng)%Imax(cgrid,ghost,myrank)
      jmin=bounds(ng)%Jmin(cgrid,ghost,myrank)
      jmax=bounds(ng)%Jmax(cgrid,ghost,myrank)
 
      ilen=imax-imin+1
      jlen=jmax-jmin+1
      klen=ubk-lbk+1
      llen=ubt-lbt+1
!
!  Check if the following attributes: "scale_factor", "add_offset", and
!  "_FillValue" are present in the input NetCDF variable:
!
!  If the "scale_value" attribute is present, the data is multiplied by
!  this factor after reading.
!  If the "add_offset" attribute is present, this value is added to the
!  data after reading.
!  If both "scale_factor" and "add_offset" attributes are present, the
!  data are first scaled before the offset is added.
!  If the "_FillValue" attribute is present, the data having this value
!  is treated as missing and it is replaced with zero. This feature it
!  is usually related with the land/sea masking.
!
      attname(1)='scale_factor'
      attname(2)='add_offset  '
      attname(3)='_FillValue  '
 
      CALL netcdf_get_fatt (ng, model, ncname, ncvarid, attname,        &
     &                      attvalue, foundit,                          &
     &                      ncid = ncid)
      IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
        status=ioerror
        RETURN
      END IF
 
      IF (.not.foundit(1)) THEN
        afactor=1.0_r8
      ELSE
        afactor=attvalue(1)
      END IF
 
      IF (.not.foundit(2)) THEN
        aoffset=0.0_r8
      ELSE
        aoffset=attvalue(2)
      END IF
 
      IF (.not.foundit(3)) THEN
        aspval=spval_check
      ELSE
        aspval=attvalue(3)
      END IF
!
!  Initialize checsum value.
!
      IF (PRESENT(checksum)) THEN
        lchecksum=.true.
        checksum=0_i8b
      ELSE
        lchecksum=.false.
      END IF
!
!-----------------------------------------------------------------------
!  Parallel I/O: Read in tile data from requested field and scale it.
!                Processing both water and land points.
!-----------------------------------------------------------------------
!
      IF (gtype.gt.0) THEN
!
!  Set offsets due the NetCDF dimensions. Recall that some output
!  variables not always start at one.
!
        SELECT CASE (abs(mytype))
          CASE (p2dvar, p3dvar)
            ioff=0
            joff=0
          CASE (r2dvar, r3dvar, w3dvar)
            ioff=1
            joff=1
          CASE (u2dvar, u3dvar)
            ioff=0
            joff=1
          CASE (v2dvar, v3dvar)
            ioff=1
            joff=0
          CASE DEFAULT
            ioff=1
            joff=1
        END SELECT
 
        IF (lbk.eq.0) THEN
          koff=1
        ELSE
          koff=0
        END IF
 
        IF (lbt.eq.0) THEN
          loff=1
        ELSE
          loff=0
        END IF
 
        npts=ilen*jlen*klen*llen
!
!  Allocate scratch work array.
!
        IF (.not.allocated(wrk)) THEN
          allocate ( wrk(npts) )
          wrk=0.0_r8
        END IF
!
!  Read in data: all parallel nodes read their own tile data.
!
        start(1)=imin+ioff
        total(1)=ilen
        start(2)=jmin+joff
        total(2)=jlen
        start(3)=lbk+koff
        total(3)=klen
        start(4)=lbt+loff
        total(4)=llen
        start(5)=tindex
        total(5)=1
 
        status=nf90_get_var(ncid, ncvarid, wrk, start, total)
!
!  Scale read data and process fill values, if any.  Compute minimum
!  and maximum values.
!
        IF (status.eq.nf90_noerr) THEN
          amin=spval
          amax=-spval
          DO i=1,npts
            IF (abs(wrk(i)).ge.abs(aspval)) THEN
              wrk(i)=0.0_r8                   ! masked with _FillValue
            ELSE
              wrk(i)=ascl*(afactor*wrk(i)+aoffset)
              amin=min(amin,wrk(i))
              amax=max(amax,wrk(i))
            END IF
          END DO
          IF ((abs(amin).ge.abs(aspval)).and.                           &
     &        (abs(amax).ge.abs(aspval))) THEN
            amin=0.0_r8                       ! the entire data is all
            amax=0.0_r8                       ! field value, _FillValue
          END IF
!
!  Set minimum and maximum values: global reduction.
!
          rbuffer(1)=amin
          op_handle(1)='MIN'
          rbuffer(2)=amax
          op_handle(2)='MAX'
          CALL mp_reduce (ng, model, 2, rbuffer, op_handle)
          amin=rbuffer(1)
          amax=rbuffer(2)
!
!  Unpack read data.
!
          ic=0
          DO l=lbt,ubt
            DO k=lbk,ubk
              DO j=jmin,jmax
                DO i=imin,imax
                  ic=ic+1
                  adat(i,j,k,l)=wrk(ic)
                END DO
              END DO
            END DO
          END DO
        ELSE
          exit_flag=2
          ioerror=status
        END IF
      END IF
 
# if defined MASKING && defined READ_WATER
!
!-----------------------------------------------------------------------
!  Parallel I/O: Read in tile data from requested field and scale it.
!                Processing water points only.
!-----------------------------------------------------------------------
!
      IF (gtype.lt.0) THEN
!
!  Set number of points to process, grid type switch, and offsets due
!  array packing into 1D array in column-major order.
!
        SELECT CASE (abs(mytype))
          CASE (p3dvar)
            ijlen=iobounds(ng)%xy_psi
            wtype=p2dvar
            ioff=0
            joff=1
          CASE (r3dvar, w3dvar)
            ijlen=iobounds(ng)%xy_rho
            wtype=r2dvar
            ioff=1
            joff=0
          CASE (u3dvar)
            ijlen=iobounds(ng)%xy_u
            wtype=u2dvar
            ioff=0
            joff=0
          CASE (v3dvar)
            ijlen=iobounds(ng)%xy_v
            wtype=v2dvar
            ioff=1
            joff=1
          CASE DEFAULT
            ijlen=iobounds(ng)%xy_rho
            wtype=r2dvar
            ioff=1
            joff=0
        END SELECT
 
        IF (lbk.eq.0) THEN
          koff=0
        ELSE
          koff=1
        END IF
 
        IF (lbt.eq.0) THEN
          loff=1
        ELSE
          loff=0
        END IF
 
        npts=ijlen*klen*llen
        ijsize=isize*jsize
!
!  Allocate scratch work arrays.
!
        IF (.not.allocated(a2d)) THEN
          allocate ( a2d(ijsize) )
        END IF
        IF (.not.allocated(wrk)) THEN
          allocate ( wrk(npts) )
          wrk=inival
        END IF
!
!  Read in data: all parallel nodes read a segment of the 1D data.
!  Recall that water points are pack in the NetCDF file in a single
!  dimension.
!
        CALL tile_bounds_1d (ng, myrank, npts, istr, iend)
 
        start(1)=istr
        total(1)=iend-istr+1
        start(2)=1
        total(2)=tindex
 
        status=nf90_get_var(ncid, ncvarid, wrk(istr:), start, total)
!
!  Global reduction of work array.  We need this because the packing
!  of the water point only affects the model tile partition.
!
        IF (status.eq.nf90_noerr) THEN
          CALL mp_collect (ng, model, npts, inival, wrk)
!
!  Scale read data and process fill values, if any.  Compute minimum
!  and maximum values.
!
          amin=spval
          amax=-spval
          DO i=1,npts
            IF (abs(wrk(i)).ge.abs(aspval)) THEN
              wrk(i)=0.0_r8             ! set _FillValue to zero
            ELSE
              wrk(i)=ascl*(afactor*wrk(i)+aoffset)
              amin=min(amin,wrk(i))
              amax=max(amax,wrk(i))
            END IF
          END DO
          IF ((abs(amin).ge.abs(aspval)).and.                           &
     &        (abs(amax).ge.abs(aspval))) THEN
            amin=0.0_r8                       ! the entire data is all
            amax=0.0_r8                       ! field value, _FillValue
          END IF
!
!  Unpack read data.  This is tricky in parallel I/O.  The cheapeast
!  thing to do is reconstruct a packed 2D global array and then select
!  the appropriate values for the tile.
!
          DO l=lbt,ubt
            lc=(l-loff)*ijlen*klen
            DO k=lbk,ubk
              kc=(k-koff)*ijlen+lc
              a2d=inival
              DO np=1,ijlen
                ij=scalars(ng)%IJwater(np,wtype)
                a2d(ij)=wrk(np+kc)
              END DO
              DO j=jmin,jmax
                jc=(j-joff)*isize
                DO i=imin,imax
                  ij=i+ioff+jc
                  adat(i,j,k,l)=a2d(ij)
                END DO
              END DO
            END DO
          END DO
        ELSE
          exit_flag=2
          ioerror=status
        END IF
      END IF
# endif
!
!-----------------------------------------------------------------------
!  Deallocate scratch work vector.
!-----------------------------------------------------------------------
!
# if defined MASKING && defined READ_WATER
      IF (allocated(a2d)) THEN
        deallocate (a2d)
      END IF
# endif
 
      IF (allocated(wrk)) THEN
        deallocate (wrk)
      END IF
!
      RETURN
      END FUNCTION nf90_fread4d
 
#else
 
!
!***********************************************************************
      FUNCTION nf90_fread4d (ng, model, ncname, ncid,                   &
     &                       ncvname, ncvarid,                          &
     &                       tindex, gtype, Vsize,                      &
     &                       LBi, UBi, LBj, UBj, LBk, UBk, LBt, UBt,    &
     &                       Ascl, Amin, Amax,                          &
# ifdef MASKING
     &                       Amask,                                     &
# endif
     &                       Adat,                                      &
     &                       checksum) RESULT (status)
!***********************************************************************
!
      USE mod_netcdf
!
# ifdef DISTRIBUTE
      USE distribute_mod, ONLY : mp_bcasti
#  ifdef INLINE_2DIO
      USE distribute_mod, ONLY : mp_scatter2d
#  else
      USE distribute_mod, ONLY : mp_scatter3d
#  endif
# endif
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, model, ncid, ncvarid, tindex, gtype
      integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt
      integer, intent(in) :: vsize(4)
!
      integer(i8b), intent(out), optional :: checksum
!
      real(dp), intent(in)  :: ascl
      real(r8), intent(out) :: amin
      real(r8), intent(out) :: amax
!
      character (len=*), intent(in) :: ncname
      character (len=*), intent(in) :: ncvname
!
# ifdef ASSUMED_SHAPE
#  ifdef MASKING
      real(r8), intent(in) :: amask(lbi:,lbj:)
#  endif
      real(r8), intent(out) :: adat(lbi:,lbj:,lbk:,lbt:)
# else
#  ifdef MASKING
      real(r8), intent(in) :: amask(lbi:ubi,lbj:ubj)
#  endif
      real(r8), intent(out) :: adat(lbi:ubi,lbj:ubj,lbk:ubk,lbt:ubt)
# endif
!
!  Local variable declarations.
!
      logical :: lchecksum
      logical, dimension(3) :: foundit
!
      integer :: i, j, k, ic, fourth, npts, nwpts, status, wtype
      integer :: is, ie, js, je
      integer :: imin, imax, jmin, jmax
      integer :: ilen, jlen, klen, ijlen, koff, loff
      integer :: cgrid, mytype, ghost
# ifdef DISTRIBUTE
      integer :: nghost
# endif
      integer, dimension(5) :: start, total
!
      real(r8) :: afactor, aoffset, aspval
 
      real(r8), dimension(3) :: attvalue
!
      real(r8), allocatable :: cwrk(:)           ! used for checksum
 
# if defined INLINE_2DIO && defined DISTRIBUTE
      real(r8), dimension(2+(Lm(ng)+2)*(Mm(ng)+2)) :: wrk
# else
      real(r8), dimension(2+(Lm(ng)+2)*(Mm(ng)+2)*(UBk-LBk+1)) :: wrk
# endif
!
      character (len=12), dimension(3) :: attname
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", nf90_fread4d"
!
!-----------------------------------------------------------------------
!  Set starting and ending indices to process.
!-----------------------------------------------------------------------
!
      status=nf90_noerr
      amin=spval
      amax=-spval
!
!  Set global (interior plus boundary) starting and ending grid cell
!  indices in the I- and J-directions according to staggered C-grid
!  classification.
!
      mytype=gtype
 
      SELECT CASE (abs(mytype))
        CASE (p2dvar, p3dvar)
          cgrid=1                                         ! PSI-points
          is=iobounds(ng)%ILB_psi
          ie=iobounds(ng)%IUB_psi
          js=iobounds(ng)%JLB_psi
          je=iobounds(ng)%JUB_psi
        CASE (r2dvar, r3dvar, w3dvar)
          cgrid=2                                         ! RHO-points
          is=iobounds(ng)%ILB_rho
          ie=iobounds(ng)%IUB_rho
          js=iobounds(ng)%JLB_rho
          je=iobounds(ng)%JUB_rho
        CASE (u2dvar, u3dvar)
          cgrid=3                                         ! U-points
          is=iobounds(ng)%ILB_u
          ie=iobounds(ng)%IUB_u
          js=iobounds(ng)%JLB_u
          je=iobounds(ng)%JUB_u
        CASE (v2dvar, v3dvar)
          cgrid=4                                         ! V-points
          is=iobounds(ng)%ILB_v
          ie=iobounds(ng)%IUB_v
          js=iobounds(ng)%JLB_v
          je=iobounds(ng)%JUB_v
        CASE DEFAULT
          cgrid=2                                         ! RHO-points
          is=iobounds(ng)%ILB_rho
          ie=iobounds(ng)%IUB_rho
          js=iobounds(ng)%JLB_rho
          je=iobounds(ng)%JUB_rho
      END SELECT
 
      ilen=ie-is+1
      jlen=je-js+1
      klen=ubk-lbk+1
      ijlen=ilen*jlen
 
#ifdef INLINE_2DIO
      IF (lbk.eq.0) THEN
        koff=0
      ELSE
        koff=1
      END IF
#endif
 
      IF (lbt.eq.0) THEN
        loff=1
      ELSE
        loff=0
      END IF
!
!  Set the tile computational I- and J-bounds (no ghost points).
!
      ghost=0
      imin=bounds(ng)%Imin(cgrid,ghost,myrank)
      imax=bounds(ng)%Imax(cgrid,ghost,myrank)
      jmin=bounds(ng)%Jmin(cgrid,ghost,myrank)
      jmax=bounds(ng)%Jmax(cgrid,ghost,myrank)
!
!  Check if the following attributes: "scale_factor", "add_offset", and
!  "_FillValue" are present in the input NetCDF variable:
!
!  If the "scale_value" attribute is present, the data is multiplied by
!  this factor after reading.
!  If the "add_offset" attribute is present, this value is added to the
!  data after reading.
!  If both "scale_factor" and "add_offset" attributes are present, the
!  data are first scaled before the offset is added.
!  If the "_FillValue" attribute is present, the data having this value
!  is treated as missing and it is replaced with zero. This feature it
!  is usually related with the land/sea masking.
!
      attname(1)='scale_factor'
      attname(2)='add_offset  '
      attname(3)='_FillValue  '
 
      CALL netcdf_get_fatt (ng, model, ncname, ncvarid, attname,        &
     &                      attvalue, foundit,                          &
     &                      ncid = ncid)
      IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
        status=ioerror
        RETURN
      END IF
 
      IF (.not.foundit(1)) THEN
        afactor=1.0_r8
      ELSE
        afactor=attvalue(1)
      END IF
 
      IF (.not.foundit(2)) THEN
        aoffset=0.0_r8
      ELSE
        aoffset=attvalue(2)
      END IF
 
      IF (.not.foundit(3)) THEN
        aspval=spval_check
      ELSE
        aspval=attvalue(3)
      END IF
 
# ifdef DISTRIBUTE
!
!  Set the number of tile ghost points, Nghost, to scatter in
!  distributed-memory applications. If Nghost=0, the ghost points
!  are not processed.  They will be processed elsewhere by the
!  appropriate call to any of the routines in "mp_exchange.F".
!
#  ifdef NO_READ_GHOST
      nghost=0
#  else
      IF (model.eq.iadm) THEN
        nghost=0
      ELSE
        nghost=nghostpoints
      END IF
#  endif
# endif
# if defined READ_WATER && defined MASKING
!
!  If processing water points only, set number of points and type
!  switch.
!
      SELECT CASE (abs(mytype))
        CASE (p3dvar)
          npts=iobounds(ng)%xy_psi
          wtype=p2dvar
        CASE (r3dvar, w3dvar)
          npts=iobounds(ng)%xy_rho
          wtype=r2dvar
        CASE (u3dvar)
          npts=iobounds(ng)%xy_u
          wtype=u2dvar
        CASE (v3dvar)
          npts=iobounds(ng)%xy_v
          wtype=v2dvar
        CASE DEFAULT
          npts=iobounds(ng)%xy_rho
          wtype=r2dvar
      END SELECT
      nwpts=(lm(ng)+2)*(mm(ng)+2)
#  if !(defined INLINE_2DIO && defined DISTRIBUTE)
      npts=npts*klen
#  endif
# endif
!
!  Initialize local array to avoid denormalized numbers. This
!  facilitates processing and debugging.
!
      wrk=0.0_r8
!
!  Initialize checsum value.
!
      IF (PRESENT(checksum)) THEN
        lchecksum=.true.
        checksum=0_i8b
      ELSE
        lchecksum=.false.
      END IF
!
!-----------------------------------------------------------------------
!  Serial I/O: Read in requested field and scale it.
!-----------------------------------------------------------------------
!
!  Proccess data as 3D slides.
!
      DO fourth=lbt,ubt
        IF (mytype.gt.0) THEN
          start(1)=1
          total(1)=ilen
          start(2)=1
          total(2)=jlen
          start(3)=1
          total(3)=klen
          start(4)=fourth+loff
          total(4)=1
          start(5)=tindex
          total(5)=1
          npts=ijlen
# if !(defined INLINE_2DIO && defined DISTRIBUTE)
          npts=npts*klen
# endif
# if defined READ_WATER && defined MASKING
        ELSE
          start(1)=1+(fourth+loff-1)*npts
          total(1)=npts
          start(2)=1
          total(2)=tindex
# endif
        END IF
# if defined INLINE_2DIO && defined DISTRIBUTE
!
!  Process data as 2D slices (level by level) to reduce memory
!  requirements.
!
      DO k=lbk,ubk
        start(3)=k-koff+1
        total(3)=1
# endif
        status=nf90_noerr
        IF (inpthread) THEN
          status=nf90_get_var(ncid, ncvarid, wrk, start, total)
          IF (status.eq.nf90_noerr) THEN
            DO i=1,npts
              IF (abs(wrk(i)).ge.abs(aspval)) THEN
                wrk(i)=0.0_r8                 ! masked with _FillValue
              ELSE
                wrk(i)=ascl*(afactor*wrk(i)+aoffset)
                amin=min(amin,wrk(i))
                amax=max(amax,wrk(i))
              END IF
            END DO
            IF ((abs(amin).ge.abs(aspval)).and.                         &
     &          (abs(amax).ge.abs(aspval))) THEN
              amin=0.0_r8                     ! the entire data is all
              amax=0.0_r8                     ! field value, _FillValue
            END IF
          END IF
          IF (lchecksum) THEN
            CALL get_hash (wrk, npts, checksum)
          END IF
        END IF
# ifdef DISTRIBUTE
        CALL mp_bcasti (ng, model, status)
# endif
        IF (founderror(status, nf90_noerr, __line__, myfile)) THEN
          exit_flag=2
          ioerror=status
          RETURN
        END IF
!
!-----------------------------------------------------------------------
!  Serial I/O: Unpack read field.
!-----------------------------------------------------------------------
 
# ifdef DISTRIBUTE
!
!  Scatter read data over the distributed memory tiles.
!
#  ifdef INLINE_2DIO
        CALL mp_scatter2d (ng, model, lbi, ubi, lbj, ubj,               &
     &                     nghost, mytype, amin, amax,                  &
#   if defined READ_WATER && defined MASKING
     &                     nwpts, scalars(ng)%IJwater(:,wtype),         &
#   endif
     &                     npts, wrk, adat(:,:,k,fourth))
      END DO
#  else
        CALL mp_scatter3d (ng, model, lbi, ubi, lbj, ubj, lbk, ubk,     &
     &                     nghost, mytype, amin, amax,                  &
#   if defined READ_WATER && defined MASKING
     &                     nwpts, scalars(ng)%IJwater(:,wtype),         &
#   endif
     &                     npts, wrk, adat(:,:,:,fourth))
#  endif
# else
!
!  Unpack data into the global array: serial, serial with partitions,
!  and shared-memory applications.
!
        IF (mytype.gt.0) THEN
          ic=0
          DO k=lbk,ubk
            DO j=js,je
              DO i=is,ie
                ic=ic+1
                adat(i,j,k,fourth)=wrk(ic)
              END DO
            END DO
          END DO
#  if defined MASKING || defined READ_WATER
        ELSE
          ic=0
          DO k=lbk,ubk
            DO j=js,je
              DO i=is,ie
                IF (amask(i,j).gt.0.0_r8) THEN
                  ic=ic+1
                  adat(i,j,k,fourth)=wrk(ic)
                ELSE
                  adat(i,j,k,fourth)=0.0_r8
                END IF
              END DO
            END DO
          END DO
#  endif
        END IF
# endif
      END DO
!
!-----------------------------------------------------------------------
!  If requested, compute data checksum value.
!-----------------------------------------------------------------------
!
      IF (lchecksum) THEN
# ifdef DISTRIBUTE
        npts=(imax-imin+1)*(jmax-jmin+1)*(ubk-lbk+1)*(ubt-lbt+1)
        IF (.not.allocated(cwrk)) allocate ( cwrk(npts) )
        cwrk=pack(adat(imin:imax, jmin:jmax, lbk:ubk, lbt:ubt), .true.)
        CALL get_hash (cwrk, npts, checksum, .true.)
# else
        npts=(ie-is+1)*(je-js+1)*(ubk-lbk+1)*(ubt-lbt+1)
        IF (.not.allocated(cwrk)) allocate ( cwrk(npts) )
        cwrk=pack(adat(is:ie, js:je, lbk:ubk, lbt:ubt), .true.)
        CALL get_hash (cwrk, npts, checksum)
# endif
        IF (allocated(cwrk)) deallocate (cwrk)
      END IF
 
!
      RETURN
      END FUNCTION nf90_fread4d
#endif
 
#if defined PIO_LIB && defined DISTRIBUTE
!
!***********************************************************************
      FUNCTION pio_fread4d (ng, model, ncname, pioFile,                 &
     &                      ncvname, pioVar,                            &
     &                      tindex, pioDesc, Vsize,                     &
     &                      LBi, UBi, LBj, UBj, LBk, UBk, LBt, UBt,     &
     &                      Ascl, Amin, Amax,                           &
# ifdef MASKING
     &                      Amask,                                      &
# endif
     &                      Adat,                                       &
     &                      checksum) RESULT (status)
!***********************************************************************
!
      USE mod_pio_netcdf
!
      USE distribute_mod, ONLY : mp_reduce
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, model, tindex
      integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt
      integer, intent(in) :: vsize(4)
!
      integer(i8b), intent(out), optional :: checksum
!
      real(dp), intent(in)  :: ascl
      real(r8), intent(out) :: amin
      real(r8), intent(out) :: amax
!
      character (len=*), intent(in) :: ncname
      character (len=*), intent(in) :: ncvname
!
# ifdef ASSUMED_SHAPE
#  ifdef MASKING
      real(r8), intent(in) :: amask(lbi:,lbj:)
#  endif
      real(r8), intent(out) :: adat(lbi:,lbj:,lbk:,lbt:)
# else
#  ifdef MASKING
      real(r8), intent(in) :: amask(lbi:ubi,lbj:ubj)
#  endif
      real(r8), intent(out) :: adat(lbi:ubi,lbj:ubj,lbk:ubk,lbt:ubt)
# endif
!
      TYPE (file_desc_t), intent(inout) :: piofile
      TYPE (io_desc_t),   intent(inout) :: piodesc
      TYPE (my_vardesc),  intent(inout) :: piovar
!
!  Local variable declarations.
!
      logical :: lchecksum
      logical, dimension(3) :: foundit
!
      integer :: i, j, k, l, npts, status
      integer :: is, ie, js, je
      integer :: imin, imax, jmin, jmax
      integer :: cgrid, ghost, dkind, gtype
 
      integer, dimension(5) :: start, total
!
      real(r8) :: afactor, aoffset, aspval, avalue
      real(r8) :: my_amin, my_amax
 
      real(r8), dimension(3) :: attvalue
      real(r8), dimension(2) :: rbuffer
!
      real(r4), pointer :: awrk4(:,:,:,:)        ! single precision
      real(r8), pointer :: awrk8(:,:,:,:)        ! double precision
      real(r8), allocatable :: cwrk(:)           ! used for checksum
!
      character (len=12), dimension(3) :: attname
      character (len= 3), dimension(2) :: op_handle
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", pio_fread4d"
!
!-----------------------------------------------------------------------
!  Set starting and ending indices to process.
!-----------------------------------------------------------------------
!
      status=pio_noerr
      amin=spval
      amax=-spval
      my_amin=spval
      my_amax=-spval
!
!  Set global (interior plus boundary) starting and ending grid cell
!  indices in the I- and J-directions according to staggered C-grid
!  classification.
!
      dkind=piovar%dkind
      gtype=piovar%gtype
!
      SELECT CASE (abs(gtype))
        CASE (p2dvar, p3dvar)
          cgrid=1                                         ! PSI-points
          is=iobounds(ng)%ILB_psi
          ie=iobounds(ng)%IUB_psi
          js=iobounds(ng)%JLB_psi
          je=iobounds(ng)%JUB_psi
        CASE (l4dvar, r2dvar, r3dvar, w3dvar)
          cgrid=2                                         ! RHO-points
          is=iobounds(ng)%ILB_rho
          ie=iobounds(ng)%IUB_rho
          js=iobounds(ng)%JLB_rho
          je=iobounds(ng)%JUB_rho
        CASE (u2dvar, u3dvar)
          cgrid=3                                         ! U-points
          is=iobounds(ng)%ILB_u
          ie=iobounds(ng)%IUB_u
          js=iobounds(ng)%JLB_u
          je=iobounds(ng)%JUB_u
        CASE (v2dvar, v3dvar)
          cgrid=4                                         ! V-points
          is=iobounds(ng)%ILB_v
          ie=iobounds(ng)%IUB_v
          js=iobounds(ng)%JLB_v
          je=iobounds(ng)%JUB_v
        CASE DEFAULT
          cgrid=2                                         ! RHO-points
          is=iobounds(ng)%ILB_rho
          ie=iobounds(ng)%IUB_rho
          js=iobounds(ng)%JLB_rho
          je=iobounds(ng)%JUB_rho
      END SELECT
!
!  Set the tile computational I- and J-bounds (no ghost points).
!
      ghost=0
      imin=bounds(ng)%Imin(cgrid,ghost,myrank)
      imax=bounds(ng)%Imax(cgrid,ghost,myrank)
      jmin=bounds(ng)%Jmin(cgrid,ghost,myrank)
      jmax=bounds(ng)%Jmax(cgrid,ghost,myrank)
!
!  Check if the following attributes: "scale_factor", "add_offset", and
!  "_FillValue" are present in the input NetCDF variable:
!
!  If the "scale_value" attribute is present, the data is multiplied by
!  this factor after reading.
!  If the "add_offset" attribute is present, this value is added to the
!  data after reading.
!  If both "scale_factor" and "add_offset" attributes are present, the
!  data are first scaled before the offset is added.
!  If the "_FillValue" attribute is present, the data having this value
!  is treated as missing and it is replaced with zero. This feature it
!  is usually related with the land/sea masking.
!
      attname(1)='scale_factor'
      attname(2)='add_offset  '
      attname(3)='_FillValue  '
 
      CALL pio_netcdf_get_fatt (ng, model, ncname, piovar%vd, attname,  &
     &                          attvalue, foundit,                      &
     &                          piofile = piofile)
      IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
        status=ioerror
        RETURN
      END IF
 
      IF (.not.foundit(1)) THEN
        afactor=1.0_r8
      ELSE
        afactor=attvalue(1)
      END IF
 
      IF (.not.foundit(2)) THEN
        aoffset=0.0_r8
      ELSE
        aoffset=attvalue(2)
      END IF
 
      IF (.not.foundit(3)) THEN
        aspval=spval_check
      ELSE
        aspval=attvalue(3)
      END IF
!
!  Initialize checsum value.
!
      IF (PRESENT(checksum)) THEN
        lchecksum=.true.
        checksum=0_i8b
      ELSE
        lchecksum=.false.
      END IF
!
!-----------------------------------------------------------------------
!  Read in requested field and scale it.
!-----------------------------------------------------------------------
!
!  Allocate and initialize local array used for reading. The local array
!  needs to be of the same precision as "A" and its IO decomposition
!  descriptor "pioDesc".
!
      IF (dkind.eq.pio_double) THEN                 ! double precision
        IF (.not.associated(awrk8)) THEN
          allocate ( awrk8(imin:imax, jmin:jmax, lbk:ubk, lbt:ubt) )
        END IF
        awrk8=0.0_r8
      ELSE                                          ! single precision
        IF (.not.associated(awrk4)) THEN
          allocate ( awrk4(imin:imax, jmin:jmax, lbk:ubk, lbt:ubt) )
        END IF
        awrk4=0.0_r4
      END IF
!
!  Set unlimited time dimension record to write, if any.
!
      IF (tindex.gt.0) THEN
        CALL pio_setframe (piofile,                                     &
     &                     piovar%vd,                                   &
     &                     int(tindex, kind=pio_offset_kind))
      END IF
!
!  Read in double precision data from NetCDF file.
!
      IF (dkind.eq.pio_double) THEN
        CALL pio_read_darray (piofile,                                  &
     &                        piovar%vd,                                &
     &                        piodesc,                                  &
     &                        awrk8(imin:,jmin:,lbk:,lbt:),             &
     &                        status)
!
        DO l=lbt,ubt
          DO k=lbk,ubk
            DO j=jmin,jmax
              DO i=imin,imax
                IF (abs(awrk8(i,j,k,l)).ge.abs(aspval)) THEN
                  adat(i,j,k,l)=0.0_r8          ! masked with _FillValue
                ELSE
                  avalue=ascl*(afactor*awrk8(i,j,k,l)+aoffset)
                  adat(i,j,k,l)=avalue
                  my_amin=min(my_amin,avalue)
                  my_amax=max(my_amax,avalue)
                END IF
              END DO
            END DO
          END DO
        END DO
        IF (associated(awrk8)) deallocate (awrk8)
!
!  Read in and load single precision data from NetCDF file.
!
      ELSE
        CALL pio_read_darray (piofile,                                  &
     &                        piovar%vd,                                &
     &                        piodesc,                                  &
     &                        awrk4(imin:,jmin:,lbk:,lbt:),             &
     &                        status)
!
        DO l=lbt,ubt
          DO k=lbk,ubk
            DO j=jmin,jmax
              DO i=imin,imax
                IF (abs(awrk4(i,j,k,l)).ge.abs(aspval)) THEN
                  adat(i,j,k,l)=0.0_r8          ! masked with _FillValue
                ELSE
                  avalue=real(ascl*(afactor*awrk4(i,j,k,l)+aoffset),r8)
                  adat(i,j,k,l)=avalue
                  my_amin=real(min(my_amin,avalue),r8)
                  my_amax=real(max(my_amax,avalue),r8)
                END IF
              END DO
            END DO
          END DO
        END DO
        IF (associated(awrk4)) deallocate (awrk4)
      END IF
!
!  If requested, compute input data checksum value.
!  (TODO: parallel gathering of checksum)
!
      IF (lchecksum) THEN
        npts=(imax-imin+1)*(jmax-jmin+1)*(ubk-lbk+1)*(ubt-lbt+1)
        IF (.not.allocated(cwrk)) allocate ( cwrk(npts) )
        cwrk=pack(adat(imin:imax, jmin:jmax, lbk:ubk, lbt:ubt), .true.)
        CALL get_hash (cwrk, npts, checksum, .true.)
        IF (allocated(cwrk)) deallocate (cwrk)
      END IF
!
!  Compute global minimum and maximum values.
!
      rbuffer(1)=my_amin
      rbuffer(2)=my_amax
      op_handle(1)='MIN'
      op_handle(2)='MAX'
      CALL mp_reduce (ng, model, 2, rbuffer, op_handle)
      amin=rbuffer(1)
      amax=rbuffer(2)
!
      IF ((abs(amin).ge.abs(spval)).and.                                &
     &    (abs(amax).ge.abs(spval))) THEN
        amin=0.0_r8                           ! the entire data is all
        amax=0.0_r8                           ! field value, _FillValue
      END IF                                  ! and was zeroth out
!
      RETURN
      END FUNCTION pio_fread4d
#endif
      END MODULE nf_fread4d_mod