set MODE VERIFY
!     
! ******************************************************************** 
! ef_sort_demo.jnl *acm* 9/98 - EF sorting and sampling demo
!                  *ACM* 5/00 -  correct the order of args to sample*

! Description: examples of the sorting external functions for FERRET.
! SORTI, SORTJ, SORTK, SORTL to obtain lists of indices of sorted data
! SAMPLEI, SAMPLEJ, SAMPLEK, SAMPLEL to sample data at a list of indices
 
CANCEL REGION

!         This demonstration/tutorial will introduce the FERRET
!         external functions for sorting data.
 
! ******************************************************************** 
MESSAGE

!         A first simple example of sorting: Sort a list of data.

USE coads_climatology
SET REGION/x=19e/y=77n/z=0/l=1:6
LIST sst

!         The sequence of calls is:  Compute the sorted indices, then sample 
!         the data using the new order of indices.  This is COADS sea surface 
!         data at a location where some of the data is missing.
! ******************************************************************** 
MESSAGE

!       Show the indices and the sorted temperatures

LET tsorted_indices = sortl(sst[l=1:6])

LET tsorted_sst = samplel(sst[l=1:6], tsorted_indices)
LIST tsorted_indices, tsorted_sst

!       We can also use the constant-array notation to pick out a couple
!       of points:

list samplel(sst, {1,5})



! ******************************************************************** 
MESSAGE

!         A more involved example:  What is the Sea Surface 
!         Temperature in a region 3 months after the strongest winds 
!         elsewhere?  Use the COADS monthly climatology data set.

CAN REGION 
SET DATA coads_climatology

!         Choose a region in the Pacific, and define the westerly_wind 
!         variable to be the monthly west winds averaged over this region.
!         Sort the westerly winds, lowest to highest. "sorted_indices3" 
!         is the indices 3 months after.

LET westerly_wind = uwnd[x=160e:180@ave,y=5s:5n@ave]
LET sst_e = sst[x=180:80w@ave,y=5s:5n@ave]

LET sorted_indices = sortl(westerly_wind)
LET sorted_indices3 = MOD((sorted_indices + 2), 12) + 1

!         We will plot sorted winds vs the SST's, 3 months later, 
!         and in another region.   First sample the westerly winds 
!         by their sorted indices.  Then order the SST's according  
!         to the sorted wind indices plus 3 months.  

LET wwe_by_wwe = samplel(westerly_wind, sorted_indices)
LET sst_by_wwe = samplel(sst_e, sorted_indices3)

!         The SST's of interest are to the east of our westerly_wind 
!         region.  Get the number of valid data in the sort. 

LET leng = wwe_by_wwe[l=@NGD]
! ******************************************************************** 
MESSAGE

!          Make scatter plots: sampled winds vs SST, and mark the
!          highest winds with a symbol.  The computation takes a moment.

PLOT/VS/TITLE="SST 3 months after High Westerly Winds"/SET_UP/l=1:`leng` wwe_by_wwe, sst_by_wwe
ppl xaxis -3.3, -1.9, 0.1
ppl yaxis, 25.4, 27.4, 0.2
ppl xlab "Sorted Westerly Winds from 160E to 180"
ppl ylab "SST from 180 to 80W 3 months after Westerly winds"
ppl plot

PLOT/VS/OVER/SYMBOL=18/NOLABELS/l=`leng-3`:`leng` wwe_by_wwe, sst_by_wwe

!         done