Home

Context Navigation

radiation_model_mod.f90 @ 4286

Last change on this file since 4286 was 4286, checked in by resler, 5 years ago

Merge branch resler into trunk

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