subroutine FrNetAdjust(infl1,infl2,infl3,outfl,knd,row,row1,nm,pk,nk) !基于监测量的控制网拟稳平差 !infl1-GNSS基线分量、水准测段高差或重力段差变化量文件 !infl2、infl3-监测网拟稳基准点经纬度坐标文件,监测网站点经纬度坐标文件(无头文件) !nm-站名所占的字符数 !knd-0拟稳基准零约束,1拟稳站点平均约束 USE f95_precision, ONLY: WP => DP USE lapack95 implicit none character*800::infl1,infl2,infl3,outfl character*8000::line character(len=20)::znm(8000) integer::nk,knd,row,row1,pk,nd,nm,i,j,k,nlon,nlat,len,sn,nstr,nn,kk,info integer::ndn(8000),nb(800000,2),km(800000) !ndn拟稳基准点的未知数序号,nb变化量基线端点序号 real*8::rec(800),obs(800000,6),zlb(8000,2),dm,qp,tmp,avln,dmav real*8::dtmn(800,4),qsm,qq(800)!dtmn拟稳基准点经纬度、权值 integer::status=0 !zlb全部站点经纬度 real*8,allocatable::BB(:),BPB(:,:),BPL(:),xx(:) !--------------------------------------------------------------------- open(unit=8,file=infl1,status="old",iostat=status) if(status/=0)goto 902 if(nk<1)goto 201 do i=1,nk read(8,'(a)') line !读取头文件 enddo 201 continue nd=0;dmav=0.d0 !拟稳基准点的个数 do while(.not.eof(8)) read(8,'(a)') line call PickRecord(line,len,rec,sn) if(sn<4.or.sn9900.d0)goto 906 nd=nd+1;dtmn(nd,1)=rec(2);dtmn(nd,2)=rec(3) !拟稳基准点经纬度 dtmn(nd,3)=rec(row1);dtmn(nd,4)=rec(pk)!拟稳基准点权值、参考值 906 continue enddo 905 close(8) dmav=dmav/tmp open(unit=8,file=infl2,status="old",iostat=status) read(8,'(a)') line call PickRecord(line,len,rec,sn) nstr=nint(rec(1)) !确定未知数个数nn,站点序号与站名、基线的对应关系 nn=0;kk=0;nb=0;avln=0.d0 !平均基线长度 do while(.not.eof(8)) read(8,'(a)') line call PickRecord(line,len,rec,sn) if(sn9900.d0)goto 401 kk=kk+1;obs(kk,1:5)=rec(2:6);obs(kk,6)=rec(row) avln=avln+obs(kk,5) if(kk==1)then nb(kk,1)=1;nb(kk,2)=2;nn=2 znm(1)=line(1:nm);znm(2)=line(nstr-nm+1:nstr) zlb(1,1:2)=rec(2:3);zlb(2,1:2)=rec(4:5) endif if(kk>1)then do i=1,nn if(dsqrt((zlb(i,1)-rec(2))**2+(zlb(i,2)-rec(3))**2)<1.d-6)nb(kk,1)=i if(dsqrt((zlb(i,1)-rec(4))**2+(zlb(i,2)-rec(5))**2)<1.d-6)nb(kk,2)=i enddo if(nb(kk,1)<1)then nn=nn+1;nb(kk,1)=nn;znm(nn)=line(1:nm);zlb(nn,1:2)=rec(2:3) endif if(nb(kk,2)<1)then nn=nn+1;nb(kk,2)=nn;znm(nn)=line(nstr-nm+1:nstr);zlb(nn,1:2)=rec(4:5) endif endif 401 continue enddo close(8) !确定nd个拟稳基准点在nn个未知数中的序号 avln=avln/dble(kk);qsm=0.d0;dmav=0.d0 !拟稳权值之和 do i=1,nd do j=1,nn if(dsqrt((dtmn(i,1)-zlb(j,1))**2+(dtmn(i,2)-zlb(j,2))**2)<1.d-3)then ndn(i)=j;qsm=qsm+dtmn(i,3);dmav=dmav+dtmn(i,3)*dtmn(i,4);goto 302 endif enddo 302 continue enddo dmav=dmav/qsm k=nn+nd+5 allocate(BB(nn),BPB(k,nn),BPL(k),xx(nn)) BPB=0.d0;BPL=0.d0;xx=0.d0 do k=1,kk BB=0.d0;BB(nb(k,1))=-1.d0;BB(nb(k,2))=1.d0 qp=avln/obs(k,5) do i=1,nn BPL(i)=BPL(i)+BB(i)*obs(k,6)*qp do j=1,i BPB(i,j)=BPB(i,j)+BB(i)*BB(j)*qp enddo enddo enddo do i=1,nn do j=1,i-1 BPB(j,i)=BPB(i,j) enddo enddo km=0 !拟稳基准条件构造 do i=1,nd if(ndn(i)>0)then BPB(nn+1,ndn(i))=1.d0;km(ndn(i))=1 endif enddo if(knd>0)BPL(nn+1)=dmav*3.d0 call gels(BPB,BPL,'N',info) xx(1:nn)=BPL(1:nn) open(unit=10,file=outfl,status="replace") open(unit=8,file=infl3,status="old",iostat=status) do while(.not.eof(8)) read(8,'(a)') line call PickRecord(line,len,rec,sn) if(sn<3)goto 202 kk=0 do i=1,nn if(dsqrt((rec(2)-zlb(i,1))**2+(rec(3)-zlb(i,2))**2)<1.d-3)then write(10,101)trim(line),xx(i),km(i);kk=kk+1;goto 202 endif enddo if(kk<1)write(10,101)trim(line),9999.0000,0 202 continue enddo close(8) close(10) deallocate(BB,BPB,BPL,xx) 902 continue 101 format(a,F10.4,i4) end