scale_landuse.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_landuse
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_atmos_grid_cartesc_index
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: landuse_setup
  public :: landuse_calc_fact
  public :: landuse_fillhalo
  public :: landuse_write

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  real(RP), public, allocatable :: landuse_fact_ocean(:,:)
  real(RP), public, allocatable :: landuse_fact_land (:,:)
  real(RP), public, allocatable :: landuse_fact_urban(:,:)
  real(RP), public, allocatable :: landuse_fact_lake (:,:)

  logical,  public, allocatable :: landuse_exists_ocean(:,:)
  logical,  public, allocatable :: landuse_exists_land (:,:)
  logical,  public, allocatable :: landuse_exists_urban(:,:)
  logical,  public, allocatable :: landuse_exists_lake (:,:)

  real(RP), public, allocatable :: landuse_frac_land (:,:)
  real(RP), public, allocatable :: landuse_frac_urban(:,:)
  real(RP), public, allocatable :: landuse_frac_lake (:,:)

  integer,  public, parameter   :: landuse_index_ocean  =  0
  integer,  public, parameter   :: landuse_index_urban  = -1
  integer,  public, parameter   :: landuse_index_lake   = -2

  integer,  public, parameter   :: landuse_pft_nmin   = -2
  integer,  public              :: landuse_pft_nmax   = 17
  integer,  public              :: landuse_pft_mosaic = 2

  real(RP), public, allocatable :: landuse_frac_pft (:,:,:)
  integer,  public, allocatable :: landuse_index_pft(:,:,:)

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: landuse_read

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  character(len=H_LONG),  private :: landuse_in_basename  = ''
  logical,                private :: landuse_in_aggregate
  logical,                private :: landuse_in_check_coordinates = .true.      
  character(len=H_LONG),  private :: landuse_out_basename = ''
  logical,                private :: landuse_out_aggregate
  character(len=H_MID),   private :: landuse_out_title    = 'SCALE-RM LANDUSE'
  character(len=H_SHORT), private :: landuse_out_dtype    = 'DEFAULT'
  logical,                private :: landuse_allocean     = .false.
  logical,                private :: landuse_allland      = .false.
  logical,                private :: landuse_allurban     = .false.
  logical,                private :: landuse_alllake      = .false.
  logical,                private :: landuse_ignore_lake  = .false.
  logical,                private :: landuse_mosaicworld  = .false.

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine landuse_setup( &
       OCEAN_do, &
       URBAN_do, &
       LAKE_do   )
    use scale_prc, only: &
       prc_abort
    use scale_file, only: &
       file_aggregate
    implicit none

    logical, intent(in) :: OCEAN_do
    logical, intent(in) :: URBAN_do
    logical, intent(in) :: LAKE_do

    namelist / param_landuse / &
       landuse_in_basename,          &
       landuse_in_aggregate,         &
       landuse_in_check_coordinates, &
       landuse_out_basename,         &
       landuse_out_aggregate,        &
       landuse_out_dtype,            &
       landuse_pft_mosaic,           &
       landuse_pft_nmax,             &
       landuse_allocean,             &
       landuse_allland,              &
       landuse_allurban,             &
       landuse_alllake,              &
       landuse_ignore_lake,          &
       landuse_mosaicworld

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

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

    landuse_in_aggregate  = file_aggregate
    landuse_out_aggregate = file_aggregate

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_landuse,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("LANDUSE_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("LANDUSE_setup",*) 'Not appropriate names in namelist PARAM_LANDUSE. Check!'
       call prc_abort
    endif
    log_nml(param_landuse)

    allocate( landuse_frac_land(ia,ja) )
    allocate( landuse_frac_urban(ia,ja) )
    allocate( landuse_frac_lake(ia,ja) )
    landuse_frac_land(:,:) = 0.0_rp
    landuse_frac_urban(:,:) = 0.0_rp
    landuse_frac_lake(:,:) = 0.0_rp

    allocate( landuse_index_pft(ia,ja,landuse_pft_mosaic) )
    allocate( landuse_frac_pft(ia,ja,landuse_pft_mosaic) )
    landuse_frac_pft(:,:,:) = 0.0_rp
    landuse_frac_pft(:,:,1) = 1.0_rp ! tentative, mosaic is off
    landuse_index_pft(:,:,:) = 1      ! default

    allocate( landuse_fact_ocean(ia,ja) )
    allocate( landuse_fact_land(ia,ja) )
    allocate( landuse_fact_urban(ia,ja) )
    allocate( landuse_fact_lake(ia,ja) )
    landuse_fact_ocean(:,:) = 0.0_rp
    landuse_fact_land(:,:) = 0.0_rp
    landuse_fact_urban(:,:) = 0.0_rp
    landuse_fact_lake(:,:) = 0.0_rp

    allocate( landuse_exists_ocean(ia,ja) )
    allocate( landuse_exists_land(ia,ja) )
    allocate( landuse_exists_urban(ia,ja) )
    allocate( landuse_exists_lake(ia,ja) )
    landuse_exists_ocean(:,:) = .false.
    landuse_exists_land(:,:) = .false.
    landuse_exists_urban(:,:) = .false.
    landuse_exists_lake(:,:) = .false.



    if    ( landuse_allocean ) then

       log_info("LANDUSE_setup",*) 'Assume all grids are ocean'

       call landuse_calc_fact

    elseif( landuse_allland ) then

       log_info("LANDUSE_setup",*) 'Assume all grids are land'
       log_info("LANDUSE_setup",*) 'Assume land PFT is 1 (bare ground)'
       landuse_frac_land(:,:) = 1.0_rp

       call landuse_calc_fact

    elseif( landuse_allurban ) then

       log_info("LANDUSE_setup",*) 'Assume all grids are land'
       log_info("LANDUSE_setup",*) 'Assume all lands are urban'
       log_info("LANDUSE_setup",*) 'Assume land PFT is urban: ', landuse_index_urban
       landuse_frac_land(:,:)   = 1.0_rp
       landuse_frac_urban(:,:)   = 1.0_rp
       landuse_index_pft(:,:,:) = landuse_index_urban

       call landuse_calc_fact

    elseif( landuse_alllake ) then

       log_info("LANDUSE_setup",*) 'Assume all grids are land'
       log_info("LANDUSE_setup",*) 'Assume all lands are lake'
       log_info("LANDUSE_setup",*) 'Assume land PFT is lake: ', landuse_index_lake
       landuse_frac_land(:,:)   = 1.0_rp
       landuse_frac_lake(:,:)   = 1.0_rp
       landuse_index_pft(:,:,:) = landuse_index_lake

       call landuse_calc_fact

    elseif( landuse_mosaicworld ) then

       log_info("LANDUSE_setup",*) 'Assume all grids have ocean, land, and urban'
!       LOG_INFO("LANDUSE_setup",*) 'Assume all grids have ocean, land, urban, and lake'
       log_info("LANDUSE_setup",*) 'Assume land PFT is 1 (bare ground)'
       landuse_frac_land(:,:) = 0.5_rp
       landuse_frac_urban(:,:) = 0.5_rp
!       LANDUSE_frac_lake (:,:) = 0.25_RP

       call landuse_calc_fact

    else ! default: read from file

       call landuse_read( ocean_do, urban_do, lake_do )
       call landuse_calc_fact

    endif

    return
  end subroutine landuse_setup

  !-----------------------------------------------------------------------------
  subroutine landuse_calc_fact
    implicit none

    real(RP) :: fact_soil
    integer  :: i, j
    !---------------------------------------------------------------------------

    log_newline
    log_info("LANDUSE_calc_fact",*) 'calculate landuse factor'

    ! make factors
    !$omp parallel do private(fact_soil)
    do j = 1, ja
    do i = 1, ia
       landuse_fact_ocean(i,j) = max( 1.0_rp - landuse_frac_land(i,j), 0.0_rp )
       landuse_fact_lake(i,j) = landuse_frac_land(i,j) * landuse_frac_lake(i,j)
       fact_soil = max( landuse_frac_land(i,j) - landuse_fact_lake(i,j), 0.0_rp )
       landuse_fact_urban(i,j) = fact_soil * landuse_frac_urban(i,j)
       landuse_fact_land(i,j) = max( fact_soil - landuse_fact_urban(i,j), 0.0_rp )

       if( landuse_fact_ocean(i,j) > 0.0_rp ) landuse_exists_ocean(i,j) = .true.
       if( landuse_fact_land(i,j) > 0.0_rp ) landuse_exists_land(i,j) = .true.
       if( landuse_fact_urban(i,j) > 0.0_rp ) landuse_exists_urban(i,j) = .true.
       if( landuse_fact_lake(i,j) > 0.0_rp ) landuse_exists_lake(i,j) = .true.
    enddo
    enddo

    return
  end subroutine landuse_calc_fact

  !-----------------------------------------------------------------------------
  subroutine landuse_fillhalo( FILL_BND )
    use scale_comm_cartesc, only: &
       comm_vars8, &
       comm_wait
    implicit none

    logical, intent(in), optional :: FILL_BND

    real(RP) :: temp(ia,ja)

    logical  :: FILL_BND_
    integer  :: p
    !---------------------------------------------------------------------------

    fill_bnd_ = .true.
    if ( present(fill_bnd) ) fill_bnd_ = fill_bnd

    call comm_vars8( landuse_frac_land(:,:), 1 )
    call comm_vars8( landuse_frac_lake(:,:), 2 )
    call comm_vars8( landuse_frac_urban(:,:), 3 )
    call comm_vars8( landuse_fact_ocean(:,:), 4 )
    call comm_vars8( landuse_fact_land(:,:), 5 )
    call comm_vars8( landuse_fact_urban(:,:), 6 )
    call comm_vars8( landuse_fact_lake(:,:), 7 )

    call comm_wait ( landuse_frac_land(:,:), 1, fill_bnd_ )
    call comm_wait ( landuse_frac_lake(:,:), 2, fill_bnd_ )
    call comm_wait ( landuse_frac_urban(:,:), 3, fill_bnd_ )
    call comm_wait ( landuse_fact_ocean(:,:), 4, fill_bnd_ )
    call comm_wait ( landuse_fact_land(:,:), 5, fill_bnd_ )
    call comm_wait ( landuse_fact_urban(:,:), 6, fill_bnd_ )
    call comm_wait ( landuse_fact_lake(:,:), 7, fill_bnd_ )

    do p = 1, landuse_pft_mosaic
       temp(:,:) = real(LANDUSE_index_PFT(:,:,p),kind=rp)

       call comm_vars8( landuse_frac_pft(:,:,p), 8+2*(p-1) )
       call comm_vars8( temp(:,:)  , 9+2*(p-1) )

       call comm_wait ( landuse_frac_pft(:,:,p), 8+2*(p-1), fill_bnd_ )
       call comm_wait ( temp(:,:)  , 9+2*(p-1), fill_bnd_ )

       landuse_index_pft(:,:,p) = aint(temp(:,:),kind=4)
    enddo

    return
  end subroutine landuse_fillhalo

  !-----------------------------------------------------------------------------
  subroutine landuse_read( &
       OCEAN_do, &
       URBAN_do, &
       LAKE_do   )
    use scale_file_cartesc, only: &
       file_cartesc_open,              &
       file_cartesc_read,              &
       file_cartesc_flush,             &
       file_cartesc_check_coordinates, &
       file_cartesc_close
    use scale_prc, only: &
       prc_abort
    implicit none

    logical, intent(in) :: OCEAN_do
    logical, intent(in) :: URBAN_do
    logical, intent(in) :: LAKE_do

    real(RP) :: temp(ia,ja)

    character(len=H_SHORT) :: varname

    integer  :: fid
    integer  :: p, i, j
    !---------------------------------------------------------------------------

    log_newline
    log_info("LANDUSE_read",*) 'Input landuse file '

    if ( landuse_in_basename /= '' ) then

       call file_cartesc_open( landuse_in_basename, fid, aggregate=landuse_in_aggregate )

       call file_cartesc_read( fid, 'FRAC_LAND',  'XY', landuse_frac_land(:,:) )
       call file_cartesc_read( fid, 'FRAC_LAKE',  'XY', landuse_frac_lake(:,:) )
       call file_cartesc_read( fid, 'FRAC_URBAN', 'XY', landuse_frac_urban(:,:) )

       call file_cartesc_flush( fid )

       do p = 1, landuse_pft_mosaic
          write(varname,'(A8,I1.1)') 'FRAC_PFT', p

          call file_cartesc_read( fid, varname, 'XY', landuse_frac_pft(:,:,p) )

          write(varname,'(A9,I1.1)') 'INDEX_PFT', p

          call file_cartesc_read( fid, varname, 'XY', temp(:,:) )
          call file_cartesc_flush( fid ) ! for non-blocking I/O

          do j = js, je
          do i = is, ie
             landuse_index_pft(i,j,p) = int(temp(i,j)+1.e-3_rp,kind=4)
          enddo
          enddo
       enddo

       if ( landuse_in_check_coordinates ) then
          call file_cartesc_check_coordinates( fid )
       end if

       call file_cartesc_close( fid )

       call landuse_fillhalo( fill_bnd=.false. )

       ! validity check
       !$omp parallel do
       do j = 1, ja
       do i = 1, ia
          if ( landuse_frac_land(i,j) < 0.0_rp .or. landuse_frac_land(i,j) > 1.0_rp ) then
             log_error("LANDUSE_read",*) 'LANDUSE_frac_land is invalid: ', i,j, landuse_frac_land(i,j)
             call prc_abort
          end if
          if ( landuse_frac_lake(i,j) < 0.0_rp .or. landuse_frac_lake(i,j) > 1.0_rp ) then
             log_error("LANDUSE_read",*) 'LANDUSE_frac_lake is invalid: ', i,j, landuse_frac_lake(i,j)
             call prc_abort
          end if
          if ( landuse_frac_urban(i,j) < 0.0_rp .or. landuse_frac_urban(i,j) > 1.0_rp ) then
             log_error("LANDUSE_read",*) 'LANDUSE_frac_urban is invalid: ', i,j, landuse_frac_urban(i,j)
             call prc_abort
          end if
       end do
       end do

!!$       if ( .not. OCEAN_do ) then
!!$          !$omp parallel do private(frac_ocean)
!!$          do j = 1, JA
!!$          do i = 1, IA
!!$             do p = LANDUSE_PFT_mosaic, 1, -1
!!$                frac_ocean = 1.0_RP - LANDUSE_frac_land(i,j)
!!$                if ( frac_ocean > LANDUSE_frac_PFT(i,j,p) ) then
!!$                   if ( p < LANDUSE_PFT_mosaic ) then
!!$                      LANDUSE_index_PFT(i,j,p+1) = LANDUSE_index_PFT(i,j,p)
!!$                      LANDUSE_frac_PFT (i,j,p+1) = LANDUSE_frac_PFT (i,j,p)
!!$                   end if
!!$                   LANDUSE_index_PFT(i,j,p) = LANDUSE_index_OCEAN
!!$                   LANDUSE_frac_PFT (i,j,p) = frac_ocean
!!$                else
!!$                   exit
!!$                end if
!!$             end do
!!$             LANDUSE_frac_PFT(i,j,:) = LANDUSE_frac_PFT(i,j,:) / sum( LANDUSE_frac_PFT(i,j,:) )
!!$          end do
!!$          end do
!!$          LANDUSE_frac_land(:,:) = 1.0_RP
!!$       end if

       if ( .not. urban_do ) then
          !$omp parallel do
          do j = 1, ja
          do i = 1, ia
             do p = landuse_pft_mosaic, 1, -1
                if ( landuse_frac_urban(i,j) > landuse_frac_pft(i,j,p) ) then
                   if ( p < landuse_pft_mosaic ) then
                      landuse_index_pft(i,j,p+1) = landuse_index_pft(i,j,p)
                      landuse_frac_pft(i,j,p+1) = landuse_frac_pft(i,j,p)
                   end if
                   landuse_index_pft(i,j,p) = landuse_index_urban
                   landuse_frac_pft(i,j,p) = landuse_frac_urban(i,j)
                end if
             end do
             landuse_frac_pft(i,j,:) = landuse_frac_pft(i,j,:) / sum( landuse_frac_pft(i,j,:) )
          end do
          end do
          landuse_frac_urban(:,:) = 0.0_rp
       end if

       if ( .not. lake_do ) then
          if ( landuse_ignore_lake .or. (.not. ocean_do) ) then
             !$omp parallel do
             do j = 1, ja
             do i = 1, ia
                do p = landuse_pft_mosaic, 1, -1
                   if ( landuse_frac_lake(i,j) > landuse_frac_pft(i,j,p) ) then
                      if ( p < landuse_pft_mosaic ) then
                         landuse_index_pft(i,j,p+1) = landuse_index_pft(i,j,p)
                         landuse_frac_pft(i,j,p+1) = landuse_frac_pft(i,j,p)
                      end if
                      landuse_index_pft(i,j,p) = landuse_index_lake
                      landuse_frac_pft(i,j,p) = landuse_frac_lake(i,j)
                   else
                      exit
                   end if
                end do
                landuse_frac_pft(i,j,:) = landuse_frac_pft(i,j,:) / sum( landuse_frac_pft(i,j,:) )
                landuse_frac_urban(i,j) = max( landuse_frac_urban(i,j) * ( 1.0_rp - landuse_frac_lake(i,j) ), 0.0_rp )
             end do
             end do
          else
             ! lake is assumed to be ocean
             !$omp parallel do
             do j = 1, ja
             do i = 1, ia
                if( landuse_frac_land(i,j) == 0.0_rp .and. landuse_frac_urban(i,j) > 0.0_rp)then
                   log_error("LANDUSE_read",*) 'Not appropriate original landuse (w/ Lake). Bug!',i,j
                   call prc_abort
                endif

                landuse_frac_land(i,j)  = min( max( landuse_frac_land(i,j) * (1.0_rp-landuse_frac_lake(i,j)), 0.0_rp), 1.0_rp)

                if( landuse_frac_land(i,j) == 0.0_rp .and. landuse_frac_urban(i,j) > 0.0_rp)then
                   log_error("LANDUSE_read",*) 'Not appropriate replaced landuse (w/o Lake). Bug!',i,j
                   call prc_abort
                endif

             end do
             end do
          end if
          landuse_frac_lake(:,:) = 0.0_rp
       end if

    else
       log_info_cont(*) 'landuse file is not specified.'
       log_info_cont(*) 'Assume all grids are ocean'
    endif

    return
  end subroutine landuse_read

  !-----------------------------------------------------------------------------
  subroutine landuse_write
    use scale_file_cartesc, only: &
       file_cartesc_create,    &
       file_cartesc_def_var,   &
       file_cartesc_enddef,    &
       file_cartesc_write_var, &
       file_cartesc_close
    implicit none

    real(RP) :: temp(ia,ja)

    integer                :: vid(7+landuse_pft_mosaic*2)
    character(len=H_SHORT) :: varname

    integer :: fid
    integer :: p
    !---------------------------------------------------------------------------

    if ( landuse_out_basename /= '' .and. landuse_out_basename /= landuse_in_basename ) then

       log_newline
       log_info("LANDUSE_write",*) 'Output landuse file '

       call landuse_fillhalo( fill_bnd=.false. )

       call file_cartesc_create( landuse_out_basename, landuse_out_title, landuse_out_dtype, & ! [IN]
                                 fid,                                                        & ! [OUT]
                                 haszcoord=.false., aggregate=landuse_out_aggregate          ) ! [IN]

       call file_cartesc_def_var( fid, 'FRAC_LAND'     , 'LAND fraction'          , '1', 'XY', landuse_out_dtype, vid(1), standard_name="land_area_fraction" )
       call file_cartesc_def_var( fid, 'FRAC_LAKE'     , 'LAKE fraction'          , '1', 'XY', landuse_out_dtype, vid(2) )
       call file_cartesc_def_var( fid, 'FRAC_URBAN'    , 'URBAN fraction'         , '1', 'XY', landuse_out_dtype, vid(3) )
       call file_cartesc_def_var( fid, 'FRAC_OCEAN_abs', 'absolute OCEAN fraction', '1', 'XY', landuse_out_dtype, vid(4) )
       call file_cartesc_def_var( fid, 'FRAC_LAND_abs' , 'absolute LAND fraction' , '1', 'XY', landuse_out_dtype, vid(5) )
       call file_cartesc_def_var( fid, 'FRAC_URBAN_abs', 'absolute URBAN fraction', '1', 'XY', landuse_out_dtype, vid(6) )
       call file_cartesc_def_var( fid, 'FRAC_LAKE_abs' , 'absolute LAKE fraction' , '1', 'XY', landuse_out_dtype, vid(7) )

       do p = 1, landuse_pft_mosaic
          write(varname,'(A8,I1.1)') 'FRAC_PFT', p
          call file_cartesc_def_var( fid, varname, 'PFT fraction', '1', 'XY', landuse_out_dtype, vid(8+2*(p-1)) )
          write(varname,'(A9,I1.1)') 'INDEX_PFT', p
          call file_cartesc_def_var( fid, varname, 'PFT index',    '1', 'XY', landuse_out_dtype, vid(9+2*(p-1)) )
       end do

       call file_cartesc_enddef( fid )

       call file_cartesc_write_var( fid, vid(1), landuse_frac_land(:,:), 'FRAC_LAND'     , 'XY' )
       call file_cartesc_write_var( fid, vid(2), landuse_frac_lake(:,:), 'FRAC_LAKE'     , 'XY' )
       call file_cartesc_write_var( fid, vid(3), landuse_frac_urban(:,:), 'FRAC_URBAN'    , 'XY' )
       call file_cartesc_write_var( fid, vid(4), landuse_fact_ocean(:,:), 'FRAC_OCEAN_abs', 'XY' )
       call file_cartesc_write_var( fid, vid(5), landuse_fact_land(:,:), 'FRAC_LAND_abs' , 'XY' )
       call file_cartesc_write_var( fid, vid(6), landuse_fact_urban(:,:), 'FRAC_URBAN_abs', 'XY' )
       call file_cartesc_write_var( fid, vid(7), landuse_fact_lake(:,:), 'FRAC_LAKE_abs' , 'XY' )

       do p = 1, landuse_pft_mosaic
          write(varname,'(A8,I1.1)') 'FRAC_PFT', p
          call file_cartesc_write_var( fid, vid(8+2*(p-1)), landuse_frac_pft(:,:,p), varname, 'XY' )
          write(varname,'(A9,I1.1)') 'INDEX_PFT', p
          temp(:,:) = real(LANDUSE_index_PFT(:,:,p),kind=rp)
          call file_cartesc_write_var( fid, vid(9+2*(p-1)), temp(:,:),               varname, 'XY' )
       end do

       call file_cartesc_close( fid )

    end if

    return
  end subroutine landuse_write

end module scale_landuse