Home

Context Navigation

radiation_model_mod.f90 @ 2963

Last change on this file since 2963 was 2963, checked in by suehring, 7 years ago
Minor revision of static input file checks, bugfix in initialization of surface-fractions in LSM; minor bugfix in initialization of albedo at window-surfaces; for clearer access of albedo and emissivity introduce index for vegetation/wall, pavement/green-wall and water/window surfaces
Property svn:keywords set to `Id`
File size: 373.6 KB

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