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

Last change on this file since 4198 was 4198, checked in by gronemeier, 21 months ago

Add check into radiation_check_parameters if rotation angle is set to 0. Using the radiation model for a rotated model domain is not allowed due to missing implementation within the radiation model.

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