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

Annotation 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 - (hide annotations) (download)
Tue Jan 30 22:10:34 2007 UTC (18 years, 5 months ago) by gmaze
Branch: MAIN
Add eg for bin2cdf conversion for both lat-lon cs grid

1 gmaze 1.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