c
c   2002/08/12/ : start to make CIP-CSL2 2D(Directional splitting)
c
c    df/dt + d(uf)/dx + d(vf)/dy = 0
c	Doswell Front Genesis 
c -----------------------------------------------------------------
        INCLUDE 'PARAMETER'
        COMMON /VALUE/ F(0:IM,0:JM),SIGX(0:IM,0:JM),SIGY(0:IM,0:JM)
     1			,RHO(0:IM,0:JM)
        COMMON /VELO/  YU(0:IM,0:JM),YV(0:IM,0:JM),VELX(0:IM,0:JM)
     1			,VELY(0:IM,0:JM)
        COMMON /PARAM/ DX,DY,DT,IEND,ISTEP,TIME
	COMMON ITIME,IOUTPUT,EXM
	DIMENSION X(0:IM),Y(0:JM)

c---------------------------------------------------------------------
c  file open
c                                                     netcdf - start >>>
      ix=im+1
      jx=jm+1
      call ncopn2d(idf,'out2d.cdf',idn,idi,ix,idj,jx)
      call ncputattc(idf,'comment','CSL2-2D-FrontGenesis')
      call ncdefss(idf,'time',idt,idn)
      call ncdefs2(idf,'F',idval,idn,idi,idj)
      call ncdefs2(idf,'RHO',idrho,idn,idi,idj)
      call ncdefs2(idf,'Vx',idyu,idn,idi,idj)
      call ncdefs2(idf,'Vy',idyv,idn,idi,idj)
      nd=1

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


c       << SET INITIAL VALUE + GRAPH >>
         CALL INITIAL
            IEND = int(16.2d0/DT)

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(idf,idi,'x',x,ix)
      call ncputo1(idf,idj,'y',y,jx)

c                                                     netcdf - start >>>
      istep=0
      time=istep*dt
      call ncputss(idf,idt,nd,time)
      call ncputs2(idf,idval,nd,F,ix,jx)
      call ncputs2(idf,idrho,nd,RHO,ix,jx)
      call ncputs2(idf,idyu,nd,VELX,ix,jx)
      call ncputs2(idf,idyv,nd,VELY,ix,jx)
            mstatus=nf_sync(idf)
      nd=nd+1
c                                                     netcdf - end   <<<
c                                                     netcdf - end   <<<
c---------------------------------------------------------------------



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

	  CALL VELOCITY
        DO 1000 ISTEP=1,IEND

	  CALL OUTPUT
          CALL BOUNDARY

c    ------ X-direction ------------
          CALL CIPCSL1DX(F,SIGX,1)
          CALL CIPCSL1DX(SIGY,RHO,2)
	  CALL NONADVX

        CALL BOUNDARY

c    ------ Y-direction ------------
          CALL CIPCSL1DY(F,SIGY,1)
          CALL CIPCSL1DY(SIGX,RHO,2)
	  CALL NONADVY

        CALL BOUNDARY

            time=time+dt
            write(6,*) 'TIME= ',time
	  CALL ERROR

c---------------------------------------------------------------------
c                                                     netcdf - start >
	if(mod(ISTEP,80).eq.0) then
      call ncputss(idf,idt,nd,time)
      call ncputs2(idf,idval,nd,F,ix,jx)
      call ncputs2(idf,idrho,nd,RHO,ix,jx)
      call ncputs2(idf,idyu,nd,VELX,ix,jx)
      call ncputs2(idf,idyv,nd,VELY,ix,jx)
            mstatus=nf_sync(idf)
      nd=nd+1
	 write(6,*) "netCDF output"
	end if
c                                                     netcdf - end   <<<
c---------------------------------------------------------------------

 1000   CONTINUE

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

c                                                     netcdf - end   <<<
                               
        STOP
        END

c  <<<<<<<<< END >>>>>>>>>>>>>>>>>>>>

C *******************************************************
        SUBROUTINE INITIAL
C  << INITIAL VALUE SETTING >
C
C *****************************************************
        INCLUDE "PARAMETER"
        INCLUDE "CSL2.ini"

        RETURN
        END


C *******************************************************
        SUBROUTINE VELOCITY
c   <<< X,Y Velocity on F >>>
C
C *****************************************************
        INCLUDE "PARAMETER"
        INCLUDE "CSL2.vel"

        RETURN
        END


C *******************************************************
        SUBROUTINE CIPCSL1DX(YR,SX,N)
C  << ADVECTION for Value : X-direction >>
C
C *****************************************************
        INCLUDE "PARAMETER"
        INCLUDE "CSL2.advx"

        RETURN
        END

C *******************************************************
        SUBROUTINE CIPCSL1DY(YR,SY,N)
C  << ADVECTION for Value : Y-direction >>
C
C *****************************************************
        INCLUDE "PARAMETER"
        INCLUDE "CSL2.advy"

        RETURN
        END

C *******************************************************
        SUBROUTINE NONADVX
C  << Non-ADVECTION for Value : X-direction >>
C
C *****************************************************
        INCLUDE "PARAMETER"
        INCLUDE "CSL2.nonx"

        RETURN
        END

C *******************************************************
        SUBROUTINE NONADVY
C  << Non-ADVECTION for Value : Y-direction >>
C
C *****************************************************
        INCLUDE "PARAMETER"
        INCLUDE "CSL2.nony"

        RETURN
        END

C *******************************************************
        SUBROUTINE BOUNDARY
C
C *****************************************************
        INCLUDE "PARAMETER"
        INCLUDE "CSL2.bou"

        RETURN
        END

C *******************************************************
        SUBROUTINE OUTPUT
C
C *****************************************************
        INCLUDE "PARAMETER"
        INCLUDE "CSL2.out"

        RETURN
        END


C *******************************************************
        SUBROUTINE ERROR
C
C *****************************************************
        INCLUDE "PARAMETER"
        INCLUDE "CSL2.err"

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









