mod_copytopo.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
module mod_copytopo
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use gtool_file_h
  use gtool_file, only: &
     fileread
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_tracer
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: copytopo

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  integer, public            :: cnvtopo_type = -1

  integer, public, parameter :: i_ignore     =  0
  integer, public, parameter :: i_gtopo30    =  1
  integer, public, parameter :: i_dem50m     =  2
  integer, public, parameter :: i_gmted2010  =  3

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: copytopo_transgrid
  private :: copytopo_setalpha
  private :: copytopo_input_data
  private :: copytopo_mix_data

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer, private, parameter :: handle = 1
  integer, private            :: itp_nh = 4

  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_transfact     = -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   

  real(RP), private, allocatable :: ctrx(:)
  real(RP), private, allocatable :: ctry(:)
  real(RP), private, allocatable :: copytopo_alpha(:,:)
  real(RP), private, allocatable :: topo_pd(:,:)

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine copytopo( &
      topo_cd          ) ! [inout]
    use scale_process, only: &
       prc_mpistop
    use scale_grid, only: &
       dx,         &
       dy,         &
       buffer_dx,  &
       buffer_dy,  &
       bufffact
    use scale_grid_nest, only: &
       nest_interp_level
    implicit none

    real(RP), intent(inout) :: topo_cd(:,:)

    namelist / param_copytopo / &
       copytopo_in_basename,    &
       copytopo_transition_dx,  &
       copytopo_transition_dy,  &
       copytopo_transfact,      &
       copytopo_fracx,          &
       copytopo_fracy,          &
       copytopo_taux,           &
       copytopo_tauy,           &
       copytopo_entire_region,  &
       copytopo_linear_h,       &
       copytopo_exp_h

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

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '+++ Module[COPYTOPO]/Categ[COPYTOPO]'

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

    if ( copytopo_transition_dx < 0.0_rp ) then
       copytopo_transition_dx = buffer_dx
    end if
    if ( copytopo_transition_dy < 0.0_rp ) then
       copytopo_transition_dy = buffer_dy
    end if
    if ( copytopo_transfact < 0.0_rp ) then
       copytopo_transfact = bufffact
    end if

    allocate( ctrx(ia) )
    allocate( ctry(ja) )
    allocate( copytopo_alpha(ia,ja) )
    allocate( topo_pd(ia,ja) )
    copytopo_alpha(:,:) = 0.0_rp

    itp_nh = int( nest_interp_level )

    ! copy topography from parent domain to transition region
    call copytopo_transgrid

    call copytopo_setalpha

    call copytopo_input_data( topo_pd ) ! (out)

    call copytopo_mix_data( topo_cd, &  ! (inout)
                            topo_pd  )  ! (in)

    return
  end subroutine copytopo

  !-----------------------------------------------------------------------------
  !-----------------------------------------------------------------------------

  !-----------------------------------------------------------------------------
  subroutine copytopo_transgrid
    use scale_process, only: &
       prc_mpistop, &
       prc_myrank
    use scale_rm_process, only: &
       prc_2drank,  &
       prc_num_x,   &
       prc_num_y
    use scale_const, only: &
       radius => const_radius, &
       eps    => const_eps
    use scale_grid, only: &
       dx,                  &
       dy,                  &
       cxg   => grid_cxg,   &
       fxg   => grid_fxg,   &
       cyg   => grid_cyg,   &
       fyg   => grid_fyg,   &
       cbfxg => grid_cbfxg, &
       cbfyg => grid_cbfyg, &
       bufffact
    implicit none

    real(RP), allocatable :: ctrxg(:)
    real(RP), allocatable :: ctryg(:)

    real(RP), allocatable :: buffx(:),  buffy(:)
    real(RP)              :: bufftotx,  bufftoty
    real(RP), allocatable :: transx(:), transy(:)
    real(RP)              :: transtotx, transtoty

    integer :: imain, ibuff, itrans
    integer :: jmain, jbuff, jtrans
    integer :: copy_is, copy_ie, copy_js, copy_je

    integer :: i, j, ii, jj
    !---------------------------------------------------------------------------

    allocate( buffx(0:iag) )
    allocate( buffy(0:jag) )
    allocate( transx(0:iag) )
    allocate( transy(0:jag) )
    allocate( ctrxg(  iag) )
    allocate( ctryg(  jag) )

    ! X-direction
    ! calculate buffer grid size
    buffx(:)  = dx
    transx(0) = dx
    bufftotx  = 0.0_rp
    transtotx = 0.0_rp

    do i = ihalo+1, iag
       if( abs(cbfxg(i) - 0.0_rp) < eps ) exit
       buffx(i) = buffx(i-1) * bufffact
       bufftotx = bufftotx + buffx(i)
    enddo
    ibuff = i - (ihalo+1)

    do i = 1, iag
       if( transtotx >= copytopo_transition_dx ) exit
       transx(i) = transx(i-1) * copytopo_transfact
       transtotx = transtotx + transx(i)
    enddo
    itrans = i - 1
    imain  = iag - 2*ibuff - 2*itrans - 2*ihalo

    if ( imain < 1 ) then
       write(*,*) 'xxx Not appropriate transition width for global domain(X).', copytopo_transition_dx
       write(*,*) '    # of buffer region (one side)', ibuff
       write(*,*) '    # of transion region (one side)', itrans
       call prc_mpistop
    endif

    ! calc transition factor (global domaim)
    ctrxg(:) = 0.0_rp
    do i = 1, ihalo+ibuff
       ctrxg(i) = 1.0_rp
    enddo

    if ( itrans > 0 ) then
       copy_is = ihalo+ibuff+1
       copy_ie = ihalo+ibuff+itrans
       do i = copy_is, copy_ie
          ctrxg(i) = (transtotx+bufftotx+fxg(ihalo    )-cxg(i)) / transtotx
       enddo
       copy_is = ihalo+ibuff+itrans+imain+1
       copy_ie = ihalo+ibuff+itrans+imain+itrans+ibuff
       do i = copy_is, copy_ie
          ctrxg(i) = (transtotx+bufftotx-fxg(iag-ihalo)+cxg(i)) / transtotx
       enddo
    endif

    copy_is = ihalo+ibuff+itrans+imain+itrans+ibuff+1
    copy_ie = ihalo+ibuff+itrans+imain+itrans+ibuff+ihalo
    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
    buffy(:)  = dy
    transy(0) = dy
    bufftoty  = 0.0_rp
    transtoty = 0.0_rp

    do j = jhalo+1, jag
       if( abs(cbfyg(j) - 0.0_rp) < eps ) exit
       buffy(j) = buffy(j-1) * bufffact
       bufftoty = bufftoty + buffy(j)
    enddo
    jbuff = j - (jhalo+1)

    do j = 1, jag
       if( transtoty >= copytopo_transition_dy ) exit
       transy(j) = transy(j-1) * copytopo_transfact
       transtoty = transtoty + transy(j)
    enddo
    jtrans = j - 1
    jmain  = jag - 2*jbuff - 2*jtrans - 2*jhalo

    if ( jmain < 1 ) then
       write(*,*) 'xxx Not appropriate transition width for global domain(Y).', copytopo_transition_dy
       write(*,*) '    # of buffer region (one side)', jbuff
       write(*,*) '    # of transion region (one side)', jtrans
       call prc_mpistop
    endif

    ! calc transition factor (global domaim)
    ctryg(:) = 0.0_rp
    do j = 1, jhalo+jbuff
       ctryg(j) = 1.0_rp
    enddo

    if ( jtrans > 0 ) then
       copy_js = jhalo+jbuff+1
       copy_je = jhalo+jbuff+jtrans
       do j = copy_js, copy_je
          ctryg(j) = (transtoty+bufftoty+fyg(jhalo    )-cyg(j)) / transtoty
       enddo
       copy_js = jhalo+jbuff+jtrans+jmain+1
       copy_je = jhalo+jbuff+jtrans+jmain+jtrans+jbuff
       do j = copy_js, copy_je
          ctryg(j) = (transtoty+bufftoty-fyg(jag-jhalo)+cyg(j)) / transtoty
       enddo
    endif

    copy_js = jhalo+jbuff+jtrans+jmain+jtrans+jbuff+1
    copy_je = jhalo+jbuff+jtrans+jmain+jtrans+jbuff+jhalo
    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)
    do i = 1, ia
       ii = i + prc_2drank(prc_myrank,1) * imax
       ctrx(i) = ctrxg(ii)
    enddo
    do j = 1, ja
       jj = j + prc_2drank(prc_myrank,2) * jmax
       ctry(j) = ctryg(jj)
    enddo

    deallocate( transx )
    deallocate( transy )

    return
  end subroutine copytopo_transgrid


  !-----------------------------------------------------------------------------
  subroutine copytopo_setalpha
    use scale_const, only: &
       eps => const_eps
    use scale_comm, only: &
       comm_vars8, &
       comm_wait

    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

    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 )
       end if

       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 )
       end if

       copytopo_alpha(i,j) = max( alpha_x1, alpha_y1 )
    enddo
    enddo

    call comm_vars8( copytopo_alpha(:,:), 1 )
    call comm_wait ( copytopo_alpha(:,:), 1 )

    return
  end subroutine copytopo_setalpha


  !-----------------------------------------------------------------------------
  subroutine copytopo_input_data( &
      topo_pd      ) ! (out)
    use scale_const, only: &
       d2r => const_d2r
    use scale_comm, only: &
       comm_vars8, &
       comm_wait
    use scale_grid_nest, only: &
       parent_imax,     &
       parent_jmax,     &
       nest_tile_num_x, &
       nest_tile_num_y, &
       nest_tile_id,    &
       nest_domain_shape
    use scale_interpolation_nest, only: &
       intrpnest_domain_compatibility,  &
       intrpnest_interp_fact_latlon,    &
       intrpnest_interp_2d
    use scale_grid_real, only: &
       lat => real_lat, &
       lon => real_lon
    implicit none

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

    real(RP)              :: dummy (1,1,1)
    real(RP), allocatable :: read2d(:,:)
    real(RP), allocatable :: lon_org (:,:)
    real(RP), allocatable :: lat_org (:,:)
    real(RP), allocatable :: topo_org(:,:)
    real(RP), allocatable :: hfact(:,:,:)
    integer,  allocatable :: igrd (:,:,:)
    integer,  allocatable :: jgrd (:,:,:)

    integer :: iall, jall  ! number of grids for whole domain
    integer :: pti, ptj    ! number of grids for a tile
    integer :: tilei, tilej
    integer :: rank
    integer :: i
    integer :: cxs, cxe, cys, cye  ! for child domain
    integer :: pxs, pxe, pys, pye  ! for parent domain
    !---------------------------------------------------------------------------

    pti  = parent_imax(handle)
    ptj  = parent_jmax(handle)
    iall = pti * nest_tile_num_x
    jall = ptj * nest_tile_num_y

    allocate( hfact( ia,   ja,  itp_nh ) )
    allocate( igrd( ia,   ja,  itp_nh ) )
    allocate( jgrd( ia,   ja,  itp_nh ) )
    allocate( lon_org( iall, jall        ) )
    allocate( lat_org( iall, jall        ) )
    allocate( topo_org( iall, jall        ) )

    do i = 1, size( nest_tile_id(:) )
       ! read data from split files
       rank = nest_tile_id(i)

       call nest_domain_shape ( tilei, tilej,   & ! [out]
                                cxs,   cxe,     & ! [out]
                                cys,   cye,     & ! [out]
                                pxs,   pxe,     & ! [out]
                                pys,   pye,     & ! [out]
                                i               ) ! [in ]
       allocate( read2d( tilei,tilej ) )

       call fileread( read2d(:,:), copytopo_in_basename, "lon",  1, rank )
       lon_org(cxs:cxe,cys:cye) = read2d(pxs:pxe,pys:pye) * d2r
       call fileread( read2d(:,:), copytopo_in_basename, "lat",  1, rank )
       lat_org(cxs:cxe,cys:cye) = read2d(pxs:pxe,pys:pye) * d2r
       call fileread( read2d(:,:), copytopo_in_basename, "TOPO", 1, rank )
       topo_org(cxs:cxe,cys:cye) = read2d(pxs:pxe,pys:pye)

       deallocate( read2d )
    end do

    call intrpnest_domain_compatibility( lon_org(:,:), lat_org(:,:), dummy(:,:,:), &
                                         lon(:,:),     lat(:,:),     dummy(:,:,:), &
                                         skip_z=.true.                             )
    !call check_domain_compatibility( lon_org(:,:), lat_org(:,:), dummy(:,:,:), &
    !                                 LON(:,:),     LAT(:,:),     dummy(:,:,:), &
    !                                 skip_z=.true.                             )

    call intrpnest_interp_fact_latlon( hfact(:,:,:),   & ! (out)
                                       igrd(:,:,:),   & ! (out)
                                       jgrd(:,:,:),   & ! (out)
                                       lat(:,:),     & ! (in)
                                       lon(:,:),     & ! (in)
                                       ia,   ja,         & ! (in)
                                       lat_org(:,:),     & ! (in)
                                       lon_org(:,:),     & ! (in)
                                       iall, jall        ) ! (in)

    call intrpnest_interp_2d( topo_pd(:,:),    &
                              topo_org(:,:),    &
                              hfact(:,:,:),  &
                              igrd(:,:,:),  &
                              jgrd(:,:,:),  &
                              ia, ja            )

    call comm_vars8( topo_pd(:,:), 1 )
    call comm_wait ( topo_pd(:,:), 1 )

    return
  end subroutine copytopo_input_data

  !-----------------------------------------------------------------------------
  subroutine copytopo_mix_data( &
      topo_cd,  &  ! (inout)
      topo_pd    ) ! (in)
    implicit none
    real(RP), intent(inout) :: topo_cd(:,:)  ! topography of child domain (mine)
    real(RP), intent(in)    :: topo_pd(:,:)  ! topography of parent domain

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

    if ( copytopo_entire_region ) then
       topo_cd(:,:) = topo_pd(:,:)
    else
       do j = 1, ja
       do i = 1, ia
          frac = copytopo_alpha(i,j)
          topo_cd(i,j) = topo_cd(i,j) * ( 1.0_rp - frac ) &
                       + topo_pd(i,j) * frac
       end do
       end do
    endif

    return
  end subroutine copytopo_mix_data

end module mod_copytopo