!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
! physics_lscond.f90 
! 診断型大規模凝結過程
!
! History
!   2005/09/19 Yamada Yukiko     create
!   2006/11/18 Ishiwatari Masaki modify
!

module physics_lscond_mod

  use type_mod,    only : REKIND, DBKIND, INTKIND, TOKEN, STRING

  implicit none

  private
  public :: physics_lscond

contains

  subroutine physics_lscond( xyz_Temp, xyz_Qvap, xy_LscRain, &
    &               xyz_DLscTempDt, xyz_DLscQvapDt,    &
    &               xyz_Press, xyr_Press, DelTimePhy    )
    !
    !== Large scale condensation
    !=== TODO
    ! * (i,j,k) における飽和比湿の変数を QvapSatIJK にしている.
    !   これは変数命名ルールからすると大変良くない.
    !   何か良い手はないものか???
    !
    use type_mod,      only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod,   only: im, jm, km
    use constants_mod, only: Cp    ,& ! 大気定圧比熱 
      &                   EL    ,& ! 水の凝結の潜熱 
      &                   EpsV  ,& ! 水と大気の分子量比ε
      &                   ES0   ,& ! ０℃飽和蒸気圧：Ｐａ
      &                   RVap  ,& ! 水蒸気気体定数
      &                   Grav     ! 重力加速度
    use physics_qvapsat_tetens, only: QvapSat, DQvapSatDTemp
    !use physics_qvapsat_nha92, only: QvapSat, DQvapSatDTemp
    use dc_trace, only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none

    real(DBKIND), intent(in) :: xyz_Press(im,jm,km)   ! 気圧 (整数)
    real(DBKIND), intent(in) :: xyr_Press(im,jm,km+1) ! 気圧 (半整数)
    real(DBKIND), intent(in) :: DelTimePhy            ! 2Δt
    real(DBKIND), intent(out) :: xy_LscRain(im,jm)    ! 降水量
    real(DBKIND), intent(out) :: xyz_DLscTempDt(im,jm,km) ! 温度変化率
    real(DBKIND), intent(out) :: xyz_DLscQvapDt(im,jm,km) ! 比湿変化率
    real(DBKIND), intent(inout) :: xyz_Temp(im,jm,km)     ! 温度
    real(DBKIND), intent(inout) :: xyz_Qvap(im,jm,km)     ! 比湿

    character(STRING),  parameter:: subname = "physics_lscond"
    real(DBKIND) :: xyz_Qvap_b(im,jm,km)             ! 調節前の比湿
    real(DBKIND) :: xyz_Temp_b(im,jm,km)             ! 調節前の温度
    real(DBKIND) :: QvapSatIJK                       ! 飽和比湿 ≒ 飽和混合比
    real(DBKIND) :: DQvapSatDTempIJK                 ! D(飽和比湿)/D(温度)
    real(DBKIND) :: DelTemp, DelQvap 
    real(DBKIND), parameter     :: CrtlRH = 1.0d0   ! 臨界相対湿度
    integer(INTKIND), parameter :: IterationMax = 3 ! イテレーション回数

    ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    integer(INTKIND)    :: i, j, k
    integer(INTKIND)    :: Iteration 

    continue

    !----------------------------------------------------------------
    !   開始処理
    !----------------------------------------------------------------
    call BeginSub(subname)

    !----------------------------------------------------------------
    !   大規模凝結
    !----------------------------------------------------------------

    !----- 調節前 Qvap の保存 -----    
    xyz_Qvap_b  = xyz_Qvap
    xyz_Temp_b  = xyz_Temp

    !----- 変数初期化 -----    
    xy_LscRain     = 0.0d0
    xyz_DLscTempDt = 0.0d0
    xyz_DLscQvapDt = 0.0d0
    
    !----- 調節 ------    
    do k = km, 1, -1
       do i = 1, im
          do j = 1, jm
             
             ! クラウジウスクラペイロンの式より飽和比湿計算
             QvapSatIJK = QvapSat( xyz_Temp(i,j,k), xyz_Press(i,j,k) ) 

             ! 飽和していたら, 温度と比湿の変化の収束計算
             if ( ( xyz_Qvap(i,j,k) / QvapSatIJK ) .GE. CrtlRH ) then

                do Iteration = 1, IterationMax

                   ! 飽和比湿計算
                   QvapSatIJK = QvapSat( xyz_Temp(i,j,k), xyz_Press(i,j,k) ) 
                   DQvapSatDTempIJK &
                     & = DQvapSatDTemp( xyz_Temp(i,j,k), xyz_Press(i,j,k) ) 
                   
                   ! 温度と比湿の変化分を Newton 法で求める
                   DelTemp = EL / Cp * ( xyz_Qvap(i,j,k) - QvapSatIJK ) &
                     &   / ( 1. + EL / Cp * DQvapSatDTempIJK )
                   DelQvap  = DQvapSatDTempIJK * DelTemp 

                   ! 調節
                   xyz_Temp(i,j,k) = xyz_Temp(i,j,k) + DelTemp
                   xyz_Qvap(i,j,k) = QvapSatIJK + DelQvap

                end do

            end if
          end do
       end do
    end do

    !----- 比湿変化率, 温度変化率, 降水量の算出 ----- 
    do k = km, 1, -1
       do i = 1, im
          do j = 1, jm

             ! 比湿変化率
             xyz_DLscQvapDt(i,j,k) = xyz_DLscQvapDt(i,j,k) & 
               &  + ( xyz_Qvap(i,j,k) - xyz_Qvap_b(i,j,k) ) / DelTimePhy

             ! 温度変化率
             xyz_DLscTempDt(i,j,k) = xyz_DLscTempDt(i,j,k) &
               & + ( xyz_Temp(i,j,k) - xyz_Temp_b(i,j,k) ) / DelTimePhy

             ! 降水量
             xy_LscRain(i,j) = xy_LscRain(i,j)  &
               & + ( xyz_Temp(i,j,k) - xyz_Temp_b(i,j,k) ) &
               & * CP / DelTimePhy                         &
               & * ( xyr_Press(i,j,k) - xyr_Press(i,j,k+1) ) /Grav
 
          end do
       end do
    end do
    
    !----------------------------------------------------------------
    !   終了処理
    !----------------------------------------------------------------
    call EndSub(subname)

  end subroutine physics_lscond

end module physics_lscond_mod