ごくらく DCRTM

はじめに

この文書は DCRTM を用いて手軽に実験を行うためのチュートリアルです.

DCRTM のビルド

DCRTM インストールガイド を参考に, DCRTM のビルドを行ってください. 「ビルドの手引き」の「makefile の作成, 編集, コンパイル」まで行ってください.

コンパイルが成功すると, dcrtm-バージョン ディレクトリ内に実行ファイルが作成されます. またこのディレクトリ内には, NAMELIST ファイル (dcrtm.nml) と, 簡単な実験を行うための鉛直プロファイルデータ(sample/USstandard_MLS_L032M03.txt)が用意されています. 以下では, この鉛直プロファイルを使って実験の手順を説明します.

実験の実行

DCRTM を用いた放射計算は, 以下の3つのステップで行います.

• 鉛直プロファイル(圧力, 温度, 組成) の準備
• opacity の計算
• フラックスの計算

鉛直プロファイル(圧力, 温度, 組成) の準備

DCRTM では入力データとして, 大気の鉛直プロファイル (圧力, 温度, 組成) を必要とします. #また, 計算しようとする大気に含まれる分子の分子量もあわせて用意します.

• 気圧 [Pa]
• 気温 [K]
• 組成(体積混合比) [1]

ここでは, 必要なテキストフォーマットになった地球中緯度夏の鉛直プロファイル (Ellingson et al., 1991) をもとに, DCRTM での計算に必要な netcdf ファイルに変換します. NAMELIST ファイル (dcrtm.nml) で以下のようにファイル名を指定します.

&initial_nml
InFilePTX       = './sample/USstandard_MLS_L032M03.txt',
                !
InFileMolWt     = './src/data/MolWt.txt',
                !
OutFileName     = './sample/USstandard_MLS_L032M03.nc'
/

上記 NAMELIST ファイルの設定ができたら, プログラムを実行します.

$ ./dcrtm_initialdata

これで, 光学特性の計算をするための鉛直プロファイルデータが作成されます.

opacity の計算

ここでは, 地球中緯度夏の長波放射の計算を行うことにします. 設定は NAMELIST ファイル (dcrtm.nml) で行います. opacity 計算で使用する NAMELIST 群 は optdep_nml です. NAMELIST を以下のように設定してください.

&optdep_nml
InFileName     = './sample/USstandard_MLS_L032M03.nc',
               !
DataBase       = 'HITRAN2012',
               ! DataBase: HITRAN1996, HITRAN2008, HITRAN2012, HITEMP2010 (,TestData)
PATHQ          = './src/optdep/parsum.dat.nc'
               !
MolWt0         = 28.96e-3,
               ! molecular weight of non-absorption gas [kg mol-1]
MinWN          = 0.0,
               ! mininum wavenumber [m-1]
MaxWN          = 300000.0,
               ! maximum wavenumber [m-1]
ResWN          = 1.0,
               ! wavenumber resolution [m-1]
Flag_ISOTOPE   = 0,
               ! 0: Earth abandance, 1: read user defined abandance
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_CUTOFF    = 0,
               ! 0: 25cm-1 all molecules, 1: read user defined cut off value (cutoff_nml)
Flag_CONT      = 0,
               ! 0: no continuum, 1: continuum (MT_CKD model), 2: continuum (Pollack+1993)
Flag_RAYLEIGH  = 1,
               ! 0: no Rayleigh scattering, 1: Rayleigh scattering
Flag_OUTPUT    = 1,
               ! 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)
                 ! if InAngle is negative, solar radiation is not calculated.
Flag_FluxModel   = 1,
                 ! 1: Toon et al. (1989)
                 !   for solar radiation: Two-Stream Equation (Eddington App.)
                 !   for planetary radiation: Two-Stream Source Function Method (Hemispheric Mean)
                 ! 2: Toon et al. (1989), Two-Stream Equation (for solar and planetary radiation)'
                 !   for solar radiation: Two-Stream Equation (Eddington App.)'
                 !   for planetary radiation: Two-Stream Equation (Hemispheric Mean)'
                 ! 3: Toon et al. (1989) + Source Function Iteration Method
                 !   for solar radiation: Two-Stream Equation (Eddington App.)
                 !   for planetary radiation: Two-Stream Source Function Iteration Method (Hemispheric Mean)
Flag_Emit        = 1,
                 ! 0: No emitting atmosphere, 1: emitting atmosphere
Flag_Particle    = 1,
                 ! 0: No particle atmosphere, 1: particles are contained
InFileParticle   = './sample/ParticleSample.nc'
                 !
/

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 では, 簡単な解析プログラムを用意しています. 用意されているのは, 鉛直フラックスを求めるためのフラックスの波数積分計算プログラム (dcrtm_vflux), 波数ごとの加熱率計算プログラム (dcrtm_heatrate), 大気上端から出て行くフラックスの輝度温度計算プログラム (dcrtm_brightt) です.

鉛直フラックスの計算をしてみます. この解析で使用する NAMELIST 群 は vflux_nml です. NAMELIST を以下のように設定してください.

&vflux_nml
InFileName     = 'Flux_USstandard_MLS_L032M03.nc',
               !
OutFileName    = 'vFlux_USstandard_MLS_L032M03.nc',
               !
/

NAMELIST が設定できたら, プログラムを実行します.

$ ./dcrtm_vflux

$ 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