slope_mod.F90
快速理解
它做什么： 提供坡面倾角/方位计算（getslopes）和坡面太阳辐照度（param_slope，Spiga & Forget 2008）。被 physiq 调用（L787/L1229/L1270）。
基本过程： 从地形势场计算坡面几何 → 按坡面几何修正太阳直接/散射通量。
关键结果： 坡面修正后的短波/长波辐射通量；与 comslope_mod（次网格坡面统计）配合。
路径
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 坡面框架。
定义的符号
| 符号 | 类型 | 行号 | 作用 |
|---|---|---|---|
slope_mod |
module | 1 | 坡面计算模块 |
theta_sl |
real, save, allocatable | 5 | 坡面倾角数组（度），0=水平，90=垂直 |
psi_sl |
real, save, allocatable | 6 | 坡面方位数组（度），0=北，90=东，180=南，270=西 |
getslopes |
subroutine | 14 | 从地形势场计算坡面倾角和方位 |
param_slope |
subroutine | 118 | 计算坡面上的太阳总辐照度 |
ini_slope_mod |
subroutine | 270 | 分配 theta_sl 和 psi_sl 数组 |
end_slope_mod |
subroutine | 283 | 释放 theta_sl 和 psi_sl 数组 |
依赖的模块
| use 模块 | only 列表 | 用途 | 待确认 |
|---|---|---|---|
geometry_mod |
longitude, latitude |
网格点经纬度（弧度），getslopes 使用 |
|
comcstfi_h |
g, rad, pi |
重力加速度、行星半径、圆周率 | |
mod_phys_lmdz_para |
is_parallel |
判断是否并行模式（getslopes 仅支持串行） |
|
mod_grid_phy_lmdz |
nbp_lon, nbp_lat |
经纬度网格点数 | |
comcstfi_h |
pi |
param_slope 使用 |
调用的关键例程
| 被调用例程 | 所在模块/文件 | 调用位置 | 作用 |
|---|---|---|---|
abort_physic |
物理基础设施 | 行 46, 181 | 并行模式或非法倾角时终止执行 |
本模块不调用其他 LMDZ 子程序。
输入
getslopes
| 输入 | 来源 | 类型/维度 | 单位 | 含义 |
|---|---|---|---|---|
ngrid |
physiq |
integer | — | 水平网格点数 |
geopot(ngrid) |
physiq（phisfi） |
real | m²/s² | 地表势场 |
param_slope
| 输入 | 来源 | 类型/维度 | 单位 | 含义 |
|---|---|---|---|---|
csza |
physiq（mu0(ig)） |
real | — | 太阳天顶角余弦 |
declin |
physiq |
real | rad | 太阳赤纬 |
rho |
physiq（sl_ra） |
rad | rad | 太阳赤经 |
latitude |
physiq（sl_lat） |
real | deg | 纬度 |
taudust |
physiq（sl_tau） |
real | — | 0.67 μm 参考波长 dust 光学厚度 |
albedo |
physiq（sl_alb） |
real | — | 地表反照率 |
theta_s |
theta_sl(ig) |
real | deg | 坡面倾角（0–90） |
psi_s |
psi_sl(ig) |
real | deg | 坡面方位（0=北, 90=东, 180=南, 270=西） |
fdir_0 |
physiq（sl_di0） |
real | W/m² | 水平面上直接辐照度 |
ftot_0 |
physiq（sl_fl0） |
real | W/m² | 水平面上总辐照度 |
输出
getslopes
| 输出 | 去向 | 类型/维度 | 单位 | 含义 |
|---|---|---|---|---|
theta_sl(ngrid) |
module 变量 | real | deg | 坡面倾角 |
psi_sl(ngrid) |
module 变量 | real | deg | 坡面方位 |
param_slope
| 输出 | 去向 | 类型/维度 | 单位 | 含义 |
|---|---|---|---|---|
ftot |
physiq（sl_flu） |
real | W/m² | 坡面上的总辐照度 |
共享状态与副作用
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）
- 并行检查（行 43–47）：
is_parallel=.true.时终止。 - 重排地形到 2D 网格（行 53–68）：将
geopot/g（→高度）从 1D 物理网格重排到(nbp_lon, nbp_lat)数组topogrid，极点用单点填充。 - 计算梯度（行 72–91）：
gradx（纬度方向）：(topogrid(i,j+1) - topogrid(i,j-1)) / (lat(j+1)-lat(j-1)) / radgrady（经度方向）：类似，但经度边界用2*pi偏移处理环绕- 极点梯度设为 0
- 计算倾角和方位（行 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"。
- 前置检查（行 179–182）：
theta_s超出 [0,90] 时终止。 - 太阳低于地平线（行 185–191）：
csza < 0.01时所有分量为 0。 - 坡-太阳方位角（行 197–202）：
a = psi_s + atan2(cos(declin)*sin(rho), ...)，处理atan2(0,0)特例。 - 坡-太阳相位角余弦（行 205–206）：
mu_s = csza*cos(theta_s) - cos(a)*sin(theta_s)*sqrt(1-csza²) - 天空可视因子（行 209）：
sigma_s = 0.5*(1+cos(theta_s)) - 直接辐射（行 212）：
fdir = fdir_0 * mu_s/csza - 反射辐射（行 215）：
fref = albedo*(1-sigma_s)*ftot_0 - 散射辐射（行 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_vectorfscat = ratio * fscat_0
- 散射分量
- 总辐射（行 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参与的主题流程
| 主题 | 参与方式 |
|---|---|
| 辐射计算 | param_slope 修正坡面上的短波辐射通量（直接+散射+反射） |
| 地表能量收支 | getslopes 提供坡面几何，physiq 用 theta_sl 修正长波辐射（天空可视因子） |
写法特点
!$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/idm1trick：行 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 — 辐射计算主题页