!>
!!    \mainpage Mach2D-5.8.2.1
!!
!!    \section Introduction
!!    This version of the code Mach2D solves the Euler equations for a bidimensional
!!    planar or axisymmetric flow over the forebody of an arbitrary smooth body.
!!


program main

   use data
   use grid
   use user
   use coefficients
   use solvers
   use postp

   implicit none

   ! Reads parameters, allocates variables and sets vectors to zero
   call initialize_variables

   !> \brief Generates the north and south boundaries of the grid based on a power-law model.
   call get_grid_boundary_power_law( nx, ny, akn, aks, la, lb, lr, rb, x, y) ! Last two are output

   ! Generates the grid according to kg option
   call set_grid(kg, nx, ny, a1, x, y) ! Last 2 are inoutput


   ! Calculates the centroids of each real volume
   call get_real_centroids_xy( kcm, nx, ny, x, y, xp, yp ) ! Output: last two


   ! Calculates the components of the metric tensor and the Jacobian
   call get_metrics(nx, ny, x, y, xp, yp             & ! Input
   ,          xe, ye, xen, yen, xk, yk, xke, yke, Jp & ! Output
   ,          Je, Jn, alphae, gamman, betae, betan )   ! Output


   ! All the radii are equal 1 if planar flow is choose. If axisymmetric flow, radii are calculated
   call get_radius(coord, nx, ny, y, yp, re, rn, rp    ) ! Output: last three


   ! Guess the initial conditions
   call get_initial_conditions( nx, ny, PF, TF, Rg, GF, MF, UF &
      , xe, ye, xk, yk, alphae, betae, betan, gamman &
      , p, pl, T, ro, roe, ron, u, v, ue, ve, un, vn, Uce, Vcn, Vcw ) ! UF and last 15 are output



   ! Calculations of the thermal properties
   call set_cp( nx, ny, cp) ! Last is output


   ! Beginning time evolution cycle

   ! Openning file of residuals
   open(rid, file = './mach2d_output/mach2d_' // trim(adjustl(sim_id)) // '_residual.dat' )

   ! Initialising CPU time

   call cpu_time(tcpu1)

   do it = 1, itmax

      ! Changing the time step
      if ( it <= It1 ) dt = dt1
      if ( it >  It1 .and. it < It2 ) dt = dt1 + (dt2-dt1)*(it-It1)/(It2-It1)
      if ( it >= It2 ) dt = dt2

      ! Changing the beta value
      if ( it <= itb1 ) beta = beta1
      if ( it >  itb1 .and. it < itb2 ) beta = beta1 + (beta2-beta1)*(it-itb1)/(itb2-itb1)
      if ( it >= itb2 ) beta = beta2

      ! Fields update
      roa = ro
      Ta  = T
      pa  = p
      ua  = u
      va  = v
      uea = ue
      vea = ve
      una = un
      vna = vn

      ! Starting the iterative cycle for the instant t

      do iti = 1, itimax

         ! Updating variables

         Ucea = Uce
         Vcna = Vcn


         ! Calculations of the thermal properties
         call set_cp( nx, ny, cp) ! Last is output


         ! Coefficients of the linear system for u (real volumes)
         call get_u_coefficients( nx, ny, dt, rp, re, rn, Jp &
            ,                        roe, ron, roa, Uce, Vcn, au ) ! Last one is output



         ! Source of the linear system for u (real volumes)
         call get_u_source( nx, ny, beta, dt, rp, re, rn, ye, yk &
            ,             Jp, roe, ron, roa, p, Uce, Vcn, ua, u, cup, bu ) ! The last two are output


         ! Coefficients and source of the linear system of u for the fictitious volumes
         call set_bcu( nx, ny, UF, xk, yk, alphae, betae, u, v, Uce, Vcw, au, bu)


         ! Coefficients of the linear system for v (real volumes)
         call get_v_coefficients( nx, ny, dt, rp, re, rn, Jp &
            ,                        roe, ron, roa, Uce, Vcn, av ) ! Last one is output


         ! Source of the linear system for v (real volumes)
         call get_v_source( nx, ny, beta, dt, rp, re, rn, xe, xk &
            ,          Jp, roe, ron, roa, p, Uce, Vcn, va, v, cvp, bv ) ! Last two are output


         ! Coefficients and source of the linear system of v (fictitious volumes)
         call set_bcv( nx, ny, xk, yk, u, v, Uce, Vcw, av, bv)


         ! Calculates SIMPLEC coefficients at the internal real faces
         call get_internal_simplec_coefficients( nx, ny, re, rn, xe, ye, xk, yk, au, av & ! Input
         ,                                                    due, dve, dun, dvn, de, dn ) ! Output


         ! Calculates the SIMPLEC coefficients at the boundary faces
         call get_boundary_simplec_coefficients( nx, ny, xe, ye &
            ,  due, dve, dun, dvn, de, dn) ! InOutput: last six


         ! Calculates g

         g = 1.d0/(Rg*T)


         ! Calculates the coefficients of the linear system for pressure correction
         ! g, ro, Uce and Vcn used in this subroutine are those calculated in the previous iteraction
         ! de and dn must be calculated with the coef. of the linear sytem for u and v from which
         ! u* and v* are obtained.
         call get_p_coefficients(nx, ny, dt, rp, re, rn, Jp, Uce, Vcn, roe, ron, g, de, dn, ap) ! Output: last one

         ! Solves the linear system for u

         ! LU decomposition
         call lu2d9(au,dl9,du9,nxy,nx,ny)

         ! Linear system solution
         call fb2d9(au,dl9,du9,ny,nx,nxy,u,tol_u,nitm_u,bu)


         ! Norm initialization
         norm = 0.d0

         ! Calculates the norm of the residual of the linear system for u
         call norm_l1_9d( nx, ny, u, bu, au, norm)

         ! Relative residual calculation
         !res_u = norm / sum(abs(bu))
         res_u = norm / norm_l1_b( nx, ny, bu)

         ! Solves the linear system for v

         ! LU decomposition
         call lu2d9(av,dl9,du9,nxy,nx,ny)

         ! Linear system solution
         call fb2d9(av,dl9,du9,ny,nx,nxy,v,tol_u,nitm_u,bv)


         ! Norm initialization
         norm = 0.d0

         ! Calculates the norm of the residual of the linear system for v
         call norm_l1_9d( nx, ny, v, bv, av, norm)

         ! Relative residual calculation
         !res_v = norm / sum(abs(bv))
         res_v = norm / norm_l1_b( nx, ny, bv)

         ! Uses velocities at nodes to calculate velocities at internal real faces
         call get_velocities_at_internal_faces(nx, ny, dt, rp, re, rn, xe, ye, xk, yk & ! Input
         ,                             xen, yen, xke, yke, Jp, cup, cvp                & ! Input
         ,                             au, av, roa, p, u, v, uea, vea, una, vna        & ! Input
         ,                             ue, ve, un, vn, Uce, Vcn)                         ! Output


         ! Calculates the velocities at boundary faces based on the boundary conditions
         call get_velocities_at_boundary_faces( nx, ny, UF, xe, ye, xk, yk, u, v, ue, ve, un, vn, Uce, Vcn)
         ! Last six are output


         ! ro, Uce and Vcn are the incorrect ones (obtained with p*);
         ! rom, Ucem and Vcnm are those of the previous iteraction
         call get_p_source(nx, ny, dt, rp, re, rn, Jp, roe, ron, roe, ron, ro, roa &
            , Ucea, Uce, Vcna, Vcn, bp) ! Output: last one


         do itm = 1, imax

            !Coefficients of the linear system of pl for the fictitious volumes
            call set_bcp( nx, ny, alphae, betae, betan, gamman, Uce, Vcw, p, ap, bp) ! Output: last two

            ! Solves the linear system for pl

            ! LU decomposition
            call lu2d5(ap,dl5,du5,nxy,nx,ny)

            ! Solution of the linear system
            call fb2d5(ap,dl5,du5,ny,nx,nxy,pl,tol_p,nitm_p,bp)

         end do

         ! Norm initialization
         norm = 0.d0

         ! Calculates the norm of the residuals
         call norm_l1_5d( nx, ny, pl, bp, ap, norm)

         ! Residual calculation
         res_p = norm

         ! Pressure correction  (SIMPLEC method)
         call get_pressure_density_correction_with_pl( nx, ny, pl, g, ro, p ) ! InOutput: last two


         ! Extrapolation of p to fictitious

         call get_p_extrapolation_to_fictitious(nx, ny, PF, alphae, betae, betan, gamman, Uce, Vcw, p)


         ! Density at nodes with the state equation
         call get_density_at_nodes( nx, ny, Rg, p, T, ro) ! ro is output


         ! Velocity correction (SIMPLEC method)
         call get_u_v_at_real_nodes_with_pl( nx, ny, xe, ye, xk & ! Input
         ,                                    yk, rp, pl, au, av & ! Input
         ,                                    u, v )               ! Input and Output


         ! Extrapolation is based on the boundary conditions
         call get_u_v_extrapolation_to_fictitious( nx, ny, UF, xk, yk, alphae, betae, Uce, Vcw, u, v)
         ! Last two are inout


         ! Velocity correction at faces (SIMPLEC method)
         call get_velocity_correction_at_faces_with_pl( nx, ny       & ! Input
         ,                           due, dve, dun, dvn, de, dn, pl  & ! Input
         ,                           ue, ve, un, vn, Uce, Vcn )        ! InOutput


         ! Calculates density at faces using the corrected density and velocities at nodes
         call get_density_at_faces( nx, ny, beta, ro, Uce, Vcn, roe, ron) ! roe and ron are output


         ! Coefficients and source of the temperature linear system
         call get_T_coefficients_and_source( nx, ny, beta, dt, rp, re, rn  &
            ,                                xe, ye, xk, yk, Jp, roe, ron  &
            ,                roa, p, pa, cp, Uce, Vcn, u, v, T, Ta, at, bt ) ! Last 2 are output


         ! Calculating the contravariant velocity Vcw on the east boundary
         call get_Vcw(nx, ny, xe, xk, ue, ve, Vcw)


         ! Coefficients of the linear system of T for the fictitious volumes
         call set_bcT( nx, ny, TF, Uce, Vcw, T, alphae, betae, betan, gamman, at, bt) ! Output: last two entries


         ! Solves the linear system for T

         ! LU decomposition
         call lu2d9(at,dl9,du9,nxy,nx,ny)

         ! Linear system solution
         call fb2d9(at,dl9,du9,ny,nx,nxy,T,tol_u,nitm_u,bt)

         ! Norm initialization
         norm = 0.d0

         ! Calculates the norm of the residual of the linear system for T
         call norm_l1_9d( nx, ny, T, bt, at, norm)

         ! Relative residual calculation
         !res_t = norm / sum(abs(bt))
         res_t = norm / norm_l1_b( nx, ny, bt)

         ! Density at nodes with the state equation
         call get_density_at_nodes( nx, ny, Rg, p, T, ro) ! ro is output


         ! Calculates density at faces using the corrected density and velocities at nodes
         call get_density_at_faces( nx, ny, beta, ro, Uce, Vcn, roe, ron) ! roe and ron are output

         res = res_u + res_v + res_T + res_p

         ! Checking for error

         if ( isnan(res) ) then

            write(*,*) "NaN found. Stopping..."

            write(lid,*) "NaN found. Stopping..."

            stop

         end if

      end do

      ! Periodic data print

      if ( it == 1 ) then

         call get_cdfi(nx, ny, coord, Rg, PF, TF, UF, rb, yk, rn, p, Cdfi)

         write(  *,"(A10, 2(1X, A23), A10)") 'it', 'Residuals', 'Cdfi', 'it_stop'

         write(rid,"(A10, 2(1X, A23), A10)") 'it', 'Residuals', 'Cdfi', 'it_stop'

         write(rid,"(I10, 2(1X,ES23.16), I10)") it, res, Cdfi, it_stop

      end if


      if ( mod( it, wlf ) == 0 ) then

         call get_cdfi(nx, ny, coord, Rg, PF, TF, UF, rb, yk, rn, p, Cdfi)

         write(*,"(I10, 2(1X,ES23.16), I10)") it, res, Cdfi, it_stop

         write(rid,"(I10, 2(1X,ES23.16), I10)") it, res, Cdfi, it_stop

      end if

      ! Updating residual mean

      res_mean = res_mean + res

      ! Checking for convergence

      if ( mod( it, nit_res ) == 0 ) then

         res_mean = res_mean / dble(nit_res)

         if ( ( res_mean < tol_res ) .and. ( res_check == 0 ) ) then

            it_stop = 2 * it

            res_check = 1

         end if

         if ( it == it_stop ) exit

         res_mean = 0.d0

      end if


   end do

   ! Finishing cpu time

   call cpu_time(tcpu2)

   tcpu = tcpu2 - tcpu1

   ! Closing file of residuals
   close(rid)


   ! Calculating the RAM memory used
   RAM = get_RAM(sim_id)


   ! Calculating gamma from cp
   call set_gamma( nx, ny, Rg, cp, gcp) ! Last is output

   ! Post processing data
   call post_processing(nx, ny, it, sem_a, sem_g, w_g, w_cam, lid &
      ,  sim_id, tcpu, RAM, res, lr, rb, x, y, xe, ye, xk, yk, Jp, Rg, cp &
      ,  gcp, Uce, Vcn, de, dn, pl, bp, xp, yp, u, v, p, T, ro, Cdfi)


end program

