riemann2d.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: riemann2d.f90                                                     #
!#                                                                           #
!# Copyright (C) 2006-2024                                                   #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!#                                                                           #
!# This program is free software; you can redistribute it and/or modify      #
!# it under the terms of the GNU General Public License as published by      #
!# the Free Software Foundation; either version 2 of the License, or (at     #
!# your option) any later version.                                           #
!#                                                                           #
!# This program is distributed in the hope that it will be useful, but       #
!# WITHOUT ANY WARRANTY; without even the implied warranty of                #
!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or        #
!# NON INFRINGEMENT.  See the GNU General Public License for more            #
!# details.                                                                  #
!#                                                                           #
!# You should have received a copy of the GNU General Public License         #
!# along with this program; if not, write to the Free Software               #
!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 #
!#                                                                           #
!#############################################################################
PROGRAM riemann2d
  USE fosite_mod
#include "tap.h"
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! simulation parameter
  INTEGER, PARAMETER :: select_test = 6   ! run this test, should be in 1..19
  LOGICAL, PARAMETER :: run_all = .true.  ! set this if you want to run all 19 tests
  REAL, PARAMETER    :: gamma = 1.4        ! ratio of specific heats
  ! mesh settings
  INTEGER, PARAMETER :: mgeo = cartesian   ! geometry of the mesh
!   INTEGER, PARAMETER :: MGEO = CYLINDRICAL
!   INTEGER, PARAMETER :: MGEO = LOGCYLINDRICAL
  INTEGER, PARAMETER :: xres = 100         ! resolution
  INTEGER, PARAMETER :: yres = 100
  INTEGER, PARAMETER :: zres = 1
  REAL, PARAMETER    :: rmin = 1.0e-4      ! inner radius for polar grids
  ! output file parameter
  INTEGER, PARAMETER :: onum = 1           ! number of output data sets
  CHARACTER(LEN=256), PARAMETER &          ! output data dir
                     :: odir = './'
  CHARACTER(LEN=256), PARAMETER &          ! output data file name
                     :: ofname = 'riemann2d'
  ! do not change these parameters
  INTEGER, PARAMETER :: num_tests = 19     ! total number of 2D setups
  REAL, PARAMETER, DIMENSION(NUM_TESTS) :: & ! stop time for each test
                    test_stoptime = (/ 0.2,0.2,0.3,0.25,0.23,0.3,0.25, &
                                       0.25,0.3,0.15,0.3,0.25,0.3,0.1, &
                                       0.2,0.2,0.3,0.2,0.3/)
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE      :: sim
  LOGICAL                         :: ok(num_tests), mpifinalize = .false.
  INTEGER                         :: ic,err
  !--------------------------------------------------------------------------!
 
  IF (run_all) THEN
tap_plan(num_tests)
  ELSE
tap_plan(1)
  END IF
 
  DO ic=1,num_tests
    IF (.NOT.run_all.AND.ic.NE.select_test) cycle ! skip if not selected
 
    IF (ALLOCATED(sim)) DEALLOCATE(sim)
    ALLOCATE(sim,stat=err)
    IF (err.NE.0) &
      CALL sim%Error("riemann2d","cannot allocate memory for Sim")
 
    CALL sim%InitFosite()
    CALL makeconfig(sim, sim%config, ic)
 
!  CALL PrintDict(config)
 
    CALL sim%Setup()
 
    ! set initial condition
    CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc, ic)
    CALL sim%Run()
    ok(ic) = .NOT.sim%aborted
 
    IF ((run_all.AND.ic.EQ.num_tests).OR.(.NOT.run_all)) mpifinalize=.true.
    CALL sim%Finalize(mpifinalize)
 
  END DO
 
  DO ic=1,num_tests
    IF (.NOT.run_all.AND.ic.NE.select_test) cycle ! skip if not selected
tap_check(ok(ic),"Simulation finished")
  END DO
 
  IF (ALLOCATED(sim)) DEALLOCATE(sim)
 
tap_done
 
CONTAINS
 
  SUBROUTINE makeconfig(Sim, config, ic)
    USE functions, ONLY : asinh
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fosite)     :: Sim
    TYPE(dict_typ),POINTER :: config
    INTEGER           :: ic
    !------------------------------------------------------------------------!
    ! Local variable declaration
    INTEGER           :: bc(6)
    TYPE(dict_typ), POINTER :: mesh, physics, boundary, datafile, &
                               timedisc, fluxes
    REAL              :: x1,x2,y1,y2,z1,z2,sc
    CHARACTER(LEN=16) :: fext
    !------------------------------------------------------------------------!
    INTENT(INOUT)     :: sim
    INTENT(IN)        :: ic
    CHARACTER(LEN=8)  :: geo_str
    !------------------------------------------------------------------------!
    IF (ic.LT.1.OR.ic.GT.num_tests) &
       CALL sim%Error("riemann2d::MakeConfig","invalid test number selected")
 
    ! mesh settings and boundary conditions
    SELECT CASE(mgeo)
    CASE(cartesian)
       sc = 1.0
       x1 = -0.5
       x2 =  0.5
       y1 = -0.5
       y2 =  0.5
       z1 = -0.0
       z2 =  0.0
       bc(west)  = no_gradients
       bc(east)  = no_gradients
       bc(south) = no_gradients
       bc(north) = no_gradients
       bc(bottom)= no_gradients
       bc(top)   = no_gradients
       geo_str = "-cart"
   CASE(cylindrical)
      sc = 1.0
      x1 = 0.0
      x2 = 0.5*sqrt(2.0)
      y1 = 0.0
      y2 = 2*pi
      z1 = -0.0
      z2 = 0.0
      bc(west)   = absorbing
      bc(east)   = no_gradients
      bc(south)  = periodic
      bc(north)  = periodic
      bc(bottom) = reflecting
      bc(top)    = reflecting
      geo_str = "-cyl"
   CASE(logcylindrical)
      sc = 0.3
      x1 = log(rmin/sc)
      x2 = log(0.5*sqrt(2.0)/sc)
      y1 = -pi
      y2 = pi
      z1 = 0.0
      z2 = 0.0
      bc(west)   = axis
      bc(east)   = no_gradients
      bc(south)  = periodic
      bc(north)  = periodic
      bc(bottom) = no_gradients
      bc(top)    = no_gradients
      geo_str = "-logcyl"
    CASE DEFAULT
      CALL sim%Error("riemann2d::MakeConfig","geometry not supported in this test")
    END SELECT
 
    ! mesh settings
    mesh => dict("meshtype" / midpoint, &
           "geometry" / mgeo, &
           "inum"     / xres, &
           "jnum"     / yres, &
           "knum"     / zres, &
           "xmin"     / x1, &
           "xmax"     / x2, &
           "ymin"     / y1, &
           "ymax"     / y2, &
           "zmin"     / z1, &
           "zmax"     / z2, &
           "output/commutator" / 1,&
           "gparam"   / sc)
 
    ! boundary conditions
    boundary => dict("western" / bc(west), &
               "eastern" / bc(east), &
               "southern" / bc(south), &
               "northern" / bc(north), &
               "bottomer" / bc(bottom),&
               "topper"   / bc(top))
 
    ! physics settings
    physics => dict("problem" / euler, &
              "gamma"   / gamma)         ! ratio of specific heats        !
 
    ! flux calculation and reconstruction method
    fluxes => dict("order"     / linear, &
             "fluxtype"  / kt, &
             "variables" / conservative, &        ! vars. to use for reconstruction!
             "limiter"   / monocent, &    ! one of: minmod, monocent,...   !
             "theta"     / 1.2)          ! optional parameter for limiter !
 
    ! time discretization settings
    timedisc => dict( &
           "method"   / modified_euler, &
           "order"    / 3, &
           "cfl"      / 0.4, &
           "stoptime" / test_stoptime(ic), &
           "dtlimit"  / 1.0e-10, &
           "maxiter"  / 10000000, &
           "output/geometrical_sources" / 1)
 
    ! initialize data input/output
    ! append test number to file names
    WRITE (fext, '(A,I2.2)') trim(geo_str) // "_", ic
    datafile => dict( &
!                "fileformat" / GNUPLOT, "filecycles" / 0, &
                "fileformat" / vtk, &
!                 "fileformat" / XDMF, &
                "filename"   / (trim(odir) // trim(ofname) // trim(fext)), &
                "count"      / onum)
 
    config => dict("mesh" / mesh, &
             "physics"  / physics, &
             "boundary" / boundary, &
             "fluxes"   / fluxes, &
             "timedisc" / timedisc, &
!             "logfile"  / logfile, &
             "datafile" / datafile)
  END SUBROUTINE makeconfig
 
  SUBROUTINE initdata(Mesh,Physics,Timedisc,ic)
  IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(physics_base) :: Physics
    CLASS(mesh_base)    :: Mesh
    CLASS(timedisc_base):: Timedisc
    INTEGER           :: ic
    !------------------------------------------------------------------------!
    ! Local variable declaration
    TYPE(marray_base) :: vcart,vcurv
    REAL              :: xmin,ymin,xmax,ymax,x0,y0
#ifdef PARALLEL
    REAL              :: xmin_all,xmax_all,ymin_all,ymax_all
    INTEGER           :: ierr
#endif
    REAL, DIMENSION(NUM_TESTS,4):: den,pre,vx,vy
    CHARACTER(LEN=64) :: teststr
    !------------------------------------------------------------------------!
    INTENT(IN)        :: mesh,physics,ic
    INTENT(INOUT)     :: timedisc
    !------------------------------------------------------------------------!
    ! minima and maxima of _cartesian_ coordinates
    xmin = minval(mesh%bccart(:,:,:,1))
    ymin = minval(mesh%bccart(:,:,:,2))
    xmax = maxval(mesh%bccart(:,:,:,1))
    ymax = maxval(mesh%bccart(:,:,:,2))
#ifdef PARALLEL
    CALL mpi_allreduce(xmin,xmin_all,1,default_mpi_real,mpi_min,mesh%comm_cart,ierr)
    xmin = xmin_all
    CALL mpi_allreduce(ymin,ymin_all,1,default_mpi_real,mpi_min,mesh%comm_cart,ierr)
    ymin = ymin_all
    CALL mpi_allreduce(xmax,xmax_all,1,default_mpi_real,mpi_max,mesh%comm_cart,ierr)
    xmax = xmax_all
    CALL mpi_allreduce(ymax,ymax_all,1,default_mpi_real,mpi_max,mesh%comm_cart,ierr)
    ymax = ymax_all
#endif
    x0 = xmin + 0.5*abs(xmax-xmin)
    y0 = ymin + 0.5*abs(ymax-ymin)
 
    vcart = marray_base(3)
    vcurv = marray_base(3)
 
    ! set defaults
    den(:,:) = 1.
    vx(:,:)  = 0.
    vy(:,:)  = 0.
    pre(:,:) = 1.
 
    IF (ic.LT.1.OR.ic.GT.num_tests) &
      CALL mesh%Error("InitData","Sorry, this 2D Riemann problem is currently not supported!")
 
    WRITE (teststr,'(A,I2)') "2D Riemann problem no. ", ic
 
    ! initial conditions
 
    ! Test configuration no. 1
    ! 1st quadrant
    den(1,1) = 1.
    pre(1,1) = 1.
    ! 2nd quadrant
    den(1,2) = .5197
    vx(1,2)  = -.7259
    pre(1,2) = .4
    ! 3rd quadrant
    den(1,3) = .1072
    vx(1,3)  = -.7259
    vy(1,3)  = -1.4045
    pre(1,3) = .0439
    ! 4th quadrant
    den(1,4) = .2579
    vy(1,4)  = -1.4045
    pre(1,4) = .15
    
    ! Test configuration no. 2
    ! 1st quadrant
    den(2,1) = 1.
    pre(2,1) = 1.
    ! 2nd quadrant
    den(2,2) = .5197
    vx(2,2) = -.7259
    pre(2,2) = .4
    ! 3rd quadrant
    den(2,3) = 1.
    vx(2,3) = -.7259
    vy(2,3) = -.7259
    pre(2,3) = 1.
    ! 4th quadrant
    den(2,4) = .5197
    vy(2,4) = -.7259
    pre(2,4) = .4
 
    ! Test configuration no. 3
    ! 1st quadrant
    den(3,1) = 1.5
    pre(3,1) = 1.5
    ! 2nd quadrant
    den(3,2) = .5323
    vx(3,2) = 1.206
    pre(3,2) = .3
    ! 3rd quadrant
    den(3,3) = 0.138
    vx(3,3) = 1.206
    vy(3,3) = 1.206
    pre(3,3) = 0.029
    ! 4th quadrant
    den(3,4) = .5323
    vy(3,4) = 1.206
    pre(3,4) = .3
 
    ! Test configuration no. 4
    ! 1st quadrant
    den(4,1) = 1.1
    pre(4,1) = 1.1
    ! 2nd quadrant
    den(4,2) = .5065
    vx(4,2) = .8939
    pre(4,2) = .35
    ! 3rd quadrant
    den(4,3) = 1.1
    vx(4,3) = .8939
    vy(4,3) = .8939
    pre(4,3) = 1.1
    ! 4th quadrant
    den(4,4) = .5065
    vy(4,4) = .8939
    pre(4,4) = .35
 
    ! 1st quadrant
    ! Test configuration no. 5
    den(5,1) = 1.
    vx(5,1) = -.75
    vy(5,1) = -.5
    pre(5,1) = 1.
    ! 2nd quadrant
    den(5,2) = 2.
    vx(5,2) = -.75
    vy(5,2) = .5
    pre(5,2) = 1.
    ! 3rd quadrant
    den(5,3) = 1.
    vx(5,3) = .75
    vy(5,3) = .5
    pre(5,3) = 1.
    ! 4th quadrant
    den(5,4) = 3.
    vx(5,4) = .75
    vy(5,4) = -.5
    pre(5,4) = 1.
 
    ! Test configuration no. 6
    ! 1st quadrant
    den(6,1) = 1.
    vx(6,1) = 0.75
    vy(6,1) = -0.5
    pre(6,1) = 1.
    ! 2nd quadrant
    den(6,2) = 2.
    vx(6,2) = 0.75
    vy(6,2) = 0.5
    pre(6,2) = 1.
    ! 3rd quadrant
    den(6,3) = 1.
    vx(6,3) = -0.75
    vy(6,3) = 0.5
    pre(6,3) = 1.
    ! 4th quadrant
    den(6,4) = 3.
    vx(6,4) = -0.75
    vy(6,4) = -0.5
    pre(6,4) = 1.
 
    ! Test configuration no. 7
    ! 1st quadrant
    den(7,1) = 1.
    vx(7,1) = 0.1
    vy(7,1) = 0.1
    pre(7,1) = 1.
    ! 2nd quadrant
    den(7,2) = 0.5197
    vx(7,2) = -0.6259
    vy(7,2) = 0.1
    pre(7,2) = 0.4
    ! 3rd quadrant
    den(7,3) = 0.8
    vx(7,3) = 0.1
    vy(7,3) = 0.1
    pre(7,3) = 0.4
    ! 4th quadrant
    den(7,4) = 0.5197
    vx(7,4) = 0.1
    vy(7,4) = -0.6259
    pre(7,4) = 0.4
 
    ! Test configuration no. 8
    ! 1st quadrant
    den(8,1) = 0.5197
    vx(8,1) = 0.1
    vy(8,1) = 0.1
    pre(8,1) = 0.4
    ! 2nd quadrant
    den(8,2) = 1.
    vx(8,2) = -0.6259
    vy(8,2) = 0.1
    pre(8,2) = 1.
    ! 3rd quadrant
    den(8,3) = 0.8
    vx(8,3) = 0.1
    vy(8,3) = 0.1
    pre(8,3) = 1.
    ! 4th quadrant
    den(8,4) = 1.
    vx(8,4) = 0.1
    vy(8,4) = -0.6259
    pre(8,4) = 1.
 
    ! Test configuration no. 9
    ! 1st quadrant
    den(9,1) = 1.
    vx(9,1) = 0.
    vy(9,1) = 0.3
    pre(9,1) = 1.
    ! 2nd quadrant
    den(9,2) = 2.
    vx(9,2) = 0.
    vy(9,2) = -0.3
    pre(9,2) = 1.
    ! 3rd quadrant
    den(9,3) = 1.039
    vx(9,3) = 0.
    vy(9,3) = -0.8133
    pre(9,3) = 0.4
    ! 4th quadrant
    den(9,4) = 0.5197
    vx(9,4) = 0.
    vy(9,4) = -0.4259
    pre(9,4) = 0.4
 
    ! Test configuration no. 10
    ! 1st quadrant
    den(10,1) = 1.
    vx(10,1) = 0.
    vy(10,1) = 0.4297
    pre(10,1) = 1.
    ! 2nd quadrant
    den(10,2) = 0.5
    vx(10,2) = 0.
    vy(10,2) = 0.6076
    pre(10,2) = 1.
    ! 3rd quadrant
    den(10,3) = 0.2281
    vx(10,3) = 0.
    vy(10,3) = -0.6076
    pre(10,3) = 0.3333
    ! 4th quadrant
    den(10,4) = 0.4562
    vx(10,4) = 0.
    vy(10,4) = -0.4297
    pre(10,4) = 0.3333
 
    ! Test configuration no. 11
    ! 1st quadrant
    den(11,1) = 1.
    vx(11,1) = 0.1
    vy(11,1) = 0.
    pre(11,1) = 1.
    ! 2nd quadrant
    den(11,2) = 0.5313
    vx(11,2) = 0.8276
    vy(11,2) = 0.
    pre(11,2) = 0.4
    ! 3rd quadrant
    den(11,3) = 0.8
    vx(11,3) = 0.1
    vy(11,3) = 0.
    pre(11,3) = 0.4
    ! 4th quadrant
    den(11,4) = 0.5313
    vx(11,4) = 0.1
    vy(11,4) = 0.7276
    pre(11,4) = 0.4
 
    ! Test configuration no. 12
    ! 1st quadrant
    den(12,1) = 0.5313
    vx(12,1) = 0.
    vy(12,1) = 0.
    pre(12,1) = 0.4
    ! 2nd quadrant
    den(12,2) = 1.
    vx(12,2) = 0.7276
    vy(12,2) = 0.
    pre(12,2) = 1.
    ! 3rd quadrant
    den(12,3) = 0.8
    vx(12,3) = 0.
    vy(12,3) = 0.
    pre(12,3) = 1.
    ! 4th quadrant
    den(12,4) = 1.
    vx(12,4) = 0.
    vy(12,4) = 0.7276
    pre(12,4) = 1.
 
    ! Test configuration no. 13
    ! 1st quadrant
    den(13,1) = 1.
    vx(13,1) = 0.
    vy(13,1) = -0.3
    pre(13,1) = 1.
    ! 2nd quadrant
    den(13,2) = 2.
    vx(13,2) = 0.
    vy(13,2) = 0.3
    pre(13,2) = 1.
    ! 3rd quadrant
    den(13,3) = 1.0625
    vx(13,3) = 0.
    vy(13,3) = 0.8145
    pre(13,3) = 0.4
    ! 4th quadrant
    den(13,4) = 0.5313
    vx(13,4) = 0.
    vy(13,4) = 0.4276
    pre(13,4) = 0.4
 
    ! Test configuration no. 14
    ! 1st quadrant
    den(14,1) = 2.
    vx(14,1) = 0.
    vy(14,1) = -0.5606
    pre(14,1) = 8.
    ! 2nd quadrant
    den(14,2) = 1.
    vx(14,2) = 0.
    vy(14,2) = -1.2172
    pre(14,2) = 8.
    ! 3rd quadrant
    den(14,3) = 0.4736
    vx(14,3) = 0.
    vy(14,3) = 1.2172
    pre(14,3) = 2.6667
    ! 4th quadrant
    den(14,4) = 0.9474
    vx(14,4) = 0.
    vy(14,4) = 1.1606
    pre(14,4) = 2.6667
 
    ! Test configuration no. 15
    ! 1st quadrant
    den(15,1) = 1.
    vx(15,1) = 0.1
    vy(15,1) = -0.3
    pre(15,1) = 1.
    ! 2nd quadrant
    den(15,2) = 0.5197
    vx(15,2) = -0.6259
    vy(15,2) = -0.3
    pre(15,2) = 0.4
    ! 3rd quadrant
    den(15,3) = 0.8
    vx(15,3) = 0.1
    vy(15,3) = -0.3
    pre(15,3) = 0.4
    ! 4th quadrant
    den(15,4) = 0.5313
    vx(15,4) = 0.1
    vy(15,4) = 0.4276
    pre(15,4) = 0.4
 
    ! Test configuration no. 16
    ! 1st quadrant
    den(16,1) = 0.5313
    vx(16,1) = 0.1
    vy(16,1) = 0.1
    pre(16,1) = 0.4
    ! 2nd quadrant
    den(16,2) = 1.02222
    vx(16,2) = -0.6179
    vy(16,2) = 0.1
    pre(16,2) = 1.
    ! 3rd quadrant
    den(16,3) = 0.8
    vx(16,3) = 0.1
    vy(16,3) = 0.1
    pre(16,3) = 1.
    ! 4th quadrant
    den(16,4) = 1.
    vx(16,4) = 0.1
    vy(16,4) = 0.8276
    pre(16,4) = 1.
 
    ! Test configuration no. 17
    ! 1st quadrant
    den(17,1) = 1.
    vx(17,1) = 0.
    vy(17,1) = -0.4
    pre(17,1) = 1.
    ! 2nd quadrant
    den(17,2) = 2.
    vx(17,2) = 0.
    vy(17,2) = -0.3
    pre(17,2) = 1.
    ! 3rd quadrant
    den(17,3) = 1.0625
    vx(17,3) = 0.
    vy(17,3) = 0.2145
    pre(17,3) = 0.4
    ! 4th quadrant
    den(17,4) = 0.5197
    vx(17,4) = 0.
    vy(17,4) = -1.1259
    pre(17,4) = 0.4
 
    ! Test configuration no. 18
    ! 1st quadrant
    den(18,1) = 1.
    vx(18,1) = 0.
    vy(18,1) = 1.
    pre(18,1) = 1.
    ! 2nd quadrant
    den(18,2) = 2.
    vx(18,2) = 0.
    vy(18,2) = -0.3
    pre(18,2) = 1.
    ! 3rd quadrant
    den(18,3) = 1.0625
    vx(18,3) = 0.
    vy(18,3) = 0.2145
    pre(18,3) = 0.4
    ! 4th quadrant
    den(18,4) = 0.5197
    vx(18,4) = 0.
    vy(18,4) = 0.2741
    pre(18,4) = 0.4
 
    ! Test configuration no. 19
    ! 1st quadrant
    den(19,1) = 1.
    vx(19,1) = 0.
    vy(19,1) = 0.3
    pre(19,1) = 1.
    ! 2nd quadrant
    den(19,2) = 2.
    vx(19,2) = 0.
    vy(19,2) = -0.3
    pre(19,2) = 1.
    ! 3rd quadrant
    den(19,3) = 1.0625
    vx(19,3) = 0.
    vy(19,3) = 0.2145
    pre(19,3) = 0.4
    ! 4th quadrant
    den(19,4) = 0.5197
    vx(19,4) = 0.
    vy(19,4) = -0.4259
    pre(19,4) = 0.4
 
    SELECT TYPE(pvar => timedisc%pvar)
    TYPE IS(statevector_euler) ! non-isothermal HD
      WHERE ( (mesh%bccart(:,:,:,1).GE.x0).AND.(mesh%bccart(:,:,:,2).GE.y0) )
        pvar%density%data3d(:,:,:)    = den(ic,1)
        pvar%pressure%data3d(:,:,:)   = pre(ic,1)
        vcart%data4d(:,:,:,1) = vx(ic,1)
        vcart%data4d(:,:,:,2) = vy(ic,1)
      ELSEWHERE ( (mesh%bccart(:,:,:,1).LT.x0).AND.(mesh%bccart(:,:,:,2).GE.y0) )
        pvar%density%data3d(:,:,:)    = den(ic,2)
        pvar%pressure%data3d(:,:,:)   = pre(ic,2)
        vcart%data4d(:,:,:,1) = vx(ic,2)
        vcart%data4d(:,:,:,2) = vy(ic,2)
      ELSEWHERE ( (mesh%bccart(:,:,:,1).LT.x0).AND.(mesh%bccart(:,:,:,2).LT.y0) )
        pvar%density%data3d(:,:,:)    = den(ic,3)
        pvar%pressure%data3d(:,:,:)   = pre(ic,3)
        vcart%data4d(:,:,:,1) = vx(ic,3)
        vcart%data4d(:,:,:,2) = vy(ic,3)
      ELSEWHERE ( (mesh%bccart(:,:,:,1).GE.x0).AND.(mesh%bccart(:,:,:,2).LT.y0) )
        pvar%density%data3d(:,:,:)    = den(ic,4)
        pvar%pressure%data3d(:,:,:)   = pre(ic,4)
        vcart%data4d(:,:,:,1) = vx(ic,4)
        vcart%data4d(:,:,:,2) = vy(ic,4)
      END WHERE
      vcart%data2d(:,3) = 0.0 ! z-velocity
      CALL mesh%geometry%Convert2Curvilinear(mesh%bcenter,vcart%data4d,vcurv%data4d)
      pvar%velocity%data2d(:,1:2) = vcurv%data2d(:,1:2)
    END SELECT
   
    CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
    CALL mesh%Info(" DATA-----> initial condition: " // trim(teststr))
 
  END SUBROUTINE initdata
 
END PROGRAM riemann2d