地球重力场位系数模型全空间全要素计算Fortran代码

[计算目标]

输入地球重力位系数模型，计算地球空间任意点处的高程异常（m）、空间异常（mGal）、扰动重力（mGal）、垂线偏差向量（s，ʺ/秒，南向、西向）、扰动重力梯度（径向，E）、水平重力梯度向量（E，北向、西向）和Laplace算子（E）的模型值。

当设置相等的重力位系数模型最小、最大计算阶数都为n时，程序计算第n阶位系数对扰动重力场元的贡献。可用于分析评价该重力位系数模型的谱域空域性质。

[测试入口程序]

    GGMallelementgrfd.f90

    输入计算点空间位置文件记录格式：点号/点名 经度（度小数） 纬度（度小数） 大地高（m）......

   输出文件reslt.txt记录：在输入文件记录的基础上，增加9列重力场元模型值，包括高程异常（m）、空间异常（mGal）、扰动重力（mGal）、垂线偏差向量（s，ʺ/秒，南向、西向）、扰动重力梯度（径向，E）、水平重力梯度向量（E，北向、西向）和Laplace算子（E）的模型值。

[主要调用模块]

（1）扰动重力场元模型值球谐综合计算模块

    RntGravFdpnm(nmin,maxn,rln,cnm,snm,gvm,GRS,pnm,dpt1,dpt2,gr)

    输入nmin,maxn－最小、最大计算阶数

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

    输入cnm,snm－模型重力场位系数

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

    输入pnm,dpt1,dpt2,gr－计算点处的连带勒让德函数及其一、二阶导数和正常重力值

    返回gvm(8)－高程异常(m)、空间异常(mGal)、扰动重力(mGal)、南向、西向垂线偏差(s)、重力梯度(E)、北向、东向水平重力梯度(E)

（2）正常重力位系数计算模块

    normdjn(GRS,djn)；cstonorm (maxn,djn,GRS,HC,HS)

（3）正常重力场元计算模块

    GNormalfd(BLH,NFD,GRS)

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

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

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

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

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

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

    LegPn_dt2(pn,dp1,dp2,n,t)

    计算勒让德函数Pn(t)及其对ψ一、二阶导数t=cosψ。

（6）正常椭球基本参数计算模块

    ELLIPSOIDPARA(GRS)

    GRS(1)-地心引力常数GM；GRS(2)-椭球长半轴；GRS(3)-地球动力学形状因子J2；GRS(4)-自转平均角速度omega；GRS(5)-地球椭球扁率1/f；GRS(6)-椭球面正常重力位U₀=W_G。

    输入参数：GRS(1)，GRS(2)，GRS(4)和GRS(3)、GRS(5)、GRS(6)其中之一，设置没有赋值的剩余GRS(3)、GRS(5)或GRS(6)为任意大于零的数。

    输出参数：GRS(1:6)

（7）椭球大地坐标变换为球坐标模块

    BLH_RLAT(GRS,BLH,RLAT)

[编译连接]

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

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

    7.2重力场位系数模型场元计算公式

    7.3缔合勒让德函数及其导数算法

提示：在地球两极处，垂线偏差西向、水平重力梯度西向无定义，计算结果无意义。

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