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

Last change on this file since 4238 was 4238, checked in by suehring, 20 months ago

Building data base: Indoor-model parameters for some building types adjusted in order to avoid unrealistically high indoor temperatures (S. Rissmann); Indoor model: Bugfix in determination of minimum facade height and in location message, avoid divisions by zero, minor optimizations; radiation model: Modify check in order to avoid equality comparisons of floating points

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