!+ Subroutine to determine scaling by conjugate gradients.
!
! Method:
!	The Polak-Ribiere form of the conjugate gradient algorithm
!	is implemented to find a scaling function for the amount
!	of absorber.
!
! Current Owner of Code: J. M. Edwards
!
! History:
!	Version		Date			Comment
!	1.0		12-04-95		First Version under RCS
!						(J. M. Edwards)
!
! Description of Code:
!   FORTRAN 77  with extensions listed in documentation.
!
!- ---------------------------------------------------------------------
      SUBROUTINE CONJUGATE_GRADIENT(IERR, L_RE_ENTRY
     &   , N_SET, LENGTH_SET, U_SET, TRANS_SET
     &   , I_TYPE_RESIDUAL, SCALING, WEIGHT_POINT
     &   , N_TERM, W_ESFT, K_ESFT
     &   , I_SCALE_FUNCTION, SCALE_VECTOR_LONG, P_ARRAY, T_ARRAY)
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/dimensions_spec.ucf'
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/prec_machine.ucf'
      INCLUDE 'incl/precision.pcf'
      INCLUDE 'incl/scale_parameters.acf'
      INCLUDE 'incl/scale_fnc.pcf'
      INCLUDE 'incl/error.pcf'
!
!     DUMMY ARGUMENTS
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      INTEGER	!, INTENT(OUT)
     &     N_SET
!             NUMBER OF SETS OF DATA
     &   , N_TERM
!             NUMBER OF ESFT TERMS
     &   , LENGTH_SET(NPD_SET_SCALE)
!             LENGTHS OF SETS
     &   , I_SCALE_FUNCTION
!             TYPE OF SCALING FUNCTION
     &   , I_TYPE_RESIDUAL
!             TYPE OF RESIDUAL USED
      LOGICAL	!, INTENT(IN)
     &     L_RE_ENTRY
!             FLAG FOR RE-ENTRY
      REAL*8	!, INTENT(IN)
     &     U_SET(NPD_AMOUNT, NPD_SET_SCALE)
!             AMOUNTS OF ABSORBER
     &   , TRANS_SET(NPD_AMOUNT, NPD_SET_SCALE)
!             TRANSMISSIONS
     &   , SCALING(NPD_AMOUNT, NPD_SET_SCALE)
!             SCALINGS REQUIRED
     &   , WEIGHT_POINT(NPD_SET_SCALE, NPD_AMOUNT)
!             WEIGHTS FOR DATA POINTS
     &   , W_ESFT(NPD_ESFT_TERM)
!             ESFT WEIGHTS
     &   , K_ESFT(NPD_ESFT_TERM)
!             ESFT EXPONENTS
     &   , P_ARRAY(NPD_SET_SCALE)
!             PRESSURE ARRAY
     &   , T_ARRAY(NPD_SET_SCALE)
!             TEMPERATURE ARRAY
      REAL*8	!, INTENT(OUT)
     &     SCALE_VECTOR_LONG(NPD_SCALE_VARIABLE, NPD_ESFT_TERM)
!             PARAMETERS FOR SCALING
!
!     LOCAL ARGUMENTS.
      INTEGER
     &     ITERATION
!             ITERATION COUNT
     &   , I
!             LOOP VARIABLE
     &   , J
!             LOOP VARIABLE
      REAL*8
     &     GRADIENT(NPD_SCALE_VARIABLE)
!             GRADIENT OF RESIDUAL
     &   , SCALE_VECTOR(NPD_SCALE_VARIABLE)
!             PARAMETERS FOR SCALING
     &   , NULL_VECTOR(NPD_SCALE_VARIABLE)
!             NULL VECTOR
     &   , RESIDUAL_OLD
!             OLD VALUE OF RESIDUAL
     &   , RESIDUAL_NEW
!             NEW VALUE OF RESIDUAL
     &   , RESIDUAL_TYPICAL
!             TYPICAL VALUE OF RESIDUAL
     &   , G(NPD_SCALE_VARIABLE)
!             REVERSED GRADIENT
     &   , H(NPD_SCALE_VARIABLE)
!             DIRECTION OF SEARCH
     &   , H_L2_NORM
!             L2 NORM OF H
     &   , DGI
!             TEMPORARY CONSTANT
     &   , GI2
!             TEMPORARY CONSTANT
     &   , SEARCH_DISTANCE
!             SEARCHING DISTANCE
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     RESIDUAL_GRADIENT, TERMINATE_SCALE, LINE_SEARCH
!     FUNCTIONS CALLED:
      REAL*8
     &     RESIDUAL_TRANS
!             RESIDUAL FUNCTION
      EXTERNAL
     &     RESIDUAL_TRANS
!
!
      INCLUDE 'incl/scale_fnc.dcf'
!
      DATA NULL_VECTOR/NPD_SCALE_VARIABLE*0.0D+00/
!
!
!
!     INITIALIZE THE SCALING VECTOR.
      DO I=1, N_SCALE_VARIABLE(I_SCALE_FUNCTION)
         SCALE_VECTOR(I)=SCALE_VECTOR_LONG(I, 1)
      ENDDO
!
!     FIND A TYPICAL VALUE OF THE RESIDUAL
      RESIDUAL_TYPICAL=RESIDUAL_TRANS(N_SET, LENGTH_SET
     &   , U_SET, TRANS_SET, WEIGHT_POINT
     &   , I_TYPE_RESIDUAL, SCALING
     &   , N_TERM, W_ESFT, K_ESFT
     &   , I_SCALE_FUNCTION, NULL_VECTOR, P_ARRAY, T_ARRAY
     &   )
!
      WRITE(IU_STDOUT, '(//A)') 'ATTEMPTING TO FIND OPTIMAL SCALING '
      WRITE(IU_STDOUT, '(A/)')
     &   'FOR PRESSURE AND TEMPERATURE USING THE CONJUGATE GRADIENT '
     &   //'ALGORITHM,'
!
      WRITE(IU_MONITOR, '(//A)') 'ATTEMPTING TO FIND OPTIMAL SCALING '
      WRITE(IU_MONITOR, '(A/)')
     &   'FOR PRESSURE AND TEMPERATURE USING THE CONJUGATE GRADIENT '
     &   //'ALGORITHM,'
!
!     INITIALIZE THE ITERATING LOOP.
      ITERATION=0
      RESIDUAL_OLD=RESIDUAL_TRANS(N_SET, LENGTH_SET
     &   , U_SET, TRANS_SET, WEIGHT_POINT
     &   , I_TYPE_RESIDUAL, SCALING
     &   , N_TERM, W_ESFT, K_ESFT
     &   , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &   )
      CALL RESIDUAL_GRADIENT(N_SET, LENGTH_SET
     &   , U_SET, TRANS_SET, WEIGHT_POINT
     &   , I_TYPE_RESIDUAL, SCALING
     &   , N_TERM, W_ESFT, K_ESFT
     &   , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &   , GRADIENT)
      DO I=1, N_SCALE_VARIABLE(I_SCALE_FUNCTION)
         G(I)=-GRADIENT(I)
         H(I)=G(I)
      ENDDO
!
100   ITERATION=ITERATION+1
!
         H_L2_NORM=0.0D+00
         DO I=1, N_SCALE_VARIABLE(I_SCALE_FUNCTION)
            H_L2_NORM=H_L2_NORM+H(I)*H(I)
         ENDDO
         H_L2_NORM=SQRT(H_L2_NORM)
!
!        IF THIS IS VERY SMALL FURTHER MINIMIZATION IS POINTLESS.
         IF (H_L2_NORM.LT.SQRT(TOL_TEST)*RESIDUAL_TYPICAL) THEN
            CALL TERMINATE_SCALE(N_SET, LENGTH_SET
     &         , U_SET, TRANS_SET, WEIGHT_POINT
     &         , I_TYPE_RESIDUAL, SCALING
     &         , N_TERM, W_ESFT, K_ESFT
     &         , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &         , SCALE_VECTOR_LONG
     &         )
            IF (.NOT.L_RE_ENTRY.AND.(ITERATION.EQ.1)) THEN
!              IF THE PROGRAM EXITS WITHOUT PERFORMING A MINIMIZATION
!              AND THIS IS NOT A RUN TO REFINE AN EXISTING SET OF
!              SCALING PARAMETERS THEY ARE SET TO 0.
               DO I=1, N_SCALE_VARIABLE(I_SCALE_FUNCTION)
                  DO J=1, N_TERM
                     SCALE_VECTOR_LONG(I, J)=0.0D+00
                  ENDDO
               ENDDO
            ENDIF
            RETURN
         ENDIF
!
         IF (ITERATION.EQ.1) THEN
            SEARCH_DISTANCE=1.0D+00
         ENDIF
!
!        SEARCH ALONG H FOR A MINIMUM.
         SEARCH_DISTANCE=1.0D+00
         CALL LINE_SEARCH(IERR
     &      , N_SET, LENGTH_SET, U_SET, TRANS_SET
     &      , WEIGHT_POINT
     &      , I_TYPE_RESIDUAL, SCALING
     &      , N_TERM, W_ESFT, K_ESFT
     &      , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &      , H, H_L2_NORM, SEARCH_DISTANCE
     &      )
         IF (IERR.NE.I_NORMAL) THEN
            IF (IERR.EQ.I_ABORT_CALCULATION) THEN
               WRITE(IU_STDOUT, '(/A)')
     &            'ATTEMPT TO FIND SCALING ABANDONED.'
               WRITE(IU_MONITOR, '(/A)')
     &            'ATTEMPT TO FIND SCALING ABANDONED.'
            ENDIF
            RETURN
         ENDIF
!
!        ASSESS THE CONDITION FOR NORMAL TERMINIATION
         RESIDUAL_NEW=RESIDUAL_TRANS(N_SET, LENGTH_SET
     &      , U_SET, TRANS_SET, WEIGHT_POINT
     &      , I_TYPE_RESIDUAL, SCALING
     &      , N_TERM, W_ESFT, K_ESFT
     &      , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &      )
         WRITE(IU_MONITOR, '(/A, I4)') 'ITERATION = ', ITERATION
         WRITE(IU_MONITOR, '(A, 1PE16.9)') 'MEAN SQUARE ERROR = '
     &      , RESIDUAL_NEW
         WRITE(IU_MONITOR, '(A, (T22, 3(1PE16.9, 2X)))')
     &      'SCALING PARAMETERS = '
     &      , (SCALE_VECTOR(I), I=1, N_SCALE_VARIABLE(I_SCALE_FUNCTION))
         IF (ABS(RESIDUAL_NEW-RESIDUAL_OLD)
     &      .LT.TOL_TEST*RESIDUAL_TYPICAL) THEN
!           WE ARE DONE.
            CALL TERMINATE_SCALE(N_SET, LENGTH_SET
     &         , U_SET, TRANS_SET, WEIGHT_POINT
     &         , I_TYPE_RESIDUAL, SCALING
     &         , N_TERM, W_ESFT, K_ESFT
     &         , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &         , SCALE_VECTOR_LONG
     &         )
            RETURN
         ELSE
            RESIDUAL_OLD=RESIDUAL_NEW
         ENDIF
!
!        PROCEED TO DETERMINE THE CONJUGATE DIRECTION
!        FIND THE GRADIENT AT THE NEW MINIMUM.
         CALL RESIDUAL_GRADIENT(N_SET, LENGTH_SET
     &      , U_SET, TRANS_SET, WEIGHT_POINT
     &      , I_TYPE_RESIDUAL, SCALING
     &      , N_TERM, W_ESFT, K_ESFT
     &      , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &      , GRADIENT)
         GI2=0.0D+00
         DGI=0.0D+00
         DO I=1, N_SCALE_VARIABLE(I_SCALE_FUNCTION)
            GI2=GI2+G(I)*G(I)
            DGI=DGI+GRADIENT(I)*(GRADIENT(I)+G(I))
         ENDDO
         IF (GI2.LT.TOL_TEST*RESIDUAL_TYPICAL) THEN
            CALL TERMINATE_SCALE(N_SET, LENGTH_SET
     &         , U_SET, TRANS_SET, WEIGHT_POINT
     &         , I_TYPE_RESIDUAL, SCALING
     &         , N_TERM, W_ESFT, K_ESFT
     &         , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &         , SCALE_VECTOR_LONG
     &         )
            RETURN
         ENDIF
!
         DO I=1, N_SCALE_VARIABLE(I_SCALE_FUNCTION)
            G(I)=-GRADIENT(I)
!           IN THE CASE OF NON-QUADRATIC MINIMA THE DIRECTION OF
!           SEARCHING SHOULD BE RE-INITIALIZED AFTER A NUMBER OF STEPS
!           EQUAL TO THE NUMBER OF VARIABLES, SINCE OTHERWISE THE
!           SEARCHING DIRECTIONS MAY DRIFT AWAY FROM THE OPTIMAL ONES.
            IF (MOD(ITERATION, N_SCALE_VARIABLE(I_SCALE_FUNCTION)).EQ.0)
     &         THEN
               SEARCH_DISTANCE=1.0D+00
               H(I)=G(I)
            ELSE
               H(I)=G(I)+(DGI/GI2)*H(I)
            ENDIF
         ENDDO
!
      IF (ITERATION.LT.NP_MAX_ITERATION_CG) THEN
         GOTO 100
      ELSE
         WRITE(IU_STDOUT, '(/A)')
     &      'THE CONJUGATE GRADIENT SCALING HAS NOT CONVERGED:'
         WRITE(IU_STDOUT, '(/A)')
     &      'CURRENT VALUES OF THE SCALING PARAMETERS WILL BE GIVEN.'
         WRITE(IU_MONITOR, '(/A)') 'TOO MANY ITERATIONS WITHOUT '
     &      //'CONVERGENCE.'
         WRITE(IU_MONITOR, '(A)') 'CURRENT DATA WILL BE USED.'
         CALL TERMINATE_SCALE(N_SET, LENGTH_SET
     &      , U_SET, TRANS_SET, WEIGHT_POINT
     &      , I_TYPE_RESIDUAL, SCALING
     &      , N_TERM, W_ESFT, K_ESFT
     &      , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &      , SCALE_VECTOR_LONG
     &      )
         RETURN
      ENDIF
!
      END
