Home > omviz > ecom_zslice.m

ecom_zslice

PURPOSE ^

ECOM_ZSLICE Returns a matrix containing a horizontal slice at a specified depth

SYNOPSIS ^

function [u,x,y,jd]=ecom_zslice(cdf,var,time,zuser)

DESCRIPTION ^

 ECOM_ZSLICE  Returns a matrix containing a horizontal slice at a specified depth
           at a given time step from an ECOMxx.CDF or POM.CDF file.
           Regions of grid that are shallower than requested value are
           returned as NaNs.

       USAGE: [u,x,y,jd]=zslice(cdf,var,time,zuser)

   where zuser is a depth in meters (e.g -10.)

 see also ZSLICEUV, ISLICE, JSLICE, KSLICE, KSLICEUV
 hint: use PSLICE or EXTCONTOUR to plot the results of ZSLICE

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [u,x,y,jd]=ecom_zslice(cdf,var,time,zuser)
0002 % ECOM_ZSLICE  Returns a matrix containing a horizontal slice at a specified depth
0003 %           at a given time step from an ECOMxx.CDF or POM.CDF file.
0004 %           Regions of grid that are shallower than requested value are
0005 %           returned as NaNs.
0006 %
0007 %       USAGE: [u,x,y,jd]=zslice(cdf,var,time,zuser)
0008 %
0009 %   where zuser is a depth in meters (e.g -10.)
0010 %
0011 % see also ZSLICEUV, ISLICE, JSLICE, KSLICE, KSLICEUV
0012 % hint: use PSLICE or EXTCONTOUR to plot the results of ZSLICE
0013 
0014 % Rich Signell  (rsignell@usgs.gov)
0015 % Changes:
0016 %  Dec 16, 1997:  Fixed problem with sigma levels not being correctly
0017 %                 centered.  Now T,S, RHO, U, V are computed in the
0018 %                 center of the vertical grids, and W, Q, and L are
0019 %                 computed at the tops of the vertical grids.  Also,
0020 %                 now if the user selects a zlevel that is below the
0021 %                 sea surface but above the top model data point, the
0022 %                 value at the top data point is returned. If the user
0023 %                 selects a value that is above the bottom but below
0024 %                 the bottom-most grid level, the bottom value is used.
0025 %                 An alternative, would be to use a log-layer formulation
0026 %                 for velocity components.
0027 
0028 if (nargin<2 | nargin>4),
0029   help zslice; return
0030 end
0031 mexcdf('setopts',0);
0032 ncid=mexcdf('open',cdf,'nowrite');
0033 if(ncid==-1),
0034   disp(['file ' cdf ' not found'])
0035   return
0036 end
0037 
0038 [name, nx]=mexcdf('diminq',ncid,'xpos');
0039 [name, ny]=mexcdf('diminq',ncid,'ypos');
0040 [name, nz]=mexcdf('diminq',ncid,'zpos');
0041 
0042 
0043 %
0044 % Acquire the grid.
0045 %
0046 % If "lon" and "lat" are present, grab them.
0047 % Otherwise, get "x" and "y".
0048 [lon_varid, rcode] = mexcdf('VARID', ncid, 'lon');
0049 [lat_varid, rcode] = mexcdf('VARID', ncid, 'lat');
0050 if ( (lon_varid >= 0) | (lat_varid >= 0) )
0051     x=ncmex('varget',ncid,'lon',[0 0],[-1 -1]);
0052     y=ncmex('varget',ncid,'lat',[0 0],[-1 -1]);
0053 else
0054     x=ncmex('varget',ncid,'x',[0 0],[-1 -1]);
0055     y=ncmex('varget',ncid,'y',[0 0],[-1 -1]);
0056 end
0057 
0058 
0059 %if(nargout>2),
0060 %  x=mexcdf('varget',ncid,'x',[0 0],[ny nx]);
0061 %  y=mexcdf('varget',ncid,'y',[0 0],[ny nx]);
0062 %end
0063 
0064 
0065 
0066 
0067 depth=mexcdf('varget',ncid,'depth',[0 0],[ny nx]);
0068 elev=mexcdf('varget',ncid,'elev',[(time-1) 0 0],[1 ny nx],1);
0069 dind=find(depth==-99999);
0070 igood=find(depth~=-99999);
0071 
0072 % use total depth (depth+elev) for sigma model
0073 tdepth=depth+elev;
0074 
0075 sigma=mexcdf('varget',ncid,'sigma',0,nz);
0076 n=length(sigma);
0077 
0078 % vertical location of grid cell centers (s,rho,u,v,etc, but not w)
0079 sigma2=sigma(1:n-1)+0.5*diff(sigma);
0080 
0081 % add on the bottom, so that interpolations between the bottom grid
0082 % cell center and the bottom can be made
0083 sigma2(n)=-1;
0084 
0085 % grab 3d slab (all values) at a given time step
0086 uz=mexcdf('varget',ncid,var,[(time-1) 0 0 0],[1 nz ny nx],1); %3d slab
0087 
0088 % get time stuff
0089 if(nargout>3)
0090 [nam,nt]=mexcdf('diminq',ncid,'time');
0091 base_date=[0 0 0 0 0 0];
0092 base_date(1:3)=mexcdf('attget',ncid,'global','base_date');
0093 t=mexcdf('varget',ncid,'time',time-1,1);
0094 jd0=julian(base_date);
0095 jd=jd0+t;
0096 end
0097 
0098 % Close file
0099 mexcdf('close',ncid);
0100 
0101 % Reshape the 3D array so each grid cell's vertical profile is a column.
0102 uz=reshape(uz,nx*ny,n).';
0103 
0104 % Assign bottom half of bottom grid cell the same
0105 % value as at the grid cell center.
0106 uz(n,:)=uz(n-1,:);   
0107 
0108 % Form an array the same shape as uz that contains the depths of the
0109 % grid cell centers.  (multiply sigma by total depth (depth +elev))
0110 switch var
0111 case {'w','kh','km'}
0112    zlev=sigma*(tdepth(:)');
0113 otherwise
0114    zlev=sigma2*(tdepth(:)');
0115 end
0116 % find all the values of zlev that are less than the requested depth
0117 % zind consists of zeros and ones.  One if the cell depth is less
0118 % than (below) the requested z level, zero if cell level is greater.
0119 zind=zlev < zuser;
0120 
0121 % find the indices IK of the cells that are just above the required level
0122 dz=diff(zind);
0123 dz(n,:)=zeros(size(dz(n-1,:)));
0124 ik=find(dz==1);
0125 
0126 % find the indices ISUM where the requested depth is between two data values
0127 isum=find(sum(dz));
0128 u=zeros(nx*ny,1);
0129 
0130 % do linear interpolation
0131 u(isum)=uz(ik)+(uz(ik+1)-uz(ik)).*(zuser-zlev(ik))./(zlev(ik+1)-zlev(ik));
0132 
0133 % find requested values that are above top data value but below surface
0134 iabove=find(zlev(1,:) < zuser & zuser < elev(:)');
0135 
0136 % assign top data value to requested values above it
0137 u(iabove)=uz(1,iabove);
0138 
0139 % find requested values that are above the free surface & mask 'em
0140 itoohigh=find(zuser > elev(:)');
0141 u(itoohigh)=uz(1,itoohigh)*nan;
0142 
0143 % put back in matrix form
0144 u=reshape(u,nx,ny);
0145 
0146 % mask values where the total depth is less than the requested depth
0147 tind=find(-tdepth > zuser);
0148 u(tind)=u(tind)*NaN;
0149 
0150 % mask land values
0151 u(dind)=u(dind)*NaN;
0152

Generated on Thu 01-Dec-2005 08:50:35 by m2html © 2003