 !************************************************************************************
!调用准备代码：省略变量声明
!读取地面等效水高球谐系数,负荷勒夫数
!给定时间BLH(3)
      rw=1.025d3; re=5.517d3
      !打开地面等效水高球谐系数文件
 	allocate(cnm((maxn+2)**2), stat=astat(1))
 	allocate(snm((maxn+2)**2), stat=astat(2))
	if (sum(astat(1:2)) /= 0) goto 908
      cnm=0.d0;snm=0.d0
      open(unit=8,file=clmfl,status="old",iostat=status)
      if(status/=0)goto 904
	do while(.not.eof(8))
	    read(8,*,end=901)n,m,cnm(n*(n+1)/2+m),snm(n*(n+1)/2+m)
      enddo
901   close(8)
      if(n<maxn)maxn=n
      !将等效水高球谐系数转换为位系数直接影响
      fk=3.d0*rw/re
      do n=1,maxn
        do m=0,n
            j=n*(n+1)/2+m
            cnm(j)=cnm(j)*fk/real(2*n+1);snm(j)=snm(j)*fk/real(2*n+1)
         enddo
	enddo
      flv=0.d0
      !打开全球负荷勒夫数文件
      open(unit=8,file=loadlovfl,status="old",iostat=status)
      if(status/=0)goto 908
      n=0
	do while(.not.eof(8))
        n=n+1
	  read(8,*,end=902)(flv(n,j),j=1,3)
      enddo
902   close(8)
501   continue
      if(maxn>20000)maxn=20000
      flv(1,1:3)=0.d0  !一阶项不参与计算，无需质心校正
	  call BLH_RLAT(CRS,BLH,rln)
!*********************************************************
!*****************************************************************
      subroutine GEqLoadTideCalc(rln,maxn,cnm,snm,flv,tdn,CRS,mode)
      !地面等效水高球谐综合法负荷形变场计算-计算空间点rln(3)的负荷形变
	  !mode-0直接影响,1间接影响,3负荷总影响=直接影响+间接影响
	  !当负荷为大气压时,对地面重力或地面扰动重力：负荷总影响=间接影响-直接影响
      !tdn(k)-1高程异常/大地水准面形变（mm）
      !2扰动重力变化（uGal）
      !3地倾斜北南（子午圈）分量变化（ms）
      !4地倾斜东西（卯酉圈）分量变化（ms）
      !5地面站点重力变化（uGal）
      !6地面站点东方向形变（mm）
      !7地面站点北方向形变（mm）
      !8地面站点大地高变化（mm）
      !9地面站点正（常）高变化（mm）
      !2015年7月21日,章传银
      implicit none
	integer::maxn,mode,nn,n,m,kk,i,j
	real*8::cnm((maxn+2)**2),snm((maxn+2)**2)
	real*8::rln(3),gm,tdn(9),tn(9),gr,CRS(5),cosml,sinml
	real*8::BLH(3),NFD(5),t,u
 	real*8::rr,rlat,rlon,ae,pi,RAD,flv(40000,3),dk,dh,dl
	real*8,allocatable::pnm(:),dpt1(:),dpt2(:)
!---------------------------------------------------------------------------
    	tdn=0.d0;gm=CRS(1);ae=CRS(2)
   	pi=datan(1.d0)*4.d0;RAD=pi/180.d0
 	allocate(pnm((maxn+2)**2),dpt1((maxn+2)**2),dpt2((maxn+2)**2))
      rr=rln(1);rlat=rln(2);rlon=rln(3)
      t=dsin(rlat*RAD);u=dcos(rlat*RAD)
      call RLAT_BLH(CRS,rln,BLH);call GNormalfd(BLH,NFD);gr=NFD(2)
      call BelPnmdt(pnm,dpt1,dpt2,maxn,t)
      do n=2,maxn  !一阶项不参与计算
         dh=flv(n,1);dl=flv(n,2);dk=flv(n,3);tn=0.d0
         do m=0,n
            j=n*(n+1)/2+m
            cosml=dcos(rlon*RAD*dble(m));sinml=dsin(rlon*RAD*dble(m))
	      tn(1)=tn(1)+(cnm(j)*cosml+snm(j)*sinml)*pnm(j+1)
	      tn(2)=tn(2)+(cnm(j)*cosml+snm(j)*sinml)*pnm(j+1)
	      tn(3)=tn(3)+(cnm(j)*cosml+snm(j)*sinml)*dpt1(j+1)
	      tn(4)=tn(4)+real(m)*(cnm(j)*sinml-snm(j)*cosml)*pnm(j+1)
	      tn(5)=tn(5)+(cnm(j)*cosml+snm(j)*sinml)*pnm(j+1)
	      tn(6)=tn(6)+real(m)*(cnm(j)*sinml-snm(j)*cosml)*pnm(j+1)
	      tn(7)=tn(7)+(cnm(j)*cosml+snm(j)*sinml)*dpt1(j+1)
	      tn(8)=tn(8)+(cnm(j)*cosml+snm(j)*sinml)*pnm(j+1)
         enddo
         tn=tn*dexp(real(n)*dlog(ae/rr))
         if(mode==0)then
 	      tdn(1)=tdn(1)+tn(1)
	      tdn(2)=tdn(2)+tn(2)*real(n+1)
	      tdn(3)=tdn(3)+tn(3)
	      tdn(4)=tdn(4)+tn(4)
	      tdn(5)=tdn(5)+tn(5)*real(n+1)
	      tdn(6)=0.d0
	      tdn(7)=0.d0
 	      tdn(8)=0.d0
         endif
         if(mode==1)then
 	      tdn(1)=tdn(1)+tn(1)*dk
	      tdn(2)=tdn(2)+tn(2)*real(n+1)*(-real(n+1)*dk)/real(n)
	      tdn(3)=tdn(3)+tn(3)*(dk-dh)
	      tdn(4)=tdn(4)+tn(4)*(dk-dh)
	      tdn(5)=tdn(5)+tn(5)*real(n+1)*(-real(n+1)*dk+2.d0*dh)/real(n)
	      tdn(6)=tdn(6)+tn(6)*dl
	      tdn(7)=tdn(7)+tn(7)*dl
 	      tdn(8)=tdn(8)+tn(8)*dh
         endif
         if(mode==2)then
 	      tdn(1)=tdn(1)+tn(1)*(1.d0+dk)
	      tdn(2)=tdn(2)+tn(2)*real(n+1)*(1.d0-real(n+1)*dk/real(n))
	      tdn(3)=tdn(3)+tn(3)*(1.d0+dk-dh)
	      tdn(4)=tdn(4)+tn(4)*(1.d0+dk-dh)
	      tdn(5)=tdn(5)+tn(5)*real(n+1)*(1.d0+(2.d0*dh-real(n+1)*dk)/real(n))
	      tdn(6)=tdn(6)+tn(6)*dl
	      tdn(7)=tdn(7)+tn(7)*dl
 	      tdn(8)=tdn(8)+tn(8)*dh
         endif
	enddo
	tdn(1)=tdn(1)*gm/rr/gr*1.d3
	tdn(2)=tdn(2)*gm/rr**2*1.0e8
	tdn(3)=tdn(3)*gm/rr**2/gr/RAD*36.d5
	tdn(4)=tdn(4)*gm/rr**2/gr/u/RAD*36.d5
	!tdn(5)=tdn(5)*gm/rr**2*1.0e8
	tdn(6)=-tdn(6)*gm/rr/gr/u*1.d3
 	tdn(7)=-tdn(7)*gm/rr/gr*1.d3
	tdn(8)=tdn(8)*gm/rr/gr*1.d3
	tdn(5)=tdn(2)-tdn(8)*0.3086d0 !!!!
      tdn(9)=tdn(8)-tdn(1)
      deallocate(pnm,dpt1,dpt2)
 	return
	end
