| Path: | main/init.f90 |
| Last Update: | Mon Aug 21 16:52:12 JST 2006 |
| Main Program : |
初期値を生成するためのプログラムです.
This procedure input/output NAMELIST#init_nml .
program init
!
!初期値を生成するためのプログラムです.
!
!----- 構造体参照モジュール -----
use type_mod, only: STRING, INTKIND, REKIND, DBKIND
use nmlfile_mod, only: nmlfile_init, nmlfile_open, nmlfile_close
use constants_mod, only: constants_init, pi, R0
!----- 格子点取得モジュール -----
use grid_3d_mod, only: grid_3d_init, im, jm, km, grid_3d_end
use grid_wavenumber_mod, only: grid_wavenumber_init, nm, grid_wavenumber_end
!----- 座標データ生成モジュール -----
use axis_type_mod, only: AXISINFO
use axis_x_mod, only: axis_x_init, axis_x_weight, axis_x_spectral, axis_x_end
use axis_y_mod, only: axis_y_init, axis_y_weight, axis_y_spectral, axis_y_end
use axis_z_mod, only: axis_z_init, axis_z_sigmahalf_manual, axis_z_end
!----- データI/Oモジュール -----
use io_gt4_out_mod, only : io_gt4_out_init , io_gt4_out_SetDims, io_gt4_out_SetVars, io_gt4_out_Put , io_gt4_out_end
!----- 時刻管理モジュール -----
use time_mod, only: time_init, tvar, ttype, tname, tunit, time_end
!----- SPMODEL モジュール -----
use spml_mod, only: spml_init, xya_wa, wa_Div_xya_xya ,wa_LaplaInv_wa, wa_xya, xya_GradLat_wa, xya_GradLon_wa
!----- コマンドライン引数解析ツール -----
use dc_args, only: ARGS, Open, Debug, Help, Strict, Close
!----- デバッグ・汎用ツール -----
use dc_trace, only: SetDebug, DbgMessage, BeginSub, EndSub, DataDump
use dc_string, only: toChar, StriEq, LChar, StrHead
use dc_message,only: MessageNotify
use gt4_history , only: HistoryAxisInquire
implicit none
character(STRING) :: condition = '' ! 初期値の種類
real(DBKIND) :: VelLonAve = 0.0d0 ! 速度経度成分平均値
real(DBKIND) :: VelLatAve = 0.0d0 ! 速度緯度成分平均値
real(DBKIND) :: VorAve = 0.0d0 ! 渦度平均値
real(DBKIND) :: DivAve = 0.0d0 ! 発散平均値
real(DBKIND) :: TempAve = 273.0d0 ! 温度平均値
real(DBKIND) :: QVapAve = 0.0d0 ! 比湿平均値
real(DBKIND) :: PsAve = 1.0d5 ! 地表面圧力平均値
logical :: VorDiv_Priority = .false. ! 渦度発散から風速を生成
! !for 'rigid body rotation'
real(DBKIND) :: VelLonMax_rbr = 1.0d2 ! Maximum of 'VelLon'
! !for 'convex of surface pressure' or 'convex of temperature'
real(DBKIND) :: LonLat_Radius_Deg = 20.0 ! 半径 (度数)
real(DBKIND) :: LonLat_Radius_Rad = 0.349 ! 半径 (ラジアン)
real(DBKIND) :: Lat_Center_Deg = 45.0 ! 緯度の中心位置 (度数)
real(DBKIND) :: Lon_Center_Deg = 100.0 ! 経度の中心位置 (度数)
real(DBKIND) :: Lat_Center_Rad = 0.785 ! 緯度の中心位置 (ラジアン)
real(DBKIND) :: Lon_Center_Rad = 1.745 ! 経度の中心位置 (ラジアン)
logical :: Rad_Priority = .false. ! ラジアン表記を優先
! !for 'convex of surface pressure'
real(DBKIND) :: PsMax = -200.0d2 ! Maximum of 'Ps'
! !for 'convex of temperature'
real(DBKIND) :: TempMax = 10.0d0 ! Maximum of 'Temp'
! !for 'agcm5.3 default value'
real(DBKIND) :: TPRTRB = 0.1 ! 温度擾乱最大値
namelist /init_nml/ condition, VelLonAve, VelLatAve, VorAve, DivAve, TempAve, QVapAve, PsAve, VorDiv_Priority, VelLonMax_rbr, LonLat_Radius_Deg, LonLat_Radius_Rad, Lat_Center_Deg, Lon_Center_Deg, Lat_Center_Rad, Lon_Center_Rad, Rad_Priority, PsMax, TempMax
!
! 初期値の種類の設定を行なう。
! 現在 condition に与えて有効なのは以下の値である。
!
! * <b>rigid body rotation</b>
! * 剛体回転流を与える。風速の最大値は VelLonMax_rbr に与える。
! * これを選択した場合、強制的に VorDiv_Priority は .false.
! に設定される。
!
! * <b> convex of surface pressure </b>
! * 地表面気圧の「山」を与える。位置とサイズは、度数で与える場合は
! Lon_Center_Dig, Lat_Center_Dig, LonLatRadius_Dig を用い、
! ラジアンで与える場合は Lon_Center_Rad, Lat_Center_Rad,
! LonLatRadius_Rad を用いる (ただし、Rad_Priority を .true.
! にする必要がある)。 最大値は PsMax で与える。
!
! * <b> convex of temperature </b>
! * 温度の「山」を与える。位置とサイズは、度数で与える場合は
! Lon_Center_Dig, Lat_Center_Dig, LonLatRadius_Dig を用い、
! ラジアンで与える場合は Lon_Center_Rad, Lat_Center_Rad,
! LonLatRadius_Rad を用いる (ただし、Rad_Priority を .true.
! にする必要がある)。 最大値は TempMax で与える。
!
! * <b> agcm5.3 default value </b>
! * 等温静止大気. 温度場に擾乱を与える. 既定値は以下.
! QVapAve = 1.0d-10
! TempAve = 250.0d0
! PsAve = 1.0d5
! TPRTRB = 0.1
!
! * その他
! * 等温無風の初期値を与える。
!
! VelLonAve, VelLonAve, VelLatAve, VorAve, DivAve, TempAve, QVapAve, PsAve
! には、それぞれ平均値を与える。
!
! デフォルトでは風速から渦度発散を生成するが、
! VorDiv_Priority を .true. にする事で、渦度発散から風速を生成する。
!
!For Generate AXES
type(AXISINFO):: x_Lon , x_Lon_Weight , y_Lat , y_Lat_Weight , z_Sigma , r_Sigma ! σレベル(半整数)座標
!For Generate Initial Data
real(DBKIND), pointer :: xyz_VelLon(:,:,:) , xyz_VelLat(:,:,:) , xyz_Vor(:,:,:) , xyz_Div(:,:,:) , xyz_Temp(:,:,:) , xyz_QVap(:,:,:) , xyz_Ps(:,:,:) ! 格子点データ(地表面気圧)
!Axes data (for Use as 3-dimensional data)
real(DBKIND), allocatable :: xyz_Lon(:,:,:) , xyz_Lat(:,:,:) , xyz_Sigma(:,:,:) , xyr_Sigma(:,:,:) ! σレベル(半整数)座標
real(DBKIND) :: RadDegFact ! ラジアンを度数に変換する係数
!for Generate Velocity from Vorticity and Divergence
real(DBKIND), allocatable :: wz_Psi(:,:) , wz_Chi(:,:) ! スペクトル(ポテンシャル)
!for nmlfile_mod
integer(INTKIND) :: nmlstat, nmlunit
logical :: nmlreadable
!for dc_args
type(ARGS) :: arg
!----- 作業用内部変数 -----
integer(INTKIND) :: i,j,k
character(STRING) :: axis_units
character(*), parameter:: version = '$Name: dcpam3-20070608 $' // '$Id: init.f90,v 1.5 2006/08/21 07:52:12 morikawa Exp $'
character(STRING), parameter:: subname = "init"
!-------------------------------------------------------------------
! Set Debug Mode
!-------------------------------------------------------------------
call Open(arg)
call Debug(arg)
call Help(arg)
call Strict(arg)
call Close(arg)
call BeginSub(subname, version=version)
!-------------------------------------------------------------------
! NAMELIST file Setting
!-------------------------------------------------------------------
call nmlfile_init('init.nml')
!----------------------------------------------------------------
! Read init_nml
!----------------------------------------------------------------
call nmlfile_open(nmlunit, nmlreadable)
if (nmlreadable) then
read(nmlunit, nml=init_nml, iostat=nmlstat)
call DbgMessage('Stat of NAMELIST init_nml Input is <%d>', i=(/nmlstat/))
write(0, nml=init_nml)
else
call DbgMessage('Not Read NAMELIST init_nml')
call MessageNotify('W', subname, 'Can not Read NAMELIST init_nml. Force Use Default Value.')
end if
call nmlfile_close
!-------------------------------------------------------------------
! Initialize Dependent Modules
!-------------------------------------------------------------------
call constants_init ! 定数の設定
call grid_3d_init ! 空間格子点数の取得
call grid_wavenumber_init ! 波数格子点数の取得
call axis_x_init ! 座標軸データの生成
call axis_y_init ! 座標軸データの生成
call axis_z_init ! 座標軸データの生成
call io_gt4_out_init ! データ出力の初期設定
call constants_init ! 物理定数の取得
call spml_init ! SPMODEL の初期設定
!-------------------------------------------------------------------
! 軸データ生成
!-------------------------------------------------------------------
!!$ call axis_x_weight(x_Lon_Weight) ! 経度座標重みデータ取得
call axis_x_spectral(x_Lon) ! 経度座標データ取得
!!$ call axis_y_weight(y_Lat_Weight) ! 緯度座標重みデータ取得
call axis_y_spectral(y_Lat) ! 緯度座標データ取得
call axis_z_sigmahalf_manual ( z_Sigma , r_Sigma ) ! 半整数σレベル座標データ取得
!-------------------------------------------------------------------
! 出力用の軸データ設定
!-------------------------------------------------------------------
call io_gt4_out_SetDims(x_Lon) ! 経度座標重みデータ取得
!!$ call io_gt4_out_SetDims(x_Lon_Weight) ! 経度座標データ取得
call io_gt4_out_SetDims(y_Lat) ! 緯度座標重みデータ取得
!!$ call io_gt4_out_SetDims(y_Lat_Weight) ! 緯度座標データ取得
call io_gt4_out_SetDims(z_Sigma) ! 整数σレベル座標データ取得
call io_gt4_out_SetDims(r_Sigma) ! 半整数σレベル座標データ取得
!-------------------------------------------------------------------
! 出力用の変数データ設定
!-------------------------------------------------------------------
call io_gt4_out_SetVars('VelLon')
call io_gt4_out_SetVars('VelLat')
call io_gt4_out_SetVars('Vor')
call io_gt4_out_SetVars('Div')
call io_gt4_out_SetVars('Temp')
call io_gt4_out_SetVars('QVap')
call io_gt4_out_SetVars('Ps')
!-------------------------------------------------------------------
! Allocate Variable
!-------------------------------------------------------------------
allocate( xyz_VelLon(im,jm,km) )
allocate( xyz_VelLat(im,jm,km) )
allocate( xyz_Vor(im,jm,km) )
allocate( xyz_Div(im,jm,km) )
allocate( xyz_Temp(im,jm,km) )
allocate( xyz_QVap(im,jm,km) )
allocate( xyz_Ps(im,jm,km) )
allocate( xyz_Lon(im,jm,km) )
allocate( xyz_Lat(im,jm,km) )
allocate( xyz_Sigma(im,jm,km) )
allocate( xyr_Sigma(im,jm,km+1) )
allocate( wz_Psi((nm+1)*(nm+1), km) )
allocate( wz_Chi((nm+1)*(nm+1), km) )
!-------------------------------------------------------------------
! Set Average Value
!-------------------------------------------------------------------
xyz_VelLon = VelLonAve
xyz_VelLat = VelLatAve
xyz_Vor = VorAve
xyz_Div = DivAve
xyz_Temp = TempAve
xyz_QVap = QVapAve
xyz_Ps = PsAve
!-------------------------------------------------------------------
! Create 3-dimentional axis data
!-------------------------------------------------------------------
call HistoryAxisInquire(x_Lon % axisinfo, units=axis_units)
if ( StrHead( 'radians', trim(LChar(axis_units)) ) .or. StrHead( 'rad.', trim(LChar(axis_units)) ) ) then
RadDegFact = 1.
else
RadDegFact = 180. / pi
end if
do i = 1, im
xyz_Lon(i,:,:) = x_Lon%a_Dim(i) / RadDegFact
end do
call HistoryAxisInquire(y_Lat % axisinfo, units=axis_units)
if ( StrHead( 'radians', trim(LChar(axis_units)) ) .or. StrHead( 'rad.', trim(LChar(axis_units)) ) ) then
RadDegFact = 1.
else
RadDegFact = 180. / pi
end if
do j = 1, jm
xyz_Lat(:,j,:) = y_Lat%a_Dim(j) / RadDegFact
end do
do k = 1, km
xyz_Sigma(:,:,k) = z_Sigma%a_Dim(k)
end do
do k = 1, km + 1
xyr_Sigma(:,:,k) = r_Sigma%a_Dim(k)
end do
!!$ call DataDump('x_Lon1', x_Lon%a_Dim, strlen=70)
!!$ call DataDump('x_Lon2', x_Lon%a_Dim * RadDegFact , strlen=70)
!!$ call DataDump('x_Lon3', xyz_Lon / RadDegFact , strlen=70)
!!$ call DataDump('xyz_Lon', xyz_Lon, strlen=70)
!!$ call DataDump('xyz_Lat', xyz_Lat, strlen=70)
!!$ call DataDump('xyz_Sigma', xyz_Sigma, strlen=70)
!-------------------------------------------------------------------
! Set Various Conditions
!-------------------------------------------------------------------
if ( StriEq( 'rigid body rotation', trim(condition) ) ) then
! 剛体回転流
xyz_VelLon = xyz_VelLon + VelLonMax_rbr * cos( xyz_Lat )
VorDiv_Priority = .false.
call MessageNotify('M', subname, 'Selected <%c>.', c1=trim(LChar(condition)) )
!!$ call DataDump('xyz_Lat', xyz_Lat, strlen=80)
!!$ call DataDump('y_Lat', y_Lat%a_Dim, strlen=80)
!!$ call DataDump('xyz_VelLon', xyz_VelLon, strlen=80)
elseif ( StriEq( 'convex of surface pressure', trim(condition) ) ) then
! 度数をラジアンに変換するため
RadDegFact = 180. / pi
if (.not. Rad_Priority) then
LonLat_Radius_Rad = LonLat_Radius_Deg / RadDegFact
Lat_Center_Rad = Lat_Center_Deg / RadDegFact
Lon_Center_Rad = Lon_Center_Deg / RadDegFact
end if
xyz_Ps = xyz_Ps + PsMax * ( 1.0d0 - sqrt( min( LonLat_Radius_Rad**2 , ( xyz_Lon - Lon_Center_Rad )**2 + ( xyz_Lat - Lat_Center_Rad )**2 ) ) / LonLat_Radius_Rad )
call MessageNotify('M', subname, 'Selected <%c>.', c1=trim(LChar(condition)) )
elseif ( StriEq( 'convex of temperature', trim(condition) ) ) then
! 度数をラジアンに変換するため
RadDegFact = 180. / pi
if (.not. Rad_Priority) then
LonLat_Radius_Rad = LonLat_Radius_Deg / RadDegFact
Lat_Center_Rad = Lat_Center_Deg / RadDegFact
Lon_Center_Rad = Lon_Center_Deg / RadDegFact
end if
xyz_Temp = xyz_Temp + TempMax * ( 1.0D0 - sqrt( min( LonLat_Radius_Rad**2 , ( xyz_Lon - Lon_Center_Rad )**2 + ( xyz_Lat - Lat_Center_Rad )**2 ) ) / LonLat_Radius_Rad )
!!! xyz_VelLon = xyz_VelLon &
!!! & + TempMax &
!!! & * ( 1.0D0 &
!!! & - sqrt( &
!!! & min( LonLat_Radius_Rad**2 , &
!!! & ( xyz_Lon - Lon_Center_Rad )**2 &
!!! & + ( xyz_Lat - Lat_Center_Rad )**2 &
!!! & ) &
!!! & ) / LonLat_Radius_Rad &
!!! & )
!!!
!!! xyz_VelLat = xyz_VelLat &
!!! & + TempMax &
!!! & * ( 1.0D0 &
!!! & - sqrt( &
!!! & min( LonLat_Radius_Rad**2 , &
!!! & ( xyz_Lon - Lon_Center_Rad )**2 &
!!! & + ( xyz_Lat - Lat_Center_Rad )**2 &
!!! & ) &
!!! & ) / LonLat_Radius_Rad &
!!! & )
!!!
!!! xyz_QVap = xyz_QVap &
!!! & + TempMax &
!!! & * ( 1.0D0 &
!!! & - sqrt( &
!!! & min( LonLat_Radius_Rad**2 , &
!!! & ( xyz_Lon - Lon_Center_Rad )**2 &
!!! & + ( xyz_Lat - Lat_Center_Rad )**2 &
!!! & ) &
!!! & ) / LonLat_Radius_Rad &
!!! & )
call MessageNotify('M', subname, 'Selected <%c>.', c1=trim(LChar(condition)) )
elseif ( StriEq( 'agcm5.3 default value', trim(condition) ) ) then
do k = 1, km
do i = 1, im
do j = 1, jm
xyz_Temp(i,j,k) = xyz_Temp(i,j,k) + SIN ( REAL( i*j*(km-k) ) / REAL( im*jm*km )*10.0 ) * TPRTRB
end do
end do
end do
call MessageNotify('M', subname, 'Selected <%c>.', c1=trim(LChar(condition)) )
else
call MessageNotify('M', subname, 'Selected Isothermal and Nowind (Default).' )
end if
!-------------------------------------------------------------------
! Generate Vorticity and Divergence from Velocity
!-------------------------------------------------------------------
if ( .not. VorDiv_Priority ) then
xyz_Vor = xya_wa( wa_Div_xya_xya( xyz_VelLat , - xyz_VelLon ) / R0 )
xyz_Div = xya_wa( wa_Div_xya_xya( xyz_VelLon , xyz_VelLat ) / R0 )
end if
!-------------------------------------------------------------------
! Generate Velocity from Vorticity and Divergence
!-------------------------------------------------------------------
if ( VorDiv_Priority ) then
wz_Psi = wa_LaplaInv_wa( wa_xya( xyz_Vor ) ) * R0**2
wz_Chi = wa_LaplaInv_wa( wa_xya( xyz_Div ) ) * R0**2
xyz_VelLon = ( xya_GradLon_wa( wz_Chi ) - xya_GradLat_wa( wz_Psi ) ) / R0
xyz_VelLat = ( xya_GradLon_wa( wz_Psi ) + xya_GradLat_wa( wz_Chi ) ) / R0
end if
!-------------------------------------------------------------------
! データの出力 (t-Δt)
!-------------------------------------------------------------------
call io_gt4_out_Put('VelLon' , real(xyz_VelLon(:,:,:), REKIND) )
call io_gt4_out_Put('VelLat' , real(xyz_VelLat(:,:,:), REKIND) )
call io_gt4_out_Put('Vor' , real(xyz_Vor(:,:,:) , REKIND) )
call io_gt4_out_Put('Div' , real(xyz_Div(:,:,:) , REKIND) )
call io_gt4_out_Put('Temp' , real(xyz_Temp(:,:,:), REKIND) )
call io_gt4_out_Put('QVap' , real(xyz_QVap(:,:,:), REKIND) )
call io_gt4_out_Put('Ps' , real(xyz_Ps(:,:,1) , REKIND) )
!-------------------------------------------------------------------
! データの出力 (t)
!-------------------------------------------------------------------
call io_gt4_out_Put('VelLon' , real(xyz_VelLon(:,:,:), REKIND) )
call io_gt4_out_Put('VelLat' , real(xyz_VelLat(:,:,:), REKIND) )
call io_gt4_out_Put('Vor' , real(xyz_Vor(:,:,:) , REKIND) )
call io_gt4_out_Put('Div' , real(xyz_Div(:,:,:) , REKIND) )
call io_gt4_out_Put('Temp' , real(xyz_Temp(:,:,:), REKIND) )
call io_gt4_out_Put('QVap' , real(xyz_QVap(:,:,:), REKIND) )
call io_gt4_out_Put('Ps' , real(xyz_Ps(:,:,1) , REKIND) )
!-------------------------------------------------------------------
! 終了処理
!-------------------------------------------------------------------
call grid_3d_end ! 空間格子点数の取得
call grid_wavenumber_end ! 波数格子点数の取得
call axis_x_end ! 座標軸データの生成
call axis_y_end ! 座標軸データの生成
call axis_z_end ! 座標軸データの生成
call io_gt4_out_end ! データ出力の初期設定
call EndSub(subname)
end program init