; NAME: read_rcmggcm
;   
; PURPOSE: reads rcmggcm* output files
;   
; Inputs:
;   file - path to data file
;   grid_file - rcm_grid file (defaults to looking in same path)
;   fix_phi - convert phi to -180 to 180
;   max_it - stops reading of data at specified rcm simulation time
;   rcm_start - rcm start time (openggcm sim seconds) to add to output time
;
; Outputs:
;   data - array containing data ([nphi,ntheta,ntime])
;   it - array of time values (rcm simulation seconds)
;   phi - phi values on grid (0-360 where 0 is noon)
;   theta - theta values on grid
;
; Author:
;   W. Douglas Cramer
;
pro read_rcmggcm, file, data, it, phi, theta, $
  gridfile=gridfile, fix_phi=fix_phi, max_it=max_it,rcm_start=rcm_start

  if (~file_test(file)) then stop,'ERROR: file not found ('+file+')'
  if (n_elements(grid_file) ne 1) then begin
    gridfile=strmid(file,0,strpos(strmid(file,0,strlen(file)-1),'/',/reverse_search))+'/rcm_grid'
  endif
  if (n_elements(max_it) eq 0) then max_it=1.e10 ;huge

  fname=''
  it=-1
  data=-1

  dim=0

  openr,fileid,file,/get_lun;,/f77_unformatted
  iheader=lonarr(20)
  readu,fileid,iheader
  close,fileid
  free_lun,fileid

  ;nx=iheader[7]
  ;ny=iheader[8]
  ;nz=iheader[9]
  ;dim=1
  ;if (ny ne 0) then dim=dim+1 else ny=1
  ;if (nz ne 0) then dim=dim+1 else nz=1
  ;nelem=nx*ny*nz
  ;tmp_data=fltarr(60+nx*ny*nz)
  read_rcm_grid,gridfile,phi,theta
  nphi=n_elements(phi)
  ntheta=n_elements(theta)
  nelem=nphi*ntheta
  tmp_data=fltarr(60+nphi*ntheta)

  fi=file_info(file)
  maxrec=fi.size/nelem/4  ;assumes 4-byte floats
  it=lonarr(maxrec)
  data=fltarr(maxrec*nphi*ntheta)

  num_rec = 0

  openr,fileid,file,/get_lun;,/f77_unformatted

  ;loop through records
  while (~eof(fileid)) do begin

    readu,fileid,tmp_data

    ;get header data
    tmp_rcm_it=long(tmp_data[1],0) ;time since rcm start
    tmp_it=long(tmp_data[5],0) ;time since openggcm start
    if (tmp_it eq tmp_rcm_it) then begin
      ;old code prior to change for openggcm time
      if (n_elements(rcm_start) ne 1) then begin
        print,'old rcm code only has rcm time, assumed rcm start=5400'
        rcm_start=5400
      endif
      tmp_it=tmp_it+rcm_start
    endif

    ;if (num_rec eq 0) then begin
      ;it=tmp_it
      ;data=tmp_data[60:*]
    ;endif else begin
      ;it=[it,tmp_it]
      ;data=[data,tmp_data[60:*]]
    ;endelse
    it[num_rec]=tmp_it
    data[num_rec*nelem:(num_rec+1)*nelem-1]=tmp_data[60:*]

    num_rec++

    if (tmp_it ge max_it) then break
  endwhile
  it=it[0:num_rec-1]
  data=data[0:num_rec*nelem-1]

  close,fileid
  free_lun,fileid
  
  ;data=reform(data,nx,ny,nz,num_rec)
  data=reform(data,nphi,ntheta,num_rec)
  data=reform(data)

  if (keyword_set(fix_phi)) then begin
    f=where(phi ge 0.0 and phi lt 360.0)
    phi=phi[f]
    data=data[f,*,*]

    f=where(phi gt 180.0)
    phi[f]-=360.0

    fs=sort(phi)
    phi=phi[fs]
    data=data[fs,*,*]

    fu=uniq(phi)
    phi=phi[fu]
    data=data[fu,*,*]
  endif else begin
    f=where(phi[1:*] lt phi)
    if (f[0] ne -1) then begin
      phi[0:f[0]]-=360.0
      for i=1,n_elements(f)-1 do phi[f[i]+1:*]+=360.0
    endif
  endelse
end

;function rcm_fun,x
;  a = 50.0
;  b = 1.0 
;  xm = 0.001
;
;  IF (x le xm) THEN begin
;    val = (b-a)/xm*x^2/2.0+a*x
;  endif ELSE begin
;    val = ((a-b)*x^2/2.0+(b-a*xm)*x+xm/2.0*(a-b))/(1.0-xm)
;  ENDelse
;  val = val/(a+b)*2.0
;
;  return,val
;end