Class: cloudphys

暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める.

  * 中島健介 (1994) で利用した定式をそのまま利用.

Methods

Cloudphys_K1969_Init Cloudphys_K1969_forcing FactorCloud2Gas

Included Modules

dc_types dc_iounit dc_message gtool_historyauto mpi_wrapper timeset gridset constants basicset composition axesset xyz_deriv_module ChemCalc MoistAdjust namelist_util DExnerDt SetMargin

Public Instance methods

Cloudphys_K1969_Init( )

Subroutine :

This procedure input/output NAMELIST#cloudphys_k1969_nml .

[Source]

  subroutine Cloudphys_K1969_Init

    !暗黙の型宣言禁止
    implicit none

    !内部変数
    integer  :: unit    !装置番号
    integer  :: l
    character(STRING) :: Planet = ""
    character(*), parameter:: module_name = 'Cloudphys_K1969_Init'

    !-----------------------------------------------------------
    ! NAMELIST から情報を取得
    !
    NAMELIST /cloudphys_k1969_nml/ Planet, FactorJ, AutoConvTime, QMixCr, FlagDExnerDtCloud, FlagDExnerDtFall, FactorFallRain, FactorCloud2Rain, FactorRain2Gas, FactorCloud2Gas

    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", module_name, "Planet = %c",  c1=trim(Planet))
      call MessageNotify( "M", module_name, "FactorJ = %f",  d=(/FactorJ/) )
      call MessageNotify( "M", module_name, "AutoConvTime = %f",  d=(/AutoConvTime/) )
      call MessageNotify( "M", module_name, "QMixCr = %f",  d=(/QMixCr/) )
      call MessageNotify( "M", module_name, "FactorFallRain = %f",  d=(/FactorFallRain/) )
      call MessageNotify( "M", module_name, "FactorCloud2Rain = %f",  d=(/FactorCloud2Rain/) )
      call MessageNotify( "M", module_name, "FactorRain2Gas = %f",  d=(/FactorRain2Gas/) )
      call MessageNotify( "M", module_name, "FactorCloud2Gas = %f",  d=(/FactorCloud2Gas/) )
      call MessageNotify( "M", module_name, "FlagDExnerDtCloud= %b", l=(/ FlagDExnerDtCloud /))
      call MessageNotify( "M", module_name, "FlagDExnerDtFall= %b", l=(/ FlagDExnerDtFall /))
     end if

    !-----------------------------------------------------------
    ! ヒストリデータ定義
    !
    call HistoryAutoAddVariable( varname='ExnerFall', dims=(/'x','y','z','t'/), longname='Fall term of Exner function', units='kg.kg-1.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='PTempCond', dims=(/'x','y','z','t'/), longname='Latent heat term of potential temperature', units='K.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='ExnerCondTemp', dims=(/'x','y','z','t'/), longname='Latent heat term of exner function (Temp)', units='K.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='ExnerCondQMix', dims=(/'x','y','z','t'/), longname='Latent heat term of exner function (QMix)', 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


   do l = 1, CondNum
      call HistoryAutoAddVariable( varname=trim(SpcWetSymbol(l))//'_Hum', dims=(/'x','y','z','t'/), longname='Humidity of ' //trim(SpcWetSymbol(l)), units='1', xtype='float')
    end do
    
  end subroutine Cloudphys_K1969_Init

Cloudphys_K1969_forcing( xyz_ExnerNl, xyzf_QMixNl, xyz_DExnerDt, xyz_PTempAl, xyzf_QMixAl )

*Subroutine* :
`xyz_ExnerNl(imin:imax, jmin:jmax, kmin:kmax)` :	`real(DP), intent(in)`
`xyzf_QMixNl(imin:imax, jmin:jmax, kmin:kmax, ncmax)` :	`real(DP), intent(in)`
`xyz_DExnerDt(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), intent(inout)`

[Source]

  subroutine Cloudphys_K1969_forcing(xyz_ExnerNl, xyzf_QMixNl, xyz_DExnerDt, xyz_PTempAl, xyzf_QMixAl)

    implicit none

    real(DP), intent(in)    :: xyz_ExnerNl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)    :: xyzf_QMixNl(imin:imax, jmin:jmax, kmin:kmax, ncmax)
    real(DP), intent(inout) :: xyz_DExnerDt(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_PTempCond(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_QMixCond(imin:imax, jmin:jmax, kmin:kmax, ncmax)
    real(DP) :: DelTime
    integer  :: l, s
    integer  :: iG, iC, iR

    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)
    real(DP) :: xyz_PressDry(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_ExnerCondTemp(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_ExnerCondQMix(imin:imax,jmin:jmax,kmin:kmax)

    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)  :: xyz_ExnerFall(imin:imax,jmin:jmax,kmin:kmax) 

    real(DP)  :: xyz_QMixSat(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP)  :: xyz_QMixHum(imin:imax,jmin:jmax,kmin:kmax) 


    !-------------------------------------------------------------
    ! 時間刻み幅. Leap-frog なので, 2 \del t
    !
    DelTime = 2.0d0 * DelTimeLong


    !-------------------------------------------------------------
    ! 雨粒の落下による混合比の時間変化を計算 
    !   QMixNl を用いて QMixAl を更新する. 
    !

    ! 初期化
    ! 落下する分としては, Nl の値を利用する.
    !
    xyzf_QMixAll = max( 0.0d0, xyzf_QMixNl + xyzf_QMixBZ )
    xyrf_FallRainFlux = 0.0d0
    xyzf_FallRain = 0.0d0
    xyzf_QMixWork = xyzf_QMixAl + xyzf_QMixBZ

    do s = 1, RainNum
      iR = IdxR(s)

      ! 雨粒終端速度
      xyz_VelZRain = - 12.2d0 * FactorJ * ( xyzf_QMixAll(:,:,:,iR) ** 0.125d0 )

      ! フラックスの計算
      !
      xyrf_FallRainFlux(:,:,:,iR) = xyr_avr_xyz( xyz_DensBZ * xyz_VelZRain * xyzf_QMixAll(:,:,:,iR) )

      ! 上端のフラックスはゼロ
      !
      xyrf_FallRainFlux(:,:,nz,iR) = 0.0d0

      ! 雨粒落下による時間変化 (DelTime をかけてある)
      !
      xyzf_FallRain(:,:,:,iR) = - xyz_dz_xyr( xyrf_FallRainFlux(:,:,:,iR) ) / xyz_DensBZ * DelTime

      ! 元々の存在量より減ることは無いようにする. 
      !
      xyzf_FallRain(:,:,:,iR) = max( - xyzf_QMixWork(:,:,:,iR), xyzf_FallRain(:,:,:,iR) )
    end do

    ! 落下を考慮することで, Al の値を更新
    ! 雨の落下を切る場合には FactorFallRain = 0.0 とする. 
    !
    xyzf_QMixAl = xyzf_QMixAl + xyzf_FallRain * FactorFallRain

    ! 落下に伴う変化を Output 
    !
    if ( FlagDExnerDtFall ) then
      xyz_ExnerFall = xyz_DExnerDt_xyzf( xyzf_FallRain / DelTime ) !時間変動なので
    else
      xyz_ExnerFall = 0.0d0
    end if 

    ! 保管
    !
    xyz_DExnerDt  = xyz_DExnerDt + xyz_ExnerFall

    call HistoryAutoPut(TimeN, 'ExnerFall', xyz_ExnerFall(1:nx,1:ny,1:nz))    
    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


    !------------------------------------------
    ! 初期値を保管 Store Initial Value
    !
    xyz_PTempOrig = xyz_PTempAl
    xyzf_QMixOrig = xyzf_QMixAl

    ! 全エクスナー関数・全圧を計算. サブルーチン内では変化しない.
    !
    xyz_ExnerAll = xyz_ExnerNl + xyz_ExnerBZ
    xyz_PressAll = PressBasis * (xyz_ExnerAll ** (CpDry / GasRDry))

    !------------------------------------------    
    ! 暖かい雨のパラメタリゼーション.
    ! * 雲<-->雨 の変換を行う.
    !
    ! 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( 0.0d0, xyzf_QMixAl + xyzf_QMixBZ )
    xyzf_Cloud2Rain = 0.0d0

    do s = 1, CloudNum

      ! 値を保管
      !
      iC = IdxC(s)
      iR = IdxR(s)

      !併合成長
      !
      xyz_AutoConv = DelTime / AutoConvTime * max( 0.0d0, ( xyzf_QMixAll(:,:,:,iC) - QMixCr) )

      !衝突合体成長
      !
      xyz_Collect = DelTime * 2.2d0 * FactorJ * xyzf_QMixAll(:,:,:,iC) * (xyzf_QMixAll(:,:,:,iR) * xyz_DensBZ) ** 0.875d0  

      !雲の変換量: 併合成長と合体衝突の和
      !  元々の変化量を上限値として設定する. 負の値となる.
      !
      xyzf_Cloud2Rain(:,:,:,iC) = - min( xyzf_QMixAll(:,:,:,iC), ( xyz_AutoConv + xyz_Collect ) )
      
      !雨の変換量. 符号は雲の変換量とは反対. 
      xyzf_Cloud2Rain(:,:,:,iR) = - xyzf_Cloud2Rain(:,:,:,iC) 
    end do

    ! 変化量を足し込む
    ! 雲から雨へ変換させない場合は FactorCloud2Rain = 0.0 とする. 
    !
    xyzf_QMixAl = xyzf_QMixWork + xyzf_Cloud2Rain * FactorCloud2Rain


    !-------------------------------------------    
    ! 暖かい雨のパラメタリゼーション.
    ! * 蒸気<-->雨 の変換を行う
    !
    ! 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_TempAll   = ( xyz_PTempAl + xyz_PTempBZ ) * ( xyz_ExnerNl + xyz_ExnerBZ )
    xyzf_QMixAll  = max( 0.0d0, xyzf_QMixAl + xyzf_QMixBZ )
    xyz_PressDry  = xyz_PressDry_xyzf_xyz( xyzf_QMixAll, xyz_PressAll )

    xyzf_Rain2Gas = 0.0d0
    xyzf_DelPTemp = 0.0d0
    xyzf_Rain2GasNH4SH = 0.0d0
    xyz_DelPTempNH4SH  = 0.0d0

    do s = 1, CondNum

       ! 値を保管
       !
       iG = IdxCG(s)
       iC = IdxCC(s)
       iR = IdxCR(s)

      !飽和蒸気圧と混合比の差(飽和度)を計算. 
      !  雨から蒸気への変換量は飽和度に比例する.
      !
      xyz_NonSaturate = max( 0.0d0, xyz_SvapPress(SpcWetID(iC), xyz_TempAll) * MolWtWet(iG) / ( MolWtDry * xyz_PressDry) - xyzf_QMixAll(:,:,:,iG) )

      !雨の変換量
      !  元々の雨粒の混合比以上に蒸発が生じないように上限値を設定
      !
      xyzf_Rain2Gas(:,:,:,iR) = - min( DelTime * 4.85d-2 * FactorJ * xyz_NonSaturate * ( xyzf_QMixAll(:,:,:,iR) * xyz_DensBZ )** 0.65d0, xyzf_QMixAll(:,:,:,iR) ) 

      !蒸気の変換量
      !  雨粒の変換量とは符号が逆となる
      !
      xyzf_Rain2Gas(:,:,:,iG) = - xyzf_Rain2Gas(:,:,:,iR) 
    
      ! xyzf_DelQMix を元に潜熱を計算
      !
      xyzf_DelPTemp(:,:,:,s) = xyz_LatentHeat( SpcWetID(iR), xyz_TempAll ) * xyzf_Rain2Gas(:,:,:,iR) / (xyz_ExnerAll * CpDry) 

    end do

    !飽和蒸気圧と混合比の差(飽和度)を計算. 
    !  雨から蒸気への変換量は飽和度に比例する.
    !  未飽和度を求めたいので, マイナスをかけ算している
    !  (DelQMixNH4SH は, NH4SH が増加する方向, すなわち飽和度を正としている)
    !

    if (IdxNH4SHr /= 0) then 
      xyz_NonSaturate = max( 0.0d0, - xyz_DelQMixNH4SH( xyz_TempAll, xyz_PressAll, xyz_PressDry, 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

    !変化量を足し算
    !雨から蒸気への変換を切る場合は FactorRain2Gas = 0.0 とする. 
    !
    xyzf_DelQMix = xyzf_Rain2Gas + xyzf_Rain2GasNH4SH 
    xyz_DelPTemp = sum(xyzf_DelPTemp, 4) + xyz_DelPTempNH4SH 

    ! 温位と混合比の計算. 雨から蒸気への変換分を追加
    !
    xyz_PTempAl = xyz_PTempWork + xyz_DelPTemp * FactorRain2Gas
    xyzf_QMixAl = xyzf_QMixWork + xyzf_DelQMix * FactorRain2Gas


    !-------------------------------------------
    ! 湿潤飽和調節
    ! * 蒸気<-->雲の変換を行う.
    !
    ! Moist adjustment.
    ! * Conversion from vapor to cloud.
    !
    xyzf_QMixAll = max( 0.0d0, xyzf_QMixAl + xyzf_QMixBZ )
    xyz_PressDry = xyz_PressDry_xyzf_xyz( xyzf_QMixAll, xyz_PressAll )

    call MoistAdjustSvapPress( xyz_PressDry, xyz_ExnerNl, xyz_PTempAl, xyzf_QMixAl, FactorCloud2Gas )
    if (IdxNH4SHr /= 0) then 
      call MoistAdjustNH4SH( xyz_PressDry, xyz_ExnerNl, xyz_PTempAl, xyzf_QMixAl, FactorCloud2Gas )
    end if

    !------------------------------------------
    ! Output
    !
    xyz_PTempCond = (xyz_PTempAl - xyz_PTempOrig) / DelTime
    xyzf_QMixCond = (xyzf_QMixAl - xyzf_QMixOrig) / DelTime

    if ( FlagDExnerDtCloud ) then
      xyz_ExnerCondTemp = xyz_DExnerDt_xyz( xyz_PTempCond )
      xyz_ExnerCondQMix = xyz_DExnerDt_xyzf( xyzf_QMixCond )
    else
      xyz_ExnerCondTemp = 0.0d0
      xyz_ExnerCondQMix = 0.0d0
    end if

    xyz_DExnerDt  = xyz_DExnerDt + xyz_ExnerCondTemp + xyz_ExnerCondQMix

    call HistoryAutoPut(TimeN, 'PTempCond', xyz_PTempCond(1:nx, 1:ny, 1:nz))
    call HistoryAutoPut(TimeN, 'ExnerCondTemp', xyz_ExnerCondTemp(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'ExnerCondQMix', xyz_ExnerCondQMix(1:nx,1:ny,1:nz))
    do l = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_Cond', xyzf_QMixCond(1:nx, 1:ny, 1:nz, l))
    end do

    !----------------------------------------------------------------
    ! 飽和蒸気圧と平衡定数
    !----------------------------------------------------------------
    xyz_TempAll   = ( xyz_PTempAl + xyz_PTempBZ ) * ( xyz_ExnerNl + xyz_ExnerBZ )
    xyzf_QMixAll  = max( 0.0d0, xyzf_QMixAl + xyzf_QMixBZ )
    xyz_PressDry  = xyz_PressDry_xyzf_xyz( xyzf_QMixAll, xyz_PressAll )

    do s = 1, CondNum
      iG = IdxCG(s)
      iC = IdxCC(s)

      !飽和蒸気圧
      xyz_QMixSat = xyz_SvapPress(SpcWetID(iC), xyz_TempAll) * MolWtWet(iC) / MolWtDry / xyz_PressDry

      !相対湿度
      xyz_QMixHum = xyzf_QMixAll(:,:,:,iG) / xyz_QMixSat * 100.0d0

      !出力
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_Hum', xyz_QMixHum(1:nx, 1:ny, 1:nz))
    end do

    ! Set Margin
    !
    call SetMargin_xyz(xyz_PTempAl)
    call SetMargin_xyzf(xyzf_QMixAl)

  end subroutine Cloudphys_K1969_forcing

FactorCloud2Gas

Variable :

FactorCloud2Gas = 1.0d0 :

real(DP), save, public

:	雲から蒸気への変換の有無蒸気へ変換させない場合は値をゼロにする.