mod_copytopo.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module mod_copytopo
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_atmos_grid_cartesc_index
  use scale_tracer
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: copytopo
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: copytopo_transgrid
  private :: copytopo_setalpha
  private :: copytopo_input_data
  private :: copytopo_mix_data
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer, private, parameter :: handle = 1
 
  character(len=H_SHORT), private :: COPYTOPO_IN_FILETYPE   = ''
  character(len=H_LONG),  private :: COPYTOPO_IN_BASENAME   = ''
  character(len=H_LONG),  private :: COPYTOPO_IN_POSTFIX    = ''
  character(len=H_SHORT), private :: COPYTOPO_IN_VARNAME    = 'topo'
  character(len=H_SHORT), private :: COPYTOPO_IN_LONNAME    = 'lon'
  character(len=H_SHORT), private :: COPYTOPO_IN_LATNAME    = 'lat'
  real(RP),               private :: COPYTOPO_TRANSITION_DX = -1.0_rp  
  real(RP),               private :: COPYTOPO_TRANSITION_DY = -1.0_rp  
  real(RP),               private :: COPYTOPO_FRACX         =  1.0_rp  
  real(RP),               private :: COPYTOPO_FRACY         =  1.0_rp  
  real(RP),               private :: COPYTOPO_taux          =  1.0_rp  
  real(RP),               private :: COPYTOPO_tauy          =  1.0_rp  
  logical,                private :: COPYTOPO_ENTIRE_REGION = .false.  
  logical,                private :: COPYTOPO_LINEAR_H      = .true.   
  real(RP),               private :: COPYTOPO_EXP_H         =  2.0_rp  
  character(len=H_SHORT), private :: COPYTOPO_INTRP_TYPE    = "LINEAR"
  integer,                private :: COPYTOPO_hypdiff_ORDER  = 4
  integer,                private :: COPYTOPO_hypdiff_NITER  = 20
 
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine copytopo( &
       TOPO_child )
    use scale_prc, only: &
       prc_abort
    implicit none
 
    real(rp), intent(inout) :: topo_child(:,:)
 
    real(rp) :: ctrx       (ia)
    real(rp) :: ctry       (ja)
    real(rp) :: alpha      (ia,ja)
    real(rp) :: topo_parent(ia,ja)
 
    namelist / param_copytopo / &
       copytopo_in_filetype,    &
       copytopo_in_basename,    &
       copytopo_in_postfix,     &
       copytopo_in_varname,     &
       copytopo_in_lonname,     &
       copytopo_in_latname,     &
       copytopo_transition_dx,  &
       copytopo_transition_dy,  &
       copytopo_fracx,          &
       copytopo_fracy,          &
       copytopo_taux,           &
       copytopo_tauy,           &
       copytopo_entire_region,  &
       copytopo_linear_h,       &
       copytopo_exp_h,          &
       copytopo_hypdiff_order,   &
       copytopo_hypdiff_niter
 
    integer :: ierr
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("COPYTOPO",*) 'Setup'
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_copytopo,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("COPYTOPO",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("COPYTOPO",*) 'Not appropriate names in namelist PARAM_COPYTOPO. Check!'
       call prc_abort
    endif
    log_nml(param_copytopo)
 
    ! copy topography from parent domain to transition region
 
    call copytopo_transgrid( ctrx(:), ctry(:) ) ! [OUT]
 
    call copytopo_setalpha ( ctrx(:), ctry(:), & ! [IN]
                             alpha(:,:)        ) ! [OUT]
 
    call copytopo_input_data( topo_parent(:,:) ) ! [OUT]
 
    call copytopo_mix_data( topo_parent(:,:), & ! [IN]
                            alpha(:,:), & ! [IN]
                            topo_child(:,:)  ) ! [INOUT]
 
    return
  end subroutine copytopo
 
  !-----------------------------------------------------------------------------
  subroutine copytopo_transgrid( &
       CTRX, &
       CTRY  )
    use scale_prc, only: &
       prc_abort, &
       prc_myrank
    use scale_prc_cartesc, only: &
       prc_2drank
    use scale_const, only: &
       eps => const_eps
    use scale_atmos_grid_cartesc, only: &
       cxg   => atmos_grid_cartesc_cxg,   &
       fxg   => atmos_grid_cartesc_fxg,   &
       cyg   => atmos_grid_cartesc_cyg,   &
       fyg   => atmos_grid_cartesc_fyg,   &
       cbfxg => atmos_grid_cartesc_cbfxg, &
       cbfyg => atmos_grid_cartesc_cbfyg
    implicit none
 
    real(rp), intent(out) :: ctrx(ia)
    real(rp), intent(out) :: ctry(ja)
 
    real(rp) :: ctrxg (  iag)
    real(rp) :: ctryg (  jag)
 
    real(rp) :: bufftotx, transtotx
    real(rp) :: bufftoty, transtoty
    integer  :: imain, ibuff, itrans
    integer  :: jmain, jbuff, jtrans
    integer  :: copy_is, copy_ie
    integer  :: copy_js, copy_je
 
    integer  :: i, j, ii, jj
    !---------------------------------------------------------------------------
 
    ! X-direction
    ! calculate buffer grid size
 
    do i = ihalo+1, iag
       if( abs(cbfxg(i)) < eps ) exit
    enddo
    ibuff = i - 1 - ihalo
    if ( ibuff == 0 ) then
       bufftotx = 0.0_rp
    else
       bufftotx = fxg(ibuff+ihalo) - fxg(ihalo)
    end if
 
    if ( copytopo_transition_dx < 0.0_rp ) then
       copytopo_transition_dx = bufftotx
    end if
    do i = ibuff+ihalo+1, iag
       if( cxg(i) - bufftotx - fxg(ihalo) >= copytopo_transition_dx ) exit
    enddo
    itrans = i - 1 - ihalo - ibuff
    if ( itrans == 0 ) then
       transtotx = 0.0_rp
    else
       transtotx = fxg(itrans+ibuff+ihalo) - fxg(ihalo) - bufftotx
    end if
 
    imain  = iag - 2*ibuff - 2*itrans - 2*ihalo
 
    if ( imain < 1 ) then
       log_error("COPYTOPO_transgrid",*) 'Not appropriate transition width for global domain(X).', copytopo_transition_dx
       log_error_cont(*) '# of buffer   region (one side) = ', ibuff
       log_error_cont(*) '# of transion region (one side) = ', itrans
       call prc_abort
    endif
 
    ! calc transition factor (global domaim)
    ctrxg(:) = 0.0_rp
 
    copy_is = 1
    copy_ie = ihalo+ibuff
    !$omp parallel do
    do i = copy_is, copy_ie
       ctrxg(i) = 1.0_rp
    enddo
 
    if ( itrans > 0 ) then
       copy_is = ihalo+ibuff+1
       copy_ie = ihalo+ibuff+itrans
       !$omp parallel do
       do i = copy_is, copy_ie
          ctrxg(i) = (transtotx+bufftotx+fxg(ihalo    )-cxg(i)) / transtotx
       enddo
       copy_is = iag - ihalo - ibuff - itrans + 1
       copy_ie = iag - ihalo - ibuff
       !$omp parallel do
       do i = copy_is, copy_ie
          ctrxg(i) = (transtotx+bufftotx-fxg(iag-ihalo)+cxg(i)) / transtotx
       enddo
    endif
 
    copy_is = iag - ihalo - ibuff + 1
    copy_ie = iag
    !$omp parallel do
    do i = copy_is, copy_ie
       ctrxg(i) = 1.0_rp
    enddo
 
    ctrxg(:) = max( min( ctrxg(:), 1.0_rp ), 0.0_rp )
 
    ! Y-direction
    ! calculate buffer grid size
    do j = jhalo+1, jag
       if( abs(cbfyg(j)) < eps ) exit
    enddo
    jbuff = j - 1 - jhalo
    if ( jbuff == 0 ) then
       bufftoty = 0.0_rp
    else
       bufftoty = fyg(jbuff+jhalo) - fyg(jhalo)
    end if
 
    if ( copytopo_transition_dy < 0.0_rp ) then
       copytopo_transition_dy = bufftoty
    end if
    do j = jbuff+jhalo+1, jag
       if( cyg(j) - bufftoty - fyg(jhalo) >= copytopo_transition_dy ) exit
    enddo
    jtrans = j - 1 - jhalo - jbuff
    if ( jtrans == 0 ) then
       transtoty = 0.0_rp
    else
       transtoty = fyg(jtrans+jbuff+jhalo) - fyg(jhalo) - bufftoty
    end if
 
    jmain  = jag - 2*jbuff - 2*jtrans - 2*jhalo
 
    if ( jmain < 1 ) then
       log_error("COPYTOPO_transgrid",*) 'Not appropriate transition width for global domain(Y).', copytopo_transition_dy
       log_error_cont(*) '# of buffer   region (one side)', jbuff
       log_error_cont(*) '# of transion region (one side)', jtrans
       call prc_abort
    endif
 
    ! calc transition factor (global domaim)
    ctryg(:) = 0.0_rp
 
    copy_js = 1
    copy_je = jhalo+jbuff
    !$omp parallel do
    do j = copy_js, copy_je
       ctryg(j) = 1.0_rp
    enddo
 
    if ( jtrans > 0 ) then
       copy_js = jhalo+jbuff+1
       copy_je = jhalo+jbuff+jtrans
       !$omp parallel do
       do j = copy_js, copy_je
          ctryg(j) = (transtoty+bufftoty+fyg(jhalo    )-cyg(j)) / transtoty
       enddo
       copy_js = jag - jhalo - jbuff - jtrans + 1
       copy_je = jag - jhalo - jbuff
       !$omp parallel do
       do j = copy_js, copy_je
          ctryg(j) = (transtoty+bufftoty-fyg(jag-jhalo)+cyg(j)) / transtoty
       enddo
    endif
 
    copy_js = jag - jhalo - jbuff + 1
    copy_je = jag
    !$omp parallel do
    do j = copy_js, copy_je
       ctryg(j) = 1.0_rp
    enddo
 
    ctryg(:) = max( min( ctryg(:), 1.0_rp ), 0.0_rp )
 
 
 
    ! horizontal coordinate (local domaim)
    !$omp parallel do &
    !$omp private(ii)
    do i = 1, ia
       ii = i + prc_2drank(prc_myrank,1) * imax
       ctrx(i) = ctrxg(ii)
    enddo
 
    !$omp parallel do &
    !$omp private(jj)
    do j = 1, ja
       jj = j + prc_2drank(prc_myrank,2) * jmax
       ctry(j) = ctryg(jj)
    enddo
 
    log_info("COPYTOPO_transgrid",*) 'COPYTOPO grid Information'
    log_info_cont(*) 'size of buffer   region (x and y; one side) = ', bufftotx, bufftoty
    log_info_cont(*) '#    of buffer   region (x and y; one side) = ', ibuff, jbuff
    log_info_cont(*) 'size of transion region (x and y; one side) = ', transtotx, transtoty
    log_info_cont(*) '#    of transion region (x and y; one side) = ', itrans, jtrans
 
 
    return
  end subroutine copytopo_transgrid
 
 
  !-----------------------------------------------------------------------------
  subroutine copytopo_setalpha( &
       CTRX, &
       CTRY, &
       alpha )
    use scale_const, only: &
       eps => const_eps
    use scale_comm_cartesc, only: &
       comm_vars8, &
       comm_wait
    implicit none
 
    real(rp), intent(in)  :: ctrx (ia)
    real(rp), intent(in)  :: ctry (ja)
    real(rp), intent(out) :: alpha(ia,ja)
 
    real(rp) :: coef_x, alpha_x1
    real(rp) :: coef_y, alpha_y1
    real(rp) :: ee1
 
    integer  :: i, j
    !---------------------------------------------------------------------------
 
    ! check invalid fraction
    copytopo_fracx = max( min( copytopo_fracx, 1.0_rp ), eps )
    copytopo_fracy = max( min( copytopo_fracy, 1.0_rp ), eps )
 
    if ( copytopo_taux <= 0.0_rp ) then ! invalid tau
       coef_x = 0.0_rp
    else
       coef_x = 1.0_rp / copytopo_taux
    endif
 
    if ( copytopo_tauy <= 0.0_rp ) then ! invalid tau
       coef_y = 0.0_rp
    else
       coef_y = 1.0_rp / copytopo_tauy
    endif
 
    !$omp parallel do &
    !$omp private(ee1,alpha_x1,alpha_y1)
    do j = 1, ja
    do i = 1, ia
       ee1 = ctrx(i)
 
       if ( ee1 <= 1.0_rp - copytopo_fracx ) then
          ee1 = 0.0_rp
       else
          ee1 = ( ee1 - 1.0_rp + copytopo_fracx ) / copytopo_fracx
       endif
 
       if ( copytopo_linear_h ) then
          alpha_x1 = coef_x * ee1
       else
          alpha_x1 = coef_x * ee1 * exp( -(1.0_rp-ee1) * copytopo_exp_h )
       endif
 
       ee1 = ctry(j)
 
       if ( ee1 <= 1.0_rp - copytopo_fracy ) then
          ee1 = 0.0_rp
       else
          ee1 = ( ee1 - 1.0_rp + copytopo_fracy ) / copytopo_fracy
       endif
 
       if ( copytopo_linear_h ) then
          alpha_y1 = coef_y * ee1
       else
          alpha_y1 = coef_y * ee1 * exp( -(1.0_rp-ee1) * copytopo_exp_h )
       endif
 
       alpha(i,j) = max( alpha_x1, alpha_y1 )
    enddo
    enddo
 
    call comm_vars8( alpha(:,:), 1 )
    call comm_wait ( alpha(:,:), 1 )
 
    return
  end subroutine copytopo_setalpha
 
 
  !-----------------------------------------------------------------------------
  subroutine copytopo_input_data( &
       TOPO_parent )
    use scale_prc, only: &
       prc_abort, &
       prc_ismaster
    use scale_const, only: &
       eps => const_eps, &
       pi  => const_pi
    use scale_interp, only: &
       interp_domain_compatibility, &
       interp_factor2d,             &
       interp_interp2d
    use scale_comm_cartesc, only: &
       comm_bcast, &
       comm_vars8, &
       comm_wait
    use scale_atmos_grid_cartesc_real, only: &
       lat => atmos_grid_cartesc_real_lat, &
       lon => atmos_grid_cartesc_real_lon
    use scale_atmos_grid_cartesc, only: &
       cx => atmos_grid_cartesc_cx, &
       cy => atmos_grid_cartesc_cy
    use scale_mapprojection, only: &
       mapprojection_lonlat2xy
    use scale_filter, only: &
       filter_hyperdiff
    use scale_landuse, only: &
       landuse_fact_ocean
    use scale_comm_cartesc_nest, only: &
       nest_interp_level => comm_cartesc_nest_interp_level
    implicit none
 
    real(rp), intent(out) :: topo_parent(:,:)
 
    integer :: ia_org, ja_org
 
    real(rp), allocatable :: lon_org (:,:), lat_org(:,:)
    real(rp), allocatable :: x_org   (:,:), y_org  (:,:)
    real(rp), allocatable :: topo_org(:,:)
 
    logical :: single
 
    ! for interpolation
    real(rp) :: dummy(1,1)
    integer  :: idx_i(ia,ja,nest_interp_level)
    integer  :: idx_j(ia,ja,nest_interp_level)
    real(rp) :: hfact(ia,ja,nest_interp_level)
    integer  :: interp_level
    logical  :: zonal, pole
 
    real(rp) :: topo_sign(ia,ja)
    real(rp) :: ocean_flag
 
    integer :: i, j
    !---------------------------------------------------------------------------
 
    select case ( copytopo_in_filetype )
    case ( 'SCALE' )
       call copytopo_get_size_scale( ia_org, ja_org ) ! (out)
       single = .false.
    case ( 'GrADS' )
       if ( prc_ismaster ) call copytopo_get_size_grads( ia_org, ja_org ) ! (out)
       single = .true.
    case ( 'WRF-ARW' )
       if ( prc_ismaster ) call copytopo_get_size_wrfarw( ia_org, ja_org ) ! (out)
       single = .true.
    case default
       log_error("COPYTOPO_input_data",*) 'COPYTOPO_IN_FILETYPE must be "SCALE" or "GrADS"'
       call prc_abort
    end select
 
    if ( single ) then
       call comm_bcast( ia_org )
       call comm_bcast( ja_org )
    end if
 
    allocate( lon_org(ia_org,ja_org) )
    allocate( lat_org(ia_org,ja_org) )
    allocate( topo_org(ia_org,ja_org) )
 
    select case ( copytopo_in_filetype )
    case ( 'SCALE' )
       call copytopo_get_data_scale( ia_org, ja_org,                           & ! (in)
                                     lon_org(:,:), lat_org(:,:), topo_org(:,:) ) ! (out)
    case ( 'GrADS' )
       if ( prc_ismaster ) &
       call copytopo_get_data_grads( ia_org, ja_org,                           & ! (in)
                                     lon_org(:,:), lat_org(:,:), topo_org(:,:) ) ! (out)
    case ( 'WRF-ARW' )
       if ( prc_ismaster ) &
       call copytopo_get_data_wrfarw( ia_org, ja_org,                           & ! (in)
                                      lon_org(:,:), lat_org(:,:), topo_org(:,:) ) ! (out)
    end select
 
    if ( single ) then
       call comm_bcast( lon_org(:,:), ia_org, ja_org )
       call comm_bcast( lat_org(:,:), ia_org, ja_org )
       call comm_bcast( topo_org(:,:), ia_org, ja_org )
    end if
 
    call interp_domain_compatibility( lon_org(:,:), lat_org(:,:), dummy(:,:),             &
                                      lon(:,:),     lat(:,:),     dummy(:,:), dummy(:,:), &
                                      skip_z=.true.                                       )
 
    select case ( copytopo_intrp_type )
    case ( "LINEAR" )
 
       allocate( x_org(ia_org,ja_org) )
       allocate( y_org(ia_org,ja_org) )
       call mapprojection_lonlat2xy( ia_org, 1, ia_org, ja_org, 1, ja_org, &
                                     lon_org(:,:), lat_org(:,:),           & ! [IN]
                                     x_org(:,:), y_org(:,:)                ) ! [OUT]
 
       zonal = ( maxval(lon_org) - minval(lon_org) ) > 2.0_rp * pi * 0.9_rp
       pole =  ( maxval(lat_org) > pi * 0.5_rp * 0.9_rp ) .or. ( minval(lat_org) < - pi * 0.5_rp * 0.9_rp )
       call interp_factor2d( ia_org, ja_org, ia, ja,     & ! [IN]
                             x_org(:,:), y_org(:,:),     & ! [IN]
                             cx(:), cy(:),               & ! [IN]
                             idx_i(:,:,:), idx_j(:,:,:), & ! [OUT]
                             hfact(:,:,:),               & ! [OUT]
                             zonal = zonal, pole = pole  ) ! [IN]
       deallocate( x_org, y_org )
       interp_level = 4
 
    case ( "DIST-WEIGHT" )
 
       call interp_factor2d( nest_interp_level,          & ! [IN]
                             ia_org, ja_org, ia, ja,     & ! [IN]
                             lon_org(:,:), lat_org(:,:), & ! [IN]
                             lon(:,:), lat(:,:),         & ! [IN]
                             idx_i(:,:,:), idx_j(:,:,:), & ! [OUT]
                             hfact(:,:,:)                ) ! [OUT]
       interp_level = nest_interp_level
 
       case default
          log_error("COPYTOPO_input_data",*) "Unsupported type of COPYTOPO_INTRP_TYPE : ", trim(copytopo_intrp_type)
          log_error_cont(*)                  'It must be "LINEAR" or "DIST-WEIGHT"'
          call prc_abort
    end select
 
    call interp_interp2d( interp_level,       & ! [IN]
                          ia_org, ja_org,     & ! [IN]
                          ia, ja,             & ! [IN]
                          idx_i(:,:,:), & ! [IN]
                          idx_j(:,:,:), & ! [IN]
                          hfact(:,:,:), & ! [IN]
                          topo_org(:,:),   & ! [IN]
                          topo_parent(:,:)    ) ! [OUT]
 
    deallocate( topo_org, lon_org, lat_org )
 
    if ( copytopo_hypdiff_niter > 1 ) then
       !$omp parallel do &
       !$omp private(ocean_flag)
       do j = 1, ja
       do i = 1, ia
          ocean_flag = ( 0.5_rp + sign( 0.5_rp, landuse_fact_ocean(i,j) - 1.0_rp + eps ) ) & ! fact_ocean==1
                     * ( 0.5_rp + sign( 0.5_rp, eps - abs(topo_parent(i,j)) ) )              ! |TOPO| > EPS
          topo_sign(i,j) = sign( 1.0_rp, topo_parent(i,j) ) * ( 1.0_rp - ocean_flag )
       end do
       end do
 
       call filter_hyperdiff( ia, isb, ieb, ja, jsb, jeb, &
                              topo_parent(:,:), copytopo_hypdiff_order, copytopo_hypdiff_niter, &
                              limiter_sign = topo_sign(:,:) )
    end if
 
 
    call comm_vars8( topo_parent(:,:), 1 )
    call comm_wait ( topo_parent(:,:), 1 )
 
 
    return
  end subroutine copytopo_input_data
 
  !-----------------------------------------------------------------------------
  subroutine copytopo_mix_data( &
       TOPO_parent, &
       alpha,       &
       TOPO_child   )
    implicit none
 
    real(rp), intent(in)    :: topo_parent(:,:) ! topography of parent domain
    real(rp), intent(in)    :: alpha      (:,:) ! dumping coefficient (0-1)
    real(rp), intent(inout) :: topo_child (:,:) ! topography of child  domain
 
    integer  :: i, j
    !---------------------------------------------------------------------------
 
    if ( copytopo_entire_region ) then
       !$omp parallel do
       do j = 1, ja
       do i = 1, ia
          topo_child(i,j) = topo_parent(i,j)
       enddo
       enddo
    else
       !$omp parallel do
       do j = 1, ja
       do i = 1, ia
          topo_child(i,j) = ( 1.0_rp-alpha(i,j) ) * topo_child(i,j) &
                          + (        alpha(i,j) ) * topo_parent(i,j)
       enddo
       enddo
    endif
 
    return
  end subroutine copytopo_mix_data
 
  !-----------------------------------------------------------------------------
  subroutine copytopo_get_size_scale( &
       IA_org, JA_org )
    use scale_comm_cartesc_nest, only: &
       nest_tile_num_x => comm_cartesc_nest_tile_num_x, &
       nest_tile_num_y => comm_cartesc_nest_tile_num_y, &
       parent_imax,       &
       parent_jmax
    implicit none
    integer, intent(out) :: ia_org
    integer, intent(out) :: ja_org
 
    ia_org = parent_imax(handle) * nest_tile_num_x
    ja_org = parent_jmax(handle) * nest_tile_num_y
 
    return
  end subroutine copytopo_get_size_scale
 
  !-----------------------------------------------------------------------------
  subroutine copytopo_get_data_scale( &
       IA_org, JA_org,   &
       LON_org, LAT_org, &
       TOPO_org          )
    use scale_const, only: &
       d2r => const_d2r
    use scale_comm_cartesc_nest, only: &
       nest_domain_shape => comm_cartesc_nest_domain_shape, &
       nest_tile_id      => comm_cartesc_nest_tile_id
    use scale_file, only: &
       file_open, &
       file_read, &
       file_close
    implicit none
    integer,  intent(in)  :: IA_org, JA_org
    real(RP), intent(out) :: LON_org (IA_org,JA_org)
    real(RP), intent(out) :: LAT_org (IA_org,JA_org)
    real(RP), intent(out) :: TOPO_org(IA_org,JA_org)
 
    real(RP), allocatable :: read2D(:,:)
 
    integer :: tilei, tilej
    integer :: cxs, cxe, cys, cye ! for child domain
    integer :: pxs, pxe, pys, pye ! for parent domain
    integer :: fid
 
    integer :: rank
    integer :: n
 
    do n = 1, size( nest_tile_id(:) )
       ! read data from split files
       rank = nest_tile_id(n)
 
       call nest_domain_shape( tilei, tilej, & ! [OUT]
                               cxs,   cxe,   & ! [OUT]
                               cys,   cye,   & ! [OUT]
                               pxs,   pxe,   & ! [OUT]
                               pys,   pye,   & ! [OUT]
                               n             ) ! [IN]
 
       allocate( read2d(tilei,tilej) )
 
       call file_open( copytopo_in_basename,          & ! [IN]
                       fid,                           & ! [OUT]
                       aggregate=.false., rankid=rank ) ! [IN]
 
       call file_read( fid, "lon",  read2d(:,:) )
       lon_org(cxs:cxe,cys:cye) = read2d(pxs:pxe,pys:pye) * d2r
       call file_read( fid, "lat",  read2d(:,:) )
       lat_org(cxs:cxe,cys:cye) = read2d(pxs:pxe,pys:pye) * d2r
       call file_read( fid, copytopo_in_varname, read2d(:,:) )
       topo_org(cxs:cxe,cys:cye) = read2d(pxs:pxe,pys:pye)
       deallocate( read2d )
 
       call file_close( fid )
 
    enddo
 
    return
  end subroutine copytopo_get_data_scale
 
  !-----------------------------------------------------------------------------
  subroutine copytopo_get_size_grads( &
       IA_org, JA_org )
    use scale_file_grads, only: &
       file_grads_open, &
       file_grads_get_shape
    implicit none
    integer, intent(out) :: IA_org, JA_org
 
    integer :: file_id
    integer :: shape(2)
 
    call file_grads_open( copytopo_in_basename, & ! (in)
                          file_id               ) ! (out)
 
    call file_grads_get_shape( file_id, copytopo_in_varname, & ! (in)
                               shape(:)                      ) ! (out)
 
    ia_org = shape(1)
    ja_org = shape(2)
 
    return
  end subroutine copytopo_get_size_grads
 
  !-----------------------------------------------------------------------------
  subroutine copytopo_get_data_grads( &
       IA_org, JA_org,   &
       LON_org, LAT_org, &
       TOPO_org          )
    use scale_const, only: &
       d2r => const_d2r
    use scale_file_grads, only: &
       file_grads_open,   &
       file_grads_varid,  &
       file_grads_isoned, &
       file_grads_read,   &
       file_grads_close
    implicit none
    integer,  intent(in)  :: IA_org, JA_org
    real(RP), intent(out) :: LON_org (IA_org,JA_org)
    real(RP), intent(out) :: LAT_org (IA_org,JA_org)
    real(RP), intent(out) :: TOPO_org(IA_org,JA_org)
 
    integer :: file_id, var_id
    real(RP) :: lon1D(IA_org), lat1D(JA_org)
 
    integer :: i, j
 
    call file_grads_open( copytopo_in_basename, & ! (in)
                          file_id               ) ! (out)
 
    ! lon
    call file_grads_varid( file_id, copytopo_in_lonname, & ! (in)
                           var_id                        ) ! (out)
    if ( file_grads_isoned( file_id, var_id ) ) then
       call file_grads_read( file_id, var_id, & ! (in)
                             lon1d(:)         ) ! (out)
       !$omp parallel do collapse(2)
       do j = 1, ja_org
       do i = 1, ia_org
          lon_org(i,j) = lon1d(i) * d2r
       end do
       end do
    else
       call file_grads_read( file_id, var_id,              & ! (in)
                             lon_org(:,:),                 & ! (out)
                             postfix = copytopo_in_postfix ) ! (in)
       !$omp parallel do collapse(2)
       do j = 1, ja_org
       do i = 1, ia_org
          lon_org(i,j) = lon_org(i,j) * d2r
       end do
       end do
    end if
 
    ! lat
    call file_grads_varid( file_id, copytopo_in_latname, & ! (in)
                           var_id                        ) ! (out)
    if ( file_grads_isoned( file_id, var_id ) ) then
       call file_grads_read( file_id, var_id, & ! (in)
                             lat1d(:)         ) ! (out)
       !$omp parallel do collapse(2)
       do j = 1, ja_org
       do i = 1, ia_org
          lat_org(i,j) = lat1d(j) * d2r
       end do
       end do
    else
       call file_grads_read( file_id, var_id,              & ! (in)
                             lat_org(:,:),                 & ! (out)
                             postfix = copytopo_in_postfix ) ! (in)
       !$omp parallel do collapse(2)
       do j = 1, ja_org
       do i = 1, ia_org
          lat_org(i,j) = lat_org(i,j) * d2r
       end do
       end do
    end if
 
    ! topo
    call file_grads_varid( file_id, copytopo_in_varname, & ! (in)
                           var_id                        ) ! (out)
    call file_grads_read( file_id, var_id,              & ! (in)
                          topo_org(:,:),                & ! (out)
                          postfix = copytopo_in_postfix ) ! (in)
 
    ! close
    call file_grads_close( file_id )
 
    return
  end subroutine copytopo_get_data_grads
 
  !-----------------------------------------------------------------------------
  subroutine copytopo_get_size_wrfarw( &
       IA_org, JA_org )
    use scale_file_grads, only: &
       file_grads_open, &
       file_grads_get_shape
    use scale_file, only: &
       file_open, &
       file_get_dimlength
    implicit none
    integer, intent(out) :: IA_org, JA_org
 
    integer :: fid
 
    call file_open( copytopo_in_basename,          & ! (in)
                    fid,                           & ! (out)
                    postfix = copytopo_in_postfix, & ! (in)
                    single = .true.                ) ! (in)
 
    call file_get_dimlength( fid, "west_east",   ia_org )
    call file_get_dimlength( fid, "south_north", ja_org )
 
    return
  end subroutine copytopo_get_size_wrfarw
 
  subroutine copytopo_get_data_wrfarw( &
       IA_org, JA_org,   &
       LON_org, LAT_org, &
       TOPO_org          )
    use scale_const, only: &
       d2r => const_d2r
    use scale_file, only: &
       file_open, &
       file_read, &
       file_close
    implicit none
    integer,  intent(in)  :: IA_org, JA_org
    real(RP), intent(out) :: LON_org (IA_org,JA_org)
    real(RP), intent(out) :: LAT_org (IA_org,JA_org)
    real(RP), intent(out) :: TOPO_org(IA_org,JA_org)
 
    integer :: fid
    integer :: i, j
 
    call file_open( copytopo_in_basename,          & ! (in)
                    fid,                           & ! (out)
                    postfix = copytopo_in_postfix, & ! (in)
                    single = .true.                ) ! (in)
 
    call file_read( fid, "XLAT", lat_org(:,:) )
 
    call file_read( fid, "XLONG", lon_org(:,:) )
 
    !$omp parallel do collapse(2)
    do j = 1, ja_org
    do i = 1, ia_org
       lat_org(i,j) = lat_org(i,j) * d2r
       lon_org(i,j) = lon_org(i,j) * d2r
    end do
    end do
 
    call file_read( fid, "HGT", topo_org(:,:) )
 
    call file_close( fid )
 
  end subroutine copytopo_get_data_wrfarw
 
end module mod_copytopo