FpOutput: flex_read_recepconc.m

File flex_read_recepconc.m, 2.2 KB (added by sureshraja1, 8 years ago)
Line 
1%----------------------------------------------------------
2% FLEXPART Routine: read_recpconc
3%----------------------------------------------------------
4% description
5%   - reads concentrations at a point(s) as defined in
6%     the file RECEPTORS
7% input
8%   - head (data.frame output from read_header.r)
9%   - optional: data units (1=conc (default), 2=pptv)
10%   - time_ret (vector of julian dates from read_header.r)
11%   - numreceptor (number of receptors)
12% output
13%   - storname - receptor names
14%   - storloc - receptor locations
15%   - concentrations at the receptors
16%Author: Suresh Raja
17%----------------------------------------------------------
18
19function [receptorname,storeloc,conc]=flex_read_recepconc(header,unit,numreceptor,fin)
20 
21  %set defaults
22  if unit==1
23     filein='receptor_conc'
24  end
25  if unit==2
26     filein='receptor_pptv'
27  end
28 
29  time_ret=1:length(header.dates);
30 
31  fin=[header.path filein]
32  fb=fopen(fin);
33     
34if fb<=0
35   fprintf('******************************\n could not open file %s !\n******************************\n',fin);
36else
37
38  %receptorname=repmat({'Sample Text'}, numreceptor, 1);
39  receptorname=[];
40  rname=[];
41  storename=[];
42  %receptorloc=repmat(0,numreceptor,2);
43  receptorloc=[];
44  storeloc=[];
45  %conc=repmat(0,length(time_ret),numreceptor);
46  conc=[];
47 
48  %# read receptor names
49   rl=fread(fb,1,'int32');
50 
51  for nr=1:numreceptor
52      rname=char(transpose(fread(fb,16,'char')));
53      %rname=strrep(rname,' ',''); -this removes the spaCES-
54      %receptorname={receptorname;rname};
55      receptorname=[receptorname;rname];
56  end
57   
58  rl=fread(fb,2,'int32');
59 
60  %# read receptor locations
61  for nr = 1:numreceptor
62      receptorloc=fread(fb,2,'float')';
63      storeloc=[storeloc;receptorloc];
64  end
65  rl=fread(fb,2,'int32');
66 
67  %#--- loop over dates
68  dat_cnt=0;
69  for i = 1:length(time_ret)
70     dat_cnt=dat_cnt+1;
71     ns_cnt=0;
72     itime=fread(fb,1,'int32');
73     rl=fread(fb,2,'int32');
74
75  %  # read all species
76     for ns = 1:header.nspec
77         ns_cnt=ns_cnt+1;
78  %    # read all receptors
79      rec_dump=fread(fb,numreceptor,'float')';
80      rl=fread(fb,2,'int32');
81      conc=[conc;rec_dump];
82     end
83
84 end
85 
86 fclose(fb)
87
88end
hosted by ZAMG