Home

Context Navigation

radiation_model_mod.f90 @ 4291

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

RTM is ON in case of biometeorology even if there is no vertical surfaces and/or trees to enable biommet-output in flat surfaces domain

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