phylmd/Radlwsw/swclr.f

module swclr_m

  IMPLICIT NONE

contains

  SUBROUTINE swclr(knu, palbp, pdsig, prayl, psec, ppizaz, pray1, pray2, &
       prefz, prj, prk, prmu0, ptauaz, ptra1, ptra2)

    ! PURPOSE.
    ! COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY IN CASE OF
    ! CLEAR-SKY COLUMN

    ! REFERENCE.
    ! SEE RADIATION'S PART OF THE ECMWF RESEARCH DEPARTMENT
    ! DOCUMENTATION, AND FOUQUART AND BONNEL (1980)

    ! AUTHOR.
    ! JEAN-JACQUES MORCRETTE *ECMWF*

    ! MODIFICATIONS.
    ! ORIGINAL : 94-11-15

    USE raddim, only: kdlon, kflev
    USE radepsi, only: zepsec
    USE radopt, only: novlp

    ! ARGUMENTS:

    INTEGER knu
    DOUBLE PRECISION palbp(kdlon, 2)
    DOUBLE PRECISION, intent(in):: pdsig(kdlon, kflev)
    DOUBLE PRECISION, intent(in):: prayl(kdlon)
    DOUBLE PRECISION psec(kdlon)

    DOUBLE PRECISION, intent(out):: ppizaz(kdlon, kflev)
    DOUBLE PRECISION pray1(kdlon, kflev + 1)
    DOUBLE PRECISION pray2(kdlon, kflev + 1)
    DOUBLE PRECISION prefz(kdlon, 2, kflev + 1)
    DOUBLE PRECISION prj(kdlon, 6, kflev + 1)
    DOUBLE PRECISION prk(kdlon, 6, kflev + 1)
    DOUBLE PRECISION prmu0(kdlon, kflev + 1)
    DOUBLE PRECISION, intent(out):: ptauaz(kdlon, kflev)
    DOUBLE PRECISION ptra1(kdlon, kflev + 1)
    DOUBLE PRECISION ptra2(kdlon, kflev + 1)

    ! LOCAL VARIABLES:
    DOUBLE PRECISION zc0i(kdlon, kflev + 1)
    DOUBLE PRECISION zclear(kdlon)
    DOUBLE PRECISION zr21(kdlon)
    DOUBLE PRECISION zss0(kdlon)
    DOUBLE PRECISION zscat(kdlon)
    DOUBLE PRECISION ztr(kdlon, 2, kflev + 1)
    INTEGER jl, jk, ja, jkl, jklp1, jaj, jkm1
    DOUBLE PRECISION zfacoa, zcorae
    DOUBLE PRECISION zmue, zgap, zww, zto, zden, zmu1, zden1
    DOUBLE PRECISION zbmu0, zbmu1, zre11
    double precision, parameter:: REPSCT = 1d-10

    !------------------------------------------------------------------
    
    ! 1. OPTICAL PARAMETERS FOR AEROSOLS AND RAYLEIGH

    DO jk = 1, kflev + 1
       DO ja = 1, 6
          DO jl = 1, kdlon
             prj(jl, ja, jk) = 0d0
             prk(jl, ja, jk) = 0d0
          END DO
       END DO
    END DO

    DO jk = 1, kflev
       DO jl = 1, kdlon
          ptauaz(jl, jk) = prayl(jl) * pdsig(jl, jk)
          ppizaz(jl, jk) = 1d0 - repsct
       END DO
    END DO

    ! 2. TOTAL EFFECTIVE CLOUDINESS ABOVE A GIVEN LEVEL

    DO jl = 1, kdlon
       zc0i(jl, kflev + 1) = 0d0
       zclear(jl) = 1d0
       zscat(jl) = 0d0
    END DO

    jk = 1
    jkl = kflev + 1 - jk
    jklp1 = jkl + 1
    DO jl = 1, kdlon
       zfacoa = 1d0
       zcorae = zfacoa * ptauaz(jl, jkl) * psec(jl)
       zr21(jl) = exp(- zcorae)
       zss0(jl) = 1d0 - zr21(jl)

       IF (novlp == 1) THEN
          ! maximum-random
          zclear(jl) = zclear(jl) * (1d0 - max(zss0(jl), zscat(jl))) / (1d0 &
               - min(zscat(jl), 1d0 - zepsec))
          zc0i(jl, jkl) = 1d0 - zclear(jl)
          zscat(jl) = zss0(jl)
       ELSE IF (novlp == 2) THEN
          ! maximum
          zscat(jl) = max(zss0(jl), zscat(jl))
          zc0i(jl, jkl) = zscat(jl)
       ELSE IF (novlp == 3) THEN
          ! random
          zclear(jl) = zclear(jl) * (1d0 - zss0(jl))
          zscat(jl) = 1d0 - zclear(jl)
          zc0i(jl, jkl) = zscat(jl)
       END IF
    END DO

    DO jk = 2, kflev
       jkl = kflev + 1 - jk
       jklp1 = jkl + 1
       DO jl = 1, kdlon
          zfacoa = 1d0
          zcorae = zfacoa * ptauaz(jl, jkl) * psec(jl)
          zr21(jl) = exp(- zcorae)
          zss0(jl) = 1d0 - zr21(jl)

          IF (novlp == 1) THEN
             ! maximum-random
             zclear(jl) = zclear(jl) * (1d0 - max(zss0(jl), zscat(jl))) &
                  / (1d0 - min(zscat(jl), 1d0 - zepsec))
             zc0i(jl, jkl) = 1d0 - zclear(jl)
             zscat(jl) = zss0(jl)
          ELSE IF (novlp == 2) THEN
             ! maximum
             zscat(jl) = max(zss0(jl), zscat(jl))
             zc0i(jl, jkl) = zscat(jl)
          ELSE IF (novlp == 3) THEN
             ! random
             zclear(jl) = zclear(jl) * (1d0 - zss0(jl))
             zscat(jl) = 1d0 - zclear(jl)
             zc0i(jl, jkl) = zscat(jl)
          END IF
       END DO
    END DO

    ! 3. REFLECTIVITY/TRANSMISSIVITY FOR PURE SCATTERING

    DO jl = 1, kdlon
       pray1(jl, kflev + 1) = 0d0
       pray2(jl, kflev + 1) = 0d0
       prefz(jl, 2, 1) = palbp(jl, knu)
       prefz(jl, 1, 1) = palbp(jl, knu)
       ptra1(jl, kflev + 1) = 1d0
       ptra2(jl, kflev + 1) = 1d0
    END DO

    DO jk = 2, kflev + 1
       jkm1 = jk - 1
       DO jl = 1, kdlon

          ! 3.1 EQUIVALENT ZENITH ANGLE

          zmue = (1d0 - zc0i(jl, jk)) * psec(jl) + zc0i(jl, jk) * 1.66d0
          prmu0(jl, jk) = 1d0 / zmue

          ! 3.2 REFLECT./TRANSMISSIVITY DUE TO RAYLEIGH AND AEROSOLS

          zgap = 0d0
          zbmu0 = 0.5d0 - 0.75d0 * zgap / zmue
          zww = ppizaz(jl, jkm1)
          zto = ptauaz(jl, jkm1)
          zden = 1d0 + (1d0 - zww + zbmu0 * zww) * zto * zmue + (1d0 - zww) &
               * (1d0 - zww + 2d0 * zbmu0 * zww) * zto * zto * zmue * zmue
          pray1(jl, jkm1) = zbmu0 * zww * zto * zmue / zden
          ptra1(jl, jkm1) = 1d0 / zden

          zmu1 = 0.5d0
          zbmu1 = 0.5d0 - 0.75d0 * zgap * zmu1
          zden1 = 1d0 + (1d0 - zww + zbmu1 * zww) * zto / zmu1 + (1d0 - zww) &
               * (1d0 - zww + 2d0 * zbmu1 * zww &
               ) * zto * zto / zmu1 / zmu1
          pray2(jl, jkm1) = zbmu1 * zww * zto / zmu1 / zden1
          ptra2(jl, jkm1) = 1d0 / zden1

          prefz(jl, 1, jk) = (pray1(jl, jkm1) + prefz(jl, 1, jkm1) &
               * ptra1(jl, jkm1)* ptra2(jl, jkm1) / (1d0 - pray2(jl, jkm1) &
               * prefz(jl, 1, jkm1)))

          ztr(jl, 1, jkm1) = (ptra1(jl, jkm1) / (1d0 - pray2(jl, jkm1) &
               * prefz(jl, 1, jkm1)))

          prefz(jl, 2, jk) = (pray1(jl, jkm1) + prefz(jl, 2, jkm1) &
               * ptra1(jl, jkm1) * ptra2(jl, jkm1))

          ztr(jl, 2, jkm1) = ptra1(jl, jkm1)

       END DO
    END DO
    DO jl = 1, kdlon
       zmue = (1d0 - zc0i(jl, 1)) * psec(jl) + zc0i(jl, 1) * 1.66d0
       prmu0(jl, 1) = 1d0 / zmue
    END DO

    ! 3.5 REFLECT./TRANSMISSIVITY BETWEEN SURFACE AND LEVEL

    IF (knu == 1) THEN
       jaj = 2
       DO jl = 1, kdlon
          prj(jl, jaj, kflev + 1) = 1d0
          prk(jl, jaj, kflev + 1) = prefz(jl, 1, kflev + 1)
       END DO

       DO jk = 1, kflev
          jkl = kflev + 1 - jk
          jklp1 = jkl + 1
          DO jl = 1, kdlon
             zre11 = prj(jl, jaj, jklp1) * ztr(jl, 1, jkl)
             prj(jl, jaj, jkl) = zre11
             prk(jl, jaj, jkl) = zre11 * prefz(jl, 1, jkl)
          END DO
       END DO
    ELSE
       DO jaj = 1, 2
          DO jl = 1, kdlon
             prj(jl, jaj, kflev + 1) = 1d0
             prk(jl, jaj, kflev + 1) = prefz(jl, jaj, kflev + 1)
          END DO

          DO jk = 1, kflev
             jkl = kflev + 1 - jk
             jklp1 = jkl + 1
             DO jl = 1, kdlon
                zre11 = prj(jl, jaj, jklp1) * ztr(jl, jaj, jkl)
                prj(jl, jaj, jkl) = zre11
                prk(jl, jaj, jkl) = zre11 * prefz(jl, jaj, jkl)
             END DO
          END DO
       END DO
    END IF

  END SUBROUTINE swclr

end module swclr_m
1	guez	178	module swclr_m
2
3	guez	81	IMPLICIT NONE
4
5	guez	178	contains
6	guez	81
7	guez	219	SUBROUTINE swclr(knu, palbp, pdsig, prayl, psec, ppizaz, pray1, pray2, &
8			prefz, prj, prk, prmu0, ptauaz, ptra1, ptra2)
9	guez	81
10	guez	178	! PURPOSE.
11			! COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY IN CASE OF
12			! CLEAR-SKY COLUMN
13	guez	81
14	guez	178	! REFERENCE.
15			! SEE RADIATION'S PART OF THE ECMWF RESEARCH DEPARTMENT
16			! DOCUMENTATION, AND FOUQUART AND BONNEL (1980)
17	guez	81
18	guez	178	! AUTHOR.
19	guez	218	! JEAN-JACQUES MORCRETTE ECMWF
20	guez	81
21	guez	178	! MODIFICATIONS.
22			! ORIGINAL : 94-11-15
23	guez	81
24	guez	218	USE raddim, only: kdlon, kflev
25			USE radepsi, only: zepsec
26			USE radopt, only: novlp
27
28			! ARGUMENTS:
29
30	guez	178	INTEGER knu
31			DOUBLE PRECISION palbp(kdlon, 2)
32	guez	218	DOUBLE PRECISION, intent(in):: pdsig(kdlon, kflev)
33			DOUBLE PRECISION, intent(in):: prayl(kdlon)
34	guez	178	DOUBLE PRECISION psec(kdlon)
35	guez	81
36	guez	218	DOUBLE PRECISION, intent(out):: ppizaz(kdlon, kflev)
37			DOUBLE PRECISION pray1(kdlon, kflev + 1)
38			DOUBLE PRECISION pray2(kdlon, kflev + 1)
39			DOUBLE PRECISION prefz(kdlon, 2, kflev + 1)
40			DOUBLE PRECISION prj(kdlon, 6, kflev + 1)
41			DOUBLE PRECISION prk(kdlon, 6, kflev + 1)
42			DOUBLE PRECISION prmu0(kdlon, kflev + 1)
43			DOUBLE PRECISION, intent(out):: ptauaz(kdlon, kflev)
44			DOUBLE PRECISION ptra1(kdlon, kflev + 1)
45			DOUBLE PRECISION ptra2(kdlon, kflev + 1)
46	guez	81
47	guez	218	! LOCAL VARIABLES:
48			DOUBLE PRECISION zc0i(kdlon, kflev + 1)
49	guez	178	DOUBLE PRECISION zclear(kdlon)
50			DOUBLE PRECISION zr21(kdlon)
51			DOUBLE PRECISION zss0(kdlon)
52			DOUBLE PRECISION zscat(kdlon)
53	guez	218	DOUBLE PRECISION ztr(kdlon, 2, kflev + 1)
54	guez	178	INTEGER jl, jk, ja, jkl, jklp1, jaj, jkm1
55	guez	218	DOUBLE PRECISION zfacoa, zcorae
56	guez	178	DOUBLE PRECISION zmue, zgap, zww, zto, zden, zmu1, zden1
57			DOUBLE PRECISION zbmu0, zbmu1, zre11
58	guez	218	double precision, parameter:: REPSCT = 1d-10
59	guez	81
60	guez	218	!------------------------------------------------------------------
61
62			! 1. OPTICAL PARAMETERS FOR AEROSOLS AND RAYLEIGH
63	guez	81
64	guez	178	DO jk = 1, kflev + 1
65			DO ja = 1, 6
66			DO jl = 1, kdlon
67	guez	218	prj(jl, ja, jk) = 0d0
68			prk(jl, ja, jk) = 0d0
69	guez	178	END DO
70			END DO
71	guez	81	END DO
72
73	guez	178	DO jk = 1, kflev
74			DO jl = 1, kdlon
75	guez	218	ptauaz(jl, jk) = prayl(jl) * pdsig(jl, jk)
76	guez	219	ppizaz(jl, jk) = 1d0 - repsct
77	guez	178	END DO
78			END DO
79	guez	81
80	guez	218	! 2. TOTAL EFFECTIVE CLOUDINESS ABOVE A GIVEN LEVEL
81	guez	81
82	guez	178	DO jl = 1, kdlon
83	guez	218	zc0i(jl, kflev + 1) = 0d0
84			zclear(jl) = 1d0
85			zscat(jl) = 0d0
86	guez	178	END DO
87	guez	81
88	guez	178	jk = 1
89	guez	81	jkl = kflev + 1 - jk
90			jklp1 = jkl + 1
91			DO jl = 1, kdlon
92	guez	218	zfacoa = 1d0
93			zcorae = zfacoa * ptauaz(jl, jkl) * psec(jl)
94			zr21(jl) = exp(- zcorae)
95			zss0(jl) = 1d0 - zr21(jl)
96	guez	81
97	guez	218	IF (novlp == 1) THEN
98			! maximum-random
99			zclear(jl) = zclear(jl) * (1d0 - max(zss0(jl), zscat(jl))) / (1d0 &
100			- min(zscat(jl), 1d0 - zepsec))
101			zc0i(jl, jkl) = 1d0 - zclear(jl)
102	guez	178	zscat(jl) = zss0(jl)
103	guez	218	ELSE IF (novlp == 2) THEN
104			! maximum
105	guez	178	zscat(jl) = max(zss0(jl), zscat(jl))
106			zc0i(jl, jkl) = zscat(jl)
107	guez	218	ELSE IF (novlp == 3) THEN
108			! random
109			zclear(jl) = zclear(jl) * (1d0 - zss0(jl))
110			zscat(jl) = 1d0 - zclear(jl)
111	guez	178	zc0i(jl, jkl) = zscat(jl)
112			END IF
113	guez	81	END DO
114
115	guez	178	DO jk = 2, kflev
116			jkl = kflev + 1 - jk
117			jklp1 = jkl + 1
118			DO jl = 1, kdlon
119	guez	218	zfacoa = 1d0
120			zcorae = zfacoa * ptauaz(jl, jkl) * psec(jl)
121			zr21(jl) = exp(- zcorae)
122			zss0(jl) = 1d0 - zr21(jl)
123	guez	81
124	guez	218	IF (novlp == 1) THEN
125			! maximum-random
126			zclear(jl) = zclear(jl) * (1d0 - max(zss0(jl), zscat(jl))) &
127			/ (1d0 - min(zscat(jl), 1d0 - zepsec))
128			zc0i(jl, jkl) = 1d0 - zclear(jl)
129	guez	178	zscat(jl) = zss0(jl)
130	guez	218	ELSE IF (novlp == 2) THEN
131			! maximum
132	guez	178	zscat(jl) = max(zss0(jl), zscat(jl))
133			zc0i(jl, jkl) = zscat(jl)
134	guez	218	ELSE IF (novlp == 3) THEN
135			! random
136			zclear(jl) = zclear(jl) * (1d0 - zss0(jl))
137			zscat(jl) = 1d0 - zclear(jl)
138	guez	178	zc0i(jl, jkl) = zscat(jl)
139			END IF
140			END DO
141			END DO
142	guez	81
143	guez	218	! 3. REFLECTIVITY/TRANSMISSIVITY FOR PURE SCATTERING
144	guez	81
145			DO jl = 1, kdlon
146	guez	218	pray1(jl, kflev + 1) = 0d0
147			pray2(jl, kflev + 1) = 0d0
148	guez	178	prefz(jl, 2, 1) = palbp(jl, knu)
149			prefz(jl, 1, 1) = palbp(jl, knu)
150	guez	218	ptra1(jl, kflev + 1) = 1d0
151			ptra2(jl, kflev + 1) = 1d0
152	guez	178	END DO
153	guez	81
154	guez	178	DO jk = 2, kflev + 1
155			jkm1 = jk - 1
156			DO jl = 1, kdlon
157	guez	81
158	guez	218	! 3.1 EQUIVALENT ZENITH ANGLE
159	guez	81
160	guez	218	zmue = (1d0 - zc0i(jl, jk)) * psec(jl) + zc0i(jl, jk) * 1.66d0
161			prmu0(jl, jk) = 1d0 / zmue
162	guez	81
163	guez	218	! 3.2 REFLECT./TRANSMISSIVITY DUE TO RAYLEIGH AND AEROSOLS
164	guez	81
165	guez	218	zgap = 0d0
166			zbmu0 = 0.5d0 - 0.75d0 * zgap / zmue
167	guez	178	zww = ppizaz(jl, jkm1)
168			zto = ptauaz(jl, jkm1)
169	guez	218	zden = 1d0 + (1d0 - zww + zbmu0 * zww) * zto * zmue + (1d0 - zww) &
170			* (1d0 - zww + 2d0 * zbmu0 * zww) * zto * zto * zmue * zmue
171			pray1(jl, jkm1) = zbmu0 * zww * zto * zmue / zden
172			ptra1(jl, jkm1) = 1d0 / zden
173	guez	81
174	guez	218	zmu1 = 0.5d0
175			zbmu1 = 0.5d0 - 0.75d0 * zgap * zmu1
176			zden1 = 1d0 + (1d0 - zww + zbmu1 * zww) * zto / zmu1 + (1d0 - zww) &
177			* (1d0 - zww + 2d0 * zbmu1 * zww &
178			) * zto * zto / zmu1 / zmu1
179			pray2(jl, jkm1) = zbmu1 * zww * zto / zmu1 / zden1
180			ptra2(jl, jkm1) = 1d0 / zden1
181	guez	81
182	guez	218	prefz(jl, 1, jk) = (pray1(jl, jkm1) + prefz(jl, 1, jkm1) &
183			* ptra1(jl, jkm1)* ptra2(jl, jkm1) / (1d0 - pray2(jl, jkm1) &
184			* prefz(jl, 1, jkm1)))
185	guez	81
186	guez	218	ztr(jl, 1, jkm1) = (ptra1(jl, jkm1) / (1d0 - pray2(jl, jkm1) &
187			* prefz(jl, 1, jkm1)))
188	guez	81
189	guez	218	prefz(jl, 2, jk) = (pray1(jl, jkm1) + prefz(jl, 2, jkm1) &
190			* ptra1(jl, jkm1) * ptra2(jl, jkm1))
191	guez	81
192	guez	178	ztr(jl, 2, jkm1) = ptra1(jl, jkm1)
193
194			END DO
195	guez	81	END DO
196	guez	178	DO jl = 1, kdlon
197	guez	218	zmue = (1d0 - zc0i(jl, 1)) * psec(jl) + zc0i(jl, 1) * 1.66d0
198			prmu0(jl, 1) = 1d0 / zmue
199	guez	178	END DO
200	guez	81
201	guez	218	! 3.5 REFLECT./TRANSMISSIVITY BETWEEN SURFACE AND LEVEL
202	guez	81
203	guez	218	IF (knu == 1) THEN
204	guez	178	jaj = 2
205			DO jl = 1, kdlon
206	guez	218	prj(jl, jaj, kflev + 1) = 1d0
207			prk(jl, jaj, kflev + 1) = prefz(jl, 1, kflev + 1)
208	guez	178	END DO
209	guez	81
210	guez	178	DO jk = 1, kflev
211			jkl = kflev + 1 - jk
212			jklp1 = jkl + 1
213			DO jl = 1, kdlon
214	guez	218	zre11 = prj(jl, jaj, jklp1) * ztr(jl, 1, jkl)
215	guez	178	prj(jl, jaj, jkl) = zre11
216	guez	218	prk(jl, jaj, jkl) = zre11 * prefz(jl, 1, jkl)
217	guez	178	END DO
218			END DO
219			ELSE
220			DO jaj = 1, 2
221			DO jl = 1, kdlon
222	guez	218	prj(jl, jaj, kflev + 1) = 1d0
223			prk(jl, jaj, kflev + 1) = prefz(jl, jaj, kflev + 1)
224	guez	178	END DO
225	guez	81
226	guez	178	DO jk = 1, kflev
227			jkl = kflev + 1 - jk
228			jklp1 = jkl + 1
229			DO jl = 1, kdlon
230	guez	218	zre11 = prj(jl, jaj, jklp1) * ztr(jl, jaj, jkl)
231	guez	178	prj(jl, jaj, jkl) = zre11
232	guez	218	prk(jl, jaj, jkl) = zre11 * prefz(jl, jaj, jkl)
233	guez	178	END DO
234			END DO
235			END DO
236			END IF
237	guez	81
238	guez	178	END SUBROUTINE swclr
239	guez	81
240	guez	178	end module swclr_m