subroutine get_image_datresol(qunit,datatype,fl)
! integrate emission flux along line of sight (LOS)
! in a 3D simulation box and get a 2D EUV image
use mod_global_parameters
use mod_constants
integer, intent(in) :: qunit
character(20), intent(in) :: datatype
type(te_fluid), intent(in) :: fl
double precision :: dx1,dx2,dx3
integer :: numX1,numX2,numX3,ix1,ix2,ix3
double precision, allocatable :: EUV(:,:),EUVs(:,:),Dpl(:,:),Dpls(:,:)
double precision, allocatable :: SXR(:,:),SXRs(:,:),wI(:,:,:)
double precision, allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:),dxIi
integer :: numXI1,numXI2,numSI,numWI
double precision :: xImin1,xImin2,xImin3,xImax1,xImax2,xImax3
integer :: iigrid,igrid,i,j
double precision, allocatable :: xIF1(:),xIF2(:),dxIF1(:),dxIF2(:)
integer :: nXIF1,nXIF2
double precision :: xIFmin1,xIFmin2,xIFmin3,xIFmax1,xIFmax2,xIFmax3
double precision :: unitv,arcsec,RHESSI_rsl
integer :: strtype1,strtype2,strtype3,nstrb1,nstrb2,nstrb3,nbb1,nbb2,&
nbb3,nuni1,nuni2,nuni3,nstr1,nstr2,nstr3,bnx1,bnx2,bnx3
double precision :: qs1,qs2,qs3,dxfirst1,dxfirst2,dxfirst3,dxmid1,dxmid2,&
dxmid3,lenstr1,lenstr2,lenstr3
numX1=domain_nx1*2**(refine_max_level-1)
numX2=domain_nx2*2**(refine_max_level-1)
numX3=domain_nx3*2**(refine_max_level-1)
! parameters for creating table
if (LOS_phi==0 .and. LOS_theta==90) then
nXIF1=domain_nx2*2**(refine_max_level-1)
nXIF2=domain_nx3*2**(refine_max_level-1)
xIFmin1=xprobmin2
xIFmax1=xprobmax2
xIFmin2=xprobmin3
xIFmax2=xprobmax3
bnx1=block_nx2
bnx2=block_nx3
nbb1=domain_nx2
nbb2=domain_nx3
strtype1=stretch_type(2)
strtype2=stretch_type(3)
nstrb1=nstretchedblocks_baselevel(2)
nstrb2=nstretchedblocks_baselevel(3)
qs1=qstretch_baselevel(2)
qs2=qstretch_baselevel(3)
if (mype==0) write(*,'(a)') ' LOS vector: [-1.00 0.00 0.00]'
if (mype==0) write(*,'(a)') ' xI1 vector: [ 0.00 1.00 0.00]'
if (mype==0) write(*,'(a)') ' xI2 vector: [ 0.00 0.00 1.00]'
else if (LOS_phi==90 .and. LOS_theta==90) then
nXIF1=domain_nx3*2**(refine_max_level-1)
nXIF2=domain_nx1*2**(refine_max_level-1)
xIFmin1=xprobmin3
xIFmax1=xprobmax3
xIFmin2=xprobmin1
xIFmax2=xprobmax1
bnx1=block_nx3
bnx2=block_nx1
nbb1=domain_nx3
nbb2=domain_nx1
strtype1=stretch_type(3)
strtype2=stretch_type(1)
nstrb1=nstretchedblocks_baselevel(3)
nstrb2=nstretchedblocks_baselevel(1)
qs1=qstretch_baselevel(3)
qs2=qstretch_baselevel(1)
if (mype==0) write(*,'(a)') ' LOS vector: [ 0.00 -1.00 0.00]'
if (mype==0) write(*,'(a)') ' xI1 vector: [-1.00 0.00 0.00]'
if (mype==0) write(*,'(a)') ' xI2 vector: [ 0.00 0.00 1.00]'
else
nXIF1=domain_nx1*2**(refine_max_level-1)
nXIF2=domain_nx2*2**(refine_max_level-1)
xIFmin1=xprobmin1
xIFmax1=xprobmax1
xIFmin2=xprobmin2
xIFmax2=xprobmax2
bnx1=block_nx1
bnx2=block_nx2
nbb1=domain_nx1
nbb2=domain_nx2
strtype1=stretch_type(1)
strtype2=stretch_type(2)
nstrb1=nstretchedblocks_baselevel(1)
nstrb2=nstretchedblocks_baselevel(2)
qs1=qstretch_baselevel(1)
qs2=qstretch_baselevel(2)
if (mype==0) write(*,'(a)') ' LOS vector: [ 0.00 0.00 -1.00]'
if (mype==0) write(*,'(a)') ' xI1 vector: [ 1.00 0.00 0.00]'
if (mype==0) write(*,'(a)') ' xI2 vector: [ 0.00 1.00 0.00]'
endif
allocate(xIF1(nXIF1),xIF2(nXIF2),dxIF1(nXIF1),dxIF2(nXIF2))
! initialize image coordinate
select case(strtype1)
case(0) ! uniform
dxIF1(:)=(xIFmax1-xIFmin1)/nXIF1
do ix1=1,nXIF1
xIF1(ix1)=xIFmin1+dxIF1(ix1)*(ix1-half)
enddo
case(1) ! uni stretch
qs1=qs1**(one/2**(refine_max_level-1))
dxfirst1=(xIFmax1-xIFmin1)*(one-qs1)/(one-qs1**nXIF1)
dxIF1(1)=dxfirst1
do ix1=2,nXIF1
dxIF1(ix1)=dxfirst1*qs1**(ix1-1)
xIF1(ix1)=dxIF1(1)/(one-qs1)*(one-qs1**(ix1-1))+half*dxIF1(ix1)
enddo
case(2) ! symm stretch
! base level, nbb = nstr + nuni + nstr
nstr1=nstrb1*bnx1/2
nuni1=nbb1-nstrb1*bnx1
lenstr1=(xIFmax1-xIFmin1)/(2.d0+nuni1*(one-qs1)/(one-qs1**nstr1))
dxfirst1=(xIFmax1-xIFmin1)/(dble(nuni1)+&
2.d0/(one-qs1)*(one-qs1**nstr1))
dxmid1=dxfirst1
! refine_max level, numXI = nstr + nuni + nstr
nstr1=nstr1*2**(refine_max_level-1)
nuni1=nuni1*2**(refine_max_level-1)
qs1=qs1**(one/2**(refine_max_level-1))
dxfirst1=lenstr1*(one-qs1)/(one-qs1**nstr1)
dxmid1=dxmid1/2**(refine_max_level-1)
! uniform center
if(nuni1 .gt. 0) then
do ix1=nstr1+1,nstr1+nuni1
dxIF1(ix1)=dxmid1
xIF1(ix1)=lenstr1+(dble(ix1)-0.5d0-nstr1)*dxIF1(ix1)+xIFmin1
enddo
endif
! left half
do ix1=nstr1,1,-1
dxIF1(ix1)=dxfirst1*qs1**(nstr1-ix1)
xIF1(ix1)=xIFmin1+lenstr1-dxIF1(ix1)*half-dxfirst1*(one-qs1**(nstr1-&
ix1))/(one-qs1)
enddo
! right half
do ix1=nstr1+nuni1+1,nXIF1
dxIF1(ix1)=dxfirst1*qs1**(ix1-nstr1-nuni1-1)
xIF1(ix1)=xIFmax1-lenstr1+dxIF1(ix1)*half+&
dxfirst1*(one-qs1**(ix1-nstr1-nuni1-1))/(one-qs1)
enddo
case default
call mpistop("unknown stretch type")
end select
select case(strtype2)
case(0) ! uniform
dxIF2(:)=(xIFmax2-xIFmin2)/nXIF2
do ix2=1,nXIF2
xIF2(ix2)=xIFmin2+dxIF2(ix2)*(ix2-half)
enddo
case(1) ! uni stretch
qs2=qs2**(one/2**(refine_max_level-1))
dxfirst2=(xIFmax2-xIFmin2)*(one-qs2)/(one-qs2**nXIF2)
dxIF2(1)=dxfirst2
do ix2=2,nXIF1
dxIF2(ix2)=dxfirst2*qs2**(ix2-1)
xIF2(ix2)=dxIF2(1)/(one-qs2)*(one-qs2**(ix2-1))+half*dxIF2(ix2)
enddo
case(2) ! symm stretch
! base level, nbb = nstr + nuni + nstr
nstr2=nstrb2*bnx2/2
nuni2=nbb2-nstrb2*bnx2
lenstr2=(xIFmax2-xIFmin2)/(2.d0+nuni2*(one-qs2)/(one-qs2**nstr2))
dxfirst2=(xIFmax2-xIFmin2)/(dble(nuni2)+&
2.d0/(one-qs2)*(one-qs2**nstr2))
dxmid2=dxfirst2
! refine_max level, numXI = nstr + nuni + nstr
nstr2=nstr2*2**(refine_max_level-1)
nuni2=nuni2*2**(refine_max_level-1)
qs2=qs2**(one/2**(refine_max_level-1))
dxfirst2=lenstr2*(one-qs2)/(one-qs2**nstr2)
dxmid2=dxmid2/2**(refine_max_level-1)
! uniform center
if(nuni2 .gt. 0) then
do ix2=nstr2+1,nstr2+nuni2
dxIF2(ix2)=dxmid2
xIF2(ix2)=lenstr2+(dble(ix2)-0.5d0-nstr2)*dxIF2(ix2)+xIFmin2
enddo
endif
! left half
do ix2=nstr2,1,-1
dxIF2(ix2)=dxfirst2*qs2**(nstr2-ix2)
xIF2(ix2)=xIFmin2+lenstr2-dxIF2(ix2)*half-dxfirst2*(one-qs2**(nstr2-&
ix2))/(one-qs2)
enddo
! right half
do ix2=nstr2+nuni2+1,nXIF2
dxIF2(ix2)=dxfirst2*qs2**(ix2-nstr2-nuni2-1)
xIF2(ix2)=xIFmax2-lenstr2+dxIF2(ix2)*half+&
dxfirst2*(one-qs2**(ix2-nstr2-nuni2-1))/(one-qs2)
enddo
case default
call mpistop("unknown stretch type")
end select
! integrate EUV flux and get cell average flux for image
if (datatype=='image_euv') then
if (SI_unit) then
unitv=unit_velocity/1.0e3 ! km/s
else
unitv=unit_velocity/1.0e5 ! km/s
endif
numWI=2
allocate(wI(nXIF1,nXIF2,numWI))
allocate(EUVs(nXIF1,nXIF2),EUV(nXIF1,nXIF2))
allocate(Dpl(nXIF1,nXIF2),Dpls(nXIF1,nXIF2))
EUVs=0.0d0
EUV=0.0d0
Dpl=0.d0
Dpls=0.d0
do iigrid=1,igridstail; igrid=igrids(iigrid);
call integrate_EUV_datresol(igrid,nXIF1,nXIF2,xIF1,xIF2,dxIF1,dxIF2,&
fl,EUVs,Dpls)
enddo
numSI=nXIF1*nXIF2
call MPI_ALLREDUCE(EUVs,EUV,numSI,MPI_DOUBLE_PRECISION, MPI_SUM,icomm,&
ierrmpi)
call MPI_ALLREDUCE(Dpls,Dpl,numSI,MPI_DOUBLE_PRECISION, MPI_SUM,icomm,&
ierrmpi)
do ix1=1,nXIF1
do ix2=1,nXIF2
if (EUV(ix1,ix2)<smalldouble) EUV(ix1,ix2)=zero
if(EUV(ix1,ix2)/=0) then
Dpl(ix1,ix2)=(Dpl(ix1,ix2)/EUV(ix1,ix2))*unitv
else
Dpl(ix1,ix2)=0.d0
endif
if (abs(Dpl(ix1,ix2))<smalldouble) Dpl(ix1,ix2)=zero
enddo
enddo
wI(:,:,1)=EUV(:,:)
wI(:,:,2)=Dpl(:,:)
call output_data(qunit,xIF1,xIF2,dxIF1,dxIF2,wI,nXIF1,nXIF2,numWI,&
datatype)
deallocate(WI,EUV,EUVs,Dpl,Dpls)
endif
! integrate SXR flux and get cell average flux for image
if (datatype=='image_sxr') then
if (SI_unit) then
arcsec=7.25d5
else
arcsec=7.25d7
endif
RHESSI_rsl=2.3d0/instrument_resolution_factor
numWI=1
allocate(wI(nXIF1,nXIF2,numWI))
allocate(SXRs(nXIF1,nXIF2),SXR(nXIF1,nXIF2))
SXRs=0.0d0
SXR=0.0d0
do iigrid=1,igridstail; igrid=igrids(iigrid);
call integrate_SXR_datresol(igrid,nXIF1,nXIF2,xIF1,xIF2,dxIF1,dxIF2,&
fl,SXRs)
enddo
numSI=nXIF1*nXIF2
call MPI_ALLREDUCE(SXRs,SXR,numSI,MPI_DOUBLE_PRECISION, MPI_SUM,icomm,&
ierrmpi)
SXR=SXR*(RHESSI_rsl*arcsec)**2 !photons cm-2,cm-2,cm-2 s-1,s-1,s-1 pixel-1,pixel-1,pixel-1
do ix1=1,nXIF1
do ix2=1,nXIF2
if (SXR(ix1,ix2)<smalldouble) SXR(ix1,ix2)=zero
enddo
enddo
wI(:,:,1)=SXR(:,:)
call output_data(qunit,xIF1,xIF2,dxIF1,dxIF2,wI,nXIF1,nXIF2,numWI,&
datatype)
deallocate(WI,SXR,SXRs)
endif
deallocate(xIF1,xIF2,dxIF1,dxIF2)
end subroutine get_image_datresol