Subroutine : |
|
xyr_VelLonFlux(im,jm,km+1) : | real(DBKIND), intent(inout)
|
xyr_VelLatFlux(im,jm,km+1) : | real(DBKIND), intent(inout)
|
xyr_TempFlux(im,jm,km+1) : | real(DBKIND), intent(inout)
|
xyr_QvapFlux(im,jm,km+1) : | real(DBKIND), intent(inout)
|
xy_SurfVelMatrix(im,jm) : | real(DBKIND), intent(out)
|
xyoo_SurfTempMatrix(im,jm,0:1,-1:1) : | real(DBKIND), intent(out)
|
xyoo_SurfQvapMatrix(im,jm,0:1,-1:1) : | real(DBKIND), intent(out)
|
xyz_VelLon(im,jm,km) : | real(DBKIND), intent(in)
|
xyz_VelLat(im,jm,km) : | real(DBKIND), intent(in)
|
xyz_Temp(im,jm,km) : | real(DBKIND), intent(in)
|
xyr_Temp(im,jm,km+1) : | real(DBKIND), intent(in)
|
xy_SurfTemp(im,jm) : | real(DBKIND), intent(in)
|
xyz_Qvap(im,jm,km) : | real(DBKIND), intent(in)
|
xyz_Press(im,jm,km) : | real(DBKIND), intent(in)
|
xyr_Press(im,jm,km+1) : | real(DBKIND), intent(in)
|
xyz_GeoPot(im,jm,km) : | real(DBKIND), intent(in)
|
xy_SurfHumidCoeff(im,jm) : | real(DBKIND), intent(in)
|
xy_SurfRoughLength(im,jm) : | real(DBKIND), intent(in)
|
xy_SurfCondition(im,jm) : | integer(INTKIND), intent(in)
|
subroutine physics_surface_main( xyr_VelLonFlux , xyr_VelLatFlux , xyr_TempFlux , xyr_QvapFlux , xy_SurfVelMatrix , xyoo_SurfTempMatrix , xyoo_SurfQvapMatrix , xyz_VelLon , xyz_VelLat , xyz_Temp , xyr_Temp , xy_SurfTemp , xyz_Qvap , xyz_Press , xyr_Press , xyz_GeoPot , xy_SurfHumidCoeff , xy_SurfRoughLength , xy_SurfCondition ) ! (in) 地表状態
!==== Dependency
use type_mod, only: REKIND, DBKIND, INTKIND, TOKEN, STRING
use grid_3d_mod, only: im, jm, km
use constants_mod, only: RAir, Cp, Grav, EL, ES0, RVap, EpsV
use physics_surface_coeff_mod, only: physics_surface_coeff
use dc_trace, only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump
implicit none
!==== Output
!
real(DBKIND), intent(out) :: xy_SurfVelMatrix(im,jm) , xyoo_SurfTempMatrix(im,jm,0:1,-1:1) , xyoo_SurfQvapMatrix(im,jm,0:1,-1:1) ! (out) 比湿陰解行列: 地表
!==== Input
!
real(DBKIND), intent(in) :: xyz_VelLon(im,jm,km) , xyz_VelLat(im,jm,km) , xyz_Temp(im,jm,km) , xyr_Temp(im,jm,km+1) , xy_SurfTemp(im,jm) , xyz_Qvap(im,jm,km) , xyz_Press(im,jm,km) , xyr_Press(im,jm,km+1) , xyz_GeoPot(im,jm,km) , xy_SurfHumidCoeff(im,jm) , xy_SurfRoughLength(im,jm) ! (in) 地表粗度長
integer(INTKIND), intent(in) :: xy_SurfCondition(im,jm) ! (in) 地表状態
!==== In/Out
!
real(DBKIND), intent(inout) :: xyr_VelLonFlux(im,jm,km+1) , xyr_VelLatFlux(im,jm,km+1) , xyr_TempFlux(im,jm,km+1) , xyr_QvapFlux(im,jm,km+1) ! (inout) 比湿フラックス
!----- 作業用内部変数 -----
character(STRING), parameter:: subname = "physics_surface_main"
! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
integer(INTKIND) :: i, j
! real(DBKIND), parameter :: RefPress = 100000.0D0 ! 参照気圧
real(DBKIND), parameter :: BasePotTemp = 300.0D0 ! 基本温位
real(DBKIND), parameter :: VelMinForRi = 0.01 ! Ri数用風最小値
real(DBKIND), parameter :: VelMinForVel = 0.01 ! 運動量用風最小値
real(DBKIND), parameter :: VelMinForTemp = 0.01 ! 熱用風最小値
real(DBKIND), parameter :: VelMinForQvap = 0.01 ! 水蒸気用風最小値
real(DBKIND), parameter :: VelMaxForVel = 1000. ! 運動量用風最小値
real(DBKIND), parameter :: VelMaxForTemp = 1000. ! 熱用風最小値
real(DBKIND), parameter :: VelMaxForQvap = 1000. ! 水蒸気用風最小値
real(DBKIND) :: xy_SurfBulkRiNum(im,jm) , xy_SurfTempTransCoeff(im,jm) , xy_SurfQvapTransCoeff(im,jm) , xy_SurfVelTransCoeff(im,jm) , xy_SurfTempBulkCoeff(im,jm) , xy_SurfQvapBulkCoeff(im,jm) , xy_SurfVelBulkCoeff(im,jm) , xy_SurfExner(im,jm) , xy_SurfVelAbs(im,jm) , xy_SurfSatQvap(im,jm) , xy_SurfDSatQvapDTemp(im,jm) ! 地表飽和比湿変化
continue
!----------------------------------------------------------------
! 開始処理
!----------------------------------------------------------------
call BeginSub(subname)
!----------------------------------------------------------------
! 地表面フラックス計算
!----------------------------------------------------------------
! ---- 1. 定数計算 ----
xy_SurfExner = ( xyz_Press(:,:,1) / xyr_Press(:,:,1) ) ** (RAir/Cp)
! ---- 2. バルクRi数 ----
do i = 1, im
do j = 1, jm
xy_SurfVelAbs(i,j) = SQRT ( xyz_VelLon(i,j,1)**2 + xyz_VelLat(i,j,1)**2 )
xy_SurfBulkRiNum(i,j) = Grav / BasePotTemp * ( xyz_Temp(i,j,1) / xy_SurfExner(i,j) - xy_SurfTemp(i,j) ) / MAX ( xy_SurfVelAbs(i,j) , VelMinForRi )**2 * xyz_GeoPot(i,j,1)
end do
end do
! ---- 3. バルク係数 ----
call physics_surface_coeff( xy_SurfVelBulkCoeff , xy_SurfTempBulkCoeff, xy_SurfQvapBulkCoeff, xy_SurfBulkRiNum , xy_SurfVelAbs , xy_SurfRoughLength , xy_SurfRoughLength , xyz_GeoPot(:,:,1) ) ! (in)
! ---- 5. 輸送係数 ----
do i = 1, im
do j = 1, jm
xy_SurfVelTransCoeff(i,j) = xy_SurfVelBulkCoeff(i,j) * xyr_Press(i,j,1) / (RAir * xyr_Temp(i,j,1) ) * MIN ( MAX (xy_SurfVelAbs(i,j) , VelMinForVel ) , VelMaxForVel )
xy_SurfTempTransCoeff(i,j) = xy_SurfTempBulkCoeff(i,j) * xyr_Press(i,j,1) / (RAir * xyr_Temp(i,j,1) ) * MIN ( MAX (xy_SurfVelAbs(i,j) , VelMinForTemp ) , VelMaxForTemp )
xy_SurfQvapTransCoeff(i,j) = xy_SurfQvapBulkCoeff(i,j) * xyr_Press(i,j,1) / (RAir * xyr_Temp(i,j,1) ) * MIN ( MAX (xy_SurfVelAbs(i,j) , VelMinForQvap ) , VelMaxForQvap )
end do
end do
! ---- 5. 飽和比湿----
xy_SurfSatQvap = EpsV * ES0 * EXP( EL / RVap * (1./273. - 1./xy_SurfTemp ) ) / xyr_Press(:,:,1)
xy_SurfDSatQvapDTemp = EL * xy_SurfSatQvap / ( RVap * xy_SurfTemp * xy_SurfTemp )
! ---- 6. フラックス----
xyr_VelLonFlux(:,:,1) = xyr_VelLonFlux(:,:,1) - xy_SurfVelTransCoeff * xyz_VelLon(:,:,1)
xyr_VelLatFlux(:,:,1) = xyr_VelLatFlux(:,:,1) - xy_SurfVelTransCoeff * xyz_VelLat(:,:,1)
xyr_TempFlux(:,:,1) = xyr_TempFlux(:,:,1) + Cp * xy_SurfTempTransCoeff * ( xy_SurfTemp - xyz_Temp(:,:,1) / xy_SurfExner )
xyr_QvapFlux(:,:,1) = xyr_QvapFlux(:,:,1) + EL * xy_SurfQvapTransCoeff * xy_SurfHumidCoeff * ( xy_SurfSatQvap - xyz_Qvap(:,:,1) )
! ---- 5. 陰解用行列 ----
xyoo_SurfTempMatrix = 0.0d0
xyoo_SurfQvapMatrix = 0.0d0
xy_SurfVelMatrix = xy_SurfVelTransCoeff
xyoo_SurfTempMatrix(:,:,1,0) = Cp * xy_SurfTempTransCoeff / xy_SurfExner
xyoo_SurfTempMatrix(:,:,0,1) = - Cp * xy_SurfTempTransCoeff / xy_SurfExner
xyoo_SurfQvapMatrix(:,:,1,0) = Cp * xy_SurfQvapTransCoeff * xy_SurfHumidCoeff
xyoo_SurfQvapMatrix(:,:,0,1) = - Cp * xy_SurfQvapTransCoeff * xy_SurfHumidCoeff
do i = 1, im
do j = 1, jm
if ( xy_SurfCondition(i,j) .GE. 1 ) then
xyoo_SurfTempMatrix(i,j,1,-1) = - Cp * xy_SurfTempTransCoeff(i,j)
xyoo_SurfTempMatrix(i,j,0,0) = Cp * xy_SurfTempTransCoeff(i,j)
xyoo_SurfQvapMatrix(i,j,1,-1) = - EL * xy_SurfQvapTransCoeff(i,j) * xy_SurfHumidCoeff(i,j) * xy_SurfDSatQvapDTemp(i,j)
xyoo_SurfQvapMatrix(i,j,0,0) = EL * xy_SurfQvapTransCoeff(i,j) * xy_SurfHumidCoeff(i,j) * xy_SurfDSatQvapDTemp(i,j)
end if
end do
end do
!----------------------------------------------------------------
! 終了処理
!----------------------------------------------------------------
call EndSub(subname)
end subroutine physics_surface_main