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

Last change on this file since 4643 was 4643, checked in by pavelkrc, 4 months ago

Bugfix: typo in index of surf%* arrays

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