fortran求相关及t检验程序：改变数据但是通过检验的站点数不变，不知道什么地方错

菜菜黑

program test
implicit none

real skt(8,12,66)
real staa(8,12,66),t(160,63),stat(160,63),rain6(160,63),rain7(160,63),rain8(160,63),rjpbfl(160,63)
real,dimension(8,12)::s1,jfangcha,s2,ave,r,tt
real,dimension(160)::jfangcha2,sumrain,averain
real summerrain(160,63),cj(17,63),cjr(63),avecj,jfcrain,jprain(63)
real jfc,starain(63),avejprain,sum1
integer i,j,imo,it,k


open(11,file='d:\paper\output\2.grd',form='binary')         !2月地温数据
open(13,file='d:\paper\data\r1606.txt')
open(14,file='d:\paper\data\r1607.txt')
open(15,file='d:\paper\data\r1608.txt')

! 读数据
do i=1,8
  do j=1,12
          do it=1,66
           read(11) skt(i,j,it)
enddo;enddo;enddo
close(11)
read(13, *) ((rain6(i,k), i=1,160), k=1,63)
read(14, *) ((rain7(i,k), i=1,160), k=1,63)
read(15, *) ((rain8(i,k), i=1,160), k=1,63)
close(13)
close(14)
close(15)


!高原地温场读取                  
     do j=1,12
          do i=1,8                  
           ave(i,j)=sum(skt(i,j,1:66))/66.0          
          enddo
         enddo
     do j=1,12
          do i=1,8   
          s2=0.0
       do it=1,66
            s2(i,j)=s2(i,j)+(skt(i,j,it)-ave(i,j))**2
         jfangcha(i,j)=sqrt((s2(i,j))/66.0)
            enddo
          enddo
         enddo

    do it=1,66
     do j=1,12
          do i=1,8  
           staa(i,j,it)=(skt(i,j,it)-ave(i,j))/jfangcha(i,j)
    enddo;enddo;enddo

!夏季降水写到一个数组中
do i=1,160
do it=1,63
    summerrain(i,it)=rain6(i,it)+rain7(i,it)+rain8(i,it)
enddo
enddo

!160站降水距平百分率
do i=1,160
    sumrain(i)=sum(summerrain(i,1:63))
        averain(i)=sumrain(i)/63.0
       
enddo
do i=1,160
  do it=1,63
   rjpbfl(i,it)=(summerrain(i,it)-averain(i))*100/averain(i)
   enddo
enddo

!挑出长江地区的站点并写入数组
cj(1:14,1:63)=summerrain(54:67,1:63)
cj(15:16,1:63)=summerrain(72:73,1:63)
cj(17,1:63)=summerrain(77,1:63)


!
do it=1,63
cjr(it)=sum(cj(1:17,it))
end do
avecj=sum(cjr(1:63))/63.0
do it=1,63
jprain(it)=(cjr(it)-avecj)/avecj
enddo

sum1=0.0
do it=1,63
sum1=sum1+(cjr(it)-avecj)**2
end do

jfc=sqrt(abs(sum1/63.0))

do it=1,63
starain(it)=(jprain(it)-avecj)/jfc
enddo
!计算相关系数
do j=1,12
do i=1,8
sum1=0.0
do it=1,63
sum1=sum1+staa(i,j,it+3)*starain(it)
end do
r(i,j)=sum1/63.0
end do
end do
!t检验
  do i=1,8
    do j=1,12
      tt(i,j)=r(i,j)*sqrt(61.0)/sqrt(1-r(i,j)**2)
    enddo
  enddo
  k=0
  do i=1,8
    do j=1,12
        if( tt(i,j)>=1.98.or.tt(i,j)<=-1.98)then
        k=k+1
        endif
    enddo
  enddo
  print*,k


!把数据写成文件
open(41,file='d:\paper\output\r2.txt')
write(41,"(f10.4)") ((r(i,j),i=1,8),j=1,12)
close(41)
open(42,file='d:\paper\output\tt2.txt')
write(42,"(f10.4)") ((tt(i,j),i=1,8),j=1,12)
close(42)


end program

lqouc

请先阅读提问的智慧

菜菜黑

lqouc 发表于 2014-5-9 18:18
请先阅读提问的智慧

你好，我想求青藏高原各月地表温度与长江中下游降水的相关，并且对其做t检验，在程序中我输出了通过检验的站点数，想通过这个数据来大致判断相关性再进行其他工作。我检查了数据的读入以及公式感觉没有问题，运行也没有错误，但是当我带入不同的月份去计算，通过的站点数都是一样的。我认为这个程序肯定是有错误的，但是找不出来，麻烦你帮我看一下。

lqouc

菜菜黑 发表于 2014-5-9 19:04
你好，我想求青藏高原各月地表温度与长江中下游降水的相关，并且对其做t检验，在程序中我输出了通过检验 ...

这样看不出什么的。
你测试了几个月？是相关系数一样还是只是过检站数一样？站数虽然一样但是也可能并不重合。
所以你这么肯定地说错误有待商榷，话说程序是你写的么？不是的话可以询问编写者。
这样的计算在grads里是极其简单的。

菜菜黑

lqouc 发表于 2014-5-9 20:10
这样看不出什么的。
你测试了几个月？是相关系数一样还是只是过检站数一样？站数虽然一样但是也可能并不 ...

测试了12月到6月，问了老师，他说这样应该是有错误的，但是他也不知道错在哪，程序是我在论坛里找的不同的程序来改了写到这个程序里的。这个计算在grads里怎么做我真的不知道，以后会好好去学，目前想先把这个程序的问题找出来，请您指点一下。

lqouc

确实感觉你的程序是有问题，sum函数直接可以用？
高原地温那段感觉循环有问题，其他的一时看不出来，你可以先debug一下出错的位置
我这两天出差，只有手机没法给你看

菜菜黑

lqouc 发表于 2014-5-10 00:02
确实感觉你的程序是有问题，sum函数直接可以用？
高原地温那段感觉循环有问题，其他的一时看不出来，你可 ...

debug还没学会怎么用谢谢你指出可能错误的地方，我再检查。

lqouc

不知道你弄出来了没有哈。就之前我说过的那部分，我觉得你的循环有问题，按照我的理解应该改成下面这样的。后面的内容可能还有错，但是不好确定了。
高原地温场读取                  
     do j=1,12
          do i=1,8                  
           ave(i,j)=sum(skt(i,j,1:66))/66.0         
          enddo
         enddo
     do j=1,12
          do i=1,8   
          s2=0.0
       do it=1,66
            s2(i,j)=s2(i,j)+(skt(i,j,it)-ave(i,j))**2
            enddo
         jfangcha(i,j)=sqrt((s2(i,j))/66.0)
          enddo
         enddo
这样的问题你要先尝试确定出错的位置。可以自学一下debug，用彭国伦的教材，最多二十分钟上手。

菜菜黑

修改了之后通过的站点变多了，但是还是所有月份都是那么多站。debug参考了彭国伦的书，不知道是自己没掌握还是什么，没找出问题。

菜菜黑

lqouc 发表于 2014-5-11 13:59
不知道你弄出来了没有哈。就之前我说过的那部分，我觉得你的循环有问题，按照我的理解应该改成下面这样的。 ...

修改了之后通过的站点变多了，但是还是所有月份都是那么多站。debug参考了彭国伦的书，不知道是自己没掌握还是什么，没找出问题。