!+ 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. This version is used when a separate scaling
!	is applied to each ESFT term.
!
! 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 M_CONJUGATE_GRADIENT(IERR
     &   , N_SET, LENGTH_SET, U_SET, TRANS_SET, WEIGHT_POINT
     &   , N_TERM, W_ESFT, K_ESFT
     &   , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY)
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/prec_machine.ucf'
      INCLUDE 'incl/precision.pcf'
      INCLUDE 'incl/scale_parameters.acf'
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/dimensions_spec.ucf'
      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
      REAL*8	!, INTENT(IN)
     &     U_SET(NPD_AMOUNT, NPD_SET_SCALE)
!             AMOUNTS OF ABSORBER
     &   , TRANS_SET(NPD_AMOUNT, NPD_SET_SCALE)
!             TRANSMISSIONS
     &   , 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(NPD_SCALE_VARIABLE, NPD_ESFT_TERM)
!             PARAMETERS FOR SCALING
!
!     LOCAL ARGUMENTS.
      INTEGER
     &     ITERATION
!             ITERATION COUNT
     &   , I
!             LOOP VARIABLE
     &   , J
!             LOOP VARIABLE
     &   , K
!             LOOP VARIABLE
      REAL*8
     &     GRADIENT(NPD_SCALE_VARIABLE, NPD_ESFT_TERM)
!             GRADIENT OF RESIDUAL
     &   , RESIDUAL_OLD
!             OLD VALUE OF RESIDUAL
     &   , RESIDUAL_NEW
!             NEW VALUE OF RESIDUAL
     &   , RESIDUAL_TYPICAL
!             TYPICAL VALUE OF RESIDUAL
     &   , G(NPD_SCALE_VARIABLE, NPD_ESFT_TERM)
!             REVERSED GRADIENT
     &   , H(NPD_SCALE_VARIABLE, NPD_ESFT_TERM)
!             DIRECTION OF SEARCH
     &   , H_L2_NORM
!             L2 NORM OF H
     &   , DGI
!             TEMPORARY CONSTANT
     &   , GI2
!             TEMPORARY CONSTANT
     &   , SEARCH_DISTANCE
!             SEARCHING DISTANCE
     &   , NULL_VECTOR(NPD_SCALE_VARIABLE, NPD_ESFT_TERM)
!             NULL VECTOR
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     M_RESIDUAL_GRADIENT, M_TERMINATE_SCALE, M_LINE_SEARCH
!
!     FUNCTIONS CALLED:
      REAL*8
     &     M_RESIDUAL_TRANS
!             RESIDUAL FUNCTION
      EXTERNAL
     &     M_RESIDUAL_TRANS
!
!
      INCLUDE 'incl/scale_fnc.dcf'
!
!
!
!     SET THE NULL VECTOR.
      DO I=1, NPD_SCALE_VARIABLE
         DO J=1, NPD_ESFT_TERM
            NULL_VECTOR(I, J)=0.0D+00
         ENDDO
      ENDDO
!
!     CALCULATE A TYPICAL VALUE OF THE RESIDUAL FOR ESTIMATION
!     OF ERRORS.
      RESIDUAL_TYPICAL=M_RESIDUAL_TRANS(N_SET, LENGTH_SET
     &   , U_SET, TRANS_SET, WEIGHT_POINT
     &   , 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 BY CONJUGATE GRADIENT '
      WRITE(IU_STDOUT, '(A/)')
     &   'MINIMIZATION OF THE ERROR IN THE RESIDUALS IN '
     &   //'TRANSMISSIVITY.'
!
      WRITE(IU_MONITOR, '(//A)') 'ATTEMPTING TO FIND OPTIMAL SCALING '
      WRITE(IU_MONITOR, '(A)')
     &   'FOR PRESSURE AND TEMPERATURE BY CONJUGATE GRADIENT '
      WRITE(IU_MONITOR, '(A/)')
     &   'MINIMIZATION OF THE ERROR IN THE RESIDUALS IN '
     &   //'TRANSMISSIVITY.'
!
!     INITIALIZE THE ITERATING LOOP.
      ITERATION=0
      RESIDUAL_OLD=M_RESIDUAL_TRANS(N_SET, LENGTH_SET
     &   , U_SET, TRANS_SET, WEIGHT_POINT
     &   , N_TERM, W_ESFT, K_ESFT
     &   , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &   )
      CALL M_RESIDUAL_GRADIENT(N_SET, LENGTH_SET
     &   , U_SET, TRANS_SET, WEIGHT_POINT
     &   , N_TERM, W_ESFT, K_ESFT
     &   , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &   , GRADIENT)
      DO K=1, N_TERM
         DO I=1, N_SCALE_VARIABLE(I_SCALE_FUNCTION)
            G(I, K)=-GRADIENT(I, K)
            H(I, K)=G(I, K)
         ENDDO
      ENDDO
!
100   ITERATION=ITERATION+1
!
         H_L2_NORM=0.0D+00
         DO K=1, N_TERM
            DO I=1, N_SCALE_VARIABLE(I_SCALE_FUNCTION)
               H_L2_NORM=H_L2_NORM+H(I, K)*H(I, K)
            ENDDO
         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 M_TERMINATE_SCALE(N_SET, LENGTH_SET
     &         , U_SET, TRANS_SET, WEIGHT_POINT
     &         , N_TERM, W_ESFT, K_ESFT
     &         , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &         )
            IF (ITERATION.EQ.1) THEN
!              NO MINIMIZATION HAS BEEN PERFORMED SO WE SET THE SCALING
!              PARAMETERS TO 0.
               DO I=1, N_SCALE_VARIABLE(I_SCALE_FUNCTION)
                  DO J=1, N_TERM
                     SCALE_VECTOR(I, J)=0.0D+00
                  ENDDO
               ENDDO
            ENDIF
            RETURN
         ENDIF
!
         IF (ITERATION.EQ.1) THEN
!           CALCULATE A TYPICAL SEARCHING DISTANCE.
            SEARCH_DISTANCE=MIN(DBLE(2.0D+00)
     &         , MAX(DBLE(0.5D+00), RESIDUAL_OLD/H_L2_NORM))
         ENDIF
!
!        SEARCH ALONG H FOR A MINIMUM.
         CALL M_LINE_SEARCH(IERR
     &      , N_SET, LENGTH_SET, U_SET, TRANS_SET
     &      , WEIGHT_POINT
     &      , 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=M_RESIDUAL_TRANS(N_SET, LENGTH_SET
     &      , U_SET, TRANS_SET, WEIGHT_POINT
     &      , 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/, (4(1X, 1PE16.9)))')
     &      'SCALING PARAMETERS = '
     &      , ((SCALE_VECTOR(I, K), I=1
     &      , N_SCALE_VARIABLE(I_SCALE_FUNCTION)), K=1, N_TERM)
         IF (ABS(RESIDUAL_NEW-RESIDUAL_OLD)
     &      .LT.TOL_TEST*RESIDUAL_TYPICAL) THEN
!           WE ARE DONE.
            CALL M_TERMINATE_SCALE(N_SET, LENGTH_SET
     &         , U_SET, TRANS_SET, WEIGHT_POINT
     &         , N_TERM, W_ESFT, K_ESFT
     &         , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &         )
            RETURN
         ELSE
            RESIDUAL_OLD=RESIDUAL_NEW
         ENDIF
!
!        PROCEED TO DETERMINE THE CONJUGATE DIRECTION
!        FIND THE GRADIENT AT THE NEW MINIMUM.
         CALL M_RESIDUAL_GRADIENT(N_SET, LENGTH_SET
     &      , U_SET, TRANS_SET, WEIGHT_POINT
     &      , 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 K=1, N_TERM
            DO I=1, N_SCALE_VARIABLE(I_SCALE_FUNCTION)
               GI2=GI2+G(I, K)*G(I, K)
               DGI=DGI+GRADIENT(I, K)*(GRADIENT(I, K)+G(I, K))
            ENDDO
         ENDDO
         IF (GI2.LT.TOL_TEST*RESIDUAL_TYPICAL) THEN
            CALL M_TERMINATE_SCALE(N_SET, LENGTH_SET
     &         , U_SET, TRANS_SET, WEIGHT_POINT
     &         , N_TERM, W_ESFT, K_ESFT
     &         , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &         )
            RETURN
         ENDIF
!
         DO K=1, N_TERM
            DO I=1, N_SCALE_VARIABLE(I_SCALE_FUNCTION)
               G(I, K)=-GRADIENT(I, K)
               IF (MOD(ITERATION, N_TERM
     &            *N_SCALE_VARIABLE(I_SCALE_FUNCTION))
     &            .EQ.0) THEN
                  H(I, K)=G(I, K)
               ELSE
                  H(I, K)=G(I, K)+(DGI/GI2)*H(I, K)
               ENDIF
            ENDDO
         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 M_TERMINATE_SCALE(N_SET, LENGTH_SET
     &      , U_SET, TRANS_SET, WEIGHT_POINT
     &      , N_TERM, W_ESFT, K_ESFT
     &      , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &      )
         RETURN
      ENDIF
!
      END
