Class cloud_utils
In: radiation/cloud_utils.f90

雲関系ルーチン

Cloud-related routines

Note that Japanese and English are described in parallel.

雲の分布を設定.

In this module, the amount of cloud or cloud optical depth are set. This module is under development and is still a preliminary version.

Procedures List

!$ ! RadiationFluxDennouAGCM :放射フラックスの計算
!$ ! ———— :————
!$ ! RadiationFluxDennouAGCM :Calculate radiation flux

NAMELIST

NAMELIST#cloud_utils_nml

Methods

Included Modules

dc_types dc_message gridset gtool_historyauto timeset constants constants_snowseaice constants0 saturate dc_iounit namelist_util

Public Instance methods

Subroutine :
ParentRoutineName :character(*), intent(in)
:
!$ logical , intent(in) :FlagIncludeSnowPhaseChange
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in)
xyz_TempB(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
xyz_QH2OVapB(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
xyz_QH2OLiqB(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
xyz_QH2OSolB(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
xyz_QH2OVap(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
xyz_QH2OLiq(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
xyz_QH2OSol(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in)
xy_Rain(0:imax-1, 1:jmax) :real(DP), intent(in)
xy_Snow(0:imax-1, 1:jmax) :real(DP), intent(in)

[Source]

  subroutine CloudUtilConsChk( ParentRoutineName, xyr_Press, xyz_TempB, xyz_QH2OVapB, xyz_QH2OLiqB, xyz_QH2OSolB, xyz_Temp , xyz_QH2OVap , xyz_QH2OLiq , xyz_QH2OSol , xy_Rain, xy_Snow )


    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime            ! $ \Delta t $ [s]

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry, LatentHeat, LatentHeatFusion
                              ! $ L $ [J kg-1] .
                              ! 融解の潜熱.
                              ! Latent heat of fusion

    character(*), intent(in) :: ParentRoutineName
!!$    logical , intent(in) :: FlagIncludeSnowPhaseChange
    real(DP), intent(in) :: xyr_Press   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in) :: xyz_TempB   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_QH2OVapB(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_QH2OLiqB(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_QH2OSolB(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_Temp    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_QH2OVap (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_QH2OLiq (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_QH2OSol (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xy_Rain     (0:imax-1, 1:jmax)
    real(DP), intent(in) :: xy_Snow     (0:imax-1, 1:jmax)

    ! Local variables
    !
    real(DP) :: xyz_DelMass(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_Val(0:imax-1, 1:jmax)
    real(DP) :: xy_SumB(0:imax-1, 1:jmax)
    real(DP) :: xy_Sum(0:imax-1, 1:jmax)
    real(DP) :: xy_Ratio(0:imax-1, 1:jmax)
    integer  :: i
    integer  :: j
    integer  :: k


    do k = 1, kmax
      xyz_DelMass(:,:,k) = ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do

    xy_Sum = 0.0_DP
    do k = kmax, 1, -1
      xy_Val =   CpDry * xyz_TempB(:,:,k) + LatentHeat * xyz_QH2OVapB(:,:,k) - LatentHeatFusion * xyz_QH2OSolB(:,:,k)
      xy_Sum = xy_Sum + xy_Val * xyz_DelMass(:,:,k)
    end do

    xy_SumB = xy_Sum

    xy_Sum = 0.0_DP
    do k = kmax, 1, -1
      xy_Val =   CpDry * xyz_Temp (:,:,k) + LatentHeat * xyz_QH2OVap (:,:,k) - LatentHeatFusion * xyz_QH2OSol (:,:,k)
      xy_Sum = xy_Sum + xy_Val * xyz_DelMass(:,:,k)
    end do
!!$    if ( FlagIncludeSnowPhaseChange ) then
      xy_Sum = xy_Sum - LatentHeatFusion * xy_Snow * 2.0_DP * DelTime
!!$    end if

    xy_Ratio = ( xy_Sum - xy_SumB ) / ( xy_Sum + 1.0d-100 )
    do j = 1, jmax
      do i = 0, imax-1
        if ( abs( xy_Ratio(i,j) ) > 1.0d-10 ) then
          call MessageNotify( 'M', module_name, '%c, Modified condensate static energy is not conserved, %f.', c1 = trim(ParentRoutineName), d = (/ xy_Ratio(i,j) /) )
        end if
      end do
    end do



    xy_Sum = 0.0_DP
    do k = kmax, 1, -1
      xy_Val = xyz_QH2OVapB(:,:,k) + xyz_QH2OLiqB(:,:,k) + xyz_QH2OSolB(:,:,k)
      xy_Sum = xy_Sum + xy_Val * xyz_DelMass(:,:,k)
    end do

    xy_SumB = xy_Sum

    xy_Sum = 0.0_DP
    do k = kmax, 1, -1
      xy_Val = xyz_QH2OVap (:,:,k) + xyz_QH2OLiq (:,:,k) + xyz_QH2OSol (:,:,k)
      xy_Sum = xy_Sum + xy_Val * xyz_DelMass(:,:,k)
    end do
    xy_Sum = xy_Sum + ( xy_Rain + xy_Snow ) * 2.0_DP * DelTime

    xy_Ratio = ( xy_Sum - xy_SumB ) / ( xy_Sum + 1.0d-100 )
    do j = 1, jmax
      do i = 0, imax-1
        if ( abs( xy_Ratio(i,j) ) > 1.0d-10 ) then
          call MessageNotify( 'M', module_name, '%c, H2O mass is not conserved, %f.', c1 = trim(ParentRoutineName), d = (/ xy_Ratio(i,j) /) )
        end if
      end do
    end do


  end subroutine CloudUtilConsChk
Subroutine :
xyz_TransCloudOneLayer(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_CloudCover(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyrr_OverlappedCloudTrans(0:imax-1, 1:jmax, 0:kmax, 0:kmax) :real(DP), intent(out)

[Source]

  subroutine CloudUtilsCalcOverlapCloudTrans( xyz_TransCloudOneLayer, xyz_CloudCover, xyrr_OverlappedCloudTrans )

    ! USE statements
    !

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimesetClockStart, TimesetClockStop

!!$    use sort, only : SortQuick

    real(DP), intent(in ) :: xyz_TransCloudOneLayer   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_CloudCover           (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out) :: xyrr_OverlappedCloudTrans(0:imax-1, 1:jmax, 0:kmax, 0:kmax)


    real(DP) :: xyz_EffCloudCover           (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_CloudCoverSorted        (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_EffCloudCoverSorted     (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_TransCloudOneLayerSorted(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: CloudCoverSortedCur
    real(DP) :: EffCloudCoverSortedCur
    real(DP) :: TransCloudOneLayerSortedCur
    integer  :: KInsPos
    integer  :: i
    integer  :: j
    integer  :: k
    integer  :: kk
    integer  :: kkk



    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. cloud_utils_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! Cloud optical depth
    !

    select case ( IDCloudOverlapType )
    case ( IDCloudOverlapTypeRandom )

      xyz_EffCloudCover = xyz_CloudCover * ( 1.0_DP - xyz_TransCloudOneLayer )

      do k = 0, kmax
        kk = k
        xyrr_OverlappedCloudTrans(:,:,k,kk) = 1.0_DP
        do kk = k+1, kmax
          xyrr_OverlappedCloudTrans(:,:,k,kk) = xyrr_OverlappedCloudTrans(:,:,k,kk-1) * ( 1.0_DP - xyz_EffCloudCover(:,:,kk) )
        end do
      end do

      do k = 0, kmax
        do kk = 0, k-1
          xyrr_OverlappedCloudTrans(:,:,k,kk) = xyrr_OverlappedCloudTrans(:,:,kk,k)
        end do
      end do

    case ( IDCloudOverlapTypeMaxOverlap )

      ! see Chou et al. (2001)

      xyz_EffCloudCover = xyz_CloudCover * ( 1.0_DP - xyz_TransCloudOneLayer )


      ! Original method (computationally expensive, probably)
      !
!!$        do k = 0, kmax
!!$          kk = k
!!$          xyrr_OverlappedCloudTrans(:,:,k,kk) = 1.0_DP
!!$          do kk = k+1, kmax
!!$
!!$            xyz_CloudCoverSorted         = xyz_CloudCover
!!$            xyz_EffCloudCoverSorted      = xyz_EffCloudCover
!!$            xyz_TransCloudOneLayerSorted = xyz_TransCloudOneLayer
!!$
!!$            call SortQuick( imax, jmax, kk-k,             &
!!$              & xyz_CloudCoverSorted        (:,:,k+1:kk), &
!!$              & xyz_EffCloudCoverSorted     (:,:,k+1:kk), &
!!$              & xyz_TransCloudOneLayerSorted(:,:,k+1:kk)  &
!!$              & )
!!$
!!$            xyrr_OverlappedCloudTrans(:,:,k,kk) = 0.0_DP
!!$            do kkk = k+1, kk
!!$              xyrr_OverlappedCloudTrans(:,:,k,kk) =         &
!!$                & xyz_EffCloudCoverSorted(:,:,kkk)          &
!!$                & + xyrr_OverlappedCloudTrans(:,:,k,kk)     &
!!$                &   * xyz_TransCloudOneLayerSorted(:,:,kkk)
!!$            end do
!!$            xyrr_OverlappedCloudTrans(:,:,k,kk) = &
!!$              & 1.0_DP - xyrr_OverlappedCloudTrans(:,:,k,kk)
!!$
!!$          end do
!!$        end do


      ! Economical method (probably)
      !
      do k = 0, kmax

!!$          do kkk = 1, kmax
!!$              xyz_CloudCoverSorted(:,:,kkk) = real( kmax-kkk ) / real(kmax)
!!$!              xyz_CloudCoverSorted(:,:,kkk) = abs( 0.55d0 - real( kmax-kkk ) / real(kmax) )
!!$          end do
!!$          ! debug output
!!$          if ( k == 0 ) then
!!$            kk = kmax
!!$            do kkk = k+1, kk
!!$              write( 6, * ) kkk, xyz_CloudCoverSorted(0,jmax/2+1,kkk)
!!$            end do
!!$          end if

        xyz_CloudCoverSorted         = xyz_CloudCover
        xyz_EffCloudCoverSorted      = xyz_EffCloudCover
        xyz_TransCloudOneLayerSorted = xyz_TransCloudOneLayer

        kk = k
        xyrr_OverlappedCloudTrans(:,:,k,kk) = 1.0_DP
        do kk = k+1, kmax


          do j = 1, jmax
            do i = 0, imax-1

              ! xyz_CloudCoverSorted(i,j,kk) is inserved in an appropriate position.
              !
              KInsPos = kk
              loop : do kkk = k+1, kk-1

                if ( xyz_CloudCoverSorted(i,j,kk) < xyz_CloudCoverSorted(i,j,kkk) ) then
                  KInsPos = kkk
                  exit loop
                end if

              end do loop

              ! values are saved
              CloudCoverSortedCur         = xyz_CloudCoverSorted        (i,j,kk)
              EffCloudCoverSortedCur      = xyz_EffCloudCoverSorted     (i,j,kk)
              TransCloudOneLayerSortedCur = xyz_TransCloudOneLayerSorted(i,j,kk)

              ! values are shifted upward to empty an array at insert position
              do kkk = kk, KInsPos+1, -1
                xyz_CloudCoverSorted        (i,j,kkk) = xyz_CloudCoverSorted        (i,j,kkk-1)
                xyz_EffCloudCoverSorted     (i,j,kkk) = xyz_EffCloudCoverSorted     (i,j,kkk-1)
                xyz_TransCloudOneLayerSorted(i,j,kkk) = xyz_TransCloudOneLayerSorted(i,j,kkk-1)
              end do
              kkk = KInsPos
              xyz_CloudCoverSorted        (i,j,kkk) = CloudCoverSortedCur
              xyz_EffCloudCoverSorted     (i,j,kkk) = EffCloudCoverSortedCur
              xyz_TransCloudOneLayerSorted(i,j,kkk) = TransCloudOneLayerSortedCur

            end do
          end do


!!$            xyz_CloudCoverSorted         = xyz_CloudCover
!!$            do kkk = 1, kmax
!!$              xyz_CloudCoverSorted(:,:,kkk) = real( kmax-kkk ) / real(kmax)
!!$            end do
!!$            xyz_EffCloudCoverSorted      = xyz_EffCloudCover
!!$            xyz_TransCloudOneLayerSorted = xyz_TransCloudOneLayer
!!$
!!$            call SortQuick( imax, jmax, kk-k,             &
!!$              & xyz_CloudCoverSorted        (:,:,k+1:kk), &
!!$              & xyz_EffCloudCoverSorted     (:,:,k+1:kk), &
!!$              & xyz_TransCloudOneLayerSorted(:,:,k+1:kk)  &
!!$              & )


!!$            ! debug output
!!$            if ( ( k == 0 ) .and. ( kk == kmax-2 ) ) then
!!$              do kkk = k+1, kk
!!$                write( 6, * ) kkk, xyz_CloudCoverSorted(0,jmax/2+1,kkk)
!!$              end do
!!$              write( 6, * ) '-----'
!!$            end if


          xyrr_OverlappedCloudTrans(:,:,k,kk) = 0.0_DP
          do kkk = k+1, kk
            xyrr_OverlappedCloudTrans(:,:,k,kk) = xyz_EffCloudCoverSorted(:,:,kkk) + xyrr_OverlappedCloudTrans(:,:,k,kk) * xyz_TransCloudOneLayerSorted(:,:,kkk)
          end do
          xyrr_OverlappedCloudTrans(:,:,k,kk) = 1.0_DP - xyrr_OverlappedCloudTrans(:,:,k,kk)

        end do
      end do



      do k = 0, kmax
        do kk = 0, k-1
          xyrr_OverlappedCloudTrans(:,:,k,kk) = xyrr_OverlappedCloudTrans(:,:,kk,k)
        end do
      end do


    end select

    ! Output effective cloud cover
    !
!!$    call HistoryAutoPut( TimeN, 'EffCloudCover', &
!!$      & 1.0_DP - xyrr_OverlappedCloudTrans(:,:,0,kmax) )


  end subroutine CloudUtilsCalcOverlapCloudTrans
Subroutine :
ArgFlagSnow :logical, intent(in)

This procedure input/output NAMELIST#cloud_utils_nml .

[Source]

  subroutine CloudUtilsInit( ArgFlagSnow )

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable

    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
    use saturate, only: SaturateInit


    ! 宣言文 ; Declaration statements
    !

    logical, intent(in) :: ArgFlagSnow


    character(STRING) :: CloudOverlapType

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号.
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT.
                              ! IOSTAT of NAMELIST read

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /cloud_utils_nml/ CloudOverlapType, CCNMixRatPerUnitMass
          !
          ! デフォルト値については初期化手続 "cloud_utils#CloudUtilsInit"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "cloud_utils#CloudUtilsInit" for the default values.
          !

    ! 実行文 ; Executable statement
    !

    if ( cloud_utils_inited ) return


    FlagSnow = ArgFlagSnow


    ! デフォルト値の設定
    ! Default values settings
    !

    CloudOverlapType    = "Random"
!!$    CloudOverlapType    = "MaxOverlap"

    CCNMixRatPerUnitMass = 1.0e8_DP


    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = cloud_utils_nml, iostat = iostat_nml )             ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if


    select case ( CloudOverlapType )
    case ( 'Random' )
      IDCloudOverlapType = IDCloudOverlapTypeRandom
    case ( 'MaxOverlap' )
      IDCloudOverlapType = IDCloudOverlapTypeMaxOverlap
    case default
      call MessageNotify( 'E', module_name, 'CloudOverlapType=<%c> is not supported.', c1 = trim(CloudOverlapType) )
    end select


    ! Initialization of modules used in this module
    !

    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
    call SaturateInit


    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
!!$    call HistoryAutoAddVariable( 'EffCloudCover', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'effective cloud cover', '1' )



    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'CloudOverlapType     = %c', c1 = trim(CloudOverlapType) )
    call MessageNotify( 'M', module_name, 'CCNMixRatPerUnitMass = %f', d = (/ CCNMixRatPerUnitMass /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    cloud_utils_inited = .true.

  end subroutine CloudUtilsInit
Subroutine :
xyz_CloudCover(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_DelCloudOptDep(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(inout)

[Source]

  subroutine CloudUtilsLocalizeCloud( xyz_CloudCover, xyz_DelCloudOptDep )

    ! USE statements
    !

    real(DP), intent(in   ) :: xyz_CloudCover    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout) :: xyz_DelCloudOptDep(0:imax-1, 1:jmax, 1:kmax)


    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. cloud_utils_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! Cloud optical depth is scaled by considering cloud cover less than 1. 

    xyz_DelCloudOptDep = xyz_DelCloudOptDep / max( xyz_CloudCover, 1.0d-3 )


  end subroutine CloudUtilsLocalizeCloud
Subroutine :
Press :real(DP), intent(in )
PressLI :real(DP), intent(in )
PressUI :real(DP), intent(in )
PRCPArea :real(DP), intent(in )
PRCPEvapArea :real(DP), intent(in )
Temp :real(DP), intent(inout)
QH2OVap :real(DP), intent(inout)
SurfRainFlux :real(DP), intent(inout)
SurfSnowFlux :real(DP), intent(inout)

[Source]

  subroutine CloudUtilsPRCPEvap1Grid( Press, PressLI, PressUI, PRCPArea, PRCPEvapArea, Temp, QH2OVap, SurfRainFlux, SurfSnowFlux )

    ! USE statements
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime            ! $ \Delta t $ [s]

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav
                              ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration

    real(DP), intent(in   ) :: Press
    real(DP), intent(in   ) :: PressLI
    real(DP), intent(in   ) :: PressUI
    real(DP), intent(in   ) :: PRCPArea
    real(DP), intent(in   ) :: PRCPEvapArea
    real(DP), intent(inout) :: Temp
    real(DP), intent(inout) :: QH2OVap
    real(DP), intent(inout) :: SurfRainFlux
    real(DP), intent(inout) :: SurfSnowFlux


    ! Local variables
    !
    real(DP) :: DelTemp
    real(DP) :: DelQH2OVap

    real(DP) :: DelMass

    real(DP) :: PRCPFlux
    real(DP) :: DelPRCPFlux
!!$    real(DP) :: DelQH2OVap
    character(STRING) :: CharPhase


    integer :: l


    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. cloud_utils_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    DelMass = ( PressLI - PressUI ) / Grav


    do l = 1, 2

      select case ( l )
      case ( 1 ) ! liquid
        CharPhase = 'liquid'
        PRCPFlux = SurfRainFlux
      case ( 2 ) ! solid
        CharPhase = 'solid'
        PRCPFlux = SurfSnowFlux
      case default
        call MessageNotify( 'E', module_name, 'Unexpected number for select case.' )
      end select

      call CloudUtilsPRCPEvap1GridCore( ( 2.0_DP * DelTime ), CharPhase, DelMass, Press, Temp, QH2OVap, PRCPFlux, PRCPArea, PRCPEvapArea, DelPRCPFlux, DelTemp, DelQH2OVap )

      select case ( l )
      case ( 1 ) ! liquid
        SurfRainFlux = PRCPFlux + DelPRCPFlux
      case ( 2 ) ! solid
        SurfSnowFlux = PRCPFlux + DelPRCPFlux
      case default
        call MessageNotify( 'E', module_name, 'Unexpected number for select case.' )
      end select
      QH2OVap    = QH2OVap + DelQH2OVap
      Temp       = Temp    + DelTemp
    end do


  end subroutine CloudUtilsPRCPEvap1Grid
Subroutine :
PressLI :real(DP), intent(in )
PressUI :real(DP), intent(in )
Temp :real(DP), intent(inout)
SurfRainFlux :real(DP), intent(inout)
SurfSnowFlux :real(DP), intent(inout)

[Source]

  subroutine CloudUtilsPRCPStepPC1Grid( PressLI, PressUI, Temp, SurfRainFlux, SurfSnowFlux )

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime            ! $ \Delta t $ [s]

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: CpDry, Grav, LatentHeatFusion
                              ! $ L $ [J kg-1] .
                              ! 融解の潜熱.
                              ! Latent heat of fusion

    ! 雪と海氷の定数の設定
    ! Setting constants of snow and sea ice
    !
    use constants_snowseaice, only: TempCondWater



    real(DP), intent(in   ) :: PressLI
    real(DP), intent(in   ) :: PressUI
    real(DP), intent(inout) :: Temp
    real(DP), intent(inout) :: SurfRainFlux
    real(DP), intent(inout) :: SurfSnowFlux


    ! 作業変数
    ! Work variables
    !
    real(DP) :: DelMass
    real(DP) :: MassMaxFreezeRate
    real(DP) :: MassFreezeRate
    real(DP) :: MassMaxMeltRate
    real(DP) :: MassMeltRate


    ! 初期化確認
    ! Initialization check
    !
    if ( .not. cloud_utils_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    DelMass = ( PressLI - PressUI ) / Grav

    ! Freezing and melting switching at temperature of TempCondWater

    MassMaxFreezeRate = CpDry * ( TempCondWater - Temp ) * DelMass / LatentHeatFusion / ( 2.0_DP * DelTime )
    if ( MassMaxFreezeRate >= 0.0_DP ) then
      ! freezing
      if ( SurfRainFlux >= MassMaxFreezeRate ) then
        MassFreezeRate = MassMaxFreezeRate
      else
        MassFreezeRate = SurfRainFlux
      end if
      SurfRainFlux = SurfRainFlux - MassFreezeRate
      SurfSnowFlux = SurfSnowFlux + MassFreezeRate
      Temp = Temp + LatentHeatFusion * MassFreezeRate * 2.0_DP * DelTime / ( CpDry * DelMass )
    else
      ! melting
      MassMaxMeltRate = - MassMaxFreezeRate
      if ( SurfSnowFlux >= MassMaxMeltRate ) then
        MassMeltRate = MassMaxMeltRate
      else
        MassMeltRate = SurfSnowFlux
      end if
      SurfRainFlux = SurfRainFlux + MassMeltRate
      SurfSnowFlux = SurfSnowFlux - MassMeltRate
      Temp = Temp - LatentHeatFusion * MassMeltRate * 2.0_DP * DelTime / ( CpDry * DelMass )
    end if


  end subroutine CloudUtilsPRCPStepPC1Grid
Subroutine :
xyz_CloudCover(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_DelCloudOptDep(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(inout)

[Source]

  subroutine CloudUtilsSmearCloudOptDep( xyz_CloudCover, xyz_DelCloudOptDep )

    ! USE statements
    !

    real(DP), intent(in   ) :: xyz_CloudCover    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout) :: xyz_DelCloudOptDep(0:imax-1, 1:jmax, 1:kmax)


    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. cloud_utils_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! Cloud optical depth is scaled by the way of Kiehl et al. (1994).

    xyz_DelCloudOptDep = xyz_DelCloudOptDep * xyz_CloudCover**1.5_DP


  end subroutine CloudUtilsSmearCloudOptDep

Private Instance methods

CCNMixRatPerUnitMass
Variable :
CCNMixRatPerUnitMass :real(DP), save
: number of CCN per atmospheric mass (kg-1)
                           CCN : Cloud Condensation Nuclei
Subroutine :
TimeStep :real(DP) , intent(in )
CharPhase :character(*), intent(in )
DelMass :real(DP) , intent(in )
Press :real(DP) , intent(in )
Temp :real(DP) , intent(in )
QH2OVap :real(DP) , intent(in )
PRCPFlux :real(DP) , intent(in )
PRCPArea :real(DP) , intent(in )
PRCPEvapArea :real(DP) , intent(in )
DelPRCPFlux :real(DP) , intent(out)
DelTemp :real(DP) , intent(out)
DelQH2OVap :real(DP) , intent(out)

[Source]

  subroutine CloudUtilsPRCPEvap1GridCore( TimeStep, CharPhase, DelMass, Press, Temp, QH2OVap, PRCPFlux, PRCPArea, PRCPEvapArea, DelPRCPFlux, DelTemp, DelQH2OVap )

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: PI                    ! $ \pi $ .
                              ! 円周率.  Circular constant

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry, GasRDry, LatentHeat, LatentHeatFusion, EpsV
                              ! $ \epsilon_v $ .
                              ! 水蒸気分子量比.
                              ! Molecular weight of water vapor

    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
    use saturate, only: CalcQVapSatOnLiq, CalcDQVapSatDTempOnLiq, CalcQVapSatOnSol, CalcDQVapSatDTempOnSol


    real(DP)    , intent(in ) :: TimeStep
    character(*), intent(in ) :: CharPhase
    real(DP)    , intent(in ) :: DelMass
    real(DP)    , intent(in ) :: Press
    real(DP)    , intent(in ) :: Temp
    real(DP)    , intent(in ) :: QH2OVap
    real(DP)    , intent(in ) :: PRCPFlux
    real(DP)    , intent(in ) :: PRCPArea
    real(DP)    , intent(in ) :: PRCPEvapArea
    real(DP)    , intent(out) :: DelPRCPFlux
    real(DP)    , intent(out) :: DelTemp
    real(DP)    , intent(out) :: DelQH2OVap


    ! Parameters for evaporation of rain
    real(DP), parameter :: DensWater            = 1.0d3
    !                            rho_w
    real(DP)            :: CCNND
    !                            number density of CCN, N0 or Nt (m-3)
    real(DP), parameter :: H2OVapDiffCoef       = 1.0d-5
    !                            Kd
    real(DP), parameter :: PRCPFallVel0         = 10.0_DP
    !                            m s-1

    real(DP) :: PRCPFallVel
    !                            m s-1
    real(DP) :: PRCPFallVelRatio

    real(DP) :: Dens
    !                           rho
    real(DP) :: DensPRCP
    !                           (rho q_r)

    real(DP) :: DelZ

    real(DP) :: FactorF
    real(DP) :: FactorG
    real(DP) :: FactorH
    real(DP) :: FactorI

    real(DP) :: LatentHeatSubl
    real(DP) :: LatentHeatLocal

    real(DP) :: VirTemp
    real(DP) :: QH2OVapSat
    real(DP) :: DQH2OVapSatDTemp

    real(DP) :: TempN
    real(DP) :: QH2OVapN
    real(DP) :: PRCPN
    real(DP) :: TempA
    real(DP) :: QH2OVapA
    real(DP) :: PRCPA

    real(DP) :: DelPRCP

    real(DP), parameter :: DelTempThreshold = 1.0e-2_DP
    integer, parameter :: ItrMax = 100
    real(DP) :: a_DelTemp(ItrMax)

    integer :: itr


    LatentHeatSubl = LatentHeat + LatentHeatFusion


    select case ( CharPhase )
    case ( 'liquid' )
      ! for liquid water
      PRCPFallVelRatio = 1.0_DP
      LatentHeatLocal  = LatentHeat
    case ( 'solid' )
      ! for solid water (ice)
      PRCPFallVelRatio = 0.5_DP
      LatentHeatLocal  = LatentHeatSubl
    case default
      call MessageNotify( 'E', module_name, 'Unexpected character for select case.' )
    end select


    VirTemp = Temp * ( 1.0_DP + ((( 1.0_DP / EpsV ) - 1.0_DP ) * QH2OVap) )
    Dens = Press / ( GasRDry * VirTemp )
    DelZ = DelMass / Dens

    PRCPFallVel = PRCPFallVel0 * PRCPFallVelRatio
!!$    DensPRCP = ( PRCPFlux / ( PRCPArea + 1.0e-10_DP ) ) / PRCPFallVel
    DensPRCP = ( max( PRCPFlux, 0.0_DP ) / ( PRCPArea + 1.0e-10_DP ) ) / PRCPFallVel

    ! cloud condensation
    CCNND = CCNMixRatPerUnitMass * Dens

    FactorF = Dens * H2OVapDiffCoef * PRCPEvapArea * ( 48.0_DP * ( PI * CCNND )**2 / DensWater * DensPRCP )**(1.0_DP/3.0_DP)
    FactorG = DelZ * TimeStep * FactorF
    FactorH = FactorG / DelMass
    FactorI = LatentHeatLocal / CpDry * FactorH

    TempN    = Temp
    QH2OVapN = QH2OVap
!!$    PRCPN    = PRCPFlux * TimeStep
    PRCPN    = max( PRCPFlux, 0.0_DP ) * TimeStep
    loop_evap : do itr = 1, ItrMax

      select case ( CharPhase )
      case ( 'liquid' )
        ! for liquid water
        QH2OVapSat       = CalcQVapSatOnLiq( TempN, Press )
        DQH2OVapSatDTemp = CalcDQVapSatDTempOnLiq( TempN, QH2OVapSat )
      case ( 'solid' )
        ! for solid water (ice)
        QH2OVapSat       = CalcQVapSatOnSol( TempN, Press )
        DQH2OVapSatDTemp = CalcDQVapSatDTempOnSol( TempN, QH2OVapSat )
      case default
        call MessageNotify( 'E', module_name, 'Unexpected character for select case.' )
      end select

      DelTemp = - FactorI * ( QH2OVapSat - QH2OVapN ) / ( 1.0_DP + FactorH + FactorI * DQH2OVapSatDTemp )

      DelQH2OVap = - CpDry * DelTemp / LatentHeatLocal
      DelPRCP    = - DelQH2OVap * DelMass

      if ( ( PRCPN + DelPRCP ) >= 0.0_DP ) then
        ! part of precipitation is evaporated
        ! nothing to do
      else
        ! all precipitation is evaporated
        DelPRCP    = - PRCPN
        DelQH2OVap = - DelPRCP / DelMass
        DelTemp    = - LatentHeatLocal * DelQH2OVap / CpDry
      end if


      if ( abs( DelTemp ) < DelTempThreshold ) exit loop_evap

      PRCPA    = PRCPN    + DelPRCP
      TempA    = TempN    + DelTemp
      QH2OVapA = QH2OVapN + DelQH2OVap

      PRCPN    = PRCPA
      TempN    = TempA
      QH2OVapN = QH2OVapA

      a_DelTemp(itr) = DelTemp
    end do loop_evap
    if ( itr > 100 ) then
      write( 6, * ) a_DelTemp
      call MessageNotify( 'E', module_name, 'Evaporation loop is not convergent, %d, %f.', i = (/ itr /), d = (/ DelTemp /) )
    end if

    DelPRCPFlux = DelPRCP / TimeStep


  end subroutine CloudUtilsPRCPEvap1GridCore
Subroutine :
CharPhase :character(*), intent(in )
DelMass :real(DP) , intent(in )
Press :real(DP) , intent(in )
QH2OVap :real(DP) , intent(in )
QH2OVapSat :real(DP) , intent(in )
VirTemp :real(DP) , intent(in )
PRCP :real(DP) , intent(in )
PRCPArea :real(DP) , intent(in )
PRCPEvapArea :real(DP) , intent(in )
DelPRCPFlux :real(DP) , intent(out)

[Source]

  subroutine CloudUtilsPRCPEvap1GridCoreExp( CharPhase, DelMass, Press, QH2OVap, QH2OVapSat, VirTemp, PRCP, PRCPArea, PRCPEvapArea, DelPRCPFlux )

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: PI                    ! $ \pi $ .
                              ! 円周率.  Circular constant

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, GasRDry
                              ! $ R $ [J kg-1 K-1].
                              ! 乾燥大気の気体定数.
                              ! Gas constant of air


    character(*), intent(in ) :: CharPhase
    real(DP)    , intent(in ) :: DelMass
    real(DP)    , intent(in ) :: Press
    real(DP)    , intent(in ) :: QH2OVap
    real(DP)    , intent(in ) :: QH2OVapSat
    real(DP)    , intent(in ) :: VirTemp
    real(DP)    , intent(in ) :: PRCP
    real(DP)    , intent(in ) :: PRCPArea
    real(DP)    , intent(in ) :: PRCPEvapArea
    real(DP)    , intent(out) :: DelPRCPFlux


    ! Parameters for evaporation of rain
    real(DP), parameter :: DensWater            = 1.0d3
    !                            rho_w
    !   Values below are from Kessler (1969)
    real(DP), parameter :: PRCPFallVelFactor0        = 130.0d0
    !                            K
    real(DP), parameter :: MedianDiameterFactor      = 3.67d0
    !                            C'
    real(DP), parameter :: PRCPDistFactor            = 1.0d7
    !                            N0
    real(DP), parameter :: PRCPEvapRatUnitDiamFactor = 2.24d3
    !                            C
    real(DP), parameter :: H2OVapDiffCoef            = 1.0d-5
    !                            Kd

    real(DP), parameter :: PRCPFallVelSimple0        = 10.0d0
    !                            m s-1

    real(DP) :: PRCPFallVelRatio
    real(DP) :: PRCPFallVelFactor

    real(DP) :: Dens0
    !                            rho_0
    real(DP) :: V00
    !                            V_{00}
    real(DP) :: PRCPEvapFactor

    real(DP) :: Dens
    !                           rho
    real(DP) :: DensPRCP
    !                           (rho q_r)
!!$    real(DP) :: RainArea
!!$    !                           a_p
!!$    real(DP) :: RainEvapArea
!!$    !                           A = max( a_p - a, 0 )
!!$    real(DP) :: xy_CloudCover   (0:imax-1, 1:jmax)
!!$    !                           a
    real(DP) :: PRCPEvapRate

    real(DP) :: DelZ


    select case ( CharPhase )
    case ( 'liquid' )
      ! for liquid water
      PRCPFallVelRatio = 1.0_DP
    case ( 'solid' )
      ! for solid water (ice)
      PRCPFallVelRatio = 0.5_DP
    case ( 'mixture' )
      ! for mixture, this is only for test
      PRCPFallVelRatio = 1.0_DP
    case default
      call MessageNotify( 'E', module_name, 'Unexpected character for select case.' )
    end select
    !

    ! Values for evaporation of rain
    Dens = Press / ( GasRDry * VirTemp )

    DelZ = DelMass / Dens


    if ( .false. ) then ! ECMWF version

      PRCPFallVelFactor = PRCPFallVelFactor0 * PRCPFallVelRatio

      ! Parameters for evaporation of rain
      Dens0 = 1013.0d2 / ( GasRDry * 300.0_DP )
      V00 = PRCPFallVelFactor * sqrt( MedianDiameterFactor ) / ( PI * DensWater * PRCPDistFactor )**(1.0d0/8.0d0)
      PRCPEvapFactor = PRCPEvapRatUnitDiamFactor * 1.429624558860304d0 * H2OVapDiffCoef * PRCPDistFactor**(7.0d0/20.0d0) / ( PI * DensWater )**(13.0d0/20.0d0)

!!$    RainArea   = RainArea
!!$    xy_CloudCover = CloudCover
!!$    xy_RainEvapArea = max( xy_RainArea - xy_CloudCover, 0.0_DP )
!!$    RainEvapArea = RainEvapArea

      DensPRCP = ( PRCP / ( PRCPArea + 1.0d-10 ) / ( V00 * sqrt( Dens0 / Dens ) ) )**(8.0d0/9.0d0)
      PRCPEvapRate = Dens * PRCPEvapArea * PRCPEvapFactor * max( QH2OVapSat - QH2OVap, 0.0_DP ) * DensPRCP**(13.0d0/20.0d0)

    else ! simple version

      V00 = PRCPFallVelSimple0 * PRCPFallVelRatio

      PRCPEvapFactor = ( 48.0_DP * ( PI * PRCPDistFactor )**2 / DensWater )**(1.0_DP/3.0_DP) * H2OVapDiffCoef

      DensPRCP = ( PRCP / ( PRCPArea + 1.0d-10 ) ) / V00
      PRCPEvapRate = Dens * PRCPEvapArea * PRCPEvapFactor * max( QH2OVapSat - QH2OVap, 0.0_DP ) * DensPRCP**(1.0_DP/3.0_DP)

    end if

    ! PRCPEvapRate (kg m-3 s-1)
    ! DelZ         (m)
    ! DelPRCPFlux  (kg m-2 s-1)
    DelPRCPFlux = PRCPEvapRate * DelZ

    DelPRCPFlux = min( DelPRCPFlux, PRCP )


  end subroutine CloudUtilsPRCPEvap1GridCoreExp
FlagSnow
Variable :
FlagSnow :logical , save
: A flag for snow
IDCloudOverlapType
Variable :
IDCloudOverlapType :integer , save
IDCloudOverlapTypeMaxOverlap
Constant :
IDCloudOverlapTypeMaxOverlap = 2 :integer , parameter
IDCloudOverlapTypeRandom
Constant :
IDCloudOverlapTypeRandom = 1 :integer , parameter
cloud_utils_inited
Variable :
cloud_utils_inited = .false. :logical, save
: 初期設定フラグ. Initialization flag
module_name
Constant :
module_name = ‘cloud_utils :character(*), parameter
: モジュールの名称. Module name
version
Constant :
version = ’$Name: $’ // ’$Id: cloud_utils.f90,v 1.7 2015/02/11 11:55:19 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version