scale_atmos_phy_sf_bulk.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "inc_openmp.h"
module scale_atmos_phy_sf_bulk
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_tracer
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_sf_bulk_setup
  public :: atmos_phy_sf_bulk

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

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_sf_bulk_setup( ATMOS_PHY_SF_TYPE )
    use scale_process, only: &
       prc_mpistop
    implicit none

    character(len=*), intent(in) :: atmos_phy_sf_type

    namelist / param_atmos_phy_sf_bulk / &
       atmos_phy_sf_beta

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

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[SURFACE FLUX] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
    if( io_l ) write(io_fid_log,*) '*** Bulk scheme'

    if ( atmos_phy_sf_type /= 'BULK' ) then
       write(*,*) 'xxx ATMOS_PHY_SF_TYPE is not BULK. Check!'
       call prc_mpistop
    endif

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_sf_bulk,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_SF_BULK. Check!'
       call prc_mpistop
    endif
    if( io_nml ) write(io_fid_nml,nml=param_atmos_phy_sf_bulk)

    return
  end subroutine atmos_phy_sf_bulk_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_sf_bulk( &
       ATM_TEMP, ATM_PRES, ATM_W, ATM_U, ATM_V,     &
       ATM_DENS,                                    &
       ATM_QTRC,                                    &
       ATM_Z1, dt,                                  &
       SFC_DENS, SFC_PRES,                          &
       SFLX_LW_dn, SFLX_SW_dn,                      &
       SFC_TEMP, SFC_albedo,                        &
       SFC_Z0M, SFC_Z0H, SFC_Z0E,                   &
       SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_LH, &
       SFLX_QTRC,                                   &
       U10, V10, T2, Q2                             )
    use scale_grid_index
    use scale_tracer
    use scale_const, only: &
       pre00 => const_pre00, &
       cpdry => const_cpdry, &
       rdry  => const_rdry
    use scale_atmos_hydrometeor, only: &
       hydrometeor_lhv => atmos_hydrometeor_lhv, &
       i_qv
    use scale_atmos_saturation, only: &
       saturation_pres2qsat_all => atmos_saturation_pres2qsat_all
    use scale_roughness, only: &
       roughness
    use scale_bulkflux, only: &
       bulkflux
    implicit none

    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_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_dens  (ia,ja)    ! density     at the lowermost layer (cell center) [kg/m3]
    real(RP), intent(in)    :: atm_qtrc  (ia,ja,qa) ! tracer      at the lowermost layer (cell center) [kg/kg]
    real(RP), intent(in)    :: atm_z1    (ia,ja)    ! height of the lowermost grid from surface (cell center) [m]
    real(DP), intent(in)    :: dt                   ! delta time
    real(RP), intent(in)    :: sfc_dens  (ia,ja)    ! density     at the surface atmosphere [kg/m3]
    real(RP), intent(in)    :: sfc_pres  (ia,ja)    ! pressure    at the surface atmosphere [Pa]
    real(RP), intent(in)    :: sflx_lw_dn(ia,ja)    ! downward longwave  radiation flux at the surface [J/m2/s]
    real(RP), intent(in)    :: sflx_sw_dn(ia,ja)    ! downward shortwave radiation flux at the surface [J/m2/s]
    real(RP), intent(in)    :: sfc_temp  (ia,ja)    ! temperature at the surface skin [K]
    real(RP), intent(in)    :: sfc_albedo(ia,ja,2)  ! surface albedo (LW/SW) (0-1)
    real(RP), intent(inout) :: sfc_z0m   (ia,ja)    ! surface roughness length (momentum) [m]
    real(RP), intent(inout) :: sfc_z0h   (ia,ja)    ! surface roughness length (heat) [m]
    real(RP), intent(inout) :: sfc_z0e   (ia,ja)    ! surface roughness length (vapor) [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_qtrc (ia,ja,qa) ! surface flux for tracer mass   (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) :: atm_qv   (ia,ja)
    real(RP) :: sfc_z0m_t(ia,ja)
    real(RP) :: sfc_z0h_t(ia,ja)
    real(RP) :: sfc_z0e_t(ia,ja)
    real(RP) :: sfc_qsat (ia,ja) ! saturatad water vapor mixing ratio [kg/kg]
    real(RP) :: sfc_qv   (ia,ja) ! water vapor mixing ratio [kg/kg]
    real(RP) :: lhv      (ia,ja)
    real(RP) :: pbl      (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) :: fracu10 ! calculation parameter for U10 [-]
    real(RP) :: fract2  ! calculation parameter for T2 [-]
    real(RP) :: fracq2  ! calculation parameter for Q2 [-]

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

    if( io_l ) write(io_fid_log,*) '*** Atmos physics  step: Surface flux(bulk)'

    call roughness( sfc_z0m_t(:,:), & ! [OUT]
                    sfc_z0h_t(:,:), & ! [OUT]
                    sfc_z0e_t(:,:), & ! [OUT]
                    sfc_z0m(:,:), & ! [IN]
                    sfc_z0h(:,:), & ! [IN]
                    sfc_z0e(:,:), & ! [IN]
                    atm_u(:,:), & ! [IN]
                    atm_v(:,:), & ! [IN]
                    atm_z1(:,:), & ! [IN]
                    dt              ) ! [IN]

    sfc_z0m(:,:) = sfc_z0m(:,:) + sfc_z0m_t(:,:) * dt
    sfc_z0h(:,:) = sfc_z0h(:,:) + sfc_z0h_t(:,:) * dt
    sfc_z0e(:,:) = sfc_z0e(:,:) + sfc_z0e_t(:,:) * dt

    call hydrometeor_lhv( lhv(:,:), atm_temp(:,:) )

    if ( i_qv > 0 ) then
       atm_qv(:,:) = atm_qtrc(:,:,i_qv)

       call saturation_pres2qsat_all( sfc_qsat(:,:), & ! [OUT]
                                      sfc_temp(:,:), & ! [IN]
                                      sfc_pres(:,:)  ) ! [IN]

       sfc_qv(:,:) = ( 1.0_rp - atmos_phy_sf_beta ) * atm_qv(:,:) + atmos_phy_sf_beta * sfc_qsat(:,:)
    else
       atm_qv(:,:) = 0.0_rp
       sfc_qv(:,:) = 0.0_rp
    end if

    sflx_qtrc(:,:,:) = 0.0_rp
    pbl(:,:)   = 100.0_rp ! tentative
    do j = js, je
    do i = is, ie
       call bulkflux( ustar,         & ! [OUT]
                      tstar,         & ! [OUT]
                      qstar,         & ! [OUT]
                      uabs,          & ! [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(i,j), & ! [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) = -atm_dens(i,j) * ustar**2 / uabs * atm_w(i,j)
       sflx_mu(i,j) = -atm_dens(i,j) * ustar**2 / uabs * atm_u(i,j)
       sflx_mv(i,j) = -atm_dens(i,j) * ustar**2 / uabs * atm_v(i,j)

       !-----< heat flux >-----
       sflx_sh(i,j) = - atm_dens(i,j) * ustar * tstar &
                    * cpdry * ( sfc_pres(i,j) / pre00 )**( rdry/cpdry )
       sflx_lh(i,j) = - atm_dens(i,j) * ustar * qstar * lhv(i,j)

       !-----< mass flux >-----
       if ( i_qv > 0 ) then
          sflx_qtrc(i,j,i_qv) = sflx_lh(i,j) / lhv(i,j)
       endif

       !-----< 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  (i,j) + 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(i,j) + ( atm_qv(i,j) - sfc_qv(i,j) ) &
                                * ( 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

end module scale_atmos_phy_sf_bulk