IGModel-SW 1.0

src/Williamson1992_TestCase/class_TestCase5.f90

説明を見る。
00001 
00009 module class_TestCase5
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: &
00032     & PI, ToRadians
00033 
00034   ! 座標変換
00035   ! Coordinate conversion
00036   !
00037   use igmcore_coordinate_conversion, only: &
00038     & orth_to_geo_pos, geo_to_orth_pos
00039 
00040   ! 正二十面体格子の管理クラス
00041   !
00042   !
00043   use IcGrid2D_FVM_Manager, only: &
00044     & IcGrid2D_FVM, &
00045     & get_EffSize_Min, get_EffSize_Max, get_glevel, get_IcRadius, check_pole, &
00046     & RC_REGIONS_NUM, NOT_POLE_FLAG, NORTH_POLE_FLAG, SOUTH_POLE_FLAG
00047 
00048   ! 正二十面体格子上に分布する物理場の管理クラス
00049   !
00050   !
00051   use Field_IcGrid2D_manager, only: &
00052     & Field_IcGrid2D, &
00053     & Field_IcGrid2D_Init, Field_IcGrid2D_Final, &
00054     & get_icgrid, paste_field2D_margin
00055 
00056 
00057   !
00058   !
00059   !
00060   use Field_Pattern_Builder, only: &
00061     & create_solid_rotation_field
00062 
00063   ! 正二十面体格子上の物理場における誤差ノルムチェック
00064   !
00065   !
00066   use Field_IcGrid2D_error_norm_check, only: &
00067     & numerical_error_norm1, numerical_error_norm2, &
00068     & numerical_error_norminfinity 
00069 
00070   !
00071   !
00072   !
00073   use IcGrid_ncWriter_mod, only: &
00074     & IcGrid_ncWriter, IcGrid_ncWriter_Init, &
00075     & open_ncFile, ncdef_simulation_parameter, ncdef_GridData, ncdef_FieldData, end_ncdef, &
00076     & write_GridData, write_FieldData, close_ncFile, &
00077     & increase_recorde_counter
00078 
00079   ! * IGModel-SW ****
00080   ! *
00081 
00082   ! パラメータ管理
00083   !
00084   !
00085   use param_manager, only: &
00086     & alpha, output_tick, Grav, Omega, &
00087     & integration_time, DelTime
00088 
00089   !
00090   !
00091   !
00092   use field_manager, only: &
00093     & diff_eval, xy_Coli
00094 
00095   !
00096   !
00097   !
00098   use conservation_analysis, only: &
00099     & init_conservation_analysis, &
00100     & calc_total_energy, calc_potential_enstrophy
00101 
00102   ! 宣言文 ; Declaration statement
00103   !
00104   implicit none
00105   private
00106   
00107   ! 公開変数
00108   ! Public variables
00109   !
00110 
00112   real(DP) :: h_0 
00113 
00116   real(DP), parameter :: hs_0 = 2000.0d0
00117 
00120   real(DP), parameter :: theta_c = PI / 6.0d0
00121 
00124   real(DP), parameter :: lambda_c = - PI / 2.0d0
00125   
00128   real(DP), parameter :: cb_R = PI / 9.0d0
00129  
00132   real(DP), parameter :: u_0 = 20.0d0
00133 
00136   real(DP) :: angular_speed
00137 
00138   ! クラス定義
00139   ! Class definition
00140   !
00141 
00142   !
00146   type, public :: TestCase5
00147 
00148     ! 宣言文 ; Declaration statement
00149     !
00150     integer :: dummy
00151   end type TestCase5
00152 
00153   !
00154   !
00155   !
00156   interface initialize_TestCase
00157     module procedure init_TestCase5
00158   end interface initialize_TestCase
00159 
00160   interface finalize_TestCase
00161     module procedure finalize_TestCase5
00162   end interface finalize_TestCase
00163 
00164   public :: initialize_TestCase, finalize_TestCase
00165   public :: set_initial_h, set_initial_hs, set_initial_v, timelevel_Updated
00166 
00167 
00168   ! 非公開変数
00169   ! Private variables
00170   !
00171   
00173   real(DP) :: init_TE
00174 
00176   real(DP) :: init_PE
00177 
00179   real(DP) :: TE
00180 
00182   real(DP) :: PE
00183 
00185   type(Field_IcGrid2D), pointer :: ini_v
00186 
00188   type(Field_IcGrid2D), pointer :: ini_h
00189 
00191   type(Field_IcGrid2D), pointer :: ini_hs
00192 
00193 contains
00194 
00195 !
00201 subroutine init_TestCase5(self, icgrid_ref )
00202   ! 宣言文 ; Declaration statement
00203   !
00204   type(TestCase5), intent(inout) :: self
00205   type(IcGrid2D_FVM), intent(in) :: icgrid_ref
00206 
00207   ! 作業変数
00208   ! Work variables
00209   !
00210   real(DP) :: ic_radius
00211 
00212   ! 実行文 ; Executable statement
00213   !
00214 
00215   ic_radius = get_IcRadius(icgrid_ref)
00216   h_0 = 5960d0!2.94e04 / g
00217   angular_speed = u_0 / ic_radius ! 2.0d0 * PI / ( 12.0d0 * 24.0d0 * 3600d0 )
00218 
00219   !
00220   allocate(ini_v)
00221   allocate(ini_h)
00222   allocate(ini_hs)
00223 
00224   call Field_IcGrid2D_Init(ini_v, icgrid_ref, 'init_v', 3)
00225   call Field_IcGrid2D_Init(ini_h, icgrid_ref, 'init_h', 1)
00226   call Field_IcGrid2D_Init(ini_hs, icgrid_ref, 'init_hs', 1)
00227 
00228   call calc_init_v()
00229   call paste_field2D_margin(ini_v)
00230   call calc_init_h()
00231   call calc_init_hs()
00232 
00233   call init_conservation_analysis(ini_h, ini_hs)
00234   init_TE = calc_total_energy(ini_h, ini_hs, ini_v)
00235   init_PE = calc_potential_enstrophy(ini_h, ini_hs, ini_v, xy_Coli, diff_eval)
00236 
00237   write(*,*) 'initTE', init_TE, init_PE
00238   open(10, file='TE_PE.dat',status='unknown')
00239 
00240   !
00241   !
00242 
00243   !
00244   call HistoryCreate( &
00245     & file='dissipation_rate.nc', title='disspation rates', &
00246     & source=' ', institution=' ', &
00247     & dims=(/ 't' /), dimsizes=(/ 0 /), &
00248     & longnames=(/ 'time' /), units=(/ 's' /), &
00249     & origin=real(0d0), interval=real(output_tick) )
00250 
00251   call HistoryAddVariable( &
00252     & varname='TE', dims=(/ 't' /), &
00253     & longname='total energy', units='m m2 s-2', xtype='double' )
00254 
00255   call HistoryAddVariable( &
00256     & varname='PE', dims=(/ 't' /), &
00257     & longname='potential enstrophy', units='s-2 m-1', xtype='double' )
00258 
00259   call HistoryAddVariable( &
00260     & varname='TE_DISSIP_RATE', dims=(/ 't' /), &
00261     & longname='dissipation rate of total energy', units='1', xtype='double' )
00262 
00263   call HistoryAddVariable( &
00264     & varname='PE_DISSIP_RATE', dims=(/ 't' /), &
00265     & longname='dissipation rate of potential enstrophy', units='1', xtype='double' )
00266 
00267   call HistoryAddAttr('TE', '+glevel', get_glevel(icgrid_ref))
00268   call HistoryAddAttr('TE', '+alpha', alpha)
00269 
00270 end subroutine init_TestCase5
00271 
00272 !
00277 subroutine finalize_TestCase5(self)
00278   ! 宣言文 ; Declaration statement
00279   !
00280   type(TestCase5), intent(inout) :: self
00281 
00282   ! 実行文 ; Executable statement
00283   !
00284   call HistoryClose()
00285 
00286 end subroutine finalize_TestCase5
00287 
00288 !
00294 subroutine set_initial_v(self, init_v)
00295   ! 宣言文 ; Declaration statement
00296   !
00297   type(TestCase5), intent(inout) :: self
00298   type(Field_IcGrid2D), intent(inout) :: init_v
00299 
00300   ! 実行文 ; Executable statement
00301   !
00302 
00303   write(*,*) 'init_Testcase5_v'
00304 
00305   init_v%val = ini_v%val
00306 
00307   call Field_IcGrid2D_Final(ini_v)
00308   deallocate(ini_v)
00309 
00310 end subroutine set_initial_v
00311 
00312 !
00318 subroutine set_initial_h(self, init_h)
00319   ! 宣言文 ; Declaration statement
00320   !
00321   type(TestCase5), intent(inout) :: self
00322   type(Field_IcGrid2D), intent(inout) :: init_h
00323 
00324   ! 実行文 ; Executable statement
00325   !
00326   init_h%val = ini_h%val
00327 
00328   call Field_IcGrid2D_Final(ini_h)
00329   deallocate(ini_h)
00330 
00331 end subroutine set_initial_h
00332 
00333 !
00342 subroutine set_initial_hs(self, init_hs)
00343 
00344   ! 宣言文 ; Declaration statement
00345   !
00346   type(TestCase5), intent(inout) :: self
00347   type(Field_IcGrid2D), intent(inout) :: init_hs
00348 
00349   ! 実行文 ; Executable statement
00350   !
00351   init_hs%val = ini_hs%val
00352 
00353 end subroutine set_initial_hs
00354 
00355 !
00367 subroutine timelevel_Updated(self, tstep, dt, v_n, h_n)
00368 
00369   ! 宣言文 ; Declaration statement
00370   !
00371   type(TestCase5), intent(inout) :: self
00372   integer, intent(in) :: tstep
00373   real(DP), intent(in) :: dt
00374   type(Field_IcGrid2D), intent(inout) :: v_n
00375   type(Field_IcGrid2D), intent(inout) :: h_n
00376 
00377   ! 作業変数
00378   ! Work variables
00379   !
00380   real(DP) :: TE
00381   real(DP) :: PE
00382 
00383   ! 実行文 ; Executable statement
00384   !
00385   
00386 
00387   !write(*,*) TE, PE
00388   if ( mod( tstep, int(output_tick / dt) ) == 0 ) then
00389     TE = calc_total_energy(h_n, ini_hs, v_n)
00390     PE = calc_potential_enstrophy(h_n, ini_hs, v_n, xy_Coli, diff_eval)
00391 
00392     call HistoryPut('TE', TE)
00393     call HistoryPut('PE', PE)
00394 
00395     call HistoryPut('TE_DISSIP_RATE', abs( TE - init_TE ) / init_TE )
00396     call HistoryPut('PE_DISSIP_RATE', abs( PE - init_PE ) / init_PE )
00397 
00398   end if
00399 
00400 end subroutine timelevel_Updated
00401 
00402 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00403 !
00404 ! 非公開手続き
00405 ! 
00406 !!!!!!!!!!!!!!!!!!!!!!!!!
00407 
00408 subroutine calc_init_h()
00409   ! 宣言文 ; Declaration statement
00410   !
00411 
00412   !
00413   ! Work Variables
00414   !
00415   integer :: rcID, i,j
00416   integer :: EMin, EMax
00417   type(IcGrid2D_FVM), pointer :: icgrid
00418   real(DP) :: s_pos(3)
00419   real(DP) :: ic_radius
00420   real(DP) :: tmp
00421 
00422   ! 実行文 ; Executable statement
00423   !
00424 
00425   icgrid => get_IcGrid(ini_h)
00426   ic_radius = get_IcRadius(icgrid)
00427   EMin = get_EffSize_Min(icgrid)
00428   EMax = get_EffSize_Max(icgrid)
00429 
00430   do rcID= 1, RC_REGIONS_NUM
00431     do j=EMin, EMax
00432       do i=EMin, EMax
00433 
00434         s_pos(:) = orth_to_geo_pos(icgrid%rcs_AGrid(rcID,i,j,:))
00435         if( check_pole(icgrid, rcID,i,j) == NOT_POLE_FLAG ) then
00436           !
00437           tmp = cos(s_pos(1)) * cos(s_pos(2)) * sin(alpha) + sin(s_pos(2)) * cos(alpha)
00438         else
00439           !
00440           tmp = sin(s_pos(2)) * cos(alpha)
00441         end if
00442 
00443         ini_h%val(rcID,i,j, :) = &
00444           & h_0 - ( ic_radius * Omega * u_0 + 0.5d0 * u_0**2 ) * tmp**2 / Grav
00445      end do
00446     end do
00447   end do
00448 
00449 end subroutine calc_init_h
00450 
00451 !
00452 !
00453 !
00454 subroutine calc_init_hs()
00455   ! 宣言文 ; Declaration statement
00456   !
00457 
00458   ! 作業変数
00459   ! Work Variables
00460   !
00461   integer :: rcID, i,j
00462   integer :: EMin, EMax
00463   type(IcGrid2D_FVM), pointer :: icgrid
00464   real(DP) :: s_pos(3)
00465   real(DP) :: ic_radius
00466   real(DP) :: r
00467 
00468   ! 実行文 ; Executable statement
00469   !
00470 
00471   write(*,*) 'init_Testcase5_hs'
00472 
00473   icgrid => get_IcGrid(ini_hs)
00474   ic_radius = get_IcRadius(icgrid)
00475   EMin = get_EffSize_Min(icgrid)
00476   EMax = get_EffSize_Max(icgrid)
00477 
00478   ! 下部境界の高度場として, 孤立した山岳地形を設定する.
00479   !
00480   do rcID= 1, RC_REGIONS_NUM
00481     do j=EMin, EMax
00482       do i=EMin, EMax
00483 
00484         s_pos(:) = orth_to_geo_pos(icgrid%rcs_AGrid(rcID,i,j,:))
00485         if( check_pole(icgrid, rcID,i,j) /= NOT_POLE_FLAG ) then
00486           s_pos(1) = 0.0d0
00487         end if
00488 
00489         r = sqrt(min(cb_R**2, (s_pos(1) - lambda_c)**2 + (s_pos(2) - theta_c)**2 ))
00490 !        write(*,*) cb_R**2, (s_pos(1) - lambda_c)**2 + (s_pos(2) - theta_c)**2
00491 
00492         ini_hs%val(rcID,i,j, :) = hs_0 * ( 1.0d0 - r / cb_R )
00493         if(rcID==1 .and. i == 1 .and. j==1 ) write(*,*) 'hs:', ini_hs%val(rcID,i,j,:), s_pos
00494      end do
00495     end do
00496   end do
00497 end subroutine calc_init_hs
00498 
00499 !
00500 !
00501 !
00502 subroutine calc_init_v()
00503 
00504   ! 実行文 ; Executable statement
00505   !
00506 
00507   write(*,*) 'init_Testcase5_v'
00508 
00509   call create_solid_rotation_field(ini_v, angular_speed, alpha)
00510 
00511 end subroutine calc_init_v
00512 
00513 end module class_TestCase5
 全て クラス ネームスペース ファイル 関数 変数