;-----------------------------------------------------------------------
function mtime_to_mins,mtime ; model time to integer minutes
return,mtime(0)*24*60+mtime(1)*60+mtime(2)
end
;-----------------------------------------------------------------------
function mtime_to_ddays,mtime ; model time to decimal days
return,float(mtime(0))+float(mtime(1))/24.+float(mtime(2))/(24.*60.)
end
;-----------------------------------------------------------------------
pro dayofyear_to_monthday, dayofyear,imonth,iday
;
; Given day of year (1-365), return month and day
;
iday = 0
imonth = 0
if (dayofyear eq 0) then dayofyear = 1 ; to allow for model day of 0
if (dayofyear lt 1 or dayofyear gt 365) then begin
  print,'>>> bad dayofyear=',dayofyear
  return
endif
;                J, F, M, A, M, J, J, A, S, O, N, D
dayspermonth = [31,28,31,30,31,30,31,31,30,31,30,31]
nmonths = n_elements(dayspermonth)
;
id = 0
for m=0,nmonths-1 do begin
  for i=1,dayspermonth[m] do begin
    id = id+1
    if (id eq dayofyear) then begin
      imonth = m+1
      iday = i
      break  
    endif
  endfor
  if (imonth ne 0) then break ; day was found above
endfor
;print,'dayofyear_to_monthday: dayofyear=',dayofyear,' imonth=',imonth,' iday=',iday
end
;-----------------------------------------------------------------------
function timeaxis,decyears,years,mtimes,deltamin,xtitle,nxticks
;
; As called from pltutlat:
;     time = timeaxis(years[0:ntime_plot-1],mtime_plot[*,0:ntime_plot-1],$
;       deltamin,xtitle,nxticks)
;
; Get julian day of first time:
;
dayofyear_to_monthday,mtimes(0,0),imonth0,iday0
julian0 = julday(imonth0,iday0,years[0],mtimes(1,0),mtimes(2,0),0)
;
; Generate x-axis coords in julian days (step size in minutes):
; (this array is returned as the function value)
;
ntimes = n_elements(years)
time = timegen(ntimes,start=julian0,step_size=deltamin,units="Minutes")
;
; Add 0.5 to times since the julian day starts at 12 pm noon
; (s.a., caldat call below)
;
time = time+0.5
;
; Construct x-axis label giving beginning,ending model times,
; delta time in miniutes, and total time in days:
;
ddays0 = mtime_to_ddays(mtimes[*,0])
ddays1 = mtime_to_ddays(mtimes[*,ntimes-1])
;total_ddays = ddays1-ddays0 ; total number of decimal days

; Decyears is passed from pltutvert:
;   for i=itime0_plot,itime1_plot do begin
;     decyears[i] = double(years[i])+double(days[i]-1)/365.+ut[i]/(365.*24.)
;   endfor

total_ddays = (decyears[ntimes-1]-decyears[0])*365.

;print,'timeaxis: ntimes=',ntimes,' decyears[0]=',decyears[0],$
;  ' decyears[ntimes-1]=',decyears[ntimes-1],' total_ddays=',total_ddays 

nxticks = 0 ; let idl select
;if (total_ddays ge 364) then nxticks = 7 ; does not appear to work
xtitle = string(format=$
  "('Model time (yyyy:ddd:hh:mm) ',i4,':',i3.3,':',i2.2,':',i2.2,' to ')",$
  years[0],mtimes[*,0])
xtitle = xtitle + string(format=$
  "(i4,':',i3.3,':',i2.2,':',i2.2,' by ',i6,' minutes ')",$
  years[ntimes-1],mtimes[*,ntimes-1],deltamin)
xtitle = xtitle + string(format="('(',f7.2,' days)')",total_ddays)
return,time
end