disturbenv.f90

Path: env/disturbenv.f90
Last Update: Mon Mar 28 10:48:53 +0900 2011

    Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved.

Subroutine DisturbEnv

  * Developer: SUGIYAMA Ko-ichiro, ODAKA Masatsugu
  * Version: $Id: disturbenv.f90,v 1.21 2011-03-28 01:48:53 sugiyama Exp $
  * Tag Name: $Name: arare4-20120911 $
  * Change History:

Overview

擾乱のデフォルト値を与えるためのルーチン.

Error Handling

Known Bugs

Note

Future Plans

 * 設定可能な擾乱のタイプを増やす
 * エラー処理に gf4f90io を利用

Required files

Methods

Included Modules

dc_message gridset fileset basicset Boundary ECCM

Public Instance methods

Subroutine :
cfgfile :character(*), intent(in)
xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax) :real(8), intent(out)
: 温位の擾乱成分
xz_Exner(DimXMin:DimXMax,DimZMin:DimZMax) :real(8), intent(out)
: エクスナー関数の擾乱成分
pz_VelX(DimXMin:DimXMax,DimZMin:DimZMax) :real(8), intent(out)
: 水平風速の擾乱成分
xr_VelZ(DimXMin:DimXMax,DimZMin:DimZMax) :real(8), intent(out)
: 鉛直風速の擾乱成分
xza_MixRt(DimXMin:DimXMax,DimZMin:DimZMax, SpcNum) :real(8), intent(out)
: 凝縮成分の混合比(擾乱成分)
xz_Km(DimXMin:DimXMax,DimZMin:DimZMax) :real(8), intent(out)
: 運動量に対する拡散係数
xz_Kh(DimXMin:DimXMax,DimZMin:DimZMax) :real(8), intent(out)
: 熱に対する拡散係数

擾乱のデフォルト値を与えるためのルーチン.

This procedure input/output NAMELIST#disturbset .

[Source]

subroutine DisturbEnv( cfgfile, xz_PotTemp, xz_Exner, pz_VelX, xr_VelZ, xza_MixRt, xz_Km, xz_Kh )
  !
  !擾乱のデフォルト値を与えるためのルーチン. 
  !
  
  !モジュール読み込み
  use dc_message, only: MessageNotify

  use gridset,    only:DimXMin, DimXMax, RegXMin, RegXMax, DimZMin, DimZMax, RegZMin, RegZMax, MarginZ, DelZ, SpcNum, XMin, XMax, ZMin, ZMax, s_X, s_Z               ! Z 座標軸(スカラー格子点)
  use fileset,    only:RandomFile        ! 乱数ファイル
  use basicset,   only: SpcWetMolFr, MolWtWet, MolWtDry, xz_TempBasicZ, xz_PressBasicZ, xz_ExnerBasicZ, xza_MixRtBasicZ   ! 基本場の混合比
  use Boundary, only:  BoundaryXCyc_xza , BoundaryZSym_xza 
  use ECCM,     only:  ECCM_MolFr


  !暗黙の型宣言禁止
  implicit none
  
  !変数定義
  character(*), intent(in) :: cfgfile
  real(8), intent(out)  :: pz_VelX(DimXMin:DimXMax,DimZMin:DimZMax)  
                                    !水平風速の擾乱成分
  real(8), intent(out)  :: xr_VelZ(DimXMin:DimXMax,DimZMin:DimZMax) 
                                    !鉛直風速の擾乱成分 
  real(8), intent(out)  :: xz_Exner(DimXMin:DimXMax,DimZMin:DimZMax)  
                                    !エクスナー関数の擾乱成分 
  real(8), intent(out)  :: xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax)  
                                    !温位の擾乱成分 
  real(8), intent(out)  :: xza_MixRt(DimXMin:DimXMax,DimZMin:DimZMax, SpcNum)  
                                    !凝縮成分の混合比(擾乱成分)
  real(8), intent(out)  :: xz_Km(DimXMin:DimXMax,DimZMin:DimZMax)
                                    !運動量に対する拡散係数
  real(8), intent(out)  :: xz_Kh(DimXMin:DimXMax,DimZMin:DimZMax)
                                    !熱に対する拡散係数
  real(8), parameter         :: Pi = 3.1415926535897932385d0 
                                    !円周率
  real(8)       :: Humidity         !相対湿度
  real(8)       :: XcRate           !擾乱の中心位置(水平方向)の領域に対する割合
  real(8)       :: XrRate           !擾乱の半径(水平方向)の領域に対する割合
  real(8)       :: ZcRate           !擾乱の中心位置(鉛直方向)の領域に対する割合
  real(8)       :: ZrRate           !擾乱の半径(鉛直方向)の領域に対する割合
  real(8)       :: Xc               !擾乱の中心位置(水平方向)
  real(8)       :: Xr               !擾乱の半径(水平方向)
  real(8)       :: Zc               !擾乱の中心位置(鉛直方向)
  real(8)       :: Zr               !擾乱の半径(鉛直方向)
  real(8)       :: Xpos             !擾乱の X 座標 [m] (Therma-Random 用)
  real(8)       :: Zpos             !擾乱の Z 座標 [m] (Therma-Random 用)
  real(8)       :: DryHeight2       !乾燥を与える高度の下限
  real(8)       :: Us               !シアーの強度
  real(8)       :: Hb               !シアーを与える高度の下限
  real(8)       :: Hu               !シアーを与える高度の上限
  real(8)       :: beta(DimXMin:DimXMax, DimZMin:DimZMax)
                                    !擾乱の最大値に対する割合
  real(8)       :: DelMax           !温位擾乱の最大値
  real(8)       :: Random           !ファイルから取得した乱数
  real(8)       :: RandomNum(DimXMin:DimXMax)
  real(8)       :: RandomNum2(DimXMin:DimXMax)
  character(20) :: Type             !温位擾乱のタイプ
  real(8)       :: xz_MolFr(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                    !モル比

  integer       :: i, k, s, n, B, U       ! DO ループ変数

  !-------------------------------------------------------------
  ! 初期化
  !-------------------------------------------------------------
  !NAMELIST ファイルの読み込み
  NAMELIST /disturbset/ Type, DelMax, XrRate, XcRate, ZrRate, ZcRate, Humidity, Xpos, Zpos, DryHeight2, Us, Hu, Hb

  open (10, FILE=cfgfile)
  read(10, NML=disturbset)
  close(10)
  
  !初期化
  pz_VelX    = 0.0d0
  xr_VelZ    = 0.0d0
  xz_Exner   = 0.0d0
  xz_PotTemp = 0.0d0
  xza_MixRt  = 0.0d0
  xz_Km      = 0.0d0
  xz_Kh      = 0.0d0

  Xr = minval( s_X, 1, s_X > (XMax - XMin) * XrRate )
  Xc = minval( s_X, 1, s_X > (XMax - XMin) * XcRate )
  Zr = minval( s_Z, 1, s_Z > (ZMax - ZMin) * ZrRate )
  Zc = minval( s_Z, 1, s_Z > (ZMax - ZMin) * ZcRate )

  !-------------------------------------------------------------
  ! 温位の擾乱
  !-------------------------------------------------------------
  select case(Type)

  case ("Thermal-KW1978")
    ! 指定された領域内に温位擾乱を与える (Klemp and Wilhelmson, 1978) 
    
    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        beta(i,k) = ( ( ( s_X(i) - Xc ) / Xr ) ** 2.0d0 + ( ( s_Z(k) - Zc ) / Zr ) ** 2.0d0 ) ** 5.0d-1
      end do
    end do
    
    where ( beta < 1.0d0 )
      xz_PotTemp = DelMax * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 ) / xz_ExnerBasicZ
    end where
    

  case ("Thermal-Gauss")
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        xz_PotTemp(i,k) = DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) / xz_ExnerBasicZ(i,k)
      end do
    end do

!    where ( sign(1.0d0, DelMax) * xz_PotTemp < DelMax * 1.0d-4 )
!      xz_PotTemp = 0.0d0 
!    end where
    
  case ("Thermal-MixRt-Gauss")
    ! 指定された座標を中心としたガウス分布の温位擾乱と混合比擾乱を与える. 
    ! ただし, 混合比擾乱の最大値は, 基本場の混合比の 0.01 倍とする.
    
    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        xz_PotTemp(i,k) = DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) / xz_ExnerBasicZ(i,k)

        do s = 1, SpcNum
          xza_MixRt(i,:,s) = maxval( xza_MixRtBasicZ(:,:,s) ) * 1.0d-1 * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) / xz_ExnerBasicZ(i,:)
        end do
      end do
    end do
    
!    where ( sign(1.0d0, DelMax) * xz_PotTemp < DelMax * 1.0d-4 )
!      xz_PotTemp = 0.0d0 
!    end where
    
  case ("Exner-Gauss")
    ! 指定された場所に, ガウシアンな温位の擾乱を与える. 

    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        xz_Exner(i,k) = DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) 
      end do
    end do

    where ( xz_Exner < DelMax * 1.0d-4)
      xz_Exner = 0.0d0 
    end where
    
  case ("Thermal-Random")
    ! 下層にランダムな擾乱を与える

    open(20,file=RandomFile)
    do i = DimXMin, DimXMax
      read(20,*) random
      RandomNum(i) = random
    end do
    close(20)

    do i = DimXMin, DimXMax
      !擾乱が全体としてはゼロとなるように調整
      RandomNum2(i) = RandomNum(i) - sum( RandomNum(RegXMin+1:RegXMax) ) / real(RegXMax - RegXMin, 8) 

!      xz_PotTemp(i, maxloc(s_Z, s_Z <= 90.0d3) - MarginZ) = &
!        & DelMax * RandomNum2(i) / xz_ExnerBasicZ(i, maxloc(s_Z, s_Z <= 100.0d3) - MarginZ)

      xz_PotTemp(i, maxloc(s_Z, s_Z <= Zpos) - MarginZ - 1) = DelMax * RandomNum2(i) / xz_ExnerBasicZ(i, maxloc(s_Z, s_Z <= Zpos) - MarginZ - 1)

    end do
 

  case ("HS2001")
    ! Hueso and Sanchez-Lavega を模した設定

    i = ( DimXMax - DimXMin - 10) / 2 
    k = minloc( s_Z, 1, s_Z > 2.5d4 ) - MarginZ
    n = int( 5.0d3 / DelZ )

    xz_PotTemp(i-n:i,k-n:k) = DelMax


  case ("SK1989")
    ! Skamarock and Klemp (1989) の Cold-bubble 実験
    
    xz_PotTemp = 0.0d0

    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax

        beta(i,k) = sqrt( ( ( s_X(i) - Xc ) / Xr ) ** 2.0d0 + ( ( s_Z(k) - Zc ) / Zr ) ** 2.0d0 )
      end do
    end do

    where ( beta < 1.0d0 )
      xz_PotTemp = 0.5d0*DelMax*(cos(pi*beta) + 1.0d0)
    end where

  end select

  !-------------------------------------------------------------
  ! 蒸気圧の設定. 
  !-------------------------------------------------------------
  if (Humidity /= 0.0d0) then 
    do i = DimXMin, DimXMax      
      call ECCM_MolFr( SpcWetMolFr(1:SpcNum), Humidity, xz_TempBasicZ(i,:), xz_PressBasicZ(i,:), xz_MolFr(i,:,:) )
    end do
    
    !気相のモル比を混合比に変換
    do s = 1, SpcNum
      xza_MixRt(:,:,s) = xz_MolFr(:,:,s) * MolWtWet(s) / MolWtDry - xza_MixRtBasicZ(:,:,s)
    end do
  end if

  
  !値が小さくなりすぎないように最低値を与える
!  where (xza_MixRt <= 1.0d-20 )
!    xza_MixRt = 1.0d-20
!  end where
  
  !境界条件
  call BoundaryXCyc_xza( xza_MixRt )
  call BoundaryZSym_xza( xza_MixRt )
!  xza_MixRt = xza_BoundaryXCyc_xza( xza_MixRt )
!  xza_MixRt = xza_BoundaryZSym_xza( xza_MixRt )
  
  !-------------------------------------------------------------!
  ! シアーの設定 (Takemi,2007)                                  !
  !-------------------------------------------------------------!
    !
    != 概要
    !* case "Takemi2007" での計算時に鉛直シアーのある風を与える時に使用する
    !* 風の与え方には, 以下のようなバリエーションがある
    !  (1) シアーを与える高度を変える
    !  (2) シアーのある風の最大風速 (U_s) を変える
    !
    !  (1) については, (a) 0 - 2.5 km, (b) 2.5 - 5.0 km, (c) 5.0 - 7.5 km の
    !  三パターンがある
    !  (2) については, Takemi (2007) では熱帯場と中緯度場の温度場毎に
    !  異なる値を設定している
    !
    !  その強度(Us)は, 以下の通り
    !  <熱帯場>   (1) 5 m/s, (2) 10 m/s, (3) 15 m/s
    !  <中緯度場> (1) 10 m/s, (2) 15 m/s, (3) 20 m/s
    !
    !* シアーの形の模式図 (Takemi, 2007)   |
    !                                     /| 7.5 km
    !                                    / |
    !                                   /  |
    !                                  / ←|
    !                                 |  /| 5.0 km
    !                                 | / |
    !                                 |/  |
    !                                  / ←|
    !                                 |  /| 2.5 km
    !                                 | / |
    !                                 |/  |
    !                                  / ←|
    !---------------------------------+------------- 0.0 km
    !                                Us (m/s)
    !

    ! シアーを与える高度を決める
    ! 傾き始めの高度の格子点 (B) と 傾きの終わりの格子点 (U) で決定

    do k = RegZMin+1, DimZMax
     if (s_Z(k).ge.Hb.and.s_Z(k).le.(Hb + DelZ)) then
       B = k
     end if

     if (s_Z(k).ge.Hu.and.s_Z(k).le.(Hu + DelZ)) then
       U = k
     end if
    end do

    ! 上記で決めた格子点BからUの間に鉛直シアーのある風を与える

    do k = RegZMin+1, U
     do i = DimXMin, DimXMax
       if (k.lt.B) then
         pz_VelX(i,k) = - Us
       else
         pz_VelX(i,k) = - Us + (Us / (s_Z(U) - s_Z(B))) * (s_Z(k) - s_Z(B))
       end if
     end do
    end do
end subroutine DisturbEnv