Home

Context Navigation

radiation_model_mod.f90 @ 2698

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