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

Last change on this file since 4400 was 4400, checked in by suehring, 20 months ago

Revision of the virtual-measurement module: data input from NetCDF file; removed binary output - instead parallel NetCDF output using the new data-output module; variable attributes added; further variables added that can be sampled, file connections added; Functions for coordinate transformation moved to basic_constants_and_equations; netcdf_data_input_mod: unused routines netcdf_data_input_att and netcdf_data_input_var removed; new routines to inquire fill values added; Preprocessing script (palm_cvd) to setup virtual-measurement input files provided; postprocessor combine_virtual_measurements deleted

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