/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_IO/read_nctiles.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/gcmfaces_IO/read_nctiles.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.7 - (hide annotations) (download)
Fri Mar 13 16:16:09 2015 UTC (10 years, 4 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65x, checkpoint65r, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, HEAD
Changes since 1.6: +42 -9 lines
contribution from M. Buckley: allow for
-specification of kk without specification of tt
-a range of contiguous indicies for kk  rather than just a single level

1 gforget 1.3 function [fld]=read_nctiles(fileName,fldName,varargin);
2 gforget 1.1 %usage: fld=read_nctiles(fileName); reads full field using fileName as field name
3     %usage: fld=read_nctiles(fileName,fldName); reads full field (all depths, all times)
4     %usage: fld=read_nctiles(fileName,fldName,tt); reads 3D or 2D field, at time index tt (all depths)
5 gforget 1.7 %usage: fld=read_nctiles(fileName,fldName,tt,kk); reads 3D field, at depth index(es) kk, at time index tt
6 gforget 1.1
7     gcmfaces_global;
8 gforget 1.7 nz=length(mygrid.RC);
9 gforget 1.1
10     if nargin==1;
11     tmp1=['/' fileName];
12 gforget 1.6 tmp2=strfind(tmp1,filesep);
13 gforget 1.1 fldName=tmp1(tmp2(end)+1:end);
14     end;
15 gforget 1.7 if nargin>2;
16     tt=varargin{1};
17     else;
18     tt=[];
19     end;
20     if nargin>3;
21     kk=varargin{2};
22     if min(kk)<1 | max(kk)>nz
23     error(['kk must be an integer between 1 and ' num2str(nz)])
24     return
25     end
26     if length(kk)>1
27     if ~isempty(find(diff(kk)~=1))
28     error(['kk must be a contigous range of indices'])
29     return
30     end
31     end
32     else;
33     kk=[];
34     end;
35 gforget 1.1
36 gforget 1.5 %A) determine map of tile indices (if not already done)
37    
38     global nctiles;
39    
40     test1=isempty(nctiles);
41     test2=0;
42     if ~test1;
43     tmp1=dir([fileName '*']);
44     test2=(length(tmp1)~=length(nctiles.no));
45     end;
46    
47     if test1|test2;
48     %build map of tile indices
49     if (length(dir([fileName '*']))==mygrid.nFaces);
50     nctiles.map=NaN*mygrid.XC;
51     for ff=1:mygrid.nFaces; nctiles.map{ff}(:)=ff; end;
52     else;
53     fileIn=sprintf('%s.%04d.nc',fileName,1);
54     nc=netcdf.open(fileIn,0);
55     vv = netcdf.inqVarID(nc,fldName);
56     [varname,xtype,dimids,natts]=netcdf.inqVar(nc,vv);
57     tileSize=zeros(1,2);
58     [tmp,tileSize(1)] = netcdf.inqDim(nc,dimids(1));
59     [tmp,tileSize(2)] = netcdf.inqDim(nc,dimids(2));
60     netcdf.close(nc);
61     nctiles.map=gcmfaces_loc_tile(tileSize(1),tileSize(2));
62     end;
63     %determine tile list
64     nctiles.no=unique(convert2vector(nctiles.map));
65     nctiles.no=nctiles.no(~isnan(nctiles.no));
66     %determine location of each tile
67     for ff=1:length(nctiles.no);
68     for gg=1:mygrid.nFaces;
69     [tmpi,tmpj]=find(nctiles.map{gg}==ff);
70     if ~isempty(tmpi);
71     nctiles.f{ff}=gg;
72     nctiles.i{ff}=[min(tmpi(:)):max(tmpi(:))];
73     nctiles.j{ff}=[min(tmpj(:)):max(tmpj(:))];
74     end;
75     end;
76     end;
77 gforget 1.2 end;
78 gforget 1.1
79 gforget 1.5 %B) the file read operation itself
80 gforget 1.1
81 gforget 1.5 for ff=1:length(nctiles.no);
82 gforget 1.2
83 gforget 1.5 %read one tile
84 gforget 1.2 fileIn=sprintf('%s.%04d.nc',fileName,ff);
85 gforget 1.1 nc=netcdf.open(fileIn,0);
86 gforget 1.5
87 gforget 1.1 vv = netcdf.inqVarID(nc,fldName);
88     [varname,xtype,dimids,natts]=netcdf.inqVar(nc,vv);
89    
90 gforget 1.2 [dimname,siz(1)] = netcdf.inqDim(nc,dimids(1));
91     [dimname,siz(2)] = netcdf.inqDim(nc,dimids(2));
92     siz=[siz length(mygrid.RC)];
93 gforget 1.1
94     if ~isempty(tt);
95 gforget 1.4 t0=tt(1)-1;
96     nt=tt(end)-tt(1)+1;
97 gforget 1.1 if length(dimids)==3;
98 gforget 1.4 start=[0 0 t0];
99     count=[siz(1) siz(2) nt];
100 gforget 1.1 elseif isempty(kk);
101 gforget 1.4 start=[0 0 0 t0];
102     count=[siz(1) siz(2) siz(3) nt];
103 gforget 1.1 else;
104 gforget 1.7 start=[0 0 kk(1)-1 t0];
105     count=[siz(1) siz(2) length(kk) nt];
106     end;
107     else
108     [dimname,nt] = netcdf.inqDim(nc,dimids(end));
109     if length(dimids)==2;
110     start=[0 0];
111     count=[siz(1) siz(2)];
112     elseif length(dimids)==3;
113     start=[0 0 0];
114     count=[siz(1) siz(2) nt];
115     elseif isempty(kk);
116     start=[0 0 0 0];
117     count=[siz(1) siz(2) siz(3) nt];
118     else;
119     start=[0 0 kk(1)-1 0];
120     count=[siz(1) siz(2) length(kk) nt];
121 gforget 1.1 end;
122 gforget 1.2 end;
123 gforget 1.7 fldTile=netcdf.getVar(nc,vv,start,count);
124     fldTile=squeeze(fldTile);
125 gforget 1.5 netcdf.close(nc);
126    
127     %initialize fld (full gcmfaces object)
128 gforget 1.2 if ff==1;
129     siz=[1 1 size(fldTile,3) size(fldTile,4)];
130     fld=NaN*repmat(mygrid.XC,siz);
131     end;
132    
133 gforget 1.5 %place tile in fld
134     fld{nctiles.f{ff}}(nctiles.i{ff},nctiles.j{ff},:,:)=fldTile;
135 gforget 1.1
136     end;%for ff=1:mygrid.nFaces;
137    
138    
139 gforget 1.7

  ViewVC Help
Powered by ViewVC 1.1.22