!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
! physics_cumulus_adjust.f90 
!
! History
!   2005/09/19 Yamada Yukiko     create
!

module physics_cumulus_adjust_mod

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

  implicit none

  private
  public :: physics_cumulus_adjust

contains

  subroutine physics_cumulus_adjust( & 
       & xyz_Temp          , & ! (inout) 温度
       & xyz_Qvap          , & ! (inout) 比湿
       & xy_CumulusRain    , & ! (out) 降水量
       & xyz_DCumulusTempDt, & ! (out) 温度変化率
       & xyz_DCumulusQvapDt, & ! (out) 比湿変化率
       & xyz_Press         , & ! (in) 気圧 (整数)
       & xyr_Press         , & ! (in) 気圧 (半整数)
       & DelTimePhy         )  ! (in) 2Δt

    !==== Dependency
    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  ,& ! 水蒸気気体定数
         &                   RAir  ,& ! 大気気体定数
         &                   Grav     ! 重力加速度
    use physics_qvapsat_tetens, only: QvapSat
    use dc_trace,      only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none

    !==== Input
    !
    real(DBKIND), intent(in) ::   &
         & xyz_Press(im,jm,km)   ,& ! (in) 気圧 (整数)
         & xyr_Press(im,jm,km+1) ,& ! (in) 気圧 (半整数)
         & DelTimePhy               ! (in) 2Δt

    !==== Output
    !
    real(DBKIND), intent(out) ::  &
         & xy_CumulusRain(im,jm)       ,& ! (out) 降水量
         & xyz_DCumulusTempDt(im,jm,km),& ! (out) 温度変化率
         & xyz_DCumulusQvapDt(im,jm,km)   ! (out) 比湿変化率

    !==== In/Out
    !
    real(DBKIND), intent(inout) ::  &
         & xyz_Temp(im,jm,km)      ,& ! (inout) 温度
         & xyz_Qvap(im,jm,km)         ! (inout) 比湿

    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "physics_cumulus_adjust"
    real(DBKIND)        ::  & 
         & xyz_Qvap_b(im,jm,km)              , & ! 調節前の比湿
         & xyz_Temp_b(im,jm,km)              , & ! 調節前の温度
         & xy_Adjust(im,jm)                  , & ! 今回調節されたか否か
         & xy_Adjust_b(im,jm)                , & ! 前回調節されたか否か
         & xyz_DPressDz(im,jm,km)            , & ! ΔＰ
         ! xyz_DelPress の方が良い (2006-10-28 石渡)
         & xyz_QvapSat(im,jm,km)             , & ! 飽和比湿
         & xyz_DDPressDDPress(im,jm,km)      , & ! ΔＰ/ΔＰ
         & xyz_DPFact(im,jm,km)              , & ! ファクター
         & TempSat                           , & ! 飽和温度
         ! (2006-10-30 石渡) これは Stability じゃないか???? 
         & DelTempSat                        , & ! 調節する飽和温度
         & DelQvap                           , & ! 調節する比湿
         & DelTempUpper                      , & ! 調節する温度 (k)
         & DelTempUnder                      , & ! 調節する温度 (k-1)
         & DQvapSatDTempUpper                , & ! D(飽和比湿)/D(温度) (k)
         & DQvapSatDTempUnder                , & ! D(飽和比湿)/D(温度) (k-1)
         & DHDTempUpper                      , & ! D(h)/D(温度) (k)
         & DHDTempUnder                      , & ! D(h)/D(温度) (k-1)
         & Adjust                ! 今回全領域において一度でも調節されたか否か

    real(DBKIND), parameter     :: CrtlRH = 0.990d0   ! 臨界相対湿度
    integer(INTKIND), parameter :: IterationMax = 10 ! イテレーション回数
    
    real(DBKIND) ::  TempSatMax(IterationMax)   ! 不安定の許容誤差
    data TempSatMax  / 0.01  , 0.02  , 0.02  , 0.05  , 0.05  , &
         &             0.10  , 0.10  , 0.20  , 0.20  , 0.40  /

    ! 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_CumulusRain     = 0.0d0
    xyz_DCumulusTempDt = 0.0d0
    xyz_DCumulusQvapDt = 0.0d0
    
    !----- ファクター計算 -----    
    do k = 1, km
       xyz_DPressDz(:,:,k) = xyr_Press(:,:,k) - xyr_Press(:,:,k+1)
    end do

    ! クラウジウスクラペイロンの式より飽和比湿計算
!    xyz_QvapSat = EpsV * ES0  &
!         & *  EXP( EL / RVap * ( 1./273. - 1./xyz_Temp ) ) / xyz_Press

    do i = 1, im
      do j = 1, jm
        do k = 1, km
          xyz_QvapSat(i,j,k) = QvapSat(xyz_Temp(i,j,k), xyz_Press(i,j,k))
        enddo
      enddo
    enddo

    do k = 2, km
       xyz_DDPressDDPress(:,:,k) = xyz_DPressDz(:,:,k) / xyz_DPressDz(:,:,k-1)
       xyz_DPFact(:,:,k) = RAir / Cp &
            & * ( xyz_Press(:,:,k-1) - xyz_Press(:,:,k) ) &
            & / xyr_Press(:,:,k) / 2.
    end do

    !----- 調節 ------    
    xy_Adjust_b = 1.0d0

    !----- 3.1 イテレーション ------    
    do Iteration = 1, IterationMax
       xy_Adjust = 0.0d0

       do k = 2, km
          do i = 1, im
             do j = 1, jm 
                if ( xy_Adjust_b(i,j) .GT. 0.9 ) then
                   
                   TempSat = xyz_Temp(i,j,k-1) - xyz_Temp(i,j,k) &
                        &   + ( xyz_QvapSat(i,j,k-1) - xyz_QvapSat(i,j,k) ) * EL/Cp &
                        &   - xyz_DPFact(i,j,k) &
                        &    * ( xyz_Temp(i,j,k-1) + xyz_Temp(i,j,k) )

                   ! 不安定であるならば
                   if ( TempSat .GT. TempSatMax(Iteration) ) then
                      ! かつ, 飽和しているならば
                      if (   ( xyz_Qvap(i,j,k)/xyz_QvapSat(i,j,k) .GE. CrtlRH ) &
                           & .AND. &
                           & ( xyz_Qvap(i,j,k-1)/xyz_QvapSat(i,j,k-1) .GE. CrtlRH ) &
                           &    ) then
                         
                         DelQvap = xyz_DPressDz(i,j,k-1) &
                              &   * (xyz_Qvap(i,j,k-1) - xyz_QvapSat(i,j,k-1) ) &
                              & +  xyz_DPressDz(i,j,k) &
                              &   * (xyz_Qvap(i,j,k) - xyz_QvapSat(i,j,k) ) 
                         
                         DQvapSatDTempUpper = EL * xyz_QvapSat(i,j,k) &
                              & / ( RVap * xyz_Temp(i,j,k) * xyz_Temp(i,j,k) )
                         
                         DQvapSatDTempUnder = EL * xyz_QvapSat(i,j,k-1) &
                              & / ( RVap * xyz_Temp(i,j,k-1) * xyz_Temp(i,j,k-1) )

                         DHDTempUpper = 1. + EL/Cp * DQvapSatDTempUpper
                         DHDTempUnder = 1. + EL/Cp * DQvapSatDTempUnder

                         DelTempSat = TempSat &
                              & + ( 1. - xyz_DPFact(i,j,k) / DHDTempUnder ) &
                              & * EL/Cp * DelQvap / xyz_DPressDz(i,j,k-1)

                         ! 温度の調節
                         DelTempUpper = DelTempSat / ( & 
                              & ( 1. + xyz_DDPressDDPress(i,j,k)  ) * DHDTempUpper &
                              & + xyz_DPFact(i,j,k) &
                              & * ( 1. - xyz_DDPressDDPress(i,j,k)  &
                              &         * DHDTempUpper / DHDTempUnder ) )

                         DelTempUnder = &
                              & - DHDTempUpper / DHDTempUnder &
                              &  * xyz_DDPressDDPress(i,j,k) * DelTempUpper &
                              & + EL/Cp * DelQvap & 
                              &  / ( xyz_DPressDz(i,j,k-1) * DHDTempUnder )

                         xyz_Temp(i,j,k)   = xyz_Temp(i,j,k) + DelTempUpper
                         xyz_Temp(i,j,k-1) = xyz_Temp(i,j,k-1) + DelTempUnder

                         ! 比湿の調節
                         xyz_Qvap(i,j,k)   = xyz_QvapSat(i,j,k) &
                              & + DQvapSatDTempUpper * DelTempUpper
                         xyz_Qvap(i,j,k-1) = xyz_QvapSat(i,j,k-1) &
                              & + DQvapSatDTempUnder * DelTempUnder
                         xyz_QvapSat(i,j,k)   = xyz_Qvap(i,j,k) 
                         xyz_QvapSat(i,j,k-1) = xyz_Qvap(i,j,k-1) 

                         ! 調節したか否か
                         xy_Adjust(i,j) = 1.
                      end if

                   end if

                end if
             end do
          end do
       end do

       Adjust = 0. 
       do i = 1, im
          do j = 1, jm
             xy_Adjust_b(i,j) = xy_Adjust(i,j)
             Adjust           = Adjust +  xy_Adjust(i,j)
          end do
       end do
       
       if ( Adjust .LT. 1. ) exit
       
    end do

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

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

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

             ! 降水量
             xy_CumulusRain(i,j) = xy_CumulusRain(i,j)  &
                  &  + ( xyz_Qvap_b(i,j,k) - xyz_Qvap(i,j,k) ) &
                  & * xyz_DPressDz(i,j,k) / Grav * EL  / DelTimePhy

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

  end subroutine physics_cumulus_adjust

end module physics_cumulus_adjust_mod




