scrum_islice

ISLICE  returns a vertical slice along i=iindex from a SCRUM file

 The variable must be 4D.

 USAGE: 
 >> [u,y,z]=islice(cdf,var,timestep,iindex,[jrange])
       u = the selected variable
       y = distance in *km* (assuming y units in netCDF file are in meters)
       z = depth in m
       iindex = I index along which slice is taken
       jrange = jmin and jmax indices along slice (optional).  If this
           argument is not supplied the default takes all the J indices
           except for the first and last, which are always "land" cells.


 see also JSLICE, KSLICE, ZSLICE, ZSLICEUV, KSLICEUV

0001 function [u,y,z] = scrum_islice(cdf,var,timestep,iindex,jrange)
0002 %ISLICE  returns a vertical slice along i=iindex from a SCRUM file
0003 %
0004 % The variable must be 4D.
0005 %
0006 % USAGE:
0007 % >> [u,y,z]=islice(cdf,var,timestep,iindex,[jrange])
0008 %       u = the selected variable
0009 %       y = distance in *km* (assuming y units in netCDF file are in meters)
0010 %       z = depth in m
0011 %       iindex = I index along which slice is taken
0012 %       jrange = jmin and jmax indices along slice (optional).  If this
0013 %           argument is not supplied the default takes all the J indices
0014 %           except for the first and last, which are always "land" cells.
0015 %
0016 %
0017 % see also JSLICE, KSLICE, ZSLICE, ZSLICEUV, KSLICEUV
0018 %
0019 
0020 
0021 if ( (nargin < 4) | ( nargin > 5) ) ,
0022   help scrum_islice; return
0023 end
0024 
0025 %
0026 % Suppress NetCDF warnings.
0027 ncmex('setopts',0);
0028 
0029 ncid=ncmex('open',cdf,'nowrite');
0030 if(ncid==-1),
0031   disp(['file ' cdf ' not found'])
0032   return
0033 end
0034 
0035 [name, ny]=ncmex('diminq',ncid,'eta_rho');
0036 [name, nz]=ncmex('diminq',ncid,'s_rho');
0037 if(exist('jrange')),
0038   jrange(1)=max(1,jrange(1));
0039   jrange(2)=min(ny-1,jrange(2));
0040   jstart=jrange(1)-1;
0041   jcount=jrange(2)-jrange(1)+1;
0042 else
0043   jstart=1
0044   jcount=ny-2;
0045 end
0046 
0047 %
0048 % ecom_islice just goes to (nz-1), don't really
0049 % know why.
0050 u = ncmex( 'varget', ncid, var, ...
0051             [(timestep-1) 0 jstart (iindex-1)],...
0052             [1 nz jcount 1] );
0053 u = squeeze(u);
0054 
0055 %
0056 % If the grid itself was asked for..
0057 if ( nargout > 1 )
0058     h = ncmex( 'varget', ncid, 'h', [jstart (iindex-1)], [jcount 1] );
0059     h = h';
0060 
0061     %
0062     % Find the land cells with mask_rho.
0063     [mask_rho_varid, status] = ncmex ( 'varid', ncid, 'mask_rho' );
0064     if ( mask_rho_varid ~= -1 )
0065         mask_rho = ncmex( 'varget', ncid, mask_rho_varid, ...
0066                        [jstart (iindex-1)], [jcount 1] );
0067         mask_rho = mask_rho';
0068         land = find ( mask_rho == 0 );
0069         h(land) = h(land) * NaN;
0070         u(land,:)=u(land,:)*NaN;
0071     end
0072     
0073 
0074     %
0075     % Get all the variables needed to compute z.
0076     % The equation is
0077     %
0078     % z = zeta * (1 + s) + hc*s + (h - hc)*C(s)
0079     [s_rho_dimid, status] = ncmex ( 'dimid', ncid, 's_rho' );
0080     if ( status == -1 )
0081         fprintf ( 2, 'Could not get s_rho dimid from %s.\n', cdf );
0082         ncmex ( 'close', ncid );
0083         return;
0084     end
0085 
0086     [dimname, s_rho_length, status] = ncmex ( 'diminq', ncid, s_rho_dimid );
0087     if ( status == -1 )
0088         fprintf ( 2, 'Could not get s_rho length from %s.\n', cdf );
0089         ncmex ( 'close', ncid );
0090         return;
0091     end
0092 
0093 
0094 
0095     %
0096     % w is defined at different locations than the others
0097     if ( strcmp(var,'w') )
0098         [sc, status] = ncmex ( 'varget', ncid, 'sc_w', [0], [-1] );
0099     else
0100         [sc, status] = ncmex ( 'varget', ncid, 'sc_r', [0], [-1] );
0101     end
0102 
0103 
0104     [zeta, status] = ncmex ( 'varget', ncid, 'zeta', [timestep jstart (iindex-1)], [1 jcount 1] );
0105     zeta = zeta';
0106 
0107     [hc, status] = ncmex ( 'varget1', ncid, 'hc', [0] );
0108     if ( status == -1 )
0109         fprintf ( 2, 'Could not get hc from %s.\n', cdf );
0110         ncmex ( 'close', ncid );
0111         return;
0112     end
0113 
0114 
0115     [Cs_r, status] = ncmex ( 'varget', ncid, 'Cs_r', [0], [-1] );
0116     if ( status == -1 )
0117         fprintf ( 'scrum_zslice:  could not get ''Cs_r'' in %s.', cdf );
0118         return;
0119     end
0120 
0121 
0122 
0123     z = zeta*(1+sc) + ones(size(h))*hc*sc + (h-hc)*Cs_r;
0124 
0125     pn = ncmex ( 'varget', ncid, 'pn', [jstart (iindex-1)], [jcount 1] );
0126     pn = pn';
0127 %    y=cumsum(pn)/1000;
0128 %    y=y*ones(1,nz);
0129     y = ncmex ( 'varget', ncid, 'y_rho', [jstart (iindex-1)], [jcount 1] );
0130     y = y';
0131     y=y*ones(1,nz);
0132     y = y/1000;
0133       ind=find(isnan(z));
0134     z(ind)=zeros(size(ind));
0135 
0136 end
0137 ncmex('close',ncid);
0138 
0139 u = u'; y = y'; z=z';
0140