!************************************************************************************ !调用准备代码:省略变量声明 !读取地面等效水高球谐系数,负荷勒夫数 !给定时间BLH(3) rw=1.025d3; re=5.517d3 !打开地面等效水高球谐系数文件 allocate(cnm((maxn+2)**2), stat=astat(1)) allocate(snm((maxn+2)**2), stat=astat(2)) if (sum(astat(1:2)) /= 0) goto 908 cnm=0.d0;snm=0.d0 open(unit=8,file=clmfl,status="old",iostat=status) if(status/=0)goto 904 do while(.not.eof(8)) read(8,*,end=901)n,m,cnm(n*(n+1)/2+m),snm(n*(n+1)/2+m) enddo 901 close(8) if(n20000)maxn=20000 flv(1,1:3)=0.d0 !一阶项不参与计算,无需质心校正 call BLH_RLAT(CRS,BLH,rln) !********************************************************* !***************************************************************** subroutine GEqLoadTideCalc(rln,maxn,cnm,snm,flv,tdn,CRS,mode) !地面等效水高球谐综合法负荷形变场计算-计算空间点rln(3)的负荷形变 !mode-0直接影响,1间接影响,3负荷总影响=直接影响+间接影响 !当负荷为大气压时,对地面重力或地面扰动重力:负荷总影响=间接影响-直接影响 !tdn(k)-1高程异常/大地水准面形变(mm) !2扰动重力变化(uGal) !3地倾斜北南(子午圈)分量变化(ms) !4地倾斜东西(卯酉圈)分量变化(ms) !5地面站点重力变化(uGal) !6地面站点东方向形变(mm) !7地面站点北方向形变(mm) !8地面站点大地高变化(mm) !9地面站点正(常)高变化(mm) !2015年7月21日,章传银 implicit none integer::maxn,mode,nn,n,m,kk,i,j real*8::cnm((maxn+2)**2),snm((maxn+2)**2) real*8::rln(3),gm,tdn(9),tn(9),gr,CRS(5),cosml,sinml real*8::BLH(3),NFD(5),t,u real*8::rr,rlat,rlon,ae,pi,RAD,flv(40000,3),dk,dh,dl real*8,allocatable::pnm(:),dpt1(:),dpt2(:) !--------------------------------------------------------------------------- tdn=0.d0;gm=CRS(1);ae=CRS(2) pi=datan(1.d0)*4.d0;RAD=pi/180.d0 allocate(pnm((maxn+2)**2),dpt1((maxn+2)**2),dpt2((maxn+2)**2)) rr=rln(1);rlat=rln(2);rlon=rln(3) t=dsin(rlat*RAD);u=dcos(rlat*RAD) call RLAT_BLH(CRS,rln,BLH);call GNormalfd(BLH,NFD);gr=NFD(2) call BelPnmdt(pnm,dpt1,dpt2,maxn,t) do n=2,maxn !一阶项不参与计算 dh=flv(n,1);dl=flv(n,2);dk=flv(n,3);tn=0.d0 do m=0,n j=n*(n+1)/2+m cosml=dcos(rlon*RAD*dble(m));sinml=dsin(rlon*RAD*dble(m)) tn(1)=tn(1)+(cnm(j)*cosml+snm(j)*sinml)*pnm(j+1) tn(2)=tn(2)+(cnm(j)*cosml+snm(j)*sinml)*pnm(j+1) tn(3)=tn(3)+(cnm(j)*cosml+snm(j)*sinml)*dpt1(j+1) tn(4)=tn(4)+real(m)*(cnm(j)*sinml-snm(j)*cosml)*pnm(j+1) tn(5)=tn(5)+(cnm(j)*cosml+snm(j)*sinml)*pnm(j+1) tn(6)=tn(6)+real(m)*(cnm(j)*sinml-snm(j)*cosml)*pnm(j+1) tn(7)=tn(7)+(cnm(j)*cosml+snm(j)*sinml)*dpt1(j+1) tn(8)=tn(8)+(cnm(j)*cosml+snm(j)*sinml)*pnm(j+1) enddo tn=tn*dexp(real(n)*dlog(ae/rr)) if(mode==0)then tdn(1)=tdn(1)+tn(1) tdn(2)=tdn(2)+tn(2)*real(n+1) tdn(3)=tdn(3)+tn(3) tdn(4)=tdn(4)+tn(4) tdn(5)=tdn(5)+tn(5)*real(n+1) tdn(6)=0.d0 tdn(7)=0.d0 tdn(8)=0.d0 endif if(mode==1)then tdn(1)=tdn(1)+tn(1)*dk tdn(2)=tdn(2)+tn(2)*real(n+1)*(-real(n+1)*dk)/real(n) tdn(3)=tdn(3)+tn(3)*(dk-dh) tdn(4)=tdn(4)+tn(4)*(dk-dh) tdn(5)=tdn(5)+tn(5)*real(n+1)*(-real(n+1)*dk+2.d0*dh)/real(n) tdn(6)=tdn(6)+tn(6)*dl tdn(7)=tdn(7)+tn(7)*dl tdn(8)=tdn(8)+tn(8)*dh endif if(mode==2)then tdn(1)=tdn(1)+tn(1)*(1.d0+dk) tdn(2)=tdn(2)+tn(2)*real(n+1)*(1.d0-real(n+1)*dk/real(n)) tdn(3)=tdn(3)+tn(3)*(1.d0+dk-dh) tdn(4)=tdn(4)+tn(4)*(1.d0+dk-dh) tdn(5)=tdn(5)+tn(5)*real(n+1)*(1.d0+(2.d0*dh-real(n+1)*dk)/real(n)) tdn(6)=tdn(6)+tn(6)*dl tdn(7)=tdn(7)+tn(7)*dl tdn(8)=tdn(8)+tn(8)*dh endif enddo tdn(1)=tdn(1)*gm/rr/gr*1.d3 tdn(2)=tdn(2)*gm/rr**2*1.0e8 tdn(3)=tdn(3)*gm/rr**2/gr/RAD*36.d5 tdn(4)=tdn(4)*gm/rr**2/gr/u/RAD*36.d5 !tdn(5)=tdn(5)*gm/rr**2*1.0e8 tdn(6)=-tdn(6)*gm/rr/gr/u*1.d3 tdn(7)=-tdn(7)*gm/rr/gr*1.d3 tdn(8)=tdn(8)*gm/rr/gr*1.d3 tdn(5)=tdn(2)-tdn(8)*0.3086d0 !!!! tdn(9)=tdn(8)-tdn(1) deallocate(pnm,dpt1,dpt2) return end