Home

Context Navigation

radiation_model_mod.f90 @ 4529

Last change on this file since 4529 was 4529, checked in by moh.hefny, 5 years ago

New features in coupling RTM-radiation model: consider PC-LW interactions, consider PC emiss. in the effective urban emiss., new algorithm for better performance, minor formatting changes

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