C *****************************************************
c
c  1D HYDRODYNAMIC EQATIONS(SHOCK TUBE PROBREM)
c
c        WITH CIP METHOD
c         << ADVECTION TO NON-ADAVECTION >>
c *****************************************************

        INCLUDE 'PARAMETER'

        COMMON /VALUE/ YR(0:IM),YU(0:IM),YP(0:IM)
        COMMON /NEXT/  YRN(0:IM),YUN(0:IM),YPN(0:IM)
        COMMON /GRAD/  GYU(0:IM),GYP(0:IM)
	COMMON /INTEG/ RHO(0:IM),RHON(0:IM)
        COMMON /VISCO/ VIS(0:IM)
        COMMON /PARAM/ DX,DT,GAMMA
        COMMON /BOUND/ BOUND
        COMMON /ADVECV/ UUU(0:IM)
        COMMON EXM,CALM,EXF,CALF

        dimension x(0:im)


c----------------------------------------------------------------------|
c  file open
c                                                     netcdf - start >>>
      ix=im+1
      call ncopn1d(idf,'out.cdf',idn,idi,ix)
      call ncputattc(idf,'comment','CSL2')
      call ncdefss(idf,'t',idt,idn)
      call ncdefs1(idf,'yr',idyr,idn,idi)
      call ncdefs1(idf,'yp',idyp,idn,idi)
      call ncdefs1(idf,'yu',idyu,idn,idi)
      call ncdefs1(idf,'gyp',idgyp,idn,idi)
      call ncdefs1(idf,'gyu',idgyu,idn,idi)
      nd=1

c                                                     netcdf - end   <<<

c       << SET INITIAL VALUE + GRAPH >>
        CALL INITIAL
c        CALL initg1d(IM)

            NT = 10*(0.4/DT)

c                                                     netcdf - start >>>
      do i=0,im
        x(i)=dx*real(i)
      enddo

      call ncputo1(idf,idi,'x',x,ix)
      call ncputos(idf,'gamma',gamma)
c                                                     netcdf - start >>>
      istep=0
      t=istep*dt
      call ncputss(idf,idt,nd,t)
      call ncputs1(idf,idyr,nd,yr,ix)
      call ncputs1(idf,idyp,nd,yp,ix)
      call ncputs1(idf,idyu,nd,yu,ix)
      call ncputs1(idf,idgyp,nd,gyp,ix)
      call ncputs1(idf,idgyu,nd,gyu,ix)
            mstatus=nf_sync(idf)
      nd=nd+1
c                                                     netcdf - end   <<<
c                                                     netcdf - end   <<<

c  <<<< STRAT CALCURATION >>>>>>

        DO 1000 ISTEP=1,NT+2
        write(6,*) "ISTEP=",ISTEP


C        << ARTIFICIAL VISCOCITY >>
            CALL VISCOCITY

C        << ADVECTION PHASE >>

            CALL VELOCITY

            CALL ADVECTIONCSL
            CALL ADVECTION(YP,GYP,UUU)
            CALL ADVECTION(YU,GYU,YU)


C        << NONADVECTION PHASE >>

            CALL VELOCITY
	    CALL BOUNDARY

	    CALL VELO1	

	    CALL PRESSURE
	    CALL NONGRAD(GYP,YP,YPN,YU)
	    CALL DENSITY
	    CALL VELO2	
	    CALL NONGRAD(GYU,YU,YUN,YU)

c        << boundary condition >>
	    CALL BOUNDARY

c         << SHIFT I STEP TO (I+1)STEP >>
	    CALL SHIFT
	
c	    call graph1dd(YR)

C        <<  ANALYTIC SOLUTION OF SHOCK TUBE >>
            CALL ANASHOCK(ISTEP,NT)
 	    CALL OUTPUT(ISTEP,NT)

c                                                     netcdf - start >>>
      t=istep*dt
      call ncputss(idf,idt,nd,t)
      call ncputs1(idf,idyr,nd,yr,ix)
      call ncputs1(idf,idyp,nd,yp,ix)
      call ncputs1(idf,idyu,nd,yu,ix)
      call ncputs1(idf,idgyp,nd,gyp,ix)
      call ncputs1(idf,idgyu,nd,gyu,ix)
            mstatus=nf_sync(idf)
      nd=nd+1
c                                                     netcdf - end   <<<


 1000   CONTINUE

c       call get

c----------------------------------------------------------------------|
c  file close
c                                                     netcdf - start >>>
      mstatus=nf_close(idf)
c                                                     netcdf - end   <<<


        STOP
        END
