read_dust_scenario_mod.F90

快速理解

它做什么： 从 datadir 读 dust_*.nc NetCDF 场景文件，插值到 GCM 网格返回参考气压处可见光尘埃柱厚度。

基本过程： 首次打开文件缓存经纬/时间/tau → 按 zday 时间插值 + 经纬双线性插值 → 输出 tau_pref_scenario。

关键结果： tau_pref_scenario(ngrid)，供 compute_dtau 生成注入目标和 dust_rad_adjust 生成辐射缩放目标。

路径

LMDZ.MARS\libf\phymars\read_dust_scenario_mod.F90

所属目录/模块

libf\phymars

文件定位

read_dust_scenario_mod.F90 定义尘埃场景读取模块 read_dust_scenario_mod（第 1 行）和例程 read_dust_scenario（第 7 行）。它在第一次调用时从 datadir 下打开由 iaervar 选择的 dust_*.nc NetCDF 场景文件，缓存经纬度、时间轴和三维柱光学厚度场；之后每次调用按当前 GCM 网格经纬度和 zday 做时间线性插值、经纬双线性插值，返回参考气压 odpref 处的可见光尘埃柱厚度 tau_pref_scenario(ngrid)。

该例程是尘埃观测场景约束的底层读表器：compute_dtau_mod 用它生成地表注入目标厚度，dust_rad_adjust_mod 用它生成仅辐射缩放目标厚度。它本身不修改 tracer 或辐射场，只把外部尘埃场景文件转换为当前物理网格上的可见光柱厚度。

定义的符号

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• 首次调用读取外部 NetCDF 文件；文件路径为 datadir/filename，datadir 可由 callphys.def 配置路径覆盖。
• firstcall、lat、lon、time、tautes、timelen/lonlen/latlen、filename、IRscenario 都是 SAVE 状态，并在第 56-59 行声明为 OpenMP THREADPRIVATE。
• 非 MESOSCALE 编译路径下，只有 OMP MASTER 读取 NetCDF，再通过 bcast 分发给其他并行进程。
• 错误处理会写标准输出并调用 abort_physic，包括找不到文件、缺少变量、时间长度不等于 year_day、经纬索引越界、负光学厚度和插值权重越界。
• 本例程不写 wiki、restart 或诊断文件；运行时副作用主要是读外部场景文件和填充线程私有缓存。

核心逻辑

1. 首次调用选择并读取场景文件：按 iaervar 映射 dust_tes.nc、dust_cold.nc、dust_warm.nc、dust_clim.nc、dust_MYxx.nc 或 dust_tes_MYxx.nc（第 70-94 行），再从 datadir 打开文件。
2. 读取维度和坐标：识别 Time 或 time 维度，读取 latitude、longitude 和时间轴长度；非 MESOSCALE 路径广播维度后分配缓存数组。
3. 读取尘埃柱厚度变量：优先读 dustop 并标记为可见光场景；若没有 dustop，依次读 cdod 或 cdodext 并标记为红外场景。cdodext 会先除以 1.3 转为红外吸收厚度。
4. 按需要重写时间轴：当 dustinjection>0 或 dustscaling_mode==2 时，要求 timelen==year_day，然后把 time(iloop) 写成 0..timelen-1，保证每日一张图的场景按整数 sol 索引。
5. 逐 GCM 格点找四邻点：把 geometry_mod 的弧度经纬度转为度，纬度改成余纬；对略小于 -180 的经度做容差修正；定位上下纬度和左右经度，跨 180 经线时用 xsup=1/xinf=lonlen。
6. 先时间插值、再水平双线性插值：对四个邻点分别在 tinf/tsup 两个时间层间线性插值，跨年时把分母写成 time(tsup)+(669-time(tinf))；随后用 dlat/dlon 做双线性插值生成 tau。
7. 单位/历史方案修正：若 IRscenario 为真，将 tau 乘 IRtoVIScoef(ig) 转为可见光厚度；旧 TES 文件再乘 1.3*(odpref/700.)，swrtype==1 时再乘 0.825/1.3。

伪代码

realday = mod(zday, 669)

if firstcall:
    filename = select_by_iaervar(iaervar)
    open datadir/filename
    read lon/lat/time dimensions
    allocate tautes(lonlen,latlen,timelen), lon, lat, time

    if variable dustop exists:
        read tautes
        IRscenario = false
    else if variable cdod exists:
        read tautes
        IRscenario = true
    else:
        read cdodext into tautes
        tautes = tautes / 1.3
        IRscenario = true

    read Time, latitude, longitude
    if dustinjection>0 or dustscaling_mode==2:
        require timelen == year_day
        time = [0, 1, ..., timelen-1]
    broadcast cached fields on non-MESOSCALE builds

for each GCM grid point ig:
    convert latitude/longitude from radians to degrees
    find enclosing scenario latitudes and longitudes
    find enclosing scenario times, with wraparound across sol 669
    interpolate tau at four horizontal neighbours in time
    bilinear-interpolate the four tau values in latitude/longitude
    if IRscenario:
        tau_pref_scenario(ig) = tau * IRtoVIScoef(ig)
    else:
        tau_pref_scenario(ig) = tau

if filename begins with "dust_tes":
    tau_pref_scenario *= 1.3 * (odpref / 700)
if swrtype == 1:
    tau_pref_scenario *= 0.825 / 1.3

参与的主题流程

写法特点

• 自由格式 Fortran 90；数据读取集中在 firstcall，后续调用复用缓存，避免每步重复打开 NetCDF。
• 输入场景变量名有兼容路径：dustop 表示可见光消光，cdod 表示红外吸收，cdodext 表示红外消光并先除以 1.3。
• 时间插值显式处理跨年首尾连接，但年长硬编码为 669.，而 dustinjection/dustscaling_mode==2 的逐日场景校验使用 year_day。
• 旧 TES 文件有两个历史修正：700 Pa 参考气压改到当前 odpref，以及红外吸收到红外消光的 1.3 因子。
• pplev 是接口输入但在当前源码可见逻辑中未直接使用；气压缩放只在旧 TES 文件尾部通过 odpref/700. 处理。

复现要点

• iaervar 必须选到实际存在于 datadir 的 dust_*.nc 文件；错误信息提示可用 callfis.def/callphys.def 的 datadir 指向数据目录。
• 场景文件必须至少有 Time 或 time 维度、latitude、longitude，以及 dustop、cdod、cdodext 之一。
• 若启用 dustinjection>0 或 dustscaling_mode==2，场景时间长度必须等于 planete_h::year_day，否则运行会中止。
• geometry_mod::longitude 进入本例程时以弧度给出，转成东经度后按场景经度表查找；跨 180 经线有专门分支。
• 红外场景必须提供合理的 IRtoVIScoef(ngrid)，否则输出的可见光柱厚度会系统性偏离。
• swrtype==1 的旧 Fouquart 路径会对最终厚度乘 0.825/1.3，复现实验时要确认短波辐射方案。

待确认

• realday=mod(zday,669.) 和跨年时间插值分母中的 669 是硬编码火星年长度；而逐日场景校验使用 year_day。如果 year_day 在某些配置中不是 669，需确认是否存在边界不一致风险。
• pplev(ngrid,nlayer+1) 当前未在可见逻辑中使用；保留在接口中可能是为与调用方签名或历史压力缩放路径兼容。
• abort_physic 在本文件中被直接调用但未显式 use；推断依赖工程中的全局接口或包含路径，复现编译环境时需保留原工程构建上下文。

相关页面

• dust_param_mod - 提供 odpref、dustscaling_mode 等尘埃配置状态
• compute_dtau_mod - 通过本例程读取目标厚度并转为 dustliftday
• dust_rad_adjust_mod - 通过本例程读取相邻日目标厚度并生成辐射缩放系数
• dust_scaling_mod - 在 dustscaling_mode==2 路径消费辐射缩放结果
• aeropacity_mod - 提供调用链中的尘埃柱厚度与 IR/VIS 转换背景
• datafile_mod - 管理 datadir 外部数据目录
• planete_h - 提供 year_day
• dust-cycle 主题
• callkeys_mod.F90 - 提供 iaervar、dustinjection、swrtype 配置入口。