subroutine LoadrDirectFunc(BLH,ewh,dtm,hd,nlat,nlon,tdn,CRS,dr) !------------------------------------------------------------- implicit none integer::nlat,nlon,i,j,ni,nj,i0,j0 real*8::dtm(nlat,nlon),ewh(nlat,nlon),BLH(3),hd(6),tdn(6),CRS(5) real*8::pi,gg,gr,RAD,L1,ds,rw,L0,mdr,dv,dwh,dr,rst(4) real*8::rln(3),XYZ(3),BLH1(3),rln1(3),XYZ1(3),sin2f,tt,ZCFD(3) real*8::BLH0(3),rln0(3),XYZ0(3),BLH2(3),rln2(3),XYZ2(3),rr,r0,r1 real*8 rlon,rlat,rlon1,rlat1,dlon,dlat,CGrdPntD2 !----------------------------------------------------------------- gg=6.67428d-11;rw=1.d3;pi=datan(1.d0)*4.d0;RAD=pi/180.d0;tdn=0.d0;BLH0=BLH BLH0(3)=CGrdPntD2(BLH(2),BLH(1),dtm,nlat,nlon,hd)!计算点正下方地面点位置 call BLH_RLAT(CRS,BLH,rln);call BLH_XYZ(CRS,BLH,XYZ) rr=rln(1);call normal(CRS,XYZ,ZCFD);gr=ZCFD(2) rlat=rln(2)*RAD;rlon=rln(3)*RAD call BLH_XYZ(CRS,BLH0,XYZ0) r0=dsqrt(XYZ0(1)**2+XYZ0(2)**2+XYZ0(3)**2) mdr=r0*hd(5)*RAD*dcos(rlat)/2.d0 !奇异点判断 ni=nint(dr/hd(6)+1.d0) nj=nint(dr/hd(5)/dcos(rlat)+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) BLH1(3)=dtm(i,j);call BLH_XYZ(CRS,BLH1,XYZ1) L0=dsqrt((XYZ1(1)-XYZ0(1))**2+(XYZ1(2)-XYZ0(2))**2+(XYZ1(3)-XYZ0(3))**2) if(L0>dr*RAD*r0)goto 9101 L1=dsqrt((XYZ1(1)-XYZ(1))**2+(XYZ1(2)-XYZ(2))**2+(XYZ1(3)-XYZ(3))**2) call BLH_RLAT(CRS,BLH1,rln1);r1=rln1(1) sin2f=L0/r0/2.d0;tt=1.d0-2.d0*sin2f**2 rlat1=rln1(2)*RAD;rlon1=rln1(3)*RAD ds=hd(5)*hd(6)*RAD**2*dcos(rlat1)*r1**2 dlat=rlat1-rlat;dlon=rlon1-rlon dwh=gg*rw*ewh(i,j);dv=dwh*ds if(L12.5d0*mdr)tdn(3)=tdn(3)-dv*r1/L1**3*(dcos(rlat)*dsin(rlat1)-dsin(rlat)*dcos(rlat1)*dcos(dlon)) if(L1>2.5d0*mdr)tdn(4)=tdn(4)-dv*r1/L1**3*dcos(rlat)*dcos(rlat1)*dsin(dlon) endif 9101 continue enddo 9100 continue enddo tdn(1)=tdn(1)/gr*1.d3 tdn(2)=tdn(2)*1.0e8 tdn(3)=tdn(3)/gr*36.d5/RAD tdn(4)=tdn(4)/gr/dcos(rlat)*36.d5/RAD tdn(5)=tdn(2);tdn(6)=-tdn(1) return end ! !****************************************************************!!!!! ! subroutine LoadDrctSgn(BLH,ewh,dtm,hd,nlat,nlon,rst,CRS,i0,j0,m) !------------------------------------------------------------- implicit none integer::nlat,nlon,i,j,ni,nj,i0,j0,m real*8::dtm(nlat,nlon),ewh(nlat,nlon),BLH(3),hd(6),rst(4),CRS(5) real*8::pi,gg,gr,RAD,L1,ds,rw,L0,mdr,dv,dwh,rv,lat,lon,dw real*8::rln(3),XYZ(3),BLH1(3),rln1(3),XYZ1(3),sin2f,tt,ZCFD(3) real*8::BLH0(3),rln0(3),XYZ0(3),BLH2(3),rln2(3),XYZ2(3),rr,r0,r1 real*8 rlon,rlat,rlon1,rlat1,dlon,dlat,CGrdPntD2 !----------------------------------------------------------------- gg=6.67428d-11;rw=1.d3;pi=datan(1.d0)*4.d0;RAD=pi/180.d0;BLH0=BLH BLH0(3)=CGrdPntD2(BLH(2),BLH(1),dtm,nlat,nlon,hd)!计算点正下方地面点位置 rv=hd(5)/dble(m);pi=datan(1.d0)*4.d0;RAD=pi/180.d0 lat=hd(3)+real(i0-1)*hd(6);lon=hd(1)+real(j0-1)*hd(5)!格网左下角经纬度 call BLH_XYZ(CRS,BLH,XYZ);call BLH_RLAT(CRS,BLH,rln) rr=rln(1);rlat=rln(2)*RAD;rlon=rln(3)*RAD call BLH_XYZ(CRS,BLH0,XYZ0) r0=dsqrt(XYZ0(1)**2+XYZ0(2)**2+XYZ0(3)**2) rst=0.d0;mdr=r0*rv*RAD*dcos(rlat)/4.d0 !奇异点判断 do i=1,m BLH1(1)=lat+(real(i)-0.5d0)*rv do j=1,m BLH1(2)=lon+(real(j)-0.5d0)*rv call BLH_XYZ(CRS,BLH1,XYZ1) BLH1(3)=CGrdPntD2(BLH1(2),BLH1(1),dtm,nlat,nlon,hd) dw=CGrdPntD2(BLH1(2),BLH1(1),ewh,nlat,nlon,hd) call BLH_RLAT(CRS,BLH1,rln1);r1=rln1(1);rlat1=rln1(2)*RAD;rlon1=rln1(3)*RAD L1=dsqrt((XYZ1(1)-XYZ(1))**2+(XYZ1(2)-XYZ(2))**2+(XYZ1(3)-XYZ(3))**2) L0=dsqrt((XYZ1(1)-XYZ0(1))**2+(XYZ1(2)-XYZ0(2))**2+(XYZ1(3)-XYZ0(3))**2) sin2f=L0/r0/2.d0;tt=1.d0-2.d0*sin2f**2 ds=rv**2*RAD**2*dcos(rlat1)*r1**2 dlat=rlat1-rlat;dlon=rlon1-rlon dwh=gg*rw*dw;dv=dwh*ds if(L1