    module grid_io
!
!  SYNOPSIS
!    use grid_io
!
!  DESCRIPTION
!   GRID_IO contains routines for reading and writing the physical model grid
!   and Gyorgyi's observation formats.
!
!  ROUTINES IMPLEMENTED
!   subroutine read_grid_header(bfile,expect,south,north,psouth,pnorth,nsol,
!      nassim,seqno)
!   subroutine read_model_grid(grid,south,north,nsol)
!   subroutine read_nature_grid(header,truthfn,seq)
!   subroutine write_grid_slice(gfile,header,south,north,grid,lower,upper,nsol,
!      psouth,pnorth)
!   subroutine output_analysis(analysis,lower,upper,nsol,afile,header)
!   subroutine read_obs_slice(ofile,expect,nobs,obsarray)
!   subroutine write_obs_slice(ofile,header,south,north,obs,nobs)
!
!  AUTHOR
!   Eric Kostelich, kostelich@asu.edu, 480-965-5006
!
!  MODIFICATION HISTORY
!  08/08/05 - added DUMP_LOC debugging parameter
!  08/06/05 - added GRIDKIND parameter to adjust internal grid precision and 
!             IOKIND parameter to adjust i/o precision
!  08/01/05 - added rescaling of grid in OUTPUT_ANALYSIS; added 
!             single-precision i/o buffers; removed GRIDKIND parameter.
!  07/31/05 - removed i/o of observations
!  01/12/05 - raw_to_model_obs removed.  This version reads and writes
!             observations in the raw_obs format only.
!  12/14/04 - altered input/output grid indices for more efficient processing
!             by pre- and postanalysis.
!  11/10/04 - rewrite to correspond to new data assimilation array indexing.
!
!
!  REQUIRED PACKAGES
!
    use grid_indexing
    use kdtree
    use ferror
    use io
    use resolution
    implicit none
!
!  DUMP_LOC is the patch location (if any) where verbose output to DUMPUNIT
!  should occur.
!
    type(patch_index),save::dump_loc=patch_index(0,0,0)  ! default none
    contains
! ---------------------------------------------------------------------
    subroutine read_grid_header(bfile,expect,south,north,psouth,pnorth,nsol, &
       nassim,seqno)
!  READ_GRID_HEADER - read the header record from a model grid file created
!  by preanalysis.  This routine opens BFILE and KEEPS it open.
!  Side effects:
!   - opens INNOSWAP and keeps it open.
!   - positions BFILE to just before the record containing the model grid.
!    READ_GRID_HEADER may be followed immediately by a call to READ_MODEL_GRID.
!
!  Arguments:
!  BFILE   : name of input file
!  EXPECT  : 4-character header string to look for.  If set to 'none',
!     then we assume that the input file contains no header data at all
!     and holds the entire global grid.
!  SOUTH  := Southernmost latitude contained in the grid file.
!  NORTH  := Northernmost latitude index contained in the grid file.
!  PSOUTH, PNORTH := Latitudes between which data assimilation patches can
!     be centered.
!  NSOL  := number of ensemble solutions in the grid file
!  NASSIM := number of model variables available for assimilation in the file.
!  SEQNO := sequence number of the grid file.
!
    character(*),intent(in)::bfile
    character(*),intent(in)::expect  ! 'none' or 'grid'
    integer,intent(out)::south,north,psouth,pnorth,nsol,nassim,seqno
!
    character(4)::ftype
    integer::info
!
    open(unit=INNOSWAP,form='unformatted',status='old',file=bfile,iostat=info)
    if(info.ne.0) call fexit(trim(bfile)//': cannot open',errcode=2)
!
!  If the file is to be the whole globe, indicated by an expected header
!  string of 'none', then the actual file contains no header record except
!  for the date in yyyymmddhh format.  Supply default values in this case.
!
    if(expect.eq.'none') then
       south=1
       north=LATG
       psouth=1
       pnorth=LATG
       nassim=NAVARS
       ftype='none'
       read(INNOSWAP) seqno,nsol
    else
       read(INNOSWAP) ftype,seqno,south,north,psouth,pnorth,nsol,nassim
       if(nassim.ne.NAVARS) then
          write(IERR,900) ftype,seqno
          write(IERR,901) south,north,psouth,pnorth
          write(IERR,902) nsol,nassim
          call fexit('read_grid_header: nassim must be 4 in this version')
       endif
    endif
900 format('read_grid_header: bad header'/'ftype: ',a,' seqno: ',i0)
901 format('read_grid_header: south: ',i0,' north: ',i0,' psouth: ',i0, &
       ' pnorth: ',i0)
902 format('read_grid_header: nsol: ',i0,' nassim: ',i0)
!
!  Sanity checks on the header
!
    if(ftype.ne.expect) then
       write(IERR,100) trim(bfile),ftype,expect
100    format(a,': header is: ',a,' expected: ',a)
       call fexit(trim(bfile)//': wrong grid file header',2)
    else if(seqno.ne.seqid()) then
       write(IERR,101) trim(bfile),seqno,seqid()
101    format(a,': sequence id is: ',i0,' expected: ',i0)
       call fexit(trim(bfile)//': unexpected sequence number',2)
    else if(south.lt.1.or.north.gt.LATG.or.south.gt.north.or.psouth.lt.south &
      .or.pnorth.gt.north) then
       write(IERR,110) trim(bfile),south,north,psouth,pnorth
       call fexit(trim(bfile)//': start/stop indices out of range',2)
    endif
110 format(a,': nonsensical values for south=',i0,' north=',i0, &
       ' psouth=',i0,' pnorth=',i0)
    return
    end subroutine read_grid_header
! ---------------------
    subroutine read_model_grid(g,south,north,nsol)
!  READ_MODEL_GRID - read a section of the background grid from a previously
!  split file (and do some basic sanity checks).
!  Exit if an error occurs.
!
!  Arguments:
!  SOUTH : Southernmost latitude expected in the grid file.
!  NORTH : Northernmost latitude index expected in the grid file.
!  NSOL : number of ensemble solutions expected in the grid file.
!  G  := Allocated grid slice containing the grid data from the file.
!     The latitude dimension goes from SOUTH to NORTH.
!  Side effects:
!  - Must be preceded by a call to READ_GRID_HEADER.
!  - Closes INNOSWAP.
!
    integer,intent(in)::south,north,nsol
    real(GRIDKIND),intent(out)::g(LONF,NAVARS,south:north,LEVS,nsol)
!
    real(IOKIND),allocatable::buf(:,:,:,:,:)
!
    if(IOKIND.ne.GRIDKIND) then
       allocate(buf(LONF,NAVARS,south:north,LEVS,nsol))
       read(INNOSWAP) buf
       g=buf
    else
       read(INNOSWAP) g
    endif
    close(INNOSWAP)
    return
    end subroutine read_model_grid
! -------------------
    subroutine write_grid_slice(gfile,header,south,north,g,lower,upper,nsol, &
       psouth,pnorth)
! WRITE_GRID_SLICE - Write a section of the global background grid
!  to the indicated output file in latitude-major format.
!  Halt the program if an error occurs.
!  GFILE   : name of output grid file
!  HEADER  : 4-character header string to use.  If set to 'none', we
!            write no header.
!  NORTH, SOUTH :   Define the boundaries of the latitudinal strip to write out.
!  ISTART  : index at which data assimilation should begin for the strip.
!    Assimilation occurs at patches centered at latitudes ISTART, ISTART+1,...
!    ISTART+YWIN-1, where YWIN is obtained from PATCH_RADIUS.
!  GRID  :=: grid of model data to output.  On return, contains converted/scaled
!    data in latitudes SOUTH to NORTH.
!  LOWER, UPPER : latitude dimensions of GRID.
!  NB      : the number of ensemble solutions to output.
!  PSOUTH, PNORTH : if present, these define the latitude strip over which
!    patches may be centered for the data assimilation (mainly useful in
!    preanalysis).  If absent, PSOUTH and PNORTH coincide with SOUTH and NORTH.
!    They are meaningful only if HEADER isn't 'NONE'.
!  Side effects:
!  - Opens and closes OUTNOSWAP.
!
    character(*),intent(in)::gfile
    character(*),intent(in)::header
    integer,intent(in)::south,north,lower,upper,nsol
    real(GRIDKIND),intent(in)::g(LONF,NAVARS,lower:upper,LEVS,nsol)
    integer,intent(in),optional::psouth,pnorth
!
!  Local variables
!
    integer::info,asouth,anorth,id
    character(4)::ftype
    real(IOKIND),allocatable::buf(:,:,:,:,:)
!
!  Sanity checks
!
    if(south.lt.lower.or.north.gt.upper.or.south.gt.north.or.nsol.lt.1.or. &
       min(south,lower).lt.1.or.max(north,upper).gt.LATG) &
          call fabort('write_grid_slice: bad arguments')
!
    open(unit=OUTNOSWAP,file=gfile,form='unformatted',action='write', &
       position='rewind',iostat=info)  ! overwrite any predecessor
    if(info.ne.0) call fexit(trim(gfile)//': cannot open',errcode=2)
!
!  Write header
!
    ftype=header
    if(ftype.ne.'none') then
       if(present(psouth)) then
          asouth=psouth
       else
          asouth=south
       endif
       if(present(pnorth)) then
          anorth=pnorth
       else
          anorth=north
       endif
       id=seqid()
       if(verbosity().gt.VERB_TERSE) &
          write(IERR,100) trim(gfile),south,north,nsol,asouth,anorth
100    format('write_grid_slice: ',a/'   south: ',i3,' north: ',i3,' nsol: ', &
          i3,' asouth: ',i3,' anorth: ',i3)
       write(OUTNOSWAP) ftype,id,south,north,asouth,anorth,nsol,NAVARS
    else
       write(OUTNOSWAP) seqid(),nsol
    endif
!
    if(IOKIND.ne.GRIDKIND) then
       allocate(buf(LONF,NAVARS,lower:upper,LEVS,nsol))
       buf=g
       write(OUTNOSWAP) buf(:,:,south:north,:,:)
    else
       write(OUTNOSWAP) g(:,:,south:north,:,:)
    endif
    close(OUTNOSWAP)
    return
    end subroutine write_grid_slice
!-------------------------------------
    subroutine output_analysis(agrid,lower,upper,nsol,afile,header)
!  OUTPUT_ANALYSIS - output the given grid.
!  AGRID  : array holding the analysis at each grid point.
!  LOWER, UPPER  : southern- and northernmost latitudes in the analysis.
!  AFILE  : must be nonblank and contain the name of the file to which
!            we output the analysis.
!  NSOL  : number of ensemble solutions to output
!  HEADER  : 4-character header to output
!  Side effects:
!  - opens and closes OUTNOSWAP (via write_grid_slice)
!
    integer,intent(in)::lower,upper,nsol
    real(GRIDKIND),intent(in)::agrid(nsol,lower:upper,LEVS,LONF,NAVARS)
    character(*),intent(in)::afile,header
!
    integer::south,north,j,k,lev
    character(4)::hdr
    real(GRIDKIND),allocatable::transgrid(:,:,:,:,:)
!
    hdr=header
    call assimilation_boundary(south,north)
    if(lower.gt.south.or.upper.lt.north) &
       call fabort('output_analysis: analysis array too small')
    if(south.eq.1.and.north.ge.LATG) then
       north=LATG
       hdr='none'
    else
       hdr=header
    endif
!
!  Copy contents of analysis grid to a grid whose dimensions match
!  what's expected for output.  This code isn't cache-friendly.
!  If it proves to be a problem then we can perform an intermediate
!  transpose to an array of shape (LONF,NAVARS,nsol,lower:upper,LEVS).
!
    allocate(transgrid(LONF,NAVARS,lower:upper,LEVS,nsol))
    do j=1,nsol
       do k=lower,upper
          do lev=1,LEVS
             transgrid(:,:,k,lev,j)=agrid(j,k,lev,:,:)
          enddo
       enddo
    enddo
    call write_grid_slice(afile,hdr,south,north,transgrid,lower,upper,nsol)
    return
    end subroutine output_analysis
! ---------------------------------------------------
    end module grid_io
