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

Last change on this file since 4783 was 4783, checked in by raasch, 20 months ago

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