conduction.F

路径

LMDZ.MARS\libf\aeronomars\conduction.F

所属目录/模块

libf\aeronomars

文件定位

conduction.F 定义 conduction_mod，是热层分子热传导求解器。子程序 conduction 用隐式三对角（Thomas 算法）在一维垂直方向上求解温度依赖的热传导方程，输出热传导温度倾向 zdtconduc（K/s）。热导率随温度按 k = Akk * T**skk（skk=0.69）变化，其中 Akk 取自 conc_mod::Akknew（随大气成分动态更新）。下边界用地表温度 tsurf，上边界用零热通量（phitop=0）。

该例程由 thermosphere_mod 在 callconduct 开关下调用，是 physiq → thermosphere → conduction 链路的一环；分子黏性模块 molvis.F 注明“Based on conduction.F”（同源算法）。

作者注：N. Descamp, F. Forget 05/1999（源码第 17 行）。

定义的符号

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• phitop（real, save）和 firstcall（logical, save）为模块内保存变量，并用 !$OMP THREADPRIVATE(phitop,firstcall)（第 67 行）声明为线程私有。
• firstcall 块（第 73–85 行）仅在首次调用执行：向标准输出打印 skk，把 phitop 置 0.0，然后翻转 firstcall。注意源码注释（第 74–77 行）说明原想同时打印 Akk,skk，但 Akk 在此阶段未定义，故只打印 skk。
• phitop 自始至终只被赋值为 0.0，从不更新——即上边界恒为零热通量，是占位的零通量边界条件。
• 无文件 I/O、无诊断输出、不修改任何模块级共享数组；pdt 虽由调用方以 inout 传入，但本例程声明为 intent(in)，只读不改。

核心逻辑

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

1. 预测温度：zt(i)=pt(ig,i)+pdt(ig,i)*ptimestep（第 89–98 行）。此时 pdt 已含上游 euvheat 注入的 zdteuv（thermosphere_mod 第 86 行先 pdt=pdt+zdteuv 再调 conduction），故热传导在线性化温度中已包含 EUV 加热。
2. 顶界面延伸：zlev(nlayer+1)=zlev(nlayer)+10000.（第 100 行）——在最高层界面之上硬编码延伸 10 km，用于定义顶层厚度。
3. 传导系数 lambda（第 102–108 行）：底层 lambda(1)=Akknew(ig,1)*tsurf**skk/zlay(1)（用地表温度、除以第 1 层全高度）；i=2,nlayer 用 lambda(i)=Akknew(ig,i)*zt(i)**skk/(zlay(i)-zlay(i-1))（用预测温度、除以层间厚度）。
4. 热容项 alpha（第 109–117 行）：muvol(i)=pplay(ig,i)/(rnew(ig,i)*zt(i))（即密度 ρ），alpha(i)=cpnew(ig,i)*(muvol(i)/ptimestep)*(zlev(i+1)-zlev(i))。顶层 alpha(nlayer) 见“待确认/复现风险”中的索引疑点。
5. Thomas 前推 C,D（第 125–143 行）：
   • 底层 den(1)=alpha(1)+lambda(2)+lambda(1)；C(1)=[lambda(1)*(tsurf-zt(1))+lambda(2)*(zt(2)-zt(1))]/den(1)；D(1)=lambda(2)/den(1)。
   • 中间层 i=2,nlayer-1：den(i)=alpha(i)+lambda(i+1)+lambda(i)*(1-D(i-1))；C(i)=[lambda(i+1)*(zt(i+1)-zt(i))+lambda(i)*(zt(i-1)-zt(i)+C(i-1))]/den(i)；D(i)=lambda(i+1)/den(i)。
   • 顶层 den(nlayer)=alpha(nlayer)+lambda(nlayer)*(1-D(nlayer-1))；C(nlayer)=[(C(nlayer-1)+zt(nlayer-1)-zt(nlayer))*lambda(nlayer)+phitop]/den(nlayer)（phitop=0 即零通量上边界）。
6. 回代求温度增量（第 151–157 行）：pdtc(nlayer)=C(nlayer)；i=nlayer-1,1,-1：pdtc(i)=C(i)+D(i)*pdtc(i+1)。
7. 输出倾向（第 164–166 行）：zdtconduc(ig,i)=pdtc(i)/ptimestep。

伪代码

conduction(ngrid,nlayer,ptimestep,pplay,pplev,pt,pdt,tsurf,zzlev,zzlay,zdtconduc):
  if firstcall:
    print skk
    phitop = 0.0
    firstcall = false

  for ig = 1..ngrid:
    zt(i)  = pt(ig,i) + pdt(ig,i)*ptimestep      # predicted T (includes EUV in pdt)
    zlay(i)= zzlay(ig,i)
    zlev(i)= zzlev(ig,i)
    zlev(nlayer+1) = zlev(nlayer) + 10000.        # artificial top, 10 km

    lambda(1) = Akknew(ig,1)*tsurf(ig)**skk / zlay(1)
    for i=2..nlayer:
      lambda(i) = Akknew(ig,i)*zt(i)**skk / (zlay(i)-zlay(i-1))

    for i=1..nlayer-1:
      muvol(i) = pplay(ig,i)/(rnew(ig,i)*zt(i))
      alpha(i) = cpnew(ig,i)*(muvol(i)/ptimestep)*(zlev(i+1)-zlev(i))
    muvol(nlayer) = pplay(ig,nlayer)/(rnew(ig,nlayer)*zt(nlayer))
    alpha(nlayer) = cpnew(ig,i)*(muvol(nlayer)/ptimestep)*(zlev(nlayer+1)-zlev(nlayer))
        # NOTE: i is still nlayer-1 here (see 待确认)

    # Thomas forward sweep
    den(1)=alpha(1)+lambda(2)+lambda(1)
    C(1)=[lambda(1)*(tsurf-zt(1))+lambda(2)*(zt(2)-zt(1))]/den(1)
    D(1)=lambda(2)/den(1)
    for i=2..nlayer-1:
      den(i)=alpha(i)+lambda(i+1)+lambda(i)*(1-D(i-1))
      C(i)=[lambda(i+1)*(zt(i+1)-zt(i))+lambda(i)*(zt(i-1)-zt(i)+C(i-1))]/den(i)
      D(i)=lambda(i+1)/den(i)
    den(nlayer)=alpha(nlayer)+lambda(nlayer)*(1-D(nlayer-1))
    C(nlayer)=[(C(nlayer-1)+zt(nlayer-1)-zt(nlayer))*lambda(nlayer)+phitop]/den(nlayer)

    # back substitution
    pdtc(nlayer)=C(nlayer)
    for i=nlayer-1..1:
      pdtc(i)=C(i)+D(i)*pdtc(i+1)

    zdtconduc(ig,i)=pdtc(i)/ptimestep

参与的主题流程

写法特点

• 固定格式 .F（续行用第 6 列 $/&）。
• phitop、firstcall 为 SAVE + THREADPRIVATE；并行复现需确认每线程都完成 firstcall 初始化。
• 热导率温度指数 skk=0.69 为 REAL, PARAMETER 硬编码（第 63 行）。
• 顶界面高度硬编码延伸 10000. m（第 100 行）。
• 局部变量 m、tmean、l 声明后未在函数体使用；入参 pplev 同样声明后未使用（见待确认）。

复现要点

• 调用前提：conduction 仅在 thermosphere_mod（callthermos 路径）内被调用，而 conc_mod::Akknew 由 update_r_cp_mu_ak 在 photochem.or.callthermos 为真时更新。因此正常运行中 Akknew 已是当步成分相关值；若试图脱离 callthermos 路径单独调用 conduction，需自行保证 Akknew/cpnew/rnew 已赋值（init_r_cp_mu 不写 Akknew）。
• 预测温度 zt 含上游 zdteuv（thermosphere_mod 先把 EUV 倾向加进 pdt 再调用本例程），复现热层倾向顺序时须保持 euvheat → conduction 的先后。
• 上边界为零热通量（phitop=0）；顶界面用硬编码 zlev(nlayer)+10000. 定义，结果依赖模式顶高度与该 10 km 假设。
• zdtconduc 在 callconduct=.false. 时由 thermosphere_mod 预置为 0，不进入 pdt。
• 顶层 alpha(nlayer) 的 cpnew 索引疑点会影响顶层热容，复现顶层热传导时按源码字面行为执行（见待确认）。

待确认

• 顶层 alpha(nlayer) 的 cpnew 索引（第 116 行）：alpha(nlayer)=cpnew(ig,i)*(muvol(nlayer)/ptimestep)*(zlev(nlayer+1)-zlev(nlayer))。此处 i 是前一个 DO i=1,nlayer-1 循环结束后的残留值 nlayer-1（Fortran 循环变量在循环结束后保留末值），故实际取 cpnew(ig,nlayer-1) 而非 cpnew(ig,nlayer)；同层 muvol/pplay/rnew/zt 均正确取 nlayer。疑似 off-by-one，待开发者确认是否应为 cpnew(ig,nlayer)。已同步记入 symbol-index.md“源码疑点与待确认符号”。
• 入参 pplev 声明为 intent(in) 但函数体未使用；是否为历史接口残留待确认。
• 局部变量 m、tmean、l 声明后未使用，疑似历史遗留。

相关页面

• aeronomars/index：aeronomars 目录总览，本文件在其“加热与传导”分类下。
• conc_mod：提供 Akknew/cpnew/rnew；update_r_cp_mu_ak 在 callthermos 路径更新这些数组。
• callkeys_mod：callconduct 开关门控本例程。
• physiq_mod：主物理调度，经 thermosphere 间接调用本例程。
• thermosphere_mod：热层物理过程总调度，直接调用方。
• molvis：分子黏性速度扩散文件页，其源码注明“Based on conduction.F”。