    module grid_obs
!  GRID_OBS - support for "real" observations.
!  
!  SYNOPSIS
!    use grid_obs
!
!  DESCRIPTION
!    GRID_OBS contains initial support for handling "real" observations.
!    The routines here take input files consisting of one raw_obs datum
!    per record and massage the values into a model_obs datum that contains
!    a subset of the information required for the data assimilation program
!    and also scales the observation values and standard errors accordingly.
!    In addition, GRID_OBS 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.
!
!  TYPES DEFINEED
!    type(raw_obs) - the "raw" observation format (after prepbufr preprocessing)
!    type(model_obs) - a condensed observation format used in the data
!      assimilation code
!    type(ensemble_stat) - statistics on the relevant ensemble solutions
!      interpolated to the observation location
!
!  ROUTINES IMPLEMENTED
!    subroutine obs_init(ofile,expect,nobs,nsol)
!    private subroutine retained(obs,nobs)
!    subroutine h_operator(nsol,nobs,kept)
!    private subroutine linear_interp(obs,n,hx,grid,ensstat)
!    subroutine apply_consistency_check(checktype)
!    subroutine set_observation_types(list,temp_units)
!    pure subroutine get_observation_types(list,num,temp_units)
!    subroutine set_qcm_threshold(q)
!    subroutine linear_interp(idx,nsol,hx,gridname,ensstat,dh)
!    subroutine gridpt_interp(x,ix,iadj,dx)
!    subroutine obs_energy_units(obs,nobs)
!    subroutine copy_from_obs(loc,yobs,hbar,rdiag,nsol,nobs,hx,dh)
!    private subroutine quality_control(obslist,es,nobs,kept,index)
!    subroutine validate(obslist,nobs,info)
!    private subroutine move_to_head(list,nobs,kept,index)
!    subroutine dump_obs(obs,msg,[optval],[raw])
!    pure subroutine kth_obs(k,obs,ierr,bg)
!    pure subroutine obs_count(kept,total)
!    subroutine read_obs_slice(ofile,expect,nobs,obsarray)
!    subroutine write_obs_slice(ofile,tag,south,north,obs,nobs)
!
!  REQUIRED PACKAGES
!    grid_model - in main source directory
!    spectral - in main source directory (for energy scalings)
!    interpolate - in main source directory (for linear interpolations)
!
!  MODIFICATION HISTORY
!  10/30/05 - Fixed "observation" quality control criterion to be the
!    maximum of 6.0 times the background or 6.0 times the observation 
!    standard error.
!  10/25/05 - Modified OBS_CONSISTENCY_CHECK to apply the quality control
!    only to the surface pressure for the time being.
!  10/10/05 - Updates to OBS_CONSISTENCY_CHECK: added the simple covariance
!    inflation (where applicable) to the background and added the "absdiff"
!    check.
!  10/08/05 - Fixes to interpolation algorithm in LINEAR_INTERP and SIGMA_FIXUP;
!     changed threshold in OBS_CONSISTENCY_CHECK from 5 to 6 times the 
!     approporiate standard error.
!  09/27/05 - Altered handling of vertical levels to vary by the interpolated
!     surface pressure of each ensemble member.  Rewrote SIGMA_FIXUP
!     and incorporated new energy unit conversion routines in spectral.f90.
!  08/23/05 - Fixed gridpt_interp to use XY_ADJUST to handle interpolations
!     correctly for patches over the poles.
!  08/07/05 - replaced KD tree with lookup that handles all 3 dimensions;
!     added code to include surface pressure observations in all patches.
!  07/25/05 - added artificial pressure sdev's for synthetic grid observations
!  07/18/05 - added DHX array.
!  06/29/05 - now converts observations to energy units.
!  06/15/05 - added support for interpolating values from a nature grid to
!     the location of each accepted observation.
!  06/02/05 - use z_grd, which is the sigma value computed from the interpolated
!     background surface pressure mean, instead of precomputed sigma from 
!     Gyorgyi's code, which assumed a fixed constant background surface 
!     pressure.
!  03/29/05 - modifications to accommodate new observation consistency checks.
!  02/04/05 - added checks to omit data with quality control marks above
!     a selected threshold (given by qcm_max).
!  01/21/05 - version to handle real data.  Added PHI, LAMBDA, and TIME_DIFF
!     components.  This version uses the background surface pressure to
!     obtain the sigma level of each datum.  However, for the moment,
!     we still use Gyorgyi's precomputed sigma levels.  Eventually,
!     changes will be required to obs_init to substitute the new z_grd
!     variable, which is computed from p_lev and the interpolated
!     background mean of the surface pressure.
!     Also added routines to select the observation types that are assimilated.
!  10/11/04 - initial version.  We assume that all observations are one of
!     the four model dynamical variables: surface pressure, virtual temperature,
!     u- or v-wind components.  Only simple linear interpolation is performed
!     from the model grid space to the observation space.
!
!  AUTHOR
!     Eric Kostelich, kostelich@asu.edu, 480-965-5006.
!
    use grid_model
    use spectral
    use interpolate
    implicit none
    integer,parameter::NUM_OBS_TYPES=4  ! number of different obs types
!-----------
!  Observation data types:
!   raw_obs is the "raw" data format derived from NCEP and other observation
!      files.
!   model_obs is the "processed" data format used in the data assimilation
!      step.  It contains the observation and its standard error in energy
!      units.
!   obs_stat contains statistical information for a given observation that
!      is used in the diagnostics package.
!-----------
!
!  ENSEMBLE_STAT is used to apply the observation consistency checks
!  and to compute the innovations in the Kalman filter.  Extra precision
!  is maintained here.
!
    type ensemble_stat
       real(D)::mean ! interpolated ensemble mean at observation location
       real(D)::sdev ! ensemble standard error (about interpolated mean)
    end type
!
!-----------
!  The RAW_OBS data type is used in the "real fake" datasets and in the
!  "actual" datasets as supplied by NCEP and filtered by Gyorgyi, with the
!  exceptions as noted before.  This observational information is passed
!  to preanalysis, which distributes the relevant subset to each processor 
!  for the data assimilation step.  The meaning of the components is as follows.
!
!  TYP         is an integer number assigned to all observations by NCEP. It
!              tells not only the type of the observed field
!              (e.g., temperature, wind, surface pressure, etc.)
!              but also the type of the instrument that took the
!              observation (e.g. radio-sonde, drop-sonde, aircraft,
!              etc.)
!
!  LEV         is equal to k, if this is the kth observation taken by
!              the same instrument (e.g., a radio-sonde)
!
!  X_grd       observational location in the x-direction in grid point units
!              (e.g. X_grd=15.8 means an observation between grid points
!              15 and 16, but four times closer to 16 than 15.  This field
!              is likely to be omitted in the future, as it can be computed
!              from LAMBDA.
!
!  Y_grd       observational location in the y-direction in grid point units.
!              This field is likely to be omitted in the future, as it can be
!              computed from PHI.
!
!  PHI         The latitude of the observation in degrees (-90 = south pole,
!              90 = north pole)
!
!  LAMBDA      The longitude of the observation in degrees;
!              0 <= LAMBDA < 360.  (0 = Greenwich, 180 = international dateline)
!
!  TIME_DIFF   difference between the actual observation time and the
!              date of the observation file in hours (e.g., when the 00Z
!              observation file is used and the data was taken at 23:30Z,
!              the value of TIME_DIFF is -0.5)
!
!  SIGMA   observational location in the z-direction (vertical
!          direction) in sigma units. Sigma is computed as the ratio
!          between the pressure at the observational location and
!          the surface pressure in the reference field (this field is
!          defined by the background mean in most cases).
!          For the moment, we use this value, instead of Z_grd, in the code;
!          Gyorgyi has precomputed it assuming a background surface pressure
!          of 1010 hPa.  Eventually, this field will be removed.
!
!  P_LEV       The pressure level in hPa; used to compute the sigma level
!              and model vertical coordinate of the observation.
!
!  VALUE       the numerical value of the observation in appropriate units.
!
!  SDEV        standard deviation of the assumed observational error for
!              the given type of observation (dependent on the type of
!              the observed field and the observing instrument)
!
!  QCM         quality control mark, the larger the value the more
!              likely that the observation has quality problem
!
!  CAT         extra information about the observation:
!              0 - surface level (mass report only)
!              1 - mandatory level (upper-air profile reports,
!                  e.g., radiosondes)
!              2 - significant temperature level (upper-air profile reports)
!              3 - winds-by-pressure level (upper-air profile reports)
!              4 - winds-by-height level (uperr-air profile reports)
!              5 - tropopause level (uperr-air profile reports)
!              6 - reports on a single level (e.g., aircraft, satellite wind,
!                  surface wind, precipitable water retrievals, etc.)
!                  [default value]
!              7 - Auxiliary levels generated via interpolation from
!                  spanning levels (upper-air profile reports)
!
!  PRSS        surface pressure in the reference field; this value is
!              used in the calculation of SIGMA
!  GRP         observation group (it may be useful when testing the
!              effects of different observation groups; e.g observations
!              by commercial aircraft)
!              ADPUPA - upper air (RAOB, PIBAL, RECCO, DROPS) reports
!                       RAOB - radio sonde
!                       PIBAL - drifting balloons for wind
!                       RECCO - reonnaissance mission
!                       DROPS - drop sondes
!              AIRCAR - Aicraft Commercial Addressing and Reporting
!                       System (ACARS) automated weather observations
!              AIRCFT - conventional (AIREP and PIREP) and ASDAR aircraft
!                       reports (all reported by the pilots)
!              SATWIND - satellite-derived wind reports
!              PROFLR - wind profiler reports
!              VADWND - VAD (NEXRAD) wind reports
!                       VAD - Velocity Azimuth Display (vertical profile
!                             of wind speed)
!                       NEXRAD - Doppler radar
!              SATEMP - NESDIS POES sounding, retrieval, radiance data
!                       (ATOVS)
!                       TOVS - Tiros N Operational Vertical Sounder
!              ADPSFC - surface land (SYNOPTIC, METAR) reports
!                       SYNOPTIC - surface observations at main synoptic
!                                   times (0Z, 06Z, 12Z, and 18Z)
!                       METAR - hourly surface land reports
!              SFCSHP - surface marine (SHIP, BOUY, C-MAN platform) reports
!              SFCBOG - mean sea-level pressure bogus reports
!              SPSSMI - DMSP SSM/I retrieval products (reprocessed wind
!                       speed, TPW)
!                       SSM - Special Sensor Microwave sounder
!                       (SSM/I is operated by the Air Force)
!              SYNDAT - synthetic tropical cyclone bogus reports
!              ERS1DA - ERS scatterometer data (reprocessed wind speed)
!                       ERS - Earth Remote Sensing (winds over oceans)
!              GOESND - NESDIS GOES sounding, retrieval, radiance data
!              QKSWND - QUICKSCAT scatterometer data (preprocessed wind
!                       speed; winds over ocean)
!
!  PRSS        surface pressure in the reference field; this value is not
!              used in this version
!
!  SID         id of the meteorological station that collected the
!              observation
!
!  FIELD       type of the observed meteorological field, the possible
!              values are:
!              P    - pressure [mb]=[hPa]
!              SP   - humidity [mg/kg]
!              T    - temperature [C]
!              H    - height [m]
!              U    - u component of wind [m/s]
!              V    - v component of wind [m/s]
!              WD   - wind direction [deg]
!              WS   - wind speed [knots]
!              TPW  - total precipitable water [mm]
!              RR   - rain rate [mm/h]
!              PW1  - precipitable water in the layer 1>sigma>0.9 [mm]
!              PW2  - precipitable water in the layer 0.9>sigma>0.7 [mm]
!              PW3  - precipitable water in the layer 0.7>sigma>0.3 [mm]
!              PW4  - precipitable water in the layer 0.3>sigma>0 [mm]
!              CTP  - cloud top pressure [mb]
!              CTT  - cloud top temperature [K]
!              TCC  - total cloud cover [%]
!
    integer,parameter::OBSKIND=R  ! real type of observations
    type raw_obs
       integer::typecode ! Code number for observation type
       integer::lev  ! serial # of obs. at given location
       real(R)::x_grd  ! observational location in x-direction
       real(R)::y_grd  ! observational location in y-direction
       real(R)::lambda  ! longitude grid location in radians
       real(R)::phi  ! latitude grid location in radians 
       real(R)::time_diff  ! drift time in hours
       real(R)::sigma  ! observational location in sigma
       real(OBSKIND)::p_lev ! pressure level in hPa
       real(OBSKIND)::value  ! observed value
       real(OBSKIND)::sdev  ! standard error of the observation
       real(OBSKIND)::qcm  ! quality control mark
       integer::cat  ! extra information about data
       real(OBSKIND)::prss   ! surface pressure
       character(8)::grp  ! observation group (e.g. ADPUPA)
       character(8)::sid  ! station id
       character(4)::field  ! observed meteorological field
    end type raw_obs
!
!  The MODEL_OBS type contains only the subset of the RAW_OBS data that
!  is most applicable to the data assimilation algorithm.  Prior to
!  the main assimilation step, each processor applies a reduction
!  of the components and a rescaling of selected observed values
!  and standard errors.  The LAMBDA and PHI components are copied from
!  those in RAW_OBS, but are not used except for listing in the diagnostics.
!
    type model_obs
       real(R)::x_grd(3)  ! longitude, latitude, and sigma level relative to
                          ! model grid index
       real(R)::lambda  ! longitude grid location in radians
       real(R)::phi  ! latitude grid location in radians 
       real(R)::time_diff ! observation time (relative to analysis time)
       real(OBSKIND)::p_lev  ! pressure level in hPa
       real(OBSKIND)::value  ! observed value, scaled to energy units
       real(OBSKIND)::sdev  ! standard deviation; scaled to energy units
       integer::typ  ! observation variable type (e.g., PRES_TYPE)
       integer::typecode ! code number for obs type (e.g., radiosondes)
       logical::kept  ! true if this obs is retained for the data assimilation
       !!real(D)::bgmean  ! duplicates BGSTAT entries for debugging
       !!real(D)::bgsdev  ! duplicates BGSTAT entries for debugging
    end type model_obs
!
!  Note: ZLEV and SURFPRES, below, could be allocatable arrays that are
!  part of the MODEL_OBS type.  I have elected to make them separate
!  arrays in part to avoid lots of separate heap allocations.
!
!  ZLEV(J,K) holds the interpolated vertical model grid coordinate for the
!  Kth observation using the surface pressure for the Jth ensemble member.
!  In addition, I'm not sure how compilers will treat the allocatable
!  members when the MODEL_OBS array elements are shuffled around as the
!  observations are processed.
!
    real,allocatable,save,private::zlev(:,:)  ! (nsol,nobs)
!
!  SURFPRES(J,K) is the interpolated surface pressure of the Jth ensemble
!  member corresponding to the Kth observation.
!
    real(D),allocatable,save,private::surfpres(:,:)  ! (nsol,nobs)
    private::raw_to_model_obs,quality_control,move_to_head,retained
!----------------------------------------------------------------------------
!  Private module data
!
!  TYPELIST contains all the observation types that we will assimilate
!  (default:  all observation types).  No leading or intervening blanks
!  are permitted in the string.
!
    character(NUM_OBS_TYPES),save,private::typelist='PTUV'
!
!  Units for temperature observations: C=Celsium, K=Kelvin, U=unknown
!  (determine from data).
!
    character,save,private::obs_temperature_units='U'
!
!  SURFHEIGHT holds the height of the model surface in meters.
!  For future use with real obs.
!
    real(D),allocatable,save,private::surfheight(:,:)
!
!  OBSLIST holds all the observations for the current latitude strip.
!
    type(model_obs),allocatable,save,private::obslist(:)
!
!  QCM_MAX gives the threshold for data acceptability; we retain only those
!  observations whose QCM < QCM_MAX.
!
    real(R),save,private::qcm_max=4.0
!
!  TOTAL_OBS_KEPT - the number of input observations whose type is one that
!  is to be assimilated.
!
    integer,save,private::total_obs_kept=0
!
!  HXBACK - set by H_OPERATOR.  The Kth column is H(X), where X is the
!  vector of ensemble background states closest to OBSLIST(K).
!  We retain extra precision here to try to offset cancelation that occurs
!  in the finite differences of the linear interpolation.
!
    real(D),allocatable,save,private::hxback(:,:)  ! (nsol,nobs)
!
!  BGSTAT - set by H_OPERATOR.  The Kth entry contains statistics like
!  the ensemble mean and standard error evaluated at the location of OBSLIST(K).
!
    type(ensemble_stat),allocatable,save,private::bgstat(:)  ! (nobs)
!
!  DHX - set by H_OPERATOR.  The Kth column is H(X(:)) - H(XBAR), where
!  X(:) is the ensemble solution interpolated to the Kth observation.
!  Extra precision is maintained here.
!
    real(D),allocatable,save,private::dhx(:,:)  ! (nsol,nobs)
!
!  Observation consistency checks
!
    integer,parameter,private::NO_CHECK=0,WHITAKER_CHECK=1,OBSERR_CHECK=2, &
       ABSDIFF_CHECK=3
    integer,save,private::consistency_checktype=NO_CHECK
    contains
!----------------------------------------------------------------------------
    subroutine copy_from_obs(loc,yobs,hbar,rdiag,nsol,nobs,hx,dh)
!  COPY_FROM_OBS - copy available observations from global grid to
!     a patch centered at grid location LOC = (ix,iy,iz).
!  NOBS  := on return, the number of available observations.
!  YOBS  := vector of observations
!  HBAR := vector of NOBS entries whose Kth entry is the Kth component of
!     H(xbar), where xbar is the ensemble mean over the patch.
!  RDIAG := vector whose Kth entry is R(K,K) (R is assumed diagonal).
!  NOBS := number of observations in the patch.
!  NSOL : number of ensemble solutions.
!  HX := If present, NSOL x NOBS matrix whose Kth column contains
!     H(x) for the Kth ensemble solution.  
!  DH := if present, (J,K)th element is the difference between H(x(j)) and
!     the mean of H(x) over the ensemble for the Kth observation.  DH is
!     the transpose of the matrix EB in Brian's writeup (and in Eq. A-22 of 
!     Szunyogh et al).
!  Note 1: The order in which variables are copied into YOBS must match
!  the ordering in copy_from_grid.
!  Note 2: We include the surface pressure in each vertical analysis.
!
    type(patch_index),intent(in)::loc
    real(D),allocatable,intent(out)::yobs(:),hbar(:),rdiag(:)
    real(D),allocatable,intent(out),optional::hx(:,:),dh(:,:)
    integer,intent(in)::nsol
    integer,intent(out)::nobs
!
!  Local variables
!
    integer::k,zlow,zhigh,xwin,ywin,idx,nobs2
    type(patch_index)::low,high
    type(obsinfo),allocatable::found(:),found2(:)
    logical::debug_dump
!
!  First look at the patch consisting of XWIN latitudes on either side of LOC.
!
    call patch_radius(loc,xwin,ywin)
    call zdiam(loc,zlow,zhigh)
    low%x=loc%x-xwin
    if(low%x.lt.0) low%x=low%x+LONF  ! wrap at 0 degrees
    low%y=max(loc%y-ywin,1)  ! no further south than the South Pole
    low%z=zlow
    high%x=loc%x+xwin
    if(high%x.gt.LONF) high%x=high%x-LONF
    high%y=min(loc%y+ywin,LATG)  ! no further north than the North Pole
    high%z=zhigh
    call search(low,high,found,nobs)
!
!  Include any surface pressure observations in the geographical region if
!  the patch doesn't otherwise include the surface.
!
    if(zlow.gt.1) then
       low%z=1  ! we assume surface pressure is at level 1 always
       high%z=1
       call search(low,high,found2,nobs2,(/PRES_TYPE/))  ! pressure obs only
       if(nobs2.gt.0) then  ! append to existing list
          call enlarge(found,nobs+nobs2)
          found(nobs+1:nobs+nobs2)=found2(1:nobs2)
          nobs=nobs+nobs2
       endif
    endif
!
!  Next check whether this patch extends over one of the poles.  If so, then
!  search the corresponding patch 180 degrees opposite.
!
    if(loc%y-ywin.lt.1.or.loc%y+ywin.gt.LATG) then  ! go across the pole
       low%x=low%x+LONF/2
       if(low%x.gt.LONF) low%x=low%x-LONF
       high%x=high%x+LONF/2
       if(high%x.gt.LONF) high%x=high%x-LONF
       low%z=zlow  ! in case we changed it looking for pressure obs
       high%z=zhigh
       if(loc%y-ywin.lt.1) then  ! narrower strip bordering appropriate pole
          low%y=1
          high%y=1-(loc%y-ywin)
       else
          high%y=LATG
          low%y=2*LATG-(loc%y+ywin)+1
       endif
       call search(low,high,found2,nobs2) 
       if(nobs2.gt.0) then  ! append to existing list
          call enlarge(found,nobs+nobs2)
          found(nobs+1:nobs+nobs2)=found2(1:nobs2)
          nobs=nobs+nobs2
       endif
!
!  Include any surface pressure obs on the opposite side of the pole.
!
       if(zlow.gt.1) then
          low%z=1
          high%z=1
          call search(low,high,found2,nobs2,(/PRES_TYPE/))
          if(nobs2.gt.0) then  ! append to existing list
             call enlarge(found,nobs+nobs2)
             found(nobs+1:nobs+nobs2)=found2(1:nobs2)
             nobs=nobs+nobs2
          endif
       endif
    endif
!
    if(nobs.eq.0) return
!
    allocate(yobs(nobs),hbar(nobs),rdiag(nobs))
    debug_dump=loc.eq.dump_loc
    if(debug_dump) open(DUMPUNIT,file='yobs.dump',status='replace')
    do k=1,nobs
       idx=found(k)%index
       yobs(k)=obslist(idx)%value
       rdiag(k)=obslist(idx)%sdev**2
       hbar(k)=bgstat(idx)%mean
       if(debug_dump) &
          write(DUMPUNIT,500) typelist(obslist(idx)%typ:obslist(idx)%typ), &
             obslist(idx)%x_grd,yobs(k),rdiag(k),found(k)%index
    enddo
500 format(a,1x,3f7.1,2es13.5,i6)
    if(debug_dump) close(DUMPUNIT)
    if(present(hx)) then
       allocate(hx(nsol,nobs))
       do k=1,nobs
          hx(:,k)=hxback(:,found(k)%index)
       enddo
    endif
    if(present(dh)) then
       allocate(dh(nsol,nobs))
       do k=1,nobs
          dh(:,k)=dhx(:,found(k)%index)
       enddo
    endif
    return
    end subroutine copy_from_obs
!----------------------------------------------------------------------------
    subroutine obs_init(ofile,expect,nobs,nsol)
!  OBS_INIT - read the observations, apply quality controls, and compute H.
!  OFILE  : name of observation file 
!  EXPECT : header to expect in ofile; 'none' if none.
!  NOBS  := number of observations that are retained and that pass 
!     applicable consistency checks.
!  NSOL  : number of ensemble solutions.
!  Side effects:
!  - reads and allocates OBSLIST and heaps it into a KD tree
!  - Precomputes H(x) for each observation Y by allocating and initializing
!    HXBACK and HXBAR.
!  Note: SET_OBSERVATION_TYPE must have been called first.
!
    character(*),intent(in)::ofile,expect
    integer,intent(out)::nobs
    integer,intent(in)::nsol
!
!  Local variables
!
    integer::j,info,kept,nobs_in
    type(obsinfo),allocatable::klist(:)
    type(raw_obs),allocatable::rawobs(:)
    type(model_obs),allocatable::obslist_temp(:)
    real::z(nsol)
    real(D)::pinterp(nsol)
    real,allocatable::zlev_temp(:,:)
    real(D),allocatable::surfpres_temp(:,:)
    logical::keep
!
!  NOTE:  The order in which the steps are performed is important, as the
!  various routines allocate and update the observation arrays and related
!  data.  Alter the order of the calls at your own risk!
!--------
!  Step 1.
!  Initialize OBSLIST with selected fields from the raw_obs input array
!  corresponding ONLY to those observations that we have decided to assimilate.
!
    call read_obs_slice(ofile,expect,nobs,rawobs)
    call raw_to_model_obs(obslist_temp,rawobs,nobs)
    deallocate(rawobs)
    nobs_in=nobs
!
!  Step 2.  Validate the observation values.  If any observation flunks this
!  validation test, then the program exits here.
!
    call validate(obslist_temp,nobs,info)
    if(info.ne.0) then
       write(IERR,100) j,info
       call fexit('obs_init: bad observation',1)
    endif
100 format('obs_init: bad observation :',i8,' info: ',i3)
!
!  Step 3.  Apply the conversion to energy units to match the units used in
!  the model grid.
!
    call obs_energy_units(obslist_temp,nobs)
!
!  Step 4A.  Obtain (x,y,z) model grid coords from the (lambda,phi,sigma)
!  entries and the background mean.
!
    where(obslist_temp(1:nobs)%kept) &
       obslist_temp(1:nobs)%kept=obslist_temp(1:nobs)%x_grd(3).ge.1.0.and. &
          obslist_temp(1:nobs)%x_grd(3).le.LEVS
    call move_to_head(obslist_temp,nobs,kept)
    nobs=kept
!
!  Step 4B.  Scale the vertical coordinate of each observation relative to
!  the surface pressure of each ensemble member.  KEEP is TRUE if the
!  vertical level associated with each ensemble member is between 1 and LEVS.
!  Note: the model surface pressure is in ln(dPa).
!
    allocate(zlev_temp(nsol,nobs),surfpres_temp(nsol,nobs))
    kept=0
    do j=1,nobs
       call sigma_fixup(obslist_temp(j),z,pinterp,nsol,keep)
       if(keep) then
          kept=kept+1
          obslist_temp(kept)=obslist_temp(j)
          zlev_temp(:,kept)=z
          surfpres_temp(:,kept)=pinterp
       endif
    enddo
    nobs=kept
    total_obs_kept=nobs  ! module variable
    write(LOGFILE,111) kept
111 format(i8,' observations retained with 1 <= sigma <= LEVS')
    call retained(obslist_temp,nobs)
!
!  Step 5.  Allocate the final observation list.
!  OBSLIST holds only those observations of the type that we will assimilate
!  and for which 1 <= sigma <= LEVS.  All other arrays related to grid model
!  interpolation are expected to be the same size as OBSLIST.
!  Life is MUCH easier if we allocate all the module arrays to be the
!  same size, so order the allocations to try to minimize holes in the heap.
!
    allocate(hxback(nsol,total_obs_kept),bgstat(total_obs_kept), &
       dhx(nsol,total_obs_kept))
    allocate(obslist(total_obs_kept))
    obslist=obslist_temp(1:total_obs_kept)
    deallocate(obslist_temp)
!
    allocate(zlev(nsol,total_obs_kept))
    zlev=zlev_temp(:,1:total_obs_kept)
    deallocate(zlev_temp)
!
    allocate(surfpres(nsol,total_obs_kept))
    surfpres=surfpres_temp(:,1:total_obs_kept)
    deallocate(surfpres_temp)
!
!  Step 5.
!  Precompute H(x) for x=background mean and x=individual ensemble members.
!  This procedure also fills in the sigma level for each observation
!  using the background mean.  H_OPERATOR also computes the interpolated
!  ensemble mean and standard error at each observation location and applies
!  observation consistency checks if requested.
!
    nobs=total_obs_kept
    call h_operator(nsol,nobs,kept)
    call retained(obslist,kept)
    nobs=kept
    total_obs_kept=nobs  ! module variable
!
!  Step 6.
!  Make the obslist into a KD tree using longitude and vertical coordinates.
!  Use the OBSINFO constructor to ensure that all components of KLIST are
!  initialized.
!
    write(LOGFILE,120) (100.0*total_obs_kept)/nobs_in
120 format(f8.2,' percent obs retention fraction')
    allocate(klist(total_obs_kept))
    do j=1,total_obs_kept
       klist(j)=obsinfo(obslist(j)%x_grd,j,obslist(j)%typ)
    enddo
    call build_kdtree(klist,total_obs_kept)
    return
    end subroutine obs_init
!----------------------------------------------------------------------------
    subroutine retained(obs,nobs)
!  RETAINED - report the number of observations by type.
!
    integer,intent(in)::nobs
    type(model_obs),intent(in)::obs(nobs)
!
    integer::j,obscount(NUM_OBS_TYPES)
!
    obscount=0
    do j=1,nobs
       obscount(obs(j)%typ)=obscount(obs(j)%typ)+1
    enddo
    write(LOGFILE,100) (TYPEOF(j),obscount(j),j=1,size(obscount))
100 format(9x,'retained observations by type'/9x,'type',3x,'count'/ &
       (9x,a1,3x,i8))
    return
    end subroutine retained
!----------------------------------------------------------------------------
    subroutine sigma_fixup(obs,z,pinterp,nsol,keep)
!  SIGMA_FIXUP - compute the appropriate vertical level of the observation 
!  relative to the interpolated surface pressure of each ensemble member.
!  For observations other than surface pressure:
!  sigma is determined based on the P_LEV value and the background mean
!  of the surface pressure (using linear interpolation from the model grid
!  to the observation location).  Z is determined by interpolating
!  the value of sigma so obtained with the fixed list of sigma levels
!  in resolution.f90.
!  NOTE: The observations and the background MUST be in consistent units
!  *before* this routine is called!
!  We discard observations whose Z value is outside the interval [1,LEVS].
!  OBS :=: the current observation (model_obs format).  The Z component
!     of the position is adjusted to 1 if the obs is surface pressure.
!  Z  := the interpolated observation height for each ensemble member
!  PINTERP := the interpolated surface pressure for each ensemble member
!  NSOL : number of ensemble solutions
!  KEEP := set to "true" if all the Z's are between 1 and LEVS.
!
    type(model_obs),intent(inout)::obs
    integer,intent(in)::nsol
    real,intent(out)::z(nsol)
    real(D),intent(out)::pinterp(nsol)
    logical,intent(out)::keep
!
!  Local variables
!
    real::x(3),dx(3),sigma(nsol),plev,ppa(nsol)
    real(D)::base(nsol),de(nsol,3),mu(0:2)
    type(patch_index)::loc,ix,iadj
    integer::j
!
    if(obs%p_lev.le.0.0) call fabort('sigma_fixup: obs p_lev <= 0')
!
!  Surface pressure observations don't have the Z coordinate defined.
!  We'll supply a value of 1.0 for our purposes.
!
    if(obs%typ.eq.PRES_TYPE) obs%x_grd(3)=1.0
!
!  Find the nearest integer-valued horizontal model grid indices that contain
!  the observation, accounting for the longitude wrap at LONF
!  (which is congruent to 0).
!
    x=obs%x_grd
    call gridpt_interp(x,ix,iadj,dx)
!
!  Interpolate the background surface pressure mean at X(1:2).
!  First the base point...
!
    loc=patch_index(ix%x,ix%y,1)
    call copy_grid_point(loc,base,PRES_TYPE,'background')
    if(ensemble_diffs_precomputed()) &
       call copy_grid_point(loc,mu(0:0),PRES_TYPE,'bmean')
!
!  ... then the nearest adjoining point to the north or south...
!  (making sure that the second index of DE matches DX!)
!
    loc=patch_index(ix%x,iadj%y,1)
    call copy_grid_point(loc,de(:,2),PRES_TYPE,'background')
    if(ensemble_diffs_precomputed()) &
       call copy_grid_point(loc,mu(2:2),PRES_TYPE,'bmean')
!
!  ... then the nearest adjoining point to the east or west.
!
    loc=patch_index(iadj%x,ix%y,1)
    call copy_grid_point(loc,de(:,1),PRES_TYPE,'background')
    if(ensemble_diffs_precomputed()) &
       call copy_grid_point(loc,mu(1:1),PRES_TYPE,'bmean')
!
!  Linearly interpolate between the adjoining points using the finite
!  difference approximation of the gradient vector.  Add back the
!  ensemble mean when the latter is precomputed, because we need the
!  interpolated surface pressure for each ensemble member.
!
    pinterp=base+dx(1)*(de(:,1)-base)+dx(2)*(de(:,2)-base)
    if(ensemble_diffs_precomputed()) &  ! interpolate respective means
       pinterp=pinterp+(mu(0)+dx(1)*(mu(1)-mu(0))+dx(2)*(mu(2)-mu(0)))
!
!  The p_lev pressure level measurement must be scaled in the same way as
!  the background surface pressure.  Here we assume energy units.
!
    if(obs%typ.eq.PRES_TYPE) then  ! model Z value is 1
       z=1.0
       keep=.true.
    else
       ppa=real(10.0*exp(energy_to_pres(pinterp)),kind(ppa))  ! to hPa
       sigma=real(obs%p_lev/ppa,kind(sigma))
       do j=1,nsol
          z(j)=index_interp(sigma(j),SIGMA_LEV,LEVS)
       enddo
       keep=all(z.ge.1.0.and.z.le.LEVS)
    endif
    return
    end subroutine sigma_fixup
!----------------------------------------------------------------------------
    subroutine h_operator(nsol,nobs,kept)
!  H_OPERATOR - implement the H operator that maps the model space to
!  the observation space.  In this version, we call it once at startup
!  to interpolate to each of the available observations.
!  For now, we do a linear interpolation of each observation to the
!  model grid.  (In the future, other interpolations might be done depending
!  on the variable type.)  In addition, we apply any requested data
!  consistency checks.
!  Arguments:
!  NSOL : number of ensemble solutions
!  NOBS : the number of entries in OBSLIST.
!  KEPT := the number of observations that are retained (e.g., which pass
!    any applicable internal consistency checks)
!
!  Side effects:
!  HXBACK := Kth col holds interpolated model grid values at the obs location
!    (i.e., H(x)) for each of the NSOL ensemble solutions corresponding to 
!    the observation type (e.g., the NSOL interpolated values of the
!    temperature if the Kth obs is a temperature observation).
!  DHX := Kth column is H(x(k))-H(xbar) at Kth obs location
!  BGSTAT := Kth entry is interpolated ensemble mean and std error at 
!    the location of the Kth obs.
!
!  Note:  Any scaling to the observations and model variables 
!  MUST have already been performed!
!
    integer,intent(in)::nsol,nobs
    integer,intent(out)::kept
!
!  Local variables
!
    integer::k
!
    if(nobs.eq.0) return
    do k=1,nobs
       call linear_interp(k,nsol,hxback(:,k),'background',bgstat(k),dhx(:,k))
    enddo
    if(consistency_checktype.eq.NO_CHECK) then
       kept=nobs
    else
       call quality_control(obslist,bgstat,nobs,kept)
    endif
    return
    end subroutine h_operator
!----------------------------------------------------------------------------
    subroutine linear_interp(idx,n,hx,gridname,ensstat,dh)
!  LINEAR_INTERP - interpolate the model grid to the model_obs OBS
!  by linear interpolation from the nearest model grid points.
!  IDX  : index of the observation (MODEL_OBS format)
!  N  : number of ensemble solutions, either ENSEMBLE_SIZE() or 1.
!     In the latter case, we assume interpolation using the ensemble mean.
!     Otherwise, iwe work with all the individual ensemble solutions.
!  HX := interpolated model grid values at the observation location (i.e., H(x))
!    for each of the N ensemble solutions corresponding to the observation
!    type (e.g., the N interpolated values of the temperature if OBS is
!    a temperature observation).
!  GRIDNAME : which model grid: 'background', 'analysis', or 'nature'.
!  ENSSTAT := Holds ensemble mean and standard error at the interpolated 
!    point corresponding to the observation (with the same scaling as applied
!    to the rest of the observations).
!  DH := if present, Kth entry is the the difference between H(x(k)) and
!    the mean of the H's, where x(k) is the kth ensemble solution.
!
    integer,intent(in)::idx,n
    real(D),intent(out)::hx(n)
    character(*),intent(in)::gridname
    type(ensemble_stat),intent(out),optional::ensstat
    real(D),intent(out),optional::dh(n)
!
!  Local variables
!
    type(model_obs)::obs
    real(D)::base(n),dg(n,3),mu(n,0:3),dmu(n),dz(n)
    type(patch_index)::ix,iadj,loc
    real(kind(obs%x_grd))::x(3),dx(3)
    integer::j,iz(n),iadz(n)
    character(5)::themean
!
!  Ensure that this is a variable type that we know how to handle.
!  For now this is 1, 2, 3, or 4 - see grid_types for labeling.
!
    if(idx.lt.1.or.idx.gt.size(obslist)) &  ! was total_obs_kept
       call fabort('linear_interp: obs index out of range')
    if(n.ne.1.and.n.ne.ensemble_size()) &
       call fabort('linear_interp: bad solution count')
    obs=obslist(idx)
    if(obs%typ.lt.1.or.obs%typ.gt.NAVARS) &
       call fabort('linear_interp: unknown obs type')
!
!  Special case:  Surface pressure and background grid.
!  We've already done the required interpolations when we
!  computed the observation height in SIGMA_FIXUP.
!
    if(obs%typ.eq.PRES_TYPE.and.gridname(1:1).eq.'b') then  ! background
       if(n.eq.1) then  ! ensemble mean
          hx(1)=sum(surfpres(:,idx))/ensemble_size()
       else
          hx=surfpres(:,idx)  ! previously interpolated background
       endif
       call ensstat_and_dh
       return
    endif
!
!  Otherwise we have a measurement that we need to interpolate according
!  to the vertical height of the observation (except for the surface pressure).
!  The height in turn depends on the surface pressure in each ensemble member
!  (except for surface pressure, for which z=1 identically).
!
    x(1:2)=obs%x_grd(1:2)  ! longitude and latitude only
    do j=1,n
       x(3)=merge(1.0,zlev(j,idx),obs%typ.eq.PRES_TYPE)
       call gridpt_interp(x,ix,iadj,dx)  ! find adjacent model grid points
       iz(j)=ix%z    ! save the varying vertical coordinate
       iadz(j)=iadj%z  ! and that of the opposite grid point
       dz(j)=abs(dx(3))
    enddo
    loc=ix  ! for the longitude and latitude indices
!
!  Interpolate the ensemble solutions.  First the base point...
!
    themean=gridname(1:1)//'mean'
    call copy_grid_z(loc,iz,base,obs%typ,gridname)
    if(n.gt.1.and.ensemble_diffs_precomputed(gridname)) &
       call copy_grid_z(loc,iz,mu(:,0),obs%typ,themean)
!
!  ...then the point to the north or south of the base point...
!  (to make this part of the code more cache-friendly - but ensure
!  that the subscripting of DG and MU is compatible with DX!)
!
    loc=patch_index(ix%x,iadj%y,ix%z)
    call copy_grid_z(loc,iz,dg(:,2),obs%typ,gridname)
    if(n.gt.1.and.ensemble_diffs_precomputed(gridname)) &
       call copy_grid_z(loc,iz,mu(:,2),obs%typ,themean)
    dg(:,2)=dg(:,2)-base
!
!  ...then the point above or below as appropriate...
!  (to make this part of the code more cache-friendly)
!
    if(obs%typ.ne.PRES_TYPE) then
       loc=ix
       call copy_grid_z(loc,iadz,dg(:,3),obs%typ,gridname)
       if(n.gt.1.and.ensemble_diffs_precomputed(gridname)) &
          call copy_grid_z(loc,iadz,mu(:,3),obs%typ,themean)
       dg(:,3)=dg(:,3)-base
    endif
!
!  ...and finally the point to the east or west.
!
    loc=patch_index(iadj%x,ix%y,ix%z)
    call copy_grid_z(loc,iz,dg(:,1),obs%typ,gridname)
    if(n.gt.1.and.ensemble_diffs_precomputed(gridname)) &
       call copy_grid_z(loc,iz,mu(:,1),obs%typ,themean)
    dg(:,1)=dg(:,1)-base
!
!  Linear interpolation.  The longitude and latitude
!  components are constant for all ensemble members; the last term
!  completes the interpolation in the vertical direction for non-pressure obs.
!
    hx=base+dg(:,1)*dx(1)+dg(:,2)*dx(2)
    if(obs%typ.ne.PRES_TYPE) hx=hx+dg(:,3)*dz
!
!  If N=1 then we're already interpolating the background or analysis
!  mean, presumably.  Otherwise add back the interpolated
!  mean when the ensemble differences are precomputed.
!  The update is mu(0) + (mu(1)-mu(0))*dx(1) + (mu(2)-mu(0))*dx(2) +
!  (mu(3)-mu(0))*dz; the latter term includes all individual
!  vertical displacements.
!
    if(n.gt.1.and.ensemble_diffs_precomputed(gridname)) then
       dmu=1.0-dx(1)-dx(2)-dz
       hx=hx+mu(:,0)*dmu+mu(:,1)*dx(1)+mu(:,2)*dx(2)
       if(obs%typ.ne.PRES_TYPE) hx=hx+mu(:,3)*dz
    endif
    call ensstat_and_dh
    return
!
!  If requested, compute statistics about the observation and the model
!  grid at the corresponding location.
!
    contains
       subroutine ensstat_and_dh
       real(D)::emean,ediff(n)
!
       if(present(ensstat).or.present(dh)) then
          emean=sum(hx)/ensemble_size()
          ediff=hx-emean
          !!obslist(idx)%bgmean=emean  !debug
          !!obslist(idx)%bgsdev=sqrt(dot_product(ediff,ediff)/n)  !debug
       endif
       if(present(dh)) dh=ediff
       if(present(ensstat)) then
          ensstat%mean=emean
          ensstat%sdev=sqrt(dot_product(ediff,ediff)/n)
       endif
       return
       end subroutine ensstat_and_dh
    end subroutine linear_interp
!----------------------------------------------------------------------------
    subroutine apply_consistency_check(checktype)
!  APPLY_CONSISTENCY_CHECK - argument determines which type of observation
!  consistency check to apply to the data.
!
    character(*),intent(in)::checktype
!
    select case(checktype)
    case('Whitaker','whitaker')
       consistency_checktype=WHITAKER_CHECK
    case('observation','obs','Observation','Obs')
       consistency_checktype=OBSERR_CHECK
    case('absdiff')
       consistency_checktype=ABSDIFF_CHECK
    case('none')
       consistency_checktype=NO_CHECK
    case default
       call fabort('apply_consistency_check: unknown check type ' &
          //trim(checktype))
    endselect
    write(LOGFILE,100) trim(checktype)
100 format('observation consistency check: ',a)
    return
    end subroutine apply_consistency_check
!----------------------------------------------------------------------------
    subroutine set_observation_types(list,temp_units)
!  SET_OBSERVATION_TYPES - determine which observation types will be
!  assimilated.
!  The current procedure uses single-letter variable names; later we'll have
!  to devise another convention.  For now, they are:
!  'P' or 'p' for surface pressure;
!  'T' or 't' for virtual temperature;
!  'U' or 'u' for u-wind components;
!  'V' or 'v' for v-wind components.
!  TEMP_UNITS is either 'C' or 'K' for observations in Celsius or Kelvin,
!     respectively.
!
    character(*),intent(in),optional::list
    character(*),intent(in),optional::temp_units
!
    integer::j,n
!
    if(present(list)) then
       n=0
       typelist=' '  ! clear all previous entries
       do j=1,min(len(typelist),len(list))
          select case(list(j:j))
          case('P','p')
             n=n+1
             typelist(n:n)='P'
          case('T','t')
             n=n+1
             typelist(n:n)='T'
          case('U','u')
             n=n+1
             typelist(n:n)='U'
          case('V','v')
             n=n+1
             typelist(n:n)='V'
          endselect
       enddo
    endif
    if(present(temp_units)) then
       if(index('CKU',temp_units(1:1)).gt.0) then
          obs_temperature_units=temp_units(1:1)
       else
          call fexit('set_observation_types: unknown temp_units value: '// &
             temp_units,2)
       endif
    endif
    return
    end subroutine set_observation_types
!----------------------------------------------------------------------------
    pure subroutine get_observation_types(list,num,temp_units)
!  GET_OBSERVATION_TYPES - obtain the number and type of the observations
!  that are being assimilated.
!  NUM  - if present, the number of observation types that are assimilated.
!  LIST - if present, a character string of length NUM_OBS_TYPES whose first
!  NUM entries on return contain a list of the observation types being
!  assimilated.
!  TEMP_UNITS - if present, either set to 'C' or 'K' for Celsius or Kelvin,
!     respectively,
!
    integer,intent(out),optional::num
    character(*),intent(out),optional::list,temp_units
!
    if(present(num)) num=len_trim(typelist)
    if(present(list)) list=typelist
    if(present(temp_units)) temp_units=obs_temperature_units
    return
    end subroutine get_observation_types
!----------------------------------------------------------------------------
    subroutine set_qcm_threshold(q)
!  SET_QCM_THRESHOLD - set the QCM_MAX threshold past which we omit data.
!
    real,intent(in)::q
!
    if(q.le.1.0) call fexit('set_qcm_threshold: bad value',2)
    qcm_max=q
    return
    end subroutine set_qcm_threshold
!----------------------------------------------------------------------------
    subroutine raw_to_model_obs(mdlobs,raw,nobs)
!  RAW_TO_MODEL_OBS - For each of the observation types that we will assimilate,
!  copy a subset of the raw_obs components to the model_obs.
!  Filter the data according to quality control mark and omit those observations
!  whose QCM >= QCM_MAX.
!  Note 1: NO SCALING is applied here; VALIDATE, below, should be
!  called to check for sensible observation values and to apply any scaling
!  consistent with the model grid.
!  Note 2: The sigma level of the observation is obtained from the P_LEV
!  and PRSS members of RAW.  We discard observations whose model Z valud is
!  less than 1 or greater than LEVS.
!  NOBS :=: On entry, the number of observations in RAW.  On return,
!     the number of observations of the types that we will assimilate,
!     which is also the allocated size of MDLOBS.
!
    integer,intent(inout)::nobs
    type(model_obs),allocatable,intent(out)::mdlobs(:)
    type(raw_obs),intent(in)::raw(nobs)
!
    integer::j,n,obscount(NUM_OBS_TYPES)
    logical::kept(nobs)
    character(4)::field
    character(NUM_OBS_TYPES)::list
    real::sigma
!
!  Count and flag the observations that are of a type that we will
!  assimilate.  Also retain only those observations whose quality control
!  marks are sufficiently low.
!
    do j=1,nobs
       field=adjustl(raw(j)%field)
       kept(j)=index(typelist,field(1:1)).ne.0.and.raw(j)%qcm.lt.qcm_max
    enddo
    n=count(kept)
    allocate(mdlobs(n))
!
!  This version handles synthetic data only.  Eventually we will alter
!  the handling of the PRSS field.  For the moment, we obtain the vertical
!  level of the observation from the SIGMA field of the raw observation.
!
    n=0
    obscount=0
    do j=1,nobs
       if(.not.kept(j)) cycle
       n=n+1
       mdlobs(n)%x_grd(1)=xgrid(raw(j)%lambda)
       mdlobs(n)%x_grd(2)=ygrid(raw(j)%phi)
       mdlobs(n)%x_grd(3)=raw(j)%sigma
       mdlobs(n)%lambda=raw(j)%lambda
       mdlobs(n)%phi=raw(j)%phi
       mdlobs(n)%time_diff=raw(j)%time_diff
       mdlobs(n)%p_lev=raw(j)%p_lev
       mdlobs(n)%value=raw(j)%value
       mdlobs(n)%sdev=raw(j)%sdev
       mdlobs(n)%typecode=raw(j)%typecode  ! code number for observation type
       mdlobs(n)%kept=.true.  ! unless it fails a later consistency check
       field=adjustl(raw(j)%field)
       select case(field)  ! which must appear in TYPELIST
       case('P')
          mdlobs(n)%typ=PRES_TYPE
       case('T')
          mdlobs(n)%typ=TEMP_TYPE
       case('U')
          mdlobs(n)%typ=UWIND_TYPE
       case('V')
          mdlobs(n)%typ=VWIND_TYPE
       endselect
       obscount(mdlobs(n)%typ)=obscount(mdlobs(n)%typ)+1
    enddo
    call get_observation_types(list)
    write(LOGFILE,100) list
100 format('raw_to_model_obs: Observation types assimilated: ',a)
    write(LOGFILE,101) qcm_max
101 format('observation quality control threshold:',f5.1)
    write(LOGFILE,102) nobs,'observations in data file'
102 format(i8,1x,a)
    write(LOGFILE,102) n,'obs retained based on type and quality control'
    nobs=n
    call retained(mdlobs,nobs)
    return
    end subroutine raw_to_model_obs
!----------------------------------------------------------------------------
    subroutine gridpt_interp(x,ix,iadj,dx)
!  GRIDPT_INTERP - given an observation at location X(1:3), which is a
!  floating-point number relative to the model grid units, find the nearest
!  point IX on the model grid (i.e., integer coordinates) and the next
!  closest adjoining model grid point.
!  X  : the (x,y,z) coordinates of the point (as a 3-vector)
!  IX  := PATCH_INDEX (integer) coordinates of the nearest corner point
!  IADJ := PATCH_INDEX (integer) coordinates of the opposite corner to IX.
!  DX  := a 3-vector giving the |X-IX| (pointwise *distances*)
!  Note 1:  We cannot yet handle the case where X(3) < 1; if we get a vertical
!  coordinate that lies "below" the first model level, we abort.
!  Note 2:  If an observation is above the top level of the model atmosphere,
!  then treat it as though it were exactly at the top of the atmosphere.
!
    real,intent(in)::x(3)
    type(patch_index),intent(out)::ix,iadj
    real,intent(out)::dx(3)
!
!  IX is the nearest integer grid point.
!
    ix=patch_index(nint(x(1)),nint(x(2)),nint(x(3)))
    dx=(/x(1)-ix%x,x(2)-ix%y,x(3)-ix%z/)
    call xy_adjust(ix,0,0,iadj)  ! handle longitude wrap if necessary
    ix=iadj
!
!  Project obs above top of model atmosphere onto the top.
!
    if(x(3).gt.LEVS) then 
       ix%z=LEVS
       dx(3)=0.0
    endif
!
!  Find the opposite corner of the smallest enclosing box, handling regions
!  over the poles and straddling the 0-degree longitude line as needed.
!
    call xy_adjust(ix,signum(dx(1)),signum(dx(2)),iadj)
    iadj%z=min(max(ix%z+signum(dx(3)),1),LEVS)
    dx=abs(dx) ! convert to distances
    return
    end subroutine gridpt_interp
!----------------------------------------------------------------------------
    subroutine quality_control(obslist,es,nobs,kept)
!  QUALITY_CONTROL - Retain only those observations for which
!  the interpolated ensemble mean is sufficiently close to the 
!  observation value.
!
!  The Whitaker criterion requires that the obs value and the background
!  mean lie within 5 ensemble standard deviations of each other when the
!  background ensemble grid is interpolated to the observation location.
!  The "obserr" criterion is similar but the obs value and mean must lie
!  with 6 times the maximum of either the obs std error or the ensemble std err.
!  The "absdiff" criterion is like the obserr criterion except that
!  the threshold is further expanded to be the larger of the 0.1 * observation
!  absolute value or the obserr threshold.
!  NOBS  :  the number of available observations.
!  OBSLIST :=: On entry, the list of observations.  On return, the KEPT field
!    is updated to reflect whether the Kth observation should be retained.
!    On return the first KEPT entries are the retained observations and
!    the remainder are not used in the data assimilation.
!  ES  :=:  On entry, Kth entry is to the interpolated ensemble mean
!     and standard error corresponding to the Kth input observation.
!     On return, entries are shuffled to match OBSLIST.  
!  KEPT := scalar giving the total number of retained observations.
!  NOTE: Take care that both the ensemble and the observation variables
!  are scaled consistently!
!
!  Side effects:
!    Shuffles HXBACK, ES, DHX, ZLEV, and SURFPRES so that their Kth
!    entries correspond to the Kth entry of OBSLIST.
!
    integer,intent(in)::nobs
    integer,intent(out)::kept
    type(model_obs),intent(inout)::obslist(nobs)
    type(ensemble_stat),intent(inout)::es(nobs)
!
!  Local variables
!
    real(OBSKIND)::threshold
    type(inflation_info)::infl
    type(patch_index)::loc
    integer::j,ix(3),vars(NAVARS),index(nobs)
!
    if(consistency_checktype.eq.NO_CHECK) return
!
    do j=1,nobs
       threshold=6.0*es(j)%sdev
       ix=nint(obslist(j)%x_grd)
       loc=patch_index(ix(1),ix(2),ix(3))
       call inflation(loc,infl,vars)
       if(any(vars.eq.obslist(j)%typ)) threshold=threshold*(1.0+infl%covar)
       if(consistency_checktype.ne.WHITAKER_CHECK) then
          threshold=max(threshold,real(6.0*obslist(j)%sdev,kind(threshold)))
       else if(consistency_checktype.eq.ABSDIFF_CHECK) then
          threshold=max(threshold, &
             real(abs(obslist(j)%value)/10.0,kind(threshold)))
       endif
       if(abs(es(j)%mean-obslist(j)%value).gt.threshold) obslist(j)%kept=.false.
    enddo
!
!  Move retained observations to the head of the list and
!  shuffle related arrays so that all related entries correspond.
!
    do j=1,nobs
       index(j)=j
    enddo
    call move_to_head(obslist,nobs,kept,index)
    if(kept.lt.nobs) then
       hxback=hxback(:,index)
       es=es(index)
       dhx=dhx(:,index)
       zlev(:,1:nobs)=zlev(:,index)
       surfpres=surfpres(:,index)
    endif
    write(LOGFILE,100) kept
100 format(i8,' observations pass quality control checks')
    return
    end subroutine quality_control
!----------------------------------------------------------------------------
    subroutine validate(obs,nobs,info)
!  VALIDATE - make a sequence of heuristic checks on the observation
!  data and compute (x,y) grid coordinates from the supplied (lambda,phi)
!  coordinates.  Set INFO to 0 if all checks are passed and to a nonzero value
!  otherwise.   
!  NOTE 1: this routine intuits the temperature units (either C or K) and
!  converts them as required to K.  Otherwise, it does not perform any other
!  scalings.
!  NOTE 2: Some of the "real fake" data sets have sigma=9999 for surface
!  pressure observations.  Replace these with sigma=1.
!
    integer,intent(in)::nobs
    type(model_obs),intent(inout)::obs(nobs)
    integer,intent(out)::info
!
    integer::j
!
    info=0
!
    do j=1,nobs
       select case(obs(j)%typ)
       case(PRES_TYPE)  ! can be either hPa or ln(dPa)
          if(obs(j)%value.le.0.0) then
             call dump_obs(obs(j),'nonpositive pressure',real(info))
             exit
          endif
          obs(j)%x_grd(3)=1.0
       case(TEMP_TYPE)  ! internal units are Kelvin; convert C to K as needed
          if(obs_temperature_units.eq.'U') then  ! intuit units from first obs
             if(obs(j)%value.le.0.0) then
                obs_temperature_units='C'
             else if(obs(j)%value.ge.150.0) then
                obs_temperature_units='K'
             else
                call fexit('validate: can''t determine temperature units',2)
             endif
          endif
          if(obs_temperature_units.eq.'C') then
             obs(j)%value=obs(j)%value+273.15
          else if(obs(j)%value.le.0.0) then
             info=12
             call dump_obs(obs(j),'nonpositive temperature',real(info))
             exit
          endif
       case(UWIND_TYPE:VWIND_TYPE)  ! in m/s
          if(abs(obs(j)%value).ge.500.0) then
             info=13
             call dump_obs(obs(j),'huge wind speed',real(info))
             exit
          endif
       case default  ! unknown type
          info=5
          call dump_obs(obs(j),'unknown obs type',real(info))
          exit
       endselect
!
       if(obs(j)%sdev.le.0.0) then
          info=4
          call dump_obs(obs(j),'nonpositive standard error',real(info))
          exit
       endif
    enddo
    return
    end subroutine validate
!----------------------------------------------------------------------------
    subroutine obs_energy_units(obs,nobs)
!  OBS_ENERGY_UNITS - rescale the observation array to energy units.
!  The components of OBS are assumed to have been validated.
!  NOTE:  The energy unit conversion is fixed by the sigma level of
!  the observation, which is independent of the vertical positioning
!  of the observation based on the individual ensemble members.
!
    integer,intent(in)::nobs
    type(model_obs),intent(inout)::obs(nobs)
!
    integer::j
!
!  Ensure that the spectral coefficients are initialized.
!
    call spectral_code_init
!
    do j=1,nobs
    !
    !  Pressure values < 10 are assumed to be in ln(dPa); if sdev <= 0.001
    !  then substitute a synthetic value (4.0e-04, to match gridded obs)
    !  The pressure level (P_LEV) field is always in hPa.
    !
       if(obs(j)%typ.eq.PRES_TYPE) then
          if(obs(j)%value.le.10.0) then  ! assume ln(dPa) already
             obs(j)%value=pres_to_energy(obs(j)%value)
             if(obs(j)%sdev.le.1.0e-03) then  ! substitute artificial value
                obs(j)%sdev=pres_to_energy(1.0)*4.0e-04  ! wq*val
             else
                obs(j)%sdev=pres_to_energy(1.0)*obs(j)%sdev  ! wq*val
             endif
          else    ! hPa -> ln(dPa) * energy scaling
             obs(j)%sdev=pres_to_energy(obs(j)%sdev/obs(j)%value)
                         ! log((v+s)/10) ~ log(v/10)+s/v 
             obs(j)%value=pres_to_energy(log(obs(j)%value/10.0))
          endif
          if(obs(j)%value.le.0.0) call fabort('obs_energy_units: P < 0')
       else
          if(obs(j)%typ.eq.TEMP_TYPE) then
             obs(j)%value=temp_to_energy(obs(j)%value,obs(j)%x_grd(3))
             obs(j)%sdev=temp_to_energy(obs(j)%sdev,obs(j)%x_grd(3))
             if(obs(j)%value.le.0.0) call fabort('obs_energy_units: T < 0')
          else
             obs(j)%value=wind_to_energy(obs(j)%value,obs(j)%x_grd(3))
             obs(j)%sdev=wind_to_energy(obs(j)%sdev,obs(j)%x_grd(3))
          endif
       endif
       if(obs(j)%sdev.le.0.0) call fabort('obs_energy_units: negative sdev')
    enddo
    return
    end subroutine obs_energy_units
!----------------------------------------------------------------------------
    subroutine move_to_head(list,nobs,kept,index)
!  MOVE_TO_HEAD - Move retained observations to the head of each list.
!  LIST :=: array of NOBS observations, shuffled as needed so that the
!     entries 1:KEPT contain retained observations and KEPT+1:NOBS holds
!     the rest.
!  INDEX :=: if present, an index array.  The elements of this array are
!     moved so that if INDEX(K) corresponds to LIST(K) on entry, then
!     the same is true on return.
!
    integer,intent(in)::nobs
    integer,intent(out)::kept
    type(model_obs),intent(inout)::list(nobs)
    integer,intent(inout),optional::index(nobs)
!
    type(model_obs),allocatable::otemp(:)  ! avoid stack overflow
    integer,allocatable::newindex(:)  ! avoid stack overflow
    integer::j,k,d
!
!  Copy the original list to a new one: retained observations at the
!  beginning, deleted ones at the end.
!
    allocate(otemp(nobs))
    k=0
    d=nobs+1
    if(present(index)) then
       allocate(newindex(nobs))
       do j=1,nobs
          if(list(j)%kept) then
             k=k+1
             otemp(k)=list(j)
             newindex(k)=index(j)
          else
             d=d-1
             otemp(d)=list(j)
             newindex(d)=index(j)
          endif
       enddo
       index=newindex
    else
       do j=1,nobs
          if(list(j)%kept) then
             k=k+1
             otemp(k)=list(j)
          else
             d=d-1
             otemp(d)=list(j)
          endif
       enddo
    endif
    kept=k
    list=otemp
    return
    end subroutine move_to_head
!----------------------------------------------------------------------------
    subroutine dump_obs(obs,msg,optval,raw)
    use io
    type(model_obs),intent(in)::obs  ! unscaled entries
    character(*),intent(in)::msg
    real,intent(in),optional::optval
    type(raw_obs),intent(in),optional::raw
!
    write(IERR,900) msg
900 format('------------'/'model_obs dump: ',a)
    write(IERR,901) 'location (x,y,z): ',obs%x_grd
901 format(a,3f11.3)
    if(present(optval)) write(IERR,901) 'optval: ',optval
    write(IERR,902) 'value, sdev: ',obs%value,obs%sdev
902 format(a,2es13.4)
    write(IERR,903) 'p_lev, type: ',obs%p_lev,obs%typ
903 format(a,e12.3,i4)
!
    if(.not.present(raw)) return
    write(IERR,910) raw%typecode,raw%lev
910 format('raw_obs dump: typecode: ',i0,' lev: ',i0)
    write(IERR,911) raw%x_grd,raw%y_grd,raw%lambda,raw%phi
911 format('x_grd, y_grd: ',2f11.4,' lambda, phi: ',2f11.4)
    write(IERR,912) raw%sigma,raw%p_lev
912 format('sigma, p_lev: ',f11.4,es12.3)
    write(IERR,913) raw%value,raw%sdev,raw%qcm
913 format('value, sdev, qcm: ',3es13.4)
    return
    end subroutine dump_obs
!----------------------------------------------------------------------------
    pure subroutine kth_obs(k,obs,ierr,bg,zmean)
!  KTH_OBS - return the Kth retained model_obs.  Convenience function.
!  K : the observation to get.
!  OBS := the Kth observation, if available.
!  IERR := set to 0 unless K > NOBS or K < 1.
!  BG := if present, and if K is in range, then the background ensemble
!     statistics computed in LINEAR_INTERP corresponding to the Kth observation.
!  ZMEAN := if present, and K is in range, then the average value of the
!    model Z coordinate, based on the surface pressure adjustment from the
!    background.
!
    integer,intent(in)::k
    type(model_obs),intent(out)::obs
    integer,intent(out)::ierr
    type(ensemble_stat),intent(out),optional::bg
    real,intent(out),optional::zmean
!
    if(k.ge.1.and.k.le.size(obslist)) then
       obs=obslist(k)
       if(present(bg)) bg=bgstat(k)
       if(present(zmean)) zmean=sum(zlev(:,k))/ensemble_size()
       ierr=0
    else
       ierr=-1
    endif
    return
    end subroutine kth_obs    
!----------------------------------------------------------------------------
    pure subroutine obs_count(kept,total)
!  OBS_COUNT - return the total number of observations read in and the
!  total number of obs that are retained.
!
    integer,intent(out),optional::kept,total
!
    if(present(kept)) kept=total_obs_kept
    if(present(total)) total=size(obslist)
    return
    end subroutine obs_count
!----------------------------------------------------------------------------
    subroutine read_obs_slice(ofile,expect,nobs,obsarray)
!  READ_OBS_SLICE - read a raw_obs file corresponding to a slice
!  of the background grid.
!  The latitudinal strip to which the observations correspond must
!  exactly match SOUTH and NORTH.
!  OFILE   : name of observations file
!  EXPECT  : 4-character header to expect ('none' => no header)
!  NOBS  := number of observations read and stored into OBSLIST after
!    passing simple validity tests.
!  Side effects:
!  - opens and closes INNOSWAP.
!  - allocates and fills OBSLIST
!  - heaps the observations into a KD tree for fast lookup
!
    character(*),intent(in)::ofile,expect
    integer,intent(out)::nobs
    type(raw_obs),allocatable,intent(out)::obsarray(:)
!
!  Local variables
!
    integer::istart,istop,iseq,info,south,north,nobsin
    character(4)::ftype
    character(80)::msg
    type(raw_obs)::obs
!
    open(unit=INNOSWAP,form='unformatted',status='old',file=ofile,iostat=info)
    if(info.ne.0) call fexit(trim(ofile)//': cannot open',2)
    if(verbosity().gt.VERB_TERSE) call finfo('reading '//trim(ofile)//'...')
!
!  Check headers for consistency with the expected latitude strip.
!  Endian issues are generally evident at this point.
!
    if(expect.ne.'none') then
       read(INNOSWAP) ftype,iseq,istart,istop,nobsin
       if(ftype.ne.'obs') then
          call finfo(trim(ofile)//': wrong header: '//ftype)
          call fexit(trim(ofile)//': not a split obs file',2)
       endif
       if(nobsin.lt.0) call fexit(trim(ofile)//': negative number of obs',2)
       call model_grid_boundary(south,north)
       if(istart.gt.south.or.istop.lt.north) then
          write(IERR,100) trim(ofile),istart,istop,south,north
100       format(a,' obs file has istart=',i0,' istop=',i0, &
             ' but program requires south=',i0,' and north=',i0)
          call fexit(trim(ofile)//': latitude limits disagree with grid',2)
       endif
    else
       read(INNOSWAP) iseq
    endif
    if(iseq.ne.seqid()) then
       write(IERR,110) trim(ofile),iseq,seqid()
110    format(a,': sequence id is: ',i0,' expected: ',i0)
       call fexit(trim(ofile)//': unexpected sequence id in input',2)
    endif
!
!  Allocate, read, and verify the observation grid.
!
    allocate(obsarray(nobsin))
    nobs=0
    do
       read(INNOSWAP,end=10) obs
       nobs=nobs+1
       if(nobs.gt.size(obsarray)) call grow(obsarray)
       obsarray(nobs)=obs
    enddo
 10 continue
    if(nobs.ne.nobsin) then
       write(msg,115) nobsin,nobs
       call fabort(msg)
    endif
115 format('read_obs_slice: expected ',i0,' obs; got ',i0)
    close(INNOSWAP)
    if(verbosity().gt.VERB_TERSE) then
       write(msg,120) 'number of observations: ',nobs
       call finfo(trim(msg))
       call finfo('done reading observations')
    endif
120 format(a,i8)
    return
    contains
! -------------------
       subroutine grow(f)
   !  GROW - enlarge F by 50% or 1 element, whichever is greater.
   !  Preserve the initial elements in F.
   !
       type(raw_obs),allocatable,intent(inout)::f(:)
   !
       type(raw_obs),allocatable::temp(:)
       integer::sz
   !
       sz=size(f)
       allocate(temp(sz))
       temp=f
       deallocate(f)
       allocate(f(max(3*sz/2,sz+1)))
       f(1:sz)=temp
       deallocate(temp)
       return
       end subroutine grow
! -------------------
    end subroutine read_obs_slice
!----------------------------------------------------------------------------
    subroutine write_obs_slice(ofile,tag,south,north,obs,nobs)
! WRITE_OBS_SLICE - write a section of the global observation grid
!  to the indicated output file in latitude-major format.
!  Halt the program if an error occurs.
!  OFILE   : name of output observation file
!  TAG  : 4-character header tag to use.  If set to 'none', we
!            write no tag other than the sequence id.
!  SOUTH  :  Southernmost latitude index of output observations.
!  NORTH  :  Northernmost latitude index of output observations.
!  OBS  :  Array containing all global observations.
!  NOBS  : Number of observations in OBS (not necessarily the number
!      written out)
!  Side effects:
!  - Opens and closes OUTNOSWAP.
!  - Halts the code if OFILE cannot be opened for writing
!
    character(*),intent(in)::ofile
    character(*),intent(in)::tag
    integer,intent(in)::south,north,nobs
    type(raw_obs),intent(in)::obs(nobs)
!
!  Local variables
!
    integer::info,j,nobsout
    character(4)::ftype
    real(kind(obs(1)%phi))::y
!
    open(unit=OUTNOSWAP,file=ofile,form='unformatted',action='write', &
       position='rewind',iostat=info)  ! overwrite any predecessor
    if(info.ne.0) call fexit(trim(ofile)//': cannot open',errcode=2)
!
    nobsout=count(obs%phi.ge.RLATS(south).and.obs%phi.le.RLATS(north))
!
    if(verbosity().gt.VERB_TERSE) call finfo('writing '//trim(ofile)//'...')
    ftype=tag
    if(ftype.ne.'none') then
       write(OUTNOSWAP) ftype,seqid(),south,north,nobsout
    else
       write(OUTNOSWAP) seqid()
    endif
!
!  Get the PHI value corresponding to the model grid latitudes defined by
!  SOUTH and NORTH.
!
    do j=1,nobs
       y=obs(j)%phi
       if(y.ge.RLATS(south).and.y.le.RLATS(north)) write(OUTNOSWAP) obs(j)
    enddo
    close(OUTNOSWAP)
    return
    end subroutine write_obs_slice
!----------------------------------------------------------------------------
    end module grid_obs
