! Description: 	reads in and plots lat and lon of data pts along TP tracks
!               reads in sea surface height anomaly, ssha, data
!		fills in the time [aviso day] variable via linear interpolation
!		converts ssha from mm to m
!		compares the filled, raw and linear interpolated ssha
!		saves to netCDF the important data:
! -----------------------------------------------------------------------
! by Rob Scot, June 2004:
! -----------------------------------------------------------------------
! MUST SPECIFY:
!  file info:
! -----------------------------------------------------------------------
LET ncycles = 360           ! some are missing
LET ntracks = 39
LET npts = 796  !214    796

DEFINE SYMBOL raw_ssha = cyc_tr_pt_cyc_ssha_day_tp_xxC.dat
DEFINE SYMBOL raw_latlon = lat_lon_pt_track_tp_xxC.dat

DEFINE SYMBOL out_file = ssha_lat_lon_raw_xxC.cdf


SET MEMORY/SIZE=240

!**********************************************************************
! Grid  info:
! -----------------------------------------------------------------------

define axis/z=1:`npts`:1 zaxis_pt
define axis/Y=1:`ntracks`:1 yaxis_track
define grid/z=zaxis_pt/y=yaxis_track grid_pt_track
file/G=grid_pt_track/var=track,pt,microlat,microlon,land/ORDER=ZY \
    /columns=5 ($raw_latlon)

! -----------------------------------------------------------------------

SET VARIABLE/bad = 32767 microlat
SET VARIABLE/bad = 32767 microlon
SET VARIABLE/TITLE="land mask"/BAD=32767 land
! convert units to degrees

LET/Title="point latitude"/units=latitude  lat_tmp = microlat[D=($raw_latlon)]/1E6
LET/Title="point lonitude"/units=longitude lon_tmp = microlon[D=($raw_latlon)]/1E6

! reveal the TP tracks in the Kuroshio Extension region!

plot/vs  lon_tmp,lat_tmp

PAUSE "enter to continue"
!**********************************************************************
! compare raw and filled data:
! ---------------------------------------------------------------------

!----------------------------------------------------------------------
! read ssha and time (in Aviso days) data, order: point, cycle, track
!----------------------------------------------------------------------
DEFINE AXIS/t=1:`ncycles`:1/UNITS="cycle index" tax_aviso_cycno
DEFINE GRID/LIKE=grid_pt_track/t=tax_aviso_cycno gr_pt_tr_cycno
FILE/VAR="index,trackno,ptno,tp_cycno,ssha_raw,aviso_day_raw,ssha_filled" \
    /COLUMNS=7 /GRID=gr_pt_tr_cycno/ORDER=ZTY ($raw_ssha)

SET VARIABLE/UNITS=mm/TITLE="sea level height anomaly"/BAD=32767 ssha_raw
SET VARIABLE/UNITS=mm/TITLE="sea level height anomaly"/BAD=32767 ssha_filled
SET VARIABLE/UNITS=days/TITLE="Days since Jan 1 1950"/BAD=32767 aviso_day_raw
SET VARIABLE/TITLE="cycle number"/BAD=32767 tp_cycno

! convert to m
LET ssha_tmp = ssha_filled[D=($raw_ssha)]/1000


SET VARIABLE/TITLE="sea level height anomaly" ssha_tmp
LET aviso_day_tmp = aviso_day_raw[D=($raw_ssha),t=@FLN]
LET cycle_tmp = tp_cycno[D=($raw_ssha),t=@FLN]

SAVE/FILE=($out_file)/FORMAT=cdf/CLOBBER ssha_tmp,aviso_day_tmp,\
       cycle_tmp,lat_tmp,lon_tmp, land[D=($raw_latlon)]


!**********************************************************************