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
  !    ckx     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
  !    iopd    initial oil pressure distribution, psi
  !
  real(dp) :: nt,ni,np,npc,visorg,dbdp,delt1,eps,dez,&
       dey,fpb,pi,ooip,qotot,ctim,eot,borig,porig,&
       delt,convf,en,to,dpmx,dpp,delmin,delmax,dpmax,&
	   dtmult,delp,om,rhostd,dhdx,rhoim,rho,rhoi,tqchg,pconst,&
	   rw,rp,rm,pwlim,pwell
  !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
  !	   icrn    specific solution method
  integer  :: stmax,mx,isol,ipe,iskip,iwrite,i,im,ip,& 
       kcy,l,k,lcysw,ii,icrn,int

  open (unit=5,file='res90dat.dat', status='old')
  open (unit=6,file='res90out.dat', status='replace')
  open (unit=7,file='plott.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),pest(mx),&
	   dex(mx),kx(mx),phi(mx),ckx(mx),volphi(mx),a2(mx),vol(mx),b(mx),&
	   gho(mx),ghomin(mx),ghoplus(mx))


  read  (5,*)		porig, borig, pwlim
  read  (5,*)		visorg, dbdp, delt1, eps, dez, rw
  read  (5,*)		dey,rhostd,dhdx
  read  (5,*)		(qo(i),i = 1,mx)
  read  (5,*)		delmin, delmax, dpmax, dtmult
  read  (5,*)		icrn, om, int
  read  (5,*)		(dex(i),i = 1,mx)
  read  (5,*)		(kx(i),i = 1,mx)
  read  (5,*)		(phi(i),i = 1,mx)
  read  (5,*)		tqchg,pconst
  
  write (6,3010)	mx, stmax, isol, ipe, iskip, iwrite, icrn
  
  write (6,3020)	porig, borig, pwlim, visorg,&
                    dbdp, delt1, eps, dez, dey, rw, delmin, delmax,&
					dpmax, dtmult,om,rhostd,dhdx,tqchg
  
  write (6,3030)	(i,qo(i), i = 1,mx)
  write (6,3080)
  write (6,3034)    (i,dex(i),i = 1,mx)
  write (6,3080)
  write (6,3035)	(i,kx(i), i = 1,mx)
  write (6,3080)
  write (6,3036)	(i,phi(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.

	if (rw > 0.d0) then
		rp = rw
		do i = 1,mx
			rm = rp
			rp = rm + dex(i)
			x(i) = dsqrt(rm*rp)
			phi(i) = pi*(rp+rm)/dez*phi(i)
			kx(i) = 2.d0*pi*dex(i)/dez/dlog(rp/rm)*kx(i)
		end do
	end if
  
			im=1
	do i=1,mx
	if (i > 1) im=i-1
		vol(i)=	dez*dey*dex(i)/fpb
		volphi(i)=vol(i)*phi(i)
		ckx(i)=	2.d0*convf/visorg/(dex(i)/kx(i)+dex(im)/kx(im))/dex(i)
		a2(i)	=	phi(i)*dbdp
	end do  
		ckx(1) = 2.d0*convf/visorg*kx(1)/dex(1)**2
  
		ooip	=	0.0
		x(1)	=	dex(1)/to
	
	
	
	pressure_points:&	
    do i=2,mx
		im = i -1
		x(i) = dex(i) + 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

	if (dabs(dhdx) > 1.D-06) then
		do k=1,4
			call flprop(mx,porig,dbdp,borig,rhostd,dhdx)
		end do
		do i=2,mx
		po(i)=po(im)-gho(i)
		end do
	end if

	ooip_in_place:&
	do i = 1,mx
		ooip = ooip + volphi(i) * borig
	end do ooip_in_place

write (6,3031) (ooip)
write (6,3080)
write (6,3032) (i,po(i),i=1,mx)
write (6,3080)
		delt=delmin
		


time_step_loop:&
    do l = 1, stmax
		
	if (int .eq. 1) then
			int = 0
		read(5,*) qo(1),qo(mx),delmin,delmax,dtmult,&
				dpmax,pconst,tqchg
		write(6,3080)
		write(6,3037) qo(1),qo(mx),delmin,delmax,dtmult,&
			    dpmax,pconst,tqchg
		write(6,3080)

		if(tqchg < 0.d0) stop
		

		qo(1) = qo(1)/vol(1)
		qo(mx) = qo(mx)/vol(mx)
		a9(1) = oxmin(1)*pconst
		a9(mx) = oxmin(mx)*pconst
		delt = delmin		
	end if


		dpmx=0.d0
		pold = po
		lcysw	=	0

		delt=delt*dtmult
		delp=dpmx*dtmult
	
	if(delp.gt.dpmax) then
		delt=delt*0.8*dpmax/delp
	end if

	if(delt.lt.delmin) delt=delmin
	if(delt.gt.delmax) delt=delmax
		if(ctim + delt > tqchg) then
		int = 1
		delt = tqchg - ctim
		end if


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

	7744 continue

		pest = po
		kcy = kcy + 1

		po = (po + pold)/to

	call flprop(mx,porig,dbdp,borig,rhostd,dhdx)
	call flocon(mx,delt,icrn,dhdx,pconst,rw)

select case(isol)

	case(1)
		call tridia(mx)

	case(2)
		po = to*po - pold
		comp = po
		call sor(mx,eps,om)

	case(3)
		im=1
		do i=1,mx
			if (i .gt. 1) im = i-1
			if (i .lt. mx) ip = i+1
			po(i)=a(i)*pold(im)+c(i)*pold(ip)+d(i)
		end do
end select
!*******************************************************************
	dpmx=dabs(po(1)-pold(1))
	ii=1
!
!	største trykkendring og tilhørende blokk
!
	do i=2,mx
		dpp=dabs(pold(i)-po(i))
		if(dpp.gt.dpmx) then
		dpmx=dpp
		ii=i
	end if
	end do

!
!	sjekk om største verdi er over maks. tillatte .
!	tidsstegslengden reduseres og melding skrives
!	og iterasjoner i ytre sløyfe startes på ny.
!
	if(dpmx.gt.dpmax) then
		
		delt=delt*0.8*dpmax/dpmx

	write(6,3033) ii, dpmx
	
		
		lcysw=0
		
		goto 7744
	
	end if


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

		lcysw = lcysw + 1        

		dpmx=dabs(po(1)-pest(1))												
!
!	største trykkendring og tilhørende blokk
!
	do i = 2, mx
		dpp=dabs(pest(i)-po(i))
		if(dpp.gt.dpmx) then
			dpmx=dpp
		end if
	end do

check_convergence_of_nonlinear_term:&        
	if (dpmx > eps .and. lcysw < 4) then
		goto 7744     
	
end if check_convergence_of_nonlinear_term

	if (rw>0 .and. pconst>0) qo(1) = ckx(1)*bo(1)*(pconst - po(1))

	
    if (rw>0 .and. pconst<=0.d0) 	pwell = qo(1)/(ckx(1)*bo(1))+po(1)   
	
	if (pconst <= 0 .and. pwell < pwlim)  pconst = pwlim
	
	
	write (6,3080)
	write (6,3038) pwell
	write (6,3080)

!       start material balance
!         
		ctim = ctim + delt
		nt = 0.
		qotot = 0.      

	call flprop(mx,porig,dbdp,borig,rhostd,dhdx)


		matbal: do i = 1,mx
		nt = nt + volphi(i) * bo(i)
	if (qo(i) .gt. 0.) ni = ni+qo(i)*vol(i)*delt
	if (qo(i) .lt. 0.) np = np-qo(i)*vol(i)*delt
	if (qo(i) .lt. 0.) qotot = qotot-qo(i)*vol(i)
	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)
		write (6,3080)
	end if

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

!write (6,"('  kkkcum .........', i4/)") kkkcum
if (l==1) write (7,"(t14,'Tid',t34,'po(1)',t54,'pwell')")
write (7,*) ctim, po(1), pwell
end do time_step_loop

!   
!     ** close files
!   
close (5)
close (6)
!   
!     ** format statements
!   
3000 format (&
'     ***** One Dimensional, One Phase Simulator *****    '/&
'                     by Livar Lie                        '/&
'          ------------------------------------           '//&
'     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  /&
'  Specific solution method, icrn ............... :', i3  //)
!   
3020 format (&
'  Original oil pressure, porig, psia ........... :', f10.0/&
'  Original oil fvf, borig, std vol / res vol ... :', f10.4/&
'  Bottom hole limit pressure, pwlim, psia ...... :', f10.0/&
'  Original oil viscosity, visorg, cp ........... :', f10.4/&
'  Oil fvf pressure derivative, dbdp,              ',      /&
'    std vol/res vol/psi ........................ :', f11.8/&
'  Time increment, delt1, days .................. :', f9.6/&
'  Pressure tolerence, eps ...................... :', f8.5 /&
'  Block width, dez, ft ......................... :', f8.1 /&
'  Block thickness, dey, ft ..................... :', f8.1 /&
'  Borehole radius, rw, ft ...................... :', f8.4 /&
'  Minimum time step length, delmin ............. :', f12.9/&
'  Maximum time step length, delmax ............. :', f8.1 /&
'  Maximum pressure change pr. time step, dpmax . :', f8.1 /&
'  Multiplicator on time step length, dtmult .... :', f8.1 /&
'  Omega - relaxation factor .................... :', f8.1 /&
'  Trykkgradient, rhostd, psi/ft ................ :', f8.2 /&
'  Vertikalfaktor, sin alpha=dh/dx .............. :', f10.7 /&
'  Time step for flow change block 1, tqchg ..... :', f8.4 //)
!
!
3030 format (&
'  block no. :', i3, '     oil rate, qo, [stb/d] ................:', f10.1)
!
3031 format (&
'  Original oil i place, ooip, [stb] .......................:', g14.6)
!
3032 format (&
'  block no. :', i3, '     initial oil pressure, iop, [psia] ....:', f10.1)
!
3033 format (&
'  block no. :', i3, '     time step reduction , dpmx, ..........:', f8.3)
!
3034 format (&
'  block no. :', i3, '     delta length (x-dir), dex, [ft] ..... :', f10.1) 
!
3035 format (&
'  block no. :', i3, '     horizontal permeability, kx, [md]  .. :', f10.1)
!
3036 format (&
'  block no. :', i3, '     porosity, phi, fraction ............. :', f10.4)
!
3037 format (&
'  block no. 1 ', /&
'  oil rate first block, qo(1), [stb/d] ........................ :', f10.4/&
'  oil rate last block, qo(mx), [stb/d] ........................ :', f10.4/&
'  minimum time step length, delmin ............................ :', f12.9/&
'  maximum time step length, delmax ............................ :', f8.1 /&
'  multiplicator on time step length, dtmult ................... :', f8.1 /&  
'  minimum pressure change pr. time step, dpmax ................ :', f8.1 /&
'  constant pressure value, pconst, [psia] ..................... :', f8.1 /&
'  time step for flow change, tqchg, [d] ....................... :', f8.4 //)
!
3038 format (&
'  Well pressure, pwell, psia ..............................:', f10.4)
!
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 (16('*****')/)
!   
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 sor							 *
!     ****************************************************************
!   
      subroutine SOR(mx,eps,om)
        use felles_matriser
        implicit none

        integer :: kkk, i, isw, im, ip, mx, kkkcum

        real(dp) :: eps, om

        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) = po(i)+om*(a(i)*po(im)+c(i)*po(ip)+d(i)-po(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 subroutine sor
!   
!     ****************************************************************
!     *                     end subroutine gsit						 *
!     ****************************************************************
!   
!
!     ****************************************************************
!     *                     subroutine flprop						 *
!     ****************************************************************
!
      subroutine flprop(mx,porig,dbdp,borig,rhostd,dhdx)
        use felles_matriser
        implicit none
        integer :: mx, i,im,ip
        real(dp) :: porig,dbdp,borig,rhostd,dhdx,rhoim,rhoi

	    bo = borig + (po - porig)*dbdp

		im=1

		if (dabs(dhdx) > 1.D-06) then
			rhoim=rhostd*bo(1)
			gho(1)=rhoim*dex(1)*0.5D0*dhdx
			do i=2,mx
				im=i-1
				rhoi=rhostd*bo(i)
				gho(i)=(rhoi*dex(i)+rhoim*dex(im))*0.5D0*dhdx
				rhoim=rhoi
			end do
		end if

        return
      end subroutine flprop
!
!     ****************************************************************
!     *                     end subroutine flprop					 *
!     ****************************************************************
!
!     ****************************************************************
!     *                     subroutine flocon						 *
!     ****************************************************************
!
      subroutine flocon(mx,delt,icrn,dhdx,pconst,rw)
        use felles_matriser
        implicit none
        integer :: mx,i,im,ip,icrn
        real(dp) :: delt,dhdx,pconst,rw


		oxp:	do i = 2, mx
           im = i-1
           oxmin(i) = ckx(i)*(bo(i)+bo(im))/2.
           oxplus(im) = oxmin(i)*dex(i)/dex(im)
        end do oxp
        
	
		if (pconst < 1.d0) oxmin(1) =0.d0

		if(icrn.eq.1) then
		impl:   do i = 1, mx
		   if (i < mx) ip = i+1
		   a9(i)=qo(i)+gho(ip)*oxplus(i)-gho(i)*oxmin(i)
			if (rw > 0 .and. pconst > 0) then
		   oxmin(1) = ckx(1)*bo(1)
		   a9(1)=oxmin(1)*pconst
		end if
		   b(i) = oxmin(i)+oxplus(i)+a2(i)/delt
           a(i) = oxmin(i)/b(i)
           c(i) = oxplus(i)/b(i)
           d(i) =(a2(i)/delt*pold(i)+a9(i))/b(i)
	    end do impl

		else if(icrn.eq.2) then
		im=1
		cn:     do i=1,mx
		   if (i .gt. 1)  im = i-1
           if (i .lt. mx) ip = i+1
		   b(i)	= oxmin(i)+oxplus(i)+2*a2(i)/delt
		   a(i) = oxmin(i)/b(i)
		   c(i) = oxplus(i)/b(i)
		   d(i) = (oxplus(i)*(pold(ip)-pold(i))+oxmin(i)*(pold(im)-pold(i))&
				  +2*a9(i)+2*a2(i)/delt*pold(i))/b(i)
		   end do cn

		else if(icrn.eq.3) then
		ekspl:  do i=1,mx
		   b(i)	= delt/a2(i)
		   a(i) = b(i)*oxmin(i)
		   c(i) = b(i)*oxplus(i)
		   d(i) = b(i)*(a9(i)-(oxplus(i)+oxmin(i))*pold(i))+pold(i)
		   end do ekspl
        end if
			return
      end subroutine flocon
! 
!     ****************************************************************
!     *                     end subroutine flocon					 *
!     ****************************************************************
!

 