subroutine get_spectrum_datresol(qunit,datatype,fl)
integer, intent(in) :: qunit
character(20), intent(in) :: datatype
type(te_fluid), intent(in) :: fl
integer :: numWL,numXS,iwL,ixS,numWI,numS
double precision :: dwLg,xSmin,xSmax,wLmin,wLmax
double precision, allocatable :: wL(:),xS(:),dwL(:),dxS(:)
double precision, allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
integer :: strtype,nstrb,nbb,nuni,nstr,bnx
double precision :: qs,dxfirst,dxmid,lenstr
integer :: iigrid,igrid,j,dir_loc
double precision :: xbmin(1:ndim),xbmax(1:ndim)
dwLg=1.d-3
numWL=4*int((spectrum_window_max-spectrum_window_min)/(4.d0*dwLg))
wLmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwLg*numWL/2
wLmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwLg*numWL/2
allocate(wL(numWL),dwL(numWL))
dwL(:)=dwLg
do iwL=1,numWL
wL(iwL)=wLmin+iwL*dwLg-half*dwLg
enddo
select case(direction_slit)
case (1)
numXS=domain_nx1*2**(refine_max_level-1)
xSmin=xprobmin1
xSmax=xprobmax1
bnx=block_nx1
nbb=domain_nx1
strtype=stretch_type(1)
nstrb=nstretchedblocks_baselevel(1)
qs=qstretch_baselevel(1)
if (mype==0) print *, 'Direction of the slit: x'
case (2)
numXS=domain_nx2*2**(refine_max_level-1)
xSmin=xprobmin2
xSmax=xprobmax2
bnx=block_nx2
nbb=domain_nx2
strtype=stretch_type(2)
nstrb=nstretchedblocks_baselevel(2)
qs=qstretch_baselevel(2)
if (mype==0) print *, 'Direction of the slit: y'
case (3)
numXS=domain_nx3*2**(refine_max_level-1)
xSmin=xprobmin3
xSmax=xprobmax3
bnx=block_nx3
nbb=domain_nx3
strtype=stretch_type(3)
nstrb=nstretchedblocks_baselevel(3)
qs=qstretch_baselevel(3)
if (mype==0) print *, 'Direction of the slit: z'
case default
call MPISTOP('Wrong direction_slit')
end select
allocate(xS(numXS),dxS(numXS),spectra(numWL,numXS),spectra_rc(numWL,&
numXS))
numWI=1
allocate(wI(numWL,numXS,numWI))
select case(strtype)
case(0) ! uniform
dxS(:)=(xSmax-xSmin)/numXS
do ixS=1,numXS
xS(ixS)=xSmin+dxS(ixS)*(ixS-half)
enddo
case(1) ! uni stretch
qs=qs**(one/2**(refine_max_level-1))
dxfirst=(xSmax-xSmin)*(one-qs)/(one-qs**numXS)
dxS(1)=dxfirst
do ixS=2,numXS
dxS(ixS)=dxfirst*qs**(ixS-1)
xS(ixS)=dxS(1)/(one-qs)*(one-qs**(ixS-1))+half*dxS(ixS)
enddo
case(2) ! symm stretch
! base level, nbb = nstr + nuni + nstr
nstr=nstrb*bnx/2
nuni=nbb-nstrb*bnx
lenstr=(xSmax-xSmin)/(2.d0+nuni*(one-qs)/(one-qs**nstr))
dxfirst=(xSmax-xSmin)/(dble(nuni)+2.d0/(one-qs)*(one-qs**nstr))
dxmid=dxfirst
! refine_max level, numXI = nstr + nuni + nstr
nstr=nstr*2**(refine_max_level-1)
nuni=nuni*2**(refine_max_level-1)
qs=qs**(one/2**(refine_max_level-1))
dxfirst=lenstr*(one-qs)/(one-qs**nstr)
dxmid=dxmid/2**(refine_max_level-1)
! uniform center
if(nuni .gt. 0) then
do ixS=nstr+1,nstr+nuni
dxS(ixS)=dxmid
xS(ixS)=lenstr+(dble(ixS)-0.5d0-nstr)*dxS(ixS)+xSmin
enddo
endif
! left half
do ixS=nstr,1,-1
dxS(ixS)=dxfirst*qs**(nstr-ixS)
xS(ixS)=xSmin+lenstr-dxS(ixS)*half-dxfirst*(one-qs**(nstr-ixS))/(one-&
qs)
enddo
! right half
do ixS=nstr+nuni+1,numXS
dxS(ixS)=dxfirst*qs**(ixS-nstr-nuni-1)
xS(ixS)=xSmax-lenstr+dxS(ixS)*half+&
dxfirst*(one-qs**(ixS-nstr-nuni-1))/(one-qs)
enddo
case default
call mpistop("unknown stretch type")
end select
if (LOS_phi==0 .and. LOS_theta==90 .and. direction_slit==2) then
! LOS->x slit->y
dir_loc=3
else if (LOS_phi==0 .and. LOS_theta==90 .and. direction_slit==3) then
! LOS->x slit->z
dir_loc=2
else if (LOS_phi==90 .and. LOS_theta==90 .and. direction_slit==1) then
! LOS->y slit->x
dir_loc=3
else if (LOS_phi==90 .and. LOS_theta==90 .and. direction_slit==3) then
! LOS->y slit->z
dir_loc=1
else if (LOS_theta==0 .and. direction_slit==1) then
! LOS->z slit->x
dir_loc=2
else if (LOS_theta==0 .and. direction_slit==2) then
! LOS->z slit->y
dir_loc=1
else
call MPISTOP('Wrong combination of LOS and slit direction!')
endif
if (dir_loc==1) then
if (location_slit>xprobmax1 .or. location_slit<xprobmin1) then
call MPISTOP('Wrong value for location_slit!')
endif
if(mype==0) write(*,'(a,f8.1,a)') ' Location of slit: x = ',&
location_slit,' Unit_length'
else if (dir_loc==2) then
if (location_slit>xprobmax2 .or. location_slit<xprobmin2) then
call MPISTOP('Wrong value for location_slit!')
endif
if(mype==0) write(*,'(a,f8.1,a)') ' Location of slit: y = ',&
location_slit,' Unit_length'
else
if (location_slit>xprobmax3 .or. location_slit<xprobmin3) then
call MPISTOP('Wrong value for location_slit!')
endif
if(mype==0) write(*,'(a,f8.1,a)') ' Location of slit: z = ',&
location_slit,' Unit_length'
endif
! find slit and do integration
spectra=zero
do iigrid=1,igridstail; igrid=igrids(iigrid);
xbmin(1)=rnode(rpxmin1_,igrid);xbmin(2)=rnode(rpxmin2_,igrid)
xbmin(3)=rnode(rpxmin3_,igrid);
xbmax(1)=rnode(rpxmax1_,igrid);xbmax(2)=rnode(rpxmax2_,igrid)
xbmax(3)=rnode(rpxmax3_,igrid);
if (location_slit>=xbmin(dir_loc) .and. location_slit<xbmax(dir_loc)) &
then
call integrate_spectra_datresol(igrid,wL,dwL,spectra,numWL,numXS,&
dir_loc,fl)
endif
enddo
numS=numWL*numXS
call MPI_ALLREDUCE(spectra,spectra_rc,numS,MPI_DOUBLE_PRECISION, MPI_SUM,&
icomm,ierrmpi)
do iwL=1,numWL
do ixS=1,numXS
if (spectra_rc(iwL,ixS)>smalldouble) then
wI(iwL,ixS,1)=spectra_rc(iwL,ixS)
else
wI(iwL,ixS,1)=zero
endif
enddo
enddo
call output_data(qunit,wL,xS,dwL,dxS,wI,numWL,numXS,numWI,datatype)
deallocate(wL,xS,dwL,dxS,spectra,spectra_rc,wI)
end subroutine get_spectrum_datresol