filtreg operators 滤波施加与特征基构造

源码范围：

LMDZ.COMMON-6.3\LMDZ.COMMON\libf\filtrez\filtreg.F
LMDZ.COMMON-6.3\LMDZ.COMMON\libf\filtrez\inifgn.F
LMDZ.COMMON-6.3\LMDZ.COMMON\libf\filtrez\jacobi.F90
LMDZ.COMMON-6.3\LMDZ.COMMON\libf\filtrez\eigen_sort.F
LMDZ.COMMON-6.3\LMDZ.COMMON\libf\filtrez\acc.F

源码行数：filtreg.F 322 行，inifgn.F 105 行，jacobi.F90 106 行，eigen_sort.F 32 行，acc.F 22 行。

Mars 运行参与度：条件经过。这些例程服务于动力经向滤波：inifgn/jacobi/acc/eigen_sort 在 inifilr 初始化阶段间接经过；串行 filtreg 在动力 tendency、耗散、pressure/energy helper 和 tracer 后处理路径中按调用点条件经过。并行实际滤波使用 dyn3dpar/filtreg_p.F，但复用同一套初始化矩阵和 coefils.h 状态。

文件职责

文件	主要符号	职责
`filtreg.F`	`filtreg`	对已初始化的纬带执行矩阵滤波或 inverse scalar filter，写回 `champ`。
`inifgn.F`	`inifgn`	从 `xprimu/xprimv` 构造 U/V 差分矩阵，求离散 Laplacian 特征值和特征函数。
`jacobi.F90`	`JACOBI`	对称矩阵 Jacobi 特征分解，最多 50 sweep。
`acc.F`	`acc`	对每个特征向量列做 L2 归一化。
`eigen_sort.F`	`eigen_sort`	按特征值从大到小重排特征值和对应特征向量列。

数据约定

`filtreg` 输入参数

SUBROUTINE filtreg(champ, nlat, nbniv, ifiltre, iaire, griscal, iter)

参数	方向	约束	含义
`champ(iip1,nlat,nbniv)`	inout	第 `iip1` 经度为周期闭合列	待滤波场，输出原地覆盖。
`nlat`	in	scalar/U 场必须为 `jjp1`，V/Z 场必须为 `jjm`	纬向维度。
`nbniv`	in	通常为 `llm` 或 1	垂直层数或字段层数。
`ifiltre`	in	只支持 `2` 和 `-2`	`2` 为普通滤波，`-2` 为 inverse filter。
`iaire`	in	`1` intensive，`2` extensive	决定使用 `sdd` 还是 `unsdd` 做前后缩放。
`griscal`	in	`.TRUE.` scalar/U；`.FALSE.` V/Z	决定纬带、矩阵族和 `nlat` 校验。
`iter`	in	只支持 `1`	`iter==2` 直接停止并提示使用旧滤波路径。

源码注释仍保留 ifiltre=±1 的 transform 说明，但当前代码遇到 ifiltre==1 或 -1 会 STOP 'Pas de transformee simple dans cette version'。因此当前可复现路径只应使用 2/-2。

共享状态

filtreg include coefils.h，读取：

jfiltnu/jfiltsu/jfiltnv/jfiltsv：北/南半球滤波纬带边界；
sddu/sddv/unsddu/unsddv：经向尺度缩放及其倒数；
coefil*, modfrst*, eignfn*：主要由初始化与 helper 使用，filtreg 本身直接使用纬带和缩放。

filtreg 同时从 filtreg_mod 读取：

ordinary scalar/U 矩阵：matriceun, matriceus
ordinary V/Z 矩阵：matricevn, matricevs
inverse scalar/U 矩阵：matrinvn, matrinvs

filtreg 内部还有 SAVE 状态：

first
sdd12(iim,4)

第一次进入时把 sddu/sddv/unsddu/unsddv 拷入 sdd12，之后复用这份缓存。若同一进程改变网格后重复初始化，filtreg_mod 的矩阵尺寸和 filtreg 的 sdd12 缓存都没有自动刷新机制。

`filtreg` 施加流程

1. 拒绝未支持模式

例程一开始执行三个硬约束：

ifiltre==1 或 ifiltre==-1：停止，当前版本没有 direct/inverse transform。
iter==2：停止，当前版本不做二次迭代滤波。
ifiltre==-2 .AND. .NOT. griscal：停止，inverse filter 只支持 scalar/U 网格。

随后若 ifiltre 不是 2 或 -2，同样停止。

2. 选择网格、纬带和缩放

若 griscal=.TRUE.：

nlat 必须等于 jjp1；
北半球纬带为 j=2..jfiltnu；
南半球纬带为 j=jfiltsu..jjm；
iaire=1 时前缩放用 sddv，后缩放用 unsddv；
iaire!=1 时前缩放用 unsddv，后缩放用 sddv。

若 griscal=.FALSE.：

nlat 必须等于 jjm；
北半球纬带为 j=1..jfiltnv；
南半球纬带为 j=jfiltsv..jjm；
iaire=1 时前缩放用 sddu，后缩放用 unsddu；
iaire!=1 时前缩放用 unsddu，后缩放用 sddu。

这里的 U/scalar 与 V/Z 交叉使用 sddu/sddv 来自离散算子的对偶关系：inifgn 用 xprimu/xprimv 形成两套差分矩阵，inifilr 再用对应特征函数构造滤波矩阵。

3. 分半球滤波

hemisph=1 处理北半球滤波带，hemisph=2 处理南半球滤波带。每个半球先对 champ(1:iim,j,l) 乘前缩放。

矩阵乘法选择如下：

条件	北半球矩阵	南半球矩阵	写入
`ifiltre==-2`	`matrinvn(:,:,j)`	`matrinvs(:,:,j-jfiltsu+1)`	`eignq(:,j-jdfil+1,:)`
`ifiltre==2 .AND. griscal`	`matriceun(:,:,j)`	`matriceus(:,:,j-jfiltsu+1)`	`eignq(:,j-jdfil+1,:)`
`ifiltre==2 .AND. .NOT.griscal`	`matricevn(:,:,j)`	`matricevs(:,:,j-jfiltsv+1)`	`eignq(:,j-jdfil+1,:)`

编译有 BLAS 时使用 SGEMM，否则使用 Fortran matmul。

4. 写回和周期闭合

普通滤波：

champ = (champ + eignq) * post_scale

inverse filter：

champ = (champ - eignq) * post_scale

每个半球处理完后设置：

champ(iip1,j,l) = champ(1,j,l)

即只对被滤波纬带重建周期闭合列；非滤波纬带原样保留。

`inifgn` 特征基构造

inifgn(dv) 是 filtreg_mod 中 inifilr 的前置 helper。它写入 coefils.h 中的 sddu/sddv/unsddu/unsddv/eignfnu/eignfnv，并返回 dv(iim) 给 inifilr 用作 scalar/V 相关特征值。

步骤：

从 xprimu/xprimv 取 dlonu/dlonv。
计算：

sddu(i)   = sqrt(xprimu(i))
sddv(i)   = sqrt(xprimv(i))
unsddu(i) = 1 / sddu(i)
unsddv(i) = 1 / sddv(i)

初始化一阶差分形态的 eignfnv：

eignfnv(1,1)   = -1
eignfnv(iim,1) =  1
eignfnv(i+1,i+1) = -1
eignfnv(i,i+1)   =  1
eignfnv(i,j) = eignfnv(i,j) / (sddu(i) * sddv(j))

设置 eignfnu(i,j) = -eignfnv(j,i)。
形成两个对称算子：

vec  = eignfnu * eignfnv
vec1 = eignfnv * eignfnu

CRAY 分支使用 MXM，普通分支使用三重循环。

对 vec 调 jacobi，得到 dv 和 eignfnv；归一化 eignfnv；按 dv 排序。
对 vec1 调 jacobi，得到 du 和 eignfnu；归一化 eignfnu；按 du 排序。

旧 IMSL EVCSF 路径已在源码中注释掉，当前主路径是 jacobi + acc + eigen_sort。

特征求解 helper

`JACOBI`

JACOBI(A,N,NP,D,V,NROT) 对称矩阵特征分解：

初始化 V 为单位阵；
D 和 B 取 A 的对角；
每一 sweep 求上三角非对角绝对值和 SM；
SM==0 时返回；
前 3 次 sweep 用 TRESH=0.2*SM/N**2，之后阈值为 0；
对超过阈值的 A(IP,IQ) 做 Jacobi rotation，并同步更新 A、D、V；
最多 50 次 sweep，若仍不收敛则 STOP 'Jacobi: 50 iterations should never happen'。

源码注释写着 50? I suspect this should be NP，但实际实现固定 50 次；这里应按源码行为记录，不推断为 bug。

`acc`

acc(vec,d,im) 对每一列 j：

d(i) = vec(i,j)**2
sum  = sqrt(SSUM(im,d,1))
vec(i,j) = vec(i,j) / sum

它依赖外部 SSUM，作用是把特征向量列归一化，而不是累加诊断。

`eigen_sort`

eigen_sort(d,v,n,np) 对 i=1..n-1 扫描 j=i+1..n，选择最大的 d(j) 与当前位置交换，并同步交换 v(:,i) 和 v(:,k)。因此排序方向是特征值降序。

典型调用点

调用点	参数形态	含义
`dyn3d_common/bernoui.F`	`filtreg(pbern, jjp1, llm, 2, 1, .true., 1)`	对 Bernoulli 函数做 scalar/intensive 普通滤波。
`dyn3d/dteta1.F`	`filtreg(dteta, jjp1, llm, 2, 2, .true., 1)`	对位温通量 tendency 做 scalar/extensive 普通滤波。
`dyn3d/caladvtrac.F`	`filtreg(finmasse, jjp1, llm, -2, 2, .TRUE., 1)`	对质量场走 scalar inverse filter，用于水 tracer 后处理准备。
`dyn3d_common/inidissip.F90`	`zh/zu` 用 scalar，`zv` 用 V/Z	用随机场估计耗散算子的响应时调用普通滤波。
`dyn3d_common/rotatf.F`	`filtreg(rot, jjm, klevel, 2, 2, .FALSE., 1)`	对 V/Z 网格旋度结果做普通滤波。
`dyn3d_common/divergf.F`	`filtreg(div, jjp1, klevel, 2, 2, .TRUE., 1)`	对 scalar divergence 结果做普通滤波。

源码中还有一些注释掉的 CALL filtreg，例如 leapfrog.F、integrd.F、advtrac.F90 和并行对应文件内的历史路径；这些不应算作当前活动调用链。

复现要点

必须先调用 inigeom 之类的几何初始化，使 xprimu/xprimv 和 comgeom.h 中的数组有效。
必须先调用 inifilr，它内部调用 inifgn 并分配/填充矩阵。
调 filtreg 时 griscal 和 nlat 必须匹配，否则直接停止。
ifiltre=-2 只能用于 griscal=.TRUE. 的 scalar/U 网格。
并行 GCM 不调用串行 filtreg.F 施加滤波，而是走 dyn3dpar/filtreg_p.F；但特征基和矩阵初始化仍来自同一套 inifgn/inifilr 逻辑。

风险和待确认

filtreg 缓存 sdd12 只在第一次调用时填充；同一进程内若更换水平网格并重复初始化，缩放缓存不会自动更新。
jacobi 50 sweep 不收敛会直接 STOP。没有错误码向上层传播。
待确认：eigen.F 旧 helper 仍在 filter-helpers.md 中列出，但本次源码确认的主路径未发现活动调用点。