C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved.            RAD_ALL.1986  
C                                                                                  RAD_ALL.1987  
C Use, duplication or disclosure of this code is subject to the                    RAD_ALL.1988  
C restrictions as set forth in the contract.                                       RAD_ALL.1989  
C                                                                                  RAD_ALL.1990  
C                Meteorological Office                                             RAD_ALL.1991  
C                London Road                                                       RAD_ALL.1992  
C                BRACKNELL                                                         RAD_ALL.1993  
C                Berkshire UK                                                      RAD_ALL.1994  
C                RG12 2SZ                                                          RAD_ALL.1995  
C                                                                                  RAD_ALL.1996  
C If no contract has been raised with this copy of the code, the use,              RAD_ALL.1997  
C duplication or disclosure of it is strictly prohibited.  Permission              RAD_ALL.1998  
C to do so must first be obtained in writing from the Head of Numerical            RAD_ALL.1999  
C Modelling at the above address.                                                  RAD_ALL.2000  
C ******************************COPYRIGHT******************************            RAD_ALL.2001  
C                                                                                  RAD_ALL.2002  
!+ Subroutine to calculate differences in source functions.                        RAD_ALL.2003  
!                                                                                  RAD_ALL.2004  
! Method:                                                                          RAD_ALL.2005  
!       Using the polynomial fit to the Planck function, values                    RAD_ALL.2006  
!       of this function at the boundaries of layers are found                     RAD_ALL.2007  
!       and differences across layers are determined. If the                       RAD_ALL.2008  
!       Planckian is being taken to vary quadratically across                      RAD_ALL.2009  
!       the layer second differences are found.                                    RAD_ALL.2010  
!                                                                                  RAD_ALL.2011  
! Current Owner of Code: J. M. Edwards                                             RAD_ALL.2012  
!                                                                                  RAD_ALL.2013  
! History:                                                                         RAD_ALL.2014  
!       Version         Date                    Comment                            RAD_ALL.2015  
!       4.0             27-07-95                Original Code                      RAD_ALL.2016  
!                                               (J. M. Edwards)                    RAD_ALL.2017  
!       4.1             14-03-95                Explicit formulation               RAD_ALL.2018  
!                                               for sea ice introduced.            RAD_ALL.2019  
!                                                                                  RAD_ALL.2020  
! Description of Code:                                                             RAD_ALL.2021  
!   FORTRAN 77  with extensions listed in documentation.                           RAD_ALL.2022  
!                                                                                  RAD_ALL.2023  
!- ---------------------------------------------------------------------           RAD_ALL.2024  
      SUBROUTINE DIFF_PLANCK_SOURCE(N_PROFILE, N_LAYER                             RAD_ALL.2025  
     &   , N_DEG_FIT, THERMAL_COEFFICIENT                                          RAD_ALL.2026  
     &   , T_REF_PLANCK, T_LEVEL, T_GROUND                                         RAD_ALL.2027  
     &   , PLANCK_SOURCE, DIFF_PLANCK, PLANCK_GROUND                               RAD_ALL.2028  
     &   , L_IR_SOURCE_QUAD, T, DIFF_PLANCK_2                                      RAD_ALL.2029  
     &   , N_FRAC_ICE_POINT, I_FRAC_ICE_POINT, ICE_FRACTION                        RAD_ALL.2030  
     &   , PLANCK_FREEZE_SEA                                                       RAD_ALL.2031  
     &   , NPD_PROFILE, NPD_LAYER, NPD_THERMAL_COEFF                               RAD_ALL.2032  
     &   )                                                                         RAD_ALL.2033  
!                                                                                  RAD_ALL.2034  
!                                                                                  RAD_ALL.2035  
      IMPLICIT NONE                                                                RAD_ALL.2036  
!                                                                                  RAD_ALL.2037  
!                                                                                  RAD_ALL.2038  
!     DUMMY ARRAY SIZES                                                            RAD_ALL.2039  
      INTEGER   !, INTENT(IN)                                                      RAD_ALL.2040  
     &     NPD_PROFILE                                                             RAD_ALL.2041  
!             MAXIMUM NUMBER OF PROFILES                                           RAD_ALL.2042  
     &   , NPD_LAYER                                                               RAD_ALL.2043  
!             MAXIMUM NUMBER OF LAYERS                                             RAD_ALL.2044  
     &   , NPD_THERMAL_COEFF                                                       RAD_ALL.2045  
!             NUMBER OF THERMAL COEFFICIENTS                                       RAD_ALL.2046  
!                                                                                  RAD_ALL.2047  
!     COMDECKS INCLUDED.                                                           RAD_ALL.2048  
C*L------------------COMDECK C_O_DG_C-----------------------------------           C_0_DG_C.2     
C ZERODEGC IS CONVERSION BETWEEN DEGREES CELSIUS AND KELVIN                        C_0_DG_C.3     
C TFS IS TEMPERATURE AT WHICH SEA WATER FREEZES                                    C_0_DG_C.4     
C TM IS TEMPERATURE AT WHICH FRESH WATER FREEZES AND ICE MELTS                     C_0_DG_C.5     
      REAL ZERODEGC,TFS,TM                                                         C_0_DG_C.6     
                                                                                   C_0_DG_C.7     
      PARAMETER(ZERODEGC=273.15,                                                   C_0_DG_C.8     
     &          TFS=271.35,                                                        C_0_DG_C.9     
     &          TM=273.15)                                                         C_0_DG_C.10    
C*----------------------------------------------------------------------           C_0_DG_C.11    
                                                                                   C_0_DG_C.12    
!                                                                                  RAD_ALL.2050  
!     DUMMY ARGUMENTS.                                                             RAD_ALL.2051  
      INTEGER   !, INTENT(IN)                                                      RAD_ALL.2052  
     &     N_PROFILE                                                               RAD_ALL.2053  
!             NUMBER OF PROFILES                                                   RAD_ALL.2054  
     &   , N_LAYER                                                                 RAD_ALL.2055  
!             NUMBER OF LAYERS                                                     RAD_ALL.2056  
     &   , N_DEG_FIT                                                               RAD_ALL.2057  
!             DEGREE OF FITTING FUNCTION                                           RAD_ALL.2058  
      LOGICAL   !, INTENT(IN)                                                      RAD_ALL.2059  
     &     L_IR_SOURCE_QUAD                                                        RAD_ALL.2060  
!             IR-SOURCE QUADRATIC                                                  RAD_ALL.2061  
      REAL      !, INTENT(IN)                                                      RAD_ALL.2062  
     &     THERMAL_COEFFICIENT(0: NPD_THERMAL_COEFF-1)                             RAD_ALL.2063  
!             COEFFICIENTS OF FIT TO PLANCK FNC                                    RAD_ALL.2064  
     &   , T_REF_PLANCK                                                            RAD_ALL.2065  
!             PLANCKIAN REFERENCE TEMPERATURE                                      RAD_ALL.2066  
     &   , T_LEVEL(NPD_PROFILE, 0: NPD_LAYER)                                      RAD_ALL.2067  
!             TEMPERATURES ON LEVELS                                               RAD_ALL.2068  
     &   , T_GROUND(NPD_PROFILE)                                                   RAD_ALL.2069  
!             TEMPERATURES AT GROUND                                               RAD_ALL.2070  
      REAL      !, INTENT(OUT)                                                     RAD_ALL.2071  
     &     PLANCK_SOURCE(NPD_PROFILE, 0: NPD_LAYER)                                RAD_ALL.2072  
!             PLANCK FUNCTION ON LEVELS                                            RAD_ALL.2073  
     &   , DIFF_PLANCK(NPD_PROFILE, NPD_LAYER)                                     RAD_ALL.2074  
!             DIFFERENCES IN PLANCKIAN FNC                                         RAD_ALL.2075  
     &   , DIFF_PLANCK_2(NPD_PROFILE, NPD_LAYER)                                   RAD_ALL.2076  
!             TWICE 2ND DIFFERENCES IN PLANCKIAN FUNCTION                          RAD_ALL.2077  
     &   , T(NPD_PROFILE, 0: NPD_LAYER)                                            RAD_ALL.2078  
!             TEMPERATURES AT CENTRES OF LAYERS                                    RAD_ALL.2079  
     &   , PLANCK_GROUND(NPD_PROFILE)                                              RAD_ALL.2080  
!             PLANCKIAN FUNCTION AT GROUND                                         RAD_ALL.2081  
!                                                                                  RAD_ALL.2082  
!     ARGUMENTS RELATING TO SEA ICE.                                               RAD_ALL.2083  
      INTEGER   !, INTENT(IN)                                                      RAD_ALL.2084  
     &     N_FRAC_ICE_POINT                                                        RAD_ALL.2085  
!             NUMBER OF POINTS WITH FRACTIONAL ICE COVER                           RAD_ALL.2086  
     &   , I_FRAC_ICE_POINT(NPD_PROFILE)                                           RAD_ALL.2087  
!             INDICES OF POINTS WITH FRACTIONAL ICE COVER                          RAD_ALL.2088  
      REAL      !, INTENT(IN)                                                      RAD_ALL.2089  
     &     ICE_FRACTION(NPD_PROFILE)                                               RAD_ALL.2090  
!             ICE FRACTION                                                         RAD_ALL.2091  
!                                                                                  RAD_ALL.2092  
      REAL      !, INTENT(OUT)                                                     RAD_ALL.2093  
     &     PLANCK_FREEZE_SEA                                                       RAD_ALL.2094  
!             PLANCK FUNCTION OVER FREEZING SEA                                    RAD_ALL.2095  
!                                                                                  RAD_ALL.2096  
!                                                                                  RAD_ALL.2097  
!     LOCAL VARIABLES.                                                             RAD_ALL.2098  
      INTEGER                                                                      RAD_ALL.2099  
     &     I                                                                       RAD_ALL.2100  
!             LOOP VARIABLE                                                        RAD_ALL.2101  
     &   , J                                                                       RAD_ALL.2102  
!             LOOP VARIABLE                                                        RAD_ALL.2103  
     &   , L                                                                       RAD_ALL.2104  
!             LOOP VARIABLE                                                        RAD_ALL.2105  
      REAL                                                                         RAD_ALL.2106  
     &     T_RATIO(NPD_PROFILE)                                                    RAD_ALL.2107  
!             TEMPERATURE RATIO                                                    RAD_ALL.2108  
!                                                                                  RAD_ALL.2109  
!     VARIABLES FOR FRACTIONAL ICE COVER.                                          RAD_ALL.2110  
      INTEGER                                                                      RAD_ALL.2111  
     &     LG                                                                      RAD_ALL.2112  
!             GATHERED LOOP VARIABLE                                               RAD_ALL.2113  
      REAL                                                                         RAD_ALL.2114  
     &     T_ICE_G(NPD_PROFILE)                                                    RAD_ALL.2115  
!             TEMPERATURE OF ICE SURFACE (GATHERED OVER SEA-ICE POINTS)            RAD_ALL.2116  
     &   , PLANCK_GROUND_G(NPD_PROFILE)                                            RAD_ALL.2117  
!             PLANCKIAN OF ICE SURFACE (GATHERED OVER SEA-ICE POINTS)              RAD_ALL.2118  
                                                                                   RAD_ALL.2119  
                                                                                   RAD_ALL.2120  
!                                                                                  RAD_ALL.2121  
!                                                                                  RAD_ALL.2122  
!                                                                                  RAD_ALL.2123  
!     CALCULATE THE CHANGE IN THE THERMAL SOURCE FUNCTION                          RAD_ALL.2124  
!     ACROSS EACH LAYER FOR THE INFRA-RED PART OF THE SPECTRUM.                    RAD_ALL.2125  
      DO L=1, N_PROFILE                                                            RAD_ALL.2126  
         T_RATIO(L)=T_LEVEL(L, 0)/T_REF_PLANCK                                     RAD_ALL.2127  
         PLANCK_SOURCE(L, 0)                                                       RAD_ALL.2128  
     &      =THERMAL_COEFFICIENT(N_DEG_FIT)                                        RAD_ALL.2129  
      ENDDO                                                                        RAD_ALL.2130  
      DO J=N_DEG_FIT-1, 0, -1                                                      RAD_ALL.2131  
         DO L=1, N_PROFILE                                                         RAD_ALL.2132  
            PLANCK_SOURCE(L, 0)                                                    RAD_ALL.2133  
     &         =PLANCK_SOURCE(L, 0)                                                RAD_ALL.2134  
     &         *T_RATIO(L)+THERMAL_COEFFICIENT(J)                                  RAD_ALL.2135  
         ENDDO                                                                     RAD_ALL.2136  
      ENDDO                                                                        RAD_ALL.2137  
      DO I=1, N_LAYER                                                              RAD_ALL.2138  
         DO L=1, N_PROFILE                                                         RAD_ALL.2139  
            T_RATIO(L)=T_LEVEL(L, I)/T_REF_PLANCK                                  RAD_ALL.2140  
            PLANCK_SOURCE(L, I)                                                    RAD_ALL.2141  
     &         =THERMAL_COEFFICIENT(N_DEG_FIT)                                     RAD_ALL.2142  
         ENDDO                                                                     RAD_ALL.2143  
         DO J=N_DEG_FIT-1, 0, -1                                                   RAD_ALL.2144  
            DO L=1, N_PROFILE                                                      RAD_ALL.2145  
               PLANCK_SOURCE(L, I)                                                 RAD_ALL.2146  
     &            =PLANCK_SOURCE(L, I)                                             RAD_ALL.2147  
     &            *T_RATIO(L)+THERMAL_COEFFICIENT(J)                               RAD_ALL.2148  
            ENDDO                                                                  RAD_ALL.2149  
         ENDDO                                                                     RAD_ALL.2150  
         DO L=1, N_PROFILE                                                         RAD_ALL.2151  
            DIFF_PLANCK(L, I)=PLANCK_SOURCE(L, I)                                  RAD_ALL.2152  
     &         -PLANCK_SOURCE(L, I-1)                                              RAD_ALL.2153  
         ENDDO                                                                     RAD_ALL.2154  
      ENDDO                                                                        RAD_ALL.2155  
!     CALCULATE THE SECOND DIFFERENCE IF REQUIRED.                                 RAD_ALL.2156  
      IF (L_IR_SOURCE_QUAD) THEN                                                   RAD_ALL.2157  
         DO I=1, N_LAYER                                                           RAD_ALL.2158  
!           USE THE SECOND DIFFERENCE FOR TEMPORARY STORAGE.                       RAD_ALL.2159  
!           OF THE PLANCKIAN AT THE MIDDLE OF THE LAYER.                           RAD_ALL.2160  
            DO L=1, N_PROFILE                                                      RAD_ALL.2161  
               T_RATIO(L)=T(L, I)/T_REF_PLANCK                                     RAD_ALL.2162  
               DIFF_PLANCK_2(L, I)                                                 RAD_ALL.2163  
     &            =THERMAL_COEFFICIENT(N_DEG_FIT)                                  RAD_ALL.2164  
            ENDDO                                                                  RAD_ALL.2165  
            DO J=N_DEG_FIT-1, 0, -1                                                RAD_ALL.2166  
               DO L=1, N_PROFILE                                                   RAD_ALL.2167  
                  DIFF_PLANCK_2(L, I)                                              RAD_ALL.2168  
     &               =DIFF_PLANCK_2(L, I)                                          RAD_ALL.2169  
     &               *T_RATIO(L)+THERMAL_COEFFICIENT(J)                            RAD_ALL.2170  
               ENDDO                                                               RAD_ALL.2171  
            ENDDO                                                                  RAD_ALL.2172  
            DO L=1, N_PROFILE                                                      RAD_ALL.2173  
               DIFF_PLANCK_2(L, I)=2.0E+00*(PLANCK_SOURCE(L, I)                    RAD_ALL.2174  
     &            +PLANCK_SOURCE(L, I-1)-2.0E+00*DIFF_PLANCK_2(L, I))              RAD_ALL.2175  
            ENDDO                                                                  RAD_ALL.2176  
         ENDDO                                                                     RAD_ALL.2177  
      ENDIF                                                                        RAD_ALL.2178  
!                                                                                  RAD_ALL.2179  
!     SOURCE AT THE SURFACE.                                                       RAD_ALL.2180  
      DO L=1, N_PROFILE                                                            RAD_ALL.2181  
         T_RATIO(L)=T_GROUND(L)/T_REF_PLANCK                                       RAD_ALL.2182  
         PLANCK_GROUND(L)=THERMAL_COEFFICIENT(N_DEG_FIT)                           RAD_ALL.2183  
      ENDDO                                                                        RAD_ALL.2184  
      DO J=N_DEG_FIT-1, 0, -1                                                      RAD_ALL.2185  
         DO L=1, N_PROFILE                                                         RAD_ALL.2186  
            PLANCK_GROUND(L)=PLANCK_GROUND(L)*T_RATIO(L)                           RAD_ALL.2187  
     &         +THERMAL_COEFFICIENT(J)                                             RAD_ALL.2188  
         ENDDO                                                                     RAD_ALL.2189  
      ENDDO                                                                        RAD_ALL.2190  
!                                                                                  RAD_ALL.2191  
!     WHERE THERE IS FRACTIONAL SEA-ICE THE FORMULATION MUST BE                    RAD_ALL.2192  
!     EXTENDED, BUT IT IS CONVENIENT TO CARRY OUT THE ABOVE OPERATION              RAD_ALL.2193  
!     AT ALL POINTS TO AVOID THE USE OF INDIRECT ADDRESSING.                       RAD_ALL.2194  
!                                                                                  RAD_ALL.2195  
!     CALCULATE THE SOURCE FUNCTION OVER OPEN SEA, ADOPTING THE MODEL'S            RAD_ALL.2196  
!     CONVENTION THAT THE TEMPAERTURE THERE IS FIXED.                              RAD_ALL.2197  
      T_RATIO(1)=TFS/T_REF_PLANCK                                                  RAD_ALL.2198  
      PLANCK_FREEZE_SEA=THERMAL_COEFFICIENT(N_DEG_FIT)                             RAD_ALL.2199  
      DO J=N_DEG_FIT-1, 0, -1                                                      RAD_ALL.2200  
         PLANCK_FREEZE_SEA=PLANCK_FREEZE_SEA*T_RATIO(1)                            RAD_ALL.2201  
     &      +THERMAL_COEFFICIENT(J)                                                RAD_ALL.2202  
      ENDDO                                                                        RAD_ALL.2203  
!                                                                                  RAD_ALL.2204  
!     DETERMINE THE TEMPERATURE OF THE ICE.                                        RAD_ALL.2205  
      DO L=1, N_FRAC_ICE_POINT                                                     RAD_ALL.2206  
         LG=I_FRAC_ICE_POINT(L)                                                    RAD_ALL.2207  
         T_ICE_G(L)=(T_GROUND(LG)                                                  RAD_ALL.2208  
     &      -TFS*(1.0E+00-ICE_FRACTION(LG)))/ICE_FRACTION(LG)                      RAD_ALL.2209  
      ENDDO                                                                        RAD_ALL.2210  
!                                                                                  RAD_ALL.2211  
!     CALCULATE THE SOURCE FUNCTION AT POINTS WITH FRACTIONAL ICE.                 RAD_ALL.2212  
      DO L=1, N_FRAC_ICE_POINT                                                     RAD_ALL.2213  
         T_RATIO(L)=T_ICE_G(L)/T_REF_PLANCK                                        RAD_ALL.2214  
         PLANCK_GROUND_G(L)=THERMAL_COEFFICIENT(N_DEG_FIT)                         RAD_ALL.2215  
      ENDDO                                                                        RAD_ALL.2216  
      DO J=N_DEG_FIT-1, 0, -1                                                      RAD_ALL.2217  
         DO L=1, N_FRAC_ICE_POINT                                                  RAD_ALL.2218  
            PLANCK_GROUND_G(L)=PLANCK_GROUND_G(L)*T_RATIO(L)                       RAD_ALL.2219  
     &         +THERMAL_COEFFICIENT(J)                                             RAD_ALL.2220  
         ENDDO                                                                     RAD_ALL.2221  
      ENDDO                                                                        RAD_ALL.2222  
!                                                                                  RAD_ALL.2223  
!     DETERMINE THE OVERALL PLANCKIAN FUNCTION OF THE SURFACE.                     RAD_ALL.2224  
      DO L=1, N_FRAC_ICE_POINT                                                     RAD_ALL.2225  
         LG=I_FRAC_ICE_POINT(L)                                                    RAD_ALL.2226  
         PLANCK_GROUND(LG)=ICE_FRACTION(LG)*PLANCK_GROUND_G(L)                     RAD_ALL.2227  
     &      +PLANCK_FREEZE_SEA*(1.0E+00-ICE_FRACTION(LG))                          RAD_ALL.2228  
      ENDDO                                                                        RAD_ALL.2229  
!                                                                                  RAD_ALL.2230  
!                                                                                  RAD_ALL.2231  
!                                                                                  RAD_ALL.2232  
      RETURN                                                                       RAD_ALL.2233  
      END     
