!*******************************************************************
      subroutine dpthresh(ns,nm,nl,nlat,nlon,h,s,cl0,cl1,dwin,thr,p,mval)

!  DESCRIPTION:
!     Computes an array of climatological quantile threshold values 
!     thr(nl,lon,lat) for daily DP data, corresponding to threshold value p.        
!     Climatology is defined over years cl0-cl1. A moving window of length
!     dwin is used to augment the sample size.
!     
!  References:
!    *Eade et al., 2012, JGR, doi:10.1029/2012JD018015
!    *Hamilton et al., 2012, JGR, doi:10.1029/2011JD016541
!
! argument :
      implicit none
      integer ::  ns,nm,nl,nlat,nlon
      integer ::  dwin          ! window size (days)
      integer :: s(ns)             ! hindcast start years
      integer :: cl0,cl1           ! climatology start/end years
      real    :: &
          p,               &    ! quantile level (e.g. 0.95)
          mval                  ! _FillValue
      real, dimension(nlon,nlat,nm,nl,ns) :: &
          h                     ! DP hindcast data
      real, dimension(nlon,nlat,nl) :: &
          thr                   ! threshold values for quantile p
! non-argument :
      integer ::  nsamp
      integer ::  i,j,k,l,iyearadd
      integer :: year(ns)          ! calendar year of data
      integer :: k0,k1             ! lead day window indices
      logical, dimension(nm,nl,ns) :: &
          hmask                 ! mask for selecting DP hindcast data
      real, dimension(:), allocatable :: sample


! 
  if (mod(dwin,2)==0) then
    dwin=dwin+1
  endif
! write(6,*) 'dwin = ',dwin
! write(6,*) 'cl0,cl1 = ',cl0,cl1

! Loop over lon,lat
      do j=1,nlat
       do i=1,nlon

        if (all(h(i,j,1,:,:)==mval)) then
           thr(i,j,:) = mval
        else

! Loop over lead times (daily)
         do k=1,nl
           hmask = .false.

!          write(6,*) 'k = ',k
           ! Define lead-day window
           k0 = max(1,k-floor(dwin/2.))
           k1 = min(nl,k+floor(dwin/2.))
!          write(6,*) '  k0,k1 = ',k0,k1

           ! Determine calendar year
           iyearadd = floor((k-62)/365.0)+1
           year = s+iyearadd
!          write(6,*) '  year = ',year
           
           ! Define lead days and start years that
           ! will contribute to sample:
           do l=1,ns
             if (year(l).ge.cl0.and.year(l).le.cl1) then
              where(h(i,j,:,k0:k1,l).ne.mval) hmask(:,k0:k1,l) = .true.
!             hmask(:,k0:k1,l) = .true.
             end if
           end do

           nsamp = count(hmask)
           allocate(sample(nsamp))
           sample = pack(h(i,j,:,:,:),hmask)
           write(6,*) '  nsamp = ',nsamp

           ! Compute thr
           thr(i,j,k) =  pquant(nsamp,sample,p)

           deallocate(sample)

         end do
        end if

       end do
      end do

   contains

!*******************************************************************
   function pquant(nx,x,y)
! DESCRIPTION:
! Given a 1D array x(nx), determine value of x that is the y-
! quantile.

   implicit none
   real :: pquant       ! Quantile value
   real :: y            ! 0 < y < 1
   real, dimension(nx), intent(inout) :: x         ! unsorted sample
   integer :: nx,yindex

   call QsortC(x)                 ! quicksort routine
   yindex = max(1,int(y*nx))      ! index of y-quantile
   yindex = min(yindex,nx)
   pquant = x(yindex)

   end function pquant
!*******************************************************************

!*******************************************************************
   recursive subroutine QsortC(A)
! DESCRIPTION:
! Recursive Fortran 95 quicksort routine
! sorts real numbers into ascending numerical order
! Author: Juli Rew, SCD Consulting (juliana@ucar.edu), 9/03
! Based on algorithm from Cormen et al., Introduction to Algorithms,
! 1997 printing
! Copied from:  http://www.fortran.com/qsort_c.f95
   real, intent(in out), dimension(:) :: A
   integer :: iq

   if(size(A) > 1) then
     call Partition(A, iq)
     call QsortC(A(:iq-1))
     call QsortC(A(iq:))
   endif
   end subroutine QsortC
!*******************************************************************

!*******************************************************************
   subroutine Partition(A, marker)
! DESCRIPTION:
! Recursive Fortran 95 quicksort routine
! sorts real numbers into ascending numerical order
! Author: Juli Rew, SCD Consulting (juliana@ucar.edu), 9/03
! Based on algorithm from Cormen et al., Introduction to Algorithms,
! 1997 printing
! Copied from:  http://www.fortran.com/qsort_c.f95
   real, intent(in out), dimension(:) :: A
   integer, intent(out) :: marker
   integer :: i, j
   real :: temp
   real :: x      ! pivot point
   x = A(1)
   i= 0
   j= size(A) + 1
 
   do
     j = j-1
     do
        if (A(j) <= x) exit
        j = j-1
     end do
     i = i+1
     do
        if (A(i) >= x) exit
        i = i+1
     end do
     if (i < j) then
        ! exchange A(i) and A(j)
        temp = A(i)
        A(i) = A(j)
        A(j) = temp
     elseif (i == j) then
        marker = i+1
        return
     else
        marker = i
        return
     endif
   end do

   end subroutine Partition
!*******************************************************************

 end subroutine dpthresh
!*******************************************************************