FpOutput: flex_header.m

File flex_header.m, 7.1 KB (added by sureshraja1, 3 years ago)

flex_header function

Line 
1function [header, fail] = flex_header(pathname, nest, readp, calcarea, datespath)
2%===========================================
3% flex_header.m
4%-------------------------------------------
5% input
6%   - pathname: directory path
7%   - readp   : read release points (0/1)
8% output
9%   - header object
10%   - fail indicator (-1=fail, 0=success)
11%-------------------------------------------
12% FLEXPART matlab import routines
13%-------------------------------------------
14% last changes: HSO, 06-Jun-2007
15%===========================================
16 
17  % anonymous helper function to skip some bytes
18  skip = @(fid) fread(fid,2,'int32');
19
20  % define the output object
21  if nargin==0
22    pathname='.';
23  elseif nargin==1
24    nest=0; 
25    readp=1; % default: do not read release points
26    calcarea=0;
27  end % if nargin
28 
29  % create file name for header
30  if (pathname(length(pathname))~='/')
31    pathname=[pathname ''];
32  end % if pathname
33
34
35  if nargin<=4
36     datespath=pathname;
37  else
38    datespath=deblank(datespath);
39  end   
40 
41  header.path=pathname;
42  if nest==0
43    fhd=fopen([header.path 'header']);
44    %fprintf('\nreading header');
45    header.nest=0;
46  else
47    fhd=fopen([header.path 'header_nest']);
48    header.nest=1;
49    fprintf('\nreading nested header');
50  end
51
52  % default settings
53  header.decayconstant=0;
54
55  % start reading if header was opened
56  if fhd~=-1
57    rl=fread(fhd,1,'int32'); 
58    header.open=1;
59    header.ibdate=fread(fhd,1,'int32');
60    header.ibtime=fread(fhd,1,'int32');
61    pivot_vec=(datevec(date));
62    %header.simulationstart=datenum([num2str(header.ibdate)  num2str(header.ibtime,'%06d')],'yyyymmddhhMMss',pivot_vec(1)-50);
63    header.simulationstart=datenum([num2str(header.ibdate)  num2str(header.ibtime,'%06d')],'yyyymmddhhMMss');
64    header.ver=char(transpose(fread(fhd,13,'char')));
65    skip(fhd);
66    header.loutstep=fread(fhd,1,'int32');
67    header.loutaver=fread(fhd,1,'int32');
68    header.loutsample=fread(fhd,1,'int32');
69    skip(fhd);
70    header.outlon0=fread(fhd,1,'float');
71    header.outlat0=fread(fhd,1,'float');
72    header.numxgrid=fread(fhd,1,'int32');
73    header.numygrid=fread(fhd,1,'int32');   
74    header.dxout=fread(fhd,1,'float');
75    header.dyout=fread(fhd,1,'float',8);
76    header.numzgrid=fread(fhd,1,'int32');
77    header.outheight=fread(fhd,header.numzgrid,'float');
78    header.outheight=cat(1,[0],header.outheight);
79    skip(fhd);
80    header.jjjjmmdd=fread(fhd,1,'int32');   
81    header.hhmmss=fread(fhd,1,'int32',8);   
82    header.nspec=fread(fhd,1,'int32')/3;
83    %%%%%%%%%%%%%%%%%%%%%%% diff verisons
84%     skip(fhd);
85    header.numpointspec=fread(fhd,1,'int32',8);
86    %%%%%%%%%%%%%%%%%%%%%%% diff verisons
87    % read the species
88    for i=1:header.nspec
89      one=fread(fhd,1,'int32'); % input skipped ??
90      wd=fread(fhd,10,'char');  % input skipped ??
91      skip(fhd);
92      one=fread(fhd,1,'int32'); % input skipped ??
93      dd=fread(fhd,10,'char');  % input skipped ??
94      skip(fhd);
95      header.numzgrid=fread(fhd,1,'int32'); 
96      header.species(i,1:10)=char(transpose(fread(fhd,10,'char')));
97      skip(fhd);
98    end % nspec
99 
100    % read the release points
101    header.numpoint=fread(fhd,1,'int32');
102      % initialise fields
103
104    header.ireleasestart=zeros(header.numpoint,1);
105    header.ireleaseend=zeros(header.numpoint,1);
106    header.kindz=zeros(header.numpoint,1);
107    header.xp1=zeros(header.numpoint,1);
108    header.yp1=zeros(header.numpoint,1);
109    header.xp2=zeros(header.numpoint,1);
110    header.yp2=zeros(header.numpoint,1);
111    header.zpoint1=zeros(header.numpoint,1);
112    header.zpoint2=zeros(header.numpoint,1);
113    header.npart=zeros(header.numpoint,1);
114%    for i=1:header.numpoint,
115%      header.compoint(i,1:45)='                                             ';
116%    end
117    header.xmass=zeros(header.numpoint,header.nspec);
118   
119    % read release info if requested
120    if readp==1
121      rl=fread(fhd,1,'int32'); 
122      for i=1:header.numpoint
123        rl=fread(fhd,1,'int32'); 
124        %header.ireleasestart(i)=fread(fhd,1,'int32');
125        %header.ireleaseend(i)=fread(fhd,1,'int32');
126        header.ireleasestart(i)=header.simulationstart+(fread(fhd,1,'int32'))/24/3600;
127        header.ireleaseend(i)=header.simulationstart+(fread(fhd,1,'int32'))/24/3600;
128        %%%%%%%%%%%%%%%% 2 diffs!
129%        header.kindz(i)=fread(fhd,1,'int32');
130        header.kindz(i)=fread(fhd,1,'int16');
131        %%%%%%%%%%%%%%%% 2 diffs!
132        skip(fhd);
133        header.xp1(i)=fread(fhd,1,'float');
134        header.yp1(i)=fread(fhd,1,'float');
135        header.xp2(i)=fread(fhd,1,'float');
136        header.yp2(i)=fread(fhd,1,'float');
137        header.zpoint1(i)=fread(fhd,1,'float');
138        header.zpoint2(i)=fread(fhd,1,'float');
139        skip(fhd);
140        header.npart(i)=fread(fhd,1,'int32');
141        header.mpart(i)=fread(fhd,1,'int32');
142        rl=fread(fhd,2,'int32');
143        header.compoint(i,:)=char(transpose(fread(fhd,45,'char')));
144        rl=fread(fhd,1,'int32');
145        for j=1:header.nspec,
146          rl=fread(fhd,1,'int32');
147          header.xmass(i,j)=fread(fhd,1,'float');
148          skip(fhd);
149          header.xmass(i,j)=fread(fhd,1,'float');
150          skip(fhd);
151          header.xmass(i,j)=fread(fhd,1,'float');
152          rl=fread(fhd,1,'int32');
153        end
154      end % numpoint
155    else
156      % skip reading points, fill data structure with zeroes
157      fseek(fhd,119*header.numpoint+(header.nspec*36)*header.numpoint+4,'cof');
158    end % readp
159
160   
161    rl=fread(fhd,1,'int32');
162    %%% check if ind_rel, ind_source is written out
163   
164    header.method=fread(fhd,1,'int32');
165    header.lsubgrid=fread(fhd,1,'int32');
166    header.lconvection=fread(fhd,1,'int32');
167    if rl==20
168    header.ind_source=fread(fhd,1,'int32');
169    header.ind_receptor=fread(fhd,1,'int32'); 
170    end
171    skip(fhd);
172    header.nageclass=fread(fhd,1,'int32');   
173    header.lage(1:header.nageclass)= fread(fhd,header.nageclass,'int32');
174    skip(fhd);
175
176    header.oroout=zeros(header.numxgrid,header.numygrid);
177    for ix=1:header.numxgrid,
178      header.oroout(ix,:)=fread(fhd,header.numygrid,'float');
179      skip(fhd);
180    end% fread
181
182    fclose(fhd);
183    fail=0;
184 
185    % try reading the dates file
186    f=([datespath 'dates']);
187    if exist(f)
188      dat=importdata([datespath 'dates']);
189      header.dates=(datenum(num2str(dat),'yyyymmddHHMMss'));
190    else
191       header.dates=NaN;
192       fprintf('Could not read dates file: %s\n',f);
193    end
194   
195    header.lonp = [header.outlon0+header.dxout/2:header.dxout:header.outlon0+header.dxout*header.numxgrid];
196    header.latp = [header.outlat0+header.dyout/2:header.dyout:header.outlat0+header.dyout*header.numygrid];
197    header.zp= (header.outheight(1:length(header.outheight)-1)+header.outheight(2:length(header.outheight)))./2;
198    %    header.lonp = [header.outlon0:header.dxout:header.outlon0+header.dxout*header.numxgrid-header.dxout];
199%    header.latp = [header.outlat0:header.dyout:header.outlat0+header.dyout*header.numygrid-header.dyout];
200   
201    if calcarea==1
202      header.area  = calculate_grid_area_eff(header);
203    end
204  else % if header file not opened
205 
206    fail=-1;
207    header.open=0; % this is the only field in header!
208
209  end % if header file opened
210
211end % of function flex_header
212
213%fin
hosted by ZAMG