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

Last change on this file since 4190 was 4190, checked in by suehring, 2 years ago

Implement external radiation forcing also for level-of-detail = 2 (horizontally 2D radiation)

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