basicenv.f90

Path: env/basicenv.f90
Last Update: Mon Feb 28 21:00:24 +0900 2011

Subroutine BasicEnv

Authors:SUGIYAMA Koichiro, ODAKA Masatsugu
Version:$Id: basicenv.f90,v 1.28 2011-02-28 12:00:24 sugiyama Exp $
Tag Name:$Name: arare4-20120911 $
Copyright:Copyright (C) GFD Dennou Club, 2006. All rights reserved.
License:See COPYRIGHT

Overview

デフォルトの基本場を設定するための変数参照型モジュール

  * BasicEnvFile_init: 基本場の値を netCDF ファイルから取得
  * BasicEnvCalc_Init: 基本場の情報を Namelist から取得して値を計算

Error Handling

Known Bugs

Note

Future Plans

Required files

Methods

BasicEnv  

Included Modules

dc_message gridset basicset moistset Boundary ECCM chemcalc HeatFlux_N1994

Public Instance methods

Subroutine :

デフォルトの基本場を設定するためのサブルーチン. 基本場を計算し, BasicSet モジュールの値を初期化する.

コンパイルの順序の問題から, 基本場の値(hogeBasicZ な変数)を 計算する部分をBasicSet モジュールから切り離している. ECCM 始め, BasicSet 自体に依存するが hogeBasicZ は use しない 外部サブルーチンを利用するためである.

[Source]

subroutine BasicEnv()
  !
  !デフォルトの基本場を設定するためのサブルーチン. 
  !基本場を計算し, BasicSet モジュールの値を初期化する. 
  !
  !コンパイルの順序の問題から, 基本場の値(hogeBasicZ な変数)を
  !計算する部分をBasicSet モジュールから切り離している. 
  !ECCM 始め, BasicSet 自体に依存するが hogeBasicZ は use しない
  !外部サブルーチンを利用するためである. 
  !

  !モジュール読み込み
  use dc_message, only: MessageNotify

  use gridset,  only: DimXMin, DimXMax, DimZMin, RegXMin, RegZMin, DimZMax, SpcNum, s_Z, s_X, DelZ            !Z 方向の格子点間隔
  use basicset, only: BasicSetArray_Init, PressBasis, GasRDry, CpDry, CvDry, MolWtDry, Grav, SpcWetMolFr, MolWtWet, EnvType, Tropopause, GasRUniv, Humidity, TempStrat, Dhight, TempSfc, SpcWetID        !
  use moistset, only: IdxCC, IdxCG, CondNum         !
  use Boundary, only: BoundaryXCyc_xz, BoundaryZSym_xz, BoundaryXCyc_xza, BoundaryZSym_xza   !  
  use ECCM,     only: ECCM_Dry, ECCM_Wet, ECCM_N1994, ECCM_Takemi2007, HUM, HUM_Takemi2007, HUM_Takemi2007_TDRY, HUM_Takemi2007_MDRY, HUM_Takemi2007ALL       !追加(Takemi2007)
  use chemcalc, only: SvapPress
  use HeatFlux_N1994, only: xza_MixRtFluxBulk
  
  !暗黙の型宣言禁止
  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) :: MolFrIni(SpcNum)
  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
  real(8) :: Weight
!  real(8) :: TakemiRH(DimZMin:DimZMax)
!  real(8) :: TakemiZtr
  real(8) :: z_llhum(DimZMin:DimZMax)
!  real(8) :: TakemiRHall(DimZMin:DimZMax)
  real(8) :: z_YamasakiRH(DimZMin:DimZMax)
  real(8) :: ObAltitude(1:88)
  real(8) :: ObRH(1:88)
!  real(8) :: TakemiRH_TDRY(DimZMin:DimZMax)
!  real(8) :: TakemiRH_MDRY(DimZMin:DimZMax)

!  integer :: i, k, s, tr, ll, bb, uu
  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_PressBasicZ     = 0.0d0
  za_MolFr          = 0.0d0


  !---------------------------------------------------------------
  ! EnvType を元に, 温度, 圧力, 組成を決める
  !---------------------------------------------------------------
  MolFrIni = SpcWetMolFr(1:SpcNum) 

  !場合分け. Dry なら乾燥断熱的に, Wet なら湿潤断熱的な初期場を決める
  select case(EnvType)
  case("Dry")
    call ECCM_Dry( MolFrIni, Humidity, z_TempBasicZ, z_PressBasicZ, za_MolFr )

  ! 対流圏界面より上の扱い
  !   圏界面より上は等温大気とする. 等温位大気から等温大気への遷移は 
  !   tanh を用いてなめらかにつなぐ
  do k = RegZMin+1, DimZMax-1     
    z_DTempDZ(k) = (z_TempBasicZ(k) - z_TempBasicZ(k-1)) / DelZ
  end do
  do k = RegZMin+2, DimZMax
    
    !重みつけの関数を用意. tanh を用いる
    Weight = ( tanh( (s_Z(k) - tropopause ) / Dhight ) + 1.0d0 ) * 5.0d-1
    
    !仮値として温度を計算する. 圏界面より上では TempStrat の等温大気に近づける
    z_TempBasicZ(k) = z_TempBasicZ(k) * ( 1.0d0 - Weight ) + TempStrat * Weight
    
    !温度減率が断熱温度減率より小さくならないように
    DTempDZ = max( z_DTempDZ(k), (z_TempBasicZ(k) - z_TempBasicZ(k-1)) / DelZ )

    !基本場の温度を決める
    z_TempBasicZ(k) = z_TempBasicZ(k-1) + DTempDZ * DelZ 
    
    !圧力を静水圧平衡から計算
    z_PressBasicZ(k) = z_PressBasicZ(k-1) * ( ( z_TempBasicZ(k-1) / z_TempBasicZ(k) ) ** (Grav * MolWtDry / ( DTempDZ * GasRUniv ) ) )
  end do
  
  case("Wet")
    call ECCM_Wet( MolFrIni, Humidity, z_TempBasicZ, z_PressBasicZ,  za_MolFr )

  ! 対流圏界面より上の扱い
  !   圏界面より上は等温大気とする. 等温位大気から等温大気への遷移は 
  !   tanh を用いてなめらかにつなぐ
  do k = RegZMin+1, DimZMax-1     
      z_DTempDZ(k) = (z_TempBasicZ(k) - z_TempBasicZ(k-1)) / DelZ
  end do
  do k = RegZMin+2, DimZMax
    
    !重みつけの関数を用意. tanh を用いる
    Weight = ( tanh( (s_Z(k) - tropopause ) / Dhight ) + 1.0d0 ) * 5.0d-1
    
    !仮値として温度を計算する. 圏界面より上では TempStrat の等温大気に近づける
    z_TempBasicZ(k) = z_TempBasicZ(k) * ( 1.0d0 - Weight ) + TempStrat * Weight
    
    !温度減率が断熱温度減率より小さくならないように
    DTempDZ = max( z_DTempDZ(k), (z_TempBasicZ(k) - z_TempBasicZ(k-1)) / DelZ )

    !基本場の温度を決める
    z_TempBasicZ(k) = z_TempBasicZ(k-1) + DTempDZ * DelZ 
    
    !圧力を静水圧平衡から計算
    z_PressBasicZ(k) = z_PressBasicZ(k-1) * ( ( z_TempBasicZ(k-1) / z_TempBasicZ(k) ) ** (Grav * MolWtDry / ( DTempDZ * GasRUniv ) ) )
  end do
!                                                          !
!-- Yamasaki(1983)の温度と相対湿度の観測値を使用する場合 --!
!                                                          !
 case("N1994")
    call ECCM_N1994( z_TempBasicZ, z_PressBasicZ, za_MolFr )
!                                                          !
!----------- Takemi(2007)の基本場を使用する場合 -----------!
!                                                          !
  case("Takemi2007")
    call ECCM_Takemi2007( z_TempBasicZ, z_PressBasicZ, za_MolFr, z_llhum )
  end select


  !確認のため出力
!  call MessageNotify( "M", "BasicEnv", "Basic State Atmospheric Profiles." )
!  do k = RegZMin+1, DimZMax-1     
!    write(*,*) "temp", k, s_Z(k), z_TempBasicZ(k), z_PressBasicZ(k)
!  end do
  
!  write(*,*) "#######################################"
!  write(*,*) "############--x--#############"
!  write(*,*) "#######################################"
!  do i = DimXMin, DimXMax
!    write(*,*) i, s_X(i)
!  end do

!  write(*,*) "#######################################"
!  write(*,*) "############--飽和蒸気圧--#############"
!  write(*,*) "#######################################"
!  do k = RegZMin+1,DimZMax
!    write(*,*) s_Z(k)/1000,SvapPress(6,z_TempBasicZ(k))
!  end do

!  call HUM(z_YamasakiRH, ObAltitude, ObRH)
!  write(*,*) "#######################################"
!  write(*,*) "#############--水蒸気圧--##############"
!  write(*,*) "#######################################"
!  do k = RegZMin+1,DimZMax
!    write(*,*) s_Z(k)/1000,SvapPress( 6,z_TempBasicZ(k)) * (z_YamasakiRH(k)/100)
!  end do

!   call HUM_Takemi2007_TDRY(TakemiRH_TDRY,TakemiZtr,tr,ll)
!   call HUM_Takemi2007_MDRY(TakemiRH_MDRY,TakemiZtr,tr,ll)
!   call HUM_Takemi2007(TakemiRH,TakemiZtr,tr,ll)
!   call HUM_Takemi2007ALL(TakemiRHall)
!
!     write(*,*)"---------Relative Humidity Check-------------"
!
!    do k = RegZMin+1, ll
!
!     write(*,*) k, s_Z(k), SvapPress(6,z_TempBasicZ(k)), z_PressBasicZ(k), za_MolFr(k,1), TakemiRHall(k)
!
!    end do
!
!    do k = RegZMin+1, DimZMax
!     write(*,*) TakemiZtr,tr,ll,s_Z(ll), TakemiRH(k),s_Z(k)
!    end do
!
!       write(*,*)"------------Takemi RH--------------"
!    do k = RegZMin+1, DimZMax
!      write(*,*) s_Z(k), TakemiRH(k)
!    end do
!
!       write(*,*)"------------Takemi RH TDRY--------------"
!    do k = RegZMin+1, DimZMax
!!      write(*,*) s_Z(k), TakemiRH_TDRY(k)
!
!    end do
!       write(*,*)"------------Takemi RH MDRY--------------"
!    do k = RegZMin+1, DimZMax
!      write(*,*) s_Z(k), TakemiRH_MDRY(k)
!    end do
!
!       write(*,*)"------------Takemi RHall--------------"
!    do k = RegZMin+1, DimZMax
!      write(*,*) s_Z(k), TakemiRHall(k)
!    end do
!
!    call HUM(z_YamasakiRH, ObAltitude, ObRH)
!
!       write(*,*)"------------Yamasaki RH--------------"
!     do k = RegZMin+1,DimZMax
!       write(*,*)s_Z(k), z_YamasakiRH(k)/100
!     end do
!
!     call Height(bb, uu)
!       write(*,*) bb, s_Z(bb), uu, s_Z(uu)
!
  ! 2 次元配列に格納
  do i = DimXMin, DimXMax
    xz_TempBasicZ(i,:)  = z_TempBasicZ  
    xz_PressBasicZ(i,:) = z_PressBasicZ  
  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

!  !値が小さくなりすぎないように最低値を与える
!  where (xza_MixRtBasicZ <= 1.0d-20 )
!    xza_MixRtBasicZ = 1.0d-20
!  end where
  
  !境界条件
  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