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

Last change on this file since 4392 was 4392, checked in by pavelkrc, 20 months ago

Merge branch resler (updated documentation, optional flux tracing, bugfix)

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