!************************************************************************************ !调用准备代码:省略变量声明 !给定时间XYZ(3),tm !CRS(5)椭球参数请自行添加 !读大气潮负荷球谐系数模型 open(unit=8,file=airtdcsfl,status="old",iostat=status) if(status/=0)goto 902 read(8,'(a)') line read(8,'(a)') line !FES2004(11) read(8,'(a)') line !FES2004(11) i=1 fes=0.d0 do while(.not.eof(8)) read(8,*,end=903)fes(i,1),fh,(fes(i,j),j=2,3),tmp,tmp,tmp,tmp,(fes(i,j),j=4,7) fes(i,1)=fes(i,1)*1.d3 !FES2004(11) i=i+1 enddo 903 close(8) nn=i-1;maxn=fes(nn,2) !读负荷勒夫数 open(unit=8,file=loadlovfl,status="old",iostat=status) if(status/=0)goto 902 n=0 do while(.not.eof(8)) n=n+1 read(8,*,end=904)(flv(n,j),j=1,3) enddo 904 close(8) flv(1,1:3)=0.d0 !一阶项不参与计算,无需质心校正 flv(:,1:2)=0.d0 !地球外部h=0,l=0 !********************************************************* !***************************************************************** subroutine satAirtdLoad(tm,XYZ,tdn,dxyz,fes,flv,nn,maxn) implicit none !输出:tdn(4)-重力位摄动(0.001m2/s2)和摄动力在球坐标系中三分量(μGal) !dxyz(3)-摄动力在空间直角坐标系中三分量(μGal) integer::knd,row,i,j,k,nn,IY,IM,ID,ih,imn,maxn real*8::tm,pi,RAD,djm,sec,,td real*8::tdn(9),cnm(900000),snm(900000),XYZ(3),BLH(3),rln(3) real*8::L0,de,rr,flv(40000,3),fes(900000,7) real*8::gr,ZCFD(4),su,dxyz(3) !--------------------------------------------------------------------------- pi=datan(1.d0)*4.d0;RAD=pi/180.d0 su=1.d8/36.d5*RAD call XYZ_BLH(CRS,XYZ,BLH) call BLH_RLAT(CRS,BLH,rln) call tmcnt(tm,IY,IM,ID,ih,imn,sec,tmzn)!北京时间tmzn=8 call CAL2JD (IY,IM,ID,djm,j) call CAL2JD (IY,IM,ID,djm,j) djm=djm+real(ih)/24.d0+real(imn)/1440.d0+real(sec)/864.d2 !GPS_MJD td=djm+2400000.5d0+19.d0/864.d2!; td=td-TAI_UTC(td)/864.d2 !UTC_MJD 处理跳秒 tdn=0.d0; cnm=0.d0; snm=0.d0 call AirLoadFluCS(td,cnm,snm,maxn,fes,nn) !!cnm,snm直接影响位系数 cnm(1:2)=0.d0;snm(1:2)=0.d0 call LTideFlu(rln,maxn,cnm,snm,flv,tdn,CRS,2)!2-计算海潮负荷总影响 call normal(CRS,XYZ,ZCFD);gr=ZCFD(2) tdn(1)=tdn(1)*gr;tmp=tdn(4);tdn(4)=-tdn(2) tdn(2)=-tmp*gr*su;tdn(3)=-tdn(3)*gr*su call rnl_XYZd(rln,tdn(2:4),dxyz) end