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

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

last commit documented

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