This subroutine in problem.f90 builds a data.curve file for each frame and prints to that file the q array, as well as the exact solution for the sod shock tube.

 SUBROUTINE ProblemAfterStep(Info)
     TYPE(InfoDef) :: Info
     INTEGER :: i
     REAL(KIND=qPREC), ALLOCATABLE, DIMENSION(:,:) :: qExact
     REAL(KIND=qPREC) :: um, s, max_speed
     REAL(KIND=qPREC), ALLOCATABLE, DIMENSION(:) :: wmiddle
     CHARACTER(50) :: S_FILENAME

     IF (MPI_ID == MASTER .OR. lParallelHDF5) THEN
        !write(*,*) 'levels%tnow=', levels(maxlevel)%step                                                                                                         
        WRITE(s_filename, '(A8,I5.5,A6)') 'out/data', current_frame, '.curve'  !write to the character var 's_filename' the concatenated output list              
        !PRINT *, 's_filename=', s_filename                                                                                                                       
        !PRINT *, 'current_time=', levels(maxlevel)%tnow !gives the current time                                                                                  
        !RETURN                                                                                                                                                   
        OPEN(UNIT=11, FILE=s_filename, status='replace') !open the filename defined above for each current frame                                                  
        write(11, *) '# rho'
        DO i=1, Info%mX(1)
           IF (Info%mX(2) /= Info%mX(1)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx, Info%q(i,1,1,1)
           IF (Info%mX(1)==Info%mX(2)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx*sqrt(2d0), Info%q(i,i,1,1)
        END DO
        write(11,*)
        write(11,*)
        write(11,*) '# vx'
        DO i=1, Info%mX(1)
           IF (Info%mX(2) /= Info%mX(1)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx, Info%q(i,1,1,ivx)/Info%q(i,1,1,1)
           IF (Info%mX(1)==Info%mX(2)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx*sqrt(2d0), Info%q(i,i,1,ivx)/Info%q(i,i,1,1)
        END DO
        write(11,*)
        write(11,*)
        write(11,*) '# P'
        DO i=1, Info%mX(1)
           IF (Info%mX(2) /= Info%mX(1)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx, (Press(Info%q(i,1,1,:)))
           IF (Info%mX(1)==Info%mX(2)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx*sqrt(2d0), (Press(Info%q(i,i,1,:)))
        END DO
        write(11,*)
        write(11,*)

        IF (lMHD) THEN
           write(11,*) '# vy'
           DO i=1, Info%mX(1)
              IF (Info%mX(2) /= Info%mX(1)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx, Info%q(i,1,1,ivy)/Info%q(i,1,1,1)
              IF (Info%mX(1)==Info%mX(2)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx*sqrt(2d0), Info%q(i,i,1,ivy)/Info%q(i,i,1,1)
           END DO
           write(11,*)
           write(11,*)
           write(11,*) '# vz'
           DO i=1, Info%mX(1)
              IF (Info%mX(2) /= Info%mX(1)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx, Info%q(i,1,1,ivz)/Info%q(i,1,1,1)
              IF (Info%mX(1)==Info%mX(2)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx*sqrt(2d0), Info%q(i,i,1,ivz)/Info%q(i,i,1,1)
           END DO
           write(11,*)
           write(11,*)
           write(11,*) '# Bx'
           DO i=1, Info%mX(1)
              IF (Info%mX(2) /= Info%mX(1)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx, Info%q(i,1,1,iBx)
              IF (Info%mX(1)==Info%mX(2)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx*sqrt(2d0), Info%q(i,i,1,iBx)
           END DO
           write(11,*)
           write(11,*)
           write(11,*) '# By'
           DO i=1, Info%mX(1)
              IF (Info%mX(2) /= Info%mX(1)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx, Info%q(i,1,1,iBy)
              IF (Info%mX(1)==Info%mX(2)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx*sqrt(2d0), Info%q(i,i,1,iBy)
           END DO
           write(11,*)
           write(11,*)
           write(11,*) '# Bz'
           DO i=1, Info%mX(1)
              IF (Info%mX(2) /= Info%mX(1)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx, Info%q(i,1,1,iBz)
              IF (Info%mX(1)==Info%mX(2)) write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx*sqrt(2d0), Info%q(i,i,1,iBz)
           END DO
           write(11,*)
           write(11,*)
        ELSE
           ALLOCATE(qExact(Info%mX(1), NrHydroVars))
           ALLOCATE(wmiddle(1:NrHydroVars))
           DO i=1, Info%mX(1)
              !       S=(Info%xBounds(1,1)+(REAL(i)-half)*Levels(Info%level)%dx-position(1))/final_time !get the exact solution for the final time                
              S=(Info%xBounds(1,1)+(REAL(i)-half)*Levels(Info%level)%dx-position(1))/levels(maxlevel)%tnow !get the exact solution for the current time           
              IF (Info%mX(1)==Info%mX(2)) S=S*sqrt(2d0)
              CALL vacuum_solve(qabove((/1,3,2/)), qbelow((/1,3,2/)), wmiddle, um, s, max_speed)
              qExact(i,:)=wmiddle
           END DO
           write(11,*) '# rho_Exact'
           DO i=1, Info%mX(1)
              write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx, qExact(i,1)
           END DO
           write(11,*)
           write(11,*)
           write(11,*) '# vx_Exact'
           DO i=1, Info%mX(1)
              write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx, qExact(i,3)
           END DO
           write(11,*)
           write(11,*)
           write(11,*) '# P_Exact'
           DO i=1, Info%mX(1)
              write(11,*) Info%xbounds(1,1)+(REAL(i)-half)*levels(Info%level)%dx, qExact(i,2)
           END DO
           CLOSE(11)
           IF (Info%mX(2) /= Info%mX(1)) write(*,*) 'dx, L2 Norm=', levels(Info%level)%dx, sum(abs(qExact(:,1)-Info%q(1:Info%mX(1),1,1,1)))/Info%mX(1)
           IF (Info%mX(1)==Info%mX(2)) write(*,*) 'dx, L2 Norm=', levels(Info%level)%dx, sum(abs(qExact(:,1)-(/(Info%q(i,i,1,1), i=1,Info%mX(1))/)))/Info%mX(1)
        END IF
     END IF

   END SUBROUTINE ProblemAfterStep