!+ Subroutine to calculate differences in source functions.
!
! Method:
!	Using the polynomial fit to the Planck function, values
!	of this function at the boundaries of layers are found
!	and differences across layers are determined. If the
!	Planckian is being taken to vary quadratically across 
!	the layer second differences are found.
!
! 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 DIFF_PLANCK_SOURCE(N_PROFILE, N_LAYER
     &   , N_DEG_FIT, THERMAL_COEFFICIENT
     &   , T_REF_PLANCK, T_LEVEL, T_GROUND
     &   , PLANCK_SOURCE, DIFF_PLANCK, PLANCK_GROUND
     &   , L_IR_SOURCE_QUAD, T, DIFF_PLANCK_2
     &   )
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_field.ucf'
      INCLUDE 'incl/dimensions_spec.ucf'
!
!     DUMMY ARGUMENTS.
      INTEGER	!, INTENT(IN)
     &     N_PROFILE
!             NUMBER OF PROFILES
     &   , N_LAYER
!             NUMBER OF LAYERS
     &   , N_DEG_FIT
!             DEGREE OF FITTING FUNCTION
      LOGICAL	!, INTENT(IN)
     &     L_IR_SOURCE_QUAD
!             IR-SOURCE QUADRATIC
      REAL*8	!, INTENT(IN)
     &     THERMAL_COEFFICIENT(0: NPD_THERMAL_COEFF-1)
!             COEFFICIENTS OF FIT TO PLANCK FNC
     &   , T_REF_PLANCK
!             PLANCKIAN REFERENCE TEMPERATURE
     &   , T_LEVEL(NPD_PROFILE, 0: NPD_LAYER)
!             TEMPERATURES ON LEVELS
     &   , T_GROUND(NPD_PROFILE)
!             TEMPERATURES AT GROUND
      REAL*8	!, INTENT(OUT)
     &     PLANCK_SOURCE(NPD_PROFILE, 0: NPD_LAYER)
!             PLANCK FUNCTION ON LEVELS
     &   , DIFF_PLANCK(NPD_PROFILE, NPD_LAYER)
!             DIFFERENCES IN PLANCKIAN FNC
     &   , DIFF_PLANCK_2(NPD_PROFILE, NPD_LAYER)
!             TWICE 2ND DIFFERENCES IN PLANCKIAN FUNCTION
     &   , T(NPD_PROFILE, 0: NPD_LAYER)
!             TEMPERATURES AT CENTRES OF LAYERS
     &   , PLANCK_GROUND(NPD_PROFILE)
!             PLANCKIAN FUNCTION AT GROUND
!
!
!     LOCAL VARIABLES.
      INTEGER
     &     I
!             LOOP VARIABLE
     &   , J
!             LOOP VARIABLE
     &   , L
!             LOOP VARIABLE
      REAL*8
     &     T_RATIO(NPD_PROFILE)
!             TEMPERATURE RATIO
!
!
!
!     CALCULATE THE CHANGE IN THE THERMAL SOURCE FUNCTION
!     ACROSS EACH LAYER FOR THE INFRA-RED PART OF THE SPECTRUM.
      DO L=1, N_PROFILE
         T_RATIO(L)=T_LEVEL(L, 0)/T_REF_PLANCK
         PLANCK_SOURCE(L, 0)
     &      =THERMAL_COEFFICIENT(N_DEG_FIT)
      ENDDO
      DO J=N_DEG_FIT-1, 0, -1
         DO L=1, N_PROFILE
            PLANCK_SOURCE(L, 0)
     &         =PLANCK_SOURCE(L, 0)
     &         *T_RATIO(L)+THERMAL_COEFFICIENT(J)
         ENDDO
      ENDDO
      DO I=1, N_LAYER
         DO L=1, N_PROFILE
            T_RATIO(L)=T_LEVEL(L, I)/T_REF_PLANCK
            PLANCK_SOURCE(L, I)
     &         =THERMAL_COEFFICIENT(N_DEG_FIT)
         ENDDO
         DO J=N_DEG_FIT-1, 0, -1
            DO L=1, N_PROFILE
               PLANCK_SOURCE(L, I)
     &            =PLANCK_SOURCE(L, I)
     &            *T_RATIO(L)+THERMAL_COEFFICIENT(J)
            ENDDO
         ENDDO
         DO L=1, N_PROFILE
            DIFF_PLANCK(L, I)=PLANCK_SOURCE(L, I)
     &         -PLANCK_SOURCE(L, I-1)
         ENDDO
      ENDDO
!     CALCULATE THE SECOND DIFFERENCE IF REQUIRED.
      IF (L_IR_SOURCE_QUAD) THEN
         DO I=1, N_LAYER
!           USE THE SECOND DIFFERENCE FOR TEMPORARY STORAGE.
!           OF THE PLANCKIAN AT THE MIDDLE OF THE LAYER.
            DO L=1, N_PROFILE
               T_RATIO(L)=T(L, I)/T_REF_PLANCK
               DIFF_PLANCK_2(L, I)
     &            =THERMAL_COEFFICIENT(N_DEG_FIT)
            ENDDO
            DO J=N_DEG_FIT-1, 0, -1
               DO L=1, N_PROFILE
                  DIFF_PLANCK_2(L, I)
     &               =DIFF_PLANCK_2(L, I)
     &               *T_RATIO(L)+THERMAL_COEFFICIENT(J)
               ENDDO
            ENDDO
            DO L=1, N_PROFILE
               DIFF_PLANCK_2(L, I)=2.0D+00*(PLANCK_SOURCE(L, I)
     &            +PLANCK_SOURCE(L, I-1)-2.0D+00*DIFF_PLANCK_2(L, I))
            ENDDO
         ENDDO
      ENDIF
!
!     SOURCE AT THE SURFACE.
      DO L=1, N_PROFILE
         T_RATIO(L)=T_GROUND(L)/T_REF_PLANCK
         PLANCK_GROUND(L)=THERMAL_COEFFICIENT(N_DEG_FIT)
      ENDDO
      DO J=N_DEG_FIT-1, 0, -1
         DO L=1, N_PROFILE
            PLANCK_GROUND(L)=PLANCK_GROUND(L)*T_RATIO(L)
     &         +THERMAL_COEFFICIENT(J)
         ENDDO
      ENDDO
!
!
!
      RETURN
      END
