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

Last change on this file since 4768 was 4768, checked in by suehring, 4 years ago

Enable 3D data output also with 64-bit precision

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