Warning, /utils/matlab/gmt/rdnctiles_oldflat.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
5e40d2151d Ed H*0001 function [res] = rdnctiles_oldflat(fall,vit,dlev)
0002
0003 % Function [res] = rdnctiles_oldflat(fall,vit,dlev)
0004 %
0005 % INPUTS
0006 % fall cell array of file names
0007 % vit struct containing variable and time information
0008 % vit.tdname : "time" dim name (DEF: 'T')
0009 % vit.tdname : "time" coord var name (DEF: 'iter')
0010 % vit.vars.(vname) : "time" values for each var
0011 %
0012 % dlev debug level
0013 %
0014 % OUTPUTS
0015 % res struct containing results
0016 %
0017 %
0018 % Ed Hill
0019 fi = 1;
0020 nc = netcdf(fall{fi},'read');
0021 allatts = ncnames(att(nc));
0022 if not(isempty( intersect(allatts,'exch3_ver') ))
0023 res.att.exch_ver = 3;
0024 elseif not(isempty( intersect(allatts,'exch2_tbasex') ))
0025 res.att.exch_ver = 2;
0026 else
0027 res.att.exch_ver = 1;
0028 end
0029 res.att.Nx = nc.('Nx')(:);
0030 res.att.Ny = nc.('Ny')(:);
0031 res.att.Nr = nc.('Nr')(:);
0032 res.att.Ne = min(res.att.Nx, res.att.Ny);
0033 nc = close(nc);
0034 for fi = 2:length(fall)
0035 nc = netcdf(fall{fi},'read');
0036 % Verify that all files belong to a domain that has the same overall
0037 % shape
0038 if ( (res.att.Nx ~= nc.('Nx')(:)) ...
0039 || (res.att.Ny ~= nc.('Ny')(:)) ...
0040 || (res.att.Nr ~= nc.('Nr')(:)) )
0041 error(['file ''' fall{fi} ''' does not belong to a domain ' ...
0042 'that has the same overall shape (Nx,Ny,Nr) as the ' ...
0043 'other files']);
0044 end
0045 nc = close(nc);
0046 end
0047
0048 for fi = 1:length(fall)
0049 if dlev > 10
0050 disp([' Opening : ' char(fall{fi}) ]);
0051 end
0052 nc = netcdf(fall{fi},'read');
0053
0054 % Get all the variables
0055 if dlev > 10
0056 fprintf(1,' reading : ');
0057 end
0058 vread = fields(vit.vars);
0059 for iv = 1:length(vread)
0060 if isempty(nc{char(vread{iv})})
0061 % disp([' warning: no var "',vread{iv},'" in "',fall{fi},'"']);
0062 continue
0063 end
0064 if dlev > 10
0065 fprintf(1,[' ' char(vread{iv}) ]);
0066 end
0067
0068 % Get the rank of the time dimension
0069 trank = 0;
0070 dnames = ncnames(dim(nc{vread{iv}}));
0071 if strcmp( ncnames(recdim(nc)), vit.tdname )
0072 m = regexp(dnames, [ '^' vit.tdname '$'] );
0073 for i = 1:length(m)
0074 if not(isempty(m{i}))
0075 trank = i;
0076 break
0077 end
0078 end
0079 end
0080
0081 % get the rank of the X,Y dims
0082 xrank = 0;
0083 yrank = 0;
0084 for i = 1:length(dnames)
0085 if strcmp(dnames{i},'X') || strcmp(dnames{i},'Xp1')
0086 xrank = i;
0087 end
0088 if strcmp(dnames{i},'Y') || strcmp(dnames{i},'Yp1')
0089 yrank = i;
0090 end
0091 end
0092
0093 % get the corresponding file-local indicies and global-assembly
0094 % indicies along the time dimension
0095 loc_times = nc{vit.tvname}(:);
15db930790 Ed H*0096 [v,ind1,ind2] = intersect( vit.vars.(vread{iv}), loc_times );
5e40d2151d Ed H*0097
0098 % only read the desired time values based on:
0099 % the local "kt" indicies and
0100 % the global "tk" indicies
0101 vdims = dim(nc{vread{iv}});
0102 sz_glob = [];
0103 indstr = '';
0104 for i = 1:length(dnames)
0105 if i > 1
0106 indstr = [ indstr ',' ];
0107 end
0108 if i == trank
0109 indstr = [ indstr 'kt' ];
0110 sz_glob = [ sz_glob length(vit.vars.(vread{iv})) ];
0111 elseif i == xrank
0112 indstr = [ indstr 'kx' ];
0113 sz_glob = [ sz_glob res.att.Nx ];
0114 elseif i == yrank
0115 indstr = [ indstr 'ky' ];
0116 sz_glob = [ sz_glob res.att.Ny ];
0117 else
0118 indstr = [ indstr ':' ];
0119 sz_glob = [ sz_glob vdims{i}(:) ];
0120 end
0121 end
0122 rindstr = fliplr(indstr);
0123 sz_glob = fliplr(sz_glob);
0124
0125 % If the global variable doesn't exist, then declare it and reserve
0126 % the necessary space
0127 if not(isfield(res,'var')) || not(isfield(res.var,vread{iv}))
0128 if length(sz_glob) == 1
0129 res.var.(vread{iv}) = zeros(sz_glob,1);
0130 else
0131 res.var.(vread{iv}) = zeros(sz_glob);
0132 end
0133 end
0134
0135 % file-local X,Y indicies
0136 kx = [ 1:nc.('sNx')(:) ];
0137 ky = [ 1:nc.('sNy')(:) ];
0138 % global X,Y indicies
23ca324d5d Ed H*0139 if res.att.exch_ver == 2
0140 xk = nc.('exch2_txglobalo')(:) - 1 + kx;
0141 yk = nc.('exch2_tyglobalo')(:) - 1 + ky;
0142 elseif res.att.exch_ver == 1
0143 xk = (nc.('bi')(:) - 1)*nc.('sNx')(:) + kx;
0144 yk = (nc.('bj')(:) - 1)*nc.('sNy')(:) + ky;
0145 end
15db930790 Ed H*0146 if length(ind1) > 0
0147 for jj = 1:length(ind1)
0148 kt = ind2(jj);
5e40d2151d Ed H*0149 eval([ 'tmpv = nc{vread{iv}}(' indstr ');' ]);
0150 sz = size(tmpv);
0151 nd = length(sz);
15db930790 Ed H*0152 tk = ind1(jj);
5e40d2151d Ed H*0153 comm = [ 'res.var.(char(vread{iv}))(' ...
0154 rindstr ') = permute(tmpv,[nd:-1:1]);' ];
0155 eval(comm);
0156 end
0157 else
0158 eval([ 'tmpv = nc{vread{iv}}(' indstr ');' ]);
0159 sz = size(tmpv);
0160 nd = length(sz);
0161 comm = [ 'res.var.(char(vread{iv}))(' ...
0162 rindstr ') = permute(tmpv,[nd:-1:1]);' ];
0163 eval(comm);
0164 end
0165 end
0166 if dlev > 10
0167 fprintf(1,'\n');
0168 end
0169
0170 nc = close(nc);
0171 end
0172
