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

Last change on this file since 4429 was 4429, checked in by raasch, 19 months ago

serial (non-MPI) test case added, several bugfixes for the serial mode

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