ecomts

 ecomts extracts a velocity time series from a 4D netCDF ECOM output file.

   USAGE:  [u,jd,depth]=ecomts(cdf,var,i,j,[zi],[tind]);

 cdf = cdf file e.g. 'ecomsi.cdf'
 (i,j) = station location in grid coordinates
  var = variable (e.g. 'salt')
  zi = z levels from which to extract time series 
    If zi=nan or missing, the default is to extract data at all sigma levels
 iind= time indices 

 w = vector time series of "var"
 jd = time vector  (digital Julian days)
 depth = depth(s) of velocity time series (meters)

 examples:  To extract salinity at all levels and all time steps
            for cell(i=30,j=47) from the file 'ecomsi.cdf': 

              [w,jd,depths]=ecomts('ecomsi.cdf','salt',30,47);

            To extract temperature at 5 and 25 m for time steps 1:10,
            for cell(i=30,j=47) from the file 'ecomsi.cdf': 

              [w,jd,depths]=ecomts('ecomsi.cdf','temp',30,47,[-5 -25],1:10);

            To extract temperature at 5 and 25 m for all time steps
            for cell(i=30,j=47) from the file 'ecomsi.cdf': 

              [w,jd,depths]=ecomts('ecomsi.cdf','temp',30,47,nan,1:10);

0001 function [w,jd,depth]=ecomts(cdf,var,i,j,zi,iind)
0002 % ecomts extracts a velocity time series from a 4D netCDF ECOM output file.
0003 %
0004 %   USAGE:  [u,jd,depth]=ecomts(cdf,var,i,j,[zi],[tind]);
0005 %
0006 % cdf = cdf file e.g. 'ecomsi.cdf'
0007 % (i,j) = station location in grid coordinates
0008 %  var = variable (e.g. 'salt')
0009 %  zi = z levels from which to extract time series
0010 %    If zi=nan or missing, the default is to extract data at all sigma levels
0011 % iind= time indices
0012 %
0013 % w = vector time series of "var"
0014 % jd = time vector  (digital Julian days)
0015 % depth = depth(s) of velocity time series (meters)
0016 %
0017 % examples:  To extract salinity at all levels and all time steps
0018 %            for cell(i=30,j=47) from the file 'ecomsi.cdf':
0019 %
0020 %              [w,jd,depths]=ecomts('ecomsi.cdf','salt',30,47);
0021 %
0022 %            To extract temperature at 5 and 25 m for time steps 1:10,
0023 %            for cell(i=30,j=47) from the file 'ecomsi.cdf':
0024 %
0025 %              [w,jd,depths]=ecomts('ecomsi.cdf','temp',30,47,[-5 -25],1:10);
0026 %
0027 %            To extract temperature at 5 and 25 m for all time steps
0028 %            for cell(i=30,j=47) from the file 'ecomsi.cdf':
0029 %
0030 %              [w,jd,depths]=ecomts('ecomsi.cdf','temp',30,47,nan,1:10);
0031 
0032 % Rich Signell (rsignell@usgs.gov)
0033 
0034 % suppress netcdf warnings
0035 mexcdf('setopts', 0);
0036 
0037 ncid=mexcdf('open',cdf,'nowrite');
0038 [nam,nz]=mexcdf('diminq',ncid,'zpos');
0039 [nam,nt]=mexcdf('diminq',ncid,'time');
0040 if(exist('iind')==1),
0041   istart=min(iind)-1;
0042   icount=min(max(iind),nt)-istart;
0043 else
0044   istart=0;
0045   icount=nt;
0046 end
0047 sigma=mexcdf('varget',ncid,'sigma',0,nz);
0048 sigma=.5*(sigma(1:(nz-1))+sigma(2:nz));
0049 depth=mexcdf('varget',ncid,'depth',[j-1 i-1],[1 1]);
0050 
0051 base_date=[0 0 0 0 0 0];
0052 base_date(1:3)=mexcdf('attget',ncid,'global','base_date');
0053 t=mexcdf('varget',ncid,'time',istart,icount);
0054 jd0=julian(base_date);
0055 jd=jd0+t;
0056 
0057 w=mexcdf('varget',ncid,var,[istart 0 j-1 i-1],[icount -1 1 1],1);
0058 w=w.';
0059 mexcdf('close',ncid);
0060 
0061 depth=sigma*depth;
0062 
0063 if(exist('zi')==1),
0064   if(~isnan(zi(1))),
0065   m=length(zi);
0066   if(min(zi)<min(depth)), disp('requested level below data!'),return,end
0067   if(max(zi)>max(depth)), disp('requested level above data!'),return,end
0068   for k=1:m,
0069     lev2=max(find(depth>zi(k)));
0070     lev1=lev2+1;
0071     frac=(zi(k)-depth(lev1))/(depth(lev2)-depth(lev1));
0072     wmod(:,k)=w(:,lev1)+frac*(w(:,lev2)-w(:,lev1));
0073   end
0074   w=wmod;
0075   depth=zi;
0076   end
0077 end