subroutine LGIntglr(BLH,ewh,hd,nlat,nlon,GF,L0,tdn,CRS,dr) !由地面等效水高变化计算负荷形变的间接影响 !ehw没有包含比例因子1/R,而Stokes系数包含了比例因子1/R !tdn(9)-大地水准面形变(mm) !扰动重力变化(uGal) !垂线偏差北南(子午圈)分量变化(ms) !垂线偏差东西(卯酉圈)分量变化(ms) !地面站点重力变化(uGal) !地面站点东方向形变(mm) !地面站点北方向形变(mm) !地面站点大地高变化(mm) !地面站点正(常)高变化(mm) !2016年7月2日,章传银 implicit none integer::nlat,nlon,nn,kk,i,j,maxn,ni,nj,i0,j0 real*8::ewh(nlat,nlon),hd(6),BLH(3),GF(9000,6),tdn(9),CRS(5),mdr,pi,RAD real*8::ae,gg,gm,em,gr,vfn(6),dl,ds,rw,dv,rr,psi,dr,r1,L1,L0,r0,rv real*8::rln(3),XYZ(3),BLH1(3),BLH0(3),rln1(3),XYZ0(3),XYZ1(3),NFD(5) real*8 rlon,rlat,rlon1,rlat1,sin2f,cos2f,sinf,sina,cosa,tmp !--------------------------------------------------------------------------- gg=6.67428d-11;gm=CRS(1);ae=CRS(2);rw=1.025d3;em=gm/gg tdn=0.d0;BLH0=BLH;BLH0(3)=0.d0;BLH1(3)=0.d0 pi=datan(1.d0)*4.d0;RAD=pi/180.d0;rv=hd(5)*RAD/4.d0 call BLH_RLAT(CRS,BLH,rln);rr=rln(1) rlat=rln(2)*RAD;rlon=rln(3)*RAD call BLH_XYZ(CRS,BLH,XYZ);call BLH_XYZ(CRS,BLH0,XYZ0) call GNormalfd(BLH,NFD);gr=NFD(2) r0=dsqrt(XYZ0(1)**2+XYZ0(2)**2+XYZ0(3)**2) mdr=r0*hd(5)*RAD*dcos(rlat)!奇异点判断 ni=nint(dr/hd(6)+1.d0) nj=nint(dr/hd(5)/dcos(rln(2)*RAD)+1.d0) !积分半径dr对应的地面格网数 i0=nint((BLH(1)-hd(3))/hd(6)+0.5d0) j0=nint((BLH(2)-hd(1))/hd(5)+0.5d0)!计算点所在的地面格网i0,j0 do i=i0-ni,i0+ni if(i<1.or.i>nlat)goto 9100 BLH1(1)=hd(3)+(real(i)-0.5d0)*hd(6) do j=j0-nj,j0+nj if(j<1.or.j>nlon)goto 9101 if(dabs(ewh(i,j))<1.d-8)goto 9101 BLH1(2)=hd(1)+(real(j)-0.5d0)*hd(5) call BLH_RLAT(CRS,BLH1,rln1);r1=rln1(1) rlat1=rln1(2)*RAD;rlon1=rln1(3)*RAD call BLH_XYZ(CRS,BLH1,XYZ1) dl=dsqrt((XYZ1(1)-XYZ(1))**2+(XYZ1(2)-XYZ(2))**2+(XYZ1(3)-XYZ(3))**2) L1=dsqrt((XYZ1(1)-XYZ0(1))**2+(XYZ1(2)-XYZ0(2))**2+(XYZ1(3)-XYZ0(3))**2) if(L1>dr*RAD*r0)goto 9101 sin2f=L1/r0/2.d0;cos2f=dsqrt(1.d0-sin2f**2);sinf=2.d0*sin2f*cos2f ds=r1**2*hd(5)*hd(6)*RAD**2*dcos(rlat1);psi=dasin(sin2f)*2.d0 call InterpGrF(GF,psi,vfn) if(psiL0)then cosa=(dcos(rlat)*dsin(rlat1)-dsin(rlat)*dcos(rlat1)*dcos(rlon1-rlon))/sinf sina=dcos(rlat1)*dsin(rlon1-rlon)/sinf tdn(3)=tdn(3)+tmp*vfn(3)*cosa;tdn(4)=tdn(4)+tmp*vfn(3)*sina tdn(6)=tdn(6)+tmp*vfn(5)*sina;tdn(7)=tdn(7)+tmp*vfn(5)*cosa endif else dv=rw*ewh(i,j)*ds;tdn(1)=tdn(1)+dv*vfn(1)/dl;tdn(2)=tdn(2)+dv*vfn(2)/dl tdn(5)=tdn(5)+dv*vfn(4)/dl;tdn(8)=tdn(8)+dv*vfn(6)/dl if(dl>L0)then cosa=(dcos(rlat)*dsin(rlat1)-dsin(rlat)*dcos(rlat1)*dcos(rlon1-rlon))/sinf sina=dcos(rlat1)*dsin(rlon1-rlon)/sinf tdn(3)=tdn(3)+dv*vfn(3)*cosa/dl;tdn(4)=tdn(4)+dv*vfn(3)*sina/dl tdn(6)=tdn(6)+dv*vfn(5)*sina/dl;tdn(7)=tdn(7)+dv*vfn(5)*cosa/dl endif endif 9101 continue enddo 9100 continue enddo tdn(1)=tdn(1)*1.d3 tdn(2)=tdn(2)*1.0e8 tdn(3)=tdn(3)/RAD*36.d5 tdn(4)=tdn(4)/RAD*36.d5 tdn(6)=tdn(6)*1.d3 tdn(7)=tdn(7)*1.d3 tdn(8)=tdn(8)*1.d3 tdn(5)=tdn(2)-tdn(8)*0.3086d0 !!!! tdn(9)=tdn(8)-tdn(1) return end