Line data Source code
1 : ! This code is part of Radiative Transfer for Energetics (RTE)
2 : !
3 : ! Contacts: Robert Pincus and Eli Mlawer
4 : ! email: rrtmgp@aer.com
5 : !
6 : ! Copyright 2015-, Atmospheric and Environmental Research,
7 : ! Regents of the University of Colorado, Trustees of Columbia University. All right reserved.
8 : !
9 : ! Use and duplication is permitted under the terms of the
10 : ! BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause
11 : ! -------------------------------------------------------------------------------------------------
12 : !
13 : !> ## Kernels for arrays of optical properties:
14 : !> - delta-scaling
15 : !> - adding two sets of properties
16 : !> - extracting subsets along the column dimension
17 : !
18 : ! -------------------------------------------------------------------------------------------------
19 :
20 : module mo_optical_props_kernels
21 : use, intrinsic :: iso_c_binding
22 : use mo_rte_kind, only: wp, wl
23 : implicit none
24 :
25 : public
26 :
27 : !> Delta-scale two-stream optical properties
28 : interface delta_scale_2str_kernel
29 : module procedure delta_scale_2str_f_k, delta_scale_2str_k
30 : end interface
31 :
32 : !> Subsetting, meaning extracting some portion of the 3D domain
33 : interface extract_subset
34 : module procedure extract_subset_dim1_3d, extract_subset_dim2_4d
35 : module procedure extract_subset_absorption_tau
36 : end interface extract_subset
37 :
38 : real(wp), parameter, private :: eps = 3.0_wp*tiny(1.0_wp)
39 : contains
40 : ! -------------------------------------------------------------------------------------------------
41 : !
42 : ! Delta-scaling is provided only for two-stream properties at present
43 : !
44 : ! -------------------------------------------------------------------------------------------------
45 : !> Delta-scale two-stream optical properties given user-provided value of \(f\) (forward scattering)
46 : !
47 0 : pure subroutine delta_scale_2str_f_k(ncol, nlay, ngpt, tau, ssa, g, f) &
48 : bind(C, name="rte_delta_scale_2str_f_k")
49 : integer, intent(in ) :: ncol, nlay, ngpt
50 : !! Array sizes
51 : real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: tau, ssa, g
52 : !! Optical depth, single-scattering albedo, asymmetry parameter
53 : real(wp), dimension(ncol, nlay, ngpt), intent(in ) :: f
54 : !! User-provided forward-scattering fraction
55 :
56 : real(wp) :: wf
57 : integer :: icol, ilay, igpt
58 :
59 0 : do igpt = 1, ngpt
60 0 : do ilay = 1, nlay
61 0 : do icol = 1, ncol
62 0 : wf = ssa(icol,ilay,igpt) * f(icol,ilay,igpt)
63 0 : tau(icol,ilay,igpt) = (1._wp - wf) * tau(icol,ilay,igpt)
64 0 : ssa(icol,ilay,igpt) = (ssa(icol,ilay,igpt) - wf) / max(eps,(1.0_wp - wf))
65 : g (icol,ilay,igpt) = (g (icol,ilay,igpt) - f(icol,ilay,igpt)) / &
66 0 : max(eps,(1._wp - f(icol,ilay,igpt)))
67 : end do
68 : end do
69 : end do
70 :
71 0 : end subroutine delta_scale_2str_f_k
72 : ! ---------------------------------
73 : !> Delta-scale assuming forward-scatternig fraction is the square of the asymmetry parameter
74 : !> i.e. \(f = g^2\)
75 : !
76 389263 : pure subroutine delta_scale_2str_k(ncol, nlay, ngpt, tau, ssa, g) &
77 : bind(C, name="rte_delta_scale_2str_k")
78 : integer, intent(in ) :: ncol, nlay, ngpt
79 : !! Array sizes
80 : real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: tau, ssa, g
81 : !! Optical depth, single-scattering albedo, asymmetry parameter
82 :
83 : real(wp) :: f, wf
84 : integer :: icol, ilay, igpt
85 :
86 43986719 : do igpt = 1, ngpt
87 4142147583 : do ilay = 1, nlay
88 65796884720 : do icol = 1, ncol
89 61655126400 : f = g (icol,ilay,igpt) * g (icol,ilay,igpt)
90 61655126400 : wf = ssa(icol,ilay,igpt) * f
91 61655126400 : tau(icol,ilay,igpt) = (1._wp - wf) * tau(icol,ilay,igpt)
92 61655126400 : ssa(icol,ilay,igpt) = (ssa(icol,ilay,igpt) - wf) / max(eps,(1.0_wp - wf))
93 65753287264 : g (icol,ilay,igpt) = (g (icol,ilay,igpt) - f) / max(eps,(1.0_wp - f))
94 : end do
95 : end do
96 : end do
97 :
98 389263 : end subroutine delta_scale_2str_k
99 : ! -------------------------------------------------------------------------------------------------
100 : !
101 : ! Addition of optical properties: the first set are incremented by the second set.
102 : !
103 : ! There are three possible representations of optical properties (scalar = optical depth only;
104 : ! two-stream = tau, single-scattering albedo, and asymmetry factor g, and
105 : ! n-stream = tau, ssa, and phase function moments p.) Thus we need nine routines, three for
106 : ! each choice of representation on the left hand side times three representations of the
107 : ! optical properties to be added.
108 : !
109 : ! There are two sets of these nine routines. In the first the two sets of optical
110 : ! properties are defined at the same spectral resolution. There is also a set of routines
111 : ! to add properties defined at lower spectral resolution to a set defined at higher spectral
112 : ! resolution (adding properties defined by band to those defined by g-point)
113 : !
114 : ! -------------------------------------------------------------------------------------------------
115 : !> increase one absorption optical depth by a second value
116 746136 : pure subroutine increment_1scalar_by_1scalar(ncol, nlay, ngpt, &
117 746136 : tau1, &
118 746136 : tau2) bind(C, name="rte_increment_1scalar_by_1scalar")
119 : integer, intent(in ) :: ncol, nlay, ngpt !! array sizes
120 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified
121 : real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2 !! optical properties to be added to original
122 :
123 : integer :: icol, ilay, igpt
124 :
125 96251544 : do igpt = 1, ngpt
126 9073759896 : do ilay = 1, nlay
127 >14999*10^7 : do icol = 1, ncol
128 >14990*10^7 : tau1(icol,ilay,igpt) = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
129 : end do
130 : end do
131 : end do
132 746136 : end subroutine increment_1scalar_by_1scalar
133 : ! ---------------------------------
134 : !> increase absorption optical depth with extinction optical depth (2-stream form)
135 0 : pure subroutine increment_1scalar_by_2stream(ncol, nlay, ngpt, &
136 0 : tau1, &
137 0 : tau2, ssa2) bind(C, name="rte_increment_1scalar_by_2stream")
138 : integer, intent(in ) :: ncol, nlay, ngpt !! array sizes
139 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified
140 : real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2 !! optical properties to be added to original
141 :
142 : integer :: icol, ilay, igpt
143 :
144 0 : do igpt = 1, ngpt
145 0 : do ilay = 1, nlay
146 0 : do icol = 1, ncol
147 0 : tau1(icol,ilay,igpt) = tau1(icol,ilay,igpt) + &
148 0 : tau2(icol,ilay,igpt) * (1._wp - ssa2(icol,ilay,igpt))
149 : end do
150 : end do
151 : end do
152 0 : end subroutine increment_1scalar_by_2stream
153 : ! ---------------------------------
154 : !> increase absorption optical depth with extinction optical depth (n-stream form)
155 0 : pure subroutine increment_1scalar_by_nstream(ncol, nlay, ngpt, &
156 0 : tau1, &
157 0 : tau2, ssa2) bind(C, name="rte_increment_1scalar_by_nstream")
158 : integer, intent(in ) :: ncol, nlay, ngpt !! array sizes
159 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified
160 : real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2 !! optical properties to be added to original
161 :
162 : integer :: icol, ilay, igpt
163 :
164 0 : do igpt = 1, ngpt
165 0 : do ilay = 1, nlay
166 0 : do icol = 1, ncol
167 0 : tau1(icol,ilay,igpt) = tau1(icol,ilay,igpt) + &
168 0 : tau2(icol,ilay,igpt) * (1._wp - ssa2(icol,ilay,igpt))
169 : end do
170 : end do
171 : end do
172 0 : end subroutine increment_1scalar_by_nstream
173 : ! ---------------------------------
174 : ! ---------------------------------
175 : !> increment two-stream optical properties \(\tau, \omega_0, g\) with absorption optical depth
176 0 : pure subroutine increment_2stream_by_1scalar(ncol, nlay, ngpt, &
177 0 : tau1, ssa1, &
178 0 : tau2) bind(C, name="rte_increment_2stream_by_1scalar")
179 : integer, intent(in ) :: ncol, nlay, ngpt !! array sizes
180 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified
181 : real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2 !! optical properties to be added to original
182 :
183 : integer :: icol, ilay, igpt
184 : real(wp) :: tau12
185 :
186 0 : do igpt = 1, ngpt
187 0 : do ilay = 1, nlay
188 0 : do icol = 1, ncol
189 0 : tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
190 0 : ssa1(icol,ilay,igpt) = tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) / max(eps,tau12)
191 0 : tau1(icol,ilay,igpt) = tau12
192 : ! g is unchanged
193 : end do
194 : end do
195 : end do
196 0 : end subroutine increment_2stream_by_1scalar
197 : ! ---------------------------------
198 : !> increment two-stream optical properties \(\tau, \omega_0, g\) with a second set
199 389263 : pure subroutine increment_2stream_by_2stream(ncol, nlay, ngpt, &
200 389263 : tau1, ssa1, g1, &
201 389263 : tau2, ssa2, g2) bind(C, name="rte_increment_2stream_by_2stream")
202 : integer, intent(in ) :: ncol, nlay, ngpt !! array sizes
203 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified
204 : real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2, g2 !! optical properties to be added to original
205 :
206 : integer :: icol, ilay, igpt
207 : real(wp) :: tau12, tauscat12
208 :
209 43986719 : do igpt = 1, ngpt
210 4142147583 : do ilay = 1, nlay
211 65796884720 : do icol = 1, ncol
212 : ! t=tau1 + tau2
213 61655126400 : tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
214 : ! w=(tau1*ssa1 + tau2*ssa2) / t
215 : tauscat12 = tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
216 61655126400 : tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt)
217 : g1(icol,ilay,igpt) = &
218 : (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * g1(icol,ilay,igpt) + &
219 : tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt) * g2(icol,ilay,igpt)) &
220 61655126400 : / max(eps,tauscat12)
221 61655126400 : ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
222 65753287264 : tau1(icol,ilay,igpt) = tau12
223 : end do
224 : end do
225 : end do
226 389263 : end subroutine increment_2stream_by_2stream
227 : ! ---------------------------------
228 : !> increment two-stream optical properties \(\tau, \omega_0, g\) with _n_-stream
229 0 : pure subroutine increment_2stream_by_nstream(ncol, nlay, ngpt, nmom2, &
230 0 : tau1, ssa1, g1, &
231 0 : tau2, ssa2, p2) bind(C, name="rte_increment_2stream_by_nstream")
232 : integer, intent(in ) :: ncol, nlay, ngpt, nmom2 !! array sizes
233 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified
234 : real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2 !! optical properties to be added to original
235 : real(wp), dimension(nmom2, &
236 : ncol,nlay,ngpt), intent(in ) :: p2 !! moments of the phase function to be added
237 :
238 : integer :: icol, ilay, igpt
239 : real(wp) :: tau12, tauscat12
240 :
241 0 : do igpt = 1, ngpt
242 0 : do ilay = 1, nlay
243 0 : do icol = 1, ncol
244 : ! t=tau1 + tau2
245 0 : tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
246 : ! w=(tau1*ssa1 + tau2*ssa2) / t
247 : tauscat12 = &
248 : tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
249 0 : tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt)
250 : g1(icol,ilay,igpt) = &
251 : (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * g1( icol,ilay,igpt)+ &
252 0 : tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt) * p2(1, icol,ilay,igpt)) / max(eps,tauscat12)
253 0 : ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
254 0 : tau1(icol,ilay,igpt) = tau12
255 : end do
256 : end do
257 : end do
258 0 : end subroutine increment_2stream_by_nstream
259 : ! ---------------------------------
260 : ! ---------------------------------
261 : !> increment _n_-stream optical properties \(\tau, \omega_0, p\) with absorption optical depth
262 0 : pure subroutine increment_nstream_by_1scalar(ncol, nlay, ngpt, &
263 0 : tau1, ssa1, &
264 0 : tau2) bind(C, name="rte_increment_nstream_by_1scalar")
265 : integer, intent(in ) :: ncol, nlay, ngpt !! array sizes
266 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified
267 : real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2 !! optical properties to be added to original
268 :
269 : integer :: icol, ilay, igpt
270 : real(wp) :: tau12
271 :
272 0 : do igpt = 1, ngpt
273 0 : do ilay = 1, nlay
274 0 : do icol = 1, ncol
275 0 : tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
276 0 : ssa1(icol,ilay,igpt) = tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) / max(eps,tau12)
277 0 : tau1(icol,ilay,igpt) = tau12
278 : ! p is unchanged
279 : end do
280 : end do
281 : end do
282 0 : end subroutine increment_nstream_by_1scalar
283 : ! ---------------------------------
284 : !> increment _n_-stream optical properties \(\tau, \omega_0, p\) with two-stream values
285 0 : pure subroutine increment_nstream_by_2stream(ncol, nlay, ngpt, nmom1, &
286 0 : tau1, ssa1, p1, &
287 0 : tau2, ssa2, g2) bind(C, name="rte_increment_nstream_by_2stream")
288 : integer, intent(in ) :: ncol, nlay, ngpt, nmom1 !! array sizes
289 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified
290 : real(wp), dimension(nmom1, &
291 : ncol,nlay,ngpt), intent(inout) :: p1 !! moments of the phase function be modified
292 : real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2, g2 !! optical properties to be added to original
293 :
294 : integer :: icol, ilay, igpt
295 : real(wp) :: tau12, tauscat12
296 0 : real(wp), dimension(nmom1) :: temp_moms ! TK
297 : integer :: imom !TK
298 :
299 0 : do igpt = 1, ngpt
300 0 : do ilay = 1, nlay
301 0 : do icol = 1, ncol
302 0 : tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
303 : tauscat12 = &
304 : tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
305 0 : tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt)
306 : !
307 : ! Here assume Henyey-Greenstein
308 : !
309 0 : temp_moms(1) = g2(icol,ilay,igpt)
310 0 : do imom = 2, nmom1
311 0 : temp_moms(imom) = temp_moms(imom-1) * g2(icol,ilay,igpt)
312 : end do
313 : p1(1:nmom1, icol,ilay,igpt) = &
314 : (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * p1(1:nmom1, icol,ilay,igpt) + &
315 0 : tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt) * temp_moms(1:nmom1) ) / max(eps,tauscat12)
316 0 : ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
317 0 : tau1(icol,ilay,igpt) = tau12
318 : end do
319 : end do
320 : end do
321 0 : end subroutine increment_nstream_by_2stream
322 : ! ---------------------------------
323 : !> increment one set of _n_-stream optical properties with another set
324 0 : pure subroutine increment_nstream_by_nstream(ncol, nlay, ngpt, nmom1, nmom2, &
325 0 : tau1, ssa1, p1, &
326 0 : tau2, ssa2, p2) bind(C, name="rte_increment_nstream_by_nstream")
327 : integer, intent(in ) :: ncol, nlay, ngpt, nmom1, nmom2 !! array sizes
328 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified
329 : real(wp), dimension(nmom1, &
330 : ncol,nlay,ngpt), intent(inout) :: p1 !! moments of the phase function be modified
331 : real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2 !! optical properties to be added to original
332 : real(wp), dimension(nmom2, &
333 : ncol,nlay,ngpt), intent(in ) :: p2 !! moments of the phase function to be added
334 :
335 : integer :: icol, ilay, igpt, mom_lim
336 : real(wp) :: tau12, tauscat12
337 :
338 0 : mom_lim = min(nmom1, nmom2)
339 0 : do igpt = 1, ngpt
340 0 : do ilay = 1, nlay
341 0 : do icol = 1, ncol
342 0 : tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
343 : tauscat12 = &
344 : tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
345 0 : tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt)
346 : !
347 : ! If op2 has more moments than op1 these are ignored;
348 : ! if it has fewer moments the higher orders are assumed to be 0
349 : !
350 : p1(1:mom_lim, icol,ilay,igpt) = &
351 0 : (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * p1(1:mom_lim, icol,ilay,igpt) + &
352 0 : tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt) * p2(1:mom_lim, icol,ilay,igpt)) / max(eps,tauscat12)
353 0 : ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
354 0 : tau1(icol,ilay,igpt) = tau12
355 : end do
356 : end do
357 : end do
358 0 : end subroutine increment_nstream_by_nstream
359 : ! -------------------------------------------------------------------------------------------------
360 : !
361 : ! Incrementing when the second set of optical properties is defined at lower spectral resolution
362 : ! (e.g. by band instead of by gpoint)
363 : !
364 : ! -------------------------------------------------------------------------------------------------
365 : !> increase one absorption optical depth defined on g-points by a second value defined on bands
366 746136 : pure subroutine inc_1scalar_by_1scalar_bybnd(ncol, nlay, ngpt, &
367 746136 : tau1, &
368 746136 : tau2, &
369 746136 : nbnd, gpt_lims) bind(C, name="rte_inc_1scalar_by_1scalar_bybnd")
370 : integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes
371 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified (defined on g-points)
372 : real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2 !! optical properties to be added to original (defined on bands)
373 : integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band
374 :
375 : integer :: ibnd, igpt
376 :
377 12684312 : do ibnd = 1, nbnd
378 108189720 : do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
379 >15001*10^7 : tau1(:,:,igpt) = tau1(:,:,igpt) + tau2(:,:,ibnd)
380 : end do
381 : end do
382 746136 : end subroutine inc_1scalar_by_1scalar_bybnd
383 : ! ---------------------------------
384 : !> increase absorption optical depth defined on g-points with extinction optical depth (2-stream form) defined on bands
385 0 : pure subroutine inc_1scalar_by_2stream_bybnd(ncol, nlay, ngpt, &
386 0 : tau1, &
387 0 : tau2, ssa2, &
388 0 : nbnd, gpt_lims) bind(C, name="rte_inc_1scalar_by_2stream_bybnd")
389 : integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes
390 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified (defined on g-points)
391 : real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands)
392 : integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band
393 :
394 : integer :: ibnd, igpt
395 :
396 0 : do ibnd = 1, nbnd
397 0 : do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
398 0 : tau1(:,:,igpt) = tau1(:,:,igpt) + tau2(:,:,ibnd) * (1._wp - ssa2(:,:,ibnd))
399 : end do
400 : end do
401 0 : end subroutine inc_1scalar_by_2stream_bybnd
402 : ! ---------------------------------
403 : !> increase absorption optical depth defined on g-points with extinction optical depth (n-stream form) defined on bands
404 0 : pure subroutine inc_1scalar_by_nstream_bybnd(ncol, nlay, ngpt, &
405 0 : tau1, &
406 0 : tau2, ssa2, &
407 0 : nbnd, gpt_lims) bind(C, name="rte_inc_1scalar_by_nstream_bybnd")
408 : integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes
409 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified (defined on g-points)
410 : real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands)
411 : integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band
412 :
413 : integer :: ibnd, igpt
414 :
415 0 : do ibnd = 1, nbnd
416 0 : do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
417 0 : tau1(:,:,igpt) = tau1(:,:,igpt) + tau2(:,:,ibnd) * (1._wp - ssa2(:,:,ibnd))
418 : end do
419 : end do
420 0 : end subroutine inc_1scalar_by_nstream_bybnd
421 :
422 : ! ---------------------------------
423 : !> increment two-stream optical properties \(\tau, \omega_0, g\) defined on g-points with absorption optical depth defined on bands
424 0 : pure subroutine inc_2stream_by_1scalar_bybnd(ncol, nlay, ngpt, &
425 0 : tau1, ssa1, &
426 0 : tau2, &
427 0 : nbnd, gpt_lims) bind(C, name="rte_inc_2stream_by_1scalar_bybnd")
428 : integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes
429 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points)
430 : real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2 !! optical properties to be added to original (defined on bands)
431 : integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band
432 :
433 : integer :: icol, ilay, igpt, ibnd
434 : real(wp) :: tau12
435 :
436 0 : do ibnd = 1, nbnd
437 0 : do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
438 0 : do ilay = 1, nlay
439 0 : do icol = 1, ncol
440 0 : tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,ibnd)
441 0 : ssa1(icol,ilay,igpt) = tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) / max(eps,tau12)
442 0 : tau1(icol,ilay,igpt) = tau12
443 : ! g is unchanged
444 : end do
445 : end do
446 : end do
447 : end do
448 0 : end subroutine inc_2stream_by_1scalar_bybnd
449 : ! ---------------------------------
450 : !> increment 2-stream optical properties defined on g-points with another set defined on bands
451 389263 : pure subroutine inc_2stream_by_2stream_bybnd(ncol, nlay, ngpt, &
452 389263 : tau1, ssa1, g1, &
453 389263 : tau2, ssa2, g2, &
454 389263 : nbnd, gpt_lims) bind(C, name="rte_inc_2stream_by_2stream_bybnd")
455 : integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes
456 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified (defined on g-points)
457 : real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2, g2 !! optical properties to be added to original (defined on bands)
458 : integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band
459 :
460 : integer :: icol, ilay, igpt, ibnd
461 : real(wp) :: tau12, tauscat12
462 :
463 5838945 : do ibnd = 1, nbnd
464 49436401 : do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
465 4147208002 : do ilay = 1, nlay
466 65796884720 : do icol = 1, ncol
467 : ! t=tau1 + tau2
468 61655126400 : tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,ibnd)
469 : ! w=(tau1*ssa1 + tau2*ssa2) / t
470 : tauscat12 = &
471 : tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
472 61655126400 : tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd)
473 : g1(icol,ilay,igpt) = &
474 : (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * g1(icol,ilay,igpt) + &
475 61655126400 : tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd) * g2(icol,ilay,ibnd)) / max(eps,tauscat12)
476 61655126400 : ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
477 65753287264 : tau1(icol,ilay,igpt) = tau12
478 : end do
479 : end do
480 : end do
481 : end do
482 389263 : end subroutine inc_2stream_by_2stream_bybnd
483 : ! ---------------------------------
484 : !> increment 2-stream optical properties defined on g-points with _n_-stream properties set defined on bands
485 0 : pure subroutine inc_2stream_by_nstream_bybnd(ncol, nlay, ngpt, nmom2, &
486 0 : tau1, ssa1, g1, &
487 0 : tau2, ssa2, p2, &
488 0 : nbnd, gpt_lims) bind(C, name="rte_inc_2stream_by_nstream_bybnd")
489 : integer, intent(in ) :: ncol, nlay, ngpt, nmom2, nbnd !! array sizes
490 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified (defined on g-points)
491 : real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands)
492 : real(wp), dimension(nmom2, &
493 : ncol,nlay,nbnd), intent(in ) :: p2 !! moments of the phase function to be added
494 : integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band
495 :
496 : integer :: icol, ilay, igpt, ibnd
497 : real(wp) :: tau12, tauscat12
498 :
499 0 : do ibnd = 1, nbnd
500 0 : do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
501 0 : do ilay = 1, nlay
502 0 : do icol = 1, ncol
503 : ! t=tau1 + tau2
504 0 : tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,ibnd)
505 : ! w=(tau1*ssa1 + tau2*ssa2) / t
506 : tauscat12 = &
507 : tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
508 0 : tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd)
509 : g1(icol,ilay,igpt) = &
510 : (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * g1( icol,ilay,igpt)+ &
511 0 : tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd) * p2(1, icol,ilay,ibnd)) / max(eps,tauscat12)
512 0 : ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
513 0 : tau1(icol,ilay,igpt) = tau12
514 : end do
515 : end do
516 : end do
517 : end do
518 0 : end subroutine inc_2stream_by_nstream_bybnd
519 : ! ---------------------------------
520 : ! ---------------------------------
521 : !> increment _n_-stream optical properties defined on g-points with absorption optical depth defined on bands
522 0 : pure subroutine inc_nstream_by_1scalar_bybnd(ncol, nlay, ngpt, &
523 0 : tau1, ssa1, &
524 0 : tau2, &
525 0 : nbnd, gpt_lims) bind(C, name="rte_inc_nstream_by_1scalar_bybnd")
526 : integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes
527 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points)
528 : real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2 !! optical properties to be added to original (defined on bands)
529 : integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band
530 :
531 : integer :: icol, ilay, igpt, ibnd
532 : real(wp) :: tau12
533 :
534 0 : do ibnd = 1, nbnd
535 0 : do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
536 0 : do ilay = 1, nlay
537 0 : do icol = 1, ncol
538 0 : tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,ibnd)
539 0 : ssa1(icol,ilay,igpt) = tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) / max(eps,tau12)
540 0 : tau1(icol,ilay,igpt) = tau12
541 : ! p is unchanged
542 : end do
543 : end do
544 : end do
545 : end do
546 0 : end subroutine inc_nstream_by_1scalar_bybnd
547 : ! ---------------------------------
548 : !> increment n-stream optical properties defined on g-points with 2-stream properties set defined on bands
549 0 : pure subroutine inc_nstream_by_2stream_bybnd(ncol, nlay, ngpt, nmom1, &
550 0 : tau1, ssa1, p1, &
551 0 : tau2, ssa2, g2, &
552 0 : nbnd, gpt_lims) bind(C, name="rte_inc_nstream_by_2stream_bybnd")
553 : integer, intent(in ) :: ncol, nlay, ngpt, nmom1, nbnd !! array sizes
554 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points)
555 : real(wp), dimension(nmom1, &
556 : ncol,nlay,ngpt), intent(inout) :: p1 !! moments of the phase function be modified
557 : real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2, g2 !! optical properties to be added to original (defined on bands)
558 : integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band
559 :
560 : integer :: icol, ilay, igpt, ibnd
561 : real(wp) :: tau12, tauscat12
562 0 : real(wp), dimension(nmom1) :: temp_moms ! TK
563 : integer :: imom !TK
564 :
565 0 : do ibnd = 1, nbnd
566 0 : do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
567 0 : do ilay = 1, nlay
568 0 : do icol = 1, ncol
569 0 : tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,ibnd)
570 : tauscat12 = &
571 : tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
572 0 : tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd)
573 : !
574 : ! Here assume Henyey-Greenstein
575 : !
576 0 : temp_moms(1) = g2(icol,ilay,ibnd)
577 0 : do imom = 2, nmom1
578 0 : temp_moms(imom) = temp_moms(imom-1) * g2(icol,ilay,ibnd)
579 : end do
580 : p1(1:nmom1, icol,ilay,igpt) = &
581 : (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * p1(1:nmom1, icol,ilay,igpt) + &
582 0 : tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd) * temp_moms(1:nmom1) ) / max(eps,tauscat12)
583 0 : ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
584 0 : tau1(icol,ilay,igpt) = tau12
585 : end do
586 : end do
587 : end do
588 : end do
589 0 : end subroutine inc_nstream_by_2stream_bybnd
590 : ! ---------------------------------
591 : !> increment _n_-stream optical properties defined on g-points with a second set defined on bands
592 0 : pure subroutine inc_nstream_by_nstream_bybnd(ncol, nlay, ngpt, nmom1, nmom2, &
593 0 : tau1, ssa1, p1, &
594 0 : tau2, ssa2, p2, &
595 0 : nbnd, gpt_lims) bind(C, name="rte_inc_nstream_by_nstream_bybnd")
596 : integer, intent(in ) :: ncol, nlay, ngpt, nmom1, nmom2, nbnd !! array sizes
597 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points)
598 : real(wp), dimension(nmom1, &
599 : ncol,nlay,ngpt), intent(inout) :: p1 !! moments of the phase function be modified
600 : real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands)
601 : real(wp), dimension(nmom2, &
602 : ncol,nlay,nbnd), intent(in ) :: p2 !! moments of the phase function to be added
603 : integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band
604 :
605 : integer :: icol, ilay, igpt, ibnd, mom_lim
606 : real(wp) :: tau12, tauscat12
607 :
608 0 : mom_lim = min(nmom1, nmom2)
609 0 : do ibnd = 1, nbnd
610 0 : do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
611 0 : do ilay = 1, nlay
612 0 : do icol = 1, ncol
613 0 : tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,ibnd)
614 : tauscat12 = &
615 : tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
616 0 : tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd)
617 : !
618 : ! If op2 has more moments than op1 these are ignored;
619 : ! if it has fewer moments the higher orders are assumed to be 0
620 : !
621 : p1(1:mom_lim, icol,ilay,igpt) = &
622 0 : (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * p1(1:mom_lim, icol,ilay,igpt) + &
623 0 : tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd) * p2(1:mom_lim, icol,ilay,ibnd)) / max(eps,tauscat12)
624 0 : ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
625 0 : tau1(icol,ilay,igpt) = tau12
626 : end do
627 : end do
628 : end do
629 : end do
630 0 : end subroutine inc_nstream_by_nstream_bybnd
631 : ! -------------------------------------------------------------------------------------------------
632 : !
633 : ! Subsetting, meaning extracting some portion of the 3D domain
634 : !
635 : ! -------------------------------------------------------------------------------------------------
636 : !>
637 : !> Extract a subset from the first dimension (normally columns) of a 3D field.
638 : !> Applicable to most variables e.g. tau, ssa, g
639 : !>
640 0 : pure subroutine extract_subset_dim1_3d(ncol, nlay, ngpt, array_in, colS, colE, array_out) &
641 : bind (C, name="rte_extract_subset_dim1_3d")
642 : integer, intent(in ) :: ncol, nlay, ngpt !! Array sizes
643 : real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: array_in !! Array to subset
644 : integer, intent(in ) :: colS, colE !! Starting and ending index
645 : real(wp), dimension(colE-colS+1,&
646 : nlay,ngpt), intent(out) :: array_out !! subset of the input array
647 :
648 : integer :: icol, ilay, igpt
649 0 : do igpt = 1, ngpt
650 0 : do ilay = 1, nlay
651 0 : do icol = colS, colE
652 0 : array_out(icol-colS+1, ilay, igpt) = array_in(icol, ilay, igpt)
653 : end do
654 : end do
655 : end do
656 :
657 0 : end subroutine extract_subset_dim1_3d
658 : ! ---------------------------------
659 : !> Extract a subset from the second dimension (normally columns) of a 4D field.
660 : !> Applicable to phase function moments, where the first dimension is the moment
661 0 : pure subroutine extract_subset_dim2_4d(nmom, ncol, nlay, ngpt, array_in, colS, colE, array_out) &
662 : bind (C, name="rte_extract_subset_dim2_4d")
663 : integer, intent(in ) :: nmom, ncol, nlay, ngpt !! Array sizes
664 : real(wp), dimension(nmom,ncol,nlay,ngpt), intent(in ) :: array_in !! Array to subset
665 : integer, intent(in ) :: colS, colE !! Starting and ending index
666 : real(wp), dimension(nmom,colE-colS+1,&
667 : nlay,ngpt), intent(out) :: array_out !! subset of the input array
668 :
669 : integer :: icol, ilay, igpt, imom
670 :
671 0 : do igpt = 1, ngpt
672 0 : do ilay = 1, nlay
673 0 : do icol = colS, colE
674 0 : do imom = 1, nmom
675 0 : array_out(imom, icol-colS+1, ilay, igpt) = array_in(imom, icol, ilay, igpt)
676 : end do
677 : end do
678 : end do
679 : end do
680 :
681 0 : end subroutine extract_subset_dim2_4d
682 : ! ---------------------------------
683 : !
684 : !> Extract the absorption optical thickness \(\tau_{abs} = 1 - \omega_0 \tau_{ext}\)
685 : !
686 0 : pure subroutine extract_subset_absorption_tau(ncol, nlay, ngpt, tau_in, ssa_in, &
687 0 : colS, colE, tau_out) &
688 : bind (C, name="rte_extract_subset_absorption_tau")
689 : integer, intent(in ) :: ncol, nlay, ngpt !! Array sizes
690 : real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau_in, ssa_in !! Optical thickness, single scattering albedo
691 : integer, intent(in ) :: colS, colE !! Starting and ending index
692 : real(wp), dimension(colE-colS+1,&
693 : nlay,ngpt), intent(out) :: tau_out !! absorption optical thickness subset
694 :
695 : integer :: icol, ilay, igpt
696 :
697 0 : do igpt = 1, ngpt
698 0 : do ilay = 1, nlay
699 0 : do icol = colS, colE
700 0 : tau_out(icol-colS+1, ilay, igpt) = &
701 0 : tau_in(icol, ilay, igpt) * (1._wp - ssa_in(icol, ilay, igpt))
702 : end do
703 : end do
704 : end do
705 :
706 0 : end subroutine extract_subset_absorption_tau
707 : end module mo_optical_props_kernels
|
