swde.F90

Go to the documentation of this file.

!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