/[MITgcm]/MITgcm_contrib/gmaze_pv/subfct/latlon2ingrid_netcdf.m
ViewVC logotype

Contents of /MITgcm_contrib/gmaze_pv/subfct/latlon2ingrid_netcdf.m

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


Revision 1.1 - (show annotations) (download)
Tue Jan 30 22:10:34 2007 UTC (18 years, 6 months ago) by gmaze
Branch: MAIN
Add eg for bin2cdf conversion for both lat-lon cs grid

1 %
2
3 function latlon2ingrid_netcdf(pathname,pathout, ...
4 stepnum,fpref,otab, ...
5 lon_c, lon_u, ...
6 lat_c, lat_v, ...
7 z_c, z_w, ...
8 subname, ...
9 lonmin,lonmax,latmin,latmax,depmin,depmax);
10
11 irow=strmatch({fpref},otab(:,1),'exact');
12 if length(irow) ~= 1
13 fprintf('Bad irow value in latlon2ingrid_netcdf2\n');
14 return
15 end
16 loc=otab{irow,3};
17 id=otab{irow,4};
18 units=otab{irow,5};
19 dimspec=otab{irow,2};
20 if strmatch(id,'unknown_id','exact')
21 id = fpref;
22 end
23 fprintf('Field %s, loc=%s, id=%s, units=%s, dimspec=%s\n',fpref,loc,id,units,dimspec);
24 wordlen=otab{irow,6};
25 if wordlen == 4
26 numfmt='float32';
27 end
28 if wordlen == 8
29 numfmt='float64';
30 end
31 %numfmt='float64';
32 %wordlen =8;
33
34 %ilo_c=min(find(lon_c >= lonmin & lon_c <= lonmax));
35 %ilo_u=min(find(lon_u >= lonmin & lon_u <= lonmax));
36 %ihi_c=max(find(lon_c >= lonmin & lon_c <= lonmax));
37 %ihi_u=max(find(lon_u >= lonmin & lon_u <= lonmax));
38 ilo_c=min(find(lon_c-180 >= lonmin & lon_c-180 <= lonmax));
39 ilo_u=min(find(lon_u-180 >= lonmin & lon_u-180 <= lonmax));
40 ihi_c=max(find(lon_c-180 >= lonmin & lon_c-180 <= lonmax));
41 ihi_u=max(find(lon_u-180 >= lonmin & lon_u-180 <= lonmax));
42 jlo_c=min(find(lat_c >= latmin & lat_c <= latmax));
43 jlo_v=min(find(lat_v >= latmin & lat_v <= latmax));
44 jhi_c=max(find(lat_c >= latmin & lat_c <= latmax));
45 jhi_v=max(find(lat_v >= latmin & lat_v <= latmax));
46 klo_w=min(find(z_w >= depmin & z_w <= depmax));
47 khi_w=max(find(z_w >= depmin & z_w <= depmax));
48 klo_c=min(find(z_c >= depmin & z_c <= depmax));
49 khi_c=max(find(z_c >= depmin & z_c <= depmax));
50
51 fnam=sprintf('%s.%10.10d.data',fpref,stepnum);
52 if loc == 'c'
53 ilo=ilo_c;
54 ihi=ihi_c;
55 jlo=jlo_c;
56 jhi=jhi_c;
57 klo=klo_c;
58 khi=khi_c;
59 lon=lon_c;
60 lat=lat_c;
61 dep=-z_c;
62 end
63 if loc == 'u'
64 ilo=ilo_u;
65 ihi=ihi_u;
66 jlo=jlo_c;
67 jhi=jhi_c;
68 klo=klo_c;
69 khi=khi_c;
70 lon=lon_u;
71 lat=lat_c;
72 dep=-z_c;
73 end
74 if loc == 'v'
75 ilo=ilo_c;
76 ihi=ihi_c;
77 jlo=jlo_v;
78 jhi=jhi_v;
79 klo=klo_c;
80 khi=khi_c;
81 lon=lon_c;
82 lat=lat_v;
83 dep=-z_c;
84 end
85 if loc == 'w'
86 ilo=ilo_c;
87 ihi=ihi_c;
88 jlo=jlo_c;
89 jhi=jhi_c;
90 klo=klo_w;
91 khi=khi_w;
92 lon=lon_c;
93 lat=lat_v;
94 dep=-z_w;
95 end
96
97 nx=1;ny=1;nz=1;
98 if strmatch(dimspec,'xyz','exact');
99 nx=length(lon);
100 ny=length(lat);
101 nz=length(dep);
102 end
103 if strmatch(dimspec,'xy','exact');
104 nx=length(lon);
105 ny=length(lat);
106 end
107
108 if klo > nz
109 klo = nz;
110 end
111 if khi > nz
112 khi = nz;
113 end
114
115 phiXYZ=zeros(ihi-ilo+1,jhi-jlo+1,khi-klo+1,'single');
116 disp(strcat('in:',pathname,fnam))
117 %[klo khi khi-klo+1]
118
119 % Read a single level (selected by k)
120 for k = klo : khi
121 fid = fopen(strcat(pathname,fnam),'r','ieee-be');
122 fseek(fid,(k-1)*nx*ny*wordlen,'bof');
123 phi = fread(fid,nx*ny,numfmt);
124 %whos phi, [k nx ny]
125 phiXY = reshape(phi,[nx ny]);
126 phiXY = phiXY(ilo:ihi,jlo:jhi);
127 phiXYZ(:,:,k) = phiXY;
128 %phiXYZ(100,100,k)
129 fclose(fid);
130 end
131
132 %%%clear phi;
133 %%%clear phiXY;
134 phiXYZ(find(phiXYZ==0))=NaN;
135
136 if subname == ' '
137 %outname=sprintf('%s.nc',id);
138 outname = sprintf('%s.nc',otab{irow,1});
139 else
140 %outname=sprintf('%s_%s.nc',subname,id);
141 outname = sprintf('%s.%s.nc',otab{irow,1},subname);
142 %outname = sprintf('%s.%s.nc',strcat(otab{irow,1},'s'),subname);
143
144 end
145 nc = netcdf(strcat(pathout,outname),'clobber');
146 %disp(strcat(pathout,outname))
147
148 nc('X')=ihi-ilo+1;
149 nc('Y')=jhi-jlo+1;
150 nc('Z')=khi-klo+1;
151
152 nc{'X'}='X';
153 nc{'Y'}='Y';
154 nc{'Z'}='Z';
155
156 nc{'X'}.uniquename='X';
157 nc{'X'}.long_name='longitude';
158 nc{'X'}.gridtype=ncint(0);
159 nc{'X'}.units='degrees_east';
160 nc{'X'}(:) = lon(ilo:ihi);
161
162 nc{'Y'}.uniquename='Y';
163 nc{'Y'}.long_name='latitude';
164 nc{'Y'}.gridtype=ncint(0);
165 nc{'Y'}.units='degrees_north';
166 nc{'Y'}(:) = lat(jlo:jhi);
167
168 nc{'Z'}.uniquename='Z';
169 nc{'Z'}.long_name='depth';
170 nc{'Z'}.gridtype=ncint(0);
171 nc{'Z'}.units='m';
172 nc{'Z'}(:) = dep(klo:khi);
173
174 ncid=id;
175 nc{ncid}={'Z' 'Y' 'X'};
176 nc{ncid}.missing_value = ncdouble(NaN);
177 nc{ncid}.FillValue_ = ncdouble(0.0);
178 nc{ncid}(:,:,:) = permute(phiXYZ,[3 2 1]);
179 nc{ncid}.units=units;
180
181 close(nc);

  ViewVC Help
Powered by ViewVC 1.1.22