mod_land_driver.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module mod_land_driver
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_land_grid_cartesc_index
  use scale_tracer
  use scale_cpl_sfc_index
 
  use scale_const, only: &
     i_lw => const_i_lw, &
     i_sw => const_i_sw
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: land_driver_setup
  public :: land_driver_calc_tendency
  public :: land_driver_update
  public :: land_surface_get
  public :: land_surface_set
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine land_driver_setup
    use scale_prc, only: &
       prc_abort
    use mod_land_admin, only: &
       land_do, &
       land_dyn_type, &
       land_sfc_type, &
       snow_type
    use scale_land_dyn_bucket, only: &
       land_dyn_bucket_setup
    use scale_land_phy_snow_ky90, only: &
       land_phy_snow_ky90_setup
    use scale_cpl_phy_sfc_skin, only: &
       cpl_phy_sfc_skin_setup
    use scale_cpl_phy_sfc_fixed_temp, only: &
       cpl_phy_sfc_fixed_temp_setup
    implicit none
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("LAND_driver_setup",*) 'Setup'
 
    if ( land_do ) then
 
       select case ( land_dyn_type )
       case ( 'BUCKET' )
          call land_dyn_bucket_setup
       case ( 'INIT' )
          ! do nothing
       case default
          log_error("LAND_driver_setup",*) 'LAND_DYN_TYPE is invalid: ', trim(land_dyn_type)
          call prc_abort
       end select
 
       select case ( land_sfc_type )
       case ( 'SKIN' )
          call cpl_phy_sfc_skin_setup
       case ( 'FIXED-TEMP' )
          call cpl_phy_sfc_fixed_temp_setup
       case default
          log_error("LAND_driver_setup",*) 'LAND_SFC_TYPE is invalid: ', trim(land_sfc_type)
          call prc_abort
       end select
 
       select case ( snow_type )
       case ( 'NONE', 'OFF' )
       case ( 'KY90' )
          log_warn("LAND_driver_setup",*) 'SNOW model is enabled'
          log_warn("LAND_driver_setup",*) 'SNOW model is on experimental stage.'
          log_warn("LAND_driver_setup",*) 'Use this with your own risk.'
          call land_phy_snow_ky90_setup
       case default
          log_error("LAND_driver_setup",*) 'SNOW_TYPE is invalid: ', trim(snow_type)
          call prc_abort
       end select
 
    end if
 
    return
  end subroutine land_driver_setup
 
  !-----------------------------------------------------------------------------
  subroutine land_driver_calc_tendency( force )
    use scale_const, only: &
       tem00 => const_tem00
    use scale_time, only: &
       dt => time_dtsec_land
    use scale_file_history, only: &
       file_history_in
    use scale_atmos_grid_cartesc_real, only: &
       real_z1 => atmos_grid_cartesc_real_z1
    use scale_topography, only: &
       tansl_x => topography_tansl_x, &
       tansl_y => topography_tansl_y
    use scale_atmos_hydrometeor, only: &
       hydrometeor_lhv => atmos_hydrometeor_lhv, &
       hydrometeor_lhs => atmos_hydrometeor_lhs, &
       atmos_hydrometeor_dry,                    &
       cv_water,                                 &
       cv_ice,                                   &
       lhf,                                      &
       i_qv
    use scale_land_grid_cartesc, only: &
       lcz => land_grid_cartesc_cz, &
       cdz => land_grid_cartesc_cdz
    use scale_land_phy_snow_ky90, only: &
       land_phy_snow_ky90
    use scale_land_phy_snow_diagnos, only: &
       land_phy_snow_diags
    use scale_cpl_phy_sfc_skin, only: &
       cpl_phy_sfc_skin
    use scale_cpl_phy_sfc_fixed_temp, only: &
       cpl_phy_sfc_fixed_temp
    use scale_bulkflux, only: &
       bulkflux_diagnose_scales
    use mod_atmos_admin, only: &
       atmos_sw_phy_ch
    use mod_atmos_phy_ch_driver, only: &
       atmos_phy_ch_driver_land_flux
    use mod_land_admin, only: &
       land_sfc_type, &
       snow_type
    use mod_land_vars, only: &
       i_waterlimit,      &
       i_watercritical,   &
       i_stomataresist,   &
       i_thermalcond,     &
       i_heatcapacity,    &
       i_waterdiff,       &
       i_alblw,           &
       i_albsw,           &
       i_z0m,             &
       i_z0h,             &
       i_z0e,             &
       snow_flag,         &
       land_property,     &
       land_temp,         &
       land_water,        &
       land_ice,          &
       land_sfc_temp,     &
       land_sfc_albedo,   &
       snow_sfc_temp,     &
       snow_swe,          &
       snow_depth,        &
       snow_dzero,        &
       snow_nosnowsec,    &
       land_temp_t,       &
       land_water_t,      &
       land_ice_t,        &
       land_sflx_gh,      &
       land_sflx_water,   &
       land_sflx_engi,    &
       land_sflx_mw,      &
       land_sflx_mu,      &
       land_sflx_mv,      &
       land_sflx_sh,      &
       land_sflx_lh,      &
       land_sflx_qtrc,    &
       land_u10,          &
       land_v10,          &
       land_t2,           &
       land_q2,           &
       land_ustar,        &
       land_tstar,        &
       land_qstar,        &
       land_wstar,        &
       land_rlmo,         &
       soil_ustar,        &
       soil_tstar,        &
       soil_qstar,        &
       soil_wstar,        &
       soil_rlmo,         &
       snow_ustar,        &
       snow_tstar,        &
       snow_qstar,        &
       snow_wstar,        &
       snow_rlmo,         &
       atmos_temp,        &
       atmos_pres,        &
       atmos_u,           &
       atmos_v,           &
       atmos_dens,        &
       atmos_qv,          &
       atmos_pbl,         &
       atmos_sfc_dens,    &
       atmos_sfc_pres,    &
       atmos_sflx_rad_dn, &
       atmos_sflx_water,  &
       atmos_sflx_engi
    use scale_landuse, only: &
       exists_land => landuse_exists_land
    implicit none
 
    logical, intent(in) :: force
 
    ! parameters
    real(rp), parameter :: beta_max = 1.0_rp
 
    ! works
    real(rp) :: snow_qvef (lia,lja)
    real(rp) :: land_wstr (lia,lja)
    real(rp) :: land_qvef (lia,lja)
    real(rp) :: land_tc_dz(lia,lja)
    real(rp) :: sflx_qv   (lia,lja)
    real(rp) :: sflx_engi (lia,lja)
    real(rp) :: lh        (lia,lja) ! latent heat of vaporization [J/kg]
    real(rp) :: atmos_w   (lia,lja)
    real(rp) :: total
 
    ! for snow
    real(rp) :: snow_albedo         (lia,lja,2)
    real(rp) :: snow_atmos_sflx_sh  (lia,lja)
    real(rp) :: snow_atmos_sflx_lh  (lia,lja)
    real(rp) :: snow_atmos_sflx_gh  (lia,lja)
    real(rp) :: snow_atmos_sflx_qv  (lia,lja)
    real(rp) :: snow_land_sflx_gh   (lia,lja)
    real(rp) :: snow_land_sflx_water(lia,lja)
    real(rp) :: snow_land_sflx_engi (lia,lja)
    real(rp) :: snow_frac           (lia,lja)
 
    real(rp) :: snow_atmos_sflx_mw  (lia,lja)
    real(rp) :: snow_atmos_sflx_mu  (lia,lja)
    real(rp) :: snow_atmos_sflx_mv  (lia,lja)
    real(rp) :: snow_u10            (lia,lja)
    real(rp) :: snow_v10            (lia,lja)
    real(rp) :: snow_t2             (lia,lja)
    real(rp) :: snow_q2             (lia,lja)
 
    ! monitor
    !real(RP) :: MONIT_WCONT0        (LIA,LJA)
    !real(RP) :: MONIT_WCONT1        (LIA,LJA)
    !real(RP) :: MONIT_ENG0          (LIA,LJA)
    !real(RP) :: MONIT_ENG1          (LIA,LJA)
    !
    !real(RP) :: MONIT_SNOW_heat     (LIA,LJA)
    !real(RP) :: MONIT_SNOW_water    (LIA,LJA)
    !real(RP) :: MONIT_LAND_heat     (LIA,LJA)
    !real(RP) :: MONIT_LAND_water    (LIA,LJA)
 
    integer :: k, i, j, iq, idir
    !---------------------------------------------------------------------------
 
    call prof_rapstart('LND_CalcTend', 1)
 
    !########## Get Surface Boundary from coupler ##########
    call land_surface_get
 
    !########## reset tendencies ##########
!OCL XFILL
    !$omp parallel do
    do j = ljs, lje
    do i = lis, lie
    do k = lks, lke
       land_temp_t(k,i,j) = 0.0_rp
       land_water_t(k,i,j) = 0.0_rp
       land_ice_t(k,i,j) = 0.0_rp
    enddo
    enddo
    enddo
!OCL XFILL
    do iq = 1, qa
    !$omp parallel do
    do j  = ljs, lje
    do i  = lis, lie
       land_sflx_qtrc(i,j,iq) = 0.0_rp
    enddo
    enddo
    enddo
 
    !$omp parallel do
    do j = ljs, lje
    do i = lis, lie
       atmos_w(i,j) = atmos_u(i,j) * tansl_x(i,j) + atmos_v(i,j) * tansl_y(i,j)
    end do
    end do
 
    if ( snow_flag ) then
       !------------------------------------------------------------------------
 
!OCL XFILL
       !$omp parallel do
       do j = ljs, lje
       do i = lis, lie
          ! This is for debug---adachi start
          !if(( int(SNOW_frac(i,j)) == 1 ).and.( abs(SNOW_SFC_TEMP(i,j)-LAND_SFC_TEMP(i,j))/=0 ))then
          !   LOG_ERROR("LAND_driver_calc_tendency",*) "Error please check SNOW_SFC_TEMP routine"
          !   call PRC_abort
          !endif
          ! This is for debug---adachi end
          snow_sfc_temp(i,j) = land_sfc_temp(i,j)
       end do
       end do
 
       select case ( snow_type )
       case ( 'KY90' )
          ! accumulation and melt of snow if there is snow
 
          !MONIT_WCONT0 = 0.0_RP
          !call monitor_snow_water(  SNOW_Depth          (:,:),   & ! [IN]
          !                          SNOW_Dzero          (:,:),   & ! [IN]
          !                          MONIT_WCONT0        (:,:)    ) ! [OUT]
 
          call land_phy_snow_ky90( lia, lis, lie, lja, ljs, lje, &
                                   atmos_sflx_water(:,:), atmos_sflx_engi(:,:),      & ! [IN]
                                   atmos_pres(:,:), atmos_temp(:,:), atmos_qv(:,:),  & ! [IN]
                                   atmos_w(:,:), atmos_u(:,:), atmos_v(:,:),         & ! [IN]
                                   atmos_sfc_dens(:,:),                              & ! [IN]
                                   atmos_sflx_rad_dn(:,:,:,:),                       & ! [IN]
                                   exists_land(:,:), dt,                             & ! [IN]
                                   snow_sfc_temp(:,:), snow_swe(:,:),                & ! [INOUT]
                                   snow_depth(:,:), snow_dzero(:,:),                 & ! [INOUT]
                                   snow_nosnowsec(:,:),                              & ! [INOUT]
                                   snow_albedo(:,:,:),                               & ! [OUT]
                                   snow_atmos_sflx_sh(:,:),                          & ! [OUT]
                                   snow_atmos_sflx_lh(:,:), snow_atmos_sflx_qv(:,:), & ! [OUT]
                                   sflx_engi(:,:),                                   & ! [OUT]
                                   snow_atmos_sflx_gh(:,:), snow_land_sflx_gh(:,:),  & ! [OUT]
                                   snow_land_sflx_water(:,:),                        & ! [OUT]
                                   snow_frac(:,:)                         ) ! [OUT]
 
!OCL XFILL
          !$omp parallel do
          do j = ljs, lje
          do i = lis, lie
          if ( exists_land(i,j) ) then
             snow_land_sflx_engi(i,j) = atmos_sflx_engi(i,j) & ! internal energy of precipitation
                                      - sflx_engi(i,j)         ! internal energy of evapolation
          end if
          enddo
          enddo
       end select
 
!OCL XFILL
       !call monitor_snow_water(  SNOW_Depth          (:,:),   & ! [IN]
       !                          SNOW_Dzero          (:,:),   & ! [IN]
       !                          MONIT_WCONT1        (:,:)    ) ! [OUT]
 
       !call monitor_land_regidual( ATMOS_SFLX_water    (:,:),   & ! [IN] ! downward at surface
       !                            ATMOS_SFLX_ENGI     (:,:),   & ! [IN] ! downward at surface
       !                            SNOW_ATMOS_SFLX_evap(:,:),   & ! [IN] ! upward   at surface
       !                            SNOW_LAND_SFLX_water(:,:),   & ! [IN] ! downward at bottom
       !                            MONIT_WCONT0        (:,:),   & ! [IN]
       !                            MONIT_WCONT1        (:,:),   & ! [IN]
       !                            MONIT_SNOW_water    (:,:)    ) ! [OUT]
 
!OCL XFILL
       !$omp parallel do
       do j = ljs, lje
       do i = lis, lie
          snow_qvef(i,j) = 1.0_rp ! tentative
       end do
       end do
 
       ! momentum fluxes and diagnostic variables above snowpack
       call land_phy_snow_diags( lia, lis, lie, lja, ljs, lje, &
                                 snow_frac(:,:),                                                & ! [IN]
                                 atmos_temp(:,:), atmos_pres(:,:),                              & ! [IN]
                                 atmos_w(:,:), atmos_u(:,:), atmos_v(:,:),                      & ! [IN]
                                 atmos_dens(:,:), atmos_qv(:,:),                                & ! [IN]
                                 real_z1(:,:), atmos_pbl(:,:),                                  & ! [IN]
                                 atmos_sfc_dens(:,:), atmos_sfc_pres(:,:), snow_sfc_temp(:,:), & ! [IN]
                                 snow_qvef(:,:),                                                & ! [IN]
                                 land_property(:,:,i_z0m),                                      & ! [IN]
                                 land_property(:,:,i_z0h),                                      & ! [IN]
                                 land_property(:,:,i_z0e),                                      & ! [IN]
                                 snow_atmos_sflx_mw(:,:),                                       & ! [OUT]
                                 snow_atmos_sflx_mu(:,:),                                       & ! [OUT]
                                 snow_atmos_sflx_mv(:,:),                                       & ! [OUT]
                                 snow_ustar(:,:), snow_tstar(:,:), snow_qstar(:,:),             & ! [OUT]
                                 snow_wstar(:,:),                                               & ! [OUT]
                                 snow_rlmo(:,:),                                                & ! [OUT]
                                 snow_u10(:,:), snow_v10(:,:),                                  & ! [OUT]
                                 snow_t2(:,:), snow_q2(:,:)                                     ) ! [OUT]
 
       call file_history_in( snow_frac(:,:),      'LAND_SNOW_frac',    'Snow fraction on land subgrid',           '1'      )
       call file_history_in( snow_albedo(:,:,i_sw), 'LAND_SNOW_ALB_SW',  'Snow surface albedo (short wave)',        '1'      )
       call file_history_in( snow_albedo(:,:,i_lw), 'LAND_SNOW_ALB_LW',  'Snow surface albedo (long wave)',         '1'      )
       call file_history_in( snow_atmos_sflx_sh(:,:),      'LAND_SNOW_SFLX_SH', 'Snow surface sensible heat flux',         'J/m2/s' )
       call file_history_in( snow_atmos_sflx_lh(:,:),      'LAND_SNOW_SFLX_LH', 'Snow surface latent heat flux',           'J/m2/s' )
       call file_history_in( snow_atmos_sflx_gh(:,:),      'LAND_SNOW_SFLX_GH', 'Snowpack received heat flux',             'J/m2/s' )
       call file_history_in( snow_atmos_sflx_mw(:,:),      'LAND_SNOW_SFLX_MW', 'Snow surface w-momentum flux',            'J/m2/s' )
       call file_history_in( snow_atmos_sflx_mu(:,:),      'LAND_SNOW_SFLX_MU', 'Snow surface u-momentum flux',            'J/m2/s' )
       call file_history_in( snow_atmos_sflx_mv(:,:),      'LAND_SNOW_SFLX_MV', 'Snow surface v-momentum flux',            'J/m2/s' )
       call file_history_in( snow_u10(:,:),      'LAND_SNOW_U10',     'Wind velocity u at 10 m on snow surface', 'm/s'    )
       call file_history_in( snow_v10(:,:),      'LAND_SNOW_V10',     'Wind velocity v at 10 m on snow surface', 'm/s'    )
       call file_history_in( snow_t2(:,:),      'LAND_SNOW_T2',      'Air temperature at 2m on snow surface',   'K'      )
       call file_history_in( snow_q2(:,:),      'LAND_SNOW_Q2',      'Specific humidity at 2m on snow surface', 'kg/kg'  )
 
       call file_history_in( snow_land_sflx_gh(:,:), 'LAND_SNOW_LAND_SFLX_GH',    'land surface ground heat flux under snow',     'J/m2/s'  )
       call file_history_in( snow_land_sflx_water(:,:), 'LAND_SNOW_LAND_SFLX_water', 'land surface water mass flux under snow',      'kg/m2/s' )
       call file_history_in( snow_land_sflx_engi(:,:), 'LAND_SNOW_LAND_SFLX_ENGI',  'land surface internal energy flux under snow', 'kg/m2/s' )
    endif
 
 
!OCL XFILL
    !$omp parallel do &
    !$omp private(total)
    do j = ljs, lje
    do i = lis, lie
    if ( exists_land(i,j) ) then
       total = land_water(lks,i,j) + land_ice(lks,i,j)
       land_wstr(i,j) = total * cdz(lks) &
                      + dt * ( atmos_sflx_water(i,j) &
                             + max( 0.0_rp, 2.0_rp * land_property(i,j,i_waterdiff) &
                                    * ( land_water(lks+1,i,j) - land_water(lks,i,j) ) / ( lcz(lks) + lcz(lks+1) ) ) )
       if ( atmos_hydrometeor_dry ) then
          land_qvef(i,j) = 0.0_rp
       else
          land_qvef(i,j) = min( total / land_property(i,j,i_watercritical), beta_max )
       end if
 
       ! eq.(12) in Merlin et al.(2011) but simplified P=0.5 used
       !sw = 0.5_RP + sign(0.5_RP,LAND_WATER(LKS,i,j)-LAND_PROPERTY(i,j,I_WaterCritical)) ! if W > Wc, sw = 1
       !LAND_QVEF(i,j) = (        sw ) * 1.0_RP &
       !               + ( 1.0_RP-sw ) * sqrt( 0.5_RP - 0.5_RP * cos( PI * LAND_WATER(LKS,i,j) / LAND_PROPERTY(i,j,I_WaterCritical) ) )
 
       land_tc_dz(i,j) = land_property(i,j,i_thermalcond) / lcz(lks)
    end if
    end do
    end do
 
 
    !------------------------------------------------------------------------
 
 
    !$omp parallel do
    do j = ljs, lje
    do i = lis, lie
    if ( exists_land(i,j) ) then
       if ( land_ice(lks,i,j) > 0.0_rp ) then
          call hydrometeor_lhs( land_sfc_temp(i,j), lh(i,j) )
       else
          call hydrometeor_lhv( land_sfc_temp(i,j), lh(i,j) )
       end if
    end if
    end do
    end do
 
 
    select case ( land_sfc_type )
    case ( 'SKIN' )
!OCL XFILL
       !$omp parallel do
       do j = ljs, lje
       do i = lis, lie
       if ( exists_land(i,j) ) then
          land_sfc_albedo(i,j,i_r_direct ,i_r_ir ) = land_property(i,j,i_alblw)
          land_sfc_albedo(i,j,i_r_diffuse,i_r_ir ) = land_property(i,j,i_alblw)
          land_sfc_albedo(i,j,i_r_direct ,i_r_nir) = land_property(i,j,i_albsw)
          land_sfc_albedo(i,j,i_r_diffuse,i_r_nir) = land_property(i,j,i_albsw)
          land_sfc_albedo(i,j,i_r_direct ,i_r_vis) = land_property(i,j,i_albsw)
          land_sfc_albedo(i,j,i_r_diffuse,i_r_vis) = land_property(i,j,i_albsw)
       end if
       end do
       end do
 
       call cpl_phy_sfc_skin( lia, lis, lie, lja, ljs, lje, &
                              atmos_temp(:,:), atmos_pres(:,:),                        & ! [IN]
                              atmos_w(:,:), atmos_u(:,:), atmos_v(:,:),                & ! [IN]
                              atmos_dens(:,:), atmos_qv(:,:),                          & ! [IN]
                              lh(:,:), real_z1(:,:), atmos_pbl(:,:),                   & ! [IN]
                              atmos_sfc_dens(:,:), atmos_sfc_pres(:,:),                & ! [IN]
                              atmos_sflx_rad_dn(:,:,:,:),                              & ! [IN]
                              land_temp(lks,:,:), land_wstr(:,:), land_qvef(:,:),      & ! [IN]
                              land_sfc_albedo(:,:,:,:),                                & ! [IN]
                              land_property(:,:,i_stomataresist),                      & ! [IN]
                              land_tc_dz(:,:),                                         & ! [IN]
                              land_property(:,:,i_z0m),                                & ! [IN]
                              land_property(:,:,i_z0h),                                & ! [IN]
                              land_property(:,:,i_z0e),                                & ! [IN]
                              exists_land(:,:), dt,                                    & ! [IN]
                              'LAND',                                                  & ! [IN]
                              land_sfc_temp(:,:),                                      & ! [INOUT]
                              land_sflx_mw(:,:), land_sflx_mu(:,:), land_sflx_mv(:,:), & ! [OUT]
                              land_sflx_sh(:,:), land_sflx_lh(:,:), sflx_qv(:,:),      & ! [OUT]
                              land_sflx_gh(:,:),                                       & ! [OUT]
                              soil_ustar(:,:), soil_tstar(:,:), soil_qstar(:,:),       & ! [OUT]
                              soil_wstar(:,:),                                         & ! [OUT]
                              soil_rlmo(:,:),                                          & ! [OUT]
                              land_u10(:,:), land_v10(:,:), land_t2(:,:), land_q2(:,:) ) ! [OUT]
 
    case ( 'FIXED-TEMP' )
!OCL XFILL
       !$omp parallel do
       do j = ljs, lje
       do i = lis, lie
       if ( exists_land(i,j) ) then
          land_sfc_temp(i,j) = land_temp(lks,i,j)
       end if
       end do
       end do
!OCL XFILL
       !$omp parallel do
       do j = ljs, lje
       do i = lis, lie
       if ( exists_land(i,j) ) then
          land_sfc_albedo(i,j,i_r_direct ,i_r_ir ) = land_property(i,j,i_alblw)
          land_sfc_albedo(i,j,i_r_diffuse,i_r_ir ) = land_property(i,j,i_alblw)
          land_sfc_albedo(i,j,i_r_direct ,i_r_nir) = land_property(i,j,i_albsw)
          land_sfc_albedo(i,j,i_r_diffuse,i_r_nir) = land_property(i,j,i_albsw)
          land_sfc_albedo(i,j,i_r_direct ,i_r_vis) = land_property(i,j,i_albsw)
          land_sfc_albedo(i,j,i_r_diffuse,i_r_vis) = land_property(i,j,i_albsw)
       end if
       end do
       end do
 
       call cpl_phy_sfc_fixed_temp( lia, lis, lie, lja, ljs, lje, &
                                    atmos_temp(:,:), atmos_pres(:,:),                        & ! [IN]
                                    atmos_w(:,:), atmos_u(:,:), atmos_v(:,:),                & ! [IN]
                                    atmos_dens(:,:), atmos_qv(:,:), lh(:,:),                 & ! [IN]
                                    real_z1(:,:), atmos_pbl(:,:),                            & ! [IN]
                                    atmos_sfc_dens(:,:), atmos_sfc_pres(:,:),                & ! [IN]
                                    atmos_sflx_rad_dn(:,:,:,:),                              & ! [IN]
                                    land_sfc_temp(:,:), land_wstr(:,:), land_qvef(:,:),      & ! [IN]
                                    land_sfc_albedo(:,:,:,:),                                & ! [IN]
                                    land_property(:,:,i_stomataresist),                      & ! [IN]
                                    land_property(:,:,i_z0m),                                & ! [IN]
                                    land_property(:,:,i_z0h),                                & ! [IN]
                                    land_property(:,:,i_z0e),                                & ! [IN]
                                    exists_land(:,:), dt,                                    & ! [IN]
                                    land_sflx_mw(:,:), land_sflx_mu(:,:), land_sflx_mv(:,:), & ! [OUT]
                                    land_sflx_sh(:,:), land_sflx_lh(:,:), sflx_qv(:,:),      & ! [OUT]
                                    land_sflx_gh(:,:),                                       & ! [OUT]
                                    soil_ustar(:,:), soil_tstar(:,:), soil_qstar(:,:),       & ! [OUT]
                                    soil_wstar(:,:),                                         & ! [OUT]
                                    soil_rlmo(:,:),                                          & ! [OUT]
                                    land_u10(:,:), land_v10(:,:),                            & ! [OUT]
                                    land_t2(:,:), land_q2(:,:)                               ) ! [OUT]
    end select
 
    !$omp parallel do
    do j = ljs, lje
    do i = lis, lie
    if ( exists_land(i,j) ) then
       if ( land_ice(lks,i,j) > 0.0_rp ) then
          sflx_engi(i,j) = ( cv_ice * land_sfc_temp(i,j) - lhf ) * sflx_qv(i,j)
       else
          sflx_engi(i,j) = cv_water * land_sfc_temp(i,j) * sflx_qv(i,j)
       end if
    end if
    end do
    end do
 
    ! LAND_SFLX_* are positive for downward
!OCL XFILL
    !$omp parallel do
    do j = ljs, lje
    do i = lis, lie
    if ( exists_land(i,j) ) then
       land_sflx_water(i,j) = atmos_sflx_water(i,j) - sflx_qv(i,j)
       land_sflx_engi(i,j) = atmos_sflx_engi(i,j) & ! internal energy of precipitation
                            - sflx_engi(i,j)         ! internal energy of evapolation or sublimation
    end if
    end do
    end do
 
    if ( snow_flag ) then
 
       call file_history_in( land_sflx_mw(:,:), 'SOIL_SFLX_MW',  'soil surface w-momentum flux (upward)',    'kg/m2/s' )
       call file_history_in( land_sflx_mu(:,:), 'SOIL_SFLX_MU',  'soil surface u-momentum flux (upward)',    'kg/m2/s' )
       call file_history_in( land_sflx_mv(:,:), 'SOIL_SFLX_MV',  'soil surface v-momentum flux (upward)',    'kg/m2/s' )
       call file_history_in( land_sflx_sh(:,:), 'SOIL_SFLX_SH',  'soil surface sensible heat flux (upward)', 'J/m2/s'  )
       call file_history_in( land_sflx_lh(:,:), 'SOIL_SFLX_LH',  'soil surface latent heat flux (upward)',   'J/m2/s'  )
       call file_history_in( land_u10(:,:), 'LAND_SOIL_U10', 'Wind velocity u at 10 m on soil surface',  'm/s'     )
       call file_history_in( land_v10(:,:), 'LAND_SOIL_V10', 'Wind velocity v at 10 m on soil surface',  'm/s'     )
       call file_history_in( land_t2(:,:), 'LAND_SOIL_T2',  'Air temperature at 2m on soil surface',    'K'       )
       call file_history_in( land_q2(:,:), 'LAND_SOIL_Q2',  'Specific humidity at 2m on soil surface',  'kg/kg'   )
 
       ! marge land surface and snow surface !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!OCL XFILL
       !$omp parallel do
       do j = ljs, lje
       do i = lis, lie
       if ( exists_land(i,j) ) then
          land_sfc_temp(i,j) = (        snow_frac(i,j) ) * snow_sfc_temp(i,j) &
                             + ( 1.0_rp-snow_frac(i,j) ) * land_sfc_temp(i,j)
 
          do idir = i_r_direct, i_r_diffuse
             land_sfc_albedo(i,j,idir,i_r_ir ) = (        snow_frac(i,j) ) * snow_albedo(i,j,i_lw)         &
                                               + ( 1.0_rp-snow_frac(i,j) ) * land_sfc_albedo(i,j,idir,i_r_ir)
             land_sfc_albedo(i,j,idir,i_r_nir) = (        snow_frac(i,j) ) * snow_albedo(i,j,i_sw)         &
                                               + ( 1.0_rp-snow_frac(i,j) ) * land_sfc_albedo(i,j,idir,i_r_nir)
             land_sfc_albedo(i,j,idir,i_r_vis) = (        snow_frac(i,j) ) * snow_albedo(i,j,i_sw)         &
                                               + ( 1.0_rp-snow_frac(i,j) ) * land_sfc_albedo(i,j,idir,i_r_vis)
          enddo
 
          ! flux to the soil
          land_sflx_gh(i,j) = (        snow_frac(i,j) ) * snow_land_sflx_gh(i,j) &
                               + ( 1.0_rp-snow_frac(i,j) ) *      land_sflx_gh(i,j)
          land_sflx_water(i,j) = (        snow_frac(i,j) ) * snow_land_sflx_water(i,j) &
                               + ( 1.0_rp-snow_frac(i,j) ) *      land_sflx_water(i,j)
          land_sflx_engi(i,j) = (        snow_frac(i,j) ) * snow_land_sflx_engi(i,j) &
                               + ( 1.0_rp-snow_frac(i,j) ) *      land_sflx_engi(i,j)
          ! flux to the atmosphere
          land_sflx_mw(i,j) = (        snow_frac(i,j) ) * snow_atmos_sflx_mw(i,j) &
                            + ( 1.0_rp-snow_frac(i,j) ) *       land_sflx_mw(i,j)
          land_sflx_mu(i,j) = (        snow_frac(i,j) ) * snow_atmos_sflx_mu(i,j) &
                            + ( 1.0_rp-snow_frac(i,j) ) *       land_sflx_mu(i,j)
          land_sflx_mv(i,j) = (        snow_frac(i,j) ) * snow_atmos_sflx_mv(i,j) &
                            + ( 1.0_rp-snow_frac(i,j) ) *       land_sflx_mv(i,j)
          land_sflx_sh(i,j) = (        snow_frac(i,j) ) * snow_atmos_sflx_sh(i,j) &
                            + ( 1.0_rp-snow_frac(i,j) ) *       land_sflx_sh(i,j)
          land_sflx_lh(i,j) = (        snow_frac(i,j) ) * snow_atmos_sflx_lh(i,j) &
                            + ( 1.0_rp-snow_frac(i,j) ) *       land_sflx_lh(i,j)
               sflx_qv(i,j) = (        snow_frac(i,j) ) * snow_atmos_sflx_qv(i,j) &
                            + ( 1.0_rp-snow_frac(i,j) ) *            sflx_qv(i,j)
          ! diagnostics
          land_u10(i,j) = (        snow_frac(i,j) ) * snow_u10(i,j) &
                        + ( 1.0_rp-snow_frac(i,j) ) * land_u10(i,j)
          land_v10(i,j) = (        snow_frac(i,j) ) * snow_v10(i,j) &
                        + ( 1.0_rp-snow_frac(i,j) ) * land_v10(i,j)
          land_t2(i,j) = (        snow_frac(i,j) ) * snow_t2(i,j) &
                        + ( 1.0_rp-snow_frac(i,j) ) * land_t2(i,j)
          land_q2(i,j) = (        snow_frac(i,j) ) * snow_q2(i,j) &
                        + ( 1.0_rp-snow_frac(i,j) ) * land_q2(i,j)
       end if
       enddo
       enddo
 
       call bulkflux_diagnose_scales( lia, lis, lie, lja, ljs, lje, &
                                      land_sflx_mw(:,:), land_sflx_mu(:,:), land_sflx_mv(:,:), & ! [IN]
                                      land_sflx_sh(:,:), sflx_qv(:,:),                         & ! [IN]
                                      atmos_sfc_dens(:,:), land_sfc_temp(:,:), atmos_pbl(:,:), & ! [IN]
                                      land_ustar(:,:), land_tstar(:,:), land_qstar(:,:),       & ! [OUT]
                                      land_wstar(:,:), land_rlmo(:,:),                         & ! [OUT]
                                      mask = exists_land(:,:)                                  ) ! [IN]
 
    end if
 
    if ( .NOT. atmos_hydrometeor_dry ) then
       !$omp parallel do
       do j = ljs, lje
       do i = lis, lie
       if ( exists_land(i,j) ) then
          land_sflx_qtrc(i,j,i_qv) = sflx_qv(i,j)
       end if
       enddo
       enddo
    end if
 
 
    ! Surface flux for chemical tracers
    if ( atmos_sw_phy_ch ) then
       call atmos_phy_ch_driver_land_flux( land_sflx_qtrc(:,:,:) ) ! [INOUT]
    endif
 
    call prof_rapend  ('LND_CalcTend', 1)
 
    !########## Set Surface Boundary to coupler ##########
    call land_surface_set( countup=.true. )
 
    return
  end subroutine land_driver_calc_tendency
 
  !-----------------------------------------------------------------------------
  subroutine land_driver_update
    use scale_time, only: &
       dt => time_dtsec_land
    use mod_land_vars, only: &
       land_property,     &
       i_waterlimit,      &
       i_thermalcond,     &
       i_heatcapacity,    &
       i_waterdiff,       &
       land_temp,         &
       land_water,        &
       land_ice,          &
       land_sflx_gh,      &
       land_sflx_water,   &
       land_sflx_engi,    &
       land_runoff,       &
       land_runoff_engi,  &
       land_sfc_temp,     &
       land_temp_t,       &
       land_water_t,      &
       land_ice_t,        &
       land_vars_check
    use scale_land_grid_cartesc, only: &
       lcdz => land_grid_cartesc_cdz
    use scale_land_dyn_bucket, only: &
       land_dyn_bucket
    use scale_landuse, only: &
       exists_land => landuse_exists_land
    use scale_time, only: &
       nowdaysec => time_nowdaysec
    use scale_file_history, only: &
       file_history_in
    use mod_land_admin, only: &
       land_dyn_type
    implicit none
 
    integer :: k, i, j
    !---------------------------------------------------------------------------
 
    call prof_rapstart('LND_Update', 2)
 
    !########## Get Surface Boundary from coupler ##########
    call land_surface_get
 
    !########## Dynamics / Update variables ##########
    select case ( land_dyn_type )
    case ( 'BUCKET' )
       call land_dyn_bucket( lkmax, lks, lke, lia, lis, lie, lja, ljs, lje, &
                             land_temp_t(:,:,:),                     & ! [IN]
                             land_water_t(:,:,:), land_ice_t(:,:,:), & ! [IN]
                             land_property(:,:,i_waterlimit),        & ! [IN]
                             land_property(:,:,i_thermalcond),       & ! [IN]
                             land_property(:,:,i_heatcapacity),      & ! [IN]
                             land_property(:,:,i_waterdiff),         & ! [IN]
                             land_sflx_gh(:,:),                      & ! [IN]
                             land_sflx_water(:,:),                   & ! [IN]
                             land_sflx_engi(:,:),                    & ! [IN]
                             exists_land(:,:), lcdz(:),              & ! [IN]
                             dt, nowdaysec,                          & ! [IN]
                             land_temp(:,:,:),                       & ! [INOUT]
                             land_water(:,:,:), land_ice(:,:,:),     & ! [INOUT]
                             land_runoff(:,:), land_runoff_engi(:,:) ) ! [OUT]
    case ( 'INIT' )
       ! Never update LAND_TEMP and LAND_WATER from initial condition
    end select
 
    !########## Negative Fixer ##########
    !$omp parallel do
    do j = ljs, lje
    do i = lis, lie
       if ( exists_land(i,j) ) then
          do k = lks, lke
             land_water(k,i,j) = max( land_water(k,i,j), 0.0_rp )
             land_ice(k,i,j) = max( land_ice(k,i,j), 0.0_rp )
          enddo
       end if
    enddo
    enddo
 
    call land_vars_check
 
    call prof_rapend  ('LND_Update', 1)
 
    return
  end subroutine land_driver_update
 
  !-----------------------------------------------------------------------------
  subroutine land_surface_get
    use mod_land_admin, only: &
       land_do
    use mod_land_vars, only: &
       atmos_temp,        &
       atmos_pres,        &
       atmos_w,           &
       atmos_u,           &
       atmos_v,           &
       atmos_dens,        &
       atmos_qv,          &
       atmos_pbl,         &
       atmos_sfc_dens,    &
       atmos_sfc_pres,    &
       atmos_sflx_rad_dn, &
       atmos_cossza,      &
       atmos_sflx_water,  &
       atmos_sflx_engi
    use mod_cpl_vars, only: &
       cpl_getatm_lnd
    implicit none
    !---------------------------------------------------------------------------
 
    call prof_rapstart('LND_SfcExch', 2)
 
    if ( land_do ) then
       call cpl_getatm_lnd( atmos_temp(:,:),     & ! [OUT]
                            atmos_pres(:,:),     & ! [OUT]
                            atmos_w(:,:),     & ! [OUT]
                            atmos_u(:,:),     & ! [OUT]
                            atmos_v(:,:),     & ! [OUT]
                            atmos_dens(:,:),     & ! [OUT]
                            atmos_qv(:,:),     & ! [OUT]
                            atmos_pbl(:,:),     & ! [OUT]
                            atmos_sfc_dens(:,:),     & ! [OUT]
                            atmos_sfc_pres(:,:),     & ! [OUT]
                            atmos_sflx_rad_dn(:,:,:,:), & ! [OUT]
                            atmos_cossza(:,:),     & ! [OUT]
                            atmos_sflx_water(:,:),     & ! [OUT]
                            atmos_sflx_engi(:,:)      ) ! [OUT]
    endif
 
    call prof_rapend  ('LND_SfcExch', 2)
 
    return
  end subroutine land_surface_get
 
  !-----------------------------------------------------------------------------
  subroutine land_surface_set( countup )
    use mod_land_admin, only: &
       land_do
    use mod_land_vars, only: &
       land_property,   &
       i_z0m,           &
       i_z0h,           &
       i_z0e,           &
       land_sfc_temp,   &
       land_sfc_albedo, &
       land_sflx_mw,    &
       land_sflx_mu,    &
       land_sflx_mv,    &
       land_sflx_sh,    &
       land_sflx_lh,    &
       land_sflx_gh,    &
       land_sflx_qtrc,  &
       land_u10,        &
       land_v10,        &
       land_t2,         &
       land_q2
    use mod_cpl_vars, only: &
       cpl_putlnd
    use scale_landuse, only: &
       exists_land => landuse_exists_land
    implicit none
 
    ! arguments
    logical, intent(in) :: countup
    !---------------------------------------------------------------------------
 
    call prof_rapstart('LND_SfcExch', 2)
 
    if ( land_do ) then
       call cpl_putlnd( land_sfc_temp(:,:),       & ! [IN]
                        land_sfc_albedo(:,:,:,:),   & ! [IN]
                        land_property(:,:,i_z0m), & ! [IN]
                        land_property(:,:,i_z0h), & ! [IN]
                        land_property(:,:,i_z0e), & ! [IN]
                        land_sflx_mw(:,:),       & ! [IN]
                        land_sflx_mu(:,:),       & ! [IN]
                        land_sflx_mv(:,:),       & ! [IN]
                        land_sflx_sh(:,:),       & ! [IN]
                        land_sflx_lh(:,:),       & ! [IN]
                        land_sflx_gh(:,:),       & ! [IN]
                        land_sflx_qtrc(:,:,:),     & ! [IN]
                        land_u10(:,:),       & ! [IN]
                        land_v10(:,:),       & ! [IN]
                        land_t2(:,:),       & ! [IN]
                        land_q2(:,:),       & ! [IN]
                        exists_land(:,:),       & ! [IN]
                        countup                     ) ! [IN]
    endif
 
    call prof_rapend  ('LND_SfcExch', 2)
 
    return
  end subroutine land_surface_set
 
end module mod_land_driver