Mars 水平耗散链路

输入范围

dyn3d\dissip.F
dyn3dpar\dissip_p.F
dyn3d\leapfrog.F
dyn3dpar\leapfrog_p.F
dyn3d_common\inidissip.F90

Mars 运行参与度：条件经过。耗散由 apdiss 门控，按 dissip_period 在非 forward 步触发；llm==1 和 NODYN 编译选项会抑制。

链路定位

本页覆盖 Mars run 中水平耗散从初始化到施加再到守恒修正的完整数据链路，侧重各阶段之间的数据流、内部算子调用和串并行差异。apdiss 门控的 leapfrog 侧细节由 leapfrog-dissipation-conservation 覆盖；算子族本身的离散格式由差分算子总览覆盖；inidissip 的配置和校验细节由 inidissip 文件页覆盖。

链路阶段

阶段 1：inidissip 配置输出

初始化阶段由 inidissip.F90 完成，它读取 conf_gcm 的耗散配置键，通过幂迭代估算各算子最大特征值，构造垂直耗散剖面，最终写入两个 COMMON 块：

comdissipn.h 携带每层耗散时间尺度 tetaudiv(llm)、tetaurot(llm)、tetah(llm) 和全局系数 cdivu、crot、cdivh。comdissnew.h 携带算子选择逻辑 lstardis、迭代次数 nitergdiv/nitergrot/niterh 和中间参数 tetagdiv/tetagrot/tetatemp。

inidissip 还计算 dissip_period 并写回 comconst_mod，供 leapfrog 门控使用。

阶段 2：apdiss 门控

leapfrog 主循环每步重置 apdiss=.FALSE.，再按以下条件设置：

purmats 模式：MOD(itau,dissip_period)==0 且 .NOT.forward。
leapfrog 模式：MOD(itau+1,dissip_period)==0 且 .NOT.forward。
llm==1（Shallow Water）强制抑制。
#ifdef NODYN 编译选项强制抑制。

当 apdiss 为真时，leapfrog 进入耗散段。

阶段 3：耗散前准备

进入 IF(apdiss) 块后，在调用 dissip 之前依次执行：

sponge：若 callsponge 为真，先调用 sponge(ucov,vcov,teta,ps,dtdiss,mode_sponge) 施加顶层阻尼，修改进入水平耗散的风场和温度状态。
动能基准：covcont 把协变风转为逆变风，enercin 计算耗散前动能 ecin0(ij,l)。
位温转温度：tpot2t 把 teta 转为 temp，供后续 conservative 修正使用。

阶段 4：dissip 内部三路施加

dissip（串行 149 行）和 dissip_p（并行 215 行）签名相同：

SUBROUTINE dissip( vcov,ucov,teta,p, dv,du,dh )
REAL,INTENT(IN)  :: vcov(ip1jm,llm), ucov(ip1jmp1,llm), teta(ip1jmp1,llm), p(ip1jmp1,llmp1)
REAL,INTENT(OUT) :: dv(ip1jm,llm), du(ip1jmp1,llm), dh(ip1jmp1,llm)

输入是协变风 ucov/vcov、位温 teta 和层界压力 p；输出是三个增量数组（不是 tendency 率，是增量）。内部按三路执行：

第一路 grad(div)：从 comdissipn 取 tetaudiv(l)，乘以 dtdiss 得到 te1dt(l)。按 lstardis 调用 gradiv2/gradiv（或 _p 并行变体），对协变风场做 nitergdiv 次迭代，输出 gdx/gdy。在极点 ij=1..iip1 处把 gdx 置零，再把 -te1dt*gdx 加到 du、-te1dt*gdy 加到 dv。

第二路 n×grad(rot)：从 comdissipn 取 tetaurot(l)，乘以 dtdiss 得到 te2dt(l)。按 lstardis 调用 nxgraro2/nxgrarot（或 _p），对协变风场做 nitergrot 次迭代，输出 grx/gry。在极点 ij=1..iip1 处把 grx 置零，再把 -te2dt*grx 加到 du、-te2dt*gry 加到 dv。

第三路 div(grad)：从 comdissipn 取 tetah(l)，乘以 dtdiss 得到 te3dt(l)。当 lstardis 为真时，先计算 deltapres(ij,l)=AMAX1(0.,p(ij,l)-p(ij,l+1))（层压差，保证非负），然后调用 divgrad2(llm,teta,deltapres,niterh,gdx)；否则调用 divgrad(llm,teta,niterh,gdx)。把 -te3dt*gdx 加到 dh。

注意 lstardis 控制 star 变体选择：star 算子（gradiv2、nxgraro2、divgrad2）是改进版，divgrad2 额外接受 deltapres 做权重。标准版（gradiv、nxgrarot、divgrad）不接受压力权重。

阶段 5：风场更新和 conservative 修正

dissip 返回后，leapfrog 执行：

风场更新：ucov=ucov+dudis、vcov=vcov+dvdis。
tendency 转换：dudis=dudis/dtdiss、dvdis=dvdis/dtdiss，从增量转为 (m/s)/s 率，供后续诊断和输出。
conservative 修正：若 dissip_conservative 为真（串行是 SAVE 逻辑量默认 .true.，并行是 PARAMETER .TRUE.），重新计算更新后的动能 ecin，把动能损失 (ecin0-ecin) 以热力形式回注：temp=temp+(ecin0-ecin)/cpdet(temp)，再通过 t2tpot 转回位温修正 ztetaec，叠加到 dtetadis。
温度更新：teta=teta+dtetadis，再 dtetadis=dtetadis/dtdiss 转成 K/s 率。
极点温度平均：用面积加权把极点 iip1 个格点的 teta 统一为极区均值。

阶段 6：leapfrog_nogcm 排除

MARS phymars/leapfrog_nogcm.F 不设置 apdiss、不调用 sponge/dissip/covcont/enercin。所有耗散相关代码已被注释删除（ED18 标注）。它保留 dudis/dvdis/dtetadis 数组和 comdissnew.h include 仅为接口兼容；当前 COMMON 源码已删除 dyn3d/leapfrog_nogcm.F。

串并行差异

并行 halo 和 band 切换

leapfrog_p 在进入耗散段前执行 halo 请求注册和发送（Register_Hallo/SendRequest/WaitRequest），确保 ucov/vcov/teta/p 的子域边界带满足 stencil 宽度。随后通过 SetDistrib(jj_Nb_dissip) 切换到耗散专用的 band 分布，并启动 timer_dissip 计时。

OpenMP 差异

串行 dissip.F 中所有循环是普通 DO，输出数组 du/dv/dh 整体置零。并行 dissip_p.F 中循环带 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)，输出数组按 ij_begin:ij_end 子域置零。极点零化条件也不同：串行无条件清除南北两极 ij=1..iip1 和 ij=ip1jm+1..ip1jmp1；并行按 pole_nord/pole_sud 逻辑量条件清除。

conservative 修正的并行 halo

并行版在 conservative 修正前需要额外一次 halo 交换：更新后的 ucov/vcov 必须同步边界才能重新计算 ecin。串行版不需要这次通信。

dissip_conservative 声明差异

串行版中 dissip_conservative 是 SAVE 逻辑量，初始化为 .true.，理论上可在运行时修改。并行版中是 PARAMETER，编译期固定为 .TRUE.。

输入输出数据流

inidissip → comdissipn{tetaudiv,tetaurot,tetah,cdivu,crot,cdivh}
          → comdissnew{lstardis,nitergdiv,nitergrot,niterh,...}
          → comconst_mod{dissip_period}

leapfrog apdiss gate ← dissip_period, itau, forward, llm

sponge(ucov,vcov,teta,ps) → 修改 ucov,vcov,teta
covcont(ucov,vcov) → ucont,vcont
enercin(vcov,ucov,vcont,ucont) → ecin0

dissip(vcov,ucov,teta,p) ← comdissipn, comdissnew, comconst_mod{dtdiss}
      ├─ gradiv/gradiv2(ucov,vcov,nitergdiv) → gdx,gdy → du,dv
      ├─ nxgrarot/nxgraro2(ucov,vcov,nitergrot) → grx,gry → du,dv
      └─ divgrad/divgrad2(teta,deltapres,niterh) → gdx → dh

ucov += dudis; vcov += dvdis
dudis /= dtdiss; dvdis /= dtdiss  → (m/s)/s

if dissip_conservative:
  covcont + enercin → ecin
  temp += (ecin0-ecin)/cpdet(temp)
  t2tpot → ztetaec
  dtetadis += ztetaec - teta

teta += dtetadis
dtetadis /= dtdiss  → K/s
极点 teta 面积加权平均

复现要点

复现耗散链路必须先完成 inidissip，否则 comdissipn 和 comdissnew 中的系数和开关未初始化。
dissip_period 决定施加频率：Mars 常用值来自 inidissip 自动计算，而非用户直接指定。
lstardis 为真时走 star 算子路径，divgrad2 需要 deltapres 做压力权重；为假时走标准路径，divgrad 不需要。
增量和 tendency 的区分：dissip 输出增量（乘过 dtdiss），leapfrog 加到风场后才除以 dtdiss 得到率。复现诊断输出时必须使用转换后的率。
conservative 修正是动能→热能的单方向注入，不会改变风场，只改变 temp/teta。
并行的 halo 交换是正确性关键：遗漏 stencil 边界同步会导致算子边界计算错误，且误差随迭代次数放大。

待确认

lstardis 在 Mars 当前 COMMON/MARS 配置下的实际值需从 conf_gcm/run.def 追踪确认。
sponge 的 callsponge 和 mode_sponge 在 Mars 运行中的默认值需在 sponge 页面补齐。
dissip_p 内部 pole_sud 条件下 ije=ij_end-iip1 的范围缩减对 dv 零化与施加的影响，需与串行 ij+ip1jm 处的零化做几何对齐确认。