nf_fread2d_xtr.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE nf_fread2d_xtr_mod
 
#ifdef GRID_EXTRACT
!
!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 2D array from an    !
!  input NetCDF file for field extraction.                             !
!                                                                      !
!  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)                     !
!     Ascl       Factor to scale field after reading (real).           !
!     Amask      Land/Sea mask, if any (real 2D array)                 !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     Amin       Field minimum value (real)                            !
!     Amax       Field maximum value (real)                            !
!     Adat       Field to read in (real 2D array)                      !
!     checksum   Field checksum value (integer*8; OPTIONAL)            !
!     Lregrid    Field regrid interpolation switch (logical; OPTIONAL) !
!     status     Result error flag (integer)                           !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_grid
      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_fread2d_xtr
        MODULE PROCEDURE nf90_fread2d_xtr
# if defined PIO_LIB && defined DISTRIBUTE
        MODULE PROCEDURE pio_fread2d_xtr
# endif
      END INTERFACE nf_fread2d_xtr
!
      CONTAINS
 
# if defined PARALLEL_IO && defined DISTRIBUTE
!
!***********************************************************************
      FUNCTION nf90_fread2d_xtr (ng, model, ncname, ncid,               &
     &                           ncvname, ncvarid,                      &
     &                           tindex, gtype, Vsize,                  &
     &                           LBi, UBi, LBj, UBj,                    &
     &                           Ascl, Amin, Amax,                      &
#  ifdef MASKING
     &                           Amask,                                 &
#  endif
     &                           Adat,                                  &
     &                           checksum,                              &
     &                           Lregrid) RESULT (status)
!***********************************************************************
!
      USE mod_netcdf
!
      USE distribute_mod, ONLY : mp_bcasti, mp_collect, mp_reduce
      USE get_bounds_mod, ONLY : tile_bounds_1d, tile_bounds_2d
      USE regrid_mod,     ONLY : regrid_nf90
!
!  Imported variable declarations.
!
      logical, intent(out), optional :: Lregrid
!
      integer, intent(in) :: ng, model, ncid, ncvarid, tindex, gtype
      integer, intent(in) :: LBi, UBi, LBj, UBj
      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:)
#  else
#   ifdef MASKING
      real(r8), intent(in) :: Amask(LBi:UBi,LBj:UBj)
#   endif
      real(r8), intent(out) :: Adat(LBi:UBi,LBj:UBj)
#  endif
!
!  Local variable declarations.
!
      logical :: Lchecksum, interpolate
!
      logical, dimension(3) :: foundit
!
      integer :: i, ic, ij, j, jc, np, MyNpts, Npts
      integer :: Imin, Imax, Isize, Jmin, Jmax, Jsize, IJsize
      integer :: Istr, Iend, Jstr, Jend
      integer :: Ioff, Joff, IJoff
      integer :: Ilen, Itile, Jlen, Jtile, IJlen
      integer :: Cgrid, MyType, ghost, status, wtype
 
      integer, dimension(3) :: start, total
!
      real(r8) :: Afactor, Aoffset, Aspval
 
      real(r8), parameter :: IniVal = 0.0_r8
 
      real(r8), dimension(2) :: rbuffer
      real(r8), dimension(3) :: AttValue
 
      real(r8), allocatable :: Awrk(:,:)
#  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_fread2d"
!
!-----------------------------------------------------------------------
!  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=xtr_iobounds(ng)%xi_psi
          jsize=xtr_iobounds(ng)%eta_psi
        CASE (r2dvar, r3dvar)
          cgrid=2
          isize=xtr_iobounds(ng)%xi_rho
          jsize=xtr_iobounds(ng)%eta_rho
        CASE (u2dvar, u3dvar)
          cgrid=3
          isize=xtr_iobounds(ng)%xi_u
          jsize=xtr_iobounds(ng)%eta_u
        CASE (v2dvar, v3dvar)
          cgrid=4
          isize=xtr_iobounds(ng)%xi_v
          jsize=xtr_iobounds(ng)%eta_v
        CASE DEFAULT
          cgrid=2
          isize=xtr_iobounds(ng)%xi_rho
          jsize=xtr_iobounds(ng)%eta_rho
      END SELECT
 
      imin=xtr_bounds(ng)%Imin(cgrid,ghost,myrank)
      imax=xtr_bounds(ng)%Imax(cgrid,ghost,myrank)
      jmin=xtr_bounds(ng)%Jmin(cgrid,ghost,myrank)
      jmax=xtr_bounds(ng)%Jmax(cgrid,ghost,myrank)
 
      ilen=imax-imin+1
      jlen=jmax-jmin+1
!
!  Determine if interpolating from coarse gridded data to model grid
!  is required.  This is only allowed for gridded 2D fields.  This is
!  convenient for atmospheric forcing datasets that are usually on
!  coarser grids. The user can provide coarser gridded data to avoid
!  very large input files.
!
      interpolate=.false.
      IF (((vsize(1).gt.0).and.(vsize(1).ne.isize)).or.                 &
     &    ((vsize(2).gt.0).and.(vsize(2).ne.jsize))) THEN
        interpolate=.true.
        ilen=vsize(1)
        jlen=vsize(2)
      END IF
      IF (PRESENT(lregrid)) THEN
        lregrid=interpolate
      END IF
!
!  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)
            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
 
        npts=ilen*jlen
!
!  Allocate scratch work arrays.
!
        IF (interpolate) THEN
          IF (.not.allocated(awrk)) THEN
            allocate ( awrk(ilen,jlen) )
            awrk=inival
          END IF
          IF (.not.allocated(wrk)) THEN
            allocate ( wrk(npts) )
            wrk=inival
          END IF
        ELSE
          IF (.not.allocated(wrk)) THEN
            allocate ( wrk(tilesize(ng)) )
            wrk=0.0_r8
          END IF
        END IF
!
!  Read in data: all parallel nodes read their own tile data.
!
        IF (interpolate) THEN
          CALL tile_bounds_2d (ng, myrank, ilen, jlen, itile, jtile,    &
     &                         istr, iend, jstr, jend)
          start(1)=istr
          total(1)=iend-istr+1
          start(2)=jstr
          total(2)=jend-jstr+1
        ELSE
          start(1)=imin+ioff
          total(1)=ilen
          start(2)=jmin+joff
          total(2)=jlen
        END IF
        mynpts=total(1)*total(2)
        start(3)=tindex
        total(3)=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,mynpts
            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. If not interpolating, the data is loaded into
!  model array.
!
          IF (interpolate) THEN
            CALL mp_collect (ng, model, npts, inival, wrk)
            ic=0
            DO j=jstr,jend
              DO i=istr,iend
                ic=ic+1
                awrk(i,j)=wrk(ic)
              END DO
            END DO
          ELSE
            ic=0
            DO j=jmin,jmax
              DO i=imin,imax
                ic=ic+1
                adat(i,j)=wrk(ic)
              END DO
            END DO
          END IF
        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. Interpolation is not
!                allowed.
!-----------------------------------------------------------------------
!
      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 (p2dvar)
            npts=xtr_iobounds(ng)%xy_psi
            wtype=p2dvar
            ioff=0
            joff=1
          CASE (r2dvar)
            npts=xtr_iobounds(ng)%xy_rho
            wtype=r2dvar
            ioff=1
            joff=0
          CASE (u2dvar)
            npts=xtr_iobounds(ng)%xy_u
            wtype=u2dvar
            ioff=0
            joff=0
          CASE (v2dvar)
            npts=xtr_iobounds(ng)%xy_v
            wtype=v2dvar
            ioff=1
            joff=1
          CASE DEFAULT
            npts=xtr_iobounds(ng)%xy_rho
            wtype=r2dvar
            ioff=1
            joff=0
        END SELECT
        ijsize=isize*jsize
!
!  Allocate scratch work arrays.
!
        IF (.not.allocated(a2d)) THEN
          allocate ( a2d(ijsize) )
          a2d=inival
        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                     ! 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
!
!  Unpack read data.  This is tricky in parallel I/O.  The cheapeast
!  thing to do is reconstruct a packed global array and then select
!  the appropriate values for the tile.
!
          DO np=1,npts
            ij=scalars(ng)%IJwater(np,wtype)
            a2d(ij)=wrk(np)
          END DO
 
          DO j=jmin,jmax
            jc=(j-joff)*isize
            DO i=imin,imax
              ij=i+ioff+jc
              adat(i,j)=a2d(ij)
            END DO
          END DO
        ELSE
          exit_flag=2
          ioerror=status
        END IF
      END IF
#  endif
!
!-----------------------------------------------------------------------
!  Parallel I/O: If interpolating from gridded data, read its associated
!                locations and interpolate.
!-----------------------------------------------------------------------
!
      IF (interpolate.and.(status.eq.nf90_noerr)) THEN
        SELECT CASE (abs(mytype))
          CASE (p2dvar, p3dvar)
            IF (spherical) THEN
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, awrk, amin, amax,           &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          imin, imax, jmin, jmax,                 &
#  ifdef MASKING
     &                          grid(ng) % pmask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % lonp,                        &
     &                          grid(ng) % latp,                        &
     &                          adat)
            ELSE
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, awrk, amin, amax,           &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          imin, imax, jmin, jmax,                 &
#  ifdef MASKING
     &                          grid(ng) % pmask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % xp,                          &
     &                          grid(ng) % yp,                          &
     &                          adat)
            END IF
          CASE (r2dvar, r3dvar)
            IF (spherical) THEN
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, awrk, amin, amax,           &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          imin, imax, jmin, jmax,                 &
#  ifdef MASKING
     &                          grid(ng) % rmask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % lonr,                        &
     &                          grid(ng) % latr,                        &
     &                          adat)
            ELSE
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, awrk, amin, amax,           &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          imin, imax, jmin, jmax,                 &
#  ifdef MASKING
     &                          grid(ng) % rmask,                       &
#  endif
    &                           grid(ng) % MyLon,                       &
     &                          grid(ng) % xr,                          &
     &                          grid(ng) % yr,                          &
     &                          adat)
            END IF
          CASE (u2dvar, u3dvar)
            IF (spherical) THEN
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, awrk, amin, amax,           &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          imin, imax, jmin, jmax,                 &
#  ifdef MASKING
     &                          grid(ng) % umask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % lonu,                        &
     &                          grid(ng) % latu,                        &
     &                          adat)
            ELSE
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, awrk, amin, amax,           &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          imin, imax, jmin, jmax,                 &
#  ifdef MASKING
     &                          grid(ng) % umask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % xu,                          &
     &                          grid(ng) % yu,                          &
     &                          adat)
            END IF
          CASE (v2dvar, v3dvar)
            IF (spherical) THEN
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, awrk, amin, amax,           &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          imin, imax, jmin, jmax,                 &
#  ifdef MASKING
     &                          grid(ng) % vmask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % lonv,                        &
     &                          grid(ng) % latv,                        &
     &                          adat)
            ELSE
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, awrk, amin, amax,           &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          imin, imax, jmin, jmax,                 &
#  ifdef MASKING
     &                          grid(ng) % vmask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % xv,                          &
     &                          grid(ng) % yv,                          &
     &                          adat)
            END IF
          CASE DEFAULT
            IF (spherical) THEN
              CALL regrid_nf90 (ng, model, ncname, ncid,                  &
     &                          ncvname, ncvarid, mytype, interpflag,     &
     &                          ilen, jlen, awrk, amin, amax,             &
     &                          lbi, ubi, lbj, ubj,                       &
     &                          imin, imax, jmin, jmax,                   &
#  ifdef MASKING
     &                          grid(ng) % rmask,                         &
#  endif
     &                          grid(ng) % MyLon,                         &
     &                          grid(ng) % lonr,                          &
     &                          grid(ng) % latr,                          &
     &                          adat)
            ELSE
              CALL regrid_nf90 (ng, model, ncname, ncid,                  &
     &                          ncvname, ncvarid, mytype, interpflag,     &
     &                          ilen, jlen, awrk, amin, amax,             &
     &                          lbi, ubi, lbj, ubj,                       &
     &                          imin, imax, jmin, jmax,                   &
#  ifdef MASKING
     &                          grid(ng) % rmask,                         &
#  endif
     &                          grid(ng) % MyLon,                         &
     &                          grid(ng) % xr,                            &
     &                          grid(ng) % yr,                            &
     &                          adat)
            END IF
        END SELECT
      END IF
!
!-----------------------------------------------------------------------
!  Deallocate scratch work arrays.
!-----------------------------------------------------------------------
!
      IF (interpolate) THEN
        IF (allocated(awrk)) THEN
          deallocate ( awrk )
        END IF
      END IF
 
#  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_fread2d_xtr
 
# else
 
!
!***********************************************************************
      FUNCTION nf90_fread2d_xtr (ng, model, ncname, ncid,               &
     &                           ncvname, ncvarid,                      &
     &                           tindex, gtype, Vsize,                  &
     &                           LBi, UBi, LBj, UBj,                    &
     &                           Ascl, Amin, Amax,                      &
#  ifdef MASKING
     &                           Amask,                                 &
#  endif
     &                           Adat,                                  &
     &                           checksum,                              &
     &                           Lregrid) RESULT (status)
!***********************************************************************
!
      USE mod_netcdf
 
#  ifdef DISTRIBUTE
!
      USE distribute_mod, ONLY : mp_bcastf, mp_bcasti, mp_scatter2d_xtr
#  endif
      USE regrid_mod,     ONLY : regrid_nf90
!
!  Imported variable declarations.
!
      logical, intent(out), optional :: Lregrid
!
      integer, intent(in) :: ng, model, ncid, ncvarid, tindex, gtype
      integer, intent(in) :: LBi, UBi, LBj, UBj
      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:)
#  else
#   ifdef MASKING
      real(r8), intent(in) :: Amask(LBi:UBi,LBj:UBj)
#   endif
      real(r8), intent(out) :: Adat(LBi:UBi,LBj:UBj)
#  endif
!
!  Local variable declarations.
!
      logical :: Lchecksum, interpolate
 
      logical, dimension(3) :: foundit
!
      integer :: i, j, ic, Npts, NWpts, status, wtype
      integer :: Is, Ie, Js, Je
      integer :: Imin, Imax, Jmin, Jmax
      integer :: Ilen, Jlen, IJlen
      integer :: Cgrid, MyType, ghost
#  ifdef DISTRIBUTE
      integer :: Nghost
#  endif
      integer, dimension(3) :: start, total
!
      real(r8) :: Afactor, Aoffset, Aspval
 
      real(r8), dimension(3) :: AttValue
 
      real(r8), allocatable :: Cwrk(:)           ! used for checksum
      real(r8), allocatable :: wrk(:)
!
      character (len=12), dimension(3) :: AttName
 
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", nf90_fread2d"
!
!-----------------------------------------------------------------------
!  Set starting and ending indices to process.
!-----------------------------------------------------------------------
!
      status=nf90_noerr
!
!  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)
          cgrid=1                                         ! PSI-points
          is=xtr_iobounds(ng)%ILB_psi
          ie=xtr_iobounds(ng)%IUB_psi
          js=xtr_iobounds(ng)%JLB_psi
          je=xtr_iobounds(ng)%JUB_psi
        CASE (r2dvar)
          cgrid=2                                         ! RHO-points
          is=xtr_iobounds(ng)%ILB_rho
          ie=xtr_iobounds(ng)%IUB_rho
          js=xtr_iobounds(ng)%JLB_rho
          je=xtr_iobounds(ng)%JUB_rho
        CASE (u2dvar)
          cgrid=3                                         ! U-points
          is=xtr_iobounds(ng)%ILB_u
          ie=xtr_iobounds(ng)%IUB_u
          js=xtr_iobounds(ng)%JLB_u
          je=xtr_iobounds(ng)%JUB_u
        CASE (v2dvar)
          cgrid=4                                         ! V-points
          is=xtr_iobounds(ng)%ILB_v
          ie=xtr_iobounds(ng)%IUB_v
          js=xtr_iobounds(ng)%JLB_v
          je=xtr_iobounds(ng)%JUB_v
        CASE DEFAULT
          cgrid=2                                         ! RHO-points
          is=xtr_iobounds(ng)%ILB_rho
          ie=xtr_iobounds(ng)%IUB_rho
          js=xtr_iobounds(ng)%JLB_rho
          je=xtr_iobounds(ng)%JUB_rho
      END SELECT
 
      ilen=ie-is+1
      jlen=je-js+1
!
!  Set the tile computational I- and J-bounds (no ghost points).
!
      ghost=0
      imin=xtr_bounds(ng)%Imin(cgrid,ghost,myrank)
      imax=xtr_bounds(ng)%Imax(cgrid,ghost,myrank)
      jmin=xtr_bounds(ng)%Jmin(cgrid,ghost,myrank)
      jmax=xtr_bounds(ng)%Jmax(cgrid,ghost,myrank)
 
#  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                       ! no ghost points exchange
#   else
      IF (model.eq.iadm) THEN
        nghost=0                     ! no ghost points exchange
      ELSE
        nghost=nghostpoints          ! do ghost points exchange
      END IF
#   endif
#  endif
!
!  Determine if interpolating from coarse gridded data to model grid
!  is required.  This is only allowed for gridded 2D fields.  This is
!  convenient for atmospheric forcing datasets that are usually on
!  coarser grids. The user can provide coarser gridded data to avoid
!  very large input files.
!
      interpolate=.false.
      IF (((vsize(1).gt.0).and.(vsize(1).ne.ilen)).or.                  &
     &    ((vsize(2).gt.0).and.(vsize(2).ne.jlen))) THEN
        interpolate=.true.
        ilen=vsize(1)
        jlen=vsize(2)
      END IF
      IF (PRESENT(lregrid)) THEN
        lregrid=interpolate
      END IF
      ijlen=ilen*jlen
!
!  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
 
#  if defined READ_WATER && defined MASKING
!
!  If processing water points only, set number of points and type
!  switch.
!
      SELECT CASE (abs(mytype))
        CASE (p2dvar)
          npts=xtr_iobounds(ng)%xy_psi
          wtype=p2dvar
        CASE (r2dvar)
          npts=xtr_iobounds(ng)%xy_rho
          wtype=r2dvar
        CASE (u2dvar)
          npts=xtr_iobounds(ng)%xy_u
          wtype=u2dvar
        CASE (v2dvar)
          npts=xtr_iobounds(ng)%xy_v
          wtype=v2dvar
        CASE DEFAULT
          npts=xtr_iobounds(ng)%xy_rho
          wtype=r2dvar
      END SELECT
      nwpts=(lm(ng)+2)*(mm(ng)+2)
#  endif
!
!  Set NetCDF dimension counters for processing requested field.
!
      IF (mytype.gt.0) THEN
        npts=ijlen
        start(1)=1
        total(1)=ilen
        start(2)=1
        total(2)=jlen
        start(3)=tindex
        total(3)=1
#  if defined READ_WATER && defined MASKING
      ELSE
        start(1)=1
        total(1)=npts
        start(2)=1
        total(2)=tindex
#  endif
      END IF
!
!  Allocate scratch work vector. The dimension of this vector is
!  unknown when interpolating input data to model grid. Notice
!  that the array length is increased by two because the minimum
!  and maximum values are appended in distributed-memory
!  communications.
!
      IF (.not.allocated(wrk)) THEN
        IF (interpolate) THEN
          allocate ( wrk(npts) )
        ELSE
          allocate ( wrk(npts+2) )
        END IF
        wrk=0.0_r8
      END IF
!
!  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.
!-----------------------------------------------------------------------
!
      status=nf90_noerr
      IF (inpthread) THEN
        status=nf90_get_var(ncid, ncvarid, wrk, start, total)
        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
        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: If not interpolating, unpack read field.
!-----------------------------------------------------------------------
!
      IF (.not.interpolate) THEN
#  ifdef DISTRIBUTE
!
!  Scatter read data over the distributed memory tiles.
!
        CALL mp_scatter2d_xtr (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)
#  else
!
!  Unpack data into the global array: serial, serial with partitions,
!  and shared-memory applications.
!
        IF (mytype.gt.0) THEN
          ic=0
          DO j=js,je
            DO i=is,ie
              ic=ic+1
              adat(i,j)=wrk(ic)
            END DO
          END DO
#   if defined MASKING || defined READ_WATER
        ELSE
          ic=0
          DO j=js,je
            DO i=is,ie
              IF (amask(i,j).gt.0.0_r8) THEN
                ic=ic+1
                adat(i,j)=wrk(ic)
              ELSE
                adat(i,j)=0.0_r8
              END IF
            END DO
          END DO
#   endif
        END IF
#  endif
#  ifdef DISTRIBUTE
      ELSE
        CALL mp_bcastf (ng, model, wrk)
#  endif
      END IF
!
!-----------------------------------------------------------------------
!  Serial I/O: If interpolating from gridded data, read its associated
!  locations and interpolate.
!-----------------------------------------------------------------------
!
      IF (interpolate) THEN
        SELECT CASE (abs(mytype))
          CASE (p2dvar, p3dvar)
            IF (spherical) THEN
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, wrk, amin, amax,            &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          is, ie, js, je,                         &
#  ifdef MASKING
     &                          amask,                                  &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % lonp,                        &
     &                          grid(ng) % latp,                        &
     &                          adat)
            ELSE
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, wrk, amin, amax,            &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          is, ie, js, je,                         &
#  ifdef MASKING
     &                          amask,                                  &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % xp,                          &
     &                          grid(ng) % yp,                          &
     &                          adat)
            END IF
          CASE (r2dvar, r3dvar)
            IF (spherical) THEN
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, wrk, amin, amax,            &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          is, ie, js, je,                         &
#  ifdef MASKING
     &                          grid(ng) % rmask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % lonr,                        &
     &                          grid(ng) % latr,                        &
     &                          adat)
            ELSE
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, wrk, amin, amax,            &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          is, ie, js, je,                         &
#  ifdef MASKING
     &                          grid(ng) % rmask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % xr,                          &
     &                          grid(ng) % yr,                          &
     &                          adat)
            END IF
          CASE (u2dvar, u3dvar)
            IF (spherical) THEN
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, wrk, amin, amax,            &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          is, ie, js, je,                         &
#  ifdef MASKING
     &                          grid(ng) % umask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % lonu,                        &
     &                          grid(ng) % latu,                        &
     &                          adat)
            ELSE
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, wrk, amin, amax,            &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          is, ie, js, je,                         &
#  ifdef MASKING
     &                          grid(ng) % umask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % xu,                          &
     &                          grid(ng) % yu,                          &
     &                          adat)
            END IF
          CASE (v2dvar, v3dvar)
            IF (spherical) THEN
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, wrk, amin, amax,            &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          is, ie, js, je,                         &
#  ifdef MASKING
     &                          grid(ng) % vmask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % lonv,                        &
     &                          grid(ng) % latv,                        &
     &                          adat)
            ELSE
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, wrk, amin, amax,            &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          is, ie, js, je,                         &
#  ifdef MASKING
     &                          grid(ng) % vmask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % xv,                          &
     &                          grid(ng) % yv,                          &
     &                          adat)
            END IF
          CASE DEFAULT
            IF (spherical) THEN
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, wrk, amin, amax,            &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          is, ie, js, je,                         &
#  ifdef MASKING
     &                          grid(ng) % rmask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % lonr,                        &
     &                          grid(ng) % latr,                        &
     &                          adat)
            ELSE
              CALL regrid_nf90 (ng, model, ncname, ncid,                &
     &                          ncvname, ncvarid, mytype, interpflag,   &
     &                          ilen, jlen, wrk, amin, amax,            &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          is, ie, js, je,                         &
#  ifdef MASKING
     &                          grid(ng) % rmask,                       &
#  endif
     &                          grid(ng) % MyLon,                       &
     &                          grid(ng) % xr,                          &
     &                          grid(ng) % yr,                          &
     &                          adat)
            END IF
        END SELECT
      END IF
!
!-----------------------------------------------------------------------
!  If requested, compute data checksum value.
!-----------------------------------------------------------------------
!
      IF (lchecksum) THEN
#  ifdef DISTRIBUTE
        npts=(imax-imin+1)*(jmax-jmin+1)
        IF (.not.allocated(cwrk)) allocate ( cwrk(npts) )
        cwrk = pack(adat(imin:imax, jmin:jmax), .true.)
        CALL get_hash (cwrk, npts, checksum, .true.)
#  else
        npts=(ie-is+1)*(je-js+1)
        IF (.not.allocated(cwrk)) allocate ( cwrk(npts) )
        cwrk = pack(adat(is:ie, js:je), .true.)
        CALL get_hash (cwrk, npts, checksum)
#  endif
        IF (allocated(cwrk)) deallocate (cwrk)
      END IF
!
!-----------------------------------------------------------------------
!  Deallocate scratch work vector.
!-----------------------------------------------------------------------
!
      IF (allocated(wrk)) THEN
        deallocate (wrk)
      END IF
!
      RETURN
      END FUNCTION nf90_fread2d_xtr
# endif
 
# if defined PIO_LIB && defined DISTRIBUTE
!
!***********************************************************************
      FUNCTION pio_fread2d_xtr (ng, model, ncname, pioFile,             &
     &                          ncvname, pioVar,                        &
     &                          tindex, pioDesc, Vsize,                 &
     &                          LBi, UBi, LBj, UBj,                     &
     &                          Ascl, Amin, Amax,                       &
#  ifdef MASKING
     &                          Amask,                                  &
#  endif
     &                          Adat,                                   &
     &                          checksum,                               &
     &                          Lregrid) RESULT (status)
!***********************************************************************
!
      USE mod_pio_netcdf
!
      USE distribute_mod, ONLY : mp_reduce
      USE regrid_mod,     ONLY : regrid_pio
!
!  Imported variable declarations.
!
      logical, intent(out), optional :: Lregrid
!
      integer, intent(in) :: ng, model, tindex
      integer, intent(in) :: LBi, UBi, LBj, UBj
      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:)
#  else
#   ifdef MASKING
      real(r8), intent(in) :: Amask(LBi:UBi,LBj:UBj)
#   endif
      real(r8), intent(out) :: Adat(LBi:UBi,LBj:UBj)
#  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, interpolate
 
      logical, dimension(3) :: foundit
!
      integer :: i, j, Npts, status
      integer :: Is, Ie, Js, Je                 ! global bounds
      integer :: Imin, Imax, Jmin, Jmax         ! tile bounds
      integer :: Ilen, Jlen, IJlen
      integer :: Cgrid, ghost, dkind, gtype
 
      integer, dimension(3) :: 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
      real(r8), allocatable :: wrk(:,:)          ! interpolating, regrid
!
      character (len=12), dimension(3) :: AttName
      character (len= 3), dimension(2) :: op_handle
 
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", pio_fread2d"
!
!-----------------------------------------------------------------------
!  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=xtr_iobounds(ng)%ILB_psi
          ie=xtr_iobounds(ng)%IUB_psi
          js=xtr_iobounds(ng)%JLB_psi
          je=xtr_iobounds(ng)%JUB_psi
        CASE (r2dvar, r3dvar)
          cgrid=2                                         ! RHO-points
          is=xtr_iobounds(ng)%ILB_rho
          ie=xtr_iobounds(ng)%IUB_rho
          js=xtr_iobounds(ng)%JLB_rho
          je=xtr_iobounds(ng)%JUB_rho
        CASE (u2dvar, u3dvar)
          cgrid=3                                         ! U-points
          is=xtr_iobounds(ng)%ILB_u
          ie=xtr_iobounds(ng)%IUB_u
          js=xtr_iobounds(ng)%JLB_u
          je=xtr_iobounds(ng)%JUB_u
        CASE (v2dvar, v3dvar)
          cgrid=4                                         ! V-points
          is=xtr_iobounds(ng)%ILB_v
          ie=xtr_iobounds(ng)%IUB_v
          js=xtr_iobounds(ng)%JLB_v
          je=xtr_iobounds(ng)%JUB_v
        CASE DEFAULT
          cgrid=2                                         ! RHO-points
          is=xtr_iobounds(ng)%ILB_rho
          ie=xtr_iobounds(ng)%IUB_rho
          js=xtr_iobounds(ng)%JLB_rho
          je=xtr_iobounds(ng)%JUB_rho
      END SELECT
!
!  Compute the global lengths of the I- and J-directions. They are
!  needed to check if regridding is required.
!
      ilen=ie-is+1
      jlen=je-js+1
!
!  Set the tile computational I- and J-bounds (no ghost points).
!
      ghost=0
      imin=xtr_bounds(ng)%Imin(cgrid,ghost,myrank)
      imax=xtr_bounds(ng)%Imax(cgrid,ghost,myrank)
      jmin=xtr_bounds(ng)%Jmin(cgrid,ghost,myrank)
      jmax=xtr_bounds(ng)%Jmax(cgrid,ghost,myrank)
!
!  Determine if interpolating from coarse gridded data to model grid
!  is required.  This is only allowed for gridded 2D fields.  This is
!  convenient for atmospheric forcing datasets that are usually on
!  coarser grids. The user can provide coarser gridded data to avoid
!  very large input files.
!
      interpolate=.false.
      IF (((vsize(1).gt.0).and.(vsize(1).ne.ilen)).or.                  &
     &    ((vsize(2).gt.0).and.(vsize(2).ne.jlen))) THEN
        interpolate=.true.
        ilen=vsize(1)          ! data and state variable are incongruent
        jlen=vsize(2)          ! horizontal interpolation is required
      END IF
      IF (PRESENT(lregrid)) THEN
        lregrid=interpolate
      END IF
      ijlen=ilen*jlen
!
!  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
!
!-----------------------------------------------------------------------
!  If interpolting, read in input data with the non-tiled PIO interface.
!-----------------------------------------------------------------------
!
      IF (interpolate) THEN
        IF (.not.allocated(wrk)) THEN
          allocate ( wrk(ilen,jlen) )
        END IF
        wrk=0.0_r8
!
        start(1)=1
        total(1)=ilen
        start(2)=1
        total(2)=jlen
        start(3)=tindex
        total(3)=1
!
!  The generic routine "pio_netcdf_get_fvar" checks if the attributes
!  "scale_factor", "add_offset", and "_FillValue" are present in the
!  input NetCDF variable and it applies its operators.
!
        CALL pio_netcdf_get_fvar (ng, model, ncname,                    &
     &                            ncvname, wrk,                         &
     &                            piofile = piofile,                    &
     &                            start = start,                        &
     &                            total = total,                        &
     &                            broadcast = .false.,                  &
     &                            min_val = amin,                       &
     &                            max_val = amax)
!
!  Multipy by user metadata scale.
!
        DO j=1,jlen
          DO i=1,ilen
            wrk(i,j)=ascl*wrk(i,j)
          END DO
        END DO
      END IF
!
!  Initialize checsum value.
!
      IF (PRESENT(checksum)) THEN
        lchecksum=.true.
        checksum=0_i8b
      ELSE
        lchecksum=.false.
      END IF
!
!-----------------------------------------------------------------------
!  If not interpolated, read in requested tiled field and scale it.
!-----------------------------------------------------------------------
!
      IF (.not.interpolate) THEN
!
!  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) )
          END IF
          awrk8=0.0_r8
        ELSE                                    ! single precision
          IF (.not.associated(awrk4)) THEN
            allocate ( awrk4(imin:imax, jmin:jmax) )
          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 and load double precision data from NetCDF file.
!
        IF (dkind.eq.pio_double) THEN
          CALL pio_read_darray (piofile,                                &
     &                          piovar%vd,                              &
     &                          piodesc,                                &
     &                          awrk8(imin:,jmin:),                     &
     &                          status)
!
          DO j=jmin,jmax
            DO i=imin,imax
              IF (abs(awrk8(i,j)).ge.abs(aspval)) THEN
                adat(i,j)=0.0_r8              ! masked with _FillValue
              ELSE
                avalue=ascl*(afactor*awrk8(i,j)+aoffset)
                adat(i,j)=avalue
                my_amin=min(my_amin,avalue)
                my_amax=max(my_amax,avalue)
              END IF
            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:),                     &
     &                          status)
!
          DO j=jmin,jmax
            DO i=imin,imax
              IF (abs(awrk4(i,j)).ge.abs(real(aspval,r4))) THEN
                adat(i,j)=0.0_r8               ! masked with _FillValue
              ELSE
                avalue=real(ascl*(afactor*awrk4(i,j)+aoffset),r8)
                adat(i,j)=avalue
                my_amin=real(min(my_amin,avalue),r8)
                my_amax=real(max(my_amax,avalue),r8)
              END IF
            END DO
          END DO
          IF (associated(awrk4)) deallocate (awrk4)
        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
      END IF
!
!-----------------------------------------------------------------------
!  If interpolating from gridded data, read its associated locations
!  and interpolate.
!-----------------------------------------------------------------------
!
      IF (interpolate) THEN
        SELECT CASE (gtype)
          CASE (p2dvar, p3dvar)
            IF (spherical) THEN
              CALL regrid_pio (ng, model, ncname, piofile,              &
     &                         ncvname, piovar, gtype, interpflag,      &
     &                         ilen, jlen, wrk, amin, amax,             &
     &                         lbi, ubi, lbj, ubj,                      &
     &                         imin, imax, jmin, jmax,                  &
#  ifdef MASKING
     &                         amask,                                   &
#  endif
     &                         grid(ng) % MyLon,                        &
     &                         grid(ng) % lonp,                         &
     &                         grid(ng) % latp,                         &
     &                         adat)
            ELSE
              CALL regrid_pio (ng, model, ncname, piofile,              &
     &                         ncvname, piovar, gtype, interpflag,      &
     &                         ilen, jlen, wrk, amin, amax,             &
     &                         lbi, ubi, lbj, ubj,                      &
     &                         imin, imax, jmin, jmax,                  &
#  ifdef MASKING
     &                         amask,                                   &
#  endif
     &                         grid(ng) % MyLon,                        &
     &                         grid(ng) % xp,                           &
     &                         grid(ng) % yp,                           &
     &                         adat)
            END IF
          CASE (r2dvar, r3dvar)
            IF (spherical) THEN
              CALL regrid_pio (ng, model, ncname, piofile,              &
     &                         ncvname, piovar, gtype, interpflag,      &
     &                         ilen, jlen, wrk, amin, amax,             &
     &                         lbi, ubi, lbj, ubj,                      &
     &                         imin, imax, jmin, jmax,                  &
#  ifdef MASKING
     &                         grid(ng) % rmask,                        &
#  endif
     &                         grid(ng) % MyLon,                        &
     &                         grid(ng) % lonr,                         &
     &                         grid(ng) % latr,                         &
     &                         adat)
            ELSE
              CALL regrid_pio (ng, model, ncname, piofile,              &
     &                         ncvname, piovar, gtype, interpflag,      &
     &                         ilen, jlen, wrk, amin, amax,             &
     &                         lbi, ubi, lbj, ubj,                      &
     &                         imin, imax, jmin, jmax,                  &
#  ifdef MASKING
     &                         grid(ng) % rmask,                        &
#  endif
     &                         grid(ng) % MyLon,                        &
     &                         grid(ng) % xr,                           &
     &                         grid(ng) % yr,                           &
     &                         adat)
            END IF
          CASE (u2dvar, u3dvar)
            IF (spherical) THEN
              CALL regrid_pio (ng, model, ncname, piofile,              &
     &                         ncvname, piovar, gtype, interpflag,      &
     &                         ilen, jlen, wrk, amin, amax,             &
     &                         lbi, ubi, lbj, ubj,                      &
     &                         imin, imax, jmin, jmax,                  &
#  ifdef MASKING
     &                         grid(ng) % umask,                        &
#  endif
     &                         grid(ng) % MyLon,                        &
     &                         grid(ng) % lonu,                         &
     &                         grid(ng) % latu,                         &
     &                         adat)
            ELSE
              CALL regrid_pio (ng, model, ncname, piofile,              &
     &                         ncvname, piovar, gtype, interpflag,      &
     &                         ilen, jlen, wrk, amin, amax,             &
     &                         lbi, ubi, lbj, ubj,                      &
     &                         imin, imax, jmin, jmax,                  &
#  ifdef MASKING
     &                         grid(ng) % umask,                        &
#  endif
     &                         grid(ng) % MyLon,                        &
     &                         grid(ng) % xu,                           &
     &                         grid(ng) % yu,                           &
     &                         adat)
            END IF
          CASE (v2dvar, v3dvar)
            IF (spherical) THEN
              CALL regrid_pio (ng, model, ncname, piofile,              &
     &                         ncvname, piovar, gtype, interpflag,      &
     &                         ilen, jlen, wrk, amin, amax,             &
     &                         lbi, ubi, lbj, ubj,                      &
     &                         imin, imax, jmin, jmax,                  &
#  ifdef MASKING
     &                         grid(ng) % vmask,                        &
#  endif
     &                         grid(ng) % MyLon,                        &
     &                         grid(ng) % lonv,                         &
     &                         grid(ng) % latv,                         &
     &                         adat)
            ELSE
              CALL regrid_pio (ng, model, ncname, piofile,              &
     &                         ncvname, piovar, gtype, interpflag,      &
     &                         ilen, jlen, wrk, amin, amax,             &
     &                         lbi, ubi, lbj, ubj,                      &
     &                         imin, imax, jmin, jmax,                  &
#  ifdef MASKING
     &                         grid(ng) % vmask,                        &
#  endif
     &                         grid(ng) % MyLon,                        &
     &                         grid(ng) % xv,                           &
     &                         grid(ng) % yv,                           &
     &                         adat)
            END IF
          CASE DEFAULT
            IF (spherical) THEN
              CALL regrid_pio (ng, model, ncname, piofile,              &
     &                         ncvname, piovar, gtype, interpflag,      &
     &                         ilen, jlen, wrk, amin, amax,             &
     &                         lbi, ubi, lbj, ubj,                      &
     &                         imin, imax, jmin, jmax,                  &
#  ifdef MASKING
     &                         grid(ng) % rmask,                        &
#  endif
     &                         grid(ng) % MyLon,                        &
     &                         grid(ng) % lonr,                         &
     &                         grid(ng) % latr,                         &
     &                         adat)
            ELSE
              CALL regrid_pio (ng, model, ncname, piofile,              &
     &                         ncvname, piovar, gtype, interpflag,      &
     &                         ilen, jlen, wrk, amin, amax,             &
     &                         lbi, ubi, lbj, ubj,                      &
     &                         imin, imax, jmin, jmax,                  &
#  ifdef MASKING
     &                         grid(ng) % rmask,                        &
#  endif
     &                         grid(ng) % MyLon,                        &
     &                         grid(ng) % xr,                           &
     &                         grid(ng) % yr,                           &
     &                         adat)
            END IF
        END SELECT
!
!  Deallocate regridding work array.
!
        IF (allocated(wrk)) THEN
          deallocate (wrk)
        END IF
      END IF
!
!  If requested, compute input data checksum value.
!
      IF (lchecksum) THEN
        npts=(imax-imin+1)*(jmax-jmin+1)
        IF (.not.allocated(cwrk)) allocate ( cwrk(npts) )
        cwrk=pack(adat(imin:imax, jmin:jmax), .true.)
        CALL get_hash (cwrk, npts, checksum, .true.)
        IF (allocated(cwrk)) deallocate (cwrk)
      END IF
!
      RETURN
      END FUNCTION pio_fread2d_xtr
# endif
#endif
      END MODULE nf_fread2d_xtr_mod