--- surfaceflux_bulk.bk.f90 2011-06-18 04:07:27.000000000 +0900 +++ surfaceflux_bulk.f90 2011-09-08 15:59:18.000000000 +0900 @@ -1,23 +1,27 @@ != Module HeatFlux ! -! Authors:: ODAKA Masatsugu -! Version:: $Id: surfaceflux_bulk.f90,v 1.8 2011-06-17 19:07:27 sugiyama Exp $ -! Tag Name:: $Name: arare5-20110623 $ +! Authors:: ODAKA Masatsugu, TAKAHASHI Yoshiyuki +! Version:: $Id: surfaceflux_bulk.f90,v 1.10 2011-09-04 14:40:00 yot Exp $ +! Tag Name:: $Name: $ ! 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 +! H = Fh/Δz_1 +! Q = Fq/Δz_1 ! -! と表される. ここで Δz_1 は最下層の格子間隔である. +! と表される. ここで Δz_1 は最下層の格子間隔である. ! ! 熱フラックス Fh と凝結物質のフラックス Fq は以下の式にしたがって ! 計算する. @@ -51,7 +55,7 @@ ! !下部境界でのフラックスの計算モジュール ! - + !モジュール読み込み use dc_types, only: DP, STRING use dc_iounit, only: FileOpen @@ -68,13 +72,15 @@ & nx, ny, nz, ncmax use axesset, only: z_dz, & !z 方向の格子点間隔 & xyz_avr_pyz, & - & xyz_avr_xqz + & xyz_avr_xqz, & + & pyz_avr_xyz, & + & xqz_avr_xyz use basicset, only: xyz_ExnerBZ, & !エクスナー関数の基本場 & xyz_PressBZ, & ! & xyz_PTempBZ, & !温位の基本場 & xyz_TempBZ, & ! & xyzf_QMixBZ !温位の基本場 - use constants,only: MolWtDry + 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 @@ -93,10 +99,19 @@ !変数定義 real(DP), save :: Bulk = 1.5d-3 !熱・運動量フラックスのバルク係数 + + + ! Variables, VelX0, VelY0, should be unified to one variable, such as Vel0, + ! since only the magnitude of the value is meaningful. (YOT, 2011/08/31) + ! real(DP), save :: VelX0 = 0.0d0 !下層での水平速度嵩上げ値 real(DP), save :: VelY0 = 0.0d0 !下層での水平速度嵩上げ値 + character(*), parameter:: module_name = 'surfaceflux_bulk' + ! モジュールの名称. + ! Module name + contains !!!------------------------------------------------------------------------!!! subroutine Surfaceflux_Bulk_init @@ -110,50 +125,82 @@ !内部変数 integer :: l, unit - !--------------------------------------------------------------- + !--------------------------------------------------------------- ! NAMELIST から情報を取得 ! NAMELIST /surfaceflux_bulk_nml/ Bulk, VelX0, VelY0 - + call FileOpen(unit, file=namelist_filename, mode='r') read(unit, NML=surfaceflux_bulk_nml) close(unit) if (myrank == 0) then - call MessageNotify( "M", "SurfaceSfc", "Bulk = %f", d=(/Bulk/)) - call MessageNotify( "M", "SurfaceSfc", "VelX0 = %f", d=(/VelX0/)) - call MessageNotify( "M", "SurfaceSfc", "VelY0 = %f", d=(/VelY0/)) + call MessageNotify( "M", module_name, "Bulk = %f", d=(/Bulk/) ) + call MessageNotify( "M", module_name, "VelX0 = %f", d=(/VelX0/)) + call MessageNotify( "M", module_name, "VelY0 = %f", d=(/VelY0/)) 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='surface flux of potential temperature', & - & units='kg.kg-1.s-1', & - & xtype='double') + & longname='potential temperature tendency by surface flux', & + & units='K.s-1', & + & xtype='float') call HistoryAutoAddVariable( & & varname='VelXSfc', & & dims=(/'x','y','z','t'/), & - & longname='flux of heating from surface', & - & units='kg.m.s-1', & - & xtype='double') + & longname='x-component velocity tendency by surface flux', & + & units='m.s-2', & + & xtype='float') call HistoryAutoAddVariable( & & varname='VelYSfc', & & dims=(/'x','y','z','t'/), & - & longname='flux of momentum from surface', & - & units='kg.m.s-1', & - & xtype='double') + & 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='Surface Flux term of ' & - & //trim(SpcWetSymbol(l))//' mixing ratio', & - & units='kg.kg-1.s-1', & - & xtype='double') + & longname=trim(SpcWetSymbol(l))//' mixing ratio tendency by surface flux', & + & units='s-1', & + & xtype='float') end do end subroutine Surfaceflux_Bulk_init @@ -171,7 +218,7 @@ !暗黙の型宣言禁止 implicit none - + !変数定義 real(DP), intent(in) :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax) !水平風速 @@ -187,114 +234,175 @@ 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) :: pyz_DVelXDt0(imin:imax,jmin:jmax,kmin:kmax) - real(DP) :: xqz_DVelYDt0(imin:imax,jmin:jmax,kmin:kmax) - real(DP) :: xyz_DPTempDt0(imin:imax,jmin:jmax,kmin:kmax) - real(DP) :: xyzf_DQMixDt0(imin:imax,jmin:jmax,kmin:kmax, ncmax) - real(DP) :: pyz_Momflux(imin:imax,jmin:jmax,kmin:kmax) + real(DP) :: py_VelXflux (imin:imax,jmin:jmax) !運動量フラックス - real(DP) :: xqz_Momflux(imin:imax,jmin:jmax,kmin:kmax) + !(strictly speaking, this value is not + !momentum flux, but is it divided by density) + real(DP) :: xq_VelYflux (imin:imax,jmin:jmax) !運動量フラックス - real(DP) :: xyz_Heatflux(imin:imax,jmin:jmax,kmin:kmax) + !(strictly speaking, this value is not + !momentum flux, but is it divided by density) + real(DP) :: xy_PTempFlux(imin:imax,jmin:jmax) !地表面熱フラックス - real(DP) :: xyzf_QMixflux(imin:imax,jmin:jmax,kmin:kmax, ncmax) + !(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) :: pyz_VelXSfc(imin:imax,jmin:jmax,kmin:kmax) - real(DP) :: xqz_VelYSfc(imin:imax,jmin:jmax,kmin:kmax) - real(DP) :: xyz_VelXSfc(imin:imax,jmin:jmax,kmin:kmax) + real(DP) :: xyz_VelX(imin:imax,jmin:jmax,kmin:kmax) !水平風速 (xyz 格子) - real(DP) :: xyz_VelYSfc(imin:imax,jmin:jmax,kmin:kmax) + real(DP) :: xyz_VelY(imin:imax,jmin:jmax,kmin:kmax) !水平風速 (xyz 格子) - real(DP) :: TempSfc - real(DP) :: PressSfc - real(DP) :: xyz_VelSfc(imin:imax,jmin:jmax,kmin:kmax) + real(DP) :: xyz_AbsVel(imin:imax,jmin:jmax,kmin:kmax) !水平風速 (xyz 格子) - real(DP) :: xyz_TempAll(imin:imax,jmin:jmax,kmin:kmax) + 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_Heatflux = 0.0d0 - xyzf_QMixflux = 0.0d0 - pyz_Momflux = 0.0d0 - xqz_Momflux = 0.0d0 - - pyz_DVelXDt0 = pyz_DVelXDt - xqz_DVelYDt0 = xqz_DVelYDt - xyz_DPTempDt0 = xyz_DPTempDt - xyzf_DQMixDt0 = xyzf_DQMixDt - - xyz_TempAll = (xyz_Exner + xyz_ExnerBZ) * (xyz_PTemp + xyz_PTempBZ) - xyzf_QMixAll = xyzf_QMix + xyzf_QMixBZ - - pyz_VelXSfc = pyz_VelX + VelX0 - xqz_VelYSfc = xqz_VelY + VelY0 - ! 地表面での値 - ! - TempSfc = xyz_TempBZ(1,1,kz) !水平一様なので, i=j=1 で代表させる - PressSfc= xyz_PressBZ(1,1,kz) !水平一様なので, i=j=1 で代表させる - - !地表面熱フラックスによる加熱率を計算 - ! * 単位は K/s - ! * エクスナー関数は基本場の値で代表させる. - ! * 格子点 xz では, 物理領域の最下端の添え字は kz = 1 - xyz_VelXSfc = xyz_avr_pyz(pyz_VelXSfc) - xyz_VelYSfc = xyz_avr_xqz(xqz_VelYSfc) - - xyz_VelSfc = SQRT( xyz_VelXSfc ** 2.0 + xyz_VelYSfc ** 2.0 ) - - xyz_Heatflux(:,:,kz) = & - & MAX( & - & 0.0d0, & - & - Bulk * xyz_VelSfc(:,:,kz) & - & * ( xyz_TempAll(:,:,kz) - TempSfc ) & - & / ( z_dz(kz) * 5.0d-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) + + ! Magnitude of wind velocities at xyz, pyz, and xqz grid points are calculated. + ! Variables, VelX0, VelY0, should be unified to one variable, such as Vel0, + ! since only the magnitude of the value is meaningful. (YOT, 2011/08/31) + ! + xyz_AbsVel = SQRT( & + & ( abs( xyz_VelX ) + VelX0 )** 2.0 & + & + ( abs( xyz_VelY ) + VelY0 )** 2.0 & + & ) + 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 - xyzf_QMixflux(:,:,kz,IdxCG(s)) = & - & MAX( & - & 0.0d0, & - & - Bulk * xyz_VelSfc(:,:,kz) & - & * ( xyzf_QMixAll(:,:,kz,s) & - & - SvapPress( SpcWetID(IdxCC(s)), TempSfc ) & - & / PressSfc * (MolWtWet(IdxCG(s)) / MolWtDry)& - & ) / ( z_dz(kz) * 5.0d-1 ) & - & ) + 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) - !地表面運動量フラックスによる変化率を計算 - ! * 単位は m/s^2 - ! * 格子点 xz では, 物理領域の最下端の添え字は kz = 1 - pyz_MomFlux(:,:,kz) = & - & MAX( & - & 0.0d0, & - & - Bulk * abs(pyz_VelXSfc(:,:,kz)) * pyz_VelX(:,:,kz) & - & / ( z_dz(kz) * 5.0d-1 ) & - & ) - - xqz_MomFlux(:,:,kz) = & - & MAX( & - & 0.0d0, & - & - Bulk * abs(xyz_VelYSfc(:,:,kz)) * xqz_VelY(:,:,kz) & - & / ( z_dz(kz) * 5.0d-1 ) & - & ) - - xyz_DPTempDt = xyz_DPTempDt0 + xyz_Heatflux - xyzf_DQMixDt = xyzf_DQMixDt0 + xyzf_Qmixflux - pyz_DVelXDt = pyz_DVelXDt0 + pyz_MomFlux - xqz_DVelYDt = xqz_DVelYDt0 + xqz_MomFlux - - call HistoryAutoPut(TimeN, 'PTempSfc', xyz_HeatFlux(1:nx,1:ny,1:nz)) - call HistoryAutoPut(TimeN, 'VelXSfc', pyz_MomFlux(1:nx,1:ny,1:nz)) - call HistoryAutoPut(TimeN, 'VelYSfc', xqz_MomFlux(1:nx,1:ny,1:nz)) + + ! Add tendency by surface flux convergence + ! + xyz_DPTempDt(:,:,kz) = xyz_DPTempDt(:,:,kz) + xy_DPTempDtBulk do s = 1, ncmax - call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_Sfc', xyzf_Qmixflux(1:nx,1:ny,1:nz,s)) + 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 + end subroutine Surfaceflux_Bulk_forcing end module Surfaceflux_bulk