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)
    
    USE raddim, only: kdlon, kflev
    USE radepsi, only: repsct, zepsec
    USE radopt, only: novlp

    ! ------------------------------------------------------------------
    ! 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
    ! ------------------------------------------------------------------
    ! * ARGUMENTS:

    INTEGER knu
    ! -OB
    logical, intent(in):: flag_aer
    DOUBLE PRECISION palbp(kdlon, 2)
    DOUBLE PRECISION pdsig(kdlon, kflev)
    DOUBLE PRECISION prayl(kdlon)
    DOUBLE PRECISION psec(kdlon)

    DOUBLE PRECISION pcgaz(kdlon, kflev)
    DOUBLE PRECISION 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 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 ztray, zgar, zratio, zff, zfacoa, zcorae
    DOUBLE PRECISION zmue, zgap, zww, zto, zden, zmu1, zden1
    DOUBLE PRECISION zbmu0, zbmu1, zre11

    ! ------------------------------------------------------------------

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


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

    DO jk = 1, kflev
       DO jl = 1, kdlon
          ptauaz(jl, jk) = 0d0
          ppizaz(jl, jk) = 0d0
          pcgaz(jl, jk) = 0d0
       END DO

       IF (flag_aer) THEN
          ! -OB
          DO jl = 1, kdlon
             ztray = prayl(jl)*pdsig(jl, jk)
             zratio = ztray/(ztray+ptauaz(jl,jk))
             zgar = pcgaz(jl, jk)
             zff = zgar*zgar
             ptauaz(jl, jk) = ztray + ptauaz(jl, jk)*(1.-ppizaz(jl,jk)*zff)
             pcgaz(jl, jk) = zgar*(1.-zratio)/(1.+zgar)
             ppizaz(jl, jk) = zratio + (1.-zratio)*ppizaz(jl, jk)*(1.-zff)/(1.- &
                  ppizaz(jl,jk)*zff)
          END DO
       ELSE
          DO jl = 1, kdlon
             ztray = prayl(jl)*pdsig(jl, jk)
             ptauaz(jl, jk) = ztray
             pcgaz(jl, jk) = 0.
             ppizaz(jl, jk) = 1. - repsct
          END DO
       END IF
    END DO

    ! ------------------------------------------------------------------

    ! *         2.    TOTAL EFFECTIVE CLOUDINESS ABOVE A GIVEN LEVEL
    ! ----------------------------------------------


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

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

       IF (novlp==1) THEN
          ! * maximum-random
          zclear(jl) = zclear(jl)*(1.0-max(zss0(jl),zscat(jl)))/ &
               (1.0-min(zscat(jl),1.-zepsec))
          zc0i(jl, jkl) = 1.0 - 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)*(1.0-zss0(jl))
          zscat(jl) = 1.0 - 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 = 1. - ppizaz(jl, jkl)*pcgaz(jl, jkl)*pcgaz(jl, jkl)
          zcorae = zfacoa*ptauaz(jl, jkl)*psec(jl)
          zr21(jl) = exp(-zcorae)
          zss0(jl) = 1. - zr21(jl)

          IF (novlp==1) THEN
             ! * maximum-random
             zclear(jl) = zclear(jl)*(1.0-max(zss0(jl),zscat(jl)))/ &
                  (1.0-min(zscat(jl),1.-zepsec))
             zc0i(jl, jkl) = 1.0 - 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)*(1.0-zss0(jl))
             zscat(jl) = 1.0 - 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) = 0.
       pray2(jl, kflev+1) = 0.
       prefz(jl, 2, 1) = palbp(jl, knu)
       prefz(jl, 1, 1) = palbp(jl, knu)
       ptra1(jl, kflev+1) = 1.
       ptra2(jl, kflev+1) = 1.
    END DO

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


          ! ------------------------------------------------------------------

          ! *         3.1  EQUIVALENT ZENITH ANGLE
          ! -----------------------


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


          ! ------------------------------------------------------------------

          ! *         3.2  REFLECT./TRANSMISSIVITY DUE TO RAYLEIGH AND AEROSOLS
          ! ----------------------------------------------------


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

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


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

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