Home

Context Navigation

radiation_model_mod.f90 @ 2696

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