nf_fwrite4d.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE nf_fwrite4d_mod
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This function writes out a generic floating point 4D array into an  !
!  output NetCDF file.                                                 !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng           Nested grid number                                  !
!     model        Calling model identifier                            !
!     ncid         NetCDF file ID                                      !
!     ifield       Field metadata index (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
!     ncvarid      NetCDF variable ID                                  !
#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 write                   !
!     gtype        Grid type. If negative, only write water points     !
#if defined PIO_LIB && defined DISTRIBUTE
!or   pioDesc      IO data decomposition descriptor, TYPE(IO_desc_t)   !
#endif
!     LBi          I-dimension Lower bound                             !
!     UBi          I-dimension Upper bound                             !
!     LBj          J-dimension Lower bound                             !
!     UBj          J-dimension Upper bound                             !
!     LBk          K-dimension Lower bound                             !
!     UBk          K-dimension Upper bound                             !
!     LBt          Time-dimension Lower bound                          !
!     UBt          Time-dimension Upper bound                          !
!     Amask        land/Sea mask, if any (real)                        !
!     Ascl         Factor to scale field before writing (real)         !
!     Adat         Field to write out (real)                           !
!     SetFillVal   Logical switch to set fill value in land areas      !
!                    (OPTIONAL)                                        !
!     ExtractField Field extraction flag (integer, OPTIONAL)           !
!                    ExtractField = 0   no extraction                  !
!                    ExtractField = 1   extraction by intrpolation     !
!                    ExtractField > 1   extraction by decimation       !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     status       Error flag (integer)                                !
!     MinValue     Minimum value (real, OPTIONAL)                      !
!     MaxValue     Maximum value (real, OPTIONAL)                      !
!                                                                      !
#ifdef POSITIVE_ZERO
!  Starting F95 zero values can be signed (-0 or +0) following the     !
!  IEEE 754 floating point standard.  This may produce different       !
!  output data in serial and parallel applications. Since comparing    !
!  serial and parallel output is essential for tracking parallel       !
!  partition  bugs, "positive zero" is enforced.                       !
!                                                                      !
#endif
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_ncparam
      USE mod_scalars
!
      USE pack_field_mod,  ONLY : pack_field4d
!
      implicit none
!
      INTERFACE nf_fwrite4d
        MODULE PROCEDURE nf90_fwrite4d
#if defined PIO_LIB && defined DISTRIBUTE
        MODULE PROCEDURE pio_fwrite4d
#endif
      END INTERFACE nf_fwrite4d
!
      CONTAINS
 
#if defined PARALLEL_IO && defined DISTRIBUTE
!
!***********************************************************************
      FUNCTION nf90_fwrite4d (ng, model, ncid, ifield,                  &
     &                        ncvarid, tindex, gtype,                   &
     &                        LBi, UBi, LBj, UBj, LBk, UBk, LBt, UBt,   &
     &                        Ascl,                                     &
# ifdef MASKING
     &                        Amask,                                    &
# endif
     &                        Adat,                                     &
     &                        SetFillVal,                               &
     &                        ExtractField,                             &
     &                        MinValue, MaxValue) RESULT (status)
!***********************************************************************
!
      USE mod_netcdf
!
# ifdef MASKING
      USE distribute_mod, ONLY : mp_collect
# endif
      USE get_bounds_mod, ONLY : tile_bounds_1d
!
!  Imported variable declarations.
!
      logical, intent(in), optional :: setfillval
!
      integer, intent(in) :: ng, model, ncid, ncvarid, tindex, gtype
      integer, intent(in) :: ifield
      integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt
!
      integer, intent(in), optional :: extractfield
!
      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), optional :: minvalue
      real(r8), intent(out), optional :: maxvalue
!
!  Local variable declarations.
!
      logical :: landfill
!
      integer :: extract_flag
      integer :: i, ic, j, jc, k, kc, l, lc, npts
      integer :: imin, imax, jmin, jmax, kmin, kmax
      integer :: ioff, joff, koff, loff
      integer :: istr, iend
      integer :: ilen, isize, jlen, jsize, ijlen, ijklen, klen, llen
      integer :: mytype
      integer :: status
 
      integer, dimension(5) :: start, total
!
      real(r8), parameter :: inival = 0.0_r8
 
      real(r8), allocatable :: awrk(:)
!
!-----------------------------------------------------------------------
!  Set starting and ending indices to process.
!-----------------------------------------------------------------------
!
      status=nf90_noerr
!
!  Set first and last grid point according to staggered C-grid
!  classification.
!
      mytype=gtype
!
      SELECT CASE (abs(mytype))
        CASE (p2dvar, p3dvar)
          imin=bounds(ng)%Istr (myrank)
          imax=bounds(ng)%IendR(myrank)
          jmin=bounds(ng)%Jstr (myrank)
          jmax=bounds(ng)%JendR(myrank)
          isize=iobounds(ng)%xi_psi
          jsize=iobounds(ng)%eta_psi
        CASE (r2dvar, r3dvar)
          imin=bounds(ng)%IstrR(myrank)
          imax=bounds(ng)%IendR(myrank)
          jmin=bounds(ng)%JstrR(myrank)
          jmax=bounds(ng)%JendR(myrank)
          isize=iobounds(ng)%xi_rho
          jsize=iobounds(ng)%eta_rho
        CASE (u2dvar, u3dvar)
          imin=bounds(ng)%Istr (myrank)
          imax=bounds(ng)%IendR(myrank)
          jmin=bounds(ng)%JstrR(myrank)
          jmax=bounds(ng)%JendR(myrank)
          isize=iobounds(ng)%xi_u
          jsize=iobounds(ng)%eta_u
        CASE (v2dvar, v3dvar)
          imin=bounds(ng)%IstrR(myrank)
          imax=bounds(ng)%IendR(myrank)
          jmin=bounds(ng)%Jstr (myrank)
          jmax=bounds(ng)%JendR(myrank)
          isize=iobounds(ng)%xi_v
          jsize=iobounds(ng)%eta_v
        CASE DEFAULT
          imin=bounds(ng)%IstrR(myrank)
          imax=bounds(ng)%IendR(myrank)
          jmin=bounds(ng)%JstrR(myrank)
          jmax=bounds(ng)%JendR(myrank)
          isize=iobounds(ng)%xi_rho
          jsize=iobounds(ng)%eta_rho
      END SELECT
!
      ilen=imax-imin+1
      jlen=jmax-jmin+1
      klen=ubk-lbk+1
      llen=ubt-lbt+1
!
!  Set switch to replace land areas with fill value, spval.
!
# ifdef MASKING
      IF (PRESENT(setfillval)) THEN
        landfill=setfillval
      ELSE
        landfill=tindex.gt.0
      END IF
# else
      landfill=.false.
# endif
!
!  If appropriate, set the field extraction flag to the provided grid
!  geometry through interpolation or decimation.
!
      IF (PRESENT(extractfield)) THEN
        extract_flag=extractfield
      ELSE
        extract_flag=0
      END IF
!
!-----------------------------------------------------------------------
!  Parallel I/O: Pack tile data into 1D array in column-major order.
# ifdef MASKING
!                Overwrite masked points with special value.
# endif
!-----------------------------------------------------------------------
!
      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
 
        IF (lbk.eq.0) THEN
          koff=1
        ELSE
          koff=0
        END IF
 
        IF (lbt.eq.0) THEN
          loff=1
        ELSE
          loff=0
        END IF
!
!  Allocate and initialize scratch work array.
!
        ijlen=ilen*jlen
        ijklen=ijlen*klen
        npts=ijklen*llen
 
        IF (.not.allocated(awrk)) THEN
          allocate ( awrk(npts) )
          awrk=inival
        END IF
!
!  Pack and scale tile data.
!
        ic=0
        DO l=lbt,ubt
          DO k=lbk,ubk
            DO j=jmin,jmax
              DO i=imin,imax
                ic=ic+1
                awrk(ic)=adat(i,j,k,l)*ascl
# ifdef POSITIVE_ZERO
                IF (abs(awrk(ic)).eq.0.0_r8) THEN
                  awrk(ic)=0.0_r8                 ! impose positive zero
                END IF
# endif
# 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
        END DO
!
!  Write out data: all parallel nodes write their own packed tile data.
!
        start(1)=imin+ioff
        total(1)=ilen
        start(2)=jmin+joff
        total(2)=jlen
        start(3)=lbk+koff
        total(3)=klen
        start(4)=lbt+loff
        total(4)=llen
        start(5)=tindex
        total(5)=1
 
        status=nf90_put_var(ncid, ncvarid, awrk, start, total)
      END IF
 
# if defined WRITE_WATER && defined MASKING
!
!-----------------------------------------------------------------------
!  Parallel I/O: Remove land points and pack tile data into 1D array.
!-----------------------------------------------------------------------
!
      IF (gtype.lt.0) THEN
!
!  Set offsets due array packing into 1D array in column-major order.
!
        SELECT CASE (abs(mytype))
          CASE (p2dvar, p3dvar)
            ioff=0
            joff=1
          CASE (r2dvar, r3dvar)
            ioff=1
            joff=0
          CASE (u2dvar, u3dvar)
            ioff=0
            joff=0
          CASE (v2dvar, v3dvar)
            ioff=1
            joff=1
          CASE DEFAULT
            ioff=1
            joff=0
        END SELECT
 
        IF (lbk.eq.0) THEN
          koff=0
        ELSE
          koff=1
        END IF
 
        IF (lbt.eq.0) THEN
          loff=0
        ELSE
          loff=1
        END IF
!
!  Allocate and initialize scratch work array.
!
        ijlen=isize*jsize
        ijklen=ijlen*klen
        npts=ijklen*llen
 
        IF (.not.allocated(awrk)) THEN
          allocate ( awrk(npts) )
          awrk=inival
        END IF
!
!  Scale and gather data from all spawned nodes. Store data into a 1D
!  global array, packed in column-major order. Flag land point with
!  special value.
!
        DO l=lbt,ubt
          lc=(l-loff)*ijklen
          DO k=lbk,ubk
            kc=(k-koff)*ijlen+lc
            DO j=jmin,jmax
              jc=(j-joff)*isize+kc
              DO i=imin,imax
                ic=i+ioff+jc
                awrk(ic)=adat(i,j,k,l)*ascl
#  ifdef POSITIVE_ZERO
                IF (abs(awrk(ic)).eq.0.0_r8) THEN
                  awrk(ic)=0.0_r8                 ! impose positive zero
                END IF
#  endif
                IF (amask(i,j).eq.0.0_r8) THEN
                  awrk(ic)=spval
                END IF
              END DO
            END DO
          END DO
        END DO
!
!  Global reduction of work array.
!
        CALL mp_collect (ng, model, npts, inival, awrk)
!
!  Remove land points.
!
        ic=0
        DO i=1,npts
          IF (awrk(i).lt.spval) THEN
            ic=ic+1
            awrk(ic)=awrk(i)
          END IF
        END DO
        npts=ic
!
!  Write out data: all parallel nodes write a section of the packed
!                  data.
!
        CALL tile_bounds_1d (ng, myrank, npts, istr, iend)
!
        start(1)=istr
        total(1)=iend-istr+1
        start(2)=tindex
        total(2)=1
 
        status=nf90_put_var(ncid, ncvarid, awrk(istr:), start, total)
      END IF
# endif
!
!-----------------------------------------------------------------------
!  Deallocate scratch work array.
!-----------------------------------------------------------------------
!
      IF (allocated(awrk)) THEN
        deallocate (awrk)
      END IF
!
      RETURN
      END FUNCTION nf90_fwrite4d
 
#else
 
!
!***********************************************************************
      FUNCTION nf90_fwrite4d (ng, model, ncid, ifield,                  &
     &                        ncvarid, tindex, gtype,                   &
     &                        LBi, UBi, LBj, UBj, LBk, UBk, LBt, UBt,   &
     &                        Ascl,                                     &
# ifdef MASKING
     &                        Amask,                                    &
# endif
     &                        Adat,                                     &
     &                        SetFillVal,                               &
     &                        ExtractField,                             &
     &                        MinValue, MaxValue) RESULT (status)
!***********************************************************************
!
      USE mod_netcdf
 
# ifdef DISTRIBUTE
!
      USE distribute_mod, ONLY : mp_bcasti
#  ifdef INLINE_2DIO
      USE distribute_mod, ONLY : mp_gather2d
#  else
      USE distribute_mod, ONLY : mp_gather3d
#  endif
# endif
# ifdef OUTPUT_STATS
      USE stats_mod,      ONLY : stats_4dfld
# endif
!
!  Imported variable declarations.
!
      logical, intent(in), optional :: setfillval
!
      integer, intent(in) :: ng, model, ncid, ncvarid, tindex, gtype
      integer, intent(in) :: ifield
      integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt
!
      integer, intent(in), optional :: extractfield
!
      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), optional :: minvalue
      real(r8), intent(out), optional :: maxvalue
!
!  Local variable declarations.
!
      logical :: landfill
!
      integer :: extract_flag
      integer :: i, fourth, npts, tile
      integer :: status
 
      integer, dimension(5) :: start, total
!
      real(r8), dimension((Lm(ng)+2)*(Mm(ng)+2)*(UBk-LBk+1)) :: awrk
 
# ifdef OUTPUT_STATS
!
      TYPE (t_stats) :: stats
# endif
!
!-----------------------------------------------------------------------
!  Set starting and ending indices to process.
!-----------------------------------------------------------------------
!
      status=nf90_noerr
!
!  Set parallel tile.
!
# ifdef DISTRIBUTE
      tile=myrank
# else
      tile=-1
# endif
!
!  Set switch to replace land areas with fill value, spval.
!
# ifdef MASKING
      IF (PRESENT(setfillval)) THEN
        landfill=setfillval
      ELSE
        landfill=tindex.gt.0
      END IF
# else
      landfill=.false.
# endif
!
!  If appropriate, set the field extraction flag to the provided grid
!  geometry through interpolation or decimation.
!
      IF (PRESENT(extractfield)) THEN
        extract_flag=extractfield
      ELSE
        extract_flag=0
      END IF
 
# ifdef OUTPUT_STATS
!
!  Initialize output variable statistics structure.
!
      stats % checksum=0_i8b
      stats % count=0
      stats % min=spval
      stats % max=-spval
      stats % avg=0.0_r8
      stats % rms=0.0_r8
# endif
!
!  If appropriate, initialize minimum and maximum processed values.
!
      IF (PRESENT(minvalue)) THEN
        minvalue=spval
        maxvalue=-spval
      END IF
!
!  Initialize local array to avoid denormalized numbers. This
!  facilitates processing and debugging.
!
      awrk=0.0_r8
!
!-----------------------------------------------------------------------
!  Pack 4D field data into 1D array. Proccess fourth dimension record
!  by record.
!-----------------------------------------------------------------------
!
      dim4_loop : DO fourth=lbt,ubt
        CALL 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)
!
!-----------------------------------------------------------------------
!  If applicable, compute output field minimum and maximum values.
!-----------------------------------------------------------------------
!
        IF (PRESENT(minvalue)) THEN
          IF (outthread) THEN
            DO i=1,npts
              IF (abs(awrk(i)).lt.spval) THEN
                minvalue=min(minvalue,awrk(i))
                maxvalue=max(maxvalue,awrk(i))
              END IF
            END DO
          END IF
        END IF
!
!-----------------------------------------------------------------------
!  Write output buffer into NetCDF file.
!-----------------------------------------------------------------------
!
        IF (outthread) THEN
          status=nf90_put_var(ncid, ncvarid, awrk, start, total)
        END IF
      END DO dim4_loop
 
# ifdef OUTPUT_STATS
!
!-----------------------------------------------------------------------
!  Compute and report output field statistics.
!-----------------------------------------------------------------------
!
      CALL stats_4dfld (ng, tile, model, gtype, stats,                  &
     &                  extract_flag,                                   &
     &                  lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt,         &
     &                  adat,                                           &
#  ifdef MASKING
     &                  fmask = amask,                                  &
#  endif
     &                  debug = .false.)
      IF (outthread) THEN
        WRITE (stdout,10) trim(vname(1,ifield)), stats%min, stats%max,  &
     &                    stats%checksum
  10    FORMAT (4x,'- ',a,':',t35,'Min = ',1p,e15.8,',  Max = ',        &
     &          1p,e15.8,',  CheckSum = ',i0)
      END IF
# endif
# ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Broadcast IO error flag to all nodes.
!-----------------------------------------------------------------------
!
      CALL mp_bcasti (ng, model, status)
# endif
!
      RETURN
      END FUNCTION nf90_fwrite4d
#endif
 
#if defined PIO_LIB && defined DISTRIBUTE
!
!***********************************************************************
      FUNCTION pio_fwrite4d (ng, model, pioFile, ifield,                &
     &                       pioVar, tindex, pioDesc,                   &
     &                       LBi, UBi, LBj, UBj, LBk, UBk, LBt, UBt,    &
     &                       Ascl,                                      &
# ifdef MASKING
     &                       Amask,                                     &
# endif
     &                       Adat, SetFillVal,                          &
     &                       MinValue, MaxValue) RESULT (status)
!***********************************************************************
!
      USE mod_pio_netcdf
!
      USE distribute_mod, ONLY : mp_reduce
# ifdef OUTPUT_STATS
      USE stats_mod,      ONLY : stats_4dfld
# endif
!
!  Imported variable declarations.
!
      logical, intent(in), optional :: setfillval
!
      integer, intent(in) :: ng, model, tindex
      integer, intent(in) :: ifield
      integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt
!
      integer, intent(in), optional :: extractfield
!
      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), optional :: minvalue
      real(r8), intent(out), optional :: maxvalue
!
      TYPE (file_desc_t), intent(inout) :: piofile
      TYPE (io_desc_t),   intent(inout) :: piodesc
      TYPE (my_vardesc),  intent(inout) :: piovar
!
!  Local variable declarations.
!
      logical :: landfill, lminmax
 
      logical, pointer :: lwater(:,:,:,:)
!
      integer :: extract_flag
      integer :: i, j, k, l, tile
      integer :: imin, imax, jmin, jmax
      integer :: cgrid, dkind, ghost, gtype
      integer :: status
 
      integer, dimension(5) :: start, total
!
      real(r8), dimension(2) :: rbuffer
 
      real(r4), pointer :: awrk4(:,:,:,:)
      real(r8), pointer :: awrk8(:,:,:,:)
 
# ifdef OUTPUT_STATS
!
      TYPE (t_stats) :: stats
# endif
!
      character (len= 3), dimension(2) :: op_handle
!
!-----------------------------------------------------------------------
!  Set starting and ending indices to process.
!-----------------------------------------------------------------------
!
      status=pio_noerr
!
!  Set first and last grid point according to staggered C-grid
!  classification. Set loops offsets.
!
      ghost=0
      dkind=piovar%dkind
      gtype=piovar%gtype
!
      SELECT CASE (abs(gtype))
        CASE (p3dvar)
          cgrid=1
        CASE (l3dvar, l4dvar, r3dvar)
          cgrid=2
        CASE (u3dvar)
          cgrid=3
        CASE (v3dvar)
          cgrid=4
        CASE DEFAULT
          cgrid=2
      END SELECT
!
      imin=bounds(ng)%Imin(cgrid,ghost,myrank)
      imax=bounds(ng)%Imax(cgrid,ghost,myrank)
      jmin=bounds(ng)%Jmin(cgrid,ghost,myrank)
      jmax=bounds(ng)%Jmax(cgrid,ghost,myrank)
!
!  Set switch to compute minimum and maximum values.
!
      IF (PRESENT(minvalue)) THEN
        lminmax=.true.
        IF (.not.associated(lwater)) THEN
          allocate ( lwater(lbi:ubi,lbj:ubj,lbk:ubk,lbt:ubt) )
          lwater=.true.
        END IF
      ELSE
        lminmax=.false.
      END IF
!
!  Set switch to replace land areas with fill value, spval.
!
# ifdef MASKING
      IF (PRESENT(setfillval)) THEN
        landfill=setfillval
      ELSE
        landfill=tindex.gt.0
      END IF
# else
      landfill=.false.
# endif
!
!  If appropriate, set the field extraction flag to the provided grid
!  geometry through interpolation or decimation.
!
      IF (PRESENT(extractfield)) THEN
        extract_flag=extractfield
      ELSE
        extract_flag=0
      END IF
 
# ifdef OUTPUT_STATS
!
!  Initialize output variable statistics structure.
!
      stats % checksum=0_i8b
      stats % count=0
      stats % min=spval
      stats % max=-spval
      stats % avg=0.0_r8
      stats % rms=0.0_r8
# endif
!
!-----------------------------------------------------------------------
!  Write out data into NetCDF file.
!-----------------------------------------------------------------------
!
!  Allocate, initialize and load data into local array used for
!  writing. Overwrite masked points with special value.
!
      IF (dkind.eq.pio_double) THEN                 ! double precision
        IF (.not.associated(awrk8)) THEN
          allocate ( awrk8(lbi:ubi,lbj:ubj,lbk:ubk,lbt:ubt) )
          awrk8=0.0_r8
        END IF
!
        DO l=lbt,ubt
          DO k=lbk,ubk
            DO j=jmin,jmax
              DO i=imin,imax
                awrk8(i,j,k,l)=adat(i,j,k,l)*ascl
# ifdef MASKING
                IF ((amask(i,j).eq.0.0_r8).and.landfill) THEN
                  awrk8(i,j,k,l)=spval
                  IF (lminmax) lwater(i,j,k,l)=.false.
                END IF
# endif
              END DO
            END DO
          END DO
        END DO
        IF (lminmax) THEN
          rbuffer(1)=minval(awrk8, mask=lwater)
          rbuffer(2)=maxval(awrk8, mask=lwater)
        END IF
      ELSE                                          ! single precision
        IF (.not.associated(awrk4)) THEN
          allocate ( awrk4(lbi:ubi,lbj:ubj,lbk:ubk,lbt:ubt) )
          awrk4=0.0_r8
        END IF
!
        DO l=lbt,ubt
          DO k=lbk,ubk
            DO j=jmin,jmax
              DO i=imin,imax
                awrk4(i,j,k,l)=real(adat(i,j,k,l)*ascl, r4)
# ifdef MASKING
                IF ((amask(i,j).eq.0.0_r8).and.landfill) THEN
                  awrk4(i,j,k,l)=real(spval, r4)
                END IF
# endif
              END DO
            END DO
          END DO
        END DO
        IF (lminmax) THEN
          rbuffer(1)=real(minval(awrk4, mask=lwater),r8)
          rbuffer(2)=real(maxval(awrk4, mask=lwater),r8)
        END IF
      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
!
!  Write out data into NetCDF.
!
      IF (dkind.eq.pio_double) THEN                 ! double precision
        CALL pio_write_darray (piofile,                                 &
     &                         piovar%vd,                               &
     &                         piodesc,                                 &
     &                         awrk8(imin:imax,jmin:jmax,               &
     &                               lbk:ubk,lbt:ubt),                  &
     &                         status)
      ELSE                                          ! single precision
        CALL pio_write_darray (piofile,                                 &
     &                         piovar%vd,                               &
     &                         piodesc,                                 &
     &                         awrk4(imin:imax,jmin:jmax,               &
     &                               lbk:ubk,lbt:ubt),                  &
     &                         status)
      END IF
 
# ifdef OUTPUT_STATS
!
!-----------------------------------------------------------------------
!  Compute and report output field statistics.
!-----------------------------------------------------------------------
!
#  ifdef DISTRIBUTE
      tile=myrank
#  else
      tile=-1
#  endif
      CALL stats_4dfld (ng, tile, model, gtype, stats,                  &
     &                  extract_flag,                                   &
     &                  lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt,         &
     &                  adat,                                           &
#  ifdef MASKING
     &                  fmask = amask,                                  &
#  endif
     &                  debug = .false.)
      IF (outthread) THEN
        WRITE (stdout,10) trim(vname(1,ifield)), stats%min, stats%max,  &
     &                    stats%checksum
  10    FORMAT (4x,'- ',a,':',t35,'Min = ',1p,e15.8,',  Max = ',        &
     &          1p,e15.8,',  CheckSum = ',i0)
      END IF
# endif
!
!-----------------------------------------------------------------------
!  If applicable, compute global minimum and maximum values.
!-----------------------------------------------------------------------
!
      IF (lminmax) THEN
        op_handle(1)='MIN'
        op_handle(2)='MAX'
        CALL mp_reduce (ng, model, 2, rbuffer, op_handle)
        minvalue=rbuffer(1)
        maxvalue=rbuffer(2)
        IF (associated(lwater)) deallocate (lwater)
      END IF
!
!  Deallocate local array.
!
      IF (dkind.eq.pio_double) THEN
        IF (associated(awrk8)) deallocate (awrk8)
      ELSE
        IF (associated(awrk4)) deallocate (awrk4)
      END IF
!
      RETURN
      END FUNCTION pio_fwrite4d
#endif
!
      END MODULE nf_fwrite4d_mod