广义Stokes/Hotine数值积分外部高程异常计算Fortran代码

[计算目标]

    由等位边界面大地高格网（m）及其面上的残差空间异常/扰动重力格网（mGal），按广义Stokes/Hotine数值积分方法，计算近地空间大地水准面及其外部残差高程异常（m）。

    等位边界面不必是大地水准面，可以是某一等高面，甚至可以由地面外部等位面（等高面）上的残差空间异常/扰动重力，计算地面残差高程异常（教科书中一般不会告诉你哦）。

    大地水准面上的高程异常，即为大地水准面差距或大地水准面（大地）高。Stokes边值问题要求，边界面必须是等位面，即空间异常/扰动重力必须位于重力等位面上。

    为实现有限半径积分，通常需采用参考重力场移去恢复法，先移去等位边界面上的模型空间异常/扰动重力，再积分得到计算点处的残差高程异常，最后恢复计算点处的模型高程异常。

    等位面可采用参考重力场模型（不大于360阶）构造，在高度不大于10千米的近地空间，可用等正（常）高面大地高格网表示。

[测试入口程序]

    StokesHotinenumintegral.f90

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

    输入参数mode－=0广义stokes积分；=1广义Hotine积分。

    输入等位边界面大地高格网文件－表示等位边界面的位置，用于计算积分面元与计算点的积分距离。

    输入与等位边界面大地高格网规格相同的空间异常（用于Stokes积分）或扰动重力（用于Hotine积分）格网文件。

    输出文件reslt.txt记录：在输入文件记录的基础上，增加1列残差高程异常积分值（m）。

[主要调用模块]

（1）广义Stokes积分算法模块

    Real*8 StokesBLH(BLH,gra,sfh,nlat,nlon,hd,dr,GRS)

    输入BLH(3)－空间计算点的经纬度（度小数）和大地高（m）。

    输入sfh(nlat,nlon)－等位边界面大地高格网，用于精确计算积分面元与计算点的积分距离。

    输入gra(nlat,nlon)－等位边界面上空间异常格网（mGal）。

    输入dr, hd(6)－积分半径（m）和格网规格参数（最小最大经度，最小最大纬度，经度间隔，纬度间隔）。

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

    返回－计算点处残差高程异常（m）。

（2）广义Hotine积分算法模块

    Real*8 HotineBLH(BLH,rga,sfh,nlat,nlon,hd,dr,GRS)

    输入rga(nlat,nlon)－等位边界面上扰动重力格网（mGal）。

    返回－计算点处残差高程异常（m）。

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

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

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

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

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

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

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

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

    RLAT_BLH(GRS, RLAT, BLH)

（6）格网内插点值算法包

    CGrdPntD(lon,lat,dt,row,col,hd)；CGrdPntD2(lon,lat,dt,row,col,hd)

    CShepard(lon,lat,dt,row,col,hd)；Gauss2D(lon,lat,dt,row,col,hd)

（7）其他辅助模块

    PickRecord(str0, kln, rec, nn)

[编译连接]

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

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

    1.4.1大地测量数据文件格式约定

    7.9.1广义Stokes与Hotine积分公式

    7.1（4）低阶勒让德函数及其一、二阶导数算法

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