filtreg_mod.F90 滤波初始化模块

源码范围：

LMDZ.COMMON-6.3\LMDZ.COMMON\libf\filtrez\filtreg_mod.F90

源码行数：559。

Mars 运行参与度：条件经过。Mars 3D GCM 和 Mars nogcm 在动力初始化后调用 inifilr；newstart、start2archive 等 Mars 工具链在构造或转换 start/archive 文件前也调用它，以便后续压力、几何或动力变量处理拥有同一套纬向滤波状态。

文件职责

filtreg_mod 是 filtrez 目录的滤波状态初始化模块。它不直接滤波字段，而是在 inifilr 中：

• 读取当前水平网格几何和 zoom/stretch 状态；
• 调用 inifgn 构造离散 Laplacian 特征值和特征函数；
• 计算 U/scalar 网格与 V/Z 网格的滤波纬带、起滤波模态和滤波系数；
• 分配并填充矩阵滤波路径使用的北半球、南半球矩阵；
• 在 CPP_PARA 且 use_filtre_fft 为真时初始化全局和本地 FFT 滤波 backend。

实际施加滤波的是 filtreg.F 和 dyn3dpar/filtreg_p.F；它们使用本模块保存的矩阵数组和 coefils.h 中的滤波边界状态。

符号清单

本文件没有 THREADPRIVATE 声明；并行路径通过模块全局数组和 coefils.h 状态共享初始化结果。

依赖和状态

USE 依赖

include 状态

调用关系

主要调用方

被调例程和后续使用位置

inifilr 调用 inifgn(eignvl)，由 filter-helpers 覆盖。矩阵路径后续由 filtreg-system 中的串行 filtreg.F 和并行 dyn3dpar/filtreg_p.F 使用。

初始化流程

1. 读取网格尺度并计算特征值

inifilr 首先设置 pi = 2*asin(1.)，把 xprimu(i) 拷贝到 dlonu(i)，调用 inifgn(eignvl) 计算离散 Laplacian 特征值。随后从 rlatu(j)-rlatu(j+1) 得到 dlatu(j)。

普通构建分支用显式循环求：

• dxmin = min(dlonu(:))
• dymin = min(dlatu(:))

CRAY 分支则用 ISMIN 找最小值位置。

2. 计算 colat0 和 lamdamax

默认：

colat0 = min(0.5, dymin / dxmin)

这表示规则网格从约 60 度纬度附近开始滤波；zoom 网格则在纬向长度和经向长度相近时开始向极区滤波。

特殊分支：

if (.not. fxyhypb .and. ysinus) then
  colat0 = 0.6
  alphax = 0.
endif

若 alphax == 1.，例程打印错误并 STOP。否则：

lamdamax = iim / (pi * colat0 * (1 - alphax))
rlamda(i) = lamdamax / sqrt(abs(eignvl(i)))   for i=2..iim

rlamda(1) 没有在源码中赋值，因为后续模态循环从 k=2 开始，零波数不进入滤波系数构造。

3. 初始化系数和纬带边界

coefilu/coefilv/coefilu2/coefilv2 先在全部 i=1..iim, j=1..jjm 上清零。

默认滤波纬带为：

jfiltnu = 2
jfiltsu = jjm
jfiltnv = 1
jfiltsv = jjm

随后例程在北半球和南半球分别扫描 scalar/U 纬度 rlatu 与 V/Z 纬度 rlatv。当同时满足：

cos(latitude) / colat0 < 1
rlamda(imx) * cos(latitude) < 1

对应 jfiltn* 或 jfilts* 被更新。源码明确不在 U/scalar 极点滤波；V 网格没有极点行，因此 jfiltnv 从 1 开始。

边界检查包括：

• jfiltnu <= jjm/2 + 1
• jfiltsu <= jjm + 1
• jfiltnv <= jjm/2
• jfiltsv <= jjm

任一越界都会打印错误并 STOP。

当 jfiltnv >= jjm/2 或 jfiltnu >= jjm/2 时，源码有一个 limit-case 修正：若南北滤波带在赤道附近相遇，则把 jfiltsv 或 jfiltsu 推到南侧下一行，避免同一纬带被两侧重复覆盖。源码注释本身也提示这里是边界情形处理。

4. 分配矩阵

只在 first_call_inifilr 为真时分配：

matriceun(iim,iim,jfiltnu)
matriceus(iim,iim,jjm-jfiltsu+1)
matricevn(iim,iim,jfiltnv)
matricevs(iim,iim,jjm-jfiltsv+1)
matrinvn(iim,iim,jfiltnu)
matrinvs(iim,iim,jjm-jfiltsu+1)

随后 first_call_inifilr 置为 .FALSE.。源码没有 DEALLOCATE 或尺寸变化后的重新分配逻辑；同一进程内若第二次调用时网格或滤波纬带尺寸不同，会复用第一次分配的数组形状。这在标准 GCM 启动的一次性初始化中不是问题，但在工具链或测试中重复初始化不同网格时需要特别注意。

5. 构造滤波系数

modfrstu(j) 和 modfrstv(j) 默认置为 iim，表示所有模态保留。然后对北/南 scalar/U 纬带和北/南 V/Z 纬带分别寻找第一个满足：

rlamda(k) * cos(latitude) < 1

的模态 k，将其记录为 modfrstu(j) 或 modfrstv(j)。从该模态到 modemax=iim，写入：

coefilu(k,j)  = cof - 1
coefilu2(k,j) = cof*cof - 1
coefilv(k,j)  = cof - 1
coefilv2(k,j) = cof*cof - 1

其中 cof = rlamda(k) * cos(latitude)。因此低模态不滤波，高模态从 modfrst* 起按纬度和特征值衰减。

6. 构造普通滤波矩阵

scalar/U 网格使用 eignfnv 参与矩阵构造：

eignft(i,k) = eignfnv(k,i) * coff
matriceu    = eignfnv * eignft

其中若 i < modfrstu(j)，coff 强制为 0。北半球写入 matriceun(:,:,j)，南半球写入 matriceus(:,:,j-jfiltsu+1)。

V/Z 网格使用 eignfnu：

eignft(i,k) = eignfnu(k,i) * coff
matricev    = eignfnu * eignft

北半球写入 matricevn(:,:,j)，南半球写入 matricevs(:,:,j-jfiltsv+1)。

乘法实现按编译条件选择：

• CRAY：MXM
• BLAS：SGEMM
• 否则：三重循环手写矩阵乘法

7. 构造 inverse filter 矩阵

inverse filter 只为 scalar/U 网格构造：

coff = coefilu(i,j) / (1 + coefilu(i,j))

同样使用 eignfnv * eignft 得到：

• 北半球 matrinvn(:,:,j)
• 南半球 matrinvs(:,:,j-jfiltsu+1)

源码没有为 V/Z 网格构造 matrinv 对应矩阵；串行和并行滤波例程也只在 iaire == 3 的 inverse 分支使用 matrinvn/matrinvs。

8. 可选 FFT 初始化

在 CPP_PARA 构建中，如果 use_filtre_fft 为真，inifilr 调用：

Init_filtre_fft(coefilu, modfrstu, jfiltnu, jfiltsu,
                coefilv, modfrstv, jfiltnv, jfiltsv)
Init_filtre_fft_loc(coefilu, modfrstu, jfiltnu, jfiltsu,
                    coefilv, modfrstv, jfiltnv, jfiltsv)

这说明 FFT backend 不取代本例程的系数和纬带判定，而是复用同一套 coefil*、modfrst* 和 jfiltn*/jfilts* 状态。

复现顺序

最小复现上下文必须先完成：

1. 编译维度和水平网格常量可用：dimensions.h, paramet.h。
2. 几何初始化已写入 comgeom.h 数组：rlatu, rlatv, xprimu, eignfnu, eignfnv 等。
3. logic_mod 与 serre_mod 已给出 fxyhypb, ysinus, alphax。
4. coefils.h 中的滤波状态数组已按编译维度存在。
5. 调用 inifilr。
6. 后续调用 filtreg 或 filtreg_p 时，按 iaire 和 ifiltre 选择普通矩阵、inverse 矩阵或 FFT 路径。

在 Mars 3D GCM 中，调用顺序是 conf_gcm/iniconst/inigeom 等初始化之后进入 inifilr，再进入动力主循环；Mars newstart/start2archive 工具链也在 inigeom 后调用。

风险和待确认

• first_call_inifilr 只控制分配，不控制重分配。若同一进程里对不同水平网格重复调用，数组形状可能与第二次计算出的 jfiltn*/jfilts* 不一致；源码没有保护这一情形。
• ysinus 特殊分支会把 alphax 改为 0，这是对 serre_mod 全局状态的写入副作用，不只是局部计算。
• rlamda(1) 未赋值但也不被滤波系数循环读取；审查时不应把它误判为运行期 bug。
• 待确认：Mars 标准 run.def 中 use_filtre_fft 的默认开启状态仍应由 conf_gcm 并行 主页面最终确认。

相关页面

• filtrez 总览
• filtreg-system
• filter-helpers
• fft-backends
• dyn3d/gcm
• dyn3dpar/gcm
• mars-start-restart-archive-boundary
• common-gcm-mars-run-entry