mod_atmos_phy_sf_driver.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module mod_atmos_phy_sf_driver
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_atmos_grid_cartesc_index
  use scale_tracer
  use scale_cpl_sfc_index
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_sf_driver_setup
  public :: atmos_phy_sf_driver_calc_tendency

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_sf_driver_setup
    use scale_prc, only: &
       prc_abort
    use scale_atmos_phy_sf_bulk, only: &
       atmos_phy_sf_bulk_setup
    use scale_atmos_phy_sf_const, only: &
       atmos_phy_sf_const_setup
    use mod_atmos_admin, only: &
       atmos_phy_sf_type, &
       atmos_sw_phy_sf
    use mod_atmos_phy_sf_vars, only: &
       sfc_z0m   => atmos_phy_sf_sfc_z0m,   &
       sfc_z0h   => atmos_phy_sf_sfc_z0h,   &
       sfc_z0e   => atmos_phy_sf_sfc_z0e,   &
       sflx_mw   => atmos_phy_sf_sflx_mw,   &
       sflx_mu   => atmos_phy_sf_sflx_mu,   &
       sflx_mv   => atmos_phy_sf_sflx_mv,   &
       sflx_sh   => atmos_phy_sf_sflx_sh,   &
       sflx_lh   => atmos_phy_sf_sflx_lh,   &
       sflx_qtrc => atmos_phy_sf_sflx_qtrc
    use mod_cpl_admin, only: &
       cpl_sw
    implicit none
    !---------------------------------------------------------------------------

    log_newline
    log_info("ATMOS_PHY_SF_driver_setup",*) 'Setup'

    if ( atmos_sw_phy_sf ) then

       if ( cpl_sw ) then
          log_info("ATMOS_PHY_SF_driver_setup",*) 'Coupler is enabled.'
       else
          ! setup library component
          select case( atmos_phy_sf_type )
          case ( 'BULK' )
             call atmos_phy_sf_bulk_setup
          case ( 'CONST' )
             call atmos_phy_sf_const_setup
          case default
             log_error("ATMOS_PHY_SF_driver_setup",*) 'invalid Surface flux type(', trim(atmos_phy_sf_type), '). CHECK!'
             call prc_abort
          end select
       endif

    else

       log_info("ATMOS_PHY_SF_driver_setup",*) 'this component is never called.'
       log_info("ATMOS_PHY_SF_driver_setup",*) 'surface fluxes are set to zero.'
       sflx_mw(:,:) = 0.0_rp
       sflx_mu(:,:) = 0.0_rp
       sflx_mv(:,:) = 0.0_rp
       sflx_sh(:,:) = 0.0_rp
       sflx_lh(:,:) = 0.0_rp
       log_info("ATMOS_PHY_SF_driver_setup",*) 'SFC_TEMP, SFC_albedo is set in ATMOS_PHY_SF_vars.'

    endif

    sflx_qtrc(:,:,:) = 0.0_rp

    return
  end subroutine atmos_phy_sf_driver_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_sf_driver_calc_tendency( update_flag )
    use scale_const, only: &
       eps    => const_eps,    &
       grav   => const_grav,   &
       karman => const_karman, &
       cpdry  => const_cpdry
    use scale_atmos_grid_cartesc_real, only: &
       cz => atmos_grid_cartesc_real_cz, &
       fz => atmos_grid_cartesc_real_fz, &
       z1 => atmos_grid_cartesc_real_z1, &
       atmos_grid_cartesc_real_area, &
       atmos_grid_cartesc_real_totarea
    use scale_topography, only: &
       topo_zsfc
    use scale_time, only: &
       dt_sf => time_dtsec_atmos_phy_sf
    use scale_statistics, only: &
       statistics_checktotal, &
       statistics_total
    use scale_atmos_bottom, only: &
       bottom_estimate => atmos_bottom_estimate
    use scale_atmos_hydrometeor, only: &
       atmos_hydrometeor_dry, &
       i_qv
    use scale_atmos_phy_sf_bulk, only: &
       atmos_phy_sf_bulk_flux
    use scale_atmos_phy_sf_const, only: &
       atmos_phy_sf_const_flux
    use mod_atmos_admin, only: &
       atmos_phy_sf_type
    use mod_atmos_vars, only: &
       dens   => dens_av, &
       rhot   => rhot_av, &
       pott,              &
       temp,              &
       pres,              &
       w,                 &
       u,                 &
       v,                 &
       qv,                &
       dens_t => dens_tp, &
       momz_t => momz_tp, &
       rhou_t => rhou_tp, &
       rhov_t => rhov_tp, &
       rhoh   => rhoh_p,  &
       rhot_t => rhot_tp, &
       rhoq_t => rhoq_tp
    use mod_atmos_phy_rd_vars, only: &
       sflx_lw_dn => atmos_phy_rd_sflx_lw_dn, &
       sflx_sw_dn => atmos_phy_rd_sflx_sw_dn
    use mod_atmos_phy_bl_vars, only: &
       pbl_zi => atmos_phy_bl_zi
    use mod_atmos_phy_sf_vars, only: &
       dens_t_sf => atmos_phy_sf_dens_t,    &
       momz_t_sf => atmos_phy_sf_momz_t,    &
       rhou_t_sf => atmos_phy_sf_rhou_t,    &
       rhov_t_sf => atmos_phy_sf_rhov_t,    &
       rhoh_sf   => atmos_phy_sf_rhoh,      &
       rhot_t_sf => atmos_phy_sf_rhot_t,    &
       rhoq_t_sf => atmos_phy_sf_rhoq_t,    &
       sfc_dens  => atmos_phy_sf_sfc_dens,  &
       sfc_pres  => atmos_phy_sf_sfc_pres,  &
       sfc_temp  => atmos_phy_sf_sfc_temp,  &
       sfc_z0m   => atmos_phy_sf_sfc_z0m,   &
       sfc_z0h   => atmos_phy_sf_sfc_z0h,   &
       sfc_z0e   => atmos_phy_sf_sfc_z0e,   &
       sflx_mw   => atmos_phy_sf_sflx_mw,   &
       sflx_mu   => atmos_phy_sf_sflx_mu,   &
       sflx_mv   => atmos_phy_sf_sflx_mv,   &
       sflx_sh   => atmos_phy_sf_sflx_sh,   &
       sflx_lh   => atmos_phy_sf_sflx_lh,   &
       sflx_gh   => atmos_phy_sf_sflx_gh,   &
       sflx_qtrc => atmos_phy_sf_sflx_qtrc, &
       u10       => atmos_phy_sf_u10,       &
       v10       => atmos_phy_sf_v10,       &
       t2        => atmos_phy_sf_t2,        &
       q2        => atmos_phy_sf_q2,        &
       l_mo      => atmos_phy_sf_l_mo
    use mod_cpl_admin, only: &
       cpl_sw
    implicit none

    logical, intent(in) :: update_flag

    real(RP) :: ATM_W   (ia,ja)
    real(RP) :: ATM_U   (ia,ja)
    real(RP) :: ATM_V   (ia,ja)
    real(RP) :: ATM_DENS(ia,ja)
    real(RP) :: ATM_TEMP(ia,ja)
    real(RP) :: ATM_PRES(ia,ja)
    real(RP) :: ATM_QV  (ia,ja)
    real(RP) :: SFLX_QV (ia,ja)
    real(RP) :: us, SFLX_PT
    real(RP) :: work

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

    if ( update_flag ) then

       ! update surface density, surface pressure
       call bottom_estimate( ka, ks,  ke,                        & ! [IN]
                             ia, isb, ieb,                       & ! [IN]
                             ja, jsb, jeb,                       & ! [IN]
                             dens(:,:,:), pres(:,:,:),           & ! [IN]
                             cz(:,:,:), topo_zsfc(:,:), z1(:,:), & ! [IN]
                             sfc_dens(:,:), sfc_pres(:,:)        ) ! [OUT]

       if ( .NOT. cpl_sw ) then

          !$omp parallel do
          do j = jsb, jeb
          do i = isb, ieb
             atm_w(i,j) = w(ks,i,j)
             atm_u(i,j) = u(ks,i,j)
             atm_v(i,j) = v(ks,i,j)
             atm_dens(i,j) = dens(ks,i,j)
             atm_temp(i,j) = temp(ks,i,j)
             atm_pres(i,j) = pres(ks,i,j)
             atm_qv(i,j) = qv(ks,i,j)
          enddo
          enddo

          select case ( atmos_phy_sf_type )
          case ( 'BULK' )

             call atmos_phy_sf_bulk_flux( ia, isb, ieb, ja, jsb, jeb,                  & ! [IN]
                                          atm_w(:,:), atm_u(:,:), atm_v(:,:),          & ! [IN]
                                          atm_temp(:,:), atm_pres(:,:), atm_qv(:,:),   & ! [IN]
                                          sfc_dens(:,:), sfc_temp(:,:), sfc_pres(:,:), & ! [IN]
                                          sfc_z0m(:,:), sfc_z0h(:,:), sfc_z0e(:,:),    & ! [IN]
                                          pbl_zi(:,:), z1(:,:),                        & ! [IN]
                                          sflx_mw(:,:), sflx_mu(:,:), sflx_mv(:,:),    & ! [OUT]
                                          sflx_sh(:,:), sflx_lh(:,:), sflx_qv(:,:),    & ! [OUT]
                                          u10(:,:), v10(:,:), t2(:,:), q2(:,:)         ) ! [OUT]

          case ( 'CONST' )

             call atmos_phy_sf_const_flux( ia, is, ie, ja, js, je,                            & ! [IN]
                                           atm_w(:,:), atm_u(:,:), atm_v(:,:), atm_temp(:,:), & ! [IN]
                                           z1(:,:), sfc_dens(:,:),                            & ! [IN]
                                           sflx_mw(:,:), sflx_mu(:,:), sflx_mv(:,:),          & ! [OUT]
                                           sflx_sh(:,:), sflx_lh(:,:), sflx_qv(:,:),          & ! [OUT]
                                           u10(:,:), v10(:,:)                                 ) ! [OUT]
             t2(:,:) = atm_temp(:,:)
             q2(:,:) = atm_qv(:,:)

          end select

          if ( .NOT. atmos_hydrometeor_dry ) then
             sflx_qtrc(:,:,i_qv) = sflx_qv(:,:)
          endif

       endif

       ! temtative
       !$omp parallel do private(us,sflx_pt)
       do j = js, je
       do i = is, ie
          us = max( 1.e-6_rp, sqrt( sqrt( sflx_mu(i,j)**2 + sflx_mv(i,j)**2 ) / dens(ks,i,j) ) ) ! frictional velocity

          sflx_pt = sflx_sh(i,j) / ( cpdry * dens(ks,i,j) ) * pott(ks,i,j) / temp(ks,i,j)
          sflx_pt = sign( max(abs(sflx_pt), eps), sflx_pt )

          l_mo(i,j) = - us**3 * pott(ks,i,j) / ( karman * grav * sflx_pt )
       enddo
       enddo

       call history_output

!OCL XFILL
       !$omp parallel do
       do j = js, je
       do i = is, ie
          momz_t_sf(i,j) = sflx_mw(i,j) / ( cz(ks+1,i,j) - cz(ks,i,j) )
          rhou_t_sf(i,j) = sflx_mu(i,j) / ( fz(ks,i,j) - fz(ks-1,i,j) )
          rhov_t_sf(i,j) = sflx_mv(i,j) / ( fz(ks,i,j) - fz(ks-1,i,j) )
          rhoh_sf(i,j) = sflx_sh(i,j) / ( fz(ks,i,j) - fz(ks-1,i,j) )
          dens_t_sf(i,j) = 0.0_rp
          rhot_t_sf(i,j) = 0.0_rp
       enddo
       enddo

       if ( .NOT. atmos_hydrometeor_dry ) then
          !$omp parallel do &
          !$omp private(work)
          do j = js, je
          do i = is, ie
             do iq = 1, qa
                work = sflx_qtrc(i,j,iq) / ( fz(ks,i,j) - fz(ks-1,i,j) )

                rhoq_t_sf(i,j,iq) = work
                dens_t_sf(i,j)    = dens_t_sf(i,j) + work * tracer_mass(iq)
                rhot_t_sf(i,j)    = rhot_t_sf(i,j) + work * tracer_mass(iq) * rhot(ks,i,j) / dens(ks,i,j)
             enddo
          enddo
          enddo
       endif

    endif

    !$omp parallel do
    do j = js, je
    do i = is, ie
       momz_t(ks,i,j) = momz_t(ks,i,j) + momz_t_sf(i,j)
       rhou_t(ks,i,j) = rhou_t(ks,i,j) + rhou_t_sf(i,j)
       rhov_t(ks,i,j) = rhov_t(ks,i,j) + rhov_t_sf(i,j)
       rhoh(ks,i,j) = rhoh(ks,i,j) + rhoh_sf(i,j)
    enddo
    enddo

    if ( .NOT. atmos_hydrometeor_dry ) then
       !$omp parallel do
       do j = js, je
       do i = is, ie
          do iq = 1, qa
             rhoq_t(ks,i,j,iq) = rhoq_t(ks,i,j,iq) + rhoq_t_sf(i,j,iq)
          enddo
          dens_t(ks,i,j)    = dens_t(ks,i,j)    + dens_t_sf(i,j)
          rhot_t(ks,i,j)    = rhot_t(ks,i,j)    + rhot_t_sf(i,j)
       enddo
       enddo
    endif

    if ( statistics_checktotal ) then

       if ( .NOT. atmos_hydrometeor_dry ) then
          do iq = 1, qa
             call statistics_total( ia, is, ie, ja, js, je, &
                                    sflx_qtrc(:,:,iq), 'SFLX_'//trim(tracer_name(iq)), &
                                    atmos_grid_cartesc_real_area(:,:),                 &
                                    atmos_grid_cartesc_real_totarea                    )
          enddo
       endif

       call statistics_total( ia, is, ie, ja, js, je, &
                              dens_t_sf(:,:), 'DENS_t_SF',       &
                              atmos_grid_cartesc_real_area(:,:), &
                              atmos_grid_cartesc_real_totarea    )
       call statistics_total( ia, is, ie, ja, js, je, &
                              momz_t_sf(:,:), 'MOMZ_t_SF',       &
                              atmos_grid_cartesc_real_area(:,:), &
                              atmos_grid_cartesc_real_totarea    )
       call statistics_total( ia, is, ie, ja, js, je, &
                              rhou_t_sf(:,:), 'RHOU_t_SF',       &
                              atmos_grid_cartesc_real_area(:,:), &
                              atmos_grid_cartesc_real_totarea    )
       call statistics_total( ia, is, ie, ja, js, je, &
                              rhov_t_sf(:,:), 'RHOV_t_SF',       &
                              atmos_grid_cartesc_real_area(:,:), &
                              atmos_grid_cartesc_real_totarea    )
       call statistics_total( ia, is, ie, ja, js, je, &
                              rhoh_sf(:,:), 'RHOH_SF',         &
                              atmos_grid_cartesc_real_area(:,:), &
                              atmos_grid_cartesc_real_totarea    )

       if ( .NOT. atmos_hydrometeor_dry ) then
          call statistics_total( ia, is, ie, ja, js, je, &
                                 rhot_t_sf(:,:), 'RHOT_t_SF',       &
                                 atmos_grid_cartesc_real_area(:,:), &
                                 atmos_grid_cartesc_real_totarea    )
          do iq = 1, qa
             call statistics_total( ia, is, ie, ja, js, je, &
                                    rhoq_t_sf(:,:,iq), trim(tracer_name(iq))//'_t_SF', &
                                    atmos_grid_cartesc_real_area(:,:),                 &
                                    atmos_grid_cartesc_real_totarea                    )
          enddo
       endif
    endif

    return
  end subroutine atmos_phy_sf_driver_calc_tendency

  !-----------------------------------------------------------------------------
  subroutine history_output
    use scale_const, only: &
       undef => const_undef
    use scale_file_history, only: &
       file_history_in
    use scale_atmos_hydrostatic, only: &
       barometric_law_mslp => atmos_hydrostatic_barometric_law_mslp
    use scale_atmos_grid_cartesc_real, only: &
       real_cz => atmos_grid_cartesc_real_cz
    use scale_atmos_grid_cartesc_metric, only: &
       rotc => atmos_grid_cartesc_metric_rotc
    use scale_atmos_hydrometeor, only: &
       atmos_hydrometeor_dry, &
       i_qv
    use mod_atmos_vars, only: &
       temp, &
       pres, &
       qv
    use mod_atmos_phy_sf_vars, only: &
       sfc_dens   => atmos_phy_sf_sfc_dens,   &
       sfc_pres   => atmos_phy_sf_sfc_pres,   &
       sfc_temp   => atmos_phy_sf_sfc_temp,   &
       sfc_albedo => atmos_phy_sf_sfc_albedo, &
       sfc_z0m    => atmos_phy_sf_sfc_z0m,    &
       sfc_z0h    => atmos_phy_sf_sfc_z0h,    &
       sfc_z0e    => atmos_phy_sf_sfc_z0e,    &
       sflx_mw    => atmos_phy_sf_sflx_mw,    &
       sflx_mu    => atmos_phy_sf_sflx_mu,    &
       sflx_mv    => atmos_phy_sf_sflx_mv,    &
       sflx_sh    => atmos_phy_sf_sflx_sh,    &
       sflx_lh    => atmos_phy_sf_sflx_lh,    &
       sflx_gh    => atmos_phy_sf_sflx_gh,    &
       sflx_qtrc  => atmos_phy_sf_sflx_qtrc,  &
       u10        => atmos_phy_sf_u10,        &
       v10        => atmos_phy_sf_v10,        &
       t2         => atmos_phy_sf_t2,         &
       q2         => atmos_phy_sf_q2
    implicit none

    real(RP) :: MSLP  (ia,ja) ! mean sea-level pressure [Pa]
    real(RP) :: Uabs10(ia,ja) ! 10m absolute wind [m/s]
    real(RP) :: U10m  (ia,ja) ! eastward 10m wind
    real(RP) :: V10m  (ia,ja) ! northward 10m wind

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

!OCL XFILL
    !$omp parallel do
    do j = js, je
    do i = is, ie
       uabs10(i,j) = sqrt( u10(i,j)**2 + v10(i,j)**2 )
       u10m(i,j) = u10(i,j) * rotc(i,j,1) - v10(i,j) * rotc(i,j,2)
       v10m(i,j) = u10(i,j) * rotc(i,j,2) + v10(i,j) * rotc(i,j,1)
    enddo
    enddo


    call barometric_law_mslp( ka, ks, ke, ia, is, ie, ja, js, je,  & ! [IN]
                              pres(:,:,:), temp(:,:,:), qv(:,:,:), & ! [IN]
                              real_cz(:,:,:),                      & ! [IN]
                              mslp(:,:)                            ) ! [OUT]

    call file_history_in( sfc_dens(:,:),                     'SFC_DENS',        'surface atmospheric density',          'kg/m3'   )
    call file_history_in( sfc_pres(:,:),                     'SFC_PRES',        'surface atmospheric pressure',         'Pa'      )
    call file_history_in( sfc_temp(:,:),                     'SFC_TEMP',        'surface skin temperature (merged)',    'K'       )
    call file_history_in( sfc_albedo(:,:,i_r_direct ,i_r_ir ), 'SFC_ALB_IR_dir' , 'surface albedo (IR; direct; merged)',  '1'       , fill_halo=.true. )
    call file_history_in( sfc_albedo(:,:,i_r_diffuse,i_r_ir ), 'SFC_ALB_IR_dif' , 'surface albedo (IR; diffuse; merged)',  '1'       , fill_halo=.true. )
    call file_history_in( sfc_albedo(:,:,i_r_direct ,i_r_nir), 'SFC_ALB_NIR_dir', 'surface albedo (NIR; direct; merged)',  '1'       , fill_halo=.true. )
    call file_history_in( sfc_albedo(:,:,i_r_diffuse,i_r_nir), 'SFC_ALB_NIR_dif', 'surface albedo (NIR; diffuse; merged)', '1'       , fill_halo=.true. )
    call file_history_in( sfc_albedo(:,:,i_r_direct ,i_r_vis), 'SFC_ALB_VIS_dir', 'surface albedo (VIS; direct; merged)',  '1'       , fill_halo=.true. )
    call file_history_in( sfc_albedo(:,:,i_r_diffuse,i_r_vis), 'SFC_ALB_VIS_dif', 'surface albedo (VIS; diffuse; merged)', '1'       , fill_halo=.true. )
    call file_history_in( sfc_z0m(:,:),                     'SFC_Z0M',         'roughness length (momentum)',           'm'       , fill_halo=.true. )
    call file_history_in( sfc_z0h(:,:),                     'SFC_Z0H',         'roughness length (heat)',               'm'       , fill_halo=.true. )
    call file_history_in( sfc_z0e(:,:),                     'SFC_Z0E',         'roughness length (vapor)',              'm'       , fill_halo=.true. )
    call file_history_in( sflx_mw(:,:),                     'MWFLX',           'w-momentum flux (merged)',              'kg/m/s2' )
    call file_history_in( sflx_mu(:,:),                     'MUFLX',           'u-momentum flux (merged)',              'kg/m/s2' )
    call file_history_in( sflx_mv(:,:),                     'MVFLX',           'v-momentum flux (merged)',              'kg/m/s2' )
    call file_history_in( sflx_sh(:,:),                     'SHFLX',           'sensible heat flux (merged)',           'W/m2'    , fill_halo=.true. )
    call file_history_in( sflx_lh(:,:),                     'LHFLX',           'latent heat flux (merged)',             'W/m2'    , fill_halo=.true. )
    call file_history_in( sflx_gh(:,:),                     'GHFLX',           'ground heat flux (merged)',             'W/m2'    , fill_halo=.true. )
    if ( .NOT. atmos_hydrometeor_dry ) then
       do iq = 1, qa
          call file_history_in( sflx_qtrc(:,:,iq), 'SFLX_'//trim(tracer_name(iq)),   &
                                'surface '//trim(tracer_name(iq))//' flux (merged)', &
                                'kg/m2/s' , fill_halo=.true.                         )
       enddo
    endif
    call file_history_in( uabs10(:,:), 'Uabs10', '10m absolute wind',       'm/s'  , fill_halo=.true. )
    call file_history_in( u10(:,:), 'U10',    '10m x-wind',              'm/s'  , fill_halo=.true. )
    call file_history_in( v10(:,:), 'V10',    '10m y-wind',              'm/s'  , fill_halo=.true. )
    call file_history_in( u10m(:,:), 'U10m',   '10m eastward wind',       'm/s'  , fill_halo=.true. )
    call file_history_in( v10m(:,:), 'V10m',   '10m northward wind',      'm/s'  , fill_halo=.true. )
    call file_history_in( t2(:,:), 'T2 ',    '2m air temperature',      'K'    , fill_halo=.true. )
    call file_history_in( q2(:,:), 'Q2 ',    '2m specific humidity',    'kg/kg', fill_halo=.true. )
    call file_history_in( mslp(:,:), 'MSLP',   'mean sea-level pressure', 'Pa'   , fill_halo=.true., standard_name='air_pressure_at_mean_sea_level' )

    return
  end subroutine history_output

end module mod_atmos_phy_sf_driver