subroutine write_image_vtuCC(qunit,xC,wC,dxC,nPiece,nC1,nC2,nWC,datatype)
! write image data to vtu
use mod_global_parameters
integer, intent(in) :: qunit
integer, intent(in) :: nPiece,nC1,nC2,nWC
double precision, intent(in) :: xC(nPiece,nC1,nC2,2),dxC(nPiece,nc1,nc2,&
2)
double precision, intent(in) :: wC(nPiece,nC1,nC2,nWC)
character(20), intent(in) :: datatype
integer :: nP1,nP2
double precision :: xP(nPiece,nC1+1,nC2+1,2)
integer :: filenr
logical :: fileopen
character (70) :: subname,wname,vname,nameL,nameS
character (len=std_len) :: filename
integer :: ixC1,ixC2,ixP,ix1,ix2,j
integer :: nc,np,icel,VTK_type
integer :: mass
double precision :: logTe
character (30) :: ion
double precision :: line_center
double precision :: spatial_rsl,spectral_rsl,sigma_PSF,wslit
nP1=nC1+1
nP2=nC2+1
np=nP1*nP2
nc=nC1*nC2
! cell corner location
do ixP=1,nPiece
do ix1=1,nP1
do ix2=1,nP2
if (ix1<nP1) xP(ixP,ix1,ix2,1)=xC(ixP,ix1,1,1)-0.5d0*dxC(ixP,ix1,1,&
1)
if (ix1==nP1) xP(ixP,ix1,ix2,1)=xC(ixP,ix1-1,1,1)+0.5d0*dxC(ixP,&
ix1-1,1,1)
if (ix2<nP2) xP(ixP,ix1,ix2,2)=xC(ixP,1,ix2,2)-0.5d0*dxC(ixP,1,ix2,&
2)
if (ix2==nP2) xP(ixP,ix1,ix2,2)=xC(ixP,1,ix2-1,2)+0.5d0*dxC(ixP,1,&
ix2-1,2)
enddo
enddo
enddo
! get information of emission line
if (datatype=='image_euv') then
call get_line_info(wavelength,ion,mass,logTe,line_center,spatial_rsl,&
spectral_rsl,sigma_PSF,wslit)
else if (datatype=='spectrum_euv') then
call get_line_info(spectrum_wl,ion,mass,logTe,line_center,spatial_rsl,&
spectral_rsl,sigma_PSF,wslit)
endif
if (mype==0) then
inquire(qunit,opened=fileopen)
if(.not.fileopen)then
! generate filename
filenr=snapshotini
if (autoconvert) filenr=snapshotnext
if (datatype=='image_euv') then
write(filename,'(a,i4.4,a)') trim(filename_euv),filenr,".vtu"
else if (datatype=='image_sxr') then
write(filename,'(a,i4.4,a)') trim(filename_sxr),filenr,".vtu"
else if (datatype=='image_whitelight') then
write(filename,'(a,i4.4,a)') trim(filename_whitelight),filenr,&
".vtu"
else if (datatype=='spectrum_euv') then
write(filename,'(a,i4.4,a)') trim(filename_spectrum),filenr,".vtu"
endif
open(qunit,file=filename,status='unknown',form='formatted')
endif
! generate xml header
write(qunit,'(a)')'<?xml version="1.0"?>'
write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
write(qunit,'(a)')'<UnstructuredGrid>'
write(qunit,'(a)')'<FieldData>'
write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
'NumberOfTuples="1" format="ascii">'
write(qunit,*) real(global_time*time_convert_factor)
write(qunit,'(a)')'</DataArray>'
if (datatype=='image_euv' .or. datatype=='spectrum_euv') then
write(qunit,'(2a)')'<DataArray type="Float32" Name="logT" ',&
'NumberOfTuples="1" format="ascii">'
write(qunit,*) real(logTe)
write(qunit,'(a)')'</DataArray>'
endif
write(qunit,'(a)')'</FieldData>'
do ixP=1,nPiece
write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
'" NumberOfCells="',nc,'">'
write(qunit,'(a)')'<CellData>'
do j=1,nWC
if (datatype=='image_euv') then
if (j==1) then
if (wavelength<100) then
write(vname,'(a,i2)') "AIA",wavelength
else if (wavelength<1000) then
write(vname,'(a,i3)') "AIA",wavelength
else
write(vname,'(a,i4)') "IRIS",wavelength
endif
endif
if (j==2 .and. dat_resolution) vname='Doppler_velocity'
else if (datatype=='image_sxr') then
if (emin_sxr<10 .and. emax_sxr<10) then
write(vname,'(a,i1,a,i1,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
else if (emin_sxr<10 .and. emax_sxr>=10) then
write(vname,'(a,i1,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
else
write(vname,'(a,i2,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
endif
else if (datatype=='image_whitelight') then
write(vname,'(a)')'whitelight'
else if (datatype=='spectrum_euv') then
if (spectrum_wl==1354) then
write(vname,'(a,i4)') "SG",spectrum_wl
else
write(vname,'(a,i3)') "EIS",spectrum_wl
endif
endif
write(qunit,'(a,a,a)')'<DataArray type="Float64" Name="',&
TRIM(vname),'" format="ascii">'
write(qunit,'(200(1pe14.6))') ((wC(ixP,ixC1,ixC2,j),ixC1=1,nC1),&
ixC2=1,nC2)
write(qunit,'(a)')'</DataArray>'
enddo
write(qunit,'(a)')'</CellData>'
write(qunit,'(a)')'<Points>'
write(qunit,'(a)'&
)'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
do ix2=1,nP2
do ix1=1,nP1
if (datatype=='image_euv' .and. dat_resolution) then
if (LOS_phi==0 .and. LOS_theta==90) then
write(qunit,'(3(1pe14.6))') 0.d0,xP(ixP,ix1,ix2,1),xP(ixP,&
ix1,ix2,2)
else if (LOS_phi==90 .and. LOS_theta==90) then
write(qunit,'(3(1pe14.6))') xP(ixP,ix1,ix2,2),0.d0,xP(ixP,&
ix1,ix2,1)
else
write(qunit,'(3(1pe14.6))') xP(ixP,ix1,ix2,1),xP(ixP,ix1,ix2,&
2),0.d0
endif
else if (datatype=='image_sxr' .and. dat_resolution) then
if (LOS_phi==0 .and. LOS_theta==90) then
write(qunit,'(3(1pe14.6))') 0.d0,xP(ixP,ix1,ix2,1),xP(ixP,&
ix1,ix2,2)
else if (LOS_phi==90 .and. LOS_theta==90) then
write(qunit,'(3(1pe14.6))') xP(ixP,ix1,ix2,2),0.d0,xP(ixP,&
ix1,ix2,1)
else
write(qunit,'(3(1pe14.6))') xP(ixP,ix1,ix2,1),xP(ixP,ix1,ix2,&
2),0.d0
endif
else
write(qunit,'(3(1pe14.6))') xP(ixP,ix1,ix2,1),xP(ixP,ix1,ix2,&
2),0.d0
endif
enddo
enddo
write(qunit,'(a)')'</DataArray>'
write(qunit,'(a)')'</Points>'
! connetivity part
write(qunit,'(a)')'<Cells>'
write(qunit,'(a)'&
)'<DataArray type="Int32" Name="connectivity" format="ascii">'
do ix2=1,nC2
do ix1=1,nC1
write(qunit,'(4(i7))') ix1-1+(ix2-1)*nP1,ix1+(ix2-1)*nP1,&
ix1-1+ix2*nP1,ix1+ix2*nP1
enddo
enddo
write(qunit,'(a)')'</DataArray>'
! offsets data array
write(qunit,'(a)'&
)'<DataArray type="Int32" Name="offsets" format="ascii">'
do icel=1,nc
write(qunit,'(i7)') icel*(2**2)
enddo
write(qunit,'(a)')'</DataArray>'
! VTK cell type data array
write(qunit,'(a)'&
)'<DataArray type="Int32" Name="types" format="ascii">'
! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
VTK_type=8
do icel=1,nc
write(qunit,'(i2)') VTK_type
enddo
write(qunit,'(a)')'</DataArray>'
write(qunit,'(a)')'</Cells>'
write(qunit,'(a)')'</Piece>'
enddo
write(qunit,'(a)')'</UnstructuredGrid>'
write(qunit,'(a)')'</VTKFile>'
close(qunit)
endif
end subroutine write_image_vtuCC