!+ Subroutine to initialize for CG-algorithm to find scaling.
!
! Method:
!	An initial guess is made for the scaling parameters and
!	scaled independent variables are defined for numerical 
!	convenience. The CG-algorithm is then called. If different
!	terms are allowed their own scalings, a second 
!	CG-minimzation is envoked.
!
! 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 SCALE_FIT(IERR, I_SCALE
     &   , N_OTHER
     &   , LENGTH_SET, U_SET, TRANS_SET, P_SET, T_SET
     &   , N_TERM, K_ESFT, W_ESFT
     &   , I_SCALE_FUNCTION, SCALE_VECTOR)
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/dimensions_spec.ucf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/prec_machine.ucf'
      INCLUDE 'incl/precision.pcf'
      INCLUDE 'incl/esft_scale.pcf'
      INCLUDE 'incl/type_residual.pcf'
      INCLUDE 'incl/tol_scale.tcf'
      INCLUDE 'incl/scale_fnc.pcf'
      INCLUDE 'incl/error.pcf'
!
!     DUMMY ARGUMENTS
      INTEGER	!, INTENT(IN)
     &     IERR
!             ERROR FLAG
     &   , I_SCALE
!             TYPE OF SCALING
     &   , I_SCALE_FUNCTION
!             TYPE OF SCALING FUNCTION
     &   , I_TYPE_RESIDUAL
!             TYPE OF RESIDUAL USED
     &   , N_OTHER
!             NUMBER OF SIMILAR SETS
     &   , LENGTH_SET(0: NPD_SET_SCALE)
!             SIZES OF SETS
     &   , LENGTH_SET_FIT(0: NPD_SET_SCALE)
!             REDUCED SIZES OF SETS
     &   , N_TERM
!             NUMBER OF ESFT TERMS
      REAL*8	!, INTENT(IN)
     &     U_SET(NPD_AMOUNT, 0: NPD_SET_SCALE)
!             AMOUNTS OF ABSORBER
     &   , TRANS_SET(NPD_AMOUNT, 0: NPD_SET_SCALE)
!             TRANSMISSIONS
     &   , P_SET(0: NPD_SET_SCALE)
!             PRESSURES
     &   , T_SET(0: NPD_SET_SCALE)
!             TEMPERATURES
     &   , K_ESFT(NPD_ESFT_TERM)
!             ESFT EXPONENTS
     &   , W_ESFT(NPD_ESFT_TERM)
!             ESFT WEIGHTS
      REAL*8	!, INTENT(OUT)
     &     SCALE_VECTOR(NPD_SCALE_VARIABLE, NPD_ESFT_TERM)
!             SCALING PARAMETERS
!
!     LOCAL VARIABLES.
      INTEGER
     &     N_POINT
!             TOTAL NUMBER OF DATA POINTS
     &   , I
!             LOOP VARIABLE
     &   , J
!             LOOP VARIABLE
      REAL*8
     &     P_ARRAY(NPD_SET_SCALE)
!             ARRAY OF PRESSURES
     &   , T_ARRAY(NPD_SET_SCALE)
!             ARRAY OF TEMPERATURES
     &   , WEIGHT_POINT(NPD_SET_SCALE, NPD_AMOUNT)
!             WEIGHTS FOR DATA POINTS
     &   , SCALING(NPD_AMOUNT, NPD_SET_SCALE)
!             SCALINGS REQUIRED
     &   , RANGE_T
!             RANGE OF TEMPERATURES ENCOUNTERED
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     FIND_SCALING_POINT, FIND_SCALING_MEAN
     &   , CONJUGATE_GRADIENT, M_CONJUGATE_GRADIENT
!
!
      INCLUDE 'incl/scale_fnc.dcf'
      INCLUDE 'incl/set_residual.ucf'
!
!
!
!     NORMALIZE THE PRESSURE AND THE TEMPERATURE TO A CONVENIENT FORM 
!     FOR FITTING AND SELECT THE SETS OF DATA NEEDED, 
!     RECORDING THEIR LENGTHS.
      IF ( (I_SCALE_FUNCTION.EQ.IP_SCALE_POWER_QUAD).OR.
     &   (I_SCALE_FUNCTION.EQ.IP_SCALE_DOPPLER_QUAD) ) THEN
         RANGE_T=0.0D+00
         DO I=1, N_OTHER
            RANGE_T=MAX(RANGE_T, ABS(T_SET(I)-T_SET(0)))
         ENDDO
      ENDIF
      DO I=1, N_OTHER
         IF (I_SCALE_FUNCTION.EQ.IP_SCALE_POWER_LAW) THEN
            P_ARRAY(I)=LOG(P_SET(I)/P_SET(0))
            T_ARRAY(I)=LOG(T_SET(I)/T_SET(0))
         ELSE IF (I_SCALE_FUNCTION.EQ.IP_SCALE_POWER_QUAD) THEN
            P_ARRAY(I)=LOG(P_SET(I)/P_SET(0))
            T_ARRAY(I)=(T_SET(I)-T_SET(0))/RANGE_T
         ELSE IF (I_SCALE_FUNCTION.EQ.IP_SCALE_DOPPLER_QUAD) THEN
            P_ARRAY(I)=P_SET(I)/P_SET(0)
            T_ARRAY(I)=(T_SET(I)-T_SET(0))/RANGE_T
         ENDIF
      ENDDO
!
!     INITIALIZE THE WEIGHTING ARRAYS.
!     HERE ALL POINTS ARE WEIGHTED EQUALLY.
      N_POINT=0
      DO I=1, N_OTHER
         N_POINT=N_POINT+LENGTH_SET(I)
      ENDDO
      DO I=1, N_OTHER
         DO J=1, LENGTH_SET(I)
            WEIGHT_POINT(I, J)=1.0D+00/DBLE(N_POINT)
         ENDDO
      ENDDO
!
!     INITIALIZE THE SCALING PARAMETERS
      IF (I_SCALE_FUNCTION.EQ.IP_SCALE_POWER_LAW) THEN
         DO I=1, N_TERM
            SCALE_VECTOR(1, I)=1.0D+00
            SCALE_VECTOR(2, I)=0.5D+00
         ENDDO
      ELSE IF (I_SCALE_FUNCTION.EQ.IP_SCALE_POWER_QUAD) THEN
         DO I=1, N_TERM
            SCALE_VECTOR(1, I)=1.0D+00
            SCALE_VECTOR(2, I)=0.0D+00
            SCALE_VECTOR(3, I)=0.0D+00
         ENDDO
      ELSE IF (I_SCALE_FUNCTION.EQ.IP_SCALE_DOPPLER_QUAD) THEN
         DO I=1, N_TERM
            SCALE_VECTOR(1, I)=1.0D+00
            SCALE_VECTOR(2, I)=-2.0D+00
            SCALE_VECTOR(3, I)=0.0D+00
            SCALE_VECTOR(4, I)=0.0D+00
         ENDDO
      ELSE
         WRITE(IU_ERR, '(/A)') '*** ERROR: AN UNSUPPORTED SCALING '
     &      //'FUNCTION HAS BEEN REQUESTED.'
         IERR=I_ERR_FATAL
         RETURN
      ENDIF
!
      WRITE(IU_STDOUT, '(//A)')
     &   '------- CACULATION OF SCALING OF AMOUNT OF ABSORBER -------'
      WRITE(IU_MONITOR, '(//A)')
     &   '------- CACULATION OF SCALING OF AMOUNT OF ABSORBER -------'
!
      IF (I_TYPE_RESIDUAL.EQ.IP_SCALE_TRANS_RESIDUAL) THEN
         DO I=1, N_OTHER
            LENGTH_SET_FIT(I)=LENGTH_SET(I)
         ENDDO
      ELSE IF (I_TYPE_RESIDUAL.EQ.IP_SCALE_FULL_SCL_RES) THEN
         DO I=1, N_OTHER
            LENGTH_SET_FIT(I)=LENGTH_SET(I)
            DO J=LENGTH_SET(I), 1, -1
!              THIS FIT CANNOT BE DONE WITH SMALL TRANSMITTANCES
!              SO THE LENGTH IS REDUCED IF SUCH ARE ENCOUNTERED.
               IF (TRANS_SET(J, I).LT.TOL_TRANSMISSION) THEN
                  LENGTH_SET_FIT(I)=J-1
               ELSE
                  CALL FIND_SCALING_POINT(IERR, N_TERM, W_ESFT, K_ESFT
     &               , U_SET(J, I), TRANS_SET(J, I)
     &               , SCALING(J, I)
     &               )
               ENDIF
            ENDDO
         ENDDO
      ELSE IF (I_TYPE_RESIDUAL.EQ.IP_SCALE_MEAN_SCL_RES) THEN
         DO I=1, N_OTHER
            LENGTH_SET_FIT(I)=LENGTH_SET(I)
            CALL FIND_SCALING_MEAN(IERR, N_TERM, W_ESFT, K_ESFT
     &         , LENGTH_SET(I), U_SET(1, I), TRANS_SET(1, I)
     &         , SCALING(1, I)
     &         )
         ENDDO
      ENDIF
!
!     PASS TO THE CONJUGATE GRADIENT ALGORITHM.
      IF (I_SCALE.EQ.IP_SCALE_TERM) THEN
         WRITE(IU_STDOUT, '(/A/)')
     &      'CALCULATE ONE SCALING FOR ALL TERMS AS A STARTING POINT.'
         WRITE(IU_MONITOR, '(/A/)')
     &      'CALCULATE ONE SCALING FOR ALL TERMS AS A STARTING POINT.'
      ENDIF
      CALL CONJUGATE_GRADIENT(IERR, .FALSE.
     &   , N_OTHER, LENGTH_SET_FIT(1), U_SET(1, 1), TRANS_SET(1, 1)
     &   , I_TYPE_RESIDUAL, SCALING, WEIGHT_POINT
     &   , N_TERM, W_ESFT, K_ESFT
     &   , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &   )
      IF ( (IERR.NE.I_NORMAL).AND.
     &     (IERR.NE.I_ABORT_CALCULATION) ) RETURN

!
      IF (I_TYPE_RESIDUAL.NE.IP_SCALE_TRANS_RESIDUAL) THEN
!        REPEAT THE CALCULATION MINIMIZING THE ERROR IN THE TRANSMISSION
!        EVEN IF THE PREVIOUS OPTIMIZATION SHOULD HAVE FAILED TO 
!        CONVERGE, WE MAY TRY ONCE MORE.
         IF (IERR.EQ.I_ABORT_CALCULATION) IERR=I_NORMAL
         WRITE(IU_STDOUT, '(///A)')
     &      'REPEATING THE CALCULATION MINIMIZING THE ERROR IN THE '
     &      //'TRANSMISSION.'
         WRITE(IU_MONITOR, '(/A)')
     &      'REPEATING THE CALCULATION MINIMIZING THE ERROR IN THE '
     &      //'TRANSMISSION.'
!        THE FULL LENGTH OF THE SETS IS USED HERE REGARDLESSLY.
         CALL CONJUGATE_GRADIENT(IERR, .TRUE.
     &      , N_OTHER, LENGTH_SET(1), U_SET(1, 1), TRANS_SET(1, 1)
     &      , IP_SCALE_TRANS_RESIDUAL, SCALING, WEIGHT_POINT
     &      , N_TERM, W_ESFT, K_ESFT
     &      , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &      )
      ENDIF
!
      IF (IERR.EQ.I_ABORT_CALCULATION) THEN
!        IF THE CALCULATION HAS FAILED WE SET ALL SCALING 
!        PARAMETERS TO 0.
         IF (I_SCALE_FUNCTION.EQ.IP_SCALE_POWER_LAW) THEN
            DO I=1, N_TERM
               SCALE_VECTOR(1, I)=0.0D+00
               SCALE_VECTOR(2, I)=0.0D+00
            ENDDO
         ELSE IF (I_SCALE_FUNCTION.EQ.IP_SCALE_POWER_QUAD) THEN
            DO I=1, N_TERM
               SCALE_VECTOR(1, I)=0.0D+00
               SCALE_VECTOR(2, I)=0.0D+00
               SCALE_VECTOR(3, I)=0.0D+00
            ENDDO
         ELSE IF (I_SCALE_FUNCTION.EQ.IP_SCALE_DOPPLER_QUAD) THEN
            DO I=1, N_TERM
               SCALE_VECTOR(1, I)=0.0D+00
!              THIS IS LATER TO BE EXPONENTIATED SO 0 IS INAPPROPRIATE.
               SCALE_VECTOR(2, I)=-1.0D+02
               SCALE_VECTOR(3, I)=0.0D+00
               SCALE_VECTOR(4, I)=0.0D+00
            ENDDO
         ENDIF
         RETURN
      ENDIF
!
      IF (I_SCALE.EQ.IP_SCALE_TERM) THEN
         WRITE(IU_STDOUT, '(/A/)')
     &      'NOW ALLOW A DIFFERENT SCALING FOR EACH TERM.'
         WRITE(IU_MONITOR, '(/A/)')
     &      'NOW ALLOW A DIFFERENT SCALING FOR EACH TERM.'
         CALL M_CONJUGATE_GRADIENT(IERR
     &      , N_OTHER, LENGTH_SET(1), U_SET(1, 1), TRANS_SET(1, 1)
     &      , WEIGHT_POINT
     &      , N_TERM, W_ESFT, K_ESFT
     &      , I_SCALE_FUNCTION, SCALE_VECTOR, P_ARRAY, T_ARRAY
     &      )
!        THERE IS NO TEST FOR CONVERGENCE HERE SINCE THE VALUES FOR JUST
!        ONE TERM CAN BE USED.
      ENDIF
!
!     NORMALIZE THE SCALING VARIABLES TO THE CORRECT FORM
!     FOR THE RADIATION CODE.
      IF (I_SCALE_FUNCTION.EQ.IP_SCALE_POWER_QUAD) THEN
         DO I=1, N_TERM
            SCALE_VECTOR(2, I)=T_SET(0)*SCALE_VECTOR(2, I)/RANGE_T
            SCALE_VECTOR(3, I)=SCALE_VECTOR(3, I)*(T_SET(0)/RANGE_T)**2
         ENDDO
      ENDIF
      IF (I_SCALE_FUNCTION.EQ.IP_SCALE_DOPPLER_QUAD) THEN
         DO I=1, N_TERM
            SCALE_VECTOR(2, I)=EXP(SCALE_VECTOR(2, I))*P_SET(0)
            SCALE_VECTOR(3, I)=T_SET(0)*SCALE_VECTOR(3, I)/RANGE_T
            SCALE_VECTOR(4, I)=SCALE_VECTOR(4, I)*(T_SET(0)/RANGE_T)**2
         ENDDO
      ENDIF
!
!
      RETURN
      END
