ts2

 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]=ts(cdf,var,sta);  to return a matrix containing time,depth
                                 information for variable "var" at station
                                 sta for all sigma levels.

   2.  [u,jd]=ts(cdf,var,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]=ts('tsepic.cdf','temp',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 depths are positive, assume that they are heights above bottom

   z are the *water depths*, elevation not included

0001 function [w,jd,z]=ts2(cdf,var,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]=ts(cdf,var,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]=ts(cdf,var,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]=ts('tsepic.cdf','temp',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 %   if depths are positive, assume that they are heights above bottom
0019 %
0020 %   z are the *water depths*, elevation not included
0021 
0022 % use Chuck's tookit
0023 nc=netcdf(cdf);
0024 
0025 % valid station?
0026 nsta=length(nc{'stations'}(:));
0027 if(ista>nsta),
0028   disp('invalid station');
0029   return
0030 end
0031 
0032 % determine if the variable is in the center of the sigma layer or
0033 % at the top
0034 
0035 switch var
0036  case {'u','v','salt','temp','conc','am'}
0037    vloc='center';
0038  case {'w','kh','km'}
0039    vloc='top';
0040  otherwise
0041    disp(['Error: Unknown variable ' var ': please modify TS.M!']);
0042    return
0043 end
0044 
0045 % get all the data at station ista
0046 uz=nc{var}(:,:,:,ista);
0047 [m,n]=size(uz);
0048 % find out the vertical coordinate (depth + elevation)
0049 
0050 sigma=nc{'sigma'}(:);
0051 
0052 sigma2=sigma(1:n-1)+0.5*diff(sigma);
0053 
0054 % add on the bottom, so that interpolations between the bottom grid
0055 % cell center and the bottom can be made
0056 sigma2(n)=-1;
0057 
0058 base_date=zeros(1,6);
0059 t=nc{'time'}(:);
0060 base_date(1:3)=nc.base_date(:);
0061 jd0=julian(base_date);
0062 jd=jd0+t/(3600*24);
0063 
0064 depth=nc{'depth'}(ista);
0065 elev=nc{'elev'}(:,:,:,ista);
0066 close(nc)
0067 tdepth=depth+elev;
0068 switch vloc
0069  case 'top'
0070    zlevt=tdepth*sigma';    % total depth including surface
0071  case 'center'
0072    zlevt=tdepth*sigma2';    % total depth including surface
0073 end
0074 %   want length(sigma) rows
0075 
0076 uz=uz.';
0077 zlevt=zlevt.';
0078   
0079 % Assign bottom half of bottom grid cell the same
0080 % value as at the grid cell center.
0081 uz(n,:)=uz(n-1,:);
0082 
0083 for k=1:length(zuser);
0084    % zlev(k) negative denotes distance below surface
0085    % zlev(k) positive denotes distance above bottom
0086      if(zuser(k)<0),
0087        zlev=zlevt;
0088      else
0089        zlev=zlevt-ones(n,1)*zlevt(n,:);
0090      end
0091      zind=zlev < zuser(k);
0092     % find the indices IK of the cells that are just above the required level
0093     dz=diff(zind);
0094     dz(n,:)=zeros(size(dz(n-1,:)));
0095     ik=find(dz==1);
0096 
0097 % find the indices ISUM where the requested depth is between two data values
0098     isum=find(sum(dz));
0099     u=zeros(m,1);
0100 % do linear interpolation
0101    u(isum)=uz(ik)+(uz(ik+1)-uz(ik)).*(zuser(k)-zlev(ik))./(zlev(ik+1)-zlev(ik));
0102 
0103 % find requested values that are above top data value but below surface
0104    if(zuser(k)<0),
0105     iabove=find(zlev(1,:) < zuser(k) & zuser(k) < elev(:)');
0106     itoohigh=find(zuser(k) > elev(:)');
0107    else
0108     iabove=find(zlev(1,:) < zuser(k) & zuser(k) < tdepth(:)');
0109     itoohigh=find(zuser(k) > tdepth(:)');
0110    end
0111 % assign top data value to requested values above it
0112    u(iabove)=uz(1,iabove);
0113 
0114 % find requested values that are above the free surface & mask 'em
0115    u(itoohigh)=uz(1,itoohigh)*nan;
0116 
0117    w(:,k)=u(:);
0118 
0119 end