molvis.F

路径

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。

定义的符号

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• 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 独立求解：

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

参与的主题流程

写法特点

• 固定格式 .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)，复现时应避免单层配置或先确认调用约束。

相关页面

• aeronomars/index：aeronomars 目录总览，本文件在“分子扩散”组中。
• conduction：同源的热层分子热传导求解器，源码注释称 molvis.F 基于该文件。
• conc_mod：提供 Akknew/cpnew/rnew，本例程用它们构造黏性扩散系数和密度。
• callkeys_mod：定义 callmolvis 开关。
• conf_phys：读取 callmolvis 并检查其依赖 callthermos。
• thermosphere_mod：热层物理过程总调度，直接调用方。