From f2a35cdb45a6be05aaa9cfd98f1f122e8e56d2a3 Mon Sep 17 00:00:00 2001 From: Oscar Theran Date: Mon, 12 Aug 2024 19:00:36 -0300 Subject: [PATCH 01/16] 12/08/2024 update --- oil1.nml | 13 ++++++++++++- src/routines.f90 | 18 ++++++++++-------- src/types.f90 | 1 + 3 files changed, 23 insertions(+), 9 deletions(-) diff --git a/oil1.nml b/oil1.nml index b1e1ed0..5901db2 100644 --- a/oil1.nml +++ b/oil1.nml @@ -11,6 +11,10 @@ ! number of defined components before scn scn_nc = 14 ! number of single carbon number components + scn_nc_ps = 20 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped / &nml_components def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" @@ -37,4 +41,11 @@ ! molecular weights of single carbon numbers plus_mw = 470.7 ! molecular weight of residual fraction -/ \ No newline at end of file +/ +&nml_density + scn_density = 0.685 0.749 0.768 0.793 0.808 0.815 0.836 0.85 0.861 0.873 0.882 0.873 0.875 0.885 + ! density of single carbon numbers + plus_density = 0.936 + ! density of residual fraction +/ + diff --git a/src/routines.f90 b/src/routines.f90 index 445acec..31339a6 100644 --- a/src/routines.f90 +++ b/src/routines.f90 @@ -544,7 +544,8 @@ subroutine calculate_volume_6plus(fluid, mw_source, characterization) cn_plus => characterization%carbon_number_plus, & z_m_plus_i => characterization%product_z_mw_plus_i, & scn_z_i => characterization%scn_z, & - scn_mw_i => characterization%scn_mw & + scn_mw_i => characterization%scn_mw, & + scn_zm => characterization%scn_zm & ) scn_density_cal = ((a_d)*(exp(- (real(fluid%scn, pr))/10._pr))) + b_d @@ -553,12 +554,12 @@ subroutine calculate_volume_6plus(fluid, mw_source, characterization) select case (mw_source) case("experimental") - volume_6plus_cal = sum((scn_z_i*scn_mw_i)/(fluid%scn_density)) + & + volume_6plus_cal = sum((scn_zm)/(fluid%scn_density)) + & sum(z_m_plus_i/plus_density_cal) plus6_density = [fluid%scn_density, plus_density_cal] case("calculated") - volume_6plus_cal = sum((scn_z_i*scn_mw_i)/(scn_density_cal)) + & + volume_6plus_cal = sum((scn_zm)/(scn_density_cal)) + & sum(z_m_plus_i/plus_density_cal) plus6_density = [scn_density_cal, plus_density_cal] end select @@ -598,6 +599,7 @@ subroutine lump (fluid, characterization) associate (& plus_zm => characterization%plus_zm, & + plus_density => characterization%plus_density, & scn_zm => characterization%scn_zm, & plus_z => characterization%plus_z, & scn_z => characterization%scn_z, & @@ -704,10 +706,10 @@ subroutine lump (fluid, characterization) print*, 'ps',i, plus_z_ps(i), plus_mw_ps(i), density_ps(i), j_ps(i), carbon_number_plus(scn_nc_new+i) end do - end associate - - !characterization%product_z_mw_plus_i = z_m_plus_i - !characterization%plus_z_i = plus_z_i + plus_density = plus_zm / sum(plus_z_ps(1:numbers_ps)* & + plus_mw_ps(1:numbers_ps)/(density_ps(1:numbers_ps))) + print*, plus_density + characterization%last_C = last_C characterization%i_last = i_last @@ -716,7 +718,7 @@ subroutine lump (fluid, characterization) !do i = 1, scn_nc_new ! print*, prev_i + i , characterization%scn_z(i) , characterization%scn_mw(i) !end do - + end associate end subroutine lump diff --git a/src/types.f90 b/src/types.f90 index afec0a3..e00cde9 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -63,6 +63,7 @@ module dtypes real(pr), allocatable :: scn_i(:) real(pr), allocatable :: plus6_density(:) real(pr) :: plus_w + real(pr) :: plus_density end type FluidDataOut end module dtypes From afa98e7657bb65de453284848b399a075b827a85 Mon Sep 17 00:00:00 2001 From: Oscar Theran Date: Fri, 16 Aug 2024 17:30:18 -0300 Subject: [PATCH 02/16] changes until 16/08/2024 --- app/main.f90 | 2 +- oil1.nml | 2 +- src/routines.f90 | 43 +++++++++++++++++++++++++++---------------- src/types.f90 | 1 + 4 files changed, 30 insertions(+), 18 deletions(-) diff --git a/app/main.f90 b/app/main.f90 index 170cbc0..b9efea0 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -9,7 +9,7 @@ program main type(FluidData) :: fluid - prueba = characterize(file='oil1.nml', mw_source="calculated", method = "plus_mw", fix_C=.true.) + prueba = characterize(file='oil1.nml', mw_source="calculated", method = "plus_mw", fix_C=.true., eos='eos') print*, prueba%n_init print*, prueba%a , prueba%b print*, prueba%C, prueba%plus_mw diff --git a/oil1.nml b/oil1.nml index 5901db2..081200a 100644 --- a/oil1.nml +++ b/oil1.nml @@ -11,7 +11,7 @@ ! number of defined components before scn scn_nc = 14 ! number of single carbon number components - scn_nc_ps = 20 + scn_nc_ps = 10 ! CN from which all SCN fractions will be lumped into the specified number of pseudos numbers_ps = 5 ! number of pseudos in which the scn fractions grouped diff --git a/src/routines.f90 b/src/routines.f90 index 31339a6..82082b3 100644 --- a/src/routines.f90 +++ b/src/routines.f90 @@ -566,7 +566,7 @@ subroutine calculate_volume_6plus(fluid, mw_source, characterization) characterization%volume_6plus_cal = volume_6plus_cal characterization%plus6_density = plus6_density - + end associate end subroutine calculate_volume_6plus @@ -596,8 +596,11 @@ subroutine lump (fluid, characterization) real(pr), allocatable :: density_ps(:) real(pr), allocatable :: w_ps(:) real(pr) :: sum_zm_last_ps + real(pr), allocatable :: moles(:) associate (& + def_nc => fluid%def_comp_nc , & + scn_mw => characterization%scn_mw, & plus_zm => characterization%plus_zm, & plus_density => characterization%plus_density, & scn_zm => characterization%scn_zm, & @@ -622,6 +625,7 @@ subroutine lump (fluid, characterization) allocate (plus_mw_ps(0)) allocate(density_ps(0)) allocate (w_ps(0)) + allocate(moles(0)) scn_nc_input = fluid%scn_nc scn_nc_new = fluid%scn_nc_ps - 6 @@ -656,7 +660,7 @@ subroutine lump (fluid, characterization) ! Lumping into Nps pseudos rec_zm = plus_zm / numbers_ps ! Recommended value for the product z*M (proportional to weight) - !for each pseudo according to pedersen + ! for each pseudo according to pedersen remain_plus_zm = plus_zm i_ps = 1._pr sum_z = 0.0_pr @@ -692,7 +696,6 @@ subroutine lump (fluid, characterization) ! at this point, Nps is the order for the last NON-Asphaltenic pseudo comp. ! (e.g. 4 if 5 is for Asphaltenes) plus_mw_ps =[plus_mw_ps, sum(z_m_plus_i(i+1:i_last))/plus_z_ps(numbers_ps)] - ! (zMp - sum(zMpi(1:i))) / zps(Nps) w_ps = [w_ps, characterization%plus_w * ((plus_z_ps(numbers_ps))* & (plus_mw_ps(numbers_ps))/(plus_zm))] carbon_number_plus(scn_nc_new + numbers_ps) = 6 + & @@ -708,28 +711,34 @@ subroutine lump (fluid, characterization) plus_density = plus_zm / sum(plus_z_ps(1:numbers_ps)* & plus_mw_ps(1:numbers_ps)/(density_ps(1:numbers_ps))) - print*, plus_density - + + moles = [moles, (fluid%def_comp_w)/(fluid%def_comp_mw)] + moles = [moles, fluid%w(def_nc+1:def_nc+scn_nc_new)/ scn_mw(1:scn_nc_new)] + moles = [moles, w_ps/plus_mw_ps ] + characterization%mol_fraction = moles/sum(moles) ! mole fractions normalize + characterization%last_C = last_C characterization%i_last = i_last - - !! esto es lo que haria en el write, luego eliminar este comentario. - !prev_i = fluid%scn(1)-1 - !do i = 1, scn_nc_new - ! print*, prev_i + i , characterization%scn_z(i) , characterization%scn_mw(i) - !end do + end associate end subroutine lump - subroutine get_critical_constants(fluid, characterization) - !! this subroutine .... + subroutine get_critical_constants(fluid, characterization, eos) + !! this subroutine doing: + !! - Calculation of Tc, Pc, omega from Pedersen correlations (adapted from code "C7plusPedersenTcPcOm") + !! - Calculation of EoS parameters (taken from CubicParam) + + implicit none type(FluidData), intent(in) :: fluid type(FluidDataOut), intent(inout) :: characterization - !implicit none - + character(len=*), intent(in) :: eos + + !select case () + + !end select @@ -737,7 +746,7 @@ subroutine get_critical_constants(fluid, characterization) end subroutine get_critical_constants - type(FluidDataOut) function characterize(file, mw_source, method, fix_C) & + type(FluidDataOut) function characterize(file, mw_source, method, fix_C, eos) & result(characterization) type(FluidData) :: fluid @@ -745,6 +754,7 @@ type(FluidDataOut) function characterize(file, mw_source, method, fix_C) & character(len=*), intent(in), optional :: method !! plus_mw or global_mw character(len=*), intent(in) :: mw_source logical, intent(in), optional :: fix_C + character(len=*), intent(in) :: eos fluid = data_from_file(file) @@ -757,6 +767,7 @@ type(FluidDataOut) function characterize(file, mw_source, method, fix_C) & allocate(characterization%scn_zm(fluid%scn_nc)) allocate(characterization%scn_i(0)) allocate(characterization%plus6_density(0)) + allocate(characterization%mol_fraction(0)) call get_c_or_m_plus(fluid=fluid, mw_source=mw_source, method=method, fix_C=fix_C, characterization= characterization) diff --git a/src/types.f90 b/src/types.f90 index e00cde9..3571c42 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -64,6 +64,7 @@ module dtypes real(pr), allocatable :: plus6_density(:) real(pr) :: plus_w real(pr) :: plus_density + real(pr), allocatable :: mol_fraction(:) end type FluidDataOut end module dtypes From 7ee745f56ba45120fdc8c39f3cd081e199ae9d21 Mon Sep 17 00:00:00 2001 From: Oscar Theran Date: Thu, 22 Aug 2024 11:57:00 -0300 Subject: [PATCH 03/16] en desarrollo constantes criticas 22/08/2024 --- src/routines.f90 | 70 ++++++++++++++++++++++++++++++++++++++++++++++-- src/types.f90 | 2 ++ 2 files changed, 70 insertions(+), 2 deletions(-) diff --git a/src/routines.f90 b/src/routines.f90 index 82082b3..24e8eb9 100644 --- a/src/routines.f90 +++ b/src/routines.f90 @@ -734,11 +734,77 @@ subroutine get_critical_constants(fluid, characterization, eos) type(FluidDataOut), intent(inout) :: characterization character(len=*), intent(in) :: eos + select case (eos) + + case("srk") + ! Coeficientes originales Tabla 5.3 Pedersen + c1 = 1.6312e2_pr + c2 = 8.6052d1 + c3 = 4.3475d-1 + c4 = -1.8774d3 + d1P = -1.3408d-1 + d2 = 2.5019 + d3 = 2.0846d2 + d4 = -3.9872d3 + d5 = 1.0d0 + e1 = 7.4310d-1 + e2 = 4.8122d-3 + e3 = 9.6707d-3 + e4 = -3.7184d-6 + A = -0.176 + B = 1.574 + C0 = 0.48 + del1 = 1.0D0 + case("pr") + ! Coeficientes originales Tabla 5.3 Pedersen + c1 = 7.34043d1 + c2 = 9.73562d1 + c3 = 6.18744d-1 + c4 = -2.05932d3 + d1P = 7.28462d-2 + d2 = 2.18811d0 + d3 = 1.63910d2 + d4 = -4.04323d3 + d5 = 0.25d0 + e1 = 3.73765d-1 + e2 = 5.49269d-3 + e3 = 1.17934d-2 + e4 = -4.93049d-6 + A = -0.26992 + B = 1.54226 + C0 = 0.37464 + del1 = 1.0D0 + sqrt(2.0) + case("rkpr") + c1 = 7.34043d1 + c2 = 9.73562d1 + c3 = 6.18744d-1 + c4 = -2.05932d3 + d1P = 2.3d0 + d2 = -0.5d0 + d3 = 1.85d2 + d4 = -4.0d3 + d5 = 0.25d0 + e1 = 3.73765d-1 + e2 = 5.49269d-3 + e3 = 1.17934d-2 + e4 = -4.93049d-6 + A = -0.26992 + B = 1.54226 + C0 = 0.37464 + end select + + if (eos == "pr" .or. eos=="srk") + + d1 = (1 + del1**2)/(1 + del1) + y = 1 + (2*(1 + del1))**(1.0d0/3) + (4/(1 + del1))**(1.0d0/3) + OMa = (3*y*y + 3*y*d1 + d1**2 + d1 - 1.0d0)/(3*y + d1 - 1.0d0)**2 + OMb = 1/(3*y + d1 - 1.0d0) + Zc = y/(3*y + d1 - 1.0d0) + endif - !select case () - !end select + diff --git a/src/types.f90 b/src/types.f90 index 3571c42..de5820e 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -67,6 +67,8 @@ module dtypes real(pr), allocatable :: mol_fraction(:) end type FluidDataOut + + end module dtypes From cdfb48fe369037627639e2523971ae6534dac766 Mon Sep 17 00:00:00 2001 From: Oscar Theran Date: Mon, 2 Sep 2024 10:47:01 -0300 Subject: [PATCH 04/16] version 02/09/2024 --- app/main.f90 | 6 +- oil1.nml | 2 +- src/critical_parameters.f90 | 88 ++++ src/module.f90 | 836 ------------------------------------ src/prueba.f90 | 466 -------------------- src/routines.f90 | 155 ++++--- src/types.f90 | 8 + 7 files changed, 175 insertions(+), 1386 deletions(-) create mode 100644 src/critical_parameters.f90 delete mode 100644 src/module.f90 delete mode 100644 src/prueba.f90 diff --git a/app/main.f90 b/app/main.f90 index b9efea0..d85eaf3 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -4,18 +4,18 @@ program main use routines implicit none - integer :: k, i type(FluidDataOut) :: prueba type(FluidData) :: fluid - prueba = characterize(file='oil1.nml', mw_source="calculated", method = "plus_mw", fix_C=.true., eos='eos') + prueba = characterize(file='oil1.nml', mw_source="calculated", method = "plus_mw", fix_C=.true., eos='PR') print*, prueba%n_init print*, prueba%a , prueba%b print*, prueba%C, prueba%plus_mw print*, prueba%a_d, prueba%b_d print*, "-------------------------------------------------------------------" - print*, prueba%nc_plus + + end program main diff --git a/oil1.nml b/oil1.nml index 081200a..5901db2 100644 --- a/oil1.nml +++ b/oil1.nml @@ -11,7 +11,7 @@ ! number of defined components before scn scn_nc = 14 ! number of single carbon number components - scn_nc_ps = 10 + scn_nc_ps = 20 ! CN from which all SCN fractions will be lumped into the specified number of pseudos numbers_ps = 5 ! number of pseudos in which the scn fractions grouped diff --git a/src/critical_parameters.f90 b/src/critical_parameters.f90 new file mode 100644 index 0000000..e8b2377 --- /dev/null +++ b/src/critical_parameters.f90 @@ -0,0 +1,88 @@ +module critical_parameters + use constants + use dtypes, only: FluidData + + implicit none + + real(pr), allocatable :: tc_def(:), pc_def(:), om_def(:) + real(pr):: c1, c2, c3, c4, d1_p, d2, d3, d4, d5, e1, e2, e3, e4 + real(pr):: a, b, c0, del_1 + +contains + + subroutine get_parameteres_for_critical(eos) + implicit none + character(len=*) :: eos + + tc_def = [126.2, 304.21, 190.564, 305.32, 369.83, 408.14, 425.12, 460.43, 469.7] + pc_def = [34.0, 73.83, 45.99, 48.72, 42.48, 36.48, 37.96, 33.81, 33.7] + om_def = [0.038, 0.224, 0.012, 0.099, 0.152, 0.181, 0.20, 0.228, 0.252] + + select case(eos) + + case("SRK") + c1 = 1.6312d2 + c2 = 8.6052d1 + c3 = 4.3475d-1 + c4 = -1.8774d3 + d1_p = -1.3408d-1 + d2 = 2.5019 + d3 = 2.0846d2 + d4 = -3.9872d3 + d5 = 1.0d0 + e1 = 7.4310d-1 + e2 = 4.8122d-3 + e3 = 9.6707d-3 + e4 = -3.7184d-6 + a = -0.176 + b = 1.574 + c0 = 0.48 + !del_1 = 1.0D0 + + case("PR") ! Coeficientes originales Tabla 5.3 Pedersen + c1 = 7.34043d1 + c2 = 9.73562d1 + c3 = 6.18744d-1 + c4 = -2.05932d3 + d1_p = 7.28462d-2 + d2 = 2.18811d0 + d3 = 1.63910d2 + d4 = -4.04323d3 + d5 = 0.25d0 + e1 = 3.73765d-1 + e2 = 5.49269d-3 + e3 = 1.17934d-2 + e4 = -4.93049d-6 + a = -0.26992 + b = 1.54226 + c0 = 0.37464 + !del_1 = 1.0D0 + sqrt(2.0) + + case("RKPR") + c1 = 7.34043d1 ! Coeficientes originales Tabla 5.3 Pedersen + c2 = 9.73562d1 + c3 = 6.18744d-1 + c4 = -2.05932d3 + d1_p = 2.3d0 + d2 = -0.5d0 + d3 = 1.85d2 + d4 = -4.0d3 + d5 = 0.25d0 + e1 = 3.73765d-1 + e2 = 5.49269d-3 + e3 = 1.17934d-2 + e4 = -4.93049d-6 + a = -0.26992 + b = 1.54226 + c0 = 0.37464 + + end select + + end subroutine get_parameteres_for_critical + + end module + + + + + diff --git a/src/module.f90 b/src/module.f90 deleted file mode 100644 index 19e8c82..0000000 --- a/src/module.f90 +++ /dev/null @@ -1,836 +0,0 @@ - -module data - - use constants, only: pr - - !use types, only: FluidData - - implicit none - integer :: nv, nv1, nNcut, nNdef - real(pr) :: Zplus - - -contains - - subroutine data_in (Nfluid,Oil,Ncut,Ndef,DefComp,zdef,rMdef,rCN,zcomp,rMW,& - Den,Plus, rMWplus, DenPlus,Nps,wat,watPlus,zM,zMp,Z6p,sumV,w, rn ) - - implicit none - - integer, parameter :: maxD=15, imax=48 - integer :: i,j!, nv1 - integer :: file_unit_1, file_unit_2, file_unit_3, Nfluid, Ncut, Ndef, Nps - integer, dimension(330) :: rCN - real(pr), dimension(maxD) :: zdef, rMdef - real(pr), dimension(imax) :: zcomp, rMW, Den, wat, zM, zMdef, w, rn,MWt - real(pr) :: rMWplus, DenPlus , watPlus, zMp, Z6p, sumV, zMtotal - character(len=30) :: infile, outfile - character(len=7) :: Oil, Plus - character(len=4) :: DefComp(maxD) - - write(*,*) 'ARE MOLECULAR WEIGTHS FROM DATAIN EXPERIMENTAL?' - write(*,*) '1: YES 2: NO' - read(*,*) nv - - write(*,*) 'ENTER THE VERSION NUMBER' - write(*,*) '1: C.M 2: CMplus' - read(*,*) nv1 - - !nv2 = nv1 - - ! open input file - write(*,*) 'ENTER INFILE' - read(*,*) infile - write(*, *) "READING INFILE: " // infile - open(newunit=file_unit_1, file =infile) - write(*,*) 'ENTER OUTFILE' - read(*,*) outfile - open(newunit=file_unit_2, file = outfile) - open(newunit=file_unit_3, file = 'TfPsatIn.txt') - ! read input file:read number fluid, single carbon number for each fluid. - read(file_unit_1,*) Nfluid - do j = 1, Nfluid - read (file_unit_1,*) Oil - read (file_unit_1,*) Ncut, Ndef - end do - - nNcut = Ncut - nNdef = Ndef - - - ! read defined compounds with molar fraction and M by fluid j and compound i. - do i = 1, Ndef - read(file_unit_1,*) DefComp(i), zdef(i), rMdef(i) - end do - ! read carbon number and molar fraction of fluid j and comp. i - do i = 1, Ncut - read(file_unit_1,*) rCN(i), zcomp(i), rMW(i),Den(i) - end do - - read(file_unit_1,*) Plus, zPlus, rMWplus, DenPlus ! plus fraction properties - read (file_unit_1,*) Nps !number of pseudos for the C20+ fraction - read (file_unit_1,*) wat(1:Ndef+Ncut),watPlus ! atmospheric oil weight fractions - - !zpl = zPlus ! added by oscar - ! There get Zi*Mi and volume from cut to plus fraction - - zMdef(1:Ndef) = zdef(1:Ndef)*rMdef(1:Ndef) - zM(1:Ncut) = zcomp(1:Ncut) * rMW(1:Ncut) - zMp = zPlus * rMWplus - zMtotal = sum(zMdef(1:Ndef)) + sum(zM(1:Ncut)) + zMp - w(1:Ndef) = zMdef(1:Ndef) / zMtotal - w(Ndef+1:Ndef+Ncut) = zM(1:Ncut) / zMtotal - w(Ndef+Ncut+1) = zMp / zMtotal - rn(1:Ndef) = zdef(1:Ndef)/zMtotal - Z6p = sum(zcomp(1:Ncut)) + zPlus - sumV = sum(zM(1:Ncut)/Den(1:Ncut)) + zMp/DenPlus - - !agregado para calcular el peso molecular global - !MWt(1:Ndef) = rMdef(1:Ndef) - !MWt(Ndef+1:Ndef+Ncut) = rMW(1:Ncut) - - !do i= 1, Ndef+Ncut - ! print*, MWt(i) - !end do - - - - - !print*, Zplus - end subroutine data_in - -end module data - -module routines_pre - - use constants, only: pr - !use dtypes, only: FluidDataOut - use data - - - - contains - - subroutine Linear_Regression(x,y,a,b,r2) - !! This subroutine computes the regression line for a data set of x, y variables. - implicit none - - integer :: i - integer :: n !! total number of data set. - real(pr), allocatable, intent(in) :: x(:) !! x: input array of length n which contains the set of independent variable - real(pr), allocatable, intent(in) :: y(:) !! y: input array of length n which contains the set of dependent variable - real(pr) :: t1,t2,t3,t4, aux1,aux2,aux3,aux4,aux5 ! internal variables - real(pr), intent(out) :: a !! a: output real variable. Slope of the regression line - real(pr), intent(out) :: b !! b: output real variable. Intercept of the regression line - real(pr), intent(out) :: r2 !! r2: output real variable. Square correlation coefficient - - n = size(x) - !Calculation of a y b --> y=a*x+b - t1=0.;t2=0.;t3=0.;t4=0.; - - do i=1,n - t1=t1+x(i)*y(i) - t2=t2+x(i) - t3=t3+y(i) - t4=t4+x(i)**2 - end do - - a=(n*t1-t2*t3)/(n*t4-t2**2) - b=(t3-a*t2)/n - !coefficient calculation of correlations r2 - aux1=0.;aux2=0.;aux3=0.;aux4=0.;aux5=0.; - - do i=1, n - aux1= aux1 + x(i)*y(i) - aux2= aux2 + x(i) - aux3= aux3 + y(i) - aux4= aux4 + x(i)**2 - aux5= aux5 + y(i)**2 - end do - - r2=(aux1-aux2*aux3/n)**2 /((aux4-aux2**2/n)*(aux5-aux3**2/n)) - - end subroutine Linear_Regression - - - subroutine Best_Linear_Regression(scn_nc,scn,scn_z,plus_z,a,b,r2,ninit) - !! This subroutine calculates the best regression line for an oil. - implicit none - - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the oil - integer, allocatable, intent(in) :: scn(:) !! set of singles cuts being considered in the oil - integer :: i, j, k , kold, nbest, xaux,cbmax ! internal variables - real, allocatable, intent(in) :: scn_z(:) !! set of corresponding mole fractions of scn cuts - real(pr), intent(in) :: plus_z - real(pr), intent(out) :: a !! output real variable. Slope of the best regression line. - real(pr), intent(out) :: b !! output real variable. Intercept of the best regression line. - real(pr), intent(out) :: r2 !! output real variable. Square correlation coefficient. - integer, intent(out) :: ninit !! minimum carbon number obtained from the best linear regression - real(pr), allocatable :: xBR(:), yBR(:) - real(pr) :: r2old,r2best, aold, bold, abest, bbest, zsum,zaux - - allocate(scn(scn_nc)) - allocate(scn_z(scn_nc)) - allocate(xBR(scn_nc)) - allocate(yBR(scn_nc)) - - k=5 - r2=0.0001d0 - r2old=0.00001d0 - r2best=0.0001d0 - - do while (r2.gt.r2old.or.r2old.lt.0.9) - kold=k - r2old=r2 - aold=a - bold=b - - if (r2.gt.r2best) then - r2best=r2 - abest=a - bbest=b - nbest=scn(scn_nc-k+2) - end if - - if (k.gt.scn_nc) then - if (r2.gt.r2best)then - ninit=scn(scn_nc-k+2) - go to 22 - else - r2=r2best - a=abest - b=bbest - ninit=nbest - go to 22 - end if - end if - - j=1 - xBR=0. - yBR=0. - - do i=scn_nc-k+1, scn_nc - xBR(j)=scn(i) - yBR(j)=scn_z(i) - j=j+1 - end do - - call LinearRegression(xBR(:k), yBR(:k), a, b, r2) - k=k+1 - end do - - r2=r2old - a=aold - b=bold - ninit=scn(scn_nc-kold+2) - - 22 continue - - zsum = 0d0 - xaux=scn(scn_nc) - - do while (zsum.lt.plus_z.and.xaux<300) - xaux=xaux+1 - zaux= exp(a*xaux+b) - zsum=zsum+zaux - end do - - cbmax=xaux - - print*, a, b, cbmax, r2, Ninit - - end subroutine Best_Linear_Regression - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine get_C(start,C,Ncut,rCN,zM,zMp,Z6p,a,b,rMp,maxC,z,zpi,zMpi,Ndef,zp,w,rn,rMWplus,zcomp,zdef) - - implicit none - - integer, parameter :: imax=48, maxD=15 - logical :: start - real(pr) :: C, Z6p, zMp, a, b, rMp, rMWplus - real(pr) :: aux, dold, dif, Cold,zp - real(pr), dimension(imax) :: zM ! inputs coming from main program - integer, dimension(imax) :: rCN ! inputs coming from main program - real(pr), dimension(imax) :: z, rM,w,rn, zcomp - integer :: Ncut, maxC, Ndef - real(pr), dimension(300):: zpi, zMpi - real(pr), dimension(maxD) :: zdef - - !write(*,*) 'ENTER THE VERSION NUMBER' - !write(*,*) '1: CMOD 2: CMM' - !read(*,*) nv - - - if(nv1==1)then - - ! Find the C value for molecular weights line - C = 13.0 ! initial guess - Start = .true. - rMp = rMWplus - - call difMpfromC(start,C,Ncut,rCN,zM,zMp,Z6p,dold,a,b,rMp,maxC,z,zpi,zMpi,zcomp,zdef,zp) - Cold = C - C = min(12.8,C-0.07) - Start = .false. - - - - call difMpfromC(start,C,Ncut,rCN,zM,zMp,Z6p,dif,a,b,rMp,maxC,z,zpi,zMpi,zcomp,zdef,zp) - - do while (abs(dif) > 0.1) - aux = C - C = C - dif*(C-Cold)/(dif-dold) - Cold = aux - dold = dif - call difMpfromC(start,C,Ncut,rCN,zM,zMp,Z6p,dif,a,b,rMp,maxC,z,zpi,zMpi,zcomp,zdef,zp) - end do - - !print*, C,maxC, rMp - - - else - - ! Find the C value for molecular weights line - C = 14.0 ! initial guess and maximum C allowed - Start = .true. - rMp = rMWplus - - call difMpfromCfull (Start,C,Ndef,Ncut,rCN,dold,a,b,rM,rMp,maxC,z,zp,zpi,zMpi,w,rn,zcomp,zdef) - Cold = C - C = 13.5 !min(12.8,C-0.07) - Start = .false. - - call difMpfromCfull (Start,C,Ndef,Ncut,rCN,dif,a,b,rM,rMp,maxC,z,zp,zpi,zMpi,w,rn, zcomp,zdef) - - do while (abs(dif) > 0.1) - aux = C - C = C - dif*(C-Cold)/(dif-dold) - Cold = aux - dold = dif - call difMpfromCfull (Start,C,Ndef,Ncut,rCN,dif,a,b,rM,rMp,maxC,z,zp,zpi,zMpi,w,rn,zcomp,zdef) - end do - - !print*, C, maxC, rMp, zp - - end if - - end subroutine get_C - - subroutine difMpfromC(start,C,Ncut,rCN,zM,zMp,Z6p,dif,a,b,rMp,maxC,z,zpi,zMpi,zcomp,zdef,zp) - - use data - - implicit none - - integer, parameter :: imax=48, maxD=15 - logical :: start - real(pr) :: C, Z6p, zMp, dif, a, b, rMp, zp,aBE,bBE,r2, alim, blim, half - real(pr) :: a60,b60, sumz, rMpcalc - real(pr), dimension(imax) :: zM ! inputs coming from main program - integer, dimension(imax) :: rCN,x ! inputs coming from main program - real(pr), dimension(imax) :: ylog, rM, z, zcomp - integer :: i, Ncut, maxC, i0, CmaxL, C60max - real(pr), dimension(300):: zpi, zMpi - real(pr), dimension(maxD) :: zdef - - ! 1 i0 = rCN(Ncut) - - i = 300 - sumz = -1 - do while (i==300.and.sumz12 y C<12. - a = alim - b = blim - else - !call LineC60max (Ncut,Zp,rCN,aBE,bBE,half,a60,b60,C60max) - if(aBE > a60)then - a = a60 - b = b60 - else - a = aBE - b = bBE - end if - end if - - sumz = 0.0d0 - i = 0 - do while (sumz a60)then - a = a60 - b = b60 - else - a = aBE - b = bBE - end if - end if - - - sumz = 0.0d0 - i = 0 - do while (sumzalim) eliminado por oscar 20/03/2024 - ! aold = a - ! a = max(1.0001*a, alim) - ! b = b - (a-aold)*19.51 - ! sumz = 0.0d0 - ! i = 0 - ! do while (sumz 0.00001) - aux = rMp - rMp = rMp - dif*(rMp-rMpold)/(dif-dold) - rMpold = aux - dold = dif - call difnewMp(Start,C,Ndef,Ncut,rCN,dif,a,b,rM,rMp,maxC,z,zp,zpi,zMpi,w,rn,zcomp,zdef,rMW,rMdef, Mglobal) - end do - - - - - !print*, C, maxC, rMp - !print*, zp , z6p - !print*, z, zp - end if - - - end subroutine GetNewMpfromC - - - subroutine difnewMp(Start,C,Ndef,Ncut,rCN,dif,a,b,rM,rMp,maxC,z,zp,zpi,zMpi,w,rn,zcomp,zdef,rMW,rMdef, Mglobal) - !! this subroutines for a C value returns the correponding Mp value - !use MassFracAndMoles, only: w,rn - use data - implicit none - - integer, parameter :: imax=48, maxD=15 - logical :: start - real(pr) :: C, dif, a, b, rMp, zp ,aBE ,bBE ,r2 , alim, blim, half - real(pr) :: a60,b60, sumz, rMpcalc, rntot, aold, Mglobal - integer, dimension(imax) :: rCN,x ! inputs coming from main program - real(pr), dimension(imax) :: ylog, rM, z, w, rn, zcomp, ztotal,rMW, MWt - integer :: i,Ndef, Ncut, maxC, i0, CmaxL, C60max - real(pr), dimension(300):: zpi, zMpi,ZM_1 - real(pr), dimension(maxD) :: zdef,rMdef - - 1 i0 = rCN(Ncut) - - if(nv == 1)then !introduce 'if' to choise where experimental data from - x = rCN - rn(Ndef+1:Ndef+Ncut) = w(Ndef+1:Ndef+Ncut)/rMW(1:Ncut) - rn(Ndef+Ncut+1) = w(Ndef+Ncut+1)/rMp - rntot = sum(rn(1:Ndef+Ncut+1)) - ztotal(1:Ndef+Ncut) = rn(1:Ndef+Ncut) / rntot - z(1:Ncut) = rn(Ndef+1:Ndef+Ncut) / rntot - Zp = rn(Ndef+Ncut+1) / rntot - !ylog = log(zcomp) - ylog = log(z) - - !calculo m global - MWt(1:Ndef) = rMdef(1:Ndef) - MWt(Ndef+1:Ndef+Ncut) = rMW(1:Ncut) - ZM_1(1:Ndef+Ncut) = ztotal(1:Ndef+Ncut)*MWt(1:Ndef+Ncut) - - !do i= 1, Ndef+Ncut - ! print*, MWt(i) - !end do - - end if - - if(nv == 2)then - rM(1:Ncut) = 84 + C*(rCN(1:Ncut)-6) - rn(Ndef+1:Ndef+Ncut) = w(Ndef+1:Ndef+Ncut)/rM(1:Ncut) - rn(Ndef+Ncut+1) = w(Ndef+Ncut+1)/rMp - rntot = sum(rn(1:Ndef+Ncut+1)) - ztotal(1:Ndef+Ncut) = rn(1:Ndef+Ncut) / rntot - z(1:Ncut) = rn(Ndef+1:Ndef+Ncut) / rntot - Zp = rn(Ndef+Ncut+1) / rntot - x = rCN - ylog = log(z) - !calculo m global - MWt(1:Ndef) = rMdef(1:Ndef) - MWt(Ndef+1:Ndef+Ncut) = rM(1:Ncut) - ZM_1(1:Ndef+Ncut) = ztotal(1:Ndef+Ncut)*MWt(1:Ndef+Ncut) - - !do i= 1, Ndef+Ncut - ! print*, MWt(i) - !end do - - end if - !print*, zp - !print*, ztotal - - call BestLinearRegression(Ncut,x,ylog,Zp,aBE,bBE,r2) - call LimitLine(Ncut,rCN,Zp,aBE,bBE,half,alim,blim,CmaxL) - call LineC60max (Ncut,Zp,rCN,aBE,bBE,half,a60,b60,C60max) - - if(aBE < alim)then !, para C se excede el valor de 14 en la mayoria de los casos. - a = alim - b = blim - else - !call LineC60max (Ncut,Zp,rCN,aBE,bBE,half,a60,b60,C60max) - if(aBE > a60)then - a = a60 - b = b60 - else - a = aBE - b = bBE - end if - end if - - - sumz = 0.0d0 - i = 0 - do while (sumzalim) - aold = a - a = max(1.0001*a, alim) - b = b - (a-aold)*19.51 - sumz = 0.0d0 - i = 0 - do while (sumz y=a*x+b - t1=0.;t2=0.;t3=0.;t4=0.; - - do i=1,n - t1=t1+x(i)*y(i) - t2=t2+x(i) - t3=t3+y(i) - t4=t4+x(i)**2 - end do - - a=(n*t1-t2*t3)/(n*t4-t2**2) - b=(t3-a*t2)/n - !coefficient calculation of correlations r2 - aux1=0.;aux2=0.;aux3=0.;aux4=0.;aux5=0.; - - do i=1, n - aux1= aux1 + x(i)*y(i) - aux2= aux2 + x(i) - aux3= aux3 + y(i) - aux4= aux4 + x(i)**2 - aux5= aux5 + y(i)**2 - end do - - r2=(aux1-aux2*aux3/n)**2 /((aux4-aux2**2/n)*(aux5-aux3**2/n)) - - end subroutine Linear_Regression - - subroutine Best_Linear_Regression(scn_nc,scn,scn_z,plus_z,a,b,r2,n_init,c_max_blr) - !! This subroutine calculates the best regression line for an oil. - implicit none - - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the oil - integer, allocatable, intent(in) :: scn(:) !! set of singles cuts being considered in the oil - integer, intent(out) :: c_max_blr - !! output CN at which plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) - real(pr), allocatable, intent(in) :: scn_z(:) !! set of corresponding mole fractions of scn cuts - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file - real(pr), intent(out) :: a !! output real variable. Slope of the best regression line. - real(pr), intent(out) :: b !! output real variable. Intercept of the best regression line. - real(pr), intent(out) :: r2 !! output real variable. Square correlation coefficient. - integer, intent(out) :: n_init !! minimum carbon number obtained from the best linear regression - integer :: i, j, k , k_old, n_best, x_aux ! internal variables - real(pr), dimension(scn_nc) :: x_blr, y_blr - real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux - - k=5 - r2=0.0001_pr - r2_old=0.00001_pr - r2_best=0.0001_pr - - do while (r2.gt.r2_old.or.r2_old.lt.0.9) - k_old=k - r2_old=r2 - a_old=a - b_old=b - - if (r2.gt.r2_best) then - r2_best=r2 - a_best=a - b_best=b - n_best=scn(scn_nc-k+2) - end if - - if (k.gt.scn_nc) then - if (r2.gt.r2_best)then - n_init=scn(scn_nc-k+2) - go to 22 - else - r2=r2_best - a=a_best - b=b_best - n_init=n_best - go to 22 - end if - end if - - j=1 - x_blr=0. - y_blr=0. - - do i=scn_nc-k+1, scn_nc - x_blr(j)=scn(i) - y_blr(j)=scn_z(i) - j=j+1 - end do - - call Linear_Regression(x_blr(:k), y_blr(:k), a, b, r2) - k=k+1 - end do - - r2=r2_old - a=a_old - b=b_old - n_init=scn(scn_nc-k_old+2) - - 22 continue - - z_sum = 0d0 - x_aux=scn(scn_nc) - - do while (z_sum.lt.plus_z.and.x_aux<300) - x_aux=x_aux+1 - z_aux= exp(a*x_aux+b) - z_sum=z_sum+z_aux - end do - - c_max_blr=x_aux - - !print*, a, b, c_max_blr, r2, n_init - - end subroutine Best_Linear_Regression - - subroutine LimitLine(scn_nc,scn,plus_z,a_blr,b_blr,a_lim,b_lim,c_max_lim,half) - !! This subroutine obtains the limit line constants for an oil. - implicit none - - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the oil - integer, allocatable, intent(in) :: scn(:) !! set of singles cuts being considered in the oil - integer, intent(out) :: c_max_lim - !! output CN at which plus_z is reached, as the summation of single z(i) from the limit distribution (blr) - real(pr), intent(in) :: a_blr !! input constant from the Best linear regression - real(pr), intent(in) :: b_blr !! input constant from the Best linear regression - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file - real(pr), intent(out) :: a_lim !! output real variable. Slope of the limit line - real(pr), intent(out) :: b_lim !! output real variable. Intercept of the limit line. - real(pr), intent(out) :: half - real(pr) :: z_lim, cross_cn, z_cross, z_aux ! cn: carbon number - integer :: x_aux - - !A and B limit calculation. - z_lim = 0.0_pr - half = 0.506_pr ! modified by oscar, the variable half was removed fron the argument - !print*, scn - do while (z_lim.lt.plus_z) - half = half + 0.002d0 - cross_cn = scn(scn_nc)+half ! typically 19.508 in first try - z_cross = exp((a_blr*cross_cn) + b_blr) - a_lim = -(z_cross)/plus_z - b_lim = log(z_cross)-(a_lim*cross_cn) - x_aux=scn(scn_nc) - z_lim = 0._pr - do while (z_lim.lt.plus_z.and.x_aux<300) - x_aux=x_aux+1 - z_aux= exp(a_lim*x_aux+b_lim) - z_lim=z_lim+z_aux - end do - continue - end do - - c_max_lim = x_aux - !print*, a_lim,b_lim,c_max_lim - - end subroutine LimitLine - - subroutine Line_C60_max (scn_nc,scn,plus_z,a_blr,b_blr,half,a_60,b_60,c_max_60) - !! This subroutine obtains the Cmax60 line constants for an oil - - implicit none - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the oil - integer, allocatable, intent(in) :: scn(:) !! set of singles cuts being considered in the oil - integer, intent(out) :: c_max_60 !!output CN at which Zp is reached, as the summation of single z(i) from the Cmax60 distribution. - real(pr), intent(in) :: a_blr !! input constant from the Best linear regression - real(pr), intent(in) :: b_blr !! input constant from the Best linear regression - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file - real(pr), intent(out) :: a_60 !! output real variable. Slope of the limit line - real(pr), intent(out) :: b_60 !! output real variable. Intercept of the C60max line. - real(pr), intent(in) :: half - integer :: x_aux - real(pr) :: z_sum, cross_cn, z_cross, z_aux, F_tol, a_tol, var_range - real(pr) :: a_old, Zp_60, F, dF_dA - - cross_cn = scn(scn_nc)+half ! typically 19.51 - z_cross= exp((a_blr*cross_cn)+b_blr) - F_tol=1_pr - a_tol=1_pr - a_60=a_blr - var_range = 60.5_pr - cross_cn - - do while (a_tol.gt.1e-6_pr.or.F_tol.gt.1e-6_pr) - a_old=a_60 - Zp_60 = (exp(a_60*var_range+log(z_cross))-z_cross)/a_60 - F = plus_z - Zp_60 - dF_dA = - (var_range*(a_60*Zp_60+z_cross)-Zp_60) / a_60 ! modified by oscar, error in the derivative in previous version is corrected - a_60=a_old-(F/dF_dA) - a_tol=abs(a_old-a_60) - F_tol=abs(F) - end do - - b_60=log(z_cross)-a_60*cross_cn - z_sum=0.0_pr - x_aux=scn(scn_nc) - - do while (z_sum.lt.plus_z) - x_aux=x_aux+1 - z_aux= exp(a_60*x_aux+b_60) - z_sum=z_sum+z_aux - end do - - c_max_60=x_aux - - !print*, a_60,b_60,c_max_60 - end subroutine Line_C60_max - - subroutine select_method(oil,mw_source,method,C,log_scn_z,plus_z,plus_mw) - - implicit none - type(FluidData) :: oil - !character(len=*), intent(in) :: file !! file name - character(len=*), intent(in), optional :: method - character(len=*), intent(in) :: mw_source - real(pr), allocatable :: scn_mw(:) !! set to molecular weights calculated by \[M = 84-C(i-6)\] - real(pr) :: C !! C constants which is used in equation \[M = 84-C(i-6)\] - real(pr), allocatable :: scn_z(:) !! set of corresponding calculated mole fractions of scn cuts - real(pr) :: sum_def_comp_z_plus !! compositions sum from (define component +1) to (plus component) - real(pr), intent(out) :: plus_z !! composition of residual fraction from input file - real(pr) :: plus_mw !! calculated molecular weight of residual fraction - real(pr), allocatable, intent(out) :: log_scn_z(:) !! logarithm of corresponding mole fractions of scn cuts - real(pr), allocatable :: def_comp_moles(:) - real(pr), allocatable :: scn_moles(:) - real(pr) :: plus_moles - real(pr) :: total_moles - - - !oil = data_from_file(file) - - allocate(scn_z(oil%scn_nc)) - allocate(scn_mw(oil%scn_nc)) - allocate(log_scn_z(oil%scn_nc)) - allocate(def_comp_moles(oil%def_comp_nc)) - allocate(scn_moles(oil%scn_nc)) - - def_comp_moles = (oil%def_comp_w)/(oil%def_comp_mw) - - select case (mw_source) - - case("experimental") - scn_moles = (oil%scn_w)/(oil%scn_mw) - plus_moles = (oil%plus_w)/(plus_mw) - total_moles = sum(def_comp_moles)+sum(scn_moles)+(plus_moles) - scn_z = scn_moles / total_moles - plus_z = plus_moles / total_moles - log_scn_z = log(scn_z) - - case("calculated") - - if (method=="global_mw")then - scn_mw = 84 + C*(oil%scn-6) - scn_z = oil%product_z_mw_scn / scn_mw - sum_def_comp_z_plus = sum(oil%scn_z) + (oil%plus_z) - plus_z = sum_def_comp_z_plus - sum(scn_z) ! new molar fraction for residual cut, based on C - plus_mw = oil%product_z_mw_plus / plus_z - log_scn_z = log(scn_z) - end if - - if (method=="plus_mw")then - scn_mw = 84 + C*(oil%scn-6) - scn_moles = (oil%scn_w)/(scn_mw) - plus_moles = (oil%plus_w)/(plus_mw) - total_moles = sum(def_comp_moles)+sum(scn_moles)+(plus_moles) - scn_z = scn_moles / total_moles - plus_z = plus_moles / total_moles - log_scn_z = log(scn_z) - endif - - end select - - end subroutine select_method - - subroutine difference_mw_plus(oil,mw_source,method,start,C,difference,plus_mw) - - !! this subroutine ... - implicit none - type(FluidData) :: oil - logical :: start - !character(len=*), intent(in) :: file !! file name - character(len=*), intent(in), optional :: method - character(len=*), intent(in) :: mw_source - integer :: c_max !! output CN at which plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) - real(pr), allocatable :: scn_z(:) !! set of corresponding calculated mole fractions of scn cuts - integer :: n_init ! minimum carbon number obtained from the best linear regression - real(pr), allocatable :: log_scn_z(:) !! logarithm of corresponding mole fractions of scn cuts - real(pr) :: C !! C constants which is used in equation \[M = 84-C(i-6)\] - real(pr) :: plus_z !! composition of residual fraction from input file - real(pr) :: plus_mw !! calculated molecular weight of residual fraction - real(pr) :: plus_mw_cal !! calculated molecular weight of residual fraction - real(pr), intent(out) :: difference !! - real(pr) :: a !! output real variable. Slope of the best regression line. - real(pr) :: b !! output real variable. Intercept of the best regression line. - real(pr) :: r2 !! output real variable. Square correlation coefficient. - real(pr), allocatable :: x_blr(:), y_blr(:) - real(pr) :: a_blr, b_blr, a_60, b_60, a_lim, b_lim, half - real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux - real(pr), dimension(300) :: plus_z_i - real(pr), dimension(300) :: product_z_mw_plus_i - !real(pr) :: plus_mw - real(pr):: sum_z - integer :: j, k , k_old, n_best, x_aux, i_0! internal variables - integer :: c_max_blr, c_max_lim, c_max_60 - integer :: i - real(pr), allocatable :: scn_mw(:) !! set to molecular weights calculated by \[M = 84-C(i-6)\] - real(pr) :: sum_def_comp_z_plus !! compositions sum from (define component +1) to (plus component) - integer, dimension(300) :: carbon_number_plus - real(pr) :: denom - - !oil = data_from_file(file) - - allocate(scn_z(oil%scn_nc)) - allocate(scn_mw(oil%scn_nc)) - allocate(log_scn_z(oil%scn_nc)) - allocate(x_blr(oil%scn_nc)) - allocate(y_blr(oil%scn_nc)) - - 1 i_0 = oil%scn(oil%scn_nc) - - call select_method(oil,mw_source,method,C,log_scn_z,plus_z,plus_mw) ! select case - call Best_Linear_Regression(oil%scn_nc,oil%scn,log_scn_z,plus_z,a_blr,b_blr,r2, & - n_init,c_max_blr) - call LimitLine(oil%scn_nc,oil%scn,plus_z,a_blr,b_blr,a_lim,b_lim,c_max_lim,half) - call Line_C60_max (oil%scn_nc,oil%scn,plus_z,a_blr,b_blr,half,a_60,b_60,c_max_60) ! add line by oscar. - - if(a_blr < a_lim)then ! added by oscar 05/12/2023.se elimina por que se agregan las restriccones para C>12 y C<12. - a = a_lim - b = b_lim - else ! this line was removed call Line_C60_max - if(a_blr > a_60)then - a = a_60 - b = b_60 - else - a = a_blr - b = b_blr - end if - end if - - sum_z = 0.0_pr - i = 0 - do while (sum_z12 .and. C<14 )then !! mayor igual a 12 o menor igual 14. - denom = sum((plus_z_i(1:i))*((carbon_number_plus(1:i)-6))) !! denominator of equation to obtain C value directly - C = (plus_z*(plus_mw-84))/denom - endif - - plus_mw_cal = sum(product_z_mw_plus_i(1:i))/plus_z - difference = plus_mw_cal-plus_mw - - end subroutine difference_mw_plus - - subroutine get_C_or_m_plus(oil,mw_source,method,start,C) - !! This subroutine... - implicit none - type(FluidData) :: oil - logical :: start - !character(len=*), intent(in) :: file !! file name - character(len=*), intent(in), optional :: method - character(len=*), intent(in) :: mw_source - integer :: c_max !! output CN at which plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) - real(pr), allocatable :: log_scn_z(:) !! logarithm of corresponding mole fractions of scn cuts - real(pr) :: C !! C constants which is used in equation \[M = 84-C(i-6)\] - real(pr) :: plus_z !! composition of residual fraction from input file - real(pr) :: plus_mw !! calculated molecular weight of residual fraction - real(pr) :: difference !! - real(pr) :: difference_old, plus_mw_old - real(pr) :: C_old - real(pr) :: aux - - !oil = data_from_file(file) - - select case (mw_source) - - case("experimental") - C = 13.5 - Start = .true. - plus_mw = oil%plus_mw - call difference_mw_plus(oil,mw_source,method,start,C,difference,plus_mw) - print*, C, plus_mw - - case("calculated") - - if(method == 'global_mw') C=13 - if(method == 'plus_mw') C=14 - - Start = .true. - plus_mw = oil%plus_mw - - call difference_mw_plus(oil,mw_source,method,start,C,difference_old,plus_mw) - C_old = C - if(method == 'global_mw') C=min(12.8, C-0.07) - if(method == 'plus_mw') C=13.5 - Start = .false. - call difference_mw_plus(oil,mw_source,method,start,C,difference,plus_mw) - - do while (abs(difference) > 0.1) - aux = C - C = C - difference*(C-C_old)/(difference-difference_old) - C_old = aux - difference_old = difference - call difference_mw_plus(oil,mw_source,method,start,C,difference,plus_mw) - end do - print*, C, plus_mw - end select - - if(C>14.or.C<12)then - - if(C>14) C=14 - if(C<12) C=12 - - Start = .true. - plus_mw = oil%plus_mw - call difference_mw_plus(oil,mw_source,method,start,C,difference_old,plus_mw) - plus_mw_old = plus_mw - plus_mw = 0.9*plus_mw - Start = .false. - call difference_mw_plus(oil,mw_source,method,start,C,difference,plus_mw) - - do while (abs(difference)>0.00001) - aux = plus_mw - plus_mw = plus_mw - difference*(plus_mw-plus_mw_old)/(difference-difference_old) - plus_mw_old = aux - difference_old = difference - call difference_mw_plus(oil,mw_source,method,start,C,difference,plus_mw) - end do - print*, C, plus_mw - end if - end subroutine get_C_or_m_plus - - end module pruebas - - - \ No newline at end of file diff --git a/src/routines.f90 b/src/routines.f90 index 24e8eb9..786b31f 100644 --- a/src/routines.f90 +++ b/src/routines.f90 @@ -589,7 +589,7 @@ subroutine lump (fluid, characterization) integer :: i, prev_i integer, dimension(15) :: j_ps integer :: i_ps ! auxiliary variables - real(pr) :: rec_zm, remain_plus_zm, sum_z, sum_zm, sum_volume_ps ! auxiliary variables + real(pr) :: rec_zm, remain_plus_zm, sum_z, sum_zm, sum_volume_ps ! auxiliary variables real(pr), allocatable :: plus_z_ps(:) real(pr), allocatable :: plus_mw_ps(:) integer :: numbers_ps @@ -597,6 +597,9 @@ subroutine lump (fluid, characterization) real(pr), allocatable :: w_ps(:) real(pr) :: sum_zm_last_ps real(pr), allocatable :: moles(:) + real(pr), allocatable :: lumped_z(:) + real(pr), allocatable :: lumped_mw(:) + real(pr), allocatable :: lumped_densities(:) associate (& def_nc => fluid%def_comp_nc , & @@ -610,6 +613,7 @@ subroutine lump (fluid, characterization) w => fluid%w, & plus_mw => characterization%plus_mw, & carbon_number_plus => characterization%carbon_number_plus, & + plus6_density => characterization%plus6_density, & C => characterization%C & ) @@ -626,6 +630,10 @@ subroutine lump (fluid, characterization) allocate(density_ps(0)) allocate (w_ps(0)) allocate(moles(0)) + allocate(lumped_z(0)) + allocate(lumped_mw(0)) + allocate(lumped_densities(0)) + scn_nc_input = fluid%scn_nc scn_nc_new = fluid%scn_nc_ps - 6 @@ -654,6 +662,10 @@ subroutine lump (fluid, characterization) characterization%product_z_mw_plus_i = z_m_plus_i characterization%plus_z_i = plus_z_i end if + + !do i= 1, scn_nc_new + ! print*, scn_z(i), scn_mw(i) + !end do numbers_ps = fluid%numbers_ps j_ps = 0 @@ -674,8 +686,7 @@ subroutine lump (fluid, characterization) j_ps(i_ps) = j_ps(i_ps) + 1 sum_z = sum_z + plus_z_i(i) sum_zm = sum_zm + z_m_plus_i(i) - sum_volume_ps = sum_volume_ps + z_m_plus_i(i) / & - characterization%plus6_density(scn_nc_new+i) + sum_volume_ps = sum_volume_ps + z_m_plus_i(i) / plus6_density(scn_nc_new+i) if (z_m_plus_i(i+1) > 2*(rec_zm - sum_zm)) then ! when adding one more would go too far plus_z_ps = [plus_z_ps, sum_z] @@ -701,9 +712,8 @@ subroutine lump (fluid, characterization) carbon_number_plus(scn_nc_new + numbers_ps) = 6 + & ((plus_mw_ps(numbers_ps) - 84) / C) sum_zm_last_ps = (sum(z_m_plus_i(i+1:i_last))) - density_ps = [density_ps, ((sum_zm_last_ps)/(sum(z_m_plus_i(i+1:i_last)) & - /(characterization%plus6_density(scn_nc_new+i+1:i_last))))] - j_ps(numbers_ps) = i_last - sum(j_ps(1: numbers_ps-1)) + density_ps = [density_ps, (sum_zm_last_ps)/sum((z_m_plus_i(i+1:i_last)) / (plus6_density(scn_nc_new+i+1:last)))] + j_ps(numbers_ps) = i_last - sum(j_ps(1: numbers_ps-1)) !! revisar manana do i = 1, numbers_ps print*, 'ps',i, plus_z_ps(i), plus_mw_ps(i), density_ps(i), j_ps(i), carbon_number_plus(scn_nc_new+i) @@ -715,10 +725,23 @@ subroutine lump (fluid, characterization) moles = [moles, (fluid%def_comp_w)/(fluid%def_comp_mw)] moles = [moles, fluid%w(def_nc+1:def_nc+scn_nc_new)/ scn_mw(1:scn_nc_new)] moles = [moles, w_ps/plus_mw_ps ] - characterization%mol_fraction = moles/sum(moles) ! mole fractions normalize + characterization%mol_fraction = moles/sum(moles) ! mole fractions normalize characterization%last_C = last_C characterization%i_last = i_last + characterization%last = last + + ! In this point we create a new array to content final Mis and densities + ! whit pseudos. + + lumped_z = [scn_z(1:scn_nc_new), plus_z_ps] + lumped_mw = [scn_mw(1:scn_nc_new), plus_mw_ps] + lumped_densities = [plus6_density(1:scn_nc_new), density_ps ] + + + characterization%lumped_z = lumped_z + characterization%lumped_mw = lumped_mw + characterization%lumped_densities = lumped_densities end associate @@ -728,87 +751,52 @@ subroutine get_critical_constants(fluid, characterization, eos) !! this subroutine doing: !! - Calculation of Tc, Pc, omega from Pedersen correlations (adapted from code "C7plusPedersenTcPcOm") !! - Calculation of EoS parameters (taken from CubicParam) - + use critical_parameters + implicit none type(FluidData), intent(in) :: fluid type(FluidDataOut), intent(inout) :: characterization character(len=*), intent(in) :: eos + real(pr), allocatable :: modified_c(:) + real(pr), allocatable :: tc(:) + real(pr), allocatable :: pc(:) + real(pr), allocatable :: om(:) + real(pr), allocatable :: m_funtion(:) + integer :: i - select case (eos) - - case("srk") - ! Coeficientes originales Tabla 5.3 Pedersen - c1 = 1.6312e2_pr - c2 = 8.6052d1 - c3 = 4.3475d-1 - c4 = -1.8774d3 - d1P = -1.3408d-1 - d2 = 2.5019 - d3 = 2.0846d2 - d4 = -3.9872d3 - d5 = 1.0d0 - e1 = 7.4310d-1 - e2 = 4.8122d-3 - e3 = 9.6707d-3 - e4 = -3.7184d-6 - A = -0.176 - B = 1.574 - C0 = 0.48 - del1 = 1.0D0 - case("pr") - ! Coeficientes originales Tabla 5.3 Pedersen - c1 = 7.34043d1 - c2 = 9.73562d1 - c3 = 6.18744d-1 - c4 = -2.05932d3 - d1P = 7.28462d-2 - d2 = 2.18811d0 - d3 = 1.63910d2 - d4 = -4.04323d3 - d5 = 0.25d0 - e1 = 3.73765d-1 - e2 = 5.49269d-3 - e3 = 1.17934d-2 - e4 = -4.93049d-6 - A = -0.26992 - B = 1.54226 - C0 = 0.37464 - del1 = 1.0D0 + sqrt(2.0) - case("rkpr") - c1 = 7.34043d1 - c2 = 9.73562d1 - c3 = 6.18744d-1 - c4 = -2.05932d3 - d1P = 2.3d0 - d2 = -0.5d0 - d3 = 1.85d2 - d4 = -4.0d3 - d5 = 0.25d0 - e1 = 3.73765d-1 - e2 = 5.49269d-3 - e3 = 1.17934d-2 - e4 = -4.93049d-6 - A = -0.26992 - B = 1.54226 - C0 = 0.37464 - end select - - if (eos == "pr" .or. eos=="srk") - - d1 = (1 + del1**2)/(1 + del1) - y = 1 + (2*(1 + del1))**(1.0d0/3) + (4/(1 + del1))**(1.0d0/3) - OMa = (3*y*y + 3*y*d1 + d1**2 + d1 - 1.0d0)/(3*y + d1 - 1.0d0)**2 - OMb = 1/(3*y + d1 - 1.0d0) - Zc = y/(3*y + d1 - 1.0d0) - + associate(& + rho => characterization%lumped_densities, & + mw => characterization%lumped_mw & + ) + + allocate(modified_c(0)) + allocate(tc(0)) + allocate(pc(0)) + allocate(om(0)) + allocate(m_funtion(0)) + + call get_parameteres_for_critical(eos) + + tc = c1*rho + c2*log(mw) + c3*mw + (c4/mw) + pc = exp(d1_p + d2*(rho**(d5)) + (d3/mw) + (d4/(mw**2))) + + if (eos == "SRK") then + m_funtion = e1 + e2*mw + e3*rho + e4*(mw**(2)) + else + if (eos == "PR" .or. eos =="RKPR") then + m_funtion = e3*rho + 0.2 + 1.8*(1-exp(-mw/220)) + ! Alternative expression without max with MW - May 2019 + endif endif + modified_c = c0 - m_funtion + om = (-b + sqrt(b**2 - 4*a*modified_c))/(2*a) + do i = 1, size(mw) + print*, mw(i), rho(i), tc(i), pc(i), om(i), m_funtion(i) + end do - - - - + end associate end subroutine get_critical_constants @@ -834,11 +822,18 @@ type(FluidDataOut) function characterize(file, mw_source, method, fix_C, eos) & allocate(characterization%scn_i(0)) allocate(characterization%plus6_density(0)) allocate(characterization%mol_fraction(0)) - + allocate(characterization%lumped_z(0)) + allocate(characterization%lumped_mw(0)) + allocate(characterization%lumped_densities(0)) + allocate(characterization%critical_temperature(0)) + allocate(characterization%critical_pressure(0)) + allocate(characterization%acentric_factor(0)) + allocate(characterization%m_funtion(0)) call get_c_or_m_plus(fluid=fluid, mw_source=mw_source, method=method, fix_C=fix_C, characterization= characterization) call density_funtion(fluid=fluid, mw_source=mw_source, characterization=characterization) call lump(fluid=fluid, characterization=characterization) + call get_critical_constants(fluid= fluid, characterization=characterization, eos=eos) end function characterize diff --git a/src/types.f90 b/src/types.f90 index de5820e..5bd309a 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -60,11 +60,19 @@ module dtypes real(pr) :: volume_6plus_cal integer :: last_C integer :: i_last + integer :: last real(pr), allocatable :: scn_i(:) real(pr), allocatable :: plus6_density(:) real(pr) :: plus_w real(pr) :: plus_density real(pr), allocatable :: mol_fraction(:) + real(pr), allocatable :: lumped_z(:) + real(pr), allocatable :: lumped_mw(:) + real(pr), allocatable :: lumped_densities(:) + real(pr), allocatable :: critical_temperature(:) + real(pr), allocatable :: critical_pressure(:) + real(pr), allocatable :: acentric_factor(:) + real(pr), allocatable :: m_funtion (:) end type FluidDataOut From 93381b6d968c3e59c6b1c2924a703fe1681e0a6c Mon Sep 17 00:00:00 2001 From: Oscar Theran Date: Wed, 4 Sep 2024 16:12:20 -0300 Subject: [PATCH 05/16] changes 04/09/2024 --- app/main.f90 | 8 +--- fort.1 | 0 src/routines.f90 | 100 ++++++++++++++++++++++++++++++++++++++++++----- src/tools/io.f90 | 22 +++++++++++ src/types.f90 | 2 + 5 files changed, 115 insertions(+), 17 deletions(-) create mode 100644 fort.1 create mode 100644 src/tools/io.f90 diff --git a/app/main.f90 b/app/main.f90 index d85eaf3..49d4272 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -8,13 +8,7 @@ program main type(FluidData) :: fluid - prueba = characterize(file='oil1.nml', mw_source="calculated", method = "plus_mw", fix_C=.true., eos='PR') - print*, prueba%n_init - print*, prueba%a , prueba%b - print*, prueba%C, prueba%plus_mw - print*, prueba%a_d, prueba%b_d - print*, "-------------------------------------------------------------------" - + prueba = characterize(file='oil1.nml', mw_source="calculated", method = "plus_mw", fix_C=.false., eos='PR') end program main diff --git a/fort.1 b/fort.1 new file mode 100644 index 0000000..e69de29 diff --git a/src/routines.f90 b/src/routines.f90 index 786b31f..afd9560 100644 --- a/src/routines.f90 +++ b/src/routines.f90 @@ -564,14 +564,14 @@ subroutine calculate_volume_6plus(fluid, mw_source, characterization) plus6_density = [scn_density_cal, plus_density_cal] end select - characterization%volume_6plus_cal = volume_6plus_cal + characterization%volume_6plus_cal = volume_6plus_cal characterization%plus6_density = plus6_density end associate end subroutine calculate_volume_6plus - subroutine lump (fluid, characterization) + subroutine lump(fluid, characterization) !! This sobroutine ... implicit none type(FluidData), intent(in) :: fluid @@ -614,7 +614,8 @@ subroutine lump (fluid, characterization) plus_mw => characterization%plus_mw, & carbon_number_plus => characterization%carbon_number_plus, & plus6_density => characterization%plus6_density, & - C => characterization%C & + C => characterization%C, & + scn_nc_new => characterization%scn_nc_new & ) last_C = fluid%scn(fluid%scn_nc) !19 @@ -715,14 +716,14 @@ subroutine lump (fluid, characterization) density_ps = [density_ps, (sum_zm_last_ps)/sum((z_m_plus_i(i+1:i_last)) / (plus6_density(scn_nc_new+i+1:last)))] j_ps(numbers_ps) = i_last - sum(j_ps(1: numbers_ps-1)) !! revisar manana - do i = 1, numbers_ps - print*, 'ps',i, plus_z_ps(i), plus_mw_ps(i), density_ps(i), j_ps(i), carbon_number_plus(scn_nc_new+i) - end do + !do i = 1, numbers_ps + ! print*, 'ps',i, plus_z_ps(i), plus_mw_ps(i), density_ps(i), j_ps(i), carbon_number_plus(scn_nc_new+i) + !end do plus_density = plus_zm / sum(plus_z_ps(1:numbers_ps)* & plus_mw_ps(1:numbers_ps)/(density_ps(1:numbers_ps))) - moles = [moles, (fluid%def_comp_w)/(fluid%def_comp_mw)] + moles = [(fluid%def_comp_w)/(fluid%def_comp_mw)] moles = [moles, fluid%w(def_nc+1:def_nc+scn_nc_new)/ scn_mw(1:scn_nc_new)] moles = [moles, w_ps/plus_mw_ps ] @@ -792,14 +793,91 @@ subroutine get_critical_constants(fluid, characterization, eos) modified_c = c0 - m_funtion om = (-b + sqrt(b**2 - 4*a*modified_c))/(2*a) - do i = 1, size(mw) - print*, mw(i), rho(i), tc(i), pc(i), om(i), m_funtion(i) - end do + !do i = 1, size(mw) + ! print*, mw(i), rho(i), tc(i), pc(i), om(i), m_funtion(i) + !end do + + characterization%critical_temperature = tc + characterization%critical_pressure = pc + characterization%acentric_factor = om + characterization%m_funtion = m_funtion end associate end subroutine get_critical_constants + + + + subroutine write_result(fluid,characterization) + use ftools__io, only: str + use critical_parameters + + implicit none + type(FluidData), intent(in) :: fluid + type(FluidDataOut), intent(inout) :: characterization + integer :: i, i_prev + character(len=*), parameter :: str_fmt_1 = "(A4,7(A15,2x))" + character(len=*), parameter :: num_fmt_1 = "(A3,1x,*(E15.5,2x))" + character(len=*), parameter :: num_fmt_2 = "(A3,1x,*(E15.5,2x))" + + associate(& + def_nc => fluid%def_comp_nc, & + def_name => fluid%def_components, & + def_mw => fluid%def_comp_mw, & + z => characterization%mol_fraction, & + mw => characterization%lumped_mw, & + rho => characterization%lumped_densities, & + scn => fluid%scn, & + scn_nc_new => characterization%scn_nc_new, & + tc => characterization%critical_temperature, & + pc => characterization%critical_pressure, & + omega => characterization%acentric_factor, & + m_funtion => characterization%m_funtion & + ) + + print*, "----------------------------------------------------------------" + print*, "Best feasible regresion parameters" + print*, "" + print*, "Init_BFR:", characterization%n_init + print*, "A:", characterization%a," ","B:", characterization%b + print*, "C:", characterization%C + print*, "MW+:", characterization%plus_mw + print*, "Cmax:", characterization%c_max + print*, "" + print*, "Density funtion parameters parameters" + print*, "" + print*, "ad:", characterization%a_d," ","bd:", characterization%b_d + print*, "----------------------------------------------------------------" + print*, "" + print*, "compositional result" + print*, "" + + + write(*, fmt=str_fmt_1) "Comp", "Z", "Mw", "Tc", "Pc", "Omega", "Rho", "M" + ! defined components + do i = 1, def_nc + write(*, fmt=num_fmt_1) def_name(i), z(i), def_mw(i), tc_def(i), pc_def(i), om_def(i) + end do + ! + i_prev = scn(1) - 1 + do i = 1, scn_nc_new + write(*, fmt=num_fmt_2) 'C' // str(i_prev + i), z(i+def_nc), mw(i), tc(i), pc(i), omega(i), rho(i), m_funtion(i) + end do + ! + do i = 1 + scn_nc_new, size(mw) + write(*, fmt=num_fmt_2) 'ps' // str(i - scn_nc_new), z(i+def_nc), mw(i), rho(i), tc(i), pc(i), omega(i), m_funtion(i) + end do + + end associate + end subroutine write_result + + + + + + + type(FluidDataOut) function characterize(file, mw_source, method, fix_C, eos) & result(characterization) @@ -809,6 +887,7 @@ type(FluidDataOut) function characterize(file, mw_source, method, fix_C, eos) & character(len=*), intent(in) :: mw_source logical, intent(in), optional :: fix_C character(len=*), intent(in) :: eos + fluid = data_from_file(file) @@ -834,6 +913,7 @@ type(FluidDataOut) function characterize(file, mw_source, method, fix_C, eos) & call density_funtion(fluid=fluid, mw_source=mw_source, characterization=characterization) call lump(fluid=fluid, characterization=characterization) call get_critical_constants(fluid= fluid, characterization=characterization, eos=eos) + call write_result(fluid=fluid,characterization=characterization) end function characterize diff --git a/src/tools/io.f90 b/src/tools/io.f90 new file mode 100644 index 0000000..7fde863 --- /dev/null +++ b/src/tools/io.f90 @@ -0,0 +1,22 @@ +module ftools__io + implicit none + + interface str + module procedure :: str_gen + end interface str + + contains + function str_gen(int_in) result(str_out) + class(*), intent(in) :: int_in + character(len=99) :: str_mid + character(len=:), allocatable :: str_out + + select type(int_in) + type is (real) + write (str_mid, *) int_in + type is (integer) + write (str_mid, *) int_in + end select + str_out = trim(adjustl(str_mid)) + end function str_gen + end module ftools__io \ No newline at end of file diff --git a/src/types.f90 b/src/types.f90 index 5bd309a..8364ed9 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -61,6 +61,7 @@ module dtypes integer :: last_C integer :: i_last integer :: last + integer :: scn_nc_new real(pr), allocatable :: scn_i(:) real(pr), allocatable :: plus6_density(:) real(pr) :: plus_w @@ -73,6 +74,7 @@ module dtypes real(pr), allocatable :: critical_pressure(:) real(pr), allocatable :: acentric_factor(:) real(pr), allocatable :: m_funtion (:) + end type FluidDataOut From da9702a692d63ca5f4f64f2701be92277b47dd97 Mon Sep 17 00:00:00 2001 From: Oscar Theran Date: Thu, 5 Sep 2024 11:53:37 -0300 Subject: [PATCH 06/16] 05/09/2024 --- app/main.f90 | 9 +++- src/critical_parameters.f90 | 2 +- src/routines.f90 | 70 ++------------------------ src/types.f90 | 99 ++++++++++++++++++++++++++++++++----- 4 files changed, 102 insertions(+), 78 deletions(-) diff --git a/app/main.f90 b/app/main.f90 index 49d4272..75dc1aa 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -5,10 +5,17 @@ program main implicit none type(FluidDataOut) :: prueba - type(FluidData) :: fluid + !type(FluidData) :: fluid + integer :: i prueba = characterize(file='oil1.nml', mw_source="calculated", method = "plus_mw", fix_C=.false., eos='PR') + write(*, *) prueba + + !do i = 1, size(prueba%mol_fraction) + ! print*, prueba%mol_fraction(i) + !end do + end program main diff --git a/src/critical_parameters.f90 b/src/critical_parameters.f90 index e8b2377..a6baee5 100644 --- a/src/critical_parameters.f90 +++ b/src/critical_parameters.f90 @@ -1,6 +1,6 @@ module critical_parameters use constants - use dtypes, only: FluidData + implicit none diff --git a/src/routines.f90 b/src/routines.f90 index afd9560..cc25321 100644 --- a/src/routines.f90 +++ b/src/routines.f90 @@ -809,69 +809,7 @@ end subroutine get_critical_constants - subroutine write_result(fluid,characterization) - use ftools__io, only: str - use critical_parameters - - implicit none - type(FluidData), intent(in) :: fluid - type(FluidDataOut), intent(inout) :: characterization - integer :: i, i_prev - character(len=*), parameter :: str_fmt_1 = "(A4,7(A15,2x))" - character(len=*), parameter :: num_fmt_1 = "(A3,1x,*(E15.5,2x))" - character(len=*), parameter :: num_fmt_2 = "(A3,1x,*(E15.5,2x))" - - associate(& - def_nc => fluid%def_comp_nc, & - def_name => fluid%def_components, & - def_mw => fluid%def_comp_mw, & - z => characterization%mol_fraction, & - mw => characterization%lumped_mw, & - rho => characterization%lumped_densities, & - scn => fluid%scn, & - scn_nc_new => characterization%scn_nc_new, & - tc => characterization%critical_temperature, & - pc => characterization%critical_pressure, & - omega => characterization%acentric_factor, & - m_funtion => characterization%m_funtion & - ) - - print*, "----------------------------------------------------------------" - print*, "Best feasible regresion parameters" - print*, "" - print*, "Init_BFR:", characterization%n_init - print*, "A:", characterization%a," ","B:", characterization%b - print*, "C:", characterization%C - print*, "MW+:", characterization%plus_mw - print*, "Cmax:", characterization%c_max - print*, "" - print*, "Density funtion parameters parameters" - print*, "" - print*, "ad:", characterization%a_d," ","bd:", characterization%b_d - print*, "----------------------------------------------------------------" - print*, "" - print*, "compositional result" - print*, "" - - - write(*, fmt=str_fmt_1) "Comp", "Z", "Mw", "Tc", "Pc", "Omega", "Rho", "M" - ! defined components - do i = 1, def_nc - write(*, fmt=num_fmt_1) def_name(i), z(i), def_mw(i), tc_def(i), pc_def(i), om_def(i) - end do - ! - i_prev = scn(1) - 1 - do i = 1, scn_nc_new - write(*, fmt=num_fmt_2) 'C' // str(i_prev + i), z(i+def_nc), mw(i), tc(i), pc(i), omega(i), rho(i), m_funtion(i) - end do - ! - do i = 1 + scn_nc_new, size(mw) - write(*, fmt=num_fmt_2) 'ps' // str(i - scn_nc_new), z(i+def_nc), mw(i), rho(i), tc(i), pc(i), omega(i), m_funtion(i) - end do - - end associate - end subroutine write_result - + @@ -890,7 +828,8 @@ type(FluidDataOut) function characterize(file, mw_source, method, fix_C, eos) & - fluid = data_from_file(file) + fluid = data_from_file(file) + characterization%input_data = fluid allocate(characterization%scn_z(fluid%scn_nc)) allocate(characterization%log_scn_z(fluid%scn_nc)) allocate(characterization%scn_mw(fluid%scn_nc)) @@ -908,12 +847,13 @@ type(FluidDataOut) function characterize(file, mw_source, method, fix_C, eos) & allocate(characterization%critical_pressure(0)) allocate(characterization%acentric_factor(0)) allocate(characterization%m_funtion(0)) + call get_c_or_m_plus(fluid=fluid, mw_source=mw_source, method=method, fix_C=fix_C, characterization= characterization) call density_funtion(fluid=fluid, mw_source=mw_source, characterization=characterization) call lump(fluid=fluid, characterization=characterization) call get_critical_constants(fluid= fluid, characterization=characterization, eos=eos) - call write_result(fluid=fluid,characterization=characterization) + end function characterize diff --git a/src/types.f90 b/src/types.f90 index 8364ed9..38daa35 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -7,20 +7,20 @@ module dtypes integer :: scn_nc !! number of single cuts being considered in the oil integer :: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified number of pseudos integer :: numbers_ps !! number of pseudos in which the scn fractions grouped. - character(len=:), allocatable :: filename - integer, allocatable :: scn (:) !! set of singles cuts being considered in the oil - character(len=15), allocatable :: def_components (:) !! set of defined components being considered in the oil - character(len=15) :: scn_plus !! name of residual fraction - real(pr), allocatable :: def_comp_z(:) !! set of corresponding mole fractions of defined components - real(pr), allocatable :: scn_z(:) !! set of corresponding mole fractions of scn cuts - real(pr), allocatable :: def_comp_mw(:) !! set of corresponding molecular weights of defined components + character(len=:), allocatable :: filename + integer, allocatable :: scn (:) !! set of singles cuts being considered in the oil + character(len=15), allocatable :: def_components (:) !! set of defined components being considered in the oil + character(len=15) :: scn_plus !! name of residual fraction + real(pr), allocatable :: def_comp_z(:) !! set of corresponding mole fractions of defined components + real(pr), allocatable :: scn_z(:) !! set of corresponding mole fractions of scn cuts + real(pr), allocatable :: def_comp_mw(:) !! set of corresponding molecular weights of defined components real(pr), allocatable :: scn_mw(:) !! set of corresponding molecular weights of scn cuts real(pr), allocatable:: product_z_mw_def_comp(:) !! product between composition and molecular weight of defind components real(pr), allocatable:: product_z_mw_scn(:) !! product between composition and molecular weight of scn fractions real(pr) :: sum_z_mw_i !! sum of the product between composition and molecular weight of the fluid's compounds real(pr), allocatable :: w(:) !! mass fractions of the fluid's compounds real(pr) :: plus_z !! composition of residual fraction - real(pr) :: plus_mw !! molecular weight of residual fraction + real(pr) :: plus_mw !! molecular weight of residual fraction real(pr) :: product_z_mw_plus !! product between composition and molecular weight of residual fraction real(pr), allocatable :: def_comp_w(:) !! mass fractions of the defined compounds real(pr), allocatable :: scn_w(:) !! !! mass fractions of the scn-s compounds @@ -31,6 +31,7 @@ module dtypes type :: FluidDataOut ! incluir aqui las variables que quiero que salgan como salida + type(FluidData) :: input_data real(pr), allocatable :: scn_z(:) !! set of corresponding mole fractions of scn cuts calculated real(pr), allocatable :: log_scn_z(:) !! set of corresponding mole fractions of scn cuts calculated real(pr) :: plus_mw !! molecular weight of residual fraction @@ -61,7 +62,7 @@ module dtypes integer :: last_C integer :: i_last integer :: last - integer :: scn_nc_new + integer :: scn_nc_new real(pr), allocatable :: scn_i(:) real(pr), allocatable :: plus6_density(:) real(pr) :: plus_w @@ -74,10 +75,86 @@ module dtypes real(pr), allocatable :: critical_pressure(:) real(pr), allocatable :: acentric_factor(:) real(pr), allocatable :: m_funtion (:) - + + contains + private + procedure, pass :: write => write_result + generic, public :: write (FORMATTED) => write end type FluidDataOut - +contains + + subroutine write_result(characterization,unit,iotype,v_list,iostat,iomsg) + use ftools__io, only: str + use critical_parameters + + implicit none + class(FluidDataOut), intent(in) :: characterization + integer :: i, i_prev + integer, intent(in) :: unit + integer, intent(out) :: iostat + character(*), optional, intent(in) :: iotype + character(*), optional, intent(inout) :: iomsg + integer, optional, intent(in) :: v_list(:) + character(len=*), parameter :: str_fmt_1 = "(A4,7(A15,2x),/)" + character(len=*), parameter :: num_fmt_1 = "(A3,1x,*(E15.5,2x),/)" + character(len=*), parameter :: num_fmt_2 = "(A3,1x,*(E15.5,2x),/)" + + associate(& + def_nc => characterization%input_data%def_comp_nc, & + def_name => characterization%input_data%def_components, & + def_mw => characterization%input_data%def_comp_mw, & + z => characterization%mol_fraction, & + mw => characterization%lumped_mw, & + rho => characterization%lumped_densities, & + scn => characterization%input_data%scn, & + scn_nc_new => characterization%scn_nc_new, & + tc => characterization%critical_temperature, & + pc => characterization%critical_pressure, & + omega => characterization%acentric_factor, & + m_funtion => characterization%m_funtion & + ) + + print*, "----------------------------------------------------------------" + print*, "Best feasible regresion parameters" + print*, "" + print*, "Init_BFR:", characterization%n_init + print*, "A:", characterization%a," ","B:", characterization%b + print*, "C:", characterization%C + print*, "MW+:", characterization%plus_mw + print*, "Cmax:", characterization%c_max + print*, "" + print*, "Density funtion parameters parameters" + print*, "" + print*, "ad:", characterization%a_d," ","bd:", characterization%b_d + print*, "----------------------------------------------------------------" + print*, "" + print*, "compositional result" + print*, "" + + + write(unit, fmt=str_fmt_1, iostat = iostat) "Comp", "Z", "Mw", "Tc", & + "Pc", "Omega", "Rho", "M" + ! defined components + do i = 1, def_nc + write(unit, fmt=num_fmt_1, iostat = iostat) def_name(i), z(i), & + def_mw(i), tc_def(i), pc_def(i), om_def(i) + end do + ! + i_prev = scn(1) - 1 + do i = 1, scn_nc_new + write(unit, fmt=num_fmt_2, iostat = iostat) 'C' // str(i_prev + i), & + z(i+def_nc), mw(i), tc(i), pc(i), omega(i), rho(i), m_funtion(i) + end do + ! + do i = 1 + scn_nc_new, size(mw) + write(unit, fmt=num_fmt_2, iostat = iostat) 'ps' // str(i - scn_nc_new), & + z(i+def_nc), mw(i), rho(i), tc(i), pc(i), omega(i), m_funtion(i) + end do + + end associate + end subroutine write_result + end module dtypes From 05cbd2a7f2e093e29b26a91655e3dd138f70dd86 Mon Sep 17 00:00:00 2001 From: Oscar Theran Date: Mon, 4 Nov 2024 10:12:44 -0300 Subject: [PATCH 07/16] start a new refactoring 04/11/2024 --- app/main.f90 | 8 ++--- app/optimize.f90 | 65 +++++++++++++++++++++++++++++++++++++++++ app/pedersen.f90 | 76 ++++++++++++++++++++++++++++++++++++++++++++++++ fpm.toml | 4 ++- src/routines.f90 | 10 ------- src/types.f90 | 47 ++++++++++++++++++------------ 6 files changed, 174 insertions(+), 36 deletions(-) create mode 100644 app/optimize.f90 create mode 100644 app/pedersen.f90 diff --git a/app/main.f90 b/app/main.f90 index 75dc1aa..a15948a 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -5,18 +5,14 @@ program main implicit none type(FluidDataOut) :: prueba - !type(FluidData) :: fluid - integer :: i - - prueba = characterize(file='oil1.nml', mw_source="calculated", method = "plus_mw", fix_C=.false., eos='PR') + prueba = characterize(file='YPF2.nml', mw_source="calculated", method = "plus_mw", fix_C=.false., eos='PR') write(*, *) prueba + print*, prueba%a , prueba%b , prueba%n_init !do i = 1, size(prueba%mol_fraction) ! print*, prueba%mol_fraction(i) !end do - - end program main diff --git a/app/optimize.f90 b/app/optimize.f90 new file mode 100644 index 0000000..ceeaacc --- /dev/null +++ b/app/optimize.f90 @@ -0,0 +1,65 @@ +module my_objective + use ForTimize, only: pr + + type :: ExpData + real(pr), allocatable :: x(:) + real(pr), allocatable :: y(:) + end type + +contains + + subroutine foo(x, F, dF, data) + real(pr), intent(in) :: x(:) + !! Vector of parameters to optimize + real(pr), intent(out) :: F + !! Value of the objetive function after optimization. + real(pr), optional, intent(out) :: dF(:) !! + !! Optional gradient. + class(*), optional, intent(in out) :: data + !! Special data that the function could use. + + real(pr) :: a, b, c + + a = x(1) + b = x(2) + c = x(3) + + F = 0 + + select type(data) + type is (ExpData) + F = sum( (data%y - (quad(data%x, a, b, c)))**2 ) + print *, F, data%x + end select + end subroutine + + elemental function quad(x, a, b, c) + real(pr), intent(in) :: x, a, b, c + real(pr) :: quad + + quad = a*x**2 + b*x + c + end function +end module + + +program oscar + use ForTimize, only: pr, minimize + use my_objective, only: foo, ExpData + + real(pr) :: x(3), F + type(ExpData) :: exp + + exp%x = [2, 5, 7] + exp%y = [4, 22, 50] + + ! Initial guess + x = [1, 4, 3] + + ! Minimize uses the Nelder-Mead algorithm as a default + call minimize(foo, x, F, data=exp) + ! call foo(x, F, data=exp) + + ! Print results + print *, x + print *, F +end program oscar diff --git a/app/pedersen.f90 b/app/pedersen.f90 new file mode 100644 index 0000000..51d939b --- /dev/null +++ b/app/pedersen.f90 @@ -0,0 +1,76 @@ +module my_objective1 + use ForTimize, only: pr + + + type :: ExpData + real(pr), allocatable :: x(:) + real(pr), allocatable :: zplus(:) + real(pr), allocatable :: mplus(:) + end type ExpData + +contains + + subroutine foo(x, F, dF, data) + real(pr), intent(in) :: x(:) + !! Vector of parameters to optimize + real(pr), intent(out) :: F + !! Value of the objetive function after optimization. + real(pr), optional, intent(out) :: dF(:) !! + !! Optional gradient. + class(*), optional, intent(in out) :: data + !! Special data that the function could use. + + real(pr) :: a, b + integer :: i + integer, parameter :: n = 181 + real(pr) :: cn(n) + real(pr) :: z_i(n) + real(pr) :: m_i(n) + real(pr) :: zplus_cal + real(pr) :: mplus_cal + + + do i = 1, n + cn(i) = 20.0d0 + i - 1 + end do + + a = x(1) + b = x(2) + + z_i = exp(a*cn+b) + m_i = 14*cn-4 + zplus_cal = sum(z_i) + mplus_cal = sum(z_i*m_i)/zplus_cal + + F = 0 + + select type(data) + type is (ExpData) + F = sum((data%zplus - (zplus_cal))**2 + (data%mplus - (mplus_cal))**2) + end select + end subroutine foo + +end module my_objective1 + + +program oscar + use ForTimize, only: pr, minimize + use my_objective1, only: foo, ExpData + + real(pr) :: x(2), F + type(ExpData) :: exp + + exp%zplus = [0.2133] + exp%mplus = [442.1] + + ! Initial guess + x = [0, 0] + + ! Minimize uses the Nelder-Mead algorithm as a default + call minimize(foo, x, F, data=exp) + + ! Print results + print *, x + print *, F + !print*, cn +end program oscar diff --git a/fpm.toml b/fpm.toml index 77a72ab..d4d70ef 100644 --- a/fpm.toml +++ b/fpm.toml @@ -14,4 +14,6 @@ library = false [fortran] implicit-typing = false implicit-external = false -source-form = "free" \ No newline at end of file +source-form = "free" +[dependencies] +ForTimize = {git='https://github.com/fedebenelli/ForTimize'} \ No newline at end of file diff --git a/src/routines.f90 b/src/routines.f90 index cc25321..107c7d8 100644 --- a/src/routines.f90 +++ b/src/routines.f90 @@ -806,16 +806,6 @@ subroutine get_critical_constants(fluid, characterization, eos) end subroutine get_critical_constants - - - - - - - - - - type(FluidDataOut) function characterize(file, mw_source, method, fix_C, eos) & result(characterization) diff --git a/src/types.f90 b/src/types.f90 index 38daa35..9f77624 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -97,8 +97,9 @@ subroutine write_result(characterization,unit,iotype,v_list,iostat,iomsg) character(*), optional, intent(inout) :: iomsg integer, optional, intent(in) :: v_list(:) character(len=*), parameter :: str_fmt_1 = "(A4,7(A15,2x),/)" - character(len=*), parameter :: num_fmt_1 = "(A3,1x,*(E15.5,2x),/)" - character(len=*), parameter :: num_fmt_2 = "(A3,1x,*(E15.5,2x),/)" + character(len=*), parameter :: num_fmt_1 = "(A3,1x,5(E15.5,2x),/)" + character(len=*), parameter :: num_fmt_2 = "(A3,1x,7(E15.5,2x),/)" + character(len=*), parameter :: str_fmt_2 = "(130('-'),/),/)" associate(& def_nc => characterization%input_data%def_comp_nc, & @@ -115,22 +116,30 @@ subroutine write_result(characterization,unit,iotype,v_list,iostat,iomsg) m_funtion => characterization%m_funtion & ) - print*, "----------------------------------------------------------------" - print*, "Best feasible regresion parameters" - print*, "" - print*, "Init_BFR:", characterization%n_init - print*, "A:", characterization%a," ","B:", characterization%b - print*, "C:", characterization%C - print*, "MW+:", characterization%plus_mw - print*, "Cmax:", characterization%c_max - print*, "" - print*, "Density funtion parameters parameters" - print*, "" - print*, "ad:", characterization%a_d," ","bd:", characterization%b_d - print*, "----------------------------------------------------------------" - print*, "" - print*, "compositional result" - print*, "" + + + write(unit, fmt=str_fmt_2, iostat = iostat) + !write(*,"(A64,/)") "----------------------------------------------------------------" + !write(*,*) "Best feasible regresion parameters" + !write(unit, fmt=str_fmt_1, iostat = iostat) "" + !write(unit, fmt=num_fmt_2, iostat = iostat) "Init_BFR:", characterization%n_init + !write(*,*) "A:", characterization%a," ","B:", characterization%b + !write(*,*) "C:", characterization%C + !write(*,*) "MW+:", characterization%plus_mw + !write(*,*) "Cmax:", characterization%c_max + !write(*,*) "" + !write(*,*) "Density funtion parameters parameters" + !write(*,*) "" + !write(*,*) "ad:", characterization%a_d," ","bd:", characterization%b_d + !write(*,*) "----------------------------------------------------------------" + !write(*,*) "" + !write(*,*) "compositional result" + !write(*,*) "" + + + + + write(unit, fmt=str_fmt_1, iostat = iostat) "Comp", "Z", "Mw", "Tc", & @@ -143,7 +152,7 @@ subroutine write_result(characterization,unit,iotype,v_list,iostat,iomsg) ! i_prev = scn(1) - 1 do i = 1, scn_nc_new - write(unit, fmt=num_fmt_2, iostat = iostat) 'C' // str(i_prev + i), & + write(unit, fmt=num_fmt_2, iostat = iostat) adjustl('C' // str(i_prev + i)), & z(i+def_nc), mw(i), tc(i), pc(i), omega(i), rho(i), m_funtion(i) end do ! From 180da8e6919008fc5c0077cb64e98fcc56160285 Mon Sep 17 00:00:00 2001 From: oscartb25 Date: Tue, 7 Jan 2025 21:08:26 -0500 Subject: [PATCH 08/16] 07 enero 2025 --- .github/workflows/docs.yml | 0 .gitignore | 16 +- PVT5.nml | 78 ++ README.md | 0 YPF1.nml | 0 app/main.f90 | 14 +- app/my_objective.mod | Bin 0 -> 3632 bytes app/optimize.f90 | 0 app/pedersen.f90 | 6 +- example/demo.f90 | 0 fpm.toml | 0 input.nml | 0 oil1.nml | 0 src/bfr_routines/bfr_routines.f90 | 226 +++++ src/characterization/characterization.f90 | 58 ++ src/constans.f90 | 0 src/critical_parameters.f90 | 88 -- .../critical_parameters.f90 | 72 ++ .../defined_critical_param.f90 | 82 ++ src/density/density_funtion.f90 | 160 ++++ src/input_data/data_from_input.mod | Bin 0 -> 3919 bytes src/{ => input_data}/data_input.f90 | 103 ++- src/lumping/lumping_method.f90 | 184 ++++ src/molecular_weigth/mw_calculations.f90 | 268 ++++++ src/routines.f90 | 852 ------------------ src/tools/io.f90 | 0 src/types.f90 | 20 +- 27 files changed, 1248 insertions(+), 979 deletions(-) mode change 100644 => 100755 .github/workflows/docs.yml mode change 100644 => 100755 .gitignore create mode 100755 PVT5.nml mode change 100644 => 100755 README.md mode change 100644 => 100755 YPF1.nml mode change 100644 => 100755 app/main.f90 create mode 100755 app/my_objective.mod mode change 100644 => 100755 app/optimize.f90 mode change 100644 => 100755 app/pedersen.f90 mode change 100644 => 100755 example/demo.f90 mode change 100644 => 100755 fpm.toml mode change 100644 => 100755 input.nml mode change 100644 => 100755 oil1.nml create mode 100755 src/bfr_routines/bfr_routines.f90 create mode 100755 src/characterization/characterization.f90 mode change 100644 => 100755 src/constans.f90 delete mode 100644 src/critical_parameters.f90 create mode 100755 src/critical_parameters/critical_parameters.f90 create mode 100755 src/critical_parameters/defined_critical_param.f90 create mode 100755 src/density/density_funtion.f90 create mode 100755 src/input_data/data_from_input.mod rename src/{ => input_data}/data_input.f90 (71%) mode change 100644 => 100755 create mode 100755 src/lumping/lumping_method.f90 create mode 100755 src/molecular_weigth/mw_calculations.f90 delete mode 100644 src/routines.f90 mode change 100644 => 100755 src/tools/io.f90 mode change 100644 => 100755 src/types.f90 diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml old mode 100644 new mode 100755 diff --git a/.gitignore b/.gitignore old mode 100644 new mode 100755 index 32b8afb..259de11 --- a/.gitignore +++ b/.gitignore @@ -86,4 +86,18 @@ infiles/YPF10 YPF10 fig77 fpm.toml - +src/tools/ftools__io.mod +YPF1.nml +YPF2.nml +YPF1.nml +YPF1.nml +fpm.toml +fpm.toml +src/molecular_weigth/molecular_weigth.mod +.gitignore +src/bfr_routines/bfr_routines.mod +src/characterization/characterization.mod +src/critical_parameters/critical_parameters.mod +.gitignore +src/density/density.mod +src/lumping/lumping.mod diff --git a/PVT5.nml b/PVT5.nml new file mode 100755 index 0000000..7ffb5b0 --- /dev/null +++ b/PVT5.nml @@ -0,0 +1,78 @@ +! Namelist based input file +! ============================================================================== +! - PVT5 +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 24 + ! number of single carbon number components + scn_nc_ps = 30 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped + density_setup = .true. + ! Logical variable to choose the density calculation method + number_plus_density = 6 + ! CN from which we considered ....... +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 + ! names of scn fractions + scn_plus= "C30+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.0026 0.0053 0.5403 0.1260 0.0764 0.0139 0.0341 0.0136 0.0168 + ! compositions of defined components + scn_z = 0.0213 0.0183 0.0205 0.0152 0.0112 0.0098 0.0079 0.0072 0.0055 0.0040 0.0031 0.0027 0.0025 0.0022 0.0017 0.0015 0.0013 0.0012 0.0010 0.0010 0.0008 0.0007 0.0006 0.0006 + ! compositions of single carbon numbers + plus_z = 0.0057 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.1711 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.1492 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.0512 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.0057 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.2 72.2 + ! molecular weights of defined components + scn_mw = 86 95.9 108.2 124 138.3 147.5 161.1 174.1 189.5 207.7 223.1 236.6 247.8 262 276.8 297.8 310.3 312.2 325.3 353.6 360.2 377.3 384.0 401.5 + ! molecular weights of single carbon numbers + plus_mw = 579.2 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 157.21 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 167.73 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 265.89 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 579.2 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + ! density of single carbon numbers + plus_density = 0.9837 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.8102 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.8214 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.8732 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.9837 + ! density of residual fraction CN 30 plus +/ + + diff --git a/README.md b/README.md old mode 100644 new mode 100755 diff --git a/YPF1.nml b/YPF1.nml old mode 100644 new mode 100755 diff --git a/app/main.f90 b/app/main.f90 old mode 100644 new mode 100755 index a15948a..98886dc --- a/app/main.f90 +++ b/app/main.f90 @@ -1,18 +1,22 @@ program main - use data_from_input, only: data_from_file, FluidData - use routines + !use data_from_input, only: data_from_file, FluidData + use characterization implicit none + type(FluidDataOut) :: prueba + integer :: i - prueba = characterize(file='YPF2.nml', mw_source="calculated", method = "plus_mw", fix_C=.false., eos='PR') + prueba = characterize(file='PVT5.nml', mw_source="experimental", method = "plus_mw", pho_method = 3 , fix_C=.true., eos='PR') write(*, *) prueba - print*, prueba%a , prueba%b , prueba%n_init - + print*, prueba%C, prueba%a, prueba%b, prueba%plus_mw, prueba%plus_z + + !print*, prueba%a , prueba%b , prueba%n_init, prueba%input_data%number_plus_density !do i = 1, size(prueba%mol_fraction) ! print*, prueba%mol_fraction(i) !end do + end program main diff --git a/app/my_objective.mod b/app/my_objective.mod new file mode 100755 index 0000000000000000000000000000000000000000..6ad6d7c2e8eb644e274bd09c78b6f7c852d9d1a7 GIT binary patch literal 3632 zcmV-04$tu)iwFP!000001KnL)bJ|E2em}pW-@p&ErPiIdC7^hz!9Zc~O!h@ZYY%-IfVz7|DbHDWI&h5wPd@=7&$d~QY*O!QVi+205#`=v^qpjL<&zu0&+yglMBiVZR)z5@--td(zLEU*_g|>lVZTHy9(2B}zK81r zsvv&GaNGuSVjpjlzJC(}>K=|{&=`CQ+$WuXBme(k3$j|W;v%b+|Z@?;rP}@1S2>YF6TknwcJgI z%W;q}WEc#_1#a9gKg_4g#b^?s%KX;Mq9Kjj@qHByUKorTQs`RfLM*lub>u0-O8>F|@3#G~w?L@Q}Y zE10onn3R|%F)d=+LPAYc>D!|KpQ8O^#QYdf<$u?QKSPx0`rC@}VxNmUzZ_WrGOce+ z1enU2?lpve)BBm?|H^xGR*yy&|Lrs&nf^}}|1;pdQbYm;cB`lx)S_+VHn41q`aR3b zv#n#Yr4yo9lsY*?Sw`Zfq3EWm+e3sVy#WjTbQnIqtoHj*#Yu;G<{ut9&?@)W40I-! z`;*aqI2a?jj4PIj*pr&23+cTrMIWFYZD&(<)kAeoGJn zgHqD2Gy@QXarE>SdV19J`PJ8!FrgEh8Z>Sn+hpq{BEw; zAtNL_=#m2MLX=TdQj8L*l4kmGsX!%BQn377OEVC?HwIy&0}o%K4Hsil=XvfJgIUmm z`NC0uTr`zRLDbe#AlEco>_S<-7iQ^yS^EmvU`s`#c^1fPoXxpKwan<0^6V&VNTHOZ zZd+B&v+LX=oKWv(ng_M+QA6e+^NxnByirVmnbfl)piHM?opZL_-)8F3wOwL+#P*5Z zBaR`ib4(f8;_9Xx-*l^+vsH2qYff|>lQO2zIDpt+Etx%-xS1Uq}`O3 zo~1L$yAao_O4E%k2eb>6#&TIAO}D4!K+e+j1~CRtH~l=RX`MTei{~QfKLW|xU9B1w zr~n}qg*7)19q{+s_CHA4|0n1Z*Zn88-G3A?GU5k^4=VU%%M((33O#sLe5eOHtsj9G zIf@({*Rf64L&tIDeh?>!nGZ)eN>?C zdUM#VH~aPDO>>6dNYBURWUxRRA;CC6b6KEy&hS17l8l{DWoJ$U!Yn539t$s+3q6PH z&H51X1`sy5#)tlV6fZ5R$Ypewlp(*|_aB0ibO@T?ZLmmqn`KRcd_~r@x&fa>sPNli z^*K!h+m;?FRZ7}cr79r>9?-2N0BaMl8iHs>P+09H7BU?PVzEBS5R)q6)Hi@4VlzW5 zZRIFX(Ndc#R6uRZq=pPjM21RJU$>=yF}c}#ad{I{avm^T zGqjjblskun-6E{!bHywQGy;bKwZdA(V^6se5Y9zjgyDK$`@WJsu(EXU)f zxt|zLDQQ$A7kY8xAFWvR=IG@T#8qn{+Ys*ID*Gqv0imk*;py1tJ6%aMDNCUd=NE?dG z9p&m&?iH)~BWh4ClGKiEoSmieNK(&FwMbHNP)B-A?MMQ3qSML{lu`gC$(^t?P;#}l zox9pL6=&OVo6P-E!^-{d8%)vbC;Z2m;?xaO!+xKHpBl~;Wc<`{#VP-LW&G6eL@9rr z2a``eypK7v-3`7lrKd(u$rGQj!YoU0g*kn*Csr(BDBI6b>xm_GjF#v?EIPGfj29i$ zQBv4yzzfe5$PTpVkS6sY-)^9q)Ps1Ng-%8<-`Y{m=*3~?5*@*l3j)ffjIOyNNH_R& z;24pp^*)6GrM9p#XoyJXtqdr2B&ET{MzI?z1L}$@14_NiB-NZJDE6R))xLhjlBCj9 z60xW}Vx2Tw$!7s#I}!ny|T$H1+SmD0yK-!3wTMMOo>WQsxo8QwkK64 z5L))6$`K6W2Q_a9DKR>@+9DWt|!h*L7*UNA z$a_PITbTsaT#+K+DRScp?Nd0KjxZCBx3W_O5*Uo-$n51^Fbwbi(Vqoc8~Rg47l_A0 z@ra30AtfgO8VC?P8lCGgZnbCSbjIu{8w@2?Wp$O5HSw(`5MLY`3NaYHr$W`5e z)yS2VQ_o^bE-hD97d@TuY|cC~>FV$kr2N!lCvWqt&_rF#@Smj~^xOl)#m^bld0bz8ycW$EDjn&h(4KT-x*Z9C{9Gv2{`)$0}HhQ=32j6e*BnbN14Xqt|0Ej>FkYk(i6yi(9mpR;+Ba z7scpuw3l9XB28|>Eymp;(I&}pGw~A@kDKB5_4-9{aTvZ9aGp~H#~KY<=n8GOlh8r6 zQ`Ae&PuafyIAQMt7L=Z{^Qu(8T>nYY5yEV~toB&dwY6gqAbMo2u4` znbO!1q*&Nc*TPqP)+S`j;IEO$8d(kQ7w_CAoFPn=ifbgy<(s=VH#tX=f9KS=Q?8N+!YVPe}^?RKjl87QYszm*4-xgqo* zi_z_Jy1Z9KnyqddoR6RPC%=joyG&H-ABQO1A9m}fNP@p6hKf%J$_aVskCy?DTrznC zO{;ROC$-+JG3}}2^f{4^r?;a4{vvIi9p5)1bt=kWzl|~owHDNM!IkKOeQmnn$#mgj zNZaZ{ujo7DLhOvbQ#4C0h6C|Juu*iKWjJSWow|t3vaz^rT8b@>gSlDTZB4crD zjz4)9V__#I3XSD*NLv~UI|Pw87S}m*qM@mpUeQ>}aL$ouXd>!WjOA1VR!uWbGZuEp zqR?3UD>D{$n4#EME{3$7v9Nz|EkECpoz+@B^3Nd0lK_h+x8cV=G#`wHm-VU#j~T=_ zpmcf{ewLr;RnxY`>r7~S02TYT1@)&n;U(OLHV!j?kWvB~b(cTRq87DnL+f2~_;ri> zx)is=O~o1{o|27MVA; y=a*x+b + t1=0.;t2=0.;t3=0.;t4=0.; + + do i=1,n + t1=t1+x(i)*y(i) + t2=t2+x(i) + t3=t3+y(i) + t4=t4+x(i)**2 + end do + + a=(n*t1-t2*t3)/(n*t4-t2**2) + b=(t3-a*t2)/n + !coefficient calculation of correlations r2 + aux1=0.;aux2=0.;aux3=0.;aux4=0.;aux5=0.; + + do i=1, n + aux1= aux1 + x(i)*y(i) + aux2= aux2 + x(i) + aux3= aux3 + y(i) + aux4= aux4 + x(i)**2 + aux5= aux5 + y(i)**2 + end do + + r2=(aux1-aux2*aux3/n)**2 /((aux4-aux2**2/n)*(aux5-aux3**2/n)) + + end subroutine Linear_Regression + + subroutine Best_Linear_Regression(scn_nc, scn, scn_z, plus_z, a, b, r2, & + n_init, c_max_blr) + !! This subroutine calculates the best regression line for an fluid. + implicit none + + integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid + integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid + integer, intent(out) :: c_max_blr !! output CN at which plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) + real(pr), intent(in) :: scn_z(:) !! set of corresponding mole fractions of scn cuts + real(pr), intent(in) :: plus_z !! composition of residual fraction from input file + real(pr), intent(out) :: a !! output real variable. Slope of the best regression line. + real(pr), intent(out) :: b !! output real variable. Intercept of the best regression line. + real(pr), intent(out) :: r2 !! output real variable. Square correlation coefficient. + integer, intent(out) :: n_init !! minimum carbon number obtained from the best linear regression + integer :: i, j, k , k_old, n_best, x_aux ! internal variables + real(pr), dimension(scn_nc) :: x_blr, y_blr ! x_blr: carbon number vector, y_blr: logarithm of mole fraction vector; for each step + real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux ! auxiliary variables + + k=5 + r2=0.0001_pr + r2_old=0.00001_pr + r2_best=0.0001_pr + + do while (r2.gt.r2_old.or.r2_old.lt.0.9) + k_old=k + r2_old=r2 + a_old=a + b_old=b + + if (r2.gt.r2_best) then + r2_best=r2 + a_best=a + b_best=b + n_best=scn(scn_nc-k+2) + end if + + if (k.gt.scn_nc) then + if (r2.gt.r2_best)then + n_init=scn(scn_nc-k+2) + go to 22 + else + r2=r2_best + a=a_best + b=b_best + n_init=n_best + go to 22 + end if + end if + + j=1 + x_blr=0. + y_blr=0. + + do i=scn_nc-k+1, scn_nc + x_blr(j)=scn(i) + y_blr(j)=scn_z(i) + j=j+1 + end do + + call Linear_Regression(x_blr(:k), y_blr(:k), a, b, r2) + k=k+1 + end do + + r2=r2_old + a=a_old + b=b_old + n_init=scn(scn_nc-k_old+2) + + 22 continue + + z_sum = 0d0 + x_aux=scn(scn_nc) + + do while (z_sum.lt.plus_z.and.x_aux<300) + x_aux=x_aux+1 + z_aux= exp(a*x_aux+b) + z_sum=z_sum+z_aux + end do + + c_max_blr=x_aux + + !print*, a, b, c_max_blr, r2, n_init + + end subroutine Best_Linear_Regression + + subroutine LimitLine(scn_nc, scn, plus_z, a_blr, b_blr, a_lim, b_lim, & + c_max_lim, half) + !! This subroutine obtains the limit line constants for an fluid. + implicit none + + integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid. + integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid. + integer, intent(out) :: c_max_lim !! maximum carbon number obtained for limit line. + !! output CN at which plus_z is reached, as the summation of single z(i) from the limit distribution (blr). + real(pr), intent(in) :: a_blr !! input constant from the Best linear regression. + real(pr), intent(in) :: b_blr !! input constant from the Best linear regression. + real(pr), intent(in) :: plus_z !! composition of residual fraction from input file. + real(pr), intent(out) :: a_lim !! output real variable. Slope of the limit line. + real(pr), intent(out) :: b_lim !! output real variable. Intercept of the limit line. + real(pr), intent(out) :: half ! variable defined in the numerical method to converge the area under the curve of the function. + real(pr) :: z_lim, cross_cn, z_cross, z_aux + integer :: x_aux ! internal variable + + !A and B limit calculation. + z_lim = 0.0_pr + half = 0.506_pr ! modified by oscar, the variable half was removed from the argument + + do while (z_lim.lt.plus_z) + half = half + 0.002d0 + cross_cn = scn(scn_nc)+half ! typically 19.508 in first try + z_cross = exp((a_blr*cross_cn) + b_blr) + a_lim = -(z_cross)/plus_z + b_lim = log(z_cross)-(a_lim*cross_cn) + x_aux=scn(scn_nc) + z_lim = 0._pr + do while (z_lim.lt.plus_z.and.x_aux<300) + x_aux=x_aux+1 + z_aux= exp(a_lim*x_aux+b_lim) + z_lim=z_lim+z_aux + end do + continue + end do + + c_max_lim = x_aux + + end subroutine LimitLine + + subroutine Line_C60_max (scn_nc, scn, plus_z, a_blr, b_blr, half, a_60, & + b_60, c_max_60) + !! This subroutine obtains the Cmax60 line constants for an fluid. + + implicit none + integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid + integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid + integer, intent(out) :: c_max_60 !!output CN at which Zp is reached, as the summation of single z(i) from the Cmax60 distribution. + real(pr), intent(in) :: a_blr !! input constant from the Best linear regression + real(pr), intent(in) :: b_blr !! input constant from the Best linear regression + real(pr), intent(in) :: plus_z !! composition of residual fraction from input file + real(pr), intent(out) :: a_60 !! output real variable. Slope of the limit line + real(pr), intent(out) :: b_60 !! output real variable. Intercept of the C60max line. + real(pr), intent(in) :: half !variable defined in the numerical method to converge the area under the curve of the function + integer :: x_aux ! auxiliary variable + real(pr) :: z_sum, cross_cn, z_cross, z_aux, F_tol, a_tol, var_range ! internal variables + real(pr) :: a_old, Zp_60, F, dF_dA ! internal variables + + cross_cn = scn(scn_nc)+half ! typically 19.51 + z_cross= exp((a_blr*cross_cn)+b_blr) + F_tol=1_pr + a_tol=1_pr + a_60=a_blr + var_range = 60.5_pr - cross_cn + + do while (a_tol.gt.1e-6_pr.or.F_tol.gt.1e-6_pr) + a_old=a_60 + Zp_60 = (exp(a_60*var_range+log(z_cross))-z_cross)/a_60 + F = plus_z - Zp_60 + dF_dA = - (var_range*(a_60*Zp_60+z_cross)-Zp_60) / a_60 ! modified by oscar, error in the derivative in previous version is corrected + a_60=a_old-(F/dF_dA) + a_tol=abs(a_old-a_60) + F_tol=abs(F) + end do + + b_60=log(z_cross)-a_60*cross_cn + z_sum=0.0_pr + x_aux=scn(scn_nc) + + do while (z_sum.lt.plus_z) + x_aux=x_aux+1 + z_aux= exp(a_60*x_aux+b_60) + z_sum=z_sum+z_aux + end do + + c_max_60=x_aux + + !print*, a_60,b_60,c_max_60 + end subroutine Line_C60_max + + +end module bfr_routines + diff --git a/src/characterization/characterization.f90 b/src/characterization/characterization.f90 new file mode 100755 index 0000000..9d42ae0 --- /dev/null +++ b/src/characterization/characterization.f90 @@ -0,0 +1,58 @@ +module characterization + + use constants + use dtypes, only: FluidData, FluidDataOut + use data_from_input, only: data_from_file, FluidData + use molecular_weigth + use density + use lumping + use critical_parameters + + +contains + + type(FluidDataOut) function characterize(file, mw_source, method, pho_method, fix_C, eos) & + result(characterization) + + type(FluidData) :: fluid + character(len=*), intent(in) :: file !! file name + character(len=*), intent(in), optional :: method !! plus_mw or global_mw + character(len=*), intent(in) :: mw_source + logical, intent(in), optional :: fix_C + character(len=*), intent(in) :: eos + integer, intent(in) :: pho_method !! selector method of density funtion. + + + + fluid = data_from_file(file) + characterization%input_data = fluid + allocate(characterization%scn_z(fluid%scn_nc)) + allocate(characterization%log_scn_z(fluid%scn_nc)) + allocate(characterization%scn_mw(fluid%scn_nc)) + allocate(characterization%carbon_number_plus(0)) + allocate(characterization%plus_z_i(0)) + allocate(characterization%product_z_mw_plus_i(0)) + allocate(characterization%scn_zm(fluid%scn_nc)) + allocate(characterization%scn_i(0)) + allocate(characterization%plus6_density(0)) + allocate(characterization%mol_fraction(0)) + allocate(characterization%lumped_z(0)) + allocate(characterization%lumped_mw(0)) + allocate(characterization%lumped_densities(0)) + allocate(characterization%critical_temperature(0)) + allocate(characterization%critical_pressure(0)) + allocate(characterization%acentric_factor(0)) + allocate(characterization%m_funtion(0)) + + + call get_c_or_m_plus(fluid=fluid, mw_source=mw_source, method=method, fix_C=fix_C, characterization= characterization) + call density_funtion(fluid=fluid, mw_source=mw_source, pho_method= pho_method, characterization=characterization) + call lump(fluid=fluid, characterization=characterization) + call get_critical_constants(fluid= fluid, characterization=characterization, eos=eos) + + + end function characterize + +end module characterization + + diff --git a/src/constans.f90 b/src/constans.f90 old mode 100644 new mode 100755 diff --git a/src/critical_parameters.f90 b/src/critical_parameters.f90 deleted file mode 100644 index a6baee5..0000000 --- a/src/critical_parameters.f90 +++ /dev/null @@ -1,88 +0,0 @@ -module critical_parameters - use constants - - - implicit none - - real(pr), allocatable :: tc_def(:), pc_def(:), om_def(:) - real(pr):: c1, c2, c3, c4, d1_p, d2, d3, d4, d5, e1, e2, e3, e4 - real(pr):: a, b, c0, del_1 - -contains - - subroutine get_parameteres_for_critical(eos) - implicit none - character(len=*) :: eos - - tc_def = [126.2, 304.21, 190.564, 305.32, 369.83, 408.14, 425.12, 460.43, 469.7] - pc_def = [34.0, 73.83, 45.99, 48.72, 42.48, 36.48, 37.96, 33.81, 33.7] - om_def = [0.038, 0.224, 0.012, 0.099, 0.152, 0.181, 0.20, 0.228, 0.252] - - select case(eos) - - case("SRK") - c1 = 1.6312d2 - c2 = 8.6052d1 - c3 = 4.3475d-1 - c4 = -1.8774d3 - d1_p = -1.3408d-1 - d2 = 2.5019 - d3 = 2.0846d2 - d4 = -3.9872d3 - d5 = 1.0d0 - e1 = 7.4310d-1 - e2 = 4.8122d-3 - e3 = 9.6707d-3 - e4 = -3.7184d-6 - a = -0.176 - b = 1.574 - c0 = 0.48 - !del_1 = 1.0D0 - - case("PR") ! Coeficientes originales Tabla 5.3 Pedersen - c1 = 7.34043d1 - c2 = 9.73562d1 - c3 = 6.18744d-1 - c4 = -2.05932d3 - d1_p = 7.28462d-2 - d2 = 2.18811d0 - d3 = 1.63910d2 - d4 = -4.04323d3 - d5 = 0.25d0 - e1 = 3.73765d-1 - e2 = 5.49269d-3 - e3 = 1.17934d-2 - e4 = -4.93049d-6 - a = -0.26992 - b = 1.54226 - c0 = 0.37464 - !del_1 = 1.0D0 + sqrt(2.0) - - case("RKPR") - c1 = 7.34043d1 ! Coeficientes originales Tabla 5.3 Pedersen - c2 = 9.73562d1 - c3 = 6.18744d-1 - c4 = -2.05932d3 - d1_p = 2.3d0 - d2 = -0.5d0 - d3 = 1.85d2 - d4 = -4.0d3 - d5 = 0.25d0 - e1 = 3.73765d-1 - e2 = 5.49269d-3 - e3 = 1.17934d-2 - e4 = -4.93049d-6 - a = -0.26992 - b = 1.54226 - c0 = 0.37464 - - end select - - end subroutine get_parameteres_for_critical - - end module - - - - - diff --git a/src/critical_parameters/critical_parameters.f90 b/src/critical_parameters/critical_parameters.f90 new file mode 100755 index 0000000..39bfcf4 --- /dev/null +++ b/src/critical_parameters/critical_parameters.f90 @@ -0,0 +1,72 @@ +module critical_parameters + + use constants + use dtypes + use defined_critical_parameters + +contains + + subroutine get_critical_constants(fluid, characterization, eos) + !! this subroutine doing: + !! - Calculation of Tc, Pc, omega from Pedersen correlations (adapted from code "C7plusPedersenTcPcOm") + !! - Calculation of EoS parameters (taken from CubicParam) + !use critical_parameters + + implicit none + type(FluidData), intent(in) :: fluid + type(FluidDataOut), intent(inout) :: characterization + character(len=*), intent(in) :: eos + real(pr), allocatable :: modified_c(:) + real(pr), allocatable :: tc(:) + real(pr), allocatable :: pc(:) + real(pr), allocatable :: om(:) + real(pr), allocatable :: m_funtion(:) + integer :: i + + associate(& + rho => characterization%lumped_densities, & + mw => characterization%lumped_mw & + ) + + allocate(modified_c(0)) + allocate(tc(0)) + allocate(pc(0)) + allocate(om(0)) + allocate(m_funtion(0)) + + call get_parameteres_for_critical(eos) + + tc = c1*rho + c2*log(mw) + c3*mw + (c4/mw) + pc = exp(d1_p + d2*(rho**(d5)) + (d3/mw) + (d4/(mw**2))) + + if (eos == "SRK") then + m_funtion = e1 + e2*mw + e3*rho + e4*(mw**(2)) + else + if (eos == "PR" .or. eos =="RKPR") then + m_funtion = e3*rho + 0.2 + 1.8*(1-exp(-mw/220)) + ! Alternative expression without max with MW - May 2019 + endif + endif + + modified_c = c0 - m_funtion + om = (-b + sqrt(b**2 - 4*a*modified_c))/(2*a) + + !do i = 1, size(mw) + ! print*, mw(i), rho(i), tc(i), pc(i), om(i), m_funtion(i) + !end do + + characterization%critical_temperature = tc + characterization%critical_pressure = pc + characterization%acentric_factor = om + characterization%m_funtion = m_funtion + + end associate + + end subroutine get_critical_constants + +end module critical_parameters + + + + + diff --git a/src/critical_parameters/defined_critical_param.f90 b/src/critical_parameters/defined_critical_param.f90 new file mode 100755 index 0000000..b00e1d1 --- /dev/null +++ b/src/critical_parameters/defined_critical_param.f90 @@ -0,0 +1,82 @@ +module defined_critical_parameters + use constants + + implicit none + + real(pr), allocatable :: tc_def(:), pc_def(:), om_def(:) + real(pr):: c1, c2, c3, c4, d1_p, d2, d3, d4, d5, e1, e2, e3, e4 + real(pr):: a, b, c0, del_1 + +contains + + subroutine get_parameteres_for_critical(eos) + implicit none + character(len=*) :: eos + + tc_def = [126.2, 304.21, 190.564, 305.32, 369.83, 408.14, 425.12, 460.43, 469.7] + pc_def = [34.0, 73.83, 45.99, 48.72, 42.48, 36.48, 37.96, 33.81, 33.7] + om_def = [0.038, 0.224, 0.012, 0.099, 0.152, 0.181, 0.20, 0.228, 0.252] + + select case(eos) + + case("SRK") + c1 = 1.6312d2 + c2 = 8.6052d1 + c3 = 4.3475d-1 + c4 = -1.8774d3 + d1_p = -1.3408d-1 + d2 = 2.5019 + d3 = 2.0846d2 + d4 = -3.9872d3 + d5 = 1.0d0 + e1 = 7.4310d-1 + e2 = 4.8122d-3 + e3 = 9.6707d-3 + e4 = -3.7184d-6 + a = -0.176 + b = 1.574 + c0 = 0.48 + !del_1 = 1.0D0 + + case("PR") ! Coeficientes originales Tabla 5.3 Pedersen + c1 = 7.34043d1 + c2 = 9.73562d1 + c3 = 6.18744d-1 + c4 = -2.05932d3 + d1_p = 7.28462d-2 + d2 = 2.18811d0 + d3 = 1.63910d2 + d4 = -4.04323d3 + d5 = 0.25d0 + e1 = 3.73765d-1 + e2 = 5.49269d-3 + e3 = 1.17934d-2 + e4 = -4.93049d-6 + a = -0.26992 + b = 1.54226 + c0 = 0.37464 + !del_1 = 1.0D0 + sqrt(2.0) + + case("RKPR") + c1 = 7.34043d1 ! Coeficientes originales Tabla 5.3 Pedersen + c2 = 9.73562d1 + c3 = 6.18744d-1 + c4 = -2.05932d3 + d1_p = 2.3d0 + d2 = -0.5d0 + d3 = 1.85d2 + d4 = -4.0d3 + d5 = 0.25d0 + e1 = 3.73765d-1 + e2 = 5.49269d-3 + e3 = 1.17934d-2 + e4 = -4.93049d-6 + a = -0.26992 + b = 1.54226 + c0 = 0.37464 + + end select + + end subroutine get_parameteres_for_critical + +end module defined_critical_parameters \ No newline at end of file diff --git a/src/density/density_funtion.f90 b/src/density/density_funtion.f90 new file mode 100755 index 0000000..5ce4e14 --- /dev/null +++ b/src/density/density_funtion.f90 @@ -0,0 +1,160 @@ +module density + use constants + use dtypes, only: FluidData, FluidDataOut + +contains + + subroutine density_funtion(fluid, mw_source, pho_method, characterization) + !! this subroutine ... + implicit none + type(FluidData) :: fluid + type(FluidDataOut) :: characterization + character(len=*), intent(in) :: mw_source + integer, intent(in) :: pho_method + !real(pr), allocatable :: volume_6plus_cal(:) + real(pr) :: volume_exp + real(pr) :: a_d_old, aux + real(pr) :: difference, difference_old + + + !volume_exp = sum((fluid%product_z_mw_scn)/(fluid%scn_density)) + & + ! (fluid%product_z_mw_plus/fluid%plus_density) + + + + + + + associate( & + a_d => characterization%a_d , & + b_d => characterization%b_d, & + volume_cal => characterization%volume_cal, & + volume_exp => characterization%volume_exp & + ) + + call density_method (fluid, pho_method,characterization) + + ! Now find constants for the Density function + a_d = -0.50_pr ! initial guess + b_d = 0.685_pr - a_d*exp(-0.60_pr) + call calculate_volume(fluid, mw_source, pho_method, characterization) + difference_old = volume_cal - volume_exp + a_d_old = a_d + a_d = -0.49_pr + b_d = 0.685_pr - a_d*exp(-0.60_pr) + call calculate_volume(fluid, mw_source, pho_method, characterization) + difference = volume_cal - volume_exp + + do while (abs(difference) > 0.001_pr) + aux = a_d + a_d = a_d - difference*(a_d - a_d_old)/(difference - difference_old) + b_d = 0.685_pr - a_d*exp(-0.60_pr) + a_d_old = aux + difference_old = difference + call calculate_volume(fluid, mw_source, pho_method, characterization) + difference = volume_cal - volume_exp + end do + + end associate + end subroutine density_funtion + + subroutine density_method (fluid, pho_method,characterization) + !! this subroutine ... + implicit none + type(FluidData) :: fluid + type(FluidDataOut) :: characterization + integer, intent(in) :: pho_method + real(pr) :: volume_exp + + select case (pho_method) + + case (1) + ! case 1 use experimental scn's densities reported to calculate + ! experimental volume for 6plus. The densities values are calculated + ! since CN plus + volume_exp = sum((fluid%product_z_mw_scn)/(fluid%scn_density)) + & + (fluid%product_z_mw_plus/fluid%plus_density) + characterization%volume_exp = volume_exp + + case (2) + !This case use scn densities reported to calculate experimental volume + ! this densities doesn't experimental, they can be reference values. + ! the densities values are calculated since 6 plus + + volume_exp = sum((fluid%product_z_mw_scn)/(fluid%scn_density)) + & + (fluid%product_z_mw_plus/fluid%plus_density) + characterization%volume_exp = volume_exp + + case (3) + ! This case use one experimental value (6plus) to calculate the density + ! distribution + + volume_exp = (fluid%plus6_z_exp*fluid%plus6_mw_exp)/(fluid%plus6_density_exp) + characterization%volume_exp = volume_exp + + end select + + + + + + + + + end subroutine density_method + + subroutine calculate_volume(fluid, mw_source, pho_method, characterization) + !! this subroutine ... + implicit none + type(FluidData), intent(in) :: fluid + type(FluidDataOut), intent(inout) :: characterization + character(len=*), intent(in) :: mw_source + integer, intent(in) :: pho_method + real(pr) :: volume_cal + real(pr), allocatable :: plus_density_cal(:), scn_density_cal(:), plus6_density(:) + + !allocate(density_plus_cal(0), density_scn_cal(0), plus6_density(0)) + + + associate(& + a_d => characterization%a_d, & + b_d => characterization%b_d, & + cn_plus => characterization%carbon_number_plus, & + z_m_plus_i => characterization%product_z_mw_plus_i, & + scn_z_i => characterization%scn_z, & + scn_mw_i => characterization%scn_mw, & + scn_zm => characterization%scn_zm & + ) + + scn_density_cal = ((a_d)*(exp(- (real(fluid%scn, pr))/10._pr))) + b_d + plus_density_cal =((a_d)*(exp(- (real(cn_plus, pr) )/10._pr))) + b_d + + + select case (pho_method) + + case(1) + volume_cal = sum((scn_zm)/(fluid%scn_density)) + & + sum(z_m_plus_i/plus_density_cal) + plus6_density = [fluid%scn_density, plus_density_cal] + + case(2) + volume_cal = sum((scn_zm)/(scn_density_cal)) + & + sum(z_m_plus_i/plus_density_cal) + plus6_density = [scn_density_cal, plus_density_cal] + + case(3) + volume_cal = sum((scn_zm)/(scn_density_cal)) + & + sum(z_m_plus_i/plus_density_cal) + plus6_density = [scn_density_cal, plus_density_cal] + + end select + + characterization%volume_cal = volume_cal + characterization%plus6_density = plus6_density + + end associate + + end subroutine calculate_volume + + +end module density diff --git a/src/input_data/data_from_input.mod b/src/input_data/data_from_input.mod new file mode 100755 index 0000000000000000000000000000000000000000..4031f7fdb67fa1d2ce1d564960eb7801e0de2e38 GIT binary patch literal 3919 zcmV-V53ukbiwFP!000001KnLoliJ7^GPC;A$H{Cl>y7FDxVi2e`rSF5x5oqRh-!zf zPtNM%Z0PIhxTl-dWo2&vn^~b!caE`Ss?U zJ5&mP9(>8P^D+3~@p@rC+v!k|cDn;>Hhi%LbUFTOJo$V4$8fx`#*6UfbT;V||3&a7 z72v<#Xf)~D{mJtbV09JRafP-$K;cs*_p<^CUV{LPWWgXGM_I!MF9HTCup4W3hXp?f^jTl7AREIJrITjTj~GEOdk9805}-SMj(z$EE7X)aDu0B4d9 z6F3PEr%pR>yX)EB|43J6_=l9!9gB%Sy5WzWF_ZTj+@zg#fpM$P1K%=4fwzw@(ybo9bG{TcVGZX_chTX1Yp4R>h$G!un*V$7v~hH(#D} zC_q+11kTU5RuBLm?^frFy>5?&WT4}%jVM~gpnG~~5rN_D55X}VoFA?EM%naWWRq2# zUP&c#sJcs41kVj<_QpA!x0he(ab+P_P&ULZPOPM>h`!4zavIch)U5>D`9M#kI)um! z@^S_#5P)R&$DRFpTCFeJ<01W6q5Y`>HU=(!b2d5WplL8`Z~}{ybi5f&6xI&gUfFB08?)M@p#} zk3H7T;d0up?N?xgW9oSusifUpmMUYwsz5q>`}KO@5Z41)?h~CZKa7U`P%S<}A$**d z9TKQ_b}D%JpKizAfME9P@w8uEF0e{}fF+^z_>Y^sfJnoEtEp)Z#`#Dx7XM@^y!JO;r^%4I&j)^hdG5 z#;-e0q8rs{4KP6*)ARFWupC*QJq)b=D5B%}^22PhTnxt+sQd5U=eceDGX>d)i-?`Y z73v<$tC)c9sxnm-s;X4gsH#)d@WFQ(IRHP{y=|`{!USLJ0=L<`y`7%9FS$(>nan*M z&h7pPbdHlvG?=J({$BRQ%FExE-YZWh_{kYqrf&xBlZgh?Xn%jZoEQumFOPpW5Xqo7 zf6dU_nqh2|_%;f&D=}3V_bS2AAdcENulD-5e}g?ja(_Vh`;_XgJ|zpQ1lyF>&Qmot z+$p^kJ-};zqom;>qI+s~ZWV@(cHu_%_Sv~y@f2Y+`84ch8L+jt1WY2q)R~fR0K<6T04qomSk{ODJK$Nl}FiL2?KvMNE$U=)4w03Jm^}rN? zD%w(k5g22O7|hTd+XG<=p*^+;LgW-Iv`~d?SeeeI1YpOq}(zNlLW-vE)XIL zH+pUz3Z9+^Ct5cA_4V?NIXVikG+wVld$R>*Y=*tj6Nwz?ivsYqyzZDR#8sd$onBEh$RuZ+IikjQd zdN1G_ua`9ox8mUjaGB=8&;H~)u$RT>sRiS~8iPpL{=|oFU+GnJ1tD4h_|a!QzQUcXT9Y3|;+J z19o^)d3f%j5O;W&h>vd=z1%)4igKI<8`<3xCU0(yEC{a_;QywWeaC-dCo`{6k493f z0%%ZJ;JO!V$^zGNaI14;)J_C(I@q2J%4!nrCH8Rr-GkM&k-MA^= z71U=aY+rk1h66ZX8d!%=#Kqhel@<`+Sfo+U0>7*VbwSyw7Ff`VOIJV0i^ zUW8>YDQWkOph;yf5-S2K#%77-P55jkiCi`lB`H7N? zYQNquLy*JMg=^Isp=s<8Er(lz&NvBSRM-bI#M z`AwM0Zc{45B2=91{r@N)H8puww{J_+xz z3LECaU?LV(uamRdL`qX1X`)osohTJ1NQDVfVS-eXFAH{P++@LiIqW)m zBw;ks6$Zj!QIp{dtWN{$(?B>4gwsGcEh!v9;q*YbFw6Q#g6W=MItZqNU^)nYnmbph{c-24ue;UB{^4xDe#@f&P`daArK)r!kLEuUgkXS#21sat zga$}xfP@A}Xc&ARGy)0z-h#VLS}t1?(e%VHK@1bbFhL9x#4te&6T~oe?j#}At-r01 z#+}pl_sU;Y?jKcRGfKm8uG^!7<@57r@A+~10qM$#!kKmhQS8rO(==QNM$;3I`6=v+Jih%( zPokPe7QI9;XVGT#2%&&aFb$7xOezTipZ_HxgrA*fk%UZ+%oc(o#f;83V=eLZcKYC3n|!X3dwQMef-I7m>fk!`n){Z1* zBigJ~(KLuq+z5CmOm5YCb=B1YqYBZ|_v$KnkJX}zVg37KFKbo0-M>{GX|0LXVnhsE zu$G8e)subPTu^n10;ICXy=QAQoG-YzNp80|H!g_$zKarkqZ)preJN>+8f#M83@OTy zyK1r>>DAE}n-$u?Zy`xFKCE|$tOHGKYaK)`Q~8|&s4mwgwVg^Gq_90`x90hZpv#(a z1KwU!Z@-kHI_e7e_>-^28gVN@BB8&p2gyAg*YsFt?gKzF;;p(5fXR%v?zWFfkF%QI z8msAYtmTf2sE~G=U>RzYQRWDSsLYa={*22opRgu z9srra@8WA88Qw_|;SWP-i$4r0dg%=x(L?fwA>VGnsf4Ukg#) z0Gc9ULyKP%@c@|g#I+HnK;4H&jKnn(?SnGotXqoYnV_3dKWxEGW+}Q!qydoyKHYzvi`0xQ zma@E&W174ks|?c=h1g=LLLto}lH!X+th;~Rb5mir;C zwpgAB>9RM}3smA>Z6i?j(wk(-mymhMj9XHKkI?xyqu`Gel5sJCKr?SH@Sm?qgzZ|u zt~VvwErf5YoN8r!-wX?%z$o3`nAZphY)^^ze&E}Z*E4~w6*5V79pL*4uYw fluid%def_comp_nc , & + scn_mw => characterization%scn_mw, & + plus_zm => characterization%plus_zm, & + plus_density => characterization%plus_density, & + scn_zm => characterization%scn_zm, & + plus_z => characterization%plus_z, & + scn_z => characterization%scn_z, & + plus_w => fluid%plus_w, & + w => fluid%w, & + plus_mw => characterization%plus_mw, & + carbon_number_plus => characterization%carbon_number_plus, & + plus6_density => characterization%plus6_density, & + C => characterization%C, & + scn_nc_new => characterization%scn_nc_new & + ) + + last_C = fluid%scn(fluid%scn_nc) !19 + i_last = characterization%nc_plus !124 + last = fluid%scn_nc + i_last !138 + + allocate (z_m_plus_i(0)) + allocate (plus_z_i(0)) + allocate (var_aux_1(0)) + allocate (var_aux_2(0)) + allocate (plus_z_ps(0)) + allocate (plus_mw_ps(0)) + allocate(density_ps(0)) + allocate (w_ps(0)) + allocate(moles(0)) + allocate(lumped_z(0)) + allocate(lumped_mw(0)) + allocate(lumped_densities(0)) + + + scn_nc_input = fluid%scn_nc + scn_nc_new = fluid%scn_nc_ps - 6 + ! scn_nc_ps : CN from which all SCN fractions will be lumped + ! into the specified number of pseudos + + characterization%plus_w = plus_w + z_m_plus_i = characterization%product_z_mw_plus_i + plus_z_i = characterization%plus_z_i + + if (scn_nc_new < scn_nc_input) then + ! Plus Fraction needs to be extended to include lower CN's + plus_zm = plus_zm + sum(scn_zm(scn_nc_new + 1 : scn_nc_input)) + plus_z = plus_z + sum(scn_z(scn_nc_new + 1 : scn_nc_input)) !extended zp + plus_w_new = plus_w + sum(w(fluid%def_comp_nc+scn_nc_new + 1 & + : fluid%def_comp_nc+scn_nc_input)) + !! use new variable for plus_w because type fluidata can't be modified. + plus_mw = plus_zm / plus_z + z_m_plus_i = scn_zm(scn_nc_new + 1 : scn_nc_input ) + z_m_plus_i = [z_m_plus_i, characterization%product_z_mw_plus_i] + plus_z_i = scn_z(scn_nc_new + 1 : scn_nc_input ) + plus_z_i = [plus_z_i, characterization%plus_z_i] + i_last = i_last + scn_nc_input - scn_nc_new + last_C = last_C - scn_nc_input + scn_nc_new + characterization%plus_w = plus_w_new + characterization%product_z_mw_plus_i = z_m_plus_i + characterization%plus_z_i = plus_z_i + end if + + !do i= 1, scn_nc_new + ! print*, scn_z(i), scn_mw(i) + !end do + + numbers_ps = fluid%numbers_ps + j_ps = 0 + ! Lumping into Nps pseudos + rec_zm = plus_zm / numbers_ps + ! Recommended value for the product z*M (proportional to weight) + ! for each pseudo according to pedersen + remain_plus_zm = plus_zm + i_ps = 1._pr + sum_z = 0.0_pr + sum_zm = 0.0_pr + sum_volume_ps = 0.0_pr + + i = 0.0_pr + + do while (i_ps < numbers_ps) + i = i + 1 + j_ps(i_ps) = j_ps(i_ps) + 1 + sum_z = sum_z + plus_z_i(i) + sum_zm = sum_zm + z_m_plus_i(i) + sum_volume_ps = sum_volume_ps + z_m_plus_i(i) / plus6_density(scn_nc_new+i) + if (z_m_plus_i(i+1) > 2*(rec_zm - sum_zm)) then + ! when adding one more would go too far + plus_z_ps = [plus_z_ps, sum_z] + plus_mw_ps = [plus_mw_ps, sum_zm/sum_z] + w_ps = [w_ps, (characterization%plus_w*(sum_zm/plus_zm))] + remain_plus_zm = remain_plus_zm - sum_zm + if (remain_plus_zm < rec_zm) numbers_ps = i_ps + 1 + carbon_number_plus(scn_nc_new + i_ps) = 6+((plus_mw_ps(i_ps)-84)/C) + density_ps = [density_ps, (sum_zm/sum_volume_ps)] + i_ps = i_ps + 1 + sum_z = 0._pr + sum_zm = 0._pr + sum_volume_ps = 0._pr + end if + end do + + plus_z_ps = [plus_z_ps, sum(plus_z_i(i+1:i_last))] + ! at this point, Nps is the order for the last NON-Asphaltenic pseudo comp. + ! (e.g. 4 if 5 is for Asphaltenes) + plus_mw_ps =[plus_mw_ps, sum(z_m_plus_i(i+1:i_last))/plus_z_ps(numbers_ps)] + w_ps = [w_ps, characterization%plus_w * ((plus_z_ps(numbers_ps))* & + (plus_mw_ps(numbers_ps))/(plus_zm))] + carbon_number_plus(scn_nc_new + numbers_ps) = 6 + & + ((plus_mw_ps(numbers_ps) - 84) / C) + sum_zm_last_ps = (sum(z_m_plus_i(i+1:i_last))) + density_ps = [density_ps, (sum_zm_last_ps)/sum((z_m_plus_i(i+1:i_last)) / (plus6_density(scn_nc_new+i+1:last)))] + j_ps(numbers_ps) = i_last - sum(j_ps(1: numbers_ps-1)) !! revisar manana + + !do i = 1, numbers_ps + ! print*, 'ps',i, plus_z_ps(i), plus_mw_ps(i), density_ps(i), j_ps(i), carbon_number_plus(scn_nc_new+i) + !end do + + plus_density = plus_zm / sum(plus_z_ps(1:numbers_ps)* & + plus_mw_ps(1:numbers_ps)/(density_ps(1:numbers_ps))) + + moles = [(fluid%def_comp_w)/(fluid%def_comp_mw)] + moles = [moles, fluid%w(def_nc+1:def_nc+scn_nc_new)/ scn_mw(1:scn_nc_new)] + moles = [moles, w_ps/plus_mw_ps ] + + characterization%mol_fraction = moles/sum(moles) ! mole fractions normalize + characterization%last_C = last_C + characterization%i_last = i_last + characterization%last = last + + ! In this point we create a new array to content final Mis and densities + ! whit pseudos. + + lumped_z = [scn_z(1:scn_nc_new), plus_z_ps] + lumped_mw = [scn_mw(1:scn_nc_new), plus_mw_ps] + lumped_densities = [plus6_density(1:scn_nc_new), density_ps ] + + + characterization%lumped_z = lumped_z + characterization%lumped_mw = lumped_mw + characterization%lumped_densities = lumped_densities + + end associate + + end subroutine lump + +end module lumping + diff --git a/src/molecular_weigth/mw_calculations.f90 b/src/molecular_weigth/mw_calculations.f90 new file mode 100755 index 0000000..de34db4 --- /dev/null +++ b/src/molecular_weigth/mw_calculations.f90 @@ -0,0 +1,268 @@ +module molecular_weigth + use constants + use dtypes, only: FluidData, FluidDataOut + use bfr_routines + +contains + + subroutine select_method(fluid, mw_source, method, characterization) + !! This subroutine defines the calculation method to perform the characterization,& + !! based on the available experimental data. If the molecular weight data is experimental, + !! define the mole fractions, molecular weights and densities directly from the input file. + !! On the other hand, if the molecular weights are assumed, they are recalculated according + !! to the methodology described by Martin et al and the molar fractions of the fluid are recalculated. + + implicit none + type(FluidData), intent(inout) :: fluid !! Derived type to save inlet or experimental information of fluid to be characterized. + type(FluidDataOut), intent(inout) :: characterization !! Derive type variable to save results of characterization methodology. + character(len=*), intent(in), optional :: method + !! this option allows choose beetwen obtaining plus_mw or global_mw reported in the experimental information . + character(len=*), intent(in) :: mw_source !! This variable indicates the source of the input data, that is, whether they are experimental or not. + real(pr) :: sum_def_comp_z_plus !! compositions sum from (define component +1) to (plus component) + real(pr), allocatable :: def_comp_moles(:) !! + real(pr), allocatable :: scn_moles(:) !! + real(pr) :: plus_moles + real(pr) :: total_moles + + associate (& + scn_nc => fluid%scn_nc, def_comp_nc => fluid%def_comp_nc, & + def_comp_w => fluid%def_comp_w, def_comp_mw => fluid%def_comp_mw, & + scn_w => fluid%scn_w, plus_w => fluid%plus_w, scn => fluid%scn, & + plus_mw => characterization%plus_mw, scn_z => characterization%scn_z, & + plus_z => characterization%plus_z, log_scn_z => characterization%log_scn_z, & + C => characterization%C, scn_mw => characterization%scn_mw, & + scn_zm => characterization%scn_zm, & + plus_zm => characterization%plus_zm & + ) + + allocate(def_comp_moles(def_comp_nc)) + allocate(scn_moles(scn_nc)) + + def_comp_moles = (def_comp_w)/(def_comp_mw) + + select case (mw_source) + + case("experimental") + scn_mw = fluid%scn_mw + scn_moles = (scn_w)/(scn_mw) + plus_moles = (plus_w)/(plus_mw) + total_moles = sum(def_comp_moles)+sum(scn_moles)+(plus_moles) + scn_z = scn_moles / total_moles + plus_z = plus_moles / total_moles + + case("calculated") + + if (method=="global_mw")then + scn_mw = 84 + C*(scn-6) + scn_z = fluid%product_z_mw_scn / scn_mw + sum_def_comp_z_plus = sum(fluid%scn_z) + (fluid%plus_z) + plus_z = sum_def_comp_z_plus - sum(scn_z) ! new molar fraction for residual cut, based on C + plus_mw = fluid%product_z_mw_plus / plus_z + end if + + if (method=="plus_mw")then + scn_mw = 84 + C*(scn-6) + scn_moles = (scn_w)/(scn_mw) + plus_moles = (plus_w)/(plus_mw) + total_moles = sum(def_comp_moles)+sum(scn_moles)+(plus_moles) + scn_z = scn_moles / total_moles + plus_z = plus_moles / total_moles + + endif + end select + scn_zm = scn_z*scn_mw + plus_zm = plus_z*plus_mw + log_scn_z = log(scn_z) + end associate + + end subroutine select_method + + subroutine difference_mw_plus(fluid, mw_source, method, start, difference, & + characterization) + + !! this subroutine ... + implicit none + + type(FluidData), intent(inout) :: fluid + type(FluidDataOut), intent(inout) :: characterization + logical :: start + character(len=*), intent(in), optional :: method + character(len=*), intent(in) :: mw_source + real(pr) :: plus_mw_cal !! calculated molecular weight of residual fraction + real(pr), intent(out) :: difference + real(pr) :: half + real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux + real(pr):: sum_z + real(pr) :: denom + !! set to molecular weights calculated by \[M = 84-characterization%C(i-6)\] + real(pr) :: sum_def_comp_z_plus + !! compositions sum from (define component +1) to (plus component) + integer :: j, k , k_old, n_best, x_aux, i_0! internal variables + integer :: i + integer, dimension(300) :: carbon_number_plus + real(pr), dimension(300) :: plus_z_i + real(pr), dimension(300) :: product_z_mw_plus_i + real(pr), dimension(300) :: scn_i + + associate (& + scn_nc => fluid%scn_nc, scn => fluid%scn,& + plus_z => characterization%plus_z, plus_mw => characterization%plus_mw, & + log_scn_z => characterization%log_scn_z, C => characterization%C, & + a_blr => characterization%a_blr, b_blr => characterization%b_blr, & + a_60 => characterization%a_60, b_60 => characterization%b_60, & + a_lim => characterization%a_lim, b_lim => characterization%b_lim, & + r2 => characterization%r2, n_init => characterization%n_init, & + c_max_blr => characterization%c_max_lim, a => characterization%a, & + c_max_lim => characterization%c_max_lim, b => characterization%b, & + c_max_60 => characterization%c_max_lim, & + c_max => characterization%c_max & + ) + + 1 i_0 = scn(scn_nc) + + call select_method(fluid, mw_source, method, characterization) ! select case + call Best_Linear_Regression(scn_nc, scn, log_scn_z, & + plus_z, a_blr, b_blr, r2, n_init, c_max_blr) + call LimitLine(scn_nc, scn, plus_z, a_blr, b_blr,a_lim,b_lim,c_max_lim,half) + call Line_C60_max (scn_nc,scn,plus_z,a_blr,b_blr,half,a_60,b_60,c_max_60) ! add line by oscar. + + if(a_blr < a_lim)then + !Added by oscar 05/12/2023.se elimina por que se agregan las restricciones + !para characterization%C>12 y characterization%C<12. + a = a_lim + b = b_lim + else ! this line was removed call Line_C60_max + if(a_blr > a_60)then + a = a_60 + b = b_60 + else + a = a_blr + b = b_blr + end if + end if + + sum_z = 0.0_pr + i = 0 + do while (sum_z < plus_z .and. i< 300) + i = i+1 + plus_z_i(i) = exp(a*(i+i_0)+b) + sum_z = sum_z + plus_z_i(i) + product_z_mw_plus_i(i) = plus_z_i(i)*(84+C*(i+i_0-6)) + carbon_number_plus(i) = i+i_0 + end do + + if (i == 300.and.sum_z < plus_z)then + if(start) C = C - 0.07 + if(.not.start) C = C - 0.01 + go to 1 + end if + + c_max = i+i_0 + plus_z_i(i) = plus_z_i(i) - (sum_z-plus_z) ! Adjustment to Zp (z20+) + product_z_mw_plus_i(i) = plus_z_i(i)*(84+C*(i+i_0-6)) + + if (mw_source == "experimental" .and. C>12 .and. C<14 ) then + ! mayor igual a 12 o menor igual 14. + denom = sum((plus_z_i(1:i))*((carbon_number_plus(1:i)-6))) + ! denominator of equation to obtain C value directly + C = (plus_z*(plus_mw-84))/denom + endif + + plus_mw_cal = sum(product_z_mw_plus_i(1:i))/plus_z + difference = plus_mw_cal-plus_mw + scn_i = [scn, carbon_number_plus(1:i)] + end associate + + characterization%plus_z_i = plus_z_i(1:i) + characterization%product_z_mw_plus_i = product_z_mw_plus_i(1:i) + characterization%carbon_number_plus = carbon_number_plus(1:i) + characterization%nc_plus = i + characterization%scn_i = scn_i + + + + end subroutine difference_mw_plus + + subroutine get_c_or_m_plus(fluid, mw_source, method, fix_C, characterization) + !! This subroutine... + implicit none + type(FluidData) :: fluid + type(FluidDataOut) :: characterization + logical :: start + character(len=*), intent(in) :: mw_source + character(len=*), intent(in), optional :: method + logical, intent(in) :: fix_C !! If C is < 12 or > 14 then fix to the + !! closest value. + !integer :: c_max !! output CN at which characterization%plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) + real(pr), allocatable :: log_scn_z(:) !! logarithm of corresponding mole fractions of scn cutsn + real(pr) :: difference !! + real(pr) :: difference_old, plus_mw_old + real(pr) :: C_old + real(pr) :: aux + + associate (C => characterization%C, plus_mw => characterization%plus_mw) + + select case (mw_source) + case("experimental") + C = 13.5 + Start = .true. + plus_mw = fluid%plus_mw + call difference_mw_plus(fluid, mw_source, method, start, difference, & + characterization) + + case("calculated") + + if(method == 'global_mw') C=13 + if(method == 'plus_mw') C=14 + + Start = .true. + plus_mw = fluid%plus_mw + call difference_mw_plus(fluid, mw_source, method, start, & + difference_old, characterization) + C_old = C + if(method == 'global_mw') C = min(12.8, C-0.07) + if(method == 'plus_mw') C = 13.5 + Start = .false. + call difference_mw_plus(fluid, mw_source, method, start, difference, & + characterization) + + do while (abs(difference) > 0.1) + aux = C + C = C - difference*(C-C_old)/(difference-difference_old) + C_old = aux + difference_old = difference + call difference_mw_plus(fluid, mw_source, method, start, & + difference, characterization) + end do + end select + + if(fix_C .and. C>14 .or. C<12)then + + + if(C>14) C=14 + if(C<12) C=12 + + Start = .true. + plus_mw = fluid%plus_mw + call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & + difference_old, characterization) + plus_mw_old = plus_mw + plus_mw = 0.9*plus_mw + Start = .false. + call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & + difference, characterization) + + do while (abs(difference)>0.00001) + aux = plus_mw + plus_mw = plus_mw - difference*(plus_mw-plus_mw_old)/(difference-difference_old) + plus_mw_old = aux + difference_old = difference + call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & + difference, characterization) + end do + end if + end associate + end subroutine get_c_or_m_plus + +end module molecular_weigth + diff --git a/src/routines.f90 b/src/routines.f90 deleted file mode 100644 index 107c7d8..0000000 --- a/src/routines.f90 +++ /dev/null @@ -1,852 +0,0 @@ -module routines - use constants - use dtypes, only: FluidData, FluidDataOut - use data_from_input, only: data_from_file, FluidData - -contains - - subroutine Linear_Regression(x, y, a, b, r2) - !! This subroutine computes the regression line for a data set of x, y variables. - implicit none - - integer :: i - integer :: n !! total number of data set. - real(pr), intent(in) :: x(:) !! x: input array of length n which contains the set of independent variable - real(pr), intent(in) :: y(:) !! y: input array of length n which contains the set of dependent variable - real(pr) :: t1,t2,t3,t4, aux1,aux2,aux3,aux4,aux5 ! internal variables - real(pr), intent(out) :: a !! a: output real variable. Slope of the regression line - real(pr), intent(out) :: b !! b: output real variable. Intercept of the regression line - real(pr), intent(out) :: r2 !! r2: output real variable. Square correlation coefficient - - n = size(x) - !Calculation of a y b --> y=a*x+b - t1=0.;t2=0.;t3=0.;t4=0.; - - do i=1,n - t1=t1+x(i)*y(i) - t2=t2+x(i) - t3=t3+y(i) - t4=t4+x(i)**2 - end do - - a=(n*t1-t2*t3)/(n*t4-t2**2) - b=(t3-a*t2)/n - !coefficient calculation of correlations r2 - aux1=0.;aux2=0.;aux3=0.;aux4=0.;aux5=0.; - - do i=1, n - aux1= aux1 + x(i)*y(i) - aux2= aux2 + x(i) - aux3= aux3 + y(i) - aux4= aux4 + x(i)**2 - aux5= aux5 + y(i)**2 - end do - - r2=(aux1-aux2*aux3/n)**2 /((aux4-aux2**2/n)*(aux5-aux3**2/n)) - - end subroutine Linear_Regression - - subroutine Best_Linear_Regression(scn_nc, scn, scn_z, plus_z, a, b, r2, & - n_init, c_max_blr) - !! This subroutine calculates the best regression line for an fluid. - implicit none - - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid - integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid - integer, intent(out) :: c_max_blr !! output CN at which plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) - real(pr), intent(in) :: scn_z(:) !! set of corresponding mole fractions of scn cuts - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file - real(pr), intent(out) :: a !! output real variable. Slope of the best regression line. - real(pr), intent(out) :: b !! output real variable. Intercept of the best regression line. - real(pr), intent(out) :: r2 !! output real variable. Square correlation coefficient. - integer, intent(out) :: n_init !! minimum carbon number obtained from the best linear regression - integer :: i, j, k , k_old, n_best, x_aux ! internal variables - real(pr), dimension(scn_nc) :: x_blr, y_blr ! x_blr: carbon number vector, y_blr: logarithm of mole fraction vector; for each step - real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux ! auxiliary variables - - k=5 - r2=0.0001_pr - r2_old=0.00001_pr - r2_best=0.0001_pr - - do while (r2.gt.r2_old.or.r2_old.lt.0.9) - k_old=k - r2_old=r2 - a_old=a - b_old=b - - if (r2.gt.r2_best) then - r2_best=r2 - a_best=a - b_best=b - n_best=scn(scn_nc-k+2) - end if - - if (k.gt.scn_nc) then - if (r2.gt.r2_best)then - n_init=scn(scn_nc-k+2) - go to 22 - else - r2=r2_best - a=a_best - b=b_best - n_init=n_best - go to 22 - end if - end if - - j=1 - x_blr=0. - y_blr=0. - - do i=scn_nc-k+1, scn_nc - x_blr(j)=scn(i) - y_blr(j)=scn_z(i) - j=j+1 - end do - - call Linear_Regression(x_blr(:k), y_blr(:k), a, b, r2) - k=k+1 - end do - - r2=r2_old - a=a_old - b=b_old - n_init=scn(scn_nc-k_old+2) - - 22 continue - - z_sum = 0d0 - x_aux=scn(scn_nc) - - do while (z_sum.lt.plus_z.and.x_aux<300) - x_aux=x_aux+1 - z_aux= exp(a*x_aux+b) - z_sum=z_sum+z_aux - end do - - c_max_blr=x_aux - - !print*, a, b, c_max_blr, r2, n_init - - end subroutine Best_Linear_Regression - - subroutine LimitLine(scn_nc, scn, plus_z, a_blr, b_blr, a_lim, b_lim, & - c_max_lim, half) - !! This subroutine obtains the limit line constants for an fluid. - implicit none - - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid. - integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid. - integer, intent(out) :: c_max_lim !! maximum carbon number obtained for limit line. - !! output CN at which plus_z is reached, as the summation of single z(i) from the limit distribution (blr). - real(pr), intent(in) :: a_blr !! input constant from the Best linear regression. - real(pr), intent(in) :: b_blr !! input constant from the Best linear regression. - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file. - real(pr), intent(out) :: a_lim !! output real variable. Slope of the limit line. - real(pr), intent(out) :: b_lim !! output real variable. Intercept of the limit line. - real(pr), intent(out) :: half ! variable defined in the numerical method to converge the area under the curve of the function. - real(pr) :: z_lim, cross_cn, z_cross, z_aux - integer :: x_aux ! internal variable - - !A and B limit calculation. - z_lim = 0.0_pr - half = 0.506_pr ! modified by oscar, the variable half was removed from the argument - - do while (z_lim.lt.plus_z) - half = half + 0.002d0 - cross_cn = scn(scn_nc)+half ! typically 19.508 in first try - z_cross = exp((a_blr*cross_cn) + b_blr) - a_lim = -(z_cross)/plus_z - b_lim = log(z_cross)-(a_lim*cross_cn) - x_aux=scn(scn_nc) - z_lim = 0._pr - do while (z_lim.lt.plus_z.and.x_aux<300) - x_aux=x_aux+1 - z_aux= exp(a_lim*x_aux+b_lim) - z_lim=z_lim+z_aux - end do - continue - end do - - c_max_lim = x_aux - - end subroutine LimitLine - - subroutine Line_C60_max (scn_nc, scn, plus_z, a_blr, b_blr, half, a_60, & - b_60, c_max_60) - !! This subroutine obtains the Cmax60 line constants for an fluid. - - implicit none - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid - integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid - integer, intent(out) :: c_max_60 !!output CN at which Zp is reached, as the summation of single z(i) from the Cmax60 distribution. - real(pr), intent(in) :: a_blr !! input constant from the Best linear regression - real(pr), intent(in) :: b_blr !! input constant from the Best linear regression - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file - real(pr), intent(out) :: a_60 !! output real variable. Slope of the limit line - real(pr), intent(out) :: b_60 !! output real variable. Intercept of the C60max line. - real(pr), intent(in) :: half !variable defined in the numerical method to converge the area under the curve of the function - integer :: x_aux ! auxiliary variable - real(pr) :: z_sum, cross_cn, z_cross, z_aux, F_tol, a_tol, var_range ! internal variables - real(pr) :: a_old, Zp_60, F, dF_dA ! internal variables - - cross_cn = scn(scn_nc)+half ! typically 19.51 - z_cross= exp((a_blr*cross_cn)+b_blr) - F_tol=1_pr - a_tol=1_pr - a_60=a_blr - var_range = 60.5_pr - cross_cn - - do while (a_tol.gt.1e-6_pr.or.F_tol.gt.1e-6_pr) - a_old=a_60 - Zp_60 = (exp(a_60*var_range+log(z_cross))-z_cross)/a_60 - F = plus_z - Zp_60 - dF_dA = - (var_range*(a_60*Zp_60+z_cross)-Zp_60) / a_60 ! modified by oscar, error in the derivative in previous version is corrected - a_60=a_old-(F/dF_dA) - a_tol=abs(a_old-a_60) - F_tol=abs(F) - end do - - b_60=log(z_cross)-a_60*cross_cn - z_sum=0.0_pr - x_aux=scn(scn_nc) - - do while (z_sum.lt.plus_z) - x_aux=x_aux+1 - z_aux= exp(a_60*x_aux+b_60) - z_sum=z_sum+z_aux - end do - - c_max_60=x_aux - - !print*, a_60,b_60,c_max_60 - end subroutine Line_C60_max - - subroutine select_method(fluid, mw_source, method, characterization) - !! This subroutine defines the calculation method to perform the characterization,& - !! based on the available experimental data. If the molecular weight data is experimental, - !! define the mole fractions, molecular weights and densities directly from the input file. - !! On the other hand, if the molecular weights are assumed, they are recalculated according - !! to the methodology described by Martin et al and the molar fractions of the fluid are recalculated. - - implicit none - type(FluidData), intent(inout) :: fluid !! Derived type to save inlet or experimental information of fluid to be characterized. - type(FluidDataOut), intent(inout) :: characterization !! Derive type variable to save results of characterization methodology. - character(len=*), intent(in), optional :: method - !! this option allows choose beetwen obtaining plus_mw or global_mw reported in the experimental information . - character(len=*), intent(in) :: mw_source !! This variable indicates the source of the input data, that is, whether they are experimental or not. - real(pr) :: sum_def_comp_z_plus !! compositions sum from (define component +1) to (plus component) - real(pr), allocatable :: def_comp_moles(:) !! - real(pr), allocatable :: scn_moles(:) !! - real(pr) :: plus_moles - real(pr) :: total_moles - - associate (& - scn_nc => fluid%scn_nc, def_comp_nc => fluid%def_comp_nc, & - def_comp_w => fluid%def_comp_w, def_comp_mw => fluid%def_comp_mw, & - scn_w => fluid%scn_w, plus_w => fluid%plus_w, scn => fluid%scn, & - plus_mw => characterization%plus_mw, scn_z => characterization%scn_z, & - plus_z => characterization%plus_z, log_scn_z => characterization%log_scn_z, & - C => characterization%C, scn_mw => characterization%scn_mw, & - scn_zm => characterization%scn_zm, & - plus_zm => characterization%plus_zm & - ) - - allocate(def_comp_moles(def_comp_nc)) - allocate(scn_moles(scn_nc)) - - def_comp_moles = (def_comp_w)/(def_comp_mw) - - select case (mw_source) - - case("experimental") - scn_mw = fluid%scn_mw - scn_moles = (scn_w)/(scn_mw) - plus_moles = (plus_w)/(plus_mw) - total_moles = sum(def_comp_moles)+sum(scn_moles)+(plus_moles) - scn_z = scn_moles / total_moles - plus_z = plus_moles / total_moles - - case("calculated") - - if (method=="global_mw")then - scn_mw = 84 + C*(scn-6) - scn_z = fluid%product_z_mw_scn / scn_mw - sum_def_comp_z_plus = sum(fluid%scn_z) + (fluid%plus_z) - plus_z = sum_def_comp_z_plus - sum(scn_z) ! new molar fraction for residual cut, based on C - plus_mw = fluid%product_z_mw_plus / plus_z - end if - - if (method=="plus_mw")then - scn_mw = 84 + C*(scn-6) - scn_moles = (scn_w)/(scn_mw) - plus_moles = (plus_w)/(plus_mw) - total_moles = sum(def_comp_moles)+sum(scn_moles)+(plus_moles) - scn_z = scn_moles / total_moles - plus_z = plus_moles / total_moles - - endif - end select - scn_zm = scn_z*scn_mw - plus_zm = plus_z*plus_mw - log_scn_z = log(scn_z) - end associate - - end subroutine select_method - - subroutine difference_mw_plus(fluid, mw_source, method, start, difference, & - characterization) - - !! this subroutine ... - implicit none - - type(FluidData), intent(inout) :: fluid - type(FluidDataOut), intent(inout) :: characterization - logical :: start - character(len=*), intent(in), optional :: method - character(len=*), intent(in) :: mw_source - real(pr) :: plus_mw_cal !! calculated molecular weight of residual fraction - real(pr), intent(out) :: difference - real(pr) :: half - real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux - real(pr):: sum_z - real(pr) :: denom - !! set to molecular weights calculated by \[M = 84-characterization%C(i-6)\] - real(pr) :: sum_def_comp_z_plus - !! compositions sum from (define component +1) to (plus component) - integer :: j, k , k_old, n_best, x_aux, i_0! internal variables - integer :: i - integer, dimension(300) :: carbon_number_plus - real(pr), dimension(300) :: plus_z_i - real(pr), dimension(300) :: product_z_mw_plus_i - real(pr), dimension(300) :: scn_i - - associate (& - scn_nc => fluid%scn_nc, scn => fluid%scn,& - plus_z => characterization%plus_z, plus_mw => characterization%plus_mw, & - log_scn_z => characterization%log_scn_z, C => characterization%C, & - a_blr => characterization%a_blr, b_blr => characterization%b_blr, & - a_60 => characterization%a_60, b_60 => characterization%b_60, & - a_lim => characterization%a_lim, b_lim => characterization%b_lim, & - r2 => characterization%r2, n_init => characterization%n_init, & - c_max_blr => characterization%c_max_lim, a => characterization%a, & - c_max_lim => characterization%c_max_lim, b => characterization%b, & - c_max_60 => characterization%c_max_lim, & - c_max => characterization%c_max & - ) - - 1 i_0 = scn(scn_nc) - - call select_method(fluid, mw_source, method, characterization) ! select case - call Best_Linear_Regression(scn_nc, scn, log_scn_z, & - plus_z, a_blr, b_blr, r2, n_init, c_max_blr) - call LimitLine(scn_nc, scn, plus_z, a_blr, b_blr,a_lim,b_lim,c_max_lim,half) - call Line_C60_max (scn_nc,scn,plus_z,a_blr,b_blr,half,a_60,b_60,c_max_60) ! add line by oscar. - - if(a_blr < a_lim)then - !Added by oscar 05/12/2023.se elimina por que se agregan las restricciones - !para characterization%C>12 y characterization%C<12. - a = a_lim - b = b_lim - else ! this line was removed call Line_C60_max - if(a_blr > a_60)then - a = a_60 - b = b_60 - else - a = a_blr - b = b_blr - end if - end if - - sum_z = 0.0_pr - i = 0 - do while (sum_z < plus_z .and. i< 300) - i = i+1 - plus_z_i(i) = exp(a*(i+i_0)+b) - sum_z = sum_z + plus_z_i(i) - product_z_mw_plus_i(i) = plus_z_i(i)*(84+C*(i+i_0-6)) - carbon_number_plus(i) = i+i_0 - end do - - if (i == 300.and.sum_z < plus_z)then - if(start) C = C - 0.07 - if(.not.start) C = C - 0.01 - go to 1 - end if - - c_max = i+i_0 - plus_z_i(i) = plus_z_i(i) - (sum_z-plus_z) ! Adjustment to Zp (z20+) - product_z_mw_plus_i(i) = plus_z_i(i)*(84+C*(i+i_0-6)) - - if (mw_source == "experimental" .and. C>12 .and. C<14 ) then - ! mayor igual a 12 o menor igual 14. - denom = sum((plus_z_i(1:i))*((carbon_number_plus(1:i)-6))) - ! denominator of equation to obtain C value directly - C = (plus_z*(plus_mw-84))/denom - endif - - plus_mw_cal = sum(product_z_mw_plus_i(1:i))/plus_z - difference = plus_mw_cal-plus_mw - scn_i = [scn, carbon_number_plus(1:i)] - end associate - - characterization%plus_z_i = plus_z_i(1:i) - characterization%product_z_mw_plus_i = product_z_mw_plus_i(1:i) - characterization%carbon_number_plus = carbon_number_plus(1:i) - characterization%nc_plus = i - characterization%scn_i = scn_i - - - - end subroutine difference_mw_plus - - subroutine get_c_or_m_plus(fluid, mw_source, method, fix_C, characterization) - !! This subroutine... - implicit none - type(FluidData) :: fluid - type(FluidDataOut) :: characterization - logical :: start - character(len=*), intent(in) :: mw_source - character(len=*), intent(in), optional :: method - logical, intent(in) :: fix_C !! If C is < 12 or > 14 then fix to the - !! closest value. - !integer :: c_max !! output CN at which characterization%plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) - real(pr), allocatable :: log_scn_z(:) !! logarithm of corresponding mole fractions of scn cutsn - real(pr) :: difference !! - real(pr) :: difference_old, plus_mw_old - real(pr) :: C_old - real(pr) :: aux - - associate (C => characterization%C, plus_mw => characterization%plus_mw) - - select case (mw_source) - case("experimental") - C = 13.5 - Start = .true. - plus_mw = fluid%plus_mw - call difference_mw_plus(fluid, mw_source, method, start, difference, & - characterization) - - case("calculated") - - if(method == 'global_mw') C=13 - if(method == 'plus_mw') C=14 - - Start = .true. - plus_mw = fluid%plus_mw - call difference_mw_plus(fluid, mw_source, method, start, & - difference_old, characterization) - C_old = C - if(method == 'global_mw') C = min(12.8, C-0.07) - if(method == 'plus_mw') C = 13.5 - Start = .false. - call difference_mw_plus(fluid, mw_source, method, start, difference, & - characterization) - - do while (abs(difference) > 0.1) - aux = C - C = C - difference*(C-C_old)/(difference-difference_old) - C_old = aux - difference_old = difference - call difference_mw_plus(fluid, mw_source, method, start, & - difference, characterization) - end do - end select - - if(fix_C .and. C>14 .or. C<12)then - - - if(C>14) C=14 - if(C<12) C=12 - - Start = .true. - plus_mw = fluid%plus_mw - call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & - difference_old, characterization) - plus_mw_old = plus_mw - plus_mw = 0.9*plus_mw - Start = .false. - call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & - difference, characterization) - - do while (abs(difference)>0.00001) - aux = plus_mw - plus_mw = plus_mw - difference*(plus_mw-plus_mw_old)/(difference-difference_old) - plus_mw_old = aux - difference_old = difference - call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & - difference, characterization) - end do - end if - end associate - end subroutine get_c_or_m_plus - - subroutine density_funtion(fluid, mw_source, characterization) - !! this subroutine ... - implicit none - type(FluidData) :: fluid - type(FluidDataOut) :: characterization - character(len=*), intent(in) :: mw_source - !real(pr), allocatable :: volume_6plus_cal(:) - real(pr) :: volume_6plus_exp - real(pr) :: a_d_old, aux - real(pr) :: difference, difference_old - - - volume_6plus_exp = sum((fluid%product_z_mw_scn)/(fluid%scn_density)) + & - (fluid%product_z_mw_plus/fluid%plus_density) - - associate( & - a_d => characterization%a_d , & - b_d => characterization%b_d, & - volume_6plus_cal => characterization%volume_6plus_cal & - ) - ! Now find constants for the Density function - a_d = -0.50_pr ! initial guess - b_d = 0.685_pr - a_d*exp(-0.60_pr) - call calculate_volume_6plus(fluid, mw_source, characterization) - difference_old = volume_6plus_cal - volume_6plus_exp - a_d_old = a_d - a_d = -0.49_pr - b_d = 0.685_pr - a_d*exp(-0.60_pr) - call calculate_volume_6plus(fluid, mw_source, characterization) - difference = volume_6plus_cal - volume_6plus_exp - - do while (abs(difference) > 0.001_pr) - aux = a_d - a_d = a_d - difference*(a_d - a_d_old)/(difference - difference_old) - b_d = 0.685_pr - a_d*exp(-0.60_pr) - a_d_old = aux - difference_old = difference - call calculate_volume_6plus(fluid, mw_source, characterization) - difference = volume_6plus_cal - volume_6plus_exp - end do - - end associate - end subroutine density_funtion - - subroutine calculate_volume_6plus(fluid, mw_source, characterization) - !! this subroutine ... - implicit none - type(FluidData), intent(in) :: fluid - type(FluidDataOut), intent(inout) :: characterization - character(len=*), intent(in) :: mw_source - real(pr) :: volume_6plus_cal - real(pr), allocatable :: plus_density_cal(:), scn_density_cal(:), plus6_density(:) - - !allocate(density_plus_cal(0), density_scn_cal(0), plus6_density(0)) - - - associate(& - a_d => characterization%a_d, & - b_d => characterization%b_d, & - cn_plus => characterization%carbon_number_plus, & - z_m_plus_i => characterization%product_z_mw_plus_i, & - scn_z_i => characterization%scn_z, & - scn_mw_i => characterization%scn_mw, & - scn_zm => characterization%scn_zm & - ) - - scn_density_cal = ((a_d)*(exp(- (real(fluid%scn, pr))/10._pr))) + b_d - plus_density_cal =((a_d)*(exp(- (real(cn_plus, pr) )/10._pr))) + b_d - - - select case (mw_source) - case("experimental") - volume_6plus_cal = sum((scn_zm)/(fluid%scn_density)) + & - sum(z_m_plus_i/plus_density_cal) - plus6_density = [fluid%scn_density, plus_density_cal] - - case("calculated") - volume_6plus_cal = sum((scn_zm)/(scn_density_cal)) + & - sum(z_m_plus_i/plus_density_cal) - plus6_density = [scn_density_cal, plus_density_cal] - end select - - characterization%volume_6plus_cal = volume_6plus_cal - characterization%plus6_density = plus6_density - - end associate - - end subroutine calculate_volume_6plus - - subroutine lump(fluid, characterization) - !! This sobroutine ... - implicit none - type(FluidData), intent(in) :: fluid - type(FluidDataOut), intent(inout) :: characterization - integer :: last_C ! CN of last single cut: typically 19 - integer :: i_last ! Number of elements in the distribution of CN+ - integer :: last ! Total elements starting after defined components. - integer :: scn_nc_input !! inicial number of single cuts being considered in the oil from data input - integer :: scn_nc_new !! new number of single cuts being considered in the oil defined as from scn_nc_ps variable. - real(pr) :: plus_w_new - real(pr), allocatable :: z_m_plus_i(:) - real(pr), allocatable :: plus_z_i(:) - real(pr), allocatable :: var_aux_1(:) - real(pr), allocatable :: var_aux_2(:) - integer :: i, prev_i - integer, dimension(15) :: j_ps - integer :: i_ps ! auxiliary variables - real(pr) :: rec_zm, remain_plus_zm, sum_z, sum_zm, sum_volume_ps ! auxiliary variables - real(pr), allocatable :: plus_z_ps(:) - real(pr), allocatable :: plus_mw_ps(:) - integer :: numbers_ps - real(pr), allocatable :: density_ps(:) - real(pr), allocatable :: w_ps(:) - real(pr) :: sum_zm_last_ps - real(pr), allocatable :: moles(:) - real(pr), allocatable :: lumped_z(:) - real(pr), allocatable :: lumped_mw(:) - real(pr), allocatable :: lumped_densities(:) - - associate (& - def_nc => fluid%def_comp_nc , & - scn_mw => characterization%scn_mw, & - plus_zm => characterization%plus_zm, & - plus_density => characterization%plus_density, & - scn_zm => characterization%scn_zm, & - plus_z => characterization%plus_z, & - scn_z => characterization%scn_z, & - plus_w => fluid%plus_w, & - w => fluid%w, & - plus_mw => characterization%plus_mw, & - carbon_number_plus => characterization%carbon_number_plus, & - plus6_density => characterization%plus6_density, & - C => characterization%C, & - scn_nc_new => characterization%scn_nc_new & - ) - - last_C = fluid%scn(fluid%scn_nc) !19 - i_last = characterization%nc_plus !124 - last = fluid%scn_nc + i_last !138 - - allocate (z_m_plus_i(0)) - allocate (plus_z_i(0)) - allocate (var_aux_1(0)) - allocate (var_aux_2(0)) - allocate (plus_z_ps(0)) - allocate (plus_mw_ps(0)) - allocate(density_ps(0)) - allocate (w_ps(0)) - allocate(moles(0)) - allocate(lumped_z(0)) - allocate(lumped_mw(0)) - allocate(lumped_densities(0)) - - - scn_nc_input = fluid%scn_nc - scn_nc_new = fluid%scn_nc_ps - 6 - ! scn_nc_ps : CN from which all SCN fractions will be lumped - ! into the specified number of pseudos - - characterization%plus_w = plus_w - z_m_plus_i = characterization%product_z_mw_plus_i - plus_z_i = characterization%plus_z_i - - if (scn_nc_new < scn_nc_input) then - ! Plus Fraction needs to be extended to include lower CN's - plus_zm = plus_zm + sum(scn_zm(scn_nc_new + 1 : scn_nc_input)) - plus_z = plus_z + sum(scn_z(scn_nc_new + 1 : scn_nc_input)) !extended zp - plus_w_new = plus_w + sum(w(fluid%def_comp_nc+scn_nc_new + 1 & - : fluid%def_comp_nc+scn_nc_input)) - !! use new variable for plus_w because type fluidata can't be modified. - plus_mw = plus_zm / plus_z - z_m_plus_i = scn_zm(scn_nc_new + 1 : scn_nc_input ) - z_m_plus_i = [z_m_plus_i, characterization%product_z_mw_plus_i] - plus_z_i = scn_z(scn_nc_new + 1 : scn_nc_input ) - plus_z_i = [plus_z_i, characterization%plus_z_i] - i_last = i_last + scn_nc_input - scn_nc_new - last_C = last_C - scn_nc_input + scn_nc_new - characterization%plus_w = plus_w_new - characterization%product_z_mw_plus_i = z_m_plus_i - characterization%plus_z_i = plus_z_i - end if - - !do i= 1, scn_nc_new - ! print*, scn_z(i), scn_mw(i) - !end do - - numbers_ps = fluid%numbers_ps - j_ps = 0 - ! Lumping into Nps pseudos - rec_zm = plus_zm / numbers_ps - ! Recommended value for the product z*M (proportional to weight) - ! for each pseudo according to pedersen - remain_plus_zm = plus_zm - i_ps = 1._pr - sum_z = 0.0_pr - sum_zm = 0.0_pr - sum_volume_ps = 0.0_pr - - i = 0.0_pr - - do while (i_ps < numbers_ps) - i = i + 1 - j_ps(i_ps) = j_ps(i_ps) + 1 - sum_z = sum_z + plus_z_i(i) - sum_zm = sum_zm + z_m_plus_i(i) - sum_volume_ps = sum_volume_ps + z_m_plus_i(i) / plus6_density(scn_nc_new+i) - if (z_m_plus_i(i+1) > 2*(rec_zm - sum_zm)) then - ! when adding one more would go too far - plus_z_ps = [plus_z_ps, sum_z] - plus_mw_ps = [plus_mw_ps, sum_zm/sum_z] - w_ps = [w_ps, (characterization%plus_w*(sum_zm/plus_zm))] - remain_plus_zm = remain_plus_zm - sum_zm - if (remain_plus_zm < rec_zm) numbers_ps = i_ps + 1 - carbon_number_plus(scn_nc_new + i_ps) = 6+((plus_mw_ps(i_ps)-84)/C) - density_ps = [density_ps, (sum_zm/sum_volume_ps)] - i_ps = i_ps + 1 - sum_z = 0._pr - sum_zm = 0._pr - sum_volume_ps = 0._pr - end if - end do - - plus_z_ps = [plus_z_ps, sum(plus_z_i(i+1:i_last))] - ! at this point, Nps is the order for the last NON-Asphaltenic pseudo comp. - ! (e.g. 4 if 5 is for Asphaltenes) - plus_mw_ps =[plus_mw_ps, sum(z_m_plus_i(i+1:i_last))/plus_z_ps(numbers_ps)] - w_ps = [w_ps, characterization%plus_w * ((plus_z_ps(numbers_ps))* & - (plus_mw_ps(numbers_ps))/(plus_zm))] - carbon_number_plus(scn_nc_new + numbers_ps) = 6 + & - ((plus_mw_ps(numbers_ps) - 84) / C) - sum_zm_last_ps = (sum(z_m_plus_i(i+1:i_last))) - density_ps = [density_ps, (sum_zm_last_ps)/sum((z_m_plus_i(i+1:i_last)) / (plus6_density(scn_nc_new+i+1:last)))] - j_ps(numbers_ps) = i_last - sum(j_ps(1: numbers_ps-1)) !! revisar manana - - !do i = 1, numbers_ps - ! print*, 'ps',i, plus_z_ps(i), plus_mw_ps(i), density_ps(i), j_ps(i), carbon_number_plus(scn_nc_new+i) - !end do - - plus_density = plus_zm / sum(plus_z_ps(1:numbers_ps)* & - plus_mw_ps(1:numbers_ps)/(density_ps(1:numbers_ps))) - - moles = [(fluid%def_comp_w)/(fluid%def_comp_mw)] - moles = [moles, fluid%w(def_nc+1:def_nc+scn_nc_new)/ scn_mw(1:scn_nc_new)] - moles = [moles, w_ps/plus_mw_ps ] - - characterization%mol_fraction = moles/sum(moles) ! mole fractions normalize - characterization%last_C = last_C - characterization%i_last = i_last - characterization%last = last - - ! In this point we create a new array to content final Mis and densities - ! whit pseudos. - - lumped_z = [scn_z(1:scn_nc_new), plus_z_ps] - lumped_mw = [scn_mw(1:scn_nc_new), plus_mw_ps] - lumped_densities = [plus6_density(1:scn_nc_new), density_ps ] - - - characterization%lumped_z = lumped_z - characterization%lumped_mw = lumped_mw - characterization%lumped_densities = lumped_densities - - end associate - - end subroutine lump - - subroutine get_critical_constants(fluid, characterization, eos) - !! this subroutine doing: - !! - Calculation of Tc, Pc, omega from Pedersen correlations (adapted from code "C7plusPedersenTcPcOm") - !! - Calculation of EoS parameters (taken from CubicParam) - use critical_parameters - - implicit none - type(FluidData), intent(in) :: fluid - type(FluidDataOut), intent(inout) :: characterization - character(len=*), intent(in) :: eos - real(pr), allocatable :: modified_c(:) - real(pr), allocatable :: tc(:) - real(pr), allocatable :: pc(:) - real(pr), allocatable :: om(:) - real(pr), allocatable :: m_funtion(:) - integer :: i - - associate(& - rho => characterization%lumped_densities, & - mw => characterization%lumped_mw & - ) - - allocate(modified_c(0)) - allocate(tc(0)) - allocate(pc(0)) - allocate(om(0)) - allocate(m_funtion(0)) - - call get_parameteres_for_critical(eos) - - tc = c1*rho + c2*log(mw) + c3*mw + (c4/mw) - pc = exp(d1_p + d2*(rho**(d5)) + (d3/mw) + (d4/(mw**2))) - - if (eos == "SRK") then - m_funtion = e1 + e2*mw + e3*rho + e4*(mw**(2)) - else - if (eos == "PR" .or. eos =="RKPR") then - m_funtion = e3*rho + 0.2 + 1.8*(1-exp(-mw/220)) - ! Alternative expression without max with MW - May 2019 - endif - endif - - modified_c = c0 - m_funtion - om = (-b + sqrt(b**2 - 4*a*modified_c))/(2*a) - - !do i = 1, size(mw) - ! print*, mw(i), rho(i), tc(i), pc(i), om(i), m_funtion(i) - !end do - - characterization%critical_temperature = tc - characterization%critical_pressure = pc - characterization%acentric_factor = om - characterization%m_funtion = m_funtion - - end associate - - end subroutine get_critical_constants - - type(FluidDataOut) function characterize(file, mw_source, method, fix_C, eos) & - result(characterization) - - type(FluidData) :: fluid - character(len=*), intent(in) :: file !! file name - character(len=*), intent(in), optional :: method !! plus_mw or global_mw - character(len=*), intent(in) :: mw_source - logical, intent(in), optional :: fix_C - character(len=*), intent(in) :: eos - - - - fluid = data_from_file(file) - characterization%input_data = fluid - allocate(characterization%scn_z(fluid%scn_nc)) - allocate(characterization%log_scn_z(fluid%scn_nc)) - allocate(characterization%scn_mw(fluid%scn_nc)) - allocate(characterization%carbon_number_plus(0)) - allocate(characterization%plus_z_i(0)) - allocate(characterization%product_z_mw_plus_i(0)) - allocate(characterization%scn_zm(fluid%scn_nc)) - allocate(characterization%scn_i(0)) - allocate(characterization%plus6_density(0)) - allocate(characterization%mol_fraction(0)) - allocate(characterization%lumped_z(0)) - allocate(characterization%lumped_mw(0)) - allocate(characterization%lumped_densities(0)) - allocate(characterization%critical_temperature(0)) - allocate(characterization%critical_pressure(0)) - allocate(characterization%acentric_factor(0)) - allocate(characterization%m_funtion(0)) - - - call get_c_or_m_plus(fluid=fluid, mw_source=mw_source, method=method, fix_C=fix_C, characterization= characterization) - call density_funtion(fluid=fluid, mw_source=mw_source, characterization=characterization) - call lump(fluid=fluid, characterization=characterization) - call get_critical_constants(fluid= fluid, characterization=characterization, eos=eos) - - - end function characterize - -end module routines - - diff --git a/src/tools/io.f90 b/src/tools/io.f90 old mode 100644 new mode 100755 diff --git a/src/types.f90 b/src/types.f90 old mode 100644 new mode 100755 index 9f77624..242f369 --- a/src/types.f90 +++ b/src/types.f90 @@ -7,6 +7,8 @@ module dtypes integer :: scn_nc !! number of single cuts being considered in the oil integer :: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified number of pseudos integer :: numbers_ps !! number of pseudos in which the scn fractions grouped. + logical :: density_setup + integer :: number_plus_density character(len=:), allocatable :: filename integer, allocatable :: scn (:) !! set of singles cuts being considered in the oil character(len=15), allocatable :: def_components (:) !! set of defined components being considered in the oil @@ -20,13 +22,26 @@ module dtypes real(pr) :: sum_z_mw_i !! sum of the product between composition and molecular weight of the fluid's compounds real(pr), allocatable :: w(:) !! mass fractions of the fluid's compounds real(pr) :: plus_z !! composition of residual fraction + real(pr) :: plus6_z_exp !! composition of residual fraction + real(pr) :: plus7_z_exp !! composition of residual fraction + real(pr) :: plus12_z_exp !! composition of residual fraction + real(pr) :: plus30_z_exp !! composition of residual fraction real(pr) :: plus_mw !! molecular weight of residual fraction + real(pr) :: plus6_mw_exp !! molecular weight of residual fraction + real(pr) :: plus7_mw_exp !! molecular weight of residual fraction + real(pr) :: plus12_mw_exp !! molecular weight of residual fraction + real(pr) :: plus30_mw_exp !! molecular weight of residual fraction real(pr) :: product_z_mw_plus !! product between composition and molecular weight of residual fraction real(pr), allocatable :: def_comp_w(:) !! mass fractions of the defined compounds real(pr), allocatable :: scn_w(:) !! !! mass fractions of the scn-s compounds real(pr) :: plus_w !! mass fractions of the plus fraction real(pr), allocatable :: scn_density(:) !! set of corresponding densities of scn cuts real(pr) :: plus_density !! experimental density of the plus fraction + real(pr) :: plus6_density_exp + real(pr) :: plus7_density_exp + real(pr) :: plus12_density_exp + real(pr) :: plus30_density_exp + end type FluidData type :: FluidDataOut @@ -58,7 +73,8 @@ module dtypes real(pr), allocatable :: product_z_mw_plus_i(:) real(pr) :: a_d !! ad constant which is used in equation \[rho_i = ad*exp(-i/10) +bd\] for density. real(pr) :: b_d !! bd constant which is used in equation \[rho_i = ad*exp(-i/10) +bd\] for density. - real(pr) :: volume_6plus_cal + real(pr) :: volume_cal + real(pr) :: volume_exp integer :: last_C integer :: i_last integer :: last @@ -86,7 +102,7 @@ module dtypes subroutine write_result(characterization,unit,iotype,v_list,iostat,iomsg) use ftools__io, only: str - use critical_parameters + use defined_critical_parameters implicit none class(FluidDataOut), intent(in) :: characterization From aba100f5528ac4d8f2225c3485f01099dbb9a041 Mon Sep 17 00:00:00 2001 From: oscartb25 Date: Tue, 7 Jan 2025 21:16:24 -0500 Subject: [PATCH 09/16] commit --- ford.md | 0 fort.1 | 0 2 files changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 ford.md delete mode 100644 fort.1 diff --git a/ford.md b/ford.md old mode 100644 new mode 100755 diff --git a/fort.1 b/fort.1 deleted file mode 100644 index e69de29..0000000 From 2ff312d1b470b9628ded3d299948c311522aeb73 Mon Sep 17 00:00:00 2001 From: oscartb25 Date: Tue, 7 Jan 2025 21:17:48 -0500 Subject: [PATCH 10/16] . --- YPF2 copy.nml | 78 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100755 YPF2 copy.nml diff --git a/YPF2 copy.nml b/YPF2 copy.nml new file mode 100755 index 0000000..d8ffef1 --- /dev/null +++ b/YPF2 copy.nml @@ -0,0 +1,78 @@ +! Namelist based input file +! ============================================================================== +! - YPF2 +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 14 + ! number of single carbon number components + scn_nc_ps = 20 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped + density_setup = .true. + ! Logical variable to choose the density calculation method + number_plus_density = 6 + ! CN from which we considered ....... +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 + ! names of scn fractions + scn_plus= "C20+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.00005 0.01963 0.06271 0.00707 0.01354 0.0045 0.01607 0.01069 0.01562 + ! compositions of defined components + scn_z = 0.03374 0.06051 0.08007 0.06365 0.05495 0.03993 0.03139 0.03392 0.03133 0.03183 0.03042 0.02984 0.03127 0.02964 + ! compositions of single carbon numbers + plus_z = 0.26762 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.26762 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.26762 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.26762 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.26762 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.1 72.1 + ! molecular weights of defined components + scn_mw = 84 96 107 121 134 147 161 175 190 206 222 237 251 263 + ! molecular weights of single carbon numbers + plus_mw = 543.35 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 465 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 500 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 600 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 700 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 0.685 0.749 0.768 0.793 0.808 0.815 0.836 0.85 0.861 0.873 0.882 0.873 0.875 0.885 + ! density of single carbon numbers + plus_density = 0.936 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.77 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.77 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.77 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.77 + ! density of residual fraction CN 30 plus +/ + + From c4f87955ec3b445b8ffa2e6577f2e1096ba35bf4 Mon Sep 17 00:00:00 2001 From: oscartb25 Date: Tue, 7 Jan 2025 21:18:31 -0500 Subject: [PATCH 11/16] . --- app/my_objective.mod | Bin 3632 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100755 app/my_objective.mod diff --git a/app/my_objective.mod b/app/my_objective.mod deleted file mode 100755 index 6ad6d7c2e8eb644e274bd09c78b6f7c852d9d1a7..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 3632 zcmV-04$tu)iwFP!000001KnL)bJ|E2em}pW-@p&ErPiIdC7^hz!9Zc~O!h@ZYY%-IfVz7|DbHDWI&h5wPd@=7&$d~QY*O!QVi+205#`=v^qpjL<&zu0&+yglMBiVZR)z5@--td(zLEU*_g|>lVZTHy9(2B}zK81r zsvv&GaNGuSVjpjlzJC(}>K=|{&=`CQ+$WuXBme(k3$j|W;v%b+|Z@?;rP}@1S2>YF6TknwcJgI z%W;q}WEc#_1#a9gKg_4g#b^?s%KX;Mq9Kjj@qHByUKorTQs`RfLM*lub>u0-O8>F|@3#G~w?L@Q}Y zE10onn3R|%F)d=+LPAYc>D!|KpQ8O^#QYdf<$u?QKSPx0`rC@}VxNmUzZ_WrGOce+ z1enU2?lpve)BBm?|H^xGR*yy&|Lrs&nf^}}|1;pdQbYm;cB`lx)S_+VHn41q`aR3b zv#n#Yr4yo9lsY*?Sw`Zfq3EWm+e3sVy#WjTbQnIqtoHj*#Yu;G<{ut9&?@)W40I-! z`;*aqI2a?jj4PIj*pr&23+cTrMIWFYZD&(<)kAeoGJn zgHqD2Gy@QXarE>SdV19J`PJ8!FrgEh8Z>Sn+hpq{BEw; zAtNL_=#m2MLX=TdQj8L*l4kmGsX!%BQn377OEVC?HwIy&0}o%K4Hsil=XvfJgIUmm z`NC0uTr`zRLDbe#AlEco>_S<-7iQ^yS^EmvU`s`#c^1fPoXxpKwan<0^6V&VNTHOZ zZd+B&v+LX=oKWv(ng_M+QA6e+^NxnByirVmnbfl)piHM?opZL_-)8F3wOwL+#P*5Z zBaR`ib4(f8;_9Xx-*l^+vsH2qYff|>lQO2zIDpt+Etx%-xS1Uq}`O3 zo~1L$yAao_O4E%k2eb>6#&TIAO}D4!K+e+j1~CRtH~l=RX`MTei{~QfKLW|xU9B1w zr~n}qg*7)19q{+s_CHA4|0n1Z*Zn88-G3A?GU5k^4=VU%%M((33O#sLe5eOHtsj9G zIf@({*Rf64L&tIDeh?>!nGZ)eN>?C zdUM#VH~aPDO>>6dNYBURWUxRRA;CC6b6KEy&hS17l8l{DWoJ$U!Yn539t$s+3q6PH z&H51X1`sy5#)tlV6fZ5R$Ypewlp(*|_aB0ibO@T?ZLmmqn`KRcd_~r@x&fa>sPNli z^*K!h+m;?FRZ7}cr79r>9?-2N0BaMl8iHs>P+09H7BU?PVzEBS5R)q6)Hi@4VlzW5 zZRIFX(Ndc#R6uRZq=pPjM21RJU$>=yF}c}#ad{I{avm^T zGqjjblskun-6E{!bHywQGy;bKwZdA(V^6se5Y9zjgyDK$`@WJsu(EXU)f zxt|zLDQQ$A7kY8xAFWvR=IG@T#8qn{+Ys*ID*Gqv0imk*;py1tJ6%aMDNCUd=NE?dG z9p&m&?iH)~BWh4ClGKiEoSmieNK(&FwMbHNP)B-A?MMQ3qSML{lu`gC$(^t?P;#}l zox9pL6=&OVo6P-E!^-{d8%)vbC;Z2m;?xaO!+xKHpBl~;Wc<`{#VP-LW&G6eL@9rr z2a``eypK7v-3`7lrKd(u$rGQj!YoU0g*kn*Csr(BDBI6b>xm_GjF#v?EIPGfj29i$ zQBv4yzzfe5$PTpVkS6sY-)^9q)Ps1Ng-%8<-`Y{m=*3~?5*@*l3j)ffjIOyNNH_R& z;24pp^*)6GrM9p#XoyJXtqdr2B&ET{MzI?z1L}$@14_NiB-NZJDE6R))xLhjlBCj9 z60xW}Vx2Tw$!7s#I}!ny|T$H1+SmD0yK-!3wTMMOo>WQsxo8QwkK64 z5L))6$`K6W2Q_a9DKR>@+9DWt|!h*L7*UNA z$a_PITbTsaT#+K+DRScp?Nd0KjxZCBx3W_O5*Uo-$n51^Fbwbi(Vqoc8~Rg47l_A0 z@ra30AtfgO8VC?P8lCGgZnbCSbjIu{8w@2?Wp$O5HSw(`5MLY`3NaYHr$W`5e z)yS2VQ_o^bE-hD97d@TuY|cC~>FV$kr2N!lCvWqt&_rF#@Smj~^xOl)#m^bld0bz8ycW$EDjn&h(4KT-x*Z9C{9Gv2{`)$0}HhQ=32j6e*BnbN14Xqt|0Ej>FkYk(i6yi(9mpR;+Ba z7scpuw3l9XB28|>Eymp;(I&}pGw~A@kDKB5_4-9{aTvZ9aGp~H#~KY<=n8GOlh8r6 zQ`Ae&PuafyIAQMt7L=Z{^Qu(8T>nYY5yEV~toB&dwY6gqAbMo2u4` znbO!1q*&Nc*TPqP)+S`j;IEO$8d(kQ7w_CAoFPn=ifbgy<(s=VH#tX=f9KS=Q?8N+!YVPe}^?RKjl87QYszm*4-xgqo* zi_z_Jy1Z9KnyqddoR6RPC%=joyG&H-ABQO1A9m}fNP@p6hKf%J$_aVskCy?DTrznC zO{;ROC$-+JG3}}2^f{4^r?;a4{vvIi9p5)1bt=kWzl|~owHDNM!IkKOeQmnn$#mgj zNZaZ{ujo7DLhOvbQ#4C0h6C|Juu*iKWjJSWow|t3vaz^rT8b@>gSlDTZB4crD zjz4)9V__#I3XSD*NLv~UI|Pw87S}m*qM@mpUeQ>}aL$ouXd>!WjOA1VR!uWbGZuEp zqR?3UD>D{$n4#EME{3$7v9Nz|EkECpoz+@B^3Nd0lK_h+x8cV=G#`wHm-VU#j~T=_ zpmcf{ewLr;RnxY`>r7~S02TYT1@)&n;U(OLHV!j?kWvB~b(cTRq87DnL+f2~_;ri> zx)is=O~o1{o|27MVA; Date: Tue, 7 Jan 2025 21:23:55 -0500 Subject: [PATCH 12/16] . --- .gitignore | 4 --- YPF2.nml | 78 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 78 insertions(+), 4 deletions(-) create mode 100755 YPF2.nml diff --git a/.gitignore b/.gitignore index 259de11..a01e08b 100755 --- a/.gitignore +++ b/.gitignore @@ -88,10 +88,6 @@ fig77 fpm.toml src/tools/ftools__io.mod YPF1.nml -YPF2.nml -YPF1.nml -YPF1.nml -fpm.toml fpm.toml src/molecular_weigth/molecular_weigth.mod .gitignore diff --git a/YPF2.nml b/YPF2.nml new file mode 100755 index 0000000..d8ffef1 --- /dev/null +++ b/YPF2.nml @@ -0,0 +1,78 @@ +! Namelist based input file +! ============================================================================== +! - YPF2 +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 14 + ! number of single carbon number components + scn_nc_ps = 20 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped + density_setup = .true. + ! Logical variable to choose the density calculation method + number_plus_density = 6 + ! CN from which we considered ....... +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 + ! names of scn fractions + scn_plus= "C20+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.00005 0.01963 0.06271 0.00707 0.01354 0.0045 0.01607 0.01069 0.01562 + ! compositions of defined components + scn_z = 0.03374 0.06051 0.08007 0.06365 0.05495 0.03993 0.03139 0.03392 0.03133 0.03183 0.03042 0.02984 0.03127 0.02964 + ! compositions of single carbon numbers + plus_z = 0.26762 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.26762 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.26762 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.26762 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.26762 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.1 72.1 + ! molecular weights of defined components + scn_mw = 84 96 107 121 134 147 161 175 190 206 222 237 251 263 + ! molecular weights of single carbon numbers + plus_mw = 543.35 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 465 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 500 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 600 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 700 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 0.685 0.749 0.768 0.793 0.808 0.815 0.836 0.85 0.861 0.873 0.882 0.873 0.875 0.885 + ! density of single carbon numbers + plus_density = 0.936 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.77 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.77 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.77 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.77 + ! density of residual fraction CN 30 plus +/ + + From 4d83f878bf091d374da7533fd830dceb2a65ef90 Mon Sep 17 00:00:00 2001 From: oscartb25 Date: Wed, 8 Jan 2025 11:33:37 -0500 Subject: [PATCH 13/16] 08 de enero 2025, version actualizada de namelist --- PVT5.nml | 6 +----- src/input_data/data_input.f90 | 28 ++++++++-------------------- 2 files changed, 9 insertions(+), 25 deletions(-) diff --git a/PVT5.nml b/PVT5.nml index 7ffb5b0..dbcaad7 100755 --- a/PVT5.nml +++ b/PVT5.nml @@ -1,6 +1,6 @@ ! Namelist based input file ! ============================================================================== -! - PVT5 +! - PVT5 from whitson dataset ! - Composition: molar ! - Molecular Weights: g/mol ! - Density: g/cm3 @@ -15,10 +15,6 @@ ! CN from which all SCN fractions will be lumped into the specified number of pseudos numbers_ps = 5 ! number of pseudos in which the scn fractions grouped - density_setup = .true. - ! Logical variable to choose the density calculation method - number_plus_density = 6 - ! CN from which we considered ....... / &nml_components def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" diff --git a/src/input_data/data_input.f90 b/src/input_data/data_input.f90 index 0e25663..a0b658a 100755 --- a/src/input_data/data_input.f90 +++ b/src/input_data/data_input.f90 @@ -5,17 +5,15 @@ module data_from_input implicit none contains - subroutine read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps, density_setup, number_plus_density) + subroutine read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps) !! Reads the setup data from input file character(len=*), intent(in) :: file !! file name integer, intent(out):: def_comp_nc !! number of defined components being considered in the oil integer, intent(out):: scn_nc !! number of single cuts being considered in the oil integer, intent(out):: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified number of pseudos integer, intent(out):: numbers_ps !! number of pseudos in which the scn fractions grouped - logical, intent(out) :: density_setup - integer, intent(out) :: number_plus_density integer :: funit - namelist /nml_setup/ def_comp_nc, scn_nc, scn_nc_ps, numbers_ps, density_setup, number_plus_density + namelist /nml_setup/ def_comp_nc, scn_nc, scn_nc_ps, numbers_ps open(newunit=funit, file=file) read(funit, nml=nml_setup) @@ -32,12 +30,10 @@ subroutine read_components(file, def_components, scn, scn_plus) integer :: scn_nc !! number of single cuts being considered in the oil integer :: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified number of pseudos integer :: numbers_ps !! number of pseudos in which the scn fractions grouped - logical :: density_setup - integer :: number_plus_density integer :: funit namelist /nml_components/ def_components, scn, scn_plus - call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps, density_setup, number_plus_density) + call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps) allocate(def_components(def_comp_nc)) allocate(scn(scn_nc)) @@ -60,12 +56,10 @@ subroutine read_composition(file, def_comp_z, scn_z, plus_z, plus6_z_exp, plus7_ integer :: scn_nc !! number of single cuts being considered in the oil integer :: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified number of pseudos integer :: numbers_ps !! number of pseudos in which the scn fractions grouped - logical :: density_setup - integer :: number_plus_density integer :: funit namelist /nml_composition/ def_comp_z, scn_z, plus_z, plus6_z_exp, plus7_z_exp, plus12_z_exp, plus30_z_exp - call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps, density_setup, number_plus_density) + call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps) allocate(def_comp_z(def_comp_nc)) allocate(scn_z(scn_nc)) @@ -88,12 +82,10 @@ subroutine read_molecular_weight(file, def_comp_mw, scn_mw, plus_mw, plus6_mw_ex integer :: scn_nc !! number of single cuts being considered in the oil integer :: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified number of pseudos integer :: numbers_ps !! number of pseudos in which the scn fractions grouped - logical :: density_setup - integer :: number_plus_density integer :: funit namelist /nml_molecular_weight/ def_comp_mw, scn_mw, plus_mw, plus6_mw_exp, plus7_mw_exp, plus12_mw_exp, plus30_mw_exp - call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps, density_setup, number_plus_density) + call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps) allocate(def_comp_mw(def_comp_nc)) allocate(scn_mw(scn_nc)) @@ -116,8 +108,6 @@ subroutine mass_fractions(file, w, product_z_mw_def_comp, product_z_mw_scn,& integer :: scn_nc !! number of single cuts being considered in the oil integer :: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified number of pseudos integer :: numbers_ps !! number of pseudos in which the scn fractions grouped - logical :: density_setup - integer :: number_plus_density integer :: funit real(pr), allocatable :: def_comp_z(:) !! set of corresponding mole fractions of defined components real(pr), allocatable :: scn_z(:) !! set of corresponding mole fractions of scn cuts @@ -138,7 +128,7 @@ subroutine mass_fractions(file, w, product_z_mw_def_comp, product_z_mw_scn,& real(pr), allocatable, intent(out) :: scn_w(:) real(pr), intent(out) :: plus_w - call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps, density_setup, number_plus_density) + call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps) allocate(product_z_mw_def_comp(def_comp_nc)) allocate(product_z_mw_scn(scn_nc)) allocate(def_comp_z(def_comp_nc)) @@ -183,12 +173,10 @@ subroutine read_density(file,scn_density,plus_density, plus6_density_exp, & integer :: scn_nc !! number of single cuts being considered in the oil integer :: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified number of pseudos integer :: numbers_ps !! number of pseudos in which the scn fractions grouped - logical :: density_setup - integer :: number_plus_density integer :: funit namelist /nml_density/ scn_density, plus_density, plus6_density_exp, plus7_density_exp, plus12_density_exp, plus30_density_exp - call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps, density_setup, number_plus_density) + call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps) allocate(scn_density(scn_nc)) open(newunit=funit, file=file) @@ -204,7 +192,7 @@ type(FluidData) function data_from_file(file) result(data) data%filename = file call read_setup(file, data%def_comp_nc, data%scn_nc, data%scn_nc_ps, & - data%numbers_ps, data%density_setup, data%number_plus_density & + data%numbers_ps & ) call read_components(file, data%def_components, data%scn, data%scn_plus) call read_composition(file, data%def_comp_z, data%scn_z, data%plus_z, & From 22ee3678bdced91054d01db0074c3e60058870e7 Mon Sep 17 00:00:00 2001 From: oscartb25 Date: Wed, 8 Jan 2025 11:34:57 -0500 Subject: [PATCH 14/16] eliminar archivo de YPF2 redundante --- YPF2 copy.nml | 78 --------------------------------------------------- 1 file changed, 78 deletions(-) delete mode 100755 YPF2 copy.nml diff --git a/YPF2 copy.nml b/YPF2 copy.nml deleted file mode 100755 index d8ffef1..0000000 --- a/YPF2 copy.nml +++ /dev/null @@ -1,78 +0,0 @@ -! Namelist based input file -! ============================================================================== -! - YPF2 -! - Composition: molar -! - Molecular Weights: g/mol -! - Density: g/cm3 -! -! -&nml_setup - def_comp_nc = 9 - ! number of defined components before scn - scn_nc = 14 - ! number of single carbon number components - scn_nc_ps = 20 - ! CN from which all SCN fractions will be lumped into the specified number of pseudos - numbers_ps = 5 - ! number of pseudos in which the scn fractions grouped - density_setup = .true. - ! Logical variable to choose the density calculation method - number_plus_density = 6 - ! CN from which we considered ....... -/ -&nml_components - def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" - ! names of defined components before scn - scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 - ! names of scn fractions - scn_plus= "C20+" - ! name of plus fraction -/ -&nml_composition - def_comp_z = 0.00005 0.01963 0.06271 0.00707 0.01354 0.0045 0.01607 0.01069 0.01562 - ! compositions of defined components - scn_z = 0.03374 0.06051 0.08007 0.06365 0.05495 0.03993 0.03139 0.03392 0.03133 0.03183 0.03042 0.02984 0.03127 0.02964 - ! compositions of single carbon numbers - plus_z = 0.26762 - ! composition of residual fraction CN 20 plus - plus6_z_exp = 0.26762 - ! composition of residual fraction CN 6 plus - plus7_z_exp = 0.26762 - ! composition of residual fraction CN 7 plus - plus12_z_exp = 0.26762 - ! composition of residual fraction CN 12 plus - plus30_z_exp = 0.26762 - ! composition of residual fraction CN 30 plus -/ -&nml_molecular_weight - def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.1 72.1 - ! molecular weights of defined components - scn_mw = 84 96 107 121 134 147 161 175 190 206 222 237 251 263 - ! molecular weights of single carbon numbers - plus_mw = 543.35 - ! molecular weight of residual fraction CN 20 plus - plus6_mw_exp = 465 - ! molecular weight of residual fraction CN 6 plus - plus7_mw_exp = 500 - ! molecular weight of residual fraction CN 7 plus - plus12_mw_exp = 600 - ! molecular weight of residual fraction CN 12 plus - plus30_mw_exp = 700 - ! molecular weight of residual fraction CN 30 plus -/ -&nml_density - scn_density = 0.685 0.749 0.768 0.793 0.808 0.815 0.836 0.85 0.861 0.873 0.882 0.873 0.875 0.885 - ! density of single carbon numbers - plus_density = 0.936 - ! density of residual fraction CN 20 plus - plus6_density_exp = 0.77 - ! density of residual fraction CN 6 plus - plus7_density_exp = 0.77 - ! density of residual fraction CN 7 plus - plus12_density_exp = 0.77 - ! density of residual fraction CN 12 plus - plus30_density_exp = 0.77 - ! density of residual fraction CN 30 plus -/ - - From 89b6037f26a8f61f82d5fcc12e7a3a69f7d490d7 Mon Sep 17 00:00:00 2001 From: oscartb25 Date: Fri, 7 Feb 2025 11:36:55 -0300 Subject: [PATCH 15/16] inputs loads --- PVT1.nml | 72 ++++++++++++++++++++++++++++++++ PVT2.nml | 73 +++++++++++++++++++++++++++++++++ PVT4.nml | 72 ++++++++++++++++++++++++++++++++ app/main.f90 | 1 + src/density/density_funtion.f90 | 2 +- 5 files changed, 219 insertions(+), 1 deletion(-) create mode 100755 PVT1.nml create mode 100755 PVT2.nml create mode 100755 PVT4.nml diff --git a/PVT1.nml b/PVT1.nml new file mode 100755 index 0000000..5340f70 --- /dev/null +++ b/PVT1.nml @@ -0,0 +1,72 @@ +! Namelist based input file +! ============================================================================== +! - PVT1 from whitson dataset +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 24 + ! number of single carbon number components + scn_nc_ps = 30 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 + ! names of scn fractions + scn_plus= "C30+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.0088 0.0052 0.6037 0.1259 0.0728 0.0146 0.0308 0.0117 0.0143 + ! compositions of defined components + scn_z = 0.0189 0.0154 0.0192 0.0121 0.0088 0.0067 0.0054 0.0048 0.0037 0.0029 0.0022 0.0018 0.0016 0.0013 0.0010 0.0008 0.0007 0.0006 0.0005 0.0004 0.0004 0.0003 0.0003 0.0002 + ! compositions of single carbon numbers + plus_z = 0.0022 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.1122 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.0933 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.0311 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.0022 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.2 72.2 + ! molecular weights of defined components + scn_mw = 86.2 96.4 108.7 123.9 138.3 146.9 160.3 174.8 190.7 207.8 223.7 241.6 247.9 261.1 267 290.8 299.7 298.8 320.4 337.4 350.7 364 377.3 390.6 + ! molecular weights of single carbon numbers + plus_mw = 561.11 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 146.72 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 157.58 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 243.55 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 561.11 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + ! density of single carbon numbers + plus_density = 0.9816 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.7969 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.8113 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.8596 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.9816 + ! density of residual fraction CN 30 plus +/ diff --git a/PVT2.nml b/PVT2.nml new file mode 100755 index 0000000..42be41b --- /dev/null +++ b/PVT2.nml @@ -0,0 +1,73 @@ +! Namelist based input file +! ============================================================================== +! - PVT2 from whitson dataset +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 24 + ! number of single carbon number components + scn_nc_ps = 30 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 + ! names of scn fractions + scn_plus= "C30+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.0021 0.0035 0.6048 0.1317 0.0604 0.0125 0.0279 0.0120 0.0148 + ! compositions of defined components + scn_z = 0.0188 0.0186 0.0233 0.0149 0.0107 0.0076 0.0057 0.0054 0.0041 0.0030 0.0024 0.0020 0.0018 0.0016 0.0012 0.0011 0.0009 0.0008 0.0007 0.0006 0.0006 0.0005 0.0004 0.0004 + ! compositions of single carbon numbers + plus_z = 0.0032 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.1303 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.1114 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.0364 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.0032 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.2 72.2 + ! molecular weights of defined components + scn_mw = 85.92 95.74 107.97 123.86 138.88 146.65 162.24 174.24 189.11 207.29 220.29 238.25 255.75 262.48 282.67 279 309.59 312.96 334.59 346.06 359.99 373.92 387.86 401.79 + ! molecular weights of single carbon numbers + plus_mw = 553.6 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 150.4 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 159.93 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 253.63 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 553.6 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + ! density of single carbon numbers + plus_density = 0.9807 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.8022 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.8146 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.8654 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.9807 + ! density of residual fraction CN 30 plus +/ + diff --git a/PVT4.nml b/PVT4.nml new file mode 100755 index 0000000..d3a0602 --- /dev/null +++ b/PVT4.nml @@ -0,0 +1,72 @@ +! Namelist based input file +! ============================================================================== +! - PVT4 from whitson dataset +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 24 + ! number of single carbon number components + scn_nc_ps = 30 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 + ! names of scn fractions + scn_plus= "C30+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.0054 0.0023 0.7261 0.1179 0.0532 0.0112 0.0176 0.011 0.0074 + ! compositions of defined components + scn_z = 0.00978 0.00778 0.01038 0.00559 0.00369 0.00249 0.0017 0.0016 0.0011 0.0009 0.0006 0.0005 0.0004 0.0003 0.00021 0.00016 0.00013 0.0001 0.00008 0.00007 0.00005 0.00004 0.00003 0.00003 + ! compositions of single carbon numbers + plus_z = 0.0002 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.0478 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.0381 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.0082 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.0002 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.2 72.2 + ! molecular weights of defined components + scn_mw = 85.71 96.64 108.98 124.03 137.52 145.38 161.53 169.94 186.01 197.43 233.32 242.3 255.76 269.22 279.99 294.34 308.69 323.04 337.39 351.74 366.08 380.43 394.78 409.13 + ! molecular weights of single carbon numbers + plus_mw = 551.2 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 126.64 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 134.38 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 217.35 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 551.2 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + ! density of single carbon numbers + plus_density = 0.9804 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.7728 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.7889 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.8414 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.9804 + ! density of residual fraction CN 30 plus +/ diff --git a/app/main.f90 b/app/main.f90 index 98886dc..bc7f772 100755 --- a/app/main.f90 +++ b/app/main.f90 @@ -11,6 +11,7 @@ program main prueba = characterize(file='PVT5.nml', mw_source="experimental", method = "plus_mw", pho_method = 3 , fix_C=.true., eos='PR') write(*, *) prueba print*, prueba%C, prueba%a, prueba%b, prueba%plus_mw, prueba%plus_z + print*, prueba%c_max !print*, prueba%a , prueba%b , prueba%n_init, prueba%input_data%number_plus_density !do i = 1, size(prueba%mol_fraction) diff --git a/src/density/density_funtion.f90 b/src/density/density_funtion.f90 index 5ce4e14..32eb4b1 100755 --- a/src/density/density_funtion.f90 +++ b/src/density/density_funtion.f90 @@ -66,7 +66,7 @@ subroutine density_method (fluid, pho_method,characterization) integer, intent(in) :: pho_method real(pr) :: volume_exp - select case (pho_method) + select case (pho_method) ! cambiar a rho case (1) ! case 1 use experimental scn's densities reported to calculate From 7dbd7372024cb02e0af88634d462845604b4d204 Mon Sep 17 00:00:00 2001 From: Oscar Theran Date: Fri, 7 Feb 2025 12:40:02 -0300 Subject: [PATCH 16/16] 2025 --- .github/workflows/docs.yml | 0 .gitignore | 0 PVT1.nml | 72 +++++++++++++++++ PVT2.nml | 73 ++++++++++++++++++ PVT4.nml | 72 +++++++++++++++++ PVT5.nml | 0 README.md | 0 YPF1.nml | 0 YPF2.nml | 0 app/main.f90 | 1 + app/optimize.f90 | 0 app/pedersen.f90 | 0 example/demo.f90 | 0 ford.md | 0 fpm.toml | 0 input.nml | 0 oil1.nml | 0 src/bfr_routines/bfr_routines.f90 | 0 src/characterization/characterization.f90 | 0 src/constans.f90 | 0 .../critical_parameters.f90 | 0 .../defined_critical_param.f90 | 0 src/density/density_funtion.f90 | 2 +- src/input_data/data_from_input.mod | Bin 3919 -> 0 bytes src/input_data/data_input.f90 | 0 src/lumping/lumping_method.f90 | 0 src/molecular_weigth/mw_calculations.f90 | 0 src/tools/io.f90 | 0 src/types.f90 | 0 29 files changed, 219 insertions(+), 1 deletion(-) mode change 100755 => 100644 .github/workflows/docs.yml mode change 100755 => 100644 .gitignore create mode 100644 PVT1.nml create mode 100644 PVT2.nml create mode 100644 PVT4.nml mode change 100755 => 100644 PVT5.nml mode change 100755 => 100644 README.md mode change 100755 => 100644 YPF1.nml mode change 100755 => 100644 YPF2.nml mode change 100755 => 100644 app/main.f90 mode change 100755 => 100644 app/optimize.f90 mode change 100755 => 100644 app/pedersen.f90 mode change 100755 => 100644 example/demo.f90 mode change 100755 => 100644 ford.md mode change 100755 => 100644 fpm.toml mode change 100755 => 100644 input.nml mode change 100755 => 100644 oil1.nml mode change 100755 => 100644 src/bfr_routines/bfr_routines.f90 mode change 100755 => 100644 src/characterization/characterization.f90 mode change 100755 => 100644 src/constans.f90 mode change 100755 => 100644 src/critical_parameters/critical_parameters.f90 mode change 100755 => 100644 src/critical_parameters/defined_critical_param.f90 mode change 100755 => 100644 src/density/density_funtion.f90 delete mode 100755 src/input_data/data_from_input.mod mode change 100755 => 100644 src/input_data/data_input.f90 mode change 100755 => 100644 src/lumping/lumping_method.f90 mode change 100755 => 100644 src/molecular_weigth/mw_calculations.f90 mode change 100755 => 100644 src/tools/io.f90 mode change 100755 => 100644 src/types.f90 diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml old mode 100755 new mode 100644 diff --git a/.gitignore b/.gitignore old mode 100755 new mode 100644 diff --git a/PVT1.nml b/PVT1.nml new file mode 100644 index 0000000..5340f70 --- /dev/null +++ b/PVT1.nml @@ -0,0 +1,72 @@ +! Namelist based input file +! ============================================================================== +! - PVT1 from whitson dataset +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 24 + ! number of single carbon number components + scn_nc_ps = 30 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 + ! names of scn fractions + scn_plus= "C30+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.0088 0.0052 0.6037 0.1259 0.0728 0.0146 0.0308 0.0117 0.0143 + ! compositions of defined components + scn_z = 0.0189 0.0154 0.0192 0.0121 0.0088 0.0067 0.0054 0.0048 0.0037 0.0029 0.0022 0.0018 0.0016 0.0013 0.0010 0.0008 0.0007 0.0006 0.0005 0.0004 0.0004 0.0003 0.0003 0.0002 + ! compositions of single carbon numbers + plus_z = 0.0022 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.1122 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.0933 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.0311 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.0022 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.2 72.2 + ! molecular weights of defined components + scn_mw = 86.2 96.4 108.7 123.9 138.3 146.9 160.3 174.8 190.7 207.8 223.7 241.6 247.9 261.1 267 290.8 299.7 298.8 320.4 337.4 350.7 364 377.3 390.6 + ! molecular weights of single carbon numbers + plus_mw = 561.11 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 146.72 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 157.58 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 243.55 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 561.11 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + ! density of single carbon numbers + plus_density = 0.9816 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.7969 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.8113 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.8596 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.9816 + ! density of residual fraction CN 30 plus +/ diff --git a/PVT2.nml b/PVT2.nml new file mode 100644 index 0000000..42be41b --- /dev/null +++ b/PVT2.nml @@ -0,0 +1,73 @@ +! Namelist based input file +! ============================================================================== +! - PVT2 from whitson dataset +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 24 + ! number of single carbon number components + scn_nc_ps = 30 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 + ! names of scn fractions + scn_plus= "C30+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.0021 0.0035 0.6048 0.1317 0.0604 0.0125 0.0279 0.0120 0.0148 + ! compositions of defined components + scn_z = 0.0188 0.0186 0.0233 0.0149 0.0107 0.0076 0.0057 0.0054 0.0041 0.0030 0.0024 0.0020 0.0018 0.0016 0.0012 0.0011 0.0009 0.0008 0.0007 0.0006 0.0006 0.0005 0.0004 0.0004 + ! compositions of single carbon numbers + plus_z = 0.0032 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.1303 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.1114 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.0364 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.0032 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.2 72.2 + ! molecular weights of defined components + scn_mw = 85.92 95.74 107.97 123.86 138.88 146.65 162.24 174.24 189.11 207.29 220.29 238.25 255.75 262.48 282.67 279 309.59 312.96 334.59 346.06 359.99 373.92 387.86 401.79 + ! molecular weights of single carbon numbers + plus_mw = 553.6 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 150.4 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 159.93 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 253.63 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 553.6 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + ! density of single carbon numbers + plus_density = 0.9807 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.8022 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.8146 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.8654 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.9807 + ! density of residual fraction CN 30 plus +/ + diff --git a/PVT4.nml b/PVT4.nml new file mode 100644 index 0000000..d3a0602 --- /dev/null +++ b/PVT4.nml @@ -0,0 +1,72 @@ +! Namelist based input file +! ============================================================================== +! - PVT4 from whitson dataset +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 24 + ! number of single carbon number components + scn_nc_ps = 30 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 + ! names of scn fractions + scn_plus= "C30+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.0054 0.0023 0.7261 0.1179 0.0532 0.0112 0.0176 0.011 0.0074 + ! compositions of defined components + scn_z = 0.00978 0.00778 0.01038 0.00559 0.00369 0.00249 0.0017 0.0016 0.0011 0.0009 0.0006 0.0005 0.0004 0.0003 0.00021 0.00016 0.00013 0.0001 0.00008 0.00007 0.00005 0.00004 0.00003 0.00003 + ! compositions of single carbon numbers + plus_z = 0.0002 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.0478 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.0381 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.0082 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.0002 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.2 72.2 + ! molecular weights of defined components + scn_mw = 85.71 96.64 108.98 124.03 137.52 145.38 161.53 169.94 186.01 197.43 233.32 242.3 255.76 269.22 279.99 294.34 308.69 323.04 337.39 351.74 366.08 380.43 394.78 409.13 + ! molecular weights of single carbon numbers + plus_mw = 551.2 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 126.64 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 134.38 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 217.35 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 551.2 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + ! density of single carbon numbers + plus_density = 0.9804 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.7728 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.7889 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.8414 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.9804 + ! density of residual fraction CN 30 plus +/ diff --git a/PVT5.nml b/PVT5.nml old mode 100755 new mode 100644 diff --git a/README.md b/README.md old mode 100755 new mode 100644 diff --git a/YPF1.nml b/YPF1.nml old mode 100755 new mode 100644 diff --git a/YPF2.nml b/YPF2.nml old mode 100755 new mode 100644 diff --git a/app/main.f90 b/app/main.f90 old mode 100755 new mode 100644 index 98886dc..bc7f772 --- a/app/main.f90 +++ b/app/main.f90 @@ -11,6 +11,7 @@ program main prueba = characterize(file='PVT5.nml', mw_source="experimental", method = "plus_mw", pho_method = 3 , fix_C=.true., eos='PR') write(*, *) prueba print*, prueba%C, prueba%a, prueba%b, prueba%plus_mw, prueba%plus_z + print*, prueba%c_max !print*, prueba%a , prueba%b , prueba%n_init, prueba%input_data%number_plus_density !do i = 1, size(prueba%mol_fraction) diff --git a/app/optimize.f90 b/app/optimize.f90 old mode 100755 new mode 100644 diff --git a/app/pedersen.f90 b/app/pedersen.f90 old mode 100755 new mode 100644 diff --git a/example/demo.f90 b/example/demo.f90 old mode 100755 new mode 100644 diff --git a/ford.md b/ford.md old mode 100755 new mode 100644 diff --git a/fpm.toml b/fpm.toml old mode 100755 new mode 100644 diff --git a/input.nml b/input.nml old mode 100755 new mode 100644 diff --git a/oil1.nml b/oil1.nml old mode 100755 new mode 100644 diff --git a/src/bfr_routines/bfr_routines.f90 b/src/bfr_routines/bfr_routines.f90 old mode 100755 new mode 100644 diff --git a/src/characterization/characterization.f90 b/src/characterization/characterization.f90 old mode 100755 new mode 100644 diff --git a/src/constans.f90 b/src/constans.f90 old mode 100755 new mode 100644 diff --git a/src/critical_parameters/critical_parameters.f90 b/src/critical_parameters/critical_parameters.f90 old mode 100755 new mode 100644 diff --git a/src/critical_parameters/defined_critical_param.f90 b/src/critical_parameters/defined_critical_param.f90 old mode 100755 new mode 100644 diff --git a/src/density/density_funtion.f90 b/src/density/density_funtion.f90 old mode 100755 new mode 100644 index 5ce4e14..32eb4b1 --- a/src/density/density_funtion.f90 +++ b/src/density/density_funtion.f90 @@ -66,7 +66,7 @@ subroutine density_method (fluid, pho_method,characterization) integer, intent(in) :: pho_method real(pr) :: volume_exp - select case (pho_method) + select case (pho_method) ! cambiar a rho case (1) ! case 1 use experimental scn's densities reported to calculate diff --git a/src/input_data/data_from_input.mod b/src/input_data/data_from_input.mod deleted file mode 100755 index 4031f7fdb67fa1d2ce1d564960eb7801e0de2e38..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 3919 zcmV-V53ukbiwFP!000001KnLoliJ7^GPC;A$H{Cl>y7FDxVi2e`rSF5x5oqRh-!zf zPtNM%Z0PIhxTl-dWo2&vn^~b!caE`Ss?U zJ5&mP9(>8P^D+3~@p@rC+v!k|cDn;>Hhi%LbUFTOJo$V4$8fx`#*6UfbT;V||3&a7 z72v<#Xf)~D{mJtbV09JRafP-$K;cs*_p<^CUV{LPWWgXGM_I!MF9HTCup4W3hXp?f^jTl7AREIJrITjTj~GEOdk9805}-SMj(z$EE7X)aDu0B4d9 z6F3PEr%pR>yX)EB|43J6_=l9!9gB%Sy5WzWF_ZTj+@zg#fpM$P1K%=4fwzw@(ybo9bG{TcVGZX_chTX1Yp4R>h$G!un*V$7v~hH(#D} zC_q+11kTU5RuBLm?^frFy>5?&WT4}%jVM~gpnG~~5rN_D55X}VoFA?EM%naWWRq2# zUP&c#sJcs41kVj<_QpA!x0he(ab+P_P&ULZPOPM>h`!4zavIch)U5>D`9M#kI)um! z@^S_#5P)R&$DRFpTCFeJ<01W6q5Y`>HU=(!b2d5WplL8`Z~}{ybi5f&6xI&gUfFB08?)M@p#} zk3H7T;d0up?N?xgW9oSusifUpmMUYwsz5q>`}KO@5Z41)?h~CZKa7U`P%S<}A$**d z9TKQ_b}D%JpKizAfME9P@w8uEF0e{}fF+^z_>Y^sfJnoEtEp)Z#`#Dx7XM@^y!JO;r^%4I&j)^hdG5 z#;-e0q8rs{4KP6*)ARFWupC*QJq)b=D5B%}^22PhTnxt+sQd5U=eceDGX>d)i-?`Y z73v<$tC)c9sxnm-s;X4gsH#)d@WFQ(IRHP{y=|`{!USLJ0=L<`y`7%9FS$(>nan*M z&h7pPbdHlvG?=J({$BRQ%FExE-YZWh_{kYqrf&xBlZgh?Xn%jZoEQumFOPpW5Xqo7 zf6dU_nqh2|_%;f&D=}3V_bS2AAdcENulD-5e}g?ja(_Vh`;_XgJ|zpQ1lyF>&Qmot z+$p^kJ-};zqom;>qI+s~ZWV@(cHu_%_Sv~y@f2Y+`84ch8L+jt1WY2q)R~fR0K<6T04qomSk{ODJK$Nl}FiL2?KvMNE$U=)4w03Jm^}rN? zD%w(k5g22O7|hTd+XG<=p*^+;LgW-Iv`~d?SeeeI1YpOq}(zNlLW-vE)XIL zH+pUz3Z9+^Ct5cA_4V?NIXVikG+wVld$R>*Y=*tj6Nwz?ivsYqyzZDR#8sd$onBEh$RuZ+IikjQd zdN1G_ua`9ox8mUjaGB=8&;H~)u$RT>sRiS~8iPpL{=|oFU+GnJ1tD4h_|a!QzQUcXT9Y3|;+J z19o^)d3f%j5O;W&h>vd=z1%)4igKI<8`<3xCU0(yEC{a_;QywWeaC-dCo`{6k493f z0%%ZJ;JO!V$^zGNaI14;)J_C(I@q2J%4!nrCH8Rr-GkM&k-MA^= z71U=aY+rk1h66ZX8d!%=#Kqhel@<`+Sfo+U0>7*VbwSyw7Ff`VOIJV0i^ zUW8>YDQWkOph;yf5-S2K#%77-P55jkiCi`lB`H7N? zYQNquLy*JMg=^Isp=s<8Er(lz&NvBSRM-bI#M z`AwM0Zc{45B2=91{r@N)H8puww{J_+xz z3LECaU?LV(uamRdL`qX1X`)osohTJ1NQDVfVS-eXFAH{P++@LiIqW)m zBw;ks6$Zj!QIp{dtWN{$(?B>4gwsGcEh!v9;q*YbFw6Q#g6W=MItZqNU^)nYnmbph{c-24ue;UB{^4xDe#@f&P`daArK)r!kLEuUgkXS#21sat zga$}xfP@A}Xc&ARGy)0z-h#VLS}t1?(e%VHK@1bbFhL9x#4te&6T~oe?j#}At-r01 z#+}pl_sU;Y?jKcRGfKm8uG^!7<@57r@A+~10qM$#!kKmhQS8rO(==QNM$;3I`6=v+Jih%( zPokPe7QI9;XVGT#2%&&aFb$7xOezTipZ_HxgrA*fk%UZ+%oc(o#f;83V=eLZcKYC3n|!X3dwQMef-I7m>fk!`n){Z1* zBigJ~(KLuq+z5CmOm5YCb=B1YqYBZ|_v$KnkJX}zVg37KFKbo0-M>{GX|0LXVnhsE zu$G8e)subPTu^n10;ICXy=QAQoG-YzNp80|H!g_$zKarkqZ)preJN>+8f#M83@OTy zyK1r>>DAE}n-$u?Zy`xFKCE|$tOHGKYaK)`Q~8|&s4mwgwVg^Gq_90`x90hZpv#(a z1KwU!Z@-kHI_e7e_>-^28gVN@BB8&p2gyAg*YsFt?gKzF;;p(5fXR%v?zWFfkF%QI z8msAYtmTf2sE~G=U>RzYQRWDSsLYa={*22opRgu z9srra@8WA88Qw_|;SWP-i$4r0dg%=x(L?fwA>VGnsf4Ukg#) z0Gc9ULyKP%@c@|g#I+HnK;4H&jKnn(?SnGotXqoYnV_3dKWxEGW+}Q!qydoyKHYzvi`0xQ zma@E&W174ks|?c=h1g=LLLto}lH!X+th;~Rb5mir;C zwpgAB>9RM}3smA>Z6i?j(wk(-mymhMj9XHKkI?xyqu`Gel5sJCKr?SH@Sm?qgzZ|u zt~VvwErf5YoN8r!-wX?%z$o3`nAZphY)^^ze&E}Z*E4~w6*5V79pL*4uYw