1 | function 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 |
|
---|
54 | if i<1
|
---|
55 | disp('i<0 ?')
|
---|
56 | disp(i);
|
---|
57 | end
|
---|
58 |
|
---|
59 | % determine nspec
|
---|
60 | nspec = 1;
|
---|
61 | nVarargs = length(varargin);
|
---|
62 | if nVarargs>0
|
---|
63 | nspec = varargin{1};
|
---|
64 | end
|
---|
65 |
|
---|
66 | % file name
|
---|
67 | if isnumeric(i)
|
---|
68 | dates=textread(strcat (directory,'/dates'));
|
---|
69 | %i=2
|
---|
70 | file=['partposit_' num2str(dates(i))];
|
---|
71 | elseif ischar(i)
|
---|
72 | file=['partposit_' i];
|
---|
73 | end
|
---|
74 | file = strcat (directory,'/',file);
|
---|
75 |
|
---|
76 | %t1 = cputime;
|
---|
77 |
|
---|
78 | % if numpart not given, calculate numpart (total) and return it
|
---|
79 | if 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
|
---|
87 | end
|
---|
88 |
|
---|
89 | % if numpart is given, read the file and return the data
|
---|
90 |
|
---|
91 | [fid,message] = fopen(file,'r');
|
---|
92 | if fid < 0
|
---|
93 | display(message)
|
---|
94 | end
|
---|
95 |
|
---|
96 | %t2 = cputime;
|
---|
97 |
|
---|
98 | % record time record
|
---|
99 | fread(fid,1,'int32');
|
---|
100 | output.outheader=fread(fid,1,'int32');
|
---|
101 | fread(fid,1,'int32');
|
---|
102 |
|
---|
103 | % if numpart < 1, return here
|
---|
104 | if (numpart < 1)
|
---|
105 | return
|
---|
106 | end
|
---|
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
|
---|
116 | nvals = 14 + nspec;
|
---|
117 |
|
---|
118 | % read data as int32 first, to get itramem and npoint
|
---|
119 | dataAsInt = fread(fid,[nvals,numpart],'int32');
|
---|
120 | fclose(fid);
|
---|
121 |
|
---|
122 | % reread file, this time as float32, to get the rest of the data
|
---|
123 | fid = fopen(file,'r');
|
---|
124 | fread(fid,3,'int32'); % skip header
|
---|
125 | dataAsFloat = fread(fid,[nvals,numpart],'float32');
|
---|
126 | fclose(fid);
|
---|
127 |
|
---|
128 | %t3 = cputime;
|
---|
129 |
|
---|
130 | % create the output structure
|
---|
131 | output.npoint = transpose(dataAsInt(2,1:numpart));
|
---|
132 | output.xyz = transpose(dataAsFloat(3:5,1:numpart));
|
---|
133 | output.itramem = transpose(dataAsInt(6,1:numpart));
|
---|
134 | output.vars = transpose(dataAsFloat(7:13,1:numpart));
|
---|
135 | output.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 |
|
---|
142 | return
|
---|
143 |
|
---|