!+ Subroutine to fit to cloud properties using conjugate gradients.
!
! Method:
!	The Polak-Ribiere from of the conjugate gradient algorithm
!	is implemented to find a fit to cloudy properties.
!
! 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_CLOUD(IERR
     &   , N_DATA, I_FIT_TYPE, N_PARAMETER
     &   , VOL_FRACTION, RADIUS_EFFECTIVE
     &   , DATA, PARM
     &   , RESIDUAL, RESIDUAL_GRADIENT
     &   )
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_spec.ucf'
      INCLUDE 'incl/cloud_fit_parm.acf'
      INCLUDE 'incl/prec_machine.ucf'
      INCLUDE 'incl/precision.pcf'
      INCLUDE 'incl/tol_conj_grad.tcf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/error.pcf'
!
!     DUMMY ARGUMENTS.
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      INTEGER	!, INTENT(IN)
     &     N_DATA
!             NUMBER OF DATA POINTS
     &   , I_FIT_TYPE
!             TYPE OF FIT REQUIRED
     &   , N_PARAMETER
!             NUMBER OF FITTING PARAMETERS
      REAL*8	!, INTENT(IN)
     &     DATA(N_DATA)
!             DATA TO BE FITTED
     &   , VOL_FRACTION(N_DATA)
!             VOLUME FRACTION OF ABSORBER
     &   , RADIUS_EFFECTIVE(N_DATA)
!             EFFECTIVE RADIUS OF PARTICLES
      REAL*8	!, INTENT(OUT)
     &     PARM(N_PARAMETER)
!             FITTING PARAMETERS
!
      REAL*8	!, INTENT(IN)
     &     RESIDUAL
!             NAME OF RESIDUAL ROUTINE
!
      EXTERNAL	
     &     RESIDUAL
!             NAME OF RESIDUAL ROUTINE
     &   , RESIDUAL_GRADIENT
!             GRADIENT OF RESIDUAL
!
!     LOCAL ARGUMENTS
      INTEGER
     &     ITERATION
!             NUMBER OF CG ITERATION
     &   , I
!             LOOP VARIABLE
      REAL*8
     &     RESIDUAL_OLD
!             FORMER VALUE OF RESIDUAL
     &   , RESIDUAL_NEW
!             NEWER VALUE OF RESIDUAL
     &   , GRADIENT(NPD_CLOUD_PARAMETER)
!             GRADIENT
     &   , G(NPD_CLOUD_PARAMETER)
!             REVERSED GRADIENT
     &   , H(NPD_CLOUD_PARAMETER)
!             SEARCH DIRECTION
     &   , H_L2_NORM
!             L2-NORM OF H
     &   , SEARCH_DISTANCE
!             ESTIMATED SEARCH DISTANCE
     &   , GI2
!             TEMPORARY CONSTANT
     &   , DGI
!             TEMPORARY CONSTANT
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     LINE_SEARCH_CLOUD
!
!
!
!     CALCULATE THE INITIAL VALUE OF THE RESIDUAL.
      RESIDUAL_OLD=RESIDUAL(N_DATA, I_FIT_TYPE, N_PARAMETER
     &   , VOL_FRACTION, RADIUS_EFFECTIVE
     &   , DATA, PARM
     &   )
!
!     INITIALIZE THE GRADIENT AND THE DIRECTION OF SEARCHING.
      CALL RESIDUAL_GRADIENT(N_DATA, I_FIT_TYPE, N_PARAMETER
     &   , VOL_FRACTION, RADIUS_EFFECTIVE
     &   , DATA, PARM
     &   , GRADIENT
     &   )
       DO I=1, N_PARAMETER
          G(I)=-GRADIENT(I)
          H(I)=G(I)
       ENDDO
!
!     INITIALIZE THE COUNT OF ITERATIONS.
      ITERATION=0
100   ITERATION=ITERATION+1
!
         H_L2_NORM=0.0D+00
         DO I=1, N_PARAMETER
            H_L2_NORM=H_L2_NORM+H(I)*H(I)
         ENDDO
         H_L2_NORM=SQRT(H_L2_NORM)
!
!        RETURN IF H IS VIRTUALLY 0.
         IF (H_L2_NORM.LT.TOL_GRADIENT) RETURN
!
         IF (ITERATION.EQ.1) THEN
!           ESTIMATE A SUITABLE SEARCHING DISTANCE.
            SEARCH_DISTANCE=1.0D+00
         ENDIF
!
!        SEARCH ALONG THE DIRECTION H.
         SEARCH_DISTANCE=1.0D+00
         CALL LINE_SEARCH_CLOUD(IERR
     &      , N_DATA, I_FIT_TYPE, N_PARAMETER
     &      , VOL_FRACTION, RADIUS_EFFECTIVE
     &      , DATA, PARM
     &      , H, H_L2_NORM, SEARCH_DISTANCE
     &      , RESIDUAL
     &      )
         IF (IERR.NE.I_NORMAL) THEN
            IF (IERR.EQ.I_ABORT_CALCULATION) THEN
!              LINE SEARCH WAS ABORTED: WE RECOVER FROM THIS.
               IERR=I_NORMAL
            ELSE
               RETURN
            ENDIF
         ENDIF
!
!        ASSESS THE CONDITION FOR NORMAL TERMINATION.
         RESIDUAL_NEW=RESIDUAL(N_DATA, I_FIT_TYPE, N_PARAMETER
     &      , VOL_FRACTION, RADIUS_EFFECTIVE
     &      , DATA, PARM
     &      )
         WRITE(IU_MONITOR, '(/A, I4)') 'ITERATION = ', ITERATION
         WRITE(IU_MONITOR, '(A, 1PE16.9)') 'MEAN SQUARE ERROR = '
     &      , RESIDUAL_NEW
         WRITE(IU_MONITOR, '(A, (2(1X, 1PE16.9)))')
     &      'SCALING PARAMETERS = '
     &      , (PARM(I), I=1, N_PARAMETER)
         IF (ABS(RESIDUAL_NEW-RESIDUAL_OLD).LT.TOL_RESIDUAL) THEN
            RETURN
         ELSE
            RESIDUAL_OLD=RESIDUAL_NEW
         ENDIF
!
!        NOW SEARCH IN A DIRECTION CONJUGATE TO THE CURRENT ONE.
         CALL RESIDUAL_GRADIENT(N_DATA, I_FIT_TYPE, N_PARAMETER
     &      , VOL_FRACTION, RADIUS_EFFECTIVE
     &      , DATA, PARM
     &      , GRADIENT
     &      )
         GI2=0.0D+00
         DGI=0.0D+00
         DO  I=1, N_PARAMETER
            GI2=GI2+G(I)*G(I)
            DGI=DGI+GRADIENT(I)*(GRADIENT(I)+G(I))
         ENDDO
         IF (GI2.LT.TOL_GRADIENT) RETURN
!
         DO  I=1, N_PARAMETER
            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_PARAMETER).EQ.0)
     &         THEN
               H(I)=G(I)
               SEARCH_DISTANCE=1.0D+00
            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)') 'TOO MANY ITERATIONS IN THE '
     &      //'CONJUGATE GRADIENT ROUTINE '
         WRITE(IU_STDOUT, '(A)') 'WITHOUT CONVERGENCE.'
         WRITE(IU_STDOUT, '(A)') 'CURRENT DATA WILL BE USED.'
         RETURN
      ENDIF
!
!
      END
