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

Last change on this file since 4452 was 4452, checked in by suehring, 3 years ago

Bugfix in calc_albedo

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