! $Id$ ! ! This module provides routines for obtaining slice data. ! module Slices_methods ! use Cdata ! implicit none ! interface assign_slices_scal module procedure assign_slices_f_scal module procedure assign_slices_sep_scal endinterface interface assign_slices_vec module procedure assign_slices_f_vec module procedure assign_slices_sep_vec endinterface interface store_slices module procedure store_slices_scal module procedure store_slices_vec endinterface interface process_slices module procedure process_slices_func module procedure process_slices_fac endinterface interface alloc_rslice module procedure alloc_rslice_scal module procedure alloc_rslice_scal_2D module procedure alloc_rslice_vec endinterface interface alloc_slice_buffers module procedure alloc_slice_buffers_scal module procedure alloc_slice_buffers_vec endinterface public :: store_slices public :: process_slices, addto_slices public :: assign_slices_vec public :: assign_slices_scal, assign_slices_f_scal public :: alloc_slice_buffers public :: nullify_slice_pointers public :: write_rslice_position public :: alloc_rslice, prep_rslice, store_dims_rslice public :: exp2d, log2d, abs2d ! type slices_scal ! real, dimension (:,:), allocatable :: xy ! real, dimension (:,:), allocatable :: xz ! real, dimension (:,:), allocatable :: xz2 ! real, dimension (:,:), allocatable :: yz ! real, dimension (:,:), allocatable :: xy2 ! real, dimension (:,:), allocatable :: xy3 ! real, dimension (:,:), allocatable :: xy4 ! real, dimension (:,:,:,:,:), allocatable :: r ! endtype slices_scal ! type slices_vec ! real, dimension (:,:,:), allocatable :: xy ! real, dimension (:,:,:), allocatable :: xz ! real, dimension (:,:,:), allocatable :: xz2 ! real, dimension (:,:,:), allocatable :: yz ! real, dimension (:,:,:), allocatable :: xy2 ! real, dimension (:,:,:), allocatable :: xy3 ! real, dimension (:,:,:), allocatable :: xy4 ! real, dimension (:,:,:,:,:,:), allocatable :: r ! endtype slices_vec private ! ! For spherical slices in Cartesian geometry: ! real, dimension(:), allocatable :: cph_slice,sph_slice,cth_slice,sth_slice integer :: ith_min,ith_max,iph_min,iph_max real, dimension(:,:,:,:,:), allocatable :: rslice_interp_weights integer, dimension(:,:,:), allocatable :: rslice_adjec_corn_inds contains !*********************************************************************** subroutine alloc_slice_buffers_scal(xy,xz,yz,xy2,xy3,xy4,xz2,r) !(slices) !type(slices_scal) :: slices real, dimension(:,:), allocatable :: xy,xz,yz,xy2,xy3,xy4,xz2 real, dimension(:,:,:,:,:), allocatable, optional :: r if (lwrite_slice_xy .and..not.allocated(xy) ) allocate(xy (nx,ny)) if (lwrite_slice_xz .and..not.allocated(xz) ) allocate(xz (nx,nz)) if (lwrite_slice_yz .and..not.allocated(yz) ) allocate(yz (ny,nz)) if (lwrite_slice_xy2.and..not.allocated(xy2)) allocate(xy2(nx,ny)) if (lwrite_slice_xy3.and..not.allocated(xy3)) allocate(xy3(nx,ny)) if (lwrite_slice_xy4.and..not.allocated(xy4)) allocate(xy4(nx,ny)) if (lwrite_slice_xz2.and..not.allocated(xz2)) allocate(xz2(nx,nz)) if (lwrite_slice_r.and.present(r)) then if (.not.allocated(r)) call alloc_rslice(r) endif endsubroutine alloc_slice_buffers_scal !*********************************************************************** subroutine alloc_slice_buffers_vec(xy,xz,yz,xy2,xy3,xy4,xz2,r,ncomp) use General, only: ioptest real, dimension(:,:,:), allocatable :: xy,xz,yz,xy2,xy3,xy4,xz2 real, dimension(:,:,:,:,:,:), allocatable, optional :: r integer, optional :: ncomp integer :: nc nc=ioptest(ncomp,3) if (lwrite_slice_xy .and..not.allocated(xy) ) allocate(xy (nx,ny,nc)) if (lwrite_slice_xz .and..not.allocated(xz) ) allocate(xz (nx,nz,nc)) if (lwrite_slice_yz .and..not.allocated(yz) ) allocate(yz (ny,nz,nc)) if (lwrite_slice_xy2.and..not.allocated(xy2)) allocate(xy2(nx,ny,nc)) if (lwrite_slice_xy3.and..not.allocated(xy3)) allocate(xy3(nx,ny,nc)) if (lwrite_slice_xy4.and..not.allocated(xy4)) allocate(xy4(nx,ny,nc)) if (lwrite_slice_xz2.and..not.allocated(xz2)) allocate(xz2(nx,nz,nc)) if (lwrite_slice_r.and.present(r)) then if (.not.allocated(r)) call alloc_rslice(r,nc) endif endsubroutine alloc_slice_buffers_vec !*********************************************************************** subroutine assign_slices_sep_scal(slices,xy,xz,yz,xy2,xy3,xy4,xz2,r) ! ! Assignment of slice pointers according to slice selection switches. ! ! 12-apr-17/MR: coded ! 8-oct-18/MR: made some pars optional ! type (slice_data), intent(OUT):: slices real, dimension(:,:), target, intent(IN) :: xy real, dimension(:,:), target, optional, intent(IN) :: xz,yz,xy2,xy3,xy4,xz2 real, dimension(:,:,:,:,:), target, optional, intent(IN) :: r if (lwrite_slice_yz ) slices%yz =>yz if (lwrite_slice_xz ) slices%xz =>xz if (lwrite_slice_xy ) slices%xy =>xy if (lwrite_slice_xy2) slices%xy2=>xy2 if (lwrite_slice_xy3) slices%xy3=>xy3 if (lwrite_slice_xy4) slices%xy4=>xy4 if (lwrite_slice_xz2) slices%xz2=>xz2 if (lwrite_slice_r) then if (present(r)) call interp_rslice(slices%r,rslice=r) endif slices%ready=.true. endsubroutine assign_slices_sep_scal !*********************************************************************** subroutine assign_slices_sep_vec(slices,xy,xz,yz,xy2,xy3,xy4,xz2,r,ncomp) ! ! Assignment of slice pointers according to slice selection switches. ! ! 12-apr-17/MR: coded ! use General, only: ioptest ! type (slice_data) , intent(INOUT):: slices real, dimension(:,:,:), target, intent(IN) :: xy,xz,yz,xy2,xy3,xy4,xz2 real, dimension(:,:,:,:,:,:), target, optional, intent(IN) :: r integer, optional , intent(IN) :: ncomp integer :: nc nc=ioptest(ncomp,3) if (slices%index>=nc) then slices%ready=.false. else slices%index=slices%index+1 if (lwrite_slice_yz ) slices%yz => yz(:,:,slices%index) if (lwrite_slice_xz ) slices%xz => xz(:,:,slices%index) if (lwrite_slice_xy ) slices%xy => xy(:,:,slices%index) if (lwrite_slice_xy2) slices%xy2=>xy2(:,:,slices%index) if (lwrite_slice_xy3) slices%xy3=>xy3(:,:,slices%index) if (lwrite_slice_xy4) slices%xy4=>xy4(:,:,slices%index) if (lwrite_slice_xz2) slices%xz2=>xz2(:,:,slices%index) if (lwrite_slice_r) then if (present(r)) & call interp_rslice(slices%r,rslice=r(:,:,:,:,:,slices%index)) endif if (slices%index<=nc) slices%ready=.true. endif endsubroutine assign_slices_sep_vec !*********************************************************************** subroutine assign_slices_f_scal(slices,f,ind1,ind2) ! ! Copying of scalar data from f-array to arrays assigned to the slice pointers ! according to slice selection switches. ! ! 12-apr-17/MR: coded ! 28-may-19/MR: added optional argument ind2 ! use General, only: ioptest type (slice_data) , intent(OUT):: slices real, dimension(mx,my,mz,*), intent(IN) :: f integer , intent(IN) :: ind1 integer , optional , intent(IN) :: ind2 integer :: ind2_ ind2_=ioptest(ind2,ind1) if (lwrite_slice_yz) slices%yz =f(ix_loc,m1:m2 ,n1:n2 ,ind1) !:ind2_) if (lwrite_slice_xz) slices%xz =f(l1:l2 ,iy_loc ,n1:n2 ,ind1) !:ind2_) if (lwrite_slice_xy) slices%xy =f(l1:l2 ,m1:m2 ,iz_loc ,ind1) !:ind2_) if (lwrite_slice_xy2) slices%xy2=f(l1:l2 ,m1:m2 ,iz2_loc,ind1) !:ind2_) if (lwrite_slice_xy3) slices%xy3=f(l1:l2 ,m1:m2 ,iz3_loc,ind1) !:ind2_) if (lwrite_slice_xy4) slices%xy4=f(l1:l2 ,m1:m2 ,iz4_loc,ind1) !:ind2_) if (lwrite_slice_xz2) slices%xz2=f(l1:l2 ,iy2_loc,n1:n2 ,ind1) !:ind2_) if (lwrite_slice_r ) call interp_rslice(slices%r,f,ind1) slices%ready=.true. endsubroutine assign_slices_f_scal !*********************************************************************** subroutine assign_slices_f_vec(slices,f,ind,ncomp) ! ! Copying of multi-component data from f-array to arrays assigned to the slice pointers ! according to slice selection switches. If ncomp is not present the quantity ! in f(:,:,:,ind:ind+2) is considered a proper vector (to be transformed in ! Yin-Yang context), otherwise a general multi-component quantity (not to be ! transformed). ! ! 12-apr-17/MR: coded ! use General, only: transform_thph_yy_other, ioptest type (slice_data) , intent(INOUT):: slices real, dimension(:,:,:,:), intent(IN) :: f integer , intent(IN) :: ind integer, optional , intent(IN) :: ncomp real, dimension(:,:,:,:), allocatable, save :: transformed integer :: nc nc=ioptest(ncomp,3) if (slices%index>=nc) then slices%ready=.false. else slices%index=slices%index+1 if (lwrite_slice_yz) then if (lyang.and..not.present(ncomp).and.slices%index>=2) then ! ! On Yang grid: transform theta and phi components of vector to Yin-grid basis. ! (phi component is saved in transformed for use in next call. ! if (slices%index==2) then if (.not.allocated(transformed)) allocate(transformed(1,ny,nz,3)) call transform_thph_yy_other(f(ix_loc:ix_loc,:,:,ind+1:ind+2),m1,m2,n1,n2,transformed) endif ! ! theta component is used immediately, phi component with next call. ! slices%yz=transformed(1,:,:,slices%index-1) else slices%yz=f(ix_loc,m1:m2,n1:n2,ind-1+slices%index) endif endif if (lwrite_slice_xz) & slices%xz =f(l1:l2,iy_loc ,n1:n2 ,ind-1+slices%index) if (lwrite_slice_xy) & slices%xy =f(l1:l2,m1:m2 ,iz_loc ,ind-1+slices%index) if (lwrite_slice_xy2) & slices%xy2=f(l1:l2,m1:m2 ,iz2_loc,ind-1+slices%index) if (lwrite_slice_xy3) & slices%xy3=f(l1:l2,m1:m2 ,iz3_loc,ind-1+slices%index) if (lwrite_slice_xy4) & slices%xy4=f(l1:l2,m1:m2 ,iz4_loc,ind-1+slices%index) if (lwrite_slice_xz2) & slices%xz2=f(l1:l2,iy2_loc,n1:n2 ,ind-1+slices%index) if (lwrite_slice_r ) call interp_rslice(slices%r,f,ind-1+slices%index) if (slices%index<=nc) slices%ready=.true. endif endsubroutine assign_slices_f_vec !*********************************************************************** subroutine store_slices_scal(pencil,xy,xz,yz,xy2,xy3,xy4,xz2,r) ! ! Stores scalar data from pencil in auxiliary arrays when needed. ! ! 12-apr-17/MR: coded ! 8-oct-18/MR: made some pars optional ! real, dimension(nx) , intent(IN) :: pencil real, dimension(:,:), intent(OUT) :: xy real, dimension(:,:), intent(OUT),optional :: xz,yz,xy2,xy3,xy4,xz2 real, dimension(:,:,:,:,:), intent(OUT), optional :: r if (n==iz_loc) xy(:,m-m1+1)=pencil if (present(xz) .and.m==iy_loc) xz(:,n-n1+1)=pencil if (present(yz) .and.lwrite_slice_yz) yz(m-m1+1,n-n1+1)=pencil(ix_loc-l1+1) if (present(xy2).and.n==iz2_loc) xy2(:,m-m1+1)=pencil if (present(xy3).and.n==iz3_loc) xy3(:,m-m1+1)=pencil if (present(xy4).and.n==iz4_loc) xy4(:,m-m1+1)=pencil if (present(xz2).and.m==iy2_loc) xz2(:,n-n1+1)=pencil if (lwrite_slice_r) then if (present(r)) call store_rslice_scal(pencil,r) endif endsubroutine store_slices_scal !*********************************************************************** subroutine store_slices_vec(pencil,xy,xz,yz,xy2,xy3,xy4,xz2,r,ncomp) ! ! Stores multi-component data from pencil in auxiliary arrays when needed. ! ! 12-apr-17/MR: coded ! use General, only: transform_thph_yy, ioptest real, dimension(:,:) , intent(IN) :: pencil real, dimension(:,:,:), intent(OUT):: xy,xz,yz,xy2,xy3,xy4,xz2 real, dimension(:,:,:,:,:,:), intent(OUT), optional :: r integer, optional , intent(IN) :: ncomp real, dimension(1,3) :: tmp !nc=ioptest(ncomp,3) ! ! Note: i[xyz][234]*_loc is the index with respect to array with ghost zones. ! if (lwrite_slice_yz) then if (lyang.and..not.present(ncomp).and.size(pencil,2)==3) then ! ! On Yang grid: transform theta and phi components of vector to Yin-grid basis. ! call transform_thph_yy(pencil(ix_loc-l1+1:ix_loc-l1+1,:),(/1,1,1/),tmp) yz(m-m1+1,n-n1+1,:)=tmp(1,:) else yz(m-m1+1,n-n1+1,:)=pencil(ix_loc-l1+1,:) endif endif ! if (m==iy_loc ) xz (:,n-n1+1,:)=pencil if (n==iz_loc ) xy (:,m-m1+1,:)=pencil if (n==iz2_loc) xy2(:,m-m1+1,:)=pencil if (n==iz3_loc) xy3(:,m-m1+1,:)=pencil if (n==iz4_loc) xy4(:,m-m1+1,:)=pencil if (m==iy2_loc) xz2(:,n-n1+1,:)=pencil if (lwrite_slice_r) then if (present(r)) call store_rslice_vec(pencil,r) endif endsubroutine store_slices_vec !*********************************************************************** function exp2d(arr) result(res) real, dimension(:,:) :: arr real, dimension(size(arr,1),size(arr,2)) :: res res=exp(arr) endfunction exp2d !*********************************************************************** function log2d(arr) result(res) real, dimension(:,:) :: arr real, dimension(size(arr,1),size(arr,2)) :: res res=alog(arr) endfunction log2d !*********************************************************************** function abs2d(arr) result(res) real, dimension(:,:) :: arr real, dimension(size(arr,1),size(arr,2)) :: res res=abs(arr) endfunction abs2d !*********************************************************************** subroutine process_slices_func(slices,func) type(slice_data), intent(INOUT):: slices interface function func(arr) result(res) real, dimension(:,:) :: arr real, dimension(size(arr,1),size(arr,2)) :: res end function func endinterface if (lwrite_slice_yz ) slices%yz =func(slices%yz) if (lwrite_slice_xz ) slices%xz =func(slices%xz) if (lwrite_slice_xy ) slices%xy =func(slices%xy) if (lwrite_slice_xy2) slices%xy2=func(slices%xy2) if (lwrite_slice_xy3) slices%xy3=func(slices%xy3) if (lwrite_slice_xy4) slices%xy4=func(slices%xy4) if (lwrite_slice_xz2) slices%xz2=func(slices%xz2) if (lwrite_slice_r ) & ! ! If rank has point, process. ! where(rslice_adjec_corn_inds(:,:,1)/=0) slices%r=func(slices%r) endsubroutine process_slices_func !*********************************************************************** subroutine process_slices_fac(slices,fac) type(slice_data), intent(INOUT):: slices real, intent(IN) :: fac if (lwrite_slice_yz ) slices%yz =fac*slices%yz if (lwrite_slice_xz ) slices%xz =fac*slices%xz if (lwrite_slice_xy ) slices%xy =fac*slices%xy if (lwrite_slice_xy2) slices%xy2=fac*slices%xy2 if (lwrite_slice_xy3) slices%xy3=fac*slices%xy3 if (lwrite_slice_xy4) slices%xy4=fac*slices%xy4 if (lwrite_slice_xz2) slices%xz2=fac*slices%xz2 if (lwrite_slice_r ) slices%r =fac*slices%r endsubroutine process_slices_fac !*********************************************************************** subroutine addto_slices(slices,pencil) ! ! Adds pencil data to scalar slices according to slice selection switches. ! ! 12-apr-17/MR: coded ! use Messages, only: not_implemented type (slice_data), intent(INOUT):: slices real, dimension(nx),intent(IN) :: pencil integer :: i if (lwrite_slice_yz ) slices%yz=slices%yz+pencil(ix_loc-l1+1) if (lwrite_slice_xz ) then do i=1,nz; slices%xz (:,i)=slices%xz (:,i)+pencil; enddo endif if (lwrite_slice_xy ) then do i=1,ny; slices%xy (:,i)=slices%xy (:,i)+pencil; enddo endif if (lwrite_slice_xy2) then do i=1,ny; slices%xy2(:,i)=slices%xy2(:,i)+pencil; enddo endif if (lwrite_slice_xy3) then do i=1,ny; slices%xy3(:,i)=slices%xy3(:,i)+pencil; enddo endif if (lwrite_slice_xy4) then do i=1,ny; slices%xy4(:,i)=slices%xy4(:,i)+pencil; enddo endif if (lwrite_slice_xz2) then do i=1,nz; slices%xz2(:,i)=slices%xz2(:,i)+pencil; enddo endif if (lwrite_slice_r) call not_implemented('addto_slices','for r-slices') endsubroutine addto_slices !*********************************************************************** subroutine nullify_slice_pointers(slices) ! ! Nullifies all pointers in slices struc. ! ! 12-apr-17/MR: coded ! type (slice_data), intent(OUT):: slices nullify(slices%yz) nullify(slices%xz) nullify(slices%xy) nullify(slices%xy2) nullify(slices%xy3) nullify(slices%xy4) nullify(slices%xz2) nullify(slices%r) endsubroutine nullify_slice_pointers !*********************************************************************** subroutine prep_rslice integer :: ith,iph,ll,mm,nn,iph_min_loc,iph_max_loc,iphm,iph_min_loc0 real :: xs,ys,zs,dth,dphi,dth0,dph0,dxs,dys,dzs,wsum,poms logical :: lfound,l0 if (allocated(cph_slice)) deallocate(cph_slice,sph_slice,cth_slice,sth_slice) allocate(cph_slice(nph_rslice),sph_slice(nph_rslice)) allocate(cth_slice(nth_rslice),sth_slice(nth_rslice)) dth=pi/nth_rslice; dphi=2*pi/nph_rslice; dth0=.5*dth; dph0=0. !pi/4. do iph=1,nph_rslice cph_slice(iph)=cos((iph-1)*dphi+dph0) sph_slice(iph)=sin((iph-1)*dphi+dph0) enddo ith_min=0; iph_min=0; l0=.true. do ith=1,nth_rslice ! loop over all parallels of the r-slice cth_slice(ith)=r_rslice*cos(dth0+(ith-1)*dth) sth_slice(ith)=r_rslice*sin(dth0+(ith-1)*dth) zs=cth_slice(ith); poms=sth_slice(ith)**2 if (z(n1-1)<=zs.and.zs<=z(n2).and.poms>dx**2+dy**2) then ! equidistant grid!!! iph_min_loc=0; iph_max_loc=0; iph_min_loc0=0; lfound=.false. do iph=1,nph_rslice ! loop over all points on a parallel xs=sth_slice(ith)*cph_slice(iph) ys=sth_slice(ith)*sph_slice(iph) if ( x(l1-1)<=xs.and.xs<=x(l2) .and. y(m1-1)<=ys.and.ys<=y(m2) ) then if (iph_min_loc==0) iph_min_loc=iph iph_max_loc=max(iph,iph_max_loc) lfound=.true. if (iph_min_loc0==0) iph_min_loc0=iph elseif (lfound) then if (iph_min_loc>0) iph_max_loc=iph_max_loc+nph_rslice iph_min_loc=0 endif !if (ith<3) & ! write(10+iproc,'(2i3,2f8.4,4i3)') ith,iph,xs,ys,iph_min_loc,iph_max_loc,iph_min_loc0 enddo if (lfound) then ! rank has points on the parallel if (iph_max_loc>nph_rslice) then if (iph_min_loc==0) then iph_min_loc=iph_min_loc0 else iph_min_loc=iph_min_loc-nph_rslice endif iph_max_loc=iph_max_loc-nph_rslice endif if (l0) then iph_min=iph_min_loc; iph_max=iph_max_loc; l0=.false. else iph_min=min(iph_min,iph_min_loc) iph_max=max(iph_max,iph_max_loc) endif if (ith_min==0) ith_min=ith; ith_max=ith endif endif !if (ith<4) & ! write(10+iproc,*) 'ith,iphmm=',ith,iph_min,iph_max enddo ! end loop over r-slice parallels lwrite_slice_r = (ith_min/=0.and.iph_max/=0) if (lwrite_slice_r) then if (allocated(rslice_adjec_corn_inds)) & deallocate(rslice_adjec_corn_inds,rslice_interp_weights) allocate(rslice_adjec_corn_inds(ith_min:ith_max,iph_min:iph_max,3)) allocate(rslice_interp_weights(ith_min:ith_max,iph_min:iph_max,2,2,2)) rslice_adjec_corn_inds=0 do ith=ith_min,ith_max ! loop over all r-slice parallels shared by the rank zs=cth_slice(ith) do iph=iph_min,iph_max if (iph<1) then ! correction for negative indices iphm=nph_rslice+iph else iphm=iph endif ! ! For each point of the r-slice on one of its parallels. ! If the rank shares it, store neighbor indices and calculate weights. ! xs=sth_slice(ith)*cph_slice(iphm) ys=sth_slice(ith)*sph_slice(iphm) do nn=n1,n2 dzs=z(nn)-zs if (dzs>=0.) then do mm=m1,m2 dys=y(mm)-ys if (dys>=0..and.dys=0..and.dxs1e-5) print*, iproc,ith,iph,wsum goto 100 endif enddo endif enddo endif enddo 100 continue enddo enddo !write(20+iproc,'(3(i3,1x)/)') (rslice_adjec_corn_inds(ith_min,ll,:), ll=iph_min,iph_max) endif endsubroutine prep_rslice !*********************************************************************** subroutine alloc_rslice_scal(slice) real, dimension(:,:,:,:,:), allocatable :: slice if (allocated(slice)) deallocate(slice) allocate(slice(ith_min:ith_max,iph_min:iph_max,-1:0,-1:0,-1:0)) slice=0. endsubroutine alloc_rslice_scal !*********************************************************************** subroutine alloc_rslice_scal_2D(slice) real, dimension(:,:), allocatable :: slice if (allocated(slice)) deallocate(slice) allocate(slice(ith_min:ith_max,iph_min:iph_max)) slice=0. endsubroutine alloc_rslice_scal_2D !*********************************************************************** subroutine alloc_rslice_vec(slice,ncomp) use General, only: ioptest real, dimension(:,:,:,:,:,:), allocatable :: slice integer, optional :: ncomp if (allocated(slice)) deallocate(slice) allocate(slice(ith_min:ith_max,iph_min:iph_max,-1:0,-1:0,-1:0,ioptest(ncomp,3))) slice=0. endsubroutine alloc_rslice_vec !*********************************************************************** subroutine store_rslice_scal(pencil,rslice) real, dimension(nx) :: pencil real, dimension(ith_min:ith_max,iph_min:iph_max,-1:0,-1:0,-1:0) :: rslice integer :: ith,iph,mdif,ndif,lind,mind,nind do ith=ith_min,ith_max do iph=iph_min,iph_max lind=rslice_adjec_corn_inds(ith,iph,1) if (lind/=0) then ! rank has point mind=rslice_adjec_corn_inds(ith,iph,2) mdif=m-mind if (mdif==-1.or.mdif==0) then ! m appropriate nind=rslice_adjec_corn_inds(ith,iph,3) ndif=n-nind if (ndif==-1.or.ndif==0) then ! n appropriate lind=lind-nghost ! take relevant part of pencil if (lind==1) then rslice(ith,iph,-1,mdif,ndif)=2.*pencil(lind)-pencil(lind+1) ! extrapolate rslice(ith,iph, 0,mdif,ndif)=pencil(lind) else rslice(ith,iph,:,mdif,ndif)=pencil(lind-1:lind) endif nind=nind-nghost; mind=mind-nghost if (mind==1) & ! m=m1: duplicate in y-direction rslice(ith,iph,:,-1,ndif)=rslice(ith,iph,:,mdif,ndif) if (nind==1) & ! n=n1: duplicate in z-direction rslice(ith,iph,:,mdif,-1)=rslice(ith,iph,:,mdif,ndif) if (nind==1.and.mind==1) & ! m=m1 and n=n1: populate 3rd edge of cubicle rslice(ith,iph,:,-1,-1)=rslice(ith,iph,:,mdif,ndif) endif endif endif enddo enddo endsubroutine store_rslice_scal !*********************************************************************** subroutine store_rslice_vec(pencil,rslice) real, dimension(nx,3) :: pencil real, dimension(ith_min:ith_max,iph_min:iph_max,-1:0,-1:0,-1:0,3) :: rslice integer :: ith,iph,mdif,ndif,lind,mind,nind do ith=ith_min,ith_max do iph=iph_min,iph_max lind=rslice_adjec_corn_inds(ith,iph,1) if (lind/=0) then ! rank has point mind=rslice_adjec_corn_inds(ith,iph,2) mdif=m-mind if (mdif==-1.or.mdif==0) then ! m appropriate nind=rslice_adjec_corn_inds(ith,iph,3) ndif=n-nind if (ndif==-1.or.ndif==0) then ! n appropriate lind=lind-nghost ! take relevant part of pencil if (lind==1) then rslice(ith,iph,-1,mdif,ndif,:)=2.*pencil(lind,:)-pencil(lind+1,:) ! extrapolate rslice(ith,iph, 0,mdif,ndif,:)=pencil(lind,:) else rslice(ith,iph,:,mdif,ndif,:)=pencil(lind-1:lind,:) endif nind=nind-nghost; mind=mind-nghost if (mind==1) & ! m=m1: duplicate in y-direction rslice(ith,iph,:,-1,ndif,:)=rslice(ith,iph,:,mdif,ndif,:) if (nind==1) & ! n=n1: duplicate in z-direction rslice(ith,iph,:,mdif,-1,:)=rslice(ith,iph,:,mdif,ndif,:) if (nind==1.and.mind==1) & ! m=m1 and n=n1: populate 3rd edge of cubicle rslice(ith,iph,:,-1,-1,:)=rslice(ith,iph,:,mdif,ndif,:) endif endif endif enddo enddo endsubroutine store_rslice_vec !*********************************************************************** subroutine interp_rslice(slice,f,ind,rslice) use Messages, only: fatal_error real, dimension(ith_min:ith_max,iph_min:iph_max) :: slice real, dimension(mx,my,mz,*), optional :: f integer, optional :: ind real, dimension(ith_min:ith_max,iph_min:iph_max,-1:0,-1:0,-1:0), optional :: rslice integer :: ith,iph,lind,mind,nind do ith=ith_min,ith_max do iph=iph_min,iph_max lind=rslice_adjec_corn_inds(ith,iph,1) if (lind/=0) then ! proc has point if (present(f)) then mind=rslice_adjec_corn_inds(ith,iph,2) nind=rslice_adjec_corn_inds(ith,iph,3) slice(ith,iph)=sum(f(lind-1:lind,mind-1:mind,nind-1:nind,ind)*rslice_interp_weights(ith,iph,:,:,:)) elseif (present(rslice)) then slice(ith,iph)=sum(rslice(ith,iph,:,:,:)*rslice_interp_weights(ith,iph,:,:,:)) else call fatal_error('interp_rslice', 'either f and ind or rslice must be provided') endif endif enddo enddo endsubroutine interp_rslice !*********************************************************************** function store_dims_rslice() result(dims) integer, dimension(6) :: dims dims=(/nth_rslice,nph_rslice,ith_min,ith_max,iph_min,iph_max/) endfunction store_dims_rslice !*********************************************************************** subroutine write_rslice_position(unit) integer :: unit integer :: iph_min_ if (iph_min==max_int) then iph_min_=0 else iph_min_=iph_min endif write (unit, '(l5,4(1x,i5)," R")') lwrite_slice_r, ith_min,ith_max, iph_min_,iph_max endsubroutine write_rslice_position !*********************************************************************** endmodule Slices_methods