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

Last change on this file since 4644 was 4644, checked in by pavelkrc, 3 months ago

Bugfix: remove invalid reordering of vertical surfaces in RTM

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