        MODULE functions    
      
        USE constants, ONLY : dp   ! double precision (i.e. accuracy)    
        IMPLICIT NONE
        
      
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !| |                                                                          | |
        !| |                                                                          | |
        !| |                 INFORMATION SUPPLIED BY THE USER:                        | | 
        !| |                                                                          | |
        !| |                                                                          | |
        !| |    * PROBLEM specification:                                              | |
        !| |                                                                          | |
        !| |        - the number of objectives:    'number_of_obj'                    | |       
        !| |        - the number of constraints:   'number_of_const'                  | |
        !| |        - the number of variables:     'user_n'                           | |       
        !| |                                                                          | |       
        !| |        - The DC components f1 and g1: 'problem1'                         | |
        !| |        - The DC components f2 and g2: 'problem2'                         | |       
        !| |        - The starting point:          'startpoint'                       | |       
        !| |                                                                          | |
        !| |    * Different PARAMETERS:                                               | |
        !| |                                                                          | |           
        !| |        - the number of variables:     'user_n'                           | |       
        !| |        - the stopping tolerance:      'user_crit_tol'                    | |
        !| |        - the extra descent parameter  'user_m_extra'                     | |               
        !| |                                                                          | |       
        !| |        MAIN ITERATION:                                                   | |  
        !| |                                                                          | |               
        !| |        - the size of bundle B_1:      'user_size_b1'                     | |       
        !| |        - the size of bundle B_2:      'user_size_b2'                     | |       
        !| |        - the descent parameter:       'user_m'                           | |       
        !| |        - the decrease parameter:      'user_r_dec'                       | |       
        !| |        - the increase parameter:      'user_r_inc'                       | |           
        !| |        - the small number:            'user_eps'                         | |       
        !| |        - the decrease parameter:      'user_c1'                          | |       
        !| |        - the decrease parameter:      'user_c2'                          | |       
        !| |        - the decrease parameter:      'user_c3'                          | |       
        !| |                                                                          | |       
        !| |        CLARKE STATIONARY ALGORITHM:                                      | | 
        !| |                                                                          | |               
        !| |        - the size of bundle B:        'user_size'                        | |       
        !| |        - the descent parameter:       'user_m_clarke'                    | |       
        !| |        - the step-length tolerance:   'user_step_tol'                    | |            
        !| |                                                                          | |       
        !| |                       .**..**..**..**..**..**.                           | |       
        !| |                              OBJECTIVES:                                 | |
        !| |                       .**..**..**..**..**..**.                           | |               
        !| |                                                                          | |       
        !| |    * Computation of the values of the DC functions f_1 and f_2:          | |
        !| |        - f1(y, problem)    the value of DC component f_1 at a point y    | |
        !| |        - f2(y, problem)    the value of DC component f_2 at a point y    | |           
        !| |                                                                          | |               
        !| |                                                                          | |               
        !| |    * Computation of the subgradient of the DC components f_1 and f_2:    | |
        !| |        - subgradient_f1(y, problem)    the subgradient of f_1 at y       | |
        !| |        - subgradient_f2(y, problem)    the subgradient of f_2 at y       | |       
        !| |                                                                          | |
        !| |                                                                          | |
        !| |                       .**..**..**..**..**..**.                           | |       
        !| |                              CONSTRAINTS:                                | |
        !| |                       .**..**..**..**..**..**.                           | |               
        !| |                                                                          | |       
        !| |    * Computation of the values of the DC functions g_1 and g_2:          | |
        !| |        - g1(y, problem)    the value of DC component g_1 at a point y    | |
        !| |        - g2(y, problem)    the value of DC component g_2 at a point y    | |           
        !| |                                                                          | |               
        !| |                                                                          | |               
        !| |    * Computation of the subgradient of the DC components g_1 and g_2:    | |
        !| |        - subgradient_g1(y, problem1)    the subgradient of g_1 at y      | |
        !| |        - subgradient_g2(y, problem2)    the subgradient of g_2 at y      | |       
        !| |                                                                          | |
        !| |                                                                          | |       
        !| |                                                                          | |
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*
        
        
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !| |                                                                                              | |
        !| |                                                                                              | |
        !| |                           OTHER FUNCTIONS:                                                   | |
        !| |                                                                                              | |
        !| |    The value of H_1:           H1(AandB_y_all)                                               | |
        !| |    The value of H_2:           H2(fg2_y_all)                                                 | |
        !| |                                                                                              | |
        !| |    The subgradient of H_1:     subgradient_H1(grad1_y_all, grad2_y_all, AandB_y_all)         | |
        !| |    The subgradient of H_2:     subgradient_H2(grad2_y_all)                                   | |
        !| |                                                                                              | |
        !| |    Index of max comp. in H_1:  index_H1(AandB_y_all)                                         | |
        !| |                                                                                              | |
        !| |    The value of A_i or B_l:         AorB_i(fg1_y_all, fg2_y_all, fg1_k_all, fg2_k_all, ind)  | |
        !| |    The subgradient of A_i of B_l:   subgradient_AorB_i(grad1_y_all, grad2_y_all, ind)        | |
        !| |                                                                                              | |
        !| |                                                                                              | |
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*       
 

        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !| |                                                                                                              | |
        !| |                                                                                                              | |
        !| |                                             TEST PROBLEMS:                                                   | |
        !| |                                                                                                              | |
        !| |    * UNCONSTRAINED PROBLEMS:                                                                                 | |
        !| |                                                                                                              | |
        !| |                     objective functions       starting point        dimension                                | |
        !| |                   -----------------------    -----------------    -------------                              | |
        !| |    Problem 1:             2, 6                      1                   2                                    | |
        !| |    Problem 2:             2, 7                      2                   2                                    | |
        !| |    Problem 3:             2, 7                      3                   2                                    | |
        !| |    Problem 4:             6, 7                      4                   2                                    | |
        !| |    Problem 5:             3, 9                      5                   4                                    | |
        !| |    Problem 6:             10, 13                    6                   10                                   | |
        !| |    Problem 7:             12, 13                    7                   10                                   | |
        !| |    Problem 8:             4, 10                     8                 varies                                 | |
        !| |    Problem 9:             10, 12                    9                 varies                                 | |
        !| |    Problem 10:            14, 15                    10                varies                                 | |
        !| |                                                                                                              | |
        !| |    Problem 11:            2, 6, 7                   11                  2                                    | |
        !| |    Problem 12:            3, 4, 9                   12                  4                                    | |
        !| |    Problem 13:            4, 10, 12                 13                varies                                 | |
        !| |    Problem 14:            4, 10, 16                 14                varies                                 | |
        !| |    Problem 15:            10, 14, 15                15                varies                                 | |
        !| |                                                                                                              | |
        !| |                                                                                                              | |
        !| |    * CONSTRAINED PROBLEMS:                                                                                   | |
        !| |                                                                                                              | |
        !| |                     objective functions       constraint functions      starting point        dimension      | |
        !| |                   -----------------------   ------------------------   -----------------    -------------    | |
        !| |    Problem 16:            2, 7                         1                      16                  2          | |
        !| |    Problem 17:            3, 9                         2                      17                  4          | |
        !| |    Problem 18:            10, 12                       3                      18                varies       | |
        !| |    Problem 19:            2, 6, 7                      1                      19                  2          | |
        !| |    Problem 20:            3, 4, 9                      2                      20                  4          | |
        !| |    Problem 21:            4, 10, 16                    3                      21                varies       | |
        !| |                                                                                                              | |
        !| |                                                                                                              | |
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*
        
        
        !--------------------------------------------------------------------------------
        ! ----------------------------------------------------------------------------- |
        ! |                  INFORMATION ABOUT PARAMETERS:                            | |
        ! ----------------------------------------------------------------------------- |
        !--------------------------------------------------------------------------------       
        
      
        !****************** GLOBAL PARAMETERS *******************************************
      
        INTEGER, PARAMETER :: number_of_obj = 3                         ! The number of objectives in the multiobjective problem
                                                                        ! IF 'number_of_obj'=1 THEN we have a single objective case.
        INTEGER, PARAMETER :: number_of_const = 1                       ! The number of constraints in the multiobjective problem
                                                                        ! IF 'number_of_obj'=0 THEN we have an unconstrained problem    
                                                                        
        INTEGER, PARAMETER :: number_of_func = number_of_obj + number_of_const       ! the sum of objective and constraint functions                                                               

        !---------------------------------------------------------------------------------------------------------------------------------
        ! DC component lists where the first 'number_of_obj' places tell the DC components f1 for the objectives
        ! and the latter 'number_of_const' places tell the DC components g1 for the constraints         
        INTEGER, DIMENSION(number_of_func), PARAMETER :: problem1 = (/ 4,10,16,3 &     ! DC components f1 of objective functions used
                                                                    & /)          ! DC components g1 of constraint functions used
        INTEGER, DIMENSION(number_of_func), PARAMETER :: problem2 = problem1      ! DC components f2 and g2 of objective and constraint functions used 
        
        INTEGER, PARAMETER :: startpoint = 21                           ! Starting point used

        !-------------------------------------------------------------------------------------------------
      
        INTEGER, PARAMETER :: user_n = 500                             ! The number of variables in the problem
        
        !-------------------------------------------------------------------------------------------------           

 
        !__________________________________________________________________________________________
        !>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>>
        !****************** PARAMETRES NEEDED ONLY IN MAIN ITERATION AlGORITHM ********************
        
        INTEGER, PARAMETER :: user_size_b1 = MIN((user_n + 5)*number_of_func,1000)  ! The biggest possible size of the bundle B_1 
                                                                                    ! If user_size_b1 <= 0 then DEFAULT value MIN((user_n+5)*number_of_obj,1000) is used 
                                                                   
        INTEGER, PARAMETER :: user_size_b2 = 3                      ! The biggest possible size of the bundle B_2 
                                                                    ! If user_size_b2 <= 0 then DEFAULT value 3 is used  
            
        REAL(KIND=dp), PARAMETER :: user_m = 0.01_dp                ! The descent parameter ('m2' in the paper)  If user_m <= 0.0_dp .OR. user_m >= 1.0_dp 
                                                                    !                                            then DEFAULT value 0.2_dp is used
                                                                   
        REAL(KIND=dp), PARAMETER :: user_m_extra = 0.1_dp           ! The extra descent parameter ('m3' in the paper)  If user_m_extra <= user_m .OR. user_m >= 1.0_dp 
                                                                    !                                                  then DEFAULT value '2*user_m' is used if it is smaller than 1.0_dp
                                                                    !                                                  if it is .NOT. then the value 0.5(user_m + 1.0_dp) is chosen
                                                               
        REAL(KIND=dp), PARAMETER :: user_r_dec = -0.99_dp           ! The decrease parameter     
        
        !If user_r_dec <= 0.0_dp .OR. user_r_dec >= 1.0_dp then DEFAULT value is used.
        !                               
        !   DEFAULT value:                          
        !     If user_n < 10:           user_r_dec = 0.75_dp    
        !     If 10 <= user_n < 300:    user_r_dec = n/(n+5)
        !     If user_n >= 300:         user_r_dec = 0.99_dp

        
        REAL(KIND=dp), PARAMETER :: user_r_inc = (10.0_dp)**(10)    ! The increase parameter     If user_r_inc <= 1.0_dp 
                                                                    !                            then DEFAULT value (10.0_dp)**10 is used
                                                                 
        REAL(KIND=dp), PARAMETER :: user_c1 = 0.5_dp                ! The decrease parameter  If user_c1 <= 0.0_dp .OR. user_c1 >= 1.0_dp 
                                                                    !                        then DEFAULT value 0.5_dp is used  
        REAL(KIND=dp), PARAMETER :: user_c2 = -0.3_dp               ! The decrease parameter  If user_c2 <= 0.0_dp .OR. user_c2 >= 1.0_dp 
                                                                    !                        then DEFAULT values are used  (formed using values: 0.5, 0.25, 0.1, 0.01, 0.001 )                                                              
        REAL(KIND=dp), PARAMETER :: user_c3 = -0.1_dp               ! The decrease parameter  If user_c3 <= 0.0_dp .OR. user_c3 >= 1.0_dp 
                                                                    !                        then DEFAULT value 0.1_dp is used                                                                  
                                                                
        !____________________________________________________________________________________________                                                       
        !>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<
        !****************** PARAMETRES NEEDED ONLY IN CLARKE STATIONARY AlGORITHM *******************
        
        
         INTEGER, PARAMETER :: user_size = (user_n + 5) * 2       ! The biggest possible bundle size for the bundle used in 'Clarke stationary' algorithm
                                                                  ! If user_size_b1 < 2 then DEFAULT value '(user_n+5)*2' is used 
         REAL(KIND=dp), PARAMETER :: user_m_clarke = 0.01_dp        ! The descent parameter ('m1' in the paper): If user_m_clarke <= 0.0_dp .OR. user_m_clarke >= 1.0_dp 
                                                                    !                        then DEFAULT value 0.01_dp is used
                                                                    
          
         REAL(KIND=dp), PARAMETER :: user_step_tol = (10.0_dp)**(-4)  ! The step-length tolerance (i.e. proximity measure in the paper): If user_set_tol <= 0.0_dp 
                                                                      !                                                                  then DEFAULT value (10.0_dp)**(-4) is used 
                                                                    
        !________________________________________________________________________________
        !>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<
        !****************** PARAMETRES NEEDED IN STOPPING CONDITIONS ********************

        ! ** The stopping tolerance **     
        REAL(KIND=dp), PARAMETER :: user_crit_tol = (10.0_dp)**(-5)     ! The stopping tolerance: If user_crit_tol <= 0.0_dp .OR. user_crit_tol >= 1.0_dp 
                                                                        ! then DEFAULT value (10.0_dp)**(-5) is used
                                                                        
        REAL(KIND=dp), PARAMETER :: user_dec_tol = 1.0_dp               ! If user_dec_tol > 0.0_dp
                                                                        ! then DEFAULT value   0.0_dp         is used when user_n<=100
                                                                        !                    -(10.0_dp)**(-5} is used when 100<user_n<250
                                                                        !                    -(10.0_dp)**(-4) is used when 250<=user_n
                                                                                
        ! ** The small number **
        REAL(KIND=dp), PARAMETER :: user_eps = 0.00005_dp               ! The enlargement parameter: If user_eps <= 0.0_dp then DEFAULT value 0.00005_dp is used
                
       
        !________________________________________________________________________________
        !>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<>><<
        
        
        
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !| |                                                                          | |
        !| |     To EXECUTE the bundle algorithm MDBDC the USER needs to DETERMINE:   | | 
        !| |                                                                          | |       
        !| |        f1(y, problem1)     - value of DC component f_1 at a point y      | |
        !| |        f2(y, problem2)     - value of DC component f_2 at a point y      | |
        !| |                                                                          | |
        !| |        subgradient_f1(y, problem1)   - subgradient of f_1 at y           | |
        !| |        subgradient_f2(y, problem2)   - subgradient of f_2 at y           | |
        !| |                                                                          | |
        !| |        g1(y, problem1)     - value of DC component g_1 at a point y      | |
        !| |        g2(y, problem2)     - value of DC component g_2 at a point y      | |
        !| |                                                                          | |
        !| |        subgradient_g1(y, problem1)   - subgradient of g_1 at y           | |
        !| |        subgradient_g2(y, problem2)   - subgradient of g_2 at y           | |       
        !| |                                                                          | |
        !| |                                                                          | |       
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*
    
        
        CONTAINS
        
        
        
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !| |                                                                          | |
        !| |                        OBJECTIVE FUNCTIONS:                              | |
        !| |                                                                          | |
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*       

        
        !********************************************************************************
        !                                                                               |
        !              FUNCTION VALUES OF THE DC COMPONENTS f_1 AND f_2                 |
        !                                                                               |
        !********************************************************************************

           FUNCTION f1(y, problem1, factor) RESULT(f)        
                !
                ! Calculates the function value of the DC component f_1 at a point 'y'.
                ! Variable 'problem1' identifies the objective function used
                ! 'factor' is a scaling coefficient for objective               
                !
                ! NOTICE: The dimension of 'y' has to be 'user_n'.
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: y    ! a point where the function value of the DC component f_1 is calculated
                INTEGER, INTENT(IN) :: problem1                 ! the objective function selected for the DC component f_1
                REAL(KIND=dp), INTENT(IN) :: factor             ! the scaling factor for objective
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp) :: f                              ! the function value of the DC component f_1 at a point 'y'
                REAL(KIND=dp), DIMENSION(4) :: a                ! help variables
                REAL(KIND=dp) :: apu, largest                   ! help variables
                REAL(KIND=dp), PARAMETER :: zero=0.0_dp, one=1.0_dp
                REAL(KIND=dp), PARAMETER :: two=2.0_dp
                REAL(KIND=dp) :: a1, a2, a3, a4, a5, a6 
                INTEGER :: i, j, ind                                    ! help variable
                
                REAL(KIND=dp), DIMENSION(user_n) :: f_i    ! only first n-1 places are used            
                
                SELECT CASE(problem1)


                   !-------------------------------------
                   !           Problem   1
                   !-------------------------------------
                   CASE(1)  
                     a1 = y(1)**4 + y(2)**2
                     a2 = (2.0_dp - y(1))**2 + (2.0_dp - y(2))**2 
                     a3 = 2.0_dp * EXP(-y(1)+y(2))

                     a4 = y(1)**2 - 2.0_dp * y(1) + y(2)**2 - 4.0_dp * y(2) + 4.0_dp
                     a5 = 2.0_dp * y(1)**2 - 5.0_dp * y(1) + y(2)**2 - 2.0_dp * y(2) 
                     a5 = a5 + 4.0_dp
                     a6 = y(1)**2 + 2.0_dp * y(2)**2 - 4.0_dp * y(2) + 1.0_dp               
                
                     f = MAX(a1,a2,a3) + a4 + a5 + a6                  

                     f = f * factor                  
                     
                   !-------------------------------------  
                   
                   !-------------------------------------
                   !           Problem   2
                   !-------------------------------------
                   CASE(2)  
                     IF ((ABS(y(1)) - y(2) ) > 0.0_dp) THEN
                        f = ABS(y(1) - 1.0_dp) + 200.0_dp * (ABS(y(1)) - y(2) )
                     ELSE
                        f = ABS(y(1) - 1.0_dp)
                     END IF
                     
                     f = f * factor                  
                   !-------------------------------------    
                  
                   !-------------------------------------
                   !           Problem   3
                   !-------------------------------------
                   CASE(3)
                     f = ABS(y(1) - 1.0_dp)
                     f = f + ABS(y(3) - 1.0_dp )                
                     f = f + 4.95_dp * ( ABS(y(2) + y(4) - 2.0_dp ) )
                     f = f + 10.1_dp * ( ABS(y(2)-1.0_dp) + ABS(y(4)-1.0_dp) )      
                
                     IF ((ABS(y(1)) - y(2) ) > 0.0_dp) THEN
                        f = f + 200.0_dp * (ABS(y(1)) - y(2) )
                     END IF 
                
                     IF ((ABS(y(3)) - y(4) ) > 0.0_dp) THEN
                        f = f + 180.0_dp * (ABS(y(3)) - y(4) )
                     END IF  

                     f = f * factor
                   !-------------------------------------                  
                   
                   !-------------------------------------
                   !           Problem   4
                   !-------------------------------------
                   CASE(4)
                     apu = ABS(y(1))
                     DO i = 2, user_n
                       IF (apu < ABS(y(i)) ) THEN 
                          apu = ABS(y(i))
                       END IF
                     END DO
                     f = user_n * apu 

                     f = f * factor
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   5
                   !-------------------------------------
                   CASE(5)
                     apu = ABS(summa(y,1)) 
                     DO j = 2, 20
                       IF (apu <= ABS(summa(y,j)) ) THEN 
                         apu = ABS(summa(y,j))
                       END IF
                     END DO
                     f = 20.0_dp * apu

                     f = f * factor                  
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   6
                   !-------------------------------------
                   CASE(6)
                     f = y(2) + 0.1_dp * (y(1)**2 + y(2)**2)  

                     IF (-y(2) <= 0.0_dp ) THEN 
                         f = f
                     ELSE
                         f = f + 10.0_dp * (-y(2))
                     END IF
                     
                     f = f * factor                  
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   7
                   !-------------------------------------
                   CASE(7)
                     f = ABS(y(1)-1.0_dp) + 200.0_dp * MAX(0.0_dp, ABS(y(1))-y(2) )                 
                
                     a1 = y(1)**2 + y(2)**2 + ABS(y(2))
                     a2 = y(1) + y(1)**2 + y(2)**2 + ABS(y(2)) - 0.5_dp
                     a3 = ABS( y(1) - y(2) ) + ABS(y(2)) - 1.0_dp
                     a4 = y(1) + y(1)**2 + y(2)**2
                
                     f = f + 10.0_dp * MAX(a1,a2,a3,a4)
                     
                     f = f * factor                  
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   8
                   !-------------------------------------                    
                   CASE(8)
                     f = 9.0_dp - 8.0_dp * y(1) -6.0_dp * y(2) - 4.0_dp * y(3)  
                     f = f + 2.0_dp * ABS(y(1)) + 2.0_dp * ABS(y(2))+ 2.0_dp * ABS(y(3)) 
                     f = f + 4.0_dp * y(1)**2 + 2.0_dp * y(2)**2 + 2.0_dp * y(3)**2 

                     a(1) = y(1) + y(2) + 2.0_dp * y(3) - 3.0_dp
                     a(2) = -y(1)
                     a(3) = -y(2)
                     a(4) = -y(3)
                
                     f = f + 10.0_dp * MAX(0.0_dp, a(1), a(2), a(3), a(4))
                     
                     f = f * factor                  
                   !------------------------------------
                    
                   !-------------------------------------
                   !           Problem   9
                   !-------------------------------------                    
                   CASE(9)
                     f = y(1)**2 + (y(1)-1.0_dp)**2 + 2.0_dp*(y(1)-2.0_dp)**2 
                     f = f + (y(1)-3.0_dp)**2 + 2.0_dp * y(2)**2 + (y(2)-1.0_dp)**2 
                     f = f + 2.0_dp*(y(2)-2.0_dp)**2 + y(3)**2 + (y(3)-1.0_dp)**2 
                     f = f + 2.0_dp*(y(3)-2.0_dp)**2 + (y(3)-3.0_dp)**2 
                     f = f + 2.0_dp*y(4)**2 + (y(4)-1.0_dp)**2 + 2.0_dp*(y(4)-2.0_dp)**2 

                     f = f * factor                  
                   !------------------------------------
                               
                   !-------------------------------------
                   !           Problem   10
                   !-------------------------------------
                   CASE(10)
                     f = 0.0_dp
                     DO i = 1, user_n
                       f =  f + y(i)**2                             
                     END DO      

                     f = f * factor
                   !-------------------------------------

                   !-------------------------------------
                   !           Problem   11
                   !-------------------------------------
                   CASE(11)
                     f = 0.0_dp
                     f = f + 4.0_dp * ABS(y(1)) 
                     f = f + 2.0_dp * ABS(y(2)) 
                     f = f + 2.0_dp * ABS(y(3)) 
                     f = f - 33.0_dp*y(1) + 16.0_dp*y(2) - 24.0_dp*y(3)
                     f = f + 100.0_dp*MAX(0.0_dp, 2.0_dp*ABS(y(2))-3.0_dp*y(1)-7.0_dp)
                     f = f + 100.0_dp*MAX(0.0_dp, ABS(y(3))-4.0_dp*y(1)-11.0_dp)
                     
                     f = f * factor                  
                   !-------------------------------------   
                   
                   !-------------------------------------
                   !           Problem   12
                   !-------------------------------------
                   CASE(12)                
                     f = 0.0_dp
                     
                     DO i = 1, user_n
                        f = f + ABS(y(i))
                     END DO
                     
                     DO i = 1, user_n
                        apu = 2.0_dp * (y(i)*y(i) -y(i) -1.0_dp)
                        IF (apu > 0.0_dp) THEN
                           f = f + 10.0_dp*apu
                        END IF
                     END DO
                     
                     f = f * factor
                 
                   !-------------------------------------   
                   
                   !-------------------------------------
                   !           Problem   13
                   !-------------------------------------
                   CASE(13)             
                    
                     f = 0.0_dp
                     
                     DO i = 1, user_n-1
                        f = f + ABS(y(i)+y(i+1))
                     END DO
                     
                     DO i = 1, user_n-2
                        f = f + ABS(y(i)+y(i+2))
                     END DO
                     
                     f = f + ABS(y(1)+y(9)) 
                     f = f + ABS(y(1)+y(10)) 
                     f = f + ABS(y(2)+y(10)) 
                     f = f + ABS(y(1)+y(5)) 
                     f = f + ABS(y(4)+y(7)) 
                     
                     apu = 0.0_dp
                     DO i = 1, user_n
                        apu = apu + y(i)
                     END DO
                     
                     IF ((apu-1.0_dp)>0.0_dp) THEN 
                        f = f + 10.0_dp * (apu-1.0_dp)
                     END IF
                     
                     DO i = 1, user_n
                        IF (-y(i) > 0.0_dp) THEN 
                            f = f + 10.0_dp*(-y(i))
                        END IF
                     END DO
                    
                     f = f * factor                 
                   !-------------------------------------   

                   !-------------------------------------
                   !           Problem   14
                   !-------------------------------------
                   CASE(14)         

                     f = 0.0_dp
                     
                     apu = 0.0_dp
                     DO j =1, user_n
                        apu = apu + y(j)/(1+j-1.0_dp)
                     END DO
                     largest = ABS(apu) 
                     
                     DO i = 2, user_n
                        apu = 0.0_dp
                        DO j = 1, user_n
                          apu = apu + y(j)/(i+j-1.0_dp)
                        END DO
                        IF (ABS(apu) > largest) THEN 
                           largest = ABS(apu)
                        END IF
                     END DO
                     
                     f = user_n * largest   

                     f = f * factor                  
                   
                   !-------------------------------------   

                   !-------------------------------------
                   !           Problem   15
                   !-------------------------------------
                   CASE(15)     

                     f = 0.0_dp
                     
                     DO i = 1, user_n-1
                         a(1) = y(i)**4 + y(i+1)**2
                         a(2) = (2.0_dp-y(i))**2 + (2.0_dp-y(i+1))**2
                         a(3) = 2.0_dp*EXP(-y(i)+y(i+1))
                         IF (a(1)>a(2)) THEN
                            IF (a(1)>a(3)) THEN
                              f_i(i) = a(1)
                            ELSE
                              f_i(i) = a(3)
                            END IF
                         ELSE
                            IF (a(2)>a(3)) THEN
                              f_i(i) = a(2)
                            ELSE
                              f_i(i) = a(3)
                            END IF      
                         END IF
                     END DO
                     
                     largest = f_i(1)
                     ind = 1
                     DO i =2, user_n-1
                         IF (f_i(i) > largest) THEN 
                            largest = f_i(i)
                            ind = i
                         END IF
                     END DO
                     
                     f = (user_n-1.0_dp)*largest
                     
                     f = f * factor
                   
                   !-------------------------------------                      
                   
               
                   !-------------------------------------
                   !              Problem 16:
                   !          Chained Crescent I
                   !-------------------------------------
                   CASE(16)
                     f = 0.0_dp
                     
                     apu = 0.0_dp
                     DO i = 1, user_n-1
                        apu = apu + y(i)**2 + (y(i+1)-1.0_dp)**2
                        apu = apu + y(i+1) -1.0_dp
                     END DO
                     
                    apu = 2.0_dp *apu
                    IF (apu > 0.0_dp) THEN
                       f = apu
                    END IF
                    
                    f = f * factor
                   !------------------------------------- 
               

                
                END SELECT
                
           END FUNCTION f1      
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
           
           FUNCTION f2(y, problem2, factor) RESULT(f)           
                !
                ! Calculates the function value of DC component f_2 at a point 'y'.
                ! Variable 'problem2' identifies the objective function used.
                ! 'factor' is a scaling coefficient for the objective
                !
                ! NOTICE: The dimension of 'y' has to be 'user_n'.
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: y    ! a point where the function value of the DC component f_2 is calculated
                INTEGER, INTENT(IN) :: problem2                 ! the objective function selected for the DC component f_2   
                REAL(KIND=dp), INTENT(IN) :: factor             ! the scaling factor for objective              
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp) :: f                          ! the function value of the DC component f_2 at a point 'y'
                REAL(KIND=dp) :: apu, largest               
                REAL(KIND=dp) :: a, b              ! help variables
                REAL(KIND=dp), PARAMETER :: zero=0.0_dp, one=1.0_dp
                REAL(KIND=dp), PARAMETER :: two=2.0_dp
                REAL(KIND=dp), DIMENSION(3) :: aa
                REAL(KIND=dp), DIMENSION(user_n-1) :: f_i
                REAL(KIND=dp), DIMENSION(user_n) :: term
                REAL(KIND=dp) :: a4, a5, a6
                INTEGER :: i, j, ind                                     ! help variable
                
                
                SELECT CASE(problem2)
                
                
                   !-------------------------------------
                   !           Problem   1
                   !-------------------------------------
                   CASE(1)
                     a4 = y(1)**2 - 2.0_dp * y(1) + y(2)**2 - 4.0_dp * y(2) + 4.0_dp
                     a5 = 2.0_dp * y(1)**2 - 5.0_dp * y(1) + y(2)**2 - 2.0_dp * y(2) 
                     a5 = a5 + 4.0_dp
                     a6 = y(1)**2 + 2.0_dp * y(2)**2 - 4.0_dp * y(2) + 1.0_dp
                
                     f = MAX( (a4 + a5) , (a4 + a6) , (a5 + a6) )                  

                     f = f * factor
                   !-------------------------------------  
                   
                   !-------------------------------------
                   !           Problem   2
                   !-------------------------------------                  
                   CASE(2)
                     f = 100.0_dp * ( ABS(y(1)) - y(2) )
                     
                     f = f * factor                  
                   !-------------------------------------                   
                   
                   !-------------------------------------
                   !           Problem   3
                   !-------------------------------------
                   CASE(3)
                     f = 4.95_dp * ( ABS(y(2) - y(4)) )
                     f = f + 90.0_dp * ( ABS(y(3)) - y(4) )
                     f = f + 100.0_dp * ( ABS(y(1)) - y(2) )

                     f = f * factor                  
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   4
                   !-------------------------------------
                   CASE(4)
                     f = 0.0_dp
                     DO i = 1, user_n
                        f = f + ABS(y(i))                                   
                     END DO  

                     f = f * factor
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   5
                   !-------------------------------------
                   CASE(5)
                     f = 0.0_dp
                     DO j = 1, 20
                        f = f + ABS( summa(y,j) )
                     END DO 

                     f = f * factor                  
                   !-------------------------------------                  
                   
                   !-------------------------------------
                   !           Problem   6
                   !-------------------------------------
                   CASE(6)
                     f = 0.0_dp
                
                     IF (y(1) <= 0.0_dp ) THEN 
                         f = f - y(1)
                     ELSE
                         f = f + y(1) 
                     END IF
                
                     IF (y(2) <= 0.0_dp ) THEN 
                         f = f - y(2)
                     ELSE
                         f = f + y(2)
                     END IF   

                     f = f * factor                  
                   !-------------------------------------                   
                   
                   !-------------------------------------
                   !           Problem   7
                   !-------------------------------------
                   CASE(7)
                     f = 10.0_dp * ( y(1)**2 + y(2)**2 + ABS(y(2)) )
                     f = f + 100.0_dp * ( ABS(y(1)) - y(2) )

                     f = f * factor                  
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   8
                   !-------------------------------------                   
                   CASE(8)
                     f = ABS( y(1) -y(2) ) + ABS( y(1) -y(3) )
                     
                     f = f * factor                  
                   !------------------------------------
                    
                   !-------------------------------------
                   !           Problem   9
                   !-------------------------------------                    
                   CASE(9)
                     f = 0.0_dp 

                     a = (y(1) -2.0_dp)**2 + y(2)**2
                     b = (y(3) -2.0_dp)**2 + y(4)**2
                
                     IF ( a >= b) THEN 
                       f = f + a
                     ELSE 
                       f = f + b
                     END IF

                     a = (y(1) -2.0_dp)**2 + (y(2)-1.0_dp)**2
                     b = (y(3) -2.0_dp)**2 + (y(4)-1.0_dp)**2
                
                     IF ( a >= b) THEN 
                       f = f + a
                     ELSE 
                       f = f + b
                     END IF

                     a = (y(1) -3.0_dp)**2 + y(2)**2
                     b = (y(3) -3.0_dp)**2 + y(4)**2
                
                     IF ( a >= b) THEN 
                       f = f + a
                     ELSE 
                       f = f + b
                     END IF

                    a = (y(1))**2 + (y(2)-2.0_dp)**2
                    b = (y(3))**2 + (y(4)-2.0_dp)**2
                
                     IF ( a >= b) THEN 
                       f = f + a
                     ELSE 
                       f = f + b
                     END IF

                     a = (y(1)-1.0_dp)**2 + (y(2)-2.0_dp)**2
                     b = (y(3)-1.0_dp)**2 + (y(4)-2.0_dp)**2
                
                     IF ( a >= b) THEN 
                       f = f + a
                     ELSE 
                       f = f + b
                     END IF            

                     f = f * factor                  
                   !------------------------------------  
                   
                   !-------------------------------------
                   !           Problem   10
                   !-------------------------------------
                   CASE(10)
                     f = 0.0_dp
                     DO i = 2, user_n
                         f = f + ABS( y(i) - y(i-1) ) 
                     END DO
                     
                     f = f * factor
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   11
                   !-------------------------------------
                   CASE(11)
                     apu = - 7.0_dp * y(1) 
                     apu = apu + 2.0_dp * ABS(y(2)) - ABS(y(3)) - 18.0_dp
                     f = 20.0_dp * apu
                     
                     f = f * factor                  
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   12
                   !-------------------------------------
                   CASE(12)                
                     f = 0.0_dp

                     DO i = 1, user_n
                        apu = (y(i)*y(i) -y(i) -1.0_dp)
                        f = f + 10.0_dp * apu
                     END DO
                     
                     term = 0.0_dp
                     DO i = 1, user_n
                        DO j = 1, user_n 
                           IF (i /= j) THEN
                               term(i) = term(i) + ABS(y(j))
                           END IF
                        END DO
                     END DO
                     
                     largest = term(1)
                     ind = 1 
                     
                     DO i = 2, user_n
                        IF (term(i) > largest) THEN 
                           largest = term(i)
                           ind = i
                        END IF
                     END DO
                     
                     f = f + largest
                     
                    f = f * factor                   
                   !-------------------------------------   
                   
                   !-------------------------------------
                   !           Problem   13
                   !-------------------------------------
                   CASE(13)             
                     f = 0.0_dp
                     
                     DO i = 1, user_n-1
                        f  = f + ABS(y(i)) + ABS(y(i+1))
                     END DO
                     
                     DO i = 1, user_n-2
                        f  = f + ABS(y(i)) + ABS(y(i+2))
                     END DO                  
                     
                     f = f + 3.0_dp*ABS(y(1))
                     f = f + ABS(y(2))
                     f = f + ABS(y(4))
                     f = f + ABS(y(5))
                     f = f + ABS(y(7))
                     f = f + ABS(y(9))
                     f = f + 2.0_dp * ABS(y(10))    
    
                     f = f * factor 
                   !-------------------------------------   

                   !-------------------------------------
                   !           Problem   14
                   !-------------------------------------
                   CASE(14)         
                  
                     f = 0.0_dp
                     
                     DO i = 1, user_n
                        apu = 0.0_dp
                        DO j = 1, user_n
                          apu = apu + y(j)/(i+j-1.0_dp)
                        END DO
                        f = f + ABS(apu)
                     END DO
                     
                    f = f * factor                   
                  
                   !-------------------------------------   

                   !-------------------------------------
                   !           Problem   15
                   !-------------------------------------
                   CASE(15)     
                     f = 0.0_dp
                     
                     DO i = 1, user_n-1
                         aa(1) = y(i)**4 + y(i+1)**2
                         aa(2) = (2.0_dp-y(i))**2 + (2.0_dp-y(i+1))**2
                         aa(3) = 2.0_dp*EXP(-y(i)+y(i+1))
                         IF (aa(1)>aa(2)) THEN
                            IF (aa(1)>aa(3)) THEN
                              f_i(i) = aa(1)
                            ELSE
                              f_i(i) = aa(3)
                            END IF
                         ELSE
                            IF (aa(2)>aa(3)) THEN
                              f_i(i) = aa(2)
                            ELSE
                              f_i(i) = aa(3)
                            END IF      
                         END IF
                     END DO

                     DO i =1, user_n-1
                         f = f + f_i(i)
                     END DO 

                     f = f * factor
                    
                   !-------------------------------------                  
                   
                   !-------------------------------------
                   !            Problem 16:
                   !         Chained Crescent I
                   !-------------------------------------
                   CASE(16)
                     f = 0.0_dp

                     DO i = 1, user_n-1
                        f = f + y(i)**2 + (y(i+1)-1.0_dp)**2
                        f = f + y(i+1) -1.0_dp
                     END DO
                     
                     f = f * factor                  
                   !-------------------------------------                  
            
                
                END SELECT              
            


           END FUNCTION f2
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
           

        !********************************************************************************
        !                                                                               |
        !                SUBGRADIENTS OF THE DC COMPONENTS f_1 AND f_2                  |
        !                                                                               |       
        !********************************************************************************       
        
           FUNCTION subgradient_f1(y, problem1, factor) RESULT(grad)           
                !
                ! Calculates a subgradient of the DC component f_1 at a point 'y'.
                ! Variable 'problem1' identifies the objective function used.
                ! 'factor' is a scaling coefficient for the objective
                !
                ! NOTICE: * The dimension of 'y' has to be 'user_n'.
                !         * The dimension of 'grad' is also 'user_n'.
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: y    ! a point where the subgradient of the DC component f_1 is calculated
                INTEGER, INTENT(IN) :: problem1                 ! the objective function selected for the DC component f_1 
                REAL(KIND=dp), INTENT(IN) :: factor             ! the scaling factor for objective              
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp), DIMENSION(SIZE(y)) :: grad       ! the subgradient of the DC component f_1 at a point 'y'
                REAL(KIND=dp), DIMENSION(4) :: a
                REAL(KIND=dp), DIMENSION(5) :: b                ! help variable
                REAL(KIND=dp), PARAMETER :: zero=0.0_dp, one=1.0_dp
                REAL(KIND=dp), PARAMETER :: two=2.0_dp
                REAL(KIND=dp) :: a1, a2, a3
                REAL(KIND=dp) :: apu                                ! help varaible
                REAL(KIND=dp) :: largest                            ! help varaible
                REAL(KIND=dp), DIMENSION(user_n) :: abs_sign        ! help varaible
                REAL(KIND=dp), DIMENSION(user_n) :: f_i             ! help varaible (only first n-1 placs are used
                INTEGER :: i, j, ind                                ! help variable
                


                SELECT CASE(problem1)
                

 
                   !-------------------------------------
                   !           Problem   1
                   !-------------------------------------                    
                   CASE(1)              
                     a1 = y(1)**4 + y(2)**2
                     a2 = (2.0_dp - y(1))**2 + (2.0_dp - y(2))**2 
                     a3 = 2.0_dp * EXP(-y(1)+y(2))              
                
                     IF (a1 >= a2) THEN
                       IF (a1 >= a3) THEN
                         grad(1) = 4.0_dp * y(1)**3 
                         grad(2) = 2.0_dp * y(2)                    
                       ELSE
                         grad(1) = -2.0_dp * EXP(-y(1)+y(2))
                         grad(2) = 2.0_dp * EXP(-y(1)+y(2))                 
                       END IF               
                     ELSE
                       IF (a2 >= a3) THEN
                         grad(1) = 2.0_dp * y(1) - 4.0_dp
                         grad(2) = 2.0_dp * y(2) - 4.0_dp                   
                       ELSE
                         grad(1) = -2.0_dp * EXP(-y(1)+y(2))
                         grad(2) = 2.0_dp * EXP(-y(1)+y(2))                 
                       END IF                   
                     END IF
                
                     grad(1) = grad(1) + 2.0_dp * y(1) - 2.0_dp
                     grad(1) = grad(1) + 4.0_dp * y(1) - 5.0_dp
                     grad(1) = grad(1) + 2.0_dp * y(1)
  
                     grad(2) = grad(2) + 2.0_dp * y(2) - 4.0_dp
                     grad(2) = grad(2) + 2.0_dp * y(2) - 2.0_dp
                     grad(2) = grad(2) + 4.0_dp * y(2) - 4.0_dp                

                     grad = grad * factor
                     
                   !-------------------------------------      
                   
                   !-------------------------------------
                   !           Problem   2
                   !-------------------------------------                    
                   CASE(2)              
                     IF ( y(1) <= 1.0_dp ) THEN
                        grad(1) = -1.0_dp
                     ELSE
                        grad(1) = 1.0_dp
                     END IF
                
                     IF( ( ABS(y(1))-y(2) ) > 0.0_dp  ) THEN
                        IF(y(1) <= 0.0_dp) THEN 
                          grad(1) = grad(1) - 200.0_dp
                        ELSE
                          grad(1) = grad(1) + 200.0_dp  
                        END IF                   
                        grad(2) = -200.0_dp
                     ELSE        
                        grad(1) = grad(1) 
                        grad(2) = 0.0_dp
                     END IF
                     
                     grad = grad * factor                    
                   !-------------------------------------                   

                   !-------------------------------------
                   !           Problem   3
                   !-------------------------------------  
                   CASE(3)
                     grad(1) = 0.0_dp
                     grad(2) = 0.0_dp
                     grad(3) = 0.0_dp
                     grad(4) = 0.0_dp
                
                     IF ( y(1) <= 1.0_dp ) THEN
                       grad(1) = -1.0_dp
                     ELSE
                       grad(1) = 1.0_dp
                     END IF
                
                     IF( ( ABS(y(1))-y(2) ) > 0.0_dp  ) THEN
                       IF(y(1) <= 0.0_dp) THEN 
                          grad(1) = grad(1) - 200.0_dp
                       ELSE
                          grad(1) = grad(1) + 200.0_dp  
                       END IF                    
                       grad(2) = -200.0_dp
                     ELSE        
                       grad(1) = grad(1) 
                       grad(2) = 0.0_dp
                     END IF

                     IF( ( ABS(y(3))-y(4) ) > 0.0_dp  ) THEN
                       IF(y(3) <= 0.0_dp) THEN 
                          grad(3) =  - 180.0_dp
                       ELSE
                          grad(3) =  180.0_dp  
                       END IF                    
                       grad(4) = -180.0_dp
                     ELSE        
                       grad(3) = 0.0_dp 
                       grad(4) = 0.0_dp
                     END IF
                
                     IF ( y(3) <= 1.0_dp ) THEN
                       grad(3) = grad(3) - 1.0_dp
                     ELSE
                       grad(3) = grad(3) + 1.0_dp
                     END IF         
                
                     IF ( y(2) <= 1.0_dp ) THEN
                       grad(2) = grad(2) - 10.1_dp  
                     ELSE
                       grad(2) = grad(2) + 10.1_dp
                     END IF         
                
                     IF ( y(4) <= 1.0_dp ) THEN
                       grad(4) = grad(4) - 10.1_dp  
                     ELSE
                       grad(4) = grad(4) + 10.1_dp
                     END IF 
                
                     IF ( (y(2) + y(4) - 2.0_dp) <= 0.0_dp ) THEN
                       grad(2) = grad(2) - 4.95_dp
                       grad(4) = grad(4) - 4.95_dp  
                     ELSE
                       grad(2) = grad(2) + 4.95_dp
                       grad(4) = grad(4) + 4.95_dp
                     END IF 
                     
                     grad = grad * factor                    
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   4
                   !-------------------------------------  
                   CASE(4)
                     grad = 0.0_dp 
                     apu = ABS(y(1))
                     ind = 1
                
                     DO i = 2, user_n
                       IF (apu < ABS(y(i))) THEN 
                         apu = ABS(y(i))
                         ind = i
                       END IF
                     END DO
                
                     IF (y(ind) <= 0.0_dp) THEN 
                         grad(ind) = - user_n
                     ELSE 
                         grad(ind) = user_n
                    END IF  

                     grad = grad * factor
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   5
                   !-------------------------------------  
                   CASE(5)
                     grad = 0.0_dp 
                     apu = ABS(summa(y,1)) 
                     ind = 1

                     DO j = 2, 20
                       IF (apu <= ABS(summa(y,j)) ) THEN 
                         apu = ABS(summa(y,j))
                         ind = j
                       END IF
                     END DO

                     DO i = 1, user_n
                       grad(i) = (0.05_dp * ind)**(i-1)
                     END DO             
                
                     IF ( summa(y,ind) <= 0.0_dp ) THEN 
                         grad = -20.0_dp * grad
                     ELSE
                         grad = 20.0_dp * grad
                     END IF 
                     
                     grad = grad * factor                    
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   6
                   !-------------------------------------  
                   CASE(6)
                     grad = 0.0_dp 
                
                     grad(1) = 0.2_dp * y(1)
                
                     grad(2) = 1.0_dp + 0.2_dp * y(2) 

                     IF ( -y(2) <= 0.0_dp ) THEN 
                         grad(2) = grad(2)
                     ELSE
                         grad(2) = grad(2) - 10.0_dp
                     END IF        

                     grad = grad * factor                    
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   7
                   !-------------------------------------  
                   CASE(7)
                     grad = 0.0_dp 
                
                     IF ( (y(1)-1.0_dp) <= 0.0_dp ) THEN 
                       grad(1) = -1.0_dp
                     ELSE
                       grad(1) = 1.0_dp
                     END IF 

                     a(1) = y(1)**2 + y(2)**2 + ABS(y(2))
                     a(2) = y(1) + y(1)**2 + y(2)**2 + ABS(y(2)) - 0.5_dp
                     a(3) = ABS( y(1) - y(2) ) + ABS(y(2)) - 1.0_dp
                     a(4) = y(1) + y(1)**2 + y(2)**2                
                
                     ind = 1
                     DO i = 2, 4
                       IF ( a(ind) < a(i) ) THEN 
                         ind = i
                       END IF 
                     END DO
                
                     IF (ind == 1) THEN
                       grad(1) = grad(1) + 20.0_dp * y(1)
                       grad(2) = grad(2) + 20.0_dp * y(2)
                       IF (y(2) <= 0.0_dp) THEN
                          grad(2) = grad(2) - 10.0_dp
                       ELSE
                          grad(2) = grad(2) + 10.0_dp
                       END IF
                     END IF
                
                     IF (ind == 2) THEN
                       grad(1) = grad(1) + 10.0_dp + 20.0_dp * y(1)
                       grad(2) = grad(2) + 20.0_dp * y(2)
                       IF (y(2) <= 0.0_dp) THEN
                          grad(2) = grad(2) - 10.0_dp
                       ELSE
                          grad(2) = grad(2) + 10.0_dp
                       END IF
                     END IF 
                
                     IF (ind == 3) THEN
                        IF ( ( y(1) - y(2) ) <= 0.0_dp ) THEN
                          grad(1) = grad(1) - 10.0_dp 
                          grad(2) = grad(2) + 10.0_dp
                        ELSE
                          grad(1) = grad(1) + 10.0_dp
                          grad(2) = grad(2) - 10.0_dp                   
                        END IF
                    
                        IF (y(2) <= 0.0_dp) THEN
                          grad(2) = grad(2) - 10.0_dp                   
                        ELSE
                          grad(2) = grad(2) +   10.0_dp                 
                        END IF
                     END IF 
                
                     IF (ind == 4) THEN
                       grad(1) = grad(1) + 10.0_dp + 20.0_dp * y(1)
                       grad(2) = grad(2) + 20.0_dp * y(2)
                     END IF
                
                     IF ( (ABS(y(1)) - y(2) ) >= 0.0_dp ) THEN
                       grad(2) = grad(2) - 200.0_dp
                       IF ( y(1) <= 0.0_dp ) THEN
                         grad(1) = grad(1) - 200.0_dp
                       ELSE
                         grad(1) = grad(1) + 200.0_dp
                       END IF 
                     END IF
                     
                     grad = grad * factor                    
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   8
                   !-------------------------------------                      
                   CASE(8)
                     grad(1) = -8.0_dp + 8.0_dp * y(1)
                     grad(2) = -6.0_dp + 4.0_dp * y(2)
                     grad(3) = -4.0_dp + 4.0_dp * y(3)

                     IF ( y(1) <= 0.0_dp ) THEN 
                         grad(1) = grad(1) - 2.0_dp
                     ELSE
                         grad(1) = grad(1) + 2.0_dp
                     END IF                 
                
                     IF ( y(2) <= 0.0_dp ) THEN 
                         grad(2) = grad(2) - 2.0_dp
                     ELSE
                         grad(2) = grad(2) + 2.0_dp
                     END IF 

                     IF ( y(3) <= 0.0_dp ) THEN 
                         grad(3) = grad(3) - 2.0_dp
                     ELSE
                         grad(3) = grad(3) + 2.0_dp
                     END IF                 
                
                     b(1) = y(1) + y(2) + 2.0_dp * y(3) - 3.0_dp
                     b(2) = -y(1)
                     b(3) = -y(2)
                     b(4) = -y(3)
                     b(5) = 0.0_dp
                
                     ind = 5
                
                     DO i = 1, 4
                       IF ( b(ind) < b(i) ) THEN 
                         ind = i
                       END IF 
                     END DO
                
                     IF (ind == 1) THEN
                       grad(1) = grad(1) + 10.0_dp 
                       grad(2) = grad(2) + 10.0_dp 
                       grad(3) = grad(3) + 20.0_dp 
                     END IF
                
                     IF (ind == 2) THEN
                       grad(1) = grad(1) - 10.0_dp 
                     END IF 
                
                     IF (ind == 3) THEN
                       grad(2) = grad(2) - 10.0_dp 
                     END IF 
                
                     IF (ind == 4) THEN
                       grad(3) = grad(3) - 10.0_dp 
                     END IF
                     
                     grad = grad * factor                    
                   !-------------------------------------

                   !-------------------------------------
                   !           Problem   9
                   !-------------------------------------                       
                   CASE(9)
                     grad(1) = 10.0_dp * y(1) - 16.0_dp
                     grad(2) = 10.0_dp * y(2) - 10.0_dp 
                     grad(3) = 10.0_dp * y(3) - 16.0_dp
                     grad(4) = 10.0_dp * y(4) - 10.0_dp  

                     grad = grad * factor                    
                   !------------------------------------
                    
                   !-------------------------------------
                   !           Problem   10
                   !-------------------------------------  
                   CASE(10)
                      grad = 2.0_dp * y    
                      grad = grad * factor
                   !------------------------------------- 

                   !-------------------------------------
                   !           Problem   11
                   !-------------------------------------
                   CASE(11)
                     grad = 0.0_dp
                     
                     IF (y(1)> 0.0_dp) THEN
                        grad(1) = grad(1) + 4.0_dp
                     ELSE
                        grad(1) = grad(1) -4.0_dp
                     END IF
                     
                     IF (y(2)> 0.0_dp) THEN
                        grad(2) = grad(2) + 2.0_dp
                     ELSE
                        grad(2) = grad(2) -2.0_dp
                     END IF

                     IF (y(3)> 0.0_dp) THEN
                        grad(3) = grad(3) + 2.0_dp
                     ELSE
                        grad(3) = grad(3) -2.0_dp
                     END IF                  
                     
                     grad(1) = grad(1) -33.0_dp 
                     grad(2) = grad(2) +16.0_dp 
                     grad(3) = grad(3) -24.0_dp 

                     apu = 2.0_dp*ABS(y(2))-3.0_dp*y(1)-7.0_dp
                     IF (apu > 0.0_dp) THEN 
                        grad(1) = grad(1) -300.0_dp
                        IF (y(2)>0.0_dp) THEN 
                           grad(2) = grad(2) +200.0_dp
                        ELSE
                           grad(2) = grad(2) -200.0_dp
                        END IF 
                     END IF
                     
                     apu = ABS(y(3))-4.0_dp*y(1)-11.0_dp
                     IF (apu > 0.0_dp) THEN 
                        grad(1) = grad(1) -400.0_dp
                        IF (y(3)>0.0_dp) THEN 
                           grad(3) = grad(3) +100.0_dp
                        ELSE
                           grad(3) = grad(3) -100.0_dp
                        END IF 
                     END IF  

                     grad = grad * factor                    
                   !-------------------------------------   

                   !-------------------------------------
                   !           Problem   12
                   !-------------------------------------
                   CASE(12)                
                     grad = 0.0_dp
                     
                     DO i = 1, user_n
                        IF (y(i) > 0.0_dp) THEN 
                           grad(i) = grad(i) +1.0_dp
                        ELSE
                           grad(i) = grad(i) -1.0_dp
                        END IF               
                     END DO                  

                     DO i = 1, user_n
                         apu = 2.0_dp * (y(i)*y(i)- y(i)-1.0_dp) 
                         IF (apu > 0.0_dp) THEN 
                            grad(i) = grad(i) + 40.0_dp * y(i) -20.0_dp
                         END IF
                     END DO
                     
                     grad = grad * factor
                   !-------------------------------------   
                   
                   !-------------------------------------
                   !           Problem   13
                   !-------------------------------------
                   CASE(13)             
                     grad = 0.0_dp
                     
                     DO i = 1, user_n-1
                        IF ( (y(i)+y(i+1)) > 0.0_dp) THEN 
                           grad(i) = grad(i) + 1.0_dp
                           grad(i+1) = grad(i+1) + 1.0_dp
                        ELSE
                           grad(i) = grad(i) - 1.0_dp
                           grad(i+1) = grad(i+1) - 1.0_dp
                        END IF
                     END DO
                     
                     DO i = 1, user_n-2
                        IF ( (y(i)+y(i+2)) > 0.0_dp) THEN 
                           grad(i) = grad(i) + 1.0_dp
                           grad(i+2) = grad(i+2) + 1.0_dp
                        ELSE
                           grad(i) = grad(i) - 1.0_dp
                           grad(i+2) = grad(i+2) - 1.0_dp
                        END IF
                     END DO
                     
                     IF ( (y(1)+y(9)) > 0.0_dp) THEN 
                        grad(1) = grad(1) + 1.0_dp
                        grad(9) = grad(9) + 1.0_dp
                     ELSE
                        grad(1) = grad(1) - 1.0_dp
                        grad(9) = grad(9) - 1.0_dp
                     END IF 
                     
                     IF ( (y(1)+y(10)) > 0.0_dp) THEN 
                        grad(1) = grad(1) + 1.0_dp
                        grad(10) = grad(10) + 1.0_dp
                     ELSE
                        grad(1) = grad(1) - 1.0_dp
                        grad(10) = grad(10) - 1.0_dp
                     END IF 
                     
                     IF ( (y(2)+y(10)) > 0.0_dp) THEN 
                        grad(2) = grad(2) + 1.0_dp
                        grad(10) = grad(10) + 1.0_dp
                     ELSE
                        grad(2) = grad(2) - 1.0_dp
                        grad(10) = grad(10) - 1.0_dp
                     END IF 
                     
                     IF ( (y(1)+y(5)) > 0.0_dp) THEN 
                        grad(1) = grad(1) + 1.0_dp
                        grad(5) = grad(5) + 1.0_dp
                     ELSE
                        grad(1) = grad(1) - 1.0_dp
                        grad(5) = grad(5) - 1.0_dp
                     END IF 
                     
                     IF ( (y(4)+y(7)) > 0.0_dp) THEN 
                        grad(4) = grad(4) + 1.0_dp
                        grad(7) = grad(7) + 1.0_dp
                     ELSE
                        grad(4) = grad(4) - 1.0_dp
                        grad(7) = grad(7) - 1.0_dp
                     END IF                      
                 
                     apu = 0.0_dp
                     DO i = 1, user_n
                        apu = apu + y(i)
                     END DO
                     
                     IF ((apu-1.0_dp)>0.0_dp) THEN 
                        grad = grad + 10.0_dp 
                     END IF
                     
                     DO i = 1, user_n
                        IF (-y(i) > 0.0_dp) THEN 
                            grad(i) = grad(i) - 10.0_dp
                        END IF
                     END DO     

                     grad = grad * factor                    
                   !-------------------------------------   

                   !-------------------------------------
                   !           Problem   14
                   !-------------------------------------
                   CASE(14)         
                     grad = 0.0_dp
                   
                     apu = 0.0_dp
                     DO j =1, user_n
                        apu = apu + y(j)/(1+j-1.0_dp)
                     END DO
                     IF (apu >= 0.0_dp) THEN 
                           abs_sign(1) = 1.0_dp
                     ELSE   
                           abs_sign(1) = -1.0_dp
                     END IF        
                     largest = ABS(apu) 
                     ind = 1
                     
                     DO i = 2, user_n
                        apu = 0.0_dp
                        DO j = 1, user_n
                          apu = apu + y(j)/(i+j-1.0_dp)
                        END DO
                        IF (apu >= 0.0_dp) THEN 
                           abs_sign(i) = 1.0_dp
                        ELSE   
                           abs_sign(i) = -1.0_dp
                        END IF
                        IF (ABS(apu) > largest) THEN 
                            largest = ABS(apu)
                            ind = i
                        END IF
                     END DO
                     
                     DO j = 1, user_n
                        grad(j) = abs_sign(ind) * user_n / (ind+j-1)
                     END DO                      
    
                     grad = grad * factor   
                   !-------------------------------------   

                   !-------------------------------------
                   !           Problem   15
                   !-------------------------------------
                   CASE(15)     
                     grad = 0.0_dp
                     
                     DO i = 1, user_n-1
                         a(1) = y(i)**4 + y(i+1)**2
                         a(2) = (2.0_dp-y(i))**2 + (2.0_dp-y(i+1))**2
                         a(3) = 2.0_dp*EXP(-y(i)+y(i+1))
                         IF (a(1)>a(2)) THEN
                            IF (a(1)>a(3)) THEN
                              f_i(i) = a(1)
                            ELSE
                              f_i(i) = a(3)
                            END IF
                         ELSE
                            IF (a(2)>a(3)) THEN
                              f_i(i) = a(2)
                            ELSE
                              f_i(i) = a(3)
                            END IF      
                         END IF
                     END DO
                     
                     largest = f_i(1)
                     ind = 1
                     DO i =2, user_n-1
                         IF (f_i(i) > largest) THEN 
                            largest = f_i(i)
                            ind = i
                         END IF
                     END DO
                     
                         a(1) = y(ind)**4 + y(ind+1)**2
                         a(2) = (2.0_dp-y(ind))**2 + (2.0_dp-y(ind+1))**2
                         a(3) = 2.0_dp*EXP(-y(ind)+y(ind+1))
                         IF (a(1)>a(2)) THEN
                            IF (a(1)>a(3)) THEN
                              grad(ind) = 4.0_dp*y(ind)**3 
                              grad(ind+1) = 2.0_dp*y(ind+1)
                            ELSE
                              grad(ind) = -2.0_dp * EXP(-y(ind)+y(ind+1)) 
                              grad(ind+1) = 2.0_dp * EXP(-y(ind)+y(ind+1)) 
                            END IF
                         ELSE
                            IF (a(2)>a(3)) THEN
                              grad(ind) = -2.0_dp * (2.0_dp -y(ind))
                              grad(ind+1) = -2.0_dp * (2.0_dp -y(ind+1))
                            ELSE
                              grad(ind) = -2.0_dp * EXP(-y(ind)+y(ind+1)) 
                              grad(ind+1) = 2.0_dp * EXP(-y(ind)+y(ind+1)) 
                            END IF      
                         END IF                  
                     
                     grad = (user_n-1.0_dp)*grad

                     grad = grad * factor
                     
                   !-------------------------------------   
                       
                   
                   !-------------------------------------
                   !             Problem 16:
                   !         Chained Crescent I
                   !-------------------------------------
                   CASE(16)
                     grad = 0.0_dp
                     
                     apu = 0.0_dp
                     DO i = 1, user_n-1
                        apu = apu + y(i)**2 + (y(i+1)-1.0_dp)**2
                        apu = apu + y(i+1) -1.0_dp
                     END DO
                     
                     apu = 2.0_dp *apu
                     IF (apu > 0.0_dp) THEN
                       DO i = 1, user_n-1
                          grad(i) = grad(i) + 4.0_dp *y(i)
                          grad(i+1) = grad(i+1) + 4.0_dp *(y(i+1)-1.0_dp) +2.0_dp
                       END DO
                     END IF
                    
                     grad = grad * factor
                     
                   !-------------------------------------
                      
 
                END SELECT                      
                
           END FUNCTION subgradient_f1      
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
           
           FUNCTION subgradient_f2(y, problem2, factor) RESULT(grad)                
                !
                ! Calculate a subgradient of the DC component f_2 at a point 'y'.
                ! Variable 'problem1' identifies the objective function used.
                ! 'factor' is a scaling coefficient for the objective
                !
                ! NOTICE: * The dimension of 'y' has to be 'user_n'.
                !         * The dimension of 'grad' is also 'user_n'.
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: y    ! a point where the subgradient of the DC component f_2 is calculated
                INTEGER, INTENT(IN) :: problem2                 ! the objective function selected for the DC component f_2
                REAL(KIND=dp), INTENT(IN) :: factor             ! the scaling factor for objective              
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp), DIMENSION(SIZE(y)) :: grad       ! the subgradient of the DC component f_2 at a point 'y'
                REAL(KIND=dp), PARAMETER :: zero=0.0_dp, one=1.0_dp
                REAL(KIND=dp), PARAMETER :: two=2.0_dp
                REAL(KIND=dp) :: a,b 
                REAL(KIND=dp) :: a4, a5, a6 
                REAL(KIND=dp) :: apu                                ! help varaible
                REAL(KIND=dp) :: largest, abs_sign                  ! help varaible
                INTEGER :: i, j, ind                                ! help variable
                
                REAL(KIND=dp), DIMENSION(3) :: aa               ! sign of the active partial objective for f1 and for f2    
                REAL(KIND=dp), DIMENSION(user_n) :: term        ! sign of the active partial objective for f1 and for f2    
                
                
                SELECT CASE(problem2)
                
                  
                   !-------------------------------------
                   !           Problem   1
                   !-------------------------------------                    
                   CASE(1)              
                     grad = 0.0_dp
                
                     a4 = y(1)**2 - 2.0_dp * y(1) + y(2)**2 - 4.0_dp * y(2) + 4.0_dp
                     a5 = 2.0_dp * y(1)**2 - 5.0_dp * y(1) + y(2)**2 - 2.0_dp * y(2) 
                     a5 = a5 + 4.0_dp
                     a6 = y(1)**2 + 2.0_dp * y(2)**2 - 4.0_dp * y(2) + 1.0_dp               
                
                     IF ( (a4 + a5) >= (a4 + a6) ) THEN
                        IF ( (a4 + a5) >= (a5 + a6)) THEN
                          grad(1) = 6.0_dp * y(1) - 7.0_dp
                          grad(2) = 4.0_dp * y(2) - 6.0_dp
                        ELSE
                          grad(1) = 6.0_dp * y(1) - 5.0_dp
                          grad(2) = 6.0_dp * y(2) - 6.0_dp
                        END IF              
                     ELSE
                        IF ( (a4 + a6) >= (a5 + a6) ) THEN
                          grad(1) = 4.0_dp * y(1) - 2.0_dp
                          grad(2) = 6.0_dp * y(2) - 8.0_dp
                        ELSE
                          grad(1) = 6.0_dp * y(1) - 5.0_dp
                          grad(2) = 6.0_dp * y(2) - 6.0_dp 
                        END IF                  
                     END IF
                     
                     grad = grad * factor                    
                   !------------------------------------- 
                   
                   !-------------------------------------
                   !           Problem   2
                   !-------------------------------------                    
                   CASE(2)              
                     IF(y(1) <= 0) THEN
                       grad(1) = -100.0_dp 
                     ELSE
                       grad(1) = 100.0_dp
                     END IF
                
                     grad(2) = -100.0_dp
                     
                     grad = grad * factor                    
                   !-------------------------------------                   
                   
                   !-------------------------------------
                   !           Problem   3
                   !-------------------------------------  
                   CASE(3)
                     IF(y(1) <= 0.0_dp) THEN
                       grad(1) = -100.0_dp 
                     ELSE
                       grad(1) = 100.0_dp
                     END IF
                
                     grad(2) = -100.0_dp
                
                     IF(y(3) <= 0.0_dp) THEN
                       grad(3) = -90.0_dp 
                     ELSE
                       grad(3) = 90.0_dp
                     END IF
                
                     grad(4) = -90.0_dp             
                
                     IF( (y(2) - y(4) ) <= 0.0_dp) THEN
                       grad(2) = grad(2) - 4.95_dp 
                       grad(4) = grad(4) + 4.95_dp
                     ELSE
                       grad(2) = grad(2) + 4.95_dp 
                       grad(4) = grad(4) - 4.95_dp                 
                     END IF         

                     grad = grad * factor                    
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   4
                   !-------------------------------------  
                   CASE(4)
                     DO i = 1, user_n
                       IF (y(i) <= 0.0_dp ) THEN 
                         grad(i) = -1.0_dp
                       ELSE
                         grad(i) = 1.0_dp
                       END IF
                     END DO
                     
                     grad = grad * factor
                   !-------------------------------------
                    
                   !-------------------------------------
                   !           Problem   5
                   !-------------------------------------  
                   CASE(5)
                     grad = 0.0_dp 

                     DO j = 1, 20
                      IF (summa(y,j) <= 0.0_dp) THEN 
                        DO i = 1, user_n
                            grad(i) = grad(i) - (0.05_dp * j)**(i-1)
                        END DO                  
                      ELSE 
                        DO i = 1, user_n
                            grad(i) = grad(i) + (0.05_dp * j)**(i-1)
                        END DO                   
                      END IF
                     END DO   

                     grad = grad * factor                    
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   6
                   !-------------------------------------  
                   CASE(6)
                     DO i = 1, user_n
                       IF (y(i) <= 0.0_dp ) THEN 
                         grad(i) = -1.0_dp
                       ELSE
                         grad(i) = 1.0_dp
                       END IF
                     END DO 

                     grad = grad * factor                    
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   7
                   !-------------------------------------  
                   CASE(7)
                     grad = 0.0_dp
                
                     IF (y(2) <= 0.0_dp ) THEN 
                         grad(2) = 20.0_dp * y(2) - 10.0_dp 
                     ELSE
                         grad(2) = 20.0_dp * y(2) + 10.0_dp 
                     END IF
                
                     grad(2) = grad(2) - 100.0_dp
                     grad(1) = 20.0_dp * y(1)

                     IF (y(1) <= 0.0_dp ) THEN 
                         grad(1) = grad(1) - 100.0_dp
                     ELSE
                         grad(1) = grad(1) + 100.0_dp
                     END IF
                     
                     grad = grad * factor                    
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   8
                   !-------------------------------------                    
                   CASE(8)
                     grad = 0.0_dp
                
                     IF ( (y(1)-y(2)) <= 0.0_dp ) THEN 
                         grad(1) = grad(1) - 1.0_dp 
                         grad(2) = grad(2) + 1.0_dp
                     ELSE
                         grad(1) = grad(1) + 1.0_dp
                         grad(2) = grad(2) - 1.0_dp 
                     END IF
                
                     IF ( (y(1)-y(3)) <= 0.0_dp ) THEN 
                         grad(1) = grad(1) - 1.0_dp 
                         grad(3) = grad(3) + 1.0_dp
                     ELSE
                         grad(1) = grad(1) + 1.0_dp
                         grad(3) = grad(3) - 1.0_dp     
                     END IF 
                     
                     grad = grad * factor                    
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   9
                   !-------------------------------------                     
                   CASE(9)
                     grad = 0.0_dp 

                     a = (y(1) -2.0_dp)**2 + y(2)**2
                     b = (y(3) -2.0_dp)**2 + y(4)**2
                
                     IF ( a >= b) THEN 
                       grad(1) = grad(1) + 2.0_dp * (y(1)-2.0_dp)
                       grad(2) = grad(2) + 2.0_dp * y(2)
                     ELSE 
                       grad(3) = grad(3) + 2.0_dp * (y(3)-2.0_dp)
                       grad(4) = grad(4) + 2.0_dp * y(4)
                     END IF

                     a = (y(1) -2.0_dp)**2 + (y(2)-1.0_dp)**2
                     b = (y(3) -2.0_dp)**2 + (y(4)-1.0_dp)**2
                
                     IF ( a >= b) THEN 
                       grad(1) = grad(1) + 2.0_dp * (y(1)-2.0_dp)
                       grad(2) = grad(2) + 2.0_dp * (y(2)-1.0_dp)
                     ELSE 
                       grad(3) = grad(3) + 2.0_dp * (y(3)-2.0_dp)
                       grad(4) = grad(4) + 2.0_dp * (y(4)-1.0_dp)
                     END IF

                     a = (y(1) -3.0_dp)**2 + y(2)**2
                     b = (y(3) -3.0_dp)**2 + y(4)**2
                
                     IF ( a >= b) THEN 
                       grad(1) = grad(1) + 2.0_dp * (y(1)-3.0_dp)
                       grad(2) = grad(2) + 2.0_dp * (y(2))
                     ELSE 
                       grad(3) = grad(3) + 2.0_dp * (y(3)-3.0_dp)
                       grad(4) = grad(4) + 2.0_dp * (y(4)) 
                     END IF

                     a = (y(1))**2 + (y(2)-2.0_dp)**2
                     b = (y(3))**2 + (y(4)-2.0_dp)**2
                
                     IF ( a >= b) THEN 
                       grad(1) = grad(1) + 2.0_dp * (y(1))
                       grad(2) = grad(2) + 2.0_dp * (y(2)-2.0_dp)
                     ELSE 
                       grad(3) = grad(3) + 2.0_dp * (y(3))
                       grad(4) = grad(4) + 2.0_dp * (y(4)-2.0_dp)
                     END IF

                     a = (y(1)-1.0_dp)**2 + (y(2)-2.0_dp)**2
                     b = (y(3)-1.0_dp)**2 + (y(4)-2.0_dp)**2
                
                     IF ( a >= b) THEN 
                       grad(1) = grad(1) + 2.0_dp * (y(1)-1.0_dp)
                       grad(2) = grad(2) + 2.0_dp * (y(2)-2.0_dp)
                     ELSE 
                       grad(3) = grad(3) + 2.0_dp * (y(3)-1.0_dp)
                       grad(4) = grad(4) + 2.0_dp * (y(4)-2.0_dp)
                     END IF           

                     grad = grad * factor                    
                   !------------------------------------
                    
                   !-------------------------------------
                   !           Problem   10
                   !-------------------------------------  
                   CASE(10)
                     grad = 0.0_dp
                
                     DO i = 2, user_n
                       IF ( y(i) - y(i-1) <= 0.0_dp) THEN
                         grad(i-1) = grad(i-1) + 1.0_dp 
                         grad(i) = grad(i) - 1.0_dp 
                       ELSE
                         grad(i-1) = grad(i-1) - 1.0_dp  
                         grad(i) = grad(i) + 1.0_dp     
                       END IF
                     END DO
                     
                     grad = grad * factor
                   !-------------------------------------
                   
                   !-------------------------------------
                   !           Problem   11
                   !-------------------------------------  
                   CASE(11)
                     grad = 0.0_dp
                
                     grad(1) = -140.0_dp
                     
                     IF (y(2)>0.0_dp) THEN 
                        grad(2) = 40.0_dp
                     ELSE
                        grad(2) = -40.0_dp
                     END IF
                     
                     IF (y(3)>0.0_dp) THEN 
                        grad(3) = -20.0_dp
                     ELSE
                        grad(3) = 20.0_dp
                     END IF       

                     grad = grad * factor                    
                   !-------------------------------------       

                   !-------------------------------------
                   !           Problem   12
                   !-------------------------------------
                   CASE(12)                
                     grad = 0.0_dp

                     DO i = 1, user_n
                        grad(i) = grad(i) +20.0_dp*y(i) - 10.0_dp
                     END DO
                     
                     term = 0.0_dp
                     DO i = 1, user_n
                        DO j = 1, user_n 
                           IF (i /= j) THEN
                               term(i) = term(i) + ABS(y(j))
                           END IF
                        END DO
                     END DO
                     
                     largest = term(1)
                     ind = 1 
                     
                     DO i = 2, user_n
                        IF (term(i) > largest) THEN 
                           largest = term(i)
                           ind = i
                        END IF
                     END DO
                     
                     DO i = 1, user_n
                        IF (i /= ind) THEN
                           IF (y(i) > 0.0_dp) THEN 
                              grad(i) = grad(i) +1.0_dp
                           ELSE
                              grad(i) = grad(i) -1.0_dp
                           END IF
                        END IF                   
                     END DO
                     
                     grad = grad * factor                    
                   !-------------------------------------   
                   
                   !-------------------------------------
                   !           Problem   13
                   !-------------------------------------
                   CASE(13)             
                     grad = 0.0_dp
                     
                     DO i = 1, user_n-1
                        IF ( y(i) > 0.0_dp) THEN
                           grad(i) = grad(i) + 1.0_dp
                        ELSE
                           grad(i) = grad(i) - 1.0_dp
                        END IF
                        IF ( y(i+1) > 0.0_dp) THEN
                           grad(i+1) = grad(i+1) + 1.0_dp
                        ELSE
                           grad(i+1) = grad(i+1) - 1.0_dp
                        END IF                      
                     END DO
                     
                     DO i = 1, user_n-2
                        IF ( y(i) > 0.0_dp) THEN
                           grad(i) = grad(i) + 1.0_dp
                        ELSE
                           grad(i) = grad(i) - 1.0_dp
                        END IF
                        IF ( y(i+2) > 0.0_dp) THEN
                           grad(i+2) = grad(i+2) + 1.0_dp
                        ELSE
                           grad(i+2) = grad(i+2) - 1.0_dp
                        END IF      
                     END DO                  
                     
                    
                     IF ( y(1) > 0.0_dp) THEN
                        grad(1) = grad(1) + 3.0_dp
                     ELSE
                        grad(1) = grad(1) - 3.0_dp
                     END IF

                     IF ( y(2) > 0.0_dp) THEN
                        grad(2) = grad(2) + 1.0_dp
                     ELSE
                        grad(2) = grad(2) - 1.0_dp
                     END IF
                     
                     IF ( y(4) > 0.0_dp) THEN
                        grad(4) = grad(4) + 1.0_dp
                     ELSE
                        grad(4) = grad(4) - 1.0_dp
                     END IF 

                     IF ( y(5) > 0.0_dp) THEN
                        grad(5) = grad(5) + 1.0_dp
                     ELSE
                        grad(5) = grad(5) - 1.0_dp
                     END IF

                     IF ( y(7) > 0.0_dp) THEN
                        grad(7) = grad(7) + 1.0_dp
                     ELSE
                        grad(7) = grad(7) - 1.0_dp
                     END IF

                     IF ( y(9) > 0.0_dp) THEN
                        grad(9) = grad(9) + 1.0_dp
                     ELSE
                        grad(9) = grad(9) - 1.0_dp
                     END IF

                     IF ( y(10) > 0.0_dp) THEN
                        grad(10) = grad(10) + 2.0_dp
                     ELSE
                        grad(10) = grad(10) - 2.0_dp
                     END IF     

                     grad = grad * factor                    
                   !-------------------------------------   

                   !-------------------------------------
                   !           Problem   14
                   !-------------------------------------
                   CASE(14)         
                     grad = 0.0_dp
                     
                     DO i = 1, user_n
                        apu = 0.0_dp
                        DO j = 1, user_n
                          apu = apu + y(j)/(i+j-1.0_dp)
                        END DO
                        IF (apu >= 0.0_dp) THEN 
                           abs_sign = 1.0_dp
                        ELSE   
                           abs_sign = -1.0_dp
                        END IF
                        DO j = 1, user_n
                          grad(j) = grad(j) + abs_sign / (i+j-1.0_dp)
                        END DO
                     END DO  

                     grad = grad * factor                    
                   !-------------------------------------   

                   !-------------------------------------
                   !           Problem   15
                   !-------------------------------------
                   CASE(15)     
                     grad = 0.0_dp
                                 
                     DO ind = 1, user_n-1
                         aa(1) = y(ind)**4 + y(ind+1)**2
                         aa(2) = (2.0_dp-y(ind))**2 + (2.0_dp-y(ind+1))**2
                         aa(3) = 2.0_dp*EXP(-y(ind)+y(ind+1))
                         IF (aa(1)>aa(2)) THEN
                            IF (aa(1)>aa(3)) THEN
                              grad(ind) = grad(ind) + 4.0_dp*y(ind)**3 
                              grad(ind+1) = grad(ind+1) + 2.0_dp*y(ind+1)
                            ELSE
                              grad(ind) = grad(ind) -2.0_dp * EXP(-y(ind)+y(ind+1)) 
                              grad(ind+1) = grad(ind+1) + 2.0_dp * EXP(-y(ind)+y(ind+1)) 
                            END IF
                         ELSE
                            IF (aa(2)>aa(3)) THEN
                              grad(ind) = grad(ind) - 2.0_dp * (2.0_dp -y(ind))
                              grad(ind+1) = grad(ind+1) - 2.0_dp * (2.0_dp -y(ind+1))
                            ELSE
                              grad(ind) = grad(ind) -2.0_dp * EXP(-y(ind)+y(ind+1)) 
                              grad(ind+1) = grad(ind+1) +2.0_dp * EXP(-y(ind)+y(ind+1)) 
                            END IF      
                         END IF 
                     END DO  
                                      
                     grad = grad * factor
                     
                   !-------------------------------------                      
                   
                   !-------------------------------------
                   !              Problem 16:
                   !         Chained Crescent I
                   !-------------------------------------
                   CASE(16)
                     grad = 0.0_dp

                     DO i = 1, user_n-1
                        grad(i) = grad(i) + 2.0_dp * y(i) 
                        grad(i+1) = grad(i+1) + 2.0_dp * (y(i+1) -1.0_dp) +1.0_dp
                     END DO

                     grad = grad * factor
                     
                   !-------------------------------------       
                   
                
                END SELECT  

           END FUNCTION subgradient_f2
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -      
           
           
           FUNCTION summa(y, j) RESULT(f)           
                !
                ! Calculates the sum used in f_1 and f_2 for parameter t_j
                !
                ! NOTICE: The dimension of 'y' has to be 'n'. 'j' needs to be an integer from interval [1,20]
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: y    ! a point where the function value of the sum is calculated
                INTEGER, INTENT(IN) ::  j                       ! used to determines the parameter t_j
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp) :: f, t, apu                      ! the function value of the sum, the parameter t_j and a help variable
                INTEGER :: i                                    ! help variable
                
                f = 0.0_dp
                apu = 1.0_dp/ user_n
                
                DO i = 1, user_n 
                    t = (0.05_dp * j )**(i-1)
                    f = f + (y(i)-apu)*t            
                END DO 
            
           END FUNCTION summa

           
        
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !| |                                                                          | |
        !| |                       CONSTRAINT FUNCTIONS:                              | |
        !| |                                                                          | |
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*   

        
           
        !********************************************************************************
        !                                                                               |
        !              FUNCTION VALUES OF THE DC COMPONENTS g_1 AND g_2                 |
        !                                                                               |
        !********************************************************************************

           FUNCTION g1(y, problem, factor) RESULT(g)      
                !
                ! For the constraint calculates the function value of the DC component g_1 at a point 'y'.
                ! Variable 'problem' identifies the constraint function used.
                ! 'factor' is a scaling coefficient for the objective
                !
                ! NOTICE: * The dimension of 'y' has to be 'user_n'.
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: y    ! a point where the function value of the DC component g_1 is calculated
                REAL(KIND=dp), INTENT(IN) :: factor             ! a scaling factor
                INTEGER, INTENT(IN) :: problem                  ! identifies the constraint 
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp) :: g       ! the function value of the DC component g_1 at a point 'y'
                REAL(KIND=dp) :: a1     
               
                g = 0.0_dp
                SELECT CASE(problem)
                   
                    
                   !-------------------------------------
                   !           Constraint   1
                   !-------------------------------------                    
                   CASE(1)              
                       
                      a1 = (y(1)+1.5_dp)**2+y(2)**2-4.0_dp
                      a1 = a1+(y(1)-1.0_dp)**2+(y(2)-1.0_dp)**2-1.0_dp
                      IF(a1>=0)THEN
                        g = a1
                      ELSE 
                        g = 0.0_dp
                      END IF
                      g = g*factor
                    
                   !------------------------------------- 
                   
                   !-------------------------------------
                   !           Constraint   2
                   !-------------------------------------                    
                   CASE(2)
                      g = 0.0_dp

                   !-------------------------------------                   
                   
                   !-------------------------------------
                   !           Constraint   3
                   !-------------------------------------  
                   CASE(3)
                   
                       g = user_n * 0.5 

                   !-------------------------------------

                
                END SELECT              
            


           END FUNCTION g1
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
           
           FUNCTION g2(y, problem, factor) RESULT(g)      
                !
                ! For the constraint calculates the function value of the DC component g_2 at a point 'y'.
                ! Variable 'problem' identifies the constraint function used.
                ! 'factor' is a scaling coefficient for the objective           
                !
                ! NOTICE: * The dimension of 'y' has to be 'user_n'.
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: y    ! a point where the function value of the DC component g_2 is calculated
                REAL(KIND=dp), INTENT(IN) :: factor             ! a scaling factor
                INTEGER, INTENT(IN) :: problem                  ! identifies the constraint 
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp) :: g       ! the function value of the DC component g_2 at a point 'y'
                REAL(KIND=dp) :: a1, a2    
                INTEGER :: i
                
                g = 0.0_dp
                
                
                SELECT CASE(problem)

                   
                   !-------------------------------------
                   !           Constraint   1
                   !-------------------------------------                    
                   CASE(1)
                   
                      g = (y(1)-1.0_dp)**2+(y(2)-1.0_dp)**2-1.0_dp
                      g = g*factor
                   
                   !------------------------------------- 
                   
                   !-------------------------------------
                   !           Constraint   2
                   !-------------------------------------                    
                    CASE(2)
                    
                      a1 = y(1)+y(2)+y(3)+y(4)-5.5_dp
                      a2 = y(1)**2+y(2)**2+y(3)**2+y(4)**2-10.0_dp
                      
                      IF (a1>=a2) THEN
                         g = a1
                      ELSE   
                         g = a2
                      END IF
                      g = g * factor                  
                      
                   !-------------------------------------                   
                   
                   !-------------------------------------
                   !           Constraint   3
                   !-------------------------------------  
                   CASE(3)
                   
                     g = 0.0_dp                 
                     
                     DO i = 2, user_n, 2
                       g = g + (y(i)-0.5_dp)**2
                     END DO
                     DO i = 1, user_n, 2
                       g = g + (y(i)+0.5_dp)**2
                     END DO
                     g = g * factor 
                   
                   !-------------------------------------
             
                
                END SELECT              
            


           END FUNCTION g2
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
           

        !********************************************************************************
        !                                                                               |
        !                SUBGRADIENTS OF THE DC COMPONENTS g_1 AND g_2                  |
        !                                                                               |       
        !********************************************************************************       
        
           FUNCTION subgradient_g1(y, problem, factor) RESULT(grad)                
                !
                ! For the constraint calculates a subgradient of the DC component g_1 at a point 'y'.
                ! Variable 'problem' identifies the constraint function used.
                ! 'factor' is a scaling coefficient for the objective               
                !
                ! NOTICE: * The dimension of 'y' has to be 'user_n'.
                !         * The dimension of 'grad' is also 'user_n'.
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: y    ! a point where the subgradient of the DC component g_1 is calculated
                REAL(KIND=dp), INTENT(IN) :: factor             ! a scaling factor
                INTEGER, INTENT(IN) :: problem                  ! identifies the constraint 
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp), DIMENSION(SIZE(y)) :: grad       ! the subgradient of the DC component g_1 at a point 'y'
                REAL(KIND=dp) :: a1
                 
                grad = 0.0_dp
 
                SELECT CASE(problem)

                   !-------------------------------------
                   !           Constraint   1
                   !-------------------------------------                    
                   CASE(1)              
                       
                       a1 = (y(1)+1.5_dp)**2+y(2)**2-4.0_dp
                       a1 = a1+(y(1)-1.0_dp)**2+(y(2)-1.0_dp)**2-1.0_dp
                       IF(a1 >= 0)THEN
                          grad(1) = 2.0_dp*(y(1)+1.5_dp)+2.0_dp*(y(1)-1.0_dp)
                          grad(2) = 2.0_dp*y(2)+2.0_dp*(y(2)-1.0_dp)
                       ELSE 
                          grad = 0.0_dp
                       END IF
                       grad = grad*factor                     
                    
                   !------------------------------------- 
                   
                   !-------------------------------------
                   !           Constraint   2
                   !-------------------------------------                    
                   CASE(2)
                       
                       grad = 0.0_dp               
                       
                   !-------------------------------------                   
                   
                   !-------------------------------------
                   !           Constraint   3
                   !-------------------------------------  
                   CASE(3)
                       
                       grad = 0.0_dp
                   
                   !-------------------------------------
                   
                 
                END SELECT  

           END FUNCTION subgradient_g1
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  
           
           FUNCTION subgradient_g2(y, problem, factor) RESULT(grad)                
                !
                ! For the constraint calculates a subgradient of the DC component g_2 at a point 'y'.
                ! Variable 'problem' identifies the constraint function used.
                ! 'factor' is a scaling coefficient for the objective               
                !
                ! NOTICE: * The dimension of 'y' has to be 'user_n'.
                !         * The dimension of 'grad' is also 'user_n'.
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: y         ! a point where the subgradient of the DC component g_2 is calculated
                REAL(KIND=dp), INTENT(IN) :: factor                  ! a scaling factor
                INTEGER, INTENT(IN) :: problem                       ! identifies the constraint 
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp), DIMENSION(SIZE(y)) :: grad       ! the subgradient of the DC component g_2 at a point 'y'
                REAL(KIND=dp) :: a1, a2    
                INTEGER :: i
                
                grad = 0.0_dp
               
                SELECT CASE(problem)

                   !-------------------------------------
                   !           Constraint   1
                   !-------------------------------------                    
                   CASE(1) 
                   
                       grad(1) = 2.0_dp*(y(1)-1.0_dp)
                       grad(2) = 2.0_dp*(y(2)-1.0_dp)
                       grad = factor*grad   
                       
                   !------------------------------------- 
                   
                   !-------------------------------------
                   !           Constraint   2
                   !-------------------------------------                    
                   CASE(2)
                   
                       a1 = y(1)+y(2)+y(3)+y(4)-5.5_dp
                       a2 = y(1)**2+y(2)**2+y(3)**2+y(4)**2-10.0_dp
                       
                       IF (a1 >= a2) THEN
                          grad = 1.0_dp
                       ELSE   
                          DO i = 1, user_n
                             grad(i) = 2.0_dp * y(i)
                          END DO
                       END IF
                       
                       grad = factor * grad                  
                   !-------------------------------------                   
                   
                   !-------------------------------------
                   !           Constraint   3
                   !-------------------------------------  
                   CASE(3)
                       grad = 0.0_dp
                     
                       DO i = 2, user_n, 2
                         grad(i) = grad(i) + 2.0_dp*(y(i)-0.5_dp)
                       END DO
                       DO i = 1, user_n, 2
                         grad(i) = grad(i) + 2.0_dp*(y(i)+0.5_dp)
                       END DO
                       grad = grad * factor 
                   
                   !-------------------------------------
       
       
                END SELECT  

           END FUNCTION subgradient_g2
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -              
           
           
        !----------------------------------------------------------------------------------   
        ! ------------------------------------------------------------------------------- |  
        ! |                                                                             | |
        ! |                            FUCTIONS NEEDED IN                               | |
        ! |                  MULTIOBJECTIVE OPTIMIZATION METHOD MDBDC                   | |
        ! |                                                                             | |
        ! ------------------------------------------------------------------------------- |
        !----------------------------------------------------------------------------------

       

        !********************************************************************************
        !                                                                               |
        !   FUNCTION VALUES OF THE DC COMPONENTS H_1 AND H_2 IN MULTIOBJECTIVE MODEL    |
        !                                                                               |
        !********************************************************************************


           FUNCTION H1(AandB_y_all) RESULT(f)        
                !
                ! Calculates the function value of the DC component H_1=max{A_i,B_l} in multiobjective model when the componets used are 'AandB_y_all' 
                !
                ! NOTICE: * The dimension of 'AandB_y_all' has to be 'number_of_func'.
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: AandB_y_all ! the values of the components in H_1
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp) :: f                              ! the function value of the DC component H_1 at a point 'y'

                INTEGER :: i, ind                               ! help variable
           
                   !-------------------------------------
                     ind = 1 
                     DO i = 2, number_of_func
                        IF (AandB_y_all(i) > AandB_y_all(ind)) THEN 
                           ind = i
                        END IF 
                     END DO
     
                     f = AandB_y_all(ind)
                   !-------------------------------------                  

           END FUNCTION H1   
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
           
           FUNCTION H2(fg2_y_all) RESULT(f)          
                !
                ! Calculates the function value of the DC component H_2 when the DC components f2_i and g2_l have values 'fg2_y_all'.
                !
                ! NOTICE: * The dimension of 'fg2_y_all' has to be 'number_of_func'.
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fg2_y_all    ! function values for all DC components f2_i and g2_l at point y
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp) :: f                                     ! the function value of the DC component H_2 at a point 'y'
                
                INTEGER :: i
               
                    !-------------------------------------             
                      f = 0.0_dp
                      DO i = 1, number_of_func
                          f = f + fg2_y_all(i)
                      END DO
                    !-------------------------------------

           END FUNCTION H2
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
           
           
           
        !********************************************************************************
        !                                                                               |
        !   SUBGRADIENTS OF THE DC COMPONENTS H_1 AND H_2 IN MULTIOBJECTIVE MODEL       |
        !                                                                               |
        !********************************************************************************          
           
            FUNCTION subgradient_H1(grad1_y_all, grad2_y_all, AandB_y_all) RESULT(grad)
                !
                ! Calculates the subgradient of the DC component H_1=max{A_i,B_l} at 'y' in multiobjective model when 
                ! the componets used at 'y' are 'AandB_y_all'
                !
                ! NOTICE: * The dimension of 'AandB_y_all' has to be 'number_of_func'.
                !         * The dimension of 'grad1_y_all' and 'grad2_y_all' has to be '(number_of_func,user_n)'.             
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: AandB_y_all       ! the values of the components A_i in H_1 at y
                REAL(KIND=dp), DIMENSION(:,:), INTENT(IN) :: grad1_y_all     ! all subgradients for the DC components f1_i and g1_l at y             
                REAL(KIND=dp), DIMENSION(:,:), INTENT(IN) :: grad2_y_all     ! all subgradients for the DC components f2_i and g2_l at y                             
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp), DIMENSION(user_n) :: grad        ! the subgradient of the DC component H_1 at a point 'y'
                REAL(KIND=dp), DIMENSION(user_n) :: grad_f2     ! the sum of the subgradients of the DC components f2_i at 'y'
                REAL(KIND=dp), DIMENSION(user_n) :: grad_g2     ! the sum of the subgradients of the DC components g2_i at 'y'  
                
                INTEGER :: i, j, ind                              ! help variable

                   !-------------------------------------
                     ind = 1 
                     DO i = 2, number_of_func
                        IF (AandB_y_all(i) > AandB_y_all(ind)) THEN 
                           ind = i
                        END IF 
                     END DO
                   !-------------------------------------
                    
                   IF (ind < number_of_obj+1) THEN  ! subgradient of A_i is calculated                  
                   !-------------------------------------
                     grad_g2 = 0.0_dp
                     DO i = number_of_obj+1, number_of_func
                       DO j = 1, user_n
                          grad_g2(j) = grad_g2(j) + grad2_y_all(i,j)
                       END DO
                     END DO
                       
                     grad = 0.0_dp
                     DO i = 1, user_n
                        grad(i) = grad1_y_all(ind,i)
                     END DO     

                     DO j = 1, number_of_obj
                        IF ( ind /= j ) THEN 
                           DO i = 1, user_n
                             grad(i) = grad(i) + grad2_y_all(j,i)
                           END DO   
                        END IF
                     END DO
                     
                     grad = grad + grad_g2
                   !-------------------------------------
                   ELSE ! subgradient of B_l is calculated
                   !-------------------------------------
                     grad_f2 = 0.0_dp
                     DO i = 1, number_of_obj
                       DO j = 1, user_n
                          grad_f2(j) = grad_f2(j) + grad2_y_all(i,j)
                       END DO
                     END DO
                       
                     grad = 0.0_dp
                     DO i = 1, user_n
                        grad(i) = grad1_y_all(ind,i)
                     END DO     

                     DO j = number_of_obj+1, number_of_func
                        IF ( ind /= j ) THEN 
                           DO i = 1, user_n
                             grad(i) = grad(i) + grad2_y_all(j,i)
                           END DO   
                        END IF
                     END DO
                     
                     grad = grad + grad_f2
                   !-------------------------------------
                   END IF                    
    
                
           END FUNCTION subgradient_H1  
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
           
           FUNCTION subgradient_H2(grad2_y_all) RESULT(grad)              
                !
                ! Calculates a subgradient of the DC component H_2 at a point 'y' in multiobjective model when 
                ! subgradients of DC componets f2_i and g2_l at y are 'grad2_y_all'.
                !
                ! NOTICE: * The dimension of 'grad_fg2_y_all' has to be '(number_of_func,user_n)'.             
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:,:), INTENT(IN) :: grad2_y_all   ! all subgradients for the DC components f2_i and g2_l at y 
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp), DIMENSION(user_n) :: grad         ! the subgradient of the DC component H_2 at a point 'y'

                INTEGER :: i , j

                   !-------------------------------------
                      grad = 0.0_dp
                      DO i = 1, number_of_func
                        DO j = 1, user_n
                          grad(j) = grad(j) + grad2_y_all(i,j)
                        END DO  
                      END DO                  
                   !------------------------------------- 

           END FUNCTION subgradient_H2
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -          
           
 
           
        !********************************************************************************
        !                                                                               |
        !           GIVES THE INDEX OF THE COMPONENT YIELDING THE VALUE OF H_1          |
        !                                                                               |       
        !********************************************************************************       
        
           FUNCTION index_H1(AandB_y_all) RESULT(ind)
                !
                ! Calculates the index of the component yielding the value of H_1=max{A_i,B_l} when
                ! components used are 'AandB_y_all'
                !
                ! NOTICE: * The dimension of 'AandB_y_all' has to be 'number_of_func'
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: AandB_y_all     ! the values of components used in H_1

                !**************************** OTHER VARIABLES **************************************
                INTEGER :: ind   ! the index of the component yielding the value of H_1                
                INTEGER :: i     ! help variable
              
                     ind = 1 
                     DO i = 2, number_of_func
                        IF (AandB_y_all(i) > AandB_y_all(ind)) THEN 
                           ind = i
                        END IF 
                     END DO
                                
                
           END FUNCTION index_H1                   
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
           
           
        
        !***************************************************************************************
        !                                                                                      |
        !    FUNCTION VALUE AND SUBGRADIENT OF COMPONENT A_i OR B_l IN THE DC COMPONENT H_1    |
        !                                                                                      |
        !***************************************************************************************


           FUNCTION AorB_i(fg1_y_all, fg2_y_all, fg1_k_all, fg2_k_all, ind) RESULT(f)        
                !
                ! Calculates the value of the component A_i or B_l in the DC component H_1 at a point 'y',
                ! when the model is formed at 'x_k'. Variable 'ind' will identify the component function used.
                !
                ! NOTICE: * The dimension of 'fg1_y_all' and 'fg2_y_all' has to be 'number_of_func'.
                !         * The dimension of 'fg1_k_all' and 'fg2_k_all' has to be 'number_of_func'.
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fg1_y_all   ! the values of the DC components f1_i and g1_l at point y
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fg2_y_all   ! the values of the DC components f2_i and g2_l at point y
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fg1_k_all   ! the values of the DC components f1_i and g1_l at point x_k
                REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fg2_k_all   ! the values of the DC components f2_i and g2_l at point x_k
                INTEGER, INTENT(IN) :: ind                             ! index of the component used
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp) :: f                                 ! the value of the component A_'ind' in the DC component H_1 at a point 'y'
                REAL(KIND=dp) :: sum_f2                            ! the sum of values of the all the DC components f2_i at a point 'y'
                REAL(KIND=dp) :: sum_g2                            ! the sum of values of the all the DC components g2_i at a point 'y'

                INTEGER :: i, j                                    ! help variable
                    
                    IF (ind < number_of_obj+1) THEN  ! A_i is calculated 
                   !-------------------------------------
                   
                       sum_g2 = 0.0_dp
                       DO i = number_of_obj+1, number_of_func
                         sum_g2 = sum_g2 + fg2_y_all(i)
                       END DO
                    
                       f = fg1_y_all(ind) 
                       f = f -(fg1_k_all(ind)-fg2_k_all(ind))
     
                       DO j = 1, number_of_obj
                           IF (ind /= j) THEN 
                              f = f + fg2_y_all(j)
                           END IF 
                       END DO    
                       f = f + sum_g2
                   !-------------------------------------  
                    ELSE  ! B_l is calculated
                   !-------------------------------------
                   
                       sum_f2 = 0.0_dp
                       DO i = 1, number_of_obj
                         sum_f2 = sum_f2 + fg2_y_all(i)
                       END DO                  
                       
                       f =  fg1_y_all(ind) 
     
                       DO j = number_of_obj+1, number_of_func
                           IF (ind /= j) THEN 
                              f = f + fg2_y_all(j)
                           END IF 
                       END DO    
                       f = f + sum_f2                   
                   
                   !-------------------------------------
                    END IF                 
                
           END FUNCTION AorB_i 
           
           !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -   
           
           FUNCTION subgradient_AorB_i(grad1_y_all, grad2_y_all, ind) RESULT(grad)
                !
                ! Calculates a subgradient of the component A_i or B_l in the DC component H_1 at a point 'y'.
                ! Variable 'ind' identifies the objective function used.
                !
                ! NOTICE: * The dimension of 'grad1_y_all' and 'grad2_y_all' has to be '(number_of_func,user_n)'.
                !
                IMPLICIT NONE
                !**************************** NEEDED FROM USER *************************************
                REAL(KIND=dp), DIMENSION(:,:), INTENT(IN) :: grad1_y_all    ! subgradients of DC components f1_i and g1_l at point y
                REAL(KIND=dp), DIMENSION(:,:), INTENT(IN) :: grad2_y_all    ! subgradients of DC components f2_i and g2_l at point y
                INTEGER, INTENT(IN) :: ind                                  ! the index of the component used
                !**************************** OTHER VARIABLES **************************************
                REAL(KIND=dp), DIMENSION(user_n) :: grad        ! the subgradient of the component 'ind' in the DC component H_1 at a point 'y'
                REAL(KIND=dp), DIMENSION(user_n) :: grad_f2     ! the sum of the subgradients of the DC components f2_i at 'y'
                REAL(KIND=dp), DIMENSION(user_n) :: grad_g2     ! the sum of the subgradients of the DC components g2_i at 'y'

                INTEGER :: i, j                                 ! help variable

                    
                   IF (ind < number_of_obj+1) THEN  ! subgradient of A_i is calculated                  
                   !-------------------------------------
                     grad_g2 = 0.0_dp
                     DO i = number_of_obj+1, number_of_func
                       DO j = 1, user_n
                          grad_g2(j) = grad_g2(j) + grad2_y_all(i,j)
                       END DO
                     END DO
                       
                     grad = 0.0_dp
                     DO i = 1, user_n
                        grad(i) = grad1_y_all(ind,i)
                     END DO     

                     DO j = 1, number_of_obj
                        IF ( ind /= j ) THEN 
                           DO i = 1, user_n
                             grad(i) = grad(i) + grad2_y_all(j,i)
                           END DO   
                        END IF
                     END DO
                     
                     grad = grad + grad_g2
                   !-------------------------------------
                   ELSE ! subgradient of B_l is calculated
                   !-------------------------------------
                     grad_f2 = 0.0_dp
                     DO i = 1, number_of_obj
                       DO j = 1, user_n
                          grad_f2(j) = grad_f2(j) + grad2_y_all(i,j)
                       END DO
                     END DO
                       
                     grad = 0.0_dp
                     DO i = 1, user_n
                        grad(i) = grad1_y_all(ind,i)
                     END DO     

                     DO j = number_of_obj+1, number_of_func
                        IF ( ind /= j ) THEN 
                           DO i = 1, user_n
                             grad(i) = grad(i) + grad2_y_all(j,i)
                           END DO   
                        END IF
                     END DO
                     
                     grad = grad + grad_f2
                   !-------------------------------------
                   END IF

           END FUNCTION subgradient_AorB_i     

   
                   
           
      END MODULE functions     