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