!= Module HeatFlux ! ! Authors:: ODAKA Masatsugu, TAKAHASHI Yoshiyuki ! Version:: $Id: surfaceflux_bulk.f90,v 1.13 2011-10-10 15:44:24 yot Exp $ ! Tag Name:: $Name: arare5-20111010 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! ! ! *** Explanation below are obsolete. (YOT, 2011/09/01) *** ! ! ! 下部境界からのフラックスによる温度と凝結成分の変化率をバルク方法に ! 基づいて計算するモジュール. これは中島 (1994) で用いられた方法である. ! ! 熱フラックスを Fh と凝結物質のフラックス Fq とすると, 温度と凝結物質 ! の変化率 H, Q は ! ! H = Fh/Δz_1 ! Q = Fq/Δz_1 ! ! と表される. ここで Δz_1 は最下層の格子間隔である. ! ! 熱フラックス Fh と凝結物質のフラックス Fq は以下の式にしたがって ! 計算する. ! ! Fh = - Cdρ|V| * (π_1θ_1 - T_sfc) ! Fh = - Cdρ|V| * (Q_1 - Q*(T_sfc)) ! ! ここで θ_1, Q_1, π_1 は最下層の温度と凝結物質の混合比および無次元 ! 圧力関数, T_sfc は下部境界の温度, Q*(T_sfc) は T_sfc で決まる飽和混 ! 合比である. ! ! バルク係数 Cd は一定とする. 無次元圧力関数 π_1 は基本場の値を用いる. ! 風速値 |V| は ! ! V = ( V^2 + V_0^2 )^(1/2) ! ! と計算する. ! !== Error Handling ! !== Bugs ! !== Note ! ! !== Future Plans ! ! module Surfaceflux_bulk ! !下部境界でのフラックスの計算モジュール ! !モジュール読み込み 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 gridset, only: imin, & !x 方向の配列の下限 & imax, & !x 方向の配列の上限 & jmin, & !y 方向の配列の下限 & jmax, & !y 方向の配列の上限 & kmin, & !z 方向の配列の下限 & kmax, & !z 方向の配列の上限 & nx, ny, nz, ncmax use axesset, only: z_dz, & !z 方向の格子点間隔 & xyz_avr_pyz, & & xyz_avr_xqz, & & pyz_avr_xyz, & & xqz_avr_xyz use basicset, only: xyz_ExnerBZ, & !エクスナー関数の基本場 & xyz_PressBZ, & ! & xyz_PTempBZ, & !温位の基本場 & xyz_TempBZ, & ! & xyzf_QMixBZ, & !温位の基本場 & xyz_DensBZ !基本場の密度 use constants,only: MolWtDry, PressBasis, TempSfc, PressSfc, CpDry, GasRDry use composition,only : IdxCC, IdxCG, SpcWetID, CondNum, MolWtWet, SpcWetSymbol use chemcalc,only : SvapPress use namelist_util, only: namelist_filename use timeset, only: TimeN use setmargin,only: SetMargin_xyz, SetMargin_xyzf, SetMargin_pyz, SetMargin_xqz !暗黙の型宣言禁止 implicit none !属性の指定 private !関数を public に設定 public surfaceflux_bulk_init public surfaceflux_bulk_forcing !変数定義 real(DP), save :: Bulk = 1.5d-3 !熱・運動量フラックスのバルク係数 real(DP), save :: Vel0 = 0.0d0 !下層での水平速度嵩上げ値 character(*), parameter:: module_name = 'surfaceflux_bulk' ! モジュールの名称. ! Module name contains !!!------------------------------------------------------------------------!!! subroutine Surfaceflux_Bulk_init ! !NAMELIST から必要な情報を読み取り, 時間関連の変数の設定を行う. ! !暗黙の型宣言禁止 implicit none !内部変数 integer :: l, unit !--------------------------------------------------------------- ! NAMELIST から情報を取得 ! NAMELIST /surfaceflux_bulk_nml/ Bulk, Vel0 call FileOpen(unit, file=namelist_filename, mode='r') read(unit, NML=surfaceflux_bulk_nml) close(unit) if (myrank == 0) then call MessageNotify( "M", module_name, "Bulk = %f", d=(/Bulk/) ) call MessageNotify( "M", module_name, "Vel0 = %f", d=(/Vel0/)) end if call HistoryAutoAddVariable( & & varname='PTempSfcFlux', & & dims=(/'x','y','t'/), & & longname='surface potential temperature flux (heat flux divided by density and specific heat)', & & units='K.m.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelXSfcFlux', & & dims=(/'x','y','t'/), & & longname='surface flux of x-component of velocity (momentum flux divided by density)', & & units='m2.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelYSfcFlux', & & dims=(/'x','y','t'/), & & longname='surface flux of y-component of velocity (momentum flux divided by density)', & & units='m2.s-2', & & xtype='float') do l = 1, ncmax call HistoryAutoAddVariable( & & varname=trim(SpcWetSymbol(l))//'_SfcFlux', & & dims=(/'x','y','t'/), & & longname='surface flux of ' & & //trim(SpcWetSymbol(l))//' mixing ratio (mass flux divided by density)', & & units='m.s-1', & & xtype='float') end do call HistoryAutoAddVariable( & & varname='PTempSfc', & & dims=(/'x','y','z','t'/), & & longname='potential temperature tendency by surface flux', & & units='K.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelXSfc', & & dims=(/'x','y','z','t'/), & & longname='x-component velocity tendency by surface flux', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelYSfc', & & dims=(/'x','y','z','t'/), & & longname='y-component velocity tendency by surface flux', & & units='m.s-2', & & xtype='float') do l = 1, ncmax call HistoryAutoAddVariable( & & varname=trim(SpcWetSymbol(l))//'_Sfc', & & dims=(/'x','y','z','t'/), & & longname=trim(SpcWetSymbol(l))//' mixing ratio tendency by surface flux', & & units='s-1', & & xtype='float') end do call HistoryAutoAddVariable( & & varname='SfcHeatFlux', & & dims=(/'x','y','t'/), & & longname='surface heat flux', & & units='W.m-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='SfcXMomFlux', & & dims=(/'x','y','t'/), & & longname='surface x-component momentum flux', & & units='kg.m-2.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='SfcYMomFlux', & & dims=(/'x','y','t'/), & & longname='surface y-component momentum flux', & & units='kg.m-2.s-1', & & xtype='float') do l = 1, ncmax call HistoryAutoAddVariable( & & varname=trim(SpcWetSymbol(l))//'_SfcMassFlux', & & dims=(/'x','y','t'/), & & longname=trim(SpcWetSymbol(l))//' surface mass flux', & & units='kg.m-2.s-1', & & xtype='float') end do end subroutine Surfaceflux_Bulk_init !!!------------------------------------------------------------------------!!! subroutine Surfaceflux_Bulk_forcing( & & pyz_VelX, xqz_VelY, xyz_PTemp, xyz_Exner, xyzf_QMix, & & pyz_DVelXDt, xqz_DVelYDt, xyz_DPTempDt, xyzf_DQMixDt & & ) ! ! 下部境界からのフラックスによる温度の変化率を, ! バルク方法に基づいて計算する. ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax) !水平風速 real(DP), intent(in) :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax) !水平風速 real(DP), intent(in) :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax) !温位の擾乱成分 real(DP), intent(in) :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax) !温位の擾乱成分 real(DP), intent(in) :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax) !温位の擾乱成分 real(DP), intent(inout):: pyz_DVelXDt(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout):: xqz_DVelYDt(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout):: xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout):: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax) real(DP) :: py_VelXflux (imin:imax,jmin:jmax) !運動量フラックス !(strictly speaking, this value is not !momentum flux, but is it divided by density) real(DP) :: xq_VelYflux (imin:imax,jmin:jmax) !運動量フラックス !(strictly speaking, this value is not !momentum flux, but is it divided by density) real(DP) :: xy_PTempFlux(imin:imax,jmin:jmax) !地表面熱フラックス !(strictly speaking, this value is not !heat flux, but is it divided by density and !specific heat) real(DP) :: xyf_QMixFlux(imin:imax,jmin:jmax,ncmax) !物質的フラックス real(DP) :: xyz_VelX(imin:imax,jmin:jmax,kmin:kmax) !水平風速 (xyz 格子) real(DP) :: xyz_VelY(imin:imax,jmin:jmax,kmin:kmax) !水平風速 (xyz 格子) real(DP) :: xyz_AbsVel(imin:imax,jmin:jmax,kmin:kmax) !水平風速 (xyz 格子) real(DP) :: pyz_AbsVel(imin:imax,jmin:jmax,kmin:kmax) !水平風速 (pyz 格子) real(DP) :: xqz_AbsVel(imin:imax,jmin:jmax,kmin:kmax) !水平風速 (xqz 格子) real(DP) :: xyz_PTempAll(imin:imax,jmin:jmax,kmin:kmax) !Total value of potential temperature real(DP) :: xyz_TempAll (imin:imax,jmin:jmax,kmin:kmax) !Total value of temperature real(DP) :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax) !Total value of mixing ratios real(DP) :: xy_DPTempDtBulk(imin:imax,jmin:jmax) !potential temperature tendency by surface flux real(DP) :: xyf_DQMixDtBulk(imin:imax,jmin:jmax, ncmax) !mixing ratio tendency by surface flux real(DP) :: py_DVelXDtBulk (imin:imax,jmin:jmax) !x-component velocity tendency by surface flux real(DP) :: xq_DVelYDtBulk (imin:imax,jmin:jmax) !y-component velocity tendency by surface flux real(DP) :: xyz_DPTempDtBulk(imin:imax,jmin:jmax,kmin:kmax) ! variable for output real(DP) :: xyzf_DQMixDtBulk(imin:imax,jmin:jmax,kmin:kmax, ncmax) ! variable for output real(DP) :: pyz_DVelXDtBulk (imin:imax,jmin:jmax,kmin:kmax) ! variable for output real(DP) :: xqz_DVelYDtBulk (imin:imax,jmin:jmax,kmin:kmax) ! variable for output real(DP) :: ExnerBZSfc ! Basic state Exner function at the surface real(DP) :: xy_PressSfc(imin:imax,jmin:jmax) ! Total pressure at the surface integer :: kz !配列添字 integer :: s !ループ変数 ! 初期化 ! kz = 1 xyz_PTempAll = xyz_PTemp + xyz_PTempBZ xyz_TempAll = (xyz_Exner + xyz_ExnerBZ) * xyz_PTempAll xyzf_QMixAll = xyzf_QMix + xyzf_QMixBZ ExnerBZSfc = (PressSfc / PressBasis) ** (GasRDry / CpDry) xy_PressSfc = PressBasis * ( ExnerBZSfc + xyz_Exner(:,:,kz) )**(CpDry / GasRDry) ! Perturbation component of Exner function at the surface is assumed ! to be same as that at the lowest layer. (YOT, 2011/09/03) ! Velocities at xyz grid points are calculated. xyz_VelX = xyz_avr_pyz(pyz_VelX) xyz_VelY = xyz_avr_xqz(xqz_VelY) xyz_AbsVel = SQRT( xyz_VelX**2 + xyz_VelY**2 + Vel0**2 ) pyz_AbsVel = pyz_avr_xyz(xyz_AbsVel) xqz_AbsVel = xqz_avr_xyz(xyz_AbsVel) ! Something like heat, mass, and momentum fluxes are calculated. ! The values below are not heat, mass, and momentum fluxes, but are those divided by ! by density and specific heat, density, and density, respectively. ! xy_PTempFlux = - Bulk * xyz_AbsVel(:,:,kz) & & * ( xyz_PTempAll(:,:,kz) - TempSfc / ( ExnerBZSfc + xyz_Exner(:,:,kz) ) ) ! xyf_QMixFlux = 0.0d0 do s = 1, CondNum xyf_QMixFlux(:,:,IdxCG(s)) = & & - Bulk * xyz_AbsVel(:,:,kz) & & * ( & & xyzf_QMixAll(:,:,kz,s) & & - SvapPress( SpcWetID(IdxCC(s)), TempSfc ) & & / ( xy_PressSfc ) & & * (MolWtWet(IdxCG(s)) / MolWtDry) & & ) end do ! py_VelXFlux = - Bulk * pyz_AbsVel(:,:,kz) * pyz_VelX(:,:,kz) ! xq_VelYFlux = - Bulk * xqz_AbsVel(:,:,kz) * xqz_VelY(:,:,kz) ! Something like heat flux and mass flux are restricted. ! This would be arbitrary treatment. xy_PTempFlux = max( 0.0d0, xy_PTempFlux ) xyf_QMixFlux = max( 0.0d0, xyf_QMixFlux ) ! Tendencies by surface fluxes (convergences of fluxes) are calculated. ! xy_DPTempDtBulk = - ( 0.0d0 - xy_PTempFlux ) / z_dz(kz) ! xyf_DQMixDtBulk = - ( 0.0d0 - xyf_QMixFlux ) / z_dz(kz) ! py_DVelXDtBulk = - ( 0.0d0 - py_VelXFlux ) / z_dz(kz) ! xq_DVelYDtBulk = - ( 0.0d0 - xq_VelYFlux ) / z_dz(kz) ! Add tendency by surface flux convergence ! xyz_DPTempDt(:,:,kz) = xyz_DPTempDt(:,:,kz) + xy_DPTempDtBulk do s = 1, ncmax xyzf_DQMixDt(:,:,kz,s) = xyzf_DQMixDt(:,:,kz,s) + xyf_DQMixDtBulk(:,:,s) end do pyz_DVelXDt (:,:,kz) = pyz_DVelXDt (:,:,kz) + py_DVelXDtBulk xqz_DVelYDt (:,:,kz) = xqz_DVelYDt (:,:,kz) + xq_DVelYDtBulk ! Output ! xyz_DPTempDtBulk = 0.0d0 xyzf_DQMixDtBulk = 0.0d0 pyz_DVelXDtBulk = 0.0d0 xqz_DVelYDtBulk = 0.0d0 xyz_DPTempDtBulk(:,:,kz) = xy_DPTempDtBulk do s = 1, ncmax xyzf_DQMixDtBulk(:,:,kz,s) = xyf_DQMixDtBulk(:,:,s) end do pyz_DVelXDtBulk (:,:,kz) = py_DVelXDtBulk xqz_DVelYDtBulk (:,:,kz) = xq_DVelYDtBulk ! call HistoryAutoPut(TimeN, 'PTempSfc', xyz_DPTempDtBulk(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelXSfc', pyz_DVelXDtBulk (1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelYSfc', xqz_DVelYDtBulk (1:nx,1:ny,1:nz)) do s = 1, ncmax call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_Sfc', & & xyzf_DQMixDtBulk(1:nx,1:ny,1:nz,s)) end do call HistoryAutoPut(TimeN, 'PTempSfcFlux', xy_PTempFlux(1:nx,1:ny)) call HistoryAutoPut(TimeN, 'VelXSfcFlux', py_VelXFlux (1:nx,1:ny)) call HistoryAutoPut(TimeN, 'VelYSfcFlux', xq_VelYFlux (1:nx,1:ny)) do s = 1, ncmax call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_SfcFlux', & & xyf_QMixFlux(1:nx,1:ny,s)) end do call HistoryAutoPut(TimeN, 'SfcHeatFlux', & & CpDry * xyz_DensBZ(1:nx,1:ny,1) * xy_PTempFlux(1:nx,1:ny) & & * ( ExnerBZSfc + xyz_Exner(1:nx,1:ny,1) ) ) call HistoryAutoPut(TimeN, 'SfcXMomFlux', & & xyz_DensBZ(1:nx,1:ny,1) * py_VelXFlux (1:nx,1:ny)) call HistoryAutoPut(TimeN, 'SfcYMomFlux', & & xyz_DensBZ(1:nx,1:ny,1) * xq_VelYFlux (1:nx,1:ny)) do s = 1, ncmax call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_SfcMassFlux', & & xyz_DensBZ(1:nx,1:ny,1) * xyf_QMixFlux(1:nx,1:ny,s)) end do ! Set Margin ! call SetMargin_xyz( xyz_DPTempDt ) call SetMargin_pyz( pyz_DVelXDt ) call SetMargin_xqz( xqz_DVelYDt ) call SetMargin_xyzf(xyzf_DQMixDt ) end subroutine Surfaceflux_Bulk_forcing end module Surfaceflux_bulk