!-- !---------------------------------------------------------------------- ! Copyright (c) 2002-2013 SPMDOEL Development Group !---------------------------------------------------------------------- !表題 wa_base_module ! ! spml/wa_base_module モジュールは球面上での流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するための ! モジュール wa_module の下部モジュールであり, スペクトル計算の ! 基本的な Fortran90 関数を提供する. ! ! 球面上の 1 層モデル用 w_base_module モジュールを多層モデル用に ! 拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに ! 対する変換が行える. ! ! 内部で ISPACK の SPPACK と SNPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK/SNPACK,SPPACK のマニュアルを参照されたい. ! ! このモジュールを使うためには前もって w_base_initial を呼んで ! 切断波数, 格子点数の設定をしておく必要がある. ! ! !履歴 2002/02/02 竹広真一 多層用に改造 ! 2002/03/30 竹広真一 モジュール名変更 ! 2002/05/25 竹広真一 格子点座標を区別すべく命名法変更 ! 2005/01/09 竹広真一 msgdmp -> MessageNotify に変更 ! 2005/07/09 竹広真一 OPENMP 版変換ルーチンに対応 ! バンク競合を避けるための作業配列追加 ! 2005/07/10 竹広真一 OpenMP セットアップのメッセージ出力 ! 2006/03/08 竹広真一 コメントを RDoc 用に修正 ! 2007/11/21 竹広真一 初期化サブルーチンメッセージ出力 ! 2008/05/31 竹広真一 初期化サブルーチン分離 ! 2008/06/22 佐々木洋平 水平方向の格子点データ配列の始点を ! 1 から 0 に変更 ! 2008/06/28 佐々木洋平 コメントを RDoc 用に微修正 ! 2008/07/07 竹広真一 配列最終次元の不定性を許容する変更 ! 2008/12/29 竹広真一 xya_wa, wa_xya コメント修正 ! 2009/01/09 竹広真一 wa_base_Initial メッセージに日付を追加 ! 2009/01/29 佐々木洋平 コメントを RDoc 用に微修正 ! 2009/07/30 竹広真一 作業領域をローカル変数に変更(for OpenMP) ! 2013/02/13 竹広真一 wa_StreamPotential2Vector, ! wa_Vector2VorDiv 導入 ! 2013/02/14 竹広真一 w_VectorCosLat2VorDiv 導入 ! 2013/02/23 竹広真一 wa_base_Finalize 導入 ! !++ module wa_base_module ! != wa_base_module ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id$ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要 ! ! spml/wa_base_module モジュールは球面上での流体運動を ! 球面調和函数を用いたスペクトル法によって数値計算するための ! モジュール wa_module の下部モジュールであり, スペクトル計算の ! 基本的な Fortran90 関数を提供する. ! ! 球面上の 1 層モデル用 w_base_module モジュールを多層モデル用に ! 拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに ! 対する変換が行える. ! ! 内部で ISPACK の SPPACK と SNPACK の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に ! ついては ISPACK/SNPACK,SPPACK のマニュアルを参照されたい. ! ! このモジュールを使うためには前もって w_base_initial を呼んで ! 切断波数, 格子点数の設定をしておく必要がある. ! use dc_message use w_base_module, only : im, jm, nm, it, t, y, ip, p, r, ia, a, openmp, np,& w_StreamPotential2Vector, & w_Vector2VorDiv, w_VectorCosLat2VorDiv implicit none integer :: km=16 ! 同時に処理する最大データ数(層の数) integer, allocatable :: ipk(:,:) ! 変換用配列(多層用) real(8), allocatable :: pk(:,:), rk(:,:) ! 変換用配列(多層用) integer :: id=65, jd=33 ! xya_work の大きさ integer :: iw ! ww, ws の大きさ logical :: wa_base_initialize = .false. ! 初期化フラグ real(8), parameter :: pi=3.14159265358979 private public km ! 層数 public wa_base_Initial ! 初期化サブルーチン public wa_base_Finalize ! 終了処理 public xya_wa, wa_xya ! 変換関数 public wa_StreamPotential2Vector ! 流線ポテンシャルから速度場計算 public wa_Vector2VorDiv ! 速度場から渦度発散を計算 public wa_VectorCosLat2VorDiv ! 速度場から渦度発散を計算 save km ! 最大データ数(層数)を記憶 save ipk, pk, rk ! 変換用配列を記憶 save id, jd, iw ! 変換用配列の大きさ save wa_base_initialize ! 初期化フラグ contains !--------------- 初期化 ----------------- subroutine wa_base_initial(k_in) ! ! スペクトル変換の最大データ数(層数)を設定する. ! ! このサブルーチンを単独で用いるのでなく, ! 上位サブルーチン wa_Initial を使用すること. ! integer,intent(in) :: k_in !(in) 最大データ数(層数) km = k_in allocate(ipk(km,((nm+1)/2+nm+1)*2)) ! 変換用配列(多層用) allocate(pk(km,((nm+1)/2+nm+1)*jm)) ! 変換用配列(多層用) allocate(rk(km,((nm+1)/2*2+3)*(nm/2+1))) ! 変換用配列(多層用) if ( im/2*2 .eq. im ) then id = im+1 else id = im endif if ( openmp ) then jd = jm else if ( jm/2*2 .eq. jm ) then jd = jm+1 else jd = jm endif if ( openmp ) then iw=km*(im+nm+1)*3*jm/2 call MessageNotify('M','wa_base_Initial', & 'OpenMP computation was set up.') else iw=km * max((nm+4)*(nm+3),jd*3*(nm+1),jd*im) endif call snkini(nm,jm,km,ip,p,r,ipk,pk,rk) wa_base_initialize = .true. call MessageNotify('M','wa_base_initial',& 'wa_base_module (2013/02/23) is initialized') end subroutine wa_base_Initial !--------------- 基本変換 ----------------- function xya_wa(wa_data,ipow,iflag) ! 球面調和関数スペクトル -> 格子点 ! ! スペクトルデータから格子データへ変換する(多層用). ! real(8), intent(in) :: wa_data(:,:) !(in) スペクトルデータ((nm+1)*(nm+1),:) ! real(8) :: xya_wa(0:im-1,1:jm,size(wa_data,2)) !(out) 格子点データ(0:im-1,1:jm,:) ! integer, intent(in), optional :: ipow !(in) 作用させる 1/cosφ の次数. 省略時は 0. integer, intent(in), optional :: iflag !(in) 変換の種類 ! 0 : 通常の正変換 ! -1 : 経度微分を作用させた正変換 ! 1 : 緯度微分を作用させた正変換 ! 2 : sinφを作用させた正変換 ! 省略時は 0. ! integer, parameter :: ipow_default = 0 integer, parameter :: iflag_default = 0 integer ipval, ifval integer k, i real(8) :: xya_work(id,jd,km) ! 変換用配列 real(8) :: q(km*((nm+1)/2+nm+1)*jm) ! 作業配列(多層用) real(8) :: ws(iw),ww(iw) ! 作業用配列(多層用) real(8) :: wv(km*(nm+4)*(nm+3)*np) ! 作業用配列(OPENMP用) logical :: first=.true. ! 初回判定スイッチ save first if (present(ipow)) then ipval = ipow else ipval = ipow_default endif if (present(iflag)) then ifval = iflag else ifval = iflag_default endif if ( size(wa_data,1) /= (nm+1)**2 ) then call MessageNotify('E','xya_wa','Size of 1st dimension invalid.') end if k=size(wa_data,2) if ( k > km ) then call MessageNotify('E','xya_wa','Size of 2nd dimension invalid.') else if ( openmp ) then if ( first ) then call MessageNotify('M','xya_wa', & 'OpenMP routine SNTSOG/SNPACK is used for spherical harmonic transformation.') endif call sntsog(nm,im,id,jm,k,wa_data,xya_work,& it,t,y,ipk(1:k,:),pk(1:k,:),rk(1:k,:),& ia,a,q,ws,ww,wv,ipval,ifval) else call snts2g(nm,im,id,jm,jd,k,wa_data, xya_work,& it,t,y,ipk(1:k,:),pk(1:k,:),rk(1:k,:),ia,a,q,ws,ww,ipval,ifval) endif do i=0,im-1 xya_wa(i,1:jm,1:k) = xya_work(i+1,1:jm,1:k) enddo first = .false. end function xya_wa function wa_xya(xya_data,ipow,iflag) ! 格子点 -> 球面調和関数スペクトル ! ! 格子データからスペクトルデータへ(正)変換する(多層用). ! real(8), intent(in) :: xya_data(0:,:,:) !(in) 格子点データ(0:im-1,1:jm,:) real(8) :: wa_xya((nm+1)*(nm+1),size(xya_data,3)) !(out) スペクトルデータ((nm+1)*(nm+1),:) integer, intent(in), optional :: ipow !(in) 変換時に同時に作用させる 1/cosφ の次数. 省略時は 0. integer, intent(in), optional :: iflag ! 変換の種類 ! 0 : 通常の正変換 ! -1 : 経度微分を作用させた正変換 ! 1 : 緯度微分を作用させた正変換 ! 2 : sinφを作用させた正変換 ! 省略時は 0. integer, parameter :: ipow_default = 0 ! スイッチデフォルト値 integer, parameter :: iflag_default = 0 ! スイッチデフォルト値 integer ipval, ifval integer i,k real(8) :: xya_work(id,jd,km) ! 変換用配列 real(8) :: q(km*((nm+1)/2+nm+1)*jm) ! 作業配列(多層用) real(8) :: ws(iw),ww(iw) ! 作業用配列(多層用) real(8) :: wv(km*(nm+4)*(nm+3)*np) ! 作業用配列(OPENMP用) logical :: first=.true. ! 初回判定スイッチ save first if (present(ipow)) then ipval = ipow else ipval = ipow_default endif if (present(iflag)) then ifval = iflag else ifval = iflag_default endif if ( size(xya_data,1) /= im ) then call MessageNotify('E','wa_xya','Size of 1st dimension invalid.') endif if ( size(xya_data,2) /= jm ) then call MessageNotify('E','wa_xya','Size of 2nd dimension invalid.') endif k = size(xya_data,3) if ( k > km ) then call MessageNotify('E','wa_xya','Size of 3rd dimension invalid.') endif do i=0,im-1 xya_work(i+1,1:jm,1:k) = xya_data(i,1:jm,1:k) enddo if ( openmp ) then if ( first ) then call MessageNotify('M','wa_xya', & 'OpenMP routine SNTGOS/SNPACK is used for spherical harmonic transformation.') endif call sntgos(nm,im,id,jm,k,xya_work,wa_xya,& it,t,y,ipk(1:k,:),pk(1:k,:),rk(1:k,:),& ia,a,q,ws,ww,wv,ipval,ifval) else call sntg2s(nm,im,id,jm,jd,k,xya_work,wa_xya,& it,t,y,ipk(1:k,:),pk(1:k,:),rk(1:k,:),ia,a,q,ws,ww,ipval,ifval) endif first = .false. end function wa_xya !----------- 速度, 渦度・発散, 流線・ポテンシャル計算 ------------- subroutine wa_StreamPotential2Vector(wa_Psi, wa_Chi, xya_U, xya_V) ! ! 流線・ポテンシャル(スペクトルデータ)から速度場(格子データ)に ! (逆)変換する(多層用) ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. ! ! u cosφ = ∂χ/∂λ - cosφ∂ψ/∂φ, ! v cosφ = cosφ∂χ/∂φ + ∂ψ/∂λ ! real(8), intent(in) :: wa_Psi(:,:) !(in) 流線関数((nm+1)*(nm+1),:) real(8), intent(in) :: wa_Chi(:,:) !(in) 速度ポテンシャル((nm+1)*(nm+1),:) real(8), intent(out) :: xya_U(0:im-1,1:jm,size(wa_Psi,2)) !(out) 速度経度成分(0:im-1,1:jm,:) real(8), intent(out) :: xya_V(0:im-1,1:jm,size(wa_Psi,2)) !(out) 速度緯度成分(0:im-1,1:jm,:) integer :: k do k=1,size(wa_Psi,2) call w_StreamPotential2Vector & ( wa_Psi(:,k), wa_Chi(:,k), xya_U(:,:,k), xya_V(:,:,k) ) enddo end subroutine wa_StreamPotential2Vector subroutine wa_Vector2VorDiv(xya_U, xya_V, wa_Vor, wa_Div) ! ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に ! (正)変換する(多層用) ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. ! ! ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ ! D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ ! real(8), intent(in) :: xya_U(0:,:,:) !(in) 速度経度成分(0:im-1,1:jm,:) real(8), intent(in) :: xya_V(0:,:,:) !(in) 速度緯度成分(0:im-1,1:jm,:) real(8), intent(out) :: wa_Vor((nm+1)*(nm+1),size(xya_U,3)) !(out) 流線関数 real(8), intent(out) :: wa_Div((nm+1)*(nm+1),size(xya_U,3)) !(out) 速度ポテンシャル integer :: k do k=1,size(xya_U,3) call w_Vector2VorDiv & (xya_U(:,:,k), xya_V(:,:,k), wa_Vor(:,k), wa_Div(:,k)) enddo end subroutine wa_Vector2VorDiv subroutine wa_VectorCosLat2VorDiv(xya_UCosLat, xya_VCosLat, wa_Vor, wa_Div) ! ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に ! (正)変換する(多層用) ! ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ. ! ! ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ ! D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ ! real(8), intent(in) :: xya_UCosLat(0:,:,:) !(in) 速度経度成分 * cos(lat) (0:im-1,1:jm,:) real(8), intent(in) :: xya_VCosLat(0:,:,:) !(in) 速度緯度成分 * cos(lat) (0:im-1,1:jm,:) real(8), intent(out) :: wa_Vor((nm+1)*(nm+1),size(xya_UCosLat,3)) !(out) 流線関数 real(8), intent(out) :: wa_Div((nm+1)*(nm+1),size(xya_UCosLat,3)) !(out) 速度ポテンシャル integer :: k do k=1,size(xya_UCosLat,3) call w_VectorCosLat2VorDiv & (xya_UCosLat(:,:,k), xya_VCosLat(:,:,k), wa_Vor(:,k), wa_Div(:,k)) enddo end subroutine wa_VectorCosLat2VorDiv !--------------- 終了処理 ----------------- subroutine wa_base_Finalize ! ! モジュールの終了処理(割り付け配列の解放)をおこなう. ! ! このサブルーチンを単独で用いるのでなく, ! 上位サブルーチン wa_Finalize を使用すること. ! if ( .not. wa_base_initialize ) then call MessageNotify('W','wa_base_Finalize',& 'wa_base_module not initialized yet') return endif deallocate(ipk) ! 変換用配列(多層用) deallocate(pk) ! 変換用配列(多層用) deallocate(rk) ! 変換用配列(多層用) call MessageNotify('M','wa_base_finalize',& 'wa_base_module (2013/02/23) is finalized') end subroutine wa_base_Finalize end module wa_base_module