LMDZ.COMMON Wiki

interp_horiz.F / iniinterp_horiz.F

路径

LMDZ.COMMON-6.3\LMDZ.COMMON\libf\dyn3d\interp_horiz.F
LMDZ.COMMON-6.3\LMDZ.COMMON\libf\dyn3d\iniinterp_horiz.F

文件定位

interp_horiz.F 和 iniinterp_horiz.F 是旧式 F77 水平重网格 helper。它们把一个 LMDZ scalar 网格上的 2D/3D intensive field 插值到另一个 LMDZ scalar 网格，方法是预先计算新旧网格 cell 的经纬边界交叠面积，再按交叠面积加权平均。

Mars 运行参与度：工具链条件经过。Mars datareadnc.F 用它把 surface.nc 的 360x180 surface 数据投到当前 GCM scalar 网格；lect_start_archive.F 用它把旧 start_archive.nc 中的 surface、soil、3D atmosphere 和 tracer 字段投到当前网格；newstart.F 的 MONS 数据加载也用它把外部数据投到 GCM 网格。普通时间推进主循环不调用这两个 helper。

定义的符号

符号 文件 行号 作用
interp_horiz interp_horiz.F 4 顶层水平插值入口，初始化交叠表后把 varo 累加到 varn。
iniinterp_horiz iniinterp_horiz.F 4 预计算新旧 scalar cell 边界、交叠列表 iik/jjk/ik/jk、交叠面积 intersec 和新网格面积 airen。

输入输出

interp_horiz(varo,varn,imo,jmo,imn,jmn,lm,rlonuo,rlatvo,rlonun,rlatvn) 的数组约定如下：

参数 方向 维度 说明
varo in (imo+1,jmo+1,lm) 旧 scalar 网格上的场。经向含极点行，纬向含周期补列。
varn out (imn+1,jmn+1,lm) 新 scalar 网格上的结果，例程内部先清零再累加。
imo,jmo in scalar 旧网格经纬参数；实际数组是 imo+1 by jmo+1。
imn,jmn in scalar 新网格经纬参数；实际数组是 imn+1 by jmn+1。
lm in scalar 垂直或字段层数；2D field 传 1，土壤或大气 3D 传对应层数。
rlonuo,rlonun in imo+1, imn+1 旧/新 scalar cell 经向边界，单位为弧度，最后一项通常是第一项加 2*pi。
rlatvo,rlatvn in jmo, jmn 旧/新 scalar cell 纬向内边界，单位为弧度；两极由例程内部补为 pi/2 和 -pi/2。

iniinterp_horiz 额外接收固定容量 kllm，并输出：

核心算法

1. 构造 cell 边界

iniinterp_horiz 先把输入的 scalar boundary 数组展开为旧网格 [a(i),b(i)] x [c(j),d(j)] 和新网格 [an(ii),bn(ii)] x [cn(jj),dn(jj)]。

经向边界有周期处理：

a(1) = -rlonuo(imo+1)
b(1) =  rlonuo(1)
a(i) =  rlonuo(i-1), b(i) = rlonuo(i) for i=2..imo+1

新网格同理。纬向边界则把北极和南极补为：

d(1) = pi/2
c(jmo+1) = -pi/2

这意味着调用方只需要提供 interior scalar-box latitude boundaries，极点由例程内部固定。

2. 计算新网格 cell 面积

每个新 cell 面积写为：

airen(ii,jj) = (bn(ii)-an(ii)) * (sin(dn(jj))-sin(cn(jj)))

源码注释称面积为 m2，但公式没有乘行星半径平方；实际是单位球面的 solid-angle 面积。由于插值使用 intersec/airen 比值，半径因子会相消。复现时不要再额外给 intersec 或 airen 单边乘半径。

3. 枚举新旧 cell 交叠

iniinterp_horiz 对每个新 cell 和旧 cell 判断纬向是否重叠、经向是否重叠。经向除正常交叠外，还检查旧 cell 平移 -2*pi 或 +2*pi 后的周期交叠，用于处理跨经度断点的 cell。

一旦重叠：

ktotal += 1
iik(ktotal)=ii, jjk(ktotal)=jj
ik(ktotal)=i,   jk(ktotal)=j
intersec = (overlap_lon_width) * (sin(overlap_north)-sin(overlap_south))

固定容量 kllm 在 interp_horiz.F 中设为 1400*200*10。源码没有检查 ktotal > kllm，因此极高分辨率或异常边界输入存在越界风险。

4. 按面积权重累加字段

interp_horiz 调 iniinterp_horiz 后先把 varn 全部清零，再遍历 k=1..ktotal：

varn(iik(k),jjk(k),l) += varo(ik(k),jk(k),l) * intersec(k) / airen(iik(k),jjk(k))

对每个 l=1..lm 独立执行同一水平映射。这个公式对 intensive field 做面积加权平均；若把 extensive total 当作 varo 传入，会改变物理含义。

5. 极点行平均

最后例程对新网格北极行 jj=1 和南极行 jj=jmn+1 做经向平均：

totn = sum_i varn(i,1,l)
tots = sum_i varn(i,jmn+1,l)
varn(:,1,l) = totn / (imn+1)
varn(:,jmn+1,l) = tots / (imn+1)

这保证极点只有一个物理值，即使数组中仍保留 imn+1 个经向位置。

Mars 调用链

调用方 用法 后续处理
datareadnc.F 读取 surface.nc 的 z0/albedo/thermal/zMOL 等 360x180 数据，补周期列后调 interp_horiz(zdataS,pfield,imd,jmd,iim,jjm,1,rlonud,rlatvd,rlonu,rlatv)。 手动恢复周期列；newstart 再用 gr_dyn_fi 把动力 scalar grid field 转到物理网格。
lect_start_archive.F 2D fields 对 phis/co2ice/tsurf/emis/tauscaling/totcloudfrac/watercap/peren_co2ice/albedo/ps 等按旧 archive 网格到新网格重投影。 多数 surface/physics 字段随后 gr_dyn_fi；ps 和 CO2 ice 后续再做总量缩放。
lect_start_archive.F soil fields 对 inertiedat/tsoil 等按 nsoilmx 或旧 soil 层数传入 lm，做逐层水平插值。 必要时先经 interp_line 做土壤垂直插值，再水平插值。
lect_start_archive.F atmosphere fields 对温度、q2、风和 tracer 先经 interp_vert 到当前垂直层，再用 interp_horiz 到当前水平网格。 风场插值后再 scal_wind，并乘 cu/cv 写回 covariant wind。
newstart.F MONS 数据 把 Hdnx 和 d21x 从 180x90 MONS 网格投到 MONS_Hdn/MONS_d21。 d21 在读入时从 g/cm2 乘 10 转为 kg/m2，再水平插值。

复现检查清单

  1. 确认旧/新经度边界数组单位为弧度，且最后一列是周期补列。
  2. 确认纬向边界数组不包含极点本身；北/南极由 iniinterp_horiz 补入。
  3. 确认 varo 和 varn 的第一维都是 im+1，并包含周期列。
  4. 对 2D field 传 lm=1；对 soil、atmosphere 或 tracer profile 传真实层数。
  5. 对需要守恒的 Mars 物理量，检查调用方是否在 interp_horiz 后另做总量缩放；lect_start_archive 对 ps 和 CO2 ice 有显式缩放。
  6. 若出现异常打印 ktotal 1 = ...，先检查边界数组顺序、kllm 容量和经度断点。

边界和风险

相关页面