    module grid_model
!
!  SYNOPSIS
!    use grid_model
!
!  DESCRIPTION
!    GRID_MODEL contains support routines for copying variables to and from
!    model grid arrays, whether the global ensemble grid, the analysis grid,
!    or the nature run grid.
!    In addition, GRID_MODEL contains the initial implementation of the H
!    operator, which for now calls only a simple linear interpolation routine
!    to interpolate the model state to each observation.  In the future,
!    additional routines will be implemented for nonlinear interpolations.
!
!  ROUTINES IMPLEMENTED
!   subroutine copy_grid_point(loc,x,typ,what)
!   subroutine copy_from_ensemble(loc,x,nsol,patchsz,typelist)
!   subroutine copy_from_mean(loc,x,patchsz,typelist)
!   subroutine copy_grid_z(loc,iz,x,typ,what)
!   private subroutine copy_from(loc,x,patchsz,grid,south,north,nsol,typelist)
!   private subroutine copy_to(loc,nc,agrid,south,north,nsol,center)
!   subroutine model_grid_init(bfile,afile,mfile,dfile,expect)
!   subroutine model_grid_fini(have_model_grid)
!   subroutine transpose_grid(inbuf,outbuf,south,north,nsol)
!   subroutine precompute_ensemble_differences
!
!  REQUIRED PACKAGES
!  assim_vars, grid_io,  interpolate - in main source directory
!  ferror - in support directory
!
!  MODIFICATION HISTORY
!  10/13/05 - changed COPY_TO to account for change in wind direction
!     when going over the poles.
!  09/22/05 - added COPY_GRID_Z.
!  08/06/05 - added GRIDKIND to manage internal types of model grids.
!  06/19/05 - removed diagnostics; added support for nature grid.
!  06/10/05 - incorporated diagnostics from previous small and seaparate
!    module.
!  01/19/05 - changes to accommodate real observations (finally).
!     Added code to get the sigma level (relative to model grid) from the
!     background mean surface pressure and the P_LEV observation value.
!  10/11/04 - further changes to accommodate real observations.
!     Installed new i/o format for the model grid.  First version of the
!     H operator (linear interpolation from nearly model grid points).
!  09/27/04 - intermediate version to prepare for actual observations.
!     Only surface pressure, virtual temperature, and u- and v-wind components
!     are assimilated in this version.  The grid indexing is changed
!     for greater cache coherency.   The model_vars type has been eliminated
!     in favor of a single large array.
!  AUTHOR
!     Eric Kostelich, kostelich@asu.edu, 480-965-5006.
!
    use grid_io
    use command_line_support,only:program_exit
    implicit none
!
!  Flags used in check_model_grid
!
    integer,parameter::NONE=0,BAD_VALUE=1,INF_OR_NAN=2
!
!  Private module data
!  -------------------
!  BACKGROUND defines the model background (forecast).
!  BMEAN contains the ensemble mean.
!  ANALYSIS contains the analysis.
!  AMEAN contains the analysis mean.
!  When the ensemble differences are precomputed, BACKGROUND eventually
!  holds the differene between the actual model background and the
!  background mean (an optimization to avoid recomputing them on each
!  patch when doing the data assimilation).
!  ANALYSIS and AMEAN are smaller grids than the
!  background and the background mean, because they omits the padding
!  necessary to center patches at every boundary point.
!
!  NOTE:  All model variables are scaled to energy units!
!  See to_energy_units and from_energy_units in spectral.f90.
!
    real(GRIDKIND),allocatable,save,private::background(:,:,:,:,:) 
                                         ! (nsol,south:north,LEVS,LONF,NVAR)
    real(GRIDKIND),allocatable,save,private::bmean(:,:,:,:)
                                         ! (south:north,LEVS,LONF,var)
    real(GRIDKIND),allocatable,save,private::analysis(:,:,:,:,:)
                                         ! (nsol,psouth:pnorth,LEVS,LONF,NVAR)
    real(GRIDKIND),allocatable,save,private::amean(:,:,:,:)
                                         ! (psouth:pnorth,LEVS,LONF,var)
    logical,save,private::did_diffs=.false.  ! true if ensemble diffs are
                                             ! precomputed
!
!  The "true" model variables from the nature run, if available;
!  natgrid is not associated otherwise.
!
    real(GRIDKIND),allocatable,save,private::natgrid(:,:,:,:) 
                                           ! (south:north,LEVS,LONF,NVAR)
!
!  Names of the output files containing the analysis and analysis mean
!
    character(MAX_FILE_NAME_LEN),save,private::analysis_filename=' '
    character(MAX_FILE_NAME_LEN),save,private::amean_filename=' '
!
!  Internal work routines
!
    private::copy_from,copy_to
    contains
!----------------------------------------------------------------------------
    subroutine model_grid_init(bfile,afile,mfile,dfile,expect)
!  MODEL_GRID_INIT - initialize the model grid by reading in the physical
!  space data output by preanalysis and ecomputing the ensemble means.
!  BFILE : background grid file
!  AFILE : output analysis grid file
!  MFILE : output analysis mean file
!  DFILE : output diagnostics file
!  EXPECT : header to expect in background grid file
!
    character(*),intent(in)::bfile,afile,expect,mfile,dfile
!
!  Local variables
!
    integer::south,north,psouth,pnorth,nb,nsol,nassim,seqno,info,k, &
       errtype,var
    type(patch_index)::loc
    character(MAX_FILE_NAME_LEN+80)::msg
    real(GRIDKIND),allocatable::buf(:,:,:,:,:),avg(:,:,:,:)
    real(D)::badval
!
    call read_grid_header(bfile,expect,south,north,psouth,pnorth, &
       nb,nassim,seqno)
    call set_model_grid_boundary(south,north,psouth,pnorth)
!
!  Sanity checks
!
    nsol=ensemble_size()
    if(nb.ne.nsol) then
       write(msg,110) trim(bfile),nb,ensemble_size()
110    format(a,': ensemble size is: ',i0,' expected: ',i0)
       call fexit(trim(msg),2)
    else if(seqno.ne.seqid()) then
       write(IERR,111) trim(bfile),seqno,seqid()
111    format(a,' has sequence number: ',i0,' expected: ',i0)
       call fexit(trim(msg),2)
    else if(nassim.ne.NAVARS) then
       write(msg,112) NAVARS,nassim
       call fexit(trim(msg),2)
112    format('model_grid_init: expected ',i3,' variables; got ',i3)
    endif
!
!  Read in the ensemble solutions and initialize the global grid.
!
    allocate(background(nsol,south:north,LEVS,LONF,NAVARS),stat=info)
    if(info.ne.0) then
       write(msg,120) 'background',south,north
       call fexit(trim(msg),2)
    endif
120 format('model_grid_init: cannot allocate ',a,' latitudes: ',i3,' to ',i3)
!
!  Temporary input buffer.  Note: we can modify the structure of
!  READ_MODEL_GRID to read one solution per call; at the risk of slightly
!  more complicated open/close file logic, we would only need BUF to
!  be large enough to hold one solution.  This could become a factor
!  when the ensemble sizes get big enough.
!
    allocate(buf(LONF,NAVARS,south:north,LEVS,nsol),stat=info)
    if(info.ne.0) then
       write(msg,120) 'buf',south,north
       call fexit(trim(msg),2)
    endif
    call read_model_grid(buf,south,north,nsol)
!
!  Convert the input buffer to the subscript order that's most efficient
!  for the rest of the data assimilation algorithm.
!
    call transpose_grid(buf,background,south,north,nsol)
!
!  Do perfunctory checks on the data to try to identify definite garbage values.
!
    call check_model_grid(background,nsol,south,north,errtype,loc,var,badval)
    if(errtype.ne.0) then
       if(errtype.eq.INF_OR_NAN) then
          write(IERR,150) trim(bfile),'inf or nan',TYPEOF(var),loc
       else
          write(IERR,150,advance='no') trim(bfile),'bad value',TYPEOF(var),loc
          write(IERR,151) badval
       endif
150    format(a,': '/a,': variable ',a1,' at lon:',i4,' lat:',i4,' lev:',i3)
151    format(' value: ',es13.4)
       call fexit('read_model_grid: bad value/infinity detected on input',2)
    endif
!
!  Get ensemble means (needed by the H operator).  Use a cache-friendly
!  algorithm, then convert to the subscript storage order needed elsewhere.
!
    allocate(avg(LONF,NAVARS,south:north,LEVS),stat=info)
    if(info.ne.0) then
       write(msg,120) 'avg',south,north
       call fexit(trim(msg),2)
    endif
    avg=buf(:,:,:,:,1)
    do k=2,nsol
       avg=avg+buf(:,:,:,:,k)
    enddo
    avg=avg/nsol
    deallocate(buf)
!
    allocate(bmean(south:north,LEVS,LONF,NAVARS),stat=info)
    if(info.ne.0) then
       write(msg,120) 'bmean',south,north
       call fexit(trim(msg),2)
    endif
    call transpose_grid(avg,bmean,south,north,1)
    deallocate(avg)
!
!  Store the names of the analysis output files.  We must allocate space
!  the analysis mean if diagnostics are requested, because we need the
!  analysis grid and mean to compare with the observations.  And of course
!  we need to allocate the analysis mean if its computation is requested.
!  Eric: potential optimization: copy background mean to analysis mean and
!  the background ensemble to the analysis ensemble; then return_background
!  isn't necessary when nobs=0 in the assimilation step.
!
    analysis_filename=afile
    if(len_trim(mfile).gt.0.or.len_trim(dfile).gt.0) then
       amean_filename=mfile
       allocate(amean(psouth:pnorth,LEVS,LONF,NAVARS),stat=info)
       if(info.ne.0) then
          write(msg,120) 'analysis mean',psouth,pnorth
          call fexit(trim(msg),2)
       endif
    endif
!
!  The analysis is a narrower latitude strip, in general, because we
!  don't need the extra padding to center patches.
!
    allocate(analysis(nsol,psouth:pnorth,LEVS,LONF,NAVARS),stat=info)
    if(info.ne.0) then
       write(msg,120) 'analysis',psouth,pnorth
       call fexit(trim(msg),2)
    endif
    return
    end subroutine model_grid_init
!----------------------------------------------------------------------------
    subroutine model_grid_fini
!  MODEL_GRID_FINI - output the analysis and optionally the analysis mean
!  and sum of squares (so we can compute global RMS error).
!
    integer::asouth,anorth,nsol
!
    nsol=ensemble_size()
    call assimilation_boundary(asouth,anorth)
    call output_analysis(analysis,asouth,anorth,nsol,analysis_filename,'aout')
    if(len_trim(amean_filename).gt.0.and.allocated(amean)) &
       call output_analysis(amean,asouth,anorth,1,amean_filename,'mean')
    return
    end subroutine model_grid_fini
!----------------------------------------------------------------------------
    subroutine copy_grid_point(loc,x,typ,what)
!  COPY_GRID_POINT - copy the ensemble solutions for the given dynamical
!  variable at the point on the model grid specified by LOC.
!  If LOC has invalid entries, or if LOC specifies a location that is
!  not available on the latitude strip on which we are working, then
!  the routine aborts with an error message.
!  LOC  : specifies the model grid location
!  X   := gets as many ensemble solution values as will fit.
!  TYP  : the dynamical variable type to obtain (one of PRES_TYPE, etc.)
!  WHAT : one of 'background', 'bmean', 'analysis', 'amean', or 'nature'.
!     (Only the first two letters are significant.  MUST be in lower case.)
!
    type(patch_index),intent(in)::loc
    integer,intent(in)::typ
    real(D),intent(out)::x(:)
    character(*),intent(in)::what
!
    integer::n
    character(80)::msg
!
    if(loc%x.lt.1.or.loc%x.gt.LONF.or.loc%z.lt.1.or.loc%z.gt.LEVS) &
       call fabort('copy_grid_point: bad index')
!
!  Ensure that LOC%Y corresponds to a latitude contained in the appropriate
!  grid.
!
    if(valid_dynamical_type(typ)) then
       select case(what(1:2))
       case('ba')
          if(loc%y.lt.lbound(background,2).or.loc%y.gt.ubound(background,2)) &
             goto 911
          n=min(size(x),ensemble_size())
          x(1:n)=background(1:n,loc%y,loc%z,loc%x,typ)
       case('bm')
          if(loc%y.lt.lbound(bmean,1).or.loc%y.gt.ubound(bmean,1)) goto 911
          x=bmean(loc%y,loc%z,loc%x,typ)
       case('an')
          if(loc%y.lt.lbound(analysis,2).or.loc%y.gt.ubound(analysis,2)) &
             goto 911
          n=min(size(x),ensemble_size())
          x(1:n)=analysis(1:n,loc%y,loc%z,loc%x,typ)
       case('am')
          if(loc%y.lt.lbound(amean,1).or.loc%y.gt.ubound(amean,1)) goto 911
          x=amean(loc%y,loc%z,loc%x,typ)
       case('na')
          if(.not.allocated(natgrid)) &
             call fabort('copy_grid_point: no nature grid')
          if(loc%y.lt.lbound(natgrid,1).or.loc%y.gt.ubound(natgrid,1)) goto 911
          x=natgrid(loc%y,loc%z,loc%x,typ)
       case default
          call fabort('copy_grid_point: unknown source')
       endselect
    else
       call fabort('copy_grid_point: bad type')
    endif
    return
!
!  Oops
!
911 continue
    write(msg,100) what,loc
100 format('copy_grid_point: bad index for ',a,': loc: ',3i4)
    call fabort(msg)
    end subroutine copy_grid_point
!----------------------------------------------------------------------------
    subroutine copy_grid_z(loc,iz,x,typ,what)
!  COPY_GRID_Z - like COPY_GRID_POINT, but the Jth value returned in X
!  is the dynamical value at (LOC%X,LOC%Y,IZ(J)).
!  LOC  : specifies the longitude and latitude of the model grid location
!  IZ  : vector of vertical positions to obtain.
!  X  := Jth entry is dynamical variable at (LOC%X,LOC%Y,IZ(J)).  Must be
!     same size as IZ.
!  TYP  : the dynamical variable type to obtain (one of PRES_TYPE, etc.)
!  WHAT : one of 'background', 'bmean', 'analysis', 'amean', or 'nature'.
!     (Only the first two letters are significant.  MUST be in lower case.)
!
    type(patch_index),intent(in)::loc
    integer,intent(in)::iz(:),typ
    real(D),intent(out)::x(:)
    character(*),intent(in)::what
!
    integer::j
    character(80)::msg
!
    if(loc%x.lt.1.or.loc%x.gt.LONF.or.any(iz.lt.1).or.any(iz.gt.LEVS)) &
       call fabort('copy_grid_z: bad index')
    if(size(iz).ne.size(x)) call fabort('copy_grid_z: array size mismatch')
!
!  Ensure that LOC%Y corresponds to a latitude contained in the appropriate
!  grid.
!
    if(valid_dynamical_type(typ)) then
       select case(what(1:2))
       case('ba')
          if(loc%y.lt.lbound(background,2).or.loc%y.gt.ubound(background,2)) &
             goto 911
          do j=1,size(x)
             x(j)=background(j,loc%y,iz(j),loc%x,typ)
          enddo
       case('bm')
          if(loc%y.lt.lbound(bmean,1).or.loc%y.gt.ubound(bmean,1)) goto 911
          do j=1,size(x)
             x(j)=bmean(loc%y,iz(j),loc%x,typ)
          enddo
       case('an')
          if(loc%y.lt.lbound(analysis,2).or.loc%y.gt.ubound(analysis,2)) &
             goto 911
          do j=1,size(x)
             x(j)=analysis(j,loc%y,iz(j),loc%x,typ)
          enddo
       case('am')
          if(loc%y.lt.lbound(amean,1).or.loc%y.gt.ubound(amean,1)) goto 911
          do j=1,size(x)
             x(j)=amean(loc%y,iz(j),loc%x,typ)
          enddo
       case('na')
          if(.not.allocated(natgrid)) call fabort('copy_grid_z: no nature grid')
          if(loc%y.lt.lbound(natgrid,1).or.loc%y.gt.ubound(natgrid,1)) goto 911
          do j=1,size(x)
             x(j)=natgrid(loc%y,iz(j),loc%x,typ)
          enddo
       case default
          call fabort('copy_grid_z: unknown source')
       endselect
    else
       call fabort('copy_grid_z: bad type')
    endif
    return
!
!  Oops
!
911 continue
    write(msg,100) what,loc
100 format('copy_grid_z: bad index for ',a,': loc: ',3i4)
    call fabort(msg)
    end subroutine copy_grid_z
!----------------------------------------------------------------------------
    subroutine copy_from_ensemble(loc,x,nsol,patchsz,typelist)
!  COPY_FROM_ENSEMBLE - copy the relevant model dynamical variables from the
!  background ensemble to an array representing a patch whose logical
!  center is LOC.  If the ensemble mean has been subtracted off, then
!  the values returned are the _differences_ from the mean, not the
!  original background values.
!
!  Arguments:
!  LOC   : (longitude,latitude,z) coordinates of center of patch
!  X     := vector whose Kth row contains the PATCHSZ dynamical variables 
!     of the Kth ensemble solution for the patch centered at (IX,IY).
!  NSOL : leading dimension of X.  Must equal actual number of ensemble
!     solutions.
!  TYPELIST := integer array whose (J,K)th element is the variable type in
!    X(J,K).  Optional argument.
!
!  Variables are copied into X from bottom to top and left to right
!  (i.e., in order of increasing latitude, then increasing longitude).
!  X(1) contains the surface pressure; X(2:NLEV+1) holds the temperature
!  at the vertical levels in question, and so on.  At each point,
!  X holds 3*NLEV+1 dynamical quantities, where NLEV is the number of
!  vertical levels in the patch.  We include the surface pressure
!  in each analysis, even in columns that aren't near the surface.
!
    integer,intent(in)::nsol,patchsz
    type(patch_index),intent(in)::loc
    real(D),intent(out)::x(nsol,patchsz)
    integer,intent(out),optional::typelist(nsol,patchsz)
!
!  Local variables
!
    integer::south,north
!
    if(patchsz.lt.patchsize(loc).or.nsol.ne.ensemble_size()) &
       call fabort('copy_from_ensemble: array X is wrong size')
    call model_grid_boundary(south,north)
    call copy_from(loc,x,patchsz,background,south,north,nsol,typelist)
    return
    end subroutine copy_from_ensemble
!----------------------------------------------------------------------------
    subroutine copy_from_mean(loc,x,patchsz,typelist)
!  COPY_FROM_MEAN - return in X the mean solution over the patch whose
!  logical center is LOC.  The ordering of the variables must be consistent
!  with that in COPY_FROM_ENSEMBLE.
!
    integer,intent(in)::patchsz
    type(patch_index),intent(in)::loc
    real(D),intent(out)::x(patchsz,1)
    integer,intent(out),optional::typelist(patchsz)
!
    integer::south,north
!
    call model_grid_boundary(south,north)
    call copy_from(loc,x,patchsz,bmean,south,north,1,typelist)
    return
    end subroutine copy_from_mean
!----------------------------------------------------------------------------
    subroutine copy_from(loc,x,patchsz,grid,south,north,nsol,typelist)
!  COPY_FROM_RG - copy variables from the model grid array into 
!  a vector that represents the dynamical variables on a patch.
!  Internal work routine.  This routine defines the ordering of the 
!  model dynamical variables in a patch.
!  Note: we always copy the surface pressure into the patch, even if it
!  does not intersect the surface.
!
    integer,intent(in)::patchsz,south,north,nsol
    type(patch_index),intent(in)::loc
    real(GRIDKIND),intent(in)::grid(nsol,south:north,LEVS,LONF,NAVARS)
    real(D),intent(out)::x(nsol,patchsz)
    integer,intent(out),optional::typelist(nsol,patchsz)
!
!  Local variables
!
    integer::offset,j,k,zlo,zhi,xwin,ywin,nlev,ix,iy
    type(patch_index)::newloc
    logical::over_pole
!
    if(patchsz.lt.patchsize(loc)) &
       call fabort('copy_from_mean: array X is too small')
!
    call patch_radius(loc,xwin,ywin)
    do j=-xwin,xwin
       do k=-ywin,ywin
          call xy_adjust(loc,j,k,newloc,zlo,zhi,over_pole)
          offset=patch_offset(loc,j,k)
          nlev=zhi-zlo+1
          if(offset+3*nlev.gt.patchsz)  &
             call fabort('copy_from: patch offset out of range')
          ix=newloc%x
          iy=newloc%y
          x(:,offset)=grid(:,iy,1,ix,PRES_TYPE)
          if(present(typelist)) typelist(:,offset)=PRES_TYPE
          x(:,offset+1:offset+nlev)=grid(:,iy,zlo:zhi,ix,TEMP_TYPE)
          if(present(typelist)) typelist(:,offset+1:offset+nlev)=TEMP_TYPE
          offset=offset+nlev
!
!  An east wind at a point centered near the pole is blowing west relative to
!  a point centered opposite.
!
          if(over_pole) then
             x(:,offset+1:offset+nlev)=-grid(:,iy,zlo:zhi,ix,UWIND_TYPE)
          else
             x(:,offset+1:offset+nlev)=grid(:,iy,zlo:zhi,ix,UWIND_TYPE)
          endif
          if(present(typelist)) typelist(:,offset+1:offset+nlev)=UWIND_TYPE
          offset=offset+nlev
          if(over_pole) then
             x(:,offset+1:offset+nlev)=-grid(:,iy,zlo:zhi,ix,VWIND_TYPE)
          else
             x(:,offset+1:offset+nlev)=grid(:,iy,zlo:zhi,ix,VWIND_TYPE)
          endif
          if(present(typelist)) typelist(:,offset+1:offset+nlev)=VWIND_TYPE
          offset=offset+nlev
       enddo
    enddo
    return
    end subroutine copy_from
!----------------------------------------------------------------------------
    subroutine copy_to_analysis(loc,nc,nsol,center,error,var,badval)
!  COPY_TO_ANALYSIS - copy the analysis at the CENTER of the patch at
!  location LOC to the analysis grid.
!  NC  : Number of variables in the center Column of the patch
!  NSOL : ensemble size
!  CENTER : Kth column contains the NC variables for the Kth ensemble solution.
!  If either of ERROR or VAR is present, then check for nonsense values
!  and flag any that are found in ERROR and note the bad variable type in VAR.
!  ERROR (if present) is set to 0 unless a problem is found.
!
    type(patch_index),intent(in)::loc
    integer,intent(in)::nc,nsol
    real(D),intent(in)::center(nsol,nc)
    integer,intent(out),optional::error,var
    real(D),intent(out),optional::badval  ! if any
!
    integer::asouth,anorth,err,which
    real(D)::val
!
    call assimilation_boundary(asouth,anorth)
    call copy_to(loc,nc,analysis,asouth,anorth,nsol,center)
    if(present(error).or.present(var)) then
       call check_grid_at(analysis,nsol,asouth,anorth,loc,err,which,val)
       if(present(error)) error=err
       if(present(var)) var=which
       if(present(badval)) badval=val
    endif
    return
    end subroutine copy_to_analysis
!----------------------------------------------------------------------------
    subroutine copy_to_mean(loc,nc,cmean)
!  COPY_TO_MEAN - store the analysis mean at the center of the patch at LOC.
!  NC  : Number of variables in the center Column of the patch
!  CMEAN : Kth entry is the mean of the Kth row of CENTER.
!
    type(patch_index),intent(in)::loc
    integer,intent(in)::nc
    real(D),intent(in)::cmean(nc)
!
    integer::asouth,anorth
!
    if(allocated(amean)) then
       call assimilation_boundary(asouth,anorth)
       call copy_to(loc,nc,amean,asouth,anorth,1,cmean)
    endif
    return
    end subroutine copy_to_mean
!----------------------------------------------------------------------------
    subroutine copy_to(loc,nc,agrid,asouth,anorth,nsol,center)
!  COPY_TO - copy CENTER to AGRID, which can be either the full analysis
!  grid or the analysis mean.  Internal work routine.
!  The order in which variables are copied back must be compatible with
!  the ordering defined in COPY_FROM.
!  LOC=(IX,IY,IZ) to an analysis grid.  
!  NC : number of dynamical variables (including the surface pressure)
!     included in the vertical column at CENTER.
!  AGRID : destination grid from latitudes ASOUTH to ANORTH.
!  NSOL  : number of ensemble solutions to copy
!  CENTER : array containing the dynamical variables at the logical center
!     of the patch.
!  Note 1:  although the surface pressure is always included in the analysis,
!  we copy it to the analysis grid only for vertical columns
!  that are intended to give a surface analysis.
!  The ordering of the variables in CENTER depends on copy_from_grid.
!  Note 2:  any deliberate type conversions (single/double) are performed
!     with intrinsic REAL.
!
    integer,intent(in)::asouth,anorth,nsol,nc
    type(patch_index),intent(in)::loc
    real(GRIDKIND),intent(inout)::agrid(nsol,asouth:anorth,LEVS,LONF,NAVARS)
    real(D),intent(in)::center(nsol,nc)
!
!  Local variables
!
    integer::off,ix,iy,iz,zlow,zhigh,nlev
!
!  Sanity checks
!
    if(loc%y.lt.asouth.or.loc%y.gt.anorth)  &
       call fabort('copy_to: bad latitude index')
    if(loc%x.lt.1.or.loc%x.gt.LONF) call fabort('copy_to: bad longitude index')
!
!  When we copy the solution at the center, we compute not only
!  the (x,y) offset to obtain the center of the grid, but we also
!  copy only the dynamical variables at the logical center of the patch.
!  (The logical center differs from the actual
!  center at patches at the surface or at the top of the atmosphere.)
!  When we have vertical localization, the logical center is a single
!  point; otherwise, it is the entire vertical column.
!
    call zdiam(loc,zlow,zhigh)
    nlev=zhigh-zlow+1
    ix=loc%x
    iy=loc%y
    if(nlev.eq.LEVS) then   ! no vertical localization
       do iz=1,LEVS
          agrid(:,iy,iz,ix,PRES_TYPE)=real(center(:,1),kind(agrid))
       enddo
       off=2
       do iz=1,LEVS
          agrid(:,iy,iz,ix,TEMP_TYPE)=real(center(:,off+iz-1),kind(agrid))
       enddo
       off=off+LEVS
       do iz=1,LEVS
          agrid(:,iy,iz,ix,UWIND_TYPE)=real(center(:,off+iz-1),kind(agrid))
       enddo
       off=off+LEVS
       do iz=1,LEVS
          agrid(:,iy,iz,ix,VWIND_TYPE)=real(center(:,off+iz-1),kind(agrid))
       enddo
       off=off+LEVS
       if(off.ne.nc+1) then  ! OFF should point exactly 1 above the top
          call fabort('copy_to: center indexing botched')
       endif
    else
       if(loc%z.eq.1) then ! include surface pressure
          do iz=1,LEVS
             agrid(:,iy,iz,ix,PRES_TYPE)=real(center(:,1),kind(agrid))
          enddo
       endif
       iz=loc%z
       off=iz-zlow+2  ! point to temperature at level iz
       agrid(:,iy,iz,ix,TEMP_TYPE)=real(center(:,off),kind(agrid))
       off=off+nlev
       agrid(:,iy,iz,ix,UWIND_TYPE)=real(center(:,off),kind(agrid))
       off=off+nlev
       agrid(:,iy,iz,ix,VWIND_TYPE)=real(center(:,off),kind(agrid))
       off=off+zhigh-iz
       if(off.gt.nc) then
          call fabort('copy_to: center indexing botched')
       endif
    endif
    return
    end subroutine copy_to
!----------------------------------------------------------------------------
    subroutine transpose_grid(inbuf,outbuf,south,north,nsol)
!  TRANSPOSE_GRID - Convert the input buffer from the subscript order that's
!  most efficient for spectral transforms to the order that's most efficient
!  for the rest of the data assimilation algorithm.
!  This part of the code is not cache-friendly.  If it becomes
!  a problem, we can do an intermediate transposition step to reduce the
!  number of cache misses before performing the scaling.
!
    integer,intent(in)::south,north,nsol
    real(GRIDKIND),intent(in)::inbuf(LONF,NAVARS,south:north,LEVS,nsol)
    real(GRIDKIND),intent(out)::outbuf(nsol,south:north,LEVS,LONF,NAVARS)
!
    integer::k,lat,lev
!
    do k=1,nsol
      do lat=south,north
         do lev=1,LEVS
             outbuf(k,lat,lev,:,:)=inbuf(:,:,lat,lev,k)
         enddo
       enddo
    enddo
    return
    end subroutine transpose_grid
!----------------------------------------------------------------------------
    subroutine precompute_ensemble_differences
!  PRECOMPUTE_ENSEMBLE_DIFFERENCES -
!  Precompute the ensemble differences: we compute BMEAN, whose
!  (Y,Z,X,V) entry is the ensemble mean of variable V at (X,Y,Z),
!  and subtract from BACKGROUND, whose (K,Y,Z,X,V)th entry becomes the 
!  difference between the Kth ensemble solution and the corresponding mean.
!
    integer::k
!
    if(did_diffs) return   ! don't do this twice
    do k=1,size(background,1)
       background(k,:,:,:,:)=background(k,:,:,:,:)-bmean(:,:,:,:)
    enddo
    did_diffs=.true.
    write(LOGFILE,'(a)') 'ensemble differences precomputed'
    return
    end subroutine precompute_ensemble_differences
!----------------------------------------------------------------------------
    logical function ensemble_diffs_precomputed(grid)
!  ENSEMBLE_DIFFS_PRECOMPUTED - return true or false according as
!  the ensemble differences have been precomputed for the indicated grid.
!  If GRID is omitted, it defaults to the background.
!
    character(*),intent(in),optional::grid
!
    if(present(grid)) then
       ensemble_diffs_precomputed=grid.eq.'background'.and.did_diffs
    else
       ensemble_diffs_precomputed=did_diffs
    endif
    return
    end function ensemble_diffs_precomputed
!----------------------------------------------------------------------------
    pure subroutine check_grid_at(grid,nsol,south,north,loc,errtype,var,badval)
!  CHECK_GRID_AT - perform heuristic checks for obvious problems with
!  the dynamical variables in GRID at LOC.  On return, ERROR is set if
!  VAR at LOC fails the test.
!  The checks are as follows:
!  1) NaN (not-a-number) value
!  2) surface pressure <= 0
!  3) temperature <= 0 
!  4) wind speed >= 500 m/s
!  ERRTYPE is set to 0 for no error, 1 for bad value, 2 for inf/NaN.
!  BADVAL is set to the first bad value (if any) that is found, provided that
!     ERRTYPE is BAD_VALUE.  The output value is undefined if ERRTYPE
!     is inf/NaN or if ERRTYPE is 0.
!  NOTE: This routine is to be called AFTER any internal scalings are applied!
!
    use inf_and_nan
    integer,intent(in)::nsol,south,north
    real(GRIDKIND),intent(in)::grid(nsol,south:north,LEVS,LONF,NAVARS)
    type(patch_index),intent(in)::loc
    integer,intent(out)::errtype
    integer,intent(out)::var
    real(D),intent(out)::badval
!
    real(D)::val
    integer::j
!
!  In energy units, values larger than this are likely trouble
!
    real(D),parameter::WINDMAX=10000.0 
!
!  Check for ridiculous values of surface pressure if the patch includes
!  the surface
!
    errtype=NONE
    var=PRES_TYPE
    if(any(isnan(grid(:,loc%y,loc%z,loc%x,PRES_TYPE)))) then
       errtype=INF_OR_NAN
    else if(loc%z.eq.1) then
       do j=1,nsol
          val=grid(j,loc%y,loc%z,loc%x,PRES_TYPE)
          if(val.le.0.0) then  ! in energy units, ln(dPa)*scale < 0 is trouble
             errtype=BAD_VALUE
             badval=val
             exit
          endif
       enddo
    endif
    if(errtype.ne.NONE) return
!
!  Check temperature
!
    errtype=NONE
    var=TEMP_TYPE
    if(any(isnan(grid(:,loc%y,loc%z,loc%x,TEMP_TYPE)))) then
       errtype=INF_OR_NAN
    else 
       do j=1,nsol
          val=grid(j,loc%y,loc%z,loc%x,TEMP_TYPE)
          if(val.le.0.0) then
             errtype=BAD_VALUE
             badval=val
             exit
          endif
       enddo
    endif
    if(errtype.ne.NONE) return
!
!  Check wind
!
    errtype=NONE
    var=UWIND_TYPE
    if(any(isnan(grid(:,loc%y,loc%z,loc%x,UWIND_TYPE)))) then
       errtype=INF_OR_NAN
    else
       do j=1,nsol
          val=abs(grid(j,loc%y,loc%z,loc%x,UWIND_TYPE))
          if(val.ge.WINDMAX) then
             errtype=BAD_VALUE
             badval=val
             exit
          endif
       enddo
    endif
    if(errtype.ne.NONE) return
!
    errtype=NONE
    var=VWIND_TYPE
    if(any(isnan(grid(:,loc%y,loc%z,loc%x,VWIND_TYPE)))) then
       errtype=INF_OR_NAN
    else
       do j=1,nsol
          val=abs(grid(j,loc%y,loc%z,loc%x,VWIND_TYPE))
          if(val.ge.WINDMAX) then
             errtype=BAD_VALUE
             badval=val
             exit
          endif
       enddo
    endif
    if(errtype.eq.NONE) var=NONE
    return
    end subroutine check_grid_at
!----------------------------------------------------------------------------
    pure subroutine check_model_grid(grid,nsol,south,north,errtype,loc,var,bad)
!  CHECK_MODEL_GRID - check whether the grid contains suspect entries.
!  If there are such entries, then toggle ERROR and
!  set LOC to the (J,K,LEV) index and variable type VAR
!  at which such an entry is found.
!
    use inf_and_nan
    integer,intent(in)::nsol,south,north
    real(GRIDKIND),intent(in)::grid(nsol,south:north,LEVS,LONF,NAVARS)
    integer,intent(out)::errtype
    type(patch_index),intent(out)::loc
    integer,intent(out)::var
    real(D),intent(out)::bad
!
    integer::ix,iy,iz
!
    do ix=1,LONF
       do iz=1,LEVS
          do iy=south,north
             loc=patch_index(ix,iy,iz)
             call check_grid_at(grid,nsol,south,north,loc,errtype,var,bad)
             if(errtype.ne.NONE) return
          enddo
       enddo
    enddo
    return
    end subroutine check_model_grid
! ----------------------------------------------------------------------
    subroutine read_nature_grid(header,truthfn,seq)
!  READ_NATURE_GRID - read the nature grid and scale its entries to be 
!  compatible with those used in the rest of the data assimilation code.
!  Side effects:
!   - allocates and defines natgrid
!
!   Arguments:
!   HEADER : character label to expect, if any, in the nature grid file
!   TRUTHFN : file name containing the nature grid
!
    character(*),intent(in)::header,truthfn
    integer,intent(in)::seq
!
    integer::psouth,pnorth,isouth,inorth,nassim,iens,seqno,south,north,info
    real(GRIDKIND),allocatable::buf(:,:,:,:,:)
!
    call read_grid_header(truthfn,header,isouth,inorth,psouth,pnorth,iens, &
       nassim,seqno)
    call model_grid_boundary(south,north)
    if(isouth.ne.south.or.inorth.ne.north) then
       write(IERR,900) isouth,south,inorth,north,psouth,pnorth
       call fexit(trim(truthfn)//': latitude boundary mismatch',2)
    endif
900 format('read_nature_grid: bad strip boundary: isouth: ',i0, &
       ' south: ',i0, ' inorth: ',i0,' north: ',i0,' psouth: ',i0, &
       ' pnorth: ',i0)
    if(iens.ne.1) then
       call fexit(trim(truthfn)//': contains more than one solution',2)
    else if(nassim.ne.NAVARS) then
       call fexit(trim(truthfn)//': number of model variables is wrong',2)
    else if(seqno.ne.seq) then
       call fexit(trim(truthfn)//': sequence number mismatch',2)
    endif
    allocate(buf(LONF,NAVARS,south:north,LEVS,1),stat=info)
    if(info.ne.0) then
       write(IERR,911) 'buf',info
       call program_exit(1)
    endif
911 format('diag_init: cannot allocate ',a,' info: ',i0)
    allocate(natgrid(south:north,LEVS,LONF,NAVARS),stat=info)
    if(info.ne.0) then
       write(IERR,911) 'natgrid',info
       call program_exit(1)
    endif
    call read_model_grid(buf,south,north,1)
    call transpose_grid(buf,natgrid,south,north,1)
    return
    end subroutine read_nature_grid
! ----------------------------------------------------------------------
    end module grid_model
