calcstormfract_mod.F90

快速理解

它做什么： 计算每格点被火箭式尘暴占据的亚格点面积分数。被 physiq 在 rdstorm 时调用（physiq_mod.F:1069）。

基本过程： 取 stormdust_mass 廓线最大值 → 用参考值归一化 → 限制在 [1e-2, 0.6]。

关键结果： totstormfract(ngrid)，供 callradite/aeropacity/rocketduststorm 把格点平均换算为风暴内部浓度。

路径

LMDZ.MARS\libf\phymars\calcstormfract_mod.F90

所属目录/模块

libf\phymars

文件定位

该文件定义火箭式尘暴（rocket dust storm）亚格点面积分数计算模块。它只提供一个例程 calcstormfract：从每个水平格点的 stormdust_mass 垂直廓线中取第 2 层到顶层的最大绝对质量混合比，用参考质量混合比 5.e-4 归一化，再把结果限制在 [1.e-2, 0.6]。输出 totstormfract(ngrid) 表示每个网格中被风暴尘埃占据的面积分数。

physiq 在辐射计算前、且 rdstorm 开启时调用本例程。后续 callradite / aeropacity / rocketduststorm 使用 totstormfract 把格点平均 stormdust 转换为风暴内部浓度或风暴内部光学厚度。

定义的符号

符号	类型	行号	作用
`calcstormfract_mod`	module	1	火箭式尘暴面积分数计算模块
`calcstormfract`	subroutine	20	根据 `stormdust_mass` 垂直廓线计算 `totstormfract`

依赖的模块

use 模块	only 列表	用途	待确认
`tracer_mod`	`igcm_stormdust_mass`	定位 `pq(:,:,igcm_stormdust_mass)` 风暴尘埃质量混合比 tracer	否

调用的关键例程

被调用例程	所在模块/文件	调用位置	作用
无	-	-	本例程只使用 Fortran intrinsic：`maxval`、`abs`、`max`、`min`

输入

输入	来源	类型/维度	单位	含义
`ngrid`	调用方 `physiq`	integer scalar	-	水平格点数
`nlayer`	调用方 `physiq`	integer scalar	-	垂直层数
`nq`	调用方 `physiq`	integer scalar	-	tracer 数；本例程只用于声明 `pq` 维度
`pq(ngrid,nlayer,nq)`	调用方 `physiq`	real array	kg/kg	advected tracer 场；本例程只读取 `igcm_stormdust_mass` 分量
`igcm_stormdust_mass`	`tracer_mod`	integer scalar	-	`stormdust_mass` tracer 索引

输出

输出	去向	类型/维度	单位	含义
`totstormfract(ngrid)`	`physiq`、`callradite`、`rocketduststorm`、`aeropacity`	real array	-	网格中火箭式尘暴占据的面积分数，限制在 0.01 到 0.6

共享状态与副作用

只读 tracer_mod 中的 igcm_stormdust_mass。
不写 module 变量、common block、文件或日志。
无 SAVE、无 THREADPRIVATE、无 firstcall 初始化。
无显式错误处理；若 nlayer < 2，切片 pq(ig,2:nlayer,...) 不具备正常物理含义。复现 1 层测试时需额外确认。

核心逻辑

设置三个硬编码参数：
- mmr_ref = 5.e-4：参考质量混合比，源码注释说明对应 tau=10 的 OMEGA 本地尘暴。
- fracmax = 0.6：最大风暴面积分数。
- fracmin = 1.e-2：最小风暴面积分数。
对每个水平格点 ig：
- 读取 pq(ig,2:nlayer,igcm_stormdust_mass)，即跳过第 1 层，从第 2 层到顶层取风暴尘埃质量混合比。
- 取该垂直切片的最大值，再取绝对值：abs(maxval(...))。
- 除以 mmr_ref 得到未限制的面积分数。
- 先用 max(..., fracmin) 施加下限 0.01。
- 再用 min(..., fracmax) 施加上限 0.6。

源码行 56 的注释写有 why from 2nd layer?，说明跳过第 1 层的理由在源码内未解释。

伪代码

calcstormfract(ngrid, nlayer, nq, pq, totstormfract):
    mmr_ref = 5.e-4
    fracmin = 1.e-2
    fracmax = 0.6

    for ig in 1..ngrid:
        storm_mmr = abs(maxval(pq[ig, 2:nlayer, igcm_stormdust_mass]))
        fraction = storm_mmr / mmr_ref
        fraction = max(fraction, fracmin)
        fraction = min(fraction, fracmax)
        totstormfract[ig] = fraction

参与的主题流程

主题	参与方式
尘埃循环	为火箭式尘暴分支提供亚格点面积分数
辐射	`aeropacity` 用 `totstormfract` 把 stormdust 光学厚度换算为风暴区域内部光学厚度
火箭式尘暴输运	`rocketduststorm` 用 `totstormfract` 在格点平均浓度和风暴内部浓度之间转换

写法特点

.F90 后缀但代码排版保留固定格式风格缩进。
纯函数式辅助例程，无分配、无 I/O、无诊断。
nq 只用于 pq 形状声明，例程主体不直接读取。
totstormfract 永远至少为 0.01，即没有 stormdust 的格点也不会返回 0；这是为了后续除以 totstormfract 时避免除零。
totstormfract 最大为 0.6，对应最多 60% 网格面积被风暴尘埃占据。
第 1 层被排除在 maxval 之外，源码只留下疑问注释，未说明物理理由。

复现要点

必须使用 stormdust_mass，不是 stormdust_number 或背景 dust_mass。
垂直范围按源码为 2:nlayer；不要改成 1:nlayer，除非另有物理依据。
下限和上限顺序为先 max(fracmin) 再 min(fracmax)。
低 stormdust 或无 stormdust 时返回 0.01，会影响 aeropacity 中 aerosol(...stormdust...)/totstormfract 以及 rocketduststorm 中 stormdust/totstormfract。
abs(maxval(...)) 只对最大值取绝对值；如果数组中全是负值，则结果是最接近 0 的负值取绝对值，而不是最大绝对值。正常 tracer 应非负，但复现异常输入时应按源码字面执行。

待确认

为什么 maxval 跳过第 1 层（源码注释直接写 why from 2nd layer?）。
abs(maxval(...)) 是否有意而不是 maxval(abs(...))；正常非负 tracer 下二者等价。
1D 或测试配置若 nlayer=1 时该切片如何处理，当前源码没有保护。