phylmd/Radlwsw/swde.f90

      SUBROUTINE SWDE (PGG,PREF,PRMUZ,PTO1,PW,
     S                 PRE1,PRE2,PTR1,PTR2)
      use dimens_m
      use dimphy
      use raddim
      IMPLICIT none
C
C     ------------------------------------------------------------------
C     PURPOSE.
C     --------
C           COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY OF A CLOUDY
C     LAYER USING THE DELTA-EDDINGTON'S APPROXIMATION.
C
C     METHOD.
C     -------
C
C          STANDARD DELTA-EDDINGTON LAYER CALCULATIONS.
C
C     REFERENCE.
C     ----------
C
C        SEE RADIATION'S PART OF THE MODEL'S DOCUMENTATION AND
C        ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE IFS
C
C     AUTHOR.
C     -------
C        JEAN-JACQUES MORCRETTE  *ECMWF*
C
C     MODIFICATIONS.
C     --------------
C        ORIGINAL : 88-12-15
C     ------------------------------------------------------------------
C* ARGUMENTS:
C
      DOUBLE PRECISION PGG(KDLON)   ! ASSYMETRY FACTOR
      DOUBLE PRECISION PREF(KDLON)  ! REFLECTIVITY OF THE UNDERLYING LAYER
      DOUBLE PRECISION PRMUZ(KDLON) ! COSINE OF SOLAR ZENITH ANGLE
      DOUBLE PRECISION PTO1(KDLON)  ! OPTICAL THICKNESS
      DOUBLE PRECISION PW(KDLON)    ! SINGLE SCATTERING ALBEDO
      DOUBLE PRECISION PRE1(KDLON)  ! LAYER REFLECTIVITY (NO UNDERLYING-LAYER REFLECTION)
      DOUBLE PRECISION PRE2(KDLON)  ! LAYER REFLECTIVITY
      DOUBLE PRECISION PTR1(KDLON)  ! LAYER TRANSMISSIVITY (NO UNDERLYING-LAYER REFLECTION)
      DOUBLE PRECISION PTR2(KDLON)  ! LAYER TRANSMISSIVITY
C
C* LOCAL VARIABLES:
C
      INTEGER jl
      DOUBLE PRECISION ZFF, ZGP, ZTOP, ZWCP, ZDT, ZX1, ZWM
      DOUBLE PRECISION ZRM2, ZRK, ZX2, ZRP, ZALPHA, ZBETA, ZARG
      DOUBLE PRECISION ZEXMU0, ZARG2, ZEXKP, ZEXKM, ZXP2P, ZXM2P, ZAP2B
      DOUBLE PRECISION ZAM2B
      DOUBLE PRECISION ZA11, ZA12, ZA13, ZA21, ZA22, ZA23
      DOUBLE PRECISION ZDENA, ZC1A, ZC2A, ZRI0A, ZRI1A
      DOUBLE PRECISION ZRI0B, ZRI1B
      DOUBLE PRECISION ZB21, ZB22, ZB23, ZDENB, ZC1B, ZC2B
      DOUBLE PRECISION ZRI0C, ZRI1C, ZRI0D, ZRI1D
C     ------------------------------------------------------------------
C
C*         1.      DELTA-EDDINGTON CALCULATIONS
C
 100  CONTINUE
C
      DO 131 JL   =   1, KDLON
C
C*         1.1     SET UP THE DELTA-MODIFIED PARAMETERS
C
 110  CONTINUE
C
      ZFF = PGG(JL)*PGG(JL)
      ZGP = PGG(JL)/(1.+PGG(JL))
      ZTOP = (1.- PW(JL) * ZFF) * PTO1(JL)
      ZWCP = (1-ZFF)* PW(JL) /(1.- PW(JL) * ZFF)
      ZDT = 2./3.
      ZX1 = 1.-ZWCP*ZGP
      ZWM = 1.-ZWCP
      ZRM2 =  PRMUZ(JL) * PRMUZ(JL)
      ZRK = SQRT(3.*ZWM*ZX1)
      ZX2 = 4.*(1.-ZRK*ZRK*ZRM2)
      ZRP=ZRK/ZX1
      ZALPHA = 3.*ZWCP*ZRM2*(1.+ZGP*ZWM)/ZX2
      ZBETA = 3.*ZWCP* PRMUZ(JL) *(1.+3.*ZGP*ZRM2*ZWM)/ZX2
CMAF      ZARG=MIN(ZTOP/PRMUZ(JL),200.)
      ZARG=MIN(ZTOP/PRMUZ(JL),2.0d+2)
      ZEXMU0=EXP(-ZARG)
CMAF      ZARG2=MIN(ZRK*ZTOP,200.)
      ZARG2=MIN(ZRK*ZTOP,2.0d+2)
      ZEXKP=EXP(ZARG2)
      ZEXKM = 1./ZEXKP
      ZXP2P = 1.+ZDT*ZRP
      ZXM2P = 1.-ZDT*ZRP
      ZAP2B = ZALPHA+ZDT*ZBETA
      ZAM2B = ZALPHA-ZDT*ZBETA
C
C*         1.2     WITHOUT REFLECTION FROM THE UNDERLYING LAYER
C
 120  CONTINUE
C
      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)
C
C*         1.3     WITH REFLECTION FROM THE UNDERLYING LAYER
C
 130  CONTINUE
C
      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)
C
 131  CONTINUE
      RETURN
      END
1	guez	24	SUBROUTINE SWDE (PGG,PREF,PRMUZ,PTO1,PW,
2			S PRE1,PRE2,PTR1,PTR2)
3			use dimens_m
4			use dimphy
5			use raddim
6			IMPLICIT none
7			C
8			C ------------------------------------------------------------------
9			C PURPOSE.
10			C --------
11			C COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY OF A CLOUDY
12			C LAYER USING THE DELTA-EDDINGTON'S APPROXIMATION.
13			C
14			C METHOD.
15			C -------
16			C
17			C STANDARD DELTA-EDDINGTON LAYER CALCULATIONS.
18			C
19			C REFERENCE.
20			C ----------
21			C
22			C SEE RADIATION'S PART OF THE MODEL'S DOCUMENTATION AND
23			C ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE IFS
24			C
25			C AUTHOR.
26			C -------
27			C JEAN-JACQUES MORCRETTE ECMWF
28			C
29			C MODIFICATIONS.
30			C --------------
31			C ORIGINAL : 88-12-15
32			C ------------------------------------------------------------------
33			C* ARGUMENTS:
34			C
35	guez	71	DOUBLE PRECISION PGG(KDLON) ! ASSYMETRY FACTOR
36			DOUBLE PRECISION PREF(KDLON) ! REFLECTIVITY OF THE UNDERLYING LAYER
37			DOUBLE PRECISION PRMUZ(KDLON) ! COSINE OF SOLAR ZENITH ANGLE
38			DOUBLE PRECISION PTO1(KDLON) ! OPTICAL THICKNESS
39			DOUBLE PRECISION PW(KDLON) ! SINGLE SCATTERING ALBEDO
40			DOUBLE PRECISION PRE1(KDLON) ! LAYER REFLECTIVITY (NO UNDERLYING-LAYER REFLECTION)
41			DOUBLE PRECISION PRE2(KDLON) ! LAYER REFLECTIVITY
42			DOUBLE PRECISION PTR1(KDLON) ! LAYER TRANSMISSIVITY (NO UNDERLYING-LAYER REFLECTION)
43			DOUBLE PRECISION PTR2(KDLON) ! LAYER TRANSMISSIVITY
44	guez	24	C
45			C* LOCAL VARIABLES:
46			C
47			INTEGER jl
48	guez	71	DOUBLE PRECISION ZFF, ZGP, ZTOP, ZWCP, ZDT, ZX1, ZWM
49			DOUBLE PRECISION ZRM2, ZRK, ZX2, ZRP, ZALPHA, ZBETA, ZARG
50			DOUBLE PRECISION ZEXMU0, ZARG2, ZEXKP, ZEXKM, ZXP2P, ZXM2P, ZAP2B
51			DOUBLE PRECISION ZAM2B
52			DOUBLE PRECISION ZA11, ZA12, ZA13, ZA21, ZA22, ZA23
53			DOUBLE PRECISION ZDENA, ZC1A, ZC2A, ZRI0A, ZRI1A
54			DOUBLE PRECISION ZRI0B, ZRI1B
55			DOUBLE PRECISION ZB21, ZB22, ZB23, ZDENB, ZC1B, ZC2B
56			DOUBLE PRECISION ZRI0C, ZRI1C, ZRI0D, ZRI1D
57	guez	24	C ------------------------------------------------------------------
58			C
59			C* 1. DELTA-EDDINGTON CALCULATIONS
60			C
61			100 CONTINUE
62			C
63			DO 131 JL = 1, KDLON
64			C
65			C* 1.1 SET UP THE DELTA-MODIFIED PARAMETERS
66			C
67			110 CONTINUE
68			C
69			ZFF = PGG(JL)*PGG(JL)
70			ZGP = PGG(JL)/(1.+PGG(JL))
71			ZTOP = (1.- PW(JL) * ZFF) * PTO1(JL)
72			ZWCP = (1-ZFF)* PW(JL) /(1.- PW(JL) * ZFF)
73			ZDT = 2./3.
74			ZX1 = 1.-ZWCP*ZGP
75			ZWM = 1.-ZWCP
76			ZRM2 = PRMUZ(JL) * PRMUZ(JL)
77			ZRK = SQRT(3.ZWMZX1)
78			ZX2 = 4.(1.-ZRKZRK*ZRM2)
79			ZRP=ZRK/ZX1
80			ZALPHA = 3.ZWCPZRM2(1.+ZGPZWM)/ZX2
81			ZBETA = 3.ZWCP PRMUZ(JL) (1.+3.ZGPZRM2ZWM)/ZX2
82			CMAF ZARG=MIN(ZTOP/PRMUZ(JL),200.)
83			ZARG=MIN(ZTOP/PRMUZ(JL),2.0d+2)
84			ZEXMU0=EXP(-ZARG)
85			CMAF ZARG2=MIN(ZRK*ZTOP,200.)
86			ZARG2=MIN(ZRK*ZTOP,2.0d+2)
87			ZEXKP=EXP(ZARG2)
88			ZEXKM = 1./ZEXKP
89			ZXP2P = 1.+ZDT*ZRP
90			ZXM2P = 1.-ZDT*ZRP
91			ZAP2B = ZALPHA+ZDT*ZBETA
92			ZAM2B = ZALPHA-ZDT*ZBETA
93			C
94			C* 1.2 WITHOUT REFLECTION FROM THE UNDERLYING LAYER
95			C
96			120 CONTINUE
97			C
98			ZA11 = ZXP2P
99			ZA12 = ZXM2P
100			ZA13 = ZAP2B
101			ZA22 = ZXP2P*ZEXKP
102			ZA21 = ZXM2P*ZEXKM
103			ZA23 = ZAM2B*ZEXMU0
104			ZDENA = ZA11 * ZA22 - ZA21 * ZA12
105			ZC1A = (ZA22ZA13-ZA12ZA23)/ZDENA
106			ZC2A = (ZA11ZA23-ZA21ZA13)/ZDENA
107			ZRI0A = ZC1A+ZC2A-ZALPHA
108			ZRI1A = ZRP*(ZC1A-ZC2A)-ZBETA
109			PRE1(JL) = (ZRI0A-ZDT*ZRI1A)/ PRMUZ(JL)
110			ZRI0B = ZC1AZEXKM+ZC2AZEXKP-ZALPHA*ZEXMU0
111			ZRI1B = ZRP(ZC1AZEXKM-ZC2AZEXKP)-ZBETAZEXMU0
112			PTR1(JL) = ZEXMU0+(ZRI0B+ZDT*ZRI1B)/ PRMUZ(JL)
113			C
114			C* 1.3 WITH REFLECTION FROM THE UNDERLYING LAYER
115			C
116			130 CONTINUE
117			C
118			ZB21 = ZA21- PREF(JL) ZXP2PZEXKM
119			ZB22 = ZA22- PREF(JL) ZXM2PZEXKP
120			ZB23 = ZA23- PREF(JL) ZEXMU0(ZAP2B - PRMUZ(JL) )
121			ZDENB = ZA11 * ZB22 - ZB21 * ZA12
122			ZC1B = (ZB22ZA13-ZA12ZB23)/ZDENB
123			ZC2B = (ZA11ZB23-ZB21ZA13)/ZDENB
124			ZRI0C = ZC1B+ZC2B-ZALPHA
125			ZRI1C = ZRP*(ZC1B-ZC2B)-ZBETA
126			PRE2(JL) = (ZRI0C-ZDT*ZRI1C) / PRMUZ(JL)
127			ZRI0D = ZC1BZEXKM + ZC2BZEXKP - ZALPHA*ZEXMU0
128			ZRI1D = ZRP * (ZC1BZEXKM - ZC2BZEXKP) - ZBETA*ZEXMU0
129			PTR2(JL) = ZEXMU0 + (ZRI0D + ZDT*ZRI1D) / PRMUZ(JL)
130			C
131			131 CONTINUE
132			RETURN
133			END