2011 長崎大学
FORTRAN 講習会
筑波大学生命環境科学研究科地球科学専攻 ょ 津秀敏
九 大学総合理 学府大気海洋環境クケゾヘ学専攻ょ 山 慶人
屋大学環境学研究科 ょ ょ ょ ょ 木 太郎
Linux
端曒 開 方
゚ハモォヴクョル → ゚ェセキモ → GNOME端曒
メルスホ 表示 く 便利
こ 端曒 カブルチを入力 る
端曒 閉 方
端曒 ” exit ” 入力 る
よく使うカブルチ mkdir
mkdir make directory 略 ある directory Windows いう ころ ネァャジ
こ あり mkdir ネァャジを くるカブルチ ある 試 端曒 よう
入力 実行 よう！
mkdir test
ls ll
ls ll 入力 るこ いるタ゛ヤェダモ 何 入 いる 確認 る 両方
試 よう！
る 先ほ 作 test いうdirectory ある
ls ll 何 違う 考え よう！
ls -a
入力 よう！
先ほ 比べ 何 違う 考え よう！
cd
change directory 略 タ゛ヤェダモを変える 移動 る 使うカブルチ
Windows ゜ベヴグ：
例え word ネ゙゜ャを見る
home → チキュベルダ → word → 自 いネ゙゜ャ
いう作業をェモッェ らやる こ 作業をこ cd いうカブルチを入力 行う
カブルチを入力 先ほ 作 タ゛ヤェダモ 移動 よう！
cd ~/test
※ここ ~ homeを意味 る
ls 入力 よう！
い 確認 ら
mkdir test2 入力
ls
入力 よう！先ほ 比べ 何 違う 確認 よう！
確認 ら
cd ..
入力 よう！
“cdょ ..” 入力 るこ 前 directory 移動 るこ る
次
cd ~/test/test2
入力 よう！
る い home らtest2 いうdirectory 移動 る！
home 戻り い場合
cd ~
入力 れ 戻れる
pwd
present working directory 略 いるタ゛ヤェダモ こ を示 カブルチ
gedit , vi
gedit ゾキケダ゠タ゛シ viを使 ハュエメヘを作成 る
gedit ベペ帳 よう ゜ベヴグ vi ょ く 実際 端曒 ハュエメヘを書く
いう゜ベヴグ vi 設定を い 不便 初心者 gedit を使 ほう
いい 思われる 慣れ らvi 楽
試 両方使 よう！
gedit 起動方法 端曒
gedit 入力
゚ェセキモ → ゾキケダ゠タ゛シ
ら開く
保 り閉 る ベペ帳 様
vi 起動方法 操作方法
端曒
vi ネ゙゜ャ
例ょ vi testfile 入力
入力 方
“ i ”を押 挿入ペヴチ ら描 始 る
これを い 書け い 注意！
保 方法
保 閉 る場合：
Escを押 後 こ 操作 より挿入ペヴチを解除 る shiftを押 らzz 入力
保 い 閉 る場合：
Escを押 後 ” :q! ” 入力
先ほ 作成 test いうdirectory 移動 cd ~/test gedit vi 適当 ネ
゙゜ャを作 よう！ mv
作成 ネ゙゜ャ を別 directory 移 カブルチ
例
先ほ 作成 ネ゙゜ャをtest らtest2 移動 る場合
test いうdirectory いるこ を”pwd” 確認 ら カブルチを入力 よう！ mv filename ~/test/test2
※filename 先ほ 自 作成 ネ゙゜ャ 前を入力
cd test2 ls
入力 移動 確認 よう！
cp
いるdirectory 別 directory あるネ゙゜ャ をカヌヴ るカブルチ
例 先ほ test2 移動 ネ゙゜ャをtest カヌヴ よう！
pwd test2 いるこ を確認 後 次 カブルチを入力
cp filename ~/test
cd .. ls
カヌヴ 確認
awk
タヴシ 処理を行うカブルチ 列 タヴシを り出 こ る
具体的 例 後ほ 明
rm
作 ネ゙゜ャ を消去 るカブルチ
rm filename
先ほ く ネ゙゜ャを消 よう！
※rm 消去 ネ゙゜ャ ガプ箱 行 復元不可能 注意
紹 基曓的 カブルチ ある こ カブルチ ある
各自 調べ よう
復習
Home fortran いうdirectoryを作成 ら 中 28dec2011 いう directoryを作成 よう！
home 戻 適当 ネ゙゜ャを作り 先ほ 作 28dec2011 いうdirectory 移動 よう！
移動 を確認 ら 移動 ネ゙゜ャを消去 よう！
FORTRAN
ハュエメプルエ 実行 流れ
vi, gedit ゾキケダ゠タ゛シを使 ハュエメヘを書 filename.f90 いう 前 保
端曒
gfortran filename.f90
f95 filename.f90 打
゠メヴ 出 けれ
ょ ./a.out 打
Fortran90 ハュエメヘ 書 方
ハュエメヘ べ 英語 ゚ャネ゙パッダ 数 記号 書く
カベルダ 日曓語 OK
ハュエメヘ 英語 ゚ャネ゙パッダ 数 デ゜ネル(-) ゚ルジヴケカ゚ _ 書
くこ 暷初 文 ゚ャネ゙パッダ るこ
1行 長 暷大 半角 132文
2行 わ る場合 行 暷後 ょ ＆ょ を けるこ
! 書く 行 れ 降 カベルダ れる
ハュエメヘ 基曓的 構成
ハュエメヘ文 program
暗黙 型宣言 無効 implicit none 宣言文 integer, real
入力文 read
実行文 計算式
出力文 write
゠ルチハュエメヘ文 end program
実際 ハュエメヘ
や り実際 ハュエメプルエ 実行 る 習得 近道！ 簡 ハュエメプル
エ らや よう！
ensyu1.f90
name 部 自 前を入れ よう！
program ensyu1 ← program文 write(*,*) ‘name’ ← write文 stop ←stop文
end program ensyu1 ← end program文 明
ハュエメヘ 前をensyu1 いう 前 る
文 や数値を書 込 (*,*) る 画面 ケェモヴル 直接出力 表示 る (6,*)
よい こ 場合 name いう文 を出力 る ’ ’ くくるこ 文 いう
こ を示 ’ ’ょ くくら 場合 変数を意味 る
こ 後 ハュエメヘを書い 実行 れ いこ を意味 る
ハュエメヘ 終わりを表 文 ハュエメヘ 1行目 るこ
りや いよう よう 書くこ を勧 る よう 省略 書くこ 可
能
write(*,*)ょ ‘name’ end
ensyu2.f90
次 ょ ヤレゟャ 次 面積を計算 るハュエメヘを作 よう！
面積ょ S=πr2 program ensyu2
implicit none real :: r,s
real,parameter :: pi=3.141592 write(*,*) ‘input r’
read(5,*) r s=pi*r**2.0 write(*,*) ’S=’,s stop
write文 write( * , * )
ィッカ 右側 装置番号
側 書式を設定 る(format文)
詳 いこ 後 明
←ょ 宣言文
←ょ read文 write文
end program ensyu2 明
宣言文：
これを書くこ よ 宣言文や り 間違いを指摘 くれる
real 変数 実数 小数 あるこ を宣言 る こ 場合 r s いう変数 実数
あるこ を宣言 いる 整数 場合 integer る
parameter文 pi ドメベヴシ値を設定 やる 小数 real いる
read文 write文：
read文 数値を 込 命 (5,*) る キヴピヴチ ら入力 値を 取
る いう命 る
面積 計算を行う ’*’ 乗算を意味 ’**’ 2乗を意味 る ここ 注意 い
数式 よう ゜カヴャ 意味 く s pi*r**2.0を 入 る いう意味
実行 らキヴピヴチ r 具体的 数値を入れ よう！
ensyu3.f90
次 If文を使 ハュエメヘを作成 よう！
if文 明
～ ら い
If (条件) then
処理ょ ←条件 真 場合
End if
例：
if (x > 60) then write(*,*) ‘OK’ end if
x 60より大 い らOK 出力 れる
～ ら い けれ い
if (条件 ) then
ょ 処理 ょ ←条件 真
else
ょ 処理 ょ ←条件 偽
end if 例：
if (x > 60) then write(*,*) ‘OK’ else
write(*,*) ‘Not Good’ end if
ういう意味 る 考え よう！
program ensyu3 implicit none integer :: y,x write(*,*) ‘input x’ read(5,*) x
if (x > 0) then y = x+10 end if
write(*,*) ‘y=’,y stop
end program ensyu3
例
条件 条件 れ 外 場合：
if (x > 80) then
write(6,*) ‘Very Good’ else if (x > 60) then write(6,*) ‘Good’ else
write(6,*) ‘Not Good’ end if
条件 条件 を 時 満 場合：
if (x > 80 .and. y > 80) then write(6,*) ‘Very Good’ end if
ら 一方 条件を満 よい場合： if ( x > 60 .or. y>60 ) then
‘OK’ end if
一致 値 場合：
if ( x == 100 ) then ‘Excelent’
end if
ensyu4.f90
do文 いう構文を使 ハュエメヘ do文 ある操作を繰り返 いう命 を る
1~100 和を求 均値を求 るハュエメヘ
! 1~100 足 算 均値を求 るハュエメヘ
program ensyu4
! 宣言文
implicit none
integer :: i ! iを整数型 宣言 real :: n ! nを実数型 宣言
! 1~100 和 計算 n=0
do i=1,100 n=n+i end do
! 画面 均値を求 出力
write(*,*) "ave=",n/100.
stop
end program ensyu4
←ょ !を けるこ カベルダ るこ
る
←ょ do文
゜ベヴグ
i=1 n=0+1=1 i=2 n=1+2=3 i=3 n=3+3=6 ン
ン ン
こ 作業をi=100 繰り返
ensyu5.f90
配列 open文を使 ハュエメヘ
program ensyu5
implicit none
! integer 変数 整数 あ を指定 例 .1,2,3,4,....
! real 変数 実数 つ 少数 あ を指定 例 .1.23 2.31 3.67...
integer i,j,mx,my
parameter(mx=3,my=3)
real,dimension(mx,my) :: field
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! open input file
open(10,file='test1.dat',status='old')
! read input file
do j=1,my
read(10,*) (field(i,j),i=1,mx)
end do
! close input file
close(10)
!
! open output file
open(20,file='result1.dat',status='new')
! write 'read data' to output file
do j=1,my
write(20,*) (field(i,j),i=1,mx)
end do
! close output file
close(20)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end program ensyu4
使用タヴシ test1.dat 1
2 3 4 5 6 7 8 9
ensyu6.f90
配列を使 計算
program ensyu6
implicit none
integer i,j,mx,my
parameter(mx=3,my=3)
real,dimension(mx,my) :: A, B, C, D, E
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! タ !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! open test1 data file
open(10,file='test1.dat',status='old')
! read test1 data file
do j=1,my
read(10,*) (A(i,j),i=1,mx)
end do
! close test1 data file
close(10)
! open test2 data file
open(20,file='test2.dat',status='old')
! read test2 data file
do j=1,my
read(20,*) (B(i,j),i=1,mx)
end do
! close test2 data file
close(20)
! open test3 file
open(30,file='test3.dat',status='old')
! read test3 file
do j=1,my
read(30,*) (C(i,j),i=1,mx)
end do
! close test3 file
11 12 13 14 15 16 17 18 19
21 22 23 24 25 26 27 28 29
31 32 33 34 35 36 37 38 39
test2.dat
test3.dat
test4.dat 使用タヴシ test1.dat 先 ほ タ ヴ シ を使用
close(30)
! open test4 file
open(40,file='test4.dat',status='old')
! read test4 file
do j=1,my
read(40,*) (D(i,j),i=1,mx)
end do
! close test4 file
close(40)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 平均を計算 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
E = (A+B+C+D)/4
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 計算結果を別 ァイ 入力 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! open output file
open(50,file='heikin.dat',status='unknown')
!60 format(' ',3f10.2)
! write 'read data' to output file
do j=1,my
write(50,*) (E(i,j),i=1,mx)
end do
! close output file
close(50)
end program ensyu6
気象タヴシ
NCEP/NCAR再解析タヴシ
NCEP/NCAR National Centers for Environmental Prediction/Natinal Center for Atmospheric Research 略
形式：net-CDF 実際 開い 見るこ いタヴシ
水 解能：2.5度×2.5度
鉛直ヤレゟャ：17層 1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, 20, 10hPa
期間：1948/01/01 - 現在
取得キ゜ダ：http://www.esrl.noaa.gov/psd/data/reanalysis/reanalysis.shtml ERA-40
European Center for Medium-Range Weather Forecasts Re-Analysis 略
形式：net-CDF
水 解能：2.5度×2.5度
鉛直ヤレゟャ：23層：1000, 925, 850, 775, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, 20, 10, 7, 5, 3, 2, 1hPa
期間：1957/09/01 – 2002/08/31
取得キ゜ダ：http://database.rish.kyoto-u.ac.jp/arch/era40/
※ 登録 必要
気象庁 ゚ベジケタヴシ
取得キ゜ダ：http://www.jma.go.jp/jma/menu/report.html Climate Prediction Center タヴシ
月別ゾヤカネェクョル タヴシ ある
取得キ゜ダ：ftp://ftp.cpc.ncep.noaa.gov/wd52dg/data/indices/tele_oldindex.nh
GrADS
描画
先ほ 紹 ビヴヘヒヴグ らNCEP/NCAR再解析タヴシ 2011 気温 daily
dataをジゞルュヴチ よう！
端曒 grads
入力 Enterキヴを押
先ほ ジゞルュヴチ ネ゙゜ャ 中身を見 る
sdfopen air.2011.nc q ctlinfoょ
実際 タヴシを開い 1000hPa面 気温 タヴシを描画 る
sdfopenょ air.2011.nc set time 01jan2011 set lev 1000
d air
これをeps いう形式 ネ゙゜ャ 保 る
print 01jan2011T.eps
GrADSを終了 る quit
作業をexecネ゙゜ャを作 一気 実行 る
gedit vi カブルチを書い 01jan2011.exec いう 前を け 保 る 01 jan2011.exec
sdfopenょ air.2011.nc set time 01jan2011 d air
print 01jan2011T.eps
GrADSを開い カブルチを実行 る exec 01jan2011.exec
GrADSを閉 端曒 convert いうカブルチを用い eps らjpg 変換 る convert 01jan2011T.eps 01jan2011T.jpg
東 大 ビヴヘヒヴグ
http://wind.gp.tohoku.ac.jp/index.php?%B8%F8%B3%AB%BE%F0%CA%F3%2FGr ADS%2FGrADS%A4%CETips
を参考 いろいろや よう！
タヴシ り出
GrADS 描画 け く タヴシをト゜ヂモ形式 り出 こ る
例 先ほ 描画 タヴシを り出 る
01jan2011T1000.exec sdfopenょ air.2011.nc set time 01jan2011 set lev 1000
set fwrite -le 01jan2011T1000.bin set gxout fwrite
d air
disable fwrite reinit
次 演習 使うタヴシを作成 る 次 作業を る
NCEP/NCAR daily dataをジゞルュヴチ る ょ air.1989.nc uwnd.1989.nc vwnd.1989.nc
exec fileを作成 実行 る mkbin1989feb.exec
sdfopen air.1989.nc
set lat 10 80
set lon 80 280
set lev 150
set time 04feb1989
set fwrite -le 04feb1989T150.bin
set gxout fwrite
d air
disable fwrite
reinit
sdfopen air.1989.nc
set lat 10 80
set lon 80 280
set lev 250
set time 04feb1989
set fwrite -le 04feb1989T250.bin
set gxout fwrite
d air
disable fwrite
reinit
sdfopen uwnd.1989.nc
set lat 7.5 82.5
set lon 77.5 282.5
set lev 200
set time 04feb1989
set fwrite -le 04feb1989U200.bin
set gxout fwrite
d uwnd
disable fwrite
reinit
sdfopen vwnd.1989.nc
set lat 7.5 82.5
set lon 77.5 282.5
set lev 200
set time 04feb1989
set fwrite -le 04feb1989V200.bin
set gxout fwrite
d vwnd
disable fwrite
reinit
気象要素 計算
カモアモドメベヴシ 計算
カモアモドメベヴシ f =2 Ω sinφ coli.f90
!! コ パ タ f を計算 グ !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1. 宣言文
! mx: 緯度方向 グ ッ 数
! my: 経度方向 グ ッ 数
! m: 領域
! omega: Ω 7.292×10^(-5) [1/s]
! lat0: 最初 緯度
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
integer i,j,mx,my,m,lat0,omega
parameter(mx=29,my=81,m=mx*my,omega=7.292*10**(-5),lat0=10)
real,dimension(mx,my) :: f
real rad,pai
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
pai=4.0*atan(1.0)
rad=pai/180
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 惑星渦度 計算
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do i=1,my
do j=1,mx
f(i,j)=2*omega*sin((lat0+2.5*(i-1))*rad)
end do
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
open(10,file='f.dat',status='unknown')
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do i=1,my
do j=1,mx
write(10,110) f(i,j)
end do
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
110 format(3f12.9)
end
温 計算
温 pt.f90
!! 温位 気圧変化を計算 グ !!
!
! 温位 気圧変化
! d θ/dp
!
! < グ 主 流 >
! 1. 宣言文
! 2. タ 込
! 3. 計算 , タ 書 込
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! 用意 タ
! 気圧 タ
! 温度 タ
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1. 宣言文
! mx: 緯度方向 グ ッ 数
! my: 経度方向 グ ッ 数
! m: 領域
! pr: 気圧
! tem: 温度
! dp: 鉛直方向 気圧面 中央差分
! 例 :200hPa 面 d θ/dp を求 たい
! →dp = pB - pA
! = 250-150=100[hPa]
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
integer i,j,mx,my,m,dp,pA,pB
parameter(mx=81,my=29,pA=150,pB=250,dp=pB-pA,m=mx*my)
real,dimension(mx,my) :: temA,temB,onniA,onniB
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 2. タ 込
! 2.1. ァイ を開い タを 込
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! < 温度 >
! 150hPa
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!open(30,file='t150.dat',status='old')
!read(30,*) ((temA(i,j),i=1,mx),j=1,my)
open(30,file='04feb1989T150.bin',status='old',form='unformatted',access='direct',recl=
m*4)
read(30,rec=1) ((temA(i,j),i=1,mx),j=1,my)
close(30)
! 250hPa
!open(40,file='t250.dat',status='old')
!read(40,*) ((temB(i,j),i=1,mx),j=1,my)
open(40,file='04feb1989T250.bin',status='old',form='unformatted',access='direct',recl=
m*4)
read(40,rec=1) ((temB(i,j),i=1,mx),j=1,my)
close(40)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 3. 温位 計算
! 150hPa
do j=1,my
do i=1,mx
onniA(i,j)=temA(i,j)*(1000/pA)**(0.2860)
end do
end do
close(50)
! 250hPa
do j=1,my
do i=1,mx
onniB(i,j)=temB(i,j)*(1000/pB)**(0.2860)
end do
end do
close(60)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 4. 温位勾配 計算
open(100,file='onnikobai.dat',status='unknown')
do j=1,my
do i=1,mx
write(100,*) (onniB(i,j)-onniA(i,j))/(dp*100)
end do
end do
close(100)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! format 文
110 format(f10.8)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end
球面 ける渦度 計算
球面 ける渦度 ζ vort.f90
! 中央差 球面 渦度を求 るハュエメヘ
program uzudo implicit none
integer,parameter :: nx=83 ! 経度方向 エモッチ数 integer,parameter :: ny=31 ! 緯度方向 エモッチ数ょ real :: u(nx,ny) ! 東西風 (m/s)
real :: v(nx,ny) ! 風 (m/s) real(8) :: vort(nx,ny) ! 相対渦度 (1/s) real :: lon(nx),lat(ny) ! 経度,緯度 (degree) integer :: i,j,m
m=nx*ny ! 込 タヴシ エモッチ数
! ネ゙゜ャを開い 東西風 風 タヴシを 込
open(10,file='feb1989U200.bin',status='old',form='unformatted',access='direct',recl=m*4 )
open(20,file='feb1989V200.bin',status='old',form='unformatted',access='direct',recl=m*4 )
read(10,rec=1) ((u(i,j),i=1,nx),j=1,ny) read(20,rec=1) ((v(i,j),i=1,nx),j=1,ny)
close(10) close(20)
!ょ 緯度 経度 計算 do i=1,nx
lon(i)=77.5+2.5*(i-1) end do
do j=1,ny
lat(j)=7.5+2.5*(j-1)
end do
! subroutine rvortをcall る
call rvort(nx,ny,u,v, & ょ ょ ょ ! input vort) ょ ょ ょ ! output
!ょ rvort.dat いうネ゙゜ャを開い 経度ょ 緯度ょ 渦度 を書 込
!ょ 緯度 ら高緯度 け タヴシを書 込
! i doャヴハ 経度方向 タヴシを書く
! j doャヴハ 緯度方向 タヴシを書く open(30,file='rvort.dat')
do j=2,ny-1 do i=2,nx-1
write(30,'(2f8.2,f12.4)') lon(i),lat(j),vort(i,j)*10**5 end do
end do close(30)
stop
end program uzudo
! 相対渦度を計算 るキノャヴスル subroutine rvort(nx,ny,u,v,vort) implicit none
integer, intent(in) :: nx,nyょ ょ ょ ょ ょ ょ ょ ょ ! 入力用引数 指定 整数型 real, intent(in) :: u(nx,ny),v(nx,ny)ょ ょ ! 入力用引数 指定 実数型 real(8), intent(out):: vort(nx,ny) ! 出力用引数 指定 実数型 real :: a,pi,dx,dy
integer :: j,i
a=6370000.0 ! 地球 半径 pi=4.0*atan(1.0) ! 周率
! 渦度 計算
do j=2,ny-1
dx=a*cos((j+3)*pi/73.)*2*pi/144 ! 経度 エモッチ間 距 dy=a*cos((j+3)*pi/73.)*pi/73 ! 緯度 エモッチ間 距
do i=2,nx-1
vort(i,j) = (v(i+1,j)-v(i-1,j))/(2.0*dx) &
-(u(i,j+1)*cos((j+4)*pi/73)-u(i,j-1)*cos((j+2)*pi/73))/(2.0*dy) end do
end do end subroutine
中央差 方
ンょ ン ン ン ン ン ン ン ン ン ン ょ ン ン ン ン ン ン ン ン ン ン ン ン ン ン ン ン ン ン ン ン ン ン ンょ ン ン ン ン ン ン ン ン ン ン ンょ ン ン ン ン ン ン ン ン ン ン
• ： 風 タヴシ
• ：東西風 タヴシ
ょ ：相対渦度 タヴシ
• タヴシを用い 中央差 相対渦度 値を求
る
述 タヴシ 方 書 方 よう 点線 赤
領域 相対渦度 タヴシを求 る
球面 ける相対渦度：ょ ζ
ハュエメヘ い
∂v：v(i+1,j)-v(i-1,j) ょ acosθλ：2.0*dx
∂ucosθ：u(i,j+1)*cos((j+4)*pi/73)-u(i,j-1)*cos((j+2)*pi/73) acosθ∂θ：2.0*dy
偏差 計算
偏差ょ X’ = anomaly.f90
program anomaly
implicit none
integer,parameter :: nx=83 ! 経度 グ ッ 数
integer,parameter :: ny=31 ! 緯度 グ ッ 数
integer,parameter :: t=31 ! 日数
integer,parameter :: m=nx*ny
real :: dta(m,t)
real :: dt(m,t)
real :: dth(m,t)
integer :: i,j
open(10,file='climate_dec_H500.bin',status='old',form='unformatted',access='direct',recl
=m*4)
open(20,file='jan2011.bin',status='old',form='unformatted',access='direct',recl=m*4)
do j=1,t
read(10,rec=j)(dta(i,j),i=1,m)
end do
do j=1,t
read(20,rec=j)(dt(i,j),i=1,m)
end do
close(10)
close(20)
open(30,file='jan2011H500anomaly.bin',status='unknown',form='unformatted',access='d
irect',recl=m*t*4)
do j=1,t
do i=1,m
dth(i,j)=dt(i,j)-dta(i,j)
end do
end do
do j=1,t
write(30,rec=j) (dth(i,j),i=1,m)
end do
stop
end program anomaly
統計ハュエメヘ
タヴシ 作成
ナケダエメヘ用
Climate Prediction Centerょ
ょ ftp://ftp.cpc.ncep.noaa.gov/wd52dg/data/indices/tele_index.nh
゚ェセケ
1950 1月～2011 9月 ゾヤカネェクョル タヴシをカヌヴ tc.dat 保
バッジヴ カヌヴ い
端曒 カブルチを実行
ょ awk ‘{print $5}’ tc.dat > wp.dat
awk 列 タヴシを り出 カブルチ こ 場合 tc.dat 5列目 タヴシを り出 ょ
ょ タヴシをwp.dat いうネ゙゜ャ 保 る
-99.90 残 いる ある 不要 geditを開い 操作 消去 る
ょ 検索 S → 置換 R → 検索 る文 列 S 欄 ”-99.90” 入力
ょ 全 置換 A を選択 る い 消去 れる
相関係数 回帰直線用
気象庁｜日曓 月 均気温偏差
http://www.data.kishou.go.jp/climate/cpdinfo/temp/list/mon_jpn.html
ょ ゚ェセケ 1950-2010 タヴシをカヌヴ monthly_anomaly_japan.dat
いうネ゙゜ャ 保 バッジヴ け い
端曒 カブルチを実行
ょ awk ‘{print $2” ” $3” ” $4” ” $5” ” $6” ” $7” ” $8” ” $9” ”
$10” ” $11” ” $12” ” $13}’ monthly_anmaly_japan.dat > monthly_japan.dat ハュエメヘを実行
change.f90
! 列 タ 変換 グ
program change
implicit none
integer :: i, j
integer,parameter :: nn=61, nnn=12, n=nn*nnn
real,dimension(nn,nnn) :: xx
!nn: 縦 数
!nnn: 横 数
open(10,file='monthly_japan.dat',status='old')
do i=1,nn
read(10,*) (xx(i,j),j=1,nnn)
end do
open(20,file='T_mon_jap.dat',status='unknown')
do i=1,nn
do j=1,nnn
write(20,100) xx(i,j)
end do
end do
100 format(f5.2)
stop
end program change
tc.dat う 1950 1 月～2010 12 月 ゾヤカネェクョル タヴシをカヌヴ tc2.dat 保
端曒 カブルチを実行
ょ awk ‘{print $5}’ tc2.dat > wp1950_2010.dat ハュエメヘを実行
wp_jt.f90
program gmt_dataset
implicit none
integer :: i
integer,parameter :: n=732
real,dimension(n) :: x,y
open(10,file='wp1950_2010.dat',status='old')
do i=1,n
read(10,*) x(i)
end do
close(10)
open(20,file='T_mon_jap.dat',status='old')
do i=1,n
read(20,*) y(i)
end do
close(20)
open(30,file='wp_T_1950_2010.dat',status='unknown')
do i=1,n
write(30,100) x(i),y(i)
end do
close(30)
100 format(f5.2,1x,f5.2)
stop
end program gmt_dataset
相関係数 sokan.f90
program correlation_coefficient implicit none
integer :: i
integer,parameter :: n=732
real :: sx,sy,rx,ry,ax,ay,ssx,ssy,ssxy,cc real,dimension(n) :: x,y
! れ れ ドメベヴシ 明
! n:タヴシ数
! sx,sy:x,y れ れ 和
! ax,ay:x,y れ れ 均
! rx,ry:x,y れ れ 偏差
! ssx,ssy,ssxy:x,y,xy 散
! cc:相関係数
! 変える ころ
! n,ネ゙゜ャ
open(10,file='wp1950_2010.dat') do i=1,n
read(10,*) x(i) end do
open(20,file='T_mon_jap.dat') do i=1,n
read(20,*) y(i) end do
! x 散を計算 sx=0.0
do i=1,n sx=sx+x(i) end do
ax=sx/n
ssx=0.0 do i=1,n rx=x(i)-ax ssx= ssx+rx*rx end do
! y 散を計算 sy=0.0
do i=1,n sy=sy+y(i)
end do
ay=sy/n
ssy=0.0 do i=1,n ry=y(i)-ay ssy= ssy+ry*ry end do
! xy 散を計算 ssxy=0.0
do i=1,n rx=x(i)-ax ry=y(i)-ay
ssxy= ssxy+rx*ry end do
! 相関係数 計算
cc=(ssxy/n)/(sqrt(ssx/n)*sqrt(ssy/n)) write(*,*) "相関係数"
write(*,100) cc
100 format(f5.2)
stop
end program correlation_coefficient
回帰直線 kaiki.f90
容：暷小二乗法 回帰直線を求 るハュエメヘ
傾 ：
,
片：
program regression_line implicit none
integer :: i
real :: sx,sy,sxy,sx2,a,b integer,parameter :: n=732 real,dimension(n) :: x,y
! 回帰直線 傾 a 片 b を計算
! a = n∑(xi*yi)-∑xi∑yi / n∑xi^2-(∑xi)^2
! b = ∑xi^2∑yi-∑(xi*yi)∑xi / n∑xi^2-(∑xi)^2
! タヴシ 込
open(10,file='wp1950_2010.dat',status='old') do i=1,n
read(10,*) x(i) end do
open(20,file='T_mon_jap.dat',status='old') do i=1,n
read(20,*) y(i) end do
! 和 計算 sx=0.0
do i=1,n sx=sx+x(i) end do
sy=0.0 do i=1,n sy=sy+y(i) end do
sxy=0.0 do i=1,n
sxy=sxy+x(i)*y(i) end do
! 2乗和 計算 sx2=0.0 do i=1,n
sx2=sx2+x(i)**2 end do
! 傾 a 片 b を求 る
a=(n*sxy-sx*sy)/(n*sx2-(sx)**2) b=(sx2*sy-sxy*sx)/(n*sx2-(sx)**2)
! a,bを画面 出力 write(*,*) "回帰曲線" write(*,*) "y=",a,"x+",b
!GMT用 タヴシを作成
open(30,file='gmt_kaiki.dat',status='unknown') write(30,*) "-4", a*(-4)+b
write(30,*) "4",a*4+b close(30)
stop
end program regression_line
ナケダエメヘ histo.f90
program histogram
! ドメベヴシ文 タヴシ数 n ナケダエメヘを書く際 横軸 間隔 d を指定
integer :: i,j,mini,x real :: dmin,dmax
integer,parameter :: n=741 real,parameter :: d=1.0 real,dimension(n) :: dat integer,dimension(8) :: rf
! タヴシ 込
open(10,file='wp.dat',status='old') do i=1,n
read(10,*) dat(i) end do
! 並び暶え
! 値 小 い ら順 書く
dmin=0.0 do j=1,n-1 dmin = dat(j) mini = j do i=1+j,n
if (dat(i)<dmin)then dmin=dat(i)
mini=i end if end do
if(mini /= j)then dat(mini) = dat(j) dat(j) = dmin end if
ハュエメヘ ゜ベヴグ
暷小値を け 何番目 タヴシ をmini
タヴシを dat(j) 入れる 値を dmin る
mini番目 外 タヴシ 中 ら暷小値dminを探 dat(j) 入れる
をくり え
end do
! 結果 出力
open(20,file='wp2.dat',status='unknown') do j=1,n
write(20,100) dat(j) end do
! 暷小値を書く do j=1,1 dmin=dat(j) write(*,110) dmin enddo
! 暷大値を書く do j=n,n dmax=dat(j) write(*,120) dmax enddo
! ネァヴブッダ文 100 format(f5.2)
110 format(' ','min=',6f6.2) 120 format(' ','max=',6f6.2)
! ナケダエメヘ 作成
x=(int(dmax)+d)-(int(dmin)-d) rf(:) = 0
do j=1,n
if(dat(j) < int(dmin))then rf(1)=rf(1)+1
end if
A
do i=0,x-3
if (int(dmin) + i*d <= dat(j) .and. dat(j) < int(dmin) + (i+1) * d ) then rf(i+2) = rf(i+2) + 1
end if end do
if (int(dmin)+(x-2)*d < dat(j)) then rf(x) = rf(x) + 1
end if end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! rf:相対度数(relative frequency) open(30,file='rf.dat')
do j = 1,x
write(30,*) int(dmin)-int(d)+int(d*(j-1)),rf(j) end do
stop
end program histogram
ナケダエメヘ 値を求 る部 ゜ベヴグ
求 暷小値 暷大値 整数部 をintを使 求 る
int(dmin)-d 横軸 暷小値 int(dmax)+d 暷大値を求 る
A ころ 整数部 int(dmin)-dを る場合をィゞルダ る C 部
整数部 int(dmax)+dを る場合を れ 外 ころをB 部 れ れ
ィゞルダ る
ょ
B
C
C ィゞルダ る
A ィゞルダ る
れ れB ィゞルダ る
れ れ 横軸 値を求 る
GMT
概要
GMT General Mapping Tools 略
エメネや地図ン 布図 作成ソネダ
地図投影法 豊富 細 い海岸線タヴシ 暷初 ら利用可
LINUX カブルチを用い 重 書 作図 いく
GrADSより手間 る反面 見栄え よい図 作成 る
地点タヴシ 格子点タヴシ く 比較的容易 図 描ける
GrADS 扱えるタヴシ る いン面倒 場合 GMT 楽
実行 方
GMT カブルチをクゟャ 書い 実行 る
具体的 filename.csh filename.sh いうネ゙゜ャを く こ カブル
チを書
csh filename.csh sh filename.sh
実行 る
実際 例
地図 描画
map1.csh
pscoast -R -180/180/-80/80 -JM15c -Ba30g15 -G150 -S255 > map1.eps
＜ 明＞
pscoast: 海岸線 を地 ッ コマ
-R-180/180/-80/80: ッ 範 をコ ショ 場合 経度 -180
180 経度 -80 80 範 を ッ いう意味
-JM15c: ッ 形式 大 をコ ショ J 後
”M”or”m” を使用 指定 15c 大 場合
横幅 15cm いう意味
-Ba30g15: ベ や目盛をコ ショ ”a30” 30 度 ベ を付
”g15” 15 グ ッ を ッ いう意味
-G150: 陸地 色をコ ショ RGB 指定 例え
”-G/255/0/0” 赤 ッ 色 濃 0 ～ 255 範 選択
-S255: 海 色をコ ショ 使用方法 ”G” 同
> map.eps:pscoast コマ 出力を map.eps
＜そ 他 使う ショ ＞
-D: 海岸線 密度をコ ショ (-Df: full; -Dh: high; -Di:
intermediate; -Dl : low)
-W(numbers): 線 太 をコ ショ
-N(numbers): 境界線を 描 コ ショ
-N1=National boundaries; -N2=State boundaries within the Americans; -N3=Marine
boundaries; -Na=All boundaries(1 ～ 3)
-P: を縦書
地 を入 う
map2.csh
pscoast -R -180/180/-80/80 -JM15c -Ba30g15 -G150 -S255 -K > map2.eps
pstext -R -JM -O << EOF >> map2.eps
135 35 18 0 1 1 Japan
EOF
-K: 作成を 使用 ショ
-O: 上書 使用 ショ
-R: 後 何 い場合 前 設定を引 い い
<< EOF: 入力 場合 使用 EOF 終了
>>: ァイ 後 出力を 使用
135 35 18 0 1 1 Japan ： x 位置 y 位置 サイ
角度 を配置 位置 を意味
シ ボ を ッ う
map3.csh
pscoast -R-180/180/-80/80 -JM15c -Ba30g15 -G150 -S255 -K > map3.eps
pstext -R -JM -O -K << EOF >> map3.eps
135 35 18 0 1 1 Japan
EOF
psxy -R -JM -Sa0.1c -O -G255/255/0 -W1 << EOF >> map3.eps
135 35
EOF
psxy: シ ボ や線を ッ 使用 コマ
-S: シ ボ タイ 大 をコ
先ほ 作成 渦度 タヴシを描画 よう！
uzudo.csh
pscoast -JQ160/20.0 -R80/280/10/80 -G100 -Ba10g10WeSn -X3.0 -Y2.0
-K > uzudo.eps
surface rvort.dat -Gvorticity.grd -R -I2.5/2.5 -T0.25
grdcontour vorticity.grd -R -JQ -C1 -A1 -W0/0/0 -O >> uzudo.eps
エメネ 描画
先ほ 作 タヴシをエメネ よう！
regression.csh ょ ょ ょ 回帰直線 相関係数を書く
echo "6.0 12.5 14 0 5 2 R=0.30, Y=0.27X-0.23" >! sokan.ymd pstext sokan.ymd -R0/8.5/0/11.0 -Jx1 -P -N -K > WPI_T.eps psxy wp_T_1950_2010.dat -Jx1.0 -R-4/4/-4/4 -Sc0.1
-Ba2.0f1.0g1.0:"WPI":/a2.0f1.0g1.0:"Temperature Anomaly (Japan)" :WSne:."Correlation Chart": -P -V -X2 -Y4 -O -K >> WPI_T.eps psxy gmt_kaiki.dat -Jx -R -W3/255/0/0 -O >> WPI_T.eps
明：
”R=0.30, Y=0.27X-0.23” 書 方 指定
”R=0.30, Y=0.27X-0.23” 書 く カ ブ ル チ ” wp_T_1950_2010.dat ” ”ょ wp1950_2010.dat ” “ T_mon_jap.dat ”を2列 書い タヴシ
相関図を描く チッダを打 カブルチ
回帰直線を描くカブルチ
histo.csh ょ ょ ナケダエメヘを作成
psxy rf.dat –JX16/16 –X3 –Y7 –R-5/4/0/300 –I1
–Ba1f1:WPI:/a50f25:Frequency:neWS –Sb2 –G0/0/255 –P –V > rf.eps
役 立 ハュエメヘ
