grid_noro1.F

路径

LMDZ.MARS\libf\dynphy_lonlat\phymars\grid_noro1.F

所属目录/模块

libf/dynphy_lonlat/phymars

文件定位

亚网格地形参数计算例程（无旋转网格）。grid_noro1 从高分辨率 US Navy（USN）高程数据计算 GCM 每个网格单元的 5 个亚网格地形统计参数：平均高程（zmea）、标准差（zstd）、坡度（zsig）、各向异性（zgam）和主方向角（zthe），用于中尺度山脉参数化方案。由 datareadnc 在处理 k=4（亚网格地形统计）时调用。

定义的符号

依赖的模块

INCLUDE 头文件：

调用的关键例程

输入

输出

共享状态与副作用

• 读取 comconst_mod::rad（只读）。
• 大量 print * 调试输出（维度校验、进度标记 OKM1/OK0/OK1、最终统计值）。
• 无文件 I/O、无模块变量写入。

核心逻辑

1. 维度校验：验证 imar <= 2200、jmar <= 1100、imdep == 360、jmdep == 180、imar+1 <= iim+1、jmar <= jjm+1，否则 ABORT。
2. 数据扩展：将 USN 高程数据从 (360,180) 扩展到 (360+2×40, 180+2) 的带边界数组 zusn，经度方向周期性延拓 iext=40 列，纬度方向外推两极。
3. 高纬滤波：按 ideltax = 1/cos(lat) 计算经向滤波窗口宽度（奇数），对高程做滑动平均得到 zusnfi，使高纬度梯度计算更各向同性。
4. 梯度相关计算：在 USN 网格上逐点计算经向梯度平方 zxtzxusn、纬向梯度平方 zytzyusn 和交叉梯度 zxtzyusn。
5. 网格聚合：对每个 GCM 目标网格单元，按面积权重（weighx × weighy）累加 USN 格点的高程、高程平方、梯度相关量。
6. 统计参数推导：
   • zmea = Σ(z × w) / Σw（加权平均高程）
   • zstd = sqrt(z²_mean - zmea²)（标准差）
   • xk, xl, xm = 梯度张量分量 → xp, xq = 特征值
   • zsig = sqrt(xq)（主方向坡度）
   • zgam = xp/xq（各向异性比）
   • zthe = atan2(xm, xl) / 2 × 57.3°（主方向角）
7. 极点处理：经度周期性补齐 imar+1 列；极点（j=1, jmar）做纬圈加权平均，zgam=1、zthe=0。

伪代码

SUBROUTINE grid_noro1(imdep, jmdep, xdata, ydata, entree,
                      imar, jmar, x, y, zmea, zstd, zsig, zgam, zthe)

  校验维度 → ABORT if 超限

  -- 数据扩展 --
  zusn(i+iext, j+1) ← entree(i,j) + 周期性延拓 + 极区填充

  -- 高纬滤波 --
  FOR j = 1..jusn:
    ideltax = 1/cos(lat)  (奇数化)
    zusnfi ← 滑动平均(entree, ideltax)

  -- 梯度相关 --
  FOR each USN point:
    zxtzxusn ← (dz/dx)²
    zytzyusn ← (dz/dy)²
    zxtzyusn ← dz/dx × dz/dy

  -- 面积加权聚合到 GCM 网格 --
  FOR each GCM cell (ii, jj):
    FOR each USN point in cell:
      面积权重 = overlap_x × overlap_y
      累加 z, z², 梯度量

  -- 推导统计参数 --
  FOR each GCM cell:
    zmea = weighted_mean(z)
    zstd = sqrt(mean(z²) - zmea²)
    xk, xl, xm = gradient tensor components
    xp, xq = eigenvalues
    zsig = sqrt(xq)
    zgam = xp/xq
    zthe = atan2(xm, xl)/2 in degrees

  -- 极点特殊处理 --
  极点做纬圈平均, zgam=1, zthe=0
  经度周期性补齐 imar+1
END SUBROUTINE

参与的主题流程

写法特点

• Fixed-form Fortran（72 列限制、c 注释）。
• 硬编码 USN 网格参数：iusn=360, jusn=180, iext=40（扩展边界）。
• 固定数组上限：a(2200), b(2200), c(1100), d(1100) 对应最大目标网格维度。
• 法语注释：源码全法语注释（F. Lott）。
• 梯度张量特征值：用 xk ± sqrt(xl² + xm²) 直接计算主特征值，避免完整矩阵对角化。
• 数值保护：xp <= 1e-8 时置 0、xq <= 1e-8 时置 1e-8、|xm| <= 1e-8 时保留符号的最小值。

复现要点

• 输入高程须为 360×180 全球等经纬度网格。
• datareadnc 调用时传入 imdep=360, jmdep=180，坐标已转弧度，高程已 ×1000 转 m。
• 目标网格维度 imar=iim, jmar=jjp1。
• comconst_mod::rad 为角度→弧度转换因子（π/180）。

待确认

• rad 的具体值（推测为 π/180）和定义位置。
• 极点处理中 zmeanor/zweinor 分母为零的可能性（极地无数据时）。

相关页面

• datareadnc — 唯一调用方，在 k=4 时调用
• avg_horiz_mod — 后续由 datareadnc 调用处理平均高程
• dynphy_lonlat/index — 所属目录概览