      subroutine RentGField(maxn,BLH,cnm,snm,gvm,CRS)
!输入maxn－模型重力场最大阶数
!输入cnm、snm－模型重力场位系数
!输出gvm(8)-高程异常(m)、空间异常(mGal)、扰动重力(mGal)、南向、西向垂线偏差(s)、重力梯度(E)、北向、东向水平重力梯度(E)，扰动位(m2/s2)
!2018年6月6日,章传银
      implicit none
	integer::maxn
	real*8::cnm((maxn+2)**2),snm((maxn+2)**2)
	real*8::cosml,sinml,tn(5),gr,pi,RAD,rln(3),BLH(3),NFD(5)
	integer n,m,kk,i,astat(6)
	real*8 rr,rlat,rlon,gm,ae,gvm(9),CRS(5),u,t,sp
	real*8,allocatable::pnm(:),dpt1(:),dpt2(:)
!---------------------------------------------------------------------------
   	gm=CRS(1);ae=CRS(2);pi=datan(1.d0)*4.d0; RAD=pi/180.d0
 	allocate(pnm((maxn+2)**2), stat=astat(1))
 	allocate(dpt1((maxn+2)**2), stat=astat(2))
 	allocate(dpt2((maxn+2)**2), stat=astat(3))
	if (sum(astat(1:3)) /= 0) goto 904
      call GNormalfd(BLH,NFD);call BLH_RLAT(CRS,BLH,rln)
	rr=rln(1);rlat=rln(2);rlon=rln(3);gr=NFD(2)
      t=dsin(rlat*RAD);u=dcos(rlat*RAD)
      call BelPnmdt(pnm,dpt1,dpt2,maxn,t)
	do i=1,9
	  gvm(i)=0.d0
	enddo
	do n=2,maxn
        tn=0.d0
	  do m=0,n
          kk=n*(n+1)/2+m
          cosml=dcos(dble(m)*rlon*RAD);sinml=dsin(dble(m)*rlon*RAD)
	    tn(1)=tn(1)+(cnm(kk)*cosml+snm(kk)*sinml)*pnm(kk+1)
	    tn(2)=tn(2)+(cnm(kk)*cosml+snm(kk)*sinml)*dpt1(kk+1)
	    tn(3)=tn(3)+dble(m)*(cnm(kk)*sinml-snm(kk)*cosml)*pnm(kk+1)
	    tn(4)=tn(4)+(cnm(kk)*cosml+snm(kk)*sinml)*dpt2(kk+1)
	    tn(5)=tn(5)+dble(m**2)*(cnm(kk)*cosml+snm(kk)*sinml)*pnm(kk+1)
        enddo
        sp=dexp(dble(n)*dlog(ae/rr))
	  gvm(1)=gvm(1)+tn(1)*sp
	  gvm(2)=gvm(2)+tn(1)*sp*dble(n-1)
	  gvm(3)=gvm(3)+tn(1)*sp*dble(n+1)
	  gvm(4)=gvm(4)+tn(2)*sp
	  gvm(5)=gvm(5)+tn(3)*sp
	  gvm(6)=gvm(6)+tn(1)*sp*dble(n+1)*dble(n+2)
	  gvm(7)=gvm(7)+tn(4)*sp
	  gvm(8)=gvm(8)+tn(5)*sp
	enddo
	gvm(9)=gm/rr*gvm(1) !m2/s2
	gvm(1)=(gm/rr/gr)*gvm(1)  !m
	gvm(2)=(gm/rr**2)*gvm(2)*1.d5 !mGal
	gvm(3)=(gm/rr**2)*gvm(3)*1.d5 !mGal
	gvm(4)=gm/rr**2/gr*gvm(4)/RAD*36.d2 !s
	gvm(5)=gm/rr**2/gr*gvm(5)/u/RAD*36.d2 !s
	gvm(6)=gm/rr**3*gvm(6)*1.d9 !E
	gvm(7)=gm/rr**3*gvm(7)*1.d9 !E
	gvm(8)=-gm/rr**3*gvm(8)/u**2*1.d9 !E
      deallocate(pnm,dpt1,dpt2)
904	return
      end
