IGModel-SW 1.0

src/Williamson1992_TestCase/class_TestCase6.f90

説明を見る。
00001 
00009 module class_TestCase6
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     & geo_to_orth_vec
00040 
00041   ! 正二十面体格子の管理クラス
00042   !
00043   !
00044   use IcGrid2D_FVM_Manager, only: &
00045     & IcGrid2D_FVM, &
00046     & get_EffSize_Min, get_EffSize_Max, get_glevel, get_IcRadius, check_pole, &
00047     & RC_REGIONS_NUM, NOT_POLE_FLAG, NORTH_POLE_FLAG, SOUTH_POLE_FLAG
00048 
00049   ! 正二十面体格子上に分布する物理場の管理クラス
00050   !
00051   !
00052   use Field_IcGrid2D_manager, only: &
00053     & Field_IcGrid2D, &
00054     & Field_IcGrid2D_Init, Field_IcGrid2D_Final, &
00055     & get_icgrid, paste_field2D_margin
00056   
00057 
00058   !
00059   !
00060   !
00061   use Field_Pattern_Builder, only: &
00062     & create_solid_rotation_field
00063 
00064 
00065   !
00066   !
00067   !
00068   use IcGrid_ncWriter_mod, only: &
00069     & IcGrid_ncWriter, IcGrid_ncWriter_Init, &
00070     & open_ncFile, ncdef_simulation_parameter, ncdef_GridData, ncdef_FieldData, end_ncdef, &
00071     & write_GridData, write_FieldData, close_ncFile, &
00072     & increase_recorde_counter
00073 
00074 
00075   ! * IGModel-SW ****
00076   ! *
00077 
00078 
00079   ! パラメータ管理
00080   !
00081   !
00082   use param_manager, only: &
00083     & alpha, output_tick, Grav, Omega, &
00084     & integration_time, DelTime
00085 
00086   !
00087   !
00088   !
00089   use field_manager, only: &
00090     & diff_eval, xy_Coli
00091 
00092   !
00093   !
00094   !
00095   use conservation_analysis, only: &
00096     & init_conservation_analysis, &
00097     & calc_total_energy, calc_potential_enstrophy
00098 
00099   ! 宣言文 ; Declaration statement
00100   !
00101   implicit none
00102   private
00103   
00104   ! 公開変数
00105   ! Public variables
00106   !
00107 
00109   real(DP) :: h_0 
00110 
00112   real(DP), parameter :: ome = 7.848d-06
00113 
00115   real(DP), parameter :: K =  7.848d-06
00116 
00118   real(DP), parameter :: R = 4.0d0
00119  
00120   ! クラス定義
00121   ! Class definition
00122   !
00123 
00124   !
00128   type, public :: TestCase6
00129 
00130     ! 宣言文 ; Declaration statement
00131     !
00132     integer :: dummy
00133   end type TestCase6
00134 
00135   !
00136   !
00137   !
00138   interface initialize_TestCase
00139     module procedure init_TestCase6
00140   end interface initialize_TestCase
00141 
00142   interface finalize_TestCase
00143     module procedure finalize_TestCase6
00144   end interface finalize_TestCase
00145 
00146   public :: initialize_TestCase, finalize_TestCase
00147   public :: set_initial_h, set_initial_hs, set_initial_v, timelevel_Updated
00148 
00149 
00150   ! 非公開変数
00151   ! Private variables
00152   !
00153   
00155   real(DP) :: init_TE
00156 
00158   real(DP) :: init_PE
00159 
00161   real(DP) :: TE
00162 
00164   real(DP) :: PE
00165 
00167   type(Field_IcGrid2D), pointer :: ini_v
00168 
00170   type(Field_IcGrid2D), pointer :: ini_h
00171 
00173   type(Field_IcGrid2D), pointer :: ini_hs
00174 
00175 contains
00176 
00177 !
00183 subroutine init_TestCase6(self, icgrid_ref )
00184   ! 宣言文 ; Declaration statement
00185   !
00186   type(TestCase6), intent(inout) :: self
00187   type(IcGrid2D_FVM), intent(in) :: icgrid_ref
00188 
00189   ! 作業変数
00190   ! Work variables
00191   !
00192 
00193   ! 実行文 ; Executable statement
00194   !
00195   h_0 = 8.0e03
00196 
00197   !
00198   allocate(ini_v)
00199   allocate(ini_h)
00200   allocate(ini_hs)
00201 
00202   call Field_IcGrid2D_Init(ini_v, icgrid_ref, 'init_v', 3)
00203   call Field_IcGrid2D_Init(ini_h, icgrid_ref, 'init_h', 1)
00204   call Field_IcGrid2D_Init(ini_hs, icgrid_ref, 'init_hs', 1)
00205 
00206   call calc_init_v()
00207   call paste_field2D_margin(ini_v)
00208   call calc_init_h()
00209   call calc_init_hs()
00210 
00211   call init_conservation_analysis(ini_h, ini_hs)
00212   init_TE = calc_total_energy(ini_h, ini_hs, ini_v)
00213   init_PE = calc_potential_enstrophy(ini_h, ini_hs, ini_v, xy_Coli, diff_eval)
00214 
00215   write(*,*) 'initTE', init_TE, init_PE
00216   open(10, file='TE_PE.dat',status='unknown')
00217 
00218   !
00219   !
00220 
00221   !
00222   call HistoryCreate( &
00223     & file='dissipation_rate.nc', title='disspation rates', &
00224     & source=' ', institution=' ', &
00225     & dims=(/ 't' /), dimsizes=(/ 0 /), &
00226     & longnames=(/ 'time' /), units=(/ 's' /), &
00227     & origin=real(0d0), interval=real(output_tick) )
00228 
00229   call HistoryAddVariable( &
00230     & varname='TE', dims=(/ 't' /), &
00231     & longname='total energy', units='m m2 s-2', xtype='double' )
00232 
00233   call HistoryAddVariable( &
00234     & varname='PE', dims=(/ 't' /), &
00235     & longname='potential enstrophy', units='s-2 m-1', xtype='double' )
00236 
00237   call HistoryAddVariable( &
00238     & varname='TE_DISSIP_RATE', dims=(/ 't' /), &
00239     & longname='dissipation rate of total energy', units='1', xtype='double' )
00240 
00241   call HistoryAddVariable( &
00242     & varname='PE_DISSIP_RATE', dims=(/ 't' /), &
00243     & longname='dissipation rate of potential enstrophy', units='1', xtype='double' )
00244 
00245   call HistoryAddAttr('TE', '+glevel', get_glevel(icgrid_ref))
00246   call HistoryAddAttr('TE', '+alpha', alpha)
00247 
00248 end subroutine init_TestCase6
00249 
00250 !
00255 subroutine finalize_TestCase6(self)
00256   ! 宣言文 ; Declaration statement
00257   !
00258   type(TestCase6), intent(inout) :: self
00259 
00260   ! 実行文 ; Executable statement
00261   !
00262   call HistoryClose()
00263 
00264 end subroutine finalize_TestCase6
00265 
00266 !
00272 subroutine set_initial_v(self, init_v)
00273   ! 宣言文 ; Declaration statement
00274   !
00275   type(TestCase6), intent(inout) :: self
00276   type(Field_IcGrid2D), intent(inout) :: init_v
00277 
00278   ! 実行文 ; Executable statement
00279   !
00280 
00281   write(*,*) 'init_TestCase6_v'
00282 
00283   init_v%val = ini_v%val
00284 
00285   call Field_IcGrid2D_Final(ini_v)
00286   deallocate(ini_v)
00287 
00288 end subroutine set_initial_v
00289 
00290 !
00296 subroutine set_initial_h(self, init_h)
00297   ! 宣言文 ; Declaration statement
00298   !
00299   type(TestCase6), intent(inout) :: self
00300   type(Field_IcGrid2D), intent(inout) :: init_h
00301 
00302   ! 実行文 ; Executable statement
00303   !
00304   init_h%val = ini_h%val
00305 
00306   call Field_IcGrid2D_Final(ini_h)
00307   deallocate(ini_h)
00308 
00309 end subroutine set_initial_h
00310 
00311 !
00320 subroutine set_initial_hs(self, init_hs)
00321 
00322   ! 宣言文 ; Declaration statement
00323   !
00324   type(TestCase6), intent(inout) :: self
00325   type(Field_IcGrid2D), intent(inout) :: init_hs
00326 
00327   ! 実行文 ; Executable statement
00328   !
00329   init_hs%val = ini_hs%val
00330 
00331 end subroutine set_initial_hs
00332 
00333 !
00345 subroutine timelevel_Updated(self, tstep, dt, v_n, h_n)
00346 
00347   ! 宣言文 ; Declaration statement
00348   !
00349   type(TestCase6), intent(inout) :: self
00350   integer, intent(in) :: tstep
00351   real(DP), intent(in) :: dt
00352   type(Field_IcGrid2D), intent(inout) :: v_n
00353   type(Field_IcGrid2D), intent(inout) :: h_n
00354 
00355   ! 作業変数
00356   ! Work variables
00357   !
00358   real(DP) :: TE
00359   real(DP) :: PE
00360 
00361   ! 実行文 ; Executable statement
00362   !
00363   
00364 
00365   !write(*,*) TE, PE
00366   if ( mod( tstep, int(output_tick / dt) ) == 0 ) then
00367     TE = calc_total_energy(h_n, ini_hs, v_n)
00368     PE = calc_potential_enstrophy(h_n, ini_hs, v_n, xy_Coli, diff_eval)
00369 
00370     call HistoryPut('TE', TE)
00371     call HistoryPut('PE', PE)
00372 
00373     call HistoryPut('TE_DISSIP_RATE', abs( TE - init_TE ) / init_TE )
00374     call HistoryPut('PE_DISSIP_RATE', abs( PE - init_PE ) / init_PE )
00375 
00376   end if
00377 
00378 end subroutine timelevel_Updated
00379 
00380 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00381 !
00382 ! 非公開手続き
00383 ! 
00384 !!!!!!!!!!!!!!!!!!!!!!!!!
00385 
00386 !
00387 !
00388 !
00389 subroutine calc_init_v()
00390 
00391   ! 実行文 ; Executable statement
00392   !
00393 
00394   !
00395   ! Work Variables
00396   !
00397   integer :: rcID, i,j
00398   integer :: EMin, EMax
00399   type(IcGrid2D_FVM), pointer :: icgrid
00400   real(DP) :: s_pos(3)
00401   real(DP) :: a
00402   real(DP) :: s_u, s_v
00403 
00404   ! 実行文 ; Executable statement
00405   !
00406 
00407   write(*,*) 'init_TestCase6_v'
00408 
00409   icgrid => get_IcGrid(ini_h)
00410   a = get_IcRadius(icgrid)
00411   EMin = get_EffSize_Min(icgrid)
00412   EMax = get_EffSize_Max(icgrid)
00413 
00414   do rcID= 1, RC_REGIONS_NUM
00415     do j=EMin, EMax
00416       do i=EMin, EMax
00417         s_pos(:) = orth_to_geo_pos(icgrid%rcs_AGrid(rcID,i,j,:))
00418         if( check_pole(icgrid, rcID,i,j) == NOT_POLE_FLAG ) then
00419           !
00420            s_u =   a * ome * cos(s_pos(2)) &
00421              &   + a * K * cos(s_pos(2))**(R-1) * (R*sin(s_pos(2))**2 - cos(s_pos(2))**2) * cos(R*s_pos(1))
00422 
00423            s_v = - a * K * R * cos(s_pos(2))**(R-1) * sin(s_pos(2)) * sin(R * s_pos(1))
00424 
00425            ini_v%val(rcID,i,j,:) = geo_to_orth_vec((/ s_u, s_v, 0.0d0 /), s_pos(:))
00426 
00427         else
00428           ini_v%val(rcID,i,j,:) = 0.0d0
00429         end if
00430 
00431       end do
00432     end do
00433   end do
00434 
00435 end subroutine calc_init_v
00436 
00437 subroutine calc_init_h()
00438   ! 宣言文 ; Declaration statement
00439   !
00440 
00441   !
00442   ! Work Variables
00443   !
00444   integer :: rcID, i,j
00445   integer :: EMin, EMax
00446   type(IcGrid2D_FVM), pointer :: icgrid
00447   real(DP) :: s_pos(3)
00448   real(DP) :: ic_radius
00449   real(DP) :: tmp
00450 
00451   ! 実行文 ; Executable statement
00452   !
00453 
00454   icgrid => get_IcGrid(ini_h)
00455   ic_radius = get_IcRadius(icgrid)
00456   EMin = get_EffSize_Min(icgrid)
00457   EMax = get_EffSize_Max(icgrid)
00458 
00459   do rcID= 1, RC_REGIONS_NUM
00460     do j=EMin, EMax
00461       do i=EMin, EMax
00462 
00463         s_pos(:) = orth_to_geo_pos(icgrid%rcs_AGrid(rcID,i,j,:))
00464         if( check_pole(icgrid, rcID,i,j) == NOT_POLE_FLAG ) then
00465           !
00466           tmp =   Grav * h_0 + ic_radius**2 * A(s_pos(2)) &
00467             &   + ic_radius**2 * B(s_pos(2)) * cos(R * s_pos(1)) &
00468             &   + ic_radius**2 * C(s_pos(2)) * cos(2.0d0 * R * s_pos(1))
00469         else
00470           !
00471           tmp = Grav * h_0
00472         end if
00473 
00474         ini_h%val(rcID,i,j, :) = tmp / Grav
00475 
00476      end do
00477     end do
00478   end do
00479 
00480 end subroutine calc_init_h
00481 
00482 !
00483 !
00484 !
00485 subroutine calc_init_hs()
00486   ! 宣言文 ; Declaration statement
00487   !
00488 
00489   ! 作業変数
00490   ! Work Variables
00491   !
00492   integer :: rcID, i,j
00493   integer :: EMin, EMax
00494   type(IcGrid2D_FVM), pointer :: icgrid
00495   real(DP) :: s_pos(3)
00496   real(DP) :: ic_radius
00497   real(DP) :: r
00498 
00499   ! 実行文 ; Executable statement
00500   !
00501 
00502   write(*,*) 'init_TestCase6_hs'
00503 
00504   ini_hs%val = 0.0d0
00505 
00506 end subroutine calc_init_hs
00507 
00508 function A(theta) result(val)
00509 
00510   real(DP), intent(in) :: theta
00511   real(DP) val
00512 
00513   val =    0.5d0 * ome * (2 * Omega + ome) * cos(theta)**2 &
00514     &    + 0.25d0 * K**2 * cos(theta) **(2*R) * ( &
00515     &       (R+1.0d0) * cos(theta)**2 + ( 2.0d0 * R**2 - R - 2.0d0 ) - 2.0d0 * R**2 * cos(theta)**(-2) &
00516     &        )
00517 
00518 end function A
00519 
00520 function B(theta) result(val)
00521 
00522   real(DP), intent(in) :: theta
00523   real(DP) val
00524 
00525   val =  2.0d0 * (Omega + ome) * K / ( ( R + 1.0d0 ) * ( R + 2.0d0 ) ) * cos(theta)**R * &
00526     &       ( ( R**2 + 2.0d0 * R + 2.0d0 ) - ( R + 1.0d0 )**2 * cos(theta)**2 )
00527 
00528 end function B
00529 
00530 function C(theta) result(val)
00531 
00532   real(DP), intent(in) :: theta
00533   real(DP) val
00534 
00535   val =  0.25d0 * K**2 * cos(theta)**(2*R) * ( ( R + 1.0d0 ) * cos(theta)**2 - ( R + 2.0d0 ) )
00536 
00537 end function C
00538 
00539 end module class_TestCase6
 全て クラス ネームスペース ファイル 関数 変数