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