简介Flexpart 10.4 的源码资源包flexpart.tar.gz面向在 Ubuntu 18.04 下开展大气粒子轨迹模拟、需要从源码编译安装模型的研究者与开发者尤其适合有一定 Linux 基础、希望理解依赖库构建逻辑的读者。整个压缩包约 328.8 MB共 718 个文件除 Fortran 源代码、makefile 和 shell 编译脚本外还包含大量 HTML 说明文档、PNG 示意图、RST/TXT 技术文档、CSV 样例数据、f90 源码、控制文件与运行日志以及用于测试的 GRIB 气象数据文件与 Python 辅助脚本目录层级完整便于对照查阅。针对源码编译前最棘手的依赖准备问题包内相关说明覆盖 libemos、HDF5、NetCDF、eccodes、grib_api、zlib、szip、jasper 等核心依赖库的作用、版本匹配与编译顺序要点能帮助读者理清 Flexpart 在 Ubuntu 18.04 环境中的完整编译链路也方便后续针对特定数据格式进行二次开发。已有 608 人学习适合需要结合系统环境完成 Flexpart 安装、调试与结果验证的研究人员。1. flexpart 是什么一个能逆向追污染源的拉格朗日粒子扩散模型flexpart 这个包在我手里是当“后悔药”用的。某模拟项目X里监测点凌晨出现一个高浓度峰值排放源清单怎么都对不上我把这套拉格朗日粒子扩散模型拿到 Ubuntu 上编译跑通用反向追踪模式把 72 小时气团轨迹追回去才定位到问题源头。它能解决的不是“风往哪吹”这种大路问题而是更细的一层给定浓度观测或排放假设反推气团从哪来、向哪去、沿途经过了哪些高度层、有没有可能被抬升到边界层以上。适合做大气污染溯源、火山灰扩散、放射性物质传输评估的从业者新手按本文也能在 Ubuntu 上把一个最小案例跑起来。全文按我实际拆包顺序展开先讲编译再讲配置最后讲跑结果和踩坑。2. 在 Ubuntu 上编译 flexpart先把三类依赖装齐再动 makefileflexpart 是 Fortran 写的整个模型本质上是一个可执行文件加一堆输入文件。编译这步卡住的人不少原因不是源码有问题而是依赖没备齐。我先后在几台 Ubuntu 20.04 和 22.04 上编过结论一致依赖装对了make 一次过依赖缺了报错信息又短又玄学。2.1 依赖选型gfortran、NetCDF、eccodes 三件套flexpart 的编译需要三样东西Fortran 编译器、NetCDF 库带 Fortran 接口、以及读写气象 GRIB 数据的底层库。下面这个表格是我实际使用的组合按“够用且不容易踩坑”的标准选的。组件选型说明编译器gfortran开源、apt 直接装ifort 也能编但没必要NetCDFlibnetcdf-dev libnetcdf-fortran-dev必须带 Fortran 头文件 netcdff.modGRIB 读取libeccodes-dev新版库老的 grib_api 在新系统上容易编不过辅助make, openblasmake 必有openblas 可选但建议装很多人在 NetCDF 上翻车只装了 libnetcdf-dev没装 libnetcdf-fortran-dev结果 make 时找不到 netcdff.mod。这两个包的区别前者是 C 接口后者才是 Fortran 接口。flexpart 的源码里大量调用nf_open、nf_create这类函数没有 Fortran 接口直接编译失败。安装命令如下sudo apt update sudo apt install -y gfortran make \ libnetcdf-dev libnetcdf-fortran-dev \ libeccodes-dev libopenblas-dev说明一下每个包的作用gfortran 提供编译器libnetcdf-fortran-dev 提供 netcdff.mod 和 Fortran 链接库libeccodes-dev 是 ECMWF 编码库的 Debian 包flexpart 读 GRIB 气象数据时靠它解析openblas 不是必须的但某些版本编译时链接数学库装上能省去报错。如果你的 Ubuntu 是 18.04可能没有 libeccodes-dev 这个包名需要从源码编译 eccodes或者换用 grib_api但我在新版本上没遇到这种情况。2.2 编译操作解压、改 makefile 三处变量、make拿到 flexpart.tar.gz 后先解压看目录结构。这一步不要跳确认包内包含 src、tools、docs 目录再开始改配置。tar -xzf flexpart.tar.gz cd flexpart ls -la解压后进入源码根目录会看到一个 makefile 或 makefile.mk。flexpart 的 makefile 逻辑很直接顶部定义编译器、NetCDF 路径、eccodes 路径然后按编译选项生成可执行文件。我一般只改下面这几个变量# 编译器 FC gfortran # NetCDF 安装根目录 NETCDF /usr # eccodes 安装根目录 ECCODES /usr # 编译选项 FFLAGS -O2 -fdefault-real-8 -ffpe-trapinvalid,zero,overflow LDFLAGS -L$(NETCDF)/lib -lnetcdff -L$(ECCODES)/lib -leccodes这里有个关键参数-fdefault-real-8。flexpart 的源码假设默认 real 是 8 字节不加这个选项所有浮点运算都会按 4 字节处理结果会出现大量 NaN而且很难排查。-ffpe-trap是让浮点异常直接暴露出来编译时不加也行运行时能帮你更快发现数据问题。NetCDF 和 eccodes 的路径在 Ubuntu 上一般就是/usr因为 apt 装的库都在 /usr/lib 和 /usr/include 里。如果你是自己编译安装到 /opt/local 这类目录要把路径改成实际的安装位置否则链接时找不到库文件。改完之后直接 makemake -j4-j4 是用 4 个并行任务编译速度快一些。第一次编译如果机器配置低建议先不加 -j避免内存不够。编译过程会输出大量 Fortran 编译信息看到 src/flexpart 这个可执行文件生成就算成了。2.3 验证编译结果先跑一个默认用例编译通过不代表能跑通。flexpart 编译完我习惯先不急着配自己的数据而是看它有没有内置测试用例。很多发行版源码树里带一个小规模示例数据放在 tests 或 examples 目录下。ls src/flexpart ./src/flexpart不带参数直接运行正常情况下会读取当前目录下的 pathnames 文件。如果当前目录没有这个文件它会报错退出提示缺少 pathnames。这个报错其实是个好消息说明可执行文件能正常启动、Fortran 运行时没问题。如果连这个都出不来比如提示找不到共享库那就是链接路径没配对用ldd src/flexpart查一下缺哪个库。验证完这一步编译阶段就算结束了。接下来真正花时间的是配置运行环境。3. 跑通第一个模拟pathnames 与 COMMAND 的配置细节编译只是热身flexpart 的使用门槛主要在配置。它不像普通软件有图形界面所有控制全靠几个文本文件。我最早接触时对着文档折腾了两天才弄清楚每个文件的职责这里把最关键的几个讲透。3.1 运行目录文件清单到底需要准备哪些文件flexpart 运行时当前目录下必须存在一组固定文件。不同版本文件数量略有差异以最常见的源码树为例至少需要下面这些文件作用pathnames告诉你各输入输出数据放在哪个目录COMMAND控制模拟主参数如起止时间、粒子数、输出间隔OUTGRID定义浓度输出的网格范围与分辨率RELEASES定义释放源的经纬度、高度、时间窗口available物种与沉降参数数据库注意运行时只认当前目录下的这几个名字所以每次做新案例我一般复制一个干净目录再改里面的文件而不是直接在源码目录里跑避免污染原始配置。常见做法是这样组织目录mkdir -p /data/flexpart/case1 cd /data/flexpart/case1 cp /data/flexpart/available ./available然后把 pathnames、COMMAND、OUTGRID、RELEASES 这四个文件按需求改好再执行 flexpart。源码树里会带模板先把模板复制过来改比自己新建文件少踩很多格式坑。3.2 pathnames路径文件里每一行的含义pathnames 是 flexpart 第一个读取的文件它按固定行序告诉模型各数据目录在哪。格式很简单但有一个隐蔽的坑Fortran 读取时对行列位置敏感路径后面不能有多余空格目录层级也不能错。我常用的一个 pathnames 长这样/data/flexpart/input/met /data/flexpart/input/geo /data/flexpart/input/oh /data/flexpart/output第一行是气象数据的根目录第二行是地理数据目录第三行是 OH 自由基浓度场数据目录第四行是输出目录。如果你的模拟不涉及 OH 化学第三行可以指向一个空目录但这一行不能省略。行末如果多了一个空格Fortran 的固定格式读取会把它当成路径的一部分运行时直接报错找不到目录。这个文件改完可以用一个很简单的命令检查行尾cat -A pathnames每行结尾应该是干净的$如果看到^M或者空格$说明文件带了 Windows 换行或多余空格。这个细节看起来不起眼但能省掉你半小时排查时间。3.3 COMMAND 和 RELEASES粒子数、释放时间、输出间隔COMMAND 文件用 Fortran namelist 格式写成控制模拟的全局行为。我一般重点关注这几个参数COMMAND PARTICLE_NUMBER 10000 LDTEFF .TRUE. LSUBGRID .TRUE. LDIRD 1 MAXAGE 72.0 LOUTSTEP 3600.0 LOUTAVER .TRUE. /PARTICLE_NUMBER 是释放的总粒子数数值越大浓度场越平滑但运行时间和输出文件体积也越大。LDIRD 是方向控制1 表示正向追踪-1 表示反向追踪做溯源时用 -1。MAXAGE 是粒子最大存活时间单位小时比如 72.0 表示追踪三天。LOUTSTEP 是浓度输出时间间隔单位秒3600 就是每小时输出一帧。LOUTAVER 决定输出的是瞬时值还是时间平均值做浓度图一般用 .TRUE.。RELEASES 文件定义释放源格式比 COMMAND 更固定。我常用的写法是一个释放组包含时间窗口和空间范围1 2019-06-01 00:00 2019-06-01 02:00 0.0 72.0 1.0E06 0.0 0.0 113.0 30.0 0.0 100.0第一行是释放组编号第二行和第三行是释放开始和结束时间第四行是释放持续时间、粒子寿命上限、总释放质量、一个标志位第五行是释放中心经度、纬度、最低高度、最高高度单位分别是度和米。这里的经纬度我写的是某模拟项目X里用过的测试点替换成你研究区域的位置即可。每次改完这几个文件我都会做一次“最小配置验证”先只放 1000 个粒子、输出间隔拉大跑通流程后再把粒子数和输出频率加上去。这个习惯帮我避免了很多“配置错了但跑了几小时才发现”的尴尬。4. 运行与输出从粒子轨迹到浓度网格结果怎么读配置写完终于可以真正运行 flexpart。这一章讲输入数据怎么对齐、运行日志怎么看、以及输出文件如何解读。4.1 气象数据与地理数据先把输入对齐flexpart 不是自带气象数据的它需要你提供气象驱动场。常见做法是去全球再分析数据集下载 GRIB 格式的气象场然后放到 pathnames 指定的目录里。下载时要注意两个事一是时间范围要覆盖模拟窗口二是变量要完整至少包含风速、温度、湿度、气压和边界层高度这些核心场。下载完先检查文件再放入目录。我一般用 eccodes 提供的工具看一眼grib_ls /data/flexpart/input/met/20190601.gribgrib_ls 会列出 GRIB 文件里的所有场包括变量名、层数、时间。这里重点看数据的 level 类型flexpart 需要多层三维场如果下的是单层场运行时会报“无法插值”之类的错误。地理数据方面flexpart 需要地表高度、土地利用类型等文件这些一般在包内自带生成工具或者用模型配套脚本从公开地表数据集转换。把气象数据放好后再检查路径ls /data/flexpart/input/met/ ls /data/flexpart/input/geo/确认目录里确实有文件不是空目录。我遇到过好多次pathnames 指对了但目录是空的flexpart 一句话不报就退出看起来像程序崩了。所以养成习惯在运行之前先手动确认每个目录非空。4.2 启动运行命令行入口与运行日志所有输入就位后在运行目录里执行可执行文件cd /data/flexpart/case1 ../src/flexpart运行时终端会滚动大量文本包括正在读取哪个 GRIB 文件、插值进度、粒子释放信息。两个 Linux 小技巧可以帮上忙../src/flexpart run.log 21 tail -f run.log把输出重定向到日志文件再用 tail 实时看进度这样即使终端断开日志也保留下来。运行时间取决于粒子数、模拟时长和气象数据覆盖范围小案例几分钟大案例几小时。如果运行几分钟后日志不动了按 Ctrl-C 中断检查日志最后几行多半是路径问题或数据时间不匹配。4.3 输出文件解读partposit、grid_conc 与 header运行结束后输出目录里会出现几类文件。最简单的是 header记录模拟基本信息比如坐标系参数、网格定义。核心输出是两类一类是粒子位置文件 partposit记录每个粒子随时间变化的坐标一类是浓度网格文件 grid_conc按你设置的网格输出浓度场。partposit 是二进制格式不能直接打开看。flexpart 源码树里带了转换工具常见做法是到 tools 目录编译 readpartpositcd flexpart/tools gfortran readpartposit.f90 -o readpartposit cd /data/flexpart/output ../flexpart/tools/readpartposit partposit partposit.txt转换出来的文本一行一个粒子列包含时间、经纬度、高度等信息。用 head 看一眼格式head -20 partposit.txtgrid_conc 文件如果编译时 NetCDF 接口正常会输出为 NetCDF 格式直接用 Python 的 xarray 读取画图很方便import xarray as xr import matplotlib.pyplot as plt ds xr.open_dataset(grid_conc_20190601_00.nc) conc ds[conc].sum(dimtime).squeeze() fig, ax plt.subplots(figsize(8, 6)) mesh conc.plot(axax, cmapYlOrRd) ax.set_title(Column-integrated concentration (ng/m3)) plt.savefig(concentration.png, dpi150)这段代码把整个模拟时段的浓度按时间维度求和得到一张累计浓度分布图。对于污染溯源场景这张图能直接告诉你高浓度覆盖了哪些区域再叠加监测点位就能对照。5. 避坑清单编译与运行阶段的五条血泪经验这一章直接写我实际踩过的坑每条按现象、原因、解决三个步骤展开照着排查能省不少时间。5.1 编译报错找不到 netcdff.mod现象执行 make 后gfortran 报错Error: Cannot open file netcdff.mod或者提示cant find -lnetcdff。原因只装了 libnetcdf-dev没装 libnetcdf-fortran-dev。前者是 C 接口的库不包含 Fortran module 文件。解决sudo apt install libnetcdf-fortran-dev装完后检查 /usr/include 或 /usr/lib/gcc 相关目录下是否存在 netcdff.mod。如果存在但 makefile 里没指定头文件路径在 FFLAGS 里加一行FFLAGS -I/usr/include从那以后我每次新机器编译前都先执行apt search netcdf看一眼 Fortran 接口包是否在列再开始动 makefile。5.2 运行秒退但无任何报错信息现象在运行目录执行./flexpart终端只输出一两行进程就退出没有明确的 error 信息。原因最常见的两个。一是 pathnames 里某个路径末尾多了一个空格Fortran 固定格式读取把空格当成了路径的一部分二是某个数据目录是空的flexpart 检查不到文件直接退出。解决先cat -A pathnames确认行尾干净再逐个ls检查目录内容最后在 makefile 里加-g -fbacktrace重新编译一遍。flexpart 的日志本身比较简略编译时加回溯选项能定位到具体是哪个文件读取失败。这个坑我连续遇到两次后来养成了“先 cat -A再 ls -l”的固定套路。5.3 模拟正常结束但浓度场全零现象粒子数设了 10 万跑了几小时日志显示正常释放、正常结束但打开 grid_conc 全是 0 或 NaN。原因OUTGRID 文件定义的网格范围没有覆盖释放源或者 COMMAND 里 LOUTSTEP 设得太大以至于每个输出间隔里粒子恰好还没进入网格。还有一种情况是释放质量设了 0RELEASES 里那一行的质量参数写成了 0.0粒子虽然没有密度但浓度算出来自然是零。解决先检查 OUTGRID 的经纬度范围确认释放点落在网格内再确认 RELEASES 里的质量参数大于 0最后把 LOUTSTEP 缩短到 600 秒做一次小规模测试。做污染溯源时网格范围宁大勿小哪怕分辨率低一点先把浓度场跑出来再局部加密。5.4 反向追踪结果在原地打转现象设置 LDIRD -1 做溯源跑出来的粒子轨迹不往外走就围着释放点附近转圈。原因反向模式对气象数据的起始时间要求与正向不同。反向追踪时粒子从观测时刻出发往回走如果气象数据没有覆盖观测时刻之前的时间窗或者 RELEASES 里的释放时间写得比气象数据还早粒子在初始时刻就找不到有效风场只能原地乱飘。解决把 RELEASES 的释放时间设为观测时刻之前至少一个 MAXAGE 小时同时确认气象数据的起始时间比释放时间还早。比如观测是 6 月 1 日 02 时MAXAGE 设 72 小时气象数据从 5 月 30 日 00 时就要有不能只覆盖 6 月 1 日当天。另外反向模式下粒子数建议加大到至少 5 万不然统计噪声会被近似成原地不动。5.5 GRIB 数据读取报错“key not found”现象程序运行到读取 GRIB 文件时报错提示某个 key 找不到。原因气象数据的格式或变量命名与 flexpart 预期不一致。flexpart 对不同来源的气象数据有严格约定它要求特定的字段名、层数、网格类型。下载数据时选了错误的参数组合就会这样。解决先看报错里提到的 key 是什么对比数据里的变量名和包内文档中列出的变量清单。最常见的是用了单层场而非多层场或网格类型是 regular 而非 reduced。另外检查 eccodes 版本新版库对 GRIB2 支持更完整旧版容易出现 key 缺失。我一般会先跑一个只用 3 小时数据的极短模拟确认数据读取正常再跑正式案例。6. 进阶技巧批量跑情景模拟与 partposit 可视化配置和排错都稳了之后日常使用中我最常用的两个进阶操作批量改释放点跑多情景以及把 partposit 轨迹画出来。批量做情景时不需要手动一个一个改 RELEASES用 bash 循环加 sed 替换就可以。RELEASES 文件里的经纬度可以用占位符代替比如把经度写成变量 LAT然后循环替换for i in 1 2 3 4 5; do mkdir -p run_case_$i sed s/LAT/$i/g template_RELEASES run_case_$i/RELEASES cd run_case_$i cp ../COMMAND ../pathnames ../OUTGRID ../available . ../src/flexpart run.log 21 cd .. done每次循环生成一个独立目录各自放一套配置跑完不互相干扰。这里 sed 命令把模板里的 LAT 替换成实际数值适合快速生成多释放点实验。需要注意 COMMAND 里输出的文件名如果指定了绝对路径每次循环要一并改掉否则结果会覆盖。partposit 转成文本后画轨迹最直接。readpartposit 转换出的列一般依次是时间、经度、纬度、高度、序号等。画三维轨迹时用 matplotlibimport matplotlib.pyplot as plt import pandas as pd df pd.read_csv(partposit.txt, delim_whitespaceTrue, headerNone, names[time, lon, lat, z, num]) sample df.sample(2000, random_state42) fig plt.figure(figsize(10, 8)) ax fig.add_subplot(111, projection3d) ax.plot(sample[lon], sample[lat], sample[z], lw0.3, colorsteelblue) ax.set_xlabel(Longitude) ax.set_ylabel(Latitude) ax.set_zlabel(Height (m)) plt.savefig(trajectory.png, dpi150)粒子数多的时候画 10 万条线会很密集所以先随机抽 2000 个粒子再画轨迹趋势不受影响。如果 partposit 列数和你手上版本不同先head -5看实际列数再改 names 列表。关于 partposit 还有一个容易被忽略的点不同版本的转换工具输出的列顺序不完全一样。我用过一次某个历史版本时间列和高度列是反的画出来的轨迹高度变化看起来非常离谱。后来我固定了“先 head 看表头再画图”的流程从那以后每次跑完模拟都强制走一遍先 cat -A 检查路径文件再 head 检查转换结果的列数最后才进画图和浓度分析。希望帮到你。本文还有配套的精品资源点击获取