!-- !---------------------------------------------------------------------- ! Copyright (c) 2002-2013 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- !表題 w_deriv_module ! ! spml/w_deriv_module モジュールは球面上での 2 次元流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するための ! モジュール w_module の下部モジュールであり, スペクトル法の ! 微分計算のための Fortran90 関数を提供する. ! ! 内部で ISPACK の SPPACK と SNPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK/SNPACK,SPPACK のマニュアルを参照されたい. ! ! このモジュールを使うためには前もって w_base_initial を呼んで ! 切断波数, 格子点数の設定をしておく必要がある. ! ! !履歴 2001/12/08 竹広真一 ! 2001/12/26 竹広真一 関数,変数名前変更 ! 2002/02/07 竹広真一 関数,変数名前再変更 ! 2002/03/30 竹広真一 関数,変数名前再再変更 ! 2002/05/25 竹広真一 格子点座標を区別すべく命名法変更 ! 2005/07/04 竹広真一 OPENMP 版変換ルーチンに対応 ! 2006/03/08 竹広真一 コメントを RDoc 用に修正 ! 2008/05/31 竹広真一 初期化ルーチンを分離 ! 2008/06/22 佐々木洋平 格子点データの配列開始点を 1 から 0 へ. ! 2008/06/23 佐々木洋平 格子点データの格納は (0:im-1, 1:jm) へ. ! 2008/07/01 佐々木洋平 コメントを RDoc 用に微修正 ! 2009/01/09 竹広真一 w_deriv_Initial メッセージに日付を追加 ! 2009/01/29 佐々木洋平 コメントを RDoc 用に微修正 ! 2009/07/30 竹広真一 作業領域をローカル変数に変更(for OpenMP) ! 2013/02/23 竹広真一 w_deriv_Finalize 導入 ! ! 制限 ! ・変換する格子点データ, スペクトルデータの配列の大きさは決めうち ! !++ module w_deriv_module ! != w_deriv_module ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id$ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要 ! ! spml/w_deriv_module モジュールは球面上での 2 次元流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するための ! モジュール w_module の下部モジュールであり, スペクトル法の ! 微分計算のための Fortran90 関数を提供する. ! ! 内部で ISPACK の SPPACK と SNPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK/SNPACK,SPPACK のマニュアルを参照されたい. ! ! このモジュールを使うためには前もって w_base_initial を呼んで ! 切断波数, 格子点数の設定をしておく必要がある. ! use dc_message, only : MessageNotify use w_base_module, only : im, jm, nm, it, t, y, ip, p, r, ia, a, & w_base_Initial, xy_w, w_xy implicit none real(8), allocatable :: rn(:,:) ! ラプラシアン演算用配列 ! ! スペクトルデータのラプラシアンを計算するための係数 ! 配列のサイズは((nm+1)*(nm+1), 2) ! ! r(L,1) には L 番目の格納位置のスペクトルに対するラプラシアン計算の ! 係数 -n(n+1) の値が格納されている. ! integer, allocatable :: irm(:,:) ! 経度微分演算用配列 ! ! スペクトルデータの経度微分を計算するための係数. ! 配列サイズは ( (nm+1)*(nm+1),2 ) である. ! ! L番目の格納位置のスペクトルが実部なら, irm(L,1)には対応する虚部の格納位置が, ! irm(L,2) には東西波数 m が格納されている. また, L番目の格納位置のスペクトル ! が虚部なら, irm(L,1)には対応する実部の格納位置が, irm(L,2)には -m が格納され ! ている. ! integer, allocatable :: ip2(:), ip3(:) ! ヤコビアン計算用配列 real(8), allocatable :: p2(:), p3(:) ! ヤコビアン計算用配列 real(8), allocatable :: r2(:), r3(:) ! ヤコビアン計算用配列 integer iw ! 作業配列の大きさ logical :: w_deriv_initialize=.false. ! 初期化フラッグ private public w_deriv_Initial ! 初期化 public w_deriv_Finalize ! 終了処理 public w_Lapla_w, w_LaplaInv_w ! ラプラシアンと逆演算 public w_DLon_w ! 経度微分 public xy_GradLon_w, xy_GradLat_w ! 勾配型微分 public w_DivLon_xy, w_DivLat_xy ! 発散型微分 public w_Div_xy_xy ! 発散型微分 public w_Jacobian_w_w ! ヤコビアン public xy_GradLambda_w, xy_GradMu_w ! 勾配型微分(λ,μ座標) public w_DivLambda_xy, w_DivMu_xy ! 発散型微分(λ,μ座標) public rn, irm ! ラプラシアン/経度微分演算用配列 save rn, irm, ip2, ip3, p2, p3, r2, r3 save iw save w_deriv_initialize ! 初期化フラグ contains !--------------- 初期化 ----------------- subroutine w_deriv_initial ! ! スペクトル微分計算に必要となる作業領域を設定する. ! ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで ! 初期設定をしなければならない. ! ! このサブルーチンを単独で用いるのでなく, ! 上位サブルーチン w_Initial を使用すること. ! allocate(rn((nm+1)*(nm+1),2)) ! ラプラシアン演算用配列 allocate(irm((nm+1)*(nm+1),2)) ! 経度微分演算用配列 call spnini(nm,rn) call spmini(nm,irm) allocate(ip2(2*((nm+1)/2+nm+1)*2)) ! ヤコビアン計算用配列 allocate(p2(2*((nm+1)/2+nm+1)*jm)) ! ヤコビアン計算用配列 allocate(r2(2*((nm+1)/2*2+3)*(nm/2+1))) ! ヤコビアン計算用配列 allocate(ip3(3*((nm+1)/2+nm+1)*2)) ! ヤコビアン計算用配列 allocate(p3(3*((nm+1)/2+nm+1)*jm)) ! ヤコビアン計算用配列 allocate(r3(3*((nm+1)/2*2+3)*(nm/2+1))) ! ヤコビアン計算用配列 call snkini(nm,jm,2,ip,p,r,ip2,p2,r2) call snkini(nm,jm,3,ip,p,r,ip3,p3,r3) iw=3*max( ((nm+1)/2*2+3)*(nm/2+2)*2, & jm*((nm+1)/2+nm+1)*2, jm*im ) w_deriv_initialize=.true. call MessageNotify('M','w_deriv_initial',& 'w_deriv_module (2013/02/23) is initialized') end subroutine w_deriv_initial !--------------- 微分計算 ----------------- function w_Lapla_w(w_data) ! ! 入力スペクトルデータにラプラシアン ! ! ▽^2 = 1/cos^2φ・∂^2/∂λ^2 + 1/cosφ・∂/∂φ(cosφ∂/∂φ) ! ! を作用する(1 層用). ! ! スペクトルデータのラプラシアンとは, 対応する格子点データに ! ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8) :: w_Lapla_w((nm+1)*(nm+1)) !(out) 入力スペクトルデータのラプラシアン real(8), intent(in) :: w_data((nm+1)*(nm+1)) !(in) 入力スペクトルデータ call spclap(nm,w_data,w_Lapla_w,rn(1,1)) end function w_Lapla_w function w_LaplaInv_w(w_data) ! ! 入力スペクトルデータに逆ラプラシアン ! ! ▽^{-2} ! =[1/cos^2φ・∂^2/∂λ^2 + 1/cosφ・∂/∂φ(cosφ∂/∂φ)]^{-1} ! ! を作用する(1 層用). ! ! スペクトルデータの逆ラプラシアンとは, 対応する格子点データに ! 逆ラプラシアンを作用させたデータのスペクトル変換のことである. ! real(8) :: w_LaplaInv_w((nm+1)*(nm+1)) !(out) スペクトルデータの逆ラプラシアン real(8), intent(in) :: w_data((nm+1)*(nm+1)) !(in) 入力スペクトルデータ call spclap(nm,w_data,w_LaplaInv_w,rn(1,2)) end function w_LaplaInv_w function w_DLon_w(w_data) ! ! スペクトルデータに経度微分 ∂/∂λ を作用させる(1 層用). ! ! スペクトルデータの経度微分とは, 対応する格子点データに ! 経度微分∂/∂λを作用させたデータのスペクトル変換のことである. ! real(8) :: w_DLon_w((nm+1)*(nm+1)) !(out) スペクトルデータの経度微分 real(8), intent(in) :: w_data((nm+1)*(nm+1)) !(in) 入力スペクトルデータ call spclam(nm,w_data,w_DLon_w,irm) end function w_DLon_w function xy_GradLon_w(w_data) ! ! スペクトルデータに勾配型経度微分 1/cosφ・∂/∂λ を ! 作用させた格子点データを返す(1 層用). ! real(8) :: xy_GradLon_w(0:im-1,1:jm) !(out) スペクトルデータを勾配型経度微分した格子点データ real(8), intent(in) :: w_data((nm+1)*(nm+1)) !(in) 入力スペクトルデータ xy_GradLon_w = xy_w(w_data,ipow=1,iflag=-1) end function xy_GradLon_w function xy_GradLat_w(w_data) ! ! スペクトルデータに勾配型緯度微分 ∂/∂φ を作用させて ! 格子点データに変換して返す(1 層用). ! real(8) :: xy_GradLat_w(0:im-1,1:jm) !(out) スペクトルデータを勾配型緯度微分した格子点データ real(8), intent(in) :: w_data((nm+1)*(nm+1)) !(in) 入力スペクトルデータ xy_GradLat_w = xy_w(w_data,ipow=1,iflag=1) end function xy_GradLat_w function w_DivLon_xy(xy_data) ! ! 格子点データに発散型経度微分 1/cosφ・∂/∂λ を作用させて ! スペクトルデータに変換して返す(1 層用). ! real(8) :: w_DivLon_xy((nm+1)*(nm+1)) !(out) 格子点データを発散型経度微分したスペクトルデータ real(8), intent(in) :: xy_data(0:im-1,1:jm) !(in) 入力格子点データ w_DivLon_xy = w_xy(xy_data,ipow=1,iflag=-1) end function w_DivLon_xy function w_DivLat_xy(xy_data) ! ! 格子点データに発散型緯度微分 1/cosφ・∂(f cosφ)/∂φ を作用させて ! スペクトルデータに変換して返す(1 層用). ! real(8) :: w_DivLat_xy((nm+1)*(nm+1)) !(out) 格子点データを発散型緯度微分したスペクトルデータ real(8), intent(in) :: xy_data(0:im-1,1:jm) !(in) 入力格子点データ w_DivLat_xy = w_xy(xy_data,ipow=1,iflag=1) end function w_DivLat_xy function w_Div_xy_xy(xy_u,xy_v) ! ! 2 つの入力格子点データをベクトル成分とする発散を計算し, ! スペクトルデータとして返す(1 層用). ! real(8) :: w_Div_xy_xy((nm+1)*(nm+1)) !(out) 2 つの入力格子点データをベクトル成分とする発散のスペクトルデータ real(8), intent(in) :: xy_u(0:im-1,1:jm) !(in) ベクトル経度成分の格子点データ real(8), intent(in) :: xy_v(0:im-1,1:jm) !(in) ベクトル緯度成分の格子点データ w_Div_xy_xy = w_Divlon_xy(xy_u) + w_Divlat_xy(xy_v) end function w_Div_xy_xy function w_Jacobian_w_w(w_a,w_b) ! 2 つのスペクトルデータにヤコビアン ! ! J(f,g) = ∂f/∂λ・∂g/∂μ - ∂g/∂λ・∂f/∂μ ! = ∂f/∂λ・1/cosφ・∂g/∂φ ! - ∂g/∂λ・1/cosφ・∂f/∂φ ! ! を作用させる(1 層用). real(8) :: w_Jacobian_w_w((nm+1)*(nm+1)) !(out) 2 つのスペクトルデータのヤコビアン real(8), intent(in) :: w_a((nm+1)*(nm+1)) !(in) 1つ目の入力スペクトルデータ real(8), intent(in) :: w_b((nm+1)*(nm+1)) !(in) 2つ目の入力スペクトルデータ real(8) :: q(3*((nm+1)/2+nm+1)*jm) ! 作業用配列 real(8) :: ws(iw),ww(iw) ! 作業用配列 call spnjcb(nm,im,im,jm,jm,w_a,w_b,w_Jacobian_w_w,& it,t,y,ip2,p2,r2,ip3,p3,r3,ia,a,q,ws,ww) end function w_Jacobian_w_w !--------------- 微分計算 (λ,μ座標系用) ----------------- function xy_GradLambda_w(w_data) ! ! スペクトルデータに勾配型経度微分 ∂/∂λ を作用する(1 層用). ! real(8) :: xy_GradLambda_w(0:im-1,1:jm) !(out) スペクトルデータを勾配型経度微分した格子点データ real(8), intent(in) :: w_data((nm+1)*(nm+1)) !(in) 入力スペクトルデータ xy_GradLambda_w = xy_w(w_data,ipow=0,iflag=-1) end function xy_GradLambda_w function xy_GradMu_w(w_data) ! ! スペクトルデータに勾配型緯度微分 (1-μ^2)∂/∂μ (μ=sinφ) ! を作用させて格子点データに変換して返す(1 層用). ! real(8) :: xy_GradMu_w(0:im-1,1:jm) !(out) スペクトルデータを勾配型緯度微分した格子点データ real(8), intent(in) :: w_data((nm+1)*(nm+1)) !(in) 入力スペクトルデータ xy_GradMu_w = xy_w(w_data,ipow=0,iflag=1) end function xy_GradMu_w function w_DivLambda_xy(xy_data) ! ! 格子点データに発散型経度微分 1/(1-μ^2)・∂/∂λ (μ=sinφ) ! を作用させてスペクトルデータに変換して返す(1 層用). ! real(8) :: w_DivLambda_xy((nm+1)*(nm+1)) !(out) 格子点データを発散型経度微分したスペクトルデータ real(8), intent(in) :: xy_data(0:im-1,1:jm) !(in) 入力格子点データ w_DivLambda_xy = w_xy(xy_data,ipow=2,iflag=-1) end function w_DivLambda_xy function w_DivMu_xy(xy_data) ! ! 格子点データに発散型緯度微分 ∂/∂μ (μ=sinφ)を作用させて ! スペクトルデータに変換して返す(1 層用). ! real(8) :: w_DivMu_xy((nm+1)*(nm+1)) !(out) 格子点データを発散型緯度微分したスペクトルデータ real(8), intent(in) :: xy_data(0:im-1,1:jm) !(in) 入力格子点データ w_DivMu_xy = w_xy(xy_data,ipow=2,iflag=1) end function w_DivMu_xy !--------------- 終了処理 ----------------- subroutine w_deriv_finalize ! ! モジュールの終了処理(割り付け配列の解放)をおこなう. ! ! このサブルーチンを単独で用いるのでなく, ! 上位サブルーチン w_Finalize を使用すること. ! if ( .not. w_deriv_initialize ) then call MessageNotify('W','w_deriv_Finalize',& 'w_deriv_module not initialized yet') return endif deallocate(rn) ! ラプラシアン演算用配列 deallocate(irm) ! 経度微分演算用配列 deallocate(ip2) ! ヤコビアン計算用配列 deallocate(p2) ! ヤコビアン計算用配列 deallocate(r2) ! ヤコビアン計算用配列 deallocate(ip3) ! ヤコビアン計算用配列 deallocate(p3) ! ヤコビアン計算用配列 deallocate(r3) ! ヤコビアン計算用配列 w_deriv_initialize = .false. call MessageNotify('M','w_deriv_Finalize',& 'w_deriv_module (2013/02/23) is finalized') end subroutine w_deriv_finalize end module w_deriv_module