disturbenv_3d.f90

Path: env/disturbenv_3d.f90
Last Update: Wed Jun 25 05:16:40 +0900 2008

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

Subroutine DisturbEnv_3d

  * Developer: SUGIYAMA Ko-ichiro, ODAKA Masatsugu
  * Version: $Id: disturbenv_3d.f90,v 1.6 2008-06-24 20:16:40 odakker2 Exp $
  * Tag Name: $Name: arare4-20111215 $
  * Change History:

Overview

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

Error Handling

Known Bugs

Note

Future Plans

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

Required files

Methods

Included Modules

dc_types dc_iounit dc_message gridset_3d fileset_3d basicset_3d xyz_module

Public Instance methods

Subroutine :
cfgfile :character(*), intent(in)
xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) :real(DP), intent(out)
: 温位の擾乱成分
xyz_Exner(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) :real(DP), intent(out)
: エクスナー関数の擾乱成分
pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) :real(DP), intent(out)
: 水平風速の擾乱成分
xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) :real(DP), intent(out)
: 水平風速の擾乱成分
xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) :real(DP), intent(out)
: 鉛直風速の擾乱成分
xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum) :real(DP), intent(out)
: 凝縮成分の混合比(擾乱成分)
xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) :real(DP), intent(out)
: 運動量に対する拡散係数
xyz_Kh(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) :real(DP), intent(out)
: 熱に対する拡散係数

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

This procedure input/output NAMELIST#disturbset .

[Source]

subroutine DisturbEnv_3d( cfgfile, xyz_PotTemp, xyz_Exner, pyz_VelX, xqz_VelY, xyr_VelZ, xyza_MixRt, xyz_Km, xyz_Kh )
  !
  !擾乱のデフォルト値を与えるためのルーチン. 
  !
  
  !モジュール読み込み
  use dc_types, only : DP
  use dc_iounit,  only: FileOpen
  use dc_message, only: MessageNotify

  use gridset_3d, only:DimXMin, DimXMax, DimYMin, DimYMax, DimZMin, DimZMax, Zmargin, SpcNum, XMin, XMax, YMin, YMax, ZMin, ZMax, x_X, y_Y, z_Z, xyz_X, xyz_Y, xyz_Z, x_dx, y_dy, z_dz
  use fileset_3d, only:RandomFile        ! 乱数ファイル
  use basicset_3d,only: SpcWetMolFr, MolWtWet, MolWtDry, xyz_TempBasicZ, xyz_PressBasicZ, xyz_ExnerBasicZ   ! 無次元圧力
!    &                  xyza_MixRtBasicZ  ! 基本場の混合比

  use xyz_module, only: BoundaryXCyc_xyz, BoundaryYCyc_xyz, BoundaryZSym_xyz 
!  use ECCM_3d,     only:  ECCM_MolFr


  !暗黙の型宣言禁止
  implicit none
  
  !変数定義
  character(*), intent(in) :: cfgfile
  real(DP), intent(out)  :: pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)  
                                    !水平風速の擾乱成分
  real(DP), intent(out)  :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)  
                                    !水平風速の擾乱成分
  real(DP), intent(out)  :: xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) 
                                    !鉛直風速の擾乱成分 
  real(DP), intent(out)  :: xyz_Exner(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)  
                                    !エクスナー関数の擾乱成分 
  real(DP), intent(out)  :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)  
                                    !温位の擾乱成分 
  real(DP), intent(out)  :: xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum)  
                                    !凝縮成分の混合比(擾乱成分)
  real(DP), intent(out)  :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                    !運動量に対する拡散係数
  real(DP), intent(out)  :: xyz_Kh(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                    !熱に対する拡散係数
  real(DP), parameter         :: Pi = 3.1415926535897932385d0 
                                    !円周率
  real(DP)       :: Humidity        !相対湿度
  real(DP)       :: XcRate          !擾乱の中心位置(水平方向)の領域に対する割合
  real(DP)       :: XrRate          !擾乱の半径(水平方向)の領域に対する割合
  real(DP)       :: YcRate          !擾乱の中心位置(水平方向)の領域に対する割合
  real(DP)       :: YrRate          !擾乱の半径(水平方向)の領域に対する割合
  real(DP)       :: ZcRate          !擾乱の中心位置(鉛直方向)の領域に対する割合
  real(DP)       :: ZrRate          !擾乱の半径(鉛直方向)の領域に対する割合
  real(DP)       :: Xc              !擾乱の中心位置(水平方向)
  real(DP)       :: Xr              !擾乱の半径(水平方向)
  real(DP)       :: Yc              !擾乱の中心位置(水平方向)
  real(DP)       :: Yr              !擾乱の半径(水平方向)
  real(DP)       :: Zc              !擾乱の中心位置(鉛直方向)
  real(DP)       :: Zr              !擾乱の半径(鉛直方向)
  real(DP)       :: Xpos            !擾乱の X 座標 [m] (Therma-Random 用)
  real(DP)       :: Ypos            !擾乱の Y 座標 [m] (Therma-Random 用)
  real(DP)       :: Zpos            !擾乱の Z 座標 [m] (Therma-Random 用)
  real(DP)       :: beta(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                    !擾乱の最大値に対する割合
  real(DP)       :: DelMax          !温位擾乱の最大値
  real(DP)       :: Random          !ファイルから取得した乱数
  real(DP)       :: RandomNum(DimXMin:DimXMax,DimYMin:DimYMax)
  real(DP)       :: RandomNum2(DimXMin:DimXMax,DimYMin:DimYMax)
  real(DP)       :: RandomMean       
  character(20)  :: Type            !温位擾乱のタイプ
!  real(DP)       :: xyza_MolFr(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                    !モル比

  integer       :: i, j, k, s, n    ! DO ループ変数
  integer       :: unit             ! 装置番号

  !-------------------------------------------------------------
  ! 初期化
  !-------------------------------------------------------------
  !NAMELIST ファイルの読み込み
  NAMELIST /disturbset/ Type, DelMax, XrRate, XcRate, YrRate, YcRate, ZrRate, ZcRate, Humidity, Xpos, Ypos, Zpos

  call FileOpen(unit, file=cfgfile, mode='r')
  read(unit, NML=disturbset)
  close(unit)

  !初期化
  pyz_VelX    = 0.0d0
  xqz_VelY    = 0.0d0
  xyr_VelZ    = 0.0d0
  xyz_Exner   = 0.0d0
  xyz_PotTemp = 0.0d0
  xyza_MixRt  = 0.0d0
  xyz_Km      = 0.0d0
  xyz_Kh      = 0.0d0

  Xr = minval( x_X, 1, x_X > (XMax - XMin) * XrRate )
  Xc = minval( x_X, 1, x_X > (XMax - XMin) * XcRate )
  Yr = minval( y_Y, 1, y_Y > (YMax - YMin) * YrRate )
  Yc = minval( y_Y, 1, y_Y > (YMax - YMin) * YcRate )
  Zr = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZrRate )
  Zc = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZcRate )

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

  case ("Thermal-KW1978")
    ! 指定された領域内に温位擾乱を与える (Klemp and Wilhelmson, 1978) 
    
    beta = ( ( ( xyz_X - Xc ) / Xr ) ** 2.0d0 + ( ( xyz_Y - Yc ) / Yr ) ** 2.0d0 + ( ( xyz_Z - Zc ) / Zr ) ** 2.0d0 ) ** 5.0d-1
    
    where ( beta < 1.0d0 )
      xyz_PotTemp = DelMax * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 ) / xyz_ExnerBasicZ
    end where
    

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

    xyz_PotTemp = DelMax * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 - ( (xyz_Z - Zc) / Yr )**2.0d0 * 5.0d-1 ) / xyz_ExnerBasicZ

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

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

    xyz_PotTemp = DelMax * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) / xyz_ExnerBasicZ

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where
  
  case ("Thermal-GaussYZ")
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = DelMax * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) / xyz_ExnerBasicZ

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where
 
  case ("Adv-XZ-X")
    ! X 方向移流のテスト
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = DelMax * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

    where ( xyz_PotTemp > DelMax  )
      xyz_PotTemp = DelMax
    end where

    pyz_VelX = 20.0d0

  case ("Adv-XZ-Z")
    ! Z 方向移流のテスト
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = DelMax * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

    where ( xyz_PotTemp > DelMax  )
      xyz_PotTemp = DelMax
    end where

    xyr_VelZ = 20.0d0

  case ("Adv-YZ-Y")
    ! Y 方向移流のテスト
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = DelMax * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

    where ( xyz_PotTemp > DelMax  )
      xyz_PotTemp = DelMax
    end where

    xqz_VelY = 20.0d0

  case ("Adv-YZ-Z")
    ! Z 方向移流のテスト
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = DelMax * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

    where ( xyz_PotTemp > DelMax  )
      xyz_PotTemp = DelMax
    end where

    xyr_VelZ = 20.0d0

  case ("Adv-XZ-XZ")
    ! XZ 方向移流のテスト
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = DelMax * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

    where ( xyz_PotTemp > DelMax  )
      xyz_PotTemp = DelMax
    end where

    pyz_VelX = 20.0d0
    xyr_VelZ = 20.0d0

  case ("Adv-YZ-YZ")
    ! YZ 方向移流のテスト
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = DelMax * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

    where ( xyz_PotTemp > DelMax  )
      xyz_PotTemp = DelMax
    end where

    xqz_VelY = 20.0d0
    xyr_VelZ = 20.0d0

  case ("Exner-Gauss")
    ! 指定された場所に, ガウシアンな圧力の擾乱を与える. 

    xyz_Exner = DelMax * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 ) * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 ) * dexp( - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( xyz_Exner < DelMax * 1.0d-4)
      xyz_Exner = 0.0d0 
    end where
    
  case ("Exner-GaussX")
    ! 指定された場所に, ガウシアンな圧力の擾乱を与える. 
    ! (X 方向一様)

    xyz_Exner = DelMax * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 ) * dexp( - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( xyz_Exner < DelMax * 1.0d-4)
      xyz_Exner = 0.0d0 
    end where

  case ("Exner-GaussY")
    ! 指定された場所に, ガウシアンな圧力の擾乱を与える. 
    ! (Y 方向一様)

    xyz_Exner = DelMax * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 ) * dexp( - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( xyz_Exner < DelMax * 1.0d-4)
      xyz_Exner = 0.0d0 
    end where

  case ("Exner-GaussZ")
    ! 指定された場所に, ガウシアンな圧力の擾乱を与える. 
    ! (Z 方向一様)

    xyz_Exner = DelMax * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 ) * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 ) 

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

    call FileOpen(unit, file=RandomFile, mode='r')
    do j = DimYMin, DimYMax
      do i = DimXMin, DimXMax
        read(unit,*) random
        RandomNum(i,j) = random
    end do
    end do
    close(unit)

    ! 乱数の平均値を計算
    RandomMean = 0.0d0

    do j = DimYMin, DimYMax
      do i = DimXMin, DimXMax
        RandomMean = RandomMean + (RandomNum(i,j)*x_dx(i)*y_dy(j))/((XMax-Xmin)*(YMax-YMin))
      end do
    end do

    do j = DimYMin, DimYMax
      do i = DimXMin, DimXMax
        !擾乱が全体としてはゼロとなるように調整
        RandomNum2(i,j) = RandomNum(i,j)  - RandomMean

        xyz_PotTemp(i,j,maxloc(z_Z, z_Z <= Zpos) - Zmargin ) = DelMax * RandomNum2(i,j) / xyz_ExnerBasicZ(i,j,maxloc(z_Z, z_Z <= Zpos) - Zmargin - 1)
        
      end do
    end do

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

    i = ( DimXMax - DimXMin - 10) / 2 
    j = ( DimYMax - DimYMin - 10) / 2 
    k = minloc( z_Z, 1, z_Z > 2.5d4 ) - Zmargin
    n = int( 5.0d3 / z_dz(1) )

    xyz_PotTemp(i-n:i,j-n:j,k-n:k) = DelMax


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

    beta = sqrt( ( ( xyz_X - Xc ) / Xr ) ** 2.0d0 + ( ( xyz_Y - Yc ) / Yr ) ** 2.0d0 + ( ( xyz_Z - Zc ) / Zr ) ** 2.0d0 )

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

  end select

  !-------------------------------------------------------------
  ! 蒸気圧の設定. 
  !-------------------------------------------------------------
!  if (Humidity /= 0.0d0) then 
!    do i = DimXMin, DimXMax      
!      do j = DimYMin, DimYMax      
!        call ECCM_MolFr( SpcWetMolFr(1:SpcNum), Humidity, &
!          &              xyz_TempBasicZ(i,j,:),           &
!          &              xyz_PressBasicZ(i,j,:), xyza_MolFr(i,j,:,:) )
!      end do
!    end do
    
    !気相のモル比を混合比に変換
!    do s = 1, SpcNum
!      xyza_MixRt(:,:,:,s) =  &
!        & xyza_MolFr(:,:,:,s) &
!        & * MolWtWet(s) / MolWtDry - xyza_MixRtBasicZ(:,:,:,s)
!    end do
!  end if
  
  
  !値が小さくなりすぎないように最低値を与える
!  where (xza_MixRt <= 1.0d-20 )
!    xza_MixRt = 1.0d-20
!  end where
  
  !境界条件
!  do s =1, SpcNum
!    call BoundaryXCyc_xyz( xyza_MixRt(:,:,:,s) )
!    call BoundaryYCyc_xyz( xyza_MixRt(:,:,:,s) )
!    call BoundaryZSym_xyz( xyza_MixRt(:,:,:,s) )
!  end do

end subroutine DisturbEnv_3d