!= Module cloudphys_k1969 ! ! Authors:: 杉山耕一朗(SUGIYAMA Ko-ichiro), 小高正嗣 (ODAKA Masatsugu), 高橋芳幸 (YOSHIYUKI Takahashi) ! Version:: $Id: cloudphys_k1969.f90,v 1.15 2011-10-10 15:43:00 yot Exp $ ! Tag Name:: $Name: arare5-20111010 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! !暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める. ! * 中島健介 (1994) で利用した定式をそのまま利用. ! !== Error Handling ! !== Bugs ! !== Note ! !== Future Plans ! ! module cloudphys_k1969 ! !暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める. ! * 中島健介 (1994) で利用した定式をそのまま利用. ! !モジュール読み込み 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 timeset, only: DelTimeLong, TimeN use gridset, only : imin, &!x 方向の配列の下限 & imax, &!x 方向の配列の上限 & jmin, &!y 方向の配列の上限 & jmax, &!y 方向の配列の上限 & kmin, &!z 方向の配列の下限 & kmax, &!z 方向の配列の上限 & nx, ny, nz, ncmax !物理領域の大きさ use constants,only: PressBasis, &!温位の基準圧力 & CpDry, &!乾燥成分の比熱 & MolWtDry, &! & GasRDry !乾燥成分の気体定数 use basicset, only: xyz_DensBZ, &!基本場の密度 & xyz_PTempBZ, &!基本場の温位 & xyz_ExnerBZ, &!基本場の無次元圧力 & xyzf_QMixBZ !基本場の混合比 use composition, only: & & MolWtWet, &! & SpcWetID, &! & SpcWetSymbol, &! & CondNum, &!凝結過程の数 & IdxCG, &!凝結過程(蒸気)の配列添え字 & IdxCC, &!凝結過程(雲)の配列添え字 & IdxCR, &!凝結過程(雲)の配列添え字 & CloudNum, &!雲の数 & RainNum, &!雨の数 & IdxC, &!雲の配列添え字 & IdxR, &!雨の配列添え字 & IdxNH3, &!NH3(蒸気)の配列添え字 & IdxH2S, &!H2S(蒸気)の配列添え字 & IdxNH4SHr !NH4SH(雨)の配列添え字 use axesset, only : xyr_avr_xyz use xyz_deriv_module,only : xyz_dz_xyr use ChemCalc, only : xyz_SvapPress, xyz_LatentHeat, ReactHeatNH4SH, xyz_DelQMixNH4SH use MoistAdjust, only: MoistAdjustSvapPress, MoistAdjustNH4SH use setmargin, only: SetMargin_xyzf, SetMargin_xyz use namelist_util, only: namelist_filename use setmargin,only: SetMargin_xyz, SetMargin_xyzf !暗黙の型宣言禁止 implicit none !属性の指定 private !関数を public にする public Cloudphys_K1969_Init public Cloudphys_K1969_forcing public Cloudphys_K1969_FallRain real(DP), save :: FactorJ = 1.0d0 !雲物理過程のパラメータ !木星では 3.0d0 !地球では 1.0d0 とする real(DP), save :: AutoConvTime = 1.0d3 !併合成長の時定数 [sec] real(DP), save :: QMixCr = 1.0d-3 !併合成長を生じる臨界混合比 [kg/kg] contains !!!=================================================================================!!! subroutine Cloudphys_K1969_Init !暗黙の型宣言禁止 implicit none !内部変数 integer :: unit !装置番号 integer :: l character(STRING) :: Planet = "" !----------------------------------------------------------- ! NAMELIST から情報を取得 !----------------------------------------------------------- ! NAMELIST から情報を取得 NAMELIST /cloudphys_k1969_nml/ Planet, FactorJ, AutoConvTime, QMixCr call FileOpen(unit, file=namelist_filename, mode='r') read(unit, NML=cloudphys_k1969_nml) close(unit) if (trim(Planet) == "Earth") then FactorJ = 1.0d0 elseif (trim(Planet) == "Jupiter") then FactorJ = 3.0d0 end if if (myrank == 0) then call MessageNotify( "M", & & "Cloudphys_K1969_Init", "Planet = %c", c1=trim(Planet)) call MessageNotify( "M", & & "Cloudphys_K1969_Init", "FactorJ = %f", d=(/FactorJ/) ) call MessageNotify( "M", & & "Cloudphys_K1969_Init", "AutoConvTime = %f", d=(/AutoConvTime/) ) call MessageNotify( "M", & & "Cloudphys_K1969_Init", "QMixCr = %f", d=(/QMixCr/) ) end if call HistoryAutoAddVariable( & & varname='PTempCond',& & dims=(/'x','y','z','t'/), & & longname='Latent heat term of potential temperature', & & units='K.s-1', & & xtype='float') do l = 1, ncmax call HistoryAutoAddVariable( & & varname=trim(SpcWetSymbol(l))//'_Cond', & & dims=(/'x','y','z','t'/), & & longname='Condensation term of ' & & //trim(SpcWetSymbol(l))//' mixing ratio', & & units='kg.kg-1.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname=trim(SpcWetSymbol(l))//'_Fall', & & dims=(/'x','y','z','t'/), & & longname='Fall Rain term of ' & & //trim(SpcWetSymbol(l))//' mixing ratio', & & units='kg.kg-1.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname=trim(SpcWetSymbol(l))//'_FallFluxAtLB', & & dims=(/'x','y','t'/), & & longname='Falling Rain Flux ' & & //trim(SpcWetSymbol(l)), & & units='kg.m-2.s-1', & & xtype='float') end do end subroutine Cloudphys_K1969_Init !!!=================================================================================!!! subroutine Cloudphys_K1969_forcing(xyz_ExnerNl, xyz_PTempAl, xyzf_QMixAl) implicit none real(DP), intent(in) :: xyz_ExnerNl(imin:imax, jmin:jmax, kmin:kmax) real(DP), intent(inout) :: xyz_PTempAl(imin:imax, jmin:jmax, kmin:kmax) real(DP), intent(inout) :: xyzf_QMixAl(imin:imax, jmin:jmax, kmin:kmax, ncmax) real(DP) :: xyz_PTempOrig(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_PTempWork(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_DelPTemp(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_Del(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyzf_QMixOrig(imin:imax, jmin:jmax, kmin:kmax, ncmax) real(DP) :: xyzf_QMixWork(imin:imax, jmin:jmax, kmin:kmax, ncmax) real(DP) :: xyzf_DelQMix(imin:imax, jmin:jmax, kmin:kmax, ncmax) real(DP) :: xyzf_Del(imin:imax, jmin:jmax, kmin:kmax, ncmax) real(DP) :: DelTime integer :: l, s real(DP) :: xyzf_Cloud2Rain(imin:imax,jmin:jmax,kmin:kmax, ncmax) !雲から雨への変換量 real(DP) :: xyz_AutoConv(imin:imax,jmin:jmax,kmin:kmax) !飽和混合比 real(DP) :: xyz_Collect(imin:imax,jmin:jmax,kmin:kmax) !規格化された潜熱 real(DP) :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax) !混合比の擾乱成分 + 平均成分 real(DP) :: xyz_TempAll(imin:imax,jmin:jmax,kmin:kmax) !温度の擾乱成分 + 平均成分 real(DP) :: xyz_PressAll(imin:imax,jmin:jmax,kmin:kmax) !全圧 real(DP) :: xyz_ExnerAll(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_NonSaturate(imin:imax,jmin:jmax,kmin:kmax) !未飽和度(飽和混合比と蒸気の混合比の差) real(DP) :: xyzf_Rain2Gas(imin:imax,jmin:jmax,kmin:kmax, ncmax) real(DP) :: xyzf_Rain2GasNH4SH(imin:imax,jmin:jmax,kmin:kmax, ncmax) real(DP) :: xyzf_DelPTemp(imin:imax,jmin:jmax,kmin:kmax, ncmax) real(DP) :: xyz_DelPTempNH4SH(imin:imax,jmin:jmax,kmin:kmax) !----------------------------------------- ! 時間刻み幅. Leap-frog なので, 2 \del t ! DelTime = 2.0d0 * DelTimeLong !------------------------------------------ ! 初期値を保管 Store Initial Value ! xyz_PTempOrig = xyz_PTempAl xyzf_QMixOrig = xyzf_QMixAl !------------------------------------------ ! 暖かい雨のパラメタリゼーション. ! * 雲<-->雨 の変換を行う. ! ! Warm rain parameterization. ! * Conversion from cloud to rain. !これまでの値を作業配列に保管 ! Previous values are stored to work area. ! xyzf_QMixWork = xyzf_QMixAl !雨への変化量を計算 ! Conversion values are calculated. ! xyzf_QMixAll = max( 1.0d-60, xyzf_QMixAl + xyzf_QMixBZ ) xyzf_Cloud2Rain = 0.0d0 do s = 1, CloudNum xyz_AutoConv = 0.0d0 xyz_Collect = 0.0d0 !併合成長 ! xyz_AutoConv = & & DelTime / AutoConvTime & & * max( 1.0d-60, ( xyzf_QMixAll(:,:,:,IdxC(s)) - QMixCr) ) !衝突合体成長 ! xyz_Collect = & & DelTime & & * 2.2d0 * FactorJ * xyzf_QMixAll(:,:,:,IdxC(s)) & & * (xyzf_QMixAll(:,:,:,IdxR(s)) * xyz_DensBZ) ** 0.875d0 !雲の変換量: 併合成長と合体衝突の和 ! 元々の変化量を上限値として設定する. 負の値となる. ! xyzf_Cloud2Rain(:,:,:,IdxC(s)) = & & - min( xyzf_QMixAll(:,:,:,IdxC(s)), ( xyz_AutoConv + xyz_Collect ) ) !雨の変換量. 符号は雲の変換量とは反対. xyzf_Cloud2Rain(:,:,:,IdxR(s)) = - xyzf_Cloud2Rain(:,:,:,IdxC(s)) end do ! 変化量を足し込む ! xyzf_QMixAl = xyzf_QMixWork + xyzf_Cloud2Rain ! Set Margin ! call SetMargin_xyzf(xyzf_QMixAl) !------------------------------------------- ! 湿潤飽和調節 ! * 蒸気<-->雲の変換を行う. ! ! Moist adjustment. ! * Conversion from vapor to cloud. ! call MoistAdjustSvapPress( & & xyz_ExnerNl, & ! (in) & xyz_PTempAl, & ! (inout) & xyzf_QMixAl & ! (inout) & ) if (IdxNH4SHr /= 0) then call MoistAdjustNH4SH( & & xyz_ExnerNl, & !(in) & xyz_PTempAl, & !(inout) & xyzf_QMixAl & !(inout) & ) end if ! Set Margin ! call SetMargin_xyz(xyz_PTempAl) call SetMargin_xyzf(xyzf_QMixAl) !------------------------------------------- ! 暖かい雨のパラメタリゼーション. ! * 蒸気<-->雨 の変換を行う ! ! Warm rain parameterization. ! * Conversion from rain to vapor. !これまでの値を作業配列に保管 ! Previous values are stored to work area. ! xyz_PTempWork = xyz_PTempAl xyzf_QMixWork = xyzf_QMixAl ! 雨から蒸気への混合比変化を求める ! * 温位の計算において, 混合比変化が必要となるため, ! 混合比変化を 1 つの配列として用意する. ! ! Conversion values are calculated. ! !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 ! xyz_ExnerAll = xyz_ExnerNl + xyz_ExnerBZ xyz_TempAll = ( xyz_PTempAl + xyz_PTempBZ ) * ( xyz_ExnerNl + xyz_ExnerBZ ) xyz_PressAll = PressBasis * ((xyz_ExnerNl + xyz_ExnerBZ) ** (CpDry / GasRDry)) xyzf_QMixAll = max( 1.0d-60, xyzf_QMixAl + xyzf_QMixBZ ) xyzf_Rain2Gas = 0.0d0 xyzf_Rain2GasNH4SH = 0.0d0 xyz_NonSaturate = 0.0d0 xyzf_DelPTemp = 0.0d0 do s = 1, CondNum !飽和蒸気圧と混合比の差(飽和度)を計算. ! 雨から蒸気への変換量は飽和度に比例する. ! xyz_NonSaturate = & & max( & & 1.0d-60, & & xyz_SvapPress(SpcWetID(IdxCC(s)), xyz_TempAll) & & * MolWtWet(IdxCG(s)) / ( MolWtDry * xyz_PressAll) & & - xyzf_QMixAll(:,:,:,IdxCG(s)) & & ) !雨の変換量 ! 元々の雨粒の混合比以上に蒸発が生じないように上限値を設定 ! xyzf_Rain2Gas(:,:,:,IdxCR(s)) = & & - min( & & DelTime * 4.85d-2 * FactorJ * xyz_NonSaturate & & * ( xyzf_QMixAll(:,:,:,IdxCR(s)) * xyz_DensBZ )** 0.65d0, & & xyzf_QMixAll(:,:,:,IdxCR(s)) & & ) !蒸気の変換量 ! 雨粒の変換量とは符号が逆となる ! xyzf_Rain2Gas(:,:,:,IdxCG(s)) = - xyzf_Rain2Gas(:,:,:,IdxCR(s)) ! xyzf_DelQMix を元に潜熱を計算 ! xyzf_DelPTemp(:,:,:,s) = & & xyz_LatentHeat( SpcWetID(IdxCR(s)), xyz_TempAll ) & & * xyzf_Rain2Gas(:,:,:,IdxCR(s)) & & / (xyz_ExnerAll * CpDry) end do !飽和蒸気圧と混合比の差(飽和度)を計算. ! 雨から蒸気への変換量は飽和度に比例する. ! 未飽和度を求めたいので, マイナスをかけ算している ! (DelQMixNH4SH は, NH4SH が増加する方向, すなわち飽和度を正としている) ! xyz_NonSaturate = 0.0d0 xyzf_Rain2GasNH4SH = 0.0d0 xyz_DelPTempNH4SH = 0.0d0 if (IdxNH4SHr /= 0) then xyz_NonSaturate = & & max( & & 1.0d-60, & & - xyz_DelQMixNH4SH( & & xyz_TempAll, xyz_PressAll, & & xyzf_QMixAll(:,:,:,IdxNH3), xyzf_QMixAll(:,:,:,IdxH2S), & & MolWtWet(IdxNH3), MolWtWet(IdxH2S) & & ) & & ) !雨の変換量 ! 元々の雨粒の混合比以上に蒸発が生じないように上限値を設定 ! xyzf_Rain2GasNH4SH(:,:,:,IdxNH4SHr) = & & - min( & & DelTime * 4.85d-2 * FactorJ * xyz_NonSaturate & & * (xyzf_QMixAll(:,:,:,IdxNH4SHr) * xyz_DensBZ) ** 0.65d0, & & xyzf_QMixAll(:,:,:,IdxNH4SHr) & & ) !蒸気の変換量 ! 雨粒の変換量とは符号が逆となる ! xyzf_Rain2GasNH4SH(:,:,:,IdxNH3) = & & - xyzf_Rain2GasNH4SH(:,:,:,IdxNH4SHr) * MolWtWet(IdxNH3) & & / MolWtWet(IdxNH4SHr) xyzf_Rain2GasNH4SH(:,:,:,IdxH2S) = & & - xyzf_Rain2GasNH4SH(:,:,:,IdxNH4SHr) * MolWtWet(IdxH2S) & & / MolWtWet(IdxNH4SHr) xyz_DelPTempNH4SH & & = ReactHeatNH4SH * xyzf_Rain2GasNH4SH(:,:,:,IdxNH4SHr) & & / (xyz_ExnerAll * CpDry) end if !変化量を足し算 ! xyzf_DelQMix = xyzf_Rain2Gas + xyzf_Rain2GasNH4SH xyz_DelPTemp = sum(xyzf_DelPTemp, 4) + xyz_DelPTempNH4SH ! 温位と混合比の計算. 雨から蒸気への変換分を追加 ! xyz_PTempAl = xyz_PTempWork + xyz_DelPTemp xyzf_QMixAl = xyzf_QMixWork + xyzf_DelQMix !------------------------------------------ ! Output ! xyz_Del = (xyz_PTempAl - xyz_PTempOrig) / DelTime xyzf_Del = (xyzf_QMixAl - xyzf_QMixOrig) / DelTime call HistoryAutoPut(TimeN, 'PTempCond', xyz_Del(1:nx, 1:ny, 1:nz) / DelTime) do l = 1, ncmax call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_Cond', xyzf_Del(1:nx, 1:ny, 1:nz, l) / DelTime) end do ! Set Margin ! call SetMargin_xyz(xyz_PTempAl) call SetMargin_xyzf(xyzf_QMixAl) end subroutine Cloudphys_K1969_forcing !!!=================================================================================!!! subroutine CloudPhys_K1969_FallRain(xyzf_QMix, xyzf_DQMixDt ) ! ! 雨粒の落下による移流を求める. ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax) !蒸気混合比(擾乱) real(DP), intent(inout) :: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax) !蒸気混合比の変化量 real(DP) :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax) !蒸気混合比(擾乱 + 平均場) real(DP) :: xyzf_FallRain(imin:imax,jmin:jmax,kmin:kmax, ncmax) !雨粒の落下効果 real(DP) :: xyz_VelZRain(imin:imax,jmin:jmax,kmin:kmax) !雨粒落下速度 real(DP) :: xyrf_FallRainFlux(imin:imax,jmin:jmax,kmin:kmax, ncmax) !雨粒落下フラックス real(DP) :: xyzf_DQMixDtOrig(imin:imax,jmin:jmax,kmin:kmax, ncmax) !蒸気混合比の変化量 integer :: s, l xyzf_QMixAll = max( 0.0d0, xyzf_QMix + xyzf_QMixBZ ) xyzf_DQMixDtOrig = xyzf_DQMixDt xyrf_FallRainFlux = 0.0d0 do s = 1, RainNum ! 雨粒終端速度 xyz_VelZRain = - 12.2d0 * FactorJ & * ( xyzf_QMixAll(:,:,:,IdxR(s)) ** 0.125d0 ) xyrf_FallRainFlux(:,:,:,IdxR(s)) = & & xyr_avr_xyz( & & xyz_DensBZ & & * xyzf_QMixAll(:,:,:,IdxR(s)) & & * xyz_VelZRain & & ) end do ! 上端のフラックスはゼロ do s = 1, RainNum xyrf_FallRainFlux(:,:,nz,IdxR(s)) = 0.0d0 end do ! 雨粒落下による時間変化率 xyzf_FallRain = 0.0d0 do s = 1, RainNum xyzf_FallRain(:,:,:,IdxR(s)) = & & - xyz_dz_xyr( xyrf_FallRainFlux(:,:,:,IdxR(s)) ) & & / xyz_DensBZ end do xyzf_DQMixDt = xyzf_DQMixDtOrig + xyzf_FallRain do l = 1, ncmax call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_Fall', xyzf_FallRain(1:nx, 1:ny, 1:nz, l)) call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_FallFluxAtLB', xyrf_FallRainFlux(1:nx, 1:ny, 0, l)) end do ! SetMargin ! call SetMargin_xyzf(xyzf_DQMixDt) end subroutine CloudPhys_K1969_FallRain end module Cloudphys_k1969