| Class | phy_implicit_sdh |
| In: |
phy_implicit/phy_implicit_sdh.f90
|
Note that Japanese and English are described in parallel.
| PhyImplSDHTendency : | 時間変化率の計算 |
| PhyImplSDHSetCalcMethodFromMatthewsIndex : | SurfType から計算法インデクスの作成 |
| PhyImplSDHInit : | 初期化 |
| —————————————- : | ———— |
| PhyImplSDHTendency : | Calculate tendency |
| PhyImplSDHSetCalcMethodFromMatthewsIndex : | Set index for calculation method |
| PhyImplSDHInit : | Initialization |
| Subroutine : | |||
| ArgFlagSlabOcean : | logical, intent(in )
| ||
| ArgFlagMajCompPhaseChange : | logical, intent(in )
|
This procedure input/output NAMELIST#phy_implicit_sdh_nml .
subroutine PhyImplSDHInit( ArgFlagSlabOcean, ArgFlagMajCompPhaseChange )
!
! 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
! 宣言文 ; Declaration statements
!
logical, intent(in ) :: ArgFlagSlabOcean
! flag for use of slab ocean
logical, intent(in ) :: ArgFlagMajCompPhaseChange
! flag for use of major component phase change
! 作業変数
! 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_nml/ SOHeatCapacity ! Slab ocean heat capacity (J m-2 K-1)
!
! デフォルト値については初期化手続 "phy_implicit#PhyImplInit"
! のソースコードを参照のこと.
!
! Refer to source codes in the initialization procedure
! "phy_implicit#PhyImplInit" for the default values.
!
! 実行文 ; Executable statement
!
if ( phy_implicit_sdh_inited ) return
! Set flag for slab ocean
FlagSlabOcean = ArgFlagSlabOcean
! Set flag for major component phase change
FlagMajCompPhaseChange = ArgFlagMajCompPhaseChange
! デフォルト値の設定
! Default values settings
!
SOHeatCapacity = 4.187d3 * 1.0d3 * 60.0_DP
! 4.187d3 (J (kg K)-1) * 1.0d3 (kg m-3) * 60.0d0 (m)
! 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_nml, iostat = iostat_nml ) ! (out)
close( unit_nml )
call NmlutilMsg( iostat_nml, module_name ) ! (in)
end if
! 印字 ; Print
!
call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
call MessageNotify( 'M', module_name, ' SOHeatCapacity = %f', d = (/ SOHeatCapacity /) )
call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )
phy_implicit_sdh_inited = .true.
end subroutine PhyImplSDHInit
| Subroutine : | |||
| xy_SurfCond(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 PhyImplSDHSetCalcMethodFromMatthewsIndex( xy_SurfCond, xy_SeaIceConc, xy_IndexCalcMethod )
!
!
!
! Set index for calculation method from Matthews' index
!
! モジュール引用 ; USE statements
!
! 雪と海氷の定数の設定
! Setting constants of snow and sea ice
!
use constants_snowseaice, only: SeaIceThreshold
! 宣言文 ; Declaration statements
!
integer , intent(in ) :: xy_SurfCond (0:imax-1, 1:jmax)
! 地表状態.
! Surface condition
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_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) >= 1 ) then
! land
xy_IndexCalcMethod(i,j) = IndexLand
else
if ( xy_SeaIceConc(i,j) > SeaIceThreshold ) then
! sea ice
xy_IndexCalcMethod(i,j) = IndexSeaIce
else if ( FlagSlabOcean ) then
! slab ocean
xy_IndexCalcMethod(i,j) = IndexSlabOcean
else
! open ocean
xy_IndexCalcMethod(i,j) = IndexOceanPresSST
end if
end if
end do
end do
end subroutine PhyImplSDHSetCalcMethodFromMatthewsIndex
| Subroutine : | |||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , 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_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_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 PhyImplSDHTendency( xy_IndexCalcMethod, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, 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_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
!
! 座標データ設定
! 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, SeaIceThreshold, SeaIceThickness, TempBelowSeaIce
! 時刻管理
! Time control
!
use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop
! バケツモデル
! bucket model
!
use Bucket_Model, only: FlagBucketModel, FlagBucketModelSnow
! 陰解法による時間積分のためのルーチン
! 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):: 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(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_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):: 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):: xy_SurfQVapSat (0:imax-1, 1:jmax)
!!$ ! 地表飽和比湿.
!!$ ! Saturated specific humidity on surface
!!$ real(DP):: xy_SurfDQVapSatDTemp (0:imax-1, 1:jmax)
!!$ ! 地表飽和比湿変化.
!!$ ! Saturated specific humidity tendency on surface
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)
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_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. * 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. * 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. * 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 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. * 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. * 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. * 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. * 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. * 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. * DelTime ) + xyr_QMixTransCoef(:,:,k-1)
xyza_QMixMtx(:,:,k, 1) = 0.0_DP
do n = 1, ncmax
do k = 1, kmax
xyzf_QMixVec(:,:,k,n) = - ( xyrf_QMixFlux(:,:,k,n) - xyrf_QMixFlux(:,:,k-1,n) )
end do
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. * 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. * 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. * 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 * SeaIceThickness / ( 2.0d0 * DelTime ) + CpDry * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,0) + SeaIceThermCondCoef / SeaIceThickness
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. * 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 ( IndexOceanPresSST )
! open ocean
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 ) / SeaIceThickness
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 ( IndexOceanPresSST )
! open ocean
xy_SurfRH(i,j) = 0.0_DP
case default
call MessageNotify( 'E', module_name, 'Unexpected Error.' )
end select
end do
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
! 温度と比湿の計算
! 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. * 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. * DelTime )
case ( IndexSeaIce )
! sea ice
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexSlabOcean )
! slab ocean
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexOceanPresSST )
! open ocean
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_DelTempSoilTempLUVec(:,:,k) / ( 2. * DelTime )
end do
!
! Calculation of tendencies of soil moisture and surface snow amount
!
if ( FlagBucketModel ) then
if ( FlagBucketModelSnow ) 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
xy_DSurfSnowDt(i,j) = - xyrf_QMixFlux(i,j,0,IndexH2OVap)
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.0d0 * DelTime )
xy_DSurfSnowDt (i,j) = - xy_SurfSnowB(i,j) / ( 2.0d0 * 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_DSoilMoistDt(i,j) = 0.0_DP
else
xy_DSurfSnowDt (i,j) = 0.0_DP
xy_DSoilMoistDt(i,j) = - xyrf_QMixFlux(i,j,0,IndexH2OVap)
end if
end if
end do
end do
else
! Evaporation is subtracted from soil moisture
!
xy_DSoilMoistDt = - xyrf_QMixFlux(:,:,0,IndexH2OVap)
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
!!$ call PhyImplSDHSnowMeltCorrection( &
!!$ & xyr_Press(:,:,0), & ! (in)
!!$ & xyr_HeatFlux, xy_SurfLatentHeatFlux, & ! (in)
!!$ & xyr_SoilHeatFlux, & ! (in)
!!$ & xyr_SoilTempTransCoef, & ! (in)
!!$ & xyr_RadSFlux, xyr_RadLFlux, & ! (in)
!!$ & xy_DeepSubSurfHeatFlux, & ! (in)
!!$ & xy_SurfTemp, xyz_SoilTemp, & ! (in)
!!$ & xy_SurfSnowB, & ! (in)
!!$ & xy_SurfMajCompIceB, & ! (in)
!!$ & xy_SurfHeatCapacity, & ! (in)
!!$ & xy_SoilHeatCap, xy_SoilHeatDiffCoef, & ! (in)
!!$ & xy_IndexCalcMethod, & ! (in)
!!$ & xyra_DelRadLFlux, & ! (in)
!!$ & xyz_Exner, xyr_Exner, & ! (in)
!!$ & xy_SurfTempTransCoef, & ! (in)
!!$ & xyza_TempMtx, xyz_TempVec, & ! (in)
!!$ & xyaa_SurfMtx, xy_SurfRH, & ! (in)
!!$ & xyaa_SoilTempMtx, xya_SoilTempVec, & ! (in)
!!$ & xyz_DTempDt, & ! (in)
!!$ & xy_DSurfTempDt, & ! (inout)
!!$ & xyz_DSoilTempDt, & ! (inout)
!!$ & xy_DSurfMajCompIceDt, & ! (inout)
!!$ & xy_DSoilMoistDt, xy_DSurfSnowDt, & ! (inout)
!!$ & xy_LatHeatFluxByMajCompIceSubl, & ! (out)
!!$ & xy_LatHeatFluxBySnowMelt & ! (out)
!!$ & )
if ( FlagMajCompPhaseChange ) then
xy_DAtmMassDt = 0.0_DP
xy_DSurfMajCompIceDt = 0.0_DP
! A dummy value
!
xy_LatHeatFluxByOtherSpc = 0.0_DP
call PhyImplSDHIceSnowPhaseChangeCorrection( IndexSpcCO2, xyr_Press(:,:,0), xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, 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
call PhyImplSDHIceSnowPhaseChangeCorrection( IndexSpcH2O, xyr_Press(:,:,0), xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, 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 )
call PhyImplSDHSeaIceCorrection( xy_IndexCalcMethod, 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 )
!!$ call PhyImplSDHSeaIceCorrection_OLD_Difficult( &
!!$ & xy_IndexCalcMethod, & ! (in)
!!$ & xy_SurfTemp, & ! (in)
!!$ & xyr_HeatFlux, & ! (in)
!!$ & xyr_SoilHeatFlux, & ! (in)
!!$ & xy_SurfTempTransCoef, & ! (in)
!!$ & xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux, & ! (in)
!!$ & xy_SurfLatentHeatFlux, & ! (in)
!!$ & xyza_TempMtx, xyz_TempVec, & ! (in)
!!$ & xyaa_SurfMtx, xy_SurfRH, & ! (in)
!!$ & xyaa_SoilTempMtx, xya_SoilTempVec, & ! (in)
!!$ & xyz_DTempDt, & ! (inout)
!!$ & xy_DSurfTempDt, & ! (inout)
!!$ & xyz_DSoilTempDt, & ! (inout)
!!$ & xy_LatHeatFluxBySeaIceMelt & ! (out)
!!$ & )
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. * DelTime )
end do
end do
! Debug routine
!
call PhyImplSDHChkConservation( xy_IndexCalcMethod, xyr_Press, xyz_Exner, xyr_Exner, xy_SurfTemp, xy_SurfHeatCapacity, xy_SoilHeatCap, xyr_HeatFlux, xyrf_QMixFlux, 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 PhyImplSDHTendency
| Variable : | |||
| FlagMajCompPhaseChange : | logical, save
|
| Subroutine : | |||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , 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)
| ||
| xyrf_QMixFlux(0:imax-1, 1:jmax, 0:kmax, 1:ncmax) : | 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 PhyImplSDHChkConservation( xy_IndexCalcMethod, xyr_Press, xyz_Exner, xyr_Exner, xy_SurfTemp, xy_SurfHeatCapacity, xy_SoilHeatCap, xyr_HeatFlux, xyrf_QMixFlux, 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, SeaIceThreshold, SeaIceThickness, TempBelowSeaIce
! 宣言文 ; Declaration statements
!
integer , intent(in):: xy_IndexCalcMethod (0:imax-1, 1:jmax)
!
! Index for calculation method
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_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_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.0d0 * DelTime ) + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0d0 * 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_SeaIceHeatCondFlux = - SeaIceThermCondCoef * ( xy_SurfTemp + xy_DSurfTempDt * ( 2.0_DP * DelTime ) - TempBelowSeaIce ) / 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_Residual = xy_SurfRadSFlux + xy_SurfRadLFlux + xy_SurfSensHeatFlux + xy_SurfLatentHeatFlux - xy_SurfSoilHeatCondFlux + xy_LatHeatFluxByMajCompIceSubl + xy_LatHeatFluxBySnowMelt
!
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, '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 ( 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, 'Slab ocean heating res. : %f.', d = (/ MaxResidual /) )
end if
!-----
! Sea ice heating
!
xy_SumAtmRate = SeaIceVolHeatCap * 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 ( 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, '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 - xyrf_QMixFlux(:,:,0,IndexH2OVap) ) - 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.0d-10 ) 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 ( 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, '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 PhyImplSDHChkConservation
| 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_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 PhyImplSDHIceSnowPhaseChangeCorrection( IndexSpc, xy_Ps, xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, 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
! バケツモデル
! bucket model
!
use Bucket_Model, only: FlagBucketModel, FlagBucketModelSnow
! CO2 相変化
! Phase change of CO2
!
use co2_phase_change, only : LatentHeatCO2Subl, CO2Saturate2D
! 宣言文 ; 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_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)
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):: SurfSolATentative
real(DP):: xy_SurfSolATentativeSave(0:imax-1, 1:jmax)
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)
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_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. FlagBucketModelSnow ) ) return
if ( kslmax == 0 ) then
call MessageNotify( 'E', module_name, 'kslmax <= 0 in PhyImplSDHSnowMeltCorrection.' )
else
xy_DSurfLiqDtSave = xy_DSurfLiqDt
xy_DSurfSolDtSave = xy_DSurfSolDt
select case ( IndexSpc )
case ( IndexSpcH2O )
xy_TempCond = TempCondWater
LatentHeatLocal = LatentHeatFusion
case ( IndexSpcCO2 )
call CO2Saturate2D( xy_Ps, xy_TempCond )
LatentHeatLocal = LatentHeatCO2Subl
case default
call MessageNotify( 'E', module_name, 'Undefined IndexSpc, %d.', i = (/ IndexSpc /) )
end select
xy_SurfSolATentativeSave = xy_SurfSolB + xy_DSurfSolDt * ( 2.0_DP * DelTime )
!----------
! A case that a part of snow/ice melts
!----------
do j = 1, jmax
do i = 0, imax-1
if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
SurfTempATentative = xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * 2.0d0 * DelTime
!!$ SurfSnowATentative = xy_SurfSnowB(i,j) &
!!$ & + xy_DSurfSnowDt(i,j) * 2.0d0 * 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
else
xy_FlagCalc(i,j) = .false.
end if
end do
end do
select case ( IndexSpc )
case ( IndexSpcCO2 )
do j = 1, jmax
do i = 0, imax-1
if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
SurfTempATentative = xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * 2.0d0 * DelTime
if ( SurfTempATentative < xy_TempCond(i,j) ) then
xy_FlagCalc(i,j) = .true.
end if
end if
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. * 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. * DelTime )
case ( IndexSeaIce )
! sea ice
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexSlabOcean )
! slab ocean
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexOceanPresSST )
! open ocean
xy_DSurfTempDt(i,j) = 0.
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. * 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.0d0 * DelTime ) + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0d0 * 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_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 melts
!----------
do j = 1, jmax
do i = 0, imax-1
if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
if ( xy_FlagCalc(i,j) ) then
SurfSolATentative = xy_SurfSolB(i,j) + xy_DSurfSolDt(i,j) * 2.0d0 * DelTime
if ( SurfSolATentative < 0.0_DP ) then
xy_FlagCalc(i,j) = .true.
else
xy_FlagCalc(i,j) = .false.
end if
else
xy_FlagCalc(i,j) = .false.
end if
else
xy_FlagCalc(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_FlagCalc(i,j) ) then
!!$ SurfSnowATentative = xy_SurfSnowB(i,j) &
!!$ & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
SurfSolATentative = xy_SurfSolATentativeSave(i,j)
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 * SurfSolATentative / ( 2.0d0 * 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. * 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. * DelTime )
case ( IndexSeaIce )
! sea ice
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexSlabOcean )
! slab ocean
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexOceanPresSST )
! open ocean
xy_DSurfTempDt(i,j) = 0.
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. * 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
SurfSolATentative = xy_SurfSolATentativeSave(i,j)
xy_LatHeatFluxBySnowMelt(i,j) = LatentHeatLocal * SurfSolATentative / ( 2.0d0 * 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
!----------
end if
!!$ ! 計算時間計測一時停止
!!$ ! Pause measurement of computation time
!!$ !
!!$ call TimesetClockStop( module_name )
end subroutine PhyImplSDHIceSnowPhaseChangeCorrection
| 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_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 PhyImplSDHIceSnowPhaseChangeCorrection_OLD_Difficult( IndexSpc, xy_Ps, xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, 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
! バケツモデル
! bucket model
!
use Bucket_Model, only: FlagBucketModel, FlagBucketModelSnow
! CO2 相変化
! Phase change of CO2
!
use co2_phase_change, only : LatentHeatCO2Subl, CO2Saturate2D
! 宣言文 ; 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_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)
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):: SurfSolATentative
real(DP):: xy_SurfSolATentativeSave(0:imax-1, 1:jmax)
real(DP) :: xy_TempMajCompCond(0:imax-1, 1:jmax)
real(DP) :: SurfMajCompIceATentative
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_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. FlagBucketModelSnow ) ) return
if ( kslmax == 0 ) then
call MessageNotify( 'E', module_name, 'kslmax <= 0 in PhyImplSDHSnowMeltCorrection.' )
else
xy_DSurfLiqDtSave = xy_DSurfLiqDt
xy_DSurfSolDtSave = xy_DSurfSolDt
select case ( IndexSpc )
case ( IndexSpcH2O )
xy_TempCond = TempCondWater
LatentHeatLocal = LatentHeatFusion
case ( IndexSpcCO2 )
call CO2Saturate2D( xy_Ps, xy_TempCond )
LatentHeatLocal = LatentHeatCO2Subl
case default
call MessageNotify( 'E', module_name, 'Undefined IndexSpc, %d.', i = (/ IndexSpc /) )
end select
xy_SurfSolATentativeSave = xy_SurfSolB + xy_DSurfSolDt * ( 2.0_DP * DelTime )
!----------
! A case that a part of snow melts
!----------
do j = 1, jmax
do i = 0, imax-1
if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
SurfTempATentative = xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * 2.0d0 * DelTime
!!$ SurfSnowATentative = xy_SurfSnowB(i,j) &
!!$ & + xy_DSurfSnowDt(i,j) * 2.0d0 * 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
else
xy_FlagCalc(i,j) = .false.
end if
end do
end do
select case ( IndexSpc )
case ( IndexSpcCO2 )
do j = 1, jmax
do i = 0, imax-1
if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
SurfTempATentative = xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * 2.0d0 * DelTime
if ( SurfTempATentative < xy_TempCond(i,j) ) then
xy_FlagCalc(i,j) = .true.
end if
end if
end do
end do
end select
xyza_TempMtx = xyza_ArgTempMtx
xyz_TempVec = xyz_ArgTempVec
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
k = 1
xyza_TempMtx(i,j,k,-1) = 0.0_DP
xyz_TempVec(i,j,k) = - ( xyr_HeatFlux(i,j,k) - xyr_HeatFlux(i,j,k-1) )
xyz_TempVec(i,j,k) = xyz_TempVec(i,j,k) + CpDry * xy_SurfTempTransCoef(i,j) * ( xy_TempCond(i,j) - xy_SurfTemp(i,j) )
end if
end do
end do
!
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, 0) = 1.0_DP
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_SurfHeatCapacity(i,j) / ( 2. * DelTime ) + CpDry * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,0) - xyr_SoilTempTransCoef(i,j,0) ) * ( xy_TempCond(i,j) - xy_SurfTemp(i,j) ) - xy_LatHeatFluxByOtherSpc(i,j)
end if
end do
end do
!
xyaa_SoilTempMtx = xyaa_ArgSoilTempMtx
xya_SoilTempVec = xya_ArgSoilTempVec
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
k = 1
xyaa_SoilTempMtx(i,j,k,-1) = 0.0_DP
xya_SoilTempVec(i,j,k) = - ( xyr_SoilHeatFlux(i,j,k) - xyr_SoilHeatFlux(i,j,k-1) )
xya_SoilTempVec(i,j,k) = xya_SoilTempVec(i,j,k) + xyr_SoilTempTransCoef(i,j,k-1) * ( xy_TempCond(i,j) - xy_SurfTemp(i,j) )
end if
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
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. * 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
xy_DSurfTempDt(i,j) = ( xy_TempCond(i,j) - xy_SurfTemp(i,j) ) / ( 2.0_DP * DelTime )
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. * DelTime )
end if
end do
end do
end do
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
xy_LatHeatFluxBySnowMelt(i,j) = xya_DelTempSoilTempLUVec(i,j,0)
!!$ xy_DSurfSnowDt(i,j) = - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
!!$ xy_DSoilMoistDt(i,j) = - xy_DSurfSnowDt(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 melts
!----------
do j = 1, jmax
do i = 0, imax-1
if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
if ( xy_FlagCalc(i,j) ) then
SurfSolATentative = xy_SurfSolB(i,j) + xy_DSurfSolDt(i,j) * 2.0d0 * DelTime
!!$ SurfSnowATentative = xy_SurfSnowATentativeSave(i,j)
!!$ if ( SurfSnowATentative < 0.0_DP ) then
!!$ SurfTempATentative = &
!!$ & xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * ( 2.0_DP * DelTime )
!!$ if ( SurfTempATentative < xy_TempCond(i,j) ) then
if ( SurfSolATentative < 0.0_DP ) then
xy_FlagCalc(i,j) = .true.
else
xy_FlagCalc(i,j) = .false.
end if
else
xy_FlagCalc(i,j) = .false.
end if
else
xy_FlagCalc(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_FlagCalc(i,j) ) then
!!$ SurfSnowATentative = xy_SurfSnowB(i,j) &
!!$ & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
SurfSolATentative = xy_SurfSolATentativeSave(i,j)
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 * SurfSolATentative / ( 2.0d0 * 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. * 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. * DelTime )
case ( IndexSeaIce )
! sea ice
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexSlabOcean )
! slab ocean
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexOceanPresSST )
! open ocean
xy_DSurfTempDt(i,j) = 0.
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. * 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
SurfSolATentative = xy_SurfSolATentativeSave(i,j)
xy_LatHeatFluxBySnowMelt(i,j) = LatentHeatLocal * SurfSolATentative / ( 2.0d0 * 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
!----------
end if
!!$ ! 計算時間計測一時停止
!!$ ! Pause measurement of computation time
!!$ !
!!$ call TimesetClockStop( module_name )
end subroutine PhyImplSDHIceSnowPhaseChangeCorrection_OLD_Difficult
| Subroutine : | |||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , 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 PhyImplSDHSeaIceCorrection( xy_IndexCalcMethod, 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, SeaIceThreshold , SeaIceVolHeatCap , SeaIceThickness , TempBelowSeaIce
! 陰解法による時間積分のためのルーチン
! Routines for time integration with implicit scheme
!
use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3
! バケツモデル
! bucket model
!
use Bucket_Model, only: FlagBucketModel, FlagBucketModelSnow
! 宣言文 ; Declaration statements
!
integer , intent(in ) :: xy_IndexCalcMethod(0:imax-1, 1:jmax)
!
! Index for calculation method
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_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. FlagBucketModelSnow ) ) return
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. * 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. * DelTime )
case ( IndexSeaIce )
! sea ice
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexSlabOcean )
! slab ocean
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexOceanPresSST )
! open ocean
xy_DSurfTempDt(i,j) = 0.
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. * 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.0d0 * DelTime ) + xyra_DelRadLFlux(:,:,0,1) * xyz_DTempDt(:,:,1) * ( 2.0d0 * 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_SeaIceHeatCondFlux = - SeaIceThermCondCoef * ( xy_SurfTemp + xy_DSurfTempDt * ( 2.0_DP * DelTime ) - TempBelowSeaIce ) / SeaIceThickness
xy_HeatingTendency = SeaIceVolHeatCap * 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 PhyImplSDHSeaIceCorrection
| Subroutine : | |||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , intent(in )
| ||
| xy_SurfTemp(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_SurfSensHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xy_SurfLatentHeatFlux(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xy_SeaIceHeatCondFlux(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in )
| ||
| xy_DSurfTempDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xy_LatHeatFluxBySeaIceMelt(0:imax-1, 1:jmax) : | real(DP), intent(out ) |
融雪による時間変化率の修正を行います.
Correction of tendencies due to melt of snow.
subroutine PhyImplSDHSeaIceCorrection_OLD( xy_IndexCalcMethod, xy_SurfTemp, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux, xy_SurfSensHeatFlux, xy_SurfLatentHeatFlux, xy_SeaIceHeatCondFlux, xyz_DTempDt, xy_DSurfTempDt, 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
! 雪と海氷の定数の設定
! Setting constants of snow and sea ice
!
use constants_snowseaice, only: TempCondWater , SeaIceThermCondCoef, SeaIceThreshold , SeaIceVolHeatCap , SeaIceThickness
! バケツモデル
! bucket model
!
use Bucket_Model, only: FlagBucketModel, FlagBucketModelSnow
! 宣言文 ; Declaration statements
!
integer , intent(in ) :: xy_IndexCalcMethod(0:imax-1, 1:jmax)
!
! Index for calculation method
real(DP), intent(in ) :: xy_SurfTemp (0:imax-1, 1:jmax)
! 地表面温度.
! Surface 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_SurfSensHeatFlux(0:imax-1, 1:jmax)
!
! Sensible heat flux used only for debug
real(DP), intent(in ):: xy_SurfLatentHeatFlux(0:imax-1, 1:jmax)
! 惑星表面潜熱フラックス.
! Latent heat flux at the surface
real(DP), intent(in ):: xy_SeaIceHeatCondFlux (0:imax-1, 1:jmax)
!
! Heat conduction flux at the bottom of sea ice
real(DP), intent(in ):: 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(out ) :: xy_LatHeatFluxBySeaIceMelt(0:imax-1, 1:jmax)
! 作業変数
! Work variables
!
real(DP) :: SurfTempATentative
real(DP) :: SurfRadLFlux
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_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. FlagBucketModelSnow ) ) return
do j = 1, jmax
do i = 0, imax-1
! old version
!!$ if ( &
!!$ & ( xy_IndexCalcMethod(i,j) == IndexSeaIce ) .and. &
!!$ & ( xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * DelTime > TempCondWater ) &
!!$ & ) then
!!$
!!$ xy_DSurfTempDt(i,j) = ( TempCondWater - xy_SurfTemp(i,j) ) / DelTime
!!$
!!$ end if
SurfTempATentative = xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * 2.0_DP * DelTime
if ( ( xy_IndexCalcMethod(i,j) == IndexSeaIce ) .and. ( SurfTempATentative > TempCondWater ) ) then
xy_DSurfTempDt(i,j) = ( TempCondWater - xy_SurfTemp(i,j) ) / ( 2.0_DP * DelTime )
SurfRadLFlux = xyr_RadLFlux(i,j,0) + xyra_DelRadLFlux(i,j,0,0) * xy_DSurfTempDt(i,j) * ( 2.0d0 * DelTime ) + xyra_DelRadLFlux(i,j,0,1) * xyz_DTempDt(i,j,1) * ( 2.0d0 * DelTime )
xy_LatHeatFluxBySeaIceMelt(i,j) = - SeaIceVolHeatCap * SeaIceThickness * xy_DSurfTempDt(i,j) - xyr_RadSFlux(i,j,0) - SurfRadLFlux - xy_SurfSensHeatFlux(i,j) - xy_SurfLatentHeatFlux(i,j) + xy_SeaIceHeatCondFlux(i,j)
if ( xy_LatHeatFluxBySeaIceMelt(i,j) < 0.0_DP ) then
call MessageNotify( 'E', module_name, 'Latent heat flux by sea ice melt is negative, %f.', d = (/ xy_LatHeatFluxBySeaIceMelt(i,j) /) )
end if
!!$ xy_LatHeatFluxBySeaIceMelt(i,j) = &
!!$ & ( SurfTempATentative - TempCondWater ) &
!!$ & * SeaIceVolHeatCap * SeaIceThickness &
!!$ & / ( 2.0_DP * DelTime )
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 PhyImplSDHSeaIceCorrection_OLD
| Subroutine : | |||
| xy_IndexCalcMethod(0:imax-1, 1:jmax) : | integer , 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_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 PhyImplSDHSeaIceCorrection_OLD_Difficult( xy_IndexCalcMethod, xy_SurfTemp, xyr_HeatFlux, xyr_SoilHeatFlux, xy_SurfTempTransCoef, 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, SeaIceThreshold , SeaIceVolHeatCap , SeaIceThickness , TempBelowSeaIce
! 陰解法による時間積分のためのルーチン
! Routines for time integration with implicit scheme
!
use phy_implicit_utils, only : PhyImplLUDecomp3, PhyImplLUSolve3
! バケツモデル
! bucket model
!
use Bucket_Model, only: FlagBucketModel, FlagBucketModelSnow
! 宣言文 ; Declaration statements
!
integer , intent(in ) :: xy_IndexCalcMethod(0:imax-1, 1:jmax)
!
! Index for calculation method
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_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
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_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. FlagBucketModelSnow ) ) return
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
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagSeaIceMelt(i,j) ) then
k = 1
xyza_TempMtx(i,j,k,-1) = 0.0_DP
xyz_TempVec(i,j,k) = - ( xyr_HeatFlux(i,j,k) - xyr_HeatFlux(i,j,k-1) )
xyz_TempVec(i,j,k) = xyz_TempVec(i,j,k) + CpDry * xy_SurfTempTransCoef(i,j) * ( xy_TempCond(i,j) - xy_SurfTemp(i,j) )
end if
end do
end do
!
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, 0) = 1.0_DP
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 ) / SeaIceThickness - ( SeaIceVolHeatCap * SeaIceThickness / ( 2.0d0 * DelTime ) + CpDry * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,0) + SeaIceThermCondCoef / SeaIceThickness ) * ( 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. * 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
xy_DSurfTempDt(i,j) = ( xy_TempCond(i,j) - xy_SurfTemp(i,j) ) / ( 2.0_DP * DelTime )
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. * DelTime )
end if
end do
end do
end do
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagSeaIceMelt(i,j) ) then
xy_LatHeatFluxBySeaIceMelt(i,j) = xya_DelTempSoilTempLUVec(i,j,0)
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 PhyImplSDHSeaIceCorrection_OLD_Difficult
| Subroutine : | |||
| 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_SurfSnowB(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfMajCompIceB(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)
| ||
| 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_DSurfMajCompIceDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xy_DSoilMoistDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xy_DSurfSnowDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax) : | real(DP), intent(out )
| ||
| xy_LatHeatFluxBySnowMelt(0:imax-1, 1:jmax) : | real(DP), intent(out )
|
融雪による時間変化率の修正を行います.
Correction of tendencies due to melt of snow.
subroutine PhyImplSDHSnowMeltCorrection( xy_Ps, xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfSnowB, xy_SurfMajCompIceB, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xy_IndexCalcMethod, xyra_DelRadLFlux, xyz_Exner, xyr_Exner, xy_SurfTempTransCoef, xyza_ArgTempMtx, xyz_ArgTempVec, xyaa_ArgSurfMtx, xy_ArgSurfRH, xyaa_ArgSoilTempMtx, xya_ArgSoilTempVec, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DSurfMajCompIceDt, xy_DSoilMoistDt, xy_DSurfSnowDt, xy_LatHeatFluxByMajCompIceSubl, 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
! バケツモデル
! bucket model
!
use Bucket_Model, only: FlagBucketModel, FlagBucketModelSnow
! CO2 相変化
! Phase change of CO2
!
use co2_phase_change, only : LatentHeatCO2Subl, CO2Saturate2D
! 宣言文 ; Declaration statements
!
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_SurfSnowB (0:imax-1, 1:jmax)
! 積雪量.
! Surface snow amount.
real(DP), intent(in):: xy_SurfMajCompIceB(0:imax-1, 1:jmax)
!
! Surface major component ice 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(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_DSurfMajCompIceDt(0:imax-1, 1:jmax)
!
! Major component ice amount tendency (kg m-2 s-1)
real(DP), intent(inout):: xy_DSoilMoistDt (0:imax-1, 1:jmax)
! 土壌温度時間変化率 (kg m-2 s-1)
! Soil temperature tendency (kg m-2 s-1)
real(DP), intent(inout):: xy_DSurfSnowDt (0:imax-1, 1:jmax)
! 積雪率時間変化率 (kg m-2 s-1)
! Surface snow amount tendency (kg m-2 s-1)
real(DP), intent(out ):: xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax)
!
! Latent heat flux by major component ice sublimation
! (variable only for debug)
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_DSoilMoistDtSave(0:imax-1, 1:jmax)
real(DP):: xy_DSurfSnowDtSave (0:imax-1, 1:jmax)
real(DP):: xy_TempCond(0:imax-1, 1:jmax)
logical :: xy_FlagCalc(0:imax-1, 1:jmax)
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):: SurfSnowATentative
real(DP):: xy_SurfSnowATentativeSave(0:imax-1, 1:jmax)
real(DP) :: xy_TempMajCompCond(0:imax-1, 1:jmax)
real(DP) :: SurfMajCompIceATentative
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_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. FlagBucketModelSnow ) ) return
if ( kslmax == 0 ) then
call MessageNotify( 'E', module_name, 'kslmax <= 0 in PhyImplSDHSnowMeltCorrection.' )
else
xy_DSoilMoistDtSave = xy_DSoilMoistDt
xy_DSurfSnowDtSave = xy_DSurfSnowDt
xy_LatHeatFluxByMajCompIceSubl = 0.0_DP
!!$ if ( FlagMajCompPhaseChange ) then
!!$ call CO2Saturate2D( &
!!$ & xy_Ps, & ! (in)
!!$ & xy_TempMajCompCond & ! (out)
!!$ & )
!!$ end if
xy_TempCond = TempCondWater
LatentHeatLocal = LatentHeatFusion
xy_SurfSnowATentativeSave = xy_SurfSnowB + xy_DSurfSnowDt * ( 2.0_DP * DelTime )
!----------
! A case that a part of snow melts
!----------
do j = 1, jmax
do i = 0, imax-1
if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
SurfTempATentative = xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * 2.0d0 * DelTime
!!$ SurfSnowATentative = xy_SurfSnowB(i,j) &
!!$ & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
SurfSnowATentative = xy_SurfSnowATentativeSave(i,j)
if ( ( SurfSnowATentative > 0.0_DP ) .and. ( SurfTempATentative > xy_TempCond(i,j) ) ) then
xy_FlagCalc(i,j) = .true.
else
xy_FlagCalc(i,j) = .false.
end if
else
xy_FlagCalc(i,j) = .false.
end if
end do
end do
xyza_TempMtx = xyza_ArgTempMtx
xyz_TempVec = xyz_ArgTempVec
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
k = 1
xyza_TempMtx(i,j,k,-1) = 0.0_DP
xyz_TempVec(i,j,k) = - ( xyr_HeatFlux(i,j,k) - xyr_HeatFlux(i,j,k-1) )
xyz_TempVec(i,j,k) = xyz_TempVec(i,j,k) + CpDry * xy_SurfTempTransCoef(i,j) * ( xy_TempCond(i,j) - xy_SurfTemp(i,j) )
end if
end do
end do
!
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, 0) = 1.0_DP
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_SurfHeatCapacity(i,j) / ( 2. * DelTime ) + CpDry * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,0) - xyr_SoilTempTransCoef(i,j,0) ) * ( xy_TempCond(i,j) - xy_SurfTemp(i,j) )
end if
end do
end do
!
xyaa_SoilTempMtx = xyaa_ArgSoilTempMtx
xya_SoilTempVec = xya_ArgSoilTempVec
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
k = 1
xyaa_SoilTempMtx(i,j,k,-1) = 0.0_DP
xya_SoilTempVec(i,j,k) = - ( xyr_SoilHeatFlux(i,j,k) - xyr_SoilHeatFlux(i,j,k-1) )
xya_SoilTempVec(i,j,k) = xya_SoilTempVec(i,j,k) + xyr_SoilTempTransCoef(i,j,k-1) * ( xy_TempCond(i,j) - xy_SurfTemp(i,j) )
end if
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
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. * 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
xy_DSurfTempDt(i,j) = ( xy_TempCond(i,j) - xy_SurfTemp(i,j) ) / ( 2.0_DP * DelTime )
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. * DelTime )
end if
end do
end do
end do
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
xy_LatHeatFluxBySnowMelt(i,j) = xya_DelTempSoilTempLUVec(i,j,0)
!!$ xy_DSurfSnowDt(i,j) = - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
!!$ xy_DSoilMoistDt(i,j) = - xy_DSurfSnowDt(i,j)
xy_DSurfSnowDt(i,j) = xy_DSurfSnowDtSave(i,j) - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
xy_DSoilMoistDt(i,j) = xy_DSoilMoistDtSave(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 melts
!----------
do j = 1, jmax
do i = 0, imax-1
if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
if ( xy_FlagCalc(i,j) ) then
SurfSnowATentative = xy_SurfSnowB(i,j) + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
!!$ SurfSnowATentative = xy_SurfSnowATentativeSave(i,j)
!!$ if ( SurfSnowATentative < 0.0_DP ) then
!!$ SurfTempATentative = &
!!$ & xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * ( 2.0_DP * DelTime )
!!$ if ( SurfTempATentative < xy_TempCond(i,j) ) then
if ( SurfSnowATentative < 0.0_DP ) then
xy_FlagCalc(i,j) = .true.
else
xy_FlagCalc(i,j) = .false.
end if
else
xy_FlagCalc(i,j) = .false.
end if
else
xy_FlagCalc(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_FlagCalc(i,j) ) then
!!$ SurfSnowATentative = xy_SurfSnowB(i,j) &
!!$ & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
SurfSnowATentative = xy_SurfSnowATentativeSave(i,j)
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) - LatentHeatLocal * SurfSnowATentative / ( 2.0d0 * 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. * 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. * DelTime )
case ( IndexSeaIce )
! sea ice
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexSlabOcean )
! slab ocean
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexOceanPresSST )
! open ocean
xy_DSurfTempDt(i,j) = 0.
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. * 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
SurfSnowATentative = xy_SurfSnowATentativeSave(i,j)
xy_LatHeatFluxBySnowMelt(i,j) = LatentHeatLocal * SurfSnowATentative / ( 2.0d0 * DelTime )
xy_DSurfSnowDt(i,j) = xy_DSurfSnowDtSave(i,j) - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
xy_DSoilMoistDt(i,j) = xy_DSoilMoistDtSave(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
!----------
end if
!!$ ! 計算時間計測一時停止
!!$ ! Pause measurement of computation time
!!$ !
!!$ call TimesetClockStop( module_name )
end subroutine PhyImplSDHSnowMeltCorrection
| Subroutine : | |||
| 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_SurfSnowB(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfMajCompIceB(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)
| ||
| xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in )
| ||
| 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_DSurfMajCompIceDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xy_DSoilMoistDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xy_DSurfSnowDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax) : | real(DP), intent(out )
| ||
| xy_LatHeatFluxBySnowMelt(0:imax-1, 1:jmax) : | real(DP), intent(out )
|
融雪による時間変化率の修正を行います.
Correction of tendencies due to melt of snow.
subroutine PhyImplSDHSnowMeltCorrection_OLD( xy_Ps, xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfSnowB, xy_SurfMajCompIceB, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xy_IndexCalcMethod, xyra_DelRadLFlux, xyz_Exner, xyr_Exner, xy_SurfTempTransCoef, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DSurfMajCompIceDt, xy_DSoilMoistDt, xy_DSurfSnowDt, xy_LatHeatFluxByMajCompIceSubl, 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
! バケツモデル
! bucket model
!
use Bucket_Model, only: FlagBucketModel, FlagBucketModelSnow
! CO2 相変化
! Phase change of CO2
!
use co2_phase_change, only : LatentHeatCO2Subl, CO2Saturate2D
! 宣言文 ; Declaration statements
!
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_SurfSnowB (0:imax-1, 1:jmax)
! 積雪量.
! Surface snow amount.
real(DP), intent(in):: xy_SurfMajCompIceB(0:imax-1, 1:jmax)
!
! Surface major component ice 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 ):: 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_DSurfMajCompIceDt(0:imax-1, 1:jmax)
!
! Major component ice amount tendency (kg m-2 s-1)
real(DP), intent(inout):: xy_DSoilMoistDt (0:imax-1, 1:jmax)
! 土壌温度時間変化率 (kg m-2 s-1)
! Soil temperature tendency (kg m-2 s-1)
real(DP), intent(inout):: xy_DSurfSnowDt (0:imax-1, 1:jmax)
! 積雪率時間変化率 (kg m-2 s-1)
! Surface snow amount tendency (kg m-2 s-1)
real(DP), intent(out ):: xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax)
!
! Latent heat flux by major component ice sublimation
! (variable only for debug)
real(DP), intent(out ):: xy_LatHeatFluxBySnowMelt(0:imax-1, 1:jmax)
!
! Latent heat flux by melt (variable only for debug)
! 作業変数
! Work variables
!
real(DP):: TempCond
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):: SurfSnowATentative
real(DP) :: xy_TempMajCompCond(0:imax-1, 1:jmax)
real(DP) :: SurfMajCompIceATentative
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_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. FlagBucketModelSnow ) ) return
if ( kslmax == 0 ) then
call MessageNotify( 'E', module_name, 'kslmax <= 0 in PhyImplSDHSnowMeltCorrection.' )
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$
!!$ TempCond = TempCondWater
!!$
!!$ SurfSnowATentative = xy_SurfSnowB(i,j) &
!!$ & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
!!$ SurfTempATentative = xy_SurfTemp(i,j) &
!!$ & + xy_DSurfTempDt(i,j) * 2.0d0 * DelTime
!!$
!!$ if ( &
!!$ & ( xy_IndexCalcMethod(i,j) == IndexLand ) .and. &
!!$ & ( SurfSnowATentative > 0.0_DP ) .and. &
!!$ & ( SurfTempATentative > TempCond ) &
!!$ & ) then
!!$
!!$ ! if all snow is melting,
!!$
!!$ LatentHeatFluxByMelt = &
!!$ & SurfSnowATentative * LatentHeatFusion / ( 2.0d0 * DelTime )
!!$
!!$ SenHeatFluxA = &
!!$ & xyr_HeatFlux(i,j,0) &
!!$ & - CpDry * xyr_Exner(i,j,0) * xy_SurfTempTransCoef(i,j) &
!!$ & * ( xyz_DTempDt(i,j,1) / xyz_Exner(i,j,1) &
!!$ & - xy_DSurfTempDt(i,j) / xyr_Exner(i,j,0) ) * ( 2.0d0 * DelTime )
!!$ ! NOTICE: LatentHeatFlux^{n+1} = LatentHeatFlux^{n-1} due to mass
!!$ ! conservation
!!$ LatHeatFluxA = xy_SurfLatentHeatFlux(i,j)
!!$
!!$ CondHeatFluxA = xy_DeepSubSurfHeatFlux(i,j)
!!$
!!$ ValueAlpha = &
!!$ & xyr_RadSFlux(i,j,0) &
!!$ & + xyr_RadLFlux(i,j,0) &
!!$ & - xyra_DelRadLFlux(i,j,0,0) * xy_SurfTemp(i,j) &
!!$ & + xyra_DelRadLFlux(i,j,0,1) * xyz_DTempDt(i,j,1) &
!!$ & * ( 2.0d0 * DelTime ) &
!!$ & + SenHeatFluxA &
!!$ & + LatHeatFluxA &
!!$ & + LatentHeatFluxByMelt
!!$
!!$
!!$ SurfTempATentative = &
!!$ & xy_SurfHeatCapacity(i,j) / ( 2.0d0 * DelTime ) &
!!$ & * xy_SurfTemp(i,j) &
!!$ & - ValueAlpha + CondHeatFluxA
!!$ SurfTempATentative = &
!!$ & SurfTempATentative &
!!$ & / ( &
!!$ & xy_SurfHeatCapacity(i,j) / ( 2.0d0 * DelTime ) &
!!$ & + xyra_DelRadLFlux(i,j,0,0) &
!!$ & )
!!$
!!$ if ( SurfTempATentative >= TempCond ) then
!!$
!!$ xy_DSurfTempDt(i,j) = &
!!$ & ( SurfTempATentative - xy_SurfTemp(i,j) ) / ( 2.0d0 * DelTime )
!!$ xy_DSoilMoistDt(i,j) = &
!!$ & xy_DSoilMoistDt(i,j) + LatentHeatFluxByMelt / LatentHeatFusion
!!$ xy_DSurfSnowDt(i,j) = &
!!$ & xy_DSurfSnowDt (i,j) - LatentHeatFluxByMelt / LatentHeatFusion
!!$
!!$ else
!!$
!!$ ! if part of snow is melting,
!!$
!!$ SurfTempATentative = TempCond
!!$
!!$ ValueAlpha = &
!!$ & xyr_RadSFlux(i,j,0) &
!!$ & + xyr_RadLFlux(i,j,0) &
!!$ & + xyra_DelRadLFlux(i,j,0,0) * SurfTempATentative &
!!$ & - xyra_DelRadLFlux(i,j,0,0) * xy_SurfTemp(i,j) &
!!$ & + xyra_DelRadLFlux(i,j,0,1) * xyz_DTempDt(i,j,1) &
!!$ & * ( 2.0d0 * DelTime ) &
!!$ & + SenHeatFluxA &
!!$ & + LatHeatFluxA &
!!$ & - CondHeatFluxA
!!$
!!$ LatentHeatFluxByMelt = &
!!$ & xy_SurfHeatCapacity(i,j) &
!!$ & * ( SurfTempATentative - xy_SurfTemp(i,j) ) &
!!$ & / ( 2.0d0 * DelTime ) &
!!$ & - ValueAlpha
!!$
!!$ xy_DSurfTempDt(i,j) = &
!!$ & ( SurfTempATentative - xy_SurfTemp(i,j) ) / ( 2.0d0 * DelTime )
!!$ xy_DSoilMoistDt(i,j) = &
!!$ & xy_DSoilMoistDt(i,j) + LatentHeatFluxByMelt / LatentHeatFusion
!!$ xy_DSurfSnowDt(i,j) = &
!!$ & xy_DSurfSnowDt (i,j) - LatentHeatFluxByMelt / LatentHeatFusion
!!$
!!$ end if
!!$
!!$ SurfSnowATentative = xy_SurfSnowB(i,j) &
!!$ & + xy_DSurfSnowDt(i,j) * 2.0_DP * DelTime
!!$ if ( SurfSnowATentative < 0.0_DP ) then
!!$ xy_DSurfSnowDt(i,j) = - xy_SurfSnowB(i,j) / ( 2.0d0 * DelTime )
!!$ end if
!!$
!!$
!!$ ! Save latent heat flux for debug
!!$ xy_LatHeatFluxBySnowMelt(i,j) = LatentHeatFluxByMelt
!!$
!!$ else
!!$
!!$ ! Save latent heat flux for debug
!!$ xy_LatHeatFluxBySnowMelt(i,j) = 0.0_DP
!!$
!!$ end if
!!$
!!$ end do
!!$ end do
else
if ( FlagMajCompPhaseChange ) then
call CO2Saturate2D( xy_Ps, xy_TempMajCompCond )
end if
do j = 1, jmax
do i = 0, imax-1
if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
SoilTempATentative = xyz_SoilTemp(i,j,1) + xyz_DSoilTempDt(i,j,1) * 2.0d0 * DelTime
SurfMajCompIceATentative = xy_SurfMajCompIceB(i,j) + xy_DSurfMajCompIceDt(i,j) * 2.0d0 * DelTime
SenHeatFluxA = xyr_HeatFlux(i,j,0) - CpDry * xyr_Exner(i,j,0) * xy_SurfTempTransCoef(i,j) * ( xyz_DTempDt(i,j,1) / xyz_Exner(i,j,1) - xy_DSurfTempDt(i,j) / xyr_Exner(i,j,0) ) * ( 2.0d0 * DelTime )
!
! NOTICE: LatentHeatFlux^{n+1} = LatentHeatFlux^{n-1} due to mass
! conservation
!
LatHeatFluxA = xy_SurfLatentHeatFlux(i,j)
!
CondHeatFluxA = xyr_SoilHeatFlux(i,j,1) - xyr_SoilTempTransCoef(i,j,1) * ( xyz_DSoilTempDt(i,j,2) - xyz_DSoilTempDt(i,j,1) ) * ( 2.0d0 * DelTime )
! Phase change of major component
!
if ( FlagMajCompPhaseChange ) then
TempCond = xy_TempMajCompCond(i,j)
LatentHeatLocal = LatentHeatCO2Subl
if ( SoilTempATentative < TempCond ) then
! Ice accumulates
SoilTempATentative = TempCond
SurfTempATentative = ( r_SSDepth(0) - r_SSDepth(1) ) * ( r_SSDepth(0) - z_SSDepth(1) ) / xy_SoilHeatDiffCoef(i,j) * ( ( xy_SoilHeatCap(i,j) / ( 2.0d0 * DelTime ) + xy_SoilHeatDiffCoef(i,j) / ( ( r_SSDepth(0) - r_SSDepth(1) ) * ( r_SSDepth(0) - z_SSDepth(1) ) ) ) * SoilTempATentative - xy_SoilHeatCap(i,j) / ( 2.0d0 * DelTime ) * xyz_SoilTemp(i,j,1) - CondHeatFluxA / ( r_SSDepth(0) - r_SSDepth(1) ) )
LatentHeatFluxByMelt = - xy_SoilHeatDiffCoef(i,j) * ( SurfTempATentative - SoilTempATentative ) / ( r_SSDepth(0) - z_SSDepth(1) ) - xyr_RadSFlux(i,j,0) - ( xyr_RadLFlux(i,j,0) + xyra_DelRadLFlux(i,j,0,0) * ( SurfTempATentative - xy_SurfTemp(i,j) ) + xyra_DelRadLFlux(i,j,0,1) * xyz_DTempDt(i,j,1) * ( 2.0d0 * DelTime ) ) - SenHeatFluxA - LatHeatFluxA
xyz_DSoilTempDt(i,j,1) = ( SoilTempATentative - xyz_SoilTemp(i,j,1) ) / ( 2.0d0 * DelTime )
xy_DSurfTempDt(i,j) = ( SurfTempATentative - xy_SurfTemp(i,j) ) / ( 2.0d0 * DelTime )
xy_DSurfMajCompIceDt(i,j) = xy_DSurfMajCompIceDt(i,j) - LatentHeatFluxByMelt / LatentHeatLocal
! ! Save latent heat flux for debug
xy_LatHeatFluxByMajCompIceSubl(i,j) = LatentHeatFluxByMelt
else if ( SurfMajCompIceATentative > 0.0_DP ) then
! Ice sublimes
! if all ice is melting,
LatentHeatFluxByMelt = SurfMajCompIceATentative * LatentHeatLocal / ( 2.0d0 * DelTime )
ValueAlpha = xyr_RadSFlux(i,j,0) + xyr_RadLFlux(i,j,0) - xyra_DelRadLFlux(i,j,0,0) * xy_SurfTemp(i,j) + xyra_DelRadLFlux(i,j,0,1) * xyz_DTempDt(i,j,1) * ( 2.0d0 * DelTime ) + SenHeatFluxA + LatHeatFluxA + 0.0_DP + LatentHeatFluxByMelt
ValueAlpha = ValueAlpha * ( r_SSDepth(0) - z_SSDepth(1) ) / xy_SoilHeatDiffCoef(i,j)
SurfTempATentative = - ( xy_SoilHeatCap(i,j) / ( 2.0d0 * DelTime ) + xy_SoilHeatDiffCoef(i,j) / ( ( r_SSDepth(0) - r_SSDepth(1) ) * ( r_SSDepth(0) - z_SSDepth(1) ) ) ) * ValueAlpha + xy_SoilHeatCap(i,j) / ( 2.0d0 * DelTime ) * xyz_SoilTemp(i,j,1) + CondHeatFluxA / ( r_SSDepth(0) - r_SSDepth(1) )
SurfTempATentative = SurfTempATentative / ( xy_SoilHeatCap(i,j) / ( 2.0d0 * DelTime ) * ( 1.0d0 + ( r_SSDepth(0) - z_SSDepth(1) ) / xy_SoilHeatDiffCoef(i,j) * xyra_DelRadLFlux(i,j,0,0) ) + xyra_DelRadLFlux(i,j,0,0) / ( r_SSDepth(0) - r_SSDepth(1) ) )
SoilTempATentative = ( 1.0d0 + ( r_SSDepth(0) - z_SSDepth(1) ) / xy_SoilHeatDiffCoef(i,j) * xyra_DelRadLFlux(i,j,0,0) ) * SurfTempATentative + ValueAlpha
if ( SoilTempATentative >= TempCond ) then
xyz_DSoilTempDt(i,j,1) = ( SoilTempATentative - xyz_SoilTemp(i,j,1) ) / ( 2.0d0 * DelTime )
xy_DSurfTempDt(i,j) = ( SurfTempATentative - xy_SurfTemp(i,j) ) / ( 2.0d0 * DelTime )
xy_DSurfMajCompIceDt(i,j) = xy_DSurfMajCompIceDt (i,j) - LatentHeatFluxByMelt / LatentHeatLocal
else
! if part of ice sublimes,
SoilTempATentative = TempCond
SurfTempATentative = ( r_SSDepth(0) - r_SSDepth(1) ) * ( r_SSDepth(0) - z_SSDepth(1) ) / xy_SoilHeatDiffCoef(i,j) * ( ( xy_SoilHeatCap(i,j) / ( 2.0d0 * DelTime ) + xy_SoilHeatDiffCoef(i,j) / ( ( r_SSDepth(0) - r_SSDepth(1) ) * ( r_SSDepth(0) - z_SSDepth(1) ) ) ) * SoilTempATentative - xy_SoilHeatCap(i,j) / ( 2.0d0 * DelTime ) * xyz_SoilTemp(i,j,1) - CondHeatFluxA / ( r_SSDepth(0) - r_SSDepth(1) ) )
LatentHeatFluxByMelt = - xy_SoilHeatDiffCoef(i,j) * ( SurfTempATentative - SoilTempATentative ) / ( r_SSDepth(0) - z_SSDepth(1) ) - xyr_RadSFlux(i,j,0) - ( xyr_RadLFlux(i,j,0) + xyra_DelRadLFlux(i,j,0,0) * ( SurfTempATentative - xy_SurfTemp(i,j) ) + xyra_DelRadLFlux(i,j,0,1) * xyz_DTempDt(i,j,1) * ( 2.0d0 * DelTime ) ) - SenHeatFluxA - LatHeatFluxA
xyz_DSoilTempDt(i,j,1) = ( SoilTempATentative - xyz_SoilTemp(i,j,1) ) / ( 2.0d0 * DelTime )
xy_DSurfTempDt(i,j) = ( SurfTempATentative - xy_SurfTemp(i,j) ) / ( 2.0d0 * DelTime )
xy_DSurfMajCompIceDt(i,j) = xy_DSurfMajCompIceDt (i,j) - LatentHeatFluxByMelt / LatentHeatLocal
end if
SurfMajCompIceATentative = xy_SurfMajCompIceB(i,j) + xy_DSurfMajCompIceDt(i,j) * 2.0d0 * DelTime
if ( SurfMajCompIceATentative < 0.0_DP ) then
! This should be unnecessary.
xy_DSurfMajCompIceDt(i,j) = - xy_SurfMajCompIceB(i,j) / ( 2.0d0 * DelTime )
end if
! ! Save latent heat flux for debug
xy_LatHeatFluxByMajCompIceSubl(i,j) = LatentHeatFluxByMelt
else
! ! Save latent heat flux for debug
xy_LatHeatFluxByMajCompIceSubl(i,j) = 0.0_DP
end if
else
xy_LatHeatFluxByMajCompIceSubl(i,j) = 0.0_DP
end if ! if ( FlagMajCompPhaseChange ) then
! Phase change of H2O
!
TempCond = TempCondWater
LatentHeatLocal = LatentHeatFusion
SurfSnowATentative = xy_SurfSnowB(i,j) + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
if ( ( SurfSnowATentative > 0.0_DP ) .and. ( SoilTempATentative > TempCond ) ) then
! if all snow is melting,
LatentHeatFluxByMelt = SurfSnowATentative * LatentHeatLocal / ( 2.0d0 * DelTime )
ValueAlpha = xyr_RadSFlux(i,j,0) + xyr_RadLFlux(i,j,0) - xyra_DelRadLFlux(i,j,0,0) * xy_SurfTemp(i,j) + xyra_DelRadLFlux(i,j,0,1) * xyz_DTempDt(i,j,1) * ( 2.0d0 * DelTime ) + SenHeatFluxA + LatHeatFluxA + xy_LatHeatFluxByMajCompIceSubl(i,j) + LatentHeatFluxByMelt
ValueAlpha = ValueAlpha * ( r_SSDepth(0) - z_SSDepth(1) ) / xy_SoilHeatDiffCoef(i,j)
SurfTempATentative = - ( xy_SoilHeatCap(i,j) / ( 2.0d0 * DelTime ) + xy_SoilHeatDiffCoef(i,j) / ( ( r_SSDepth(0) - r_SSDepth(1) ) * ( r_SSDepth(0) - z_SSDepth(1) ) ) ) * ValueAlpha + xy_SoilHeatCap(i,j) / ( 2.0d0 * DelTime ) * xyz_SoilTemp(i,j,1) + CondHeatFluxA / ( r_SSDepth(0) - r_SSDepth(1) )
SurfTempATentative = SurfTempATentative / ( xy_SoilHeatCap(i,j) / ( 2.0d0 * DelTime ) * ( 1.0d0 + ( r_SSDepth(0) - z_SSDepth(1) ) / xy_SoilHeatDiffCoef(i,j) * xyra_DelRadLFlux(i,j,0,0) ) + xyra_DelRadLFlux(i,j,0,0) / ( r_SSDepth(0) - r_SSDepth(1) ) )
SoilTempATentative = ( 1.0d0 + ( r_SSDepth(0) - z_SSDepth(1) ) / xy_SoilHeatDiffCoef(i,j) * xyra_DelRadLFlux(i,j,0,0) ) * SurfTempATentative + ValueAlpha
if ( SoilTempATentative >= TempCond ) then
xyz_DSoilTempDt(i,j,1) = ( SoilTempATentative - xyz_SoilTemp(i,j,1) ) / ( 2.0d0 * DelTime )
xy_DSurfTempDt(i,j) = ( SurfTempATentative - xy_SurfTemp(i,j) ) / ( 2.0d0 * DelTime )
xy_DSoilMoistDt(i,j) = xy_DSoilMoistDt(i,j) + LatentHeatFluxByMelt / LatentHeatLocal
xy_DSurfSnowDt(i,j) = xy_DSurfSnowDt (i,j) - LatentHeatFluxByMelt / LatentHeatLocal
else
! if part of snow is melting,
SoilTempATentative = TempCond
SurfTempATentative = ( r_SSDepth(0) - r_SSDepth(1) ) * ( r_SSDepth(0) - z_SSDepth(1) ) / xy_SoilHeatDiffCoef(i,j) * ( ( xy_SoilHeatCap(i,j) / ( 2.0d0 * DelTime ) + xy_SoilHeatDiffCoef(i,j) / ( ( r_SSDepth(0) - r_SSDepth(1) ) * ( r_SSDepth(0) - z_SSDepth(1) ) ) ) * SoilTempATentative - xy_SoilHeatCap(i,j) / ( 2.0d0 * DelTime ) * xyz_SoilTemp(i,j,1) - CondHeatFluxA / ( r_SSDepth(0) - r_SSDepth(1) ) )
LatentHeatFluxByMelt = - xy_SoilHeatDiffCoef(i,j) * ( SurfTempATentative - SoilTempATentative ) / ( r_SSDepth(0) - z_SSDepth(1) ) - xyr_RadSFlux(i,j,0) - ( xyr_RadLFlux(i,j,0) + xyra_DelRadLFlux(i,j,0,0) * ( SurfTempATentative - xy_SurfTemp(i,j) ) + xyra_DelRadLFlux(i,j,0,1) * xyz_DTempDt(i,j,1) * ( 2.0d0 * DelTime ) ) - SenHeatFluxA - LatHeatFluxA
xyz_DSoilTempDt(i,j,1) = ( SoilTempATentative - xyz_SoilTemp(i,j,1) ) / ( 2.0d0 * DelTime )
xy_DSurfTempDt(i,j) = ( SurfTempATentative - xy_SurfTemp(i,j) ) / ( 2.0d0 * DelTime )
xy_DSoilMoistDt(i,j) = xy_DSoilMoistDt(i,j) + LatentHeatFluxByMelt / LatentHeatLocal
xy_DSurfSnowDt(i,j) = xy_DSurfSnowDt (i,j) - LatentHeatFluxByMelt / LatentHeatLocal
end if
SurfSnowATentative = xy_SurfSnowB(i,j) + xy_DSurfSnowDt(i,j) * 2.0_DP * DelTime
if ( SurfSnowATentative < 0.0_DP ) then
! This should be unnecessary.
xy_DSurfSnowDt(i,j) = - xy_SurfSnowB(i,j) / ( 2.0d0 * DelTime )
end if
! Save latent heat flux for debug
xy_LatHeatFluxBySnowMelt(i,j) = LatentHeatFluxByMelt
else
! Save latent heat flux for debug
xy_LatHeatFluxBySnowMelt(i,j) = 0.0_DP
end if
else ! if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
xy_LatHeatFluxByMajCompIceSubl(i,j) = 0.0_DP
xy_LatHeatFluxBySnowMelt (i,j) = 0.0_DP
end if ! if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
end do
end do
end if
!!$ ! 計算時間計測一時停止
!!$ ! Pause measurement of computation time
!!$ !
!!$ call TimesetClockStop( module_name )
end subroutine PhyImplSDHSnowMeltCorrection_OLD
| Subroutine : | |||
| 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_SurfSnowB(0:imax-1, 1:jmax) : | real(DP), intent(in)
| ||
| xy_SurfMajCompIceB(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)
| ||
| xyza_ArgTempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1) : | 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 )
| ||
| 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_DSurfMajCompIceDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xy_DSoilMoistDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xy_DSurfSnowDt(0:imax-1, 1:jmax) : | real(DP), intent(inout)
| ||
| xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax) : | real(DP), intent(out )
| ||
| xy_LatHeatFluxBySnowMelt(0:imax-1, 1:jmax) : | real(DP), intent(out )
|
融雪による時間変化率の修正を行います.
Correction of tendencies due to melt of snow.
subroutine PhyImplSDHSnowMeltCorrection_OLD2( xy_Ps, xyr_HeatFlux, xy_SurfLatentHeatFlux, xyr_SoilHeatFlux, xyr_SoilTempTransCoef, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfSnowB, xy_SurfMajCompIceB, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xy_IndexCalcMethod, xyra_DelRadLFlux, xyz_Exner, xyr_Exner, xy_SurfTempTransCoef, xyza_ArgTempMtx, xyaa_ArgSurfMtx, xy_ArgSurfRH, xyaa_ArgSoilTempMtx, xyz_DTempDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DSurfMajCompIceDt, xy_DSoilMoistDt, xy_DSurfSnowDt, xy_LatHeatFluxByMajCompIceSubl, 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
! バケツモデル
! bucket model
!
use Bucket_Model, only: FlagBucketModel, FlagBucketModelSnow
! CO2 相変化
! Phase change of CO2
!
use co2_phase_change, only : LatentHeatCO2Subl, CO2Saturate2D
! 宣言文 ; Declaration statements
!
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_SurfSnowB (0:imax-1, 1:jmax)
! 積雪量.
! Surface snow amount.
real(DP), intent(in):: xy_SurfMajCompIceB(0:imax-1, 1:jmax)
!
! Surface major component ice 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 ):: 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(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_DSurfMajCompIceDt(0:imax-1, 1:jmax)
!
! Major component ice amount tendency (kg m-2 s-1)
real(DP), intent(inout):: xy_DSoilMoistDt (0:imax-1, 1:jmax)
! 土壌温度時間変化率 (kg m-2 s-1)
! Soil temperature tendency (kg m-2 s-1)
real(DP), intent(inout):: xy_DSurfSnowDt (0:imax-1, 1:jmax)
! 積雪率時間変化率 (kg m-2 s-1)
! Surface snow amount tendency (kg m-2 s-1)
real(DP), intent(out ):: xy_LatHeatFluxByMajCompIceSubl(0:imax-1, 1:jmax)
!
! Latent heat flux by major component ice sublimation
! (variable only for debug)
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_TempCond(0:imax-1, 1:jmax)
logical :: xy_FlagCalc(0:imax-1, 1:jmax)
real(DP):: xyza_TempMtx(0:imax-1, 1:jmax, 1:kmax, -1:1)
! 温度陰解行列.
! Implicit matrix 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):: 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):: SurfSnowATentative
!!$ real(DP):: xy_SurfSnowATentativeSave(0:imax-1, 1:jmax)
real(DP) :: xy_TempMajCompCond(0:imax-1, 1:jmax)
real(DP) :: SurfMajCompIceATentative
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_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. FlagBucketModelSnow ) ) return
if ( kslmax == 0 ) then
call MessageNotify( 'E', module_name, 'kslmax <= 0 in PhyImplSDHSnowMeltCorrection.' )
else
xy_LatHeatFluxByMajCompIceSubl = 0.0_DP
!!$ if ( FlagMajCompPhaseChange ) then
!!$ call CO2Saturate2D( &
!!$ & xy_Ps, & ! (in)
!!$ & xy_TempMajCompCond & ! (out)
!!$ & )
!!$ end if
xy_TempCond = TempCondWater
LatentHeatLocal = LatentHeatFusion
!!$ xy_SurfSnowATentativeSave = xy_SurfSnowB + xy_DSurfSnowDt * ( 2.0_DP * DelTime )
!----------
! A case that all snow melts
!----------
do j = 1, jmax
do i = 0, imax-1
if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
SurfTempATentative = xy_SurfTemp(i,j) + xy_DSurfTempDt(i,j) * 2.0d0 * DelTime
SurfSnowATentative = xy_SurfSnowB(i,j) + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
!!$ SurfSnowATentative = xy_SurfSnowATentativeSave(i,j)
if ( ( SurfSnowATentative > 0.0_DP ) .and. ( SurfTempATentative > xy_TempCond(i,j) ) ) then
xy_FlagCalc(i,j) = .true.
else
xy_FlagCalc(i,j) = .false.
end if
else
xy_FlagCalc(i,j) = .false.
end if
end do
end do
xyza_TempMtx = xyza_ArgTempMtx
xyaa_SoilTempMtx = xyaa_ArgSoilTempMtx
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
!!$ SurfSnowATentative = xy_SurfSnowATentativeSave(i,j)
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) - LatentHeatLocal * SurfSnowATentative / ( 2.0d0 * DelTime )
end if
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) = - ( xyr_SoilHeatFlux(:,:,k) - xyr_SoilHeatFlux(:,:,k-1) )
end do
k = 0
xya_DelTempSoilTempLUVec(:,:,k) = xy_SurfRH
do k = 1, kmax
xya_DelTempSoilTempLUVec(:,:,k) = - ( xyr_HeatFlux(:,:,k) - xyr_HeatFlux(:,:,k-1) )
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. * 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. * DelTime )
case ( IndexSeaIce )
! sea ice
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexSlabOcean )
! slab ocean
xy_DSurfTempDt(i,j) = xya_DelTempSoilTempLUVec(i,j,0) / ( 2. * DelTime )
case ( IndexOceanPresSST )
! open ocean
xy_DSurfTempDt(i,j) = 0.
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. * 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
!!$ SurfSnowATentative = xy_SurfSnowATentativeSave(i,j)
xy_LatHeatFluxBySnowMelt(i,j) = LatentHeatLocal * SurfSnowATentative / ( 2.0d0 * DelTime )
xy_DSurfSnowDt(i,j) = - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
xy_DSoilMoistDt(i,j) = - xy_DSurfSnowDt(i,j)
else
xy_LatHeatFluxBySnowMelt(i,j) = 0.0_DP
end if
end do
end do
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ 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_SurfTemp(i,j) /) )
!!$ end if
!!$ end do
!!$ end do
!----------
! A case that a part of snow melts
!----------
do j = 1, jmax
do i = 0, imax-1
if ( xy_IndexCalcMethod(i,j) == IndexLand ) then
if ( xy_FlagCalc(i,j) ) then
!!$ SurfSnowATentative = xy_SurfSnowB(i,j) &
!!$ & + xy_DSurfSnowDt(i,j) * 2.0d0 * DelTime
!!$ SurfSnowATentative = xy_SurfSnowATentativeSave(i,j)
!!$ if ( SurfSnowATentative < 0.0_DP ) then
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.
else
xy_FlagCalc(i,j) = .false.
end if
else
xy_FlagCalc(i,j) = .false.
end if
else
xy_FlagCalc(i,j) = .false.
end if
end do
end do
xyza_TempMtx = xyza_ArgTempMtx
xyaa_SoilTempMtx = xyaa_ArgSoilTempMtx
xyaa_SurfMtx = xyaa_ArgSurfMtx
xy_SurfRH = xy_ArgSurfRH
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
k = 1
xyza_TempMtx(i,j,k,-1) = 0.0_DP
xyaa_SurfMtx(i,j,0, 0) = 1.0_DP
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_SurfHeatCapacity(i,j) / ( 2. * DelTime ) + CpDry * xy_SurfTempTransCoef(i,j) + xyra_DelRadLFlux(i,j,0,0) - xyr_SoilTempTransCoef(i,j,0) ) * ( xy_TempCond(i,j) - xy_SurfTemp(i,j) )
k = 1
xyaa_SoilTempMtx(:,:,k,-1) = 0.0_DP
end if
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) = - ( xyr_SoilHeatFlux(:,:,k) - xyr_SoilHeatFlux(:,:,k-1) )
end do
k = 1
xya_DelTempSoilTempLUVec(:,:,-k) = xya_DelTempSoilTempLUVec(:,:,-k) + xyr_SoilTempTransCoef(:,:,k-1) * ( xy_TempCond - xy_SurfTemp )
k = 0
xya_DelTempSoilTempLUVec(:,:,k) = xy_SurfRH
do k = 1, kmax
xya_DelTempSoilTempLUVec(:,:,k) = - ( xyr_HeatFlux(:,:,k) - xyr_HeatFlux(:,:,k-1) )
end do
k = 1
xya_DelTempSoilTempLUVec(:,:,k) = xya_DelTempSoilTempLUVec(:,:,k) + CpDry * xy_SurfTempTransCoef * ( xy_TempCond - xy_SurfTemp )
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. * 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
xy_DSurfTempDt(i,j) = ( xy_TempCond(i,j) - xy_SurfTemp(i,j) ) / ( 2.0_DP * DelTime )
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. * DelTime )
end if
end do
end do
end do
do j = 1, jmax
do i = 0, imax-1
if ( xy_FlagCalc(i,j) ) then
xy_LatHeatFluxBySnowMelt(i,j) = xya_DelTempSoilTempLUVec(i,j,0)
xy_DSurfSnowDt(i,j) = - xy_LatHeatFluxBySnowMelt(i,j) / LatentHeatLocal
xy_DSoilMoistDt(i,j) = - xy_DSurfSnowDt(i,j)
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
end if
end do
end do
!----------
end if
!!$ ! 計算時間計測一時停止
!!$ ! Pause measurement of computation time
!!$ !
!!$ call TimesetClockStop( module_name )
end subroutine PhyImplSDHSnowMeltCorrection_OLD2
| Constant : | |||
| module_name = ‘phy_implicit_sdh‘ : | character(*), parameter
|
| Variable : | |||
| phy_implicit_sdh_inited = .false. : | logical, save
|
| Constant : | |||
| version = ’$Name: dcpam5-20120226 $’ // ’$Id: phy_implicit_sdh.f90,v 1.10 2012-01-20 00:25:40 yot Exp $’ : | character(*), parameter
|