ごくらく DCRTM
はじめに
この文書は dcrtm を用いて手軽に実験を行うためのチュートリアルです.
DCRTM のビルド
DCRTM インストールガイド を参考に, dcrtm のビルドを行ってください. 「ビルドの手引き」の「makefile の作成, 編集, コンパイル」まで行ってください.
コンパイルが成功すると, dcrtm-バージョン ディレクトリ内に実行ファイルが作成されます. またこのディレクトリ内には, NAMELIST ファイル (dcrtm.nml) と, 簡単な実験を行うための鉛直プロファイルデータ(sample/USstandard_MLS_L032M03.nc)が用意されています. 以下では, この鉛直プロファイルを使って実験の手順を説明します.
実験の実行
dcrtm を用いた放射計算は, 以下の3つのステップで行います.
- 鉛直プロファイル(圧力, 温度, 組成) の準備
- opacity の計算
- フラックスの計算
鉛直プロファイル(圧力, 温度, 組成) の準備
dcrtm では入力データとして, 大気の鉛直プロファイル (圧力, 温度, 組成) を必要とします. また, 計算しようとする大気に含まれる分子の分子量もあわせて用意します.
- 気圧 [Pa]
- 気温 [K]
- 組成(体積混合比) [1]
- 分子量 [kg mol-1]
ここでは, 地球中緯度夏の鉛直プロファイル (Ellingson et al., 1991) を使用します.
opacity の計算
ここでは, 地球中緯度夏の長波放射の計算を行うことにします. 設定は NAMELIST ファイル (dcrtm.nml) で行います. opacity 計算で使用する NAMELIST 群 は optdep_nml です. NAMELIST を以下のように設定してください.
&optdep_nml
InFileName = './sample/USstandard_MLS_L032M03.nc',
!
DataBase = 'HITRAN2012',
! DataBase: HITRAN1996, HITRAN2008, HITRAN2012, HITEMP2010
PATHQ = './src/optdep/parsum.dat'
!
MinWN = 0.0,
! mininum wavenumber [m-1]
MaxWN = 300000.0,
! maximum wavenumber [m-1]
ResWN = 1.0,
! wavenumber resolution [m-1]
Flag_LINE = 1,
! 0: no line calculation, 1: line calculation
Flag_LINESHAPE = 1
! 1: voigt (all), 2: voigt & CO2 sub-Lorentzian
! 3: Lorentz (all), 4: Lorentz & CO2 sub-Lorentzian
Flag_CONT = 0,
! 0: no continuum, 1: continuum (MT_CKD model)
MolWt0 = 28.96e-3,
! molecular weight of non-absorption gas [kg mol-1]
Flag_OUTPUT = 3,
! 1: optical depth only, 2, coefficient only, 3. optical depth and coefficient
OutFileOptDep = 'OptDep_USstandard_MLS_L032M03.nc',
!
OutFileOptProp = 'OptProp_USstandard_MLS_L032M03.nc',
!
/
DataBase については, 実験に使用する HITRAN バージョンを選んでください. 今回の設定は, 計算波数 0 - 300000 [m-1] (0 - 3000 [cm-1]), 波数解像度 1.0 [m-1] (0.01 [cm-1]) です. (公開版では連続吸収の効果は計算されません. Flag_CONT = 0 として実行してください.) なお, Flag_OUTPUT を, 2 とすると, 入力データで与えた温度, 圧力, 組成での吸収係数, 散乱係数を出力します. (3 とすると両方出力されます)
!! 注意 !! この計算は, それなりに時間がかかります. (計算環境にもよりますが, 1時間半くらい??)
初めて, dcrtm を実行する方は, 波数領域を小さくして実行してください. (例えば MaxWN = 10.0)
上記 NAMELIST ファイルの設定ができたら, プログラムを実行します.
$ ./dcrtm_optdep
計算終了までしばらく時間がかかることがあります. 計算が終了すると, opacity の計算結果が netcdf ファイルで出力されます.
$ ls OptDep_USstandard_MLS_L032M03.nc OptProp_USstandard_MLS_L032M03.nc
出力された netcdf データのヘッダの部分を見ると以下のようになっています.
$ ncdump -h OptDep_USstandard_MLS_L032M03.nc
netcdf OptDep_USstandard_MLS_L032M03 {
dimensions:
r_Press = 33 ;
w_WaveNum = 300001 ;
z_Press = 32 ;
variables:
float r_Press(r_Press) ;
r_Press:long_name = "boundary pressure" ;
r_Press:units = "Pa" ;
float w_WaveNum(w_WaveNum) ;
w_WaveNum:long_name = "wavenumber" ;
w_WaveNum:units = "m-1" ;
float z_Press(z_Press) ;
z_Press:long_name = "mid-layer pressure" ;
z_Press:units = "Pa" ;
double r_Temp(r_Press) ;
r_Temp:long_name = "temperature" ;
r_Temp:units = "K" ;
double rw_OptDep(w_WaveNum, r_Press) ;
rw_OptDep:long_name = "optical depth from TOA" ;
rw_OptDep:units = "1" ;
double zw_SingleScatA(w_WaveNum, z_Press) ;
zw_SingleScatA:long_name = "single scattering albedo" ;
zw_SingleScatA:units = "1" ;
rw_OptDep は, 大気上端からある高度までの波数ごとの光学的厚さです. zw_SingleScatA は, それぞれの大気層での 1 次散乱アルベドです.
$ ncdump -h OptProp_USstandard_MLS_L032M03.nc
netcdf OptProp_USstandard_MLS_L032M03 {
dimensions:
r_Press = 33 ;
m_MolNum = 4 ;
w_WaveNum = 300001 ;
variables:
float r_Press(r_Press) ;
r_Press:long_name = "boundary pressure" ;
r_Press:units = "Pa" ;
float m_MolNum(m_MolNum) ;
m_MolNum:long_name = "molecular number" ;
m_MolNum:units = "1" ;
float w_WaveNum(w_WaveNum) ;
w_WaveNum:long_name = "wavenumber" ;
w_WaveNum:units = "m-1" ;
double m_MolWt(m_MolNum) ;
m_MolWt:long_name = "molecular weight" ;
m_MolWt:units = "kg mol-1" ;
double r_Temp(r_Press) ;
r_Temp:long_name = "temperature" ;
r_Temp:units = "K" ;
double rm_MixRatio(m_MolNum, r_Press) ;
rm_MixRatio:long_name = "volume mixing ratio" ;
rm_MixRatio:units = "1" ;
double rmw_LAbsorpK(w_WaveNum, m_MolNum, r_Press) ;
rmw_LAbsorpK:long_name = "line absorption coefficient" ;
rmw_LAbsorpK:units = "m2 kg-1" ;
double rmw_CAbsorpK(w_WaveNum, m_MolNum, r_Press) ;
rmw_CAbsorpK:long_name = "continuum absorption coefficient" ;
rmw_CAbsorpK:units = "m2 kg-1" ;
double rmw_TAbsorpK(w_WaveNum, m_MolNum, r_Press) ;
rmw_TAbsorpK:long_name = "line & continuum absorption coefficient" ;
rmw_TAbsorpK:units = "m2 kg-1" ;
double rmw_RayScatK(w_WaveNum, m_MolNum, r_Press) ;
rmw_RayScatK:long_name = "rayleigh scattering coefficient" ;
rmw_RayScatK:units = "m2 kg-1" ;
フラックスの計算
フラックスの計算は, 上記で計算した opacity の計算結果をもとに行います. フラックス計算で使用する NAMELIST 群 は flux_nml です. NAMELIST を以下のように設定してください.
&flux_nml
Flag_InputFormat = 1
! 1: optical depth and single scattering albedo
! 2: absorption, scattering coefficient
InFileName = 'OptDep_USstandard_MLS_L032M03.nc'
!
OutFileFlux = 'Flux_USstandard_MLS_L032M03.nc'
!
SurfAlbedo = 0.0
! surface albedo
SolTemp = 5800.0
! temperature of central star [K]
SolConst = 1366.0
! solar constant [W m-2]
InAngle = 2.0
! sec (angle of incidence)
/
Flag_InputFormat で, opacity データの選択を行います. opacity 計算の際, 光学的厚さ, 一次散乱アルベドを出力した場合は 1, 吸収散乱係数を出力した際は 2 を選択してください. フラックスの計算では, 地表面アルベドと入射太陽光の設定を行います. 今回の設定は, 地表面アルベド = 0, 太陽光は 5800 [K] の黒体スペクトル(太陽定数 1366[W m-2]) で入射角度が 60 度 (sec(入射角)= 2) としています.
上記 NAMELIST ファイルの設定ができたら, プログラムを実行します.
$ ./dcrtm_flux
計算終了までしばらく時間がかかることがあります. 計算が終了すると, フラックスの計算結果が netcdf ファイルで出力されます.
$ ls Flux_USstandard_MLS_L032M03.nc
出力された netcdf データのヘッダの部分を見ると以下のようになっています.
$ ncdump -h Flux_USstandard_MLS_L032M03.nc
netcdf Flux_USstandard_MLS_L032M03 {
dimensions:
r_Press = 33 ;
w_WaveNum = 300001 ;
variables:
float r_Press(r_Press) ;
r_Press:long_name = "pressure" ;
r_Press:units = "Pa" ;
float w_WaveNum(w_WaveNum) ;
w_WaveNum:long_name = "wavenumber" ;
w_WaveNum:units = "m-1" ;
double r_Temp(r_Press) ;
r_Temp:long_name = "temperature" ;
r_Temp:units = "K" ;
double rw_RadFluxUp(w_WaveNum, r_Press) ;
rw_RadFluxUp:long_name = "upward flux(planet)" ;
rw_RadFluxUp:units = "W m-2 m" ;
double rw_RadFluxDn(w_WaveNum, r_Press) ;
rw_RadFluxDn:long_name = "downward flux(planet)" ;
rw_RadFluxDn:units = "W m-2 m" ;
double rw_SolFluxUp(w_WaveNum, r_Press) ;
rw_SolFluxUp:long_name = "upward flux(solar)" ;
rw_SolFluxUp:units = "W m-2 m" ;
double rw_SolFluxDn(w_WaveNum, r_Press) ;
rw_SolFluxDn:long_name = "downward flux(solar)" ;
rw_SolFluxDn:units = "W m-2 m" ;
rw_RadFluxUp, rw_RadFluxDn は惑星大気の熱放射に由来する上向き, 下向きフラックスです. rw_SolFluxUp, rw_SolFluxDn は入射太陽光に由来する上向き, 下向きフラックスです.
結果の可視化
DCRTM では入出力するファイルとして Gtool4 NetCDF 規約 に基づいた NetCDF データを扱います.
数値実験の結果を解析・可視化するためには, NetCDF データを取り扱うことのできる解析・可視化ツールが必要です. ここでは, 電脳 Ruby プロジェクト から提供される Gphys を使った可視化の例を紹介します.
可視化ツールのインストールについては, 電脳Ruby謹製品 インストールガイド を参照してください.
フラックススペクトルの可視化
惑星大気の熱放射に由来する上向き, 下向きフラックスについて, スペクトルを可視化します.
大気上端での上向きフラックス (惑星大気から宇宙空間へ射出される放射フラックス) は, Gphys のコマンドを以下のように使って描画することができます.
$ gpview Flux_USstandard_MLS_L032M03.nc@rw_RadFluxUp,r_Press=0 --range 0:0.005
rw_RadFluxUp は上向きフラックス, r_Press=0 は, 大気上端 (圧力 = 0 [Pa]) での値を表しています. --range 0:0.005 は, フラックスの値の範囲を指定しています.
地表面 (r_Press=101300 [Pa]) での, 上向き放射フラックスは, 下記のようにして描画することができます.
$ gpview Flux_USstandard_MLS_L032M03.nc@rw_RadFluxUp,r_Press=101300 --range 0:0.005
地表面温度でのプランク放射となっています.
地表面 (r_Press=101300 [Pa]) での, 下向き放射フラックスは, 下記のようにして描画することができます.
$ gpview Flux_USstandard_MLS_L032M03.nc@rw_RadFluxDn,r_Press=101300 --range 0:0.005
鉛直フラックスの可視化
波数方向に積分した鉛直フラックスを計算するために, DCRTM では, 解析プログラムを用意しています.
プログラムは, src/analysis/sum_verticalprof.f90 です. 以下のようにしてコンパイルします.
$ gt5frt src/analysis/sum_verticalprof.f90
この解析で使用する NAMELIST 群 は vflux_nml です. NAMELIST を以下のように設定してください.
&vflux_nml
InFileName = 'Flux_USstandard_MLS_L032M03.nc',
!
OutFileName = 'vFlux_USstandard_MLS_L032M03.nc',
!
/
NAMELIST が設定できたら, プログラムを実行します.
$ ./a.out (実行ファイルが a.out の場合) $ ls vFlux_USstandard_MLS_L032M03.nc
出力された netcdf データのヘッダの部分を見ると以下のようになっています.
$ ncdump -h Flux_USstandard_MLS_L032M03.nc
netcdf vFlux_USstandard_MLS_L032M03 {
dimensions:
r_Press = 33 ;
variables:
float r_Press(r_Press) ;
r_Press:long_name = "pressure" ;
r_Press:units = "Pa" ;
double r_RadFluxUp(r_Press) ;
r_RadFluxUp:long_name = "upward flux(planet)" ;
r_RadFluxUp:units = "W m-2" ;
double r_RadFluxDn(r_Press) ;
r_RadFluxDn:long_name = "upward flux(planet)" ;
r_RadFluxDn:units = "W m-2" ;
double r_SolFluxUp(r_Press) ;
r_SolFluxUp:long_name = "upward flux(solar)" ;
r_SolFluxUp:units = "W m-2" ;
double r_SolFluxDn(r_Press) ;
r_SolFluxDn:long_name = "upward flux(solar)" ;
r_SolFluxDn:units = "W m-2" ;
double r_RadNetFlux(r_Press) ;
r_RadNetFlux:long_name = "net upward flux(planet)" ;
r_RadNetFlux:units = "W m-2" ;
double r_SolNetFlux(r_Press) ;
r_SolNetFlux:long_name = "net upward flux(solar)" ;
r_SolNetFlux:units = "W m-2" ;
r_RadNetFlux, r_SolNetFlux はそれぞれ, 惑星大気の熱放射, 太陽光に由来する放射の正味上向きフラックスを表しています.
惑星大気の熱放射について, 正味上向き, 上向き, 下向き放射を可視化します.
$ gpview --overplot=3 vFlux_USstandard_MLS_L032M03.nc@r_RadNetFlux vFlux_USstandard_MLS_L032M03.nc@r_RadFluxUp vFlux_USstandard_MLS_L032M03.nc@r_RadFluxDn --ex --itr 2 --range 0:450
