Home

Context Navigation

radiation_model_mod.f90 @ 4493

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

Merge brach resler

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