| Class | ECCM |
| In: |
moist/eccm.f90
|
| Subroutine : |
subroutine ECCM_Init( )
!暗黙の型宣言禁止
implicit none
!変数定義
integer :: s
integer :: n1
!-----------------------------------------------------------
! 雲粒と気体の ID の組を作る
!-----------------------------------------------------------
!初期化
LoopNum = 0
!化学種の中から雲粒を作るものを選び, その配列添え字と分子量を保管.
SelectCloud: do s = 1, SpcNum
! NH4SH については無視
if ( trim(SpcWetSymbol(s)) == 'NH4SH-s-Cloud' ) then
cycle SelectCloud
end if
!'Cloud' という文字列が含まれるものの個数を数える
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
end do SelectCloud
!-----------------------------------------------------------
! 硫化アンモニウム, およびアンモニアと硫化水素の ID を取得
!-----------------------------------------------------------
NH3Num = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('NH3-g'))
H2SNum = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('H2S-g'))
!-----------------------------------------------------------
! 確認
!-----------------------------------------------------------
if ( LoopNum == 0 ) then
write(*,*) "ECCM: CloudNum = 0, please comment out of MoistAdjust"
! stop
end if
write(*,*) "ECCM_Init, LoopNum: ", LoopNum
write(*,*) "ECCM_Init, CloudNum: ", CloudNum
write(*,*) "ECCM_Init, GasNum: ", GasNum
write(*,*) "ECCM_Init, NH3Num: ", NH3Num
write(*,*) "ECCM_Init, H2SNum: ", H2SNum
end subroutine ECCM_Init
| Subroutine : | |
| a_MolFrIni(1:SpcNum) : | real(8), intent(in) |
| Humidity : | real(8), intent(in) |
| z_Temp(DimZMin:DimZMax) : | real(8), intent(in) |
| z_Press(DimZMin:DimZMax) : | real(8), intent(in) |
| za_MolFr(DimZMin:DimZMax, 1:SpcNum) : | real(8), intent(out) |
与えられた温度に対し, 気塊が断熱的に上昇した時に実現される モル比のプロファイルを求める
subroutine ECCM_MolFr( a_MolFrIni, Humidity, z_Temp, z_Press, za_MolFr )
!
! 与えられた温度に対し, 気塊が断熱的に上昇した時に実現される
! モル比のプロファイルを求める
!
!暗黙の型宣言禁止
implicit none
real(8), intent(in) :: a_MolFrIni(1:SpcNum)
real(8), intent(in) :: Humidity
real(8), intent(in) :: z_Temp(DimZMin:DimZMax)
real(8), intent(in) :: z_Press(DimZMin:DimZMax)
real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum)
real(8) :: DelMolFr
integer :: k, s
!-----------------------------------------------------------
! 配列の初期化
!-----------------------------------------------------------
do s = 1, SpcNum
za_MolFr(:,s) = a_MolFrIni(s)
end do
!-----------------------------------------------------------
! 断熱減率 dT/dz の計算.
!-----------------------------------------------------------
do k = RegZMin, DimZMax
za_MolFr(k,:) = za_MolFr(k-1,:)
!------------------------------------------------------------
!NH4SH 以外の化学種の平衡条件
!------------------------------------------------------------
do s = 1, LoopNum
!モル比を求める
!モル比は前のステップでのモル比を超えることはない
za_MolFr(k,GasNum(s)) = min( za_MolFr(k-1,GasNum(s)), SvapPress( SpcWetID(CloudNum(s)), z_Temp(k) ) * Humidity / z_Press(k) )
end do
!------------------------------------------------------------
!NH4SH の平衡条件
!------------------------------------------------------------
if ( NH3Num /= 0 ) then
!モル比の変化.
!とりあえず NH4SH に対する飽和比は 1.0 とする(手抜き...).
DelMolFr = max ( DelMolFrNH4SH( z_Temp(k), z_Press(k), za_MolFr(k,NH3Num), za_MolFr(k,H2SNum), Humidity ), 0.0d0 )
za_MolFr(k,NH3Num) = za_MolFr(k,NH3Num) - DelMolFr
za_MolFr(k,H2SNum) = za_MolFr(k,H2SNum) - DelMolFr
write(*,*) k, z_Temp(k), za_MolFr(k,NH3Num), za_MolFr(k,H2SNum)
end if
end do
end subroutine ECCM_MolFr
| Subroutine : | |
| xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax) : | real(8), intent(in) |
| xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) : | real(8), intent(in) |
| xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) : | real(8), intent(in) |
| xz_Stab(DimXMin:DimXMax, DimZMin:DimZMax) : | real(8), intent(out) |
| xz_StabTemp(DimXMin:DimXMax, DimZMin:DimZMax) : | real(8), intent(out) |
| xz_StabMolWt(DimXMin:DimXMax, DimZMin:DimZMax) : | real(8), intent(out) |
subroutine ECCM_Stab( xz_PotTemp, xz_Exner, xza_MixRt, xz_Stab, xz_StabTemp, xz_StabMolWt )
use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax, SpcNum !
use basicset,only: MolWtDry, MolWtWet, CpDry, Grav, xz_ExnerBasicZ, xz_PotTempBasicZ, xza_MixRtBasicZ
use storeset,only: StoreStab
use average, only: xz_avr_xr
use differentiate_center2, only: xr_dz_xz
implicit none
real(8), intent(in) :: xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax)
real(8), intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
real(8), intent(in) :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
real(8), intent(out) :: xz_Stab(DimXMin:DimXMax, DimZMin:DimZMax)
real(8), intent(out) :: xz_StabTemp(DimXMin:DimXMax, DimZMin:DimZMax)
real(8), intent(out) :: xz_StabMolWt(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: xza_MolFrAll(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum)
real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: xz_MolWtWet(DimXMin:DimXMax, DimZMin:DimZMax)
integer :: i, k, s
xz_TempAll = (xz_PotTemp + xz_PotTempBasicZ) * (xz_Exner + xz_ExnerBasicZ)
do s = 1, SpcNum
xza_MolFrAll(:,:,s) = (xza_MixRt(:,:,s) + xza_MixRtBasicZ(:,:,s)) * MolWtDry / MolWtWet(s)
end do
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
xz_MolWtWet(i,k) = dot_product( MolWtWet(1:3), xza_MolFrAll(i,k,1:3) )
end do
end do
xz_StabTemp = Grav / xz_TempAll * ( xz_avr_xr( xr_dz_xz( xz_TempAll ) ) + Grav / CpDry )
xz_StabMolWt = - Grav * xz_avr_xr( xr_dz_xz( xz_MolWtWet ) ) / MolWtDry
xz_Stab = Grav / xz_TempAll * ( xz_avr_xr( xr_dz_xz( xz_TempAll ) ) + Grav / CpDry ) - Grav * xz_avr_xr( xr_dz_xz( xz_MolWtWet ) ) / MolWtDry
where (xz_Stab < 1.0d-7)
xz_Stab = 1.0d-7
end where
call StoreStab( xz_Stab, xz_StabTemp, xz_StabMolWt)
end subroutine ECCM_Stab
| Subroutine : | |
| MolFrIni(1:SpcNum) : | real(8), intent(in) |
| z_Temp(DimZMin:DimZMax) : | real(8), intent(out) |
| z_Press(DimZMin:DimZMax) : | real(8), intent(out) |
subroutine ECCM_Temp_Press( MolFrIni, z_Temp, z_Press )
!暗黙の型宣言禁止
implicit none
real(8), intent(in) :: MolFrIni(1:SpcNum)
real(8), intent(out):: z_Temp(DimZMin:DimZMax)
real(8), intent(out):: z_Press(DimZMin:DimZMax)
real(8) :: z_MolFr(DimZMin:DimZMax, 0:SpcNum)
!モル分率
real(8) :: z_DTempDZ(DimZMin:DimZMax)
real(8) :: MolFr(0:SpcNum)
integer :: k
real(8) :: Temp1, Press1, DTempDZ1
real(8) :: Temp2, Press2, DTempDZ2
real(8) :: Temp3, Press3, DTempDZ3
real(8) :: Temp4, Press4, DTempDZ4
!-------------------------------------------------------------
! 配列の初期化
!-------------------------------------------------------------
!初期モル比
z_MolFr = 0.0d0
z_MolFr(RegZMin, 1:SpcNum) = MolFrIni(1:SpcNum)
z_MolFr(RegZMin, 0) = 1.0d0 - sum(MolFrIni)
z_MolFr(RegZMin-1, 1:SpcNum) = MolFrIni(1:SpcNum)
z_MolFr(RegZMin-1, 0) = 1.0d0 - sum(MolFrIni)
!地表面での温度(RegZMin は, 高度 DelZ / 2 に相当)
z_Temp = 1.0d-60
z_Temp(RegZMin) = TempSfc - Grav / CpDry * ( DelZ * 5.0d-1 )
!地表面での圧力(RegZMin は, 高度 DelZ / 2 に相当)
z_Press = 1.0d-60
z_Press(RegZMin) = PressSfc *((TempSfc / z_Temp(RegZMin)) ** ( - CpDryMol / GasRUniv ))
!-----------------------------------------------------------
! 断熱減率 dT/dz の計算.
!-----------------------------------------------------------
DtDz: do k = RegZMin, DimZMax-1
!初期化
z_MolFr(k,:) = z_MolFr(k-1,:)
!==============================================================
!湿潤断熱減率をルンゲクッタ法を用いて計算
Temp1 = z_Temp(k)
Press1 = z_Press(k)
MolFr = z_MolFr(k-1,:)
call ECCM_DTempDZ( Temp1, Press1, z_MolFr(k,:), DTempDZ1 )
Temp2 = Temp1 + DTempDZ1 * DelZ * 5.0d-1
if (Temp2 < 0.0d0) exit DtDz
Press2 = Press1 * ( ( Temp1 / Temp2 ) ** (Grav * MolWtDry / (GasRUniv * DTempDZ1) ) )
call ECCM_DTempDZ( Temp2, Press2, MolFr, DTempDZ2 )
Temp3 = Temp1 + DTempDZ2 * DelZ * 5.0d-1
if (Temp3 < 0.0d0) exit DtDz
Press3 = Press1 * ( ( Temp1 / Temp3 ) ** (Grav * MolWtDry / (GasRUniv * DTempDZ2) ) )
call ECCM_DTempDZ( Temp3, Press3, MolFr, DTempDZ3 )
Temp4 = Temp1 + DTempDZ3 * DelZ
if (Temp4 < 0.0d0) exit DtDz
Press4 = Press1 * ( ( Temp1 / Temp4 ) ** (Grav * MolWtDry / (GasRUniv * DTempDZ3) ) )
call ECCM_DTempDZ( Temp4, Press4, MolFr, DTempDZ4 )
z_DTempDZ(k) = ( + DTempDZ1 + DTempDZ2 * 2.0d0 + DTempDZ3 * 2.0d0 + DTempDZ4 ) / 6.0d0
!==============================================================
!温度を計算
z_Temp(k+1) = z_Temp(k) + z_DTempDz(k) * DelZ
!圧力を静水圧平衡から計算
z_Press(k+1) = z_Press(k) * ( ( z_Temp(k) / z_Temp(k+1)) ** (Grav * MolWtDry / ( z_DTempDZ(k) * GasRUniv ) ) )
!モル比 (DimZMax でも値を持つようにするため)
z_MolFr(k+1,:) = z_MolFr(k,:)
z_MolFr(k+1,0) = 1.0d0 - sum( z_MolFr(k,1:SpcNum) )
if (z_Temp(k+1) <= 0.0d0 ) exit DtDz
end do DtDz
end subroutine ECCM_Temp_Press