Home

Context Navigation

radiation_model_mod.f90 @ 4245

Last change on this file since 4245 was 4245, checked in by pavelkrc, 4 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-4244
/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: 545.4 KB

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