Home

Context Navigation

radiation_model_mod.f90 @ 4392

Last change on this file since 4392 was 4392, checked in by pavelkrc, 5 years ago

Merge branch resler (updated documentation, optional flux tracing, bugfix)

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

/palm/branches/chemistry/SOURCE/radiation_model_mod.f90	2047-3190,3218-3297
/palm/branches/forwind/SOURCE/radiation_model_mod.f90	1564-1913
/palm/branches/mosaik_M2/radiation_model_mod.f90	2360-3471
/palm/branches/palm4u/SOURCE/radiation_model_mod.f90	2540-2692
/palm/branches/radiation/SOURCE/radiation_model_mod.f90	2081-3493
/palm/branches/rans/SOURCE/radiation_model_mod.f90	2078-3128
/palm/branches/resler/SOURCE/radiation_model_mod.f90	2023-4391
/palm/branches/salsa/SOURCE/radiation_model_mod.f90	2503-3460
/palm/branches/fricke/SOURCE/radiation_model_mod.f90	942-977
/palm/branches/hoffmann/SOURCE/radiation_model_mod.f90	989-1052
/palm/branches/letzel/masked_output/SOURCE/radiation_model_mod.f90	296-409
/palm/branches/suehring/radiation_model_mod.f90	423-666

File size: 554.7 KB

Line
1	!> @file radiation_model_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 2015-2020 Institute of Computer Science of the
18	! Czech Academy of Sciences, Prague
19	! Copyright 2015-2019 Czech Technical University in Prague
20	! Copyright 1997-2020 Leibniz Universitaet Hannover
21	!------------------------------------------------------------------------------!
22	!
23	! Current revisions:
24	! ------------------
25	!
26	!
27	! Former revisions:
28	! -----------------
29	! $Id: radiation_model_mod.f90 4392 2020-01-31 16:14:57Z pavelkrc $
30	! - Add debug tracing of large radiative fluxes (option trace_fluxes_above)
31	! - Print exact counts of SVF and CSF if debut_output is enabled
32	! - Update descriptions of RTM 3.0 and related comments
33	!
34	! 4360 2020-01-07 11:25:50Z suehring
35	! Renamed pc_heating_rate, pc_transpiration_rate, pc_transpiration_rate to
36	! pcm_heating_rate, pcm_latent_rate, pcm_transpiration_rate
37	!
38	! 4340 2019-12-16 08:17:03Z Giersch
39	! Albedo indices for building_surface_pars are now declared as parameters to
40	! prevent an error if the gfortran compiler with -Werror=unused-value is used
41	!
42	! 4291 2019-11-11 12:36:54Z moh.hefny
43	! Enabled RTM in case of biometeorology even if there is no vertical
44	! surfaces or 3D vegetation in the domain
45	!
46	! 4286 2019-10-30 16:01:14Z resler
47	! - Fix wrong treating of time_rad during interpolation in radiation_model_mod
48	! - Fix wrong checks of time_rad from dynamic driver in radiation_model_mod
49	! - Add new directional model of human body for MRT: ellipsoid
50	!
51	! 4271 2019-10-23 10:46:41Z maronga
52	! Bugfix: missing parentheses in calculation of snow albedo
53	!
54	! 4245 2019-09-30 08:40:37Z pavelkrc
55	! Initialize explicit per-surface albedos from building_surface_pars
56	!
57	! 4238 2019-09-25 16:06:01Z suehring
58	! Modify check in order to avoid equality comparisons of floating points
59	!
60	! 4227 2019-09-10 18:04:34Z gronemeier
61	! implement new palm_date_time_mod
62	!
63	! 4226 2019-09-10 17:03:24Z suehring
64	! - Netcdf input routine for dimension length renamed
65	! - Define time variable for external radiation input relative to time_utc_init
66	!
67	! 4210 2019-09-02 13:07:09Z suehring
68	! - Revise steering of splitting diffuse and direct radiation
69	! - Bugfixes in checks
70	! - Optimize mapping of radiation components onto 2D arrays, avoid unnecessary
71	! operations
72	!
73	! 4208 2019-09-02 09:01:07Z suehring
74	! Bugfix in accessing albedo_pars in the clear-sky branch
75	! (merge from branch resler)
76	!
77	! 4198 2019-08-29 15:17:48Z gronemeier
78	! Prohibit execution of radiation model if rotation_angle is not zero
79	!
80	! 4197 2019-08-29 14:33:32Z suehring
81	! Revise steering of surface albedo initialization when albedo_pars is provided
82	!
83	! 4190 2019-08-27 15:42:37Z suehring
84	! Implement external radiation forcing also for level-of-detail = 2
85	! (horizontally 2D radiation)
86	!
87	! 4188 2019-08-26 14:15:47Z suehring
88	! Minor adjustment in error message
89	!
90	! 4187 2019-08-26 12:43:15Z suehring
91	! - Take external radiation from root domain dynamic input if not provided for
92	! each nested domain
93	! - Combine MPI_ALLREDUCE calls to reduce mpi overhead
94	!
95	! 4182 2019-08-22 15:20:23Z scharf
96	! Corrected "Former revisions" section
97	!
98	! 4179 2019-08-21 11:16:12Z suehring
99	! Remove debug prints
100	!
101	! 4178 2019-08-21 11:13:06Z suehring
102	! External radiation forcing implemented.
103	!
104	! 4168 2019-08-16 13:50:17Z suehring
105	! Replace function get_topography_top_index by topo_top_ind
106	!
107	! 4157 2019-08-14 09:19:12Z suehring
108	! Give informative message on raytracing distance only by core zero
109	!
110	! 4148 2019-08-08 11:26:00Z suehring
111	! Comments added
112	!
113	! 4134 2019-08-02 18:39:57Z suehring
114	! Bugfix in formatted write statement
115	!
116	! 4127 2019-07-30 14:47:10Z suehring
117	! Remove unused pch_index (merge from branch resler)
118	!
119	! 4089 2019-07-11 14:30:27Z suehring
120	! - Correct level 2 initialization of spectral albedos in rrtmg branch, long- and
121	! shortwave albedos were mixed-up.
122	! - Change order of albedo_pars so that it is now consistent with the defined
123	! order of albedo_pars in PIDS
124	!
125	! 4069 2019-07-01 14:05:51Z Giersch
126	! Masked output running index mid has been introduced as a local variable to
127	! avoid runtime error (Loop variable has been modified) in time_integration
128	!
129	! 4067 2019-07-01 13:29:25Z suehring
130	! Bugfix, pass dummy string to MPI_INFO_SET (J. Resler)
131	!
132	! 4039 2019-06-18 10:32:41Z suehring
133	! Bugfix for masked data output
134	!
135	! 4008 2019-05-30 09:50:11Z moh.hefny
136	! Bugfix in check variable when a variable's string is less than 3
137	! characters is processed. All variables now are checked if they
138	! belong to radiation
139	!
140	! 3992 2019-05-22 16:49:38Z suehring
141	! Bugfix in rrtmg radiation branch in a nested run when the lowest prognistic
142	! grid points in a child domain are all inside topography
143	!
144	! 3987 2019-05-22 09:52:13Z kanani
145	! Introduce alternative switch for debug output during timestepping
146	!
147	! 3943 2019-05-02 09:50:41Z maronga
148	! Missing blank characteer added.
149	!
150	! 3900 2019-04-16 15:17:43Z suehring
151	! Fixed initialization problem
152	!
153	! 3885 2019-04-11 11:29:34Z kanani
154	! Changes related to global restructuring of location messages and introduction
155	! of additional debug messages
156	!
157	! 3881 2019-04-10 09:31:22Z suehring
158	! Output of albedo and emissivity moved from USM, bugfixes in initialization
159	! of albedo
160	!
161	! 3861 2019-04-04 06:27:41Z maronga
162	! Bugfix: factor of 4.0 required instead of 3.0 in calculation of rad_lw_out_change_0
163	!
164	! 3859 2019-04-03 20:30:31Z maronga
165	! Added some descriptions
166	!
167	! 3847 2019-04-01 14:51:44Z suehring
168	! Implement check for dt_radiation (must be > 0)
169	!
170	! 3846 2019-04-01 13:55:30Z suehring
171	! unused variable removed
172	!
173	! 3814 2019-03-26 08:40:31Z pavelkrc
174	! Change zenith(0:0) and others to scalar.
175	! Code review.
176	! Rename exported nzu, nzp and related variables due to name conflict
177	!
178	! 3771 2019-02-28 12:19:33Z raasch
179	! rrtmg preprocessor for directives moved/added, save attribute added to temporary
180	! pointers to avoid compiler warnings about outlived pointer targets,
181	! statement added to avoid compiler warning about unused variable
182	!
183	! 3769 2019-02-28 10:16:49Z moh.hefny
184	! removed unused variables and subroutine radiation_radflux_gridbox
185	!
186	! 3767 2019-02-27 08:18:02Z raasch
187	! unused variable for file index removed from rrd-subroutines parameter list
188	!
189	! 3760 2019-02-21 18:47:35Z moh.hefny
190	! Bugfix: initialized simulated_time before calculating solar position
191	! to enable restart option with reading in SVF from file(s).
192	!
193	! 3754 2019-02-19 17:02:26Z kanani
194	! (resler, pavelkrc)
195	! Bugfixes: add further required MRT factors to read/write_svf,
196	! fix for aggregating view factors to eliminate local noise in reflected
197	! irradiance at mutually close surfaces (corners, presence of trees) in the
198	! angular discretization scheme.
199	!
200	! 3752 2019-02-19 09:37:22Z resler
201	! added read/write number of MRT factors to the respective routines
202	!
203	! 3705 2019-01-29 19:56:39Z suehring
204	! Make variables that are sampled in virtual measurement module public
205	!
206	! 3704 2019-01-29 19:51:41Z suehring
207	! Some interface calls moved to module_interface + cleanup
208	!
209	! 3667 2019-01-10 14:26:24Z schwenkel
210	! Modified check for rrtmg input files
211	!
212	! 3655 2019-01-07 16:51:22Z knoop
213	! nopointer option removed
214	!
215	! 1496 2014-12-02 17:25:50Z maronga
216	! Initial revision
217	!
218	!
219	! Description:
220	! ------------
221	!> Radiation models and interfaces:
222	!> constant, simple and RRTMG models, interface to external radiation model
223	!> Radiative Transfer Model (RTM) version 3.0 for modelling of radiation
224	!> interactions within urban canopy or other surface layer in complex terrain
225	!> Integrations of RTM with other PALM-4U modules:
226	!> integration with RRTMG, USM, LSM, PCM, BIO modules
227	!>
228	!> @todo move variable definitions used in radiation_init only to the subroutine
229	!> as they are no longer required after initialization.
230	!> @todo Output of full column vertical profiles used in RRTMG
231	!> @todo Output of other rrtm arrays (such as volume mixing ratios)
232	!> @todo Optimize radiation_tendency routines
233	!>
234	!> @note Many variables have a leading dummy dimension (0:0) in order to
235	!> match the assume-size shape expected by the RRTMG model.
236	!------------------------------------------------------------------------------!
237	MODULE radiation_model_mod
238
239	USE arrays_3d, &
240	ONLY: dzw, hyp, nc, pt, p, q, ql, u, v, w, zu, zw, exner, d_exner
241
242	USE basic_constants_and_equations_mod, &
243	ONLY: c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, sigma_sb, &
244	barometric_formula
245
246	USE calc_mean_profile_mod, &
247	ONLY: calc_mean_profile
248
249	USE control_parameters, &
250	ONLY: biometeorology, cloud_droplets, coupling_char, &
251	debug_output, debug_output_timestep, debug_string, &
252	dt_3d, &
253	dz, dt_spinup, end_time, &
254	humidity, &
255	initializing_actions, io_blocks, io_group, &
256	land_surface, large_scale_forcing, &
257	latitude, longitude, lsf_surf, &
258	message_string, plant_canopy, pt_surface, &
259	rho_surface, simulated_time, spinup_time, surface_pressure, &
260	read_svf, write_svf, &
261	time_since_reference_point, urban_surface, varnamelength
262
263	USE cpulog, &
264	ONLY: cpu_log, log_point, log_point_s
265
266	USE grid_variables, &
267	ONLY: ddx, ddy, dx, dy
268
269	USE indices, &
270	ONLY: nnx, nny, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, &
271	nzb, nzt, topo_top_ind
272
273	USE, INTRINSIC :: iso_c_binding
274
275	USE kinds
276
277	USE bulk_cloud_model_mod, &
278	ONLY: bulk_cloud_model, microphysics_morrison, na_init, nc_const, sigma_gc
279
280	#if defined ( __netcdf )
281	USE NETCDF
282	#endif
283
284	USE netcdf_data_input_mod, &
285	ONLY: albedo_type_f, &
286	albedo_pars_f, &
287	building_type_f, &
288	building_surface_pars_f, &
289	pavement_type_f, &
290	vegetation_type_f, &
291	water_type_f, &
292	char_fill, &
293	char_lod, &
294	check_existence, &
295	close_input_file, &
296	get_attribute, &
297	get_dimension_length, &
298	get_variable, &
299	inquire_num_variables, &
300	inquire_variable_names, &
301	input_file_dynamic, &
302	input_pids_dynamic, &
303	num_var_pids, &
304	pids_id, &
305	open_read_file, &
306	real_1d_3d, &
307	vars_pids
308
309	USE palm_date_time_mod, &
310	ONLY: date_time_str_len, get_date_time, &
311	hours_per_day, seconds_per_hour
312
313	USE plant_canopy_model_mod, &
314	ONLY: lad_s, &
315	pcm_heating_rate, &
316	pcm_transpiration_rate, &
317	pcm_latent_rate, &
318	plant_canopy_transpiration, &
319	pcm_calc_transpiration_rate
320
321	USE pegrid
322
323	#if defined ( __rrtmg )
324	USE parrrsw, &
325	ONLY: naerec, nbndsw
326
327	USE parrrtm, &
328	ONLY: nbndlw
329
330	USE rrtmg_lw_init, &
331	ONLY: rrtmg_lw_ini
332
333	USE rrtmg_sw_init, &
334	ONLY: rrtmg_sw_ini
335
336	USE rrtmg_lw_rad, &
337	ONLY: rrtmg_lw
338
339	USE rrtmg_sw_rad, &
340	ONLY: rrtmg_sw
341	#endif
342	USE statistics, &
343	ONLY: hom
344
345	USE surface_mod, &
346	ONLY: ind_pav_green, ind_veg_wall, ind_wat_win, &
347	surf_lsm_h, surf_lsm_v, surf_type, surf_usm_h, surf_usm_v, &
348	vertical_surfaces_exist
349
350	IMPLICIT NONE
351
352	CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
353
354	!
355	!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
356	CHARACTER(37), DIMENSION(0:33), PARAMETER :: albedo_type_name = (/ &
357	'user defined ', & ! 0
358	'ocean ', & ! 1
359	'mixed farming, tall grassland ', & ! 2
360	'tall/medium grassland ', & ! 3
361	'evergreen shrubland ', & ! 4
362	'short grassland/meadow/shrubland ', & ! 5
363	'evergreen needleleaf forest ', & ! 6
364	'mixed deciduous evergreen forest ', & ! 7
365	'deciduous forest ', & ! 8
366	'tropical evergreen broadleaved forest', & ! 9
367	'medium/tall grassland/woodland ', & ! 10
368	'desert, sandy ', & ! 11
369	'desert, rocky ', & ! 12
370	'tundra ', & ! 13
371	'land ice ', & ! 14
372	'sea ice ', & ! 15
373	'snow ', & ! 16
374	'bare soil ', & ! 17
375	'asphalt/concrete mix ', & ! 18
376	'asphalt (asphalt concrete) ', & ! 19
377	'concrete (Portland concrete) ', & ! 20
378	'sett ', & ! 21
379	'paving stones ', & ! 22
380	'cobblestone ', & ! 23
381	'metal ', & ! 24
382	'wood ', & ! 25
383	'gravel ', & ! 26
384	'fine gravel ', & ! 27
385	'pebblestone ', & ! 28
386	'woodchips ', & ! 29
387	'tartan (sports) ', & ! 30
388	'artifical turf (sports) ', & ! 31
389	'clay (sports) ', & ! 32
390	'building (dummy) ' & ! 33
391	/)
392	!
393	!-- Indices of radiation-related input attributes in building_surface_pars
394	!-- (other are in urban_surface_mod)
395	INTEGER(iwp), PARAMETER :: ind_s_alb_b_wall = 19 !< index for Broadband albedo of wall fraction
396	INTEGER(iwp), PARAMETER :: ind_s_alb_l_wall = 20 !< index for Longwave albedo of wall fraction
397	INTEGER(iwp), PARAMETER :: ind_s_alb_s_wall = 21 !< index for Shortwave albedo of wall fraction
398	INTEGER(iwp), PARAMETER :: ind_s_alb_b_win = 22 !< index for Broadband albedo of window fraction
399	INTEGER(iwp), PARAMETER :: ind_s_alb_l_win = 23 !< index for Longwave albedo of window fraction
400	INTEGER(iwp), PARAMETER :: ind_s_alb_s_win = 24 !< index for Shortwave albedo of window fraction
401	INTEGER(iwp), PARAMETER :: ind_s_alb_b_green = 24 !< index for Broadband albedo of green fraction
402	INTEGER(iwp), PARAMETER :: ind_s_alb_l_green = 25 !< index for Longwave albedo of green fraction
403	INTEGER(iwp), PARAMETER :: ind_s_alb_s_green = 26 !< index for Shortwave albedo of green fraction
404
405	INTEGER(iwp) :: albedo_type = 9999999_iwp, & !< Albedo surface type
406	dots_rad = 0_iwp !< starting index for timeseries output
407
408	LOGICAL :: unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
409	constant_albedo = .FALSE., & !< flag parameter indicating whether the albedo may change depending on zenith
410	force_radiation_call = .FALSE., & !< flag parameter for unscheduled radiation calls
411	lw_radiation = .TRUE., & !< flag parameter indicating whether longwave radiation shall be calculated
412	radiation = .FALSE., & !< flag parameter indicating whether the radiation model is used
413	sun_up = .TRUE., & !< flag parameter indicating whether the sun is up or down
414	sw_radiation = .TRUE., & !< flag parameter indicating whether shortwave radiation shall be calculated
415	sun_direction = .FALSE., & !< flag parameter indicating whether solar direction shall be calculated
416	average_radiation = .FALSE., & !< flag to set the calculation of radiation averaging for the domain
417	radiation_interactions = .FALSE., & !< flag to activiate RTM (TRUE only if vertical urban/land surface and trees exist)
418	surface_reflections = .TRUE., & !< flag to switch the calculation of radiation interaction between surfaces.
419	!< When it switched off, only the effect of buildings and trees shadow
420	!< will be considered. However fewer SVFs are expected.
421	radiation_interactions_on = .TRUE. !< namelist flag to force RTM activiation regardless to vertical urban/land surface and trees
422
423	REAL(wp) :: albedo = 9999999.9_wp, & !< NAMELIST alpha
424	albedo_lw_dif = 9999999.9_wp, & !< NAMELIST aldif
425	albedo_lw_dir = 9999999.9_wp, & !< NAMELIST aldir
426	albedo_sw_dif = 9999999.9_wp, & !< NAMELIST asdif
427	albedo_sw_dir = 9999999.9_wp, & !< NAMELIST asdir
428	decl_1, & !< declination coef. 1
429	decl_2, & !< declination coef. 2
430	decl_3, & !< declination coef. 3
431	dt_radiation = 0.0_wp, & !< radiation model timestep
432	emissivity = 9999999.9_wp, & !< NAMELIST surface emissivity
433	lon = 0.0_wp, & !< longitude in radians
434	lat = 0.0_wp, & !< latitude in radians
435	net_radiation = 0.0_wp, & !< net radiation at surface
436	skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
437	sky_trans, & !< sky transmissivity
438	time_radiation = 0.0_wp, & !< time since last call of radiation code
439	trace_fluxes_above = -1.0_wp !< NAMELIST option for debug tracing of large radiative fluxes (W/m2;W/m3)
440
441	INTEGER(iwp) :: day_of_year !< day of the current year
442
443	REAL(wp) :: cos_zenith !< cosine of solar zenith angle, also z-coordinate of solar unit vector
444	REAL(wp) :: d_hours_day !< 1 / hours-per-day
445	REAL(wp) :: d_seconds_hour !< 1 / seconds-per-hour
446	REAL(wp) :: second_of_day !< second of the current day
447	REAL(wp) :: sun_dir_lat !< y-coordinate of solar unit vector
448	REAL(wp) :: sun_dir_lon !< x-coordinate of solar unit vector
449
450	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_net_av !< average of net radiation (rad_net) at surface
451	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_lw_in_xy_av !< average of incoming longwave radiation at surface
452	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_lw_out_xy_av !< average of outgoing longwave radiation at surface
453	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_xy_av !< average of incoming shortwave radiation at surface
454	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_out_xy_av !< average of outgoing shortwave radiation at surface
455
456	REAL(wp), PARAMETER :: emissivity_atm_clsky = 0.8_wp !< emissivity of the clear-sky atmosphere
457	!
458	!-- Land surface albedos for solar zenith angle of 60degree after Briegleb (1992)
459	!-- (broadband, longwave, shortwave ): bb, lw, sw,
460	REAL(wp), DIMENSION(0:2,1:33), PARAMETER :: albedo_pars = RESHAPE( (/&
461	0.06_wp, 0.06_wp, 0.06_wp, & ! 1
462	0.19_wp, 0.28_wp, 0.09_wp, & ! 2
463	0.23_wp, 0.33_wp, 0.11_wp, & ! 3
464	0.23_wp, 0.33_wp, 0.11_wp, & ! 4
465	0.25_wp, 0.34_wp, 0.14_wp, & ! 5
466	0.14_wp, 0.22_wp, 0.06_wp, & ! 6
467	0.17_wp, 0.27_wp, 0.06_wp, & ! 7
468	0.19_wp, 0.31_wp, 0.06_wp, & ! 8
469	0.14_wp, 0.22_wp, 0.06_wp, & ! 9
470	0.18_wp, 0.28_wp, 0.06_wp, & ! 10
471	0.43_wp, 0.51_wp, 0.35_wp, & ! 11
472	0.32_wp, 0.40_wp, 0.24_wp, & ! 12
473	0.19_wp, 0.27_wp, 0.10_wp, & ! 13
474	0.77_wp, 0.65_wp, 0.90_wp, & ! 14
475	0.77_wp, 0.65_wp, 0.90_wp, & ! 15
476	0.82_wp, 0.70_wp, 0.95_wp, & ! 16
477	0.08_wp, 0.08_wp, 0.08_wp, & ! 17
478	0.17_wp, 0.17_wp, 0.17_wp, & ! 18
479	0.17_wp, 0.17_wp, 0.17_wp, & ! 19
480	0.30_wp, 0.30_wp, 0.30_wp, & ! 20
481	0.17_wp, 0.17_wp, 0.17_wp, & ! 21
482	0.17_wp, 0.17_wp, 0.17_wp, & ! 22
483	0.17_wp, 0.17_wp, 0.17_wp, & ! 23
484	0.17_wp, 0.17_wp, 0.17_wp, & ! 24
485	0.17_wp, 0.17_wp, 0.17_wp, & ! 25
486	0.17_wp, 0.17_wp, 0.17_wp, & ! 26
487	0.17_wp, 0.17_wp, 0.17_wp, & ! 27
488	0.17_wp, 0.17_wp, 0.17_wp, & ! 28
489	0.17_wp, 0.17_wp, 0.17_wp, & ! 29
490	0.17_wp, 0.17_wp, 0.17_wp, & ! 30
491	0.17_wp, 0.17_wp, 0.17_wp, & ! 31
492	0.17_wp, 0.17_wp, 0.17_wp, & ! 32
493	0.17_wp, 0.17_wp, 0.17_wp & ! 33
494	/), (/ 3, 33 /) )
495
496	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
497	rad_lw_cs_hr, & !< longwave clear sky radiation heating rate (K/s)
498	rad_lw_cs_hr_av, & !< average of rad_lw_cs_hr
499	rad_lw_hr, & !< longwave radiation heating rate (K/s)
500	rad_lw_hr_av, & !< average of rad_sw_hr
501	rad_lw_in, & !< incoming longwave radiation (W/m2)
502	rad_lw_in_av, & !< average of rad_lw_in
503	rad_lw_out, & !< outgoing longwave radiation (W/m2)
504	rad_lw_out_av, & !< average of rad_lw_out
505	rad_sw_cs_hr, & !< shortwave clear sky radiation heating rate (K/s)
506	rad_sw_cs_hr_av, & !< average of rad_sw_cs_hr
507	rad_sw_hr, & !< shortwave radiation heating rate (K/s)
508	rad_sw_hr_av, & !< average of rad_sw_hr
509	rad_sw_in, & !< incoming shortwave radiation (W/m2)
510	rad_sw_in_av, & !< average of rad_sw_in
511	rad_sw_out, & !< outgoing shortwave radiation (W/m2)
512	rad_sw_out_av !< average of rad_sw_out
513
514
515	!
516	!-- Variables and parameters used in RRTMG only
517	#if defined ( __rrtmg )
518	CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
519
520
521	!
522	!-- Flag parameters to be passed to RRTMG (should not be changed until ice phase in clouds is allowed)
523	INTEGER(iwp), PARAMETER :: rrtm_idrv = 1, & !< flag for longwave upward flux calculation option (0,1)
524	rrtm_inflglw = 2, & !< flag for lw cloud optical properties (0,1,2)
525	rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
526	rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
527	rrtm_inflgsw = 2, & !< flag for sw cloud optical properties (0,1,2)
528	rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
529	rrtm_liqflgsw = 1 !< flag for sw liquid droplet specifications
530
531	!
532	!-- The following variables should be only changed with care, as this will
533	!-- require further setting of some variables, which is currently not
534	!-- implemented (aerosols, ice phase).
535	INTEGER(iwp) :: nzt_rad, & !< upper vertical limit for radiation calculations
536	rrtm_icld = 0, & !< cloud flag (0: clear sky column, 1: cloudy column)
537	rrtm_iaer = 0 !< aerosol option flag (0: no aerosol layers, for lw only: 6 (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented)
538
539	INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
540
541	LOGICAL :: snd_exists = .FALSE. !< flag parameter to check whether a user-defined input files exists
542	LOGICAL :: sw_exists = .FALSE. !< flag parameter to check whether that required rrtmg sw file exists
543	LOGICAL :: lw_exists = .FALSE. !< flag parameter to check whether that required rrtmg lw file exists
544
545
546	REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
547
548	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd, & !< hypostatic pressure from sounding data (hPa)
549	rrtm_tsfc, & !< dummy array for storing surface temperature
550	t_snd !< actual temperature from sounding data (hPa)
551
552	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rrtm_ccl4vmr, & !< CCL4 volume mixing ratio (g/mol)
553	rrtm_cfc11vmr, & !< CFC11 volume mixing ratio (g/mol)
554	rrtm_cfc12vmr, & !< CFC12 volume mixing ratio (g/mol)
555	rrtm_cfc22vmr, & !< CFC22 volume mixing ratio (g/mol)
556	rrtm_ch4vmr, & !< CH4 volume mixing ratio
557	rrtm_cicewp, & !< in-cloud ice water path (g/m2)
558	rrtm_cldfr, & !< cloud fraction (0,1)
559	rrtm_cliqwp, & !< in-cloud liquid water path (g/m2)
560	rrtm_co2vmr, & !< CO2 volume mixing ratio (g/mol)
561	rrtm_emis, & !< surface emissivity (0-1)
562	rrtm_h2ovmr, & !< H2O volume mixing ratio
563	rrtm_n2ovmr, & !< N2O volume mixing ratio
564	rrtm_o2vmr, & !< O2 volume mixing ratio
565	rrtm_o3vmr, & !< O3 volume mixing ratio
566	rrtm_play, & !< pressure layers (hPa, zu-grid)
567	rrtm_plev, & !< pressure layers (hPa, zw-grid)
568	rrtm_reice, & !< cloud ice effective radius (microns)
569	rrtm_reliq, & !< cloud water drop effective radius (microns)
570	rrtm_tlay, & !< actual temperature (K, zu-grid)
571	rrtm_tlev, & !< actual temperature (K, zw-grid)
572	rrtm_lwdflx, & !< RRTM output of incoming longwave radiation flux (W/m2)
573	rrtm_lwdflxc, & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
574	rrtm_lwuflx, & !< RRTM output of outgoing longwave radiation flux (W/m2)
575	rrtm_lwuflxc, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
576	rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
577	rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
578	rrtm_lwhr, & !< RRTM output of longwave radiation heating rate (K/d)
579	rrtm_lwhrc, & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
580	rrtm_swdflx, & !< RRTM output of incoming shortwave radiation flux (W/m2)
581	rrtm_swdflxc, & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
582	rrtm_swuflx, & !< RRTM output of outgoing shortwave radiation flux (W/m2)
583	rrtm_swuflxc, & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
584	rrtm_swhr, & !< RRTM output of shortwave radiation heating rate (K/d)
585	rrtm_swhrc, & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
586	rrtm_dirdflux, & !< RRTM output of incoming direct shortwave (W/m2)
587	rrtm_difdflux !< RRTM output of incoming diffuse shortwave (W/m2)
588
589	REAL(wp), DIMENSION(1) :: rrtm_aldif, & !< surface albedo for longwave diffuse radiation
590	rrtm_aldir, & !< surface albedo for longwave direct radiation
591	rrtm_asdif, & !< surface albedo for shortwave diffuse radiation
592	rrtm_asdir !< surface albedo for shortwave direct radiation
593
594	!
595	!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
596	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rad_lw_cs_in, & !< incoming clear sky longwave radiation (W/m2) (not used)
597	rad_lw_cs_out, & !< outgoing clear sky longwave radiation (W/m2) (not used)
598	rad_sw_cs_in, & !< incoming clear sky shortwave radiation (W/m2) (not used)
599	rad_sw_cs_out, & !< outgoing clear sky shortwave radiation (W/m2) (not used)
600	rrtm_lw_tauaer, & !< lw aerosol optical depth
601	rrtm_lw_taucld, & !< lw in-cloud optical depth
602	rrtm_sw_taucld, & !< sw in-cloud optical depth
603	rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
604	rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
605	rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
606	rrtm_sw_tauaer, & !< sw aerosol optical depth
607	rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
608	rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
609	rrtm_sw_ecaer !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
610
611	#endif
612	!
613	!-- Parameters of urban and land surface models
614	INTEGER(iwp) :: nz_urban !< number of layers of urban surface (will be calculated)
615	INTEGER(iwp) :: nz_plant !< number of layers of plant canopy (will be calculated)
616	INTEGER(iwp) :: nz_urban_b !< bottom layer of urban surface (will be calculated)
617	INTEGER(iwp) :: nz_urban_t !< top layer of urban surface (will be calculated)
618	INTEGER(iwp) :: nz_plant_t !< top layer of plant canopy (will be calculated)
619	!-- parameters of urban and land surface models
620	INTEGER(iwp), PARAMETER :: nzut_free = 3 !< number of free layers above top of of topography
621	INTEGER(iwp), PARAMETER :: ndsvf = 2 !< number of dimensions of real values in SVF
622	INTEGER(iwp), PARAMETER :: idsvf = 2 !< number of dimensions of integer values in SVF
623	INTEGER(iwp), PARAMETER :: ndcsf = 1 !< number of dimensions of real values in CSF
624	INTEGER(iwp), PARAMETER :: idcsf = 2 !< number of dimensions of integer values in CSF
625	INTEGER(iwp), PARAMETER :: kdcsf = 4 !< number of dimensions of integer values in CSF calculation array
626	INTEGER(iwp), PARAMETER :: id = 1 !< position of d-index in surfl and surf
627	INTEGER(iwp), PARAMETER :: iz = 2 !< position of k-index in surfl and surf
628	INTEGER(iwp), PARAMETER :: iy = 3 !< position of j-index in surfl and surf
629	INTEGER(iwp), PARAMETER :: ix = 4 !< position of i-index in surfl and surf
630	INTEGER(iwp), PARAMETER :: im = 5 !< position of surface m-index in surfl and surf
631	INTEGER(iwp), PARAMETER :: nidx_surf = 5 !< number of indices in surfl and surf
632
633	INTEGER(iwp), PARAMETER :: nsurf_type = 10 !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
634
635	INTEGER(iwp), PARAMETER :: iup_u = 0 !< 0 - index of urban upward surface (ground or roof)
636	INTEGER(iwp), PARAMETER :: idown_u = 1 !< 1 - index of urban downward surface (overhanging)
637	INTEGER(iwp), PARAMETER :: inorth_u = 2 !< 2 - index of urban northward facing wall
638	INTEGER(iwp), PARAMETER :: isouth_u = 3 !< 3 - index of urban southward facing wall
639	INTEGER(iwp), PARAMETER :: ieast_u = 4 !< 4 - index of urban eastward facing wall
640	INTEGER(iwp), PARAMETER :: iwest_u = 5 !< 5 - index of urban westward facing wall
641
642	INTEGER(iwp), PARAMETER :: iup_l = 6 !< 6 - index of land upward surface (ground or roof)
643	INTEGER(iwp), PARAMETER :: inorth_l = 7 !< 7 - index of land northward facing wall
644	INTEGER(iwp), PARAMETER :: isouth_l = 8 !< 8 - index of land southward facing wall
645	INTEGER(iwp), PARAMETER :: ieast_l = 9 !< 9 - index of land eastward facing wall
646	INTEGER(iwp), PARAMETER :: iwest_l = 10 !< 10- index of land westward facing wall
647
648	INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1/) !< surface normal direction x indices
649	INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/) !< surface normal direction y indices
650	INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/) !< surface normal direction z indices
651	REAL(wp), DIMENSION(0:nsurf_type) :: facearea !< area of single face in respective
652	!< direction (will be calc'd)
653
654
655	!-- indices and sizes of urban and land surface models
656	INTEGER(iwp) :: startland !< start index of block of land and roof surfaces
657	INTEGER(iwp) :: endland !< end index of block of land and roof surfaces
658	INTEGER(iwp) :: nlands !< number of land and roof surfaces in local processor
659	INTEGER(iwp) :: startwall !< start index of block of wall surfaces
660	INTEGER(iwp) :: endwall !< end index of block of wall surfaces
661	INTEGER(iwp) :: nwalls !< number of wall surfaces in local processor
662
663	!-- indices needed for RTM netcdf output subroutines
664	INTEGER(iwp), PARAMETER :: nd = 5
665	CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
666	INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint_u = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
667	INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint_l = (/ iup_l, isouth_l, inorth_l, iwest_l, ieast_l /)
668	INTEGER(iwp), DIMENSION(0:nd-1) :: dirstart
669	INTEGER(iwp), DIMENSION(0:nd-1) :: dirend
670
671	!-- indices and sizes of urban and land surface models
672	INTEGER(iwp), DIMENSION(:,:), POINTER :: surfl !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x, m]
673	INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surfl_linear !< dtto (linearly allocated array)
674	INTEGER(iwp), DIMENSION(:,:), POINTER :: surf !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x, m]
675	INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surf_linear !< dtto (linearly allocated array)
676	INTEGER(iwp) :: nsurfl !< number of all surfaces in local processor
677	INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: nsurfs !< array of number of all surfaces in individual processors
678	INTEGER(iwp) :: nsurf !< global number of surfaces in index array of surfaces (nsurf = proc nsurfs)
679	INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surfstart !< starts of blocks of surfaces for individual processors in array surf (indexed from 1)
680	!< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1]
681
682	!-- block variables needed for calculation of the plant canopy model inside the urban surface model
683	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pct !< top layer of the plant canopy
684	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch !< heights of the plant canopy
685	INTEGER(iwp) :: npcbl = 0 !< number of the plant canopy gridboxes in local processor
686	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pcbl !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i]
687	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinsw !< array of absorbed sw radiation for local plant canopy box
688	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinswdir !< array of absorbed direct sw radiation for local plant canopy box
689	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinswdif !< array of absorbed diffusion sw radiation for local plant canopy box
690	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinlw !< array of absorbed lw radiation for local plant canopy box
691
692	!-- configuration parameters (they can be setup in PALM config)
693	LOGICAL :: raytrace_mpi_rma = .TRUE. !< use MPI RMA to access LAD and gridsurf from remote processes during raytracing
694	LOGICAL :: rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for
695	!< reflected radiation (as opposed to all mutually visible pairs)
696	LOGICAL :: plant_lw_interact = .TRUE. !< whether plant canopy interacts with LW radiation (in addition to SW)
697	INTEGER(iwp) :: mrt_nlevels = 0 !< number of vertical boxes above surface for which to calculate MRT
698	LOGICAL :: mrt_skip_roof = .TRUE. !< do not calculate MRT above roof surfaces
699	LOGICAL :: mrt_include_sw = .TRUE. !< should MRT calculation include SW radiation as well?
700	INTEGER(wp) :: mrt_geom = 1 !< method for MRT direction weights simulating a sphere or a human body
701	REAL(wp), DIMENSION(2) :: mrt_geom_params = (/ .12_wp, .88_wp /) !< parameters for the selected method
702	INTEGER(iwp) :: nrefsteps = 3 !< number of reflection steps to perform
703	REAL(wp), PARAMETER :: ext_coef = 0.6_wp !< extinction coefficient (a.k.a. alpha)
704	INTEGER(iwp), PARAMETER :: rad_version_len = 10 !< length of identification string of rad version
705	CHARACTER(rad_version_len), PARAMETER :: rad_version = 'RAD v. 3.0' !< identification of version of binary svf and restart files
706	INTEGER(iwp) :: raytrace_discrete_elevs = 40 !< number of discretization steps for elevation (nadir to zenith)
707	INTEGER(iwp) :: raytrace_discrete_azims = 80 !< number of discretization steps for azimuth (out of 360 degrees)
708	REAL(wp) :: max_raytracing_dist = -999.0_wp !< maximum distance for raytracing (in metres)
709	REAL(wp) :: min_irrf_value = 1e-6_wp !< minimum potential irradiance factor value for raytracing
710	REAL(wp), DIMENSION(1:30) :: svfnorm_report_thresh = 1e21_wp !< thresholds of SVF normalization values to report
711	INTEGER(iwp) :: svfnorm_report_num !< number of SVF normalization thresholds to report
712
713	!-- radiation related arrays to be used in radiation_interaction routine
714	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_dir !< direct sw radiation
715	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_diff !< diffusion sw radiation
716	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_lw_in_diff !< diffusion lw radiation
717
718	!-- parameters required for RRTMG lower boundary condition
719	REAL(wp) :: albedo_urb !< albedo value retuned to RRTMG boundary cond.
720	REAL(wp) :: emissivity_urb !< emissivity value retuned to RRTMG boundary cond.
721	REAL(wp) :: t_rad_urb !< temperature value retuned to RRTMG boundary cond.
722
723	!-- type for calculation of svf
724	TYPE t_svf
725	INTEGER(iwp) :: isurflt !<
726	INTEGER(iwp) :: isurfs !<
727	REAL(wp) :: rsvf !<
728	REAL(wp) :: rtransp !<
729	END TYPE
730
731	!-- type for calculation of csf
732	TYPE t_csf
733	INTEGER(iwp) :: ip !<
734	INTEGER(iwp) :: itx !<
735	INTEGER(iwp) :: ity !<
736	INTEGER(iwp) :: itz !<
737	INTEGER(iwp) :: isurfs !< Idx of source face / -1 for sky
738	REAL(wp) :: rcvf !< Canopy view factor for faces /
739	!< canopy sink factor for sky (-1)
740	END TYPE
741
742	!-- arrays storing the values of USM
743	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: svfsurf !< svfsurf[:,isvf] = index of target and source surface for svf[isvf]
744	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: svf !< array of shape view factors+direct irradiation factors for local surfaces
745	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins !< array of sw radiation falling to local surface after i-th reflection
746	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl !< array of lw radiation for local surface after i-th reflection
747
748	REAL(wp), DIMENSION(:), ALLOCATABLE :: skyvf !< array of sky view factor for each local surface
749	REAL(wp), DIMENSION(:), ALLOCATABLE :: skyvft !< array of sky view factor including transparency for each local surface
750	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dsitrans !< dsidir[isvfl,i] = path transmittance of i-th
751	!< direction of direct solar irradiance per target surface
752	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dsitransc !< dtto per plant canopy box
753	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dsidir !< dsidir[:,i] = unit vector of i-th
754	!< direction of direct solar irradiance
755	INTEGER(iwp) :: ndsidir !< number of apparent solar directions used
756	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: dsidir_rev !< dsidir_rev[ielev,iazim] = i for dsidir or -1 if not present
757
758	INTEGER(iwp) :: nmrtbl !< No. of local grid boxes for which MRT is calculated
759	INTEGER(iwp) :: nmrtf !< number of MRT factors for local processor
760	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: mrtbl !< coordinates of i-th local MRT box - surfl[:,i] = [z, y, x]
761	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: mrtfsurf !< mrtfsurf[:,imrtf] = index of target MRT box and source surface for mrtf[imrtf]
762	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtf !< array of MRT factors for each local MRT box
763	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtft !< array of MRT factors including transparency for each local MRT box
764	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtsky !< array of sky view factor for each local MRT box
765	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtskyt !< array of sky view factor including transparency for each local MRT box
766	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mrtdsit !< array of direct solar transparencies for each local MRT box
767	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtinsw !< mean SW radiant flux for each MRT box
768	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtinlw !< mean LW radiant flux for each MRT box
769	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrt !< mean radiant temperature for each MRT box
770	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtinsw_av !< time average mean SW radiant flux for each MRT box
771	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtinlw_av !< time average mean LW radiant flux for each MRT box
772	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrt_av !< time average mean radiant temperature for each MRT box
773
774	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw !< array of sw radiation falling to local surface including radiation from reflections
775	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw !< array of lw radiation falling to local surface including radiation from reflections
776	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdir !< array of direct sw radiation falling to local surface
777	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif !< array of diffuse sw radiation from sky and model boundary falling to local surface
778	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif !< array of diffuse lw radiation from sky and model boundary falling to local surface
779
780	!< Outward radiation is only valid for nonvirtual surfaces
781	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsl !< array of reflected sw radiation for local surface in i-th reflection
782	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutll !< array of reflected + emitted lw radiation for local surface in i-th reflection
783	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfouts !< array of reflected sw radiation for all surfaces in i-th reflection
784	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutl !< array of reflected + emitted lw radiation for all surfaces in i-th reflection
785	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlg !< global array of incoming lw radiation from plant canopy
786	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw !< array of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
787	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
788	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfemitlwl !< array of emitted lw radiation for local surface used to calculate effective surface temperature for radiation model
789
790	!-- block variables needed for calculation of the plant canopy model inside the urban surface model
791	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: csfsurf !< csfsurf[:,icsf] = index of target surface and csf grid index for csf[icsf]
792	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: csf !< array of plant canopy sink fators + direct irradiation factors (transparency)
793	REAL(wp), DIMENSION(:,:,:), POINTER :: sub_lad !< subset of lad_s within urban surface, transformed to plain Z coordinate
794	REAL(wp), DIMENSION(:), POINTER :: sub_lad_g !< sub_lad globalized (used to avoid MPI RMA calls in raytracing)
795	REAL(wp) :: prototype_lad !< prototype leaf area density for computing effective optical depth
796	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nzterr, plantt !< temporary global arrays for raytracing
797	INTEGER(iwp) :: plantt_max
798
799	!-- arrays and variables for calculation of svf and csf
800	TYPE(t_svf), DIMENSION(:), POINTER :: asvf !< pointer to growing svc array
801	TYPE(t_csf), DIMENSION(:), POINTER :: acsf !< pointer to growing csf array
802	TYPE(t_svf), DIMENSION(:), POINTER :: amrtf !< pointer to growing mrtf array
803	TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET :: asvf1, asvf2 !< realizations of svf array
804	TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET :: acsf1, acsf2 !< realizations of csf array
805	TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET :: amrtf1, amrtf2 !< realizations of mftf array
806	INTEGER(iwp) :: nsvfla !< dimmension of array allocated for storage of svf in local processor
807	INTEGER(iwp) :: ncsfla !< dimmension of array allocated for storage of csf in local processor
808	INTEGER(iwp) :: nmrtfa !< dimmension of array allocated for storage of mrt
809	INTEGER(iwp) :: msvf, mcsf, mmrtf!< mod for swapping the growing array
810	INTEGER(iwp), PARAMETER :: gasize = 100000_iwp !< initial size of growing arrays
811	REAL(wp), PARAMETER :: grow_factor = 1.4_wp !< growth factor of growing arrays
812	INTEGER(iwp) :: nsvfl !< number of svf for local processor
813	INTEGER(iwp) :: ncsfl !< no. of csf in local processor
814	!< needed only during calc_svf but must be here because it is
815	!< shared between subroutines calc_svf and raytrace
816	INTEGER(iwp), DIMENSION(:,:,:,:), POINTER :: gridsurf !< reverse index of local surfl[d,k,j,i] (for case rad_angular_discretization)
817	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: gridpcbl !< reverse index of local pcbl[k,j,i]
818	INTEGER(iwp), PARAMETER :: nsurf_type_u = 6 !< number of urban surf types (used in gridsurf)
819
820	!-- temporary arrays for calculation of csf in raytracing
821	INTEGER(iwp) :: maxboxesg !< max number of boxes ray can cross in the domain
822	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: boxes !< coordinates of gridboxes being crossed by ray
823	REAL(wp), DIMENSION(:), ALLOCATABLE :: crlens !< array of crossing lengths of ray for particular grid boxes
824	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: lad_ip !< array of numbers of process where lad is stored
825	#if defined( __parallel )
826	INTEGER(kind=MPI_ADDRESS_KIND), &
827	DIMENSION(:), ALLOCATABLE :: lad_disp !< array of displaycements of lad in local array of proc lad_ip
828	INTEGER(iwp) :: win_lad !< MPI RMA window for leaf area density
829	INTEGER(iwp) :: win_gridsurf !< MPI RMA window for reverse grid surface index
830	#endif
831	REAL(wp), DIMENSION(:), ALLOCATABLE :: lad_s_ray !< array of received lad_s for appropriate gridboxes crossed by ray
832	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: target_surfl
833	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: rt2_track
834	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rt2_track_lad
835	REAL(wp), DIMENSION(:), ALLOCATABLE :: rt2_track_dist
836	REAL(wp), DIMENSION(:), ALLOCATABLE :: rt2_dist
837
838	!-- arrays for time averages
839	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfradnet_av !< average of net radiation to local surface including radiation from reflections
840	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw_av !< average of sw radiation falling to local surface including radiation from reflections
841	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw_av !< average of lw radiation falling to local surface including radiation from reflections
842	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdir_av !< average of direct sw radiation falling to local surface
843	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif_av !< average of diffuse sw radiation from sky and model boundary falling to local surface
844	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif_av !< average of diffuse lw radiation from sky and model boundary falling to local surface
845	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswref_av !< average of sw radiation falling to surface from reflections
846	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwref_av !< average of lw radiation falling to surface from reflections
847	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw_av !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
848	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw_av !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
849	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins_av !< average of array of residua of sw radiation absorbed in surface after last reflection
850	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl_av !< average of array of residua of lw radiation absorbed in surface after last reflection
851	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinlw_av !< Average of pcbinlw
852	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinsw_av !< Average of pcbinsw
853	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinswdir_av !< Average of pcbinswdir
854	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinswdif_av !< Average of pcbinswdif
855	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinswref_av !< Average of pcbinswref
856
857
858	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
859	!-- Energy balance variables
860	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
861	!-- parameters of the land, roof and wall surfaces
862	REAL(wp), DIMENSION(:), ALLOCATABLE :: albedo_surf !< albedo of the surface
863	REAL(wp), DIMENSION(:), ALLOCATABLE :: emiss_surf !< emissivity of the wall surface
864	!
865	!-- External radiation. Depending on the given level of detail either a 1D or
866	!-- a 3D array will be allocated.
867	TYPE( real_1d_3d ) :: rad_lw_in_f !< external incoming longwave radiation, from observation or model
868	TYPE( real_1d_3d ) :: rad_sw_in_f !< external incoming shortwave radiation, from observation or model
869	TYPE( real_1d_3d ) :: rad_sw_in_dif_f !< external incoming shortwave radiation, diffuse part, from observation or model
870	TYPE( real_1d_3d ) :: time_rad_f !< time dimension for external radiation, from observation or model
871
872	INTERFACE radiation_check_data_output
873	MODULE PROCEDURE radiation_check_data_output
874	END INTERFACE radiation_check_data_output
875
876	INTERFACE radiation_check_data_output_ts
877	MODULE PROCEDURE radiation_check_data_output_ts
878	END INTERFACE radiation_check_data_output_ts
879
880	INTERFACE radiation_check_data_output_pr
881	MODULE PROCEDURE radiation_check_data_output_pr
882	END INTERFACE radiation_check_data_output_pr
883
884	INTERFACE radiation_check_parameters
885	MODULE PROCEDURE radiation_check_parameters
886	END INTERFACE radiation_check_parameters
887
888	INTERFACE radiation_clearsky
889	MODULE PROCEDURE radiation_clearsky
890	END INTERFACE radiation_clearsky
891
892	INTERFACE radiation_constant
893	MODULE PROCEDURE radiation_constant
894	END INTERFACE radiation_constant
895
896	INTERFACE radiation_control
897	MODULE PROCEDURE radiation_control
898	END INTERFACE radiation_control
899
900	INTERFACE radiation_3d_data_averaging
901	MODULE PROCEDURE radiation_3d_data_averaging
902	END INTERFACE radiation_3d_data_averaging
903
904	INTERFACE radiation_data_output_2d
905	MODULE PROCEDURE radiation_data_output_2d
906	END INTERFACE radiation_data_output_2d
907
908	INTERFACE radiation_data_output_3d
909	MODULE PROCEDURE radiation_data_output_3d
910	END INTERFACE radiation_data_output_3d
911
912	INTERFACE radiation_data_output_mask
913	MODULE PROCEDURE radiation_data_output_mask
914	END INTERFACE radiation_data_output_mask
915
916	INTERFACE radiation_define_netcdf_grid
917	MODULE PROCEDURE radiation_define_netcdf_grid
918	END INTERFACE radiation_define_netcdf_grid
919
920	INTERFACE radiation_header
921	MODULE PROCEDURE radiation_header
922	END INTERFACE radiation_header
923
924	INTERFACE radiation_init
925	MODULE PROCEDURE radiation_init
926	END INTERFACE radiation_init
927
928	INTERFACE radiation_parin
929	MODULE PROCEDURE radiation_parin
930	END INTERFACE radiation_parin
931
932	INTERFACE radiation_rrtmg
933	MODULE PROCEDURE radiation_rrtmg
934	END INTERFACE radiation_rrtmg
935
936	#if defined( __rrtmg )
937	INTERFACE radiation_tendency
938	MODULE PROCEDURE radiation_tendency
939	MODULE PROCEDURE radiation_tendency_ij
940	END INTERFACE radiation_tendency
941	#endif
942
943	INTERFACE radiation_rrd_local
944	MODULE PROCEDURE radiation_rrd_local
945	END INTERFACE radiation_rrd_local
946
947	INTERFACE radiation_wrd_local
948	MODULE PROCEDURE radiation_wrd_local
949	END INTERFACE radiation_wrd_local
950
951	INTERFACE radiation_interaction
952	MODULE PROCEDURE radiation_interaction
953	END INTERFACE radiation_interaction
954
955	INTERFACE radiation_interaction_init
956	MODULE PROCEDURE radiation_interaction_init
957	END INTERFACE radiation_interaction_init
958
959	INTERFACE radiation_presimulate_solar_pos
960	MODULE PROCEDURE radiation_presimulate_solar_pos
961	END INTERFACE radiation_presimulate_solar_pos
962
963	INTERFACE radiation_calc_svf
964	MODULE PROCEDURE radiation_calc_svf
965	END INTERFACE radiation_calc_svf
966
967	INTERFACE radiation_write_svf
968	MODULE PROCEDURE radiation_write_svf
969	END INTERFACE radiation_write_svf
970
971	INTERFACE radiation_read_svf
972	MODULE PROCEDURE radiation_read_svf
973	END INTERFACE radiation_read_svf
974
975
976	SAVE
977
978	PRIVATE
979
980	!
981	!-- Public functions / NEEDS SORTING
982	PUBLIC radiation_check_data_output, radiation_check_data_output_pr, &
983	radiation_check_data_output_ts, &
984	radiation_check_parameters, radiation_control, &
985	radiation_header, radiation_init, radiation_parin, &
986	radiation_3d_data_averaging, &
987	radiation_data_output_2d, radiation_data_output_3d, &
988	radiation_define_netcdf_grid, radiation_wrd_local, &
989	radiation_rrd_local, radiation_data_output_mask, &
990	radiation_calc_svf, radiation_write_svf, &
991	radiation_interaction, radiation_interaction_init, &
992	radiation_read_svf, radiation_presimulate_solar_pos
993
994
995	!
996	!-- Public variables and constants / NEEDS SORTING
997	PUBLIC albedo, albedo_type, decl_1, decl_2, decl_3, dots_rad, dt_radiation,&
998	emissivity, force_radiation_call, lat, lon, mrt_geom, &
999	mrt_geom_params, &
1000	mrt_include_sw, mrt_nlevels, mrtbl, mrtinsw, mrtinlw, nmrtbl, &
1001	rad_net_av, radiation, radiation_scheme, rad_lw_in, &
1002	rad_lw_in_av, rad_lw_out, rad_lw_out_av, &
1003	rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in, &
1004	rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr, &
1005	rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, solar_constant, &
1006	skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
1007	cos_zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon, &
1008	idir, jdir, kdir, id, iz, iy, ix, &
1009	iup_u, inorth_u, isouth_u, ieast_u, iwest_u, &
1010	iup_l, inorth_l, isouth_l, ieast_l, iwest_l, &
1011	nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf, &
1012	idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct, &
1013	radiation_interactions, startwall, startland, endland, endwall, &
1014	skyvf, skyvft, radiation_interactions_on, average_radiation, &
1015	rad_sw_in_diff, rad_sw_in_dir
1016
1017
1018	#if defined ( __rrtmg )
1019	PUBLIC radiation_tendency, rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
1020	#endif
1021
1022	CONTAINS
1023
1024
1025	!------------------------------------------------------------------------------!
1026	! Description:
1027	! ------------
1028	!> This subroutine controls the calls of the radiation schemes
1029	!------------------------------------------------------------------------------!
1030	SUBROUTINE radiation_control
1031
1032
1033	IMPLICIT NONE
1034
1035
1036	IF ( debug_output_timestep ) CALL debug_message( 'radiation_control', 'start' )
1037
1038
1039	SELECT CASE ( TRIM( radiation_scheme ) )
1040
1041	CASE ( 'constant' )
1042	CALL radiation_constant
1043
1044	CASE ( 'clear-sky' )
1045	CALL radiation_clearsky
1046
1047	CASE ( 'rrtmg' )
1048	CALL radiation_rrtmg
1049
1050	CASE ( 'external' )
1051	!
1052	!-- During spinup apply clear-sky model
1053	IF ( time_since_reference_point < 0.0_wp ) THEN
1054	CALL radiation_clearsky
1055	ELSE
1056	CALL radiation_external
1057	ENDIF
1058
1059	CASE DEFAULT
1060
1061	END SELECT
1062
1063	IF ( debug_output_timestep ) CALL debug_message( 'radiation_control', 'end' )
1064
1065	END SUBROUTINE radiation_control
1066
1067	!------------------------------------------------------------------------------!
1068	! Description:
1069	! ------------
1070	!> Check data output for radiation model
1071	!------------------------------------------------------------------------------!
1072	SUBROUTINE radiation_check_data_output( variable, unit, i, ilen, k )
1073
1074
1075	USE control_parameters, &
1076	ONLY: data_output, message_string
1077
1078	IMPLICIT NONE
1079
1080	CHARACTER (LEN=*) :: unit !<
1081	CHARACTER (LEN=*) :: variable !<
1082
1083	INTEGER(iwp) :: i, k
1084	INTEGER(iwp) :: ilen
1085	CHARACTER(LEN=varnamelength) :: var !< TRIM(variable)
1086
1087	var = TRIM(variable)
1088
1089	IF ( len(var) < 3_iwp ) THEN
1090	unit = 'illegal'
1091	RETURN
1092	ENDIF
1093
1094	IF ( var(1:3) /= 'rad' .AND. var(1:3) /= 'rtm' ) THEN
1095	unit = 'illegal'
1096	RETURN
1097	ENDIF
1098
1099	!-- first process diractional variables
1100	IF ( var(1:12) == 'rtm_rad_net_' .OR. var(1:13) == 'rtm_rad_insw_' .OR. &
1101	var(1:13) == 'rtm_rad_inlw_' .OR. var(1:16) == 'rtm_rad_inswdir_' .OR. &
1102	var(1:16) == 'rtm_rad_inswdif_' .OR. var(1:16) == 'rtm_rad_inswref_' .OR. &
1103	var(1:16) == 'rtm_rad_inlwdif_' .OR. var(1:16) == 'rtm_rad_inlwref_' .OR. &
1104	var(1:14) == 'rtm_rad_outsw_' .OR. var(1:14) == 'rtm_rad_outlw_' .OR. &
1105	var(1:14) == 'rtm_rad_ressw_' .OR. var(1:14) == 'rtm_rad_reslw_' ) THEN
1106	IF ( .NOT. radiation ) THEN
1107	message_string = 'output of "' // TRIM( var ) // '" require'&
1108	// 's radiation = .TRUE.'
1109	CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
1110	ENDIF
1111	unit = 'W/m2'
1112	ELSE IF ( var(1:7) == 'rtm_svf' .OR. var(1:7) == 'rtm_dif' .OR. &
1113	var(1:9) == 'rtm_skyvf' .OR. var(1:9) == 'rtm_skyvft' .OR. &
1114	var(1:12) == 'rtm_surfalb_' .OR. var(1:13) == 'rtm_surfemis_' ) THEN
1115	IF ( .NOT. radiation ) THEN
1116	message_string = 'output of "' // TRIM( var ) // '" require'&
1117	// 's radiation = .TRUE.'
1118	CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
1119	ENDIF
1120	unit = '1'
1121	ELSE
1122	!-- non-directional variables
1123	SELECT CASE ( TRIM( var ) )
1124	CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out', &
1125	'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out' )
1126	IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' ) THEN
1127	message_string = '"output of "' // TRIM( var ) // '" requi' // &
1128	'res radiation = .TRUE. and ' // &
1129	'radiation_scheme = "rrtmg"'
1130	CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
1131	ENDIF
1132	unit = 'K/h'
1133
1134	CASE ( 'rad_net', 'rrtm_aldif', 'rrtm_aldir', 'rrtm_asdif', &
1135	'rrtm_asdir', 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', &
1136	'rad_sw_out*')
1137	IF ( i == 0 .AND. ilen == 0 .AND. k == 0) THEN
1138	! Workaround for masked output (calls with i=ilen=k=0)
1139	unit = 'illegal'
1140	RETURN
1141	ENDIF
1142	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
1143	message_string = 'illegal value for data_output: "' // &
1144	TRIM( var ) // '" & only 2d-horizontal ' // &
1145	'cross sections are allowed for this value'
1146	CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
1147	ENDIF
1148	IF ( .NOT. radiation .OR. radiation_scheme /= "rrtmg" ) THEN
1149	IF ( TRIM( var ) == 'rrtm_aldif*' .OR. &
1150	TRIM( var ) == 'rrtm_aldir*' .OR. &
1151	TRIM( var ) == 'rrtm_asdif*' .OR. &
1152	TRIM( var ) == 'rrtm_asdir*' ) &
1153	THEN
1154	message_string = 'output of "' // TRIM( var ) // '" require'&
1155	// 's radiation = .TRUE. and radiation_sch'&
1156	// 'eme = "rrtmg"'
1157	CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
1158	ENDIF
1159	ENDIF
1160
1161	IF ( TRIM( var ) == 'rad_net*' ) unit = 'W/m2'
1162	IF ( TRIM( var ) == 'rad_lw_in*' ) unit = 'W/m2'
1163	IF ( TRIM( var ) == 'rad_lw_out*' ) unit = 'W/m2'
1164	IF ( TRIM( var ) == 'rad_sw_in*' ) unit = 'W/m2'
1165	IF ( TRIM( var ) == 'rad_sw_out*' ) unit = 'W/m2'
1166	IF ( TRIM( var ) == 'rad_sw_in' ) unit = 'W/m2'
1167	IF ( TRIM( var ) == 'rrtm_aldif*' ) unit = ''
1168	IF ( TRIM( var ) == 'rrtm_aldir*' ) unit = ''
1169	IF ( TRIM( var ) == 'rrtm_asdif*' ) unit = ''
1170	IF ( TRIM( var ) == 'rrtm_asdir*' ) unit = ''
1171
1172	CASE ( 'rtm_rad_pc_inlw', 'rtm_rad_pc_insw', 'rtm_rad_pc_inswdir', &
1173	'rtm_rad_pc_inswdif', 'rtm_rad_pc_inswref')
1174	IF ( .NOT. radiation ) THEN
1175	message_string = 'output of "' // TRIM( var ) // '" require'&
1176	// 's radiation = .TRUE.'
1177	CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
1178	ENDIF
1179	unit = 'W'
1180
1181	CASE ( 'rtm_mrt', 'rtm_mrt_sw', 'rtm_mrt_lw' )
1182	IF ( i == 0 .AND. ilen == 0 .AND. k == 0) THEN
1183	! Workaround for masked output (calls with i=ilen=k=0)
1184	unit = 'illegal'
1185	RETURN
1186	ENDIF
1187
1188	IF ( .NOT. radiation ) THEN
1189	message_string = 'output of "' // TRIM( var ) // '" require'&
1190	// 's radiation = .TRUE.'
1191	CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
1192	ENDIF
1193	IF ( mrt_nlevels == 0 ) THEN
1194	message_string = 'output of "' // TRIM( var ) // '" require'&
1195	// 's mrt_nlevels > 0'
1196	CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
1197	ENDIF
1198	IF ( TRIM( var ) == 'rtm_mrt_sw' .AND. .NOT. mrt_include_sw ) THEN
1199	message_string = 'output of "' // TRIM( var ) // '" require'&
1200	// 's rtm_mrt_sw = .TRUE.'
1201	CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
1202	ENDIF
1203	IF ( TRIM( var ) == 'rtm_mrt' ) THEN
1204	unit = 'K'
1205	ELSE
1206	unit = 'W m-2'
1207	ENDIF
1208
1209	CASE DEFAULT
1210	unit = 'illegal'
1211
1212	END SELECT
1213	ENDIF
1214
1215	END SUBROUTINE radiation_check_data_output
1216
1217
1218	!------------------------------------------------------------------------------!
1219	! Description:
1220	! ------------
1221	!> Set module-specific timeseries units and labels
1222	!------------------------------------------------------------------------------!
1223	SUBROUTINE radiation_check_data_output_ts( dots_max, dots_num )
1224
1225
1226	INTEGER(iwp), INTENT(IN) :: dots_max
1227	INTEGER(iwp), INTENT(INOUT) :: dots_num
1228
1229	!
1230	!-- Next line is just to avoid compiler warning about unused variable.
1231	IF ( dots_max == 0 ) CONTINUE
1232
1233	!
1234	!-- Temporary solution to add LSM and radiation time series to the default
1235	!-- output
1236	IF ( land_surface .OR. radiation ) THEN
1237	IF ( TRIM( radiation_scheme ) == 'rrtmg' ) THEN
1238	dots_num = dots_num + 15
1239	ELSE
1240	dots_num = dots_num + 11
1241	ENDIF
1242	ENDIF
1243
1244
1245	END SUBROUTINE radiation_check_data_output_ts
1246
1247	!------------------------------------------------------------------------------!
1248	! Description:
1249	! ------------
1250	!> Check data output of profiles for radiation model
1251	!------------------------------------------------------------------------------!
1252	SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, &
1253	dopr_unit )
1254
1255	USE arrays_3d, &
1256	ONLY: zu
1257
1258	USE control_parameters, &
1259	ONLY: data_output_pr, message_string
1260
1261	USE indices
1262
1263	USE profil_parameter
1264
1265	USE statistics
1266
1267	IMPLICIT NONE
1268
1269	CHARACTER (LEN=*) :: unit !<
1270	CHARACTER (LEN=*) :: variable !<
1271	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
1272
1273	INTEGER(iwp) :: var_count !<
1274
1275	SELECT CASE ( TRIM( variable ) )
1276
1277	CASE ( 'rad_net' )
1278	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )&
1279	THEN
1280	message_string = 'data_output_pr = ' // &
1281	TRIM( data_output_pr(var_count) ) // ' is' // &
1282	'not available for radiation = .FALSE. or ' //&
1283	'radiation_scheme = "constant"'
1284	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1285	ELSE
1286	dopr_index(var_count) = 99
1287	dopr_unit = 'W/m2'
1288	hom(:,2,99,:) = SPREAD( zw, 2, statistic_regions+1 )
1289	unit = dopr_unit
1290	ENDIF
1291
1292	CASE ( 'rad_lw_in' )
1293	IF ( ( .NOT. radiation) .OR. radiation_scheme == 'constant' ) &
1294	THEN
1295	message_string = 'data_output_pr = ' // &
1296	TRIM( data_output_pr(var_count) ) // ' is' // &
1297	'not available for radiation = .FALSE. or ' //&
1298	'radiation_scheme = "constant"'
1299	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1300	ELSE
1301	dopr_index(var_count) = 100
1302	dopr_unit = 'W/m2'
1303	hom(:,2,100,:) = SPREAD( zw, 2, statistic_regions+1 )
1304	unit = dopr_unit
1305	ENDIF
1306
1307	CASE ( 'rad_lw_out' )
1308	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) &
1309	THEN
1310	message_string = 'data_output_pr = ' // &
1311	TRIM( data_output_pr(var_count) ) // ' is' // &
1312	'not available for radiation = .FALSE. or ' //&
1313	'radiation_scheme = "constant"'
1314	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1315	ELSE
1316	dopr_index(var_count) = 101
1317	dopr_unit = 'W/m2'
1318	hom(:,2,101,:) = SPREAD( zw, 2, statistic_regions+1 )
1319	unit = dopr_unit
1320	ENDIF
1321
1322	CASE ( 'rad_sw_in' )
1323	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) &
1324	THEN
1325	message_string = 'data_output_pr = ' // &
1326	TRIM( data_output_pr(var_count) ) // ' is' // &
1327	'not available for radiation = .FALSE. or ' //&
1328	'radiation_scheme = "constant"'
1329	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1330	ELSE
1331	dopr_index(var_count) = 102
1332	dopr_unit = 'W/m2'
1333	hom(:,2,102,:) = SPREAD( zw, 2, statistic_regions+1 )
1334	unit = dopr_unit
1335	ENDIF
1336
1337	CASE ( 'rad_sw_out')
1338	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )&
1339	THEN
1340	message_string = 'data_output_pr = ' // &
1341	TRIM( data_output_pr(var_count) ) // ' is' // &
1342	'not available for radiation = .FALSE. or ' //&
1343	'radiation_scheme = "constant"'
1344	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1345	ELSE
1346	dopr_index(var_count) = 103
1347	dopr_unit = 'W/m2'
1348	hom(:,2,103,:) = SPREAD( zw, 2, statistic_regions+1 )
1349	unit = dopr_unit
1350	ENDIF
1351
1352	CASE ( 'rad_lw_cs_hr' )
1353	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
1354	THEN
1355	message_string = 'data_output_pr = ' // &
1356	TRIM( data_output_pr(var_count) ) // ' is' // &
1357	'not available for radiation = .FALSE. or ' //&
1358	'radiation_scheme /= "rrtmg"'
1359	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
1360	ELSE
1361	dopr_index(var_count) = 104
1362	dopr_unit = 'K/h'
1363	hom(:,2,104,:) = SPREAD( zu, 2, statistic_regions+1 )
1364	unit = dopr_unit
1365	ENDIF
1366
1367	CASE ( 'rad_lw_hr' )
1368	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
1369	THEN
1370	message_string = 'data_output_pr = ' // &
1371	TRIM( data_output_pr(var_count) ) // ' is' // &
1372	'not available for radiation = .FALSE. or ' //&
1373	'radiation_scheme /= "rrtmg"'
1374	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
1375	ELSE
1376	dopr_index(var_count) = 105
1377	dopr_unit = 'K/h'
1378	hom(:,2,105,:) = SPREAD( zu, 2, statistic_regions+1 )
1379	unit = dopr_unit
1380	ENDIF
1381
1382	CASE ( 'rad_sw_cs_hr' )
1383	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
1384	THEN
1385	message_string = 'data_output_pr = ' // &
1386	TRIM( data_output_pr(var_count) ) // ' is' // &
1387	'not available for radiation = .FALSE. or ' //&
1388	'radiation_scheme /= "rrtmg"'
1389	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
1390	ELSE
1391	dopr_index(var_count) = 106
1392	dopr_unit = 'K/h'
1393	hom(:,2,106,:) = SPREAD( zu, 2, statistic_regions+1 )
1394	unit = dopr_unit
1395	ENDIF
1396
1397	CASE ( 'rad_sw_hr' )
1398	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
1399	THEN
1400	message_string = 'data_output_pr = ' // &
1401	TRIM( data_output_pr(var_count) ) // ' is' // &
1402	'not available for radiation = .FALSE. or ' //&
1403	'radiation_scheme /= "rrtmg"'
1404	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
1405	ELSE
1406	dopr_index(var_count) = 107
1407	dopr_unit = 'K/h'
1408	hom(:,2,107,:) = SPREAD( zu, 2, statistic_regions+1 )
1409	unit = dopr_unit
1410	ENDIF
1411
1412
1413	CASE DEFAULT
1414	unit = 'illegal'
1415
1416	END SELECT
1417
1418
1419	END SUBROUTINE radiation_check_data_output_pr
1420
1421
1422	!------------------------------------------------------------------------------!
1423	! Description:
1424	! ------------
1425	!> Check parameters routine for radiation model
1426	!------------------------------------------------------------------------------!
1427	SUBROUTINE radiation_check_parameters
1428
1429	USE control_parameters, &
1430	ONLY: land_surface, message_string, rotation_angle, urban_surface
1431
1432	USE netcdf_data_input_mod, &
1433	ONLY: input_pids_static
1434
1435	IMPLICIT NONE
1436
1437	!
1438	!-- In case no urban-surface or land-surface model is applied, usage of
1439	!-- a radiation model make no sense.
1440	IF ( .NOT. land_surface .AND. .NOT. urban_surface ) THEN
1441	message_string = 'Usage of radiation module is only allowed if ' // &
1442	'land-surface and/or urban-surface model is applied.'
1443	CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
1444	ENDIF
1445
1446	IF ( radiation_scheme /= 'constant' .AND. &
1447	radiation_scheme /= 'clear-sky' .AND. &
1448	radiation_scheme /= 'rrtmg' .AND. &
1449	radiation_scheme /= 'external' ) THEN
1450	message_string = 'unknown radiation_scheme = '// &
1451	TRIM( radiation_scheme )
1452	CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
1453	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
1454	#if ! defined ( __rrtmg )
1455	message_string = 'radiation_scheme = "rrtmg" requires ' // &
1456	'compilation of PALM with pre-processor ' // &
1457	'directive -D__rrtmg'
1458	CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
1459	#endif
1460	#if defined ( __rrtmg ) && ! defined( __netcdf )
1461	message_string = 'radiation_scheme = "rrtmg" requires ' // &
1462	'the use of NetCDF (preprocessor directive ' // &
1463	'-D__netcdf'
1464	CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
1465	#endif
1466
1467	ENDIF
1468	!
1469	!-- Checks performed only if data is given via namelist only.
1470	IF ( .NOT. input_pids_static ) THEN
1471	IF ( albedo_type == 0 .AND. albedo == 9999999.9_wp .AND. &
1472	radiation_scheme == 'clear-sky') THEN
1473	message_string = 'radiation_scheme = "clear-sky" in combination'//&
1474	'with albedo_type = 0 requires setting of'// &
1475	'albedo /= 9999999.9'
1476	CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
1477	ENDIF
1478
1479	IF ( albedo_type == 0 .AND. radiation_scheme == 'rrtmg' .AND. &
1480	( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
1481	.OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp&
1482	) ) THEN
1483	message_string = 'radiation_scheme = "rrtmg" in combination' // &
1484	'with albedo_type = 0 requires setting of ' // &
1485	'albedo_lw_dif /= 9999999.9' // &
1486	'albedo_lw_dir /= 9999999.9' // &
1487	'albedo_sw_dif /= 9999999.9 and' // &
1488	'albedo_sw_dir /= 9999999.9'
1489	CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
1490	ENDIF
1491	ENDIF
1492	!
1493	!-- Parallel rad_angular_discretization without raytrace_mpi_rma is not implemented
1494	#if defined( __parallel )
1495	IF ( rad_angular_discretization .AND. .NOT. raytrace_mpi_rma ) THEN
1496	message_string = 'rad_angular_discretization can only be used ' // &
1497	'together with raytrace_mpi_rma or when ' // &
1498	'no parallelization is applied.'
1499	CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
1500	ENDIF
1501	#endif
1502
1503	IF ( cloud_droplets .AND. radiation_scheme == 'rrtmg' .AND. &
1504	average_radiation ) THEN
1505	message_string = 'average_radiation = .T. with radiation_scheme'// &
1506	'= "rrtmg" in combination cloud_droplets = .T.'// &
1507	'is not implementd'
1508	CALL message( 'check_parameters', 'PA0560', 1, 2, 0, 6, 0 )
1509	ENDIF
1510
1511	!
1512	!-- Incialize svf normalization reporting histogram
1513	svfnorm_report_num = 1
1514	DO WHILE ( svfnorm_report_thresh(svfnorm_report_num) < 1e20_wp &
1515	.AND. svfnorm_report_num <= 30 )
1516	svfnorm_report_num = svfnorm_report_num + 1
1517	ENDDO
1518	svfnorm_report_num = svfnorm_report_num - 1
1519	!
1520	!-- Check for dt_radiation
1521	IF ( dt_radiation <= 0.0 ) THEN
1522	message_string = 'dt_radiation must be > 0.0'
1523	CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 )
1524	ENDIF
1525	!
1526	!-- Check rotation angle
1527	!> @todo Remove this limitation
1528	IF ( rotation_angle /= 0.0 ) THEN
1529	message_string = 'rotation of the model domain is not considered in the radiation ' // &
1530	'model.&Using rotation_angle /= 0.0 is not allowed in combination ' // &
1531	'with the radiation model at the moment!'
1532	CALL message( 'check_parameters', 'PA0675', 1, 2, 0, 6, 0 )
1533	ENDIF
1534
1535	END SUBROUTINE radiation_check_parameters
1536
1537
1538	!------------------------------------------------------------------------------!
1539	! Description:
1540	! ------------
1541	!> Initialization of the radiation model and Radiative Transfer Model
1542	!------------------------------------------------------------------------------!
1543	SUBROUTINE radiation_init
1544
1545	IMPLICIT NONE
1546
1547	INTEGER(iwp) :: i !< running index x-direction
1548	INTEGER(iwp) :: is !< running index for input surface elements
1549	INTEGER(iwp) :: ioff !< offset in x between surface element reference grid point in atmosphere and actual surface
1550	INTEGER(iwp) :: j !< running index y-direction
1551	INTEGER(iwp) :: joff !< offset in y between surface element reference grid point in atmosphere and actual surface
1552	INTEGER(iwp) :: k !< running index z-direction
1553	INTEGER(iwp) :: l !< running index for orientation of vertical surfaces
1554	INTEGER(iwp) :: m !< running index for surface elements
1555	INTEGER(iwp) :: ntime = 0 !< number of available external radiation timesteps
1556	#if defined( __rrtmg )
1557	INTEGER(iwp) :: ind_type !< running index for subgrid-surface tiles
1558	#endif
1559	LOGICAL :: radiation_input_root_domain !< flag indicating the existence of a dynamic input file for the root domain
1560
1561
1562	IF ( debug_output ) CALL debug_message( 'radiation_init', 'start' )
1563	!
1564	!-- Activate radiation_interactions according to the existence of vertical surfaces and/or trees
1565	! or if biometeorology output is required for flat surfaces.
1566	!-- The namelist parameter radiation_interactions_on can override this behavior.
1567	!-- (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in
1568	!-- init_surface_arrays.)
1569	IF ( radiation_interactions_on ) THEN
1570	IF ( vertical_surfaces_exist .OR. plant_canopy .OR. biometeorology ) THEN
1571	radiation_interactions = .TRUE.
1572	average_radiation = .TRUE.
1573	ELSE
1574	radiation_interactions_on = .FALSE. !< reset namelist parameter: no interactions
1575	!< calculations necessary in case of flat surface
1576	ENDIF
1577	ELSEIF ( vertical_surfaces_exist .OR. plant_canopy .OR. biometeorology ) THEN
1578	message_string = 'radiation_interactions_on is set to .FALSE. although ' // &
1579	'vertical surfaces and/or trees or biometeorology exist ' // &
1580	'is ON. The model will run without RTM (no shadows, no ' // &
1581	'radiation reflections)'
1582	CALL message( 'init_3d_model', 'PA0348', 0, 1, 0, 6, 0 )
1583	ENDIF
1584	!
1585	!-- Precalculate some time constants
1586	d_hours_day = 1.0_wp / REAL( hours_per_day, KIND = wp )
1587	d_seconds_hour = 1.0_wp / seconds_per_hour
1588
1589	!
1590	!-- If required, initialize radiation interactions between surfaces
1591	!-- via sky-view factors. This must be done before radiation is initialized.
1592	IF ( radiation_interactions ) CALL radiation_interaction_init
1593	!
1594	!-- Allocate array for storing the surface net radiation
1595	IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_net ) .AND. &
1596	surf_lsm_h%ns > 0 ) THEN
1597	ALLOCATE( surf_lsm_h%rad_net(1:surf_lsm_h%ns) )
1598	surf_lsm_h%rad_net = 0.0_wp
1599	ENDIF
1600	IF ( .NOT. ALLOCATED ( surf_usm_h%rad_net ) .AND. &
1601	surf_usm_h%ns > 0 ) THEN
1602	ALLOCATE( surf_usm_h%rad_net(1:surf_usm_h%ns) )
1603	surf_usm_h%rad_net = 0.0_wp
1604	ENDIF
1605	DO l = 0, 3
1606	IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_net ) .AND. &
1607	surf_lsm_v(l)%ns > 0 ) THEN
1608	ALLOCATE( surf_lsm_v(l)%rad_net(1:surf_lsm_v(l)%ns) )
1609	surf_lsm_v(l)%rad_net = 0.0_wp
1610	ENDIF
1611	IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_net ) .AND. &
1612	surf_usm_v(l)%ns > 0 ) THEN
1613	ALLOCATE( surf_usm_v(l)%rad_net(1:surf_usm_v(l)%ns) )
1614	surf_usm_v(l)%rad_net = 0.0_wp
1615	ENDIF
1616	ENDDO
1617
1618
1619	!
1620	!-- Allocate array for storing the surface longwave (out) radiation change
1621	IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_lw_out_change_0 ) .AND. &
1622	surf_lsm_h%ns > 0 ) THEN
1623	ALLOCATE( surf_lsm_h%rad_lw_out_change_0(1:surf_lsm_h%ns) )
1624	surf_lsm_h%rad_lw_out_change_0 = 0.0_wp
1625	ENDIF
1626	IF ( .NOT. ALLOCATED ( surf_usm_h%rad_lw_out_change_0 ) .AND. &
1627	surf_usm_h%ns > 0 ) THEN
1628	ALLOCATE( surf_usm_h%rad_lw_out_change_0(1:surf_usm_h%ns) )
1629	surf_usm_h%rad_lw_out_change_0 = 0.0_wp
1630	ENDIF
1631	DO l = 0, 3
1632	IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_lw_out_change_0 ) .AND. &
1633	surf_lsm_v(l)%ns > 0 ) THEN
1634	ALLOCATE( surf_lsm_v(l)%rad_lw_out_change_0(1:surf_lsm_v(l)%ns) )
1635	surf_lsm_v(l)%rad_lw_out_change_0 = 0.0_wp
1636	ENDIF
1637	IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_lw_out_change_0 ) .AND. &
1638	surf_usm_v(l)%ns > 0 ) THEN
1639	ALLOCATE( surf_usm_v(l)%rad_lw_out_change_0(1:surf_usm_v(l)%ns) )
1640	surf_usm_v(l)%rad_lw_out_change_0 = 0.0_wp
1641	ENDIF
1642	ENDDO
1643
1644	!
1645	!-- Allocate surface arrays for incoming/outgoing short/longwave radiation
1646	IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_sw_in ) .AND. &
1647	surf_lsm_h%ns > 0 ) THEN
1648	ALLOCATE( surf_lsm_h%rad_sw_in(1:surf_lsm_h%ns) )
1649	ALLOCATE( surf_lsm_h%rad_sw_out(1:surf_lsm_h%ns) )
1650	ALLOCATE( surf_lsm_h%rad_sw_dir(1:surf_lsm_h%ns) )
1651	ALLOCATE( surf_lsm_h%rad_sw_dif(1:surf_lsm_h%ns) )
1652	ALLOCATE( surf_lsm_h%rad_sw_ref(1:surf_lsm_h%ns) )
1653	ALLOCATE( surf_lsm_h%rad_sw_res(1:surf_lsm_h%ns) )
1654	ALLOCATE( surf_lsm_h%rad_lw_in(1:surf_lsm_h%ns) )
1655	ALLOCATE( surf_lsm_h%rad_lw_out(1:surf_lsm_h%ns) )
1656	ALLOCATE( surf_lsm_h%rad_lw_dif(1:surf_lsm_h%ns) )
1657	ALLOCATE( surf_lsm_h%rad_lw_ref(1:surf_lsm_h%ns) )
1658	ALLOCATE( surf_lsm_h%rad_lw_res(1:surf_lsm_h%ns) )
1659	surf_lsm_h%rad_sw_in = 0.0_wp
1660	surf_lsm_h%rad_sw_out = 0.0_wp
1661	surf_lsm_h%rad_sw_dir = 0.0_wp
1662	surf_lsm_h%rad_sw_dif = 0.0_wp
1663	surf_lsm_h%rad_sw_ref = 0.0_wp
1664	surf_lsm_h%rad_sw_res = 0.0_wp
1665	surf_lsm_h%rad_lw_in = 0.0_wp
1666	surf_lsm_h%rad_lw_out = 0.0_wp
1667	surf_lsm_h%rad_lw_dif = 0.0_wp
1668	surf_lsm_h%rad_lw_ref = 0.0_wp
1669	surf_lsm_h%rad_lw_res = 0.0_wp
1670	ENDIF
1671	IF ( .NOT. ALLOCATED ( surf_usm_h%rad_sw_in ) .AND. &
1672	surf_usm_h%ns > 0 ) THEN
1673	ALLOCATE( surf_usm_h%rad_sw_in(1:surf_usm_h%ns) )
1674	ALLOCATE( surf_usm_h%rad_sw_out(1:surf_usm_h%ns) )
1675	ALLOCATE( surf_usm_h%rad_sw_dir(1:surf_usm_h%ns) )
1676	ALLOCATE( surf_usm_h%rad_sw_dif(1:surf_usm_h%ns) )
1677	ALLOCATE( surf_usm_h%rad_sw_ref(1:surf_usm_h%ns) )
1678	ALLOCATE( surf_usm_h%rad_sw_res(1:surf_usm_h%ns) )
1679	ALLOCATE( surf_usm_h%rad_lw_in(1:surf_usm_h%ns) )
1680	ALLOCATE( surf_usm_h%rad_lw_out(1:surf_usm_h%ns) )
1681	ALLOCATE( surf_usm_h%rad_lw_dif(1:surf_usm_h%ns) )
1682	ALLOCATE( surf_usm_h%rad_lw_ref(1:surf_usm_h%ns) )
1683	ALLOCATE( surf_usm_h%rad_lw_res(1:surf_usm_h%ns) )
1684	surf_usm_h%rad_sw_in = 0.0_wp
1685	surf_usm_h%rad_sw_out = 0.0_wp
1686	surf_usm_h%rad_sw_dir = 0.0_wp
1687	surf_usm_h%rad_sw_dif = 0.0_wp
1688	surf_usm_h%rad_sw_ref = 0.0_wp
1689	surf_usm_h%rad_sw_res = 0.0_wp
1690	surf_usm_h%rad_lw_in = 0.0_wp
1691	surf_usm_h%rad_lw_out = 0.0_wp
1692	surf_usm_h%rad_lw_dif = 0.0_wp
1693	surf_usm_h%rad_lw_ref = 0.0_wp
1694	surf_usm_h%rad_lw_res = 0.0_wp
1695	ENDIF
1696	DO l = 0, 3
1697	IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_sw_in ) .AND. &
1698	surf_lsm_v(l)%ns > 0 ) THEN
1699	ALLOCATE( surf_lsm_v(l)%rad_sw_in(1:surf_lsm_v(l)%ns) )
1700	ALLOCATE( surf_lsm_v(l)%rad_sw_out(1:surf_lsm_v(l)%ns) )
1701	ALLOCATE( surf_lsm_v(l)%rad_sw_dir(1:surf_lsm_v(l)%ns) )
1702	ALLOCATE( surf_lsm_v(l)%rad_sw_dif(1:surf_lsm_v(l)%ns) )
1703	ALLOCATE( surf_lsm_v(l)%rad_sw_ref(1:surf_lsm_v(l)%ns) )
1704	ALLOCATE( surf_lsm_v(l)%rad_sw_res(1:surf_lsm_v(l)%ns) )
1705
1706	ALLOCATE( surf_lsm_v(l)%rad_lw_in(1:surf_lsm_v(l)%ns) )
1707	ALLOCATE( surf_lsm_v(l)%rad_lw_out(1:surf_lsm_v(l)%ns) )
1708	ALLOCATE( surf_lsm_v(l)%rad_lw_dif(1:surf_lsm_v(l)%ns) )
1709	ALLOCATE( surf_lsm_v(l)%rad_lw_ref(1:surf_lsm_v(l)%ns) )
1710	ALLOCATE( surf_lsm_v(l)%rad_lw_res(1:surf_lsm_v(l)%ns) )
1711
1712	surf_lsm_v(l)%rad_sw_in = 0.0_wp
1713	surf_lsm_v(l)%rad_sw_out = 0.0_wp
1714	surf_lsm_v(l)%rad_sw_dir = 0.0_wp
1715	surf_lsm_v(l)%rad_sw_dif = 0.0_wp
1716	surf_lsm_v(l)%rad_sw_ref = 0.0_wp
1717	surf_lsm_v(l)%rad_sw_res = 0.0_wp
1718
1719	surf_lsm_v(l)%rad_lw_in = 0.0_wp
1720	surf_lsm_v(l)%rad_lw_out = 0.0_wp
1721	surf_lsm_v(l)%rad_lw_dif = 0.0_wp
1722	surf_lsm_v(l)%rad_lw_ref = 0.0_wp
1723	surf_lsm_v(l)%rad_lw_res = 0.0_wp
1724	ENDIF
1725	IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_sw_in ) .AND. &
1726	surf_usm_v(l)%ns > 0 ) THEN
1727	ALLOCATE( surf_usm_v(l)%rad_sw_in(1:surf_usm_v(l)%ns) )
1728	ALLOCATE( surf_usm_v(l)%rad_sw_out(1:surf_usm_v(l)%ns) )
1729	ALLOCATE( surf_usm_v(l)%rad_sw_dir(1:surf_usm_v(l)%ns) )
1730	ALLOCATE( surf_usm_v(l)%rad_sw_dif(1:surf_usm_v(l)%ns) )
1731	ALLOCATE( surf_usm_v(l)%rad_sw_ref(1:surf_usm_v(l)%ns) )
1732	ALLOCATE( surf_usm_v(l)%rad_sw_res(1:surf_usm_v(l)%ns) )
1733	ALLOCATE( surf_usm_v(l)%rad_lw_in(1:surf_usm_v(l)%ns) )
1734	ALLOCATE( surf_usm_v(l)%rad_lw_out(1:surf_usm_v(l)%ns) )
1735	ALLOCATE( surf_usm_v(l)%rad_lw_dif(1:surf_usm_v(l)%ns) )
1736	ALLOCATE( surf_usm_v(l)%rad_lw_ref(1:surf_usm_v(l)%ns) )
1737	ALLOCATE( surf_usm_v(l)%rad_lw_res(1:surf_usm_v(l)%ns) )
1738	surf_usm_v(l)%rad_sw_in = 0.0_wp
1739	surf_usm_v(l)%rad_sw_out = 0.0_wp
1740	surf_usm_v(l)%rad_sw_dir = 0.0_wp
1741	surf_usm_v(l)%rad_sw_dif = 0.0_wp
1742	surf_usm_v(l)%rad_sw_ref = 0.0_wp
1743	surf_usm_v(l)%rad_sw_res = 0.0_wp
1744	surf_usm_v(l)%rad_lw_in = 0.0_wp
1745	surf_usm_v(l)%rad_lw_out = 0.0_wp
1746	surf_usm_v(l)%rad_lw_dif = 0.0_wp
1747	surf_usm_v(l)%rad_lw_ref = 0.0_wp
1748	surf_usm_v(l)%rad_lw_res = 0.0_wp
1749	ENDIF
1750	ENDDO
1751	!
1752	!-- Fix net radiation in case of radiation_scheme = 'constant'
1753	IF ( radiation_scheme == 'constant' ) THEN
1754	IF ( ALLOCATED( surf_lsm_h%rad_net ) ) &
1755	surf_lsm_h%rad_net = net_radiation
1756	IF ( ALLOCATED( surf_usm_h%rad_net ) ) &
1757	surf_usm_h%rad_net = net_radiation
1758	!
1759	!-- Todo: weight with inclination angle
1760	DO l = 0, 3
1761	IF ( ALLOCATED( surf_lsm_v(l)%rad_net ) ) &
1762	surf_lsm_v(l)%rad_net = net_radiation
1763	IF ( ALLOCATED( surf_usm_v(l)%rad_net ) ) &
1764	surf_usm_v(l)%rad_net = net_radiation
1765	ENDDO
1766	! radiation = .FALSE.
1767	!
1768	!-- Calculate orbital constants
1769	ELSE
1770	decl_1 = SIN(23.45_wp * pi / 180.0_wp)
1771	decl_2 = 2.0_wp * pi / 365.0_wp
1772	decl_3 = decl_2 * 81.0_wp
1773	lat = latitude * pi / 180.0_wp
1774	lon = longitude * pi / 180.0_wp
1775	ENDIF
1776
1777	IF ( radiation_scheme == 'clear-sky' .OR. &
1778	radiation_scheme == 'constant' .OR. &
1779	radiation_scheme == 'external' ) THEN
1780	!
1781	!-- Allocate arrays for incoming/outgoing short/longwave radiation
1782	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
1783	ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
1784	ENDIF
1785	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
1786	ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
1787	ENDIF
1788
1789	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
1790	ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
1791	ENDIF
1792	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
1793	ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
1794	ENDIF
1795
1796	!
1797	!-- Allocate average arrays for incoming/outgoing short/longwave radiation
1798	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
1799	ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
1800	ENDIF
1801	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
1802	ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
1803	ENDIF
1804
1805	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
1806	ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
1807	ENDIF
1808	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
1809	ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
1810	ENDIF
1811	!
1812	!-- Allocate arrays for broadband albedo, and level 1 initialization
1813	!-- via namelist paramter, unless not already allocated.
1814	IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) ) THEN
1815	ALLOCATE( surf_lsm_h%albedo(0:2,1:surf_lsm_h%ns) )
1816	surf_lsm_h%albedo = albedo
1817	ENDIF
1818	IF ( .NOT. ALLOCATED(surf_usm_h%albedo) ) THEN
1819	ALLOCATE( surf_usm_h%albedo(0:2,1:surf_usm_h%ns) )
1820	surf_usm_h%albedo = albedo
1821	ENDIF
1822
1823	DO l = 0, 3
1824	IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) ) THEN
1825	ALLOCATE( surf_lsm_v(l)%albedo(0:2,1:surf_lsm_v(l)%ns) )
1826	surf_lsm_v(l)%albedo = albedo
1827	ENDIF
1828	IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) ) THEN
1829	ALLOCATE( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns) )
1830	surf_usm_v(l)%albedo = albedo
1831	ENDIF
1832	ENDDO
1833	!
1834	!-- Level 2 initialization of broadband albedo via given albedo_type.
1835	!-- Only if albedo_type is non-zero. In case of urban surface and
1836	!-- input data is read from ASCII file, albedo_type will be zero, so that
1837	!-- albedo won't be overwritten.
1838	DO m = 1, surf_lsm_h%ns
1839	IF ( surf_lsm_h%albedo_type(ind_veg_wall,m) /= 0 ) &
1840	surf_lsm_h%albedo(ind_veg_wall,m) = &
1841	albedo_pars(0,surf_lsm_h%albedo_type(ind_veg_wall,m))
1842	IF ( surf_lsm_h%albedo_type(ind_pav_green,m) /= 0 ) &
1843	surf_lsm_h%albedo(ind_pav_green,m) = &
1844	albedo_pars(0,surf_lsm_h%albedo_type(ind_pav_green,m))
1845	IF ( surf_lsm_h%albedo_type(ind_wat_win,m) /= 0 ) &
1846	surf_lsm_h%albedo(ind_wat_win,m) = &
1847	albedo_pars(0,surf_lsm_h%albedo_type(ind_wat_win,m))
1848	ENDDO
1849	DO m = 1, surf_usm_h%ns
1850	IF ( surf_usm_h%albedo_type(ind_veg_wall,m) /= 0 ) &
1851	surf_usm_h%albedo(ind_veg_wall,m) = &
1852	albedo_pars(0,surf_usm_h%albedo_type(ind_veg_wall,m))
1853	IF ( surf_usm_h%albedo_type(ind_pav_green,m) /= 0 ) &
1854	surf_usm_h%albedo(ind_pav_green,m) = &
1855	albedo_pars(0,surf_usm_h%albedo_type(ind_pav_green,m))
1856	IF ( surf_usm_h%albedo_type(ind_wat_win,m) /= 0 ) &
1857	surf_usm_h%albedo(ind_wat_win,m) = &
1858	albedo_pars(0,surf_usm_h%albedo_type(ind_wat_win,m))
1859	ENDDO
1860
1861	DO l = 0, 3
1862	DO m = 1, surf_lsm_v(l)%ns
1863	IF ( surf_lsm_v(l)%albedo_type(ind_veg_wall,m) /= 0 ) &
1864	surf_lsm_v(l)%albedo(ind_veg_wall,m) = &
1865	albedo_pars(0,surf_lsm_v(l)%albedo_type(ind_veg_wall,m))
1866	IF ( surf_lsm_v(l)%albedo_type(ind_pav_green,m) /= 0 ) &
1867	surf_lsm_v(l)%albedo(ind_pav_green,m) = &
1868	albedo_pars(0,surf_lsm_v(l)%albedo_type(ind_pav_green,m))
1869	IF ( surf_lsm_v(l)%albedo_type(ind_wat_win,m) /= 0 ) &
1870	surf_lsm_v(l)%albedo(ind_wat_win,m) = &
1871	albedo_pars(0,surf_lsm_v(l)%albedo_type(ind_wat_win,m))
1872	ENDDO
1873	DO m = 1, surf_usm_v(l)%ns
1874	IF ( surf_usm_v(l)%albedo_type(ind_veg_wall,m) /= 0 ) &
1875	surf_usm_v(l)%albedo(ind_veg_wall,m) = &
1876	albedo_pars(0,surf_usm_v(l)%albedo_type(ind_veg_wall,m))
1877	IF ( surf_usm_v(l)%albedo_type(ind_pav_green,m) /= 0 ) &
1878	surf_usm_v(l)%albedo(ind_pav_green,m) = &
1879	albedo_pars(0,surf_usm_v(l)%albedo_type(ind_pav_green,m))
1880	IF ( surf_usm_v(l)%albedo_type(ind_wat_win,m) /= 0 ) &
1881	surf_usm_v(l)%albedo(ind_wat_win,m) = &
1882	albedo_pars(0,surf_usm_v(l)%albedo_type(ind_wat_win,m))
1883	ENDDO
1884	ENDDO
1885
1886	!
1887	!-- Level 3 initialization at grid points where albedo type is zero.
1888	!-- This case, albedo is taken from file. In case of constant radiation
1889	!-- or clear sky, only broadband albedo is given.
1890	IF ( albedo_pars_f%from_file ) THEN
1891	!
1892	!-- Horizontal surfaces
1893	DO m = 1, surf_lsm_h%ns
1894	i = surf_lsm_h%i(m)
1895	j = surf_lsm_h%j(m)
1896	IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill ) THEN
1897	surf_lsm_h%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
1898	surf_lsm_h%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
1899	surf_lsm_h%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
1900	ENDIF
1901	ENDDO
1902	DO m = 1, surf_usm_h%ns
1903	i = surf_usm_h%i(m)
1904	j = surf_usm_h%j(m)
1905	IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill ) THEN
1906	surf_usm_h%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
1907	surf_usm_h%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
1908	surf_usm_h%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
1909	ENDIF
1910	ENDDO
1911	!
1912	!-- Vertical surfaces
1913	DO l = 0, 3
1914
1915	ioff = surf_lsm_v(l)%ioff
1916	joff = surf_lsm_v(l)%joff
1917	DO m = 1, surf_lsm_v(l)%ns
1918	i = surf_lsm_v(l)%i(m) + ioff
1919	j = surf_lsm_v(l)%j(m) + joff
1920	IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill ) THEN
1921	surf_lsm_v(l)%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
1922	surf_lsm_v(l)%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
1923	surf_lsm_v(l)%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
1924	ENDIF
1925	ENDDO
1926
1927	ioff = surf_usm_v(l)%ioff
1928	joff = surf_usm_v(l)%joff
1929	DO m = 1, surf_usm_v(l)%ns
1930	i = surf_usm_v(l)%i(m) + ioff
1931	j = surf_usm_v(l)%j(m) + joff
1932	IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill ) THEN
1933	surf_usm_v(l)%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
1934	surf_usm_v(l)%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
1935	surf_usm_v(l)%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
1936	ENDIF
1937	ENDDO
1938	ENDDO
1939
1940	ENDIF
1941	!
1942	!-- Read explicit albedo values from building surface pars. If present,
1943	!-- they override all less specific albedo values and force a albedo_type
1944	!-- to zero in order to take effect.
1945	IF ( building_surface_pars_f%from_file ) THEN
1946	DO m = 1, surf_usm_h%ns
1947	i = surf_usm_h%i(m)
1948	j = surf_usm_h%j(m)
1949	k = surf_usm_h%k(m)
1950	!
1951	!-- Iterate over surfaces in column, check height and orientation
1952	DO is = building_surface_pars_f%index_ji(1,j,i), &
1953	building_surface_pars_f%index_ji(2,j,i)
1954	IF ( building_surface_pars_f%coords(4,is) == -surf_usm_h%koff .AND. &
1955	building_surface_pars_f%coords(1,is) == k ) THEN
1956
1957	IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /= &
1958	building_surface_pars_f%fill ) THEN
1959	surf_usm_h%albedo(ind_veg_wall,m) = &
1960	building_surface_pars_f%pars(ind_s_alb_b_wall,is)
1961	surf_usm_h%albedo_type(ind_veg_wall,m) = 0
1962	ENDIF
1963
1964	IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /= &
1965	building_surface_pars_f%fill ) THEN
1966	surf_usm_h%albedo(ind_wat_win,m) = &
1967	building_surface_pars_f%pars(ind_s_alb_b_win,is)
1968	surf_usm_h%albedo_type(ind_wat_win,m) = 0
1969	ENDIF
1970
1971	IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /= &
1972	building_surface_pars_f%fill ) THEN
1973	surf_usm_h%albedo(ind_pav_green,m) = &
1974	building_surface_pars_f%pars(ind_s_alb_b_green,is)
1975	surf_usm_h%albedo_type(ind_pav_green,m) = 0
1976	ENDIF
1977
1978	EXIT ! surface was found and processed
1979	ENDIF
1980	ENDDO
1981	ENDDO
1982
1983	DO l = 0, 3
1984	DO m = 1, surf_usm_v(l)%ns
1985	i = surf_usm_v(l)%i(m)
1986	j = surf_usm_v(l)%j(m)
1987	k = surf_usm_v(l)%k(m)
1988	!
1989	!-- Iterate over surfaces in column, check height and orientation
1990	DO is = building_surface_pars_f%index_ji(1,j,i), &
1991	building_surface_pars_f%index_ji(2,j,i)
1992	IF ( building_surface_pars_f%coords(5,is) == -surf_usm_v(l)%joff .AND. &
1993	building_surface_pars_f%coords(6,is) == -surf_usm_v(l)%ioff .AND. &
1994	building_surface_pars_f%coords(1,is) == k ) THEN
1995
1996	IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /= &
1997	building_surface_pars_f%fill ) THEN
1998	surf_usm_v(l)%albedo(ind_veg_wall,m) = &
1999	building_surface_pars_f%pars(ind_s_alb_b_wall,is)
2000	surf_usm_v(l)%albedo_type(ind_veg_wall,m) = 0
2001	ENDIF
2002
2003	IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /= &
2004	building_surface_pars_f%fill ) THEN
2005	surf_usm_v(l)%albedo(ind_wat_win,m) = &
2006	building_surface_pars_f%pars(ind_s_alb_b_win,is)
2007	surf_usm_v(l)%albedo_type(ind_wat_win,m) = 0
2008	ENDIF
2009
2010	IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /= &
2011	building_surface_pars_f%fill ) THEN
2012	surf_usm_v(l)%albedo(ind_pav_green,m) = &
2013	building_surface_pars_f%pars(ind_s_alb_b_green,is)
2014	surf_usm_v(l)%albedo_type(ind_pav_green,m) = 0
2015	ENDIF
2016
2017	EXIT ! surface was found and processed
2018	ENDIF
2019	ENDDO
2020	ENDDO
2021	ENDDO
2022	ENDIF
2023	!
2024	!-- Initialization actions for RRTMG
2025	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
2026	#if defined ( __rrtmg )
2027	!
2028	!-- Allocate albedos for short/longwave radiation, horizontal surfaces
2029	!-- for wall/green/window (USM) or vegetation/pavement/water surfaces
2030	!-- (LSM).
2031	ALLOCATE ( surf_lsm_h%aldif(0:2,1:surf_lsm_h%ns) )
2032	ALLOCATE ( surf_lsm_h%aldir(0:2,1:surf_lsm_h%ns) )
2033	ALLOCATE ( surf_lsm_h%asdif(0:2,1:surf_lsm_h%ns) )
2034	ALLOCATE ( surf_lsm_h%asdir(0:2,1:surf_lsm_h%ns) )
2035	ALLOCATE ( surf_lsm_h%rrtm_aldif(0:2,1:surf_lsm_h%ns) )
2036	ALLOCATE ( surf_lsm_h%rrtm_aldir(0:2,1:surf_lsm_h%ns) )
2037	ALLOCATE ( surf_lsm_h%rrtm_asdif(0:2,1:surf_lsm_h%ns) )
2038	ALLOCATE ( surf_lsm_h%rrtm_asdir(0:2,1:surf_lsm_h%ns) )
2039
2040	ALLOCATE ( surf_usm_h%aldif(0:2,1:surf_usm_h%ns) )
2041	ALLOCATE ( surf_usm_h%aldir(0:2,1:surf_usm_h%ns) )
2042	ALLOCATE ( surf_usm_h%asdif(0:2,1:surf_usm_h%ns) )
2043	ALLOCATE ( surf_usm_h%asdir(0:2,1:surf_usm_h%ns) )
2044	ALLOCATE ( surf_usm_h%rrtm_aldif(0:2,1:surf_usm_h%ns) )
2045	ALLOCATE ( surf_usm_h%rrtm_aldir(0:2,1:surf_usm_h%ns) )
2046	ALLOCATE ( surf_usm_h%rrtm_asdif(0:2,1:surf_usm_h%ns) )
2047	ALLOCATE ( surf_usm_h%rrtm_asdir(0:2,1:surf_usm_h%ns) )
2048
2049	!
2050	!-- Allocate broadband albedo (temporary for the current radiation
2051	!-- implementations)
2052	IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) ) &
2053	ALLOCATE( surf_lsm_h%albedo(0:2,1:surf_lsm_h%ns) )
2054	IF ( .NOT. ALLOCATED(surf_usm_h%albedo) ) &
2055	ALLOCATE( surf_usm_h%albedo(0:2,1:surf_usm_h%ns) )
2056
2057	!
2058	!-- Allocate albedos for short/longwave radiation, vertical surfaces
2059	DO l = 0, 3
2060
2061	ALLOCATE ( surf_lsm_v(l)%aldif(0:2,1:surf_lsm_v(l)%ns) )
2062	ALLOCATE ( surf_lsm_v(l)%aldir(0:2,1:surf_lsm_v(l)%ns) )
2063	ALLOCATE ( surf_lsm_v(l)%asdif(0:2,1:surf_lsm_v(l)%ns) )
2064	ALLOCATE ( surf_lsm_v(l)%asdir(0:2,1:surf_lsm_v(l)%ns) )
2065
2066	ALLOCATE ( surf_lsm_v(l)%rrtm_aldif(0:2,1:surf_lsm_v(l)%ns) )
2067	ALLOCATE ( surf_lsm_v(l)%rrtm_aldir(0:2,1:surf_lsm_v(l)%ns) )
2068	ALLOCATE ( surf_lsm_v(l)%rrtm_asdif(0:2,1:surf_lsm_v(l)%ns) )
2069	ALLOCATE ( surf_lsm_v(l)%rrtm_asdir(0:2,1:surf_lsm_v(l)%ns) )
2070
2071	ALLOCATE ( surf_usm_v(l)%aldif(0:2,1:surf_usm_v(l)%ns) )
2072	ALLOCATE ( surf_usm_v(l)%aldir(0:2,1:surf_usm_v(l)%ns) )
2073	ALLOCATE ( surf_usm_v(l)%asdif(0:2,1:surf_usm_v(l)%ns) )
2074	ALLOCATE ( surf_usm_v(l)%asdir(0:2,1:surf_usm_v(l)%ns) )
2075
2076	ALLOCATE ( surf_usm_v(l)%rrtm_aldif(0:2,1:surf_usm_v(l)%ns) )
2077	ALLOCATE ( surf_usm_v(l)%rrtm_aldir(0:2,1:surf_usm_v(l)%ns) )
2078	ALLOCATE ( surf_usm_v(l)%rrtm_asdif(0:2,1:surf_usm_v(l)%ns) )
2079	ALLOCATE ( surf_usm_v(l)%rrtm_asdir(0:2,1:surf_usm_v(l)%ns) )
2080	!
2081	!-- Allocate broadband albedo (temporary for the current radiation
2082	!-- implementations)
2083	IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) ) &
2084	ALLOCATE( surf_lsm_v(l)%albedo(0:2,1:surf_lsm_v(l)%ns) )
2085	IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) ) &
2086	ALLOCATE( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns) )
2087
2088	ENDDO
2089	!
2090	!-- Level 1 initialization of spectral albedos via namelist
2091	!-- paramters. Please note, this case all surface tiles are initialized
2092	!-- the same.
2093	IF ( surf_lsm_h%ns > 0 ) THEN
2094	surf_lsm_h%aldif = albedo_lw_dif
2095	surf_lsm_h%aldir = albedo_lw_dir
2096	surf_lsm_h%asdif = albedo_sw_dif
2097	surf_lsm_h%asdir = albedo_sw_dir
2098	surf_lsm_h%albedo = albedo_sw_dif
2099	ENDIF
2100	IF ( surf_usm_h%ns > 0 ) THEN
2101	IF ( surf_usm_h%albedo_from_ascii ) THEN
2102	surf_usm_h%aldif = surf_usm_h%albedo
2103	surf_usm_h%aldir = surf_usm_h%albedo
2104	surf_usm_h%asdif = surf_usm_h%albedo
2105	surf_usm_h%asdir = surf_usm_h%albedo
2106	ELSE
2107	surf_usm_h%aldif = albedo_lw_dif
2108	surf_usm_h%aldir = albedo_lw_dir
2109	surf_usm_h%asdif = albedo_sw_dif
2110	surf_usm_h%asdir = albedo_sw_dir
2111	surf_usm_h%albedo = albedo_sw_dif
2112	ENDIF
2113	ENDIF
2114
2115	DO l = 0, 3
2116
2117	IF ( surf_lsm_v(l)%ns > 0 ) THEN
2118	surf_lsm_v(l)%aldif = albedo_lw_dif
2119	surf_lsm_v(l)%aldir = albedo_lw_dir
2120	surf_lsm_v(l)%asdif = albedo_sw_dif
2121	surf_lsm_v(l)%asdir = albedo_sw_dir
2122	surf_lsm_v(l)%albedo = albedo_sw_dif
2123	ENDIF
2124
2125	IF ( surf_usm_v(l)%ns > 0 ) THEN
2126	IF ( surf_usm_v(l)%albedo_from_ascii ) THEN
2127	surf_usm_v(l)%aldif = surf_usm_v(l)%albedo
2128	surf_usm_v(l)%aldir = surf_usm_v(l)%albedo
2129	surf_usm_v(l)%asdif = surf_usm_v(l)%albedo
2130	surf_usm_v(l)%asdir = surf_usm_v(l)%albedo
2131	ELSE
2132	surf_usm_v(l)%aldif = albedo_lw_dif
2133	surf_usm_v(l)%aldir = albedo_lw_dir
2134	surf_usm_v(l)%asdif = albedo_sw_dif
2135	surf_usm_v(l)%asdir = albedo_sw_dir
2136	ENDIF
2137	ENDIF
2138	ENDDO
2139
2140	!
2141	!-- Level 2 initialization of spectral albedos via albedo_type.
2142	!-- Please note, for natural- and urban-type surfaces, a tile approach
2143	!-- is applied so that the resulting albedo is calculated via the weighted
2144	!-- average of respective surface fractions.
2145	DO m = 1, surf_lsm_h%ns
2146	!
2147	!-- Spectral albedos for vegetation/pavement/water surfaces
2148	DO ind_type = 0, 2
2149	IF ( surf_lsm_h%albedo_type(ind_type,m) /= 0 ) THEN
2150	surf_lsm_h%aldif(ind_type,m) = &
2151	albedo_pars(1,surf_lsm_h%albedo_type(ind_type,m))
2152	surf_lsm_h%asdif(ind_type,m) = &
2153	albedo_pars(2,surf_lsm_h%albedo_type(ind_type,m))
2154	surf_lsm_h%aldir(ind_type,m) = &
2155	albedo_pars(1,surf_lsm_h%albedo_type(ind_type,m))
2156	surf_lsm_h%asdir(ind_type,m) = &
2157	albedo_pars(2,surf_lsm_h%albedo_type(ind_type,m))
2158	surf_lsm_h%albedo(ind_type,m) = &
2159	albedo_pars(0,surf_lsm_h%albedo_type(ind_type,m))
2160	ENDIF
2161	ENDDO
2162
2163	ENDDO
2164	!
2165	!-- For urban surface only if albedo has not been already initialized
2166	!-- in the urban-surface model via the ASCII file.
2167	IF ( .NOT. surf_usm_h%albedo_from_ascii ) THEN
2168	DO m = 1, surf_usm_h%ns
2169	!
2170	!-- Spectral albedos for wall/green/window surfaces
2171	DO ind_type = 0, 2
2172	IF ( surf_usm_h%albedo_type(ind_type,m) /= 0 ) THEN
2173	surf_usm_h%aldif(ind_type,m) = &
2174	albedo_pars(1,surf_usm_h%albedo_type(ind_type,m))
2175	surf_usm_h%asdif(ind_type,m) = &
2176	albedo_pars(2,surf_usm_h%albedo_type(ind_type,m))
2177	surf_usm_h%aldir(ind_type,m) = &
2178	albedo_pars(1,surf_usm_h%albedo_type(ind_type,m))
2179	surf_usm_h%asdir(ind_type,m) = &
2180	albedo_pars(2,surf_usm_h%albedo_type(ind_type,m))
2181	surf_usm_h%albedo(ind_type,m) = &
2182	albedo_pars(0,surf_usm_h%albedo_type(ind_type,m))
2183	ENDIF
2184	ENDDO
2185
2186	ENDDO
2187	ENDIF
2188
2189	DO l = 0, 3
2190
2191	DO m = 1, surf_lsm_v(l)%ns
2192	!
2193	!-- Spectral albedos for vegetation/pavement/water surfaces
2194	DO ind_type = 0, 2
2195	IF ( surf_lsm_v(l)%albedo_type(ind_type,m) /= 0 ) THEN
2196	surf_lsm_v(l)%aldif(ind_type,m) = &
2197	albedo_pars(1,surf_lsm_v(l)%albedo_type(ind_type,m))
2198	surf_lsm_v(l)%asdif(ind_type,m) = &
2199	albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_type,m))
2200	surf_lsm_v(l)%aldir(ind_type,m) = &
2201	albedo_pars(1,surf_lsm_v(l)%albedo_type(ind_type,m))
2202	surf_lsm_v(l)%asdir(ind_type,m) = &
2203	albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_type,m))
2204	surf_lsm_v(l)%albedo(ind_type,m) = &
2205	albedo_pars(0,surf_lsm_v(l)%albedo_type(ind_type,m))
2206	ENDIF
2207	ENDDO
2208	ENDDO
2209	!
2210	!-- For urban surface only if albedo has not been already initialized
2211	!-- in the urban-surface model via the ASCII file.
2212	IF ( .NOT. surf_usm_v(l)%albedo_from_ascii ) THEN
2213	DO m = 1, surf_usm_v(l)%ns
2214	!
2215	!-- Spectral albedos for wall/green/window surfaces
2216	DO ind_type = 0, 2
2217	IF ( surf_usm_v(l)%albedo_type(ind_type,m) /= 0 ) THEN
2218	surf_usm_v(l)%aldif(ind_type,m) = &
2219	albedo_pars(1,surf_usm_v(l)%albedo_type(ind_type,m))
2220	surf_usm_v(l)%asdif(ind_type,m) = &
2221	albedo_pars(2,surf_usm_v(l)%albedo_type(ind_type,m))
2222	surf_usm_v(l)%aldir(ind_type,m) = &
2223	albedo_pars(1,surf_usm_v(l)%albedo_type(ind_type,m))
2224	surf_usm_v(l)%asdir(ind_type,m) = &
2225	albedo_pars(2,surf_usm_v(l)%albedo_type(ind_type,m))
2226	surf_usm_v(l)%albedo(ind_type,m) = &
2227	albedo_pars(0,surf_usm_v(l)%albedo_type(ind_type,m))
2228	ENDIF
2229	ENDDO
2230
2231	ENDDO
2232	ENDIF
2233	ENDDO
2234	!
2235	!-- Level 3 initialization at grid points where albedo type is zero.
2236	!-- This case, spectral albedos are taken from file if available
2237	IF ( albedo_pars_f%from_file ) THEN
2238	!
2239	!-- Horizontal
2240	DO m = 1, surf_lsm_h%ns
2241	i = surf_lsm_h%i(m)
2242	j = surf_lsm_h%j(m)
2243	!
2244	!-- Spectral albedos for vegetation/pavement/water surfaces
2245	DO ind_type = 0, 2
2246	IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill ) &
2247	surf_lsm_h%albedo(ind_type,m) = &
2248	albedo_pars_f%pars_xy(0,j,i)
2249	IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill ) &
2250	surf_lsm_h%aldir(ind_type,m) = &
2251	albedo_pars_f%pars_xy(1,j,i)
2252	IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill ) &
2253	surf_lsm_h%aldif(ind_type,m) = &
2254	albedo_pars_f%pars_xy(1,j,i)
2255	IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill ) &
2256	surf_lsm_h%asdir(ind_type,m) = &
2257	albedo_pars_f%pars_xy(2,j,i)
2258	IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill ) &
2259	surf_lsm_h%asdif(ind_type,m) = &
2260	albedo_pars_f%pars_xy(2,j,i)
2261	ENDDO
2262	ENDDO
2263	!
2264	!-- For urban surface only if albedo has not been already initialized
2265	!-- in the urban-surface model via the ASCII file.
2266	IF ( .NOT. surf_usm_h%albedo_from_ascii ) THEN
2267	DO m = 1, surf_usm_h%ns
2268	i = surf_usm_h%i(m)
2269	j = surf_usm_h%j(m)
2270	!
2271	!-- Broadband albedos for wall/green/window surfaces
2272	DO ind_type = 0, 2
2273	IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )&
2274	surf_usm_h%albedo(ind_type,m) = &
2275	albedo_pars_f%pars_xy(0,j,i)
2276	ENDDO
2277	!
2278	!-- Spectral albedos especially for building wall surfaces
2279	IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill ) THEN
2280	surf_usm_h%aldir(ind_veg_wall,m) = &
2281	albedo_pars_f%pars_xy(1,j,i)
2282	surf_usm_h%aldif(ind_veg_wall,m) = &
2283	albedo_pars_f%pars_xy(1,j,i)
2284	ENDIF
2285	IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill ) THEN
2286	surf_usm_h%asdir(ind_veg_wall,m) = &
2287	albedo_pars_f%pars_xy(2,j,i)
2288	surf_usm_h%asdif(ind_veg_wall,m) = &
2289	albedo_pars_f%pars_xy(2,j,i)
2290	ENDIF
2291	!
2292	!-- Spectral albedos especially for building green surfaces
2293	IF ( albedo_pars_f%pars_xy(3,j,i) /= albedo_pars_f%fill ) THEN
2294	surf_usm_h%aldir(ind_pav_green,m) = &
2295	albedo_pars_f%pars_xy(3,j,i)
2296	surf_usm_h%aldif(ind_pav_green,m) = &
2297	albedo_pars_f%pars_xy(3,j,i)
2298	ENDIF
2299	IF ( albedo_pars_f%pars_xy(4,j,i) /= albedo_pars_f%fill ) THEN
2300	surf_usm_h%asdir(ind_pav_green,m) = &
2301	albedo_pars_f%pars_xy(4,j,i)
2302	surf_usm_h%asdif(ind_pav_green,m) = &
2303	albedo_pars_f%pars_xy(4,j,i)
2304	ENDIF
2305	!
2306	!-- Spectral albedos especially for building window surfaces
2307	IF ( albedo_pars_f%pars_xy(5,j,i) /= albedo_pars_f%fill ) THEN
2308	surf_usm_h%aldir(ind_wat_win,m) = &
2309	albedo_pars_f%pars_xy(5,j,i)
2310	surf_usm_h%aldif(ind_wat_win,m) = &
2311	albedo_pars_f%pars_xy(5,j,i)
2312	ENDIF
2313	IF ( albedo_pars_f%pars_xy(6,j,i) /= albedo_pars_f%fill ) THEN
2314	surf_usm_h%asdir(ind_wat_win,m) = &
2315	albedo_pars_f%pars_xy(6,j,i)
2316	surf_usm_h%asdif(ind_wat_win,m) = &
2317	albedo_pars_f%pars_xy(6,j,i)
2318	ENDIF
2319
2320	ENDDO
2321	ENDIF
2322	!
2323	!-- Vertical
2324	DO l = 0, 3
2325	ioff = surf_lsm_v(l)%ioff
2326	joff = surf_lsm_v(l)%joff
2327
2328	DO m = 1, surf_lsm_v(l)%ns
2329	i = surf_lsm_v(l)%i(m)
2330	j = surf_lsm_v(l)%j(m)
2331	!
2332	!-- Spectral albedos for vegetation/pavement/water surfaces
2333	DO ind_type = 0, 2
2334	IF ( albedo_pars_f%pars_xy(0,j+joff,i+ioff) /= &
2335	albedo_pars_f%fill ) &
2336	surf_lsm_v(l)%albedo(ind_type,m) = &
2337	albedo_pars_f%pars_xy(0,j+joff,i+ioff)
2338	IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /= &
2339	albedo_pars_f%fill ) &
2340	surf_lsm_v(l)%aldir(ind_type,m) = &
2341	albedo_pars_f%pars_xy(1,j+joff,i+ioff)
2342	IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /= &
2343	albedo_pars_f%fill ) &
2344	surf_lsm_v(l)%aldif(ind_type,m) = &
2345	albedo_pars_f%pars_xy(1,j+joff,i+ioff)
2346	IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /= &
2347	albedo_pars_f%fill ) &
2348	surf_lsm_v(l)%asdir(ind_type,m) = &
2349	albedo_pars_f%pars_xy(2,j+joff,i+ioff)
2350	IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /= &
2351	albedo_pars_f%fill ) &
2352	surf_lsm_v(l)%asdif(ind_type,m) = &
2353	albedo_pars_f%pars_xy(2,j+joff,i+ioff)
2354	ENDDO
2355	ENDDO
2356	!
2357	!-- For urban surface only if albedo has not been already initialized
2358	!-- in the urban-surface model via the ASCII file.
2359	IF ( .NOT. surf_usm_v(l)%albedo_from_ascii ) THEN
2360	ioff = surf_usm_v(l)%ioff
2361	joff = surf_usm_v(l)%joff
2362
2363	DO m = 1, surf_usm_v(l)%ns
2364	i = surf_usm_v(l)%i(m)
2365	j = surf_usm_v(l)%j(m)
2366	!
2367	!-- Broadband albedos for wall/green/window surfaces
2368	DO ind_type = 0, 2
2369	IF ( albedo_pars_f%pars_xy(0,j+joff,i+ioff) /= &
2370	albedo_pars_f%fill ) &
2371	surf_usm_v(l)%albedo(ind_type,m) = &
2372	albedo_pars_f%pars_xy(0,j+joff,i+ioff)
2373	ENDDO
2374	!
2375	!-- Spectral albedos especially for building wall surfaces
2376	IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /= &
2377	albedo_pars_f%fill ) THEN
2378	surf_usm_v(l)%aldir(ind_veg_wall,m) = &
2379	albedo_pars_f%pars_xy(1,j+joff,i+ioff)
2380	surf_usm_v(l)%aldif(ind_veg_wall,m) = &
2381	albedo_pars_f%pars_xy(1,j+joff,i+ioff)
2382	ENDIF
2383	IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /= &
2384	albedo_pars_f%fill ) THEN
2385	surf_usm_v(l)%asdir(ind_veg_wall,m) = &
2386	albedo_pars_f%pars_xy(2,j+joff,i+ioff)
2387	surf_usm_v(l)%asdif(ind_veg_wall,m) = &
2388	albedo_pars_f%pars_xy(2,j+joff,i+ioff)
2389	ENDIF
2390	!
2391	!-- Spectral albedos especially for building green surfaces
2392	IF ( albedo_pars_f%pars_xy(3,j+joff,i+ioff) /= &
2393	albedo_pars_f%fill ) THEN
2394	surf_usm_v(l)%aldir(ind_pav_green,m) = &
2395	albedo_pars_f%pars_xy(3,j+joff,i+ioff)
2396	surf_usm_v(l)%aldif(ind_pav_green,m) = &
2397	albedo_pars_f%pars_xy(3,j+joff,i+ioff)
2398	ENDIF
2399	IF ( albedo_pars_f%pars_xy(4,j+joff,i+ioff) /= &
2400	albedo_pars_f%fill ) THEN
2401	surf_usm_v(l)%asdir(ind_pav_green,m) = &
2402	albedo_pars_f%pars_xy(4,j+joff,i+ioff)
2403	surf_usm_v(l)%asdif(ind_pav_green,m) = &
2404	albedo_pars_f%pars_xy(4,j+joff,i+ioff)
2405	ENDIF
2406	!
2407	!-- Spectral albedos especially for building window surfaces
2408	IF ( albedo_pars_f%pars_xy(5,j+joff,i+ioff) /= &
2409	albedo_pars_f%fill ) THEN
2410	surf_usm_v(l)%aldir(ind_wat_win,m) = &
2411	albedo_pars_f%pars_xy(5,j+joff,i+ioff)
2412	surf_usm_v(l)%aldif(ind_wat_win,m) = &
2413	albedo_pars_f%pars_xy(5,j+joff,i+ioff)
2414	ENDIF
2415	IF ( albedo_pars_f%pars_xy(6,j+joff,i+ioff) /= &
2416	albedo_pars_f%fill ) THEN
2417	surf_usm_v(l)%asdir(ind_wat_win,m) = &
2418	albedo_pars_f%pars_xy(6,j+joff,i+ioff)
2419	surf_usm_v(l)%asdif(ind_wat_win,m) = &
2420	albedo_pars_f%pars_xy(6,j+joff,i+ioff)
2421	ENDIF
2422	ENDDO
2423	ENDIF
2424	ENDDO
2425
2426	ENDIF
2427	!
2428	!-- Read explicit albedo values from building surface pars. If present,
2429	!-- they override all less specific albedo values and force a albedo_type
2430	!-- to zero in order to take effect.
2431	IF ( building_surface_pars_f%from_file ) THEN
2432	DO m = 1, surf_usm_h%ns
2433	i = surf_usm_h%i(m)
2434	j = surf_usm_h%j(m)
2435	k = surf_usm_h%k(m)
2436	!
2437	!-- Iterate over surfaces in column, check height and orientation
2438	DO is = building_surface_pars_f%index_ji(1,j,i), &
2439	building_surface_pars_f%index_ji(2,j,i)
2440	IF ( building_surface_pars_f%coords(4,is) == -surf_usm_h%koff .AND. &
2441	building_surface_pars_f%coords(1,is) == k ) THEN
2442
2443	IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /= &
2444	building_surface_pars_f%fill ) THEN
2445	surf_usm_h%albedo(ind_veg_wall,m) = &
2446	building_surface_pars_f%pars(ind_s_alb_b_wall,is)
2447	surf_usm_h%albedo_type(ind_veg_wall,m) = 0
2448	ENDIF
2449
2450	IF ( building_surface_pars_f%pars(ind_s_alb_l_wall,is) /= &
2451	building_surface_pars_f%fill ) THEN
2452	surf_usm_h%aldir(ind_veg_wall,m) = &
2453	building_surface_pars_f%pars(ind_s_alb_l_wall,is)
2454	surf_usm_h%aldif(ind_veg_wall,m) = &
2455	building_surface_pars_f%pars(ind_s_alb_l_wall,is)
2456	surf_usm_h%albedo_type(ind_veg_wall,m) = 0
2457	ENDIF
2458
2459	IF ( building_surface_pars_f%pars(ind_s_alb_s_wall,is) /= &
2460	building_surface_pars_f%fill ) THEN
2461	surf_usm_h%asdir(ind_veg_wall,m) = &
2462	building_surface_pars_f%pars(ind_s_alb_s_wall,is)
2463	surf_usm_h%asdif(ind_veg_wall,m) = &
2464	building_surface_pars_f%pars(ind_s_alb_s_wall,is)
2465	surf_usm_h%albedo_type(ind_veg_wall,m) = 0
2466	ENDIF
2467
2468	IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /= &
2469	building_surface_pars_f%fill ) THEN
2470	surf_usm_h%albedo(ind_wat_win,m) = &
2471	building_surface_pars_f%pars(ind_s_alb_b_win,is)
2472	surf_usm_h%albedo_type(ind_wat_win,m) = 0
2473	ENDIF
2474
2475	IF ( building_surface_pars_f%pars(ind_s_alb_l_win,is) /= &
2476	building_surface_pars_f%fill ) THEN
2477	surf_usm_h%aldir(ind_wat_win,m) = &
2478	building_surface_pars_f%pars(ind_s_alb_l_win,is)
2479	surf_usm_h%aldif(ind_wat_win,m) = &
2480	building_surface_pars_f%pars(ind_s_alb_l_win,is)
2481	surf_usm_h%albedo_type(ind_wat_win,m) = 0
2482	ENDIF
2483
2484	IF ( building_surface_pars_f%pars(ind_s_alb_s_win,is) /= &
2485	building_surface_pars_f%fill ) THEN
2486	surf_usm_h%asdir(ind_wat_win,m) = &
2487	building_surface_pars_f%pars(ind_s_alb_s_win,is)
2488	surf_usm_h%asdif(ind_wat_win,m) = &
2489	building_surface_pars_f%pars(ind_s_alb_s_win,is)
2490	surf_usm_h%albedo_type(ind_wat_win,m) = 0
2491	ENDIF
2492
2493	IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /= &
2494	building_surface_pars_f%fill ) THEN
2495	surf_usm_h%albedo(ind_pav_green,m) = &
2496	building_surface_pars_f%pars(ind_s_alb_b_green,is)
2497	surf_usm_h%albedo_type(ind_pav_green,m) = 0
2498	ENDIF
2499
2500	IF ( building_surface_pars_f%pars(ind_s_alb_l_green,is) /= &
2501	building_surface_pars_f%fill ) THEN
2502	surf_usm_h%aldir(ind_pav_green,m) = &
2503	building_surface_pars_f%pars(ind_s_alb_l_green,is)
2504	surf_usm_h%aldif(ind_pav_green,m) = &
2505	building_surface_pars_f%pars(ind_s_alb_l_green,is)
2506	surf_usm_h%albedo_type(ind_pav_green,m) = 0
2507	ENDIF
2508
2509	IF ( building_surface_pars_f%pars(ind_s_alb_s_green,is) /= &
2510	building_surface_pars_f%fill ) THEN
2511	surf_usm_h%asdir(ind_pav_green,m) = &
2512	building_surface_pars_f%pars(ind_s_alb_s_green,is)
2513	surf_usm_h%asdif(ind_pav_green,m) = &
2514	building_surface_pars_f%pars(ind_s_alb_s_green,is)
2515	surf_usm_h%albedo_type(ind_pav_green,m) = 0
2516	ENDIF
2517
2518	EXIT ! surface was found and processed
2519	ENDIF
2520	ENDDO
2521	ENDDO
2522
2523	DO l = 0, 3
2524	DO m = 1, surf_usm_v(l)%ns
2525	i = surf_usm_v(l)%i(m)
2526	j = surf_usm_v(l)%j(m)
2527	k = surf_usm_v(l)%k(m)
2528	!
2529	!-- Iterate over surfaces in column, check height and orientation
2530	DO is = building_surface_pars_f%index_ji(1,j,i), &
2531	building_surface_pars_f%index_ji(2,j,i)
2532	IF ( building_surface_pars_f%coords(5,is) == -surf_usm_v(l)%joff .AND. &
2533	building_surface_pars_f%coords(6,is) == -surf_usm_v(l)%ioff .AND. &
2534	building_surface_pars_f%coords(1,is) == k ) THEN
2535
2536	IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /= &
2537	building_surface_pars_f%fill ) THEN
2538	surf_usm_v(l)%albedo(ind_veg_wall,m) = &
2539	building_surface_pars_f%pars(ind_s_alb_b_wall,is)
2540	surf_usm_v(l)%albedo_type(ind_veg_wall,m) = 0
2541	ENDIF
2542
2543	IF ( building_surface_pars_f%pars(ind_s_alb_l_wall,is) /= &
2544	building_surface_pars_f%fill ) THEN
2545	surf_usm_v(l)%aldir(ind_veg_wall,m) = &
2546	building_surface_pars_f%pars(ind_s_alb_l_wall,is)
2547	surf_usm_v(l)%aldif(ind_veg_wall,m) = &
2548	building_surface_pars_f%pars(ind_s_alb_l_wall,is)
2549	surf_usm_v(l)%albedo_type(ind_veg_wall,m) = 0
2550	ENDIF
2551
2552	IF ( building_surface_pars_f%pars(ind_s_alb_s_wall,is) /= &
2553	building_surface_pars_f%fill ) THEN
2554	surf_usm_v(l)%asdir(ind_veg_wall,m) = &
2555	building_surface_pars_f%pars(ind_s_alb_s_wall,is)
2556	surf_usm_v(l)%asdif(ind_veg_wall,m) = &
2557	building_surface_pars_f%pars(ind_s_alb_s_wall,is)
2558	surf_usm_v(l)%albedo_type(ind_veg_wall,m) = 0
2559	ENDIF
2560
2561	IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /= &
2562	building_surface_pars_f%fill ) THEN
2563	surf_usm_v(l)%albedo(ind_wat_win,m) = &
2564	building_surface_pars_f%pars(ind_s_alb_b_win,is)
2565	surf_usm_v(l)%albedo_type(ind_wat_win,m) = 0
2566	ENDIF
2567
2568	IF ( building_surface_pars_f%pars(ind_s_alb_l_win,is) /= &
2569	building_surface_pars_f%fill ) THEN
2570	surf_usm_v(l)%aldir(ind_wat_win,m) = &
2571	building_surface_pars_f%pars(ind_s_alb_l_win,is)
2572	surf_usm_v(l)%aldif(ind_wat_win,m) = &
2573	building_surface_pars_f%pars(ind_s_alb_l_win,is)
2574	surf_usm_v(l)%albedo_type(ind_wat_win,m) = 0
2575	ENDIF
2576
2577	IF ( building_surface_pars_f%pars(ind_s_alb_s_win,is) /= &
2578	building_surface_pars_f%fill ) THEN
2579	surf_usm_v(l)%asdir(ind_wat_win,m) = &
2580	building_surface_pars_f%pars(ind_s_alb_s_win,is)
2581	surf_usm_v(l)%asdif(ind_wat_win,m) = &
2582	building_surface_pars_f%pars(ind_s_alb_s_win,is)
2583	surf_usm_v(l)%albedo_type(ind_wat_win,m) = 0
2584	ENDIF
2585
2586	IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /= &
2587	building_surface_pars_f%fill ) THEN
2588	surf_usm_v(l)%albedo(ind_pav_green,m) = &
2589	building_surface_pars_f%pars(ind_s_alb_b_green,is)
2590	surf_usm_v(l)%albedo_type(ind_pav_green,m) = 0
2591	ENDIF
2592
2593	IF ( building_surface_pars_f%pars(ind_s_alb_l_green,is) /= &
2594	building_surface_pars_f%fill ) THEN
2595	surf_usm_v(l)%aldir(ind_pav_green,m) = &
2596	building_surface_pars_f%pars(ind_s_alb_l_green,is)
2597	surf_usm_v(l)%aldif(ind_pav_green,m) = &
2598	building_surface_pars_f%pars(ind_s_alb_l_green,is)
2599	surf_usm_v(l)%albedo_type(ind_pav_green,m) = 0
2600	ENDIF
2601
2602	IF ( building_surface_pars_f%pars(ind_s_alb_s_green,is) /= &
2603	building_surface_pars_f%fill ) THEN
2604	surf_usm_v(l)%asdir(ind_pav_green,m) = &
2605	building_surface_pars_f%pars(ind_s_alb_s_green,is)
2606	surf_usm_v(l)%asdif(ind_pav_green,m) = &
2607	building_surface_pars_f%pars(ind_s_alb_s_green,is)
2608	surf_usm_v(l)%albedo_type(ind_pav_green,m) = 0
2609	ENDIF
2610
2611	EXIT ! surface was found and processed
2612	ENDIF
2613	ENDDO
2614	ENDDO
2615	ENDDO
2616	ENDIF
2617
2618	!
2619	!-- Calculate initial values of current (cosine of) the zenith angle and
2620	!-- whether the sun is up
2621	CALL get_date_time( time_since_reference_point, &
2622	day_of_year=day_of_year, &
2623	second_of_day=second_of_day )
2624	CALL calc_zenith( day_of_year, second_of_day )
2625	!
2626	!-- Calculate initial surface albedo for different surfaces
2627	IF ( .NOT. constant_albedo ) THEN
2628	#if defined( __netcdf )
2629	!
2630	!-- Horizontally aligned natural and urban surfaces
2631	CALL calc_albedo( surf_lsm_h )
2632	CALL calc_albedo( surf_usm_h )
2633	!
2634	!-- Vertically aligned natural and urban surfaces
2635	DO l = 0, 3
2636	CALL calc_albedo( surf_lsm_v(l) )
2637	CALL calc_albedo( surf_usm_v(l) )
2638	ENDDO
2639	#endif
2640	ELSE
2641	!
2642	!-- Initialize sun-inclination independent spectral albedos
2643	!-- Horizontal surfaces
2644	IF ( surf_lsm_h%ns > 0 ) THEN
2645	surf_lsm_h%rrtm_aldir = surf_lsm_h%aldir
2646	surf_lsm_h%rrtm_asdir = surf_lsm_h%asdir
2647	surf_lsm_h%rrtm_aldif = surf_lsm_h%aldif
2648	surf_lsm_h%rrtm_asdif = surf_lsm_h%asdif
2649	ENDIF
2650	IF ( surf_usm_h%ns > 0 ) THEN
2651	surf_usm_h%rrtm_aldir = surf_usm_h%aldir
2652	surf_usm_h%rrtm_asdir = surf_usm_h%asdir
2653	surf_usm_h%rrtm_aldif = surf_usm_h%aldif
2654	surf_usm_h%rrtm_asdif = surf_usm_h%asdif
2655	ENDIF
2656	!
2657	!-- Vertical surfaces
2658	DO l = 0, 3
2659	IF ( surf_lsm_v(l)%ns > 0 ) THEN
2660	surf_lsm_v(l)%rrtm_aldir = surf_lsm_v(l)%aldir
2661	surf_lsm_v(l)%rrtm_asdir = surf_lsm_v(l)%asdir
2662	surf_lsm_v(l)%rrtm_aldif = surf_lsm_v(l)%aldif
2663	surf_lsm_v(l)%rrtm_asdif = surf_lsm_v(l)%asdif
2664	ENDIF
2665	IF ( surf_usm_v(l)%ns > 0 ) THEN
2666	surf_usm_v(l)%rrtm_aldir = surf_usm_v(l)%aldir
2667	surf_usm_v(l)%rrtm_asdir = surf_usm_v(l)%asdir
2668	surf_usm_v(l)%rrtm_aldif = surf_usm_v(l)%aldif
2669	surf_usm_v(l)%rrtm_asdif = surf_usm_v(l)%asdif
2670	ENDIF
2671	ENDDO
2672
2673	ENDIF
2674
2675	!
2676	!-- Allocate 3d arrays of radiative fluxes and heating rates
2677	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
2678	ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2679	rad_sw_in = 0.0_wp
2680	ENDIF
2681
2682	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
2683	ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2684	ENDIF
2685
2686	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
2687	ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2688	rad_sw_out = 0.0_wp
2689	ENDIF
2690
2691	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
2692	ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2693	ENDIF
2694
2695	IF ( .NOT. ALLOCATED ( rad_sw_hr ) ) THEN
2696	ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2697	rad_sw_hr = 0.0_wp
2698	ENDIF
2699
2700	IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) ) THEN
2701	ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2702	rad_sw_hr_av = 0.0_wp
2703	ENDIF
2704
2705	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) ) THEN
2706	ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2707	rad_sw_cs_hr = 0.0_wp
2708	ENDIF
2709
2710	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) ) THEN
2711	ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2712	rad_sw_cs_hr_av = 0.0_wp
2713	ENDIF
2714
2715	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
2716	ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2717	rad_lw_in = 0.0_wp
2718	ENDIF
2719
2720	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
2721	ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2722	ENDIF
2723
2724	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
2725	ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2726	rad_lw_out = 0.0_wp
2727	ENDIF
2728
2729	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
2730	ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2731	ENDIF
2732
2733	IF ( .NOT. ALLOCATED ( rad_lw_hr ) ) THEN
2734	ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2735	rad_lw_hr = 0.0_wp
2736	ENDIF
2737
2738	IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) ) THEN
2739	ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2740	rad_lw_hr_av = 0.0_wp
2741	ENDIF
2742
2743	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) ) THEN
2744	ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2745	rad_lw_cs_hr = 0.0_wp
2746	ENDIF
2747
2748	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) ) THEN
2749	ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2750	rad_lw_cs_hr_av = 0.0_wp
2751	ENDIF
2752
2753	ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2754	ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2755	rad_sw_cs_in = 0.0_wp
2756	rad_sw_cs_out = 0.0_wp
2757
2758	ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2759	ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2760	rad_lw_cs_in = 0.0_wp
2761	rad_lw_cs_out = 0.0_wp
2762
2763	!
2764	!-- Allocate 1-element array for surface temperature
2765	!-- (RRTMG anticipates an array as passed argument).
2766	ALLOCATE ( rrtm_tsfc(1) )
2767	!
2768	!-- Allocate surface emissivity.
2769	!-- Values will be given directly before calling rrtm_lw.
2770	ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
2771
2772	!
2773	!-- Initialize RRTMG, before check if files are existent
2774	INQUIRE( FILE='rrtmg_lw.nc', EXIST=lw_exists )
2775	IF ( .NOT. lw_exists ) THEN
2776	message_string = 'Input file rrtmg_lw.nc' // &
2777	'&for rrtmg missing. ' // &
2778	'&Please provide <jobname>_lsw file in the INPUT directory.'
2779	CALL message( 'radiation_init', 'PA0583', 1, 2, 0, 6, 0 )
2780	ENDIF
2781	INQUIRE( FILE='rrtmg_sw.nc', EXIST=sw_exists )
2782	IF ( .NOT. sw_exists ) THEN
2783	message_string = 'Input file rrtmg_sw.nc' // &
2784	'&for rrtmg missing. ' // &
2785	'&Please provide <jobname>_rsw file in the INPUT directory.'
2786	CALL message( 'radiation_init', 'PA0584', 1, 2, 0, 6, 0 )
2787	ENDIF
2788
2789	IF ( lw_radiation ) CALL rrtmg_lw_ini ( c_p )
2790	IF ( sw_radiation ) CALL rrtmg_sw_ini ( c_p )
2791
2792	!
2793	!-- Set input files for RRTMG
2794	INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists)
2795	IF ( .NOT. snd_exists ) THEN
2796	rrtm_input_file = "rrtmg_lw.nc"
2797	ENDIF
2798
2799	!
2800	!-- Read vertical layers for RRTMG from sounding data
2801	!-- The routine provides nzt_rad, hyp_snd(1:nzt_rad),
2802	!-- t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
2803	!-- rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
2804	CALL read_sounding_data
2805
2806	!
2807	!-- Read trace gas profiles from file. This routine provides
2808	!-- the rrtm_ arrays (1:nzt_rad+1)
2809	CALL read_trace_gas_data
2810	#endif
2811	ENDIF
2812	!
2813	!-- Initializaion actions exclusively required for external
2814	!-- radiation forcing
2815	IF ( radiation_scheme == 'external' ) THEN
2816	!
2817	!-- Open the radiation input file. Note, for child domain, a dynamic
2818	!-- input file is often not provided. In order to do not need to
2819	!-- duplicate the dynamic input file just for the radiation input, take
2820	!-- it from the dynamic file for the parent if not available for the
2821	!-- child domain(s). In this case this is possible because radiation
2822	!-- input should be the same for each model.
2823	INQUIRE( FILE = TRIM( input_file_dynamic ), &
2824	EXIST = radiation_input_root_domain )
2825
2826	IF ( .NOT. input_pids_dynamic .AND. &
2827	.NOT. radiation_input_root_domain ) THEN
2828	message_string = 'In case of external radiation forcing ' // &
2829	'a dynamic input file is required. If no ' // &
2830	'dynamic input for the child domain(s) is ' // &
2831	'provided, at least one for the root domain ' // &
2832	'is needed.'
2833	CALL message( 'radiation_init', 'PA0315', 1, 2, 0, 6, 0 )
2834	ENDIF
2835	#if defined( __netcdf )
2836	!
2837	!-- Open dynamic input file for child domain if available, else, open
2838	!-- dynamic input file for the root domain.
2839	IF ( input_pids_dynamic ) THEN
2840	CALL open_read_file( TRIM( input_file_dynamic ) // &
2841	TRIM( coupling_char ), &
2842	pids_id )
2843	ELSEIF ( radiation_input_root_domain ) THEN
2844	CALL open_read_file( TRIM( input_file_dynamic ), &
2845	pids_id )
2846	ENDIF
2847
2848	CALL inquire_num_variables( pids_id, num_var_pids )
2849	!
2850	!-- Allocate memory to store variable names and read them
2851	ALLOCATE( vars_pids(1:num_var_pids) )
2852	CALL inquire_variable_names( pids_id, vars_pids )
2853	!
2854	!-- Input time dimension.
2855	IF ( check_existence( vars_pids, 'time_rad' ) ) THEN
2856	CALL get_dimension_length( pids_id, ntime, 'time_rad' )
2857
2858	ALLOCATE( time_rad_f%var1d(0:ntime-1) )
2859	!
2860	!-- Read variable
2861	CALL get_variable( pids_id, 'time_rad', time_rad_f%var1d )
2862
2863	time_rad_f%from_file = .TRUE.
2864	ENDIF
2865	!
2866	!-- Input shortwave downwelling.
2867	IF ( check_existence( vars_pids, 'rad_sw_in' ) ) THEN
2868	!
2869	!-- Get _FillValue attribute
2870	CALL get_attribute( pids_id, char_fill, rad_sw_in_f%fill, &
2871	.FALSE., 'rad_sw_in' )
2872	!
2873	!-- Get level-of-detail
2874	CALL get_attribute( pids_id, char_lod, rad_sw_in_f%lod, &
2875	.FALSE., 'rad_sw_in' )
2876	!
2877	!-- Level-of-detail 1 - radiation depends only on time_rad
2878	IF ( rad_sw_in_f%lod == 1 ) THEN
2879	ALLOCATE( rad_sw_in_f%var1d(0:ntime-1) )
2880	CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var1d )
2881	rad_sw_in_f%from_file = .TRUE.
2882	!
2883	!-- Level-of-detail 2 - radiation depends on time_rad, y, x
2884	ELSEIF ( rad_sw_in_f%lod == 2 ) THEN
2885	ALLOCATE( rad_sw_in_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )
2886
2887	CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var3d, &
2888	nxl, nxr, nys, nyn, 0, ntime-1 )
2889
2890	rad_sw_in_f%from_file = .TRUE.
2891	ELSE
2892	message_string = '"rad_sw_in" has no valid lod attribute'
2893	CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
2894	ENDIF
2895	ENDIF
2896	!
2897	!-- Input longwave downwelling.
2898	IF ( check_existence( vars_pids, 'rad_lw_in' ) ) THEN
2899	!
2900	!-- Get _FillValue attribute
2901	CALL get_attribute( pids_id, char_fill, rad_lw_in_f%fill, &
2902	.FALSE., 'rad_lw_in' )
2903	!
2904	!-- Get level-of-detail
2905	CALL get_attribute( pids_id, char_lod, rad_lw_in_f%lod, &
2906	.FALSE., 'rad_lw_in' )
2907	!
2908	!-- Level-of-detail 1 - radiation depends only on time_rad
2909	IF ( rad_lw_in_f%lod == 1 ) THEN
2910	ALLOCATE( rad_lw_in_f%var1d(0:ntime-1) )
2911	CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var1d )
2912	rad_lw_in_f%from_file = .TRUE.
2913	!
2914	!-- Level-of-detail 2 - radiation depends on time_rad, y, x
2915	ELSEIF ( rad_lw_in_f%lod == 2 ) THEN
2916	ALLOCATE( rad_lw_in_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )
2917
2918	CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var3d, &
2919	nxl, nxr, nys, nyn, 0, ntime-1 )
2920
2921	rad_lw_in_f%from_file = .TRUE.
2922	ELSE
2923	message_string = '"rad_lw_in" has no valid lod attribute'
2924	CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
2925	ENDIF
2926	ENDIF
2927	!
2928	!-- Input shortwave downwelling, diffuse part.
2929	IF ( check_existence( vars_pids, 'rad_sw_in_dif' ) ) THEN
2930	!
2931	!-- Read _FillValue attribute
2932	CALL get_attribute( pids_id, char_fill, rad_sw_in_dif_f%fill, &
2933	.FALSE., 'rad_sw_in_dif' )
2934	!
2935	!-- Get level-of-detail
2936	CALL get_attribute( pids_id, char_lod, rad_sw_in_dif_f%lod, &
2937	.FALSE., 'rad_sw_in_dif' )
2938	!
2939	!-- Level-of-detail 1 - radiation depends only on time_rad
2940	IF ( rad_sw_in_dif_f%lod == 1 ) THEN
2941	ALLOCATE( rad_sw_in_dif_f%var1d(0:ntime-1) )
2942	CALL get_variable( pids_id, 'rad_sw_in_dif', &
2943	rad_sw_in_dif_f%var1d )
2944	rad_sw_in_dif_f%from_file = .TRUE.
2945	!
2946	!-- Level-of-detail 2 - radiation depends on time_rad, y, x
2947	ELSEIF ( rad_sw_in_dif_f%lod == 2 ) THEN
2948	ALLOCATE( rad_sw_in_dif_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )
2949
2950	CALL get_variable( pids_id, 'rad_sw_in_dif', &
2951	rad_sw_in_dif_f%var3d, &
2952	nxl, nxr, nys, nyn, 0, ntime-1 )
2953
2954	rad_sw_in_dif_f%from_file = .TRUE.
2955	ELSE
2956	message_string = '"rad_sw_in_dif" has no valid lod attribute'
2957	CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
2958	ENDIF
2959	ENDIF
2960	!
2961	!-- Finally, close the input file and deallocate temporary arrays
2962	DEALLOCATE( vars_pids )
2963
2964	CALL close_input_file( pids_id )
2965	#endif
2966	!
2967	!-- Make some consistency checks.
2968	IF ( .NOT. rad_sw_in_f%from_file .OR. &
2969	.NOT. rad_lw_in_f%from_file ) THEN
2970	message_string = 'In case of external radiation forcing ' // &
2971	'both, rad_sw_in and rad_lw_in are required.'
2972	CALL message( 'radiation_init', 'PA0195', 1, 2, 0, 6, 0 )
2973	ENDIF
2974
2975	IF ( .NOT. time_rad_f%from_file ) THEN
2976	message_string = 'In case of external radiation forcing ' // &
2977	'dimension time_rad is required.'
2978	CALL message( 'radiation_init', 'PA0196', 1, 2, 0, 6, 0 )
2979	ENDIF
2980
2981	CALL get_date_time( 0.0_wp, second_of_day=second_of_day )
2982
2983	IF ( end_time - spinup_time > time_rad_f%var1d(ntime-1) ) THEN
2984	message_string = 'External radiation forcing does not cover ' // &
2985	'the entire simulation time.'
2986	CALL message( 'radiation_init', 'PA0314', 1, 2, 0, 6, 0 )
2987	ENDIF
2988	!
2989	!-- Check for fill values in radiation
2990	IF ( ALLOCATED( rad_sw_in_f%var1d ) ) THEN
2991	IF ( ANY( rad_sw_in_f%var1d == rad_sw_in_f%fill ) ) THEN
2992	message_string = 'External radiation array "rad_sw_in" ' // &
2993	'must not contain any fill values.'
2994	CALL message( 'radiation_init', 'PA0197', 1, 2, 0, 6, 0 )
2995	ENDIF
2996	ENDIF
2997
2998	IF ( ALLOCATED( rad_lw_in_f%var1d ) ) THEN
2999	IF ( ANY( rad_lw_in_f%var1d == rad_lw_in_f%fill ) ) THEN
3000	message_string = 'External radiation array "rad_lw_in" ' // &
3001	'must not contain any fill values.'
3002	CALL message( 'radiation_init', 'PA0198', 1, 2, 0, 6, 0 )
3003	ENDIF
3004	ENDIF
3005
3006	IF ( ALLOCATED( rad_sw_in_dif_f%var1d ) ) THEN
3007	IF ( ANY( rad_sw_in_dif_f%var1d == rad_sw_in_dif_f%fill ) ) THEN
3008	message_string = 'External radiation array "rad_sw_in_dif" ' //&
3009	'must not contain any fill values.'
3010	CALL message( 'radiation_init', 'PA0199', 1, 2, 0, 6, 0 )
3011	ENDIF
3012	ENDIF
3013
3014	IF ( ALLOCATED( rad_sw_in_f%var3d ) ) THEN
3015	IF ( ANY( rad_sw_in_f%var3d == rad_sw_in_f%fill ) ) THEN
3016	message_string = 'External radiation array "rad_sw_in" ' // &
3017	'must not contain any fill values.'
3018	CALL message( 'radiation_init', 'PA0197', 1, 2, 0, 6, 0 )
3019	ENDIF
3020	ENDIF
3021
3022	IF ( ALLOCATED( rad_lw_in_f%var3d ) ) THEN
3023	IF ( ANY( rad_lw_in_f%var3d == rad_lw_in_f%fill ) ) THEN
3024	message_string = 'External radiation array "rad_lw_in" ' // &
3025	'must not contain any fill values.'
3026	CALL message( 'radiation_init', 'PA0198', 1, 2, 0, 6, 0 )
3027	ENDIF
3028	ENDIF
3029
3030	IF ( ALLOCATED( rad_sw_in_dif_f%var3d ) ) THEN
3031	IF ( ANY( rad_sw_in_dif_f%var3d == rad_sw_in_dif_f%fill ) ) THEN
3032	message_string = 'External radiation array "rad_sw_in_dif" ' //&
3033	'must not contain any fill values.'
3034	CALL message( 'radiation_init', 'PA0199', 1, 2, 0, 6, 0 )
3035	ENDIF
3036	ENDIF
3037	!
3038	!-- Currently, 2D external radiation input is not possible in
3039	!-- combination with topography where average radiation is used.
3040	IF ( ( rad_lw_in_f%lod == 2 .OR. rad_sw_in_f%lod == 2 .OR. &
3041	rad_sw_in_dif_f%lod == 2 ) .AND. average_radiation ) THEN
3042	message_string = 'External radiation with lod = 2 is currently '//&
3043	'not possible with average_radiation = .T..'
3044	CALL message( 'radiation_init', 'PA0670', 1, 2, 0, 6, 0 )
3045	ENDIF
3046	!
3047	!-- All radiation input should have the same level of detail. The sum
3048	!-- of lods divided by the number of available radiation arrays must be
3049	!-- 1 (if all are lod = 1) or 2 (if all are lod = 2).
3050	IF ( REAL( MERGE( rad_lw_in_f%lod, 0, rad_lw_in_f%from_file ) + &
3051	MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) + &
3052	MERGE( rad_sw_in_dif_f%lod, 0, rad_sw_in_dif_f%from_file ),&
3053	KIND = wp ) / &
3054	( MERGE( 1.0_wp, 0.0_wp, rad_lw_in_f%from_file ) + &
3055	MERGE( 1.0_wp, 0.0_wp, rad_sw_in_f%from_file ) + &
3056	MERGE( 1.0_wp, 0.0_wp, rad_sw_in_dif_f%from_file ) ) &
3057	/= 1.0_wp .AND. &
3058	REAL( MERGE( rad_lw_in_f%lod, 0, rad_lw_in_f%from_file ) + &
3059	MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) + &
3060	MERGE( rad_sw_in_dif_f%lod, 0, rad_sw_in_dif_f%from_file ),&
3061	KIND = wp ) / &
3062	( MERGE( 1.0_wp, 0.0_wp, rad_lw_in_f%from_file ) + &
3063	MERGE( 1.0_wp, 0.0_wp, rad_sw_in_f%from_file ) + &
3064	MERGE( 1.0_wp, 0.0_wp, rad_sw_in_dif_f%from_file ) ) &
3065	/= 2.0_wp ) THEN
3066	message_string = 'External radiation input should have the same '//&
3067	'lod.'
3068	CALL message( 'radiation_init', 'PA0673', 1, 2, 0, 6, 0 )
3069	ENDIF
3070
3071	ENDIF
3072	!
3073	!-- Perform user actions if required
3074	CALL user_init_radiation
3075
3076	!
3077	!-- Calculate radiative fluxes at model start
3078	SELECT CASE ( TRIM( radiation_scheme ) )
3079
3080	CASE ( 'rrtmg' )
3081	CALL radiation_rrtmg
3082
3083	CASE ( 'clear-sky' )
3084	CALL radiation_clearsky
3085
3086	CASE ( 'constant' )
3087	CALL radiation_constant
3088
3089	CASE ( 'external' )
3090	!
3091	!-- During spinup apply clear-sky model
3092	IF ( time_since_reference_point < 0.0_wp ) THEN
3093	CALL radiation_clearsky
3094	ELSE
3095	CALL radiation_external
3096	ENDIF
3097
3098	CASE DEFAULT
3099
3100	END SELECT
3101
3102	!
3103	!-- Find all discretized apparent solar positions for radiation interaction.
3104	IF ( radiation_interactions ) CALL radiation_presimulate_solar_pos
3105
3106	!
3107	!-- If required, read or calculate and write out the SVF
3108	IF ( radiation_interactions .AND. read_svf) THEN
3109	!
3110	!-- Read sky-view factors and further required data from file
3111	CALL radiation_read_svf()
3112
3113	ELSEIF ( radiation_interactions .AND. .NOT. read_svf) THEN
3114	!
3115	!-- calculate SFV and CSF
3116	CALL radiation_calc_svf()
3117	ENDIF
3118
3119	IF ( radiation_interactions .AND. write_svf) THEN
3120	!
3121	!-- Write svf, csf svfsurf and csfsurf data to file
3122	CALL radiation_write_svf()
3123	ENDIF
3124
3125	!
3126	!-- Adjust radiative fluxes. In case of urban and land surfaces, also
3127	!-- call an initial interaction.
3128	IF ( radiation_interactions ) THEN
3129	CALL radiation_interaction
3130	ENDIF
3131
3132	IF ( debug_output ) CALL debug_message( 'radiation_init', 'end' )
3133
3134	RETURN !todo: remove, I don't see what we need this for here
3135
3136	END SUBROUTINE radiation_init
3137
3138
3139	!------------------------------------------------------------------------------!
3140	! Description:
3141	! ------------
3142	!> A simple clear sky radiation model
3143	!------------------------------------------------------------------------------!
3144	SUBROUTINE radiation_external
3145
3146	IMPLICIT NONE
3147
3148	INTEGER(iwp) :: l !< running index for surface orientation
3149	INTEGER(iwp) :: t !< index of current timestep
3150	INTEGER(iwp) :: tm !< index of previous timestep
3151
3152	LOGICAL :: horizontal !< flag indicating treatment of horinzontal surfaces
3153
3154	REAL(wp) :: fac_dt !< interpolation factor
3155	REAL(wp) :: second_of_day_init !< second of the day at model start
3156
3157	TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine
3158
3159	!
3160	!-- Calculate current zenith angle
3161	CALL get_date_time( time_since_reference_point, &
3162	day_of_year=day_of_year, &
3163	second_of_day=second_of_day )
3164	CALL calc_zenith( day_of_year, second_of_day )
3165	!
3166	!-- Interpolate external radiation on current timestep
3167	IF ( time_since_reference_point <= 0.0_wp ) THEN
3168	t = 0
3169	tm = 0
3170	fac_dt = 0
3171	ELSE
3172	CALL get_date_time( 0.0_wp, second_of_day=second_of_day_init )
3173	t = 0
3174	DO WHILE ( time_rad_f%var1d(t) <= time_since_reference_point )
3175	t = t + 1
3176	ENDDO
3177
3178	tm = MAX( t-1, 0 )
3179
3180	fac_dt = ( time_since_reference_point &
3181	- time_rad_f%var1d(tm) + dt_3d ) &
3182	/ ( time_rad_f%var1d(t) - time_rad_f%var1d(tm) )
3183	fac_dt = MIN( 1.0_wp, fac_dt )
3184	ENDIF
3185	!
3186	!-- Call clear-sky calculation for each surface orientation.
3187	!-- First, horizontal surfaces
3188	horizontal = .TRUE.
3189	surf => surf_lsm_h
3190	CALL radiation_external_surf
3191	surf => surf_usm_h
3192	CALL radiation_external_surf
3193	horizontal = .FALSE.
3194	!
3195	!-- Vertical surfaces
3196	DO l = 0, 3
3197	surf => surf_lsm_v(l)
3198	CALL radiation_external_surf
3199	surf => surf_usm_v(l)
3200	CALL radiation_external_surf
3201	ENDDO
3202
3203	CONTAINS
3204
3205	SUBROUTINE radiation_external_surf
3206
3207	USE control_parameters
3208
3209	IMPLICIT NONE
3210
3211	INTEGER(iwp) :: i !< grid index along x-dimension
3212	INTEGER(iwp) :: j !< grid index along y-dimension
3213	INTEGER(iwp) :: k !< grid index along z-dimension
3214	INTEGER(iwp) :: m !< running index for surface elements
3215
3216	REAL(wp) :: lw_in !< downwelling longwave radiation, interpolated value
3217	REAL(wp) :: sw_in !< downwelling shortwave radiation, interpolated value
3218	REAL(wp) :: sw_in_dif !< downwelling diffuse shortwave radiation, interpolated value
3219
3220	IF ( surf%ns < 1 ) RETURN
3221	!
3222	!-- level-of-detail = 1. Note, here it must be distinguished between
3223	!-- averaged radiation and non-averaged radiation for the upwelling
3224	!-- fluxes.
3225	IF ( rad_sw_in_f%lod == 1 ) THEN
3226
3227	sw_in = ( 1.0_wp - fac_dt ) * rad_sw_in_f%var1d(tm) &
3228	+ fac_dt * rad_sw_in_f%var1d(t)
3229
3230	lw_in = ( 1.0_wp - fac_dt ) * rad_lw_in_f%var1d(tm) &
3231	+ fac_dt * rad_lw_in_f%var1d(t)
3232	!
3233	!-- Limit shortwave incoming radiation to positive values, in order
3234	!-- to overcome possible observation errors.
3235	sw_in = MAX( 0.0_wp, sw_in )
3236	sw_in = MERGE( sw_in, 0.0_wp, sun_up )
3237
3238	surf%rad_sw_in = sw_in
3239	surf%rad_lw_in = lw_in
3240
3241	IF ( average_radiation ) THEN
3242	surf%rad_sw_out = albedo_urb * surf%rad_sw_in
3243
3244	surf%rad_lw_out = emissivity_urb * sigma_sb * t_rad_urb**4 &
3245	+ ( 1.0_wp - emissivity_urb ) * surf%rad_lw_in
3246
3247	surf%rad_net = surf%rad_sw_in - surf%rad_sw_out &
3248	+ surf%rad_lw_in - surf%rad_lw_out
3249
3250	surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb &
3251	* sigma_sb &
3252	* t_rad_urb**3
3253	ELSE
3254	DO m = 1, surf%ns
3255	k = surf%k(m)
3256	surf%rad_sw_out(m) = ( surf%frac(ind_veg_wall,m) * &
3257	surf%albedo(ind_veg_wall,m) &
3258	+ surf%frac(ind_pav_green,m) * &
3259	surf%albedo(ind_pav_green,m) &
3260	+ surf%frac(ind_wat_win,m) * &
3261	surf%albedo(ind_wat_win,m) ) &
3262	* surf%rad_sw_in(m)
3263
3264	surf%rad_lw_out(m) = ( surf%frac(ind_veg_wall,m) * &
3265	surf%emissivity(ind_veg_wall,m) &
3266	+ surf%frac(ind_pav_green,m) * &
3267	surf%emissivity(ind_pav_green,m) &
3268	+ surf%frac(ind_wat_win,m) * &
3269	surf%emissivity(ind_wat_win,m) &
3270	) &
3271	* sigma_sb &
3272	* ( surf%pt_surface(m) * exner(k) )**4
3273
3274	surf%rad_lw_out_change_0(m) = &
3275	( surf%frac(ind_veg_wall,m) * &
3276	surf%emissivity(ind_veg_wall,m) &
3277	+ surf%frac(ind_pav_green,m) * &
3278	surf%emissivity(ind_pav_green,m) &
3279	+ surf%frac(ind_wat_win,m) * &
3280	surf%emissivity(ind_wat_win,m) &
3281	) * 4.0_wp * sigma_sb &
3282	* ( surf%pt_surface(m) * exner(k) )**3
3283	ENDDO
3284
3285	ENDIF
3286	!
3287	!-- If diffuse shortwave radiation is available, store it on
3288	!-- the respective files.
3289	IF ( rad_sw_in_dif_f%from_file ) THEN
3290	sw_in_dif= ( 1.0_wp - fac_dt ) * rad_sw_in_dif_f%var1d(tm) &
3291	+ fac_dt * rad_sw_in_dif_f%var1d(t)
3292
3293	IF ( ALLOCATED( rad_sw_in_diff ) ) rad_sw_in_diff = sw_in_dif
3294	IF ( ALLOCATED( rad_sw_in_dir ) ) rad_sw_in_dir = sw_in &
3295	- sw_in_dif
3296	!
3297	!-- Diffuse longwave radiation equals the total downwelling
3298	!-- longwave radiation
3299	IF ( ALLOCATED( rad_lw_in_diff ) ) rad_lw_in_diff = lw_in
3300	ENDIF
3301	!
3302	!-- level-of-detail = 2
3303	ELSE
3304
3305	DO m = 1, surf%ns
3306	i = surf%i(m)
3307	j = surf%j(m)
3308	k = surf%k(m)
3309
3310	surf%rad_sw_in(m) = ( 1.0_wp - fac_dt ) &
3311	* rad_sw_in_f%var3d(tm,j,i) &
3312	+ fac_dt * rad_sw_in_f%var3d(t,j,i)
3313	!
3314	!-- Limit shortwave incoming radiation to positive values, in
3315	!-- order to overcome possible observation errors.
3316	surf%rad_sw_in(m) = MAX( 0.0_wp, surf%rad_sw_in(m) )
3317	surf%rad_sw_in(m) = MERGE( surf%rad_sw_in(m), 0.0_wp, sun_up )
3318
3319	surf%rad_lw_in(m) = ( 1.0_wp - fac_dt ) &
3320	* rad_lw_in_f%var3d(tm,j,i) &
3321	+ fac_dt * rad_lw_in_f%var3d(t,j,i)
3322	!
3323	!-- Weighted average according to surface fraction.
3324	surf%rad_sw_out(m) = ( surf%frac(ind_veg_wall,m) * &
3325	surf%albedo(ind_veg_wall,m) &
3326	+ surf%frac(ind_pav_green,m) * &
3327	surf%albedo(ind_pav_green,m) &
3328	+ surf%frac(ind_wat_win,m) * &
3329	surf%albedo(ind_wat_win,m) ) &
3330	* surf%rad_sw_in(m)
3331
3332	surf%rad_lw_out(m) = ( surf%frac(ind_veg_wall,m) * &
3333	surf%emissivity(ind_veg_wall,m) &
3334	+ surf%frac(ind_pav_green,m) * &
3335	surf%emissivity(ind_pav_green,m) &
3336	+ surf%frac(ind_wat_win,m) * &
3337	surf%emissivity(ind_wat_win,m) &
3338	) &
3339	* sigma_sb &
3340	* ( surf%pt_surface(m) * exner(k) )**4
3341
3342	surf%rad_lw_out_change_0(m) = &
3343	( surf%frac(ind_veg_wall,m) * &
3344	surf%emissivity(ind_veg_wall,m) &
3345	+ surf%frac(ind_pav_green,m) * &
3346	surf%emissivity(ind_pav_green,m) &
3347	+ surf%frac(ind_wat_win,m) * &
3348	surf%emissivity(ind_wat_win,m) &
3349	) * 4.0_wp * sigma_sb &
3350	* ( surf%pt_surface(m) * exner(k) )**3
3351
3352	surf%rad_net(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m) &
3353	+ surf%rad_lw_in(m) - surf%rad_lw_out(m)
3354	!
3355	!-- If diffuse shortwave radiation is available, store it on
3356	!-- the respective files.
3357	IF ( rad_sw_in_dif_f%from_file ) THEN
3358	IF ( ALLOCATED( rad_sw_in_diff ) ) &
3359	rad_sw_in_diff(j,i) = ( 1.0_wp - fac_dt ) &
3360	* rad_sw_in_dif_f%var3d(tm,j,i) &
3361	+ fac_dt * rad_sw_in_dif_f%var3d(t,j,i)
3362	!
3363	!-- dir = sw_in - sw_in_dif.
3364	IF ( ALLOCATED( rad_sw_in_dir ) ) &
3365	rad_sw_in_dir(j,i) = surf%rad_sw_in(m) - &
3366	rad_sw_in_diff(j,i)
3367	!
3368	!-- Diffuse longwave radiation equals the total downwelling
3369	!-- longwave radiation
3370	IF ( ALLOCATED( rad_lw_in_diff ) ) &
3371	rad_lw_in_diff(j,i) = surf%rad_lw_in(m)
3372	ENDIF
3373
3374	ENDDO
3375
3376	ENDIF
3377	!
3378	!-- Store radiation also on 2D arrays, which are still used for
3379	!-- direct-diffuse splitting. Note, this is only required
3380	!-- for horizontal surfaces, which covers all x,y position.
3381	IF ( horizontal ) THEN
3382	DO m = 1, surf%ns
3383	i = surf%i(m)
3384	j = surf%j(m)
3385
3386	rad_sw_in(0,j,i) = surf%rad_sw_in(m)
3387	rad_lw_in(0,j,i) = surf%rad_lw_in(m)
3388	rad_sw_out(0,j,i) = surf%rad_sw_out(m)
3389	rad_lw_out(0,j,i) = surf%rad_lw_out(m)
3390	ENDDO
3391	ENDIF
3392
3393	END SUBROUTINE radiation_external_surf
3394
3395	END SUBROUTINE radiation_external
3396
3397	!------------------------------------------------------------------------------!
3398	! Description:
3399	! ------------
3400	!> A simple clear sky radiation model
3401	!------------------------------------------------------------------------------!
3402	SUBROUTINE radiation_clearsky
3403
3404	IMPLICIT NONE
3405
3406	INTEGER(iwp) :: l !< running index for surface orientation
3407
3408	LOGICAL :: horizontal !< flag indicating treatment of horinzontal surfaces
3409
3410	REAL(wp) :: pt1 !< potential temperature at first grid level or mean value at urban layer top
3411	REAL(wp) :: pt1_l !< potential temperature at first grid level or mean value at urban layer top at local subdomain
3412	REAL(wp) :: ql1 !< liquid water mixing ratio at first grid level or mean value at urban layer top
3413	REAL(wp) :: ql1_l !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain
3414
3415	TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine
3416
3417	!
3418	!-- Calculate current zenith angle
3419	CALL get_date_time( time_since_reference_point, &
3420	day_of_year=day_of_year, &
3421	second_of_day=second_of_day )
3422	CALL calc_zenith( day_of_year, second_of_day )
3423
3424	!
3425	!-- Calculate sky transmissivity
3426	sky_trans = 0.6_wp + 0.2_wp * cos_zenith
3427
3428	!
3429	!-- Calculate value of the Exner function at model surface
3430	!
3431	!-- In case averaged radiation is used, calculate mean temperature and
3432	!-- liquid water mixing ratio at the urban-layer top.
3433	IF ( average_radiation ) THEN
3434	pt1 = 0.0_wp
3435	IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = 0.0_wp
3436
3437	pt1_l = SUM( pt(nz_urban_t,nys:nyn,nxl:nxr) )
3438	IF ( bulk_cloud_model .OR. cloud_droplets ) ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) )
3439
3440	#if defined( __parallel )
3441	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
3442	CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
3443	IF ( ierr /= 0 ) THEN
3444	WRITE(9,*) 'Error MPI_AllReduce1:', ierr, pt1_l, pt1
3445	FLUSH(9)
3446	ENDIF
3447
3448	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
3449	CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
3450	IF ( ierr /= 0 ) THEN
3451	WRITE(9,*) 'Error MPI_AllReduce2:', ierr, ql1_l, ql1
3452	FLUSH(9)
3453	ENDIF
3454	ENDIF
3455	#else
3456	pt1 = pt1_l
3457	IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = ql1_l
3458	#endif
3459
3460	IF ( bulk_cloud_model .OR. cloud_droplets ) pt1 = pt1 + lv_d_cp / exner(nz_urban_t) * ql1
3461	!
3462	!-- Finally, divide by number of grid points
3463	pt1 = pt1 / REAL( ( nx + 1 ) * ( ny + 1 ), KIND=wp )
3464	ENDIF
3465	!
3466	!-- Call clear-sky calculation for each surface orientation.
3467	!-- First, horizontal surfaces
3468	horizontal = .TRUE.
3469	surf => surf_lsm_h
3470	CALL radiation_clearsky_surf
3471	surf => surf_usm_h
3472	CALL radiation_clearsky_surf
3473	horizontal = .FALSE.
3474	!
3475	!-- Vertical surfaces
3476	DO l = 0, 3
3477	surf => surf_lsm_v(l)
3478	CALL radiation_clearsky_surf
3479	surf => surf_usm_v(l)
3480	CALL radiation_clearsky_surf
3481	ENDDO
3482
3483	CONTAINS
3484
3485	SUBROUTINE radiation_clearsky_surf
3486
3487	IMPLICIT NONE
3488
3489	INTEGER(iwp) :: i !< index x-direction
3490	INTEGER(iwp) :: j !< index y-direction
3491	INTEGER(iwp) :: k !< index z-direction
3492	INTEGER(iwp) :: m !< running index for surface elements
3493
3494	IF ( surf%ns < 1 ) RETURN
3495
3496	!
3497	!-- Calculate radiation fluxes and net radiation (rad_net) assuming
3498	!-- homogeneous urban radiation conditions.
3499	IF ( average_radiation ) THEN
3500
3501	k = nz_urban_t
3502
3503	surf%rad_sw_in = solar_constant * sky_trans * cos_zenith
3504	surf%rad_sw_out = albedo_urb * surf%rad_sw_in
3505
3506	surf%rad_lw_in = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k+1))**4
3507
3508	surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4 &
3509	+ (1.0_wp - emissivity_urb) * surf%rad_lw_in
3510
3511	surf%rad_net = surf%rad_sw_in - surf%rad_sw_out &
3512	+ surf%rad_lw_in - surf%rad_lw_out
3513
3514	surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb * sigma_sb &
3515	* (t_rad_urb)**3
3516
3517	!
3518	!-- Calculate radiation fluxes and net radiation (rad_net) for each surface
3519	!-- element.
3520	ELSE
3521
3522	DO m = 1, surf%ns
3523	i = surf%i(m)
3524	j = surf%j(m)
3525	k = surf%k(m)
3526
3527	surf%rad_sw_in(m) = solar_constant * sky_trans * cos_zenith
3528
3529	!
3530	!-- Weighted average according to surface fraction.
3531	!-- ATTENTION: when radiation interactions are switched on the
3532	!-- calculated fluxes below are not actually used as they are
3533	!-- overwritten in radiation_interaction.
3534	surf%rad_sw_out(m) = ( surf%frac(ind_veg_wall,m) * &
3535	surf%albedo(ind_veg_wall,m) &
3536	+ surf%frac(ind_pav_green,m) * &
3537	surf%albedo(ind_pav_green,m) &
3538	+ surf%frac(ind_wat_win,m) * &
3539	surf%albedo(ind_wat_win,m) ) &
3540	* surf%rad_sw_in(m)
3541
3542	surf%rad_lw_out(m) = ( surf%frac(ind_veg_wall,m) * &
3543	surf%emissivity(ind_veg_wall,m) &
3544	+ surf%frac(ind_pav_green,m) * &
3545	surf%emissivity(ind_pav_green,m) &
3546	+ surf%frac(ind_wat_win,m) * &
3547	surf%emissivity(ind_wat_win,m) &
3548	) &
3549	* sigma_sb &
3550	* ( surf%pt_surface(m) * exner(nzb) )**4
3551
3552	surf%rad_lw_out_change_0(m) = &
3553	( surf%frac(ind_veg_wall,m) * &
3554	surf%emissivity(ind_veg_wall,m) &
3555	+ surf%frac(ind_pav_green,m) * &
3556	surf%emissivity(ind_pav_green,m) &
3557	+ surf%frac(ind_wat_win,m) * &
3558	surf%emissivity(ind_wat_win,m) &
3559	) * 4.0_wp * sigma_sb &
3560	* ( surf%pt_surface(m) * exner(nzb) )** 3
3561
3562
3563	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
3564	pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i)
3565	surf%rad_lw_in(m) = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k))**4
3566	ELSE
3567	surf%rad_lw_in(m) = emissivity_atm_clsky * sigma_sb * (pt(k,j,i) * exner(k))**4
3568	ENDIF
3569
3570	surf%rad_net(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m) &
3571	+ surf%rad_lw_in(m) - surf%rad_lw_out(m)
3572
3573	ENDDO
3574
3575	ENDIF
3576
3577	!
3578	!-- Fill out values in radiation arrays. Note, this is only required
3579	!-- for horizontal surfaces, which covers all x,y position.
3580	IF ( horizontal ) THEN
3581	DO m = 1, surf%ns
3582	i = surf%i(m)
3583	j = surf%j(m)
3584	rad_sw_in(0,j,i) = surf%rad_sw_in(m)
3585	rad_sw_out(0,j,i) = surf%rad_sw_out(m)
3586	rad_lw_in(0,j,i) = surf%rad_lw_in(m)
3587	rad_lw_out(0,j,i) = surf%rad_lw_out(m)
3588	ENDDO
3589	ENDIF
3590
3591	END SUBROUTINE radiation_clearsky_surf
3592
3593	END SUBROUTINE radiation_clearsky
3594
3595
3596	!------------------------------------------------------------------------------!
3597	! Description:
3598	! ------------
3599	!> This scheme keeps the prescribed net radiation constant during the run
3600	!------------------------------------------------------------------------------!
3601	SUBROUTINE radiation_constant
3602
3603
3604	IMPLICIT NONE
3605
3606	INTEGER(iwp) :: l !< running index for surface orientation
3607
3608	LOGICAL :: horizontal !< flag indicating treatment of horinzontal surfaces
3609
3610	REAL(wp) :: pt1 !< potential temperature at first grid level or mean value at urban layer top
3611	REAL(wp) :: pt1_l !< potential temperature at first grid level or mean value at urban layer top at local subdomain
3612	REAL(wp) :: ql1 !< liquid water mixing ratio at first grid level or mean value at urban layer top
3613	REAL(wp) :: ql1_l !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain
3614
3615	TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine
3616
3617	!
3618	!-- In case averaged radiation is used, calculate mean temperature and
3619	!-- liquid water mixing ratio at the urban-layer top.
3620	IF ( average_radiation ) THEN
3621	pt1 = 0.0_wp
3622	IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = 0.0_wp
3623
3624	pt1_l = SUM( pt(nz_urban_t,nys:nyn,nxl:nxr) )
3625	IF ( bulk_cloud_model .OR. cloud_droplets ) ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) )
3626
3627	#if defined( __parallel )
3628	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
3629	CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
3630	IF ( ierr /= 0 ) THEN
3631	WRITE(9,*) 'Error MPI_AllReduce3:', ierr, pt1_l, pt1
3632	FLUSH(9)
3633	ENDIF
3634	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
3635	CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
3636	IF ( ierr /= 0 ) THEN
3637	WRITE(9,*) 'Error MPI_AllReduce4:', ierr, ql1_l, ql1
3638	FLUSH(9)
3639	ENDIF
3640	ENDIF
3641	#else
3642	pt1 = pt1_l
3643	IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = ql1_l
3644	#endif
3645	IF ( bulk_cloud_model .OR. cloud_droplets ) pt1 = pt1 + lv_d_cp / exner(nz_urban_t+1) * ql1
3646	!
3647	!-- Finally, divide by number of grid points
3648	pt1 = pt1 / REAL( ( nx + 1 ) * ( ny + 1 ), KIND=wp )
3649	ENDIF
3650
3651	!
3652	!-- First, horizontal surfaces
3653	horizontal = .TRUE.
3654	surf => surf_lsm_h
3655	CALL radiation_constant_surf
3656	surf => surf_usm_h
3657	CALL radiation_constant_surf
3658	horizontal = .FALSE.
3659	!
3660	!-- Vertical surfaces
3661	DO l = 0, 3
3662	surf => surf_lsm_v(l)
3663	CALL radiation_constant_surf
3664	surf => surf_usm_v(l)
3665	CALL radiation_constant_surf
3666	ENDDO
3667
3668	CONTAINS
3669
3670	SUBROUTINE radiation_constant_surf
3671
3672	IMPLICIT NONE
3673
3674	INTEGER(iwp) :: i !< index x-direction
3675	INTEGER(iwp) :: ioff !< offset between surface element and adjacent grid point along x
3676	INTEGER(iwp) :: j !< index y-direction
3677	INTEGER(iwp) :: joff !< offset between surface element and adjacent grid point along y
3678	INTEGER(iwp) :: k !< index z-direction
3679	INTEGER(iwp) :: koff !< offset between surface element and adjacent grid point along z
3680	INTEGER(iwp) :: m !< running index for surface elements
3681
3682	IF ( surf%ns < 1 ) RETURN
3683
3684	!-- Calculate homogenoeus urban radiation fluxes
3685	IF ( average_radiation ) THEN
3686
3687	surf%rad_net = net_radiation
3688
3689	surf%rad_lw_in = emissivity_atm_clsky * sigma_sb * (pt1 * exner(nz_urban_t+1))**4
3690
3691	surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4 &
3692	+ ( 1.0_wp - emissivity_urb ) & ! shouldn't be this a bulk value -- emissivity_urb?
3693	* surf%rad_lw_in
3694
3695	surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb * sigma_sb &
3696	* t_rad_urb**3
3697
3698	surf%rad_sw_in = ( surf%rad_net - surf%rad_lw_in &
3699	+ surf%rad_lw_out ) &
3700	/ ( 1.0_wp - albedo_urb )
3701
3702	surf%rad_sw_out = albedo_urb * surf%rad_sw_in
3703
3704	!
3705	!-- Calculate radiation fluxes for each surface element
3706	ELSE
3707	!
3708	!-- Determine index offset between surface element and adjacent
3709	!-- atmospheric grid point
3710	ioff = surf%ioff
3711	joff = surf%joff
3712	koff = surf%koff
3713
3714	!
3715	!-- Prescribe net radiation and estimate the remaining radiative fluxes
3716	DO m = 1, surf%ns
3717	i = surf%i(m)
3718	j = surf%j(m)
3719	k = surf%k(m)
3720
3721	surf%rad_net(m) = net_radiation
3722
3723	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
3724	pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i)
3725	surf%rad_lw_in(m) = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k))**4
3726	ELSE
3727	surf%rad_lw_in(m) = emissivity_atm_clsky * sigma_sb * &
3728	( pt(k,j,i) * exner(k) )**4
3729	ENDIF
3730
3731	!
3732	!-- Weighted average according to surface fraction.
3733	surf%rad_lw_out(m) = ( surf%frac(ind_veg_wall,m) * &
3734	surf%emissivity(ind_veg_wall,m) &
3735	+ surf%frac(ind_pav_green,m) * &
3736	surf%emissivity(ind_pav_green,m) &
3737	+ surf%frac(ind_wat_win,m) * &
3738	surf%emissivity(ind_wat_win,m) &
3739	) &
3740	* sigma_sb &
3741	* ( surf%pt_surface(m) * exner(nzb) )**4
3742
3743	surf%rad_sw_in(m) = ( surf%rad_net(m) - surf%rad_lw_in(m) &
3744	+ surf%rad_lw_out(m) ) &
3745	/ ( 1.0_wp - &
3746	( surf%frac(ind_veg_wall,m) * &
3747	surf%albedo(ind_veg_wall,m) &
3748	+ surf%frac(ind_pav_green,m) * &
3749	surf%albedo(ind_pav_green,m) &
3750	+ surf%frac(ind_wat_win,m) * &
3751	surf%albedo(ind_wat_win,m) ) &
3752	)
3753
3754	surf%rad_sw_out(m) = ( surf%frac(ind_veg_wall,m) * &
3755	surf%albedo(ind_veg_wall,m) &
3756	+ surf%frac(ind_pav_green,m) * &
3757	surf%albedo(ind_pav_green,m) &
3758	+ surf%frac(ind_wat_win,m) * &
3759	surf%albedo(ind_wat_win,m) ) &
3760	* surf%rad_sw_in(m)
3761
3762	ENDDO
3763
3764	ENDIF
3765
3766	!
3767	!-- Fill out values in radiation arrays. Note, this is only required
3768	!-- for horizontal surfaces, which covers all x,y position.
3769	IF ( horizontal ) THEN
3770	DO m = 1, surf%ns
3771	i = surf%i(m)
3772	j = surf%j(m)
3773	rad_sw_in(0,j,i) = surf%rad_sw_in(m)
3774	rad_sw_out(0,j,i) = surf%rad_sw_out(m)
3775	rad_lw_in(0,j,i) = surf%rad_lw_in(m)
3776	rad_lw_out(0,j,i) = surf%rad_lw_out(m)
3777	ENDDO
3778	ENDIF
3779
3780	END SUBROUTINE radiation_constant_surf
3781
3782
3783	END SUBROUTINE radiation_constant
3784
3785	!------------------------------------------------------------------------------!
3786	! Description:
3787	! ------------
3788	!> Header output for radiation model
3789	!------------------------------------------------------------------------------!
3790	SUBROUTINE radiation_header ( io )
3791
3792
3793	IMPLICIT NONE
3794
3795	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
3796
3797
3798
3799	!
3800	!-- Write radiation model header
3801	WRITE( io, 3 )
3802
3803	IF ( radiation_scheme == "constant" ) THEN
3804	WRITE( io, 4 ) net_radiation
3805	ELSEIF ( radiation_scheme == "clear-sky" ) THEN
3806	WRITE( io, 5 )
3807	ELSEIF ( radiation_scheme == "rrtmg" ) THEN
3808	WRITE( io, 6 )
3809	IF ( .NOT. lw_radiation ) WRITE( io, 10 )
3810	IF ( .NOT. sw_radiation ) WRITE( io, 11 )
3811	ELSEIF ( radiation_scheme == "external" ) THEN
3812	WRITE( io, 14 )
3813	ENDIF
3814
3815	IF ( albedo_type_f%from_file .OR. vegetation_type_f%from_file .OR. &
3816	pavement_type_f%from_file .OR. water_type_f%from_file .OR. &
3817	building_type_f%from_file ) THEN
3818	WRITE( io, 13 )
3819	ELSE
3820	IF ( albedo_type == 0 ) THEN
3821	WRITE( io, 7 ) albedo
3822	ELSE
3823	WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
3824	ENDIF
3825	ENDIF
3826	IF ( constant_albedo ) THEN
3827	WRITE( io, 9 )
3828	ENDIF
3829
3830	WRITE( io, 12 ) dt_radiation
3831
3832
3833	3 FORMAT (//' Radiation model information:'/ &
3834	' ----------------------------'/)
3835	4 FORMAT (' --> Using constant net radiation: net_radiation = ', F6.2, &
3836	// 'W/m**2')
3837	5 FORMAT (' --> Simple radiation scheme for clear sky is used (no clouds,',&
3838	' default)')
3839	6 FORMAT (' --> RRTMG scheme is used')
3840	7 FORMAT (/' User-specific surface albedo: albedo =', F6.3)
3841	8 FORMAT (/' Albedo is set for land surface type: ', A)
3842	9 FORMAT (/' --> Albedo is fixed during the run')
3843	10 FORMAT (/' --> Longwave radiation is disabled')
3844	11 FORMAT (/' --> Shortwave radiation is disabled.')
3845	12 FORMAT (' Timestep: dt_radiation = ', F6.2, ' s')
3846	13 FORMAT (/' Albedo is set individually for each xy-location, according ', &
3847	'to given surface type.')
3848	14 FORMAT (' --> External radiation forcing is used')
3849
3850
3851	END SUBROUTINE radiation_header
3852
3853
3854	!------------------------------------------------------------------------------!
3855	! Description:
3856	! ------------
3857	!> Parin for &radiation_parameters for radiation model and RTM
3858	!------------------------------------------------------------------------------!
3859	SUBROUTINE radiation_parin
3860
3861
3862	IMPLICIT NONE
3863
3864	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
3865
3866	NAMELIST /radiation_par/ albedo, albedo_lw_dif, albedo_lw_dir, &
3867	albedo_sw_dif, albedo_sw_dir, albedo_type, &
3868	constant_albedo, dt_radiation, emissivity, &
3869	lw_radiation, max_raytracing_dist, &
3870	min_irrf_value, mrt_geom, mrt_geom_params, &
3871	mrt_include_sw, mrt_nlevels, &
3872	mrt_skip_roof, net_radiation, nrefsteps, &
3873	plant_lw_interact, rad_angular_discretization,&
3874	radiation_interactions_on, radiation_scheme, &
3875	raytrace_discrete_azims, &
3876	raytrace_discrete_elevs, raytrace_mpi_rma, &
3877	trace_fluxes_above, &
3878	skip_time_do_radiation, surface_reflections, &
3879	svfnorm_report_thresh, sw_radiation, &
3880	unscheduled_radiation_calls
3881
3882
3883	NAMELIST /radiation_parameters/ albedo, albedo_lw_dif, albedo_lw_dir, &
3884	albedo_sw_dif, albedo_sw_dir, albedo_type, &
3885