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 ZSLICE0001 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