! optimizes process described in chapter 12 ! employing the modified model program chapt12 implicit none integer :: option, maxnx, i, j integer :: nContour ! useful in contour plots parameter(maxnx = 100) real :: DT, T, P ! these are the design variables real :: Tw1, T1, Tw2, Tc real :: ec, Tf, Xf, F, Pf, kexp real :: T0, EF, Qf, Fw, Af, Uf, Qe, Fr, Ue, Ae, Qc, Fc, Ac, Uc, Er, X, Fp real :: Cw, Ce, Ca, b1, b2, c1, c2, e, ty real :: inp(22), des(3), modl(14) real :: lowerF, higherF, deltaF, eps, x1l, x1u, x2l, x2u, x1(maxnx), x2(maxnx) ! useful in contour plots real :: objective common /arrays/ inp, des, modl external objective call read_data(DT, T, Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, P, Pf, kexp, & Uf, Ue, Uc, Cw, Ce, Ca, b1, b2, c1, c2, e, ty) 1 write(*,*) 'Options: ' write(*,*) '0 = solve without optimizing,' write(*,*) '1 = optimize using direct search,' write(*,*) '2 = optimize using steepest descent,' write(*,*) '3 = plot contour lines.' read (*,*) option if(option == 0) then call model(Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc, & P, T, DT, & Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X) call results(P, T, DT, T0, EF, Qf, Fw, Af, Qe, Fr, Ae, Qc, Fc, Ac, Er, X, Fp, & Cw, Ce, Ca, b1, b2, c1, c2, e, ty) stop endif ! call put_data(Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc, & Ca, Cw, Ce, ty, e, b1, b2, c1, c2, inp, & P, T, DT, des, & Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X, modl) ! if(option == 1) then call direct_search(inp, des, modl) else if(option == 2) then call steepest_descent(inp, des, modl) else if(option == 3) then write(*,*) 'low function value?' read(*,*) lowerF write(*,*) 'high function value?' read(*,*) higherF write(*,*) 'how many contour lines?' read(*,*) nContour deltaF=(higherF-lowerF)/(nContour) open(unit=20,file="contour.out") Do i=0,nContour-1 x1 = 0. x2 = 0. call contour(lowerF+i*deltaF, 1.e-3, -90., 20., 100, 0., 25., 100, x1, x2, maxnx) write(20,'("# data set No ",i2)') i+1 Do j = 1, maxnx write(20,*) j, x1(j), x2(j), lowerF+i*deltaF End Do write(20,'(" ")') write(20,'(" ")') End Do close(unit=20) else write(*,*) 'wrong option!' goto 1 endif end program chapt12 ! solves set of equations that comprises our model, ! given a tuple of design variables assigned specific values subroutine model(Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc, & ! input data P, T, DT, & ! design variables Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X) ! other model variables implicit none real :: DT, T, P ! these are the design variables real :: Te, Tw1, T1, Tw2, Tc real :: ec, M, Tf, Xf, F, Pf, kexp real :: EF, T0, Qf, enthalpy, Fw, Af, Uf, DTL, Qe, Fr real :: Ae, Ue, Qc, Fc, Ac, Uc, Er, X, Fp real :: R, c2k, Cpw, DHr, Cpr, X_sat, mw_air, mw_ac parameter(R = 8.3144, Cpw = 4.18, DHr = 490., Cpr = 2.45, c2k = 273.15) parameter(mw_air = 28.8, mw_ac = 58.) external enthalpy, DTL, X_sat ! ! Designer*s options; eqs 15, 16 and 17 Te = T - DT T1 = Tw1 + DT Tc = Tw2 + DT ! ! F Compressor M = (mw_air + Xf * mw_ac) / ( 1. + Xf) if(P < Pf) P = Pf ! safeguard against decompression EF = (1.0/ec)*(R/M)*(c2k+Tf)*(1.0+Xf)*F*((P/Pf)**kexp-1)/kexp ! eq 1 T0 = (c2k+Tf)*(P/Pf)**kexp-c2k ! eq 2 ! ! F Heat Exchanger Qf = F * (enthalpy(P, T0, Xf) - enthalpy(P, T1, Xf)) ! eq 3 Fw = Qf / Cpw / (Tw2 - Tw1) ! eq 4 Af = Qf / Uf / DTL(T1, Tw1, T0, Tw2) ! eq 5 ! ! E Heat Exchanger Qe = F * (enthalpy(P, T1, Xf) - enthalpy(P, T, Xf)) ! eq 6 Fr = Qe / (DHr - Cpr * (Tc - Te)) ! eq 7 Ae = Qe / Ue / DTL(T1, Te, T, Te) ! eq 8 ! ! C Heat Exchanger Qc = Fr*(DHr*(Tc+c2k)/(Te+c2k)-Cpr*(Tc-Te)) ! eq 9 Fc = Qc / Cpw / (Tw2 - Tw1) ! eq 10 Ac = Qc / Uc / DTL(Tc, Tw2, Tc, Tw1) ! eq 11 ! ! R compressor Er = (1.0/ec) * DHr * Fr * (Tc - Te) / (Te+c2k) ! eq 12 ! ! Separation vessel X = X_sat(P, T) ! eq 13 Fp = F * (Xf - X) ! eq 14 return end subroutine model !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function enthalpy(P, T, X) ! calculates the enthalpy of dry air-acetone mixture implicit none real :: enthalpy, P, T, X real :: Xv, Xl, P0ac, P0eth, X_sat real :: Cpa, Cpv, Cpl, DH0 parameter(Cpa = 1.01, Cpv = 1.18, Cpl = 2.18, DH0 = 502.0) external X_sat ! vapor and liquid concentrations Xv = amin1(X_sat(P, T), X) Xl = X - Xv enthalpy = Cpa * T + Xv * (DH0 + Cpv * T) + Xl * Cpl * T return end function enthalpy !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function X_sat(P, T) ! acetone saturation concentration (kg ac / kg dry air) implicit none real :: P, T, X_sat real :: P0, A, B, C, lambda, c2k parameter(A = 10.0311, B = 2940.5, C = 237.20) parameter(lambda = 2.011) ! MW acetone / MW air - check this! ! Vapor pressure from Antoine equation - T in Celsius P0 = exp(A - B /(C + T)) ! acetone ! Saturation concentration (kg ac / kg dry air) X_sat = lambda * P0 / (P - P0) if(X_sat < 0.) X_sat = 1000. return end function X_sat !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function DTL(t1, t2, t3, t4) implicit none real :: DTL, t1, t2, t3, t4 ! DTL = (t1-t2-t3+t4)/alog((t1-t2)/(t3-t4)) ! return end function DTL !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_data(DT, T, Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, P, Pf, kexp, & Uf, Ue, Uc, Cw, Ce, Ca, b1, b2, c1, c2, e, ty) ! implicit none real :: DT, T, P ! these are design variables real :: Tw1, T1, Tw2, Tc real :: ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc, Cw, Ce, Ca, b1, b2, c1, c2, e, ty ! open(unit=1,file='acetone.in') read(1,*) read(1,*) read(1,*) P ! design variable: assign initial value read(1,*) T ! design variable: assign initial value read(1,*) DT ! design variable: assign initial value read(1,*) Tw1 read(1,*) T1 read(1,*) Tw2 read(1,*) Tc read(1,*) ec read(1,*) Tf read(1,*) Xf read(1,*) F read(1,*) Pf read(1,*) kexp read(1,*) Uf read(1,*) Ue read(1,*) Uc read(1,*) Cw read(1,*) Ce read(1,*) Ca read(1,*) e read(1,*) b1 read(1,*) b2 read(1,*) c1 read(1,*) c2 read(1,*) ty close(1) return end subroutine read_data !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine results(P, T, DT, T0, EF, Qf, Fw, Af, Qe, Fr, Ae, Qc, Fc, Ac, Er, X, Fp, & Cw, Ce, Ca, b1, b2, c1, c2, e, ty) implicit none real :: P, T, DT, T0, EF, Qf, Fw, Af, Qe, Fr, Ae, Qc, Fc, Ac, Er, X, Fp real :: Cw, Ce, Ca, b1, b2, c1, c2, e, ty real :: fix_cost, oper_cost, tot_cost, sales, profit ! open(unit=10,file='acetone.out') write(10,1) P write(10,2) T write(10,3) DT 1 format(t1,'Design variables'/'operation pressure : ',f8.3,' bar') 2 format('Acetone condensation temperature : ',f8.3,' C') 3 format('Heat exchangers DT : ',f8.3,' C') write(10,4) Fp 4 format(/'Acetone produced : ',e8.3,' kg/s') write(10,5) 5 format(/'Equipment required'/t1,'Type',t30,'Size',t45,'Power',t60,'Cost') write(10,6) Af, Qf, b1*Af**b2 write(10,7) Ae, Qe, b1*Ae**b2 write(10,8) Ac, Qc, b1*Ac**b2 write(10,9) Ef, c1*Ef**c2 write(10,10) Er, c1*Er**c2 fix_cost = b1*Af**b2+b1*Ae**b2+b1*Ac**b2+c1*Ef**c2+c1*Er**c2 write(10,11) fix_cost 6 format(t1,'F heat exchanger',t30,f5.1,' m^2',t45,f5.1,' kW',t60,f5.1,' × 10^6 δρχ') 7 format(t1,'E heat exchanger',t30,f5.1,' m^2',t45,f5.1,' kW',t60,f5.1,' × 10^6 δρχ') 8 format(t1,'C heat exchanger',t30,f5.1,' m^2',t45,f5.1,' kW',t60,f5.1,' × 10^6 δρχ') 9 format(t1,'F compressor',t45,f5.1,' kW',t60,f5.1,' × 10^6 δρχ') 10 format(t1,'R compressor',t45,f5.1,' kW',t60,f5.1,' × 10^6 δρχ') 11 format(/t1,'total fixed cost',t60,f5.1,' × 10^6 δρχ') write(10,12) 12 format(/'Facilities required'/t1,'Type',t45,'Quantity',t60,'Cost') write(10,13) Fw+Fc, 3600.*ty*(Cw/1000.)*(Fw+Fc)/1.e+6 ! X 3600 because rate is in s, div by 1000 for cost is per m^3 write(10,14) Er+Ef, ty*Ce*(Er+Ef)/1.e+6 ! no 3600 factor because cost is per kWh oper_cost = ty*(3600.*(Cw/1000.)*(Fw+Fc)+Ce*(Er+Ef))/1.e+6 write(10,15) oper_cost 13 format(t1,'Cooling water',t45,f5.1,' kg/s',t60,f5.1,' × 10^6 δρχ') 14 format(t1,'Power (F+R)',t45,f5.1,' kW',t60,f5.1,' × 10^6 δρχ') 15 format(/t1,'total operating cost',t60,f5.1,' × 10^6 δρχ') tot_cost = e*fix_cost + oper_cost write(10,16) 16 format(/t1,'Annual balance') sales = 3600*ty*Fp*Ca/1.e+6 profit = sales - tot_cost write(10,17) sales write(10,18) oper_cost write(10,19) e*fix_cost write(10,20) profit close(unit=10) 17 format(t1,'Acetone sales',t60,f5.1,' × 10^6 δρχ') 18 format(t1,'Operating cost',t60,f5.1,' × 10^6 δρχ') 19 format(t1,'Fixed cost',t60,f5.1,' × 10^6 δρχ') 20 format(t1,'Profit',t60,f5.1,' × 10^6 δρχ') return end subroutine results !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function objective(Fp, Fw, Fc, Af, Ae, Ac, Ef, Er, Ca, Cw, Ce, ty, e, b1, b2, c1, c2) implicit none real :: objective real :: Fp, Ca, Cw, Fw, Fc, Ce, Ef, Er, ty, e, b1, b2, Af, Ae, Ac, c1, c2 real :: Ceq, Cop, Cac Ceq = e*(b1*Af**b2+b1*Ae**b2+b1*Ac**b2+c1*Ef**c2+c1*Er**c2) Cop = ty*(3600.*(Cw/1000.)*(Fw+Fc)+Ce*(Er+Ef))/1.e+6 Cac = 3600*ty*Fp*Ca/1.e+6 objective = Cac - Cop - Ceq ! return end function objective !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine put_data(Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc, & Ca, Cw, Ce, ty, e, b1, b2, c1, c2, inp, & P, T, DT, des, & Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X, modl) implicit none real :: Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc real :: Ca, Cw, Ce, ty, e, b1, b2, c1, c2, inp(22) real :: P, T, DT, des(3) real :: Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X, modl(14) inp(1) = Tw1 inp(2) = T1 inp(3) = Tw2 inp(4) = Tc inp(5) = ec inp(6) = Tf inp(7) = Xf inp(8) = F inp(9) = Pf inp(10) = kexp inp(11) = Uf inp(12) = Ue inp(13) = Uc inp(14) = Ca inp(15) = Cw inp(16) = Ce inp(17) = ty inp(18) = e inp(19) = b1 inp(20) = b2 inp(21) = c1 inp(22) = c2 ! des(1) = P des(2) = T des(3) = DT ! modl(1) = Fp modl(2) = Fr modl(3) = Fw modl(4) = Fc modl(5) = Af modl(6) = Ae modl(7) = Ac modl(8) = Ef modl(9) = Er modl(10) = T0 modl(11) = Qf modl(12) = Qe modl(13) = Qc modl(14) = X return end subroutine put_data !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine get_data(Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc, & Ca, Cw, Ce, ty, e, b1, b2, c1, c2, inp, & P, T, DT, des, & Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X, modl) implicit none real :: Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc real :: Ca, Cw, Ce, ty, e, b1, b2, c1, c2, inp(22) real :: P, T, DT, des(3) real :: Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X, modl(14) Tw1 = inp(1) T1 = inp(2) Tw2 = inp(3) Tc = inp(4) ec = inp(5) Tf = inp(6) Xf = inp(7) F = inp(8) Pf = inp(9) kexp = inp(10) Uf = inp(11) Ue = inp(12) Uc = inp(13) Ca = inp(14) Cw = inp(15) Ce = inp(16) ty = inp(17) e = inp(18) b1 = inp(19) b2 = inp(20) c1 = inp(21) c2 = inp(22) ! P = des(1) T = des(2) DT = des(3) ! Fp = modl(1) Fr = modl(2) Fw = modl(3) Fc = modl(4) Af = modl(5) Ae = modl(6) Ac = modl(7) Ef = modl(8) Er = modl(9) T0 = modl(10) Qf = modl(11) Qe = modl(12) Qc = modl(13) X = modl(14) return end subroutine get_data !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine direct_search(inp, des, modl) implicit none integer :: i, imax, iter, itermax, j, k, l, m, n real :: inp(22), des(3), modl(14), FUNC real :: Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc real :: Ca, Cw, Ce, ty, e, b1, b2, c1, c2 real :: P, T, DT real :: Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X real :: bak, y0, ymax, dx(3), tol, ddx(3), y, ini_dx parameter(itermax = 1000, ini_dx = 0.1, tol = 1.e-3) external FUNC ! dx = ini_dx * des ddx = tol * des open(unit=11,file="direct.out") write(*,1) 1 format(t1,'Iteration', t17, 'Function', t40, 'Variables') do1: Do iter = 1, itermax y0 = FUNC(inp, des, modl) ymax = y0 imax = 0 Do i = 1, 3 bak = des(i) des(i) = bak - dx(i) ! up y = FUNC(inp, des, modl) if(y > ymax) then ymax = y imax = 2*i-1 end if des(i) = bak + dx(i) ! down y = FUNC(inp, des, modl) if(y > ymax) then ymax = y imax = 2*i end if des(i) = bak ! restore original value End Do if(imax == 0) dx = dx/2 Do i=1,3 if(abs(dx(i)) <= abs(ddx(i))) exit do1 End Do if(imax > 0) then j = (imax+1)/2 des(j) = des(j) + (-1)**imax*dx(j) end if write(*,2) iter, ymax, des write(11,2) iter, ymax, des End Do do1 close(unit=11) 2 format(t1,i5,t15,f8.2,t30,3(f8.3,3x)) ! call get_data(Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc, & Ca, Cw, Ce, ty, e, b1, b2, c1, c2, inp, & P, T, DT, des, & Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X, modl) call results(P, T, DT, T0, EF, Qf, Fw, Af, Qe, Fr, Ae, Qc, Fc, Ac, Er, X, Fp, & Cw, Ce, Ca, b1, b2, c1, c2, e, ty) return end subroutine direct_search !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine steepest_descent(inp, des, modl) implicit none integer :: i, iter, itermax real :: inp(22), des(3), modl(14), FUNC real :: Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc real :: Ca, Cw, Ce, ty, e, b1, b2, c1, c2 real :: P, T, DT real :: Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X real :: ini_dx, lambda, dx, tol, bak, dFunc(3), dFuncNorm, bakDes(3), y, ynew parameter(itermax = 1000, ini_dx = 1.e+0, dx = 1.e-6, tol = 1.e-5) external FUNC ! open(unit=12,file="steepest.out") lambda = ini_dx * minval(abs(des)) write(*,1) 1 format(t1,'Iteration', t17, 'Function', t40, 'Variables') do1: Do iter = 1, itermax y = FUNC(inp, des, modl) Do i = 1, 3 bak = des(i) des(i) = bak + dx * abs(bak) dFunc(i) = (FUNC(inp, des, modl) - y) / abs(bak) / dx des(i) = bak End Do ! ! form unit gradient vector and move along dFuncNorm = sqrt(dot_product(dFunc, dFunc)) dFunc = dFunc / dFuncNorm bakDes = des des = des + lambda * dFunc ynew = FUNC(inp, des, modl) ! ! if not improved restore previous point and try smaller step if(ynew <= y) then des = bakDes lambda = lambda/2 ! termination criterion: if(lambda < tol) exit do1 else y = ynew endif write(*,2) iter, y, des write(12,2) iter, y, des End Do do1 close(unit=12) 2 format(t1,i5,t15,f8.2,t30,3(f8.3,3x)) ! call get_data(Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc, & Ca, Cw, Ce, ty, e, b1, b2, c1, c2, inp, & P, T, DT, des, & Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X, modl) call results(P, T, DT, T0, EF, Qf, Fw, Af, Qe, Fr, Ae, Qc, Fc, Ac, Er, X, Fp, & Cw, Ce, Ca, b1, b2, c1, c2, e, ty) return end subroutine steepest_descent !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function FUNC(inp, des, modl) implicit none real :: inp(22), des(3), modl(14), FUNC, objective real :: Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc real :: Ca, Cw, Ce, ty, e, b1, b2, c1, c2 real :: P, T, DT real :: Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X external objective call get_data(Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc, & Ca, Cw, Ce, ty, e, b1, b2, c1, c2, inp, & P, T, DT, des, & Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X, modl) call model(Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc, & P, T, DT, & Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X) FUNC = objective(Fp, Fw, Fc, Af, Ae, Ac, Ef, Er, Ca, Cw, Ce, ty, e, b1, b2, c1, c2) call put_data(Tw1, T1, Tw2, Tc, ec, Tf, Xf, F, Pf, kexp, Uf, Ue, Uc, & Ca, Cw, Ce, ty, e, b1, b2, c1, c2, inp, & P, T, DT, des, & Fp, Fr, Fw, Fc, Af, Ae, Ac, Ef, Er, T0, Qf, Qe, Qc, X, modl) return end function FUNC !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine contour(const, eps, x1l, x1u, n1, x2l, x2u, n2, x1, x2, nx) implicit none integer :: j1, j2, jac, itest integer :: n1, n2, nx,npt real :: const, eps, x1l, x1u, x2l, x2u, x1(nx), x2(nx) real :: dx1, dx2, xa, xb, ya, yb, x1s, x2a, x2b, x0, y0 real :: my_function external my_function ! npt=0 dx1 = (x1u-x1l)/(n1-1) dx2 = (x2u-x2l)/(n2-1) Do jac = 1,2 Do j1 = 1, n1 x1s = x1l*(2-jac)+x1u*(jac-1)-(-1)**jac*dx1*(j1-1) itest = 0 Do j2 = 1, n2-1 if(itest /= 1) then if(j2 /= 1) then x2a = x2b ya = yb else x2a = x2l*(2-jac)+x2u*(jac-1)-(-1)**jac*dx2*(j2-1) ya = my_function(x1s, x2a) - const endif x2b = x2a -(-1)**jac*dx2 yb = my_function(x1s, x2b) - const if(ya*yb <= 0.) then xa = x2a xb = x2b do while(abs(xa-xb) >= eps) x0 = (xa+xb)/2. y0 = my_function(x1s,x0)-const if(ya*y0 <= 0.) xb = x0 if(ya*y0 <= 0.) yb = y0 if(yb*y0 <= 0.) xa = x0 if(yb*y0 <= 0.) ya = y0 end do npt = npt+1 if(npt > nx) then print*, 'return before exceeding array size!' return endif itest = 1 x1(npt) = x1s x2(npt) = xa endif endif End Do End Do End Do return end subroutine contour !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function my_function(x, y) ! useful in drawing the contour plots. It*s somekind of interface between the contour subroutine ! and the objective function, to make the latter appear as dependent on two variables only. implicit none real :: x, y, my_function, FUNC real :: inp(22), des(3), modl(14), bak(3) common /arrays/ inp, des, modl external FUNC ! bak = des des(1) = 1. ! set pressure equal to 1 bar, so we restrict ourselves to 2 dimensions des(2) = x des(3) = y my_function = FUNC(inp, des, modl) des = bak ! return end function my_function