      subroutine GNormalfd(CRS,BLH,NFD)
!输入CRS(5)-gm,a,j2,w,1/f
!输入BLH(3)-(纬度,经度,大地高),单位:度小数,m
!输出NFD(5)-正常重力位,正常重力,正常重力梯度,正常重力线方向,正常梯度方向
      implicit none
 	real*8::BLH(3),NFD(5),val(5),rln(3),pn(40),dp1(40),dp2(40)
      integer::i,j,n
	real*8::CRS(5),djn(40),pi,RAD,gm,ae,ww,rr,t,u,vdn(5),atr
!---------------------------------------------------------------------
      pi=datan(1.d0)*4.d0;RAD=pi/180.d0
	gm=CRS(1);ae=CRS(2);ww=CRS(4);vdn=0.d0
      call BLH_RLAT(CRS,BLH,rln)
      rr=rln(1);t=dsin(rln(2)*RAD);u=dcos(rln(2)*RAD)
      call LegPn_dt2(pn,dp1,dp2,40,t)
      do n=1,20
        atr=dexp(dble(2*n)*dlog(ae/rr))*djn(2*n)
        vdn(1)=vdn(1)+atr*pn(2*n)
        vdn(2)=vdn(2)+atr*pn(2*n)*dble(2*n+1)
        vdn(3)=vdn(3)+atr*dp1(2*n)
        vdn(4)=vdn(4)+atr*pn(2*n)*dble((2*n+1)*(n+1))
        vdn(5)=vdn(5)+atr*dp2(2*n)
      enddo
      val(1)=gm/rr*(1.d0-vdn(1))+(ww*rr*u)**2/2.d0
      val(2)=-gm/rr**2*(1.d0-vdn(2))+ww**2*rr*u**2 !径向正常重力
      val(3)=-gm/rr**2*vdn(3)+ww**2*rr*u*t         !南向正常重力
      val(4)=-2.d0*gm/rr**3*(1.d0-vdn(4))+(ww*u)**2
      val(5)=-gm/rr**3*vdn(5)+ww**2*(2.d0*t**2-1.d0)
      NFD(4)=datan(val(3)/val(2))/RAD*36.d2
      NFD(5)=datan(val(5)/val(4))/RAD*36.d2
      val(2)=dsqrt(val(2)**2+val(3)**2)
      val(4)=dsqrt(val(4)**2+val(5)**2)
      NFD(1)=val(1);NFD(2)=val(2);NFD(3)=val(4)
      end
