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

Last change on this file since 4210 was 4210, checked in by suehring, 2 years ago

Radiation: revise steering of splitting diffuse and direct radiation; revise some checks; optimize mapping of radiation components onto 2D arrays, avoid unnecessary operations

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