; NAME: read_jrrle_multi
;   
; PURPOSE: reads openggcm jrrle-format output files (diags2/combined).
;   
; Inputs:
;   file - path to data file
;   search_var - array of variables to retrieve (e.g., ['bx',...])
;                (defaults to all; narrowing search can speed up)
;
; Outputs:
;   data - array of structures containing data (e.g., data[].bx)
;   it - array of time values (simulation seconds)
;   time - array of time values (julday)
;   info - data line from file (not usually needed)
;
; Notes:
;   some code borrowed from probe/read_mhd2d.f
;
; Author:
;   W. Douglas Cramer
;
pro read_jrrle_multi, file, data, it, time, info, search_var=search_var

  if (~file_test(file)) then begin
    it=-1
    print,'ERROR: file not found ('+file+')'
  endif
  num_vars = n_elements(search_var)
  if (num_vars eq 0) then begin
    find_all=1 
    search_var=['']
  endif else begin
    find_all=0
  endelse
  num_vars_found = 0

  type=''
  line=''
  var=''
  it=0UL
  nx=0UL
  ny=0UL
  nz=0UL
  did=''
  n=0UL
  zmin=0.0
  zmax=0.0
  rid=0.0
  it=0UL
  cid=''
  i3 = intarr(64)
  a1 = fltarr(64)

  openr,fileid,file,/get_lun

  ;loop through variables
  repeat begin
    repeat begin
      readf,fileid,type  ; FIELD-2D-1, etc.
      ;print,type
    endrep until (strmid(type,0,5) eq 'FIELD' || eof(fileid))

    readf,fileid,var 
    var=strtrim(var,2)
    var_found = find_all or (where(var eq search_var) ge 0)
    ;print,var,var_found

    readf,fileid,line  ; time=...
    ;dummy=strsplit(line,/extract)
    ;modtime1=double(dummy[1])
    ;modtime2=double(dummy[2])
    ;acttime=dummy[3]

    dim = 0
    nrec=0UL
    switch (type) of
      'FIELD-1D-1': begin
                      dim = 1
                      readf,fileid,it,nx
                      nrec=nx
                      break
                    end
      'FIELD-2D-1': begin
                      dim = 2
                      readf,fileid,it,nx,ny
                      nrec=nx*ny
                      break
                    end
      'FIELD-3D-1': begin
                      dim = 3
                      readf,fileid,it,nx,ny,nz
                      nrec=nx*ny*nz
                      break
                    end
      else:         begin
                      close,fileid
                      free_lun,fileid
                      stop,'unknown field type ','<'+type+'>'
                    end
    endswitch

    ;rdn2 stuff
    repeat begin
      readf,fileid,did,n,zmin,zmax,rid,it,cid,format='(a4,i8,3e14.7,i8,a8)'
    endrep until (did eq 'WRN2')
        
    if (n ne nrec) then begin
      close,fileid
      free_lun,fileid
      stop,'ERROR: n ne nrec ',n,nrec
    endif
  
    if (var_found) then begin
      dat = fltarr(nrec)
    endif
  
    if (zmin eq zmax) then begin
      if (var_found) then dat[*] = zmin
    endif else begin
      q = 0UL
      dzi = (zmax-zmin)/4410.0
     
      for k=1UL,n,64 do begin
        nk = min([63,n-k])

        wrndec,fileid,rec1,nn
        if (nn ne nk) then begin
          close,fileid
          free_lun,fileid
          stop,'ERROR: nn ne nk ',nn,nk,n,k
        endif

        wrndec,fileid,rec2,nn
        if (nn ne nk) then begin
          close,fileid
          free_lun,fileid
          stop,'ERROR: nn ne nk ',nn,nk,n,k
        endif

;        if (var_found) then begin
;          for i=0,nk do begin
;            rec1[i] -= 33
;            rec2[i] -= 33
;            sig = 1.
;  
;            if (rec1[i] ge 47) then begin
;              sig = -1.
;              rec1[i] -= 47
;            endif
;
;            i3[i] = rec2[i] + 94*rec1[i]       
;            a1[i] = float(i3[i])
;            a1[i] = dzi*a1[i] + zmin
;            a1[i] = sig*exp(a1[i])
;
;            dat[q] = a1[i]
;            q += 1
;          endfor
;        endif
        if (var_found) then begin
          rec1[0:nk] -= 33
          rec2[0:nk] -= 33
          sig = fltarr(nk+1) + 1.
  
          f = where(rec1 ge 47)
          if (f[0] ne -1) then begin
            sig[f] = -1.
            rec1[f] -= 47
          endif

          i3 = rec2 + 94*rec1       
          a1 = float(i3)
          a1 = dzi*a1 + zmin
          a1 = sig*exp(a1)

          dat[q:q+nk] = a1
          q += nk + 1
        endif
      endfor
      if (var_found) then begin
        ;replace any illegal characters in variable name
        varlen = strlen(var)
        if (stregex(strmid(var,0,1),'[a-zA-Z]') eq -1) then strput,var,'Z',0
        for i=1,varlen-1 do begin
          if (stregex(strmid(var,i,1),'[a-zA-Z0-9_$]') eq -1) then begin
            strput,var,'_',i
          endif
        endfor
 
        if (n_elements(vars_found) eq 0) then begin
          vars_found=[var]
          info=line;[line]
          case dim of
            1: data = create_struct(var, reform(dat, nx))
            2: data = create_struct(var, reform(dat, nx, ny))
            3: data = create_struct(var, reform(dat, nx, ny, nz))
          endcase
        endif else begin
          vars_found=[vars_found,var]
          ;info=[info,line]
          case dim of
            1: data = create_struct(data, var, reform(dat, nx))
            2: data = create_struct(data, var, reform(dat, nx, ny))
            3: data = create_struct(data, var, reform(dat, nx, ny, nz))
          endcase
        endelse
        
        ;extract time from info
        dummy=strsplit(info,' ',/extract)
        if (n_elements(dummy) ge 4) then begin
          tvals=strsplit(dummy[3],':',/extract)
          tsec=float(tvals[5])
          tvals=fix(tvals[0:4])
          time=julday(tvals[1],tvals[2],tvals[0],tvals[3],tvals[4],tsec) 
        endif else begin
          time=0.0
        endelse

        num_vars_found = n_elements(vars_found)
      endif
    endelse

  endrep until (num_vars_found eq num_vars || eof(fileid))

  close,fileid
  free_lun,fileid

  ;if (num_vars_found ne num_vars) then begin
    ;dat = dat[0:nrec*num_vars_found-1]
  ;endif
end

pro wrndec, fileid, record, num
  record = intarr(64)
  dat = intarr(70)
  indata = ''

  readf,fileid,indata,format='(A)'
  hexval = byte(indata)
  len = n_elements(hexval)

  i = 0
  j = 0
  repeat begin
    ;reads,strmid(indata,j,1),hexval,format='(Z)'
    ;if (hexval[j] eq 32) then break
    if (hexval[j] gt 127) then begin
      nrep = hexval[j] - 170
      j += 1
      ;reads,strmid(indata,j,1),hexval,format='(Z)'
      for q=1,nrep do begin
        i += 1
        dat[i] = hexval[j]
      endfor
    endif else begin
      i += 1
      dat[i] = hexval[j]
    endelse
    j += 1
  endrep until (j gt min([66,len-1]))

  if (j gt 66) then begin
    close,fileid
    free_lun,fileid
    stop,'cant find end of encoded record'
  endif

  num = i-2
  if (num gt 63) then begin
    close,fileid
    free_lun,fileid
    stop,'num gt 63 ',indata
  endif

  checksum = 0
  for i=0,num do begin
    record[i] = dat[i+2]
    checksum += record[i]
  endfor

  if (dat[1] ne 33+(checksum mod 92)) then begin
    close,fileid
    free_lun,fileid
    stop,'checksum error', dat[1], checksum, 33+(checksum mod 92)
  endif

  return
end