FpOutput: readpart10.m

File readpart10.m, 4.6 KB (added by ignacio, 4 months ago)

matlab function to read particle output

Line 
1function output=readpart10(directory,i,numpart,varargin)
2% reads particle files from FLEXPART into structure output
3% inputs:
4% directory [string]: path to part files
5%       e.g. directory='/Users/userid/flexpart/FLEXPART8.1/output/'
6% i [integer or string]: order of particle file
7%       e.g. 1 for the firstdate in 'dates' file or 
8%       '19800103120000' to read file='partposit_19800103120000'
9% numpart [integer or string]: number of particles to be read.
10%       If unknown input any string. This will determine numpart from the size
11%       of the file and return total nr of particles in the file instead of
12%       the output given below.
13% varargin [integer]: optional number of species nspec. Default is 1.
14%
15% output:
16%     outheader: [double] time record in seconds since simulation start
17%        npoint: [numpart x 1 double] release point ID-number for each particle
18%           xyz: [numpart x 3 double] position coordinates
19%       itramem: [numpart x 1 double] relase times of each particle
20%          vars: [numpart x 7 double] particle parameters
21%         xmass: [numpart x nspec double] particle masses
22
23%xyz (from flexpart fortran partout subroutine):
24%1  xlon    =   longitude
25%2  ylat    =   latitude
26%3  ztra1(i)=   particle heigth in m
27%vars:
28%1  topo    =   topography heigth in m
29%2  pvi     =   potential vorticity
30%3  qvi     =   humodity
31%4  rhoi    =   air density
32%5  hmixi   =   'PBL' heigth in m
33%6  tri     =   tropopause heigth in m
34%7  tti     =   temperature
35%xmass  =   particle mass (for nspec species)
36
37% NB: If numpart < 1, only outheader is returned (output.outheader)
38
39% Written by Ignacio Pisso, January 2008
40% Modified
41% 5/2011: added output.outheader IP
42% 5/2013: modified to take more than one species into account
43% aded varargin for more than 1 species. Defalut is 1. 4 arguments required
44
45% Modified by Eivind G Waersted
46% 19/12-2014: made npoint (release point ID) be returned for each trajectory
47% 12/1-2015: vectorized the file-reading and made the function compute numpart directly from
48% the size of the file (NB: it now reads large files ~100 times faster than before)
49% 13/1-2015: handle the case that numpart < 1, and added closing of the file.
50% 2018-08-09 adapted for distribution with flexpart version 10.3 IP
51
52%t0 = cputime;
53
54if i<1
55    disp('i<0 ?')
56    disp(i);
57end
58
59% determine nspec
60nspec = 1;
61nVarargs = length(varargin);
62if nVarargs>0
63   nspec = varargin{1};
64end
65
66% file name
67if isnumeric(i)
68    dates=textread(strcat (directory,'/dates'));
69    %i=2
70    file=['partposit_' num2str(dates(i))];
71elseif ischar(i)
72    file=['partposit_' i];
73end
74file = strcat (directory,'/',file);
75
76%t1 = cputime;
77
78% if numpart not given, calculate numpart (total) and return it
79if ischar(numpart)
80    finfo = dir(file);
81    nbytes = finfo.bytes;
82    % assuming the file has the exact format written in partoutput.f90 of Flexpart v. 9.02,
83    % the size of the file is exactly 4 * [3 + (numpart+1)*(14+nspec)] bytes.
84    numpart = round((nbytes/4+3)/(14+nspec)-1);
85    output = numpart;
86    return
87end
88
89% if numpart is given, read the file and return the data
90
91[fid,message] = fopen(file,'r');
92if fid < 0
93    display(message)
94end
95
96%t2 = cputime;
97
98% record time record
99fread(fid,1,'int32');
100output.outheader=fread(fid,1,'int32');
101fread(fid,1,'int32');
102
103% if numpart < 1, return here
104if (numpart < 1)
105  return
106end
107
108% FLEXPART writes the output as:
109%           write(unitpartout) npoint(i),xlon,ylat,ztra1(i),
110%      +    itramem(i),topo,pvi,qvi,rhoi,hmixi,tri,tti,
111%      +    (xmass1(i,j),j=1,nspec)
112
113% as 4byte numbers, with an empty 4byte value before and afterwards
114
115% -> thus, the number of 4byte values read for each trajectory is 14+nspec
116nvals = 14 + nspec;
117
118% read data as int32 first, to get itramem and npoint
119dataAsInt = fread(fid,[nvals,numpart],'int32');
120fclose(fid);
121
122% reread file, this time as float32, to get the rest of the data
123fid = fopen(file,'r');
124fread(fid,3,'int32'); % skip header
125dataAsFloat = fread(fid,[nvals,numpart],'float32');
126fclose(fid);
127
128%t3 = cputime;
129
130% create the output structure
131output.npoint = transpose(dataAsInt(2,1:numpart));
132output.xyz = transpose(dataAsFloat(3:5,1:numpart));
133output.itramem = transpose(dataAsInt(6,1:numpart));
134output.vars = transpose(dataAsFloat(7:13,1:numpart));
135output.xmass = transpose(dataAsFloat(14:(13+nspec),1:numpart));
136
137%t4 = cputime;
138
139% Print how long time the function used on each part of the process
140%fprintf('Performance:\n %15s: %.2f s\n %15s: %.2f s\n %15s: %.2f s\n %15s: %.2f s\n','Inputs',t1-t0,'Allocations',t2-t1,'Reading',t3-t2,'Structure',t4-t3)
141
142return
143
hosted by ZAMG