scale_land_phy_slab.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
module scale_land_phy_slab
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_debug
  use scale_grid_index
  use scale_land_grid_index
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: land_phy_slab_setup
  public :: land_phy_slab

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  logical, private :: LAND_PHY_UPDATE_BOTTOM_TEMP  = .false. ! Is LAND_TEMP  updated in the lowest level?
  logical, private :: LAND_PHY_UPDATE_BOTTOM_WATER = .false. ! Is LAND_WATER updated in the lowest level?

  real(RP), private :: WATER_DENSCS

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine land_phy_slab_setup( LAND_TYPE )
    use scale_process, only: &
       prc_mpistop
    use scale_const, only: &
       dwatr => const_dwatr, &
       cl    => const_cl
    implicit none

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

    namelist / param_land_phy_slab / &
       land_phy_update_bottom_temp,   &
       land_phy_update_bottom_water

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

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[SLAB] / Categ[LAND PHY] / Origin[SCALElib]'

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

    water_denscs = dwatr * cl

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '*** Update soil temperature of bottom layer? : ', land_phy_update_bottom_temp
    if( io_l ) write(io_fid_log,*) '*** Update soil moisture    of bottom layer? : ', land_phy_update_bottom_water

    return
  end subroutine land_phy_slab_setup

  !-----------------------------------------------------------------------------
  subroutine land_phy_slab( &
       TEMP_t,       &
       WATER_t,      &
       TEMP,         &
       WATER,        &
       WaterLimit,   &
       ThermalCond,  &
       HeatCapacity, &
       WaterDiff,    &
       SFLX_GH,      &
       SFLX_prec,    &
       SFLX_evap,    &
       CDZ,          &
       dt            )
    use scale_const, only: &
       dwatr => const_dwatr
    use scale_landuse, only: &
       landuse_fact_land
    use scale_matrix, only: &
       matrix_solver_tridiagonal
    implicit none

    ! arguments
    real(RP), intent(out) :: temp_t      (lkmax,ia,ja)
    real(RP), intent(out) :: water_t     (lkmax,ia,ja)

    real(RP), intent(in)  :: temp        (lkmax,ia,ja)
    real(RP), intent(in)  :: water       (lkmax,ia,ja)
    real(RP), intent(in)  :: waterlimit  (ia,ja)
    real(RP), intent(in)  :: thermalcond (ia,ja)
    real(RP), intent(in)  :: heatcapacity(ia,ja)
    real(RP), intent(in)  :: waterdiff   (ia,ja)
    real(RP), intent(in)  :: sflx_gh     (ia,ja)
    real(RP), intent(in)  :: sflx_prec   (ia,ja)
    real(RP), intent(in)  :: sflx_evap   (ia,ja)
    real(RP), intent(in)  :: cdz         (lkmax)
    real(DP), intent(in)  :: dt

    ! work
    real(RP) :: temp1 (lkmax,ia,ja)
    real(RP) :: water1(lkmax,ia,ja)

    real(RP) :: land_denscs(lkmax,ia,ja)
    real(RP) :: thermaldiff(lkmax,ia,ja)

!    real(RP) :: RUNOFF(IA,JA)

    real(RP) :: u(lkmax,ia,ja)
    real(RP) :: m(lkmax,ia,ja)
    real(RP) :: l(lkmax,ia,ja)
    real(RP) :: v(lkmax,ia,ja)

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

    if( io_l ) write(io_fid_log,*) '*** Land physics step: Slab'

    ! Solve diffusion of soil moisture (tridiagonal matrix)
    do j = js, je
    do i = is, ie
      l(lks,i,j) = 0.0_rp
      u(lks,i,j) = -2.0_rp * waterdiff(i,j) / ( cdz(lks) * ( cdz(lks) + cdz(lks+1) ) ) * dt
      l(lke,i,j) = -2.0_rp * waterdiff(i,j) / ( cdz(lke) * ( cdz(lke) + cdz(lke-1) ) ) * dt
      u(lke,i,j) = 0.0_rp

      m(lks,i,j) = 1.0_rp - l(lks,i,j) - u(lks,i,j)
      m(lke,i,j) = 1.0_rp - l(lke,i,j) - u(lke,i,j)
    end do
    end do

    do j = js, je
    do i = is, ie
    do k = lks+1, lke-1
      l(k,i,j) = -2.0_rp * waterdiff(i,j) / ( cdz(k) * ( cdz(k) + cdz(k-1) ) ) * dt
      u(k,i,j) = -2.0_rp * waterdiff(i,j) / ( cdz(k) * ( cdz(k) + cdz(k+1) ) ) * dt
      m(k,i,j) = 1.0_rp - l(k,i,j) - u(k,i,j)
    end do
    end do
    end do

    ! input from atmosphere
    do j = js, je
    do i = is, ie
      v(lks,i,j) = water(lks,i,j) + ( sflx_prec(i,j) - sflx_evap(i,j) ) / ( cdz(lks) * dwatr ) * dt
    end do
    end do

    do j = js, je
    do i = is, ie
    do k = lks+1, lke
      v(k,i,j) = water(k,i,j)
    end do
    end do
    end do

    call matrix_solver_tridiagonal( lkmax,         & ! [IN]
                                    ia, is, ie,    & ! [IN]
                                    ja, js, je,    & ! [IN]
                                    u(:,:,:), & ! [IN]
                                    m(:,:,:), & ! [IN]
                                    l(:,:,:), & ! [IN]
                                    v(:,:,:), & ! [IN]
                                    water1(:,:,:)  ) ! [OUT]

    if ( .not. land_phy_update_bottom_water ) then
      do j = js, je
      do i = is, ie
        water1(lke,i,j) = water(lke,i,j)
      end do
      end do
    endif

    ! runoff of soil moisture (vertical sum)
    do j = js, je
    do i = is, ie
!      RUNOFF(i,j) = 0.0_RP
      do k = lks, lke
!        RUNOFF(i,j) = RUNOFF(i,j) + max( WATER1(k,i,j) - WaterLimit(i,j), 0.0_RP ) * CDZ(k) * DWATR
        water1(k,i,j) = min( water1(k,i,j), waterlimit(i,j) )
      end do
    end do
    end do

    ! estimate thermal diffusivity
    do j = js, je
    do i = is, ie
    do k = lks, lke
      land_denscs(k,i,j) = ( 1.0_rp - waterlimit(i,j) ) * heatcapacity(i,j) + water_denscs * water1(k,i,j)
      thermaldiff(k,i,j) = thermalcond(i,j) / land_denscs(k,i,j)
    end do
    end do
    end do

    ! Solve diffusion of soil temperature (tridiagonal matrix)
    do j = js, je
    do i = is, ie
      l(lks,i,j) = 0.0_rp
      u(lks,i,j) = -2.0_rp * thermaldiff(lks,i,j) / ( cdz(lks) * ( cdz(lks) + cdz(lks+1) ) ) * dt
      l(lke,i,j) = -2.0_rp * thermaldiff(lke,i,j) / ( cdz(lke) * ( cdz(lke) + cdz(lke-1) ) ) * dt
      u(lke,i,j) = 0.0_rp

      m(lks,i,j) = 1.0_rp - l(lks,i,j) - u(lks,i,j)
      m(lke,i,j) = 1.0_rp - l(lke,i,j) - u(lke,i,j)
    end do
    end do

    do j = js, je
    do i = is, ie
    do k = lks+1, lke-1
      l(k,i,j) = -2.0_rp * thermaldiff(k,i,j) / ( cdz(k) * ( cdz(k) + cdz(k-1) ) ) * dt
      u(k,i,j) = -2.0_rp * thermaldiff(k,i,j) / ( cdz(k) * ( cdz(k) + cdz(k+1) ) ) * dt
      m(k,i,j) = 1.0_rp - l(k,i,j) - u(k,i,j)
    end do
    end do
    end do

    ! input from atmosphere
    do j = js, je
    do i = is, ie
      v(lks,i,j) = temp(lks,i,j) - sflx_gh(i,j) / ( land_denscs(lks,i,j) * cdz(lks) ) * dt
    end do
    end do

    do j = js, je
    do i = is, ie
    do k = lks+1, lke
      v(k,i,j) = temp(k,i,j)
    end do
    end do
    end do

    call matrix_solver_tridiagonal( lkmax,        & ! [IN]
                                    ia, is, ie,   & ! [IN]
                                    ja, js, je,   & ! [IN]
                                    u(:,:,:), & ! [IN]
                                    m(:,:,:), & ! [IN]
                                    l(:,:,:), & ! [IN]
                                    v(:,:,:), & ! [IN]
                                    temp1(:,:,:)  ) ! [OUT]

    if ( .not. land_phy_update_bottom_temp ) then
      do j = js, je
      do i = is, ie
        temp1(lke,i,j) = temp(lke,i,j)
      end do
      end do
    endif

    ! calculate tendency
    do j = js, je
    do i = is, ie
    do k = lks, lke
      if( landuse_fact_land(i,j) > 0.0_rp ) then
        temp_t(k,i,j) = ( temp1(k,i,j) - temp(k,i,j) ) / dt
        water_t(k,i,j) = ( water1(k,i,j) - water(k,i,j) ) / dt
      else
        temp_t(k,i,j) = 0.0_rp
        water_t(k,i,j) = 0.0_rp
      end if
    end do
    end do
    end do

    return
  end subroutine land_phy_slab

end module scale_land_phy_slab