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

Last change on this file since 4271 was 4271, checked in by maronga, 21 months ago

bugfix in calculation of snow albedo

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