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())
!        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