/[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.2 - (hide annotations) (download)
Wed Sep 19 15:24:41 2007 UTC (17 years, 10 months ago) by gmaze
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +8 -1 lines
General Update

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

  ViewVC Help
Powered by ViewVC 1.1.22