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

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

Urban-surface model and indoor model: Default building parameters updated and extended; revision of summer- and wintertime-dependent indoor parameters (responsible S. Rissmann); extension of albedo parameters to different building facades and roofs; first and second wall layer initialized with individual building properties rather than with with the same properties (heat capacity and conductivity)

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