| Class | sltt_extarr |
| In: |
sltt/sltt_extarr.f90
|
| Subroutine : | |
| 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(in ) |
| 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) : | real(DP), intent(out) |
subroutine SLTTExtArrExt( y_ExtLatS, y_ExtLatN, x_SinLonS, x_CosLonS, x_SinLonN, x_CosLonN, xyzf_QMix, xyz_U, xyz_V, xyzf_ExtDQMixDLatS, xyzf_ExtDQMixDLatN, xyzf_ExtQMixS, xyzf_ExtQMixN, xyz_ExtUS, xyz_ExtUN, xyz_ExtVS, xyz_ExtVN )
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)
real(DP), allocatable :: xyza_SendBuf(:,:,:,:,:)
real(DP), allocatable :: xyza_RecvBuf(:,:,:,:,:)
real(DP), allocatable :: xyz_USPrev (:,:,:)
real(DP), allocatable :: xyz_UNPrev (:,:,:)
real(DP), allocatable :: xyz_USNext (:,:,:)
real(DP), allocatable :: xyz_UNNext (:,:,:)
real(DP), allocatable :: xyz_VSPrev (:,:,:)
real(DP), allocatable :: xyz_VNPrev (:,:,:)
real(DP), allocatable :: xyz_VSNext (:,:,:)
real(DP), allocatable :: xyz_VNNext (:,:,:)
real(DP), allocatable :: xyzf_QMixSPrev(:,:,:,:)
real(DP), allocatable :: xyzf_QMixNPrev(:,:,:,:)
real(DP), allocatable :: xyzf_QMixSNext(:,:,:,:)
real(DP), allocatable :: xyzf_QMixNNext(:,:,:,:)
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, 2+ncmax, 4) )
allocate( xyza_RecvBuf(0:imax-1, 1:dtjw, 1:kmax, 2+ncmax, 4) )
if ( myrank > 0 ) then
allocate( xyz_USPrev (0:imax-1, 1:dtjw, 1:kmax) )
allocate( xyz_UNPrev (0:imax-1, 1:dtjw, 1:kmax) )
allocate( xyz_VSPrev (0:imax-1, 1:dtjw, 1:kmax) )
allocate( xyz_VNPrev (0:imax-1, 1:dtjw, 1:kmax) )
allocate( xyzf_QMixSPrev(0:imax-1, 1:dtjw, 1:kmax, ncmax) )
allocate( xyzf_QMixNPrev(0:imax-1, 1:dtjw, 1:kmax, ncmax) )
end if
if ( myrank < (nprocs-1) ) then
allocate( xyz_USNext (0:imax-1, 1:dtjw, 1:kmax) )
allocate( xyz_UNNext (0:imax-1, 1:dtjw, 1:kmax) )
allocate( xyz_VSNext (0:imax-1, 1:dtjw, 1:kmax) )
allocate( xyz_VNNext (0:imax-1, 1:dtjw, 1:kmax) )
allocate( xyzf_QMixSNext(0:imax-1, 1:dtjw, 1:kmax, ncmax) )
allocate( xyzf_QMixNNext(0:imax-1, 1:dtjw, 1:kmax, ncmax) )
end if
! y_Array(1:dtjw) values are transfered (y_SPrev).
!
if ( myrank < (nprocs-1) ) then
xyza_SendBuf(:,:,:,1 ,1) = xyz_U (:,1:dtjw,:)
xyza_SendBuf(:,:,:,2 ,1) = xyz_V (:,1:dtjw,:)
xyza_SendBuf(:,:,:,3:2+ncmax,1) = xyzf_QMix(:,1:dtjw,:,:)
idest = myrank + 1
call MPIWrapperISend( idest, imax, dtjw, kmax, 2+ncmax, xyza_SendBuf(:,:,:,:,1), a_ireq_send(1) )
end if
if ( myrank > 0 ) then
idep = myrank - 1
call MPIWrapperIRecv( idep, imax, dtjw, kmax, 2+ncmax, xyza_RecvBuf(:,:,:,:,1), a_ireq_recv(1) )
end if
!
! y_Array(jmax/2+1-dtjw:jmax/2) values are transfered (y_SNext).
!
if ( myrank > 0 ) then
xyza_SendBuf(:,:,:,1 ,2) = xyz_U (:,jmax/2+1-dtjw:jmax/2,:)
xyza_SendBuf(:,:,:,2 ,2) = xyz_V (:,jmax/2+1-dtjw:jmax/2,:)
xyza_SendBuf(:,:,:,3:2+ncmax,2) = xyzf_QMix(:,jmax/2+1-dtjw:jmax/2,:,:)
idest = myrank - 1
call MPIWrapperISend( idest, imax, dtjw, kmax, 2+ncmax, xyza_SendBuf(:,:,:,:,2), a_ireq_send(2) )
end if
if ( myrank < (nprocs-1) ) then
idep = myrank + 1
call MPIWrapperIRecv( idep, imax, dtjw, kmax, 2+ncmax, xyza_RecvBuf(:,:,:,:,2), a_ireq_recv(2) )
end if
!
! y_Array(jmax/2+1:jmax/2+dtjw) values are transfered (y_NNext).
!
if ( myrank > 0 ) then
xyza_SendBuf(:,:,:,1 ,3) = xyz_U (:,jmax/2+1:jmax/2+dtjw,:)
xyza_SendBuf(:,:,:,2 ,3) = xyz_V (:,jmax/2+1:jmax/2+dtjw,:)
xyza_SendBuf(:,:,:,3:2+ncmax,3) = xyzf_QMix(:,jmax/2+1:jmax/2+dtjw,:,:)
idest = myrank - 1
call MPIWrapperISend( idest, imax, dtjw, kmax, 2+ncmax, xyza_SendBuf(:,:,:,:,3), a_ireq_send(3) )
end if
if ( myrank < (nprocs-1) ) then
idep = myrank + 1
call MPIWrapperIRecv( idep, imax, dtjw, kmax, 2+ncmax, xyza_RecvBuf(:,:,:,:,3), a_ireq_recv(3) )
end if
!
! y_Array(jmax+1-dtjw:jmax) values are transfered (y_NPrev).
!
if ( myrank < (nprocs-1) ) then
xyza_SendBuf(:,:,:,1 ,4) = xyz_U (:,jmax+1-dtjw:jmax,:)
xyza_SendBuf(:,:,:,2 ,4) = xyz_V (:,jmax+1-dtjw:jmax,:)
xyza_SendBuf(:,:,:,3:2+ncmax,4) = xyzf_QMix(:,jmax+1-dtjw:jmax,:,:)
idest = myrank + 1
call MPIWrapperISend( idest, imax, dtjw, kmax, 2+ncmax, xyza_SendBuf(:,:,:,:,4), a_ireq_send(4) )
end if
if ( myrank > 0 ) then
idep = myrank - 1
call MPIWrapperIRecv( idep, imax, dtjw, kmax, 2+ncmax, xyza_RecvBuf(:,:,:,:,4), a_ireq_recv(4) )
end if
if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_send(1) )
if ( myrank > 0 ) call MPIWrapperWait( a_ireq_recv(1) )
if ( myrank > 0 ) call MPIWrapperWait( a_ireq_send(2) )
if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_recv(2) )
if ( myrank > 0 ) call MPIWrapperWait( a_ireq_send(3) )
if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_recv(3) )
if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_send(4) )
if ( myrank > 0 ) call MPIWrapperWait( a_ireq_recv(4) )
! y_Array(1:dtjw) values are transfered (y_SPrev).
if ( myrank > 0 ) then
xyz_USPrev = xyza_RecvBuf(:,:,:,1 ,1)
xyz_VSPrev = xyza_RecvBuf(:,:,:,2 ,1)
xyzf_QMixSPrev = xyza_RecvBuf(:,:,:,3:2+ncmax,1)
end if
! y_Array(jmax/2+1-dtjw:jmax/2) values are transfered (y_SNext).
if ( myrank < (nprocs-1) ) then
xyz_USNext = xyza_RecvBuf(:,:,:,1 ,2)
xyz_VSNext = xyza_RecvBuf(:,:,:,2 ,2)
xyzf_QMixSNext = xyza_RecvBuf(:,:,:,3:2+ncmax,2)
end if
! y_Array(jmax/2+1:jmax/2+dtjw) values are transfered (y_NNext).
if ( myrank < (nprocs-1) ) then
xyz_UNNext = xyza_RecvBuf(:,:,:,1 ,3)
xyz_VNNext = xyza_RecvBuf(:,:,:,2 ,3)
xyzf_QMixNNext = xyza_RecvBuf(:,:,:,3:2+ncmax,3)
end if
! y_Array(jmax+1-dtjw:jmax) values are transfered (y_NPrev).
if ( myrank > 0 ) then
xyz_UNPrev = xyza_RecvBuf(:,:,:,1 ,4)
xyz_VNPrev = xyza_RecvBuf(:,:,:,2 ,4)
xyzf_QMixNPrev = xyza_RecvBuf(:,:,:,3:2+ncmax,4)
end if
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_USNext(i,dtjw-(j-1),k)
xyz_ExtVS(i,1-j,k) = xyz_VSNext(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_QMixSNext(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_USPrev(i,j,k)
xyz_ExtVS(i,jmax/2+j,k) = xyz_VSPrev(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_QMixSPrev(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_UNPrev(i,dtjw-(j-1),k)
xyz_ExtVN(i,1-j,k) = xyz_VNPrev(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_QMixNPrev(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_UNNext(i,j,k)
xyz_ExtVN(i,jmax/2+j,k) = xyz_VNNext(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_QMixNNext(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( xyz_USPrev )
deallocate( xyz_UNPrev )
deallocate( xyz_VSPrev )
deallocate( xyz_VNPrev )
deallocate( xyzf_QMixSPrev )
deallocate( xyzf_QMixNPrev )
end if
if ( myrank < (nprocs-1) ) then
deallocate( xyz_USNext )
deallocate( xyz_UNNext )
deallocate( xyz_VSNext )
deallocate( xyz_VNNext )
deallocate( xyzf_QMixSNext )
deallocate( xyzf_QMixNNext )
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 : | |||
| 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 : | real(DP), intent(in )
| ||
| 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) | ||
| PoleMethod : | character(*), intent(in )
|
subroutine SLTTExtArrExt2( x_SinLonS, x_CosLonS, x_SinLonN, x_CosLonN, xyzf_QMix, pm, xyzf_ExtQMixS, xyzf_ExtQMixN, PoleMethod )
! メッセージ出力
! 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_QMixSPrev(:,:,:,:)
real(DP), allocatable :: xyzf_QMixNPrev(:,:,:,:)
real(DP), allocatable :: xyzf_QMixSNext(:,:,:,:)
real(DP), allocatable :: xyzf_QMixNNext(:,:,:,:)
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_QMixSPrev(0:imax-1, 1:dtjw, 1:kmax, ncmax) )
allocate( xyzf_QMixNPrev(0:imax-1, 1:dtjw, 1:kmax, ncmax) )
end if
if ( myrank < (nprocs-1) ) then
allocate( xyzf_QMixSNext(0:imax-1, 1:dtjw, 1:kmax, ncmax) )
allocate( xyzf_QMixNNext(0:imax-1, 1:dtjw, 1:kmax, ncmax) )
end if
! y_Array(1:dtjw) values are transfered (y_SPrev).
!
if ( myrank < (nprocs-1) ) then
xyza_SendBuf(:,:,:,:,1) = xyzf_QMix(:,1:dtjw,:,:)
idest = myrank + 1
call MPIWrapperISend( idest, imax, dtjw, kmax, ncmax, xyza_SendBuf(:,:,:,:,1), a_ireq_send(1) )
end if
if ( myrank > 0 ) then
idep = myrank - 1
call MPIWrapperIRecv( idep, imax, dtjw, kmax, ncmax, xyza_RecvBuf(:,:,:,:,1), a_ireq_recv(1) )
end if
!
! y_Array(jmax/2+1-dtjw:jmax/2) values are transfered (y_SNext).
!
if ( myrank > 0 ) then
xyza_SendBuf(:,:,:,:,2) = xyzf_QMix(:,jmax/2+1-dtjw:jmax/2,:,:)
idest = myrank - 1
call MPIWrapperISend( idest, imax, dtjw, kmax, ncmax, xyza_SendBuf(:,:,:,:,2), a_ireq_send(2) )
end if
if ( myrank < (nprocs-1) ) then
idep = myrank + 1
call MPIWrapperIRecv( idep, imax, dtjw, kmax, ncmax, xyza_RecvBuf(:,:,:,:,2), a_ireq_recv(2) )
end if
!
! y_Array(jmax/2+1:jmax/2+dtjw) values are transfered (y_NNext).
!
if ( myrank > 0 ) then
xyza_SendBuf(:,:,:,:,3) = xyzf_QMix(:,jmax/2+1:jmax/2+dtjw,:,:)
idest = myrank - 1
call MPIWrapperISend( idest, imax, dtjw, kmax, ncmax, xyza_SendBuf(:,:,:,:,3), a_ireq_send(3) )
end if
if ( myrank < (nprocs-1) ) then
idep = myrank + 1
call MPIWrapperIRecv( idep, imax, dtjw, kmax, ncmax, xyza_RecvBuf(:,:,:,:,3), a_ireq_recv(3) )
end if
!
! y_Array(jmax+1-dtjw:jmax) values are transfered (y_NPrev).
!
if ( myrank < (nprocs-1) ) then
xyza_SendBuf(:,:,:,:,4) = xyzf_QMix(:,jmax+1-dtjw:jmax,:,:)
idest = myrank + 1
call MPIWrapperISend( idest, imax, dtjw, kmax, ncmax, xyza_SendBuf(:,:,:,:,4), a_ireq_send(4) )
end if
if ( myrank > 0 ) then
idep = myrank - 1
call MPIWrapperIRecv( idep, imax, dtjw, kmax, ncmax, xyza_RecvBuf(:,:,:,:,4), a_ireq_recv(4) )
end if
if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_send(1) )
if ( myrank > 0 ) call MPIWrapperWait( a_ireq_recv(1) )
if ( myrank > 0 ) call MPIWrapperWait( a_ireq_send(2) )
if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_recv(2) )
if ( myrank > 0 ) call MPIWrapperWait( a_ireq_send(3) )
if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_recv(3) )
if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_send(4) )
if ( myrank > 0 ) call MPIWrapperWait( a_ireq_recv(4) )
! y_Array(1:dtjw) values are transfered (y_SPrev).
if ( myrank > 0 ) then
xyzf_QMixSPrev = xyza_RecvBuf(:,:,:,:,1)
end if
! y_Array(jmax/2+1-dtjw:jmax/2) values are transfered (y_SNext).
if ( myrank < (nprocs-1) ) then
xyzf_QMixSNext = xyza_RecvBuf(:,:,:,:,2)
end if
! y_Array(jmax/2+1:jmax/2+dtjw) values are transfered (y_NNext).
if ( myrank < (nprocs-1) ) then
xyzf_QMixNNext = xyza_RecvBuf(:,:,:,:,3)
end if
! y_Array(jmax+1-dtjw:jmax) values are transfered (y_NPrev).
if ( myrank > 0 ) then
xyzf_QMixNPrev = xyza_RecvBuf(:,:,:,:,4)
end if
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_LatSNext(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_QMixSNext(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_QMixSPrev(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_QMixNPrev(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_QMixNNext(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_QMixSPrev )
deallocate( xyzf_QMixNPrev )
end if
if ( myrank < (nprocs-1) ) then
deallocate( xyzf_QMixSNext )
deallocate( xyzf_QMixNNext )
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 : | |
| 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(in ) |
| 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) : | real(DP), intent(out) |
subroutine SLTTExtArrInit( x_LonS, y_LatS, x_LonN, y_LatN, x_ExtLonS, y_ExtLatS, x_ExtLonN, y_ExtLatN )
!
! 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(4)
integer :: a_ireq_recv(4)
real(DP), allocatable :: y_LatSPrev(:)
real(DP), allocatable :: y_LatNPrev(:)
real(DP), allocatable :: y_LatSNext(:)
real(DP), allocatable :: y_LatNNext(:)
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_LatSPrev(1:dtjw) )
allocate( y_LatNPrev(1:dtjw) )
end if
if ( myrank < (nprocs-1) ) then
allocate( y_LatSNext(1:dtjw) )
allocate( y_LatNNext(1:dtjw) )
end if
! y_Lat(1:dtjw) values are transfered (y_LatSPrev).
!
if ( myrank < (nprocs-1) ) then
idest = myrank + 1
call MPIWrapperISend( idest, dtjw, y_Lat(1:dtjw), a_ireq_send(1) )
end if
if ( myrank > 0 ) then
idep = myrank - 1
call MPIWrapperIRecv( idep, dtjw, y_LatSPrev, a_ireq_recv(1) )
end if
!
! y_Lat(jmax/2+1-dtjw:jmax/2) values are transfered (y_LatSNext).
!
if ( myrank > 0 ) then
idest = myrank - 1
call MPIWrapperISend( idest, dtjw, y_Lat(jmax/2+1-dtjw:jmax/2), a_ireq_send(2) )
end if
if ( myrank < (nprocs-1) ) then
idep = myrank + 1
call MPIWrapperIRecv( idep, dtjw, y_LatSNext, a_ireq_recv(2) )
end if
!
! y_Lat(jmax/2+1:jmax/2+dtjw) values are transfered (y_LatNNext).
!
if ( myrank > 0 ) then
idest = myrank - 1
call MPIWrapperISend( idest, dtjw, y_Lat(jmax/2+1:jmax/2+dtjw), a_ireq_send(3) )
end if
if ( myrank < (nprocs-1) ) then
idep = myrank + 1
call MPIWrapperIRecv( idep, dtjw, y_LatNNext, a_ireq_recv(3) )
end if
!
! y_Lat(jmax+1-dtjw:jmax) values are transfered (y_LatNPrev).
!
if ( myrank < (nprocs-1) ) then
idest = myrank + 1
call MPIWrapperISend( idest, dtjw, y_Lat(jmax+1-dtjw:jmax), a_ireq_send(4) )
end if
if ( myrank > 0 ) then
idep = myrank - 1
call MPIWrapperIRecv( idep, dtjw, y_LatNPrev, a_ireq_recv(4) )
end if
if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_send(1) )
if ( myrank > 0 ) call MPIWrapperWait( a_ireq_recv(1) )
if ( myrank > 0 ) call MPIWrapperWait( a_ireq_send(2) )
if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_recv(2) )
if ( myrank > 0 ) call MPIWrapperWait( a_ireq_send(3) )
if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_recv(3) )
if ( myrank < (nprocs-1) ) call MPIWrapperWait( a_ireq_send(4) )
if ( myrank > 0 ) call MPIWrapperWait( a_ireq_recv(4) )
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_LatSNext(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_LatSPrev(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_LatNPrev(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_LatNNext(j)
end do
end if
if ( myrank > 0 ) then
deallocate( y_LatSPrev )
deallocate( y_LatNPrev )
end if
if ( myrank < (nprocs-1) ) then
deallocate( y_LatSNext )
deallocate( y_LatNNext )
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
| Subroutine : |
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
| Constant : | |||
| version = ’$Name: dcpam5-20140204-4 $’ // ’$Id: sltt_extarr.f90,v 1.2 2013-01-27 11:26:14 yot Exp $’ : | character(*), parameter
|