Program Main

  character(*), parameter:: InputFile='input.txt',OutputFile='data.plt' ! names of input and output files
  character MeshFile*30        ! name of file with computational mesh
  integer, parameter:: IO = 12 ! input-output unit
  real,allocatable,dimension(:,:):: X,Y,P,CellVolume ! scalar arrays
  real,allocatable,dimension(:,:,:):: CellCenter,IFaceCenter,IFaceVector,JFaceCenter,JFaceVector ! vector arrays
  real,allocatable,dimension(:,:,:):: GradP, GradPExact, GradPError, V
  real,allocatable,dimension(:,:):: DivV, DivVExact, DivVError
  real,allocatable,dimension(:,:):: LapP, LapPExact, LapPError

  INTEGER Niter, it, mode

!===  READ INPUT FILE ===
  WRITE(*,*) 'Read input file: ', InputFile
  OPEN(IO,FILE=InputFile)
  READ(IO,*) MeshFile  ! read name of file with computational mesh
  READ(IO,*) Niter
  READ(IO,*) mode
  CLOSE(IO)

!===   READ NODES NUMBER (NI,NJ) FROM FILE WITH MESH ===
  WRITE(*,*) 'Read nodes number from file: ', MeshFile
  OPEN(IO,FILE = MeshFile)
  READ(IO,*) NI,NJ
  WRITE(*,*) 'NI, NJ = ',NI,NJ

!=== ALLOCATE ALL ARRAYS ===
  WRITE(*,*) 'Allocate arrays'
  allocate(X(NI,NJ)) ! mesh nodes X-coordinates
  allocate(Y(NI,NJ)) ! mesh nodes Y-coordinates
  allocate(P(0:NI,0:NJ))   ! Pressure
  allocate(CellVolume(NI-1,NJ-1))   ! Cell Volumes
  allocate(CellCenter(0:NI,0:NJ,2)) ! Cell Centers
  allocate(IFaceCenter( NI,NJ-1,2)) ! Face Centers for I-faces
  allocate(IFaceVector( NI,NJ-1,2)) ! Face Vectors for I-faces
  allocate(JFaceCenter( NI-1,NJ,2)) ! Face Centers for J-faces
  allocate(JFaceVector( NI-1,NJ,2)) ! Face Vectors for I-faces
  allocate(GradP(0:NI,0:NJ,2))  !Gradient
  allocate(GradPExact(0:NI,0:NJ,2))  !GradientP Exact
  allocate(GradPError(0:NI,0:NJ,2))  !GradientP Error

  allocate(V(0:NI,0:NJ,2))  !V
  allocate(DivV(0:NI,0:NJ))  !DivV
  allocate(DivVExact(0:NI,0:NJ))  !DivVExact
  allocate(DivVError(0:NI,0:NJ))  !DivVError

  allocate(LapP(0:NI,0:NJ))  !LapP
  allocate(LapPExact(0:NI,0:NJ))  !LapPExact
  allocate(LapPError(0:NI,0:NJ))  !LapPError


!===  READ GRID ===
  WRITE(*,*) 'Read mesh from file: ', MeshFile
  READ(IO,*) ((X(I,J),Y(I,J),I=1,NI),J=1,NJ)
  CLOSE(IO)

!=== CALCULATE METRIC ===
  WRITE(*,*) 'Calculate metric'
  Call B_CalcMetric(NI,NJ,X,Y,CellCenter,CellVolume,IFaceCenter,IFaceVector,JFaceCenter,JFaceVector)

!=== INITIATE FIELDS ===
  WRITE(*,*) 'Initiate fields'
  DO  J = 0,NJ
    DO  I = 0,NI
      P(I,J) = Pressure(CellCenter(I,J,1),CellCenter(i,j,2))
      call CalcGradPExact(CellCenter(I,J,1),CellCenter(i,j,2),GradPExact(I,J,:))
      call Velocity(CellCenter(I,J,1),CellCenter(I,J,2),V(I,J,:))
      IF (mode.eq.0) THEN
        DivVExact(I,J)=DivVelocityExact(CellCenter(I,J,1),CellCenter(I,J,2))
        ELSE IF (mode.ge.1) THEN
        DivVExact(I,J)=DivVelocityPExact(CellCenter(I,J,1),CellCenter(I,J,2))
    END IF
    LapPExact(I,J)=RLaplacianP(CellCenter(I,J,1),CellCenter(I,J,2))
    ENDDO
  ENDDO

!=== CALCULATE GRADIENT ===
  WRITE(*,*) 'Calculate derivatives'

  GradP=0.0
  Do it =1, Niter

  Call B_CalcGradient(NI,NJ,P,GradP,CellCenter,CellVolume,IFaceCenter,IFaceVector,JFaceCenter,JFaceVector)
  GradPError=ABS(GradPExact - GradP)/GradPExact
  write(*,*) 'Iteration ', it, 'maximum error', maxval(GradPError(1:NI-1,1:NJ-1,:))

  end do

  GradPError=ABS(GradPExact - GradP)/GradPExact

  write(*,*) 'maximum x-error', maxval(GradPError(1:NI-1,1:NJ-1,1))
  write(*,*) 'maximum y-error', maxval(GradPError(1:NI-1,1:NJ-1,2))

!=== CALCULATE DIVERGENCE ===
  WRITE(*,*) 'Calculate Divergence'
  DivV=0.0
  Call B_CalcDivergence(NI,NJ,V,P,GradP,DivV,CellCenter,CellVolume,IFaceCenter,IFaceVector,JFaceCenter,JFaceVector,mode)
  DivVError=ABS(DivVExact-DivV)/DivVExact
  write(*,*) 'maximum error', maxval(DivVError(1:NI-1,1:NJ-1))

!=== CALCULATE LAPLACIAN ===
  WRITE(*,*) 'Calculate Laplacian'
  LapP=0.0
  Call B_CalcLaplacian(NI,NJ,P,GradP,LapP,CellCenter,CellVolume,IFaceCenter,IFaceVector,JFaceCenter,JFaceVector)
  LapPError=ABS(LapPExact-LapP)/LapPExact
  write(*,*) 'maximum error', maxval(LapPError(1:NI-1,1:NJ-1))

!=== OUTPUT FIELDS ===
  WRITE(*,*) 'Output fields to file: ', OutputFile
  Open(IO,FILE=OutputFile)
  Call B_OutputFields(IO,NI,NJ,X,Y,P,GradP,GradPError,V,DivV,DivVError,LapP,LapPError)
  Close(IO)

END PROGRAM Main