module resultados

! objetivo: calcular a solução numérica e 
!           apresentar gráficos e listagens dos resultados

!-------------------------------------------------

use coeficientes
use solvers

implicit none

contains

!-------------------------------------------------

  subroutine solucao_numerica
  

    real*8 :: hm ! tamanho médio dos volumes de controle

    open(8,file=trim(adjustl(caso))//'.itn')

    write(8,1) 
    1 format( t1,'Iteração', t12,'fmr (kg/s)',  t36,'Fsr (N)', &
             t60,'uc (m/s)', t84,'vc (m/s)',   t108,'pc (Pa)'  )

    open(11,file='r.dat')

    write(10,3) 
    3 format(//,t3,'Iteração',t14,'M?dia da norma L1 dos resíduos / norma na itera??o 1',/)

    call cpu_time ( tcpu1 )

    do it = ito, itmax    ! ciclo iterativo externo

       norma = 0.0d0

       ! resolve u

       call coeficientes_velocidade
       call fonte_velocidade_u
       call ccu
       select case ( solver )
          case (1) 
             call jacobi (nx,ny,nxy,nitm_u,a,b,u)
          case (2) 
             call gauss_seidel (nx,ny,nxy,nitm_u,a,b,u)
          case (3) 
             call sor (nx,ny,nxy,nitm_u,w,a,b,u)
          case (4) 
             call adi (nx,ny,nxy,nitm_u,a,b,u)
          case (5)
             call lu2d5(a,dl,du,nxy,nx,ny)
             call fb2d5(a,dl,du,ny,nx,nxy,u,tolu,nitm_u,b)
       end select
       call norma_L1 (nx,ny,nxy,a,b,u,norma)

       ! resolve v

       call fonte_velocidade_v
       call ccv
       select case ( solver )
          case (1) 
             call jacobi (nx,ny,nxy,nitm_u,a,b,v)
          case (2) 
             call gauss_seidel (nx,ny,nxy,nitm_u,a,b,v)
          case (3) 
             call sor (nx,ny,nxy,nitm_u,w,a,b,v)
          case (4) 
             call adi (nx,ny,nxy,nitm_u,a,b,v)
          case (5)
             call lu2d5(a,dl,du,nxy,nx,ny)
             call fb2d5(a,dl,du,ny,nx,nxy,v,tolu,nitm_u,b)
       end select
       call norma_L1 (nx,ny,nxy,a,b,v,norma)

       ! calcula ue, vn

       call media_coef_x_vel ! NHT/94

       do im = 1, imax     ! ciclo iterativo da massa

          ! resolve pl

          call coeficientes_pressao
          call ccp
          select case ( solver )
             case (1) 
                call jacobi (nx,ny,nxy,nitm_p,a,b,pl)
             case (2) 
                call gauss_seidel (nx,ny,nxy,nitm_p,a,b,pl)
             case (3) 
                call sor (nx,ny,nxy,nitm_p,w,a,b,pl)
             case (4) 
                call adi (nx,ny,nxy,nitm_p,a,b,pl)
             case (5)
                call lu2d5(a,dl,du,nxy,nx,ny)
                call fb2d5(a,dl,du,ny,nx,nxy,pl,tolp,nitm_p,b)
          end select
       
          ! correção de p, u, v, ue, vn com pl

          call atualiza_pressao
          call atualiza_velocidades_nodais
          call atualiza_velocidades_medias

       end do

       call norma_L1 (nx,ny,nxy,a,b,pl,norma)

       if ( it == 1 ) norma_1 = norma

       if ( it==1 .or. it==itmax .or. norma/norma_1<=tole .or. mod(it,w_r)==0 ) then
          write(10,4) it, norma/norma_1
          write(11,4) it, norma/norma_1
          4 format(2x,i8,1pe15.6)
          call globais
          write(8,2) it, fmr, Fsr, uc, vc, pc
          2 format(i8,5(1pe24.15))
       end if

       if ( norma/norma_1 <= tole ) exit

       ! atualização para novo incremento de tempo

       ua  = u
       va  = v
       uea = ue
       vna = vn

    end do

    call cpu_time ( tcpu2 )

    tcpu = tcpu2-tcpu1 + tcpuo

    close(8)

    open(99,file=trim(adjustl(caso))//'.bin',form='unformatted')
    write(99) itmax,norma_1,tcpu,u,v,p,pl,ue,vn
    close(99)

    write(10,6) norma_1
    6 format(/,5x,'Média da norma L1 dos res?duos na iteração 1 =',1pe15.6)

    close(11)

    call curva_residuo

    deallocate ( a, ua, va, uea, vna, cup, cvp )

    call perfil_u

    call perfil_v

    call chain_retangulo

    call perfil_chain_retangulo

    call chain_trapezio
    
    call extremos
    
    call esc_psi_p2_2d
    
    call esc_min_u_p2_2d
    
    call esc_min_v_p2_2d
    
    call esc_max_v_p2_2d
    
    if (nxy >= 25) then
       call esc_psi_p4_2d
       call esc_min_u_p4_2d
       call esc_min_v_p4_2d
       call esc_max_v_p4_2d
    end if
    
    if (nxy >= 49) then
       call esc_psi_p6_2d
       call esc_min_u_p6_2d
       call esc_min_v_p6_2d
       call esc_max_v_p6_2d
    end if
    
    
    if ( w_cam == 1 ) then

       write(10,10)
       10 format(//,10x,'*** campos ***', //, &
                    10x,'na horizontal: contorno oeste para leste', //, &
                    10x,'na vertical  : contorno sul   para norte'  )

       write(10,11)
       11 format(//,30x,'velocidades u (nodais)')
       call print1(u,ny,nx,nxy)

       write(10,12)
       12 format(//,30x,'velocidades v (nodais)')
       call print1(v,ny,nx,nxy)

       write(10,13)
       13 format(//,30x,'velocidades ue (face leste)')
       call print1(ue,ny,nx,nxy)

       write(10,14)
       14 format(//,30x,'velocidades vn (face norte)')
       call print1(vn,ny,nx,nxy)

       write(10,15)
       15 format(//,30x,'press?es (nodais)')
       call print1(p,ny,nx,nxy)

       write(10,16)
       16 format(//,30x,'varia??es de pressão (nodais)')
       call print1(pl,ny,nx,nxy)

       write(10,17)
       17 format(//,30x,'fontes de plinha (nodais)')
       call print1(b,ny,nx,nxy)

       write(10,18)
       18 format(//,30x,'função de corrente (integral retângulo)')
       call print1(psi_r,ny,nx,nxy)

       write(10,19)
       19 format(//,30x,'função de corrente (integral trapézio)')
       call print1(psi_t,ny,nx,nxy)

    end if

    hm = dsqrt(1.0d0/((nx-2)*(ny-2)))

    write(10,20) hm, uc, umin, y_umin, vc, vmin, x_vmin, vmax, x_vmax, &
                 pc, fmr, fmt, dudy, Fsr, Fnr, Fnr2, Fnt, &
                 psi_r_min, xmin_psi_r, ymin_psi_r, &
                 psi_t_min, xmin_psi_t, ymin_psi_t

    20 format(//,5x,'h (m) =',1pe24.15,//, &
                 5x,'uc (m/s) =',1pe24.15,//, &
                 5x,'umin (m/s) =',1pe24.15,/, &
                 5x,'y_umin (m) =',1pe24.15,//, &
                 5x,'vc (m/s) =',1pe24.15,//, &
                 5x,'vmin (m/s) =',1pe24.15,/, &
                 5x,'x_vmin (m) =',1pe24.15,//, &
                 5x,'vmax (m/s) =',1pe24.15,/, &
                 5x,'x_vmax (m) =',1pe24.15,//, &
                 5x,'pc (Pa) =',1pe24.15,//, &
                 5x,'fmr (kg/s) =',1pe24.15,/, &
                 5x,'fmt (kg/s) =',1pe24.15,//, &
                 5x,'dudy (1/s) =',1pe24.15,//, &
                 5x,'Fsr (N) =',1pe24.15,//, &
                 5x,'Fnr (N)  =',1pe24.15,/, &
                 5x,'Fnr2 (N) =',1pe24.15,/, &
                 5x,'Fnt (N)  =',1pe24.15,//, &
                 5x,'Psi_r_min (m2/s) =',1pe24.15,/, &
                 5x,'xmin_psi_r (m) =',1pe24.15,/,&
                 5x,'ymin_psi_r (m) =',1pe24.15,//,&
                 5x,'Psi_t_min (m2/s) =',1pe24.15,/, &
                 5x,'xmin_psi_t (m) =',1pe24.15,/,&
                 5x,'ymin_psi_t (m) =',1pe24.15)

    open(77,file=trim(adjustl(caso))//'.Richardson_3p0')

    write(77,78) hm, uc, umin, y_umin, vc, vmin, x_vmin, vmax, x_vmax, &
                 pc, fmr, fmt, dudy, Fsr, Fnr, Fnr2, Fnt, &
                 psi_r_min, xmin_psi_r, ymin_psi_r, &
                 psi_t_min, xmin_psi_t, ymin_psi_t

    78 format(  2x, 1pe24.15, 2x, 'h (m)',/, &
              ' 1', 1pe24.15, 2x, 'uc (m/s)',/, &
              ' 2', 1pe24.15, 2x, 'umin (m/s)',/, &
              ' 3', 1pe24.15, 2x, 'y_umin (m)',/, &
              ' 4', 1pe24.15, 2x, 'vc (m/s)',/, &
              ' 5', 1pe24.15, 2x, 'vmin (m/s)',/, &
              ' 6', 1pe24.15, 2x, 'x_vmin (m)',/, &
              ' 7', 1pe24.15, 2x, 'vmax (m/s)',/, &
              ' 8', 1pe24.15, 2x, 'x_vmax (m)',/, &
              ' 9', 1pe24.15, 2x, 'pc (Pa)',/, &
              '10', 1pe24.15, 2x, 'fmr (kg/s)',/, &
              '11', 1pe24.15, 2x, 'fmt (kg/s)',/, &
              '12', 1pe24.15, 2x, 'dudy (1/s)',/, &
              '13', 1pe24.15, 2x, 'Fsr (N)',/, &
              '14', 1pe24.15, 2x, 'Fnr (N)',/, &
              '15', 1pe24.15, 2x, 'Fnr2 (N)',/, &
              '16', 1pe24.15, 2x, 'Fnt (N)',/, &
              '17', 1pe24.15, 2x, 'Psi_r_min (m2/s)',/, &
              '18', 1pe24.15, 2x, 'xmin_psi_r (m)',/,&
              '19', 1pe24.15, 2x, 'ymin_psi_r (m)',/,&
              '20', 1pe24.15, 2x, 'Psi_t_min (m2/s)',/, &
              '21', 1pe24.15, 2x, 'xmin_psi_t (m)',/,&
              '22', 1pe24.15, 2x, 'ymin_psi_t (m)')

    call perfil_u_interesse

    call perfil_v_interesse

    write(77,780) psi(no_psi_min), x_psi_min, y_psi_min, u_min_2d, x_umin_2d, y_umin_2d, &
                 v_min_2d, x_vmin_2d, y_vmin_2d, v_max_2d, x_vmax_2d, y_vmax_2d 
    
   780 format('53', 1pe24.15, 2x, 'Psi_t_min (m2/s) (2D)',/, &
              '54', 1pe24.15, 2x, 'xmin_psi_t (m) (2D)',/, &
              '55', 1pe24.15, 2x, 'ymin_psi_t (m) (2D)',/, &
              '56', 1pe24.15, 2x, 'umin (m/s) (2D)',/, &
              '57', 1pe24.15, 2x, 'x_umin (m)(2D)',/, &
              '58', 1pe24.15, 2x, 'y_umin (m)(2D)',/, &
              '59', 1pe24.15, 2x, 'vmin (m/s) (2D)',/, &
              '60', 1pe24.15, 2x, 'x_vmin (m) (2D)',/, &
              '61', 1pe24.15, 2x, 'y_vmin (m) (2D)',/, &
              '62', 1pe24.15, 2x, 'vmax (m/s) (2D)',/, &
              '63', 1pe24.15, 2x, 'x_vmax (m) (2D)',/, &
              '64', 1pe24.15, 2x, 'y_vmax (m) (2D)')
    
    
    close(77)

    write(10,31) tcpu
    31 format(//,f14.2,' = tempo de CPU acumulado (s)')

    note_caso = 'notepad '// trim(adjustl(caso)) // '.itn'
    dos = system(note_caso) ! lista arquivo de itera??es

    note_caso = 'notepad '// trim(adjustl(caso)) // '.Richardson_3p0'
    dos = system(note_caso) ! lista arquivo das vari?veis de interesse

  end subroutine solucao_numerica

!-------------------------------------------------

  subroutine norma_L1 (nx,ny,nxy,a,b,x,norma)

    integer,intent(in) :: nx, ny, nxy

    real*8,intent(in),dimension(nxy,5) :: a ! coeficientes
    real*8,intent(in),dimension(nxy)   :: b ! fonte
    real*8,intent(in),dimension(nxy)   :: x ! solução

    real*8,intent(inout) :: norma

    np = 1 ! canto inferior-esquerdo 
    norma = norma + dabs(x(np)*a(np,3) + a(np,4)*x(np+1) &
          + a(np,5)*x(np+nx) - b(np))

    ! contorno inferior
    j = 1
    do i = 2, nx
       np = (j-1)*nx + i
       norma = norma + dabs(x(np)*a(np,3) &
             + a(np,2)*x(np-1) + a(np,4)*x(np+1) &
             + a(np,5)*x(np+nx) - b(np))
    end do

    ! pontos internos e contornos esquerdo e direito
    do i = 1, nx
       do j = 2, ny-1
          np = (j-1)*nx + i
          norma = norma + dabs(x(np)*a(np,3) + a(np,1)*x(np-nx) &
                + a(np,2)*x(np-1) + a(np,4)*x(np+1) &
                + a(np,5)*x(np+nx) - b(np))
       end do
     end do

    ! contorno superior
    j = ny
    do i = 1, nx-1
       np = (j-1)*nx + i
       norma = norma + dabs(x(np)*a(np,3) + a(np,1)*x(np-nx) &
             + a(np,2)*x(np-1) + a(np,4)*x(np+1) - b(np))
    end do

    ! canto superior-direito
    np = nxy
    norma = norma + dabs(x(np)*a(np,3) + a(np,1)*x(np-nx) &
          + a(np,2)*x(np-1) - b(np))

  end subroutine norma_L1

!-------------------------------------------------

  subroutine globais

    real*8 :: x, y, fp, glp, uexe, uexeo ! auxiliares

    i = nx/2
    j = ny/2
    np = (j-1)*nx + i

    uc = 5.0d-1*(ue(np)+ue(np+nx))

    vc = 5.0d-1*(vn(np)+vn(np+1))

    pc = 2.5d-1*(p(np)+p(np+1)+p(np+nx)+p(np+nx+1))

    ! Fnr2 = força viscosa na direções x da tampa (contorno norte) sobre o
             ! fluido com UDS-2, uP e integral retângulo 
    j = ny-1
    Fnr2 = 0.0d0
    do i = 2, nx-1
       np = (j-1)*nx + i
       if ( tipo == 1 ) then
          Fnr2 = Fnr2 + (8*uex(np+nx) + u(np-nx) - 9*u(np))
       else
          Fnr2 = Fnr2 + (8*Ut + u(np-nx) - 9*u(np))
       end if
    end do
    Fnr2 = Fnr2 * mi * dx / ( 3 * dy )

    ! Fnr = força viscosa na direção x da tampa (contorno norte) sobre o
            ! fluido com UDS, uP e integral ret?ngulo 
    j = ny-1
    Fnr = 0.0d0
    do i = 2, nx-1
       np = (j-1)*nx + i
       if ( tipo == 1 ) then
          Fnr = Fnr + uex(np+nx) - u(np)
       else
          Fnr = Fnr + Ut - u(np)
       end if
    end do
    Fnr = 2.0d0 * Fnr * mi * dx / dy

    ! Fnt = força viscosa na direção x da tampa (contorno norte) sobre o
            ! fluido com UDS, ue e integral trap?zio 
    j = ny-1
    Fnt = 0.0d0
    i = 2
    np = (j-1)*nx + i
    if ( tipo == 1 ) then
       x = (i-1)*dx
       y = 1.0d0
       fp = (x**4) - 2.0d0*(x**3) + (x**2)
       glp = 4.0d0*(y**3) - 2.0d0*y
       uexe = 8.0d0*fp*glp
       Fnt = Fnt + uexe - ue(np)
       do i = 3, nx-2
          np = (j-1)*nx + i
          uexeo = uexe
          x = (i-1)*dx
          fp = (x**4) - 2.0d0*(x**3) + (x**2)
          uexe = 8.0d0*fp*glp
          Fnt = Fnt + uexeo - ue(np-1) + uexe - ue(np)
       end do
       i = Nx-2
       np = (j-1)*nx + i
       Fnt = Fnt + uexe - ue(np)
    else
       Fnt = Fnt + Ut - ue(np)
       do i = 3, nx-2
          np = (j-1)*nx + i
          Fnt = Fnt + 2*Ut - ue(np-1) - ue(np)
       end do
       i = Nx-2
       np = (j-1)*nx + i
       Fnt = Fnt + Ut - ue(np)
    end if
    Fnt = Fnt * dx * mi / dy

    ! Fsr = força viscosa na direçao x da parede do contorno sul sobre o
            ! fluido com DDS, uP e integral retângulo 
    j = 2
    Fsr = 0.0d0
    do i = 2, nx-1
       np = (j-1)*nx + i
       Fsr = Fsr + u(np)
    end do
    Fsr = - 2 * Fsr * mi * dx / dy

    ! derivada du/dy com ue e UDS em x=0.5 e y=1 (1/s); válido para nx par 
    i  = nx/2
    j  = ny-1
    np = (j-1)*nx + i
    dudy = 2.0d0 * (1.0d0 - ue(np)) / dy

    ! fmr = fluxo de massa em y=1/2 e x[0;1/2] com vn e integral retângulo
    j = ny/2
    fmr = 0.0d0
    do i = 2, nx/2
       np = (j-1)*nx + i
       fmr = fmr + vn(np)
    end do
    fmr = fmr * dx

    ! fmt = fluxo de massa em y=1/2 e x[0;1/2] com vn e integral trapézio
    j = ny/2
    fmt = 0.0d0
    i = 2
    np = (j-1)*nx + i
    fmt = fmt + vn(np) / 2
    do i = 2, nx/2-1
       np = (j-1)*nx + i
       fmt = fmt + vn(np) + vn(np+1)
    end do
    i = Nx/2
    np = (j-1)*nx + i
    fmt = fmt + ( vn(np) + vc ) / 2
    fmt = fmt * dx / 2

  end subroutine globais

!-------------------------------------------------

  subroutine print1 ( b, nj, ni, nij )

    integer,intent(in) :: nij, ni, nj
    real*8,intent(in),dimension(nij) :: b ! incógnita a imprimir

    integer :: k, ll, ipart, npi, nps, l

    k=0
    ll=8
    ipart = 1

    2 write(10,11)ipart

    do 1 j=1,nj
       np=(j-1)*ni + k
       npi=np+1
       nps=np+ll
       if(ni.lt.8) nps=np+ni
       write(10,10)(b(l),l=npi,nps)
    1 continue

    if(ni.lt.8) goto 3
    if(nps.eq.nij) goto 3

    k=k+8
    ll=ni-k
    if(ll.gt.8)ll=8
    ipart = ipart + 1

    goto 2

    10 format(8(1pe15.6))
    11 format(/,2x,'parte ',i3)
  
    3 continue

  end subroutine print1

! -----------------------------------------------

  subroutine curva_residuo

    integer :: k ! auxiliar
    real*8 :: uj, y

    open(7,file='Stokes_1p6_r.gnu')

    do k = 1, 7
       read(7,*)
    end do
    write(7,3) head
    3 format("set title '",a62,/, &
             "replot")
    close(7)

    dos = system('wgnuplot Stokes_1p6_r.gnu')

  end subroutine curva_residuo

! -----------------------------------------------

  subroutine perfil_u

    integer :: k ! auxiliar
    real*8 :: uj, y
    allocate ( vet_u_min_1d(ny-2) ) ! vetor contendo o perfil de interesse


    open(11,file='u.dat')

    write(10,1)
    1 format(//,5x,'velocidade ue em x=1/2 (todos os nós)',//, &
                t4,'ue', t19,'y')

    umin = Ut
    i = nx/2
    do j = 1, ny
       np = (j-1)*nx + i
       
       if ( j == 1 ) then
          uj = 0.0d0
          y  = 0.0d0
       else
          uj = ue(np)
          y  = (j-2)*dy + dy/2.0d0
          if ( uj < umin ) then
             umin = uj
             y_umin = y
          end if
       end if
       if ( j == ny) then
          uj = Ut
          y  = H
       end if
       write(10,2) uj, y
       write(11,2) uj, y
       if ( ( j > 1 ) .and. ( j < ny ) ) then  
             vet_u_min_1d(j-1)= uj !salva perfil em vetor
       end if
       
    end do

    close(11)

    if ( tipo == 1 ) then
       open(7,file='Stokes_1p6_u_Shih.gnu')
       do k = 1, 8
          read(7,*)
       end do
    else
       open(7,file='Stokes_1p6_u_Ghia.gnu')
       do k = 1, 10
          read(7,*)
       end do
    end if
    write(7,3) head
    3 format("set title '",a62,/, &
             "replot")

    close(7)

    if ( tipo == 1 ) then
       dos = system('wgnuplot Stokes_1p6_u_Shih.gnu')
    else
       dos = system('wgnuplot Stokes_1p6_u_Ghia.gnu')
    end if

    2 format(2(1pe15.6))
    
    ! Salva dados para entrada em interp_1
    
    itp = minval([minloc(vet_u_min_1d)-1,(ny-2)-minloc(vet_u_min_1d)])
    
    1301 format(A<len(trim(adjustl(caso)))>)
    
    write(caso1,1301) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso1)),"_umin_1d",".txt"
   
    itp1 = itp ! minloc(vet_u_min_1d)-1  !indica quantos pontos serão coletados abaixo e acima de u_min_1d
    
    itp2 =  minloc(vet_u_min_1d) - itp1  !localiza o primeiro ponto a ser coletado (abaixo e u_min_1d)
    
    open(1302,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(1302,1303) arquivo_caso
    
    1303 format (3x, A40, 19x, 'nome do caso') 
        
    write (1302, 1304) 2*itp1(1)+1
        
    1304 format (1x, i6,48x, '= N, número de nós centrados em u_min_1d')
   				            
	y = 0.0q0 !inicializa a coordenada y do perfil
	
	do j = itp2(1), itp2(1)+2*itp1(1) 
	   
	   y =  (j-1)*dy + dy/2.0d0
	  
	   write (1302, 1305) y, vet_u_min_1d(j),  j
	    
	   1305 format (1pe48.32, 1pe48.32,2x, '(y, u(1/2,y)) do nó ', 1x, i3)
				   
	end do
    
    close(1302)

  end subroutine perfil_u

! -----------------------------------------------

  subroutine perfil_v

    integer :: k ! auxiliar
    real*8 :: vi, x
    allocate ( vet_v_min_1d(nx-2), vet_v_max_1d(nx-2)  ) ! vetores contendo os perfis de interesse

    open(11,file='v.dat')
    1 format(//,5x,'velocidade vn em y=1/2 (todos os nós)',//, &
                t4,'x', t19,'vn')

    write(10,1)

    vmin = Ut
    vmax = 0.0d0
    j = ny/2
    do i = 1, nx
       np = (j-1)*nx + i
       
       if ( i == 1 ) then
          vi = 0.0d0
          x  = 0.0d0
       else
          vi = vn(np)
          x  = (i-2)*dx + dx/2.0d0
          if ( vi < vmin ) then
             vmin = vi
             x_vmin = x
          end if
          if ( vi > vmax ) then
             vmax = vi
             x_vmax = x
          end if
       end if
       if ( i == nx) then
          vi = 0.0d0
          x  = L
       end if
       write(10,2) x, vi
       write(11,2) x, vi
       
       
       if ( ( i > 1 ) .and. ( i < nx ) ) then  
             vet_v_min_1d(i-1)= vi !salva perfil em vetor
             vet_v_max_1d(i-1)= vi !salva perfil em vetor
       end if
              
    
    end do

    close(11)

    if ( tipo == 1 ) then
       open(7,file='Stokes_1p6_v_Shih.gnu')
       do k = 1, 7
          read(7,*)
       end do
    else
       open(7,file='Stokes_1p6_v_Ghia.gnu')
       do k = 1, 9
          read(7,*)
       end do
    end if
    write(7,3) head
    3 format("set title '",a62,/, &
             "replot")

    close(7)

    if ( tipo == 1 ) then
       dos = system('wgnuplot Stokes_1p6_v_Shih.gnu')
    else
       dos = system('wgnuplot Stokes_1p6_v_Ghia.gnu')
    end if

    2 format(2(1pe15.6))
    
    ! Salva dados para entrada em interp_1 (min v) 
    
    itp = minval([minloc(vet_v_min_1d)-1,(nx-2)-minloc(vet_v_min_1d)])
    
    1301 format(A<len(trim(adjustl(caso)))>)
    
    write(caso1,1301) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso1)),"_vmin_1d",".txt"
   
    itp1 = itp ! indica quantos pontos serão coletados abaixo e acima de v_min_1d
    
    itp2 =  minloc(vet_v_min_1d) - itp1  !localiza o primeiro ponto a ser coletado (abaixo e v_min_1d)
    
    open(1302,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(1302,1303) arquivo_caso
    
    1303 format (3x, A40, 19x, 'nome do caso') 
        
    write (1302, 1304) 2*itp1(1)+1 
        
    1304 format (1x, i6,48x, '= N, número de nós centrados em v_min_1d')
   				            
	x = 0.0q0 !inicializa a coordenada y do perfil
	
	do i = itp2(1), itp2(1)+2*itp1(1)

	      x  = (i-1)*dx + dx/2.0d0
	  
	      write (1302, 1305) x, vet_v_min_1d(i),  i
	    
		1305 format (1pe48.32, 1pe48.32,2x, '(v, v(x,1/2)) do nó ', 1x, i3)
				   
	end do
    
    close(1302)
    
    ! Salva dados para entrada em interp_1 (max v) 
    
    itp = minval([maxloc(vet_v_max_1d)-1,(nx-2)-maxloc(vet_v_max_1d)])
    
    1306 format(A<len(trim(adjustl(caso)))>)
    
    write(caso1,1306) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso1)),"_vmax_1d",".txt"
   
    itp1 = itp ! indica quantos pontos serão coletados abaixo e acima de v_min_1d
    
    itp2 =  maxloc(vet_v_max_1d) - itp1  !localiza o primeiro ponto a ser coletado (abaixo e v_min_1d)
    
    open(1307,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(1307,1308) arquivo_caso
    
    1308 format (3x, A40, 19x, 'nome do caso') 
        
    write (1307, 1309) 2*itp1(1)+1 
        
    1309 format (1x, i6,48x, '= N, número de nós centrados em v_max_1d')
   				            
	x = 0.0q0 !inicializa a coordenada y do perfil
	
	do i = itp2(1), itp2(1)+2*itp1(1)

	      x  = (i-1)*dx + dx/2.0d0
	  
	      write (1307, 1310) x, vet_v_max_1d(i),  i
	    
		1310 format (1pe48.32, 1pe48.32,2x, '(v, v(x,1/2)) do nó ', 1x, i3)
				   
	end do
    
    close(1307)

  end subroutine perfil_v

! -----------------------------------------------

  subroutine perfil_u_interesse
    
    
    integer :: w_int, k ! auxiliares
    real*8  :: uj, y

    if ( base > (ny-2) ) return

    write(10,1)
    1 format(//,5x,'m?dia de ue em x=1/2 (s? pontos de interesse)',//, &
                t4,'u', t28,'y')

    w_int = (ny-2) / base

    k = 22
    i = nx/2
    do j = 2, ny - w_int
       np = (j-1)*nx + i
       if ( mod((j-1),w_int)==0 .and. j/=ny-1 ) then
          uj = ( ue(np) + ue(np+nx) ) / 2
          y  = (j-1)*dy
          k  = k + 1
          write(10,2) uj, y
          write(77,3) k, uj, y
       end if
    end do

    2 format(2(1pe24.15))

    3 format( i2, 1pe24.15, 2x, 'u em x=1/2 e y=', 1pe11.5 )

  end subroutine perfil_u_interesse

! -----------------------------------------------

  subroutine perfil_v_interesse

    integer :: w_int, k ! auxiliares
    real*8  :: vi, x

    if ( base > (nx-2) ) return

    write(10,1)
    1 format(//,5x,'média de vn em y=1/2 (s? pontos de interesse)',//, &
                t4,'x', t28,'v')

    w_int = (nx-2) / base

    k = 37
    j = ny/2
    do i = 2, nx - w_int
       np = (j-1)*nx + i
       if ( mod((i-1),w_int)==0 .and. i/=nx-1 ) then
          vi = ( vn(np) + vn(np+1) ) / 2
          x  = (i-1)*dx
          k  = k + 1
          write(10,2) x, vi
          write(77,3) k, vi, x
       end if
    end do

    2 format(2(1pe24.15))

    3 format( i2, 1pe24.15, 2x, 'v em y=1/2 e x=', 1pe11.5 )

  end subroutine perfil_v_interesse

! -----------------------------------------------

  subroutine chain_retangulo 

    integer :: aux(1)

    allocate ( psi_r(nxy) )

    psi_r = 0.0d0

    do i = 1, nx-1
       do j = 2, ny-1
          np = (j-1)*nx + i
          psi_r(np) = psi_r(np-nx) + ue(np) * dy
       end do
    end do

    psi_r_min = minval(psi_r)

    aux = minloc(psi_r)
    xmin_psi_r = ( aux(1) - (aux(1)/nx)*nx - 1 ) * dx
    ymin_psi_r = ( aux(1)/nx ) * dy

  end subroutine chain_retangulo

! -----------------------------------------------

  subroutine chain_trapezio 

    real*8 :: x, y, fp, glp, uexe ! auxiliares
    integer :: aux(1)

    allocate ( psi_t(nxy), psi(nxy) )

    psi_t = 0.0d0

    do i = 1, nx-1

       j = 2
       np = (j-1)*nx + i
       psi_t(np) = psi_t(np-nx) + (2*ue(np) + ue(np+nx)) * dy / 4

       do j = 3, ny-2
          np = (j-1)*nx + i
          psi_t(np) = psi_t(np-nx) + (2*ue(np) + ue(np-nx) + ue(np+nx) ) * dy / 4
       end do

       j = ny-1
       np = (j-1)*nx + i
       if ( tipo == 1 ) then
          x = (i-1)*dx
          y = 1.0d0
          fp = (x**4) - 2.0d0*(x**3) + (x**2)
          glp = 4.0d0*(y**3) - 2.0d0*y
          uexe = 8.0d0*fp*glp
          psi_t(np) = psi_t(np-nx) + (2*ue(np) + ue(np-nx) + uexe ) * dy / 4
       else 
          psi_t(np) = psi_t(np-nx) + (2*ue(np) + ue(np-nx) + Ut ) * dy / 4
       end if

    end do

    psi_t_min = minval(psi_t)

    aux = minloc(psi_t)
    no_psi_min = aux(1)
    xmin_psi_t = ( aux(1) - (aux(1)/nx)*nx - 1 ) * dx
    ymin_psi_t = ( aux(1)/nx ) * dy
    x_psi_min = xmin_psi_t
    y_psi_min = ymin_psi_t
    
    psi = psi_t
    
  end subroutine chain_trapezio

! -----------------------------------------------

  subroutine perfil_chain_retangulo

    real*8  :: x

    write(10,1)
    1 format(//,5x,'Psi(x) ret?ngulo em y=1 (ideal = zero)',//, &
                t4,'x', t19,'Psi_r(x)')

    j = ny-1
    do i = 1, nx-1
       np = (j-1)*nx + i
       x  = (i-1)*dx
       write(10,2) x, Psi_r(np)
    end do

    2 format(2(1pe15.6))

  end subroutine perfil_chain_retangulo


!****************************************************************************************
 
 subroutine esc_psi_p2_2d
 
 ! Salva dados para entrada em interp_2d, p = 2, psi
    
    380 format(A<len(trim(adjustl(caso)))>)
    
    write(caso1,380) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso1)),"_psi_p2_2d",".txt" 
    
    open(185,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(185,129) arquivo_caso
    
    129 format (3x, A40, 19x, 'nome do caso') 
        
    write (185, 130) 9
   
    130 format (1x, i6,48x, '= N, número de nós para interpolação 2D, p = 2')				         
	
	write (185, 131) x_psi_min - dx, y_psi_min - dy, psi(no_psi_min - nx - 1), 1
	
	write(185, 131)  x_psi_min, y_psi_min - dy, psi(no_psi_min - nx), 2  
	
	write(185, 131) x_psi_min + dx, y_psi_min - dy, psi(no_psi_min - nx + 1), 3
	
	write(185, 131) x_psi_min - dx, y_psi_min, psi(no_psi_min - 1), 4
	
	write(185, 131) x_psi_min, y_psi_min, psi(no_psi_min), 5

	write(185, 131) x_psi_min + dx, y_psi_min, psi(no_psi_min + 1), 6

	write(185, 131) x_psi_min - dx, y_psi_min + dy, psi(no_psi_min + nx - 1), 7 
	
	write(185, 131) x_psi_min, y_psi_min + dy, psi(no_psi_min + nx), 8 
	
	write(185, 131) x_psi_min + dx, y_psi_min + dy, psi(no_psi_min + nx + 1), 9
	    
    131 format (1pe48.32, 1pe48.32, 1pe48.32, 2x, '(x, y, psi) do nó ', 1x, i3)		   
    
    close(185)
    
 end subroutine esc_psi_p2_2d
 
 
 !****************************************************************************************
 
 subroutine esc_psi_p4_2d
 
 ! Salva dados para entrada em interp_2d, p = 4, psi
    
    3800 format(A<len(trim(adjustl(caso)))>)
    
    write(caso,3800) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso)),"_psi_p4_2d",".txt" 
    
    open(1850,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(1850,1290) arquivo_caso
    
    1290 format (3x, A40, 19x, 'nome do caso') 
        
    write (1850, 1300) 25
   
    1300 format (1x, i6,48x, '= N, número de nós para interpolação 2D, p = 4')				         
	
    write (1850, 1310) x_psi_min - 2.0q0*dx, y_psi_min - 2.0q0*dy, psi(no_psi_min - 2*nx - 2), 1
	
	write(1850, 1310)  x_psi_min-dx, y_psi_min - 2.0q0*dy, psi(no_psi_min - 2*nx - 1), 2  
	
	write(1850, 1310) x_psi_min, y_psi_min - 2.0q0*dy, psi(no_psi_min - 2*nx), 3
	
	write(1850, 1310) x_psi_min + dx, y_psi_min - 2.0q0*dy, psi(no_psi_min - 2*nx + 1), 4
	
	write(1850, 1310) x_psi_min + 2.0q0*dx, y_psi_min - 2.0q0*dy, psi(no_psi_min - 2*nx + 2), 5

	write(1850, 1310) x_psi_min - 2.0q0*dx, y_psi_min - dy, psi(no_psi_min - nx - 2), 6 

	write(1850, 1310) x_psi_min - dx, y_psi_min - dy, psi(no_psi_min - nx - 1), 7 
	
	write(1850, 1310) x_psi_min, y_psi_min - dy, psi(no_psi_min - nx), 8 
	
	write(1850, 1310) x_psi_min + dx, y_psi_min - dy, psi(no_psi_min - nx + 1), 9
	
	write(1850, 1310) x_psi_min + 2.0q0*dx, y_psi_min - dy, psi(no_psi_min - nx + 2), 10
	
	write(1850, 1310) x_psi_min - 2.0q0*dx, y_psi_min, psi(no_psi_min - 2), 11
	
	write(1850, 1310) x_psi_min - dx, y_psi_min, psi(no_psi_min - 1), 12
	
	write(1850, 1310) x_psi_min, y_psi_min, psi(no_psi_min), 13
	
	write(1850, 1310) x_psi_min + dx, y_psi_min, psi(no_psi_min + 1), 14
	
	write(1850, 1310) x_psi_min + 2.0q0*dx, y_psi_min, psi(no_psi_min + 2), 15
	
	write(1850, 1310) x_psi_min - 2.0q0*dx, y_psi_min + dy, psi(no_psi_min + nx - 2), 16
	
	write(1850, 1310) x_psi_min - dx, y_psi_min + dy, psi(no_psi_min + nx - 1), 17
	
	write(1850, 1310) x_psi_min, y_psi_min + dy, psi(no_psi_min + nx), 18
	
	write(1850, 1310) x_psi_min + dx, y_psi_min + dy, psi(no_psi_min + nx + 1), 19
	
	write(1850, 1310) x_psi_min + 2.0q0*dx, y_psi_min + dy, psi(no_psi_min + nx + 2), 20
	
	write(1850, 1310) x_psi_min - 2.0q0*dx, y_psi_min + 2.0q0*dy, psi(no_psi_min + 2*nx - 2), 21
	
	write(1850, 1310) x_psi_min - dx, y_psi_min + 2.0q0*dy, psi(no_psi_min + 2*nx - 1), 22
	
	write(1850, 1310) x_psi_min, y_psi_min + 2.0q0*dy, psi(no_psi_min + 2*nx), 23
	
	write(1850, 1310) x_psi_min + dx, y_psi_min + 2.0q0*dy, psi(no_psi_min + 2*nx + 1), 24
	
	write(1850, 1310) x_psi_min + 2.0q0*dx, y_psi_min + 2.0q0*dy, psi(no_psi_min + 2*nx + 2), 25
	    
    1310 format (1pe48.32, 1pe48.32, 1pe48.32, 2x, 'psi(x, y, psi) do nó ', 1x, i3)		   
    
    close(1850)
    
 end subroutine esc_psi_p4_2d
 
  !****************************************************************************************
 
 subroutine esc_psi_p6_2d
 
 ! Salva dados para entrada em interp_2d, p = 6, psi
    
    33800 format(A<len(trim(adjustl(caso)))>)
    
    write(caso,33800) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso)),"_psi_p6_2d",".txt" 
    
    open(11850,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(11850,11290) arquivo_caso
    
    11290 format (3x, A40, 19x, 'nome do caso') 
        
    write (11850, 11300) 49
   
    11300 format (1x, i6,48x, '= N, número de nós para interpolação 2D, p = 6')	
    	         
	        
    write (11850, 11310) x_psi_min - 3.0q0*dx, y_psi_min - 3.0q0*dy, psi(no_psi_min - 3*nx - 3), 1
	
	write (11850, 11310) x_psi_min - 2.0q0*dx, y_psi_min - 3.0q0*dy, psi(no_psi_min - 3*nx - 2), 2
	
	write (11850, 11310) x_psi_min - dx, y_psi_min - 3.0q0*dy, psi(no_psi_min - 3*nx - 1), 3
	
	write (11850, 11310) x_psi_min, y_psi_min - 3.0q0*dy, psi(no_psi_min - 3*nx), 4
	
	write (11850, 11310) x_psi_min + dx, y_psi_min - 3.0q0*dy, psi(no_psi_min - 3*nx +1), 5 
	
	write (11850, 11310) x_psi_min + 2.0q0*dx, y_psi_min - 3.0q0*dy, psi(no_psi_min - 3*nx +2), 6 
	
	write (11850, 11310) x_psi_min + 3.0q0*dx, y_psi_min - 3.0q0*dy, psi(no_psi_min - 3*nx +3), 7 
	
    write (11850, 11310) x_psi_min - 3.0q0*dx, y_psi_min - 2.0q0*dy, psi(no_psi_min - 2*nx - 3), 8
	
	write (11850, 11310) x_psi_min - 2.0q0*dx, y_psi_min - 2.0q0*dy, psi(no_psi_min - 2*nx - 2), 9
	
	write (11850, 11310) x_psi_min - dx, y_psi_min - 2.0q0*dy, psi(no_psi_min - 2*nx - 1), 10
	
	write (11850, 11310) x_psi_min, y_psi_min - 2.0q0*dy, psi(no_psi_min - 2*nx), 11
	
	write (11850, 11310) x_psi_min + dx, y_psi_min - 2.0q0*dy, psi(no_psi_min - 2*nx +1), 12 
	
	write (11850, 11310) x_psi_min + 2.0q0*dx, y_psi_min - 2.0q0*dy, psi(no_psi_min - 2*nx +2), 13 
	
	write (11850, 11310) x_psi_min + 3.0q0*dx, y_psi_min - 2.0q0*dy, psi(no_psi_min - 2*nx +3), 14 
	
	write (11850, 11310) x_psi_min - 3.0q0*dx, y_psi_min - dy, psi(no_psi_min - nx - 3), 15
	
	write (11850, 11310) x_psi_min - 2.0q0*dx, y_psi_min - dy, psi(no_psi_min - nx - 2), 16
	
	write (11850, 11310) x_psi_min - dx, y_psi_min - dy, psi(no_psi_min - nx - 1), 17
	
	write (11850, 11310) x_psi_min, y_psi_min - dy, psi(no_psi_min - nx), 18
	
	write (11850, 11310) x_psi_min + dx, y_psi_min - dy, psi(no_psi_min - nx +1), 19 
	
	write (11850, 11310) x_psi_min + 2.0q0*dx, y_psi_min - dy, psi(no_psi_min - nx +2), 20 
	
	write (11850, 11310) x_psi_min + 3.0q0*dx, y_psi_min - dy, psi(no_psi_min - nx +3), 21
	
	write (11850, 11310) x_psi_min - 3.0q0*dx, y_psi_min, psi(no_psi_min - 3), 22
	
	write (11850, 11310) x_psi_min - 2.0q0*dx, y_psi_min, psi(no_psi_min - 2), 23
	
	write (11850, 11310) x_psi_min - dx, y_psi_min, psi(no_psi_min - 1), 24
	
	write (11850, 11310) x_psi_min, y_psi_min, psi(no_psi_min), 25
	
	write (11850, 11310) x_psi_min + dx, y_psi_min, psi(no_psi_min + 1), 26 
	
	write (11850, 11310) x_psi_min + 2.0q0*dx, y_psi_min, psi(no_psi_min + 2), 27 
	
	write (11850, 11310) x_psi_min + 3.0q0*dx, y_psi_min, psi(no_psi_min + 3), 28
	
	write (11850, 11310) x_psi_min - 3.0q0*dx, y_psi_min+dy, psi(no_psi_min + nx - 3), 29
	
	write (11850, 11310) x_psi_min - 2.0q0*dx, y_psi_min+dy, psi(no_psi_min + nx - 2), 30
	
	write (11850, 11310) x_psi_min - dx, y_psi_min+dy, psi(no_psi_min + nx - 1), 31
	
	write (11850, 11310) x_psi_min, y_psi_min+dy, psi(no_psi_min + nx), 32
	
	write (11850, 11310) x_psi_min + dx, y_psi_min+dy, psi(no_psi_min + nx + 1), 33
	
	write (11850, 11310) x_psi_min + 2.0q0*dx, y_psi_min+dy, psi(no_psi_min + nx + 2), 34
	
	write (11850, 11310) x_psi_min + 3.0q0*dx, y_psi_min+dy, psi(no_psi_min + nx + 3), 35
	
	write (11850, 11310) x_psi_min - 3.0q0*dx, y_psi_min+2.0q0*dy, psi(no_psi_min + 2*nx - 3), 36
	
	write (11850, 11310) x_psi_min - 2.0q0*dx, y_psi_min+2.0q0*dy, psi(no_psi_min + 2*nx - 2), 37
	
	write (11850, 11310) x_psi_min - dx, y_psi_min+2.0q0*dy, psi(no_psi_min + 2*nx - 1), 38
	
	write (11850, 11310) x_psi_min, y_psi_min+2.0q0*dy, psi(no_psi_min + 2*nx), 39
	
	write (11850, 11310) x_psi_min + dx, y_psi_min+2.0q0*dy, psi(no_psi_min + 2*nx + 1), 40
	
	write (11850, 11310) x_psi_min + 2.0q0*dx, y_psi_min+2.0q0*dy, psi(no_psi_min + 2*nx + 2), 41
	
	write (11850, 11310) x_psi_min + 3.0q0*dx, y_psi_min+2.0q0*dy, psi(no_psi_min + 2*nx + 3), 42
	
	write (11850, 11310) x_psi_min - 3.0q0*dx, y_psi_min+3.0q0*dy, psi(no_psi_min + 3*nx - 3), 43
	
	write (11850, 11310) x_psi_min - 2.0q0*dx, y_psi_min+3.0q0*dy, psi(no_psi_min + 3*nx - 2), 44
	
	write (11850, 11310) x_psi_min - dx, y_psi_min+3.0q0*dy, psi(no_psi_min + 3*nx - 1), 45
	
	write (11850, 11310) x_psi_min, y_psi_min+3.0q0*dy, psi(no_psi_min + 3*nx), 46
	
	write (11850, 11310) x_psi_min + dx, y_psi_min+3.0q0*dy, psi(no_psi_min + 3*nx + 1), 47
	
	write (11850, 11310) x_psi_min + 2.0q0*dx, y_psi_min+3.0q0*dy, psi(no_psi_min + 3*nx + 2), 48
	
	write (11850, 11310) x_psi_min + 3.0q0*dx, y_psi_min+3.0q0*dy, psi(no_psi_min + 3*nx + 3), 49
	    
    11310 format (1pe48.32, 1pe48.32, 1pe48.32, 2x, 'psi(x, y, psi) do nó ', 1x, i3)		   
    
    close(11850)
    
 end subroutine esc_psi_p6_2d
 
!************************************* ********************************************************
	
subroutine extremos

!Obtém extremos 2D da malha (valores nodais)

x1 = 0.5q0*dx

y1 = x1  !0.5q0*dy

u_min_2d = u(2) !inicializa valor mínimo da componente na direção horizontal do vetor velocidade (2D)

v_min_2d = v(2) !inicializa valor mínimo da componente na direção horizontal do vetor velocidade (2D)

v_max_2d = v(2) !inicializa valor mínimo da componente na direção horizontal do vetor velocidade (2D)

	do j = 2,ny-1

		yy = y1 +(j-2)*dy

        k1 = (j-1)*nx 
        
		do i = 2, nx-1
				 
				xx = x1 +(i-2)*dx     
								
				np = k1 + i
				
									
					if ( u_min_2d > u(np) ) then
					
					     u_min_2d = u(np) ! valor mínimo de u no domínio 2D
					     
					     no_umin_2d = np  ! nó correspondente
					     
					     x_umin_2d = xx !abcissa do nó correspondente
					     
					     y_umin_2d = yy  !ordenada do nó correspondente
					
					end if	
					
					if ( v_min_2d > v(np) ) then
					
					     v_min_2d = v(np) ! valor mínimo de v no domínio 2D
					     
					     no_vmin_2d = np  ! nó correspondente
					     
					     x_vmin_2d = xx !abcissa do nó correspondente
					     
					     y_vmin_2d = yy  !ordenada do nó correspondente
					
					end if	
					
					if ( v_max_2d < v(np) ) then
					
					     v_max_2d = v(np) ! valor máximo de v no domínio 2D
					     
					     no_vmax_2d = np  ! nó correspondente
					     
					     x_vmax_2d = xx !abcissa do nó correspondente
					     
					     y_vmax_2d = yy  !ordenada do nó correspondente
					
					end if	
						
		    end do
										
	    end do

 
 end subroutine extremos

!*************************************************************************************************************
 
 subroutine esc_min_u_p2_2d
 
 ! Salva dados para entrada em interp_2d, p = 2, u
    
    350 format(A<len(trim(adjustl(caso)))>)
    
    write(caso4,350) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso)),"_min_u_p2_2d",".txt" 
    
    open(145,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(145,120) arquivo_caso
    
    120 format (3x, A40, 19x, 'nome do caso') 
        
    write (145, 121) 9
   
    121 format (1x, i6,48x, '= N, número de nós para interpolação, 2D, p = 2')
	
	write (145, 122) x_umin_2d - dx, y_umin_2d - dy, u(no_umin_2d - nx - 1), 1
	
	write(145, 122)  x_umin_2d, y_umin_2d - dy, u(no_umin_2d - nx), 2  
	
	write(145, 122) x_umin_2d + dx, y_umin_2d - dy, u(no_umin_2d - nx + 1), 3
	
	write(145, 122) x_umin_2d - dx, y_umin_2d, u(no_umin_2d - 1), 4
	
	write(145, 122) x_umin_2d, y_umin_2d, u(no_umin_2d), 5

	write(145, 122) x_umin_2d + dx, y_umin_2d, u(no_umin_2d + 1), 6

	write(145, 122) x_umin_2d - dx, y_umin_2d + dy, u(no_umin_2d +nx - 1), 7 
	
	write(145, 122) x_umin_2d, y_umin_2d + dy, u(no_umin_2d + nx), 8 
	
	write(145, 122) x_umin_2d + dx, y_umin_2d + dy, u(no_umin_2d + nx + 1), 9
	    
    122 format (1pe48.32, 1pe48.32, 1pe48.32, 2x, '(x, y, u) do nó ', 1x, i3)		   
    
    close(145)
    
 end subroutine esc_min_u_p2_2d
 
 !****************************************************************************************
 
 subroutine esc_min_u_p4_2d
 
 ! Salva dados para entrada em interp_2d, p = 4, u
    
    3500 format(A<len(trim(adjustl(caso)))>)
    
    write(caso,3500) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso)),"_min_u_p4_2d",".txt" 
    
    open(1450,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(1450,1200) arquivo_caso
    
    1200 format (3x, A40, 19x, 'nome do caso') 
        
    write (1450, 1210) 25
   
    1210 format (1x, i6,48x, '= N, número de nós para interpolação')
	
	write (1450, 1220) x_umin_2d - 2.0q0*dx, y_umin_2d - 2.0q0*dy, u(no_umin_2d - 2*nx - 2), 1
	
	write(1450, 1220)  x_umin_2d-dx, y_umin_2d - 2.0q0*dy, u(no_umin_2d - 2*nx - 1), 2  
	
	write(1450, 1220) x_umin_2d, y_umin_2d - 2.0q0*dy, u(no_umin_2d - 2*nx), 3
	
	write(1450, 1220) x_umin_2d + dx, y_umin_2d - 2.0q0*dy, u(no_umin_2d - 2*nx + 1), 4
	
	write(1450, 1220) x_umin_2d + 2.0q0*dx, y_umin_2d - 2.0q0*dy, u(no_umin_2d - 2*nx + 2), 5

	write(1450, 1220) x_umin_2d - 2.0q0*dx, y_umin_2d - dy, u(no_umin_2d - nx - 2), 6 

	write(1450, 1220) x_umin_2d - dx, y_umin_2d - dy, u(no_umin_2d - nx - 1), 7 
	
	write(1450, 1220) x_umin_2d, y_umin_2d - dy, u(no_umin_2d - nx), 8 
	
	write(1450, 1220) x_umin_2d + dx, y_umin_2d - dy, u(no_umin_2d - nx + 1), 9
	
	write(1450, 1220) x_umin_2d + 2.0q0*dx, y_umin_2d - dy, u(no_umin_2d - nx + 2), 10
	
	write(1450, 1220) x_umin_2d - 2.0q0*dx, y_umin_2d, u(no_umin_2d - 2), 11
	
	write(1450, 1220) x_umin_2d - dx, y_umin_2d, u(no_umin_2d - 1), 12
	
	write(1450, 1220) x_umin_2d, y_umin_2d, u(no_umin_2d), 13
	
	write(1450, 1220) x_umin_2d + dx, y_umin_2d, u(no_umin_2d + 1), 14
	
	write(1450, 1220) x_umin_2d + 2.0q0*dx, y_umin_2d, u(no_umin_2d + 2), 15
	
	write(1450, 1220) x_umin_2d - 2.0q0*dx, y_umin_2d + dy, u(no_umin_2d + nx - 2), 16
	
	write(1450, 1220) x_umin_2d - dx, y_umin_2d + dy, u(no_umin_2d + nx - 1), 17
	
	write(1450, 1220) x_umin_2d, y_umin_2d + dy, u(no_umin_2d + nx), 18
	
	write(1450, 1220) x_umin_2d + dx, y_umin_2d + dy, u(no_umin_2d + nx + 1), 19
	
	write(1450, 1220) x_umin_2d + 2.0q0*dx, y_umin_2d + dy, u(no_umin_2d + nx + 2), 20
	
	write(1450, 1220) x_umin_2d - 2.0q0*dx, y_umin_2d + 2.0q0*dy, u(no_umin_2d + 2*nx - 2), 21
	
	write(1450, 1220) x_umin_2d - dx, y_umin_2d + 2.0q0*dy, u(no_umin_2d + 2*nx - 1), 22
	
	write(1450, 1220) x_umin_2d, y_umin_2d + 2.0q0*dy, u(no_umin_2d + 2*nx), 23
	
	write(1450, 1220) x_umin_2d + dx, y_umin_2d + 2.0q0*dy, u(no_umin_2d + 2*nx + 1), 24
	
	write(1450, 1220) x_umin_2d + 2.0q0*dx, y_umin_2d + 2.0q0*dy, u(no_umin_2d + 2*nx + 2), 25
	    
    1220 format (1pe48.32, 1pe48.32, 1pe48.32, 2x, 'u(x, y, u) do nó ', 1x, i3)		   
    
    close(1450)
    
 end subroutine esc_min_u_p4_2d
 
 !****************************************************************************************
 
 subroutine esc_min_u_p6_2d
 
 ! Salva dados para entrada em interp_2d, p = 6, u
    
    33500 format(A<len(trim(adjustl(caso)))>)
    
    write(caso,33500) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso)),"_min_u_p6_2d",".txt" 
    
    open(11450,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(11450,11200) arquivo_caso
    
    11200 format (3x, A40, 19x, 'nome do caso') 
        
    write (11450, 11210) 49
   
    11210 format (1x, i6,48x, '= N, número de nós para interpolação')
	
	write (11450, 11220) x_umin_2d - 3.0q0*dx, y_umin_2d - 3.0q0*dy, u(no_umin_2d - 3*nx - 3), 1
	
	write (11450, 11220) x_umin_2d - 2.0q0*dx, y_umin_2d - 3.0q0*dy, u(no_umin_2d - 3*nx - 2), 2
	
	write (11450, 11220) x_umin_2d - dx, y_umin_2d - 3.0q0*dy, u(no_umin_2d - 3*nx - 1), 3
	
	write (11450, 11220) x_umin_2d, y_umin_2d - 3.0q0*dy, u(no_umin_2d - 3*nx), 4
	
	write (11450, 11220) x_umin_2d + dx, y_umin_2d - 3.0q0*dy, u(no_umin_2d - 3*nx +1), 5 
	
	write (11450, 11220) x_umin_2d + 2.0q0*dx, y_umin_2d - 3.0q0*dy, u(no_umin_2d - 3*nx +2), 6 
	
	write (11450, 11220) x_umin_2d + 3.0q0*dx, y_umin_2d - 3.0q0*dy, u(no_umin_2d - 3*nx +3), 7 
	
    write (11450, 11220) x_umin_2d - 3.0q0*dx, y_umin_2d - 2.0q0*dy, u(no_umin_2d - 2*nx - 3), 8
	
	write (11450, 11220) x_umin_2d - 2.0q0*dx, y_umin_2d - 2.0q0*dy, u(no_umin_2d - 2*nx - 2), 9
	
	write (11450, 11220) x_umin_2d - dx, y_umin_2d - 2.0q0*dy, u(no_umin_2d - 2*nx - 1), 10
	
	write (11450, 11220) x_umin_2d, y_umin_2d - 2.0q0*dy, u(no_umin_2d - 2*nx), 11
	
	write (11450, 11220) x_umin_2d + dx, y_umin_2d - 2.0q0*dy, u(no_umin_2d - 2*nx +1), 12 
	
	write (11450, 11220) x_umin_2d + 2.0q0*dx, y_umin_2d - 2.0q0*dy, u(no_umin_2d - 2*nx +2), 13 
	
	write (11450, 11220) x_umin_2d + 3.0q0*dx, y_umin_2d - 2.0q0*dy, u(no_umin_2d - 2*nx +3), 14 
	
	write (11450, 11220) x_umin_2d - 3.0q0*dx, y_umin_2d - dy, u(no_umin_2d - nx - 3), 15
	
	write (11450, 11220) x_umin_2d - 2.0q0*dx, y_umin_2d - dy, u(no_umin_2d - nx - 2), 16
	
	write (11450, 11220) x_umin_2d - dx, y_umin_2d - dy, u(no_umin_2d - nx - 1), 17
	
	write (11450, 11220) x_umin_2d, y_umin_2d - dy, u(no_umin_2d - nx), 18
	
	write (11450, 11220) x_umin_2d + dx, y_umin_2d - dy, u(no_umin_2d - nx +1), 19 
	
	write (11450, 11220) x_umin_2d + 2.0q0*dx, y_umin_2d - dy, u(no_umin_2d - nx +2), 20 
	
	write (11450, 11220) x_umin_2d + 3.0q0*dx, y_umin_2d - dy, u(no_umin_2d - nx +3), 21
	
	write (11450, 11220) x_umin_2d - 3.0q0*dx, y_umin_2d, u(no_umin_2d  - 3), 22
	
	write (11450, 11220) x_umin_2d - 2.0q0*dx, y_umin_2d, u(no_umin_2d  - 2), 23
	
	write (11450, 11220) x_umin_2d - dx, y_umin_2d, u(no_umin_2d - 1), 24
	
	write (11450, 11220) x_umin_2d, y_umin_2d, u(no_umin_2d), 25
	
	write (11450, 11220) x_umin_2d + dx, y_umin_2d, u(no_umin_2d + 1), 26 
	
	write (11450, 11220) x_umin_2d + 2.0q0*dx, y_umin_2d, u(no_umin_2d + 2), 27 
	
	write (11450, 11220) x_umin_2d + 3.0q0*dx, y_umin_2d, u(no_umin_2d + 3), 28
	
	write (11450, 11220) x_umin_2d - 3.0q0*dx, y_umin_2d+dy, u(no_umin_2d + nx - 3), 29
	
	write (11450, 11220) x_umin_2d - 2.0q0*dx, y_umin_2d+dy, u(no_umin_2d + nx - 2), 30
	
	write (11450, 11220) x_umin_2d - dx, y_umin_2d+dy, u(no_umin_2d + nx - 1), 31
	
	write (11450, 11220) x_umin_2d, y_umin_2d+dy, u(no_umin_2d + nx), 32
	
	write (11450, 11220) x_umin_2d + dx, y_umin_2d+dy, u(no_umin_2d + nx + 1), 33
	
	write (11450, 11220) x_umin_2d + 2.0q0*dx, y_umin_2d+dy, u(no_umin_2d + nx + 2), 34
	
	write (11450, 11220) x_umin_2d + 3.0q0*dx, y_umin_2d+dy, u(no_umin_2d + nx + 3), 35
	
	write (11450, 11220) x_umin_2d - 3.0q0*dx, y_umin_2d+2.0q0*dy, u(no_umin_2d + 2*nx - 3), 36
	
	write (11450, 11220) x_umin_2d - 2.0q0*dx, y_umin_2d+2.0q0*dy, u(no_umin_2d + 2*nx - 2), 37
	
	write (11450, 11220) x_umin_2d - dx, y_umin_2d+2.0q0*dy, u(no_umin_2d + 2*nx - 1), 38
	
	write (11450, 11220) x_umin_2d, y_umin_2d+2.0q0*dy, u(no_umin_2d + 2*nx), 39
	
	write (11450, 11220) x_umin_2d + dx, y_umin_2d+2.0q0*dy, u(no_umin_2d + 2*nx + 1), 40
	
	write (11450, 11220) x_umin_2d + 2.0q0*dx, y_umin_2d+2.0q0*dy, u(no_umin_2d + 2*nx + 2), 41
	
	write (11450, 11220) x_umin_2d + 3.0q0*dx, y_umin_2d+2.0q0*dy, u(no_umin_2d + 2*nx + 3), 42
	
	write (11450, 11220) x_umin_2d - 3.0q0*dx, y_umin_2d+3.0q0*dy, u(no_umin_2d + 3*nx - 3), 43
	
	write (11450, 11220) x_umin_2d - 2.0q0*dx, y_umin_2d+3.0q0*dy, u(no_umin_2d + 3*nx - 2), 44
	
	write (11450, 11220) x_umin_2d - dx, y_umin_2d+3.0q0*dy, u(no_umin_2d + 3*nx - 1), 45
	
	write (11450, 11220) x_umin_2d, y_umin_2d+3.0q0*dy, u(no_umin_2d + 3*nx), 46
	
	write (11450, 11220) x_umin_2d + dx, y_umin_2d+3.0q0*dy, u(no_umin_2d + 3*nx + 1), 47
	
	write (11450, 11220) x_umin_2d + 2.0q0*dx, y_umin_2d+3.0q0*dy, u(no_umin_2d + 3*nx + 2), 48
	
	write (11450, 11220) x_umin_2d + 3.0q0*dx, y_umin_2d+3.0q0*dy, u(no_umin_2d + 3*nx + 3), 49
	    
    11220 format (1pe48.32, 1pe48.32, 1pe48.32, 2x, 'u(x, y, u) do nó ', 1x, i3)		   
    
    close(11450)
    
 end subroutine esc_min_u_p6_2d

!*************************************************************************************************************
 
 subroutine esc_min_v_p2_2d
 
 ! Salva dados para entrada em interp_2d, p = 2, min v
    
    350 format(A<len(trim(adjustl(caso)))>)
    
    write(caso4,350) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso)),"_min_v_p2_2d",".txt" 
    
    open(145,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(145,120) arquivo_caso
    
    120 format (3x, A40, 19x, 'nome do caso') 
        
    write (145, 121) 9
   
    121 format (1x, i6,48x, '= N, número de nós para interpolação, 2D, p = 2')
	
	write (145, 122) x_vmin_2d - dx, y_vmin_2d - dy, v(no_vmin_2d - nx - 1), 1
	
	write(145, 122)  x_vmin_2d, y_vmin_2d - dy, v(no_vmin_2d - nx), 2  
	
	write(145, 122) x_vmin_2d + dx, y_vmin_2d - dy, v(no_vmin_2d - nx + 1), 3
	
	write(145, 122) x_vmin_2d - dx, y_vmin_2d, v(no_vmin_2d - 1), 4
	
	write(145, 122) x_vmin_2d, y_vmin_2d, v(no_vmin_2d), 5

	write(145, 122) x_vmin_2d + dx, y_vmin_2d, v(no_vmin_2d + 1), 6

	write(145, 122) x_vmin_2d - dx, y_vmin_2d + dy, v(no_vmin_2d +nx - 1), 7 
	
	write(145, 122) x_vmin_2d, y_vmin_2d + dy, v(no_vmin_2d + nx), 8 
	
	write(145, 122) x_vmin_2d + dx, y_vmin_2d + dy, v(no_vmin_2d + nx + 1), 9
	    
    122 format (1pe48.32, 1pe48.32, 1pe48.32, 2x, '(x, y, v) do nó ', 1x, i3)		   
    
    close(145)
    
 end subroutine esc_min_v_p2_2d
 
 !****************************************************************************************
 
 subroutine esc_min_v_p4_2d
 
 ! Salva dados para entrada em interp_2d, p = 4, min v
    
    3500 format(A<len(trim(adjustl(caso)))>)
    
    write(caso,3500) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso)),"_min_v_p4_2d",".txt" 
    
    open(1450,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(1450,1200) arquivo_caso
    
    1200 format (3x, A40, 19x, 'nome do caso') 
        
    write (1450, 1210) 25
   
    1210 format (1x, i6,48x, '= N, número de nós para interpolação')
	
	write (1450, 1220) x_vmin_2d - 2.0q0*dx, y_vmin_2d - 2.0q0*dy, v(no_vmin_2d - 2*nx - 2), 1
	
	write(1450, 1220)  x_vmin_2d-dx, y_vmin_2d - 2.0q0*dy, v(no_vmin_2d - 2*nx - 1), 2  
	
	write(1450, 1220) x_vmin_2d, y_vmin_2d - 2.0q0*dy, v(no_vmin_2d - 2*nx), 3
	
	write(1450, 1220) x_vmin_2d + dx, y_vmin_2d - 2.0q0*dy, v(no_vmin_2d - 2*nx + 1), 4
	
	write(1450, 1220) x_vmin_2d + 2.0q0*dx, y_vmin_2d - 2.0q0*dy, v(no_vmin_2d - 2*nx + 2), 5

	write(1450, 1220) x_vmin_2d - 2.0q0*dx, y_vmin_2d - dy, v(no_vmin_2d - nx - 2), 6 

	write(1450, 1220) x_vmin_2d - dx, y_vmin_2d - dy, v(no_vmin_2d - nx - 1), 7 
	
	write(1450, 1220) x_vmin_2d, y_vmin_2d - dy, v(no_vmin_2d - nx), 8 
	
	write(1450, 1220) x_vmin_2d + dx, y_vmin_2d - dy, v(no_vmin_2d - nx + 1), 9
	
	write(1450, 1220) x_vmin_2d + 2.0q0*dx, y_vmin_2d - dy, v(no_vmin_2d - nx + 2), 10
	
	write(1450, 1220) x_vmin_2d - 2.0q0*dx, y_vmin_2d, v(no_vmin_2d - 2), 11
	
	write(1450, 1220) x_vmin_2d - dx, y_vmin_2d, v(no_vmin_2d - 1), 12
	
	write(1450, 1220) x_vmin_2d, y_vmin_2d, v(no_vmin_2d), 13
	
	write(1450, 1220) x_vmin_2d + dx, y_vmin_2d, v(no_vmin_2d + 1), 14
	
	write(1450, 1220) x_vmin_2d + 2.0q0*dx, y_vmin_2d, v(no_vmin_2d + 2), 15
	
	write(1450, 1220) x_vmin_2d - 2.0q0*dx, y_vmin_2d + dy, v(no_vmin_2d + nx - 2), 16
	
	write(1450, 1220) x_vmin_2d - dx, y_vmin_2d + dy, v(no_vmin_2d + nx - 1), 17
	
	write(1450, 1220) x_vmin_2d, y_vmin_2d + dy, v(no_vmin_2d + nx), 18
	
	write(1450, 1220) x_vmin_2d + dx, y_vmin_2d + dy, v(no_vmin_2d + nx + 1), 19
	
	write(1450, 1220) x_vmin_2d + 2.0q0*dx, y_vmin_2d + dy, v(no_vmin_2d + nx + 2), 20
	
	write(1450, 1220) x_vmin_2d - 2.0q0*dx, y_vmin_2d + 2.0q0*dy, v(no_vmin_2d + 2*nx - 2), 21
	
	write(1450, 1220) x_vmin_2d - dx, y_vmin_2d + 2.0q0*dy, v(no_vmin_2d + 2*nx - 1), 22
	
	write(1450, 1220) x_vmin_2d, y_vmin_2d + 2.0q0*dy, v(no_vmin_2d + 2*nx), 23
	
	write(1450, 1220) x_vmin_2d + dx, y_vmin_2d + 2.0q0*dy, v(no_vmin_2d + 2*nx + 1), 24
	
	write(1450, 1220) x_vmin_2d + 2.0q0*dx, y_vmin_2d + 2.0q0*dy, v(no_vmin_2d + 2*nx + 2), 25
	    
    1220 format (1pe48.32, 1pe48.32, 1pe48.32, 2x, 'v(x, y, v) do nó ', 1x, i3)		   
    
    close(1450)
    
 end subroutine esc_min_v_p4_2d
 
 !****************************************************************************************
 
 subroutine esc_min_v_p6_2d
 
 ! Salva dados para entrada em interp_2d, p = 6, min v
    
    33500 format(A<len(trim(adjustl(caso)))>)
    
    write(caso,33500) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso)),"_min_v_p6_2d",".txt" 
    
    open(11450,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(11450,11200) arquivo_caso
    
    11200 format (3x, A40, 19x, 'nome do caso') 
        
    write (11450, 11210) 49
   
    11210 format (1x, i6,48x, '= N, número de nós para interpolação')
	
	write (11450, 11220) x_vmin_2d - 3.0q0*dx, y_vmin_2d - 3.0q0*dy, v(no_vmin_2d - 3*nx - 3), 1
	
	write (11450, 11220) x_vmin_2d - 2.0q0*dx, y_vmin_2d - 3.0q0*dy, v(no_vmin_2d - 3*nx - 2), 2
	
	write (11450, 11220) x_vmin_2d - dx, y_vmin_2d - 3.0q0*dy, v(no_vmin_2d - 3*nx - 1), 3
	
	write (11450, 11220) x_vmin_2d, y_vmin_2d - 3.0q0*dy, v(no_vmin_2d - 3*nx), 4
	
	write (11450, 11220) x_vmin_2d + dx, y_vmin_2d - 3.0q0*dy, v(no_vmin_2d - 3*nx +1), 5 
	
	write (11450, 11220) x_vmin_2d + 2.0q0*dx, y_vmin_2d - 3.0q0*dy, v(no_vmin_2d - 3*nx +2), 6 
	
	write (11450, 11220) x_vmin_2d + 3.0q0*dx, y_vmin_2d - 3.0q0*dy, v(no_vmin_2d - 3*nx +3), 7 
	
    write (11450, 11220) x_vmin_2d - 3.0q0*dx, y_vmin_2d - 2.0q0*dy, v(no_vmin_2d - 2*nx - 3), 8
	
	write (11450, 11220) x_vmin_2d - 2.0q0*dx, y_vmin_2d - 2.0q0*dy, v(no_vmin_2d - 2*nx - 2), 9
	
	write (11450, 11220) x_vmin_2d - dx, y_vmin_2d - 2.0q0*dy, v(no_vmin_2d - 2*nx - 1), 10
	
	write (11450, 11220) x_vmin_2d, y_vmin_2d - 2.0q0*dy, v(no_vmin_2d - 2*nx), 11
	
	write (11450, 11220) x_vmin_2d + dx, y_vmin_2d - 2.0q0*dy, v(no_vmin_2d - 2*nx +1), 12 
	
	write (11450, 11220) x_vmin_2d + 2.0q0*dx, y_vmin_2d - 2.0q0*dy, v(no_vmin_2d - 2*nx +2), 13 
	
	write (11450, 11220) x_vmin_2d + 3.0q0*dx, y_vmin_2d - 2.0q0*dy, v(no_vmin_2d - 2*nx +3), 14 
	
	write (11450, 11220) x_vmin_2d - 3.0q0*dx, y_vmin_2d - dy, v(no_vmin_2d - nx - 3), 15
	
	write (11450, 11220) x_vmin_2d - 2.0q0*dx, y_vmin_2d - dy, v(no_vmin_2d - nx - 2), 16
	
	write (11450, 11220) x_vmin_2d - dx, y_vmin_2d - dy, v(no_vmin_2d - nx - 1), 17
	
	write (11450, 11220) x_vmin_2d, y_vmin_2d - dy, v(no_vmin_2d - nx), 18
	
	write (11450, 11220) x_vmin_2d + dx, y_vmin_2d - dy, v(no_vmin_2d - nx +1), 19 
	
	write (11450, 11220) x_vmin_2d + 2.0q0*dx, y_vmin_2d - dy, v(no_vmin_2d - nx +2), 20 
	
	write (11450, 11220) x_vmin_2d + 3.0q0*dx, y_vmin_2d - dy, v(no_vmin_2d - nx +3), 21
	
	write (11450, 11220) x_vmin_2d - 3.0q0*dx, y_vmin_2d, v(no_vmin_2d  - 3), 22
	
	write (11450, 11220) x_vmin_2d - 2.0q0*dx, y_vmin_2d, v(no_vmin_2d  - 2), 23
	
	write (11450, 11220) x_vmin_2d - dx, y_vmin_2d, v(no_vmin_2d - 1), 24
	
	write (11450, 11220) x_vmin_2d, y_vmin_2d, v(no_vmin_2d), 25
	
	write (11450, 11220) x_vmin_2d + dx, y_vmin_2d, v(no_vmin_2d + 1), 26 
	
	write (11450, 11220) x_vmin_2d + 2.0q0*dx, y_vmin_2d, v(no_vmin_2d + 2), 27 
	
	write (11450, 11220) x_vmin_2d + 3.0q0*dx, y_vmin_2d, v(no_vmin_2d + 3), 28
	
	write (11450, 11220) x_vmin_2d - 3.0q0*dx, y_vmin_2d+dy, v(no_vmin_2d + nx - 3), 29
	
	write (11450, 11220) x_vmin_2d - 2.0q0*dx, y_vmin_2d+dy, v(no_vmin_2d + nx - 2), 30
	
	write (11450, 11220) x_vmin_2d - dx, y_vmin_2d+dy, v(no_vmin_2d + nx - 1), 31
	
	write (11450, 11220) x_vmin_2d, y_vmin_2d+dy, v(no_vmin_2d + nx), 32
	
	write (11450, 11220) x_vmin_2d + dx, y_vmin_2d+dy, v(no_vmin_2d + nx + 1), 33
	
	write (11450, 11220) x_vmin_2d + 2.0q0*dx, y_vmin_2d+dy, v(no_vmin_2d + nx + 2), 34
	
	write (11450, 11220) x_vmin_2d + 3.0q0*dx, y_vmin_2d+dy, v(no_vmin_2d + nx + 3), 35
	
	write (11450, 11220) x_vmin_2d - 3.0q0*dx, y_vmin_2d+2.0q0*dy, v(no_vmin_2d + 2*nx - 3), 36
	
	write (11450, 11220) x_vmin_2d - 2.0q0*dx, y_vmin_2d+2.0q0*dy, v(no_vmin_2d + 2*nx - 2), 37
	
	write (11450, 11220) x_vmin_2d - dx, y_vmin_2d+2.0q0*dy, v(no_vmin_2d + 2*nx - 1), 38
	
	write (11450, 11220) x_vmin_2d, y_vmin_2d+2.0q0*dy, v(no_vmin_2d + 2*nx), 39
	
	write (11450, 11220) x_vmin_2d + dx, y_vmin_2d+2.0q0*dy, v(no_vmin_2d + 2*nx + 1), 40
	
	write (11450, 11220) x_vmin_2d + 2.0q0*dx, y_vmin_2d+2.0q0*dy, v(no_vmin_2d + 2*nx + 2), 41
	
	write (11450, 11220) x_vmin_2d + 3.0q0*dx, y_vmin_2d+2.0q0*dy, v(no_vmin_2d + 2*nx + 3), 42
	
	write (11450, 11220) x_vmin_2d - 3.0q0*dx, y_vmin_2d+3.0q0*dy, v(no_vmin_2d + 3*nx - 3), 43
	
	write (11450, 11220) x_vmin_2d - 2.0q0*dx, y_vmin_2d+3.0q0*dy, v(no_vmin_2d + 3*nx - 2), 44
	
	write (11450, 11220) x_vmin_2d - dx, y_vmin_2d+3.0q0*dy, v(no_vmin_2d + 3*nx - 1), 45
	
	write (11450, 11220) x_vmin_2d, y_vmin_2d+3.0q0*dy, v(no_vmin_2d + 3*nx), 46
	
	write (11450, 11220) x_vmin_2d + dx, y_vmin_2d+3.0q0*dy, v(no_vmin_2d + 3*nx + 1), 47
	
	write (11450, 11220) x_vmin_2d + 2.0q0*dx, y_vmin_2d+3.0q0*dy, v(no_vmin_2d + 3*nx + 2), 48
	
	write (11450, 11220) x_vmin_2d + 3.0q0*dx, y_vmin_2d+3.0q0*dy, v(no_vmin_2d + 3*nx + 3), 49
	    
    11220 format (1pe48.32, 1pe48.32, 1pe48.32, 2x, 'v(x, y, v) do nó ', 1x, i3)		   
    
    close(11450)
    
 end subroutine esc_min_v_p6_2d
 
 !*************************************************************************************************************
 
 subroutine esc_max_v_p2_2d
 
 ! Salva dados para entrada em interp_2d, p = 2, max v
    
    350 format(A<len(trim(adjustl(caso)))>)
    
    write(caso4,350) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso)),"_max_v_p2_2d",".txt" 
    
    open(145,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(145,120) arquivo_caso
    
    120 format (3x, A40, 19x, 'nome do caso') 
        
    write (145, 121) 9
   
    121 format (1x, i6,48x, '= N, número de nós para interpolação, 2D, p = 2')
	
	write (145, 122) x_vmax_2d - dx, y_vmax_2d - dy, v(no_vmax_2d - nx - 1), 1
	
	write(145, 122)  x_vmax_2d, y_vmax_2d - dy, v(no_vmax_2d - nx), 2  
	
	write(145, 122) x_vmax_2d + dx, y_vmax_2d - dy, v(no_vmax_2d - nx + 1), 3
	
	write(145, 122) x_vmax_2d - dx, y_vmax_2d, v(no_vmax_2d - 1), 4
	
	write(145, 122) x_vmax_2d, y_vmax_2d, v(no_vmax_2d), 5

	write(145, 122) x_vmax_2d + dx, y_vmax_2d, v(no_vmax_2d + 1), 6

	write(145, 122) x_vmax_2d - dx, y_vmax_2d + dy, v(no_vmax_2d +nx - 1), 7 
	
	write(145, 122) x_vmax_2d, y_vmax_2d + dy, v(no_vmax_2d + nx), 8 
	
	write(145, 122) x_vmax_2d + dx, y_vmax_2d + dy, v(no_vmax_2d + nx + 1), 9
	    
    122 format (1pe48.32, 1pe48.32, 1pe48.32, 2x, '(x, y, v) do nó ', 1x, i3)		   
    
    close(145)
    
 end subroutine esc_max_v_p2_2d
 
 !****************************************************************************************
 
 subroutine esc_max_v_p4_2d
 
 ! Salva dados para entrada em interp_2d, p = 4, max v
    
    3500 format(A<len(trim(adjustl(caso)))>)
    
    write(caso,3500) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso)),"_max_v_p4_2d",".txt" 
    
    open(1450,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(1450,1200) arquivo_caso
    
    1200 format (3x, A40, 19x, 'nome do caso') 
        
    write (1450, 1210) 25
   
    1210 format (1x, i6,48x, '= N, número de nós para interpolação')
	
	write (1450, 1220) x_vmax_2d - 2.0q0*dx, y_vmax_2d - 2.0q0*dy, v(no_vmax_2d - 2*nx - 2), 1
	
	write(1450, 1220)  x_vmax_2d-dx, y_vmax_2d - 2.0q0*dy, v(no_vmax_2d - 2*nx - 1), 2  
	
	write(1450, 1220) x_vmax_2d, y_vmax_2d - 2.0q0*dy, v(no_vmax_2d - 2*nx), 3
	
	write(1450, 1220) x_vmax_2d + dx, y_vmax_2d - 2.0q0*dy, v(no_vmax_2d - 2*nx + 1), 4
	
	write(1450, 1220) x_vmax_2d + 2.0q0*dx, y_vmax_2d - 2.0q0*dy, v(no_vmax_2d - 2*nx + 2), 5

	write(1450, 1220) x_vmax_2d - 2.0q0*dx, y_vmax_2d - dy, v(no_vmax_2d - nx - 2), 6 

	write(1450, 1220) x_vmax_2d - dx, y_vmax_2d - dy, v(no_vmax_2d - nx - 1), 7 
	
	write(1450, 1220) x_vmax_2d, y_vmax_2d - dy, v(no_vmax_2d - nx), 8 
	
	write(1450, 1220) x_vmax_2d + dx, y_vmax_2d - dy, v(no_vmax_2d - nx + 1), 9
	
	write(1450, 1220) x_vmax_2d + 2.0q0*dx, y_vmax_2d - dy, v(no_vmax_2d - nx + 2), 10
	
	write(1450, 1220) x_vmax_2d - 2.0q0*dx, y_vmax_2d, v(no_vmax_2d - 2), 11
	
	write(1450, 1220) x_vmax_2d - dx, y_vmax_2d, v(no_vmax_2d - 1), 12
	
	write(1450, 1220) x_vmax_2d, y_vmax_2d, v(no_vmax_2d), 13
	
	write(1450, 1220) x_vmax_2d + dx, y_vmax_2d, v(no_vmax_2d + 1), 14
	
	write(1450, 1220) x_vmax_2d + 2.0q0*dx, y_vmax_2d, v(no_vmax_2d + 2), 15
	
	write(1450, 1220) x_vmax_2d - 2.0q0*dx, y_vmax_2d + dy, v(no_vmax_2d + nx - 2), 16
	
	write(1450, 1220) x_vmax_2d - dx, y_vmax_2d + dy, v(no_vmax_2d + nx - 1), 17
	
	write(1450, 1220) x_vmax_2d, y_vmax_2d + dy, v(no_vmax_2d + nx), 18
	
	write(1450, 1220) x_vmax_2d + dx, y_vmax_2d + dy, v(no_vmax_2d + nx + 1), 19
	
	write(1450, 1220) x_vmax_2d + 2.0q0*dx, y_vmax_2d + dy, v(no_vmax_2d + nx + 2), 20
	
	write(1450, 1220) x_vmax_2d - 2.0q0*dx, y_vmax_2d + 2.0q0*dy, v(no_vmax_2d + 2*nx - 2), 21
	
	write(1450, 1220) x_vmax_2d - dx, y_vmax_2d + 2.0q0*dy, v(no_vmax_2d + 2*nx - 1), 22
	
	write(1450, 1220) x_vmax_2d, y_vmax_2d + 2.0q0*dy, v(no_vmax_2d + 2*nx), 23
	
	write(1450, 1220) x_vmax_2d + dx, y_vmax_2d + 2.0q0*dy, v(no_vmax_2d + 2*nx + 1), 24
	
	write(1450, 1220) x_vmax_2d + 2.0q0*dx, y_vmax_2d + 2.0q0*dy, v(no_vmax_2d + 2*nx + 2), 25
	    
    1220 format (1pe48.32, 1pe48.32, 1pe48.32, 2x, 'v(x, y, v) do nó ', 1x, i3)		   
    
    close(1450)
    
 end subroutine esc_max_v_p4_2d
 
 !****************************************************************************************
 
 subroutine esc_max_v_p6_2d
 
 ! Salva dados para entrada em interp_2d, p = 6, max v
    
    33500 format(A<len(trim(adjustl(caso)))>)
    
    write(caso,33500) caso
    
    write(arquivo_caso, '(A,A,A)') trim(adjustl(caso)),"_max_v_p6_2d",".txt" 
    
    open(11450,file=arquivo_caso)          !Arquivo que contém os dados para interpolação
    
    write(11450,11200) arquivo_caso
    
    11200 format (3x, A40, 19x, 'nome do caso') 
        
    write (11450, 11210) 49
   
    11210 format (1x, i6,48x, '= N, número de nós para interpolação')
	
	write (11450, 11220) x_vmax_2d - 3.0q0*dx, y_vmax_2d - 3.0q0*dy, v(no_vmax_2d - 3*nx - 3), 1
	
	write (11450, 11220) x_vmax_2d - 2.0q0*dx, y_vmax_2d - 3.0q0*dy, v(no_vmax_2d - 3*nx - 2), 2
	
	write (11450, 11220) x_vmax_2d - dx, y_vmax_2d - 3.0q0*dy, v(no_vmax_2d - 3*nx - 1), 3
	
	write (11450, 11220) x_vmax_2d, y_vmax_2d - 3.0q0*dy, v(no_vmax_2d - 3*nx), 4
	
	write (11450, 11220) x_vmax_2d + dx, y_vmax_2d - 3.0q0*dy, v(no_vmax_2d - 3*nx +1), 5 
	
	write (11450, 11220) x_vmax_2d + 2.0q0*dx, y_vmax_2d - 3.0q0*dy, v(no_vmax_2d - 3*nx +2), 6 
	
	write (11450, 11220) x_vmax_2d + 3.0q0*dx, y_vmax_2d - 3.0q0*dy, v(no_vmax_2d - 3*nx +3), 7 
	
    write (11450, 11220) x_vmax_2d - 3.0q0*dx, y_vmax_2d - 2.0q0*dy, v(no_vmax_2d - 2*nx - 3), 8
	
	write (11450, 11220) x_vmax_2d - 2.0q0*dx, y_vmax_2d - 2.0q0*dy, v(no_vmax_2d - 2*nx - 2), 9
	
	write (11450, 11220) x_vmax_2d - dx, y_vmax_2d - 2.0q0*dy, v(no_vmax_2d - 2*nx - 1), 10
	
	write (11450, 11220) x_vmax_2d, y_vmax_2d - 2.0q0*dy, v(no_vmax_2d - 2*nx), 11
	
	write (11450, 11220) x_vmax_2d + dx, y_vmax_2d - 2.0q0*dy, v(no_vmax_2d - 2*nx +1), 12 
	
	write (11450, 11220) x_vmax_2d + 2.0q0*dx, y_vmax_2d - 2.0q0*dy, v(no_vmax_2d - 2*nx +2), 13 
	
	write (11450, 11220) x_vmax_2d + 3.0q0*dx, y_vmax_2d - 2.0q0*dy, v(no_vmax_2d - 2*nx +3), 14 
	
	write (11450, 11220) x_vmax_2d - 3.0q0*dx, y_vmax_2d - dy, v(no_vmax_2d - nx - 3), 15
	
	write (11450, 11220) x_vmax_2d - 2.0q0*dx, y_vmax_2d - dy, v(no_vmax_2d - nx - 2), 16
	
	write (11450, 11220) x_vmax_2d - dx, y_vmax_2d - dy, v(no_vmax_2d - nx - 1), 17
	
	write (11450, 11220) x_vmax_2d, y_vmax_2d - dy, v(no_vmax_2d - nx), 18
	
	write (11450, 11220) x_vmax_2d + dx, y_vmax_2d - dy, v(no_vmax_2d - nx +1), 19 
	
	write (11450, 11220) x_vmax_2d + 2.0q0*dx, y_vmax_2d - dy, v(no_vmax_2d - nx +2), 20 
	
	write (11450, 11220) x_vmax_2d + 3.0q0*dx, y_vmax_2d - dy, v(no_vmax_2d - nx +3), 21
	
	write (11450, 11220) x_vmax_2d - 3.0q0*dx, y_vmax_2d, v(no_vmax_2d  - 3), 22
	
	write (11450, 11220) x_vmax_2d - 2.0q0*dx, y_vmax_2d, v(no_vmax_2d  - 2), 23
	
	write (11450, 11220) x_vmax_2d - dx, y_vmax_2d, v(no_vmax_2d - 1), 24
	
	write (11450, 11220) x_vmax_2d, y_vmax_2d, v(no_vmax_2d), 25
	
	write (11450, 11220) x_vmax_2d + dx, y_vmax_2d, v(no_vmax_2d + 1), 26 
	
	write (11450, 11220) x_vmax_2d + 2.0q0*dx, y_vmax_2d, v(no_vmax_2d + 2), 27 
	
	write (11450, 11220) x_vmax_2d + 3.0q0*dx, y_vmax_2d, v(no_vmax_2d + 3), 28
	
	write (11450, 11220) x_vmax_2d - 3.0q0*dx, y_vmax_2d+dy, v(no_vmax_2d + nx - 3), 29
	
	write (11450, 11220) x_vmax_2d - 2.0q0*dx, y_vmax_2d+dy, v(no_vmax_2d + nx - 2), 30
	
	write (11450, 11220) x_vmax_2d - dx, y_vmax_2d+dy, v(no_vmax_2d + nx - 1), 31
	
	write (11450, 11220) x_vmax_2d, y_vmax_2d+dy, v(no_vmax_2d + nx), 32
	
	write (11450, 11220) x_vmax_2d + dx, y_vmax_2d+dy, v(no_vmax_2d + nx + 1), 33
	
	write (11450, 11220) x_vmax_2d + 2.0q0*dx, y_vmax_2d+dy, v(no_vmax_2d + nx + 2), 34
	
	write (11450, 11220) x_vmax_2d + 3.0q0*dx, y_vmax_2d+dy, v(no_vmax_2d + nx + 3), 35
	
	write (11450, 11220) x_vmax_2d - 3.0q0*dx, y_vmax_2d+2.0q0*dy, v(no_vmax_2d + 2*nx - 3), 36
	
	write (11450, 11220) x_vmax_2d - 2.0q0*dx, y_vmax_2d+2.0q0*dy, v(no_vmax_2d + 2*nx - 2), 37
	
	write (11450, 11220) x_vmax_2d - dx, y_vmax_2d+2.0q0*dy, v(no_vmax_2d + 2*nx - 1), 38
	
	write (11450, 11220) x_vmax_2d, y_vmax_2d+2.0q0*dy, v(no_vmax_2d + 2*nx), 39
	
	write (11450, 11220) x_vmax_2d + dx, y_vmax_2d+2.0q0*dy, v(no_vmax_2d + 2*nx + 1), 40
	
	write (11450, 11220) x_vmax_2d + 2.0q0*dx, y_vmax_2d+2.0q0*dy, v(no_vmax_2d + 2*nx + 2), 41
	
	write (11450, 11220) x_vmax_2d + 3.0q0*dx, y_vmax_2d+2.0q0*dy, v(no_vmax_2d + 2*nx + 3), 42
	
	write (11450, 11220) x_vmax_2d - 3.0q0*dx, y_vmax_2d+3.0q0*dy, v(no_vmax_2d + 3*nx - 3), 43
	
	write (11450, 11220) x_vmax_2d - 2.0q0*dx, y_vmax_2d+3.0q0*dy, v(no_vmax_2d + 3*nx - 2), 44
	
	write (11450, 11220) x_vmax_2d - dx, y_vmax_2d+3.0q0*dy, v(no_vmax_2d + 3*nx - 1), 45
	
	write (11450, 11220) x_vmax_2d, y_vmax_2d+3.0q0*dy, v(no_vmax_2d + 3*nx), 46
	
	write (11450, 11220) x_vmax_2d + dx, y_vmax_2d+3.0q0*dy, v(no_vmax_2d + 3*nx + 1), 47
	
	write (11450, 11220) x_vmax_2d + 2.0q0*dx, y_vmax_2d+3.0q0*dy, v(no_vmax_2d + 3*nx + 2), 48
	
	write (11450, 11220) x_vmax_2d + 3.0q0*dx, y_vmax_2d+3.0q0*dy, v(no_vmax_2d + 3*nx + 3), 49
	    
    11220 format (1pe48.32, 1pe48.32, 1pe48.32, 2x, 'v(x, y, v) do nó ', 1x, i3)		   
    
    close(11450)
    
 end subroutine esc_max_v_p6_2d


end module resultados