LMDZ.MARS Wiki

conduction.F

快速理解

它做什么： 热层分子热传导求解器。被 thermosphere_mod 在 callconduct 时调用。

基本过程： 隐式三对角求解温度依赖的热传导方程（k=Akk*T^0.69）→ 下边界 tsurf、上边界零通量。

关键结果： 热传导温度倾向 zdtconduc（K/s），累加回 pdt。

路径

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 行）。

定义的符号

符号 类型 行号 作用
conduction_mod module 1（END MODULE 172） 包装热传导例程的模块
conduction subroutine 7（END SUBROUTINE 170） 隐式三对角求解垂直分子热传导倾向

依赖的模块

use 模块 only 列表 用途 待确认
conc_mod Akknew, rnew, cpnew 取逐格点逐层热导率系数 Akknew、比气体常数 rnew、定压比热 cpnew 构造传导系数与热容项 -

调用的关键例程

被调用例程 所在模块/文件 调用位置 作用
无外部例程 - - 本例程只做局部数组运算和 Thomas 递推，不调用其他子程序（仅用 Fortran intrinsics）

输入

输入 来源 类型/维度 单位 含义
ngrid thermosphere_mod integer - 大气列数
nlayer thermosphere_mod integer - 垂直层数
ptimestep thermosphere_mod real s 物理时间步
pplay thermosphere_mod real (ngrid,nlayer) Pa 层中气压，用于数密度/密度 muvol=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 已累计温度倾向；用于构造预测温度 zt=pt+pdt*ptimestep
tsurf thermosphere_mod real (ngrid) K 地表温度，作为下边界
zzlay thermosphere_mod real (ngrid,nlayer) m 层中高度
zzlev thermosphere_mod real (ngrid,nlayer+1) m 层界面高度
Akknew / rnew / cpnew conc_mod real (ngrid,nlayer) 热导率系数 / J/kg/K / J/kg/K 成分依赖的热力学属性（由 update_r_cp_mu_ak 在 photochem/callthermos 路径更新）

输出

输出 去向 类型/维度 单位 含义
zdtconduc thermosphere_mod（第 93 行 pdt=pdt+zdtconduc） real (ngrid,nlayer) K/s 分子热传导温度倾向；callconduct=.false. 时由调用方预置为 0

共享状态与副作用

核心逻辑

逐列 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

参与的主题流程

主题 参与方式
热层加热/冷却 thermosphere_mod 在 callconduct 下调用，把热传导倾向累加回 pdt，与 EUV 加热、分子黏性/扩散并列
大气成分热力学 通过 conc_mod::Akknew/cpnew/rnew 读取随 tracer 组成动态更新的热导率、比热和气体常数

写法特点

复现要点

待确认

相关页面