!+ Subroutine to fit parametrization of scattering by ice in a band.
!
! Method:
!	For each quantity in the parametrization a numerical scaling
!	of the data is performed to improve conditioning and a 
!	conjugate gradient algorithm is called to find the coefficients 
!	in the parametrization.
!
! 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 ICE_CLOUD_FIT_BAND(IERR
     &   , N_DATA, I_FIT_TYPE
     &   , DENSITY_PARTICLE
     &   , VOL_FRACTION, RADIUS_EFFECTIVE
     &   , SCATTERING, ABSORPTION, ASYMMETRY
     &   , ICE_PARAMETER_LIST)
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADERS.
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/dimensions_spec.ucf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/opt_variable.pcf'
      INCLUDE 'incl/ice_cloud_parametrization.pcf'
      INCLUDE 'incl/ice_cloud_parametrization.ncf'
      INCLUDE 'incl/error.pcf'
!
!     DUMMY ARGUMENTS.
      INTEGER 	!, INTENT(OUT)
     &     IERR
!             SUBROUTINE ERROR FLAG
      INTEGER 	!, INTENT(IN)
     &     N_DATA
!             NUMBER OF DATA POINTS
     &   , I_FIT_TYPE
!             TYPE OF FITTING REQUIRED
      REAL*8	!, INTENT(IN)
     &     DENSITY_PARTICLE
!             DENSITY OF PARTICLES
     &   , VOL_FRACTION(N_DATA)
!             VOLUME FRACTION OF ABSORBER
     &   , RADIUS_EFFECTIVE(N_DATA)
!             EFFECTIVE RADIUS OF DROPS
     &   , SCATTERING(N_DATA)
!             SCATTERING IN BAND
     &   , ABSORPTION(N_DATA)
!             ABSORPTION IN BAND
     &   , ASYMMETRY(N_DATA)
!             ASYMMETRY AVERAGED ACROSS BAND
      REAL*8	!, INTENT(OUT)
     &     ICE_PARAMETER_LIST(NPD_CLOUD_PARAMETER)
!             PARAMETERS OF CLOUD SCHEME
!
!     LOCAL VARIABLES.
      INTEGER
     &     I
!             LOOP VARIABLE
     &   , N_PARAMETER
!             NUMBER OF FITTING PARAMETERS
      REAL*8
     &     EXTINCTION(NPD_MIE_BLOCK)
!             TOTAL EXTINCTION
     &   , SS_ALBEDO(NPD_MIE_BLOCK)
!             SINGLE SCATTERING ALBEDO
     &   , RADIUS_SCALED(NPD_MIE_BLOCK)
!             SCALED EFFECTIVE RADIUS
     &   , EXTINCTION_SCALED(NPD_MIE_BLOCK)
!             DIMENSIONLESS EXTINCTION
     &   , PARAMETER_LOCAL(NPD_CLOUD_PARAMETER
     &        , NPD_CLOUD_OPTICAL_VARIABLE)
!             PARAMETERS EACH OPTICAL PROP.
     &   , RADIUS_MEAN
!             MEAN EFFECTIVE RADIUS
     &   , MEAN_EXTINCTION
!             MEAN EXTINCTION/UNIT VOL.
     &   , SS_COALBEDO
!             COALBEDO
!
!     FUNCTIONS CALLED:
      REAL*8
     &     ICE_EXTINCTION_RESIDUAL
!             RESIDUAL ROUTINE FOR EXTINCTION
     &   , ICE_SS_ALBEDO_RESIDUAL
!             RESIDUAL ROUTINE FOR SCATTERING
     &   , ICE_ASYMMETRY_RESIDUAL
!             RESIDUAL ROUTINE FOR ASYMMETRY
     &   , ICE_SCATTERING_RESIDUAL
!             RESIDUAL ROUTINE FOR SCATTERING
     &   , ICE_ABSORPTION_RESIDUAL
!             RESIDUAL ROUTINE FOR SCATTERING
      EXTERNAL
     &     ICE_EXTINCTION_RESIDUAL, ICE_SS_ALBEDO_RESIDUAL
     &   , ICE_ASYMMETRY_RESIDUAL, ICE_SCATTERING_RESIDUAL
     &   , ICE_ABSORPTION_RESIDUAL
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     CONJUGATE_GRADIENT_CLOUD, BUILD_ICE_PARAMETER_LIST
     &   , ICE_EXTINCTION_GRADIENT, ICE_SS_ALBEDO_GRADIENT
     &   , ICE_ASYMMETRY_GRADIENT, ICE_SCATTERING_GRADIENT
     &   , ICE_ABSORPTION_GRADIENT
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/initial.pcf'
      INCLUDE 'incl/ice_cloud_parm_data.dcf'
!
!
!
!     SEPARATE THE PARAMETRIZATION SCHEMES ACCORDING TO THE VARIABLES TO
!     BE FITTED.
      IF (I_FIT_TYPE.EQ.IP_SLINGO_SCHRECKER_ICE) THEN
!        FIT THE EXTINCTION, SINGLE SCATTERING ALBEDO AND ASYMMETRY.
         DO I=1, N_DATA
            EXTINCTION(I)=ABSORPTION(I)+SCATTERING(I)
            SS_ALBEDO(I)=SCATTERING(I)/EXTINCTION(I)
         ENDDO
         DO I=1, N_ICE_CLOUD_FIT_PARAMETER(I_FIT_TYPE, IP_EXTINCTION)
            PARAMETER_LOCAL(I, IP_EXTINCTION)=0.0D+00
         ENDDO
         DO I=1, N_ICE_CLOUD_FIT_PARAMETER(I_FIT_TYPE, IP_SS_ALBEDO)
            PARAMETER_LOCAL(I, IP_SS_ALBEDO)=0.0D+00
         ENDDO
         DO I=1, N_ICE_CLOUD_FIT_PARAMETER(I_FIT_TYPE, IP_ASYMMETRY)
            PARAMETER_LOCAL(I, IP_ASYMMETRY)=0.0D+00
         ENDDO
!
!        SCALE THE QUANTITIES FOR FITTING.
         RADIUS_MEAN=0.0D+00
         SS_COALBEDO=0.0D+00
         MEAN_EXTINCTION=0.0D+00
         DO I=1, N_DATA
            RADIUS_MEAN=RADIUS_MEAN+RADIUS_EFFECTIVE(I)
            SS_COALBEDO=MAX(SS_COALBEDO, 1.0D+00-SS_ALBEDO(I))
            MEAN_EXTINCTION=MEAN_EXTINCTION
     &         +EXTINCTION(I)/VOL_FRACTION(I)
         ENDDO
         RADIUS_MEAN=RADIUS_MEAN/DBLE(N_DATA)
         MEAN_EXTINCTION=MEAN_EXTINCTION/DBLE(N_DATA)
         DO I=1, N_DATA
            EXTINCTION_SCALED(I)=EXTINCTION(I)/(MEAN_EXTINCTION
     &         *VOL_FRACTION(I))
            RADIUS_SCALED(I)=RADIUS_EFFECTIVE(I)/RADIUS_MEAN
            SS_ALBEDO(I)=(1.0D+00-SS_ALBEDO(I))/SS_COALBEDO
         ENDDO
!
         N_PARAMETER
     &      =N_ICE_CLOUD_FIT_PARAMETER(I_FIT_TYPE, IP_EXTINCTION)
!
         WRITE(IU_MONITOR, '(/A ,/A)')
     &      '------------------------------------'
     &      , 'FITTING THE EXTINCTION.'
!
         CALL CONJUGATE_GRADIENT_CLOUD(IERR
     &      , N_DATA, I_FIT_TYPE, N_PARAMETER
     &      , VOL_FRACTION, RADIUS_SCALED
     &      , EXTINCTION_SCALED, PARAMETER_LOCAL(1, IP_EXTINCTION)
     &      , ICE_EXTINCTION_RESIDUAL, ICE_EXTINCTION_GRADIENT
     &      )
         IF (IERR.NE.I_NORMAL) THEN
            IF (IERR.EQ.I_ABORT_CALCULATION) THEN
!              RECOVER FORM MILD ERRORS OR RETURN FROM FATAL ONES.
               IERR=I_NORMAL
            ELSE
               RETURN
            ENDIF
         ENDIF
!
!        RESCALE PARAMETERS AND CONVERT DIVIDE BY THE DENSITY OF THE
!        PARTICLES TO ALLOW FOR THE SUBSEQUENT USE OF 
!        MASS MIXING RATIOS.
         PARAMETER_LOCAL(1, IP_EXTINCTION)
     &      =MEAN_EXTINCTION*PARAMETER_LOCAL(1, IP_EXTINCTION)
     &      /DENSITY_PARTICLE
         PARAMETER_LOCAL(2, IP_EXTINCTION)=RADIUS_MEAN
     &      *MEAN_EXTINCTION*PARAMETER_LOCAL(2, IP_EXTINCTION)
     &      /DENSITY_PARTICLE
!
!
!        THE QUANTITY PASSED TO THE SUBROUTINE IS ACTUALLY THE SCALED
!        CO-ALBEDO.
         N_PARAMETER=N_ICE_CLOUD_FIT_PARAMETER(I_FIT_TYPE, IP_SS_ALBEDO)
!
         WRITE(IU_MONITOR, '(/A ,/A)')
     &      '------------------------------------'
     &      , 'FITTING THE COALBEDO.'
!
         CALL CONJUGATE_GRADIENT_CLOUD(IERR
     &      , N_DATA, I_FIT_TYPE, N_PARAMETER
     &      , VOL_FRACTION, RADIUS_SCALED
     &      , SS_ALBEDO, PARAMETER_LOCAL(1, IP_SS_ALBEDO)
     &      , ICE_SS_ALBEDO_RESIDUAL, ICE_SS_ALBEDO_GRADIENT
     &      )
         IF (IERR.NE.I_NORMAL) THEN
            IF (IERR.EQ.I_ABORT_CALCULATION) THEN
!              RECOVER FORM MILD ERRORS OR RETURN FROM FATAL ONES.
               IERR=I_NORMAL
            ELSE
               RETURN
            ENDIF
         ENDIF
!
!        RESCALE PARAMETERS
         PARAMETER_LOCAL(1, IP_SS_ALBEDO)=SS_COALBEDO
     &      *PARAMETER_LOCAL(1, IP_SS_ALBEDO)
         PARAMETER_LOCAL(2, IP_SS_ALBEDO)=SS_COALBEDO
     &      *PARAMETER_LOCAL(2, IP_SS_ALBEDO)/RADIUS_MEAN
!
!
         N_PARAMETER=N_ICE_CLOUD_FIT_PARAMETER(I_FIT_TYPE, IP_ASYMMETRY)
!
         WRITE(IU_MONITOR, '(/A ,/A)')
     &      '------------------------------------'
     &      , 'FITTING THE ASYMMETRY.'
!
         CALL CONJUGATE_GRADIENT_CLOUD(IERR
     &      , N_DATA, I_FIT_TYPE, N_PARAMETER
     &      , VOL_FRACTION, RADIUS_SCALED
     &      , ASYMMETRY, PARAMETER_LOCAL(1, IP_ASYMMETRY)
     &      , ICE_ASYMMETRY_RESIDUAL, ICE_ASYMMETRY_GRADIENT
     &      )
         IF (IERR.NE.I_NORMAL) THEN
!           RECOVER FORM MILD ERRORS OR RETURN FROM FATAL ONES.
            IF (IERR.EQ.I_ABORT_CALCULATION) THEN
               IERR=I_NORMAL
            ELSE
               RETURN
            ENDIF
         ENDIF
!        RESCALE PARAMETERS
         PARAMETER_LOCAL(2, IP_ASYMMETRY)
     &      =PARAMETER_LOCAL(2, IP_ASYMMETRY)/RADIUS_MEAN
!
!        ASSIGN THE LOCAL ARRAYS OF PARAMETERS TO THE MAIN ARRAY.
         CALL BUILD_ICE_PARAMETER_LIST(I_FIT_TYPE
     &      , PARAMETER_LOCAL, ICE_PARAMETER_LIST)
!
!
!     CIRRUS PARAMETERIZATION BASED iON ADT.
      ELSE IF (I_FIT_TYPE.EQ.IP_CIRRUS_ADT) THEN
!        FIT THE SINGLE SCATTERING ALBEDO AND ASYMMETRY.
         DO I=1, N_DATA
            SS_ALBEDO(I)=SCATTERING(I)/(ABSORPTION(I)+SCATTERING(I))
         ENDDO
         DO I=1, N_ICE_CLOUD_FIT_PARAMETER(I_FIT_TYPE, IP_SS_ALBEDO)
            PARAMETER_LOCAL(I, IP_SS_ALBEDO)=0.0D+00
         ENDDO
         DO I=1, N_ICE_CLOUD_FIT_PARAMETER(I_FIT_TYPE, IP_ASYMMETRY)
            PARAMETER_LOCAL(I, IP_ASYMMETRY)=0.0D+00
         ENDDO
!
!        SCALE THE QUANTITIES FOR FITTING.
         RADIUS_MEAN=0.0D+00
         DO I=1, N_DATA
            RADIUS_MEAN=RADIUS_MEAN+RADIUS_EFFECTIVE(I)
         ENDDO
         RADIUS_MEAN=RADIUS_MEAN/DBLE(N_DATA)
         DO I=1, N_DATA
            RADIUS_SCALED(I)=RADIUS_EFFECTIVE(I)/RADIUS_MEAN
         ENDDO
!
!        THE EXTINCTION IS NOT FITTED IN THIS CASE.
!
         N_PARAMETER=N_ICE_CLOUD_FIT_PARAMETER(I_FIT_TYPE, IP_SS_ALBEDO)
!
         WRITE(IU_MONITOR, '(/A ,/A)')
     &      '------------------------------------'
     &      , 'FITTING THE ALBEDO.'
!
         CALL CONJUGATE_GRADIENT_CLOUD(IERR
     &      , N_DATA, I_FIT_TYPE, N_PARAMETER
     &      , VOL_FRACTION, RADIUS_SCALED
     &      , SS_ALBEDO, PARAMETER_LOCAL(1, IP_SS_ALBEDO)
     &      , ICE_SS_ALBEDO_RESIDUAL, ICE_SS_ALBEDO_GRADIENT
     &      )
         IF (IERR.NE.I_NORMAL) THEN
            IF (IERR.EQ.I_ABORT_CALCULATION) THEN
!              RECOVER FORM MILD ERRORS OR RETURN FROM FATAL ONES.
               IERR=I_NORMAL
            ELSE
               RETURN
            ENDIF
         ENDIF
!
!        RESCALE PARAMETERS
         PARAMETER_LOCAL(2, IP_SS_ALBEDO)
     &      =PARAMETER_LOCAL(2, IP_SS_ALBEDO)/RADIUS_MEAN
!
!
         N_PARAMETER=N_ICE_CLOUD_FIT_PARAMETER(I_FIT_TYPE, IP_ASYMMETRY)
!
         WRITE(IU_MONITOR, '(/A ,/A)')
     &      '------------------------------------'
     &      , 'FITTING THE ASYMMETRY.'
!
         CALL CONJUGATE_GRADIENT_CLOUD(IERR
     &      , N_DATA, I_FIT_TYPE, N_PARAMETER
     &      , VOL_FRACTION, RADIUS_SCALED
     &      , ASYMMETRY, PARAMETER_LOCAL(1, IP_ASYMMETRY)
     &      , ICE_ASYMMETRY_RESIDUAL, ICE_ASYMMETRY_GRADIENT
     &      )
         IF (IERR.NE.I_NORMAL) THEN
!           RECOVER FORM MILD ERRORS OR RETURN FROM FATAL ONES.
            IF (IERR.EQ.I_ABORT_CALCULATION) THEN
               IERR=I_NORMAL
            ELSE
               RETURN
            ENDIF
         ENDIF
!
!        ASSIGN THE LOCAL ARRAYS OF PARAMETERS TO THE MAIN ARRAY.
         CALL BUILD_ICE_PARAMETER_LIST(I_FIT_TYPE
     &      , PARAMETER_LOCAL, ICE_PARAMETER_LIST)
!
      ELSE
         WRITE(IU_ERR, '(/A)')
     &      '***ERROR: ILLEGAL CLOUD PARAMETRIZATION SELECTED.'
         IERR=I_ERR_FATAL
         RETURN
      ENDIF
!
!
      RETURN
      END
