module backup
!
!  SUBROUTINES
!
!  1) save_main_fields
!  2) load_backup
!  3) save_backup
!
implicit none

contains
!------------------------------------------------------------------------------------------
   subroutine save_main_fields(nx, ny, u, v, T, p, ue, un, ve, vn, sim_id)
      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) :: u   ! Cartesian velocity
      real(8), dimension (nx*ny), intent(in) :: v   ! Cartesian velocity
      real(8), dimension (nx*ny), intent(in) :: T   ! Temperature
      real(8), dimension (nx*ny), intent(in) :: p   ! Pressure at center o 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) :: un  ! Cartesian velocity u at center of north face
      real(8), dimension (nx*ny), intent(in) :: ve  ! Cartesian velocity v at center of east face
      real(8), dimension (nx*ny), intent(in) :: vn  ! Cartesian velocity v at center of north face
      character*100,  intent(in) :: sim_id ! Simulation identification
      


      ! Auxiliary variables

      integer :: i, j, np, nps, npw

      real(8), dimension(nx) :: Tbs, Tbn, pbs, pbn

      real(8), dimension(ny) :: Tbw, Tbe, pbw, pbe



      ! p and T on the south boundary

      j = 2

      do i = 2, nx-1

         np  = nx * (j-1) + i

         nps  = np - nx

         Tbs(i) = ( T(np) + T(nps) ) / 2.d0

         pbs(i) = ( p(np) + p(nps) ) / 2.d0

      end do

      ! p and T on the north boundary

      j = ny

      do i = 2, nx-1

         np  = nx * (j-1) + i

         nps  = np - nx

         Tbn(i) = ( T(np) + T(nps) ) / 2.d0

         pbn(i) = ( p(np) + p(nps) ) / 2.d0

      end do

      ! p and T on the west boundary

      i = 2

      do j = 2, ny-1

         np  = nx * (j-1) + i

         npw  = np - 1

         Tbw(j) = ( T(np) + T(npw) ) / 2.d0

         pbw(j) = ( p(np) + p(npw) ) / 2.d0

      end do

      ! p and T on the east boundary

      i = nx

      do j = 2, ny-1

         np  = nx * (j-1) + i

         npw  = np - 1

         Tbe(j) = ( T(np) + T(npw) ) / 2.d0

         pbe(j) = ( p(np) + p(npw) ) / 2.d0

      end do




      ! Saving fields to a file for further comparisons

      open(15,file='./mach2d-output/mainfields_' // trim(adjustl(sim_id)) // '.dat')

      write(15,*) u
      write(15,*) v
      write(15,*) T
      write(15,*) p
      write(15,*) ue
      write(15,*) un
      write(15,*) ve
      write(15,*) vn
      write(15,*) Tbs
      write(15,*) Tbn
      write(15,*) Tbw
      write(15,*) Tbe
      write(15,*) pbs
      write(15,*) pbn
      write(15,*) pbw
      write(15,*) pbe

      close(15)


   end subroutine save_main_fields
!------------------------------------------------------------------------------------------
   subroutine load_backup( nx, ny, ito, norm1, tcpu, p, pin, pl, T, ro, roe, ron &
         ,                     u, v, ue, ve, un, vn, Uce, Vcn, sim_id) ! 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
      character*100,  intent(in) :: sim_id ! Simulation identification
      
      ! Auxiliary variables

      integer :: nx_s, ny_s, i, j, aux_i, np
      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,*) nx_s,ny_s
         read(10,*) ito, norm1, tcpu  
         ! Carregar a pressão na entrada
         do j = 1, ny_s
            read(10,*) pin(j)
         end do  
         ! CARREGAR TODOS OS OUTROS CAMPOS
         do j=1, ny_s
            aux_i = nx_s * (j-1)
            do i=1, nx_s
              np = aux_i + i
              read(10,*) p(np),pl(np),T(np),ro(np),roe(np),ron(np),u(np),v(np),ue(np),&
                         ve(np),un(np),vn(np),Uce(np),Vcn(np)
            enddo
         end do  

         close(10)

         ito = ito + 1

      else

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

      end if

   end subroutine load_backup
!------------------------------------------------------------------------------------------
   subroutine save_backup( nx, ny, it, norm1, tcpu, p, pin, pl, T, ro, roe, ron &
         ,                     u, v, ue, ve, un, vn, Uce, Vcn, sim_id) ! 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
      character*100,  intent(in) :: sim_id ! Simulation identification

        integer :: np, i, j, aux_i
        open(10, file = './mach2d-output/mach2d-backup-'// trim(adjustl(sim_id)))

        write(10,*) nx, ny
        write(10,*) it, norm1, tcpu
 
         ! Salvar a pressão na entrada
        do j = 1, ny
            write(10,'(1PE22.15)') pin(j)
        end do
 
        ! SALVAR TODOS OS OUTROS CAMPOS
         do j=1, ny
            aux_i = nx * (j-1)
            do i=1,nx
                np = aux_i + i
                write(10,'(14(1PE23.15))') p(np),pl(np),T(np),ro(np),roe(np),ron(np),u(np),v(np),ue(np),&
                            ve(np),un(np),vn(np),Uce(np),Vcn(np)
            enddo
        end do      
     
        close(10)

   end subroutine save_backup
!------------------------------------------------------------------------------------------
end module backup
