!-- !---------------------------------------------------------------------- ! Copyright (c) 2001-2013 SPMODEL Development Group, All rights reserved. !---------------------------------------------------------------------- ! !表題 esc_module ! ! spml/esc_module モジュールは 2 次元チャネル領域での流体運動を ! スペクトル法により数値計算するための Fortran90 関数を提供する. ! ! 内部で ISPACK/C2PACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法については ! ISPACK/C2PACK のマニュアルを参照されたい. ! !履歴 2001/10/07 竹広真一 ! 2001/12/26 竹広真一 関数, 変数命名法変更 ! 2002/03/25 竹広真一 モジュール名変更 ! 2002/08/19 竹広真一 格子データ添字を gg -> xy に変更 ! 2002/08/20 竹広真一 積分・平均関数追加 ! 2005/03/15 竹広真一 xy -> yx に接頭子を変更 ! 2006/03/06 竹広真一 コメントを RDoc 用に修正 ! 2006/03/19 竹広真一 変数・手続き群の要約をコメントに追加 ! 2007/11/12 竹広真一 ec_LaplaInv_ec 関数追加 ! 2007/11/21 竹広真一 初期化サブルーチンメッセージ出力 ! 2009/01/09 竹広真一 esc_Initial メッセージに日付を追加 ! 2009/01/29 佐々木洋平 コメントを RDoc 用に修正 ! 2009/07/30 竹広真一 作業領域をローカル変数に変更(for OpenMP) ! 2013/09/26 佐々木洋平 tj,ti を倍精度に変更 ! !++ module esc_module ! != esc_module ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id$ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要 ! ! spml/esc_module モジュールは 2 次元チャネル領域での流体運動を ! スペクトル法により数値計算するための Fortran90 関数を提供する. ! ! 内部で ISPACK/C2PACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法については ! ISPACK/C2PACK のマニュアルを参照されたい. ! !== 関数・変数の名前と型について ! !=== 命名法 ! ! * 関数名の先頭 (es_, ec_, yx_, x_, y_) は, 返す値の形を示している. ! es_,ec_ :: スペクトルデータ(Y 方向 SIN 展開, COS 展開) ! yx_ :: 2 次元格子点データ ! x_ :: X 方向 1 次元格子点データ ! y_ :: Y 方向 1 次元格子点データ ! ! * 関数名の間の文字列(Dx, Dy, Lapla, LaplaInv, Jacobian)は, ! その関数の作用を表している. ! ! * 関数名の最後 (_es_es, _es_ec, _es, _ec, _yx, _x, _y) は, ! 入力変数の形スペクトルデータおよび格子点データであることを示している. ! _es :: スペクトルデータ(Y 方向 SIN 型) ! _ec :: スペクトルデータ(Y 方向 COS 型) ! _es_es :: 2 つのスペクトルデータ ! _es_ec :: 2 つのスペクトルデータ ! _yx :: 2 次元格子点データ, ! _x :: X 方向 1 次元格子点データ ! _y :: Y 方向 1 次元格子点データ ! !=== 各データの種類の説明 ! ! * yx : 2 次元格子点データ. ! * 変数の種類と次元は real(8), dimension(0:jm,0:im-1). ! * im, jm はそれぞれ X, Y 座標の格子点数であり, ! サブルーチン esc_Initial にてあらかじめ設定しておく. ! * 第 1 次元が Y 座標の格子点位置番号, 第 2 次元が X 座標の格子点位置番号 ! である (X, Y の順ではない)ことに注意. ! ! * es : X 方向フーリエ型, Y 方向 SIN 型スペクトルデータ. ! * 変数の種類と次元は real(8), dimension(-km:km,lm). ! * km, lm はそれぞれ X, Y 方向の最大波数であり, ! サブルーチン esc_Initial にて ! あらかじめ設定しておく. ! * スペクトルデータの格納のされ方については... ! ! * ec : X 方向フーリエ型, Y 方向 COS 型スペクトルデータ. ! * 変数の種類と次元は real(8), dimension(-km:km,0:lm). ! * km, lm はそれぞれ ! X, Y 方向の最大波数であり, ! サブルーチン esc_Initial にてあらかじめ設定しておく. ! * スペクトルデータの格納のされ方については... ! ! * x, y : X, Y 方向 1 次元格子点データ. ! * 変数の種類と次元はそれぞれ real(8), dimension(0:im-1) ! および real(8), dimension(0:jm). ! ! * es_, ec_ で始まる関数が返す値はスペクトルデータに同じ. ! ! * yx_ で始まる関数が返す値は 2 次元格子点データに同じ. ! ! * x_, y_ で始まる関数が返す値は 1 次元格子点データに同じ. ! ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに ! 微分などを作用させたデータをスペクトル変換したものことである. ! !== 変数・手続き群の要約 ! !==== 初期化 ! ! esc_initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定 ! !==== 座標変数 ! ! x_X, y_Y :: 格子点座標(X,Y座標)を格納した 1 次元配列 ! x_X_Weight, y_Y_Weight :: 重み座標を格納した 1 次元配列 ! yx_X, yx_Y :: 格子点データの XY 座標(X,Y)(格子点データ型 2 次元配列) ! !==== 基本変換 ! ! yx_es, yx_ec :: スペクトルデータから格子データへの変換 ! es_yx, ec_yx :: 格子データからスペクトルデータへの変換 ! !==== 微分 ! ! es_Lapla_es, ec_Lapla_ec :: スペクトルデータにラプラシアンを作用させる ! es_LaplaInv_es, ec_LaplaInv_ec :: スペクトルデータに ! ラプラシアンの逆変換を作用させる ! es_Dx_es, ec_Dx_ec :: スペクトルデータに X 微分を作用させる ! ec_Dy_es, es_Dy_ec :: スペクトルデータに Y 微分を作用させる ! es_Jacobian_es_es, ec_Jacobian_es_ec :: 2 つのスペクトルデータから ! ヤコビアンを計算する ! !==== 積分・平均 ! ! IntYX_yx, AvrYX_yx :: 2 次元格子点データの全領域積分および平均 ! y_IntX_yx, y_AvrX_yx :: 2 次元格子点データの X 方向積分および平均 ! IntX_x, AvrX_x :: 1 次元(X)格子点データの X 方向積分および平均 ! x_IntY_yx, x_AvrY_yx :: 2 次元格子点データの Y 方向積分および平均 ! IntY_y, AvrY_y :: 1 次元(Y)格子点データの Y 方向積分および平均 ! ! use dc_message, only : MessageNotify use dc_types, only: DP implicit none private public esc_Initial ! 初期化ルーチン public yx_es, yx_ec, es_yx, ec_yx ! 基本変換 public es_Lapla_es, es_LaplaInv_es, es_Dx_es, ec_Dy_es ! 微分 public ec_Lapla_ec, ec_LaplaInv_ec, ec_Dx_ec, es_Dy_ec ! 微分 public es_Jacobian_es_es, ec_Jacobian_es_ec ! 非線形計算 public IntYX_yx, y_IntX_yx, x_IntY_yx, IntX_x, IntY_y ! 積分 public AvrYX_yx, y_AvrX_yx, x_AvrY_yx, AvrX_x, AvrY_y ! 平均 public x_X, y_Y, x_X_Weight, y_Y_Weight, yx_X, yx_Y ! 座標変数 integer :: im=32, jm=8 ! 格子点の設定(X,Y) integer :: km=10, lm=5 ! 切断波数の設定(X,Y) real(DP) :: xl=2.0, yl=1.0 ! 領域の大きさ integer, dimension(5) :: itj real(DP), dimension(:),allocatable :: tj integer, dimension(5) :: iti real(DP), dimension(:),allocatable :: ti real(DP), dimension(:), allocatable :: x_X ! 格子点座標(X)を格納した 1 次元配列 real(DP), dimension(:), allocatable :: y_Y ! 格子点座標(X)を格納した 1 次元配列 real(DP), dimension(:), allocatable :: x_X_Weight ! 重み座標(X)を格納した 1 次元配列. X 方向の格子点の間隔が格納してある. real(DP), dimension(:), allocatable :: y_Y_Weight ! 重み座標(Y)を格納した 1 次元配列. Y 方向の格子点の間隔が格納してある. real(DP), dimension(:,:), allocatable :: yx_X ! 各格子点(i,j)の位置の X 座標を格納した格子データ. real(DP), dimension(:,:), allocatable :: yx_Y ! 各格子点(i,j)の位置の Y 座標を格納した格子データ. real(DP), dimension(:), allocatable :: wg, ws, wgj real(DP), dimension(:,:), allocatable :: yx_work,es_work,ec_work real(DP), parameter :: pi=3.1415926535897932385D0 save im, jm, km, lm, itj, tj, iti, ti, xl, yl save x_X, y_Y, x_X_Weight, y_Y_Weight, yx_X, yx_Y contains !--------------- 初期化 ----------------- subroutine esc_Initial(i,j,k,l,xmin,xmax,ymin,ymax) ! ! スペクトル変換の格子点数, 波数, 領域の大きさを設定する. ! ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで ! 初期設定をしなければならない. ! integer,intent(in) :: i !(in) 格子点数(X) integer,intent(in) :: j !(in) 格子点数(Y) integer,intent(in) :: k !(in) 切断波数(X) integer,intent(in) :: l !(in) 切断波数(Y) real(DP),intent(in) :: xmin, xmax !(in) X 座標範囲 real(DP),intent(in) :: ymin, ymax !(in) Y 座標範囲 integer :: ii, jj im = i ; jm = j km = k ; lm = l xl = xmax-xmin ; yl = ymax-ymin allocate(tj(jm*6),ti(im*2)) call c2init(jm,im,itj,tj,iti,ti) allocate(x_X(0:im-1), x_X_Weight(0:im-1)) allocate(y_Y(0:jm), y_Y_Weight(0:jm)) allocate(yx_X(0:jm,0:im-1), yx_Y(0:jm,0:im-1)) do ii=0,im-1 x_X(ii) = xmin + xl/im*ii enddo x_X_Weight = xl/im do jj=0,jm y_Y(jj) = ymin + yl/jm*jj enddo y_Y_Weight(0) = yl/(2*jm) y_Y_Weight(1:jm-1) = yl/jm y_Y_Weight(jm) = yl/(2*jm) yx_X = spread(x_X,1,jm+1) yx_Y = spread(y_Y,2,im) call MessageNotify('M','esc_initial', & 'esc_module (2009/07/30) is initialized') end subroutine esc_Initial !--------------- 基本変換 ----------------- function yx_es(es) ! ! SIN(Y)型スペクトルデータから格子データへ変換する. ! real(DP), dimension(0:jm,0:im-1) :: yx_es !(out) 格子点データ real(DP), dimension(-km:km,lm), intent(in) :: es !(in) SIN(Y)型スペクトルデータ real(DP) :: wg((jm+1)*im) ! 作業領域 call c2s2ga(lm,km,jm,im,es,yx_es,wg,itj,tj,iti,ti,1) end function yx_es function yx_ec(ec) ! ! COS(Y)型スペクトルデータから格子データへ変換する. ! real(DP), dimension(0:jm,0:im-1) :: yx_ec !(out) 格子点データ real(DP), dimension(-km:km,0:lm), intent(in) :: ec !(in) COS(Y)型スペクトルデータ real(DP) :: wg((jm+1)*im) ! 作業領域 call c2s2ga(lm,km,jm,im,ec,yx_ec,wg,itj,tj,iti,ti,2) end function yx_ec function es_yx(yx) ! ! 格子データから SIN(Y)型スペクトルデータへ変換する. ! real(DP), dimension(-km:km,lm) :: es_yx !(out) SIN(Y)型スペクトルデータ real(DP), dimension(0:jm,0:im-1), intent(in) :: yx !(in) 格子点データ real(DP) :: wg((jm+1)*im) real(DP) :: yx_work(0:jm,0:im-1) ! 作業領域 yx_work = yx call c2g2sa(lm,km,jm,im,yx_work,es_yx,wg,itj,tj,iti,ti,1) end function es_yx function ec_yx(yx) ! ! 格子データから COS(Y)型スペクトルデータへ変換する. ! real(DP), dimension(-km:km,0:lm) :: ec_yx !(in) COS(Y)型スペクトルデータ real(DP), dimension(0:jm,0:im-1), intent(in) :: yx !(out) 格子点データ real(DP) :: wg((jm+1)*im) real(DP) :: yx_work(0:jm,0:im-1) ! 作業領域 yx_work = yx call c2g2sa(lm,km,jm,im,yx_work,ec_yx,wg,itj,tj,iti,ti,2) end function ec_yx !--------------- 微分計算 ----------------- function es_Lapla_es(es) ! ! SIN(Y)型入力スペクトルデータにラプラシアン(∂xx+∂yy)を作用する. ! ! スペクトルデータのラプラシアンとは, 対応する格子点データに ! ラプラシアンを作用させたデータのスペクトル変換のことである. ! ! 実際にはスペクトルデータに全波数 (k**2 + l**2) をかける ! 計算を行っている. ! real(DP), dimension(-km:km,lm) :: es_Lapla_es !(out) スペクトルデータのラプラシアン real(DP), dimension(-km:km,lm), intent(in) :: es !(in) SIN(Y)型入力スペクトルデータ integer k,l do l=1,lm do k=-km,km es_Lapla_es(k,l) = -((2*pi*k/xl)**2+(pi*l/yl)**2)*es(k,l) enddo enddo end function es_Lapla_es function ec_Lapla_ec(ec) ! ! COS(Y)型入力スペクトルデータにラプラシアン(∂xx+∂yy)を作用する. ! ! スペクトルデータのラプラシアンとは, 対応する格子点データに ! ラプラシアンを作用させたデータのスペクトル変換のことである. ! ! 実際にはスペクトルデータに全波数 (k**2 + l**2) をかける ! 計算を行っている. ! real(DP), dimension(-km:km,0:lm) :: ec_Lapla_ec !(out) スペクトルデータのラプラシアン real(DP), dimension(-km:km,0:lm), intent(in) :: ec !(in) COS(Y)型入力スペクトルデータ integer k,l do l=0,lm do k=-km,km ec_Lapla_ec(k,l) = -((2*pi*k/xl)**2+(pi*l/yl)**2)*ec(k,l) enddo enddo end function ec_Lapla_ec function es_LaplaInv_es(es) ! スペクトル SINY に作用する逆 \lapla 演算子 ! ! 入力スペクトルデータに逆ラプラシアン(∂xx+∂yy)**(-1)を作用する. ! ! スペクトルデータの逆ラプラシアンとは, 対応する格子点データに ! 逆ラプラシアンを作用させたデータのスペクトル変換のことである. ! ! 実際にはスペクトルデータに全波数 (k**2 + l**2) で割る ! 計算を行っている. ! real(DP), dimension(-km:km,lm) :: es_LaplaInv_es !(out) スペクトルデータの逆ラプラシアン real(DP), dimension(-km:km,lm), intent(in) :: es !(in) SIN(Y)型スペクトルデータ integer k,l do l=1,lm do k=-km,km es_LaplaInv_es(k,l) = -es(k,l)/((2*pi*k/xl)**2+(pi*l/yl)**2) enddo enddo end function es_LaplaInv_es function ec_LaplaInv_ec(ec) ! スペクトル COSY に作用する逆 \lapla 演算子 ! ! 入力スペクトルデータに逆ラプラシアン(∂xx+∂yy)**(-1)を作用する. ! ! スペクトルデータの逆ラプラシアンとは, 対応する格子点データに ! 逆ラプラシアンを作用させたデータのスペクトル変換のことである. ! ! 実際にはスペクトルデータに全波数 (k**2 + l**2) で割る ! 計算を行っている. ! real(DP), dimension(-km:km,0:lm) :: ec_LaplaInv_ec !(out) スペクトルデータの逆ラプラシアン real(DP), dimension(-km:km,0:lm), intent(in) :: ec !(in) SIN(Y)型スペクトルデータ integer k,l if ( abs(ec(0,0)) >= epsilon(0.0D0) ) then call MessageNotify('W','ec_LaplaInv_ec', & '0-0 component of input data is not zero.') call MessageNotify('W','ec_LaplaInv_ec', & '0-0 component of output set to zero.') endif do l=1,lm do k=-km,km ec_LaplaInv_ec(k,l) = -ec(k,l)/((2*pi*k/xl)**2+(pi*l/yl)**2) enddo enddo do k=1,km ec_LaplaInv_ec(k,0) = -ec(k,0)/(2*pi*k/xl)**2 ec_LaplaInv_ec(-k,0) = -ec(-k,0)/(2*pi*k/xl)**2 enddo ec_LaplaInv_ec(0,0) = 0.0 end function ec_LaplaInv_ec function es_Dx_es(es) ! ! SIN(Y)型入力スペクトルデータに X 微分(∂x)を作用する. ! ! スペクトルデータの X 微分とは, 対応する格子点データに X 微分を ! 作用させたデータのスペクトル変換のことである. ! ! 実際にはスペクトルデータに X 方向波数 k をかけて ! sin(kx) <-> cos(kx) 成分に入れ換える計算を行っている. ! real(DP), dimension(-km:km,lm) :: es_Dx_es !(out) スペクトルデータの X 微分 real(DP), dimension(-km:km,lm), intent(in) :: es !(in) SIN(Y)型入力スペクトルデータ integer k,l do l=1,lm do k=-km,km es_Dx_es(k,l) = (-2*pi*k/xl)*es(-k,l) enddo enddo end function es_Dx_es function ec_Dx_ec(ec) ! ! COS(Y)型入力スペクトルデータに X 微分(∂x)を作用する. ! ! スペクトルデータの X 微分とは, 対応する格子点データに X 微分を ! 作用させたデータのスペクトル変換のことである. ! ! 実際にはスペクトルデータに X 方向波数 k をかけて ! sin(kx) <-> cos(kx) 成分に入れ換える計算を行っている. ! real(DP), dimension(-km:km,0:lm) :: ec_Dx_ec !(out) スペクトルデータの X 微分 real(DP), dimension(-km:km,0:lm), intent(in) :: ec !(in) COS(Y)型入力スペクトルデータ integer k,l do l=0,lm do k=-km,km ec_Dx_ec(k,l) = (-2*pi*k/xl)*ec(-k,l) enddo enddo end function ec_Dx_ec function ec_Dy_es(es) ! スペクトル SINY に作用する y 微分演算子 ! ! SIN(Y)型入力スペクトルデータに Y 微分(∂y)を作用する. ! ! スペクトルデータの X 微分とは, 対応する格子点データに Y 微分を ! 作用させたデータのスペクトル変換のことである. ! ! 実際にはスペクトルデータに X 方向波数 l をかけている. ! real(DP), dimension(-km:km,0:lm) :: ec_Dy_es !(out) スペクトルデータの Y 微分, COS(Y)型. real(DP), dimension(-km:km,lm), intent(in) :: es !(in) SIN(Y)型入力スペクトルデータ integer k,l do k=-km,km ec_Dy_es(k,0) = 0.0 enddo do l=1,lm do k=-km,km ec_Dy_es(k,l) = (pi*l/yl)*es(k,l) enddo enddo end function ec_Dy_es function es_Dy_ec(ec) ! スペクトル COSY に作用する y 微分演算子 ! ! COS(Y)型入力スペクトルデータに Y 微分(∂y)を作用する. ! ! スペクトルデータの X 微分とは, 対応する格子点データに Y 微分を ! 作用させたデータのスペクトル変換のことである. ! ! 実際にはスペクトルデータに X 方向波数 l をかけている. ! real(DP), dimension(-km:km,lm) :: es_Dy_ec !(out) スペクトルデータの Y 微分, SIN(Y)型. real(DP), dimension(-km:km,0:lm), intent(in) :: ec !(in) COS(Y)型入力スペクトルデータ integer k,l do l=1,lm do k=-km,km es_Dy_ec(k,l) = -(pi*l/yl)*ec(k,l) enddo enddo end function es_Dy_ec !------------------- 非線形項計算 ---------------------- function es_Jacobian_es_es(es_a,es_b) !スペクトル SINY に作用するヤコビアン ! ! 2 つのスペクトルデータからヤコビアン ! ! J(A,B)=(∂xA)(∂yB)-(∂yA)(∂xB) ! ! を計算する. ! ! 2 つのスペクトルデータのヤコビアンとは, 対応する 2 つの ! 格子点データのヤコビアンのスペクトル変換のことである. ! real(DP), dimension(-km:km,lm) :: es_Jacobian_es_es !(out) 2 つのスペクトルデータのヤコビアン real(DP), dimension(-km:km,lm), intent(in) :: es_A,es_B !(in) 2つの SIN(Y)型入力スペクトルデータ real(DP) :: ws((2*km+1)*(lm+1)) real(DP) :: wgj((jm+1)*im*3) real(DP) :: es_work(-km:km,lm) ! 作業領域 integer k,l call c2ajcb(lm,km,jm,im,es_A,es_B,es_work,ws,wgj,itj,tj,iti,ti) do l=1,lm do k=-km,km es_Jacobian_es_es(k,l) = (2*pi/xl)*(pi/yl)*es_work(k,l) enddo enddo end function es_Jacobian_es_es function ec_Jacobian_es_ec(es,ec) ! スペクトル COS(Y) に作用するヤコビアン ! ! 2 つのスペクトルデータからヤコビアン ! ! J(A,B)=(∂xA)(∂yB)-(∂yA)(∂xB) ! ! を計算する. ! ! 2 つのスペクトルデータのヤコビアンとは, 対応する 2 つの ! 格子点データのヤコビアンのスペクトル変換のことである. ! real(DP), dimension(-km:km,0:lm) :: ec_Jacobian_es_ec !(out) 2 つのスペクトルデータのヤコビアン real(DP), dimension(-km:km,lm), intent(in) :: es !(in) 1つ目の SIN(Y)型入力スペクトルデータ real(DP), dimension(-km:km,0:lm), intent(in) :: ec !(in) 2つ目の COS(Y)型入力スペクトルデータ real(DP) :: ws((2*km+1)*(lm+1)) real(DP) :: wgj((jm+1)*im*3) real(DP) :: ec_work(-km:km,0:lm) ! 作業領域 integer k,l call c2ajcc(lm,km,jm,im,es,ec,ec_work,ws,wgj,itj,tj,iti,ti) do l=0,lm do k=-km,km ec_Jacobian_es_ec(k,l) = (2*pi/xl)*(pi/yl)*ec_work(k,l) enddo enddo end function ec_Jacobian_es_ec !--------------- 積分計算 ----------------- function IntYX_yx(yx) ! ! 2 次元格子点データの全領域積分および平均. ! real(DP), dimension(0:jm,0:im-1) :: yx !(in) 2 次元格子点データ real(DP) :: IntYX_yx !(out) 積分値 integer :: i, j IntYX_yx = 0.0d0 do i=0,im-1 do j=0,jm IntYX_yx = IntYX_yx + yx(j,i) * y_Y_Weight(j) * x_X_Weight(i) enddo enddo end function IntYX_yx function y_IntX_yx(yx) ! X 方向積分 ! ! 1 次元(X)格子点データの X 方向積分 ! real(DP), dimension(0:jm,0:im-1) :: yx !(in) 2 次元格子点データ real(DP), dimension(0:jm) :: y_IntX_yx !(out) X 方向に積分された 1 次元(Y)格子点データ integer :: i y_IntX_yx = 0.0d0 do i=0,im-1 y_IntX_yx(:) = y_IntX_yx(:) + yx(:,i) * x_X_Weight(i) enddo end function y_IntX_yx function x_IntY_yx(yx) ! ! 2 次元格子点データの Y 方向積分 ! real(DP), dimension(0:jm,0:im-1) :: yx !(in) 2 次元格子点データ real(DP), dimension(0:im-1) :: x_IntY_yx !(out) Y 方向に積分された 1 次元(X)格子点データ integer :: j x_IntY_yx = 0.0d0 do j=0,jm x_IntY_yx(:) = x_IntY_yx(:) + yx(j,:) * y_Y_Weight(j) enddo end function x_IntY_yx function IntX_x(x) ! X 方向積分 ! ! 1 次元(X)格子点データの X 方向積分 ! real(DP), dimension(0:im-1) :: x !(in) 1 次元(X)格子点データ real(DP) :: IntX_x !(out) 積分値 IntX_x = sum(x*x_X_Weight) end function IntX_x function IntY_y(y) ! Y 方向積分 ! ! 1 次元(Y)格子点データの Y 方向積分 ! real(DP), dimension(0:jm) :: y !(in) 1 次元(Y)格子点データ real(DP) :: IntY_y !(out) 積分値 IntY_y = sum(y*y_Y_Weight) end function IntY_y !--------------- 平均計算 ----------------- function AvrYX_yx(yx) ! 全領域平均 ! ! 2 次元格子点データの全領域平均 ! real(DP), dimension(0:jm,0:im-1) :: yx !(in) 2 次元格子点データ real(DP) :: AvrYX_yx !(out) 平均値 AvrYX_yx = IntYX_yx(yx)/(sum(x_X_weight)*sum(y_Y_weight)) end function AvrYX_yx function y_AvrX_yx(yx) ! X 方向平均 ! ! 1 次元(X)格子点データの X 方向平均 ! real(DP), dimension(0:jm,0:im-1) :: yx !(in) 2 次元格子点データ real(DP), dimension(0:jm) :: y_AvrX_yx !(out) X 方向に平均された 1 次元(Y)格子点データ y_AvrX_yx = y_IntX_yx(yx)/sum(x_X_weight) end function y_AvrX_yx function x_AvrY_yx(yx) ! Y 方向平均 ! ! 2 次元格子点データの Y 方向平均 ! real(DP), dimension(0:jm,0:im-1) :: yx !(in) 2 次元格子点データ real(DP), dimension(0:im-1) :: x_AvrY_yx !(out) Y 方向に平均された 1 次元(X)格子点データ x_AvrY_yx = x_IntY_yx(yx)/sum(y_Y_weight) end function x_AvrY_yx function AvrX_x(x) ! X 方向平均 ! ! 1 次元(X)格子点データの X 方向平均 ! real(DP), dimension(0:im-1) :: x !(in) 1 次元(X)格子点データ real(DP) :: AvrX_x !(out) 平均値 AvrX_x = IntX_x(x)/sum(x_X_weight) end function AvrX_x function AvrY_y(y) ! Y 方向平均 ! ! 1 次元(Y)格子点データの Y 方向平均 ! real(DP), dimension(0:jm) :: y !(in) 1 次元(Y)格子点データ real(DP) :: AvrY_y ! 平均値 AvrY_y = IntY_y(y)/sum(y_Y_weight) end function AvrY_y end module esc_module