!= Module Radiation ! ! Authors:: YAMASHITA Tatsuya, SUGIYAMA Ko-ichiro ! Version:: $Id: radiation_heatbalance.f90,v 1.6 2011-10-04 05:16:27 sugiyama Exp $ ! Tag Name:: $Name: arare5-20111010 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! !モデルの放射過程を計算するためのパッケージ型モジュール !具体的には以下の項を計算するための関数を格納する. ! * 一様冷却 ! !== Error Handling ! !== Bugs ! !== Note ! ! !== Future Plans ! ! module Radiation_HeatBalance ! !モデルの放射過程を計算するためのパッケージ型モジュール !具体的には以下の項を計算するための関数を格納する. ! * 一様冷却 ! !モジュール読み込み use dc_types, only: DP, STRING use dc_iounit, only: FileOpen use dc_message, only: MessageNotify use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut use mpi_wrapper,only: myrank use namelist_util, only: namelist_filename use timeset, only: TimeN use gridSet, only: imin, & !x 方向の配列の下限 & imax, & !x 方向の配列の上限 & jmin, & !y 方向の配列の下限 & jmax, & !y 方向の配列の上限 & kmin, & !z 方向の配列の下限 & kmax, & !z 方向の配列の上限 & nx, ny, nz use axesset, only: z_Z, z_dz !Z 座標軸(スカラー格子点) use constants,only: DayTime, & ! 1 日の長さ [s] & PressSfc, & !地表面圧力 & CpDry, & !定圧比熱 & CvDry, & !定積比熱 & GasRDry !気体定数 use basicset, only: xyz_ExnerBZ, & !エクスナー関数の基本場 & xyz_VelSoundBZ, &! 音速 & xyz_VPTempBZ, &! 仮温位 & xyz_PTempBZ !温位の基本場 use setmargin,only: SetMargin_xyz !暗黙の型宣言禁止 implicit none !private 属性 private !関数を public にする. public Radiation_heatbalance_init public Radiation_heatbalance_forcing !変数定義 real(DP),save :: RadCoolRate = 0.0d0 !一様放射加熱率 [K/day] integer, save :: IdxHeatUp = 0.0d0 !加熱域上限の鉛直座標に対応する整数値 integer, save :: IdxHeatDown = 0.0d0 !冷却域上限の鉛直座標に対応する整数値 integer, save :: IdxCoolUp = 0.0d0 !加熱域上限の鉛直座標に対応する整数値 integer, save :: IdxCoolDown = 0.0d0 !冷却域上限の鉛直座標に対応する整数値 real(DP),save, allocatable :: xyz_RadHeightHeat(:,:,:) !放射加熱が存在する領域 real(DP),save, allocatable :: xyz_RadHeightCool(:,:,:) !放射加熱が存在する領域 contains !!!------------------------------------------------------------------------!!! subroutine Radiation_heatbalance_init ! !NAMELIST から放射強制の設定を取得 ! !暗黙の型宣言禁止 implicit none !入力変数 real(8) :: HeightHeatUp = 0.0d0 !放射強制を与える鉛直領域の上限 real(8) :: HeightHeatDown = 0.0d0 !放射強制を与える鉛直領域の上限 real(8) :: HeightCoolUp = 0.0d0 !放射強制を与える鉛直領域の上限 real(8) :: HeightCoolDown = 0.0d0 !放射強制を与える鉛直領域の上限 integer :: k !ループ変数 integer :: unit ! NAMELIST から情報を取得 NAMELIST /radiation_heatbalance_nml/ & & RadCoolRate, HeightHeatUp, HeightHeatDown, HeightCoolUp, HeightCoolDown call FileOpen(unit, file=namelist_filename, mode='r') read(unit, NML=radiation_heatbalance_nml) close(unit) do k = kmin, kmax-1 if ( z_Z(k) <= HeightHeatUp .AND. HeightHeatUp < z_Z(k+1) ) then IdxHeatUp = k elseif ( z_Z(k) <= HeightHeatDown .AND. HeightHeatDown < z_Z(k+1) ) then IdxHeatDown = k elseif ( z_Z(k) <= HeightCoolUp .AND. HeightCoolUp < z_Z(k+1) ) then IdxCoolUp = k elseif ( z_Z(k) <= HeightCoolDown .AND. HeightCoolDown < z_Z(k+1) ) then IdxCoolDown = k end if end do allocate( xyz_RadHeightHeat(imin:imax, jmin:jmax, kmin:kmax) ) allocate( xyz_RadHeightCool(imin:imax, jmin:jmax, kmin:kmax) ) xyz_RadHeightHeat = 1.0d0 xyz_RadHeightHeat(:,:,IdxHeatDown:IdxHeatUp) = 1.0d0 xyz_RadHeightCool = 0.0d0 xyz_RadHeightCool(:,:,IdxCoolDown:IdxCoolUp) = 1.0d0 ! Output ! if (myrank == 0) then call MessageNotify( "M", & & "Radiation", "RadCoolRate = %f", d=(/RadCoolRate/)) call MessageNotify( "M", & & "Radiation", "HeightHeatUp = %f", d=(/HeightHeatUp/)) call MessageNotify( "M", & & "Radiation", "HeightHeatDown= %f", d=(/HeightHeatDown/)) call MessageNotify( "M", & & "Radiation", "HeightCoolUp = %f", d=(/HeightCoolUp/)) call MessageNotify( "M", & & "Radiation", "HeightCoolDown= %f", d=(/HeightCoolDown/)) end if call HistoryAutoAddVariable( & & varname='PTempRad', & & dims=(/'x','y','z','t'/), & & longname='Radiation term of potential temperature', & & units='K.s-1"', & & xtype='float') call HistoryAutoAddVariable( & & varname='ExnerRad', & & dims=(/'x','y','z','t'/), & & longname='Radiation term of exner function', & & units='s-1"', & & xtype='float') end subroutine Radiation_heatbalance_init !!!------------------------------------------------------------------------!!! subroutine radiation_heatbalance_forcing(xyz_Exner, xyz_PTemp, xyz_DPTempDt, xyz_DExnerDt) ! ! 温位の放射強制項. ! 地表面から RadHeight で指定された高度までの間で空間的に一様な加熱, ! RadHeight から RadHeight2 までの間で空間的な一様な冷却を与える. ! 放射強制全体として加熱と冷却がバランスするように時々刻々一様加熱率 ! および一様冷却率を変化させる. ! 冷却の振幅を与え, 加熱の振幅を調節する. !暗黙の型宣言を禁止 implicit none !変数定義 real(DP), intent(in) :: xyz_Exner(imin:imax, jmin:jmax, kmin:kmax) real(DP), intent(in) :: xyz_PTemp(imin:imax, jmin:jmax, kmin:kmax) real(DP), intent(inout) :: xyz_DPTempDt(imin:imax, jmin:jmax, kmin:kmax) real(DP), intent(inout) :: xyz_DExnerDt(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_DPTempDt0(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_DExnerDt0(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_Rad(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_RadPI(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_DensSum(imin:imax, jmin:jmax, kmin:kmax) !密度(基本場成分と擾乱成分の和) real(DP) :: HeatSum !ある時刻での単位時間当たりの加熱量 real(DP) :: CoolSum !ある時刻での単位時間当たりの冷却量 real(DP) :: RadHeatRate integer :: i, j, k ! 初期化 ! HeatSum = 0.0d0 CoolSum = 0.0d0 xyz_DPTempDt0 = xyz_DPTempDt xyz_DExnerDt0 = xyz_DExnerDt ! 密度の算出 xyz_DensSum = PressSfc & & * (xyz_ExnerBZ + xyz_Exner)**( CvDry /GasRDry ) & & / (GasRDry * (xyz_PTempBZ + xyz_PTemp ) ) ! ある時刻での (加熱量/加熱率) を計算 do k = IdxHeatDown, IdxHeatUp do j = 1, ny do i = 1, nx HeatSum = HeatSum + z_dz(k) * CpDry * xyz_DensSum(i,j,k) end do end do end do ! ある時刻での (冷却量/冷却率) を計算 do k = IdxCoolDown, IdxCoolUp do j = 1, ny do i = 1, nx CoolSum = CoolSum + z_dz(k) * CpDry * xyz_DensSum(i,j,k) end do end do end do ! 加熱率の算出 ! RadCoolRate が負値であることに注意. RadHeatRate = - RadCoolRate * CoolSum / HeatSum xyz_Rad = & & xyz_RadHeightHeat * RadHeatRate / ( xyz_ExnerBZ * DayTime ) & & + xyz_RadHeightCool * RadCoolRate / ( xyz_ExnerBZ * DayTime ) xyz_DPTempDt = xyz_DPTempDt0 + xyz_Rad ! 圧力変化 ! xyz_RadPI = (xyz_VelSoundBZ ** 2.0d0) * xyz_Rad & & / (CpDry * xyz_VPTempBZ ** 2.0d0) xyz_DExnerDt = xyz_DExnerDt0 + xyz_RadPI ! output ! call HistoryAutoPut(TimeN, 'PTempRad', xyz_Rad(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'ExnerRad', xyz_RadPI(1:nx,1:ny,1:nz)) ! Set Margin ! call SetMargin_xyz(xyz_DPTempDt) call SetMargin_xyz(xyz_DExnerDt) end subroutine radiation_heatbalance_forcing end module Radiation_HeatBalance