pack_field.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE pack_field_mod
!
!git $Id$
!=======================================================================
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                            Hernan G. Arango   !
!=======================================================================
!                                                                      !
!  These routines pack output variables data into a 1D array. If       !
!  distributed memory, it gathers tiled data from all processors       !
!  and serializes.                                                     !
!                                                                      !
!  Then, if applicable, they extract fields at different geometry by   !
!  decimation or interpolation.                                        !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_scalars
!
#ifdef DISTRIBUTE
      USE distribute_mod,   ONLY : mp_collect,                          &
     &                             mp_gather2d,                         &
# ifdef grid_extract
     &                             mp_gather2d_xtr,                     &
# endif
     &                             mp_gather3d
#endif
#ifdef GRID_EXTRACT
# ifdef ADJUST_BOUNDARY
      USE extract_field_mod, ONLY: extract_boundary
# endif
      USE extract_field_mod, ONLY: extract_field
#endif
!
      implicit none
!
#ifdef ADJUST_BOUNDARY
      PUBLIC :: pack_boundary2d
      PUBLIC :: pack_boundary3d
#endif
      PUBLIC :: pack_field2d
      PUBLIC :: pack_field3d
      PUBLIC :: pack_field4d
!
      CONTAINS
 
#ifdef ADJUST_BOUNDARY
!
      SUBROUTINE pack_boundary2d (ng, model, tile,                      &
     &                            gtype, ncvname, tindex,               &
     &                            Extract_Flag,                         &
     &                            LBij, UBij, Nrec,                     &
     &                            Bscl, Bdat,                           &
     &                            start, total, Npts, Bwrk)
!
!=======================================================================
!                                                                      !
!  Packs 2D boundary data into 1D array to be written into output      !
!  NetCDF file elsewhere.                                              !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng           Nested grid number (integer)                        !
!     model        Calling model identifier (integer)                  !
!     tile         Domain partition (integer)                          !
!     gtype        Staggered C-grid type (integer)                     !
!     ncvname      NetCDF variable name (string)                       !
!     tindex       Time record index to process (integer)              !
!     Extract_Flag Extraction flag interpolation/decimation (integer)  !
!     LBij         Donor IJ-dimension Lower bound (integer)            !
!     UBij         Donor IJ-dimension Upper bound (integer)            !
!     Nrec         Number of boundary records (integer)                !
!     Bscl         Factor to scale boundary data before writing (real) !
!     Bdat         2D boundary data to process (real)                  !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     start        Start index where the first of the data values will !
!                    be written along each dimension (integer)         !
!     total        Number of data values to be written along each      !
!                    dimension (integer)                               !
!     Npts         Number of points processed in Bwrk (integer)        !
!     Bwrk         Packed 2D boundary data (real)                      !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      integer, intent(in)   :: ng, model, tile
      integer, intent(in)   :: gtype, tindex
      integer, intent(in)   :: extract_flag
      integer, intent(in)   :: lbij, ubij, nrec
      integer, intent(in)   :: npts
      integer, intent(out)  :: start(:), total(:)
!
      real(dp), intent(in)  :: bscl
!
# ifdef ASSUMED_SHAPE
      real(r8), intent(in)  :: bdat(lbij:,:,:)
# else
      real(r8), intent(in)  :: bdat(lbij:lbij,4,nrec)
# endif
      real(r8), intent(out) :: bwrk(:)
!
      character (len=*), intent(in) :: ncvname
!
!  Local variable declarations.
!
      logical, dimension(4) :: bounded
!
      integer :: bc, i, ib, ic, ir, j, rc
      integer :: iorj, imin, imax, jmin, jmax
!
      real(r8), parameter :: bspv = 0.0_r8
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", pack_boundary2d"
!
!------------------------------------------------------------------------
!  Pack 2D boundary data.
!------------------------------------------------------------------------
!
!  Set first and last grid point according to staggered C-grid
!  classification.
!
# ifdef DISTRIBUTE
      SELECT CASE (gtype)
!
!       ---------------------
        CASE (p2dvar, p3dvar)               ! PSI points field
!       ---------------------
!
          imin=bounds(ng)%Istr (tile)
          imax=bounds(ng)%Iend (tile)
          jmin=bounds(ng)%Jstr (tile)
          jmax=bounds(ng)%Jend (tile)
!
!       ---------------------
        CASE (r2dvar, r3dvar)               ! RHO points field
!       ---------------------
!
          imin=bounds(ng)%IstrR(tile)
          imax=bounds(ng)%IendR(tile)
          jmin=bounds(ng)%JstrR(tile)
          jmax=bounds(ng)%JendR(tile)
!
!       ---------------------
        CASE (u2dvar, u3dvar)               ! U points field
!       ---------------------
!
          imin=bounds(ng)%Istr (tile)
          imax=bounds(ng)%IendR(tile)
          jmin=bounds(ng)%JstrR(tile)
          jmax=bounds(ng)%JendR(tile)
!
!       ---------------------
        CASE (v2dvar, v3dvar)               ! V points field
!       ---------------------
!
          imin=bounds(ng)%IstrR(tile)
          imax=bounds(ng)%IendR(tile)
          jmin=bounds(ng)%Jstr (tile)
          jmax=bounds(ng)%JendR(tile)
!
!       ---------------------
        CASE DEFAULT                        ! RHO points field
!       ---------------------
!
          imin=bounds(ng)%IstrR(tile)
          imax=bounds(ng)%IendR(tile)
          jmin=bounds(ng)%JstrR(tile)
          jmax=bounds(ng)%JendR(tile)
      END SELECT
# else
      imin=lbij
      imax=ubij
      jmin=lbij
      jmax=ubij
# endif
!
      iorj=iobounds(ng)%IorJ
!
!  Set switch to process boundary data by their associated tiles.
!
      bounded(iwest )=domain(ng)%Western_Edge (tile)
      bounded(ieast )=domain(ng)%Eastern_Edge (tile)
      bounded(isouth)=domain(ng)%Southern_Edge(tile)
      bounded(inorth)=domain(ng)%Northern_Edge(tile)
!
!-----------------------------------------------------------------------
!  Pack and scale output data.
!-----------------------------------------------------------------------
!
      bwrk=bspv
!
      DO ir=1,nrec
        rc=(ir-1)*iorj*4
        DO ib=1,4
          IF (bounded(ib)) THEN
            bc=(ib-1)*iorj+rc
            IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN
              DO j=jmin,jmax
                ic=1+(j-lbij)+bc
                bwrk(ic)=bdat(j,ib,ir)*bscl
# ifdef POSITIVE_ZERO
                IF (abs(bwrk(ic)).eq.0.0_r8) THEN
                  bwrk(ic)=0.0_r8                 ! impose positive zero
                END IF
# endif
              END DO
            ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN
              DO i=imin,imax
                ic=1+(i-lbij)+bc
                bwrk(ic)=bdat(i,ib,ir)*bscl
# ifdef POSITIVE_ZERO
                IF (abs(bwrk(ic)).eq.0.0_r8) THEN
                  bwrk(ic)=0.0_r8                 ! impose positive zero
                END IF
# endif
              END DO
            END IF
          END IF
        END DO
      END DO
 
# ifdef DISTRIBUTE
!
!  If distributed-memory set-up, collect data from all spawned
!  processes.
!
      CALL mp_collect (ng, model, npts, bspv, bwrk)
# endif
!
!-----------------------------------------------------------------------
!  If there is no extracting data, set the start and total vectors
!  needed by the NetCDF library for writing.
!-----------------------------------------------------------------------
!
      IF (extract_flag.le.0) THEN
        start(1)=1
        total(1)=iorj
        start(2)=1
        total(2)=4
        start(3)=1
        total(3)=nrec
        start(4)=tindex
        total(4)=1
 
# ifdef GRID_EXTRACT
!
!-----------------------------------------------------------------------
!  Otherwise, extract field by decimation or interpolation.
!-----------------------------------------------------------------------
!
      ELSE IF (extract_flag.ge.1) THEN
        CALL extract_boundary (ng, model, tile,                         &
     &                         gtype, ncvname, tindex,                  &
     &                         extract_flag,                            &
     &                         imin, imax, jmin, jmax, nrec,            &
     &                         start, total,                            &
     &                         npts, bwrk)
# endif
      END IF
!
      RETURN
      END SUBROUTINE pack_boundary2d
!
      SUBROUTINE pack_boundary3d (ng, model, tile,                      &
     &                            gtype, ncvname, tindex,               &
     &                            Extract_Flag,                         &
     &                            LBij, UBij, LBk, UBk, Nrec,           &
     &                            Bscl, Bdat,                           &
     &                            start, total, Npts, Bwrk)
!
!=======================================================================
!                                                                      !
!  Packs 3D boundary data into 1D array to be written into output      !
!  NetCDF file elsewhere.                                              !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng           Nested grid number (integer)                        !
!     model        Calling model identifier (integer)                  !
!     tile         Domain partition (integer)                          !
!     gtype        Staggered C-grid type (integer)                     !
!     ncvname      NetCDF variable name (string)                       !
!     tindex       Time record index to process (integer)              !
!     Extract_Flag Extraction flag interpolation/decimation (integer)  !
!     LBij         Donor IJ-dimension Lower bound (integer)            !
!     UBij         Donor IJ-dimension Upper bound (integer)            !
!     LBk          Donor K-dimension lower bound (integer)             !
!     UBk          Donor K-dimension upper bound (integer)             !
!     Nrec         Number of boundary records (integer)                !
!     Bscl         Factor to scale boundary data before writing (real) !
!     Bdat         3D boundary data to process (real)                  !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     start        Start index where the first of the data values will !
!                    be written along each dimension (integer)         !
!     total        Number of data values to be written along each      !
!                    dimension (integer)                               !
!     Npts         Number of points processed in Bwrk (integer)        !
!     Bwrk         Packed 3D boundary data (real)                      !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      integer, intent(in)   :: ng, model, tile
      integer, intent(in)   :: gtype, tindex
      integer, intent(in)   :: extract_flag
      integer, intent(in)   :: lbij, ubij, lbk, ubk, nrec
      integer, intent(in)   :: npts
      integer, intent(out)  :: start(:), total(:)
!
      real(dp), intent(in)  :: bscl
!
# ifdef ASSUMED_SHAPE
      real(r8), intent(in)  :: bdat(lbij:,lbk:,:,:)
# else
      real(r8), intent(in)  :: bdat(lbij:lbij,lbk:ubk,4,nrec)
# endif
      real(r8), intent(out) :: bwrk(:)
!
      character (len=*), intent(in) :: ncvname
!
!  Local variable declarations.
!
      logical, dimension(4) :: bounded
!
      integer :: imin, imax, jmin, jmax
      integer :: iorj, ijklen, klen
      integer :: bc, i, ib, ic, ir, j, k, kc, rc
!
      real(r8), parameter :: bspv = 0.0_r8
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", pack_boundary3d"
!
!------------------------------------------------------------------------
!  Pack 3D boundary data.
!------------------------------------------------------------------------
!
!  Set first and last grid point according to staggered C-grid
!  classification.
!
# ifdef DISTRIBUTE
      SELECT CASE (gtype)
!
!       ---------------------
        CASE (p2dvar, p3dvar)               ! PSI points field
!       ---------------------
!
          imin=bounds(ng)%Istr (tile)
          imax=bounds(ng)%Iend (tile)
          jmin=bounds(ng)%Jstr (tile)
          jmax=bounds(ng)%Jend (tile)
!
!       ---------------------
        CASE (r2dvar, r3dvar)               ! RHO points field
!       ---------------------
!
          imin=bounds(ng)%IstrR(tile)
          imax=bounds(ng)%IendR(tile)
          jmin=bounds(ng)%JstrR(tile)
          jmax=bounds(ng)%JendR(tile)
!
!       ---------------------
        CASE (u2dvar, u3dvar)               ! U points field
!       ---------------------
!
          imin=bounds(ng)%Istr (tile)
          imax=bounds(ng)%IendR(tile)
          jmin=bounds(ng)%JstrR(tile)
          jmax=bounds(ng)%JendR(tile)
!
!       ---------------------
        CASE (v2dvar, v3dvar)               ! V points field
!       ---------------------
!
          imin=bounds(ng)%IstrR(tile)
          imax=bounds(ng)%IendR(tile)
          jmin=bounds(ng)%Jstr (tile)
          jmax=bounds(ng)%JendR(tile)
!
!       ---------------------
        CASE DEFAULT                        ! RHO points field
!       ---------------------
!
          imin=bounds(ng)%IstrR(tile)
          imax=bounds(ng)%IendR(tile)
          jmin=bounds(ng)%JstrR(tile)
          jmax=bounds(ng)%JendR(tile)
      END SELECT
# else
      imin=lbij
      imax=ubij
      jmin=lbij
      jmax=ubij
# endif
!
      iorj=iobounds(ng)%IorJ
      klen=ubk-lbk+1
      ijklen=iorj*klen
!
!  Set switch to process boundary data by their associated tiles.
!
      bounded(iwest )=domain(ng)%Western_Edge (tile)
      bounded(ieast )=domain(ng)%Eastern_Edge (tile)
      bounded(isouth)=domain(ng)%Southern_Edge(tile)
      bounded(inorth)=domain(ng)%Northern_Edge(tile)
!
!-----------------------------------------------------------------------
!  Pack and scale output data.
!-----------------------------------------------------------------------
!
      bwrk=bspv
!
      DO ir=1,nrec
        rc=(ir-1)*ijklen*4
        DO ib=1,4
          IF (bounded(ib)) THEN
            bc=(ib-1)*ijklen+rc
            IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN
              DO k=lbk,ubk
                kc=(k-lbk)*iorj+bc
                DO j=jmin,jmax
                  ic=1+(j-lbij)+kc
                  bwrk(ic)=bdat(j,k,ib,ir)*bscl
# ifdef POSITIVE_ZERO
                  IF (abs(bwrk(ic)).eq.0.0_r8) THEN
                    bwrk(ic)=0.0_r8               ! impose positive zero
                  END IF
# endif
                END DO
              END DO
            ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN
              DO k=lbk,ubk
                kc=(k-lbk)*iorj+bc
                DO i=imin,imax
                  ic=1+(i-lbij)+kc
                  bwrk(ic)=bdat(i,k,ib,ir)*bscl
# ifdef POSITIVE_ZERO
                  IF (abs(bwrk(ic)).eq.0.0_r8) THEN
                    bwrk(ic)=0.0_r8               ! impose positive zero
                  END IF
# endif
                END DO
              END DO
            END IF
          END IF
        END DO
      END DO
 
# ifdef DISTRIBUTE
!
!  If distributed-memory set-up, collect data from all spawned
!  processes.
!
      CALL mp_collect (ng, model, npts, bspv, bwrk)
# endif
!
!-----------------------------------------------------------------------
!  If there is no extracting data, set the start and total vectors
!  needed by the NetCDF library for writing.
!-----------------------------------------------------------------------
!
      IF (extract_flag.le.0) THEN
        start(1)=1
        total(1)=iorj
        start(2)=1
        total(2)=klen
        start(3)=1
        total(3)=4
        start(4)=1
        total(4)=nrec
        start(5)=tindex
        total(5)=1
 
# ifdef GRID_EXTRACT
!
!-----------------------------------------------------------------------
!  Otherwise, extract field by decimation or interpolation.
!-----------------------------------------------------------------------
!
      ELSE IF (extract_flag.ge.1) THEN
        CALL extract_boundary (ng, model, tile,                         &
     &                         gtype, ncvname, tindex,                  &
     &                         extract_flag,                            &
     &                         imin, imax, jmin, jmax, lbk, ubk, nrec,  &
     &                         start, total,                            &
     &                         npts, bwrk)
# endif
      END IF
!
      RETURN
      END SUBROUTINE pack_boundary3d
#endif
!
      SUBROUTINE pack_field2d (ng, model, tile,                         &
     &                         gtype, ifield, tindex,                   &
     &                         LandFill, Extract_Flag,                  &
     &                         LBi, UBi, LBj, UBj,                      &
#ifdef MASKING
     &                         Amask,                                   &
#endif
     &                         Ascl, Adat,                              &
     &                         start, total, Npts, Awrk)
!
!=======================================================================
!                                                                      !
!  Packs 2D field data into 1D array to be written into output NetCDF  !
!  file elsewhere.                                                     !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng           Nested grid number (integer)                        !
!     model        Calling model identifier (integer)                  !
!     tile         Domain partition (integer)                          !
!     gtype        Staggered C-grid type (integer)                     !
!     ifield       Field metadata index (integer)                      !
!     tindex       Time record index to process (integer)              !
!     LandFill     Switch to set fill value in land areas (logical)    !                                         !
!     Extract_Flag Extraction flag interpolation/decimation (integer)  !
!     LBi          Donor field I-dimension Lower bound (integer)       !
!     UBi          Donor field I-dimension Upper bound (integer)       !
!     LBj          Donor field J-dimension Lower bound (integer)       !
!     UBj          Donor field J-dimension Upper bound (integer)       !
#ifdef MASKING
!     Amask        Land/Sea mask, if any (real 2D array)               !
#endif
!     Ascl         Factor to scale field after reading (real)          !
!     Adat         2D field to process (real)                          !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     start        Start index where the first of the data values will !
!                    be written along each dimension (integer)         !
!     total        Number of data values to be written along each      !
!                    dimension (integer)                               !
!     Npts         Number of points processed in Awrk.                 !
!     Awrk         Packed 2D field data (real)                         !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      logical, intent(in)  :: landfill
!
      integer, intent(in)  :: ng, model, tile
      integer, intent(in)  :: gtype, ifield, tindex
      integer, intent(in)  :: extract_flag
      integer, intent(in)  :: lbi, ubi, lbj, ubj
      integer, intent(out) :: npts
      integer, intent(out) :: start(:), total(:)
!
      real(dp), intent(in) :: ascl
!
#ifdef ASSUMED_SHAPE
# ifdef MASKING
      real(r8), intent(in) :: amask(lbi:,lbj:)
# endif
      real(r8), intent(in) :: adat(lbi:,lbj:)
#else
# ifdef MASKING
      real(r8), intent(in) :: amask(lbi:ubi,lbj:ubj)
# endif
      real(r8), intent(in) :: adat(lbi:ubi,lbj:ubj)
#endif
      real(r8), intent(out) :: awrk(:)
!
!  Local variable declarations.
!
      integer :: i, ic, j
      integer :: imin, imax, jmin, jmax
      integer :: isize, jsize, ijsize, mytype
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", pack_field2d"
 
!-----------------------------------------------------------------------
!  Set starting and ending indices to process.
!-----------------------------------------------------------------------
!
!  Set first and last grid point according to staggered C-grid
!  classification. Set loops offsets.
#ifdef GRID_EXTRACT
!
!  Here, the ELSE statements (Extract_Flag < 0) indicate the processing
!  of data already on the extracted grid, like geometry variables
!  requiring no extraction since they are read during initialization.
#endif
!
      mytype=gtype
!
      SELECT CASE (abs(mytype))
!
!       ---------------------
        CASE (p2dvar, p3dvar)               ! PSI points field
!       ---------------------
!
          IF (extract_flag.ge.0) THEN
            imin=iobounds(ng)%ILB_psi
            imax=iobounds(ng)%IUB_psi
            jmin=iobounds(ng)%JLB_psi
            jmax=iobounds(ng)%JUB_psi
#ifdef GRID_EXTRACT
          ELSE                              ! Already on extract grid
            imin=xtr_iobounds(ng)%ILB_psi
            imax=xtr_iobounds(ng)%IUB_psi
            jmin=xtr_iobounds(ng)%JLB_psi
            jmax=xtr_iobounds(ng)%JUB_psi
#endif
          END IF
!
!       ---------------------
        CASE (r2dvar, r3dvar)               ! RHO points field
!       ---------------------
!
          IF (extract_flag.ge.0) THEN
            imin=iobounds(ng)%ILB_rho
            imax=iobounds(ng)%IUB_rho
            jmin=iobounds(ng)%JLB_rho
            jmax=iobounds(ng)%JUB_rho
#ifdef GRID_EXTRACT
          ELSE                              ! Already on extract grid
            imin=xtr_iobounds(ng)%ILB_rho
            imax=xtr_iobounds(ng)%IUB_rho
            jmin=xtr_iobounds(ng)%JLB_rho
            jmax=xtr_iobounds(ng)%JUB_rho
#endif
          END IF
!
!       ---------------------
        CASE (u2dvar, u3dvar)               ! U points field
!       ---------------------
!
          IF (extract_flag.ge.0) THEN
            imin=iobounds(ng)%ILB_u
            imax=iobounds(ng)%IUB_u
            jmin=iobounds(ng)%JLB_u
            jmax=iobounds(ng)%JUB_u
#ifdef GRID_EXTRACT
          ELSE                              ! Already on extract grid
            imin=xtr_iobounds(ng)%ILB_u
            imax=xtr_iobounds(ng)%IUB_u
            jmin=xtr_iobounds(ng)%JLB_u
            jmax=xtr_iobounds(ng)%JUB_u
#endif
          END IF
!
!       ---------------------
        CASE (v2dvar, v3dvar)               ! V points field
!       ---------------------
!
          IF (extract_flag.ge.0) THEN
            imin=iobounds(ng)%ILB_v
            imax=iobounds(ng)%IUB_v
            jmin=iobounds(ng)%JLB_v
            jmax=iobounds(ng)%JUB_v
#ifdef GRID_EXTRACT
          ELSE                              ! Already on extract grid
            imin=xtr_iobounds(ng)%ILB_v
            imax=xtr_iobounds(ng)%IUB_v
            jmin=xtr_iobounds(ng)%JLB_v
            jmax=xtr_iobounds(ng)%JUB_v
#endif
          END IF
!
!       ---------------------
        CASE DEFAULT                        ! RHO points field
!       ---------------------
!
          IF (extract_flag.ge.0) THEN
            imin=iobounds(ng)%ILB_rho
            imax=iobounds(ng)%IUB_rho
            jmin=iobounds(ng)%JLB_rho
            jmax=iobounds(ng)%JUB_rho
#ifdef GRID_EXTRACT
          ELSE                              ! Already on extract grid
            imin=xtr_iobounds(ng)%ILB_rho
            imax=xtr_iobounds(ng)%IUB_rho
            jmin=xtr_iobounds(ng)%JLB_rho
            jmax=xtr_iobounds(ng)%JUB_rho
#endif
          END IF
      END SELECT
!
      isize=imax-imin+1
      jsize=jmax-jmin+1
      ijsize=isize*jsize
 
#ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  If distributed-memory set-up, collect tile data from all spawned
!  nodes and store it into a global scratch 1D array, packed in column-
!  major order.
# ifdef MASKING
#  ifdef WRITE_WATER
!  Remove land points and pack water points into 1D-array.
#  else
!  Overwrite masked points with special value.
#  endif
# endif
!-----------------------------------------------------------------------
!
      IF (extract_flag.ge.0) THEN
        CALL mp_gather2d (ng, model, lbi, ubi, lbj, ubj,                &
     &                    tindex, gtype, ascl,                          &
# ifdef MASKING
     &                    amask,                                        &
# endif
     &                    adat, npts, awrk, landfill)
# ifdef GRID_EXTRACT
      ELSE                                  ! Already on extract grid
        CALL mp_gather2d_xtr (ng, model, lbi, ubi, lbj, ubj,            &
     &                        tindex, gtype, ascl,                      &
#  ifdef MASKING
     &                        amask,                                    &
#  endif
     &                        adat, npts, awrk, landfill)
# endif
      END IF
 
#else
!
!-----------------------------------------------------------------------
!  If serial or shared-memory applications and serial output, pack data
!  into a global 1D array in column-major order.
# ifdef MASKING
#  ifdef WRITE_WATER
!  Remove land points and pack water points into 1D-array.
#  else
!  Overwrite masked points with special value.
#  endif
# endif
!-----------------------------------------------------------------------
!
      IF (gtype.gt.0) THEN
        ic=0
        DO j=jmin,jmax
          DO i=imin,imax
            ic=ic+1
            awrk(ic)=adat(i,j)*ascl
# ifdef MASKING
            IF ((amask(i,j).eq.0.0_r8).and.landfill) THEN
              awrk(ic)=spval
            END IF
# endif
          END DO
        END DO
        npts=ic
# ifdef MASKING
      ELSE
        ic=0
        DO j=jmin,jmax
          DO i=imin,imax
            IF (amask(i,j).gt.0.0_r8) THEN
              ic=ic+1
              awrk(ic)=adat(i,j)*ascl
            END IF
          END DO
        END DO
        npts=ic
# endif
      END IF
#endif
!
!-----------------------------------------------------------------------
!  If there is no extracting data, set the start and total vectors
!  needed by the NetCDF library for writing.
!-----------------------------------------------------------------------
!
      IF (extract_flag.le.0) THEN
        IF (gtype.gt.0) THEN
          start(1)=1
          total(1)=isize
          start(2)=1
          total(2)=jsize
          start(3)=tindex
          total(3)=1
#ifdef MASKING
        ELSE
          start(1)=1
          total(1)=npts
          start(2)=tindex
          total(2)=1
#endif
        END IF
 
#ifdef GRID_EXTRACT
!
!-----------------------------------------------------------------------
!  Otherwise, extract field by decimation or interpolation.
!-----------------------------------------------------------------------
!
      ELSE IF (extract_flag.ge.1) THEN
        CALL extract_field (ng, model, tile,                            &
     &                      gtype, ifield, tindex,                      &
     &                      extract_flag,                               &
     &                      imin, imax, jmin, jmax,                     &
     &                      start, total,                               &
     &                      npts, awrk)
#endif
      END IF
!
      RETURN
      END SUBROUTINE pack_field2d
!
      SUBROUTINE pack_field3d (ng, model, tile,                         &
     &                         gtype, ifield, tindex,                   &
     &                         LandFill, Extract_Flag,                  &
     &                         LBi, UBi, LBj, UBj, LBk, UBk,            &
#ifdef MASKING
     &                         Amask,                                   &
#endif
     &                         Ascl, Adat,                              &
     &                         start, total, Npts, Awrk)
!
!=======================================================================
!                                                                      !
!  Packs 3D field data into 1D array to be written into output NetCDF  !
!  file elsewhere.                                                     !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng           Nested grid number (integer)                        !
!     model        Calling model identifier (integer)                  !
!     tile         Domain partition (integer)                          !
!     gtype        Staggered C-grid type (integer)                     !
!     ifield       Field metadata index (integer)                      !
!     tindex       Time record index to process (integer)              !
!     LandFill     Switch to set fill value in land areas (logical)    !
!     Extract_Flag Extraction flag interpolation/decimation (integer)  !
!     LBi          Field I-dimension Lower bound (integer)             !
!     UBi          Field I-dimension Upper bound (integer)             !
!     LBj          Field J-dimension Lower bound (integer)             !
!     UBj          Field J-dimension Upper bound (integer)             !
!     LBk          Field K-dimension Lower bound (integer)             !
!     UBk          Field K-dimension Upper bound (integer)             !
#ifdef MASKING
!     Amask        Land/Sea mask, if any (real 2D array)               !
#endif
!     Ascl         Factor to scale field after reading (real)          !
!     Adat         3D field to process (real)                          !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     start        Start index where the first of the data values will !
!                    be written along each dimension (integer)         !
!     total        Number of data values to be written along each      !
!                    dimension (integer)                               !
!     Npts         Number of points processed in Awrk.                 !
!     Awrk         Packed 3D field data (real)                         !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      logical, intent(in)  :: landfill
!
      integer, intent(in)  :: ng, model, tile
      integer, intent(in)  :: gtype, ifield, tindex
      integer, intent(in)  :: extract_flag
      integer, intent(in)  :: lbi, ubi, lbj, ubj, lbk, ubk
      integer, intent(out) :: npts
      integer, intent(out) :: start(:), total(:)
!
      real(dp), intent(in) :: ascl
!
#ifdef ASSUMED_SHAPE
# ifdef MASKING
      real(r8), intent(in) :: amask(lbi:,lbj:)
# endif
      real(r8), intent(in) :: adat(lbi:,lbj:,lbk:)
#else
# ifdef MASKING
      real(r8), intent(in) :: amask(lbi:ubi,lbj:ubj)
# endif
      real(r8), intent(in) :: adat(lbi:ubi,lbj:ubj,lbk:ubk)
#endif
      real(r8), intent(out) :: awrk(:)
!
!  Local variable declarations.
!
      integer :: i, j, k, ic
      integer :: imin, imax, jmin, jmax, koff
      integer :: isize, jsize, ksize, ijsize, mytype
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", pack_field3d"
!
!------------------------------------------------------------------------
!  Pack 3D field data.
!------------------------------------------------------------------------
!
!  Set first and last grid point according to staggered C-grid
!  classification. Set loops offsets.
!
      mytype=gtype
!
      SELECT CASE (abs(mytype))
!
!       ---------------------
        CASE (p3dvar)                       ! PSI points field
!       ---------------------
!
          imin=iobounds(ng)%ILB_psi
          imax=iobounds(ng)%IUB_psi
          jmin=iobounds(ng)%JLB_psi
          jmax=iobounds(ng)%JUB_psi
!
!       ---------------------
        CASE (r3dvar)                       ! RHO points field
!       ---------------------
!
          imin=iobounds(ng)%ILB_rho
          imax=iobounds(ng)%IUB_rho
          jmin=iobounds(ng)%JLB_rho
          jmax=iobounds(ng)%JUB_rho
!
!       ---------------------
        CASE (u3dvar)                       ! U points field
!       ---------------------
!
          imin=iobounds(ng)%ILB_u
          imax=iobounds(ng)%IUB_u
          jmin=iobounds(ng)%JLB_u
          jmax=iobounds(ng)%JUB_u
!
!       ---------------------
        CASE (v3dvar)                       ! V points field
!       ---------------------
!
          imin=iobounds(ng)%ILB_v
          imax=iobounds(ng)%IUB_v
          jmin=iobounds(ng)%JLB_v
          jmax=iobounds(ng)%JUB_v
!
!       ---------------------
        CASE DEFAULT                        ! RHO points field
!       ---------------------
!
          imin=iobounds(ng)%ILB_rho
          imax=iobounds(ng)%IUB_rho
          jmin=iobounds(ng)%JLB_rho
          jmax=iobounds(ng)%JUB_rho
      END SELECT
!
      isize=imax-imin+1
      jsize=jmax-jmin+1
      ksize=ubk-lbk+1
      ijsize=isize*jsize
!
      IF (lbk.eq.0) THEN
        koff=0
      ELSE
        koff=1
      END IF
 
#ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  If distributed-memory set-up, collect tile data from all spawned
!  nodes and store it into a global scratch 1D array, packed in column-
!  major order.
# ifdef MASKING
#  ifdef WRITE_WATER
!  Remove land points and pack water points into 1D-array.
#  else
!  Overwrite masked points with special value.
#  endif
# endif
!-----------------------------------------------------------------------
!
      CALL mp_gather3d (ng, model, lbi, ubi, lbj, ubj, lbk, ubk,        &
     &                  tindex, gtype, ascl,                            &
# ifdef MASKING
     &                  amask,                                          &
# endif
     &                  adat, npts, awrk, landfill)
#else
!
!-----------------------------------------------------------------------
!  If serial or shared-memory applications and serial output, pack data
!  into a global 1D array in column-major order.
# ifdef MASKING
#  ifdef WRITE_WATER
!  Remove land points and pack water points into 1D-array.
#  else
!  Overwrite masked points with special value.
#  endif
# endif
!-----------------------------------------------------------------------
!
      IF (gtype.gt.0) THEN
        ic=0
        DO k=lbk,ubk
          DO j=jmin,jmax
            DO i=imin,imax
              ic=ic+1
              awrk(ic)=adat(i,j,k)*ascl
# ifdef MASKING
              IF ((amask(i,j).eq.0.0_r8).and.landfill) THEN
                awrk(ic)=spval
              END IF
# endif
            END DO
          END DO
        END DO
        npts=ic
# ifdef MASKING
      ELSE
        ic=0
        DO k=lbk,ubk
          DO j=jmin,jmax
            DO i=imin,imax
              IF (amask(i,j).gt.0.0_r8) THEN
                ic=ic+1
                awrk(ic)=adat(i,j,k)*ascl
              END IF
            END DO
          END DO
        END DO
        npts=ic
# endif
      END IF
#endif
!
!-----------------------------------------------------------------------
!  If there is no extracting data, set the start and total vectors
!  needed by the NetCDF library for writing.
!-----------------------------------------------------------------------
!
      IF (extract_flag.le.0) THEN
        IF (gtype.gt.0) THEN
          start(1)=1
          total(1)=isize
          start(2)=1
          total(2)=jsize
          start(3)=1
          total(3)=ksize
          start(4)=tindex
          total(4)=1
#ifdef MASKING
        ELSE
          start(1)=1
          total(1)=npts
          start(2)=tindex
          total(2)=1
#endif
        END IF
 
#ifdef GRID_EXTRACT
!
!-----------------------------------------------------------------------
!  Otherwise, extract field by decimation or interpolation.
!-----------------------------------------------------------------------
!
      ELSE IF (extract_flag.ge.1) THEN
        CALL extract_field (ng, model, tile,                            &
     &                      gtype, ifield, tindex,                      &
     &                      extract_flag,                               &
     &                      imin, imax, jmin, jmax, lbk, ubk,           &
     &                      start, total,                               &
     &                      npts, awrk)
#endif
      END IF
 
#ifdef POSITIVE_ZERO
!
!-----------------------------------------------------------------------
!  Impose positive zero. Recall that IEEE allows both negative and
!  positive zeros, which affects bit-by-bit comparison betwee different
!  tile partitions.
!-----------------------------------------------------------------------
!
      DO ic=1,npts
        IF (abs(awrk(ic)).eq.0.0_r8) THEN
          awrk(ic)=0.0_r8
        END IF
      END DO
#endif
!
      RETURN
      END SUBROUTINE pack_field3d
!
      SUBROUTINE pack_field4d (ng, model, tile,                         &
     &                         gtype, ifield, tindex,                   &
     &                         LandFill, Extract_Flag,                  &
     &                         LBi, UBi, LBj, UBj, LBk, UBk, LBt, UBt,  &
     &                         fourth,                                  &
#ifdef MASKING
     &                         Amask,                                   &
#endif
     &                         Ascl, Adat,                              &
     &                         start, total, Npts, Awrk)
!
!=======================================================================
!                                                                      !
!  Packs 4D field data into 1D array to be written into output NetCDF  !
!  file elsewhere. The field data is processed by 3D slices to reduce  !
!  memory requirements.                                                !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng           Nested grid number (integer)                        !
!     model        Calling model identifier (integer)                  !
!     tile         Domain partition (integer)                          !
!     gtype        Staggered C-grid type (integer)                     !
!     ifield       Field metadata index (integer)                      !
!     tindex       Time record index to process (integer)              !
!     LandFill     Switch to set fill value in land areas (logical)    !
!     Extract_Flag Extraction flag interpolation/decimation (integer)  !
!     LBi          Field I-dimension Lower bound (integer)             !
!     UBi          Field I-dimension Upper bound (integer)             !
!     LBj          Field J-dimension Lower bound (integer)             !
!     UBj          Field J-dimension Upper bound (integer)             !
!     LBk          Field K-dimension Lower bound (integer)             !
!     UBk          Field K-dimension Upper bound (integer)             !
!     LBt          Field fourth-dimension Lower bound (integer)        !
!     UBt          Field fourth-dimension Upper bound (integer)        !
!     fourth       Fourth dimension index to process (integer)         !
#ifdef MASKING
!     Amask        Land/Sea mask, if any (real 2D array)               !
#endif
!     Ascl         Factor to scale field after reading (real)          !
!     Adat         4D field to process (real)                          !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     start        Start index where the first of the data values will !
!                    be written along each dimension (integer)         !
!     total        Number of data values to be written along each      !
!                    dimension (integer)                               !
!     Npts         Number of points processed in Awrk.                 !
!     Awrk         Packed 4D field data (real)                         !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      logical, intent(in)  :: landfill
!
      integer, intent(in)  :: ng, model, tile
      integer, intent(in)  :: gtype, ifield, tindex
      integer, intent(in)  :: extract_flag
      integer, intent(in)  :: lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt
      integer, intent(in)  :: fourth
      integer, intent(out) :: npts
      integer, intent(out) :: start(:), total(:)
!
      real(dp), intent(in) :: ascl
!
#ifdef ASSUMED_SHAPE
# ifdef MASKING
      real(r8), intent(in) :: amask(lbi:,lbj:)
# endif
      real(r8), intent(in) :: adat(lbi:,lbj:,lbk:,lbt:)
#else
# ifdef MASKING
      real(r8), intent(in) :: amask(lbi:ubi,lbj:ubj)
# endif
      real(r8), intent(in) :: adat(lbi:ubi,lbj:ubj,lbk:ubk,lbt:ubt)
#endif
      real(r8), intent(out) :: awrk(:)
!
!  Local variable declarations.
!
      integer :: i, j, k, ic
      integer :: imin, imax, jmin, jmax, kmin, kmax, koff, loff
      integer :: isize, jsize, ksize, ijsize, mytype
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", pack_field4d"
!
!------------------------------------------------------------------------
!  Pack 4D field data by 3D slices.
!------------------------------------------------------------------------
!
!  Set first and last grid point according to staggered C-grid
!  classification. Set loops offsets.
!
      mytype=gtype
!
      SELECT CASE (abs(mytype))
!
!       ---------------------
        CASE (p2dvar, p3dvar)               ! PSI points field
!       ---------------------
!
          imin=iobounds(ng)%ILB_psi
          imax=iobounds(ng)%IUB_psi
          jmin=iobounds(ng)%JLB_psi
          jmax=iobounds(ng)%JUB_psi
!
!       ---------------------
        CASE (r2dvar, r3dvar)               ! RHO points field
!       ---------------------
!
          imin=iobounds(ng)%ILB_rho
          imax=iobounds(ng)%IUB_rho
          jmin=iobounds(ng)%JLB_rho
          jmax=iobounds(ng)%JUB_rho
!
!       ---------------------
        CASE (u2dvar, u3dvar)               ! U points field
!       ---------------------
!
          imin=iobounds(ng)%ILB_u
          imax=iobounds(ng)%IUB_u
          jmin=iobounds(ng)%JLB_u
          jmax=iobounds(ng)%JUB_u
!
!       ---------------------
        CASE (v2dvar, v3dvar)               ! V points field
!       ---------------------
!
          imin=iobounds(ng)%ILB_v
          imax=iobounds(ng)%IUB_v
          jmin=iobounds(ng)%JLB_v
          jmax=iobounds(ng)%JUB_v
!
!       ---------------------
        CASE DEFAULT                        ! RHO points field
!       ---------------------
!
          imin=iobounds(ng)%ILB_rho
          imax=iobounds(ng)%IUB_rho
          jmin=iobounds(ng)%JLB_rho
          jmax=iobounds(ng)%JUB_rho
      END SELECT
!
      isize=imax-imin+1
      jsize=jmax-jmin+1
      ksize=ubk-lbk+1
      ijsize=isize*jsize
!
      IF (lbk.eq.0) THEN
        koff=0
      ELSE
        koff=1
      END IF
!
      IF (lbt.eq.0) THEN
        loff=1
      ELSE
        loff=0
      END IF
 
#ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  If distributed-memory set-up, collect tile data from all spawned
!  nodes and store it into a global scratch 1D array, packed in column-
!  major order.
# ifdef MASKING
#  ifdef WRITE_WATER
!  Remove land points and pack water points into 1D-array.
#  else
!  Overwrite masked points with special value.
#  endif
# endif
!-----------------------------------------------------------------------
!
!  Process the data by 3D slices.
!
      CALL mp_gather3d (ng, model, lbi, ubi, lbj, ubj, lbk, ubk,        &
     &                  tindex, gtype, ascl,                            &
# ifdef MASKING
     &                  amask,                                          &
# endif
     &                  adat(:,:,:,fourth), npts, awrk, landfill)
#else
!
!-----------------------------------------------------------------------
!  If serial or shared-memory applications and serial output, pack data
!  into a global 1D array in column-major order.
# ifdef MASKING
#  ifdef WRITE_WATER
!  Remove land points and pack water points into 1D-array.
#  else
!  Overwrite masked points with special value.
#  endif
# endif
!-----------------------------------------------------------------------
!
!  The data is processed in 3D slices.
!
      IF (gtype.gt.0) THEN
        ic=0
        npts=ijsize*ksize
        DO k=lbk,ubk
          DO j=jmin,jmax
            DO i=imin,imax
              ic=ic+1
              awrk(ic)=adat(i,j,k,fourth)*ascl
# ifdef MASKING
              IF ((amask(i,j).eq.0.0_r8).and.landfill) THEN
                awrk(ic)=spval
              END IF
# endif
            END DO
          END DO
        END DO
# ifdef MASKING
      ELSE
        npts=0
        DO k=lbk,ubk
          DO j=jmin,jmax
            DO i=imin,imax
              IF (amask(i,j).gt.0.0_r8) THEN
                npts=npts+1
                awrk(npts)=adat(i,j,k,fourth)*ascl
              END IF
            END DO
          END DO
        END DO
# endif
      END IF
#endif
!
!-----------------------------------------------------------------------
!  If there is no extracting data, set the start and total vectors
!  needed by the NetCDF library for writing.
!-----------------------------------------------------------------------
!
      IF (extract_flag.le.0) THEN
        IF (gtype.gt.0) THEN
          start(1)=1
          total(1)=isize
          start(2)=1
          total(2)=jsize
          start(3)=1
          total(3)=ksize
          start(4)=fourth+loff
          total(4)=1
          start(5)=tindex
          total(5)=1
#ifdef MASKING
        ELSE
          start(1)=1+(fourth+loff-1)*npts
          total(1)=npts
          start(2)=tindex
          total(2)=1
#endif
        END IF
 
#ifdef GRID_EXTRACT
!
!-----------------------------------------------------------------------
!  Otherwise, extract field by decimation or interpolation.
!-----------------------------------------------------------------------
!
      ELSE IF (extract_flag.ge.1) THEN
        CALL extract_field (ng, model, tile,                            &
     &                      gtype, ifield, tindex,                      &
     &                      extract_flag,                               &
     &                      imin, imax, jmin, jmax, lbk, ubk,           &
     &                      start, total,                               &
     &                      npts, awrk)
#endif
      END IF
 
#ifdef POSITIVE_ZERO
!
!-----------------------------------------------------------------------
!  Impose positive zero. Recall that IEEE allows both negative and
!  positive zeros, which affects bit-by-bit comparison betwee different
!  tile partitions.
!-----------------------------------------------------------------------
!
      DO ic=1,npts
        IF (abs(awrk(ic)).eq.0.0_r8) THEN
          awrk(ic)=0.0_r8
        END IF
      END DO
#endif
!
      RETURN
      END SUBROUTINE pack_field4d
!
      END MODULE pack_field_mod