!************************************************************************************
!调用准备代码：省略变量声明
!由地面高程模型dtm和地面大地高模型hgt，计算地形质量补偿密度mu及其一、二阶导数ux,ux,uxx,uyy
      open(unit=8,file=infl1,status="old",iostat=status)
      if(status/=0)goto 902
      read(8,'(a)') line
      call PickRecord(line,len,rec,sn)
      hd(1:6)=rec(1:6)
	nlat=nint((hd(4)-hd(3))/hd(6))
	nlon=nint((hd(2)-hd(1))/hd(5))
	hd(5)=(hd(2)-hd(1))/real(nlon)
	hd(6)=(hd(4)-hd(3))/real(nlat)
 	allocate(dtm(nlat,nlon), stat=astat(1))
 	allocate(hgt(nlat,nlon), stat=astat(2))
 	allocate(mu(nlat,nlon), stat=astat(3))
 	allocate(uxx(nlat,nlon), stat=astat(4))
 	allocate(uyy(nlat,nlon), stat=astat(5))
 	allocate(ux(nlat,nlon), stat=astat(6))
 	allocate(uy(nlat,nlon), stat=astat(7))
	if (sum(astat(1:7)) /= 0) then
          close(8);goto 902
	endif
 	do i=1,nlat
	   read(8,*,end=903)(dtm(i,j),j=1,nlon)
         do j=1,nlon
           if(dtm(i,j)<0)dtm(i,j)=0.d0
         enddo
      enddo
903   close(8)
      open(unit=10,file=infl2,status="old",iostat=status)
      if(status/=0)goto 904
      read(10,'(a)') line
      call PickRecord(line,len,rec,sn)
      hd1(1:6)=rec(1:6)
      if(sum(hd1-hd)>1.d-5)then  !格网规格不同
         close(10);dp(11)=1.d0;goto 904
      endif
 	do i=1,nlat
	   read(10,*,end=905)(hgt(i,j),j=1,nlon)
      enddo
905   close(10)
 	do i=1,nlat
        BLH(1)=hd(3)+(real(i)-0.5d0)*hd(6)
        do j=1,nlon
	    BLH(2)=hd(1)+(real(j)-0.5d0)*hd(5)
          BLH(3)=CGrdPntD2(BLH(2),BLH(1),hgt,nlat,nlon,hd)
          call BLH_RLAT(CRS,BLH,rln);rr=rln(1);zch=dtm(i,j)
	    mu(i,j)=tp*zch*(1.d0+zch/rr+zch**2/rr**2/3.d0)*1.d-5  !值太大，超过系统缺省值9900.0，统一乘1.d-5,在结果补回去
        enddo
      enddo
      call HGradientTwo(mu,hgt,nlat,nlon,hd,uxx,uyy)
      call HGradientOne(mu,hgt,nlat,nlon,hd,ux,uy)
!*********************************************************
      real*8 function LRlfBLHdPotent(BLH,dtm,hgt,uxx,uyy,nlat,nlon,hd,dr,CRS)
      !按严密球面积分公式计算地球外空间点BLH(3)的扰动位局部地形影响
	  !dtm(nlat,nlon),hgt(nlat,nlon)地面高程模型和地面大地高模型
	  !uxx(nlat,nlon),uyy(nlat,nlon)地形补偿密度二阶水平梯度北方向、东方向分量
      !dr-积分半径m
!-------------------------------------------------------------
      implicit none
	integer::i,j,nlat,nlon,i0,j0,ni,nj,astat(5)
	real*8::dr,mu(nlat,nlon),hgt(nlat,nlon),uxx(nlat,nlon),uyy(nlat,nlon)
	real*8::hd(6),gp,pi,RAD,ds,mdr,tt,r0,r1,rst
	real*8::BLH(3),XYZ(3),BLH0(3),XYZ0(3),BLH1(3),XYZ1(3),rln1(3),BLH2(3),XYZ2(3)
	real*8 CGrdPntD2,HCmpnsSinguDp,L0,L1,up,u1
	real*8::CRS(5),ge,re,tp
      integer::kfn !奇异积分核函数精化参数2*kfn
!-----------------------------------------------------------------
      ge=6.67428d-11!万有引力常数
      re=6371000.79d0!地球平均半径
      tp=2.67d3!地形密度
	  kfn=2 !可变
      pi=datan(1.d0)*4.d0;RAD=pi/180.d0;
      BLH0=BLH;BLH0(3)=CGrdPntD2(BLH(2),BLH(1),hgt,nlat,nlon,hd)!计算点正下方地面点位置
      call BLH_XYZ(CRS,BLH,XYZ);call BLH_XYZ(CRS,BLH0,XYZ0)
      r0=dsqrt(XYZ0(1)**2+XYZ0(2)**2+XYZ0(3)**2)
      up=CGrdPntD2(BLH(2),BLH(1),mu,nlat,nlon,hd)
      rst=0.d0;mdr=r0*hd(5)*RAD*dcos(BLH(1)*RAD)/2.d0 !奇异点判断
      ni=nint(dr/r0/RAD/hd(6)+1.d0) !积分半径dr对应的地面格网数
      nj=nint(dr/r0/RAD/hd(5)/dcos(BLH(1)*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)
          L0=dsqrt((XYZ1(1)-XYZ0(1))**2+(XYZ1(2)-XYZ0(2))**2+(XYZ1(3)-XYZ0(3))**2)
          if(L0>dr)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)
          if(L1<mdr)then!计算奇异积分
             rst=rst+HCmpnsSinguDp(BLH,mu,hgt,uxx,uyy,nlat,nlon,hd,i,j,kfn,CRS)
             goto 9101 
          endif
          u1=mu(i,j)
          tt=1.d0-2.d0*(L0/r1/2.d0)**2
          ds=hd(5)*hd(6)*RAD**2*dcos(rln1(2)*RAD)
          rst=rst+(u1-up)*ge/L1*ds*r1**2
9101      continue
	  enddo
9100    continue
	enddo
	HCmpnsBLHdPotent=rst*1.d5 !补回1.d5
9002	return
      end
!--------------------------------------------------------------------------------
      real*8 function HCmpnsSinguDp(BLH,mu,hgt,uxx,uyy,nlat,nlon,hd,i0,j0,m,CRS)
      !细化核函数，计算BLH点的扰动位局部地形补偿的奇异积分
      !m-核函数细化为m*m
      !i0,j0-奇异点格网位置
!-------------------------------------------------------------
      implicit none
	integer::m,i,j,nlat,nlon,i0,j0,ni,nj
	real*8::mu(nlat,nlon),hgt(nlat,nlon),uxx(nlat,nlon),uyy(nlat,nlon)
	real*8::hd(6),gp,pi,RAD,ds,mdr,tt,r0,rr,r1,rst,rv,tx,ty
	real*8::BLH(3),XYZ(3),BLH0(3),XYZ0(3),BLH1(3),XYZ1(3),rln1(3),BLH2(3),XYZ2(3)
	real*8 CGrdPntD2,zch,zch1,L0,L1,lon,lat,up,u1
	real*8::CRS(5),ge,re,tp
!-----------------------------------------------------------------
      ge=6.67428d-11!万有引力常数
      re=6371000.79d0!地球平均半径
      tp=2.67d3!地形密度
      rv=hd(5)/dble(m);pi=datan(1.d0)*4.d0;RAD=pi/180.d0;gp=ge*tp
      lat=hd(3)+real(i0-1)*hd(6);lon=hd(1)+real(j0-1)*hd(5)!格网左下角经纬度
      BLH0=BLH;BLH0(3)=CGrdPntD2(BLH(2),BLH(1),hgt,nlat,nlon,hd)!计算点正下方地面点位置
      call BLH_XYZ(CRS,BLH,XYZ);call BLH_XYZ(CRS,BLH0,XYZ0)
      rr=dsqrt(XYZ(1)**2+XYZ(2)**2+XYZ(3)**2)
      r0=dsqrt(XYZ0(1)**2+XYZ0(2)**2+XYZ0(3)**2)
      up=CGrdPntD2(BLH(2),BLH(1),mu,nlat,nlon,hd)
      rst=0.d0;mdr=r0*rv*RAD*dcos(BLH(1)*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)
          L0=dsqrt((XYZ1(1)-XYZ0(1))**2+(XYZ1(2)-XYZ0(2))**2+(XYZ1(3)-XYZ0(3))**2)
          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
          u1=CGrdPntD2(BLH1(2),BLH1(1),mu,nlat,nlon,hd)
          tt=1.d0-2.d0*(L0/r1/2.d0)**2
          rst=rst+(u1-up)*ge/L1*ds*r1**2
9101      continue
          if(L1<mdr)then !计算奇异积分
            tx=CGrdPntD2(BLH1(2),BLH1(1),uxx,nlat,nlon,hd)
            ty=CGrdPntD2(BLH1(2),BLH1(1),uyy,nlat,nlon,hd)
            rst=rst+(tx+ty)*ge*ds*r1**5*dsqrt(ds/pi)/rr**2/6.d0
          endif
	  enddo
9100    continue
	enddo
	HCmpnsSinguDp=rst
9002	return
      end
