    module diagnostics
!
!  SYNOPSIS
!    use diagnostics
!
!  DESCRIPTION
!  Implements diagnostics for Brian's version of the data assimilation method.
!
!  ROUTINES IMPLEMENTED
!   subroutine diag_init(header,diagfn,truthfn,backfn,seq,obsfn,analfn)
!   subroutine output_patch_diag(loc,pd)
!   subroutine compute_rms
!   private subroutine report_rms
!   subroutine summary_diagnostics(nsol,have_natgrid)
!   private subroutine standard_units
!
!  REQUIRED PACKAGES
!   grid_obs - in current directory
!   precision, ferror, io, command_line_support - in support directory
!
!  MODIFICATION HISTORY
!  09/24/05 -  Output results on observations is converted to physical units.
!  08/08/05 -  Added optional reporting on the singular values of the
!              matrix (k-1)I+Eb^T*inv(R)*Eb.
!  07/20/05 -  Replaced "BV" field of diagnostics with "EDIM" to reflect
!              current terminology.
!  07/16/05 -  Updated list of error codes from assimilation method;
!              adjusted entries of ERROR_REASON type.
!  06/26/05 -  Added reporting for enhanced variance inflation.
!  06/18/05 -  Consolidated observation diagnostics back here.  Added
!              interpolation to nature grid so we can compute obs minus truth.
!  01/21/05  - Moved observation diagnostics to grid_obs.
!              Diagnostic output no longer is appended to previous log files
!              of the same name; instead, all previous content is overwritten.
!  12/16/04  - bug fixes to 11/19/04 version.
!  11/19/04  - modifications to deal with new grid i/o formats.  Still need
!              a way to select particular observations and report the
!              interpolated model grid values there.
!  06/22/04  - add computation and output of maximum deviation from truth
!              at patch center.
!  06/19/04  - replace 'cond' with 'rcond' everywhere, since we're computing
!              the reciprocal condition number.  Also print out ensemble
!              spread and additional data when rcond < 1e-4.
!  06/17/04  - support for both Brian's method and the original scheme.
!  06/07/04  - modified to work with Brian's assimilation scheme.
!              Removed code to compute explained variance (since we don't
!              have a Q matrix any more).  BV dimension is computed using
!              the eigenvalues of P^a (which we need anyway to get the
!              symmetric square root).
!  04/19/04  - added output of background ensemble mean of p,t,u,v
!  12/05/03  - moved physical scalings to separate module
!  11/11/03  - revision to supply conversion coefficients to physical units.
!  10/19/03  - revision to handle vertical localization and to output
!              traces of background and analysis covariance matrices.
!  08/31/03  - revision to handle latitudinal strips and patches at poles.
!  08/11/03  - no longer require observations at every point
!            and add 'es' lines to estimate ensemble spread.
!
    use grid_obs
    use command_line_support
    implicit none
!
    character(*),parameter,private::VERSION='09/24/05'
!
!  The "true" model variables from the nature run, if available;
!  natgrid is not associated otherwise.
!
    real(D),allocatable,save,private::natgrid(:,:,:,:) 
                                           ! (south:north,LEVS,LONF,NVAR)
!
!  Determine which diagnostics we're going to compute.
!  By default, we do NOT compute and display ensemble statistics for each
!  observation.
!
    logical,save,private::obs_statistics=.false.
!
    type patch_diag
       real(D)::edim=0.0  ! E-dimension of the patch
       real(D)::rcondpb=0.0  ! recip condition numbe of Pb~ (Ed only)
       real(D)::rcondice=0.0  ! recip cond number of (I+C*Eb) (both methods)
       real(D),allocatable::iceb_sing(:)  ! singular vals of (I+C*Eb) (both)
       real(D)::smax=0.0  ! largest eigenval of F in enhanced inflation
       integer::effrank=0  ! effective rank of F in enhanced inflation
       integer::info=0  ! valid return codes are nonnegative
       integer::nobs=0  ! number of observations in the patch
       logical::brian_failed=.false.  ! true if Brian's method failed
    end type patch_diag
!
!  Error return values for assim()
!
!    INFLATION_FAIL :  enhanced variance inflation failed.
!    SYM_SQRT_FAIL :  the computation of the matrix Y failed (singular)
!    PA_NO_EV      :  computation of eigenvectors of the covar mtx failed
!    P_NEG_EV     :   covariance matrix has nontrivial negative eigenvalue
!    NO_OBS_UPDATE :  No observations available to complete the update
!    INFLATE_FAIL  :  Enhanced variance inflation failed
!    P_SINGULAR    :  Either Pb or Pa is numerically singular
!
    integer,parameter::INFLATION_FAIL=1
    integer,parameter::SYM_SQRT_FAIL=2
    integer,parameter::PA_NO_EV=3
    integer,parameter::P_NEG_EV=4
    integer,parameter::NO_OBS_UPDATE=5
    integer,parameter::PA_SINGULAR=6
    integer,parameter::GARBAGE_VALUES=7
    integer,parameter::PB_NO_EV=8
    integer,parameter::OTHER_REASON=9
    character(26),parameter::errcodes(OTHER_REASON)=(/ &
       "enhanced inflation failed ", &
       "can't compute sym sqrt    ", &
       "can't compute analysis Pa ", &
       "covar matrix asymmetric   ", &
       "no observations in patch  ", &
       "singular covariance matrix", &
       "svd of P^b failed         ", & 
       "method generates garbage  ", &
       "other unspecified problem "/)
!
!  Local routines for reporting results in physical units
!
    private::standard_units,report_rms
    contains
! -----------------------------------------------------------------------------
    subroutine diag_init(header,diagfn,truthfn,backfn,seq,obsfn,analfn)
!  DIAG_INIT - initialize diagnostics package for the data
!  assimilation algorithm.
!  HEADER : 4-character string to expect when reading the "true" atmospheric
!           state file (if any).  Use 'none' of no header is expected.
!  DIAGFN : name of file to which diagnostic output should go.
!  TRUTHFN : name of file (if any) containing the "true" atmospheric state
!  BACKFN  : name of the file (if any) containing the background.
!  SEQ     : sequence id (if any) of the background file
!  OBSFN   : name of the file (if any) containing the observations
!  ANALFN  : name of the file (if any) containing the analysis
!
!  Side effects/requirements:
!  - READIN must have been called first
!  - kills the program if problems arise trying to read the "true" grid file.
!  - opens DIAGUNIT and prints info on background, truth, and analysis files
!
!  NOTE:  this interface is incomplete.  Ultimately we should have
!  switches to control how much information about the observations is
!  written (e.g., 'none', 'all', or 'rejected') and whether all the singular
!  values of I+Eb^T*inv(R)*Eb are written out.  For the moment we'll just
!  determine in MAINLOOP whether the singular values are written out,
!  and we'll dump the data on all the observations.
!
    use grid_io
    implicit none
    character(*),intent(in)::header,diagfn
    character(*),intent(in)::truthfn,backfn,obsfn,analfn
    integer,intent(in)::seq
!
!  Local variables
!  PSOUTH, PNORTH define latitude strip at which assimilation patches can
!  be centered (not important here for the nature grid but are checked
!  for consistency).  INORTH, ISOUTH define the latitude strip containing
!  the nature grid; they must coincide with the latitude strip for each
!  ensemble model grid.
!
    if(len_trim(diagfn).eq.0) &
       call fexit('diag_init: no diagnostic output file specified',2)
    open(unit=DIAGUNIT,file=diagfn,action='write')
!
!  Get the "true" nature (sub)grid and check for compatible boundaries
!
    write(DIAGUNIT,100) VERSION
100 format('#diagnostics package: version ',a)
    if(len_trim(truthfn).gt.0) call read_nature_grid(header,truthfn,seq)
101 format('#',a,1x,a)
    if(len_trim(backfn).gt.0) write(DIAGUNIT,101) 'background',trim(backfn)
    write(DIAGUNIT,102) 'seqid',seq
102 format('#',a,1x,i0)
    if(len_trim(obsfn).gt.0) write(DIAGUNIT,101) 'observations',trim(obsfn)
    if(len_trim(analfn).gt.0) write(DIAGUNIT,101) 'analysis',trim(analfn)
!
!  Explain what the entries mean
!
    write(DIAGUNIT,110) 'loc gives (ix iy iz thickness inflation)'
    write(DIAGUNIT,110) '  where inflation=simple,enhanced,Snyder-Zhang)'
    write(DIAGUNIT,110) 'output values are in physical units'
    write(DIAGUNIT,110) 'patch gives nobs edim rcondice [rcondpb], where'
    write(DIAGUNIT,110) '  nobs is number of observations in the patch'
    write(DIAGUNIT,110) '  edim is E-dimension over entire patch'
    write(DIAGUNIT,110) '  rcondice is recip cond number of I+Eb^T*inv(R)*Eb'
    write(DIAGUNIT,110) '  rcondpb is reciprocal condition number of Pb^'
    write(DIAGUNIT,110) 'SV is singular vals of I+Eb^T*inv(R)*Eb'
110 format('#',a)
    return
    end subroutine diag_init
! -----------------------------------------------------------------------------
    subroutine output_patch_diag(loc,pd)
!  OUTPUT_DIAGS - print out the diagnostics.
!  Side effects: writes to DIAGUNIT
!
    type(patch_index),intent(in)::loc
    type(patch_diag),intent(in)::pd
!
    integer::low,high
    type(inflation_info)::inflate
!
!  If there are no observations, then no useful diagnostics are computed.
!
    if(pd%nobs.eq.0) return
!
!  If there was an error in the data assimilation step other than no
!  observations, display it here.
!
    call zdiam(loc,low,high)
    call inflation(loc,inflate)
    write(DIAGUNIT,100) loc,high-low+1,inflate
100 format('loc ',4i4,3f9.3)
    if(pd%brian_failed) then  ! we used the original formulation
       write(DIAGUNIT,101) pd%nobs,pd%edim,pd%rcondice,pd%rcondpb
    else
       write(DIAGUNIT,101) pd%nobs,pd%edim,pd%rcondice
    endif
101 format('p ',i4,f9.3,2es11.3)
    if(pd%info.ne.0) write(DIAGUNIT,102) 'failed:',errcodes(pd%info)
102 format(a,1x,a)
    if(pd%effrank.gt.0) write(DIAGUNIT,103) pd%effrank,pd%smax
103 format('F',i4,es11.3)
    if(allocated(pd%iceb_sing)) write(DIAGUNIT,104) pd%iceb_sing
104 format('SV',7es11.4/('> ',7es11.4))
    return
    end subroutine output_patch_diag
! ----------------------------------------------------------------------
    subroutine compute_rms
!  COMPUTE_RMS - compute the sum of squares of the difference
!  between the analysis grid (AGRID) and the nature grid over the strip
!  processed by this assimilation run at selected atmospheric levels.
!  (The surface pressure is accessed only at level 1.)
!  This routine assumes that the nature grid is available!  Calling it without
!  a nature grid WILL crash the program.
!  Results are computed and displayed to DIAGUNIT.
!
    integer::levlist(3)=(/1,5,15/)
    integer::asouth,anorth
!
    call assimilation_boundary(asouth,anorth)
    write(DIAGUNIT,100)
100 format('#rms var lat lev value')
!
!  Surface pressure first - level 1 only.
!
    call report_rms(asouth,anorth,levlist(1:1),PRES_TYPE)
!
!  Next temperature at the selected atmospheric levels
!
    call report_rms(asouth,anorth,levlist,TEMP_TYPE)
!
!  Finally, the U component of the wind.  We won't repeat the computation
!  for V since I expect the results to be comparable for both.
!
    call report_rms(asouth,anorth,levlist,UWIND_TYPE)
    return
    end subroutine compute_rms
!----------------------------------------------------------------------------
    subroutine report_rms(asouth,anorth,levlist,typ)
!  REPORT_RMS - compute and report the RMS between the analysis and the truth
!  for the specified variable TYP between the indicated latitudes for the
!  specified model atmospheric levels.
!  Results go to DIAGUNIT.
!
    integer,intent(in)::asouth,anorth,typ,levlist(:)
!
    real(D)::aval(1),natval(1)  ! must be arrays, despite logical scalar quality
    real(D)::rms(asouth:anorth,size(levlist))
    real::v
    type(patch_index)::loc
    integer::ix,iy,k,lev
!
    rms=0.0
    do ix=1,LONF
       do k=1,size(levlist)
          lev=levlist(k)
          do iy=asouth,anorth
             loc=patch_index(ix,iy,lev)
             call copy_grid_point(loc,aval,typ,'amean')
             call copy_grid_point(loc,natval,typ,'nature')
             v=real(aval(1),R)
             call standard_units(v,typ,real(lev))
             aval(1)=v
             v=real(natval(1),R)
             call standard_units(v,typ,real(lev))
             natval(1)=v
             rms(iy,k)=rms(iy,k)+(aval(1)-natval(1))**2
          enddo
       enddo
    enddo
    rms=sqrt(rms/LONF)
    do k=1,size(levlist)
       lev=levlist(k)
       do iy=asouth,anorth
          write(DIAGUNIT,100) TYPEOF(typ),iy,lev,rms(iy,k)
       enddo
    enddo
    return
100 format('rms ',a,2i4,es11.3)
    end subroutine report_rms
!---------------------------------------------------------------------------
    subroutine summary_diagnostics(nsol,have_natgrid)
!  SUMMARY_DIAGNOSTICS - compute and display statistics for each
!  observation that appears in the final analysis grid.  Output goes
!  to the diagnostics file.   We report on all observations that are
!  of types that are assimilated.  (For example, if we don't assimilate
!  surface pressure, then no statistics about pressure observations and
!  the associated interpolated model grid points are output.)
!  Lines beginning with 'o' are observations that pass all applicable
!  consistency checks; 'r' are those that fail some consistency test.
!  IN ADDITION, we compute the global RMS at selected atmospheric levels.
!
!  HAVE_NATGRID : if true, then the nature grid is available for comparison.
!
    implicit none
!
    integer,intent(in)::nsol  ! needed for analysis sdev computation
    logical,intent(in)::have_natgrid
!
    type(ensemble_stat)::bgstat  ! background statistics
    type(ensemble_stat)::anstat  ! analysis statistics
    integer::asouth,anorth,k,total_kept,total_obs,ierr
    character::mnemonic
    real::phisouth,phinorth,v,s,zmean
    real(D)::hx(nsol)
    type(model_obs)::obs
!
    if(have_natgrid) call compute_rms
!
    call assimilation_boundary(asouth,anorth)
    if(asouth.le.1) then  ! assimilation goes to south pole
       phisouth=-90.0
    else
       phisouth=RLATS(asouth)
    endif
    if(anorth.ge.LATG) then  ! assimilation goes to north pole
       phinorth=90.0
    else
       phinorth=RLATS(anorth)
    endif
!
!  Statistics on retained observations
!
    write(DIAGUNIT,100) 'lambda','phi','z','value','sdev', &
        'back mean','back sdev','anal mean','anal sdev','nature val'
100 format('#o|r',2x,2a8,a7,9a11)
!
    call obs_count(total_kept,total_obs)
    do k=1,total_kept
       call kth_obs(k,obs,ierr,bgstat,zmean)
       if(ierr.eq.0.and.obs%phi.ge.phisouth.and.obs%phi.lt.phinorth) then
          call linear_interp(k,nsol,hx,'analysis',anstat)
          call standard_units(obs%value,obs%typ,zmean,obs%sdev)
          v=real(bgstat%mean,R)
          s=real(bgstat%sdev,R)
          call standard_units(v,obs%typ,zmean,s)
          bgstat%mean=v
          bgstat%sdev=s
          v=real(anstat%mean,R)
          s=real(anstat%sdev,R)
          call standard_units(v,obs%typ,zmean,s)
          anstat%mean=v
          anstat%sdev=s
          mnemonic=TYPEOF(obs%typ)
          if(have_natgrid) then
             call linear_interp(k,1,hx,'nature')
             v=real(hx(1),R)
             call standard_units(v,obs%typ,zmean)
             hx(1)=v
             write(DIAGUNIT,110) 'o',obs%typecode,mnemonic, &  ! types
                obs%lambda,obs%phi,obs%x_grd(3), &  ! location
                obs%value,obs%sdev,bgstat%mean,bgstat%sdev, &
                anstat%mean,anstat%sdev,hx(1)
          else
             write(DIAGUNIT,110) 'o',obs%typecode,mnemonic, &  ! types
                obs%lambda,obs%phi,obs%x_grd(3), &  ! location
                obs%value,obs%sdev,bgstat%mean,bgstat%sdev, &
                anstat%mean,anstat%sdev
          endif
       endif
    enddo
110 format(a1,i4,a2,2f8.3,f7.3,7es11.3)
!
!  Statistics on rejected observations
!
    do k=total_kept+1,total_obs
       call kth_obs(k,obs,ierr,bgstat,zmean)
       if(ierr.eq.0.and.obs%phi.ge.phisouth.and.obs%phi.lt.phinorth) then
          call standard_units(obs%value,obs%typ,obs%x_grd(3),obs%sdev)
          v=real(bgstat%mean,R)
          s=real(bgstat%sdev,R)
          call standard_units(v,obs%typ,zmean,s)
          bgstat%mean=v
          bgstat%sdev=s
          mnemonic=TYPEOF(obs%typ)
          if(have_natgrid) then
             call linear_interp(k,1,hx,'nature')
             v=real(hx(1),R)
             call standard_units(v,obs%typ,zmean)
             write(DIAGUNIT,110) 'r',obs%typecode,mnemonic, &  ! types
                obs%lambda,obs%phi,obs%x_grd(3), &  ! location
                obs%value,obs%sdev,bgstat%mean,bgstat%sdev,v
          else
             write(DIAGUNIT,110) 'r',obs%typecode,mnemonic, &  ! types
                obs%lambda,obs%phi,obs%x_grd(3), &  ! location
                obs%value,obs%sdev,bgstat%mean,bgstat%sdev
          endif
       endif
    enddo
    return
    end subroutine summary_diagnostics
! ----------------------------------------------------------------------
    pure subroutine standard_units(val,typ,z,sdev)
!  STANDARD_UNITS - return the relevant value of the arguments in
!  "standard" physical units; see spectral.f90 for the definition of 
!  "standard".
!
    real(R),intent(inout)::val
    integer,intent(in)::typ
    real(R),intent(in)::z
    real(R),intent(inout),optional::sdev
!
    select case(typ)
    case(PRES_TYPE)  ! ln(dPa) -> hPa; Taylor gives log(v/10+s)~log(v/10)+10s/v
       val=10.0*exp(energy_to_pres(val))
       if(present(sdev)) sdev=energy_to_pres(sdev)*val/10.0
    case(TEMP_TYPE)
       val=energy_to_temp(val,z)
       if(present(sdev)) sdev=energy_to_temp(sdev,z)
    case(UWIND_TYPE,VWIND_TYPE)
       val=energy_to_wind(val,z)
       if(present(sdev)) sdev=energy_to_wind(sdev,z)
    case default
       val=0.0
       if(present(sdev)) sdev=0.0
    endselect
    return
    end subroutine standard_units
! ----------------------------------------------------------------------
    end module diagnostics
