Home

Context Navigation

radiation_model_mod.f90 @ 4227

Last change on this file since 4227 was 4227, checked in by gronemeier, 5 years ago

implement new palm_date_time_mod; replaced namelist parameters time_utc_init and day_of_year_init by origin_date_time

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