subroutine dcpam_ape_physics( Dims, Vars_a )
!
! Physics main module
! DCPAM-APE Physics Main module
!
use dycore_type_mod, only: DYCORE_VARS, DYCORE_DIMS, STRING, INTKIND, REKIND, DBKIND
use grid_3d_mod, only: im, jm, km
use time_mod, only: DelTime
use spml_mod, only: wa_Div_xya_xya, xya_wa, wa_xya, xy_w, w_xy, xya_GradLon_wa,xya_GradLat_wa, wa_LaplaInv_wa
use constants_mod, only: R0
use io_gt4_out_mod , only: io_gt4_out_Put
use dc_trace, only: BeginSub, EndSub, DbgMessage
use physics_interpolate_mod, only: physics_interpolate_temp, physics_interpolate_geopot
use physics_negq_mod, only: physics_negq
use physics_lscond_mod , only: physics_lscond
use physics_cumulus_adjust_mod , only: physics_cumulus_adjust
use physics_dryadjust_mod , only: physics_dryadjust
use physics_ground_mod, only: physics_ground
use physics_radiation_main_mod, only: physics_radiation_main, physics_radiation_deltemp
use physics_verdiff_main_mod, only: physics_verdiff_main
use physics_surface_main_mod, only: physics_surface_main
use physics_implicit_mod, only: physics_implicit_init, physics_implicit_integrate, physics_implicit_fluxcorrection
use physics_ps_correction_mod, only: physics_ps_correction
use phys_vfilter_conserve_mod, only: phys_vfilter_conserve
use dycore_grid_mod, only: nm
implicit none
type(DYCORE_DIMS), intent(in) :: Dims ! 次元データ全種
type(DYCORE_VARS), intent(inout):: Vars_a ! 格子点データ全種(t+Δt)
character(STRING), parameter:: subname = "dcpam_ape_physics"
real(DBKIND) :: xyr_Temp(im,jm,km+1) ! 温度 (半整数)
real(DBKIND) :: xyz_Press(im,jm,km) ! 気圧 (整数)
real(DBKIND) :: xyr_Press(im,jm,km+1) ! 気圧 (半整数)
real(DBKIND) :: xyz_GeoPot(im,jm,km) ! ジオポテンシャル(整数)
real(DBKIND) :: xyr_GeoPot(im,jm,km+1) ! ジオポテンシャル(半整数)
real(DBKIND) :: xyz_DNegQvap1Dt(im,jm,km) ! 負の比湿除去 (比湿変化率)
real(DBKIND) :: xyz_DNegQvap2Dt(im,jm,km) ! 負の比湿除去 (比湿変化率)
real(DBKIND) :: xyz_DLscTempDt(im,jm,km) ! 大規模凝結による温度変化率
real(DBKIND) :: xyz_DLscQvapDt(im,jm,km) ! 大規模凝結による比湿変化率
real(DBKIND) :: xy_LscRain(im,jm) ! 大規模凝結による降水量
real(DBKIND) :: xyz_DCumulusTempDt(im,jm,km) ! 積雲スキームによる温度変化率
real(DBKIND) :: xyz_DCumulusQvapDt(im,jm,km) ! 積雲スキームによる比湿変化率
real(DBKIND) :: xy_CumulusRain(im,jm) ! 積雲スキームによる降水量
real(DBKIND) :: xyz_DDryTempDt(im,jm,km) ! 乾燥対流調節による温度変化率
real(DBKIND) :: xy_Rain(im,jm) ! cumulus + lsc 降水量
real(DBKIND) :: xy_Ps_b(im,jm) ! 地表面気圧
real(DBKIND) :: xyz_DRadLTempDt(im,jm,km) ! 長波加熱率
real(DBKIND) :: xyz_DRadSTempDt(im,jm,km) ! 短波加熱率
real(DBKIND) :: xy_SurfTemp(im,jm) ! 地表面温度
real(DBKIND) :: xy_SurfAlbedo(im,jm) ! 地表アルベド
real(DBKIND) :: xy_SurfHumidCoeff(im,jm) ! 地表湿潤度
real(DBKIND) :: xy_SurfRoughLength(im,jm) ! 地表粗度長
real(DBKIND) :: xy_SurfHeatCapacity(im,jm) ! 地表熱容量
real(DBKIND) :: xy_GroundTempFlux(im,jm) ! 地中熱フラックス
real(DBKIND) :: xyr_RadLFlux(im,jm,km+1) ! 長波フラックス
real(DBKIND) :: xyo_SurfRadLMatrix(im,jm,-1:1) ! Ｔ陰解行列：地表
real(DBKIND) :: xyro_DelRadLFlux(im,jm,km+1,0:1) ! 長波地表温度変化
real(DBKIND) :: xyr_RadSFlux(im,jm,km+1) ! 日射フラックス
real(DBKIND) :: xyr_VelLonFlux(im,jm,km+1) ! 速度経度成分フラックス
real(DBKIND) :: xyr_VelLatFlux(im,jm,km+1) ! 速度緯度成分フラックス
real(DBKIND) :: xyr_TempFlux(im,jm,km+1) ! 温度フラックス
real(DBKIND) :: xyr_QvapFlux(im,jm,km+1) ! 比湿フラックス
real(DBKIND) :: xyzo_VelMatrix(im,jm,km,-1:1) ! 速度陰解行列
real(DBKIND) :: xyzo_TempMatrix(im,jm,0:km,-1:1) ! 温度陰解行列
real(DBKIND) :: xyzo_QvapMatrix(im,jm,km,-1:1) ! 比湿陰解行列
real(DBKIND) :: xy_SurfVelMatrix(im,jm) ! 速度陰解行列: 地表
real(DBKIND) :: xyoo_SurfTempMatrix(im,jm,0:1,-1:1) ! 温度陰解行列: 地表
real(DBKIND) :: xyoo_SurfQvapMatrix(im,jm,0:1,-1:1) ! 比湿陰解行列: 地表
real(DBKIND) :: xyz_DVerdiffVelLonDt(im,jm,km) ! 経度成分 鉛直拡散加速度
real(DBKIND) :: xyz_DVerdiffVelLatDt(im,jm,km) ! 緯度成分 鉛直拡散加速度
real(DBKIND) :: xyz_DVerdiffTempDt(im,jm,km) ! 鉛直拡散加熱率
real(DBKIND) :: xyz_DVerdiffSurfTempDt(im,jm) ! 地表面 鉛直拡散加熱率
real(DBKIND) :: xyz_DVerdiffQvapDt(im,jm,km) ! 鉛直拡散加湿率
integer(INTKIND) :: xy_SurfCondition(im,jm) ! 地表状態
integer(INTKIND) :: k
real(DBKIND) :: wz_Psi_a((nm+1)*(nm+1), km)
real(DBKIND) :: wz_Chi_a((nm+1)*(nm+1), km)
real(DBKIND) :: xyz_DTempDtVertFiltCons(im,jm,km) ! 温度変化率(鉛直フィルター)
real(DBKIND) :: xyz_DVelLonDtVertFiltCons(im,jm,km) ! Ｕ変化率(鉛直フィルター)
real(DBKIND) :: xyz_DVelLatDtVertFilstCons(im,jm,km) ! Ｖ変化率(鉛直フィルター)
continue
! 開始処理
call BeginSub(subname)
! 1. 物理過程演算用の変数を初期化
xyz_DNegQvap1Dt = 0.0d0
xyz_DNegQvap2Dt = 0.0d0
xyz_DLscTempDt = 0.0d0
xyz_DLscQvapDt = 0.0d0
xy_LscRain = 0.0d0
xyz_DCumulusTempDt = 0.0d0
xyz_DCumulusQvapDt = 0.0d0
xy_CumulusRain = 0.0d0
xyz_DDryTempDt = 0.0d0
xy_Rain = 0.0d0
xy_Ps_b = Vars_a%xy_Ps
! 2. 地表条件設定
call physics_ground( xy_SurfTemp , xy_SurfAlbedo, xy_SurfHumidCoeff, xy_SurfRoughLength, xy_SurfCondition, xy_SurfHeatCapacity, xy_GroundTempFlux, Dims%y_Lat%a_Dim ) ! (in) 経度座標
!----------------------------------------------------------------
! 3. 温度半整数 sigma 補間, 気圧と高度の算出 (1)
!----------------------------------------------------------------
!----- 温度半整数 sigma 補間 -----
call physics_interpolate_temp( xyr_Temp , Vars_a%xyz_Temp , Dims%z_Sigma%a_Dim , Dims%r_Sigma%a_Dim ) ! (in) σレベル(半整数)座標
!----- 気圧と高度の算出 -----
call physics_interpolate_geopot( xyz_Press , xyr_Press , xyz_GeoPot , xyr_GeoPot , Vars_a%xy_Ps , Vars_a%xyz_Temp , xyr_Temp , Dims%z_Sigma%a_Dim , Dims%r_Sigma%a_Dim ) ! (in) σレベル(半整数)座標
!----------------------------------------------------------------
! 4. 負の水蒸気除去(1)
!----------------------------------------------------------------
call physics_negq( Vars_a%xyz_Qvap , xyz_DNegQvap1Dt , xyr_Press , 2.0d0*DelTime ) ! (in) 2Δt
!----------------------------------------------------------------
! 5. 湿潤過程 (積雲)
!----------------------------------------------------------------
call physics_cumulus_adjust( Vars_a%xyz_Temp , Vars_a%xyz_Qvap , xy_CumulusRain , xyz_DCumulusTempDt , xyz_DCumulusQvapDt , xyz_Press , xyr_Press , 2.0d0*DelTime ) ! (in) 2Δt
!----------------------------------------------------------------
! 6. 湿潤過程 (大規模凝結)
!----------------------------------------------------------------
call physics_lscond( Vars_a%xyz_Temp , Vars_a%xyz_Qvap , xy_LscRain , xyz_DLscTempDt , xyz_DLscQvapDt , xyz_Press , xyr_Press , 2.0d0*DelTime ) ! (in) 2Δt
!----------------------------------------------------------------
! 7. 負の水蒸気除去(2)
!----------------------------------------------------------------
call physics_negq( Vars_a%xyz_Qvap , xyz_DNegQvap1Dt , xyr_Press , 2.0d0*DelTime ) ! (in) 2Δt
!----------------------------------------------------------------
! 8. 温度半整数 sigma 補間, 気圧と高度の算出 (2)
!----------------------------------------------------------------
!----- Ps の計算しなおし -----
do k = 1, km
Vars_a%xy_Ps(:,:) = Vars_a%xy_Ps(:,:) + ( xyz_DLscQvapDt(:,:,k) + xyz_DCumulusQvapDt(:,:,k) + xyz_DNegQvap1Dt(:,:,k) ) * ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) * 2.0d0 * DelTime
end do
!----- 温度半整数 sigma 補間 -----
call physics_interpolate_temp( xyr_Temp , Vars_a%xyz_Temp , Dims%z_Sigma%a_Dim , Dims%r_Sigma%a_Dim ) ! (in) σレベル(半整数)座標
!----- 気圧と高度の算出 -----
call physics_interpolate_geopot( xyz_Press , xyr_Press , xyz_GeoPot , xyr_GeoPot , Vars_a%xy_Ps , Vars_a%xyz_Temp , xyr_Temp , Dims%z_Sigma%a_Dim , Dims%r_Sigma%a_Dim ) ! (in) σレベル(半整数)座標
!----------------------------------------------------------------
! 9. 陰解配列初期化
!----------------------------------------------------------------
call physics_implicit_init( xyr_VelLonFlux , xyr_VelLatFlux , xyr_TempFlux , xyr_QvapFlux , xyzo_VelMatrix , xyzo_TempMatrix , xyzo_QvapMatrix , xyr_Press , 2.0d0*DelTime , xy_SurfHeatCapacity , xy_SurfCondition ) ! (in) 地表状態
!----------------------------------------------------------------
! 10. 放射 flux
!----------------------------------------------------------------
call physics_radiation_main( xyr_RadLFlux , xyo_SurfRadLMatrix , xyro_DelRadLFlux , xyr_RadSFlux , Vars_a%xyz_Temp , xy_SurfTemp , Vars_a%xyz_Qvap , xyr_Press , Dims%x_Lon%a_Dim , Dims%y_Lat%a_Dim , xy_SurfAlbedo ) ! (in) 地表アルベド
!----------------------------------------------------------------
! 11. 鉛直拡散 flux
!----------------------------------------------------------------
call physics_verdiff_main( xyr_VelLonFlux , xyr_VelLatFlux , xyr_TempFlux , xyr_QvapFlux , xyzo_VelMatrix , xyzo_TempMatrix , xyzo_QvapMatrix , Vars_a%xyz_VelLon , Vars_a%xyz_VelLat , Vars_a%xyz_Temp , xyr_Temp , Vars_a%xyz_Qvap , xyz_Press , xyr_Press , xyz_GeoPot , xyr_GeoPot ) ! (in) 高度 (半整数)
!----------------------------------------------------------------
! 12. 地表 flux
!----------------------------------------------------------------
call physics_surface_main( xyr_VelLonFlux , xyr_VelLatFlux , xyr_TempFlux , xyr_QvapFlux , xy_SurfVelMatrix , xyoo_SurfTempMatrix , xyoo_SurfQvapMatrix , Vars_a%xyz_VelLon , Vars_a%xyz_VelLat , Vars_a%xyz_Temp , xyr_Temp , xy_SurfTemp , Vars_a%xyz_Qvap , xyz_Press , xyr_Press , xyz_GeoPot , xy_SurfHumidCoeff , xy_SurfRoughLength , xy_SurfCondition ) ! (in) 地表状態
!----------------------------------------------------------------
! 13. 時間変化率の計算 (implicit)
!----------------------------------------------------------------
call physics_implicit_integrate( xyz_DVerdiffVelLonDt , xyz_DVerdiffVelLatDt , xyz_DVerdiffTempDt , xyz_DVerdiffSurfTempDt , xyz_DVerdiffQvapDt , xyr_VelLonFlux , xyr_VelLatFlux , xyr_TempFlux , xyr_RadSFlux(:,:,1) , xyr_RadLFlux(:,:,1) , xy_GroundTempFlux , xyr_QvapFlux , xyzo_VelMatrix , xyzo_TempMatrix , xyzo_QvapMatrix , xy_SurfVelMatrix , xyoo_SurfTempMatrix , xyoo_SurfQvapMatrix , xyo_SurfRadLMatrix , 2.0d0*DelTime , xy_SurfCondition ) ! (in) 地表状態
! (2006-7-30 石渡) 追加
! フラックス補正
call physics_implicit_fluxcorrection ( xyr_VelLonFlux, xyr_VelLatFlux, xyr_TempFlux, xyr_QvapFlux, xyz_DVerdiffVelLonDt, xyz_DVerdiffVelLatDt, xyz_DVerdiffTempDt, xyz_DVerdiffSurfTempDt, xyz_DVerdiffQvapDt, xyzo_VelMatrix, xyzo_TempMatrix, xyzo_QvapMatrix, xy_SurfVelMatrix, xyoo_SurfTempMatrix, xyoo_SurfQvapMatrix, 2.0d0*DelTime )
!----------------------------------------------------------------
! 14. 放射による温度変化率
!----------------------------------------------------------------
do k = 1, km+1
xyr_RadLFlux(:,:,k) = xyr_RadLFlux(:,:,k) + (xyz_DVerdiffSurfTempDt(:,:) * xyro_DelRadLFlux(:,:,k,0) + xyz_DVerdiffTempDt(:,:,1) * xyro_DelRadLFlux(:,:,k,1) ) * 2.0d0*DelTime
end do
call physics_radiation_deltemp( xyz_DRadLTempDt , xyz_DRadSTempDt , xyr_RadLFlux , xyr_RadSFlux , xyr_Press ) ! (in) 圧力 (半整数)
!----------------------------------------------------------------
! 15. 温度変化分の足し込み
! (2006-7-28 石渡) フラックス補正してないからここで
! まとめてやっている.
!----------------------------------------------------------------
Vars_a%xyz_Temp = Vars_a%xyz_Temp + ( xyz_DRadLTempDt + xyz_DRadSTempDt ) * 2.0d0*DelTime
Vars_a%xyz_Temp = Vars_a%xyz_Temp + ( xyz_DVerdiffTempDt ) * 2.0d0* DelTime
Vars_a%xyz_Qvap = Vars_a%xyz_Qvap + ( xyz_DVerdiffQvapDt ) * 2.0d0* DelTime
Vars_a%xyz_VelLon = Vars_a%xyz_VelLon + ( xyz_DVerdiffVelLonDt ) * 2.0d0* DelTime
Vars_a%xyz_VelLat = Vars_a%xyz_VelLat + ( xyz_DVerdiffVelLatDt ) * 2.0d0* DelTime
! (2006-7-28 石渡) 追加.
call physics_integrate_surftemp( xy_SurfTemp, xyr_RadLFlux, xyro_DelRadLFlux, xyz_DVerdiffSurfTempDt, DelTime )
!----------------------------------------------------------------
! 16. 温度半整数 sigma 補間, 気圧と高度の算出 (3)
!----------------------------------------------------------------
!----- Ps の計算しなおし -----
do k = 1, km
Vars_a%xy_Ps(:,:) = Vars_a%xy_Ps(:,:) + xyz_DVerdiffQvapDt(:,:,k) * ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) * 2.0d0 * DelTime
end do
!----- 温度半整数 sigma 補間 -----
call physics_interpolate_temp( xyr_Temp , Vars_a%xyz_Temp , Dims%z_Sigma%a_Dim , Dims%r_Sigma%a_Dim ) ! (in) σレベル(半整数)座標
!----- 気圧と高度の算出 -----
call physics_interpolate_geopot( xyz_Press , xyr_Press , xyz_GeoPot , xyr_GeoPot , Vars_a%xy_Ps , Vars_a%xyz_Temp , xyr_Temp , Dims%z_Sigma%a_Dim , Dims%r_Sigma%a_Dim ) ! (in) σレベル(半整数)座標
!----------------------------------------------------------------
! 17. 乾燥対流調節
!----------------------------------------------------------------
call physics_dryadjust( Vars_a%xyz_Temp , xyz_DDryTempDt , xyz_Press , xyr_Press , 2.0d0*DelTime ) ! (in) 2Δt
!----------------------------------------------------------------
! 鉛直フィルター
! 2006-12-14 石渡追加
! デフォルトでコメントにしておく.
! 使用する時にはコメントを外すこと.
!----------------------------------------------------------------
! call phys_vfilter_conserve( &
! & Vars_a%xyz_Temp , & !(inout)
! & Vars_a%xyz_VelLon , & !(inout)
! & Vars_a%xyz_VelLat , & !(inout)
! & xyz_DTempDtVertFiltCons , & !(out)
! & xyz_DVelLonDtVertFiltCons , & !(out)
! & xyz_DVelLatDtVertFilstCons, & !(out)
! & xyr_Temp , & !(in)
! & xyr_Press , & !(in)
! & 2.0d0*DelTime ) !(in)
!----------------------------------------------------------------
! 18. 負の水蒸気除去(3)
!----------------------------------------------------------------
call physics_negq( Vars_a%xyz_Qvap , xyz_DNegQvap2Dt , xyr_Press , 2.0d0*DelTime ) ! (in) 2Δt
!----------------------------------------------------------------
! 19. Ps の計算しなおし
!----------------------------------------------------------------
do k = 1, km
Vars_a%xy_Ps(:,:) = Vars_a%xy_Ps(:,:) + xyz_DNegQvap2Dt(:,:,k) * ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) * 2.0d0 * DelTime
end do
! 凝結・降水による表面気圧変化
! このサブルーチンを使用する時は以下のコメントを外すこと.
! call physics_ps_correction( &
! & Vars_a%xyz_Qvap , & !(inout)
! & Vars_a%xy_Ps, & !(inout)
! & xyr_QvapFlux , & !(in)
! & xyz_DCumulusQvapDt, & !(in)
! & xyr_Press , & !(in)
! & 2.0d0 * DelTime ) !(in)
!----------------------------------------------------------------
! 20. 変数出力
!----------------------------------------------------------------
! call io_gt4_out_Put( 'GeoPot', real(xyz_GeoPot(:,:,:), DBKIND) )
! call io_gt4_out_Put( 'Press', real(xyz_Press(:,:,:), DBKIND) )
call io_gt4_out_Put( 'SurfTemp', real(xy_SurfTemp(:,:), DBKIND) )
call io_gt4_out_Put( 'DNegQvapDt', real( (xyz_DNegQvap1Dt(:,:,:) + xyz_DNegQvap2Dt(:,:,:)), DBKIND) )
call io_gt4_out_Put( 'DLscTempDt', real(xyz_DLscTempDt(:,:,:), DBKIND) )
call io_gt4_out_Put( 'DLscQvapDt', real(xyz_DLscQvapDt(:,:,:), DBKIND) )
call io_gt4_out_Put( 'LscRain', real(xy_LscRain(:,:), DBKIND) )
call io_gt4_out_Put( 'DCumulusTempDt', real(xyz_DCumulusTempDt(:,:,:), DBKIND) )
call io_gt4_out_Put( 'DCumulusQvapDt', real(xyz_DCumulusQvapDt(:,:,:), DBKIND) )
call io_gt4_out_Put( 'CumulusRain', real(xy_CumulusRain(:,:), DBKIND) )
call io_gt4_out_Put( 'DDryTempDt', real(xyz_DDryTempDt(:,:,:), DBKIND) )
xy_Rain = xy_LscRain + xy_CumulusRain
call io_gt4_out_Put( 'Rain', real(xy_Rain(:,:), DBKIND) )
call io_gt4_out_Put( 'DPsDt', real( (Vars_a%xy_Ps(:,:) - xy_Ps_b(:,:)), DBKIND) )
call io_gt4_out_Put( 'DRadLTempDt', real(xyz_DRadLTempDt(:,:,:), DBKIND) )
call io_gt4_out_Put( 'DRadSTempDt', real(xyz_DRadSTempDt(:,:,:), DBKIND) )
call io_gt4_out_Put('DVerdiffVelLonDt', real(xyz_DVerdiffVelLonDt(:,:,:), DBKIND) )
call io_gt4_out_Put('DVerdiffVelLatDt', real(xyz_DVerdiffVelLatDt(:,:,:), DBKIND) )
call io_gt4_out_Put('DVerdiffTempDt', real(xyz_DVerdiffTempDt(:,:,:), DBKIND) )
call io_gt4_out_Put('DVerdiffQvapDt', real(xyz_DVerdiffQvapDt(:,:,:), DBKIND) )
!-------------------------------------------------------------------
! 21. Generate Vorticity and Divergence from Velocity
!-------------------------------------------------------------------
Vars_a%xyz_Vor = xya_wa( wa_Div_xya_xya( Vars_a%xyz_VelLat,-Vars_a%xyz_VelLon )/R0)
Vars_a%xyz_Div = xya_wa( wa_Div_xya_xya( Vars_a%xyz_VelLon, Vars_a%xyz_VelLat ) /R0)
! wz_Psi_a = wa_LaplaInv_wa( wa_xya( Vars_a%xyz_Vor ) ) * R0**2
! wz_Chi_a = wa_LaplaInv_wa( wa_xya( Vars_a%xyz_Div ) ) * R0**2
! Vars_a%xyz_VelLon = ( xya_GradLon_wa( wz_Chi_a ) &
! & - xya_GradLat_wa( wz_Psi_a ) ) / R0
! Vars_a%xyz_VelLat = ( xya_GradLon_wa( wz_Psi_a ) &
! & + xya_GradLat_wa( wz_Chi_a ) ) / R0
! Vars_a%xyz_Temp = xya_wa( wa_xya(Vars_a%xyz_Temp) )
! Vars_a%xyz_Qvap = xya_wa( wa_xya(Vars_a%xyz_Qvap) )
! Vars_a%xy_Ps = xy_w( w_xy(Vars_a%xy_Ps) )
!----------------------------------------------------------------
! 終了処理
!----------------------------------------------------------------
call EndSub(subname)
end subroutine dcpam_ape_physics