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

Last change on this file since 4360 was 4360, checked in by suehring, 16 months ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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