program utdelt_sim
	use felles_matriser
	implicit none
      
!list of real variables in double precision:
!    kx      permeabilitet i x-retning, md
!    nt      oljevolum tilstede ved et gitt tidspunkt, stb
!    ni      oljevolum injisert, stb
!    np      oljevolum produsert, stb
!    npc     oljevolum produsert i % 
!    visorg  opprinnelig viskositet, cp
!    dbdp    derivert av bo mhp trykk, stb/rb/psi
!    delt1   de 12 første tidssteg er delt1/(1.5**(12-l)), hvor l er antall
!            tidssteg, days
!    eps     trykktoleranse, iterasjon på ulineære ledd _og_ i iterativ
!            løsningsrutine, psi        
!    dex     delta x, ft
!    dey     delta y, ft
!    dez     delta z, ft
!    fpb     cubic feet per barrel, 5.61...
!    pi      3.14...
!    vol     dex*dey*dez/fpb, rb
!    volphi  vol*phi, rb
!    ckbak   konstantleddet i strømningskoeffisientene
!    a2      phi*dbdp
!    ooip    original oil in place, stb
!    pest    estimated value of po(1), psi
!    qotot   total oil rate, stb/day
!    ctim    kumulativ simulert tid, days
!    eot     total oil error
!    borig   original value of bo, stb/rb
!    phi     porosity
!    porig   original oil pressure, psi
!    delt    time step size, days
!    convf   conversion factor from Darcy units to oil-field units
!    en      1.d0
!    to      2.d0
!    

	real(dp) :: kx,nt,ni,np,npc,visorg,dbdp,delt1,eps,dez,&
				dey,fpb,pi,vol,volphi,ckbak,a2,ooip,&
				pest,qotot,ctim,eot,borig,phi,dex,porig,&
    			delt,convf,en,to

!list of integers
!    stmax   number of time steps
!    mx      number of numerical blocks
!    isol    specify solution method
!    ipe     trigger, not in use
!    iskip   trigger, not in use
!    iwrite  trigger, if equal to 1, intermediat pressures are written
!    i       loop counter
!    im      i-1
!    ip      i+1
!    lcysw   number of iterations per time step on non-linear terms
!    kcy     sum of lcysw
!    l       time step counter
!    k       loop counter
!

	integer  :: stmax,mx,isol,ipe,iskip,iwrite,i,im,ip,&
                kcy,l,k,lcysw

        open (unit=5,file='res90dat.dat', status='old')
        open (unit=6,file='res90out.dat', status='replace')
 
!     ** read and write input data.

        write (6,3000)
        read  (5,*)       mx, stmax, isol, ipe, iskip, iwrite

allocate (&
            a(mx),c(mx),d(mx),al(mx),s(mx),comp(mx),po(mx),&
            pold(mx),oxplus(mx),oxmin(mx),a9(mx),qo(mx),bo(mx),x(mx)&
         )

      write (6,3010)	mx, stmax, isol, ipe, iskip, iwrite
      read  (5,*)		kx, phi, dex, porig, borig
      read  (5,*)		visorg, dbdp, delt1, eps, dez
      read  (5,*)		dey
      write (6,3020)	kx, phi, dex, porig, borig, visorg, &
						dbdp, delt1, eps, dez, dey
      read (5,*)		(qo(i),i = 1,mx)

      write (6,3030)	(i,qo(i), i = 1,mx)
      write (6,3080)
!   
!     ** begin calculations; initialize variables
!   
      fpb	=	5.6146d0
      pi	=	3.1415926d0
	  convf	=	0.00632827

	  en	=	1.d0
	  to	=	2.d0

	  oxplus(mx)	= 0.
      oxmin(1)		= 0.
      delt          = 0.
      ni			= 0.
      np			= 0.
      ctim			= 0.
      kcy			= 0

      vol	=	dez*dey*dex/fpb
      volphi=	vol*phi
      ckbak =	convf*kx/dex/visorg/dex
      a2	=	phi*dbdp

      ooip	=	0.0
      x(1)	=	dex/to
pressure_points:&	
	do i	= 2,mx
		im = i -1
		x(i) = dex + x(im)
	end do pressure_points

      write (6,3040) (i,x(i),i=1,mx)
      write (6,3080)

	  po	= porig
	  qo	= qo /vol
	  a9	= qo
	  pold	= porig
	   
ooip_in_place:&
	do i = 1,mx
		ooip = ooip + volphi * borig
	end do ooip_in_place

time_step_loop:&
	do l = 1, stmax
		
		pold = po

		lcysw	=	0
		if (l <= 12) delt = delt1/1.5**(12-l)

		write (6,3070) 
!   
!   cc start iteration loop on unlinear variables i.e. bo
!  

7744	continue
 
		pest = po(1)
		kcy = kcy + 1

		po = (po + pold)/to

        call flprop(mx,porig,dbdp,borig)
		call flocon(mx,delt,a2,ckbak)


		select case(isol)

            case(1)
				call tridia(mx)

            case(2)
				po = to*po - pold
                comp = po
                call gsit(mx,eps)
        end select

        if (iwrite == 1) then
            write(6,3050) (i,po(i), i = 1,mx)
            write(6,3080)
        end if
        
        lcysw = lcysw + 1        

check_convergence_of_nonlinear_term:&        
        if (&
                dabs(po(1)-pest) > eps &
            .and. &
                lcysw < 4) then

					goto 7744     
					  
        end if check_convergence_of_nonlinear_term
        
!   
!       start material balance
!         
        ctim = ctim + delt

        nt = 0.
        qotot = 0.      

        call flprop(mx,porig,dbdp,borig)

      
matbal:&       
        do i = 1,mx
           nt = nt + volphi * bo(i)
           if (qo(i) .gt. 0.) ni = ni+qo(i)*vol*delt
           if (qo(i) .lt. 0.) np = np-qo(i)*vol*delt
           if (qo(i) .lt. 0.) qotot = qotot-qo(i)*vol
        end do matbal           

        eot = (-ni+np+nt-ooip)/(ooip+ni)*100.
        npc = np/ooip*100.

        if (iwrite /=  1) then 
           write (6,3050) (i,po(i),i=1,mx)
        end if

        write (6,3060)l,ctim,np,delt,kcy,ni,nt,qotot,eot,npc

	end do time_step_loop

!   
!     ** close files
!   
      close (5)
      close (6)
!   
!     ** format statements
!   
 3000 format (&
     '     ***** one dimensional, one phase simulator *****    '/&
     '           ------------------------------------          '//&
     '     basic assumptions of this model include:            '//&
     '      - horizontal flow of a compressible fluid          '/&
     '      - pressure dependent volume factor                 '/&
     '      - pressure independent viscosity                   '/&
     '      - homogeneous rock properties, that is the same    '/&
     '        porosity and permeability in all blocks.         '/&
     '      - flowing fluid is undersaturated oil              '//)
!   
 3010 format (&
     '  number of blocks in x-direction, mx .......... :', i3  /&
     '  maximum number of time steps, stmax .......... :', i3  /&
     '  solution method used, isol ................... :', i3  /&
     '  ipe .......................................... :', i3  /&
     '  iskip ........................................ :', i3  /&
     '  write option, iwrite ......................... :', i3  //)
!   
 3020 format (&
     '  horizontal permeability, kx, md  ............. :', f12.4/&
     '  porosity, phi, fraction ...................... :', f8.4 /&
     '  delta length (x-dir), dex, ft ................ :', f10.2/&
     '  original oil pressure, porig, psia ........... :', f10.1/&
     '  original oil fvf, borig, std vol / res vol ... :', f10.4/&
     '  original oil viscosity, visorg, cp ........... :', f10.4/&
     '  oil fvf pressure derivative, dbdp,              ',      /&
     '    std vol/res vol/psi ........................ :', f10.8/&
     '  time increment, delt1, days .................. :', f10.6/&
     '  pressure tolerence, eps ...................... :', f8.5 /&
     '  block width, dez, ft ......................... :', f8.1 /&
     '  block thickness, dey, ft ..................... :', f8.1 //)
!   
 3030 format (&
     '  block no. :', i3, '     oil rate, qo, stb/d .... ', f10.1)
!   
 3040 format (&
     '  block no. :', i3, '     pressure point, x(i), ft ', f8.1)
!   
 3050 format (&
     '  block no. :', i3, '     oil pressure, psia ..... ', f12.3)
!   
 3060 format (2x,'steps=',i4,t22,'time=',d12.4,t43,'np=',d12.4/2x,&
    'delt=',d12.4,t22,'cum cycle=',i4,t43,'ni=',d12.4/2x,&
    'nt=',d12.4,t22,'qotot=',d12.4,t43,'cum oil error=',d12.4/2x,&
    'npc=',d12.4)
!   
 3070 format (25('*****')/)
!   
 3080 format (/)
!      stop
      end
!   
!     ****************************************************************
!     *                  end of main program
!     ****************************************************************
!   
!   
!     ****************************************************************
!     *                  subroutine tridia
!     ****************************************************************
!   
subroutine tridia(mx)
  use felles_matriser
		implicit none

		integer :: i, mx, k

        al(1) = 1.
		s(1) = d(1)

upper_trian:	do i = 2,mx
                    al(i) =  1.   - a(i)/al(i-1)*c(i-1)
                    s(i) = d(i) + a(i)/al(i-1)*s(i-1)
				end do upper_trian

				po(mx) = s(mx)/al(mx)	!last pressure

back_subst:     do k = 2,mx
                    i = mx - k + 1
                    po(i) = (s(i)+c(i)*po(i+1)) / al(i)
                end do back_subst

      return
end subroutine tridia
!   
!     ****************************************************************
!     *                end of subroutine tridia
!     ****************************************************************
!   
!   
!     ****************************************************************
!     *                     subroutine gsit
!     ****************************************************************
!   
subroutine gsit(mx,eps)
	use felles_matriser
	implicit none

	integer :: kkk, i, isw, im, ip, mx, kkkcum
	real(dp) :: eps

itera:	do kkk = 1,300
			kkkcum = kkk
			isw = 0
			im = 1

pressures:	do i = 1,mx
				if (i .gt. 1) im = i-1
				if (i .lt. mx) ip = i+1

				po(i) = a(i)*po(im)+c(i)*po(ip)+d(i)

				if (isw == 0) then
					if (dabs(comp(i)-po(i)) > eps) isw = 1
				end if

				comp(i) = po(i)

			end do pressures

			if (isw .eq. 0) exit

		end do itera

		write (6,"('  kkkcum .........', i4/)") kkkcum

	return
end
!   
!     ****************************************************************
!     *                     end subroutine gsit
!     ****************************************************************
!   
!
!     ****************************************************************
!     *                     subroutine flprop
!     ****************************************************************
!
subroutine flprop(mx,porig,dbdp,borig)
	use felles_matriser
	implicit none
		integer :: mx, i
		real(dp) :: porig,dbdp,borig

		bo = borig + (po - porig)*dbdp

	return
end subroutine flprop
!
!     ****************************************************************
!     *                     end subroutine flprop
!     ****************************************************************
!
!     ****************************************************************
!     *                     subroutine flocon
!     ****************************************************************
!
subroutine flocon(mx,delt,a2,ckbak)
	use felles_matriser
	implicit none
		integer :: mx,i,im
		real(dp) :: delt,a2,ckbak,b
oxp:	do i = 2, mx
			im = i-1
            oxmin(i) = ckbak*(bo(i)+bo(im))/2.
            oxplus(im) = oxmin(i)
		end do oxp
acd:	do i = 1, mx
			b    = oxmin(i)+oxplus(i)+a2/delt
			a(i) = oxmin(i)/b
			c(i) = oxplus(i)/b
            d(i) =(a2/delt*pold(i)+a9(i))/b
		end do acd
return
end subroutine flocon
!
!     ****************************************************************
!     *                     end subroutine flocon
!     ****************************************************************
!


