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

Last change on this file since 4340 was 4340, checked in by Giersch, 4 years ago

Topography closed channel flow with symmetric boundaries implemented, ID tag in radiation module corrected

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