% % 表題 SPMODEL サンプルプログラムドキュメント % gravcurrent1.f90 % % 履歴 2004/05/27 小高 正嗣 % 2004/06/08 小高 正嗣 % 2026/07/08 竹広 真一 for et_module % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \documentclass[a4j,12pt,dvipdfmx]{jarticle} % % プリアンブル % \usepackage{Dennou6} % dennou-sty-6 を使用 \Dtitle[{\small 2 次元ブシネスクモデル:重力流の計算}] % ヘッダ部タイトル, モデル名と同じ { {\large SPMODEL サンプルプログラム} \\ % ここは変更しない 水路領域 \\ 2 次元ブシネスク方程式モデル:\\ % モデル名 重力流の計算 \\ {\large gravcurrent1.f90} % プログラム名 } % \Dauthor[小高 正嗣, 竹広 真一]{小高 正嗣, 竹広 真一} % 著者 % \Ddate{2004 年 6 月 8 日} % 日付 \Ddate{2026 年 7 月 8 日} % 日付 \Dnoparindent % 段落の字下げをしない % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % 本文開始 % % 基本的な章節構成は以下の通り. % % 1. はじめに % 2. 支配方程式系 % 3. 離散化 % 4. 利用モジュールとその他の設定 % 5. 数値実験 % 6. 参考文献 % -. 謝辞 % % 必要に応じて構成を変更すること. % \begin{document} \maketitle % タイトルの作成 \tableofcontents % 目次の作成 \pagebreak %---------------------------------------------------------------------- \Dparskip %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{概要} SPMODEL サンプルプログラム『gravcurrent1.f90』に用いられている基礎方程 式と境界条件, および, このプログラムを用いた数値実験の方法について解説 する. 基礎方程式は 2 次元のブシネスク方程式系である. 計算はスペクトル法を用いて行い, 展開関数は水平方向にはフーリエ級数, 鉛直方向にはチェビシェフ多項式を用いる. スペクトル変換と逆変換および微分演算には, SPMODEL ライブラリ(spml)を用いている. 数値実験では重力流の計算を行った.a {\bf プログラム名} \\ % プログラム名 {\footnotesize gravcurrent1.f90 } {\bf プログラム取得元}\\ % プログラム取得先 {\footnotesize http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash 2d-channel-et\slash boussinesq\slash gravcurrent\slash SIGEN.htm } {\bf SPMODEL サンプルプログラム目次}\\ % サンプルプログラム目次(変更しない) {\footnotesize http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash sample.htm } {\bf SPMODEL の使い方}\\ % SPMODEL の使い方(変更しない) {\footnotesize http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\/ } %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{支配方程式系} ここでは支配方程式系と境界条件を記す. \subsection{支配方程式系} 支配方程式系は 2 次元のブシネスク方程式系である. \begin{eqnarray} && \DP{u}{t}+u\DP{u}{x}+v\DP{u}{y} = -\frac{1}{\rho _{0}}\DP{p}{x} + \nu \Dlapla u, \\ && \DP{v}{t}+u\DP{v}{x}+v\DP{v}{y} = -\frac{1}{\rho _{0}}\DP{p}{y} + \frac{\rho}{\rho _{0}}g + \nu \Dlapla v, \\ && \DP{u}{x}+\DP{v}{y} = 0,\\ && \DP{\rho}{t}+u\DP{\rho}{x}+v\DP{\rho}{y} = \kappa \Dlapla \rho. \\ \end{eqnarray} 各記号の定義は\Dtabref{変数, 物理定数の定義}に表す. 渦度 $\zeta$ と流線関数 $\psi$ \begin{eqnarray} &&\zeta =\DP{v}{x}-\DP{u}{y}, \\ &&v=\DP{\psi}{x}, \quad u=-\DP{\psi}{y}, \end{eqnarray} を導入し, 渦度と密度の予報方程式として支配方程式系を表すと以下のようになる. \begin{eqnarray} && \DP{\zeta }{t} + J(\psi, \zeta) + \frac{g}{\rho _{0}}\DP{\rho }{x} = \nu \Dlapla \zeta , \Deqlab{渦度方程式}\\ && \DP{\rho}{t} + J(\psi, \rho) = \kappa \Dlapla \rho. \Deqlab{密度の予報方程式} \end{eqnarray} ここで $J(A, B)$ はヤコビアン \[ J(A,B) = \DP{A}{x}\DP{B}{y} - \DP{A}{y}\DP{B}{x} \] である. \begin{table}[h] \begin{center} \begin{tabular}{c l l } \hline 記号 &\qquad\qquad& 変数/物理定数 \\ \hline $x$ && 水平座標 \\ $z$ && 鉛直座標 \\ $t$ && 時間 \\ $u$ && $x$ 方向速度 \\ $v$ && $y$ 方向速度 \\ $\rho_{0}$ && 平均密度 \\ $\rho $ && 密度 \\ $\psi$ && 流線関数 \\ $\zeta$ && 渦度 \\ $g$ && 重力加速度 \\ $\nu$ && 動粘性係数 \\ $\kappa$ && 拡散係数 \\ \hline \end{tabular} \caption{変数, 物理定数の定義} \Dtablab{変数, 物理定数の定義} \end{center} \end{table} \subsection{境界条件} 境界条件は水平に周期境界条件, 鉛直方向には $z=0, y_{m}$ に置いた剛体壁 面で $v=0$, 応力なし, 熱フラックスなしとする. すなわち水平境界条件は 計算領域をそれぞれ $x_{m}$ とすると \begin{equation} \zeta(x+x_{m},z) = \zeta(x,z) \end{equation} などと表される. 鉛直境界条件は $z=0, y_{m}$ において, \begin{eqnarray} &&\psi = \mbox{Const.}, \\ &&\zeta = 0, \\ &&\DP{\rho}{y} = 0 \end{eqnarray} である. 境界で与える $\psi $ の値は簡単のため $0$ とした. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % \newpage % 改ページ % \section{離散化} % この節では方程式系の空間離散化および使用した時間積分法について説明し, % プログラム内で実際に用いられている方程式を記述する. % \subsection{空間離散化} % 支配方程式\Deqref{渦度方程式}, \Deqref{密度の予報方程式}の離散表現は以 % 下のようになる. % \begin{eqnarray} % && \DP{\zeta _{i,j}}{t} + % \DP{\psi_{i,j}}{x}\DP{\zeta _{i,j}}{y} - % \DP{\psi_{i,j}}{y}\DP{\zeta _{i,j}}{x} % + \frac{g}{\rho _{0}}\DP{\rho _{i,j}}{x} % = \nu \Dlapla \zeta_{i,j} , \Deqlab{離散化された渦度方程式}\\ % && \DP{\rho_{i,j}}{t} + % \DP{\psi_{i,j}}{x}\DP{\rho _{i,j}}{y} - % \DP{\psi_{i,j}}{y}\DP{\rho _{i,j}}{x} % = \kappa \Dlapla \rho_{i,j}. \Deqlab{離散化された密度の予報方程式} % \end{eqnarray} % ここで添字 $i,j$ は格子点 $(x_{i}, y_{j})$ 上の値であることを表す. % \subsection{空間方向のスペクトル表現} % 空間離散化した支配方程式\Deqref{離散化された渦度方程式}, \Deqref{離散 % 化された密度の予報方程式}をスペクトル法を用いて表現する. スペクトル % 展開は水平方向にフーリエ級数, 鉛直方法には境界条件にあわせてフーリエ正 % 弦級数またはフーリエ余弦級数を用いて行う. 非線形項を扱う場合は, 先に格 % 子点上での非線形項の値を計算し, その値のスペクトルを求める方法(変換法) % を用いる. 浮力項についても同様に扱う. 以下では $k,l$ をそれぞれ$x,z$ % 方向波数, $K,L$ を切断波数, $I,J$ を格子点数とする. % $\zeta_{i,j}, \psi _{i,j}, \rho _{i,j}$ はスペクトル逆変換によって以下 % のように展開される. % \begin{eqnarray} % \zeta_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L} % \exp\left(\frac{2\pi i k x_{i}}{x_{m}}\right) % \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)\hat \zeta_{k,l}, \\ % \psi_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L} % \exp\left(\frac{2\pi i k x_{i}}{x_{m}}\right) % \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)\hat \psi_{k,l}, \\ % \rho_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L} % \exp\left(\frac{2\pi i k x_{i}}{x_{m}}\right) % \cos\left(\frac{\pi l y_{j}}{y_{m}}\right)\hat \rho_{k,l}. % \end{eqnarray} % スペクトル係数 $\hat \zeta_{i,j}, \hat \psi _{i,j}, \hat \rho _{i,j}$ は % スペクトル変換によって以下のように与えられる. % \begin{eqnarray} % \hat \zeta_{k,l} &=& \frac{1}{I}\frac{1}{J} % \sum _{i=0}^{I-1}\sum _{j=0}^{J-1} % \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right) % \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)\zeta_{i,j}, \\ % \hat \psi_{k,l} &=& \frac{1}{I}\frac{1}{J} % \sum _{i=0}^{I-1}\sum _{j=0}^{J-1} % \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right) % \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)\psi_{i,j}, \\ % \hat \rho_{k,l} &=& \frac{1}{I}\frac{1}{J} % \sum _{i=0}^{I-1}\sum _{j=0}^{J-1} % \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right) % \cos\left(\frac{\pi l y_{j}}{y_{m}}\right)\rho_{i,j}. % \end{eqnarray} % 以上を用いると支配方程式\Deqref{離散化された渦度方程式}, \Deqref{離散 % 化された密度の予報方程式}のスペクトル表現は以下のようになる. % \begin{eqnarray} % \DP{\hat \zeta _{k,l}(t)}{t} &=& - % \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1} % \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right) % \sin\left(\frac{\pi l y_{j}}{y_{m}}\right) % \left(\DP{\psi}{x}\right)_{i,j} % \left(\DP{\zeta}{y}\right)_{i,j} \nonumber \\ % && + \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1} % \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right) % \sin\left(\frac{\pi l y_{j}}{y_{m}}\right) % \left(\DP{\psi}{y}\right)_{i,j} % \left(\DP{\zeta}{x}\right)_{i,j} \nonumber \\ % && - \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1} % \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right) % \sin\left(\frac{\pi l y_{j}}{y_{m}}\right) % \frac{g}{\rho _{0}}\left(\DP{\rho}{x}\right)_{i,j} \nonumber \\ % && - \nu \left[ % \left(\frac{2\pi k}{x_{m}}\right)^{2}+ % \left(\frac{\pi l}{y_{m}}\right)^{2}\right]\hat \zeta _{k,l}, % \Deqlab{渦度方程式のスペクトル表現}\\ % \DP{\hat \rho _{k,l}(t)}{t} &=& - % \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1} % \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right) % \cos\left(\frac{\pi l y_{j}}{y_{m}}\right) % \left(\DP{\psi}{x}\right)_{i,j} % \left(\DP{\rho}{y}\right)_{i,j} \nonumber \\ % && + \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1} % \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right) % \cos\left(\frac{\pi l y_{j}}{y_{m}}\right) % \left(\DP{\psi}{y}\right)_{i,j} % \left(\DP{\rho}{x}\right)_{i,j} \nonumber \\ % && - \kappa \left[ % \left(\frac{2\pi k}{x_{m}}\right)^{2}+ % \left(\frac{\pi l}{y_{m}}\right)^{2}\right]\hat \rho _{k,l}. % \Deqlab{密度の予報方程式のスペクトル表現} % \end{eqnarray} % ここで % \begin{eqnarray} % \left(\DP{\zeta}{x}\right)_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L} % \exp\left(\frac{2\pi i k }{x_{m}}\right) % \sin\left(\frac{\pi l }{y_{m}}\right) % \frac{2\pi i k}{x_{m}}\hat \zeta _{k,l}, \\ % \left(\DP{\zeta}{y}\right)_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L} % \exp\left(\frac{2\pi i k }{x_{m}}\right) % \cos\left(\frac{\pi l }{y_{m}}\right) % \frac{\pi l}{y_{m}}\hat \zeta _{k,l}, \\ % \left(\DP{\psi}{x}\right)_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L} % \exp\left(\frac{2\pi i k }{x_{m}}\right) % \sin\left(\frac{\pi l }{y_{m}}\right) % \frac{2\pi i k}{x_{m}}\hat \psi _{k,l}, \\ % \left(\DP{\psi}{y}\right)_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L} % \exp\left(\frac{2\pi i k }{x_{m}}\right) % \cos\left(\frac{\pi l }{y_{m}}\right) % \frac{\pi l}{y_{m}}\hat \psi _{k,l}, \\ % \left(\DP{\rho}{x}\right)_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L} % \exp\left(\frac{2\pi i k }{x_{m}}\right) % \cos\left(\frac{\pi l }{y_{m}}\right) % \frac{2\pi i k}{x_{m}}\hat \rho _{k,l}, \\ % \left(\DP{\rho}{y}\right)_{i,j} &=& - \sum _{k=-K}^{K}\sum _{l=0}^{L} % \exp\left(\frac{2\pi i k }{x_{m}}\right) % \sin\left(\frac{\pi l }{y_{m}}\right) % \frac{\pi l}{y_{m}}\hat \rho _{k,l} % \end{eqnarray} % である. % \subsection{時間積分} % ここでは時間積分法について記述し, プログラム内で実際に用いられている方 % 程式を記述する. 以下では $\Delta t$ を時間格子間隔, 時刻 $\tau \Delta t$ % における $\hat \zeta_{k,l}$ の値を $\hat \zeta_{k,l}^{\tau}$ 等と表す. % 時間方向の離散化は Euler スキームを用いて行う. 時空間方向に離散化され % た方程式は以下のように表される. % \begin{eqnarray} % &&\hat \zeta _{k,l}^{\tau +1} = \hat \zeta _{k,l}^{\tau} + \Delta t \hat F_{k,l}^{\tau}, \\ % &&\hat \rho _{k,l}^{\tau +1} = \hat \rho _{k,l}^{\tau} + \Delta t \hat G_{k,l}^{\tau}, % \end{eqnarray} % \begin{eqnarray} % \hat F_{k,l}^{\tau} &=& % - \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1} % \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right) % \sin\left(\frac{\pi l y_{j}}{y_{m}}\right) % \left(\DP{\psi}{x}\right)_{i,j}^{\tau} % \left(\DP{\zeta}{y}\right)_{i,j}^{\tau} \nonumber \\ % && + \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1} % \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right) % \sin\left(\frac{\pi l y_{j}}{y_{m}}\right) % \left(\DP{\psi}{y}\right)_{i,j}^{\tau} % \left(\DP{\zeta}{x}\right)_{i,j}^{\tau} \nonumber \\ % && - \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1} % \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right) % \sin\left(\frac{\pi l y_{j}}{y_{m}}\right) % \frac{g}{\rho _{0}}\left(\DP{\rho}{x}\right)_{i,j}^{\tau} \nonumber \\ % && - \nu \left[ % \left(\frac{2\pi k}{x_{m}}\right)^{2}+ % \left(\frac{\pi l}{y_{m}}\right)^{2}\right]\hat \zeta _{k,l}^{\tau}, \\ % \hat G_{i,j}^{\tau} &=& - % \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1} % \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right) % \cos\left(\frac{\pi l y_{j}}{y_{m}}\right) % \left(\DP{\psi}{x}\right)_{i,j}^{\tau} % \left(\DP{\rho}{y}\right)_{i,j}^{\tau} \nonumber \\ % && + \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1} % \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right) % \cos\left(\frac{\pi l y_{j}}{y_{m}}\right) % \left(\DP{\psi}{y}\right)_{i,j}^{\tau} % \left(\DP{\rho}{x}\right)_{i,j}^{\tau} \nonumber \\ % && - \kappa \left[ % \left(\frac{2\pi k}{x_{m}}\right)^{2}+ % \left(\frac{\pi l}{y_{m}}\right)^{2}\right]\hat \rho _{k,l}^{\tau}. % \end{eqnarray} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{使用モジュールとその他の設定} スペクトル変換と逆変換, 微分演算は SPMODEL ライブラリ (spml) の et\_module に含まれる関数を用いて行う. フーリエ変換とチェビシェフ変換, それらの逆変換の際の数値積分は台形公式を用いて行う. spml が下位で使用する ISPACK の仕様から, 格子点数 $I,J$ は偶数で, かつ $I/2, J/2=2^{a}3^{b}5^{c}$ ($a, b, c$ は 0 または整数) でなければならな い. 非線形項の計算によって生じるエリアジングを防ぐため, 格子点数 $I,J$ と切断波数 $K,L$ は $I>3K,J>3K/2$ を満たすように与える. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{数値実験} 数値実験では重力流の室内実験を模した計算を行う. 初期条件として以下のよ うな密度分布を持つ静止した流体を考える. \begin{equation} \rho = \rho_{0} + \frac{\delta \rho}{2} \left[\mbox{tanh}\left(\frac{x-x_{1}}{\sigma}\right) - \mbox{tanh}\left(\frac{x-x_{2}}{\sigma}\right) \right] \end{equation} パラメータは\Dtabref{使用したパラメータの値}にまとめた値を用いる. 格子点数 $I, J$ と切断波数 $K, L$ はそれぞれ $I=128, J=64, K=L=42$ と する. 時間格子間隔 $\Delta t$ は $10^{-4}$ sec, 計算ステップ数は 300,000 ステップである. \Dfigref{密度分布の時間変化}に計算された密度分布の時間変化を示した. 室 内実験との対比がしやすいように $0\leq x\leq 60$ cm までの範囲を示してあ る. \begin{table}[h] \begin{center} \begin{tabular}{c l l } \hline パラメータ &\qquad\qquad& 数値 \\ \hline $\rho_{0}$ && 1000 kgm${}^{-3}$ \\ $\delta \rho $ && 1 kgm${}^{-3}$ \\ $g$ && 9.8 msec${}^{-2}$ \\ $\nu$ && 10${}^{-7}$ m${}^{2}$sec${}^{-1}$ \\ $\kappa$ && 10${}^{-11}$ m${}^{2}$sec${}^{-1}$ \\ $x_{m}$ && 0.2 m \\ $y_{m}$ && 1.2 m \\ $x_{1}$ && $x_{m}/2$ $+$ 0.01 m \\ $x_{2}$ && $x_{m}/2$ $-$ 0.01 m \\ $\sigma$ && 0.02 m \\\hline \end{tabular} \caption{使用したパラメータの値} \Dtablab{使用したパラメータの値} \end{center} \end{table} \begin{figure}[p] \begin{center} \includegraphics[height=0.8\textheight]{./figs/gravcurrent_tdevelop.png} \caption{密度分布の時間変化. 左上から右下へと順に $t=0.0$ から 7.0 sec 毎の結果. $0\leq x\leq 60$ cm までの範囲を示した. 寒色から暖色へと色が変化するにつれて密度が高い.} \Dfiglab{密度分布の時間変化} \end{center} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \newpage \section{参考文献} \begin{description} \item 竹広真一, 石岡圭一, 森川靖大, 小高正嗣, 石渡正樹, 林祥介, SPMODEL 開発グループ, 2004: 階層的地球流体力学スペクトルモデル集 (SPMODEL), http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash , 地球流体電脳倶楽部. \end{description} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section*{謝辞} % 謝辞. このまま利用する. \markright{謝辞} % ヘッダに表示 \addcontentsline{toc}{section}{謝辞} % 目次に表示 本資源は, 地球流体電脳倶楽部のインターネット上での学術知識の集積と活用 の実験の一環として \begin{center} http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash \end{center} において公開されているものである (\copyright 地球流体電脳倶楽部スペクトルモデルプロジェクト spmodel@gfd-dennou.org 2002. ). 本資源は, 著作者の諸権利に抵触しない(迷惑をかけない)限りにおいて自由に 利用していただいて構わない. なお, 利用する際には今一度自ら内容を確かめ ることをお願いする(無保証無責任原則). 本資源に含まれる元資源提供者(図等の版元等を含む)からは, 直接的な形で の WEB 上での著作権または使用許諾を得ていない場合があるが, 勝手ながら, 「未来の教育」のための実験という学術目的であることをご理解いただけるも のと信じ, 学術標準の引用手順を守ることで諸手続きを略させていただいて いる. 本資源の利用者には, この点を理解の上, 注意して扱っていただける ようお願いする. 万一, 不都合のある場合には \begin{center} spmodel@gfd-dennou.org \end{center} まで連絡していただければ幸いである. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \end{document} %%% Local Variables: %%% mode: latex %%% TeX-master: t %%% End: