program phy_radiation_flux_test
use phy_radiation_flux, only: PHYRADFLX, PhyRadFluxCreate, PhyRadFluxClose, PhyRadFluxPutLine, PhyRadFluxInitialized, RadiationFlux, RadiationDTempDt
use constants, only: CONST, Create, Get
use dc_test, only: AssertEqual, AssertGreaterThan, AssertLessThan
use dc_types, only: DP, STRING, TOKEN
use dc_string, only: StoA, PutLine, Printf, toChar
use dc_message, only: MessageNotify
use dc_args, only: ARGS, DCArgsOpen, DCArgsHelpMsg, DCArgsOption, DCArgsDebug, DCArgsHelp, DCArgsStrict, DCArgsClose
use dc_date_types, only: DC_DIFFTIME
use dc_date, only: DCDiffTimeCreate, EvalSec, EvalMin
use gt4_history, only: HistoryGet, HistoryCreate, HistoryAddVariable, HistoryAddAttr, HistoryPut, HistoryClose
implicit none
!---------------------------------------------------------
! 実験の表題, モデルの名称, 所属機関名
! Title of a experiment, name of model, sub-organ
!---------------------------------------------------------
character(*), parameter:: title = 'phy_radiation_flux_test $Name: dcpam4-20080624-1 $ :: ' // 'Test program of "phy_radiation_flux" module'
character(*), parameter:: source = 'dcpam4 ' // '(See http://www.gfd-dennou.org/library/dcpam)'
character(*), parameter:: institution = 'GFD Dennou Club (See http://www.gfd-dennou.org)'
!-------------------------------------------------------------------
! 格子点数・最大全波数
! Grid points and maximum truncated wavenumber
!-------------------------------------------------------------------
integer, parameter:: imax = 32
! 経度格子点数.
! Number of grid points in longitude
integer, parameter:: jmax = 16
! 緯度格子点数.
! Number of grid points in latitude
integer, parameter:: kmax = 12
! 鉛直層数.
! Number of vertical level
!-----------------------------------------------------------------
! 物理定数等の設定
! Configure physical constants etc.
!-----------------------------------------------------------------
type(CONST):: const_earth
real(DP):: PI ! $ \pi $ . 円周率. Circular constant
real(DP):: Grav ! $ g $ . 重力加速度. Gravitational acceleration
real(DP):: Cp ! $ C_p $ . 大気定圧比熱. Specific heat of air at constant pressure
real(DP):: StB ! $ \sigma_{SB} $ . ステファンボルツマン定数. Stefan-Boltzmann constant
!---------------------------------------------------------
! 軸データ
! Axes data
!---------------------------------------------------------
real(DP):: x_Lon (0:imax-1)
! 経度. Longitude
real(DP):: y_Lat (0:jmax-1)
! 緯度. Latitude
real(DP):: z_Sigma (0:kmax-1)
! $ \sigma $ レベル (整数).
! Full $ \sigma $ level
real(DP):: r_Sigma (0:kmax)
! $ \sigma $ レベル (半整数).
! Half $ \sigma $ level
!---------------------------------------------------------
! 時刻
! Time
!---------------------------------------------------------
type(DC_DIFFTIME):: delta_time
! $ \Delta t$ . 時刻. Time
character(TOKEN):: time_range
!---------------------------------------------------------
! 物理量
! Physical values
!---------------------------------------------------------
real(DP):: xy_SurfTemp (0:imax-1, 0:jmax-1)
! 地表面温度.
! Surface temperature
real(DP):: xy_SurfAlbedo (0:imax-1, 0:jmax-1)
! 地表アルベド.
! Surface albedo
real(DP):: xyz_Temp (0:imax-1, 0:jmax-1, 0:kmax-1)
! $ T $ . 温度. Temperature
real(DP):: xyz_QVap (0:imax-1, 0:jmax-1, 0:kmax-1)
! $ q $ . 比湿. Specific humidity
real(DP):: xyr_Press (0:imax-1, 0:jmax-1, 0:kmax)
! $ p_s $ . 地表面気圧 (半整数レベル).
! Surface pressure (half level)
!-----------------------------------
! フラックス
! Fluxes
real(DP):: xyr_RadLFlux (0:imax-1, 0:jmax-1, 0:kmax)
! 長波フラックス.
! Long wave flux
real(DP):: xya_SurfRadLMatrix (0:imax-1, 0:jmax-1, -1:1)
! $ T $ 陰解行列: 地表.
! $ T $ implicit matrix: surface
real(DP):: xyra_DelRadLFlux (0:imax-1, 0:jmax-1, 0:kmax, 0:1)
! 長波地表温度変化.
! Surface temperature tendency with long wave
real(DP):: xyr_RadSFlux (0:imax-1, 0:jmax-1, 0:kmax)
! 短波 (日射) フラックス.
! Short wave (insolation) flux
real(DP):: xyr_RadLFluxAns (0:imax-1, 0:jmax-1, 0:kmax)
! 長波フラックス.
! Long wave flux
real(DP):: xya_SurfRadLMatrixAns (0:imax-1, 0:jmax-1, -1:1)
! $ T $ 陰解行列: 地表.
! $ T $ implicit matrix: surface
real(DP):: xyra_DelRadLFluxAns (0:imax-1, 0:jmax-1, 0:kmax, 0:1)
! 長波地表温度変化.
! Surface temperature tendency with long wave
real(DP):: xyr_RadSFluxAns (0:imax-1, 0:jmax-1, 0:kmax)
! 短波 (日射) フラックス.
! Short wave (insolation) flux
!-----------------------------------
! 温度変換率
! Temperature tendency
real(DP):: xyz_DRadLTempDt (0:imax-1, 0:jmax-1, 0:kmax-1)
! 長波加熱率.
! Temperature tendency with long wave
real(DP):: xyz_DRadSTempDt (0:imax-1, 0:jmax-1, 0:kmax-1)
! 短波加熱率.
! Temperature tendency with short wave
real(DP):: xyz_DRadLTempDtAns (0:imax-1, 0:jmax-1, 0:kmax-1)
! 長波加熱率.
! Temperature tendency with long wave
real(DP):: xyz_DRadSTempDtAns (0:imax-1, 0:jmax-1, 0:kmax-1)
! 短波加熱率.
! Temperature tendency with short wave
!---------------------------------------------------------
! 作業変数
! Work variables
!---------------------------------------------------------
type(ARGS):: arg ! コマンドライン引数.
! Command line options
logical:: OPT_namelist ! -N, --namelist オプションの有無.
! Existence of '-N', '--namelist' option
character(STRING):: VAL_namelist
! -N, --namelist オプションの値.
! Value of '-N', '--namelist' option
type(PHYRADFLX):: phy_radflx00, phy_radflx01, phy_radflx02, phy_radflx03
logical:: err
integer:: i ! DO ループ用作業変数
! Work variables for DO loop
character(*), parameter:: subname = 'phy_radiation_flux_test'
continue
!---------------------------------------------------------
! コマンドライン引数の処理
! Command line options handling
!---------------------------------------------------------
call cmdline_optparse ! これは内部サブルーチン. This is an internal subroutine
!---------------------------------------------------------
! 物理定数の準備
! Prepare physical constants
!---------------------------------------------------------
call Create( const_earth ) ! (inout)
call Get( constant = const_earth, PI = PI, Grav = Grav, Cp = Cp, StB = StB ) ! (out)
!---------------------------------------------------------
! 緯度データの準備
! Prepare latitude data
!---------------------------------------------------------
call HistoryGet( file = 'phy_radiation_flux_test00.nc', varname = 'lon', array = x_Lon ) ! (out)
x_Lon = x_Lon * PI / 180.0_DP
call HistoryGet( file = 'phy_radiation_flux_test00.nc', varname = 'lat', array = y_Lat ) ! (out)
y_Lat = y_Lat * PI / 180.0_DP
call HistoryGet( file = 'phy_radiation_flux_test00.nc', varname = 'sig', array = z_Sigma ) ! (out)
call HistoryGet( file = 'phy_radiation_flux_test00.nc', varname = 'sigm', array = r_Sigma ) ! (out)
!---------------------------------------------------------
! 地表面データの準備
! Prepare surface data
!---------------------------------------------------------
xy_SurfTemp = 273.0_DP
xy_SurfAlbedo = 0.3_DP
!---------------------------------------------------------
! 時刻管理
! Time control
!---------------------------------------------------------
call DCDiffTimeCreate( diff = delta_time, value = 10.0_DP, unit = 'min' ) ! (in)
!---------------------------------------------------------
! 初期設定テスト
! Initialization test
!---------------------------------------------------------
call PhyRadFluxCreate( phy_radflx = phy_radflx00, imax = imax, jmax = jmax, kmax = kmax, y_Lat = y_Lat, Grav = Grav, Cp = Cp, StB = StB, xy_SurfTemp = xy_SurfTemp, xy_SurfAlbedo = xy_SurfAlbedo, current_time_value = 0.0_DP, current_time_unit = 'sec', delta_time_value = EvalMin( delta_time ), delta_time_unit = 'min', delta_time_value_RadLong = EvalMin( delta_time ) * 3, delta_time_unit_RadLong = 'min', delta_time_value_RadShort = EvalMin( delta_time ) * 2, delta_time_unit_RadShort = 'min' ) ! (in)
call AssertEqual( 'initialization test 1', answer = .true., check = PhyRadFluxInitialized(phy_radflx00) )
call PhyRadFluxPutLine( phy_radflx = phy_radflx00 ) ! (in)
!---------------------------------------------------------
! RadiationFlux, RadiationDelTemp テスト
! RadiationFlux, RadiationDelTemp test
!---------------------------------------------------------
call HistoryGet( file = 'phy_radiation_flux_test00.nc', varname = 'Temp', array = xyz_Temp ) ! (out)
call HistoryGet( file = 'phy_radiation_flux_test00.nc', varname = 'QVap', array = xyz_QVap ) ! (out)
call HistoryGet( file = 'phy_radiation_flux_test00.nc', varname = 'PressM', array = xyr_Press ) ! (out)
!!$ call output_init ! これは内部サブルーチン. This is an internal subroutine
do i = 1, 4
call RadiationFlux( phy_radflx = phy_radflx00, xyz_Temp = xyz_Temp, xyz_QVap = xyz_QVap, xyr_Press = xyr_Press, xyr_RadLFlux = xyr_RadLFlux, xya_SurfRadLMatrix = xya_SurfRadLMatrix, xyra_DelRadLFlux = xyra_DelRadLFlux, xyr_RadSFlux = xyr_RadSFlux ) ! (out)
call Printf(fmt = '')
call MessageNotify('M', subname, 'current_time=%d min', i = (/ i * 10 /) )
call PhyRadFluxPutLine( phy_radflx = phy_radflx00 ) ! (in)
call Printf(fmt = '')
call RadiationDTempDt( phy_radflx = phy_radflx00, xyr_RadLFlux = xyr_RadLFlux, xyr_RadSFlux = xyr_RadSFlux, xyr_Press = xyr_Press, xyz_DRadLTempDt = xyz_DRadLTempDt, xyz_DRadSTempDt = xyz_DRadSTempDt ) ! (out)
xyz_Temp = xyz_Temp * ( 1.0_DP + EvalSec ( delta_time ) * ( xyz_DRadLTempDt + xyz_DRadLTempDt ) )
xyz_QVap = xyz_QVap * 5.0_DP
!!$ call output_put ! これは内部サブルーチン. This is an internal subroutine
time_range = 'time=' // toChar( EvalMin( delta_time ) * ( i - 1 ) )
call HistoryGet( file = 'phy_radiation_flux_test01.nc', varname = 'RadLFlux', range = time_range, array = xyr_RadLFluxAns ) ! (out)
call AssertEqual( 'RadLFlux test 1-' // toChar( i ), answer = xyr_RadLFluxAns, check = xyr_RadLFlux, significant_digits = 13, ignore_digits = -15 )
call HistoryGet( file = 'phy_radiation_flux_test01.nc', varname = 'SurfRadLMatrix', range = time_range, array = xya_SurfRadLMatrixAns ) ! (out)
call AssertEqual( 'SurfRadLMatrix test 1-' // toChar( i ), answer = xya_SurfRadLMatrixAns, check = xya_SurfRadLMatrix, significant_digits = 13, ignore_digits = -15 )
call HistoryGet( file = 'phy_radiation_flux_test01.nc', varname = 'DelRadLFlux', range = time_range, array = xyra_DelRadLFluxAns ) ! (out)
call AssertEqual( 'DelRadLFlux test 1-' // toChar( i ), answer = xyra_DelRadLFluxAns, check = xyra_DelRadLFlux, significant_digits = 13, ignore_digits = -15 )
call HistoryGet( file = 'phy_radiation_flux_test01.nc', varname = 'RadSFlux', range = time_range, array = xyr_RadSFluxAns ) ! (out)
call AssertEqual( 'RadSFlux test 1-' // toChar( i ), answer = xyr_RadSFluxAns, check = xyr_RadSFlux, significant_digits = 14, ignore_digits = -15 )
call HistoryGet( file = 'phy_radiation_flux_test01.nc', varname = 'DRadLTempDt', range = time_range, array = xyz_DRadLTempDtAns ) ! (out)
call AssertEqual( 'DRadLTempDt test 1-' // toChar( i ), answer = xyz_DRadLTempDtAns, check = xyz_DRadLTempDt, significant_digits = 14, ignore_digits = -15 )
call HistoryGet( file = 'phy_radiation_flux_test01.nc', varname = 'DRadSTempDt', range = time_range, array = xyz_DRadSTempDtAns ) ! (out)
call AssertEqual( 'DRadSTempDt test 1-' // toChar( i ), answer = xyz_DRadSTempDtAns, check = xyz_DRadSTempDt, significant_digits = 14, ignore_digits = -15 )
end do
!!$ call output_close ! これは内部サブルーチン. This is an internal subroutine
!---------------------------------------------------------
! NAMELIST 読み込みテスト
! NAMELIST loading test
!---------------------------------------------------------
call PhyRadFluxCreate( phy_radflx = phy_radflx01, imax = imax, jmax = jmax, kmax = kmax, y_Lat = y_Lat, Grav = Grav, Cp = Cp, StB = StB, xy_SurfTemp = xy_SurfTemp, xy_SurfAlbedo = xy_SurfAlbedo, current_time_value = 0.0_DP, current_time_unit = 'sec', delta_time_value = EvalMin( delta_time ), delta_time_unit = 'min', nmlfile = VAL_namelist ) ! (in)
call AssertEqual( 'initialization test 2', answer = .true., check = PhyRadFluxInitialized(phy_radflx01) )
call PhyRadFluxPutLine( phy_radflx = phy_radflx01 ) ! (in)
!---------------------------------------------------------
! 終了処理テスト
! Termination test
!---------------------------------------------------------
call PhyRadFluxClose( phy_radflx = phy_radflx00 ) ! (inout)
call AssertEqual( 'termination test 1', answer = .false., check = PhyRadFluxInitialized(phy_radflx00) )
call PhyRadFluxPutLine( phy_radflx = phy_radflx00 ) ! (in)
call PhyRadFluxClose( phy_radflx = phy_radflx02, err = err ) ! (out)
call AssertEqual( 'termination test 2', answer = .true., check = err )
!---------------------------------------------------------
! 係数をゼロとした際の無効化チェック
! Check nullification by zeroing coefficient
!---------------------------------------------------------
call PhyRadFluxCreate( phy_radflx = phy_radflx03, imax = imax, jmax = jmax, kmax = kmax, y_Lat = y_Lat, Grav = Grav, Cp = Cp, StB = StB, xy_SurfTemp = xy_SurfTemp, xy_SurfAlbedo = xy_SurfAlbedo, LongBandNumber = 1, LongAbsorpCoeffQVap = (/ 0.0_DP /), LongAbsorpCoeffDryAir = (/ 0.0_DP /), LongBandWeight = (/ 1.0_DP /), ShortBandNumber = 1 , ShortAbsorpCoeffQVap = (/ 0.0_DP /), ShortAbsorpCoeffDryAir = (/ 0.0_DP /), ShortBandWeight = (/ 1.0_DP /), SolarCoeff = 0.0_DP ) ! (in)
call RadiationFlux( phy_radflx = phy_radflx03, xyz_Temp = xyz_Temp, xyz_QVap = xyz_QVap, xyr_Press = xyr_Press, xyr_RadLFlux = xyr_RadLFlux, xya_SurfRadLMatrix = xya_SurfRadLMatrix, xyra_DelRadLFlux = xyra_DelRadLFlux, xyr_RadSFlux = xyr_RadSFlux ) ! (out)
call RadiationDTempDt( phy_radflx = phy_radflx03, xyr_RadLFlux = xyr_RadLFlux, xyr_RadSFlux = xyr_RadSFlux, xyr_Press = xyr_Press, xyz_DRadLTempDt = xyz_DRadLTempDt, xyz_DRadSTempDt = xyz_DRadSTempDt ) ! (out)
xyz_DRadLTempDtAns = 0.0_DP
call AssertEqual( 'DRadLTempDt test 2', answer = xyz_DRadLTempDtAns, check = xyz_DRadLTempDt, significant_digits = 15, ignore_digits = -15 )
xyz_DRadSTempDtAns = 0.0_DP
call AssertEqual( 'DRadSTempDt test 2', answer = xyz_DRadSTempDtAns, check = xyz_DRadSTempDt, significant_digits = 15, ignore_digits = -15 )
call PhyRadFluxClose( phy_radflx = phy_radflx03 ) ! (out)
contains
subroutine output_init
call HistoryCreate( file = 'phy_radiation_flux_test01.nc', title = 'Radiation flux module test', source = 'dcpam4 (see http://www.gfd-dennou.org) ' , institution = 'GFD Dennou Club', dims = StoA('lon', 'lat', 'sig', 'sigm', 'd1', 'd2', 'time'), dimsizes = (/imax, jmax, kmax, kmax + 1, 3, 2, 0/), longnames = StoA('longitude', 'latitude', 'sigma at layer midpoints', 'sigma at layer end-points (half level)', 'dimension for SurfRadLMatrix', 'dimension for DelRadLFlux', 'time'), units = StoA('degree_east', 'degree_north', '1', '1', '1', '1', 'min'), origin = 0.0 , interval = real( EvalMin( delta_time ) ) ) ! (in)
call HistoryAddAttr( varname = 'lon', attrname = 'standard_name', value = 'longitude' ) ! (in)
call HistoryAddAttr( varname = 'lat', attrname = 'standard_name', value = 'latitude' ) ! (in)
call HistoryAddAttr( varname = 'sig', attrname = 'standard_name', value = 'atmosphere_sigma_coordinate' ) ! (in)
call HistoryAddAttr( varname = 'sigm', attrname = 'standard_name', value = 'atmosphere_sigma_coordinate' ) ! (in)
call HistoryAddAttr( varname = 'time', attrname = 'standard_name', value = 'time' ) ! (in)
call HistoryAddAttr( varname = 'sig', attrname = 'positive', value = 'down' ) ! (in)
call HistoryAddAttr( varname = 'sigm', attrname = 'positive', value = 'down' ) ! (in)
call HistoryPut( varname = 'lon', array = x_Lon / PI * 180.0_DP ) ! (in)
call HistoryPut( varname = 'lat', array = y_Lat / PI * 180.0_DP ) ! (in)
call HistoryPut( varname = 'sig', array = z_Sigma ) ! (in)
call HistoryPut( varname = 'sigm', array = r_Sigma ) ! (in)
call HistoryPut( varname = 'd1', array = (/-1, 0, 1/) ) ! (in)
call HistoryPut( varname = 'd2', array = (/0, 1/) ) ! (in)
call HistoryAddVariable( varname = 'RadLFlux', dims = StoA('lon', 'lat', 'sigm', 'time'), longname = 'long wave flux', units = '?', xtype = 'double' ) ! (in)
call HistoryAddVariable( varname = 'RadSFlux', dims = StoA('lon', 'lat', 'sigm', 'time'), longname = 'short wave flux', units = '?', xtype = 'double' ) ! (in)
call HistoryAddVariable( varname = 'SurfRadLMatrix', dims = StoA('lon', 'lat', 'd1', 'time'), longname = 'T implicit matrix: surface', units = '?', xtype = 'double' ) ! (in)
call HistoryAddVariable( varname = 'DelRadLFlux', dims = StoA('lon', 'lat', 'sigm', 'd2', 'time'), longname = 'surface temperature tendency with long wave', units = '?', xtype = 'double' ) ! (in)
call HistoryAddVariable( varname = 'DRadLTempDt', dims = StoA('lon', 'lat', 'sig', 'time'), longname = 'temperature tendency with long wave', units = 'K s-1', xtype = 'double' ) ! (in)
call HistoryAddVariable( varname = 'DRadSTempDt', dims = StoA('lon', 'lat', 'sig', 'time'), longname = 'temperature tendency with short wave', units = 'K s-1', xtype = 'double' ) ! (in)
end subroutine output_init
subroutine output_put
call HistoryPut( varname = 'RadLFlux', array = xyr_RadLFlux ) ! (in)
call HistoryPut( varname = 'SurfRadLMatrix', array = xya_SurfRadLMatrix ) ! (in)
call HistoryPut( varname = 'DelRadLFlux', array = xyra_DelRadLFlux ) ! (in)
call HistoryPut( varname = 'RadSFlux', array = xyr_RadSFlux ) ! (in)
call HistoryPut( varname = 'DRadLTempDt', array = xyz_DRadLTempDt ) ! (in)
call HistoryPut( varname = 'DRadSTempDt', array = xyz_DRadSTempDt ) ! (in)
end subroutine output_put
subroutine output_close
call HistoryClose
end subroutine output_close
subroutine cmdline_optparse
!
! コマンドライン引数の処理を行います
!
! Handle command line options
!
call DCArgsOpen( arg = arg ) ! (out)
call DCArgsHelpMsg( arg = arg, category = 'Title', msg = title ) ! (in)
call DCArgsHelpMsg( arg = arg, category = 'Usage', msg = './' // trim(subname) // ' [Options]' ) ! (in)
call DCArgsHelpMsg( arg = arg, category = 'Source', msg = source ) ! (in)
call DCArgsHelpMsg( arg = arg, category = 'Institution', msg = institution ) ! (in)
call DCArgsOption( arg = arg, options = StoA('-N', '--namelist'), flag = OPT_namelist, value = VAL_namelist, help = "Namelist filename") ! (in)
call DCArgsDebug( arg = arg ) ! (inout)
call DCArgsHelp( arg = arg ) ! (inout)
call DCArgsStrict( arg = arg ) ! (inout)
call DCArgsClose( arg = arg ) ! (inout)
end subroutine cmdline_optparse
end program phy_radiation_flux_test