Home

Context Navigation

radiation_model_mod.f90 @ 3046

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