!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
! physics_cumulus_adjust_muchwater.f90 
!
! History
!   2006/10/27 ISHIWATARI Masaki
!              imported from AGCM5.3 momoko version(p2cumaG1n.F)

module physics_cumulus_adjust_muchwater_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  , & ! (inout)
    & xyz_DCumulusTempDt, & ! (inout)
    & xyz_DCumulusQvapDt, & ! (inout)
    & xyz_Press, & !(in)
    & xyr_Press,  & !(in)
    &  DelTimePhy  )  !(in)
    !
    != 湿潤対流調節スキーム(水少ない近似を緩和したバージョン)
    !==TODO
    ! * NAMELIST によるパラメータ読み込み

    use type_mod,      only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod,   only: im, jm, km
    use constants_mod, only: Cp    ,& ! 大気定圧比熱 
         &                   EL    ,& ! 水の凝結の潜熱 
         &                   RAir  ,& ! 大気気体定数
         &                   Grav     ! 重力加速度
    use physics_qvapsat_nha92, only: QvapSat, DQvapSatDTemp
    use dc_trace,      only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none

    real(DBKIND), intent(inout) :: xyz_Temp ( im*jm, km ) ! 温度  Ｔ
    real(DBKIND), intent(inout) :: xyz_Qvap ( im*jm, km ) ! 比湿  ｑ
    real(DBKIND), intent(inout) :: xy_CumulusRain ( im*jm ) ! 降水量
    real(DBKIND), intent(inout) :: xyz_DCumulusTempDt ( im*jm,km ) 
                                   ! 凝結調節加熱率
    real(DBKIND), intent(inout) :: xyz_DCumulusQvapDt ( im*jm, km )
                                   ! 凝結調節比湿変化
    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

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

    NAMELIST  /NMMADJ/ TempSatMax, CrtlRH

    real(DBKIND) :: xy_Adjust ( im*jm )     ! 今回調節されたか否か
    real(DBKIND) :: xy_Adjust_b ( im*jm )   ! 前回調節されたか否か
    real(DBKIND) :: xyz_Temp_b  ( im*jm, km )   ! 調節前の温度
    real(DBKIND) :: xyz_Qvap_b  ( im*jm, km )   ! 調節前の比湿
    real(DBKIND) :: xyz_DPressDz  ( im*jm, km )   ! ΔＰ
    real(DBKIND) :: xyz_DDPressDDPress  ( im*jm, km )   ! ファクター ΔＰ/ΔＰ
    ! REAL       DPPMK ( im*jm, km )   ! ファクター
    real(DBKIND) :: xyz_QvapSat  ( im*jm, km )   ! 飽和比湿

    integer(INTKIND) :: ij, k, i, j ! ループ変数
    integer(INTKIND) :: Iteration ! 調節の繰り返し回数
    real(DBKIND) :: ELF, AKAPPA
    real(DBKIND) :: ST
    real(DBKIND) :: QEXE
    real(DBKIND) :: DQvapSatDTempUpper
    real(DBKIND) :: DQvapSatDTempLower
    real(DBKIND) :: DelTempUpper
    real(DBKIND) :: DelTempLower
    real(DBKIND) :: STEXE
    real(DBKIND) :: GammaUpper
    real(DBKIND) :: GammaLower
    real(DBKIND) :: ADJPTS
    real(DBKIND) :: MM, TP, TM, M1, M2, QP, QM, PP1, PP2, B, C
    character(STRING),  parameter:: subname = "physics_cumulus_adjust"

    continue

    !  開始処理
    call BeginSub(subname)

    ! 調節前温度, 比湿のセーブ >
    xyz_Temp_b = xyz_Temp
    xyz_Qvap_b = xyz_Qvap

    ! ファクター計算
    ELF    = EL/Cp
    AKAPPA = RAir / Cp

    do k = 1, km
      xyz_DPressDz ( :,k ) = xyr_Press( :,k ) - xyr_Press( :,k+1 )
    enddo

    do k = 2, km
      xyz_DDPressDDPress ( :,k ) = xyz_DPressDz( :,k ) / xyz_DPressDz( :,k-1 )
    enddo

    do k = 1, km
      do ij = 1, im*jm
        xyz_QvapSat  ( ij,k ) = QvapSat( xyz_Temp( ij,k ), xyz_Press( ij,k ) )
      enddo
    enddo

    ! 3. 調節
    xy_Adjust_b = 0.0d0
    do i = 1, im
      do j = 1, jm
        xy_Adjust_b( im*(j-1)+i ) = 1.  
      enddo
    enddo

    ! 3.1 イテレーション
    do Iteration = 1, IterationMax
      xy_Adjust = 0.0d0
      do k = 2, km
        do ij = 1, im*jm
          if ( xy_Adjust_b( ij ) .GT. 0.9   ) then
            ! (2006-11-13 石渡) xyz_QvapSat はループの先頭で計算し直さないと
            !                   ダメじゃない!?
            MM = 1.0 - 0.5 * xyz_QvapSat(ij,k-1) - 0.5 * xyz_QvapSat(ij,k)
            TP = xyz_Temp(ij,k-1) + xyz_Temp(ij,k)
            TM = xyz_Temp(ij,k-1) - xyz_Temp(ij,k)
            M1 = 1.0 - xyz_QvapSat(ij,k-1)
            M2 = 1.0 - xyz_QvapSat(ij,k  )
            QP = xyz_QvapSat(ij,k-1) + xyz_QvapSat(ij,k)
            QM = xyz_QvapSat(ij,k-1) - xyz_QvapSat(ij,k)
            PP1 = xyz_Press(ij,k-1)/xyr_Press(ij,k)
            PP2 = xyz_Press(ij,k  )/xyr_Press(ij,k)
            ST = AKAPPA/4.0d0 * TP**2 * MM*(PP1*M1 - PP2*M2) &
              & - 0.5d0*MM*TP*TM &
              & + ELF/2.0d0 * (QP*MM*TM - TP*QM)

            if ( ST .GT. TempSatMax(Iteration) ) then ! 不安定の場合
              if ( ( xyz_Qvap(ij,k  ) .GE. CrtlRH*xyz_QvapSat(ij,k) ) &
                &  .AND.( xyz_Qvap(ij,k-1).GE. CrtlRH*xyz_QvapSat(ij,k-1) ) &
                & ) then  ! 飽和
                QEXE   = xyz_DPressDz( ij,k-1 ) &
                  & *( xyz_Qvap( ij,k-1 )-xyz_QvapSat( ij,k-1 ) ) &
                  & + xyz_DPressDz( ij,k   ) &
                  & *( xyz_Qvap( ij,k   )-xyz_QvapSat( ij,k   ) )

                DQvapSatDTempUpper = DQvapSatDTemp( xyz_Temp(ij,k  ), xyz_Press(ij,k  ) )
                DQvapSatDTempLower = DQvapSatDTemp( xyz_Temp(ij,k-1), xyz_Press(ij,k-1) )
                GammaUpper   = 1.0d0 + ELF * DQvapSatDTempUpper
                GammaLower  = 1.0d0 + ELF * DQvapSatDTempLower

                B = - AKAPPA/4.0d0 &
                  &   * ( - TP**2 * MM * PP1 *DQvapSatDTempLower &
                  &       + (PP1*M1 - PP2*M2) &
                  &          * ( - 0.5d0 * TP**2 * DQvapSatDTempLower &
                  &              + 2.0d0 * TP * MM )   ) &
                  & + 0.5d0 * ( MM*(TP + TM) - 0.5d0*TP*TM*DQvapSatDTempLower ) &
                  & - ELF/2.0d0 &
                  &   * (   QP*MM - 0.5d0*QP*TM*DQvapSatDTempLower &
                  &       + MM * TM * DQvapSatDTempLower ) &
                  & + ELF/2.0d0 * ( TP*DQvapSatDTempLower + QM )
                C = - AKAPPA/4.0d0 &
                  &   * (   TP**2 * MM * PP2 * DQvapSatDTempUpper &
                  &       + (PP1*M1 - PP2*M2) &
                  &         * ( - 0.5d0 * TP**2 * DQvapSatDTempUpper + 2.0d0*TP*MM)  ) &
                  & + 0.5d0 &
                  &   * (   MM*(-TP+TM) &
                  &       - 0.5d0*TP*TM*DQvapSatDTempUpper ) &
                  & - ELF/2.0d0 &
                  &   * ( - QP*MM - 0.5d0*QP*TM*DQvapSatDTempUpper &
                  &       + MM * TM * DQvapSatDTempUpper ) &
                  & + ELF/2.0d0 * ( - TP*DQvapSatDTempUpper + QM )

                STEXE =   ST &
                  & - B * ELF * QEXE &
                  &   / ( xyz_DPressDz( ij,k-1 ) * GammaLower )

                ! 温度の調節
                DelTempUpper = STEXE &
                 &   /( C + B * ( - xyz_DDPressDDPress( ij,k )*GammaUpper/GammaLower ) )

                DelTempLower &
                  & = - GammaUpper/GammaLower*xyz_DDPressDDPress( ij,k )*DelTempUpper &
                  &   + ELF * QEXE / ( xyz_DPressDz( ij,k-1 ) * GammaLower )

                xyz_Temp( ij,k   ) = xyz_Temp( ij,k   ) + DelTempUpper
                xyz_Temp( ij,k-1 ) = xyz_Temp( ij,k-1 ) + DelTempLower

                ! 比湿の調節
                xyz_Qvap(ij,k  ) = xyz_QvapSat( ij,k ) + DQvapSatDTempUpper * DelTempUpper
                xyz_Qvap(ij,k-1) = xyz_QvapSat( ij,k-1 ) + DQvapSatDTempLower * DelTempLower

                xyz_QvapSat( ij,k   ) = xyz_Qvap( ij,k   )
                xyz_QvapSat( ij,k-1 ) = xyz_Qvap( ij,k-1 )

                xy_Adjust( ij ) = 1.0d0
              endif
            endif
          endif
        enddo
      enddo
      ADJPTS = 0.0d0  
      do ij = 1, im*jm
        xy_Adjust_b( ij ) = xy_Adjust( ij )
        ADJPTS = ADJPTS + xy_Adjust( ij )
      enddo
      if ( ADJPTS .LT. 1.0d0 ) exit
    enddo

    ! 降水量, 温度変化
    do k = 1, km
      do ij = 1, im*jm
        xyz_DCumulusTempDt ( ij,k ) = xyz_DCumulusTempDt ( ij,k ) &
           & + ( xyz_Temp ( ij,k ) - xyz_Temp_b( ij,k ) ) / DelTimePhy
        xyz_DCumulusQvapDt  ( ij,k ) = xyz_DCumulusQvapDt ( ij,k ) &
           & + ( xyz_Qvap ( ij,k ) - xyz_Qvap_b( ij,k ) ) / DelTimePhy
        xy_CumulusRain   ( ij   ) = xy_CumulusRain  ( ij ) & 
           & + ( xyz_Qvap_b ( ij,k ) - xyz_Qvap ( ij,k ) ) &
           &   * xyz_DPressDz( ij,k ) / Grav * EL / DelTimePhy
      enddo
    enddo

    ! 終了処理
    call EndSub(subname)

  end subroutine physics_cumulus_adjust

end module physics_cumulus_adjust_muchwater_mod
