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

Last change on this file since 4441 was 4441, checked in by suehring, 19 months ago

Change order of dimension in surface arrays %frac, %emissivity and %albedo to allow for better vectorization in the radiation interactions; Set back turbulent length scale to 8 x grid spacing in the parametrized mode for the synthetic turbulence generator (was accidentally changed in last commit)

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