!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2006. All rights reserved.
!---------------------------------------------------------------------
!= Module MoistBuoyancy
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: moistbuoyancy.f90,v 1.1.1.1 2006/04/25 03:43:58 deepconv Exp $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview
!
!
!== Error Handling
!
!== Bugs
!
!== Note
!
!
!== Future Plans
!
!

module MoistBuoyancy
  !
  !
  !

  !$B%b%8%e!<%kFI$_9~$_(B 
  use gridset, only:  SpcNum,            &! $B2=3X<o$N?t(B
    &                 DimXMin,           &! x $BJ}8~$NG[Ns$N2<8B(B
    &                 DimXMax,           &! x $BJ}8~$NG[Ns$N>e8B(B
    &                 DimZMin,           &! z $BJ}8~$NG[Ns$N2<8B(B
    &                 DimZMax,           &! z $BJ}8~$NG[Ns$N>e8B(B
    &                 DelX,              &! x $BJ}8~$N3J;RE@4V3V(B
    &                 DelZ                ! z $BJ}8~$N3J;RE@4V3V(B
  use basicset, only: CpDry,             &!$B4%Ag@.J,$NHfG.(B
    &                 MolWtDry,          &!$B4%Ag@.J,$NJ,;RNL(B
    &                 Grav,              &!$B=ENO2CB.EY(B
    &                 SpcWetID,          &!$B6E=L@.J,$N(B ID
    &                 SpcWetSymbol,      &!$B6E=L@.J,$NL>A0(B
    &                 MolWtWet,          &!$B6E=L@.J,$NJ,;RNL(B
    &                 xz_PotTempBasicZ,  &!$B4pK\>l$N290L(B
    &                 xz_ExnerBasicZ,    &!$B%(%/%9%J!<4X?t$N4pK\>l(B
    &                 xza_MixRtBasicZ      !$B4pK\>l$N:.9gHf(B
  use ChemCalc, only: xz_LatentHeatPerMass       !$B@xG.(B
  use ChemData, only: ChemData_OneSpcID
  use average,  only: xr_avr_xz, xz_avr_xr
  use differentiate_center4, only: xr_dz_xz

  !$B0EL[$N7?@k8@6X;_(B
  implicit none
  
  !$BB0@-$N;XDj(B
  private
  
  !$B4X?t$r(B public $B$K@_Dj(B
  public MoistBuoy_Init
  public xz_MoistBuoyKm
  public xr_MoistBuoy

  !$BJQ?t$NDj5A(B
  integer     :: LoopNum      = 0
  integer     :: GasNum(10)   = 0
  integer     :: CloudNum(10) = 0
  real(8)     :: Cm           = 2.0d-1
  real(8)     :: MixLen       = 0.0d0
  
  save GasNum, CloudNum, LoopNum
  save Cm, MixLen

contains

!!!------------------------------------------------------------------!!!
  subroutine MoistBuoy_Init( )

    !$B0EL[$N7?@k8@6X;_(B
    implicit none

    !$BJQ?tDj5A(B
    integer                  :: s
    integer                  :: n1

    !-----------------------------------------------------------
    ! $B:.9g5wN%(B
    !-----------------------------------------------------------
    MixLen = sqrt(DelX * DelZ) 
    
    !-----------------------------------------------------------
    ! $B1@N3$H5$BN$N(B ID $B$NAH$r:n$k(B
    !-----------------------------------------------------------
    !$B2=3X<o$NCf$+$i1@N3$r:n$k$b$N$rA*$S(B, $B$=$NG[NsE:$(;z$HJ,;RNL$rJ]4I(B.
    LoopNum = 0
    SelectCloud: do s = 1, SpcNum
      
!      !'Cloud' $B$H$$$&J8;zNs$,4^$^$l$k$b$N$N8D?t$r?t$($k(B
!      n1 = index(SpcWetSymbol(s), '-Cloud' )
!      if (n1 /= 0) then
!        LoopNum          = LoopNum + 1
!        CloudNum(LoopNum)= s
!        GasNum(LoopNum)  = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID(SpcWetSymbol(s)(1:n1-3) // '-g'))
!      end if

      !'-g' $B$H$$$&J8;zNs$,4^$^$l$k$b$N$N8D?t$r?t$($k(B
      n1 = index(SpcWetSymbol(s), '-g' )
      if (n1 /= 0) then
        LoopNum          = LoopNum + 1
!        CloudNum(LoopNum)= s
        GasNum(LoopNum)  = s
      end if
    end do SelectCloud
    
    !-----------------------------------------------------------
    ! $B3NG'(B
    !-----------------------------------------------------------
    write(*,*) "MoistAdjust_Init, LoopNum:  ", LoopNum
!    write(*,*) "MoistAdjust_Init, CloudNum: ", CloudNum
    write(*,*) "MoistAdjust_Init, GasNum:   ", GasNum    

  end subroutine MoistBuoy_Init


!!!------------------------------------------------------------------!!!
  function xz_MoistBuoyKm(xz_Exner,  xz_PotTemp, xza_MixRt)
    !
    !$BMpN.%(%M%k%.!<$NIbNO9`$r7W;;(B
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none
    
    !$BJQ?tDj5A(B
    real(8), intent(in)  :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
                                               !$B%(%/%9%J!<4X?t$N>qMp(B
    real(8), intent(in)  :: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
                                               !$B290L>qMp(B
    real(8), intent(in)  :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                               !$B6E=L@.J,$N:.9gHf(B
    real(8)              :: xz_MoistBuoyKm(DimXMin:DimXMax, DimZMin:DimZMax)
                                               !$BIbNO9`(B
    real(8)              :: xz_EquivPotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
                                               !$BAjEv290L(B
    real(8)              :: xz_Temp(DimXMin:DimXMax, DimZMin:DimZMax)
                                               !$B29EY(B
    real(8)              :: xz_LatentHeatPerMassAll(DimXMin:DimXMax, DimZMin:DimZMax, LoopNum)
                                               !$B@xG.$NAmOB(B
    real(8)              :: xza_MixRtPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax, LoopNum)
                                               !$B:.9gHf(B / $BJ,;RNL(B
    real(8)              :: xza_MixRtBasicZPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax, LoopNum)
                                               !$B4pK\>l$N:.9gHf(B / $BJ,;RNL(B
    integer              :: s

    !$B=i4|2=(B
    xz_Temp                = 0.0d0
    xz_LatentHeatPerMassAll       = 0.0d0
    xza_MixRtPerMolWt       = 0.0d0
    xza_MixRtBasicZPerMolWt = 0.0d0

    !$B:n6HG[Ns$N=i4|2=(B. $B5$BN$N$_MxMQ(B
    xz_Temp = (xz_Exner + xz_ExnerBasicZ) * (xz_PotTemp + xz_PotTempBasicZ)
    do s = 1, LoopNum
      xza_MixRtPerMolWt(:,:,GasNum(s)) = xza_MixRt(:,:,GasNum(s)) / MolWtWet(GasNum(s))

      xza_MixRtBasicZPerMolWt(:,:,GasNum(s)) = &
        & xza_MixRtBasicZ(:,:,GasNum(s)) / MolWtWet(GasNum(s))

      xz_LatentHeatPerMassAll(:,:,GasNum(s)) = &
        &  (xza_MixRt(:,:,GasNum(s)) + xza_MixRtBasicZ(:,:,GasNum(s))) &
        &   * xz_LatentHeatPerMass(SpcWetID(CloudNum(s)), xz_Temp)
    end do
  
    !$BAjEv290L$N7W;;(B
    xz_EquivPotTemp = 0.0d0 
    xz_EquivPotTemp =                                           & 
      &   sum(xz_LatentHeatPerMassAll, 3) / (CpDry * xz_ExnerBasicZ)   &
      & + xz_PotTempBasicZ                                      &
      &   * (                                                   &
      &       + sum(xza_MixRtPerMolWt, 3)                        &
      &           / ( +  1.0d0 / MolWtDry                       &
      &               + sum(xza_MixRtBasicZPerMolWt, 3) )        &
      &       - sum(xza_MixRt, 3)                                &
      &          / ( 1.0d0 + sum(xza_MixRtBasicZ, 3))            &
      &      )
    
    !$BIbNO9`$N7W;;(B
    xz_MoistBuoyKm = 0.0d0
    xz_MoistBuoyKm = &
      &  - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) &
      &       * xz_avr_xr( xr_dz_xz( xz_EquivPotTemp ) ) &
      &       / ( 2.0d0 * xz_PotTempBasicZ )
    
  end function xz_MoistBuoyKm
  

!!!------------------------------------------------------------------------!!!
  function xr_MoistBuoy(xza_MixRt)
    !
    ! $B1tD>J}8~$N1?F0J}Dx<0$K8=$l$kIbNO9`$N$&$A(B, 
    ! $BJ,;RNL$N8z2L$@$1$r5a$a$k(B
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none
    
    !$BJQ?tDj5A(B
    real(8), intent(in)  :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                               !$B6E=L@.J,$N:.9gHf(B
    real(8)              :: xr_MoistBuoy(DimXMin:DimXMax, DimZMin:DimZMax)
                                               !$BIbNO9`(B($BJ,;RNL8z2L(B)
    real(8)              :: xr_MoistBuoy2(DimXMin:DimXMax, DimZMin:DimZMax)
                                               !$BIbNO9`(B($BJ,;RNL8z2L(B)
    real(8)              :: xza_MixRtPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax, LoopNum)
                                               !$B:.9gHf(B / $BJ,;RNL(B
    real(8)              :: xza_MixRtBasicZPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax, LoopNum)
                                               !$B4pK\>l$N:.9gHf(B / $BJ,;RNL(B
    integer              :: s
    
    !$B=i4|2=(B
    xr_MoistBuoy           = 0.0d0
    xza_MixRtPerMolWt       = 0.0d0
    xza_MixRtBasicZPerMolWt = 0.0d0

    !$B:n6HG[Ns$N=i4|2=(B. $B5$BN$N$_MxMQ(B
    do s = 1, LoopNum
      xza_MixRtPerMolWt(:,:,GasNum(s)) = xza_MixRt(:,:,GasNum(s)) / MolWtWet(GasNum(s))

      xza_MixRtBasicZPerMolWt(:,:,GasNum(s)) = &
        & xza_MixRtBasicZ(:,:,GasNum(s)) / MolWtWet(GasNum(s))
    end do
    
    !$BIbNO9`$N7W;;(B
    xr_MoistBuoy =                                           &
      & + Grav * xr_avr_xz( sum(xza_MixRtPerMolWt, 3) )        &
      &    / ( 1.0d0 / MolWtDry                               &
      &        + xr_avr_xz( sum(xza_MixRtBasicZPerMolWt, 3) )) &
      & - Grav * xr_avr_xz( sum(xza_MixRt, 3) )                &
      &    / ( 1.0d0 + xr_avr_xz( sum(xza_MixRtBasicZ, 3) ) )

  end function xr_MoistBuoy

  
end module MoistBuoyancy
