Home

Context Navigation

radiation_model_mod.f90 @ 2809

Last change on this file since 2809 was 2809, checked in by schwenkel, 7 years ago
Bugfix for gfortran: Replace the function C_SIZEOF with STORAGE_SIZE
Property svn:keywords set to `Id`
File size: 375.2 KB

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