source: palm/trunk/SOURCE/radiation_model_mod.f90

Last change on this file was 4899, checked in by raasch, 3 months ago

small adjustments regarding r4897, bugfix for r4898 (missing preprocessor directive added)

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