slope_mod.F90

路径

LMDZ.MARS\libf\phymars\slope_mod.F90

所属目录/模块

libf/phymars

文件定位

slope_mod.F90 提供火星地形坡面的两个核心功能：(1) 从地形势场计算坡面倾角和方位（getslopes），(2) 计算坡面上的太阳总辐照度（param_slope，Spiga & Forget 2008 方法）。这两个子程序由 physiq_mod 在 callslope=.true. 时调用，用于修正地表短波/长波辐射通量。模块还管理 theta_sl 和 psi_sl 两个 module 级数组的分配/释放。

与 comslope_mod（次网格坡面统计分布和共享状态）配合使用：slope_mod 提供几何计算，comslope_mod 提供多 bin 坡面框架。

定义的符号

依赖的模块

调用的关键例程

本模块不调用其他 LMDZ 子程序。

输入

getslopes

param_slope

输出

getslopes

param_slope

共享状态与副作用

• theta_sl 和 psi_sl：module 级 SAVE + !$OMP THREADPRIVATE 数组（行 5–8），由 ini_slope_mod 分配、end_slope_mod 释放、getslopes 写入、physiq 读取。
• getslopes 仅串行模式：并行时调用 abort_physic 终止（行 43–47）。
• param_slope 纯计算：无副作用，仅通过参数返回结果。

核心逻辑

getslopes（行 14–114）

1. 并行检查（行 43–47）：is_parallel=.true. 时终止。
2. 重排地形到 2D 网格（行 53–68）：将 geopot/g（→高度）从 1D 物理网格重排到 (nbp_lon, nbp_lat) 数组 topogrid，极点用单点填充。
3. 计算梯度（行 72–91）：
   • gradx（纬度方向）：(topogrid(i,j+1) - topogrid(i,j-1)) / (lat(j+1)-lat(j-1)) / rad
   • grady（经度方向）：类似，但经度边界用 2*pi 偏移处理环绕
   • 极点梯度设为 0
4. 计算倾角和方位（行 96–112）：
   • theta_val = atan(sqrt(gradx² + grady²))
   • psi_val：从梯度分量通过 atan2 和象限修正推导方位角，转为度并取模 360

param_slope（行 118–266）

参考：Spiga & Forget, "Fast and accurate estimation of irradiance on Martian slopes"。

1. 前置检查（行 179–182）：theta_s 超出 [0,90] 时终止。
2. 太阳低于地平线（行 185–191）：csza < 0.01 时所有分量为 0。
3. 坡-太阳方位角（行 197–202）：a = psi_s + atan2(cos(declin)*sin(rho), ...)，处理 atan2(0,0) 特例。
4. 坡-太阳相位角余弦（行 205–206）：mu_s = csza*cos(theta_s) - cos(a)*sin(theta_s)*sqrt(1-csza²)
5. 天空可视因子（行 209）：sigma_s = 0.5*(1+cos(theta_s))
6. 直接辐射（行 212）：fdir = fdir_0 * mu_s/csza
7. 反射辐射（行 215）：fref = albedo*(1-sigma_s)*ftot_0
8. 散射辐射（行 218–253）：
   • 散射分量 fscat_0 = ftot_0 - fdir_0
   • 构建散射向量 s_vector、几何向量 g_vector、耦合矩阵 mat_T = mat_M + csza*mat_N
   • 矩阵系数按 csza >= 0.5 和 < 0.5 分两组（经验拟合参数）
   • 低倾角（theta_s ≤ 5°）用线性过渡避免数值问题
   • ratio = s_vector · mat_T · g_vector
   • fscat = ratio * fscat_0
9. 总辐射（行 257）：ftot = fdir + fref + fscat

ini_slope_mod / end_slope_mod（行 270–290）

分配/释放 theta_sl 和 psi_sl，由 phys_state_var_init_mod 调用。

伪代码

! getslopes
abort if parallel
topogrid = reshape(geopot/g, [nbp_lon, nbp_lat])  ! 极点用单点填充
for interior points:
  gradx = (topo(j+1)-topo(j-1)) / (lat(j+1)-lat(j-1)) / rad
  grady = (topo(i+1)-topo(i-1)) / (lon(i+1)-lon(i-1)) / rad  ! 经度环绕修正
poles: gradx=grady=0
theta_sl = atan(sqrt(gradx²+grady²))
psi_sl = azimuth_from_gradient(gradx, grady) → [0,360)°

! param_slope
if csza < 0.01: ftot=0; return
a = psi_s + atan2(sun_geometry)                    ! 坡-太阳方位
mu_s = csza*cos(θ) - cos(a)*sin(θ)*sin(SZA)       ! 坡-太阳相位角余弦
sigma_s = 0.5*(1+cos(θ))                           ! 天空可视因子
fdir = fdir_0 * mu_s/csza                          ! 直接
fref = albedo*(1-sigma_s)*ftot_0                   ! 反射
fscat_0 = ftot_0 - fdir_0                          ! 散射（水平面）
mat_T = mat_M + csza*mat_N                         ! 耦合矩阵
ratio = s_vector · mat_T · g_vector                ! 散射比
fscat = ratio * fscat_0                            ! 散射（坡面）
ftot = fdir + fref + fscat

参与的主题流程

写法特点

• !$OMP THREADPRIVATE：theta_sl 和 psi_sl 为线程私有（行 8）。
• 仅串行：getslopes 明确要求串行模式（行 43–47），不支持并行。
• 经度环绕处理：grady 边界点用 2*pi 偏移处理 -180°/180° 环绕（行 77–80）。
• 经验拟合矩阵：param_slope 的 mat_M/mat_N 系数来自 Spiga & Forget (2008) 的拟合结果，按 csza 阈值分两组。
• 低角度线性过渡：theta_s ≤ 5°（0.0872664626 rad）时用线性插值避免数值不稳定（行 241–245）。
• id2/idm1 trick：行 41 注释说这是为了在 testphys1d 调试模式下编译通过。

复现要点

• callslope=.true.（callphys.def）时才调用 getslopes 和 param_slope。
• getslopes 仅在 3D 串行模式下工作；1D 和并行模式不调用。
• param_slope 的拟合矩阵系数是固定的（来自文献），不可调。
• 长波辐射修正使用 sky_slope（comslope_mod 中的天空可视因子），而非本模块。

待确认

• 行 199 a = deg2rad*psi_s 在 atan2(0,0) 时的处理（注释说"some compiler need specifying value"，待确认：哪些编译器有此问题）。
• param_slope 中 mu_s <= 0 时设为 0 但不终止（行 206），后续 fdir = fdir_0*mu_s/csza 会得到 0（待确认：是否为有意设计，允许背面坡面）。
• getslopes 在极点处 gradx=grady=0（行 88–91），意味着极点倾角恒为 0（待确认：是否因极点地形数据不可靠）。

相关页面

• comslope_mod — 次网格坡面参数化共享状态模块，管理多 bin 坡面框架
• vdifc-water-surface-exchange — 涉及坡面投影的地表通量
• co2condens — CO2 凝结中使用 subslope_dist 的坡面循环
• water-cycle — 水循环主题页
• co2-cycle — CO2 循环主题页
• radiation — 辐射计算主题页