! -*- mode: f90; coding: utf-8 -*- !------------------------------------------------------------------------- ! Copyright (c) 2019 SPMODEL Development Group. All rights reserved. ! ! 履歴 2019/03/16 竹広真一 ! 2019/10/04 佐々木洋平 ! 2019/10/09 佐々木洋平 ! 2020/11/20 竹広真一 !------------------------------------------------------------------------- !> @author Shin-ichi Takehiro, Youhei SASAKI !> @copyright Copyright (C) SPMODEL Development Group, 2019. !> License is MIT/X11. see [COPYRIGHT](@ref COPYRIGHT) in detail !> @brief @ref w_module の下部モジュール: 微分 !> @warning !> 本モジュールを直接呼ぶ事は想定されていないことに注意されたい． !> ユーザは @ref w_module を呼ぶこと． !> module w_module_deriv_mint use dc_message, only : MessageNotify use dc_types, only: DP use w_module_base_mint, only : & & im, jm, mm, nn, mi, & & w_base_Initial, xy_w, w_xy implicit none private public w_deriv_Initial public w_deriv_Finalize public w_Lapla_w public w_LaplaInv_w public w_DLon_w public xy_GradLon_w public xy_GradLat_w public w_DivLon_xy public w_DivLat_xy public w_Div_xy_xy public w_Jacobian_w_w public xy_GradLambda_w public xy_GradMu_w public w_DivLambda_xy public w_DivMu_xy public rn !--------------------------------------------------------------------- !> ラプラシアン演算用配列(後方互換性のため) !> !> * スペクトルデータのラプラシアンを計算するための係数 !> 配列のサイズは((nn+1)*(nn+1), 2) !> * r(L,1) には L 番目の格納位置のスペクトルに対するラプラシアン計算の !> 係数 -n(n+1) の値が格納されている. !> real(DP), allocatable :: rn(:,:) !--------------------------------------------------------------------- !> ラプラシアン演算用配列 !> !> * スペクトルデータのラプラシアンを計算するための係数 !> 配列のサイズは((mm+1)*(mm+1)*2) !> real(DP), allocatable :: D(:) !> 初期化フラグ logical :: w_deriv_initialize=.false. save D, rn contains !---------------------------------------------------------- !> スペクトル微分計算に必要となる作業領域を設定する. !> !> 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで !> 初期設定をしなければならない. !> !> @warning !> このサブルーチンを単独で用いるのでなく, !> 上位サブルーチン !> [w_Initial](@ref w_module#w_initial) を使用すること. !> subroutine w_deriv_initial allocate(D((MM/MI+1)*(2*(NN+1)-MM/MI*MI)*2_8)) allocate(rn((MM/MI+1)*(2*(NN+1)-MM/MI*MI),2_8)) call sxinmd(mm,nn,D,mi) rn = reshape(D,(/(MM/MI+1)*(2*(NN+1)-MM/MI*MI),2_8/)) w_deriv_initialize=.true. call MessageNotify('M','w_deriv_initial',& 'w_module_deriv_mint (2020/11/18) is initialized') end subroutine w_deriv_initial !---------------------------------------------------------- !> 入力スペクトルデータにラプラシアン !> \f{align*}{ !> ▽^2 = 1/\cos^2φ・∂^2/∂λ^2 + 1/\cosφ・∂/∂φ(\cosφ∂/∂φ) !> \f} !> を作用する(1 層用). !> function w_Lapla_w(w_data) !> 入力スペクトルデータ real(DP), intent(in) :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) !> 入力スペクトルデータのラプラシアン real(DP) :: w_Lapla_w((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) call sxclam(mm,nn,w_data,w_Lapla_w,D,1_8,mi) end function w_Lapla_w !---------------------------------------------------------- !> 入力スペクトルデータに逆ラプラシアン !> \f{align*}{ !> ▽^{-2} !> =[1/\cos^2φ・∂^2/∂λ^2 + 1/\cosφ・∂/∂φ(\cosφ∂/∂φ)]^{-1} !> \f} !> を作用する(1 層用). !> function w_LaplaInv_w(w_data) !> 入力スペクトルデータ real(DP), intent(in) :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) !> スペクトルデータの逆ラプラシアン real(DP) :: w_LaplaInv_w((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) call sxclam(mm,nn,w_data,w_LaplaInv_w,D,2_8,mi) end function w_LaplaInv_w !---------------------------------------------------------- !> スペクトルデータに経度微分 \f$∂/∂λ\f$ を作用させる(1 層用). !> function w_DLon_w(w_data) !> 入力スペクトルデータ real(DP), intent(in) :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) !> スペクトルデータの経度微分 real(DP) :: w_DLon_w((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) call sxcsmx(mm,nn,w_data,w_DLon_w,mi) end function w_DLon_w !------------------------------------------------------------- !> スペクトルデータに勾配型経度微分 \f$1/\cosφ・∂/∂λ\f$ を !> 作用させた格子点データを返す(1 層用). !> function xy_GradLon_w(w_data) !> 入力スペクトルデータ real(DP), intent(in) :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) !> スペクトルデータを勾配型経度微分した格子点データ real(DP) :: xy_GradLon_w(0:im-1,1:jm) xy_GradLon_w = xy_w(w_data,ipow=1,iflag=-1) end function xy_GradLon_w !--------------------------------------------------------------- !> スペクトルデータに勾配型緯度微分 \f$∂/∂φ\f$ を作用させて !> 格子点データに変換して返す(1 層用). !> function xy_GradLat_w(w_data) !> 入力スペクトルデータ real(DP), intent(in) :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) !> スペクトルデータを勾配型緯度微分した格子点データ real(DP) :: xy_GradLat_w(0:im-1,1:jm) xy_GradLat_w = xy_w(w_data,ipow=1,iflag=1) end function xy_GradLat_w !--------------------------------------------------------------------- !> 格子点データに発散型経度微分 \f$1/\cosφ・∂/∂λ\f$ を作用させて !> スペクトルデータに変換して返す(1 層用). !> function w_DivLon_xy(xy_data) !> 入力格子点データ real(DP), intent(in) :: xy_data(0:im-1,1:jm) !> 格子点データを発散型経度微分したスペクトルデータ real(DP) :: w_DivLon_xy((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) w_DivLon_xy = w_xy(xy_data,ipow=1,iflag=-1) end function w_DivLon_xy !------------------------------------------------------------------- !> 格子点データに発散型緯度微分 \f$1/\cosφ・∂(f \cosφ)/∂φ\f$ !> を作用させてスペクトルデータに変換して返す(1 層用). !> function w_DivLat_xy(xy_data) !> 入力格子点データ real(DP), intent(in) :: xy_data(0:im-1,1:jm) !> 格子点データを発散型緯度微分したスペクトルデータ real(DP) :: w_DivLat_xy((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) w_DivLat_xy = w_xy(xy_data,ipow=1,iflag=1) end function w_DivLat_xy !---------------------------------------------------------- !> 2 つの入力格子点データをベクトル成分とする発散を計算し, !> スペクトルデータとして返す(1 層用). !> function w_Div_xy_xy(xy_u,xy_v) !> ベクトル経度成分の格子点データ real(DP), intent(in) :: xy_u(0:im-1,1:jm) !> ベクトル緯度成分の格子点データ real(DP), intent(in) :: xy_v(0:im-1,1:jm) !> 2 つの入力格子点データをベクトル成分とする発散のスペクトルデータ real(DP) :: w_Div_xy_xy((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) w_Div_xy_xy = w_Divlon_xy(xy_u) + w_Divlat_xy(xy_v) end function w_Div_xy_xy !------------------------------------------------------------------------- !> 2 つのスペクトルデータにヤコビアン !> \f{align*}{ !> J(f,g) & = ∂f/∂λ・∂g/∂μ - ∂g/∂λ・∂f/∂μ \\ !> & = ∂f/∂λ・1/\cos φ・∂g/∂φ !> - ∂g/∂λ・1/\cos φ・∂f/∂φ !> \f} !> を作用させる(1 層用). !> function w_Jacobian_w_w(w_a,w_b) !> 1つ目の入力スペクトルデータ real(DP), intent(in) :: w_a((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) !> 2つ目の入力スペクトルデータ real(DP), intent(in) :: w_b((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) !> 2 つのスペクトルデータのヤコビアン real(DP) :: w_Jacobian_w_w((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) w_Jacobian_w_w = & & w_xy( xy_w(w_DLon_w(w_a))*xy_w(w_b,ipow=2,iflag=1) & & - xy_w(w_DLon_w(w_b))*xy_w(w_a,ipow=2,iflag=1) ) end function w_Jacobian_w_w !---------------------------------------------------------- !> スペクトルデータに勾配型経度微分 ∂/∂λ を作用する(1 層用). !> function xy_GradLambda_w(w_data) !> 入力スペクトルデータ real(DP), intent(in) :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) !> スペクトルデータを勾配型経度微分した格子点データ real(DP) :: xy_GradLambda_w(0:im-1,1:jm) xy_GradLambda_w = xy_w(w_data,ipow=0,iflag=-1) end function xy_GradLambda_w !---------------------------------------------------------- !> スペクトルデータに勾配型緯度微分 !> \f$(1-μ^2)∂/∂μ ,\, (μ=\sinφ)\f$ !> を作用させて格子点データに変換して返す(1 層用). !> function xy_GradMu_w(w_data) !> 入力スペクトルデータ real(DP), intent(in) :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) !> スペクトルデータを勾配型緯度微分した格子点データ real(DP) :: xy_GradMu_w(0:im-1,1:jm) xy_GradMu_w = xy_w(w_data,ipow=0,iflag=1) end function xy_GradMu_w !---------------------------------------------------------- !> 格子点データに発散型経度微分 !> \f$1/(1-μ^2)・∂/∂λ, \, (μ=\sinφ)\f$ !> を作用させてスペクトルデータに変換して返す(1 層用). !> function w_DivLambda_xy(xy_data) !> 入力格子点データ real(DP), intent(in) :: xy_data(0:im-1,1:jm) !> 格子点データを発散型経度微分したスペクトルデータ real(DP) :: w_DivLambda_xy((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) w_DivLambda_xy = w_xy(xy_data,ipow=2,iflag=-1) end function w_DivLambda_xy !---------------------------------------------------------- !> 格子点データに発散型緯度微分 !> \f$∂/∂μ, \, (μ=\sinφ)\f$ !> を作用させてスペクトルデータに変換して返す(1 層用). !> function w_DivMu_xy(xy_data) !> 入力格子点データ real(DP), intent(in) :: xy_data(0:im-1,1:jm) !> 格子点データを発散型緯度微分したスペクトルデータ real(DP) :: w_DivMu_xy((mm/mi+1)*(2*(nn+1)-mm/mi*mi)) w_DivMu_xy = w_xy(xy_data,ipow=2,iflag=1) end function w_DivMu_xy !---------------------------------------------------------- !> モジュールの終了処理(割り付け配列の解放)をおこなう. !> !> @warning !> このサブルーチンを単独で用いるのでなく, !> 上位サブルーチン !> [w_Finalize](@ref w_module#w_finalize) を使用すること. !> subroutine w_deriv_finalize if ( .not. w_deriv_initialize ) then call MessageNotify('W','w_deriv_Finalize',& 'w_module_deriv_mint not initialized yet') return endif deallocate(D) ! ラプラシアン演算用配列 deallocate(rn) ! ラプラシアン演算用配列 w_deriv_initialize = .false. call MessageNotify('M','w_deriv_finalize',& 'w_module_deriv_mint (2020/11/18) is finalized') end subroutine w_deriv_finalize end module w_module_deriv_mint