!= Subroutine BasicEnv ! ! Authors:: SUGIYAMA Koichiro, ODAKA Masatsugu ! Version:: $Id: basicenv_mmc.f90,v 1.1 2011-03-30 13:53:06 sugiyama Exp $ ! Tag Name:: $Name: arare4-20120911 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! !デフォルトの基本場を設定するための変数参照型モジュール ! * BasicEnvFile_init: 基本場の値を netCDF ファイルから取得 ! * BasicEnvCalc_Init: 基本場の情報を Namelist から取得して値を計算 ! !== Error Handling ! !== Known Bugs ! !== Note ! !== Future Plans ! subroutine BasicEnv_mmc() ! !デフォルトの基本場を設定するためのサブルーチン. !基本場を計算し, BasicSet モジュールの値を初期化する. ! !コンパイルの順序の問題から, 基本場の値(hogeBasicZ な変数)を !計算する部分をBasicSet モジュールから切り離している. !ECCM 始め, BasicSet 自体に依存するが hogeBasicZ は use しない !外部サブルーチンを利用するためである. ! !モジュール読み込み use dc_message, only: MessageNotify use gridset, only: DimXMin, &!配列の X 方向の下限 & RegXMin, &!配列の X 方向の下限 & DimXMax, &!配列の X 方向の上限 & RegXMax, &!配列の X 方向の上限 & DimZMin, &!配列の Z 方向の下限 & RegZMin, &!配列の Z 方向の下限 & DimZMax, &!配列の Z 方向の上限 & RegZMax, &!配列の X 方向の上限 & SpcNum, &!凝縮成分の数 & s_Z, &!スカラー格子点での高度 & DelZ !Z 方向の格子点間隔 use basicset, only: BasicSetArray_Init, &! & PressBasis, &!温位の基準圧力 & GasRDry, &!乾燥成分の定圧比熱 & CpDry, &!乾燥成分の定圧比熱 & CvDry, &!乾燥成分の定積比熱 & MolWtDry, &!乾燥成分の分子量 & Grav, &!重力加速度 & SpcWetMolFr, &!凝縮成分の初期モル比 & MolWtWet, &!凝縮成分の分子量 & EnvType, &!基本場の設定 & Tropopause, &!対流圏界面高度 & GasRUniv, &! & Humidity, &!基本場の相対湿度 & TempStrat, &!成層圏の温度 [k] & Dhight, & !重み関数のパラメータ [m] & TempSfc, & & PressSfc use Boundary, only: BoundaryXCyc_xz, &! & BoundaryZSym_xz, &! & BoundaryXCyc_xza, &! & BoundaryZSym_xza ! use ECCM, only: ECCM_Dry, &! & ECCM_Wet use cloudset, only: SatRtWetAdia use ChemData, only: ChemData_SVapPress_AntoineA, & & ChemData_SVapPress_AntoineB !暗黙の型宣言禁止 implicit none !変数の定義 real(8) :: xz_DensBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_PressBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_ExnerBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_TempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_PotTempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_VelSoundBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xza_MixRtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) real(8) :: xz_EffMolWtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: z_TempBasicZ(DimZMin:DimZMax) real(8) :: z_PressBasicZ(DimZMin:DimZMax) real(8) :: xza_MolFr(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) real(8) :: za_MolFr(DimZMin:DimZMax, SpcNum) real(8) :: xza_MixRtDivMolWt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) real(8) :: z_DTempDZ(DimZMin:DimZMax) real(8) :: z_MolWtMean(DimZMin:DimZMax) real(8) :: DTempDZ(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: Weight1(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: Weight2(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: LapseRateRatio ! 乾燥断熱減率に対する下層での温度減率の比 ! 乾燥断熱よりも小さく設定する. real(8) :: xz_TempAdia(DimXMin:DimXMax, DimZMin:DimZMax) ! 乾燥断熱線に沿った温度 real(8) :: xz_TempSat(DimXMin:DimXMax, DimZMin:DimZMax) ! 凝結線に沿った温度 real(8) :: xz_TempIso(DimXMin:DimXMax, DimZMin:DimZMax) ! 等温線に沿った温度 real(8) :: xz_PressAdia(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_PressSat(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_PressIso(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_TempWork(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_PressWork(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_Z(DimXMin:DimXMax, DimZMin:DimZMax) ! 2D 座標 real(8) :: Temp_0, Temp_1, Press_0, Press_1 ! 乾燥断熱線と湿潤断熱線とが交わる高度を反復法で求める ! 際に用いる作業変数 real(8) :: Work ! 作業変数 real(8) :: LCL ! 下層の温度分布と湿潤断熱線が交わる高度 real(8) :: LTP ! 湿潤断熱線と等温線が交わる高度 integer :: i, k, s !--------------------------------------------------------------- ! 配列の初期化 !--------------------------------------------------------------- xz_PressBasicZ = 0.0d0 xz_ExnerBasicZ = 0.0d0 xz_TempBasicZ = 0.0d0 xz_PotTempBasicZ = 0.0d0 xz_VelSoundBasicZ = 0.0d0 xza_MixRtBasicZ = 0.0d0 xz_EffMolWtBasicZ = 0.0d0 ! z_TempBasicZ = 0.0d0 z_TempBasicZ = TempSfc z_PressBasicZ = 0.0d0 za_MolFr = 0.0d0 ! 座標の初期化 do k = DimZMin, DimZMax xz_Z(:,k) = s_Z(k) end do ! 北守修論計算の基本場設定(極冠周縁での典型的温度プロファイル)を改良したもの. ! * 温度勾配の不連続を解消する為, 重み付き関数でなめらかに接続. ! * 基本場の温度分布が超断熱, 過飽和とならないよう調整. ! 乾燥断熱減率に対する下層の温度減率の比 LapseRateRatio = 0.8d0 xz_TempAdia = TempSfc - LapseRateRatio * Grav * xz_Z / CpDry xz_TempIso = 135.0d0 !--- 乾燥断熱線, 湿潤断熱線, 等温線が交わる高度を計算し, !--- 各領域で成り立つ式を用いて温度, 圧力を計算 !--- 乾燥断熱線と湿潤断熱線が交わる高度(LCL)を反復法で計算 Press_0 = PressSfc Temp_0 = TempSfc do Temp_1 = ChemData_SVapPress_AntoineB(12) / & & (ChemData_SVapPress_AntoineA(12) - dlog(Press_0/SatRtWetAdia)) Press_1 = PressSfc * & & (Temp_1/TempSfc)**(CpDry / (LapseRateRatio * GasRDry)) if (abs(Temp_1 - Temp_0) < epsilon(0.0d0)) then LCL = (TempSfc * CpDry) / (Grav * LapseRateRatio) & & * (1.0d0 - (Press_1/PressSfc)**( GasRDry/CpDry )) exit else Temp_0 = Temp_1 Press_0 = Press_1 end if end do !--- 湿潤断熱線と等温線が交わる高度(LTP)を計算 LTP = LCL + GasRDry * ChemData_SVapPress_AntoineB(12) & & / Grav * dlog(Temp_1/xz_TempIso(1,1)) !--- LCL, LTP の値を表示. write(*,*) 'LCL', LCL write(*,*) 'LTP', LTP !--- 温度, (圧力), 無次元圧力分布を計算する xz_TempSat = & & Temp_1 * exp(-Grav * (xz_Z - LCL) & & /(GasRDry*ChemData_SVapPress_AntoineB(12))) LCL = 3.0d3 Weight1 = ( tanh( (xz_Z - LCL ) / Dhight ) + 1.0d0 ) * 5.0d-1 Weight2 = ( tanh( (xz_Z - LTP ) / Dhight ) + 1.0d0 ) * 5.0d-1 ! 温度の不連続をなくす為に, tanh を用いてなめらかにつなぐ. xz_TempWork = & & xz_TempAdia * ( 1.0d0 - Weight1 ) & & + xz_TempSat * Weight1 xz_TempBasicZ = & & xz_TempWork * ( 1.0d0 - Weight2 ) + xz_TempIso * Weight2 & & - 0.92d0 ! 圧力を静水圧平衡から計算 do i = DimXMin, DimXMax DTempDZ(i,RegZMin) = & & (xz_TempBasicZ(i,RegZMin) - TempSfc) & & / ( 0.5d0 * DelZ) xz_PressBasicZ(i,RegZMin) = & & Presssfc * & & ( ( Tempsfc / xz_TempBasicZ(i,RegZMin) ) & & ** (Grav / ( DTempDZ(i,RegZMin) * GasRDry ) ) ) end do do k = RegZMin+1, DimZMax-1 do i = DimXMin, DimXMax ! 局所的な温度減率 DTempDZ(i,k) = (xz_TempBasicZ(i,k) - xz_TempBasicZ(i,k-1)) / DelZ xz_PressBasicZ(i,k) = & & xz_PressBasicZ(i,k-1) * & & ( ( xz_TempBasicZ(i,k-1) / xz_TempBasicZ(i,k) ) & & ** (Grav / ( DTempDZ(i,k) * GasRDry ) ) ) end do end do !確認のため出力 call MessageNotify( "M", "BasicEnv", "Basic State Atmospheric Profiles." ) do k = RegZMin+1, DimZMax-1 write(*,*) "temp", k, s_Z(k), xz_TempBasicZ(1,k), xz_PressBasicZ(1,k) end do !境界条件 call BoundaryXCyc_xz( xz_TempBasicZ ) call BoundaryZSym_xz( xz_TempBasicZ ) call BoundaryXCyc_xz( xz_PressBasicZ ) call BoundaryZSym_xz( xz_PressBasicZ ) !--------------------------------------------------------------- ! 混合比 !--------------------------------------------------------------- !水平方向には一様 do i = DimXMin, DimXMax xza_MolFr(i,:,:) = za_MolFr end do !気相のモル比を混合比に変換 do s = 1, SpcNum xza_MixRtBasicZ(:,:,s) = xza_MolFr(:,:,s) * MolWtWet(s) / MolWtDry end do !境界条件 call BoundaryXCyc_xza( xza_MixRtBasicZ ) call BoundaryZSym_xza( xza_MixRtBasicZ ) !--------------------------------------------------------------- ! 分子量の効果 !--------------------------------------------------------------- do s = 1, SpcNum xza_MixRtDivMolWt(:,:,s) = xza_MixRtBasicZ(:,:,s) / MolWtWet(s) end do xz_EffMolWtBasicZ = & & (1.0d0 + sum(xza_MixRtBasicZ,3) ) & & / ( MolWtDry * ((1.0d0 / MolWtDry) + sum(xza_MixRtDivMolWt,3)) ) !境界条件 call BoundaryXCyc_xz( xz_EffMolWtBasicZ ) call BoundaryZSym_xz( xz_EffMolWtBasicZ ) !--------------------------------------------------------------- ! 温位 !--------------------------------------------------------------- xz_PotTempBasicZ = & & xz_TempBasicZ * (PressBasis / xz_PressBasicZ) ** (GasRDry / CpDry) !境界条件 call BoundaryXCyc_xz( xz_PotTempBasicZ ) call BoundaryZSym_xz( xz_PotTempBasicZ ) !--------------------------------------------------------------- ! エクスナー関数 !--------------------------------------------------------------- xz_ExnerBasicZ = xz_TempBasicZ / xz_PotTempBasicZ !境界条件 call BoundaryXCyc_xz( xz_ExnerBasicZ ) call BoundaryZSym_xz( xz_ExnerBasicZ ) !--------------------------------------------------------------- ! 密度 !--------------------------------------------------------------- xz_DensBasicZ = & & PressBasis * (xz_ExnerBasicZ ** (CvDry / GasRDry)) & & / (GasRDry * xz_PotTempBasicZ / xz_EffMolWtBasicZ) !境界条件 call BoundaryXCyc_xz( xz_DensBasicZ ) call BoundaryZSym_xz( xz_DensBasicZ ) !--------------------------------------------------------------- ! 音速 !--------------------------------------------------------------- xz_VelSoundBasicZ = & & sqrt ( & & CpDry * GasRDry * xz_ExnerBasicZ * xz_PotTempBasicZ & & / (CvDry * xz_EffMolWtBasicZ) & & ) !境界条件 call BoundaryXCyc_xz( xz_VelSoundBasicZ ) call BoundaryZSym_xz( xz_VelSoundBasicZ ) !---------------------------------------------------------- ! BasicSet モジュールに値を設定 !---------------------------------------------------------- call BasicSetArray_Init( & & xz_PressBasicZ, xz_ExnerBasicZ, xz_TempBasicZ, & & xz_PotTempBasicZ, xz_DensBasicZ, xz_VelSoundBasicZ, & & xza_MixRtBasicZ, xz_EffMolWtBasicZ ) end subroutine BasicEnv_mmc