!= 配列の分割／拡張
!
!= Division and expansion of arrays
!
! Authors::   Hiroki KASHIMURA, Yoshiyuki O. TAKAHASHI
! Version::
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2013. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!
module sltt_extarr
  != 配列の分割／拡張
  !
  != Division and expansion of arrays


  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision.
    &                 STRING     ! 文字列.       Strings.

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only:       &
    &                imax, & ! 経度格子点数.
                             ! Number of grid points in longitude
    &                jmax, & ! 緯度格子点数.
                             ! Number of grid points in latitude
    &                kmax !, & ! 鉛直層数.
                             ! Number of vertical level
  ! 組成に関わる配列の設定
  ! Settings of array for atmospheric composition
  !
  use composition, only : ncmax

  implicit none

  private

  !
  ! local variables
  !
!!$  real(DP), save, allocatable :: a_LQIFUVNP(:)
!!$  real(DP), save, allocatable :: a_LQIFUVSP(:)

  real(DP), save, allocatable :: a_LCIFUVNP(:)
  real(DP), save, allocatable :: a_LCIFUVSP(:)

  real(DP), save              :: SumSinSq


  real(DP), save              :: sum_sinsq
  real(DP), save, allocatable :: asinlon1( : ), acoslon1( : )


  public :: SLTTExtArrInit
  public :: SLTTExtArrExt
  public :: SLTTExtArrExt2


  character(*), parameter:: module_name = 'sltt_extarr'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: sltt_extarr.f90,v 1.3 2014/02/18 02:59:19 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

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

contains

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

  subroutine SLTTExtArrInit(                      &
    & x_LonS, y_LatS, x_LonN, y_LatN,             & ! (in )
    & x_ExtLonS, y_ExtLatS, x_ExtLonN, y_ExtLatN  & ! (out)
    & )

    !
    ! MPI
    !
    use mpi_wrapper   , only : myrank, nprocs, &
      & MPIWrapperISend, MPIWrapperIRecv, MPIWrapperWait

    use constants0, only : PI
    use axesset   , only : y_Lat
    use sltt_const, only : PIx2, PIH, &
      & dtjw, iexmin, iexmax, jexmins, jexmaxs, jexminn, jexmaxn

    real(DP), intent(in ) :: x_LonS   (   0:imax-1    )
    real(DP), intent(in ) :: y_LatS   (   1:jmax/2    )
    real(DP), intent(in ) :: x_LonN   (   0:imax-1    )
    real(DP), intent(in ) :: y_LatN   (   1:jmax/2    )
    real(DP), intent(out) :: x_ExtLonS(iexmin :iexmax )
    real(DP), intent(out) :: y_ExtLatS(jexmins:jexmaxs)
    real(DP), intent(out) :: x_ExtLonN(iexmin :iexmax )
    real(DP), intent(out) :: y_ExtLatN(jexminn:jexmaxn)


    !
    ! local variables
    !
    integer  :: idest
    integer  :: idep
    integer  :: a_ireq_send_ss
    integer  :: a_ireq_send_sn
    integer  :: a_ireq_send_ns
    integer  :: a_ireq_send_nn
    integer  :: a_ireq_recv_ss
    integer  :: a_ireq_recv_sn
    integer  :: a_ireq_recv_ns
    integer  :: a_ireq_recv_nn
    ! SS : Southern hemisphere, Southward array
    ! SN : Southern hemisphere, Northward array
    ! NS : Northern hemisphere, Southward array
    ! NN : Northern hemisphere, Northward array
    real(DP), allocatable :: y_LatSN   (:)
    real(DP), allocatable :: y_LatNS   (:)
    real(DP), allocatable :: y_LatSS   (:)
    real(DP), allocatable :: y_LatNN   (:)

    real(DP) :: h, theta, thetasq

    integer  :: i, j, k, m

    logical, save :: sw_fs
    data sw_fs /.true./


    if( .not. sw_fs ) return
    sw_fs = .false.


    x_ExtLonS(-2) = 0.0_DP - ( x_LonS(2) - x_LonS(0) )
    x_ExtLonS(-1) = 0.0_DP - ( x_LonS(1) - x_LonS(0) )
    do i = 0, imax-1
      x_ExtLonS(i) = x_LonS(i)
    end do
    x_ExtLonS(imax-1+1) = PIx2
    x_ExtLonS(imax-1+2) = PIx2 + ( x_LonS(1) - x_LonS(0) )
    x_ExtLonS(imax-1+3) = PIx2 + ( x_LonS(2) - x_LonS(0) )
    !
    x_ExtLonN(-2) = 0.0_DP - ( x_LonN(2) - x_LonN(0) )
    x_ExtLonN(-1) = 0.0_DP - ( x_LonN(1) - x_LonN(0) )
    do i = 0, imax-1
      x_ExtLonN(i) = x_LonN(i)
    end do
    x_ExtLonN(imax-1+1) = PIx2
    x_ExtLonN(imax-1+2) = PIx2 + ( x_LonN(1) - x_LonN(0) )
    x_ExtLonN(imax-1+3) = PIx2 + ( x_LonN(2) - x_LonN(0) )


    !====================================================================

    if ( myrank > 0          ) then
      allocate( y_LatSN   (1:dtjw) )
      allocate( y_LatNS   (1:dtjw) )
    end if
    if ( myrank < (nprocs-1) ) then
      allocate( y_LatSS   (1:dtjw) )
      allocate( y_LatNN   (1:dtjw) )
    end if


    ! y_Lat(1:dtjw) values are transfered (y_LatSN   ).
    !
    if ( myrank < (nprocs-1) ) then
      idest = myrank + 1
      call MPIWrapperISend(               &
        & idest, dtjw,                    & ! (in)
        & y_Lat(1:dtjw),                  & ! (in)
        & a_ireq_send_sn                  & ! (out)
        & )
    end if
    if ( myrank > 0 ) then
      idep = myrank - 1
      call MPIWrapperIRecv( &
        & idep, dtjw,                     & ! (in)
        & y_LatSN   ,                     & ! (out)
        & a_ireq_recv_sn                  & ! (out)
        & )
    end if
    !
    ! y_Lat(jmax/2+1-dtjw:jmax/2) values are transfered (y_LatSS   ).
    !
    if ( myrank > 0 ) then
      idest = myrank - 1
      call MPIWrapperISend(               &
        & idest, dtjw,                    & ! (in)
        & y_Lat(jmax/2+1-dtjw:jmax/2),    & ! (in)
        & a_ireq_send_ss                  & ! (out)
        & )
    end if
    if ( myrank < (nprocs-1) ) then
      idep = myrank + 1
      call MPIWrapperIRecv(               &
        & idep, dtjw,                     & ! (in)
        & y_LatSS   ,                     & ! (out)
        & a_ireq_recv_ss                  & ! (out)
        & )
    end if
    !
    ! y_Lat(jmax/2+1:jmax/2+dtjw) values are transfered (y_LatNN   ).
    !
    if ( myrank > 0 ) then
      idest = myrank - 1
      call MPIWrapperISend(               &
        & idest, dtjw,                    & ! (in)
        & y_Lat(jmax/2+1:jmax/2+dtjw),    & ! (in)
        & a_ireq_send_nn                  & ! (out)
        & )
    end if
    if ( myrank < (nprocs-1) ) then
      idep = myrank + 1
      call MPIWrapperIRecv(               &
        & idep, dtjw,                     & ! (in)
        & y_LatNN   ,                     & ! (out)
        & a_ireq_recv_nn                  & ! (out)
        & )
    end if
    !
    ! y_Lat(jmax+1-dtjw:jmax) values are transfered (y_LatNS   ).
    !
    if ( myrank < (nprocs-1) ) then
      idest = myrank + 1
      call MPIWrapperISend(               &
        & idest, dtjw,                    & ! (in)
        & y_Lat(jmax+1-dtjw:jmax),        & ! (in)
        & a_ireq_send_ns                  & ! (out)
        & )
    end if
    if ( myrank > 0 ) then
      idep = myrank - 1
      call MPIWrapperIRecv( &
        & idep, dtjw,                     & ! (in)
        & y_LatNS   ,                     & ! (out)
        & a_ireq_recv_ns                  & ! (out)
        & )
    end if

    if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_send_sn )
    if ( myrank > 0          ) call MPIWrapperWait( a_ireq_recv_sn )
    if ( myrank > 0          ) call MPIWrapperWait( a_ireq_send_ss )
    if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_recv_ss )
    if ( myrank > 0          ) call MPIWrapperWait( a_ireq_send_nn )
    if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_recv_nn )
    if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_send_ns )
    if ( myrank > 0          ) call MPIWrapperWait( a_ireq_recv_ns )


    do j = 1, jmax/2
      y_ExtLatS(j) = y_LatS(j)
      y_ExtLatN(j) = y_LatN(j)
    end do


    ! southern edge of southern array
    if( myrank == (nprocs-1) ) then
      y_ExtLatS(0) = -PIH
      do j = 1, dtjw
        y_ExtLatS(0-j) = -PIH - ( y_LatS(j) - ( -PIH ) )
      end do
    else
      do j = 1, dtjw
        y_ExtLatS(1-j) = y_LatSS   (dtjw-(j-1))
      end do
    end if
    ! northern edge of southern array
    if ( myrank == 0 ) then
      do j = 1, dtjw
        y_ExtLatS(jmax/2+j) = y_LatN(j)
      end do
    else
      do j = 1, dtjw
        y_ExtLatS(jmax/2+j) = y_LatSN   (j)
      end do
    end if
    !
    ! southern edge of northern array
    if ( myrank == 0 ) then
      do j = 1, dtjw
        y_ExtLatN(-j+1) = y_LatS(jmax/2-(j-1))
      end do
    else
      do j = 1, dtjw
        y_ExtLatN(-j+1) = y_LatNS   (dtjw-(j-1))
      end do
    end if
    ! northern edge of northern array
    if( myrank == (nprocs-1) ) then
      y_ExtLatN(jmax/2+1) = PIH
      do j = 1, dtjw
        y_ExtLatN(jmax/2+1+j) = PIH + ( PIH - y_LatN(jmax/2+1-j) )
      end do
    else
      do j = 1, dtjw
        y_ExtLatN(jmax/2+j) = y_LatNN   (j)
      end do
    end if

    if ( myrank > 0          ) then
      deallocate( y_LatSN    )
      deallocate( y_LatNS    )
    end if
    if ( myrank < (nprocs-1) ) then
      deallocate( y_LatSS    )
      deallocate( y_LatNN    )
    end if



!!$    allocate( a_LQIFUVNP(jmax/2-1:jmax/2+2) )
!!$    allocate( a_LQIFUVSP(      -1:2       ) )

    allocate( a_LCIFUVNP(jmax/2-1:jmax/2+3) )
    allocate( a_LCIFUVSP(      -2:2       ) )



    !
    ! calculation of factors for Lagrange cubic interpolation used 
    ! to estimate mixing ratios at poles
    !
    if( myrank == (nprocs-1) ) then

!!$      !
!!$      ! calculation of factors for Lagrange quadratic interpolation used 
!!$      ! to estimate wind velocities at south pole (old)
!!$      !
!!$      a_LQIFUVSP(-1) =                                                              &
!!$        &     ( y_ExtLatS( 0) - y_ExtLatS( 1) ) * ( y_ExtLatS( 0) - y_ExtLatS( 2) )   &
!!$        & / ( ( y_ExtLatS(-1) - y_ExtLatS( 1) ) * ( y_ExtLatS(-1) - y_ExtLatS( 2) ) )
!!$      a_LQIFUVSP( 0) = 1.0d100
!!$      a_LQIFUVSP( 1) =                                                              &
!!$        &     ( y_ExtLatS( 0) - y_ExtLatS(-1) ) * ( y_ExtLatS( 0) - y_ExtLatS( 2) )   &
!!$        & / ( ( y_ExtLatS( 1) - y_ExtLatS(-1) ) * ( y_ExtLatS( 1) - y_ExtLatS( 2) ) )
!!$      a_LQIFUVSP( 2) =                                                              &
!!$        &     ( y_ExtLatS( 0) - y_ExtLatS(-1) ) * ( y_ExtLatS( 0) - y_ExtLatS( 1) )   &
!!$        & / ( ( y_ExtLatS( 2) - y_ExtLatS(-1) ) * ( y_ExtLatS( 2) - y_ExtLatS( 1) ) )
!!$      !
!!$      ! calculation of factors for Lagrange quadratic interpolation used 
!!$      ! to estimate wind velocities at north pole (old)
!!$      !
!!$      a_LQIFUVNP(jmax/2-1) =                                                                              &
!!$        &     ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2  ) ) * ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2+2) ) &
!!$        & / ( ( y_ExtLatN(jmax/2-1) - y_ExtLatN(jmax/2  ) ) * ( y_ExtLatN(jmax/2-1) - y_ExtLatN(jmax/2+2) ) )
!!$      a_LQIFUVNP(jmax/2  ) =                                                                              &
!!$        &     ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2-1) ) * ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2+2) ) &
!!$        & / ( ( y_ExtLatN(jmax/2  ) - y_ExtLatN(jmax/2-1) ) * ( y_ExtLatN(jmax/2  ) - y_ExtLatN(jmax/2+2) ) )
!!$      a_LQIFUVNP(jmax/2+1) = 1.0d100
!!$      a_LQIFUVNP(jmax/2+2) =                                                                              &
!!$        &     ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2-1) ) * ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2  ) ) &
!!$        & / ( ( y_ExtLatN(jmax/2+2) - y_ExtLatN(jmax/2-1) ) * ( y_ExtLatN(jmax/2+2) - y_ExtLatN(jmax/2  ) ) )


      !
      ! calculation of factors for Lagrange cubic interpolation used 
      ! to estimate wind velocities at south pole
      !
      a_LCIFUVSP(-2) =                          &
        &     ( y_ExtLatS( 0) - y_ExtLatS(-1) ) &
        &   * ( y_ExtLatS( 0) - y_ExtLatS( 1) ) &
        &   * ( y_ExtLatS( 0) - y_ExtLatS( 2) ) &
        & / ( ( y_ExtLatS(-2) - y_ExtLatS(-1) ) &
        &   * ( y_ExtLatS(-2) - y_ExtLatS( 1) ) &
        &   * ( y_ExtLatS(-2) - y_ExtLatS( 2) ) )
      a_LCIFUVSP(-1) =                          &
        &     ( y_ExtLatS( 0) - y_ExtLatS(-2) ) &
        &   * ( y_ExtLatS( 0) - y_ExtLatS( 1) ) &
        &   * ( y_ExtLatS( 0) - y_ExtLatS( 2) ) &
        & / ( ( y_ExtLatS(-1) - y_ExtLatS(-2) ) &
        &   * ( y_ExtLatS(-1) - y_ExtLatS( 1) ) &
        &   * ( y_ExtLatS(-1) - y_ExtLatS( 2) ) )
      a_LCIFUVSP( 0) = 1.0d100
      a_LCIFUVSP( 1) =                          &
        &     ( y_ExtLatS( 0) - y_ExtLatS(-2) ) &
        &   * ( y_ExtLatS( 0) - y_ExtLatS(-1) ) &
        &   * ( y_ExtLatS( 0) - y_ExtLatS( 2) ) &
        & / ( ( y_ExtLatS( 1) - y_ExtLatS(-2) ) &
        &   * ( y_ExtLatS( 1) - y_ExtLatS(-1) ) &
        &   * ( y_ExtLatS( 1) - y_ExtLatS( 2) ) )
      a_LCIFUVSP( 2) =                          &
        &     ( y_ExtLatS( 0) - y_ExtLatS(-2) ) &
        &   * ( y_ExtLatS( 0) - y_ExtLatS(-1) ) &
        &   * ( y_ExtLatS( 0) - y_ExtLatS( 1) ) &
        & / ( ( y_ExtLatS( 2) - y_ExtLatS(-2) ) &
        &   * ( y_ExtLatS( 2) - y_ExtLatS(-1) ) &
        &   * ( y_ExtLatS( 2) - y_ExtLatS( 1) ) )
      !
      ! calculation of factors for Lagrange cubic interpolation used 
      ! to estimate wind velocities at north pole
      !
      a_LCIFUVNP(jmax/2-1) =                                &
        &     ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2  ) ) &
        &   * ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2+2) ) &
        &   * ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2+3) ) &
        & / ( ( y_ExtLatN(jmax/2-1) - y_ExtLatN(jmax/2  ) ) &
        &   * ( y_ExtLatN(jmax/2-1) - y_ExtLatN(jmax/2+2) ) &
        &   * ( y_ExtLatN(jmax/2-1) - y_ExtLatN(jmax/2+3) ) )
      a_LCIFUVNP(jmax/2  ) =                                &
        &     ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2-1) ) &
        &   * ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2+2) ) &
        &   * ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2+3) ) &
        & / ( ( y_ExtLatN(jmax/2  ) - y_ExtLatN(jmax/2-1) ) &
        &   * ( y_ExtLatN(jmax/2  ) - y_ExtLatN(jmax/2+2) ) &
        &   * ( y_ExtLatN(jmax/2  ) - y_ExtLatN(jmax/2+3) ) )
      a_LCIFUVNP(jmax/2+1) = 1.0d100
      a_LCIFUVNP(jmax/2+2) =                                &
        &     ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2-1) ) &
        &   * ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2  ) ) &
        &   * ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2+3) ) &
        & / ( ( y_ExtLatN(jmax/2+2) - y_ExtLatN(jmax/2-1) ) &
        &   * ( y_ExtLatN(jmax/2+2) - y_ExtLatN(jmax/2  ) ) &
        &   * ( y_ExtLatN(jmax/2+2) - y_ExtLatN(jmax/2+3) ) )
      a_LCIFUVNP(jmax/2+3) =                                &
        &     ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2-1) ) &
        &   * ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2  ) ) &
        &   * ( y_ExtLatN(jmax/2+1) - y_ExtLatN(jmax/2+2) ) &
        & / ( ( y_ExtLatN(jmax/2+3) - y_ExtLatN(jmax/2-1) ) &
        &   * ( y_ExtLatN(jmax/2+3) - y_ExtLatN(jmax/2  ) ) &
        &   * ( y_ExtLatN(jmax/2+3) - y_ExtLatN(jmax/2+2) ) )
    end if


    !
    ! variable for estimating polar values
    !
    if( myrank == ( nprocs-1 ) ) then
      SumSinSq = imax / 2
    else
      SumSinSq = 1.0d100
    end if


  end subroutine SLTTExtArrInit

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

  ! Under development (yot, 2013/01/27)
  subroutine SLTTExtArrMkJMAXTable

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

    !
    ! MPI
    !
    use mpi_wrapper, only : nprocs, myrank, &
      & MPIWrapperISend, MPIWrapperIRecv, MPIWrapperWait

    !
    ! local variables
    !
    integer          , save, allocatable :: a_TblJMAX(:)
    integer          , save, allocatable :: a_NSSepTblRank(:)
    integer          , save, allocatable :: a_NSSepTblNGrid(:)
    character(STRING), save, allocatable :: a_NSSepTblNS(:)

    integer              :: a_SendBuf(1)
    integer, allocatable :: aa_RecvBuf(:,:)
    integer, allocatable :: a_iReqSend(:)
    integer, allocatable :: a_iReqRecv(:)

    integer              :: n
    integer              :: l


    allocate( a_TblJMAX(0:nprocs-1) )

    ! Make a table containing jmax in all processes
    !
    allocate( aa_RecvBuf(1,0:nprocs-1) )
    allocate( a_iReqSend(0:nprocs-1) )
    allocate( a_iReqRecv(0:nprocs-1) )
    !
    a_SendBuf = jmax
    do n = 0, nprocs-1
      if ( n == myrank ) cycle
      call MPIWrapperISend( n, 1, a_SendBuf      , a_iReqSend(n) )
      call MPIWrapperIRecv( n, 1, aa_RecvBuf(:,n), a_iReqRecv(n) )
    end do
    do n = 0, nprocs-1
      if ( n == myrank ) cycle
      call MPIWrapperWait( a_iReqSend(n) )
      call MPIWrapperWait( a_iReqRecv(n) )
    end do
    !
    aa_RecvBuf(:,myrank) = a_SendBuf
    do n = 0, nprocs-1
      a_TblJMAX(n) = aa_RecvBuf(1,n)
    end do
    !
    deallocate( aa_RecvBuf )
    deallocate( a_iReqSend )
    deallocate( a_iReqRecv )


    ! Table is checked
    do n = 0, nprocs-1
      if ( mod( a_TblJMAX(n), 2 ) /= 0 ) then
        call MessageNotify( 'E', module_name, 'Unexpected jmax in process %d, %d.', i = (/ n, a_TblJMAX(n) /) )
      end if
    end do


    allocate( a_NSSepTblRank (1:nprocs*2) )
    allocate( a_NSSepTblNGrid(1:nprocs*2) )
    allocate( a_NSSepTblNS   (1:nprocs*2) )

    ! a_NSSepTblRank : rank included in a North-South separate table
    l = 1
    do n = nprocs-1, 0, -1
      a_NSSepTblRank(l) = n
      l = l + 1
    end do
    do n = 0, nprocs-1
      a_NSSepTblRank(l) = n
      l = l + 1
    end do
    ! a_NSSepTblNGrid : number of grid included in a North-South separate table
    do l = 1, nprocs*2
      a_NSSepTblNGrid(l) = a_TblJMAX(a_NSSepTblRank(l)) / 2
    end do
    ! a_NSSepTblNS : Symbol representing south- or north-array
    !              : included in a North-South separate table
    do l = 1, nprocs*2
      if ( l <= nprocs ) then
        a_NSSepTblNS(l) = 'S'
      else
        a_NSSepTblNS(l) = 'N'
      end if
    end do

!!$    a_NSSepTblNToBeSent
!!$    a_NSSepTblNToBeRecv

    deallocate( a_NSSepTblRank )
    deallocate( a_NSSepTblNGrid )
    deallocate( a_TblJMAX )

  end subroutine SLTTExtArrMkJMAXTable

  !**************************************************************************
  ! Making new array(3+jm+3,1+km+2) from array(jm,km). 
  ! The array(3+jm+3,1+km+2) includes values at opposite hemisphere 
  ! in both latitudinal and longitudinal direction. 
  ! This routine uses int1d subroutine in interpolate2d.f file 
  ! to estimate values at poles, and fft3 subroutine in fft3.f file 
  ! to constrain polar values. 
  !**************************************************************************
  ! 01/09/08 :
  ! 01/09/09 : The size of array is changed from (2+jm+3,1+km+2) 
  !          : to (3+jm+3,1+km+2). 
  ! 01/09/10 : The name of subroutine which interpolates in one-dimensional 
  !          : space is changed from inte1d to int1d. 
  ! 01/09/10 : The on/off switch of shape-preserving constraint is added 
  !          : to int1d function. 
  ! 01/09/10 : The on/off switch of shape-preserving constraint is added. 
  !**************************************************************************
  ! jm,km
  ! lat
  ! lon
  ! arr
  ! ex_lat
  ! ex_lon
  ! ex_arr
  ! sw      : sw=0,  Shape-preserving constraint is NOT applied. 
  !         : sw!=0, Shape-preserving constraint is applied. 
  !**************************************************************************

  subroutine SLTTExtArrExt(                       &
    & y_ExtLatS, y_ExtLatN,                       & ! (in)
    & x_SinLonS, x_CosLonS, x_SinLonN, x_CosLonN, & ! (in)
    & xyzf_QMix, xyz_U, xyz_V,                    & ! (in)
    & xyzf_ExtDQMixDLatS, xyzf_ExtDQMixDLatN,     & ! (in)
    & xyzf_ExtQMixS, xyzf_ExtQMixN,               & ! (out)
    & xyz_ExtUS, xyz_ExtUN,                       & ! (out)
    & xyz_ExtVS, xyz_ExtVN                        & ! (out)
    & )

    use mpi_wrapper, only : nprocs, myrank
!!$    use mpi_wrapper, only : nprocs, myrank, &
!!$      & MPIWrapperISend, MPIWrapperIRecv, MPIWrapperWait
    use sltt_const , only : &
      & dtjw, iexmin, iexmax, jexmins, jexmaxs, jexminn, jexmaxn
    use sltt_lagint, only : SLTTIrrHerIntQui1DNonUni

    real(DP), intent(in ) :: y_ExtLatS(jexmins:jexmaxs)
    real(DP), intent(in ) :: y_ExtLatN(jexminn:jexmaxn)
    real(DP), intent(in ) :: x_SinLonS(0:imax-1)
    real(DP), intent(in ) :: x_CosLonS(0:imax-1)
    real(DP), intent(in ) :: x_SinLonN(0:imax-1)
    real(DP), intent(in ) :: x_CosLonN(0:imax-1)
    real(DP), intent(in ) :: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    real(DP), intent(in ) :: xyz_U    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_V    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyzf_ExtDQMixDLatS(iexmin:iexmax, jexmins:jexmaxs, 1:kmax, 1:ncmax)
    real(DP), intent(in ) :: xyzf_ExtDQMixDLatN(iexmin:iexmax, jexminn:jexmaxn, 1:kmax, 1:ncmax)
    real(DP), intent(out) :: xyzf_ExtQMixS(iexmin:iexmax, jexmins:jexmaxs, 1:kmax, 1:ncmax)
    real(DP), intent(out) :: xyzf_ExtQMixN(iexmin:iexmax, jexminn:jexmaxn, 1:kmax, 1:ncmax)
    real(DP), intent(out) :: xyz_ExtUS    (iexmin:iexmax, jexmins:jexmaxs, 1:kmax)
    real(DP), intent(out) :: xyz_ExtUN    (iexmin:iexmax, jexminn:jexmaxn, 1:kmax)
    real(DP), intent(out) :: xyz_ExtVS    (iexmin:iexmax, jexmins:jexmaxs, 1:kmax)
    real(DP), intent(out) :: xyz_ExtVN    (iexmin:iexmax, jexminn:jexmaxn, 1:kmax)


      !
      ! local variables
      !
!!$      !
!!$      ! variables for data transfer using MPI
!!$      !
!!$      real(DP)                :: sendbuf( imax, dtjw, kmax, ntrc+2 )
!!$      real(DP)                :: recvbuf( imax, dtjw, kmax, ntrc+2 )
!!$      integer            :: irec_send, irec_recv
!!$      integer            :: istatus( mpi_status_size )
!!$      integer            :: ierr

    !
    ! variables for estimation of values at poles
    !
    real(DP)                :: Ave
    real(DP)                :: SumC
    real(DP)                :: SumS
!!$    real(DP)                :: &
!!$      hcpif_q_dataz( imax*kmax, 6 ), &
!!$      d1_1d        ( imax*kmax )   , d2_1d( imax*kmax )
!!$    integer            :: mm
!!$
!!$    integer            :: i, j, k, m, n, nt
!!$    integer            :: ii
!!$
!!$
!!$    imaxh = imax / 2
!!$    mm    = imax * kmax


      !
      ! initialization for debug
      !
!!$      ex_gq( :, :, :, : ) = 1.0d100




    integer  :: idest
    integer  :: idep
!!$    integer  :: a_ireq_send(4)
!!$    integer  :: a_ireq_recv(4)

    ! SS : Southern hemisphere, Southward array
    ! SN : Southern hemisphere, Northward array
    ! NS : Northern hemisphere, Southward array
    ! NN : Northern hemisphere, Northward array
    real(DP), allocatable :: xyz_USN       (:,:,:)
    real(DP), allocatable :: xyz_UNS       (:,:,:)
    real(DP), allocatable :: xyz_USS       (:,:,:)
    real(DP), allocatable :: xyz_UNN       (:,:,:)
    real(DP), allocatable :: xyz_VSN       (:,:,:)
    real(DP), allocatable :: xyz_VNS       (:,:,:)
    real(DP), allocatable :: xyz_VSS       (:,:,:)
    real(DP), allocatable :: xyz_VNN       (:,:,:)
    real(DP), allocatable :: xyzf_QMixSN   (:,:,:,:)
    real(DP), allocatable :: xyzf_QMixNS   (:,:,:,:)
    real(DP), allocatable :: xyzf_QMixSS   (:,:,:,:)
    real(DP), allocatable :: xyzf_QMixNN   (:,:,:,:)

    real(DP) :: h, theta, thetasq

    integer  :: i
    integer  :: j
    integer  :: k
    integer  :: n
    integer  :: ii


    !====================================================================

!!$    if ( myrank > 0          ) then
      allocate( xyz_USN       (0:imax-1, 1:dtjw, 1:kmax) )
      allocate( xyz_UNS       (0:imax-1, 1:dtjw, 1:kmax) )
      allocate( xyz_VSN       (0:imax-1, 1:dtjw, 1:kmax) )
      allocate( xyz_VNS       (0:imax-1, 1:dtjw, 1:kmax) )
      allocate( xyzf_QMixSN   (0:imax-1, 1:dtjw, 1:kmax, ncmax) )
      allocate( xyzf_QMixNS   (0:imax-1, 1:dtjw, 1:kmax, ncmax) )
!!$    end if
!!$    if ( myrank < (nprocs-1) ) then
      allocate( xyz_USS       (0:imax-1, 1:dtjw, 1:kmax) )
      allocate( xyz_UNN       (0:imax-1, 1:dtjw, 1:kmax) )
      allocate( xyz_VSS       (0:imax-1, 1:dtjw, 1:kmax) )
      allocate( xyz_VNN       (0:imax-1, 1:dtjw, 1:kmax) )
      allocate( xyzf_QMixSS   (0:imax-1, 1:dtjw, 1:kmax, ncmax) )
      allocate( xyzf_QMixNN   (0:imax-1, 1:dtjw, 1:kmax, ncmax) )
!!$    end if


    call SLTTExtArrExtMPICore(                                   &
      & xyzf_QMix,                                                      & ! (in)
      & xyzf_QMixSN   , xyzf_QMixSS   , xyzf_QMixNS   , xyzf_QMixNN   , & ! (out)
      & xyz_U, xyz_V,                                                   & ! (in)
      & xyz_USN   , xyz_USS   , xyz_UNS   , xyz_UNN   ,                 & ! (out)
      & xyz_VSN   , xyz_VSS   , xyz_VNS   , xyz_VNN                     & ! (out)
      & )



    do k = 1, kmax
      do j = 1, jmax/2
        do i = 0, imax-1
          xyz_ExtUS(i,j,k) = xyz_U(i,j       ,k)
          xyz_ExtUN(i,j,k) = xyz_U(i,j+jmax/2,k)
          xyz_ExtVS(i,j,k) = xyz_V(i,j       ,k)
          xyz_ExtVN(i,j,k) = xyz_V(i,j+jmax/2,k)
        end do
      end do
    end do
    do n = 1, ncmax
      do k = 1, kmax
        do j = 1, jmax/2
          do i = 0, imax-1
            xyzf_ExtQMixS(i,j,k,n) = xyzf_QMix(i,j       ,k,n)
            xyzf_ExtQMixN(i,j,k,n) = xyzf_QMix(i,j+jmax/2,k,n)
          end do
        end do
      end do
    end do


    ! southern edge of southern array
    if( myrank == (nprocs-1) ) then
      ! values at south pole
      do k = 1, kmax
        do j = 1, dtjw
          do i = 0, imax/2-1
            ii = i + imax/2
            xyz_ExtUS(ii,0-j,k) = - xyz_ExtUS(i,j,k)
            xyz_ExtVS(ii,0-j,k) = - xyz_ExtVS(i,j,k)
          end do
          do i = imax/2, imax-1
            ii = i - imax/2
            xyz_ExtUS(ii,0-j,k) = - xyz_ExtUS(i,j,k)
            xyz_ExtVS(ii,0-j,k) = - xyz_ExtVS(i,j,k)
          end do
        end do
      end do
      do k = 1, kmax
        do i = 0, imax-1
          ! quadratic interpolation (old)
!!$          xyz_ExtUS(i,0,k) =                       &
!!$            &   a_LQIFUVSP(-1) * xyz_ExtUS(i,-1,k) &
!!$            & + a_LQIFUVSP( 1) * xyz_ExtUS(i, 1,k) &
!!$            & + a_LQIFUVSP( 2) * xyz_ExtUS(i, 2,k)
!!$          xyz_ExtVS(i,0,k) =                       &
!!$            &   a_LQIFUVSP(-1) * xyz_ExtVS(i,-1,k) &
!!$            & + a_LQIFUVSP( 1) * xyz_ExtVS(i, 1,k) &
!!$            & + a_LQIFUVSP( 2) * xyz_ExtVS(i, 2,k)
          ! cubic interpolation
          xyz_ExtUS(i,0,k) =                       &
            &   a_LCIFUVSP(-2) * xyz_ExtUS(i,-2,k) &
            & + a_LCIFUVSP(-1) * xyz_ExtUS(i,-1,k) &
            & + a_LCIFUVSP( 1) * xyz_ExtUS(i, 1,k) &
            & + a_LCIFUVSP( 2) * xyz_ExtUS(i, 2,k)
          xyz_ExtVS(i,0,k) =                       &
            &   a_LCIFUVSP(-2) * xyz_ExtVS(i,-2,k) &
            & + a_LCIFUVSP(-1) * xyz_ExtVS(i,-1,k) &
            & + a_LCIFUVSP( 1) * xyz_ExtVS(i, 1,k) &
            & + a_LCIFUVSP( 2) * xyz_ExtVS(i, 2,k)
        end do
      end do
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax/2-1
              ii = i + imax/2
              xyzf_ExtQMixS(ii,0-j,k,n) = xyzf_ExtQMixS(i,j,k,n)
            end do
            do i = imax/2, imax-1
              ii = i - imax/2
              xyzf_ExtQMixS(ii,0-j,k,n) = xyzf_ExtQMixS(i,j,k,n)
            end do
          end do
        end do
        do k = 1, kmax
          do i = 0, imax-1
          ! quadratic interpolation (old)
!!$            xyzf_ExtQMixS(i,0,k,n) =                         &
!!$              &   a_LQIFUVSP(-1) * xyzf_ExtQMixS(i,-1,k,n) &
!!$              & + a_LQIFUVSP( 1) * xyzf_ExtQMixS(i, 1,k,n) &
!!$              & + a_LQIFUVSP( 2) * xyzf_ExtQMixS(i, 2,k,n)
!!$          ! cubic interpolation
!!$            xyzf_ExtQMixS(i,0,k,n) =                       &
!!$              &   a_LCIFUVSP(-2) * xyzf_ExtQMixS(i,-2,k,n) &
!!$              & + a_LCIFUVSP(-1) * xyzf_ExtQMixS(i,-1,k,n) &
!!$              & + a_LCIFUVSP( 1) * xyzf_ExtQMixS(i, 1,k,n) &
!!$              & + a_LCIFUVSP( 2) * xyzf_ExtQMixS(i, 2,k,n)
!!$
            ! Polar value estimated by 1D Hermite Quintic Interpolation
            xyzf_ExtQMixS(i,0,k,n) =            &
              & SLTTIrrHerIntQui1DNonUni        &
              & (                               &
              &   xyzf_ExtQMixS(i,-2,k,n),      &
              &   xyzf_ExtQMixS(i,-1,k,n),      &
              &   xyzf_ExtQMixS(i, 1,k,n),      &
              &   xyzf_ExtQMixS(i, 2,k,n),      &
              &   xyzf_ExtDQMixDLatS(i,-1,k,n), &
              &   xyzf_ExtDQMixDLatS(i, 1,k,n), &
              &   y_ExtLatS(-2)-y_ExtLatS(-1),  &
              &   y_ExtLatS( 1)-y_ExtLatS(-1),  &
              &   y_ExtLatS( 2)-y_ExtLatS(-1),  &
              &   y_ExtLatS( 0)-y_ExtLatS(-1)   &
              & )

          end do
        end do
      end do

      ! only wavenumber 1 component is retained for zonal and meridional 
      ! wind velocities
      j = 0
      do k = 1, kmax
        SumC = 0.0_DP
        SumS = 0.0_DP
        do i = 0, imax-1
          SumC = SumC + xyz_ExtUS(i,j,k) * x_CosLonS(i)
          SumS = SumS + xyz_ExtUS(i,j,k) * x_SinLonS(i)
        end do
        SumC = SumC / SumSinSq
        SumS = SumS / SumSinSq
        do i = 0, imax-1
          xyz_ExtUS(i,j,k) = SumC * x_CosLonS(i) + SumS * x_SinLonS(i)
        end do
      end do
      do k = 1, kmax
        SumC = 0.0_DP
        SumS = 0.0_DP
        do i = 0, imax-1
          SumC = SumC + xyz_ExtVS(i,j,k) * x_CosLonS(i)
          SumS = SumS + xyz_ExtVS(i,j,k) * x_SinLonS(i)
        end do
        SumC = SumC / SumSinSq
        SumS = SumS / SumSinSq
        do i = 0, imax-1
          xyz_ExtVS(i,j,k) = SumC * x_CosLonS(i) + SumS * x_SinLonS(i)
        end do
      end do
      ! zonal average is set for mixing ratio
      j = 0
      do n = 1, ncmax
        do k = 1, kmax
          Ave = 0.0_DP
          do i = 0, imax-1
            Ave = Ave + xyzf_ExtQMixS(i,j,k,n)
          end do
          Ave = Ave / dble( imax )
          do i = 0, imax-1
            xyzf_ExtQMixS(i,j,k,n) = Ave
          end do
        end do
      end do

    else

      do k = 1, kmax
        do j = 1, dtjw
          do i = 0, imax-1
            xyz_ExtUS(i,1-j,k) = xyz_USS   (i,dtjw-(j-1),k)
            xyz_ExtVS(i,1-j,k) = xyz_VSS   (i,dtjw-(j-1),k)
          end do
        end do
      end do
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax-1
              xyzf_ExtQMixS(i,1-j,k,n) = xyzf_QMixSS   (i,dtjw-(j-1),k,n)
            end do
          end do
        end do
      end do
    end if
    ! northern edge of southern array
    if ( myrank == 0 ) then
      do k = 1, kmax
        do j = 1, dtjw
          do i = 0, imax-1
            xyz_ExtUS(i,jmax/2+j,k) = xyz_ExtUN(i,j,k)
            xyz_ExtVS(i,jmax/2+j,k) = xyz_ExtVN(i,j,k)
          end do
        end do
      end do
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax-1
              xyzf_ExtQMixS(i,jmax/2+j,k,n) = xyzf_ExtQMixN(i,j,k,n)
            end do
          end do
        end do
      end do
    else
      do k = 1, kmax
        do j = 1, dtjw
          do i = 0, imax-1
            xyz_ExtUS(i,jmax/2+j,k) = xyz_USN   (i,j,k)
            xyz_ExtVS(i,jmax/2+j,k) = xyz_VSN   (i,j,k)
          end do
        end do
      end do
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax-1
              xyzf_ExtQMixS(i,jmax/2+j,k,n) = xyzf_QMixSN   (i,j,k,n)
            end do
          end do
        end do
      end do
    end if
    !
    ! southern edge of northern array
    if ( myrank == 0 ) then
      do k = 1, kmax
        do j = 1, dtjw
          do i = 0, imax-1
            xyz_ExtUN(i,1-j,k) = xyz_ExtUS(i,jmax/2-(j-1),k)
            xyz_ExtVN(i,1-j,k) = xyz_ExtVS(i,jmax/2-(j-1),k)
          end do
        end do
      end do
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax-1
              xyzf_ExtQMixN(i,1-j,k,n) = xyzf_ExtQMixS(i,jmax/2-(j-1),k,n)
            end do
          end do
        end do
      end do
    else
      do k = 1, kmax
        do j = 1, dtjw
          do i = 0, imax-1
            xyz_ExtUN(i,1-j,k) = xyz_UNS   (i,dtjw-(j-1),k)
            xyz_ExtVN(i,1-j,k) = xyz_VNS   (i,dtjw-(j-1),k)
          end do
        end do
      end do
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax-1
              xyzf_ExtQMixN(i,1-j,k,n) = xyzf_QMixNS   (i,dtjw-(j-1),k,n)
            end do
          end do
        end do
      end do
    end if
    ! northern edge of northern array
    if( myrank == (nprocs-1) ) then
      ! values at north pole

      do k = 1, kmax
        do j = 1, dtjw
          do i = 0, imax/2-1
            ii = i + imax/2
            xyz_ExtUN(ii,jmax/2+1+j,k) = - xyz_ExtUN(i,jmax/2+1-j,k)
            xyz_ExtVN(ii,jmax/2+1+j,k) = - xyz_ExtVN(i,jmax/2+1-j,k)
          end do
          do i = imax/2, imax-1
            ii = i - imax/2
            xyz_ExtUN(ii,jmax/2+1+j,k) = - xyz_ExtUN(i,jmax/2+1-j,k)
            xyz_ExtVN(ii,jmax/2+1+j,k) = - xyz_ExtVN(i,jmax/2+1-j,k)
          end do
        end do
      end do
      do k = 1, kmax
        do i = 0, imax-1
          ! quadratic interpolation (old)
!!$          xyz_ExtUN(i,jmax/2+1,k) =                              &
!!$            &   a_LQIFUVNP(jmax/2-1) * xyz_ExtUN(i,jmax/2-1,k) &
!!$            & + a_LQIFUVNP(jmax/2  ) * xyz_ExtUN(i,jmax/2  ,k) &
!!$            & + a_LQIFUVNP(jmax/2+2) * xyz_ExtUN(i,jmax/2+2,k)
!!$          xyz_ExtVN(i,jmax/2+1,k) =                              &
!!$            &   a_LQIFUVNP(jmax/2-1) * xyz_ExtVN(i,jmax/2-1,k) &
!!$            & + a_LQIFUVNP(jmax/2  ) * xyz_ExtVN(i,jmax/2  ,k) &
!!$            & + a_LQIFUVNP(jmax/2+2) * xyz_ExtVN(i,jmax/2+2,k)
          ! qcubic interpolation
          xyz_ExtUN(i,jmax/2+1,k) =                            &
            &   a_LCIFUVNP(jmax/2-1) * xyz_ExtUN(i,jmax/2-1,k) &
            & + a_LCIFUVNP(jmax/2  ) * xyz_ExtUN(i,jmax/2  ,k) &
            & + a_LCIFUVNP(jmax/2+2) * xyz_ExtUN(i,jmax/2+2,k) &
            & + a_LCIFUVNP(jmax/2+3) * xyz_ExtUN(i,jmax/2+3,k)
          xyz_ExtVN(i,jmax/2+1,k) =                            &
            &   a_LCIFUVNP(jmax/2-1) * xyz_ExtVN(i,jmax/2-1,k) &
            & + a_LCIFUVNP(jmax/2  ) * xyz_ExtVN(i,jmax/2  ,k) &
            & + a_LCIFUVNP(jmax/2+2) * xyz_ExtVN(i,jmax/2+2,k) &
            & + a_LCIFUVNP(jmax/2+3) * xyz_ExtVN(i,jmax/2+3,k)
        end do
      end do
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax/2-1
              ii = i + imax/2
              xyzf_ExtQMixN(ii,jmax/2+1+j,k,n) = xyzf_ExtQMixN(i,jmax/2+1-j,k,n)
            end do
            do i = imax/2, imax-1
              ii = i - imax/2
              xyzf_ExtQMixN(ii,jmax/2+1+j,k,n) = xyzf_ExtQMixN(i,jmax/2+1-j,k,n)
            end do
          end do
        end do
        do k = 1, kmax
          do i = 0, imax-1
!!$          ! quadratic interpolation (old)
!!$            xyzf_ExtQMixN(i,jmax/2+1,k,n) =                              &
!!$              &   a_LQIFUVNP(jmax/2-1) * xyzf_ExtQMixN(i,jmax/2-1,k,n) &
!!$              & + a_LQIFUVNP(jmax/2  ) * xyzf_ExtQMixN(i,jmax/2  ,k,n) &
!!$              & + a_LQIFUVNP(jmax/2+2) * xyzf_ExtQMixN(i,jmax/2+2,k,n)
!!$          ! cubic interpolation
!!$            xyzf_ExtQMixN(i,jmax/2+1,k,n) =                            &
!!$              &   a_LCIFUVNP(jmax/2-1) * xyzf_ExtQMixN(i,jmax/2-1,k,n) &
!!$              & + a_LCIFUVNP(jmax/2  ) * xyzf_ExtQMixN(i,jmax/2  ,k,n) &
!!$              & + a_LCIFUVNP(jmax/2+2) * xyzf_ExtQMixN(i,jmax/2+2,k,n) &
!!$              & + a_LCIFUVNP(jmax/2+3) * xyzf_ExtQMixN(i,jmax/2+3,k,n)
!!$
            ! Polar value estimated by 1D Hermite Quintic Interpolation
            xyzf_ExtQMixN(i,jmax/2+1,k,n) =                &
              & SLTTIrrHerIntQui1DNonUni                   &
              & (                                          &
              &   xyzf_ExtQMixN(i,jmax/2-1,k,n),           &
              &   xyzf_ExtQMixN(i,jmax/2  ,k,n),           &
              &   xyzf_ExtQMixN(i,jmax/2+2,k,n),           &
              &   xyzf_ExtQMixN(i,jmax/2+3,k,n),           &
              &   xyzf_ExtDQMixDLatN(i,jmax/2  ,k,n),      &
              &   xyzf_ExtDQMixDLatN(i,jmax/2+2,k,n),      &
              &   y_ExtLatN(jmax/2-1)-y_ExtLatN(jmax/2  ), &
              &   y_ExtLatN(jmax/2+2)-y_ExtLatN(jmax/2  ), &
              &   y_ExtLatN(jmax/2+3)-y_ExtLatN(jmax/2  ), &
              &   y_ExtLatN(jmax/2+1)-y_ExtLatN(jmax/2  )  &
              & )
          end do
        end do
      end do

      ! only wavenumber 1 component is retained for zonal and meridional 
      ! wind velocities
      j = jmax/2+1
      do k = 1, kmax
        SumC = 0.0_DP
        SumS = 0.0_DP
        do i = 0, imax-1
          SumC = SumC + xyz_ExtUN(i,j,k) * x_CosLonN(i)
          SumS = SumS + xyz_ExtUN(i,j,k) * x_SinLonN(i)
        end do
        SumC = SumC / SumSinSq
        SumS = SumS / SumSinSq
        do i = 0, imax-1
          xyz_ExtUN(i,j,k) = SumC * x_CosLonN(i) + SumS * x_SinLonN(i)
        end do
      end do
      do k = 1, kmax
        SumC = 0.0_DP
        SumS = 0.0_DP
        do i = 0, imax-1
          SumC = SumC + xyz_ExtVN(i,j,k) * x_CosLonN(i)
          SumS = SumS + xyz_ExtVN(i,j,k) * x_SinLonN(i)
        end do
        SumC = SumC / SumSinSq
        SumS = SumS / SumSinSq
        do i = 0, imax-1
          xyz_ExtVN(i,j,k) = SumC * x_CosLonN(i) + SumS * x_SinLonN(i)
        end do
      end do
      ! zonal average is set for mixing ratio
      j = jmax/2+1
      do n = 1, ncmax
        do k = 1, kmax
          Ave = 0.0_DP
          do i = 0, imax-1
            Ave = Ave + xyzf_ExtQMixN(i,j,k,n)
          end do
          Ave = Ave / dble( imax )
          do i = 0, imax-1
            xyzf_ExtQMixN(i,j,k,n) = Ave
          end do
        end do
      end do

    else

      do k = 1, kmax
        do j = 1, dtjw
          do i = 0, imax-1
            xyz_ExtUN(i,jmax/2+j,k) = xyz_UNN   (i,j,k)
            xyz_ExtVN(i,jmax/2+j,k) = xyz_VNN   (i,j,k)
          end do
        end do
      end do
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax-1
              xyzf_ExtQMixN(i,jmax/2+j,k,n) = xyzf_QMixNN   (i,j,k,n)
            end do
          end do
        end do
      end do
    end if

!!$    if ( myrank > 0          ) then
      deallocate( xyz_USN        )
      deallocate( xyz_UNS        )
      deallocate( xyz_VSN        )
      deallocate( xyz_VNS        )
      deallocate( xyzf_QMixSN    )
      deallocate( xyzf_QMixNS    )
!!$    end if
!!$    if ( myrank < (nprocs-1) ) then
      deallocate( xyz_USS        )
      deallocate( xyz_UNN        )
      deallocate( xyz_VSS        )
      deallocate( xyz_VNN        )
      deallocate( xyzf_QMixSS    )
      deallocate( xyzf_QMixNN    )
!!$    end if



    !===========================================
    ! set values at longitudinal edge
    !-------------------------------------------
    do k = 1, kmax
      do j = jexmins, jexmaxs
        do i = iexmin, 0-1
          xyz_ExtUS(i,j,k) = xyz_ExtUS(imax+i,j,k)
          xyz_ExtVS(i,j,k) = xyz_ExtVS(imax+i,j,k)
        end do
        do i = imax-1+1, iexmax
          xyz_ExtUS(i,j,k) = xyz_ExtUS(i-imax,j,k)
          xyz_ExtVS(i,j,k) = xyz_ExtVS(i-imax,j,k)
        end do
      end do
      do j = jexminn, jexmaxn
        do i = iexmin, 0-1
          xyz_ExtUN(i,j,k) = xyz_ExtUN(imax+i,j,k)
          xyz_ExtVN(i,j,k) = xyz_ExtVN(imax+i,j,k)
        end do
        do i = imax-1+1, iexmax
          xyz_ExtUN(i,j,k) = xyz_ExtUN(i-imax,j,k)
          xyz_ExtVN(i,j,k) = xyz_ExtVN(i-imax,j,k)
        end do
      end do
    end do
    do n = 1, ncmax
      do k = 1, kmax
        do j = jexmins, jexmaxs
          do i = iexmin, 0-1
            xyzf_ExtQMixS(i,j,k,n) = xyzf_ExtQMixS(imax+i,j,k,n)
          end do
          do i = imax-1+1, iexmax
            xyzf_ExtQMixS(i,j,k,n) = xyzf_ExtQMixS(i-imax,j,k,n)
          end do
        end do
        do j = jexminn, jexmaxn
          do i = iexmin, 0-1
            xyzf_ExtQMixN(i,j,k,n) = xyzf_ExtQMixN(imax+i,j,k,n)
          end do
          do i = imax-1+1, iexmax
            xyzf_ExtQMixN(i,j,k,n) = xyzf_ExtQMixN(i-imax,j,k,n)
          end do
        end do
      end do
    end do


  end subroutine SLTTExtArrExt

  !--------------------------------------------------------------------------------------
  subroutine SLTTExtArrExt2(                      &
    & x_SinLonS, x_CosLonS, x_SinLonN, x_CosLonN, & ! (in)
    & xyzf_QMix, pm,                              & ! (in)
    & xyzf_ExtQMixS, xyzf_ExtQMixN,               & ! (out)
    & PoleMethod                                  & ! (in)
    & )

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

    use mpi_wrapper, only : nprocs, myrank, &
      & MPIWrapperISend, MPIWrapperIRecv, MPIWrapperWait
    use sltt_const , only : dtjw, iexmin, iexmax, jexmins, jexmaxs, jexminn, jexmaxn

    real(DP), intent(in ) :: x_SinLonS(0:imax-1)
    real(DP), intent(in ) :: x_CosLonS(0:imax-1)
    real(DP), intent(in ) :: x_SinLonN(0:imax-1)
    real(DP), intent(in ) :: x_CosLonN(0:imax-1)
    real(DP), intent(in ) :: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    real(DP), intent(in ) :: pm   ! 配列拡張する際、極ごえ後に符号が変わる場合は -1.0を与える。そうでない場合は1.0を与える。
    real(DP), intent(out) :: xyzf_ExtQMixS(iexmin:iexmax, jexmins:jexmaxs, 1:kmax, 1:ncmax)
    real(DP), intent(out) :: xyzf_ExtQMixN(iexmin:iexmax, jexminn:jexmaxn, 1:kmax, 1:ncmax)
    character(*), intent(in ) :: PoleMethod
                                   ! "Mean"  : Longitudinal mean
                                   ! "Wave1" : Only wave #1 is retained.

    !
    ! local variables
    !
    !
    ! variables for estimation of values at poles
    !
    real(DP)                :: Ave
    real(DP)                :: SumC
    real(DP)                :: SumS


!!$    integer  :: idest
!!$    integer  :: idep
!!$    integer  :: a_ireq_send(4)
!!$    integer  :: a_ireq_recv(4)
!!$    real(DP), allocatable :: xyza_SendBuf(:,:,:,:,:)
!!$    real(DP), allocatable :: xyza_RecvBuf(:,:,:,:,:)
    real(DP), allocatable :: xyzf_QMixSN   (:,:,:,:)
    real(DP), allocatable :: xyzf_QMixNS   (:,:,:,:)
    real(DP), allocatable :: xyzf_QMixSS   (:,:,:,:)
    real(DP), allocatable :: xyzf_QMixNN   (:,:,:,:)

    real(DP) :: h, theta, thetasq

    integer  :: i
    integer  :: j
    integer  :: k
    integer  :: n
    integer  :: ii


    !====================================================================

!!$    allocate( xyza_SendBuf(0:imax-1, 1:dtjw, 1:kmax, ncmax, 4) )
!!$    allocate( xyza_RecvBuf(0:imax-1, 1:dtjw, 1:kmax, ncmax, 4) )
!!$    if ( myrank > 0          ) then
      allocate( xyzf_QMixSN   (0:imax-1, 1:dtjw, 1:kmax, ncmax) )
      allocate( xyzf_QMixNS   (0:imax-1, 1:dtjw, 1:kmax, ncmax) )
!!$    end if
!!$    if ( myrank < (nprocs-1) ) then
      allocate( xyzf_QMixSS   (0:imax-1, 1:dtjw, 1:kmax, ncmax) )
      allocate( xyzf_QMixNN   (0:imax-1, 1:dtjw, 1:kmax, ncmax) )
!!$    end if



    call SLTTExtArrExtMPICore(                                    &
      & xyzf_QMix,                                                      & ! (in)
      & xyzf_QMixSN   , xyzf_QMixSS   , xyzf_QMixNS   , xyzf_QMixNN     & ! (out)
    & )


    do n = 1, ncmax
      do k = 1, kmax
        do j = 1, jmax/2
          do i = 0, imax-1
            xyzf_ExtQMixS(i,j,k,n) = xyzf_QMix(i,j       ,k,n)
            xyzf_ExtQMixN(i,j,k,n) = xyzf_QMix(i,j+jmax/2,k,n)
          end do
        end do
      end do
    end do


    ! southern edge of southern array
    if( myrank == (nprocs-1) ) then
      ! values at south pole

      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax/2-1
              ii = i + imax/2
              xyzf_ExtQMixS(ii,0-j,k,n) = pm * xyzf_ExtQMixS(i,j,k,n)
            end do
            do i = imax/2, imax-1
              ii = i - imax/2
              xyzf_ExtQMixS(ii,0-j,k,n) = pm * xyzf_ExtQMixS(i,j,k,n)
            end do
          end do
        end do
        do k = 1, kmax
          do i = 0, imax-1
!!$            xyzf_ExtQMixS(i,0,k,n) =                       & !南極点の値
!!$              &   a_LQIFUVSP(-1) * xyzf_ExtQMixS(i,-1,k,n) &
!!$              & + a_LQIFUVSP( 1) * xyzf_ExtQMixS(i, 1,k,n) &
!!$              & + a_LQIFUVSP( 2) * xyzf_ExtQMixS(i, 2,k,n)
            xyzf_ExtQMixS(i,0,k,n) =                       & !南極点の値
              &   a_LCIFUVSP(-2) * xyzf_ExtQMixS(i,-2,k,n) &
              & + a_LCIFUVSP(-1) * xyzf_ExtQMixS(i,-1,k,n) &
              & + a_LCIFUVSP( 1) * xyzf_ExtQMixS(i, 1,k,n) &
              & + a_LCIFUVSP( 2) * xyzf_ExtQMixS(i, 2,k,n)
          end do
        end do
      end do

      select case ( PoleMethod )
      case ( "Mean" )
        ! Longitudinal mean
        j = 0
        do n = 1, ncmax
          do k = 1, kmax
            Ave = 0.0_DP
            do i = 0, imax-1
              Ave = Ave + xyzf_ExtQMixS(i,j,k,n)
            end do
            Ave = Ave / dble( imax )
            do i = 0, imax-1
              xyzf_ExtQMixS(i,j,k,n) = Ave !!南極点の値を各iで統一する。
            end do
          end do
        end do
      case ( "Wave1" )
        ! Only wave #1 is retained.
        j = 0
        do n = 1, ncmax
          do k = 1, kmax
            SumC = 0.0_DP
            SumS = 0.0_DP
            do i = 0, imax-1
              SumC = SumC + xyzf_ExtQMixS(i,j,k,n) * x_CosLonS(i)
              SumS = SumS + xyzf_ExtQMixS(i,j,k,n) * x_SinLonS(i)
            end do
            SumC = SumC / SumSinSq
            SumS = SumS / SumSinSq
            do i = 0, imax-1
              xyzf_ExtQMixS(i,j,k,n) = SumC * x_CosLonS(i) + SumS * x_SinLonS(i)
            end do
          end do
        end do
      case default
        call MessageNotify( 'E', module_name, 'PoleMethod of %c is inappropriate.', c1 = trim(PoleMethod) )
      end select

    else
!!$      do j = 1, jew
!!$        y_ExtLatS(1-j) = y_LatSS   (jew-(j-1))
!!$      end do
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax-1
              xyzf_ExtQMixS(i,1-j,k,n) = xyzf_QMixSS   (i,dtjw-(j-1),k,n)
            end do
          end do
        end do
      end do
    end if
    ! northern edge of southern array
    if ( myrank == 0 ) then
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax-1
              xyzf_ExtQMixS(i,jmax/2+j,k,n) = xyzf_ExtQMixN(i,j,k,n)
            end do
          end do
        end do
      end do
    else
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax-1
              xyzf_ExtQMixS(i,jmax/2+j,k,n) = xyzf_QMixSN   (i,j,k,n)
            end do
          end do
        end do
      end do
    end if
    !
    ! southern edge of northern array
    if ( myrank == 0 ) then
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax-1
              xyzf_ExtQMixN(i,-j+1,k,n) = xyzf_ExtQMixS(i,jmax/2-(j-1),k,n)
            end do
          end do
        end do
      end do
    else
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax-1
              xyzf_ExtQMixN(i,-j+1,k,n) = xyzf_QMixNS   (i,dtjw-(j-1),k,n)
            end do
          end do
        end do
      end do
    end if

    ! northern edge of northern array
    if( myrank == (nprocs-1) ) then
      ! values at north pole

      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax/2-1
              ii = i + imax/2
              xyzf_ExtQMixN(ii,jmax/2+1+j,k,n) = pm * xyzf_ExtQMixN(i,jmax/2+1-j,k,n)
            end do
            do i = imax/2, imax-1
              ii = i - imax/2
              xyzf_ExtQMixN(ii,jmax/2+1+j,k,n) = pm * xyzf_ExtQMixN(i,jmax/2+1-j,k,n)
            end do
          end do
        end do
        do k = 1, kmax
          do i = 0, imax-1
!!$            xyzf_ExtQMixN(i,jmax/2+1,k,n) =                            & !北極点での値
!!$              &   a_LQIFUVNP(jmax/2-1) * xyzf_ExtQMixN(i,jmax/2-1,k,n) &
!!$              & + a_LQIFUVNP(jmax/2  ) * xyzf_ExtQMixN(i,jmax/2  ,k,n) &
!!$              & + a_LQIFUVNP(jmax/2+2) * xyzf_ExtQMixN(i,jmax/2+2,k,n)
            xyzf_ExtQMixN(i,jmax/2+1,k,n) =                            & !北極点での値
              &   a_LCIFUVNP(jmax/2-1) * xyzf_ExtQMixN(i,jmax/2-1,k,n) &
              & + a_LCIFUVNP(jmax/2  ) * xyzf_ExtQMixN(i,jmax/2  ,k,n) &
              & + a_LCIFUVNP(jmax/2+2) * xyzf_ExtQMixN(i,jmax/2+2,k,n) &
              & + a_LCIFUVNP(jmax/2+3) * xyzf_ExtQMixN(i,jmax/2+3,k,n)
          end do
        end do
      end do

      select case ( PoleMethod )
      case ( "Mean" )
        ! Longitudinal mean
        j = jmax/2+1
        do n = 1, ncmax
          do k = 1, kmax
            Ave = 0.0_DP
            do i = 0, imax-1
              Ave = Ave + xyzf_ExtQMixN(i,j,k,n)
            end do
            Ave = Ave / dble( imax )
            do i = 0, imax-1
              xyzf_ExtQMixN(i,j,k,n) = Ave !値を各iで統一する。
            end do
          end do
        end do
      case ( "Wave1" )
        ! Only wave #1 is retained.
        j = jmax/2+1
        do n = 1, ncmax
          do k = 1, kmax
            SumC = 0.0_DP
            SumS = 0.0_DP
            do i = 0, imax-1
              SumC = SumC + xyzf_ExtQMixN(i,j,k,n) * x_CosLonN(i)
              SumS = SumS + xyzf_ExtQMixN(i,j,k,n) * x_SinLonN(i)
            end do
            SumC = SumC / SumSinSq
            SumS = SumS / SumSinSq
            do i = 0, imax-1
              xyzf_ExtQMixN(i,j,k,n) = SumC * x_CosLonN(i) + SumS * x_SinLonN(i)
            end do
          end do
        end do
      case default
        call MessageNotify( 'E', module_name, 'PoleMethod of %c is inappropriate.', c1 = trim(PoleMethod) )
      end select

    else
      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, dtjw
            do i = 0, imax-1
              xyzf_ExtQMixN(i,jmax/2+j,k,n) = xyzf_QMixNN   (i,j,k,n)
            end do
          end do
        end do
      end do
    end if

!!$    deallocate( xyza_SendBuf )
!!$    deallocate( xyza_RecvBuf )
!!$    if ( myrank > 0          ) then
      deallocate( xyzf_QMixSN    )
      deallocate( xyzf_QMixNS    )
!!$    end if
!!$    if ( myrank < (nprocs-1) ) then
      deallocate( xyzf_QMixSS    )
      deallocate( xyzf_QMixNN    )
!!$    end if



    !===========================================
    ! set values at longitudinal edge
    !-------------------------------------------
    do n = 1, ncmax
      do k = 1, kmax
        do j = jexmins, jexmaxs
          do i = iexmin, 0-1
            xyzf_ExtQMixS(i,j,k,n) = xyzf_ExtQMixS(imax+i,j,k,n)
          end do
          do i = imax-1+1, iexmax
            xyzf_ExtQMixS(i,j,k,n) = xyzf_ExtQMixS(i-imax,j,k,n)
          end do
        end do
        do j = jexminn, jexmaxn
          do i = iexmin, 0-1
            xyzf_ExtQMixN(i,j,k,n) = xyzf_ExtQMixN(imax+i,j,k,n)
          end do
          do i = imax-1+1, iexmax
            xyzf_ExtQMixN(i,j,k,n) = xyzf_ExtQMixN(i-imax,j,k,n)
          end do
        end do
      end do
    end do


  end subroutine SLTTExtArrExt2

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

  subroutine SLTTExtArrExtMPICore(                                    &
    & xyzf_QMix,                                                      & ! (in)
    & xyzf_QMixSN   , xyzf_QMixSS   , xyzf_QMixNS   , xyzf_QMixNN   , & ! (out)
    & xyz_U, xyz_V,                                                   & ! (in)
    & xyz_USN   , xyz_USS   , xyz_UNS   , xyz_UNN   ,                 & ! (out)
    & xyz_VSN   , xyz_VSS   , xyz_VNS   , xyz_VNN                     & ! (out)
    & )

    use dc_message, only: MessageNotify
    use mpi_wrapper, only : nprocs, myrank, &
      & MPIWrapperISend, MPIWrapperIRecv, MPIWrapperWait
    use sltt_const , only : &
      & dtjw, iexmin, iexmax, jexmins, jexmaxs, jexminn, jexmaxn

    real(DP), intent(in ) :: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    real(DP), intent(out) :: xyzf_QMixSN   (0:imax-1, 1:dtjw, 1:kmax, ncmax)
    real(DP), intent(out) :: xyzf_QMixSS   (0:imax-1, 1:dtjw, 1:kmax, ncmax)
    real(DP), intent(out) :: xyzf_QMixNS   (0:imax-1, 1:dtjw, 1:kmax, ncmax)
    real(DP), intent(out) :: xyzf_QMixNN   (0:imax-1, 1:dtjw, 1:kmax, ncmax)
    real(DP), intent(in ), optional :: xyz_U    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ), optional :: xyz_V    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out), optional :: xyz_USN       (0:imax-1, 1:dtjw, 1:kmax)
    real(DP), intent(out), optional :: xyz_USS       (0:imax-1, 1:dtjw, 1:kmax)
    real(DP), intent(out), optional :: xyz_UNS       (0:imax-1, 1:dtjw, 1:kmax)
    real(DP), intent(out), optional :: xyz_UNN       (0:imax-1, 1:dtjw, 1:kmax)
    real(DP), intent(out), optional :: xyz_VSN       (0:imax-1, 1:dtjw, 1:kmax)
    real(DP), intent(out), optional :: xyz_VSS       (0:imax-1, 1:dtjw, 1:kmax)
    real(DP), intent(out), optional :: xyz_VNN       (0:imax-1, 1:dtjw, 1:kmax)
    real(DP), intent(out), optional :: xyz_VNS       (0:imax-1, 1:dtjw, 1:kmax)


    !
    ! local variables
    !

    logical  :: FlagIncludeUV
    integer  :: narrsize
    integer  :: idest
    integer  :: idep
    integer  :: a_ireq_send_ss
    integer  :: a_ireq_send_sn
    integer  :: a_ireq_send_ns
    integer  :: a_ireq_send_nn
    integer  :: a_ireq_recv_ss
    integer  :: a_ireq_recv_sn
    integer  :: a_ireq_recv_ns
    integer  :: a_ireq_recv_nn
    real(DP), allocatable :: xyz_SendBufSS(:,:,:,:)
    real(DP), allocatable :: xyz_SendBufSN(:,:,:,:)
    real(DP), allocatable :: xyz_SendBufNS(:,:,:,:)
    real(DP), allocatable :: xyz_SendBufNN(:,:,:,:)
    real(DP), allocatable :: xyz_RecvBufSS(:,:,:,:)
    real(DP), allocatable :: xyz_RecvBufSN(:,:,:,:)
    real(DP), allocatable :: xyz_RecvBufNS(:,:,:,:)
    real(DP), allocatable :: xyz_RecvBufNN(:,:,:,:)


    if ( present( xyz_U ) .and. &
      &  present( xyz_V ) .and. &
      &  present( xyz_USN    ) .and. &
      &  present( xyz_USS    ) .and. &
      &  present( xyz_UNS    ) .and. &
      &  present( xyz_UNN    ) .and. &
      &  present( xyz_VSN    ) .and. &
      &  present( xyz_VSS    ) .and. &
      &  present( xyz_VNN    ) .and. &
      &  present( xyz_VNS    ) ) then
      FlagIncludeUV = .true.
    else
      if ( ( .not. present( xyz_U ) ) .and. &
        &  ( .not. present( xyz_V ) ) .and. &
        &  ( .not. present( xyz_USN    ) ) .and. &
        &  ( .not. present( xyz_USS    ) ) .and. &
        &  ( .not. present( xyz_UNS    ) ) .and. &
        &  ( .not. present( xyz_UNN    ) ) .and. &
        &  ( .not. present( xyz_VSN    ) ) .and. &
        &  ( .not. present( xyz_VSS    ) ) .and. &
        &  ( .not. present( xyz_VNN    ) ) .and. &
        &  ( .not. present( xyz_VNS    ) ) ) then
        FlagIncludeUV = .false.
      else
        call MessageNotify( 'E', module_name, 'Argument is inappropriate in SLTTExtArrExtMPICore.' )
      end if
    end if


    !====================================================================

    if ( FlagIncludeUV ) then
      narrsize = 2 + ncmax
    else
      narrsize = ncmax
    end if

    allocate( xyz_SendBufSS(0:imax-1, 1:dtjw, 1:kmax, narrsize) )
    allocate( xyz_SendBufSN(0:imax-1, 1:dtjw, 1:kmax, narrsize) )
    allocate( xyz_SendBufNS(0:imax-1, 1:dtjw, 1:kmax, narrsize) )
    allocate( xyz_SendBufNN(0:imax-1, 1:dtjw, 1:kmax, narrsize) )
    allocate( xyz_RecvBufSS(0:imax-1, 1:dtjw, 1:kmax, narrsize) )
    allocate( xyz_RecvBufSN(0:imax-1, 1:dtjw, 1:kmax, narrsize) )
    allocate( xyz_RecvBufNS(0:imax-1, 1:dtjw, 1:kmax, narrsize) )
    allocate( xyz_RecvBufNN(0:imax-1, 1:dtjw, 1:kmax, narrsize) )

    ! y_Array(1:dtjw) values are transfered (y_SN   ).
    !
    if ( myrank < (nprocs-1) ) then
      xyz_SendBufSN(:,:,:,1:ncmax) = xyzf_QMix(:,1:dtjw,:,:)
      if ( FlagIncludeUV ) then
        xyz_SendBufSN(:,:,:,ncmax+1) = xyz_U(:,1:dtjw,:)
        xyz_SendBufSN(:,:,:,ncmax+2) = xyz_V(:,1:dtjw,:)
      end if
      idest = myrank + 1
      call MPIWrapperISend(                  &
        & idest, imax, dtjw, kmax, narrsize, & ! (in)
        & xyz_SendBufSN(:,:,:,:),           & ! (in)
        & a_ireq_send_sn                     & ! (out)
        & )
    end if
    if ( myrank > 0 ) then
      idep = myrank - 1
      call MPIWrapperIRecv(                 &
        & idep, imax, dtjw, kmax, narrsize, & ! (in)
        & xyz_RecvBufSN(:,:,:,:),          & ! (out)
        & a_ireq_recv_sn                    & ! (out)
        & )
    end if
    !
    ! y_Array(jmax/2+1-dtjw:jmax/2) values are transfered (y_SS   ).
    !
    if ( myrank > 0 ) then
      xyz_SendBufSS(:,:,:,1:ncmax) = xyzf_QMix(:,jmax/2+1-dtjw:jmax/2,:,:)
      if ( FlagIncludeUV ) then
        xyz_SendBufSS(:,:,:,ncmax+1) = xyz_U(:,jmax/2+1-dtjw:jmax/2,:)
        xyz_SendBufSS(:,:,:,ncmax+2) = xyz_V(:,jmax/2+1-dtjw:jmax/2,:)
      end if
      idest = myrank - 1
      call MPIWrapperISend(                  &
        & idest, imax, dtjw, kmax, narrsize, & ! (in)
        & xyz_SendBufSS(:,:,:,:),           & ! (in)
        & a_ireq_send_ss                     & ! (out)
        & )
    end if
    if ( myrank < (nprocs-1) ) then
      idep = myrank + 1
      call MPIWrapperIRecv(                 &
        & idep, imax, dtjw, kmax, narrsize, & ! (in)
        & xyz_RecvBufSS(:,:,:,:),          & ! (out)
        & a_ireq_recv_ss                    & ! (out)
        & )
    end if
    !
    ! y_Array(jmax/2+1:jmax/2+dtjw) values are transfered (y_NN   ).
    !
    if ( myrank > 0 ) then
      xyz_SendBufNN(:,:,:,1:ncmax) = xyzf_QMix(:,jmax/2+1:jmax/2+dtjw,:,:)
      if ( FlagIncludeUV ) then
        xyz_SendBufNN(:,:,:,ncmax+1) = xyz_U(:,jmax/2+1:jmax/2+dtjw,:)
        xyz_SendBufNN(:,:,:,ncmax+2) = xyz_V(:,jmax/2+1:jmax/2+dtjw,:)
      end if
      idest = myrank - 1
      call MPIWrapperISend(                  &
        & idest, imax, dtjw, kmax, narrsize, & ! (in)
        & xyz_SendBufNN(:,:,:,:),           & ! (in)
        & a_ireq_send_nn                     & ! (out)
        & )
    end if
    if ( myrank < (nprocs-1) ) then
      idep = myrank + 1
      call MPIWrapperIRecv(                 &
        & idep, imax, dtjw, kmax, narrsize, & ! (in)
        & xyz_RecvBufNN(:,:,:,:),          & ! (out)
        & a_ireq_recv_nn                    & ! (out)
        & )
    end if
    !
    ! y_Array(jmax+1-dtjw:jmax) values are transfered (y_NS   ).
    !
    if ( myrank < (nprocs-1) ) then
      xyz_SendBufNS(:,:,:,1:ncmax) = xyzf_QMix(:,jmax+1-dtjw:jmax,:,:)
      if ( FlagIncludeUV ) then
        xyz_SendBufNS(:,:,:,ncmax+1) = xyz_U(:,jmax+1-dtjw:jmax,:)
        xyz_SendBufNS(:,:,:,ncmax+2) = xyz_V(:,jmax+1-dtjw:jmax,:)
      end if
      idest = myrank + 1
      call MPIWrapperISend(                  &
        & idest, imax, dtjw, kmax, narrsize, & ! (in)
        & xyz_SendBufNS(:,:,:,:),           & ! (in)
        & a_ireq_send_ns                     & ! (out)
        & )
    end if
    if ( myrank > 0 ) then
      idep = myrank - 1
      call MPIWrapperIRecv(                 &
        & idep, imax, dtjw, kmax, narrsize, & ! (in)
        & xyz_RecvBufNS(:,:,:,:),          & ! (out)
        & a_ireq_recv_ns                    & ! (out)
        & )
    end if

    if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_send_sn )
    if ( myrank > 0          ) call MPIWrapperWait( a_ireq_recv_sn )
    if ( myrank > 0          ) call MPIWrapperWait( a_ireq_send_ss )
    if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_recv_ss )
    if ( myrank > 0          ) call MPIWrapperWait( a_ireq_send_nn )
    if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_recv_nn )
    if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_send_ns )
    if ( myrank > 0          ) call MPIWrapperWait( a_ireq_recv_ns )

    ! y_Array(1:dtjw) values are transfered (y_SN   ).
    if ( myrank > 0 ) then
      xyzf_QMixSN    = xyz_RecvBufSN(:,:,:,1:ncmax)
      if ( FlagIncludeUV ) then
        xyz_USN    = xyz_RecvBufSN(:,:,:,ncmax+1)
        xyz_VSN    = xyz_RecvBufSN(:,:,:,ncmax+2)
      end if
    end if
    ! y_Array(jmax/2+1-dtjw:jmax/2) values are transfered (y_SS   ).
    if ( myrank < (nprocs-1) ) then
      xyzf_QMixSS    = xyz_RecvBufSS(:,:,:,1:ncmax)
      if ( FlagIncludeUV ) then
        xyz_USS    = xyz_RecvBufSS(:,:,:,ncmax+1)
        xyz_VSS    = xyz_RecvBufSS(:,:,:,ncmax+2)
      end if
    end if
    ! y_Array(jmax/2+1:jmax/2+dtjw) values are transfered (y_NN   ).
    if ( myrank < (nprocs-1) ) then
      xyzf_QMixNN    = xyz_RecvBufNN(:,:,:,1:ncmax)
      if ( FlagIncludeUV ) then
        xyz_UNN    = xyz_RecvBufNN(:,:,:,ncmax+1)
        xyz_VNN    = xyz_RecvBufNN(:,:,:,ncmax+2)
      end if
    end if
    ! y_Array(jmax+1-dtjw:jmax) values are transfered (y_NS   ).
    if ( myrank > 0 ) then
      xyzf_QMixNS    = xyz_RecvBufNS(:,:,:,1:ncmax)
      if ( FlagIncludeUV ) then
        xyz_UNS    = xyz_RecvBufNS(:,:,:,ncmax+1)
        xyz_VNS    = xyz_RecvBufNS(:,:,:,ncmax+2)
      end if
    end if

    deallocate( xyz_SendBufSS )
    deallocate( xyz_SendBufSN )
    deallocate( xyz_SendBufNS )
    deallocate( xyz_SendBufNN )
    deallocate( xyz_RecvBufSS )
    deallocate( xyz_RecvBufSN )
    deallocate( xyz_RecvBufNS )
    deallocate( xyz_RecvBufNN )


  end subroutine SLTTExtArrExtMPICore

end module sltt_extarr
