solarlong.F

路径

LMDZ.MARS\libf\phymars\solarlong.F

所属目录 / 模块

libf\phymars

文件定位

solarlong.F 定义 solarlong_mod 模块和 solarlong(pday, psollong) 子程序，把以北半球春分为零点的火星日期 pday 转换为太阳经度 Ls。它读取 planete_h 中由 iniorbit 派生的轨道参数，先由日期相对近日点得到平均近点角，再用 Newton 迭代解 Kepler 方程，最后加上近日点太阳经度 lsperi 并把结果折回 [0,2*pi]。

定义的符号

符号	类型	行号	作用
`solarlong_mod`	module	1	火星日期到太阳经度转换模块。
`solarlong`	subroutine	7	输入 sol 日期 `pday`，输出弧度制太阳经度 `psollong`。

依赖的模块

use 模块	only 列表	用途	待确认
`planete_h`	`lsperi`, `peri_day`, `year_day`, `e_elips`	提供近日点日期、年长、轨道偏心率和近日点太阳经度。	-
`comcstfi_h`	`pi`	构造 `2*pi` 周期和角度归一化。	-

调用的关键例程

被调用例程	所在模块 / 来源	调用位置	作用
`nint`	Fortran intrinsic	33	把日期相对近日点的年周期间隔折到最近周期。
`sin`, `cos`, `abs`, `sqrt`, `tan`, `atan`	Fortran intrinsic	34-52	Newton 迭代和偏近点角到真近点角转换。

下游调用点

调用方	源码位置	用途
`physiq_mod.F`	`use solarlong_mod` 在 65 行；调用在 879、881 行	在主物理步中由 `zday` 或 `day_ini` 计算 `comsaison_h::zls`，供后续轨道几何、辐射和季节依赖过程使用。
`dyn1d\testphys1d.F90`	170 行注释调用	历史或测试路径注释中展示由 `day` 计算 `zls` 的用法。
`dynphy_lonlat\phymars\lect_start_archive.F`	745 行注释调用	归档读取路径中保留按时间列表计算太阳经度的注释调用。

输入

输入	来源	类型 / 维度	单位	含义
`pday`	调用方	real scalar	sol	火星日期；注释说明 `pday=0` 对应北半球春分。
`peri_day`	`planete_h`	real saved scalar	sol	近日点日期。
`year_day`	`planete_h`	real saved scalar	sol	火星年长。
`e_elips`	`planete_h`	real saved scalar	-	轨道偏心率。
`lsperi`	`planete_h`	real saved scalar	rad	近日点太阳经度。

输出

输出	去向	类型 / 维度	单位	含义
`psollong`	调用方，通常为 `comsaison_h::zls`	real scalar	rad	太阳经度 `Ls`，经简单边界处理后落在 `[0,2*pi]` 附近。

共享状态与副作用

本模块没有 SAVE 状态、文件 I/O 或标准输出。
计算结果依赖 planete_h 的当前轨道状态；调用前需要确保 iniorbit 已初始化 e_elips/lsperi。
Newton 循环最多 10 次，收敛阈值为 abs(zdx)<=1.e-7。

核心逻辑

计算相对近日点的年周期比例 zz=(pday-peri_day)/year_day。
用 zanom=2*pi*(zz-nint(zz)) 得到平均近点角，并取绝对值 xref。
以 xref + e_elips*sin(xref) 为初值，解 zx0 - e*sin(zx0) = xref。
若 zanom<0，把偏近点角 zx0 改为负号。
用椭圆轨道关系把偏近点角转换为真近点角 zteta。
psollong=zteta+lsperi，若小于 0 加 2*pi，若大于 2*pi 减 2*pi。

伪代码

zz = (pday - peri_day) / year_day
zanom = 2*pi * (zz - nearest_integer(zz))
xref = abs(zanom)

zx0 = xref + e_elips*sin(xref)
repeat up to 10:
    zdx = -(zx0 - e_elips*sin(zx0) - xref) / (1 - e_elips*cos(zx0))
    if abs(zdx) <= 1e-7: exit
    zx0 = zx0 + zdx
zx0 = zx0 + zdx
if zanom < 0: zx0 = -zx0

zteta = 2*atan(sqrt((1+e_elips)/(1-e_elips))*tan(zx0/2))
psollong = zteta + lsperi
wrap once into [0, 2*pi]

参与的主题流程

主题	参与方式
季节推进	`physiq` 由当前 `zday` 计算 `zls`，作为季节相位的主入口。
轨道 / 太阳几何	`solarlong` 输出的 `zls` 下游进入 `orbite`，计算日火距离和太阳赤纬。
辐射与水/CO2 循环	`zls` 影响辐射、极冠、尘暴、坡面和季节性循环路径的判定。

写法特点

固定格式 Fortran 文件，使用模块封装。
Newton 收敛后无论是否因阈值退出，都会再执行一次 zx0=zx0+zdx；复现时应保留这一字面行为。
结果只做一次加/减 2*pi 的边界处理，而不是通用 modulo 归一化；当前输入通常接近一个火星年内日期，足够覆盖主路径。

复现要点

调用前必须先初始化 planete_h 的轨道派生量，通常由 tabfi 或 1D 初始化路径调用 iniorbit 完成。
pday 单位是 sol，不是秒；year_day/peri_day 必须使用同一单位。
psollong 输出是弧度；若页面或配置使用 degree，需要显式转换。

待确认

结果归一化只处理一次越界；若调用方传入远离当前年周期的极端 pday，是否仍满足预期，需要结合上游日期范围确认。