!--------------------------------------------------------------------- ! Copyright (C) GFD Dennou Club, 2004. All rights reserved. !--------------------------------------------------------------------- !=begin != Module BasicSet ! ! * Developer: SUGIYAMA Ko-ichiro ! * Version: $Id: basicset.f90,v 1.2.2.8 2006/01/07 12:32:05 kitamo Exp $ ! * Tag Name: $Name: $ ! * Change History: ! !== Overview ! !基本場の値を取得する. ! * basicset_file_init: 基本場の値を netCDF ファイルから取得 ! * basicset_calc_init: 基本場の情報を Namelist から取得して値を計算 ! !== Error Handling ! !== Known Bugs ! !== Note ! !== Future Plans ! !=end module BasicSet !== 暗黙の型宣言禁止 implicit none !== save 属性 save !=begin !== Public Interface real(8), allocatable :: ss_CpBasicZ(:,:) !比熱 real(8), allocatable :: ss_DensBasicZ(:,:) !密度 real(8), allocatable :: ss_ExnerBasicZ(:,:) !無次元圧力 real(8), allocatable :: ss_PotTempBasicZ(:,:) !温位 real(8), allocatable :: ss_VelSoundBasicZ(:,:) !音速 !=end contains !=begin !== Procedure Interface ! !=== Initialize module ! !基本場の値を netCDF ファイルから取得する. !netCDF から値を取得するために gt4f90io を利用する. ! subroutine BasicSet_file_init(basicfile) !=== Dependency use dc_trace, only: BeginSub, EndSub use gt4_history use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax !=== Input character(*), intent(in) :: basicfile !=end character(20) :: name !変数名 call BeginSub("basicset_file_init", & & fmt="%c", & & c1="Get basic state variables from netCDF file.") !--- 配列の割り当て allocate( & & ss_CpBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & ss_ExnerBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & ss_PotTempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & ss_DensBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & ss_VelSoundBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) ) !=begin !=== Get a Value from netCDF File name = "CpBasicZ" call HistoryGet( basicfile, name, ss_CpBasicZ ) name = "ExnerBasicZ" call HistoryGet( basicfile, name, ss_ExnerBasicZ ) name = "PotTempBasicZ" call HistoryGet( basicfile, name, ss_PotTempBasicZ ) name = "DensBasicZ" call HistoryGet( basicfile, name, ss_DensBasicZ ) name = "VelSoundBasicZ" call HistoryGet( basicfile, name, ss_VelSoundBasicZ ) !=end call EndSub("basicset_file_init") end subroutine BasicSet_file_init !=begin !== Procedure Interface ! !=== Initialize module and acquire NAMELIST ! !典型的な基本場を計算するためのプログラム basic の初期値を !読み込み, 基本場の値を計算するためのモジュール ! !典型的な基本場のを計算するためのサブルーチン. !タイプとしては以下が設定可能. ! * 温度一様 or 温位一様 ! * 圧力一様 or 静水圧平衡 ! subroutine basicset_calc_init(cfgfile) !=== Dependency use dc_trace, only: BeginSub, EndSub use dc_message, only: MessageNotify use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax, RegZMin, s_Z, DelZ use physset, only: Grav, TempSfc, PressSfc, MolWtDry, CpDry, GasRUniv, GasR use cloudset, only: SatPressA, SatPressB, SatRatioCr !=== Input character(*), intent(in) :: cfgfile !設定ファイル(NAMELIST) !=end !=== Work character(80) :: TempType character(80) :: PressType real(8) :: ss_TempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) ! 平均温度場 real(8) :: ss_CvBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: ss_Z(DimXMin:DimXMax, DimZMin:DimZMax) ! 2-D 座標 !! real(8) :: GasR ! real(8) :: s_TempAdiabat(DimZMin:DimZMax) ! 等温時の温位分布 real(8) :: s_TempSaturate(DimZMin:DimZMax) ! 等温時の温位分布 real(8) :: MidHeight ! 等温と等温位が切り替わる高度(Odaka1998にて使用) real(8) :: TempMidHeight ! 高度 MidHeight での温度 real(8) :: ss_TempAdia(DimXMin:DimXMax, DimZMin:DimZMax) ! 乾燥断熱線に沿った温度 real(8) :: ss_TempSat(DimXMin:DimXMax, DimZMin:DimZMax) ! 凝結線に沿った温度 real(8) :: ss_TempIso(DimXMin:DimXMax, DimZMin:DimZMax) ! 等温線に沿った温度 integer :: k,i !=begin !=== NAMELIST NAMELIST /basicset/ TempType, PressType !=end call BeginSub("basicset_calc_init", & & fmt="%c", & & c1="calculate the basic state variables.") open (10, FILE=cfgfile) read(10, NML=basicset) close(10) !==== 確認 ! write(*,*) "TempType", TempType ! write(*,*) "PressType", PressType !=== 変数の初期化 allocate( ss_CpBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & ss_DensBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & ss_ExnerBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & ss_PotTempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & ss_VelSoundBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) ) ss_CpBasicZ = 0.0d0 ss_DensBasicZ = 0.0d0 ss_ExnerBasicZ = 0.0d0 ss_PotTempBasicZ = 0.0d0 ss_VelSoundBasicZ = 0.0d0 TempMidHeight = 0.0d0 MidHeight = 0.0d0 !--- 座標を初期化 do k = DimZMin, DimZMax ss_Z(:,k) = s_Z(k) end do !-------------------------------------------- !--- 気体定数を決める !-------------------------------------------- !! GasR = GasRUniv / MolWtDry !-------------------------------------------- !--- 比熱を決める !--- 比熱は乾燥大気を念頭に一様とする. !-------------------------------------------- ss_CpBasicZ = CpDry / MolWtDry ss_CvBasicZ = ss_CpBasicZ - GasR !-------------------------------------------- !--- 温度/温位, 無次元圧力の基本場 !-------------------------------------------- if (TempType == "uniform" .AND. PressType == "uniform") then !--- 温度場は TempSfc で一様 ss_TempBasicZ = TempSfc !--- 地表面圧力で一定. 定義より値は 1 となる. ss_ExnerBasicZ = 1.0d0 !--- 温位 ss_PotTempBasicZ = ss_TempBasicZ / ss_ExnerBasicZ elseif (TempType == "uniform" .AND. PressType == "hydrostatic") then !--- 温度場は TempSfc で一様 ss_TempBasicZ = TempSfc !--- 無次元圧力と温位を連立して解く. ss_ExnerBasicZ = & & dexp( - Grav * ss_Z / ( ss_CpBasicZ * ss_TempBasicZ ) ) !--- 無次元圧力を利用して温位を計算 ss_PotTempBasicZ = ss_TempBasicZ / ss_ExnerBasicZ elseif (TempType == "adiabat" .AND. PressType == "uniform") then !--- 地表面での温位(= TempSfc)で一定 ss_PotTempBasicZ = TempSfc !--- 地表面圧力で一定. 定義より値は 1 となる. ss_ExnerBasicZ = 1.0d0 elseif (TempType == "adiabat" .AND. PressType == "hydrostatic") then !--- 地表面での温位(= TempSfc)で一定 ss_PotTempBasicZ = TempSfc !--- 無次元圧力を計算. 温位一定 ss_ExnerBasicZ = & & - Grav * ss_Z / ( ss_CpBasicZ * ss_PotTempBasicZ ) + 1.0d0 elseif (TempType == "odaka1998" .AND. PressType == "hydrostatic") then !--- 温度分布を与える. 各高度毎に等温線と乾燥断熱線の高い方を選ぶ. ss_TempBasicZ = & & max(TempSfc - 25.0d0, TempSfc - Grav * ss_Z / ss_CpBasicZ) !--- 与えた温度分布から静水圧平衡の式を用いて圧力を求める ss_ExnerBasicZ = ss_HydroStaticExnerZ(ss_TempBasicZ) !--- 温度, 圧力から温位を求める ss_PotTempBasicZ = ss_TempBasicZ / ss_ExnerBasicZ elseif (TempType == "adia-sat-iso" .AND. PressType == "hydrostatic") then ss_TempAdia = TempSfc - Grav * ss_Z / ss_CpBasicZ ss_TempIso = 135.0d0 !--- 初期化. とりあえず乾燥断熱分布を与えておく ss_TempBasicZ = ss_TempAdia ss_ExnerBasicZ = ss_HydroStaticExnerZ(ss_TempBasicZ) !--- 下から圧力, 温度場を求め直す do k = RegZMin+1, DimZMax do i = DimXMin, DimXMax ss_TempSat(i,k) = SatPressB & & / ( & & SatPressA & & - log(PressSfc & & * ss_ExnerBasicZ(i,k)**(ss_CpBasicZ(i,k)/GasR) & & / SatRatioCr & & ) & & ) if (ss_TempBasicZ(i,k) < ss_TempSat(i,k) ) then do ss_TempSat(i,k) = SatPressB & & / ( & & SatPressA & & - log( PressSfc & & * ss_ExnerBasicZ(i,k)**(ss_CpBasicZ(i,k)/GasR) & & / SatRatioCr) & & ) if ( abs(ss_TempSat(i,k) - ss_TempBasicZ(i,k)) > epsilon(1.0d0) ) then ss_TempbasicZ(i,k) = ss_TempbasicZ(i,k) + (ss_TempSat(i,k) - ss_TempBasicZ(i,k)) ss_ExnerBasicZ(i,k) = ss_ExnerBasicZ(i,k-1) * & & exp(- Grav / ss_CpBasicZ(i,k-1) * 0.5d0 * DelZ & & * (1.0d0 / ss_TempBasicZ(i,k-1) + 1.0d0 / ss_TempBasicZ(i,k))) else exit end if end do end if ss_TempBasicZ(i,k) = max(ss_TempBasicZ(i,k), ss_TempIso(i,k)) end do end do !--- 上の層の圧力を計算 !! ss_ExnerBasicZ(:,k) = ss_ExnerBasicZ(:,k-1) & !! & * exp(- Grav / ss_CpBasicZ(:,k) * DelZ & ! & * 1.0d0 / ss_TempBasicZ(:,k-1) & ! & ) ! !--- 凝結温度を求める ! ! ノイズによる凝結を防ぐため過飽和度を 1 % 小さくして計算 ! ss_TempSat(:,k) = SatPressB & ! & / ( & ! & SatPressA & ! & - log(PressSfc & ! & * ss_ExnerBasicZ(:,k)**(ss_CpBasicZ(:,k)/GasR) & ! & / SatRatioCr & ! & ) & !! & - log(PressSfc) !! & - log(ss_ExnerBasicZ(:,k)) * ss_CpBasicZ(:,k) /GasR & !! & + log(SatRatio) & ! & ) ! !--- 上の層の温度を決める ! ss_TempBasicZ(:,k) = & ! & max(ss_TempAdia(:,k), ss_TempIso(:,k), ss_TempSat(:,k)) ! end do !--- 温度, 圧力から温位を求める ss_PotTempBasicZ = ss_TempBasicZ / ss_ExnerBasicZ else call MessageNotify("Message", "basicset_calc_init", & & "unknown TempType, or PressType") end if !-------------------------------------------- !---- 密度 !-------------------------------------------- ss_DensBasicZ = & & PressSfc * ( ss_ExnerBasicZ ** (ss_CvBasicZ / GasR) ) & & / ( GasR * ss_PotTempBasicZ ) !-------------------------------------------- !--- 平均場の音速 !-------------------------------------------- ss_VelSoundBasicZ = & & ss_CpBasicZ * GasR * ss_ExnerBasicZ * ss_PotTempBasicZ & & / ss_CvBasicZ ss_VelSoundBasicZ = sqrt( ss_VelSoundBasicZ ) !==== 確認 ! write(*,*) "ss_CpBasicZ ", maxval( ss_CpBasicZ ) ! write(*,*) "ss_DensBasicZ ", maxval( ss_DensBasicZ ) ! write(*,*) "ss_ExnerBasicZ ", maxval( ss_ExnerBasicZ ) ! write(*,*) "ss_PotTempBasicZ ", maxval( ss_PotTempBasicZ ) ! write(*,*) "ss_VelSoundBasicZ ", maxval( ss_VelSoundBasicZ ) call EndSub("basicset_calc_init") end subroutine basicset_calc_init function ss_HydroStaticExnerZ(ss_TempZ) !=begin !== Overviews ! 静水圧平衡を仮定して, 温度分布から無次元圧力分布を求める. ! 地表面から上端まで静水圧平衡の式を鉛直積分することで計算する. ! !== Dependancy use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax, RegZMin, DelZ use physset, only: PressSfc, TempSfc, Grav !==Input real(8), intent(in) :: ss_TempZ(DimXMin:DimXMax, DimZMin:DimZMax) ! 温度 !==Output real(8) :: ss_HydroStaticExnerZ(DimXMin:DimXMax, DimZMin:DimZMax) ! 無次元圧力 !=end !==Work integer :: k ! 初期化. 地表面より下の圧力は計算しないのでこのまま 1 が入る. ss_HydroStaticExnerZ = 1.0d0 ! 最下層の圧力を計算 ss_HydroStaticExnerZ(:,RegZMin+1) = & & exp(- Grav / ss_CpBasicZ(:,RegZMin+1) * 0.25d0 * DelZ & & * (1.0d0 / TempSfc + 1.0d0 / ss_TempZ(:,RegZMin+1) ) & & ) ! ひとつ下の格子点の圧力から次の格子点の圧力を求める do k = RegZMin+2, DimZMax ss_HydroStaticExnerZ(:,k) = ss_HydroStaticExnerZ(:,k-1) & & * exp(- Grav / ss_CpBasicZ(:,k) * 0.5d0 * DelZ & & * (1.0d0 / ss_TempZ(:,k-1) + 1.0d0 / ss_TempZ(:,k)) & & ) end do ! 下側マージン部分の圧力を計算 do k = RegZMin, DimZMin, -1 ss_HydroStaticExnerZ(:,k) = ss_HydroStaticExnerZ(:,k+1) & & * exp(Grav / ss_CpBasicZ(:,k) * 0.5d0 * DelZ & & * (1.0d0 / ss_TempZ(:,k+1) + 1.0d0 / ss_TempZ(:,k)) & & ) end do end function ss_HydroStaticExnerZ end module BasicSet