module obs_hx
!  OBS_HX - H operator and supporting routines.
!  
!  SYNOPSIS
!    use obs_hx
!
!  DESCRIPTION
!    OBS_HX contains support for handling "real" observations.
!    The code here computes H(x) for a subset of the available ensemble
!    solutions.
!
!  PUBLIC ROUTINES IMPLEMENTED
!    subroutine compute_hx(obslist,nobs,nsol,igrid,oint)
!
!  REQUIRED PACKAGES/!    obs_types, interpolate, grid_model - in main source directory
!
!  MODIFICATION HISTORY
!  10/19/06 - New implementation of the model linear interpolation code
!    to handle reduced-resolution grids.
!  09/06/06 - Version for MPI-based code that handles a global ensemble grid
!    for a limited number of ensemble solutions.
!
!  IMPLEMENTATION NOTES
!  09/06/06 - This version does not support the SSI surface pressure fixup,
!    as the mean background temperature is required.  Here we are computing
!    H(x) for only a subset of the ensemble members, so we don't have the mean.
!
!  AUTHOR
!     Eric Kostelich, kostelich@asu.edu, 480-965-5006.
!
    use interpolate
    use obs_types
    use grid_model
    use io
    use resolution
!    use surface
    implicit none

    integer,parameter,private::QUIKSCAT=285


!
    private::find_z,sigma_fixup,sp_fixup_whitaker,h_operator, &
       four_d,three_d,two_d,one_d
    contains
!----------------------------------------------------------------------------
    subroutine compute_hx(obslist,nobs,nsol,igrid,oint)
!  COMPUTE_HX - Compute H(x) for each observation in OBSLIST.
!  OBSLIST :=: On input, the list of NOBS available observations.  On return,
!     some surface pressure measurements may have been moved to the
!     model surface; all others are unchanged.
!  NSOL : number of ensemble solutions in the model grid.
!  IGRID : grid identifier.  We assume that the geopotential height is available
!    from it.
!  OINT := on return, all components have been allocated, as follows:
!     HX(:,K) holds the computed value of H(x) for the Kth
!     observation, one value for each available ensemble solution.
!     ZLEV(:,K) holds the sigma levels for the Kth observation,
!     one per ensemble solution.
!     KEEP(K) is the number of ensemble members for which
!     the Kth observation satisfies criteria for retention.  (So if KEEP(K) <
!     NSOL, we will throw it out - but KEEP(K)=NSOL does not assure
!     retention, because in general we see only a limited number of
!     ensemble solutions at this point).
!  Note: SET_OBSERVATION_TYPES must have been called first.
!
    integer,intent(in)::nobs,nsol,igrid
    type(model_obs),intent(inout)::obslist(nobs)
    type(hx_info),intent(out)::oint
!
!  Local variables
!
    integer::j,mode
!
    if(.not.all(valid_observation_type(obslist%ityp))) &
       call fabort('compute_hx: unknown observation type in obslist')
    oint%nobs=nobs
    allocate(oint%hx(nsol,nobs),oint%zlev(nsol,nobs),oint%keep(nobs))
!
!  If we're running in 3-d mode, then we project all the observations onto the
!  analysis time.
!
    call atime_limit(mode=mode)
    if(mode.eq.3) obslist%time_diff=0.0
!
!  SIGMA_FIXUP sets the vertical coordinate of each observation 
!  relative to the surface pressure of each ensemble member.
!  Each observation is retained only if the vertical level associated 
!  with each ensemble member is between 1 and LEVS
!
    do j=1,nobs
!  ML Begin scat test code; July 9,2012  
       if(obslist(j)%ncepcode.eq.QUIKSCAT) &         
          write(LOGFILE,*) 'in compute_hx: scat debug, just before call to sigma_fixup, j = ',j
!  ML End scat test code  
       call sigma_fixup(obslist(j),nsol,igrid,oint%zlev(:,j),oint%keep(j))
       if(oint%keep(j).eq.nsol) &
          call h_operator(obslist(j),nsol,igrid,oint%zlev(:,j),oint%hx(:,j))
       if(any(isnan(oint%hx(:,j)))) then  ! debug
          call fabort('NaN in H(x) at observation '//j)
       endif
    enddo
    return
    end subroutine compute_hx
!----------------------------------------------------------------------------
    subroutine h_operator(obs,nsol,igrid,zlev,hx)
!  H_OPERATOR - interpolate the model grid to the observation OBS.
!  by linear interpolation from the nearest model grid points.
!  OBS :=: the observation (in MODEL_OBS format).  It must be one that we
!     have not obviously rejected (i.e., ZLEV must be defined and sensible)
!     and it must be one that we know how to handle.  Units of measurement
!     may be adjusted to conform with background grid units.
!  NSOL  : number of ensemble solutions in GRID.
!  GRID : the background model grid to use for interpolation.
!  ZLEV :  previously computed model Z values from SIGMA_FIXUP.
!  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
!  For the moment, this routine can handle only conventional observations
!  (surface pressure, temperature, wind).  We will extend it later to
!  handle satellite radiances and other observation types.
!
    type(model_obs),intent(inout)::obs
    integer,intent(in)::nsol,igrid
    real,intent(in)::zlev(nsol)
    real(R),intent(out)::hx(nsol)
!
!  Local variables
!
    real(R)::tv1(nsol),w1(nsol),tv(nsol)! model virtual temperature and mixing ratio
    real(R)::u1(nsol),v1(nsol),t1(nsol) ! model U/V and T at ZLEV=1
    real(R)::prs1(nsol), speed1(nsol)   ! model pressure & wind speed at ZLEV=1
    real(R)::pinterp(nsol),rq(nsol)           ! model surface pressure
    real(R)::dp1(nsol), zlev1(nsol)           ! local workspace
    real(R)::dz1(nsol), z0(nsol), z1(nsol)    ! for reduction of
    real(R)::speed10(nsol)                    ! lowest sigma-level
    real(R)::pave(nsol), rhoave(nsol)         ! model winds to
    real(R)::theta1(nsol), thetaV1(nsol)
    real(R)::expon(nsol),fracspd(nsol)
    real(R)::p0                               ! 10 meters

    integer::gtyp
    real::hr
  
    real,parameter::C0=3.7       ! empirical parameter for roughness length eqn
    real,parameter::C1=-1.165    ! empirical parameter for roughness length eqn
    real,parameter::A=0.018      ! Charnock constant
    real,parameter::K=0.4        ! von Karman constant 

!
!  The background grid is maintained in Kelvin, so convert the temperature obs.
!
    if(obs%ityp.eq.OB_TEMP) obs%value=obs%value+KELVIN
!
!  If we have a sensible temperature measurement (type codes 130, 131, 133),
!  then interpolate the background virtual temperature to the observation
!  location and convert to sensible temperature.  This requires two 
!  steps: one to get the interpolated virtual temperature and another to get
!  the interpolated mixing ratio.
!
    hr=obs%time_diff
    if(obs%ityp.eq.OB_TEMP.and.any(obs%ncepcode.eq.(/130,131,133/))) then
       call four_d(obs%x,zlev,hr,nsol,igrid,G_TEMP,tv)  ! virtual temperature
       call four_d(obs%x,zlev,hr,nsol,igrid,G_MIXING,rq)  ! mixing ratio
!       rq=rq/1000.0  ! g/kg -> g/g
       rq=rq/(1.0-rq) ! g/kg -> g/g  !should be corrected but it is not high impact, leav it!
       hx=tv*(rq+1.0)/(rq/0.622+1.0)  ! convert virtual to sensible T

!  If we have a scatterometer measurement (type code 285), then reduce the lowest
!  sigma level wind to the scatterometer height of 10 meters above the sea surface.
!
    else if(obs%ncepcode.eq.QUIKSCAT) then
       zlev1(:)=1.0
!
!  Find lowest sigma-level U/V/T and Psfc at obs location & time
!
       call four_d(obs%x,zlev1,hr,nsol,igrid,G_UWIND, u1)
       call four_d(obs%x,zlev1,hr,nsol,igrid,G_VWIND, v1)
       call four_d(obs%x,zlev1,hr,nsol,igrid,G_TEMP,  tv1)
       call four_d(obs%x,zlev1,hr,nsol,igrid,G_SURFP, pinterp)
       call four_d(obs%x,zlev1,hr,nsol,igrid,G_MIXING,w1)
!
!  Convert surface pressure from units of ln(dPa) to Pa.
!  We need the units in Pa for thermodynamic calculations that follow.
!  Convert mixing ratio from g/kg -> kg/kg.
!
       pinterp=1000.0*exp(pinterp)
       w1=w1/1000.0
!
!  Compute lowest sigma level..
!
!  .. wind speed
       speed1=sqrt(u1*u1 + v1*v1)
!
!  .. pressure
       prs1=SIGMA_LEV(1)*pinterp
!
!  Delta pressure across lowest sigma layer (top minus bottom, or
!  zlev=1 minus zlev=0; the result will be a negative number).
       dp1=prs1 - pinterp 
!
!  Lowest sigma layer average pressure.
       pave=(pinterp + prs1)/2.0 
!
!  Compute potential temperature with the Exner function at zlev=1
!
       t1=tv1/(1.0 + 0.61*w1)
       p0=100000.0
       theta1=t1*(p0/pinterp)**(RD/CP)
!
!  Compute virtual potential temperature at zlev=1
!
       thetaV1 = theta1*(1.0 + 0.61*w1)
!
!  Solve for lowest sigma layer average density, rhoave,
!  assuming a well-mixed lowest sigma layer. (In a well-mixed,
!  lowest sigma layer, thetaV is conserved and w is 
!  nearly constant.)  So substitute the eqn of state
!  for Tave, Tave = pave/(rhoave*RD), on the R.H.S of 
!  the thetaV definition and solve for rhoave.
!
       rhoave=(pave/(thetaV1*RD))*(1.0 + 0.61*w1)*(p0/pave)**(RD/CP)

!  Use hydrostatic eqn to determine height of lowest model level.
       dz1=(-1.0)*dp1/(rhoave*GRAV)
!
!  Set height of lowest model sigma level.  
       z1=dz1
!
!  Empirical calculation to determine z0, surface roughness, 
!  from Hoffman (2011; arXiv:1107.1416v1 [physics.ao-ph]).
!  Assumes neutral stability in the boundary layer.
!
       expon=C0 + C1*log10((A*K**2*speed1**2)/(GRAV*z1)) 
       z0=z1**((-1.0)*expon)   
!
!  Reduce lowest sigma-level winds to 10 meters  
       speed10 = speed1*((log10(10.0/z0))/(log10(z1/z0)))
       fracspd = speed10/speed1
!
!  Apply fractional reduction in windspeed to wind components.
!
       select case(obs%ityp)
       case(OB_UWIND)
          hx = u1 * fracspd
       case(OB_VWIND)
          hx = v1 * fracspd
       endselect

       write(LOGFILE,*) 'h_operator: scat debug, pinterp=',pinterp
       write(LOGFILE,*) 'h_operator: scat debug, prs1   =',prs1
       write(LOGFILE,*) 'h_operator: scat debug, w1     =',w1
       write(LOGFILE,*) 'h_operator: scat debug, t1     =',t1
       write(LOGFILE,*) 'h_operator: scat debug, theta1 =',theta1
       write(LOGFILE,*) 'h_operator: scat debug, thetaV1=',thetaV1
       write(LOGFILE,*) 'h_operator: scat debug, t1     =',t1
       write(LOGFILE,*) 'h_operator: scat debug, z1     =',z1
       write(LOGFILE,*) 'h_operator: scat debug, z0     =',z0
       write(LOGFILE,*) 'h_operator: scat debug, speed1 =',speed1
       write(LOGFILE,*) 'h_operator: scat debug, speed10=',speed10
       write(LOGFILE,*) 'h_operator: scat debug, fracspd=',fracspd

!  Otherwise we have a measurement that we need to interpolate according
!  to space and time.
!
    else
       select case(obs%ityp)  ! map the obs type to the grid type
       case(OB_SURFP)
          gtyp=G_SURFP
       case(OB_TEMP)
          gtyp=G_TEMP
       case(OB_UWIND)
          gtyp=G_UWIND
       case(OB_VWIND)
          gtyp=G_VWIND
       case default  ! other obs types can be added here later
          call fabort('h_operator: unknown observation type')
       endselect
       call four_d(obs%x,zlev,hr,nsol,igrid,gtyp,hx)
    endif
    return
    end subroutine h_operator
!----------------------------------------------------------------------------
    subroutine sigma_fixup(obs,nsol,igrid,zlev,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).  ZLEV 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 ZLEV value is outside the interval [1,LEVS].
!  IGRID : model grid identifier.
!  NSOL : number of available ensemble solutions.
!  OBS :=: the current observation (MODEL_OBS format).  The Z component
!     of the position is adjusted to 1 if the obs is surface pressure.
!     For observation types other than the surface pressure, the Z component
!     is set to the mean of the interpolated observation heights.
!  ZLEV  := the interpolated observation height for each ensemble member.
!  KEEP := The number of ensemble members for which 1 <= Z <= LEVS and
!    for which the observation is within the forecast time range.
!
    type(model_obs),intent(inout)::obs
    integer,intent(in)::nsol,igrid
    real,intent(out)::zlev(nsol)
    integer,intent(out)::keep
!
!  Local variables
!
    real::hr,early,late  ! time measurements in hours
    real(R)::pinterp(nsol)

   type(grid_index)::near
    integer,parameter::S_LANDSEA_MASK=1 ! land-sea mask (0=water, 1=land, 2=ice)
    real(R)::surftype(nsol)  ! surface type (land, water, ice, etc.)

!
    real(R),dimension(nsol)::tv1,w1,prs1,pave,t1,theta1,zlev1
    real(R),dimension(nsol)::thetaV1,qsigma,rhoave,dp
    real(R)::p0,dz


!  For measurements not taken at the surface (CAT field nonzero), P_LEV is
!  the pressure level in hPa and must be nonzero.
!
    if(obs%surface.eq.0.and.obs%p_lev.le.0.0) then
       call fabort('sigma_fixup: obs p_lev <= 0')
    endif
    if(obs%ityp.ne.OB_SURFP.and.obs%prss.le.0.0) &
       call fabort('sigma_fixup: obs prss <= 0')
!
!  If the observation precedes or follows any available background forecast,
!  then discard it.  If we're running in 3-d mode, then
!  EARLY and LATE give the observation time limits around the analysis time.
!
    call atime_limit(early,late)
    hr=obs%time_diff   ! time differential of observation in hours
    if(hr.lt.early.or.hr.gt.late) then  ! outside range of forecasts
       keep=0
       return
    endif
!
!  Surface pressure observations.
!
    if(obs%ityp.eq.OB_SURFP) then
       call sp_fixup_whitaker(igrid,obs,keep)
       if(keep.gt.0) then
          zlev=1.0  ! surface pressure is at the model surface by definition
          keep=nsol  ! and it's either OK for all solutions or not
          obs%x(3)=1.0  ! by definition for surface pressure
       endif

  else if(obs%ncepcode.eq.QUIKSCAT) then
!
!  Find the sigma level corresponding to 10 meters from the surface,
!  using the hydrostatic equation and assuming a well-mixed
!  lowest sigma layer, virtual potential temperature is conserved,
!  allowing for the layer mean density to be computed.
!
!  Find lowest sigma-level Tv, Psfc and mixing ratio at 
!  obs location & time.
!
       zlev1(:)=1.0
       call four_d(obs%x,zlev1,hr,nsol,igrid,G_TEMP,  tv1)
       call four_d(obs%x,zlev1,hr,nsol,igrid,G_SURFP, pinterp)
       call four_d(obs%x,zlev1,hr,nsol,igrid,G_MIXING,w1)
!    dimension(nsol)::tv1, pinterp, w1, prs1, pave, t1, theta1, thetaV1, dp
!    real:: p0
!  Convert surface pressure from units of ln(dPa) to Pa.
!  We need the units in Pa for thermodynamic calculations that follow.
!
       pinterp=1000.0*exp(pinterp)
       w1=w1/1000.0
!
!  Compute lowest sigma level pressure.
!
       prs1=SIGMA_LEV(1)*pinterp
!
!  Lowest sigma layer average pressure.
!
       pave=(pinterp + prs1)/2.0 
!
!  Compute potential temperature with the Exner function at zlev=1
!
       t1=tv1/(1.0 + 0.61*w1)
       p0=100000.0
       theta1=t1*(p0/prs1)**(RD/CP)
!
!  Compute virtual potential temperature at zlev=1
!
       thetaV1 = theta1*(1.0 + 0.61*w1)
!
!  Solve for lowest sigma layer average density, rhoave,
!  assuming a well-mixed lowest sigma layer. (In a well-mixed,
!  lowest sigma layer, thetaV is conserved and w is 
!  nearly constant.)  So substitute the eqn of state
!  for Tave on the R.H.S of the thetaV definition, i.e., 
!  T = pave/(rhoave*RD), and solve for rhoave.
!
       rhoave=(pave/(thetaV1*RD))*(1.0 + 0.61*w1)*(p0/pave)**(RD/CP)

!  Use hydrostatic eqn to determine the delta pressure from the 
!  surface to 10 meters above the surface.
!
       dz=10.0
       dp=(-1.0)*dz*rhoave*GRAV
!
!   Compute the sigma value, given dp.
!
       qsigma = (dp + pinterp)/pinterp
!
!  Compute sigma index.
!
       zlev = (1.0-qsigma)/(1.0-SIGMA_LEV(1))
!
!  Assign observation z location to average sigma index.
!
       obs%x(3)=sum(zlev)/nsol

       obs%sdev=1.8 ! assign scat obs expected error (valid for scalar windspeed 3-20 m/s)

!  This QUIKSCAT part is taken from Eric's latest global code of Nov. 2010 (4d   
!       if(obs%ityp.eq.OB_UWIND.or.obs%ityp.eq.OB_VWIND) then
!          near=nearest_grid_index(obs%x)
!          call from_sfc_grid(near,S_LANDSEA_MASK,nsol,surftype)
!          keep=merge(nsol,0,surftype(1).eq.0.0)  ! water
!       else
!          zlev=0.0
!          keep=0  ! Quikscat measures wind only, as far as I know.
!       endif`

       keep=nsol

       write(LOGFILE,*) 'in sigma_fixup: scat debug,    zlev = ', zlev
       write(LOGFILE,*) 'x value prss p_lev surface =', &
            obs%x,obs%value,obs%prss,obs%p_lev,obs%surface




    else if(obs%surface.ne.0) then  ! skip other surface obs for now
       keep=0

!  The following part before *** is also taken from Eric's up to date code
! It seems a good idea to include this with the QUIKSCAT data

!  ****************
!  We don't assimilate the wind above the highest latitude line (i.e.,
!  within the polar circles).  I'm not convinced that the interpolation
!  routines handle it correctly, and in any case it will be
!  converted to pseudo-wind for the model (and multiplied by the cosine
!  of the latitude).
!
! Although Dasa later in this code or earlier? says that we have no polar
! problem with the regional


!    else if((obs%ityp.eq.OB_UWIND.or.obs%ityp.eq.OB_VWIND).and. &
!      abs(obs%x(2)).gt.rlats(LATG)) then  ! inside polar circle
!       zlev=0.0
!       keep=0



   !**********




    else
!
!  Interpolate the background surface pressure (in ln(dPa)) for all other types.
!  We need the units in hPa for computing the Z levels.
!  ZLEV is undefined but unused in this particular call to FOUR_D.
!
       call four_d(obs%x,zlev,hr,nsol,igrid,G_SURFP,pinterp)
       pinterp=10.0*exp(pinterp)
       call find_z(obs,pinterp,zlev,nsol)
!  ML Begin scat test code; July 9, 2012
!       if(obs%ncepcode.eq.QUIKSCAT) then
!          write(LOGFILE,*) 'in sigma_fixup: scat debug,    zlev = ', zlev
!          write(LOGFILE,*) 'in sigma_fixup: scat debug, pinterp = ', pinterp
!          write(LOGFILE,*) 'x value prss p_lev surface = ',obs%x,obs%value,obs%prss,obs%p_lev,obs%surface
!       endif
!  ML End scat test code
       keep=count(zlev.ge.1.0.and.zlev.le.LEVS)  !original filter
!       keep=count(zlev.gt.0.0.and.zlev.le.LEVS) !scat test filter
!if(obs%ityp.eq.OB_VWIND) print *,obs%x,obs%value,obs%prss,obs%p_lev
    endif
    return
    end subroutine sigma_fixup
!----------------------------------------------------------------------------
    subroutine sp_fixup_whitaker(igrid,obs,keep)
!  SP_FIXUP_WHITAKER - adjust the value of a surface pressure observation to
!  the model surface using Jeff Whitaker's simplied approach.
!  The standard deviation of the observational error is also modified.
!  Arguments:
!  IGRID : model grid identifier (only the geopotential height is needed).
!  OBS  :=:  the surface pressure observation, which is in ln(dPa).
!     The observation value and associated surface temperature are adjusted
!     as described below.
!  KEEP  := set to 1 if the observation is to be retained and 0 otherwise.
!  Note: no model bias correction is applicable here.
!
    integer,intent(in)::igrid
    type(model_obs),intent(inout)::obs
    integer,intent(out)::keep
!
!  Local variables
!
    real::prssk
    real(R)::height(1),gzm,dgz,tml
!
    real,parameter::GRAD=-0.0065  ! Standard temperature gradient
!
!  We'll generate garbage if the surface pressure isn't in ln(dPa).
!
    if(obs%value.gt.5.0.or.obs%value.lt.3.5) &
       call fabort('sp_fixup_whitaker: obs value out of range')
!
!  Interpolate (in 2D) the model orography to the observation location.
!  Construct a fake observation for this purpose.
!
    call two_d(obs%x,1,0,igrid,G_HEIGHT,height,lev=1)
    gzm=height(1)
!
!  Calculate the difference between the heights of the real and model
!  orographies.  When dgz is positive (negative), the observation is above
!  (below) the model surface.
!
    dgz=obs%p_lev-gzm
    if(abs(dgz).gt.600.0) then  ! discard obs if the change is more than 600 m
       keep=0
       return
    else if(abs(dgz).lt.10.0) then  ! do nothing if change is small enough
       keep=1
       return
    endif
!
!  Calculate the distance between the first sigma level and the model surface.
!
    prssk=obs%prss+KELVIN  ! convert temperature to Kelvin
!
!  TML is the mean temperature of the layer, for use in the hydrostatic 
!  balance equation.
!
    tml=prssk+(0.5*GRAD)*dgz
!
!  Adjust the observed pressure value (which is ln(dPa)) and its standard
!  error.  The factor of 0.003 assumes that there is an error of 3 K/km
!  in the lapse rate.
!
    obs%value=real(obs%value+(GRAV/RD)*dgz/tml,kind(obs%value))
    obs%sdev=real(obs%sdev+(GRAV/RD)*(dgz/prssk**2)*(0.5*dgz*0.003_D), &
       kind(obs%sdev))
    keep=1
    return
    end subroutine sp_fixup_whitaker
!----------------------------------------------------------------------------
    subroutine find_z(obs,pinterp,z,nsol)
!  FIND_Z - find the Z level for each ensemble solution, given the
!  observation pressure level and previously interpolated ensemble surface
!  pressures.  Utility routine.
!  OBS :=: the observation.  The third component of X_GRD is computed.
!  PINTERP : Kth component is the surface pressure (in hPa) 
!     of the Kth ensemble solution interpolated to the observation location.
!  Z  := Kth component is the Z level for the Kth ensemble solution.
!  NSOL : number of ensemble solutions.
!
    type(model_obs),intent(inout)::obs
    integer,intent(in)::nsol
    real(R),intent(in)::pinterp(nsol)
    real,intent(out)::z(nsol)
!
    real::sigma(nsol)
    integer::j
!

! if(obs%nceptype.eq.QUIKSCAT) then 
!    TBD
!    calculate obs p_lev for each ensemble member given z=10m
!    sigma = real(my_p_lev/pinterp,kind(sigma))
! else
    sigma=real(obs%p_lev/pinterp,kind(sigma))
! endif
    if(any(sigma.le.0.0)) call fabort('find_z: negative SIGMA value')
!
    do j=1,nsol
       z(j)=index_interp(log(sigma(j)),LOG_SIGMA_LEV,LEVS)
    enddo
    obs%x(3)=sum(z)/nsol
    return
    end subroutine find_z
!----------------------------------------------------------------------------
    subroutine four_d(x,zlev,t,nsol,igrid,gtyp,value)
!  FOUR_D - internal work routine to perform linear space and time
!  interpolation.  All variables, including pressure, are interpolated in
!  their "standard" model units (i.e., we interpolate ln(dPa) for pressure).
!  X  :  (longitude, latitude) coord to interpolate to (in degrees).
!  ZLEV  : Kth entry is the model vertical level of the observation relative
!     to the model grid of the Kth ensemble solution.
!  T  : observation time in hours relative to the assimilation time.
!  NSOL  : the number of ensemble solutions.
!  IGRID  : the model grid to use.
!  GTYP  : model grid variable type to be interpolated.
!  VALUE := the interpolated values; the Kth component is the space/time
!     interpolation of the Kth ensemble solution.
!  Internal work routine.
!  Algorithm:  Determine the times T1 and T2 that bracket the observation
!  time (MINUTES).  Perform 3-d interpolation at each of T1 and T2,
!  then linearly interpolate to the time T defined by MINUTES.
!
    integer,intent(in)::nsol,igrid,gtyp
    real,intent(in)::x(2),t,zlev(nsol)
    real(R),intent(out)::value(nsol)
!
!  Local variables.
!
    real(R)::vlo(nsol),vhi(nsol)
    real::hr,dt,tlo,thi
    integer::lo,hi,mode,ierr
!
!  If we're running in 3-d mode, then we copy background values from the only 
!  available background grid, which is at the analysis time.  The special
!  index value LO=HI=0 indicates this grid (regardless of mode).
!  If LO != HI, then we need to do a time interpolation, so we must copy
!  the relevant variables from two different grids.
!
    call atime_limit(mode=mode)
    if(mode.eq.3) then
       lo=0
       hi=0
       hr=0.0
    else
       hr=t
       call time_index(hr,lo,hi,tlo,thi,ierr)
       if(ierr.ne.0) call fabort('four_d: obs time out of range')
    endif
!
!  Interpolate the background as needed to the vertical height of 
!  the observation.
!
    if(gtyp.eq.G_HEIGHT.or.gtyp.eq.G_SURFP) then
       if(mode.eq.3.or.t.eq.0.0) then
          call two_d(x,nsol,0,igrid,gtyp,value,lev=1)
       else if(lo.eq.hi) then  ! observation time is exactly at a forecast time
          call two_d(x,nsol,lo,igrid,gtyp,value,lev=1)
       else
          call two_d(x,nsol,lo,igrid,gtyp,vlo,lev=1)
          call two_d(x,nsol,hi,igrid,gtyp,vhi,lev=1)
          dt=(hr-tlo)/(thi-tlo)
          value=(1.0-dt)*vlo+dt*vhi
       endif
    else
       if(mode.eq.3.or.t.eq.0.0) then
          call three_d(x,nsol,0,igrid,gtyp,zlev,value)
       else if(lo.eq.hi) then  ! observation time is exactly at a forecast time
          call three_d(x,nsol,lo,igrid,gtyp,zlev,value)
       else
          call three_d(x,nsol,lo,igrid,gtyp,zlev,vlo)
          call three_d(x,nsol,hi,igrid,gtyp,zlev,vhi)
          dt=(hr-tlo)/(thi-tlo)
          value=(1.0-dt)*vlo+dt*vhi
       endif
    endif
    return
    end subroutine four_d
!----------------------------------------------------------------------------
    subroutine three_d(x,nsol,idt,igrid,gtyp,zlev,value)
!  THREE_D - Do a three-dimensional interpolation.
!  X  :  (longitude, latitude) coord to interpolate to (in degrees).
!  NSOL  : the number of ensemble solutions.
!  IDT  : time identifier for the grid.
!  IGRID  : the model grid to use.
!  GTYP  : model grid variable type to be interpolated.
!  ZLEV  : Kth entry is the model vertical level of the observation relative
!     to the model grid of the Kth ensemble solution.
!  VALUE := the interpolated values; the Kth component is the space/time
!     interpolation of the Kth ensemble solution.
!  Internal work routine.
!  Algorithm:  For each ensemble solution, find the model levels ZLO and ZHI
!  that bracket the observation.  Perform 2-d interpolation at each of the
!  two levels, then linearly interpolate to ZLEV.
!
    integer,intent(in)::nsol,idt,igrid,gtyp
    real,intent(in)::x(2),zlev(nsol)
    real(R),intent(out)::value(nsol)
!
!  Local variables
!
    integer::iz(nsol)
    real::dz(nsol)
    real(R)::zfloor(nsol),zceil(nsol)
!
    iz=max(1,floor(zlev))
    dz=zlev-iz
    call two_d(x,nsol,idt,igrid,gtyp,zfloor,iz)
    iz=min(LEVS,ceiling(zlev))
    call two_d(x,nsol,idt,igrid,gtyp,zceil,iz)
!
!  The vertical spacing between model grid points is 1 by definition, which
!  is why we don't need to normalize DZ.
!
    value=(1.0-dz)*zfloor+dz*zceil
    return
    end subroutine three_d
!----------------------------------------------------------------------------
    subroutine two_d(x,nsol,idt,igrid,gtyp,value,iz,lev)
!  TWO_D - perform a 2-d linear interpolation on the model grid.
!  X := (longitude,latitude) coordinates for which the interpolation is desired.
!  NSOL : number of ensemble solutions.
!  IDT : grid time index.
!  IGRID : grid identifier (G_BACKGROUND, etc.)
!  GTYP : model grid variable type to be interpolated.
!  VALUE := Kth entry holds the Kth interpolated ensemble solution.
!  LEV :  (optional) the vertical model level at which all interpolations
!        are done.
!  IZ  : (optional) Kth entry is the vertical model level at which the
!      interpolation is done for the Kth ensemble member.
!  One of LEV or IZ must be specified.
!  Internal work routine.
!  Algorithm: Find the model latitude values that bracket the observation.
!  Then find the pair of grid points along each latitude line that brackets
!  the longitude of the observation.  Do a 1-d linear interpolation to the
!  longitude point of the observation on each model latitude line.  Then
!  interpolate linearly in the latitudinal direction.
!
    real,intent(in)::x(2)
    integer,intent(in)::nsol,idt,igrid,gtyp
    real(R),intent(out)::value(nsol)
    integer,intent(in),optional::lev,iz(nsol)
!
!  Local variables
!
    real::xopp,dlo,dy,arclen
    real(R)::vlo(nsol),vhi(nsol)
    type(integer_range)::latrange,lonrange
    integer::lat
!    logical::polar
!
    if(x(1).lt.0.0.or.x(1).gt.360.0.or.abs(x(2)).gt.90.0) &
       call fabort('two_d: obs location out of bounds')
!
!  First find the model latitude strip bracketing the observation.
!  in regional model find latitude and longitude range at once
!    call latitude_range(x(2),x(2),latrange)
    call latlon_range(x(1),x(1),x(2),x(2),latrange,lonrange)
!
!  Linearly interpolate between the longitude points on the southern boundary
!  of the box containing the observation.
!
    call one_d(x(1),lonrange,latrange%lo,nsol,idt,igrid,gtyp,vlo,iz,lev)
!
!  Likewise (in the ordinary case) for the northern boundary.
!  If we're near the poles, go to the point 180 degrees opposite; then
!  interpolate along the arc between these two points.
!
!    dlo=rlats(LATG)  ! at T62, DLO ~ 88.5
!   in regional no problems with poles
!    polar=abs(x(2)).gt.dlo  ! inside a polar circle?
!    if(.not.polar) then  ! usual case
       xopp=x(1)
       dlo=rlats(lonrange%lo,latrange%lo)
       arclen=rlats(lonrange%lo,latrange%hi)-dlo
       dy=x(2)-dlo
       lat=latrange%hi
!    else  ! go 180 degrees across through the pole
!       xopp=x(1)+sign(180.0,180.0-x(1))
!       arclen=2.0*(90.0-dlo)
!       dy=abs(x(2))-dlo
!       lat=latrange%lo
!    endif
    call one_d(xopp,lonrange,lat,nsol,idt,igrid,gtyp,vhi,iz,lev)
!
!  If a V-wind observation lies over the pole with respect to the center
!  of the region, then change the sign.  (If, at the observation location, the
!  wind is blowing toward the pole, and at a point opposite the wind is
!  blowing away from the pole at the same speed, then the best interpolated
!  value is one that puts the wind going over the pole in the appropriate
!  direction at the same speed.  If the sign isn't changed, then the
!  interpolated grid value is zero.)  Note that this consideration does not
!  apply to the U-wind; if, say, the wind is blowing parallel to the entire
!  polarmost latitude circle at constant speed, then the best interpolated
!  value at the pole is zero.
!
!    if(polar.and.gtyp.eq.G_VWIND) vhi=sign(vhi,vlo)
    dy=dy/arclen
    value=(1.0-dy)*vlo+dy*vhi
    return
    end subroutine two_d
!----------------------------------------------------------------------------
    subroutine one_d(lambda,lonrange,lat,nsol,idt,igrid,gtyp,value,iz,lev)
!  ONE_D - linear interpolation between two grid points on the
!  latitude circle LAT.
!  LAMBDA : the longitude (in degrees) at which the interpolation is desired.
!  LAT  : integer model latitude index at which interpolation is desired.
!  NSOL : number of ensemble solutions.
!  IDT : grid time index.
!  IGRID : grid identifier (G_BACKGROUND, etc.)
!  GTYP : model grid variable type to be interpolated.
!  VALUE  := Kth entry holds the Kth interpolated ensemble solution.
!  IZ  : (optional) Kth entry is the vertical model level at which the
!      interpolation is done for the Kth ensemble member.
!  LEV :  (optional) the vertical model level at which all interpolations
!        are done.
!  One of LEV or IZ must be specified.
!  Internal work routine.
!  Algorithm: Find the model longitude points on the model latitude LAT
!  that bracket the observation.  Interpolate the grid point values at
!  each endpoint linearly to the longitude LAMBDA.
!  Internal work routine.
!
    real,intent(in)::lambda
    integer,intent(in)::lat,nsol,idt,igrid,gtyp
    real(R),intent(out)::value(nsol)
    integer,intent(in),optional::lev,iz(nsol)
    type(integer_range),intent(in)::lonrange
!
!  Local variables
!
    type(grid_index)::left,right
    real::xleft,xright,dx
    real(R)::vleft(nsol),vright(nsol)
!
!    if(lambda.lt.0.0.or.lambda.gt.360.0.or.lat.lt.1.or.lat.gt.LATG) &
!       print *,'=============  ',lambda,lat
    if(lambda.lt.0.0.or.lambda.gt.360.0.or.lat.lt.1.or.lat.gt.LATG) &
       call fabort('one_d: location out of bounds')
    if(.not.present(lev).and..not.present(iz)) &
       call fabort('one_d: one of LEV or IZ must be given')
!
!  Determine the bracketing model grid points and their locations in degrees.
!
!    call longitude_range(lambda,lambda,lat,lonrange)
    if(present(iz)) then  ! vertical coord comes from IZ
       left=grid_index(lonrange%lo,lat,1)  ! COPY_FROM_GRID ignores the "1"
       right=grid_index(lonrange%hi,lat,1)
    else  ! constant level from LEV
       left=grid_index(lonrange%lo,lat,lev)
       right=grid_index(lonrange%hi,lat,lev)
    endif
    xleft=rlons(left%x,lat)
    xright=rlons(right%x,lat)
    dx=(lambda-xleft)/(xright-xleft)
!
!  Linearly interpolate model endpoint values over the interval.
!
    call copy_from_grid(left,gtyp,idt,igrid,vleft,nsol,iz)
    call copy_from_grid(right,gtyp,idt,igrid,vright,nsol,iz)
    value=real((1.0-dx)*vleft+dx*vright,kind(value))
    return
    end subroutine one_d
!----------------------------------------------------------------------------
    end module obs_hx