Home > omviz > ecom_zsliceuv.m

ecom_zsliceuv

PURPOSE ^

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

SYNOPSIS ^

function [w,x,y,jd]=zsliceuv(cdf,time,zuser)

DESCRIPTION ^

  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 ZSLICE

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 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);

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