;
pro wrdata_1d, array,id1,fname,unit
printf,unit,fname
nvline = 5	; nvline values per line
nlines = id1/nvline
nextra = id1 mod nvline
if nextra gt 0 then nlines = nlines+1
format = "(5x,'+ '," + strcompress(string(nvline),/remove_all) + $
  "(e12.4,','))"
fline = fltarr(nvline)
for i=1,nlines-1 do begin
  for ii=0,nvline-1 do fline(ii) = array((i-1)*nvline+ii) 
  printf,unit,format=format,fline
end
if nextra eq 0 then begin
  format = "(5x,'+ '," + strcompress(string(nvline-1),/remove_all)+$
    "(e12.4,','),e12.4,'/')"
  for ii=0,nvline-1 do fline(ii) = array((nlines-1)*nvline+ii) 
  printf,unit,format=format,fline
endif else if nextra eq 1 then begin
  format = "(5x,'+ ',e12.4,'/')"
  printf,unit,format=format,array(id1-1)
endif else begin
  format = "(5x,'+ '," + strcompress(string(nextra-1),/remove_all)+$
    "(e12.4,','),e12.4,'/')"
  fline = fltarr(nextra)
  for ii=0,nextra-1 do fline(ii) = array(id1-nextra+ii)
  printf,unit,format=format,fline
endelse
return
end
;
pro wrdata, array,name,luout=luout
;
; Given a floating point array of 1 or 2 dimensions write fortran 
;   style data statements to luout using "name" as array name.
; On input, luout is already open and positioned for write.
;
if n_elements(luout) gt 0 then $
  unit = luout $
else $
  unit = -1
arraysize = size(array)
ndim = arraysize(0)
if ndim ne 1 and ndim ne 2 then begin
  print,'>>> wrdata: array must be 1 or 2 dimensions: ndim=',ndim
  return
endif
if ndim eq 1 then begin
  id1 = arraysize(1)
  itype = arraysize(2)
  if itype ne 4 then begin
    print,'>>> wrdata: array should be floating point: itype=',itype
    return 
  endif
  fname = '      DATA ' + name + ' /'
  wrdata_1d,array,id1,fname,luout
endif else begin
  id1 = arraysize(1)
  id2 = arraysize(2)
  for i=0,id2-1 do begin
    fname = '      DATA (' + name + "(i," + $
      strcompress(string(i+1),/remove_all) + "),i=1," + $
      strcompress(string(id1),/remove_all) + ") /"
    wrdata_1d,array(*,i),id1,fname,luout
  endfor
endelse
return
end