scale_gridtrans.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
module scale_gridtrans
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: gtrans_setup

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  real(RP), public, allocatable :: gtrans_mapf (:,:,:,:)
  real(RP), public, allocatable :: gtrans_rotc (:,:,:)

  real(RP), public, allocatable :: gtrans_gsqrt(:,:,:,:)
  real(RP), public, allocatable :: gtrans_j13g (:,:,:,:)
  real(RP), public, allocatable :: gtrans_j23g (:,:,:,:)
  real(RP), public              :: gtrans_j33g

  real(RP), public, allocatable :: gtrans_limyz(:,:,:,:)
  real(RP), public, allocatable :: gtrans_limxz(:,:,:,:)
  real(RP), public, allocatable :: gtrans_limxy(:,:,:,:)

  integer,  public :: i_xyz = 1 ! at (x,y,z)
  integer,  public :: i_xyw = 2 ! at (x,y,w)
  integer,  public :: i_uyw = 3 ! at (u,y,w)
  integer,  public :: i_xvw = 4 ! at (x,v,w)
  integer,  public :: i_uyz = 5 ! at (u,y,z)
  integer,  public :: i_xvz = 6 ! at (x,v,z)
  integer,  public :: i_uvz = 7 ! at (u,v,z)

  integer,  public :: i_xy  = 1 ! at (x,y)
  integer,  public :: i_uy  = 2 ! at (u,y)
  integer,  public :: i_xv  = 3 ! at (x,v)
  integer,  public :: i_uv  = 4 ! at (u,v)

  integer,  public :: i_fyz = 1 ! y-z face limiting x-flux
  integer,  public :: i_fxz = 2 ! x-z face limiting y-flux
  integer,  public :: i_fxy = 3 ! x-y face limiting z-flux

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: gtrans_mapfactor
  private :: gtrans_terrainfollowing
  private :: gtrans_thin_wall
  private :: gtrans_step_mountain
  private :: gtrans_write

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  character(len=H_LONG),  private :: gtrans_out_basename  = ''
  character(len=H_MID),   private :: gtrans_out_title     = 'SCALE-RM GEOMETRICS'
  character(len=H_SHORT), private :: gtrans_out_dtype     = 'DEFAULT'

  character(len=H_SHORT), private :: gtrans_topo_type     = 'TERRAINFOLLOWING'
  integer,                private :: gtrans_thinwall_xdiv = 50
  integer,                private :: gtrans_thinwall_ydiv = 50
  logical,                private :: debug                = .false.
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine gtrans_setup
    use scale_process, only: &
       prc_mpistop
    implicit none

    namelist / param_gtrans / &
       gtrans_out_basename,  &
       gtrans_out_dtype,     &
       gtrans_topo_type,     &
       gtrans_thinwall_xdiv, &
       gtrans_thinwall_ydiv, &
       debug

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

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[GRIDTRANS] / Categ[ATMOS-RM GRID] / Origin[SCALElib]'

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

    allocate( gtrans_mapf(ia,ja,2,4) )

    allocate( gtrans_rotc(ia,ja,2) )

    allocate( gtrans_gsqrt(ka,ia,ja,7) )
    allocate( gtrans_j13g(ka,ia,ja,7) )
    allocate( gtrans_j23g(ka,ia,ja,7) )

    gtrans_gsqrt(:,:,:,:) = 1.0_rp
    gtrans_j13g(:,:,:,:) = 0.0_rp
    gtrans_j23g(:,:,:,:) = 0.0_rp
    gtrans_j33g = 1.0_rp

    allocate( gtrans_limyz(ka,ia,ja,7) )
    allocate( gtrans_limxz(ka,ia,ja,7) )
    allocate( gtrans_limxy(ka,ia,ja,7) )
    gtrans_limyz(:,:,:,:) = 1.0_rp
    gtrans_limxz(:,:,:,:) = 1.0_rp
    gtrans_limxy(:,:,:,:) = 1.0_rp

    ! calc metrics for orthogonal curvelinear coordinate
    call gtrans_mapfactor

    ! calc coeficient for rotaion of velocity vector
    call gtrans_rotcoef

    ! calc metrics for terrain-following,step-mountain,thin-wall coordinate
    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '*** Terrain coordinate type : ', trim(gtrans_topo_type)
    select case(gtrans_topo_type)
    case('TERRAINFOLLOWING')
      if( io_l ) write(io_fid_log,*) '*** => Terrain-following method'
      call gtrans_terrainfollowing
    case('STEPMOUNTAIN')
      if( io_l ) write(io_fid_log,*) '*** => Step-mountain method'
      call gtrans_thin_wall
      call gtrans_step_mountain
    case('THINWALL')
      if( io_l ) write(io_fid_log,*) '*** => Thin-wall approximation method'
      call gtrans_thin_wall
    case default
       write(*,*) 'xxx Unsupported GTRANS_TOPO_type. STOP'
       call prc_mpistop
    end select

    ! output metrics (for debug)
    call gtrans_write

    return
  end subroutine gtrans_setup

  !-----------------------------------------------------------------------------
  subroutine gtrans_mapfactor
    use scale_mapproj, only: &
       mprj_mapfactor
    use scale_grid_real, only: &
       real_calc_areavol, &
       real_lat,          &
       real_latx,         &
       real_laty,         &
       real_latxy
    implicit none
    !---------------------------------------------------------------------------

    call mprj_mapfactor( real_lat(1:ia,1:ja), gtrans_mapf(:,:,1,i_xy), gtrans_mapf(:,:,2,i_xy))
    call mprj_mapfactor( real_latx(1:ia,1:ja), gtrans_mapf(:,:,1,i_uy), gtrans_mapf(:,:,2,i_uy))
    call mprj_mapfactor( real_laty(1:ia,1:ja), gtrans_mapf(:,:,1,i_xv), gtrans_mapf(:,:,2,i_xv))
    call mprj_mapfactor( real_latxy(1:ia,1:ja), gtrans_mapf(:,:,1,i_uv), gtrans_mapf(:,:,2,i_uv))

    call real_calc_areavol( gtrans_mapf(:,:,:,i_xy) )

    return
  end subroutine gtrans_mapfactor

  !-----------------------------------------------------------------------------
  subroutine gtrans_rotcoef
    use scale_mapproj, only: &
       mprj_rotcoef
    use scale_grid_real, only: &
       real_lon,  &
       real_lat
    implicit none
    !---------------------------------------------------------------------------

    call mprj_rotcoef( gtrans_rotc(:,:,:), & ! [OUT]
                       real_lon(:,:),   & ! [IN]
                       real_lat(:,:)    ) ! [IN]

    return
  end subroutine gtrans_rotcoef

  !-----------------------------------------------------------------------------
  subroutine gtrans_terrainfollowing
    use scale_grid, only: &
       grid_rcdz, &
       grid_rcdx, &
       grid_rcdy, &
       grid_rfdz, &
       grid_rfdx, &
       grid_rfdy
    use scale_grid_real, only: &
       real_cz, &
       real_fz
    use scale_comm, only: &
       comm_vars8, &
       comm_wait
    implicit none

    real(RP) :: REAL_CZ_U (  KA,IA,JA)
    real(RP) :: REAL_CZ_V (  KA,IA,JA)
    real(RP) :: REAL_CZ_UV(  KA,IA,JA)
    real(RP) :: REAL_FZ_U (0:KA,IA,JA)
    real(RP) :: REAL_FZ_V (0:KA,IA,JA)
    real(RP) :: REAL_FZ_UV(0:KA,IA,JA)

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

    ! calc Z-coordinate height at staggered position
    do j = 1, ja
    do i = 1, ia-1
    do k = 1, ka
       real_cz_u(k,i,j) = 0.5_rp * ( real_cz(k,i+1,j) + real_cz(k,i,j) )
    enddo
    enddo
    enddo

    do j = 1, ja
    do i = 1, ia-1
    do k = 0, ka
       real_fz_u(k,i,j) = 0.5_rp * ( real_fz(k,i+1,j) + real_fz(k,i,j) )
    enddo
    enddo
    enddo

    do j = 1, ja-1
    do i = 1, ia
    do k = 1, ka
       real_cz_v(k,i,j) = 0.5_rp * ( real_cz(k,i,j+1) + real_cz(k,i,j) )
    enddo
    enddo
    enddo

    do j = 1, ja-1
    do i = 1, ia
    do k = 0, ka
       real_fz_v(k,i,j) = 0.5_rp * ( real_fz(k,i,j+1) + real_fz(k,i,j) )
    enddo
    enddo
    enddo

    do j = 1, ja-1
    do i = 1, ia-1
    do k = 1, ka
       real_cz_uv(k,i,j) = 0.25_rp * ( real_cz(k,i+1,j+1) + real_cz(k,i+1,j) &
                                    + real_cz(k,i  ,j+1) + real_cz(k,i  ,j) )
    enddo
    enddo
    enddo

    do j = 1, ja-1
    do i = 1, ia-1
    do k = 0, ka
       real_fz_uv(k,i,j) = 0.25_rp * ( real_fz(k,i+1,j+1) + real_fz(k,i+1,j) &
                                    + real_fz(k,i  ,j+1) + real_fz(k,i  ,j) )
    enddo
    enddo
    enddo

    ! G^1/2
    do j = js, je
    do i = is, ie
       ! at (x,y,z)
       do k = 1, ka
          gtrans_gsqrt(k,i,j,i_xyz) = ( real_fz(k,i,j) - real_fz(k-1,i,j) ) * grid_rcdz(k)
       enddo

       ! at (x,y,w)
       do k = 1, ka-1
          gtrans_gsqrt(k,i,j,i_xyw) = ( real_cz(k+1,i,j) - real_cz(k,i,j) ) * grid_rfdz(k)
       enddo
       gtrans_gsqrt(ka,i,j,i_xyw) = gtrans_gsqrt(ka-1,i,j,i_xyw)

       ! at (u,y,w)
       do k = 1, ka-1
          gtrans_gsqrt(k,i,j,i_uyw) = ( real_cz_u(k+1,i,j) - real_cz_u(k,i,j) ) * grid_rfdz(k)
       enddo
       gtrans_gsqrt(ka,i,j,i_uyw) = gtrans_gsqrt(ka-1,i,j,i_uyw)

       ! at (x,v,w)
       do k = 1, ka-1
          gtrans_gsqrt(k,i,j,i_xvw) = ( real_cz_v(k+1,i,j) - real_cz_v(k,i,j) ) * grid_rfdz(k)
       enddo
       gtrans_gsqrt(ka,i,j,i_xvw) = gtrans_gsqrt(ka-1,i,j,i_xvw)

       ! at (u,y,z)
       do k = 1, ka
          gtrans_gsqrt(k,i,j,i_uyz) = ( real_fz_u(k,i,j) - real_fz_u(k-1,i,j) ) * grid_rcdz(k)
       enddo

       ! at (x,v,z)
       do k = 1, ka
          gtrans_gsqrt(k,i,j,i_xvz) = ( real_fz_v(k,i,j) - real_fz_v(k-1,i,j) ) * grid_rcdz(k)
       enddo

       ! at (u,v,z)
       do k = 1, ka
          gtrans_gsqrt(k,i,j,i_uvz) = ( real_fz_uv(k,i,j) - real_fz_uv(k-1,i,j) ) * grid_rcdz(k)
       enddo
    enddo
    enddo

    call comm_vars8( gtrans_gsqrt(:,:,:,1), 1 )
    call comm_vars8( gtrans_gsqrt(:,:,:,2), 2 )
    call comm_vars8( gtrans_gsqrt(:,:,:,3), 3 )
    call comm_vars8( gtrans_gsqrt(:,:,:,4), 4 )
    call comm_vars8( gtrans_gsqrt(:,:,:,5), 5 )
    call comm_vars8( gtrans_gsqrt(:,:,:,6), 6 )
    call comm_vars8( gtrans_gsqrt(:,:,:,7), 7 )
    call comm_wait ( gtrans_gsqrt(:,:,:,1), 1 )
    call comm_wait ( gtrans_gsqrt(:,:,:,2), 2 )
    call comm_wait ( gtrans_gsqrt(:,:,:,3), 3 )
    call comm_wait ( gtrans_gsqrt(:,:,:,4), 4 )
    call comm_wait ( gtrans_gsqrt(:,:,:,5), 5 )
    call comm_wait ( gtrans_gsqrt(:,:,:,6), 6 )
    call comm_wait ( gtrans_gsqrt(:,:,:,7), 7 )

    ! Jacobian * G^1/2
    do j = js, je
    do i = is, ie
    do k = 1,  ka
       gtrans_j13g(k,i,j,i_xyz) = -( real_cz_u(k,i  ,j) - real_cz_u(k,i-1,j) ) * grid_rcdx(i)
       gtrans_j13g(k,i,j,i_xyw) = -( real_fz_u(k,i  ,j) - real_fz_u(k,i-1,j) ) * grid_rcdx(i)
       gtrans_j13g(k,i,j,i_uyw) = -( real_fz(k,i+1,j) - real_fz(k,i  ,j) ) * grid_rfdx(i)
       gtrans_j13g(k,i,j,i_xvw) = -( real_fz_uv(k,i  ,j) - real_fz_uv(k,i-1,j) ) * grid_rcdx(i)
       gtrans_j13g(k,i,j,i_uyz) = -( real_cz(k,i+1,j) - real_cz(k,i  ,j) ) * grid_rfdx(i)
       gtrans_j13g(k,i,j,i_xvz) = -( real_cz_uv(k,i  ,j) - real_cz_uv(k,i-1,j) ) * grid_rcdx(i)
       gtrans_j13g(k,i,j,i_uvz) = -( real_cz_v(k,i+1,j) - real_cz_v(k,i  ,j) ) * grid_rfdx(i)
    enddo
    enddo
    enddo

    do j = js, je
    do i = is, ie
    do k = 1,  ka
       gtrans_j23g(k,i,j,i_xyz) = -( real_cz_v(k,i,j  ) - real_cz_v(k,i,j-1) ) * grid_rcdy(j)
       gtrans_j23g(k,i,j,i_xyw) = -( real_fz_v(k,i,j  ) - real_fz_v(k,i,j-1) ) * grid_rcdy(j)
       gtrans_j23g(k,i,j,i_xvw) = -( real_fz(k,i,j+1) - real_fz(k,i,j  ) ) * grid_rfdy(j)
       gtrans_j23g(k,i,j,i_uyw) = -( real_fz_uv(k,i,j  ) - real_fz_uv(k,i,j-1) ) * grid_rcdy(j)
       gtrans_j23g(k,i,j,i_xvz) = -( real_cz(k,i,j+1) - real_cz(k,i,j  ) ) * grid_rfdy(j)
       gtrans_j23g(k,i,j,i_uyz) = -( real_cz_uv(k,i,j  ) - real_cz_uv(k,i,j-1) ) * grid_rcdy(j)
       gtrans_j23g(k,i,j,i_uvz) = -( real_cz_u(k,i,j+1) - real_cz_u(k,i,j  ) ) * grid_rfdy(j)
    enddo
    enddo
    enddo

    gtrans_j33g = 1.0_rp ! - 1 / G^1/2 * G^1/2

    call comm_vars8( gtrans_j13g(:,:,:,i_xyz),  8 )
    call comm_vars8( gtrans_j13g(:,:,:,i_xyw),  9 )
    call comm_vars8( gtrans_j13g(:,:,:,i_uyw), 10 )
    call comm_vars8( gtrans_j13g(:,:,:,i_xvw), 11 )
    call comm_vars8( gtrans_j13g(:,:,:,i_uyz), 12 )
    call comm_vars8( gtrans_j13g(:,:,:,i_xvz), 13 )
    call comm_vars8( gtrans_j13g(:,:,:,i_uvz), 14 )
    call comm_wait ( gtrans_j13g(:,:,:,i_xyz),  8 )
    call comm_wait ( gtrans_j13g(:,:,:,i_xyw),  9 )
    call comm_wait ( gtrans_j13g(:,:,:,i_uyw), 10 )
    call comm_wait ( gtrans_j13g(:,:,:,i_xvw), 11 )
    call comm_wait ( gtrans_j13g(:,:,:,i_uyz), 12 )
    call comm_wait ( gtrans_j13g(:,:,:,i_xvz), 13 )
    call comm_wait ( gtrans_j13g(:,:,:,i_uvz), 14 )

    call comm_vars8( gtrans_j23g(:,:,:,i_xyz), 15 )
    call comm_vars8( gtrans_j23g(:,:,:,i_xyw), 16 )
    call comm_vars8( gtrans_j23g(:,:,:,i_uyw), 17 )
    call comm_vars8( gtrans_j23g(:,:,:,i_xvw), 18 )
    call comm_vars8( gtrans_j23g(:,:,:,i_uyz), 19 )
    call comm_vars8( gtrans_j23g(:,:,:,i_xvz), 20 )
    call comm_vars8( gtrans_j23g(:,:,:,i_uvz), 21 )
    call comm_wait ( gtrans_j23g(:,:,:,i_xyz), 15 )
    call comm_wait ( gtrans_j23g(:,:,:,i_xyw), 16 )
    call comm_wait ( gtrans_j23g(:,:,:,i_uyw), 17 )
    call comm_wait ( gtrans_j23g(:,:,:,i_xvw), 18 )
    call comm_wait ( gtrans_j23g(:,:,:,i_uyz), 19 )
    call comm_wait ( gtrans_j23g(:,:,:,i_xvz), 20 )
    call comm_wait ( gtrans_j23g(:,:,:,i_uvz), 21 )

    return
  end subroutine gtrans_terrainfollowing

  !-----------------------------------------------------------------------------
  subroutine gtrans_thin_wall
    use scale_process, only: &
       prc_mpistop
    use scale_grid, only: &
       grid_cz, &
       grid_cx, &
       grid_cy, &
       grid_fz, &
       grid_fx, &
       grid_fy
    use scale_topography, only : &
       topo_zsfc
    use scale_comm, only : &
       comm_vars8, &
       comm_wait
    implicit none

    real(RP) :: TOPO_ZsfcALL(2*IA,2*JA)
    real(RP) :: TOPO_ZsfcXY (IA,JA)
    real(RP) :: TOPO_ZsfcUY (IA,JA)
    real(RP) :: TOPO_ZsfcXV (IA,JA)
    real(RP) :: TOPO_ZsfcUV (IA,JA)

    real(RP) :: GTRANS_QLIM (2*KA,2*IA,2*JA,3)
    real(RP) :: QDZ(2*KA)
    real(RP) :: QDX(2*IA)
    real(RP) :: QDY(2*JA)
    real(RP) :: AQAF (3)
    real(RP) :: XSLOPE, YSLOPE
    real(RP) :: Ztop
    real(RP) :: DX_piece, DY_piece
    real(RP) :: DX, DY, DZ

    integer  :: I_QLIMtoLIM(3,7)

    integer  :: iii, jjj, n
    integer  :: k, i, j, kk, ii, jj
    !---------------------------------------------------------------------------

    ! calc absolute height at staggered position
    ! at (x,y)
    do j = 1, ja
    do i = 1, ia
      topo_zsfcxy(i,j) = topo_zsfc(i,j)
    enddo
    enddo
    ! at (u,y)
    do j = 1, ja
    do i = 1, ia-1
      topo_zsfcuy(i,j) = 0.5_rp * ( topo_zsfc(i,j) + topo_zsfc(i+1,j) )
    enddo
    enddo
    ! at (x,v)
    do j = 1, ja-1
    do i = 1, ia
      topo_zsfcxv(i,j) = 0.5_rp * ( topo_zsfc(i,j) + topo_zsfc(i,j+1) )
    enddo
    enddo
    ! at (u,v)
    do j = 1, ja-1
    do i = 1, ia-1
      topo_zsfcuv(i,j) = 0.25_rp * ( topo_zsfc(i  ,j  ) + topo_zsfc(i  ,j+1) &
                                   + topo_zsfc(i  ,j+1) + topo_zsfc(i+1,j+1) )
    enddo
    enddo

    ! reset topography
    topo_zsfc(:,:) = 0.d0

    call comm_vars8( topo_zsfcxy(:,:), 1 )
    call comm_vars8( topo_zsfcuy(:,:), 2 )
    call comm_vars8( topo_zsfcxv(:,:), 3 )
    call comm_vars8( topo_zsfcuv(:,:), 4 )
    call comm_wait ( topo_zsfcxy(:,:), 1 )
    call comm_wait ( topo_zsfcuy(:,:), 2 )
    call comm_wait ( topo_zsfcxv(:,:), 3 )
    call comm_wait ( topo_zsfcuv(:,:), 4 )

    ! all height
    do j = 1, ja
    do i = 1, ia
       ii = (i-1) * 2 + 1
       jj = (j-1) * 2 + 1

       topo_zsfcall(ii  ,jj  ) = topo_zsfcxy(i,j)
       topo_zsfcall(ii+1,jj  ) = topo_zsfcuy(i,j)
       topo_zsfcall(ii  ,jj+1) = topo_zsfcxv(i,j)
       topo_zsfcall(ii+1,jj+1) = topo_zsfcuv(i,j)
    enddo
    enddo

    ! length of control volume of quarter cell
    do k = 1, ka
       kk = (k-1) * 2 + 1

       qdz(kk  ) = grid_cz(k) - grid_fz(k-1)
       qdz(kk+1) = grid_fz(k) - grid_cz(k  )
    enddo

    do i = 1, ia-1
       ii = (i-1) * 2 + 1

       qdx(ii  ) = grid_fx(i  ) - grid_cx(i)
       qdx(ii+1) = grid_cx(i+1) - grid_fx(i)
    enddo

    do j = 1, ja-1
       jj = (j-1) * 2 + 1

       qdy(jj  ) = grid_fy(j  ) - grid_cy(j)
       qdy(jj+1) = grid_cy(j+1) - grid_fy(j)
    enddo

    ! quarter flux limiter
    do jj = 1, 2*(ja-1)
    do ii = 1, 2*(ia-1)
    do kk = ks, 2*ke
       dx_piece = qdx(ii) / real(GTRANS_ThinWall_XDIV,kind=rp)
       dy_piece = qdy(jj) / real(GTRANS_ThinWall_YDIV,kind=rp)

       aqaf(1:3) = 0.0_rp
       ztop = sum(qdz(ks:kk))

       !--- y-z face ---
       yslope = ( topo_zsfcall(ii,jj+1) - topo_zsfcall(ii,jj) ) / qdy(jj)

       do jjj = 1, gtrans_thinwall_ydiv
          dy = ( real(jjj,kind=RP) - 0.5_RP ) * DY_piece
          dz = ztop - topo_zsfcall(ii,jj) - yslope * dy

          if ( dz > 0.0_rp ) then
             if ( dz < qdz(kk) ) then
                aqaf(i_fyz) = aqaf(i_fyz) + dz      * dy_piece
             else
                aqaf(i_fyz) = aqaf(i_fyz) + qdz(kk) * dy_piece
             endif
          endif
       enddo

       !--- x-z face ---
       xslope = ( topo_zsfcall(ii+1,jj) - topo_zsfcall(ii,jj) ) / qdx(ii)

       do iii = 1, gtrans_thinwall_xdiv
          dx = ( real(iii,kind=RP) - 0.5_RP ) * DX_piece
          dz = ztop - topo_zsfcall(ii,jj) + xslope * dx

          if ( dz > 0.0_rp ) then
             if ( dz < qdz(kk) ) then
                aqaf(i_fxz) = aqaf(i_fxz) + dz      * dx_piece
             else
                aqaf(i_fxz) = aqaf(i_fxz) + qdz(kk) * dx_piece
             endif
          endif
       enddo

       !--- x-y face ---
       do jjj = 1, gtrans_thinwall_ydiv
       do iii = 1, gtrans_thinwall_xdiv
          dx = ( real(iii,kind=RP) - 0.5_RP ) * DX_piece
          dy = ( real(jjj,kind=RP) - 0.5_RP ) * DY_piece
          dz = ztop - topo_zsfcall(ii,jj) - xslope * dx - yslope * dy

          if ( dz > 0.0_rp ) then
             aqaf(i_fxy) = aqaf(i_fxy) + dx_piece * dy_piece
          endif
       enddo
       enddo

       gtrans_qlim(kk,ii,jj,i_fyz) = aqaf(i_fyz) / ( qdy(jj) * qdz(kk) )
       gtrans_qlim(kk,ii,jj,i_fxz) = aqaf(i_fxz) / ( qdx(ii) * qdz(kk) )
       gtrans_qlim(kk,ii,jj,i_fxy) = aqaf(i_fxy) / ( qdy(jj) * qdx(ii) )
    enddo
    enddo
    enddo

    ! index i,j,k
    i_qlimtolim(1:3,i_xyz) = (/ 1, 1, 1 /)
    i_qlimtolim(1:3,i_xyw) = (/ 1, 1, 0 /)
    i_qlimtolim(1:3,i_uyw) = (/ 0, 1, 0 /)
    i_qlimtolim(1:3,i_xvw) = (/ 1, 0, 0 /)
    i_qlimtolim(1:3,i_uyz) = (/ 0, 1, 1 /)
    i_qlimtolim(1:3,i_xvz) = (/ 1, 0, 1 /)
    i_qlimtolim(1:3,i_uvz) = (/ 0, 0, 1 /)

    do n = 1, 7
       do j = js, je
       do i = is, ie
       do k = ks, ke
          ii = (i-1) * 2 + 1 - i_qlimtolim(1,n)
          jj = (j-1) * 2 + 1 - i_qlimtolim(2,n)
          kk = (k-1) * 2 + 1 - i_qlimtolim(3,n)

         gtrans_limyz(k,i,j,n) = 0.25_rp * ( gtrans_qlim(kk  ,ii,jj,i_fyz) + gtrans_qlim(kk  ,ii  ,jj+1,i_fyz) &
                                           + gtrans_qlim(kk+1,ii,jj,i_fyz) + gtrans_qlim(kk+1,ii  ,jj+1,i_fyz) )
         gtrans_limxz(k,i,j,n) = 0.25_rp * ( gtrans_qlim(kk  ,ii,jj,i_fxz) + gtrans_qlim(kk  ,ii+1,jj  ,i_fxz) &
                                           + gtrans_qlim(kk+1,ii,jj,i_fxz) + gtrans_qlim(kk+1,ii+1,jj  ,i_fxz) )
         gtrans_limxy(k,i,j,n) = 0.25_rp * ( gtrans_qlim(kk  ,ii,jj,i_fxy) + gtrans_qlim(kk  ,ii+1,jj+1,i_fxy) &
                                           + gtrans_qlim(kk  ,ii,jj,i_fxy) + gtrans_qlim(kk  ,ii+1,jj  ,i_fxy) )

         if ( gtrans_limyz(k,i,j,n) > 1.d0 ) then
            write(*,*) 'xxx Facter miss! Check!'
            write(*,*) k,i,j,n,gtrans_limyz(k,i,j,n)
            call prc_mpistop
         endif
       enddo
       enddo
       enddo
    enddo

    call comm_vars8( gtrans_limyz(:,:,:,i_xyz),  1 )
    call comm_vars8( gtrans_limyz(:,:,:,i_xyw),  2 )
    call comm_vars8( gtrans_limyz(:,:,:,i_uyw),  3 )
    call comm_vars8( gtrans_limyz(:,:,:,i_xvw),  4 )
    call comm_vars8( gtrans_limyz(:,:,:,i_uyz),  5 )
    call comm_vars8( gtrans_limyz(:,:,:,i_xvz),  6 )
    call comm_vars8( gtrans_limyz(:,:,:,i_uvz),  7 )
    call comm_wait ( gtrans_limyz(:,:,:,i_xyz),  1 )
    call comm_wait ( gtrans_limyz(:,:,:,i_xyw),  2 )
    call comm_wait ( gtrans_limyz(:,:,:,i_uyw),  3 )
    call comm_wait ( gtrans_limyz(:,:,:,i_xvw),  4 )
    call comm_wait ( gtrans_limyz(:,:,:,i_uyz),  5 )
    call comm_wait ( gtrans_limyz(:,:,:,i_xvz),  6 )
    call comm_wait ( gtrans_limyz(:,:,:,i_uvz),  7 )

    call comm_vars8( gtrans_limxz(:,:,:,i_xyz),  8 )
    call comm_vars8( gtrans_limxz(:,:,:,i_xyw),  9 )
    call comm_vars8( gtrans_limxz(:,:,:,i_uyw), 10 )
    call comm_vars8( gtrans_limxz(:,:,:,i_xvw), 11 )
    call comm_vars8( gtrans_limxz(:,:,:,i_uyz), 12 )
    call comm_vars8( gtrans_limxz(:,:,:,i_xvz), 13 )
    call comm_vars8( gtrans_limxz(:,:,:,i_uvz), 14 )
    call comm_wait ( gtrans_limxz(:,:,:,i_xyz),  8 )
    call comm_wait ( gtrans_limxz(:,:,:,i_xyw),  9 )
    call comm_wait ( gtrans_limxz(:,:,:,i_uyw), 10 )
    call comm_wait ( gtrans_limxz(:,:,:,i_xvw), 11 )
    call comm_wait ( gtrans_limxz(:,:,:,i_uyz), 12 )
    call comm_wait ( gtrans_limxz(:,:,:,i_xvz), 13 )
    call comm_wait ( gtrans_limxz(:,:,:,i_uvz), 14 )

    call comm_vars8( gtrans_limxy(:,:,:,i_xyz), 15 )
    call comm_vars8( gtrans_limxy(:,:,:,i_xyw), 16 )
    call comm_vars8( gtrans_limxy(:,:,:,i_uyw), 17 )
    call comm_vars8( gtrans_limxy(:,:,:,i_xvw), 18 )
    call comm_vars8( gtrans_limxy(:,:,:,i_uyz), 19 )
    call comm_vars8( gtrans_limxy(:,:,:,i_xvz), 20 )
    call comm_vars8( gtrans_limxy(:,:,:,i_uvz), 21 )
    call comm_wait ( gtrans_limxy(:,:,:,i_xyz), 15 )
    call comm_wait ( gtrans_limxy(:,:,:,i_xyw), 16 )
    call comm_wait ( gtrans_limxy(:,:,:,i_uyw), 17 )
    call comm_wait ( gtrans_limxy(:,:,:,i_xvw), 18 )
    call comm_wait ( gtrans_limxy(:,:,:,i_uyz), 19 )
    call comm_wait ( gtrans_limxy(:,:,:,i_xvz), 20 )
    call comm_wait ( gtrans_limxy(:,:,:,i_uvz), 21 )

    return
  end subroutine gtrans_thin_wall

  !-----------------------------------------------------------------------------
  subroutine gtrans_step_mountain
    use scale_comm, only: &
       comm_vars8, &
       comm_wait
    implicit none

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

    do n = 1, 7
       do j = js, je
       do i = is, ie
       do k = ks, ke

         if ( gtrans_limyz(k,i,j,n) > 0.0_rp ) then
            gtrans_limyz(k,i,j,n) = 1.0_rp
         endif

         if ( gtrans_limxz(k,i,j,n) > 0.0_rp ) then
            gtrans_limxz(k,i,j,n) = 1.0_rp
         endif

         if ( gtrans_limxy(k,i,j,n) > 0.0_rp ) then
            gtrans_limxy(k,i,j,n) = 1.0_rp
         endif
       enddo
       enddo
       enddo
    enddo

    call comm_vars8( gtrans_limyz(:,:,:,i_xyz),  1 )
    call comm_vars8( gtrans_limyz(:,:,:,i_xyw),  2 )
    call comm_vars8( gtrans_limyz(:,:,:,i_uyw),  3 )
    call comm_vars8( gtrans_limyz(:,:,:,i_xvw),  4 )
    call comm_vars8( gtrans_limyz(:,:,:,i_uyz),  5 )
    call comm_vars8( gtrans_limyz(:,:,:,i_xvz),  6 )
    call comm_vars8( gtrans_limyz(:,:,:,i_uvz),  7 )
    call comm_wait ( gtrans_limyz(:,:,:,i_xyz),  1 )
    call comm_wait ( gtrans_limyz(:,:,:,i_xyw),  2 )
    call comm_wait ( gtrans_limyz(:,:,:,i_uyw),  3 )
    call comm_wait ( gtrans_limyz(:,:,:,i_xvw),  4 )
    call comm_wait ( gtrans_limyz(:,:,:,i_uyz),  5 )
    call comm_wait ( gtrans_limyz(:,:,:,i_xvz),  6 )
    call comm_wait ( gtrans_limyz(:,:,:,i_uvz),  7 )

    call comm_vars8( gtrans_limxz(:,:,:,i_xyz),  8 )
    call comm_vars8( gtrans_limxz(:,:,:,i_xyw),  9 )
    call comm_vars8( gtrans_limxz(:,:,:,i_uyw), 10 )
    call comm_vars8( gtrans_limxz(:,:,:,i_xvw), 11 )
    call comm_vars8( gtrans_limxz(:,:,:,i_uyz), 12 )
    call comm_vars8( gtrans_limxz(:,:,:,i_xvz), 13 )
    call comm_vars8( gtrans_limxz(:,:,:,i_uvz), 14 )
    call comm_wait ( gtrans_limxz(:,:,:,i_xyz),  8 )
    call comm_wait ( gtrans_limxz(:,:,:,i_xyw),  9 )
    call comm_wait ( gtrans_limxz(:,:,:,i_uyw), 10 )
    call comm_wait ( gtrans_limxz(:,:,:,i_xvw), 11 )
    call comm_wait ( gtrans_limxz(:,:,:,i_uyz), 12 )
    call comm_wait ( gtrans_limxz(:,:,:,i_xvz), 13 )
    call comm_wait ( gtrans_limxz(:,:,:,i_uvz), 14 )

    call comm_vars8( gtrans_limxy(:,:,:,i_xyz), 15 )
    call comm_vars8( gtrans_limxy(:,:,:,i_xyw), 16 )
    call comm_vars8( gtrans_limxy(:,:,:,i_uyw), 17 )
    call comm_vars8( gtrans_limxy(:,:,:,i_xvw), 18 )
    call comm_vars8( gtrans_limxy(:,:,:,i_uyz), 19 )
    call comm_vars8( gtrans_limxy(:,:,:,i_xvz), 20 )
    call comm_vars8( gtrans_limxy(:,:,:,i_uvz), 21 )
    call comm_wait ( gtrans_limxy(:,:,:,i_xyz), 15 )
    call comm_wait ( gtrans_limxy(:,:,:,i_xyw), 16 )
    call comm_wait ( gtrans_limxy(:,:,:,i_uyw), 17 )
    call comm_wait ( gtrans_limxy(:,:,:,i_xvw), 18 )
    call comm_wait ( gtrans_limxy(:,:,:,i_uyz), 19 )
    call comm_wait ( gtrans_limxy(:,:,:,i_xvz), 20 )
    call comm_wait ( gtrans_limxy(:,:,:,i_uvz), 21 )

    return
  end subroutine gtrans_step_mountain

  !-----------------------------------------------------------------------------
  subroutine gtrans_write
    use scale_const, only: &
       const_radius
    use scale_vector, only: &
       vectr_distance
    use scale_fileio, only: &
       fileio_write
    use scale_grid_real, only: &
       real_basepoint_lon, &
       real_basepoint_lat, &
       real_lon,           &
       real_lat
    use scale_mapproj, only: &
       mprj_lonlat2xy
    implicit none

    real(RP) :: check_X_XY(IA,JA)
    real(RP) :: check_Y_XY(IA,JA)
    real(RP) :: distance  (IA,JA)

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

    if ( gtrans_out_basename /= '' ) then

       if( io_l ) write(io_fid_log,*)
       if( io_l ) write(io_fid_log,*) '*** Output metrics file ***'

       call fileio_write( gtrans_mapf(:,:,1,i_xy),       gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'MAPF_X_XY', 'Map factor x-dir at XY', 'NIL', 'XY', gtrans_out_dtype  ) ! [IN]
       call fileio_write( gtrans_mapf(:,:,2,i_xy),       gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'MAPF_Y_XY', 'Map factor y-dir at XY', 'NIL', 'XY', gtrans_out_dtype  ) ! [IN]
       call fileio_write( gtrans_mapf(:,:,1,i_uy),       gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'MAPF_X_UY', 'Map factor x-dir at UY', 'NIL', 'UY', gtrans_out_dtype  ) ! [IN]
       call fileio_write( gtrans_mapf(:,:,2,i_uy),       gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'MAPF_Y_UY', 'Map factor y-dir at UY', 'NIL', 'UY', gtrans_out_dtype  ) ! [IN]
       call fileio_write( gtrans_mapf(:,:,1,i_xv),       gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'MAPF_X_XV', 'Map factor x-dir at XV', 'NIL', 'XV', gtrans_out_dtype  ) ! [IN]
       call fileio_write( gtrans_mapf(:,:,2,i_xv),       gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'MAPF_Y_XV', 'Map factor y-dir at XV', 'NIL', 'XV', gtrans_out_dtype  ) ! [IN]
       call fileio_write( gtrans_mapf(:,:,1,i_uv),       gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'MAPF_X_UV', 'Map factor x-dir at UV', 'NIL', 'UV', gtrans_out_dtype  ) ! [IN]
       call fileio_write( gtrans_mapf(:,:,2,i_uv),       gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'MAPF_Y_UV', 'Map factor y-dir at UV', 'NIL', 'UV', gtrans_out_dtype  ) ! [IN]

       call fileio_write( gtrans_rotc(:,:,1),            gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'ROTC_COS',  'Rotation factor (cos)',  'NIL', 'XY', gtrans_out_dtype  ) ! [IN]
       call fileio_write( gtrans_rotc(:,:,2),            gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'ROTC_SIN',  'Rotation factor (sin)',  'NIL', 'XY', gtrans_out_dtype  ) ! [IN]

       call fileio_write( gtrans_rotc(:,:,1),            gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'ROTC_COS',  'Rotation factor (cos)',  'NIL', 'XY', gtrans_out_dtype  ) ! [IN]

       do j = 1, ja
       do i = 1, ia
          call mprj_lonlat2xy( real_lon(i,j), real_lat(i,j), check_x_xy(i,j), check_y_xy(i,j) )
       enddo
       enddo

       call fileio_write( check_x_xy(:,:),     gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'X_XY', 'x at XY for check', 'NIL', 'XY', gtrans_out_dtype  ) ! [IN]
       call fileio_write( check_y_xy(:,:),     gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'Y_XY', 'y at XY for check', 'NIL', 'XY', gtrans_out_dtype  ) ! [IN]

       do j = 1, ja
       do i = 1, ia
          call vectr_distance( const_radius,       & ! [IN]
                               real_basepoint_lon, & ! [IN]
                               real_basepoint_lat, & ! [IN]
                               real_lon(i,j),      & ! [IN]
                               real_lat(i,j),      & ! [IN]
                               distance(i,j)       ) ! [OUT]
       enddo
       enddo

       call fileio_write( distance(:,:),               gtrans_out_basename, gtrans_out_title, & ! [IN]
                          'distance', 'distance from basepoint', 'm', 'XY', gtrans_out_dtype  ) ! [IN]

    endif

    return
  end subroutine gtrans_write

end module scale_gridtrans