Home

Context Navigation

radiation_model_mod.f90 @ 4531

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

update test cases and bugfix in radiation module in non-parallel case

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

/palm/branches/chemistry/SOURCE/radiation_model_mod.f90	2047-3190,3218-3297
/palm/branches/forwind/SOURCE/radiation_model_mod.f90	1564-1913
/palm/branches/mosaik_M2/radiation_model_mod.f90	2360-3471
/palm/branches/palm4u/SOURCE/radiation_model_mod.f90	2540-2692
/palm/branches/radiation/SOURCE/radiation_model_mod.f90	2081-3493
/palm/branches/rans/SOURCE/radiation_model_mod.f90	2078-3128
/palm/branches/resler/SOURCE/radiation_model_mod.f90	2023-4492
/palm/branches/salsa/SOURCE/radiation_model_mod.f90	2503-3460
/palm/branches/fricke/SOURCE/radiation_model_mod.f90	942-977
/palm/branches/hoffmann/SOURCE/radiation_model_mod.f90	989-1052
/palm/branches/letzel/masked_output/SOURCE/radiation_model_mod.f90	296-409
/palm/branches/suehring/radiation_model_mod.f90	423-666

File size: 566.2 KB

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