program main
  !
  !============================================================
  !                  MODULE DEPENDENCY
  !     -------------------------------------------------
  !
  !     |-<- module solvers -<- modules mtdma2d5 and mtdma2d9
  !     |
  !     |-<- module user -<- module coefficients
  !     |    |
  !     |    |->- gridboundary.dat
  !     |
  !     |-<- module grid -<- gridboundary.dat
  !     |
  !     |-<- module data -<- input_file_parameters
  !     |
  !     |-<- module blomax
  !     |
  !     |-<- module post-processing
  !     |    |
  !     |    |->- output data
  !     |
  !     |->- main program
  !
  !============================================================
  !
  !
  !
  !$ use omp_lib ! OpenMP library for parallel computation
  !
  use data
  use grid
  use solvers
  use user
  use blomax
  use postp

  implicit none

  call get_numerical_solution

contains

  subroutine get_numerical_solution
    implicit none
    integer :: it, im


    ! Read input data from file
    ! Allocate and Initialize variables
    call initialize_variables


    ! Generates the boundary nodes of the domain. Also gives ig and Sg.
    call get_boundary_nodes( nx, ny, ig, Sg, rcg) ! Output: last three entries


    ! Read boundary grid from a file and generates the grid
    call set_grid(kg, nx, ny, a1, x, y)


    ! Centroids of control volumes
    call get_real_centroids_xy( 1, nx, ny, x, y, xp, yp ) ! Output: last two entries


    ! Metrics calculation
    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


    ! Radius
    call get_radius(coord, nx, ny, y, yp, radius, re, rn, rp ) ! Output: last four entries


    ! Wall temperature
    if ( ccTw == 1 ) call set_wall_temperature( nx, Twall )

    roe = 0.d0
    ron = 0.d0

    ! Flow initialization (1D isentropic flow)
    call  get_initial_guess( nx, ny, ig, modvis, beta, po, T0, gamma  & ! Input
         ,                   Rg, Sg, ye, yk, radius, rn, x, xp        & ! Input
         ,                   M1D, p1D, T1D, u1D, p, T, u, v, ue       & ! Output
         ,                   un, Uce, Vcn, uin, vin, pin, Tin, Mw     & ! Output
         ,                   fm1D, Fd1D, Fpv1D, de, dn, ro, roe, ron  & ! Output
         ,                   au, ap, b ) ! Output

    av = au
    at = au

    ! Baldwin-Lomax module initialization

    if ( modtur == 1 ) then

        call initialize_blomax(nx, ny, x, y)

    end if

    if ( reload == 1 ) then
        call load_backup( nx, ny, ito, norm1, tcpuo, p, pin, pl, T, ro, roe, ron &
        ,                     u, v, ue, ve, un, vn, Uce, Vcn) ! Output: all, except nx and ny
    end if

    !
    ! rid = residual file id
    !
    if ( reload == 1 ) then
        open(rid, file = "./mach2d-output/mach2d-" // trim(adjustl(sim_id)) // "-residual.dat", POSITION='APPEND')
    else
        open(rid, file = "./mach2d-output/mach2d-" // trim(adjustl(sim_id)) // "-residual.dat")
    end if

    ! Starts cpu time
    call cpu_time(tcpu1)

    do it = ito, itmax

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

       ! Selecting beta
       if ( it <= itb1 ) beta = beta1
       if ( it >  itb1 .and. it < itb2 ) beta = beta1 + (beta2-beta1)*dble(it-itb1)/dble(itb2-itb1)
       if ( it >= itb2 ) beta = beta2

       norm = 0.d0

       ! Calculating mean variables used in the eddy viscosity calculation

        if ( modtur == 1 ) then

            ro_av  = 0.5d0 * ( ro + roa )
            ron_av = 0.5d0 * ( ron + rona )
            u_av   = 0.5d0 * ( u + ua )
            v_av   = 0.5d0 * ( v + va )

        end if

       ! Updating variables for the new time step

       ua  = u
       va  = v
       uea = ue
       vea = ve
       una = un
       vna = vn
       pa  = p
       Ta  = T
       roa = ro
       rona= ron

       Ucem = Uce
       Vcnm = Vcn

       pina = pin

       ! Updating conditions at the nozzle entrance
       call get_uin_vin_pin_Tin_Mw( nx, ny, gamma, Rg, po, T0, u & ! Input
            ,                       uin, vin, pin, Tin, Mw)        ! Output

       call get_plin_and_p_fictitious( nx, ny, pina, pin, plin, p) ! Output: last two entries


       ! Calculation of the gas properties
       call set_cp_and_gamma( nx, ny, Rg, cp, gcp) ! Output: last two entries

       call set_laminar_viscosity_at_nodes( nx, ny, vlp) ! Output: last entry

       ! CAUTION: In order to apply correctly this function, laminar viscosity vlp must be known
       ! in real and fictitious volumes (except CV of corners SW, SE, NW and NE)
       call get_laminar_viscosity_at_faces( nx, ny, vlp, vle, vln) ! Output: last two entries

        if ( modtur == 1 ) then

            ! Calculates the dimensional eddy viscosity at center of real volumes
            call get_eddy_viscosity_blomax(nx, ny, x, y, xp, yp, xk, yk, Jp, Jn, vln &
            ,                         ro_av, ron_av, u_av, v_av, vtp) ! Last one is output

            ! Extrapolation of the eddy viscosity to the fictitious volumes
            call get_eddy_viscosity_fictitious(nx, ny, vtp)

            vlp = vlp + vtp

            ! CAUTION: In order to apply correctly this function, laminar viscosity vlp must be known
            ! in real and fictitious volumes (except CV of corners SW, SE, NW and NE)
            call get_laminar_viscosity_at_faces( nx, ny, vlp, vle, vln) ! Output: last two entries

        end if

       call set_thermal_conductivity_at_nodes( nx, ny, kp) ! Output: last entry

        if ( modtur == 1 ) then

            kp = kp + vtp * cp / 0.9d0

        end if

       ! CAUTION: In order to apply correctly this function, thermal conductivity kp must be known
       ! in real and fictitious volumes (except CV of corners SW, SE, NW and NE)
       call get_thermal_conductivity_at_faces( nx, ny, kp, ke, kn) ! Output: last two entries


       !======================= BEGIN PARALLEL COMPUTATION ===================================

       !$OMP PARALLEL DEFAULT(SHARED)

       !$OMP SECTIONS

       !$OMP SECTION
       ! Calculating u coefficients and source
       call get_u_coefficients( nx, ny, modvis, dt, rp, re, rn, Jp, Je, Jn &
            ,                   ye, yk, alphae, betae, betan, gamman       &
            ,                   vle, vln, roe, ron, roa, Uce, Vcn, au )    ! Output: last entry

       call get_u_source( nx, ny, modvis, beta, dt, rp, re, rn    &
            ,              xe, ye, xk, yk, xke, yke, xen, yen     &
            ,              Jp, Je, Jn, roe, ron, roa, p, vle, vln &
            ,              Uce, Vcn, ua, u, v, cup, sup, bu )      ! Output: last three entries

       call set_bcu( nx, ny, modvis, x, xp, u, au, bu) ! Output: last two entries

       !$OMP SECTION
       ! Calculating v coefficients and source
       call get_v_coefficients( nx, ny, coord, modvis, dt, rp, re, rn, Jp, Je, Jn &
            ,                          xe, xk, alphae, betae, betan, gamman       &
            ,                          vle, vln, vlp, roe, ron, roa, Uce, Vcn, av ) ! Output: last entry

       call get_v_source( nx, ny, coord, modvis, beta, dt, rp, re, rn   &
            ,                    xe, ye, xk, yk, xke, yke, xen, yen     &
            ,                    Jp, Je, Jn, roe, ron, roa, p, vle, vln &
            ,                         Uce, Vcn, va, u, v, cvp, svp, b ) ! Output: last three entries

       call set_bcv( nx, ny, modvis, x, xp, v, av, b) ! Output: last two entries

       !$OMP END SECTIONS
       !$OMP END PARALLEL

       !======================== END PARALLEL COMPUTATION ====================================

        ! Calculating SIMPLEC coefficients

        call get_internal_simplec_coefficients( nx, ny, re, rn, xe, ye, xk, yk, au, av & ! Input
            ,                                              due, dve, dun, dvn, de, dn ) ! Output

        call get_boundary_simplec_coefficients( nx, ny, de, dn) ! InOutput: de, dn

        ! Calculating 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 &
            ,                       ro, 1.d0/(Rg*T), de, dn, ap) ! Output: last one


        ! Solving the linear system for u

       call tdma2d9(nx, ny, nitm_u, au, bu, u)

       call norm_l1_9d( nx, ny, u, bu, au, norm)

       ! Solving the linear system for v

       call tdma2d9(nx, ny, nitm_u, av, b, v)

       call norm_l1_9d( nx, ny, v, b, av, norm)


       ! Calculating incorrect velocities at faces (Uce* and Vcn*)

       call get_velocities_at_faces(nx, ny, dt, rp, re, rn, xe, ye, xk, yk     & ! Input
            ,                       xen, yen, xke, yke, Jp, cup, cvp, sup, svp & ! Input
            ,                       au, av, roa, p, u, v, uea, vea, una, vna   & ! Input
            ,                       ue, ve, un, vn, Uce, Vcn)                    ! Output

       call get_Uce_Vcn_at_boundary_faces( nx, ny, ye, u, Uce, Vcn) ! Output: Uce, Vcn

       ! Calculating source of the equation for pressure correction

       ! 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, beta, dt, rp, re, rn, Jp, ro, ro, roa, Ucem, Uce, Vcnm, Vcn, bp) ! Output: last one

       ! Mass correction cycle

       do im = 1, imax

          call set_bcp( nx, ny, x, xp, plin, pl, ap, bp) ! Output: last two entries

          call tdma2d5(nx, ny, nitm_p, ap, bp, pl)

       end do

       call norm_l1_5d( nx, ny, pl, bp, ap, norm)

       call get_pressure_density_correction_with_pl( nx, ny, pl, 1.d0/(Rg*T), ro, p ) ! InOutput: last two entries

       call get_u_v_at_real_nodes_with_pl( nx, ny, xe, ye, xk, yk, rp, pl, au, av, u, v ) ! InOutput: last two entries

       call get_u_v_at_fictitious_nodes_with_pl( nx, ny, modvis, x, xp, u, v) ! InOutput: last two entries

       call get_velocities_at_internal_faces_with_pl( nx, ny, xe, ye, xk, yk   & ! Input
            ,                                         due, dve, dun, dvn, pl   & ! Input
            ,                                         ue, ve, un, vn, Uce, Vcn ) ! InOutput

       call get_Uce_Vcn_at_boundary_faces_with_pl( nx, ny, pl, de, Uce) ! InOutput: Last entry

       call get_density_at_faces( nx, ny, beta, ro, Uce, Vcn, roe, ron) ! Output: last two entries

       ! Solving the linear system for temperature

       call get_T_coefficients_and_source( nx, ny, coord, modvis, beta, dt, rp, re, rn       &
            ,            xe, ye, xk, yk, alphae, betae, betan, gamman, Jp, Je, Jn, roe, ron  &
            ,            roa, p, pa, cp, vlp, ke, kn, Uce, Vcn, u, v, T, Ta, at, b ) ! Output: last two entries

       call set_bcT( nx, ny, ccTw, x, xp, Tin, Twall, T, at, b) ! Output: last two entries

       call tdma2d9(nx, ny, nitm_u, at, b, T)

       call norm_l1_9d( nx, ny, T, b, at, norm)

       call get_density_at_nodes( nx, ny, Rg, p, T, ro) ! ro is output

       call get_density_at_faces( nx, ny, beta, ro, Uce, Vcn, roe, ron) ! roe and ron are output


        ! Saving backup data
        if ( mod(it,wbkp) == 0 ) then

            ! Finishes cpu time
            call cpu_time(tcpu2)

            call save_backup( nx, ny, it, norm1, tcpuo + tcpu2 - tcpu1, p, pin, pl, T, ro, roe, ron &
                ,                     u, v, ue, ve, un, vn, Uce, Vcn) ! Output: none

            !write(*,*)  "Backup"!, tcpu, tcpu1, tcpu2, tcpu2 - tcpu1, tcpu + tcpu2 - tcpu1

        end if

       ! Printing results to listing file
       if ( it == ito ) then

          if ( it == 1 ) norm1 = norm

          ! lid - listing file id
          write(lid,*)
          write(lid,"(A7,4(1X,A23))") "it", "norm / norm1", "fmi (kg/s)", "fme/fmi", "Fd (N)"
          write(lid,*)

          write(*,*)
          write(*,"(A7,4(1X,A23))") "it", "norm / norm1", "fmi (kg/s)", "fme/fmi", "Fd (N)"
          write(*,*)

          ! rid = residual file id
          write(rid,"('#', A4,(1X,A14))") "it", "norm / norm1"

       end if

        ! Printing results to listing file
       if ( mod(it,wlf) == 0 ) then

          call get_mass_flow_rate_and_thrust( nx, ny, re, roe, u, Uce, fmi, fme, Fd) ! Output: last three entries

          write(lid,"(I7,4(1X,ES23.15))") it, norm / norm1, fmi, fme/fmi, Fd
          write(*  ,"(I7,4(1X,ES23.15))") it, norm / norm1, fmi, fme/fmi, Fd
          write(rid,"(I5,1X,ES14.7)") it, norm / norm1
          call flush(lid)
          call flush(rid)
       end if

    end do

    close(rid)

    ! Finishes cpu time
    call cpu_time(tcpu2)

    tcpu = tcpuo + tcpu2 - tcpu1

    call save_backup( nx, ny, itmax, norm1, tcpu, p, pin, pl, T, ro, roe, ron &
        ,                     u, v, ue, ve, un, vn, Uce, Vcn) ! Output: none

    call set_laminar_viscosity_at_nodes( nx, ny, vlp) ! Output: last entry
    call set_thermal_conductivity_at_nodes( nx, ny, kp) ! Output: last entry

    call get_mass_flow_rate_and_thrust( nx, ny, re, roe, u, Uce, fmi, fme, Fd) ! Output: last three entries

    write(lid,*)
    write(lid,*) "MAIN RESULTS"
    write(lid,*)
    write(lid,"(ES23.16,' .....: ',A)") fmi, "fmi  - Mass flow rate in the entrance (kg/s) "
    write(lid,"(ES23.16,' .....: ',A)") fme, "fme  - Mass flow rate in the exit (kg/s)"
    write(lid,"(ES23.16,' .....: ',A)") Fd,  "Fd   - Dynamic thrust (N) "
    write(lid,"(ES23.16,' .....: ',A)") tcpu,"tcpu - CPU time (s) "

    call post_processing(nx, ny, sem_a, sem_g, w_g, w_cam, ccTw, ig, lid, modtur &
       , sim_id, x, y, rp, re, xe, ye, xk, yk, Jp, Sg, rcg, gamma, pr, go, po, Rg &
       , tcpuo, tcpu1, tcpu2, tcpu, Fpv1D, Fd1D, fm1D, fme, cp, vlp, vtp, kp, gcp &
       , Twall, u1D, p1D, T1D, M1D, uin, vin, pin, Tin &
       , roe, Uce, Vcn, de, dn, pl, bp, xp, yp, u, v, p, T, ro)

  end subroutine get_numerical_solution

  subroutine get_mass_flow_rate_and_thrust( nx, ny, re, roe, u, Uce, fmi, fme, Fd)
    implicit none
    integer, intent(in) :: nx  ! Number of volumes in csi direction (real+fictitious)
    integer, intent(in) :: ny  ! Number of volumes in eta direction (real+fictitious)
    real(8), dimension (nx*ny), intent(in) :: re   ! Radius of the center of east face of volume P
    real(8), dimension (nx*ny), intent(in) :: roe  ! Absolute density at east face
    real(8), dimension (nx*ny), intent(in) :: u    ! Cartesian velocity of the last iteraction
    real(8), dimension (nx*ny), intent(in) :: Uce  ! Contravariant velocity U at east face
    real(8), intent(out) :: fmi ! Mass flow rate at entrance
    real(8), intent(out) :: fme ! Mass flow rate at exit
    real(8), intent(out) :: Fd  ! Thrust

    real(8), parameter :: pi = acos(-1.d0)

    integer :: i, j, np

    fmi = 0.d0

    i = 1

    do j = 2, ny-1

       np = nx * (j-1) + i

       fmi = fmi + re(np) * roe(np) * Uce(np)

    end do

    fmi = fmi * 2.d0 * pi



    fme = 0.d0

    Fd = 0.d0

    i = nx-1

    do j = 2, ny-1

       np = nx * (j-1) + i

       fme = fme + re(np) * roe(np) * Uce(np)

       Fd = Fd + re(np) * roe(np) * Uce(np) * ( u(np) + u(np+1) ) / 2.d0

    end do

    fme = fme * 2.d0 * pi

    Fd = Fd * 2.d0 * pi

  end subroutine get_mass_flow_rate_and_thrust

  subroutine save_backup( nx, ny, it, norm1, tcpu, p, pin, pl, T, ro, roe, ron &
        ,                     u, v, ue, ve, un, vn, Uce, Vcn) ! Output: none
        implicit none
        integer, intent(in) :: nx  ! Number of volumes in csi direction (real+fictitious)
        integer, intent(in) :: ny  ! Number of volumes in eta direction (real+fictitious)
        integer, intent(in) :: it  ! Iteraction number when backup was made
        real(8), intent(in) :: norm1 ! Norm of the linear systems of the first iteraction
        real(8), intent(in) :: tcpu  ! Elapsed cpu time
        real(8), dimension (nx*ny), intent(in) :: p   ! Pressure at center of volume P
        real(8), dimension (ny),    intent(in) :: pin ! Pressure in the entrance
        real(8), dimension (nx*ny), intent(in) :: pl  ! Pressure correction at center of volume P
        real(8), dimension (nx*ny), intent(in) :: T   ! Temperature at center of volume P
        real(8), dimension (nx*ny), intent(in) :: ro  ! Specific mass at center of volume P
        real(8), dimension (nx*ny), intent(in) :: roe ! Absolute density at east face
        real(8), dimension (nx*ny), intent(in) :: ron ! Absolute density at north face
        real(8), dimension (nx*ny), intent(in) :: u   ! Cartesian velocity at center of volume P
        real(8), dimension (nx*ny), intent(in) :: v   ! Cartesian velocity at center of volume P
        real(8), dimension (nx*ny), intent(in) :: ue  ! Cartesian velocity u at center of east face
        real(8), dimension (nx*ny), intent(in) :: ve  ! Cartesian velocity v at center of east face
        real(8), dimension (nx*ny), intent(in) :: un  ! Cartesian velocity u at center of north face
        real(8), dimension (nx*ny), intent(in) :: vn  ! Cartesian velocity v at center of north face
        real(8), dimension (nx*ny), intent(in) :: Uce ! Contravariant velocity U at east face
        real(8), dimension (nx*ny), intent(in) :: Vcn ! Contravariant velocity V at north face

        open(10, file = './mach2d-output/mach2d-backup-'// trim(adjustl(sim_id)))

        write(10,*) it, norm1, tcpu, p, pin, pl, T, ro, roe, ron, u, v, ue, ve, un, vn, Uce, Vcn

        close(10)

    end subroutine save_backup

    subroutine load_backup( nx, ny, ito, norm1, tcpu, p, pin, pl, T, ro, roe, ron &
        ,                     u, v, ue, ve, un, vn, Uce, Vcn) ! Output: all, except nx and ny
        implicit none
        integer, intent(in) :: nx  ! Number of volumes in csi direction (real+fictitious)
        integer, intent(in) :: ny  ! Number of volumes in eta direction (real+fictitious)
        integer, intent(out) :: ito ! Iteraction number when backup was made
        real(8), intent(out) :: norm1 ! Norm of the linear systems of the first iteraction
        real(8), intent(out) :: tcpu  ! Elapsed cpu time
        real(8), dimension (nx*ny), intent(out) :: p   ! Pressure at center of volume P
        real(8), dimension (ny),    intent(out) :: pin ! Pressure in the entrance
        real(8), dimension (nx*ny), intent(out) :: pl  ! Pressure correction at center of volume P
        real(8), dimension (nx*ny), intent(out) :: T   ! Temperature at center of volume P
        real(8), dimension (nx*ny), intent(out) :: ro  ! Specific mass at center of volume P
        real(8), dimension (nx*ny), intent(out) :: roe ! Absolute density at east face
        real(8), dimension (nx*ny), intent(out) :: ron ! Absolute density at north face
        real(8), dimension (nx*ny), intent(out) :: u   ! Cartesian velocity at center of volume P
        real(8), dimension (nx*ny), intent(out) :: v   ! Cartesian velocity at center of volume P
        real(8), dimension (nx*ny), intent(out) :: ue  ! Cartesian velocity u at center of east face
        real(8), dimension (nx*ny), intent(out) :: ve  ! Cartesian velocity v at center of east face
        real(8), dimension (nx*ny), intent(out) :: un  ! Cartesian velocity u at center of north face
        real(8), dimension (nx*ny), intent(out) :: vn  ! Cartesian velocity v at center of north face
        real(8), dimension (nx*ny), intent(out) :: Uce ! Contravariant velocity U at east face
        real(8), dimension (nx*ny), intent(out) :: Vcn ! Contravariant velocity V at north face

        ! Auxiliary variables

        logical :: lexist

        inquire(file = './mach2d-output/mach2d-backup-'// trim(adjustl(sim_id)), exist = lexist)

        if ( lexist ) then

            open(10, file = './mach2d-output/mach2d-backup-'// trim(adjustl(sim_id)))

            read(10,*) ito, norm1, tcpu, p, pin, pl, T, ro, roe, ron, u, v, ue, ve, un, vn, Uce, Vcn

            close(10)

            ito = ito + 1

        else

            write(*,*) 'load_backup: Backup file does not exist.'

        end if

    end subroutine load_backup
end program
