;---------------------------------------------------------------------------
; This file contains following modules...
;    
;     PROC_STRING  : Breaks the string parameter  & stores in the structure. 
;     PROC_ZDEF    : Stores the level info in the appropriate structure.
;     PROC_DATE    : Stores the date info in the appropriate structure.
;     CHECK_TDEF   : processes and validates the date.
;     PROC_TIME    : stores the time info in the appropriate structure. 
;     PROC_VARS    : Processes variable information.
;     READ_STRUCT  : Reads the structure from one to another.
;     FORM_STRUCT  : Forms the huge attribute structure to be return.
;     READ_HDR     : Reads the header line by line and stores all the info in
;                    a huge structure.
;    
; LAST UPDATE: 21 July 2000    
;    
;     MODIFICATIONS: 
;       10 March 2000 DMA, CIRES/ETL   Close and free_lun statements added.
;       10 July 2000  DMA, CIRES/ETL   Set dat_type properly (not alway FLOAT) 
;       12 July 2000  DMA, CIRES/ETL   Added /nodata and /noscale flags 
;---------------------------------------------------------------------------

; Processing the space variables
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function proc_string, string, struct

common datastruct, x_linear, y_linear, z_levels, t_linear, $
mask,datatype,dataset,ttl, undefined, var_struct,header_info, no_of_var


space=' '
comp = StrParse(string,space,list)

if((comp eq 3) and (strupcase(list(1)) ne 'LINEAR')) then begin

	print,'Syntax error in space variable definition'
	return,0

endif else begin

	struct.gridpts=fix(list(0))
	struct.maptype=list(1)
	struct.startpt=float(list(2))
	struct.incr=float(list(3))
endelse

return,1
end
;;;;;;;;;;;; Processing  of zdef ;;;;;;;;;;;;;;;;;

function proc_zdef,string,struct

space=' '
comp=StrParse(string,space,list)

struct.no_levels=fix(list(0))
struct.maptype ='LEVELS'
if strupcase(list(1)) eq 'LINEAR' then begin
        st_pt=fix(list(2))
        inc=fix(list(3))
	for loop=0, struct.no_levels-1,1 do begin
            struct.lev(loop)= st_pt+(inc * loop)
        endfor
endif else begin
    if strupcase(list(1)) eq 'LEVELS' then begin
	for loop=0,struct.no_levels-1,1 do begin
           struct.lev(loop)=fix(list(loop+2))
        endfor
    endif else begin
        print,'Wrong format in z-definition'
        return,0
    endelse
endelse

return,1
end
        
;;;;;;;;;;;; Processing date ;;;;;;;;;;;;;;;;;;;

pro proc_date, string, flag


common datastruct, x_linear, y_linear, z_levels, t_linear, $
mask,datatype,dataset,ttl, undefined, var_struct,header_info, no_of_vars


temp=bytarr(10)
len = strlen(string)
temp_string = byte(string)
i=0

if (len eq 5) or (len eq 7) then begin
	num1=1
endif else begin
	while (temp_string(i) le 57) and (temp_string(i) ge 48)  do begin  
		;loop untill an alphabet is hit
		temp(i) = temp_string(i)
		i = i+1
	endwhile
	num1 = fix(string(temp))
endelse

x=i
while (temp_string(i) ge 65)and (temp_string(i) le 122)  do begin        
; loop untill a number is hit
	temp(i-x) = temp_string(i)
        if ( (i+1) lt len ) then begin
             i = i+1
        endif else begin
             goto, endloop1
        endelse
endwhile  
endloop1: str = string(temp)

if flag eq 'date' then begin
	x=i
	while (temp_string(i) le 57) and (temp_string(i) ge 48) do begin   
		temp(i-x) = temp_string(i)
		if ( (i+1) lt len ) then begin
			i = i+1
		endif else begin
			goto, endloop2
		endelse
	endwhile
endloop2: num3 = fix(string(temp))
	  t_linear.day = num1
	  t_linear.mon = strupcase(str)
	  t_linear.year = num3
endif else begin
	t_linear.incr_by = num1
	t_linear.incr_unit = strupcase(str)
endelse

end

;;;;;;;;;;;;;;;; Check for erros in time def ;;;;;;;;;;;;;;;;;;;;;

function check_tdef

common datastruct, x_linear, y_linear, z_levels, t_linear, $
mask,datatype,dataset,ttl, undefined, var_struct,header_info, no_of_vars


;help, /structure, t_linear
;;;; Time checking....

if ( (t_linear.hour lt 0) or (t_linear.hour gt 24) ) then return,'hours'
if ( (t_linear.minutes lt 0) or (t_linear.minutes gt 59) ) $
    then return,'minutes'

;;;; Date checking....

if ( t_linear.year mod 4 ne 0 ) then leap_feb = 28 else leap_feb = 29

if ( ((t_linear.mon ne 'JAN') or $
     (t_linear.day lt 1) or (t_linear.day gt 31)) and $
     ((t_linear.mon ne 'FEB') or $
     (t_linear.day lt 1) or (t_linear.day gt leap_feb)) and $
     ((t_linear.mon ne 'MAR') or $
     (t_linear.day lt 1) or (t_linear.day gt 31)) and $
     ((t_linear.mon ne 'APR') or $
     (t_linear.day lt 1) or (t_linear.day gt 30)) and $
     ((t_linear.mon ne 'MAY') or $
     (t_linear.day lt 1) or (t_linear.day gt 31)) and $
     ((t_linear.mon ne 'JUN') or $
     (t_linear.day lt 1) or (t_linear.day gt 30)) and $
     ((t_linear.mon ne 'JUL') or $
     (t_linear.day lt 1) or (t_linear.day gt 31)) and $
     ((t_linear.mon ne 'AUG') or $
     (t_linear.day lt 1) or (t_linear.day gt 31)) and $
     ((t_linear.mon ne 'SEP') or $
     (t_linear.day lt 1) or (t_linear.day gt 30)) and $
     ((t_linear.mon ne 'OCT') or $
     (t_linear.day lt 1) or (t_linear.day gt 31)) and $
     ((t_linear.mon ne 'NOV') or $
     (t_linear.day lt 1) or (t_linear.day gt 30)) and $
     ((t_linear.mon ne 'DEC') or $
     (t_linear.day lt 1) or (t_linear.day gt 31)) )  then return,'date'

if (t_linear.year lt 1950) or (t_linear.year gt 2049) or $
   (t_linear.year lt 50) then return,'year'

if (t_linear.incr_by lt 0) or (t_linear.incr_by gt 99) then $
                 return,'tdef increments'


;;;; Time increment checking...

if ((t_linear.incr_unit ne 'MN') and (t_linear.incr_unit ne 'HR') and $
   (t_linear.incr_unit ne 'DY') and (t_linear.incr_unit ne 'MO') and $
   (t_linear.incr_unit ne 'YR')) then return,'tdef increment units'

return,1

end
;;;;;;;;;;;;; Processing of time ;;;;;;;;;;;;;;;;

function proc_time, value

common datastruct, x_linear, y_linear, z_levels, t_linear, $
mask,datatype,dataset,ttl, undefined, var_struct,header_info, no_of_vars


space =' '
components = StrParse(value,space,list)
t_linear.no_times= fix(list(0))
t_linear.maptype = list(1)

pos = strpos(list(2), ':')
zpos = strpos(list(2), 'Z')
if pos ne -1 then begin
        t_linear.hour = fix(strmid(list(2), 0, pos))
	t_linear.minutes = fix( strmid(list(2), pos+1, $
                                (strpos(list(2),'Z') -1)) )
endif else begin
	if (zpos ne -1) then begin
	        t_linear.hour = fix( strmid(list(2), pos+1, $
                                strpos(list(2),'Z')) )
	endif else begin
        	t_linear.hour = 0
	endelse
endelse


pos = strpos(list(2), 'Z')
temp = strmid( list(2), pos+1, strlen(list(2)) )
proc_date, temp, 'date'
temp = list(3)
proc_date, temp,'incr'

ret = check_tdef()
if (ret  ne '1') then begin
	print, 'Wrong format in ', ret
	return,0
endif else return,1
end

;;;;;;;;;;;; Processing of variables ;;;;;;;;;;;;;;;
function proc_vars, string, st_index 
common datastruct, x_linear, y_linear, z_levels, t_linear, $
mask,datatype,dataset,ttl, undefined, var_struct,header_info, no_of_vars


space=' ' & temp=''
no=StrParse(string,space,list)
if (no lt 3) then begin
	print,'Wrong format in variable definition'
	return,0
endif else begin
	var_struct(st_index).var_name = list(0)
	var_struct(st_index).var_scale = float(list(1))
        valid=0
        ON_IOERROR, not_int
        i=fix(list[1])
	var_struct(st_index).var_offset = float(list(2))
        valid=1
        not_int: IF NOT valid THEN BEGIN
                   ind=2
                   var_struct(st_index).var_offset=0.0
                 ENDIF ELSE BEGIN
                   ind=3
                 ENDELSE
	var_struct(st_index).var_units = strtrim(list(ind))
	for i=ind+1, no, 1 do begin
		temp = temp + ' ' + list(i)
	endfor
	var_struct(st_index).subtitle = strtrim(temp , 2)
endelse
return,1

end

;;;;;;;;;;;; read structure routine ;;;;;;;;;;;;;;
pro read_struct,dest,scr 


if n_tags(dest) ne n_tags(scr) then print,'Reading wrong structure'
for i=0, n_tags(scr)-1 do begin
	temp=scr.(i)             ;Direct assignment is not possible with struct
	dest.(i)=temp
endfor

end

;;;;;;;;;;;; Forming the structure ;;;;;;;;;;;;;;;;

pro form_struct, attr_struct, flag, noscale=noscale 

; Note: data order reversed if flag is 1, else not reversed 

common datastruct, x_linear, y_linear, z_levels, t_linear, $
mask,datatype,dataset,ttl, undefined, var_struct,header_info, no_of_vars

temp_struct=attr_struct.x_lnr
read_struct,temp_struct,x_linear
attr_struct.x_lnr=temp_struct

temp_struct=attr_struct.y_lnr
read_struct,temp_struct,y_linear
if ( temp_struct.incr LT 0 and flag eq 1 ) then begin
  ;attr_struct.order_info=1
  temp_struct.startpt=temp_struct.startpt $
                      + temp_struct.incr*(temp_struct.gridpts-1)
  temp_struct.incr = -temp_struct.incr
endif 
attr_struct.y_lnr=temp_struct

temp_struct=attr_struct.t_lnr
read_struct,temp_struct,t_linear
attr_struct.t_lnr=temp_struct

temp_struct=attr_struct.z_lvl
read_struct,temp_struct,z_levels
attr_struct.z_lvl=temp_struct

for i=0,19 do begin
	temp_struct=attr_struct.var_strc(i)
	read_struct,temp_struct,var_struct(i)
	attr_struct.var_strc(i)=temp_struct
endfor

attr_struct.dat_set=dataset

if ( keyword_set (noscale) ) then begin 
  case datatype of
    'INTEGER' : attr_struct.dat_type='INTEGER' 
    'LONG'    : attr_struct.dat_type='LONG' 
    'FLOAT'   : attr_struct.dat_type='FLOAT' 
    'BYTE'    : attr_struct.dat_type='BYTE' 
  endcase
endif else begin 
  attr_struct.dat_type='FLOAT'
endelse

attr_struct.msk_val=mask
attr_struct.ttl_str=ttl
attr_struct.unk=undefined
attr_struct.h_flag=header_info
attr_struct.no_vars=no_of_vars
attr_struct.scl_ftr = 1.0
if attr_struct.var_strc(0).var_units eq '0' then begin
	attr_struct.units=''
endif else begin
	attr_struct.units=attr_struct.var_strc(0).var_units
endelse

if attr_struct.x_lnr.gridpts EQ 1 then begin
  attr_struct.x_ext=strtrim(string(attr_struct.x_lnr.startpt,format='(f10.2)'))
endif
if attr_struct.y_lnr.gridpts EQ 1 then begin
  attr_struct.y_ext=strtrim(string(attr_struct.y_lnr.startpt,format='(f10.2)'))
endif
if attr_struct.z_lvl.no_levels GT 1 then begin
   attr_struct.z_ext=strtrim(string(attr_struct.z_lvl.lev(0)),2)
endif else begin
   attr_struct.z_ext=strtrim(string(1),2)
endelse

attr_struct.v_ext=attr_struct.var_strc(0).var_name
attr_struct.t_ext=strtrim(string(attr_struct.t_lnr.day),2)+' '+attr_struct.t_lnr.mon+' '+strtrim(string(attr_struct.t_lnr.year),2)
end

;-------------------------------------------------------------------------------
;+
; NAME:
;       READ_HDR
; PURPOSE:
;       Read the grads header and store the info in the grads header in the 
;       attribute structure.
; CATEGORY:
; CALLING SEQUENCE:
;       var = read_hdr()           -> Displays description of read_hdr
;       var = read_hdr('header file')   
; INPUTS:
;       The function parameter should be a scalar string
; COMMON BLOCKS:
;       datastruct -- It is local to this file and should not be used outside.
; KEYWORD PARAMETERS:
; OUTPUTS:
;       A structure containing the information from the grads header.
; SIDE EFFECTS:
;      If the data array is stored in the reverse order (N to S) then the
;      array is reversed and reformed to be of the order (S to N). 
; AUTHOR:  J. Kirani  - June 5, 1995 (jvk@cdc.noaa.gov)
;
;-------------------------------------------------------------------------------

;;;;;;;;;;;; main file read function ;;;;;;;;;;;;;;

function read_hdr, filename, swap=swap, nodata=nodata, noscale=noscale

common datastruct, x_linear, y_linear, z_levels, t_linear, $
mask,datatype,dataset,ttl, undefined, var_struct,header_info, no_of_vars

if n_params() eq 0 then begin
  print,' Calling Sequence is struct=read_hdr(header filename)'
  print,' Keywords:'
  print,'   /swap    -> Swaps bytes for data files from linux machines'
  print,'   /nodata  -> Read header info only, do not load data'
  print,'   /noscale -> Do not reset scales and offsets for each variable'
  goto, stop
endif

header_info = 0 
rd_order = 0
openr, d1, filename, /get_lun, error=err 
if (err ne 0) then begin
        print,'Non existent file'
	print,'Calling Sequence is struct=read_hdr(header filename)'
	goto,stop
endif

;;;;;;;;;; local variables ;;;;;;;;;;;;;;;;;;;;;;
header_info = 1
var='' & data='' & space = ' '
;dataset = '' & ttl='' & undefined='' & line='' & undefined=-999.0
dataset = '' & ttl='' & line='' & undefined=-999.0
mask=-998.0 & datatype='FLOAT' & scale_factor=1.0 & offset=0.0

;;;;;;;;;;;;; Structure templates ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
; NOTE: The number of variables and number of levels that can be stored 
;       in the attribute structure is limited to 200.
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

x_linear = {lin, gridpts:0, maptype:'', startpt:0.0, incr:0.0}
y_linear = {lin, gridpts:0, maptype:'', startpt:0.0, incr:0.0}
t_linear = {time, no_times:0, maptype:'', hour:0, minutes:0, $
	    day:0, mon:'', year:0, incr_by:0, incr_unit:''}

z_levels = {levl, no_levels:1, maptype:'', lev:intarr(100)}
var_struct=replicate({v_st, var_name:'', var_scale:1.0, var_offset:0.0, $
                      var_units:'', subtitle:''}, 200)

while not eof(d1) do begin
	readf,d1,'A',line
	no_of_elements = StrParse(line,space,list)
	var = strupcase(list(0))

;;;;;;;;; Getting rid of inline comments

	for i=1, no_of_elements, 1 do begin
                if strpos(list(i),'*') eq 0 then begin
                  j=i
                  while j le no_of_elements do begin
                       list(j) = ''
                       j=j+1
                  endwhile
                endif 
		data = data + ' ' + list(i)
	endfor
	data = strtrim(data,2)
	case var of
	     'DSET'  : dataset = data 
  	     'TITLE' : ttl = data
	     'UNDEF' : undefined = float(data)
	     'XDEF'  : if(proc_string(data,x_linear) eq 0) then goto, stop
	     'YDEF'  : if(proc_string(data,y_linear) eq 0) then goto, stop
	     'ZDEF'  : if(proc_zdef(data,z_levels) eq 0) then goto, stop

	     'TDEF'  : if(proc_time(data) eq 0) then goto, stop
	     'VARS'  :begin 
			no_of_vars=fix(data)
		        for i=0,no_of_vars-1 do begin 
		        	readf,d1,'A',line
		        	ret = proc_vars(line, i)
		        	if (ret eq 0) then goto, stop
	                endfor	
		        readf,d1,'A',line
		        if (strtrim(strupcase(line),2) ne 'ENDVARS') then begin
                                print, line,ret
				print,'Incorrect end of header file'
				goto,stop
		       endif
		      end

	     'MDEF'  : mask = float(data) 
	     'DTYPE' : datatype = (strupcase(data))
	     'SCALE' : scale_factor=float(data)
	     'OFFSET': offset=float(data)
	     else   : 
	endcase
	data = ''
endwhile

; Override scale and offset factors for each variable if SCALE or OFFSET is set
; and if /noscale keyword is not provided

if ( not keyword_set(noscale) ) then begin  
  if (scale_factor NE 1.0) then begin
    for i=0,no_of_vars-1 do begin 
      var_struct[i].var_scale=fix(1.0/scale_factor)
    endfor
  endif
  if (offset NE 0.0) then begin
    for i=0,no_of_vars-1 do begin 
      var_struct[i].var_offset=offset
    endfor
  endif
endif

close, d1
free_lun, d1

;;;;;;;;;;;;; Reading the binary file into an array ;;;;;;;

if ( not keyword_set(nodata) ) then begin  

  ; Extract the path from the filename 
  slashloc = strpos(filename,'/',/reverse_search)
  path = strmid(filename,0,slashloc)

  s=findfile(dataset,count=cnt)
  if (cnt EQ 0) then begin
    ; Get the binary data file name (minus path)
    slashloc = strpos(dataset,'/',/reverse_search)
    binname = strmid(dataset,slashloc+1,strlen(dataset)-slashloc)

    ; Build the complete binary file name
    dataset = path+'/'+binname
  endif

  ;if (path eq "") then begin
  ;  dataset = binname
  ;endif

  if (n_elements(swap) GT 0) then begin
    if (strpos(dataset,'.gz') EQ (strlen(dataset)-3)) then begin
      openr, d, dataset, /get_lun, error=err, /compress, /swap_endian
    endif else begin
      openr, d, dataset, /get_lun, error=err, /swap_endian
    endelse
  endif else begin
    if (strpos(dataset,'.gz') EQ (strlen(dataset)-3)) then begin
      openr, d, dataset, /get_lun, error=err, /compress
    endif else begin
      openr, d, dataset, /get_lun, error=err
    endelse
  endelse
  if (err ne 0) then begin
     print,'Cant open the binary data file ',dataset
  endif else begin
     print,'reading bin file'
     case datatype of
        'INTEGER' : temp_arr=intarr(x_linear.gridpts,y_linear.gridpts,$
                    z_levels.no_levels,t_linear.no_times, no_of_vars)

        'LONG'    : temp_arr=lonarr(x_linear.gridpts,y_linear.gridpts,$
                    z_levels.no_levels,t_linear.no_times, no_of_vars)

	'FLOAT'   : temp_arr=fltarr(x_linear.gridpts,y_linear.gridpts,$
                    z_levels.no_levels, t_linear.no_times, no_of_vars)

	'BYTE'    : temp_arr=bytarr(x_linear.gridpts,y_linear.gridpts,$
                    z_levels.no_levels, t_linear.no_times, no_of_vars)
     endcase

     if ( keyword_set (noscale) ) then begin 
       attr_struct = {x_lnr:{lin}, y_lnr:{lin}, t_lnr:{time}, z_lvl:{levl}, $
                      var_strc:replicate({v_st},200), dat_set:'', dat_type:'', $
                      msk_val:0.0,units:'',scl_ftr:1.0, $
                      ttl_str:'', unk:0.0, h_flag:0, order_info:0, no_vars:0,$
                      x_ext:'', y_ext:'', z_ext:'', t_ext:'', v_ext:'',$
                      a:temp_arr} 
     endif else begin
       attr_struct = {x_lnr:{lin}, y_lnr:{lin}, t_lnr:{time}, z_lvl:{levl}, $
                      var_strc:replicate({v_st},200), dat_set:'', dat_type:'', $
                      msk_val:0.0,units:'',scl_ftr:1.0, $
                      ttl_str:'', unk:0.0, h_flag:0, order_info:0, no_vars:0, $
                      x_ext:'', y_ext:'', z_ext:'', t_ext:'', v_ext:'', $
                      a:fltarr(x_linear.gridpts,y_linear.gridpts, $
                      z_levels.no_levels, t_linear.no_times, no_of_vars)}
     endelse
  endelse

  on_ioerror, End_of_file
  readu, d, temp_arr
  free_lun, d
  undex=where((temp_arr EQ undefined),ucnt)
  mndex=where((temp_arr EQ mask),mcnt)
  temp_arr=float(temp_arr)
  for i=0,no_of_vars-1 do begin 
    if ( not keyword_set ( noscale ) ) then begin 
      temp_arr[*,*,*,*,i]=temp_arr[*,*,*,*,i] / var_struct[i].var_scale $
                        - var_struct[i].var_offset
      var_struct[i].var_scale = 1.0
      var_struct[i].var_offset = 0.0
    endif
  endfor
  if (ucnt GT 0) then temp_arr(undex)=undefined
  if (mcnt GT 0) then temp_arr(mndex)=mask

; The following code reorients the data to be from South to North if neccessary

  if y_linear.incr lt 0 then begin
     print,'NOTE: DATA ARRAY IS REORIENTED TO BE FROM SOUTH TO NORTH'
     size_out = size(temp_arr)
     if (size_out(0) GT 2) then begin
     s3=size_out(size_out(0)+2)/(size_out(1)* size_out(2))
     temp_arr=reverse(reform(temp_arr,size_out(1),size_out(2),s3),2)
     if (size_out(0) EQ 4) then $
         temp_arr=reform(temp_arr,size_out(1),size_out(2),size_out(3),$
                         size_out(4))
     if (size_out(0) EQ 5) then $
       temp_arr=reform(temp_arr,size_out(1),size_out(2),size_out(3),$
                       size_out(4),size_out(5))
     endif else begin
      temp_arr=reverse(temp_arr,2)
     endelse
  endif
  attr_struct.a = temp_arr
  print,'Bin file read completed'
  if ( keyword_set (noscale) ) then form_struct, attr_struct, 1, /noscale $ 
  else form_struct, attr_struct, 1 

; If the nodata flag was passed, return attr_struct without data

endif else begin 

  attr_struct = {x_lnr:{lin},y_lnr:{lin},$
                 t_lnr:{time},z_lvl:{levl},$
                 var_strc:replicate({v_st},200),dat_set:'',dat_type:'',$
                 msk_val:0.0,units:'',scl_ftr:1.0,$
                 ttl_str:'',unk:0.0,h_flag:0,order_info:0,no_vars:0,$
                 x_ext:'',y_ext:'',z_ext:'',t_ext:'',v_ext:'',$
                 a:fltarr(x_linear.gridpts,y_linear.gridpts,$
                 z_levels.no_levels, t_linear.no_times, no_of_vars)}
  form_struct, attr_struct, 0, /noscale 

endelse

return, attr_struct
End_of_file:    if !err_string ne '' then print,'Bin file read incomplete' & print, 'Float data type expected' 

stop : if n_elements(d) ne 0 then begin 
         close, d 
         free_lun, d
       endif
return,0
end