scale_atmos_phy_sf_const.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "inc_openmp.h"
module scale_atmos_phy_sf_const
  !-----------------------------------------------------------------------------
  !
  !++ 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_const_setup
  public :: atmos_phy_sf_const

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer,  private            :: ATMOS_PHY_SF_FLG_MOM_FLUX = 0 ! application type for momentum flux
                                                                ! 0: Bulk coefficient  is constant
                                                                ! 1: Friction velocity is constant

  real(RP), private, parameter :: ATMOS_PHY_SF_U_maxM      =  100.0_rp ! maximum limit of absolute velocity for momentum [m/s]
  real(RP), private            :: ATMOS_PHY_SF_U_minM      =    0.0_rp ! minimum limit of absolute velocity for momentum [m/s]
  real(RP), private, parameter :: ATMOS_PHY_SF_Cm_max      = 2.5e-3_rp ! maximum limit of bulk coefficient for momentum [NIL]
  real(RP), private            :: ATMOS_PHY_SF_Cm_min      = 1.0e-5_rp ! minimum limit of bulk coefficient for momentum [NIL]

  real(RP), private            :: ATMOS_PHY_SF_Const_Ustar =   0.25_rp ! constant friction velocity [m/s]
  real(RP), private            :: ATMOS_PHY_SF_Const_Cm    = 0.0011_rp ! constant bulk coefficient for momentum [NIL]
  real(RP), private            :: ATMOS_PHY_SF_Const_SH    =   15.0_rp ! constant surface sensible heat flux [W/m2]
  real(RP), private            :: ATMOS_PHY_SF_Const_LH    =  115.0_rp ! constant surface latent   heat flux [W/m2]

  logical,  private            :: ATMOS_PHY_SF_FLG_SH_DIURNAL = .false. ! diurnal modulation for sensible heat flux?
  real(RP), private            :: ATMOS_PHY_SF_Const_FREQ     = 24.0_rp ! frequency of sensible heat flux modulation [hour]

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_sf_const_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_const / &
       atmos_phy_sf_flg_mom_flux,   &
       atmos_phy_sf_u_minm,         &
       atmos_phy_sf_cm_min,         &
       atmos_phy_sf_const_ustar,    &
       atmos_phy_sf_const_cm,       &
       atmos_phy_sf_const_sh,       &
       atmos_phy_sf_const_lh,       &
       atmos_phy_sf_flg_sh_diurnal, &
       atmos_phy_sf_const_freq

    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,*) '*** Constant flux'

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

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_sf_const,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_CONST. Check!'
       call prc_mpistop
    endif
    if( io_nml ) write(io_fid_nml,nml=param_atmos_phy_sf_const)

    return
  end subroutine atmos_phy_sf_const_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_sf_const( &
       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: &
       pi    => const_pi
    use scale_atmos_hydrometeor, only: &
       hydrometeor_lhv => atmos_hydrometeor_lhv, &
       i_qv
    use scale_time, only: &
       time_nowsec
    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_uabs(ia,ja) ! absolute velocity at z1 [m/s]

    real(RP) :: cm(ia,ja)       ! bulk coefficient for momentum
    real(RP) :: r10

    real(RP) :: modulation
    real(RP) :: uabs_lim
    real(RP) :: lhv(ia,ja)
    integer  :: i, j
    !---------------------------------------------------------------------------

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

    do j = js, je
    do i = is, ie
       atm_uabs(i,j) = sqrt( atm_w(i,j)*atm_w(i,j) &
                           + atm_u(i,j)*atm_u(i,j) &
                           + atm_v(i,j)*atm_v(i,j) ) ! at cell center
    enddo
    enddo

    if   ( atmos_phy_sf_flg_mom_flux == 0 ) then ! Bulk coefficient is constant
       do j = js, je
       do i = is, ie
          cm(i,j) = atmos_phy_sf_const_cm
       enddo
       enddo
    elseif( atmos_phy_sf_flg_mom_flux == 1 ) then ! Friction velocity is constant
       do j = js, je
       do i = is, ie
          cm(i,j) = ( atmos_phy_sf_const_ustar / atm_uabs(i,j) )**2
          cm(i,j) = min( max( cm(i,j), atmos_phy_sf_cm_min ), atmos_phy_sf_cm_max )
       enddo
       enddo
    endif

    !-----< momentum >-----

    do j = js, je
    do i = is, ie
       uabs_lim = min( max( atm_uabs(i,j), atmos_phy_sf_u_minm ), atmos_phy_sf_u_maxm )

       sflx_mw(i,j) = -cm(i,j) * uabs_lim * sfc_dens(i,j) * atm_w(i,j)
       sflx_mu(i,j) = -cm(i,j) * uabs_lim * sfc_dens(i,j) * atm_u(i,j)
       sflx_mv(i,j) = -cm(i,j) * uabs_lim * sfc_dens(i,j) * atm_v(i,j)
    enddo
    enddo

    !-----< heat flux >-----

    if ( atmos_phy_sf_flg_sh_diurnal ) then
       modulation = sin( 2.0_rp * pi * time_nowsec / 3600.0_rp / atmos_phy_sf_const_freq )
    else
       modulation = 1.0_rp
    endif

    do j = js, je
    do i = is, ie
       sflx_sh(i,j) = atmos_phy_sf_const_sh * modulation
       sflx_lh(i,j) = atmos_phy_sf_const_lh
    enddo
    enddo

    !-----< mass flux >-----
   call hydrometeor_lhv( lhv(:,:), atm_temp(:,:) )

    sflx_qtrc(:,:,:) = 0.0_rp
    if ( i_qv > 0 ) then
       do j = js, je
       do i = is, ie
          sflx_qtrc(i,j,i_qv) = sflx_lh(i,j) / lhv(i,j)
       enddo
       enddo
    end if

    !-----< U10, T2, q2 >-----

    do j = js, je
    do i = is, ie
       r10 = 10.0_rp / atm_z1(i,j)

       u10(i,j) = r10 * atm_u(i,j)
       v10(i,j) = r10 * atm_v(i,j)
    enddo
    enddo

    do j = js, je
    do i = is, ie
       t2(i,j) = atm_temp(i,j)
    enddo
    enddo

    if ( i_qv > 0 ) then
       do j = js, je
       do i = is, ie
          q2(i,j) = atm_qtrc(i,j,i_qv)
       enddo
       enddo
    else
       do j = js, je
       do i = is, ie
          q2(i,j) = 0.0_rp
       enddo
       enddo
    end if

    return
  end subroutine atmos_phy_sf_const

end module scale_atmos_phy_sf_const