!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

 module io_binary 2,13

!BOP
! !MODULE: io_binary
! !DESCRIPTION:
!  This module provides a binary parallel input/output interface
!  for writing fields.
!
! !REVISION HISTORY:
!  SVN:$Id: io_binary.F90 19531 2009-11-19 19:05:44Z njn01 $

! !USES:

   use POP_CommMod
   use kinds_mod
   use domain_size
   use domain
   use constants
   use communicate
   use broadcast
   use gather_scatter
   use exit_mod
   use io_types

#ifdef USEPIO
   use piolib_mod, only:  &                             ! _EXTERNAL
        PIO_OpenFile, PIO_CloseFile, PIO_CreateFile, &
        PIO_write_darray, PIO_read_darray, &
        PIO_numToRead, PIO_numToWrite
   use pio_types, only: IO_desc_t, Var_desc_t, iotype_direct_pbinary,iosystem_desc_t          ! _EXTERNAL
   use pio
#endif

   implicit none
   private
   save

! !PUBLIC MEMBER FUNCTIONS:

   public :: open_read_binary,    &
             open_binary,         &
             close_binary,        &
             define_field_binary, &
             read_field_binary,   &
             write_field_binary

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  module variables
!
!-----------------------------------------------------------------------

   character (1), parameter :: &
      attrib_separator=':'  ! character separating attrib name and value

!-----------------------------------------------------------------------
!
!  binary interfaces for i/o routines
!
!-----------------------------------------------------------------------


   interface write_array 12
      module procedure write_int_2d
                       write_real4_2d, &
                       write_real8_2d, &
                       write_int_3d,   &
                       write_real4_3d, &
                       write_real8_3d
   end interface


   interface read_array 12
     module procedure read_int_2d
                      read_real4_2d, &
                      read_real8_2d, &
                      read_int_3d,   &
                      read_real4_3d, &
                      read_real8_3d
   end interface

!EOC
!***********************************************************************

 contains

!***********************************************************************
!BOP
! !IROUTINE: open_read_binary
! !INTERFACE:


 subroutine open_read_binary(data_file) 1,10

! !DESCRIPTION:
!  This routine opens a binary data file (and header file if it
!  exists) for reading.  It also reads global file attributes.
!
! !REVISION HISTORY:
!  same as module

! !INPUT/OUTPUT PARAMETERS:

   type (datafile), intent (inout)  :: &
      data_file       ! file to be opened

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   character (char_len) :: &
      path,          &! filename to read
      work_line,     &! temporary to use for parsing file lines
      comp_line,     &! temporary to use for parsing file lines
      data_type,     &! temporary to use for data type
      att_name        ! attribute name

   integer (i4) ::  &
      n,            &! loop index
      att_ival,     &! value for integer attribute
      hdr_error,    &! error for reading header file
      cindx1,cindx2  ! indices for unpacking character strings

   logical (log_kind) ::     &
      att_lval       ! value for logical attribute

   real (r4) ::     &
      att_rval       ! value for real attribute

   real (r8) ::     &
      att_dval       ! value for double attribute

   logical (log_kind) :: &
      header_exists,     &! flag to check existence of header file
      lmatch              ! flag to use for attribute search

   integer (i4) :: iostat,ierr
!-----------------------------------------------------------------------
!
!  set the readonly flag in the data file descriptor
!
!-----------------------------------------------------------------------

   data_file%readonly = .true.

!-----------------------------------------------------------------------
!
!  if unit not assigned already, assign a unit (id) to this file
!
!-----------------------------------------------------------------------

   if (data_file%id(1) == 0) then
     call get_unit(data_file%id(1))
     call get_unit(data_file%id(2))
   endif

!-----------------------------------------------------------------------
!
!  check to see whether a header file exists
!    if so, open it
!    if not, set the unit to -unit so it can be released later
!
!-----------------------------------------------------------------------

   header_exists = .false.

   if (my_task == master_task) then
      path = trim(data_file%full_name)/&
             &/'.hdr'

      inquire(file=path, exist=header_exists)

      if (header_exists) then
         open (unit=data_file%id(2), file=path, status='old')
      else
         write(stdout,*) 'WARNING: Input header file does not exist'
         write(stdout,*) 'for file: ',trim(path)
         write(stdout,*) &
                  'Assuming fields will be read in sequential order'
      endif
   endif

   call broadcast_scalar(header_exists, master_task)

   if (.not. header_exists) then
      data_file%id(2) = -data_file%id(2)
      data_file%current_record = 1
   endif

!-----------------------------------------------------------------------
!
!  now open data file
!
!-----------------------------------------------------------------------

#ifdef USEPIOREAD
   iostat = PIO_OpenFile(PIO_desc,data_file%file,iotype_direct_pbinary, trim(data_file%full_name))  
#else
   if (my_task < data_file%num_iotasks) then
      path = trim(data_file%full_name)
      open (unit=data_file%id(1),action='read',status='unknown', &
            file=trim(path), form='unformatted',access='direct', &
            recl=data_file%record_length)
   endif
#endif

!-----------------------------------------------------------------------
!
!  if the header file exists read file attributes from header file
!
!-----------------------------------------------------------------------

   if (header_exists) then

      work_line = char_blank
      hdr_error = 0

      if (my_task == master_task) then
         !*** first read until hit beginning of global attribute
         !*** definition section

         do while (hdr_error == 0 .and. work_line(1:7) /= '&GLOBAL')
            read(data_file%id(2),'(a80)',iostat=hdr_error) comp_line
            work_line = adjustl(comp_line)
         end do
      end if

      !*** now read until hit end of global attribute definition

      att_loop: do while (hdr_error == 0)

         work_line = char_blank
         if (my_task == master_task) then
            read(data_file%id(2),'(a80)',iostat=hdr_error) comp_line
            work_line = adjustl(comp_line)
         endif

         call broadcast_scalar(work_line, master_task)
         call broadcast_scalar(hdr_error, master_task)

         if (work_line(1:1) == '/') exit att_loop

         !*** find location of separator in the string
         !*** to determine name of attribute and extract
         !*** attribute name from beginning of string

         cindx1 = index(work_line,attrib_separator)
         att_name = char_blank
         att_name(1:cindx1-1) = work_line(1:cindx1-1)

         comp_line = work_line
         do n=1,cindx1
            comp_line(n:n) = ' '
         end do
         work_line = adjustl(comp_line)

         !*** now find location of separator in the string
         !*** to determine data type of attribute and extract
         !*** data type from beginning of string

         cindx1 = index(work_line,attrib_separator)
         data_type = char_blank
         data_type(1:cindx1-1) = work_line(1:cindx1-1)

         comp_line = work_line
         do n=1,cindx1
            comp_line(n:n) = ' '
         end do
         work_line = adjustl(comp_line)

         !*** check for standard file attributes

         select case(trim(att_name))
         case ('title','TITLE')
            data_file%title = trim(work_line)
         case ('history','HISTORY')
            data_file%history = trim(work_line)
         case ('conventions','CONVENTIONS')
            data_file%conventions = trim(work_line)
         case default

            !*** if additional attributes exist, add them as
            !*** additional attributes
            !*** note that if they are already defined, the
            !*** call to add_attrib will overwrite value

            select case (trim(data_type))
            case('char','CHAR')
               call add_attrib_file(data_file, trim(att_name), &
                                               trim(work_line))
            case('log','LOG','logical','LOGICAL')
               read(work_line,*) att_lval 
               call add_attrib_file(data_file, trim(att_name), att_lval)
            case('int','INT','i4','I4')
               read(work_line,*) att_ival 
               call add_attrib_file(data_file, trim(att_name), att_ival)
            case('r4','R4','REAL','real','float','FLOAT')
               read(work_line,*) att_rval 
               call add_attrib_file(data_file, trim(att_name), att_rval)
            case('r8','R8','dbl','DBL','double','DOUBLE')
               read(work_line,*) att_dval 
               call add_attrib_file(data_file, trim(att_name), att_dval)
            end select

         end select
      end do att_loop

   endif ! header exists

!-----------------------------------------------------------------------
!EOC

 end subroutine open_read_binary

!***********************************************************************
!BOP
! !IROUTINE: open_binary
! !INTERFACE:


 subroutine open_binary(data_file) 1,2

! !DESCRIPTION:
!  This routine opens a binary data file and header file for writing
!  and writes global file attributes to the header file.
!
! !REVISION HISTORY:
!  same as module

! !INPUT/OUTPUT PARAMETERS:

   type (datafile), intent (inout)  :: &
      data_file              ! file to open

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   character (255) :: path

   character (char_len_long) :: &
      work_line  ! temp space for manipulating strings 

   character (5), parameter :: &
      hdr_fmt = '(a80)'

   integer (i4) ::    &
      n,              &! loop index
      cindx1, cindx2, &! indices for character strings
      hdr_unit         ! unit number for header file
    
   integer(i4) :: ierr

!-----------------------------------------------------------------------
!
!  assign units to file if not assigned already
!
!-----------------------------------------------------------------------

   if (data_file%id(1) == 0) then
      call get_unit(data_file%id(1))
      call get_unit(data_file%id(2))
   endif

!-----------------------------------------------------------------------
!
!  open data and header files for writing
!
!-----------------------------------------------------------------------

   !*** open header file from one task only

   if (my_task == master_task) then
      path = trim(data_file%full_name)/&
                                       &/'.hdr'
      open (unit=data_file%id(2), file=path, status='unknown')
      hdr_unit = data_file%id(2)
   end if

   path = trim(data_file%full_name)

#ifdef USEPIO
   ierr = PIO_CreateFile(PIO_desc,data_file%file,iotype_direct_pbinary, path)
#else
   !*** open data file from all io tasks

   if (my_task < data_file%num_iotasks) then
      open(data_file%id(1), file=path, access='direct',      &
           form='unformatted', recl=data_file%record_length, &
           status='unknown')
   endif
#endif

   !*** initialize record number

   data_file%current_record = 1

!-----------------------------------------------------------------------
!
!  write attributes to header file
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then

      work_line = char_blank
      work_line(1:7) = '&GLOBAL'
      write(hdr_unit,hdr_fmt) work_line

      work_line = char_blank
      work_line(1:5) = 'title'
      work_line(6:6) = attrib_separator
      work_line(7:10) = 'char'
      work_line(11:11) = attrib_separator
      write(work_line(12:),*) trim(data_file%title)
      write(hdr_unit,hdr_fmt) work_line

      work_line = char_blank
      work_line(1:7) = 'history'
      work_line(8:8) = attrib_separator
      work_line(9:12) = 'char'
      work_line(13:13) = attrib_separator
      write(work_line(14:),*) trim(data_file%history)
      write(hdr_unit,hdr_fmt) work_line

      work_line = char_blank
      work_line(1:11) = 'conventions'
      work_line(12:12) = attrib_separator
      work_line(11:14) = 'char'
      work_line(15:15) = attrib_separator
      write(work_line(16:),*) trim(data_file%conventions)
      write(hdr_unit,hdr_fmt) work_line

      !*** if additional attributes are defined in the
      !*** file definition, write these as well

      if (associated(data_file%add_attrib_cval)) then
         do n=1,size(data_file%add_attrib_cval)

            work_line = char_blank
            cindx1 = 1
            cindx2 = len_trim(data_file%add_attrib_cname(n))
            work_line(cindx1:cindx2) = trim(data_file%add_attrib_cname(n))
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            cindx2 = cindx1 + 3
            work_line(cindx1:cindx2) = 'char'
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            write(work_line(cindx1:),*) trim(data_file%add_attrib_cval(n))
            write(hdr_unit,hdr_fmt) work_line

         end do
      endif ! cval

      if (associated(data_file%add_attrib_lval)) then
         do n=1,size(data_file%add_attrib_lval)

            work_line = char_blank
            cindx1 = 1
            cindx2 = len_trim(data_file%add_attrib_lname(n))
            work_line(cindx1:cindx2)=trim(data_file%add_attrib_lname(n))
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            cindx2 = cindx1 + 2
            work_line(cindx1:cindx2) = 'log'
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            write(work_line(cindx1:),*) data_file%add_attrib_lval(n)
            write(hdr_unit,hdr_fmt) work_line

         end do
      endif ! lval

      if (associated(data_file%add_attrib_ival)) then
         do n=1,size(data_file%add_attrib_ival)

            work_line = char_blank
            cindx1 = 1
            cindx2 = len_trim(data_file%add_attrib_iname(n))
            work_line(cindx1:cindx2)=trim(data_file%add_attrib_iname(n))
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            cindx2 = cindx1 + 2
            work_line(cindx1:cindx2) = 'int'
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            write(work_line(cindx1:),*) data_file%add_attrib_ival(n)
            write(hdr_unit,hdr_fmt) work_line

         end do
      endif ! ival

      if (associated(data_file%add_attrib_rval)) then
         do n=1,size(data_file%add_attrib_rval)

            work_line = char_blank
            cindx1 = 1
            cindx2 = len_trim(data_file%add_attrib_rname(n))
            work_line(cindx1:cindx2) = &
                                trim(data_file%add_attrib_rname(n))
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            cindx2 = cindx1 + 1
            work_line(cindx1:cindx2) = 'r4'
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            write(work_line(cindx1:),*) data_file%add_attrib_rval(n)
            write(hdr_unit,hdr_fmt) work_line

         end do
      endif ! rval

      if (associated(data_file%add_attrib_dval)) then
         do n=1,size(data_file%add_attrib_dval)

            work_line = char_blank
            cindx1 = 1
            cindx2 = len_trim(data_file%add_attrib_dname(n))
            work_line(cindx1:cindx2)=trim(data_file%add_attrib_dname(n))
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            cindx2 = cindx1 + 1
            work_line(cindx1:cindx2) = 'r8'
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            write(work_line(cindx1:),*) data_file%add_attrib_dval(n)
            write(hdr_unit,hdr_fmt) work_line

         end do
      endif ! dval

   endif ! master task

!-----------------------------------------------------------------------
!EOC

 end subroutine open_binary

!***********************************************************************
!BOP
! !IROUTINE: close_binary
! !INTERFACE:


 subroutine close_binary(data_file) 1

! !DESCRIPTION:
!  This routine closes an open binary data file and the associated
!  header file if it exists.
!
! !REVISION HISTORY:
!  same as module

! !INPUT/OUTPUT PARAMETERS:

   type (datafile), intent (inout)  :: data_file

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  close a data file
!
!-----------------------------------------------------------------------
#ifdef USEPIO
   if(data_file%readonly) then
#ifdef USEPIOREAD
      call PIO_CloseFile(data_file%file)
#else
      if (my_task < data_file%num_iotasks) close(data_file%id(1))
#endif
   else
      call PIO_CloseFile(data_file%file)
   endif
#else
   if (my_task < data_file%num_iotasks) close(data_file%id(1))
#endif

   if (my_task == master_task .and. data_file%id(2) > 0) &
      close(data_file%id(2))

!-----------------------------------------------------------------------
!EOC

 end subroutine close_binary

!***********************************************************************
!BOP
! !IROUTINE: define_field_binary
! !INTERFACE:


 subroutine define_field_binary(data_file, io_field) 1,12

! !DESCRIPTION:
!  This routine defines an io field.  When writing a file, this means
!  computing the starting record number for the field and then writing
!  the field information to the header file.  When reading a file,
!  the define routine will attempt to fill an io field with information
!  from the header file; if a header file does not exist, it will leave
!  the io field mostly undefined and will rely on an input record from
!  the read routine to locate the field in the file (this latter is
!  provided for compatibility with forcing and grid files).
!
! !REVISION HISTORY:
!  same as module

! !INPUT/OUTPUT PARAMETERS:

   type (datafile), intent (inout)  :: &
      data_file

   type (io_field_desc), intent (inout) :: &
      io_field

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   character (80) :: &
      work_line,     &! workspace for manipulating input string
      comp_line,     &! comparison string
      att_name,      &! attribute name string
      ctype           ! attribute data type

   integer (i4) ::  &
      n,            &! loop index
      att_ival,     &! temp value for integer attribute
      hdr_error,    &! io error status for header file
      cindx1,       &! character string index
      cindx2,       &! character string index
      unit           ! unit for header file

   logical (log_kind) :: &
      att_lval   ! temp value for logical attribute

   real (r4) :: &
      att_rval   ! temp value for real attribute

   real (r8) :: &
      att_dval   ! temp value for double attribute

!-----------------------------------------------------------------------
!
!  for output files, set id to current record number and write
!  field attributes to header file.
!
!-----------------------------------------------------------------------

   if (.not. data_file%readonly) then

      !*** define id as current record

      io_field%id = data_file%current_record

      !*** increment record counter for next field

      if (io_field%nfield_dims == 2) then

         data_file%current_record = data_file%current_record + 1

      else if (io_field%nfield_dims == 3) then

         if (associated(io_field%field_i_3d)) then
            n = size(io_field%field_i_3d,dim=3)
         else if (associated(io_field%field_r_3d)) then
            n = size(io_field%field_r_3d,dim=3)
         else if (associated(io_field%field_d_3d)) then
            n = size(io_field%field_d_3d,dim=3)
         endif

         data_file%current_record = data_file%current_record + n

      endif

      !***
      !*** now write attributes to header file
      !***

      unit = data_file%id(2)

      if (my_task == master_task) then

         !*** first write variable name header

         write(unit,'(a1,a)') '&',trim(io_field%short_name)

         !*** write standard attributes to header file

         if (io_field%long_name /= char_blank)                           &
            write(unit,'(a9,a1,a4,a1,a)') 'long_name', attrib_separator, &
                                          'char',      attrib_separator, &
                                          trim(io_field%long_name)

         if (io_field%units /= char_blank)                           &
            write(unit,'(a5,a1,a4,a1,a)') 'units', attrib_separator, &
                                          'char' , attrib_separator, &
                                          trim(io_field%units)

         if (io_field%coordinates /= char_blank)                           &
            write(unit,'(a5,a1,a4,a1,a)') 'coordinates', attrib_separator, &
                                          'char' , attrib_separator,       &
                                          trim(io_field%coordinates)

         if (io_field%grid_loc /= '    ')                                &
            write(unit,'(a8,a1,a4,a1,a4)') 'grid_loc', attrib_separator, &
                                           'char'    , attrib_separator, &
                                           io_field%grid_loc

         if (io_field%missing_value /= undefined) then
            work_line = char_blank
            work_line( 1:13) = 'missing_value'
            work_line(14:14) = attrib_separator
            work_line(15:16) = 'r8'
            work_line(17:17) = attrib_separator
            write(work_line(18:),*) io_field%missing_value
            write(unit,'(a)') trim(work_line)
         endif

         if (any(io_field%valid_range(:) /= undefined)) then
            work_line = char_blank
            work_line( 1:11) = 'valid_range'
            work_line(12:12) = attrib_separator
            work_line(13:14) = 'r8'
            work_line(15:15) = attrib_separator
            write(work_line(16:),*) io_field%valid_range(:)
            write(unit,'(a)') trim(work_line)
         endif

         work_line = char_blank
         work_line(1:2) = 'id'
         work_line(3:3) = attrib_separator
         work_line(4:6) = 'int'
         work_line(7:7) = attrib_separator
         write(work_line(8:),*) io_field%id
         write(unit,'(a)') trim(work_line)

         !**** number of dimensions and dimension info

         work_line = char_blank
         work_line( 1:11) = 'nfield_dims'
         work_line(12:12) = attrib_separator
         work_line(13:15) = 'int'
         work_line(16:16) = attrib_separator
         write(work_line(17:),*) io_field%nfield_dims
         write(unit,'(a)') trim(work_line)

         !**** not sure what to do about dimension info

         ! write field_dim stuff

         !*** write additional attributes to header file

         if (associated(io_field%add_attrib_cval)) then
            do n=1,size(io_field%add_attrib_cval)

               work_line = char_blank
               cindx1 = len_trim(io_field%add_attrib_cname(n))
               cindx2 = cindx1 + 7
               work_line(1:cindx1) = trim(io_field%add_attrib_cname(n))
               work_line(cindx1+1:cindx1+1) = attrib_separator 
               work_line(cindx1+2:cindx2-2) = 'char'
               work_line(cindx2-1:cindx2-1) = attrib_separator 
               cindx1 = cindx2 + len_trim(io_field%add_attrib_cval(n))
               work_line(cindx2:cindx2) = trim(io_field%add_attrib_cval(n))
               write(unit,'(a)') trim(work_line)
            end do
         endif

         if (associated(io_field%add_attrib_lval)) then
            do n=1,size(io_field%add_attrib_lval)
   
               work_line = char_blank
               cindx1 = len_trim(io_field%add_attrib_lname(n))
               cindx2 = cindx1 + 6
               work_line(1:cindx1) = trim(io_field%add_attrib_lname(n))
               work_line(cindx1+1:cindx1+1) = attrib_separator 
               work_line(cindx1+2:cindx2-2) = 'log'
               work_line(cindx2-1:cindx2-1) = attrib_separator 
               write(work_line(cindx2:),*) io_field%add_attrib_lval(n)
               write(unit,'(a)') trim(work_line)
            end do
         endif

         if (associated(io_field%add_attrib_ival)) then
            do n=1,size(io_field%add_attrib_ival)
   
               work_line = char_blank
               cindx1 = len_trim(io_field%add_attrib_iname(n))
               cindx2 = cindx1 + 6
               work_line(1:cindx1) = trim(io_field%add_attrib_iname(n))
               work_line(cindx1+1:cindx1+1) = attrib_separator 
               work_line(cindx1+2:cindx2-2) = 'int'
               work_line(cindx2-1:cindx2-1) = attrib_separator 
               write(work_line(cindx2:),*) io_field%add_attrib_ival(n)
               write(unit,'(a)') trim(work_line)
            end do
         endif

         if (associated(io_field%add_attrib_rval)) then
            do n=1,size(io_field%add_attrib_rval)
   
               work_line = char_blank
               cindx1 = len_trim(io_field%add_attrib_rname(n))
               cindx2 = cindx1 + 5
               work_line(1:cindx1) = trim(io_field%add_attrib_rname(n))
               work_line(cindx1+1:cindx1+1) = attrib_separator 
               work_line(cindx1+2:cindx2-2) = 'r4'
               work_line(cindx2-1:cindx2-1) = attrib_separator 
               write(work_line(cindx2:),*) io_field%add_attrib_rval(n)
               write(unit,'(a)') trim(work_line)
            end do
         endif

         if (associated(io_field%add_attrib_dval)) then
            do n=1,size(io_field%add_attrib_dval)
   
               work_line = char_blank
               cindx1 = len_trim(io_field%add_attrib_dname(n))
               cindx2 = cindx1 + 5
               work_line(1:cindx1) = trim(io_field%add_attrib_dname(n))
               work_line(cindx1+1:cindx1+1) = attrib_separator 
               work_line(cindx1+2:cindx2-2) = 'r8'
               work_line(cindx2-1:cindx2-1) = attrib_separator 
               write(work_line(cindx2:),*) io_field%add_attrib_dval(n)
               write(unit,'(a)') work_line
               write(unit,'(a)') trim(work_line)
            end do
         endif

         !*** end variable attribute section with a delimiter

         write(unit,'(a1)') '/'

      endif ! master task

!-----------------------------------------------------------------------
!
!  attempt to define an io field from an input header file
!
!-----------------------------------------------------------------------

   else ! this is an input file

      unit = data_file%id(2)
      if (unit <= 0) then  ! no header file, assume fields are defined 
                           ! in the order they exist in input file
                           ! set id as current record and increment

         io_field%id = data_file%current_record

         if (associated(io_field%field_i_2d)) then
            data_file%current_record = data_file%current_record + 1
         else if (associated(io_field%field_i_3d)) then
            data_file%current_record = data_file%current_record + &
                                       size(io_field%field_i_3d, dim=3)
         else if (associated(io_field%field_r_2d)) then
            data_file%current_record = data_file%current_record + 1
         else if (associated(io_field%field_r_3d)) then
            data_file%current_record = data_file%current_record + &
                                       size(io_field%field_r_3d, dim=3)
         else if (associated(io_field%field_d_2d)) then
            data_file%current_record = data_file%current_record + 1
         else if (associated(io_field%field_d_3d)) then
            data_file%current_record = data_file%current_record + &
                                       size(io_field%field_d_3d, dim=3)
         else
            call exit_POP(sigAbort, &
             'define: No known binary field descriptor associated')
         end if

         
      else  ! header exists: read all attributes from header file
         
         if (my_task == master_task) then
            hdr_error = 0
            rewind (unit)
            cindx1 = len_trim(io_field%short_name) + 1
            comp_line(1:1) = '&'
            comp_line(2:cindx1) = trim(io_field%short_name)

            srch_loop: do while (hdr_error == 0) ! look for field in file

               read(unit,'(a80)',iostat=hdr_error) work_line
               work_line = adjustl(work_line)
               if (work_line(1:cindx1) == comp_line(1:cindx1)) &
                  exit srch_loop
            end do srch_loop
         endif ! master_task

         call broadcast_scalar(hdr_error, master_task)
         if (hdr_error /= 0) call exit_POP(sigAbort,&
            'could not find field in binary header file')

         !*** found io_field definition - now extract field attributes

         att_loop: do

            if (my_task == master_task) then
               read(unit,'(a80)',iostat=hdr_error) work_line
               work_line = adjustl(work_line)
            endif

            call broadcast_scalar(hdr_error, master_task)
            call broadcast_scalar(work_line, master_task)

            !*** exit if error reading or reached end of definition

            if (hdr_error /= 0 .or. work_line(1:1) == '/') &
               exit att_loop

            !*** extract attribute name

            cindx1 = index(work_line,attrib_separator)
            att_name = char_blank
            att_name(1:cindx1-1) = work_line(1:cindx1-1)

            !*** strip attribute name from input line

            comp_line = work_line
            do n=1,cindx1
               comp_line(n:n) = ' '
            end do
            work_line = adjustl(comp_line)

            !*** extract attribute data type

            cindx1 = index(work_line,attrib_separator)
            ctype = char_blank
            ctype(1:cindx1-1) = work_line(1:cindx1-1)

            !*** strip data type from input line

            comp_line = work_line
            do n=1,cindx1
               comp_line(n:n) = ' '
            end do
            work_line = adjustl(comp_line)

            !*** if this attribute matches a standard attribute
            !*** set attribute value.  otherwise, define an
            !*** additional attribute

            select case (trim(att_name))

            case('long_name','LONG_NAME')
               io_field%long_name = char_blank
               io_field%long_name = trim(work_line)

            case('units','UNITS')
               io_field%units = char_blank
               io_field%units = trim(work_line)

            case('coordinates','COORDINATES')
               io_field%coordinates = char_blank
               io_field%coordinates = trim(work_line)

            case('grid_loc','GRID_LOC')
               io_field%grid_loc = char_blank
               io_field%grid_loc = trim(work_line)

            case('missing_value','MISSING_VALUE')
               read(work_line,*) io_field%missing_value

            case('valid_range','VALID_RANGE')
               read(work_line,*) io_field%valid_range(:)

            case('id','ID')
               read(work_line,*) io_field%id

            case('field_dim','FIELD_DIM')
               !**** do not know what to do about dimensions yet

            case('nfield_dims','NFIELD_DIMS')
               read(work_line,*) io_field%nfield_dims

            case default

               !*** add attribute to io field

               select case (ctype)
               case ('char','CHAR','character','CHARACTER')
                  call add_attrib_io_field(io_field, trim(att_name), &
                                                     trim(work_line))

               case ('log','LOG','logical','LOGICAL')
                  read(work_line,*) att_lval
                  call add_attrib_io_field(io_field, trim(att_name), &
                                                     att_lval)

               case ('int','INT','i4','I4','integer','INTEGER')
                  read(work_line,*) att_ival
                  call add_attrib_io_field(io_field, trim(att_name), &
                                                     att_ival)

               case ('r4','R4','real','REAL','float','FLOAT')
                  read(work_line,*) att_rval
                  call add_attrib_io_field(io_field, trim(att_name), &
                                                     att_rval)

               case ('r8','R8','dbl','DBL','double','DOUBLE')
                  read(work_line,*) att_dval
                  call add_attrib_io_field(io_field, trim(att_name), &
                                                     att_dval)

               case default
                  call exit_POP(sigAbort, &
                                'define_io_field: unknown data type')
               end select

            end select ! att_name

         end do att_loop

         if (hdr_error /= 0) call exit_POP(sigAbort, &
            'define_io_field: error reading attribute from header')

      endif ! unit > 0 (header exists)

   endif !readonly

!-----------------------------------------------------------------------
!EOC

 end subroutine define_field_binary

!***********************************************************************
!BOP
! !IROUTINE: write_field_binary
! !INTERFACE:


 subroutine write_field_binary(data_file, io_field) 1,14

! !INPUT PARAMETERS:

   type (datafile), intent (inout)  :: &
      data_file             ! file to which data will be written

   type (io_field_desc), intent (inout) :: &
      io_field              ! field to be written

! !DESCRIPTION:
!  This routine writes a binary field to the data file.
!
! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (i4) :: &
      record        ! starting record for writing array

!-----------------------------------------------------------------------
!
!  extract the starting record from the io field
!
!-----------------------------------------------------------------------

   record = io_field%id

!-----------------------------------------------------------------------
!
!  write the io field
!
!-----------------------------------------------------------------------

#ifdef USEPIO
   if (associated(io_field%field_i_2d)) then
      call write_array(data_file,io_field%varDesc, io_field%ioDesc, io_field%field_i_2d,record)
   else if (associated(io_field%field_i_3d)) then
      call write_array(data_file,io_field%varDesc, io_field%ioDesc, io_field%field_i_3d,record)
   else if (associated(io_field%field_r_2d)) then
      call write_array(data_file,io_field%varDesc, io_field%ioDesc, io_field%field_r_2d,record)
   else if (associated(io_field%field_r_3d)) then
      call write_array(data_file,io_field%varDesc, io_field%ioDesc, io_field%field_r_3d,record)
   else if (associated(io_field%field_d_2d)) then
      call write_array(data_file,io_field%varDesc, io_field%ioDesc, io_field%field_d_2d,record)
   else if (associated(io_field%field_d_3d)) then
      call write_array(data_file,io_field%varDesc, io_field%ioDesc, io_field%field_d_3d,record)
   else
      call exit_POP(sigAbort, &
                   'write: No known binary field descriptor associated')
   end if
#else
   if (associated(io_field%field_i_2d)) then
      call write_array(data_file,io_field%field_i_2d,record)
   else if (associated(io_field%field_i_3d)) then
      call write_array(data_file,io_field%field_i_3d,record)
   else if (associated(io_field%field_r_2d)) then
      call write_array(data_file,io_field%field_r_2d,record)
   else if (associated(io_field%field_r_3d)) then
      call write_array(data_file,io_field%field_r_3d,record)
   else if (associated(io_field%field_d_2d)) then
      call write_array(data_file,io_field%field_d_2d,record)
   else if (associated(io_field%field_d_3d)) then
      call write_array(data_file,io_field%field_d_3d,record)
   else
      call exit_POP(sigAbort, &
                   'write: No known binary field descriptor associated')
   end if
#endif

!-----------------------------------------------------------------------
!EOC

 end subroutine write_field_binary

!***********************************************************************
!BOP
! !IROUTINE: read_field_binary
! !INTERFACE:


 subroutine read_field_binary(data_file, io_field) 1,14

! !DESCRIPTION:
!  This routine reads a field from a binary input file.
!
! !REVISION HISTORY:
!  same as module

! !INPUT/OUTPUT PARAMETERS:

   type (datafile), intent (inout)  :: &
      data_file             ! file from which to read data

   type (io_field_desc), intent (inout) :: &
      io_field              ! data to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  if the field location parameters undefined, assume center scalar
!
!-----------------------------------------------------------------------

   if (io_field%field_loc == field_loc_unknown) then
      io_field%field_loc  = field_loc_center
      io_field%field_type = field_type_scalar
   endif

!-----------------------------------------------------------------------
!
!  read an io field
!
!-----------------------------------------------------------------------

#ifdef USEPIOREAD
   if (associated(io_field%field_i_2d)) then
      call read_array(data_file,io_field%varDesc, io_field%ioDesc, io_field%field_i_2d, &
                io_field%id, io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_i_3d)) then
      call read_array(data_file,io_field%varDesc, io_field%ioDesc, io_field%field_i_3d, &
                io_field%id, io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_r_2d)) then
      call read_array(data_file,io_field%varDesc, io_field%ioDesc, io_field%field_r_2d, &
                io_field%id, io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_r_3d)) then
      call read_array(data_file,io_field%varDesc, io_field%ioDesc, io_field%field_r_3d, &
                io_field%id, io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_d_2d)) then
      call read_array(data_file, io_field%varDesc, io_field%ioDesc, io_field%field_d_2d, &
                io_field%id, io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_d_3d)) then
      call read_array(data_file, io_field%varDesc, io_field%ioDesc, io_field%field_d_3d, &
                io_field%id, io_field%field_loc, io_field%field_type)
   else
      call exit_POP(sigAbort,'read_field: field not associated')
   end if
#else
   if (associated(io_field%field_i_2d)) then
      call read_array(data_file,io_field%field_i_2d,io_field%id, &
                                io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_i_3d)) then
      call read_array(data_file,io_field%field_i_3d,io_field%id, &
                                io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_r_2d)) then
      call read_array(data_file,io_field%field_r_2d,io_field%id, &
                                io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_r_3d)) then
      call read_array(data_file,io_field%field_r_3d,io_field%id, &
                                io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_d_2d)) then
      call read_array(data_file,io_field%field_d_2d,io_field%id, &
                                io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_d_3d)) then
      call read_array(data_file,io_field%field_d_3d,io_field%id, &
                                io_field%field_loc, io_field%field_type)
   else
      call exit_POP(sigAbort,'read_field: field not associated')
   end if
#endif

!-----------------------------------------------------------------------
!EOC

 end subroutine read_field_binary

!***********************************************************************
!BOP
! !IROUTINE: read_int_2d
! !INTERFACE:

#ifdef USEPIOREAD

 subroutine read_int_2d(data_file, varDesc, ioDesc,  INT2D, start_record, & 1,1
                        field_loc, field_type)
#else

 subroutine read_int_2d(data_file, INT2D, start_record, & 1,1
                        field_loc, field_type)
#endif

! !DESCRIPTION:
!  Reads a 2-d horizontal slice of integers from a binary file
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   type (datafile), intent(inout) :: &
      data_file                   ! info about data file

#ifdef USEPIOREAD
   type (VAR_desc_t), intent(inout) :: varDesc
   type (IO_desc_t),  intent(inout) :: ioDesc
#endif

   integer (i4), intent(in) :: &
      start_record,            &! starting record of field in file
      field_loc,               &! loc of field on horiz grid
      field_type                ! type of field (scalar, vector, angle)

! !INPUT/OUTPUT PARAMETERS:

   integer (i4), dimension(:,:,:), intent(inout) :: &
      INT2D                        ! array to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (i4), dimension(:,:), allocatable ::  &
      IOBUFI                ! local global-sized buffer

#ifdef USEPIOREAD
   integer(i4) :: llen
   integer(i4) :: iostat
   integer(i4), dimension(:),allocatable :: lbuf_i4

!DBG      write(*,*) 'read_int_2d'
      llen = PIO_numToRead(ioDesc)
      allocate(lbuf_i4(llen))
      varDesc%rec = start_record
      call PIO_read_darray(data_file%file,varDesc,ioDesc,lbuf_i4,iostat)
      call Copyto2D(lbuf_i4,INT2d)
!NEWPIO      call update_ghost_cells(INT2d,bndy_clinic,field_loc,field_type)
      deallocate(lbuf_i4)
#else

!-----------------------------------------------------------------------
!
!  read in global 2-d slice from one processor
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then
      allocate(IOBUFI(nx_global,ny_global))
      read(data_file%id(1),rec=start_record) IOBUFI
   endif

!-----------------------------------------------------------------------
!
!  send chunks to processors who own them
!
!-----------------------------------------------------------------------

   call scatter_global(INT2D, IOBUFI, master_task, distrb_clinic, &
                       field_loc, field_type)

   if (my_task == master_task) deallocate(IOBUFI)
#endif

!-----------------------------------------------------------------------
!EOC

 end subroutine read_int_2d

!***********************************************************************
!BOP
! !IROUTINE: read_real4_2d
! !INTERFACE:

#ifdef USEPIOREAD

 subroutine read_real4_2d(data_file, varDesc, ioDesc, REAL2D, start_record, &,1
                          field_loc, field_type)
#else

 subroutine read_real4_2d(data_file, REAL2D, start_record, &,1
                          field_loc, field_type)
#endif

! !DESCRIPTION:
!  Reads a 2-d horizontal slice of reals from a binary file
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   type (datafile), intent(inout) :: &
      data_file                   ! info about data file

#ifdef USEPIOREAD
   type (Var_desc_t), intent(inout) :: varDesc
   type (IO_desc_t), intent(inout)  :: ioDesc
#endif

   integer (i4), intent(in) :: &
      start_record,            &! starting record of field in file
      field_loc,               &! loc of field on horiz grid
      field_type                ! type of field (scalar, vector, angle)

! !INPUT/OUTPUT PARAMETERS:

   real (r4), dimension(:,:,:), intent(inout) :: &
      REAL2D                    ! array to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   real (r4), dimension(:,:), allocatable ::  &
      IOBUFR                ! local global-sized buffer

#ifdef USEPIOREAD
   integer(i4) :: llen
   integer(i4) :: iostat
   real(r4), dimension(:),allocatable :: lbuf_r4

!      write(*,*) 'read_real4_2d'
      llen = PIO_numToRead(ioDesc)
      allocate(lbuf_r4(llen))
      varDesc%rec = start_record
      call PIO_read_darray(data_file%file,varDesc,ioDesc, lbuf_r4,iostat)
      call Copyto2D(lbuf_r4,REAL2d)
!NEWPIO      call update_ghost_cells(REAL2d,bndy_clinic,field_loc,field_type)
      deallocate(lbuf_r4)

#else

!-----------------------------------------------------------------------
!
!  read in global 2-d slice from one processor
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then
      allocate(IOBUFR(nx_global,ny_global))
      read(data_file%id(1),rec=start_record) IOBUFR
   endif

!-----------------------------------------------------------------------
!
!  send chunks to processors who own them
!
!-----------------------------------------------------------------------

   call scatter_global(REAL2D, IOBUFR, master_task, distrb_clinic, &
                       field_loc, field_type)

   if (my_task == master_task) deallocate(IOBUFR)
#endif

!-----------------------------------------------------------------------
!EOC

 end subroutine read_real4_2d

!***********************************************************************
!BOP
! !IROUTINE: read_real8_2d
! !INTERFACE:

#ifdef USEPIOREAD

 subroutine read_real8_2d(data_file, varDesc,ioDesc,  DBL2D, start_record, &,1
                          field_loc, field_type)
#else

 subroutine read_real8_2d(data_file, DBL2D, start_record, &,1
                          field_loc, field_type)
#endif

! !DESCRIPTION:
!  Reads a 2-d horizontal slice of 64-bit reals from a binary file
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   type (datafile), intent(inout) :: &
      data_file                   ! info about data file

#ifdef USEPIOREAD
   type (VAR_desc_t), intent(inout) :: varDesc
   type (IO_desc_t), intent(inout) :: ioDesc
#endif

   integer (i4), intent(in) :: &
      start_record,            &! starting record of field in file
      field_loc,               &! loc of field on horiz grid
      field_type                ! type of field (scalar, vector, angle)

! !INPUT/OUTPUT PARAMETERS:

   real (r8), dimension(:,:,:), intent(inout) :: &
      DBL2D                        ! array to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   real (r8), dimension(:,:), allocatable ::  &
      IOBUFD                ! local global-sized buffer

#ifdef USEPIOREAD
   integer(i4) :: llen
   integer(i4) :: iostat
   real(r8), dimension(:),allocatable :: lbuf_r8

!      write(*,*) 'read_real8_2d'
      llen = PIO_numToRead(ioDesc)
      allocate(lbuf_r8(llen))
      varDesc%rec = start_record
      call PIO_read_darray(data_file%file,varDesc,ioDesc, lbuf_r8,iostat)
      call Copyto2D(lbuf_r8,DBL2d)
      deallocate(lbuf_r8)
!NEWPIO      call update_ghost_cells(DBL2d,bndy_clinic,field_loc,field_type)
#else

!-----------------------------------------------------------------------
!
!  read in global 2-d slice from one processor
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then
      allocate(IOBUFD(nx_global,ny_global))
      read(data_file%id(1),rec=start_record) IOBUFD
   endif

!-----------------------------------------------------------------------
!
!  send chunks to processors who own them
!
!-----------------------------------------------------------------------

   call scatter_global(DBL2D, IOBUFD, master_task, distrb_clinic, &
                       field_loc, field_type)

   if (my_task == master_task) deallocate(IOBUFD)
#endif

!-----------------------------------------------------------------------
!EOC

 end subroutine read_real8_2d

!***********************************************************************
!BOP
! !IROUTINE: read_int_3d
! !INTERFACE:

#ifdef USEPIOREAD

 subroutine read_int_3d(data_file, varDesc, ioDesc, INT3D, start_record, &,1
                        field_loc, field_type)
#else

 subroutine read_int_3d(data_file, INT3D, start_record, &,1
                        field_loc, field_type)
#endif

! !DESCRIPTION:
!  Reads a 3-d integer array from a binary file
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   type (datafile), intent(inout) :: &
      data_file                   ! info on input data file

#ifdef USEPIOREAD
   type (VAR_desc_t), intent(inout) :: varDesc
   type (IO_desc_t), intent(inout) :: ioDesc
#endif

   integer (i4), intent(in) :: &
      start_record,            &! starting record of field in file
      field_loc,               &! loc of field on horiz grid
      field_type                ! type of field (scalar, vector, angle)

! !INPUT/OUTPUT PARAMETERS:

   integer (i4), dimension(:,:,:,:), intent(inout) ::  &
      INT3D  ! array to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (i4) k,n,krecs,klvl,nz

#ifdef USEPIOREAD
   integer(i4) :: llen
   integer(i4), allocatable :: lbuf_i4(:)
   integer(i4) :: iostat
#endif

   integer (i4), dimension(:,:), allocatable ::  &
      IOBUFI                ! global-sized array buffer

!-----------------------------------------------------------------------
!
!  determine the number of records each i/o process must read to
!  get all the records
!
!-----------------------------------------------------------------------

   nz = size(INT3D, DIM=3)

#ifdef USEPIOREAD

!      write(*,*) 'read_int_3d'
   llen = PIO_numToRead(ioDesc)
   allocate(lbuf_i4(llen))
   do k=1,nz
      varDesc%rec = start_record + (k-1)
      call PIO_read_darray(data_file%file,varDesc,ioDesc,lbuf_i4,iostat)
      call Copyto2D(lbuf_i4,INT3d(:,:,k,:))
   enddo
   deallocate(lbuf_i4)
!NEWPIO   call update_ghost_cells(INT3d,bndy_clinic,field_loc,field_type)
#else

   if (mod(nz,data_file%num_iotasks) == 0) then
      krecs = nz/data_file%num_iotasks
   else
      krecs = nz/data_file%num_iotasks + 1
   endif

   if (my_task < data_file%num_iotasks) &
      allocate(IOBUFI(nx_global,ny_global))

!-----------------------------------------------------------------------
!
!  each i/o process reads a horizontal slab from a record and sends
!  the data to processors who own it. read and distribute num_iotasks
!  records at a time to keep the messages from getting mixed up.
!
!-----------------------------------------------------------------------

   do k=1,krecs

      if (my_task < data_file%num_iotasks) then
         klvl = (k-1)*data_file%num_iotasks + my_task + 1
         if (klvl <= nz) then
            read(data_file%id(1), rec=start_record+klvl-1) IOBUFI
         endif
      endif

      do n=1,data_file%num_iotasks
         klvl = (k-1)*data_file%num_iotasks + n
         if (klvl <= nz) then
            call scatter_global(INT3D(:,:,klvl,:), IOBUFI, n-1, &
                                distrb_clinic, field_loc, field_type)
         endif
      end do

   end do

   if (my_task < data_file%num_iotasks) deallocate(IOBUFI)
#endif

!-----------------------------------------------------------------------
!EOC

 end subroutine read_int_3d

!***********************************************************************
!BOP
! !IROUTINE: read_real4_3d
! !INTERFACE:

#ifdef USEPIOREAD

 subroutine read_real4_3d(data_file, varDesc, ioDesc, REAL3D, start_record, &,1
                          field_loc, field_type)
#else

 subroutine read_real4_3d(data_file, REAL3D, start_record, &,1
                          field_loc, field_type)
#endif

! !DESCRIPTION:
!  Reads a 3-d real array from a binary file
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   type (datafile), intent(inout) :: &
      data_file                   ! info on input data file

#ifdef USEPIOREAD
   type (VAR_desc_t), intent(inout) :: varDesc
   type (IO_desc_t), intent(inout) :: ioDesc
#endif

   integer (i4), intent(in) :: &
      start_record,            &! starting record of field in file
      field_loc,               &! loc of field on horiz grid
      field_type                ! type of field (scalar, vector, angle)

! !INPUT/OUTPUT PARAMETERS:

   real (r4), dimension(:,:,:,:), intent(inout) ::  &
      REAL3D  ! array to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (i4) k,n,krecs,klvl,nz

#ifdef USEPIOREAD
   integer(i4) :: llen
   integer(i4) :: iostat
   real(r4), allocatable :: lbuf_r4(:)
#endif

   real (r4), dimension(:,:), allocatable ::  &
      IOBUFR                ! global-sized array buffer

!-----------------------------------------------------------------------
!
!  determine the number of records each i/o process must read to
!  get all the records
!
!-----------------------------------------------------------------------

   nz = size(REAL3D, DIM=3)

#ifdef USEPIOREAD
!      write(*,*) 'read_real4_3d'
   llen = PIO_numToRead(ioDesc)
   allocate(lbuf_r4(llen))
   do k=1,nz
      varDesc%rec = start_record + (k-1)
      call PIO_read_darray(data_file%file,varDesc,ioDesc, lbuf_r4,iostat)
      call Copyto2D(lbuf_r4,REAL3D(:,:,k,:))
   enddo
   deallocate(lbuf_r4)
!NEWPIO   call update_ghost_cells(REAL3D,bndy_clinic,field_loc,field_type)
#else

   if (mod(nz,data_file%num_iotasks) == 0) then
      krecs = nz/data_file%num_iotasks
   else
      krecs = nz/data_file%num_iotasks + 1
   endif

   if (my_task < data_file%num_iotasks) &
      allocate(IOBUFR(nx_global,ny_global))

!-----------------------------------------------------------------------
!
!  each i/o process reads a horizontal slab from a record and sends
!  the data to processors who own it. read and distribute num_iotasks
!  records at a time to keep the messages from getting mixed up.
!
!-----------------------------------------------------------------------

   do k=1,krecs

      if (my_task < data_file%num_iotasks) then
         klvl = (k-1)*data_file%num_iotasks + my_task + 1
         if (klvl <= nz) then
            read(data_file%id(1), rec=start_record+klvl-1) IOBUFR
         endif
      endif

      do n=1,data_file%num_iotasks
         klvl = (k-1)*data_file%num_iotasks + n
         if (klvl <= nz) then
            call scatter_global(REAL3D(:,:,klvl,:), IOBUFR, n-1, &
                                distrb_clinic, field_loc, field_type)
         endif
      end do

   end do

   if (my_task < data_file%num_iotasks) deallocate(IOBUFR)
#endif

!-----------------------------------------------------------------------
!EOC

 end subroutine read_real4_3d

!***********************************************************************
!BOP
! !IROUTINE: read_real8_3d
! !INTERFACE:

#ifdef USEPIOREAD

 subroutine read_real8_3d(data_file, varDesc, ioDesc, DBL3D, start_record, &,1
                          field_loc, field_type)
#else

 subroutine read_real8_3d(data_file, DBL3D, start_record, &,1
                          field_loc, field_type)
#endif

! !DESCRIPTION:
!  Reads a 3-d 64-bit real array from a binary file
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   type (datafile), intent(inout) :: &
      data_file                   ! info on input data file

#ifdef USEPIOREAD
   type (VAR_desc_t), intent(inout) :: varDesc
   type (IO_desc_t), intent(inout) :: ioDesc
#endif

   integer (i4), intent(in) :: &
      start_record,            &! starting record of field in file
      field_loc,               &! loc of field on horiz grid
      field_type                ! type of field (scalar, vector, angle)

! !INPUT/OUTPUT PARAMETERS:

   real (r8), dimension(:,:,:,:), intent(inout) ::  &
      DBL3D  ! array to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (i4) k,n,krecs,klvl,nz

#ifdef USEPIOREAD
   integer(i4)  :: llen,iostat
   real (r8), allocatable, dimension(:) :: lbuf_r8
   real (r8), dimension(:,:,:), allocatable :: DBL2D
#endif

   real (r8), dimension(:,:), allocatable ::  &
      IOBUFD                ! global-sized array buffer

!-----------------------------------------------------------------------
!
!  determine the number of records each i/o process must read to
!  get all the records
!
!-----------------------------------------------------------------------

   nz = size(DBL3D, DIM=3)

#ifdef USEPIOREAD
!      write(*,*) 'read_real8_3d'
   llen = PIO_numToRead(ioDesc)
   allocate(lbuf_r8(llen))
   do k=1,nz
      varDesc%rec = start_record + (k-1)
      call PIO_read_darray(data_file%file,varDesc,ioDesc,lbuf_r8,iostat)
      call Copyto2D(lbuf_r8,DBL3d(:,:,k,:))
   enddo
   deallocate(lbuf_r8)
!NEWPIO   call update_ghost_cells(DBL3d,bndy_clinic,field_loc,field_type)
#else

   if (mod(nz,data_file%num_iotasks) == 0) then
      krecs = nz/data_file%num_iotasks
   else
      krecs = nz/data_file%num_iotasks + 1
   endif

   if (my_task < data_file%num_iotasks) &
      allocate(IOBUFD(nx_global,ny_global))

!-----------------------------------------------------------------------
!
!  each i/o process reads a horizontal slab from a record and sends
!  the data to processors who own it. read and distribute num_iotasks
!  records at a time to keep the messages from getting mixed up.
!
!-----------------------------------------------------------------------

   do k=1,krecs

      if (my_task < data_file%num_iotasks) then
         klvl = (k-1)*data_file%num_iotasks + my_task + 1
         if (klvl <= nz) then
            read(data_file%id(1), rec=start_record+klvl-1) IOBUFD
         endif
      endif

      do n=1,data_file%num_iotasks
         klvl = (k-1)*data_file%num_iotasks + n
         if (klvl <= nz) then
            call scatter_global(DBL3D(:,:,klvl,:), IOBUFD, n-1, &
                                distrb_clinic, field_loc, field_type)
         endif
      end do

   end do

   if (my_task < data_file%num_iotasks) deallocate(IOBUFD)
#endif

!-----------------------------------------------------------------------
!EOC

 end subroutine read_real8_3d

!***********************************************************************
!BOP
! !IROUTINE: write_int_2d
! !INTERFACE:

#ifdef USEPIO

 subroutine write_int_2d(data_file,varDesc, ioDesc, INT2D,start_record) 1,1
#else

 subroutine write_int_2d(data_file,INT2D,start_record) 1,1
#endif

! !DESCRIPTION:
!  Writes a 2-d slab of integers to a binary file.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   type (datafile), intent(inout) :: &
      data_file                   ! file information

#ifdef USEPIO
   type (VAR_desc_t), intent(inout) :: varDesc
   type (IO_desc_t), intent(inout) :: ioDesc
#endif

   integer (i4), intent(in) ::  &
      start_record              ! starting record of array in file

   integer (i4), dimension(:,:,:), intent(in) ::  &
      INT2D                 ! array to be written

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (i4), dimension(:,:), allocatable ::  &
      IOBUFI                ! local global-sized buffer

#ifdef USEPIO
   integer(i4) :: llen
   integer(i4), dimension(:), allocatable :: lbuf_i4
   integer(i4) :: iostat

     llen = PIO_numToWrite(ioDesc)
     allocate(lbuf_i4(llen))
     call Copyto1D(INT2D,lbuf_i4)
     varDesc%rec = start_record
     call PIO_write_darray(data_file%file,varDesc,ioDesc,lbuf_i4,iostat)
     deallocate(lbuf_i4)
#else
!-----------------------------------------------------------------------
!
!  receive chunks from processors who own them
!
!-----------------------------------------------------------------------

   if (my_task == master_task) allocate(IOBUFI(nx_global,ny_global))

   call gather_global(IOBUFI, INT2D, master_task, distrb_clinic)

!-----------------------------------------------------------------------
!
!  write global 2-d slice from one processor
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then
      write(data_file%id(1),rec=start_record) IOBUFI
      deallocate(IOBUFI)
   endif
#endif

!-----------------------------------------------------------------------
!EOC
 end subroutine write_int_2d

!***********************************************************************
!BOP
! !IROUTINE: write_real4_2d
! !INTERFACE:
#ifdef USEPIO

 subroutine write_real4_2d(data_file,varDesc,ioDesc, REAL2D,start_record),1
#else

 subroutine write_real4_2d(data_file,REAL2D,start_record),1
#endif

! !DESCRIPTION:
!  Writes a 2-d slab of reals to a binary file.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   type (datafile), intent(inout) :: &
      data_file                   ! file information

#ifdef USEPIO
   type (VAR_desc_t), intent(inout) :: varDesc
   type (IO_desc_t), intent(inout) :: ioDesc
#endif

   integer (i4), intent(in) ::  &
      start_record              ! starting record of array in file

   real (r4), dimension(:,:,:), intent(in) ::  &
      REAL2D                 ! array to be written

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   real (r4), dimension(:,:), allocatable ::  &
      IOBUFR                ! local global-sized buffer

!-----------------------------------------------------------------------
!
!  receive chunks from processors who own them
!
!-----------------------------------------------------------------------
#ifdef USEPIO
   integer(i4) :: llen
   real(r4), dimension(:), allocatable :: lbuf_r4
   integer(i4) :: iostat

     llen = PIO_numToWrite(ioDesc)
     allocate(lbuf_r4(llen))
     lbuf_r4 = c0
!DBG     write(*,*) 'write_array(2D,R4): IAM: ',my_task,' rec: ',start_record
     call Copyto1D(REAL2D,lbuf_r4)
     varDesc%rec = start_record
     call PIO_write_darray(data_file%file,varDesc,ioDesc,lbuf_r4,iostat)
     deallocate(lbuf_r4)
#else

   if (my_task == master_task) allocate(IOBUFR(nx_global,ny_global))
   call gather_global(IOBUFR, REAL2D, master_task, distrb_clinic)

!-----------------------------------------------------------------------
!
!  write global 2-d slice from one processor
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then
      write(data_file%id(1),rec=start_record) IOBUFR
      deallocate(IOBUFR)
   endif
#endif

!-----------------------------------------------------------------------
!EOC
 end subroutine write_real4_2d

!***********************************************************************
!BOP
! !IROUTINE: write_real8_2d
! !INTERFACE:

#ifdef USEPIO

 subroutine write_real8_2d(data_file,varDesc, ioDesc, DBL2D,start_record),1
#else

 subroutine write_real8_2d(data_file,DBL2D,start_record),1
#endif

! !DESCRIPTION:
!  Writes a 2-d slab of doubles to a binary file.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   type (datafile), intent(inout) :: &
      data_file                   ! file information

#ifdef USEPIO
   type (VAR_desc_t), intent(inout) :: varDesc
   type (IO_desc_t), intent(inout) :: ioDesc
#endif

   integer (i4), intent(in) ::  &
      start_record              ! starting record of array in file

   real (r8), dimension(:,:,:), intent(in) ::  &
      DBL2D                 ! array to be written

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   real (r8), dimension(:,:), allocatable ::  &
      IOBUFD                ! local global-sized buffer

#ifdef USEPIO
   integer(i4) :: llen
   real(r8), dimension(:), allocatable :: lbuf_r8
   integer(i4) :: iostat
#endif
!-----------------------------------------------------------------------
!
!  receive chunks from processors who own them
!
!-----------------------------------------------------------------------
#ifdef USEPIO
     llen = PIO_numToWrite(ioDesc)
     allocate(lbuf_r8(llen))
     lbuf_r8 = c0
     call Copyto1D_real8(DBL2D,lbuf_r8)
     varDesc%rec = start_record
!     rsum = MAXVAL(lbuf_r8);
!     write(*,*) 'before PIO_write_darray IAM: ',my_task,' rec: ',start_record,rsum
     call PIO_write_darray(data_file%file,varDesc,ioDesc,lbuf_r8,iostat)
     deallocate(lbuf_r8)
#else

   if (my_task == master_task) allocate(IOBUFD(nx_global,ny_global))
   call gather_global(IOBUFD, DBL2D, master_task, distrb_clinic)

!-----------------------------------------------------------------------
!
!  write global 2-d slice from one processor
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then
      write(data_file%id(1),rec=start_record) IOBUFD
      deallocate(IOBUFD)
   endif
#endif

!-----------------------------------------------------------------------
!EOC
 end subroutine write_real8_2d

!***********************************************************************
!BOP
! !IROUTINE: write_int_3d
! !INTERFACE:

#ifdef USEPIO

 subroutine write_int_3d(data_file,varDesc,ioDesc,INT3D,start_record),1
#else

 subroutine write_int_3d(data_file,INT3D,start_record),1
#endif

! !DESCRIPTION:
!  Writes a 3-d integer array as a series of 2-d slabs to a binary
!  output file.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   type (datafile), intent(inout) :: &
      data_file                   ! info on output data file

#ifdef USEPIO
   type (VAR_desc_t), intent(inout) :: varDesc
   type (IO_desc_t), intent(inout) :: ioDesc
#endif

   integer (i4), intent(inout) :: &
      start_record             ! starting position of array in file

   integer (i4), dimension(:,:,:,:), intent(in) :: &
      INT3D                      ! array to be written

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (i4) ::  &
      k, n,   &  ! dummy counters
      krecs,  &  ! number of records each iotask must write
      klvl,   &  ! k index corresponding to current record
      nz         ! size of 3rd dimension of 3-d array to be written

   integer (i4), dimension(:,:), allocatable ::  &
      IOBUFI                ! global-sized buffer array

#ifdef USEPIO
   integer(i4) :: llen
   integer(i4), allocatable :: lbuf_i4(:)
   integer(i4) :: iostat
#endif
!-----------------------------------------------------------------------
!
!  determine the number of records each i/o process must write to
!  get all the records
!
!-----------------------------------------------------------------------

   nz = size(INT3D, DIM=3)

#ifdef USEPIO
     llen = PIO_numToWrite(ioDesc)
     allocate(lbuf_i4(llen))
     do k=1,nz
         varDesc%rec = start_record + (k-1)
         call Copyto1D(INT3D(:,:,k,:),lbuf_i4)
         call PIO_write_darray(data_file%file,varDesc,ioDesc,lbuf_i4,iostat)
     enddo
     deallocate(lbuf_i4)
#else
   if (mod(nz,data_file%num_iotasks) == 0) then
      krecs = nz/data_file%num_iotasks
   else
      krecs = nz/data_file%num_iotasks + 1
   endif

   if (my_task < data_file%num_iotasks) &
      allocate(IOBUFI(nx_global,ny_global))

!-----------------------------------------------------------------------
!
!  gather and write num_iotasks records at a time to keep the
!  messages from getting mixed up.
!
!-----------------------------------------------------------------------

   do k=1,krecs

!-----------------------------------------------------------------------
!
!     gather a global slice on each iotask
!
!-----------------------------------------------------------------------

      do n=1,data_file%num_iotasks
         klvl = (k-1)*data_file%num_iotasks + n
         if (klvl <= nz) then
            call gather_global(IOBUFI, INT3D(:,:,klvl,:), &
                               n-1, distrb_clinic)
         endif
      end do

!-----------------------------------------------------------------------
!
!     iotasks wait for the gather to be complete and write global
!     slabs to the file
!
!-----------------------------------------------------------------------

      if (my_task < data_file%num_iotasks) then
         klvl = (k-1)*data_file%num_iotasks + my_task + 1
         if (klvl <= nz) then
            write(data_file%id(1), rec=start_record+klvl-1) IOBUFI
         endif
      endif

   end do

!-----------------------------------------------------------------------

   if (my_task < data_file%num_iotasks) deallocate(IOBUFI)
#endif

!-----------------------------------------------------------------------
!EOC

 end subroutine write_int_3d

!***********************************************************************
!BOP
! !IROUTINE: write_real4_3d
! !INTERFACE:

#ifdef USEPIO

 subroutine write_real4_3d(data_file,varDesc,ioDesc, REAL3D,start_record),2
#else

 subroutine write_real4_3d(data_file,REAL3D,start_record),2
#endif

! !DESCRIPTION:
!  Writes a 3-d real array as a series of 2-d slabs to a binary
!  output file.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   type (datafile), intent(inout) :: &
      data_file                   ! info on output data file

#ifdef USEPIO
   type (VAR_desc_t), intent(inout) :: varDesc
   type (IO_desc_t), intent(inout) :: ioDesc
#endif

   integer (i4), intent(in) :: &
      start_record             ! starting position of array in file

   real (r4), dimension(:,:,:,:), intent(in) :: &
      REAL3D                      ! array to be written

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (i4) ::  &
      k, n,   &  ! dummy counters
      krecs,  &  ! number of records each iotask must write
      klvl,   &  ! k index corresponding to current record
      nz         ! size of 3rd dimension of 3-d array to be written

#ifdef USEPIO
    integer(i4) :: llen,ierr
    real(r4), allocatable :: lbuf_r4(:)
    integer(i4) :: iostat
#endif

   real (r4), dimension(:,:), allocatable ::  &
      IOBUFR                ! global-sized buffer array

!-----------------------------------------------------------------------
!
!  determine the number of records each i/o process must write to
!  get all the records
!
!-----------------------------------------------------------------------

   nz = size(REAL3D, DIM=3)

#ifdef USEPIO
     llen = PIO_numToWrite(ioDesc)
     allocate(lbuf_r4(llen))
     varDesc%rec = start_record
     do k=1,nz
         lbuf_r4 = c0
         varDesc%rec = start_record + (k-1)
!DBG         write(*,*) 'write_array(3D,R4): IAM: ',my_task,' rec: ',varDesc%rec
         call Copyto1D(REAL3D(:,:,k,:),lbuf_r4)
         call PIO_write_darray(data_file%file,varDesc,ioDesc,lbuf_r4,iostat)
         call POP_barrier
     enddo
     deallocate(lbuf_r4)
#else
   if (mod(nz,data_file%num_iotasks) == 0) then
      krecs = nz/data_file%num_iotasks
   else
      krecs = nz/data_file%num_iotasks + 1
   endif

   if (my_task < data_file%num_iotasks) &
      allocate(IOBUFR(nx_global,ny_global))

!-----------------------------------------------------------------------
!
!  gather and write num_iotasks records at a time to keep the
!  messages from getting mixed up.
!
!-----------------------------------------------------------------------

   do k=1,krecs

!-----------------------------------------------------------------------
!
!     gather a global slice on each iotask
!
!-----------------------------------------------------------------------

      do n=1,data_file%num_iotasks
         klvl = (k-1)*data_file%num_iotasks + n
         if (klvl <= nz) then
            call gather_global(IOBUFR, REAL3D(:,:,klvl,:), &
                               n-1, distrb_clinic)
         endif
      end do

!-----------------------------------------------------------------------
!
!     iotasks wait for the gather to be complete and write global
!     slabs to the file
!
!-----------------------------------------------------------------------

      if (my_task < data_file%num_iotasks) then
         klvl = (k-1)*data_file%num_iotasks + my_task + 1
         if (klvl <= nz) then
            write(data_file%id(1), rec=start_record+klvl-1) IOBUFR
         endif
      endif

   end do

!-----------------------------------------------------------------------

   if (my_task < data_file%num_iotasks) deallocate(IOBUFR)
#endif

!-----------------------------------------------------------------------
!EOC

 end subroutine write_real4_3d

!***********************************************************************
!BOP
! !IROUTINE: write_real8_3d
! !INTERFACE:

#ifdef USEPIO

 subroutine write_real8_3d(data_file,varDesc, ioDesc, DBL3D,start_record),1
#else

 subroutine write_real8_3d(data_file,DBL3D,start_record),1
#endif

! !DESCRIPTION:
!  Writes a 3-d 64-bit real array as a series of 2-d slabs to a binary
!  output file.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   type (datafile), intent(inout) :: &
      data_file                   ! info on output data file

#ifdef USEPIO
   type (VAR_desc_t), intent(inout) :: varDesc
   type (IO_desc_t), intent(inout) :: ioDesc
#endif

   integer (i4), intent(in) :: &
      start_record             ! starting position of array in file

   real (r8), dimension(:,:,:,:), intent(in) :: &
      DBL3D                      ! array to be written

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (i4) ::  &
      k, n,   &  ! dummy counters
      krecs,  &  ! number of records each iotask must write
      klvl,   &  ! k index corresponding to current record
      nz         ! size of 3rd dimension of 3-d array to be written

#ifdef USEPIO
   integer(i4) :: llen
   real(r8), allocatable :: lbuf_r8(:)
   integer(i4) :: iostat
#endif

   real (r8), dimension(:,:), allocatable ::  &
      IOBUFD                ! global-sized buffer array

!-----------------------------------------------------------------------
!
!  determine the number of records each i/o process must write to
!  get all the records
!
!-----------------------------------------------------------------------

   nz = size(DBL3D, DIM=3)

#ifdef USEPIO
!DBG     write(*,*) 'write_array(3D,R8): IAM: rec: ',data_file%file%comp_rank,start_record
     llen = PIO_numToWrite(ioDesc)
     allocate(lbuf_r8(llen))
     varDesc%rec = start_record
     do k=1,nz
         lbuf_r8 = c0
!DBG         if(my_task == 0) then
!DBG             write(*,*) 'IAM: ',my_task,' write_array(3D,R8) [slice]: ',k
!DBG         endif
         varDesc%rec = start_record + (k-1)
!JMD         call Copyto1D(DBL3D(:,:,k,:),lbuf_r8)
         call Copyto1D(DBL3D,k,lbuf_r8)
         call PIO_write_darray(data_file%file,varDesc,ioDesc,lbuf_r8,iostat)
     enddo
     deallocate(lbuf_r8)
#else

   if (mod(nz,data_file%num_iotasks) == 0) then
      krecs = nz/data_file%num_iotasks
   else
      krecs = nz/data_file%num_iotasks + 1
   endif

   if (my_task < data_file%num_iotasks) &
      allocate(IOBUFD(nx_global,ny_global))

!-----------------------------------------------------------------------
!
!  gather and write num_iotasks records at a time to keep the
!  messages from getting mixed up.
!
!-----------------------------------------------------------------------

   do k=1,krecs

!-----------------------------------------------------------------------
!
!     gather a global slice on each iotask
!
!-----------------------------------------------------------------------

      do n=1,data_file%num_iotasks
         klvl = (k-1)*data_file%num_iotasks + n
         if (klvl <= nz) then
            call gather_global(IOBUFD, DBL3D(:,:,klvl,:), &
                               n-1, distrb_clinic)
         endif
      end do

!-----------------------------------------------------------------------
!
!     iotasks wait for the gather to be complete and write global
!     slabs to the file
!
!-----------------------------------------------------------------------

      if (my_task < data_file%num_iotasks) then
         klvl = (k-1)*data_file%num_iotasks + my_task + 1
         if (klvl <= nz) then
            write(data_file%id(1), rec=start_record+klvl-1) IOBUFD
         endif
      endif

   end do

!-----------------------------------------------------------------------

   if (my_task < data_file%num_iotasks) deallocate(IOBUFD)
#endif

!-----------------------------------------------------------------------
!EOC

 end subroutine write_real8_3d

!***********************************************************************

 end module io_binary

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||