source: palm/trunk/SOURCE/radiation_model_mod.f90 @ 3424

Last change on this file since 3424 was 3424, checked in by gronemeier, 4 years ago

bugfix: add rad_lw_in, rad_lw_out, rad_sw_out to radiation_check_data_output

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