IGModel-SW 1.0

src/Williamson1992_TestCase/class_TestCase2.f90

説明を見る。
00001 
00009 module class_TestCase2
00010   
00011   ! モジュール引用 ; Use statements
00012   !
00013 
00014   ! 種類型パラメータ
00015   ! Kind type parameter 
00016   !
00017   use dc_types, only: DP,  & ! 倍精度実数型. Double precision. 
00018     &                 TOKEN
00019 
00020   !
00021   !
00022   !
00023   use gtool_history
00024 
00025   ! * IGMBaseLib *
00026   ! *
00027 
00028   ! 数学
00029   ! Mathematic
00030   !
00031   use igmcore_math, only: PI
00032 
00033   ! 座標変換
00034   ! Coordinate conversion
00035   !
00036   use igmcore_coordinate_conversion, only: &
00037     & orth_to_geo_pos, geo_to_orth_pos
00038 
00039   ! 正二十面体格子の管理クラス
00040   !
00041   !
00042   use IcGrid2D_FVM_Manager, only: &
00043     & IcGrid2D_FVM, &
00044     & get_glevel, get_IcRadius, get_EffSize_Min, get_EffSize_Max, check_pole, &
00045     & RC_REGIONS_NUM, NOT_POLE_FLAG
00046   
00047   ! 正二十面体格子上に分布する物理場の管理クラス
00048   !
00049   !
00050   use Field_IcGrid2D_manager, only: &
00051     & Field_IcGrid2D, &
00052     & Field_IcGrid2D_Init, get_icgrid
00053   
00054   !
00055   !
00056   !
00057   use Field_Pattern_Builder, only: &
00058     & create_solid_rotation_field
00059 
00060   ! 正二十面体格子上の物理場における誤差ノルムチェック
00061   !
00062   !
00063   use Field_IcGrid2D_error_norm_check, only: &
00064     & numerical_error_norm1, numerical_error_norm2, &
00065     & numerical_error_norminfinity 
00066 
00067 
00068   !
00069   !
00070   !
00071   use IcGrid_ncWriter_mod, only: &
00072     & IcGrid_ncWriter, IcGrid_ncWriter_Init, &
00073     & open_ncFile, ncdef_simulation_parameter, ncdef_GridData, ncdef_FieldData, end_ncdef, &
00074     & write_GridData, write_FieldData, close_ncFile, &
00075     & increase_recorde_counter
00076 
00077 
00078   ! * IGModel-SW ****
00079   ! *
00080 
00081 
00082   ! パラメータ管理
00083   !
00084   !
00085   use param_manager, only: &
00086     & alpha, output_tick, Grav, Omega, &
00087     & integration_time, DelTime
00088 
00089 
00090   ! 宣言文 ; Declaration statement
00091   !
00092   implicit none
00093   private
00094   
00095   ! 公開変数
00096   ! Public variables
00097   !
00098   
00101   real(DP) :: h_0
00102 
00105   real(DP) :: u_0
00106 
00109   real(DP) :: angular_speed
00110 
00111   ! クラス定義
00112   ! Class definition
00113   !
00114 
00115   !
00119   type, public :: TestCase2
00120 
00121     ! 宣言文 ; Declaration statement
00122     !
00123 
00124     ! 非公開メンバ変数
00125     ! Private member variable
00126     !
00127 
00129     type(Field_IcGrid2D) :: true_h
00130     ! f2003: type(Field_IcGrid2D), private :: true_h
00131 
00132     
00134     type(Field_IcGrid2D) :: error_h
00135     ! f2003: type(Field_IcGrid2D), private :: error_h
00136     
00138     type(Field_IcGrid2D) :: true_v
00139     ! f2003: type(Field_IcGrid2D), private :: true_v
00140 
00141     
00143     integer :: errorh_ncVarID
00144     ! f2003: integer, private :: errorh_ncVarID
00145 
00146   end type TestCase2
00147 
00148   !
00149   !
00150   !
00151   interface initialize_TestCase
00152     module procedure init_TestCase2
00153   end interface initialize_TestCase
00154 
00155   interface finalize_TestCase
00156     module procedure finalize_TestCase2
00157   end interface finalize_TestCase
00158 
00159   public :: initialize_TestCase, finalize_TestCase
00160   public :: set_initial_h, set_initial_hs, set_initial_v, timelevel_Updated
00161 
00162 
00163   ! 非公開変数
00164   ! Private variables
00165   !
00166   
00168   character(TOKEN) :: filename = 'error_norm.dat'
00169   
00171   type(IcGrid_ncWriter), save :: writer
00172 
00173 contains
00174 
00175 !
00181 subroutine init_TestCase2(self, icgrid_ref )
00182   ! 宣言文 ; Declaration statement
00183   !
00184   type(TestCase2), intent(inout) :: self
00185   type(IcGrid2D_FVM), intent(in) :: icgrid_ref
00186 
00187   ! 作業変数
00188   ! Work variables
00189   !
00190   real(DP) :: ic_radius
00191 
00192   ! 実行文 ; Executable statement
00193   !
00194 
00195   ic_radius = get_IcRadius(icgrid_ref)
00196   
00197  
00198   u_0 = 2.0d0 * PI * ic_radius / ( 12.0d0 * 24.0d0 * 3600d0 )
00199   h_0 = 2.94e04 / Grav
00200   angular_speed = 2.0d0 * PI / ( 12.0d0 * 24.0d0 * 3600d0 )
00201 
00202   !
00203   call Field_IcGrid2D_Init(self%true_h, icgrid_ref, 'true_h', 1)
00204   call Field_IcGrid2D_Init(self%error_h, icgrid_ref, 'error_h', 1)
00205 
00206   !
00207   call Field_IcGrid2D_Init(self%true_v, icgrid_ref, 'true_v', 3)
00208 
00209   ! 数値誤差ノルムを書きだすための, NetCDF ファイルを用意する.  
00210   
00211   call HistoryCreate( &
00212     & file='error_norm.nc', title='error norm check', &
00213     & source=' ', institution=' ', &
00214     & dims=(/ 't' /), dimsizes=(/ 0 /), &
00215     & longnames=(/ 'time' /), units=(/ 's' /), &
00216     & origin=real(0d0), interval=real(output_tick) )
00217   
00218   call HistoryAddVariable( &
00219     & varname='norm1', dims=(/ 't' /), &
00220     & longname='norm1', units=' ', xtype='double' )
00221 
00222   call HistoryAddVariable( &
00223     & varname='norm2', dims=(/ 't' /), &
00224     & longname='norm2', units=' ', xtype='double' )
00225 
00226   call HistoryAddVariable( &
00227     & varname='norm_infinity', dims=(/ 't' /), &
00228     & longname='norm_infinity', units=' ', xtype='double' )
00229 
00230   call HistoryAddAttr('norm1', '+glevel', get_glevel(icgrid_ref))
00231   call HistoryAddAttr('norm1', '+alpha', alpha)
00232   
00233   !
00234   !
00235   call IcGrid_ncWriter_Init(writer, icgrid_ref)
00236   call open_ncFile(writer, 'h_error_check.nc')
00237   call ncdef_Simulation_Parameter(writer, integration_time, DelTime, output_tick)
00238 
00239   ! 格子点情報の定義 
00240   call ncdef_GridData(writer)
00241   self%errorh_ncVarID = ncdef_FieldData(writer, self%error_h)
00242 
00243   ! 定義の終了. 
00244   call end_ncdef(writer)
00245 
00246   ! 格子情報を netCDF に書きこむ. 
00247   !
00248   call write_GridData(writer)
00249 
00250 end subroutine init_TestCase2
00251 
00252 !
00257 subroutine finalize_TestCase2(self)
00258   ! 宣言文 ; Declaration statement
00259   !
00260   type(TestCase2), intent(inout) :: self
00261 
00262   ! 実行文 ; Executable statement
00263   !
00264 
00265   !
00266   call HistoryClose()
00267   call close_ncFile(writer)
00268 
00269 end subroutine finalize_TestCase2
00270 
00271 !
00277 subroutine set_initial_v(self, init_v)
00278   ! 宣言文 ; Declaration statement
00279   !
00280   type(TestCase2), intent(inout) :: self
00281   type(Field_IcGrid2D), intent(inout) :: init_v
00282 
00283   ! 実行文 ; Executable statement
00284   !
00285 
00286   write(*,*) 'init_Testcase2_v'
00287 
00288   call create_solid_rotation_field(self%true_v, angular_speed, alpha)
00289   init_v%val(:,:,:,:) =  self%true_v%val(:,:,:,:)
00290 
00291 end subroutine set_initial_v
00292 
00293 !
00299 subroutine set_initial_h(self, init_h)
00300   ! 宣言文 ; Declaration statement
00301   !
00302   type(TestCase2), intent(inout) :: self
00303   type(Field_IcGrid2D), intent(inout) :: init_h
00304 
00305   !
00306   ! Work Variables
00307   !
00308   integer :: rcID, i,j
00309   integer :: EMin, EMax
00310   type(IcGrid2D_FVM), pointer :: icgrid
00311   real(DP) :: s_pos(3)
00312   real(DP) :: ic_radius
00313   real(DP) :: tmp
00314 
00315   ! 実行文 ; Executable statement
00316   !
00317   
00318   write(*,*) 'init_Testcase2_h'
00319 
00320   icgrid => get_IcGrid(init_h)
00321   ic_radius = get_IcRadius(icgrid)
00322   EMin = get_EffSize_Min(icgrid)
00323   EMax = get_EffSize_Max(icgrid)
00324   
00325   do rcID= 1, RC_REGIONS_NUM
00326     do j=EMin, EMax
00327       do i=EMin, EMax
00328 
00329         s_pos(:) = orth_to_geo_pos(icgrid%rcs_AGrid(rcID,i,j,:))
00330         if( check_pole(icgrid, rcID,i,j) == NOT_POLE_FLAG ) then
00331           !
00332           tmp = dcos(s_pos(1)) * dcos(s_pos(2)) * dsin(alpha) + dsin(s_pos(2)) * dcos(alpha)
00333         else
00334           !
00335           tmp = dsin(s_pos(2)) * dcos(alpha)
00336 
00337         end if
00338         
00339         self%true_h%val(rcID,i,j, :) = &
00340           & h_0 - ( ic_radius * Omega * u_0 + 0.5d0 * u_0**2 ) * tmp**2 / Grav       
00341 
00342      end do
00343     end do
00344   end do
00345 
00346   init_h%val(:,:,:,:) = self%true_h%val(:,:,:,:)
00347 
00348 end subroutine set_initial_h
00349 
00350 !
00356 subroutine set_initial_hs(self, init_hs)
00357 
00358   ! 宣言文 ; Declaration statement
00359   !
00360   type(TestCase2), intent(inout) :: self
00361   type(Field_IcGrid2D), intent(inout) :: init_hs
00362 
00363   ! 実行文 ; Executable statement
00364   !
00365 
00366   write(*,*) 'Testcase2: initialize height surface field'  
00367   init_hs%val(:,:,:,:) = 0.0d0
00368 
00369 end subroutine set_initial_hs
00370 
00371 !
00383 subroutine timelevel_Updated(self, tstep, dt, v_n, h_n)
00384 
00385   ! 宣言文 ; Declaration statement
00386   !
00387   type(TestCase2), intent(inout) :: self
00388   integer, intent(in) :: tstep
00389   real(DP), intent(in) :: dt
00390   type(Field_IcGrid2D), intent(inout) :: v_n
00391   type(Field_IcGrid2D), intent(inout) :: h_n
00392 
00393 
00394   ! 実行文 ; Executable statement
00395   !
00396   
00397   if ( mod(tstep * dt, output_tick) == 0 ) then
00398   !
00399     call eval_numcal_h_solution( tstep*dt, self%true_h, h_n)
00400     !
00401     self%error_h%val(:,:,:,:) = self%true_h%val(:,:,:,:) - h_n%val(:,:,:,:)
00402     call write_FieldData(writer, self%errorh_ncVarID, self%error_h)
00403     call increase_recorde_counter(writer)
00404 
00405   end if
00406 
00407 end subroutine timelevel_Updated
00408 
00409 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00410 !
00411 ! 非公開手続き
00412 ! 
00413 !!!!!!!!!!!!!!!!!!!!!!!!!
00414 
00415 !
00421 subroutine eval_numcal_h_solution(t, true_h, h_n)
00422 
00423   ! 宣言文 ; Declaration statement
00424   !
00425   real(DP), intent(in) :: t
00426   type(Field_IcGrid2D), intent(inout) :: true_h
00427   type(Field_IcGrid2D), intent(in) :: h_n
00428   
00429   ! 作業変数
00430   ! Work variables
00431   !
00432   real(DP) norm1Val, norm2Val, normInfityVal
00433 
00434   ! 実行文 ; Executable statement
00435   !
00436   
00437   !
00438   norm1Val = numerical_error_norm1(true_h, h_n)
00439   norm2Val = numerical_error_norm2(true_h, h_n)
00440   normInfityVal = numerical_error_norminfinity(true_h, h_n, maxmin_info_flag=.false.)
00441 
00442   call HistoryPut('norm1', norm1Val)
00443   call HistoryPut('norm2', norm2Val)
00444   call HistoryPut('norm_infinity', normInfityVal)
00445 
00446 end subroutine eval_numcal_h_solution
00447 
00448 end module class_TestCase2
 全て クラス ネームスペース ファイル 関数 変数