IGModel-SW 1.0
|
00001 00012 module motion_equation 00013 00014 ! モジュール引用 ; Use statements 00015 ! 00016 00017 !* gtool **** 00018 !* 00019 00020 ! 種類型パラメータ 00021 ! Kind type parameter 00022 ! 00023 use dc_types, only: DP ! 倍精度実数型. Double precision. 00024 00025 !* IGMBaseLib **** 00026 !* 00027 00028 ! 数学 00029 ! Mathematica 00030 ! 00031 use igmcore_math, only: PI 00032 00033 ! 線形代数 00034 ! Linear algebra 00035 ! 00036 use igmcore_linear_algebra, only: & 00037 & dot, cross 00038 00039 ! 座標変換 00040 ! Coordinate conversion 00041 ! 00042 use igmcore_coordinate_conversion, only: & 00043 & orth_to_geo_pos, geo_to_orth_vec 00044 00045 00046 ! 正二十面体格子上の物理場の管理 00047 ! Manager of a physical field on icosahedral grid 00048 ! 00049 use Field_IcGrid2D_Manager, only: & 00050 & Field_IcGrid2D, & 00051 & Field_IcGrid2D_Init, & 00052 & get_icgrid, paste_field2D_margin 00053 00054 ! 正二十面体格子の管理 00055 ! Manager of icoshedral grid 00056 ! 00057 use IcGrid2D_FVM_Manager, only: & 00058 & IcGrid2D_FVM, & 00059 & get_EffSize_Min, get_EffSize_Max, check_pole, & 00060 & RC_REGIONS_NUM, NOT_POLE_FLAG 00061 00062 ! 正二十面体格子上の物理場に作用する微分演算を提供するクラス 00063 ! Class providing some differential operators acting on a physcial field on icosahedral grid. 00064 ! 00065 use Derivate_Field_IcGrid2D_Manager, only: & 00066 & Derivate_Field_IcGrid2D, & 00067 & vertical_curl_op, gradient_op 00068 00069 !* IGModel-SW **** 00070 !* 00071 00072 ! シミュレーションパラメータの管理 00073 ! Manager of simulation parameter 00074 ! 00075 use param_manager, only: & 00076 & Grav, Omega, alpha 00077 00078 ! OpenMP 00079 ! 00080 ! 00081 use omp_lib 00082 00083 ! 宣言文 ; Declaration statements 00084 ! 00085 implicit none 00086 private 00087 00088 ! 公開手続き 00089 ! Public statements 00090 ! 00091 public :: motion_equation_Init, calc_motion_eq_dvdt 00092 00093 ! 非公開変数 00094 ! Private statements 00095 ! 00096 00099 type(Field_IcGrid2D), save :: xy_Zeta 00100 00103 type(Field_IcGrid2D), save :: mecha_energy 00104 00107 type(Field_IcGrid2D), save :: grad_mecha_energy 00108 00109 contains 00110 00111 ! 00122 subroutine motion_equation_Init( & 00123 & icgrid & ! (inout) 00124 & ) 00125 00126 ! 宣言文 ; Declaration statements 00127 ! 00128 type(IcGrid2D_FVM), intent(inout) :: icgrid 00129 00130 ! 実行文 ; Executable statement 00131 ! 00132 00133 ! 渦度場, 運動エネルギーの場, 運動エネルギーの勾配場を管理する, 構造型 Field_IcGrid2D の変数を初期化する. 00134 ! 00135 00136 call Field_IcGrid2D_Init( & 00137 & self=xy_Zeta, & 00138 & icgrid=icgrid, name='relative vorticity field', rank=1 ) 00139 00140 call Field_IcGrid2D_Init( & 00141 & self=mecha_energy, & 00142 & icgrid=icgrid, name='mechanical_energy_field', rank=1 ) 00143 00144 call Field_IcGrid2D_Init( & 00145 & self=grad_mecha_energy, & 00146 & icgrid=icgrid, name='gradient of mechanical energy field', rank=3 ) 00147 00148 end subroutine motion_equation_Init 00149 00150 ! 00186 subroutine calc_motion_eq_dvdt( & 00187 & DVelDtN, & ! (inout) 00188 & xy_VelN, xy_HeightN, xy_Coli, & ! (in) 00189 & diff_eval & ! (inout) 00190 & ) 00191 00192 ! 宣言文 ; Declaration statements 00193 ! 00194 type(Field_IcGrid2D), intent(inout) :: DVelDtN 00195 type(Field_IcGrid2D), intent(in) :: xy_VelN 00196 type(Field_IcGrid2D), intent(in) :: xy_HeightN 00197 type(Field_IcGrid2D), intent(in) :: xy_Coli 00198 type(Derivate_Field_IcGrid2D), intent(inout) :: diff_eval 00199 00200 00201 ! 作業変数 00202 ! Work variables 00203 ! 00204 integer :: EMax, EMin 00205 integer :: rcID,i,j 00206 real(DP) :: s_pos(3) 00207 real(DP) :: unit_k(3) 00208 type(IcGrid2D_FVM), pointer :: icgrid 00209 00210 ! 実行文 ; Executable statements 00211 ! 00212 00213 icgrid => get_icgrid(xy_VelN) 00214 EMin = get_EffSize_Min(icgrid) 00215 EMax = get_EffSize_Max(icgrid) 00216 00217 00218 ! 速度場に回転演算子を作用させて, 相対渦度場を計算する. 00219 ! Curl operator acts on the velocity field to calculate the relative vorticity field. 00220 call vertical_curl_op(diff_eval, xy_VelN, xy_Zeta ) 00221 00222 ! ジオポテンシャルと運動エネルギーの和の空間勾配場を計算する. 00223 ! Calculates the spatial gradient geopotential and kinetic energy. 00224 ! 00225 00226 ! まず, 各格子点において力学的エネルギーを計算する. 00227 ! First, calculates the mechanical energy at each grid point. 00228 do rcID=1, RC_REGIONS_NUM 00229 !$omp parallel do 00230 do j=EMin, EMax 00231 do i=EMin, EMax 00232 mecha_energy%val(rcID,i,j,:) = & 00233 & Grav * xy_HeightN%val(rcID,i,j,:) & 00234 & + 0.5d0 * dot( xy_VelN%val(rcID,i,j,:), xy_VelN%val(rcID,i,j,:) ) 00235 end do 00236 end do 00237 end do 00238 00239 ! 袖領域の貼りあわせ 00240 call paste_field2D_margin(mecha_energy) 00241 00242 ! 力学的エネルギーの場に, 勾配演算子を作用させる. 00243 ! Gradient operator acts on the mechanical energy field. 00244 call gradient_op(diff_eval, mecha_energy, grad_mecha_energy) 00245 00246 00247 ! 運動方程式の半離散式の右辺を評価し, タイムレベル n における速度変化率を計算する. 00248 ! Evaluates the right-hand side of the semi-discreted equation of motion, 00249 ! and calculates the rate of change of velocity at time level n with respect to time. 00250 ! 00251 00252 do rcID=1, RC_REGIONS_NUM 00253 !$omp parallel do private(s_pos, unit_k) 00254 do j=EMin, EMax 00255 do i=EMin, EMax 00256 00257 ! 格子点位置を極座標として取得する. 00258 s_pos(:) = orth_to_geo_pos( icgrid%rcs_AGrid(rcID,i,j,:) ) 00259 00260 ! 格子点上におけるコリオリパラメータと鉛直方向の単位ベクトル(デカルト座標系表記)を求める. 00261 if( check_pole(icgrid, rcID,i,j) == NOT_POLE_FLAG ) then 00262 unit_k(:) = local_k(s_pos(:)) 00263 else 00264 ! 両極 00265 if ( rcID <= 5 ) then 00266 unit_k(:) = (/ 0.0d0, 0.0d0, 1.0d0 /) 00267 s_pos(2) = PI / 2.0d0 00268 else 00269 unit_k(:) = (/ 0.0d0, 0.0d0, -1.0d0 /) 00270 s_pos(2) = - PI / 2.0d0 00271 end if 00272 end if 00273 00274 ! 速度変化率を計算する. 00275 DVelDtN%val(rcID,i,j,:) = & 00276 & - ( xy_Zeta%val(rcID,i,j,1) + xy_Coli%val(rcID,i,j,1) ) * cross( unit_k(:), xy_VelN%val(rcID,i,j,:) ) & 00277 & - grad_mecha_energy%val(rcID,i,j,:) 00278 00279 end do 00280 end do 00281 end do 00282 00283 end subroutine calc_motion_eq_dvdt 00284 00285 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00286 ! 00287 ! 非公開手続き 00288 ! Private procedure 00289 ! 00290 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00291 00292 ! 00305 function local_k( & 00306 & s_pos & ! (in) 00307 & ) result(vec) 00308 00309 ! 宣言文 ; Declaration statements 00310 ! 00311 real(DP), intent(in) :: s_pos(3) 00312 real(DP) :: vec(3) 00313 00314 ! 作業変数 00315 ! Work variables 00316 ! 00317 00318 ! 地理座標系における鉛直方向の単位ベクトル. 00319 real(DP), parameter :: s_k(3) = (/ 0.0d0, 0.0d0, 1.0d0 /) 00320 00321 ! 実行文 00322 ! 00323 00324 ! 地理座標系における鉛直方向の単位ベクトルを, 直交座標系へと変換する. 00325 vec = geo_to_orth_vec( s_k, s_pos) 00326 00327 end function local_k 00328 00329 end module motion_equation