!--------------------------------------------------------------------- ! Copyright (C) GFD Dennou Club, 2005. All rights reserved. !--------------------------------------------------------------------- !=begin != Program init ! ! * Developers: Morikawa Yasuhiro ! * Version: $Id: init.f90,v 1.3 2006/07/24 18:37:35 morikawa Exp $ ! * Tag Name: $Name: dcpam3-20060725 $ ! * Change History: ! !== Overview ! !Generate Various Initial Data ! !== Error Handling ! !== Known Bugs ! !== Note ! !== Future Plans ! ! !=end program init !=begin !== Dependency !----- 構造体参照モジュール ----- 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_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 !=end implicit none !=begin !== NAMELIST ! !初期値の種類の設定を行なう。 !現在 condition に与えて有効なのは以下の値である。 ! ! * (({ rigid body rotation })) ! * 剛体回転流を与える。風速の最大値は VelLonMax_rbr に与える。 ! * これを選択した場合、強制的に VorDiv_Priority は .false. ! に設定される。 ! ! * (({ convex of surface pressure })) ! * 地表面気圧の「山」を与える。位置とサイズは、度数で与える場合は ! Lon_Center_Dig, Lat_Center_Dig, LonLatRadius_Dig を用い、 ! ラジアンで与える場合は Lon_Center_Rad, Lat_Center_Rad, ! LonLatRadius_Rad を用いる (ただし、Rad_Priority を .true. ! にする必要がある)。 最大値は PsMax で与える。 ! ! * (({ convex of temperature })) ! * 温度の「山」を与える。位置とサイズは、度数で与える場合は ! Lon_Center_Dig, Lat_Center_Dig, LonLatRadius_Dig を用い、 ! ラジアンで与える場合は Lon_Center_Rad, Lat_Center_Rad, ! LonLatRadius_Rad を用いる (ただし、Rad_Priority を .true. ! にする必要がある)。 最大値は TempMax で与える。 ! ! * (({ agcm5.3 default value })) ! * 等温静止大気. 温度場に擾乱を与える. 既定値は以下. ! QVapAve = 1.0d-10 ! TempAve = 250.0d0 ! PsAve = 1.0d5 ! TPRTRB = 0.1 ! ! * その他 ! * 等温無風の初期値を与える。 ! !VelLonAve, VelLonAve, VelLatAve, VorAve, DivAve, TempAve, QVapAve, PsAve !には、それぞれ平均値を与える。 ! !デフォルトでは風速から渦度発散を生成するが、 !VorDiv_Priority を .true. にする事で、渦度発散から風速を生成する。 ! 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.174 ! 半径 (ラジアン) !=10.0度 real(DBKIND) :: Lat_Center_Deg = 45.0 ! 緯度の中心位置 (度数) real(DBKIND) :: Lon_Center_Deg = 100.0 ! 経度の中心位置 (度数) real(DBKIND) :: Lat_Center_Rad = 0.0 ! 緯度の中心位置 (ラジアン) !=0.0度 real(DBKIND) :: Lon_Center_Rad = 3.14 ! 経度の中心位置 (ラジアン)!=180度 logical :: Rad_Priority = .true. ! ラジアン表記を優先 ! !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 , & ! Maximum of 'Ps' & TempMax ! Maximum of 'Temp' !=end !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 !----- 作業用内部変数 ----- integer(INTKIND) :: i,j,k character(STRING) :: axis_units character(STRING), parameter:: version = & & '$Id: init.f90,v 1.3 2006/07/24 18:37:35 morikawa Exp $' character(STRING), parameter:: tagname = '$Name: dcpam3-20060725 $' character(STRING), parameter:: subname = "init" !------------------------------------------------------------------- ! Set Debug Mode !------------------------------------------------------------------- ! call SetDebug call BeginSub(subname) !---------------------------------------------------------------- ! Version identifier !---------------------------------------------------------------- call DbgMessage('%c :: %c', c1=trim(version), c2=trim(tagname)) !------------------------------------------------------------------- ! 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_Temp(:,:,1) = 280.000 xyz_Temp(:,:,2) = 277.368 xyz_Temp(:,:,3) = 274.737 xyz_Temp(:,:,4) = 272.105 xyz_Temp(:,:,5) = 269.474 xyz_Temp(:,:,6) = 266.842 xyz_Temp(:,:,7) = 264.211 xyz_Temp(:,:,8) = 261.579 xyz_Temp(:,:,9) = 258.947 xyz_Temp(:,:,10) = 256.316 xyz_Temp(:,:,11) = 253.684 xyz_Temp(:,:,12) = 251.053 xyz_Temp(:,:,13) = 248.421 xyz_Temp(:,:,14) = 245.789 xyz_Temp(:,:,15) = 243.158 xyz_Temp(:,:,16) = 240.526 xyz_Temp(:,:,17) = 237.895 xyz_Temp(:,:,18) = 235.263 xyz_Temp(:,:,19) = 232.632 xyz_Temp(:,:,20) = 230.000 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