startarchive2icosa — ICOSA 初始场转换链

路径

LMDZ.MARS\util\startarchive2icosa

所属目录/模块

util

目录定位

startarchive2icosa 是把 LMDZ.MARS 的 start_archive.nc（lon-lat 网格初始场）转换为 Dynamico ICOSA 正二十面体网格初始场的离线工具链。它解决的问题是：LMDZ GCM 的 restart 文件使用规则经纬网格，而 Dynamico 使用非结构化的正二十面体网格，两者不能直接互用。

整个工具链分三步：（1）start_archive2icosa 主程序通过 XIOS 把 lon-lat 场重映射到 ICOSA 网格并输出 *_prefinalize.nc；（2）icosa_finalize.bash 用 NCO 工具修补维度名、补充缺失变量；（3）rearrange_startphy 按经纬度匹配重新排列物理场的列顺序，使其与 ICOSA restart 文件一致。

README 明确标注此工具仍在开发中，用户需自行调整硬编码变量。

目录文件总览

定义的符号

start_archive2icosa.f90

rearrange_startphy.f90

依赖的模块

输入

start_archive2icosa 主程序

rearrange_startphy 程序

输出

start_archive2icosa 主程序

icosa_finalize.bash 后处理

核心逻辑

start_archive2icosa.f90

1. MPI+XIOS 初始化：MPI_INIT → xios_initialize 获取通信子 → xios_context_initialize("test")
2. 读地形文件维度：打开 surface_nc4.nc，从 zMOL 变量获取 src_topo_ni_glo/src_topo_nj_glo，按 MPI rank 沿纬度分割
3. 读 start_archive 维度：打开 start_archive_nc4.nc，从 temp 获取 src_ni_glo/src_nj_glo/src_nlev；从 tsoil 获取 src_nsoil/src_nslope；硬编码 src_nq=11
4. 读坐标和控制参数：从 start_archive 读 rlonv/rlatu（弧度→度）、ap/bp；从 startphy_icosa_ref.nc 读 controle（gravity=controle(7)、kappa=controle(9)）；preff 硬编码为 610 Pa
5. 读目标网格：从 start_icosa_ref.nc 读 lon_mesh/lat_mesh/bounds_*_mesh，按 MPI rank 分配列
6. XIOS 网格/轴定义：设置 src_domain_regular（rectilinear）、src_topo_domain_regular、src_domain_regular_clean（去掉经度最后一列的冗余点）、dst_domain_unstructured（ICOSA）；定义 lev/lev_p1/nq/soil_layers/subslope 轴
7. 时间循环（单步 src_nt=1）：
   • 地形：xios_recv_field("zMOL") → xios_send_field("topo", zMOL*gravity*1000)（km→m²/s²）
   • 地表压力：xios_recv_field("src_ps") → 计算层界面压力 src_ap(l)+src_bp(l)*ps → xios_send_field("ps_clean")
   • 地表温度/发射率：接收 subslope 维度 → 去掉最后一列 → 发送
   • q2（湍流动能）：接收 src_q2surf 和 src_q2atm，拼成 lev_p1 维度
   • 温度→theta_rhodz：T * (ΔP/g) / ((P_mean/preff)^κ)
   • 风场：u/v 直接去掉最后一列发送
   • 11 个 tracer（循环 i=1,src_nq）：按硬编码的 if/elseif 分支映射到 co2_clean、dust_number_clean、dust_mass_clean、ccn_number_clean、ccn_mass_clean、h2o_ice_clean、h2o_vap_clean、stormdust_*、topdust_*，同时处理地表 tracer 值
   • 土壤温度/热惯量：tsoil（含 subslope）和 inertiedat
   • SSO 参数：ZMEA/ZSTD/ZSIG/ZGAM/ZTHE → *_clean
   • 反照率/粗糙度：albedodat/z0 → *_clean
8. 收尾：rank 0 向 startfi_prefinalize.nc 写 controle，向 start_icosa_prefinalize.nc 写 iteration=0 → xios_finalize → MPI_FINALIZE

rearrange_startphy.f90

1. 从 startphy_icosa_ref.nc 读所有列的 latitude/longitude
2. 从 startfi.nc 读 latitude/longitude
3. 逐列匹配：对 startfi.nc 的第 i 列，在 ref 中找 lon/lat 相对误差 ≤10⁻⁵ 的第 j 列，记录 cell_index(i)=j
4. 按维度分类遍历所有变量：1D（physical_points）、2D（nvertex×physical_points、subslope×physical_points、lev_p1×physical_points、subsurface_layers×physical_points）、3D（subsurface_layers×subslope×physical_points），按 cell_index 重排并原地写回

伪代码

program start_archive2icosa:
  MPI_INIT + xios_initialize
  read surface_nc4.nc: get topo grid (zMOL)
  read start_archive_nc4.nc: get atm grid (temp, tsoil), ap/bp, lon/lat
  read startphy_icosa_ref.nc: get controle (gravity, kappa)
  read start_icosa_ref.nc: get ICOSA mesh (lon_mesh, lat_mesh, bounds)
  set XIOS domains (rectilinear src + unstructured dst) and axes
  time loop (1 step):
    recv topo → send topo (convert km→m²/s²)
    recv ps → compute layer pressures → send ps_clean
    recv tsurf/emis → strip last lon → send *_clean
    recv q2surf+q2atm → stack → send q2_clean
    recv temp → compute theta_rhodz → send theta_rhodz_clean
    recv u/v → strip last lon → send u_clean/v_clean
    for i=1..11: recv tracer[i] → send to named field (co2/dust/ccn/h2o/stormdust/topdust)
    recv tsoil/inertiedat → send tsoil_clean/inertiedat_clean
    recv SSO (ZMEA/ZSTD/ZSIG/ZGAM/ZTHE) → send *_clean
    recv albedodat/z0 → send *_clean
  rank 0: write controle to startfi_prefinalize, iteration to start_icosa_prefinalize
  finalize XIOS + MPI

program rearrange_startphy:
  read ref lon/lat from startphy_icosa_ref.nc
  read current lon/lat from startfi.nc
  for each cell i: find matching ref cell j by lon/lat (relative error ≤ 1e-5)
  for each variable: rearrange by cell_index and write back in place

参与的主题流程

写法特点

• 硬编码文件名：5 个输入/输出文件名（start_archive_nc4.nc、surface_nc4.nc 等）在源码中以 PARAMETER 硬编码，改文件名需改源码
• 硬编码 tracer 映射：11 个 tracer 的 if/elseif 分支（行 432–498）将索引 1–11 硬映射到 co2/dust_number/dust_mass/ccn_number/ccn_mass/h2o_ice/h2o_vap/stormdust_number/stormdust_mass/topdust_number/topdust_mass，与 tracer.def_no_start 的顺序一致
• 硬编码常数：src_nq=11（行 26）、preff=610（行 216，注释说明 restart 文件中缺此值）
• 去掉最后一列经度：src_domain_regular_clean 维度为 ni_glo-1，所有 *_clean 字段取 src_field(1:src_ni-1,...) 以去除 lon-lat 网格的周期性冗余点
• theta_rhodz 计算：行 406–410 将温度转为 theta * rho * dz 形式：T * (ΔP/g) / ((P_mean/preff)^κ)
• 单步时间循环：src_nt=1，XIOS 只做一次场传输
• MPI 沿纬度分割：源和目标网格都按 latitude 方向均分给 MPI 进程
• 列重排暴力匹配：rearrange_startphy 用 O(n²) 双重循环逐列匹配 lon/lat，相对误差阈值 10⁻⁵

复现要点

• 必须先用 Dynamico 跑一次（startphy_file=.false.）生成 start.nc/startfi.nc，重命名为 start_icosa_ref.nc/startphy_icosa_ref.nc
• start_archive.nc 必须含 SSO 参数（需用修改版 start2archive_SSO 生成）
• 输入 NetCDF 必须转为 NetCDF4 格式（nccopy -k 4）
• compile 脚本假设 XIOS 位于 ../../../XIOS/，需根据实际环境调整
• icosa_finalize.bash 中有多处硬编码路径（../trunk/XIOS/arch.env、restartfi.nc），需自行适配
• rearrange_startphy 的 lon/lat 匹配要求两者都以弧度存储（源码行 84 做了 deg→rad 转换）
• controle(3)（day_ini）被强制设为 0（行 208），因为 ref 文件是 restart
• 使用生成的 start 文件运行 Dynamico 时，run_icosa.def 需加 etat0_start_file_colocated = true

待确认

• src_nq=11 硬编码是否与所有配置的 tracer 数一致。当前 tracer.def_no_start 正好 11 个，但其他 tracer 配置可能不同
• preff=610 硬编码是否在所有 Dynamico 版本中一致（注释说 restart 文件中缺此值）
• rearrange_startphy 的 O(n²) 匹配在大网格上是否可接受（README 提到高分辨率地形插值需要大量内存）
• icosa_finalize.bash 中 ncks -A -v soildepth restartfi.nc 等操作假设 restartfi.nc 存在于当前目录，来源和时机未说明
• start_archive2icosa 中 src_topo_iend=src_topo_ibegin+src_ni-1（行 101）使用了 src_ni（atm 维度）而非 src_topo_ni，两者在当前逻辑下相等（都等于 src_ni_glo），但若源和地形分辨率不同会出错

复现风险

• 文件名全部硬编码，任何重命名都会导致程序失败
• tracer 映射硬编码为 11 个且顺序固定，与 initracer.F90 的字符串匹配无关——如果 tracer 配置或顺序变化，程序不会报错但会把错误的物理量映射到错误的输出字段
• icosa_finalize.bash 中 ncks -A 从 restartfi.nc 补充变量，如果该文件不存在或维度不匹配会静默失败
• 高分辨率地形插值需要大量内存（README 建议 2 节点 56 核但只用 28 进程）
• README 标注"currently under development"，用户需自行检查和调整

相关页面

• util/index.md — util 后处理工具总览
• zrecast.md — 同属 util 的垂直坐标重投影工具