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

Last change on this file since 4227 was 4227, checked in by gronemeier, 20 months ago

implement new palm_date_time_mod; replaced namelist parameters time_utc_init and day_of_year_init by origin_date_time

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