set_pio.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE set_pio_mod
#if defined PIO_LIB && defined DISTRIBUTE
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This module contains several routines initializes and configures    !
!  Parallel-IO (PIO) in ROMS.  It uses the SCORPIO library which is    !
!  based on the PIO library developed at NCAR.                         !
!                                                                      !
!=======================================================================
!
      USE pio
!
      USE mod_kinds
      USE mod_param
      USE mod_parallel
      USE mod_pio_netcdf
      USE mod_iounits
# if defined PROPAGATOR && defined CHECKPOINTING
      USE mod_storage
# endif
!
# if defined SSH_TIDES || defined UV_TIDES
      USE mod_stepping, ONLY : ntc
# endif
      USE strings_mod,  ONLY : founderror
!
      implicit none
!
      PUBLIC  :: initialize_pio
      PUBLIC  :: finalize_pio
      PUBLIC  :: field_iodecomp
# if defined PROPAGATOR && defined CHECKPOINTING
      PUBLIC  :: state_iodecomp
# endif
      PUBLIC  :: set_iodecomp
# ifdef ASYNCHRONOUS_SCORPIO
      PUBLIC  :: set_pio_async
# endif
!
      CONTAINS
!
      SUBROUTINE initialize_pio
!
!***********************************************************************
!                                                                      !
!  Initializes the PIO subsystem. It sets PIO decomposition for all    !
!  ROMS variables.                                                     !
!                                                                      !
!***********************************************************************
!
!  Local variable declarations.
!
# ifdef ASYNCHRONOUS_PIO
      logical, allocatable :: lranks(:)
!
      integer :: computesize
 
      integer, allocatable :: compute_comm(:,:),  io_comm(:)
      integer, allocatable :: compute_ranks(:,:), io_ranks(:)
# endif
      integer :: myerror
      integer :: i, ic, ng
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", initialize_pio"
!
!-----------------------------------------------------------------------
!  Initialize PIO and get IO system descriptor.  It uses collective
!  communicatios.
!-----------------------------------------------------------------------
!
      IF (.not.allocated(piosystem)) THEN
        allocate ( piosystem(npiocomps,ngrids) )
      END IF
!
!  Set PIO internal level of debug information. The default value is 0,
!  allowed values 0-6.
!
      IF (pio_debug.gt.0) THEN
        CALL pio_setdebuglevel (pio_debug)
      END IF
 
# if defined ASYNCHRONOUS_PIO
!
!-----------------------------------------------------------------------
!  If NCAR/UNIDATA PIO library, set the ranks of the computational and
!  IO processes with respect the initial (peer) communicator. Notice
!  that the communicator is natively split inside "PIO_init". Therefore,
!  the I/O processes do not return from the call. Instead go to an
!  internal loop and wait to receive further instructions from the
!  computational processes. I think that this is a better strategy.
!-----------------------------------------------------------------------
!
!  Set peer communicator as the initial ROMS communicator.
!
      peer_comm_world=ocn_comm_world
      peersize=numthreads
      peerrank=myrank
      computesize=peersize-pio_numiotasks
!
!  Allocate and initialize local arrays.
!
      IF (.not.allocated(compute_comm)) THEN
        allocate ( compute_comm(npiocomps,ngrids) )
        compute_comm=mpi_comm_null
      END IF
      IF (.not.allocated(io_comm)) THEN
        allocate ( io_comm(ngrids) )
        io_comm=mpi_comm_null
      END IF
!
      IF (.not.allocated(lranks)) THEN
        allocate ( lranks(0:peersize-1) )
        lranks(0:peersize-1)=.true.
      END IF
      IF (.not.allocated(compute_ranks)) THEN
        allocate ( compute_ranks(computesize,npiocomps) )
        compute_ranks=-1
      END IF
      IF (.not.allocated(io_ranks)) THEN
        allocate ( io_ranks(pio_numiotasks) )
        io_ranks=-1
      END IF
!
!  Set the ranks of the dedicated I/O processes with respect the peer
!  communicator.
!
      ic=pio_base
      io_ranks(1)=ic
      lranks(ic)=.false.
      DO i=2,pio_numiotasks
        IF ((ic+pio_stride).le.(peersize-1)) THEN
          ic=ic+pio_stride
          io_ranks(i)=ic
          lranks(ic)=.false.
        END IF
      END DO
      WRITE (cioranks,'(*(i0,1x))') io_ranks
!
!  Set the ranks of the computational processes with respect the peer
!  communicator.
!
      ic=0
      DO i=0,peersize-1
        IF (lranks(i)) THEN
          ic=ic+1
          compute_ranks(ic,ipioroms)=i
        END IF
      END DO
      WRITE (ccompranks,'(*(i0,1x))') compute_ranks
!
!  Create a new IO system for asynchronous or synchronous I/O.  The
!  asynchronous I/O is only possible PIO type ("io_pio") files.
!
!!    IF ((inp_lib.eq.io_pio).and.(out_lib.eq.io_pio)) THEN
      IF (out_lib.eq.io_pio) THEN
        DO ng=1,ngrids
          CALL pio_init (piosystem(ipioroms:,ng),                       &
     &                   peer_comm_world,                               &
     &                   (/computesize/),                               &
     &                   compute_ranks,                                 &
     &                   io_ranks,                                      &
     &                   pio_rearranger,                                &
     &                   compute_comm(ipioroms:,ng),                    &
     &                   io_comm(ng))
        END DO
!
!  Initialize ROMS kernel communicators.
!
        ng=1
        ocn_comm_world=compute_comm(ipioroms,ng)
        io_comm_world =io_comm(ng)
!
!  Reset ROMS communicator parameters.
!
        CALL mpi_comm_size (ocn_comm_world, numthreads, myerror)
        CALL mpi_comm_rank (ocn_comm_world, myrank, myerror)
!
        master=myrank.eq.mymaster
#  ifdef PARALLEL_IO
        inpthread=.true.
        outthread=.true.
#  else
        IF (myrank.eq.0) THEN
          inpthread=.true.
          outthread=.true.
        ELSE
          inpthread=.false.
          outthread=.false.
        END IF
#  endif
!
!  Deallocate local arrays.
!
        IF (allocated(compute_comm))  deallocate (compute_comm)
        IF (allocated(io_comm))       deallocate (io_comm)
        IF (allocated(lranks))        deallocate (lranks)
        IF (allocated(compute_ranks)) deallocate (compute_ranks)
        IF (allocated(io_ranks))      deallocate (io_ranks)
!
!  Otherwise, do synchronous I/O.
!
      ELSE
        DO ng=1,ngrids
          CALL pio_init (myrank,                                        &
     &                   ocn_comm_world,                                &
     &                   pio_numiotasks,                                &
     &                   pio_aggregator,                                &
     &                   pio_stride,                                    &
     &                   pio_rearranger,                                &
     &                   piosystem(ipioroms,ng),                        &
     &                   base = pio_base)
        END DO
      END IF
 
# elif defined ASYNCHRONOUS_SCORPIO
!
!-----------------------------------------------------------------------
!  If SCORPIO library, set the ranks of the computational and I/O
!  IO processes with respect the initial (peer) communicator. Unlike
!  the NCAR/UNIDATA version, there is not a version to split the
!  communicator internally in the call to "PIO_init".  We need to
!  split the communicator previously by calling "set_pio_async".
!  (HGA: this option does not work yet).
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
!!      IF ((inp_lib.eq.io_pio).and.(out_lib.eq.io_pio)) THEN
        IF (out_lib.eq.io_pio) THEN
          CALL pio_init (piosystem(ipioroms:,ng),                       &
     &                   peer_comm_world,                               &
     &                   (/ocn_comm_world/),                            &
     &                   io_comm_world,                                 &
     &                   pio_rearranger)
        ELSE
          CALL pio_init (myrank,                                        &
     &                   ocn_comm_world,                                &
     &                   pio_numiotasks,                                &
     &                   pio_aggregator,                                &
     &                   pio_stride,                                    &
     &                   pio_rearranger,                                &
     &                   piosystem(ipioroms,ng),                        &
     &                   base = pio_base)
        END IF
      END DO
 
# else
!
!-----------------------------------------------------------------------
!  Initialize synchronous PIO system.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        CALL pio_init (myrank,                                          &
     &                 ocn_comm_world,                                  &
     &                 pio_numiotasks,                                  &
     &                 pio_aggregator,                                  &
     &                 pio_stride,                                      &
     &                 pio_rearranger,                                  &
     &                 piosystem(ipioroms,ng),                          &
     &                 base = pio_base)
      END DO
# endif
!
!-----------------------------------------------------------------------
!  Set PIO rearrangement communication options.
!-----------------------------------------------------------------------
!
      lpioinitialized=.true.
!
!  The rearranger communication type "pio_rearr_comm" has two choices:
!
!    PIO_rearr_comm_p2p               Point to point (send/recive)
!    PIO_rearr_comm_coll              Collective (gather/scatter)
!
!  The rearranger communication flow control direction "pio_rearr_fcd"
!  has four choices:
!
!    PIO_rearr_comm_fc_2d_enable      COMM to IO processes and viceversa
!    PIO_rearr_comm_fc_1d_comp2io     COMM to IO processes only
!    PIO_rearr_comm_fc_1d_io2comp     IO to COMM processes only
!    PIO_rearr_comm_fc_2d_disable     Disable flow control
!
!  Compute to IO (C2I) processes:
!
!    pio_rearr_C2I_HS                 Enable handshake (true/false)
!    pio_rearr_C2I_iS                 Enable Isends (true/false)
!    pio_rearr_C2I_PR                 Maximum pending requests
!
!  IO to compute (I2C) processes:
!
!    pio_rearr_I2C_HS                 Enable handshake (true/false)
!    pio_rearr_I2C_iS                 Enable Isends (true/false)
!    pio_rearr_I2C_PR                 Maximum pending requests
!
!  Use PIO_REARR_COMM_UNLIMITED_PEND_REQ for unlimited number of
!  requests.
!
      DO ng=1,ngrids
        myerror=pio_set_rearr_opts(piosystem(ipioroms,ng),              &
                                   pio_rearr_comm,                      &
     &                             pio_rearr_fcd,                       &
     &                             pio_rearr_c2i_hs,                    &
     &                             pio_rearr_c2i_is,                    &
     &                             pio_rearr_c2i_pr,                    &
     &                             pio_rearr_i2c_hs,                    &
     &                             pio_rearr_i2c_is,                    &
     &                             pio_rearr_i2c_pr)
        IF (founderror(myerror, pio_noerr, __line__, myfile)) RETURN
      END DO
!
      RETURN
      END SUBROUTINE initialize_pio
!
      SUBROUTINE finalize_pio
!
!***********************************************************************
!                                                                      !
!  Finalizes the PIO subsystem.  It frees all the storage memory       !
!  associated with the IO decomposition.                               !
!                                                                      !
!***********************************************************************
!
!  Local variable declarations.
!
      integer :: i, ng, status
!
!-----------------------------------------------------------------------
!  Deallocate storage memory associated with IO decomposition.
!-----------------------------------------------------------------------
!
      IF (lpioinitialized) THEN
!
!  Single precision decomposition descriptors.
!
        DO ng=1,ngrids
          DO i=1,npiocomps
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_p2dvar(ng))
# ifdef ADJUST_BOUNDARY
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_r2dobc(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_u2dobc(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_v2dobc(ng))
# endif
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_r2dvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_u2dvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_v2dvar(ng))
# if defined SSH_TIDES || defined UV_TIDES
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_rtides(ng))
# endif
# ifdef SOLVE3D
#  ifdef SEDIMENT
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_b3dvar(ng))
#  endif
#  if defined DIAGNOSTICS_BIO && defined ECOSIM
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_l3dvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_l4dvar(ng))
#  endif
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_p3dvar(ng))
# ifdef ADJUST_BOUNDARY
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_r3dobc(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_u3dobc(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_v3dobc(ng))
# endif
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_r3dvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_u3dvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_v3dvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_w3dvar(ng))
# endif
          END DO
        END DO
!
!  Double precision decomposition descriptors.
!
        DO ng=1,ngrids
          DO i=1,npiocomps
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_p2dvar(ng))
# ifdef ADJUST_BOUNDARY
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_r2dobc(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_u2dobc(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_v2dobc(ng))
# endif
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_r2dvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_u2dvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_v2dvar(ng))
# if defined SSH_TIDES || defined UV_TIDES
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_rtides(ng))
# endif
# ifdef SOLVE3D
#  ifdef SEDIMENT
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_b3dvar(ng))
#  endif
#  if defined DIAGNOSTICS_BIO && defined ECOSIM
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_l3dvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_l4dvar(ng))
#  endif
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_p3dvar(ng))
# ifdef ADJUST_BOUNDARY
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_r3dobc(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_u3dobc(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_v3dobc(ng))
# endif
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_r3dvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_u3dvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_v3dvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_w3dvar(ng))
# endif
          END DO
        END DO
!
!  Special restart and harmonics single precision decomposition
!  descriptors.
!
        DO ng=1,ngrids
          DO i=1,npiocomps
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_rubar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_rvbar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_rzeta(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_ubar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_vbar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_zeta(ng))
# ifdef SOLVE3D
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_ruvel(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_rvvel(ng))
#  if defined GLS_MIXING || defined MY25_MIXING
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_tkevar(ng))
#  endif
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_trcvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_uvel(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_vvel(ng))
# endif
# if defined AVERAGES  && defined AVERAGES_DETIDE && \
    (defined ssh_tides || defined uv_tides)
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_r2dhar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_u2dhar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_v2dhar(ng))
#  ifdef SOLVE3D
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_r3dhar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_u3dhar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_v3dhar(ng))
#  endif
# endif
          END DO
        END DO
!
!  Special restart and harmonics double precision decomposition
!  descriptors.
!
        DO ng=1,ngrids
          DO i=1,npiocomps
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_rubar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_rvbar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_rzeta(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_ubar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_vbar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_zeta(ng))
# ifdef SOLVE3D
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_ruvel(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_rvvel(ng))
#  if defined GLS_MIXING || defined MY25_MIXING
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_tkevar(ng))
#  endif
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_trcvar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_uvel(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_vvel(ng))
# endif
# if defined AVERAGES  && defined AVERAGES_DETIDE && \
    (defined ssh_tides || defined uv_tides)
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_r2dhar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_u2dhar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_v2dhar(ng))
#  ifdef SOLVE3D
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_r3dhar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_u3dhar(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_v3dhar(ng))
#  endif
# endif
# if defined PROPAGATOR && defined CHECKPOINTING
!
!  I/O decomposition descriptors for GST single and double precision
!  state propagator data.
!
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_bvec(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_resid(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_sp_sworkd(ng))
!
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_bvec(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_resid(ng))
            CALL pio_freedecomp (piosystem(i,ng), iodesc_dp_sworkd(ng))
# endif
          END DO
        END DO
!
!-----------------------------------------------------------------------
!  Shut down and clean up any memory associated with the PIO library.
!-----------------------------------------------------------------------
!
        DO ng=1,ngrids
          DO i=1,npiocomps
            CALL pio_finalize (piosystem(i,ng), status)
          END DO
        END DO
      END IF
!
      RETURN
      END SUBROUTINE finalize_pio
!
      SUBROUTINE field_iodecomp (ng, ioSystem, ioType, ioDesc,          &
     &                           gtype, ndims, LBk, UBk, LBt, UBt)
!
!***********************************************************************
!                                                                      !
!  Sets the IO decomposition descriptor for ROMS field variable types. !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng           Nested grid number (integer)                        !
!     ioSystem     PIO system descriptor (TYPE IOSystem_desc_t)        !
!     ioType       PIO kind variable type (integer)                    !
!     gtype        Variable C-grid type (integer)                      !
!     ndims        Number of state variable dimensions (integer)       !
!     LBk          K- or 3rd-dimension Lower bound (integer, OPTIONAL) !
!     UBk          K- or 3rd-dimension Upper bound (integer, OPTIONAL) !
!     LBt          T- or 4th-dimension Lower bound (integer, OPTIONAL) !
!     UBt          T- or 4th-dimension Upper bound (integer, OPTIONAL) !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     ioDesc       IO decomposition descriptor (TYPE io_desc_t)        !
!                                                                      !
!***********************************************************************
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, iotype, gtype, ndims
 
      integer, intent(in), optional :: lbk, ubk
      integer, intent(in), optional :: lbt, ubt
!
      TYPE (iosystem_desc_t), intent(in) :: iosystem
      TYPE (io_desc_t), intent(out) :: iodesc
!
!  Local variable declarations.
!
      logical :: lboundary
!
      integer :: cgrid, ghost
      integer :: i, ic, j, jc, k, kc, l, lc, np
      integer :: is, ie, js, je
      integer :: imin, imax, jmin, jmax
      integer :: ioff, joff, koff, loff
      integer :: ilen, isize, jlen, jsize, klen, ksize, llen, lsize
      integer :: ijlen, ijklen
      integer :: my_size
!
      integer(PIO_Offset_kind), allocatable :: map_decomp(:)
!
!-----------------------------------------------------------------------
!  Set the PIO computational decomposition for ROMS C-type variables
!  and array rank. It is based on variable kind type and its mapping
!  from storage order to memory order.
!-----------------------------------------------------------------------
!
      lboundary=.false.
 
!  Get GLOBAL lower and upper bounds for each variable type in input
!  or ouput NetCDF files.
!
      SELECT CASE (gtype)
        CASE (r2dobc, u2dobc, v2dobc)
          lboundary=.true.
          cgrid=2
          is=0
          ie=iobounds(ng) % IorJ
          js=1
          je=4
          ioff=1
          joff=0
        CASE (r3dobc, u3dobc, v3dobc)
          lboundary=.true.
          cgrid=2
          is=0
          ie=iobounds(ng) % IorJ
          js=1
          je=4
          ioff=1
          joff=0
        CASE (p2dvar, p3dvar)
          cgrid=1
          is=iobounds(ng) % ILB_psi
          ie=iobounds(ng) % IUB_psi
          js=iobounds(ng) % JLB_psi
          je=iobounds(ng) % JUB_psi
          ioff=0
          joff=1
        CASE (r2dvar, b3dvar, l3dvar, l4dvar, r3dvar)
          cgrid=2
          is=iobounds(ng) % ILB_rho
          ie=iobounds(ng) % IUB_rho
          js=iobounds(ng) % JLB_rho
          je=iobounds(ng) % JUB_rho
          ioff=1
          joff=0
        CASE (u2dvar, u3dvar)
          cgrid=3
          is=iobounds(ng) % ILB_u
          ie=iobounds(ng) % IUB_u
          js=iobounds(ng) % JLB_u
          je=iobounds(ng) % JUB_u
          ioff=0
          joff=0
        CASE (v2dvar, v3dvar)
          cgrid=4
          is=iobounds(ng) % ILB_v
          ie=iobounds(ng) % IUB_v
          js=iobounds(ng) % JLB_v
          je=iobounds(ng) % JUB_v
          ioff=1
          joff=1
        CASE (w3dvar)
          cgrid=2
          is=iobounds(ng) % ILB_rho
          ie=iobounds(ng) % IUB_rho
          js=iobounds(ng) % JLB_rho
          je=iobounds(ng) % JUB_rho
          ioff=1
          joff=0
      END SELECT
!
!  Get GLOBAL length for each variable dimension.
!
      ilen=ie-is+1
      jlen=je-js+1
      ijlen=ilen*jlen
!
      IF (PRESENT(lbk)) THEN
        IF (lbk.eq.0) THEN
          koff=0
        ELSE
          koff=1
        END IF
        klen=ubk-lbk+1
        ksize=klen
        ijklen=ijlen*klen
      END IF
!
      IF (PRESENT(lbt)) THEN
        IF (lbt.eq.0) THEN
          loff=0
        ELSE
          loff=1
        END IF
        llen=ubt-lbt+1
        lsize=llen
      END IF
!
!  Starting/ending I- and J-indices for each decomposition tile
!  according to C-grid locatation, excluding ghost points.
!
      IF (lboundary) THEN
        imin=is
        imax=ie
        jmin=js
        jmax=je
      ELSE
        ghost=0
        imin=bounds(ng) % Imin(cgrid,ghost,myrank)
        imax=bounds(ng) % Imax(cgrid,ghost,myrank)
        jmin=bounds(ng) % Jmin(cgrid,ghost,myrank)
        jmax=bounds(ng) % Jmax(cgrid,ghost,myrank)
      END IF
!
!  Allocate 1D array for mapping of the storage order of the variable to
!  its memory order.
!
      isize=imax-imin+1
      jsize=jmax-jmin+1
!
      IF (ndims.eq.2) THEN
        my_size=isize*jsize
      ELSE IF (ndims.eq.3) THEN
        my_size=isize*jsize*ksize
      ELSE IF (ndims.eq.4) THEN
        my_size=isize*jsize*ksize*lsize
      END IF
!
      IF (.not.ALLOCATED(map_decomp)) THEN
        allocate ( map_decomp(my_size) )
      END IF
      map_decomp=0_pio_offset_kind
!
!  Set variable decomposition mapping.
!
      IF (ndims.eq.2) THEN
        np=0
        DO j=jmin,jmax
          jc=(j-joff)*ilen
          DO i=imin,imax
            np=np+1
            ic=i+ioff+jc
            map_decomp(np)=ic
          END DO
        END DO
      ELSE IF (ndims.eq.3) THEN
        np=0
        DO k=lbk,ubk
          kc=(k-koff)*ijlen
          DO j=jmin,jmax
            jc=(j-joff)*ilen+kc
            DO i=imin,imax
              np=np+1
              ic=i+ioff+jc
              map_decomp(np)=ic
            END DO
          END DO
        END DO
      ELSE IF (ndims.eq.4) THEN
        np=0
        DO l=lbt,ubt
          lc=(l-loff)*ijklen
          DO k=lbk,ubk
            kc=(k-koff)*ijlen+lc
            DO j=jmin,jmax
              jc=(j-joff)*ilen+kc
              DO i=imin,imax
                np=np+1
                ic=i+ioff+jc
                map_decomp(np)=ic
              END DO
            END DO
          END DO
        END DO
      END IF
!
!  Set IO decomposition descriptor
!
      IF (ndims.eq.2) THEN
        CALL pio_initdecomp (iosystem, iotype, (/ilen,jlen/),           &
     &                       map_decomp, iodesc)
      ELSE IF (ndims.eq.3) THEN
        CALL pio_initdecomp (iosystem, iotype, (/ilen,jlen,klen/),      &
     &                       map_decomp, iodesc)
      ELSE IF (ndims.eq.4) THEN
        CALL pio_initdecomp (iosystem, iotype, (/ilen,jlen,klen,llen/), &
     &                       map_decomp, iodesc)
      END IF
!
!  Deallocate.
!
      IF (allocated(map_decomp)) deallocate (map_decomp)
!
      RETURN
      END SUBROUTINE field_iodecomp
!
      SUBROUTINE set_iodecomp
!
!***********************************************************************
!                                                                      !
!  Sets the IO decomposition descriptors for ROMS input and output     !
!  variables.  They are used for the mapping between computational     !
!  and I/O processes.                                                  !
!                                                                      !
!***********************************************************************
!
!  Local variable declarations.
!
      integer :: ng
!
!-----------------------------------------------------------------------
!  Allocate I/O decomposition descriptors.
!-----------------------------------------------------------------------
!
!  I/O decomposition descriptors for single precision data.
!
      allocate ( iodesc_sp_p2dvar(ngrids) )
# ifdef ADJUST_WSTRESS
      allocate ( iodesc_sp_u2dfrc(ngrids) )
      allocate ( iodesc_sp_v2dfrc(ngrids) )
# endif
# ifdef ADJUST_BOUNDARY
      allocate ( iodesc_sp_r2dobc(ngrids) )
      allocate ( iodesc_sp_u2dobc(ngrids) )
      allocate ( iodesc_sp_v2dobc(ngrids) )
# endif
      allocate ( iodesc_sp_r2dvar(ngrids) )
      allocate ( iodesc_sp_u2dvar(ngrids) )
      allocate ( iodesc_sp_v2dvar(ngrids) )
# if defined SSH_TIDES || defined UV_TIDES
      allocate ( iodesc_sp_rtides(ngrids) )
# endif
# ifdef SOLVE3D
#  ifdef SEDIMENT
      allocate ( iodesc_sp_b3dvar(ngrids) )
#  endif
#  if defined DIAGNOSTICS_BIO && defined ECOSIM
      allocate ( iodesc_sp_l3dvar(ngrids) )
      allocate ( iodesc_sp_l4dvar(ngrids) )
#  endif
      allocate ( iodesc_sp_p3dvar(ngrids) )
#  ifdef ADJUST_STFLUX
      allocate ( iodesc_sp_r2dfrc(ngrids) )
#  endif
#  ifdef ADJUST_BOUNDARY
      allocate ( iodesc_sp_r3dobc(ngrids) )
      allocate ( iodesc_sp_u3dobc(ngrids) )
      allocate ( iodesc_sp_v3dobc(ngrids) )
#  endif
      allocate ( iodesc_sp_r3dvar(ngrids) )
      allocate ( iodesc_sp_u3dvar(ngrids) )
      allocate ( iodesc_sp_v3dvar(ngrids) )
      allocate ( iodesc_sp_w3dvar(ngrids) )
# endif
!
! I/O decomposition descriptors for double precision data.
!
      allocate ( iodesc_dp_p2dvar(ngrids) )
# ifdef ADJUST_WSTRESS
      allocate ( iodesc_dp_u2dfrc(ngrids) )
      allocate ( iodesc_dp_v2dfrc(ngrids) )
# endif
# ifdef ADJUST_BOUNDARY
      allocate ( iodesc_dp_r2dobc(ngrids) )
      allocate ( iodesc_dp_u2dobc(ngrids) )
      allocate ( iodesc_dp_v2dobc(ngrids) )
# endif
      allocate ( iodesc_dp_r2dvar(ngrids) )
      allocate ( iodesc_dp_u2dvar(ngrids) )
      allocate ( iodesc_dp_v2dvar(ngrids) )
# if defined SSH_TIDES || defined UV_TIDES
      allocate ( iodesc_dp_rtides(ngrids) )
# endif
# ifdef SOLVE3D
#  ifdef SEDIMENT
      allocate ( iodesc_dp_b3dvar(ngrids) )
#  endif
#  if defined DIAGNOSTICS_BIO && defined ECOSIM
      allocate ( iodesc_dp_l3dvar(ngrids) )
      allocate ( iodesc_dp_l4dvar(ngrids) )
#  endif
      allocate ( iodesc_dp_p3dvar(ngrids) )
#  ifdef ADJUST_STFLUX
      allocate ( iodesc_dp_r2dfrc(ngrids) )
#  endif
#  ifdef ADJUST_BOUNDARY
      allocate ( iodesc_dp_r3dobc(ngrids) )
      allocate ( iodesc_dp_u3dobc(ngrids) )
      allocate ( iodesc_dp_v3dobc(ngrids) )
#  endif
      allocate ( iodesc_dp_r3dvar(ngrids) )
      allocate ( iodesc_dp_u3dvar(ngrids) )
      allocate ( iodesc_dp_v3dvar(ngrids) )
      allocate ( iodesc_dp_w3dvar(ngrids) )
# endif
!
!  I/O decomposition descriptors for special single precision
!  restart and harmonics data.
!
      allocate ( iodesc_sp_rubar(ngrids) )
      allocate ( iodesc_sp_rvbar(ngrids) )
      allocate ( iodesc_sp_rzeta(ngrids) )
      allocate ( iodesc_sp_ubar(ngrids) )
      allocate ( iodesc_sp_vbar(ngrids) )
      allocate ( iodesc_sp_zeta(ngrids) )
# ifdef SOLVE3D
      allocate ( iodesc_sp_ruvel(ngrids) )
      allocate ( iodesc_sp_rvvel(ngrids) )
#  if defined GLS_MIXING || defined MY25_MIXING
      allocate ( iodesc_sp_tkevar(ngrids) )
#  endif
      allocate ( iodesc_sp_trcvar(ngrids) )
      allocate ( iodesc_sp_uvel(ngrids) )
      allocate ( iodesc_sp_vvel(ngrids) )
# endif
# if defined AVERAGES  && defined AVERAGES_DETIDE && \
    (defined ssh_tides || defined uv_tides)
!
      allocate ( iodesc_sp_r2dvar(ngrids) )
      allocate ( iodesc_sp_u2dvar(ngrids) )
      allocate ( iodesc_sp_v2dvar(ngrids) )
#  ifdef SOLVE3D
      allocate ( iodesc_sp_r3dvar(ngrids) )
      allocate ( iodesc_sp_u3dvar(ngrids) )
      allocate ( iodesc_sp_v3dvar(ngrids) )
#  endif
# endif
!
!  I/O decomposition descriptors for special double precison
!  restart and harmonics data.
!
      allocate ( iodesc_dp_rubar(ngrids) )
      allocate ( iodesc_dp_rvbar(ngrids) )
      allocate ( iodesc_dp_rzeta(ngrids) )
      allocate ( iodesc_dp_ubar(ngrids) )
      allocate ( iodesc_dp_vbar(ngrids) )
      allocate ( iodesc_dp_zeta(ngrids) )
# ifdef SOLVE3D
      allocate ( iodesc_dp_ruvel(ngrids) )
      allocate ( iodesc_dp_rvvel(ngrids) )
#  if defined GLS_MIXING || defined MY25_MIXING
      allocate ( iodesc_dp_tkevar(ngrids) )
#  endif
      allocate ( iodesc_dp_trcvar(ngrids) )
      allocate ( iodesc_dp_uvel(ngrids) )
      allocate ( iodesc_dp_vvel(ngrids) )
# endif
# if defined AVERAGES  && defined AVERAGES_DETIDE && \
    (defined ssh_tides || defined uv_tides)
!
      allocate ( iodesc_dp_r2dvar(ngrids) )
      allocate ( iodesc_dp_u2dvar(ngrids) )
      allocate ( iodesc_dp_v2dvar(ngrids) )
#  ifdef SOLVE3D
      allocate ( iodesc_dp_r3dvar(ngrids) )
      allocate ( iodesc_dp_u3dvar(ngrids) )
      allocate ( iodesc_dp_v3dvar(ngrids) )
#  endif
# endif
# if defined PROPAGATOR && defined CHECKPOINTING
!
!  I/O decomposition descriptors for GST single and double precision
!  state propagator data. Its values are set in routine "wpoints_tile".
!
      allocate ( iodesc_sp_bvec(ngrids) )
      allocate ( iodesc_sp_resid(ngrids) )
      allocate ( iodesc_sp_sworkd(ngrids) )
!
      allocate ( iodesc_dp_bvec(ngrids) )
      allocate ( iodesc_dp_resid(ngrids) )
      allocate ( iodesc_dp_sworkd(ngrids) )
# endif
!
!-----------------------------------------------------------------------
!  Set the PIO computational decomposition for ROMS C-type variables
!  and array rank. It is based on variable kind type and its mapping
!  from storage order to memory order.
!-----------------------------------------------------------------------
!
!  Set I/O decomposition descriptors for single precision data
!
      DO ng=1,ngrids
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_p2dvar(ng),                      &
     &                       p2dvar, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_r2dvar(ng),                      &
     &                       r2dvar, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_u2dvar(ng),                      &
     &                       u2dvar, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_v2dvar(ng),                      &
     &                       v2dvar, 2)
# if defined SSH_TIDES || defined UV_TIDES
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_rtides(ng),                      &
     &                       r2dvar+4, 3, 1, ntc(ng))
# endif
# ifdef SOLVE3D
#  ifdef SEDIMENT
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_b3dvar(ng),                      &
     &                       b3dvar, 3, 1, nbed)
#  endif
#  if defined DIAGNOSTICS_BIO && defined ECOSIM
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_l3dvar(ng),                      &
     &                       l3dvar, 3, 1, ndbands)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_l4dvar(ng),                      &
     &                       l4dvar, 4, 1, n(ng), ndbands)
#  endif
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_p3dvar(ng),                      &
     &                       p3dvar, 3, 1, n(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_r3dvar(ng),                      &
     &                       r3dvar, 3, 1, n(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_u3dvar(ng),                      &
     &                       u3dvar, 3, 1, n(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_v3dvar(ng),                      &
     &                       v3dvar, 3, 1, n(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_w3dvar(ng),                      &
     &                       w3dvar, 3, 0, n(ng))
# endif
# if defined ADJUST_STFLUX && defined DISTRIBUTE
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_r2dfrc(ng),                      &
     &                       r2dvar, 3, 1, nfrec(ng))
# endif
# ifdef ADJUST_WSTRESS
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_u2dfrc(ng),                      &
     &                       u2dvar, 3, 1, nfrec(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_v2dfrc(ng),                      &
     &                       v2dvar, 3, 1, nfrec(ng))
# endif
# ifdef ADJUST_BOUNDARY
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_r2dobc(ng),                      &
     &                       r2dobc, 3, 1, nbrec(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_u2dobc(ng),                      &
     &                       u2dobc, 3, 1, nbrec(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_v2dobc(ng),                      &
     &                       v2dobc, 3, 1, nbrec(ng))
#  ifdef SOLVE3D
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_r3dobc(ng),                      &
     &                       r3dobc, 4, 1, n(ng), 1, nbrec(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_u3dobc(ng),                      &
     &                       u3dobc, 4, 1, n(ng), 1, nbrec(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_v3dobc(ng),                      &
     &                       v3dobc, 4, 1, n(ng), 1, nbrec(ng))
#  endif
# endif
      END DO
!
!  Set IO decomposition descriptors for double precision data.
!
      DO ng=1,ngrids
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_p2dvar(ng),                      &
     &                       p2dvar, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_r2dvar(ng),                      &
     &                       r2dvar, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_u2dvar(ng),                      &
     &                       u2dvar, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_v2dvar(ng),                      &
     &                       v2dvar, 2)
# if defined SSH_TIDES || defined UV_TIDES
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_rtides(ng),                      &
     &                       r2dvar+4, 3, 1, ntc(ng))
# endif
# ifdef SOLVE3D
#  ifdef SEDIMENT
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_b3dvar(ng),                      &
     &                       b3dvar, 3, 1, nbed)
#  endif
#  if defined DIAGNOSTICS_BIO && defined ECOSIM
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_l3dvar(ng),                      &
     &                       l3dvar, 3, 1, ndbands)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_l4dvar(ng),                      &
     &                       l4dvar, 4, 1, n(ng), ndbands)
#  endif
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_p3dvar(ng),                      &
     &                       p3dvar, 3, 1, n(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_r3dvar(ng),                      &
     &                       r3dvar, 3, 1, n(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_u3dvar(ng),                      &
     &                       u3dvar, 3, 1, n(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_v3dvar(ng),                      &
     &                       v3dvar, 3, 1, n(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_w3dvar(ng),                      &
     &                       w3dvar, 3, 0, n(ng))
# endif
# if defined ADJUST_STFLUX && defined DISTRIBUTE
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_r2dfrc(ng),                      &
     &                       r2dvar, 3, 1, nfrec(ng))
# endif
# ifdef ADJUST_WSTRESS
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_u2dfrc(ng),                      &
     &                       u2dvar, 3, 1, nfrec(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_v2dfrc(ng),                      &
     &                       v2dvar, 3, 1, nfrec(ng))
# endif
# ifdef ADJUST_BOUNDARY
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_r2dobc(ng),                      &
     &                       r2dobc, 3, 1, nbrec(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_u2dobc(ng),                      &
     &                       u2dobc, 3, 1, nbrec(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_v2dobc(ng),                      &
     &                       v2dobc, 3, 1, nbrec(ng))
#  ifdef SOLVE3D
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_r3dobc(ng),                      &
     &                       r3dobc, 4, 1, n(ng), 1, nbrec(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_u3dobc(ng),                      &
     &                       u3dobc, 4, 1, n(ng), 1, nbrec(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_v3dobc(ng),                      &
     &                       v3dobc, 4, 1, n(ng), 1, nbrec(ng))
#  endif
# endif
      END DO
!
!  Set I/O decomposition descriptors for special single precision
!  restart and harmonics data.
!
      DO ng=1,ngrids
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_rubar(ng),                       &
     &                       u2dvar, 3, 1, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_rvbar(ng),                       &
     &                       v2dvar, 3, 1, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_rzeta(ng),                       &
     &                       r2dvar, 3, 1, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_ubar(ng),                        &
     &                       u2dvar, 3, 1, 3)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_vbar(ng),                        &
     &                       v2dvar, 3, 1, 3)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_zeta(ng),                        &
     &                       r2dvar, 3, 1, 3)
# ifdef SOLVE3D
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_ruvel(ng),                       &
     &                       u3dvar, 4, 0, n(ng), 1, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_rvvel(ng),                       &
     &                       v3dvar, 4, 0, n(ng), 1, 2)
#  if defined GLS_MIXING || defined MY25_MIXING
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_tkevar(ng),                      &
     &                       r3dvar, 4, 0, n(ng), 1, 2)
#  endif
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_trcvar(ng),                      &
     &                       r3dvar, 4, 1, n(ng), 1, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_uvel(ng),                        &
     &                       u3dvar, 4, 1, n(ng), 1, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_vvel(ng),                        &
     &                       v3dvar, 4, 1, n(ng), 1, 2)
# endif
 
# if defined AVERAGES  && defined AVERAGES_DETIDE && \
    (defined ssh_tides || defined uv_tides)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_r2dhar(ng),                      &
     &                       r2dvar, 3, 0, 2*ntc(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_u2dhar(ng),                      &
     &                       u2dvar, 3, 0, 2*ntc(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_v2dhar(ng),                      &
     &                       v2dvar, 3, 0, 2*ntc(ng))
#  ifdef SOLVE3D
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_r3dhar(ng),                      &
     &                       r3dvar, 4, 1, n(ng), 0, 2*ntc(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_u3dhar(ng),                      &
     &                       u3dvar, 4, 1, n(ng), 0, 2*ntc(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_real,      &
     &                       iodesc_sp_v3dhar(ng),                      &
     &                       v3dvar, 4, 1, n(ng), 0, 2*ntc(ng))
#  endif
# endif
      END DO
!
!  Set I/O decomposition descriptors for special double precision
!  restart and harmonics data.
!
      DO ng=1,ngrids
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_rubar(ng),                       &
     &                       u2dvar, 3, 1, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_rvbar(ng),                       &
     &                       v2dvar, 3, 1, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_rzeta(ng),                       &
     &                       r2dvar, 3, 1, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_ubar(ng),                        &
     &                       u2dvar, 3, 1, 3)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_vbar(ng),                        &
     &                       v2dvar, 3, 1, 3)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_zeta(ng),                        &
     &                       r2dvar, 3, 1, 3)
# ifdef SOLVE3D
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_ruvel(ng),                       &
     &                       u3dvar, 4, 0, n(ng), 1, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_rvvel(ng),                       &
     &                       v3dvar, 4, 0, n(ng), 1, 2)
#  if defined GLS_MIXING || defined MY25_MIXING
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_tkevar(ng),                      &
     &                       r3dvar, 4, 0, n(ng), 1, 2)
#  endif
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_trcvar(ng),                      &
     &                       r3dvar, 4, 1, n(ng), 1, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_uvel(ng),                        &
     &                       u3dvar, 4, 1, n(ng), 1, 2)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_vvel(ng),                        &
     &                       v3dvar, 4, 1, n(ng), 1, 2)
# endif
 
# if defined AVERAGES  && defined AVERAGES_DETIDE && \
    (defined ssh_tides || defined uv_tides)
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_r2dhar(ng),                      &
     &                       r2dvar, 3, 0, 2*ntc(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_u2dhar(ng),                      &
     &                       u2dvar, 3, 0, 2*ntc(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_v2dhar(ng),                      &
     &                       v2dvar, 3, 0, 2*ntc(ng))
#  ifdef SOLVE3D
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_r3dhar(ng),                      &
     &                       r3dvar, 4, 1, n(ng), 0, 2*ntc(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_u3dhar(ng),                      &
     &                       u3dvar, 4, 1, n(ng), 0, 2*ntc(ng))
        CALL field_iodecomp (ng, piosystem(ipioroms,ng), pio_double,    &
     &                       iodesc_dp_v3dhar(ng),                      &
     &                       v3dvar, 4, 1, n(ng), 0, 2*ntc(ng))
#  endif
# endif
      END DO
!
      RETURN
      END SUBROUTINE set_iodecomp
 
# ifdef ASYNCHRONOUS_SCORPIO
!
      SUBROUTINE set_pio_async
!
!***********************************************************************
!                                                                      !
!  If SCORPIO library, splits the distributed-memory communicator to   !
!  allow asynchronous I/O with dedicated processes.                    !
!                                                                      !
!***********************************************************************
!
!  Local variable declarations.
!
      logical :: lsplit
      logical, allocatable :: lranks(:)
!
      integer :: computesize, lstr, myerror, serror
      integer :: iomaster, ocnmaster
      integer :: i, ic, last, ng, npets, tag
      integer :: grp_compute, grp_initial, grp_io
      integer :: petrange(3,1)                        ! one triplet only
!
      integer, allocatable :: compute_ranks(:,:), io_ranks(:)
!
      character (len=MPI_MAX_ERROR_STRING) :: string
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", set_pio_async"
!
!-----------------------------------------------------------------------
!  Split distributed-memory communicator into computational and
!  dedicated I/O processes.
!-----------------------------------------------------------------------
!
!  Set peer communicator as the initial ROMS communicator.  After
!  splitting, the OCN_COMM_WORLD must have the value of MPI_COMM_NULL
!  on those processes dedicated to I/O.
!
      peer_comm_world=ocn_comm_world
      ocn_comm_world=mpi_comm_null
      io_comm_world=mpi_comm_null
!
      peersize=numthreads
      peerrank=myrank
      computesize=peersize-pio_numiotasks
!
!  Allocate and initialize local arrays.
!
      IF (.not.allocated(lranks)) THEN
        allocate ( lranks(0:peersize-1) )
        lranks(0:peersize-1)=.true.
      END IF
      IF (.not.allocated(compute_ranks)) THEN
        allocate ( compute_ranks(computesize,npiocomps) )
        compute_ranks=-1
      END IF
      IF (.not.allocated(io_ranks)) THEN
        allocate ( io_ranks(pio_numiotasks) )
        io_ranks=-1
      END IF
!
!  Set the ranks of the dedicated I/O processes with respect the peer
!  communicator.
!
      ic=pio_base
      io_ranks(1)=ic
      lranks(ic)=.false.
      DO i=2,pio_numiotasks
        IF ((ic+pio_stride).le.(peersize-1)) THEN
          ic=ic+pio_stride
          io_ranks(i)=ic
          lranks(ic)=.false.
        END IF
      END DO
      iomaster=io_ranks(1)
      WRITE (cioranks,'(*(i0,1x))') io_ranks
!
!  Set the ranks of the computational processes with respect the peer
!  communicator.
!
      ic=0
      DO i=0,peersize-1
        IF (lranks(i)) THEN
          ic=ic+1
          compute_ranks(ic,ipioroms)=i
        END IF
      END DO
      ocnmaster=compute_ranks(1,ipioroms)
      WRITE (ccompranks,'(*(i0,1x))') compute_ranks(:,ipioroms)
!
!  Set switch to split communicator.
!
      IF ((inp_lib.eq.io_pio).or.(out_lib.eq.io_pio)) THEN
        lsplit=.true.
      ELSE
        lsplit=.false.
      END IF
!
!  Check for correct parallel decompositon between computational and
!  I/O processes or when split communicator is unnecessary.
!
      DO ng=1,ngrids
        IF (lsplit) THEN
          npets=pio_numiotasks+ntilei(ng)*ntilej(ng)
          IF (npets.ne.peersize) THEN
            IF (peerrank.eq.0)                                          &
              print 10, ' asynchonous ', ng,                            &
     &                  'NumIOtasks+NtileI(ng)*NtileJ(ng) = ', npets,   &
     &                  'Peer Communicator Size           = ', peersize,&
     &                  'NumIOtasks+NtileI(ng)*NtileJ(ng) = ', peersize
            exit_flag=6
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END IF
        ELSE
          npets=ntilei(ng)*ntilej(ng)
          IF (npets.ne.peersize) THEN
            IF (peerrank.eq.0)                                          &
              print 10, ' synchonous ', ng,                             &
     &                  'NtileI(ng)*NtileJ(ng) = ', npets,              &
     &                  'Communicator Size     = ', peersize,           &
     &                  'NtileI(ng)*NtileJ(ng) = ', npets
            exit_flag=6
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END IF
        END IF
      END DO
!
!-----------------------------------------------------------------------
!  Split initial communicator into disjointed computational and I/O
!  subgroups.
!-----------------------------------------------------------------------
!
!  Set dedicated range I/O processes triplets. All ranks must be valid
!  in group and all computed ranks must be distinct.
!  Hmm, the grouping below works only for base=0. That is weid!
!
      petrange(1,1)=pio_base                              ! first rank
      petrange(2,1)=io_ranks(pio_numiotasks)              ! last rank
      petrange(3,1)=pio_stride                            ! rank stride
!
!  Split initial communicator into disjointed subgroups.
!
      IF (lsplit) THEN
!
!  Create a group associated with the initial communicator.
!
        CALL mpi_comm_group (peer_comm_world, grp_initial, myerror)
        IF (myerror.ne.mpi_success) THEN
          CALL mpi_error_string (myerror, string, lstr, serror)
          IF (peerrank.eq.0)                                            &
             print 20, 'MPI_COMM_GROUP',                                &
     &                 'Comm = PEER_COMM_WORLD',                        &
     &                 peerrank, myerror, trim(string)
          exit_flag=2
        END IF
!
!  Set group of processes the I/O group. Include the processes range
!  defined in PetRange.
!
        CALL mpi_group_range_incl (grp_initial, 1, petrange,            &
     &                             grp_io, myerror)
        IF (myerror.ne.mpi_success) THEN
          CALL mpi_error_string (myerror, string, lstr, serror)
          IF (peerrank.eq.iomaster)                                     &
            print 20, 'MPI_GROUP_RANGE_INCL',                           &
     &                'Group = GRP_IO',                                 &
     &                iomaster, myerror, trim(string)
          exit_flag=2
        END IF
!
!  Set group of process in the computational group. Exclude the
!  precesses range defined in PetRange.
!
        CALL mpi_group_range_excl (grp_initial, 1, petrange,            &
     &                             grp_compute, myerror)
        IF (myerror.ne.mpi_success) THEN
          CALL mpi_error_string (myerror, string, lstr, serror)
          IF (peerrank.eq.ocnmaster)                                    &
     &      print 20, 'MPI_GROUP_RANGE_EXCL',                           &
     &                'Group = GRP_COMPUTE',                            &
     &                ocnmaster, myerror, trim(string)
          exit_flag=2
        END IF
!
!  Create the computational group associated with the initial
!  communicator.
!
        CALL mpi_comm_create (peer_comm_world, grp_compute,             &
     &                        ocn_comm_world, myerror)
        IF (myerror.ne.mpi_success) THEN
          CALL mpi_error_string (myerror, string, lstr, serror)
          IF (peerrank.eq.ocnmaster)                                    &
     &      print 20, 'MPI_COMM_CREATE',                                &
     &                'Comm = OCN_COMM_WORLD',                          &
     &                ocnmaster, myerror, trim(string)
          exit_flag=2
        END IF
!
!  Create dedicated I/O group associated the initial communicator.
!  Notice that at least
!
        CALL mpi_comm_create (peer_comm_world, grp_io,                  &
     &                        io_comm_world, myerror)
        IF (myerror.ne.mpi_success) THEN
          CALL mpi_error_string (myerror, string, lstr, serror)
          IF (peerrank.eq.iomaster)                                     &
     &      print 20, 'MPI_COMM_CREATE',                                &
     &                'Comm = IO_COMM_WORLD',                           &
     &                iomaster, myerror, trim(string)
          exit_flag=2
        END IF
!
!  Create inter-communicator from two existing intra-communicators.
!  Notice that at least one selected process from each group has the
!  ability to communicate with the selected member of the other group.
!
        tag=1
        IF (ocn_comm_world.ne.mpi_comm_null) THEN
          CALL mpi_intercomm_create (ocn_comm_world, 0, peer_comm_world,&
     &                               pio_base, tag, inter_comm_world,   &
     &                               myerror)
          IF (myerror.ne.mpi_success) THEN
            CALL mpi_error_string (myerror, string, lstr, serror)
            IF (peerrank.eq.ocnmaster)                                  &
              print 20, 'MPI_ITERCOMM_CREATE',                          &
     &                  'Comm = INTER_COMM_WORLD',                      &
     &                  ocnmaster, myerror, trim(string)
            exit_flag=2
          END IF
        ELSE IF (io_comm_world.ne.mpi_comm_null) THEN
          IF (pio_base.eq.0) THEN
            IF (pio_stride.gt.1) THEN
              CALL mpi_intercomm_create (io_comm_world, 0,              &
     &                                   peer_comm_world, 1, tag,       &
     &                                   inter_comm_world, myerror)
              IF (myerror.ne.mpi_success) THEN
                CALL mpi_error_string (myerror, string, lstr, serror)
                IF (peerrank.eq.iomaster)                               &
                  print 20, 'MPI_ITERCOMM_CREATE',                      &
     &                      'Comm = INTER_COMM_WORLD',                  &
     &                      iomaster, myerror, trim(string)
                exit_flag=2
              END IF
            ELSE
              CALL mpi_intercomm_create (io_comm_world, 0,              &
     &                                   peer_comm_world,               &
     &                                   pio_numiotasks, tag,           &
     &                                   inter_comm_world, myerror)
              IF (myerror.ne.mpi_success) THEN
                CALL mpi_error_string (myerror, string, lstr, serror)
                IF (peerrank.eq.iomaster)                               &
                  print 20, 'MPI_ITERCOMM_CREATE',                      &
     &                      'Comm = INTER_COMM_WORLD',                  &
     &                      iomaster, myerror, trim(string)
                exit_flag=2
              END IF
            END IF
          ELSE
            CALL mpi_intercomm_create (io_comm_world, 0,                &
     &                                 peer_comm_world, 0, tag,         &
     &                                 inter_comm_world, myerror)
            IF (myerror.ne.mpi_success) THEN
              CALL mpi_error_string (myerror, string, lstr, serror)
              IF (peerrank.eq.iomaster)                                 &
     &          print 20, 'MPI_ITERCOMM_CREATE',                        &
     &                    'Comm = INTER_COMM_WORLD',                    &
     &                    iomaster, myerror, trim(string)
              exit_flag=2
            END IF
          END IF
        ELSE
          exit_flag=2
          IF (peerrank.eq.0) print 30
        END IF
      ELSE
        ocn_comm_world=peer_comm_world                    ! not splitted
      END IF
!
!-----------------------------------------------------------------------
!  Initialize parallel control switches.
!-----------------------------------------------------------------------
!
      IF (lsplit) THEN
        CALL mpi_comm_size (ocn_comm_world, numthreads, myerror)
        CALL mpi_comm_rank (ocn_comm_world, myrank, myerror)
!
        master=myrank.eq.mymaster
 
# ifdef PARALLEL_IO
        inpthread=.true.
        outthread=.true.
# else
        IF (myrank.eq.0) THEN
          inpthread=.true.
          outthread=.true.
        ELSE
          inpthread=.false.
          outthread=.false.
        END IF
# endif
      END IF
!
!  Deallocate local arrays.
!
      IF (allocated(lranks))        deallocate (lranks)
      IF (allocated(compute_ranks)) deallocate (compute_ranks)
      IF (allocated(io_ranks))      deallocate (io_ranks)
!
 10   FORMAT (/,' SET_PIO_ASYNC - Wrong number of processors for',      &
     &        a,'I/O in grid: ',i0,/,17x,a,i0,/,17x,a,i0,               &
     &        ', but we need: ',/,17x,a,i0,/)
 20   FORMAT (/,' SET_PIO_ASYNC - error during ',a,                     &
     &        ' call, ',a,', Rank = ',i0,' Error = ',i0,/,22x,a)
 30   FORMAT (/,' SET_PIO_ASYNC - no option to create ',                &
     &        'inter-communicator.')
!
      RETURN
      END SUBROUTINE set_pio_async
# endif
# if defined PROPAGATOR && defined CHECKPOINTING
!
      SUBROUTINE state_iodecomp (ng, ioSystem, ioType, ioDesc,          &
     &                           ioVname, ndims)
!
!***********************************************************************
!                                                                      !
!  Sets the IO decomposition descriptor for ROMS packed state variable !
!  types.                                                              !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng           Nested grid number (integer)                        !
!     ioSystem     PIO system descriptor (TYPE IOSystem_desc_t)        !
!     ioType       PIO kind variable type (integer)                    !
!     ioVname      State variable name (string)                        !
!     ndims        Number of state variable dimensions (integer)       !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     ioDesc       IO decomposition descriptor (TYPE io_desc_t)        !
!                                                                      !
!***********************************************************************
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, iotype, ndims
!
      character (len=*) :: iovname
!
      TYPE (iosystem_desc_t), intent(in) :: iosystem
      TYPE (io_desc_t), intent(out) :: iodesc
!
!  Local variable declarations.
!
      integer :: is, ie, isize, js, je, jsize, my_size
      integer :: imin, imax, ioff, jmin, jmax, joff
      integer :: i, ic, j, jc, np
 
      integer(PIO_Offset_kind), allocatable :: map_decomp(:)
!
!-----------------------------------------------------------------------
!  Set the PIO computational decomposition for ROMS packed state
!  variables.
!-----------------------------------------------------------------------
!
      SELECT CASE (trim(iovname))
        CASE ('Bvec')
          is=nstr(ng)
          ie=nend(ng)
          js=1
          je=ncv
          ioff=0
          joff=1
        CASE ('resid')
          is=nstr(ng)
          ie=nend(ng)
        CASE ('SworkD')
          is=myrank*3*nstate(ng)+1
          ie=min(is+3*nstate(ng)-1, 3*mstate(ng))
      END SELECT
!
!  Starting/ending I- and J-indices for each decomposition tile.
!
      IF (ndims.eq.1) THEN
        imin=is
        imax=ie
        isize=ie-is+1
        my_size=isize
      ELSE IF (ndims.eq.2) THEN
        imin=is
        imax=ie
        jmin=js
        jmax=je
        isize=ie-is+1
        jsize=je-js+1
        my_size=isize*jsize
      END IF
!
!  Allocate 1D array for mapping of the storage order of the variable to
!  its memory order.
!
      IF (.not.ALLOCATED(map_decomp)) THEN
        allocate ( map_decomp(my_size) )
      END IF
      map_decomp=0_pio_offset_kind
!
!  Set variable decomposition mapping.
!
      IF (ndims.eq.1) THEN
        np=0
        DO i=imin,imax
          np=np+1
          map_decomp(np)=i
        END DO
      ELSE IF (ndims.eq.2) THEN
        np=0
        DO j=jmin,jmax
          jc=(j-joff)*isize
          DO i=imin,imax
            np=np+1
            ic=i+ioff+jc
            map_decomp(np)=ic
          END DO
        END DO
      END IF
!
!  Set IO decomposition descriptor
!
      IF (ndims.eq.1) THEN
        CALL pio_initdecomp (iosystem, iotype, (/isize/),               &
     &                       map_decomp, iodesc)
      ELSE IF (ndims.eq.2) THEN
        CALL pio_initdecomp (iosystem, iotype, (/isize,jsize/),         &
     &                       map_decomp, iodesc)
      END IF
!
!  Deallocate.
!
      IF (ALLOCATED(map_decomp)) deallocate (map_decomp)
!
      RETURN
      END SUBROUTINE state_iodecomp
# endif
#endif
      END MODULE set_pio_mod