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_LONG), private :: copytopo_in_basename   = ''
  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  
  integer,               private :: copytopo_filter_order  = 2
  integer,               private :: copytopo_filter_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_basename,    &
       copytopo_transition_dx,  &
       copytopo_transition_dy,  &
       copytopo_fracx,          &
       copytopo_fracy,          &
       copytopo_taux,           &
       copytopo_tauy,           &
       copytopo_entire_region,  &
       copytopo_linear_h,       &
       copytopo_exp_h,          &
       copytopo_filter_order,   &
       copytopo_filter_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_const, only: &
       eps => const_eps, &
       d2r => const_d2r
    use scale_file, only: &
       file_open, &
       file_read, &
       file_close
    use scale_interp, only: &
       interp_domain_compatibility, &
       interp_factor2d,             &
       interp_interp2d
    use scale_comm_cartesc, only: &
       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_comm_cartesc_nest, only: &
       nest_interp_level => comm_cartesc_nest_interp_level, &
       nest_domain_shape => comm_cartesc_nest_domain_shape, &
       nest_tile_num_x   => comm_cartesc_nest_tile_num_x,   &
       nest_tile_num_y   => comm_cartesc_nest_tile_num_y,   &
       nest_tile_id      => comm_cartesc_nest_tile_id,      &
       parent_imax,       &
       parent_jmax
    use scale_filter, only: &
       filter_hyperdiff
    use scale_landuse, only: &
       landuse_fact_ocean
    implicit none

    real(RP), intent(out) :: TOPO_parent(:,:)

    real(RP), allocatable :: LON_org (:,:)
    real(RP), allocatable :: LAT_org (:,:)
    real(RP), allocatable :: TOPO_org(:,:)
    real(RP), allocatable :: read2D  (:,:)

    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)

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

    integer :: IA_org, JA_org     ! number of grids for whole domain
    integer :: tilei, tilej
    integer :: cxs, cxe, cys, cye ! for child domain
    integer :: pxs, pxe, pys, pye ! for parent domain
    integer :: rank

    real(RP) :: ocean_flag

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

    ia_org = parent_imax(handle) * nest_tile_num_x
    ja_org = parent_jmax(handle) * nest_tile_num_y

    allocate( lon_org(ia_org,ja_org) )
    allocate( lat_org(ia_org,ja_org) )
    allocate( topo_org(ia_org,ja_org) )

    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, "TOPO", read2d(:,:) )
       topo_org(cxs:cxe,cys:cye) = read2d(pxs:pxe,pys:pye)
       deallocate( read2d )

       call file_close( fid )

    enddo

    call interp_domain_compatibility( lon_org(:,:), lat_org(:,:), dummy(:,:),             &
                                      lon(:,:),     lat(:,:),     dummy(:,:), dummy(:,:), &
                                      skip_z=.true.                                       )

    call interp_factor2d( nest_interp_level, & ! [IN]
                          ia_org, ja_org,    & ! [IN]
                          lon_org(:,:),      & ! [IN]
                          lat_org(:,:),      & ! [IN]
                          ia, ja,            & ! [IN]
                          lon(:,:),      & ! [IN]
                          lat(:,:),      & ! [IN]
                          idx_i(:,:,:),    & ! [OUT]
                          idx_j(:,:,:),    & ! [OUT]
                          hfact(:,:,:)     ) ! [OUT]

    call interp_interp2d( nest_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]

    if ( copytopo_filter_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, is, ie, ja, js, je, &
                              topo_parent(:,:), copytopo_filter_order, copytopo_filter_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

end module mod_copytopo