      SUBROUTINE PLANCK3(BF, BST, URBF, BF0, URBF0, DBF, TFULL, GTT, 
     1                   IB, ITILE,
     2                   IL1, IL2, ILG, LAY, LEV, NTILE)
C
C     * MAY 01,2012 - M.LAZARE. NEW VERSION FOR GCM16:
C     *                         - ASSUME ISOTHERMAL LAPSE RATE FOR
C     *                           MOON LAYER INSTEAD OF EXTRAPOLATION.
C     * DEC 05,2007 - M.LAZARE. PREVIOUS VERSION PLANCK2 FOR GCM15G/H/I:
C     *                         - XX NOW INTERNAL WORK ARRAY.
C     * JUN 19/2003 - J.LI.     PREVIOUS VERSION PLANCK UP TO GCM15F.
C----------------------------------------------------------------------C
C     CALCULATION OF PLANCK FUNCTION IN VALID RANGE 120 - 360 K        C
C                                                                      C
C     BF:    BLACKBODY INTENSITY INTEGRATED OVER EACH BAND AT EACH     C
C            LEVEL IN UNITS W / M^2 / SR.                              C
C     BS:    THE BLACKBODY INTENSITY AT EACH TILED SURFACE.             C
C     BF0:   THE BLACKBODY INTENSITY AT THE TOA (ASSUME 210 K).        C
C     TFULL: TEMPERATURE AT EACH LEVEL                                 C
C     GT:    TEMPERATURE AT GROUND                                     C
C     U:     1 / DIFFUSIVITY FACTOR                                    C
C     URBF:  U TIMES THE DIFFERENCE OF LOG(BF) FOR TWO NEIGHBOR LEVELS C
C            USED FOR EXPONENTIAL SOURCE FUNCTION (LI, 2002 JAS P3302) C
C     URBF0: U TIMES THE DIFFERENCE OF LOG(BF) FOR TOA AND LEVEL 1     C
C     DBF:   DIFFERENCE OF BF FOR TWO NEIGHBOR LEVELS USED FOR LINEAR  C
C            SOURCE FUNCTION (LI, 2002 JAS P3302)                      C
C                                                                      C
C     0.0040816327 = 1 / 245 (245 THE STANDARD TEMPERATURE FOR POLY.   C
C     FIT)                                                             C
C----------------------------------------------------------------------C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C
      REAL    BF(ILG,LEV), BF0(ILG), URBF(ILG,LAY), URBF0(ILG),
     1        DBF(ILG,LAY), TFULL(ILG,LEV)
      REAL    BST(ILG,NTILE), GTT(ILG,NTILE)
      INTEGER ITILE(ILG,NTILE)
C
C     * INTERNAL WORK ARRAY.
C
      REAL XX(ILG,LAY)
      REAL    XP(6,9)
      DATA  U / 0.60653066 /, RTSTAND / 0.0040816327 /
      DATA ((XP(I,J), I = 1, 6), J = 1, 9) /
     1  -2.9876423E+00,    1.3660089E+01,   -1.2944461E+01,
     1   1.1775748E+01,   -1.9236798E+01,    2.3584435E+01,
     2  -1.6414103E+00,    1.1898535E+01,   -1.1262182E+01,
     2   1.0236863E+01,   -1.6677772E+01,    2.0423136E+01,
     3   6.5215205E-01,    9.2657366E+00,   -8.5872301E+00,
     3   7.6765044E+00,   -1.2287254E+01,    1.4990547E+01,
     4   1.5442143E+00,    7.2253228E+00,   -6.7811515E+00,
     4   6.1572299E+00,   -9.8725011E+00,    1.1997278E+01,
     5   1.2777580E+00,    6.1257638E+00,   -5.7906013E+00,
     5   5.3296782E+00,   -8.7529282E+00,    1.0741367E+01,
     6   2.1005257E+00,    5.2376301E+00,   -4.8915631E+00,
     6   4.5030997E+00,   -7.3199981E+00,    8.9204038E+00,
     7   2.9091223E+00,    3.9860795E+00,   -3.5829565E+00,
     7   3.2692193E+00,   -5.1799711E+00,    6.2157752E+00,
     8   2.7856424E+00,    2.8179582E+00,   -2.3780464E+00,
     8   2.1432949E+00,   -3.4540206E+00,    4.1814100E+00,
     9   2.4623332E+00,    1.8731841E+00,   -1.3659538E+00,
     9   1.1484948E+00,   -1.5975564E+00,    1.7791135E+00 /
C
      DO 100 I = IL1, IL2
        DT          =  TFULL(I,1) * RTSTAND - 1.0
        XXT         =  XP(1,IB) + DT * (XP(2,IB) + DT * (XP(3,IB) +
     2                            DT * (XP(4,IB) + DT * (XP(5,IB) +
     3                            DT *  XP(6,IB) ))))
C
        DT          =  TFULL(I,1) * RTSTAND - 1.0
        XX0         =  XP(1,IB) + DT * (XP(2,IB) + DT * (XP(3,IB) +
     1                            DT * (XP(4,IB) + DT * (XP(5,IB) +
     2                            DT *  XP(6,IB) ))))
        DT          =  TFULL(I,2) * RTSTAND - 1.
        XX(I,1)     =  XP(1,IB) + DT * (XP(2,IB) + DT * (XP(3,IB) +
     1                            DT * (XP(4,IB) + DT * (XP(5,IB) +
     2                            DT *  XP(6,IB) ))))
C
        BF0(I)      =  EXP(XXT)
        URBF0(I)    =  U * (XX0 - XXT)
        BF(I,1)     =  EXP(XX0)
        BF(I,2)     =  EXP(XX(I,1))
        DBF(I,1)    =  BF(I,2) - BF(I,1)
        URBF(I,1)   =  U * (XX(I,1) - XX0)
  100 CONTINUE
      DO M = 1, NTILE
         DO I = IL1, IL2
            IF (ITILE(I,M) .GT. 0) THEN
               DT       =  GTT(I,M) * RTSTAND - 1.0
               BST(I,M) =  EXP( XP(1,IB) +
     1                     DT * (XP(2,IB) + DT * (XP(3,IB) +
     2                     DT * (XP(4,IB) + DT * (XP(5,IB) +
     3                     DT *  XP(6,IB) )))) )  
            END IF
         END DO ! I
      END DO ! M
C
      DO 205 K = 2, LAY
        KM1 = K - 1
        KP1 = K + 1
        DO 200 I = IL1, IL2
          DT        =  TFULL(I,KP1) * RTSTAND - 1.0
          XX(I,K)   =  XP(1,IB) + DT * (XP(2,IB) + DT * (XP(3,IB) +
     1                            DT * (XP(4,IB) + DT * (XP(5,IB) +
     2                            DT *  XP(6,IB) ))))
C
          BF(I,KP1) =  EXP(XX(I,K))
          DBF(I,K)  =  BF(I,KP1) - BF(I,K)
          URBF(I,K) =  U * (XX(I,K) - XX(I,KM1))
  200     CONTINUE
  205   CONTINUE
C
      RETURN
      END
