!--------------------------------------------------------------------- ! 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) ! ! == TODO ! * 水と乾燥空気の比熱が違っていても良いかどうか ! 確認しないといけない. ! * NAMELIST によるパラメータ読み込み 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) ! != 湿潤対流調節スキーム(水少ない近似を緩和したバージョン) 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