LMDZ.MARS Wiki

vlz_fi.F

路径

LMDZ.MARS\libf\phymars\vlz_fi.F

所属目录/模块

libf/phymars

文件定位

源码依据（第 12-14 行注释）：vlz_fi 是 LMDZ.MARS 中用于物理过程（沉降）的垂直方向 "pseudo upstream" 平流格式。

它实现 Van Leer 坡度限制的半拉格朗日输运：给定一个示踪物混合比场 q 和每个层间界面在本时间步内被穿越的大气质量 w，本例程计算各界面的示踪物通量 wq 并原地更新 q。

在系统中，下落速度公式由调用者（newsedim_mod.F）提供——本例程只负责"已知 w 后如何保守地输运 q"。它是尘埃沉降、水冰沉降和 HDO 同位素沉降共享的底层数值内核。

定义的符号

符号 类型 行号 作用
vlz_fi_mod module 1 容器模块，无模块级变量
vlz_fi subroutine 7 Van Leer 坡度限制垂直平流：已知 w，计算通量 wq 并更新 q

依赖的模块

use 模块 only 列表 用途 待确认
（无显式 USE 语句） — 本文件不 USE 任何模块 —

外部调用：abort_physic（第 162 行），来自 LMDZ.COMMON/libf/dyn3d/abort_gcm.F（通过 libf/phymars/dyn1d/abort_gcm.F 的符号链接引入）。推断：编译时由链接器解析，属于无显式模块接口的外部过程调用。

调用的关键例程

被调用例程 所在模块/文件 调用位置 作用
abort_physic abort_gcm.F（LMDZ.COMMON） 第 162 行 向上平流分支中遇到不可能的气柱质量配置时中止

输入

输入 来源 类型/维度 单位 含义
ngrid 调用者 integer（in） — 大气柱数（水平格点数）
nlay 调用者 integer（in） — 大气层数
masse 调用者 real(ngrid,nlay)（in） kg·m⁻²（推断） 每层大气质量 δ(P)/g；注释写 "mass of atmospheric layer delta(P)/g"（第 25 行）
pente_max 调用者 real（in） — 坡度限制最大斜率，推荐值 2（第 26 行注释）；两个调用者均传 2.

输出

输出 去向 类型/维度 单位 含义
q inout，调用者回收 real(ngrid,nlay)（inout） kg/kg 示踪物混合比，被原地更新为平流后的值
w inout，可能被裁剪 real(ngrid,nlay)（inout） kg·m⁻² 本时间步穿过各界面的大气质量；扩展模式中若穿越质量超过可用气柱，w 会被截断为 Mtot（第 113 行）
wq out，调用者使用 real(ngrid,nlay+1)（out） kg·m⁻²（推断；源码注释写 kg，调用者注释写 kg.m-2） 各界面示踪物通量增量；wq(ij,nlay+1)=0（模式顶无通量），wq(ij,1) 为落到地表的量

共享状态与副作用

核心逻辑

按执行顺序（行号为源码行号）：

阶段 1：梯度计算与坡度限制（第 45-68 行）

  1. 原始梯度（第 45-50 行）：对 l=2..nlay 逐层逐列计算 dzqw(ij,l) = q(ij,l-1) - q(ij,l) 方向约定：沿气压增大方向（向下），与 w 正方向一致。

  2. Van Leer 坡度限制（第 52-63 行）：对内部层 l=2..nlay-1：

    • 若 dzqw(ij,l) 与 dzqw(ij,l+1) 同号（单调区域）： dzq = 0.5 * (dzqw(ij,l) + dzqw(ij,l+1))（中心差分）
    • 否则（极值点附近）：dzq = 0（一阶迎风，消除振荡）
    • 限幅：|dzq| ≤ pente_max * min(|dzqw(ij,l)|, |dzqw(ij,l+1)|) 使用 sign(min(abs(dzq), dzqmax), dzq) 实现（第 61 行）。

  3. 边界条件（第 65-68 行）： dzq(ij,1) = 0 且 dzq(ij,nlay) = 0——模式顶和底界面不施加二阶修正。

阶段 2：向下通量计算（w > 0，沉降方向）（第 76-120 行）

顶部边界：wq(ij,nlay+1) = 0（模式顶无通量，第 77 行）。

对每层 l=1..nlay、每列 ij=1..ngrid，当 w(ij,l) > 0 时：

常规格式（w ≤ masse，即穿越质量不超过本层质量，第 89-92 行）：

sigw = w(ij,l) / masse(ij,l)
wq(ij,l) = w(ij,l) * (q(ij,l) + 0.5*(1-sigw)*dzq(ij,l))

这是标准的 Van Leer 格式：迎风值加二阶修正项，修正幅度随 sigw 增大而减小（穿越越多，越接近纯迎风）。

扩展格式（w > masse，穿越质量超过一层，第 96-116 行）：

阶段 3：向上通量计算（w < 0，非沉降方向）（第 122-169 行）

对沉降场景被跳过：第 124 行 goto 99 直接跳到阶段 4。

向上分支的数值格式是向下分支的镜像，但循环方向相反（l=nlay-1..1），且边界条件为：wq(ij,1) = 0（地表无向上通量，第 128 行）。

异常处理（第 159-163 行）：若扩展循环到达底层仍无法覆盖穿越质量，打印错误信息并调用 abort_physic 中止。

阶段 4：保守更新 q（第 171-187 行）

对每层 l=1..nlay、每列 ij=1..ngrid：

  1. 非负保护（第 181-183 行）：若净通量会导致 q 变负，截断 if (wq(l+1) - wq(l)) < -(masse(l)*q(l)) 则将 wq(l+1) = wq(l) - masse(l)*q(l)。

  2. 更新混合比（第 185 行）： q(ij,l) = q(ij,l) + (wq(ij,l+1) - wq(ij,l)) / masse(ij,l)

注释掉的行（第 175-178 行）展示了同时更新 masse 的做法（用于非沉降的完整输运场景），沉降模式下不使用——大气质量不因沉降而改变。

伪代码

subroutine vlz_fi(ngrid, nlay, q, pente_max, masse, w, wq):

  #--- 阶段 1: 坡度限制梯度 ---
  for l=2..nlay, ij=1..ngrid:
    dzqw(ij,l) = q(ij,l-1) - q(ij,l)
    adzqw(ij,l) = |dzqw(ij,l)|

  for l=2..nlay-1, ij=1..ngrid:
    if sign(dzqw(l)) == sign(dzqw(l+1)):    # 单调区域
      dzq = 0.5*(dzqw(l) + dzqw(l+1))       # 中心差分
    else:
      dzq = 0                                 # 极值点：一阶迎风
    dzq = sign(min(|dzq|, pente_max*min(|dzqw(l)|,|dzqw(l+1)|)), dzq)  # 限幅

  dzq(:,1) = 0; dzq(:,nlay) = 0               # 边界条件

  #--- 阶段 2: 向下通量 (w > 0) ---
  wq(:,nlay+1) = 0                            # 顶无通量

  for l=1..nlay, ij=1..ngrid:
    if w(ij,l) > 0:
      if w(ij,l) <= masse(ij,l):              # 常规：穿越 ≤ 1 层
        sigw = w / masse
        wq(ij,l) = w * (q(ij,l) + 0.5*(1-sigw)*dzq(ij,l))
      else:                                    # 扩展：穿越 > 1 层
        m = l; Mtot = masse(l); MQtot = masse(l)*q(l)
        while w > Mtot + masse(m+1):
          m += 1; Mtot += masse(m); MQtot += masse(m)*q(m)
        if m < nlay:
          sigw = (w - Mtot) / masse(m+1)
          wq = MQtot + (w-Mtot)*(q(m+1) + 0.5*(1-sigw)*dzq(m+1))
        else:
          w(ij,l) = Mtot; wq(ij,l) = MQtot    # 截断：气柱不够

  #--- 阶段 3: 向上通量 (w < 0) --- 沉降时跳过 (goto 99) ---

  #--- 阶段 4: 保守更新 ---
  for l=1..nlay, ij=1..ngrid:
    if wq(l+1) - wq(l) < -(masse(l)*q(l)):    # 非负保护
      wq(l+1) = wq(l) - masse(l)*q(l)
    q(ij,l) = q(ij,l) + (wq(l+1) - wq(l)) / masse(l)

参与的主题流程

主题 参与方式
尘埃循环 newsedim_mod.F:216 — call vlz_fi(ngrid,nlay,pqi,2.,masse,w,wq)，每个尘埃示踪物的重力沉降通量计算与混合比更新
水循环 同上，callsedim_mod.F 第 593/599 行通过 newsedim 间接调用
HDO 同位素沉降 callsedim_mod.F:631 — call vlz_fi(ngrid,nlay,Ratio,2.,masseq,w,wq)，直接调用，输运的是 HDO/H2O 比值 Ratio 而非混合比；w 来自母体 H2O 冰的沉降通量（第 627 行注释："hdo is transported by h2o"）
CO2 云 co2cloud_mod.F90 通过 newsedim 间接调用（CO2 冰与 CCN 的微时间步沉降）

写法特点

复现要点

待确认

相关页面