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

Last change on this file since 4286 was 4286, checked in by resler, 20 months ago

Merge branch resler into trunk

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