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

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

Revise steering of surface albedo initialization when albedo_pars is provided

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