;|-----------------------------------------------------------------|
;|                                                                 |
;|    Subroutine WATS                                              |
;|                                                                 |
;|    extracts basic parameters from Wdata.                        |
;|                                                                 |
;|-----------------------------------------------------------------|


pro wats


@gia_var.inc


; Now you can use Neu_Data and Columns in this code

neu_orbitnum=intarr(wcolumns)
nflag = intarr(wcolumns)
neu_vel=fltarr(wcolumns)
sc_vel=fltarr(wcolumns)
neu_flag_row = 2
neu_vel_row = 11
;sc_vel_row = 7

neu_orbitnum = fix(wdata(13,0:wcolumns-1))
neu_vel = wdata(neu_vel_row,0:wcolumns-1)
;sc_vel = wdata(sc_vel_row,0:wcolumns-1)
nflag = fix(wdata(neu_flag_row,0:wcolumns-1))/100
Orbitnumber=fix(idata(4,0))

; This next section does the same thing (keep it for future ref)

;neu_vel=EXTRAC(wdata, neu_vel_row, 0, 1, wcolumns)
;nflag=fix(EXTRAC(wdata, neu_flag_row, 0, 1, wcolumns)) /100
;neu_orbitnum=fix(EXTRAC(wdata, 13, 0, 1, wcolumns))

; Get Horizontal and Vertical Components of Neutral Velocity

tneu_hvel = FLTARR(wcolumns)
tneu_vvel = FLTARR(wcolumns)
tneu_schvel = FLTARR(wcolumns)
tneu_scvvel = FLTARR(wcolumns)
tneu_hut = FLTARR(wcolumns)
tneu_vut = FLTARR(wcolumns)
tneu_hglat = FLTARR(wcolumns)
tneu_hilat = FLTARR(wcolumns)
tneu_vglat = FLTARR(wcolumns)
tneu_vilat = FLTARR(wcolumns)

max1=-1.e3
min1=+1.e3
;tmax=1.
;tmin=1.e10

hcnt =1
vcnt =1

;a=""
for j=0,(wcolumns-1) do begin

;neutral_vel(j)=wdata(7,j)
;nflag(j)=fix((wdata(2,j))/100)



 if (neu_orbitnum(j) eq orbitnumber) then begin

   ; Set boundaries
        if neu_vel(j) ne 9999. then begin
		if neu_vel(j) gt max1 then max1=neu_vel(j)
		if neu_vel(j) lt min1 then min1=neu_vel(j)
		;if wdata(1,j) gt tmax then tmax=wdata(1,j)
		;if wdata(1,j) lt tmin then tmin=wdata(1,j)

     ;Horizontal neutral velocity

		if ((nflag(j) eq 3) or (nflag(j) eq 4)) and (neu_vel(j) gt -9000) then begin
        		tneu_hvel(hcnt)=neu_vel(j)
        		tneu_schvel(hcnt)=sc_vel(j)
        		tneu_hglat(hcnt)=wdata(15,j) 
        		tneu_hilat(hcnt)= wdata(20,j)
        		tneu_hut(hcnt)=wdata(1,j)    ;/3600000.
        		if abs(tneu_hglat(hcnt)) lt abs(tneu_hglat(hcnt-1)) then begin 
        			tneu_hvel(hcnt) =- tneu_hvel(hcnt) 
        			endif
        		;print,hcnt
      			hcnt=hcnt +1
      			endif 
     
     ;Vertical neutral velocity

     		if (nflag(j) eq 5) or (nflag(j) eq 6) then begin
        		tneu_vvel(vcnt)=neu_vel(j)
        		;tneu_scvvel(vcnt)=sc_vel(j)
       			tneu_vglat(vcnt)=wdata(15,j) 
		        tneu_vilat(vcnt)= wdata(20,j)
		        tneu_vut(vcnt)=wdata(1,j)   ;/3600000.
           		if abs(tneu_vglat(vcnt)) lt abs(tneu_vglat(vcnt-1)) then begin 
               			tneu_vvel(vcnt) =- tneu_vvel(vcnt) 
           			endif
        		;print,vcnt
        		vcnt=vcnt +1
     			endif 
		endif


 endif

endfor

;
;clearing up the zeros at the end
;sometimes the neutral data is not available so the if statement is necessary
;

if vcnt gt 1 then begin
  neu_hvel=tneu_hvel(1:hcnt-1)
  neu_schvel=tneu_schvel(1:hcnt-1)
  neu_hut=tneu_hut(1:hcnt-1)
  neu_hilat=tneu_hilat(1:hcnt-1)
  neu_hglat=tneu_hglat(1:hcnt-1)
endif else if vcnt eq 1 then begin 
  neu_hvel= tneu_hvel
  neu_schvel= tneu_hvel
  neu_hut = tneu_hut
  neu_hglat = tneu_hglat
  neu_hilat = tneu_hilat
  ;hcnt = wcolumns
endif 

if vcnt gt 1 then begin
  neu_vvel=tneu_vvel(1:vcnt-1)
  neu_scvvel=tneu_scvvel(1:vcnt-1)
  neu_vut=tneu_vut(1:vcnt-1)
  neu_vilat=tneu_vilat(1:vcnt-1)
  neu_vglat=tneu_vglat(1:vcnt-1)
endif else if vcnt eq 1 then begin
  neu_hvel= tneu_hvel
  neu_scvvel=tneu_scvvel
  neu_hut = tneu_hut
  neu_vilat=tneu_vilat
  neu_vglat=tneu_vglat
  ;hcnt = wcolumns
endif

end;--------------------------------------------------------------------------------------