IGModel-SW 1.0

src/motion_equation.f90

説明を見る。
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
 全て クラス ネームスペース ファイル 関数 変数