molvis.F

快速理解

它做什么： 热层风速的垂直分子黏性扩散求解器。被 thermosphere_mod 在 callmolvis 时调用。

基本过程： 复用 conduction 同源三对角算法 → 热导率换算为分子黏性系数 → 隐式求解。

关键结果： 纬向/经向速度倾向 zdumolvis/zdvmolvis，累加回 pdu/pdv。

路径

LMDZ.MARS\libf\aeronomars\molvis.F

所属目录 / 模块

libf\aeronomars

文件定位

molvis.F 定义 molvis_mod，提供热层风速的垂直分子黏性扩散求解器。子程序 molvis 对一个风速分量 pvel 做一维垂直隐式三对角求解，输出同维度的速度 tendency zdvelmolvis（m/s/s）。thermosphere_mod.F 在 callmolvis 开关为真时分别把纬向风 pu 和经向风 pv 传入本例程，得到 zdumolvis 与 zdvmolvis，再累加回 pdu 和 pdv。

源码注释说明本文件基于 conduction.F，但这里求解的是速度增量而不是温度增量：热导率系数 Akknew 会先通过 fac = 0.25*(9*Cp - 5*Cv) 换算为分子黏性相关系数，其中 Cv = cpnew-rnew，Cp 和 R 来自 conc_mod::cpnew/rnew。下边界固定为地表速度 velsurf=0.0，上边界通量 phitop=0.0。

定义的符号

符号	类型	行号	作用
`molvis_mod`	module	1（`END MODULE` 184）	包装分子黏性风速扩散例程
`molvis`	subroutine	7（`END SUBROUTINE` 182）	对一个风速分量求解垂直分子黏性扩散 tendency
`firstcall`	saved local logical	75	每个 OpenMP 线程首次调用时打印 `skk`，随后置为 `.false.`
`skk`	local parameter	71	温度指数，固定为 `0.69`
`velsurf`	local parameter	73	下边界速度，固定为 `0.0`

依赖的模块

use 模块	only 列表	用途	待确认
`conc_mod`	`cpnew`, `Akknew`, `rnew`	读取当前成分相关的定压比热、热导率系数和气体常数，用于把热导率系数换算为分子黏性系数并计算密度	-

调用的关键例程

被调用例程	所在模块 / 文件	调用位置	作用
无外部例程	-	-	本例程只做局部数组运算和 Thomas 前推/回代；除 Fortran intrinsic 外不调用外部子程序

输入

输入	来源	类型/维度	单位	含义
`ngrid`	`thermosphere_mod`	integer	-	大气列数
`nlayer`	`thermosphere_mod`	integer	-	垂直层数
`ptimestep`	`thermosphere_mod`	real	s	物理时间步
`pplay`	`thermosphere_mod`	real `(ngrid,nlayer)`	Pa	层中压力，用于 `rho=p/(R*T)`
`pplev`	`thermosphere_mod`	real `(ngrid,nlayer+1)`	Pa	层界面压力；声明后未在例程体使用
`pt`	`thermosphere_mod`	real `(ngrid,nlayer)`	K	当前层中温度
`pdt`	`thermosphere_mod`	real `(ngrid,nlayer)`	K/s	已累计温度 tendency，用于构造预测温度 `zt=pt+pdt*ptimestep`
`pvel`	`thermosphere_mod`	real `(ngrid,nlayer)`	m/s	单个风速分量；调用方分别传入 `pu` 和 `pv`
`tsurf`	`thermosphere_mod`	real `(ngrid)`	K	地表温度，参与底层 `lambda(1)` 的温度依赖系数
`zzlev`	`thermosphere_mod`	real `(ngrid,nlayer+1)`	m	层界面高度；本例程只复制 1..`nlayer`，然后把顶界面设为 `zlev(nlayer)+10000.`
`zzlay`	`thermosphere_mod`	real `(ngrid,nlayer)`	m	层中高度，用于计算层间距离
`Akknew`, `cpnew`, `rnew`	`conc_mod`	real `(ngrid,nlayer)`	mixed	成分相关热导率系数、定压比热和比气体常数

输出

输出	去向	类型/维度	单位	含义
`zdvelmolvis`	`thermosphere_mod`	real `(ngrid,nlayer)`	m/s/s	分子黏性造成的速度 tendency；调用方把它分别接到 `zdumolvis` 或 `zdvmolvis`

共享状态与副作用

firstcall 是 SAVE + THREADPRIVATE 的局部状态；每个 OpenMP 线程首次调用会向标准输出打印 molvis: coeff of molecular viscosity skk 和 skk=0.69。
phitop 是普通局部变量，每次调用在第 93 行设为 0.0，作为顶边界零动量通量。
velsurf=0.0 是硬编码底边界速度，底层方程使用 lambda(1)*(velsurf-zvel(1))，等价于无滑移地表速度条件。
本例程不写文件、不修改 conc_mod 数组，也不修改输入的 pdt/pvel；真正的风场 tendency 累加发生在调用方 thermosphere_mod.F 第 102-103 行。

核心逻辑

逐列 ig=1,ngrid 独立求解：

构造预测温度与当前风速：zt(l)=pt(ig,l)+pdt(ig,l)*ptimestep，zvel(l)=pvel(ig,l)，并复制 zzlay/zzlev 的前 nlayer 个高度。
把顶层界面重设为 zlev(nlayer+1)=zlev(nlayer)+10000.，即忽略传入的 zzlev(ig,nlayer+1)，用硬编码 10 km 顶层厚度。
计算黏性扩散系数 lambda。底层用 tsurf(ig)**skk/zlay(1)；上方各层用 zt(l)**skk/(zlay(l)-zlay(l-1))。每层都先用 fac=0.25*(9*cpnew-5*(cpnew-rnew)) 把 Akknew 换算为黏性相关系数。
计算惯性项 alpha(l)=rho(l)/ptimestep * dz(l)，其中 rho=pplay/(rnew*zt)；与 conduction.F 不同，alpha 不再乘 cpnew，因为这里扩散的是速度而不是温度。
做 Thomas 前推。底层包含地表速度边界 velsurf=0，中间层按相邻层风速差和上一层 C/D 递推，顶层用 phitop=0 的零通量边界。
自顶向下回代得到 pdvelm(l)，最后输出 zdvelmolvis(ig,l)=pdvelm(l)/ptimestep。

伪代码

molvis(..., pvel, ..., zdvelmolvis):
  if firstcall:
    print skk
    firstcall = false

  phitop = 0.0
  for each column ig:
    for l = 1..nlayer:
      zt(l) = pt(ig,l) + pdt(ig,l) * ptimestep
      zvel(l) = pvel(ig,l)
      zlay(l) = zzlay(ig,l)
      zlev(l) = zzlev(ig,l)
    zlev(nlayer+1) = zlev(nlayer) + 10000.

    fac = 0.25 * (9*cpnew(ig,1) - 5*(cpnew(ig,1)-rnew(ig,1)))
    lambda(1) = Akknew(ig,1) * tsurf(ig)**skk / zlay(1) / fac
    for l = 2..nlayer:
      fac = 0.25 * (9*cpnew(ig,l) - 5*(cpnew(ig,l)-rnew(ig,l)))
      lambda(l) = Akknew(ig,l) / fac * zt(l)**skk / (zlay(l)-zlay(l-1))

    for l = 1..nlayer:
      rho = pplay(ig,l) / (rnew(ig,l) * zt(l))
      alpha(l) = rho / ptimestep * (zlev(l+1)-zlev(l))

    den(1) = alpha(1) + lambda(2) + lambda(1)
    C(1) = (lambda(1)*(0-zvel(1)) + lambda(2)*(zvel(2)-zvel(1))) / den(1)
    D(1) = lambda(2) / den(1)
    for l = 2..nlayer-1:
      den(l) = alpha(l) + lambda(l+1) + lambda(l)*(1-D(l-1))
      C(l) = (lambda(l+1)*(zvel(l+1)-zvel(l))
              + lambda(l)*(zvel(l-1)-zvel(l)+C(l-1))) / den(l)
      D(l) = lambda(l+1) / den(l)
    den(nlayer) = alpha(nlayer) + lambda(nlayer)*(1-D(nlayer-1))
    C(nlayer) = ((C(nlayer-1)+zvel(nlayer-1)-zvel(nlayer))*lambda(nlayer)
                 + phitop) / den(nlayer)

    pdvelm(nlayer) = C(nlayer)
    for l = nlayer-1..1:
      pdvelm(l) = C(l) + D(l) * pdvelm(l+1)

    zdvelmolvis(ig,l) = pdvelm(l) / ptimestep

参与的主题流程

主题	参与方式
热层动量扩散	`thermosphere_mod` 在 `callmolvis` 下对 `pu` 和 `pv` 各调用一次本例程，再把输出累加回 `pdu/pdv`
高层大气成分热力学	通过 `conc_mod::Akknew/cpnew/rnew` 读取随 tracer 组成更新的热导率系数、比热和气体常数

写法特点

固定格式 .F，续行混用第 6 列 & 和 $。
算法结构与 conduction.F 相似，但 alpha 不乘 cpnew，底边界从地表温度改为地表速度 velsurf=0。
pplev 声明为输入但未使用；zzlev(ig,nlayer+1) 也没有被使用，顶界面由 zlev(nlayer)+10000. 覆盖。
局部变量 m、tmean 声明后未使用。

复现要点

callmolvis 默认在 conf_phys.F 中设为 .false.，通过 getin_p("callmolvis",callmolvis) 覆盖；若 callthermos=.false. 但 callmolvis=.true.，conf_phys 会终止运行。
thermosphere_mod 的顺序是 EUV、热传导、分子黏性、分子扩散；因此 molvis 看到的 pdt 已可能包含上游 EUV 与热传导 tendency。
正常运行前需要 conc_mod::Akknew/cpnew/rnew 已由 update_r_cp_mu_ak 更新；若脱离 callthermos 路径单独调用，需要自行保证这些数组已分配并赋值。
顶层厚度硬编码为 10 km，底边界速度硬编码为 0 m/s；这两个边界条件会直接影响最高层和最低层的动量 tendency。

待确认

pplev 作为入参但未使用，可能是沿用 conduction 接口的历史残留。
zzlev(ig,nlayer+1) 被接口传入但不用于顶层厚度，实际顶界面强制为 zzlev(ig,nlayer)+10000.；是否应使用真实顶界面高度需开发者确认。
方程从 den(1)=alpha(1)+lambda(2)+lambda(1) 开始，隐含 nlayer>=2；单层垂直网格会访问 lambda(2)，复现时应避免单层配置或先确认调用约束。