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