Skip to content

GitLab

fvcom-toolbox
fvcom-toolbox
Commit 378f9be5 authored by Pierre Cazenave's avatar Pierre Cazenave
Browse files
Options

Merge branch 'dev'

parents 59350f4e 9687b1fc
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 1200 additions and 84 deletions
examples/example_fvcom_inputs.m 0 → 100644
View file @ 378f9be5
This diff is collapsed.
fvcom_prepro/add_obc_nodes_list.m
View file @ 378f9be5
function [Mobj] = add_obc_nodes_list(Mobj,Nlist,ObcName,ObcType,plotFig)
function [Mobj] = add_obc_nodes_list(Mobj,Nlist,ObcName,ObcType,plotFig)
% Add a set of obc nodes comprising a single obc boundary to Mesh structure
% Add a set of obc nodes comprising a single obc boundary to Mesh structure
% Using a list of nodes
%
% [Mobj] = add_obc_nodes_list(Mobj,Nlist,ObcName,ObcType)
......@@ -21,7 +21,7 @@ function [Mobj] = add_obc_nodes_list(Mobj,Nlist,ObcName,ObcType,plotFig)
% EXAMPLE USAGE
% Mobj = add_obc_nodes_list(Mobj,Nlist,'OpenOcean')
%
% Author(s):
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
% Pierre Cazenave (Plymouth Marine Laboratory)
% Karen Amoudry (National Oceanography Centre, Liverpool)
......@@ -31,13 +31,14 @@ function [Mobj] = add_obc_nodes_list(Mobj,Nlist,ObcName,ObcType,plotFig)
% 2013-01-02 KJA bug fix: amended usage of 'unique' in line 53 to
% prevent it from sorting the values it returns. Amended by Pierre to
% support pre-2012 versions of MATLAB whilst giving the same result.
%
% 2015-02-23 Output number of nodes if the verbose flag is set.
%
%==========================================================================
subname = 'add_obc_nodes';
subname = 'add_obc_nodes_list';
global ftbverbose
if(ftbverbose)
fprintf('\n')
fprintf(['begin : ' subname '\n'])
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
% Do we want a figure showing how we're getting along?
......@@ -46,7 +47,7 @@ if nargin == 4
end
%--------------------------------------------------------------------------
% Get a unique list and make sure they are in the range of node numbers
% Get a unique list and make sure they are in the range of node numbers
%--------------------------------------------------------------------------
% Make this works in versions of MATLAB older than 2012a (newer versions
% can just use unique(A, 'stable'), but checking versions is a pain).
......@@ -54,12 +55,12 @@ end
Nlist = Nlist(sort(Nidx));
if max(Nlist) > Mobj.nVerts
fprintf('your open boundary node number exceed the total number of nodes in the domain\n');
fprintf('Your open boundary node number exceed the total number of nodes in the domain\n');
error('stopping...')
end
%--------------------------------------------------------------------------
% Plot the mesh
% Plot the mesh
%--------------------------------------------------------------------------
if plotFig == 1
if strcmpi(Mobj.nativeCoords(1:3), 'car')
......@@ -71,25 +72,28 @@ if plotFig == 1
end
figure
patch('Vertices',[x,y],'Faces',Mobj.tri,...
'Cdata',Mobj.h,'edgecolor','k','facecolor','interp');
hold on;
patch('Vertices', [x, y] , 'Faces', Mobj.tri, ...
'Cdata', Mobj.h, 'edgecolor', 'k', 'facecolor', 'interp')
hold on
whos Nlist
plot(x(Nlist),y(Nlist),'ro');
axis('equal','tight')
plot(x(Nlist), y(Nlist), 'ro');
axis('equal', 'tight')
title('open boundary nodes');
end
% add to mesh object
npts = numel(Nlist);
Mobj.nObs = Mobj.nObs + 1;
Mobj.nObcNodes(Mobj.nObs) = npts;
Mobj.nObcNodes(Mobj.nObs) = npts;
Mobj.obc_nodes(Mobj.nObs,1:npts) = Nlist;
Mobj.obc_name{Mobj.nObs} = ObcName;
Mobj.obc_type(Mobj.nObs) = ObcType;
if ftbverbose
fprintf('found %d open boundary nodes', npts)
end
if(ftbverbose)
fprintf(['end : ' subname '\n'])
if ftbverbose
fprintf('\nend : %s\n', subname)
end
fvcom_prepro/add_stations_list.m
View file @ 378f9be5
......@@ -10,7 +10,12 @@ function [Mobj] = add_stations_list(Mobj,Positions,Names,Dist)
% grid node to those supplied will be used in the output file.
%
% INPUT
% Mobj = Matlab mesh object
% Mobj = Matlab mesh object with the following fields:
% - nativeCoords = native grid coordinate type.
% - x, y = grid node positions
% - xc, yc = grid element positions
% - h = water depth at the nodes
% - tri = grid triangulation table
% Positions = 2xn array of the XY positions of the stations
% Names = Cell array of the names of the stations defined in Positions
% Dist = Maximum distance from a station for a node to be included
......@@ -21,8 +26,14 @@ function [Mobj] = add_stations_list(Mobj,Positions,Names,Dist)
% ensure Dist is in those units.
%
% OUTPUT:
% Mobj = Matlab mesh object with an additional cell array containing id,
% x, y, nodelist, depth and station name.
% Mobj = Matlab mesh object with an additional cell array (stations)
% containing:
% {id, x, y, nodelist, depth, 'station name', elem}
% where id is the station ID (from 1 number of valid stations, x and y
% are the coordinates of the station, nodelist is the indices of the
% model grid nodes, depth is the depth at the grid node, station name is
% a string of the name from the input file and eleme is the grid element
% ID closest to each station.
%
% EXAMPLE USAGE
% Mobj = add_stations_list(Mobj, [-5.54, 50.103; -3.0865, 58.441], ...
......@@ -31,17 +42,18 @@ function [Mobj] = add_stations_list(Mobj,Positions,Names,Dist)
% Author(s):
% Pierre Cazenave (Plymouth Marine Laboratory)
%
%
% Revision history
% 2012-11-30 First version.
% 2015-01-14 Add support for exporting the element ID closest to the
% station of interest (as well as the node ID which was already there).
% Fix some formatting issues.
%
%==========================================================================
subname = 'add_stations_list';
global ftbverbose
if(ftbverbose)
fprintf('\n')
fprintf(['begin : ' subname '\n'])
end;
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
%--------------------------------------------------------------------------
% Check the inputs
......@@ -74,19 +86,20 @@ else
end
inc = 1;
% out = cell(nPos, 1);
out = cell(1); % don't preallocate as we don't know how many we'll have
for s=1:nPos
for s = 1:nPos
[node, dist] = find_nearest_pt(Positions(s, cols(1)), Positions(s, cols(2)), Mobj);
[~, elem] = min(abs(sqrt((Mobj.xc - Positions(s, cols(1))).^2 + Mobj.yc - Positions(s, cols(2))).^2));
if dist >= Dist
% Skip out for this station
if(ftbverbose)
if ftbverbose
fprintf('Skipping station %s (%g, %g). Nodal distance from station position = %f\n', Names{s}, Positions(s, 1), Positions(s, 2), dist)
end
continue
end
out{inc} = {inc, Positions(s, cols(1)), Positions(s, cols(2)), node, Mobj.h(node), Names{s}};
out{inc} = {inc, Positions(s, cols(1)), Positions(s, cols(2)), node, Mobj.h(node), Names{s}, elem};
inc = inc + 1;
end
......@@ -98,3 +111,7 @@ else
fprintf('No stations found within the model domain.\n')
end
end
if ftbverbose
fprintf('\nend : %s\n', subname)
end
\ No newline at end of file
fvcom_prepro/find_nearest_pt.m
View file @ 378f9be5
function [Point,Distance] = find_nearest_pt(xloc,yloc,Mobj)
% Find nearest point in Mesh structure to (x,y)
%
% function [Point,Distance] = find_nearest_pt(xloc,yloc,Mobj)
......@@ -8,11 +7,14 @@ function [Point,Distance] = find_nearest_pt(xloc,yloc,Mobj)
% Find nearest point to (xloc,yloc) in the domain of Mobj
% using native coordinates of Mobj
%
% INPUT:
% INPUT:
% xloc = x location of point (in native Mobj coordinates)
% yloc = y location of point (in native Mobj coordinates)
% Mobj = Mesh object
%
% Mobj = Mesh object with the following fields:
% - nativeCoords = grid type (cartesian or spherical)
% - x, y and/or lon, lat = coordinates (dependent on
% nativeCoords).
%
% OUTPUT:
% Point = index of nearest vertex in the mesh
% Distance = Distance from x,y to Point in Mobj native coordinates
......@@ -22,40 +24,47 @@ function [Point,Distance] = find_nearest_pt(xloc,yloc,Mobj)
%
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history
%
% 2015-01-14 Tidy up the code a bit and add extra information to the
% help.
%
%==============================================================================
% global ftbverbose
subname = 'find_nearest_pt';
%fprintf('\n')
%fprintf(['begin : ' subname '\n'])
% if ftbverbose
% fprintf('\nbegin : %s\n', subname)
% end
%------------------------------------------------------------------------------
% Parse input arguments
%------------------------------------------------------------------------------
if(exist('xloc')*exist('yloc')*exist('Mobj') == 0)
error('arguments to find_nearest_pt are missing')
end;
if ~exist('xloc', 'var') || ~exist('yloc', 'var') || ~exist('Mobj', 'var')
error('arguments to %s are missing', subname)
end
%------------------------------------------------------------------------------
% Set native coordinates
%------------------------------------------------------------------------------
if(lower(Mobj.nativeCoords(1:3)) == 'car')
if strcmpi(Mobj.nativeCoords, 'cartesian')
x = Mobj.x;
y = Mobj.y;
else
elseif strcmpi(Mobj.nativeCoords, 'spherical')
x = Mobj.lon;
y = Mobj.lat;
end;
else
error('Unrecognised coordinate type.')
end
%------------------------------------------------------------------------------
% Find the nearest point
%------------------------------------------------------------------------------
radvec = sqrt( (xloc-x).^2 + (yloc-y).^2);
[Distance,Point] = min(radvec);
[Distance, Point] = min(sqrt((xloc - x).^2 + (yloc - y).^2));
%fprintf(['end : ' subname '\n'])
% if ftbverbose
% fprintf(['end : ' subname '\n'])
% end
fvcom_prepro/get_NCEP_forcing.m
View file @ 378f9be5
......@@ -579,7 +579,17 @@ for aa = 1:length(fields)
% [~,length3] = netcdf.inqDim(ncid,dimids(3))
start=[min(index_lon)-1,min(index_lat)-1,min(data_time_idx)-1];
count=[length(index_lon),length(index_lat),length(data_time_idx)];
data1.(fields{aa}).(fields{aa}).(fields{aa}) = netcdf.getVar(ncid,varid,start,count,'double');
% The air data was failing with a three long start and count
% array, so try first without (to retain original behaviour for
% other potentially unaffected variables) but fall back to
% getting only the first level (start = 0, count = 1).
try
data1.(fields{aa}).(fields{aa}).(fields{aa}) = netcdf.getVar(ncid,varid,start,count,'double');
catch
start=[min(index_lon)-1,min(index_lat)-1,0,min(data_time_idx)-1];
count=[length(index_lon),length(index_lat),1,length(data_time_idx)];
data1.(fields{aa}).(fields{aa}).(fields{aa}) = netcdf.getVar(ncid,varid,start,count,'double');
end
else
eval(['data1.(fields{aa}) = loaddap(''', ncep.(fields{aa}),'?',...
......
fvcom_prepro/read_fvcom_obc.m 0 → 100644
View file @ 378f9be5
function Mobj = read_fvcom_obc(Mobj, obcfile)
% Read fvcom open boundary file
%
% Mobj = function read_fvcom_obc(obcfile)
%
% DESCRIPTION:
% Read FVCOM open boundary node specification file.
%
% INPUT:
% Mobj : existing mesh object into which obc data is added
% obcfile : fvcom open boundary file
%
% OUTPUT:
% Mobj.read_obc_nodes : open boundary node cell array (length = number
% of open boundaries).
%
% EXAMPLE USAGE
% Mobj = read_fvcom_obc('tst_obc.dat')
%
% Author(s):
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history
% 2015-01-12 Create function.
%
%==========================================================================
subname = 'read_fvcom_obc';
global ftbverbose
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
have_strings = false;
%--------------------------------------------------------------------------
% read in the FVCOM open boundary nodes
%--------------------------------------------------------------------------
fid = fopen(obcfile,'r');
assert(fid >= 0, 'file: %s does not exist.', obcfile)
reading = true;
while reading
lin = fgetl(fid);
if lin ~= -1 % EOF is -1
lin = strtrim(lin);
switch lin(1:2)
case 'OB'
tmp = regexp(lin, ' = ', 'split');
nObcNodes = str2double(tmp(end));
clear tmp
read_obc_nodes = cell(1, 1);
end
else
reading = false;
end
end
fclose(fid);
if nObcNodes > 0
have_strings = true;
fprintf('nObcNodes = %d\n', nObcNodes)
end
% We have to assume we have only a single open boundary as the _obc.dat
% file has no way of indicating how many open boundaries there are.
fid = fopen(obcfile,'r');
S = textscan(fid, '%f%f%f%*s%*s%[^\n\r]', 'HeaderLines', 1, ...
'Delimiter', {'\t',' '}, 'MultipleDelimsAsOne', true);
fclose(fid);
read_obc_nodes{1} = S{:, 2};
if have_strings
Mobj.have_strings = have_strings;
Mobj.read_obc_nodes = read_obc_nodes;
end
if ftbverbose
fprintf('\nend : %s\n', subname)
end
\ No newline at end of file
fvcom_prepro/read_sms_mesh.m
View file @ 378f9be5
......@@ -44,6 +44,7 @@ function [Mobj] = read_sms_mesh(varargin)
% 2013-12-11 Closed the sms_2dm file using fclose (ROM).
% 2014-04-10 Fix bugs when not using bathymetry (i.e. only reading the
% grid data in).
% 2015-03-19 Add spherical coordinates on element centres.
%
%==============================================================================
......@@ -348,6 +349,10 @@ Mobj.tri = tri;
% Make a depth array for the element centres.
Mobj.hc = nodes2elems(h, Mobj);
% Add element spherical coordinates too.
Mobj.lonc = nodes2elems(lon, Mobj);
Mobj.latc = nodes2elems(lat, Mobj);
if ftbverbose
fprintf('end : %s\n', subname)
end
fvcom_prepro/write_FVCOM_bath.m
View file @ 378f9be5
......@@ -27,7 +27,7 @@ subname = 'write_FVCOM_bath';
global ftbverbose
if(ftbverbose)
fprintf('\n'); fprintf(['begin : ' subname '\n']);
end;
end
%------------------------------------------------------------------------------
% Parse input arguments
......@@ -45,21 +45,21 @@ if(lower(Mobj.nativeCoords(1:3)) == 'car')
else
x = Mobj.lon;
y = Mobj.lat;
end;
end
if(Mobj.have_bath)
if(ftbverbose); fprintf('writing FVCOM bathymetry file %s\n',filename); end;
fid = fopen(filename,'w');
fprintf(fid,'Node Number = %d\n',Mobj.nVerts);
for i=1:Mobj.nVerts
fprintf(fid,'%f %f %f\n',x(i),y(i),Mobj.h(i));
end;
end
fclose(fid);
else
error('can''t write bathymetry to file, mesh object has no bathymetry')
end;
end
if(ftbverbose)
fprintf(['end : ' subname '\n'])
end;
end
fvcom_prepro/write_FVCOM_nested_forcing.m 0 → 100644
View file @ 378f9be5
This diff is collapsed.
fvcom_prepro/write_FVCOM_probes.m
View file @ 378f9be5
......@@ -12,7 +12,9 @@ function write_FVCOM_probes(nml_file, interval, probes)
% interval - output interval (in seconds)
% probe - struct of structs with fields whose names are the probe title.
% Each probe struct must have the following fields:
% node - unstructured grid node number.
% node - unstructured grid node ID (for all scalars except u
% and v).
% elem - unstructured grid element ID (for u and v only).
% file - file name for the current probe's output. If the
% 'variable' field is a cell array of multiple variables, each
% output file name will have the variable name appended (e.g.
......@@ -39,6 +41,7 @@ function write_FVCOM_probes(nml_file, interval, probes)
% EXAMPLE USAGE:
% probes.newlyn.file = 'newlyn_elev.dat';
% probes.newlyn.node = 1045;
% probes.newlyn.elem = 467;
% probes.newlyn.description = 'Elevation at Newlyn';
% probes.newlyn.variable = 'el';
% probes.newlyn.longname = 'Surface elevation (m)';
......@@ -65,6 +68,10 @@ function write_FVCOM_probes(nml_file, interval, probes)
%
% Revision history:
% 2014-02-10 - First version.
% 2015-01-14 - Fix export of the u and v locations to use the closest
% element instead of the node ID. Using the node ID for the u and v data
% will yield velocity time series miles away from the actual location of
% interest.
%
%==========================================================================
......@@ -103,7 +110,11 @@ for p = 1:np
file = fullfile(pathstr, name);
fprintf(f, '&NML_PROBE\n');
fprintf(f, ' PROBE_INTERVAL = "seconds=%.1f",\n', interval);
fprintf(f, ' PROBE_LOCATION = %i,\n', probes.(fnames{p}).node);
if any(strcmpi(vname, {'u', 'v'}))
fprintf(f, ' PROBE_LOCATION = %i,\n', probes.(fnames{p}).elem);
else
fprintf(f, ' PROBE_LOCATION = %i,\n', probes.(fnames{p}).node);
end
fprintf(f, ' PROBE_TITLE = "%s",\n', file);
% If we've got something which is vertically resolved, output the
% vertical levels.
......
fvcom_prepro/write_FVCOM_restart.m
View file @ 378f9be5
......@@ -330,11 +330,11 @@ for ii = 1:numvars
ss = 0:1 / (nt - 1):1; % scale from 0 to 1.
% Use the first modelled time step.
startdata = squeeze(data(:, :, 1));
td = indata.(fnames{vv}) - startdata;
for tt = 1:nt
if tt == 1
sfvdata(:, :, 1) = startdata;
else
td = indata.(fnames{vv}) - startdata;
sfvdata(:, :, tt) = startdata + (ss(tt) .* td);
end
end
......
utilities/read_netCDF_FVCOM.m
View file @ 378f9be5
function [data,selection] =read_netCDF_FVCOM(varargin)
function [data,selection] = read_netCDF_FVCOM(varargin)
% Function to extract data from a Netcdf file output from FVCOM.
%
% data = read_netCDF_FVCOM(varargin)
......@@ -74,10 +74,17 @@ function [data,selection] =read_netCDF_FVCOM(varargin)
% 'node_idx', node_idx, ...
% 'varnames', vars);
%
% BUGS:
% - When loading all times with the argument pair:
% 'time', -1
% the returned time series is nt - 1 (where nt is the number of time
% steps in the netCDF file). Not sure where this is broken, but
% probably around line 377.
%
% Author(s):
% Ricardo Torres - Plymouth Marine Laboratory 2012
% Hakeem Johnson - CH2M
% Pierre Cazenave - Plymouth Marine Laboratory
% Ricardo Torres - Plymouth Marine Laboratory 2012
% Hakeem Johnson - CH2M
% Pierre Cazenave - Plymouth Marine Laboratory
%
% Revision history:
% v0 March 2012
......@@ -191,18 +198,18 @@ try
Itime.Data(2) = netcdf.getVar(nc,Itime.ID,last_entry-1,1,'int32');
Itime2.Data(1) = netcdf.getVar(nc,Itime.ID+1,0,1,'int32');
Itime2.Data(2) = netcdf.getVar(nc,Itime.ID+1,last_entry-1,1,'int32');
[start_d(1),end_d(1)] = deal(double(Itime.Data(1))+time_offset,double(Itime.Data(end))+time_offset);
[start_d(2),end_d(2)] = deal(double(Itime2.Data(1)),double(Itime2.Data(end)));
start_date=sum(start_d.*[1 1/(24*60*60*1000)]); %hkj missing 1000 inserted
end_date = sum(end_d.*[1 1/(24*60*60*1000)]); %hkj missing 1000 inserted
var_time = netcdf.getVar(nc,Itime.ID,[0],[10],'double')+time_offset+...
netcdf.getVar(nc,Itime.ID+1,0,10,'double')./(24*600*6000) ;
DeltaT=median(diff(var_time));
var_time = start_date:DeltaT:(end_date-DeltaT);
catch me
if ftbverbose
warning('No ''Itime'' and/or ''Itime2'' variables, using less precise ''time'' instead.\n(%s)\n', me.message)
......@@ -277,7 +284,9 @@ for aa=1:length(varnames)
% Extract number of dimensions, lengths and names of all variables
%----------------------------------------------------------------------
% tic
fprintf('Processing variable %s: ', varnames{aa})
if ftbverbose
fprintf('Processing variable %s: ', varnames{aa})
end
% Tidy up the previous iteration's variables so we don't get confused.
clear dimName dimLength
......@@ -289,7 +298,7 @@ for aa=1:length(varnames)
else
netcdf.close(nc)
varargout{1} = 0;
fprintf('\n')
if ftbverbose; fprintf('\n'); end
error('Variable %s NOT found in file. Stopping. Check input variable names.', varnames{aa})
end
varID=netcdf.inqVarID(nc,vars{varidx(aa)});
......@@ -302,12 +311,18 @@ for aa=1:length(varnames)
if ftbverbose
if dd == 1
if length(dimids) == 1
fprintf('%i dimension: %s ', dimens, dimName{dd})
if ftbverbose
fprintf('%i dimension: %s ', dimens, dimName{dd})
end
else
fprintf('%i dimensions: %s ', dimens, dimName{dd})
if ftbverbose
fprintf('%i dimensions: %s ', dimens, dimName{dd})
end
end
else
fprintf('%s ', dimName{dd})
if ftbverbose
fprintf('%s ', dimName{dd})
end
end
end
end
......@@ -371,7 +386,7 @@ for aa=1:length(varnames)
start.(dimName{dd})=RestrictDims.idx{dimidx(dd)}(1)-1;
count.(dimName{dd})=length(RestrictDims.idx{dimidx(dd)});
stride.(dimName{dd})=timestrd;
else
for ss=1:length(RestrictDims.idx{dimidx(dd)})
start.(dimName{dd})(ss)=RestrictDims.idx{dimidx(dd)}(ss)-1;
......@@ -415,14 +430,14 @@ for aa=1:length(varnames)
read_start(find(~do_restrict))=start.(cc_names{find(~do_restrict)});
read_count(find(~do_restrict))=count.(cc_names{find(~do_restrict)});
read_stride(find(~do_restrict))=stride.(cc_names{find(~do_restrict)});
for cc=1:length(start.(cc_names{find(do_restrict)}))
read_start(find(do_restrict))=start.(cc_names{find(do_restrict)})(cc);
read_count(find(do_restrict))=count.(cc_names{find(do_restrict)})(cc);
read_stride(find(do_restrict))=stride.(cc_names{find(do_restrict)});
var_dump=netcdf.getVar(nc,varID,read_start,read_count,read_stride,'double');
eval([varnames{aa},'(cc,:)=var_dump;'])
clear var_dump
end
......@@ -470,7 +485,7 @@ for aa=1:length(varnames)
read_start(tt)=start.(cc_names{tt});
read_count(tt)=count.(cc_names{tt});
read_stride(tt)=stride.(cc_names{tt});
end
% check if time is one of them
......@@ -487,7 +502,7 @@ for aa=1:length(varnames)
read_count(find(do_restrict))=count.(cc_names{find(do_restrict)})(cc);
read_stride(find(do_restrict))=stride.(cc_names{find(do_restrict)});
var_dump=netcdf.getVar(nc,varID,read_start,read_count,read_stride,'double');
switch dimName(find(do_restrict))
case 'node' | 'nele'
eval([varnames{aa},'(cc,:,:)=var_dump;'])
......@@ -564,7 +579,7 @@ for aa=1:length(varnames)
read_count(2)=count.(cc_names{2})(cc);
read_stride(2)=stride.(cc_names{2});
var_dump=netcdf.getVar(nc,varID,read_start,read_count,read_stride,'double');
eval([varnames{aa},'(kk,cc)=var_dump;'])
clear var_dump
end
......@@ -575,7 +590,7 @@ for aa=1:length(varnames)
eval(['selection.',varnames{aa},'.start=start;'])
eval(['selection.',varnames{aa},'.count=count;'])
case 3 % three dimension to restrict!
% but the variable can have more than 2 dimensions
switch dimens
......@@ -587,7 +602,7 @@ for aa=1:length(varnames)
min(sum(count.(cc_names{2})),dimLength(2)),...
min(sum(count.(cc_names{3})),dimLength(3))];
end
eval([varnames{aa},'=nan(rr);'])
% check if time is one of them
if isempty(find(dimidx==5));disp('This won''t work, try again');return;end
......@@ -596,7 +611,7 @@ for aa=1:length(varnames)
read_start(do_time)=start.(cc_names{do_time});
read_count(do_time)=count.(cc_names{do_time});
read_stride(do_time)=stride.(cc_names{do_time});
% search for the non_restrictive variable
% cc=1
% while ~(length( start.(cc_names{cc}))==1);cc=cc+1;end
......@@ -608,17 +623,17 @@ for aa=1:length(varnames)
else
do_one=do_other(2);do_two=do_other(1);
end
for cc=1:length(start.(cc_names{do_one}))
read_start(do_one)=start.(cc_names{do_one})(cc);
read_count(do_one)=count.(cc_names{do_one})(cc);
read_stride(do_one)=stride.(cc_names{do_one});
for pp=1:length(start.(cc_names{do_two}))
read_start(do_two)=start.(cc_names{do_two})(pp);
read_count(do_two)=count.(cc_names{do_two})(pp);
read_stride(do_two)=stride.(cc_names{do_two});
var_dump=netcdf.getVar(nc,varID,read_start,read_count,read_stride,'double');
eval([varnames{aa},'(pp,cc,:)=var_dump;'])
end
......@@ -627,7 +642,7 @@ for aa=1:length(varnames)
eval(['selection.',varnames{aa},'.start=start;'])
eval(['selection.',varnames{aa},'.count=count;'])
case 0 % there are NO dimensions to restrict and 3 dimensions haven't been coded yet!!
for nn=1:length(cc_names)
read_start(nn)=start.(cc_names{nn});
read_count(nn)=count.(cc_names{nn});
......