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 : !>## Numeric calculations for radiative transfer solvers
14 : !> - Emission/absorption (no-scattering) calculations
15 : !> - solver for multi-angle Gaussian quadrature
16 : !> - solver for a single angle, calling
17 : !> - source function computation (linear-in-tau)
18 : !> - transport
19 : !> - Extinction-only calculation (direct solar beam)
20 : !> - Two-stream calculations:
21 : !> solvers for LW and SW with different boundary conditions and source functions
22 : !> - source function calculation for LW, SW
23 : !> - two-stream calculations for LW, SW (using different assumtions about phase function)
24 : !> - transport (adding)
25 : !> - Application of boundary conditions
26 : !
27 : ! -------------------------------------------------------------------------------------------------
28 : module mo_rte_solver_kernels
29 : use, intrinsic :: iso_c_binding
30 : use mo_rte_kind, only: wp, wl
31 : use mo_rte_util_array,only: zero_array
32 : implicit none
33 : private
34 :
35 : public :: lw_solver_noscat, lw_solver_2stream, &
36 : sw_solver_noscat, sw_solver_2stream
37 :
38 : real(wp), parameter :: pi = acos(-1._wp)
39 : contains
40 : ! -------------------------------------------------------------------------------------------------
41 : !
42 : ! Top-level longwave kernels
43 : !
44 : ! -------------------------------------------------------------------------------------------------
45 : !
46 : !> LW fluxes, no scattering, mu (cosine of integration angle) specified by column
47 : !> Does radiation calculation at user-supplied angles; converts radiances to flux
48 : !> using user-supplied weights
49 : !
50 : ! ---------------------------------------------------------------
51 1492272 : subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, &
52 1492272 : tau, lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, &
53 1492272 : incident_flux, &
54 1492272 : flux_up, flux_dn, &
55 1492272 : do_broadband, broadband_up, broadband_dn, &
56 1492272 : do_Jacobians, sfc_srcJac, flux_upJac, &
57 1492272 : do_rescaling, ssa, g)
58 : integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points
59 : logical(wl), intent(in ) :: top_at_1
60 : real(wp), dimension(ncol, ngpt), intent(in ) :: D ! secant of propagation angle []
61 : real(wp), intent(in ) :: weight ! quadrature weight
62 : real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau ! Absorption optical thickness []
63 : real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lay_source ! Planck source at layer average temperature [W/m2]
64 : ! Planck source at layer edge for radiation in increasing/decreasing ilay direction
65 : ! lev_source_dec applies the mapping in layer i to the Planck function at layer i
66 : ! lev_source_inc applies the mapping in layer i to the Planck function at layer i+1
67 : real(wp), dimension(ncol,nlay, ngpt), target, &
68 : intent(in ) :: lev_source_inc, lev_source_dec
69 : real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis ! Surface emissivity []
70 : real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src ! Surface source function [W/m2]
71 : real(wp), dimension(ncol, ngpt), intent(in ) :: incident_flux! Boundary condition for flux [W/m2]
72 : real(wp), dimension(ncol,nlay+1,ngpt), target, & ! Fluxes [W/m2]
73 : intent( out) :: flux_up, flux_dn
74 : !
75 : ! Optional variables - arrays aren't referenced if corresponding logical == False
76 : !
77 : logical(wl), intent(in ) :: do_broadband
78 : real(wp), dimension(ncol,nlay+1 ), intent( out) :: broadband_up, broadband_dn ! Spectrally-integrated fluxes [W/m2]
79 : logical(wl), intent(in ) :: do_Jacobians
80 : real(wp), dimension(ncol ,ngpt), intent(in ) :: sfc_srcJac ! surface temperature Jacobian of surface source function [W/m2/K]
81 : real(wp), dimension(ncol,nlay+1 ), intent( out) :: flux_upJac ! surface temperature Jacobian of Radiances [W/m2-str / K]
82 : logical(wl), intent(in ) :: do_rescaling
83 : real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: ssa, g ! single-scattering albedo, asymmetry parameter
84 : ! ------------------------------------
85 : ! Local variables, no g-point dependency
86 : !
87 : integer :: icol, ilay, igpt
88 : integer :: top_level, sfc_level
89 2984544 : real(wp), dimension(ncol,nlay) :: tau_loc, & ! path length (tau/mu)
90 2984544 : trans ! transmissivity = exp(-tau)
91 2984544 : real(wp), dimension(ncol,nlay) :: source_dn, source_up
92 2984544 : real(wp), dimension(ncol ) :: sfc_albedo
93 :
94 1492272 : real(wp), dimension(:,:,:), pointer :: lev_source_up, lev_source_dn ! Mapping increasing/decreasing indicies to up/down
95 :
96 : real(wp), parameter :: pi = acos(-1._wp)
97 : ! loc_fluxes hold a single g-point flux if fluxes are being integrated instead of returned
98 : ! with spectral detail
99 : real(wp), dimension(ncol,nlay+1), &
100 2984544 : target :: loc_flux_up, loc_flux_dn
101 : ! gpt_fluxes point to calculations for the current g-point
102 1492272 : real(wp), dimension(:,:), pointer :: gpt_flux_up, gpt_flux_dn
103 : ! -------------------------------------------------------------------------------------------------
104 : ! Optionally, use an approximate treatment of scattering using rescaling
105 : ! Implemented based on the paper
106 : ! Tang G, et al, 2018: https://doi.org/10.1175/JAS-D-18-0014.1
107 : ! a) relies on rescaling of the optical parameters based on asymetry factor and single scattering albedo
108 : ! scaling can be computed by scaling_1rescl
109 : ! b) adds adustment term based on cloud properties (lw_transport_1rescl)
110 : ! adustment terms is computed based on solution of the Tang equations
111 : ! for "linear-in-tau" internal source (not in the paper)
112 : !
113 : ! Used when approximating scattering
114 : !
115 : real(wp) :: ssal, wb, scaleTau
116 2984544 : real(wp), dimension(ncol,nlay ) :: An, Cn
117 2984544 : real(wp), dimension(ncol,nlay+1) :: gpt_flux_Jac
118 : ! ------------------------------------
119 : ! Which way is up?
120 : ! Level Planck sources for upward and downward radiation
121 : ! When top_at_1, lev_source_up => lev_source_dec
122 : ! lev_source_dn => lev_source_inc, and vice-versa
123 1492272 : if(top_at_1) then
124 1492272 : top_level = 1
125 1492272 : sfc_level = nlay+1
126 1492272 : lev_source_up => lev_source_dec
127 1492272 : lev_source_dn => lev_source_inc
128 : else
129 0 : top_level = nlay+1
130 0 : sfc_level = 1
131 0 : lev_source_up => lev_source_inc
132 0 : lev_source_dn => lev_source_dec
133 : end if
134 :
135 : !
136 : ! Integrated fluxes need zeroing
137 : !
138 1492272 : if(do_broadband) then
139 746136 : call zero_array(ncol, nlay+1, broadband_up )
140 746136 : call zero_array(ncol, nlay+1, broadband_dn )
141 : end if
142 1492272 : if(do_Jacobians) &
143 0 : call zero_array(ncol, nlay+1, flux_upJac )
144 :
145 192503088 : do igpt = 1, ngpt
146 191010816 : if(do_broadband) then
147 95505408 : gpt_flux_up => loc_flux_up
148 95505408 : gpt_flux_dn => loc_flux_dn
149 : else
150 95505408 : gpt_flux_up => flux_up (:,:,igpt)
151 95505408 : gpt_flux_dn => flux_dn (:,:,igpt)
152 : end if
153 : !
154 : ! Transport is for intensity
155 : ! convert flux at top of domain to intensity assuming azimuthal isotropy
156 : !
157 3189436416 : gpt_flux_dn(:,top_level) = incident_flux(:,igpt)/(2._wp * pi * weight)
158 : !
159 : ! Optical path and transmission, used in source function and transport calculations
160 : !
161 191010816 : if (do_rescaling) then
162 : !
163 : ! The scaling and scaleTau terms are independent of propagation
164 : ! angle D and could be pre-computed if several values of D are used
165 : ! We re-compute them here to keep not have to localize memory use
166 : !
167 0 : do ilay = 1, nlay
168 0 : do icol = 1, ncol
169 0 : ssal = ssa(icol, ilay, igpt)
170 :
171 : ! w is the layer single scattering albedo
172 : ! b is phase function parameter (Eq.13 of the paper)
173 : ! for the similarity principle scaling scheme
174 : ! b = (1-g)/2 (where g is phase function avergae cosine)
175 0 : wb = ssal*(1._wp - g(icol, ilay, igpt)) * 0.5_wp
176 :
177 : ! scaleTau=1-w(1-b) is a scaling factor of the optical thickness representing
178 : ! the radiative transfer equation in a nonscattering form Eq(14) of the paper
179 0 : scaleTau = (1._wp - ssal + wb)
180 :
181 : ! Cn = 0.5*wb/(1-w(1-b)) is parameter of Eq.21-22 of the Tang paper
182 : ! Tang paper, p.2222 advises to replace 0.5 with 0.4 based on simulations
183 0 : Cn(icol,ilay) = 0.4_wp*wb/scaleTau
184 :
185 : ! Eqs.15, 18ab and 19 of the paper,
186 : ! rescaling of the optical depth multiplied by path length
187 0 : tau_loc(icol,ilay) = tau(icol,ilay,igpt)*D(icol,igpt)*scaleTau
188 : end do
189 0 : trans (:,ilay) = exp(-tau_loc(:,ilay))
190 0 : An(:,ilay) = (1._wp-trans(:,ilay)**2)
191 : end do
192 : else
193 18146027520 : do ilay = 1, nlay
194 >29980*10^7 : tau_loc(:,ilay) = tau(:,ilay,igpt)*D(:,igpt)
195 >29999*10^7 : trans (:,ilay) = exp(-tau_loc(:,ilay))
196 : end do
197 : end if
198 : !
199 : ! Source function for diffuse radiation
200 : !
201 : call lw_source_noscat(ncol, nlay, &
202 : lay_source(:,:,igpt), lev_source_up(:,:,igpt), lev_source_dn(:,:,igpt), &
203 191010816 : tau_loc, trans, source_dn, source_up)
204 : !
205 : ! Transport down
206 : !
207 191010816 : call lw_transport_noscat_dn(ncol, nlay, top_at_1, trans, source_dn, gpt_flux_dn)
208 : !
209 : ! Surface albedo, surface source function, reflection and emission
210 : !
211 3189436416 : sfc_albedo(:) = 1._wp - sfc_emis(:,igpt)
212 0 : gpt_flux_up (:,sfc_level) = gpt_flux_dn(:,sfc_level)*sfc_albedo(:) + &
213 6187862016 : sfc_emis(:,igpt) * sfc_src (:,igpt)
214 191010816 : if(do_Jacobians) &
215 0 : gpt_flux_Jac(:,sfc_level) = sfc_emis(:,igpt) * sfc_srcJac(:,igpt)
216 : !
217 : ! Transport up, or up and down again if using rescaling
218 : !
219 191010816 : if(do_rescaling) then
220 : call lw_transport_1rescl(ncol, nlay, top_at_1, trans, &
221 : source_dn, source_up, &
222 : gpt_flux_up, gpt_flux_dn, An, Cn, &
223 0 : do_Jacobians, gpt_flux_Jac) ! Standing in for Jacobian, i.e. rad_up_Jac(:,:,igpt), rad_dn_Jac(:,:,igpt))
224 : else
225 : call lw_transport_noscat_up(ncol, nlay, top_at_1, trans, source_up, gpt_flux_up, &
226 191010816 : do_Jacobians, gpt_flux_Jac)
227 : end if
228 :
229 191010816 : if(do_broadband) then
230 >15159*10^7 : broadband_up(:,:) = broadband_up(:,:) + gpt_flux_up(:,:)
231 >15159*10^7 : broadband_dn(:,:) = broadband_dn(:,:) + gpt_flux_dn(:,:)
232 : else
233 : !
234 : ! Convert intensity to flux assuming azimuthal isotropy and quadrature weight
235 : !
236 >15159*10^7 : gpt_flux_dn(:,:) = 2._wp * pi * weight * gpt_flux_dn(:,:)
237 >15159*10^7 : gpt_flux_up(:,:) = 2._wp * pi * weight * gpt_flux_up(:,:)
238 : end if
239 : !
240 : ! Only broadband-integrated Jacobians are provided
241 : !
242 191010816 : if(do_Jacobians) &
243 1492272 : flux_upJac(:,:) = flux_upJac(:,:) + gpt_flux_Jac(:,:)
244 : end do ! g point loop
245 :
246 1492272 : if(do_broadband) then
247 1184326056 : broadband_up(:,:) = 2._wp * pi * weight* broadband_up(:,:)
248 1184326056 : broadband_dn(:,:) = 2._wp * pi * weight* broadband_dn(:,:)
249 : end if
250 1492272 : if(do_Jacobians) &
251 0 : flux_upJac(:,:) = 2._wp * pi * weight * flux_upJac(:,:)
252 :
253 1492272 : end subroutine lw_solver_noscat_oneangle
254 : ! -------------------------------------------------------------------------------------------------
255 : !
256 : !> LW transport, no scattering, multi-angle quadrature
257 : !> Users provide a set of weights and quadrature angles
258 : !> Routine sums over single-angle solutions for each sets of angles/weights
259 : !
260 : ! ---------------------------------------------------------------
261 1492272 : subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, &
262 1492272 : nmus, Ds, weights, &
263 1492272 : tau, &
264 1492272 : lay_source, lev_source_inc, lev_source_dec, &
265 1492272 : sfc_emis, sfc_src, &
266 1492272 : inc_flux, &
267 1492272 : flux_up, flux_dn, &
268 1492272 : do_broadband, broadband_up, broadband_dn, &
269 1492272 : do_Jacobians, sfc_srcJac, flux_upJac, &
270 1492272 : do_rescaling, ssa, g) bind(C, name="rte_lw_solver_noscat")
271 : integer, intent(in ) :: ncol, nlay, ngpt
272 : !! Number of columns, layers, g-points
273 : logical(wl), intent(in ) :: top_at_1
274 : !! ilay = 1 is the top of the atmosphere?
275 : integer, intent(in ) :: nmus
276 : !! number of quadrature angles
277 : real(wp), dimension (ncol, ngpt, &
278 : nmus), intent(in ) :: Ds
279 : !! quadrature secants
280 : real(wp), dimension(nmus), intent(in ) :: weights
281 : !! quadrature weights
282 : real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau
283 : !! Absorption optical thickness []
284 : real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lay_source
285 : !! Planck source at layer average temperature [W/m2]
286 : real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lev_source_inc
287 : !! Planck source at layer edge for radiation in increasing ilay direction [W/m2]
288 : real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lev_source_dec
289 : !! Planck source at layer edge for radiation in decreasing ilay direction [W/m2]
290 : real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis
291 : !! Surface emissivity []
292 : real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src
293 : !! Surface source function [W/m2]
294 : real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux
295 : !! Incident diffuse flux, probably 0 [W/m2]
296 : real(wp), dimension(ncol,nlay+1,ngpt), target, &
297 : intent( out) :: flux_up, flux_dn
298 : !! Fluxes [W/m2]
299 : !
300 : ! Optional variables - arrays aren't referenced if corresponding logical == False
301 : !
302 : logical(wl), intent(in ) :: do_broadband
303 : real(wp), dimension(ncol,nlay+1 ), target, &
304 : intent( out) :: broadband_up, broadband_dn
305 : !! Spectrally-integrated fluxes [W/m2]
306 : logical(wl), intent(in ) :: do_Jacobians
307 : !! compute Jacobian with respect to surface temeprature?
308 : real(wp), dimension(ncol ,ngpt), intent(in ) :: sfc_srcJac
309 : !! surface temperature Jacobian of surface source function [W/m2/K]
310 : real(wp), dimension(ncol,nlay+1 ), target, &
311 : intent( out) :: flux_upJac
312 : !! surface temperature Jacobian of Radiances [W/m2-str / K]
313 : logical(wl), intent(in ) :: do_rescaling
314 : !! Approximate treatment of scattering (10.1175/JAS-D-18-0014.1)
315 : real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: ssa, g
316 : !! single-scattering albedo, asymmetry parameter
317 : ! ------------------------------------
318 : !
319 : ! Local variables - used for a single quadrature angle
320 : !
321 1492272 : real(wp), dimension(:,:,:), pointer :: this_flux_up, this_flux_dn
322 1492272 : real(wp), dimension(:,:), pointer :: this_broadband_up, this_broadband_dn, this_flux_upJac
323 :
324 : integer :: imu
325 : ! ------------------------------------
326 : !
327 : ! For the first angle output arrays store total flux
328 : !
329 : call lw_solver_noscat_oneangle(ncol, nlay, ngpt, &
330 0 : top_at_1, Ds(:,:,1), weights(1), tau, &
331 : lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, &
332 : inc_flux, &
333 : flux_up, flux_dn, &
334 : do_broadband, broadband_up, broadband_dn, &
335 : do_Jacobians, sfc_srcJac, flux_upJac, &
336 1492272 : do_rescaling, ssa, g)
337 : !
338 : ! For more than one angle use local arrays
339 : !
340 1492272 : if(nmus > 1) then
341 0 : if(do_broadband) then
342 0 : allocate(this_broadband_up(ncol,nlay+1), this_broadband_dn(ncol,nlay+1))
343 : ! Spectrally-resolved fluxes won't be filled in so can point to caller-supplied memory
344 0 : this_flux_up => flux_up
345 0 : this_flux_dn => flux_dn
346 : else
347 0 : allocate(this_flux_up(ncol,nlay+1,ngpt), this_flux_dn(ncol,nlay+1,ngpt))
348 : ! Spectrally-integrated fluxes won't be filled in so can point to caller-supplied memory
349 0 : this_broadband_up => broadband_up
350 0 : this_broadband_dn => broadband_dn
351 : end if
352 0 : if(do_Jacobians) then
353 0 : allocate(this_flux_upJac(ncol,nlay+1))
354 : else
355 0 : this_flux_upJac => flux_upJac
356 : end if
357 : end if
358 1492272 : do imu = 2, nmus
359 : call lw_solver_noscat_oneangle(ncol, nlay, ngpt, &
360 0 : top_at_1, Ds(:,:,imu), weights(imu), tau, &
361 : lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, &
362 : inc_flux, &
363 : this_flux_up, this_flux_dn, &
364 : do_broadband, this_broadband_up, this_broadband_dn, &
365 : do_Jacobians, sfc_srcJac, this_flux_upJac, &
366 0 : do_rescaling, ssa, g)
367 0 : if(do_broadband) then
368 0 : broadband_up(:,:) = broadband_up(:,:) + this_broadband_up(:,:)
369 0 : broadband_dn(:,:) = broadband_dn(:,:) + this_broadband_dn(:,:)
370 : else
371 0 : flux_up (:,:,:) = flux_up (:,:,:) + this_flux_up (:,:,:)
372 0 : flux_dn (:,:,:) = flux_dn (:,:,:) + this_flux_dn (:,:,:)
373 : end if
374 0 : if (do_Jacobians) &
375 1492272 : flux_upJac(:,:) = flux_upJac(:,: ) + this_flux_upJac(:,: )
376 : end do
377 1492272 : if(nmus > 1) then
378 0 : if( do_broadband) deallocate(this_broadband_up, this_broadband_dn)
379 0 : if(.not. do_broadband) deallocate(this_flux_up, this_flux_dn)
380 0 : if( do_Jacobians) deallocate(this_flux_upJac)
381 : end if
382 1492272 : end subroutine lw_solver_noscat
383 : ! -------------------------------------------------------------------------------------------------
384 : !
385 : !> Longwave two-stream calculation:
386 : !> - combine RRTMGP-specific sources at levels
387 : !> - compute layer reflectance, transmittance
388 : !> - compute total source function at levels using linear-in-tau
389 : !> - transport
390 : !
391 : ! -------------------------------------------------------------------------------------------------
392 0 : subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, &
393 0 : tau, ssa, g, &
394 0 : lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, &
395 0 : inc_flux, &
396 0 : flux_up, flux_dn) bind(C, name="rte_lw_solver_2stream")
397 : integer, intent(in ) :: ncol, nlay, ngpt
398 : !! Number of columns, layers, g-points
399 : logical(wl), intent(in ) :: top_at_1
400 : !! ilay = 1 is the top of the atmosphere?
401 : real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau, ssa, g
402 : !! Optical thickness, single-scattering albedo, asymmetry parameter []
403 : real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lay_source
404 : !! Planck source at layer average temperature [W/m2]
405 : real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lev_source_inc
406 : !! Planck source at layer edge for radiation in increasing ilay direction [W/m2]
407 : real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lev_source_dec
408 : !! Planck source at layer edge for radiation in decreasing ilay direction [W/m2]
409 : real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis
410 : !! Surface emissivity []
411 : real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src
412 : !! Surface source function [W/m2]
413 : real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux
414 : !! Incident diffuse flux, probably 0 [W/m2]
415 : real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: flux_up, flux_dn
416 : !! Fluxes [W/m2]
417 : ! ----------------------------------------------------------------------
418 : integer :: igpt, top_level
419 0 : real(wp), dimension(ncol,nlay ) :: Rdif, Tdif, gamma1, gamma2
420 0 : real(wp), dimension(ncol ) :: sfc_albedo
421 0 : real(wp), dimension(ncol,nlay+1) :: lev_source
422 0 : real(wp), dimension(ncol,nlay ) :: source_dn, source_up
423 0 : real(wp), dimension(ncol ) :: source_sfc
424 : ! ------------------------------------
425 0 : top_level = nlay+1
426 0 : if(top_at_1) top_level = 1
427 0 : do igpt = 1, ngpt
428 : !
429 : ! RRTMGP provides source functions at each level using the spectral mapping
430 : ! of each adjacent layer. Combine these for two-stream calculations
431 : !
432 : call lw_combine_sources(ncol, nlay, top_at_1, &
433 : lev_source_inc(:,:,igpt), lev_source_dec(:,:,igpt), &
434 0 : lev_source)
435 : !
436 : ! Cell properties: reflection, transmission for diffuse radiation
437 : ! Coupling coefficients needed for source function
438 : !
439 : call lw_two_stream(ncol, nlay, &
440 : tau (:,:,igpt), ssa(:,:,igpt), g(:,:,igpt), &
441 0 : gamma1, gamma2, Rdif, Tdif)
442 : !
443 : ! Source function for diffuse radiation
444 : !
445 : call lw_source_2str(ncol, nlay, top_at_1, &
446 : sfc_emis(:,igpt), sfc_src(:,igpt), &
447 : lay_source(:,:,igpt), lev_source, &
448 : gamma1, gamma2, Rdif, Tdif, tau(:,:,igpt), &
449 0 : source_dn, source_up, source_sfc)
450 : !
451 : ! Transport
452 : !
453 0 : sfc_albedo(1:ncol) = 1._wp - sfc_emis(:,igpt)
454 : !
455 : ! Boundary condition
456 : !
457 0 : flux_dn(:,top_level,igpt) = inc_flux(:,igpt)
458 : call adding(ncol, nlay, top_at_1, &
459 : sfc_albedo, &
460 : Rdif, Tdif, &
461 : source_dn, source_up, source_sfc, &
462 0 : flux_up(:,:,igpt), flux_dn(:,:,igpt))
463 : end do
464 :
465 0 : end subroutine lw_solver_2stream
466 : ! -------------------------------------------------------------------------------------------------
467 : !
468 : ! Top-level shortwave kernels
469 : !
470 : ! -------------------------------------------------------------------------------------------------
471 : !
472 : ! !> Extinction-only shortwave solver i.e. solar direct beam
473 : !
474 : ! -------------------------------------------------------------------------------------------------
475 0 : pure subroutine sw_solver_noscat(ncol, nlay, ngpt, top_at_1, &
476 0 : tau, mu0, inc_flux_dir, flux_dir) bind(C, name="rte_sw_solver_noscat")
477 : integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points
478 : !! Number of columns, layers, g-points
479 : logical(wl), intent(in ) :: top_at_1
480 : !! ilay = 1 is the top of the atmosphere?
481 : real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau
482 : !! Absorption optical thickness []
483 : real(wp), dimension(ncol,nlay ), intent(in ) :: mu0
484 : !! cosine of solar zenith angle
485 : real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux_dir
486 : !! Direct beam incident flux [W/m2]
487 : real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: flux_dir
488 : !! Direct-beam flux, spectral [W/m2]
489 :
490 : integer :: ilev, igpt
491 :
492 : ! ------------------------------------
493 : ! Indexing into arrays for upward and downward propagation depends on the vertical
494 : ! orientation of the arrays (whether the domain top is at the first or last index)
495 : ! We write the loops out explicitly so compilers will have no trouble optimizing them.
496 :
497 : ! Downward propagation
498 0 : if(top_at_1) then
499 : ! For the flux at this level, what was the previous level, and which layer has the
500 : ! radiation just passed through?
501 : ! layer index = level index - 1
502 : ! previous level is up (-1)
503 0 : do igpt = 1, ngpt
504 0 : flux_dir(:, 1,igpt) = inc_flux_dir(:,igpt) * mu0(:,1)
505 0 : do ilev = 2, nlay+1
506 0 : flux_dir(:,ilev,igpt) = flux_dir(:,ilev-1,igpt) * exp(-tau(:,ilev-1,igpt)/mu0(:,ilev-1))
507 : end do
508 : end do
509 : else
510 : ! layer index = level index
511 : ! previous level is up (+1)
512 0 : do igpt = 1, ngpt
513 0 : flux_dir(:,nlay+1,igpt) = inc_flux_dir(:,igpt) * mu0(:,nlay)
514 0 : do ilev = nlay, 1, -1
515 0 : flux_dir(:,ilev,igpt) = flux_dir(:,ilev+1,igpt) * exp(-tau(:,ilev,igpt)/mu0(:,ilev))
516 : end do
517 : end do
518 : end if
519 0 : end subroutine sw_solver_noscat
520 : ! -------------------------------------------------------------------------------------------------
521 : !
522 : !> Shortwave two-stream calculation:
523 : !> compute layer reflectance, transmittance
524 : !> compute solar source function for diffuse radiation
525 : !> transport
526 : !
527 : ! -------------------------------------------------------------------------------------------------
528 778526 : subroutine sw_solver_2stream (ncol, nlay, ngpt, top_at_1, &
529 778526 : tau, ssa, g, mu0, &
530 778526 : sfc_alb_dir, sfc_alb_dif, &
531 778526 : inc_flux_dir, &
532 778526 : flux_up, flux_dn, flux_dir, &
533 778526 : has_dif_bc, inc_flux_dif, &
534 778526 : do_broadband, broadband_up, &
535 778526 : broadband_dn, broadband_dir) bind(C, name="rte_sw_solver_2stream")
536 : integer, intent(in ) :: ncol, nlay, ngpt
537 : !! Number of columns, layers, g-points
538 : logical(wl), intent(in ) :: top_at_1
539 : !! ilay = 1 is the top of the atmosphere?
540 : real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau, ssa, g
541 : !! Optical thickness, single-scattering albedo, asymmetry parameter []
542 : real(wp), dimension(ncol,nlay ), intent(in ) :: mu0
543 : !! cosine of solar zenith angle
544 : real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_alb_dir, sfc_alb_dif
545 : !! Spectral surface albedo for direct and diffuse radiation
546 : real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux_dir
547 : !! Direct beam incident flux
548 : real(wp), dimension(ncol,nlay+1,ngpt), target, &
549 : intent( out) :: flux_up, flux_dn, flux_dir
550 : !! Fluxes [W/m2]
551 : logical(wl), intent(in ) :: has_dif_bc
552 : !! Is a boundary condition for diffuse flux supplied?
553 : real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux_dif
554 : !! Boundary condition for diffuse flux [W/m2]
555 : logical(wl), intent(in ) :: do_broadband
556 : !! Provide broadband-integrated, not spectrally-resolved, fluxes?
557 : real(wp), dimension(ncol,nlay+1 ), intent( out) :: broadband_up, broadband_dn, broadband_dir
558 : !! Broadband integrated fluxes
559 : ! -------------------------------------------
560 : integer :: igpt, top_level, top_layer
561 1557052 : real(wp), dimension(ncol,nlay ) :: Rdif, Tdif
562 1557052 : real(wp), dimension(ncol,nlay ) :: source_up, source_dn
563 1557052 : real(wp), dimension(ncol ) :: source_srf
564 : ! loc_fluxes hold a single g-point flux if fluxes are being integrated instead of returned
565 : ! with spectral detail
566 : real(wp), dimension(ncol,nlay+1), &
567 1557052 : target :: loc_flux_up, loc_flux_dn, loc_flux_dir
568 : ! gpt_fluxes point to calculations for the current g-point
569 778526 : real(wp), dimension(:,:), pointer :: gpt_flux_up, gpt_flux_dn, gpt_flux_dir
570 : ! ------------------------------------
571 778526 : if(top_at_1) then
572 : top_level = 1
573 : top_layer = 1
574 : else
575 0 : top_level = nlay+1
576 0 : top_layer = nlay
577 : end if
578 : !
579 : ! Integrated fluxes need zeroing
580 : !
581 778526 : if(do_broadband) then
582 389263 : call zero_array(ncol, nlay+1, broadband_up )
583 389263 : call zero_array(ncol, nlay+1, broadband_dn )
584 389263 : call zero_array(ncol, nlay+1, broadband_dir)
585 : end if
586 :
587 87973438 : do igpt = 1, ngpt
588 87194912 : if(do_broadband) then
589 43597456 : gpt_flux_up => loc_flux_up
590 43597456 : gpt_flux_dn => loc_flux_dn
591 43597456 : gpt_flux_dir => loc_flux_dir
592 : else
593 43597456 : gpt_flux_up => flux_up (:,:,igpt)
594 43597456 : gpt_flux_dn => flux_dn (:,:,igpt)
595 43597456 : gpt_flux_dir => flux_dir(:,:,igpt)
596 : end if
597 : !
598 : ! Boundary conditions direct beam...
599 : !
600 1399006112 : gpt_flux_dir(:,top_level) = inc_flux_dir(:,igpt) * mu0(:,top_layer)
601 : !
602 : ! ... and diffuse field, using 0 if no BC is provided
603 : !
604 87194912 : if(has_dif_bc) then
605 0 : gpt_flux_dn(:,top_level) = inc_flux_dif(:,igpt)
606 : else
607 1399006112 : gpt_flux_dn(:,top_level) = 0._wp
608 : end if
609 : !
610 : ! Cell properties: transmittance and reflectance for diffuse radiation
611 : ! Direct-beam and source for diffuse radiation
612 : !
613 : call sw_dif_and_source(ncol, nlay, top_at_1, mu0, sfc_alb_dir(:,igpt), &
614 : tau(:,:,igpt), ssa(:,:,igpt), g(:,:,igpt), &
615 : Rdif, Tdif, source_dn, source_up, source_srf, &
616 87194912 : gpt_flux_dir)
617 : !
618 : ! Transport
619 : !
620 : call adding(ncol, nlay, top_at_1, &
621 : sfc_alb_dif(:,igpt), Rdif, Tdif, &
622 87194912 : source_dn, source_up, source_srf, gpt_flux_up, gpt_flux_dn)
623 : !
624 : ! adding() computes only diffuse flux; flux_dn is total
625 : !
626 87973438 : if(do_broadband) then
627 66496387776 : broadband_up (:,:) = broadband_up (:,:) + gpt_flux_up (:,:)
628 66496387776 : broadband_dn (:,:) = broadband_dn (:,:) + gpt_flux_dn (:,:) + gpt_flux_dir(:,:)
629 66496387776 : broadband_dir(:,:) = broadband_dir(:,:) + gpt_flux_dir(:,:)
630 : else
631 >13299*10^7 : gpt_flux_dn(:,:) = gpt_flux_dn (:,:) + gpt_flux_dir(:,:)
632 : end if
633 : end do
634 778526 : end subroutine sw_solver_2stream
635 : ! -------------------------------------------------------------------------------------------------
636 : !
637 : ! Lower-level longwave kernels
638 : !
639 : ! -------------------------------------------------------------------------------------------------
640 : !
641 : ! Compute LW source function for upward and downward emission at levels using linear-in-tau assumption
642 : ! See Clough et al., 1992, doi: 10.1029/92JD01419, Eq 13
643 : !
644 : ! ---------------------------------------------------------------
645 191010816 : subroutine lw_source_noscat(ncol, nlay, lay_source, lev_source_up, lev_source_dn, tau, trans, &
646 191010816 : source_dn, source_up)
647 : integer, intent(in) :: ncol, nlay
648 : real(wp), dimension(ncol, nlay), intent(in) :: lay_source, & ! Planck source at layer center
649 : lev_source_up, & ! Planck source at levels (layer edges),
650 : lev_source_dn, & ! increasing/decreasing layer index
651 : tau, & ! Optical path (tau/mu)
652 : trans ! Transmissivity (exp(-tau))
653 : real(wp), dimension(ncol, nlay), intent(out):: source_dn, source_up
654 : ! Source function at layer edges
655 : ! Down at the bottom of the layer, up at the top
656 : ! --------------------------------
657 : integer :: icol, ilay
658 : real(wp) :: fact
659 : real(wp), parameter :: tau_thresh = sqrt(sqrt(epsilon(tau)))
660 : ! ---------------------------------------------------------------
661 18146027520 : do ilay = 1, nlay
662 >29999*10^7 : do icol = 1, ncol
663 : !
664 : ! Weighting factor. Use 2nd order series expansion when rounding error (~tau^2)
665 : ! is of order epsilon (smallest difference from 1. in working precision)
666 : ! Thanks to Peter Blossey
667 : ! Updated to 3rd order series and lower threshold based on suggestion from Dmitry Alexeev (Nvidia)
668 : !
669 >28185*10^7 : if(tau(icol, ilay) > tau_thresh) then
670 >21727*10^7 : fact = (1._wp - trans(icol,ilay))/tau(icol,ilay) - trans(icol,ilay)
671 : else
672 64576670882 : fact = tau(icol, ilay) * (0.5_wp + tau(icol, ilay) * (- 1._wp/3._wp + tau(icol, ilay) * 1._wp/8._wp ) )
673 : end if
674 : !
675 : ! Equation below is developed in Clough et al., 1992, doi:10.1029/92JD01419, Eq 13
676 : !
677 : source_dn(icol,ilay) = (1._wp - trans(icol,ilay)) * lev_source_dn(icol,ilay) + &
678 >28185*10^7 : 2._wp * fact * (lay_source(icol,ilay) - lev_source_dn(icol,ilay))
679 : source_up(icol,ilay) = (1._wp - trans(icol,ilay)) * lev_source_up(icol,ilay ) + &
680 >29980*10^7 : 2._wp * fact * (lay_source(icol,ilay) - lev_source_up(icol,ilay))
681 : end do
682 : end do
683 191010816 : end subroutine lw_source_noscat
684 : ! -------------------------------------------------------------------------------------------------
685 : !
686 : ! Longwave no-scattering transport - separate routines for up and down
687 : !
688 : ! -------------------------------------------------------------------------------------------------
689 191010816 : subroutine lw_transport_noscat_dn(ncol, nlay, top_at_1, &
690 191010816 : trans, source_dn, radn_dn)
691 : integer, intent(in ) :: ncol, nlay ! Number of columns, layers, g-points
692 : logical(wl), intent(in ) :: top_at_1 !
693 : real(wp), dimension(ncol,nlay ), intent(in ) :: trans ! transmissivity = exp(-tau)
694 : real(wp), dimension(ncol,nlay ), intent(in ) :: source_dn ! Diffuse radiation emitted by the layer
695 : real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_dn ! Radiances [W/m2-str] Top level must contain incident flux boundary condition
696 :
697 : ! ---------------------------------------------------
698 : ! Local variables
699 : integer :: ilev
700 : ! ---------------------------------------------------
701 191010816 : if(top_at_1) then
702 : !
703 : ! Top of domain is index 1
704 : !
705 18146027520 : do ilev = 2, nlay+1
706 >29999*10^7 : radn_dn(:,ilev) = trans(:,ilev-1)*radn_dn(:,ilev-1) + source_dn(:,ilev-1)
707 : end do
708 : else
709 : !
710 : ! Top of domain is index nlay+1
711 : !
712 0 : do ilev = nlay, 1, -1
713 0 : radn_dn(:,ilev) = trans(:,ilev )*radn_dn(:,ilev+1) + source_dn(:,ilev)
714 : end do
715 : end if
716 191010816 : end subroutine lw_transport_noscat_dn
717 : ! -------------------------------------------------------------------------------------------------
718 191010816 : subroutine lw_transport_noscat_up(ncol, nlay, top_at_1, &
719 191010816 : trans, source_up, radn_up, do_Jacobians, radn_upJac)
720 : integer, intent(in ) :: ncol, nlay ! Number of columns, layers, g-points
721 : logical(wl), intent(in ) :: top_at_1 !
722 : real(wp), dimension(ncol,nlay ), intent(in ) :: trans ! transmissivity = exp(-tau)
723 : real(wp), dimension(ncol,nlay ), intent(in ) :: source_up ! Diffuse radiation emitted by the layer
724 : real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_up ! Radiances [W/m2-str] Top level must contain incident flux boundary condition
725 : logical(wl), intent(in ) :: do_Jacobians
726 : real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_upJac ! surface temperature Jacobian of Radiances [W/m2-str / K]
727 :
728 : ! ---------------------------------------------------
729 : ! Local variables
730 : integer :: ilev
731 : ! ---------------------------------------------------
732 191010816 : if(top_at_1) then
733 : !
734 : ! Top of domain is index 1
735 : !
736 : ! Upward propagation
737 18146027520 : do ilev = nlay, 1, -1
738 >31776*10^7 : radn_up (:,ilev) = trans(:,ilev )*radn_up (:,ilev+1) + source_up(:,ilev)
739 17955016704 : if(do_Jacobians) &
740 191010816 : radn_upJac(:,ilev) = trans(:,ilev )*radn_upJac(:,ilev+1)
741 : end do
742 : else
743 : !
744 : ! Top of domain is index nlay+1
745 : !
746 : ! Upward propagation
747 0 : do ilev = 2, nlay+1
748 0 : radn_up (:,ilev) = trans(:,ilev-1) * radn_up (:,ilev-1) + source_up(:,ilev-1)
749 0 : if(do_Jacobians) &
750 0 : radn_upJac(:,ilev) = trans(:,ilev-1) * radn_upJac(:,ilev-1)
751 : end do
752 : end if
753 191010816 : end subroutine lw_transport_noscat_up
754 : ! -------------------------------------------------------------------------------------------------
755 : ! Upward and (second) downward transport for re-scaled longwave solution
756 : ! adds adjustment factor based on cloud properties
757 : !
758 : ! implementation notice:
759 : ! the adjustmentFactor computation can be skipped where Cn <= epsilon
760 : ! -------------------------------------------------------------------------------------------------
761 0 : subroutine lw_transport_1rescl(ncol, nlay, top_at_1, &
762 0 : trans, source_dn, source_up, &
763 0 : radn_up, radn_dn, An, Cn,&
764 0 : do_Jacobians, radn_up_Jac)
765 : integer, intent(in ) :: ncol, nlay ! Number of columns, layers, g-points
766 : logical(wl), intent(in ) :: top_at_1 !
767 : real(wp), dimension(ncol,nlay ), intent(in ) :: trans ! transmissivity = exp(-tau)
768 : real(wp), dimension(ncol,nlay ), intent(in ) :: source_dn, &
769 : source_up ! Diffuse radiation emitted by the layer
770 : real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_up ! Radiances [W/m2-str]
771 : real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_dn !Top level must contain incident flux boundary condition
772 : real(wp), dimension(ncol,nlay), intent(in ) :: An, Cn
773 : logical(wl), intent(in ) :: do_Jacobians
774 : real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_up_Jac ! Surface temperature Jacobians [W/m2-str/K]
775 : !
776 : ! We could in principle compute a downwelling Jacobian too, but it's small
777 : ! (only a small proportion of LW is scattered) and it complicates code and the API,
778 : ! so we will not
779 : !
780 :
781 : ! Local variables
782 : integer :: ilev, icol
783 : ! ---------------------------------------------------
784 : real(wp) :: adjustmentFactor
785 0 : if(top_at_1) then
786 : !
787 : ! Top of domain is index 1
788 : !
789 : ! Upward propagation
790 : ! adjustment factor is obtained as a solution of 18b of the Tang paper
791 : ! eqvivalent to Eq.20 of the Tang paper but for linear-in-tau source
792 0 : do ilev = nlay, 1, -1
793 0 : do icol=1,ncol
794 0 : adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_dn(icol,ilev) - &
795 0 : trans(icol,ilev)*source_dn(icol,ilev) - source_up(icol,ilev) )
796 0 : radn_up (icol,ilev) = trans(icol,ilev)*radn_up(icol,ilev+1) + source_up(icol,ilev) + &
797 0 : adjustmentFactor
798 : end do
799 0 : if(do_Jacobians) &
800 0 : radn_up_Jac(:,ilev) = trans(:,ilev)*radn_up_Jac(:,ilev+1)
801 : end do
802 : ! Downward propagation
803 : ! radn_dn_Jac(:,1) = 0._wp
804 : ! adjustment factor is obtained as a solution of 19 of the Tang paper
805 : ! eqvivalent to Eq.21 of the Tang paper but for linear-in-tau source
806 0 : do ilev = 1, nlay
807 : ! radn_dn_Jac(:,ilev+1) = trans(:,ilev)*radn_dn_Jac(:,ilev)
808 0 : do icol=1,ncol
809 0 : adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_up(icol,ilev) - &
810 0 : trans(icol,ilev)*source_up(icol,ilev) - source_dn(icol,ilev) )
811 0 : radn_dn(icol,ilev+1) = trans(icol,ilev)*radn_dn(icol,ilev) + source_dn(icol,ilev) + &
812 0 : adjustmentFactor
813 : ! adjustmentFactor = Cn(icol,ilev)*An(icol,ilev)*radn_up_Jac(icol,ilev)
814 : ! radn_dn_Jac(icol,ilev+1) = radn_dn_Jac(icol,ilev+1) + adjustmentFactor
815 : enddo
816 : end do
817 : else
818 : !
819 : ! Top of domain is index nlay+1
820 : !
821 : ! Upward propagation
822 : ! adjustment factor is obtained as a solution of 18b of the Tang paper
823 : ! eqvivalent to Eq.20 of the Tang paper but for linear-in-tau source
824 0 : do ilev = 1, nlay
825 0 : radn_up (:,ilev+1) = trans(:,ilev) * radn_up (:,ilev) + source_up(:,ilev)
826 0 : do icol=1,ncol
827 0 : adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_dn(icol,ilev+1) - &
828 0 : trans(icol,ilev)*source_dn(icol,ilev) - source_up(icol,ilev) )
829 : radn_up(icol,ilev+1) = trans(icol,ilev)*radn_up(icol,ilev) + source_up(icol,ilev) + &
830 0 : adjustmentFactor
831 : enddo
832 0 : if(do_Jacobians) &
833 0 : radn_up_Jac(:,ilev+1) = trans(:,ilev) * radn_up_Jac(:,ilev)
834 : end do
835 :
836 : ! Downward propagation
837 : ! adjustment factor is obtained as a solution of 19 of the Tang paper
838 : ! eqvivalent to Eq.21 of the Tang paper but for linear-in-tau source
839 : ! radn_dn_Jac(:,nlay+1) = 0._wp
840 0 : do ilev = nlay, 1, -1
841 : ! radn_dn_Jac(:,ilev) = trans(:,ilev)*radn_dn_Jac(:,ilev+1)
842 0 : do icol=1,ncol
843 0 : adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_up(icol,ilev) - &
844 0 : trans(icol,ilev)*source_up(icol,ilev) - source_dn(icol,ilev) )
845 0 : radn_dn(icol,ilev) = trans(icol,ilev)*radn_dn(icol,ilev+1) + source_dn(icol,ilev) + &
846 0 : adjustmentFactor
847 : ! adjustmentFactor = Cn(icol,ilev)*An(icol,ilev)*radn_up_Jac(icol,ilev)
848 : ! radn_dn_Jac(icol,ilev) = radn_dn_Jac(icol,ilev) + adjustmentFactor
849 : enddo
850 : end do
851 : end if
852 0 : end subroutine lw_transport_1rescl
853 : ! -------------------------------------------------------------------------------------------------
854 : !
855 : ! Longwave two-stream solutions to diffuse reflectance and transmittance for a layer
856 : ! with optical depth tau, single scattering albedo w0, and asymmetery parameter g.
857 : !
858 : ! Equations are developed in Meador and Weaver, 1980,
859 : ! doi:10.1175/1520-0469(1980)037<0630:TSATRT>2.0.CO;2
860 : !
861 : ! -------------------------------------------------------------------------------------------------
862 0 : pure subroutine lw_two_stream(ncol, nlay, tau, w0, g, &
863 0 : gamma1, gamma2, Rdif, Tdif)
864 : integer, intent(in) :: ncol, nlay
865 : real(wp), dimension(ncol,nlay), intent(in) :: tau, w0, g
866 : real(wp), dimension(ncol,nlay), intent(out) :: gamma1, gamma2, Rdif, Tdif
867 :
868 : ! -----------------------
869 : integer :: i, j
870 :
871 : ! Variables used in Meador and Weaver
872 0 : real(wp) :: k(ncol)
873 :
874 : ! Ancillary variables
875 0 : real(wp) :: RT_term(ncol)
876 0 : real(wp) :: exp_minusktau(ncol), exp_minus2ktau(ncol)
877 :
878 : real(wp), parameter :: LW_diff_sec = 1.66 ! 1./cos(diffusivity angle)
879 : ! ---------------------------------
880 0 : do j = 1, nlay
881 0 : do i = 1, ncol
882 : !
883 : ! Coefficients differ from SW implementation because the phase function is more isotropic
884 : ! Here we follow Fu et al. 1997, doi:10.1175/1520-0469(1997)054<2799:MSPITI>2.0.CO;2
885 : ! and use a diffusivity sec of 1.66
886 : !
887 0 : gamma1(i,j)= LW_diff_sec * (1._wp - 0.5_wp * w0(i,j) * (1._wp + g(i,j))) ! Fu et al. Eq 2.9
888 0 : gamma2(i,j)= LW_diff_sec * 0.5_wp * w0(i,j) * (1._wp - g(i,j)) ! Fu et al. Eq 2.10
889 : ! Eq 18; k = SQRT(gamma1**2 - gamma2**2), limited below to avoid div by 0.
890 : ! k = 0 for isotropic, conservative scattering; this lower limit on k
891 : ! gives relative error with respect to conservative solution
892 : ! of < 0.1% in Rdif down to tau = 10^-9
893 0 : k(i) = sqrt(max((gamma1(i,j) - gamma2(i,j)) * (gamma1(i,j) + gamma2(i,j)), 1.e-12_wp))
894 : end do
895 :
896 : ! Written to encourage vectorization of exponential
897 0 : exp_minusktau(1:ncol) = exp(-tau(1:ncol,j)*k(1:ncol))
898 :
899 : !
900 : ! Diffuse reflection and transmission
901 : !
902 0 : do i = 1, ncol
903 0 : exp_minus2ktau(i) = exp_minusktau(i) * exp_minusktau(i)
904 :
905 : ! Refactored to avoid rounding errors when k, gamma1 are of very different magnitudes
906 : RT_term(i) = 1._wp / (k (i ) * (1._wp + exp_minus2ktau(i)) + &
907 0 : gamma1(i,j) * (1._wp - exp_minus2ktau(i)) )
908 :
909 : ! Equation 25
910 0 : Rdif(i,j) = RT_term(i) * gamma2(i,j) * (1._wp - exp_minus2ktau(i))
911 :
912 : ! Equation 26
913 0 : Tdif(i,j) = RT_term(i) * 2._wp * k(i) * exp_minusktau(i)
914 : end do
915 :
916 : end do
917 0 : end subroutine lw_two_stream
918 : ! -------------------------------------------------------------------------------------------------
919 : !
920 : ! Source function combination
921 : ! RRTMGP provides two source functions at each level
922 : ! using the spectral mapping from each of the adjascent layers.
923 : ! Need to combine these for use in two-stream calculation.
924 : !
925 : ! -------------------------------------------------------------------------------------------------
926 0 : subroutine lw_combine_sources(ncol, nlay, top_at_1, &
927 0 : lev_src_inc, lev_src_dec, lev_source)
928 : integer, intent(in ) :: ncol, nlay
929 : logical(wl), intent(in ) :: top_at_1
930 : real(wp), dimension(ncol, nlay ), intent(in ) :: lev_src_inc, lev_src_dec
931 : real(wp), dimension(ncol, nlay+1), intent(out) :: lev_source
932 :
933 : integer :: icol, ilay
934 : ! ---------------------------------------------------------------
935 0 : ilay = 1
936 0 : do icol = 1,ncol
937 0 : lev_source(icol, ilay) = lev_src_dec(icol, ilay)
938 : end do
939 0 : do ilay = 2, nlay
940 0 : do icol = 1,ncol
941 0 : lev_source(icol, ilay) = sqrt(lev_src_dec(icol, ilay) * &
942 0 : lev_src_inc(icol, ilay-1))
943 : end do
944 : end do
945 0 : ilay = nlay+1
946 0 : do icol = 1,ncol
947 0 : lev_source(icol, ilay) = lev_src_inc(icol, ilay-1)
948 : end do
949 :
950 0 : end subroutine lw_combine_sources
951 : ! ---------------------------------------------------------------
952 : !
953 : ! Compute LW source function for upward and downward emission at levels using linear-in-tau assumption
954 : ! This version straight from ECRAD
955 : ! Source is provided as W/m2-str; factor of pi converts to flux units
956 : !
957 : ! ---------------------------------------------------------------
958 0 : subroutine lw_source_2str(ncol, nlay, top_at_1, &
959 0 : sfc_emis, sfc_src, &
960 0 : lay_source, lev_source, &
961 0 : gamma1, gamma2, rdif, tdif, tau, source_dn, source_up, source_sfc) &
962 : bind (C, name="rte_lw_source_2str")
963 : integer, intent(in) :: ncol, nlay
964 : logical(wl), intent(in) :: top_at_1
965 : real(wp), dimension(ncol ), intent(in) :: sfc_emis, sfc_src
966 : real(wp), dimension(ncol, nlay), intent(in) :: lay_source, & ! Planck source at layer center
967 : tau, & ! Optical depth (tau)
968 : gamma1, gamma2,& ! Coupling coefficients
969 : rdif, tdif ! Layer reflectance and transmittance
970 : real(wp), dimension(ncol, nlay+1), target, &
971 : intent(in) :: lev_source ! Planck source at layer edges
972 : real(wp), dimension(ncol, nlay), intent(out) :: source_dn, source_up
973 : real(wp), dimension(ncol ), intent(out) :: source_sfc ! Source function for upward radation at surface
974 :
975 : integer :: icol, ilay
976 : real(wp) :: Z, Zup_top, Zup_bottom, Zdn_top, Zdn_bottom
977 0 : real(wp), dimension(:), pointer :: lev_source_bot, lev_source_top
978 : ! ---------------------------------------------------------------
979 0 : do ilay = 1, nlay
980 0 : if(top_at_1) then
981 0 : lev_source_top => lev_source(:,ilay)
982 0 : lev_source_bot => lev_source(:,ilay+1)
983 : else
984 0 : lev_source_top => lev_source(:,ilay+1)
985 0 : lev_source_bot => lev_source(:,ilay)
986 : end if
987 0 : do icol = 1, ncol
988 0 : if (tau(icol,ilay) > 1.0e-8_wp) then
989 : !
990 : ! Toon et al. (JGR 1989) Eqs 26-27
991 : !
992 0 : Z = (lev_source_bot(icol)-lev_source_top(icol)) / (tau(icol,ilay)*(gamma1(icol,ilay)+gamma2(icol,ilay)))
993 0 : Zup_top = Z + lev_source_top(icol)
994 0 : Zup_bottom = Z + lev_source_bot(icol)
995 0 : Zdn_top = -Z + lev_source_top(icol)
996 0 : Zdn_bottom = -Z + lev_source_bot(icol)
997 0 : source_up(icol,ilay) = pi * (Zup_top - rdif(icol,ilay) * Zdn_top - tdif(icol,ilay) * Zup_bottom)
998 0 : source_dn(icol,ilay) = pi * (Zdn_bottom - rdif(icol,ilay) * Zup_bottom - tdif(icol,ilay) * Zdn_top)
999 : else
1000 0 : source_up(icol,ilay) = 0._wp
1001 0 : source_dn(icol,ilay) = 0._wp
1002 : end if
1003 : end do
1004 : end do
1005 0 : do icol = 1, ncol
1006 0 : source_sfc(icol) = pi * sfc_emis(icol) * sfc_src(icol)
1007 : end do
1008 0 : end subroutine lw_source_2str
1009 : ! -------------------------------------------------------------------------------------------------
1010 : !
1011 : ! Lower-level shortwave kernels
1012 : !
1013 : ! -------------------------------------------------------------------------------------------------
1014 : !
1015 : ! Two-stream solutions to diffuse reflectance and transmittance for a layer
1016 : ! with optical depth tau, single scattering albedo w0, and asymmetery parameter g.
1017 : ! Direct reflectance and transmittance used to compute direct beam source for diffuse radiation
1018 : ! in layers and at surface; report direct beam as a byproduct
1019 : ! Computing the direct-beam source for diffuse radiation at the same time as R and T for
1020 : ! direct radiation reduces memory traffic and use.
1021 : !
1022 : ! Equations are developed in Meador and Weaver, 1980,
1023 : ! doi:10.1175/1520-0469(1980)037<0630:TSATRT>2.0.CO;2
1024 : !
1025 : ! -------------------------------------------------------------------------------------------------
1026 87194912 : pure subroutine sw_dif_and_source(ncol, nlay, top_at_1, mu0, sfc_albedo, &
1027 87194912 : tau, w0, g, &
1028 87194912 : Rdif, Tdif, source_dn, source_up, source_sfc, flux_dn_dir) bind (C, name="rte_sw_source_dir")
1029 : integer, intent(in ) :: ncol, nlay
1030 : logical(wl), intent(in ) :: top_at_1
1031 : real(wp), dimension(ncol ), intent(in ) :: sfc_albedo ! surface albedo for direct radiation
1032 : real(wp), dimension(ncol,nlay ), intent(in ) :: tau, w0, g, mu0
1033 : real(wp), dimension(ncol,nlay ), intent( out) :: Rdif, Tdif, source_dn, source_up
1034 : real(wp), dimension(ncol ), intent( out) :: source_sfc ! Source function for upward radation at surface
1035 : real(wp), dimension(ncol,nlay+1), target, &
1036 : intent(inout) :: flux_dn_dir ! Direct beam flux
1037 :
1038 : ! -----------------------
1039 : integer :: i, j
1040 :
1041 : ! Variables used in Meador and Weaver
1042 : real(wp) :: gamma1, gamma2, gamma3, gamma4, alpha1, alpha2
1043 :
1044 :
1045 : ! Ancillary variables
1046 : real(wp), parameter :: min_k = 1.e4_wp * epsilon(1._wp) ! Suggestion from Chiel van Heerwaarden
1047 : real(wp) :: k, exp_minusktau, k_mu, k_gamma3, k_gamma4
1048 : real(wp) :: RT_term, exp_minus2ktau
1049 : real(wp) :: Rdir, Tdir, Tnoscat
1050 87194912 : real(wp), pointer, dimension(:) :: dir_flux_inc, dir_flux_trans
1051 : integer :: lay_index
1052 : real(wp) :: tau_s, w0_s, g_s, mu0_s
1053 : ! ---------------------------------
1054 :
1055 8283516640 : do j = 1, nlay
1056 8196321728 : if(top_at_1) then
1057 8196321728 : lay_index = j
1058 8196321728 : dir_flux_inc => flux_dn_dir(:,lay_index )
1059 8196321728 : dir_flux_trans => flux_dn_dir(:,lay_index+1)
1060 : else
1061 0 : lay_index = nlay-j+1
1062 0 : dir_flux_inc => flux_dn_dir(:,lay_index+1)
1063 0 : dir_flux_trans => flux_dn_dir(:,lay_index )
1064 : end if
1065 :
1066 >13159*10^7 : do i = 1, ncol
1067 : !
1068 : ! Scalars
1069 : !
1070 >12331*10^7 : tau_s = tau(i, lay_index)
1071 >12331*10^7 : w0_s = w0 (i, lay_index)
1072 >12331*10^7 : g_s = g (i, lay_index)
1073 >12331*10^7 : mu0_s = mu0(i, lay_index)
1074 : !
1075 : ! Zdunkowski Practical Improved Flux Method "PIFM"
1076 : ! (Zdunkowski et al., 1980; Contributions to Atmospheric Physics 53, 147-66)
1077 : !
1078 >12331*10^7 : gamma1 = (8._wp - w0_s * (5._wp + 3._wp * g_s)) * .25_wp
1079 >12331*10^7 : gamma2 = 3._wp *(w0_s * (1._wp - g_s)) * .25_wp
1080 : !
1081 : ! Direct reflect and transmission
1082 : !
1083 : ! Eq 18; k = SQRT(gamma1**2 - gamma2**2), limited below to avoid div by 0.
1084 : ! k = 0 for isotropic, conservative scattering; this lower limit on k
1085 : ! gives relative error with respect to conservative solution
1086 : ! of < 0.1% in Rdif down to tau = 10^-9
1087 >12331*10^7 : k = sqrt(max((gamma1 - gamma2) * (gamma1 + gamma2), min_k))
1088 >12331*10^7 : exp_minusktau = exp(-tau_s*k)
1089 >12331*10^7 : exp_minus2ktau = exp_minusktau * exp_minusktau
1090 :
1091 : ! Refactored to avoid rounding errors when k, gamma1 are of very different magnitudes
1092 : RT_term = 1._wp / (k * (1._wp + exp_minus2ktau) + &
1093 >12331*10^7 : gamma1 * (1._wp - exp_minus2ktau) )
1094 : ! Equation 25
1095 >12331*10^7 : Rdif(i,lay_index) = RT_term * gamma2 * (1._wp - exp_minus2ktau)
1096 :
1097 : ! Equation 26
1098 >12331*10^7 : Tdif(i,lay_index) = RT_term * 2._wp * k * exp_minusktau
1099 :
1100 : !
1101 : ! On a round earth, where mu0 can increase with depth in the atmosphere,
1102 : ! levels with mu0 <= 0 have no direct beam and hence no source for diffuse light
1103 : !
1104 >13150*10^7 : if(mu0_s > 0._wp) then
1105 >12331*10^7 : k_mu = k * mu0_s
1106 : !
1107 : ! Equation 14, multiplying top and bottom by exp(-k*tau)
1108 : ! and rearranging to avoid div by 0.
1109 : !
1110 : RT_term = w0_s * RT_term/merge(1._wp - k_mu*k_mu, &
1111 : epsilon(1._wp), &
1112 >12331*10^7 : abs(1._wp - k_mu*k_mu) >= epsilon(1._wp))
1113 : !
1114 : ! Zdunkowski Practical Improved Flux Method "PIFM"
1115 : ! (Zdunkowski et al., 1980; Contributions to Atmospheric Physics 53, 147-66)
1116 : !
1117 >12331*10^7 : gamma3 = (2._wp - 3._wp * mu0_s * g_s ) * .25_wp
1118 >12331*10^7 : gamma4 = 1._wp - gamma3
1119 >12331*10^7 : alpha1 = gamma1 * gamma4 + gamma2 * gamma3 ! Eq. 16
1120 >12331*10^7 : alpha2 = gamma1 * gamma3 + gamma2 * gamma4 ! Eq. 17
1121 :
1122 : !
1123 : ! Transmittance of direct, unscattered beam.
1124 : !
1125 >12331*10^7 : k_gamma3 = k * gamma3
1126 >12331*10^7 : k_gamma4 = k * gamma4
1127 >12331*10^7 : Tnoscat = exp(-tau_s/mu0_s)
1128 : Rdir = RT_term * &
1129 : ((1._wp - k_mu) * (alpha2 + k_gamma3) - &
1130 : (1._wp + k_mu) * (alpha2 - k_gamma3) * exp_minus2ktau - &
1131 >12331*10^7 : 2.0_wp * (k_gamma3 - alpha2 * k_mu) * exp_minusktau * Tnoscat)
1132 : !
1133 : ! Equation 15, multiplying top and bottom by exp(-k*tau),
1134 : ! multiplying through by exp(-tau/mu0) to
1135 : ! prefer underflow to overflow
1136 : ! Omitting direct transmittance
1137 : !
1138 : Tdir = -RT_term * &
1139 : ((1._wp + k_mu) * (alpha1 + k_gamma4) * Tnoscat - &
1140 : (1._wp - k_mu) * (alpha1 - k_gamma4) * exp_minus2ktau * Tnoscat - &
1141 >12331*10^7 : 2.0_wp * (k_gamma4 + alpha1 * k_mu) * exp_minusktau)
1142 : ! Final check that energy is not spuriously created, by recognizing that
1143 : ! the beam can either be reflected, penetrate unscattered to the base of a layer,
1144 : ! or penetrate through but be scattered on the way - the rest is absorbed
1145 : ! Makes the equations safer in single precision. Credit: Robin Hogan, Peter Ukkonen
1146 >12331*10^7 : Rdir = max(0.0_wp, min(Rdir, (1.0_wp - Tnoscat ) ))
1147 >12331*10^7 : Tdir = max(0.0_wp, min(Tdir, (1.0_wp - Tnoscat - Rdir) ))
1148 :
1149 >12331*10^7 : source_up(i,lay_index) = Rdir * dir_flux_inc(i)
1150 >12331*10^7 : source_dn(i,lay_index) = Tdir * dir_flux_inc(i)
1151 >12331*10^7 : dir_flux_trans(i) = Tnoscat * dir_flux_inc(i)
1152 : else
1153 0 : source_up(i,lay_index) = 0._wp
1154 0 : source_dn(i,lay_index) = 0._wp
1155 0 : dir_flux_trans(i) = 0._wp
1156 : end if
1157 : end do
1158 : end do
1159 1399006112 : source_sfc(:) = dir_flux_trans(:)*sfc_albedo(:)
1160 :
1161 87194912 : end subroutine sw_dif_and_source
1162 : ! ---------------------------------------------------------------
1163 : !
1164 : ! Transport of diffuse radiation through a vertically layered atmosphere.
1165 : ! Equations are after Shonk and Hogan 2008, doi:10.1175/2007JCLI1940.1 (SH08)
1166 : ! This routine is shared by longwave and shortwave
1167 : !
1168 : ! -------------------------------------------------------------------------------------------------
1169 87194912 : subroutine adding(ncol, nlay, top_at_1, &
1170 87194912 : albedo_sfc, &
1171 87194912 : rdif, tdif, &
1172 87194912 : src_dn, src_up, src_sfc, &
1173 87194912 : flux_up, flux_dn)
1174 : integer, intent(in ) :: ncol, nlay
1175 : logical(wl), intent(in ) :: top_at_1
1176 : real(wp), dimension(ncol ), intent(in ) :: albedo_sfc
1177 : real(wp), dimension(ncol,nlay ), intent(in ) :: rdif, tdif
1178 : real(wp), dimension(ncol,nlay ), intent(in ) :: src_dn, src_up
1179 : real(wp), dimension(ncol ), intent(in ) :: src_sfc
1180 : real(wp), dimension(ncol,nlay+1), intent( out) :: flux_up
1181 : ! intent(inout) because top layer includes incident flux
1182 : real(wp), dimension(ncol,nlay+1), intent(inout) :: flux_dn
1183 : ! ------------------
1184 : integer :: ilev
1185 174389824 : real(wp), dimension(ncol,nlay+1) :: albedo, & ! reflectivity to diffuse radiation below this level
1186 : ! alpha in SH08
1187 174389824 : src ! source of diffuse upwelling radiation from emission or
1188 : ! scattering of direct beam
1189 : ! G in SH08
1190 174389824 : real(wp), dimension(ncol,nlay ) :: denom ! beta in SH08
1191 : ! ------------------
1192 : !
1193 : ! Indexing into arrays for upward and downward propagation depends on the vertical
1194 : ! orientation of the arrays (whether the domain top is at the first or last index)
1195 : ! We write the loops out explicitly so compilers will have no trouble optimizing them.
1196 : !
1197 87194912 : if(top_at_1) then
1198 87194912 : ilev = nlay + 1
1199 : ! Albedo of lowest level is the surface albedo...
1200 1399006112 : albedo(:,ilev) = albedo_sfc(:)
1201 : ! ... and source of diffuse radiation is surface emission
1202 1399006112 : src(:,ilev) = src_sfc(:)
1203 :
1204 : !
1205 : ! From bottom to top of atmosphere --
1206 : ! compute albedo and source of upward radiation
1207 : !
1208 8283516640 : do ilev = nlay, 1, -1
1209 >13150*10^7 : denom(:, ilev) = 1._wp/(1._wp - rdif(:,ilev)*albedo(:,ilev+1)) ! Eq 10
1210 8196321728 : albedo(:,ilev) = rdif(:,ilev) + &
1211 >13150*10^7 : tdif(:,ilev)*tdif(:,ilev) * albedo(:,ilev+1) * denom(:,ilev) ! Equation 9
1212 : !
1213 : ! Equation 11 -- source is emitted upward radiation at top of layer plus
1214 : ! radiation emitted at bottom of layer,
1215 : ! transmitted through the layer and reflected from layers below (tdiff*src*albedo)
1216 : !
1217 : src(:,ilev) = src_up(:, ilev) + &
1218 : tdif(:,ilev) * denom(:,ilev) * &
1219 >13159*10^7 : (src(:,ilev+1) + albedo(:,ilev+1)*src_dn(:,ilev))
1220 : end do
1221 :
1222 : ! Eq 12, at the top of the domain upwelling diffuse is due to ...
1223 1399006112 : ilev = 1
1224 : flux_up(:,ilev) = flux_dn(:,ilev) * albedo(:,ilev) + & ! ... reflection of incident diffuse and
1225 1399006112 : src(:,ilev) ! emission from below
1226 :
1227 : !
1228 : ! From the top of the atmosphere downward -- compute fluxes
1229 : !
1230 8283516640 : do ilev = 2, nlay+1
1231 8196321728 : flux_dn(:,ilev) = (tdif(:,ilev-1)*flux_dn(:,ilev-1) + & ! Equation 13
1232 8196321728 : rdif(:,ilev-1)*src(:,ilev) + &
1233 >14789*10^7 : src_dn(:,ilev-1)) * denom(:,ilev-1)
1234 : flux_up(:,ilev) = flux_dn(:,ilev) * albedo(:,ilev) + & ! Equation 12
1235 >13159*10^7 : src(:,ilev)
1236 : end do
1237 : else
1238 0 : ilev = 1
1239 : ! Albedo of lowest level is the surface albedo...
1240 0 : albedo(:,ilev) = albedo_sfc(:)
1241 : ! ... and source of diffuse radiation is surface emission
1242 0 : src(:,ilev) = src_sfc(:)
1243 :
1244 : !
1245 : ! From bottom to top of atmosphere --
1246 : ! compute albedo and source of upward radiation
1247 : !
1248 0 : do ilev = 1, nlay
1249 0 : denom(:, ilev ) = 1._wp/(1._wp - rdif(:,ilev)*albedo(:,ilev)) ! Eq 10
1250 0 : albedo(:,ilev+1) = rdif(:,ilev) + &
1251 0 : tdif(:,ilev)*tdif(:,ilev) * albedo(:,ilev) * denom(:,ilev) ! Equation 9
1252 : !
1253 : ! Equation 11 -- source is emitted upward radiation at top of layer plus
1254 : ! radiation emitted at bottom of layer,
1255 : ! transmitted through the layer and reflected from layers below (tdiff*src*albedo)
1256 : !
1257 : src(:,ilev+1) = src_up(:, ilev) + &
1258 : tdif(:,ilev) * denom(:,ilev) * &
1259 0 : (src(:,ilev) + albedo(:,ilev)*src_dn(:,ilev))
1260 : end do
1261 :
1262 : ! Eq 12, at the top of the domain upwelling diffuse is due to ...
1263 0 : ilev = nlay+1
1264 : flux_up(:,ilev) = flux_dn(:,ilev) * albedo(:,ilev) + & ! ... reflection of incident diffuse and
1265 0 : src(:,ilev) ! scattering by the direct beam below
1266 :
1267 : !
1268 : ! From the top of the atmosphere downward -- compute fluxes
1269 : !
1270 0 : do ilev = nlay, 1, -1
1271 0 : flux_dn(:,ilev) = (tdif(:,ilev)*flux_dn(:,ilev+1) + & ! Equation 13
1272 0 : rdif(:,ilev)*src(:,ilev) + &
1273 0 : src_dn(:, ilev)) * denom(:,ilev)
1274 : flux_up(:,ilev) = flux_dn(:,ilev) * albedo(:,ilev) + & ! Equation 12
1275 0 : src(:,ilev)
1276 :
1277 : end do
1278 : end if
1279 87194912 : end subroutine adding
1280 : end module mo_rte_solver_kernels
|
