function [t,t0,lon,lat] = thredds_coverage(url)
% USAGE: thredds_coverage(url)
%        thredds_coverage(ncfile)
%
% simple script to return metadata on geospatial and time coverage
% for a thredds catalog.xml

try
  lon=nc_varget(url,'lon');
  lat=nc_varget(url,'lat');
catch
  lon=nc_varget(url,'lon_rho');
  lat=nc_varget(url,'lat_rho');
end

try
  t=nc_varget(url,'time');
  t=nc_varget(url,'time');
  tunits = nc_attget(url,'time','units');
  t0=datenum(parsetnc(tunits));
catch
  t=nc_varget(url,'ocean_time');
  tunits = nc_attget(url,'ocean_time','units');
  t0=datenum(parsetnc(tunits));
end
switch tunits(1:3)
  case 'day'
    fac = 1;
  case 'sec'
    fac = 1/86400;
  otherwise
    error(['Did not recognize time units ' tunits])
end
t = fac*t;

disp(' ')
disp(['          <timeCoverage>'])
disp(['             <start>' datestr(t0+t(1),'yyyy-mm-dd HH:MM:SS') ' UTC</start>'])
disp(['             <end>' datestr(t0+t(end),'yyyy-mm-dd HH:MM:SS') ' UTC</end>'])
disp(['             <resolution>' sprintf('%6.4f',diff(t(1:2))) ' days</resolution>'])
disp('          </timeCoverage>')

disp(' ')
disp(['             <northsouth>'])
disp(['                <start>' sprintf('%5.2f',min(lat(:))) '</start>'])
disp(['                <size>' sprintf('%5.2f',max(lat(:))-min(lat(:))) '</size>'])
disp(['                <units>degrees_north</units>'])
disp('             </northsouth>')


disp(' ')
disp(['             <eastwest>'])
disp(['               <start>' sprintf('%5.2f',min(lon(:))) '</start>'])
disp(['               <size>' sprintf('%5.2f',max(lon(:))-min(lon(:))) '</size>'])
disp(['               <units>degrees_east</units>'])
disp('             </eastwest>')

