tsvel2

 TS Reads a scalar variable at a particular station from TSEPIC.CDF files.
 Returns data from all sigma layers, or intepolates to specified depths.

 USAGE:
   1.  [u,jd,z]=tsvel2(cdf,sta);  to return a matrix containing time,depth
                                 information for variable "var" at station
                                 sta for all sigma levels.

   2.  [u,jd]=tsvel2(cdf,sta,depths);
                                 to return a matrix containing time,depth
                                 information for variable "var" at station
                                 sta for depths specified in vector "depths".

       Example: [t,jd]=tsvel2('tsepic.cdf',2,[-2 -10 -24]);
              Returns a matrix containing temperature at Station 2 at
              2, 10 and 24 m below surface.  jd is the julian day vector.

   If the User specifies positive depths, they denote distance above bottom
   instead of distance below the free surface

0001 function [w,jd,z]=tsvel2(cdf,ista,zuser);
0002 % TS Reads a scalar variable at a particular station from TSEPIC.CDF files.
0003 % Returns data from all sigma layers, or intepolates to specified depths.
0004 %
0005 % USAGE:
0006 %   1.  [u,jd,z]=tsvel2(cdf,sta);  to return a matrix containing time,depth
0007 %                                 information for variable "var" at station
0008 %                                 sta for all sigma levels.
0009 %
0010 %   2.  [u,jd]=tsvel2(cdf,sta,depths);
0011 %                                 to return a matrix containing time,depth
0012 %                                 information for variable "var" at station
0013 %                                 sta for depths specified in vector "depths".
0014 %
0015 %       Example: [t,jd]=tsvel2('tsepic.cdf',2,[-2 -10 -24]);
0016 %              Returns a matrix containing temperature at Station 2 at
0017 %              2, 10 and 24 m below surface.  jd is the julian day vector.
0018 %
0019 %   If the User specifies positive depths, they denote distance above bottom
0020 %   instead of distance below the free surface
0021 %
0022 
0023 % use Chuck's tookit
0024 nc=netcdf(cdf);
0025 
0026 % valid station?
0027 nsta=length(nc{'stations'}(:));
0028 if(ista>nsta),
0029   disp('invalid station');
0030   return
0031 end
0032 
0033 % determine if the variable is in the center of the sigma layer or
0034 % at the top
0035 
0036 vloc='center';
0037 
0038 % get all the data at station ista
0039 uz=nc{'u'}(:,:,:,ista);
0040 uz=uz+sqrt(-1)*nc{'v'}(:,:,:,ista);
0041 
0042 [m,n]=size(uz);
0043 % find out the vertical coordinate (depth + elevation)
0044 
0045 sigma=nc{'sigma'}(:);
0046 
0047 sigma2=sigma(1:n-1)+0.5*diff(sigma);
0048 
0049 % add on the bottom, so that interpolations between the bottom grid
0050 % cell center and the bottom can be made
0051 sigma2(n)=-1;
0052 
0053 base_date=zeros(1,6);
0054 t=nc{'time'}(:);
0055 base_date(1:3)=nc.base_date(:);
0056 jd0=julian(base_date);
0057 jd=jd0+t/(3600*24);
0058 
0059 depth=nc{'depth'}(ista);
0060 elev=nc{'elev'}(:,:,:,ista);
0061 close(nc)
0062 tdepth=depth+elev;
0063 switch vloc
0064  case 'top'
0065    zlevt=tdepth*sigma';    % total depth including surface
0066  case 'center'
0067    zlevt=tdepth*sigma2';    % total depth including surface
0068 end
0069 %   want length(sigma) rows
0070 
0071 uz=uz.';
0072 zlevt=zlevt.';
0073   
0074 % Assign bottom half of bottom grid cell the same
0075 % value as at the grid cell center.
0076 uz(n,:)=uz(n-1,:);
0077 
0078 for k=1:length(zuser);
0079    % zlev(k) negative denotes distance below surface
0080    % zlev(k) positive denotes distance above bottom
0081      if(zuser(k)<0),
0082        zlev=zlevt;
0083      else
0084        zlev=zlevt-ones(n,1)*zlevt(n,:);
0085      end
0086      zind=zlev < zuser(k);
0087     % find the indices IK of the cells that are just above the required level
0088     dz=diff(zind);
0089     dz(n,:)=zeros(size(dz(n-1,:)));
0090     ik=find(dz==1);
0091 
0092 % find the indices ISUM where the requested depth is between two data values
0093     isum=find(sum(dz));
0094     u=zeros(m,1);
0095 % do linear interpolation
0096    u(isum)=uz(ik)+(uz(ik+1)-uz(ik)).*(zuser(k)-zlev(ik))./(zlev(ik+1)-zlev(ik));
0097 
0098 % find requested values that are above top data value but below surface
0099    if(zuser(k)<0),
0100     iabove=find(zlev(1,:) < zuser(k) & zuser(k) < elev(:)');
0101     itoohigh=find(zuser(k) > elev(:)');
0102    else
0103     iabove=find(zlev(1,:) < zuser(k) & zuser(k) < tdepth(:)');
0104     itoohigh=find(zuser(k) > tdepth(:)');
0105    end
0106 % assign top data value to requested values above it
0107    u(iabove)=uz(1,iabove);
0108 
0109 % find requested values that are above the free surface & mask 'em
0110    u(itoohigh)=uz(1,itoohigh)*nan;
0111 
0112    w(:,k)=u(:);
0113 
0114 end