!************************************************************************************
!调用准备代码：省略变量声明
!读取海潮潮高球谐系数,负荷勒夫数
!给定时间IY,IM,ID,ih,imn,sec
      !读FES2004海潮模型文件oloadcsfl
      open(unit=8,file=oloadcsfl,status="old",iostat=status)
      if(status/=0)goto 902
      read(8,'(a)') line
      read(8,'(a)') line   !FES2004(11)
      read(8,'(a)') line   !FES2004(11)
      i=0;fes=0.d0
      do while(.not.eof(8))  
	  read(8,*,end=903)st(1),fh,(st(j),j=2,3),tmp,tmp,tmp,tmp,(st(j),j=4,7)
        if(st(2)<mm+1)then
          i=i+1;st(1)=st(1)*1.d3;fes(i,1:7)=st(1:7);
        endif
	enddo
903   close(8)
      nn=i
	!读负荷勒夫数
      open(unit=8,file=loadlovfl,status="old",iostat=status)
      if(status/=0)goto 902
      n=0
	do while(.not.eof(8))
        n=n+1
	  read(8,*,end=904)(flv(n,j),j=1,3)
      enddo
904   close(8)
      if(mode==1)flv(1,1:3)=0.d0  !一阶项不参与计算，无需质心校正
	  call BLH_RLAT(CRS,BLH,rln)
      call CAL2JD (IY,IM,ID,djm,j)
      djm=djm+real(ih)/24.d0+real(imn)/1440.d0+real(sec)/864.d2 !GPS_MJD
	  dt=djm+2400000.5d0+19.d0/864.d2!; td=td-TAI_UTC(td)/864.d2 !UTC_MJD 处理跳秒
	  call OLoadDFluCS(td,cnm,snm,maxn,fes,nn)   !!cnm,snm直接影响位系数
!*********************************************************
!*****************************************************************
      subroutine LTideFlu(rln,maxn,cnm,snm,flv,tdn,CRS,mode)
      !由地面等效水高球谐系数变化（cnm,snm）计算空间点rln(3)的负荷形变
	  !不同时间时间，由该时间的cnm、snm体现！调用OLoadDFlu.f90实现
	  !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
