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

質問

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

回答No.6

f272 · ベストアンサー

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))
とすることもできます。

お礼

ありがとうございました。
少し混乱しています。下のように書きましたが、さらに値がおかしくなってしまいました。

do doy = 1,365

ave = 0.0
var = 0.0

do i=1,ndata(doy)
ave = ave + ndatax(doy,i)
end do

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

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

stdevplus(doy) = temp(doy) + stdev(doy)
stdevminus(doy) = temp(doy) - stdev(doy)
! stdevplus(doy)/2 = temp(doy) + stdev(doy)/2
! stdevminus(doy)/2 = temp(doy) - stdev(doy)/2

write(6,*) doy, temp(doy), ndata(doy),stdevplus(doy),stdevminus(doy)
! write(11,*) doy,',',temp(i),',',ndata(i)
end do

回答No.7

f272

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

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

お礼

申し訳ありませんでした。僕のグラフ化の操作ミスでした。

元データでの力技でだしたグラフと一致しました。ありがとうございました。
色々とご指導いただき、参考にさせていただきます。
まだまだ創造しにくいもので困っています。

次に気温ではなく、日照時間などのプログラムをつくるつもりなのですが、これは日にちを旬別に平均しなければなりません。少し自力でやってみまして、わからないところがあればまた質問しようと思っております。もし見かけることがありましたら、ぜひまたアドバイスいただければ幸いです。

本当にありがとうございました。

回答No.5

f272

(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
として使用する。

お礼

ありがとうございました。
(３)について次回から参考にさせていただきます。

実は、まだうまくいかなくて、おそらく平均の方は合っているような気がするのですが、標準偏差が明らかに違います。どうかアドバイスください。

program sapporo_kikouchi

INTEGER :: sum, no, point
INTEGER :: year, mon, day, data
INTEGER :: doy
REAL,dimension(365) :: temp
REAL,dimension(365) :: stdev
REAL,dimension(365) :: stdevplus,stdevminus
INTEGER,dimension(365) :: ndata
REAL,dimension(365,33) :: ndatax
REAL :: lon, lat
REAL :: ave, var
CHARACTER*4 yyyy
CHARACTER*5 sssss
ndata(:)=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)
ndata(doy) = ndata(doy) + 1
if( ndata(doy) > 33 ) then
write(6,*)"Error"
stop
end if
ndatax(doy,ndata(doy)) = data/10.0

!　ndata(doy)は一年のうち何日目か,個数を表している

enddo

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

do doy = 1,365

ave = 0.0
var = 0.0

do i=1,ndata(doy)

ave = ave + ndatax(doy,i)

if( ndata(doy) == 0 ) then
temp(doy) = -99999.9
else
temp(doy) = ave/ndata(doy)
var = var + ( ndatax(doy,i) - temp(doy))**2
stdev(doy) = sqrt(var/(ndata(doy)-1))
endif

end do

stdevplus(doy) = temp(doy) + stdev(doy)
stdevminus(doy) = temp(doy) - stdev(doy)
! stdevplus(doy)/2 = temp(doy) + stdev(doy)/2
! stdevminus(doy)/2 = temp(doy) - stdev(doy)/2

write(6,*) doy, temp(doy), ndata(doy),stdevplus(doy),stdevminus(doy)
! write(11,*) doy,',',temp(i),',',ndata(i)
end do

stop
end program

---------------------サブルーチン文省略----------------------------------------------

回答No.4

f272

気がついたところ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)
のはずです。

お礼

すごいです。できました。ありがとうございます。
文字制限の関係からエラーメッセージを載せれなかったのですが、次から工夫して載せるようにします。
ただ、2つの疑問がでました。
私の一番初めの質問した際のプログラムと今回のプログラムだと、値に誤差が出ています。エクセルで確認したところ、後者の方の値が間違っていました。
どうしてなのかでしょうか。。
また、標準偏差をつけてみたのですが、値が大きく違いました。これは、私のミスだと思いますが。
一応、今のプログラムを下に載せて起きます。もし、なにかお気づきになれば、アドバイスいただけると幸いです。

INTEGER :: sum, no, point
INTEGER :: year, mon, day, data
INTEGER :: doy
REAL,dimension(365) :: temp
REAL,dimension(365) :: stdev
REAL,dimension(365) :: stdevplus,stdevminus
INTEGER,dimension(365) :: ndata
INTEGER,dimension(365,33) :: ndatax
REAL :: lon, lat
REAL :: ave, var
CHARACTER*4 yyyy
CHARACTER*5 sssss
ndata(:)=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)
ndata(doy) = ndata(doy) + 1
if( ndata(doy) > 33 ) then
write(6,*)"Error"
stop
end if
ndatax(doy,ndata(doy)) = data/10.0

!　ndata(doy)は一年のうち何日目か,個数を表している

enddo

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

do doy = 1,365

ave = 0.0
var = 0.0

do i=1,ndata(doy)

ave = ave + ndatax(doy,i)
var = var + ( ndatax(doy,i) - temp(doy))**2
end do

if( ndata(doy) == 0 ) then
temp(doy) = -99999.9
else
temp(doy) = ave/ndata(doy)
stdev(doy) = sqrt(var/(ndata(doy)-1))
endif

stdevplus(doy) = temp(doy) + stdev(doy)
stdevminus(doy) = temp(doy) - stdev(doy)
　　write(6,*) doy, temp(doy), ndata(doy),stdevplus(doy),stdevminus(doy)
end do
-------------------subroutine文省略--------------------------------------------

回答No.3

f272

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回が最大数ではなくなります。
その場合は，色々とやり方はあるだろうが，それが必要になった時にやり方を学んでください。

お礼

連絡が遅くなってしまい申し訳ありませんでした。色々考えていたのですが、うまくまわりません。したにプログラムを載せます。ポカミスでしょうか？
real dimension(365、33) :: ndatax
としています。データは1976~2008の33年間。1月1日は33個以上にいくことは恐らくないと考えています。

program sapporo_kikouchi

INTEGER :: sum, no, point
INTEGER :: year, mon, day, data
INTEGER :: doy
REAL dimension(365) :: temp, ndata
REAL dimension(365,33) :: ndatax
REAL :: lon, lat
REAL :: ave
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)
ndata(doy) = ndata(doy) + 1
ndatax(doy,ndata(doy)) = data/10.0

!　ndata(doy)は一年のうち何日目か,個数を表している

enddo

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

do doy = 1,365

ave = 0.0

do i=1,ndata(doy)

ave = ave + ndatax(doy,i)

end do

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

write(6,*) i, temp(doy), ndata(doy)
! 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

回答No.2

f272

データが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

で計算出来るでしょう。

お礼

お返事が遅くなってしまい本当に申し訳ありませんでした。
ご丁寧に本当にありがとうございます。なんとなくわかってきました。
今日、色々試してみまして、再度お礼or質問追加させていただきます。
取り急ぎ、一度感謝をお伝えいたします。

回答No.1

f272

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

補足

program sapporo_kikouchi

INTEGER :: sum, no, point
INTEGER :: year, mon, day, data
INTEGER :: doy
REAL,dimension(365) :: temp, ndata, ndatax
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)
ndatax(doy,i) = data/10.0
temp(doy) = temp(doy) + ndatax(doy,i)
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

ndatax(i) = ndatax(i) - temp(i)

write(6,*) i, temp(i), ndata(i), nadatax(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

と書きましたが、nadatax(doy,i) = data/10.0 の部分でエラーが出ます。

まだ2次元配列がうまく理解できません。どう考えれいいのか。どう書けばいいのか。。
も少しアドバイスください。