Line data Source code
1 : ! path: $Source: /storm/rc1/cvsroot/rc/rrtmg_sw/src/rrtmg_sw_reftra.f90,v $
2 : ! author: $Author: mike $
3 : ! revision: $Revision: 1.2 $
4 : ! created: $Date: 2007/08/23 20:40:14 $
5 :
6 : module rrtmg_sw_reftra
7 :
8 : ! --------------------------------------------------------------------------
9 : ! | |
10 : ! | Copyright 2002-2007, Atmospheric & Environmental Research, Inc. (AER). |
11 : ! | This software may be used, copied, or redistributed as long as it is |
12 : ! | not sold and this copyright notice is reproduced on each copy made. |
13 : ! | This model is provided as is without any express or implied warranties. |
14 : ! | (http://www.rtweb.aer.com/) |
15 : ! | |
16 : ! --------------------------------------------------------------------------
17 :
18 : ! ------- Modules -------
19 :
20 : use shr_kind_mod, only: r8 => shr_kind_r8
21 :
22 : use rrsw_tbl, only : tblint, bpade, od_lo, exp_tbl
23 :
24 : implicit none
25 :
26 : contains
27 :
28 : ! --------------------------------------------------------------------
29 8601600 : subroutine reftra_sw(nlayers, ncol, lrtchk, pgg, prmuz, ptau, pw, &
30 8601600 : pref, prefd, ptra, ptrad)
31 : ! --------------------------------------------------------------------
32 :
33 : ! Purpose: computes the reflectivity and transmissivity of a clear or
34 : ! cloudy layer using a choice of various approximations.
35 : !
36 : ! Interface: *rrtmg_sw_reftra* is called by *rrtmg_sw_spcvrt*
37 : !
38 : ! Description:
39 : ! explicit arguments :
40 : ! --------------------
41 : ! inputs
42 : ! ------
43 : ! lrtchk = .t. for all layers in clear profile
44 : ! lrtchk = .t. for cloudy layers in cloud profile
45 : ! = .f. for clear layers in cloud profile
46 : ! pgg = asymmetry factor
47 : ! prmuz = cosine solar zenith angle
48 : ! ptau = optical thickness
49 : ! pw = single scattering albedo
50 : !
51 : ! outputs
52 : ! -------
53 : ! pref : collimated beam reflectivity
54 : ! prefd : diffuse beam reflectivity
55 : ! ptra : collimated beam transmissivity
56 : ! ptrad : diffuse beam transmissivity
57 : !
58 : !
59 : ! Method:
60 : ! -------
61 : ! standard delta-eddington, p.i.f.m., or d.o.m. layer calculations.
62 : ! kmodts = 1 eddington (joseph et al., 1976)
63 : ! = 2 pifm (zdunkowski et al., 1980)
64 : ! = 3 discrete ordinates (liou, 1973)
65 : !
66 : !
67 : ! Modifications:
68 : ! --------------
69 : ! Original: J-JMorcrette, ECMWF, Feb 2003
70 : ! Revised for F90 reformatting: MJIacono, AER, Jul 2006
71 : ! Revised to add exponential lookup table: MJIacono, AER, Aug 2007
72 : !
73 : ! ------------------------------------------------------------------
74 :
75 : ! ------- Declarations ------
76 :
77 : ! ------- Input -------
78 :
79 : integer, intent(in) :: nlayers
80 : integer, intent(in) :: ncol
81 : logical, intent(in) :: lrtchk(ncol,nlayers) ! Logical flag for reflectivity and
82 : ! and transmissivity calculation;
83 : ! Dimensions: (nlayers)
84 :
85 : real(kind=r8), intent(in) :: pgg(ncol,nlayers) ! asymmetry parameter
86 : ! Dimensions: (nlayers)
87 : real(kind=r8), intent(in) :: ptau(ncol,nlayers) ! optical depth
88 : ! Dimensions: (nlayers)
89 : real(kind=r8), intent(in) :: pw(ncol,nlayers) ! single scattering albedo
90 : ! Dimensions: (nlayers)
91 : real(kind=r8), intent(in) :: prmuz(ncol) ! cosine of solar zenith angle
92 :
93 : ! ------- Output -------
94 :
95 : real(kind=r8), intent(inout) :: pref(ncol,nlayers+1) ! direct beam reflectivity
96 : ! Dimensions: (nlayers+1)
97 : real(kind=r8), intent(inout) :: prefd(ncol,nlayers+1) ! diffuse beam reflectivity
98 : ! Dimensions: (nlayers+1)
99 : real(kind=r8), intent(inout) :: ptra(ncol,nlayers+1) ! direct beam transmissivity
100 : ! Dimensions: (nlayers+1)
101 : real(kind=r8), intent(inout) :: ptrad(ncol,nlayers+1) ! diffuse beam transmissivity
102 : ! Dimensions: (nlayers+1)
103 :
104 : ! ------- Local -------
105 :
106 : integer :: jk, icol, kmodts
107 : integer :: itind
108 :
109 : real(kind=r8) :: tblind
110 : real(kind=r8) :: za, za1, za2
111 : real(kind=r8) :: zbeta, zdend, zdenr, zdent
112 : real(kind=r8) :: ze1, ze2, zem1, zem2, zemm, zep1, zep2
113 : real(kind=r8) :: zg, zg3, zgamma1, zgamma2, zgamma3, zgamma4, zgt
114 : real(kind=r8) :: zr1, zr2, zr3, zr4, zr5
115 : real(kind=r8) :: zrk, zrk2, zrkg, zrm1, zrp, zrp1, zrpp
116 : real(kind=r8) :: zsr3, zt1, zt2, zt3, zt4, zt5, zto1
117 : real(kind=r8) :: zw, zwcrit, zwo
118 : real(kind=r8) :: temp1
119 : real(kind=r8), parameter :: eps = 1.e-08_r8
120 :
121 : ! ------------------------------------------------------------------
122 :
123 : ! Initialize
124 :
125 8601600 : zsr3=sqrt(3._r8)
126 8601600 : zwcrit=0.9999995_r8
127 8601600 : kmodts=2
128 :
129 481689600 : do jk=1, nlayers
130 3887923200 : do icol = 1,ncol
131 3879321600 : if (.not.lrtchk(icol,jk)) then
132 1557510521 : pref(icol,jk) =0._r8
133 1557510521 : ptra(icol,jk) =1._r8
134 1557510521 : prefd(icol,jk)=0._r8
135 1557510521 : ptrad(icol,jk)=1._r8
136 : else
137 1848723079 : zto1=ptau(icol,jk)
138 1848723079 : zw =pw(icol,jk)
139 1848723079 : zg =pgg(icol,jk)
140 :
141 : ! General two-stream expressions
142 :
143 1848723079 : zg3= 3._r8 * zg
144 : if (kmodts == 1) then
145 : zgamma1= (7._r8 - zw * (4._r8 + zg3)) * 0.25_r8
146 : zgamma2=-(1._r8 - zw * (4._r8 - zg3)) * 0.25_r8
147 : zgamma3= (2._r8 - zg3 * prmuz(icol) ) * 0.25_r8
148 : else if (kmodts == 2) then
149 1848723079 : zgamma1= (8._r8 - zw * (5._r8 + zg3)) * 0.25_r8
150 1848723079 : zgamma2= 3._r8 *(zw * (1._r8 - zg )) * 0.25_r8
151 1848723079 : zgamma3= (2._r8 - zg3 * prmuz(icol) ) * 0.25_r8
152 : else if (kmodts == 3) then
153 : zgamma1= zsr3 * (2._r8 - zw * (1._r8 + zg)) * 0.5_r8
154 : zgamma2= zsr3 * zw * (1._r8 - zg ) * 0.5_r8
155 : zgamma3= (1._r8 - zsr3 * zg * prmuz(icol) ) * 0.5_r8
156 : end if
157 1848723079 : zgamma4= 1._r8 - zgamma3
158 :
159 : ! Recompute original s.s.a. to test for conservative solution
160 :
161 1848723079 : temp1 = 1._r8 - 2._r8 * zg
162 1848723079 : zwo= zw * (temp1 + zg**2)/(temp1 + zg**2 * zw)
163 1848723079 : if (zwo >= zwcrit) then
164 : ! Conservative scattering
165 :
166 21899621 : za = zgamma1 * prmuz(icol)
167 21899621 : za1 = za - zgamma3
168 21899621 : zgt = zgamma1 * zto1
169 :
170 : ! Homogeneous reflectance and transmittance,
171 : ! collimated beam
172 :
173 21899621 : ze1 = min ( zto1 / prmuz(icol) , 500._r8)
174 21899621 : ze2 = exp( -ze1 )
175 :
176 : ! Use exponential lookup table for transmittance, or expansion of
177 : ! exponential for low tau
178 : ! if (ze1 .le. od_lo) then
179 : ! ze2 = 1._r8 - ze1 + 0.5_r8 * ze1 * ze1
180 : ! else
181 : ! tblind = ze1 / (bpade + ze1)
182 : ! itind = tblint * tblind + 0.5_r8
183 : ! ze2 = exp_tbl(itind)
184 : ! endif
185 : !
186 :
187 21899621 : pref(icol,jk) = (zgt - za1 * (1._r8 - ze2)) / ( 1._r8 + zgt)
188 21899621 : ptra(icol,jk) = 1._r8 - pref(icol,jk)
189 :
190 : ! isotropic incidence
191 :
192 21899621 : prefd(icol,jk) = zgt / ( 1._r8 + zgt)
193 21899621 : ptrad(icol,jk) = 1._r8 - prefd(icol,jk)
194 :
195 : ! This is applied for consistency between total (delta-scaled) and direct (unscaled)
196 : ! calculations at very low optical depths (tau < 1.e-4) when the exponential lookup
197 : ! table returns a transmittance of 1.0.
198 21899621 : if (ze2 .eq. 1.0_r8) then
199 0 : pref(icol,jk) = 0.0_r8
200 0 : ptra(icol,jk) = 1.0_r8
201 0 : prefd(icol,jk) = 0.0_r8
202 0 : ptrad(icol,jk) = 1.0_r8
203 : endif
204 :
205 : else
206 : ! Non-conservative scattering
207 :
208 1826823458 : za1 = zgamma1 * zgamma4 + zgamma2 * zgamma3
209 1826823458 : za2 = zgamma1 * zgamma3 + zgamma2 * zgamma4
210 1826823458 : zrk = sqrt ( zgamma1**2 - zgamma2**2)
211 1826823458 : zrp = zrk * prmuz(icol)
212 1826823458 : zrp1 = 1._r8 + zrp
213 1826823458 : zrm1 = 1._r8 - zrp
214 1826823458 : zrk2 = 2._r8 * zrk
215 1826823458 : zrpp = 1._r8 - zrp*zrp
216 1826823458 : zrkg = zrk + zgamma1
217 1826823458 : zr1 = zrm1 * (za2 + zrk * zgamma3)
218 1826823458 : zr2 = zrp1 * (za2 - zrk * zgamma3)
219 1826823458 : zr3 = zrk2 * (zgamma3 - za2 * prmuz(icol) )
220 1826823458 : zr4 = zrpp * zrkg
221 1826823458 : zr5 = zrpp * (zrk - zgamma1)
222 1826823458 : zt1 = zrp1 * (za1 + zrk * zgamma4)
223 1826823458 : zt2 = zrm1 * (za1 - zrk * zgamma4)
224 1826823458 : zt3 = zrk2 * (zgamma4 + za1 * prmuz(icol) )
225 1826823458 : zt4 = zr4
226 1826823458 : zt5 = zr5
227 1826823458 : zbeta = (zgamma1 - zrk) / zrkg !- zr5 / zr4
228 :
229 : ! Homogeneous reflectance and transmittance
230 :
231 1826823458 : ze1 = min ( zrk * zto1, 500._r8)
232 1826823458 : ze2 = min ( zto1 / prmuz(icol) , 500._r8)
233 : !
234 : ! Original
235 : ! zep1 = exp( ze1 )
236 : ! zem1 = exp(-ze1 )
237 : ! zep2 = exp( ze2 )
238 : ! zem2 = exp(-ze2 )
239 : !
240 : ! Revised original, to reduce exponentials
241 1826823458 : zep1 = exp( ze1 )
242 1826823458 : zem1 = 1._r8 / zep1
243 1826823458 : zep2 = exp( ze2 )
244 1826823458 : zem2 = 1._r8 / zep2
245 : !
246 : ! Use exponential lookup table for transmittance, or expansion of
247 : ! exponential for low tau
248 : ! if (ze1 .le. od_lo) then
249 : ! zem1 = 1._r8 - ze1 + 0.5_r8 * ze1 * ze1
250 : ! zep1 = 1._r8 / zem1
251 : ! else
252 : ! tblind = ze1 / (bpade + ze1)
253 : ! itind = tblint * tblind + 0.5_r8
254 : ! zem1 = exp_tbl(itind)
255 : ! zep1 = 1._r8 / zem1
256 : ! endif
257 : !
258 : ! if (ze2 .le. od_lo) then
259 : ! zem2 = 1._r8 - ze2 + 0.5_r8 * ze2 * ze2
260 : ! zep2 = 1._r8 / zem2
261 : ! else
262 : ! tblind = ze2 / (bpade + ze2)
263 : ! itind = tblint * tblind + 0.5_r8
264 : ! zem2 = exp_tbl(itind)
265 : ! zep2 = 1._r8 / zem2
266 : ! endif
267 :
268 : ! collimated beam
269 :
270 1826823458 : zdenr = zr4*zep1 + zr5*zem1
271 1826823458 : zdent = zt4*zep1 + zt5*zem1
272 1826823458 : if (zdenr .ge. -eps .and. zdenr .le. eps) then
273 1 : pref(icol,jk) = eps
274 1 : ptra(icol,jk) = zem2
275 : else
276 1826823457 : pref(icol,jk) = zw * (zr1*zep1 - zr2*zem1 - zr3*zem2) / zdenr
277 1826823457 : ptra(icol,jk) = zem2 - zem2 * zw * (zt1*zep1 - zt2*zem1 - zt3*zep2) / zdent
278 : endif
279 :
280 : ! diffuse beam
281 :
282 1826823458 : zemm = zem1*zem1
283 1826823458 : zdend = 1._r8 / ( (1._r8 - zbeta*zemm ) * zrkg)
284 1826823458 : prefd(icol,jk) = zgamma2 * (1._r8 - zemm) * zdend
285 1826823458 : ptrad(icol,jk) = zrk2*zem1*zdend
286 :
287 : endif
288 :
289 : endif
290 :
291 : enddo
292 : end do
293 8601600 : end subroutine reftra_sw
294 :
295 : end module rrtmg_sw_reftra
296 :
|
