!OPTIONS XOPT(HSFUN)
SUBROUTINE SWDE &
  &( KIDIA, KFDIA, KLON &
  &, PGG  , PREF , PRMUZ, PTO1, PW &
  &, PRE1 , PRE2 , PTR1 , PTR2 &
  &)

!**** *SWDE* - DELTA-EDDINGTON IN A CLOUDY LAYER

!     PURPOSE.
!     --------
!           COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY OF A CLOUDY
!     LAYER USING THE DELTA-EDDINGTON'S APPROXIMATION.

!**   INTERFACE.
!     ----------
!          *SWDE* IS CALLED BY *SWR*, *SWNI*


!        EXPLICIT ARGUMENTS :
!        --------------------
! PGG    : (KLON)             ; ASSYMETRY FACTOR
! PREF   : (KLON)             ; REFLECTIVITY OF THE UNDERLYING LAYER
! PRMUZ  : (KLON)             ; COSINE OF SOLAR ZENITH ANGLE
! PTO1   : (KLON)             ; OPTICAL THICKNESS
! PW     : (KLON)             ; SINGLE SCATTERING ALBEDO
!     ==== OUTPUTS ===
! PRE1   : (KLON)             ; LAYER REFLECTIVITY ASSUMING NO
!                             ; REFLECTION FROM UNDERLYING LAYER
! PTR1   : (KLON)             ; LAYER TRANSMISSIVITY ASSUMING NO
!                             ; REFLECTION FROM UNDERLYING LAYER
! PRE2   : (KLON)             ; LAYER REFLECTIVITY ASSUMING
!                             ; REFLECTION FROM UNDERLYING LAYER
! PTR2   : (KLON)             ; LAYER TRANSMISSIVITY ASSUMING
!                             ; REFLECTION FROM UNDERLYING LAYER

!        IMPLICIT ARGUMENTS :   NONE
!        --------------------

!     METHOD.
!     -------

!          STANDARD DELTA-EDDINGTON LAYER CALCULATIONS.

!     EXTERNALS.
!     ----------

!          NONE

!     REFERENCE.
!     ----------

!        SEE RADIATION'S PART OF THE MODEL'S DOCUMENTATION AND
!        ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE IFS

!     AUTHOR.
!     -------
!        JEAN-JACQUES MORCRETTE  *ECMWF*

!     MODIFICATIONS.
!     --------------
!        ORIGINAL : 88-12-15
!                   96-05-30 Michel Deque (security in EXP()) 
!                   08-03-28 Hubert Gallee(lower/upper limit on PTR1,2)
    
!     ------------------------------------------------------------------


!     ------------------------------------------------------------------

!*       0.1   ARGUMENTS
!              ---------


#include "tsmbkind.h"

IMPLICIT NONE


!     DUMMY INTEGER SCALARS
INTEGER_M :: KFDIA
INTEGER_M :: KIDIA
INTEGER_M :: KLON

REAL_B :: PGG(KLON),PREF(KLON),PRMUZ(KLON),PTO1(KLON),PW(KLON)
REAL_B :: PRE1(KLON),PRE2(KLON),PTR1(KLON),PTR2(KLON)

!     LOCAL INTEGER SCALARS
INTEGER_M :: JL

!     LOCAL REAL SCALARS
REAL_B :: ZA11, ZA12, ZA13, ZA21, ZA22, ZA23, ZALPHA,&
          &ZAM2B, ZAP2B, ZARG, ZARG2, ZB21, ZB22, ZB23, &
          &ZBETA, ZC1A, ZC1B, ZC2A, ZC2B, ZDENA, ZDENB, &
          &ZDT, ZEXKM, ZEXKP, ZEXMU0, ZFF, ZGP, ZRI0A, &
          &ZRI0B, ZRI0C, ZRI0D, ZRI1A, ZRI1B, ZRI1C, &
          &ZRI1D, ZRK, ZRM2, ZRP, ZTOP, ZWCP, ZWM, ZX1, &
          &ZX2, ZXM2P, ZXP2P



!     ------------------------------------------------------------------

!*         1.      DELTA-EDDINGTON CALCULATIONS


DO JL   =   KIDIA,KFDIA

!*         1.1     SET UP THE DELTA-MODIFIED PARAMETERS


  ZFF = PGG(JL)*PGG(JL)
  ZGP = PGG(JL)/(_ONE_+PGG(JL))
  ZTOP = (_ONE_- PW(JL) * ZFF) * PTO1(JL)
  ZWCP = (1-ZFF)* PW(JL) /(_ONE_- PW(JL) * ZFF)
  ZDT = _TWO_/3._JPRB
  ZX1 = _ONE_-ZWCP*ZGP
  ZWM = _ONE_-ZWCP
  ZRM2 =  PRMUZ(JL) * PRMUZ(JL)
  ZRK = SQRT(3._JPRB*ZWM*ZX1)
  ZX2 = 4._JPRB*(_ONE_-ZRK*ZRK*ZRM2)
  ZRP=ZRK/ZX1
  ZALPHA = 3._JPRB*ZWCP*ZRM2*(_ONE_+ZGP*ZWM)/ZX2
  ZBETA = 3._JPRB*ZWCP* PRMUZ(JL) *(_ONE_+3._JPRB*ZGP*ZRM2*ZWM)/ZX2
!     ZARG=MIN(ZTOP/PRMUZ(JL),200.)
  ZARG=MAX(-085._JPRB,MIN(ZTOP/PRMUZ(JL),085._JPRB))
  ZEXMU0=EXP(-ZARG)
  ZARG2=MIN(ZRK*ZTOP,085._JPRB)
  ZEXKP=EXP(ZARG2)
  ZEXKM = _ONE_/ZEXKP
  ZXP2P = _ONE_+ZDT*ZRP
  ZXM2P = _ONE_-ZDT*ZRP
  ZAP2B = ZALPHA+ZDT*ZBETA
  ZAM2B = ZALPHA-ZDT*ZBETA

!*         1.2     WITHOUT REFLECTION FROM THE UNDERLYING LAYER


  ZA11 = ZXP2P
  ZA12 = ZXM2P
  ZA13 = ZAP2B
  ZA22 = ZXP2P*ZEXKP
  ZA21 = ZXM2P*ZEXKM
  ZA23 = ZAM2B*ZEXMU0
  ZDENA = ZA11 * ZA22 - ZA21 * ZA12
  ZC1A = (ZA22*ZA13-ZA12*ZA23)/ZDENA
  ZC2A = (ZA11*ZA23-ZA21*ZA13)/ZDENA
  ZRI0A = ZC1A+ZC2A-ZALPHA
  ZRI1A = ZRP*(ZC1A-ZC2A)-ZBETA
  PRE1(JL) = (ZRI0A-ZDT*ZRI1A)/ PRMUZ(JL)
  ZRI0B = ZC1A*ZEXKM+ZC2A*ZEXKP-ZALPHA*ZEXMU0
  ZRI1B = ZRP*(ZC1A*ZEXKM-ZC2A*ZEXKP)-ZBETA*ZEXMU0
  PTR1(JL) = ZEXMU0+(ZRI0B+ZDT*ZRI1B)/ PRMUZ(JL)

  PRE1(JL) = MAX(_ZERO_,MIN(PRE1(JL),_ONE_))   ! lower/upper limit (Hubert Gallee, LGGE, 28-03-2008)
  PTR1(JL) = MAX(_ZERO_,MIN(PTR1(JL),_ONE_))   ! lower/upper limit (Hubert Gallee, LGGE, 28-03-2008)

!*         1.3     WITH REFLECTION FROM THE UNDERLYING LAYER


  ZB21 = ZA21- PREF(JL) *ZXP2P*ZEXKM
  ZB22 = ZA22- PREF(JL) *ZXM2P*ZEXKP
  ZB23 = ZA23- PREF(JL) *ZEXMU0*(ZAP2B - PRMUZ(JL) )
  ZDENB = ZA11 * ZB22 - ZB21 * ZA12
  ZC1B = (ZB22*ZA13-ZA12*ZB23)/ZDENB
  ZC2B = (ZA11*ZB23-ZB21*ZA13)/ZDENB
  ZRI0C = ZC1B+ZC2B-ZALPHA
  ZRI1C = ZRP*(ZC1B-ZC2B)-ZBETA
  PRE2(JL) = (ZRI0C-ZDT*ZRI1C) / PRMUZ(JL)
  ZRI0D = ZC1B*ZEXKM + ZC2B*ZEXKP - ZALPHA*ZEXMU0
  ZRI1D = ZRP * (ZC1B*ZEXKM - ZC2B*ZEXKP) - ZBETA*ZEXMU0
  PTR2(JL) = ZEXMU0 + (ZRI0D + ZDT*ZRI1D) / PRMUZ(JL)

  PRE2(JL) = MAX(_ZERO_,MIN(PRE2(JL),_ONE_))   ! lower/upper limit (Hubert Gallee, LGGE, 28-03-2008)
  PTR2(JL) = MAX(_ZERO_,MIN(PTR2(JL),_ONE_))   ! lower/upper limit (Hubert Gallee, LGGE, 28-03-2008)

ENDDO
RETURN
END SUBROUTINE SWDE