subroutine integrate_EUV_datresol(igrid,nXIF1,nXIF2,xIF1,xIF2,dxIF1,dxIF2,&
fl,EUV,Dpl)
use mod_global_parameters
integer, intent(in) :: igrid,nXIF1,nXIF2
double precision, intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
double precision, intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
type(te_fluid), intent(in) :: fl
double precision, intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
integer :: ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,ixO1,ixO2,&
ixO3,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,ix1,ix2,ix3,i,j
double precision :: xbmin1,xbmin2,xbmin3,xbmax1,xbmax2,xbmax3,xd1,xd2,&
xd3
double precision, allocatable :: flux(:,:,:),v(:,:,:),rho(:,:,:)
double precision, allocatable :: dxb1(:,:,:),dxb2(:,:,:),dxb3(:,:,:)
double precision, allocatable :: EUVg(:,:),Fvg(:,:),xg1(:),xg2(:),&
dxg1(:),dxg2(:)
integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg1,&
iXg2,iXg3
double precision :: EUVt,Fvt,xcmin1,xcmin2,xcmin3,xcmax1,xcmax2,xcmax3,&
xgmin1,xgmin2,xgmin3,xgmax1,xgmax2,xgmax3,r2
integer :: ixPmin1,ixPmin2,ixPmin3,ixPmax1,ixPmax2,ixPmax3,ixP1,ixP2,&
ixP3
integer :: direction_LOS
if (LOS_phi==0 .and. LOS_theta==90) then
direction_LOS=1
else if (LOS_phi==90 .and. LOS_theta==90) then
direction_LOS=2
else
direction_LOS=3
endif
ixOmin1=ixmlo1
ixOmin2=ixmlo2
ixOmin3=ixmlo3
ixOmax1=ixmhi1
ixOmax2=ixmhi2
ixOmax3=ixmhi3
ixImin1=ixglo1
ixImin2=ixglo2
ixImin3=ixglo3
ixImax1=ixghi1
ixImax2=ixghi2
ixImax3=ixghi3
xbmin1=rnode(rpxmin1_,igrid)
xbmin2=rnode(rpxmin2_,igrid)
xbmin3=rnode(rpxmin3_,igrid)
xbmax1=rnode(rpxmax1_,igrid)
xbmax2=rnode(rpxmax2_,igrid)
xbmax3=rnode(rpxmax3_,igrid)
allocate(flux(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
v(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
rho(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3))
allocate(dxb1(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
dxb2(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
dxb3(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3))
dxb1(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)=ps(igrid)%dx(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3,1)
dxb2(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)=ps(igrid)%dx(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3,2)
dxb3(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)=ps(igrid)%dx(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3,3)
! get local EUV flux and velocity
call get_EUV(wavelength,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,ps(igrid)%w,&
ps(igrid)%x,fl,flux)
flux(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)=flux(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)/instrument_resolution_factor**2 !adjust flux due to artifical change of resolution
call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixImin1,ixImin2,ixImin3,ixImax1,&
ixImax2,ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,rho)
v(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3)=-&
ps(igrid)%w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
iw_mom(direction_LOS))/rho(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)
deallocate(rho)
! grid parameters
levelg=ps(igrid)%level
rft=2**(refine_max_level-levelg)
! fine table for storing EUV flux of current grid
select case(direction_LOS)
case(1)
nXg1=ixImax2*rft
nXg2=ixImax3*rft
case(2)
nXg1=ixImax3*rft
nXg2=ixImax1*rft
case(3)
nXg1=ixImax1*rft
nXg2=ixImax2*rft
end select
allocate(EUVg(nXg1,nXg2),Fvg(nXg1,nXg2),xg1(nXg1),xg2(nXg2),dxg1(nXg1),&
dxg2(nXg2))
EUVg=zero
Fvg=zero
xg1=zero
xg2=zero
! integrate for different direction
select case(direction_LOS)
case(1)
do ix2=ixOmin2,ixOmax2
iXgmin1=(ix2-1)*rft+1
iXgmax1=ix2*rft
do ix3=ixOmin3,ixOmax3
iXgmin2=(ix3-1)*rft+1
iXgmax2=ix3*rft
EUVt=0.d0
Fvt=0.d0
do ix1=ixOmin1,ixOmax1
EUVt=EUVt+flux(ix1,ix2,ix3)*dxb1(ix1,ix2,ix3)*unit_length
Fvt=Fvt+flux(ix1,ix2,ix3)*dxb1(ix1,ix2,ix3)*unit_length*v(ix1,&
ix2,ix3)
enddo
EUVg(iXgmin1:iXgmax1,iXgmin2:iXgmax2)=EUVt
Fvg(iXgmin1:iXgmax1,iXgmin2:iXgmax2)=Fvt
enddo
enddo
case(2)
do ix3=ixOmin3,ixOmax3
iXgmin1=(ix3-1)*rft+1
iXgmax1=ix3*rft
do ix1=ixOmin1,ixOmax1
iXgmin2=(ix1-1)*rft+1
iXgmax2=ix1*rft
EUVt=0.d0
Fvt=0.d0
do ix2=ixOmin2,ixOmax2
EUVt=EUVt+flux(ix1,ix2,ix3)*dxb2(ix1,ix2,ix3)*unit_length
Fvt=Fvt+flux(ix1,ix2,ix3)*dxb2(ix1,ix2,ix3)*unit_length*v(ix1,&
ix2,ix3)
enddo
EUVg(iXgmin1:iXgmax1,iXgmin2:iXgmax2)=EUVt
Fvg(iXgmin1:iXgmax1,iXgmin2:iXgmax2)=Fvt
enddo
enddo
case(3)
do ix1=ixOmin1,ixOmax1
iXgmin1=(ix1-1)*rft+1
iXgmax1=ix1*rft
do ix2=ixOmin2,ixOmax2
iXgmin2=(ix2-1)*rft+1
iXgmax2=ix2*rft
EUVt=0.d0
Fvt=0.d0
do ix3=ixOmin3,ixOmax3
EUVt=EUVt+flux(ix1,ix2,ix3)*dxb3(ix1,ix2,ix3)*unit_length
Fvt=Fvt+flux(ix1,ix2,ix3)*dxb3(ix1,ix2,ix3)*unit_length*v(ix1,&
ix2,ix3)
enddo
EUVg(iXgmin1:iXgmax1,iXgmin2:iXgmax2)=EUVt
Fvg(iXgmin1:iXgmax1,iXgmin2:iXgmax2)=Fvt
enddo
enddo
end select
if (SI_unit) then
EUVg=EUVg*1.d2
Fvg=Fvg*1.d2
endif
! mapping grid data to global table
! index ranges in local table
select case(direction_LOS)
case(1)
iXgmin1=(ixOmin2-1)*rft+1
iXgmax1=ixOmax2*rft
iXgmin2=(ixOmin3-1)*rft+1
iXgmax2=ixOmax3*rft
case(2)
iXgmin1=(ixOmin3-1)*rft+1
iXgmax1=ixOmax3*rft
iXgmin2=(ixOmin1-1)*rft+1
iXgmax2=ixOmax1*rft
case(3)
iXgmin1=(ixOmin1-1)*rft+1
iXgmax1=ixOmax1*rft
iXgmin2=(ixOmin2-1)*rft+1
iXgmax2=ixOmax2*rft
end select
! index ranges in global table & mapping
select case(direction_LOS)
case(1)
ixPmin1=(node(pig2_,igrid)-1)*rft*block_nx2+1
ixPmax1=node(pig2_,igrid)*rft*block_nx2
ixPmin2=(node(pig3_,igrid)-1)*rft*block_nx3+1
ixPmax2=node(pig3_,igrid)*rft*block_nx3
case(2)
ixPmin1=(node(pig3_,igrid)-1)*rft*block_nx3+1
ixPmax1=node(pig3_,igrid)*rft*block_nx3
ixPmin2=(node(pig1_,igrid)-1)*rft*block_nx1+1
ixPmax2=node(pig1_,igrid)*rft*block_nx1
case(3)
ixPmin1=(node(pig1_,igrid)-1)*rft*block_nx1+1
ixPmax1=node(pig1_,igrid)*rft*block_nx1
ixPmin2=(node(pig2_,igrid)-1)*rft*block_nx2+1
ixPmax2=node(pig2_,igrid)*rft*block_nx2
end select
xg1(iXgmin1:iXgmax1)=xIF1(ixPmin1:ixPmax1)
xg2(iXgmin2:iXgmax2)=xIF2(ixPmin2:ixPmax2)
dxg1(iXgmin1:iXgmax1)=dxIF1(ixPmin1:ixPmax1)
dxg2(iXgmin2:iXgmax2)=dxIF2(ixPmin2:ixPmax2)
EUV(ixPmin1:ixPmax1,ixPmin2:ixPmax2)=EUV(ixPmin1:ixPmax1,&
ixPmin2:ixPmax2)+EUVg(iXgmin1:iXgmax1,iXgmin2:iXgmax2)
Dpl(ixPmin1:ixPmax1,ixPmin2:ixPmax2)=Dpl(ixPmin1:ixPmax1,&
ixPmin2:ixPmax2)+FVg(iXgmin1:iXgmax1,iXgmin2:iXgmax2)
deallocate(flux,v,dxb1,dxb2,dxb3,EUVg,Fvg,xg1,xg2,dxg1,dxg2)
end subroutine integrate_EUV_datresol