phylmd/Radlwsw/swclr.f

module swclr_m

  IMPLICIT NONE

contains

  SUBROUTINE swclr(knu, flag_aer, palbp, pdsig, prayl, psec, pcgaz, 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
    ! -OB
    logical, intent(in):: flag_aer
    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):: pcgaz(kdlon, kflev)
    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
          pcgaz(jl, jk) = 0d0
          ptauaz(jl, jk) = prayl(jl) * pdsig(jl, jk)
       END DO

       IF (flag_aer) THEN
          ! -OB
          DO jl = 1, kdlon
             ppizaz(jl, jk) = 1d0
          END DO
       ELSE
          DO jl = 1, kdlon
             ppizaz(jl, jk) = 1d0 - repsct
          END DO
       END IF
    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	217	SUBROUTINE swclr(knu, flag_aer, palbp, pdsig, prayl, psec, pcgaz, ppizaz, &
8			pray1, pray2, 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			! -OB
32	guez	217	logical, intent(in):: flag_aer
33	guez	178	DOUBLE PRECISION palbp(kdlon, 2)
34	guez	218	DOUBLE PRECISION, intent(in):: pdsig(kdlon, kflev)
35			DOUBLE PRECISION, intent(in):: prayl(kdlon)
36	guez	178	DOUBLE PRECISION psec(kdlon)
37	guez	81
38	guez	218	DOUBLE PRECISION, intent(out):: pcgaz(kdlon, kflev)
39			DOUBLE PRECISION, intent(out):: ppizaz(kdlon, kflev)
40			DOUBLE PRECISION pray1(kdlon, kflev + 1)
41			DOUBLE PRECISION pray2(kdlon, kflev + 1)
42			DOUBLE PRECISION prefz(kdlon, 2, kflev + 1)
43			DOUBLE PRECISION prj(kdlon, 6, kflev + 1)
44			DOUBLE PRECISION prk(kdlon, 6, kflev + 1)
45			DOUBLE PRECISION prmu0(kdlon, kflev + 1)
46			DOUBLE PRECISION, intent(out):: ptauaz(kdlon, kflev)
47			DOUBLE PRECISION ptra1(kdlon, kflev + 1)
48			DOUBLE PRECISION ptra2(kdlon, kflev + 1)
49	guez	81
50	guez	218	! LOCAL VARIABLES:
51			DOUBLE PRECISION zc0i(kdlon, kflev + 1)
52	guez	178	DOUBLE PRECISION zclear(kdlon)
53			DOUBLE PRECISION zr21(kdlon)
54			DOUBLE PRECISION zss0(kdlon)
55			DOUBLE PRECISION zscat(kdlon)
56	guez	218	DOUBLE PRECISION ztr(kdlon, 2, kflev + 1)
57	guez	178	INTEGER jl, jk, ja, jkl, jklp1, jaj, jkm1
58	guez	218	DOUBLE PRECISION zfacoa, zcorae
59	guez	178	DOUBLE PRECISION zmue, zgap, zww, zto, zden, zmu1, zden1
60			DOUBLE PRECISION zbmu0, zbmu1, zre11
61	guez	218	double precision, parameter:: REPSCT = 1d-10
62	guez	81
63	guez	218	!------------------------------------------------------------------
64
65			! 1. OPTICAL PARAMETERS FOR AEROSOLS AND RAYLEIGH
66	guez	81
67	guez	178	DO jk = 1, kflev + 1
68			DO ja = 1, 6
69			DO jl = 1, kdlon
70	guez	218	prj(jl, ja, jk) = 0d0
71			prk(jl, ja, jk) = 0d0
72	guez	178	END DO
73			END DO
74	guez	81	END DO
75
76	guez	178	DO jk = 1, kflev
77			DO jl = 1, kdlon
78	guez	217	pcgaz(jl, jk) = 0d0
79	guez	218	ptauaz(jl, jk) = prayl(jl) * pdsig(jl, jk)
80	guez	178	END DO
81	guez	81
82	guez	217	IF (flag_aer) THEN
83	guez	178	! -OB
84			DO jl = 1, kdlon
85	guez	218	ppizaz(jl, jk) = 1d0
86	guez	178	END DO
87			ELSE
88			DO jl = 1, kdlon
89	guez	218	ppizaz(jl, jk) = 1d0 - repsct
90	guez	178	END DO
91	guez	217	END IF
92	guez	178	END DO
93	guez	81
94	guez	218	! 2. TOTAL EFFECTIVE CLOUDINESS ABOVE A GIVEN LEVEL
95	guez	81
96	guez	178	DO jl = 1, kdlon
97	guez	218	zc0i(jl, kflev + 1) = 0d0
98			zclear(jl) = 1d0
99			zscat(jl) = 0d0
100	guez	178	END DO
101	guez	81
102	guez	178	jk = 1
103	guez	81	jkl = kflev + 1 - jk
104			jklp1 = jkl + 1
105			DO jl = 1, kdlon
106	guez	218	zfacoa = 1d0
107			zcorae = zfacoa * ptauaz(jl, jkl) * psec(jl)
108			zr21(jl) = exp(- zcorae)
109			zss0(jl) = 1d0 - zr21(jl)
110	guez	81
111	guez	218	IF (novlp == 1) THEN
112			! maximum-random
113			zclear(jl) = zclear(jl) * (1d0 - max(zss0(jl), zscat(jl))) / (1d0 &
114			- min(zscat(jl), 1d0 - zepsec))
115			zc0i(jl, jkl) = 1d0 - zclear(jl)
116	guez	178	zscat(jl) = zss0(jl)
117	guez	218	ELSE IF (novlp == 2) THEN
118			! maximum
119	guez	178	zscat(jl) = max(zss0(jl), zscat(jl))
120			zc0i(jl, jkl) = zscat(jl)
121	guez	218	ELSE IF (novlp == 3) THEN
122			! random
123			zclear(jl) = zclear(jl) * (1d0 - zss0(jl))
124			zscat(jl) = 1d0 - zclear(jl)
125	guez	178	zc0i(jl, jkl) = zscat(jl)
126			END IF
127	guez	81	END DO
128
129	guez	178	DO jk = 2, kflev
130			jkl = kflev + 1 - jk
131			jklp1 = jkl + 1
132			DO jl = 1, kdlon
133	guez	218	zfacoa = 1d0
134			zcorae = zfacoa * ptauaz(jl, jkl) * psec(jl)
135			zr21(jl) = exp(- zcorae)
136			zss0(jl) = 1d0 - zr21(jl)
137	guez	81
138	guez	218	IF (novlp == 1) THEN
139			! maximum-random
140			zclear(jl) = zclear(jl) * (1d0 - max(zss0(jl), zscat(jl))) &
141			/ (1d0 - min(zscat(jl), 1d0 - zepsec))
142			zc0i(jl, jkl) = 1d0 - zclear(jl)
143	guez	178	zscat(jl) = zss0(jl)
144	guez	218	ELSE IF (novlp == 2) THEN
145			! maximum
146	guez	178	zscat(jl) = max(zss0(jl), zscat(jl))
147			zc0i(jl, jkl) = zscat(jl)
148	guez	218	ELSE IF (novlp == 3) THEN
149			! random
150			zclear(jl) = zclear(jl) * (1d0 - zss0(jl))
151			zscat(jl) = 1d0 - zclear(jl)
152	guez	178	zc0i(jl, jkl) = zscat(jl)
153			END IF
154			END DO
155			END DO
156	guez	81
157	guez	218	! 3. REFLECTIVITY/TRANSMISSIVITY FOR PURE SCATTERING
158	guez	81
159			DO jl = 1, kdlon
160	guez	218	pray1(jl, kflev + 1) = 0d0
161			pray2(jl, kflev + 1) = 0d0
162	guez	178	prefz(jl, 2, 1) = palbp(jl, knu)
163			prefz(jl, 1, 1) = palbp(jl, knu)
164	guez	218	ptra1(jl, kflev + 1) = 1d0
165			ptra2(jl, kflev + 1) = 1d0
166	guez	178	END DO
167	guez	81
168	guez	178	DO jk = 2, kflev + 1
169			jkm1 = jk - 1
170			DO jl = 1, kdlon
171	guez	81
172	guez	218	! 3.1 EQUIVALENT ZENITH ANGLE
173	guez	81
174	guez	218	zmue = (1d0 - zc0i(jl, jk)) * psec(jl) + zc0i(jl, jk) * 1.66d0
175			prmu0(jl, jk) = 1d0 / zmue
176	guez	81
177	guez	218	! 3.2 REFLECT./TRANSMISSIVITY DUE TO RAYLEIGH AND AEROSOLS
178	guez	81
179	guez	218	zgap = 0d0
180			zbmu0 = 0.5d0 - 0.75d0 * zgap / zmue
181	guez	178	zww = ppizaz(jl, jkm1)
182			zto = ptauaz(jl, jkm1)
183	guez	218	zden = 1d0 + (1d0 - zww + zbmu0 * zww) * zto * zmue + (1d0 - zww) &
184			* (1d0 - zww + 2d0 * zbmu0 * zww) * zto * zto * zmue * zmue
185			pray1(jl, jkm1) = zbmu0 * zww * zto * zmue / zden
186			ptra1(jl, jkm1) = 1d0 / zden
187	guez	81
188	guez	218	zmu1 = 0.5d0
189			zbmu1 = 0.5d0 - 0.75d0 * zgap * zmu1
190			zden1 = 1d0 + (1d0 - zww + zbmu1 * zww) * zto / zmu1 + (1d0 - zww) &
191			* (1d0 - zww + 2d0 * zbmu1 * zww &
192			) * zto * zto / zmu1 / zmu1
193			pray2(jl, jkm1) = zbmu1 * zww * zto / zmu1 / zden1
194			ptra2(jl, jkm1) = 1d0 / zden1
195	guez	81
196	guez	218	prefz(jl, 1, jk) = (pray1(jl, jkm1) + prefz(jl, 1, jkm1) &
197			* ptra1(jl, jkm1)* ptra2(jl, jkm1) / (1d0 - pray2(jl, jkm1) &
198			* prefz(jl, 1, jkm1)))
199	guez	81
200	guez	218	ztr(jl, 1, jkm1) = (ptra1(jl, jkm1) / (1d0 - pray2(jl, jkm1) &
201			* prefz(jl, 1, jkm1)))
202	guez	81
203	guez	218	prefz(jl, 2, jk) = (pray1(jl, jkm1) + prefz(jl, 2, jkm1) &
204			* ptra1(jl, jkm1) * ptra2(jl, jkm1))
205	guez	81
206	guez	178	ztr(jl, 2, jkm1) = ptra1(jl, jkm1)
207
208			END DO
209	guez	81	END DO
210	guez	178	DO jl = 1, kdlon
211	guez	218	zmue = (1d0 - zc0i(jl, 1)) * psec(jl) + zc0i(jl, 1) * 1.66d0
212			prmu0(jl, 1) = 1d0 / zmue
213	guez	178	END DO
214	guez	81
215	guez	218	! 3.5 REFLECT./TRANSMISSIVITY BETWEEN SURFACE AND LEVEL
216	guez	81
217	guez	218	IF (knu == 1) THEN
218	guez	178	jaj = 2
219			DO jl = 1, kdlon
220	guez	218	prj(jl, jaj, kflev + 1) = 1d0
221			prk(jl, jaj, kflev + 1) = prefz(jl, 1, kflev + 1)
222	guez	178	END DO
223	guez	81
224	guez	178	DO jk = 1, kflev
225			jkl = kflev + 1 - jk
226			jklp1 = jkl + 1
227			DO jl = 1, kdlon
228	guez	218	zre11 = prj(jl, jaj, jklp1) * ztr(jl, 1, jkl)
229	guez	178	prj(jl, jaj, jkl) = zre11
230	guez	218	prk(jl, jaj, jkl) = zre11 * prefz(jl, 1, jkl)
231	guez	178	END DO
232			END DO
233			ELSE
234			DO jaj = 1, 2
235			DO jl = 1, kdlon
236	guez	218	prj(jl, jaj, kflev + 1) = 1d0
237			prk(jl, jaj, kflev + 1) = prefz(jl, jaj, kflev + 1)
238	guez	178	END DO
239	guez	81
240	guez	178	DO jk = 1, kflev
241			jkl = kflev + 1 - jk
242			jklp1 = jkl + 1
243			DO jl = 1, kdlon
244	guez	218	zre11 = prj(jl, jaj, jklp1) * ztr(jl, jaj, jkl)
245	guez	178	prj(jl, jaj, jkl) = zre11
246	guez	218	prk(jl, jaj, jkl) = zre11 * prefz(jl, jaj, jkl)
247	guez	178	END DO
248			END DO
249			END DO
250			END IF
251	guez	81
252	guez	178	END SUBROUTINE swclr
253	guez	81
254	guez	178	end module swclr_m