全空间大地测量全要素固体潮效应算法Fortran代码

[计算目标]

由给定经纬度大地高和时刻，预报地面及地球外部空间计算点10种大地测量要素的固体潮效应。改善IERS2010站位移固体潮算法，实现全空间大地测量全要素固体潮效应统一解析计算。天体的地球引潮位计算时，月球取6阶、太阳取3阶和太阳系地球外部行星取2阶。

兼容IERS2010协议地球重力位与地面站点位移固体潮效应的数值标准与算法，顾及体潮勒夫数的纬度依赖和频率相关性，严密实现地面及固体地球外部全要素几何和物理大地测量固体潮效应算法统一，以维持各种大地测量要素固体潮效应之间的解析关系。

[输出结果]

大地测量固体潮效应tdn(14)

tdn(1:14)存放10种大地测量要素固体潮效应：高程异常（大地水准面mm）tdn(1)、地面重力⊙（μGal）tdn(2)、扰动重力（μGal）tdn(3)、地倾斜⊙（SW南向/西向mas）tdn(4) tdn(5)、垂线偏差（SW南向/西向mas）tdn(6) tdn(7)、水平位移⊙（EN东向/北向mm）tdn(8) tdn(9)、地面径向⊙（大地高mm）tdn(10)、地面正（常）高⊙（mm）tdn(11)、扰动重力梯度（径向10μE）tdn(12)和水平重力梯度（NW北向/西向10μE）tdn(13) tdn(14)。

计算点可以位于地面、低空、卫星、海洋或水下空间，上述标注⊙的要素，当且仅当计算点位与地球固连时有效。

[3个输入地球物理模型文件]

（1）JPL月球和行星星历文件JEPH.440。JPL月球和行星星历DE440/LE440起止时间1850年至2250年。

（2）IERS地球定向参数EOP产品。IERSEOP_C04格式文件IERSeopc04.dat。

（3）勒夫数频率相关系数文件frqadjlovekhl.txt。通过改造IERS2010协议标准中的表6.5a、6.5b、6.5c、7.2、7.3a和7.3b，增加Cartwright-Tayler相应分潮平衡潮高全球最大振幅（10⁻⁵m）生成。

[测试入口程序]

    Solidtidaleffects.f90

    输出文件reslt.txt记录：长整型ETideLoad格式时间，相对开始时间的天数（实数），tdn(1:14)

[6个核心模块]

（1）大地测量全要素固体潮效应计算模块

    GravFdSldTd(td, BLH, tdn, eoput1, GRS, val, NUM, frd, TRStCRS)

    输入：td, BLH(3)－计算时刻（儒略日）td (JD)= MJD+2400000.5 + 32.184/86400，计算点的纬度、经度（度小数）和大地高（m）

    输入：eoput1(6)－计算历元时刻td的地球定向参数EOP[x(") y(") UT1-UT(s) LOD(s) dX(") dY(")]

    输入：GRS(6)－gm, ae, j2, omega, 1/f, 缺省值

    输入：frd(119,8)－勒夫数频率相关系数

    输入：NUM－=645 DE440，=156 DE405

    输入：TRStCRS(3,3)－ITRS到GCRS转换矩阵

    返回：tdn(1:14)－全要素固体潮效应。

（2）N体(含日月)理论固体潮位系数计算模块

    NormTide(td, val, nn, cnm, snm, rln, GRS, TRStCRS)

    输入：rln(3)－计算点的球坐标（ITRS）

    返回：cnm, snm－全部引潮天体的地球位系数直接影响，即全部天体的引潮位（生成位TGP）对应的位系数变化。

（3）用标称勒夫数的全要素固体潮效应计算模块

    solidtdnormal(rln, 6, cnm, snm, tdn, GRS)

    输入：cnm, snm-理论固体潮位系数

    返回：tdn(1:14)-用标称标称勒夫数（代码内部赋值）计算的全要素固体潮效应，直接影响与间接影响之和。

（4）位勒夫数频率相关性位系数校正模块

    Freqdploven(td,cnm,snm,frd,eop)

    输入：frd(119,8)－勒夫数频率相关系数

    返回：cnm, snm-勒夫数频率相关性导致的位系数变化

（5）全要素固体潮效应勒夫数频率相关性校正模块

    Tdlovefreqadj(rln, cnm, snm, ftd, GRS)

    输入：cnm, snm-勒夫数频率相关性导致的位系数变化

    返回：ftd(1:14)-全要素固体潮效应勒夫数频率相关性校正。

（6）扰动地球重力场球谐综合计算模块

    RentGField(maxn, rln, cnm, snm, gvm, GRS)

    返回：gvm(9)－扰动位（m2/s2），高程异常（mm）, 空间异常、扰动重力（μGal），垂线偏差（南向、西向，mas），径向重力梯度（10E），水平重力梯度（北向、西向，10μE）。

[主要辅助模块]

（7）JPL 行星星历读取软件包

    ReadDE405.f90 !改化，适合DE405/DE440。

（8）正常重力场计算模块

    normdjn(GRS, djn)；GNormalfd(BLH, NFD, GRS)

    返回NFD(5)－正常重力位，正常重力，正常重力梯度，正常重力线方向，正常梯度方向

（9）规格化连带勒让德函数及其一、二阶导数模块

    BelPnmdt(pnm, dpt1, dpt2, maxn, t)

    计算规格化连带勒让德函数Pnm及其对θ一、二阶导数dpt1, dpt2，t=cosθ。

    规格化Pnm采用改进的Belikov递推算法，一、二阶导数采用非奇异递推算法。

（10）勒让德函数及导数计算模块

    PlmBar_d(p, dp, lmax, rlat)

    计算勒让德函数p(t)及其对ψ一阶导数dp，t=cosψ。

（11）大地坐标形式变换包

    BLH_RLAT(GRS, BLH, RLAT)；BLH_XYZ(GRS, BLH, XYZ)

    RLAT_BLH(GRS, RLAT, BLH)

（12）时间系统转换包

    CAL2JD (IY0, IM0, ID0, DJM, J)；JD2CAL(DJ1, DJ2, IY, IM, ID, FD, J)

（13）IAU SOFA2000库

    TRStoCRS(tm, eoput1, TRStCRS)；GMST2000(UTA, UTB, TTA, TTB)；

    RBPN2000(X, Y, S, RBPN)；……

（14）其他辅助模块

    PickRecord(str0, kln, rec, nn)；tmcnt(tm, iyr, imo, idy, ihr, imn, sec)

    mjdtotm(mjd0, ltm)；tmtostr(tm, tmstr)

[编译连接]

    Fortran固定格式代码，任何fortran编译器，无需任何外部连接库。

[算法公式]参见ETideLoad4.5参考说明书

    8.1地面及地球外部大地测量固体潮效应

    8.2.3规格化缔合勒让德函数及对θ导数

    8.3.3勒让德函数及对ψ一、二阶导数

附件rar压缩包：visual studio_intel fortran 集成环境测试项目、DOS可执行测试程序、全部地球物理模型和测试输入输出数据。