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

Last change on this file since 4598 was 4587, checked in by pavelkrc, 10 months ago

Marked RTM version 3.1, bugfixes

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