mucorr.F

路径

LMDZ.MARS\libf\phymars\mucorr.F

所属目录 / 模块

libf\phymars

文件定位

mucorr.F 定义 mucorr_mod，只提供 mucorr 一个例程。它在不显式计算日变化循环时，根据太阳赤纬和各格点纬度，计算日平均等效太阳入射角余弦 pmu 以及日照时间比例 pfract。当前源码树中，physiq_mod.F 在 diurnal=.false. 分支调用它，把输出写入共享太阳几何数组 mu0/fract，随后供辐射和相关地表/尘埃流程读取。

本文件不 use 其他模块，也不调用其他自定义例程；核心依赖是 Fortran 内置三角函数、平方根和反正切。源码末尾还对 pmu 做一次“大气厚度/行星半径”球面修正，但活动代码使用硬编码常数 1224. 和 35.，而不是同一段中已计算但被注释掉的 alph=phaut/prad 公式。

定义的符号

依赖的模块

调用的关键例程

上游调用点

输入

输出

共享状态与副作用

• 本文件没有 module 变量、SAVE 变量、common block、文件 I/O 或诊断输出。
• mucorr 只通过 intent(out) 参数写出 pmu/pfract。
• 上游 physiq_mod.F 把输出放入 comsaison_h::mu0/fract 后，下游短波辐射、NIR CO2 加热、坡面/尘埃相关流程会读取这些太阳几何量。

核心逻辑

1. 计算 pi，把 pdeclin 记为 z，预先得到 cz=cos(z) 和 sz=sin(z)。
2. 对每个纬度点，计算 sin(phi)、cos(phi)、tan(phi)，其中 cos(phi) 小于等于 1.e-9 时被钳制到 1.e-9。
3. 由 t=-tan(phi)*sin(declin)/cos(declin) 和 a=1-t*t 推导日照半角 tp；当 a<0 时先保留原始 ap，再把用于平方根的 a 钳制为 0。
4. 用 pmu=(sin(phi)*sin(declin)*t)/pi + cos(phi)*cos(declin)*sin(t)/pi 和 pfract=t/pi 计算日平均量。
5. 若原始 ap<0，进入极昼/极夜分支，把 pmu 设为 sin(phi)*sin(declin)、pfract 设为 1。
6. 将非正 pmu 钳制为 0，再用 pfract 归一化 pmu；若归一化后 pmu 为 0，则把 pfract 也置为 0。
7. 计算 alph=phaut/prad，但活动修正公式对每个格点执行 pmu=sqrt(1224.*pmu*pmu+1.)/35.；基于 alph 的公式保留在注释中。

伪代码

mucorr(npts, declin, lat, pmu, pfract, phaut, prad):
    pi = 4 * atan(1)
    cz = cos(declin)
    sz = sin(declin)

    for each grid point j:
        phi = lat[j]
        cphi = max(cos(phi), 1.e-9)
        sphi = sin(phi)
        t = -(sphi / cphi) * sz / cz
        a = 1 - t*t
        ap = a

        if t == 0:
            daylight_half_angle = pi / 2
        else:
            daylight_half_angle = atan_branch(sqrt(max(a, 0)) / t)

        pmu[j] = daily_mean_cosine(daylight_half_angle, sphi, sz, cphi, cz)
        pfract[j] = daylight_half_angle / pi

        if ap < 0:
            pmu[j] = sphi * sz
            pfract[j] = 1

        pmu[j] = max(pmu[j], 0)
        pmu[j] = pmu[j] / pfract[j]
        if pmu[j] == 0:
            pfract[j] = 0

    alph = phaut / prad
    for each grid point j:
        pmu[j] = sqrt(1224 * pmu[j]^2 + 1) / 35

参与的主题流程

写法特点

• 固定格式 Fortran .F 文件，但内容是一个 module 加一个 contained subroutine。
• tz 局部变量声明后未在当前源码中使用。
• alph=phaut/prad 被计算，但活动代码没有使用 alph；实际生效的是硬编码 sqrt(1224.*pmu^2+1.)/35.。
• 极昼/极夜判断依赖 ap=1-t*t 的原始符号；后续用于平方根的 a 会被截断到非负。

复现要点

• plat 和 pdeclin 在源码中按弧度参与三角函数；不要按角度制复现。
• cphi 有下限 1.e-9，但 cz=cos(pdeclin) 没有显式下限；移植时要保留源码行为，或把新增保护标为行为变更。
• 活动球面修正必须使用硬编码 1224. 和 35.；注释中的 alph 公式不是当前运行路径。
• pmu 先除以 pfract，之后才在 pmu==0 时把 pfract 置零；如果改写为先清零 pfract，会改变边界条件行为。

待确认

• 源码没有注释 1224. 和 35. 的来源；从形式看它们替代了 phaut/prad 参数化修正，但本页只按源码事实记录，不推断常数物理来源。
• cos(pdeclin)=0 附近的数值行为没有源码保护；需结合实际轨道参数范围确认是否可能触发极端分母。

相关页面

• phymars/index：phymars 目录总览。
• comsaison_h：mu0/fract/local_time 共享太阳几何数组。
• nirco2abs：近红外 CO2 吸收；说明非日变化时单独日内积分太阳角。
• radiation：辐射主题页，mu0/fract 是辐射入口参数之一。