phylmd/Radlwsw/swde.f

      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	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	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	C
45	C* LOCAL VARIABLES:
46	C
47	INTEGER jl
48	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	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