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