Home

Context Navigation

radiation_model_mod.f90 @ 4226

Last change on this file since 4226 was 4226, checked in by suehring, 4 years ago

Offline nesting: data input modularized; time variable is defined relative to time_utc_init, so that input data is correctly mapped to actual model time; checks rephrased and new checks for the time dimension added; Netcdf input: routine to retrieve dimension length renamed

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-4166
/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: 526.8 KB

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