Home

Context Navigation

radiation_model_mod.f90 @ 2930

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