!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
!
!= Module dcpam_ape_physics_mod
!
!   * Developers: Morikawa Yasuhiro, Yamada Yukiko
!   * Version: $Id: dcpam_ape_physics.f90,v 1.11 2007/08/27 14:45:43 momoko Exp $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview
!
!DCPAM-APE Physics Main module
!
!== Error Handling
!
!== Known Bugs
!
!== Note
! * The module for corrections of Surface pressure and specific humidity 
!   due to the change of atmospheric mass can be used if you nedd.
!   When you use the module, delete the comment at the part of
!   call physics_ps_correction( ....
!
!== Future Plans
! * xy_SurfTemp はメインで管理するようにしてメインに返すようにすること.
!
!== Histrory
!  2007-5-12 M.Ishiwatari 放射モジュールのiterface を変更. xyr_Temp 追加.
!  2007-5-17 M.Ishiwatari 変数 xyr_RadLFluxCorrected を新設
!

module dcpam_ape_physics_mod
  use type_mod,      only : STRING
  use grid_3d_mod,   only: im, jm, km
  implicit none
  private
  public :: dcpam_ape_physics       ! subroutines
  character(STRING),parameter:: version = &
       & '$Id: dcpam_ape_physics.f90,v 1.11 2007/08/27 14:45:43 momoko Exp $'
  character(STRING),parameter:: tagname = '$Name:  $'
  logical,save:: dcpam_ape_physics_initialized = .false.


contains

  subroutine dcpam_ape_physics( Dims, Vars_a )
    !
    ! Physics main module
    ! DCPAM-APE Physics Main module
    !
    use dycore_type_mod, only: DYCORE_VARS, DYCORE_DIMS, STRING, INTKIND, &
         &                     REKIND, DBKIND
    use grid_3d_mod,     only: im, jm, km
    use time_mod,    only: DelTime
    use spml_mod,    only: wa_Div_xya_xya, xya_wa, wa_xya, xy_w, w_xy, &
         & xya_GradLon_wa,xya_GradLat_wa, wa_LaplaInv_wa
    use constants_mod, only: R0
    use io_gt4_out_mod , only: io_gt4_out_Put
    use dc_trace,        only: BeginSub, EndSub, DbgMessage
    use physics_interpolate_mod, only: physics_interpolate_temp, &
         &                             physics_interpolate_geopot
    use physics_negq_mod, only: physics_negq
    use physics_lscond_mod, only: physics_lscond

    ! switch of cumulus parameterization
    use physics_cumulus_adjust_mod , only: physics_cumulus_adjust
    !use physics_cumulus_adjust_muchwater_mod , only: physics_cumulus_adjust
    !use physics_cumulus_none_mod , only: physics_cumulus_adjust

    use physics_dryadjust_mod , only: physics_dryadjust
    use physics_ground_mod, only: physics_ground
    use physics_radiation_main_mod, only: physics_radiation_main, &
         & physics_radiation_deltemp
    use physics_verdiff_main_mod, only: physics_verdiff_main
    use physics_surface_main_mod, only: physics_surface_main
    use physics_implicit_mod, only: physics_implicit_init, &
      & physics_implicit_integrate, &
      & physics_implicit_fluxcorrection

    ! switch of surface pressure correction 
    use physics_ps_correction_none_mod, only: physics_ps_correction
    !use physics_ps_correction_mod, only: physics_ps_correction

    ! switch of vertical filter
    use phys_vfilter_none_mod, only: phys_vfilter
    !use phys_vfilter_conserve_mod, only: phys_vfilter

    use dycore_grid_mod, only: nm

    implicit none
    type(DYCORE_DIMS), intent(in)   :: Dims   ! 次元データ全種
    type(DYCORE_VARS), intent(inout):: Vars_a ! 格子点データ全種(t+Δt)
    character(STRING),  parameter:: subname = "dcpam_ape_physics"
    real(DBKIND) :: xyr_Temp(im,jm,km+1)    ! 温度 (半整数)
    real(DBKIND) :: xyz_Press(im,jm,km)      ! 気圧 (整数)
    real(DBKIND) :: xyr_Press(im,jm,km+1)    ! 気圧 (半整数)
    real(DBKIND) :: xyz_GeoPot(im,jm,km)     ! ジオポテンシャル(整数)
    real(DBKIND) :: xyr_GeoPot(im,jm,km+1)   ! ジオポテンシャル(半整数)
    real(DBKIND) :: xyz_DNegQvap1Dt(im,jm,km) ! 負の比湿除去 (比湿変化率)
    real(DBKIND) :: xyz_DNegQvap2Dt(im,jm,km) ! 負の比湿除去 (比湿変化率)
    real(DBKIND) :: xyz_DLscTempDt(im,jm,km) ! 大規模凝結による温度変化率
    real(DBKIND) :: xyz_DLscQvapDt(im,jm,km) ! 大規模凝結による比湿変化率
    real(DBKIND) :: xy_LscRain(im,jm)        ! 大規模凝結による降水量
    real(DBKIND) :: xyz_DCumulusTempDt(im,jm,km) ! 積雲スキームによる温度変化率
    real(DBKIND) :: xyz_DCumulusQvapDt(im,jm,km) ! 積雲スキームによる比湿変化率
    real(DBKIND) :: xy_CumulusRain(im,jm)    ! 積雲スキームによる降水量
    real(DBKIND) :: xyz_DDryTempDt(im,jm,km) ! 乾燥対流調節による温度変化率
    real(DBKIND) :: xy_Rain(im,jm)           ! cumulus + lsc 降水量
    real(DBKIND) :: xy_Ps_b(im,jm)              ! 地表面気圧

    real(DBKIND) :: xyz_DRadLTempDt(im,jm,km)  ! 長波加熱率
    real(DBKIND) :: xyz_DRadSTempDt(im,jm,km)   ! 短波加熱率
! (2007-5-17 石渡) 緊急措置
!    real(DBKIND) :: xy_SurfTemp(im,jm)          ! 地表面温度
    real(DBKIND),allocatable,save :: xy_SurfTemp(:,:)          ! 地表面温度
    real(DBKIND) :: xy_SurfAlbedo(im,jm)        ! 地表アルベド
    real(DBKIND) :: xy_SurfHumidCoeff(im,jm)    ! 地表湿潤度
    real(DBKIND) :: xy_SurfRoughLength(im,jm)   ! 地表粗度長
    real(DBKIND) :: xy_SurfHeatCapacity(im,jm)  ! 地表熱容量
    real(DBKIND) :: xy_GroundTempFlux(im,jm)     ! 地中熱フラックス
    real(DBKIND) :: xyr_RadLFlux(im,jm,km+1)           ! 長波フラックス
    real(DBKIND) :: xyr_RadLFluxCorrected(im,jm,km+1) ! 補正した長波フラックス
    real(DBKIND) :: xyo_SurfRadLMatrix(im,jm,-1:1)     ! Ｔ陰解行列：地表
    real(DBKIND) :: xyro_DelRadLFlux(im,jm,km+1,0:1)   ! 長波地表温度変化
    real(DBKIND) :: xyr_RadSFlux(im,jm,km+1)           ! 日射フラックス
    real(DBKIND) :: xyr_VelLonFlux(im,jm,km+1)         ! 速度経度成分フラックス
    real(DBKIND) :: xyr_VelLatFlux(im,jm,km+1) ! 速度緯度成分フラックス
    real(DBKIND) :: xyr_TempFlux(im,jm,km+1)           ! 温度フラックス
    real(DBKIND) :: xyr_QvapFlux(im,jm,km+1)           ! 比湿フラックス
    real(DBKIND) :: xyzo_VelMatrix(im,jm,km,-1:1)      ! 速度陰解行列
    real(DBKIND) :: xyzo_TempMatrix(im,jm,0:km,-1:1)   ! 温度陰解行列
    real(DBKIND) :: xyzo_QvapMatrix(im,jm,km,-1:1)     ! 比湿陰解行列
    real(DBKIND) :: xy_SurfVelMatrix(im,jm)            !  速度陰解行列: 地表
    real(DBKIND) :: xyoo_SurfTempMatrix(im,jm,0:1,-1:1) ! 温度陰解行列: 地表
    real(DBKIND) :: xyoo_SurfQvapMatrix(im,jm,0:1,-1:1) ! 比湿陰解行列: 地表
    real(DBKIND) :: xyz_DVerdiffVelLonDt(im,jm,km) ! 経度成分 鉛直拡散加速度
    real(DBKIND) :: xyz_DVerdiffVelLatDt(im,jm,km) ! 緯度成分 鉛直拡散加速度
    real(DBKIND) :: xyz_DVerdiffTempDt(im,jm,km)       ! 鉛直拡散加熱率
    real(DBKIND) :: xyz_DVerdiffSurfTempDt(im,jm)      ! 地表面 鉛直拡散加熱率
    real(DBKIND) :: xyz_DVerdiffQvapDt(im,jm,km)       ! 鉛直拡散加湿率
    integer(INTKIND) :: xy_SurfCondition(im,jm)    ! 地表状態
    integer(INTKIND) :: k
    real(DBKIND) :: wz_Psi_a((nm+1)*(nm+1), km) 
    real(DBKIND) :: wz_Chi_a((nm+1)*(nm+1), km)
    real(DBKIND) :: xyz_DTempDtVertFiltCons(im,jm,km) ! 温度変化率(鉛直フィルター)
    real(DBKIND) :: xyz_DVelLonDtVertFiltCons(im,jm,km) ! Ｕ変化率(鉛直フィルター)
    real(DBKIND) :: xyz_DVelLatDtVertFilstCons(im,jm,km) ! Ｖ変化率(鉛直フィルター)

    continue


    !   開始処理
    call BeginSub(subname)
  
    ! (2007-5-17 石渡) 緊急措置で追加.
    if (dcpam_ape_physics_initialized) then
    else
      allocate( xy_SurfTemp(im,jm) )
      dcpam_ape_physics_initialized = .true.
    endif


    ! 1. 物理過程演算用の変数を初期化
    xyz_DNegQvap1Dt = 0.0d0
    xyz_DNegQvap2Dt = 0.0d0
    xyz_DLscTempDt = 0.0d0
    xyz_DLscQvapDt = 0.0d0
    xy_LscRain     = 0.0d0
    xyz_DCumulusTempDt = 0.0d0
    xyz_DCumulusQvapDt = 0.0d0
    xy_CumulusRain = 0.0d0
    xyz_DDryTempDt = 0.0d0
    xy_Rain        = 0.0d0

    xy_Ps_b = Vars_a%xy_Ps

    !  2. 地表条件設定
    call physics_ground( &
      & xy_SurfTemp        , &  !(inout) 地表温度
      & xy_SurfAlbedo, &        !(out) 地表アルベド
      & xy_SurfHumidCoeff, &    !(out) 地表湿潤度
      & xy_SurfRoughLength, &   !(out) 地表粗度長
      & xy_SurfCondition, &     !(out) 地表状態. 
      & xy_SurfHeatCapacity, &  !(out) 地表熱容量
      & xy_GroundTempFlux, &    !(out) 地中熱フラックス
      &  Dims%y_Lat%a_Dim     ) ! (in) 経度座標

    !----------------------------------------------------------------
    !  3. 温度半整数 sigma 補間, 気圧と高度の算出 (1)
    !----------------------------------------------------------------

    !----- 温度半整数 sigma 補間 -----
    call physics_interpolate_temp( &
      & xyr_Temp              , & ! (out) 温度 (半整数)
      & Vars_a%xyz_Temp       , & ! (in) 温度 (整数) (t+Δt)
      & Dims%z_Sigma%a_Dim    , & ! (in) σレベル(整数)座標
      & Dims%r_Sigma%a_Dim    )   ! (in) σレベル(半整数)座標

    !----- 気圧と高度の算出 -----
    call physics_interpolate_geopot( &
      & xyz_Press            , & ! (out) 圧力 (整数)
      & xyr_Press            , & ! (out) 圧力 (半整数)
      & xyz_GeoPot           , & ! (out) ジオポテンシャル(整数)
      & xyr_GeoPot           , & ! (out) ジオポテンシャル(半整数)
      & Vars_a%xy_Ps         , & ! (in) 地表面気圧 (t+Δt)
      & Vars_a%xyz_Temp      , & ! (in) 温度 (整数) (t+Δt)
      & xyr_Temp             , & ! (in) 温度 (半整数)
      & Dims%z_Sigma%a_Dim   , & ! (in) σレベル(整数)座標
      & Dims%r_Sigma%a_Dim   )   ! (in) σレベル(半整数)座標

    !----------------------------------------------------------------
    !  4. 負の水蒸気除去(1)
    !----------------------------------------------------------------
    call physics_negq( &
      & Vars_a%xyz_Qvap      , & ! (inout) 比湿 (整数) (t+Δt)
      & xyz_DNegQvap1Dt      , & ! (inout) 比湿変化 (整数)
      & xyr_Press            , & ! (in) 圧力 (半整数)
      & 2.0d0*DelTime  )         ! (in) 2Δt

    !----------------------------------------------------------------
    !  5. 湿潤過程 (積雲)
    !----------------------------------------------------------------
    call physics_cumulus_adjust( &
      & Vars_a%xyz_Temp     , & ! (inout) 温度
      & Vars_a%xyz_Qvap     , & ! (inout) 比湿
      & xy_CumulusRain      , & ! (out) 降水量
      & xyz_DCumulusTempDt  , & ! (out) 温度変化率
      & xyz_DCumulusQvapDt  , & ! (out) 比湿変化率
      & xyz_Press           , & ! (in) 気圧 (整数)
      & xyr_Press           , & ! (in) 気圧 (半整数)
      & 2.0d0*DelTime )         ! (in) 2Δt

    !----------------------------------------------------------------
    !  6. 湿潤過程 (大規模凝結)
    !----------------------------------------------------------------
    call physics_lscond( &
      & Vars_a%xyz_Temp      , & ! (inout) 温度 (整数) (t+Δt)
      & Vars_a%xyz_Qvap      , & ! (inout) 比湿 (整数) (t+Δt)
      & xy_LscRain           , & ! (out) 降水量
      & xyz_DLscTempDt       , & ! (out) 温度変化率
      & xyz_DLscQvapDt       , & ! (out) 比湿変化率
      & xyz_Press            , & ! (in) 圧力 (整数)
      & xyr_Press            , & ! (in) 圧力 (半整数)
      & 2.0d0*DelTime  )         ! (in) 2Δt

    !----------------------------------------------------------------
    !  7. 負の水蒸気除去(2)
    !----------------------------------------------------------------
    call physics_negq( &
      & Vars_a%xyz_Qvap      , & ! (inout) 比湿 (整数) (t+Δt)
      & xyz_DNegQvap1Dt       , & ! (inout) 比湿変化 (整数)
      & xyr_Press            , & ! (in) 圧力 (半整数)
      & 2.0d0*DelTime  )         ! (in) 2Δt

    !----------------------------------------------------------------
    !  8. 温度半整数 sigma 補間, 気圧と高度の算出 (2)
    !----------------------------------------------------------------

    !----- Ps の計算しなおし -----
!   (2007-7-10 石渡)
!    試しにやめる.
!    そもそも, agcm5 ではこんなことやってないよね?
!    それに式正しいか???
!    do k = 1, km
!       Vars_a%xy_Ps(:,:) = Vars_a%xy_Ps(:,:) &
!            & + ( xyz_DLscQvapDt(:,:,k) &
!            &     + xyz_DCumulusQvapDt(:,:,k) &
!            &     + xyz_DNegQvap1Dt(:,:,k)     ) &
!            &  * ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) &
!            &  * 2.0d0 * DelTime
!    end do

    !----- 温度半整数 sigma 補間 -----
    call physics_interpolate_temp( &
      & xyr_Temp              , & ! (out) 温度 (半整数)
      & Vars_a%xyz_Temp       , & ! (in) 温度 (整数) (t+Δt)
      & Dims%z_Sigma%a_Dim    , & ! (in) σレベル(整数)座標
      & Dims%r_Sigma%a_Dim    )   ! (in) σレベル(半整数)座標

    !----- 気圧と高度の算出 -----
    call physics_interpolate_geopot( &
      & xyz_Press            , & ! (out) 圧力 (整数)
      & xyr_Press            , & ! (out) 圧力 (半整数)
      & xyz_GeoPot           , & ! (out) ジオポテンシャル(整数)
      & xyr_GeoPot           , & ! (out) ジオポテンシャル(半整数)
      & Vars_a%xy_Ps         , & ! (in) 地表面気圧 (t+Δt)
      & Vars_a%xyz_Temp      , & ! (in) 温度 (整数) (t+Δt)
      & xyr_Temp             , & ! (in) 温度 (半整数)
      & Dims%z_Sigma%a_Dim   , & ! (in) σレベル(整数)座標
      & Dims%r_Sigma%a_Dim   )   ! (in) σレベル(半整数)座標

    !----------------------------------------------------------------
    !  9. 陰解配列初期化
    !----------------------------------------------------------------
    call physics_implicit_init( &
      & xyr_VelLonFlux       , & ! (out) 速度経度成分フラックス
      & xyr_VelLatFlux       , & ! (out) 速度緯度成分フラックス
      & xyr_TempFlux         , & ! (out) 温度フラックス
      & xyr_QvapFlux         , & ! (out) 比湿フラックス
      & xyzo_VelMatrix       , & ! (out) 速度陰解行列
      & xyzo_TempMatrix      , & ! (out) 温度陰解行列
      & xyzo_QvapMatrix      , & ! (out) 比湿陰解行列
      & xyr_Press            , & ! (in) 気圧 (半整数)
      & 2.0d0*DelTime        , & ! (in) ２Δt
      & xy_SurfHeatCapacity  , & ! (in) 地表熱容量
      & xy_SurfCondition       ) ! (in) 地表状態

    !----------------------------------------------------------------
    !  10. 放射 flux
    !----------------------------------------------------------------
    call physics_radiation_main( &
      & xyr_RadLFlux              , & ! (out) 長波フラックス
      & xyo_SurfRadLMatrix        , & ! (out) Ｔ陰解行列：地表
      & xyro_DelRadLFlux          , & ! (out) 長波地表温度変化
      & xyr_RadSFlux              , & ! (out) 日射フラックス
      & Vars_a%xyz_Temp           , & ! (in) 温度
      & xyr_Temp                  , & ! (in) 温度(半整数)
      & xy_SurfTemp               , & ! (in) 地表面温度
      & Vars_a%xyz_Qvap           , & ! (in) 比湿
      & xyr_Press                 , & ! (in) 圧力
      & Dims%x_Lon%a_Dim          , & ! (in) 緯度
      & Dims%y_Lat%a_Dim          , & ! (in) 経度
      & xy_SurfAlbedo             )   ! (in) 地表アルベド

    !----------------------------------------------------------------
    !  11. 鉛直拡散 flux
    !----------------------------------------------------------------
    call physics_verdiff_main( &
      & xyr_VelLonFlux       , & ! (out) 速度経度成分フラックス
      & xyr_VelLatFlux       , & ! (out) 速度緯度成分フラックス
      & xyr_TempFlux         , & ! (out) 温度フラックス
      & xyr_QvapFlux         , & ! (out) 比湿フラックス
      & xyzo_VelMatrix       , & ! (inout) 速度陰解行列
      & xyzo_TempMatrix      , & ! (inout) 温度陰解行列
      & xyzo_QvapMatrix      , & ! (inout) 比湿陰解行列
      & Vars_a%xyz_VelLon    , & ! (in) 速度経度成分
      & Vars_a%xyz_VelLat    , & ! (in) 速度緯度成分
      & Vars_a%xyz_Temp      , & ! (in) 温度 (整数)
      & xyr_Temp             , & ! (in) 温度 (半整数)
      & Vars_a%xyz_Qvap      , & ! (in) 比湿 (整数)
      & xyz_Press            , & ! (in) 気圧 (整数)
      & xyr_Press            , & ! (in) 気圧 (半整数)
      & xyz_GeoPot           , & ! (in) 高度 (整数)
      & xyr_GeoPot             ) ! (in) 高度 (半整数)

    !----------------------------------------------------------------
    !  12. 地表 flux
    !----------------------------------------------------------------
    call physics_surface_main( &
      & xyr_VelLonFlux       , & ! (inout) 速度経度成分フラックス
      & xyr_VelLatFlux       , & ! (inout) 速度緯度成分フラックス
      & xyr_TempFlux         , & ! (inout) 温度フラックス
      & xyr_QvapFlux         , & ! (inout) 比湿フラックス
      & xy_SurfVelMatrix     , & ! (out) 速度陰解行列: 地表
      & xyoo_SurfTempMatrix  , & ! (out) 温度陰解行列: 地表
      & xyoo_SurfQvapMatrix  , & ! (out) 比湿陰解行列: 地表
      & Vars_a%xyz_VelLon    , & ! (in) 速度経度成分
      & Vars_a%xyz_VelLat    , & ! (in) 速度緯度成分
      & Vars_a%xyz_Temp      , & ! (in) 温度 (整数)
      & xyr_Temp             , & ! (in) 温度 (半整数)
      & xy_SurfTemp          , & ! (in) 地表面温度
      & Vars_a%xyz_Qvap      , & ! (in) 比湿 (整数)
      & xyz_Press            , & ! (in) 気圧 (整数)
      & xyr_Press            , & ! (in) 気圧 (半整数)
      & xyz_GeoPot           , & ! (in) 高度 (整数)
      & xy_SurfHumidCoeff    , & ! (in) 地表湿潤度
      & xy_SurfRoughLength   , & ! (in) 地表粗度長
      & xy_SurfCondition       ) ! (in) 地表状態

    !----------------------------------------------------------------
    !  13. 時間変化率の計算 (implicit)
    !----------------------------------------------------------------
    call physics_implicit_integrate( &
      & xyz_DVerdiffVelLonDt    , & ! (out) 経度成分 鉛直拡散加速度
      & xyz_DVerdiffVelLatDt    , & ! (out) 緯度成分 鉛直拡散加速度
      & xyz_DVerdiffTempDt      , & ! (out) 鉛直拡散加熱率
      & xyz_DVerdiffSurfTempDt  , & ! (out) 地表面 鉛直拡散加熱率
      & xyz_DVerdiffQvapDt      , & ! (out) 鉛直拡散加湿率
      & xyr_VelLonFlux          , & ! (in) 速度経度成分フラックス
      & xyr_VelLatFlux          , & ! (in) 速度緯度成分フラックス
      & xyr_TempFlux            , & ! (in) 温度フラックス
      & xyr_RadSFlux(:,:,1)     , & ! (in) 日射フラックス
      & xyr_RadLFlux(:,:,1)     , & ! (in) 長波フラックス
      & xy_GroundTempFlux       , & ! (in) 地中熱フラックス
      & xyr_QvapFlux            , & ! (in) 比湿フラックス
      & xyzo_VelMatrix          , & ! (in) 速度陰解行列
      & xyzo_TempMatrix         , & ! (in) 温度陰解行列
      & xyzo_QvapMatrix         , & ! (in) 比湿陰解行列
      & xy_SurfVelMatrix        , & ! (in) 速度陰解行列: 地表
      & xyoo_SurfTempMatrix     , & ! (in) 温度陰解行列: 地表
      & xyoo_SurfQvapMatrix     , & ! (in) 比湿陰解行列: 地表
      & xyo_SurfRadLMatrix      , & ! (in) Ｔ陰解行列：地表
      & 2.0d0*DelTime           , & ! (in) ２Δt
      & xy_SurfCondition          ) ! (in) 地表状態

    ! (2006-7-30 石渡) 追加
    ! フラックス補正
    call physics_implicit_fluxcorrection ( &
      & xyr_VelLonFlux, & !(inout)
      & xyr_VelLatFlux, & !(inout)
      & xyr_TempFlux, & !(inout)
      & xyr_QvapFlux, & !(inout)
      & xyz_DVerdiffVelLonDt, & !(in)
      & xyz_DVerdiffVelLatDt, & !(in)
      & xyz_DVerdiffTempDt, & !(in)
      & xyz_DVerdiffSurfTempDt, & !(in)
      & xyz_DVerdiffQvapDt, & !(in)
      & xyzo_VelMatrix, & !(in)
      & xyzo_TempMatrix, & !(in)
      & xyzo_QvapMatrix, & !(in)
      & xy_SurfVelMatrix, & !(in)
      & xyoo_SurfTempMatrix, & !(in)
      & xyoo_SurfQvapMatrix, & !(in)
      & 2.0d0*DelTime  & !(in)
      & )

    !----------------------------------------------------------------
    !  14. 放射による温度変化率
    !----------------------------------------------------------------
    do k = 1, km+1
! (2007-5-17 石渡 変更)
!       xyr_RadLFlux(:,:,k) = xyr_RadLFlux(:,:,k) &
!            &   + (xyz_DVerdiffSurfTempDt(:,:) * xyro_DelRadLFlux(:,:,k,0)  &
!            &     + xyz_DVerdiffTempDt(:,:,1) * xyro_DelRadLFlux(:,:,k,1) )&
!            &      * 2.0d0*DelTime
       xyr_RadLFluxCorrected(:,:,k) = xyr_RadLFlux(:,:,k) &
         &     + (xyz_DVerdiffSurfTempDt(:,:) * xyro_DelRadLFlux(:,:,k,0)  &
         &        + xyz_DVerdiffTempDt(:,:,1) * xyro_DelRadLFlux(:,:,k,1) )&
         &      * 2.0d0*DelTime
    end do

! (2007-5-17 石渡 変更)
!    call physics_radiation_deltemp( &
!      &  xyz_DRadLTempDt         , & ! (out) 長波加熱率
!      &  xyz_DRadSTempDt         , & ! (out) 短波加熱率
!      &  xyr_RadLFlux            , & ! (in) 補正後の長波フラックス
!      &  xyr_RadSFlux            , & ! (in) 短波フラックス
!      &  xyr_Press                 )   ! (in) 圧力 (半整数)
    call physics_radiation_deltemp( &
      &  xyz_DRadLTempDt         , & ! (out) 長波加熱率
      &  xyz_DRadSTempDt         , & ! (out) 短波加熱率
      &  xyr_RadLFluxCorrected   , & ! (in) 補正後の長波フラックス
      &  xyr_RadSFlux            , & ! (in) 短波フラックス
      &  xyr_Press                 )   ! (in) 圧力 (半整数)

    !----------------------------------------------------------------
    !  15. 温度変化分の足し込み
    !      (2006-7-28 石渡) フラックス補正してないからここで
    !          まとめてやっている.
    !----------------------------------------------------------------
    ! (2007-5-17 石渡) まだちょっと自信がないけど, 
    ! 2.0d0*DelTime ではなく, DelTime じゃないのか?

! (2007-5-17 石渡) やまだ由バージョン
!    Vars_a%xyz_Temp = Vars_a%xyz_Temp               &
!         & +  ( xyz_DRadLTempDt + xyz_DRadSTempDt ) * 2.0d0*DelTime
!
!    Vars_a%xyz_Temp = Vars_a%xyz_Temp               &
!         & +  ( xyz_DVerdiffTempDt ) * 2.0d0* DelTime
!
!    Vars_a%xyz_Qvap = Vars_a%xyz_Qvap               &
!         & +  ( xyz_DVerdiffQvapDt ) * 2.0d0* DelTime
!
!    Vars_a%xyz_VelLon = Vars_a%xyz_VelLon               &
!         & +  ( xyz_DVerdiffVelLonDt ) * 2.0d0* DelTime
!    
!    Vars_a%xyz_VelLat = Vars_a%xyz_VelLat               &
!         & +  ( xyz_DVerdiffVelLatDt ) * 2.0d0* DelTime

! (2007-5-17 石渡) 私は次が正しいと思う.
    Vars_a%xyz_Temp = Vars_a%xyz_Temp               &
      & +  ( xyz_DRadLTempDt + xyz_DRadSTempDt ) * DelTime

    Vars_a%xyz_Temp = Vars_a%xyz_Temp               &
      & +  ( xyz_DVerdiffTempDt ) * DelTime

    Vars_a%xyz_Qvap = Vars_a%xyz_Qvap               &
      & +  ( xyz_DVerdiffQvapDt ) * DelTime

    Vars_a%xyz_VelLon = Vars_a%xyz_VelLon               &
      & +  ( xyz_DVerdiffVelLonDt ) * DelTime
    
    Vars_a%xyz_VelLat = Vars_a%xyz_VelLat               &
      & +  ( xyz_DVerdiffVelLatDt ) * DelTime

    ! (2006-7-28 石渡) 追加.
    call physics_integrate_surftemp( & 
      &  xy_SurfTemp, & !(inout)
      &  xyr_RadLFlux, & !(inout)
      &  xyro_DelRadLFlux, & !(inout)
      &  xyz_DVerdiffSurfTempDt, & !(in)
      &  DelTime & !(in)
      &  )

    !----------------------------------------------------------------
    !  16. 温度半整数 sigma 補間, 気圧と高度の算出 (3)
    !----------------------------------------------------------------

    !----- Ps の計算しなおし -----
!   (2007-7-10 石渡)
!    試しにやめる.
!    そもそも, agcm5 ではこんなことやってないよね?
!    それに式正しいか???
!    do k = 1, km
!       Vars_a%xy_Ps(:,:) = Vars_a%xy_Ps(:,:) &
!            & +  xyz_DVerdiffQvapDt(:,:,k)  &
!            &   * ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) &
!            &   * 2.0d0 * DelTime
!    end do

    !----- 温度半整数 sigma 補間 -----
    call physics_interpolate_temp( &
      & xyr_Temp              , & ! (out) 温度 (半整数)
      & Vars_a%xyz_Temp       , & ! (in) 温度 (整数) (t+Δt)
      & Dims%z_Sigma%a_Dim    , & ! (in) σレベル(整数)座標
      & Dims%r_Sigma%a_Dim    )   ! (in) σレベル(半整数)座標

    !----- 気圧と高度の算出 -----
    call physics_interpolate_geopot( &
      & xyz_Press            , & ! (out) 圧力 (整数)
      & xyr_Press            , & ! (out) 圧力 (半整数)
      & xyz_GeoPot           , & ! (out) ジオポテンシャル(整数)
      & xyr_GeoPot           , & ! (out) ジオポテンシャル(半整数)
      & Vars_a%xy_Ps         , & ! (in) 地表面気圧 (t+Δt)
      & Vars_a%xyz_Temp      , & ! (in) 温度 (整数) (t+Δt)
      & xyr_Temp             , & ! (in) 温度 (半整数)
      & Dims%z_Sigma%a_Dim   , & ! (in) σレベル(整数)座標
      & Dims%r_Sigma%a_Dim   )   ! (in) σレベル(半整数)座標

    !----------------------------------------------------------------
    !  17. 乾燥対流調節
    !----------------------------------------------------------------
    call physics_dryadjust( &
      & Vars_a%xyz_Temp       , & ! (inout) 温度
      & xyz_DDryTempDt        , & ! (out) 温度変化率
      & xyz_Press             , & ! (in) 気圧 (整数)
      & xyr_Press             , & ! (in) 気圧 (半整数)
      & 2.0d0*DelTime  )          ! (in) 2Δt

    ! 鉛直フィルター 
    !   オン/オフ の切り換えは use 文を変更しておこなう.
    call phys_vfilter( & 
      & Vars_a%xyz_Temp   , & !(inout)
      & Vars_a%xyz_VelLon , & !(inout)
      & Vars_a%xyz_VelLat , & !(inout)
      & xyr_Temp   , & !(in) 
      & xyr_Press  , & !(in)
      & 2.0d0*DelTime )   !(in)

    !----------------------------------------------------------------
    !  18. 負の水蒸気除去(3)
    !----------------------------------------------------------------
    call physics_negq( &
      & Vars_a%xyz_Qvap      , & ! (inout) 比湿 (整数) (t+Δt)
      & xyz_DNegQvap2Dt       , & ! (inout) 比湿変化 (整数)
      & xyr_Press            , & ! (in) 圧力 (半整数)
      & 2.0d0*DelTime  )         ! (in) 2Δt

    !----------------------------------------------------------------
    !  19. Ps の計算しなおし 
    !----------------------------------------------------------------
!   (2007-7-10 石渡)
!    試しにやめる.
!    そもそも, agcm5 ではこんなことやってないよね?
!    それに式正しいか???
!    do k = 1, km
!       Vars_a%xy_Ps(:,:) = Vars_a%xy_Ps(:,:) &
!            & + xyz_DNegQvap2Dt(:,:,k) &
!            &  * ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) &
!            &  * 2.0d0 * DelTime
!    end do

    ! 凝結・降水による表面気圧変化
    ! 使用する/しないの切り換えには use 文を書き換えること.
!    call physics_ps_correction( &
!      & Vars_a%xyz_Qvap   , & !(inout)
!      & Vars_a%xy_Ps, & !(inout)
!      & xyr_QvapFlux , & !(in)
!      & xyz_DCumulusQvapDt, & !(in)
!      & xyr_Press  , &  !(in)
!      & 2.0d0 * DelTime         ) !(in) 


    !----------------------------------------------------------------
    !  20. 変数出力
    !----------------------------------------------------------------
    call io_gt4_out_Put(  'GeoPot', real(xyz_GeoPot(:,:,:), DBKIND)  )
    call io_gt4_out_Put(  'Press', real(xyz_Press(:,:,:), DBKIND)  )
    call io_gt4_out_Put(  'SurfTemp', real(xy_SurfTemp(:,:), DBKIND)  )

    call io_gt4_out_Put(  'DNegQvapDt',  &
         & real( (xyz_DNegQvap1Dt(:,:,:) + xyz_DNegQvap2Dt(:,:,:)), DBKIND)  )
    call io_gt4_out_Put(  'DLscTempDt',  &
         &                     real(xyz_DLscTempDt(:,:,:), DBKIND)  )
    call io_gt4_out_Put(  'DLscQvapDt',  &
         &                     real(xyz_DLscQvapDt(:,:,:), DBKIND)  )
    call io_gt4_out_Put(  'LscRain', real(xy_LscRain(:,:), DBKIND)  )

    call io_gt4_out_Put(  'DCumulusTempDt',  &
         &                     real(xyz_DCumulusTempDt(:,:,:), DBKIND)  )

    call io_gt4_out_Put(  'DCumulusQvapDt',  &
         &                     real(xyz_DCumulusQvapDt(:,:,:), DBKIND)  )

    call io_gt4_out_Put(  'CumulusRain', real(xy_CumulusRain(:,:), DBKIND)  )

    call io_gt4_out_Put(  'DDryTempDt',  &
         &                     real(xyz_DDryTempDt(:,:,:), DBKIND)  )


    ! 2007-05-16 M. Ishiwatari 追加
    call io_gt4_out_Put(  'TempFlux',  &
         &                     real(xyr_TempFlux(:,:,:), DBKIND)  )
    call io_gt4_out_Put(  'LatentEnergyFlux',  &
         &                     real(xyr_QvapFlux(:,:,:), DBKIND)  )
    call io_gt4_out_Put(  'RadSFlux',  &
         &                     real(xyr_RadSFlux(:,:,:), DBKIND)  )
! (2007-5-17 石渡 変更)
!    call io_gt4_out_Put(  'RadLFlux',  &
!         &                     real(xyr_RadLFlux(:,:,:), DBKIND)  )
    call io_gt4_out_Put(  'RadLFlux',  &
         &                     real(xyr_RadLFluxCorrected(:,:,:), DBKIND)  )

    xy_Rain = xy_LscRain + xy_CumulusRain
    call io_gt4_out_Put(  'Rain', real(xy_Rain(:,:), DBKIND)  )

    call io_gt4_out_Put(  'DPsDt', & 
         & real( (Vars_a%xy_Ps(:,:) - xy_Ps_b(:,:)), DBKIND)  )

    call io_gt4_out_Put(  'DRadLTempDt', real(xyz_DRadLTempDt(:,:,:), DBKIND)  )
    call io_gt4_out_Put(  'DRadSTempDt', real(xyz_DRadSTempDt(:,:,:), DBKIND)  )
    call io_gt4_out_Put('DVerdiffVelLonDt',  &
         &                     real(xyz_DVerdiffVelLonDt(:,:,:), DBKIND) )
    call io_gt4_out_Put('DVerdiffVelLatDt', & 
         &                     real(xyz_DVerdiffVelLatDt(:,:,:), DBKIND) )
    call io_gt4_out_Put('DVerdiffTempDt', real(xyz_DVerdiffTempDt(:,:,:), DBKIND) )
    call io_gt4_out_Put('DVerdiffQvapDt', real(xyz_DVerdiffQvapDt(:,:,:), DBKIND) )

    !-------------------------------------------------------------------
    !  21. Generate Vorticity and Divergence from Velocity
    !-------------------------------------------------------------------
    
    Vars_a%xyz_Vor = &
      &  xya_wa( wa_Div_xya_xya( Vars_a%xyz_VelLat,-Vars_a%xyz_VelLon )/R0)
    Vars_a%xyz_Div = &
      & xya_wa( wa_Div_xya_xya( Vars_a%xyz_VelLon, Vars_a%xyz_VelLat ) /R0)
    
!    wz_Psi_a = wa_LaplaInv_wa(  wa_xya( Vars_a%xyz_Vor )  ) * R0**2
!    wz_Chi_a = wa_LaplaInv_wa(  wa_xya( Vars_a%xyz_Div )  ) * R0**2

!    Vars_a%xyz_VelLon = (  xya_GradLon_wa( wz_Chi_a ) &
!         &                - xya_GradLat_wa( wz_Psi_a )  ) / R0

!    Vars_a%xyz_VelLat = (  xya_GradLon_wa( wz_Psi_a ) &
!         &                + xya_GradLat_wa( wz_Chi_a )  ) / R0

!    Vars_a%xyz_Temp = xya_wa( wa_xya(Vars_a%xyz_Temp) )
!    Vars_a%xyz_Qvap = xya_wa( wa_xya(Vars_a%xyz_Qvap) )
!    Vars_a%xy_Ps = xy_w( w_xy(Vars_a%xy_Ps) )


    !----------------------------------------------------------------
    !   終了処理
    !----------------------------------------------------------------
    call EndSub(subname)
    
  end subroutine dcpam_ape_physics

end module dcpam_ape_physics_mod
