module satthin 16,1
!$$$ subprogram documenation block
!                .      .    .                                       .
! subprogram:    satthin
!   prgmmr: treadon          org: np20                date: 2002-10-17
!
! abstract:  This module contains code which may be used to selectively
!            thin satellite data.
!
! program history log:
!   2002-10-17 treadon
!   2004-05-28  kleist, subroutine call update
!   2005-03-28  wu, remove unused variable mlath
!   2005-10-06  treadon - add routine destroy_sfc
!   2005-11-29  parrish - add routines getsfc_nmm_binary, getsfc_nmm_netcdf,
!                         getsfc_mass_binary, getsfc_mass_netcdf
!   2005-12-08  treadon - rename global getsfc to getsfc_global, add getsfc
!                         to branch to appropriate getsfc* routine
!   2006-02-01  parrish - remove all getsfc routines, and create new getsfc
!                         which makes full sfc fields from guess_grids.
!                         rename surface fields to sst_full,sno_full, etc.
!                         to eliminate conflict with some fields with same
!                         names in guess_grids.
!                         also add new sfc fields u10_full,v10_full
!   2006-05-23  parrish - add new sfc field zs_full, the model terrain height
!   2006-07-28  derber  - use r1000 from constants
!
! Subroutines Included:
!   sub makegvals      - set up for superob weighting
!   sub makegrids      - set up thinning grids
!   sub getsfc         - create full horizontal fields of surface arrays
!   sub map2tgrid      - map observation to location on thinning grid
!   sub destroygrids   - deallocate thinning grid arrays
!   sub destroy_sfc    - deallocate full horizontal surface arrays
!   sub indexx         - sort array into ascending order
!
! Variable Definitions:
!   def mlat           - number of latitudes in thinning grid
!   def mlon           - number of longitudes in thinning grid
!   def superp         - maximum number of data types
!   def maxthin        - maximum number of obs to retain in thinning grid box
!   def itxmax         - total number of points in thinning grid
!   def istart_val     - starting location on thinning grid for superobs
!   def icount         - observation counter
!   def ibest_obs      - index of "best" quality obs in thinning grid box
!   def isli_full      - snow/land/ice mask
!   def glat           - latitudes on thinning grid
!   def super_val      - super obs factor across data types
!   def super_val1     - super obs factors summed over all mpi tasks (data types)
!   def glon           - longitudes on thinning grid
!   def hll            - (i,j) index of point on thinning grid
!   def sli_full       - 0=sea/1=land/2=ice mask
!   def sst_full       - skin temperature
!   def sno_full       - snow-ice mask
!   def u10_full       - 10 meter zonal wind
!   def v10_full       - 10 meter meridional wind
!   def zs_full        - model terrain elevation
!   def score_crit     - "best" quality obs score in thinning grid box
!   def use_all        - parameter for turning satellite thinning algorithm off
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  implicit none

  integer(i_kind) mlat,superp,maxthin,itxmax
  integer(i_kind),dimension(0:51):: istart_val
  
  integer(i_kind),allocatable,dimension(:):: mlon,icount,ibest_obs
  integer(i_kind),allocatable,dimension(:,:):: isli_full

  real(r_kind) rlat_min,rlat_max,rlon_min,rlon_max,dlat_grid,dlon_grid
  real(r_kind),allocatable,dimension(:):: glat
  real(r_kind),allocatable,dimension(:):: super_val,super_val1
  real(r_kind),allocatable,dimension(:,:):: glon,hll,sli_full,sst_full,sno_full
  real(r_kind),allocatable,dimension(:,:):: u10_full,v10_full,zs_full
  real(r_kind),allocatable,dimension(:):: score_crit
  logical use_all

contains


  subroutine makegvals 1,4
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    makegvals                            
!     prgmmr:    derber      org: np23                date: 2002-10-17
!
! abstract:  This routine allocates and initializes arrays
!            used for superobs weighting.
!
! program history log:
!   2002-10-17  derber
!   2004-06-22  treadon - update documentation
!   2004-12-09  treadon - allocate thinning grids consistent with analysis domain
!   2006-07-28  derber  - use r1000 from constants
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$
    use constants, only: deg2rad,rearth_equator,zero,two,pi,half,one,quarter,&
       rad2deg,r1000
    use obsmod, only: dmesh,dthin,ndat
    use gridmod, only: regional,region_lat,region_lon,nlat,nlon,txy2ll
    implicit none

	real(r_kind),parameter:: r90   = 90.0_r_kind
    real(r_kind),parameter:: r360  = 360.0_r_kind
    real(r_kind),parameter:: r999  = 999.0_r_kind

    logical odd
    integer(i_kind) i,ii,j,k
    integer(i_kind) mlonx,icnt,mlony,mlonj
    real(r_kind) delonx,delat,dgv,dx,dy
    real(r_kind) twopi,dlon_g,dlat_g,dlon_e,dlat_e
    real(r_kind) factor,factors,delon
    real(r_kind) rkm2dg,glatm,glatx

!   Initialize variables, set constants
    maxthin=0
    do i=1,ndat
       maxthin=max(maxthin,abs(dthin(i)))
    end do
    istart_val=0
    twopi  = two*pi
    rkm2dg = r360/(twopi*rearth_equator)*r1000

!   Set lat,lon limits for analysis grid.
    rlat_min = -r90
    rlat_max = r90
    rlon_min = zero
    rlon_max = r360
    if (regional) then
       rlat_min =  r999
       rlat_max = -r999
       rlon_min =  r999
       rlon_max = -r999
       do j=1,nlon
          dlon_g=j
          do i=1,nlat
             dlat_g=i
             call txy2ll(dlon_g,dlat_g,dlon_e,dlat_e)
             dlat_e=dlat_e*rad2deg
             dlon_e=dlon_e*rad2deg
             rlat_min = min(rlat_min,dlat_e)
             rlat_max = max(rlat_max,dlat_e)
             rlon_min = min(rlon_min,dlon_e)
             rlon_max = max(rlon_max,dlon_e)
          end do
       end do
       if (rlon_min < zero) rlon_min = rlon_min + r360
       if (rlon_max < zero) rlon_max = rlon_max + r360
    endif
    dlat_grid = rlat_max - rlat_min
    dlon_grid = rlon_max - rlon_min

    do ii=1,maxthin

!      Set up dimensions for thinning grids
       istart_val(ii+1)=istart_val(ii)
       if(abs(dmesh(ii)) > one)then
          dx    = dmesh(ii)*rkm2dg
          dy    = dx
          mlat  = dlat_grid/dy + half
          mlonx = dlon_grid/dx + half
          delat = dlat_grid/mlat
          delonx= dlon_grid/mlonx	
          dgv   = delat*half
          mlat=max(2,mlat);   mlonx=max(2,mlonx)

          icnt=0
          do j = 1,mlat
             glatx = rlat_min + (j-1)*delat
             glatx = glatx*deg2rad
             glatm = glatx + dgv*deg2rad
             
             factor = abs(cos(abs(glatm)))
             if (dmesh(ii)>zero) then
                mlonj = nint(mlonx*factor)
                mlony = max(2,mlonj)
             else
                delon = factor*dmesh(ii)
                delon = min(delon,r360)
                mlony = dlon_grid/delon
             endif
             do i = 1,mlony
                icnt=icnt+1
                istart_val(ii+1)=istart_val(ii+1)+1
             enddo

          enddo
       end if
    end do
    superp=istart_val(maxthin+1)
    
!   Allocate and initialize arrays for superobs weighthing
    allocate(super_val(0:superp),super_val1(0:superp))
    do i=0,superp
       super_val(i)=zero
    end do
    
    return
  end subroutine makegvals



  subroutine makegrids(rmesh) 9,1
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    makegrids                            
!     prgmmr:    treadon     org: np23                date: 2002-10-17
!
! abstract:  This routine sets up dimensions for and allocates
!            thinning grids.
!
! program history log:
!   2002-10-17  treadon
!   2004-06-22  treadon - update documentation
!   2004-12-09  treadon - allocate thinning grids consistent with analysis domain
!
!   input argument list:
!     rmesh - mesh size (km) of thinning grid.  If (rmesh <= one), 
!             then no thinning of the data will occur.  Instead,
!             all data will be used without thinning.
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$
    use constants, only: rearth_equator,two,deg2rad,zero,half,one,quarter,pi

    implicit none

    real(r_kind),parameter:: r360 = 360.0_r_kind
    logical odd
    integer(i_kind) i,j,k
    integer(i_kind) mlonx,mlonj
    real(r_kind) delonx,delat,dgv,halfpi,dx,dy
    real(r_kind) twopi
    real(r_kind) factor,factors,delon,rmesh
    real(r_kind) rkm2dg,glatm


!   If there is to be no thinning, simply return to calling routine
    use_all=.false.
    if(abs(rmesh) <= one)then
      use_all=.true.
      itxmax=2.e6
      return
    end if

!   Set constants
    halfpi = half*pi
    twopi  = two*pi
    rkm2dg = r360/(twopi*rearth_equator)*1.e3_r_kind

!   Set up dimensions and allocate arrays for thinning grids
    if (rmesh<zero) rkm2dg=one
    dx    = rmesh*rkm2dg
    dy    = dx
    mlat  = dlat_grid/dy + half
    mlonx = dlon_grid/dx + half
	delat = dlat_grid/mlat
	delonx= dlon_grid/mlonx
    dgv  = delat*half
    mlat=max(2,mlat);   mlonx=max(2,mlonx)

    allocate(mlon(mlat),glat(mlat),glon(mlonx,mlat),hll(mlonx,mlat))


!   Set up thinning grid lon & lat.  The lon & lat represent the location of the
!   lower left corner of the thinning grid box.
    itxmax=0
    do j = 1,mlat
       glat(j) = rlat_min + (j-1)*delat
       glat(j) = glat(j)*deg2rad
       glatm = glat(j) + dgv*deg2rad

       factor = abs(cos(abs(glatm)))
       if (rmesh>zero) then
          mlonj   = nint(mlonx*factor)	
          mlon(j) = max(2,mlonj)
          delon = dlon_grid/mlon(j)
       else
          delon = factor*rmesh
          delon = min(delon,r360)
          mlon(j) = dlon_grid/delon
       endif

       glat(j) = min(max(-halfpi,glat(j)),halfpi)
       do i = 1,mlon(j)
          itxmax=itxmax+1
          hll(i,j)=itxmax
          glon(i,j) = rlon_min + (i-1)*delon
          glon(i,j) = glon(i,j)*deg2rad
          glon(i,j) = min(max(zero,glon(i,j)),twopi)
       enddo

    end do


!   Allocate  and initialize arrays
    allocate(icount(itxmax))
    allocate(ibest_obs(itxmax))
    allocate(score_crit(itxmax))

    do j=1,itxmax
          icount(j)=0
    end do
    do j=1,itxmax
       score_crit(j) =9.99e6_r_kind
    end do

    return
  end subroutine makegrids


  subroutine getsfc(mype) 1,12
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    getsfc
!     prgmmr:    parrish     org: np23                date: 2006-02-02
!
! abstract:  This routine converts subdomain surface fields in
!            guess_grids to full horizontal fields for use in 
!            reading of observations.
!
! program history log:
!   2002-10-17  parrish
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$
    use kinds, only: r_kind,i_kind
    use gridmod, only:  nlat,nlon,lat2,lon2,lat1,lon1,&
     ltosi,ltosj,iglobal,itotsub,ijn,displs_g
    use guess_grids, only: ntguessig,isli,sfct,sno,fact10,ges_u,ges_v,ges_z
    use mpimod, only: mpi_comm_world,ierror,mpi_rtype,strip
    use constants, only: zero

    implicit none

    integer(i_kind) mype

! Local variables
    real(r_kind),dimension(lat1*lon1):: zsm
    real(r_kind),dimension(itotsub):: work1
    real(r_kind),dimension(lat2,lon2):: work2
    real(r_kind),dimension(nlat,nlon):: fact10_full
    integer(i_kind) mm1,i,j,k,it

    mm1=mype+1
    it=ntguessig

    allocate(isli_full(nlat,nlon),sli_full(nlat,nlon))
    allocate(sst_full(nlat,nlon),sno_full(nlat,nlon))
    allocate(u10_full(nlat,nlon),v10_full(nlat,nlon))
	allocate(zs_full(nlat,nlon))

! isli_full,sli_full

    do j=1,lon1*lat1
      zsm(j)=zero
    end do

    do j=1,lon2
      do i=1,lat2
        work2(i,j)=isli(i,j,it)
      end do
    end do
    call strip(work2,zsm,1)

    call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,&
       work1,ijn,displs_g,mpi_rtype,&
       mpi_comm_world,ierror)

    do k=1,iglobal
      i=ltosi(k) ; j=ltosj(k)
      sli_full(i,j)=work1(k)
      isli_full(i,j)=nint(work1(k))
    end do

! sst_full

    do j=1,lon1*lat1
      zsm(j)=zero
    end do

    call strip(sfct(1,1,it),zsm,1)

    call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,&
       work1,ijn,displs_g,mpi_rtype,&
       mpi_comm_world,ierror)

    do k=1,iglobal
      i=ltosi(k) ; j=ltosj(k)
      sst_full(i,j)=work1(k)
    end do

! sno_full

    do j=1,lon1*lat1
      zsm(j)=zero
    end do

    call strip(sno(1,1,it),zsm,1)

    call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,&
       work1,ijn,displs_g,mpi_rtype,&
       mpi_comm_world,ierror)

    do k=1,iglobal
      i=ltosi(k) ; j=ltosj(k)
      sno_full(i,j)=work1(k)
    end do

! fact10_full

    do j=1,lon1*lat1
      zsm(j)=zero
    end do

    call strip(fact10(1,1,it),zsm,1)

    call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,&
       work1,ijn,displs_g,mpi_rtype,&
       mpi_comm_world,ierror)

    do k=1,iglobal
      i=ltosi(k) ; j=ltosj(k)
      fact10_full(i,j)=work1(k)
    end do

! u10_full

    do j=1,lon1*lat1
      zsm(j)=zero
    end do

    call strip(ges_u(1,1,1,it),zsm,1)

    call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,&
       work1,ijn,displs_g,mpi_rtype,&
       mpi_comm_world,ierror)

    do k=1,iglobal
      i=ltosi(k) ; j=ltosj(k)
      u10_full(i,j)=work1(k)
    end do

! v10_full

    do j=1,lon1*lat1
      zsm(j)=zero
    end do

    call strip(ges_v(1,1,1,it),zsm,1)

    call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,&
       work1,ijn,displs_g,mpi_rtype,&
       mpi_comm_world,ierror)

    do k=1,iglobal
      i=ltosi(k) ; j=ltosj(k)
      v10_full(i,j)=work1(k)
    end do

! zs_full

    do j=1,lon1*lat1
      zsm(j)=zero
    end do

    call strip(ges_z(1,1,it),zsm,1)

    call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,&
       work1,ijn,displs_g,mpi_rtype,&
       mpi_comm_world,ierror)

    do k=1,iglobal
      i=ltosi(k) ; j=ltosj(k)
      zs_full(i,j)=work1(k)
    end do

  end subroutine getsfc


  subroutine map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse) 9,4
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    map2tgrid
!     prgmmr:    derber      org: np2                 date: 2006-05-03
!
! abstract:  This routine maps observations to the thinning grid.
!
! program history log:
!   2006-05-03  derber (created from map2grids)
!   2006-09-13  treadon - set itx=1 for the case use_all=.true.
!
!   input argument list:
!     dlat_earth - earth relative observation latitude (radians)
!     dlon_earth - earth relative observation longitude (radians)
!     crit1      - quality indicator for observation (smaller = better)
!     ithin      - number of obs to retain per thinning grid box
!
!   output argument list:
!     iobs  - observation counter
!     itx   - combined (i,j) index of observation on thinning grid
!     itt   - superobs thinning counter
!     iobsout- location for observation to be put
!     iuse  - .true. if observation should be used
!     
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$
    use kinds, only: r_kind,i_kind
    use constants, only: one, half
    implicit none
    logical,intent(out):: iuse
    integer(i_kind),intent(in):: ithin
    integer(i_kind),intent(out):: itt,itx
    real(r_kind),intent(in):: dlat_earth,dlon_earth,crit1
    real(r_kind),intent(out):: dist1

    integer(i_kind) ix,iy
    real(r_kind) dlat1,dlon1,dx,dy,dxx,dyy


!   If using all data (no thinning), simply return to calling routine
    if(use_all)then
       iuse=.true.
       itt=0
       itx=1
       return
    end if

!   Compute (i,j) indices of coarse mesh grid (grid number 1) which 
!   contains the current observation.
    dlat1=dlat_earth
    dlon1=dlon_earth

    call grdcrd(dlat1,1,glat,mlat,1)
    iy=int(dlat1)
    dy=dlat1-iy
    iy=max(1,min(iy,mlat))

    call grdcrd(dlon1,1,glon(1,iy),mlon(iy),1)
    ix=int(dlon1)
    dx=dlon1-ix
    ix=max(1,min(ix,mlon(iy)))

    dxx=half-min(dx,one-dx)
    dyy=half-min(dy,one-dy)
    dist1=dxx*dxx+dyy*dyy+half
    itx=hll(ix,iy)
    itt=istart_val(ithin)+itx
    if(ithin == 0) itt=0

!   Increment obs counter on coarse mesh grid.  Also accumulate observation
!   score and distance functions

!   ratio=1.e9
!   if ( dx > zero ) ratio=dy/dx
!   dista=sin(two*atan(ratio))
!   distb=sin(pi*dx)                !dista+distb is max at grid box center
!   dist1=one - quarter*(dista + distb)  !dist1 is min at grid box center and 
                                    !ranges from 1 (at corners)to 
                                    !.5 (at center of box)
    iuse=.true.

    if(dist1*crit1 > score_crit(itx) .and. icount(itx) == 0)iuse=.false.

    return
  end subroutine map2tgrid

  subroutine checkob(dist1,crit1,itx,iuse) 7,1
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    checkob
!     prgmmr:    checkob     org: np23                date: 2002-10-17
!
! abstract:  intermediate ob checking routine to see if we should not use.
!
! program history log:
!   2006-05-03  derber
!
!   input argument list:
!     dist1  - quality indicator for distance (smaller = better)
!     crit1      - quality indicator for observation (smaller = better)
!     itx   - combined (i,j) index of observation on thinning grid
!     iuse  - .true. if observation should be used
!
!   output argument list:
!     iuse  - .true. if observation should be used
!     
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$
    use kinds, only: r_kind,i_kind
    implicit none
    logical,intent(inout):: iuse
    integer(i_kind),intent(in):: itx
    real(r_kind),intent(in):: dist1,crit1

!   If using all data (no thinning), simply return to calling routine
    if(use_all .or. .not. iuse .or. icount(itx)==0)return

    if(crit1*dist1 > score_crit(itx))iuse=.false.

    return
  end subroutine checkob

  subroutine finalcheck(dist1,crit1,iobs,itx,iobsout,iuse,sis) 9,1
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    finalcheck
!     prgmmr:    derber     org: np23                date: 2002-10-17
!
! abstract:  This routine performs the final criterion check for sat 
!              obs and increments counters
!
! program history log:
!   2006-05-03  derber
!
!   input argument list:
!     dist1  - quality indicator for distance (smaller = better)
!     crit1  - quality indicator for observation (smaller = better)
!     itx    - combined (i,j) index of observation on thinning grid
!     iobs   - observation counter
!     iuse   - .true. if observation should be used
!     sis    - sensor/instrument/satellite 
!
!   output argument list:
!     iobs   - observation counter
!     iobsout- location for observation to be put
!     iuse   - .true. if observation should be used
!     
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$
    use kinds, only: r_kind,i_kind
    implicit none
    logical,intent(inout):: iuse
    integer(i_kind),intent(inout):: iobs,iobsout
    integer(i_kind),intent(in):: itx
    real(r_kind),intent(in):: dist1,crit1
    character(20),intent(in):: sis

    real(r_kind) crit

    if(.not. iuse)return

!   If using all data (no thinning), simply return to calling routine

    if(use_all)then
      if(iobs < itxmax)then
       iobs=iobs+1
       iobsout=iobs	
      else
       iuse = .false.
       write(6,*)' ndata > maxobs when reading data for ',sis,itxmax
      end if
      return
    end if

    crit=crit1*dist1

    if(icount(itx) == 0)then

!   Increment obs counter

      if(iobs < itxmax)then
       iobs=iobs+1
       iobsout=iobs
       score_crit(itx)= crit
       ibest_obs(itx) = iobs
       icount(itx)=icount(itx)+1
      else
       iuse = .false.
       write(6,*)' ndata > maxobs when reading data for ',sis,itxmax
      end if

    else if(crit < score_crit(itx))then
      score_crit(itx)= crit
      iobsout=ibest_obs(itx)
      icount(itx)=icount(itx)+1
    else
      iuse = .false.
    end if


    return
  end subroutine finalcheck



  subroutine destroygrids 9
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    destroygrids
!     prgmmr:    treadon     org: np23                date: 2002-10-17
!
! abstract:  This deallocate arrays used in thinning of satellite data.
!
! program history log:
!   2002-10-17  treadon
!   2004-06-22  treadon - update documentation
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$
    implicit none


!   If using all data (no thinning), the arrays following the
!   return are never allocated --> therefore nothing to deallocate
    if(use_all) return

    deallocate(mlon,glat,glon,hll)
    deallocate(icount)
    deallocate(ibest_obs)
    deallocate(score_crit)
    return
  end subroutine destroygrids


  subroutine destroy_sfc
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    destroy_sfc
!     prgmmr:    treadon     org: np23                date: 2005-10-06
!
! abstract:  This deallocate surface arrays
!
! program history log:
!   2005-10-06  treadon
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$
    implicit none

    deallocate(sst_full,sli_full,sno_full,isli_full)
	deallocate(u10_full,v10_full,zs_full)

    return
  end subroutine destroy_sfc


  subroutine indexx(n,arr,indx) 1,2
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    indexx
!     prgmmr:    treadon     org: np23                date: 2002-10-17
!
! abstract:  This routine indexes a sort key array, arr, such that
!            arr(indx(i),i=1,n) is the sort key in ascending order.
!
! program history log:
!   2002-10-17  treadon
!   2004-06-22  treadon - update documentation
!
!   input argument list:
!     n    - size of sort key array
!     arr  - sort key array
!
!   output argument list:
!     indx - index array
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$
    use kinds, only: r_kind,r_double,i_kind
    implicit none

    integer(i_kind) maxblock
    parameter (maxblock=1000)

    integer(i_kind) n,indx(n)
    real(r_kind) arr(n)

#ifdef ibm_sp
    real(r_kind),dimension(maxblock)::work
    real(r_double) oned
    real(r_kind) one

    one=1.; oned=1.
    if (digits(one)<digits(oned)) then
      call ssorts(arr,1,n,indx,work,maxblock)
    else
      call dsorts(arr,1,n,indx,work,maxblock)
    endif
#else
    integer(i_kind) m,nstack
    parameter (m=7,nstack=500)
    integer(i_kind) i,indxt,ir,itemp,j,jstack,k,l,istack(nstack)
    real(r_kind) a
    
    do j=1,n
       indx(j)=j
    end do
    jstack=0
    l=1
    ir=n
    
1   continue
    
    if(ir-l.lt.m)then
       do j=l+1,ir
          indxt=indx(j)
          a=arr(indxt)
          do i=j-1,l,-1
             if(arr(indx(i)).le.a)goto 2
             indx(i+1)=indx(i)
          end do
          i=l-1
2         continue
          indx(i+1)=indxt
       end do
       if(jstack.eq.0)return
       ir=istack(jstack)
       l=istack(jstack-1)
       jstack=jstack-2
       
    else
       k=(l+ir)/2
       itemp=indx(k)
       indx(k)=indx(l+1)
       indx(l+1)=itemp
       if(arr(indx(l)).gt.arr(indx(ir)))then
          itemp=indx(l)
          indx(l)=indx(ir)
          indx(ir)=itemp
       endif
       if(arr(indx(l+1)).gt.arr(indx(ir)))then
          itemp=indx(l+1)
          indx(l+1)=indx(ir)
          indx(ir)=itemp
       endif
       if(arr(indx(l)).gt.arr(indx(l+1)))then
          itemp=indx(l)
          indx(l)=indx(l+1)
          indx(l+1)=itemp
       endif
       i=l+1
       j=ir
       indxt=indx(l+1)
       a=arr(indxt)
3      continue
       i=i+1
       if(arr(indx(i)).lt.a)goto 3
       
4      continue
       j=j-1
       if(arr(indx(j)).gt.a)goto 4
       if(j.lt.i)goto 5
       itemp=indx(i)
       indx(i)=indx(j)
       indx(j)=itemp
       goto 3
       
5      continue
       indx(l+1)=indx(j)
       indx(j)=indxt
       jstack=jstack+2
       if(jstack.gt.nstack)then
          write(6,*)'INDEXX:  nstack=',nstack,' too small in indexx'
          call stop2(32)
       endif
       if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
       else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
       endif
    endif
    goto 1
#endif
  end subroutine indexx

end module satthin