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

Last change on this file since 4337 was 4337, checked in by Giersch, 17 months ago

Albedo indices for building_surface_pars are now declared as parameters to prevent an error if the gfortran compiler with -Werror=unused-value is used

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