flusv.F
快速理解
它做什么: 红外散射通量求解器,求解 n 层间的向上/向下红外通量。被 lwdiff 调用。
基本过程: 假设层内黑体亮度随光学厚度线性变化 → 两半球通量方法(Toon et al. 1988)。
关键结果: 向上/向下红外通量,累加回 PFLUC。
路径
LMDZ.MARS\libf\phymars\flusv.F
所属目录/模块
libf\phymars
文件定位
flusv.F 定义 flusv_mod,提供长波红外散射通量求解器。主例程 flusv 在给定每层单次散射反照率、不对称因子、光学厚度、层界黑体亮度和地表发射率后,求解 n 层之间的向上/向下红外通量。源码注释说明它假设每层内黑体亮度 B 随光学厚度 tau 线性变化,顶层向下通量为零,层号从大气顶到地面。
当前只读源码搜索显示外部调用点在 lwdiff.F:lwdiff 对 CO2 15 微米带外的红外谱带 iir=3,nir 组装 NDD=nlaylte*2 个扩散子层,然后以 nsf=0 调用 flusv 计算 hemispheric constant 通量,再把结果累加回 PFLUC。
定义的符号
| 符号 | 类型 | 行号 | 作用 |
|---|---|---|---|
flusv_mod |
module | 1 | 包含红外散射通量求解例程和内部三对角求解器。 |
flusv |
subroutine | 7 | 组装双流红外散射边值问题,求解每个层界的向上/向下通量。 |
sys3v |
subroutine | 274 | 向量化求解三对角线性系统,供 flusv 的 2*n 个积分常数求解使用。 |
依赖的模块
| use 模块 | only 列表 | 用途 | 待确认 |
|---|---|---|---|
dimradmars_mod |
ndlo2, ndlon, nflev |
定义辐射数组列数和工作数组层数。flusv 与 sys3v 都直接用这些维度声明实参数组和临时数组。 |
否 |
调用的关键例程
| 被调用例程 | 所在模块/文件 | 调用位置 | 作用 |
|---|---|---|---|
sys3v |
flusv.F |
flusv 第 162 行 |
求解 flusv 组装出的 2*n 维三对角系统,返回每层两个积分常数 y(:,2*i-1:2*i)。 |
输入
| 输入 | 来源 | 类型/维度 | 单位 | 含义 |
|---|---|---|---|---|
KDLON |
调用方 | INTEGER |
- | 本次向量化计算的水平列数;循环范围为 1:KDLON。 |
nsf |
调用方 | INTEGER |
- | 通量算法选择。0 只用 hemispheric constant 直接公式;>0 会进入 source function 角向积分分支。 |
n |
调用方 | INTEGER |
layer | 扩散层数;flusv 组装大小为 2*n 的线性系统。 |
omega |
调用方 | REAL(NDLO2,n) |
- | 每层单次散射反照率。 |
g |
调用方 | REAL(NDLO2,n) |
- | 每层散射不对称因子。 |
tau |
调用方 | REAL(NDLO2,n) |
- | 每层光学厚度。 |
emis |
调用方 | REAL(NDLO2) |
- | 地表发射率;底边界同时使用 emis*pi*bsol 和 (1-emis) 反射项。 |
bh |
调用方 | REAL(NDLO2,n+1) |
blackbody luminance | 层界黑体亮度,bh(:,n+1) 是剖面对应的地面值。 |
bsol |
调用方 | REAL(NDLO2) |
blackbody luminance | 实际地表黑体亮度,可不同于 bh(:,n+1)。 |
输出
| 输出 | 去向 | 类型/维度 | 单位 | 含义 |
|---|---|---|---|---|
fah |
调用方 | REAL(NDLO2,n+1) |
flux | 每个层顶的向上红外通量,fah(:,n+1) 为地面向上通量。 |
fdh |
调用方 | REAL(NDLO2,n+1) |
flux | 每个层顶的向下红外通量,fdh(:,1) 在顶边界为零,fdh(:,n+1) 为地面向下通量。 |
y |
sys3v 调用方 |
REAL(NDLO2,n) |
- | 三对角系统解;在 flusv 中对应每层两个积分常数。 |
共享状态与副作用
flusv.F 不定义保存型 module 变量,不读写文件,不打印日志,也不读取配置开关。副作用限于写入输出数组 fah/fdh,以及 sys3v 写入输出数组 y。所有临时状态都在例程局部数组中。
核心逻辑
- 对每层、每个向量列计算双流散射系数:
beta=(1-g)/2,gama1=2*(1-omega*(1-beta)),gama2=2*omega*beta,alambda=sqrt(gama1**2-gama2**2),并由grgama=(gama1-alambda)/gama2得到指数解的耦合系数。 - 把每层黑体亮度表示成
b0+b1*tau。当tau>1.E-3时用层顶和层底黑体亮度差分给出b1;当光学厚度太小时强制等温层,令b0=(bh(i)+bh(i+1))/2、b1=0,避免dB/dtau过大。 - 计算每层边界组合系数
e1/e2/e3/e4、非齐次项cah/cab/cdh/cdb,以及 source-function 分支需要的grg/grh/grj/grk、alpha1/alpha2、sigma1/sigma2。 - 组装
2*n维三对角系统:第 1 行施加大气顶向下通量为零;中间行连接相邻层界;第2*n行施加地表发射和反射边界条件。 - 调用
sys3v(KDLON,2*n,a,b,d,e,y)求出每层两个积分常数。 - 用
y(:,2*i-1)和y(:,2*i)更新grg/grh/grj/grk。 - 若
nsf==0,直接用 hemispheric constant 公式填充所有层界的fah/fdh。 - 若
nsf>0,先从大气顶向下沿 8 个固定高斯角度积分 source function 得到fdh,再用地表反射条件得到fah(:,n+1),最后从地面向上积分得到各层fah。 sys3v从系统底部向上做消元,保存as/ds,再从顶部向下回代生成y。
伪代码
for each layer i and column iv:
derive two-stream scattering coefficients from omega, g, tau
if tau is large enough:
linearly interpolate blackbody luminance across optical depth
else:
use isothermal layer approximation
build layer transfer/source coefficients
assemble top boundary, internal interface rows, and ground boundary
solve the 2*n tridiagonal system with sys3v
if nsf == 0:
compute upward/downward interface fluxes directly from solved constants
else:
integrate downward intensities over 8 quadrature angles
apply ground emission plus reflected downward flux
integrate upward intensities over 8 quadrature angles
参与的主题流程
| 主题 | 参与方式 |
|---|---|
| 辐射计算 | 长波辐射辅助求解器;lwdiff.F 在 CO2 15 微米带外红外谱带中调用它计算散射通量。 |
| phymars 核心物理模块目录 | 属于 libf\phymars 的辐射辅助文件。 |
写法特点
- 文件为固定格式 Fortran,模块内含两个子例程。
nq=8、x(8)、w(8)是 source-function 分支的固定高斯积分节点和权重。- 工作数组按
NDLON、NDLO2、nflev静态声明;线性系统数组维度为4*nflev,层系数数组维度为2*nflev。 - 小光学厚度保护是核心数值稳定写法:
tau<=1.E-3时不再用差分b1=(bh(i+1)-bh(i))/tau。 sys3v是专用三对角求解器,没有显式 pivoting 或零分母保护。
复现要点
KDLON必须不超过NDLON/NDLO2所能容纳的列数。n必须与nflev静态工作数组匹配;flusv会访问2*n的系统行和n+1个层界。nsf会改变算法路径:当前lwdiff.F调用传入0,因此运行的是直接 hemispheric constant 通量公式,不进入 8 点 source-function 积分。tau<=1.E-3的等温层近似是防止dB/dtau爆大的源码保护,复现实验中不要无意移除。- 复现风险:
grgama=(gama1-alambda)/gama2没有保护gama2=0的情况;当omega*beta为零时存在除零风险。 - 复现风险:source-function 分支含
alambda*x(j)-1分母,源码没有显式奇异保护。 - 复现风险:
sys3v的b(iv,n)与b(iv,i)-d(iv,i)*as(iv,i+1)分母没有 pivot 或零分母检查。
待确认
bh/bsol/fah/fdh的物理单位由调用链传入,源码注释只称black body luminance和 flux;本页不额外推断具体单位。nsf>0的 source-function 分支在当前只读搜索中未见外部调用;是否仍用于某些配置或历史路径需要更大范围运行配置核验。