Home

Context Navigation

radiation_model_mod.f90 @ 4238

Last change on this file since 4238 was 4238, checked in by suehring, 5 years ago

Building data base: Indoor-model parameters for some building types adjusted in order to avoid unrealistically high indoor temperatures (S. Rissmann); Indoor model: Bugfix in determination of minimum facade height and in location message, avoid divisions by zero, minor optimizations; radiation model: Modify check in order to avoid equality comparisons of floating points

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.4 KB

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