Home

Context Navigation

radiation_model_mod.f90 @ 4495

Last change on this file since 4495 was 4495, checked in by raasch, 5 years ago

restart data handling with MPI-IO added, first part

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