c ****************************************************
c    2D HYDRODYNAMICS(TYPE A)
c       <<<< Reflection of shock  >>>>
c *****************************************************

	INCLUDE 'PARAMETER'
	COMMON /VALUE/ YR(0:IM,0:JM),YU(0:IM,0:JM),YV(0:IM,0:JM)
     1                ,YP(0:IM,0:JM)
	COMMON /NEXTVAL/ YRN(0:IM,0:JM),YUN(0:IM,0:JM),
     1                   YVN(0:IM,0:JM),YPN(0:IM,0:JM)
	COMMON /VELO/ UIJ(0:IM,0:JM),VIJ(0:IM,0:JM),
     1                UU(0:IM,0:JM),VV(0:IM,0:JM)
	COMMON /GRADX/ GRX(0:IM,0:JM),GUX(0:IM,0:JM),GVX(0:IM,0:JM)
     1               ,GPX(0:IM,0:JM)
	COMMON /GRADY/ GRY(0:IM,0:JM),GUY(0:IM,0:JM),GVY(0:IM,0:JM)
     1               ,GPY(0:IM,0:JM)
	COMMON /PARAM/ DX,DY,GAMMA,DT
        COMMON /GR/ DEN(0:IM),DENTAI(0:IM)
	COMMON /WALL/ IW,IS

	COMMON /SHOCK/ OMA,THETA,DEN1,YU1,YV1,PRE1,
     1                  DEN2,YU2,YV2,PRE2

	COMMON /VISC/ XVISC(0:IM,0:JM),YVISC(0:IM,0:JM),
     1                CVX,CVY
        COMMON /IOUT/ ITIME,IOUTPUT,ISTEP,IEND,TIME
        DIMENSION X(0:IM),Y(0:JM)

c---------------------------------------------------------------------
c  file open : 2-dimension
c                                                     netcdf - start >>>
      ix=im+1
      jx=jm+1
      call ncopn2d(idf2d,'out2d.cdf',idn,idi1,ix,idj1,jx)
      call ncputattc(idf2d,'comment','CCUP-data2D')
      call ncdefss(idf2d,'time',idt,idn)
      call ncdefs2(idf2d,'Vx',idvx,idn,idi1,idj1)
      call ncdefs2(idf2d,'Vy',idvy,idn,idi1,idj1)
      call ncdefs2(idf2d,'den',idden,idn,idi1,idj1)
      call ncdefs2(idf2d,'pre',idpre,idn,idi1,idj1)
      nd=1

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


	CALL INISHOCK
	CALL INITIAL

c---------------------------------------------------------------------
c                                                     netcdf - start >>>
c    ---- coordinate -----
      do i=0,im
        x(i)=dx*real(i)
      enddo
      do j=0,jm
        y(j)=dy*real(j)
      enddo
            TIME = 0.0
c    ----------------------
      call ncputo1(idf2d,idi1,'x',x,ix)
      call ncputo1(idf2d,idj1,'y',y,jx)
c                                                     netcdf - start >>>
      istep=0
      time=istep*dt
      call ncputss(idf2d,idt,nd,time)
      call ncputs2(idf2d,idden,nd,yr,ix,jx)
      call ncputs2(idf2d,idpre,nd,yp,ix,jx)
      call ncputs2(idf2d,idvx,nd,yu,ix,jx)
      call ncputs2(idf2d,idvy,nd,yv,ix,jx)
            mstatus1=nf_sync(idf2d)

      nd=nd+1
         write(6,*) "netCDF output"
c                                                     netcdf - end   <<<
c---------------------------------------------------------------------

	  write(6,*) "nx,ny=",IM,JM 

C     <<<<<<< CALCULATION >>>>>>>>
	DO 1000 ISTEP=1,IEND+10

	  write(6,*) "    TIME=",DT*(ISTEP-1.0)
	  write(6,*) "cv=",CVX,"  MACH number=",OMA
	  write(6,*) "nx,ny=",IM,JM 

	
	    CALL INFLOW
	    CALL BOUNDDENSE(YR,GRX,GRY)
	    CALL BOUNDPRESS(YP,GPX,GPY)
	    CALL BOUNDU(YU,GUX,GUY)
	    CALL BOUNDV(YV,GVX,GVY)
              CALL VISCOCITY
              CALL VELOCITY
	    CALL BOUNDUU(UU,YU)
	    CALL BOUNDVV(VV,YV)
	    CALL BOUNDUIJVIJ(UIJ,VIJ)


C          << ADVECTION >>

	     CALL ADVECTION(YR,GRX,GRY,UIJ,VIJ)
	     CALL ADVECTION(YP,GPX,GPY,UIJ,VIJ)
	     CALL ADVECTION(YU,GUX,GUY,YU,VV)
	     CALL ADVECTION(YV,GVX,GVY,UU,YV)


	
	    CALL INFLOW
	    CALL BOUNDDENSE(YR,GRX,GRY)
	    CALL BOUNDPRESS(YP,GPX,GPY)
	    CALL BOUNDU(YU,GUX,GUY)
	    CALL BOUNDV(YV,GVX,GVY)
            CALL VELOCITY
	    CALL BOUNDUU(UU,YU)
	    CALL BOUNDVV(VV,YV)
	    CALL BOUNDUIJVIJ(UIJ,VIJ)


C          << NONADVECTION >>	  

	    CALL NONDEN
	    CALL NONVELX
	    CALL NONVELY
	    CALL NONPRESS

	    CALL GRADXX(YR,YRN,GRX,GRY,YU,VV)
	    CALL GRADXX(YP,YPN,GPX,GPY,YU,VV)
	    CALL GRADYY(YR,YRN,GRX,GRY,UU,YV)
	    CALL GRADYY(YP,YPN,GPX,GPY,UU,YV)
	    CALL GRAVELX(YU,YUN,GUX,GUY,VIJ,VV)
	    CALL GRAVELY(YV,YVN,GVX,GVY,UIJ,UU)


	    CALL BOUNDDENSE(YRN,GRX,GRY)
	    CALL BOUNDPRESS(YPN,GPX,GPY)
	    CALL BOUNDU(YUN,GUX,GUY)
	    CALL BOUNDV(YVN,GVX,GVY)


C------- OUTPUT
	   	CALL OUTPUT

C------- REPLACE N+1 -> N
            CALL SHIFT

c---------------------------------------------------------------------
c                                                     netcdf - start >>>
        if(mod(ISTEP,IOUTPUT).eq.0) then
      time=istep*dt
      call ncputss(idf2d,idt,nd,time)
      call ncputs2(idf2d,idden,nd,yr,ix,jx)
      call ncputs2(idf2d,idpre,nd,yp,ix,jx)
      call ncputs2(idf2d,idvx,nd,yu,ix,jx)
      call ncputs2(idf2d,idvy,nd,yv,ix,jx)
            mstatus1=nf_sync(idf2d)
      nd=nd+1
         write(6,*) "netCDF output"
        end if
c                                                     netcdf - end   <<<
c---------------------------------------------------------------------

            IF(ISTEP.EQ.IEND) THEN
                write(6,*) "calc end"
            END IF

 1000    CONTINUE

c----------------------------------------------------------------------|
c  file close
c                                                     netcdf - start >>>
      mstatus1=nf_close(idf2d)

c                                                     netcdf - end   <<<


        STOP
        END
C **************************************************************
	SUBROUTINE INISHOCK
C
C **************************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.inish'

	RETURN
	END


C **************************************************************
	SUBROUTINE INITIAL
C
C **************************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.ini'

	RETURN
	END

C *****************************************************
	SUBROUTINE VISCOCITY
C
C *****************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.vis'

	RETURN
	END

C **************************************************
	SUBROUTINE VELOCITY
C
C **************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.vel'

	RETURN
	END

C *****************************************************************
	SUBROUTINE NONDEN
C
C ****************************************************************

	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.den'

	RETURN
	END
C ***********************************************************
	SUBROUTINE NONVELX
C
C ***********************************************************

	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.velx'

 	RETURN	
	END

C *************************************************************
	SUBROUTINE NONVELY
C
C ***********************************************************

	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.vely'

 	RETURN	
	END

C ******************************************************
	SUBROUTINE NONPRESS
C
C ******************************************************0
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.pre'

	RETURN
	END

C **********************************************************
	SUBROUTINE GRADXX(Y,YN,GX,GY,U,V)
C
C *********************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.graxx'

	RETURN
	END

C **********************************************************
	SUBROUTINE GRADYY(Y,YN,GX,GY,U,V)
C
C *********************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.grayy'

	RETURN
	END
	
C **********************************************************
        SUBROUTINE GRAVELX(YU,YUN,GUX,GUY,VIJ,VV)
C
C **********************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.gravx'

	RETURN
	END

C **********************************************************
        SUBROUTINE GRAVELY(YV,YVN,GVX,GVY,UIJ,UU)
C
C **********************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.gravy'

	RETURN
	END

C ****************************************************************
	SUBROUTINE ADVECTION(F,GX,GY,U,V)
C
C ****************************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.adv'

	RETURN
	END

C ********************************************************
	SUBROUTINE SHIFT
C
C *********************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.shf'

	RETURN
	END

C ************************************************************
        SUBROUTINE OUTPUT
c         Gnuplot 
C
C ************************************************************
        INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.out'

        RETURN
        END
	
C ************************************************************
	SUBROUTINE INFLOW
C
C         <<<< INFLOW BOUNDARY >>>>
C ***********************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.inflow'

	RETURN
	END

C ************************************************************
	SUBROUTINE BOUNDDENSE(Y,GX,GY)
C
C         <<<< DENSITY,PRESSUER BOUNDARY >>>>
C ***********************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.bouden'

	RETURN
	END

c ************************************************************
        SUBROUTINE BOUNDPRESS(Y,GX,GY)
C
C         <<<< DENSITY,PRESSUER BOUNDARY >>>>
C ***********************************************************
        INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.boupre'

        RETURN
        END

C ************************************************************
	SUBROUTINE BOUNDU(Y,GX,GY)
C
C         <<<< X-VELOCITY BOUNDARY >>>>
C ***********************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.bouu'

	RETURN
	END

C ************************************************************
	SUBROUTINE BOUNDV(Y,GX,GY)
C
C         <<<< Y-VELOCITY BOUNDARY >>>>
C ***********************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.bouv'

	RETURN
	END

C ************************************************************
	SUBROUTINE BOUNDUU(UU,YU)
C
C         <<<< Y-VELOCITY BOUNDARY >>>>
C ***********************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.bouuu'

	RETURN
	END

C ************************************************************
	SUBROUTINE BOUNDVV(VV,YV)
C
C         <<<< Y-VELOCITY BOUNDARY >>>>
C ***********************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.bouvv'

	RETURN
	END

C ************************************************************
	SUBROUTINE BOUNDUIJVIJ(UIJ,VIJ)
C
C         <<<< Y-VELOCITY BOUNDARY >>>>
C ***********************************************************
	INCLUDE 'PARAMETER'
	INCLUDE 'CIP2D.bouS'

	RETURN
	END


c ------------------ A series of nedCDF file routine ----------
c======================================================================|
      subroutine ncopn2d(idf,file,idn,idi,ix,idj,jx)
c======================================================================|
c----------------------------------------------------------------------|
        INCLUDE "PARAMETER"
      parameter(mrank=3)
      character*(*) file
c----------------------------------------------------------------------|
      mstatus = nf_create(file,0,idf)

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

      mstatus = nf_enddef(idf)

      return
      end

c======================================================================|
      subroutine ncdefs2(idf,name,idd,idn,idi,idj)
c======================================================================|
c----------------------------------------------------------------------|
        INCLUDE "PARAMETER"
      parameter(mrank=3)
      integer mdims(mrank)
      character*(*) name
c----------------------------------------------------------------------|

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

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

      return
      end

c======================================================================|
      subroutine ncputs2(idf,idd,nd,data_dbl,ix,jx)
c======================================================================|
c----------------------------------------------------------------------|
        INCLUDE "PARAMETER"
      parameter(mrank=3)
      dimension data_dbl(ix,jx)
      integer mstart(mrank),mcount(mrank)
      real data(ix,jx)
c----------------------------------------------------------------------|

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

      mcount(1)=ix
      mcount(2)=jx
      mcount(3)=1

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

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

      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 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 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












