scale_atmos_phy_sf_bulk.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_atmos_phy_sf_bulk
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_sf_bulk_setup
  public :: atmos_phy_sf_bulk_flux

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  real(RP), private :: atmos_phy_sf_bulk_beta = 1.0_rp ! evaporation efficiency (0-1)

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_sf_bulk_setup
    use scale_prc, only: &
       prc_abort
    implicit none

    namelist / param_atmos_phy_sf_bulk / &
       atmos_phy_sf_bulk_beta

    integer :: ierr
    !---------------------------------------------------------------------------

    log_newline
    log_info("ATMOS_PHY_SF_bulk_setup",*) 'Setup'
    log_info("ATMOS_PHY_SF_bulk_setup",*) 'Bulk scheme'

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_sf_bulk,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("ATMOS_PHY_SF_bulk_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_PHY_SF_bulk_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_SF_BULK. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_phy_sf_bulk)

    return
  end subroutine atmos_phy_sf_bulk_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_sf_bulk_flux( &
       IA, IS, IE, JA, JS, JE, &
       ATM_W, ATM_U, ATM_V,          &
       ATM_TEMP, ATM_PRES, ATM_QV,   &
       SFC_DENS, SFC_TEMP, SFC_PRES, &
       SFC_Z0M, SFC_Z0H, SFC_Z0E,    &
       PBL, ATM_Z1,                  &
       SFLX_MW, SFLX_MU, SFLX_MV,    &
       SFLX_SH, SFLX_LH, SFLX_QV,    &
       U10, V10, T2, Q2              )
    use scale_const, only: &
       cpdry  => const_cpdry, &
       epsvap => const_epsvap
    use scale_atmos_hydrometeor, only: &
       hydrometeor_lhv => atmos_hydrometeor_lhv
    use scale_atmos_saturation, only: &
       saturation_psat_all => atmos_saturation_psat_all
    use scale_bulkflux, only: &
       bulkflux
    implicit none
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE

    real(RP), intent(in) :: ATM_W   (ia,ja) ! velocity w  at the lowermost layer (cell center) [m/s]
    real(RP), intent(in) :: ATM_U   (ia,ja) ! velocity u  at the lowermost layer (cell center) [m/s]
    real(RP), intent(in) :: ATM_V   (ia,ja) ! velocity v  at the lowermost layer (cell center) [m/s]
    real(RP), intent(in) :: ATM_TEMP(ia,ja) ! temperature at the lowermost layer (cell center) [K]
    real(RP), intent(in) :: ATM_PRES(ia,ja) ! pressure    at the lowermost layer (cell center) [Pa]
    real(RP), intent(in) :: ATM_QV  (ia,ja) ! qv          at the lowermost layer (cell center) [kg/kg]
    real(RP), intent(in) :: SFC_DENS(ia,ja) ! density     at the surface atmosphere [kg/m3]
    real(RP), intent(in) :: SFC_TEMP(ia,ja) ! temperature at the surface skin [K]
    real(RP), intent(in) :: SFC_PRES(ia,ja) ! pressure    at the surface atmosphere [Pa]
    real(RP), intent(in) :: SFC_Z0M (ia,ja) ! surface roughness length (momentum) [m]
    real(RP), intent(in) :: SFC_Z0H (ia,ja) ! surface roughness length (heat) [m]
    real(RP), intent(in) :: SFC_Z0E (ia,ja) ! surface roughness length (vapor) [m]
    real(RP), intent(in) :: PBL     (ia,ja) ! depth of the PBL [m]
    real(RP), intent(in) :: ATM_Z1  (ia,ja) ! height of the lowermost grid from surface (cell center) [m]

    real(RP), intent(out) :: SFLX_MW(ia,ja) ! surface flux for z-momentum    (area center)   [m/s*kg/m2/s]
    real(RP), intent(out) :: SFLX_MU(ia,ja) ! surface flux for x-momentum    (area center)   [m/s*kg/m2/s]
    real(RP), intent(out) :: SFLX_MV(ia,ja) ! surface flux for y-momentum    (area center)   [m/s*kg/m2/s]
    real(RP), intent(out) :: SFLX_SH(ia,ja) ! surface flux for sensible heat (area center)   [J/m2/s]
    real(RP), intent(out) :: SFLX_LH(ia,ja) ! surface flux for latent   heat (area center)   [J/m2/s]
    real(RP), intent(out) :: SFLX_QV(ia,ja) ! surface flux for qv            (area center)   [kg/m2/s]
    real(RP), intent(out) :: U10    (ia,ja) ! velocity u        at 10m height
    real(RP), intent(out) :: V10    (ia,ja) ! velocity v        at 10m height
    real(RP), intent(out) :: T2     (ia,ja) ! temperature t     at  2m height
    real(RP), intent(out) :: Q2     (ia,ja) ! water vapor q     at  2m height

    real(RP) :: SFC_PSAT (ia,ja) ! saturatad water vapor pressure [Pa]
    real(RP) :: LHV      (ia,ja)

    real(RP) :: Ustar    ! friction velocity [m]
    real(RP) :: Tstar    ! friction temperature [K]
    real(RP) :: Qstar    ! friction mixing rate [kg/kg]
    real(RP) :: Uabs     ! modified absolute velocity [m/s]
    real(RP) :: Ra       ! Aerodynamic resistance (=1/Ce) [1/s]
    real(RP) :: SFC_QSAT ! saturatad water vapor mixing ratio [kg/kg]
    real(RP) :: SFC_QV   ! water vapor mixing ratio [kg/kg]

    real(RP) :: FracU10 ! calculation parameter for U10 [-]
    real(RP) :: FracT2  ! calculation parameter for T2 [-]
    real(RP) :: FracQ2  ! calculation parameter for Q2 [-]

    integer  :: i, j
    !---------------------------------------------------------------------------

    log_progress(*) 'atmosphere / physics / surface flux / bulk'

    ! ToDo consider ATM_TEMP is appropriate
    call hydrometeor_lhv( ia, is, ie, ja, js, je, &
                          atm_temp(:,:), & ! [IN]
                          lhv(:,:)       ) ! [OUT]

    call saturation_psat_all( ia, is, ie, ja, js, je, &
                              sfc_temp(:,:), & ! [IN]
                              sfc_psat(:,:)  ) ! [OUT]

#ifndef __GFORTRAN__
    !$omp parallel do default(none) &
    !$omp private(SFC_QSAT,SFC_QV,Ustar,Tstar,Qstar,Uabs,Ra,FracU10,FracT2,FracQ2) &
    !$omp shared (IS,IE,JS,JE,EPSvap,ATMOS_PHY_SF_BULK_beta,CPdry,LHV,bulkflux,    &
    !$omp         ATM_TEMP,ATM_PRES,ATM_QV,ATM_W,ATM_U,ATM_V,ATM_Z1,               &
    !$omp         SFC_DENS,SFC_TEMP,SFC_PRES,SFC_PSAT,SFC_Z0M,SFC_Z0H,SFC_Z0E,PBL, &
    !$omp         SFLX_MW,SFLX_MU,SFLX_MV,SFLX_SH,SFLX_LH,SFLX_QV,U10,V10,T2,Q2)
#else
    !$omp parallel do default(shared) &
    !$omp private(SFC_QSAT,SFC_QV,Ustar,Tstar,Qstar,Uabs,Ra,FracU10,FracT2,FracQ2)
#endif
    do j = js, je
    do i = is, ie
       ! qdry = 1 - psat
       sfc_qsat = epsvap * sfc_psat(i,j) / ( sfc_pres(i,j) - ( 1.0_rp-epsvap ) * sfc_psat(i,j) )

       sfc_qv = ( 1.0_rp - atmos_phy_sf_bulk_beta ) * atm_qv(i,j) + atmos_phy_sf_bulk_beta * sfc_qsat

       call bulkflux( ustar,         & ! [OUT]
                      tstar,         & ! [OUT]
                      qstar,         & ! [OUT]
                      uabs,          & ! [OUT]
                      ra,            & ! [OUT]
                      fracu10,       & ! [OUT]
                      fract2,        & ! [OUT]
                      fracq2,        & ! [OUT]
                      atm_temp(i,j), & ! [IN]
                      sfc_temp(i,j), & ! [IN]
                      atm_pres(i,j), & ! [IN]
                      sfc_pres(i,j), & ! [IN]
                      atm_qv(i,j), & ! [IN]
                      sfc_qv       , & ! [IN]
                      atm_u(i,j), & ! [IN]
                      atm_v(i,j), & ! [IN]
                      atm_z1(i,j), & ! [IN]
                      pbl(i,j), & ! [IN]
                      sfc_z0m(i,j), & ! [IN]
                      sfc_z0h(i,j), & ! [IN]
                      sfc_z0e(i,j)  ) ! [IN]

       !-----< momentum >-----
       sflx_mw(i,j) = -sfc_dens(i,j) * ustar * ustar / uabs * atm_w(i,j)
       sflx_mu(i,j) = -sfc_dens(i,j) * ustar * ustar / uabs * atm_u(i,j)
       sflx_mv(i,j) = -sfc_dens(i,j) * ustar * ustar / uabs * atm_v(i,j)

       !-----< heat flux >-----
       sflx_sh(i,j) = -sfc_dens(i,j) * ustar * tstar * cpdry
       sflx_lh(i,j) = -sfc_dens(i,j) * ustar * qstar * lhv(i,j)

       !-----< mass flux >-----
       sflx_qv(i,j) = sflx_lh(i,j) / lhv(i,j)

       !-----< U10, T2, q2 >-----
       !U10(i,j) = FracU10 * ATM_U(i,j)
       !V10(i,j) = FracU10 * ATM_V(i,j)
       !T2 (i,j) = ( 1.0_RP - FracT2 ) * SFC_TEMP(i,j) + FracT2 * ATM_TEMP(i,j)
       !Q2 (i,j) = ( 1.0_RP - FracQ2 ) * SFC_QV        + FracQ2 * ATM_QV  (i,j)

       u10(i,j) = atm_u(i,j) * log( 10.0_rp / sfc_z0m(i,j) ) / log( atm_z1(i,j) / sfc_z0m(i,j) )
       v10(i,j) = atm_v(i,j) * log( 10.0_rp / sfc_z0m(i,j) ) / log( atm_z1(i,j) / sfc_z0m(i,j) )
       t2(i,j) = sfc_temp(i,j) + ( atm_temp(i,j) - sfc_temp(i,j) ) &
                                * ( log(      2.0_rp / sfc_z0m(i,j) ) * log(      2.0_rp / sfc_z0h(i,j) ) ) &
                                / ( log( atm_z1(i,j) / sfc_z0m(i,j) ) * log( atm_z1(i,j) / sfc_z0h(i,j) ) )
       q2(i,j) = sfc_qv        + ( atm_qv(i,j) - sfc_qv        ) &
                                * ( log(      2.0_rp / sfc_z0m(i,j) ) * log(      2.0_rp / sfc_z0e(i,j) ) ) &
                                / ( log( atm_z1(i,j) / sfc_z0m(i,j) ) * log( atm_z1(i,j) / sfc_z0e(i,j) ) )
    enddo
    enddo

    return
  end subroutine atmos_phy_sf_bulk_flux

end module scale_atmos_phy_sf_bulk