module func_findmin
    contains
    Integer FUNCTION findmin(LCC,maxnfe)
        use global, only : DP
        implicit none
        integer :: maxnfe,i,mindex
        real(DP) :: LCC(maxnfe),minLCC
        
        
        minLCC = minval(LCC)
        do i = 1,maxnfe
            if(LCC(i) == minLCC)then
                mindex = i
                goto 10
            endif
        enddo
        if(maxnfe == 0) mindex = 1
     10   findmin = mindex
        
     end function findmin
end module func_findmin

module opt_v_m
    use HyGSHP_calc    
    use HyGSHP_calc_down  
    use HyGSHP_calc_conv
    contains
    subroutine opt_v(DIRx_print,DIRf,nspeeds,opt1,opt2)
    use global
    use subplex_global
    use func_findmin
    use global_print
    implicit none 
    
    real(DP),intent(in),optional :: opt1    !indicates that only the control set points are optimized
    integer,intent(in),optional :: opt2 !doing multiple optimizations in 1, add a number to the file name so you get multiple runs/files
    character(len = 1) :: opt2c
    integer :: i,runs,k,shrink,ii,nspeeds,jj
    real(DP) :: in(10),LCC,totallength,size_tower,dt1,tc1,tc2,tpc,dtpc,dt2,size_boil
    
    integer :: niwmax,nwmax,nxmax
    parameter (niwmax=300,nwmax=300,nxmax=300)
    integer :: iflag,iwork(niwmax),maxnfe,mdcont,mdsing,mduser,mode, &
               nf1,nf2,nfe,nfinc,n,nfe_tot         
    real(DP) :: fx,scales(nxmax),scl, tol,tol1,tol2,tolfac,work(nwmax), &
                x(nwmax),junk,lower(nxmax),upper(nxmax),lower_orig(nxmax),upper_orig(nxmax)
    real(DP) :: lower_in(nxmax),upper_in(nxmax),x_in(nwmax),scales_in(nxmax), &
                x_in_conv(nwmax),x_in_hybrid(nwmax)           
    real(DP) :: tol1_main,tol1_starter,dfmin_lim_main,dfmin_lim_starter,dfmin_lim
    integer :: count_lim_main,count_lim_starter,count_lim
    
    character*26 flgdsc(-2:2)
    external subplx
    Integer :: mindex,mindex2 
    data flgdsc /'invalid input', &
                  'maxnfe exceeded', &
                  'tol satisfied', &
                  'limit of machine precision', &
                  'fstop reached'/
    Integer :: DIRmaxf
    real(DP) :: DIRx(128),opt_time,DIRf_main,DIRf_starter,DIRf_min,DIRx_hold(128), &
                opt_time_tot,DIRx_print_hold(10,5),DIRf_hold(5)
    real(DP), intent(out) :: DIRx_print(10),DIRf
    integer :: t1,t2,clock_rate,clock_max      
    integer :: yr_orig,m_or_s,s,r_furnace,opt_repeat
    real(DP) :: GHX_next, CCCT_next, ratio_opt 
    

    namelist /optimization_settings_1/ scales_in,x_in_conv,lower_in,upper_in,x_in_hybrid
    namelist /optimization_settings/ tol2,tolfac,nf1,nf2,nfinc,tol1_main,tol1_starter, &
                                     count_lim_main,count_lim_starter,dfmin_lim_main,dfmin_lim_starter
    
    delt_orig = delt
    yr_orig = yr
    xkeep = 0
    x = 0
    lower = 0
    upper = 0
    scales = 0
    xpass = 0
    opt_repeat = 0
    nfe_tot = 0
    opt_time_tot = 0
    DIRx_print_hold = 0
    DIRf_hold = 0
    jj = 0
                
    if(.not.present(opt1))then      
        DIRx_print = 0
    endif   
    
    if(present(opt2))then
        call file_names(opt2)
    else
        call file_names
    endif    
    
    open(unit = 21, file = trim(file1),recl = 2140000)
    if(conv == 0)then
        write(21,15) "totallength","size_towerkw","Qset_ratio","ct1","ct2", &
                     "ghx1","tpc","dtpc","LCC","Imbalance_pen", &
                     "energy_pen"
        write(21,15) "[m]","[kW]","[-]","[C]","[C]", &
                     "[C]","[C]","[C]","[$]","[$]", &
                     "[$]"
        15 format(a,',',a,',',a,',',a,',',a,',', &
                  a,',',a,',',a,',',a,',',a,',', &
                  a) 
                  
!        write(21,15) "totallength","size_towerkw","ct1","ct2","ghx1", &
!                     "tpc","dtpc","LCC","exit_early","PChours", &
!                     "Imbalance_pen","Qset_ratio","Qset","ghx_load_tot","ghx_ratio", &
!                     "energy_pen","m","b"
!        write(21,15) "[m]","[kW]","[C]","[C]","[C]", &
!                     "[C]","[C]","[$]","[]","[hrs]", &
!                     "[$]","[ratio]","[kJ/hr]","[C]","[]", &
!                     "[$]","[]","[]"
!        15 format(a,',',a,',',a,',',a,',',a,',', &
!                  a,',',a,',',a,',',a,',',a,',', &
!                  a,',',a,',',a,',',a,',',a,',', &
!                  a,',',a,',',a)                   
        
!        write(*,15) "totallength","size_towerkw","ct1","ct2","ghx1","tpc","dtpc","LCC"
!        write(*,16) "[m]","[kW]","[C]","[C]","[C]","[C]","[C]","[$]" 
!        16 format(2x,a,13x,a,13x,a,5x,a,5x,a,5x,a,5x,a,6x,a)   
    else
        write(21,17) "size_towerkw","Furnace","ct1","ct2","LCC"
        write(21,17) "[kW]","[kBH]","[C]","[C]","[$]"
        17 format(a,',',a,',',a,',',a,',',a)     
        
!        write(*,17) "size_towerkw","size_furnace","ct1","ct2","LCC"
!        write(*,18) "[kW]","[kBH]","[C]","[C]","[$]"
!        18 format(2x,a,13x,a,13x,a,5x,a,5x,a)   
    endif
    
    
    open(unit = 1,file = "..\..\namelist_files\optim_settings.txt")
    read(1,nml = optimization_settings_1) 
    close(1)

    open(unit = 1,file = "..\..\namelist_files\optim_settings.txt")
    read(1,nml = optimization_settings) 
    close(1)

    if(conv == 1)then
        x_in = x_in_conv
    else
        x_in = x_in_hybrid
    endif

    n = var_vect(1)
    if(present(opt1)) x_in(1:8) = DIRx_print(1:8)
    if(furnace)then
        !x_in(8) = 0.9
        lower_in(8) = 0
        upper_in(8) = 1
        scales_in(8) = 0.1
    endif
    
    if(nspeeds > 2)then
        lower_in(9) = 0
        upper_in(9) = 5
        lower_in(10) = -5
        upper_in(10) = 5
        scales_in(9) = 0.2
        scales_in(10) = 0.5
        x(9) = 0.1
        x(10) = -0.5        
    endif
    
    do ii = 2,n+1
        x(ii-1) = x_in(var_vect(ii))   
        lower(ii-1) = lower_in(var_vect(ii))    
        upper(ii-1) = upper_in(var_vect(ii))
        scales(ii-1) = scales_in(var_vect(ii))     
        var_vect_pl(ii-1) = var_vect_text(var_vect(ii)) 
    enddo
    if(furnace == 1)then
        do ii = 2,n+1
            if(var_vect(ii) == 8)then
                r_furnace = ii-1                    
            endif
        enddo
    endif
                
    if(conv)then      
        C_CCCT = x(1)        
    endif
    
    lower_orig = lower
    upper_orig = upper

    maxnfe = nf1
    
!    write(*,*) 3
!    read(*,*) junk
      
    ! First call to subplx so continuation mode is off.
    mdcont = 0
    ! Using default options so user options mode is off.
    mduser = 0
    ! Using optimization so single-step mode is off.
    mdsing = 0
    
    20   continue
    tol1 = tol1_main
    tol = tol1
    tol2 = tol1
    tolfac = tol1
    if(strat > 8) opt_method = 1
        
    mode = 4*mdsing + 2*mduser + mdcont   
    if((opt_method == 2).and.(yr > 0))then 
        s = 1   
        double = 0            
        yr = yr_opt
        first_run_opt = .false.
        delt = 1
        qcool(1:8760*yr_opt + 1) = qcool_opt
        qheat(1:8760*yr_opt + 1) = qheat_opt
        tdb(1:8760*yr_opt + 1) = tdb_opt
        twb(1:8760*yr_opt + 1) = twb_opt
        
        tol1 = tol1_starter
        count_lim = count_lim_starter
        dfmin_lim = dfmin_lim_starter
        tol = tol1
        tol2 = tol1
        tolfac = tol1

        if(opt_starter == 2)then    !Chose the direct method to do the 1yr optimization
            print *,'***********************************'
            print *,'****** DIRECT minimization ********'
            print *,'***********************************'
            if(present(opt1))then  
                xpass(sens_var) = opt1  
                L_hyb_m = x_in(1)
                C_CCCT = x_in(2)
                Qset_ash = x_in(8)
            endif
            select_file = 1
            call main_direct(L_hyb_m,C_CCCT,DIRx,DIRf,DIRmaxf,opt_time,0,s,lower,upper,x_in) 
            DIRf_starter = DIRf            
        elseif(opt_starter == 1)then !Chose the subplex method to do the 1yr optimization          
            print *,'************************************'
            print *,'****** SUBPLEX minimization ********'
            print *,'************************************'
    
            if(present(opt1))then
                xpass(sens_var) = opt1  
            else
                x(1) = L_hyb_m
                x(2) = C_CCCT
                x(8) = Qset_ash
            endif
            select_file = 1
            if(conv == 1)then
                call subplx(HyGSHP_sub_conv,n,tol,maxnfe,mode,scales,x,fx,nfe, &
                            work,iwork,iflag,lower,upper,opt_time,count_lim,dfmin_lim)  
            elseif(downstream == 0)then
                call subplx(HyGSHP_sub,n,tol,maxnfe,mode,scales,x,fx,nfe, &
                            work,iwork,iflag,lower,upper,opt_time,count_lim,dfmin_lim)  
            else
                call subplx(HyGSHP_sub_down,n,tol,maxnfe,mode,scales,x,fx,nfe, &
                            work,iwork,iflag,lower,upper,opt_time,count_lim,dfmin_lim)
            endif 
            
            mindex = findmin(xkeep(:,n+1),nfe)   
            DIRx(1:n) = xkeep(mindex,1:n)  
            DIRf_starter = xkeep(mindex,n+1)
            
            do ii = 2,n+1
                DIRx_print(var_vect(ii)) = DIRx(ii - 1)
            enddo
            DIRx_print(9) = xkeep(mindex,n+2)
            DIRx_print(10) = xkeep(mindex,n+3)
            
            if(present(opt1)) DIRx_print(sens_var) = opt1
            
            do i = 1,n
                write(subplex_log,230) i,DIRx(i)
            enddo
            
            write(subplex_log,200) xkeep(mindex,n+1)
            write(subplex_log,210) nfe
            write(subplex_log,220) opt_time      
            
            DIRmaxf = nfe           
        endif
        first_run_opt = .false.
        
        write(*,*)
        write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
        write(*,*) '!Initial Optimization Complete, Starting Main Optimization!'
        write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
        write(*,*)
        
        write(21,*) '!Initial Optimization Complete, Starting Main Optimization!'
        
        DIRx_hold = DIRx
        double = 1
        yr = yr_orig
        delt = delt_orig
        qcool = qcool_orig
        qheat = qheat_orig
        tdb = tdb_orig
        twb = twb_orig
        
        first_run_opt = .false.
        tol1 = tol1_main
        count_lim = count_lim_main
        dfmin_lim = dfmin_lim_main
        tol = tol1
        tol2 = tol1
        tolfac = tol1   
        lower = lower_orig
        upper = upper_orig
        
        do ii = 2,n+1
            if(var_vect(ii) <= 2)then
                lower(var_vect(ii)) = DIRx(var_vect(ii))*(0.5)
                upper(var_vect(ii)) = DIRx(var_vect(ii))*(4)
            endif
        enddo
        
        s = 0
500     jj = jj + 1                         
        if(opt_main == 2)then  !Chose the direct method to do the main optimization
            print *,'********************************************'
            print *,'******  DIRECT minimization of fun  ********'
            print *,'********************************************'
            GHX_next = DIRx(1)
            CCCT_next = DIRx(2)  
            if(present(opt1)) xpass(sens_var) = opt1  
            select_file = 2
            call main_direct(GHX_next,CCCT_next,DIRx,DIRf,DIRmaxf,opt_time,1,s,lower,upper,x_in)
            do ii = 2,n+1
                DIRx_print(var_vect(ii)) = DIRx(ii - 1)
            enddo
            DIRx_print(9) = DIRx(n+1)
            DIRx_print(10) = DIRx(n+2)
            
            if(strat >= 9)then
                DIRx_print(1) = xpass(1)
                DIRx_print(2) = xpass(2)
            endif
            if(present(opt1)) DIRx_print(sens_var) = opt1 
            Dirf_main = dirf
        elseif(opt_main == 1)then !Chose the subplex method to do the main optimization    
            print *,'********************************************'
            print *,'******  subplx minimization of fun  ********'
            print *,'********************************************'
                        
            x(1:n) = DIRx(1:n)   
            if(present(opt1)) xpass(sens_var) = opt1   
            select_file = 2
            if(conv == 1)then
                call subplx(HyGSHP_sub_conv,n,tol,maxnfe,mode,scales,x,fx,nfe, &
                            work,iwork,iflag,lower,upper,opt_time,count_lim,dfmin_lim)  
            elseif(downstream == 0)then
                call subplx(HyGSHP_sub,n,tol,maxnfe,mode,scales,x,fx,nfe, &
                            work,iwork,iflag,lower,upper,opt_time,count_lim,dfmin_lim)  
            else
                call subplx(HyGSHP_sub_down,n,tol,maxnfe,mode,scales,x,fx,nfe, &
                            work,iwork,iflag,lower,upper,opt_time,count_lim,dfmin_lim)
            endif 
            mindex = findmin(xkeep(:,n+1),nfe)   
            DIRx(1:n) = xkeep(mindex,1:n)  
            do ii = 2,n+1
                DIRx_print(var_vect(ii)) = DIRx(ii - 1)
            enddo
            DIRx_print(9) = xkeep(mindex,n+2)
            DIRx_print(10) = xkeep(mindex,n+3)
            
            if(strat >= 9)then
                DIRx_print(1) = xpass(1)
                DIRx_print(2) = xpass(2)
            endif
            if(present(opt1)) DIRx_print(sens_var) = opt1 
            !DIRf = xkeep(mindex,n+1)
            DIRf_main = xkeep(mindex,n+1)
            DIRmaxf = nfe     
            
            do i = 1,n
                write(subplex_log,230) i,DIRx(i)
            enddo
            
            write(subplex_log,200) xkeep(mindex,n+1)
            write(subplex_log,210) nfe
            write(subplex_log,220) opt_time     
            
            if(PCon_sum(mindex,1) == 0)then
                write(subplex_log,*) 'There is no PC'
            else
                write(subplex_log,*) 'There is PC'
            endif             
        endif  
        first_run_opt = .false.  
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        dirx_print_hold(:,jj) = DIRx_print
        dirf_hold(jj) = dirf_main        
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ratio_opt = dirf_main/dirf_starter
        if( ((ratio_opt > 5).and.(opt_repeat < 3)).or. &
            ((DIRmaxf <= 10).and.(ratio_opt > 2).and.(opt_repeat < 3)).or. &
            ((dirf_main >= 1e8).and.(opt_repeat < 3)))then
            opt_repeat = opt_repeat + 1
            nfe_tot = nfe_tot + DIRmaxf
            opt_time_tot = opt_time_tot + opt_time
            do ii = 2,n+1
                if((var_vect(ii) > 2).and.(var_vect(ii) /= 8))then
                    DIRx(ii-1) = DIRx_hold(ii-1)-(1*opt_repeat)
                elseif(var_vect(ii) == 8)then
                    DIRx(ii-1) = (1-opt_repeat/10)*DIRx_hold(ii-1)
                elseif(var_vect(ii) == 1)then
                    DIRx(ii-1) = DIRx_hold(ii-1) - 400*opt_repeat 
                elseif(var_vect(ii) == 2)then
                    DIRx(ii-1) = DIRx_hold(ii-1) + 100*opt_repeat    
                endif
            enddo
            write(21,*) "First optimization attempt did not produce a good result, try again"   
            if(opt_repeat == 3)then
                write(*,*) '***********************************************************'
                write(*,*) 'Notice: Main optimization had difficulty converging to an optimum'
                write(*,*) 'Look carefully at the resulting design'
                write(*,*) '***********************************************************'
                write(21,*) '***********************************************************'
                write(21,*) 'Notice: Main optimization had difficulty converging to an optimum'
                write(21,*) 'Look carefully at the resulting design'
                write(21,*) '***********************************************************'           
            endif          
            goto 500                 
        endif
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        if(opt_repeat /= 0)then
            DIRmaxf = DIRmaxf + nfe_tot
            opt_time = opt_time + opt_time_tot
            mindex2 = findmin(dirf_hold(1:jj),jj)
            DIRx_print = Dirx_print_hold(:,mindex2)
            DIRf_main = DIRf_hold(mindex2)          
        endif
        
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       
                 
    elseif(opt_method == 1)then
        s = 0
        tol1 = tol1_main
        count_lim = count_lim_main
        dfmin_lim = dfmin_lim_main
        tol = tol1
        tol2 = tol1
        tolfac = tol1
502     jj = jj + 1 
        if(opt_main == 2)then 
            if(guess == 1)then                   
50              if(strat >= 9)then
                    xpass(1) = x_in(1)
                    xpass(2) = x_in(2) 
                    if(furnace == 1) xpass(r_furnace) = x_in(8)  
                elseif(strat < 9)then 
                    x(1) = x_in(1)
                    x(2) = x_in(2)   
                    if(furnace == 1) x(r_furnace) = x_in(8)                 
                endif
                if(present(opt1)) xpass(sens_var) = opt1   
                select_file = 2   
                if(conv)then                               
                    call main_direct(x(1),x(1),DIRx,DIRf,DIRmaxf,opt_time,1,s,lower,upper,x_in)
                else
                    call main_direct(x(1),x(2),DIRx,DIRf,DIRmaxf,opt_time,1,s,lower,upper,x_in)
                endif
            else
                if(strat >= 9)then
                    xpass(1) = L_hyb_m
                    xpass(2) = C_CCCT     
                    if(furnace == 1) xpass(r_furnace) = Qset_ash               
                endif  
                if(present(opt1)) xpass(sens_var) = opt1  
                select_file = 2
                call main_direct(L_hyb_m,C_CCCT,DIRx,DIRf,DIRmaxf,opt_time,1,s,lower,upper,x_in)
            endif
            do ii = 2,n+1
                DIRx_print(var_vect(ii)) = DIRx(ii - 1)
            enddo
            DIRx_print(9) = DIRx(n+1)
            DIRx_print(10) = DIRx(n+2)
            
            if(strat >= 9)then
                DIRx_print(1) = xpass(1)
                DIRx_print(2) = xpass(2)
            endif
            if(present(opt1)) DIRx_print(sens_var) = opt1 
            !DIRf = xkeep(mindex,n+1)
             
        elseif(opt_main == 1)then
            if((guess /= 1).and.(strat < 9))then
                x(1) = L_hyb_m
                x(2) = C_CCCT
                if(furnace == 1) x(r_furnace) = Qset_ash
            elseif((guess == 1).and.(strat < 9))then
                x(1) = x_in(1)
                x(2) = x_in(2)   
                if(furnace == 1) x(r_furnace) = x_in(8) 
            endif
            if((strat >= 9).and.(guess /= 1))then
                xpass(1) = L_hyb_m
                xpass(2) = C_CCCT
                if(furnace == 1) xpass(r_furnace) = Qset_ash
            elseif((strat >= 9).and.(guess == 1))then
                xpass(1) = x_in(1)
                xpass(2) = x_in(2)   
                if(furnace == 1) xpass(r_furnace) = x_in(8)               
            endif  
            if(present(opt1)) xpass(sens_var) = opt1  
            select_file = 2
            if(conv == 1)then
                call subplx(HyGSHP_sub_conv,n,tol,maxnfe,mode,scales,x,fx,nfe, &
                            work,iwork,iflag,lower,upper,opt_time,count_lim,dfmin_lim)  
            elseif(downstream == 0)then
                call subplx(HyGSHP_sub,n,tol,maxnfe,mode,scales,x,fx,nfe, &
                            work,iwork,iflag,lower,upper,opt_time,count_lim,dfmin_lim)  
            else
                call subplx(HyGSHP_sub_down,n,tol,maxnfe,mode,scales,x,fx,nfe, &
                            work,iwork,iflag,lower,upper,opt_time,count_lim,dfmin_lim)
            endif 
            mindex = findmin(xkeep(:,n+1),nfe)   
            DIRx(1:n) = xkeep(mindex,1:n)   
            do ii = 2,n+1
                DIRx_print(var_vect(ii)) = DIRx(ii - 1)
            enddo
            DIRx_print(9) = xkeep(mindex,n+2)
            DIRx_print(10) = xkeep(mindex,n+3)
            
            if(strat >= 9)then
                DIRx_print(1) = xpass(1)
                DIRx_print(2) = xpass(2)
            endif
            if(present(opt1)) DIRx_print(sens_var) = opt1
            DIRf = xkeep(mindex,n+1)
            DIRmaxf = nfe    
            
            do i = 1,n
                write(subplex_log,230) i,DIRx(i)
            enddo
            
            write(subplex_log,200) xkeep(mindex,n+1)
            write(subplex_log,210) nfe
            write(subplex_log,220) opt_time     
            
            if(PCon_sum(mindex,1) == 0)then
                write(subplex_log,*) 'There is no PC'
            else
                write(subplex_log,*) 'There is PC'
            endif       
        endif    
        dirf_main = dirf
        
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        dirx_print_hold(:,jj) = DIRx_print
        dirf_hold(jj) = dirf_main        
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        if( ((DIRmaxf <= 10).and.(opt_repeat < 3)).or. &
            ((dirf_main >= 1e8).and.(opt_repeat < 3))    )then
            opt_repeat = opt_repeat + 1
            nfe_tot = nfe_tot + DIRmaxf
            opt_time_tot = opt_time_tot + opt_time
            do ii = 2,n+1
                if((var_vect(ii) > 2).and.(var_vect(ii) /= 8))then
                    DIRx(ii-1) = DIRx_hold(ii-1)-(1*opt_repeat)
                elseif(var_vect(ii) == 8)then
                    DIRx(ii-1) = (1-opt_repeat/10)*DIRx_hold(ii-1)
                elseif(var_vect(ii) == 1)then
                    DIRx(ii-1) = DIRx_hold(ii-1) - 400*opt_repeat 
                elseif(var_vect(ii) == 2)then
                    DIRx(ii-1) = DIRx_hold(ii-1) + 100*opt_repeat    
                endif
            enddo
            write(21,*) "First optimization attempt did not produce a good result, try again"        
            if(opt_repeat == 3)then
                write(*,*) '***********************************************************'
                write(*,*) 'Notice: Main optimization had difficulty converging to an optimum'
                write(*,*) 'Look carefully at the resulting design'
                write(*,*) '***********************************************************'
                write(21,*) '***********************************************************'
                write(21,*) 'Notice: Main optimization had difficulty converging to an optimum'
                write(21,*) 'Look carefully at the resulting design'
                write(21,*) '***********************************************************'           
            endif    
            goto 502                 
        endif
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        if(opt_repeat /= 0)then
            DIRmaxf = DIRmaxf + nfe_tot
            opt_time = opt_time + opt_time_tot
            mindex2 = findmin(dirf_hold(1:jj),jj)
            DIRx_print = Dirx_print_hold(:,mindex2)
            DIRf_main = DIRf_hold(mindex2)          
        endif        
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
        
        
    elseif(opt_method == 3)then
    !Run the DIRECT method until the change in the optimal value is less than the specified
    !amount, then run the SUBPLEX method
        s = 0
        tol1 = tol1_main
        count_lim = count_lim_main
        dfmin_lim = dfmin_lim_main
        tol = tol1
        tol2 = tol1
        tolfac = tol1
        print *,'********************************************'
        print *,'******  DIRECT minimization of fun  ********'
        print *,'********************************************'
        if(present(opt1))then
            xpass(sens_var) = opt1  
            L_hyb_m = x_in(1)
            C_CCCT = x_in(2)
            qset_ash = x_in(8)
        endif
        select_file = 1
        call main_direct(L_hyb_m,C_CCCT,DIRx,DIRf,DIRmaxf,opt_time,0,s,lower,upper,x_in)  
        dirf_starter = dirf
        first_run_opt = .false.
        
        write(*,*)
        write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
        write(*,*) '!Starting refined optimization                            !'
        write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
        write(*,*)
        
        write(21,*) '!Initial Optimization Complete, Starting Final Optimization!'
        
503     jj = jj + 1 
        print *,'********************************************'
        print *,'******  subplx minimization of fun  ********'
        print *,'********************************************'
        
        lower = lower_orig
        upper = upper_orig
        
        do ii = 2,n+1
            if(var_vect(ii) <= 2)then
                lower(var_vect(ii)) = DIRx(var_vect(ii))*(0.5)
                upper(var_vect(ii)) = DIRx(var_vect(ii))*(4)
            endif
        enddo
                    
501     x(1:n) = DIRx(1:n)   
        if(present(opt1)) xpass(sens_var) = opt1   
        select_file = 2
        if(conv == 1)then
            call subplx(HyGSHP_sub_conv,n,tol,maxnfe,mode,scales,x,fx,nfe, &
                        work,iwork,iflag,lower,upper,opt_time,count_lim,dfmin_lim)  
        elseif(downstream == 0)then
            call subplx(HyGSHP_sub,n,tol,maxnfe,mode,scales,x,fx,nfe, &
                        work,iwork,iflag,lower,upper,opt_time,count_lim,dfmin_lim)  
        else
            call subplx(HyGSHP_sub_down,n,tol,maxnfe,mode,scales,x,fx,nfe, &
                        work,iwork,iflag,lower,upper,opt_time,count_lim,dfmin_lim)
        endif 
        mindex = findmin(xkeep(:,n+1),nfe)   
        DIRx(1:n) = xkeep(mindex,1:n)  
        do ii = 2,n+1
            DIRx_print(var_vect(ii)) = DIRx(ii - 1)
        enddo
        DIRx_print(9) = xkeep(mindex,n+2)
        DIRx_print(10) = xkeep(mindex,n+3)
        
        if(strat >= 9)then
            DIRx_print(1) = xpass(1)
            DIRx_print(2) = xpass(2)
        endif
        if(present(opt1)) DIRx_print(sens_var) = opt1 
        DIRf_main = xkeep(mindex,n+1)
        DIRmaxf = nfe   
        
        do i = 1,n
            write(subplex_log,230) i,DIRx(i)
        enddo
        
        write(subplex_log,200) xkeep(mindex,n+1)
        write(subplex_log,210) nfe
        write(subplex_log,220) opt_time     
        
        if(PCon_sum(mindex,1) == 0)then
            write(subplex_log,*) 'There is no PC'
        else
            write(subplex_log,*) 'There is PC'
        endif    
        
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        dirx_print_hold(:,jj) = DIRx_print
        dirf_hold(jj) = dirf_main        
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ratio_opt = dirf_main/dirf_starter
        if( ((ratio_opt > 5).and.(opt_repeat < 3)).or. &
            ((DIRmaxf <= 10).and.(ratio_opt > 2).and.(opt_repeat < 3)).or. &
            ((dirf_main >= 1e8).and.(opt_repeat < 3)))then
            opt_repeat = opt_repeat + 1
            nfe_tot = nfe_tot + DIRmaxf
            opt_time_tot = opt_time_tot + opt_time
            do ii = 2,n+1
                if((var_vect(ii) > 2).and.(var_vect(ii) /= 8))then
                    DIRx(ii-1) = DIRx_hold(ii-1)-(1*opt_repeat)
                elseif(var_vect(ii) == 8)then
                    DIRx(ii-1) = (1-opt_repeat/10)*DIRx_hold(ii-1)
                elseif(var_vect(ii) == 1)then
                    DIRx(ii-1) = DIRx_hold(ii-1) - 400*opt_repeat 
                elseif(var_vect(ii) == 2)then
                    DIRx(ii-1) = DIRx_hold(ii-1) + 100*opt_repeat    
                endif
            enddo
            write(21,*) "First optimization attempt did not produce a good result, try again"        
            if(opt_repeat == 3)then
                write(*,*) '***********************************************************'
                write(*,*) 'Notice: Main optimization had difficulty converging to an optimum'
                write(*,*) 'Look carefully at the resulting design'
                write(*,*) '***********************************************************'
                write(21,*) '***********************************************************'
                write(21,*) 'Notice: Main optimization had difficulty converging to an optimum'
                write(21,*) 'Look carefully at the resulting design'
                write(21,*) '***********************************************************'           
            endif    
            goto 503                 
        endif
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        if(opt_repeat /= 0)then
            DIRmaxf = DIRmaxf + nfe_tot
            opt_time = opt_time + opt_time_tot
            mindex2 = findmin(dirf_hold(1:jj),jj)
            DIRx_print = Dirx_print_hold(:,mindex2)
            DIRf_main = DIRf_hold(mindex2)          
        endif        
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!                        
    endif      
    
    first_run_opt = .false.
    
       write (6,1020) maxnfe,tol,fx,DIRmaxf,iflag
    1020   format (t5,i6,t15,e10.2,t30,e13.5,t45, &
           i6,t60,i5)
           
    if (iflag .eq. -1) then
    if (maxnfe .ge. nf2) go to 30
        maxnfe = maxnfe+nfinc
    else if (iflag .eq. 0) then
        if (tol .le. tol2) go to 30
            tol = tol*tolfac
    else
        go to 30
    end if

    mdcont = 1

    go to 20

    30   continue
    print *, '****************************************'
    print *, '******* optimization results ***********'
    print *, '****************************************'
    1025   format (' iflag =',i3,5x,'(',a26,')')
    write (6,1030) DIRmaxf,fx
    1030   format (/1x,'nfe=',i7,5x,'fx=',e16.8/)
    print *, 'x ='
    write (6,1005) (x(i),i=1,n)
    1005 format (t3,5e15.5)
    1010   format (///t5,'maxnfe',t15,'tol',t30, &
           'fx',t45,'nfe',t60,'iflag'/)
    !+-----------------------------------------------------------------------+
    !| Print optimization results.                                           |
    !+-----------------------------------------------------------------------+
    open(unit = 23,file = trim(file2),recl = 2140000)
    write(23,100)
    write(23,120)
    write(23,130) DIRx_print(1)
    write(23,140) DIRx_print(2)
    write(23,250) (1-DIRx_print(8))*maxval(qheat)*0.0009478*abs((furnace+conv > 0))
    write(23,240) DIRx_print(8)    
    write(23,150) DIRx_print(3) 
    write(23,160) DIRx_print(4)
    write(23,170) DIRx_print(5)
    write(23,180) DIRx_print(6)
    write(23,190) DIRx_print(7)   
    write(23,200) DIRf_main
    write(23,210) DIRmaxf
    write(23,220) opt_time
    write(23,270) DIRx_print(9) 
    write(23,280) DIRx_print(10) 
    100   FORMAT('-------------- Final result ------------------')
    120   FORMAT('Design of optimal system    : ')
    130   FORMAT('      Length of GHX (m)    : ',f20.6)
    140   FORMAT('      Size of CT (kW)    : ',f12.7)
    250   FORMAT('      Size of Furnace (kBTU/h)    : ',f12.7)
    150   FORMAT('      CT1 (°C)    : ',f12.7)
    160   FORMAT('      CT2 (°C)    : ',f12.7)
    170   FORMAT('      GHX1 (°C)    : ',f12.7)
    180   FORMAT('      TPC (°C)    : ',f12.7)
    190   FORMAT('      DTPC (°C)    : ',f12.7)
    240   FORMAT('      Qset_ratio      : ',f12.7)
    200   FORMAT('Minimum LCC ($)   : ',f12.2)
    270   FORMAT('First Cost ($)   : ',f12.2)
    280   FORMAT('Total Fuel Cost ($)   : ',f15.2)
    210   FORMAT('Number of function-evaluations : ',I5)
    220   FORMAT('Time needed             : ',f20.6,' seconds.')    
    230   FORMAT('Final value of x',i2,': ',f20.6)   
    close(unit = 23)    

    close(unit = 21)    
    close(subplex_log)

end subroutine opt_v

end module opt_v_m

