サンプル（このままでは動かない）
入力データ：3次元データbuf(経度, 緯度, 時間)
出力: 線形トレンド, trd(経度, 緯度)
trdの単位はbufの時間間隔に依存するので、適宜bufの実際の時間刻みをかけてやる必要がある。
real,allocatable :: buf(:,:,:)
real,allocatable :: sum(:,:),ave(:,:),sdv(:,:),trd(:,:),cnt(:,:)
real,allocatable :: P(:),Q(:) !FOR LINEAR REGRESSION
allocate(sum(im,jm),ave(im,jm),sdv(im,jm),trd(im,jm),buf(im,jm,nm))
allocate(cnt(im,jm))
print '(A)','INITIALIZE'
buf(:,:,:)=0.0
sum(:,:) =0.0
ave(:,:) =0.0
sdv(:,:) =0.0
trd(:,:) =0.0
cnt(:,:) =0.0
print '(A)','SUM'
TIME_LOOP: do n=ns,ne
LAT_LOOP: do j=1,jm
LON_LOOP: do i=1,im
if(lon(i)>=wlon.and.lon(i)<=elon.and.lat(j)>=slat.and.lat(j)<=nlat)then
if(sdata(i,j,n)>0.0)then
buf(i,j,n)=sdata(i,j,n)
sum(i,j) =sum(i,j)+sdata(i,j,n)
cnt(i,j) =cnt(i,j)+1.0
endif !sdata
endif !lon lat
end do LON_LOOP
end do LAT_LOOP
end do TIME_LOOP
print '(A)','AVE'
do j=1,jm
do i=1,im
if(cnt(i,j)>10.0)then
ave(i,j)=sum(i,j)/cnt(i,j)
else
ave(i,j)=-999.
endif
end do !i
end do !j
print '(A)','SDV'
do j=1,jm
do i=1,im
if(cnt(i,j)>10.0)then
do n=ns,ne
if(buf(i,j,n)>0.0)then
sdv(i,j)=sdv(i,j)+(buf(i,j,n)-ave(i,j))**2
end if
end do !n
sdv(i,j)=sqrt(sdv(i,j)/(cnt(i,j)-1.0))
else
sdv(i,j)=-999.
end if
end do !i
end do !j
print '(A)','TREND'
do j=1,jm
do i=1,im
if(cnt(i,j)>10.0)then !cnt
ND=0
do n=ns,ne
if(buf(i,j,n)>0.0)then !buf
ND=ND+1
endif !buf
enddo !n
allocate(P(ND),Q(ND))
ntemp=0
do n=ns,ne !n
if(buf(i,j,n)>0.0)then !buf
ntemp=ntemp+1
P(ntemp)=float(ntemp)
Q(ntemp)=buf(i,j,n)
endif !buf
enddo !n
CALL LSQM(P,Q,ND,A,B)
trd(i,j)=B
deallocate(P,Q)
else !cnt
trd(i,j)=-999.
end if !cnt
end do !i
end do !j
!***********************************************************************
!* LSQM.F
!* F:\0.Computer\TOOLBOX\BASIC_STATISTICS\回帰直線
!* CODED BY A.MANDA 1997.8
!* TASK
!C Linear Regression
!* PARAMETERS
!* (1)P : INPUT DATA [X-COOD.] (R)(INPUT)
!* (2)Q : INPUT DATA [Y-COOD.] (R)(INPUT)
!* (3)N : TOTAL DATA NUMBER (I)(INPUT)
!* (4)A : Y-SECTION OF LINE (R)(OUTPUT)
!* (5)B : GRADIENT OF LINE (R)(OUTPUT)
!* METHOD
!* SLAVE SUBROUTINE(S)
!* REMARK
!* REFERENCE(S)
!***********************************************************************
SUBROUTINE LSQM(P,Q,N,A,B)
!* INPUT
integer,intent(in)::N
REAL,intent(in):: P(N),Q(N)
!* OUTPUT
REAL,intent(inout):: A,B
TP=0.0
TQ=0.0
DO 10 I=1,N
TP=TP+P(I)
TQ=TQ+Q(I)
10 CONTINUE
PM=TP/FLOAT(N)
QM=TQ/FLOAT(N)
TPQ=0.0
TPP=0.0
DO 20 I=1,N
TPQ=TPQ+(P(I)-PM)*(Q(I)-QM)
TPP=TPP+(P(I)-PM)**2
20 CONTINUE
B=TPQ/TPP
A=QM-B*PM
RETURN
END SUBROUTINE LSQM