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

Last change on this file since 4517 was 4517, checked in by raasch, 14 months ago

added restart with MPI-IO for reading local arrays

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