! This program solves the lid driven cavity flow problem. The user provides 
! the domain size in x and y directions, the mesh size, the mesh inflation factor,
! the preferred solution technique (Central Difference or Upwind), Reynold's Number 
! and lid velocity. Given this information the solution will generate the velocity
! field in both the x and y directions for the control volume.
! Author: Ahmad Osman, Student Number: 994700025.

module Dynamic_Array !Allows allocatable arrays to be passed into subroutines as an argument
    
    real*8, allocatable, dimension(:,:) :: grid_i_u, grid_f_u, grid_i_v, grid_f_v, grid_i_pc, grid_f_pc, grid_p
    !Initial grids contain initial guesses and are populated with the latest unconverged values
    !Final grids contain the computed values and are compared with the initial grids where values from them
    !are placed if they haven't converged allowing final grid to be overwritten until there's convergence
    !There are grids for x velocity (u), y velocity (v), pressure correction (pc) and pressure
    
    real*8, allocatable, dimension(:,:) :: ap, ae, aw, an, as, Fe, Fw, Fn, Fs, De, Dw, Dn, Ds, ue, uw, vn, vs
    !Arrays for each of the velocity/momentum coefficients ap, ae, aw, an and as
    !Arrays for the components of the coefficients F's and D's as well as face velocities used for F's
    
    real*8, allocatable, dimension(:,:) :: P_e, P_w, P_n, P_s, P_p, dx, dy, mdot, Sxp, Syp, Sp
    !Arrays for each of the pressure correction coefficients
    !Arrays for dx and dy across each control volume from face to face
    !Arrays for mdot and source terms used in the descritized Navier-Stokes equation (Sp and Su=>[Sxp, Syp])
   
    real*8, allocatable, dimension(:,:,:) ::  distance_grid, coordinate_grid
    !Distance grid contains non-uniform dx and dy values according to the inflation factor
    !Coordinate grid contains x and y coordinates for where the temperatures occur (CV centers)
    
    real*8, allocatable, dimension(:) :: Bj, Dj, aj, Cj, X, A_j, C_jP
    !Bj contains the first coefficients of the linear equations(bottom diagonal)
    !Dj contains the second coefficients of the linear equations(middle diagonal)
    !aj contains the third coefficients of the linear equations(top diagonal)
    !Cj contains elements '[B]' in the matrix equation '[A][x]=[B]' (right side of eqn)
    !X contains the computed results from the TDMA
    !A_j contains the first modified coefficients used in TDMA
    !C_jP contains the second modified coefficients used in TDMA
    
end module

program Cavity_Flows
    
    use Dynamic_Array !Contains initialization for allocatable arrays
    
    implicit none
    
    integer :: mesh_res_x, mesh_res_y, sol_type, box, EB, WB
    !One less than the number of nodes (how fine to discretize) in x and y
    !The choice of solution type, the presence of a box, an east border and a west border
    
    real*8 :: Lx, Ly, inf, Re, Lv, Iv, box_size
    !Respectively: size of global control volume Lx and Ly, inflation factor, 
    !Reynold's Number, lid velocity, inlet velocity and box size
    
    call Read_Input_File(Lx, Ly, inf, mesh_res_x, mesh_res_y, sol_type, Re, Lv, Iv, EB, WB, box, box_size)
    !Reads in and isolates all the variables for the problem from an input file
    
    call Generate_Geometry(mesh_res_x, mesh_res_y, Lx, Ly, inf, Lv, Iv)
    !Generates initial velocity and pressure guesses as well as boundary conditions, 
    !non-uniform values of dx and dy from point to point and from face to face according to inflation factor,
    !and generates a coordinate grid for the location of each velocity or pressure result
    
    call Solver(sol_type, mesh_res_x, mesh_res_y, Re, Lv, Iv, EB, WB, box, box_size)
    !Solves the velocity field using one of two solution methods (Central Difference or Upwind)

    call Output_Results(mesh_res_x, mesh_res_y)
    !Outputs the final array in a formatted text file readable by visualization software
    
    deallocate(grid_i_u, grid_f_u, grid_i_v, grid_f_v, grid_i_pc, grid_f_pc, grid_p, dx, dy, &
        ap, ae, aw, an, as, Fe, Fw, Fn, Fs, De, Dw, Dn, Ds, ue, uw, vn, vs, P_e, P_w, P_n, P_s, P_p,&
        distance_grid, coordinate_grid, Bj, Dj, aj, Cj, X, A_j, C_jP, Sxp, Syp, Sp, mdot)
    !Removes all the memory allocated for the various arrays

end program Cavity_Flows

subroutine Read_Input_File(Lx, Ly, inf, mesh_res_x, mesh_res_y, sol_type, Re, Lv, Iv, EB, WB, box, box_size)
    
    implicit none
    
    integer :: status, var_loc, mesh_res_x, mesh_res_y, sol_type, box, EB, WB
    !Assign variables for file status, to index strings, and store the mesh resolution,
    !solution type, and box & border presence parameters from the file
    
    real*8 :: Lx, Ly, inf, Re, Lv, Iv, box_size
    !Assign variables for each of the remaining expected input parameters from the file:
    !Size of the domain in x and y, inflation factor, lid velocity, inlet velocity and box size
    
    character(100) :: line, str, num_char
    !Assign variables for each line in the file, and to seperately store the character and number parts of the string
    
    open(1111, file="Input.txt")
    !Opens the input file with initial data for the problem 
          
    do !Loops through each line in the file until the end of the file is reached
        read(1111,'(a)', iostat=status) line !Reads the line including all spaces storing it in 'line'
        if (status < 0) then !Exits loop when end of file is reached
            exit
        else
            var_loc=index(line, "=") !Checks for what index the equal sign is at
            str=line(1:var_loc) !Stores the string up to the equal sign
            num_char=line(var_loc+1:len(line)) !Stores the string after the equal sign
            !Each if statement checks for one of the expected strings
            !If expect string is found it stores and converts the character after the equal sign into a number variable
            if (str .eq. "Size of Domain Lx =") then
                read(num_char,*)Lx
            else if (str .eq. "Size of Domain Ly =") then
                read(num_char,*)Ly
            else if (str .eq. "Mesh Resolution X =") then
                read(num_char,*)mesh_res_x
            else if (str .eq. "Mesh Resolution Y =") then
                read(num_char,*)mesh_res_y         
            else if (str .eq. "Non-Uniform Mesh Inflation Factor (Enter 1 for Uniform) =") then
                read(num_char,*)inf   
            else if (str .eq. "Solution Technique (Central:1, Upwind:2) =") then
                read(num_char,*)sol_type 
            else if (str .eq. "Reynold's Number =") then
                read(num_char,*)Re  
            else if (str .eq. "Lid Velocity =") then
                read(num_char,*)Lv            
            else if (str .eq. "Inlet Velocity =") then
                read(num_char,*)Iv          
            else if (str .eq. "Remove East Boundary (No:0, Yes:1) =") then
                read(num_char,*)EB        
            else if (str .eq. "Remove West Boundary (No:0, Yes:1) =") then
                read(num_char,*)WB         
            else if (str .eq. "Corner Box (No:0, Yes:1) =") then
                read(num_char,*)box
            else if (str .eq. "Size of Corner Box (Percentage of CV) =") then
                read(num_char,*)box_size
            end if
        end if
    end do
    
    close(1111)!Closes the opened file
    
    return
    
end subroutine Read_Input_File

subroutine Generate_Geometry(mesh_res_x, mesh_res_y, Lx, Ly, inf, Lv, Iv)
    
    use Dynamic_Array !Contains initialization for allocatable arrays
    
    implicit none
    
    integer :: mesh_res_x, mesh_res_y, i, j, count 
    !Variables for mesh resolution, rows, columns and an increment variable
    
    real*8 :: Lx, Ly, inf, Lv, Iv
    !Variables for domain size, inflation factor, lid velocity, and inlet velocity
    
    allocate(grid_i_u(1:mesh_res_y+1,1:mesh_res_x+1))
    allocate(grid_i_v(1:mesh_res_y+1,1:mesh_res_x+1))
    allocate(grid_i_pc(1:mesh_res_y+1,1:mesh_res_x+1))
    allocate(dx(1:mesh_res_y+1,1:mesh_res_x+1))
    allocate(dy(1:mesh_res_y+1,1:mesh_res_x+1))
    allocate(distance_grid(1:mesh_res_y,1:mesh_res_x,1:2))
    allocate(coordinate_grid(1:mesh_res_y+1,1:mesh_res_x+1,1:2))
    !Memory allocation for each of the arrays to be generated
    
    !Initializes the initial array with boundary velocities and guesses for velocities
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            if (i == mesh_res_y+1) then !North
                grid_i_u(i, j)=Lv !Lid velocity if present
                grid_i_v(i, j)=0
            else if (i == 1) then !South
                grid_i_u(i, j)=0
                grid_i_v(i, j)=0
            else if (j == 1) then !West
                grid_i_u(i, j)=Iv !Inlet velocity if present
                grid_i_v(i, j)=0
            else if (j == mesh_res_x+1) then !East
                grid_i_u(i, j)=0
                grid_i_v(i, j)=0
            else
                grid_i_u(i, j)=0
                grid_i_v(i, j)=0
            end if
        end do
    end do
    
    !Initializes guesses for pressure correction
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            grid_i_pc(i, j)=0
        end do
    end do
    
    if (inf == 1) then !Generates matrices for the size of dx and dy (uniform) for each control volume
        do i = 1, mesh_res_y
            do j = 1, mesh_res_x
                distance_grid(i,j,1)= Lx/mesh_res_x !Fill dx (point to point)
                distance_grid(i,j,2)= Ly/mesh_res_y !Fill dy (point to point)
            end do
        end do
    else !Generates matrices for the size of dx and dy (non-uniform) for each control volume
        do i = 1, mesh_res_y !Loop for non-uniform dx values
            count = 1 !Initialize count used to reverse the inflated sizes backwards for the second half
            do j = 1, mesh_res_x
                if (j >= 1 .and. j <= mesh_res_x/2) then
                    distance_grid(i,j,1)=((1-inf)/(1-inf**(mesh_res_x/2)))*(Lx/2)*(inf**(j-1)) !Fill dx forwards
                else if (j > mesh_res_x/2) then
                    distance_grid(i,j,1)= distance_grid(i, j-count, 1) !Fill dx backwards
                    count = count+2 !Count to subtract the j increment and the go to the previous element (1+1=2)
                end if
            end do
        end do
        do j = 1, mesh_res_x !Loop for non-uniform dy values
            count = 1 !Initialize count used to reverse the inflated sizes backwards for the second half
            do i = 1, mesh_res_y
                if (i >= 1 .and. i <= mesh_res_y/2) then
                    distance_grid(i,j,2)=((1-inf)/(1-inf**(mesh_res_y/2)))*(Ly/2)*(inf**(i-1)) !Fill dy forwards
                else if (i > mesh_res_y/2) then
                    distance_grid(i,j,2)= distance_grid(i-count, j, 2) !Fill dy backwards
                    count = count+2 !Count to subtract the j increment and the go to the previous element (1+1=2)
                end if
            end do
        end do
    end if
    
    !Generate matrix for corresponding X coordinate for each control volume
    do i = 1, mesh_res_y
        coordinate_grid(i,1,1)=0 !Initialize for the left column of X coordinates
            do j = 2, mesh_res_x+1
                coordinate_grid(i,j,1)=coordinate_grid(i,j-1,1)+distance_grid(i,j-1,1) !Generate each row of the X coordinate matrix
                coordinate_grid(mesh_res_y+1,j,1)=coordinate_grid(mesh_res_y,j,1) !Sets the last row of X coordinates to the row before
            end do
    end do

    !Generate matrix for corresponding Y coordinate for each control volume    
    do j = 1, mesh_res_x
        coordinate_grid(1,j,2)=0 !Initialize for the bottom row of Y coordinates
            do i = 2, mesh_res_y+1
                coordinate_grid(i,j,2)=coordinate_grid(i-1,j,2)+distance_grid(i-1,j,2) !Generate each column of the Y coordinate matrix
                coordinate_grid(i,mesh_res_x+1,2)=coordinate_grid(i,mesh_res_x,2) !Sets the last column of Y coordinates to the column before
            end do
    end do 
  
    !Generate matrix for corresponding dx value (face to face) for each control volume    
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            if (j == 1 .and. i == mesh_res_y+1) then
                dx(i, j)=distance_grid(i-1, j, 1)/2
            else if (j == mesh_res_x+1 .and. i == mesh_res_y+1) then
                dx(i, j)=distance_grid(i-1, j-1, 1)/2
            else if (i == mesh_res_y+1) then
                dx(i, j)=(distance_grid(i-1, j-1, 1)+distance_grid(i-1, j, 1))/2
            else if (j == 1) then
                dx(i, j)=distance_grid(i, j, 1)/2
            else if (j == mesh_res_x+1) then
                dx(i, j)=distance_grid(i, j-1, 1)/2
            else
                dx(i, j)=(distance_grid(i, j-1, 1)+distance_grid(i, j, 1))/2
            end if
        end do
    end do
    
    !Generate matrix for corresponding dy value (face to face) for each control volume    
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            if (i == 1 .and. j == mesh_res_x+1) then
                dy(i, j)=distance_grid(i, j-1, 2)/2 
            else if (i == mesh_res_y+1 .and. j == mesh_res_x+1) then
                dy(i, j)=distance_grid(i-1, j-1, 2)/2
            else if (j == mesh_res_x+1) then
                dy(i, j)=(distance_grid(i-1, j-1, 2)+distance_grid(i, j-1, 2))/2
            else if (i == 1) then
                dy(i, j)=distance_grid(i, j, 2)/2
            else if (i == mesh_res_y+1) then
                dy(i, j)=distance_grid(i-1, j, 2)/2
            else
                dy(i, j)=(distance_grid(i-1, j, 2)+distance_grid(i, j, 2))/2
            end if
        end do
    end do
    
    return
    
end subroutine Generate_Geometry

subroutine Solver(sol_type, mesh_res_x, mesh_res_y, Re, Lv, Iv, EB, WB, box, box_size)
    !Sets boundary conditions, initializes final velocity and pressure arrays,
    !Calls face velocity solver, calculates F and D coefficients,
    !Calls solution type subroutine to return 'a' coefficients using either central difference or upwind
    !Solves u and v momentum equations using initial conditions and the TDMA subroutine to convergence
    !Recalculates face velocities by calling that subroutine and updating F coefficients,
    !Calculates pressure correction coefficients and then solves for pressure correction using the TDMA
    !Corrects all velocity and final pressure values using the pressure correction, also corrects F's with them
    !Iterates from updating 'a' coefficients using new F's to the previous step 
    !until the value of mdot is minimized and also applies boundary conditions and the box if required
    
    use Dynamic_Array !Contains initialization for allocatable array
    
    implicit none

    real*8 :: Re, Lv, Iv, P_ref
    !Parameters for the Reynold's number, lid velocity, inlet velocity and reference pressure

    real*8 :: sum, avg, dif, fix, relax, box_size, compare
    !Sum of the differences between each element in the initial and final arrays
    !Average (avg) difference for all the elements in the initial and final arrays
    !The difference (dif) between two elements in the same postion in the initial and final array
    !The value for fixing (fix) certian values generally at the boundaries
    !Relaxation factor to ensure convergence
    !Box size to determine the size of the box relative to the domain
    !Compare variable to determine which control volume the box ends at
    
    integer :: sol_type, mesh_res_x, mesh_res_y, box, box_node_x, box_node_y, EB, WB
    !Respectively: stores the solution type, the mesh resolution in x and y directions, presence of box,
    !box length and width, presence of east and west boundaries
    
    integer :: num, count, count2, i, j
    !Respectively: counts the number of elements in initial/final array
    !counts the number of passes for inner and outer loops, row and column variables
            
    allocate(Bj(1:mesh_res_y+1)) !Bj = -a_n
    allocate(Dj(1:mesh_res_y+1)) !Dj = a_p = (a_n + a_s + a_w + a_e - s_p)
    allocate(aj(1:mesh_res_y+1)) !aj = -a_s
    allocate(Cj(1:mesh_res_y+1)) !Cj= a_w*T_w + a_e*T_e + s_u
    allocate(X(1:mesh_res_y+1)) !X = grid_f_u(column), X contains elements '[x]' in the matrix equation '[A][x]=[B]'
    allocate(A_j(1:mesh_res_y+1)) !A_j contains the first modified coefficients used in TDMA
    allocate(C_jP(1:mesh_res_y+1)) !C_jP contains the second modified coefficients used in TDMA
    allocate(grid_f_u(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains the final u velocity grid
    allocate(grid_f_v(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains the final v velocity grid
    allocate(grid_f_pc(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains the final pressure correction grid
    allocate(grid_p(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains the final pressure grid
    allocate(Fe(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains Fe coefficients for each node
    allocate(Fw(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains Fw coefficients for each node
    allocate(Fn(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains Fn coefficients for each node
    allocate(Fs(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains Fs coefficients for each node
    allocate(De(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains De coefficients for each node
    allocate(Dw(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains Dw coefficients for each node
    allocate(Dn(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains Dn coefficients for each node
    allocate(Ds(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains Ds coefficients for each node
    allocate(ap(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains ap coefficients for each node
    allocate(ae(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains ae coefficients for each node
    allocate(aw(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains aw coefficients for each node
    allocate(an(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains an coefficients for each node
    allocate(as(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains as coefficients for each node
    allocate(ue(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains ue coefficients for each node
    allocate(uw(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains uw coefficients for each node
    allocate(vn(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains vn coefficients for each node
    allocate(vs(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains vs coefficients for each node
    allocate(P_e(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains P'e coefficients for each node
    allocate(P_w(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains P'w coefficients for each node
    allocate(P_n(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains P'n coefficients for each node
    allocate(P_s(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains P's coefficients for each node
    allocate(P_p(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains P'p coefficients for each node
    allocate(mdot(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains mdot values for each node
    allocate(Sxp(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains pressure source term for u momentum equation
    allocate(Syp(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains pressure source term for v momentum equation
    allocate(Sp(1:mesh_res_y+1,1:mesh_res_x+1)) !Contains source term for ap
    
    open(4444, file="Iterations_to_Convergence.txt")!For outputing the number of iterations until convergence
    
    fix = 1e38 !Used to fix boundary velocities and box velocities
    relax = 0.7 !Relaxation factor used to help convergence
    
    !Sets the node at which the box ends along the x-axis  
    compare = 1
    i = 1
    do j = 1, mesh_res_x+1
        if(compare > 0) then
            compare = box_size - coordinate_grid(i, j, 1)/coordinate_grid(i, mesh_res_x+1, 1)
            box_node_x = j-1
        end if
    end do
    
    !Sets the node at which the box ends along the y-axis
    compare = 1
    j = 1
    do i = 1, mesh_res_y+1
        if (compare > 0) then
            compare = box_size - coordinate_grid(i, j, 2)/coordinate_grid(mesh_res_y+1, j, 2)
            box_node_y = i-1
        end if
    end do
    
    !Initializes the final array with boundary velocities and guesses for velocities
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            if (i == mesh_res_y+1) then !North
                grid_f_u(i, j)=Lv !Lid velocity if present
                grid_f_v(i, j)=0
            else if (i == 1) then !South
                grid_f_u(i, j)=0
                grid_f_v(i, j)=0
            else if (j == 1) then !West
                grid_f_u(i, j)=Iv !Inlet velocity if present
                grid_f_v(i, j)=0
            else if (j == mesh_res_x+1) then !East
                grid_f_u(i, j)=0
                grid_f_v(i, j)=0
            else !Inner Domain
                grid_f_u(i, j)=0
                grid_f_v(i, j)=0
            end if
        end do
    end do
    
    !Initializes the final arrays with guesses for pressure, pressure correction and ap
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            grid_f_pc(i, j)=0
            grid_p(i, j)=0
            ap(i, j)=1    
        end do
    end do

    call u_v_Solver (mesh_res_x, mesh_res_y, Lv, EB, WB)!(step 1 of SIMPLE Collocated Algorithm)
    
    !Populate arrays for De, Dw, Dn, and Ds
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
                !==========================================================!
            if (i == 1 .and. j == 1 ) then !South West corner
                !==========================================================!
                De(i, j)=dy(i, j)/(Re*distance_grid(i, j, 1))
                Dw(i, j)=0 !No West diffusion
                Dn(i, j)=dx(i, j)/(Re*distance_grid(i, j, 2))
                Ds(i, j)=0 !No South diffusion
                !==========================================================!
            else if (i == 1 .and. j == mesh_res_x+1) then !South East corner
                !==========================================================!
                De(i, j)=0 !No East diffusion
                Dw(i, j)=dy(i, j)/(Re*distance_grid(i, j-1, 1))
                Dn(i, j)=dx(i, j)/(Re*distance_grid(i, j-1, 2))
                Ds(i, j)=0 !No South diffusion
                !==========================================================!
            else if (i == mesh_res_y+1 .and. j == 1) then !North West corner
                !==========================================================!
                De(i, j)=dy(i, j)/(Re*distance_grid(i-1, j, 1))
                Dw(i, j)=0 !No West diffusion
                Dn(i, j)=0 !No North diffusion
                Ds(i, j)=dx(i, j)/(Re*distance_grid(i-1, j, 2))
                !==========================================================!
            else if (i == mesh_res_y+1 .and. j == mesh_res_x+1) then !North East corner
                !==========================================================!
                De(i, j)=0 !No East diffusion
                Dw(i, j)=dy(i, j)/(Re*distance_grid(i-1, j-1, 1))
                Dn(i, j)=0 !No North diffusion
                Ds(i, j)=dx(i, j)/(Re*distance_grid(i-1, j-1, 2))
                !==========================================================!
            else if (i == 1 .and. j /= 1 .and. j /= mesh_res_x+1) then !South border without corners
                !==========================================================!
                De(i, j)=dy(i, j)/(Re*distance_grid(i, j, 1))
                Dw(i, j)=dy(i, j)/(Re*distance_grid(i, j-1, 1))
                Dn(i, j)=dx(i, j)/(Re*distance_grid(i, j, 2))
                Ds(i, j)=0 !No South diffusion
                !==========================================================!
            else if (j == 1 .and. i /= 1 .and. 1 /= mesh_res_y+1) then !West border without corners
                !==========================================================!
                De(i, j)=dy(i, j)/(Re*distance_grid(i, j, 1))
                Dw(i, j)=0 !No West diffusion
                Dn(i, j)=dx(i, j)/(Re*distance_grid(i, j, 2))
                Ds(i, j)=dx(i, j)/(Re*distance_grid(i-1, j, 2))
                !==========================================================!
            else if (i == mesh_res_y+1 .and. j /= 1 .and. j /= mesh_res_x+1) then !North border without corners
                !==========================================================!
                De(i, j)=dy(i, j)/(Re*distance_grid(i, j, 1))
                Dw(i, j)=dy(i, j)/(Re*distance_grid(i, j-1, 1))
                Dn(i, j)=0 !No North diffusion
                Ds(i, j)=dx(i, j)/(Re*distance_grid(i-1, j, 2))
                !==========================================================!
            else if (j == mesh_res_x+1 .and. i /= 1 .and. i /= mesh_res_y+1) then !East border without corners
                !==========================================================!
                De(i, j)=0 !No East diffusion
                Dw(i, j)=dy(i, j)/(Re*distance_grid(i, j-1, 1))
                Dn(i, j)=dx(i, j)/(Re*distance_grid(i, j-1, 2))
                Ds(i, j)=dx(i, j)/(Re*distance_grid(i-1, j-1, 2))
                !==========================================================!
            else !Remaining control volume (Center)
                !==========================================================!
                De(i, j)=dy(i, j)/(Re*distance_grid(i, j, 1))
                Dw(i, j)=dy(i, j)/(Re*distance_grid(i, j-1, 1))
                Dn(i, j)=dx(i, j)/(Re*distance_grid(i, j, 2))
                Ds(i, j)=dx(i, j)/(Re*distance_grid(i-1, j, 2))
                !==========================================================!
            end if 
        end do
    end do
    
    !Populate arrays for Fe, Fw, Fn, Fs
    do i = 1, mesh_res_y+1 
        do j = 1, mesh_res_x+1
            !==========================================================!
            Fe(i, j)=ue(i, j)*dy(i,j)
            Fw(i, j)=uw(i, j)*dy(i,j)
            Fn(i, j)=vn(i, j)*dx(i,j)
            Fs(i, j)=vs(i, j)*dx(i,j)
            !==========================================================!
        end do
    end do
        
    count2 = 0 !Initialize iteration count for major loop

    do while (count2 < 400)     !>>>>>>> MAJOR ITERATION LOOP <<<<<<<!
        
        call Solution_Method(sol_type, mesh_res_x, mesh_res_y)
        
        !Calculate the pressure source terms for the u velocity
        do i = 1, mesh_res_y+1
            do j = 1, mesh_res_x+1
                !==========================================================!
                if (j == 1) then !West Border
                    Sxp(i, j)=(0.5*(grid_p(i, j+1)-grid_p(i, j)))*dy(i, j)
                else if (j == mesh_res_x+1) then !East Border
                    Sxp(i, j)=(0.5*grid_p(i, j)-0.5*grid_p(i, j-1))*dy(i, j)
                else
                    Sxp(i, j)=(0.5*grid_p(i, j+1)-0.5*grid_p(i, j-1))*dy(i, j)
                end if
                !==========================================================!
            end do
        end do 
        
        !Calculate the pressure source terms for the v velocity
        do i = 1, mesh_res_y+1
            do j = 1, mesh_res_x+1
                !==========================================================!
                if (i == 1) then !South Border
                    Syp(i, j)=(0.5*(grid_p(i+1, j)-grid_p(i, j)))*dx(i, j)
                else if (i == mesh_res_y+1) then !North Border
                    Syp(i, j)=(0.5*(grid_p(i, j)-grid_p(i-1, j)))*dx(i, j)
                else
                    Syp(i, j)=(0.5*(grid_p(i+1, j)-grid_p(i-1, j)))*dx(i, j)
                end if
                !==========================================================!
            end do
        end do 
        
        !Apply relaxation to ap
        do i = 1, mesh_res_y+1
            do j = 1, mesh_res_x+1
                !==========================================================!
                ap(i, j)=ap(i, j)/relax
                !==========================================================!
            end do
        end do     
        
        !Apply relaxation to Sxp and Syp
        do i = 2, mesh_res_y
            do j = 2, mesh_res_x
                !==========================================================!
                Sxp(i, j)=Sxp(i, j)-(1-relax)*ap(i, j)*grid_f_u(i, j)
                Syp(i, j)=Syp(i, j)-(1-relax)*ap(i, j)*grid_f_v(i, j)
                !==========================================================!
            end do
        end do 
        
        !Re-initialize the final array with boundary velocities and guesses for velocities
        do i = 1, mesh_res_y+1
            do j = 1, mesh_res_x+1
                !==========================================================!
                if (i == mesh_res_y+1) then !North
                    grid_f_u(i, j)=Lv !Lid velocity if present
                    grid_f_v(i, j)=0
                else if (i == 1) then !South
                    grid_f_u(i, j)=0
                    grid_f_v(i, j)=0
                else if (j == 1) then !West
                    if (WB == 1) then !If West border not present then do not overwrite solver answers
                        grid_f_u(i, j)=grid_f_u(i, j)
                        grid_f_v(i, j)=grid_f_v(i, j)
                    else if (WB == 0) then !If West border or inlet present fix velocity values
                        grid_f_u(i, j)=Iv !Inlet velocity if present
                        grid_f_v(i, j)=0
                    end if
                else if (j == mesh_res_x+1) then !East
                    if (EB == 1) then !If East border not present then do not overwrite solver answers
                        grid_f_u(i, j)=grid_f_u(i, j)
                        grid_f_v(i, j)=grid_f_v(i, j)
                    else if (EB == 0) then !If East border present fix velocity values
                        grid_f_u(i, j)=0
                        grid_f_v(i, j)=0
                    end if
                end if
                !==========================================================!
            end do
        end do
        
        !Set velocities inside the box to 0
        if (box == 1) then
            do i = 1, box_node_x
                do j = 1, box_node_y
                    grid_f_u(i, j)=0
                    grid_f_v(i, j)=0
                end do
            end do
        end if
        
        !Assign values to Sp, Sxp and Syp to fix border velocities
        do i = 1, mesh_res_y+1
            do j = 1, mesh_res_x+1
                if (i == 1 .or. i == mesh_res_y+1 .or. j == 1 .or. j == mesh_res_x+1) then
                    if (EB == 1 .and. WB == 1) then !Remove East and West boundary fix
                        if (i == 1 .or. i == mesh_res_y+1) then !Fix only North and South
                            Sp(i, j)=fix
                            Sxp(i, j)=Sxp(i, j)-fix*grid_f_u(i, j)
                            Syp(i, j)=Syp(i, j)-fix*grid_f_v(i, j)
                        end if
                    else if (EB == 1) then !Remove East boundary fix
                        if (i == 1 .or. i == mesh_res_y+1 .or. j == 1) then !Dont fix East
                            Sp(i, j)=fix
                            Sxp(i, j)=Sxp(i, j)-fix*grid_f_u(i, j)
                            Syp(i, j)=Syp(i, j)-fix*grid_f_v(i, j)
                        end if
                    else if (WB == 1) then !Remove West boundary fix
                        if (i == 1 .or. i == mesh_res_y+1 .or. j == mesh_res_x+1) then !Dont fix West
                            Sp(i, j)=fix
                            Sxp(i, j)=Sxp(i, j)-fix*grid_f_u(i, j)
                            Syp(i, j)=Syp(i, j)-fix*grid_f_v(i, j)
                        end if
                    else !Fix all boundaries
                        Sp(i, j)=fix
                        Sxp(i, j)=Sxp(i, j)-fix*grid_f_u(i, j)
                        Syp(i, j)=Syp(i, j)-fix*grid_f_v(i, j)
                    end if
                else
                    Sp(i, j)=0
                end if
            end do
        end do
        
        !Modify the coefficients in the box to fix velocities
        if (box == 1) then
            do i = 1, box_node_x
                do j = 1, box_node_y
                    Sp(i, j)=fix
                    Sxp(i, j)=Sxp(i, j)-fix*grid_f_u(i, j)
                    Syp(i, j)=Syp(i, j)-fix*grid_f_v(i, j)
                end do
            end do
        end if
        
        !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>Solution for u Velocity<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<!(step 2 of SIMPLE Collocated Algorithm)
        
        count=0 !Initializes the counting variable for number of scans
        avg=10 !Initializes the average difference so the loop can be entered
        
        do while (avg > 1e-11) !(count < 100)!Keeps looping until average difference between the initial and final elements is less than stated
            do j=1, mesh_res_x+1
                do i=1, mesh_res_y+1
                    if (j == 1) then !West border
                        aj(i)=-an(i, j)
                        !-a_n
                        Dj(i)=ap(i,j)+Sp(i, j)
                        !a_p = an + as + aw + ae + Fe - Fw + Fn - Fs
                        Bj(i)=-as(i, j)
                        !-a_s
                        Cj(i)=ae(i, j)*grid_f_u(i, j+1)-Sxp(i, j)!+fix*grid_f_u(i, j)
                        !Cj= aw*uw + ae*ue - (0.5*Pe - 0.5*Pp)
                    !==========================================================!
                    else if (j == mesh_res_x+1) then !East border
                        aj(i)=-an(i, j)
                        !-a_n
                        Dj(i)=ap(i,j)+Sp(i, j)
                        !a_p = an + as + aw + ae + Fe - Fw + Fn - Fs
                        Bj(i)=-as(i, j)
                        !-a_s
                        Cj(i)=aw(i, j)*grid_f_u(i, j-1)-Sxp(i, j)!+fix*grid_f_u(i, j)
                        !Cj= aw*uw + ae*ue - (0.5*Pp - 0.5*Pw)
                    !==========================================================!
                    else
                        aj(i)=-an(i, j)
                        !-a_n
                        Dj(i)=ap(i,j)+Sp(i, j)
                        !a_p = an + as + aw + ae + Fe - Fw + Fn - Fs
                        Bj(i)=-as(i, j)
                        !-a_s
                        Cj(i)=aw(i, j)*grid_f_u(i, j-1)+ae(i, j)*grid_f_u(i, j+1)-Sxp(i, j)
                        !Cj= aw*uw + ae*ue - (0.5*Pe - 0.5*Pw)
                    !==========================================================!
                    end if
                end do
                !==========================================================!
                call Tridiagonal_Matrix_Solver(mesh_res_y)
                !==========================================================!
                do i=2, mesh_res_y
                    grid_f_u(i, j)=X(i)
                end do
                !==========================================================!
            end do
            !==========================================================!
            sum=0 !Initializes the differences between each element in the initial and final arrays
            num=0 !Initializes the number of elements in initial/final array
            !==========================================================!
            !Finds the difference between each element, adds the difference to the sum of all the differences
            !Increments the number of elements differenced
            !Overwrites the initial grid after the compare which then allows the final grid to be overwritten in the next pass
            do i = 2, mesh_res_y
                do j = 2, mesh_res_x
                    dif=abs(grid_f_u(i, j)-grid_i_u(i,j))
                    sum=sum + dif
                    num=num + 1
                    grid_i_u(i, j)=grid_f_u(i, j)
                end do
            end do
            !==========================================================!
            avg=abs(sum/num) !Computes the absolute value of the average difference
            !print*, "U Avg:", avg
            count=count+1 !Increments the number of full iterations
            !print*, "U Iteration Number:", count
            !==========================================================!
        end do
        
        write (4444, *) (count) !Prints the number of iterations to a file
        write (4444, *) "u"
        
        !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>Solution for v Velocity<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<!(step 2 of SIMPLE Collocated Algorithm)
        
        count=0 !Initializes the counting variable for number of scans
        avg=10 !Initializes the average difference so the loop can be entered
        
        do while (avg > 1e-11) !Keeps looping until average difference between the initial and final elements is less than stated
            do j=1, mesh_res_x+1
                !==========================================================!
                do i=1, mesh_res_y+1
                    if (j == 1) then !West border
                        aj(i)=-an(i, j)
                        !-a_n
                        Dj(i)=ap(i,j)+Sp(i, j)
                        !a_p = [a_n + a_s + a_w + a_e - sp(-(sp) => Fs)]+Fe-Fw+Fn-Fs
                        Bj(i)=-as(i, j)
                        !-a_s
                        Cj(i)=ae(i, j)*grid_f_v(i, j+1)-Syp(i, j)!+fix*grid_f_v(i, j)
                        !Cj= ae*ve - (0.5*Pn - 0.5*Pp)
                    !==========================================================!
                    else if (j == mesh_res_y+1) then !East border
                        aj(i)=-an(i, j)
                        !-a_n
                        Dj(i)=ap(i,j)+Sp(i, j)
                        !a_p = [a_n + a_s + a_w + a_e - sp(-(sp) => Fs)]+Fe-Fw+Fn-Fs
                        Bj(i)=-as(i, j)
                        !-a_s
                        Cj(i)=aw(i, j)*grid_f_v(i, j-1)-Syp(i, j)!+fix*grid_f_v(i, j)
                        !Cj= aw*vw + ae*ve - (0.5*Pp - 0.5*Ps)
                    !==========================================================!
                    else
                        aj(i)=-an(i, j)
                        !-a_n
                        Dj(i)=ap(i,j)+Sp(i, j)
                        !a_p = [a_n + a_s + a_w + a_e - sp(-(sp) => Fs)]+Fe-Fw+Fn-Fs
                        Bj(i)=-as(i, j)
                        !-a_s
                        Cj(i)=aw(i, j)*grid_f_v(i, j-1)+ae(i, j)*grid_f_v(i, j+1)-Syp(i, j)
                        !Cj= aw*vw + ae*ve - (0.5*Pn - 0.5*Ps)
                    end if
                end do
                !==========================================================!
                call Tridiagonal_Matrix_Solver(mesh_res_y)
                !==========================================================!
                do i=2, mesh_res_y
                    grid_f_v(i, j)=X(i)
                end do
                !==========================================================!
            end do
            !==========================================================!
            sum=0 !Initializes the differences between each element in the initial and final arrays
            num=0 !Initializes the number of elements in initial/final array
            !==========================================================!
            !Finds the difference between each element, adds the difference to the sum of all the differences
            !Increments the number of elements differenced
            !Overwrites the initial grid after the compare which then allows the final grid to be overwritten in the next pass
            do i = 2, mesh_res_y
                do j = 2, mesh_res_x
                    dif=abs(grid_f_v(i, j)-grid_i_v(i,j))
                    sum=sum + dif
                    num=num + 1
                    grid_i_v(i, j)=grid_f_v(i, j)
                end do
            end do
            !==========================================================!
            avg=abs(sum/num) !Computes the absolute value of the average difference
            !print*, "V Avg:", avg
            count=count+1 !Increments the number of full iterations
            !print*, "V Iteration Number:", count
            !==========================================================!
        end do
            
        write (4444, *) (count) !Prints the number of iterations to a file
        write (4444, *) "v"
        
        call u_v_Solver (mesh_res_x, mesh_res_y, Lv, EB, WB)

        !Update arrays for Fe, Fw, Fn, Fs using updated u's and v's (step 3 of SIMPLE Collocated Algorithm)
        do i = 1, mesh_res_y+1
            do j = 1, mesh_res_x+1
                !==========================================================!
                Fe(i, j)=ue(i, j)*dy(i,j)
                Fw(i, j)=uw(i, j)*dy(i,j)
                Fn(i, j)=vn(i, j)*dx(i,j)
                Fs(i, j)=vs(i, j)*dx(i,j)
                !==========================================================!
            end do
        end do
        
        !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>Solution for P' Pressure Correction<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<!(step 4 of SIMPLE Collocated Algorithm)
        
        !Populate arrays for P'e, P'w, P'n, P's, P'p
        do i = 1, mesh_res_y+1
            do j = 1, mesh_res_x+1
                !==========================================================!
                if (i == 1 .and. j == 1) then !South West corner
                !==========================================================!   
                    P_e(i, j)=(dy(i, j)**2)/((ap(i, j+1)+ap(i, j))/2)
                    P_w(i, j)=0
                    P_n(i, j)=(dx(i, j)**2)/((ap(i+1, j)+ap(i, j))/2)
                    P_s(i, j)=0
                    P_p(i, j)=P_e(i, j)+P_w(i, j)+P_n(i, j)+P_s(i, j)
                !==========================================================!
                else if (i == 1 .and. j == mesh_res_x+1) then !South East corner
                !==========================================================!   
                    P_e(i, j)=0
                    P_w(i, j)=(dy(i, j)**2)/((ap(i, j-1)+ap(i, j))/2)
                    P_n(i, j)=(dx(i, j)**2)/((ap(i+1, j)+ap(i, j))/2)
                    P_s(i, j)=0
                    P_p(i, j)=P_e(i, j)+P_w(i, j)+P_n(i, j)+P_s(i, j)
                !==========================================================!
                else if (i == mesh_res_y+1 .and. j == 1) then !North West corner
                !==========================================================!   
                    P_e(i, j)=(dy(i, j)**2)/((ap(i, j+1)+ap(i, j))/2)
                    P_w(i, j)=0
                    P_n(i, j)=0
                    P_s(i, j)=(dx(i, j)**2)/((ap(i-1, j)+ap(i, j))/2)
                    P_p(i, j)=P_e(i, j)+P_w(i, j)+P_n(i, j)+P_s(i, j)
                !==========================================================!
                else if (i == mesh_res_y+1 .and. j == mesh_res_x+1) then !North East corner
                !==========================================================!   
                    P_e(i, j)=0
                    P_w(i, j)=(dy(i, j)**2)/((ap(i, j-1)+ap(i, j))/2)
                    P_n(i, j)=0
                    P_s(i, j)=(dx(i, j)**2)/((ap(i-1, j)+ap(i, j))/2)
                    P_p(i, j)=P_e(i, j)+P_w(i, j)+P_n(i, j)+P_s(i, j)
                !==========================================================!
                else if (i == 1) then !South Border
                !==========================================================!   
                    P_e(i, j)=(dy(i, j)**2)/((ap(i, j+1)+ap(i, j))/2)
                    P_w(i, j)=(dy(i, j)**2)/((ap(i, j-1)+ap(i, j))/2)
                    P_n(i, j)=(dx(i, j)**2)/((ap(i+1, j)+ap(i, j))/2)
                    P_s(i, j)=0
                    P_p(i, j)=P_e(i, j)+P_w(i, j)+P_n(i, j)+P_s(i, j)
                !==========================================================!
                else if (i == mesh_res_y+1) then !North Border
                !==========================================================!   
                    P_e(i, j)=(dy(i, j)**2)/((ap(i, j+1)+ap(i, j))/2)
                    P_w(i, j)=(dy(i, j)**2)/((ap(i, j-1)+ap(i, j))/2)
                    P_n(i, j)=0
                    P_s(i, j)=(dx(i, j)**2)/((ap(i-1, j)+ap(i, j))/2)
                    P_p(i, j)=P_e(i, j)+P_w(i, j)+P_n(i, j)+P_s(i, j)
                !==========================================================!
                else if (j == 1) then !West Border
                !==========================================================!   
                    P_e(i, j)=(dy(i, j)**2)/((ap(i, j+1)+ap(i, j))/2)
                    P_w(i, j)=0
                    P_n(i, j)=(dx(i, j)**2)/((ap(i+1, j)+ap(i, j))/2)
                    P_s(i, j)=(dx(i, j)**2)/((ap(i-1, j)+ap(i, j))/2)
                    P_p(i, j)=P_e(i, j)+P_w(i, j)+P_n(i, j)+P_s(i, j)
                !==========================================================!
                else if (j == mesh_res_x+1) then !East Border
                !==========================================================!   
                    P_e(i, j)=0
                    P_w(i, j)=(dy(i, j)**2)/((ap(i, j-1)+ap(i, j))/2)
                    P_n(i, j)=(dx(i, j)**2)/((ap(i+1, j)+ap(i, j))/2)
                    P_s(i, j)=(dx(i, j)**2)/((ap(i-1, j)+ap(i, j))/2)
                    P_p(i, j)=P_e(i, j)+P_w(i, j)+P_n(i, j)+P_s(i, j)
                !==========================================================!
                else !Inner Pressure Corrections
                !==========================================================!   
                    P_e(i, j)=(dy(i, j)**2)/((ap(i, j+1)+ap(i, j))/2)
                    P_w(i, j)=(dy(i, j)**2)/((ap(i, j-1)+ap(i, j))/2)
                    P_n(i, j)=(dx(i, j)**2)/((ap(i+1, j)+ap(i, j))/2)
                    P_s(i, j)=(dx(i, j)**2)/((ap(i-1, j)+ap(i, j))/2)
                    P_p(i, j)=P_e(i, j)+P_w(i, j)+P_n(i, j)+P_s(i, j)
                !==========================================================!
                end if
            end do
        end do
        
        count=0 !Initializes the counting variable for number of scans
        avg=10 !Initializes the average difference so the loop can be entered
        
        do while (avg > 1e-11) !Keeps looping until average difference between the initial and final elements is less than stated
            do j=1, mesh_res_x+1
                !==========================================================!
                do i=1, mesh_res_y+1
                    if (j == 1) then !West border
                        aj(i)=-P_n(i, j)
                        !-a_n
                        Dj(i)=P_p(i,j)
                        !a_p = an + as + aw + ae
                        Bj(i)=-P_s(i, j)
                        !-a_s
                        Cj(i)=P_e(i, j)*grid_f_pc(i, j+1)&
                            +(uw(i, j)-ue(i, j))*dy(i, j)+(vs(i, j)-vn(i, j))*dx(i, j)
                        !Cj= Pw*pc_w + Pe*pc_e + mdot
                    else if (j == mesh_res_x+1) then !East border
                        aj(i)=-P_n(i, j)
                        !-a_n
                        Dj(i)=P_p(i,j)
                        !a_p = an + as + aw + ae
                        Bj(i)=-P_s(i, j)
                        !-a_s
                        Cj(i)=P_w(i, j)*grid_f_pc(i, j-1)&
                            +(uw(i, j)-ue(i, j))*dy(i, j)+(vs(i, j)-vn(i, j))*dx(i, j)
                        !Cj= Pw*pc_w + Pe*pc_e + mdot
                    else
                        aj(i)=-P_n(i, j)
                        !-a_n
                        Dj(i)=P_p(i,j)
                        !a_p = an + as + aw + ae
                        Bj(i)=-P_s(i, j)
                        !-a_s
                        Cj(i)=P_w(i, j)*grid_f_pc(i, j-1)+P_e(i, j)*grid_f_pc(i, j+1)&
                            +(uw(i, j)-ue(i, j))*dy(i, j)+(vs(i, j)-vn(i, j))*dx(i, j)
                        !Cj= Pw*pc_w + Pe*pc_e + mdot
                    end if
                end do
                !==========================================================!
                call Tridiagonal_Matrix_Solver(mesh_res_y)
                !==========================================================!
                do i=1, mesh_res_y+1
                    grid_f_pc(i, j)=X(i)
                end do
                !==========================================================!
            end do
            !==========================================================!
            sum=0 !Initializes the differences between each element in the initial and final arrays
            num=0 !Initializes the number of elements in initial/final array
            !==========================================================!
            !Finds the difference between each element, adds the difference to the sum of all the differences
            !Increments the number of elements differenced
            !Overwrites the initial grid after the compare which then allows the final grid to be overwritten in the next pass
            do i = 1, mesh_res_y+1
                do j = 1, mesh_res_x+1
                    dif=abs(grid_f_pc(i, j)-grid_i_pc(i,j))
                    sum=sum + dif
                    num=num + 1
                    grid_i_pc(i, j)=grid_f_pc(i, j)
                end do
            end do
            !==========================================================!
            avg=abs(sum/num) !Computes the absolute value of the average difference
            !print*, "P Avg:", avg
            count=count+1 !Increments the number of full iterations
            !print*, "P Iteration Number:", count
            !==========================================================!
        end do
        
        write (4444, *) (count) !Prints the number of iterations to a file
        write (4444, *) "P"
        
        !Correct u, v, and F's using pressure corrections (step 5 & 6 of SIMPLE Collocated Algorithm)
        do i=2, mesh_res_y
            do j=2, mesh_res_x
                grid_f_u(i, j)=(dy(i, j)/ap(i, j))*(grid_f_pc(i, j-1)-grid_f_pc(i, j+1))*0.5+grid_f_u(i, j)
                grid_f_v(i, j)=(dx(i, j)/ap(i, j))*(grid_f_pc(i-1, j)-grid_f_pc(i+1, j))*0.5+grid_f_v(i, j)
                Fe(i, j)=Fe(i, j)+(dy(i, j)**2/((ap(i, j+1)+ap(i, j))/2))*(grid_f_pc(i, j)-grid_f_pc(i, j+1))
                Fw(i, j)=Fw(i, j)+(dy(i, j)**2/((ap(i, j-1)+ap(i, j))/2))*(grid_f_pc(i, j-1)-grid_f_pc(i, j))
                Fn(i, j)=Fn(i, j)+(dx(i, j)**2/((ap(i+1, j)+ap(i, j))/2))*(grid_f_pc(i, j)-grid_f_pc(i+1, j))
                Fs(i, j)=Fs(i, j)+(dx(i, j)**2/((ap(i-1, j)+ap(i, j))/2))*(grid_f_pc(i-1, j)-grid_f_pc(i, j))
            end do
        end do
        
        !Sets box velocities to 0
        if (box == 1) then
            do i=1, mesh_res_y
                do j=1, mesh_res_x
                    if (i >= 1 .and. i <= box_node_x .and. j >= 1 .and. j <= box_node_y) then
                        grid_f_u(i, j)= 0
                        grid_f_v(i, j)= 0
                    end if
                end do
            end do
        end if
        
        !Correct P using pressure corrections (step 5 of SIMPLE Collocated Algorithm)
        P_ref = grid_f_pc((mesh_res_y+2)/2, (mesh_res_x+2)/2)
        do i=1, mesh_res_y+1
            do j=1, mesh_res_x+1
                grid_p(i, j)=0.1*(grid_f_pc(i, j)-P_ref)+grid_p(i, j)
            end do
        end do 

        !==========================================================!
        sum=0 !Initializes the sum of the mdot for each element
        num=0 !Initializes the number of elements in final array
        !==========================================================!
        !Finds the difference between each element, adds the difference to the sum of all the differences
        !Increments the number of elements differenced
        !Overwrites the initial grid after the compare which then allows the final grid to be overwritten in the next pass
        do i = 1, mesh_res_y+1
            do j = 1, mesh_res_x+1
                mdot(i, j)=(uw(i, j)-ue(i, j))*dy(i, j)+(vs(i, j)-vn(i, j))*dx(i, j)
                sum = sum + abs(mdot(i, j))
                num = num + 1
            end do
        end do
        !==========================================================!
        avg=abs(sum/num) !Computes the absolute value of the average difference
        count2=count2+1 !Increments the number of full iterations
        print*, "Major Iteration Number:", count2
        print*, "Average Mass Change:", avg
        !==========================================================!
    
    end do
    
    close (4444) !Closes the file
    
    pause
    
    return

end subroutine Solver

subroutine Solution_Method(sol_type, mesh_res_x, mesh_res_y)
    
    implicit none
    
    integer :: mesh_res_x, mesh_res_y, sol_type !mesh resolution and solution type
    
    if (sol_type == 1) then
        call Difference_Scheme(mesh_res_x, mesh_res_y)
    else if (sol_type == 2) then
        call Upwind_Scheme(mesh_res_x, mesh_res_y)
    end if
    
end subroutine Solution_Method

subroutine Difference_Scheme(mesh_res_x, mesh_res_y)
    
    use Dynamic_Array !Contains initialization for allocatable array
    
    implicit none
        
    integer :: mesh_res_x, mesh_res_y, i, j
            
    !Populate arrays for ae, aw, an, as, and ap (Updates coefficients when new F's are calculated)
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            !==========================================================!
            ae(i, j)=De(i, j)-0.5*Fe(i, j)
            aw(i, j)=Dw(i, j)+0.5*Fw(i, j)
            an(i, j)=Dn(i, j)-0.5*Fn(i, j)
            as(i, j)=Ds(i, j)+0.5*Fs(i, j)
            ap(i, j)=ae(i, j)+aw(i, j)+an(i, j)+as(i, j)!+Fe(i, j)-Fw(i, j)+Fn(i, j)-Fs(i, j)
            !==========================================================!
        end do
    end do
        
end subroutine Difference_Scheme

subroutine Upwind_Scheme(mesh_res_x, mesh_res_y)
    
    use Dynamic_Array !Contains initialization for allocatable array
        
    implicit none
    
    integer :: mesh_res_x, mesh_res_y, i, j
    
    !Populate arrays for ae, aw, an, as, and ap (Updates coefficients when new F's are calculated)
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            !==========================================================!
            ae(i, j)=De(i, j)+max(0,-Fe(i, j))
            aw(i, j)=Dw(i, j)+max(0,Fw(i, j))
            an(i, j)=Dn(i, j)+max(0,-Fn(i, j))
            as(i, j)=Ds(i, j)+max(0,Fs(i, j))
            ap(i, j)=ae(i, j)+aw(i, j)+an(i, j)+as(i, j)!+Fe(i, j)-Fw(i, j)+Fn(i, j)-Fs(i, j)
            !==========================================================!
        end do
    end do

end subroutine Upwind_Scheme

subroutine u_v_Solver (mesh_res_x, mesh_res_y, Lv, EB, WB)
    !Solves for face velocities u and v at the faces for each control volume (step 1 of SIMPLE Collocated Algorithm)
    
    use Dynamic_Array !Contains initialization for allocatable arrays
    
    implicit none
    
    real*8 :: Lv !lid velocity
    integer :: mesh_res_x, mesh_res_y, EB, WB !mesh resolution
    integer :: i, j !row and column variables
    
     !Populate arrays for ue, uw, vn, and vs
     do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
                !==========================================================!
            if (i == 1 .and. j == 1 ) then !South-West corner
                !==========================================================!
                ue(i, j)=0 !Horizontal BC along South Boundary
                uw(i, j)=0 !No West face at west nodes
                vn(i, j)=0 !Vertical BC along West Boundary
                vs(i, j)=0 !No South face at south nodes
                !==========================================================!
            else if (i == 1 .and. j == mesh_res_x+1) then !South-East corner
                !==========================================================!
                ue(i, j)=0 !No South face at south nodes
                uw(i, j)=0 !Horizontal BC along South Boundary
                vn(i, j)=0 !Vertical BC along East Boundary
                vs(i, j)=0 !No South face at south nodes
                !==========================================================!
            else if (i == mesh_res_y+1 .and. j == 1) then !North-West corner
                !==========================================================!
                ue(i, j)=Lv !Horizontal BC along North Boundary
                uw(i, j)=0 !No West face at west nodes
                vn(i, j)=0 !No North face at north nodes
                vs(i, j)=0 !Vertical BC along West Boundary
                !==========================================================!
            else if (i == mesh_res_y+1 .and. j == mesh_res_x+1) then !North-East corner
                !==========================================================!
                ue(i, j)=0 !No East face at east nodes
                uw(i, j)=Lv !Horizontal BC along North Boundary
                vn(i, j)=0 !No North face at north nodes
                vs(i, j)=0 !Vertical BC along East Boundary
                !==========================================================!
            else if (i == 1 .and. j /= 1 .and. j /= mesh_res_x+1) then !South border without corner points
                !==========================================================!
                ue(i, j)=0 !Horizontal BC along South Boundary
                uw(i, j)=0 !Horizontal BC along South Boundary
                vn(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i+1, j))&
                    +0.5*((dx(i+1, j)/ap(i+1, j))*0.5*(grid_p(i+2, j)-grid_p(i, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i, j)))& !PN-PS to PN-PP
                    -(dx(i+1, j)/((ap(i+1, j)+ap(i, j))*0.5))*(grid_p(i+1, j)-grid_p(i, j))
                vs(i, j)=0 !No South face at south nodes
                !==========================================================!
            else if (j == 1 .and. i /= 1 .and. i /= mesh_res_y+1) then !West border without corner points
                !==========================================================!
                ue(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j+1))&
                    +0.5*((dy(i, j+1)/ap(i, j+1))*0.5*(grid_p(i, j+2)-grid_p(i, j))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j)))& !PE-PW to PE-PP
                    -(dy(i, j+1)/((ap(i, j+1)+ap(i, j))*0.5))*(grid_p(i, j+1)-grid_p(i, j))
                if (WB == 1) then
                    uw(i, j)=grid_f_u(i, j) !Outlet/Inlet Face Velocity
                else
                    uw(i, j)=0 !No West face at west nodes
                end if
                vn(i, j)=0 !Vertical BC along West Boundary
                vs(i, j)=0 !Vertical BC along West Boundary
                !==========================================================!
            else if (i == mesh_res_y+1 .and. j /= 1 .and. j /= mesh_res_x+1) then !North border without corner points
                !==========================================================!
                ue(i, j)=Lv !Horizontal BC along North Boundary
                uw(i, j)=Lv !Horizontal BC along North Boundary
                vn(i, j)=0 !No North face at north nodes
                vs(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i-1, j))&
                    +0.5*((dx(i-1, j)/ap(i-1, j))*0.5*(grid_p(i, j)-grid_p(i-2, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i, j)-grid_p(i-1, j)))& !PN-PS to PP-PS
                    -(dx(i-1, j)/((ap(i-1, j)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i-1, j))
                !==========================================================!
            else if (j == mesh_res_x+1 .and. i /= 1 .and. i /= mesh_res_y+1) then !East border without corner points
                !==========================================================!
                if (EB == 1) then
                    ue(i, j)=grid_f_u(i, j) !Outlet/Inlet Face Velocity
                else
                    ue(i, j)=0 !No East face at east nodes
                end if
                uw(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j-1))&
                    +0.5*((dy(i, j-1)/ap(i, j-1))*0.5*(grid_p(i, j)-grid_p(i, j-2))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j)-grid_p(i, j-1)))& !PE-PW to PP-PW
                    -(dy(i, j-1)/((ap(i, j-1)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i, j-1))
                vn(i, j)=0 !Vertical BC along East Boundary
                vs(i, j)=0 !Vertical BC along East Boundary
                !==========================================================!
            else if (i == 2 .and. j > 2 .and. j < mesh_res_x) then !1 row above South border without corners
                !==========================================================!
                ue(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j+1))&
                    +0.5*((dy(i, j+1)/ap(i, j+1))*0.5*(grid_p(i, j+2)-grid_p(i, j))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j+1)/((ap(i, j+1)+ap(i, j))*0.5))*(grid_p(i, j+1)-grid_p(i, j))
                uw(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j-1))&
                    +0.5*((dy(i, j-1)/ap(i, j-1))*0.5*(grid_p(i, j)-grid_p(i, j-2))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j-1)/((ap(i, j-1)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i, j-1))
                vn(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i+1, j))&
                    +0.5*((dx(i+1, j)/ap(i+1, j))*0.5*(grid_p(i+2, j)-grid_p(i, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i+1, j)/((ap(i+1, j)+ap(i, j))*0.5))*(grid_p(i+1, j)-grid_p(i, j))
                vs(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i-1, j))&
                    +0.5*((dx(i-1, j)/ap(i-1, j))*0.5*(grid_p(i, j)-grid_p(i-1, j))& !PP-PSS to PP-PS
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i-1, j)/((ap(i-1, j)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i-1, j))
                !==========================================================!
            else if (j == 2 .and. i > 2 .and. i < mesh_res_y) then !1 column right of West border without corners
                !==========================================================!
                ue(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j+1))&
                    +0.5*((dy(i, j+1)/ap(i, j+1))*0.5*(grid_p(i, j+2)-grid_p(i, j))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j+1)/((ap(i, j+1)+ap(i, j))*0.5))*(grid_p(i, j+1)-grid_p(i, j))
                uw(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j-1))&
                    +0.5*((dy(i, j-1)/ap(i, j-1))*0.5*(grid_p(i, j)-grid_p(i, j-1))& !PP-PWW to PP-PW
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j-1)/((ap(i, j-1)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i, j-1))
                vn(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i+1, j))&
                    +0.5*((dx(i+1, j)/ap(i+1, j))*0.5*(grid_p(i+2, j)-grid_p(i, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i+1, j)/((ap(i+1, j)+ap(i, j))*0.5))*(grid_p(i+1, j)-grid_p(i, j))
                vs(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i-1, j))&
                    +0.5*((dx(i-1, j)/ap(i-1, j))*0.5*(grid_p(i, j)-grid_p(i-2, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i-1, j)/((ap(i-1, j)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i-1, j))
                !==========================================================!
            else if (i == mesh_res_y .and. j > 2 .and. j < mesh_res_x) then !1 row below North border without corners
                !==========================================================!
                ue(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j+1))&
                    +0.5*((dy(i, j+1)/ap(i, j+1))*0.5*(grid_p(i, j+2)-grid_p(i, j))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j+1)/((ap(i, j+1)+ap(i, j))*0.5))*(grid_p(i, j+1)-grid_p(i, j))
                uw(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j-1))&
                    +0.5*((dy(i, j-1)/ap(i, j-1))*0.5*(grid_p(i, j)-grid_p(i, j-2))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j-1)/((ap(i, j-1)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i, j-1))
                vn(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i+1, j))&
                    +0.5*((dx(i+1, j)/ap(i+1, j))*0.5*(grid_p(i+1, j)-grid_p(i, j))& !PNN-PP to PN-PP
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i+1, j)/((ap(i+1, j)+ap(i, j))*0.5))*(grid_p(i+1, j)-grid_p(i, j))
                vs(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i-1, j))&
                    +0.5*((dx(i-1, j)/ap(i-1, j))*0.5*(grid_p(i, j)-grid_p(i-2, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i-1, j)/((ap(i-1, j)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i-1, j))
                !==========================================================!
            else if (j == mesh_res_x .and. i > 2 .and. i < mesh_res_y) then !1 column left of East border without corners
                !==========================================================!
                ue(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j+1))&
                    +0.5*((dy(i, j+1)/ap(i, j+1))*0.5*(grid_p(i, j+1)-grid_p(i, j))& !PEE-PP to PE-PP
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j+1)/((ap(i, j+1)+ap(i, j))*0.5))*(grid_p(i, j+1)-grid_p(i, j))
                uw(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j-1))&
                    +0.5*((dy(i, j-1)/ap(i, j-1))*0.5*(grid_p(i, j)-grid_p(i, j-2))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j-1)/((ap(i, j-1)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i, j-1))
                vn(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i+1, j))&
                    +0.5*((dx(i+1, j)/ap(i+1, j))*0.5*(grid_p(i+2, j)-grid_p(i, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i+1, j)/((ap(i+1, j)+ap(i, j))*0.5))*(grid_p(i+1, j)-grid_p(i, j))
                vs(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i-1, j))&
                    +0.5*((dx(i-1, j)/ap(i-1, j))*0.5*(grid_p(i, j)-grid_p(i-2, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i-1, j)/((ap(i-1, j)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i-1, j))
                !==========================================================!
            else if (i == 2 .and. j == 2) then !Inner South-West corner
                !==========================================================!                
                ue(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j+1))&
                    +0.5*((dy(i, j+1)/ap(i, j+1))*0.5*(grid_p(i, j+2)-grid_p(i, j))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j+1)/((ap(i, j+1)+ap(i, j))*0.5))*(grid_p(i, j+1)-grid_p(i, j))
                uw(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j-1))&
                    +0.5*((dy(i, j-1)/ap(i, j-1))*0.5*(grid_p(i, j)-grid_p(i, j-1))& !PP-PWW to PP-PW
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j-1)/((ap(i, j-1)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i, j-1))
                vn(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i+1, j))&
                    +0.5*((dx(i+1, j)/ap(i+1, j))*0.5*(grid_p(i+2, j)-grid_p(i, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i+1, j)/((ap(i+1, j)+ap(i, j))*0.5))*(grid_p(i+1, j)-grid_p(i, j))
                vs(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i-1, j))&
                    +0.5*((dx(i-1, j)/ap(i-1, j))*0.5*(grid_p(i, j)-grid_p(i-1, j))& !PP-PSS to PP-PS
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i-1, j)/((ap(i-1, j)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i-1, j))
                !==========================================================!
            else if (i == 2 .and. j == mesh_res_x) then !Inner South-East corner
                !==========================================================!
                ue(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j+1))&
                    +0.5*((dy(i, j+1)/ap(i, j+1))*0.5*(grid_p(i, j+1)-grid_p(i, j))& !PEE-PP to PE-PP
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j+1)/((ap(i, j+1)+ap(i, j))*0.5))*(grid_p(i, j+1)-grid_p(i, j))
                uw(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j-1))&
                    +0.5*((dy(i, j-1)/ap(i, j-1))*0.5*(grid_p(i, j)-grid_p(i, j-2))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j-1)/((ap(i, j-1)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i, j-1))
                vn(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i+1, j))&
                    +0.5*((dx(i+1, j)/ap(i+1, j))*0.5*(grid_p(i+2, j)-grid_p(i, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i+1, j)/((ap(i+1, j)+ap(i, j))*0.5))*(grid_p(i+1, j)-grid_p(i, j))
                vs(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i-1, j))&
                    +0.5*((dx(i-1, j)/ap(i-1, j))*0.5*(grid_p(i, j)-grid_p(i-1, j))& !PP-PSS to PP-PS
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i-1, j)/((ap(i-1, j)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i-1, j))
                !==========================================================!
            else if (i == mesh_res_y .and. j == 2) then !Inner North-West corner
                !==========================================================!
                ue(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j+1))&
                    +0.5*((dy(i, j+1)/ap(i, j+1))*0.5*(grid_p(i, j+2)-grid_p(i, j))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j+1)/((ap(i, j+1)+ap(i, j))*0.5))*(grid_p(i, j+1)-grid_p(i, j))
                uw(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j-1))&
                    +0.5*((dy(i, j-1)/ap(i, j-1))*0.5*(grid_p(i, j)-grid_p(i, j-1))& !PP-PWW to PP-PW
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j-1)/((ap(i, j-1)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i, j-1))
                vn(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i+1, j))&
                    +0.5*((dx(i+1, j)/ap(i+1, j))*0.5*(grid_p(i+1, j)-grid_p(i, j))& !PNN-PP to PN-PP
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i+1, j)/((ap(i+1, j)+ap(i, j))*0.5))*(grid_p(i+1, j)-grid_p(i, j))
                vs(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i-1, j))&
                    +0.5*((dx(i-1, j)/ap(i-1, j))*0.5*(grid_p(i, j)-grid_p(i-2, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i-1, j)/((ap(i-1, j)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i-1, j))
                !==========================================================!
            else if (i == mesh_res_y .and. j == mesh_res_x) then !Inner North-East corner
                !==========================================================!
                ue(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j+1))&
                    +0.5*((dy(i, j+1)/ap(i, j+1))*0.5*(grid_p(i, j+1)-grid_p(i, j))& !PEE-PP to PE-PP
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j+1)/((ap(i, j+1)+ap(i, j))*0.5))*(grid_p(i, j+1)-grid_p(i, j))
                uw(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j-1))&
                    +0.5*((dy(i, j-1)/ap(i, j-1))*0.5*(grid_p(i, j)-grid_p(i, j-2))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j-1)/((ap(i, j-1)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i, j-1))
                vn(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i+1, j))&
                    +0.5*((dx(i+1, j)/ap(i+1, j))*0.5*(grid_p(i+1, j)-grid_p(i, j))& !PNN-PP to PN-PP
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i+1, j)/((ap(i+1, j)+ap(i, j))*0.5))*(grid_p(i+1, j)-grid_p(i, j))
                vs(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i-1, j))&
                    +0.5*((dx(i-1, j)/ap(i-1, j))*0.5*(grid_p(i, j)-grid_p(i-2, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i-1, j)/((ap(i-1, j)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i-1, j))
                !==========================================================!
            else !Remaining control volume (Center)
                !==========================================================!
                ue(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j+1))&
                    +0.5*((dy(i, j+1)/ap(i, j+1))*0.5*(grid_p(i, j+2)-grid_p(i, j))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j+1)/((ap(i, j+1)+ap(i, j))*0.5))*(grid_p(i, j+1)-grid_p(i, j))
                uw(i, j)=0.5*(grid_f_u(i, j)+grid_f_u(i, j-1))&
                    +0.5*((dy(i, j-1)/ap(i, j-1))*0.5*(grid_p(i, j)-grid_p(i, j-2))&
                    +(dy(i, j)/ap(i, j))*0.5*(grid_p(i, j+1)-grid_p(i, j-1)))&
                    -(dy(i, j-1)/((ap(i, j-1)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i, j-1))
                vn(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i+1, j))&
                    +0.5*((dx(i+1, j)/ap(i+1, j))*0.5*(grid_p(i+2, j)-grid_p(i, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i+1, j)/((ap(i+1, j)+ap(i, j))*0.5))*(grid_p(i+1, j)-grid_p(i, j))
                vs(i, j)=0.5*(grid_f_v(i, j)+grid_f_v(i-1, j))&
                    +0.5*((dx(i-1, j)/ap(i-1, j))*0.5*(grid_p(i, j)-grid_p(i-2, j))&
                    +(dx(i, j)/ap(i, j))*0.5*(grid_p(i+1, j)-grid_p(i-1, j)))&
                    -(dx(i-1, j)/((ap(i-1, j)+ap(i, j))*0.5))*(grid_p(i, j)-grid_p(i-1, j))
                !==========================================================!
            end if
        end do
    end do
    
end subroutine u_v_Solver

subroutine Tridiagonal_Matrix_Solver(mesh_res_y)
    
    use Dynamic_Array !Contains initialization for allocatable arrays
    
    implicit none
    
    integer :: mesh_res_y, i 
    !Represents the number of linear equations - 1 (mesh_res)
    !i is used to count
    
    real*8 :: M
    !Used to store the denominator of terms A_j and C_jP
    
    !Initializes the modified coefficients (forward elimination)
    A_j(1) = aj(1)/Dj(1)
    C_jP(1) = Cj(1)/Dj(1)
    
    !Evaluates the modified coefficients (forward elimination)
    do i = 2, mesh_res_y+1
        M = Dj(i)-A_j(i-1)*Bj(i)
        A_j(i)=aj(i)/M
        C_jP(i)=(Cj(i)-C_jP(i-1)*Bj(i))/M
    end do
    
    !Initializes the final value of [x] (back substitution)
    X(mesh_res_y+1)= C_jP(mesh_res_y+1)
    
    !Evaluates the values of [x] (back substitution)
    do i = mesh_res_y, 1, -1
        X(i)=C_jP(i)-A_j(i)*X(i+1)
    end do
    
end subroutine Tridiagonal_Matrix_Solver

subroutine Output_Results(mesh_res_x, mesh_res_y)
    
    use Dynamic_Array !Contains initialization for allocatable array
    
    implicit none
    
    integer :: mesh_res_x, mesh_res_y, i, j
    !Assigns variable for mesh size and array row and column variables
    
    open(2222, file="Output.txt")
    !Opens or creates the output file 
    
    !Writes each element of the final u velocity array into the output file
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            !Writes everything for the row on the same line
            write (2222, "(f7.4, ', ')", advance="no") (grid_f_u(i, j))
        end do
        !Creates a new line for the next row
        write (2222, *) " "
    end do
    
    !Writes each element of the final v velocity array into the output file
    write (2222, *) " "
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            !Writes everything for the row on the same line
            write (2222, "(f7.4, ', ')", advance="no") (grid_f_v(i, j))
        end do
        !Creates a new line for the next row
        write (2222, *) " "
    end do
    
    !Writes each element of the final pressure correction array into the output file
    write (2222, *) " "
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            !Writes everything for the row on the same line
            write (2222, "(f9.3, ', ')", advance="no") (grid_f_pc(i, j))
        end do
        !Creates a new line for the next row
        write (2222, *) " "
    end do
    
    !Writes each element of the final pressure array into the output file
    write (2222, *) " "
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            !Writes everything for the row on the same line
            write (2222, "(f9.3, ', ')", advance="no") (grid_p(i, j))
        end do
        !Creates a new line for the next row
        write (2222, *) " "
    end do
    
    !Writes each dx into the output file
    write (2222, *) " "
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            !Writes everything for the row on the same line
            write (2222, "(f7.4, ', ')", advance="no") (dx(i, j))
        end do
        !Creates a new line for the next row
        write (2222, *) " "
    end do
    
    !Writes each dy into the output file
    write (2222, *) " "
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            !Writes everything for the row on the same line
            write (2222, "(f7.4, ', ')", advance="no") (dy(i, j))
        end do
        !Creates a new line for the next row
        write (2222, *) " "
    end do
    
    !Writes the x values in the coordinate grid
    write (2222, *) " "
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            write (2222, "(f10.5, ', ')", advance="no") (coordinate_grid(i, j, 1))
        end do
        write (2222, *) " "
    end do
    
    !Writes the y values in the coordinate grid
    write (2222, *) " "
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            write (2222, "(f10.5, ', ')", advance="no") (coordinate_grid(i, j, 2))
        end do
        write (2222, *) " "
    end do
    
    !Writes the dx values in the distance grid
    write (2222, *) " "
    do i = 1, mesh_res_y
        do j = 1, mesh_res_x
            write (2222, "(f10.5, ', ')", advance="no") (distance_grid(i, j, 1))
        end do
        write (2222, *) " "
    end do
    
    !Writes the dy values in the distance grid
    write (2222, *) " "
    do i = 1, mesh_res_y
        do j = 1, mesh_res_x
            write (2222, "(f10.5, ', ')", advance="no") (distance_grid(i, j, 2))
        end do
        write (2222, *) " "
    end do
    
    !Writes the x coordinates across horizontal line, v velocities
    write (2222, *) " "
    i = (mesh_res_y/2)+1
        do j = 1, mesh_res_x+1
            write (2222, "(f10.5, ', ')", advance="no") (coordinate_grid(i, j, 1))
            write (2222, "(f10.5, ', ')", advance="no") (grid_f_v(i, j))
            write (2222, *) " "
        end do
    
    !Writes the y coordinates across vertical line, u velocities
    write (2222, *) " "
    j = (mesh_res_x/2)+1
        do i = 1, mesh_res_y+1
            write (2222, "(f10.5, ', ')", advance="no") coordinate_grid(i, j, 2)
            write (2222, "(f10.5, ', ')", advance="no") grid_f_u(i, j)
            write (2222, *) " "
        end do
    
    close(2222)
    !Closes the opened file  
    
    open (2223, file="Tecplot_Output.dat")
    !Opens or creates output file for tecplot
    
    !Writes header information for tecplot to use coordinate and velocity data
    write (2223, *) "Title = Ahmad's Cavity Flow Solver"
    write (2223, *) "Variables = x, y, u, v"
    write (2223, 7) mesh_res_x+1, mesh_res_y+1
    7 format("Zone, I=", I5," J=", I4," F=point") !Formats previous line's output
    do i = 1, mesh_res_y+1
        do j = 1, mesh_res_x+1
            write (2223, 8) coordinate_grid(i, j, 1),coordinate_grid(i, j, 2), grid_f_u(i, j), grid_f_v(i, j)
            !Writes out the x coordinate, y coordinate, u velocity and v velocity
            8 format(1x, 4(f19.8,1x)) !Formats previous line's output
        end do
    end do
    
    close (2223)
    !Closes the opened file
    
end subroutine Output_Results