source: LMDZ6/trunk/libf/phylmd/rrtm/swde.F90 @ 3513

!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())
!        M.Hamrud      01-Oct-2003 CY28 Cleaning
!        Modified: 03-10-10 Deborah Salmond and Marta Janiskova Optimisation
!        Modified: 03-12-13 John Hague - MASS Vector Fns
!        Y.Seity 06-09-09 : add modset from O.Thouron (MesoNH) under NOVLP tests
!     ------------------------------------------------------------------

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

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

USE PARKIND1  ,ONLY : JPIM     ,JPRB
USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK

USE YOERDU   , ONLY : REPLOG
USE YOMJFH   , ONLY : N_VMASS
!++MODIFCODE
USE YOERAD   , ONLY : NOVLP 
!--MODIFCODE
IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN)    :: KLON 
INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PGG(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PREF(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PRMUZ(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTO1(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PW(KLON) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PRE1(KLON) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PRE2(KLON) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTR1(KLON) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTR2(KLON) 
REAL(KIND=JPRB) :: ZTMP   (4,KFDIA-KIDIA+1)
REAL(KIND=JPRB) :: ZTMP2  (KFDIA-KIDIA+1+N_VMASS)
REAL(KIND=JPRB) :: ZTMP3  (KFDIA-KIDIA+1+N_VMASS)
REAL(KIND=JPRB) :: ZZARG  (KFDIA-KIDIA+1+N_VMASS)
REAL(KIND=JPRB) :: ZZARG2 (KFDIA-KIDIA+1+N_VMASS)

INTEGER(KIND=JPIM) :: JL, JLL, JLEN

REAL(KIND=JPRB) :: ZA11, ZA12, ZA13, ZA21, ZA22, ZA23, ZALPHA,&
 & ZAM2B, ZAP2B,  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  


REAL(KIND=JPRB) :: MINJ, MAXJ, X, Y
REAL(KIND=JPRB) :: ZPRMUZ,ZIDENA,ZIDENB,ZRR
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!     STATEMENT DUNCTIONS
MINJ(X,Y) = Y - 0.5_JPRB*(ABS(X-Y)-(X-Y))
MAXJ(X,Y) = Y + 0.5_JPRB*(ABS(X-Y)+(X-Y))

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

!*         1.      DELTA-EDDINGTON CALCULATIONS

IF (LHOOK) CALL DR_HOOK('SWDE',0,ZHOOK_HANDLE)

ZDT = 2.0_JPRB/3._JPRB

  DO JL   =   KIDIA,KFDIA
    JLL=JL-KIDIA+1
    ZPRMUZ=1.0_JPRB/PRMUZ(JL)
!++MODIFCODE
  IF (NOVLP >= 5) THEN !MESONH_VERSION
   ZGP = PGG(JL)
   ZTOP = PTO1(JL)
   ZWCP = PW(JL)  
  ELSE !ECMWF VERSION
    ZFF = PGG(JL)*PGG(JL)
    ZGP = PGG(JL)/(1.0_JPRB+PGG(JL))
    ZTOP = (1.0_JPRB- PW(JL) * ZFF) * PTO1(JL)
    ZWCP = (1-ZFF)* PW(JL) /(1.0_JPRB- PW(JL) * ZFF)
  ENDIF
!--MODIFCODE
    ZX1 = 1.0_JPRB-ZWCP*ZGP
    ZWM = 1.0_JPRB-ZWCP
    ZRM2 =  PRMUZ(JL) * PRMUZ(JL)
    ZRK = SQRT(MAXJ(REPLOG,3._JPRB*ZWM*ZX1))
    ZX2 = (1.0_JPRB-ZRK*ZRK*ZRM2)*(4._JPRB/3._JPRB)
    ZRR = 1.0_JPRB/ZX2
    ZRP=ZRK/ZX1
    ZALPHA = ZWCP*ZRM2*(1.0_JPRB+ZGP*ZWM)*ZRR 
    ZBETA = ZWCP* PRMUZ(JL) *(1.0_JPRB+3._JPRB*ZGP*ZRM2*ZWM)*ZRR 
    ZZARG(JLL)   = -MAXJ( -200._JPRB, MINJ( ZTOP*ZPRMUZ, 200._JPRB) )
    ZZARG2(JLL)  = MINJ( ZRK*ZTOP, 200._JPRB)
    ZTMP(1,JLL) = ZPRMUZ
    ZTMP(2,JLL) = ZALPHA
    ZTMP(3,JLL) = ZBETA
    ZTMP(4,JLL) = ZRP
  ENDDO

  IF(N_VMASS /= 0 ) THEN  !USING VECTOR MASS
    JLEN=KFDIA-KIDIA+N_VMASS-MOD(KFDIA-KIDIA,N_VMASS)
    IF(KFDIA-KIDIA+1 /= JLEN) THEN
      ZZARG  (KFDIA-KIDIA+2:JLEN)=1.0_JPRB
      ZZARG2 (KFDIA-KIDIA+2:JLEN)=1.0_JPRB
    ENDIF
! Commente par MPL le 21.11.08
!   CALL VEXP(ZTMP2,ZZARG, JLEN)
!   CALL VEXP(ZTMP3,ZZARG2,JLEN)
  ELSE
    DO JL   =   KIDIA,KFDIA
      JLL=JL-KIDIA+1
      ZTMP2(JLL) = EXP(ZZARG(JLL))
      ZTMP3(JLL) = EXP(ZZARG2(JLL))
    ENDDO
  ENDIF

  DO JL   =   KIDIA,KFDIA
    JLL=JL-KIDIA+1
    ZEXMU0 = ZTMP2(JLL)
    ZEXKP  = ZTMP3(JLL)
    ZPRMUZ = ZTMP(1,JLL)
    ZALPHA = ZTMP(2,JLL)
    ZBETA  = ZTMP(3,JLL)
    ZRP    = ZTMP(4,JLL)
    ZEXKM = 1.0_JPRB/ZEXKP
    ZXP2P = 1.0_JPRB+ZDT*ZRP
    ZXM2P = 1.0_JPRB-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
    ZIDENA=1.0_JPRB/ZDENA
    ZC1A = (ZA22*ZA13-ZA12*ZA23)*ZIDENA
    ZC2A = (ZA11*ZA23-ZA21*ZA13)*ZIDENA
    ZRI0A = ZC1A+ZC2A-ZALPHA
    ZRI1A = ZRP*(ZC1A-ZC2A)-ZBETA
    PRE1(JL) = (ZRI0A-ZDT*ZRI1A)*ZPRMUZ
    ZRI0B = ZC1A*ZEXKM+ZC2A*ZEXKP-ZALPHA*ZEXMU0
    ZRI1B = ZRP*(ZC1A*ZEXKM-ZC2A*ZEXKP)-ZBETA*ZEXMU0
    PTR1(JL) = ZEXMU0+(ZRI0B+ZDT*ZRI1B)*ZPRMUZ

!*         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
    ZIDENB= 1.0_JPRB/ZDENB
    ZC1B = (ZB22*ZA13-ZA12*ZB23)*ZIDENB
    ZC2B = (ZA11*ZB23-ZB21*ZA13)*ZIDENB
    ZRI0C = ZC1B+ZC2B-ZALPHA
    ZRI1C = ZRP*(ZC1B-ZC2B)-ZBETA
    PRE2(JL) = (ZRI0C-ZDT*ZRI1C) * ZPRMUZ
    ZRI0D = ZC1B*ZEXKM + ZC2B*ZEXKP - ZALPHA*ZEXMU0
    ZRI1D = ZRP * (ZC1B*ZEXKM - ZC2B*ZEXKP) - ZBETA*ZEXMU0
    PTR2(JL) = ZEXMU0 + (ZRI0D + ZDT*ZRI1D) * ZPRMUZ
  ENDDO

IF (LHOOK) CALL DR_HOOK('SWDE',1,ZHOOK_HANDLE)
END SUBROUTINE SWDE