      subroutine StableVdeform(infl1,infl2,infl3,outfl,dp)
      !高程变化格网地面稳定性变化检测
      !infl1,infl2,infl3-高程变化率、地面高程变化梯度和地形坡度格网文件
	  !dp(1~3):地面高程变化率权,地面高程变化梯度权,地形坡度权
	  !dp(4~6):地面高程变化率指数,地面高程变化梯度指数,地形坡度指数
	  !dp(7~8):地面高程变化率时序平均值和标准差
	  !dp(9~10):地面高程变化梯度时序平均值和标准差
      implicit none
	character*800::infl1,infl2,infl3,outfl
 	character*8000::line,line0
	integer i,j,sn,len,nlon,nlat,knd,mode,astat(8)
	real*8::q1,q2,q3,hd1(6),hd2(6),hd3(6),rec(800),dp(40),x1,x2,x3,qq,rdv(4),dhd,dhd1
	real*8,allocatable::dchg(:,:),grad(:,:),dtm(:,:),rst(:,:)
	integer::status=0
!---------------------------------------------------------------------
      qq=dp(1)+dp(2)+dp(3);q1=dp(1)/qq;q2=dp(2)/qq;q3=dp(3)/qq
      open(unit=8,file=infl1,status="old",iostat=status)
      if(status/=0)goto 902
      open(unit=10,file=infl2,status="old",iostat=status)
      if(status/=0)then
         close(8);goto 902
      endif
      open(unit=12,file=infl3,status="old",iostat=status)
      if(status/=0)then
         close(8);close(10);goto 902
      endif
      read(8,'(a)') line0
      call PickRecord(line0,len,rec,sn)
      hd1(1:6)=rec(1:6)
      read(10,'(a)') line
      call PickRecord(line,len,rec,sn)
      hd2(1:6)=rec(1:6)
      read(12,'(a)') line
      call PickRecord(line,len,rec,sn)
      hd3(1:6)=rec(1:6)
      dhd=0.d0
      do i=1,6
        dhd=dhd+dabs(hd1(i)-hd2(i))
      enddo
      dhd1=0.d0
      do i=1,6
        dhd1=dhd1+dabs(hd1(i)-hd3(i))
      enddo
      if(dhd>1.d-5 .or. dhd1>1.d-5)then  !格网规格不同
         close(8);close(10);close(12);dp(21)=1.d0;goto 902
      endif
	nlat=nint((hd1(4)-hd1(3))/hd1(6))
	nlon=nint((hd1(2)-hd1(1))/hd1(5))
	hd1(5)=(hd1(2)-hd1(1))/real(nlon)
	hd1(6)=(hd1(4)-hd1(3))/real(nlat)
      hd2=hd1;hd3=hd1
 	allocate(dchg(nlat,nlon), stat=astat(1))
 	allocate(grad(nlat,nlon), stat=astat(2))
 	allocate(dtm(nlat,nlon), stat=astat(3))
      allocate(rst(nlat,nlon), stat=astat(4))
	if (sum(astat(1:4)) /= 0) then
          close(8);close(10);close(12);goto 902
	endif
 	do i=1,nlat
	   read(8,*,end=905)(dchg(i,j),j=1,nlon)
	   read(10,*,end=905)(grad(i,j),j=1,nlon)
	   read(12,*,end=905)(dtm(i,j),j=1,nlon)
      enddo
905   close(8);close(10);close(12)
      call StatGrid(dtm,nlat,nlon,rdv)
      dp(11)=rdv(1);dp(12)=rdv(2)
 	do i=1,nlat
        do j=1,nlon
          x1=(dchg(i,j)-dp(7))/dp(8)
          x2=(grad(i,j)-dp(9))/dp(10)
          x3=(dtm(i,j)-dp(11))/dp(12)
          x1=x1*(dabs(x1)+1.d-16)**(dp(4)-1.d0)
          x2=x2*(dabs(x2)+1.d-16)**(dp(5)-1.d0)
          x3=x3*(dabs(x3)+1.d-16)**(dp(6)-1.d0)
          rst(i,j)=q1*x1+q2*x2+q3*x3
        enddo
      enddo
      open(unit=12,file=outfl,status="replace")
      write(12,'(a)')trim(line0)
 	do i=1,nlat
	   write(12,'(15F12.4)')(rst(i,j),j=1,nlon)
      enddo
      close(12)
904   deallocate(dchg,grad,dtm,rst)
902   continue
      end
