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

Last change on this file since 4313 was 4291, checked in by moh.hefny, 19 months ago

RTM is ON in case of biometeorology even if there is no vertical surfaces and/or trees to enable biommet-output in flat surfaces domain

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