module m_gpsStats
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:	 module m_gpsStats
!   prgmmr:	 j guo <jguo@nasa.gov>
!      org:	 NASA/GSFC, Global Modeling and Assimilation Office, 610.3
!     date:	 2015-08-28
!
! abstract: a modular wrapper of (gsp_allhead, gps_alltail), and genstats_gps()
!
! program history log:
!   2007-06-22  cucurull - modify gps_all_ob_type structure
!   2015-08-28  j guo    - created this module on top of genstats_gps();
!                        . completed with type/data components from obsmod;
!                        . changed code where this module needs to be used;
!                        . added this document block;
!                        . for earlier history log, see the history log section
!                          inside ::genstats_gps() below.
!   2016-05-18  j guo    - Made the type private, since this is only a single
!                          instance module object, with its components defined
!                          as module variables (gps_allhead, and gsp_alltail).
!                        . Removed old interface names, which are no longer used
!                          in this version of GSI.
!                        . Edited the in-file documentation.
!
!   input argument list: see Fortran 90 style document below
!
!   output argument list: see Fortran 90 style document below
!
! attributes:
!   language: Fortran 90 and/or above
!   machine:
!
!$$$  end subprogram documentation block

! module interface:

  use m_gpsNode, only: gps_ob_type => gpsNode
  use kinds , only: r_kind,i_kind
  implicit none
  private	! except
  
        ! Data Structure:
  !public:: gps_all_ob_type      ! currently not required to be public

    type gps_all_ob_type
      type(gps_all_ob_type),pointer :: llpoint => NULL()
      type(gps_ob_type),pointer :: mmpoint => NULL()
      real(r_kind)    :: ratio_err                        
      real(r_kind)    :: obserr                        
      real(r_kind)    :: dataerr                       
      real(r_kind)    :: pg                       
      real(r_kind)    :: b                      
      real(r_kind)    :: loc                    
      real(r_kind)    :: type               

      real(r_kind),dimension(:),pointer :: rdiag => NULL()
      integer(i_kind) :: kprof
      logical         :: luse          !  flag indicating if ob is used in pen.

      logical         :: muse          !  flag indicating if ob is used in pen.
      character(8)    :: cdiag

      integer(i_kind) :: idv,iob	      ! device id and obs index for sorting
      real   (r_kind) :: elat, elon      ! earth lat-lon for redistribution
      !real   (r_kind) :: dlat, dlon      ! earth lat-lon for redistribution
    end type gps_all_ob_type

    type gps_all_ob_head
      integer(i_kind):: n_alloc=0
      type(gps_all_ob_type),pointer :: head => NULL()
    end type gps_all_ob_head

        ! Data objects:

  public:: gps_allhead
  public:: gps_alltail

        ! interfaces:

  public:: gpsStats_create
  public:: gpsStats_destroy

        interface gpsStats_create ; module procedure create_; end interface
        interface gpsStats_destroy; module procedure destroy_genstats_gps; end interface

  public:: gpsStats_genStats

        interface gpsStats_genStats; module procedure genstats_gps; end interface

        ! Synopsis:
        !       - []_create: allocated in gsimod::gsimain_initialize().  It was
        !         done through ::create_obsmod_vars(), and now next to it.
        !
        !       - externally built, node-by-node, in setupbend() or setupref().
        !         This is the reason why the ADT can not be defined as "private".
        !
        !       - []_genStats: used to update m_rhs::[ab]work, in setuprhsall(),
        !         through ::genstats_gps();
        !
        !       - []_destroy: deallocated within genstats_gps(), when it it
        !         finished.
        !
        ! The use of []_create()/[]_destroy() pair is obviously not symmetric.  It
        ! would cleaner if they are be moved to the level of setuprhsall(),
        ! where m_rhs::[ab]work are computed.  e.g.,
        !
        !       if(init_pass) call gpsStats_create()
        !       ...
        !       if(last_pass) then
        !         call gpsStats_genstats(bwork,awork,...)
        !         call gpsStats_destroy()
        !       endif
        !
        ! As it is now, the second half has been done, but the first half
        ! stayed at where it has been, for later.

! Most implementations of this module, are snap-shots of gps_all_ob_type, from
! obsmod, its original home.  These implementations include, the type
! definition (gps_all_ob_type and gsp_app_ob_head), instantiation of the type
! (gps_allhead and gsp_alltail), and interfaces for the operations of this
! object ([]_create, []_genstats, []_destroy).  The part of implementation for
! the node-growing, is remained in setupbend() and setupref() as is.

  type(gps_all_ob_head),dimension(:),pointer :: gps_allhead => null()
  type(gps_all_ob_head),dimension(:),pointer :: gps_alltail => null()

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  character(len=*),parameter :: myname='genstats_gps'
contains

  subroutine create_()
    use gsi_4dvar, only: nobs_bins
    implicit none

    ALLOCATE(gps_allhead(nobs_bins))
    ALLOCATE(gps_alltail(nobs_bins))
  end subroutine create_

  subroutine destroy_genstats_gps()
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    destroy_genstats_gps
!     prgmmr:    treadon     org: np20                date: 2005-12-21
!
! abstract:  deallocate arrays holding gps information
!
! program history log:
!   2005-12-21  treadon
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$  end documentation block
    use gsi_4dvar, only: nobs_bins
    use kinds, only: i_kind
    implicit none

    integer(i_kind):: istatus,ii

    do ii=1,nobs_bins
       gps_alltail(ii)%head => gps_allhead(ii)%head
       do while (associated(gps_alltail(ii)%head))
          gps_allhead(ii)%head => gps_alltail(ii)%head%llpoint
          deallocate(gps_alltail(ii)%head,stat=istatus)
          if (istatus/=0) write(6,*)'DESTROY_GENSTATS_GPS: deallocate error for gps_all, istatus=',istatus
          gps_alltail(ii)%head => gps_allhead(ii)%head
       end do
    end do

    return
  end subroutine destroy_genstats_gps

subroutine genstats_gps(bwork,awork,toss_gps_sub,conv_diagsave,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    genstats_gps    generate statistics for gps observations
!   prgmmr: treadon          org: np20                date: 2005-12-21
!
! abstract:  For gps observations, this routine
!              a) collects statistics for runtime diagnostic output
!              f) adjusts observation error ratio based on superobs factor
!
! program history log:
!   2005-12-21  treadon
!   2006-02-24  derber  - modify to take advantage of convinfo module
!   2006-07-28  derber  - modify to use new inner loop obs data structure
!   2006-09-20  cucurull - replace superobs factor for obs in a top (non-full) layer 
!   2007-03-01  treadon - add array toss_gps
!   2007-03-19  tremolet - binning of observations
!   2007-06-21 cucurull - add conv_diagsave and mype in argument list; 
!                         modify qc and output for diagnostic file based on toss_gps
!                         print out diagnostic files if requested
!                         add wgtlim and huge_single in constants module
!   2008-02-27 cucurull - modify diagnostics output file
!   2008-04-14 treadon  - compute super_gps within this routine
!   2008-06-04 safford  - rm unused vars and uses
!   2008-09-05 lueken   - merged ed's changes into q1fy09 code
!   2008-25-08 todling  - adapt obs-binning change to GSI-May2008
!   2009-02-05 cucurull - modify latitude range four statistics output
!   2009-10-22     shen - add high_gps
!   2010-04-09 cucurull - fix several bugs for high_gps (diag information, counters, 
!                       - consider failure of gross check, obs-binning structures, QC for CL profiles)
!                       - reorganize high_gps structure
!                       - modify dimension of diagnostic structure
!   2010-07-23 treadon  - add ratio_error=zero to reqional QC block, replace (izero,ione) with (0,1),
!                         remove _i_kind suffix from integer constants, clean up use statements
!   2010-08-17 treadon  - convert high_gps from m to km one time only; break out regional
!                         QC as separate if/then block (global will bypass); replace 
!                         ratio_errors_reg with logical toss
!   2010-10-25 cucurull - add quality control options for C/NOFS satellite
!   2011-01-18 cucurull - increase the size of nreal and mreal by one element to 
!                         add gps_dtype information
!   2012-10-16 cucurull - increase the size of nreal and mreal by one element to
!                         add qrefges information, replace qcfail=5 by 4, add regional QC for MetOpB
!                         add dtype, dobs to distinguish use of toss_gps between ref/bending, add SR QC for obs
!   2014-01-28  todling - write sensitivity slot indicator (ioff) to header of diagfile
!   2014-12-13 derber   - minor optimization modifications
!   2015-07-28 cucurull - add QC for regional bending angle assimilation
!   2015-08-28  guo     - wrapped as a  module (m_gpsStats)
!                         moved the call to obsmod::destroy_genstats_gps() to
!                         where this routine was used (setuprhsall()), with its
!                         new module interface name. gpsStats_destroy().
!   2016-11-29 shlyaeva - increase the size of nreal for saving linearized Hx for EnKF
!   2016-12-09 mccarty  - add ncdiag writing support
!
!   input argument list:
!     toss_gps_sub  - array of qc'd profile heights
!     conv_diagsave - logical to save innovation diagnostics
!     mype          - mpi task id
!
!   output argument list:
!     bwork    - array containing information about obs-ges statistics
!     awork    - array containing information for data counts and gross checks
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind,r_single
  use obsmod, only: nprof_gps,lobsdiag_forenkf
  use obsmod, only: obs_diag,lobsdiagsave,luse_obsdiag
  use obsmod, only: binary_diag,netcdf_diag,dirname,ianldate
  use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, &
                          nc_diag_write, nc_diag_data2d
  use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_get_dim, nc_diag_read_close
  use gridmod, only: nsig,regional
  use constants, only: tiny_r_kind,half,wgtlim,one,two,zero,five,four
  use qcmod, only: npres_print,ptop,pbot
  use mpimod, only: ierror,mpi_comm_world,mpi_rtype,mpi_sum,mpi_max
  use jfunc, only: jiter,miter,jiterstart
  use gsi_4dvar, only: nobs_bins
  use convinfo, only: nconvtype
  use state_vectors, only: nsdim
  implicit none

! Declare passed variables
  logical                                          ,intent(in):: conv_diagsave
  integer(i_kind)                                  ,intent(in) :: mype
  real(r_kind),dimension(100+7*nsig)               ,intent(inout):: awork
  real(r_kind),dimension(npres_print,nconvtype,5,3),intent(inout):: bwork
  real(r_kind),dimension(max(1,nprof_gps))         ,intent(in):: toss_gps_sub

! Declare local parameters
  real(r_kind),parameter:: ten   = 10.0_r_kind
  real(r_kind),parameter:: six   = 6.0_r_kind
  real(r_kind),parameter:: r1em3 = 1.0e-3_r_kind
  real(r_kind),parameter:: r20   = 20.0_r_kind
  real(r_kind),parameter:: scale = 100.0_r_kind

! Declare local variables
  logical:: luse,muse,toss,save_jacobian
  integer(i_kind):: k,jsig,icnt,khgt,kprof,ikx,nn,j,nchar,nreal,mreal,ii,ioff
  real(r_kind):: pressure,arg,wgross,wgt,term,cg_gps,valqc,elev,satid,dtype,dobs
  real(r_kind):: ress,val,ratio_errors,val2
  real(r_kind):: exp_arg,data_ikx,data_rinc,cg_term,rat_err2,elat
  real(r_kind):: wnotgross,data_ipg,data_ier,data_ib,factor,super_gps_up,rhgt
  real(r_kind),dimension(nsig,max(1,nprof_gps)):: super_gps_sub,super_gps
  real(r_kind),dimension(max(1,nprof_gps)):: toss_gps
  real(r_kind),dimension(max(1,nprof_gps)):: high_gps,high_gps_sub
  real(r_kind),dimension(max(1,nprof_gps)):: dobs_height,dobs_height_sub

  real(r_single),allocatable,dimension(:,:)::sdiag
  character(8),allocatable,dimension(:):: cdiag

  real(r_single), dimension(nsdim) :: dhx_dx_array
  
  type(obs_diag), pointer :: obsptr => NULL()

  integer(i_kind) :: nnz, nind
  type(gps_ob_type), pointer:: gpsptr
  type(gps_all_ob_type), pointer:: gps_allptr
  
  save_jacobian = conv_diagsave .and. jiter==jiterstart .and. lobsdiag_forenkf

!*******************************************************************************
! Check to see if there are any profiles to process.  If none, return.
  if (nprof_gps==0) then
     if (mype==0) write(6,*)'GENSTATS_GPS:  no profiles to process (nprof_gfs=',nprof_gps,'), EXIT routine'
     return
  endif

! Reduce sub-domain specific QC'd profile height cutoff values to
! maximum global value for each profile
  toss_gps=zero
  call mpi_allreduce(toss_gps_sub,toss_gps,nprof_gps,mpi_rtype,mpi_max,&
       mpi_comm_world,ierror)

! If netcdf diag, initialize it
  if (conv_diagsave.and.netcdf_diag) call init_netcdf_diag_

! Get height of maximum bending angle
  dobs_height_sub = zero
  DO ii=1,nobs_bins
     gps_allptr => gps_allhead(ii)%head
     do while (associated(gps_allptr))

!       Load local work variables
        kprof        = gps_allptr%kprof
        dtype        = gps_allptr%rdiag(20)
        dobs         = gps_allptr%rdiag(17)

        if (dtype == one .and. toss_gps(kprof) > zero .and. dobs == toss_gps(kprof)) then
           dobs_height_sub(kprof) = gps_allptr%rdiag(7)
        endif

        gps_allptr => gps_allptr%llpoint

!    End loop over observations
     end do

! End of loop over time bins
  END DO

! Reduce sub-domain specific QC'd profile height to maximum global value for each profile
  dobs_height=zero
  call mpi_allreduce(dobs_height_sub,dobs_height,nprof_gps,mpi_rtype,mpi_max,&
       mpi_comm_world,ierror)


! Compute superobs factor on sub-domains using global QC'd profile height
  super_gps_sub=zero
  high_gps_sub = zero
  DO ii=1,nobs_bins
     gps_allptr => gps_allhead(ii)%head
     do while (associated(gps_allptr))

!       Load local work variables
        ratio_errors = gps_allptr%ratio_err
        data_ier     = gps_allptr%obserr
        luse         = gps_allptr%luse
        kprof        = gps_allptr%kprof
        dtype        = gps_allptr%rdiag(20)

!       Accumulate superobs factors and get highest good gps obs within a profile

        if (dtype == zero) then ! refractivity
          rhgt         = gps_allptr%loc
          if (rhgt >toss_gps(kprof)) then
             if(ratio_errors*data_ier>tiny_r_kind) then
                elev         = gps_allptr%rdiag(7)
                high_gps_sub(kprof)=max(high_gps_sub(kprof),elev)
                if(luse) then
                   khgt         = gps_allptr%loc
                   k=min(max(1,khgt),nsig)
                   super_gps_sub(k,kprof)=super_gps_sub(k,kprof)+one
                endif
             endif
          endif

        else ! bending angle
            dobs         = gps_allptr%rdiag(17)
            if(toss_gps(kprof) == zero .or. (toss_gps(kprof) > zero .and. dobs < toss_gps(kprof))) then ! will not fail SR from obs qc
              elev         = gps_allptr%rdiag(7)
              if(elev > dobs_height(kprof)) then
                if(ratio_errors*data_ier>tiny_r_kind) then
                   high_gps_sub(kprof)=max(high_gps_sub(kprof),elev)
                   if(luse) then
                      khgt         = gps_allptr%loc
                      k=min(max(1,khgt),nsig)
                      super_gps_sub(k,kprof)=super_gps_sub(k,kprof)+one
                   endif
                endif
              endif
            endif
        endif


        gps_allptr => gps_allptr%llpoint
        
!    End loop over observations
     end do

! End of loop over time bins     
  END DO

  super_gps = zero
  high_gps = zero
! Reduce sub-domain specifc superobs factors to global values for each profile
  call mpi_allreduce(super_gps_sub,super_gps,nsig*nprof_gps,mpi_rtype,mpi_sum,&
       mpi_comm_world,ierror)

! Reduce sub-domain specific high_gps values to global values for each profile
  call mpi_allreduce(high_gps_sub,high_gps,nprof_gps,mpi_rtype,mpi_max,&
       mpi_comm_world,ierror)

! Convert high_gps from meters to kilometers
  high_gps = r1em3*high_gps
  

! If generating diagnostic output, need to determine dimension of output arrays.
  nreal=0
  ioff =nreal
  if (conv_diagsave) then
     icnt = zero
     DO ii=1,nobs_bins
        gps_allptr => gps_allhead(ii)%head
        do while (associated(gps_allptr))
           luse         = gps_allptr%luse
           if(luse)icnt=icnt+1
           gps_allptr => gps_allptr%llpoint
        end do
     END DO
     if(icnt > 0)then
        nreal =22
        ioff  =nreal
        if (lobsdiagsave) nreal=nreal+4*miter+1
        if (save_jacobian) then
          nnz   = 3*nsig
          nind  = 3
          nreal = nreal + 2*nind + nnz + 2
        endif
        allocate(cdiag(icnt),sdiag(nreal,icnt))
     end if
  endif



! Loop over data to apply final qc, superobs factors, accumulate
! statistics and (optionally) load diagnostic output arrays
  icnt=0
  DO ii=1,nobs_bins
     gps_allptr => gps_allhead(ii)%head
     do while (associated(gps_allptr))

!       Load local work variables
        ratio_errors = gps_allptr%ratio_err
        data_ier     = gps_allptr%obserr
        luse         = gps_allptr%luse
        muse         = gps_allptr%muse
        khgt         = gps_allptr%loc
        kprof        = gps_allptr%kprof
        dtype        = gps_allptr%rdiag(20)
        gpsptr       => gps_allptr%mmpoint
        if(muse .and. associated(gpsptr) .and. luse_obsdiag)then
           obsptr       => gpsptr%diags
        endif

!       Determine model level to which observation is mapped to
        k=min(max(1,khgt),nsig)
 
!       Normalize ratio_errors by superobs factor.  Update ratio_error 
!       term used in minimization
        super_gps_up=zero
 
        if (super_gps(k,kprof)>tiny_r_kind) then
           do j=min(k+1,nsig),nsig
              super_gps_up = max(super_gps_up,super_gps(j,kprof))
           enddo

           if (super_gps_up >tiny_r_kind) then
              factor = one / sqrt(super_gps(k,kprof))
           else
              factor = one / sqrt(max(super_gps(k-1,kprof),super_gps(k,kprof)))
           endif
           ratio_errors = ratio_errors * factor
           if(conv_diagsave .and. luse) then
              if(gps_allptr%rdiag(16) >tiny_r_kind) gps_allptr%rdiag(16)=ratio_errors*data_ier
           endif
 
!          Adjust error ratio for observations used in inner loop
           if (associated(gpsptr)) then
              gpsptr%raterr2 = ratio_errors **2
              if(associated(obsptr) .and. luse_obsdiag)then
                 obsptr%wgtjo=(ratio_errors*data_ier)**2
              end if
           endif
        endif


!       For given profile, check if observation level is below level at 
!       which profile data is tossed.   If so, set error parameter to 
!       zero (effectively tossing the obs).
 
        rhgt = gps_allptr%loc
        mreal = 22
        if(dtype == zero) then !refractivity
          if (rhgt<=toss_gps(kprof)) then
             if(ratio_errors*data_ier > tiny_r_kind) then ! obs was good
                if (luse) then
                   if(conv_diagsave) then
                      gps_allptr%rdiag(10) = four
                      gps_allptr%rdiag(12) = -one
                      gps_allptr%rdiag(16) = zero
                      if(lobsdiagsave) gps_allptr%rdiag(mreal+jiter) = -one
                   endif
                   elat         = gps_allptr%rdiag(3)
                   if(elat > r20) then
                      awork(22) = awork(22)+one
                   else if(elat< -r20)then
                      awork(23) = awork(23)+one
                   else
                      awork(24) = awork(24)+one
                   end if
                endif
             endif
             ratio_errors = zero
             if (associated(gpsptr)) then
                gpsptr%raterr2 = ratio_errors **2
                if(associated(obsptr) .and. luse_obsdiag)then
                   obsptr%wgtjo=zero
                   obsptr%muse(jiter)=.false.
                end if
             endif
          endif
        else
          elev         = gps_allptr%rdiag(7)
          dobs         = gps_allptr%rdiag(17)
          if  (toss_gps(kprof) > zero .and. (dobs == toss_gps(kprof) .or. elev < dobs_height(kprof))) then ! SR from obs
              if(ratio_errors*data_ier > tiny_r_kind) then ! obs was good
                 if (luse) then
                    if(conv_diagsave) then
                      gps_allptr%rdiag(10) = four
                      gps_allptr%rdiag(12) = -one
                      gps_allptr%rdiag(16) = zero
                      if(lobsdiagsave) gps_allptr%rdiag(mreal+jiter) = -one
                    endif
                    elat         = gps_allptr%rdiag(3)
                    if(elat > r20) then
                      awork(22) = awork(22)+one
                    else if(elat< -r20)then
                      awork(23) = awork(23)+one
                    else
                      awork(24) = awork(24)+one
                    end if
                 endif
              endif
              ratio_errors = zero
              if (associated(gpsptr)) then
                 gpsptr%raterr2 = ratio_errors **2
                 if(associated(obsptr) .and. luse_obsdiag)then
                    obsptr%wgtjo=zero
                    obsptr%muse(jiter)=.false.
                 end if
              endif
          endif
        endif



!       Regional QC.  Remove obs if highest good obs in
!       profile is below platform specific threshold height.
        if(regional) then
           toss=.false.
           if(ratio_errors*data_ier > tiny_r_kind) then
             if(dtype==zero) then !refractivity
                satid        = gps_allptr%rdiag(1)
                if((satid==41).or.(satid==722).or.(satid==723).or.(satid==4).or.(satid==786).or.(satid==3)) then
                   if ((high_gps(kprof)) < ten)  toss=.true.
                else ! OL
                   if ((high_gps(kprof)) < five) toss=.true.
                endif
             else !bending angle
                if ((high_gps(kprof)) <= six)  toss=.true.
             endif
           endif
           if (toss) then
              if (luse) then
                 if(conv_diagsave) then
                    gps_allptr%rdiag(10) = four
                    gps_allptr%rdiag(12) = -one
                    gps_allptr%rdiag(16) = zero
                    if(lobsdiagsave) gps_allptr%rdiag(mreal+jiter) = -one
                 endif
                 elat         = gps_allptr%rdiag(3)
                 if(elat > r20) then
                    awork(22) = awork(22)+one
                 else if(elat< -r20)then
                    awork(23) = awork(23)+one
                 else
                    awork(24) = awork(24)+one
                 end if
              end if
              ratio_errors = zero
              if (associated(gpsptr)) then
                 gpsptr%raterr2 = ratio_errors **2
                 if(associated(obsptr) .and. luse_obsdiag)then
                    obsptr%wgtjo=zero
                    obsptr%muse(jiter)=.false.
                 end if
              endif
           endif
        endif ! regional

!       Compute penalty terms
        if (ratio_errors*data_ier <= tiny_r_kind) muse = .false.
        if(luse)then
           val          = gps_allptr%dataerr
           val2     = val*val
           exp_arg  = -half*val2
           rat_err2 = ratio_errors**2
           data_ipg     = gps_allptr%pg
           if (data_ipg > tiny_r_kind) then
              data_ib      = gps_allptr%b
              cg_gps=cg_term/data_ib
              wnotgross= one-data_ipg
              wgross   = data_ipg*cg_gps
              arg      = exp(exp_arg)
              term     = log(wnotgross*arg+wgross)
              wgt      = wnotgross*arg/(wnotgross*arg+wgross)
           else
              term = exp_arg
              wgt  = one
           endif
           if(conv_diagsave) gps_allptr%rdiag(13) = wgt/wgtlim
           valqc = -two*rat_err2*term
         

!          Accumulate statistics for obs belonging to this task
!          based on interface (not mid-point) level
           val2=val2*rat_err2
           if(muse)then
              if(wgt < wgtlim) awork(21) = awork(21)+one

!             Accumulate values for penalty and data count
              jsig=max(1,khgt)
              awork(jsig+3*nsig+100)=awork(jsig+3*nsig+100)+valqc
              awork(jsig+5*nsig+100)=awork(jsig+5*nsig+100)+one
              awork(jsig+6*nsig+100)=awork(jsig+6*nsig+100)+val2
              nn=1
           else
              nn=2
              if(ratio_errors*data_ier >=tiny_r_kind)nn=3
           endif

           data_ikx     = gps_allptr%type
           ikx          = nint(data_ikx)
           pressure     = gps_allptr%rdiag(6)
           data_rinc    = gps_allptr%rdiag(5)*scale
!          Loop over pressure level groupings and obs to accumulate
!          statistics as a function of observation type.
           do k = 1,npres_print
              if(pressure>ptop(k) .and. pressure<=pbot(k))then
                 ress=data_rinc
                 
                 bwork(k,ikx,1,nn)  = bwork(k,ikx,1,nn)+one           ! count
                 bwork(k,ikx,2,nn)  = bwork(k,ikx,2,nn)+ress          ! (o-g)
                 bwork(k,ikx,3,nn)  = bwork(k,ikx,3,nn)+ress*ress     ! (o-g)**2
                 bwork(k,ikx,4,nn)  = bwork(k,ikx,4,nn)+val2          ! penalty
                 bwork(k,ikx,5,nn)  = bwork(k,ikx,5,nn)+valqc         ! nonlin qc penalty
                 
              end if
           end do
        end if

!       Transfer diagnostic information to output arrays
        if(conv_diagsave .and. luse) then
           icnt=icnt+1
           cdiag(icnt) = gps_allptr%cdiag
           do j=1,nreal
              sdiag(j,icnt)= gps_allptr%rdiag(j)
           enddo
        endif


        if (conv_diagsave .and. netcdf_diag .and. luse) call contents_netcdf_diag_
        
        gps_allptr => gps_allptr%llpoint

!    End loop over observations
     end do

! End of loop over time bins
  END DO

! If requested, write information to diagnostic file
  if(conv_diagsave) then
    if (netcdf_diag) call nc_diag_write
    if (binary_diag .and. icnt > 0)then
        nchar = 1
        write(7)'gps',nchar,nreal,icnt,mype,ioff
        write(7)cdiag,sdiag
        deallocate(cdiag,sdiag)
    endif
  endif


! Destroy arrays holding gps data
  call destroy_genstats_gps
contains

subroutine init_netcdf_diag_
  integer(i_kind) ncd_fileid, ncd_nobs
  character(len=80) string
  character(len=128) diag_conv_file
  logical append_diag  
  logical,parameter::verbose=.false.

     write(string,900) jiter
900  format('conv_gps_',i2.2,'.nc4')
     diag_conv_file=trim(dirname) // trim(string)

     inquire(file=diag_conv_file, exist=append_diag)

     if (append_diag) then
        call nc_diag_read_init(diag_conv_file,ncd_fileid)
        ncd_nobs = nc_diag_read_get_dim(ncd_fileid,'nobs')
        call nc_diag_read_close(diag_conv_file)

        if (ncd_nobs > 0) then
           if(verbose) print *,'file ' // trim(diag_conv_file) // ' exists.  Appending.  nobs,mype=',ncd_nobs,mype
        else
           if(verbose) print *,'file ' // trim(diag_conv_file) // ' exists but contains no obs.  Not appending. nobs,mype=',ncd_nobs,mype
           append_diag = .false. ! if there are no obs in existing file, then do not try to append
        endif
     end if

     call nc_diag_init(diag_conv_file, append=append_diag)

     if (.not. append_diag) then ! don't write headers on append - the module will break?
        call nc_diag_header("date_time",ianldate )
        call nc_diag_header("Number_of_state_vars", nsdim          )
     endif
end subroutine init_netcdf_diag_

subroutine contents_binary_diag_
end subroutine contents_binary_diag_

subroutine contents_netcdf_diag_
  use sparsearr, only: sparr2, readarray, fullarray
  integer(i_kind),dimension(miter) :: obsdiag_iuse
  integer(i_kind)                  :: obstype, obssubtype
  type(sparr2) :: dhx_dx

! Observation class
  character(7),parameter     :: obsclass = '    gps'

           call nc_diag_metadata("Station_ID",                            gps_allptr%cdiag             )
           call nc_diag_metadata("Observation_Class",                     obsclass                     )
           obstype    = gps_allptr%rdiag(1) 
           obssubtype = gps_allptr%rdiag(2)
           call nc_diag_metadata("Observation_Type",                      obstype                      )
           call nc_diag_metadata("Observation_Subtype",                   obssubtype                   )
           call nc_diag_metadata("Latitude",                              sngl(gps_allptr%rdiag(3))    )
           call nc_diag_metadata("Longitude",                             sngl(gps_allptr%rdiag(4))    )
           call nc_diag_metadata("Incremental_Bending_Angle",             sngl(gps_allptr%rdiag(5))    )
           call nc_diag_metadata("Pressure",                              sngl(gps_allptr%rdiag(6))    )
           call nc_diag_metadata("Height",                                sngl(gps_allptr%rdiag(7))    )
           call nc_diag_metadata("Time",                                  sngl(gps_allptr%rdiag(8))    )
           call nc_diag_metadata("Model_Elevation",                       sngl(gps_allptr%rdiag(9))    )
           call nc_diag_metadata("Setup_QC_Mark",                         sngl(gps_allptr%rdiag(10))   )
           call nc_diag_metadata("Prep_Use_Flag",                         sngl(gps_allptr%rdiag(11))   )
           call nc_diag_metadata("Analysis_Use_Flag",                     sngl(gps_allptr%rdiag(12))   )

           call nc_diag_metadata("Nonlinear_QC_Rel_Wgt",                  sngl(gps_allptr%rdiag(13))   )
           call nc_diag_metadata("Errinv_Input",                          sngl(gps_allptr%rdiag(14))   )
           call nc_diag_metadata("Errinv_Adjust",                         sngl(gps_allptr%rdiag(15))   )
           call nc_diag_metadata("Errinv_Final",                          sngl(gps_allptr%rdiag(16))   )
           call nc_diag_metadata("Observation",                           sngl(gps_allptr%rdiag(17))   )
           call nc_diag_metadata("Obs_Minus_Forecast_adjusted",           sngl(gps_allptr%rdiag(17))*sngl(gps_allptr%rdiag(5)) )
           call nc_diag_metadata("Obs_Minus_Forecast_unadjusted",         sngl(gps_allptr%rdiag(17))*sngl(gps_allptr%rdiag(5)) )
           call nc_diag_metadata("GPS_Type",                              sngl(gps_allptr%rdiag(20))   )
           call nc_diag_metadata("Temperature_at_Obs_Location",           sngl(gps_allptr%rdiag(18))   )
           call nc_diag_metadata("Specific_Humidity_at_Obs_Location",     sngl(gps_allptr%rdiag(21))   )

           if (save_jacobian) then
              call readarray(dhx_dx, gps_allptr%rdiag(ioff+1:nreal))
              call fullarray(dhx_dx, dhx_dx_array)
              call nc_diag_data2d("Observation_Operator_Jacobian", dhx_dx_array)
           endif



!           call nc_diag_data2d("T_Jacobian",                              gps_allptr%mmpoint%jac_t             )
           if (lobsdiagsave) then
              print *,'ERROR: OBSDIAGSAVE SKIPPED IN NCDIAG DEVELOPMENT.  STOPPING.'
              call stop2(55)
!              do jj=1,miter
!                 if (gps_allptr%diags%muse(jj)) then
!                    obsdiag_iuse(jj) =  one
!                 else
!                    obsdiag_iuse(jj) = -one
!                 endif
!              enddo
!
!              call nc_diag_data2d("ObsDiagSave_iuse",     obsdiag_iuse                             )
!              call nc_diag_data2d("ObsDiagSave_nldepart", gps_allptr%diags%nldepart )
!              call nc_diag_data2d("ObsDiagSave_tldepart", gps_allptr%diags%tldepart )
!              call nc_diag_data2d("ObsDiagSave_obssen",   gps_allptr%diags%obssen   )
           endif
end subroutine contents_netcdf_diag_
end subroutine genstats_gps

end module m_gpsStats