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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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