C *****************************************************
c
c  1D HYDRODYNAMIC EQATIONS(SHOCK TUBE PROBREM)
c
c        WITH CIP-CSL2 METHOD
c    2002-08-03 Yoichi Ogata made
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','1D-shocktube: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,'rho',idrho,idn,idi)
      call ncdefs1(idf,'gyp',idgyp,idn,idi)
      call ncdefs1(idf,'gyu',idgyu,idn,idi)
      nd=1

c                                                     netcdf - end   <<<
c---------------------------------------------------------------------


c       << SET INITIAL VALUE + GRAPH >>
        CALL INITIAL

            NT = 10*(0.4/DT)


c---------------------------------------------------------------------
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,idrho,nd,rho,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---------------------------------------------------------------------




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
            CALL OUTPUT(ISTEP,NT)


c---------------------------------------------------------------------
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,idrho,nd,rho,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---------------------------------------------------------------------

 1000   CONTINUE

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

c                                                     netcdf - end   <<<


        STOP
        END

C *******************************************************
        SUBROUTINE INITIAL
C  << INITIAL VALUE SETTING >
C
C *****************************************************

        INCLUDE 'PARAMETER'
        INCLUDE "CSL2.ini"

	RETURN
	END

C ******************************************************
        SUBROUTINE ADVECTION(YY,GG,UU)
C   <<<  ADVECTION PHASE for  pressure, velocity : normal CIP>>
C
C *******************************************************
        INCLUDE 'PARAMETER'
        INCLUDE "CIP.adv"

        RETURN
        END

C *******************************************************
        SUBROUTINE ADVECTIONCSL
C  << ADVECTION for density : CIP-CSL2>>
C
C *****************************************************
        INCLUDE 'PARAMETER'
        INCLUDE "CSL2.adv"

        RETURN
        END

C ***********************************************************
        SUBROUTINE VISCOCITY

C       <<  ARTIFICIAL VISCOCITY >>
C
C **********************************************************
        INCLUDE 'PARAMETER'
        INCLUDE "CSL2.vis"

        RETURN
        END

C ******************************************************
        SUBROUTINE VELOCITY
C        << ADVECTION VELOCITY FOR DENSITY AND ENERGY >>
C
C ******************************************************
        INCLUDE 'PARAMETER'

        COMMON /VALUE/ YR(0:IM),YU(0:IM),YP(0:IM)
        COMMON /PARAM/ DX,DT,GAMMA
        COMMON /ADVECV/ UUU(0:IM)

        DO 2000 I=1,NX

           UUU(I)=(YU(I)+YU(I-1))*0.5

 2000   CONTINUE

        RETURN
        END

C *****************************************************
	SUBROUTINE PRESSURE
C         
C ****************************************************
        INCLUDE 'PARAMETER'
        INCLUDE "CSL2.pre"

	RETURN
	END 


C *****************************************************
	SUBROUTINE VELO1
C         
C ****************************************************
        INCLUDE 'PARAMETER'
        INCLUDE "CSL2.vel1"

	RETURN
	END 

C *****************************************************
	SUBROUTINE VELO2
C         
C ****************************************************
        INCLUDE 'PARAMETER'
        INCLUDE "CSL2.vel2"

	RETURN
	END 

C *****************************************************
	SUBROUTINE DENSITY 
C         
C ****************************************************
        INCLUDE 'PARAMETER'
        INCLUDE "CSL2.del"

	RETURN
	END 


C *****************************************************
	SUBROUTINE SHIFT
C         
C ****************************************************
        INCLUDE 'PARAMETER'
        INCLUDE "CSL2.shf"

	RETURN
	END 


c **************************************************************
        SUBROUTINE BOUNDARY
C   <<  BOUNDARY CONDITION >>
C
C **************************************************************
        INCLUDE 'PARAMETER'
        INCLUDE "CSL2.bou"

        RETURN
        END

C **********************************************************
        SUBROUTINE NONGRAD(GG,YY,YYN,YU)

C     << NONADVECTION OF GRADIENT >>
C
C **********************************************************
        INCLUDE 'PARAMETER'
        INCLUDE "CSL2.grad"

        RETURN
        END

C ***************************************************
	SUBROUTINE OUTPUT(ISTEP,NT)
C          <<< for gnuplot >>>
C *************************************************
        INCLUDE 'PARAMETER'
        INCLUDE "CSL2.out"

	return
	end

c -------------- A series of NetCDF data routine : 
c	Do not change,,,
c -------------------------------------------------

c======================================================================|
      subroutine ncopn1d(idf,file,idn,idi,ix)
c======================================================================|
c----------------------------------------------------------------------|
      INCLUDE 'PARAMETER'
      character*(*) file
c----------------------------------------------------------------------|

      mstatus = nf_create(file,0,idf)

      mstatus = nf_def_dim(idf,'ix',ix,idi)
      mstatus = nf_def_dim(idf,'nx', 0,idn)

      mstatus = nf_enddef(idf)

      return
      end

c======================================================================|
      subroutine ncputattc(idf,name,data_char)
c======================================================================|
c----------------------------------------------------------------------|
      INCLUDE 'PARAMETER'
      character*(*) name
      character*(*) data_char
c----------------------------------------------------------------------|

      mstatus = nf_redef(idf)
      mstatus = nf_put_att_text(idf,0,name,len(data_char),data_char)
      mstatus = nf_enddef(idf)

      return
      end

c======================================================================|
      subroutine ncdefss(idf,name,idd,idn)
c======================================================================|
c----------------------------------------------------------------------|
      INCLUDE 'PARAMETER'
      character*(*) name
c----------------------------------------------------------------------|

      mstatus = nf_redef(idf)
      mstatus = nf_def_var(idf,name,5,1,idn,idd)
      mstatus = nf_enddef(idf)

      return
      end

c======================================================================|
      subroutine ncdefs1(idf,name,idd,idn,idi)
c======================================================================|
c----------------------------------------------------------------------|
      INCLUDE 'PARAMETER'
      parameter(mrank=2)
      integer mdims(mrank)
      character*(*) name
c----------------------------------------------------------------------|

      mdims(1)=idi
      mdims(2)=idn

      mstatus = nf_redef(idf)
      mstatus = nf_def_var(idf,name,5,mrank,mdims,idd)
      mstatus = nf_enddef(idf)

      return
      end

c======================================================================|
      subroutine ncputo1(idf,idi,name,data_dbl,ix)
c======================================================================|
c----------------------------------------------------------------------|
      INCLUDE 'PARAMETER'
      character*(*) name
      dimension data_dbl(ix)
      real data(ix)
c----------------------------------------------------------------------|

      do i=1,ix
        data(i)=data_dbl(i)
      enddo

      mstatus = nf_redef(idf)
      mstatus = nf_def_var(idf,name,5,1,idi,idd)
      mstatus = nf_enddef(idf)

      mstatus=nf_put_var_real(idf,idd,data)

      return
      end


c======================================================================|
      subroutine ncputos(idf,name,data_dbl)
c======================================================================|
c----------------------------------------------------------------------|
      INCLUDE 'PARAMETER'
      character*(*) name
c----------------------------------------------------------------------|

      mstatus = nf_redef(idf)
      mstatus = nf_def_var(idf,name,5,0,0,idd)
      mstatus = nf_enddef(idf)

      mstatus=nf_put_var_real(idf,idd,real(data_dbl))

      return
      end

c======================================================================|
      subroutine ncputss(idf,idd,nd,data_dbl)
c======================================================================|
c----------------------------------------------------------------------|
      INCLUDE 'PARAMETER'
      real data
c----------------------------------------------------------------------|

      data=data_dbl
      mstatus=nf_put_var1_real(idf,idd,nd,data)

      return
      end

c======================================================================|
      subroutine ncputs1(idf,idd,nd,data_dbl,ix)
c======================================================================|
c----------------------------------------------------------------------|
      INCLUDE 'PARAMETER'
      parameter(mrank=2)
      dimension data_dbl(ix)
      integer mstart(mrank),mcount(mrank)
      real data(ix)
c----------------------------------------------------------------------|

      do i=1,ix
         data(i)=data_dbl(i)
      enddo

      mcount(1)=ix
      mcount(2)=1

      mstart(1)=1
      mstart(2)=nd

      mstatus=nf_put_vara_real(idf,idd,mstart,mcount,data)

      return
      end



