| Class | phy_implicit_sdh_V3 |
| In: |
phy_implicit/phy_implicit_sdh_V3.f90
|
Note that Japanese and English are described in parallel.
| PhyImplSDHTendency : | 時間変化率の計算 |
| PhyImplSDHSetMethodFromMatthews : | SurfType から計算法インデクスの作成 |
| PhyImplSDHInit : | 初期化 |
| ——————————- : | ———— |
| PhyImplSDHTendency : | Calculate tendency |
| PhyImplSDHSetMethodFromMatthews : | Set index for calculation method |
| PhyImplSDHInit : | Initialization |
| Subroutine : | |||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , intent(in )
| ||
| xy_SurfSnowFlux(0:imax-1, 1:jmax) : | real(DP), intent(in ) | ||
| xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(inout) |
Set index for calculation method from Matthews’ index
subroutine PhyImplSDHV3CorSOTempBySnowMelt( xy_IndexCalcMethod, xy_SurfSnowFlux, xy_SurfTemp )
!
!
!
! Set index for calculation method from Matthews' index
!
! モジュール引用 ; USE statements
!
! 時刻管理
! Time control
!
use timeset, only: DelTime ! $ \Delta t $ [s]
! 物理定数設定
! Physical constants settings
!
use constants, only: LatentHeatFusion
! $ L $ [J kg-1] .
! 融解の潜熱.
! Latent heat of fusion
! 雪と海氷の定数の設定
! Setting constants of snow and sea ice
!
use constants_snowseaice, only: SeaIceVolHeatCap, SeaIceThickness
! 宣言文 ; Declaration statements
!
integer , intent(in ) :: xy_IndexCalcMethod(0:imax-1, 1:jmax)
!
! Index for calculation method
real(DP), intent(in ) :: xy_SurfSnowFlux (0:imax-1, 1:jmax)
real(DP), intent(inout) :: xy_SurfTemp (0:imax-1, 1:jmax)
! 作業変数
! Work variables
!
real(DP) :: xy_SeaIceThickness(0:imax-1, 1:jmax)
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. phy_implicit_sdh_V3_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
!!$ if ( .not. FlagSSModel ) then
!!$ call MessageNotify( 'E', module_name, 'FlagSSModel has to be true.' )
!!$ end if
! FlagBucketModel は関係ないよね?
! SSModel 強制にした時点で, 水蒸気は地面と分離したから.
!!$ if ( .not. FlagBucketModel ) then
!!$ call MessageNotify( 'E', module_name, 'FlagBucketModel has to be true.' )
!!$ end if
!
! Set sea ice thickness
!
xy_SeaIceThickness = SeaIceThickness
!
! Set index for calculation method
!
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
case ( IndexSeaIce )
xy_SurfTemp(i,j) = xy_SurfTemp(i,j) + LatentHeatFusion * xy_SurfSnowFlux(i,j) * 2.0_DP * DelTime / ( SeaIceVolHeatCap * xy_SeaIceThickness(i,j) )
case ( IndexSlabOcean )
xy_SurfTemp(i,j) = xy_SurfTemp(i,j) + LatentHeatFusion * xy_SurfSnowFlux(i,j) * 2.0_DP * DelTime / SOHeatCapacity
case ( IndexPresTs )
case default
call MessageNotify( 'E', module_name, 'This index is inappropriate.' )
end select
end do
end do
end subroutine PhyImplSDHV3CorSOTempBySnowMelt
| Subroutine : | |||
| ArgFlagBucketModel : | logical , intent(in )
| ||
| ArgFlagSnow : | logical , intent(in )
| ||
| ArgFlagSlabOcean : | logical , intent(in )
| ||
| ArgFlagMajCompPhaseChange : | logical , intent(in )
| ||
| CondMajCompName : | character(*), intent(in ) |
This procedure input/output NAMELIST#phy_implicit_sdh_V3_nml .
subroutine PhyImplSDHV3Init( ArgFlagBucketModel, ArgFlagSnow, ArgFlagSlabOcean, ArgFlagMajCompPhaseChange, CondMajCompName )
!
! phy_implicit モジュールの初期化を行います.
! NAMELIST#phy_implicit_nml の読み込みはこの手続きで行われます.
!
! "phy_implicit" module is initialized.
! "NAMELIST#phy_implicit_nml" is loaded in this procedure.
!
! モジュール引用 ; USE statements
!
! NAMELIST ファイル入力に関するユーティリティ
! Utilities for NAMELIST file input
!
use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid
! ファイル入出力補助
! File I/O support
!
use dc_iounit, only: FileOpen
! 種別型パラメタ
! Kind type parameter
!
use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output
! 文字列操作
! Character handling
!
use dc_string, only: StoA
! 飽和比湿の算出
! Evaluate saturation specific humidity
!
use saturate, only: SaturateInit
! 主成分相変化
! Phase change of atmospheric major component
!
use saturate_major_comp, only : SaturateMajorCompInit
! 地表面フラックスユーティリティ
! Surface flux utility routines
!
use surface_flux_util, only : SurfaceFluxUtilInit
! 宣言文 ; Declaration statements
!
logical , intent(in ) :: ArgFlagBucketModel
! flag for use of bucket model
logical , intent(in ) :: ArgFlagSnow
! flag for treating snow
logical , intent(in ) :: ArgFlagSlabOcean
! flag for use of slab ocean
logical , intent(in ) :: ArgFlagMajCompPhaseChange
! flag for use of major component phase change
character(*), intent(in ) :: CondMajCompName
! 作業変数
! Work variables
!
integer:: unit_nml ! NAMELIST ファイルオープン用装置番号.
! Unit number for NAMELIST file open
integer:: iostat_nml ! NAMELIST 読み込み時の IOSTAT.
! IOSTAT of NAMELIST read
! NAMELIST 変数群
! NAMELIST group name
!
namelist /phy_implicit_sdh_V3_nml/ SOHeatCapacity, NumMaxItr, TempItrCrit, FlagSublimation
!
! デフォルト値については初期化手続 "phy_implicit#PhyImplInit"
! のソースコードを参照のこと.
!
! Refer to source codes in the initialization procedure
! "phy_implicit#PhyImplInit" for the default values.
!
! 実行文 ; Executable statement
!
if ( phy_implicit_sdh_V3_inited ) return
! Set flag for bucket model
FlagBucketModel = ArgFlagBucketModel
! Set flag for treating snow
FlagSnow = ArgFlagSnow
! Set flag for slab ocean
FlagSlabOcean = ArgFlagSlabOcean
! Set flag for major component phase change
FlagMajCompPhaseChange = ArgFlagMajCompPhaseChange
! デフォルト値の設定
! Default values settings
!
SOHeatCapacity = 4.187e3_DP * 1.0e3_DP * 60.0_DP
! 4.187d3 (J (kg K)-1) * 1.0d3 (kg m-3) * 60.0d0 (m)
NumMaxItr = 50
TempItrCrit = 0.05_DP
FlagSublimation = .false.
! NAMELIST の読み込み
! NAMELIST is input
!
if ( trim(namelist_filename) /= '' ) then
call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)
rewind( unit_nml )
read( unit_nml, nml = phy_implicit_sdh_V3_nml, iostat = iostat_nml ) ! (out)
close( unit_nml )
call NmlutilMsg( iostat_nml, module_name ) ! (in)
end if
! Initialization of modules used in this model
!
! 飽和比湿の算出
! Evaluate saturation specific humidity
!
call SaturateInit
if ( FlagMajCompPhaseChange ) then
! 主成分相変化
! Phase change of atmospheric major component
!
call SaturateMajorCompInit( CondMajCompName )
end if
! 地表面フラックスユーティリティ
! Surface flux utility routines
!
call SurfaceFluxUtilInit
! 印字 ; Print
!
call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
call MessageNotify( 'M', module_name, ' SOHeatCapacity = %f', d = (/ SOHeatCapacity /) )
call MessageNotify( 'M', module_name, ' NumMaxItr = %d', i = (/ NumMaxItr /) )
call MessageNotify( 'M', module_name, ' TempItrCrit = %f', d = (/ TempItrCrit /) )
call MessageNotify( 'M', module_name, ' FlagSublimation = %b', l = (/ FlagSublimation /) )
call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )
phy_implicit_sdh_V3_inited = .true.
end subroutine PhyImplSDHV3Init
| Subroutine : | |||
| xy_SurfCond(0:imax-1, 1:jmax) : | integer , intent(in )
| ||
| xy_SurfType(0:imax-1, 1:jmax) : | integer , intent(in )
| ||
| xy_SeaIceConc(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , intent(out)
|
Set index for calculation method from Matthews’ index
subroutine PhyImplSDHV3SetMethodMatthews( xy_SurfCond, xy_SurfType, xy_SeaIceConc, xy_IndexCalcMethod )
!
!
!
! Set index for calculation method from Matthews' index
!
! モジュール引用 ; USE statements
!
! 雪, 氷の割合
! snow/ice fraction
!
use snowice_frac, only : SeaIceAboveThreshold
! 宣言文 ; Declaration statements
!
integer , intent(in ) :: xy_SurfCond (0:imax-1, 1:jmax)
!
! Surface condition
integer , intent(in ) :: xy_SurfType (0:imax-1, 1:jmax)
! 土地利用.
! Surface index
real(DP), intent(in ) :: xy_SeaIceConc (0:imax-1, 1:jmax)
! 海氷密度 (0 <= xy_SeaIceConc <= 1)
! Sea ice concentration (0 <= xy_SeaIceConc <= 1)
integer , intent(out) :: xy_IndexCalcMethod(0:imax-1, 1:jmax)
!
! Index for calculation method
! 作業変数
! Work variables
!
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. phy_implicit_sdh_V3_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
!!$ if ( .not. FlagSSModel ) then
!!$ call MessageNotify( 'E', module_name, 'FlagSSModel has to be true.' )
!!$ end if
! FlagBucketModel は関係ないよね?
! SSModel 強制にした時点で, 水蒸気は地面と分離したから.
!!$ if ( .not. FlagBucketModel ) then
!!$ call MessageNotify( 'E', module_name, 'FlagBucketModel has to be true.' )
!!$ end if
!
! Set index for calculation method
!
do j = 1, jmax
do i = 0, imax-1
if ( xy_SurfCond(i,j) == 0 ) then
if ( xy_SurfType(i,j) >= 1 ) then
! land
xy_IndexCalcMethod(i,j) = IndexLandWithPresTs
else
! prescribed surface temperature
xy_IndexCalcMethod(i,j) = IndexPresTs
end if
else
if ( xy_SurfType(i,j) >= 1 ) then
! land
xy_IndexCalcMethod(i,j) = IndexLand
else
if ( SeaIceAboveThreshold( xy_SeaIceConc(i,j) ) ) then
! sea ice
xy_IndexCalcMethod(i,j) = IndexSeaIce
else if ( FlagSlabOcean ) then
! slab ocean
xy_IndexCalcMethod(i,j) = IndexSlabOcean
else
! prescribed surface temperature
!!$ xy_IndexCalcMethod(i,j) = IndexPresTs
call MessageNotify( 'E', module_name, 'Unexpected behavior.' )
end if
end if
end if
end do
end do
end subroutine PhyImplSDHV3SetMethodMatthews
| Subroutine : | |||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , intent(in)
| ||
| xy_BucketFlagOceanGrid(0:imax-1, 1:jmax) : | logical , intent(in)
| ||
| xy_SnowFrac(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xy_SurfH2OVapFlux(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_SurfLatentHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xyr_SoilHeatFlux(0:imax-1, 1:jmax, 0:kslmax) : | real(DP), intent(in)
| ||
| xyr_RadSFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_RadLFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_DeepSubSurfHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_TempB(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_SoilTemp(0:imax-1, 1:jmax, 1:kslmax) : | real(DP), intent(in)
| ||
| xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in) | ||
| xy_SurfHumidCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SoilHeatCap(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xy_SoilHeatDiffCoef(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(in)
| ||
| xyr_Press(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyz_Exner(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_Exner(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_VirTemp(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyz_Height(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_VelDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_QMixDiffCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_SurfVelTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfTempTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfQVapTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_SoilTempTransCoef(0:imax-1, 1:jmax, 0:kslmax) : | real(DP), intent(in)
| ||
| xy_SurfMajCompIceB(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SoilMoistB(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfSnowB(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_DUDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_DVDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(out)
| ||
| xy_DSurfTempDt(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xyz_DSoilTempDt(0:imax-1, 1:jmax, 1:kslmax) : | real(DP), intent(out)
| ||
| xy_DPsDt(0:imax-1, 1:jmax) : | real(DP), intent(out) | ||
| xy_DSurfMajCompIceDt(0:imax-1, 1:jmax) : | real(DP), intent(out) | ||
| xy_DSoilMoistDt(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_DSurfSnowDt(0:imax-1, 1:jmax) : | real(DP), intent(out)
|
時間変化率の計算を行います.
Calculate tendencies.
subroutine PhyImplSDHV3Tendency( xy_IndexCalcMethod, xy_BucketFlagOceanGrid, xy_SnowFrac, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xyz_TempB, xy_SurfTemp, xyz_SoilTemp, xyzf_QMixB, xy_SurfHumidCoef, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xyra_DelRadLFlux, xyr_Press, xyz_Exner, xyr_Exner, xyr_VirTemp, xyz_Height, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xyr_SoilTempTransCoef, xy_SurfMajCompIceB, xy_SoilMoistB, xy_SurfSnowB, xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DPsDt, xy_DSurfMajCompIceDt, xy_DSoilMoistDt, xy_DSurfSnowDt )
!
! 時間変化率の計算を行います.
!
! Calculate tendencies.
!
! モジュール引用 ; USE statements
!
! MPI 関連ルーチン
! MPI related routines
!
use mpi_wrapper, only: myrank, MPIWrapperFindMaxVal
! 座標データ設定
! Axes data settings
!
use axesset, only: r_SSDepth, z_SSDepth ! subsurface grid at midpoint of layer
! 物理定数設定
! Physical constants settings
!
use constants, only: Grav, CpDry, GasRDry, LatentHeat
! 雪と海氷の定数の設定
! Setting constants of snow and sea ice
!
use constants_snowseaice, only: TempCondWater, SnowVolHeatCap, SnowDens, SnowMaxThermDepth, SeaIceVolHeatCap , SeaIceThermCondCoef, SeaIceThickness, TempBelowSeaIce
! 飽和比湿の算出
! Evaluation of saturation specific humidity
!
use saturate, only: xy_CalcQVapSatOnLiq, xy_CalcQVapSatOnSol, xy_CalcDQVapSatDTempOnLiq, xy_CalcDQVapSatDTempOnSol
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! バケツモデル
! Bucket model
!
use Bucket_Model, only : BucketModEvapAndLatentHeatFlux
! 地表面フラックスユーティリティ
! Surface flux utility routines
!
use surface_flux_util, only : SurfaceFluxUtilLimitFlux
! 陰解法による時間積分のためのルーチン
! Routines for time integration with implicit scheme
!
use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3
! 宣言文 ; Declaration statements
!
integer , intent(in):: xy_IndexCalcMethod(0:imax-1, 1:jmax)
!
! Index for calculation method
logical , intent(in):: xy_BucketFlagOceanGrid(0:imax-1, 1:jmax)
!
! Flag for ocean grid point used in bucket model
real(DP), intent(in):: xy_SnowFrac (0:imax-1, 1:jmax)
!
! Snow fraction
real(DP), intent(in):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
! 東西方向運動量フラックス.
! Eastward momentum flux
real(DP), intent(in):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
! 南北方向運動量フラックス.
! Northward momentum flux
real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
! 比湿フラックス.
! Specific humidity flux
real(DP), intent(out):: xy_SurfH2OVapFlux(0:imax-1, 1:jmax)
! 惑星表面水蒸気フラックス.
! Water vapor flux at the surface
real(DP), intent(out):: xy_SurfLatentHeatFlux(0:imax-1, 1:jmax)
! 惑星表面潜熱フラックス.
! Latent heat flux at the surface
real(DP), intent(in):: xyr_SoilHeatFlux (0:imax-1, 1:jmax, 0:kslmax)
! 土壌中の熱フラックス (W m-2)
! Heat flux in sub-surface soil (W m-2)
real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
! 短波 (日射) フラックス.
! Shortwave (insolation) flux
real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
! 長波フラックス.
! Longwave flux
real(DP), intent(in):: xy_DeepSubSurfHeatFlux (0:imax-1, 1:jmax)
! 地中熱フラックス.
! "Deep subsurface heat flux"
! Heat flux at the bottom of surface/soil layer.
real(DP), intent(in):: xyz_TempB(0:imax-1, 1:jmax, 1:kmax)
! 温度.
! Temperature
real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
! 地表面温度.
! Surface temperature
real(DP), intent(in):: xyz_SoilTemp (0:imax-1, 1:jmax, 1:kslmax)
! 土壌温度 (K)
! Soil temperature (K)
real(DP), intent(in):: xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
!
!
real(DP), intent(in):: xy_SurfHumidCoef (0:imax-1, 1:jmax)
! 地表湿潤度.
! Surface humidity coefficient
real(DP), intent(in ):: xy_SoilHeatCap (0:imax-1, 1:jmax)
! 土壌熱容量 (J K-1 kg-1)
! Specific heat of soil (J K-1 kg-1)
real(DP), intent(in ):: xy_SoilHeatDiffCoef (0:imax-1, 1:jmax)
! 土壌熱伝導係数 (J m-3 K-1)
! Heat conduction coefficient of soil (J m-3 K-1)
real(DP), intent(in):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
! 長波地表温度変化.
! Surface temperature tendency with longwave
real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air pressure (half level)
real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
! Exner 関数 (整数レベル).
! Exner function (full level)
real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
! Exner 関数 (半整数レベル).
! Exner function (half level)
real(DP), intent(in):: xyr_VirTemp (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{T}_v $ . 仮温度 (半整数レベル).
! Virtual temperature (half level)
real(DP), intent(in):: xyz_Height (0:imax-1, 1:jmax, 1:kmax)
! 高度 (整数レベル).
! Height (full level)
real(DP), intent(in):: xyr_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(in):: xyr_TempDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in):: xyr_QMixDiffCoef (0:imax-1, 1:jmax, 0:kmax)
! 拡散係数：比湿.
! Diffusion coefficient: specific humidity
real(DP), intent(in):: xy_SurfVelTransCoef (0:imax-1, 1:jmax)
! 輸送係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in):: xy_SurfQVapTransCoef (0:imax-1, 1:jmax)
! 輸送係数：比湿.
! Transfer coefficient: specific humidity
real(DP), intent(in):: xyr_SoilTempTransCoef (0:imax-1, 1:jmax, 0:kslmax)
! 輸送係数：土壌温度.
! Transfer coefficient: soil temperature
real(DP), intent(in):: xy_SurfMajCompIceB (0:imax-1, 1:jmax)
!
! Surface major component ice amount.
real(DP), intent(in):: xy_SoilMoistB (0:imax-1, 1:jmax)
! 土壌水分.
! Soil moisture.
real(DP), intent(in):: xy_SurfSnowB (0:imax-1, 1:jmax)
! 積雪量.
! Surface snow amount.
real(DP), intent(out):: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{u}{t} $ . 東西風速変化.
! Eastward wind tendency
real(DP), intent(out):: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{v}{t} $ . 南北風速変化.
! Northward wind tendency
real(DP), intent(out):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{T}{t} $ . 温度変化.
! Temperature tendency
real(DP), intent(out):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ \DP{q}{t} $ . 質量混合比変化.
! Mass mixing ratio tendency
real(DP), intent(out):: xy_DSurfTempDt (0:imax-1, 1:jmax)
! 地表面温度変化率 (K s-1)
! Surface temperature tendency (K s-1)
real(DP), intent(out):: xyz_DSoilTempDt (0:imax-1, 1:jmax, 1:kslmax)
! $ \DP{Tg}{t} $ . 土壌温度変化 (K s-1)
! Temperature tendency (K s-1)
real(DP), intent(out):: xy_DPsDt (0:imax-1, 1:jmax)
real(DP), intent(out):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)
real(DP), intent(out):: xy_DSoilMoistDt (0:imax-1, 1:jmax)
! 土壌温度時間変化率 (kg m-2 s-1)
! Soil temperature tendency (kg m-2 s-1)
real(DP), intent(out):: xy_DSurfSnowDt (0:imax-1, 1:jmax)
! 積雪率時間変化率 (kg m-2 s-1)
! Surface snow amount tendency (kg m-2 s-1)
! 作業変数
! Work variables
!
real(DP) :: xy_SurfHeatCapacity (0:imax-1, 1:jmax)
! 地表熱容量.
! Surface heat capacity
real(DP) :: xy_SeaIceThickness(0:imax-1, 1:jmax)
!
! Sea ice thickness
real(DP) :: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：運動量.
! Transfer coefficient: velocity
real(DP) :: xyr_TempTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP) :: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：質量.
! Transfer coefficient: mass of constituents
real(DP):: xyza_UVMtx (0:imax-1, 1:jmax, 1:kmax, -1:1)
! 速度陰解行列.
! Implicit matrix about velocity
real(DP):: xyz_UVec (0:imax-1, 1:jmax, 1:kmax)
! 速度陰解ベクトル.
! Implicit vector about velocity
real(DP):: xyz_VVec (0:imax-1, 1:jmax, 1:kmax)
! 速度陰解ベクトル.
! Implicit vector about velocity
real(DP):: xyza_TempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
! 温度陰解行列.
! Implicit matrix about temperature
real(DP):: xyz_TempVec(0:imax-1, 1:jmax, 1:kmax)
! 温度陰解ベクトル.
! Implicit vector about temperature
real(DP):: xyza_QMixMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
! 質量混合比陰解行列.
! Implicit matrix about mass mixing ratio
real(DP):: xyzf_QMixVec(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! 質量混合比陰解ベクトル.
! Implicit vector about mass mixing ratio
real(DP):: xyaa_SurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
! 惑星表面エネルギー収支用陰解行列
! Implicit matrix for surface energy balance
real(DP):: xy_SurfRH(0:imax-1,1:jmax)
real(DP):: xyza_UVLUMtx (0:imax-1, 1:jmax, 1:kmax,-1:1)
! LU 行列.
! LU matrix
!!$ real(DP):: xyza_TempQVapLUMtx (0:imax-1, 1:jmax, -kmax:kmax, -1:1)
!!$ ! LU 行列.
!!$ ! LU matrix
!!$ real(DP):: xyz_DelTempQVap (0:imax-1, 1:jmax, -kmax:kmax)
!!$ ! $ T q $ の時間変化.
!!$ ! Tendency of $ T q $
!!$
!!$ real(DP):: xyza_TempLUMtx (0:imax-1, 1:jmax, 0:kmax, -1:1)
!!$ ! LU 行列.
!!$ ! LU matrix
!!$ real(DP):: xyz_DelTempLUVec (0:imax-1, 1:jmax, 0:kmax)
!!$ ! $ T q $ の時間変化.
!!$ ! Tendency of $ T q $
real(DP):: xyza_QMixLUMtx (0:imax-1, 1:jmax, 1:kmax, -1:1)
! LU 行列.
! LU matrix
real(DP):: xyz_DelQMixLUVec (0:imax-1, 1:jmax, 1:kmax)
! $ q $ の時間変化.
! Tendency of $ q $
real(DP):: xyaa_SoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
! 土壌温度拡散方程式の行列
! Matrix for diffusion equation of soil temperature
real(DP):: xya_SoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
! 土壌温度拡散方程式のベクトル
! Vector for diffusion equation of soil temperature
real(DP):: xyaa_TempSoilTempLUMtx (0:imax-1, 1:jmax, -kslmax:kmax, -1:1)
! LU 行列.
! LU matrix
real(DP):: xya_DelTempSoilTempLUVec (0:imax-1, 1:jmax, -kslmax:kmax)
! $ T, Tg $ の時間変化.
! Tendency of $ T $ and $ Tg |
real(DP):: SurfSnowATentative
! 積雪量の仮の値 (kg m-2)
! pseudo value of surface snow amount (kg m-2)
real(DP):: xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax)
!
! Latent heat flux by major component ice sublimation
! (variable only for debug)
real(DP):: xy_LatHeatFluxBySnowMelt(0:imax-1, 1:jmax)
!
! Latent heat flux by melt
! (variable only for debug)
real(DP):: xy_LatHeatFluxBySeaIceMelt(0:imax-1, 1:jmax)
!
! Latent heat flux by sea ice melt
! (variable only for debug)
real(DP):: xy_LatHeatFluxByOtherSpc(0:imax-1, 1:jmax)
real(DP):: xy_DAtmMassDt(0:imax-1, 1:jmax)
real(DP):: xy_SurfQVapSatOnLiq(0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfQVapSatOnSol (0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfQVapSat (0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfDQVapSatDTempOnLiq (0:imax-1, 1:jmax)
! 地表飽和比湿変化.
! Saturated specific humidity tendency on surface
real(DP):: xy_SurfDQVapSatDTempOnSol (0:imax-1, 1:jmax)
! 地表飽和比湿変化.
! Saturated specific humidity tendency on surface
real(DP):: xy_SurfDQVapSatDTemp (0:imax-1, 1:jmax)
! 地表飽和比湿変化.
! Saturated specific humidity tendency on surface
real(DP):: xy_SurfSoilHeatFlux(0:imax-1, 1:jmax)
! 惑星表面土壌熱伝導フラックス.
! Soil heat conduction flux at the surface
real(DP):: xyz_TempSave (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyz_TempA (0:imax-1, 1:jmax, 1:kmax)
real(DP):: MaxTempInc
real(DP):: xy_SurfTempSave(0:imax-1, 1:jmax)
real(DP):: xy_SurfTempA (0:imax-1, 1:jmax)
real(DP):: MaxSurfTempInc
real(DP):: a_LocalMax (2)
real(DP):: a_GlobalMax(2)
integer:: iitr
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: l ! 行列用 DO ループ用作業変数
! Work variables for DO loop of matrices
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. phy_implicit_sdh_V3_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )
! Check kskmax
!
if ( kslmax < 1 ) then
call MessageNotify( 'E', module_name, 'kslmax is less than 1.' )
end if
!!$ if ( .not. FlagSSModel ) then
!!$ call MessageNotify( 'E', module_name, 'FlagSSModel has to be true.' )
!!$ end if
! FlagBucketModel は関係ないよね?
! SSModel 強制にした時点で, 水蒸気は地面と分離したから.
!!$ if ( .not. FlagBucketModel ) then
!!$ call MessageNotify( 'E', module_name, 'FlagBucketModel has to be true.' )
!!$ end if
! Set heat capacity, tentatively here.
!
!!$ xy_SurfHeatCapacity = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! J K-1 m-3 kg m-2 (kg m-3)-1 = J K-1 m-5 m3 = J K-1 m-2
xy_SurfHeatCapacity(i,j) = SnowVolHeatCap * min( max( xy_SurfSnowB(i,j) / SnowDens, 0.0_DP ), SnowMaxThermDepth )
case default
xy_SurfHeatCapacity(i,j) = 0.0_DP
end select
end do
end do
!
! Set sea ice thickness
!
xy_SeaIceThickness = SeaIceThickness
! Set coefficients for soil heat flux
!
!!$ xyr_SoilTempTransCoef = xyr_BareSoilTempTransCoef
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ if ( xy_SurfSnowB(i,j) >
!!$ xyr_SoilTempTransCoef = xyr_BareSoilTempTransCoef
! 輸送係数の計算
! Calculate transfer coefficient
!
xyr_VelTransCoef (:,:,0) = 0.0_DP
xyr_VelTransCoef (:,:,kmax) = 0.0_DP
xyr_TempTransCoef(:,:,0) = 0.0_DP
xyr_TempTransCoef(:,:,kmax) = 0.0_DP
xyr_QMixTransCoef(:,:,0) = 0.0_DP
xyr_QMixTransCoef(:,:,kmax) = 0.0_DP
do k = 1, kmax-1
xyr_VelTransCoef(:,:,k) = xyr_VelDiffCoef(:,:,k) * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
xyr_TempTransCoef(:,:,k) = xyr_TempDiffCoef(:,:,k) * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
xyr_QMixTransCoef(:,:,k) = xyr_QMixDiffCoef(:,:,k) * xyr_Press(:,:,k) / ( GasRDry * xyr_VirTemp(:,:,k) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
end do
! Calculation for momentum diffusion
!
call PhyImplSDHV3TendencyMomCore( xyr_MomFluxX, xyr_MomFluxY, xyr_Press, xyr_VelTransCoef, xy_SurfVelTransCoef, xyz_DUDt, xyz_DVDt )
! Calculation for thermal diffusion
!
! 飽和比湿の計算
! Calculate saturated specific humidity
!
xy_SurfQVapSatOnLiq = xy_CalcQVapSatOnLiq ( xy_SurfTemp, xyr_Press(:,:,0) )
xy_SurfQVapSatOnSol = xy_CalcQVapSatOnSol ( xy_SurfTemp, xyr_Press(:,:,0) )
xy_SurfQVapSat = ( 1.0_DP - xy_SnowFrac ) * xy_SurfQVapSatOnLiq + xy_SnowFrac * xy_SurfQVapSatOnSol
xy_SurfDQVapSatDTempOnLiq = xy_CalcDQVapSatDTempOnLiq( xy_SurfTemp, xy_SurfQVapSatOnLiq )
xy_SurfDQVapSatDTempOnSol = xy_CalcDQVapSatDTempOnSol( xy_SurfTemp, xy_SurfQVapSatOnSol )
xy_SurfDQVapSatDTemp = ( 1.0_DP - xy_SnowFrac ) * xy_SurfDQVapSatDTempOnLiq + xy_SnowFrac * xy_SurfDQVapSatDTempOnSol
! Initialization
!
xyzf_DQMixDt = 0.0_DP
xy_DSurfTempDt = 0.0_DP
xyz_DSoilTempDt = 0.0_DP
xyz_TempSave = xyz_TempB
xy_SurfTempSave = xy_SurfTemp
! iteration
iitr = 1
loop_itr : do
! Tendencies of atmospheric and surface temperatures, and atmospheric
! water vapor are solved with a fixed surface heat conduction flux at the surface.
! Obtained tendencies of surface temperature and atmospheric water vapor will be
! used below to estimate surface water vapor flux.
!
xy_SurfSoilHeatFlux = xyr_SoilHeatFlux(:,:,0) - xyr_SoilTempTransCoef(:,:,0) * ( xyz_DSoilTempDt(:,:,1) - xy_DSurfTempDt ) * ( 2.0_DP * DelTime )
!
call PhyImplSDHV3TendencyHeatTQCore( xy_IndexCalcMethod, xy_SeaIceThickness, xy_SnowFrac, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfSoilHeatFlux, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xy_SurfHumidCoef, xy_SurfHeatCapacity, xyra_DelRadLFlux, xyr_Press, xyz_Exner, xyr_Exner, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt )
! Tendencies of atmospheric, surface, and subsurface temperatures, and
! atmospheric water vapor are solved with a fixed surface water vapor flux.
!
n = IndexH2OVap
xy_SurfH2OVapFlux = xyrf_QMixFlux(:,:,0,n) - xy_SurfHumidCoef * xy_SurfQVapTransCoef * ( xyzf_DQMixDt(:,:,1,n) - xy_SurfDQVapSatDTemp * xy_DSurfTempDt ) * 2.0_DP * DelTime
! Limit surface flux not to be negative atmospheric content
! IMPORTANT : Now, only the water vapor flux is restricted.
call SurfaceFluxUtilLimitFlux( ( 2.0_DP * DelTime ), xyzf_QMixB, xyr_Press, xy_SurfH2OVapFlux )
! Calculation of latent heat flux
xy_SurfLatentHeatFlux = LatentHeat * xy_SurfH2OVapFlux
!
if ( FlagBucketModel ) then
! バケツモデルのための地表面フラックス修正
! Modification of surface flux for bucket model
!
call BucketModEvapAndLatentHeatFlux( xy_BucketFlagOceanGrid, xy_SoilMoistB, xy_SurfSnowB, xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux )
end if
!
call PhyImplSDHV3TendencyHeatCore( xy_IndexCalcMethod, xy_SeaIceThickness, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfHumidCoef, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xyra_DelRadLFlux, xyr_Press, xyz_Exner, xyr_Exner, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xyr_SoilTempTransCoef, xy_SurfMajCompIceB, xy_SoilMoistB, xy_SurfSnowB, xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DPsDt, xy_DSurfMajCompIceDt, xy_DSoilMoistDt, xy_DSurfSnowDt )
! Check
xyz_TempA = xyz_TempB + xyz_DTempDt * 2.0_DP * DelTime
xy_SurfTempA = xy_SurfTemp + xy_DSurfTempDt * DelTime
MaxTempInc = maxval( abs( xyz_TempA - xyz_TempSave ) )
MaxSurfTempInc = maxval( abs( xy_SurfTempA - xy_SurfTempSave ) )
a_LocalMax(1) = MaxTempInc
a_LocalMax(2) = MaxSurfTempInc
call MPIWrapperFindMaxVal( 2, a_LocalMax, a_GlobalMax )
MaxTempInc = a_GlobalMax(1)
MaxSurfTempInc = a_GlobalMax(2)
!!$ if ( myrank == 0 ) then
!!$ call MessageNotify( 'M', module_name, &
!!$ & 'Itr: %d : dT = %f, dTs = %f', &
!!$ & i = (/iitr/), d = (/ MaxTempInc, MaxSurfTempInc /) )
!!$ end if
if ( ( MaxTempInc <= TempItrCrit ) .and. ( MaxSurfTempInc <= TempItrCrit ) ) then
exit loop_itr
end if
xyz_TempSave = xyz_TempA
xy_SurfTempSave = xy_SurfTempA
iitr = iitr + 1
if ( iitr > NumMaxItr ) then
if ( NumMaxItr > 2 ) then
if ( myrank == 0 ) then
call MessageNotify( 'M', module_name, 'Too many iterations, Itr: %d : dT = %f, dTs = %f', i = (/iitr/), d = (/ MaxTempInc, MaxSurfTempInc /) )
end if
end if
exit loop_itr
end if
end do loop_itr
!!$ if ( myrank == 0 ) then
!!$ call MessageNotify( 'M', module_name, &
!!$ & 'Itr: %d : dT = %f, dTs = %f', &
!!$ & i = (/iitr/), d = (/ MaxTempInc, MaxSurfTempInc /) )
!!$ end if
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine PhyImplSDHV3Tendency
| Variable : | |||
| FlagMajCompPhaseChange : | logical, save
|
| Constant : | |||
| IndexLandWithPresTs = 14 : | integer, parameter
|
| Subroutine : | |||||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , intent(in)
| ||||
| xy_SeaIceThickness(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||||
| xyr_Press(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||||
| xyz_Exner(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||||
| xyr_Exner(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||||
| xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||||
| xy_SurfHeatCapacity(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||||
| xy_SoilHeatCap(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||||
| xy_SurfH2OVapFlux(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||||
| xy_SurfLatentHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||||
| xyr_SoilHeatFlux(0:imax-1, 1:jmax, 0:kslmax) : | real(DP), intent(in)
| ||||
| xy_SurfTempTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||||
| xyr_SoilTempTransCoef(0:imax-1, 1:jmax, 0:kslmax) : | real(DP), intent(in)
| ||||
| xyr_RadSFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||||
| xyr_RadLFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||||
| xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(in)
| ||||
| xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||||
| xy_LatHeatFluxBySnowMelt(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||||
| xy_LatHeatFluxBySeaIceMelt(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||||
| xy_DeepSubSurfHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||||
| xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in )
| ||||
| xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in )
| ||||
| xy_DSurfTempDt(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||||
| xyz_DSoilTempDt(0:imax-1, 1:jmax, 1:kslmax) : | real(DP), intent(in )
| ||||
| xy_DSoilMoistDt(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||||
| xy_DSurfSnowDt(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||||
| xy_DPsDt(0:imax-1, 1:jmax) : | real(DP), intent(in ) | ||||
| xy_DSurfMajCompIceDt(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
A part of conservation of energy is checked.
subroutine PhyImplSDHV3ChkConservation( xy_IndexCalcMethod, xy_SeaIceThickness, xyr_Press, xyz_Exner, xyr_Exner, xy_SurfTemp, xy_SurfHeatCapacity, xy_SoilHeatCap, xyr_HeatFlux, xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xy_SurfTempTransCoef, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux, xy_LatHeatFluxByMajCompIceSubl, xy_LatHeatFluxBySnowMelt, xy_LatHeatFluxBySeaIceMelt, xy_DeepSubSurfHeatFlux, xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DSoilMoistDt, xy_DSurfSnowDt, xy_DPsDt, xy_DSurfMajCompIceDt )
!
!
!
! A part of conservation of energy is checked.
!
! モジュール引用 ; USE statements
!
! 座標データ設定
! Axes data settings
!
use axesset, only: r_SSDepth, z_SSDepth ! subsurface grid at midpoint of layer
! 時刻管理
! Time control
!
use timeset, only: DelTime ! $ \Delta t $ [s]
! 物理定数設定
! Physical constants settings
!
use constants, only: Grav, CpDry, GasRDry
! $ R $ [J kg-1 K-1].
! 乾燥大気の気体定数.
! Gas constant of air
! 雪と海氷の定数の設定
! Setting constants of snow and sea ice
!
use constants_snowseaice, only: TempCondWater, SeaIceVolHeatCap , SeaIceThermCondCoef, TempBelowSeaIce
! 宣言文 ; Declaration statements
!
integer , intent(in):: xy_IndexCalcMethod (0:imax-1, 1:jmax)
!
! Index for calculation method
real(DP), intent(in):: xy_SeaIceThickness(0:imax-1, 1:jmax)
!
! Sea ice thickness
real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air pressure (half level)
real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
! Exner 関数 (整数レベル).
! Exner function (full level)
real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
! Exner 関数 (半整数レベル).
! Exner function (half level)
real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
! 地表面温度.
! Surface temperature
real(DP), intent(in):: xy_SurfHeatCapacity (0:imax-1, 1:jmax)
! 地表熱容量.
! Surface heat capacity
real(DP), intent(in):: xy_SoilHeatCap (0:imax-1, 1:jmax)
! 土壌熱容量 (J K-1 kg-1)
! Specific heat of soil (J K-1 kg-1)
real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
!!$ real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
!!$ ! 比湿フラックス.
!!$ ! Specific humidity flux
real(DP), intent(in):: xy_SurfH2OVapFlux(0:imax-1, 1:jmax)
! 惑星表面水蒸気フラックス.
! Water vapor flux at the surface
real(DP), intent(in):: xy_SurfLatentHeatFlux(0:imax-1, 1:jmax)
! 惑星表面潜熱フラックス.
! Latent heat flux at the surface
real(DP), intent(in):: xyr_SoilHeatFlux (0:imax-1, 1:jmax, 0:kslmax)
! 土壌中の熱フラックス (W m-2)
! Heat flux in sub-surface soil (W m-2)
real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in):: xyr_SoilTempTransCoef (0:imax-1, 1:jmax, 0:kslmax)
! 輸送係数：土壌温度.
! Transfer coefficient: soil temperature
real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
! 短波 (日射) フラックス.
! Shortwave (insolation) flux
real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
! 長波フラックス.
! Longwave flux
real(DP), intent(in):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
! 長波地表温度変化.
! Surface temperature tendency with longwave
real(DP), intent(in ):: xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax)
!
! Latent heat flux by major component ice sublimation
! (variable only for debug)
real(DP), intent(in ):: xy_LatHeatFluxBySnowMelt (0:imax-1, 1:jmax)
!
! Latent heat flux by melt (variable only for debug)
real(DP), intent(in ):: xy_LatHeatFluxBySeaIceMelt (0:imax-1, 1:jmax)
!
! Latent heat flux by sea ice melt (variable only for debug)
real(DP), intent(in):: xy_DeepSubSurfHeatFlux (0:imax-1, 1:jmax)
! 地中熱フラックス.
! "Deep subsurface heat flux"
! Heat flux at the bottom of surface/soil layer.
real(DP), intent(in ):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{T}{t} $ . 温度変化.
! Temperature tendency
real(DP), intent(in ):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ \DP{q}{t} $ . 質量混合比変化.
! Mass mixing ratio tendency
real(DP), intent(in ):: xy_DSurfTempDt (0:imax-1, 1:jmax)
! 地表面温度変化率 (K s-1)
! Surface temperature tendency (K s-1)
real(DP), intent(in ):: xyz_DSoilTempDt (0:imax-1, 1:jmax, 1:kslmax)
! $ \DP{Tg}{t} $ . 土壌温度変化 (K s-1)
! Temperature tendency (K s-1)
real(DP), intent(in ):: xy_DSoilMoistDt (0:imax-1, 1:jmax)
! 土壌温度時間変化率 (kg m-2 s-1)
! Soil temperature tendency (kg m-2 s-1)
real(DP), intent(in ):: xy_DSurfSnowDt (0:imax-1, 1:jmax)
! 積雪率時間変化率 (kg m-2 s-1)
! Surface snow amount tendency (kg m-2 s-1)
real(DP), intent(in ):: xy_DPsDt (0:imax-1, 1:jmax)
real(DP), intent(in ):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)
! 作業変数
! Work variables
!
real(DP) :: xy_SurfRadSFlux (0:imax-1, 1:jmax)
real(DP) :: xy_SurfRadLFlux (0:imax-1, 1:jmax)
real(DP) :: xy_SurfSensHeatFlux(0:imax-1, 1:jmax)
real(DP) :: xy_SurfSoilHeatCondFlux(0:imax-1, 1:jmax)
real(DP) :: xy_SeaIceHeatCondFlux (0:imax-1, 1:jmax)
real(DP) :: xy_Residual (0:imax-1, 1:jmax)
real(DP) :: xy_SumAtmRate (0:imax-1, 1:jmax)
real(DP) :: MaxResidual
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. phy_implicit_sdh_V3_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
xy_SurfRadSFlux = xyr_RadSFlux(:,:,0)
xy_SurfRadLFlux = xyr_RadLFlux(:,:,0) + xyra_DelRadLFlux(:,:,0,0) * xy_DSurfTempDt * ( 2.0_DP * DelTime ) + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0_DP * DelTime )
xy_SurfSoilHeatCondFlux = xyr_SoilHeatFlux(:,:,0) - xyr_SoilTempTransCoef(:,:,0) * ( xyz_DSoilTempDt(:,:,1) - xy_DSurfTempDt ) * ( 2.0_DP * DelTime )
xy_SurfSensHeatFlux = xyr_HeatFlux(:,:,0) - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef * ( xyz_DTempDt(:,:,1) / xyz_Exner(:,:,1) - xy_DSurfTempDt / xyr_Exner(:,:,0) ) * ( 2.0_DP * DelTime )
xy_SeaIceHeatCondFlux = - SeaIceThermCondCoef * ( xy_SurfTemp + xy_DSurfTempDt * ( 2.0_DP * DelTime ) - TempBelowSeaIce ) / xy_SeaIceThickness
!-----
! Atmospheric heating
!
xy_SumAtmRate = 0.0_DP
do k = kmax, 1, -1
xy_SumAtmRate = xy_SumAtmRate + CpDry * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav * xyz_DTempDt(:,:,k)
end do
!
xy_Residual = - ( 0.0_DP - xy_SurfSensHeatFlux ) - xy_SumAtmRate
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
!!$ select case ( xy_IndexCalcMethod(i,j) )
!!$ case ( IndexSeaIce )
!!$ case default
!!$ MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
!!$ end select
end do
end do
if ( MaxResidual > 1.0d-10 ) then
call MessageNotify( 'M', module_name, 'Atm. sensible heating res. : %f.', d = (/ MaxResidual /) )
end if
!-----
! Land surface
!
xy_SumAtmRate = xy_SurfHeatCapacity * xy_DSurfTempDt
!
xy_Residual = xy_SurfRadSFlux + xy_SurfRadLFlux + xy_SurfSensHeatFlux + xy_SurfLatentHeatFlux - xy_SurfSoilHeatCondFlux + xy_LatHeatFluxByMajCompIceSubl + xy_LatHeatFluxBySnowMelt + xy_SumAtmRate
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
case ( IndexSeaIce )
! sea ice
case ( IndexSlabOcean )
! slab ocean
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
if ( MaxResidual > 1.0e-10_DP ) then
call MessageNotify( 'M', module_name, 'Land surf. heat budget res.: %f.', d = (/ MaxResidual /) )
end if
!-----
! Soil heating
!
xy_SumAtmRate = 0.0_DP
do k = 1, kslmax
xy_SumAtmRate = xy_SumAtmRate + xy_SoilHeatCap * ( r_SSDepth(k-1) - r_SSDepth(k) ) * xyz_DSoilTempDt(:,:,k)
end do
!
xy_Residual = - ( xy_SurfSoilHeatCondFlux - xy_DeepSubSurfHeatFlux ) - xy_SumAtmRate
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
case ( IndexSeaIce )
! sea ice
case ( IndexSlabOcean )
! slab ocean
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
if ( MaxResidual > 1.0e-10_DP ) then
call MessageNotify( 'M', module_name, 'Soil heating res. : %f.', d = (/ MaxResidual /) )
end if
!-----
! Slab ocean heating
!
xy_SumAtmRate = SOHeatCapacity * xy_DSurfTempDt
!
xy_Residual = - ( xy_SurfRadSFlux + xy_SurfRadLFlux + xy_SurfSensHeatFlux + xy_SurfLatentHeatFlux ) - xy_SumAtmRate
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
case ( IndexSeaIce )
! sea ice
case ( IndexSlabOcean )
! slab ocean
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
if ( MaxResidual > 1.0e-10_DP ) then
call MessageNotify( 'M', module_name, 'Slab ocean heating res. : %f.', d = (/ MaxResidual /) )
end if
!-----
! Sea ice heating
!
xy_SumAtmRate = SeaIceVolHeatCap * xy_SeaIceThickness * xy_DSurfTempDt
!
xy_Residual = - ( ( xy_SurfRadSFlux + xy_SurfRadLFlux + xy_SurfSensHeatFlux + xy_SurfLatentHeatFlux ) - xy_SeaIceHeatCondFlux ) - xy_LatHeatFluxBySeaIceMelt - xy_SumAtmRate
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
case ( IndexSeaIce )
! sea ice
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
case ( IndexSlabOcean )
! slab ocean
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
if ( MaxResidual > 1.0e-10_DP ) then
call MessageNotify( 'M', module_name, 'Sea ice heating res. : %f.', d = (/ MaxResidual /) )
end if
!-----
! Atmospheric moistening
!
xy_SumAtmRate = 0.0_DP
do k = kmax, 1, -1
xy_SumAtmRate = xy_SumAtmRate + ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav * xyzf_DQMixDt(:,:,k,IndexH2OVap)
end do
!
xy_Residual = - ( 0.0_DP - xy_SurfH2OVapFlux ) - xy_SumAtmRate
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
end do
end do
if ( MaxResidual > 1.0e-10_DP ) then
call MessageNotify( 'M', module_name, 'Atm. moistening res. : %f.', d = (/ MaxResidual /) )
end if
!-----
! Land water budget
!
xy_SumAtmRate = 0.0_DP
do k = kmax, 1, -1
xy_SumAtmRate = xy_SumAtmRate + ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav * xyzf_DQMixDt(:,:,k,IndexH2OVap)
end do
!
xy_Residual = xy_DSoilMoistDt + xy_DSurfSnowDt + xy_SumAtmRate
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
case ( IndexSeaIce )
! sea ice
case ( IndexSlabOcean )
! slab ocean
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
if ( MaxResidual > 1.0e-10_DP ) then
call MessageNotify( 'M', module_name, 'Land water budget res. : %f.', d = (/ MaxResidual /) )
end if
!-----
! Atmospheric mass budget
!
xy_Residual = xy_DPsDt / Grav + xy_DSurfMajCompIceDt
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
end do
end do
if ( MaxResidual > 1.0e-10_DP ) then
call MessageNotify( 'M', module_name, 'Atm. mass budget res. : %f.', d = (/ MaxResidual /) )
end if
end subroutine PhyImplSDHV3ChkConservation
| Subroutine : | |||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , intent(in)
| ||
| xy_SeaIceThickness(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SnowFrac(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_Press(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyz_Exner(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_Exner(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xy_SurfSoilHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfTempTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfQVapTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfHumidCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfHeatCapacity(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_RadSFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_RadLFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(in)
| ||
| xy_DeepSubSurfHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in )
| ||
| xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(in )
| ||
| xy_DSurfTempDt(0:imax-1, 1:jmax) : | real(DP), intent(in )
|
A part of conservation of energy is checked.
subroutine PhyImplSDHV3ChkConservationTQ( xy_IndexCalcMethod, xy_SeaIceThickness, xy_SnowFrac, xyr_Press, xyz_Exner, xyr_Exner, xy_SurfTemp, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfSoilHeatFlux, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xy_SurfHumidCoef, xy_SurfHeatCapacity, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux, xy_DeepSubSurfHeatFlux, xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt )
!
!
!
! A part of conservation of energy is checked.
!
! モジュール引用 ; USE statements
!
! 座標データ設定
! Axes data settings
!
use axesset, only: r_SSDepth, z_SSDepth ! subsurface grid at midpoint of layer
! 時刻管理
! Time control
!
use timeset, only: DelTime ! $ \Delta t $ [s]
! 物理定数設定
! Physical constants settings
!
use constants, only: Grav, CpDry, GasRDry, LatentHeat
! 雪と海氷の定数の設定
! Setting constants of snow and sea ice
!
use constants_snowseaice, only: TempCondWater, SeaIceVolHeatCap , SeaIceThermCondCoef, TempBelowSeaIce
! 飽和比湿の算出
! Evaluation of saturation specific humidity
!
use saturate, only: xy_CalcQVapSatOnLiq, xy_CalcQVapSatOnSol, xy_CalcDQVapSatDTempOnLiq, xy_CalcDQVapSatDTempOnSol
! 宣言文 ; Declaration statements
!
integer , intent(in):: xy_IndexCalcMethod (0:imax-1, 1:jmax)
!
! Index for calculation method
real(DP), intent(in):: xy_SeaIceThickness(0:imax-1, 1:jmax)
!
! Sea ice thickness
real(DP), intent(in):: xy_SnowFrac(0:imax-1, 1:jmax)
!
! Snow fraction
real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air pressure (half level)
real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
! Exner 関数 (整数レベル).
! Exner function (full level)
real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
! Exner 関数 (半整数レベル).
! Exner function (half level)
real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
! 地表面温度.
! Surface temperature
real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
! 比湿フラックス.
! Specific humidity flux
real(DP), intent(in):: xy_SurfSoilHeatFlux(0:imax-1, 1:jmax)
! 惑星表面土壌熱伝導フラックス.
! Soil heat conduction flux at the surface
real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in):: xy_SurfQVapTransCoef (0:imax-1, 1:jmax)
! 輸送係数：比湿.
! Transfer coefficient: specific humidity
real(DP), intent(in):: xy_SurfHumidCoef (0:imax-1, 1:jmax)
! 地表湿潤度.
! Surface humidity coefficient
real(DP), intent(in):: xy_SurfHeatCapacity (0:imax-1, 1:jmax)
! 地表熱容量.
! Surface heat capacity
real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
! 短波 (日射) フラックス.
! Shortwave (insolation) flux
real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
! 長波フラックス.
! Longwave flux
real(DP), intent(in):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
! 長波地表温度変化.
! Surface temperature tendency with longwave
real(DP), intent(in):: xy_DeepSubSurfHeatFlux (0:imax-1, 1:jmax)
! 地中熱フラックス.
! "Deep subsurface heat flux"
! Heat flux at the bottom of surface/soil layer.
real(DP), intent(in ):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{T}{t} $ . 温度変化.
! Temperature tendency
real(DP), intent(in ):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ \DP{q}{t} $ . 質量混合比変化.
! Mass mixing ratio tendency
real(DP), intent(in ):: xy_DSurfTempDt (0:imax-1, 1:jmax)
! 地表面温度変化率 (K s-1)
! Surface temperature tendency (K s-1)
! 作業変数
! Work variables
!
real(DP):: xy_SurfQVapSatOnLiq(0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfQVapSatOnSol (0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfQVapSat (0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfDQVapSatDTempOnLiq (0:imax-1, 1:jmax)
! 地表飽和比湿変化.
! Saturated specific humidity tendency on surface
real(DP):: xy_SurfDQVapSatDTempOnSol (0:imax-1, 1:jmax)
! 地表飽和比湿変化.
! Saturated specific humidity tendency on surface
real(DP):: xy_SurfDQVapSatDTemp (0:imax-1, 1:jmax)
! 地表飽和比湿変化.
! Saturated specific humidity tendency on surface
real(DP) :: xy_SurfRadSFlux (0:imax-1, 1:jmax)
real(DP) :: xy_SurfRadLFlux (0:imax-1, 1:jmax)
real(DP) :: xy_SurfSensHeatFlux (0:imax-1, 1:jmax)
real(DP) :: xy_SurfH2OVapFlux (0:imax-1, 1:jmax)
real(DP) :: xy_SurfLatentHeatFlux (0:imax-1, 1:jmax)
real(DP) :: xy_SurfSoilHeatCondFlux(0:imax-1, 1:jmax)
real(DP) :: xy_SeaIceHeatCondFlux (0:imax-1, 1:jmax)
real(DP) :: xy_Residual (0:imax-1, 1:jmax)
real(DP) :: xy_SumAtmRate (0:imax-1, 1:jmax)
real(DP) :: MaxResidual
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. phy_implicit_sdh_V3_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 飽和比湿の計算
! Calculate saturated specific humidity
!
xy_SurfQVapSatOnLiq = xy_CalcQVapSatOnLiq ( xy_SurfTemp, xyr_Press(:,:,0) )
xy_SurfQVapSatOnSol = xy_CalcQVapSatOnSol ( xy_SurfTemp, xyr_Press(:,:,0) )
xy_SurfQVapSat = ( 1.0_DP - xy_SnowFrac ) * xy_SurfQVapSatOnLiq + xy_SnowFrac * xy_SurfQVapSatOnSol
xy_SurfDQVapSatDTempOnLiq = xy_CalcDQVapSatDTempOnLiq( xy_SurfTemp, xy_SurfQVapSatOnLiq )
xy_SurfDQVapSatDTempOnSol = xy_CalcDQVapSatDTempOnSol( xy_SurfTemp, xy_SurfQVapSatOnSol )
xy_SurfDQVapSatDTemp = ( 1.0_DP - xy_SnowFrac ) * xy_SurfDQVapSatDTempOnLiq + xy_SnowFrac * xy_SurfDQVapSatDTempOnSol
xy_SurfRadSFlux = xyr_RadSFlux(:,:,0)
xy_SurfRadLFlux = xyr_RadLFlux(:,:,0) + xyra_DelRadLFlux(:,:,0,0) * xy_DSurfTempDt * ( 2.0_DP * DelTime ) + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0_DP * DelTime )
xy_SurfSoilHeatCondFlux = xy_SurfSoilHeatFlux
xy_SurfSensHeatFlux = xyr_HeatFlux(:,:,0) - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef * ( xyz_DTempDt(:,:,1) / xyz_Exner(:,:,1) - xy_DSurfTempDt / xyr_Exner(:,:,0) ) * ( 2.0_DP * DelTime )
xy_SeaIceHeatCondFlux = - SeaIceThermCondCoef * ( xy_SurfTemp + xy_DSurfTempDt * ( 2.0_DP * DelTime ) - TempBelowSeaIce ) / xy_SeaIceThickness
n = IndexH2OVap
xy_SurfH2OVapFlux = xyrf_QMixFlux(:,:,0,n) - xy_SurfHumidCoef * xy_SurfQVapTransCoef * ( xyzf_DQMixDt(:,:,1,n) - xy_SurfDQVapSatDTemp * xy_DSurfTempDt ) * 2.0_DP * DelTime
xy_SurfLatentHeatFlux = LatentHeat * xy_SurfH2OVapFlux
!-----
! Atmospheric heating
!
xy_SumAtmRate = 0.0_DP
do k = kmax, 1, -1
xy_SumAtmRate = xy_SumAtmRate + CpDry * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav * xyz_DTempDt(:,:,k)
end do
!
xy_Residual = - ( 0.0_DP - xy_SurfSensHeatFlux ) - xy_SumAtmRate
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
!!$ select case ( xy_IndexCalcMethod(i,j) )
!!$ case ( IndexSeaIce )
!!$ case default
!!$ MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
!!$ end select
end do
end do
if ( MaxResidual > 1.0e-10_DP ) then
call MessageNotify( 'M', module_name, 'Atm. sensible heating res. : %f.', d = (/ MaxResidual /) )
end if
!-----
! Land surface
!
xy_SumAtmRate = xy_SurfHeatCapacity * xy_DSurfTempDt
!
xy_Residual = xy_SurfRadSFlux + xy_SurfRadLFlux + xy_SurfSensHeatFlux + xy_SurfLatentHeatFlux - xy_SurfSoilHeatCondFlux + xy_SumAtmRate
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
case ( IndexSeaIce )
! sea ice
case ( IndexSlabOcean )
! slab ocean
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
if ( MaxResidual > 1.0e-10_DP ) then
call MessageNotify( 'M', module_name, 'Land surf. heat budget res.: %f.', d = (/ MaxResidual /) )
end if
!!$ !-----
!!$ ! Soil heating
!!$ !
!!$ xy_SumAtmRate = 0.0_DP
!!$ do k = 1, kslmax
!!$ xy_SumAtmRate = xy_SumAtmRate &
!!$ & + xy_SoilHeatCap * ( r_SSDepth(k-1) - r_SSDepth(k) ) &
!!$ & * xyz_DSoilTempDt(:,:,k)
!!$ end do
!!$ !
!!$ xy_Residual = &
!!$ & - ( xy_SurfSoilHeatCondFlux - xy_DeepSubSurfHeatFlux ) &
!!$ & - xy_SumAtmRate
!!$ !
!!$ MaxResidual = 0.0_DP
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ select case ( xy_IndexCalcMethod(i,j) )
!!$ case ( IndexLand )
!!$ ! land
!!$ MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
!!$ case ( IndexSeaIce )
!!$ ! sea ice
!!$ case ( IndexSlabOcean )
!!$ ! slab ocean
!!$ case ( IndexOceanPresSST )
!!$ ! open ocean
!!$ case default
!!$ call MessageNotify( 'E', module_name, 'Unexpected Error.' )
!!$ end select
!!$ end do
!!$ end do
!!$ if ( MaxResidual > 1.0d-10 ) then
!!$ call MessageNotify( 'M', module_name, &
!!$ & 'Soil heating res. : %f.', d = (/ MaxResidual /) )
!!$ end if
!-----
! Slab ocean heating
!
xy_SumAtmRate = SOHeatCapacity * xy_DSurfTempDt
!
xy_Residual = - ( xy_SurfRadSFlux + xy_SurfRadLFlux + xy_SurfSensHeatFlux + xy_SurfLatentHeatFlux ) - xy_SumAtmRate
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
case ( IndexSeaIce )
! sea ice
case ( IndexSlabOcean )
! slab ocean
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
if ( MaxResidual > 1.0e-10_DP ) then
call MessageNotify( 'M', module_name, 'Slab ocean heating res. : %f.', d = (/ MaxResidual /) )
end if
!-----
! Sea ice heating
!
xy_SumAtmRate = SeaIceVolHeatCap * xy_SeaIceThickness * xy_DSurfTempDt
!
xy_Residual = - ( ( xy_SurfRadSFlux + xy_SurfRadLFlux + xy_SurfSensHeatFlux + xy_SurfLatentHeatFlux ) - xy_SeaIceHeatCondFlux ) - xy_SumAtmRate
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
case ( IndexSeaIce )
! sea ice
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
case ( IndexSlabOcean )
! slab ocean
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
if ( MaxResidual > 1.0e-10_DP ) then
call MessageNotify( 'M', module_name, 'Sea ice heating res. : %f.', d = (/ MaxResidual /) )
end if
!-----
! Atmospheric moistening
!
xy_SumAtmRate = 0.0_DP
do k = kmax, 1, -1
xy_SumAtmRate = xy_SumAtmRate + ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav * xyzf_DQMixDt(:,:,k,IndexH2OVap)
end do
!
xy_Residual = - ( 0.0_DP - xy_SurfH2OVapFlux ) - xy_SumAtmRate
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
end do
end do
if ( MaxResidual > 1.0e-10_DP ) then
call MessageNotify( 'M', module_name, 'Atm. moistening res. : %f.', d = (/ MaxResidual /) )
end if
!-----
! Land water budget
!
xy_SumAtmRate = 0.0_DP
do k = kmax, 1, -1
xy_SumAtmRate = xy_SumAtmRate + ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav * xyzf_DQMixDt(:,:,k,IndexH2OVap)
end do
!
!!$ xy_Residual = xy_DSoilMoistDt + xy_DSurfSnowDt + xy_SumAtmRate
xy_Residual = xy_SumAtmRate - xy_SurfH2OVapFlux
!
MaxResidual = 0.0_DP
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
case ( IndexSeaIce )
! sea ice
case ( IndexSlabOcean )
! slab ocean
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
if ( MaxResidual > 1.0e-10_DP ) then
call MessageNotify( 'M', module_name, 'Land water budget res. : %f.', d = (/ MaxResidual /) )
end if
!!$ !-----
!!$ ! Atmospheric mass budget
!!$ !
!!$ xy_Residual = xy_DPsDt / Grav + xy_DSurfMajCompIceDt
!!$ !
!!$ MaxResidual = 0.0_DP
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ MaxResidual = max( MaxResidual, abs( xy_Residual(i,j) ) )
!!$ end do
!!$ end do
!!$ if ( MaxResidual > 1.0d-10 ) then
!!$ call MessageNotify( 'M', module_name, &
!!$ & 'Atm. mass budget res. : %f.', d = (/ MaxResidual /) )
!!$ end if
end subroutine PhyImplSDHV3ChkConservationTQ
| Subroutine : | |||
| IndexSpc : | integer , intent(in) | ||
| xy_Ps(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_SurfLatentHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_SoilHeatFlux(0:imax-1, 1:jmax, 0:kslmax) : | real(DP), intent(in)
| ||
| xyr_SoilTempTransCoef(0:imax-1, 1:jmax, 0:kslmax) : | real(DP), intent(in)
| ||
| xyr_RadSFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_RadLFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_DeepSubSurfHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_SoilTemp(0:imax-1, 1:jmax, 1:kslmax) : | real(DP), intent(in)
| ||
| xy_SurfLiqB(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfSolB(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfHeatCapacity(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SoilHeatCap(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xy_SoilHeatDiffCoef(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , intent(in )
| ||
| xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(in)
| ||
| xyz_Exner(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_Exner(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_SurfTempTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_LatHeatFluxByOtherSpc(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xyza_ArgTempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1) : | real(DP), intent(in )
| ||
| xyz_ArgTempVec(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in )
| ||
| xyaa_ArgSurfMtx(0:imax-1, 1:jmax, 0:0, -1:1) : | real(DP), intent(in )
| ||
| xy_ArgSurfRH(0:imax-1,1:jmax) : | real(DP), intent(in ) | ||
| xyaa_ArgSoilTempMtx(0:imax-1, 1:jmax, 1:kslmax,-1:1) : | real(DP), intent(in )
| ||
| xya_ArgSoilTempVec(0:imax-1, 1:jmax, 1:kslmax) : | real(DP), intent(in )
| ||
| xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||
| xy_DSurfTempDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xyz_DSoilTempDt(0:imax-1, 1:jmax, 1:kslmax) : | real(DP), intent(inout)
| ||
| xy_DSurfLiqDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xy_DSurfSolDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xy_LatHeatFluxBySnowMelt(0:imax-1, 1:jmax) : | real(DP), intent(out )
|
融雪による時間変化率の修正を行います.
Correction of tendencies due to melt of snow.
subroutine PhyImplSDHV3IceSnowPhaseChgCor( IndexSpc, xy_Ps, xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfLiqB, xy_SurfSolB, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xy_IndexCalcMethod, xyra_DelRadLFlux, xyz_Exner, xyr_Exner, xy_SurfTempTransCoef, xy_LatHeatFluxByOtherSpc, xyza_ArgTempMtx, xyz_ArgTempVec, xyaa_ArgSurfMtx, xy_ArgSurfRH, xyaa_ArgSoilTempMtx, xya_ArgSoilTempVec, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DSurfLiqDt, xy_DSurfSolDt, xy_LatHeatFluxBySnowMelt )
!
! 融雪による時間変化率の修正を行います.
!
! Correction of tendencies due to melt of snow.
!
! モジュール引用 ; USE statements
!
! 座標データ設定
! Axes data settings
!
use axesset, only: r_SSDepth, z_SSDepth ! subsurface grid at midpoint of layer
! 物理定数設定
! Physical constants settings
!
use constants, only: CpDry, LatentHeatFusion
! $ L $ [J kg-1] .
! 融解の潜熱.
! Latent heat of fusion
! 雪と海氷の定数の設定
! Setting constants of snow and sea ice
!
use constants_snowseaice, only: TempCondWater
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 陰解法による時間積分のためのルーチン
! Routines for time integration with implicit scheme
!
use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3
! 主成分相変化
! Phase change of atmospheric major component
!
use saturate_major_comp, only : SaturateMajorCompCalcCondTemp, SaturateMajorCompInqLatentHeat
! 宣言文 ; Declaration statements
!
integer , intent(in):: IndexSpc
real(DP), intent(in):: xy_Ps(0:imax-1, 1:jmax)
!
! Surface pressure
real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(in):: xy_SurfLatentHeatFlux(0:imax-1, 1:jmax)
! 惑星表面潜熱フラックス.
! Latent heat flux at the surface
real(DP), intent(in):: xyr_SoilHeatFlux (0:imax-1, 1:jmax, 0:kslmax)
! 土壌中の熱フラックス (W m-2)
! Heat flux in sub-surface soil (W m-2)
real(DP), intent(in):: xyr_SoilTempTransCoef (0:imax-1, 1:jmax, 0:kslmax)
! 輸送係数：土壌温度.
! Transfer coefficient: soil temperature
real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
! 短波 (日射) フラックス.
! Shortwave (insolation) flux
real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
! 長波フラックス.
! Longwave flux
real(DP), intent(in):: xy_DeepSubSurfHeatFlux (0:imax-1, 1:jmax)
! 地中熱フラックス.
! "Deep subsurface heat flux"
! Heat flux at the bottom of surface/soil layer.
real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
! 地表面温度.
! Surface temperature
real(DP), intent(in):: xyz_SoilTemp (0:imax-1, 1:jmax, 1:kslmax)
! 土壌温度 (K)
! Soil temperature (K)
real(DP), intent(in):: xy_SurfLiqB (0:imax-1, 1:jmax)
!
! Surface liquid amount
real(DP), intent(in):: xy_SurfSolB (0:imax-1, 1:jmax)
! 積雪量.
! Surface snow amount.
real(DP), intent(in):: xy_SurfHeatCapacity (0:imax-1, 1:jmax)
! 地表熱容量.
! Surface heat capacity
real(DP), intent(in ):: xy_SoilHeatCap (0:imax-1, 1:jmax)
! 土壌熱容量 (J K-1 kg-1)
! Specific heat of soil (J K-1 kg-1)
real(DP), intent(in ):: xy_SoilHeatDiffCoef (0:imax-1, 1:jmax)
! 土壌熱伝導係数 (J m-3 K-1)
! Heat conduction coefficient of soil (J m-3 K-1)
integer , intent(in ) :: xy_IndexCalcMethod(0:imax-1, 1:jmax)
!
! Index for calculation method
real(DP), intent(in):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
! 長波地表温度変化.
! Surface temperature tendency with longwave
real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
! Exner 関数 (整数レベル).
! Exner function (full level)
real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
! Exner 関数 (半整数レベル).
! Exner function (half level)
real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in ):: xyza_ArgTempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
! 温度陰解行列.
! Implicit matrix about temperature
real(DP), intent(in ):: xyz_ArgTempVec(0:imax-1, 1:jmax, 1:kmax)
! 温度陰解ベクトル.
! Implicit vector about temperature
real(DP), intent(in ):: xyaa_ArgSurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
! 惑星表面エネルギー収支用陰解行列
! Implicit matrix for surface energy balance
real(DP), intent(in ):: xy_ArgSurfRH(0:imax-1,1:jmax)
real(DP), intent(in ):: xyaa_ArgSoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
! 土壌温度拡散方程式の行列
! Matrix for diffusion equation of soil temperature
real(DP), intent(in ):: xya_ArgSoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
! 土壌温度拡散方程式のベクトル
! Vector for diffusion equation of soil temperature
real(DP), intent(in ):: xy_LatHeatFluxByOtherSpc(0:imax-1, 1:jmax)
!
! Latent heat flux by other specie
real(DP), intent(inout):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{T}{t} $ . 温度変化.
! Temperature tendency
real(DP), intent(inout):: xy_DSurfTempDt (0:imax-1, 1:jmax)
! 地表面温度変化率 (K s-1)
! Surface temperature tendency (K s-1)
real(DP), intent(inout):: xyz_DSoilTempDt (0:imax-1, 1:jmax, 1:kslmax)
! $ \DP{Tg}{t} $ . 土壌温度変化 (K s-1)
! Temperature tendency (K s-1)
real(DP), intent(inout):: xy_DSurfLiqDt (0:imax-1, 1:jmax)
! 土壌温度時間変化率 (kg m-2 s-1)
! Soil temperature tendency (kg m-2 s-1)
real(DP), intent(inout):: xy_DSurfSolDt (0:imax-1, 1:jmax)
! 積雪率時間変化率 (kg m-2 s-1)
! Surface snow amount tendency (kg m-2 s-1)
real(DP), intent(out ):: xy_LatHeatFluxBySnowMelt(0:imax-1, 1:jmax)
!
! Latent heat flux by melt (variable only for debug)
! 作業変数
! Work variables
!
real(DP):: xy_DSurfLiqDtSave(0:imax-1, 1:jmax)
real(DP):: xy_DSurfSolDtSave (0:imax-1, 1:jmax)
real(DP):: xy_TempCond(0:imax-1, 1:jmax)
logical :: xy_FlagCalc(0:imax-1, 1:jmax)
integer :: xy_IndexMeltOrFreeze(0:imax-1, 1:jmax)
integer, parameter :: IndexOthers = 0
integer, parameter :: IndexMelt = 1
integer, parameter :: IndexFreeze = 2
real(DP):: xyza_TempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
! 温度陰解行列.
! Implicit matrix about temperature
real(DP):: xyz_TempVec(0:imax-1, 1:jmax, 1:kmax)
! 温度陰解ベクトル.
! Implicit vector about temperature
real(DP):: xyaa_SurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
! 惑星表面エネルギー収支用陰解行列
! Implicit matrix for surface energy balance
real(DP):: xy_SurfRH(0:imax-1,1:jmax)
real(DP):: xyaa_SoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
! 土壌温度拡散方程式の行列
! Matrix for diffusion equation of soil temperature
real(DP):: xya_SoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
! 土壌温度拡散方程式のベクトル
! Vector for diffusion equation of soil temperature
real(DP):: xyaa_TempSoilTempLUMtx (0:imax-1, 1:jmax, -kslmax:kmax, -1:1)
! LU 行列.
! LU matrix
real(DP):: xya_DelTempSoilTempLUVec (0:imax-1, 1:jmax, -kslmax:kmax)
! $ T, Tg $ の時間変化.
! Tendency of $ T $ and $ Tg |
real(DP):: LatentHeatLocal
real(DP):: LatentHeatFluxByMelt
real(DP):: SenHeatFluxA
real(DP):: LatHeatFluxA
real(DP):: CondHeatFluxA
real(DP):: ValueAlpha
real(DP):: SurfTempATentative
real(DP):: SoilTempATentative
real(DP):: SurfLiqATentative
real(DP):: xy_SurfLiqATentativeSave(0:imax-1, 1:jmax)
real(DP):: SurfSolATentative
real(DP):: xy_SurfSolATentativeSave(0:imax-1, 1:jmax)
real(DP):: DelSurfSol
real(DP) :: xy_TempMajCompCond(0:imax-1, 1:jmax)
real(DP) :: SurfMajCompIceATentative
real(DP) :: xy_SurfRadSFlux (0:imax-1, 1:jmax)
real(DP) :: xy_SurfRadLFlux (0:imax-1, 1:jmax)
real(DP) :: xy_SurfSoilHeatCondFlux(0:imax-1, 1:jmax)
real(DP) :: xy_SurfSensHeatFlux (0:imax-1, 1:jmax)
real(DP) :: xy_SeaIceHeatCondFlux (0:imax-1, 1:jmax)
real(DP) :: xy_HeatingTendency (0:imax-1, 1:jmax)
real(DP) :: LatHeatFluxBySnowMelt
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: l
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. phy_implicit_sdh_V3_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
!!$ ! 計算時間計測開始
!!$ ! Start measurement of computation time
!!$ !
!!$ call TimesetClockStart( module_name )
!
! check flag of snow melt
!
if ( ( .not. FlagBucketModel ) .or. ( .not. FlagSnow ) ) then
xy_LatHeatFluxBySnowMelt = 0.0_DP
return
end if
if ( kslmax == 0 ) then
call MessageNotify( 'E', module_name, 'kslmax <= 0 in PhyImplSDHV3IceSnowPhaseChgCor.' )
else
xy_DSurfLiqDtSave = xy_DSurfLiqDt
xy_DSurfSolDtSave = xy_DSurfSolDt
select case ( IndexSpc )
case ( IndexSpcMajComp )
call SaturateMajorCompCalcCondTemp( xy_Ps, xy_TempCond )
LatentHeatLocal = SaturateMajorCompInqLatentHeat()
case ( IndexSpcH2O )
xy_TempCond = TempCondWater
LatentHeatLocal = LatentHeatFusion
case default
call MessageNotify( 'E', module_name, 'Undefined IndexSpc, %d.', i = (/ IndexSpc /) )
end select
xy_SurfLiqATentativeSave = xy_SurfLiqB + xy_DSurfLiqDt * ( 2.0_DP * DelTime )
xy_SurfSolATentativeSave = xy_SurfSolB + xy_DSurfSolDt * ( 2.0_DP * DelTime )
!----------
! A case that a part of snow/ice melt or soil moisture freeze
!----------
! Melt
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
SurfTempATentative = xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
SurfSolATentative = xy_SurfSolATentativeSave(i,j)
if ( ( SurfSolATentative > 0.0_DP ) .and. ( SurfTempATentative > xy_TempCond(i,j) ) ) then
xy_FlagCalc (i,j) = .true.
else
xy_FlagCalc (i,j) = .false.
end if
case default
xy_FlagCalc (i,j) = .false.
end select
end do
end do
! Freeze
select case ( IndexSpc )
case ( IndexSpcMajComp )
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
SurfTempATentative = xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
if ( SurfTempATentative < xy_TempCond(i,j) ) then
xy_FlagCalc (i,j) = .true.
end if
end select
end do
end do
!!$ case ( IndexSpcH2O )
!!$
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ select case ( xy_IndexCalcMethod(i,j) )
!!$ case ( IndexLand )
!!$ SurfTempATentative = xy_SurfTemp(i,j) &
!!$ & + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
!!$ SurfLiqATentative = xy_SurfLiqATentativeSave(i,j)
!!$ if ( &
!!$ & ( SurfLiqATentative > 0.0_DP ) .and. &
!!$ & ( SurfTempATentative < xy_TempCond(i,j) ) &
!!$ & ) then
!!$ xy_FlagCalc (i,j) = .true.
!!$ end if
!!$ end select
!!$
!!$ end do
!!$ end do
!!$
end select
xyza_TempMtx = xyza_ArgTempMtx
xyz_TempVec = xyz_ArgTempVec
!
xyaa_SurfMtx = xyaa_ArgSurfMtx
xy_SurfRH = xy_ArgSurfRH
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
xyaa_SurfMtx(i,j,0, 0) = 1.0_DP
xyaa_SurfMtx(i,j,0, 1) = 0.0_DP
xy_SurfRH (i,j) = xy_TempCond(i,j) - xy_SurfTemp(i,j)
end if
end do
end do
!
xyaa_SoilTempMtx = xyaa_ArgSoilTempMtx
xya_SoilTempVec = xya_ArgSoilTempVec
! 温度の計算
! Calculate temperature and specific humidity
!
do l = -1, 1
do k = 1, kslmax
xyaa_TempSoilTempLUMtx(:,:,-k,-l) = xyaa_SoilTempMtx(:,:,k,l)
end do
k = 0
xyaa_TempSoilTempLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
do k = 1, kmax
xyaa_TempSoilTempLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
end do
end do
call PhyImplLUDecomp3( xyaa_TempSoilTempLUMtx, imax * jmax, kmax + 1 + kslmax )
do k = 1, kslmax
xya_DelTempSoilTempLUVec(:,:,-k) = xya_SoilTempVec(:,:,k)
end do
k = 0
xya_DelTempSoilTempLUVec(:,:,k) = xy_SurfRH
do k = 1, kmax
xya_DelTempSoilTempLUVec(:,:,k) = xyz_TempVec(:,:,k)
end do
call PhyImplLUSolve3( xya_DelTempSoilTempLUVec, xyaa_TempSoilTempLUMtx, 1, imax * jmax , kmax + 1 + kslmax )
do k = 1, kslmax
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
xyz_DSoilTempDt(i,j,k) = xya_DelTempSoilTempLUVec(i,j,-k) / ( 2.0_DP * DelTime )
case default
xyz_DSoilTempDt(i,j,k) = 0.0_DP
end select
end if
end do
end do
end do
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$ case ( IndexSeaIce )
!!$ ! sea ice
!!$ xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$ case ( IndexSlabOcean )
!!$ ! slab ocean
!!$ xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$ case ( IndexOceanPresSST )
!!$ ! open ocean
!!$ xy_DSurfTempDt(i,j) = 0.0_DP
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end if
end do
end do
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
xyz_DTempDt(i,j,k) = xya_DelTempSoilTempLUVec(i,j,k) / ( 2.0_DP * DelTime )
end if
end do
end do
end do
!----------
! Surface fluxes used below
!----------
xy_SurfRadSFlux = xyr_RadSFlux(:,:,0)
xy_SurfRadLFlux = xyr_RadLFlux(:,:,0) + xyra_DelRadLFlux(:,:,0,0) * xy_DSurfTempDt * ( 2.0_DP * DelTime ) + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0_DP * DelTime )
xy_SurfSoilHeatCondFlux = xyr_SoilHeatFlux(:,:,0) - xyr_SoilTempTransCoef(:,:,0) * ( xyz_DSoilTempDt(:,:,1) - xy_DSurfTempDt ) * ( 2.0d0 * DelTime )
xy_SurfSensHeatFlux = xyr_HeatFlux(:,:,0) - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef * ( xyz_DTempDt(:,:,1) / xyz_Exner(:,:,1) - xy_DSurfTempDt / xyr_Exner(:,:,0) ) * ( 2.0_DP * DelTime )
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
xy_LatHeatFluxBySnowMelt(i,j) = - xy_SurfRadSFlux(i,j) - xy_SurfRadLFlux(i,j) - xy_SurfSensHeatFlux(i,j) - xy_SurfLatentHeatFlux(i,j) + xy_SurfSoilHeatCondFlux(i,j) - xy_LatHeatFluxByOtherSpc(i,j) - xy_SurfHeatCapacity(i,j) * xy_DSurfTempDt(i,j)
xy_DSurfSolDt(i,j) = xy_DSurfSolDtSave(i,j) - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
xy_DSurfLiqDt(i,j) = xy_DSurfLiqDtSave(i,j) + xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
!!$ if ( xy_SurfSnowB(i,j) + xy_DSurfSnowDt(i,j) * ( 2.0_DP * DelTime ) < 0.0_DP ) then
!!$ call MessageNotify( 'M', module_name, &
!!$ & 'Surface snow amount is negative %f, %f.', &
!!$ & d = (/ xy_SurfSnowB(i,j) + xy_DSurfSnowDt(i,j) * ( 2.0_DP * DelTime ), xy_SurfSnowB(i,j) /) )
!!$ end if
else
xy_LatHeatFluxBySnowMelt(i,j) = 0.0_DP
end if
end do
end do
!----------
! A case that all snow melt or soil moisture freeze
!----------
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
if ( xy_FlagCalc(i,j) ) then
SurfSolATentative = xy_SurfSolB(i,j) + xy_DSurfSolDt(i,j) * 2.0_DP * DelTime
if ( SurfSolATentative < 0.0_DP ) then
xy_FlagCalc(i,j) = .true.
xy_IndexMeltOrFreeze(i,j) = IndexMelt
else
xy_FlagCalc(i,j) = .false.
xy_IndexMeltOrFreeze(i,j) = IndexOthers
end if
!!$ select case ( IndexSpc )
!!$ case ( IndexSpcH2O )
!!$ SurfLiqATentative = xy_SurfLiqB(i,j) &
!!$ & + xy_DSurfLiqDt(i,j) * 2.0_DP * DelTime
!!$ if ( SurfLiqATentative < 0.0_DP ) then
!!$ xy_FlagCalc(i,j) = .true.
!!$ xy_IndexMeltOrFreeze(i,j) = IndexFreeze
!!$ end if
!!$ end select
else
xy_FlagCalc(i,j) = .false.
xy_IndexMeltOrFreeze(i,j) = IndexOthers
end if
case default
xy_FlagCalc(i,j) = .false.
xy_IndexMeltOrFreeze(i,j) = IndexOthers
end select
end do
end do
xyza_TempMtx = xyza_ArgTempMtx
xyz_TempVec = xyz_ArgTempVec
!
xyaa_SurfMtx = xyaa_ArgSurfMtx
xy_SurfRH = xy_ArgSurfRH
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
!!$ SurfSnowATentative = xy_SurfSnowB(i,j) &
!!$ & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
select case ( xy_IndexMeltOrFreeze(i,j) )
case ( IndexMelt )
! all ice/snow melt
DelSurfSol = xy_SurfSolATentativeSave(i,j)
case ( IndexFreeze )
! all soil moisture freeze (= negative melt of ice/snow)
DelSurfSol = - xy_SurfLiqATentativeSave(i,j)
end select
xy_SurfRH(i,j) = - xyr_RadSFlux(i,j,0) - xyr_RadLFlux(i,j,0) - xyr_HeatFlux(i,j,0) - xy_SurfLatentHeatFlux(i,j) + xyr_SoilHeatFlux(i,j,0) - xy_LatHeatFluxByOtherSpc(i,j) - LatentHeatLocal * DelSurfSol / ( 2.0_DP * DelTime )
end if
end do
end do
!
xyaa_SoilTempMtx = xyaa_ArgSoilTempMtx
xya_SoilTempVec = xya_ArgSoilTempVec
! 温度の計算
! Calculate temperature and specific humidity
!
do l = -1, 1
do k = 1, kslmax
xyaa_TempSoilTempLUMtx(:,:,-k,-l) = xyaa_SoilTempMtx(:,:,k,l)
end do
k = 0
xyaa_TempSoilTempLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
do k = 1, kmax
xyaa_TempSoilTempLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
end do
end do
!
call PhyImplLUDecomp3( xyaa_TempSoilTempLUMtx, imax * jmax, kmax + 1 + kslmax )
!
do k = 1, kslmax
xya_DelTempSoilTempLUVec(:,:,-k) = xya_SoilTempVec(:,:,k)
end do
k = 0
xya_DelTempSoilTempLUVec(:,:,k) = xy_SurfRH
do k = 1, kmax
xya_DelTempSoilTempLUVec(:,:,k) = xyz_TempVec(:,:,k)
end do
!
call PhyImplLUSolve3( xya_DelTempSoilTempLUVec, xyaa_TempSoilTempLUMtx, 1, imax * jmax , kmax + 1 + kslmax )
do k = 1, kslmax
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
xyz_DSoilTempDt(i,j,k) = xya_DelTempSoilTempLUVec(i,j,-k) / ( 2.0_DP * DelTime )
case default
xyz_DSoilTempDt(i,j,k) = 0.0_DP
end select
end if
end do
end do
end do
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$ case ( IndexSeaIce )
!!$ ! sea ice
!!$ xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$ case ( IndexSlabOcean )
!!$ ! slab ocean
!!$ xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$ case ( IndexOceanPresSST )
!!$ ! open ocean
!!$ xy_DSurfTempDt(i,j) = 0.0_DP
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end if
end do
end do
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
xyz_DTempDt(i,j,k) = xya_DelTempSoilTempLUVec(i,j,k) / ( 2.0_DP * DelTime )
end if
end do
end do
end do
!
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
!!$ SurfSnowATentative = xy_SurfSnowB(i,j) &
!!$ & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
select case ( xy_IndexMeltOrFreeze(i,j) )
case ( IndexMelt )
! all ice/snow melt
DelSurfSol = xy_SurfSolATentativeSave(i,j)
case ( IndexFreeze )
! all soil moisture freeze (= negative melt of ice/snow)
DelSurfSol = - xy_SurfLiqATentativeSave(i,j)
end select
xy_LatHeatFluxBySnowMelt(i,j) = LatentHeatLocal * DelSurfSol / ( 2.0_DP * DelTime )
xy_DSurfSolDt(i,j) = xy_DSurfSolDtSave(i,j) - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
xy_DSurfLiqDt(i,j) = xy_DSurfLiqDtSave(i,j) + xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
end if
end do
end do
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ if ( xy_FlagCalc(i,j) ) then
!!$ if ( xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * ( 2.0_DP * DelTime ) < xy_TempCond(i,j) ) then
!!$ call MessageNotify( 'M', module_name, &
!!$ & 'Surface temperature is lower than condensation temperature, %f < %f.', &
!!$ & d = (/ xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * ( 2.0_DP * DelTime ), xy_TempCond(i,j) /) )
!!$ end if
!!$ end if
!!$ end do
!!$ end do
!----------
! Calculation for a land point with prescribed temperature
!
!----------
! Surface fluxes used below
!----------
!!$ xy_SurfRadSFlux = xyr_RadSFlux(:,:,0)
!!$ xy_SurfRadLFlux = xyr_RadLFlux(:,:,0) &
!!$ & + xyra_DelRadLFlux(:,:,0,0) * xy_DSurfTempDt * ( 2.0_DP * DelTime ) &
!!$ & + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0_DP * DelTime )
!!$ xy_SurfSoilHeatCondFlux = xyr_SoilHeatFlux(:,:,0) &
!!$ & - xyr_SoilTempTransCoef(:,:,0) &
!!$ & * ( xyz_DSoilTempDt(:,:,1) - xy_DSurfTempDt ) * ( 2.0d0 * DelTime )
!!$ xy_SurfSensHeatFlux = &
!!$ & xyr_HeatFlux(:,:,0) &
!!$ & - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef &
!!$ & * ( xyz_DTempDt(:,:,1) / xyz_Exner(:,:,1) &
!!$ & - xy_DSurfTempDt / xyr_Exner(:,:,0) ) * ( 2.0_DP * DelTime )
!!$ !
!!$ xy_LatHeatFluxBySnowMelt = &
!!$ & - xy_SurfRadSFlux &
!!$ & - xy_SurfRadLFlux &
!!$ & - xy_SurfSensHeatFlux &
!!$ & - xy_SurfLatentHeatFlux &
!!$ & + xy_SurfSoilHeatCondFlux &
!!$ & - xy_LatHeatFluxByOtherSpc &
!!$ & - xy_SurfHeatCapacity * xy_DSurfTempDt
!
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLandWithPresTs )
SurfSolATentative = xy_SurfSolATentativeSave(i,j)
LatHeatFluxBySnowMelt = - xy_SurfRadSFlux(i,j) - xy_SurfRadLFlux(i,j) - xy_SurfSensHeatFlux(i,j) - xy_SurfLatentHeatFlux(i,j) - xy_LatHeatFluxByOtherSpc(i,j) - xy_SurfHeatCapacity(i,j) * xy_DSurfTempDt(i,j)
if ( SurfSolATentative > 0.0_DP ) then
! Ice exists on the ground.
! Calculation is performed only when freezing and melting
xy_FlagCalc(i,j) = .true.
xy_LatHeatFluxBySnowMelt(i,j) = LatHeatFluxBySnowMelt
!!$ else if ( xy_LatHeatFluxBySnowMelt(i,j) < 0.0_DP ) then
else if ( LatHeatFluxBySnowMelt < 0.0_DP ) then
! Ice does not exist on the ground.
! Calculation is performed only when freezing
xy_FlagCalc(i,j) = .true.
xy_LatHeatFluxBySnowMelt(i,j) = LatHeatFluxBySnowMelt
else
xy_FlagCalc(i,j) = .false.
end if
case default
xy_FlagCalc(i,j) = .false.
end select
end do
end do
!
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
SurfSolATentative = xy_SurfSolATentativeSave(i,j) - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal * ( 2.0_DP * DelTime )
if ( SurfSolATentative < 0.0_DP ) then
xy_LatHeatFluxBySnowMelt(i,j) = LatentHeatLocal * SurfSolATentative / ( 2.0_DP * DelTime )
end if
!!$ xy_DSurfSolDt(i,j) = xy_DSurfSolDtSave(i,j) &
!!$ & - SurfSolATentative / ( 2.0_DP * DelTime )
!!$ xy_DSurfLiqDt(i,j) = xy_DSurfLiqDtSave(i,j) &
!!$ & + SurfSolATentative / ( 2.0_DP * DelTime )
xy_DSurfSolDt(i,j) = xy_DSurfSolDtSave(i,j) - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
xy_DSurfLiqDt(i,j) = xy_DSurfLiqDtSave(i,j) + xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
end if
end do
end do
end if
!!$ ! 計算時間計測一時停止
!!$ ! Pause measurement of computation time
!!$ !
!!$ call TimesetClockStop( module_name )
end subroutine PhyImplSDHV3IceSnowPhaseChgCor
| Subroutine : | |||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , intent(in )
| ||
| xy_SeaIceThickness(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_Exner(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_Exner(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_SoilHeatFlux(0:imax-1, 1:jmax, 0:kslmax) : | real(DP), intent(in)
| ||
| xy_SurfTempTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_SoilTempTransCoef(0:imax-1, 1:jmax, 0:kslmax) : | real(DP), intent(in)
| ||
| xyr_RadSFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in )
| ||
| xyr_RadLFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in )
| ||
| xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(in )
| ||
| xy_SurfLatentHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xyza_ArgTempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1) : | real(DP), intent(in )
| ||
| xyz_ArgTempVec(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in )
| ||
| xyaa_ArgSurfMtx(0:imax-1, 1:jmax, 0:0, -1:1) : | real(DP), intent(in )
| ||
| xy_ArgSurfRH(0:imax-1,1:jmax) : | real(DP), intent(in ) | ||
| xyaa_ArgSoilTempMtx(0:imax-1, 1:jmax, 1:kslmax,-1:1) : | real(DP), intent(in )
| ||
| xya_ArgSoilTempVec(0:imax-1, 1:jmax, 1:kslmax) : | real(DP), intent(in )
| ||
| xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(inout)
| ||
| xy_DSurfTempDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xyz_DSoilTempDt(0:imax-1, 1:jmax, 1:kslmax) : | real(DP), intent(inout)
| ||
| xy_LatHeatFluxBySeaIceMelt(0:imax-1, 1:jmax) : | real(DP), intent(out ) |
融雪による時間変化率の修正を行います.
Correction of tendencies due to melt of snow.
subroutine PhyImplSDHV3SeaIceCorrection( xy_IndexCalcMethod, xy_SeaIceThickness, xyz_Exner, xyr_Exner, xy_SurfTemp, xyr_HeatFlux, xyr_SoilHeatFlux, xy_SurfTempTransCoef, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux, xy_SurfLatentHeatFlux, xyza_ArgTempMtx, xyz_ArgTempVec, xyaa_ArgSurfMtx, xy_ArgSurfRH, xyaa_ArgSoilTempMtx, xya_ArgSoilTempVec, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_LatHeatFluxBySeaIceMelt )
!
! 融雪による時間変化率の修正を行います.
!
! Correction of tendencies due to melt of snow.
!
! モジュール引用 ; USE statements
!
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! 物理定数設定
! Physical constants settings
!
use constants, only: CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 雪と海氷の定数の設定
! Setting constants of snow and sea ice
!
use constants_snowseaice, only: TempCondWater , SeaIceThermCondCoef, SeaIceVolHeatCap , TempBelowSeaIce
! 陰解法による時間積分のためのルーチン
! Routines for time integration with implicit scheme
!
use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3
! 宣言文 ; Declaration statements
!
integer , intent(in ) :: xy_IndexCalcMethod(0:imax-1, 1:jmax)
!
! Index for calculation method
real(DP), intent(in):: xy_SeaIceThickness(0:imax-1, 1:jmax)
!
! Sea ice thickness
real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
! Exner 関数 (整数レベル).
! Exner function (full level)
real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
! Exner 関数 (半整数レベル).
! Exner function (half level)
real(DP), intent(in ) :: xy_SurfTemp (0:imax-1, 1:jmax)
! 地表面温度.
! Surface temperature
real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(in):: xyr_SoilHeatFlux (0:imax-1, 1:jmax, 0:kslmax)
! 土壌中の熱フラックス (W m-2)
! Heat flux in sub-surface soil (W m-2)
real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in):: xyr_SoilTempTransCoef (0:imax-1, 1:jmax, 0:kslmax)
! 輸送係数：土壌温度.
! Transfer coefficient: soil temperature
real(DP), intent(in ):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
! 短波 (日射) フラックス.
! Shortwave (insolation) flux
real(DP), intent(in ):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
! 長波フラックス.
! Longwave flux
real(DP), intent(in ):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
! 長波地表温度変化.
! Surface temperature tendency with longwave
real(DP), intent(in ):: xy_SurfLatentHeatFlux(0:imax-1, 1:jmax)
! 惑星表面潜熱フラックス.
! Latent heat flux at the surface
real(DP), intent(in ):: xyza_ArgTempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
! 温度陰解行列.
! Implicit matrix about temperature
real(DP), intent(in ):: xyz_ArgTempVec(0:imax-1, 1:jmax, 1:kmax)
! 温度陰解ベクトル.
! Implicit vector about temperature
real(DP), intent(in ):: xyaa_ArgSurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
! 惑星表面エネルギー収支用陰解行列
! Implicit matrix for surface energy balance
real(DP), intent(in ):: xy_ArgSurfRH(0:imax-1,1:jmax)
real(DP), intent(in ):: xyaa_ArgSoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
! 土壌温度拡散方程式の行列
! Matrix for diffusion equation of soil temperature
real(DP), intent(in ):: xya_ArgSoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
! 土壌温度拡散方程式のベクトル
! Vector for diffusion equation of soil temperature
real(DP), intent(inout):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{T}{t} $ . 温度変化.
! Temperature tendency
real(DP), intent(inout) :: xy_DSurfTempDt (0:imax-1, 1:jmax)
! 地表面温度変化率 (K s-1)
! Surface temperature tendency (K s-1)
real(DP), intent(inout):: xyz_DSoilTempDt (0:imax-1, 1:jmax, 1:kslmax)
! $ \DP{Tg}{t} $ . 土壌温度変化 (K s-1)
! Temperature tendency (K s-1)
real(DP), intent(out ) :: xy_LatHeatFluxBySeaIceMelt(0:imax-1, 1:jmax)
! 作業変数
! Work variables
!
real(DP):: xyza_TempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
! 温度陰解行列.
! Implicit matrix about temperature
real(DP):: xyz_TempVec(0:imax-1, 1:jmax, 1:kmax)
! 温度陰解ベクトル.
! Implicit vector about temperature
real(DP):: xyaa_SurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
! 惑星表面エネルギー収支用陰解行列
! Implicit matrix for surface energy balance
real(DP):: xy_SurfRH(0:imax-1,1:jmax)
real(DP):: xyaa_SoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
! 土壌温度拡散方程式の行列
! Matrix for diffusion equation of soil temperature
real(DP):: xya_SoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
! 土壌温度拡散方程式のベクトル
! Vector for diffusion equation of soil temperature
real(DP):: xyaa_TempSoilTempLUMtx (0:imax-1, 1:jmax, -kslmax:kmax, -1:1)
! LU 行列.
! LU matrix
real(DP):: xya_DelTempSoilTempLUVec (0:imax-1, 1:jmax, -kslmax:kmax)
! $ T, Tg $ の時間変化.
! Tendency of $ T $ and $ Tg |
logical :: xy_FlagSeaIceMelt(0:imax-1, 1:jmax)
real(DP) :: xy_TempCond(0:imax-1, 1:jmax)
real(DP) :: SurfTempATentative
real(DP) :: xy_SurfRadSFlux (0:imax-1, 1:jmax)
real(DP) :: xy_SurfRadLFlux (0:imax-1, 1:jmax)
real(DP) :: xy_SurfSoilHeatCondFlux(0:imax-1, 1:jmax)
real(DP) :: xy_SurfSensHeatFlux (0:imax-1, 1:jmax)
real(DP) :: xy_SeaIceHeatCondFlux (0:imax-1, 1:jmax)
real(DP) :: xy_HeatingTendency (0:imax-1, 1:jmax)
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: l
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. phy_implicit_sdh_V3_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
!!$ ! 計算時間計測開始
!!$ ! Start measurement of computation time
!!$ !
!!$ call TimesetClockStart( module_name )
!
! check flag of snow melt
!
if ( ( .not. FlagBucketModel ) .or. ( .not. FlagSnow ) ) then
xy_LatHeatFluxBySeaIceMelt = 0.0_DP
return
end if
xy_TempCond = TempCondWater
do j = 1, jmax
do i = 0, imax-1
SurfTempATentative = xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
if ( ( xy_IndexCalcMethod(i,j) == IndexSeaIce ) .and. ( SurfTempATentative > xy_TempCond(i,j) ) ) then
xy_FlagSeaIceMelt(i,j) = .true.
else
xy_FlagSeaIceMelt(i,j) = .false.
end if
end do
end do
xyza_TempMtx = xyza_ArgTempMtx
xyz_TempVec = xyz_ArgTempVec
!
xyaa_SurfMtx = xyaa_ArgSurfMtx
xy_SurfRH = xy_ArgSurfRH
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagSeaIceMelt(i,j) ) then
xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
xyaa_SurfMtx(i,j,0, 0) = 1.0_DP
xyaa_SurfMtx(i,j,0, 1) = 0.0_DP
xy_SurfRH (i,j) = xy_TempCond(i,j) - xy_SurfTemp(i,j)
end if
end do
end do
!
xyaa_SoilTempMtx = xyaa_ArgSoilTempMtx
xya_SoilTempVec = xya_ArgSoilTempVec
! 温度の計算
! Calculate temperature and specific humidity
!
do l = -1, 1
do k = 1, kslmax
xyaa_TempSoilTempLUMtx(:,:,-k,-l) = xyaa_SoilTempMtx(:,:,k,l)
end do
k = 0
xyaa_TempSoilTempLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
do k = 1, kmax
xyaa_TempSoilTempLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
end do
end do
call PhyImplLUDecomp3( xyaa_TempSoilTempLUMtx, imax * jmax, kmax + 1 + kslmax )
do k = 1, kslmax
xya_DelTempSoilTempLUVec(:,:,-k) = xya_SoilTempVec(:,:,k)
end do
k = 0
xya_DelTempSoilTempLUVec(:,:,k) = xy_SurfRH
do k = 1, kmax
xya_DelTempSoilTempLUVec(:,:,k) = xyz_TempVec(:,:,k)
end do
call PhyImplLUSolve3( xya_DelTempSoilTempLUVec, xyaa_TempSoilTempLUMtx, 1, imax * jmax , kmax + 1 + kslmax )
do k = 1, kslmax
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagSeaIceMelt(i,j) ) then
select case ( xy_IndexCalcMethod(i,j) )
!!$ case ( IndexLand )
!!$ xyz_DSoilTempDt(i,j,k) = &
!!$ & xya_DelTempSoilTempLUVec(i,j,-k) / ( 2.0_DP * DelTime )
case default
xyz_DSoilTempDt(i,j,k) = 0.0_DP
end select
end if
end do
end do
end do
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagSeaIceMelt(i,j) ) then
select case ( xy_IndexCalcMethod(i,j) )
!!$ case ( IndexLand )
!!$ ! land
!!$ xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
case ( IndexSeaIce )
! sea ice
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$ case ( IndexSlabOcean )
!!$ ! slab ocean
!!$ xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
!!$ case ( IndexOceanPresSST )
!!$ ! open ocean
!!$ xy_DSurfTempDt(i,j) = 0.0_DP
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end if
end do
end do
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagSeaIceMelt(i,j) ) then
xyz_DTempDt(i,j,k) = xya_DelTempSoilTempLUVec(i,j,k) / ( 2.0_DP * DelTime )
end if
end do
end do
end do
xy_SurfRadSFlux = xyr_RadSFlux(:,:,0)
xy_SurfRadLFlux = xyr_RadLFlux(:,:,0) + xyra_DelRadLFlux(:,:,0,0) * xy_DSurfTempDt * ( 2.0_DP * DelTime ) + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0_DP * DelTime )
xy_SurfSoilHeatCondFlux = xyr_SoilHeatFlux(:,:,0) - xyr_SoilTempTransCoef(:,:,0) * ( xyz_DSoilTempDt(:,:,1) - xy_DSurfTempDt ) * ( 2.0_DP * DelTime )
xy_SurfSensHeatFlux = xyr_HeatFlux(:,:,0) - CpDry * xyr_Exner(:,:,0) * xy_SurfTempTransCoef * ( xyz_DTempDt(:,:,1) / xyz_Exner(:,:,1) - xy_DSurfTempDt / xyr_Exner(:,:,0) ) * ( 2.0_DP * DelTime )
xy_SeaIceHeatCondFlux = - SeaIceThermCondCoef * ( xy_SurfTemp + xy_DSurfTempDt * ( 2.0_DP * DelTime ) - TempBelowSeaIce ) / xy_SeaIceThickness
xy_HeatingTendency = SeaIceVolHeatCap * xy_SeaIceThickness * xy_DSurfTempDt
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagSeaIceMelt(i,j) ) then
xy_LatHeatFluxBySeaIceMelt(i,j) = - ( ( xy_SurfRadSFlux(i,j) + xy_SurfRadLFlux(i,j) + xy_SurfSensHeatFlux(i,j) + xy_SurfLatentHeatFlux(i,j) ) - xy_SeaIceHeatCondFlux(i,j) ) - xy_HeatingTendency(i,j)
else
xy_LatHeatFluxBySeaIceMelt(i,j) = 0.0_DP
end if
end do
end do
!!$ ! 計算時間計測一時停止
!!$ ! Pause measurement of computation time
!!$ !
!!$ call TimesetClockStop( module_name )
end subroutine PhyImplSDHV3SeaIceCorrection
| Subroutine : | |||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , intent(in)
| ||
| xy_SeaIceThickness(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xy_SurfH2OVapFlux(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfLatentHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_SoilHeatFlux(0:imax-1, 1:jmax, 0:kslmax) : | real(DP), intent(in)
| ||
| xyr_RadSFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_RadLFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_DeepSubSurfHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_SoilTemp(0:imax-1, 1:jmax, 1:kslmax) : | real(DP), intent(in)
| ||
| xy_SurfHumidCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfHeatCapacity(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SoilHeatCap(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xy_SoilHeatDiffCoef(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(in)
| ||
| xyr_Press(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyz_Exner(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_Exner(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_SurfVelTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfTempTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfQVapTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_SoilTempTransCoef(0:imax-1, 1:jmax, 0:kslmax) : | real(DP), intent(in)
| ||
| xy_SurfMajCompIceB(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SoilMoistB(0:imax-1, 1:jmax) : | real(DP), intent(in) | ||
| xy_SurfSnowB(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(out)
| ||
| xy_DSurfTempDt(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xyz_DSoilTempDt(0:imax-1, 1:jmax, 1:kslmax) : | real(DP), intent(out)
| ||
| xy_DPsDt(0:imax-1, 1:jmax) : | real(DP), intent(out) | ||
| xy_DSurfMajCompIceDt(0:imax-1, 1:jmax) : | real(DP), intent(out) | ||
| xy_DSoilMoistDt(0:imax-1, 1:jmax) : | real(DP), intent(out)
| ||
| xy_DSurfSnowDt(0:imax-1, 1:jmax) : | real(DP), intent(out)
|
時間変化率の計算を行います.
Calculate tendencies.
subroutine PhyImplSDHV3TendencyHeatCore( xy_IndexCalcMethod, xy_SeaIceThickness, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfHumidCoef, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xyra_DelRadLFlux, xyr_Press, xyz_Exner, xyr_Exner, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xyr_SoilTempTransCoef, xy_SurfMajCompIceB, xy_SoilMoistB, xy_SurfSnowB, xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DPsDt, xy_DSurfMajCompIceDt, xy_DSoilMoistDt, xy_DSurfSnowDt )
!
! 時間変化率の計算を行います.
!
! Calculate tendencies.
!
! モジュール引用 ; USE statements
!
! 座標データ設定
! Axes data settings
!
use axesset, only: r_SSDepth, z_SSDepth ! subsurface grid at midpoint of layer
! 物理定数設定
! Physical constants settings
!
use constants, only: Grav, CpDry, GasRDry
! $ R $ [J kg-1 K-1].
! 乾燥大気の気体定数.
! Gas constant of air
! 雪と海氷の定数の設定
! Setting constants of snow and sea ice
!
use constants_snowseaice, only: TempCondWater, SeaIceVolHeatCap , SeaIceThermCondCoef, TempBelowSeaIce
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! 陰解法による時間積分のためのルーチン
! Routines for time integration with implicit scheme
!
use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3
! 宣言文 ; Declaration statements
!
integer , intent(in):: xy_IndexCalcMethod (0:imax-1, 1:jmax)
!
! Index for calculation method
real(DP), intent(in):: xy_SeaIceThickness(0:imax-1, 1:jmax)
!
! Sea ice thickness
real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
! 比湿フラックス.
! Specific humidity flux
real(DP), intent(in):: xy_SurfH2OVapFlux(0:imax-1, 1:jmax)
! 惑星表面水蒸気フラックス.
! Water vapor flux at the surface
real(DP), intent(in):: xy_SurfLatentHeatFlux(0:imax-1, 1:jmax)
! 惑星表面潜熱フラックス.
! Latent heat flux at the surface
real(DP), intent(in):: xyr_SoilHeatFlux (0:imax-1, 1:jmax, 0:kslmax)
! 土壌中の熱フラックス (W m-2)
! Heat flux in sub-surface soil (W m-2)
real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
! 短波 (日射) フラックス.
! Shortwave (insolation) flux
real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
! 長波フラックス.
! Longwave flux
real(DP), intent(in):: xy_DeepSubSurfHeatFlux (0:imax-1, 1:jmax)
! 地中熱フラックス.
! "Deep subsurface heat flux"
! Heat flux at the bottom of surface/soil layer.
real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
! 地表面温度.
! Surface temperature
real(DP), intent(in):: xyz_SoilTemp (0:imax-1, 1:jmax, 1:kslmax)
! 土壌温度 (K)
! Soil temperature (K)
real(DP), intent(in):: xy_SurfHumidCoef (0:imax-1, 1:jmax)
! 地表湿潤度.
! Surface humidity coefficient
real(DP), intent(in):: xy_SurfHeatCapacity (0:imax-1, 1:jmax)
! 地表熱容量.
! Surface heat capacity
real(DP), intent(in ):: xy_SoilHeatCap (0:imax-1, 1:jmax)
! 土壌熱容量 (J K-1 kg-1)
! Specific heat of soil (J K-1 kg-1)
real(DP), intent(in ):: xy_SoilHeatDiffCoef (0:imax-1, 1:jmax)
! 土壌熱伝導係数 (J m-3 K-1)
! Heat conduction coefficient of soil (J m-3 K-1)
real(DP), intent(in):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
! 長波地表温度変化.
! Surface temperature tendency with longwave
real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air pressure (half level)
real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
! Exner 関数 (整数レベル).
! Exner function (full level)
real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
! Exner 関数 (半整数レベル).
! Exner function (half level)
real(DP), intent(in):: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：運動量.
! Transfer coefficient: velocity
real(DP), intent(in):: xyr_TempTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：質量.
! Transfer coefficient: mass of constituents
real(DP), intent(in):: xy_SurfVelTransCoef (0:imax-1, 1:jmax)
! 輸送係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in):: xy_SurfQVapTransCoef (0:imax-1, 1:jmax)
! 輸送係数：比湿.
! Transfer coefficient: specific humidity
real(DP), intent(in):: xyr_SoilTempTransCoef (0:imax-1, 1:jmax, 0:kslmax)
! 輸送係数：土壌温度.
! Transfer coefficient: soil temperature
real(DP), intent(in):: xy_SurfMajCompIceB (0:imax-1, 1:jmax)
!
! Surface major component ice amount.
real(DP), intent(in):: xy_SoilMoistB(0:imax-1, 1:jmax)
real(DP), intent(in):: xy_SurfSnowB (0:imax-1, 1:jmax)
! 積雪量.
! Surface snow amount.
real(DP), intent(out):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{T}{t} $ . 温度変化.
! Temperature tendency
real(DP), intent(out):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ \DP{q}{t} $ . 質量混合比変化.
! Mass mixing ratio tendency
real(DP), intent(out):: xy_DSurfTempDt (0:imax-1, 1:jmax)
! 地表面温度変化率 (K s-1)
! Surface temperature tendency (K s-1)
real(DP), intent(out):: xyz_DSoilTempDt (0:imax-1, 1:jmax, 1:kslmax)
! $ \DP{Tg}{t} $ . 土壌温度変化 (K s-1)
! Temperature tendency (K s-1)
real(DP), intent(out):: xy_DPsDt (0:imax-1, 1:jmax)
real(DP), intent(out):: xy_DSurfMajCompIceDt(0:imax-1, 1:jmax)
real(DP), intent(out):: xy_DSoilMoistDt (0:imax-1, 1:jmax)
! 土壌温度時間変化率 (kg m-2 s-1)
! Soil temperature tendency (kg m-2 s-1)
real(DP), intent(out):: xy_DSurfSnowDt (0:imax-1, 1:jmax)
! 積雪率時間変化率 (kg m-2 s-1)
! Surface snow amount tendency (kg m-2 s-1)
! 作業変数
! Work variables
!
real(DP):: xyza_TempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
! 温度陰解行列.
! Implicit matrix about temperature
real(DP):: xyz_TempVec(0:imax-1, 1:jmax, 1:kmax)
! 温度陰解ベクトル.
! Implicit vector about temperature
real(DP):: xyza_QMixMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
! 質量混合比陰解行列.
! Implicit matrix about mass mixing ratio
real(DP):: xyzf_QMixVec(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! 質量混合比陰解ベクトル.
! Implicit vector about mass mixing ratio
real(DP):: xyaa_SurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
! 惑星表面エネルギー収支用陰解行列
! Implicit matrix for surface energy balance
real(DP):: xy_SurfRH(0:imax-1,1:jmax)
!!$ real(DP):: xyza_TempQVapLUMtx (0:imax-1, 1:jmax, -kmax:kmax, -1:1)
!!$ ! LU 行列.
!!$ ! LU matrix
!!$ real(DP):: xyz_DelTempQVap (0:imax-1, 1:jmax, -kmax:kmax)
!!$ ! $ T q $ の時間変化.
!!$ ! Tendency of $ T q $
!!$
!!$ real(DP):: xyza_TempLUMtx (0:imax-1, 1:jmax, 0:kmax, -1:1)
!!$ ! LU 行列.
!!$ ! LU matrix
!!$ real(DP):: xyz_DelTempLUVec (0:imax-1, 1:jmax, 0:kmax)
!!$ ! $ T q $ の時間変化.
!!$ ! Tendency of $ T q $
real(DP):: xyza_QMixLUMtx (0:imax-1, 1:jmax, 1:kmax, -1:1)
! LU 行列.
! LU matrix
real(DP):: xyz_DelQMixLUVec (0:imax-1, 1:jmax, 1:kmax)
! $ q $ の時間変化.
! Tendency of $ q $
real(DP):: xyaa_SoilTempMtx (0:imax-1, 1:jmax, 1:kslmax,-1:1)
! 土壌温度拡散方程式の行列
! Matrix for diffusion equation of soil temperature
real(DP):: xya_SoilTempVec (0:imax-1, 1:jmax, 1:kslmax)
! 土壌温度拡散方程式のベクトル
! Vector for diffusion equation of soil temperature
real(DP):: xyaa_TempSoilTempLUMtx (0:imax-1, 1:jmax, -kslmax:kmax, -1:1)
! LU 行列.
! LU matrix
real(DP):: xya_DelTempSoilTempLUVec (0:imax-1, 1:jmax, -kslmax:kmax)
! $ T, Tg $ の時間変化.
! Tendency of $ T $ and $ Tg |
real(DP):: SurfSnowATentative
! 積雪量の仮の値 (kg m-2)
! pseudo value of surface snow amount (kg m-2)
real(DP):: xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax)
!
! Latent heat flux by major component ice sublimation
! (variable only for debug)
real(DP):: xy_LatHeatFluxBySnowMelt(0:imax-1, 1:jmax)
!
! Latent heat flux by melt
! (variable only for debug)
real(DP):: xy_LatHeatFluxBySeaIceMelt(0:imax-1, 1:jmax)
!
! Latent heat flux by sea ice melt
! (variable only for debug)
real(DP):: xy_SurfMajCompLiqB (0:imax-1, 1:jmax)
real(DP):: xy_LatHeatFluxByOtherSpc(0:imax-1, 1:jmax)
real(DP):: xy_DAtmMassDt(0:imax-1, 1:jmax)
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: l ! 行列用 DO ループ用作業変数
! Work variables for DO loop of matrices
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. phy_implicit_sdh_V3_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
!!$ ! 計算時間計測開始
!!$ ! Start measurement of computation time
!!$ !
!!$ call TimesetClockStart( module_name )
!!$ if ( .not. FlagSSModel ) then
!!$ call MessageNotify( 'E', module_name, 'FlagSSModel has to be true.' )
!!$ end if
! FlagBucketModel は関係ないよね?
! SSModel 強制にした時点で, 水蒸気は地面と分離したから.
!!$ if ( .not. FlagBucketModel ) then
!!$ call MessageNotify( 'E', module_name, 'FlagBucketModel has to be true.' )
!!$ end if
! 陰解法のための行列作成
! Create matrices for implicit scheme
!
! 鉛直拡散スキームの輸送係数から陰解行列の計算 (温度)
! Calculate implicit matrices from transfer coefficient of vertical diffusion scheme (temperature)
!
k = 1
xyza_TempMtx(:,:,k,-1) = - CpDry * xy_SurfTempTransCoef(:,:)
xyza_TempMtx(:,:,k, 0) = - CpDry * ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k ) * xy_SurfTempTransCoef(:,:) + CpDry * xyr_Exner(:,:,k ) / xyz_Exner(:,:,k ) * xyr_TempTransCoef(:,:,k )
xyza_TempMtx(:,:,k, 1) = - CpDry * xyr_Exner(:,:,k ) / xyz_Exner(:,:,k+1) * xyr_TempTransCoef(:,:,k )
do k = 2, kmax-1
xyza_TempMtx(:,:,k,-1) = - CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k-1) * xyr_TempTransCoef(:,:,k-1)
xyza_TempMtx(:,:,k, 0) = - CpDry * ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k ) * xyr_TempTransCoef(:,:,k-1) + CpDry * xyr_Exner(:,:,k ) / xyz_Exner(:,:,k ) * xyr_TempTransCoef(:,:,k )
xyza_TempMtx(:,:,k, 1) = - CpDry * xyr_Exner(:,:,k ) / xyz_Exner(:,:,k+1) * xyr_TempTransCoef(:,:,k )
end do
k = kmax
xyza_TempMtx(:,:,k,-1) = - CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k-1) * xyr_TempTransCoef(:,:,k-1)
xyza_TempMtx(:,:,k, 0) = - CpDry * ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k ) * xyr_TempTransCoef(:,:,k-1)
xyza_TempMtx(:,:,k, 1) = 0.0_DP
do k = 1, kmax
xyz_TempVec(:,:,k) = - ( xyr_HeatFlux(:,:,k) - xyr_HeatFlux(:,:,k-1) )
end do
! 鉛直拡散スキームの輸送係数から陰解行列の計算 (比湿)
! Calculate implicit matrices from transfer coefficient of vertical diffusion scheme (specific humidity)
!
k = 1
xyza_QMixMtx(:,:,k,-1) = 0.0_DP
xyza_QMixMtx(:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xyr_QMixTransCoef(:,:,k )
xyza_QMixMtx(:,:,k, 1) = - xyr_QMixTransCoef(:,:,k )
do k = 2, kmax-1
xyza_QMixMtx(:,:,k,-1) = - xyr_QMixTransCoef(:,:,k-1)
xyza_QMixMtx(:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xyr_QMixTransCoef(:,:,k-1) + xyr_QMixTransCoef(:,:,k )
xyza_QMixMtx(:,:,k, 1) = - xyr_QMixTransCoef(:,:,k )
end do
k = kmax
xyza_QMixMtx(:,:,k,-1) = - xyr_QMixTransCoef(:,:,k-1)
xyza_QMixMtx(:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xyr_QMixTransCoef(:,:,k-1)
xyza_QMixMtx(:,:,k, 1) = 0.0_DP
do n = 1, ncmax
if ( n == IndexH2OVap ) then
do k = 1, 1
xyzf_QMixVec(:,:,k,n) = - ( xyrf_QMixFlux(:,:,k,n) - xy_SurfH2OVapFlux )
end do
do k = 1+1, kmax
xyzf_QMixVec(:,:,k,n) = - ( xyrf_QMixFlux(:,:,k,n) - xyrf_QMixFlux(:,:,k-1,n) )
end do
else
do k = 1, kmax
xyzf_QMixVec(:,:,k,n) = - ( xyrf_QMixFlux(:,:,k,n) - xyrf_QMixFlux(:,:,k-1,n) )
end do
end if
end do
! 土壌温度計算用の輸送係数から陰解行列の計算 (土壌温度)
! Calculate implicit matrices by using transfer coefficient (soil temperature)
!
if ( kslmax /= 0 ) then ! xyr_SoilTempMtx is not used when kslmax = 0.
do k = 1, kslmax-1
xyaa_SoilTempMtx(:,:,k,-1) = - xyr_SoilTempTransCoef(:,:,k-1)
xyaa_SoilTempMtx(:,:,k, 0) = xy_SoilHeatCap(:,:) * ( r_SSDepth(k) - r_SSDepth(k-1) ) / ( 2.0_DP * DelTime ) + xyr_SoilTempTransCoef(:,:,k-1) + xyr_SoilTempTransCoef(:,:,k )
xyaa_SoilTempMtx(:,:,k, 1) = - xyr_SoilTempTransCoef(:,:,k )
end do
k = kslmax
xyaa_SoilTempMtx(:,:,k,-1) = - xyr_SoilTempTransCoef(:,:,k-1)
xyaa_SoilTempMtx(:,:,k, 0) = xy_SoilheatCap(:,:) * ( r_SSDepth(k) - r_SSDepth(k-1) ) / ( 2.0_DP * DelTime ) + xyr_SoilTempTransCoef(:,:,k-1)
xyaa_SoilTempMtx(:,:,k, 1) = 0.0_DP
end if
do k = 1, kslmax
xya_SoilTempVec (:,:,k) = - ( xyr_SoilHeatFlux(:,:,k) - xyr_SoilHeatFlux(:,:,k-1) )
end do
! 地表面過程の輸送係数から陰解行列の計算
! Calculate implicit matrices from transfer coefficient of surface process
!
do i = 0, imax-1
do j = 1, jmax
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
xyaa_SurfMtx(i,j,0,-1) = xyr_SoilTempTransCoef(i,j,0)
xyaa_SurfMtx(i,j,0, 0) = xy_SurfHeatCapacity(i,j) / ( 2.0_DP * DelTime ) + CpDry * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,0) - xyr_SoilTempTransCoef(i,j,0)
xyaa_SurfMtx(i,j,0, 1) = - CpDry * xyr_Exner(i,j,0) / xyz_Exner(i,j,1) * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,1)
case ( IndexSeaIce )
! sea ice
xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
xyaa_SurfMtx(i,j,0, 0) = SeaIceVolHeatCap * xy_SeaIceThickness(i,j) / ( 2.0_DP * DelTime ) + CpDry * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,0) + SeaIceThermCondCoef / xy_SeaIceThickness(i,j)
xyaa_SurfMtx(i,j,0, 1) = - CpDry * xyr_Exner(i,j,0) / xyz_Exner(i,j,1) * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,1)
case ( IndexSlabOcean )
! slab ocean ocean
xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
xyaa_SurfMtx(i,j,0, 0) = SOHeatCapacity / ( 2.0_DP * DelTime ) + CpDry * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,0)
xyaa_SurfMtx(i,j,0, 1) = - CpDry * xyr_Exner(i,j,0) / xyz_Exner(i,j,1) * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,1)
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
xyaa_SurfMtx(i,j,0, 0) = 1.0_DP
xyaa_SurfMtx(i,j,0, 1) = 0.0_DP
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
xy_SurfRH(i,j) = - xyr_RadSFlux(i,j,0) - xyr_RadLFlux(i,j,0) - xyr_HeatFlux(i,j,0) - xy_SurfLatentHeatFlux(i,j) + xyr_SoilHeatFlux(i,j,0)
case ( IndexSeaIce )
! sea ice
xy_SurfRH(i,j) = - xyr_RadSFlux(i,j,0) - xyr_RadLFlux(i,j,0) - xyr_HeatFlux(i,j,0) - xy_SurfLatentHeatFlux(i,j) - SeaIceThermCondCoef * ( xy_SurfTemp(i,j) - TempBelowSeaIce ) / xy_SeaIceThickness(i,j)
case ( IndexSlabOcean )
! slab ocean
xy_SurfRH(i,j) = - xyr_RadSFlux(i,j,0) - xyr_RadLFlux(i,j,0) - xyr_HeatFlux(i,j,0) - xy_SurfLatentHeatFlux(i,j) !&
! & + xy_DeepSubSurfHeatFlux(i,j)
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
xy_SurfRH(i,j) = 0.0_DP
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
! 温度と比湿の計算
! Calculate temperature and specific humidity
!
do l = -1, 1
do k = 1, kslmax
xyaa_TempSoilTempLUMtx(:,:,-k,-l) = xyaa_SoilTempMtx(:,:,k,l)
end do
k = 0
xyaa_TempSoilTempLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
do k = 1, kmax
xyaa_TempSoilTempLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
end do
end do
call PhyImplLUDecomp3( xyaa_TempSoilTempLUMtx, imax * jmax, kmax + 1 + kslmax )
do k = 1, kslmax
xya_DelTempSoilTempLUVec(:,:,-k) = xya_SoilTempVec(:,:,k)
end do
k = 0
xya_DelTempSoilTempLUVec(:,:,k) = xy_SurfRH
do k = 1, kmax
xya_DelTempSoilTempLUVec(:,:,k) = xyz_TempVec(:,:,k)
end do
call PhyImplLUSolve3( xya_DelTempSoilTempLUVec, xyaa_TempSoilTempLUMtx, 1, imax * jmax , kmax + 1 + kslmax )
do k = 1, kslmax
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
xyz_DSoilTempDt(i,j,k) = xya_DelTempSoilTempLUVec(i,j,-k) / ( 2.0_DP * DelTime )
case default
xyz_DSoilTempDt(i,j,k) = 0.0_DP
end select
end do
end do
end do
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
case ( IndexSeaIce )
! sea ice
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
case ( IndexSlabOcean )
! slab ocean
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2.0_DP * DelTime )
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
xy_DSurfTempDt(i,j) = 0.0_DP
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
do k = 1, kmax
xyz_DTempDt(:,:,k) = xya_DelTempSoilTempLUVec(:,:,k) / ( 2.0_DP * DelTime )
end do
!
! Calculation of tendencies of soil moisture and surface snow amount
!
if ( FlagBucketModel ) then
if ( FlagSnow ) then
! Evaporation is subtracted from surface snow and soil moisture
!
do j = 1, jmax
do i = 0, imax-1
!!$ if ( xyrf_QMixFlux(i,j,0,IndexH2OVap) >= 0.0_DP ) then
if ( xy_SurfH2OVapFlux(i,j) >= 0.0_DP ) then
!!$ xy_DSurfSnowDt(i,j) = - xyrf_QMixFlux(i,j,0,IndexH2OVap)
xy_DSurfSnowDt(i,j) = - xy_SurfH2OVapFlux(i,j)
SurfSnowATentative = xy_SurfSnowB(i,j) + xy_DSurfSnowDt(i,j) * 2.0_DP * DelTime
if ( SurfSnowATentative < 0.0_DP ) then
xy_DSoilMoistDt(i,j) = SurfSnowATentative / ( 2.0_DP * DelTime )
xy_DSurfSnowDt (i,j) = - xy_SurfSnowB(i,j) / ( 2.0_DP * DelTime )
else
xy_DSoilMoistDt(i,j) = 0.0_DP
end if
else
if ( xy_SurfSnowB(i,j) > 0.0_DP ) then
!!$ xy_DSurfSnowDt (i,j) = - xyrf_QMixFlux(i,j,0,IndexH2OVap)
xy_DSurfSnowDt (i,j) = - xy_SurfH2OVapFlux(i,j)
xy_DSoilMoistDt(i,j) = 0.0_DP
else
xy_DSurfSnowDt (i,j) = 0.0_DP
!!$ xy_DSoilMoistDt(i,j) = - xyrf_QMixFlux(i,j,0,IndexH2OVap)
xy_DSoilMoistDt(i,j) = - xy_SurfH2OVapFlux(i,j)
end if
end if
end do
end do
else
! Evaporation is subtracted from soil moisture
!
!!$ xy_DSoilMoistDt = - xyrf_QMixFlux(:,:,0,IndexH2OVap)
xy_DSoilMoistDt = - xy_SurfH2OVapFlux
xy_DSurfSnowDt = 0.0_DP
end if
else
xy_DSoilMoistDt = 0.0_DP
xy_DSurfSnowDt = 0.0_DP
end if
! Temporarily set
!
xy_DSurfMajCompIceDt = 0.0_DP
if ( FlagMajCompPhaseChange ) then
xy_DAtmMassDt = 0.0_DP
xy_DSurfMajCompIceDt = 0.0_DP
! Dummy values
!
xy_SurfMajCompLiqB = 0.0_DP
xy_LatHeatFluxByOtherSpc = 0.0_DP
call PhyImplSDHV3IceSnowPhaseChgCor( IndexSpcMajComp, xyr_Press(:,:,0), xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfMajCompLiqB, xy_SurfMajCompIceB, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xy_IndexCalcMethod, xyra_DelRadLFlux, xyz_Exner, xyr_Exner, xy_SurfTempTransCoef, xy_LatHeatFluxByOtherSpc, xyza_TempMtx, xyz_TempVec, xyaa_SurfMtx, xy_SurfRH, xyaa_SoilTempMtx, xya_SoilTempVec, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DAtmMassDt, xy_DSurfMajCompIceDt, xy_LatHeatFluxByMajCompIceSubl )
else
xy_DAtmMassDt = 0.0_DP
xy_LatHeatFluxByMajCompIceSubl = 0.0_DP
end if
xy_DPsDt = xy_DAtmMassDt * Grav
xy_LatHeatFluxByOtherSpc = xy_LatHeatFluxByMajCompIceSubl
if ( FlagSublimation ) then
! If sublimation is considered, the melt of snow/ice is not calculated.
xy_LatHeatFluxBySnowMelt = 0.0_DP
else
! Else, the melt of snow/ice is calculated.
call PhyImplSDHV3IceSnowPhaseChgCor( IndexSpcH2O, xyr_Press(:,:,0), xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SoilMoistB, xy_SurfSnowB, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xy_IndexCalcMethod, xyra_DelRadLFlux, xyz_Exner, xyr_Exner, xy_SurfTempTransCoef, xy_LatHeatFluxByOtherSpc, xyza_TempMtx, xyz_TempVec, xyaa_SurfMtx, xy_SurfRH, xyaa_SoilTempMtx, xya_SoilTempVec, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DSoilMoistDt, xy_DSurfSnowDt, xy_LatHeatFluxBySnowMelt )
end if
call PhyImplSDHV3SeaIceCorrection( xy_IndexCalcMethod, xy_SeaIceThickness, xyz_Exner, xyr_Exner, xy_SurfTemp, xyr_HeatFlux, xyr_SoilHeatFlux, xy_SurfTempTransCoef, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux, xy_SurfLatentHeatFlux, xyza_TempMtx, xyz_TempVec, xyaa_SurfMtx, xy_SurfRH, xyaa_SoilTempMtx, xya_SoilTempVec, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_LatHeatFluxBySeaIceMelt )
do l = -1, 1
do k = 1, kmax
xyza_QMixLUMtx(:,:,k,l) = xyza_QMixMtx(:,:,k,l)
end do
end do
call PhyImplLUDecomp3( xyza_QMixLUMtx, imax * jmax, kmax )
do n = 1, ncmax
do k = 1, kmax
xyz_DelQMixLUVec(:,:,k) = xyzf_QMixVec(:,:,k,n)
end do
call PhyImplLUSolve3( xyz_DelQMixLUVec, xyza_QMixLUMtx, 1, imax * jmax , kmax )
do k = 1, kmax
xyzf_DQMixDt(:,:,k,n) = xyz_DelQMixLUVec(:,:,k) / ( 2.0_DP * DelTime )
end do
end do
! Debug routine
!
call PhyImplSDHV3ChkConservation( xy_IndexCalcMethod, xy_SeaIceThickness, xyr_Press, xyz_Exner, xyr_Exner, xy_SurfTemp, xy_SurfHeatCapacity, xy_SoilHeatCap, xyr_HeatFlux, xy_SurfH2OVapFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xy_SurfTempTransCoef, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux, xy_LatHeatFluxByMajCompIceSubl, xy_LatHeatFluxBySnowMelt, xy_LatHeatFluxBySeaIceMelt, xy_DeepSubSurfHeatFlux, xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DSoilMoistDt, xy_DSurfSnowDt, xy_DPsDt, xy_DSurfMajCompIceDt )
!!$ ! 計算時間計測一時停止
!!$ ! Pause measurement of computation time
!!$ !
!!$ call TimesetClockStop( module_name )
end subroutine PhyImplSDHV3TendencyHeatCore
| Subroutine : | |||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , intent(in)
| ||
| xy_SeaIceThickness(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SnowFrac(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_HeatFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) : | real(DP), intent(in)
| ||
| xy_SurfSoilHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyr_RadSFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_RadLFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_DeepSubSurfHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfHumidCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfHeatCapacity(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(in)
| ||
| xyr_Press(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyz_Exner(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in)
| ||
| xyr_Exner(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_SurfVelTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfTempTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfQVapTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) : | real(DP), intent(out)
| ||
| xy_DSurfTempDt(0:imax-1, 1:jmax) : | real(DP), intent(out)
|
時間変化率の計算を行います.
Calculate tendencies.
subroutine PhyImplSDHV3TendencyHeatTQCore( xy_IndexCalcMethod, xy_SeaIceThickness, xy_SnowFrac, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfSoilHeatFlux, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xy_SurfHumidCoef, xy_SurfHeatCapacity, xyra_DelRadLFlux, xyr_Press, xyz_Exner, xyr_Exner, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt )
!
! 時間変化率の計算を行います.
!
! Calculate tendencies.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: Grav, CpDry, GasRDry, LatentHeat
! 雪と海氷の定数の設定
! Setting constants of snow and sea ice
!
use constants_snowseaice, only: SeaIceVolHeatCap , SeaIceThermCondCoef, TempBelowSeaIce
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! 飽和比湿の算出
! Evaluation of saturation specific humidity
!
use saturate, only: xy_CalcQVapSatOnLiq, xy_CalcQVapSatOnSol, xy_CalcDQVapSatDTempOnLiq, xy_CalcDQVapSatDTempOnSol
! 陰解法による時間積分のためのルーチン
! Routines for time integration with implicit scheme
!
use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3
! 宣言文 ; Declaration statements
!
integer , intent(in):: xy_IndexCalcMethod (0:imax-1, 1:jmax)
!
! Index for calculation method
real(DP), intent(in):: xy_SeaIceThickness(0:imax-1, 1:jmax)
!
! Sea ice thickness
real(DP), intent(in):: xy_SnowFrac (0:imax-1, 1:jmax)
!
! Snow fraction
real(DP), intent(in):: xyr_HeatFlux (0:imax-1, 1:jmax, 0:kmax)
! 熱フラックス.
! Heat flux
real(DP), intent(in):: xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax)
! 比湿フラックス.
! Specific humidity flux
real(DP), intent(in):: xy_SurfSoilHeatFlux(0:imax-1, 1:jmax)
! 惑星表面土壌熱伝導フラックス.
! Soil heat conduction flux at the surface
real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
! 短波 (日射) フラックス.
! Shortwave (insolation) flux
real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
! 長波フラックス.
! Longwave flux
real(DP), intent(in):: xy_DeepSubSurfHeatFlux (0:imax-1, 1:jmax)
! 地中熱フラックス.
! "Deep subsurface heat flux"
! Heat flux at the bottom of surface/soil layer.
real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
! 地表面温度.
! Surface temperature
real(DP), intent(in):: xy_SurfHumidCoef (0:imax-1, 1:jmax)
! 地表湿潤度.
! Surface humidity coefficient
real(DP), intent(in):: xy_SurfHeatCapacity (0:imax-1, 1:jmax)
! 地表熱容量.
! Surface heat capacity
real(DP), intent(in):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
! 長波地表温度変化.
! Surface temperature tendency with longwave
real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air pressure (half level)
real(DP), intent(in):: xyz_Exner (0:imax-1, 1:jmax, 1:kmax)
! Exner 関数 (整数レベル).
! Exner function (full level)
real(DP), intent(in):: xyr_Exner (0:imax-1, 1:jmax, 0:kmax)
! Exner 関数 (半整数レベル).
! Exner function (half level)
real(DP), intent(in):: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：運動量.
! Transfer coefficient: velocity
real(DP), intent(in):: xyr_TempTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in):: xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：質量.
! Transfer coefficient: mass of constituents
real(DP), intent(in):: xy_SurfVelTransCoef (0:imax-1, 1:jmax)
! 輸送係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(in):: xy_SurfTempTransCoef (0:imax-1, 1:jmax)
! 輸送係数：温度.
! Transfer coefficient: temperature
real(DP), intent(in):: xy_SurfQVapTransCoef (0:imax-1, 1:jmax)
! 輸送係数：比湿.
! Transfer coefficient: specific humidity
real(DP), intent(out):: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{T}{t} $ . 温度変化.
! Temperature tendency
real(DP), intent(out):: xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ \DP{q}{t} $ . 質量混合比変化.
! Mass mixing ratio tendency
real(DP), intent(out):: xy_DSurfTempDt (0:imax-1, 1:jmax)
! 地表面温度変化率 (K s-1)
! Surface temperature tendency (K s-1)
! 作業変数
! Work variables
!
real(DP):: xyza_TempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
! 温度陰解行列.
! Implicit matrix about temperature
real(DP):: xyz_TempVec(0:imax-1, 1:jmax, 1:kmax)
! 温度陰解ベクトル.
! Implicit vector about temperature
real(DP):: xyza_QMixMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
! 質量混合比陰解行列.
! Implicit matrix about mass mixing ratio
real(DP):: xyz_QMixVec(0:imax-1, 1:jmax, 1:kmax)
! 質量混合比陰解ベクトル.
! Implicit vector about mass mixing ratio
real(DP):: xyaa_SurfMtx(0:imax-1, 1:jmax, 0:0, -1:1)
! 惑星表面エネルギー収支用陰解行列
! Implicit matrix for surface energy balance
real(DP):: xy_SurfRH(0:imax-1,1:jmax)
real(DP):: xy_SurfQVapSatOnLiq(0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfQVapSatOnSol (0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfQVapSat (0:imax-1, 1:jmax)
! 地表飽和比湿.
! Saturated specific humidity on surface
real(DP):: xy_SurfDQVapSatDTempOnLiq (0:imax-1, 1:jmax)
! 地表飽和比湿変化.
! Saturated specific humidity tendency on surface
real(DP):: xy_SurfDQVapSatDTempOnSol (0:imax-1, 1:jmax)
! 地表飽和比湿変化.
! Saturated specific humidity tendency on surface
real(DP):: xy_SurfDQVapSatDTemp (0:imax-1, 1:jmax)
! 地表飽和比湿変化.
! Saturated specific humidity tendency on surface
real(DP):: xyaa_TempQVapLUMtx (0:imax-1, 1:jmax, -kmax:kmax, -1:1)
! LU 行列.
! LU matrix
real(DP):: xya_DelTempQVapLUVec (0:imax-1, 1:jmax, -kmax:kmax)
! $ T, Qv $ の時間変化.
! Tendency of $ T $ and $ Qv $
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: l ! 行列用 DO ループ用作業変数
! Work variables for DO loop of matrices
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. phy_implicit_sdh_V3_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
!!$ ! 計算時間計測開始
!!$ ! Start measurement of computation time
!!$ !
!!$ call TimesetClockStart( module_name )
!!$ if ( .not. FlagSSModel ) then
!!$ call MessageNotify( 'E', module_name, 'FlagSSModel has to be true.' )
!!$ end if
! FlagBucketModel は関係ないよね?
! SSModel 強制にした時点で, 水蒸気は地面と分離したから.
!!$ if ( .not. FlagBucketModel ) then
!!$ call MessageNotify( 'E', module_name, 'FlagBucketModel has to be true.' )
!!$ end if
! 飽和比湿の計算
! Calculate saturated specific humidity
!
xy_SurfQVapSatOnLiq = xy_CalcQVapSatOnLiq ( xy_SurfTemp, xyr_Press(:,:,0) )
xy_SurfQVapSatOnSol = xy_CalcQVapSatOnSol ( xy_SurfTemp, xyr_Press(:,:,0) )
xy_SurfQVapSat = ( 1.0_DP - xy_SnowFrac ) * xy_SurfQVapSatOnLiq + xy_SnowFrac * xy_SurfQVapSatOnSol
xy_SurfDQVapSatDTempOnLiq = xy_CalcDQVapSatDTempOnLiq( xy_SurfTemp, xy_SurfQVapSatOnLiq )
xy_SurfDQVapSatDTempOnSol = xy_CalcDQVapSatDTempOnSol( xy_SurfTemp, xy_SurfQVapSatOnSol )
xy_SurfDQVapSatDTemp = ( 1.0_DP - xy_SnowFrac ) * xy_SurfDQVapSatDTempOnLiq + xy_SnowFrac * xy_SurfDQVapSatDTempOnSol
! 陰解法のための行列作成
! Create matrices for implicit scheme
!
! 鉛直拡散スキームの輸送係数から陰解行列の計算 (温度)
! Calculate implicit matrices from transfer coefficient of vertical diffusion scheme (temperature)
!
k = 1
xyza_TempMtx(:,:,k,-1) = - CpDry * xy_SurfTempTransCoef(:,:)
xyza_TempMtx(:,:,k, 0) = - CpDry * ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k ) * xy_SurfTempTransCoef(:,:) + CpDry * xyr_Exner(:,:,k ) / xyz_Exner(:,:,k ) * xyr_TempTransCoef(:,:,k )
xyza_TempMtx(:,:,k, 1) = - CpDry * xyr_Exner(:,:,k ) / xyz_Exner(:,:,k+1) * xyr_TempTransCoef(:,:,k )
do k = 2, kmax-1
xyza_TempMtx(:,:,k,-1) = - CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k-1) * xyr_TempTransCoef(:,:,k-1)
xyza_TempMtx(:,:,k, 0) = - CpDry * ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k ) * xyr_TempTransCoef(:,:,k-1) + CpDry * xyr_Exner(:,:,k ) / xyz_Exner(:,:,k ) * xyr_TempTransCoef(:,:,k )
xyza_TempMtx(:,:,k, 1) = - CpDry * xyr_Exner(:,:,k ) / xyz_Exner(:,:,k+1) * xyr_TempTransCoef(:,:,k )
end do
k = kmax
xyza_TempMtx(:,:,k,-1) = - CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k-1) * xyr_TempTransCoef(:,:,k-1)
xyza_TempMtx(:,:,k, 0) = - CpDry * ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + CpDry * xyr_Exner(:,:,k-1) / xyz_Exner(:,:,k ) * xyr_TempTransCoef(:,:,k-1)
xyza_TempMtx(:,:,k, 1) = 0.0_DP
do k = 1, kmax
xyz_TempVec(:,:,k) = - ( xyr_HeatFlux(:,:,k) - xyr_HeatFlux(:,:,k-1) )
end do
! 鉛直拡散スキームの輸送係数から陰解行列の計算 (比湿)
! Calculate implicit matrices from transfer coefficient of vertical diffusion scheme (specific humidity)
!
k = 1
xyza_QMixMtx(:,:,k,-1) = - xy_SurfHumidCoef(:,:) * xy_SurfQVapTransCoef(:,:) * xy_SurfDQVapSatDTemp(:,:)
xyza_QMixMtx(:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xy_SurfHumidCoef(:,:) * xy_SurfQVapTransCoef(:,:) + xyr_QMixTransCoef(:,:,k )
xyza_QMixMtx(:,:,k, 1) = - xyr_QMixTransCoef(:,:,k )
do k = 2, kmax-1
xyza_QMixMtx(:,:,k,-1) = - xyr_QMixTransCoef(:,:,k-1)
xyza_QMixMtx(:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xyr_QMixTransCoef(:,:,k-1) + xyr_QMixTransCoef(:,:,k )
xyza_QMixMtx(:,:,k, 1) = - xyr_QMixTransCoef(:,:,k )
end do
k = kmax
xyza_QMixMtx(:,:,k,-1) = - xyr_QMixTransCoef(:,:,k-1)
xyza_QMixMtx(:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xyr_QMixTransCoef(:,:,k-1)
xyza_QMixMtx(:,:,k, 1) = 0.0_DP
n = IndexH2OVap
do k = 1, kmax
xyz_QMixVec(:,:,k) = - ( xyrf_QMixFlux(:,:,k,n) - xyrf_QMixFlux(:,:,k-1,n) )
end do
! 地表面過程の輸送係数から陰解行列の計算
! Calculate implicit matrices from transfer coefficient of surface process
!
do i = 0, imax-1
do j = 1, jmax
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
! for thermal diffusion in soil
!!$ xyaa_SurfMtx(i,j,0,-1) = &
!!$ & xyr_SoilTempTransCoef(i,j,0)
! for water vapor diffusion in atmosphere
xyaa_SurfMtx(i,j,0,-1) = - LatentHeat * xy_SurfHumidCoef(i,j) * xy_SurfQVapTransCoef(i,j)
xyaa_SurfMtx(i,j,0, 0) = xy_SurfHeatCapacity(i,j) / ( 2.0_DP * DelTime ) + CpDry * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,0) + LatentHeat * xy_SurfHumidCoef(i,j) * xy_SurfQVapTransCoef(i,j) * xy_SurfDQVapSatDTemp(i,j)
xyaa_SurfMtx(i,j,0, 1) = - CpDry * xyr_Exner(i,j,0) / xyz_Exner(i,j,1) * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,1)
case ( IndexSeaIce )
! sea ice
!!$ xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
! for water vapor diffusion in atmosphere
xyaa_SurfMtx(i,j,0,-1) = - LatentHeat * xy_SurfHumidCoef(i,j) * xy_SurfQVapTransCoef(i,j)
xyaa_SurfMtx(i,j,0, 0) = SeaIceVolHeatCap * xy_SeaIceThickness(i,j) / ( 2.0_DP * DelTime ) + CpDry * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,0) + SeaIceThermCondCoef / xy_SeaIceThickness(i,j) + LatentHeat * xy_SurfHumidCoef(i,j) * xy_SurfQVapTransCoef(i,j) * xy_SurfDQVapSatDTemp(i,j)
xyaa_SurfMtx(i,j,0, 1) = - CpDry * xyr_Exner(i,j,0) / xyz_Exner(i,j,1) * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,1)
case ( IndexSlabOcean )
! slab ocean ocean
!!$ xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
! for water vapor diffusion in atmosphere
xyaa_SurfMtx(i,j,0,-1) = - LatentHeat * xy_SurfHumidCoef(i,j) * xy_SurfQVapTransCoef(i,j)
xyaa_SurfMtx(i,j,0, 0) = SOHeatCapacity / ( 2.0_DP * DelTime ) + CpDry * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,0) + LatentHeat * xy_SurfHumidCoef(i,j) * xy_SurfQVapTransCoef(i,j) * xy_SurfDQVapSatDTemp(i,j)
xyaa_SurfMtx(i,j,0, 1) = - CpDry * xyr_Exner(i,j,0) / xyz_Exner(i,j,1) * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,1)
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
xyaa_SurfMtx(i,j,0,-1) = 0.0_DP
xyaa_SurfMtx(i,j,0, 0) = 1.0_DP
xyaa_SurfMtx(i,j,0, 1) = 0.0_DP
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
n = IndexH2OVap
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
xy_SurfRH(i,j) = - xyr_RadSFlux(i,j,0) - xyr_RadLFlux(i,j,0) - xyr_HeatFlux(i,j,0) - LatentHeat * xyrf_QMixFlux(i,j,0,n) + xy_SurfSoilHeatFlux(i,j)
case ( IndexSeaIce )
! sea ice
xy_SurfRH(i,j) = - xyr_RadSFlux(i,j,0) - xyr_RadLFlux(i,j,0) - xyr_HeatFlux(i,j,0) - LatentHeat * xyrf_QMixFlux(i,j,0,n) - SeaIceThermCondCoef * ( xy_SurfTemp(i,j) - TempBelowSeaIce ) / xy_SeaIceThickness(i,j)
case ( IndexSlabOcean )
! slab ocean
xy_SurfRH(i,j) = - xyr_RadSFlux(i,j,0) - xyr_RadLFlux(i,j,0) - xyr_HeatFlux(i,j,0) - LatentHeat * xyrf_QMixFlux(i,j,0,n) !&
! & + xy_DeepSubSurfHeatFlux(i,j)
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
xy_SurfRH(i,j) = 0.0_DP
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
! 温度と比湿の計算
! Calculate temperature and specific humidity
!
do l = -1, 1
do k = 1, kmax
xyaa_TempQVapLUMtx(:,:,-k,-l) = xyza_QMixMtx(:,:,k,l)
end do
k = 0
xyaa_TempQVapLUMtx(:,:, k, l) = xyaa_SurfMtx(:,:,0,l)
do k = 1, kmax
xyaa_TempQVapLUMtx(:,:, k, l) = xyza_TempMtx(:,:,k,l)
end do
end do
call PhyImplLUDecomp3( xyaa_TempQVapLUMtx, imax * jmax, kmax + 1 + kmax )
do k = 1, kmax
xya_DelTempQVapLUVec(:,:,-k) = xyz_QMixVec(:,:,k)
end do
k = 0
xya_DelTempQVapLUVec(:,:,k) = xy_SurfRH
do k = 1, kmax
xya_DelTempQVapLUVec(:,:,k) = xyz_TempVec(:,:,k)
end do
call PhyImplLUSolve3( xya_DelTempQVapLUVec, xyaa_TempQVapLUMtx, 1, imax * jmax , kmax + 1 + kmax )
n = IndexH2OVap
do k = 1, kmax
xyzf_DQMixDt(:,:,k,n) = xya_DelTempQVapLUVec(:,:,-k) / ( 2.0_DP * DelTime )
end do
do j = 1, jmax
do i = 0, imax-1
select case ( xy_IndexCalcMethod(i,j) )
case ( IndexLand )
! land
xy_DSurfTempDt(i,j) = xya_DelTempQVapLUVec(i,j,0) / ( 2.0_DP * DelTime )
case ( IndexSeaIce )
! sea ice
xy_DSurfTempDt(i,j) = xya_DelTempQVapLUVec(i,j,0) / ( 2.0_DP * DelTime )
case ( IndexSlabOcean )
! slab ocean
xy_DSurfTempDt(i,j) = xya_DelTempQVapLUVec(i,j,0) / ( 2.0_DP * DelTime )
case ( IndexPresTs, IndexLandWithPresTs )
! prescribed surface temperature
xy_DSurfTempDt(i,j) = 0.
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
end do
do k = 1, kmax
xyz_DTempDt(:,:,k) = xya_DelTempQVapLUVec(:,:,k) / ( 2.0_DP * DelTime )
end do
call PhyImplSDHV3ChkConservationTQ( xy_IndexCalcMethod, xy_SeaIceThickness, xy_SnowFrac, xyr_Press, xyz_Exner, xyr_Exner, xy_SurfTemp, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfSoilHeatFlux, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xy_SurfHumidCoef, xy_SurfHeatCapacity, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux, xy_DeepSubSurfHeatFlux, xyz_DTempDt, xyzf_DQMixDt, xy_DSurfTempDt )
!!$ ! 計算時間計測一時停止
!!$ ! Pause measurement of computation time
!!$ !
!!$ call TimesetClockStop( module_name )
end subroutine PhyImplSDHV3TendencyHeatTQCore
| Subroutine : | |||
| xyr_MomFluxX(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_MomFluxY(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_Press(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xyr_VelTransCoef(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in)
| ||
| xy_SurfVelTransCoef(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xyz_DUDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
| ||
| xyz_DVDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(out)
|
時間変化率の計算を行います.
Calculate tendencies.
subroutine PhyImplSDHV3TendencyMomCore( xyr_MomFluxX, xyr_MomFluxY, xyr_Press, xyr_VelTransCoef, xy_SurfVelTransCoef, xyz_DUDt, xyz_DVDt )
!
! 時間変化率の計算を行います.
!
! Calculate tendencies.
!
! モジュール引用 ; USE statements
!
! 物理定数設定
! Physical constants settings
!
use constants, only: Grav, CpDry, GasRDry
! $ R $ [J kg-1 K-1].
! 乾燥大気の気体定数.
! Gas constant of air
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! 陰解法による時間積分のためのルーチン
! Routines for time integration with implicit scheme
!
use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3
! 宣言文 ; Declaration statements
!
real(DP), intent(in):: xyr_MomFluxX (0:imax-1, 1:jmax, 0:kmax)
! 東西方向運動量フラックス.
! Eastward momentum flux
real(DP), intent(in):: xyr_MomFluxY (0:imax-1, 1:jmax, 0:kmax)
! 南北方向運動量フラックス.
! Northward momentum flux
real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
! $ \hat{p} $ . 気圧 (半整数レベル).
! Air pressure (half level)
real(DP), intent(in):: xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax)
! 輸送係数：運動量.
! Transfer coefficient: velocity
real(DP), intent(in):: xy_SurfVelTransCoef (0:imax-1, 1:jmax)
! 輸送係数：運動量.
! Diffusion coefficient: velocity
real(DP), intent(out):: xyz_DUDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{u}{t} $ . 東西風速変化.
! Eastward wind tendency
real(DP), intent(out):: xyz_DVDt (0:imax-1, 1:jmax, 1:kmax)
! $ \DP{v}{t} $ . 南北風速変化.
! Northward wind tendency
! 作業変数
! Work variables
!
real(DP):: xyza_UVMtx (0:imax-1, 1:jmax, 1:kmax, -1:1)
! 速度陰解行列.
! Implicit matrix about velocity
real(DP):: xyz_UVec (0:imax-1, 1:jmax, 1:kmax)
! 速度陰解ベクトル.
! Implicit vector about velocity
real(DP):: xyz_VVec (0:imax-1, 1:jmax, 1:kmax)
! 速度陰解ベクトル.
! Implicit vector about velocity
real(DP):: xyza_UVLUMtx (0:imax-1, 1:jmax, 1:kmax,-1:1)
! LU 行列.
! LU matrix
integer:: i ! 経度方向に回る DO ループ用作業変数
! Work variables for DO loop in longitude
integer:: j ! 緯度方向に回る DO ループ用作業変数
! Work variables for DO loop in latitude
integer:: k ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: l ! 行列用 DO ループ用作業変数
! Work variables for DO loop of matrices
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. phy_implicit_sdh_V3_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
!!$ ! 計算時間計測開始
!!$ ! Start measurement of computation time
!!$ !
!!$ call TimesetClockStart( module_name )
!!$ if ( .not. FlagSSModel ) then
!!$ call MessageNotify( 'E', module_name, 'FlagSSModel has to be true.' )
!!$ end if
! FlagBucketModel は関係ないよね?
! SSModel 強制にした時点で, 水蒸気は地面と分離したから.
!!$ if ( .not. FlagBucketModel ) then
!!$ call MessageNotify( 'E', module_name, 'FlagBucketModel has to be true.' )
!!$ end if
! 陰解法のための行列作成
! Create matrices for implicit scheme
!
! 鉛直拡散スキームの輸送係数から陰解行列の計算 (速度)
! Calculate implicit matrices from transfer coefficient of vertical diffusion scheme (velocity)
!
k = 1
xyza_UVMtx (:,:,k,-1) = 0.0_DP
xyza_UVMtx (:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xy_SurfVelTransCoef(:,:) + xyr_VelTransCoef(:,:,k )
xyza_UVMtx (:,:,k, 1) = - xyr_VelTransCoef(:,:,k)
do k = 2, kmax-1
xyza_UVMtx (:,:,k,-1) = - xyr_VelTransCoef(:,:,k-1)
xyza_UVMtx (:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xyr_VelTransCoef(:,:,k-1) + xyr_VelTransCoef(:,:,k )
xyza_UVMtx (:,:,k, 1) = - xyr_VelTransCoef(:,:,k)
end do
k = kmax
xyza_UVMtx (:,:,k,-1) = - xyr_VelTransCoef(:,:,k-1)
xyza_UVMtx (:,:,k, 0) = - ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) ) / Grav / ( 2.0_DP * DelTime ) + xyr_VelTransCoef(:,:,k-1)
xyza_UVMtx (:,:,k, 1) = 0.0_DP
do k = 1, kmax
xyz_UVec(:,:,k) = - ( xyr_MomFluxX(:,:,k) - xyr_MomFluxX(:,:,k-1) )
xyz_VVec(:,:,k) = - ( xyr_MomFluxY(:,:,k) - xyr_MomFluxY(:,:,k-1) )
end do
! 東西風速, 南北風速の計算
! Calculate eastward and northward wind
!
xyza_UVLUMtx = xyza_UVMtx
call PhyImplLUDecomp3( xyza_UVLUMtx, imax * jmax, kmax ) ! (in)
do k = 1, kmax
xyz_DUDt(:,:,k) = xyz_UVec(:,:,k)
xyz_DVDt(:,:,k) = xyz_VVec(:,:,k)
end do
call PhyImplLUSolve3( xyz_DUDt, xyza_UVLUMtx, 1, imax * jmax, kmax ) ! (in)
call PhyImplLUSolve3( xyz_DVDt, xyza_UVLUMtx, 1, imax * jmax, kmax ) ! (in)
do k = 1, kmax
xyz_DUDt(:,:,k) = xyz_DUDt(:,:,k) / ( 2. * DelTime )
xyz_DVDt(:,:,k) = xyz_DVDt(:,:,k) / ( 2. * DelTime )
end do
!!$ ! 計算時間計測一時停止
!!$ ! Pause measurement of computation time
!!$ !
!!$ call TimesetClockStop( module_name )
end subroutine PhyImplSDHV3TendencyMomCore
| Constant : | |||
| module_name = ‘phy_implicit_sdh_V3‘ : | character(*), parameter
|
| Variable : | |||
| phy_implicit_sdh_V3_inited = .false. : | logical, save
|
| Constant : | |||
| version = ’$Name: $’ // ’$Id: phy_implicit_sdh_V3.f90,v 1.1 2015/01/29 12:05:01 yot Exp $’ : | character(*), parameter
|