Fortranについての質問です。下のプログラムは、ある地点（今回は1

Question

Fortranについての質問です。下のプログラムは、ある地点（今回は14161~14163）の33年間(1976~2008）の気温の平均を欠損値を考慮してだしているはずです。うるう年の判定はおまけで書いています。
この平均値に標準偏差もつけるプログラムにしたいのですが、どう書いていいのか手詰まり中です。どうかアドバイスください。

---------データの一部（14162_temp1997.csv)-------------------------------------

14162 1997 1 1 141.332 43.0583 -10
14162 1997 1 2 141.332 43.0583 18
14162 1997 1 3 141.332 43.0583 -2

---------プログラム------------------------------------------------------------

program sapporo_kikouchi

INTEGER :: sum, no, point
     INTEGER :: year, mon, day, data
     INTEGER :: doy
     REAL,dimension(365) :: temp, ndata
     REAL :: lon, lat
     CHARACTER*4  yyyy
     CHARACTER*5  sssss
     ndata(:)=0.0
     temp(:)=0.0

do ispot=14161,14163
     write(sssss,"(i5)") ispot
   do iwork=1976, 2008
     write(yyyy,"(i4.4)") iwork

open(50, file=''//sssss//'_temp'//yyyy//'.csv', status='old',iostat=io)
      if (io < 0) cycle
!     write(6,*) ispot iwork

do i = 1,366

read(50,*,iostat=io)  id,year,mon,day,lon,lat,data
        if(io < 0) exit
        if(mon==2 .AND. day==29) then
           cycle
        endif

call date2doy(year,mon,day,doy)
        temp(doy) = temp(doy) + data/10.0
        ndata(doy) = ndata(doy) + 1

!　doyは一年のうち何日目かを表している

enddo

close(50)
  enddo   !!! end of year loop
  enddo

do i=1,365

if( ndata(i) == 0 ) then
      temp(i) = -99999.9
    else
      temp(i)=temp(i)/ndata(i)
    endif

write(6,*) i, temp(i), ndata(i)
!   write(11,*) i,',',temp(i),',',ndata(i)

end do

stop
end program

subroutine date2doy(iy,im,id,idoy)
 INTEGER,dimension(12) :: nday
 INTEGER :: uruu  !!uruu=1: うるう年、uruu=0: 通常の年
 
 uruu=0
 DATA nday /31,28,31,30,31,30,31,31,30,31,30,31/

if(mod(iy,4)==0 .AND. mod(iy,100)/=0) then
   uruu=1
 endif
  if(mod(iy,1000)==0) then
   uruu=1
 endif

!! うるう年も無視する
 itotal = 0
 if( im /= 1 )then
 do m=1, im-1
   itotal = itotal + nday(m)
 enddo
 endif
 idoy = id + itotal

! write(6,*) iy,im,id, idoy

return
 end

f272 · Accepted Answer

1: do doy = 1,365
2: ave = 0.0
3: var = 0.0
4: do i=1,ndata(doy)
5: ave = ave + ndatax(doy,i)
6: temp(doy) = ave/ndata(doy)
7: var = var + ( ndatax(doy,i) - temp(doy))**2
8: stdev(doy) = sqrt(var/(ndata(doy)-1))
9: end do
10: write(6,*) doy, temp(doy), ndata(doy),stdevplus(doy),stdevminus(doy)
11: end do
これはあなたのプログラムを抜き出したものですが，平均値に関係しているのは
2で初期化して，4から9のループの中で5で加算し，6で割り算をしています。
まあ，これでも計算はうまくいきますが，6で割り算はループのの外でやる方が効率的でしょう。
しかし標準偏差の計算はひどいことになっています。
3で初期化して，4から9のループの中で7で加算し，8で割り算をしています。
ところがループの中で加算すべき物はデータ引く平均値の2乗であるのに，7はループの中ですから加算する時点では平均値の計算途中であって，まだ平均値になっていません。それを加算したところで計算がうまくいくはずがないのです。それで前回
> 標準偏差を計算するときには，すでに平均値が分かっていなければなりません。
> ところが，あなたのプログラムでは
> temp(doy) = ave/ndata(doy)
> で平均値を求めるところよりも前で，標準偏差の計算らしき部分に入っています。これでは計算がおかしく> なって当然でしょう。
と言ったのに理解されなかったようです。それではどうするかと言えば，平均値を計算するためのループが終わってから，あらためて標準偏差を計算するためのループをまわしてください。
!    ave = 0.0
!    var = 0.0
!    do i = 1,ndata(doy)
!      ave = ave + ndatax(doy,i)
!    end do
!    temp(doy) = ave/ndata(doy)
!    do i = 1,ndata(doy)
!      var = var + (ndatax(doy,i)-temp(doy))**2
!    end do
!    stdev(doy) = sqrt(var/(ndata(doy)-1))
こんな感じですね。明示的にループをまわさずにsum関数を使って
!   temp(doy) = sum(ndatax(doy,:ndata(doy)))/ndata(doy)
!   stdev(doy) = sqrt(sum((ndatax(doy,:ndata(doy))-temp(doy))**2)/(ndata(doy)-1))
とすることもできます。

f272 · Answer

データを読み込んだらそれを
data/10.0
という形で使っておしまいにするのではなく，例えば
datax(doy,i)=data/10.0
のように配列に保存しておけば，ループを抜けて平均値を出した後でも計算に利用できるよね。

f272 · Answer

データが10個あってその標準偏差を求めるのであれば
real dimension(10) :: x1
と定義して，このx1(1:10)について計算すればいいわけです。
この問題の場合は，求める標準偏差は各doy毎にあるわけだから
real dimension(doyの最大数,各doy毎のデータの個数の最大値) :: x2
と定義するわけです。
doyの最大数は簡単で365ですね。
各doy毎のデータの個数の最大値は，読み込まれたデータを見なければ分からないのですが，何とか考えてみましょう。
データは
do ispot=14161,14163
do iwork=1976, 2008
do i = 1,366
この3重ループの中で読んでいますが，このうちiは1から365までのdoyのどれかに対応しますから，各doy毎のデータの個数とは関係ありません。
ispotとiworkのループは合計で(14163-14161+1)*(2008-1976+1)=3*33=99回あって，各doy毎のデータの個数は絶対にこれを超えることがないことが分かるでしょう。
ということで
real dimension(365,99) :: ndatax
と定義します。変数名もあなたの使おうとしているものに変えました。

次にデータを配列に代入する部分です。
3重ループの中で，まず
read(50,*,iostat=io) id,year,mon,day,lon,lat,data
で読んで
call date2doy(year,mon,day,doy)
でdoyを出してから，配列に代入するわけですが
ndatax(doy,?)=data/10.0
の?はどうすればよいでしょうか？各doy毎に何番目のデータかを表していますから，簡単ですね。
ちゃんとndata(doy)で個数を数えていますね。これが使えるでしょう。

ndata(doy) = ndata(doy) + 1
ndatax(doy,ndata(doy)) = data/10.0
temp(doy) = temp(doy) + ndatax(doy,ndata(doy))

でうまくいくことが理解できますか？

これでデータがちゃんとndatax(doy,:)に確保できました。個数もndata(doy)で分かります。
そして各doy毎に平均値tempを出すには

do doy=1,365
ave=0.0
do i=1,ndata(doy)
ave=ave+ndatax(doy,i)
end do
temp(doy)=ave/ndata(doy)
end do

で出来ますし，平均値tempが分かれば，各doy毎の標準偏差stdevは

do doy=1,365
var=0.0
do i=1,ndata(doy)
var=var+(ndatax(doy,i)-temp(doy))**2
end do
stdev(doy)=sqrt(var/(ndata(doy)-1))
end do

で計算出来るでしょう。

f272 · Answer

ndata(doy) = ndata(doy) + 1
ndatax(doy,ndata(doy)) = data/10.0
temp(doy) = temp(doy) + ndatax(doy,ndata(doy))

でうまくいくことが理解できますか？

と書いたが
temp(doy) = temp(doy) + ndatax(doy,ndata(doy))
は，ここで加算しなくても大丈夫だった。
do doy=1,365
ave=0.0
do i=1,ndata(doy)
ave=ave+ndatax(doy,i)
end do
temp(doy)=ave/ndata(doy)
end do
で加算していますからね。

それから
real dimension(365,99) :: ndatax
と言ったが，もちろんこれは一つのcsvファイルの中に同じ日付が重複していないことを仮定しています。
もし重複しているのであれば99回が最大数ではなくなります。
その場合は，色々とやり方はあるだろうが，それが必要になった時にやり方を学んでください。

f272 · Answer

気がついたところ1
> 色々考えていたのですが、うまくまわりません。
ではなくて，どのようなエラーになるのかをちゃんと記述すること。

気がついたところ2
> real dimension(365、33) :: ndatax
ではなくて，
real,dimension(365、33) :: ndatax
とすること。このままではコンパイルエラーになる。

気がついたところ3
> 1月1日は33個以上にいくことは恐らくないと考えています。
ではなくて，ちゃんと確認する。
ndata(doy) = ndata(doy) + 1
の次の行として
if (ndata(doy)>=33) then
  print*,"Error: dimension overflow"
  stop
endif
としておけば確実です。

気がついたところ4
> write(6,*) i, temp(doy), ndata(doy)
ではなくて
write(6,*) doy, temp(doy), ndata(doy)
のはずです。

f272 · Answer

(1)
私の一番初めの質問した際のプログラムと今回のプログラムだと、値に誤差が出ています。エクセルで確認したところ、後者の方の値が間違っていました。

これは
INTEGER,dimension(365,33) :: ndatax
のせいでしょう。この変数は
ndatax(doy,ndata(doy)) = data/10.0
という使い方をしていることからも分かるように実数型の変数であるべき。

(2)
また、標準偏差をつけてみたのですが、値が大きく違いました。これは、私のミスだと思いますが。

標準偏差を計算するときには，すでに平均値が分かっていなければなりません。
ところが，あなたのプログラムでは
temp(doy) = ave/ndata(doy)
で平均値を求めるところよりも前で，標準偏差の計算らしき部分に入っています。これでは計算がおかしくなって当然でしょう。

(3)
その他，プログラムを作る上で気を付けた方が良いこと。
A. fortranの変数名について，整数型はi,j,k,l,m,nで始ます変数名にすること，また実数型であればそれ以外の文字で始めること。これは初期のfortranからの伝統です。
B. implicit none 文を使用して必ず変数は宣言すること。
C. 読みやすさのため、プログラム の構造上ひとまとまりとなる文をそのほかの文に比べて，字下げすること。
D. 例えば
REAL,dimension(365) :: temp
や
do ispot=14161,14163
のようにプログラムの中に数字をそのまま書きこまないこと。これは
integer,parameter :: ndays_per_year = 365
のようにして，この変数を使う。もっと言えばそれを
module array_size
integer,parameter :: ndays_per_year = 365
end module array_size
のようにモジュールで定義しておいて，メインプログラムでは
use array_size , only : ndays_per_year
として使用する。

f272 · Answer

> さらに値がおかしくなってしまいました。

どうして値がおかしいことが分かるの？全然状況がつかめません。

Fortranについての質問です。下のプログラムは、ある地点（今回は1

データを読み込んだらそれを

この回答への補足

データが10個あってその標準偏差を求めるのであれば

ndata(doy) = ndata(doy) + 1

気がついたところ1

(1)

1: do doy = 1,365

> さらに値がおかしくなってしまいました。

似たような質問が見つかりました

関連するカテゴリからQ&Aを探す

デイリーランキングこのカテゴリの人気デイリーQ&Aランキング

マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング