      real*8 function AnVMdeflBLH(BLH,kai,nta,hgt,kxy,nlat,nlon,hd,dr,knd)
      !按严密球面积分公式计算垂线偏差逆远算积分
      !dr-积分半径m
!-------------------------------------------------------------
      implicit none
	integer::knd,i,j,nlat,nlon,i0,j0,ni,nj,astat(5)
	real*8::dr,kai(nlat,nlon),nta(nlat,nlon),hgt(nlat,nlon),kxy(nlat,nlon)
	real*8::hd(6),pi,RAD,ds,mdr,gr,tt,rr,r1,rst,tmp,qp
	real*8::BLH(3),XYZ(3),rln(3),ZCFD(3),BLH1(3),XYZ1(3),rln1(3)
	real*8 rlon,rlat,rlon1,rlat1,dlon,dlat
	real*8 sin2f,cos2f,sinf,sina,cosa
	real*8 CGrdPntD2,AnVMdeflSng,L1
	real*8::CRS(5)
      integer::kfn !奇异积分核函数精化参数2*kfn
!-----------------------------------------------------------------
      rst=0.d0;pi=datan(1.d0)*4.d0;RAD=pi/180.d0;kfn=2
      call BLH_XYZ(CRS,BLH,XYZ);call BLH_RLAT(CRS,BLH,rln);rr=rln(1)
	call normal(CRS,XYZ,ZCFD);gr=ZCFD(2)
      rlat=rln(2)*RAD;rlon=rln(3)*RAD
      mdr=rr*hd(5)*RAD*dcos(rln(2)*RAD)/2.d0 !奇异点判断
      ni=nint(dr/rr/RAD/hd(6)+1.d0) !积分半径dr对应的地面格网数
      nj=nint(dr/rr/RAD/hd(5)/dcos(rln(2)*RAD)+1.d0)
	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
	    BLH1(2)=hd(1)+(real(j)-0.5d0)*hd(5)
          BLH1(3)=hgt(i,j);call BLH_XYZ(CRS,BLH1,XYZ1)
          L1=dsqrt((XYZ1(1)-XYZ(1))**2+(XYZ1(2)-XYZ(2))**2+(XYZ1(3)-XYZ(3))**2)
          if(L1>dr)goto 9101
          if(L1<mdr)then!计算奇异积分
             rst=rst+AnVMdeflSng(BLH,kai,nta,hgt,kxy,nlat,nlon,hd,knd,i,j,kfn)
             goto 9101 
          endif
          call BLH_RLAT(CRS,BLH1,rln1);r1=rln1(1)
          ds=hd(5)*hd(6)*RAD**2*dcos(rln1(2)*RAD)
          tt=1.d0-2.d0*(L1/r1/2.d0)**2
          rlat1=rln1(2)*RAD;rlon1=rln1(3)*RAD
          sin2f=dsqrt((1.d0-tt)/2.d0);cos2f=dsqrt(1.d0-sin2f**2);sinf=2.d0*sin2f*cos2f
	    cosa=(dcos(rlat)*dsin(rlat1)-dsin(rlat)*dcos(rlat1)*dcos(rlon1-rlon))/sinf
	    sina=dcos(rlat1)*dsin(rlon1-rlon)/sinf
          qp=-(kai(i,j)*cosa+nta(i,j)*sina)*ds/pi/4.d0
          if(knd==0)rst=rst-qp*r1*cos2f/sin2f
          if(knd==1)then
	      tmp=3.d0/sinf-1.d0/sinf/sin2f-sin2f/cos2f-2.d0*cos2f/sin2f
            rst=rst+qp*gr*tmp
          endif
          if(knd==2)rst=rst+qp*gr*(3.d0/sinf-1.d0/sinf/sin2f-sin2f/cos2f)
9101      continue
	  enddo
9100    continue
	enddo
	AnVMdeflBLH=rst
9002	return
      end
!--------------------------------------------------------------------------------
      real*8 function AnVMdeflSng(BLH,kai,nta,hgt,kxy,nlat,nlon,hd,knd,i0,j0,m)
      !细化核函数，计算BLH点的垂线偏差逆远算的奇异积分
      !m-核函数细化为m*m
      !i0,j0-奇异点格网位置
!-------------------------------------------------------------
      implicit none
	integer::knd,m,i,j,nlat,nlon,i0,j0
	real*8::kai(nlat,nlon),nta(nlat,nlon),hgt(nlat,nlon),kxy(nlat,nlon)
	real*8::hd(6),gr,pi,RAD,ds,mdr,tt,rr,r1,rst,rv,tx,lon,lat
	real*8::BLH(3),XYZ(3),rln(3),ZCFD(3),BLH1(3),XYZ1(3),rln1(3)
	real*8 rlon,rlat,rlon1,rlat1,dlon,dlat,tmp,qp
	real*8 sin2f,cos2f,sinf,sina,cosa,CGrdPntD2,L1
	real*8::CRS(5)
!-----------------------------------------------------------------
      rst=0.d0;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)
	call normal(CRS,XYZ,ZCFD);gr=ZCFD(2)
      rlat=rln(2)*RAD;rlon=rln(3)*RAD
      mdr=rr*rv*RAD*dcos(rln(2)*RAD)/2.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
          BLH1(3)=CGrdPntD2(BLH1(2),BLH1(1),hgt,nlat,nlon,hd)
          call BLH_XYZ(CRS,BLH1,XYZ1)
          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)
          ds=rv**2*RAD**2*dcos(rln1(2)*RAD)
          if(L1<mdr)goto 9101
          call BLH_RLAT(CRS,BLH1,rln1);r1=rln1(1)
          tt=1.d0-2.d0*(L1/r1/2.d0)**2
          rlat1=rln1(2)*RAD;rlon1=rln1(3)*RAD
          sin2f=dsqrt((1.d0-tt)/2.d0);cos2f=dsqrt(1.d0-sin2f**2);sinf=2.d0*sin2f*cos2f
	    cosa=(dcos(rlat)*dsin(rlat1)-dsin(rlat)*dcos(rlat1)*dcos(rlon1-rlon))/sinf
	    sina=dcos(rlat1)*dsin(rlon1-rlon)/sinf
          qp=-(kai(i,j)*cosa+nta(i,j)*sina)*ds/pi/4.d0
          if(knd==0)rst=rst-qp*r1*cos2f/sin2f
          if(knd==1)then
	      tmp=3.d0/sinf-1.d0/sinf/sin2f-sin2f/cos2f-2.d0*cos2f/sin2f
            rst=rst+qp*gr*tmp
          endif
          if(knd==2)rst=rst+qp*gr*(3.d0/sinf-1.d0/sinf/sin2f-sin2f/cos2f)
9101      continue
	  enddo
9100    continue
	enddo
	AnVMdeflSng=rst
9002	return
      end
