physiq_mod.F

路径

LMDZ.MARS\libf\phymars\physiq_mod.F

所属目录/模块

libf/phymars
MODULE physiq_mod

文件定位

physiq_mod.F 是 Mars physics package 的主时间步调度文件。它定义 physiq_mod 模块和唯一入口 physiq，由动力-物理接口 call_physiq 以及 1D 测试程序直接调用，接收动力核给出的压力、温度、风、tracer 和垂直质量通量，返回物理过程对风、温度、tracer 与表压的 tendency。源码注释在 physiq_mod.F:163-170 明确说明本文件组织 LMD Mars GCM 的物理参数化。

本文件不是某个单一物理方案，而是把 firstcall 初始化、restart 读取/写出、太阳几何、辐射、尘埃输运、垂直扩散、热羽流、对流、重力波、水循环、CO2 云/凝结、化学、热层、地表/土壤更新和诊断输出串成一个物理时间步。下游单个过程的内部算法见相关文件页；本页重点记录 physiq 如何按源码顺序连接这些过程、累加 tendency、更新共享状态并产生副作用。

定义的符号

day_ini/icount/time_phys 与两个检查开关都声明为 OpenMP THREADPRIVATE，说明每个 OpenMP 线程保存自己的局部记忆状态。

依赖的模块

physiq 的 use 列表非常长，源码位于 physiq_mod.F:16-158。为了复现调度逻辑，按用途分组如下；每一组都在调用链或共享状态更新中被直接使用。

直接调用方

调用的关键例程

下表按源码执行顺序列出主路径调用点。write_output 和 wstats 字段很多，本页只记录分段入口，具体字段见“诊断与输出”。

输入

physiq 的参数方向在 physiq_mod.F:265-291 声明。

pplev 和 pplay 是动力输入，源码在 physiq_mod.F:445-447 明确要求 physics 不得修改它们，因此 physiq 复制到工作数组 zplev/zplay，CO2 凝结后的压力坐标更新只改工作副本。

输出

这些输出在每步开始统一清零（physiq_mod.F:844-848），随后各物理过程按源码顺序做加法累积。地表 tracer、地表温度、土壤、XIOS/diagfi/stats/restart 等副作用不通过这五个输出返回，而是直接写共享模块或文件。

共享状态与副作用

初始化、restart 与诊断行为

firstcall 初始化

IF (firstcall) 从 physiq_mod.F:602 开始。非 MESOSCALE 路径调用 phyetat0("startfi.nc",...) 读取物理初始状态；MESOSCALE 路径假定模块变量已由外部设置，只打印检查信息、设置 day_ini=pday 并填充一些默认值。若非 MESOSCALE 且 startphy_file 为假，源码走 academic 初始化：给土壤层深度、tsurf/tsoil 和热惯量设置默认值，并把 day_ini=pday。

firstcall 之后会检查 pday 与 day_ini 是否同步（713-718），初始化 tracer 和地表（725, 729），按 callsoil 初始化土壤（733-746），初始化热层/电离层、平均热力属性、NLTE/NIR 表（763-775），写 restart 初始控制信息（789/795），准备 topmons、尘埃凝并、ATKE 和 XIOS 输出（804-823）。

每步初始化

每个物理步都会更新 XIOS 时间（837-839），清零所有输出 tendency、地表 tendency、通量和若干诊断数组（842-856），计算坡面到网格平均状态（857-858），设置 IRtoVIScoef=2.6 初值（860-862），计算 zday、local_time 和 zls（870-882），并把动力压力输入复制到 zplev/zplay（884-888）。

restart 写出

非 MESOSCALE 路径在诊断前计算 write_restart：ecritstart>0 且 MODULO(icount*iphysiq,ecritstart)==0 时写多时间 restart，lastcall 时一定写（2710-2721）。写出时间 ztime_fin 对 DYNAMICO/unstructured 与 LMDZ 规则不同（2723-2749），最后调用 physdem1("restartfi.nc",...) 写 tsurf/tsoil/inertiesoil/albedo/emis/q2/qsurf/qsoil/tauscaling/totcloudfrac/wstar/watercap/perennial_co2ice 等状态（2754-2758）。

诊断与输出

诊断先构造更新后工作态 zt/zu/zv/zq（2565 之后；水/CO2 调用点页也记录了 zq=pq+pdq*ptimestep 的诊断用法），再计算密度、位温、辐射通量和各主题柱积分。wstats 在非 MESOSCALE 路径累计 ps/tsurf/radiation/temp/u/v/w/rho/pressure/q2/water/co2/dust/chemistry 等字段，并在 lastcall.and.callstats 时调用 mkstats（2973-3248）。MESOSCALE 路径只填 comm_wrf 缓存（3260-3297）。常规 diagfi/XIOS 输出通过大量 write_output 调用完成（3324-4088），末尾还输出 CO2 总量守恒诊断 co2conservation（4110-4134）并发送 XIOS 固定字段（4140-4145）。

核心流程

1. 步首清零与时间/网格准备：输出 tendency 和地表增量全部清零，写 local_time/zls，复制动力压力为 zplev/zplay，可选修正动力 tracer 质量非守恒，并计算 zzlay/zzlev/zh/pw 等工作变量。
2. 辐射准备与辐射倾向：callrad 分支计算轨道几何与太阳角，按 callnlte/nltemodel/callnirco2/callslope/CLFvarying 调用 NLTE、LTE 辐射、坡面通量和近红外 CO2 加热，最终把 dtrad 累加到 pdt。
3. 尘埃与地形相关过程：按 rdstorm -> topflows -> dust_coagulation -> dustinjection 的顺序更新尘埃/风暴尘埃/topdust tendency 或 dustliftday；然后地形重力波 calldrag_noro 把 zdugw/zdvgw/zdtgw 累加到风温 tendency。
4. 边界层、热羽流、对流与非地形 GW：vdifc 用辐射+地热通量、当前累计温度 tendency 和地表状态计算垂直扩散 tendency；热羽流、干对流和非地形 GW 随后继续往 pdu/pdv/pdt/pdq 累加。
5. tracer 物理：水、CO2 云、尘卷风、沉降、化学、热层：水云和 CO2 云分别在第 9a/9a bis 段处理，callsedim 处理尘埃和水冰沉降，calchim 可读取云过程输出，thermosphere 在高层分支追加 tendency。
6. CO2 凝结与坐标更新：co2condens 被源码标注为最后一个大气物理过程（2296-2298），输出 pdpsrf、风温/tracer/地表 tendency。非 MESOSCALE 路径随后用 pdpsrf 更新 ps/zplay/zplev/zzlay/zzlev（2357-2403）。
7. 地表/土壤和非负保护：累加 qsurf += ptimestep*dqsurf，tsurf += ptimestep*zdtsurf，按水冰/CO2 霜条件调整反照率/发射率，再调用 soil；随后对风暴尘埃/背景尘埃做非负保护并更新 watercap。
8. 输出、restart、场检查和计数：计算诊断工作态，按 stats/MESOSCALE/diagfi/XIOS 写出，末尾可检查输出场并 icount=icount+1。

伪代码

physiq(args):
  pdq = 0
  if firstcall:
    read check_physics_inputs/outputs
    initialize radiation arrays and optional XIOS
    if not MESOSCALE:
      phyetat0(startfi.nc) -> tsurf/tsoil/albedo/emis/qsurf/q2/...
      if no startphy_file: build academic defaults
    else:
      check externally prepared state and set MESOSCALE defaults
    require pday == day_ini
    initracer; surfini
    if callsoil: optionally waterice_tifeedback, then soil(firstcall)
    initialize thermosphere/NLTE/NIR/mean gas properties/topmons/coagulation/ATKE/output

  if check_physics_inputs: check pt/pu/pv/pplev/pq

  zero pdu/pdv/pdt/pdq/pdpsrf, zdtsurf, dqsurf, diagnostic fluxes
  compute mesh averages, zday, local_time, zls
  copy pplev/pplay into zplev/zplay; ps = zplev(:,1)
  if call_mass_fixer_dyn: pdq += tracer_mass_fixer_dyn(...)
  if photochem or callthermos: update_r_cp_mu_ak(...)
  compute zzlay/zzlev, potential temperature, vertical velocity
  compute initial CO2 total for conservation diagnostic

  if callrad:
    orbite -> dist_sol/declin
    solang or mucorr -> mu0/fract
    optional NLTE cooling and daily nlthermeq
    callradite -> radiative tendencies and fluxes
    optional clear-sky callradite for CLFvarying
    optional slope radiation and NIR CO2 absorption
    dtrad = blendrad(...) or sum SW/LW/NIR/NLTE
    pdt += dtrad

  if rdstorm: rocketduststorm; add pdqrds to dust/stormdust pdq
  if topflows: topmons; add pdqtop to topdust and dust pdq
  if dust_coagulation: dust_coagulation_main; add pdqcoag
  if dustinjection > 0: compute_dtau -> dustliftday
  if calllott: calldrag_noro; add wind/temp tendencies
  if calldifv: vdifc; add vertical diffusion tendencies and surface fluxes
  if calltherm: calltherm_interface; add thermals tendencies
  if calladj: convadj; add dry-convection tendencies
  if calllott_nonoro: nonoro_gwd_ran; optional nonoro_gwd_mix

  if water: watercloud; inject water/HDO/CCN/dust tendencies with guards
  if co2clouds: co2cloud; inject CO2-cloud tendencies with guards
  if callddevil: dustdevil; add dustdevil tendencies
  if sedimentation: callsedim; add atmospheric and surface sedimentation tendencies
  if photochem: surfacearea; calchim; add chemical tendencies
  if callthermos: thermosphere; add high-atmosphere tendencies

  if tituscap: geticecover -> qsurf_tmp(co2)
  if callcond:
    co2condens -> zdtc/zdtsurfc/zduc/zdvc/zdqc/zdqsc/pdpsrf
    if ngrid == 1: scale selected outputs by CO2cond_ps
    add CO2 condensation tendencies
    if not MESOSCALE: update ps, zplay, zplev, zzlay, zzlev

  qsurf += ptimestep*dqsurf
  tsurf += ptimestep*zdtsurf
  if water: update frost albedo/emissivity, watercap/refill_watercap
  if callsoil: waterice_tifeedback if enabled, then soil
  protect dust/stormdust from negative updated values

  build updated diagnostic state zt/zu/zv/zq
  write optional restartfi, stats, MESOSCALE comm_wrf, diagfi/XIOS fields
  compute and write CO2 conservation diagnostic
  if check_physics_outputs: check zt/zu/zv/zplev/zq
  icount += 1

参与的主题流程

写法特点

• 固定格式 Fortran + 大型调度例程：文件主体是一个 4000 行级别的固定格式 Fortran 子程序，局部数组承担跨物理段的工作态。
• 加法式 tendency 总线：每个物理过程返回局部 tendency，随后显式 pdt += ...、pdq += ...、pdu/pdv += ...。复现时要保留源码顺序。
• 动力输入不原地修改：pplev/pplay 被复制为 zplev/zplay，CO2 凝结后的压力坐标只更新工作副本。
• firstcall 与每步初始化严格分离：表文件、restart、土壤/NLTE/NIR/ATKE/XIOS 初始化只在 firstcall；tendency、时间、太阳几何和工作压力每步重置。
• 大量 compile-time 分支：MESOSCALE、CPP_XIOS、DUSTSTORM、MESOINI 会改变 restart、输出、共享缓存和某些初始化行为。
• 成对非负保护：水云、CO2 云和尘埃相关质量/数浓度 tracer 多处使用 where 把 mass/number 成对钳到 -pq/ptimestep + 1.e-30。
• 诊断基于更新后工作态：水、CO2、尘埃、化学等柱积分和输出使用 zq = pq + pdq*ptimestep 这类更新后副本，而不是裸 pq。

复现要点

• 按源码顺序累加 pdu/pdv/pdt/pdq/pdpsrf/dqsurf/zdtsurf；不要把辐射、水云、沉降、CO2 凝结或土壤更新重排。
• co2condens 必须位于最后一个大气物理过程位置；它之后的 ps/zplay/zplev/zzlay/zzlev 更新依赖 pdpsrf。
• callradite -> watercloud -> callsedim 的顺序决定水冰半径 rice 和沉降输入；相关细节见 physiq-water-cycle-callpoints。
• co2cloud -> co2condens 之间通过 zdqssed_co2 与 zcondicea_co2microp 传递 CO2 云沉降/凝结信息；相关细节见 physiq-co2-cycle-callpoints。
• 对 pplev/pplay 的任何后处理都应写到 zplev/zplay；源码注释明确禁止 physics 修改动力输入。
• 在非 MESOSCALE 与 MESOSCALE 分支中 restart 和输出路径不同；普通 GCM 写 restartfi.nc/stats/diagfi/XIOS，MESOSCALE 写 comm_wrf 缓存。
• 若开启 check_physics_inputs 或 check_physics_outputs，失败会经 check_fields_mod 打印诊断并中止；默认值在本文件为 .false.。

待确认

• physiq_mod.F:2957 的 icetotco2(ig) = icetot(ig) + ... 疑似把水冰柱积分 icetot 用作 CO2 冰累加基值；已有 physiq-co2-cycle-callpoints 和 co2-cycle 标注为疑似 bug，仍需上游版本或开发者确认。
• callradite_mod.F 已有 callradite_mod 文件页，lwmain_mod.F 已有 lwmain_mod 文件页，swmain_mod.F 已有 swmain_mod 文件页，vdifc_mod.F 已有 vdifc_mod 文件页。本页只记录它们在 physiq 的调用位置，不展开其完整内部算法。
• calchim 对 zdqcloud/zdqscloud 的具体非均相化学使用需要在 aeronomars 化学文件页中继续确认。
• MESOSCALE 编译路径跳过部分普通 GCM restart/坐标更新逻辑；本页按源码分支记录，未复现外部 WRF/mesoscale 驱动如何消费所有 comm_* 字段。

相关页面

• phymars 目录
• callphysiq_mod / dynphy_lonlat 目录
• physiq water-cycle callpoints
• physiq CO2-cycle callpoints
• radiation
• dust-cycle
• water-cycle
• co2-cycle
• callkeys_mod
• comsaison_h
• check_fields
• comm_wrf
• wstats_mod
• write_output_mod