Class | cloud_utils |
In: |
radiation/cloud_utils.f90
|
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.
!$ ! RadiationFluxDennouAGCM : | 放射フラックスの計算 |
!$ ! ———— : | ———— |
!$ ! RadiationFluxDennouAGCM : | Calculate radiation flux |
Subroutine : | |||||
ParentRoutineName : | character(*), 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) : | real(DP), intent(in) |
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) |
subroutine CloudUtilsCalcOverlapCloudTrans( xyz_TransCloudOneLayer, xyz_CloudCover, xyrr_OverlappedCloudTrans ) ! USE statements ! ! ヒストリデータ出力 ! History data output ! use gtool_historyauto, only: HistoryAutoPut ! 時刻管理 ! Time control ! use timeset, only: TimeN, EndTime, 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 .
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) |
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) |
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) :: QH2OVapSat 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) |
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) |
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
Variable : | |||
CCNMixRatPerUnitMass : | real(DP), save
|
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) |
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) :: PRCPFallVelRatio real(DP) :: PRCPFallVelFactor real(DP) :: PRCPFallVel ! m s-1 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) :: QH2OVapSatA 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 ! 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 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) |
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
Variable : | |||
cloud_utils_inited = .false. : | logical, save
|
Constant : | |||
version = ’$Name: dcpam5-20150217 $’ // ’$Id: cloud_utils.f90,v 1.7 2015/02/11 11:55:19 yot Exp $’ : | character(*), parameter
|