subroutine permanentTide(BLH,tdn,CRS) !计算永久潮汐直接影响和间接影响之和 !tdn(9)-大地水准面形变(mm) !扰动重力变化(uGal) !地倾斜北南(子午圈)分量变化(ms) !地倾斜东西(卯酉圈)分量变化(ms) !地面站点重力变化(uGal) !地面站点东方向形变(mm) !地面站点北方向形变(mm) !地面站点大地高变化(mm) !地面站点正(常)高变化(mm) !2015年7月21日,章传银 implicit none integer::nn,n,m,kk,i,j real*8::BLH(3),rln(3),gm,tdn(9),tn(9),gr,CRS(5) real*8::cnm(40),snm(40),cosml,sinml,XYZ(3),ZCFD(4) real*8::pnm(40),pnm1(40),din(4),knm(7),hnm(7),lnm(7) real*8::p(40),dp(40) real*8::rr,rlat,rlon,ae,pi,RAD !--------------------------------------------------------------------------- pi=datan(1.d0)*4.d0;RAD=pi/180.d0 gm=CRS(1);ae=CRS(2); call BLH_RLAT(CRS,BLH,rln); call RLAT_XYZ(rln,XYZ) call normal(CRS,XYZ,ZCFD);gr=ZCFD(2) knm=0.d0; hnm=0.d0; lnm=0.d0 knm(1)=0.29525;knm(2)=0.29470;knm(3)=0.29801 knm(4)=0.093;knm(5)=0.093;knm(6)=0.093;knm(7)=0.094 hnm(1)=0.6078d0;hnm(2)=0.6078d0;hnm(3)=0.6078d0 hnm(4)=0.292d0;hnm(5)=0.292d0;hnm(6)=0.292d0;hnm(7)=0.292d0 lnm(1)=0.0847d0;lnm(2)=0.0847d0;lnm(3)=0.0847d0 lnm(4)=0.015d0;lnm(5)=0.015d0;lnm(6)=0.015d0;lnm(7)=0.015d0 tdn=0.d0; rr=rln(1);rlat=rln(2);rlon=rln(3) ! 计算全部Legendre函数pnm,数组pnm的下标比cnm和snm大3 call PlmBar_d(pnm, pnm1, 3, rlat) tn=0.d0 do m=0,2 j=m+1 cosml=dcos(real(m)*rlon*RAD);sinml=dsin(real(m)*rlon*RAD) tn(1)=tn(1)+(cnm(j)*cosml+snm(j)*sinml)*pnm(j+3)*(1.d0+knm(j)) tn(2)=tn(2)-(cnm(j)*cosml+snm(j)*sinml)*pnm(j+3)*(real(n)-real(n+1)*knm(j)) tn(3)=tn(3)-(cnm(j)*cosml+snm(j)*sinml)*pnm1(j+3)*(1.d0+knm(j)-hnm(j)) tn(4)=tn(4)-real(m)*(cnm(j)*sinml-snm(j)*cosml)*pnm(j+3)*(1.d0+knm(j)-hnm(j)) tn(5)=tn(5)-(cnm(j)*cosml+snm(j)*sinml)*pnm(j+3)*(real(n)- > real(n+1)*knm(j)+2.d0*hnm(j)) tn(6)=tn(6)+real(m)*(cnm(j)*sinml-snm(j)*cosml)*pnm(j+3)*lnm(j) tn(7)=tn(7)+(cnm(j)*cosml+snm(j)*sinml)*pnm1(j+3)*lnm(j) tn(8)=tn(8)+(cnm(j)*cosml+snm(j)*sinml)*pnm(j+3)*hnm(j) enddo tdn(1)=tn(1)*gm/rr/gr*1.d3 tdn(2)=tn(2)*gm/rr**2*1.0e8 tdn(3)=tn(3)*gm/rr**2/gr/RAD*36.d5 tdn(4)=tn(4)*gm/rr**2/gr/dcos(rlat*RAD)/RAD*36.d5 tdn(5)=tn(5)*gm/rr**2*1.0e8 tdn(6)=tn(6)*gm/rr/gr/dcos(rlat*RAD)*1.d3 tdn(7)=tn(7)*gm/rr/gr*1.d3 tdn(8)=tn(8)*gm/rr/gr*1.d3 tdn(9)=tdn(8)-tdn(1) return end ! !****************************************************************************** ! subroutine PlmBar_d(p,dp,lmax,rlat) implicit none integer::lmax,m,n,kk real*8 ::p((lmax+2)*(lmax+3)/2+5),dp((lmax+2)*(lmax+3)/2+5) real*8 ::rleg(lmax+5),dleg(lmax+5),rlat !--------------------------------------------------------------------------- ! 计算全部Legendre函数pnm,数组pnm的下标比cnm和snm大3 do n=1,lmax+1 call Legendre(lmax,n-1,rlat,rleg,dleg) do m=n,lmax+1 kk=m*(m-1)/2+n p(kk)=rleg(m) dp(kk)=dleg(m) enddo enddo end